GFS Operational Physics Documentation  Revision: 81451
radiation_aerosols.f
Go to the documentation of this file.
1 
4 
5 ! ========================================================== !!!!!
6 ! 'module_radiation_aerosols' description !!!!!
7 ! ========================================================== !!!!!
8 ! !
9 ! this module contains climatological atmospheric aerosol schemes for!
10 ! radiation computations. !
11 ! !
12 ! in the module, the externally callable subroutines are : !
13 ! !
14 ! 'aer_init' -- initialization !
15 ! inputs: !
16 ! ( NLAY, me ) !
17 ! outputs: !
18 ! ( none ) !
19 ! !
20 ! 'aer_update' -- updating aerosol data !
21 ! inputs: !
22 ! ( iyear, imon, me ) !
23 ! outputs: !
24 ! ( none ) !
25 ! !
26 ! 'setaer' -- mapping aeros profile, compute aeros opticals !
27 ! inputs: !
28 ! (prsi,prsl,prslk,tvly,rhlay,slmsk,tracer,xlon,xlat, !
29 ! IMAX,NLAY,NLP1, lsswr,lslwr, !
30 ! outputs: !
31 ! (aerosw,aerolw,tau_gocart) !
32 !! (aerosw,aerolw,aerodp) !
33 ! !
34 ! !
35 ! external modules referenced: !
36 ! 'module physparam' in 'physparam.f' !
37 ! 'module physcons' in 'physcons.f' !
38 ! 'module module_radsw_parameters' in 'radsw_xxxx#_param.f' !
39 ! 'module module_radlw_parameters' in 'radlw_xxxx#_param.f' !
40 ! 'module module_radlw_cntr_para' in 'radsw_xxxx#_param.f' !
41 ! !
42 ! output variable definitions: !
43 ! aerosw(IMAX,NLAY,NBDSW,1) - aerosols optical depth for sw !
44 ! aerosw(IMAX,NLAY,NBDSW,2) - aerosols single scat albedo for sw !
45 ! aerosw(IMAX,NLAY,NBDSW,3) - aerosols asymmetry parameter for sw!
46 ! !
47 ! aerolw(IMAX,NLAY,NBDLW,1) - aerosols optical depth for lw !
48 ! aerolw(IMAX,NLAY,NBDLW,2) - aerosols single scattering albedo !
49 ! aerolw(IMAX,NLAY,NBDLW,3) - aerosols asymetry parameter !
50 ! !
51 ! !
52 ! program history: !
53 ! apr 2003 --- y.-t. hou created !
54 ! nov 04, 2003 --- y.-t. hou modified version !
55 ! apr 15, 2005 --- y.-t. hou modified module structure !
56 ! jul 2006 --- y.-t. hou add volcanic forcing !
57 ! feb 2007 --- y.-t. hou add generalized spectral band !
58 ! interpolation for sw aerosol optical properties !
59 ! mar 2007 --- y.-t. hou add generalized spectral band !
60 ! interpolation for lw aerosol optical properties !
61 ! aug 2007 --- y.-t. hou change clim-aer vert domain !
62 ! from pressure reference to sigma reference !
63 ! sep 2007 --- y.-t. hou moving temporary allocatable !
64 ! module variable arrays to subroutine dynamically !
65 ! allocated arrays (eliminate deallocate calls) !
66 ! jan 2010 --- sarah lu add gocart option !
67 ! may 2010 --- sarah lu add geos4-gocart climo !
68 ! jul 2010 -- s. moorthi - merged NEMS version with new GFS !
69 ! version !
70 ! oct 23, 2010 --- Hsin-mu lin modified subr setclimaer to !
71 ! interpolate the 5 degree aerosol data to small domain based on!
72 ! the nearby 4 points instead of previous nearby assignment by !
73 ! using the 5 degree data. this process will eliminate the dsw !
74 ! jagged edges in the east conus where aerosol effect are lagre.!
75 ! dec 2010 --- y.-t. hou modified and optimized bi-linear!
76 ! horizontal interpolation in subr setclimaer. added safe guard !
77 ! measures in lat/lon indexing and added sea/land mask variable !
78 ! slmsk as input field to help aerosol profile selection. !
79 ! jan 2011 --- y.-t. hou divided the program into two !
80 ! separated interchangeable modules: a climatology aerosol !
81 ! module, and a gocart aerosol scheme module. the stratospheric !
82 ! volcanic aerosol part is still within the two driver modules, !
83 ! and may also become a separate one in the further development.!
84 ! unified in/out argument list for both clim and gocart types of!
85 ! schemes and added vertically integrated aer-opt-dep, aerodp, !
86 ! to replace tau_gocart as optional output for various species. !
87 ! aug 2012 --- y.-t. hou changed the initialization subr !
88 ! 'aerinit' into two parts: 'aer_init' is called at the start !
89 ! of run to set up module parameters; and 'aer_update' is !
90 ! called within the time loop to check and update data sets. !
91 ! nov 2012 --- y.-t. hou modified control parameters thru!
92 ! module 'physparam'. !
93 ! jan 2013 --- sarah lu and y.-t. hou reintegrate both !
94 ! opac-clim and gocart schemes into one module to make the !
95 ! program best utilize common components. added aerosol model !
96 ! scheme selection control variable iaer_mdl to the namelist. !
97 ! Aug 2013 --- s. moorthi - merge sarah's gocart changes with!
98 ! yutai's changes !
99 ! 13Feb2014 --- Sarah lu - compute aod at 550nm !
100 ! !
101 ! references for opac climatological aerosols: !
102 ! hou et al. 2002 (ncep office note 441) !
103 ! hess et al. 1998 - bams v79 831-844 !
104 ! !
105 ! references for gocart interactive aerosols: !
106 ! chin et al., 2000 - jgr, v105, 24671-24687 !
107 ! !
108 ! references for stratosperic volcanical aerosols: !
109 ! sato et al. 1993 - jgr, v98, d12, 22987-22994 !
110 ! !
111 ! !
112 !!!!! ========================================================== !!!!!
113 !!!!! end descriptions !!!!!
114 !!!!! ========================================================== !!!!!
115 
116 
117 
141 !========================================!
142  module module_radiation_aerosols !
143 !........................................!
144 !
145  use physparam,only : iaermdl, iaerflg, lavoflg, lalwflg, laswflg, &
146  & lalw1bd, aeros_file, ivflip, kind_phys &
147  &, kind_io4, kind_io8
148  use physcons, only : con_pi, con_rd, con_g, con_t0c, con_c, &
149  & con_boltz, con_plnk, con_amd
150 
151  use module_iounitdef, only : niaercm
152  use module_radsw_parameters, only : nbdsw, wvnsw1=>wvnum1, &
153  & nswstr, wvnsw2=>wvnum2
154  use module_radlw_parameters, only : nbdlw, wvnlw1, wvnlw2
155 !
156  use funcphys, only : fpkap
157  use gfs_phy_tracer_config, only : gfs_phy_tracer, trcindx
158 !
159  implicit none
160 !
161  private
162 
163 ! --- version tag and last revision date
164  character(40), parameter :: &
165  & VTAGAER='NCEP-Radiation_aerosols v5.2 Jan 2013 '
166 ! & VTAGAER='NCEP-Radiation_aerosols v5.1 Nov 2012 '
167 ! & VTAGAER='NCEP-Radiation_aerosols v5.0 Aug 2012 '
168 
169 ! --- general use parameter constants:
171  integer, parameter, public :: nf_aesw = 3
173  integer, parameter, public :: nf_aelw = 3
175  integer, parameter, public :: nlwstr = 1
177  integer, parameter, public :: nspc = 5
179  integer, parameter, public :: nspc1 = nspc + 1
180 
181  real (kind=kind_phys), parameter :: f_zero = 0.0
182  real (kind=kind_phys), parameter :: f_one = 1.0
183 
184 ! --- module control parameters set in subroutine "aer_init"
187  integer, save :: nswbnd = nbdsw
190  integer, save :: nlwbnd = nbdlw
192  integer, save :: nswlwbd = nbdsw+nbdlw
193 
194 ! --------------------------------------------------------------------- !
195 ! section-1 : module variables for spectral band interpolation !
196 ! similar to gfdl-sw treatment (2000 version) !
197 ! --------------------------------------------------------------------- !
198 
199 ! --- parameter constants:
201  integer, parameter, public :: nwvsol = 151
202 
204  integer, parameter, public :: nwvtot = 57600
206  integer, parameter, public :: nwvtir = 4000
207 
209  integer, dimension(NWVSOL), save :: nwvns0
210 
211  data nwvns0 / 100, 11, 14, 18, 24, 33, 50, 83, 12, 12, &
212  & 13, 15, 15, 17, 18, 20, 21, 24, 26, 30, 32, 37, 42, &
213  & 47, 55, 64, 76, 91, 111, 139, 179, 238, 333, 41, 42, 45, &
214  & 46, 48, 51, 53, 55, 58, 61, 64, 68, 71, 75, 79, 84, &
215  & 89, 95, 101, 107, 115, 123, 133, 142, 154, 167, 181, 197, 217, &
216  & 238, 263, 293, 326, 368, 417, 476, 549, 641, 758, 909, 101, 103, &
217  & 105, 108, 109, 112, 115, 117, 119, 122, 125, 128, 130, 134, 137, &
218  & 140, 143, 147, 151, 154, 158, 163, 166, 171, 175, 181, 185, 190, &
219  & 196, 201, 207, 213, 219, 227, 233, 240, 248, 256, 264, 274, 282, &
220  & 292, 303, 313, 325, 337, 349, 363, 377, 392, 408, 425, 444, 462, &
221  & 483, 505, 529, 554, 580, 610, 641, 675, 711, 751, 793, 841, 891, &
222  & 947,1008,1075,1150,1231,1323,1425,1538,1667,1633,14300 /
223 
225  real (kind=kind_phys), dimension(NWVSOL), save :: s0intv
226 
227  data s0intv( 1: 50) / &
228  & 1.60000e-6, 2.88000e-5, 3.60000e-5, 4.59200e-5, 6.13200e-5, &
229  & 8.55000e-5, 1.28600e-4, 2.16000e-4, 2.90580e-4, 3.10184e-4, &
230  & 3.34152e-4, 3.58722e-4, 3.88050e-4, 4.20000e-4, 4.57056e-4, &
231  & 4.96892e-4, 5.45160e-4, 6.00600e-4, 6.53600e-4, 7.25040e-4, &
232  & 7.98660e-4, 9.11200e-4, 1.03680e-3, 1.18440e-3, 1.36682e-3, &
233  & 1.57560e-3, 1.87440e-3, 2.25500e-3, 2.74500e-3, 3.39840e-3, &
234  & 4.34000e-3, 5.75400e-3, 7.74000e-3, 9.53050e-3, 9.90192e-3, &
235  & 1.02874e-2, 1.06803e-2, 1.11366e-2, 1.15830e-2, 1.21088e-2, &
236  & 1.26420e-2, 1.32250e-2, 1.38088e-2, 1.44612e-2, 1.51164e-2, &
237  & 1.58878e-2, 1.66500e-2, 1.75140e-2, 1.84450e-2, 1.94106e-2 /
238  data s0intv( 51:100) / &
239  & 2.04864e-2, 2.17248e-2, 2.30640e-2, 2.44470e-2, 2.59840e-2, &
240  & 2.75940e-2, 2.94138e-2, 3.13950e-2, 3.34800e-2, 3.57696e-2, &
241  & 3.84054e-2, 4.13490e-2, 4.46880e-2, 4.82220e-2, 5.22918e-2, &
242  & 5.70078e-2, 6.19888e-2, 6.54720e-2, 6.69060e-2, 6.81226e-2, &
243  & 6.97788e-2, 7.12668e-2, 7.27100e-2, 7.31610e-2, 7.33471e-2, &
244  & 7.34814e-2, 7.34717e-2, 7.35072e-2, 7.34939e-2, 7.35202e-2, &
245  & 7.33249e-2, 7.31713e-2, 7.35462e-2, 7.36920e-2, 7.23677e-2, &
246  & 7.25023e-2, 7.24258e-2, 7.20766e-2, 7.18284e-2, 7.32757e-2, &
247  & 7.31645e-2, 7.33277e-2, 7.36128e-2, 7.33752e-2, 7.28965e-2, &
248  & 7.24924e-2, 7.23307e-2, 7.21050e-2, 7.12620e-2, 7.10903e-2 /
249  data s0intv(101:151) / 7.12714e-2, &
250  & 7.08012e-2, 7.03752e-2, 7.00350e-2, 6.98639e-2, 6.90690e-2, &
251  & 6.87621e-2, 6.52080e-2, 6.65184e-2, 6.60038e-2, 6.47615e-2, &
252  & 6.44831e-2, 6.37206e-2, 6.24102e-2, 6.18698e-2, 6.06320e-2, &
253  & 5.83498e-2, 5.67028e-2, 5.51232e-2, 5.48645e-2, 5.12340e-2, &
254  & 4.85581e-2, 4.85010e-2, 4.79220e-2, 4.44058e-2, 4.48718e-2, &
255  & 4.29373e-2, 4.15242e-2, 3.81744e-2, 3.16342e-2, 2.99615e-2, &
256  & 2.92740e-2, 2.67484e-2, 1.76904e-2, 1.40049e-2, 1.46224e-2, &
257  & 1.39993e-2, 1.19574e-2, 1.06386e-2, 1.00980e-2, 8.63808e-3, &
258  & 6.52736e-3, 4.99410e-3, 4.39350e-3, 2.21676e-3, 1.33812e-3, &
259  & 1.12320e-3, 5.59000e-4, 3.60000e-4, 2.98080e-4, 7.46294e-5 /
260 
261 ! --------------------------------------------------------------------- !
262 ! section-2 : module variables for stratospheric volcanic aerosols !
263 ! from historical data (sato et al. 1993) !
264 ! --------------------------------------------------------------------- !
265 
266 ! --- parameter constants:
268  integer, parameter :: minvyr = 1850
270  integer, parameter :: maxvyr = 1999
271 
273  integer, allocatable, save :: ivolae(:,:,:)
274 
275 ! --- static control variables:
277  integer :: kyrstr
279  integer :: kyrend
281  integer :: kyrsav
283  integer :: kmonsav
284 
285 ! --------------------------------------------------------------------- !
286 ! section-3 : module variables for opac climatological aerosols !
287 ! optical properties (hess et al. 1989) !
288 ! --------------------------------------------------------------------- !
289 
290 ! --- parameters and constants:
292  integer, parameter :: nxc = 5
294  integer, parameter :: nae = 7
296  integer, parameter :: ndm = 5
298  integer, parameter :: imxae = 72
300  integer, parameter :: jmxae = 37
302  integer, parameter :: naerbnd=61
304  integer, parameter :: nrhlev =8
306  integer, parameter :: ncm1 = 6
308  integer, parameter :: ncm2 = 4
309  integer, parameter :: ncm = ncm1+ncm2
310 
312  real (kind=kind_phys), dimension(NRHLEV), save :: rhlev
313  data rhlev(:) / 0.0, 0.5, 0.7, 0.8, 0.9, 0.95, 0.98, 0.99 /
314 
315 ! --- the following arrays are for climatological data that are
316 ! allocated and read in subroutine 'clim_aerinit'.
317 ! - global aerosol distribution:
318 ! haer (NDM,NAE) - scale height of aerosols (km)
319 ! prsref(NDM,NAE) - ref pressure lev (sfc to toa) in mb (100Pa)
320 ! sigref(NDM,NAE) - ref sigma lev (sfc to toa)
321 
323  real (kind=kind_phys), save, dimension(NDM,NAE) :: haer
325  real (kind=kind_phys), save, dimension(NDM,NAE) :: prsref
327  real (kind=kind_phys), save, dimension(NDM,NAE) :: sigref
328 
329 ! --- the following arrays are allocate and setup in subr 'clim_aerinit'
330 ! - for relative humidity independent aerosol optical properties:
331 ! species : insoluble (inso); soot (soot);
332 ! mineral nuc mode (minm); mineral acc mode (miam);
333 ! mineral coa mode (micm); mineral transport(mitr).
334 ! extrhi(NCM1,NSWLWBD) - extinction coefficient for sw+lw spectral band
335 ! scarhi(NCM1,NSWLWBD) - scattering coefficient for sw+lw spectral band
336 ! ssarhi(NCM1,NSWLWBD) - single scattering albedo for sw+lw spectral band
337 ! asyrhi(NCM1,NSWLWBD) - asymmetry parameter for sw+lw spectral band
338 ! - for relative humidity dependent aerosol optical properties:
339 ! species : water soluble (waso); sea salt acc mode(ssam);
340 ! sea salt coa mode(sscm); sulfate droplets (suso).
341 ! rh level: 00%, 50%, 70%, 80%, 90%, 95%, 98%, 99%
342 ! extrhd(NRHLEV,NCM2,NSWLWBD) - extinction coefficient for sw+lw band
343 ! scarhd(NRHLEV,NCM2,NSWLWBD) - scattering coefficient for sw+lw band
344 ! ssarhd(NRHLEV,NCM2,NSWLWBD) - single scattering albedo for sw+lw band
345 ! asyrhd(NRHLEV,NCM2,NSWLWBD) - asymmetry parameter for sw+lw band
346 ! - for stratospheric aerosols optical properties:
347 ! extstra(NSWLWBD) - extinction coefficient for sw+lw band
348 
349  real (kind=kind_phys), allocatable, save, dimension(:,:) :: &
350  & extrhi, scarhi, ssarhi, asyrhi
351  real (kind=kind_phys), allocatable, save, dimension(:,:,:) :: &
352  & extrhd, scarhd, ssarhd, asyrhd
353  real (kind=kind_phys), allocatable, save, dimension(:) :: &
354  & extstra
355 
356 ! --- the following arrays are calculated in subr 'clim_aerinit'
357 ! - for topospheric aerosol profile distibution:
358 ! kprfg ( IMXAE*JMXAE) - aeros profile index
359 ! idxcg (NXC*IMXAE*JMXAE) - aeros component index
360 ! cmixg (NXC*IMXAE*JMXAE) - aeros component mixing ratio
361 ! denng ( 2 *IMXAE*JMXAE) - aerosols number density
362 
364 
366  real (kind=kind_phys), dimension(NXC,IMXAE,JMXAE), save :: cmixg
368  real (kind=kind_phys), dimension( 2 ,IMXAE,JMXAE), save :: denng
370  integer, dimension(NXC,IMXAE,JMXAE), save :: idxcg
372  integer, dimension( IMXAE,JMXAE), save :: kprfg
373 
374 ! --------------------------------------------------------------------- !
375 ! section-4 : module variables for gocart aerosol optical properties !
376 ! --------------------------------------------------------------------- !
377 
379 
380 ! --- parameters and constants:
381 ! - KCM, KCM1, KCM2 are determined from subroutine 'set_aerspc'
383  integer, parameter :: kaerbnd=61
385  integer, parameter :: krhlev =36
386 !* integer, parameter :: KCM1 = 8 ! num of rh independent aer !species
387 !* integer, parameter :: KCM2 = 5 ! num of rh dependent aer species
388 !* integer, parameter :: KCM = KCM1 + KCM2
390  integer, save :: kcm1 = 0
392  integer, save :: kcm2 = 0
394  integer, save :: kcm
395 
396  real (kind=kind_phys), dimension(KRHLEV) :: rhlev_grt &
397  data rhlev_grt (:)/ .00, .05, .10, .15, .20, .25, .30, .35, &
398  & .40, .45, .50, .55, .60, .65, .70, .75, .80, .81, .82, &
399  & .83, .84, .85, .86, .87, .88, .89, .90, .91, .92, .93, &
400  & .94, .95, .96, .97, .98, .99 /
401 
402 ! --- the following arrays are allocate and setup in subr 'gocrt_aerinit'
403 ! ------ gocart aerosol specification ------
404 ! => transported aerosol species:
405 ! DU (5-bins)
406 ! SS (4 bins for climo mode and 5 bins for fcst mode)
407 ! SU (dms, so2, so4, msa)
408 ! OC (phobic, philic) and BC (phobic, philic)
409 ! => species and lumped species for aerosol optical properties
410 ! DU (5-bins, with 4 sub-groups in the submicron bin )
411 ! SS (ssam for submicron, sscm for coarse mode)
412 ! SU (so4)
413 ! OC (phobic, philic) and BC (phobic, philic)
414 ! => specification used for aerosol optical properties luts
415 ! DU (8 bins)
416 ! SS (ssam, sscm)
417 ! SU (suso)
418 ! OC (waso) and BC (soot)
419 !
420 ! - spectral band structure:
421 ! iendwv_grt(KAERBND) - ending wavenumber (cm**-1) for each band
422 ! - relative humidity independent aerosol optical properties:
423 ! ===> species : dust (8 bins)
424 ! rhidext0_grt(KAERBND,KCM1) - extinction coefficient
425 ! rhidssa0_grt(KAERBND,KCM1) - single scattering albedo
426 ! rhidasy0_grt(KAERBND,KCM1) - asymmetry parameter
427 ! - relative humidity dependent aerosol optical properties:
428 ! ===> species : soot, suso, waso, ssam, sscm
429 ! rhdpext0_grt(KAERBND,KRHLEV,KCM2) - extinction coefficient
430 ! rhdpssa0_grt(KAERBND,KRHLEV,KCM2) - single scattering albedo
431 ! rhdpasy0_grt(KAERBND,KRHLEV,KCM2) - asymmetry parameter
432 
434  integer, allocatable, dimension(:) :: iendwv_grt
435 ! relative humidity independent aerosol optical properties:
436 !! species : dust (8 bins)
437 
440 
442  real (kind=kind_phys),allocatable, dimension(:,:) :: rhidext0_grt
444  real (kind=kind_phys),allocatable, dimension(:,:) :: rhidssa0_grt
446  real (kind=kind_phys), allocatable, dimension(:,:) :: rhidasy0_grt
447 !
448 ! relative humidity dependent aerosol optical properties:
449 ! species : soot, suso, waso, ssam, sscm
450 
453 
455  real (kind=kind_phys),allocatable,dimension(:,:,:) :: rhdpext0_grt
457  real (kind=kind_phys),allocatable,dimension(:,:,:) :: rhdpssa0_grt
459  real (kind=kind_phys),allocatable,dimension(:,:,:) :: rhdpasy0_grt
460 
461 ! - relative humidity independent aerosol optical properties:
462 ! extrhi_grt(KCM1,NSWLWBD) - extinction coefficient for sw+lw spectral band
463 ! ssarhi_grt(KCM1,NSWLWBD) - single scattering albedo for sw+lw spectral band
464 ! asyrhi_grt(KCM1,NSWLWBD) - asymmetry parameter for sw+lw spectral band
465 ! - relative humidity dependent aerosol optical properties:
466 ! extrhd_grt(KRHLEV,KCM2,NSWLWBD) - extinction coefficient for sw+lw band
467 ! ssarhd_grt(KRHLEV,KCM2,NSWLWBD) - single scattering albedo for sw+lw band
468 ! asyrhd_grt(KRHLEV,KCM2,NSWLWBD) - asymmetry parameter for sw+lw band
469 
471 
473  real (kind=kind_phys),allocatable,save,dimension(:,:) :: &
474  & extrhi_grt
475 
476  real (kind=kind_phys),allocatable,save,dimension(:,:) :: &
477  & ssarhi_grt
478 
479  real (kind=kind_phys),allocatable,save,dimension(:,:) :: &
480  & asyrhi_grt
481 
483 
485  real (kind=kind_phys),allocatable,save,dimension(:,:,:) :: &
486  & extrhd_grt
487 
488  real (kind=kind_phys),allocatable,save,dimension(:,:,:) :: &
489  & ssarhd_grt
490 
491  real (kind=kind_phys),allocatable,save,dimension(:,:,:) :: &
492  & asyrhd_grt
493 
495 
496 ! --------------------------------------------------------------------- !
497 ! section-5 : module variables for gocart aerosol climo data set !
498 ! --------------------------------------------------------------------- !
499 ! This version only supports geos3-gocart data set (Jan 2010)
500 ! Modified to support geos4-gocart data set (May 2010)
501 !
502 ! geos3-gocart vs geos4-gocart
503 ! (1) Use the same module variables
504 ! IMXG,JMXG,KMXG,NMXG,psclmg,dmclmg,geos_rlon,geos_rlat
505 ! (2) Similarity between geos3 and geos 4:
506 ! identical lat/lon grids and aerosol specification;
507 ! direction of vertical index is bottom-up (sfc to toa)
508 ! (3) Difference between geos3 and geos4
509 ! vertical coordinate (sigma for geos3/hybrid_sigma_pressure for geos4)
510 ! aerosol units (mass concentration for geos3/mixing ratio for geos4)
511 
513  integer, parameter :: imxg = 144
515  integer, parameter :: jmxg = 91
517  integer, parameter :: kmxg = 30
518 !* integer, parameter :: NMXG = 12
520  integer, save :: nmxg
521 
522  real (kind=kind_phys), parameter :: dltx = 360.0 / float(imxg)
523  real (kind=kind_phys), parameter :: dlty = 180.0 / float(jmxg-1)
524 
525 ! --- the following arrays are allocated and setup in 'rd_gocart_clim'
526 ! - geos-gocart climo data (input dataset)
527 ! psclmg - pressure in cb IMXG*JMXG*KMXG
528 ! dmclmg - aerosol dry mass in g/m3 IMXG*JMXG*KMXG*NMXG
529 ! or aerosol mixing ratio in mol/mol or Kg/Kg
530 
532  real (kind=kind_phys),allocatable, save:: psclmg(:,:,:)
534  real (kind=kind_phys),allocatable, save:: dmclmg(:,:,:,:)
535 
536 ! - geos-gocart lat/lon arrays
538  real (kind=kind_phys), allocatable, save, dimension(:):: geos_rlon
540  real (kind=kind_phys), allocatable, save, dimension(:):: geos_rlat
541 
544  character*4, save :: gocart_climo = 'xxxx'
545 
547  real (kind=kind_io4), allocatable :: molwgt(:)
548 
549 ! ---------------------------------------------------------------------
550 ! !
551 ! section-6 : module variables for gocart aerosol scheme options
552 ! !
553 ! ---------------------------------------------------------------------
554 ! !
555 
557  logical, save :: lgrtint = .true.
558 
560 ! logical, save :: lckprnt = .true.
561  logical, save :: lckprnt = .false.
562 
563 ! --- the following index/flag/weight are set up in 'set_aerspc'
564 
566  real (kind=kind_phys), save :: ctaer = f_zero ! user specified wgt
567 
569  logical, save :: get_fcst = .true.
571  logical, save :: get_clim = .true.
572 
573 ! ------ gocart aerosol specification ------
574 ! => transported aerosol species:
575 ! DU (5-bins)
576 ! SS (4 bins for climo mode and 5 bins for fcst mode)
577 ! SU (dms, so2, so4, msa)
578 ! OC (phobic, philic) and BC (phobic, philic)
579 ! => species and lumped species for aerosol optical properties
580 ! DU (5-bins, with 4 sub-groups in the submicron bin )
581 ! SS (ssam for submicron, sscm for coarse mode)
582 ! SU (so4)
583 ! OC (phobic, philic) and BC (phobic, philic)
584 ! => specification used for aerosol optical properties luts
585 ! DU (8 bins)
586 ! SS (ssam, sscm)
587 ! SU (suso)
588 ! OC (waso) and BC (soot)
589 !
590 
593  integer, save :: isoot, iwaso, isuso, issam, isscm
594 
595 ! - index for rh independent aerosol optical properties (1st
596 ! dimension for extrhi_grt, ssarhi_grt, and asyrhi_grt) is
597 ! not needed ===> hardwired to 8-bin dust
598 
599 ! - index for gocart aerosol species to be included in the
600 ! calculations of aerosol optical properties (ext, ssa, asy)
603  type gocart_index_type
604 ! dust
605  integer :: dust1, dust2, dust3, dust4, dust5
606 ! sea salt
607  integer :: ssam, sscm
608 ! sulfate
609  integer :: suso
610 ! oc
611  integer :: waso_phobic, waso_philic
612 ! bc
613  integer :: soot_phobic, soot_philic
614  endtype
615  type(gocart_index_type), save :: dm_indx
616 
618  type tracer_index_type
619 ! dust
620  integer :: du001, du002, du003, du004, du005
621 ! sea salt
622  integer :: ss001, ss002, ss003, ss004, ss005
623 ! sulfate
624  integer :: so4
625 ! oc
626  integer :: ocphobic, ocphilic
627 ! bc
628  integer :: bcphobic, bcphilic
629  endtype
630  type(tracer_index_type), save :: dmfcs_indx
631 
632 ! - grid components to be included in the aeropt calculations
634  integer, save :: num_gridcomp = 0
636  character, allocatable , save :: gridcomp(:)*2
637 
639  integer, parameter :: max_num_gridcomp = 5
641  character*2 :: max_gridcomp(max_num_gridcomp)
642  data max_gridcomp /'DU', 'BC', 'OC', 'SU', 'SS'/
643 
644 ! GOCART code modification end here (Sarah Lu)
645 ! ------------------------!
646 ! =======================================================================
647 
648 !! --- the following are for diagnostic purpose to output aerosol optical depth
649 ! aod from 10 components are grouped into 5 major different species:
650 ! 1:dust (inso,minm,miam,micm,mitr); 2:black carbon (soot)
651 ! 3:water soluble (waso); 4:sulfate (suso); 5:sea salt (ssam,sscm)
652 !
653 ! idxspc (NCM) - index conversion array
654 ! lspcaod - logical flag for aod from individual species
655 !
657  integer, dimension(NCM) :: idxspc
658  data idxspc / 1, 2, 1, 1, 1, 1, 3, 5, 5, 4 /
659 !
660 ! - wvn550 is the wavenumber (1/cm) of wavelenth 550nm for diagnostic aod output
661 ! nv_aod is the sw spectral band covering wvn550 (comp in aer_init)
662 !
664  real (kind=kind_phys), parameter :: wvn550 = 1.0e4/0.55
666  integer, save :: nv_aod = 1
667 
668 ! --- public interfaces
669 
670  public aer_init, aer_update, setaer
671 
672 
673 ! =================
674  contains
675 ! =================
676 
684 !-----------------------------------
685  subroutine aer_init &
686  & ( nlay, me ) ! --- inputs
687 ! --- outputs: ( to module variables )
688 
689 ! ================================================================== !
690 ! !
691 ! aer_init is the initialization program to set up necessary !
692 ! parameters and working arrays. !
693 ! !
694 ! inputs: !
695 ! NLAY - number of model vertical layers (not used) !
696 ! me - print message control flag !
697 ! !
698 ! outputs: (to module variables) !
699 ! !
700 ! external module variables: (in physparam) !
701 ! iaermdl - tropospheric aerosol model scheme flag !
702 ! =0 opac-clim; =1 gocart-clim, =2 gocart-prognostic !
703 ! lalwflg - logical lw aerosols effect control flag !
704 ! =t compute lw aerosol optical prop !
705 ! laswflg - logical sw aerosols effect control flag !
706 ! =t compute sw aerosol optical prop !
707 ! lavoflg - logical stratosphere volcanic aerosol control flag !
708 ! =t include volcanic aerosol effect !
709 ! lalw1bd = logical lw aeros propty 1 band vs multi-band cntl flag !
710 ! =t use 1 broad band optical property !
711 ! =f use multi bands optical property !
712 ! !
713 ! module constants: !
714 ! NWVSOL - num of wvnum regions where solar flux is constant !
715 ! NWVTOT - total num of wave numbers used in sw spectrum !
716 ! NWVTIR - total num of wave numbers used in the ir region !
717 ! NSWBND - total number of sw spectral bands !
718 ! NLWBND - total number of lw spectral bands !
719 ! !
720 ! usage: call aer_init !
721 ! !
722 ! subprograms called: clim_aerinit, gcrt_aerinit, !
723 ! wrt_aerlog, set_volcaer, set_spectrum, !
724 ! !
725 ! ================================================================== !
726 
727 ! --- inputs:
728  integer, intent(in) :: NLAY, me
729 
730 ! --- output: ( none )
731 
732 ! --- locals:
733  real (kind=kind_phys), dimension(NWVTOT) :: solfwv ! one wvn sol flux
734  real (kind=kind_phys), dimension(NWVTIR) :: eirfwv ! one wvn ir flux
735 !
736 !===> ... begin here
737 !
738  kyrstr = 1
739  kyrend = 1
740  kyrsav = 1
741  kmonsav = 1
742 
744 
745  if ( me == 0 ) then
746 
747  call wrt_aerlog ! write aerosol param info to log file
748 ! --- inputs: (in scope variables)
749 ! --- outputs: ( none )
750 
751  endif
752 
753  if ( iaerflg == 0 ) return ! return without any aerosol calculations
754 
755 ! --- ... in sw, aerosols optical properties are computed for each radiation
756 ! spectral band; while in lw, optical properties can be calculated
757 ! for either only one broad band or for each of the lw radiation bands
758 
759  if ( laswflg ) then
760  nswbnd = nbdsw
761  else
762  nswbnd = 0
763  endif
764 
765  if ( lalwflg ) then
766  if ( lalw1bd ) then
767  nlwbnd = 1
768  else
769  nlwbnd = nbdlw
770  endif
771  else
772  nlwbnd = 0
773  endif
774 
775  nswlwbd = nswbnd + nlwbnd
776 
777  if ( iaerflg /= 100 ) then
778 
781 
782  call set_spectrum
783 ! --- inputs: (module constants)
784 ! --- outputs: (in-scope variables)
785 
787 
788  if ( iaermdl == 0 ) then ! opac-climatology scheme
789 
790  call clim_aerinit &
791 ! --- inputs:
792  & ( solfwv, eirfwv, me &
793 ! --- outputs:
794  & )
795 
796 ! elseif ( iaermdl == 1 ) then ! gocart-climatology scheme
797 ! elseif ( iaermdl==1 .or. iaermdl==2 ) then ! gocart-clim/prog scheme
798 
799 ! call gcrt_climinit
800 
801 ! elseif ( iaermdl == 2 ) then ! gocart-prognostic scheme
802 
803 ! call gcrt_aerinit
804 
805  else
806  if ( me == 0 ) then
807  print *,' !!! ERROR in aerosol model scheme selection', &
808  & ' iaermdl =',iaermdl
809  stop
810  endif
811  endif
812 
813  endif ! end if_iaerflg_block
814 
817 
818  if ( lavoflg ) then
819 
820  call set_volcaer
821 ! --- inputs: (module variables)
822 ! --- outputs: (module variables)
823 
824  endif ! end if_lavoflg_block
825 
826 
827 ! =================
828  contains
829 ! =================
830 
832 !--------------------------------
833  subroutine wrt_aerlog
834 !................................
835 ! --- inputs: (in scope variables)
836 ! --- outputs: ( none )
837 
838 ! ================================================================== !
839 ! !
840 ! subprogram : wrt_aerlog !
841 ! !
842 ! write aerosol parameter configuration to run log file. !
843 ! !
844 ! ==================== defination of variables =================== !
845 ! !
846 ! external module variables: (in physparam) !
847 ! iaermdl - aerosol scheme flag: 0:opac-clm; 1:gocart-clim; !
848 ! 2:gocart-prog !
849 ! iaerflg - aerosol effect control flag: 3-digits (volc,lw,sw) !
850 ! lalwflg - toposphere lw aerosol effect: =f:no; =t:yes !
851 ! laswflg - toposphere sw aerosol effect: =f:no; =t:yes !
852 ! lavoflg - stratospherer volcanic aeros effect: =f:no; =t:yes !
853 ! !
854 ! outputs: ( none ) !
855 ! !
856 ! subroutines called: none !
857 ! !
858 ! usage: call wrt_aerlog !
859 ! !
860 ! ================================================================== !
861 
862 ! --- inputs: ( none )
863 ! --- output: ( none )
864 ! --- locals:
865 
866 !
867 !===> ... begin here
868 !
869  print *, vtagaer ! print out version tag
870 
871  if ( iaermdl == 0 ) then
872  print *,' - Using OPAC-seasonal climatology for tropospheric', &
873  & ' aerosol effect'
874  elseif ( iaermdl == 1 ) then
875  print *,' - Using GOCART-climatology for tropospheric', &
876  & ' aerosol effect'
877  elseif ( iaermdl == 2 ) then
878  print *,' - Using GOCART-prognostic aerosols for tropospheric', &
879  & ' aerosol effect'
880  else
881  print *,' !!! ERROR in selection of aerosol model scheme', &
882  & ' IAER_MDL =',iaermdl
883  stop
884  endif ! end_if_iaermdl_block
885 
886  print *,' IAER=',iaerflg,' LW-trop-aer=',lalwflg, &
887  & ' SW-trop-aer=',laswflg,' Volc-aer=',lavoflg
888 
889  if ( iaerflg <= 0 ) then ! turn off all aerosol effects
890  print *,' - No tropospheric/volcanic aerosol effect included'
891  print *,' Input values of aerosol optical properties to' &
892  & ,' both SW and LW radiations are set to zeros'
893  else
894  if ( iaerflg >= 100 ) then ! incl stratospheric volcanic aerosols
895  print *,' - Include stratospheric volcanic aerosol effect'
896  else ! no stratospheric volcanic aerosols
897  print *,' - No stratospheric volcanic aerosol effect'
898  endif
899 
900  if ( laswflg ) then ! chcek for sw effect
901  print *,' - Compute multi-band aerosol optical' &
902  & ,' properties for SW input parameters'
903  else
904  print *,' - No SW radiation aerosol effect, values of' &
905  & ,' aerosol properties to SW input are set to zeros'
906  endif
907 
908  if ( lalwflg ) then ! check for lw effect
909  if ( lalw1bd ) then
910  print *,' - Compute 1 broad-band aerosol optical' &
911  & ,' properties for LW input parameters'
912  else
913  print *,' - Compute multi-band aerosol optical' &
914  & ,' properties for LW input parameters'
915  endif
916  else
917  print *,' - No LW radiation aerosol effect, values of' &
918  & ,' aerosol properties to LW input are set to zeros'
919  endif
920  endif ! end if_iaerflg_block
921 !
922  return
923 !................................
924  end subroutine wrt_aerlog
925 !--------------------------------
926 
931 !--------------------------------
932  subroutine set_spectrum
933 !................................
934 ! --- inputs: (module constants)
935 ! --- outputs: (in-scope variables)
936 
937 ! ================================================================== !
938 ! !
939 ! subprogram : set_spectrum !
940 ! !
941 ! define the one wavenumber solar fluxes based on toa solar spectral!
942 ! distrobution, and define the one wavenumber ir fluxes based on !
943 ! black-body emission distribution at a predefined temperature. !
944 ! !
945 ! ==================== defination of variables =================== !
946 ! !
956 ! !
957 ! subroutines called: none !
958 ! !
959 ! usage: call set_spectrum !
960 ! !
961 ! ================================================================== !
962 
963 ! --- inputs: (module constants)
964 ! integer :: NWVTOT, NWVTIR
965 
966 ! --- output: (in-scope variables)
967 ! real (kind=kind_phys), dimension(NWVTOT) :: solfwv ! one wvn sol flux
968 ! real (kind=kind_phys), dimension(NWVTIR) :: eirfwv ! one wvn ir flux
969 
970 ! --- locals:
971  real (kind=kind_phys) :: soltot, tmp1, tmp2, tmp3
972 
973  integer :: nb, nw, nw1, nw2, nmax, nmin
974 !
975 !===> ... begin here
976 !
977 ! nmax = min( NWVTOT, nint( maxval(wvnsw2) ))
978 ! nmin = max( 1, nint( minval(wvnsw1) ))
979 
980 ! --- check print
981 ! print *,' MINWVN, MAXWVN = ',nmin, nmax
982 ! --- ... define the one wavenumber solar fluxes based on toa solar
983 ! spectral distribution
984 
985 ! soltot1 = f_zero
986 ! soltot = f_zero
987  do nb = 1, nwvsol
988  if ( nb == 1 ) then
989  nw1 = 1
990  else
991  nw1 = nw1 + nwvns0(nb-1)
992  endif
993 
994  nw2 = nw1 + nwvns0(nb) - 1
995 
996  do nw = nw1, nw2
997  solfwv(nw) = s0intv(nb)
998 ! soltot1 = soltot1 + s0intv(nb)
999 ! if ( nw >= nmin .and. nw <= nmax ) then
1000 ! soltot = soltot + s0intv(nb)
1001 ! endif
1002  enddo
1003  enddo
1004 
1005 ! --- ... define the one wavenumber ir fluxes based on black-body
1006 ! emission distribution at a predefined temperature
1007 
1008  tmp1 = (con_pi + con_pi) * con_plnk * con_c* con_c
1009  tmp2 = con_plnk * con_c / (con_boltz * con_t0c)
1010 
1011 !$omp parallel do private(nw,tmp3)
1012  do nw = 1, nwvtir
1013  tmp3 = 100.0 * nw
1014  eirfwv(nw) = (tmp1 * tmp3**3) / (exp(tmp2*tmp3) - 1.0)
1015  enddo
1016 !
1017  return
1018 !................................
1019  end subroutine set_spectrum
1020 !--------------------------------
1021 
1022 
1024 !-----------------------------
1025  subroutine set_volcaer
1026 !.............................
1027 ! --- inputs: ( none )
1028 ! --- outputs: (module variables)
1029 
1030 ! ================================================================== !
1031 ! !
1032 ! subprogram : set_volcaer !
1033 ! !
1034 ! this is the initialization progrmam for stratospheric volcanic !
1035 ! aerosols. !
1036 ! !
1037 ! subroutines called: none !
1038 ! !
1039 ! usage: call set_volcaer !
1040 ! !
1041 ! ================================================================== !
1042 
1043 ! --- inputs: (none)
1044 
1045 ! --- output: (module variables)
1046 ! integer :: ivolae(:,:,:)
1047 
1048 ! --- locals:
1049 !
1050 !===> ... begin here
1051 !
1052 ! --- allocate data space
1053 
1054  if ( .not. allocated(ivolae) ) then
1055  allocate ( ivolae(12,4,10) ) ! for 12-mon,4-lat_zone,10-year
1056  endif
1057 !
1058  return
1059 !................................
1060  end subroutine set_volcaer
1061 !--------------------------------
1062 !
1063 !...................................
1064  end subroutine aer_init
1065 !-----------------------------------
1066 !!@}
1067 
1068 
1079 !-----------------------------------
1080  subroutine clim_aerinit &
1081  & ( solfwv, eirfwv, me & ! --- inputs
1082  & ) ! --- outputs
1083 
1084 ! ================================================================== !
1085 ! !
1086 ! clim_aerinit is the opac-climatology aerosol initialization program !
1087 ! to set up necessary parameters and working arrays. !
1088 ! !
1089 ! inputs: !
1090 ! solfwv(NWVTOT) - solar flux for each individual wavenumber (w/m2)!
1091 ! eirfwv(NWVTIR) - ir flux(273k) for each individual wavenum (w/m2)!
1092 ! me - print message control flag !
1093 ! !
1094 ! outputs: (to module variables) !
1095 ! !
1096 ! external module variables: (in physparam) !
1097 ! iaerflg - abc 3-digit integer aerosol flag (abc:volc,lw,sw) !
1098 ! a: =0 use background stratospheric aerosol !
1099 ! =1 incl stratospheric vocanic aeros (MINVYR-MAXVYR) !
1100 ! b: =0 no topospheric aerosol in lw radiation !
1101 ! =1 include tropspheric aerosols for lw radiation !
1102 ! c: =0 no topospheric aerosol in sw radiation !
1103 ! =1 include tropspheric aerosols for sw radiation !
1104 ! lalwflg - logical lw aerosols effect control flag !
1105 ! =t compute lw aerosol optical prop !
1106 ! laswflg - logical sw aerosols effect control flag !
1107 ! =t compute sw aerosol optical prop !
1108 ! lalw1bd = logical lw aeros propty 1 band vs multi-band cntl flag !
1109 ! =t use 1 broad band optical property !
1110 ! =f use multi bands optical property !
1111 ! !
1112 ! module constants: !
1113 ! NWVSOL - num of wvnum regions where solar flux is constant !
1114 ! NWVTOT - total num of wave numbers used in sw spectrum !
1115 ! NWVTIR - total num of wave numbers used in the ir region !
1116 ! NSWBND - total number of sw spectral bands !
1117 ! NLWBND - total number of lw spectral bands !
1118 ! NAERBND - number of bands for climatology aerosol data !
1119 ! NCM1 - number of rh independent aeros species !
1120 ! NCM2 - number of rh dependent aeros species !
1121 ! !
1122 ! usage: call clim_aerinit !
1123 ! !
1124 ! subprograms called: set_aercoef, optavg !
1125 ! !
1126 ! ================================================================== !
1127 
1128 ! --- inputs:
1129  real (kind=kind_phys), dimension(:) :: solfwv ! one wvn sol flux
1130  real (kind=kind_phys), dimension(:) :: eirfwv ! one wvn ir flux
1131 
1132  integer, intent(in) :: me
1133 
1134 ! --- output: ( none )
1135 
1136 ! --- locals:
1137  real (kind=kind_phys), dimension(NAERBND,NCM1) :: &
1138  & rhidext0, rhidsca0, rhidssa0, rhidasy0
1139  real (kind=kind_phys), dimension(NAERBND,NRHLEV,NCM2):: &
1140  & rhdpext0, rhdpsca0, rhdpssa0, rhdpasy0
1141  real (kind=kind_phys), dimension(NAERBND) :: straext0
1142 
1143  real (kind=kind_phys), dimension(NSWBND,NAERBND) :: solwaer
1144  real (kind=kind_phys), dimension(NSWBND) :: solbnd
1145  real (kind=kind_phys), dimension(NLWBND,NAERBND) :: eirwaer
1146  real (kind=kind_phys), dimension(NLWBND) :: eirbnd
1147 
1148  integer, dimension(NSWBND) :: nv1, nv2
1149  integer, dimension(NLWBND) :: nr1, nr2
1150 !
1151 !===> ... begin here
1152 !
1153 ! --- ... invoke tropospheric aerosol initialization
1154 
1156  call set_aercoef
1157 ! --- inputs: (in-scope variables, module constants)
1158 ! --- outputs: (module variables)
1159 
1160 
1161 ! =================
1162  contains
1163 ! =================
1164 
1170 !--------------------------------
1171  subroutine set_aercoef
1172 !................................
1173 ! --- inputs: (in-scope variables, module constants)
1174 ! --- outputs: (module variables)
1175 
1176 ! ================================================================== !
1177 ! !
1178 ! subprogram : set_aercoef !
1179 ! !
1180 ! this is the initialization progrmam for climatological aerosols !
1181 ! !
1182 ! the program reads and maps the pre-tabulated aerosol optical !
1183 ! spectral data onto corresponding sw radiation spectral bands. !
1184 ! !
1185 ! ==================== defination of variables =================== !
1186 ! !
1187 ! inputs: (in-scope variables, module constants) !
1188 ! solfwv(:) - real, solar flux for individual wavenumber (w/m2) !
1189 ! eirfwv(:) - real, lw flux(273k) for individual wavenum (w/m2) !
1190 ! me - integer, select cpu number as print control flag !
1191 ! !
1192 ! outputs: (to the module variables) !
1193 ! !
1194 ! external module variables: (in physparam) !
1195 ! lalwflg - module control flag for lw trop-aer: =f:no; =t:yes !
1196 ! laswflg - module control flag for sw trop-aer: =f:no; =t:yes !
1197 ! aeros_file- external aerosol data file name !
1198 ! !
1199 ! internal module variables: !
1200 ! IMXAE - number of longitude points in global aeros data set !
1201 ! JMXAE - number of latitude points in global aeros data set !
1202 ! wvnsw1,wvnsw2 (NSWSTR:NSWEND) !
1203 ! - start/end wavenumbers for each of sw bands !
1204 ! wvnlw1,wvnlw2 ( 1:NBDLW) !
1205 ! - start/end wavenumbers for each of lw bands !
1206 ! NSWLWBD - total num of bands (sw+lw) for aeros optical properties!
1207 ! NSWBND - number of sw spectral bands actually invloved !
1208 ! NLWBND - number of lw spectral bands actually invloved !
1209 ! NIAERCM - unit number for reading input data set !
1210 ! extrhi - extinction coef for rh-indep aeros NCM1*NSWLWBD!
1211 ! scarhi - scattering coef for rh-indep aeros NCM1*NSWLWBD!
1212 ! ssarhi - single-scat-alb for rh-indep aeros NCM1*NSWLWBD!
1213 ! asyrhi - asymmetry factor for rh-indep aeros NCM1*NSWLWBD!
1214 ! extrhd - extinction coef for rh-dep aeros NRHLEV*NCM2*NSWLWBD!
1215 ! scarhd - scattering coef for rh-dep aeros NRHLEV*NCM2*NSWLWBD!
1216 ! ssarhd - single-scat-alb for rh-dep aeros NRHLEV*NCM2*NSWLWBD!
1217 ! asyrhd - asymmetry factor for rh-dep aeros NRHLEV*NCM2*NSWLWBD!
1218 ! !
1219 ! major local variables: !
1220 ! for handling spectral band structures !
1221 ! iendwv - ending wvnum (cm**-1) for each band NAERBND !
1222 ! for handling optical properties of rh independent species (NCM1) !
1223 ! 1. insoluble (inso); 2. soot (soot); !
1224 ! 3. mineral nuc mode (minm); 4. mineral acc mode (miam); !
1225 ! 5. mineral coa mode (micm); 6. mineral transport(mitr). !
1226 ! rhidext0 - extinction coefficient NAERBND*NCM1 !
1227 ! rhidsca0 - scattering coefficient NAERBND*NCM1 !
1228 ! rhidssa0 - single scattering albedo NAERBND*NCM1 !
1229 ! rhidasy0 - asymmetry parameter NAERBND*NCM1 !
1230 ! for handling optical properties of rh ndependent species (NCM2) !
1231 ! 1. water soluble (waso); 2. sea salt acc mode(ssam); !
1232 ! 3. sea salt coa mode(sscm); 4. sulfate droplets (suso). !
1233 ! rh level (NRHLEV): 00%, 50%, 70%, 80%, 90%, 95%, 98%, 99% !
1234 ! rhdpext0 - extinction coefficient NAERBND,NRHLEV,NCM2!
1235 ! rhdpsca0 - scattering coefficient NAERBND,NRHLEV,NCM2!
1236 ! rhdpssa0 - single scattering albedo NAERBND,NRHLEV,NCM2!
1237 ! rhdpasy0 - asymmetry parameter NAERBND,NRHLEV,NCM2!
1238 ! for handling optical properties of stratospheric bkgrnd aerosols !
1239 ! straext0 - extingction coefficients NAERBND !
1240 ! !
1241 ! usage: call set_aercoef !
1242 ! !
1243 ! subprograms called: optavg !
1244 ! !
1245 ! ================================================================== !
1246 !
1247 ! --- inputs: ( none )
1248 ! --- output: ( none )
1249 
1250 ! --- locals:
1251  integer, dimension(NAERBND) :: iendwv
1252 
1253  integer :: i, j, k, m, mb, ib, ii, id, iw, iw1, iw2
1254 
1255  real (kind=kind_phys) :: sumsol, sumir
1256 
1257  logical :: file_exist
1258  character :: cline*80
1259 !
1260 !===> ... begin here
1261 !
1264 
1265  inquire (file=aeros_file, exist=file_exist)
1266 
1267  if ( file_exist ) then
1268  close (niaercm)
1269  open (unit=niaercm,file=aeros_file,status='OLD', &
1270  & action='read',form='FORMATTED')
1271  rewind(niaercm)
1272  else
1273  print *,' Requested aerosol data file "',aeros_file, &
1274  & '" not found!'
1275  print *,' *** Stopped in subroutine aero_init !!'
1276  stop
1277  endif ! end if_file_exist_block
1278 
1279 ! --- ... skip monthly global distribution
1280 
1281  do m = 1, 12
1282  read (niaercm,12) cline
1283  12 format(a80/)
1284 
1285  do j = 1, jmxae
1286  do i = 1, imxae
1287  read(niaercm,*) id
1288  enddo
1289  enddo
1290  enddo ! end do_m_block
1291 
1292 ! --- ... aloocate and input aerosol optical data
1293 
1294  if ( .not. allocated( extrhi ) ) then
1295  allocate ( extrhi( ncm1,nswlwbd) )
1296  allocate ( scarhi( ncm1,nswlwbd) )
1297  allocate ( ssarhi( ncm1,nswlwbd) )
1298  allocate ( asyrhi( ncm1,nswlwbd) )
1299  allocate ( extrhd(nrhlev,ncm2,nswlwbd) )
1300  allocate ( scarhd(nrhlev,ncm2,nswlwbd) )
1301  allocate ( ssarhd(nrhlev,ncm2,nswlwbd) )
1302  allocate ( asyrhd(nrhlev,ncm2,nswlwbd) )
1303  allocate ( extstra( nswlwbd) )
1304  endif
1305 
1307  read(niaercm,21) cline
1308  21 format(a80)
1309  read(niaercm,22) iendwv(:)
1310  22 format(13i6)
1311 
1313  read(niaercm,21) cline
1314  read(niaercm,24) haer(:,:)
1315  24 format(20f4.1)
1316 
1318  read(niaercm,21) cline
1319  read(niaercm,26) prsref(:,:)
1320  26 format(10f7.2)
1321 
1323  read(niaercm,21) cline
1324  read(niaercm,28) rhidext0(:,:)
1325  28 format(8e10.3)
1326 
1328  read(niaercm,21) cline
1329  read(niaercm,28) rhidsca0(:,:)
1330 
1332  read(niaercm,21) cline
1333  read(niaercm,28) rhidssa0(:,:)
1334 
1336  read(niaercm,21) cline
1337  read(niaercm,28) rhidasy0(:,:)
1338 
1340  read(niaercm,21) cline
1341  read(niaercm,28) rhdpext0(:,:,:)
1342 
1344  read(niaercm,21) cline
1345  read(niaercm,28) rhdpsca0(:,:,:)
1346 
1348  read(niaercm,21) cline
1349  read(niaercm,28) rhdpssa0(:,:,:)
1350 
1352  read(niaercm,21) cline
1353  read(niaercm,28) rhdpasy0(:,:,:)
1354 
1356  read(niaercm,21) cline
1357  read(niaercm,28) straext0(:)
1358 
1359  close (niaercm)
1360 
1363 
1364  sigref(:,:) = 0.001 * prsref(:,:)
1365 
1368 
1369  if ( laswflg ) then
1370  solbnd(:) = f_zero
1371 !$omp parallel do private(i,j)
1372  do j=1,naerbnd
1373  do i=1,nswbnd
1374  solwaer(i,j) = f_zero
1375  enddo
1376  enddo
1377 
1378 !$omp parallel do private(ib,mb,ii,iw1,iw2,iw,nv_aod,sumsol)
1379  do ib = 1, nswbnd
1380  mb = ib + nswstr - 1
1381  ii = 1
1382  iw1 = nint(wvnsw1(mb))
1383  iw2 = nint(wvnsw2(mb))
1384 
1385  if ( wvnsw2(mb) >= wvn550 .and. wvn550 >= wvnsw1(mb) ) then
1386  nv_aod = ib ! sw band number covering 550nm wavelenth
1387  endif
1388 
1389  lab_swdowhile : do while ( iw1 > iendwv(ii) )
1390  if ( ii == naerbnd ) exit lab_swdowhile
1391  ii = ii + 1
1392  enddo lab_swdowhile
1393 
1394  sumsol = f_zero
1395  nv1(ib) = ii
1396 
1397  do iw = iw1, iw2
1398  solbnd(ib) = solbnd(ib) + solfwv(iw)
1399  sumsol = sumsol + solfwv(iw)
1400 
1401  if ( iw == iendwv(ii) ) then
1402  solwaer(ib,ii) = sumsol
1403 
1404  if ( ii < naerbnd ) then
1405  sumsol = f_zero
1406  ii = ii + 1
1407  endif
1408  endif
1409  enddo
1410 
1411  if ( iw2 /= iendwv(ii) ) then
1412  solwaer(ib,ii) = sumsol
1413  endif
1414 
1415  nv2(ib) = ii
1416 ! frcbnd(ib) = solbnd(ib) / soltot
1417  enddo ! end do_ib_block for sw
1418  endif ! end if_laswflg_block
1419 
1422 
1423  if ( lalwflg ) then
1424  eirbnd(:) = f_zero
1425 !$omp parallel do private(i,j)
1426  do j=1,naerbnd
1427  do i=1,nlwbnd
1428  eirwaer(i,j) = f_zero
1429  enddo
1430  enddo
1431 
1432 !$omp parallel do private(ib,ii,iw1,iw2,iw,mb,sumir)
1433  do ib = 1, nlwbnd
1434  ii = 1
1435  if ( nlwbnd == 1 ) then
1436 ! iw1 = 250 ! corresponding 40 mu
1437  iw1 = 400 ! corresponding 25 mu
1438  iw2 = 2500 ! corresponding 4 mu
1439  else
1440  mb = ib + nlwstr - 1
1441  iw1 = nint(wvnlw1(mb))
1442  iw2 = nint(wvnlw2(mb))
1443  endif
1444 
1445  lab_lwdowhile : do while ( iw1 > iendwv(ii) )
1446  if ( ii == naerbnd ) exit lab_lwdowhile
1447  ii = ii + 1
1448  enddo lab_lwdowhile
1449 
1450  sumir = f_zero
1451  nr1(ib) = ii
1452 
1453  do iw = iw1, iw2
1454  eirbnd(ib) = eirbnd(ib) + eirfwv(iw)
1455  sumir = sumir + eirfwv(iw)
1456 
1457  if ( iw == iendwv(ii) ) then
1458  eirwaer(ib,ii) = sumir
1459 
1460  if ( ii < naerbnd ) then
1461  sumir = f_zero
1462  ii = ii + 1
1463  endif
1464  endif
1465  enddo
1466 
1467  if ( iw2 /= iendwv(ii) ) then
1468  eirwaer(ib,ii) = sumir
1469  endif
1470 
1471  nr2(ib) = ii
1472  enddo ! end do_ib_block for lw
1473  endif ! end if_lalwflg_block
1474 
1477 
1478  call optavg
1479 ! --- inputs: (in-scope variables, module variables)
1480 ! --- outputs: (module variables)
1481 
1482 ! --- check print
1483 ! do ib = 1, NSWBND
1484 ! print *,' After optavg, for sw band:',ib
1485 ! print *,' extrhi:', extrhi(:,ib)
1486 ! print *,' scarhi:', scarhi(:,ib)
1487 ! print *,' ssarhi:', ssarhi(:,ib)
1488 ! print *,' asyrhi:', asyrhi(:,ib)
1489 ! mb = ib + NSWSTR - 1
1490 ! print *,' wvnsw1,wvnsw2 :',wvnsw1(mb),wvnsw2(mb)
1491 ! do i = 1, NRHLEV
1492 ! print *,' extrhd for rhlev:',i
1493 ! print *,extrhd(i,:,ib)
1494 ! print *,' scarhd for rhlev:',i
1495 ! print *,scarhd(i,:,ib)
1496 ! print *,' ssarhd for rhlev:',i
1497 ! print *,ssarhd(i,:,ib)
1498 ! print *,' asyrhd for rhlev:',i
1499 ! print *,asyrhd(i,:,ib)
1500 ! enddo
1501 ! print *,' extstra:', extstra(ib)
1502 ! enddo
1503 ! print *,' wvnlw1 :',wvnlw1
1504 ! print *,' wvnlw2 :',wvnlw2
1505 ! do ib = 1, NLWBND
1506 ! ii = NSWBND + ib
1507 ! print *,' After optavg, for lw band:',ib
1508 ! print *,' extrhi:', extrhi(:,ii)
1509 ! print *,' scarhi:', scarhi(:,ii)
1510 ! print *,' ssarhi:', ssarhi(:,ii)
1511 ! print *,' asyrhi:', asyrhi(:,ii)
1512 ! do i = 1, NRHLEV
1513 ! print *,' extrhd for rhlev:',i
1514 ! print *,extrhd(i,:,ii)
1515 ! print *,' scarhd for rhlev:',i
1516 ! print *,scarhd(i,:,ii)
1517 ! print *,' ssarhd for rhlev:',i
1518 ! print *,ssarhd(i,:,ii)
1519 ! print *,' asyrhd for rhlev:',i
1520 ! print *,asyrhd(i,:,ii)
1521 ! enddo
1522 ! print *,' extstra:', extstra(ii)
1523 ! enddo
1524 !
1525  return
1526 !................................
1527  end subroutine set_aercoef
1528 !--------------------------------
1529 !! @}
1530 
1535 !--------------------------------
1536  subroutine optavg
1537 !................................
1538 ! --- inputs: (in-scope variables, module variables
1539 ! --- outputs: (module variables)
1540 
1541 ! ==================================================================== !
1542 ! !
1543 ! subprogram: optavg !
1544 ! !
1545 ! compute mean aerosols optical properties over each sw radiation !
1546 ! spectral band for each of the species components. This program !
1547 ! follows gfdl's approach for thick cloud opertical property in !
1548 ! sw radiation scheme (2000). !
1549 ! !
1550 ! ==================== defination of variables =================== !
1551 ! !
1552 ! major input variables: !
1553 ! nv1,nv2 (NSWBND) - start/end spectral band indices of aerosol data !
1554 ! for each sw radiation spectral band !
1555 ! nr1,nr2 (NLWBND) - start/end spectral band indices of aerosol data !
1556 ! for each ir radiation spectral band !
1557 ! solwaer (NSWBND,NAERBND) !
1558 ! - solar flux weight over each sw radiation band !
1559 ! vs each aerosol data spectral band !
1560 ! eirwaer (NLWBND,NAERBND) !
1561 ! - ir flux weight over each lw radiation band !
1562 ! vs each aerosol data spectral band !
1563 ! solbnd (NSWBND) - solar flux weight over each sw radiation band !
1564 ! eirbnd (NLWBND) - ir flux weight over each lw radiation band !
1565 ! NSWBND - total number of sw spectral bands !
1566 ! NLWBND - total number of lw spectral bands !
1567 ! !
1568 ! external module variables: (in physparam) !
1569 ! laswflg - control flag for sw spectral region !
1570 ! lalwflg - control flag for lw spectral region !
1571 ! !
1572 ! output variables: (to module variables) !
1573 ! !
1574 ! ================================================================== !
1575 
1576 ! --- inputs:
1577 ! --- output:
1578 
1579 ! --- locals:
1580  real (kind=kind_phys) :: sumk, sums, sumok, sumokg, sumreft, &
1581  & sp, refb, reft, rsolbd, rirbd
1582 
1583  integer :: ib, nb, ni, nh, nc
1584 !
1585 !===> ... begin here
1586 !
1587 ! --- ... loop for each sw radiation spectral band
1588 
1589  if ( laswflg ) then
1590 
1591 !$omp parallel do private(nb,nc,sumk,sums,sumok,sumokg,sumreft)
1592 !$omp+private(ni,nh,sp,reft,refb,rsolbd)
1593  do nb = 1, nswbnd
1594  rsolbd = f_one / solbnd(nb)
1595 
1596 ! --- for rh independent aerosol species
1597 
1598  do nc = 1, ncm1 ! --- for rh independent aerosol species
1599  sumk = f_zero
1600  sums = f_zero
1601  sumok = f_zero
1602  sumokg = f_zero
1603  sumreft = f_zero
1604 
1605  do ni = nv1(nb), nv2(nb)
1606  sp = sqrt( (f_one - rhidssa0(ni,nc)) &
1607  & / (f_one - rhidssa0(ni,nc)*rhidasy0(ni,nc)) )
1608  reft = (f_one - sp) / (f_one + sp)
1609  sumreft = sumreft + reft*solwaer(nb,ni)
1610 
1611  sumk = sumk + rhidext0(ni,nc)*solwaer(nb,ni)
1612  sums = sums + rhidsca0(ni,nc)*solwaer(nb,ni)
1613  sumok = sumok + rhidssa0(ni,nc)*solwaer(nb,ni) &
1614  & * rhidext0(ni,nc)
1615  sumokg = sumokg + rhidssa0(ni,nc)*solwaer(nb,ni) &
1616  & * rhidext0(ni,nc)*rhidasy0(ni,nc)
1617  enddo
1618 
1619  refb = sumreft * rsolbd
1620 
1621  extrhi(nc,nb) = sumk * rsolbd
1622  scarhi(nc,nb) = sums * rsolbd
1623  asyrhi(nc,nb) = sumokg / (sumok + 1.0e-10)
1624  ssarhi(nc,nb) = 4.0*refb &
1625  & / ( (f_one+refb)**2 - asyrhi(nc,nb)*(f_one-refb)**2 )
1626  enddo ! end do_nc_block for rh-ind aeros
1627 
1628 
1629  do nc = 1, ncm2 ! --- for rh dependent aerosols species
1630  do nh = 1, nrhlev
1631  sumk = f_zero
1632  sums = f_zero
1633  sumok = f_zero
1634  sumokg = f_zero
1635  sumreft = f_zero
1636 
1637  do ni = nv1(nb), nv2(nb)
1638  sp = sqrt( (f_one - rhdpssa0(ni,nh,nc)) &
1639  & / (f_one - rhdpssa0(ni,nh,nc)*rhdpasy0(ni,nh,nc)) )
1640  reft = (f_one - sp) / (f_one + sp)
1641  sumreft = sumreft + reft*solwaer(nb,ni)
1642 
1643  sumk = sumk + rhdpext0(ni,nh,nc)*solwaer(nb,ni)
1644  sums = sums + rhdpsca0(ni,nh,nc)*solwaer(nb,ni)
1645  sumok = sumok + rhdpssa0(ni,nh,nc)*solwaer(nb,ni) &
1646  & * rhdpext0(ni,nh,nc)
1647  sumokg = sumokg + rhdpssa0(ni,nh,nc)*solwaer(nb,ni) &
1648  & * rhdpext0(ni,nh,nc)*rhdpasy0(ni,nh,nc)
1649  enddo
1650 
1651  refb = sumreft * rsolbd
1652 
1653  extrhd(nh,nc,nb) = sumk * rsolbd
1654  scarhd(nh,nc,nb) = sums * rsolbd
1655  asyrhd(nh,nc,nb) = sumokg / (sumok + 1.0e-10)
1656  ssarhd(nh,nc,nb) = 4.0*refb &
1657  & / ( (f_one+refb)**2 - asyrhd(nh,nc,nb)*(f_one-refb)**2 )
1658  enddo ! end do_nh_block
1659  enddo ! end do_nc_block for rh-dep aeros
1660 
1661 ! --- for stratospheric background aerosols
1662 
1663  sumk = f_zero
1664  do ni = nv1(nb), nv2(nb)
1665  sumk = sumk + straext0(ni)*solwaer(nb,ni)
1666  enddo
1667 
1668  extstra(nb) = sumk * rsolbd
1669 
1670 ! --- check print
1671 ! if ( nb > 6 .and. nb < 10) then
1672 ! print *,' in optavg for sw band',nb
1673 ! print *,' nv1, nv2:',nv1(nb),nv2(nb)
1674 ! print *,' solwaer:',solwaer(nb,nv1(nb):nv2(nb))
1675 ! print *,' extrhi:', extrhi(:,nb)
1676 ! do i = 1, NRHLEV
1677 ! print *,' extrhd for rhlev:',i
1678 ! print *,extrhd(i,:,nb)
1679 ! enddo
1680 ! print *,' sumk, rsolbd, extstra:',sumk,rsolbd,extstra(nb)
1681 ! endif
1682 
1683  enddo ! end do_nb_block for sw
1684  endif ! end if_laswflg_block
1685 
1686 ! --- ... loop for each lw radiation spectral band
1687 
1688  if ( lalwflg ) then
1689 
1690 !$omp parallel do private(nb,ib,nc,rirbd,sumk,sums,sumok,sumokg,sumreft)
1691 !$omp+private(ni,nh,sp,reft,refb,rsolbd)
1692  do nb = 1, nlwbnd
1693 
1694  ib = nswbnd + nb
1695  rirbd = f_one / eirbnd(nb)
1696 
1697  do nc = 1, ncm1 ! --- for rh independent aerosol species
1698  sumk = f_zero
1699  sums = f_zero
1700  sumok = f_zero
1701  sumokg = f_zero
1702  sumreft = f_zero
1703 
1704  do ni = nr1(nb), nr2(nb)
1705  sp = sqrt( (f_one - rhidssa0(ni,nc)) &
1706  & / (f_one - rhidssa0(ni,nc)*rhidasy0(ni,nc)) )
1707  reft = (f_one - sp) / (f_one + sp)
1708  sumreft = sumreft + reft*eirwaer(nb,ni)
1709 
1710  sumk = sumk + rhidext0(ni,nc)*eirwaer(nb,ni)
1711  sums = sums + rhidsca0(ni,nc)*eirwaer(nb,ni)
1712  sumok = sumok + rhidssa0(ni,nc)*eirwaer(nb,ni) &
1713  & * rhidext0(ni,nc)
1714  sumokg = sumokg + rhidssa0(ni,nc)*eirwaer(nb,ni) &
1715  & * rhidext0(ni,nc)*rhidasy0(ni,nc)
1716  enddo
1717 
1718  refb = sumreft * rirbd
1719 
1720  extrhi(nc,ib) = sumk * rirbd
1721  scarhi(nc,ib) = sums * rirbd
1722  asyrhi(nc,ib) = sumokg / (sumok + 1.0e-10)
1723  ssarhi(nc,ib) = 4.0*refb &
1724  & / ( (f_one+refb)**2 - asyrhi(nc,ib)*(f_one-refb)**2 )
1725  enddo ! end do_nc_block for rh-ind aeros
1726 
1727  do nc = 1, ncm2 ! --- for rh dependent aerosols species
1728  do nh = 1, nrhlev
1729  sumk = f_zero
1730  sums = f_zero
1731  sumok = f_zero
1732  sumokg = f_zero
1733  sumreft = f_zero
1734 
1735  do ni = nr1(nb), nr2(nb)
1736  sp = sqrt( (f_one - rhdpssa0(ni,nh,nc)) &
1737  & / (f_one - rhdpssa0(ni,nh,nc)*rhdpasy0(ni,nh,nc)) )
1738  reft = (f_one - sp) / (f_one + sp)
1739  sumreft = sumreft + reft*eirwaer(nb,ni)
1740 
1741  sumk = sumk + rhdpext0(ni,nh,nc)*eirwaer(nb,ni)
1742  sums = sums + rhdpsca0(ni,nh,nc)*eirwaer(nb,ni)
1743  sumok = sumok + rhdpssa0(ni,nh,nc)*eirwaer(nb,ni) &
1744  & * rhdpext0(ni,nh,nc)
1745  sumokg = sumokg + rhdpssa0(ni,nh,nc)*eirwaer(nb,ni) &
1746  & * rhdpext0(ni,nh,nc)*rhdpasy0(ni,nh,nc)
1747  enddo
1748 
1749  refb = sumreft * rirbd
1750 
1751  extrhd(nh,nc,ib) = sumk * rirbd
1752  scarhd(nh,nc,ib) = sums * rirbd
1753  asyrhd(nh,nc,ib) = sumokg / (sumok + 1.0e-10)
1754  ssarhd(nh,nc,ib) = 4.0*refb &
1755  & / ( (f_one+refb)**2 - asyrhd(nh,nc,ib)*(f_one-refb)**2 )
1756  enddo ! end do_nh_block
1757  enddo ! end do_nc_block for rh-dep aeros
1758 
1759 ! --- for stratospheric background aerosols
1760 
1761  sumk = f_zero
1762  do ni = nr1(nb), nr2(nb)
1763  sumk = sumk + straext0(ni)*eirwaer(nb,ni)
1764  enddo
1765 
1766  extstra(ib) = sumk * rirbd
1767 
1768 ! --- check print
1769 ! if ( nb >= 1 .and. nb < 5) then
1770 ! print *,' in optavg for ir band:',nb
1771 ! print *,' nr1, nr2:',nr1(nb),nr2(nb)
1772 ! print *,' eirwaer:',eirwaer(nb,nr1(nb):nr2(nb))
1773 ! print *,' extrhi:', extrhi(:,ib)
1774 ! do i = 1, NRHLEV
1775 ! print *,' extrhd for rhlev:',i
1776 ! print *,extrhd(i,:,ib)
1777 ! enddo
1778 ! print *,' sumk, rirbd, extstra:',sumk,rirbd,extstra(ib)
1779 ! endif
1780 
1781  enddo ! end do_nb_block for lw
1782  endif ! end if_lalwflg_block
1783 !
1784  return
1785 !................................
1786  end subroutine optavg
1787 !--------------------------------
1788 !
1789 !...................................
1790  end subroutine clim_aerinit
1791 !-----------------------------------
1792 !!@}
1793 
1794 
1803 !-----------------------------------
1804  subroutine aer_update &
1805  & ( iyear, imon, me ) ! --- inputs:
1806 ! --- outputs: ( to module variables )
1807 
1808 ! ================================================================== !
1809 ! !
1810 ! aer_update checks and update time varying climatology aerosol !
1811 ! data sets. !
1812 ! !
1813 ! inputs: size !
1814 ! iyear - 4-digit calender year 1 !
1815 ! imon - month of the year 1 !
1816 ! me - print message control flag 1 !
1817 ! !
1818 ! outputs: ( none ) !
1819 ! !
1820 ! external module variables: (in physparam) !
1821 ! lalwflg - control flag for tropospheric lw aerosol !
1822 ! laswflg - control flag for tropospheric sw aerosol !
1823 ! lavoflg - control flag for stratospheric volcanic aerosol !
1824 ! !
1825 ! usage: call aero_update !
1826 ! !
1827 ! subprograms called: trop_update, volc_update !
1828 ! !
1829 ! ================================================================== !
1830 
1831 ! --- inputs:
1832  integer, intent(in) :: iyear, imon, me
1833 
1834 ! --- output: ( none )
1835 
1836 ! --- locals: ( none )
1837 !
1838 !===> ... begin here
1839 !
1840  if ( imon < 1 .or. imon > 12 ) then
1841  print *,' ***** ERROR in specifying requested month !!! ', &
1842  & 'imon=', imon
1843  print *,' ***** STOPPED in subroutinte aer_update !!!'
1844  stop
1845  endif
1846 
1848  if ( lalwflg .or. laswflg ) then
1849  call trop_update
1850  endif
1851 
1853  if ( lavoflg ) then
1854  call volc_update
1855  endif
1856 
1857 
1858 ! =================
1859  contains
1860 ! =================
1861 
1864 !--------------------------------
1865  subroutine trop_update
1866 !................................
1867 ! --- inputs: (in scope variables, module variables)
1868 ! --- outputs: (module variables)
1869 
1870 ! ================================================================== !
1871 ! !
1872 ! subprogram : trop_update !
1873 ! !
1874 ! updates the monthly global distribution of aerosol profiles in !
1875 ! five degree horizontal resolution. !
1876 ! !
1877 ! ==================== defination of variables =================== !
1878 ! !
1879 ! inputs: (in-scope variables, module constants) !
1880 ! imon - integer, month of the year !
1881 ! me - integer, print message control flag !
1882 ! !
1883 ! outputs: (module variables) !
1884 ! !
1885 ! external module variables: (in physparam) !
1886 ! aeros_file - external aerosol data file name !
1887 ! !
1888 ! internal module variables: !
1889 ! kprfg ( IMXAE*JMXAE) - aeros profile index !
1890 ! idxcg (NXC*IMXAE*JMXAE) - aeros component index !
1891 ! cmixg (NXC*IMXAE*JMXAE) - aeros component mixing ratio !
1892 ! denng ( 2 *IMXAE*JMXAE) - aerosols number density !
1893 ! !
1894 ! NIAERCM - unit number for input data set !
1895 ! !
1896 ! subroutines called: none !
1897 ! !
1898 ! usage: call trop_update !
1899 ! !
1900 ! ================================================================== !
1901 
1902 ! --- inputs: ( none )
1903 ! --- output: ( none )
1904 
1905 ! --- locals:
1906 ! real (kind=kind_io8) :: cmix(NXC), denn, tem
1907  real (kind=kind_phys) :: cmix(nxc), denn, tem
1908  integer :: idxc(nxc), kprf
1909 
1910  integer :: i, id, j, k, m, nc
1911  logical :: file_exist
1912 
1913  character :: cline*80, ctyp*3
1914 !
1915 !===> ... begin here
1916 !
1917 ! --- ... reading climatological aerosols data
1918 
1919  inquire (file=aeros_file, exist=file_exist)
1920 
1921  if ( file_exist ) then
1922  close(niaercm)
1923  open (unit=niaercm,file=aeros_file,status='OLD', &
1924  & action='read',form='FORMATTED')
1925  rewind(niaercm)
1926 
1927  if ( me == 0 ) then
1928  print *,' Opened aerosol data file: ',aeros_file
1929  endif
1930  else
1931  print *,' Requested aerosol data file "',aeros_file, &
1932  & '" not found!'
1933  print *,' *** Stopped in subroutine trop_update !!'
1934  stop
1935  endif ! end if_file_exist_block
1936 
1937 !$omp parallel do private(i,j,m)
1938  do j = 1, jmxae
1939  do i = 1, imxae
1940  do m = 1, nxc
1941  idxcg(m,i,j) = 0
1942  cmixg(m,i,j) = f_zero
1943  enddo
1944  enddo
1945  enddo
1946 
1947 !$omp parallel do private(i,j)
1948  do j = 1, jmxae
1949  do i = 1, imxae
1950  denng(1,i,j) = f_zero
1951  denng(2,i,j) = f_zero
1952  enddo
1953  enddo
1954 
1955 ! --- ... loop over 12 month global distribution
1956 
1957  lab_do_12mon : do m = 1, 12
1958 
1959  read(niaercm,12) cline
1960  12 format(a80/)
1961 
1962  if ( m /= imon ) then
1963 ! if ( me == 0 ) print *,' *** Skipped ',cline
1964 
1965  do j = 1, jmxae
1966  do i = 1, imxae
1967  read(niaercm,*) id
1968  enddo
1969  enddo
1970  else
1971  if ( me == 0 ) print *,' --- Reading ',cline
1972 
1973  do j = 1, jmxae
1974  do i = 1, imxae
1975  read(niaercm,14) (idxc(k),cmix(k),k=1,nxc),kprf,denn,nc, &
1976  & ctyp
1977  14 format(5(i2,e11.4),i2,f8.2,i3,1x,a3)
1978 
1979  kprfg(i,j) = kprf
1980  denng(1,i,j) = denn ! num density of 1st layer
1981  if ( kprf >= 6 ) then
1982  denng(2,i,j) = cmix(nxc) ! num density of 2dn layer
1983  else
1984  denng(2,i,j) = f_zero
1985  endif
1986 
1987  tem = f_one
1988  do k = 1, nxc-1
1989  idxcg(k,i,j) = idxc(k) ! component index
1990  cmixg(k,i,j) = cmix(k) ! component mixing ratio
1991  tem = tem - cmix(k)
1992  enddo
1993  idxcg(nxc,i,j) = idxc(nxc)
1994  cmixg(nxc,i,j) = tem ! to make sure all add to 1.
1995  enddo
1996  enddo
1997 
1998  close (niaercm)
1999  exit lab_do_12mon
2000  endif ! end if_m_block
2001 
2002  enddo lab_do_12mon
2003 
2004 ! -- check print
2005 
2006 ! print *,' IDXCG :'
2007 ! print 16,idxcg
2008 ! 16 format(40i3)
2009 ! print *,' CMIXG :'
2010 ! print 17,cmixg
2011 ! print *,' DENNG :'
2012 ! print 17,denng
2013 ! print *,' KPRFG :'
2014 ! print 17,kprfg
2015 ! 17 format(8e16.9)
2016 !
2017  return
2018 !................................
2019  end subroutine trop_update
2020 !--------------------------------
2021 
2022 
2025 !--------------------------------
2026  subroutine volc_update
2027 !................................
2028 ! --- inputs: (in scope variables, module variables)
2029 ! --- outputs: (module variables)
2030 
2031 ! ================================================================== !
2032 ! !
2033 ! subprogram : volc_update !
2034 ! !
2035 ! searches historical volcanic data sets to find and read in !
2036 ! monthly 45-degree lat-zone band data of optical depth. !
2037 ! !
2038 ! ==================== defination of variables =================== !
2039 ! !
2040 ! inputs: (in-scope variables, module constants) !
2041 ! iyear - integer, 4-digit calender year 1 !
2042 ! imon - integer, month of the year 1 !
2043 ! me - integer, print message control flag 1 !
2044 ! NIAERCM - integer, unit number for input data set 1 !
2045 ! !
2046 ! outputs: (module variables) !
2047 ! ivolae - integer, monthly, 45-deg lat-zone volc odp 12*4*10 !
2048 ! kyrstr - integer, starting year of data in the input file !
2049 ! kyrend - integer, ending year of data in the input file !
2050 ! kyrsav - integer, the year of data in use in the input file !
2051 ! kmonsav - integer, the month of data in use in the input file !
2052 ! !
2053 ! subroutines called: none !
2054 ! !
2055 ! usage: call volc_aerinit !
2056 ! !
2057 ! ================================================================== !
2058 
2059 ! --- inputs: (in-scope variables, module constants)
2060 ! integer :: iyear, imon, me, NIAERCM
2061 
2062 ! --- output: (module variables)
2063 ! integer :: ivolae(:,:,:), kyrstr, kyrend, kyrsav, kmonsav
2064 
2065 ! --- locals:
2066  integer :: i, j, k
2067  logical :: file_exist
2068 
2069  character :: cline*80, volcano_file*32
2070  data volcano_file / 'volcanic_aerosols_1850-1859.txt ' /
2071 !
2072 !===> ... begin here
2073 !
2074  kmonsav = imon
2075 
2076  if ( kyrstr<=iyear .and. iyear<=kyrend ) then ! use previously input data
2077  kyrsav = iyear
2078  return
2079  else ! need to input new data
2080  kyrsav = iyear
2081  kyrstr = iyear - mod(iyear,10)
2082  kyrend = kyrstr + 9
2083 
2084 ! --- check print
2085 ! print *,' kyrstr, kyrend, kyrsav, kmonsav =', &
2086 ! & kyrstr,kyrend,kyrsav,kmonsav
2087 
2088  if ( iyear < minvyr .or. iyear > maxvyr ) then
2089 ! if ( .not. allocated(ivolae) ) then
2090 ! allocate ( ivolae(12,4,10) ) ! for 12-mon,4-lat_zone,10-year
2091 ! endif
2092  ivolae(:,:,:) = 1 ! set as lowest value
2093  if ( me == 0 ) then
2094  print *,' Request volcanic date out of range,', &
2095  & ' optical depth set to lowest value'
2096  endif
2097  else
2098  write(volcano_file(19:27),60) kyrstr,kyrend
2099  60 format(i4.4,'-',i4.4)
2100 
2101  inquire (file=volcano_file, exist=file_exist)
2102  if ( file_exist ) then
2103  close(niaercm)
2104  open (unit=niaercm,file=volcano_file,status='OLD', &
2105  & action='read',form='FORMATTED')
2106 
2107  read(niaercm,62) cline
2108  62 format(a80)
2109 
2110 ! --- check print
2111  if ( me == 0 ) then
2112  print *,' Opened volcanic data file: ',volcano_file
2113  print *, cline
2114  endif
2115 
2116  do k = 1, 10
2117  do j = 1, 4
2118  read(niaercm,64) (ivolae(i,j,k),i=1,12)
2119  64 format(12i5)
2120  enddo
2121  enddo
2122 
2123  close (niaercm)
2124  else
2125  print *,' Requested volcanic data file "', &
2126  & volcano_file,'" not found!'
2127  print *,' *** Stopped in subroutine VOLC_AERINIT !!'
2128  stop
2129  endif ! end if_file_exist_block
2130 
2131  endif ! end if_iyear_block
2132  endif ! end if_kyrstr_block
2133 
2134 ! --- check print
2135  if ( me == 0 ) then
2136  k = mod(kyrsav,10) + 1
2137  print *,' CHECK: Sample Volcanic data used for month, year:', &
2138  & imon, iyear
2139  print *, ivolae(kmonsav,:,k)
2140  endif
2141 !
2142  return
2143 !................................
2144  end subroutine volc_update
2145 !--------------------------------
2146 !
2147 !...................................
2148  end subroutine aer_update
2149 !-----------------------------------
2150 !! @}
2151 
2152 
2179 !-----------------------------------
2180  subroutine setaer &
2181  & ( prsi,prsl,prslk,tvly,rhlay,slmsk,tracer,xlon,xlat, & ! --- inputs
2182  & imax,nlay,nlp1, lsswr,lslwr, &
2183  & aerosw,aerolw & ! --- outputs
2184  &, aerodp &
2185  & )
2186 
2187 ! ================================================================== !
2188 ! !
2189 ! setaer computes aerosols optical properties !
2190 ! !
2191 ! inputs: size !
2192 ! prsi - pressure at interface mb IMAX*NLP1 !
2193 ! prsl - layer mean pressure mb IMAX*NLAY !
2194 ! prslk - exner function = (p/p0)**rocp IMAX*NLAY !
2195 ! tvly - layer virtual temperature k IMAX*NLAY !
2196 ! rhlay - layer mean relative humidity IMAX*NLAY !
2197 ! slmsk - sea/land mask (sea:0,land:1,sea-ice:2) IMAX !
2198 ! tracer - aerosol tracer concentration IMAX*NLAY*NTRAC !
2199 ! xlon - longitude of given points in radiance IMAX !
2200 ! ok for both 0->2pi or -pi->+pi ranges !
2201 ! xlat - latitude of given points in radiance IMAX !
2202 ! default to pi/2 -> -pi/2, otherwise see in-line comment!
2203 ! IMAX - horizontal dimension of arrays 1 !
2204 ! NLAY,NLP1-vertical dimensions of arrays 1 !
2205 ! lsswr,lslwr !
2206 ! - logical flags for sw/lw radiation calls 1 !
2207 ! !
2208 ! outputs: !
2209 ! aerosw - aeros opt properties for sw IMAX*NLAY*NBDSW*NF_AESW!
2210 ! (:,:,:,1): optical depth !
2211 ! (:,:,:,2): single scattering albedo !
2212 ! (:,:,:,3): asymmetry parameter !
2213 ! aerolw - aeros opt properties for lw IMAX*NLAY*NBDLW*NF_AELW!
2214 ! (:,:,:,1): optical depth !
2215 ! (:,:,:,2): single scattering albedo !
2216 ! (:,:,:,3): asymmetry parameter !
2217 ! tau_gocart - 550nm aeros opt depth IMAX*NLAY*MAX_NUM_GRIDCOMP!
2218 !! aerodp - vertically integrated optical depth IMAX*NSPC1 !
2219 ! !
2220 ! external module variable: (in physparam) !
2221 ! iaerflg - aerosol effect control flag (volc,lw,sw, 3-dig) !
2222 ! laswflg - tropospheric aerosol control flag for sw radiation !
2223 ! =f: no sw aeros calc. =t: do sw aeros calc. !
2224 ! lalwflg - tropospheric aerosol control flag for lw radiation !
2225 ! =f: no lw aeros calc. =t: do lw aeros calc. !
2226 ! lavoflg - control flag for stratospheric vocanic aerosols !
2227 ! =t: add volcanic aerosols to the background aerosols !
2228 ! ivflip - control flag for direction of vertical index !
2229 ! =0: index from toa to surface !
2230 ! =1: index from surface to toa !
2231 ! !
2232 ! internal module variable: (set by subroutine aer_init) !
2233 ! ivolae - stratosphere volcanic aerosol optical depth (fac 1.e4) !
2234 ! 12*4*10 !
2235 ! usage: call setaer !
2236 ! !
2237 ! subprograms called: aer_property !
2238 ! !
2239 ! ================================================================== !
2240 
2241 ! --- inputs:
2242  integer, intent(in) :: IMAX, NLAY, NLP1
2243 
2244  real (kind=kind_phys), dimension(:,:), intent(in) :: prsi, prsl, &
2245  & prslk, tvly, rhlay
2246  real (kind=kind_phys), dimension(:), intent(in) :: xlon, xlat, &
2247  & slmsk
2248  real (kind=kind_phys), dimension(:,:,:),intent(in):: tracer
2249 
2250  logical, intent(in) :: lsswr, lslwr
2251 
2252 
2253 ! --- outputs:
2254  real (kind=kind_phys), dimension(:,:,:,:), intent(out) :: &
2255  & aerosw, aerolw
2256 
2257  real (kind=kind_phys), dimension(:,:) , intent(out) :: aerodp
2258 
2259 ! --- locals:
2260  real (kind=kind_phys), parameter :: psrfh = 5.0 ! ref press (mb) for upper bound
2261 
2262  real (kind=kind_phys), dimension(IMAX) :: alon,alat,volcae,rdelp
2263 ! real (kind=kind_phys), dimension(IMAX) :: sumodp
2264  real (kind=kind_phys) :: prsln(nlp1),hz(imax,nlp1),dz(imax,nlay)
2265  real (kind=kind_phys) :: tmp1, tmp2, psrfl
2266 
2267  integer :: kcutl(imax), kcuth(imax)
2268  integer :: i, i1, j, k, m, mb, kh, kl
2269 
2270  logical :: laddsw=.false., laersw=.false.
2271  logical :: laddlw=.false., laerlw=.false.
2272 
2273 ! --- conversion constants
2274  real (kind=kind_phys), parameter :: rdg = 180.0 / con_pi
2275  real (kind=kind_phys), parameter :: rovg = 0.001 * con_rd / con_g
2276 
2277 !===> ... begin here
2278 
2279  do m = 1, nf_aesw
2280  do j = 1, nbdsw
2281  do k = 1, nlay
2282  do i = 1, imax
2283  aerosw(i,k,j,m) = f_zero
2284  enddo
2285  enddo
2286  enddo
2287  enddo
2288 
2289  do m = 1, nf_aelw
2290  do j = 1, nbdlw
2291  do k = 1, nlay
2292  do i = 1, imax
2293  aerolw(i,k,j,m) = f_zero
2294  enddo
2295  enddo
2296  enddo
2297  enddo
2298 
2299 ! sumodp = f_zero
2300  do i = 1, imax
2301  do k = 1, nspc1
2302  aerodp(i,k) = f_zero
2303  enddo
2304  enddo
2305 
2306 
2307  if ( .not. (lsswr .or. lslwr) ) then
2308  return
2309  endif
2310 
2311  if ( iaerflg <= 0 ) then
2312  return
2313  endif
2314 
2315  laersw = lsswr .and. laswflg
2316  laerlw = lslwr .and. lalwflg
2317 
2319 
2320  do i = 1, imax
2321  alon(i) = xlon(i) * rdg
2322  if (alon(i) < f_zero) alon(i) = alon(i) + 360.0
2323  alat(i) = xlat(i) * rdg ! if xlat in pi/2 -> -pi/2 range
2324 ! alat(i) = 90.0 - xlat(i)*rdg ! if xlat in 0 -> pi range
2325  enddo
2326 
2328 
2329  if ( laswflg .or. lalwflg ) then
2330 
2331  lab_do_imax : do i = 1, imax
2332 
2333  lab_if_flip : if (ivflip == 1) then ! input from sfc to toa
2334 
2335  do k = 1, nlay
2336  prsln(k) = log(prsi(i,k))
2337  enddo
2338  prsln(nlp1)= log(prsl(i,nlay))
2339 
2340  do k = nlay, 1, -1
2341  dz(i,k) = rovg * (prsln(k) - prsln(k+1)) * tvly(i,k)
2342  enddo
2343  dz(i,nlay) = 2.0 * dz(i,nlay)
2344 
2345  hz(i,1) = f_zero
2346  do k = 1, nlay
2347  hz(i,k+1) = hz(i,k) + dz(i,k)
2348  enddo
2349 
2350  else lab_if_flip ! input from toa to sfc
2351 
2352  prsln(1) = log(prsl(i,1))
2353  do k = 2, nlp1
2354  prsln(k) = log(prsi(i,k))
2355  enddo
2356 
2357  do k = 1, nlay
2358  dz(i,k) = rovg * (prsln(k+1) - prsln(k)) * tvly(i,k)
2359  enddo
2360  dz(i,1) = 2.0 * dz(i,1)
2361 
2362  hz(i,nlp1) = f_zero
2363  do k = nlay, 1, -1
2364  hz(i,k) = hz(i,k+1) + dz(i,k)
2365  enddo
2366 
2367  endif lab_if_flip
2368 
2369  enddo lab_do_imax
2370 
2371 
2381 
2382 !SARAH
2383 ! if ( iaerflg == 1 ) then ! use opac aerosol climatology
2384  if ( iaermdl == 0 ) then ! use opac aerosol climatology
2385 
2386  call aer_property &
2387 ! --- inputs:
2388  & ( prsi,prsl,prslk,tvly,rhlay,dz,hz,tracer, &
2389  & alon,alat,slmsk, laersw,laerlw, &
2390  & imax,nlay,nlp1, &
2391 ! & IMAX,NLAY,NLP1,NSPC1, &
2392 ! --- outputs:
2393  & aerosw,aerolw,aerodp &
2394  & )
2395 
2396 ! --- check print
2397 ! do m = 1, NBDSW
2398 ! print *,' *** CHECK AEROSOLS PROPERTIES FOR SW BAND =',m, &
2399 ! & ' ***'
2400 ! do k = 1, 10
2401 ! print *,' LEVEL :',k
2402 ! print *,' TAUAER:',aerosw(:,k,m,1)
2403 ! print *,' SSAAER:',aerosw(:,k,m,2)
2404 ! print *,' ASYAER:',aerosw(:,k,m,3)
2405 ! enddo
2406 ! enddo
2407 ! print *,' *** CHECK AEROSOLS OPTICAL DEPTH FOR 550nm REGION'
2408 ! print *, aerodp(:,1)
2409 ! if ( laod_out ) then
2410 ! do m = 1, NSPC1
2411 ! print *,' *** CHECK AEROSOLS OPTICAL DEPTH FOR SPECIES:', &
2412 ! & m
2413 ! print *, aerodp(:,m)
2414 ! sumodp(:) = sumodp(:) + aerodp(:,m)
2415 ! enddo
2416 
2417 !
2418 ! print *,' *** CHECK AEROSOLS OPTICAL DEPTH FOR ALL SPECIES:'
2419 ! print *, sumodp(:)
2420 ! endif
2421 ! do m = 1, NBDLW
2422 ! print *,' *** CHECK AEROSOLS PROPERTIES FOR LW BAND =',m, &
2423 ! & ' ***'
2424 ! do k = 1, 10
2425 ! print *,' LEVEL :',k
2426 ! print *,' TAUAER:',aerolw(:,k,m,1)
2427 ! print *,' SSAAER:',aerolw(:,k,m,2)
2428 ! print *,' ASYAER:',aerolw(:,k,m,3)
2429 ! enddo
2430 ! enddo
2431 ! SARAH
2432 ! elseif ( iaerflg == 2 ) then ! use gocart aerosol scheme
2433  elseif ( iaermdl == 1 ) then ! use gocart aerosol scheme
2434 
2435  call setgocartaer &
2436 
2437 ! --- inputs:
2438  & ( alon,alat,prslk,rhlay,dz,hz,nswlwbd, &
2439  & prsl,tvly,tracer, &
2440  & imax,nlay,nlp1, ivflip, lsswr,lslwr, &
2441 ! --- outputs:
2442  & aerosw,aerolw &
2443  & )
2444 
2445  endif ! end if_iaerflg_block
2446 
2447  endif ! end if_laswflg_or_lalwflg_block
2448 
2456 ! --- ... stratosphere volcanic forcing
2457 
2458  if ( lavoflg ) then
2459 
2460  if ( iaerflg == 100 ) then
2461  laddsw = lsswr
2462  laddlw = lslwr
2463  else
2464  laddsw = lsswr .and. laswflg
2465  laddlw = lslwr .and. lalwflg
2466  endif
2467 
2468  i1 = mod(kyrsav, 10) + 1
2469 
2470 ! --- select data in 4 lat bands, interpolation at the boundaires
2471 
2472  do i = 1, imax
2473  if ( alat(i) > 46.0 ) then
2474  volcae(i) = 1.0e-4 * ivolae(kmonsav,1,i1)
2475  else if ( alat(i) > 44.0 ) then
2476  volcae(i) = 5.0e-5 &
2477  & * (ivolae(kmonsav,1,i1) + ivolae(kmonsav,2,i1))
2478  else if ( alat(i) > 1.0 ) then
2479  volcae(i) = 1.0e-4 * ivolae(kmonsav,2,i1)
2480  else if ( alat(i) > -1.0 ) then
2481  volcae(i) = 5.0e-5 &
2482  & * (ivolae(kmonsav,2,i1) + ivolae(kmonsav,3,i1))
2483  else if ( alat(i) >-44.0 ) then
2484  volcae(i) = 1.0e-4 * ivolae(kmonsav,3,i1)
2485  else if ( alat(i) >-46.0 ) then
2486  volcae(i) = 5.0e-5 &
2487  & * (ivolae(kmonsav,3,i1) + ivolae(kmonsav,4,i1))
2488  else
2489  volcae(i) = 1.0e-4 * ivolae(kmonsav,4,i1)
2490  endif
2491  enddo
2492 
2493  if ( ivflip == 0 ) then ! input data from toa to sfc
2494 
2495 ! --- find lower boundary of stratosphere
2496 
2497  do i = 1, imax
2498 
2499  tmp1 = abs( alat(i) )
2500  if ( tmp1 > 70.0 ) then ! polar, fixed at 25000pa (250mb)
2501  psrfl = 250.0
2502  elseif ( tmp1 < 20.0 ) then ! tropic, fixed at 15000pa (150mb)
2503  psrfl = 150.0
2504  else ! mid-lat, interpolation
2505  psrfl = 110.0 + 2.0*tmp1
2506  endif
2507 
2508  kcuth(i) = nlay - 1
2509  kcutl(i) = 2
2510  rdelp(i) = f_one / prsi(i,2)
2511 
2512  lab_do_kcuth0 : do k = 2, nlay-2
2513  if ( prsi(i,k) >= psrfh ) then
2514  kcuth(i) = k - 1
2515  exit lab_do_kcuth0
2516  endif
2517  enddo lab_do_kcuth0
2518 
2519  lab_do_kcutl0 : do k = 2, nlay-2
2520  if ( prsi(i,k) >= psrfl ) then
2521  kcutl(i) = k - 1
2522  rdelp(i) = f_one / (prsi(i,k) - prsi(i,kcuth(i)))
2523  exit lab_do_kcutl0
2524  endif
2525  enddo lab_do_kcutl0
2526  enddo
2527 
2528 ! --- sw: add volcanic aerosol optical depth to the background value
2529 
2530  if ( laddsw ) then
2531  do m = 1, nbdsw
2532  mb = nswstr + m - 1
2533 
2534  if ( wvnsw1(mb) > 20000 ) then ! range of wvlth < 0.5mu
2535  tmp2 = 0.74
2536  elseif ( wvnsw2(mb) < 20000 ) then ! range of wvlth > 0.5mu
2537  tmp2 = 1.14
2538  else ! range of wvlth in btwn
2539  tmp2 = 0.94
2540  endif
2541  tmp1 = (0.275e-4 * (wvnsw2(mb)+wvnsw1(mb))) ** tmp2
2542 
2543  do i = 1, imax
2544  kh = kcuth(i)
2545  kl = kcutl(i)
2546  do k = kh, kl
2547  tmp2 = tmp1 * ((prsi(i,k+1) - prsi(i,k)) * rdelp(i))
2548  aerosw(i,k,m,1) = aerosw(i,k,m,1) + tmp2*volcae(i)
2549  enddo
2550 
2551 ! --- smoothing profile at boundary if needed
2552 
2553  if ( aerosw(i,kl,m,1) > 10.*aerosw(i,kl+1,m,1) ) then
2554  tmp2 = aerosw(i,kl,m,1) + aerosw(i,kl+1,m,1)
2555  aerosw(i,kl ,m,1) = 0.8 * tmp2
2556  aerosw(i,kl+1,m,1) = 0.2 * tmp2
2557  endif
2558  enddo ! end do_i_block
2559  enddo ! end do_m_block
2560 
2561 ! --- check print
2562 ! do i = 1, IMAX
2563 ! print *,' LEV PRESS TAU FOR PROFILE:',i, &
2564 ! & ' KCUTH, KCUTL =',kcuth(i),kcutl(i)
2565 ! kh = kcuth(i) - 1
2566 ! kl = kcutl(i) + 10
2567 ! do k = kh, kl
2568 ! write(6,71) k, prsl(i,k), aerosw(i,k,1,1)
2569 ! 71 format(i3,2e11.4)
2570 ! enddo
2571 ! enddo
2572 
2573  endif ! end if_laddsw_block
2574 
2575 ! --- lw: add volcanic aerosol optical depth to the background value
2576 
2577  if ( laddlw ) then
2578  if ( nlwbnd == 1 ) then
2579 
2580  tmp1 = (0.55 / 11.0) ** 1.2
2581  do i = 1, imax
2582  kh = kcuth(i)
2583  kl = kcutl(i)
2584  do k = kh, kl
2585  tmp2 = tmp1 * ((prsi(i,k+1) - prsi(i,k)) * rdelp(i)) &
2586  & * volcae(i)
2587  do m = 1, nbdlw
2588  aerolw(i,k,m,1) = aerolw(i,k,m,1) + tmp2
2589  enddo
2590  enddo
2591  enddo ! end do_i_block
2592 
2593  else
2594 
2595  do m = 1, nbdlw
2596  tmp1 = (0.275e-4 * (wvnlw2(m) + wvnlw1(m))) ** 1.2
2597 
2598  do i = 1, imax
2599  kh = kcuth(i)
2600  kl = kcutl(i)
2601  do k = kh, kl
2602  tmp2 = tmp1 * ((prsi(i,k+1)-prsi(i,k)) * rdelp(i))
2603  aerolw(i,k,m,1) = aerolw(i,k,m,1) + tmp2*volcae(i)
2604  enddo
2605  enddo ! end do_i_block
2606  enddo ! end do_m_block
2607 
2608  endif ! end if_NLWBND_block
2609  endif ! end if_laddlw_block
2610 
2611  else ! input data from sfc to toa
2612 
2613 ! --- find lower boundary of stratosphere
2614 
2615  do i = 1, imax
2616 
2617  tmp1 = abs( alat(i) )
2618  if ( tmp1 > 70.0 ) then ! polar, fixed at 25000pa (250mb)
2619  psrfl = 250.0
2620  elseif ( tmp1 < 20.0 ) then ! tropic, fixed at 15000pa (150mb)
2621  psrfl = 150.0
2622  else ! mid-lat, interpolation
2623  psrfl = 110.0 + 2.0*tmp1
2624  endif
2625 
2626  kcuth(i) = 2
2627  kcutl(i) = nlay - 1
2628  rdelp(i) = f_one / prsi(i,nlay-1)
2629 
2630  lab_do_kcuth1 : do k = nlay-1, 2, -1
2631  if ( prsi(i,k) >= psrfh ) then
2632  kcuth(i) = k
2633  exit lab_do_kcuth1
2634  endif
2635  enddo lab_do_kcuth1
2636 
2637  lab_do_kcutl1 : do k = nlay, 2, -1
2638  if ( prsi(i,k) >= psrfl ) then
2639  kcutl(i) = k
2640  rdelp(i) = f_one / (prsi(i,k) - prsi(i,kcuth(i)+1))
2641  exit lab_do_kcutl1
2642  endif
2643  enddo lab_do_kcutl1
2644  enddo
2645 
2646 ! --- sw: add volcanic aerosol optical depth to the background value
2647 
2648  if ( laddsw ) then
2649  do m = 1, nbdsw
2650  mb = nswstr + m - 1
2651 
2652  if ( wvnsw1(mb) > 20000 ) then ! range of wvlth < 0.5mu
2653  tmp2 = 0.74
2654  elseif ( wvnsw2(mb) < 20000 ) then ! range of wvlth > 0.5mu
2655  tmp2 = 1.14
2656  else ! range of wvlth in btwn
2657  tmp2 = 0.94
2658  endif
2659  tmp1 = (0.275e-4 * (wvnsw2(mb)+wvnsw1(mb))) ** tmp2
2660 
2661  do i = 1, imax
2662  kh = kcuth(i)
2663  kl = kcutl(i)
2664  do k = kl, kh
2665  tmp2 = tmp1 * ((prsi(i,k) - prsi(i,k+1)) * rdelp(i))
2666  aerosw(i,k,m,1) = aerosw(i,k,m,1) + tmp2*volcae(i)
2667  enddo
2668 
2669 ! --- smoothing profile at boundary if needed
2670 
2671  if ( aerosw(i,kl,m,1) > 10.*aerosw(i,kl-1,m,1) ) then
2672  tmp2 = aerosw(i,kl,m,1) + aerosw(i,kl-1,m,1)
2673  aerosw(i,kl ,m,1) = 0.8 * tmp2
2674  aerosw(i,kl-1,m,1) = 0.2 * tmp2
2675  endif
2676  enddo ! end do_i_block
2677  enddo ! end do_m_block
2678 
2679 ! --- check print
2680 ! do i = 1, IMAX
2681 ! print *,' LEV PRESS TAU FOR PROFILE:',i, &
2682 ! & ' KCUTH, KCUTL =',kcuth(i),kcutl(i)
2683 ! kh = kcuth(i) + 1
2684 ! kl = kcutl(i) - 10
2685 ! do k = kh, kl, -1
2686 ! write(6,71) NLP1-k,prsl(i,k),aerosw(i,k,1,1)
2687 ! enddo
2688 ! enddo
2689 
2690  endif ! end if_laddsw_block
2691 
2692 ! --- lw: add volcanic aerosol optical depth to the background value
2693 
2694  if ( laddlw ) then
2695  if ( nlwbnd == 1 ) then
2696 
2697  tmp1 = (0.55 / 11.0) ** 1.2
2698  do i = 1, imax
2699  kh = kcuth(i)
2700  kl = kcutl(i)
2701  do k = kl, kh
2702  tmp2 = tmp1 * ((prsi(i,k) - prsi(i,k+1)) * rdelp(i)) &
2703  & * volcae(i)
2704  do m = 1, nbdlw
2705  aerolw(i,k,m,1) = aerolw(i,k,m,1) + tmp2
2706  enddo
2707  enddo
2708  enddo ! end do_i_block
2709 
2710  else
2711 
2712  do m = 1, nbdlw
2713  tmp1 = (0.275e-4 * (wvnlw2(m) + wvnlw1(m))) ** 1.2
2714 
2715  do i = 1, imax
2716  kh = kcuth(i)
2717  kl = kcutl(i)
2718  do k = kl, kh
2719  tmp2 = tmp1 * ((prsi(i,k)-prsi(i,k+1)) * rdelp(i))
2720  aerolw(i,k,m,1) = aerolw(i,k,m,1) + tmp2*volcae(i)
2721  enddo
2722  enddo ! end do_i_block
2723  enddo ! end do_m_block
2724 
2725  endif ! end if_NLWBND_block
2726  endif ! end if_laddlw_block
2727 
2728  endif ! end if_ivflip_block
2729 
2730  endif ! end if_lavoflg_block
2731 !
2732  return
2733 !...................................
2734  end subroutine setaer
2735 !-----------------------------------
2737 
2738 
2767 !-----------------------------------
2768  subroutine aer_property &
2769  & ( prsi,prsl,prslk,tvly,rhlay,dz,hz,tracer, & ! --- inputs:
2770  & alon,alat,slmsk, laersw,laerlw, &
2771  & imax,nlay,nlp1, &
2772  & aerosw,aerolw,aerodp & ! --- outputs:
2773  & )
2774 
2775 ! ================================================================== !
2776 ! !
2777 ! aer_property maps the 5 degree global climatological aerosol data !
2778 ! set onto model grids, and compute aerosol optical properties for sw !
2779 ! and lw radiations. !
2780 ! !
2781 ! inputs: !
2782 ! prsi - pressure at interface mb IMAX*NLP1 !
2783 ! prsl - layer mean pressure (not used) IMAX*NLAY !
2784 ! prslk - exner function=(p/p0)**rocp (not used) IMAX*NLAY !
2785 ! tvly - layer virtual temperature (not used) IMAX*NLAY !
2786 ! rhlay - layer mean relative humidity IMAX*NLAY !
2787 ! dz - layer thickness m IMAX*NLAY !
2788 ! hz - level high m IMAX*NLP1 !
2789 ! tracer - aer tracer concentrations (not used) IMAX*NLAY*NTRAC!
2790 ! alon, alat IMAX !
2791 ! - longitude and latitude of given points in degree !
2792 ! slmsk - sea/land mask (sea:0,land:1,sea-ice:2) IMAX !
2793 ! laersw,laerlw 1 !
2794 ! - logical flag for sw/lw aerosol calculations !
2795 ! IMAX - horizontal dimension of arrays 1 !
2796 ! NLAY,NLP1-vertical dimensions of arrays 1 !
2797 !! NSPC - num of species for optional aod output fields 1 !
2798 ! !
2799 ! outputs: !
2800 ! aerosw - aeros opt properties for sw IMAX*NLAY*NBDSW*NF_AESW!
2801 ! (:,:,:,1): optical depth !
2802 ! (:,:,:,2): single scattering albedo !
2803 ! (:,:,:,3): asymmetry parameter !
2804 ! aerolw - aeros opt properties for lw IMAX*NLAY*NBDLW*NF_AELW!
2805 ! (:,:,:,1): optical depth !
2806 ! (:,:,:,2): single scattering albedo !
2807 ! (:,:,:,3): asymmetry parameter !
2808 !! aerodp - vertically integrated aer-opt-depth IMAX*NSPC+1 !
2809 ! !
2810 ! module parameters and constants: !
2811 ! NSWBND - total number of actual sw spectral bands computed !
2812 ! NLWBND - total number of actual lw spectral bands computed !
2813 ! NSWLWBD - total number of sw+lw bands computed !
2814 ! !
2815 ! external module variables: (in physparam) !
2816 ! ivflip - control flag for direction of vertical index !
2817 ! =0: index from toa to surface !
2818 ! =1: index from surface to toa !
2819 ! !
2820 ! module variable: (set by subroutine aer_init) !
2821 ! kprfg - aerosols profile index IMXAE*JMXAE !
2822 ! 1:ant 2:arc 3:cnt 4:mar 5:des 6:marme 7:cntme !
2823 ! idxcg - aerosols component index NXC*IMXAE*JMXAE !
2824 ! 1:inso 2:soot 3:minm 4:miam 5:micm !
2825 ! 6:mitr 7:waso 8:ssam 9:sscm 10:suso !
2826 ! cmixg - aerosols component mixing ratio NXC*IMXAE*JMXAE !
2827 ! denng - aerosols number density 2 *IMXAE*JMXAE !
2828 ! 1:for domain-1 2:domain-2 (prof marme/cntme only) !
2829 ! !
2830 ! usage: call aer_property !
2831 ! !
2832 ! subprograms called: radclimaer !
2833 ! !
2834 ! ================================================================== !
2835 
2836 ! --- inputs:
2837  integer, intent(in) :: IMAX, NLAY, NLP1
2838 ! integer, intent(in) :: IMAX, NLAY, NLP1, NSPC
2839  logical, intent(in) :: laersw, laerlw
2840 
2841  real (kind=kind_phys), dimension(:,:), intent(in) :: prsi, prsl, &
2842  & prslk, tvly, rhlay, dz, hz
2843  real (kind=kind_phys), dimension(:), intent(in) :: alon, alat, &
2844  & slmsk
2845  real (kind=kind_phys), dimension(:,:,:),intent(in):: tracer
2846 
2847 ! --- outputs:
2848  real (kind=kind_phys), dimension(:,:,:,:), intent(out) :: &
2849  & aerosw, aerolw
2850  real (kind=kind_phys), dimension(:,:) , intent(out) :: aerodp
2851 
2852 ! --- locals:
2853  real (kind=kind_phys), dimension(NCM) :: cmix
2854  real (kind=kind_phys), dimension( 2) :: denn
2855  real (kind=kind_phys), dimension(NSPC) :: spcodp
2856 
2857  real (kind=kind_phys), dimension(NLAY) :: delz, rh1, dz1
2858  integer, dimension(NLAY) :: idmaer
2859 
2860  real (kind=kind_phys), dimension(NLAY,NSWLWBD):: tauae,ssaae,asyae
2861 !test real (kind=kind_phys), dimension(IMAX,NLAY) :: aersav
2862 
2863  real (kind=kind_phys) :: tmp1, tmp2, rps, dtmp, h1
2864  real (kind=kind_phys) :: wi, wj, w11, w12, w21, w22
2865 
2866  integer :: i, ii, i1, i2, i3, j1, j2, j3, k, m, m1, &
2867  & kp, kpa, kpi, kpj
2868 
2869 ! --- conversion constants
2870  real (kind=kind_phys), parameter :: dltg = 360.0 / float(imxae)
2871  real (kind=kind_phys), parameter :: hdlt = 0.5 * dltg
2872  real (kind=kind_phys), parameter :: rdlt = 1.0 / dltg
2873 
2874 !
2875 !===> ... begin here
2876 !
2880 
2881  i1 = 1
2882  i2 = 2
2883  j1 = 1
2884  j2 = 2
2885 
2886  lab_do_imax : do i = 1, imax
2887 
2888 ! --- map grid in longitude direction, lon from 0 to 355 deg resolution
2889 
2890 ! print *,' Seeking lon index for point i =',i
2891  i3 = i1
2892  lab_do_imxae : do while ( i3 <= imxae )
2893  tmp1 = dltg * (i3 - 1)
2894  dtmp = alon(i) - tmp1
2895 ! print *,' alon, i3, tlon, dlon =',alon(i),i3,tmp1,dtmp
2896 
2897  if ( dtmp > dltg ) then
2898  i3 = i3 + 1
2899  if ( i3 > imxae ) then
2900  print *,' ERROR! In setclimaer alon>360. ipt =',i, &
2901  & ', dltg,alon,tlon,dlon =',dltg,alon(i),tmp1,dtmp
2902  stop
2903  endif
2904  elseif ( dtmp >= f_zero ) then
2905  i1 = i3
2906  i2 = mod(i3,imxae) + 1
2907  wi = dtmp * rdlt
2908  if ( dtmp <= hdlt ) then
2909  kpi = i3
2910  else
2911  kpi = i2
2912  endif
2913 ! print *,' found i1, i2, wi =',i1,i2,wi
2914  exit lab_do_imxae
2915  else
2916  i3 = i3 - 1
2917  if ( i3 < 1 ) then
2918  print *,' ERROR! In setclimaer alon< 0. ipt =',i, &
2919  & ', dltg,alon,tlon,dlon =',dltg,alon(i),tmp1,dtmp
2920  stop
2921  endif
2922  endif
2923  enddo lab_do_imxae
2924 
2925 ! --- map grid in latitude direction, lat from 90n to 90s in 5 deg resolution
2926 
2927 ! print *,' Seeking lat index for point i =',i
2928  j3 = j1
2929  lab_do_jmxae : do while ( j3 <= jmxae )
2930  tmp2 = 90.0 - dltg * (j3 - 1)
2931  dtmp = tmp2 - alat(i)
2932 ! print *,' alat, j3, tlat, dlat =',alat(i),j3,tmp2,dtmp
2933 
2934  if ( dtmp > dltg ) then
2935  j3 = j3 + 1
2936  if ( j3 >= jmxae ) then
2937  print *,' ERROR! In setclimaer alat<-90. ipt =',i, &
2938  & ', dltg,alat,tlat,dlat =',dltg,alat(i),tmp2,dtmp
2939  stop
2940  endif
2941  elseif ( dtmp >= f_zero ) then
2942  j1 = j3
2943  j2 = j3 + 1
2944  wj = dtmp * rdlt
2945  if ( dtmp <= hdlt ) then
2946  kpj = j3
2947  else
2948  kpj = j2
2949  endif
2950 ! print *,' found j1, j2, wj =',j1,j2,wj
2951  exit lab_do_jmxae
2952  else
2953  j3 = j3 - 1
2954  if ( j3 < 1 ) then
2955  print *,' ERROR! In setclimaer alat>90. ipt =',i, &
2956  & ', dltg,alat,tlat,dlat =',dltg,alat(i),tmp2,dtmp
2957  stop
2958  endif
2959  endif
2960  enddo lab_do_jmxae
2961 
2964 
2965  kp = kprfg(kpi,kpj) ! nearest typical aeros profile as default
2966  kpa = max( kprfg(i1,j1),kprfg(i1,j2),kprfg(i2,j1),kprfg(i2,j2) )
2967  h1 = haer(1,kp)
2968  denn(2) = f_zero
2969  ii = 1
2970 
2971  if ( kp /= kpa ) then
2972  if ( kpa == 6 ) then ! if ocean prof with mineral aeros overlay
2973  ii = 2 ! need 2 types of densities
2974  if ( slmsk(i) > f_zero ) then ! but actually a land/sea-ice point
2975  kp = 7 ! reset prof index to land
2976  h1 = 0.5*(haer(1,6) + haer(1,7)) ! use a transition scale hight
2977  else
2978  kp = kpa
2979  h1 = haer(1,6)
2980  endif
2981  elseif ( kpa == 7 ) then ! if land prof with mineral aeros overlay
2982  ii = 2 ! need 2 types of densities
2983  if ( slmsk(i) <= f_zero ) then ! but actually an ocean point
2984  kp = 6 ! reset prof index to ocean
2985  h1 = 0.5*(haer(1,6) + haer(1,7)) ! use a transition scale hight
2986  else
2987  kp = kpa
2988  h1 = haer(1,7)
2989  endif
2990  else ! lower atmos without mineral aeros overlay
2991 ! h1 = 0.5*(haer(1,kp) + haer(1,kpa)) ! use a transition scale hight
2992  h1 = haer(1,kpa)
2993  kp = kpa
2994  endif
2995  endif
2996 
2998 
2999  w11 = (f_one-wi) * (f_one-wj)
3000  w12 = (f_one-wi) * wj
3001  w21 = wi * (f_one-wj)
3002  w22 = wi * wj
3003 
3004 ! --- check print
3005 ! print *,' Grid pt', i,', alon, alat =',alon(i),alat(i), &
3006 ! & ', tlon, tlat =',tmp1,tmp2
3007 ! print *,' lon grid index i1, i2 =',i1,i2,', weight wi =',wi
3008 ! print *,' lat grid index j1, j2 =',j1,j2,', weight wj =',wj
3009 ! print *,' bi-linear weights w11,w21,w12,w22 =',w11,w21,w12,w22
3010 ! print *,' kp,kpa,slmsk,h1 =',kp,m1,slmsk(i),h1
3011 
3014 
3015  do m = 1, ii ! ii=1 for domain 1; =2 for domain 2.
3016  denn(m) = w11*denng(m,i1,j1) + w12*denng(m,i1,j2) &
3017  & + w21*denng(m,i2,j1) + w22*denng(m,i2,j2)
3018  enddo ! end_do_m_loop
3019 
3021 
3022  cmix(:) = f_zero
3023  do m = 1, nxc
3024  ii = idxcg(m,i1,j1)
3025  if ( ii > 0 ) then
3026  cmix(ii) = cmix(ii) + w11*cmixg(m,i1,j1)
3027  endif
3028  ii = idxcg(m,i1,j2)
3029  if ( ii > 0 ) then
3030  cmix(ii) = cmix(ii) + w12*cmixg(m,i1,j2)
3031  endif
3032  ii = idxcg(m,i2,j1)
3033  if ( ii > 0 ) then
3034  cmix(ii) = cmix(ii) + w21*cmixg(m,i2,j1)
3035  endif
3036  ii = idxcg(m,i2,j2)
3037  if ( ii > 0 ) then
3038  cmix(ii) = cmix(ii) + w22*cmixg(m,i2,j2)
3039  endif
3040  enddo ! end_do_m_loop
3041 
3042 ! --- check print
3043 ! print *,' denn =',denn(:)
3044 ! print *,' cmix =',cmix(:)
3045 
3048 
3049  do k = 1, nlay
3050  rh1(k) = rhlay(i,k)
3051  dz1(k) = dz (i,k)
3052  enddo
3053 
3054  lab_if_flip : if (ivflip == 1) then ! input from sfc to toa
3055 
3056  if ( prsi(i,1) > 100.0 ) then
3057  rps = f_one / prsi(i,1)
3058  else
3059  print *,' !!! (1) Error in subr radiation_aerosols:', &
3060  & ' unrealistic surface pressure =', i,prsi(i,1)
3061  stop
3062  endif
3063 
3064  ii = 1
3065  do k = 1, nlay
3066  if (prsi(i,k+1)*rps < sigref(ii,kp)) then
3067  ii = ii + 1
3068  if (ii == 2 .and. prsref(2,kp) == prsref(3,kp)) then
3069  ii = 3
3070  endif
3071  endif
3072  idmaer(k) = ii
3073 
3074  if ( ii > 1 ) then
3075  tmp1 = haer(ii,kp)
3076  else
3077  tmp1 = h1
3078  endif
3079 
3080  if (tmp1 > f_zero) then
3081  tmp2 = f_one / tmp1
3082  delz(k) = tmp1 * (exp(-hz(i,k)*tmp2)-exp(-hz(i,k+1)*tmp2))
3083  else
3084  delz(k) = dz1(k)
3085  endif
3086  enddo
3087 
3088  else lab_if_flip ! input from toa to sfc
3089 
3090  if ( prsi(i,nlp1) > 100.0 ) then
3091  rps = 1.0 / prsi(i,nlp1)
3092  else
3093  print *,' !!! (2) Error in subr radiation_aerosols:', &
3094  & ' unrealistic surface pressure =', i,prsi(i,nlp1)
3095  endif
3096 
3097  ii = 1
3098  do k = nlay, 1, -1
3099  if (prsi(i,k)*rps < sigref(ii,kp)) then
3100  ii = ii + 1
3101  if (ii == 2 .and. prsref(2,kp) == prsref(3,kp)) then
3102  ii = 3
3103  endif
3104  endif
3105  idmaer(k) = ii
3106 
3107  if ( ii > 1 ) then
3108  tmp1 = haer(ii,kp)
3109  else
3110  tmp1 = h1
3111  endif
3112 
3113  if (tmp1 > f_zero) then
3114  tmp2 = f_one / tmp1
3115  delz(k) = tmp1 * (exp(-hz(i,k+1)*tmp2)-exp(-hz(i,k)*tmp2))
3116  else
3117  delz(k) = dz1(k)
3118  endif
3119  enddo
3120 
3121  endif lab_if_flip
3122 
3123 ! --- check print
3124 
3125 ! print *,' in setclimaer, profile:',i
3126 ! print *,' rh :',rh1
3127 ! print *,' dz :',dz1
3128 ! print *,' delz :',delz
3129 ! print *,' idmaer:',idmaer
3130 
3133 
3134  call radclimaer
3135 ! --- inputs: (in-scope variables)
3136 ! --- outputs: (in-scope variables)
3137 
3138  if ( laersw ) then
3139 
3140  do m = 1, nbdsw
3141  do k = 1, nlay
3142  aerosw(i,k,m,1) = tauae(k,m)
3143  aerosw(i,k,m,2) = ssaae(k,m)
3144  aerosw(i,k,m,3) = asyae(k,m)
3145  enddo
3146  enddo
3147 
3148 ! --- total aod (optional)
3149  do k = 1, nlay
3150  aerodp(i,1) = aerodp(i,1) + tauae(k,nv_aod)
3151  enddo
3152 
3153 ! --- for diagnostic output (optional)
3154 ! if ( lspcaod ) then
3155  do m = 1, nspc
3156  aerodp(i,m+1) = spcodp(m)
3157  enddo
3158 ! endif
3159 
3160  endif ! end if_larsw_block
3161 
3162  if ( laerlw ) then
3163 
3164  if ( nlwbnd == 1 ) then
3165  m1 = nswbnd + 1
3166  do m = 1, nbdlw
3167  do k = 1, nlay
3168  aerolw(i,k,m,1) = tauae(k,m1)
3169  aerolw(i,k,m,2) = ssaae(k,m1)
3170  aerolw(i,k,m,3) = asyae(k,m1)
3171  enddo
3172  enddo
3173  else
3174  do m = 1, nbdlw
3175  m1 = nswbnd + m
3176  do k = 1, nlay
3177  aerolw(i,k,m,1) = tauae(k,m1)
3178  aerolw(i,k,m,2) = ssaae(k,m1)
3179  aerolw(i,k,m,3) = asyae(k,m1)
3180  enddo
3181  enddo
3182  endif
3183 
3184  endif ! end if_laerlw_block
3185 
3186  enddo lab_do_imax
3187 
3188 ! =================
3189  contains
3190 ! =================
3191 
3196 !--------------------------------
3197  subroutine radclimaer
3198 !................................
3199 
3200 ! --- inputs: (in scope variables)
3201 ! --- outputs: (in scope variables)
3202 
3203 ! ================================================================== !
3204 ! !
3205 ! compute aerosols optical properties in NSWLWBD bands. there are !
3206 ! seven different vertical profile structures. in the troposphere, !
3207 ! aerosol distribution at each grid point is composed from up to !
3208 ! six components out of a total of ten different substances. !
3209 ! !
3210 ! ref: wmo report wcp-112 (1986) !
3211 ! !
3212 ! input variables: !
3213 ! cmix - mixing ratioes of aerosol components - NCM !
3214 ! denn - aerosol number densities - 2 !
3215 ! rh1 - relative humidity - NLAY !
3216 ! delz - effective layer thickness km NLAY !
3217 ! idmaer - aerosol domain index - NLAY !
3218 ! NXC - number of different aerosol components- 1 !
3219 ! NLAY - vertical dimensions - 1 !
3220 ! !
3221 ! output variables: !
3222 ! tauae - optical depth - NLAY*NSWLWBD!
3223 ! ssaae - single scattering albedo - NLAY*NSWLWBD!
3224 ! asyae - asymmetry parameter - NLAY*NSWLWBD!
3225 !! aerodp - vertically integrated aer-opt-depth - IMAX*NSPC+1 !
3226 ! !
3227 ! ================================================================== !
3228 !
3229  real (kind=kind_phys) :: crt1, crt2
3230  parameter(crt1=30.0, crt2=0.03333)
3231 
3232 ! --- inputs:
3233 ! --- outputs:
3234 
3235 ! --- locals:
3236  real (kind=kind_phys) :: cm, hd, hdi, sig0u, sig0l, ratio, tt0, &
3237  & ex00, sc00, ss00, as00, ex01, sc01, ss01, as01, tt1, &
3238  & ex02, sc02, ss02, as02, ex03, sc03, ss03, as03, tt2, &
3239  & ext1, sca1, ssa1, asy1, drh0, drh1, rdrh
3240 
3241  integer :: ih1, ih2, kk, idom, icmp, ib, ii, ic, ic1
3242  integer :: idx
3243 
3244 !===> ... begin here
3245 
3246  spcodp = f_zero
3247 
3248 !===> ... loop over vertical layers from top to surface
3249 
3250  lab_do_layer : do kk = 1, nlay
3251 
3252 ! --- linear interp coeffs for rh-dep species
3253 
3254  ih2 = 1
3255  do while ( rh1(kk) > rhlev(ih2) )
3256  ih2 = ih2 + 1
3257  if ( ih2 > nrhlev ) exit
3258  enddo
3259  ih1 = max( 1, ih2-1 )
3260  ih2 = min( nrhlev, ih2 )
3261 
3262  drh0 = rhlev(ih2) - rhlev(ih1)
3263  drh1 = rh1(kk) - rhlev(ih1)
3264  if ( ih1 == ih2 ) then
3265  rdrh = f_zero
3266  else
3267  rdrh = drh1 / drh0
3268  endif
3269 
3270 ! --- assign optical properties in each domain
3271 
3272  idom = idmaer(kk)
3273 
3274  lab_if_idom : if (idom == 5) then
3275 ! --- 5th domain - upper stratosphere assume no aerosol
3276 
3277  do ib = 1, nswlwbd
3278  tauae(kk,ib) = f_zero
3279  if ( ib <= nswbnd ) then
3280  ssaae(kk,ib) = 0.99
3281  asyae(kk,ib) = 0.696
3282  else
3283  ssaae(kk,ib) = 0.5
3284  asyae(kk,ib) = 0.3
3285  endif
3286  enddo
3287 
3288  elseif (idom == 4) then lab_if_idom
3289 ! --- 4th domain - stratospheric layers
3290 
3291  do ib = 1, nswlwbd
3292  tauae(kk,ib) = extstra(ib) * delz(kk)
3293  if ( ib <= nswbnd ) then
3294  ssaae(kk,ib) = 0.99
3295  asyae(kk,ib) = 0.696
3296  else
3297  ssaae(kk,ib) = 0.5
3298  asyae(kk,ib) = 0.3
3299  endif
3300  enddo
3301 
3302 ! --- compute aod from individual species' contribution (optional)
3303  idx = idxspc(10) ! for sulfate
3304  spcodp(idx) = spcodp(idx) + tauae(kk,nv_aod)
3305 
3306  elseif (idom == 3) then lab_if_idom
3307 ! --- 3rd domain - free tropospheric layers
3308 ! 1:inso 0.17e-3; 2:soot 0.4; 7:waso 0.59983; n:730
3309 
3310  do ib = 1, nswlwbd
3311  ex01 = extrhi(1,ib)
3312  sc01 = scarhi(1,ib)
3313  ss01 = ssarhi(1,ib)
3314  as01 = asyrhi(1,ib)
3315 
3316  ex02 = extrhi(2,ib)
3317  sc02 = scarhi(2,ib)
3318  ss02 = ssarhi(2,ib)
3319  as02 = asyrhi(2,ib)
3320 
3321  ex03 = extrhd(ih1,1,ib) &
3322  & + rdrh * (extrhd(ih2,1,ib) - extrhd(ih1,1,ib))
3323  sc03 = scarhd(ih1,1,ib) &
3324  & + rdrh * (scarhd(ih2,1,ib) - scarhd(ih1,1,ib))
3325  ss03 = ssarhd(ih1,1,ib) &
3326  & + rdrh * (ssarhd(ih2,1,ib) - ssarhd(ih1,1,ib))
3327  as03 = asyrhd(ih1,1,ib) &
3328  & + rdrh * (asyrhd(ih2,1,ib) - asyrhd(ih1,1,ib))
3329 
3330  ext1 = 0.17e-3*ex01 + 0.4*ex02 + 0.59983*ex03
3331  sca1 = 0.17e-3*sc01 + 0.4*sc02 + 0.59983*sc03
3332  ssa1 = 0.17e-3*ss01*ex01 + 0.4*ss02*ex02 + 0.59983*ss03*ex03
3333  asy1 = 0.17e-3*as01*sc01 + 0.4*as02*sc02 + 0.59983*as03*sc03
3334 
3335  tauae(kk,ib) = ext1 * 730.0 * delz(kk)
3336  ssaae(kk,ib) = min(f_one, ssa1/ext1)
3337  asyae(kk,ib) = min(f_one, asy1/sca1)
3338 
3339 ! --- compute aod from individual species' contribution (optional)
3340  if ( ib==nv_aod ) then
3341  spcodp(1) = spcodp(1) + 0.17e-3*ex01*730.0*delz(kk) ! dust (inso) #1
3342  spcodp(2) = spcodp(2) + 0.4 *ex02*730.0*delz(kk) ! black carbon #2
3343  spcodp(3) = spcodp(3) + 0.59983*ex03*730.0*delz(kk) ! water soluble #7
3344  endif
3345 
3346  enddo
3347 
3348  elseif (idom == 1) then lab_if_idom
3349 ! --- 1st domain - mixing layer
3350 
3351  lab_do_ib : do ib = 1, nswlwbd
3352  ext1 = f_zero
3353  sca1 = f_zero
3354  ssa1 = f_zero
3355  asy1 = f_zero
3356 
3357  lab_do_icmp : do icmp = 1, ncm
3358  ic = icmp
3359  idx = idxspc(icmp)
3360 
3361  cm = cmix(icmp)
3362  lab_if_cm : if ( cm > f_zero ) then
3363 
3364  lab_if_ic : if ( ic <= ncm1 ) then ! component withour rh dep
3365  tt0 = cm * extrhi(ic,ib)
3366  ext1 = ext1 + tt0
3367  sca1 = sca1 + cm * scarhi(ic,ib)
3368  ssa1 = ssa1 + cm * ssarhi(ic,ib) * extrhi(ic,ib)
3369  asy1 = asy1 + cm * asyrhi(ic,ib) * scarhi(ic,ib)
3370  else lab_if_ic ! component with rh dep
3371  ic1 = ic - ncm1
3372 
3373  ex00 = extrhd(ih1,ic1,ib) &
3374  & + rdrh * (extrhd(ih2,ic1,ib) - extrhd(ih1,ic1,ib))
3375  sc00 = scarhd(ih1,ic1,ib) &
3376  & + rdrh * (scarhd(ih2,ic1,ib) - scarhd(ih1,ic1,ib))
3377  ss00 = ssarhd(ih1,ic1,ib) &
3378  & + rdrh * (ssarhd(ih2,ic1,ib) - ssarhd(ih1,ic1,ib))
3379  as00 = asyrhd(ih1,ic1,ib) &
3380  & + rdrh * (asyrhd(ih2,ic1,ib) - asyrhd(ih1,ic1,ib))
3381 
3382  tt0 = cm * ex00
3383  ext1 = ext1 + tt0
3384  sca1 = sca1 + cm * sc00
3385  ssa1 = ssa1 + cm * ss00 * ex00
3386  asy1 = asy1 + cm * as00 * sc00
3387  endif lab_if_ic
3388 
3389 ! --- compute aod from individual species' contribution (optional)
3390  if ( ib==nv_aod ) then
3391  spcodp(idx) = spcodp(idx) + tt0*denn(1)*delz(kk) ! idx for dif species
3392  endif
3393 
3394  endif lab_if_cm
3395  enddo lab_do_icmp
3396 
3397  tauae(kk,ib) = ext1 * denn(1) * delz(kk)
3398  ssaae(kk,ib) = min(f_one, ssa1/ext1)
3399  asyae(kk,ib) = min(f_one, asy1/sca1)
3400  enddo lab_do_ib
3401 
3402  elseif (idom == 2) then lab_if_idom
3403 ! --- 2nd domain - mineral transport layers
3404 
3405  do ib = 1, nswlwbd
3406  tauae(kk,ib) = extrhi(6,ib) * denn(2) * delz(kk)
3407  ssaae(kk,ib) = ssarhi(6,ib)
3408  asyae(kk,ib) = asyrhi(6,ib)
3409  enddo
3410 
3411 ! --- compute aod from individual species' contribution (optional)
3412  spcodp(1) = spcodp(1) + tauae(kk,nv_aod) ! dust
3413 
3414  else lab_if_idom
3415 ! --- domain index out off range, assume no aerosol
3416 
3417  do ib = 1, nswlwbd
3418  tauae(kk,ib) = f_zero
3419  ssaae(kk,ib) = f_one
3420  asyae(kk,ib) = f_zero
3421  enddo
3422 
3423 ! write(6,19) kk,idom
3424 ! 19 format(/' *** ERROR in sub AEROS: domain index out' &
3425 ! &, ' of range! K, IDOM =',3i5,' ***')
3426 ! stop 19
3427 
3428  endif lab_if_idom
3429 
3430  enddo lab_do_layer
3431 
3432 !
3433 !===> ... smooth profile at domain boundaries
3434 !
3435  if ( ivflip == 0 ) then ! input from toa to sfc
3436 
3437  do ib = 1, nswlwbd
3438  do kk = 2, nlay
3439  if ( tauae(kk,ib) > f_zero ) then
3440  ratio = tauae(kk-1,ib) / tauae(kk,ib)
3441  else
3442  ratio = f_one
3443  endif
3444 
3445  tt0 = tauae(kk,ib) + tauae(kk-1,ib)
3446  tt1 = 0.2 * tt0
3447  tt2 = tt0 - tt1
3448 
3449  if ( ratio > crt1 ) then
3450  tauae(kk,ib) = tt1
3451  tauae(kk-1,ib) = tt2
3452  endif
3453 
3454  if ( ratio < crt2 ) then
3455  tauae(kk,ib) = tt2
3456  tauae(kk-1,ib) = tt1
3457  endif
3458  enddo ! do_kk_loop
3459  enddo ! do_ib_loop
3460 
3461  else ! input from sfc to toa
3462 
3463  do ib = 1, nswlwbd
3464  do kk = nlay-1, 1, -1
3465  if ( tauae(kk,ib) > f_zero ) then
3466  ratio = tauae(kk+1,ib) / tauae(kk,ib)
3467  else
3468  ratio = f_one
3469  endif
3470 
3471  tt0 = tauae(kk,ib) + tauae(kk+1,ib)
3472  tt1 = 0.2 * tt0
3473  tt2 = tt0 - tt1
3474 
3475  if ( ratio > crt1 ) then
3476  tauae(kk,ib) = tt1
3477  tauae(kk+1,ib) = tt2
3478  endif
3479 
3480  if ( ratio < crt2 ) then
3481  tauae(kk,ib) = tt2
3482  tauae(kk+1,ib) = tt1
3483  endif
3484  enddo ! do_kk_loop
3485  enddo ! do_ib_loop
3486 
3487  endif
3488 
3489 !
3490  return
3491 !................................
3492  end subroutine radclimaer
3493 !--------------------------------
3494 !
3495 !...................................
3496  end subroutine aer_property
3497 !-----------------------------------
3498 
3500 ! =======================================================================
3501 ! GOCART code modification starts here (Sarah lu) ---------------------!
3502 !!
3503 !! gocart_init : set_aerspc, rd_gocart_clim, rd_gocart_luts, optavg_grt
3504 !! setgocartaer: aeropt_grt, map_aermr
3505 
3528 !-----------------------------------
3529  subroutine gocart_init &
3530  & ( nwvtot,solfwv,soltot,nwvtir,eirfwv, & ! --- inputs:
3531  & nbdsw,nlwbnd,nswlwbd,imon,me,raddt,fdaer & ! --- outputs: ( none )
3532  & )
3533 
3534 ! ================================================================== !
3535 ! !
3536 ! subprogram : gocart_init !
3537 ! !
3538 ! this is the initialization program for gocart aerosols !
3539 ! !
3540 ! - determine weight and index for aerosol composition/luts !
3541 ! - read in monthly global distribution of gocart aerosols !
3542 ! - read and map the tabulated aerosol optical spectral data !
3543 ! onto corresponding sw/lw radiation spectral bands. !
3544 ! !
3545 ! ==================== defination of variables =================== !
3546 ! !
3547 ! inputs: !
3548 ! NWVTOT - total num of wave numbers used in sw spectrum !
3549 ! solfwv(NWVTOT) - solar flux for each individual wavenumber (w/m2)!
3550 ! soltot - total solar flux for the spectrual range (w/m2)!
3551 ! NWVTIR - total num of wave numbers used in the ir region !
3552 ! eirfwv(NWVTIR) - ir flux(273k) for each individual wavenum (w/m2)!
3553 ! NBDSW - num of bands calculated for sw aeros opt prop !
3554 ! NLWBND - num of bands calculated for lw aeros opt prop !
3555 ! NSWLWBD - total num of bands calc for sw+lw aeros opt prop!
3556 ! imon - month of the year !
3557 ! me - print message control flag !
3558 ! !
3559 ! outputs: (to the module variables) !
3560 ! !
3561 ! module variables: !
3562 ! NBDSW - total number of sw spectral bands !
3563 ! wvnum1,wvnum2 (NSWSTR:NSWEND) !
3564 ! - start/end wavenumbers for each of sw bands !
3565 ! NBDLW - total number of lw spectral bands !
3566 ! wvnlw1,wvnlw2 (NBDLW) !
3567 ! - start/end wavenumbers for each of lw bands !
3568 ! NSWLWBD - total number of sw+lw bands used in this version !
3569 ! extrhi_grt - extinction coef for rh-indep aeros KCM1*NSWLWBD !
3570 ! ssarhi_grt - single-scat-alb for rh-indep aeros KCM1*NSWLWBD !
3571 ! asyrhi_grt - asymmetry factor for rh-indep aeros KCM1*NSWLWBD !
3572 ! extrhd_grt - extinction coef for rh-dep aeros KRHLEV*KCM2*NSWLWBD!
3573 ! ssarhd_grt - single-scat-alb for rh-dep aeros KRHLEV*KCM2*NSWLWBD!
3574 ! asyrhd_grt - asymmetry factor for rh-dep aerosKRHLEV*KCM2*NSWLWBD!
3575 ! ctaer - merging coefficients for fcst/clim fields !
3576 ! get_fcst - option to get fcst aerosol fields !
3577 ! get_clim - option to get clim aerosol fields !
3578 ! dm_indx - index for aer spec to be included in aeropt calculations !
3579 ! dmfcs_indx - index for prognostic aerosol fields !
3580 ! psclmg - geos3/4-gocart pressure IMXG*JMXG*KMXG !
3581 ! dmclmg - geos3-gocart aerosol dry mass IMXG*JMXG*KMXG*NMXG!
3582 ! or geos4-gocart aerosol mixing ratio !
3583 ! !
3584 ! usage: call gocart_init !
3585 ! !
3586 ! subprograms called: set_aerspc, rd_gocart_clim, !
3587 ! rd_gocart_luts, optavg_grt !
3588 ! !
3589 ! ================================================================== !
3590 
3591  implicit none
3592 
3593 ! --- inputs:
3594  integer, intent(in) :: NWVTOT,NWVTIR,NBDSW,NLWBND,NSWLWBD,imon,me
3595 
3596  real (kind=kind_phys), intent(in) :: raddt, fdaer
3597 
3598  real (kind=kind_phys), intent(in) :: solfwv(:),soltot, eirfwv(:)
3599 
3600 ! --- output: ( none )
3601 
3602 ! --- locals:
3603 
3604  real (kind=kind_phys), dimension(NBDSW,KAERBND) :: solwaer
3605  real (kind=kind_phys), dimension(NBDSW) :: solbnd
3606  real (kind=kind_phys), dimension(NLWBND,KAERBND) :: eirwaer
3607  real (kind=kind_phys), dimension(NLWBND) :: eirbnd
3608  real (kind=kind_phys) :: sumsol, sumir
3609 
3610  integer, dimension(NBDSW) :: nv1, nv2
3611  integer, dimension(NLWBND) :: nr1, nr2
3612 
3613  integer :: i, mb, ib, ii, iw, iw1, iw2
3614 
3615 !===> ... begin here
3616 
3617 !--------------------------------------------------------------------------
3618 ! (1) determine aerosol specification index and merging coefficients
3619 !--------------------------------------------------------------------------
3620 
3621  if ( .not. lgrtint ) then
3622 
3623 ! --- ... already done aerspc initialization, continue
3624 
3625  continue
3626 
3627  else
3628 
3629 ! --- ... set aerosol specification index and merging coefficients
3630 
3631  call set_aerspc(raddt,fdaer)
3632 ! --- inputs: (in scope variables)
3633 ! --- outputs: (in scope variables)
3634 
3635  endif ! end if_lgrtinit_block
3636 
3637 !
3638 !--------------------------------------------------------------------------
3639 ! (2) read gocart climatological data
3640 !--------------------------------------------------------------------------
3641 
3642 ! --- ... read gocart climatological data, if needed
3643 
3644  if ( get_clim ) then
3645 
3646  call rd_gocart_clim
3647 ! --- inputs: (in scope variables)
3648 ! --- outputs: (in scope variables)
3649 
3650  endif
3651 
3652 !
3653 !--------------------------------------------------------------------------
3654 ! (3) read and map the tabulated aerosol optical spectral data
3655 ! onto corresponding radiation spectral bands
3656 !--------------------------------------------------------------------------
3657 
3658  if ( .not. lgrtint ) then
3659 
3660 ! --- ... already done optical property interpolation, exit
3661 
3662  return
3663 
3664  else
3665 
3666 ! --- ... reset lgrtint
3667 
3668  lgrtint = .false.
3669 
3670 ! --- ... read tabulated aerosol optical input data
3671  call rd_gocart_luts
3672 ! --- inputs: (in scope variables)
3673 ! --- outputs: (in scope variables)
3674 
3675 ! --- ... compute solar flux weights and interval indices for mapping
3676 ! spectral bands between sw radiation and aerosol data
3677 
3678  solbnd(:) = f_zero
3679  solwaer(:,:) = f_zero
3680 
3681  nv_aod = 1
3682 
3683  do ib = 1, nbdsw
3684  mb = ib + nswstr - 1
3685  ii = 1
3686  iw1 = nint(wvnsw1(mb))
3687  iw2 = nint(wvnsw2(mb))
3688 !
3689 ! --- locate the spectral band for 550nm (for aod diag)
3690 !
3691  if (10000./iw1 >= 0.55 .and. &
3692  & 10000./iw2 <= 0.55 ) then
3693  nv_aod = ib
3694  endif
3695 
3696  lab_swdowhile : do while ( iw1 > iendwv_grt(ii) )
3697  if ( ii == kaerbnd ) exit lab_swdowhile
3698  ii = ii + 1
3699  enddo lab_swdowhile
3700 
3701  sumsol = f_zero
3702  nv1(ib) = ii
3703 
3704  do iw = iw1, iw2
3705  solbnd(ib) = solbnd(ib) + solfwv(iw)
3706  sumsol = sumsol + solfwv(iw)
3707 
3708  if ( iw == iendwv_grt(ii) ) then
3709  solwaer(ib,ii) = sumsol
3710 
3711  if ( ii < kaerbnd ) then
3712  sumsol = f_zero
3713  ii = ii + 1
3714  endif
3715  endif
3716  enddo
3717 
3718  if ( iw2 /= iendwv_grt(ii) ) then
3719  solwaer(ib,ii) = sumsol
3720  endif
3721 
3722  nv2(ib) = ii
3723 
3724  if((me==0) .and. lckprnt) print *,'RAD-nv1,nv2:', &
3725  & ib,iw1,iw2,nv1(ib),iendwv_grt(nv1(ib)), &
3726  & nv2(ib),iendwv_grt(nv2(ib)), &
3727  & 10000./iw1, 10000./iw2
3728  enddo ! end do_ib_block for sw
3729 
3730 ! --- check the spectral range for the nv_550 band
3731  if((me==0) .and. lckprnt) then
3732  mb = nv_aod + nswstr - 1
3733  iw1 = nint(wvnsw1(mb))
3734  iw2 = nint(wvnsw2(mb))
3735  print *,'RAD-nv_aod:', &
3736  & nv_aod, iw1, iw2, 10000./iw1, 10000./iw2
3737  endif
3738 !
3739 ! --- ... compute ir flux weights and interval indices for mapping
3740 ! spectral bands between lw radiation and aerosol data
3741 
3742  eirbnd(:) = f_zero
3743  eirwaer(:,:) = f_zero
3744 
3745  do ib = 1, nlwbnd
3746  ii = 1
3747  if ( nlwbnd == 1 ) then
3748  iw1 = 400 ! corresponding 25 mu
3749  iw2 = 2500 ! corresponding 4 mu
3750  else
3751  iw1 = nint(wvnlw1(ib))
3752  iw2 = nint(wvnlw2(ib))
3753  endif
3754 
3755  lab_lwdowhile : do while ( iw1 > iendwv_grt(ii) )
3756  if ( ii == kaerbnd ) exit lab_lwdowhile
3757  ii = ii + 1
3758  enddo lab_lwdowhile
3759 
3760  sumir = f_zero
3761  nr1(ib) = ii
3762 
3763  do iw = iw1, iw2
3764  eirbnd(ib) = eirbnd(ib) + eirfwv(iw)
3765  sumir = sumir + eirfwv(iw)
3766 
3767  if ( iw == iendwv_grt(ii) ) then
3768  eirwaer(ib,ii) = sumir
3769 
3770  if ( ii < kaerbnd ) then
3771  sumir = f_zero
3772  ii = ii + 1
3773  endif
3774  endif
3775  enddo
3776 
3777  if ( iw2 /= iendwv_grt(ii) ) then
3778  eirwaer(ib,ii) = sumir
3779  endif
3780 
3781  nr2(ib) = ii
3782 
3783  if(me==0 .and. lckprnt) print *,'RAD-nr1,nr2:', &
3784  & ib,iw1,iw2,nr1(ib),iendwv_grt(nr1(ib)), &
3785  & nr2(ib),iendwv_grt(nr2(ib)), &
3786  & 10000./iw1, 10000./iw2
3787  enddo ! end do_ib_block for lw
3788 
3789 ! --- compute spectral band mean properties for each species
3790 
3791  call optavg_grt
3792 ! --- inputs: (in scope variables)
3793 ! --- outputs: (in scope variables)
3794 
3795  if(me==0 .and. lckprnt) then
3796  print *, 'RAD -After optavg_grt, sw band info'
3797  do ib = 1, nbdsw
3798  mb = ib + nswstr - 1
3799  print *,'RAD -wvnsw1,wvnsw2: ',ib,wvnsw1(mb),wvnsw2(mb)
3800  print *,'RAD -lamda1,lamda2: ',ib,10000./wvnsw1(mb), &
3801  & 10000./wvnsw2(mb)
3802  print *,'RAD -extrhi_grt:', extrhi_grt(:,ib)
3803 ! do i = 1, KRHLEV
3804  do i = 1, krhlev, 10
3805  print *, 'RAD -extrhd_grt:',i,rhlev_grt(i), &
3806  & extrhd_grt(i,:,ib)
3807  enddo
3808  enddo
3809  print *, 'RAD -After optavg_grt, lw band info'
3810  do ib = 1, nlwbnd
3811  ii = nbdsw + ib
3812  print *,'RAD -wvnlw1,wvnlw2: ',ib,wvnlw1(ib),wvnlw2(ib)
3813  print *,'RAD -lamda1,lamda2: ',ib,10000./wvnlw1(ib), &
3814  & 10000./wvnlw2(ib)
3815  print *,'RAD -extrhi_grt:', extrhi_grt(:,ii)
3816 ! do i = 1, KRHLEV
3817  do i = 1, krhlev, 10
3818  print *, 'RAD -extrhd_grt:',i,rhlev_grt(i), &
3819  & extrhd_grt(i,:,ii)
3820  enddo
3821  enddo
3822  endif
3823 
3824 ! --- ... dealoocate input data arrays no longer needed
3825  deallocate ( iendwv_grt )
3826  if ( allocated(rhidext0_grt) ) then
3827  deallocate ( rhidext0_grt )
3828  deallocate ( rhidssa0_grt )
3829  deallocate ( rhidasy0_grt )
3830  endif
3831  if ( allocated(rhdpext0_grt) ) then
3832  deallocate ( rhdpext0_grt )
3833  deallocate ( rhdpssa0_grt )
3834  deallocate ( rhdpasy0_grt )
3835  endif
3836 
3837  endif ! end if_lgrtinit_block
3838 
3839 ! =================
3840  contains
3841 ! =================
3842 
3845 !-----------------------------
3846  subroutine set_aerspc(raddt,fdaer)
3847 !.............................
3848 ! --- inputs: (in scope variables)
3849 ! --- outputs: (in scope variables)
3850 
3851 ! ==================================================================== !
3852 ! !
3853 ! subprogram: set_aerspc !
3854 ! !
3855 ! determine merging coefficients ctaer; !
3856 ! set up aerosol specification: num_gridcomp, gridcomp, dm_indx, !
3857 ! dmfcs_indx, isoot, iwaso, isuso, issam, isscm !
3858 ! !
3859 ! Aerosol optical properties (ext, ssa, asy) are determined from !
3860 ! NMGX (<=12) aerosol species !
3861 ! ==> DU: dust1 (4 sub-micron bins), dust2, dust3, dust4, dust5 !
3862 ! BC: soot_phobic, soot_philic !
3863 ! OC: waso_phobic, waso_philic !
3864 ! SU: suso (=so4) !
3865 ! SS: ssam (accumulation mode), sscm (coarse mode) !
3866 ! !
3867 ! The current version only supports prognostic aerosols (from GOCART !
3868 ! in-line calculations) and climo aerosols (from GEOS-GOCART runs) !
3869 ! !
3870 ! ================================================================== !
3871 !
3872  implicit none
3873 
3874 ! --- inputs:
3875  real (kind=kind_phys), intent(in) :: raddt, fdaer
3876 ! --- output:
3877 
3878 ! --- local:
3879 ! real (kind=kind_phys) :: raddt
3880  integer :: i, indxr
3881  character*2 :: tp, gridcomp_tmp(max_num_gridcomp)
3882 
3883 !! ===> determine ctaer (user specified weight for fcst fields)
3884 ! raddt = min(fhswr,fhlwr) / 24.
3885  if( fdaer >= 99999. ) ctaer = f_one
3886  if((fdaer>0.).and.(fdaer<99999.)) ctaer=exp(-raddt/fdaer)
3887 
3888  if(me==0 .and. lckprnt) then
3889  print *, 'RAD -raddt, fdaer,ctaer: ', raddt, fdaer, ctaer
3890  if (ctaer == f_one ) then
3891  print *, 'LU -aerosol fields determined from fcst'
3892  elseif (ctaer == f_zero) then
3893  print *, 'LU -aerosol fields determined from clim'
3894  else
3895  print *, 'LU -aerosol fields determined from fcst/clim'
3896  endif
3897  endif
3898 
3899 !! ===> determine get_fcst and get_clim
3900 !! if fcst is chosen (ctaer == f_one ), set get_clim to F
3901 !! if clim is chosen (ctaer == f_zero), set get_fcst to F
3902  if ( ctaer == f_one ) get_clim = .false.
3903  if ( ctaer == f_zero ) get_fcst = .false.
3904 
3905 !! ===> determine aerosol species to be included in the calculations
3906 !! of aerosol optical properties (ext, ssa, asy)
3907 
3908 !* If climo option is chosen, the aerosol composition is hardwired
3909 !* to full package. If not, the composition is determined from
3910 !* tracer_config on-the-fly (full package or subset)
3911  lab_if_fcst : if ( get_fcst ) then
3912 
3913 !! use tracer_config to determine num_gridcomp and gridcomp
3914  if ( gfs_phy_tracer%doing_GOCART ) then
3915  if ( gfs_phy_tracer%doing_DU ) then
3916  num_gridcomp = num_gridcomp + 1
3917  gridcomp_tmp(num_gridcomp) = 'DU'
3918  endif
3919  if ( gfs_phy_tracer%doing_SU ) then
3920  num_gridcomp = num_gridcomp + 1
3921  gridcomp_tmp(num_gridcomp) = 'SU'
3922  endif
3923  if ( gfs_phy_tracer%doing_SS ) then
3924  num_gridcomp = num_gridcomp + 1
3925  gridcomp_tmp(num_gridcomp) = 'SS'
3926  endif
3927  if ( gfs_phy_tracer%doing_OC ) then
3928  num_gridcomp = num_gridcomp + 1
3929  gridcomp_tmp(num_gridcomp) = 'OC'
3930  endif
3931  if ( gfs_phy_tracer%doing_BC ) then
3932  num_gridcomp = num_gridcomp + 1
3933  gridcomp_tmp(num_gridcomp) = 'BC'
3934  endif
3935 !
3936  if ( num_gridcomp > 0 ) then
3937  allocate ( gridcomp(num_gridcomp) )
3938  gridcomp(1:num_gridcomp) = gridcomp_tmp(1:num_gridcomp)
3939  else
3940  print *,'ERROR: prognostic aerosols not found,abort',me
3941  stop 1000
3942  endif
3943 
3944  else ! gfs_phy_tracer%doing_GOCART=F
3945 
3946  print *,'ERROR: prognostic aerosols option off, abort',me
3947  stop 1001
3948 
3949  endif ! end_if_gfs_phy_tracer%doing_GOCART_if_
3950 
3951  else lab_if_fcst
3952 
3953 !! set to full package (max_num_gridcomp and max_gridcomp)
3954  num_gridcomp = max_num_gridcomp
3955  allocate ( gridcomp(num_gridcomp) )
3956  gridcomp(1:num_gridcomp) = max_gridcomp(1:num_gridcomp)
3957 
3958  endif lab_if_fcst
3959 
3960 !!
3961 !! Aerosol specification is determined as such:
3962 !! A. For radiation-aerosol feedback, the specification is based on the aeropt
3963 !! routine from Mian Chin and Hongbin Yu (hydrophobic and hydrophilic for
3964 !! OC/BC; submicron and supermicron for SS, 8-bins (with 4 subgroups for the
3965 !! the submicron bin) for DU, and SO4 for SU)
3966 !! B. For transport, the specification is determined from GOCART in-line module
3967 !! C. For LUTS, (waso, soot, ssam, sscm, suso, dust) is used, based on the
3968 !! the OPAC climo aerosol scheme (implemented by Yu-Tai Hou)
3969 
3970 !!=== <A> determine dm_indx and NMXG
3971  indxr = 0
3972  dm_indx%waso_phobic = -999 ! OC
3973  dm_indx%soot_phobic = -999 ! BC
3974  dm_indx%ssam = -999 ! SS
3975  dm_indx%suso = -999 ! SU
3976  dm_indx%dust1 = -999 ! DU
3977  do i = 1, num_gridcomp
3978  tp = gridcomp(i)
3979  select case ( tp )
3980  case ( 'OC') ! consider hydrophobic and hydrophilic
3981  dm_indx%waso_phobic = indxr + 1
3982  dm_indx%waso_philic = indxr + 2
3983  indxr = indxr + 2
3984  case ( 'BC') ! consider hydrophobic and hydrophilic
3985  dm_indx%soot_phobic = indxr + 1
3986  dm_indx%soot_philic = indxr + 2
3987  indxr = indxr + 2
3988  case ( 'SS') ! consider submicron and supermicron
3989  dm_indx%ssam = indxr + 1
3990  dm_indx%sscm = indxr + 2
3991  indxr = indxr + 2
3992  case ( 'SU') ! consider SO4 only
3993  dm_indx%suso = indxr + 1
3994  indxr = indxr + 1
3995  case ( 'DU') ! consider all 5 bins
3996  dm_indx%dust1 = indxr + 1
3997  dm_indx%dust2 = indxr + 2
3998  dm_indx%dust3 = indxr + 3
3999  dm_indx%dust4 = indxr + 4
4000  dm_indx%dust5 = indxr + 5
4001  indxr = indxr + 5
4002  case default
4003  print *,'ERROR: aerosol species not supported, abort',me
4004  stop 1002
4005  end select
4006  enddo
4007 !!
4008  nmxg = indxr ! num of gocart aer spec for opt cal
4009 !!
4010 
4011 !!=== <B> determine dmfcs_indx
4012 !! SS: 5-bins are considered for transport while only two groups
4013 !! (accumulation/coarse modes) are considered for radiation
4014 !! DU: 5-bins are considered for transport while 8 bins (with the
4015 !! submicorn bin exptended to 4 bins) are considered for radiation
4016 !! SU: DMS, SO2, and MSA are not considered for radiation
4017 
4018  if ( get_fcst ) then
4019  if ( gfs_phy_tracer%doing_OC ) then
4020  dmfcs_indx%ocphobic = trcindx('ocphobic', gfs_phy_tracer)
4021  dmfcs_indx%ocphilic = trcindx('ocphilic', gfs_phy_tracer)
4022  endif
4023  if ( gfs_phy_tracer%doing_BC ) then
4024  dmfcs_indx%bcphobic = trcindx('bcphobic', gfs_phy_tracer)
4025  dmfcs_indx%bcphilic = trcindx('bcphilic', gfs_phy_tracer)
4026  endif
4027  if ( gfs_phy_tracer%doing_SS ) then
4028  dmfcs_indx%ss001 = trcindx('ss001', gfs_phy_tracer)
4029  dmfcs_indx%ss002 = trcindx('ss002', gfs_phy_tracer)
4030  dmfcs_indx%ss003 = trcindx('ss003', gfs_phy_tracer)
4031  dmfcs_indx%ss004 = trcindx('ss004', gfs_phy_tracer)
4032  dmfcs_indx%ss005 = trcindx('ss005', gfs_phy_tracer)
4033  endif
4034  if ( gfs_phy_tracer%doing_SU ) then
4035  dmfcs_indx%so4 = trcindx('so4', gfs_phy_tracer)
4036  endif
4037  if ( gfs_phy_tracer%doing_DU ) then
4038  dmfcs_indx%du001 = trcindx('du001', gfs_phy_tracer)
4039  dmfcs_indx%du002 = trcindx('du002', gfs_phy_tracer)
4040  dmfcs_indx%du003 = trcindx('du003', gfs_phy_tracer)
4041  dmfcs_indx%du004 = trcindx('du004', gfs_phy_tracer)
4042  dmfcs_indx%du005 = trcindx('du005', gfs_phy_tracer)
4043  endif
4044  endif
4045 
4046 !!
4047 !!=== <C> determin KCM, KCM1, KCM2
4048 !! DU: submicron bin (dust1) contains 4 sub-groups (e.g., hardwire
4049 !! 8 bins for aerosol optical properties luts)
4050 !! OC/BC: while hydrophobic aerosols are rh-independent, the luts
4051 !! for hydrophilic aerosols are used (e.g., use the coeff
4052 !! corresponding to rh=0)
4053 !!
4054  indxr = 1
4055  isoot = -999
4056  iwaso = -999
4057  isuso = -999
4058  issam = -999
4059  isscm = -999
4060  do i = 1, num_gridcomp
4061  tp = gridcomp(i)
4062  if ( tp /= 'DU' ) then
4063  select case ( tp )
4064  case ( 'OC ')
4065  iwaso = indxr
4066  case ( 'BC ')
4067  isoot = indxr
4068  case ( 'SU ')
4069  isuso = indxr
4070  case ( 'SS ')
4071  issam = indxr
4072  isscm = indxr + 1
4073  end select
4074  if ( tp /= 'SS' ) then
4075  indxr = indxr + 1
4076  else
4077  indxr = indxr + 2
4078  endif
4079  else
4080  kcm1 = 8 ! num of rh independent aer species
4081  endif
4082  enddo
4083  kcm2 = indxr - 1 ! num of rh dependent aer species
4084  kcm = kcm1 + kcm2 ! total num of aer species
4085 
4086 !!
4087 !! check print starts here
4088  if( me == 0 .and. lckprnt) then
4089  print *, 'RAD -num_gridcomp:', num_gridcomp
4090  print *, 'RAD -gridcomp :', gridcomp(:)
4091  print *, 'RAD -NMXG:', nmxg
4092  print *, 'RAD -dm_indx ===> '
4093  print *, 'RAD -aerspc: dust1=', dm_indx%dust1
4094  print *, 'RAD -aerspc: dust2=', dm_indx%dust2
4095  print *, 'RAD -aerspc: dust3=', dm_indx%dust3
4096  print *, 'RAD -aerspc: dust4=', dm_indx%dust4
4097  print *, 'RAD -aerspc: dust5=', dm_indx%dust5
4098  print *, 'RAD -aerspc: ssam=', dm_indx%ssam
4099  print *, 'RAD -aerspc: sscm=', dm_indx%sscm
4100  print *, 'RAD -aerspc: suso=', dm_indx%suso
4101  print *, 'RAD -aerspc: waso_phobic=',dm_indx%waso_phobic
4102  print *, 'RAD -aerspc: waso_philic=',dm_indx%waso_philic
4103  print *, 'RAD -aerspc: soot_phobic=',dm_indx%soot_phobic
4104  print *, 'RAD -aerspc: soot_philic=',dm_indx%soot_philic
4105 
4106  print *, 'RAD -KCM1 =', kcm1
4107  print *, 'RAD -KCM2 =', kcm2
4108  print *, 'RAD -KCM =', kcm
4109  if ( kcm2 > 0 ) then
4110  print *, 'RAD -aerspc: issam=', issam
4111  print *, 'RAD -aerspc: isscm=', isscm
4112  print *, 'RAD -aerspc: isuso=', isuso
4113  print *, 'RAD -aerspc: iwaso=', iwaso
4114  print *, 'RAD -aerspc: isoot=', isoot
4115  endif
4116 
4117  if ( get_fcst ) then
4118  print *, 'RAD -dmfcs_indx ===> '
4119  print *, 'RAD -trc_du001=',dmfcs_indx%du001
4120  print *, 'RAD -trc_du002=',dmfcs_indx%du002
4121  print *, 'RAD -trc_du003=',dmfcs_indx%du003
4122  print *, 'RAD -trc_du004=',dmfcs_indx%du004
4123  print *, 'RAD -trc_du005=',dmfcs_indx%du005
4124  print *, 'RAD -trc_so4 =',dmfcs_indx%so4
4125  print *, 'RAD -trc_ocphobic=',dmfcs_indx%ocphobic
4126  print *, 'RAD -trc_ocphilic=',dmfcs_indx%ocphilic
4127  print *, 'RAD -trc_bcphobic=',dmfcs_indx%bcphobic
4128  print *, 'RAD -trc_bcphilic=',dmfcs_indx%bcphilic
4129  print *, 'RAD -trc_ss001=',dmfcs_indx%ss001
4130  print *, 'RAD -trc_ss002=',dmfcs_indx%ss002
4131  print *, 'RAD -trc_ss003=',dmfcs_indx%ss003
4132  print *, 'RAD -trc_ss004=',dmfcs_indx%ss004
4133  print *, 'RAD -trc_ss005=',dmfcs_indx%ss005
4134  endif
4135  endif
4136 !! check print ends here
4137 
4138  return
4139 ! !
4140  end subroutine set_aerspc
4141 
4142 !-----------------------------------
4145 !-----------------------------
4146  subroutine rd_gocart_luts
4147 !.............................
4148 ! --- inputs: (in scope variables)
4149 ! --- outputs: (in scope variables)
4150 
4151 ! ==================================================================== !
4152 ! subprogram: rd_gocart_luts !
4153 ! read input gocart aerosol optical data from Mie code calculations !
4154 ! !
4155 ! Remarks (Quanhua (Mark) Liu, JCSDA, June 2008) !
4156 ! The LUT is for NCEP selected 61 wave numbers and 6 aerosols !
4157 ! (dust, soot, suso, waso, ssam, and sscm) and 36 aerosol effective !
4158 ! size in microns. !
4159 ! !
4160 ! The LUT is computed using Mie code with a logorithm size !
4161 ! distribution for each of 36 effective sizes. The standard deviation !
4162 ! sigma of the size, and min/max size follows Chin et al. 2000 !
4163 ! For each effective size, it corresponds a relative humidity value. !
4164 ! !
4165 ! The LUT contains the density, sigma, relative humidity, mean mode !
4166 ! radius, effective size, mass extinction coefficient, single !
4167 ! scattering albedo, asymmetry factor, and phase function !
4168 ! !
4169 ! ================================================================== !
4170 !
4171  implicit none
4172 
4173 ! --- inputs:
4174 ! --- output:
4175 
4176 ! --- locals:
4177  INTEGER, PARAMETER :: NP = 100, np2 = 2*np, nwave=100, &
4178  & naero=6, n_p=36
4179  INTEGER :: NW, NS, nH, n_bin
4180  real (kind=kind_io8), Dimension( NP2 ) :: Angle, Cos_Angle, &
4181  & Cos_Weight
4182  real (kind=kind_io8), Dimension(n_p,nAero) :: RH, rm, reff
4183  real (kind=kind_io8), Dimension(nWave,n_p,nAero) :: &
4184  & ext0, sca0, asy0
4185  real (kind=kind_io8), Dimension(NP2,n_p,nWave,nAero) :: ph0
4186  real (kind=kind_io8) :: wavelength(nwave), density(naero), &
4187  & sigma(nAero), wave,n_fac,PI,t1,s1,g1
4188  CHARACTER(len=80) :: AerosolName(naero)
4189  INTEGER :: i, j, k, l, ij
4190 
4191  character :: aerosol_file*30
4192  logical :: file_exist
4193  integer :: indx_dust(8) ! map 36 dust bins to gocart size bins
4194 
4195  data aerosol_file /"NCEP_AEROSOL.bin"/
4196  data aerosolname/ ' Dust ', ' Soot ', ' SUSO ', ' WASO ', &
4197  & ' SSAM ', ' SSCM '/
4198 
4199 !! 8 dust bins
4200 !! 1 2 3 4 5 6 7 8
4201 !! .1-.18, .18-.3, .3-.6, 0.6-1.0, 1.0-1.8, 1.8-3, 3-6, 6-10 <-- def
4202 !! 0.1399 0.2399 0.4499 0.8000 1.3994 2.3964 4.4964 7.9887 <-- reff
4203  data indx_dust/4, 8, 12, 18, 21, 24, 30, 36/
4204 
4205  pi = acos(-1.d0)
4206 
4207 ! -- allocate aerosol optical data
4208  if ( .not. allocated( iendwv_grt ) ) then
4209  allocate ( iendwv_grt(kaerbnd) )
4210  endif
4211  if (.not. allocated(rhidext0_grt) .and. kcm1 > 0 ) then
4212  allocate ( rhidext0_grt(kaerbnd,kcm1))
4213  allocate ( rhidssa0_grt(kaerbnd,kcm1))
4214  allocate ( rhidasy0_grt(kaerbnd,kcm1))
4215  endif
4216  if (.not. allocated(rhdpext0_grt) .and. kcm2 > 0 ) then
4217  allocate ( rhdpext0_grt(kaerbnd,krhlev,kcm2))
4218  allocate ( rhdpssa0_grt(kaerbnd,krhlev,kcm2))
4219  allocate ( rhdpasy0_grt(kaerbnd,krhlev,kcm2))
4220  endif
4221 
4222 ! -- read luts
4223  inquire (file = aerosol_file, exist = file_exist)
4224 
4225  if ( file_exist ) then
4226  if(me==0 .and. lckprnt) print *,'RAD -open :',aerosol_file
4227  close (niaercm)
4228  open (unit=niaercm,file=aerosol_file,status='OLD', &
4229  & action='read',form='UNFORMATTED')
4230  else
4231  print *,' Requested aerosol data file "',aerosol_file, &
4232  & '" not found!', me
4233  print *,' *** Stopped in subroutine RD_GOCART_LUTS !!'
4234  stop 1003
4235  endif ! end if_file_exist_block
4236 
4237  READ(niaercm) (cos_angle(i),i=1,np)
4238  READ(niaercm) (cos_weight(i),i=1,np)
4239  READ(niaercm)
4240  READ(niaercm)
4241  READ(niaercm) nw,ns
4242  READ(niaercm)
4243  READ(niaercm) (wavelength(i),i=1,nw)
4244 
4245 ! --- check nAero and NW
4246  if (nw /= kaerbnd) then
4247  print *, "Incorrect spectral band, abort ", nw
4248  stop 1004
4249  endif
4250 
4251 ! --- convert wavelength to wavenumber
4252  do i = 1, kaerbnd
4253  iendwv_grt(i) = 10000. / wavelength(i)
4254  if(me==0 .and. lckprnt) print *,'RAD -wn,lamda:', &
4255  & i,iendwv_grt(i),wavelength(i)
4256  enddo
4257 
4258  DO j = 1, naero
4259  if(me==0 .and. lckprnt) print *,'RAD -read LUTs:', &
4260  & j,aerosolname(j)
4261  READ(niaercm)
4262  READ(niaercm)
4263  READ(niaercm) n_bin, density(j), sigma(j)
4264  READ(niaercm)
4265  READ(niaercm) (rh(i,j),i=1, n_bin)
4266  READ(niaercm)
4267  READ(niaercm) (rm(i,j),i=1, n_bin)
4268  READ(niaercm)
4269  READ(niaercm) (reff(i,j),i=1, n_bin)
4270 
4271 ! --- check n_bin
4272  if (n_bin /= krhlev ) then
4273  print *, "Incorrect rh levels, abort ", n_bin
4274  stop 1005
4275  endif
4276 
4277 ! --- read luts
4278  DO k = 1, nw
4279  READ(niaercm) wave,(ext0(k,l,j),l=1,n_bin)
4280  READ(niaercm) (sca0(k,l,j),l=1,n_bin)
4281  READ(niaercm) (asy0(k,l,j),l=1,n_bin)
4282  READ(niaercm) (ph0(1:np2,l,k,j),l=1,n_bin)
4283  END DO
4284 
4285 ! --- map luts input to module variables
4286  if (aerosolname(j) == ' Dust ' ) then
4287  if ( kcm1 > 0) then
4288  do i = 1, kcm1
4289  rhidext0_grt(1:kaerbnd,i)=ext0(1:kaerbnd,indx_dust(i),j)
4290  rhidssa0_grt(1:kaerbnd,i)=sca0(1:kaerbnd,indx_dust(i),j)
4291  rhidasy0_grt(1:kaerbnd,i)=asy0(1:kaerbnd,indx_dust(i),j)
4292  enddo
4293  endif
4294  else
4295  if ( kcm2 > 0) then
4296  if (aerosolname(j) == ' Soot ') ij = isoot
4297  if (aerosolname(j) == ' SUSO ') ij = isuso
4298  if (aerosolname(j) == ' WASO ') ij = iwaso
4299  if (aerosolname(j) == ' SSAM ') ij = issam
4300  if (aerosolname(j) == ' SSCM ') ij = isscm
4301  if ( ij .ne. -999 ) then
4302  rhdpext0_grt(1:kaerbnd,1:krhlev,ij) = &
4303  & ext0(1:kaerbnd,1:krhlev,j)
4304  rhdpssa0_grt(1:kaerbnd,1:krhlev,ij) = &
4305  & sca0(1:kaerbnd,1:krhlev,j)
4306  rhdpasy0_grt(1:kaerbnd,1:krhlev,ij) = &
4307  & asy0(1:kaerbnd,1:krhlev,j)
4308  endif ! if_ij
4309  endif ! if_KCM2
4310  endif
4311  END DO
4312 
4313  return
4314 !...................................
4315  end subroutine rd_gocart_luts
4316 !-----------------------------------
4317 ! !
4322 !-----------------------------
4323  subroutine optavg_grt
4324 !.............................
4325 ! --- inputs: (in scope variables)
4326 ! --- outputs: (in scope variables)
4327 
4328 ! ==================================================================== !
4329 ! !
4330 ! subprogram: optavg_grt !
4331 ! !
4332 ! compute mean aerosols optical properties over each sw/lw radiation !
4333 ! spectral band for each of the species components. This program !
4334 ! follows gfdl's approach for thick cloud opertical property in !
4335 ! sw radiation scheme (2000). !
4336 ! !
4337 ! ==================== defination of variables =================== !
4338 ! !
4339 ! input arguments: !
4340 ! nv1,nv2 (NBDSW) - start/end spectral band indices of aerosol data !
4341 ! for each sw radiation spectral band !
4342 ! nr1,nr2 (NLWBND) - start/end spectral band indices of aerosol data !
4343 ! for each ir radiation spectral band !
4344 ! solwaer (NBDSW,KAERBND) !
4345 ! - solar flux weight over each sw radiation band !
4346 ! vs each aerosol data spectral band !
4347 ! eirwaer (NLWBND,KAERBND) !
4348 ! - ir flux weight over each lw radiation band !
4349 ! vs each aerosol data spectral band !
4350 ! solbnd (NBDSW) - solar flux weight over each sw radiation band !
4351 ! eirbnd (NLWBND) - ir flux weight over each lw radiation band !
4352 ! NBDSW - total number of sw spectral bands !
4353 ! NLWBND - total number of lw spectral bands !
4354 ! NSWLWBD - total number of sw+lw spectral bands !
4355 ! !
4356 ! output arguments: (to module variables) !
4357 ! !
4358 ! ================================================================== !
4359 !
4360  implicit none
4361 
4362 ! --- inputs:
4363 ! --- output:
4364 
4365 ! --- locals:
4366  real (kind=kind_phys) :: sumk, sumok, sumokg, sumreft, &
4367  & sp, refb, reft, rsolbd, rirbd
4368 
4369  integer :: ib, nb, ni, nh, nc
4370 !
4371 !===> ... begin here
4372 
4373 ! --- ... allocate aerosol optical data
4374  if (.not. allocated(extrhd_grt) .and. kcm2 > 0 ) then
4375  allocate ( extrhd_grt(krhlev,kcm2,nswlwbd) )
4376  allocate ( ssarhd_grt(krhlev,kcm2,nswlwbd) )
4377  allocate ( asyrhd_grt(krhlev,kcm2,nswlwbd) )
4378  endif
4379  if (.not. allocated(extrhi_grt) .and. kcm1 > 0 ) then
4380  allocate ( extrhi_grt(kcm1,nswlwbd) )
4381  allocate ( ssarhi_grt(kcm1,nswlwbd) )
4382  allocate ( asyrhi_grt(kcm1,nswlwbd) )
4383  endif
4384 !
4385 ! --- ... loop for each sw radiation spectral band
4386 
4387  do nb = 1, nbdsw
4388  rsolbd = f_one / solbnd(nb)
4389 
4390 ! --- for rh independent aerosol species
4391 
4392  lab_rhi: if (kcm1 > 0 ) then
4393  do nc = 1, kcm1
4394  sumk = f_zero
4395  sumok = f_zero
4396  sumokg = f_zero
4397  sumreft = f_zero
4398 
4399  do ni = nv1(nb), nv2(nb)
4400  sp = sqrt( (f_one - rhidssa0_grt(ni,nc)) &
4401  & / (f_one - rhidssa0_grt(ni,nc)*rhidasy0_grt(ni,nc)) )
4402  reft = (f_one - sp) / (f_one + sp)
4403  sumreft = sumreft + reft*solwaer(nb,ni)
4404 
4405  sumk = sumk + rhidext0_grt(ni,nc)*solwaer(nb,ni)
4406  sumok = sumok + rhidssa0_grt(ni,nc)*solwaer(nb,ni) &
4407  & * rhidext0_grt(ni,nc)
4408  sumokg = sumokg + rhidssa0_grt(ni,nc)*solwaer(nb,ni) &
4409  & * rhidext0_grt(ni,nc)*rhidasy0_grt(ni,nc)
4410  enddo
4411 
4412  refb = sumreft * rsolbd
4413 
4414  extrhi_grt(nc,nb) = sumk * rsolbd
4415  asyrhi_grt(nc,nb) = sumokg / (sumok + 1.0e-10)
4416  ssarhi_grt(nc,nb) = 4.0*refb &
4417  & / ( (f_one+refb)**2 - asyrhi_grt(nc,nb)*(f_one-refb)**2 )
4418 
4419  enddo ! end do_nc_block for rh-ind aeros
4420  endif lab_rhi
4421 
4422 ! --- for rh dependent aerosols species
4423 
4424  lab_rhd: if (kcm2 > 0 ) then
4425  do nc = 1, kcm2
4426  do nh = 1, krhlev
4427  sumk = f_zero
4428  sumok = f_zero
4429  sumokg = f_zero
4430  sumreft = f_zero
4431 
4432  do ni = nv1(nb), nv2(nb)
4433  sp = sqrt( (f_one - rhdpssa0_grt(ni,nh,nc)) &
4434  & /(f_one-rhdpssa0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)))
4435  reft = (f_one - sp) / (f_one + sp)
4436  sumreft = sumreft + reft*solwaer(nb,ni)
4437 
4438  sumk = sumk + rhdpext0_grt(ni,nh,nc)*solwaer(nb,ni)
4439  sumok = sumok + rhdpssa0_grt(ni,nh,nc)*solwaer(nb,ni) &
4440  & * rhdpext0_grt(ni,nh,nc)
4441  sumokg = sumokg + rhdpssa0_grt(ni,nh,nc)*solwaer(nb,ni) &
4442  & * rhdpext0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)
4443  enddo
4444 
4445  refb = sumreft * rsolbd
4446 
4447  extrhd_grt(nh,nc,nb) = sumk * rsolbd
4448  asyrhd_grt(nh,nc,nb) = sumokg / (sumok + 1.0e-10)
4449  ssarhd_grt(nh,nc,nb) = 4.0*refb &
4450  & /((f_one+refb)**2 - asyrhd_grt(nh,nc,nb)*(f_one-refb)**2)
4451  enddo ! end do_nh_block
4452  enddo ! end do_nc_block for rh-dep aeros
4453  endif lab_rhd
4454 
4455  enddo ! end do_nb_block for sw
4456 
4457 ! --- ... loop for each lw radiation spectral band
4458 
4459  do nb = 1, nlwbnd
4460 
4461  ib = nbdsw + nb
4462  rirbd = f_one / eirbnd(nb)
4463 
4464 ! --- for rh independent aerosol species
4465 
4466  lab_rhi_lw: if (kcm1 > 0 ) then
4467  do nc = 1, kcm1
4468  sumk = f_zero
4469  sumok = f_zero
4470  sumokg = f_zero
4471  sumreft = f_zero
4472 
4473  do ni = nr1(nb), nr2(nb)
4474  sp = sqrt( (f_one - rhidssa0_grt(ni,nc)) &
4475  & / (f_one - rhidssa0_grt(ni,nc)*rhidasy0_grt(ni,nc)) )
4476  reft = (f_one - sp) / (f_one + sp)
4477  sumreft = sumreft + reft*eirwaer(nb,ni)
4478 
4479  sumk = sumk + rhidext0_grt(ni,nc)*eirwaer(nb,ni)
4480  sumok = sumok + rhidssa0_grt(ni,nc)*eirwaer(nb,ni) &
4481  & * rhidext0_grt(ni,nc)
4482  sumokg = sumokg + rhidssa0_grt(ni,nc)*eirwaer(nb,ni) &
4483  & * rhidext0_grt(ni,nc)*rhidasy0_grt(ni,nc)
4484  enddo
4485 
4486  refb = sumreft * rirbd
4487 
4488  extrhi_grt(nc,ib) = sumk * rirbd
4489  asyrhi_grt(nc,ib) = sumokg / (sumok + 1.0e-10)
4490  ssarhi_grt(nc,ib) = 4.0*refb &
4491  & / ( (f_one+refb)**2 - asyrhi_grt(nc,ib)*(f_one-refb)**2 )
4492  enddo ! end do_nc_block for rh-ind aeros
4493  endif lab_rhi_lw
4494 
4495 ! --- for rh dependent aerosols species
4496 
4497  lab_rhd_lw: if (kcm2 > 0 ) then
4498  do nc = 1, kcm2
4499  do nh = 1, krhlev
4500  sumk = f_zero
4501  sumok = f_zero
4502  sumokg = f_zero
4503  sumreft = f_zero
4504 
4505  do ni = nr1(nb), nr2(nb)
4506  sp = sqrt( (f_one - rhdpssa0_grt(ni,nh,nc)) &
4507  & /(f_one - rhdpssa0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)) )
4508  reft = (f_one - sp) / (f_one + sp)
4509  sumreft = sumreft + reft*eirwaer(nb,ni)
4510 
4511  sumk = sumk + rhdpext0_grt(ni,nh,nc)*eirwaer(nb,ni)
4512  sumok = sumok + rhdpssa0_grt(ni,nh,nc)*eirwaer(nb,ni) &
4513  & * rhdpext0_grt(ni,nh,nc)
4514  sumokg = sumokg+ rhdpssa0_grt(ni,nh,nc)*eirwaer(nb,ni) &
4515  & * rhdpext0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)
4516  enddo
4517 
4518  refb = sumreft * rirbd
4519 
4520  extrhd_grt(nh,nc,ib) = sumk * rirbd
4521  asyrhd_grt(nh,nc,ib) = sumokg / (sumok + 1.0e-10)
4522  ssarhd_grt(nh,nc,ib) = 4.0*refb &
4523  & /((f_one+refb)**2 - asyrhd_grt(nh,nc,ib)*(f_one-refb)**2 )
4524  enddo ! end do_nh_block
4525  enddo ! end do_nc_block for rh-dep aeros
4526  endif lab_rhd_lw
4527 
4528  enddo ! end do_nb_block for lw
4529 
4530 !
4531  return
4532 !................................
4533  end subroutine optavg_grt
4534 !--------------------------------
4535 !
4541 !-----------------------------------
4542  subroutine rd_gocart_clim
4543 !...................................
4544 ! --- inputs: (in scope variables)
4545 ! --- outputs: (in scope variables)
4546 
4547 ! ================================================================== !
4548 ! !
4549 ! subprogram: rd_gocart_clim !
4550 ! !
4551 ! 1. read in aerosol dry mass and surface pressure from GEOS3-GOCART !
4552 ! C3.1 2000 monthly data set !
4553 ! or aerosol mixing ratio and surface pressure from GEOS4-GOCART !
4554 ! 2000-2007 averaged monthly data set !
4555 ! 2. compute goes lat/lon array (for horizontal mapping) !
4556 ! !
4557 ! ==================== defination of variables =================== !
4558 ! !
4559 ! inputs arguments: !
4560 ! imon - month of the year !
4561 ! me - print message control flag !
4562 ! !
4563 ! outputs arguments: (to the module variables) !
4564 ! psclmg - pressure (sfc to toa) cb IMXG*JMXG*KMXG !
4565 ! dmclmg - aerosol dry mass/mixing ratio IMXG*JMXG*KMXG*NMXG !
4566 ! geos_rlon - goes longitude deg IMXG !
4567 ! geos_rlat - goes latitude deg JMXG !
4568 ! !
4569 ! usage: call rd_gocart_clim !
4570 ! !
4571 ! program history: !
4572 ! 05/18/2010 --- Lu Add the option to read GEOS4-GOCART climo !
4573 ! ================================================================== !
4574 !
4575  implicit none
4576 
4577 ! --- inputs:
4578 ! --- output:
4579 
4580 ! --- locals:
4581  integer, parameter :: MAXSPC = 5
4582  real (kind=kind_io4), parameter :: PINT = 0.01
4583  real (kind=kind_io4), parameter :: EPSQ = 0.0
4584 
4585  integer :: i, j, k, numspci, ii
4586  integer :: icmp, nrecl, nt1, nt2, nn(maxspc)
4587  character :: ymd*6, yr*4, mn*2, tp*2, &
4588  & fname*30, fin*30, aerosol_file*40
4589  logical :: file_exist
4590 
4591  real (kind=kind_io4), dimension(KMXG) :: sig
4592  real (kind=kind_io4), dimension(IMXG,JMXG) :: ps
4593  real (kind=kind_io4), dimension(IMXG,JMXG,KMXG) :: temp
4594  real (kind=kind_io4), dimension(IMXG,JMXG,KMXG,MAXSPC):: buff
4595  real (kind=kind_phys) :: pstmp
4596 
4597 ! Add the following variables for GEOS4-GOCART
4598  real (kind=kind_io4), dimension(KMXG):: hyam, hybm
4599  real (kind=kind_io4) :: p0
4600 
4601  data yr /'2000'/ !!<=== use 2000 as the climo proxy
4602 
4603 !* sigma_coordinate for GEOS3-GOCART
4604 !* P(i,j,k) = PINT + SIG(k) * (PS(i,j) - PINT)
4605  data sig / &
4606  & 9.98547e-01,9.94147e-01,9.86350e-01,9.74300e-01,9.56950e-01, &
4607  & 9.33150e-01,9.01750e-01,8.61500e-01,8.11000e-01,7.50600e-01, &
4608  & 6.82900e-01,6.10850e-01,5.37050e-01,4.63900e-01,3.93650e-01, &
4609  & 3.28275e-01,2.69500e-01,2.18295e-01,1.74820e-01,1.38840e-01, &
4610  & 1.09790e-01,8.66900e-02,6.84150e-02,5.39800e-02,4.25750e-02, &
4611  & 3.35700e-02,2.39900e-02,1.36775e-02,5.01750e-03,5.30000e-04 /
4612 
4613 !* hybrid_sigma_pressure_coordinate for GEOS4-GOCART
4614 !* p(i,j,k) = a(k)*p0 + b(k)*ps(i,j)
4615  data hyam/ &
4616  & 0, 0.0062694, 0.02377049, 0.05011813, 0.08278809, 0.1186361, &
4617  & 0.1540329, 0.1836373, 0.2043698, 0.2167788, 0.221193, &
4618  & 0.217729, 0.2062951, 0.1865887, 0.1615213, 0.1372958, &
4619  & 0.1167039, 0.09920014, 0.08432171, 0.06656809, 0.04765031, &
4620  & 0.03382346, 0.0237648, 0.01435208, 0.00659734, 0.002826232, &
4621  & 0.001118959, 0.0004086494, 0.0001368611, 3.750308e-05/
4622 
4623  data hybm / &
4624  & 0.992555, 0.9642, 0.90556, 0.816375, 0.703815, 0.576585, &
4625  & 0.44445, 0.324385, 0.226815, 0.149165, 0.089375, &
4626  & 0.045865, 0.017485, 0.00348, 0, 0, 0, 0, 0, &
4627  & 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 /
4628 
4629  data p0 /1013.25 /
4630 
4631 !===> ... begin here
4632 
4633 ! --- allocate and initialize gocart climatological data
4634  if ( .not. allocated (dmclmg) ) then
4635  allocate ( dmclmg(imxg,jmxg,kmxg,nmxg) )
4636  allocate ( psclmg(imxg,jmxg,kmxg) )
4637  allocate ( molwgt(nmxg) )
4638  endif
4639 
4640  dmclmg(:,:,:,:) = f_zero
4641  psclmg(:,:,:) = f_zero
4642  molwgt(:) = f_zero
4643 
4644 ! --- allocate and initialize geos lat and lon arrays
4645  if ( .not. allocated ( geos_rlon )) then
4646  allocate (geos_rlon(imxg))
4647  allocate (geos_rlat(jmxg))
4648  endif
4649 
4650  geos_rlon(:) = f_zero
4651  geos_rlat(:) = f_zero
4652 
4653 ! --- compute geos lat and lon arrays
4654  do i = 1, imxg
4655  geos_rlon(i) = -180. + (i-1)* dltx
4656  end do
4657  do j = 2, jmxg-1
4658  geos_rlat(j) = -90. + (j-1)* dlty
4659  end do
4660  geos_rlat(1) = -89.5
4661  geos_rlat(jmxg) = 89.5
4662 
4663 ! --- determine whether GEOS3 or GEOS4 data set is provided
4664  if ( gocart_climo == 'xxxx' ) then
4665  gocart_climo='0000'
4666 ! check geos3-gocart climo
4667  aerosol_file = '200001.PS.avg'
4668  inquire (file = aerosol_file, exist = file_exist)
4669  if ( file_exist ) gocart_climo='ver3'
4670 ! check geos4-gocart climo
4671  aerosol_file = 'gocart_climo_2000x2007_ps_01.bin'
4672  inquire (file = aerosol_file, exist = file_exist)
4673  if ( file_exist ) gocart_climo='ver4'
4674  endif
4675 !
4676 !
4677 ! --- read ps (sfc pressure) and compute 3d pressure field (psclmg)
4678 !
4679  write(mn,'(i2.2)') imon
4680  ymd = yr//mn
4681  aerosol_file = 'null'
4682  if ( gocart_climo == 'ver3' ) then
4683  aerosol_file = ymd//'.PS.avg'
4684  elseif ( gocart_climo == 'ver4' ) then
4685  aerosol_file = 'gocart_climo_2000x2007_ps_'//mn//'.bin'
4686  endif
4687 !
4688  inquire (file = aerosol_file, exist = file_exist)
4689  lab_if_ps : if ( file_exist ) then
4690 
4691  close(niaercm)
4692  if ( gocart_climo == 'ver3' ) then
4693  nrecl = 4 * (imxg * jmxg)
4694  open(niaercm, file=trim(aerosol_file), &
4695  & action='read',access='direct',recl=nrecl)
4696  read(niaercm, rec=1) ps
4697  do j = 1, jmxg
4698  do i = 1, imxg
4699  do k = 1, kmxg
4700  pstmp = pint + sig(k) * (ps(i,j) - pint)
4701  psclmg(i,j,k) = 0.1 * pstmp ! convert mb to cb
4702  enddo
4703  enddo
4704  enddo
4705 
4706  elseif ( gocart_climo == 'ver4' ) then
4707  open(niaercm, file=trim(aerosol_file), &
4708  & action='read',status='old', form='unformatted')
4709  read(niaercm) ps(:,:)
4710  do j = 1, jmxg
4711  do i = 1, imxg
4712  do k = 1, kmxg
4713  pstmp = hyam(k)*p0 + hybm(k)*ps(i,j)
4714  psclmg(i,j,k) = 0.1 * pstmp ! convert mb to cb
4715  enddo
4716  enddo
4717  enddo
4718 
4719  endif ! ---- end if_gocart_climo
4720 
4721  else lab_if_ps
4722 
4723  print *,' *** Requested aerosol data file "', &
4724  & trim(aerosol_file), '" not found!'
4725  print *,' *** Stopped in RD_GOCART_CLIM ! ', me
4726  stop 1006
4727  endif lab_if_ps
4728 !
4729 ! --- read aerosol dry mass (g/m3) or mixing ratios (mol/mol,kg/kg)
4730 !
4731  lab_do_icmp : do icmp = 1, num_gridcomp
4732 
4733  tp = gridcomp(icmp)
4734 
4735 ! determine aerosol_file
4736  aerosol_file = 'null'
4737  if ( gocart_climo == 'ver3' ) then
4738  if(tp == 'DU') fname='.DU.STD.tv20.g.avg'
4739  if(tp == 'SS') fname='.SS.STD.tv17.g.avg'
4740  if(tp == 'SU') fname='.SU.STD.tv15.g.avg'
4741  if(tp == 'OC') fname='.CC.STD.tv15.g.avg'
4742  if(tp == 'BC') fname='.CC.STD.tv15.g.avg'
4743  aerosol_file=ymd//trim(fname)
4744  elseif ( gocart_climo == 'ver4' ) then
4745  fin = 'gocart_climo_2000x2007_'
4746  if(tp == 'DU') fname=trim(fin)//'du_'
4747  if(tp == 'SS') fname=trim(fin)//'ss_'
4748  if(tp == 'SU') fname=trim(fin)//'su_'
4749  if(tp == 'OC') fname=trim(fin)//'cc_'
4750  if(tp == 'BC') fname=trim(fin)//'cc_'
4751  aerosol_file=trim(fname)//mn//'.bin'
4752  endif
4753 
4754  numspci = 4
4755  if(tp == 'DU') numspci = 5
4756  inquire (file=trim(aerosol_file), exist = file_exist)
4757  lab_if_aer: if ( file_exist ) then
4758 !
4759  close(niaercm)
4760  if ( gocart_climo == 'ver3' ) then
4761  nrecl = 4 * numspci * (imxg * jmxg * kmxg + 3)
4762  open (niaercm, file=trim(aerosol_file), &
4763  & action='read',access='direct', recl=nrecl)
4764  read(niaercm,rec=1)(nt1,nt2,nn(i),buff(:,:,:,i),i=1,numspci)
4765 
4766  elseif ( gocart_climo == 'ver4' ) then
4767  open (niaercm, file=trim(aerosol_file), &
4768  & action='read',status='old', form='unformatted')
4769  do i = 1, numspci
4770  do k = 1, kmxg
4771  read(niaercm) temp(:,:,k)
4772  buff(:,:,k,i) = temp(:,:,k)
4773  enddo
4774  enddo
4775  endif
4776 
4777 !!===> fill dmclmg with working array buff
4778  select case ( tp )
4779 
4780 ! fill in DU from DU: du1, du2, du3, du4, du5
4781  case ('DU' )
4782  if ( dm_indx%dust1 /= -999) then
4783  do ii = 1, 5
4784  dmclmg(:,:,:,dm_indx%dust1+ii-1) = buff(:,:,:,ii)
4785  enddo
4786  else
4787  print *, 'ERROR: invalid DU index, abort! ',me
4788  stop 1007
4789  endif
4790 
4791 ! fill in BC from CC: bc_phobic, oc_phobic, bc_philic, oc_philic
4792  case ('BC' )
4793  if ( dm_indx%soot_phobic /= -999) then
4794  dmclmg(:,:,:,dm_indx%soot_phobic)=buff(:,:,:,1)
4795  dmclmg(:,:,:,dm_indx%soot_philic)=buff(:,:,:,3)
4796  molwgt(dm_indx%soot_phobic) = 12.
4797  molwgt(dm_indx%soot_philic) = 12.
4798  else
4799  print *, 'ERROR: invalid BC index, abort! ',me
4800  stop 1008
4801  endif
4802 
4803 ! fill in SU from SU: dms, so2, so4, msa
4804  case ('SU' )
4805  if ( dm_indx%suso /= -999) then
4806  dmclmg(:,:,:,dm_indx%suso) = buff(:,:,:,3)
4807  molwgt(dm_indx%suso) = 96.
4808  else
4809  print *, 'ERROR: invalid SU index, abort! ',me
4810  stop 1009
4811  endif
4812 
4813 ! fill in OC from CC: bc_phobic, oc_phobic, bc_philic, oc_philic
4814  case ('OC' )
4815  if ( dm_indx%waso_phobic /= -999) then
4816  dmclmg(:,:,:,dm_indx%waso_phobic) = 1.4*buff(:,:,:,2)
4817  dmclmg(:,:,:,dm_indx%waso_philic) = 1.4*buff(:,:,:,4)
4818  molwgt(dm_indx%waso_phobic) = 12.
4819  molwgt(dm_indx%waso_philic) = 12.
4820  else
4821  print *, 'ERROR: invalid OC index, abort! ',me
4822  stop 1010
4823  endif
4824 
4825 ! fill in SS from SS: ss1, ss2, ss3, ss4
4826  case ('SS' )
4827  if ( dm_indx%ssam /= -999) then
4828  dmclmg(:,:,:,dm_indx%ssam) = buff(:,:,:,1)
4829  dmclmg(:,:,:,dm_indx%sscm) = buff(:,:,:,2) + &
4830  & buff(:,:,:,3)+buff(:,:,:,4)
4831  else
4832  print *, 'ERROR: invalid SS index, abort! ',me
4833  stop 1011
4834  endif
4835 
4836  case default
4837 
4838  print *, 'ERROR: invalid aerosol species, abort ',tp
4839  stop 1012
4840 
4841  end select
4842 
4843  else lab_if_aer
4844  print *,' *** Requested aerosol data file "',aerosol_file, &
4845  & '" not found!'
4846  print *,' *** Stopped in RD_GOCART_CLIM ! ', me
4847  stop 1013
4848  endif lab_if_aer
4849 
4850  enddo lab_do_icmp
4851 
4852  return
4853 !...................................
4854  end subroutine rd_gocart_clim
4855 !-----------------------------------
4856 !
4857 !...................................
4858  end subroutine gocart_init
4859 !-----------------------------------
4860 !! @}
4861 
4891 !-----------------------------------
4892  subroutine setgocartaer &
4893  & ( alon,alat,prslk,rhlay,dz,hz,nswlwbd, & ! --- inputs:
4894  & prsl,tvly,trcly, &
4895  & imax,nlay,nlp1, ivflip, lsswr,lslwr, &
4896  & aerosw,aerolw & ! --- outputs:
4897  & )
4898 
4899 
4900 ! ================================================================== !
4901 ! !
4902 ! setgocartaer computes sw + lw aerosol optical properties for gocart !
4903 ! aerosol species (merged from fcst and clim fields) !
4904 ! !
4905 ! inputs: !
4906 ! alon, alat IMAX !
4907 ! - longitude and latitude of given points in degree !
4908 ! prslk - pressure cb IMAX*NLAY !
4909 ! rhlay - layer mean relative humidity IMAX*NLAY !
4910 ! dz - layer thickness m IMAX*NLAY !
4911 ! hz - level high m IMAX*NLP1 !
4912 ! NSWLWBD - total number of sw+ir bands for aeros opt prop 1 !
4913 ! prsl - layer mean pressure mb IMAX*NLAY !
4914 ! tvly - layer mean virtual temperature k IMAX*NLAY !
4915 ! trcly - layer mean specific tracer g/g IMAX*NLAY*NTRAC!
4916 ! IMAX - horizontal dimension of arrays 1 !
4917 ! NLAY,NLP1-vertical dimensions of arrays 1 !
4918 ! ivflip - control flag for direction of vertical index 1 !
4919 ! =0: index from toa to surface !
4920 ! =1: index from surface to toa !
4921 ! lsswr,lslwr !
4922 ! - logical flag for sw/lw radiation calls 1 !
4923 ! !
4924 ! outputs: !
4925 ! aerosw - aeros opt properties for sw IMAX*NLAY*NBDSW*NF_AESW!
4926 ! (:,:,:,1): optical depth !
4927 ! (:,:,:,2): single scattering albedo !
4928 ! (:,:,:,3): asymmetry parameter !
4929 ! aerolw - aeros opt properties for lw IMAX*NLAY*NBDLW*NF_AELW!
4930 ! (:,:,:,1): optical depth !
4931 ! (:,:,:,2): single scattering albedo !
4932 ! (:,:,:,3): asymmetry parameter !
4933 ! tau_gocart - 550nm aeros opt depth IMAX*NLAY*MAX_NUM_GRIDCOMP!
4934 ! !
4935 ! module parameters and constants: !
4936 ! NBDSW - total number of sw bands for aeros opt prop 1 !
4937 ! NLWBND - total number of ir bands for aeros opt prop 1 !
4938 ! !
4939 ! module variable: (set by subroutine gocart_init) !
4940 ! dmclmg - aerosols dry mass/mixing ratios IMXG*JMXG*KMXG*NMXG !
4941 ! psclmg - pressure cb IMXG*JMXG*KMXG !
4942 ! !
4943 ! usage: call setgocartaer !
4944 ! !
4945 ! subprograms called: map_aermr, aeropt_grt !
4946 ! !
4947 ! ================================================================== !
4948 !
4949  implicit none
4950 
4951 ! --- inputs:
4952  integer, intent(in) :: IMAX,NLAY,NLP1,ivflip,NSWLWBD
4953  logical, intent(in) :: lsswr, lslwr
4954 
4955  real (kind=kind_phys), dimension(:,:), intent(in) :: prslk, &
4956  & prsl, rhlay, tvly, dz, hz
4957  real (kind=kind_phys), dimension(:), intent(in) :: alon, alat
4958  real (kind=kind_phys), dimension(:,:,:), intent(in) :: trcly
4959 
4960 ! --- outputs:
4961  real (kind=kind_phys), dimension(:,:,:,:), intent(out) :: &
4962  & aerosw, aerolw
4963 
4964 ! --- locals:
4965  real (kind=kind_phys), dimension(NLAY) :: rh1, dz1
4966  real (kind=kind_phys), dimension(NLAY,NSWLWBD)::tauae,ssaae,asyae
4967  real (kind=kind_phys), dimension(NLAY,max_num_gridcomp) :: &
4968  & tauae_gocart
4969 
4970  real (kind=kind_phys) :: tmp1, tmp2
4971 
4972  integer :: i, i1, i2, j1, j2, k, m, m1, kp
4973 
4974 ! prognostic aerosols on gfs grids
4975  real (kind=kind_phys), dimension(:,:,:),allocatable:: aermr,dmfcs
4976 
4977 ! aerosol (dry mass) on gfs grids/levels
4978  real (kind=kind_phys), dimension(:,:), allocatable :: &
4979  & dmanl,dmclm, dmclmx
4980  real (kind=kind_phys), dimension(KMXG) :: pstmp, pkstr
4981  real (kind=kind_phys) :: ptop, psfc, tem, plv, tv, rho
4982 
4983 ! --- conversion constants
4984  real (kind=kind_phys), parameter :: hdltx = 0.5 * dltx
4985  real (kind=kind_phys), parameter :: hdlty = 0.5 * dlty
4986 
4987 !===> ... begin here
4988 !
4989  if ( .not. allocated(dmanl) ) then
4990  allocate ( dmclmx(kmxg,nmxg) )
4991  allocate ( dmanl(nlay,nmxg) )
4992  allocate ( dmclm(nlay,nmxg) )
4993 
4994  allocate ( aermr(imax,nlay,nmxg) )
4995  allocate ( dmfcs(imax,nlay,nmxg) )
4996  endif
4997 !
5000  dmfcs(:,:,:) = f_zero
5001  lab_if_fcst : if ( get_fcst ) then
5002 
5003  call map_aermr
5004 ! --- inputs: (in scope variables)
5005 ! --- outputs: (in scope variables)
5006 
5007  endif lab_if_fcst
5008 !
5010  lab_do_imax : do i = 1, imax
5011 
5012  dmclm(:,:) = f_zero
5013 
5014  lab_if_clim : if ( get_clim ) then
5015 ! --- map grid in longitude direction
5016  i2 = 1
5017  j2 = 1
5018  tmp1 = alon(i)
5019  if (tmp1 > 180.) tmp1 = tmp1 - 360.0
5020  lab_do_imxg : do i1 = 1, imxg
5021  tmp2 = geos_rlon(i1)
5022  if (tmp2 > 180.) tmp2 = tmp2 - 360.0
5023  if (abs(tmp1-tmp2) <= hdltx) then
5024  i2 = i1
5025  exit lab_do_imxg
5026  endif
5027  enddo lab_do_imxg
5028 
5029 ! --- map grid in latitude direction
5030  lab_do_jmxg : do j1 = 1, jmxg
5031  if (abs(alat(i)-geos_rlat(j1)) <= hdlty) then
5032  j2 = j1
5033  exit lab_do_jmxg
5034  endif
5035  enddo lab_do_jmxg
5036 
5037 ! --- update local arrays pstmp and dmclmx
5038  pstmp(:)= psclmg(i2,j2,:)*1000.0 ! cb to Pa
5039  dmclmx(:,:) = dmclmg(i2,j2,:,:)
5040 
5041 ! --- map geos-gocart climo (dmclmx) to gfs level (dmclm)
5042  pkstr(:)=fpkap(pstmp(:))
5043  psfc = pkstr(1) ! pressure at sfc
5044  ptop = pkstr(kmxg) ! pressure at toa
5045 
5046 ! --- map grid in verical direction (follow how ozone is mapped
5047 ! in radiation_gases routine)
5048  do k = 1, nlay
5049  kp = k ! from sfc to toa
5050  if(ivflip==0) kp = nlay - k + 1 ! from toa to sfc
5051  tmp1 = prslk(i,kp)
5052 
5053  do m1 = 1, kmxg - 1 ! from sfc to toa
5054  if(tmp1 > pkstr(m1+1) .and. tmp1 <= pkstr(m1)) then
5055  tmp2 = f_one / (pkstr(m1)-pkstr(m1+1))
5056  tem = (pkstr(m1) - tmp1) * tmp2
5057  dmclm(kp,:) = tem * dmclmx(m1+1,:)+ &
5058  & (f_one-tem) * dmclmx(m1,:)
5059  endif
5060  enddo
5061 
5062 !* if(tmp1 > psfc) dmclm(kp,:) = dmclmx(1,:)
5063 !* if(tmp1 < ptop) dmclm(kp,:) = dmclmx(KMXG,:)
5064 
5065  enddo
5066  endif lab_if_clim
5067 !
5068 ! --- compute fcst/clim merged aerosol loading (dmanl) and the
5069 ! radiation optical properties (aerosw, aerolw)
5070 !
5071  do k = 1, nlay
5072 
5073 ! --- map global to local arrays (rh1 and dz1)
5074  rh1(k) = rhlay(i,k)
5075  dz1(k) = dz (i,k)
5076 
5077 ! --- convert from mixing ratio to dry mass (g/m3)
5078  plv = 100. * prsl(i,k) ! convert pressure from mb to Pa
5079  tv = tvly(i,k) ! virtual temp in K
5080  rho = plv / (con_rd * tv) ! air density in kg/m3
5081  if ( get_fcst ) then
5082  do m = 1, nmxg ! mixing ratio (g/g)
5083  dmfcs(i,k,m) = max(1000.*(rho*aermr(i,k,m)),f_zero)
5084  enddo ! m_do_loop
5085  endif
5086  if ( get_clim .and. (gocart_climo == 'ver4') ) then
5087  do m = 1, nmxg
5088  dmclm(k,m)=1000.*dmclm(k,m)*rho !mixing ratio (g/g)
5089  if ( molwgt(m) /= 0. ) then !mixing ratio (mol/mol)
5090  dmclm(k,m)=dmclm(k,m) * (molwgt(m)/con_amd)
5091  endif
5092  enddo ! m_do_loop
5093  endif
5094 
5095 
5096 ! --- determine dmanl from dmclm and dmfcs
5097  do m = 1, nmxg
5098  dmanl(k,m)= ctaer*dmfcs(i,k,m) + &
5099  & ( f_one-ctaer)*dmclm(k,m)
5100  enddo
5101  enddo
5102 
5105 
5106  call aeropt_grt
5107 ! --- inputs: (in scope variables)
5108 ! --- outputs: (in scope variables)
5109 
5110  if ( lsswr ) then
5111 
5112  if ( laswflg ) then
5113 
5114  do m = 1, nbdsw
5115  do k = 1, nlay
5116  aerosw(i,k,m,1) = tauae(k,m)
5117  aerosw(i,k,m,2) = ssaae(k,m)
5118  aerosw(i,k,m,3) = asyae(k,m)
5119  enddo
5120  enddo
5121 
5122  else
5123 
5124  aerosw(:,:,:,:) = f_zero
5125 
5126  endif
5127 
5128  endif ! end if_lsswr_block
5129 
5130  if ( lslwr ) then
5131 
5132  if ( lalwflg ) then
5133 
5134  if ( nlwbnd == 1 ) then
5135  m1 = nbdsw + 1
5136  do m = 1, nbdlw
5137  do k = 1, nlay
5138  aerolw(i,k,m,1) = tauae(k,m1)
5139  aerolw(i,k,m,2) = ssaae(k,m1)
5140  aerolw(i,k,m,3) = asyae(k,m1)
5141  enddo
5142  enddo
5143  else
5144  do m = 1, nbdlw
5145  m1 = nbdsw + m
5146  do k = 1, nlay
5147  aerolw(i,k,m,1) = tauae(k,m1)
5148  aerolw(i,k,m,2) = ssaae(k,m1)
5149  aerolw(i,k,m,3) = asyae(k,m1)
5150  enddo
5151  enddo
5152  endif
5153 
5154  else
5155 
5156  aerolw(:,:,:,:) = f_zero
5157 
5158  endif
5159  endif ! end if_lslwr_block
5160 
5161  enddo lab_do_imax
5162 
5163 ! =================
5164  contains
5165 ! =================
5166 
5170 !-----------------------------
5171  subroutine map_aermr
5172 !.............................
5173 ! --- inputs: (in scope variables)
5174 ! --- outputs: (in scope variables)
5175 
5176 ! ==================================================================== !
5177 ! !
5178 ! subprogram: map_aermr !
5179 ! !
5180 ! map input tracer fields (trcly) to local tracer array (aermr) !
5181 ! !
5182 ! ==================== defination of variables =================== !
5183 ! !
5184 ! input arguments: !
5185 ! IMAX - horizontal dimension of arrays 1 !
5186 ! NLAY - vertical dimensions of arrays 1 !
5187 ! trcly - layer tracer mass mixing ratio g/g IMAX*NLAY*NTRAC!
5188 ! output arguments: (to module variables) !
5189 ! aermr - layer aerosol mass mixing ratio g/g IMAX*NLAY*NMXG !
5190 ! !
5191 ! note: !
5192 ! NTRAC is the number of tracers excluding water vapor !
5193 ! NMXG is the number of prognostic aerosol species !
5194 ! ================================================================== !
5195 !
5196  implicit none
5197 
5198 ! --- inputs:
5199 ! --- output:
5200 
5201 ! --- local:
5202  integer :: i, indx, ii
5203  character :: tp*2
5204 
5205 ! initialize
5206  aermr(:,:,:) = f_zero
5207  ii = 1 !! <---- trcly does not contain q
5208 
5209 ! ==> DU: du1 (submicron bins), du2, du3, du4, du5
5210  if( gfs_phy_tracer%doing_DU ) then
5211  aermr(:,:,dm_indx%dust1) = trcly(:,:,dmfcs_indx%du001-ii)
5212  aermr(:,:,dm_indx%dust2) = trcly(:,:,dmfcs_indx%du002-ii)
5213  aermr(:,:,dm_indx%dust3) = trcly(:,:,dmfcs_indx%du003-ii)
5214  aermr(:,:,dm_indx%dust4) = trcly(:,:,dmfcs_indx%du004-ii)
5215  aermr(:,:,dm_indx%dust5) = trcly(:,:,dmfcs_indx%du005-ii)
5216  endif
5217 
5218 ! ==> OC: oc_phobic, oc_philic
5219  if( gfs_phy_tracer%doing_OC ) then
5220  aermr(:,:,dm_indx%waso_phobic) = &
5221  & trcly(:,:,dmfcs_indx%ocphobic-ii)
5222  aermr(:,:,dm_indx%waso_philic) = &
5223  & trcly(:,:,dmfcs_indx%ocphilic-ii)
5224  endif
5225 
5226 ! ==> BC: bc_phobic, bc_philic
5227  if( gfs_phy_tracer%doing_BC ) then
5228  aermr(:,:,dm_indx%soot_phobic) = &
5229  & trcly(:,:,dmfcs_indx%bcphobic-ii)
5230  aermr(:,:,dm_indx%soot_philic) = &
5231  & trcly(:,:,dmfcs_indx%bcphilic-ii)
5232  endif
5233 
5234 ! ==> SS: ss1, ss2 (submicron bins), ss3, ss4, ss5
5235  if( gfs_phy_tracer%doing_SS ) then
5236  aermr(:,:,dm_indx%ssam) = trcly(:,:,dmfcs_indx%ss001-ii) &
5237  & + trcly(:,:,dmfcs_indx%ss002-ii)
5238  aermr(:,:,dm_indx%sscm) = trcly(:,:,dmfcs_indx%ss003-ii) &
5239  & + trcly(:,:,dmfcs_indx%ss004-ii) &
5240  & + trcly(:,:,dmfcs_indx%ss005-ii)
5241  endif
5242 
5243 ! ==> SU: so4
5244  if( gfs_phy_tracer%doing_SU ) then
5245  aermr(:,:,dm_indx%suso) = trcly(:,:,dmfcs_indx%so4-ii)
5246  endif
5247 
5248  return
5249 !...................................
5250  end subroutine map_aermr
5251 !-----------------------------------
5252 
5253 
5258 !-----------------------------------
5259  subroutine aeropt_grt
5260 !...................................
5261 ! --- inputs: (in scope variables)
5262 ! --- outputs: (in scope variables)
5263 
5264 ! ================================================================== !
5265 ! !
5266 ! subprogram: aeropt_grt !
5267 ! !
5268 ! compute aerosols optical properties in NSWLWBD sw/lw bands. !
5269 ! Aerosol distribution at each grid point is composed from up to !
5270 ! NMXG aerosol species (from NUM_GRIDCOMP components). !
5271 ! !
5272 ! input variables: !
5273 ! dmanl - aerosol dry mass g/m3 NLAY*NMXG !
5274 ! rh1 - relative humidity % NLAY !
5275 ! dz1 - layer thickness km NLAY !
5276 ! NLAY - vertical dimensions - 1 !
5277 ! ivflip - control flag for direction of vertical index !
5278 ! =0: index from toa to surface !
5279 ! =1: index from surface to toa !
5280 ! !
5281 ! output variables: !
5282 ! tauae - aerosol optical depth - NLAY*NSWLWBD !
5283 ! ssaae - aerosol single scattering albedo - NLAY*NSWLWBD !
5284 ! asyae - aerosol asymmetry parameter - NLAY*NSWLWBD !
5285 ! !
5286 ! ================================================================== !
5287 !
5288  implicit none
5289 
5290 ! --- inputs:
5291 ! --- outputs:
5292 
5293 ! --- locals:
5294  real (kind=kind_phys) :: aerdm
5295  real (kind=kind_phys) :: ext1, ssa1, asy1, ex00, ss00, as00, &
5296  & ex01, ss01, as01, exint
5297  real (kind=kind_phys) :: tau, ssa, asy, &
5298  & sum_tau, sum_ssa, sum_asy
5299 
5300 ! --- subgroups for sub-micron dust
5301 ! --- corresponds to 0.1-0.18, 0.18-0.3, 0.3-0.6, 0.6-1.0 micron
5302 
5303  real (kind=kind_phys) :: fd(4)
5304  data fd / 0.01053,0.08421,0.25263,0.65263 /
5305 
5306  character :: tp*2
5307  integer :: icmp, n, kk, ib, ih2, ih1, ii, ij, ijk
5308  real (kind=kind_phys) :: drh0, drh1, rdrh
5309 
5310  real (kind=kind_phys) :: qmin
5311  data qmin / 1.e-20 /
5312 
5313 !===> ... begin here
5314 
5315 ! --- initialize (assume no aerosols)
5316  tauae = f_zero
5317  ssaae = f_one
5318  asyae = f_zero
5319 
5320  tauae_gocart = f_zero
5321 
5322 !===> ... loop over vertical layers
5323 !
5324  lab_do_layer : do kk = 1, nlay
5325 
5326 ! --- linear interp coeffs for rh-dep species
5327 
5328  ih2 = 1
5329  do while ( rh1(kk) > rhlev_grt(ih2) )
5330  ih2 = ih2 + 1
5331  if ( ih2 > krhlev ) exit
5332  enddo
5333  ih1 = max( 1, ih2-1 )
5334  ih2 = min( krhlev, ih2 )
5335 
5336  drh0 = rhlev_grt(ih2) - rhlev_grt(ih1)
5337  drh1 = rh1(kk) - rhlev_grt(ih1)
5338  if ( ih1 == ih2 ) then
5339  rdrh = f_zero
5340  else
5341  rdrh = drh1 / drh0
5342  endif
5343 
5344 ! --- loop through sw/lw spectral bands
5345 
5346  lab_do_ib : do ib = 1, nswlwbd
5347  sum_tau = f_zero
5348  sum_ssa = f_zero
5349  sum_asy = f_zero
5350 
5351 ! --- loop through aerosol grid components
5352  lab_do_icmp : do icmp = 1, num_gridcomp
5353  ext1 = f_zero
5354  ssa1 = f_zero
5355  asy1 = f_zero
5356 
5357  tp = gridcomp(icmp)
5358 
5359  select case ( tp )
5360 
5361 ! -- dust aerosols: no humidification effect
5362  case ( 'DU')
5363  do n = 1, kcm1
5364 
5365  if (n <= 4) then
5366  aerdm = dmanl(kk,dm_indx%dust1) * fd(n)
5367  else
5368  aerdm = dmanl(kk,dm_indx%dust1+n-4 )
5369  endif
5370 
5371  if (aerdm < qmin) aerdm = f_zero
5372  ex00 = extrhi_grt(n,ib)*(1000.*dz1(kk))*aerdm
5373  ss00 = ssarhi_grt(n,ib)
5374  as00 = asyrhi_grt(n,ib)
5375  ext1 = ext1 + ex00
5376  ssa1 = ssa1 + ex00 * ss00
5377  asy1 = asy1 + ex00 * ss00 * as00
5378 
5379  enddo
5380 
5381 ! -- suso aerosols: with humidification effect
5382  case ( 'SU')
5383  ij = isuso
5384  exint = extrhd_grt(ih1,ij,ib) &
5385  & + rdrh*(extrhd_grt(ih2,ij,ib) - extrhd_grt(ih1,ij,ib))
5386  ss00 = ssarhd_grt(ih1,ij,ib) &
5387  & + rdrh*(ssarhd_grt(ih2,ij,ib) - ssarhd_grt(ih1,ij,ib))
5388  as00 = asyrhd_grt(ih1,ij,ib) &
5389  & + rdrh*(asyrhd_grt(ih2,ij,ib) - asyrhd_grt(ih1,ij,ib))
5390 
5391  aerdm = dmanl(kk, dm_indx%suso)
5392  if (aerdm < qmin) aerdm = f_zero
5393  ex00 = exint*(1000.*dz1(kk))*aerdm
5394  ext1 = ex00
5395  ssa1 = ex00 * ss00
5396  asy1 = ex00 * ss00 * as00
5397 
5398 ! -- seasalt aerosols: with humidification effect
5399  case ( 'SS')
5400  do n = 1, 2
5401  ij = issam + (n-1)
5402  exint = extrhd_grt(ih1,ij,ib) &
5403  & + rdrh*(extrhd_grt(ih2,ij,ib) - extrhd_grt(ih1,ij,ib))
5404  ss00 = ssarhd_grt(ih1,ij,ib) &
5405  & + rdrh*(ssarhd_grt(ih2,ij,ib) - ssarhd_grt(ih1,ij,ib))
5406  as00 = asyrhd_grt(ih1,ij,ib) &
5407  & + rdrh*(asyrhd_grt(ih2,ij,ib) - asyrhd_grt(ih1,ij,ib))
5408 
5409  aerdm = dmanl(kk, dm_indx%ssam+n-1)
5410  if (aerdm < qmin) aerdm = f_zero
5411  ex00 = exint*(1000.*dz1(kk))*aerdm
5412  ext1 = ext1 + ex00
5413  ssa1 = ssa1 + ex00 * ss00
5414  asy1 = asy1 + ex00 * ss00 * as00
5415 
5416  enddo
5417 
5418 ! -- organic carbon/black carbon:
5419 ! using 'waso' and 'soot' for hydrophilic OC and BC
5420 ! using 'waso' and 'soot' at RH=0 for hydrophobic OC and BC
5421  case ( 'OC', 'BC')
5422  if(tp == 'OC') then
5423  ii = dm_indx%waso_phobic
5424  ij = iwaso
5425  else
5426  ii = dm_indx%soot_phobic
5427  ij = isoot
5428  endif
5429 
5430 ! --- hydrophobic
5431  aerdm = dmanl(kk, ii)
5432  if (aerdm < qmin) aerdm = f_zero
5433  ex00 = extrhd_grt(1,ij,ib)*(1000.*dz1(kk))*aerdm
5434  ss00 = ssarhd_grt(1,ij,ib)
5435  as00 = asyrhd_grt(1,ij,ib)
5436 ! --- hydrophilic
5437  aerdm = dmanl(kk, ii+1)
5438  if (aerdm < qmin) aerdm = f_zero
5439  exint = extrhd_grt(ih1,ij,ib) &
5440  & + rdrh*(extrhd_grt(ih2,ij,ib) - extrhd_grt(ih1,ij,ib))
5441  ex01 = exint*(1000.*dz1(kk))*aerdm
5442  ss01 = ssarhd_grt(ih1,ij,ib) &
5443  & + rdrh*(ssarhd_grt(ih2,ij,ib) - ssarhd_grt(ih1,ij,ib))
5444  as01 = asyrhd_grt(ih1,ij,ib) &
5445  & + rdrh*(asyrhd_grt(ih2,ij,ib) - asyrhd_grt(ih1,ij,ib))
5446 
5447  ext1 = ex00 + ex01
5448  ssa1 = (ex00 * ss00) + (ex01 * ss01)
5449  asy1 = (ex00 * ss00 * as00) + (ex01 * ss01 * as01)
5450 
5451  end select
5452 
5453 ! --- determine tau, ssa, asy for each grid component
5454  tau = ext1
5455  if (ext1 > f_zero) ssa=min(f_one,ssa1/ext1)
5456  if (ssa1 > f_zero) asy=min(f_one,asy1/ssa1)
5457 
5458 ! --- save tau at 550 nm for each grid component
5459  if ( ib == nv_aod ) then
5460  do ijk = 1, max_num_gridcomp
5461  if ( tp == max_gridcomp(ijk) ) then
5462  tauae_gocart(kk,ijk) = tau
5463  endif
5464  enddo
5465  endif
5466 
5467 ! --- update sum_tau, sum_ssa, sum_asy
5468  sum_tau = sum_tau + tau
5469  sum_ssa = sum_ssa + tau * ssa
5470  sum_asy = sum_asy + tau * ssa * asy
5471 
5472  enddo lab_do_icmp
5473 
5474 
5475 ! --- determine total tau, ssa, asy for aerosol mixture
5476  tauae(kk,ib) = sum_tau
5477  if (sum_tau > f_zero) ssaae(kk,ib) = sum_ssa / sum_tau
5478  if (sum_ssa > f_zero) asyae(kk,ib) = sum_asy / sum_ssa
5479 
5480  enddo lab_do_ib
5481 
5482  enddo lab_do_layer
5483 
5484 !
5485  return
5486 !...................................
5487  end subroutine aeropt_grt
5488 !--------------------------------
5489 
5490 !................................
5491  end subroutine setgocartaer
5492 !--------------------------------
5493 !! @}
5494 !
5495 ! GOCART code modification end here (Sarah Lu) ------------------------!
5496 ! =======================================================================
5497 
5498 !..........................................!
5499  end module module_radiation_aerosols !
5500 !==========================================!
module_radiation_aerosols::ndm
integer, parameter ndm
num of atmos aerosols domains
Definition: radiation_aerosols.f:296
module_radiation_aerosols::nspc
integer, parameter, public nspc
num of species for output aod (opnl)
Definition: radiation_aerosols.f:177
module_radiation_aerosols::kcm2
integer, save kcm2
num of rh dep aerosols (set in subr set_aerspc)
Definition: radiation_aerosols.f:392
physparam::iaermdl
integer, save iaermdl
aerosol model scheme control flag
Definition: physparam.f:146
radclimaer
subroutine radclimaer
This subroutine computes aerosols optical properties in NSWLWBD bands. there are seven different vert...
Definition: radiation_aerosols.f:3198
physcons::con_pi
real(kind=kind_phys), parameter con_pi
pi
Definition: physcons.f:48
module_radiation_aerosols::nlwstr
integer, parameter, public nlwstr
starting band number in ir region
Definition: radiation_aerosols.f:175
module_radiation_aerosols::lckprnt
logical, save lckprnt
logical parameter for gocart debug print control
Definition: radiation_aerosols.f:561
physcons::con_g
real(kind=kind_phys), parameter con_g
gravity ( )
Definition: physcons.f:59
module_radiation_aerosols::nswlwbd
integer, save nswlwbd
total number of bands for sw+lw aerosols
Definition: radiation_aerosols.f:192
physparam::lalwflg
logical, save lalwflg
LW aerosols effect control flag.
Definition: physparam.f:150
module_radiation_aerosols::kmxg
integer, parameter kmxg
num of vertical layers in geos dataset
Definition: radiation_aerosols.f:517
module_radiation_aerosols::kyrend
integer kyrend
ending year of data in the input file
Definition: radiation_aerosols.f:279
module_radiation_aerosols::get_clim
logical, save get_clim
option to get clim gocart aerosol field
Definition: radiation_aerosols.f:571
module_radiation_aerosols::nxc
integer, parameter nxc
num of max componets in a profile
Definition: radiation_aerosols.f:292
physparam::iaerflg
integer, save iaerflg
aerosol effect control flag
Definition: physparam.f:148
module_radiation_aerosols::kyrsav
integer kyrsav
the year of data in use in the input file
Definition: radiation_aerosols.f:281
module_radiation_aerosols::geos_rlat
real(kind=kind_phys), dimension(:), allocatable, save geos_rlat
geos-gocart latitude arrays
Definition: radiation_aerosols.f:540
module_radiation_aerosols::ncm2
integer, parameter ncm2
num of rh dependent aeros species
Definition: radiation_aerosols.f:308
module_radiation_aerosols::rhdpext0_grt
real(kind=kind_phys), dimension(:,:,:), allocatable rhdpext0_grt
extinction coefficient
Definition: radiation_aerosols.f:455
module_radiation_aerosols::krhlev
integer, parameter krhlev
num of rh levels for rh-dep components
Definition: radiation_aerosols.f:385
module_radiation_aerosols::isoot
integer, save isoot
index for rh dependent aerosol optical properties (2nd dimension for extrhd_grt, ssarhd_grt, and asyrhd_grt)
Definition: radiation_aerosols.f:593
optavg
subroutine optavg
This subroutine computes mean aerosols optical properties over each SW radiation spectral band for ea...
Definition: radiation_aerosols.f:1537
aeropt_grt
subroutine aeropt_grt
This subroutine computes aerosols optical properties in NSWLWBD SW/LW bands. Aerosol distribution at ...
Definition: radiation_aerosols.f:5260
module_radiation_aerosols::naerbnd
integer, parameter naerbnd
num of bands for clim aer data (opac)
Definition: radiation_aerosols.f:302
module_radiation_aerosols::iendwv_grt
integer, dimension(:), allocatable iendwv_grt
spectral band structure: ending wavenumber ( ) for each band
Definition: radiation_aerosols.f:434
module_radiation_aerosols::kcm1
integer, save kcm1
num of rh indep aerosols (set in subr set_aerspc)
Definition: radiation_aerosols.f:390
module_radiation_aerosols::psclmg
real(kind=kind_phys), dimension(:,:,:), allocatable, save psclmg
pressure in cb
Definition: radiation_aerosols.f:532
module_radiation_aerosols::geos_rlon
real(kind=kind_phys), dimension(:), allocatable, save geos_rlon
geos-gocart longitude arrays
Definition: radiation_aerosols.f:538
module_radiation_aerosols::gocart_climo
character *4, save gocart_climo
control flag for gocart climo data set: xxxx as default; ver3 for geos3; ver4 for geos4; 0000 for unk...
Definition: radiation_aerosols.f:544
module_radiation_aerosols::rhdpssa0_grt
real(kind=kind_phys), dimension(:,:,:), allocatable rhdpssa0_grt
single scattering albedo
Definition: radiation_aerosols.f:457
module_radiation_aerosols::nf_aelw
integer, parameter, public nf_aelw
num of output fields for LW rad
Definition: radiation_aerosols.f:173
module_radiation_aerosols::jmxae
integer, parameter jmxae
num of lat-points in glb aeros data set
Definition: radiation_aerosols.f:300
module_radiation_aerosols::idxcg
integer, dimension(nxc, imxae, jmxae), save idxcg
aeros component index
Definition: radiation_aerosols.f:370
module_radiation_aerosols::nwvsol
integer, parameter, public nwvsol
num of wvnum regions where solar flux is constant
Definition: radiation_aerosols.f:201
module_radiation_aerosols::num_gridcomp
integer, save num_gridcomp
number of aerosol grid components
Definition: radiation_aerosols.f:634
module_radiation_aerosols::get_fcst
logical, save get_fcst
option to get fcst gocart aerosol field
Definition: radiation_aerosols.f:569
physcons::con_t0c
real(kind=kind_phys), parameter con_t0c
temp at 0C (K)
Definition: physcons.f:96
module_radiation_aerosols::kprfg
integer, dimension( imxae, jmxae), save kprfg
aeros profile index
Definition: radiation_aerosols.f:372
module_radiation_aerosols::nwvtot
integer, parameter, public nwvtot
total num of wvnum included
Definition: radiation_aerosols.f:204
module_radiation_aerosols::nwvns0
integer, dimension(nwvsol), save nwvns0
number of wavenumbers in each region where the solar flux is constant
Definition: radiation_aerosols.f:209
module_radiation_aerosols::max_num_gridcomp
integer, parameter max_num_gridcomp
default full-package setting
Definition: radiation_aerosols.f:639
module_radiation_aerosols::prsref
real(kind=kind_phys), dimension(ndm, nae), save prsref
ref pressure lev (sfc to toa) in mb (100Pa)
Definition: radiation_aerosols.f:325
module_radiation_aerosols::setgocartaer
subroutine setgocartaer
This subroutine computes SW + LW aerosol optical properties for gocart aerosol species (merged from f...
Definition: radiation_aerosols.f:4893
module_radiation_aerosols::nrhlev
integer, parameter nrhlev
num of rh levels for rh-dep components
Definition: radiation_aerosols.f:304
physparam::lalw1bd
logical, parameter lalw1bd
=t: use 1 broad-band LW aeros properties =f: use multi bands aeros properites
Definition: physparam.f:121
rd_gocart_clim
subroutine rd_gocart_clim
This subroutine:
Definition: radiation_aerosols.f:4543
module_radiation_aerosols::ivolae
integer, dimension(:,:,:), allocatable, save ivolae
monthly, 45-deg lat-zone aerosols data set in subroutine 'aer_init'
Definition: radiation_aerosols.f:273
module_radiation_aerosols::maxvyr
integer, parameter maxvyr
upper limit (year) data available
Definition: radiation_aerosols.f:270
module_radiation_aerosols::gocart_index_type
index for gocart aerosol species to be included in the calculations of aerosol optical properties (ex...
Definition: radiation_aerosols.f:603
physcons::con_c
real(kind=kind_phys), parameter con_c
speed of light ( )
Definition: physcons.f:123
physparam::aeros_file
character, save aeros_file
external aerosols data file: aerosol.dat
Definition: physparam.f:156
module_radiation_aerosols::max_gridcomp
character *2, dimension(max_num_gridcomp) max_gridcomp
data max_gridcomp /'DU', 'BC', 'OC', 'SU', 'SS'/
Definition: radiation_aerosols.f:641
module_radsw_parameters
This module contains SW band parameters set up.
Definition: radsw_param.f:65
module_radiation_aerosols::rhdpasy0_grt
real(kind=kind_phys), dimension(:,:,:), allocatable rhdpasy0_grt
asymmetry parameter
Definition: radiation_aerosols.f:459
physcons::con_plnk
real(kind=kind_phys), parameter con_plnk
planck constant ( )
Definition: physcons.f:125
module_radiation_aerosols::kaerbnd
integer, parameter kaerbnd
num of bands for aer data (gocart)
Definition: radiation_aerosols.f:383
module_radiation_aerosols::lgrtint
logical, save lgrtint
logical parameter for gocart initialization control
Definition: radiation_aerosols.f:557
module_radiation_aerosols::gocart_init
subroutine gocart_init
The initialization program for gocart aerosols.
Definition: radiation_aerosols.f:3530
physparam::laswflg
logical, save laswflg
SW aerosols effect control flag.
Definition: physparam.f:152
optavg_grt
subroutine optavg_grt
This subroutine computes mean aerosols optical properties over each SW/LW radiation spectral band for...
Definition: radiation_aerosols.f:4324
trop_update
subroutine trop_update
This subroutine updates the monthly global distribution of aerosol profiles in five degree horizontal...
Definition: radiation_aerosols.f:1866
set_volcaer
subroutine set_volcaer
The initialization program for stratospheric volcanic aerosols.
Definition: radiation_aerosols.f:1026
module_radiation_aerosols::wvn550
real(kind=kind_phys), parameter wvn550
the wavenumber ( ) of wavelength 550nm for diagnostic aod output
Definition: radiation_aerosols.f:664
module_radiation_aerosols::s0intv
real(kind=kind_phys), dimension(nwvsol), save s0intv
solar flux in each wvnumb region where it is constant
Definition: radiation_aerosols.f:225
map_aermr
subroutine map_aermr
This subroutine maps input tracer fields (trcly) to local tracer array (aermr).
Definition: radiation_aerosols.f:5172
module_radiation_aerosols::ncm1
integer, parameter ncm1
num of rh independent aeros species
Definition: radiation_aerosols.f:306
physcons::con_rd
real(kind=kind_phys), parameter con_rd
gas constant air ( )
Definition: physcons.f:76
volc_update
subroutine volc_update
This subroutine searches historical volcanic data sets to find and read in monthly 45-degree lat-zone...
Definition: radiation_aerosols.f:2027
module_radlw_parameters
This module contains LW band parameters set up.
Definition: radlw_param.f:64
module_radiation_aerosols::aer_init
subroutine, public aer_init
The initialization program is to set up necessary parameters and working arrays.
Definition: radiation_aerosols.f:686
module_radiation_aerosols::nwvtir
integer, parameter, public nwvtir
total num of wvnum in ir range
Definition: radiation_aerosols.f:206
module_radiation_aerosols::imxae
integer, parameter imxae
num of lon-points in glb aeros data set
Definition: radiation_aerosols.f:298
wrt_aerlog
subroutine wrt_aerlog
This subroutine writes aerosol parameter configuration to run log file.
Definition: radiation_aerosols.f:834
set_aercoef
subroutine set_aercoef
The initialization program for climatological aerosols. The program reads and maps the pre-tabulated ...
Definition: radiation_aerosols.f:1172
module_radiation_aerosols::rhidasy0_grt
real(kind=kind_phys), dimension(:,:), allocatable rhidasy0_grt
asymmetry parameter
Definition: radiation_aerosols.f:446
module_radiation_aerosols::kmonsav
integer kmonsav
the month of data in use in the input file
Definition: radiation_aerosols.f:283
module_radiation_aerosols::setaer
subroutine, public setaer
This subroutine computes aerosols optical properties.
Definition: radiation_aerosols.f:2181
module_radiation_aerosols::nswbnd
integer, save nswbnd
number of actual bands for sw aerosols; calculated according to laswflg setting
Definition: radiation_aerosols.f:187
module_radiation_aerosols::rhidssa0_grt
real(kind=kind_phys), dimension(:,:), allocatable rhidssa0_grt
single scattering albedo
Definition: radiation_aerosols.f:444
module_radiation_aerosols::kyrstr
integer kyrstr
starting year of data in the input file
Definition: radiation_aerosols.f:277
module_radiation_aerosols::nv_aod
integer, save nv_aod
the sw spectral band covering wvn550 (comp in aer_init)
Definition: radiation_aerosols.f:666
module_radiation_aerosols::clim_aerinit
subroutine clim_aerinit
This subroutine is the opac-climatology aerosol initialization program to set up necessary parameters...
Definition: radiation_aerosols.f:1081
set_spectrum
subroutine set_spectrum
This subroutine defines the one wavenumber solar fluxes based on toa solar spectral distribution...
Definition: radiation_aerosols.f:933
module_radiation_aerosols::imxg
integer, parameter imxg
num of lon-points in geos dataset
Definition: radiation_aerosols.f:513
module_radiation_aerosols::nmxg
integer, save nmxg
to be determined by set_aerspc
Definition: radiation_aerosols.f:520
module_radiation_aerosols::nspc1
integer, parameter, public nspc1
total+species
Definition: radiation_aerosols.f:179
module_radiation_aerosols::sigref
real(kind=kind_phys), dimension(ndm, nae), save sigref
ref sigma lev (sfc to toa)
Definition: radiation_aerosols.f:327
module_radiation_aerosols::aer_property
subroutine aer_property
This subroutine maps the 5 degree global climatological aerosol data set onto model grids...
Definition: radiation_aerosols.f:2769
module_radiation_aerosols::ctaer
real(kind=kind_phys), save ctaer
merging coefficients for fcst/clim; determined from fdaer
Definition: radiation_aerosols.f:566
module_radiation_aerosols::rhlev
real(kind=kind_phys), dimension(nrhlev), save rhlev
predefined relative humidity levels
Definition: radiation_aerosols.f:312
module_radiation_aerosols::cmixg
real(kind=kind_phys), dimension(nxc, imxae, jmxae), save cmixg
aeros component mixing ratio
Definition: radiation_aerosols.f:366
module_radiation_aerosols::nlwbnd
integer, save nlwbnd
number of actual bands for lw aerosols; calculated according to lalwflg and lalw1bd settings ...
Definition: radiation_aerosols.f:190
module_radiation_aerosols::tracer_index_type
index for gocart aerosols from prognostic tracer fields
Definition: radiation_aerosols.f:618
module_radiation_aerosols::kcm
integer, save kcm
=KCM1+KCM2 (set in subr set_aerspc)
Definition: radiation_aerosols.f:394
module_radiation_aerosols::idxspc
integer, dimension(ncm) idxspc
index conversion array:data idxspc / 1, 2, 1, 1, 1, 1, 3, 5, 5, 4 /
Definition: radiation_aerosols.f:657
physparam::ivflip
integer, save ivflip
vertical profile indexing flag
Definition: physparam.f:222
module_radiation_aerosols::jmxg
integer, parameter jmxg
num of lat-points in geos dataset
Definition: radiation_aerosols.f:515
module_radiation_aerosols::aer_update
subroutine, public aer_update
This subroutine checks and updates time varying climatology aerosol data sets.
Definition: radiation_aerosols.f:1805
module_radiation_aerosols::molwgt
real(kind=kind_io4), dimension(:), allocatable molwgt
molecular wght of gocart aerosol species
Definition: radiation_aerosols.f:547
module_radiation_aerosols::denng
real(kind=kind_phys), dimension(2,imxae, jmxae), save denng
aeros number density
Definition: radiation_aerosols.f:368
rd_gocart_luts
subroutine rd_gocart_luts
This subroutine reads input gocart aerosol optical data from Mie code calculations.
Definition: radiation_aerosols.f:4147
module_radiation_aerosols::dmclmg
real(kind=kind_phys), dimension(:,:,:,:), allocatable, save dmclmg
aerosol dry mass in g/m3 or aerosol mixing ration in mol/mol or Kg/Kg
Definition: radiation_aerosols.f:534
physcons::con_boltz
real(kind=kind_phys), parameter con_boltz
boltzmann constant ( )
Definition: physcons.f:127
module_radiation_aerosols::rhidext0_grt
real(kind=kind_phys), dimension(:,:), allocatable rhidext0_grt
extinction coefficient
Definition: radiation_aerosols.f:442
physparam::lavoflg
logical, save lavoflg
stratospheric volcanic effect flag
Definition: physparam.f:154
module_radiation_aerosols::gridcomp
character, dimension(:), allocatable, save gridcomp
aerosol grid components
Definition: radiation_aerosols.f:636
module_radiation_aerosols::nae
integer, parameter nae
num of aerosols profile structures
Definition: radiation_aerosols.f:294
module_radiation_aerosols::minvyr
integer, parameter minvyr
lower limit (year) data available
Definition: radiation_aerosols.f:268
set_aerspc
subroutine set_aerspc(raddt, fdaer)
This subroutine determines merging coefficients ctaer; setup aerosol specification.
Definition: radiation_aerosols.f:3847
physcons::con_amd
real(kind=kind_phys), parameter con_amd
molecular wght of dry air ( )
Definition: physcons.f:136
module_radiation_aerosols::haer
real(kind=kind_phys), dimension(ndm, nae), save haer
scale height of aerosols (km)
Definition: radiation_aerosols.f:323