GFS Operational Physics Documentation  gsm/branches/DTC/phys-doc-all phys-doc-all R82971
radiation_gases.f
Go to the documentation of this file.
1 
6 
7 ! ========================================================== !!!!!
8 ! 'module_radiation_gases' description !!!!!
9 ! ========================================================== !!!!!
10 ! !
11 ! set up ozone climatological profiles and other constant gas !
12 ! profiles, such as co2, ch4, n2o, o2, and those of cfc gases. All !
13 ! data are entered as mixing ratio by volume, except ozone which is !
14 ! mass mixing ratio (g/g). !
15 ! !
16 ! in the module, the externally callabe subroutines are : !
17 ! !
18 ! 'gas_init' -- initialization !
19 ! input: !
20 ! ( me ) !
21 ! output: !
22 ! ( none ) !
23 ! !
24 ! 'gas_update' -- read in data and update with time !
25 ! input: !
26 ! ( iyear, imon, iday, ihour, loz1st, ldoco2, me ) !
27 ! output: !
28 ! ( none ) !
29 ! !
30 ! 'getozn' -- setup climatological ozone profile !
31 ! input: !
32 ! ( prslk,xlat, !
33 ! IMAX, LM ) !
34 ! output: !
35 ! ( o3mmr ) !
36 ! !
37 ! 'getgases' -- setup constant gas profiles for LW and SW !
38 ! input: !
39 ! ( plvl, xlon, xlat, !
40 ! IMAX, LMAX ) !
41 ! output: !
42 ! ( gasdat ) !
43 ! !
44 ! external modules referenced: !
45 ! 'module machine' in 'machine.f' !
46 ! 'module funcphys' in 'funcphys.f' !
47 ! 'module physcons' in 'physcons.f !
48 ! 'module module_iounitdef' in 'iounitdef.f' !
49 ! !
50 ! unit used for radiative active gases: !
51 ! ozone : mass mixing ratio (g/g) !
52 ! co2 : volume mixing ratio (p/p) !
53 ! n2o : volume mixing ratio (p/p) !
54 ! ch4 : volume mixing ratio (p/p) !
55 ! o2 : volume mixing ratio (p/p) !
56 ! co : volume mixing ratio (p/p) !
57 ! cfc11 : volume mixing ratio (p/p) !
58 ! cfc12 : volume mixing ratio (p/p) !
59 ! cfc22 : volume mixing ratio (p/p) !
60 ! ccl4 : volume mixing ratio (p/p) !
61 ! cfc113: volume mixing ratio (p/p) !
62 ! !
63 ! !
64 ! program history: !
65 ! may 2003 - y-t hou create rad_module.f that collectively !
66 ! combines several radiation computation supporting !
67 ! programs into fortran 90 module structure (gases !
68 ! and aerosols, etc.) !
69 ! apr 2004 - y-t hou modified to add astronomy and surface !
70 ! module components. !
71 ! feb 2005 - y-t hou rewrite the component modules into !
72 ! separate individule modules for thier corresponding !
73 ! tasks. here as radiation_gases.f !
74 ! mar 2006 - y-t hou add initialization subroutine to co2 and !
75 ! other gases. historical 2-d co2 data are added. !
76 ! sep 2008 - y-t hou add parameter ictm to control the input !
77 ! data time at the model initial condition. !
78 ! oct 2008 - y-t hou modify the initialization code to add the !
79 ! option of superimposing climatology seasonal cycle !
80 ! to the initial condition data (currently co2 only) !
81 ! nov 2008 - y-t hou fix bugs in superimposing climatology !
82 ! seasonal cycle calculations !
83 ! aug 2011 - y-t hou fix a bug in subr getgases doing vertical !
84 ! co2 mapping. (for iflip=0 case, not affact opr). !
85 ! aug 2012 - y-t hou modified subr getozn. moved the if-first !
86 ! block to subr gas_init to ensure threading safe in !
87 ! climatology ozone applications. (not affect gfs) !
88 ! also changed the initialization subr into two parts:!
89 ! 'gas_init' is called at the start of run to set up !
90 ! module parameters; and 'gas_update' is called within!
91 ! the time loop to check and update data sets. defined!
92 ! the climatology ozone parameters k1oz,k2oz,facoz as !
93 ! module variables and are set in subr 'gas_update' !
94 ! nov 2012 - y-t hou modified control parameters thru module !
95 ! 'physparam'. !
96 ! jan 2013 - z. janjic/y. hou modified ilon (longitude index) !
97 ! computing formula in subroutine getgases to work !
98 ! properly for models with either of 0->360 or !
99 ! -180->180 zonal grid directions. !
100 ! !
101 ! !
102 !!!!! ========================================================== !!!!!
103 !!!!! end descriptions !!!!!
104 !!!!! ========================================================== !!!!!
105 
106 
115 !========================================!
116  module module_radiation_gases !
117 !........................................!
118 !
119  use physparam, only : ico2flg, ictmflg, ioznflg, ivflip, &
120  & co2dat_file, co2gbl_file, &
121  & co2usr_file, co2cyc_file, &
122  & kind_phys, kind_io4
123  use funcphys, only : fpkapx
124  use physcons, only : con_pi
125  use ozne_def, only : jmr => latsozc, loz => levozc, &
126  & blte => blatc, dlte=> dphiozc, &
127  & timeozc => timeozc
128  use module_iounitdef, only : nio3clm, nico2cn
129 !
130  implicit none
131 !
132  private
133 
134 ! --- version tag and last revision date
135  character(40), parameter :: &
136  & VTAGGAS='NCEP-Radiation_gases v5.1 Nov 2012 '
137 ! & VTAGGAS='NCEP-Radiation_gases v5.0 Aug 2012 '
138 
140 
142  integer, parameter, public :: nf_vgas = 10
144  integer, parameter :: imxco2 = 24
146  integer, parameter :: jmxco2 = 12
148  integer, parameter :: minyear = 1957
149 
151  real (kind=kind_phys), parameter :: resco2=15.0
153  real (kind=kind_phys), parameter :: raddeg=180.0/con_pi
155  real (kind=kind_phys), parameter :: prsco2=788.0
157  real (kind=kind_phys), parameter :: hfpi =0.5*con_pi
158 
160 
161  real (kind=kind_phys), parameter :: co2vmr_def = 350.0e-6
162  real (kind=kind_phys), parameter :: n2ovmr_def = 0.31e-6
163  real (kind=kind_phys), parameter :: ch4vmr_def = 1.50e-6
164  real (kind=kind_phys), parameter :: o2vmr_def = 0.209
165  real (kind=kind_phys), parameter :: covmr_def = 1.50e-8
167  real (kind=kind_phys), parameter :: f11vmr_def = 3.520e-10
169  real (kind=kind_phys), parameter :: f12vmr_def = 6.358e-10
171  real (kind=kind_phys), parameter :: f22vmr_def = 1.500e-10
173  real (kind=kind_phys), parameter :: cl4vmr_def = 1.397e-10
175  real (kind=kind_phys), parameter :: f113vmr_def= 8.2000e-11
176 
177 ! --- ozone seasonal climatology parameters defined in module ozne_def
178 ! - 4x5 ozone data parameter
179 ! integer, parameter :: JMR=45, LOZ=17
180 ! real (kind=kind_phys), parameter :: blte=-86.0, dlte=4.0
181 ! - geos ozone data
182 ! integer, parameter :: JMR=18, LOZ=17
183 ! real (kind=kind_phys), parameter :: blte=-85.0, dlte=10.0
184 
185 ! --- module variables to be set in subroutin gas_init and/or gas_update
186 
188 
189  real (kind=kind_phys), allocatable :: pkstr(:), o3r(:,:,:)
190  integer :: k1oz = 0, k2oz = 0
191  real (kind=kind_phys) :: facoz = 0.0
192 
194 
195  real (kind=kind_phys), allocatable :: co2vmr_sav(:,:,:)
196  real (kind=kind_phys), allocatable :: co2cyc_sav(:,:,:)
197 
198  real (kind=kind_phys) :: co2_glb = co2vmr_def
199  real (kind=kind_phys) :: gco2cyc(12)
200  data gco2cyc(:) / 12*0.0 /
201 
202  integer :: kyrsav = 0
203  integer :: kmonsav = 1
204 
205 ! --- public interfaces
206 
207  public gas_init, gas_update, getgases, getozn
208 
209 
210 ! =================
211  contains
212 ! =================
213 
217 !-----------------------------------
218  subroutine gas_init &
219  & ( me )! --- inputs:
220 ! --- outputs: ( none )
221 
222 ! =================================================================== !
223 ! !
224 ! gas_init sets up ozone, co2, etc. parameters. if climatology ozone !
225 ! then read in monthly ozone data. !
226 ! !
227 ! inputs: dimemsion !
228 ! me - print message control flag 1 !
229 ! !
230 ! outputs: (to the module variables) !
231 ! ( none ) !
232 ! !
233 ! external module variables: (in physparam) !
234 ! ico2flg - co2 data source control flag !
235 ! =0: use prescribed co2 global mean value !
236 ! =1: use input global mean co2 value (co2_glb) !
237 ! =2: use input 2-d monthly co2 value (co2vmr_sav) !
238 ! ictmflg - =yyyy#, data ic time/date control flag !
239 ! = -2: same as 0, but superimpose seasonal cycle !
240 ! from climatology data set. !
241 ! = -1: use user provided external data for the fcst !
242 ! time, no extrapolation. !
243 ! = 0: use data at initial cond time, if not existed!
244 ! then use latest, without extrapolation. !
245 ! = 1: use data at the forecast time, if not existed!
246 ! then use latest and extrapolate to fcst time.!
247 ! =yyyy0: use yyyy data for the forecast time, no !
248 ! further data extrapolation. !
249 ! =yyyy1: use yyyy data for the fcst. if needed, do !
250 ! extrapolation to match the fcst time. !
251 ! ioznflg - ozone data control flag !
252 ! =0: use climatological ozone profile !
253 ! >0: use interactive ozone profile !
254 ! ivflip - vertical profile indexing flag !
255 ! co2usr_file- external co2 user defined data table !
256 ! co2cyc_file- external co2 climotology monthly cycle data table !
257 ! !
258 ! internal module variables: !
259 ! pkstr, o3r - arrays for climatology ozone data !
260 ! !
261 ! usage: call gas_init !
262 ! !
263 ! subprograms called: none !
264 ! !
265 ! =================================================================== !
266 !
267  implicit none
268 
269 ! --- inputs:
270  integer, intent(in) :: me
271 
272 ! --- output: ( none )
273 
274 ! --- locals:
275  real (kind=kind_phys), dimension(IMXCO2,JMXCO2) :: co2dat
276  real (kind=kind_phys) :: co2g1, co2g2
277  real (kind=kind_phys) :: pstr(loz)
278  real (kind=kind_io4) :: o3clim4(jmr,loz,12), pstr4(loz)
279 
280  integer :: imond(12), ilat(jmr,12)
281  integer :: i, j, k, iyr, imo
282  logical :: file_exist, lextpl
283  character :: cline*100, cform*8
284  data cform / '(24f7.2)' / !! data format in IMXCO2*f7.2
285 !
286 !===> ... begin here
287 !
288  if ( me == 0 ) print *, vtaggas ! print out version tag
289 
290  kyrsav = 0
291  kmonsav = 1
292 
293 ! --- ... climatology ozone data section
294 
295  if ( ioznflg > 0 ) then
296  if ( me == 0 ) then
297  print *,' - Using interactive ozone distribution'
298  endif
299  else
300  if ( timeozc /= 12 ) then
301  print *,' - Using climatology ozone distribution'
302  print *,' timeozc=',timeozc, ' is not monthly mean', &
303  & ' - job aborting in subroutin gas_init!!!'
304  stop
305  endif
306 
307  allocate (pkstr(loz), o3r(jmr,loz,12))
308  rewind nio3clm
309 
310  if ( loz == 17 ) then ! For the operational ozone climatology
311  do k = 1, loz
312  read (nio3clm,15) pstr4(k)
313  15 format(f10.3)
314  enddo
315 
316  do imo = 1, 12
317  do j = 1, jmr
318  read (nio3clm,16) imond(imo), ilat(j,imo), &
319  & (o3clim4(j,k,imo),k=1,10)
320  16 format(i2,i4,10f6.2)
321  read (nio3clm,20) (o3clim4(j,k,imo),k=11,loz)
322  20 format(6x,10f6.2)
323  enddo
324  enddo
325  else ! For newer ozone climatology
326  read (nio3clm)
327  do k = 1, loz
328  read (nio3clm) pstr4(k)
329  enddo
330 
331  do imo = 1, 12
332  do k = 1, loz
333  read (nio3clm) (o3clim4(j,k,imo),j=1,jmr)
334  enddo
335  enddo
336  endif ! end if_LOZ_block
337 !
338  do imo = 1, 12
339  do k = 1, loz
340  do j = 1, jmr
341  o3r(j,k,imo) = o3clim4(j,k,imo) * 1.655e-6
342  enddo
343  enddo
344  enddo
345 
346  do k = 1, loz
347  pstr(k) = pstr4(k)
348  enddo
349 
350  if ( me == 0 ) then
351  print *,' - Using climatology ozone distribution'
352  print *,' Found ozone data for levels pstr=', &
353  & (pstr(k),k=1,loz)
354 ! print *,' O3=',(o3r(15,k,1),k=1,LOZ)
355  endif
356 
357  do k = 1, loz
358  pkstr(k) = fpkapx(pstr(k)*100.0)
359  enddo
360  endif ! end if_ioznflg_block
361 
362 ! --- ... co2 data section
363 
364  co2_glb = co2vmr_def
365 
366  lab_ico2 : if ( ico2flg == 0 ) then
367 
368  if ( me == 0 ) then
369  print *,' - Using prescribed co2 global mean value=', &
370  & co2vmr_def
371  endif
372 
373  else lab_ico2
374 
375  lab_ictm : if ( ictmflg == -1 ) then ! input user provided data
376 
377  inquire (file=co2usr_file, exist=file_exist)
378  if ( .not. file_exist ) then
379  print *,' Can not find user CO2 data file: ',co2usr_file, &
380  & ' - Stopped in subroutine gas_init !!'
381  stop
382  else
383  close (nico2cn)
384  open(nico2cn,file=co2usr_file,form='formatted',status='old')
385  rewind nico2cn
386  read (nico2cn, 25) iyr, cline, co2g1, co2g2
387  25 format(i4,a94,f7.2,16x,f5.2)
388  co2_glb = co2g1 * 1.0e-6
389 
390  if ( ico2flg == 1 ) then
391  if ( me == 0 ) then
392  print *,' - Using co2 global annual mean value from', &
393  & ' user provided data set:',co2usr_file
394  print *, iyr,cline(1:94),co2g1,' GROWTH RATE =', co2g2
395  endif
396  elseif ( ico2flg == 2 ) then
397  allocate ( co2vmr_sav(imxco2,jmxco2,12) )
398 
399  do imo = 1, 12
400  read (nico2cn,cform) co2dat
401 !check print cform, co2dat
402 
403  do j = 1, jmxco2
404  do i = 1, imxco2
405  co2vmr_sav(i,j,imo) = co2dat(i,j) * 1.0e-6
406  enddo
407  enddo
408  enddo
409 
410  if ( me == 0 ) then
411  print *,' - Using co2 monthly 2-d data from user', &
412  & ' provided data set:',co2usr_file
413  print *, iyr,cline(1:94),co2g1,' GROWTH RATE =', co2g2
414 
415  print *,' CHECK: Sample of selected months of CO2 data'
416  do imo = 1, 12, 3
417  print *,' Month =',imo
418  print *, co2vmr_sav(1,:,imo)
419  enddo
420  endif
421  else
422  print *,' ICO2=',ico2flg,' is not a valid selection', &
423  & ' - Stoped in subroutine gas_init!!!'
424  stop
425  endif ! endif_ico2flg_block
426 
427  close (nico2cn)
428  endif ! endif_file_exist_block
429 
430  else lab_ictm ! input from observed data
431 
432  if ( ico2flg == 1 ) then
433  if ( me == 0 ) then
434  print *,' - Using observed co2 global annual mean value'
435  endiF
436  elseif ( ico2flg == 2 ) then
437  allocate ( co2vmr_sav(imxco2,jmxco2,12) )
438 
439  if ( me == 0 ) then
440  print *,' - Using observed co2 monthly 2-d data'
441  endif
442  else
443  print *,' ICO2=',ico2flg,' is not a valid selection', &
444  & ' - Stoped in subroutine gas_init!!!'
445  stop
446  endif
447 
448  if ( ictmflg == -2 ) then
449  inquire (file=co2cyc_file, exist=file_exist)
450  if ( .not. file_exist ) then
451  if ( me == 0 ) then
452  print *,' Can not find seasonal cycle CO2 data: ', &
453  & co2cyc_file,' - Stopped in subroutine gas_init !!'
454  endif
455  stop
456  else
457  allocate( co2cyc_sav(imxco2,jmxco2,12) )
458 
459 ! --- ... read in co2 2-d seasonal cycle data
460  close (nico2cn)
461  open (nico2cn,file=co2cyc_file,form='formatted', &
462  & status='old')
463  rewind nico2cn
464  read (nico2cn, 35) cline, co2g1, co2g2
465  35 format(a98,f7.2,16x,f5.2)
466  read (nico2cn,cform) co2dat ! skip annual mean part
467 
468  if ( me == 0 ) then
469  print *,' - Superimpose seasonal cycle to mean CO2 data'
470  print *,' Opened CO2 climatology seasonal cycle data',&
471  & ' file: ',co2cyc_file
472 !check print *, cline(1:98), co2g1, co2g2
473  endif
474 
475  do imo = 1, 12
476  read (nico2cn,45) cline, gco2cyc(imo)
477  45 format(a58,f7.2)
478 !check print *, cline(1:58),gco2cyc(imo)
479  gco2cyc(imo) = gco2cyc(imo) * 1.0e-6
480 
481  read (nico2cn,cform) co2dat
482 !check print cform, co2dat
483  do j = 1, jmxco2
484  do i = 1, imxco2
485  co2cyc_sav(i,j,imo) = co2dat(i,j) * 1.0e-6
486  enddo
487  enddo
488  enddo
489 
490  close (nico2cn)
491  endif ! endif_file_exist_block
492  endif
493 
494  endif lab_ictm
495  endif lab_ico2
496 
497  return
498 !
499 !...................................
500  end subroutine gas_init
501 !-----------------------------------
502 
514 !-----------------------------------
515  subroutine gas_update &
516  & ( iyear, imon, iday, ihour, loz1st, ldoco2, me )! --- inputs
517 ! --- outputs: ( none )
518 
519 ! =================================================================== !
520 ! !
521 ! gas_update reads in 2-d monthly co2 data set for a specified year. !
522 ! data are in a 15 degree lat/lon horizontal resolution. !
523 ! !
524 ! inputs: dimemsion !
525 ! iyear - year of the requested data for fcst 1 !
526 ! imon - month of the year 1 !
527 ! iday - day of the month 1 !
528 ! ihour - hour of the day 1 !
529 ! loz1st - clim ozone 1st time update control flag 1 !
530 ! ldoco2 - co2 update control flag 1 !
531 ! me - print message control flag 1 !
532 ! !
533 ! outputs: (to the module variables) !
534 ! ( none ) !
535 ! !
536 ! external module variables: (in physparam) !
537 ! ico2flg - co2 data source control flag !
538 ! =0: use prescribed co2 global mean value !
539 ! =1: use input global mean co2 value (co2_glb) !
540 ! =2: use input 2-d monthly co2 value (co2vmr_sav) !
541 ! ictmflg - =yyyy#, data ic time/date control flag !
542 ! = -2: same as 0, but superimpose seasonal cycle !
543 ! from climatology data set. !
544 ! = -1: use user provided external data for the fcst !
545 ! time, no extrapolation. !
546 ! = 0: use data at initial cond time, if not existed!
547 ! then use latest, without extrapolation. !
548 ! = 1: use data at the forecast time, if not existed!
549 ! then use latest and extrapolate to fcst time.!
550 ! =yyyy0: use yyyy data for the forecast time, no !
551 ! further data extrapolation. !
552 ! =yyyy1: use yyyy data for the fcst. if needed, do !
553 ! extrapolation to match the fcst time. !
554 ! ioznflg - ozone data control flag !
555 ! =0: use climatological ozone profile !
556 ! >0: use interactive ozone profile !
557 ! ivflip - vertical profile indexing flag !
558 ! co2dat_file- external co2 2d monthly obsv data table !
559 ! co2gbl_file- external co2 global annual mean data table !
560 ! !
561 ! internal module variables: !
562 ! co2vmr_sav - monthly co2 volume mixing ratio IMXCO2*JMXCO2*12 !
563 ! co2cyc_sav - monthly cycle co2 vol mixing ratio IMXCO2*JMXCO2*12 !
564 ! co2_glb - global annual mean co2 mixing ratio !
565 ! gco2cyc - global monthly mean co2 variation 12 !
566 ! k1oz,k2oz,facoz !
567 ! - climatology ozone parameters 1 !
568 ! !
569 ! usage: call gas_update !
570 ! !
571 ! subprograms called: none !
572 ! !
573 ! =================================================================== !
574 !
575  implicit none
576 
577 ! --- inputs:
578  integer, intent(in) :: iyear, imon, iday, ihour, me
579 
580  logical, intent(in) :: loz1st, ldoco2
581 
582 ! --- output: ( none )
583 
584 ! --- locals:
585  real (kind=kind_phys), dimension(IMXCO2,JMXCO2) :: co2dat, co2ann
586  real (kind=kind_phys) :: co2g1, co2g2, rate
587 
588  integer :: i, id, j, l, iyr, imo, iyr1, iyr2, jyr, idyr
589  integer, save :: mdays(13), midmon=15, midm=15, midp=45
590 ! --- number of days in a month
591  data mdays / 31,28,31,30,31,30,31,31,30,31,30,31,30 /
592 
593  logical :: file_exist, lextpl, change
594  character :: cline*100, cform*8, cfile1*26
595  data cform / '(24f7.2)' / !! data format in IMXCO2*f7.2
596 !
597 !===> ... begin here
598 !
600 
601  if ( ioznflg == 0 ) then
602  midmon = mdays(imon)/2 + 1
603  change = loz1st .or. ( (iday==midmon) .and. (ihour==0) )
604 !
605  if ( change ) then
606  if ( iday < midmon ) then
607  k1oz = mod(imon+10, 12) + 1
608  midm = mdays(k1oz)/2 + 1
609  k2oz = imon
610  midp = mdays(k1oz) + midmon
611  else
612  k1oz = imon
613  midm = midmon
614  k2oz = mod(imon, 12) + 1
615  midp = mdays(k2oz)/2 + 1 + mdays(k1oz)
616  endif
617  endif
618 !
619  if (iday < midmon) then
620  id = iday + mdays(k1oz)
621  else
622  id = iday
623  endif
624 
625  facoz = float(id - midm) / float(midp - midm)
626  endif
627 
629 
630  if ( ico2flg == 0 ) return ! use prescribed global mean co2 data
631  if ( ictmflg ==-1 ) return ! use user provided co2 data
632  if ( .not. ldoco2 ) return ! no need to update co2 data
633 
634  if ( ictmflg < 0 ) then ! use user provided external data
635  lextpl = .false. ! no time extrapolation
636  idyr = iyear ! use the model year
637  else ! use historically observed data
638  lextpl = ( mod(ictmflg,10) == 1 ) ! flag for data extrapolation
639  idyr = ictmflg / 10 ! year of data source used
640  if ( idyr == 0 ) idyr = iyear ! not specified, use model year
641  endif
642 
643 ! --- ... auto select co2 2-d data table for required year
644 
645  kmonsav = imon
646  if ( kyrsav == iyear ) return
647  kyrsav = iyear
648  iyr = iyear
649 
650 ! --- ... for data earlier than MINYEAR (1957), the data are in
651 ! the form of semi-yearly global mean values. otherwise,
652 ! data are monthly mean in horizontal 2-d map.
653 
654  lab_if_idyr : if ( idyr < minyear .and. ictmflg > 0 ) then
655 
656  if ( me == 0 ) then
657  print *,' Requested CO2 data year',iyear,' earlier than', &
658  & minyear
659  print *,' Which is the earliest monthly observation', &
660  & ' data available.'
661  print *,' Thus, historical global mean data is used'
662  endif
663 
664 ! --- ... check to see if requested co2 data file existed
665 
666  inquire (file=co2gbl_file, exist=file_exist)
667  if ( .not. file_exist ) then
668  print *,' Requested co2 data file "',co2gbl_file, &
669  & '" not found - Stopped in subroutine gas_update!!'
670  stop
671  else
672  close(nico2cn)
673  open (nico2cn,file=co2gbl_file,form='formatted',status='old')
674  rewind nico2cn
675 
676  read (nico2cn, 24) iyr1, iyr2, cline
677  24 format(i4,4x,i4,a48)
678 
679  if ( me == 0 ) then
680  print *,' Opened co2 data file: ',co2gbl_file
681 !check print *, iyr1, iyr2, cline(1:48)
682  endif
683 
684  if ( idyr < iyr1 ) then
685  iyr = iyr1
686 !check if ( me == 0 ) then
687 ! print *,' Using earlist available co2 data, year=',iyr1
688 !check endif
689  endif
690 
691  i = iyr2
692  lab_dowhile1 : do while ( i >= iyr1 )
693 ! read (NICO2CN,26) jyr, co2g1, co2g2
694 ! 26 format(i4,4x,2f7.2)
695  read (nico2cn, *) jyr, co2g1, co2g2
696 
697  if ( i == iyr .and. iyr == jyr ) then
698  co2_glb = (co2g1+co2g2) * 0.5e-6
699  if ( ico2flg == 2 ) then
700  do j = 1, jmxco2
701  do i = 1, imxco2
702  co2vmr_sav(i,j,1:6) = co2g1 * 1.0e-6
703  co2vmr_sav(i,j,7:12) = co2g2 * 1.0e-6
704  enddo
705  enddo
706  endif
707 
708  if ( me == 0 ) print *,' Co2 data for year',iyear, &
709  & co2_glb
710  exit lab_dowhile1
711  else
712 !check if ( me == 0 ) print *,' Skip co2 data for year',i
713  i = i - 1
714  endif
715  enddo lab_dowhile1
716 
717  close ( nico2cn )
718  endif ! end if_file_exist_block
719 
720  else lab_if_idyr
721 
722 ! --- ... set up input data file name
723 
724  cfile1 = co2dat_file
725  write(cfile1(19:22),34) idyr
726  34 format(i4.4)
727 
728 ! --- ... check to see if requested co2 data file existed
729 
730  inquire (file=cfile1, exist=file_exist)
731  if ( .not. file_exist ) then
732 
733  lab_if_ictm : if ( ictmflg > 10 ) then ! specified year of data not found
734  if ( me == 0 ) then
735  print *,' Specified co2 data for year',idyr, &
736  & ' not found !! Need to change namelist ICTM !!'
737  print *,' *** Stopped in subroutine gas_update !!'
738  endif
739  stop
740  else lab_if_ictm ! looking for latest available data
741  if ( me == 0 ) then
742  print *,' Requested co2 data for year',idyr, &
743  & ' not found, check for other available data set'
744  endif
745 
746  lab_dowhile2 : do while ( iyr >= minyear )
747  iyr = iyr - 1
748  write(cfile1(19:22),34) iyr
749 
750  inquire (file=cfile1, exist=file_exist)
751  if ( me == 0 ) then
752  print *,' Looking for CO2 file ',cfile1
753  endif
754 
755  if ( file_exist ) then
756  exit lab_dowhile2
757  endif
758  enddo lab_dowhile2
759 
760  if ( .not. file_exist ) then
761  if ( me == 0 ) then
762  print *,' Can not find co2 data source file'
763  print *,' *** Stopped in subroutine gas_update !!'
764  endif
765  stop
766  endif
767  endif lab_if_ictm
768  endif ! end if_file_exist_block
769 
770 ! --- ... read in co2 2-d data for the requested month
771 
772  close(nico2cn)
773  open (nico2cn,file=cfile1,form='formatted',status='old')
774  rewind nico2cn
775  read (nico2cn, 36) iyr, cline, co2g1, co2g2
776  36 format(i4,a94,f7.2,16x,f5.2)
777 
778  if ( me == 0 ) then
779  print *,' Opened co2 data file: ',cfile1
780  print *, iyr, cline(1:94), co2g1,' GROWTH RATE =', co2g2
781  endif
782 
783 ! --- ... add growth rate if needed
784  if ( lextpl ) then
785 ! rate = co2g2 * (iyear - iyr) ! rate from early year
786 ! rate = 1.60 * (iyear - iyr) ! avg rate over long period
787  rate = 2.00 * (iyear - iyr) ! avg rate for recent period
788  else
789  rate = 0.0
790  endif
791 
792  co2_glb = (co2g1 + rate) * 1.0e-6
793  if ( me == 0 ) then
794  print *,' Global annual mean CO2 data for year', &
795  & iyear, co2_glb
796  endif
797 
798  if ( ictmflg == -2 ) then ! need to calc ic time annual mean first
799 
800  if ( ico2flg == 1 ) then
801  if ( me==0 ) then
802  print *,' CHECK: Monthly deviations of climatology ', &
803  & 'to be superimposed on global annual mean'
804  print *, gco2cyc
805  endif
806  elseif ( ico2flg == 2 ) then
807  co2ann(:,:) = 0.0
808 
809  do imo = 1, 12
810  read (nico2cn,cform) co2dat
811 !check print cform, co2dat
812 
813  do j = 1, jmxco2
814  do i = 1, imxco2
815  co2ann(i,j) = co2ann(i,j) + co2dat(i,j)
816  enddo
817  enddo
818  enddo
819 
820  do j = 1, jmxco2
821  do i = 1, imxco2
822  co2ann(i,j) = co2ann(i,j) * 1.0e-6 / float(12)
823  enddo
824  enddo
825 
826  do imo = 1, 12
827  do j = 1, jmxco2
828  do i = 1, imxco2
829  co2vmr_sav(i,j,imo) = co2ann(i,j)+co2cyc_sav(i,j,imo)
830  enddo
831  enddo
832  enddo
833 
834  if ( me==0 ) then
835  print *,' CHECK: Sample of 2-d annual mean of CO2 ', &
836  & 'data used for year:',iyear
837  print *, co2ann(1,:)
838  print *,' CHECK: AFTER adding seasonal cycle, Sample ', &
839  & 'of selected months of CO2 data for year:',iyear
840  do imo = 1, 12, 3
841  print *,' Month =',imo
842  print *, co2vmr_sav(1,:,imo)
843  enddo
844  endif
845  endif ! endif_icl2flg_block
846 
847  else ! no need to calc ic time annual mean first
848 
849  if ( ico2flg == 2 ) then ! directly save monthly data
850  do imo = 1, 12
851  read (nico2cn,cform) co2dat
852 !check print cform, co2dat
853 
854  do j = 1, jmxco2
855  do i = 1, imxco2
856  co2vmr_sav(i,j,imo) = (co2dat(i,j) + rate) * 1.0e-6
857  enddo
858  enddo
859  enddo
860 
861  if ( me == 0 ) then
862  print *,' CHECK: Sample of selected months of CO2 ', &
863  & 'data used for year:',iyear
864  do imo = 1, 12, 3
865  print *,' Month =',imo
866  print *, co2vmr_sav(1,:,imo)
867  enddo
868  endif
869  endif ! endif_ico2flg_block
870 
871  do imo = 1, 12
872  gco2cyc(imo) = 0.0
873  enddo
874  endif ! endif_ictmflg_block
875  close ( nico2cn )
876 
877  endif lab_if_idyr
878 
879  return
880 !
881 !...................................
882  end subroutine gas_update
883 !-----------------------------------
884 !! @}
885 
907 !-----------------------------------
908  subroutine getgases &
909  & ( plvl, xlon, xlat, & ! --- inputs
910  & imax, lmax, &
911  & gasdat & ! --- outputs
912  & )
913 
914 ! =================================================================== !
915 ! !
916 ! getgases set up global distribution of radiation absorbing gases !
917 ! in volume mixing ratio. currently only co2 has the options from !
918 ! observed values, all other gases are asigned to the climatological !
919 ! values. !
920 ! !
921 ! inputs: !
922 ! plvl(IMAX,LMAX+1)- pressure at model layer interfaces (mb) !
923 ! xlon(IMAX) - grid longitude in radians, ok both 0->2pi or !
924 ! -pi -> +pi arrangements !
925 ! xlat(IMAX) - grid latitude in radians, default range to !
926 ! pi/2 -> -pi/2, otherwise see in-line comment !
927 ! IMAX, LMAX - horiz, vert dimensions for output data !
928 ! !
929 ! outputs: !
930 ! gasdat(IMAX,LMAX,NF_VGAS) - gases volume mixing ratioes !
931 ! (:,:,1) - co2 !
932 ! (:,:,2) - n2o !
933 ! (:,:,3) - ch4 !
934 ! (:,:,4) - o2 !
935 ! (:,:,5) - co !
936 ! (:,:,6) - cfc11 !
937 ! (:,:,7) - cfc12 !
938 ! (:,:,8) - cfc22 !
939 ! (:,:,9) - ccl4 !
940 ! (:,:,10) - cfc113 !
941 ! !
942 ! external module variables: (in physparam) !
943 ! ico2flg - co2 data source control flag !
944 ! =0: use prescribed co2 global mean value !
945 ! =1: use input global mean co2 value (co2_glb) !
946 ! =2: use input 2-d monthly co2 value (co2vmr_sav) !
947 ! ivflip - vertical profile indexing flag !
948 ! !
949 ! internal module variables used: !
950 ! co2vmr_sav - saved monthly co2 concentration from sub gas_update !
951 ! co2_glb - saved global annual mean co2 value from gas_update !
952 ! gco2cyc - saved global seasonal variation of co2 climatology !
953 ! in 12-month form !
954 ! ** note: for lower atmos co2vmr_sav may have clim monthly deviations !
955 ! superimposed on init-cond co2 value, while co2_glb only !
956 ! contains the global mean value, thus needs to add the !
957 ! monthly dglobal mean deviation gco2cyc at upper atmos. for !
958 ! ictmflg/=-2, this value will be zero. !
959 ! !
960 ! usage: call getgases !
961 ! !
962 ! subprograms called: none !
963 ! !
964 ! =================================================================== !
965 !
966  implicit none
967 
968 ! --- input:
969  integer, intent(in) :: IMAX, LMAX
970  real (kind=kind_phys), intent(in) :: plvl(:,:), xlon(:), xlat(:)
971 
972 ! --- output:
973  real (kind=kind_phys), intent(out) :: gasdat(:,:,:)
974 
975 ! --- local:
976  integer :: i, k, ilat, ilon
977 
978  real (kind=kind_phys) :: xlon1, xlat1, tmp
979 
980 !===> ... begin here
981 
982 ! --- ... assign default values
983 
984  do k = 1, lmax
985  do i = 1, imax
986  gasdat(i,k,1) = co2vmr_def
987  gasdat(i,k,2) = n2ovmr_def
988  gasdat(i,k,3) = ch4vmr_def
989  gasdat(i,k,4) = o2vmr_def
990  gasdat(i,k,5) = covmr_def
991  gasdat(i,k,6) = f11vmr_def
992  gasdat(i,k,7) = f12vmr_def
993  gasdat(i,k,8) = f22vmr_def
994  gasdat(i,k,9) = cl4vmr_def
995  gasdat(i,k,10)= f113vmr_def
996  enddo
997  enddo
998 
999 ! --- ... co2 section
1000 
1001  if ( ico2flg == 1 ) then
1002 ! --- use obs co2 global annual mean value only
1003 
1004  do k = 1, lmax
1005  do i = 1, imax
1006  gasdat(i,k,1) = co2_glb + gco2cyc(kmonsav)
1007  enddo
1008  enddo
1009 
1010  elseif ( ico2flg == 2 ) then
1011 ! --- use obs co2 monthly data with 2-d variation at lower atmos
1012 ! otherwise use global mean value
1013 
1014  tmp = raddeg / resco2
1015  do i = 1, imax
1016  xlon1 = xlon(i)
1017  if ( xlon1 < 0.0 ) xlon1 = xlon1 + con_pi ! if xlon in -pi->pi, convert to 0->2pi
1018  xlat1 = hfpi - xlat(i) ! if xlat in pi/2 -> -pi/2 range
1019 !note xlat1 = xlat(i) ! if xlat in 0 -> pi range
1020 
1021  ilon = min( imxco2, int( xlon1*tmp + 1 ))
1022  ilat = min( jmxco2, int( xlat1*tmp + 1 ))
1023 
1024  if ( ivflip == 0 ) then ! index from toa to sfc
1025  do k = 1, lmax
1026  if ( plvl(i,k) >= prsco2 ) then
1027  gasdat(i,k,1) = co2vmr_sav(ilon,ilat,kmonsav)
1028  else
1029  gasdat(i,k,1) = co2_glb + gco2cyc(kmonsav)
1030  endif
1031  enddo
1032  else ! index from sfc to toa
1033  do k = 1, lmax
1034  if ( plvl(i,k+1) >= prsco2 ) then
1035  gasdat(i,k,1) = co2vmr_sav(ilon,ilat,kmonsav)
1036  else
1037  gasdat(i,k,1) = co2_glb + gco2cyc(kmonsav)
1038  endif
1039  enddo
1040  endif
1041  enddo
1042  endif
1043 
1044 !
1045  return
1046 !...................................
1047  end subroutine getgases
1048 !-----------------------------------
1049 
1058 !-----------------------------------
1059  subroutine getozn &
1060  & ( prslk,xlat, & ! --- inputs
1061  & imax, lm, &
1062  & o3mmr & ! --- outputs
1063  & )
1064 
1065 ! =================================================================== !
1066 ! !
1067 ! getozn sets up climatological ozone profile for radiation calculation!
1068 ! !
1069 ! this code is originally written By Shrinivas Moorthi !
1070 ! !
1071 ! inputs: !
1072 ! prslk (IMAX,LM) - exner function = (p/p0)**rocp !
1073 ! xlat (IMAX) - latitude in radians, default to pi/2 -> -pi/2 !
1074 ! range, otherwise see in-line comment !
1075 ! IMAX, LM - horizontal and vertical dimensions !
1076 ! !
1077 ! outputs: !
1078 ! o3mmr (IMAX,LM) - output ozone profile in mass mixing ratio (g/g)!
1079 ! !
1080 ! module variables: !
1081 ! k1oz, k2oz - ozone data interpolation indices !
1082 ! facoz - ozone data interpolation factor !
1083 ! ivflip - control flag for direction of vertical index !
1084 ! !
1085 ! usage: call getozn !
1086 ! !
1087 ! =================================================================== !
1088 !
1089  implicit none
1090 
1091 ! --- inputs:
1092  integer, intent(in) :: IMAX, LM
1093 
1094  real (kind=kind_phys), intent(in) :: prslk(:,:), xlat(:)
1095 
1096 ! --- outputs:
1097  real (kind=kind_phys), intent(out) :: o3mmr(:,:)
1098 
1099 ! --- locals:
1100  real (kind=kind_phys) :: o3i(imax,loz), wk1(imax), deglat, elte, &
1101  & tem, tem1, tem2, tem3, tem4, temp
1102  integer :: i, j, k, l, j1, j2, ll
1103 !
1104 !===> ... begin here
1105 !
1106  elte = blte + (jmr-1)*dlte
1107 
1108  do i = 1, imax
1109  deglat = xlat(i) * raddeg ! if xlat in pi/2 -> -pi/2 range
1110 ! deglat = 90.0 - xlat(i)*raddeg ! if xlat in 0 -> pi range
1111 
1112  if (deglat > blte .and. deglat < elte) then
1113  tem1 = (deglat - blte) / dlte + 1
1114  j1 = tem1
1115  j2 = j1 + 1
1116  tem1 = tem1 - j1
1117  elseif (deglat <= blte) then
1118  j1 = 1
1119  j2 = 1
1120  tem1 = 1.0
1121  elseif (deglat >= elte) then
1122  j1 = jmr
1123  j2 = jmr
1124  tem1 = 1.0
1125  endif
1126 
1127  tem2 = 1.0 - tem1
1128  do j = 1, loz
1129  tem3 = tem2*o3r(j1,j,k1oz) + tem1*o3r(j2,j,k1oz)
1130  tem4 = tem2*o3r(j1,j,k2oz) + tem1*o3r(j2,j,k2oz)
1131  o3i(i,j) = tem4*facoz + tem3*(1.0 - facoz)
1132  enddo
1133  enddo
1134 
1135  do l = 1, lm
1136  ll = l
1137  if (ivflip == 1) ll = lm -l + 1
1138 
1139  do i = 1, imax
1140  wk1(i) = prslk(i,ll)
1141  enddo
1142 
1143  do k = 1, loz-1
1144  temp = 1.0 / (pkstr(k+1) - pkstr(k))
1145 
1146  do i = 1, imax
1147  if (wk1(i) > pkstr(k) .and. wk1(i) <= pkstr(k+1)) then
1148  tem = (pkstr(k+1) - wk1(i)) * temp
1149  o3mmr(i,ll) = tem * o3i(i,k) + (1.0 - tem) * o3i(i,k+1)
1150  endif
1151  enddo
1152  enddo
1153 
1154  do i = 1, imax
1155  if (wk1(i) > pkstr(loz)) o3mmr(i,ll) = o3i(i,loz)
1156  if (wk1(i) < pkstr(1)) o3mmr(i,ll) = o3i(i,1)
1157  enddo
1158  enddo
1159 !
1160  return
1161 !...................................
1162  end subroutine getozn
1163 !-----------------------------------
1164 
1165 !
1166 !........................................!
1167  end module module_radiation_gases !
1168 !========================================!
physparam::co2dat_file
character, save co2dat_file
external co2 2d monthly obsv data table: co2historicaldata_2004.txt
Definition: physparam.f:222
physcons::con_pi
real(kind=kind_phys), parameter con_pi
pi
Definition: physcons.f:48
module_radiation_gases::imxco2
integer, parameter imxco2
input co2 dat lon points
Definition: radiation_gases.f:144
module_radiation_gases::resco2
real(kind=kind_phys), parameter resco2
horizontal resolution in degree
Definition: radiation_gases.f:151
physparam::co2gbl_file
character, save co2gbl_file
external co2 global annual mean data tb: co2historicaldata_glob.txt
Definition: physparam.f:224
module_radiation_gases::gas_update
subroutine, public gas_update(iyear, imon, iday, ihour, loz1st, ldoco2, me)
This subroutine reads in 2-d monthly co2 data set for a specified year. Data are in a 15 degree lat/l...
Definition: radiation_gases.f:517
physparam::co2usr_file
character, save co2usr_file
external co2 user defined data table: co2userdata.txt
Definition: physparam.f:226
module_radiation_gases::f22vmr_def
real(kind=kind_phys), parameter f22vmr_def
aer 2003 value
Definition: radiation_gases.f:171
module_radiation_gases::cl4vmr_def
real(kind=kind_phys), parameter cl4vmr_def
aer 2003 value
Definition: radiation_gases.f:173
physparam::ico2flg
integer, save ico2flg
co2 data source control flag =0:prescribed value(380 ppmv) =1:yearly global averaged annual mean ...
Definition: physparam.f:205
module_radiation_gases::prsco2
real(kind=kind_phys), parameter prsco2
pressure limitation for 2-d co2 (mb)
Definition: radiation_gases.f:155
module_radiation_gases::f12vmr_def
real(kind=kind_phys), parameter f12vmr_def
aer 2003 value
Definition: radiation_gases.f:169
module_radiation_gases::f11vmr_def
real(kind=kind_phys), parameter f11vmr_def
aer 2003 value
Definition: radiation_gases.f:167
module_radiation_gases::gas_init
subroutine, public gas_init(me)
This subroutine sets up ozone, co2, etc. parameters. If climatology ozone then read in monthly ozone ...
Definition: radiation_gases.f:220
module_radiation_gases::getgases
subroutine, public getgases(plvl, xlon, xlat, IMAX, LMAX, gasdat )
This subroutine sets up global distribution of radiation absorbing gases in volume mixing ratio...
Definition: radiation_gases.f:913
module_radiation_gases::getozn
subroutine, public getozn(prslk, xlat, IMAX, LM, o3mmr )
This subroutine sets up climatological ozone profile for radiation calculation. This code is original...
Definition: radiation_gases.f:1064
physparam::ioznflg
integer, save ioznflg
ozone data source control flag =0:use seasonal climatology ozone data >0:use prognostic ozone sch...
Definition: physparam.f:220
module_radiation_gases::minyear
integer, parameter minyear
earlist year 2-d co2 data available
Definition: radiation_gases.f:148
physparam::co2cyc_file
character, save co2cyc_file
external co2 clim monthly cycle data tb: co2monthlycyc.txt
Definition: physparam.f:228
physparam::ictmflg
integer, save ictmflg
controls external data at initial time and data usage during forecast time =-2:as in 0...
Definition: physparam.f:215
module_radiation_gases::jmxco2
integer, parameter jmxco2
input co2 data lat points
Definition: radiation_gases.f:146
module_radiation_gases::raddeg
real(kind=kind_phys), parameter raddeg
rad->deg conversion
Definition: radiation_gases.f:153
module_radiation_gases::f113vmr_def
real(kind=kind_phys), parameter f113vmr_def
gfdl 1999 value
Definition: radiation_gases.f:175
module_radiation_gases::nf_vgas
integer, parameter, public nf_vgas
number of gas species
Definition: radiation_gases.f:142
physparam::ivflip
integer, save ivflip
vertical profile indexing flag
Definition: physparam.f:289
module_radiation_gases::hfpi
real(kind=kind_phys), parameter hfpi
half of pi
Definition: radiation_gases.f:157