Radiation Scheme in CCPP
Functions/Subroutines | Variables
module_radiation_astronomy Module Reference

This module sets up astronomy quantities for solar radiation calculations.

Functions/Subroutines

subroutine, public sol_update
 This subroutine computes solar parameters at forecast time. More...
 
subroutine solar
 This subroutine computes radius vector, declination and right ascension of sun, and equation of time. More...
 
subroutine, public coszmn
 This subroutine computes mean cos solar zenith angle over SW calling interval. More...
 
subroutine prtime
 This subroutine prints out forecast date, time, and astronomy quantities. More...
 
subroutine, public sol_init
 This subroutine initializes astronomy process, and set up module constants. More...
 

Variables

character(40), parameter vtagast ='NCEP-Radiation_astronomy v5.2 Jan 2013 '
 
real(kind=kind_phys), parameter degrad = 180.0/con_pi
 
real(kind=kind_phys), parameter tpi = 2.0 * con_pi
 
real(kind=kind_phys), parameter hpi = 0.5 * con_pi
 
real(kind=kind_phys), parameter f12 = 12.0
 
real(kind=kind_phys), parameter f3600 = 3600.0
 
real(kind=kind_phys), parameter czlimt = 0.0001
 
real(kind=kind_phys), parameter pid12 = (2.0*asin(1.0))/f12
 
real(kind=kind_phys), public solc0 = con_solr
 
integer isolflg = 10
 
character(26) solar_fname = ' '
 
real(kind=kind_phys) sollag =0.0
 
real(kind=kind_phys) sindec =0.0
 
real(kind=kind_phys) cosdec =0.0
 
real(kind=kind_phys) anginc =0.0
 
real(kind=kind_phys), dimension(12) smon_sav
 
integer iyr_sav =0
 
integer nstp =6
 

Function/Subroutine Documentation

subroutine, public module_radiation_astronomy::coszmn ( )
Parameters
[in]xlonreal, (IM), grids' longitudes in radians, work both on zonal, 0->2pi and -pi->+pi arrangements
[in]sinlatreal, (IM), sine of the corresponding latitudes
[in]coslatreal, (IM), cosine of the corresponding latitudes
[in]solhrreal, time after 00z in hours
[in]IMinteger, num of grids in horizontal dimension
[in]meinteger, print message control flag
[out]coszenreal, (IM), average of cosz for daytime only in sw call interval
[out]coszdgreal, (IM), average of cosz over entire sw call interval

Definition at line 788 of file radiation_astronomy.f.

References anginc, cosdec, czlimt, f12, nstp, pid12, sindec, and sollag.

Referenced by module_radiation_driver::grrad().

788 ! --- inputs:
789  & ( xlon,sinlat,coslat,solhr, im, me, &
790 ! --- outputs:
791  & coszen, coszdg &
792  & )
793 
794 ! =================================================================== !
795 ! !
796 ! coszmn computes mean cos solar zenith angle over sw calling interval !
797 ! !
798 ! inputs: !
799 ! xlon (IM) - grids' longitudes in radians, work both on zonal !
800 ! 0->2pi and -pi->+pi arrangements !
801 ! sinlat(IM) - sine of the corresponding latitudes !
802 ! coslat(IM) - cosine of the corresponding latitudes !
803 ! solhr - time after 00z in hours !
804 ! IM - num of grids in horizontal dimension !
805 ! me - print message control flag !
806 ! !
807 ! outputs: !
808 ! coszen(IM) - average of cosz for daytime only in sw call interval
809 ! coszdg(IM) - average of cosz over entire sw call interval !
810 ! !
811 ! module variables: !
812 ! sollag - equation of time !
813 ! sindec - sine of the solar declination angle !
814 ! cosdec - cosine of the solar declination angle !
815 ! anginc - solar angle increment per iteration for cosz calc !
816 ! nstp - total number of zenith angle iterations !
817 ! !
818 ! usage: call comzmn !
819 ! !
820 ! external subroutines called: none !
821 ! !
822 ! =================================================================== !
823 !
824  implicit none
825 
826 ! --- inputs:
827  integer, intent(in) :: im, me
828 
829  real (kind=kind_phys), intent(in) :: sinlat(:), coslat(:), &
830  & xlon(:), solhr
831 
832 ! --- outputs:
833  real (kind=kind_phys), intent(out) :: coszen(:), coszdg(:)
834 
835 ! --- locals:
836  real (kind=kind_phys) :: coszn, cns, ss, cc, solang, rstp
837 
838  integer :: istsun(im), i, it, j, lat
839 
840 !===> ... begin here
841 
842  solang = pid12 * (solhr - f12) ! solar angle at present time
843  rstp = 1.0 / float(nstp)
844 
845  do i = 1, im
846  coszen(i) = 0.0
847  istsun(i) = 0
848  enddo
849 
850  do it = 1, nstp
851  cns = solang + float(it-1)*anginc + sollag
852 
853  do i = 1, im
854  ss = sinlat(i) * sindec
855  cc = coslat(i) * cosdec
856 
857  coszn = ss + cc * cos(cns + xlon(i))
858  coszen(i) = coszen(i) + max(0.0, coszn)
859  if (coszn > czlimt) istsun(i) = istsun(i) + 1
860  enddo
861  enddo
862 
863 ! --- ... compute time averages
864 
865  do i = 1, im
866  coszdg(i) = coszen(i) * rstp
867  if (istsun(i) > 0) coszen(i) = coszen(i) / istsun(i)
868  enddo
869 !
870  return
871 !...................................

