GFS Operational Physics Documentation  gsm/branches/DTC/phys-doc-all phys-doc-all R82971
radiation_astronomy.f
Go to the documentation of this file.
1 
4 
5 ! ========================================================== !!!!!
6 ! 'module_radiation_astronomy' description !!!!!
7 ! ========================================================== !!!!!
8 ! !
9 ! set up astronomy quantities for solar radiation calculations. !
10 ! !
11 ! in module 'module_radiation_astronomy', externally accessable !
12 ! subroutines are listed below: !
13 ! !
14 ! 'sol_init' -- initialization !
15 ! input: !
16 ! ( me ) !
17 ! output: !
18 ! ( none ) !
19 ! !
20 ! 'sol_update' -- update astronomy related quantities !
21 ! input: !
22 ! ( jdate,kyear,deltsw,deltim,lsol_chg, me ) !
23 ! output: !
24 ! ( slag,sdec,cdec,solcon ) !
25 ! !
26 ! 'coszmn' -- compute cosin of zenith angles !
27 ! input: !
28 ! ( xlon,sinlat,coslat,solhr,IM, me ) !
29 ! output: !
30 ! ( coszen,coszdg ) !
31 ! !
32 ! !
33 ! external modules referenced: !
34 ! 'module physparam' in 'physparam.f' !
35 ! 'module physcons' in 'physcons.f' !
36 ! !
37 ! program history log: !
38 ! - a collection of programs to track solar-earth position !
39 ! may 1977 --- ray orzol (gfdl) created program compjd to !
40 ! computes julian day and fraction from year,month,dayand,time!
41 ! jun 1977 --- robert white (gfdl) created program cdate to !
42 ! computes calendar month, day, year from julian day value. !
43 ! jul 1977 --- robert white (gfdl) created program solar to !
44 ! computes radius vector, declination and right ascension of !
45 ! sun, equation of time, hour angle, fractional daylight, and !
46 ! latitudinal mean zenith angle. !
47 ! fall 1988 --- hualu pan, updated to limit the iterations in !
48 ! newton method and also ccr reduced to avoid non-convergence.!
49 ! jul 1989 --- kenneth campana modified subr solar and created !
50 ! subr zenith for computations of effective mean cosz and !
51 ! daylight fraction. !
52 ! oct 1990 --- yu-tai hou created subr coszmn to replace !
53 ! the latitudinal mean cosz by time mean cosz at grid points. !
54 ! may 1998 --- mark iredell y2k compliance !
55 ! dec 2003 --- yu-tai hou combined compjd and fcstim and !
56 ! rewritten programs in fortran 90 compatable modular form. !
57 ! feb 2006 --- yu-tai hou add 11-yr solar constant cycle !
58 ! mar 2009 --- yu-tai hou modified solinit for climate !
59 ! hindcast situation responding to ic time. !
60 ! aug 2012 --- yu-tai hou modified coszmn to allows sw !
61 ! radiation calling interval less than 1 hr limit and linked !
62 ! model time step with numb of cosz evaluations. also changed !
63 ! the initialization subroutine 'solinit' into two parts: !
64 ! 'sol_init' is called at the start of run to set up module !
65 ! parameters; and 'sol_update' is called within the time !
66 ! loop to check and update data sets. !
67 ! nov 2012 --- yu-tai hou modified control parameters thru !
68 ! model 'physparam'. !
69 ! jan 2013 --- yu-tai hou modified to include new solar !
70 ! constant tables (noaa_a0, noaa_an, cmip_an, cmip_mn) !
71 ! !
72 !!!!! ========================================================== !!!!!
73 !!!!! end descriptions !!!!!
74 !!!!! ========================================================== !!!!!
75 
76 
77 
84 !========================================!
85  module module_radiation_astronomy !
86 !........................................!
87 !
88  use physparam, only : isolar, solar_file, kind_phys
89  use physcons, only : con_solr, con_solr_old, con_pi
90  use module_iounitdef, only : niradsf
91 !
92  implicit none
93 !
94  private
95 
96 ! --- version tag and last revision date
97  character(40), parameter :: &
98  & VTAGAST='NCEP-Radiation_astronomy v5.2 Jan 2013 '
99 ! & VTAGAST='NCEP-Radiation_astronomy v5.1 Nov 2012 '
100 
102  real (kind=kind_phys), parameter :: degrad = 180.0/con_pi
103  real (kind=kind_phys), parameter :: tpi = 2.0 * con_pi
104  real (kind=kind_phys), parameter :: hpi = 0.5 * con_pi
105  real (kind=kind_phys), parameter :: f12 = 12.0
106  real (kind=kind_phys), parameter :: f3600 = 3600.0
107  real (kind=kind_phys), parameter :: czlimt = 0.0001 ! ~ cos(89.99427)
108 ! real (kind=kind_phys), parameter :: pid12 = con_pi/f12 ! angle per hour
109  real (kind=kind_phys), parameter :: pid12 = (2.0*asin(1.0))/f12
110 
112  real (kind=kind_phys), public :: solc0 = con_solr
113  integer :: isolflg = 10
114  character(26) :: solar_fname = ' '
115 
117 
119  real (kind=kind_phys) :: sollag=0.0
121  real (kind=kind_phys) :: sindec=0.0
123  real (kind=kind_phys) :: cosdec=0.0
125  real (kind=kind_phys) :: anginc=0.0
127  real (kind=kind_phys) :: smon_sav(12)
128  data smon_sav(1:12) / 12*con_solr /
129 
131  integer :: iyr_sav =0
133  integer :: nstp =6
134 
135  public sol_init, sol_update, coszmn
136 
137 
138 ! =================
139  contains
140 ! =================
141 
145  subroutine sol_init &
146  & ( me ) ! --- inputs
147 ! --- outputs: ( none )
148 
149 ! =================================================================== !
150 ! !
151 ! initialize astronomy process, set up module constants. !
152 ! !
153 ! inputs: !
154 ! me - print message control flag !
155 ! !
156 ! outputs: (to module variable) !
157 ! ( none ) !
158 ! !
159 ! external module variable: (in physparam) !
160 ! isolar - = 0: use the old fixed solar constant in "physcon" !
161 ! =10: use the new fixed solar constant in "physcon" !
162 ! = 1: use noaa ann-mean tsi tbl abs-scale with cyc apprx !
163 ! = 2: use noaa ann-mean tsi tbl tim-scale with cyc apprx !
164 ! = 3: use cmip5 ann-mean tsi tbl tim-scale with cyc apprx!
165 ! = 4: use cmip5 mon-mean tsi tbl tim-scale with cyc apprx!
166 ! solar_file- external solar constant data table !
167 ! !
168 ! internal module variable: !
169 ! isolflg - internal solar constant scheme control flag !
170 ! solc0 - solar constant (w/m**2) !
171 ! solar_fname-file name for solar constant table assigned based on !
172 ! the scheme control flag, isolflg. !
173 ! !
174 ! usage: call sol_init !
175 ! !
176 ! subprograms called: none !
177 ! !
178 ! =================================================================== !
179 !
180  implicit none
181 
182 ! --- input:
183  integer, intent(in) :: me
184 
185 ! --- output: ( none )
186 
187 ! --- local:
188  logical :: file_exist
189 !
190 !===> ... begin here
191 !
192  if ( me == 0 ) print *, vtagast !print out version tag
193 
194 ! --- initialization
195  isolflg = isolar
196  solc0 = con_solr
197  solar_fname = solar_file
198  iyr_sav = 0
199  nstp = 6
200 
201  if ( isolar == 0 ) then
202  solc0 = con_solr_old
203  if ( me == 0 ) then
204  print *,' - Using old fixed solar constant =', solc0
205  endif
206  elseif ( isolar == 10 ) then
207  if ( me == 0 ) then
208  print *,' - Using new fixed solar constant =', solc0
209  endif
210  elseif ( isolar == 1 ) then ! noaa ann-mean tsi in absolute scale
211  solar_fname(15:26) = 'noaa_a0.txt'
212 
213  if ( me == 0 ) then
214  print *,' - Using NOAA annual mean TSI table in ABS scale', &
215  & ' with cycle approximation (old values)!'
216  endif
217 
218  inquire (file=solar_fname, exist=file_exist)
219  if ( .not. file_exist ) then
220  isolflg = 10
221 
222  if ( me == 0 ) then
223  print *,' Requested solar data file "',solar_fname, &
224  & '" not found!'
225  print *,' Using the default solar constant value =',solc0,&
226  & ' reset control flag isolflg=',isolflg
227  endif
228  endif
229  elseif ( isolar == 2 ) then ! noaa ann-mean tsi in tim scale
230  solar_fname(15:26) = 'noaa_an.txt'
231 
232  if ( me == 0 ) then
233  print *,' - Using NOAA annual mean TSI table in TIM scale', &
234  & ' with cycle approximation (new values)!'
235  endif
236 
237  inquire (file=solar_fname, exist=file_exist)
238  if ( .not. file_exist ) then
239  isolflg = 10
240 
241  if ( me == 0 ) then
242  print *,' Requested solar data file "',solar_fname, &
243  & '" not found!'
244  print *,' Using the default solar constant value =',solc0,&
245  & ' reset control flag isolflg=',isolflg
246  endif
247  endif
248  elseif ( isolar == 3 ) then ! cmip5 ann-mean tsi in tim scale
249  solar_fname(15:26) = 'cmip_an.txt'
250 
251  if ( me == 0 ) then
252  print *,' - Using CMIP5 annual mean TSI table in TIM scale', &
253  & ' with cycle approximation'
254  endif
255 
256  inquire (file=solar_fname, exist=file_exist)
257  if ( .not. file_exist ) then
258  isolflg = 10
259 
260  if ( me == 0 ) then
261  print *,' Requested solar data file "',solar_fname, &
262  & '" not found!'
263  print *,' Using the default solar constant value =',solc0,&
264  & ' reset control flag isolflg=',isolflg
265  endif
266  endif
267  elseif ( isolar == 4 ) then ! cmip5 mon-mean tsi in tim scale
268  solar_fname(15:26) = 'cmip_mn.txt'
269 
270  if ( me == 0 ) then
271  print *,' - Using CMIP5 monthly mean TSI table in TIM scale', &
272  & ' with cycle approximation'
273  endif
274 
275  inquire (file=solar_fname, exist=file_exist)
276  if ( .not. file_exist ) then
277  isolflg = 10
278 
279  if ( me == 0 ) then
280  print *,' Requested solar data file "',solar_fname, &
281  & '" not found!'
282  print *,' Using the default solar constant value =',solc0,&
283  & ' reset control flag isolflg=',isolflg
284  endif
285  endif
286  else ! selection error
287  isolflg = 10
288 
289  if ( me == 0 ) then
290  print *,' - !!! ERROR in selection of solar constant data', &
291  & ' source, ISOL =',isolar
292  print *,' Using the default solar constant value =',solc0, &
293  & ' reset control flag isolflg=',isolflg
294  endif
295  endif ! end if_isolar_block
296 !
297  return
298 !...................................
299  end subroutine sol_init
300 !-----------------------------------
301 
302 
319 !-----------------------------------
320  subroutine sol_update &
321  & ( jdate,kyear,deltsw,deltim,lsol_chg, me, & ! --- inputs
322  & slag, sdec, cdec, solcon & ! --- outputs
323  & )
325 ! =================================================================== !
326 ! !
327 ! sol_update computes solar parameters at forecast time !
328 ! !
329 ! inputs: !
330 ! jdate(8)- ncep absolute date and time at fcst time !
331 ! (yr, mon, day, t-zone, hr, min, sec, mil-sec) !
332 ! kyear - usually kyear=jdate(1). if not, it is for hindcast mode,!
333 ! and it is usually the init cond time and serves as the !
334 ! upper limit of data can be used. !
335 ! deltsw - time duration in seconds per sw calculation !
336 ! deltim - timestep in seconds !
337 ! lsol_chg- logical flags for change solar constant !
338 ! me - print message control flag !
339 ! !
340 ! outputs: !
341 ! slag - equation of time in radians !
342 ! sdec, cdec - sin and cos of the solar declination angle !
343 ! solcon - sun-earth distance adjusted solar constant (w/m2) !
344 ! !
345 ! !
346 ! module variable: !
347 ! solc0 - solar constant (w/m**2) not adjusted by earth-sun dist !
348 ! isolflg - solar constant control flag !
349 ! = 0: use the old fixed solar constant !
350 ! =10: use the new fixed solar constant !
351 ! = 1: use noaa ann-mean tsi tbl abs-scale with cycle apprx !
352 ! = 2: use noaa ann-mean tsi tbl tim-scale with cycle apprx !
353 ! = 3: use cmip5 ann-mean tsi tbl tim-scale with cycle apprx!
354 ! = 4: use cmip5 mon-mean tsi tbl tim-scale with cycle apprx!
355 ! solar_fname-external solar constant data table !
356 ! sindec - sine of the solar declination angle !
357 ! cosdec - cosine of the solar declination angle !
358 ! anginc - solar angle increment per iteration for cosz calc !
359 ! nstp - total number of zenith angle iterations !
360 ! smon_sav- saved monthly solar constants (isolflg=4 only) !
361 ! iyr_sav - saved year of data previously used !
362 ! !
363 ! usage: call sol_update !
364 ! !
365 ! subprograms called: solar, prtime !
366 ! !
367 ! external functions called: iw3jdn !
368 ! !
369 ! =================================================================== !
370 !
371  implicit none
372 
373 ! --- input:
374  integer, intent(in) :: jdate(:), kyear, me
375  logical, intent(in) :: lsol_chg
376 
377  real (kind=kind_phys), intent(in) :: deltsw, deltim
378 
379 ! --- output:
380  real (kind=kind_phys), intent(out) :: slag, sdec, cdec, solcon
381 
382 ! --- locals:
383  real (kind=kind_phys), parameter :: hrday = 1.0/24.0 ! frc day/hour
384  real (kind=kind_phys), parameter :: minday= 1.0/1440.0 ! frc day/minute
385  real (kind=kind_phys), parameter :: secday= 1.0/86400.0 ! frc day/second
386 
387  real (kind=kind_phys) :: smean, solc1, dtswh, smon(12)
388  real (kind=kind_phys) :: fjd, fjd1, dlt, r1, alp
389 
390  integer :: jd, jd1, iyear, imon, iday, ihr, imin, isec
391  integer :: iw3jdn
392  integer :: i, iyr, iyr1, iyr2, jyr, nn, nswr, icy1, icy2, icy
393 
394  logical :: file_exist
395  character :: cline*60
396 !
397 !===> ... begin here
398 !
399 ! --- ... forecast time
400  iyear = jdate(1)
401  imon = jdate(2)
402  iday = jdate(3)
403  ihr = jdate(5)
404  imin = jdate(6)
405  isec = jdate(7)
406 
407  if ( lsol_chg ) then ! get solar constant from data table
408 
409  if ( iyr_sav == iyear ) then ! same year, no new reading necessary
410  if ( isolflg==4 ) then
411  solc0 = smon_sav(imon)
412  endif
413  else ! need to read in new data
414  iyr_sav = iyear
415 
416 ! --- ... check to see if the solar constant data file existed
417 
418  inquire (file=solar_fname, exist=file_exist)
419  if ( .not. file_exist ) then
420  print *,' !!! ERROR! Can not find solar constant file!!!'
421  stop
422  else
423  iyr = iyear
424 
425  close(niradsf)
426  open (niradsf,file=solar_fname,form='formatted', &
427  & status='old')
428  rewind niradsf
429 
430  read (niradsf, * ) iyr1,iyr2,icy1,icy2,smean,cline(1:60)
431 ! read (NIRADSF, 24) iyr1,iyr2,icy1,icy2,smean,cline
432 ! 24 format(4i5,f8.2,a60)
433 
434  if ( me == 0 ) then
435  print *,' Updating solar constant with cycle approx'
436  print *,' Opened solar constant data file: ',solar_fname
437 !check print *, iyr1, iyr2, icy1, icy2, smean, cline
438  endif
439 
440 ! --- ... check if there is a upper year limit put on the data table
441 
442 ! if ( iyear /= kyear ) then
443 ! icy = icy1 - iyr1 + 1 ! range of the earlest cycle in data table
444 ! if ( kyear-iyr1 < icy ) then ! need data range at least icy years
445  ! to perform cycle approximation
446 ! if ( me == 0 ) then
447 ! print *,' *** the requested year',iyear,' and upper',&
448 ! & 'limit',kyear,' do not fit the range of data ', &
449 ! & 'table of iyr1, iyr2 =',iyr1,iyr2
450 ! print *,' USE FIXED SOLAR CONSTANT=',con_solr
451 ! endif
452 ! solc0 = con_solr
453 ! isolflg = 10
454 
455 ! elseif ( kyear < iyr2 ) then
456 
457 ! --- ... because the usage limit put on the historical data table,
458 ! skip those unused data records at first
459 
460 ! i = iyr2
461 ! Lab_dowhile0 : do while ( i > kyear )
462 ! read (NIRADSF,26) jyr, solc1
463 ! 26 format(i4,f10.4)
464 ! read (NIRADSF,*) jyr, solc1
465 ! i = i - 1
466 ! enddo Lab_dowhile0
467 
468 ! iyr2 = kyear ! next record will serve the upper limit
469 
470 ! endif ! end if_kyear_block
471 ! endif ! end if_iyear_block
472 
473 ! --- ... checking the cycle range
474 
475  if ( iyr < iyr1 ) then
476  icy = icy1 - iyr1 + 1 ! range of the earlest cycle in data table
477  lab_dowhile1 : do while ( iyr < iyr1 )
478  iyr = iyr + icy
479  enddo lab_dowhile1
480 
481  if ( me == 0 ) then
482  print *,' *** Year',iyear,' out of table range!', &
483  & iyr1, iyr2
484  print *,' Using the closest-cycle year (',iyr,')'
485  endif
486  elseif ( iyr > iyr2 ) then
487  icy = iyr2 - icy2 + 1 ! range of the latest cycle in data table
488  lab_dowhile2 : do while ( iyr > iyr2 )
489  iyr = iyr - icy
490  enddo lab_dowhile2
491 
492  if ( me == 0 ) then
493  print *,' *** Year',iyear,' out of table range!', &
494  & iyr1, iyr2
495  print *,' Using the closest-cycle year (',iyr,')'
496  endif
497  endif
498 
499 ! --- ... locate the right record for the year of data
500 
501  if ( isolflg < 4 ) then ! use annual mean data tables
502  i = iyr2
503  lab_dowhile3 : do while ( i >= iyr1 )
504 ! read (NIRADSF,26) jyr, solc1
505 ! 26 format(i4,f10.4)
506  read (niradsf,*) jyr, solc1
507 
508  if ( i == iyr .and. iyr == jyr ) then
509  solc0 = smean + solc1
510 
511  if (me == 0) then
512  print *,' CHECK: Solar constant data used for year',&
513  & iyr, solc1, solc0
514  endif
515  exit lab_dowhile3
516  else
517 !check if(me == 0) print *,' Skip solar const data for yr',i
518  i = i - 1
519  endif
520  enddo lab_dowhile3
521  elseif ( isolflg == 4 ) then ! use monthly mean data tables
522  i = iyr2
523  lab_dowhile4 : do while ( i >= iyr1 )
524 ! read (NIRADSF,26) jyr, smon(:)
525 ! 26 format(i4,12f10.4)
526  read (niradsf,*) jyr, smon(1:12)
527 
528  if ( i == iyr .and. iyr == jyr ) then
529  do nn = 1, 12
530  smon_sav(nn) = smean + smon(nn)
531  enddo
532  solc0 = smean + smon(imon)
533 
534  if (me == 0) then
535  print *,' CHECK: Solar constant data used for year',&
536  & iyr,' and month',imon
537  endif
538  exit lab_dowhile4
539  else
540 !check if(me == 0) print *,' Skip solar const data for yr',i
541  i = i - 1
542  endif
543  enddo lab_dowhile4
544  endif ! end if_isolflg_block
545 
546  close ( niradsf )
547  endif ! end if_file_exist_block
548 
549  endif ! end if_iyr_sav_block
550  endif ! end if_lsol_chg_block
551 
552 ! --- ... calculate forecast julian day and fraction of julian day
553 
554  jd1 = iw3jdn(iyear,imon,iday)
555 
556 ! --- ... unlike in normal applications, where day starts from 0 hr,
557 ! in astronomy applications, day stats from noon.
558 
559  if (ihr < 12) then
560  jd1 = jd1 - 1
561  fjd1= 0.5 + float(ihr)*hrday + float(imin)*minday &
562  & + float(isec)*secday
563  else
564  fjd1= float(ihr - 12)*hrday + float(imin)*minday &
565  & + float(isec)*secday
566  endif
567 
568  fjd1 = fjd1 + jd1
569 
570  jd = int(fjd1)
571  fjd = fjd1 - jd
572 
574  call solar &
575 ! --- inputs:
576  & ( jd, fjd, &
577 ! --- outputs:
578  & r1, dlt, alp &
579  & )
580 
581 ! --- ... calculate sun-earth distance adjustment factor appropriate to date
582  solcon = solc0 / (r1*r1)
583 
584  slag = sollag
585  sdec = sindec
586  cdec = cosdec
587 
588 ! --- ... diagnostic print out
589 
590  if (me == 0) then
591 
593  call prtime &
594 ! --- inputs:
595  & ( jd, fjd, dlt, alp, r1, solcon )
596 ! --- outputs: ( none )
597 
598  endif
599 
600 ! --- ... setting up calculation parameters used by subr coszmn
601 
602  nswr = nint(deltsw / deltim) ! number of mdl t-step per sw call
603  dtswh = deltsw / f3600 ! time length in hours
604 
605  if ( deltsw >= f3600 ) then ! for longer sw call interval
606  nn = max(6, min(12, nint(f3600/deltim) )) ! num of calc per hour
607  nstp = nint(dtswh) * nn + 1 ! num of calc per sw call
608  else ! for shorter sw sw call interval
609  nstp = max(2, min(20, nswr)) + 1
610 ! nn = nint( float(nstp-1)/dtswh )
611  endif
612 
613  anginc = pid12 * dtswh / float(nstp-1) ! solar angle inc during each calc step
614 
615  if ( me == 0 ) then
616  print *,' for cosz calculations: nswr,deltim,deltsw,dtswh =', &
617  & nswr,deltim,deltsw,dtswh,' anginc,nstp =',anginc,nstp
618  endif
619 
620 ! if (me == 0) print*,'in sol_update completed sr solar'
621 !
622  return
623 !...................................
624  end subroutine sol_update
625 !-----------------------------------
626 !! @}
627 
628 
631 !-----------------------------------
632  subroutine solar &
633  & ( jd, fjd, & ! --- inputs
634  & r1, dlt, alp & ! --- outputs
635  & )
637 ! =================================================================== !
638 ! !
639 ! solar computes radius vector, declination and right ascension of !
640 ! sun, and equation of time. !
641 ! !
642 ! inputs: !
643 ! jd - julian day !
644 ! fjd - fraction of the julian day !
645 ! !
646 ! outputs: !
647 ! r1 - earth-sun radius vector !
648 ! dlt - declination of sun in radians !
649 ! alp - right ascension of sun in radians !
650 ! !
651 ! module variables: !
652 ! sollag - equation of time in radians !
653 ! sindec - sine of declination angle !
654 ! cosdec - cosine of declination angle !
655 ! !
656 ! usage: call solar !
657 ! !
658 ! external subroutines called: none !
659 ! !
660 ! =================================================================== !
661 !
662  implicit none
663 
664 ! --- inputs:
665  real (kind=kind_phys), intent(in) :: fjd
666  integer, intent(in) :: jd
667 
668 ! --- outputs:
669  real (kind=kind_phys), intent(out) :: r1, dlt, alp
670 
671 ! --- locals:
672  real (kind=kind_phys), parameter :: cyear = 365.25 ! days of year
673  real (kind=kind_phys), parameter :: ccr = 1.3e-6 ! iteration limit
674  real (kind=kind_phys), parameter :: tpp = 1.55 ! days between epoch and
675  ! perihelion passage of 1900
676  real (kind=kind_phys), parameter :: svt6 = 78.035 ! days between perihelion passage
677  ! and march equinox of 1900
678  integer, parameter :: jdor = 2415020 ! jd of epoch which is january
679  ! 0, 1900 at 12 hours ut
680 
681  real (kind=kind_phys) :: dat, t1, year, tyear, ec, angin, ador, &
682  & deleqn, sni, tini, er, qq, e1, ep, cd, eq, date, em, &
683  & cr, w1, tst, sun
684 
685  integer :: jdoe, iter
686 
687 !===> ... begin here
688 
689 ! --- ... computes time in julian centuries after epoch
690 
691  t1 = float(jd - jdor) / 36525.0
692 
693 ! --- ... computes length of anomalistic and tropical years (minus 365 days)
694 
695  year = 0.25964134e0 + 0.304e-5 * t1
696  tyear= 0.24219879e0 - 0.614e-5 * t1
697 
698 ! --- ... computes orbit eccentricity and angle of earth's inclination from t
699 
700  ec = 0.01675104e0 - (0.418e-4 + 0.126e-6 * t1) * t1
701  angin= 23.452294e0 - (0.0130125e0 + 0.164e-5 * t1) * t1
702 
703  ador = jdor
704  jdoe = ador + (svt6 * cyear) / (year - tyear)
705 
706 ! --- ... deleqn is updated svt6 for current date
707 
708  deleqn= float(jdoe - jd) * (year - tyear) / cyear
709  year = year + 365.0
710  sni = sin( angin / degrad )
711  tini = 1.0 / tan( angin / degrad )
712  er = sqrt( (1.0 + ec) / (1.0 - ec) )
713  qq = deleqn * tpi / year
714 
715 ! --- ... determine true anomaly at equinox
716 
717  e1 = 1.0
718  cd = 1.0
719  iter = 0
720 
721  lab_do_1 : do while ( cd > ccr )
722 
723  ep = e1 - (e1 - ec*sin(e1) - qq) / (1.0 - ec*cos(e1))
724  cd = abs(e1 - ep)
725  e1 = ep
726  iter = iter + 1
727 
728  if (iter > 10) then
729  write(6,*) ' ITERATION COUNT FOR LOOP 32 =', iter
730  write(6,*) ' E, EP, CD =', e1, ep, cd
731  exit lab_do_1
732  endif
733 
734  enddo lab_do_1
735 
736  eq = 2.0 * atan( er * tan( 0.5*e1 ) )
737 
738 ! --- ... date is days since last perihelion passage
739 
740  dat = float(jd - jdor) - tpp + fjd
741  date = mod(dat, year)
742 
743 ! --- ... solve orbit equations by newton's method
744 
745  em = tpi * date / year
746  e1 = 1.0
747  cr = 1.0
748  iter = 0
749 
750  lab_do_2 : do while ( cr > ccr )
751 
752  ep = e1 - (e1 - ec*sin(e1) - em) / (1.0 - ec*cos(e1))
753  cr = abs(e1 - ep)
754  e1 = ep
755  iter = iter + 1
756 
757  if (iter > 10) then
758  write(6,*) ' ITERATION COUNT FOR LOOP 31 =', iter
759  exit lab_do_2
760  endif
761 
762  enddo lab_do_2
763 
764  w1 = 2.0 * atan( er * tan( 0.5*e1 ) )
765 
766  r1 = 1.0 - ec*cos(e1)
767 
768  sindec = sni * sin(w1 - eq)
769  cosdec = sqrt( 1.0 - sindec*sindec )
770 
771  dlt = asin( sindec )
772  alp = asin( tan(dlt)*tini )
773 
774  tst = cos( w1 - eq )
775  if (tst < 0.0) alp = con_pi - alp
776  if (alp < 0.0) alp = alp + tpi
777 
778  sun = tpi * (date - deleqn) / year
779  if (sun < 0.0) sun = sun + tpi
780  sollag = sun - alp - 0.03255e0
781 !
782  return
783 !...................................
784  end subroutine solar
785 !-----------------------------------
786 
787 
800 !-----------------------------------
801  subroutine coszmn &
802  & ( xlon,sinlat,coslat,solhr, im, me, & ! --- inputs
803  & coszen, coszdg & ! --- outputs
804  & )
806 ! =================================================================== !
807 ! !
808 ! coszmn computes mean cos solar zenith angle over sw calling interval !
809 ! !
810 ! inputs: !
811 ! xlon (IM) - grids' longitudes in radians, work both on zonal !
812 ! 0->2pi and -pi->+pi arrangements !
813 ! sinlat(IM) - sine of the corresponding latitudes !
814 ! coslat(IM) - cosine of the corresponding latitudes !
815 ! solhr - time after 00z in hours !
816 ! IM - num of grids in horizontal dimension !
817 ! me - print message control flag !
818 ! !
819 ! outputs: !
820 ! coszen(IM) - average of cosz for daytime only in sw call interval
821 ! coszdg(IM) - average of cosz over entire sw call interval !
822 ! !
823 ! module variables: !
824 ! sollag - equation of time !
825 ! sindec - sine of the solar declination angle !
826 ! cosdec - cosine of the solar declination angle !
827 ! anginc - solar angle increment per iteration for cosz calc !
828 ! nstp - total number of zenith angle iterations !
829 ! !
830 ! usage: call comzmn !
831 ! !
832 ! external subroutines called: none !
833 ! !
834 ! =================================================================== !
835 !
836  implicit none
837 
838 ! --- inputs:
839  integer, intent(in) :: IM, me
840 
841  real (kind=kind_phys), intent(in) :: sinlat(:), coslat(:), &
842  & xlon(:), solhr
843 
844 ! --- outputs:
845  real (kind=kind_phys), intent(out) :: coszen(:), coszdg(:)
846 
847 ! --- locals:
848  real (kind=kind_phys) :: coszn, cns, ss, cc, solang, rstp
849 
850  integer :: istsun(im), i, it, j, lat
851 
852 !===> ... begin here
853 
854  solang = pid12 * (solhr - f12) ! solar angle at present time
855  rstp = 1.0 / float(nstp)
856 
857  do i = 1, im
858  coszen(i) = 0.0
859  istsun(i) = 0
860  enddo
861 
862  do it = 1, nstp
863  cns = solang + float(it-1)*anginc + sollag
864 
865  do i = 1, im
866  ss = sinlat(i) * sindec
867  cc = coslat(i) * cosdec
868 
869  coszn = ss + cc * cos(cns + xlon(i))
870  coszen(i) = coszen(i) + max(0.0, coszn)
871  if (coszn > czlimt) istsun(i) = istsun(i) + 1
872  enddo
873  enddo
874 
875 ! --- ... compute time averages
876 
877  do i = 1, im
878  coszdg(i) = coszen(i) * rstp
879  if (istsun(i) > 0) coszen(i) = coszen(i) / istsun(i)
880  enddo
881 !
882  return
883 !...................................
884  end subroutine coszmn
885 !-----------------------------------
886 
887 
890 !-----------------------------------
891  subroutine prtime &
892  & ( jd, fjd, dlt, alp, r1, solc & ! --- inputs
893  & ) ! --- outputs: ( none )
895 ! =================================================================== !
896 ! !
897 ! prtime prints out forecast date, time, and astronomy quantities. !
898 ! !
899 ! inputs: !
900 ! jd - forecast julian day !
901 ! fjd - forecast fraction of julian day !
902 ! dlt - declination angle of sun in radians !
903 ! alp - right ascension of sun in radians !
904 ! r1 - earth-sun radius vector in meter !
905 ! solc - solar constant in w/m^2 !
906 ! !
907 ! outputs: ( none ) !
908 ! !
909 ! module variables: !
910 ! sollag - equation of time in radians !
911 ! !
912 ! usage: call prtime !
913 ! !
914 ! external subroutines called: w3fs26 !
915 ! !
916 ! =================================================================== !
917 !
918  implicit none
919 
920 ! --- inputs:
921  integer, intent(in) :: jd
922 
923  real (kind=kind_phys), intent(in) :: fjd, dlt, alp, r1, solc
924 
925 ! --- outputs: ( none )
926 
927 ! --- locals:
928  real (kind=kind_phys), parameter :: sixty = 60.0
929 
930  character(LEN=1), parameter :: sign = '-'
931  character(LEN=1), parameter :: sigb = ' '
932 
933  character(LEN=1) :: dsig
934  character(LEN=4) :: month(12)
935 
936  data month / 'JAN.','FEB.','MAR.','APR.','MAY ','JUNE', &
937  & 'JULY','AUG.','SEP.','OCT.','NOV ','DEC.' /
938 
939  integer :: iday, imon, iyear, ihr, ltd, ltm, &
940  & ihalp, iyy, jda, mfjd, idaywk, idayyr
941  real (kind=kind_phys) :: xmin, dltd, dltm, dlts, halp, ymin, &
942  & asec, eqt, eqsec
943 
944 !===> ... begin here
945 
946 ! --- ... get forecast hour and minute from fraction of julian day
947 
948  if (fjd >= 0.5) then
949  jda = jd + 1
950  mfjd= nint( fjd*1440.0 )
951  ihr = mfjd / 60 - 12
952  xmin= float(mfjd) - (ihr + 12)*sixty
953  else
954  jda = jd
955  mfjd= nint( fjd*1440.0 )
956  ihr = mfjd / 60 + 12
957  xmin= float(mfjd) - (ihr - 12)*sixty
958  endif
959 
960 ! --- ... get forecast year, month, and day from julian day
961 
962  call w3fs26(jda, iyear,imon,iday, idaywk,idayyr)
963 
964 ! -- ... compute solar parameters
965 
966  dltd = degrad * dlt
967  ltd = dltd
968  dltm = sixty * (abs(dltd) - abs(float(ltd)))
969  ltm = dltm
970  dlts = sixty * (dltm - float(ltm))
971 
972  if ((dltd < 0.0) .and. (ltd == 0.0)) then
973  dsig = sign
974  else
975  dsig = sigb
976  endif
977 
978  halp = 6.0 * alp / hpi
979  ihalp= halp
980  ymin = abs(halp - float(ihalp)) * sixty
981  iyy = ymin
982  asec = (ymin - float(iyy)) * sixty
983 
984  eqt = 228.55735 * sollag
985  eqsec= sixty * eqt
986 
987  print 101, iday, month(imon), iyear, ihr, xmin, jd, fjd
988  101 format('0 FORECAST DATE',9x,i3,a5,i6,' AT',i3,' HRS',f6.2,' MINS'/&
989  & ' JULIAN DAY',12x,i8,2x,'PLUS',f11.6)
990 
991  print 102, r1, halp, ihalp, iyy, asec
992  102 format(' RADIUS VECTOR',9x,f10.7/' RIGHT ASCENSION OF SUN', &
993  & f12.7,' HRS, OR',i4,' HRS',i4,' MINS',f6.1,' SECS')
994 
995  print 103, dltd, dsig, ltd, ltm, dlts, eqt, eqsec, sollag, solc
996  103 format(' DECLINATION OF THE SUN',f12.7,' DEGS, OR ',a1,i3, &
997  & ' DEGS',i4,' MINS',f6.1,' SECS'/' EQUATION OF TIME',6x, &
998  & f12.7,' MINS, OR',f10.2,' SECS, OR',f9.6,' RADIANS'/ &
999  & ' SOLAR CONSTANT',8x,f12.7,' (DISTANCE AJUSTED)'//)
1000 
1001 !
1002  return
1003 !...................................
1004  end subroutine prtime
1005 !-----------------------------------
1006 
1007 !
1008 !...........................................!
1009  end module module_radiation_astronomy !
1010 !===========================================!
real(kind=kind_phys), parameter con_pi
pi
Definition: physcons.f:48
integer iyr_sav
saved year of data used
real(kind=kind_phys) anginc
solar angle increment per interation of cosz calc
real(kind=kind_phys), parameter con_solr
solar constant ( )-nasa-sorce Tim(2008)
Definition: physcons.f:68
integer nstp
total number of zenith angle iterations
real(kind=kind_phys) sollag
equation of time
subroutine, public coszmn(xlon, sinlat, coslat, solhr, IM, me, coszen, coszdg )
This subroutine computes mean cos solar zenith angle over SW calling interval.
subroutine solar(jd, fjd, r1, dlt, alp )
This subroutine computes radius vector, declination and right ascension of sun, and equation of time...
real(kind=kind_phys) cosdec
cosine of the solar declination angle
real(kind=kind_phys) sindec
sine of the solar declination angle
subroutine, public sol_init(me)
This subroutine initializes astronomy process, and set up module constants.
integer, save isolar
solar constant scheme control flag =0:fixed value=1366.0 (old standard) =10:fixed value=1360...
Definition: physparam.f:153
real(kind=kind_phys), dimension(12) smon_sav
saved monthly solar constants (isolflg=4 only)
subroutine prtime(jd, fjd, dlt, alp, r1, solc )
This subroutine prints out forecast date, time, and astronomy quantities.
subroutine, public sol_update(jdate, kyear, deltsw, deltim, lsol_chg, me, slag, sdec, cdec, solcon )
This subroutine computes solar parameters at forecast time.
real(kind=kind_phys), parameter con_solr_old
solar constant ( )-liu(2002)
Definition: physcons.f:66
character, save solar_file
external solar constant data table,solarconstant_noaa_a0.txt
Definition: physparam.f:156