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