Here is the caller graph for this function:

subroutine module_radiation_astronomy::prtime ( )
private

Definition at line 879 of file radiation_astronomy.f.

References degrad, f12, hpi, and sollag.

Referenced by sol_update().

879 ! --- inputs:
880  & ( jd, fjd, dlt, alp, r1, solc &
881 ! --- outputs: ( none )
882  & )
883 
884 ! =================================================================== !
885 ! !
886 ! prtime prints out forecast date, time, and astronomy quantities. !
887 ! !
888 ! inputs: !
889 ! jd - forecast julian day !
890 ! fjd - forecast fraction of julian day !
891 ! dlt - declination angle of sun in radians !
892 ! alp - right ascension of sun in radians !
893 ! r1 - earth-sun radius vector in meter !
894 ! solc - solar constant in w/m^2 !
895 ! !
896 ! outputs: ( none ) !
897 ! !
898 ! module variables: !
899 ! sollag - equation of time in radians !
900 ! !
901 ! usage: call prtime !
902 ! !
903 ! external subroutines called: w3fs26 !
904 ! !
905 ! =================================================================== !
906 !
907  implicit none
908 
909 ! --- inputs:
910  integer, intent(in) :: jd
911 
912  real (kind=kind_phys), intent(in) :: fjd, dlt, alp, r1, solc
913 
914 ! --- outputs: ( none )
915 
916 ! --- locals:
917  real (kind=kind_phys), parameter :: sixty = 60.0
918 
919  character(LEN=1), parameter :: sign = '-'
920  character(LEN=1), parameter :: sigb = ' '
921 
922  character(LEN=1) :: dsig
923  character(LEN=4) :: month(12)
924 
925  data month / 'JAN.','FEB.','MAR.','APR.','MAY ','JUNE', &
926  & 'JULY','AUG.','SEP.','OCT.','NOV ','DEC.' /
927 
928  integer :: iday, imon, iyear, ihr, ltd, ltm, &
929  & ihalp, iyy, jda, mfjd, idaywk, idayyr
930  real (kind=kind_phys) :: xmin, dltd, dltm, dlts, halp, ymin, &
931  & asec, eqt, eqsec
932 
933 !===> ... begin here
934 
935 ! --- ... get forecast hour and minute from fraction of julian day
936 
937  if (fjd >= 0.5) then
938  jda = jd + 1
939  mfjd= nint( fjd*1440.0 )
940  ihr = mfjd / 60 - 12
941  xmin= float(mfjd) - (ihr + 12)*sixty
942  else
943  jda = jd
944  mfjd= nint( fjd*1440.0 )
945  ihr = mfjd / 60 + 12
946  xmin= float(mfjd) - (ihr - 12)*sixty
947  endif
948 
949 ! --- ... get forecast year, month, and day from julian day
950 
951  call w3fs26(jda, iyear,imon,iday, idaywk,idayyr)
952 
953 ! -- ... compute solar parameters
954 
955  dltd = degrad * dlt
956  ltd = dltd
957  dltm = sixty * (abs(dltd) - abs(float(ltd)))
958  ltm = dltm
959  dlts = sixty * (dltm - float(ltm))
960 
961  if ((dltd < 0.0) .and. (ltd == 0.0)) then
962  dsig = sign
963  else
964  dsig = sigb
965  endif
966 
967  halp = 6.0 * alp / hpi
968  ihalp= halp
969  ymin = abs(halp - float(ihalp)) * sixty
970  iyy = ymin
971  asec = (ymin - float(iyy)) * sixty
972 
973  eqt = 228.55735 * sollag
974  eqsec= sixty * eqt
975 
976  print 101, iday, month(imon), iyear, ihr, xmin, jd, fjd
977  101 format('0 FORECAST DATE',9x,i3,a5,i6,' AT',i3,' HRS',f6.2,' MINS'/&
978  & ' JULIAN DAY',12x,i8,2x,'PLUS',f11.6)
979 
980  print 102, r1, halp, ihalp, iyy, asec
981  102 format(' RADIUS VECTOR',9x,f10.7/' RIGHT ASCENSION OF SUN', &
982  & f12.7,' HRS, OR',i4,' HRS',i4,' MINS',f6.1,' SECS')
983 
984  print 103, dltd, dsig, ltd, ltm, dlts, eqt, eqsec, sollag, solc
985  103 format(' DECLINATION OF THE SUN',f12.7,' DEGS, OR ',a1,i3, &
986  & ' DEGS',i4,' MINS',f6.1,' SECS'/' EQUATION OF TIME',6x, &
987  & f12.7,' MINS, OR',f10.2,' SECS, OR',f9.6,' RADIANS'/ &
988  & ' SOLAR CONSTANT',8x,f12.7,' (DISTANCE AJUSTED)'//)
989 
990 !
991  return
992 !...................................

Here is the caller graph for this function:

subroutine, public module_radiation_astronomy::sol_init ( )
Parameters
[in]meprint message control flag
[out]NONE

External Module Variables


physparam::isolar
- = 0: use the old fixed solar constant in "physcon"
- =10: use the new fixed solar constant in "physcon"
- = 1: use noaa ann-mean tsi tbl abs-scale with cyc apprx
- = 2: use noaa ann-mean tsi tbl tim-scale with cyc apprx
- = 3: use cmip5 ann-mean tsi tbl tim-scale with cyc apprx
- = 4: use cmip5 mon-mean tsi tbl tim-scale with cyc apprx

Definition at line 139 of file radiation_astronomy.f.

References physcons::con_solr, physcons::con_solr_old, physparam::isolar, isolflg, iyr_sav, nstp, physparam::solar_file, solar_fname, solc0, and vtagast.

Referenced by module_radiation_driver::radinit().

139 ! --- inputs:
140  & ( me )
141 ! --- outputs: ( none )
142 
143 ! =================================================================== !
144 ! !
145 ! initialize astronomy process, set up module constants. !
146 ! !
147 ! inputs: !
148 ! me - print message control flag !
149 ! !
150 ! outputs: (to module variable) !
151 ! ( none ) !
152 ! !
153 ! external module variable: (in physparam) !
154 ! isolar - = 0: use the old fixed solar constant in "physcon" !
155 ! =10: use the new fixed solar constant in "physcon" !
156 ! = 1: use noaa ann-mean tsi tbl abs-scale with cyc apprx !
157 ! = 2: use noaa ann-mean tsi tbl tim-scale with cyc apprx !
158 ! = 3: use cmip5 ann-mean tsi tbl tim-scale with cyc apprx!
159 ! = 4: use cmip5 mon-mean tsi tbl tim-scale with cyc apprx!
160 ! solar_file- external solar constant data table !
161 ! !
162 ! internal module variable: !
163 ! isolflg - internal solar constant scheme control flag !
164 ! solc0 - solar constant (w/m**2) !
165 ! solar_fname-file name for solar constant table assigned based on !
166 ! the scheme control flag, isolflg. !
167 ! !
168 ! usage: call sol_init !
169 ! !
170 ! subprograms called: none !
171 ! !
172 ! =================================================================== !
173 !
174  implicit none
175 
176 ! --- input:
177  integer, intent(in) :: me
178 
179 ! --- output: ( none )
180 
181 ! --- local:
182  logical :: file_exist
183 !
184 !===> ... begin here
185 !
186  if ( me == 0 ) print *, vtagast !print out version tag
187 
188 ! --- initialization
189  isolflg = isolar
190  solc0 = con_solr
191  solar_fname = solar_file
192  iyr_sav = 0
193  nstp = 6
194 
195  if ( isolar == 0 ) then
196  solc0 = con_solr_old
197  if ( me == 0 ) then
198  print *,' - Using old fixed solar constant =', solc0
199  endif
200  elseif ( isolar == 10 ) then
201  if ( me == 0 ) then
202  print *,' - Using new fixed solar constant =', solc0
203  endif
204  elseif ( isolar == 1 ) then ! noaa ann-mean tsi in absolute scale
205  solar_fname(15:26) = 'noaa_a0.txt'
206 
207  if ( me == 0 ) then
208  print *,' - Using NOAA annual mean TSI table in ABS scale', &
209  & ' with cycle approximation (old values)!'
210  endif
211 
212  inquire (file=solar_fname, exist=file_exist)
213  if ( .not. file_exist ) then
214  isolflg = 10
215 
216  if ( me == 0 ) then
217  print *,' Requested solar data file "',solar_fname, &
218  & '" not found!'
219  print *,' Using the default solar constant value =',solc0,&
220  & ' reset control flag isolflg=',isolflg
221  endif
222  endif
223  elseif ( isolar == 2 ) then ! noaa ann-mean tsi in tim scale
224  solar_fname(15:26) = 'noaa_an.txt'
225 
226  if ( me == 0 ) then
227  print *,' - Using NOAA annual mean TSI table in TIM scale', &
228  & ' with cycle approximation (new values)!'
229  endif
230 
231  inquire (file=solar_fname, exist=file_exist)
232  if ( .not. file_exist ) then
233  isolflg = 10
234 
235  if ( me == 0 ) then
236  print *,' Requested solar data file "',solar_fname, &
237  & '" not found!'
238  print *,' Using the default solar constant value =',solc0,&
239  & ' reset control flag isolflg=',isolflg
240  endif
241  endif
242  elseif ( isolar == 3 ) then ! cmip5 ann-mean tsi in tim scale
243  solar_fname(15:26) = 'cmip_an.txt'
244 
245  if ( me == 0 ) then
246  print *,' - Using CMIP5 annual mean TSI table in TIM scale', &
247  & ' with cycle approximation'
248  endif
249 
250  inquire (file=solar_fname, exist=file_exist)
251  if ( .not. file_exist ) then
252  isolflg = 10
253 
254  if ( me == 0 ) then
255  print *,' Requested solar data file "',solar_fname, &
256  & '" not found!'
257  print *,' Using the default solar constant value =',solc0,&
258  & ' reset control flag isolflg=',isolflg
259  endif
260  endif
261  elseif ( isolar == 4 ) then ! cmip5 mon-mean tsi in tim scale
262  solar_fname(15:26) = 'cmip_mn.txt'
263 
264  if ( me == 0 ) then
265  print *,' - Using CMIP5 monthly mean TSI table in TIM scale', &
266  & ' with cycle approximation'
267  endif
268 
269  inquire (file=solar_fname, exist=file_exist)
270  if ( .not. file_exist ) then
271  isolflg = 10
272 
273  if ( me == 0 ) then
274  print *,' Requested solar data file "',solar_fname, &
275  & '" not found!'
276  print *,' Using the default solar constant value =',solc0,&
277  & ' reset control flag isolflg=',isolflg
278  endif
279  endif
280  else ! selection error
281  isolflg = 10
282 
283  if ( me == 0 ) then
284  print *,' - !!! ERROR in selection of solar constant data', &
285  & ' source, ISOL =',isolar
286  print *,' Using the default solar constant value =',solc0, &
287  & ' reset control flag isolflg=',isolflg
288  endif
289  endif ! end if_isolar_block
290 !
291  return
292 !...................................

Here is the caller graph for this function:

subroutine, public module_radiation_astronomy::sol_update ( )
Parameters
[in]jdateinteger, 1, ncep absolute date and time at fcst time (yr, mon, day, t-zone, hr, min, sec, mil-sec)
[in]kyearinteger, 1, usually kyear=jdate(1). if not, it is for hindcast mode, and it is usually the init cond time and serves as the upper limit of data can be used.
[in]deltswreal, 1, time duration in seconds per sw calculation
[in]deltimreal, 1, timestep in seconds
[in]lsol_chglogical flags for change solar constant
[in]meinteger, print message control flag
[out]slagreal, equation of time in radians
[out]sdec,cdecreal, sin and cos of the solar declination angle
[out]solconreal, sun-earth distance adjusted solar constant \((w/m^2)\)

Definition at line 310 of file radiation_astronomy.f.

References anginc, cosdec, f3600, isolflg, iyr_sav, nstp, pid12, prtime(), sindec, smon_sav, solar(), solar_fname, solc0, and sollag.

Referenced by module_radiation_driver::radupdate().

310 ! --- inputs:
311  & ( jdate,kyear,deltsw,deltim,lsol_chg, me, &
312 ! --- outputs:
313  & slag, sdec, cdec, solcon &
314  & )
315 
316 ! =================================================================== !
317 ! !
318 ! sol_update computes solar parameters at forecast time !
319 ! !
320 ! inputs: !
321 ! jdate(8)- ncep absolute date and time at fcst time !
322 ! (yr, mon, day, t-zone, hr, min, sec, mil-sec) !
323 ! kyear - usually kyear=jdate(1). if not, it is for hindcast mode,!
324 ! and it is usually the init cond time and serves as the !
325 ! upper limit of data can be used. !
326 ! deltsw - time duration in seconds per sw calculation !
327 ! deltim - timestep in seconds !
328 ! lsol_chg- logical flags for change solar constant !
329 ! me - print message control flag !
330 ! !
331 ! outputs: !
332 ! slag - equation of time in radians !
333 ! sdec, cdec - sin and cos of the solar declination angle !
334 ! solcon - sun-earth distance adjusted solar constant (w/m2) !
335 ! !
336 ! !
337 ! module variable: !
338 ! solc0 - solar constant (w/m**2) not adjusted by earth-sun dist !
339 ! isolflg - solar constant control flag !
340 ! = 0: use the old fixed solar constant !
341 ! =10: use the new fixed solar constant !
342 ! = 1: use noaa ann-mean tsi tbl abs-scale with cycle apprx !
343 ! = 2: use noaa ann-mean tsi tbl tim-scale with cycle apprx !
344 ! = 3: use cmip5 ann-mean tsi tbl tim-scale with cycle apprx!
345 ! = 4: use cmip5 mon-mean tsi tbl tim-scale with cycle apprx!
346 ! solar_fname-external solar constant data table !
347 ! sindec - sine of the solar declination angle !
348 ! cosdec - cosine of the solar declination angle !
349 ! anginc - solar angle increment per iteration for cosz calc !
350 ! nstp - total number of zenith angle iterations !
351 ! smon_sav- saved monthly solar constants (isolflg=4 only) !
352 ! iyr_sav - saved year of data previously used !
353 ! !
354 ! usage: call sol_update !
355 ! !
356 ! subprograms called: solar, prtime !
357 ! !
358 ! external functions called: iw3jdn !
359 ! !
360 ! =================================================================== !
361 !
362  implicit none
363 
364 ! --- input:
365  integer, intent(in) :: jdate(:), kyear, me
366  logical, intent(in) :: lsol_chg
367 
368  real (kind=kind_phys), intent(in) :: deltsw, deltim
369 
370 ! --- output:
371  real (kind=kind_phys), intent(out) :: slag, sdec, cdec, solcon
372 
373 ! --- locals:
374  real (kind=kind_phys), parameter :: hrday = 1.0/24.0 ! frc day/hour
375  real (kind=kind_phys), parameter :: minday= 1.0/1440.0 ! frc day/minute
376  real (kind=kind_phys), parameter :: secday= 1.0/86400.0 ! frc day/second
377 
378  real (kind=kind_phys) :: smean, solc1, dtswh, smon(12)
379  real (kind=kind_phys) :: fjd, fjd1, dlt, r1, alp
380 
381  integer :: jd, jd1, iyear, imon, iday, ihr, imin, isec
382  integer :: iw3jdn
383  integer :: i, iyr, iyr1, iyr2, jyr, nn, nswr, icy1, icy2, icy
384 
385  logical :: file_exist
386  character :: cline*60
387 !
388 !===> ... begin here
389 !
390 ! --- ... forecast time
391  iyear = jdate(1)
392  imon = jdate(2)
393  iday = jdate(3)
394  ihr = jdate(5)
395  imin = jdate(6)
396  isec = jdate(7)
397 
398  if ( lsol_chg ) then ! get solar constant from data table
399 
400  if ( iyr_sav == iyear ) then ! same year, no new reading necessary
401  if ( isolflg==4 ) then
402  solc0 = smon_sav(imon)
403  endif
404  else ! need to read in new data
405  iyr_sav = iyear
406 
407 ! --- ... check to see if the solar constant data file existed
408 
409  inquire (file=solar_fname, exist=file_exist)
410  if ( .not. file_exist ) then
411  print *,' !!! ERROR! Can not find solar constant file!!!'
412  stop
413  else
414  iyr = iyear
415 
416  close(niradsf)
417  open (niradsf,file=solar_fname,form='formatted', &
418  & status='old')
419  rewind niradsf
420 
421  read (niradsf, * ) iyr1,iyr2,icy1,icy2,smean,cline(1:60)
422 ! read (NIRADSF, 24) iyr1,iyr2,icy1,icy2,smean,cline
423 ! 24 format(4i5,f8.2,a60)
424 
425  if ( me == 0 ) then
426  print *,' Updating solar constant with cycle approx'
427  print *,' Opened solar constant data file: ',solar_fname
428 !check print *, iyr1, iyr2, icy1, icy2, smean, cline
429  endif
430 
431 ! --- ... check if there is a upper year limit put on the data table
432 
433 ! if ( iyear /= kyear ) then
434 ! icy = icy1 - iyr1 + 1 ! range of the earlest cycle in data table
435 ! if ( kyear-iyr1 < icy ) then ! need data range at least icy years
436  ! to perform cycle approximation
437 ! if ( me == 0 ) then
438 ! print *,' *** the requested year',iyear,' and upper',&
439 ! & 'limit',kyear,' do not fit the range of data ', &
440 ! & 'table of iyr1, iyr2 =',iyr1,iyr2
441 ! print *,' USE FIXED SOLAR CONSTANT=',con_solr
442 ! endif
443 ! solc0 = con_solr
444 ! isolflg = 10
445 
446 ! elseif ( kyear < iyr2 ) then
447 
448 ! --- ... because the usage limit put on the historical data table,
449 ! skip those unused data records at first
450 
451 ! i = iyr2
452 ! Lab_dowhile0 : do while ( i > kyear )
453 ! read (NIRADSF,26) jyr, solc1
454 ! 26 format(i4,f10.4)
455 ! read (NIRADSF,*) jyr, solc1
456 ! i = i - 1
457 ! enddo Lab_dowhile0
458 
459 ! iyr2 = kyear ! next record will serve the upper limit
460 
461 ! endif ! end if_kyear_block
462 ! endif ! end if_iyear_block
463 
464 ! --- ... checking the cycle range
465 
466  if ( iyr < iyr1 ) then
467  icy = icy1 - iyr1 + 1 ! range of the earlest cycle in data table
468  lab_dowhile1 : do while ( iyr < iyr1 )
469  iyr = iyr + icy
470  enddo lab_dowhile1
471 
472  if ( me == 0 ) then
473  print *,' *** Year',iyear,' out of table range!', &
474  & iyr1, iyr2
475  print *,' Using the closest-cycle year (',iyr,')'
476  endif
477  elseif ( iyr > iyr2 ) then
478  icy = iyr2 - icy2 + 1 ! range of the latest cycle in data table
479  lab_dowhile2 : do while ( iyr > iyr2 )
480  iyr = iyr - icy
481  enddo lab_dowhile2
482 
483  if ( me == 0 ) then
484  print *,' *** Year',iyear,' out of table range!', &
485  & iyr1, iyr2
486  print *,' Using the closest-cycle year (',iyr,')'
487  endif
488  endif
489 
490 ! --- ... locate the right record for the year of data
491 
492  if ( isolflg < 4 ) then ! use annual mean data tables
493  i = iyr2
494  lab_dowhile3 : do while ( i >= iyr1 )
495 ! read (NIRADSF,26) jyr, solc1
496 ! 26 format(i4,f10.4)
497  read (niradsf,*) jyr, solc1
498 
499  if ( i == iyr .and. iyr == jyr ) then
500  solc0 = smean + solc1
501 
502  if (me == 0) then
503  print *,' CHECK: Solar constant data used for year',&
504  & iyr, solc1, solc0
505  endif
506  exit lab_dowhile3
507  else
508 !check if(me == 0) print *,' Skip solar const data for yr',i
509  i = i - 1
510  endif
511  enddo lab_dowhile3
512  elseif ( isolflg == 4 ) then ! use monthly mean data tables
513  i = iyr2
514  lab_dowhile4 : do while ( i >= iyr1 )
515 ! read (NIRADSF,26) jyr, smon(:)
516 ! 26 format(i4,12f10.4)
517  read (niradsf,*) jyr, smon(1:12)
518 
519  if ( i == iyr .and. iyr == jyr ) then
520  do nn = 1, 12
521  smon_sav(nn) = smean + smon(nn)
522  enddo
523  solc0 = smean + smon(imon)
524 
525  if (me == 0) then
526  print *,' CHECK: Solar constant data used for year',&
527  & iyr,' and month',imon
528  endif
529  exit lab_dowhile4
530  else
531 !check if(me == 0) print *,' Skip solar const data for yr',i
532  i = i - 1
533  endif
534  enddo lab_dowhile4
535  endif ! end if_isolflg_block
536 
537  close ( niradsf )
538  endif ! end if_file_exist_block
539 
540  endif ! end if_iyr_sav_block
541  endif ! end if_lsol_chg_block
542 
543 ! --- ... calculate forecast julian day and fraction of julian day
544 
545  jd1 = iw3jdn(iyear,imon,iday)
546 
547 ! --- ... unlike in normal applications, where day starts from 0 hr,
548 ! in astronomy applications, day stats from noon.
549 
550  if (ihr < 12) then
551  jd1 = jd1 - 1
552  fjd1= 0.5 + float(ihr)*hrday + float(imin)*minday &
553  & + float(isec)*secday
554  else
555  fjd1= float(ihr - 12)*hrday + float(imin)*minday &
556  & + float(isec)*secday
557  endif
558 
559  fjd1 = fjd1 + jd1
560 
561  jd = int(fjd1)
562  fjd = fjd1 - jd
563 
564  call solar &
565 ! --- inputs:
566  & ( jd, fjd, &
567 ! --- outputs:
568  & r1, dlt, alp &
569  & )
570 
571 ! --- ... calculate sun-earth distance adjustment factor appropriate to date
572  solcon = solc0 / (r1*r1)
573 
574  slag = sollag
575  sdec = sindec
576  cdec = cosdec
577 
578 ! --- ... diagnostic print out
579 
580  if (me == 0) then
581 
582  call prtime &
583 ! --- inputs:
584  & ( jd, fjd, dlt, alp, r1, solcon )
585 ! --- outputs: ( none )
586 
587  endif
588 
589 ! --- ... setting up calculation parameters used by subr coszmn
590 
591  nswr = nint(deltsw / deltim) ! number of mdl t-step per sw call
592  dtswh = deltsw / f3600 ! time length in hours
593 
594  if ( deltsw >= f3600 ) then ! for longer sw call interval
595  nn = max(6, min(12, nint(f3600/deltim) )) ! num of calc per hour
596  nstp = nint(dtswh) * nn + 1 ! num of calc per sw call
597  else ! for shorter sw sw call interval
598  nstp = max(2, min(20, nswr)) + 1
599 ! nn = nint( float(nstp-1)/dtswh )
600  endif
601 
602  anginc = pid12 * dtswh / float(nstp-1) ! solar angle inc during each calc step
603 
604  if ( me == 0 ) then
605  print *,' for cosz calculations: nswr,deltim,deltsw,dtswh =', &
606  & nswr,deltim,deltsw,dtswh,' anginc,nstp =',anginc,nstp
607  endif
608 
609 ! if (me == 0) print*,'in sol_update completed sr solar'
610 !
611  return
612 !...................................

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine module_radiation_astronomy::solar ( )
private

Definition at line 620 of file radiation_astronomy.f.

References physcons::con_pi, cosdec, degrad, sindec, sollag, and tpi.

Referenced by sol_update().

620 ! --- inputs:
621  & ( jd, fjd, &
622 ! --- outputs:
623  & r1, dlt, alp &
624  & )
625 
626 ! =================================================================== !
627 ! !
628 ! solar computes radius vector, declination and right ascension of !
629 ! sun, and equation of time. !
630 ! !
631 ! inputs: !
632 ! jd - julian day !
633 ! fjd - fraction of the julian day !
634 ! !
635 ! outputs: !
636 ! r1 - earth-sun radius vector !
637 ! dlt - declination of sun in radians !
638 ! alp - right ascension of sun in radians !
639 ! !
640 ! module variables: !
641 ! sollag - equation of time in radians !
642 ! sindec - sine of declination angle !
643 ! cosdec - cosine of declination angle !
644 ! !
645 ! usage: call solar !
646 ! !
647 ! external subroutines called: none !
648 ! !
649 ! =================================================================== !
650 !
651  implicit none
652 
653 ! --- inputs:
654  real (kind=kind_phys), intent(in) :: fjd
655  integer, intent(in) :: jd
656 
657 ! --- outputs:
658  real (kind=kind_phys), intent(out) :: r1, dlt, alp
659 
660 ! --- locals:
661  real (kind=kind_phys), parameter :: cyear = 365.25 ! days of year
662  real (kind=kind_phys), parameter :: ccr = 1.3e-6 ! iteration limit
663  real (kind=kind_phys), parameter :: tpp = 1.55 ! days between epoch and
664  ! perihelion passage of 1900
665  real (kind=kind_phys), parameter :: svt6 = 78.035 ! days between perihelion passage
666  ! and march equinox of 1900
667  integer, parameter :: jdor = 2415020 ! jd of epoch which is january
668  ! 0, 1900 at 12 hours ut
669 
670  real (kind=kind_phys) :: dat, t1, year, tyear, ec, angin, ador, &
671  & deleqn, sni, tini, er, qq, e1, ep, cd, eq, date, em, &
672  & cr, w1, tst, sun
673 
674  integer :: jdoe, iter
675 
676 !===> ... begin here
677 
678 ! --- ... computes time in julian centuries after epoch
679 
680  t1 = float(jd - jdor) / 36525.0
681 
682 ! --- ... computes length of anomalistic and tropical years (minus 365 days)
683 
684  year = 0.25964134e0 + 0.304e-5 * t1
685  tyear= 0.24219879e0 - 0.614e-5 * t1
686 
687 ! --- ... computes orbit eccentricity and angle of earth's inclination from t
688 
689  ec = 0.01675104e0 - (0.418e-4 + 0.126e-6 * t1) * t1
690  angin= 23.452294e0 - (0.0130125e0 + 0.164e-5 * t1) * t1
691 
692  ador = jdor
693  jdoe = ador + (svt6 * cyear) / (year - tyear)
694 
695 ! --- ... deleqn is updated svt6 for current date
696 
697  deleqn= float(jdoe - jd) * (year - tyear) / cyear
698  year = year + 365.0
699  sni = sin( angin / degrad )
700  tini = 1.0 / tan( angin / degrad )
701  er = sqrt( (1.0 + ec) / (1.0 - ec) )
702  qq = deleqn * tpi / year
703 
704 ! --- ... determine true anomaly at equinox
705 
706  e1 = 1.0
707  cd = 1.0
708  iter = 0
709 
710  lab_do_1 : do while ( cd > ccr )
711 
712  ep = e1 - (e1 - ec*sin(e1) - qq) / (1.0 - ec*cos(e1))
713  cd = abs(e1 - ep)
714  e1 = ep
715  iter = iter + 1
716 
717  if (iter > 10) then
718  write(6,*) ' ITERATION COUNT FOR LOOP 32 =', iter
719  write(6,*) ' E, EP, CD =', e1, ep, cd
720  exit lab_do_1
721  endif
722 
723  enddo lab_do_1
724 
725  eq = 2.0 * atan( er * tan( 0.5*e1 ) )
726 
727 ! --- ... date is days since last perihelion passage
728 
729  dat = float(jd - jdor) - tpp + fjd
730  date = mod(dat, year)
731 
732 ! --- ... solve orbit equations by newton's method
733 
734  em = tpi * date / year
735  e1 = 1.0
736  cr = 1.0
737  iter = 0
738 
739  lab_do_2 : do while ( cr > ccr )
740 
741  ep = e1 - (e1 - ec*sin(e1) - em) / (1.0 - ec*cos(e1))
742  cr = abs(e1 - ep)
743  e1 = ep
744  iter = iter + 1
745 
746  if (iter > 10) then
747  write(6,*) ' ITERATION COUNT FOR LOOP 31 =', iter
748  exit lab_do_2
749  endif
750 
751  enddo lab_do_2
752 
753  w1 = 2.0 * atan( er * tan( 0.5*e1 ) )
754 
755  r1 = 1.0 - ec*cos(e1)
756 
757  sindec = sni * sin(w1 - eq)
758  cosdec = sqrt( 1.0 - sindec*sindec )
759 
760  dlt = asin( sindec )
761  alp = asin( tan(dlt)*tini )
762 
763  tst = cos( w1 - eq )
764  if (tst < 0.0) alp = con_pi - alp
765  if (alp < 0.0) alp = alp + tpi
766 
767  sun = tpi * (date - deleqn) / year
768  if (sun < 0.0) sun = sun + tpi
769  sollag = sun - alp - 0.03255e0
770 !
771  return
772 !...................................

Here is the caller graph for this function:

Variable Documentation

real (kind=kind_phys) module_radiation_astronomy::anginc =0.0
private

Definition at line 110 of file radiation_astronomy.f.

Referenced by coszmn(), and sol_update().

110  real (kind=kind_phys) :: anginc=0.0 ! solar angle incrmt per iteration for cosz calc
real (kind=kind_phys) module_radiation_astronomy::cosdec =0.0
private

Definition at line 109 of file radiation_astronomy.f.

Referenced by coszmn(), sol_update(), and solar().

109  real (kind=kind_phys) :: cosdec=0.0 ! cosine of the solar declination angle
real (kind=kind_phys), parameter module_radiation_astronomy::czlimt = 0.0001
private

Definition at line 97 of file radiation_astronomy.f.

Referenced by coszmn().

97  real (kind=kind_phys), parameter :: czlimt = 0.0001 ! ~ cos(89.99427)
real (kind=kind_phys), parameter module_radiation_astronomy::degrad = 180.0/con_pi
private

Definition at line 92 of file radiation_astronomy.f.

Referenced by prtime(), and solar().

92  real (kind=kind_phys), parameter :: degrad = 180.0/con_pi
real (kind=kind_phys), parameter module_radiation_astronomy::f12 = 12.0
private

Definition at line 95 of file radiation_astronomy.f.

Referenced by coszmn(), and prtime().

95  real (kind=kind_phys), parameter :: f12 = 12.0
real (kind=kind_phys), parameter module_radiation_astronomy::f3600 = 3600.0
private

Definition at line 96 of file radiation_astronomy.f.

Referenced by sol_update().

96  real (kind=kind_phys), parameter :: f3600 = 3600.0
real (kind=kind_phys), parameter module_radiation_astronomy::hpi = 0.5 * con_pi
private

Definition at line 94 of file radiation_astronomy.f.

Referenced by prtime().

94  real (kind=kind_phys), parameter :: hpi = 0.5 * con_pi
integer module_radiation_astronomy::isolflg = 10
private

Definition at line 103 of file radiation_astronomy.f.

Referenced by sol_init(), and sol_update().

103  integer :: isolflg = 10
integer module_radiation_astronomy::iyr_sav =0
private

Definition at line 114 of file radiation_astronomy.f.

Referenced by sol_init(), and sol_update().

114  integer :: iyr_sav =0 ! saved year of data used
integer module_radiation_astronomy::nstp =6
private

Definition at line 115 of file radiation_astronomy.f.

Referenced by coszmn(), sol_init(), and sol_update().

115  integer :: nstp =6 ! total number of zenith angle iterations
real (kind=kind_phys), parameter module_radiation_astronomy::pid12 = (2.0*asin(1.0))/f12
private

Definition at line 99 of file radiation_astronomy.f.

Referenced by coszmn(), and sol_update().

99  real (kind=kind_phys), parameter :: pid12 = (2.0*asin(1.0))/f12
real (kind=kind_phys) module_radiation_astronomy::sindec =0.0
private

Definition at line 108 of file radiation_astronomy.f.

Referenced by coszmn(), sol_update(), and solar().

108  real (kind=kind_phys) :: sindec=0.0 ! sineof the solar declination angle
real (kind=kind_phys), dimension(12) module_radiation_astronomy::smon_sav
private

Definition at line 111 of file radiation_astronomy.f.

Referenced by sol_update().

111  real (kind=kind_phys) :: smon_sav(12) ! saved monthly solar constants (isolflg=4 only)
character(26) module_radiation_astronomy::solar_fname = ' '
private

Definition at line 104 of file radiation_astronomy.f.

Referenced by sol_init(), and sol_update().

104  character(26) :: solar_fname = ' '
real (kind=kind_phys), public module_radiation_astronomy::solc0 = con_solr

Definition at line 102 of file radiation_astronomy.f.

Referenced by sol_init(), and sol_update().

102  real (kind=kind_phys), public :: solc0 = con_solr
real (kind=kind_phys) module_radiation_astronomy::sollag =0.0
private

Definition at line 107 of file radiation_astronomy.f.

Referenced by coszmn(), prtime(), sol_update(), and solar().

107  real (kind=kind_phys) :: sollag=0.0 ! equation of time
real (kind=kind_phys), parameter module_radiation_astronomy::tpi = 2.0 * con_pi
private

Definition at line 93 of file radiation_astronomy.f.

Referenced by solar().

93  real (kind=kind_phys), parameter :: tpi = 2.0 * con_pi
character(40), parameter module_radiation_astronomy::vtagast ='NCEP-Radiation_astronomy v5.2 Jan 2013 '
private

Definition at line 87 of file radiation_astronomy.f.

Referenced by sol_init().

87  character(40), parameter :: &
88  & VTAGAST='NCEP-Radiation_astronomy v5.2 Jan 2013 '