FV3 Bundle
Zeeman_Utility.f90
Go to the documentation of this file.
1 !
2 MODULE zeeman_utility
3  ! Description:
4  ! containing routines to load geomagnetic field Lookup table, compute
5  ! geomagnetic field and its angles relative to the wave propagation
6  ! direction.
7  !
8  ! History:
9  ! ---- -------
10  ! 11/07/2007 created by Y. Han, NOAA/JCSDA
11  ! 11/24/2009 modified for CRTM
12  !
13 
14  ! -----------------
15  ! Environment setup
16  ! -----------------
17  ! Module use
18  USE type_kinds, ONLY: fp
19  USE message_handler, ONLY: success, failure, display_message
20  USE file_utility, ONLY: get_lun, file_exists
21 
22  ! Disable implicit typing
23  IMPLICIT NONE
24 
25  ! ------------
26  ! Visibilities
27  ! ------------
28  PRIVATE
29  PUBLIC :: load_bfield_lut ! function to load Geomagnetic field LUT
30  PUBLIC :: compute_bfield ! subroutine to calculate the Geomagnetic field
31  ! using a LUT and two cosine angles of the field
32  ! relative to the wave propagation direction k.
33  PUBLIC :: compute_kb_angles ! subroutine to compute two cosine angles of the
34  ! geomagnetic field relative to the wave
35  ! propagation direction k.
36 
37 
38  INTERFACE compute_bfield
39  MODULE PROCEDURE compute_bfield_f1
40  MODULE PROCEDURE compute_bfield_f2
41  MODULE PROCEDURE compute_bfield_f3
42  END INTERFACE compute_bfield
43 
44  ! Array for Earth's magnetic field
45  INTEGER, PARAMETER :: n_lat = 91, n_lon = 181
46  INTEGER, SAVE :: bfield(3, n_lat, n_lon)
47 
48  Real(fp), parameter :: degrees_to_radians = 3.141592653589793238462643_fp/180.0_fp
49 
50 CONTAINS
51 
52  !--------------------------------------------------------------------------------
53  !
54  ! NAME:
55  ! load_bfield_lut
56  !
57  ! PURPOSE:
58  ! Function to the geomagnetic filed LUT
59  !
60  ! CALLING SEQUENCE:
61  ! Error_Status = load_bfield_lut(filename_LUT)
62  !
63  ! INPUT ARGUMENT:
64  !
65  ! filename_LUT: file name for the file containing the LUT of
66  ! the Earth magnetic field.
67  ! UNITS: N/A
68  ! TYPE: character string
69  ! DIMENSION: Scale
70  ! ATTRIBUTES: INTENT(IN)
71  !----------------------------------------------------------------------------
72  Function load_bfield_lut(filename_LUT) Result(Error_Status)
73  Character(*), Intent(in) :: filename_lut
74 
75  Integer :: error_status
76 
77  ! local
78  CHARACTER(*), PARAMETER :: routine_name = 'load_bfield_lut'
79  Integer :: file_id, io_status
80  Character (len=80) :: message
81  Integer :: i, j, iskip
82 
83  error_status = success
84 
85  IF ( .NOT. file_exists( trim(filename_lut) ) )THEN
86  error_status = failure
87  message = 'File '//trim(filename_lut)//' not found.'
88  CALL display_message( routine_name, &
89  trim( message ), &
90  error_status)
91  END IF
92 
93  ! Find a free logical unit number for file access
94  file_id = get_lun()
95 
96  OPEN( file_id, file = filename_lut, &
97  & status = 'OLD', &
98  & iostat = io_status )
99 
100  IF ( io_status /= 0 ) THEN
101  error_status = failure
102  WRITE( message, '( "Error opening ", a, ". IOSTAT = ", i5 )' ) &
103  & trim(filename_lut), io_status
104  CALL display_message( routine_name, &
105  trim( message ), &
106  error_status)
107  RETURN
108  END IF
109 
110  READ(file_id, *, iostat = io_status)iskip, iskip, &
111  ((bfield(:, i, j), j=1,n_lon), i=1,n_lat)
112 
113  If(io_status /= 0) Then
114  error_status = failure
115  Write( message, '( "Error reading data from ", a, ". IOSTAT = ", i5 )' ) &
116  & trim(filename_lut), io_status
117  CALL display_message( routine_name, &
118  trim( message ), &
119  error_status)
120  Return
121  Endif
122 
123  Close( unit = file_id )
124 
125  End Function load_bfield_lut
126 
127  !--------------------------------------------------------------------------------
128  !
129  ! NAME:
130  ! compute_bfield
131  !
132  ! PURPOSE:
133  ! Subroutine to calculate the Geomagnetic field using a LUT and
134  ! two cosine angles of the field relative to the wave propagation
135  ! direction k.
136  !
137  ! CALLING SEQUENCE:
138  !
139  ! CALL compute_bfield(latitude, longitude, & ! Inputs
140  ! Bx, By, Bz, Be) ! Outputs
141  !
142  ! OR
143  !
144  ! CALL compute_bfield(latitude, & ! input
145  ! longitude, & ! input
146  ! sensor_zenang, & ! input
147  ! sensor_aziang, & ! input
148  ! Be, & ! output
149  ! cos_bkang, & ! output
150  ! cos_baziang) ! output
151  !
152  ! OR
153  !
154  ! CALL compute_bfield(latitude, & ! input
155  ! longitude, & ! input
156  ! sensor_zenang, & ! input
157  ! sensor_relative_aziang, & ! input
158  ! Julian_day, & ! input
159  ! utc_time, & ! input
160  ! Be, & ! output
161  ! cos_bkang, & ! output
162  ! cos_baziang) ! output
163  !
164  !
165  ! INPUT ARGUMENTS:
166  !
167  !
168  ! latitude : Latitude, 0 - 180 (0 North Pole; 180 - South Pole
169  ! UNITS: Degree
170  ! TYPE: real(fp)
171  ! DIMENSION: Scale
172  ! ATTRIBUTES: INTENT(IN)
173  !
174  ! longitude: longitude, 0 - 360 East
175  ! UNITS: Degree
176  ! TYPE: real(fp)
177  ! DIMENSION: Scale
178  ! ATTRIBUTES: INTENT(IN)
179  !
180  ! sensor_zenang: sensor zenith angle
181  ! UNITS: Degree
182  ! TYPE: real(fp)
183  ! DIMENSION: Scale
184  ! ATTRIBUTES: INTENT(IN)
185  !
186  ! sensor_aziang: sensor azimuth angle (starts from East,
187  ! positive counterclockwise)
188  ! UNITS: Degree
189  ! TYPE: real(fp)
190  ! DIMENSION: Scale
191  ! ATTRIBUTES: INTENT(IN)
192  !
193  ! sensor_relative_aziang: sensor azimuth angle relative to the sun azimuth
194  ! angle.
195  ! Sun_azimuth_angle - from north, East positive
196  ! sensor_relative_aziang = 90 - (sun_azimuth_angle + Sensor_aziang)
197  ! UNITS: degree
198  ! TYPE: real(fp)
199  ! DIMENSION: Scale
200  ! ATTRIBUTES: INTENT(IN)
201  !
202  ! Julian_day: Julian_Day 1=Jan 1, 365=Dec31 (366 leap year)
203  ! UNITS: day
204  ! TYPE: integer(JPIM)
205  ! DIMENSION: Scale
206  ! ATTRIBUTES: INTENT(IN)
207  !
208  ! utc_time: Universal_Time 0.00-23.999,(GMT,Z time)
209  ! UNITS: hour
210  ! TYPE: real(fp)
211  ! DIMENSION: Scale
212  ! ATTRIBUTES: INTENT(IN)
213  !
214  ! OUTPUT ARGUMENTS:
215  ! Bx: Magetic field East component
216  ! UNITS: Gauss
217  ! TYPE: real(fp)
218  ! DIMENSION: Scalar
219  ! ATTRIBUTES: INTENT(OUT)
220  !
221  ! By: Magetic field North component
222  ! UNITS: Gauss
223  ! TYPE: real(fp)
224  ! DIMENSION: Scalar
225  ! ATTRIBUTES: INTENT(OUT)
226  !
227  ! Bz: Magetic field zenith component (positive upward)
228  ! UNITS: Gauss
229  ! TYPE: real(fp)
230  ! DIMENSION: Scalar
231  ! ATTRIBUTES: INTENT(OUT)
232  !
233  ! Be: Magetic field strength (sqrt(BxBx + ByBy + BzBz))
234  ! UNITS: Gauss
235  ! TYPE: real(fp)
236  ! DIMENSION: Scalar
237  ! ATTRIBUTES: INTENT(OUT)
238  !
239  ! cos_bkang: cosine of the angle between the magnetic field Be
240  ! vector and the wave propagation direction k.
241  ! UNITS: N/A
242  ! TYPE: real(fp)
243  ! DIMENSION: Scalar
244  ! ATTRIBUTES: INTENT(OUT)
245  !
246  ! cos_baziang: cosine of the azimuth angle of the Be vector in the
247  ! (v, h, k) coordinates system, where v, h and k comprise
248  ! a right-hand orthogonal system, similar to the (x, y, z)
249  ! Catesian coordinates. The h vector is normal to the
250  ! plane containing the k and z vectors, where k points
251  ! to the wave propagation direction and z points
252  ! to the zenith. h = (z cross k)/|z cross k|. The
253  ! azimuth angle is the angle on the (v, h) plane
254  ! from the positive v axis to the projected line of the
255  ! Be vector on this plane, positive counterclockwise.
256  ! UNITS: N/A
257  ! TYPE: real(fp)
258  ! DIMENSION: Scalar
259  ! ATTRIBUTES: INTENT(OUT)
260  !
261  !--------------------------------------------------------------------------------
262 
263  SUBROUTINE compute_bfield_f1(latitude, longitude, & ! Inputs
264  Bx, By, Bz, Be) ! Outputs
265  REAL(fp), INTENT(IN) :: latitude, longitude
266  REAL(fp), INTENT(OUT) :: Bx, By, Bz, Be
267 
268  ! Local
269  REAL(fp) :: x2, w1_lat, w1_lon
270  INTEGER :: idx_lat, idx_lon
271  REAL(fp), PARAMETER :: dlat = 2.0_fp, dlon = 2.0_fp
272 
273  idx_lat = int(latitude/dlat)+1
274  IF(idx_lat >= n_lat)idx_lat = n_lat-1
275  idx_lon = int(longitude/dlat)+1
276  IF(idx_lon >= n_lon)idx_lon = n_lon-1
277 
278  x2 = REAL(idx_lat, fp)*dlat
279  w1_lat = (x2 - latitude)/dlat
280 
281  x2 = REAL(idx_lon, fp)*dlat
282  w1_lon = (x2 - longitude)/dlat
283 
284  bx = bfield_component(1, bfield, w1_lat, w1_lon, idx_lat, idx_lon)
285  by = bfield_component(2, bfield, w1_lat, w1_lon, idx_lat, idx_lon)
286  bz = bfield_component(3, bfield, w1_lat, w1_lon, idx_lat, idx_lon)
287 
288  be = sqrt(bx*bx+by*by+bz*bz)
289 
290  CONTAINS
291 
292  FUNCTION bfield_component(comp, Bfield, w1_lat, w1_lon, &
293  idx_lat, idx_lon) Result(B)
294 
295  INTEGER, INTENT(IN) :: comp
296  INTEGER, INTENT(IN) :: bfield(:,:,:)
297  REAL(fp), INTENT(IN) :: w1_lat, w1_lon
298  INTEGER, INTENT(IN) :: idx_lat, idx_lon
299 
300  REAL(fp) :: b
301 
302  REAL(fp), PARAMETER :: scale = 0.001_fp
303  REAL(fp) :: w2_lat, w2_lon
304 
305  w2_lat = 1.0_fp - w1_lat
306  w2_lon = 1.0_fp - w1_lon
307  b = (w1_lon*(w1_lat*REAL(BField(comp,idx_lat, idx_lon), fp) + &
308  w2_lat*REAL(BField(comp,idx_lat+1, idx_lon), fp)) &
309  + w2_lon*(w1_lat*REAL(BField(comp,idx_lat, idx_lon+1),fp) + &
310  w2_lat*REAL(BField(comp,idx_lat+1, idx_lon+1),fp)))*scale
311 
312  END FUNCTION bfield_component
313 
314  END SUBROUTINE compute_bfield_f1
315 
316  Subroutine compute_bfield_f2(latitude, longitude, sensor_zenang, sensor_aziang, &
317  Be, cos_bkang, cos_baziang)
318 
319  !subroutine arguments:
320  Real(fp), Intent(in) :: latitude
321  Real(fp), Intent(in) :: longitude
322  Real(fp), Intent(in) :: sensor_zenang
323  Real(fp), Intent(in) :: sensor_aziang
324  Real(fp), Intent(out) :: Be
325  Real(fp), Intent(out) :: cos_bkang
326  Real(fp), Intent(out) :: cos_baziang
327 
328  ! local variable
329  Real(fp) :: Bx, By, Bz
330 
331  ! get Earth magnetic filed from LUT
332  Call compute_bfield_f1(latitude, longitude, & ! inputs
333  bx, by, bz, be) ! outputs
334 
335  ! compute the cosines of the angle between the magnetic field Be and
336  ! propagation direction and Be's azimuth angle
337  Call compute_kb_angles(bx, by, bz, sensor_zenang, sensor_aziang, & ! Input
338  cos_bkang, cos_baziang) ! Output
339 
340  End Subroutine compute_bfield_f2
341 
342  Subroutine compute_bfield_f3(latitude, longitude, sensor_zenang, &
343  sensor_relative_aziang, julian_day, utc_time, &
344  Be, cos_bkang, cos_baziang)
345  !subroutine arguments:
346  Real(fp), Intent(in) :: latitude
347  Real(fp), Intent(in) :: longitude
348  Real(fp), Intent(in) :: sensor_zenang
349  Real(fp), Intent(in) :: sensor_relative_aziang
350  Integer, Intent(in) :: julian_day
351  Real(fp), Intent(in) :: utc_time
352  Real(fp), Intent(out) :: Be
353  Real(fp), Intent(out) :: cos_bkang
354  Real(fp), Intent(out) :: cos_baziang
355 
356  ! Local
357  Real(fp) :: Bx, By, Bz, Solar_ZA, Solar_AZ, lat, lon, sensor_aziang
358 
359  ! compute the cosines of the angle between the magnetic field Be and
360  ! propagation direction and Be's azimuth angle
361  Call compute_bfield_f1(latitude, longitude, & ! inputs
362  bx, by, bz, be) ! outputs
363 
364  lat = 90.0_fp - latitude
365  lon = longitude
366  If(lon > 180.0_fp)lon = lon - 360.0_fp
367  ! Compute Solar azimuth angle
368  Call solar_za_az(lat,lon, Real(julian_day, fp), utc_time, &
369  Solar_ZA, Solar_AZ)
370  ! Compute satellite azimuth angle (starts from East, positive counterclockwise)
371  sensor_aziang = sensor_relative_aziang + solar_az
372  sensor_aziang = 90.0_fp - sensor_aziang
373 
374  ! compute for cos_bkangle, cos_baziangle
375  Call compute_kb_angles(bx, by, bz, sensor_zenang, sensor_aziang, & ! Input
376  cos_bkang, cos_baziang) ! Output
377 
378 
379  End Subroutine compute_bfield_f3
380 
381 
382 !--------------------------------------------------------------------------------
383 !
384 ! NAME:
385 ! Compute_kb_Angles
386 !
387 ! PURPOSE:
388 ! Subroutine to calculate the cosine of the angle between the Geomagnetic
389 ! field (Be) and wave propagation direction (k) and the cosine of the
390 ! azimuth angle of the Be vector in the (v, h, k) coordinates system (see
391 ! more detailed description below)
392 !
393 !
394 ! CALLING SEQUENCE:
395 ! CALL Compute_kb_Angles(Bx, By, Bz, & ! Input
396 ! sensor_zenang, sensor_aziang, & ! Input
397 ! cos_bk_Angle, cos_baziang) ! Output
398 ! INPUT ARGUMENTS:
399 !
400 ! Bx: Magetic field East component
401 ! UNITS: Gauss
402 ! TYPE: Real(fp)
403 ! DIMENSION: Scalar
404 ! ATTRIBUTES: INTENT(IN)
405 !
406 ! By: Magetic field North component
407 ! UNITS: Gauss
408 ! TYPE: Real(fp)
409 ! DIMENSION: Scalar
410 ! ATTRIBUTES: INTENT(IN)
411 !
412 ! Bz: Magetic field zenith component (positive upward)
413 ! UNITS: Gauss
414 ! TYPE: Real(fp)
415 ! DIMENSION: Scalar
416 ! ATTRIBUTES: INTENT(IN)
417 !
418 ! sensor_zenang : sensor zenith angle
419 ! UNITS: Degree
420 ! TYPE: Real(fp)
421 ! DIMENSION: Scalar
422 ! ATTRIBUTES: INTENT(IN)
423 !
424 ! sensor_aziang : sensor zenith angle defined as the
425 ! angle from the East towards North.
426 ! UNITS: Degree
427 ! TYPE: Real(fp)
428 ! DIMENSION: Scalar
429 ! ATTRIBUTES: INTENT(IN)
430 !
431 ! OUTPUT ARGUMENTS:
432 ! cos_bkang: cosine of the angle between the magnetic field Be
433 ! vector and the wave propagation direction k.
434 ! UNITS: N/A
435 ! TYPE: real(fp)
436 ! DIMENSION: Scalar
437 ! ATTRIBUTES: INTENT(OUT)
438 !
439 ! cos_baziang: cosine of the azimuth angle of the Be vector in the
440 ! (v, h, k) coordinates system, where v, h and k comprise
441 ! a right-hand orthogonal system, similar to the (x, y, z)
442 ! Catesian coordinates. The h vector is normal to the
443 ! plane containing the k and z vectors, where k points
444 ! to the wave propagation direction and z points
445 ! to the zenith. h = (z cross k)/|z cross k|. The
446 ! azimuth angle is the angle on the (v, h) plane
447 ! from the positive v axis to the projected line of the
448 ! Be vector on this plane, positive counterclockwise.
449 ! UNITS: N/A
450 ! TYPE: Real(fp)
451 ! DIMENSION: Scalar
452 ! ATTRIBUTES: INTENT(OUT)
453 !
454 !--------------------------------------------------------------------------------
455 
456  SUBROUTINE compute_kb_angles(Bx, By, Bz, & ! Input
457  sensor_zenang, sensor_aziang, & ! Input
458  cos_bkang, cos_baziang) ! Output
459  REAL(fp), INTENT(IN) :: bx, by, bz
460  REAL(fp), INTENT(IN) :: sensor_zenang, sensor_aziang
461  REAL(fp), INTENT(OUT) :: cos_bkang, cos_baziang
462 
463  ! Local
464  REAL(fp) :: b, b_v, b_h, b_p, kx, ky, kz, &
465  sin_senza, cos_senaz, sin_senaz, cos_senza
466 
467  sin_senza = sin(sensor_zenang*degrees_to_radians)
468  cos_senza = cos(sensor_zenang*degrees_to_radians)
469  sin_senaz = sin(sensor_aziang*degrees_to_radians)
470  cos_senaz = cos(sensor_aziang*degrees_to_radians)
471 
472  ! compute k directional vector from satellite's zenith and azimuth angles
473  kx = sin_senza*cos_senaz
474  ky = sin_senza*sin_senaz
475  kz = cos_senza
476 
477  ! compute consine of the angle between the magnetic field B and k
478  b = sqrt(bx*bx+by*by+bz*bz)
479  cos_bkang = (kx*bx + ky*by + kz*bz)/b
480 
481  ! Project the B vector on the V and H plane: B_v - B component on the V
482  ! axis; B_h - B component on the H axis.
483  b_v = bx*kx*kz + by*ky*kz - bz*(ky*ky + kx*kx) ;
484  b_h = -bx*ky + by*kx
485 
486  ! compute the cosine of the azimuth angle
487  b_p = sqrt(b_v*b_v + b_h*b_h)
488  If(b_p /= 0.0_fp)Then
489  cos_baziang = b_v / b_p
490  Else
491  cos_baziang = 0.0 ! not defined (take an arbitrary value)
492  Endif
493 
494  END SUBROUTINE compute_kb_angles
495 
496  Subroutine solar_za_az(latitude, &
497  longitude, &
498  julian_day, &
499  universal_time, &
500  ZA, &
501  Az)
502 !
503 !************************************************************************
504 !* *
505 !* Module Name: Solar_Az_Za *
506 !* *
507 !* Language: Fortran Library: *
508 !* Version.Rev: 1.0 22 Feb 91 Programmer: Kleespies *
509 !* 1.1 28 Feb 91 Kleespies *
510 !* Put equation of time into hour angle. *
511 !*
512 !* updated to f95, 4, September 2007, Y. Han *
513 !* *
514 !* Calling Seq: Call Solar_Az_ZA( *
515 !* & latitude, *
516 !* & longitude, *
517 !* & julian_day, *
518 !* & universal_time, *
519 !* & ZA, *
520 !* & Az) *
521 !* *
522 !* *
523 !* Description: Computes solar azimuth and zenith angle for a *
524 !* given place and time. Solar azimuth is the angle *
525 !* positive east of north form a point to the center of *
526 !* the solar disc. Solar zenith angle is the angle *
527 !* in degrees from zenith to the center of the solar disk. *
528 !* The algorithms are taken from "Introduction to Solar *
529 !* Radiation" by Muhamad Iqbal, Academic Press, 1983. *
530 !* Note that lat=0,lon=0, 21Mar 12z will not give sun *
531 !* overhead, because that is not the way things work. *
532 !* *
533 !* Input Args: R*4 Latitude, +NH, -SH degrees *
534 !* R*4 Longitude,-W, +E *
535 !* R*4 Julian_Day 1=Jan 1, 365=Dec31 (366 leap yr)*
536 !* R*4 Universal_Time 0.00-23.99,(GMT,Z time) *
537 !* *
538 !* Output Args: R*4 ZA Solar Zenith Angle *
539 !* R*4 AZ Solar Azmuth Angle *
540 !* *
541 !* Common Blks: none *
542 !* Include: none *
543 !* Externals: none *
544 !* Data Files: none *
545 !* *
546 !* Restrictions: Accurate to within .1 deg. *
547 !* No checking made to the validity of input *
548 !* parameters. *
549 !* Since solar zenith angle is a conic angle, *
550 !* it has no sign. *
551 !* No correction made for refraction. *
552 !* *
553 !* Error Codes: none *
554 !* *
555 !* Error Messages: *
556 !* *
557 !************************************************************************
558 !
559  implicit none
560 
561  real(fp), intent(in) :: latitude,longitude,julian_day,universal_time
562  real(fp), intent(out) :: ZA, Az
563 
564  real(fp) :: local_sun_time,solar_elevation,equation_of_time
565  real(fp) :: cosza,cosaz
566  real(fp) :: hour_angle,day_angle
567  real(fp) :: solar_declination
568  real(fp) :: rlatitude, rlongitude
569  real(fp), parameter :: DEGREES_TO_RADIANS = 3.141592653589793238462643_fp/180.0_fp
570  real(fp), parameter :: one_eighty_over_pi = 1.0/degrees_to_radians
571  real(fp), parameter :: threesixty_over_24 = 15.0
572  real(fp), parameter :: threesixty_over_365 = 0.98630137
573  real(fp), parameter :: min_declination = -23.433
574  real(fp), parameter :: day_offset = 10.0 ! original equation had this nine
575 
576 !* Compute day angle
577 
578  day_angle = threesixty_over_365*(julian_day-1.0)*degrees_to_radians
579 
580  rlatitude = latitude * degrees_to_radians
581  rlongitude = longitude * degrees_to_radians
582 
583 !* Compute equation of Time
584 
585  equation_of_time = &
586  ( 0.000075 &
587  + 0.001868*cos(day_angle) &
588  - 0.032077*sin(day_angle) &
589  - 0.014615*cos(2.*day_angle) &
590  - 0.040890*sin(2.*day_angle) )*229.18/60 ! in hours
591 
592 !* Compute local sun time
593  local_sun_time = universal_time &
594  + equation_of_time &
595  + longitude/threesixty_over_24
596 
597 !* Compute solar declination
598 
599  solar_declination = &
600  ( 0.006918 &
601  - 0.399912 * cos(day_angle) &
602  + 0.070257 * sin(day_angle) &
603  - 0.006758 * cos(2*day_angle) &
604  + 0.000907 * sin(2*day_angle) &
605  - 0.002697 * cos(3*day_angle) &
606  + 0.001480 * sin(3*day_angle) )
607 
608 !* Compute hour angle
609  hour_angle = threesixty_over_24*mod(local_sun_time+12.0_fp , 24.0_fp)
610 
611 !* Compute solar zenith angle
612  cosza = sin(rlatitude)*sin(solar_declination) &
613  + cos(rlatitude)*cos(solar_declination)*cos(hour_angle*degrees_to_radians)
614 
615  za = acos(cosza)*one_eighty_over_pi
616 
617 !* Compute solar azimuth angle
618  solar_elevation = 90.0 - za
619 
620  If(solar_elevation .eq. 90.0) Then
621  az = 180.0 ! handle arbitrary case
622  Else
623  cosaz = (sin(solar_elevation*degrees_to_radians)*sin(rlatitude) - &
624  sin(solar_declination)) / &
625  (cos(solar_elevation*degrees_to_radians)*cos(rlatitude))
626 
627  If(cosaz .lt. -1.0) cosaz = -1.0
628  If(cosaz .gt. 1.0) cosaz = 1.0
629 
630  az = acos(cosaz)*one_eighty_over_pi
631 
632 ! The above formula produces azimuth positive east, zero south.
633 ! We want positive east, zero north.
634 
635 !
636  If (az .ge. 0.0) Then
637  az = 180.0 - az
638  Else
639  az = -180.0 + az
640  EndIf
641 
642  If(hour_angle .lt. 180.0) az = - az
643  If(az .lt. 0) az = 360.0 + az
644 
645  EndIf
646 
647  end Subroutine solar_za_az
648 
649 End MODULE zeeman_utility
zeeman_utility::n_lon
integer, parameter n_lon
Definition: Zeeman_Utility.f90:45
message_handler::failure
integer, parameter, public failure
Definition: Message_Handler.f90:30
zeeman_utility
Definition: Zeeman_Utility.f90:2
type_kinds::fp
integer, parameter, public fp
Definition: Type_Kinds.f90:124
zeeman_utility::n_lat
integer, parameter n_lat
Definition: Zeeman_Utility.f90:45
zeeman_utility::load_bfield_lut
integer function, public load_bfield_lut(filename_LUT)
Definition: Zeeman_Utility.f90:73
zeeman_utility::compute_bfield_f1
subroutine compute_bfield_f1(latitude, longitude, Bx, By, Bz, Be)
Definition: Zeeman_Utility.f90:265
zeeman_utility::solar_za_az
subroutine solar_za_az(latitude, longitude, julian_day, universal_time, ZA, Az)
Definition: Zeeman_Utility.f90:502
file_utility::get_lun
integer function, public get_lun()
Definition: File_Utility.f90:85
message_handler
Definition: Message_Handler.f90:4
zeeman_utility::bfield
integer, dimension(3, n_lat, n_lon), save bfield
Definition: Zeeman_Utility.f90:46
zeeman_utility::compute_bfield_f2
subroutine compute_bfield_f2(latitude, longitude, sensor_zenang, sensor_aziang, Be, cos_bkang, cos_baziang)
Definition: Zeeman_Utility.f90:318
message_handler::display_message
recursive subroutine, public display_message(Routine_Name, Message, Error_State, Message_Log)
Definition: Message_Handler.f90:118
zeeman_utility::compute_kb_angles
subroutine, public compute_kb_angles(Bx, By, Bz, sensor_zenang, sensor_aziang, cos_bkang, cos_baziang)
Definition: Zeeman_Utility.f90:459
zeeman_utility::compute_bfield_f3
subroutine compute_bfield_f3(latitude, longitude, sensor_zenang, sensor_relative_aziang, julian_day, utc_time, Be, cos_bkang, cos_baziang)
Definition: Zeeman_Utility.f90:345
file_utility::file_exists
Definition: File_Utility.f90:51
crtm_parameters::degrees_to_radians
real(fp), parameter, public degrees_to_radians
Definition: CRTM_Parameters.f90:88
zeeman_utility::compute_bfield
Definition: Zeeman_Utility.f90:38
type_kinds
Definition: Type_Kinds.f90:12
bfield_component
real(fp) function bfield_component(comp, Bfield, w1_lat, w1_lon, idx_lat, idx_lon)
Definition: Zeeman_Utility.f90:294
message_handler::success
integer, parameter, public success
Definition: Message_Handler.f90:27
file_utility
Definition: File_Utility.f90:26