File: C:\NOAA\NEMS_11731\src\atmos\phys\radiation_aerosols.f

1     !!!!!  ==========================================================  !!!!!
2     !!!!!            'module_radiation_aerosols' description           !!!!!
3     !!!!!  ==========================================================  !!!!!
4     !                                                                      !
5     !   this module contains different atmospheric aerosols schemes for    !
6     !   radiation computations.                                            !
7     !                                                                      !
8     !   in the module, the externally callable subroutines are :           !
9     !                                                                      !
10     !      'aerinit'    -- initialization, input aerosol data, etc.        !
11     !         inputs:                                                      !
12     !           (iyear, imon, IAER, me)                                    !
13     !         outputs:                                                     !
14     !           (none)                                                     !
15     !                                                                      !
16     !      'setaer'     -- mapping aeros profile, compute aeros opticals   !
17     !         inputs:                                                      !
18     !           (xlon,xlat,prsi,prsl,tlay,qlay,rhlay,                      !
19     !            prslk,oz,                                                 !
20     !            IMAX,NLAY,NLP1,iflip,lsswr,lslwr)                         !
21     !         outputs:                                                     !
22     !           (aerosw,aerolw,tau_gocart)                                 !
23     !                                                                      !
24     !                                                                      !
25     !   internal subroutine called:                                        !
26     !       clim_aerinit, setclimaer   - for opac climatological aerosols  !
27     !                                                                      !
28     !       gocart_init,  setgocart    - for gocart aerosols               !
29     !                                                                      !
30     !                                                                      !
31     !   external modules referenced:                                       !
32     !                                                                      !
33     !       'module machine'                 in 'machine.f'                !
34     !       'module physcons'                in 'physcons.f'               !
35     !       'module module_radsw_parameters' in 'radsw_xxxx#_param.f'      !
36     !       'module module_radlw_parameters' in 'radlw_xxxx#_param.f'      !
37     !       'module module_radlw_cntr_para'  in 'radsw_xxxx#_param.f'      !
38     !                                                                      !
39     !   output variable definitions:                                       !
40     !       aerosw(IMAX,NLAY,NBDSW,1) - aerosols optical depth for sw      !
41     !       aerosw(IMAX,NLAY,NBDSW,2) - aerosols single scat albedo for sw !
42     !       aerosw(IMAX,NLAY,NBDSW,3) - aerosols asymmetry parameter for sw!
43     !                                                                      !
44     !       aerolw(IMAX,NLAY,NBDLW,1) - aerosols optical depth for lw      !
45     !       aerolw(IMAX,NLAY,NBDLW,2) - aerosols single scattering albedo  !
46     !       aerolw(IMAX,NLAY,NBDLW,3) - aerosols asymetry parameter        !
47     !                                                                      !
48     !                                                                      !
49     !   program history:                                                   !
50     !     apr     2003  ---  y.-t. hou     created                         !
51     !     nov 04, 2003  ---  y.-t. hou     modified version                !
52     !     apr 15, 2005  ---  y.-t. hou     modified module structure       !
53     !     jul     2006  ---  y.-t. hou     add volcanic forcing            !
54     !     feb     2007  ---  y.-t. hou     add generalized spectral band   !
55     !                   interpolation for sw aerosol optical properties    !
56     !     mar     2007  ---  y.-t. hou     add generalized spectral band   !
57     !                   interpolation for lw aerosol optical properties    !
58     !     aug     2007  ---  y.-t. hou     change clim-aer vert domain     !
59     !                   from pressure reference to sigma reference         !
60     !     sep     2007  ---  y.-t. hou     moving temporary allocatable    !
61     !                   module variable arrays to subroutine dynamically   !
62     !                   allocated arrays (eliminate deallocate calls)      !
63     !     jan     2010  ---  sarah lu      add gocart option               !
64     !     may     2010  ---  sarah lu      add geos4-gocart climo          !
65     !     jul     2010  --   s. moorthi - merged NEMS version with new GFS !
66     !                        version                                       !
67     !                                                                      !
68     !                                                                      !
69     !   references for opac climatological aerosols:                       !
70     !     hou et al. 2002  (ncep office note 441)                          !
71     !     hess et al. 1998 - bams v79 831-844                              !
72     !                                                                      !
73     !   references for gocart interactive aerosols:                        !
74     !     chin et al., 2000 - jgr, v105, 24671-24687                       !
75     !                                                                      !
76     !                                                                      !
77     !   references for stratosperic volcanical aerosols:                   !
78     !     sato et al. 1993 - jgr, v98, d12, 22987-22994                    !
79     !                                                                      !
80     !                                                                      !
81     !!!!!  ==========================================================  !!!!!
82     !!!!!                       end descriptions                       !!!!!
83     !!!!!  ==========================================================  !!!!!
84     
85     
86     
87     !========================================!
88           module module_radiation_aerosols   !
89     !........................................!
90     !
91           use machine,      only : kind_io8, kind_phys, kind_io4
92           use physcons,     only : con_pi, con_rd, con_fvirt, con_g,        &
93          &                         con_t0c, con_c, con_boltz, con_plnk,     &
94          &                         con_amd
95     
96           use module_iounitdef,        only : NIAERCM
97           use module_radsw_parameters, only : NBDSW, NSWSTR, wvnum1, wvnum2
98           use module_radlw_parameters, only : NBDLW, wvnlw1, wvnlw2
99           use module_radlw_cntr_para,  only : iaerlw
100           use funcphys,                only : fpkap
101     !     use namelist_physics_def,    only : fhlwr, fhswr, fdaer
102           use gfs_phy_tracer_config,   only : gfs_phy_tracer, trcindx
103     !
104           implicit   none
105     !
106           private
107     
108     !  ---  general use parameter constants:
109           integer, parameter, public :: NF_AESW = 3     ! num of output fields for sw rad
110           integer, parameter, public :: NF_AELW = 3     ! num of output fields for lw rad
111     
112           real (kind=kind_phys), parameter :: f_zero = 0.0
113           real (kind=kind_phys), parameter :: f_one  = 1.0
114     
115     !  ---  module control parameters set in subroutine "aerinit"
116           integer, save :: iaerflg = 1       ! choice of tropospheric aerosol schemes
117                                              ! =0 no aer; =1 opac; =2 gocart
118           logical, public, save :: lalwflg = .false. ! lw aerosols effect
119                                              ! =t compute lw aerosol optical prop
120           logical, public, save :: laswflg = .false. ! sw aerosols effect
121                                              ! =t compute sw aerosol optical prop
122           integer, save :: NBDIR   = NBDLW   ! num of actual bands for lw aerosols
123                                              ! calculated according to iaerlw setting
124           integer, save :: NBDSWLW = NBDSW+NBDLW
125                                              ! total num of bands for sw+lw aerosols
126           logical, save :: lvolflg = .false. ! volcanic forcing
127                                              ! =t include stratos volcanic forcing
128     
129     ! --------------------------------------------------------------------- !
130     !   section-1 : module variables for spectral band interpolation        !
131     !               similar to gfdl-sw treatment (2000 version)             !
132     ! --------------------------------------------------------------------- !
133     
134     !  ---  parameter constants:
135           integer, parameter, public :: NWVSOL  = 151   ! num of wvnum regions where solar
136                                                         ! flux is constant
137           integer, parameter, public :: NWVTOT  = 57600 ! total num of wvnum included
138           integer, parameter, public :: NWVTIR  = 4000  ! total num of wvnum in ir range
139     
140     !  ---  number of wavenumbers in each region where the solar flux is constant
141           integer, dimension(NWVSOL) :: nwvns0
142     
143           data nwvns0   / 100,  11,  14,  18,  24,  33,  50,  83,  12,  12, &
144          &  13,  15,  15,  17,  18,  20,  21,  24,  26,  30,  32,  37,  42, &
145          &  47,  55,  64,  76,  91, 111, 139, 179, 238, 333,  41,  42,  45, &
146          &  46,  48,  51,  53,  55,  58,  61,  64,  68,  71,  75,  79,  84, &
147          &  89,  95, 101, 107, 115, 123, 133, 142, 154, 167, 181, 197, 217, &
148          & 238, 263, 293, 326, 368, 417, 476, 549, 641, 758, 909, 101, 103, &
149          & 105, 108, 109, 112, 115, 117, 119, 122, 125, 128, 130, 134, 137, &
150          & 140, 143, 147, 151, 154, 158, 163, 166, 171, 175, 181, 185, 190, &
151          & 196, 201, 207, 213, 219, 227, 233, 240, 248, 256, 264, 274, 282, &
152          & 292, 303, 313, 325, 337, 349, 363, 377, 392, 408, 425, 444, 462, &
153          & 483, 505, 529, 554, 580, 610, 641, 675, 711, 751, 793, 841, 891, &
154          & 947,1008,1075,1150,1231,1323,1425,1538,1667,1633,14300 /
155     
156     !  ---  solar flux (w/m**2) in each wvnumb region where it is constant
157           real (kind=kind_phys), dimension(NWVSOL) :: s0intv
158     
159           data  s0intv(  1: 50)       /                                     &
160          &     1.60000E-6, 2.88000E-5, 3.60000E-5, 4.59200E-5, 6.13200E-5,  &
161          &     8.55000E-5, 1.28600E-4, 2.16000E-4, 2.90580E-4, 3.10184E-4,  &
162          &     3.34152E-4, 3.58722E-4, 3.88050E-4, 4.20000E-4, 4.57056E-4,  &
163          &     4.96892E-4, 5.45160E-4, 6.00600E-4, 6.53600E-4, 7.25040E-4,  &
164          &     7.98660E-4, 9.11200E-4, 1.03680E-3, 1.18440E-3, 1.36682E-3,  &
165          &     1.57560E-3, 1.87440E-3, 2.25500E-3, 2.74500E-3, 3.39840E-3,  &
166          &     4.34000E-3, 5.75400E-3, 7.74000E-3, 9.53050E-3, 9.90192E-3,  &
167          &     1.02874E-2, 1.06803E-2, 1.11366E-2, 1.15830E-2, 1.21088E-2,  &
168          &     1.26420E-2, 1.32250E-2, 1.38088E-2, 1.44612E-2, 1.51164E-2,  &
169          &     1.58878E-2, 1.66500E-2, 1.75140E-2, 1.84450E-2, 1.94106E-2 /
170           data  s0intv( 51:100)       /                                     &
171          &     2.04864E-2, 2.17248E-2, 2.30640E-2, 2.44470E-2, 2.59840E-2,  &
172          &     2.75940E-2, 2.94138E-2, 3.13950E-2, 3.34800E-2, 3.57696E-2,  &
173          &     3.84054E-2, 4.13490E-2, 4.46880E-2, 4.82220E-2, 5.22918E-2,  &
174          &     5.70078E-2, 6.19888E-2, 6.54720E-2, 6.69060E-2, 6.81226E-2,  &
175          &     6.97788E-2, 7.12668E-2, 7.27100E-2, 7.31610E-2, 7.33471E-2,  &
176          &     7.34814E-2, 7.34717E-2, 7.35072E-2, 7.34939E-2, 7.35202E-2,  &
177          &     7.33249E-2, 7.31713E-2, 7.35462E-2, 7.36920E-2, 7.23677E-2,  &
178          &     7.25023E-2, 7.24258E-2, 7.20766E-2, 7.18284E-2, 7.32757E-2,  &
179          &     7.31645E-2, 7.33277E-2, 7.36128E-2, 7.33752E-2, 7.28965E-2,  &
180          &     7.24924E-2, 7.23307E-2, 7.21050E-2, 7.12620E-2, 7.10903E-2 /
181           data  s0intv(101:151)       /                        7.12714E-2,  &
182          &     7.08012E-2, 7.03752E-2, 7.00350E-2, 6.98639E-2, 6.90690E-2,  &
183          &     6.87621E-2, 6.52080E-2, 6.65184E-2, 6.60038E-2, 6.47615E-2,  &
184          &     6.44831E-2, 6.37206E-2, 6.24102E-2, 6.18698E-2, 6.06320E-2,  &
185          &     5.83498E-2, 5.67028E-2, 5.51232E-2, 5.48645E-2, 5.12340E-2,  &
186          &     4.85581E-2, 4.85010E-2, 4.79220E-2, 4.44058E-2, 4.48718E-2,  &
187          &     4.29373E-2, 4.15242E-2, 3.81744E-2, 3.16342E-2, 2.99615E-2,  &
188          &     2.92740E-2, 2.67484E-2, 1.76904E-2, 1.40049E-2, 1.46224E-2,  &
189          &     1.39993E-2, 1.19574E-2, 1.06386E-2, 1.00980E-2, 8.63808E-3,  &
190          &     6.52736E-3, 4.99410E-3, 4.39350E-3, 2.21676E-3, 1.33812E-3,  &
191          &     1.12320E-3, 5.59000E-4, 3.60000E-4, 2.98080E-4, 7.46294E-5  /
192     
193     ! --------------------------------------------------------------------- !
194     !   section-2 : module variables for stratospheric volcanic aerosols    !
195     !               from historical data (sato et al. 1993)                 !
196     ! --------------------------------------------------------------------- !
197     
198     !  ---  parameter constants:
199           integer, parameter :: MINVYR = 1850    ! lower lim (year) data available
200           integer, parameter :: MAXVYR = 1999    ! upper lim (year) data available
201     
202     !  ---  monthly, 45-deg lat-zone aerosols data set in subroutine 'aerinit'
203           integer, allocatable :: ivolae(:,:,:)
204     
205     !  ---  static control variables:
206           integer, save :: kyrstr  = 0
207           integer, save :: kyrend  = 0
208           integer, save :: kyrsav  = 0
209           integer, save :: kmonsav = 0
210     
211     ! --------------------------------------------------------------------- !
212     !   section-3 : module variables for opac climatological aerosols       !
213     !               optical properties (hess et al. 1989)                   !
214     ! --------------------------------------------------------------------- !
215     
216     !  ---  parameters and constants:
217           integer, parameter :: NXC = 5    ! num of max componets in a profile
218           integer, parameter :: NAE = 7    ! num of aerosols profile structures
219           integer, parameter :: NDM = 5    ! num of atmos aerosols domains
220           integer, parameter :: IMXAE = 72 ! num of lon-points in glb aeros data set
221           integer, parameter :: JMXAE = 37 ! num of lat-points in glb aeros data set
222           integer, parameter :: NAERBND=61 ! num of bands for clim aer data (opac)
223           integer, parameter :: NRHLEV =8  ! num of rh levels for rh-dep components
224           integer, parameter :: NCM1 = 6   ! num of rh independent aeros species
225           integer, parameter :: NCM2 = 4   ! num of rh dependent aeros species
226           integer, parameter :: NCM  = NCM1+NCM2
227     
228           real (kind=kind_phys), dimension(NRHLEV) :: rhlev
229           data  rhlev (:) / 0.0, 0.5, 0.7, 0.8, 0.9, 0.95, 0.98, 0.99 /
230     
231     !  ---  the following arrays are for climatological data that are
232     !           allocated and read in subroutine 'clim_aerinit'.
233     !   - global aerosol distribution:
234     !      haer(NDM,NAE)    - scale height of aerosols (km)
235     !      prsref(NDM,NAE)  - ref pressure lev (sfc to toa) in mb (100Pa)
236     !      sigref(NDM,NAE)  - ref sigma lev (sfc to toa)
237     
238           real (kind=kind_phys), allocatable, save, dimension(:,:) :: haer
239           real (kind=kind_phys), allocatable, save, dimension(:,:) :: prsref
240     
241     !  ---  the following arrays are allocate and setup in subr 'clim_aerinit'
242     !   - for relative humidity independent aerosol optical properties:
243     !      species : insoluble        (inso); soot             (soot);
244     !                mineral nuc mode (minm); mineral acc mode (miam);
245     !                mineral coa mode (micm); mineral transport(mitr).
246     !      extrhi(NCM1,NBDSWLW) - extinction coefficient for sw+lw spectral band
247     !      scarhi(NCM1,NBDSWLW) - scattering coefficient for sw+lw spectral band
248     !      ssarhi(NCM1,NBDSWLW) - single scattering albedo for sw+lw spectral band
249     !      asyrhi(NCM1,NBDSWLW) - asymmetry parameter for sw+lw spectral band
250     !   - for relative humidity dependent aerosol optical properties:
251     !      species : water soluble    (waso); sea salt acc mode(ssam);
252     !                sea salt coa mode(sscm); sulfate droplets (suso).
253     !      rh level: 00%, 50%, 70%, 80%, 90%, 95%, 98%, 99%
254     !      extrhd(NRHLEV,NCM2,NBDSWLW) - extinction coefficient for sw+lw band
255     !      scarhd(NRHLEV,NCM2,NBDSWLW) - scattering coefficient for sw+lw band
256     !      ssarhd(NRHLEV,NCM2,NBDSWLW) - single scattering albedo for sw+lw band
257     !      asyrhd(NRHLEV,NCM2,NBDSWLW) - asymmetry parameter for sw+lw band
258     !   - for stratospheric aerosols optical properties:
259     !      extstra(NBDSWLW)            - extinction coefficient for sw+lw band
260     !   - for topospheric aerosol profile distibution:
261     !      kprfg (    IMXAE*JMXAE)   - aeros profile index
262     !      idxcg (NXC*IMXAE*JMXAE)   - aeros component index
263     !      cmixg (NXC*IMXAE*JMXAE)   - aeros component mixing ratio
264     !      denng (NXC*IMXAE*JMXAE)   - aerosols number density
265     
266           real (kind=kind_phys), allocatable, save, dimension(:,:)   ::     &
267          &       extrhi, scarhi, ssarhi, asyrhi
268           real (kind=kind_phys), allocatable, save, dimension(:,:,:) ::     &
269          &       extrhd, scarhd, ssarhd, asyrhd
270           real (kind=kind_phys), allocatable, save, dimension(:)     ::     &
271          &       extstra
272     
273           real (kind=kind_phys),allocatable,save:: cmixg(:,:,:),denng(:,:,:)
274           integer,              allocatable,save:: kprfg(:,:),  idxcg(:,:,:)
275     
276     !  ---  logical parameter for clim opac optic prop input control
277           logical, save :: lclmin = .true.
278     
279     ! =======================================================================
280     ! GOCART code modification starts here (Sarah lu)  ---------------------!
281     
282     ! --------------------------------------------------------------------- !
283     !   section-4 : module variables for gocart aerosol optical properties  !
284     ! --------------------------------------------------------------------- !
285     
286     !  ---  parameters and constants:
287     !   - KCM, KCM1, KCM2 are determined from subroutine 'set_aerspc'
288           integer, parameter :: KAERBND=61 ! num of bands for aer data (gocart)
289           integer, parameter :: KRHLEV =36 ! num of rh levels for rh-dep components
290     !*    integer, parameter :: KCM1 = 8   ! num of rh independent aer species
291     !*    integer, parameter :: KCM2 = 5   ! num of rh dependent aer species
292     !*    integer, parameter :: KCM  = KCM1 + KCM2
293           integer, save      :: KCM1 = 0   ! num of rh indep aerosols (set in subr set_aerspc)
294           integer, save      :: KCM2 = 0   ! num of rh dep aerosols   (set in subr set_aerspc)
295           integer, save      :: KCM        ! =KCM1+KCM2               (set in subr set_aerspc)
296     
297           real (kind=kind_phys), dimension(KRHLEV) :: rhlev_grt
298           data  rhlev_grt (:)/ .00, .05, .10, .15, .20, .25, .30, .35,
299          &                     .40, .45, .50, .55, .60, .65, .70, .75,
300          &                     .80, .81, .82, .83, .84, .85, .86, .87,
301          &                     .88, .89, .90, .91, .92, .93, .94, .95,
302          &                     .96, .97, .98, .99/
303     
304     !  --- tabulated aerosol optical properties (input dataset)
305     !  --- allocated and read in subroutine 'rd_gocart_luts'
306     !
307     !   - spectral band structure:
308     !      iendwv_grt(KAERBND)      - ending wavenumber (cm**-1) for each band
309     !   - relative humidity independent aerosol optical properties:
310     !   ===> species : dust (8 bins)
311     !      rhidext0_grt(KAERBND,KCM1) - extinction coefficient
312     !      rhidssa0_grt(KAERBND,KCM1) - single scattering albedo
313     !      rhidasy0_grt(KAERBND,KCM1) - asymmetry parameter
314     !   - relative humidity dependent aerosol optical properties:
315     !   ===> species : soot, suso, waso, ssam, sscm
316     !      rhdpext0_grt(KAERBND,KRHLEV,KCM2) - extinction coefficient
317     !      rhdpssa0_grt(KAERBND,KRHLEV,KCM2) - single scattering albedo
318     !      rhdpasy0_grt(KAERBND,KRHLEV,KCM2) - asymmetry parameter
319     
320           integer,               allocatable, dimension(:) :: iendwv_grt
321           real (kind=kind_phys), allocatable, dimension(:,:)  ::            &
322          &                       rhidext0_grt, rhidssa0_grt, rhidasy0_grt
323           real (kind=kind_phys), allocatable, dimension(:,:,:)::            &
324          &                       rhdpext0_grt, rhdpssa0_grt, rhdpasy0_grt
325     
326     !  --- aerosol optical properties mapped onto scheme-specific spectral bands
327     !  --- allocated and computed in subroutine 'optavg_grt'
328     
329     !   - relative humidity independent aerosol optical properties:
330     !      extrhi_grt(KCM1,NBDSWLW) - extinction coefficient for sw+lw spectral band
331     !      ssarhi_grt(KCM1,NBDSWLW) - single scattering albedo for sw+lw spectral band
332     !      asyrhi_grt(KCM1,NBDSWLW) - asymmetry parameter for sw+lw spectral band
333     !   - relative humidity dependent aerosol optical properties:
334     !      extrhd_grt(KRHLEV,KCM2,NBDSWLW) - extinction coefficient for sw+lw band
335     !      ssarhd_grt(KRHLEV,KCM2,NBDSWLW) - single scattering albedo for sw+lw band
336     !      asyrhd_grt(KRHLEV,KCM2,NBDSWLW) - asymmetry parameter for sw+lw band
337     
338           real (kind=kind_phys), allocatable, save, dimension(:,:)   ::     &
339          &       extrhi_grt, ssarhi_grt, asyrhi_grt
340           real (kind=kind_phys), allocatable, save, dimension(:,:,:) ::     &
341          &       extrhd_grt, ssarhd_grt, asyrhd_grt
342     
343     ! --------------------------------------------------------------------- !
344     !   section-5 : module variables for gocart aerosol climo data set      !
345     ! --------------------------------------------------------------------- !
346     
347     !
348     !     This version only supports geos3-gocart data set (Jan 2010)
349     !     Modified to support geos4-gocart data set        (May 2010)
350     !
351     !  geos3-gocart vs geos4-gocart
352     !  (1) Use the same module variables
353     !      IMXG,JMXG,KMXG,NMXG,psclmg,dmclmg,geos_rlon,geos_rlat
354     !  (2) Similarity between geos3 and geos 4:
355     !      identical lat/lon grids and aerosol specification;
356     !      direction of vertical index is bottom-up (sfc to toa)
357     !  (3) Difference between geos3 and geos4
358     !      vertical coordinate (sigma for geos3/hybrid_sigma_pressure for geos4)
359     !      aerosol units (mass concentration for geos3/mixing ratio for geos4)
360     
361           integer, parameter :: IMXG = 144 ! num of lon-points in geos dataset
362           integer, parameter :: JMXG = 91  ! num of lat-points in geos dataset
363           integer, parameter :: KMXG = 30  ! num of vertical layers in geos dataset
364     !*    integer, parameter :: NMXG = 12  ! num of gocart aer spec for opt calc
365           integer, save      :: NMXG       ! to be determined by set_aerspc
366     
367           real (kind=kind_phys), parameter :: dltx = 360.0 / float(IMXG)
368           real (kind=kind_phys), parameter :: dlty = 180.0 / float(JMXG-1)
369     
370     !  --- the following arrays are allocated and setup in 'rd_gocart_clim'
371     !   - geos-gocart climo data (input dataset)
372     !     psclmg  - pressure in cb                   IMXG*JMXG*KMXG
373     !     dmclmg  - aerosol dry mass in g/m3         IMXG*JMXG*KMXG*NMXG
374     !               or aerosol mixing ratio in mol/mol or Kg/Kg
375     
376           real (kind=kind_phys),allocatable, save:: psclmg(:,:,:),          &
377          &                                          dmclmg(:,:,:,:)
378     
379     !   - geos-gocart lat/lon arrays
380           real (kind=kind_phys), allocatable, save, dimension(:)   ::       &
381          &                       geos_rlon, geos_rlat
382     
383     !   - control flag for gocart climo data set
384     !     xxxx as default; ver3 for geos3; ver4 for geos4; 0000 for unknown data
385           character*4, save  :: gocart_climo = 'xxxx'
386     
387     !   - molecular wght of gocart aerosol species
388           real (kind=kind_io4), allocatable :: molwgt(:)
389     
390     ! --------------------------------------------------------------------- !
391     !   section-6 : module variables for gocart aerosol scheme options      !
392     ! --------------------------------------------------------------------- !
393     
394     !  ---  logical parameter for gocart initialization control
395           logical, save :: lgrtint = .true.
396     
397     !  ---  logical parameter for gocart debug print control
398     !     logical, save :: lckprnt = .true.
399           logical, save :: lckprnt = .false.
400     
401     !  --- the following index/flag/weight are set up in 'set_aerspc'
402     
403     !   - merging coefficients for fcst/clim; determined from fdaer
404           real (kind=kind_phys), save :: ctaer = f_zero   ! user specified wgt
405     
406     !   - option to get fcst or clim gocart aerosol fields
407           logical, save :: get_fcst = .true.    ! option to get fcst field
408           logical, save :: get_clim = .true.    ! option to get clim field
409     
410     !  ------  gocart aerosol specification    ------
411     !  =>  transported aerosol species:
412     !      DU (5-bins)
413     !      SS (4 bins for climo mode and 5 bins for fcst mode)
414     !      SU (dms, so2, so4, msa)
415     !      OC (phobic, philic) and BC (phobic, philic)
416     !  =>  species and lumped species for aerosol optical properties
417     !      DU (5-bins, with 4 sub-groups in the submicron bin )
418     !      SS (ssam for submicron, sscm for coarse mode)
419     !      SU (so4)
420     !      OC (phobic, philic) and BC (phobic, philic)
421     !  =>  specification used for aerosol optical properties luts
422     !      DU (8 bins)
423     !      SS (ssam, sscm)
424     !      SU (suso)
425     !      OC (waso) and BC (soot)
426     !
427     
428     !   - index for rh dependent aerosol optical properties (2nd
429     !     dimension for extrhd_grt, ssarhd_grt, and asyrhd_grt)
430           integer, save :: isoot, iwaso, isuso, issam, isscm
431     
432     !   - index for rh independent aerosol optical properties (1st
433     !     dimension for extrhi_grt, ssarhi_grt, and asyrhi_grt) is
434     !     not needed ===> hardwired to 8-bin dust
435     
436     !   - index for gocart aerosol species to be included in the
437     !     calculations of aerosol optical properties (ext, ssa, asy)
438           type  gocart_index_type
439              integer :: dust1, dust2, dust3, dust4, dust5,    ! dust
440          &              ssam,  sscm,                          ! sea salt
441          &              suso,                                 ! sulfate
442          &              waso_phobic, waso_philic,             ! oc
443          &              soot_phobic, soot_philic              ! bc
444           endtype
445           type (gocart_index_type), save :: dm_indx
446     
447     !   - index for gocart aerosols from prognostic tracer fields
448           type  tracer_index_type
449              integer :: du001, du002, du003, du004, du005,    ! dust
450          &              ss001, ss002, ss003, ss004, ss005,    ! sea salt
451          &              so4,                                  ! sulfate
452          &              ocphobic, ocphilic,                   ! oc
453          &              bcphobic, bcphilic                    ! bc
454           endtype
455           type (tracer_index_type), save :: dmfcs_indx
456     
457     !   - grid components to be included in the aeropt calculations
458           integer, save                  :: num_gridcomp = 0  ! number of aerosol grid components
459           character, allocatable , save  :: gridcomp(:)*2     ! aerosol grid components
460     
461     !   - spectral band for 550nm (for diag)
462           integer, public, save       :: nv_aod = 1
463     
464     !  --- default full-package setting
465           integer, parameter          :: max_num_gridcomp = 5
466           character*2                 :: max_gridcomp(max_num_gridcomp)
467           data max_gridcomp  /'DU', 'BC', 'OC', 'SU', 'SS'/
468     
469     ! GOCART code modification end here (Sarah Lu)  ------------------------!
470     ! =======================================================================
471     
472     !  ---  public interfaces
473     
474           public aerinit, setaer
475     
476     
477     ! =================
478           contains
479     ! =================
480     
481     !-----------------------------------
482           subroutine aerinit                                                &
483     !...................................
484     
485     !  ---  inputs:
486          &     ( iyear, imon, IAER, me, raddt, fdaer )
487     !  ---  outputs: ( none )
488     
489     !  ==================================================================  !
490     !                                                                      !
491     !  aerinit does invoke aerosol initialization programs based on the    !
492     !  selection of schemes.                                               !
493     !                                                                      !
494     !  inputs:                                                             !
495     !     iyear   - 4-digit calender year                 1                !
496     !     imon    - month of the year                     1                !
497     !     IAER    - 3-digit aerosol flag (volc,lw,sw)     1                !
498     !             (a) opac - vol option implemented by Y-T Hou             !
499     !             (b) gocart option implemented by C-H Lu                  !
500     !               =  0: turn all aeros effects off (sw,lw,volc)          !
501     !               =  1: use clim tropspheric aerosol for sw only         !
502     !               = 10: use clim tropspheric aerosol for lw only         !
503     !               = 11: use clim tropspheric aerosol for both sw and lw  !
504     !               =100: volc aerosol only for both sw and lw             !
505     !               =101: volc and clim trops aerosol for sw only          !
506     !               =110: volc and clim trops aerosol for lw only          !
507     !               =111: volc and clim trops aerosol for both sw and lw   !
508     !               =  2: gocart tropspheric aerosol for sw only           !
509     !               = 20: gocart tropspheric aerosol for lw only           !
510     !               = 22: gocart tropspheric aerosol for both sw and lw    !
511     !               =102: volc and gocart trops aerosol for sw only        !
512     !               =120: volc and gocart trops aerosol for lw only        !
513     !               =122: volc and gocart trops aerosol for both sw and lw !
514     !     me      - print message control flag            1                !
515     !                                                                      !
516     !  outputs: (to the module variables)                                  !
517     !    ( none )                                                          !
518     !                                                                      !
519     !  module variables:                                                   !
520     !     * opac climo *                                                   !
521     !     kprfg   - aerosols profile index                IMXAE*JMXAE      !
522     !     idxcg   - aerosols component index              NXC*IMXAE*JMXAE  !
523     !     cmixg   - aerosols component mixing ratio       NXC*IMXAE*JMXAE  !
524     !     denng   - aerosols number density               NXC*IMXAE*JMXAE  !
525     !     * geos-gocart climo (from sfc to toa) *                          !
526     !     psclmg  - pressure in cb                      IMXG*JMXG*KMXG     !
527     !     dmclmg  - aerosol dry mass in g/m3            IMXG*JMXG*KMXG*NMXG!
528     !               or aerosol mixing ratio in mol/mol or Kg/Kg            !
529     !                                                                      !
530     !     ivolae  - stratosphere volcanic aerosol optical depth (fac 1.e4) !
531     !                                                     12*4*10          !
532     !                                                                      !
533     !  usage:    call aerinit                                              !
534     !                                                                      !
535     !  subprograms called:  clim_aerinit                                   !
536     !                       gocart_init                                    !
537     !                                                                      !
538     !  ==================================================================  !
539     !
540           implicit none
541     
542     !  ---  inputs:
543           integer,  intent(in) :: iyear, imon, IAER, me
544           real (kind=kind_phys), intent(in) :: raddt, fdaer
545     
546     !  ---  output: ( none )
547     
548     !  ---  locals:
549           real (kind=kind_phys), dimension(NWVTOT) :: solfwv        ! one wvn sol flux
550           real (kind=kind_phys), dimension(NWVTIR) :: eirfwv        ! one wvn ir flux
551           real (kind=kind_phys) :: soltot, tmp1, tmp2, tmp3
552     
553           integer :: nb, ni, nw, nw1, nw2, nmax, nmin
554     
555           character            :: cline*80, ctyp*3, volcano_file*32
556     
557           integer :: i, j, k, nc, iy, id, ilw, isw
558           logical :: file_exist
559     
560           data volcano_file / 'volcanic_aerosols_1850-1859.txt ' /
561     
562     !===>  ...  begin here
563     
564           isw = mod(IAER,10)                ! trop-aer scheme for sw
565           iy  = IAER / 10
566           ilw = mod(iy , 10)                ! trop-aer scheme for lw
567     
568           iaerflg = max( isw, ilw )         ! flag for trop-aer scheme selection
569           lalwflg = ilw  > 0                ! flag for lw trop-aer properties
570           laswflg = isw  > 0                ! flag for sw trop-aer properties
571           lvolflg = IAER >= 100             ! flag for stratospheric volcanic aer
572     
573     !  --- ...  in sw, aerosols optical properties are computed for each radiation
574     !           spectral band; while in lw, optical properties can be calculated
575     !           for either only one broad band or for each of the lw radiation bands
576     
577           if ( iaerlw == 1 ) then
578             NBDIR = NBDLW
579           else
580             NBDIR = 1
581           endif
582           NBDSWLW = NBDSW + NBDIR
583     
584           if ( lckprnt .and. (me==0) )                                      &
585          &  write(0,*)'RAD_NBDSW,NBDIR,NBDSWLW:', NBDSW,NBDIR,NBDSWLW
586     
587     !  --- ...  define the one wavenumber solar fluxes based on toa solar
588     !           spectral distribution
589     
590           nmax = min( NWVTOT, nint( maxval(wvnum2) ))
591           nmin = max( 1,      nint( minval(wvnum1) ))
592     
593     !     print *,' MINWVN, MAXWVN = ',nmin, nmax
594     
595     !     soltot1 = f_zero
596           soltot  = f_zero
597           do nb = 1, NWVSOL
598             if ( nb == 1 ) then
599               nw1 = 1
600             else
601               nw1 = nw1 + nwvns0(nb-1)
602             endif
603     
604             nw2 = nw1 + nwvns0(nb) - 1
605     
606             do nw = nw1, nw2
607               solfwv(nw) = s0intv(nb)
608     !         soltot1 = soltot1 + s0intv(nb)
609               if ( nw >= nmin .and. nw <= nmax ) then
610                 soltot = soltot + s0intv(nb)
611               endif
612             enddo
613           enddo
614     
615     !  --- ...  define the one wavenumber ir fluxes based on black-body
616     !           emission distribution at a predefined temperature
617     
618           tmp1 = 2.0 * con_pi * con_plnk * (con_c**2)
619           tmp2 = con_plnk * con_c / (con_boltz * con_t0c)
620     
621           do ni = 1, NWVTIR
622             tmp3 = 100.0 * ni
623             eirfwv(ni) = (tmp1 * tmp3**3) / (exp(tmp2*tmp3) - 1.0)
624           enddo
625     
626     !  --- ...  write aerosol parameter configuration to output logs
627     
628           if ( me == 0 ) then
629     
630             write(0,*)'   IAER=',IAER,'  iaerflg=',iaerflg,'  LW-trop-aer=' &
631          &         ,lalwflg,'  SW-trop-aer=',laswflg,'  Volc-aer=',lvolflg
632     
633             if ( IAER <= 0 ) then        ! turn off all aerosol effects
634     
635               print *,' - No tropospheric/volcanic aerosol effect included'
636               print *,'      Input values of aerosol optical properties to' &
637          &           ,' both SW and LW radiations are set to zeros'
638     
639             elseif ( IAER == 100 ) then  ! only stratospheric volcanic aerosols
640     
641               print *,' - Include onle volcanic aerosols in both SW and LW' &
642          &           ,' for year, month =', iyear, imon
643     
644             else
645     
646               if ( IAER < 100 ) then       ! no stratospheric volcanic aerosols
647                 print *,' - No stratospheric volcanic aerosol effect'
648               else                         ! include stratospheric volcanic aerosols
649                 print *,' - Include stratospheric volcanic aerosol effect'  &
650          &             ,' for year, month =', iyear, imon
651               endif
652     
653               if ( iaerflg == 1 ) then      ! opac tropospheric climatological
654     
655                 print *,' - Using OPAC climatology for tropospheric aerosol'
656     
657               elseif ( iaerflg == 2 ) then  ! gocart tropospheric climatological
658     
659                 print *,' - Using GOCART scheme for tropospheric aerosol'
660     
661               endif                         ! end if_iaerflg_block
662             
663               if ( laswflg ) then           ! shcek for sw effect
664                 print *,'      Compute aerosol optical properties for SW'   &
665          &             ,' input parameters'
666               else
667                 print *,'      No SW radiation aerosol effect, values of'   &
668          &             ,' aerosol properties to SW input are set to zeros'
669               endif                         ! end if_laswflg_block
670     
671               if ( lalwflg ) then           ! check for lw effect
672                 print *,'      Compute aerosol optical properties for LW'   &
673          &             ,' input parameters'
674               else
675                 print *,'      No LW radiation aerosol effect, values of'   &
676          &             ,' aerosol properties to LW input are set to zeros'
677               endif                         ! end if_lalwflg_block
678     
679             endif     ! end if_IAER_block
680     
681           endif       ! end if_me_block
682     
683     !  --- ...  tropospheric aerosol initialization
684     
685           if ( IAER == 0 ) then 
686     
687             return
688     
689           elseif ( IAER /= 100 ) then
690     
691     !* imon_check applies to both opac and gocart aerosol schemes          !chlu
692     !*      if ( iaerflg == 1 ) then      ! opac tropospheric climatology  !chlu
693     
694               if ( imon < 1 .or. imon > 12 ) then
695                 print *,' ***** ERROR in specifying requested month!!! ',   &
696          &              'imon=', imon
697                 print *,' ***** STOPPED in subroutinte AERINIT!!!'
698                 stop
699               endif
700     
701             if ( iaerflg == 1 ) then      ! opac tropospheric climatology
702               call clim_aerinit                                             &
703     !  ---  inputs:
704          &     ( NWVTOT,solfwv,soltot,NWVTIR,eirfwv,                        &
705          &       NBDSW,NBDIR,NBDSWLW, imon, me                              &
706     !  ---  outputs:  (none)
707          &     )
708     
709     
710             elseif ( iaerflg == 2 ) then  ! gocart prognostic aerosols
711               call gocart_init                                              &
712     !  ---  inputs:
713          &     ( NWVTOT,solfwv,soltot,NWVTIR,eirfwv,                        &
714          &       NBDSW,NBDIR,NBDSWLW,imon,me,raddt,fdaer                    &
715     !  ---  outputs:  (none)
716          &     )
717     
718             else
719               print *,'   ERROR in aerosols specification! IAER =',iaer
720               print *,'     iaerflg, lvolflg =',iaerflg,lvolflg
721               print *,'   *** Stopped in subroutine AERINIT !!'
722               stop
723             endif                          ! end if_iaerflg_block
724     
725           endif    ! end if_IAER_block
726     
727     !  --- ...  stratosperic volcanic aerosol initialization
728     
729           if ( lvolflg ) then
730     
731     !  ---  allocate data space
732     
733             if ( .not. allocated(ivolae) ) then
734               allocate ( ivolae(12,4,10) )   ! for 12-mon,4-lat_zone,10-year
735             endif
736     
737             kmonsav = imon
738     
739             if ( kyrstr<=iyear .and. iyear<=kyrend ) then   ! use previously input data
740               kyrsav = iyear
741               return
742             else                                            ! need to input new data
743               kyrsav = iyear
744               kyrstr = iyear - mod(iyear,10)
745               kyrend = kyrstr + 9
746               if ( lckprnt .and. (me==0) )  then
747               print *,'  kyrstr, kyrend, kyrsav, kmonsav =',                &
748          &            kyrstr,kyrend,kyrsav,kmonsav
749               endif
750     
751               if ( iyear < MINVYR .or. iyear > MAXVYR ) then
752                 ivolae(:,:,:) = 1            ! set as lowest value
753                 if ( me == 0 ) then
754                   print *,'   Request volcanic date out of range,',         &
755          &                ' optical depth set to lowest value'
756                 endif
757               else
758                 write(volcano_file(19:27),60) kyrstr,kyrend
759       60        format(i4.4,'-',i4.4)
760     
761                 inquire (file=volcano_file, exist=file_exist)
762                 if ( file_exist ) then
763                   open (unit=NIAERCM,file=volcano_file,status='OLD',        &
764          &              form='FORMATTED')
765     
766                   read(NIAERCM,62) cline
767       62          format(a80)
768     
769     !  ---  check print
770                   if ( me == 0 ) then
771                     print *,'   Opened volcanic data file: ',volcano_file
772                     print *, cline
773                   endif
774     
775                   do k = 1, 10
776                     do j = 1, 4
777                       read(NIAERCM,64) ivolae(:,j,k)
778     !                 read(NIAERCM,64) (ivolae(i,j,k),i=1,12)
779       64              format(12i5)
780                     enddo
781                   enddo
782     
783                   close (NIAERCM)
784                 else
785                   print *,'   Requested volcanic data file "',              &
786          &                volcano_file,'" not found!'
787                   print *,'   *** Stopped in subroutine AERINIT !!'
788                   stop
789                 endif              ! end if_file_exist_block
790               endif                ! end if_iyear_block
791     
792             endif              ! end if_kyrstr_block
793     
794             if ( me == 0 ) then
795               iy = mod(kyrsav,10) + 1
796               print *,' CHECK: Sample Volcanic data used for month, year:', &
797          &             imon, iyear
798               print *,  ivolae(kmonsav,:,iy)
799             endif
800           endif                ! end if_lvolflg_block
801     !
802           return
803     !...................................
804           end subroutine aerinit
805     !-----------------------------------
806     
807     
808     
809     !-----------------------------------
810           subroutine setaer                                                 &
811     !...................................
812     
813     !  ---  inputs:
814          &     ( xlon,xlat,prsi,prsl,tlay,qlay,rhlay,                       &
815          &       prslk, ozlay,                                              &
816          &       IMAX,NLAY,NLP1, iflip, lsswr,lslwr,                        &
817     !  ---  outputs:
818          &       aerosw,aerolw                                              &
819          &,      tau_gocart                                                 &
820          &     )
821     
822     !  ==================================================================  !
823     !                                                                      !
824     !  setaer computes aerosols optical properties from different global   !
825     !  aerosols data sets based on scheme selection.                       !
826     !                                                                      !
827     !  inputs:                                                             !
828     !     xlon, xlat                                             IMAX      !
829     !             - longitude and latitude of given points in radiance     !
830     !     prsi    - pressure at interface              mb   IMAX*NLP1      !
831     !     prsl    - layer mean pressure                mb   IMAX*NLAY      !
832     !     tlay    - layer mean temperature             k    IMAX*NLAY      !
833     !     qlay    - layer mean specific humidity       g/g  IMAX*NLAY      !
834     !     rhlay   - layer mean relative humidity            IMAX*NLAY      !
835     !     prslk   - pressure                           cb   IMAX*NLAY      !
836     !     ozlay   - layer tracer mass mixing ratio   g/g   IMAX*NLAY*NTRAC !
837     !     IMAX    - horizontal dimension of arrays                  1      !
838     !     NLAY,NLP1-vertical dimensions of arrays                   1      !
839     !     iflip   - control flag for direction of vertical index    1      !
840     !               =0: index from toa to surface                          !
841     !               =1: index from surface to toa                          !
842     !     lsswr,lslwr                                                      !
843     !             - logical flags for sw/lw radiation calls         1      !
844     !                                                                      !
845     !  outputs:                                                            !
846     !     aerosw - aeros opt properties for sw      IMAX*NLAY*NBDSW*NF_AESW!
847     !               (:,:,:,1): optical depth                               !
848     !               (:,:,:,2): single scattering albedo                    !
849     !               (:,:,:,3): asymmetry parameter                         !
850     !     aerolw - aeros opt properties for lw      IMAX*NLAY*NBDLW*NF_AELW!
851     !               (:,:,:,1): optical depth                               !
852     !               (:,:,:,2): single scattering albedo                    !
853     !               (:,:,:,3): asymmetry parameter                         !
854     !     tau_gocart - 550nm aeros opt depth     IMAX*NLAY*MAX_NUM_GRIDCOMP!
855     !                                                                      !
856     !                                                                      !
857     !  module variable: (set by subroutine aerinit)                        !
858     !     iaerflg - control flag for tropospheric aerosols selection       !
859     !               =0: do not calc tropospheric aerosol optical properties!
860     !               =1: use opac tropospheric aerosol climatology          !
861     !               =2: use gocart tropospheric interactive aerosol        !
862     !     lalwflg - control flag for lw radiation aerosols effect          !
863     !               =f: do not calc lw aerosol optical properties          !
864     !               =t: calculate lw aerosol optical properties            !
865     !     laswflg - control flag for sw radiation aerosols effect          !
866     !               =f: do not calc sw aerosol optical properties          !
867     !               =t: calculate sw aerosol optical properties            !
868     !     lvolflg - control flag for stratospheric vocanic aerosols        !
869     !               =t: add volcanic aerosols to the background aerosols   !
870     !               =f: do not add volcanic aerosols                       !
871     !                                                                      !
872     !     ivolae  - stratosphere volcanic aerosol optical depth (fac 1.e4) !
873     !                                                     12*4*10          !
874     !                                                                      !
875     !  usage:    call setaer                                               !
876     !                                                                      !
877     !  subprograms called:  setclimaer                                     !
878     !                       setgocartaer                                   !
879     !                                                                      !
880     !  ==================================================================  !
881     !
882           implicit none
883     
884     !  ---  inputs:
885           integer, intent(in) :: IMAX,NLAY,NLP1, iflip
886     
887           real (kind=kind_phys), dimension(:,:), intent(in) :: prsi, prsl,  &
888          &       tlay, qlay, rhlay
889           real (kind=kind_phys), dimension(:),   intent(in) :: xlon, xlat
890           logical, intent(in) :: lsswr, lslwr
891     
892     !     Added for gocart coupling (Sarah Lu)
893           real (kind=kind_phys), dimension(:,:), intent(in) :: prslk
894           real (kind=kind_phys), dimension(:,:,:), intent(in) :: ozlay
895     
896     !  ---  outputs:
897           real (kind=kind_phys), dimension(:,:,:,:), intent(out) ::         &
898          &       aerosw, aerolw
899     
900     !     Added for aerosol diag (Sarah Lu)
901           real (kind=kind_phys), dimension(:,:,:), intent(out) :: tau_gocart
902     
903     !  ---  locals:
904           real (kind=kind_phys), dimension(IMAX) :: alon, alat, volcae, delp
905           real (kind=kind_phys) :: prsln(NLP1),hz(IMAX,NLP1),dz(IMAX,NLAY)
906           real (kind=kind_phys) :: tmp1, tmp2, psrfh, psrfl
907     
908           integer               :: kcutl(IMAX), kcuth(IMAX)
909           integer               :: i, i1, k, m, mb, kp, kh, kl
910     
911           logical               :: laddsw = .false.
912           logical               :: laddlw = .false.
913     
914     !  ---  conversion constants
915           real (kind=kind_phys), parameter :: rdg  = 180.0 / con_pi
916           real (kind=kind_phys), parameter :: rovg = 0.001 * con_rd / con_g
917     
918     
919     !===>  ...  begin here
920     
921           if ( .not. (lsswr .or. lslwr) ) then
922             return
923           endif
924     
925           if ( .not. laswflg ) then
926             aerosw = f_zero
927           endif
928     
929           if ( .not. lalwflg ) then
930             aerolw = f_zero
931           endif
932     
933           if ( (.not.lvolflg) .and. (iaerflg==0) ) then
934             return
935           endif
936     
937     !  ---  ...  convert lat/lon from radiance to degree
938     
939           do i = 1, IMAX
940             alon(i) = xlon(i) * rdg
941             if (alon(i) < f_zero) alon(i) = alon(i) + 360.0
942             alat(i) = xlat(i) * rdg
943           enddo
944     
945     !  ---  ...  compute level height and layer thickness
946     
947           lab_do_IMAX : do i = 1, IMAX
948     
949             lab_if_flip : if (iflip == 1) then        ! input from sfc to toa
950     
951               do k = 1, NLAY
952                 prsln(k) = log(prsi(i,k))
953               enddo
954               prsln(NLP1)= log(prsl(i,NLAY))
955     
956               do k = NLAY, 1, -1
957                 dz(i,k) = rovg * (prsln(k) - prsln(k+1))                    &
958          &              * tlay(i,k) * (f_one + con_fvirt*qlay(i,k))
959               enddo
960               dz(i,NLAY)  = 2.0 * dz(i,NLAY)
961     
962               hz(i,1) = f_zero
963               do k = 1, NLAY
964                 hz(i,k+1) = hz(i,k) + dz(i,k)
965               enddo
966     
967             else  lab_if_flip                         ! input from toa to sfc
968     
969               prsln(1) = log(prsl(i,1))
970               do k = 2, NLP1
971                 prsln(k) = log(prsi(i,k))
972               enddo
973     
974               do k = 1, NLAY
975                 dz(i,k) = rovg * (prsln(k+1) - prsln(k))                    &
976          &              * tlay(i,k) * (f_one + con_fvirt*qlay(i,k))
977               enddo
978               dz(i,1) = 2.0 * dz(i,1)
979     
980               hz(i,NLP1) = f_zero
981               do k = NLAY, 1, -1
982                 hz(i,k) = hz(i,k+1) + dz(i,k)
983               enddo
984     
985             endif  lab_if_flip
986     
987           enddo  lab_do_IMAX
988     
989     !  ---  ...  calculate sw aerosol optical properties for the corresponding
990     !            frequency bands based on scheme selection
991     
992           tau_gocart(:,:,:) = f_zero
993           if ( iaerflg == 1 ) then      ! use opac aerosol climatology
994     
995             call setclimaer                                                 &
996     !  ---  inputs:
997          &     ( alon,alat,prsi,rhlay,dz,hz,NBDSWLW,                        &
998          &       IMAX,NLAY,NLP1, iflip, lsswr,lslwr,                        &
999     !  ---  outputs:
1000          &       aerosw,aerolw                                              &
1001          &     )
1002     
1003           elseif ( iaerflg == 2 )   then    ! use gocart aerosol scheme
1004     
1005             call setgocartaer                                               &
1006     
1007     !  ---  inputs:
1008          &     ( alon,alat,prslk,rhlay,dz,hz,NBDSWLW,                       &
1009          &       prsl,tlay,qlay,ozlay,                                      &
1010          &       IMAX,NLAY,NLP1, iflip, lsswr,lslwr,                        &
1011     !  ---  outputs:
1012          &       aerosw,aerolw                                              &
1013          &,      tau_gocart                                                 &
1014          &     )
1015     
1016           endif     ! end if_iaerflg_block
1017     
1018     !  ---  check print
1019           if ( lckprnt )  then
1020             do m = 1, NBDSW
1021              print *,'  ***  CHECK AEROSOLS PROPERTIES FOR SW BAND =',m,     &
1022          &          ' ***'
1023              do k = 1, 10
1024                print *,'  LEVEL :',k
1025                print *,'  TAUAER:',aerosw(:,k,m,1)
1026                print *,'  SSAAER:',aerosw(:,k,m,2)
1027                print *,'  ASYAER:',aerosw(:,k,m,3)
1028              enddo
1029             enddo
1030             do m = 1, NBDIR
1031              print *,'  ***  CHECK AEROSOLS PROPERTIES FOR LW BAND =',m,     &
1032          &          ' ***'
1033              do k = 1, 10
1034                print *,'  LEVEL :',k
1035                print *,'  TAUAER:',aerolw(:,k,m,1)
1036                print *,'  SSAAER:',aerolw(:,k,m,2)
1037                print *,'  ASYAER:',aerolw(:,k,m,3)
1038              enddo
1039             enddo
1040           endif
1041     
1042     
1043     !  ---  ...  stratosphere volcanic forcing
1044     
1045           if ( lvolflg ) then
1046     
1047             laddsw = lsswr .and. (laswflg .or. iaerflg==0)
1048             laddlw = lslwr .and. (lalwflg .or. iaerflg==0)
1049     
1050             i1 = mod(kyrsav, 10) + 1
1051     
1052     !  ---  select data in 4 lat bands, interpolation at the boundaires
1053     
1054             do i = 1, IMAX
1055               if      ( alat(i) > 46.0 ) then
1056                 volcae(i) = 1.0e-4 * ivolae(kmonsav,1,i1)
1057               else if ( alat(i) > 44.0 ) then
1058                 volcae(i) = 5.0e-5                                          &
1059          &                * (ivolae(kmonsav,1,i1) + ivolae(kmonsav,2,i1))
1060               else if ( alat(i) >  1.0 ) then
1061                 volcae(i) = 1.0e-4 * ivolae(kmonsav,2,i1)
1062               else if ( alat(i) > -1.0 ) then
1063                 volcae(i) = 5.0e-5                                          &
1064          &                * (ivolae(kmonsav,2,i1) + ivolae(kmonsav,3,i1))
1065               else if ( alat(i) >-44.0 ) then
1066                 volcae(i) = 1.0e-4 * ivolae(kmonsav,3,i1)
1067               else if ( alat(i) >-46.0 ) then
1068                 volcae(i) = 5.0e-5                                          &
1069          &                * (ivolae(kmonsav,3,i1) + ivolae(kmonsav,4,i1))
1070               else
1071                 volcae(i) = 1.0e-4 * ivolae(kmonsav,4,i1)
1072               endif
1073             enddo
1074     
1075             if ( iflip == 0 ) then          ! input data from toa to sfc
1076     
1077               psrfh = 5.0                        ! ref press for upper bound
1078     
1079     !  ---  find lower boundary of stratosphere
1080     
1081               do i = 1, IMAX
1082     
1083                 tmp1 = abs( alat(i) )
1084                 if ( tmp1 > 70.0 ) then          ! polar, fixed at 250mb
1085                   psrfl = 250.0
1086                 elseif ( tmp1 < 20.0 ) then      ! tropic, fixed at 150mb
1087                   psrfl = 150.0
1088                 else                             ! mid-lat, interp
1089                   psrfl = 110.0 + 2.0*tmp1
1090                 endif
1091     
1092                 kcuth(i) = NLAY - 1
1093                 kcutl(i) = 2
1094                 delp(i) = prsi(i,2)
1095     
1096                 lab_do_kcuth0 : do k = 2, NLAY-2
1097                   if ( prsi(i,k) >= psrfh ) then
1098                     kcuth(i) = k - 1
1099                     exit lab_do_kcuth0
1100                   endif
1101                 enddo  lab_do_kcuth0
1102     
1103                 lab_do_kcutl0 : do k = 2, NLAY-2
1104                   if ( prsi(i,k) >= psrfl ) then
1105                     kcutl(i) = k - 1
1106                     delp(i) = prsi(i,k) - prsi(i,kcuth(i))
1107                     exit lab_do_kcutl0
1108                   endif
1109                 enddo  lab_do_kcutl0
1110               enddo
1111     
1112     !  ---  sw: add volcanic aerosol optical depth to the background value
1113     
1114               if ( laddsw ) then
1115                 do m = 1, NBDSW
1116                   mb = NSWSTR + m - 1
1117     
1118                   if     ( wvnum1(mb) > 20000 ) then   ! range of wvlth < 0.5mu
1119                     tmp2 = 0.74
1120                   elseif ( wvnum2(mb) < 20000 ) then   ! range of wvlth > 0.5mu
1121                     tmp2 = 1.14
1122                   else                                 ! range of wvlth in btwn
1123                     tmp2 = 0.94
1124                   endif
1125                   tmp1 = (0.275e-4 * (wvnum2(mb)+wvnum1(mb))) ** tmp2
1126     
1127                   do i = 1, IMAX
1128                     kh = kcuth(i)
1129                     kl = kcutl(i)
1130                     do k = kh, kl
1131                       tmp2 = tmp1 * ((prsi(i,k+1) - prsi(i,k)) / delp(i))
1132                       aerosw(i,k,m,1) = aerosw(i,k,m,1) + tmp2*volcae(i)
1133                     enddo
1134     
1135     !  ---  smoothing profile at boundary if needed
1136     
1137                     if ( aerosw(i,kl,m,1) > 10.*aerosw(i,kl+1,m,1) ) then
1138                       tmp2 = aerosw(i,kl,m,1) + aerosw(i,kl+1,m,1)
1139                       aerosw(i,kl  ,m,1) = 0.8 * tmp2
1140                       aerosw(i,kl+1,m,1) = 0.2 * tmp2
1141                     endif
1142                   enddo    ! end do_i_block
1143                 enddo      ! end do_m_block
1144     
1145     !  ---  check print
1146     
1147     !           do i = 1, IMAX
1148     !             print *,' LEV  PRESS    TAUSAV    NEWTAU    FOR PROFILE:',&
1149     !    &                i,'  KCUTH, KCUTL =',kcuth(i),kcutl(i)
1150     !             kh = kcuth(i) - 1
1151     !             kl = kcutl(i) + 10
1152     !             do k = kh, kl
1153     !               write(6,71) k, prsl(i,k), aersav(i,k), aerosw(i,k,1,1)
1154     ! 71            format(i3,f9.3,2e11.4)
1155     !             enddo
1156     !           enddo
1157               endif        ! end if_laddsw_block
1158     
1159     !  ---  lw: add volcanic aerosol optical depth to the background value
1160     
1161               if ( laddlw ) then
1162                 if ( NBDIR == 1 ) then
1163     
1164                   tmp1 = (0.55 / 11.0) ** 1.2
1165                   do i = 1, IMAX
1166                     kh = kcuth(i)
1167                     kl = kcutl(i)
1168                     do k = kh, kl
1169                       tmp2 = tmp1 * ((prsi(i,k+1) - prsi(i,k)) / delp(i))
1170                       aerolw(i,k,1,1) = aerolw(i,k,1,1) + tmp2*volcae(i)
1171                     enddo
1172                   enddo    ! end do_i_block
1173     
1174                 else
1175     
1176                   do m = 1, NBDIR
1177                     tmp1 = (0.275e-4 * (wvnlw2(m) + wvnlw1(m))) ** 1.2
1178     
1179                     do i = 1, IMAX
1180                       kh = kcuth(i)
1181                       kl = kcutl(i)
1182                       do k = kh, kl
1183                         tmp2 = tmp1 * ((prsi(i,k+1)-prsi(i,k)) / delp(i))
1184                         aerolw(i,k,m,1) = aerolw(i,k,m,1) + tmp2*volcae(i)
1185                       enddo
1186                     enddo    ! end do_i_block
1187                   enddo      ! end do_m_block
1188     
1189                 endif      ! end if_NBDIR_block
1190               endif        ! end if_laddlw_block
1191     
1192             else                            ! input data from sfc to toa
1193     
1194               psrfh = 5.0                        ! ref press for upper bound
1195     
1196     !  ---  find lower boundary of stratosphere
1197     
1198               do i = 1, IMAX
1199     
1200                 tmp1 = abs( alat(i) )
1201                 if ( tmp1 > 70.0 ) then          ! polar, fixed at 250mb
1202                   psrfl = 250.0
1203                 elseif ( tmp1 < 20.0 ) then      ! tropic, fixed at 150mb
1204                   psrfl = 150.0
1205                 else                             ! mid-lat, interp
1206                   psrfl = 110.0 + 2.0*tmp1
1207                 endif
1208     
1209                 kcuth(i) = 2
1210                 kcutl(i) = NLAY - 1
1211                 delp(i) = prsi(i,NLAY-1)
1212     
1213                 lab_do_kcuth1 : do k = NLAY-1, 2, -1
1214                   if ( prsi(i,k) >= psrfh ) then
1215                     kcuth(i) = k
1216                     exit lab_do_kcuth1
1217                   endif
1218                 enddo  lab_do_kcuth1
1219     
1220                 lab_do_kcutl1 : do k = NLAY, 2, -1
1221                   if ( prsi(i,k) >= psrfl ) then
1222                     kcutl(i) = k
1223                     delp(i) = prsi(i,k) - prsi(i,kcuth(i)+1)
1224                     exit lab_do_kcutl1
1225                   endif
1226                 enddo  lab_do_kcutl1
1227               enddo
1228     
1229     !  ---  sw: add volcanic aerosol optical depth to the background value
1230     
1231               if ( laddsw ) then
1232                 do m = 1, NBDSW
1233                   mb = NSWSTR + m - 1
1234     
1235                   if     ( wvnum1(mb) > 20000 ) then   ! range of wvlth < 0.5mu
1236                     tmp2 = 0.74
1237                   elseif ( wvnum2(mb) < 20000 ) then   ! range of wvlth > 0.5mu
1238                     tmp2 = 1.14
1239                   else                                 ! range of wvlth in btwn
1240                     tmp2 = 0.94
1241                   endif
1242                   tmp1 = (0.275e-4 * (wvnum2(mb)+wvnum1(mb))) ** tmp2
1243     
1244                   do i = 1, IMAX
1245                     kh = kcuth(i)
1246                     kl = kcutl(i)
1247                     do k = kl, kh
1248                       tmp2 = tmp1 * ((prsi(i,k) - prsi(i,k+1)) / delp(i))
1249                       aerosw(i,k,m,1) = aerosw(i,k,m,1) + tmp2*volcae(i)
1250                     enddo
1251     
1252     !  ---  smoothing profile at boundary if needed
1253     
1254                     if ( aerosw(i,kl,m,1) > 10.*aerosw(i,kl-1,m,1) ) then
1255                       tmp2 = aerosw(i,kl,m,1) + aerosw(i,kl-1,m,1)
1256                       aerosw(i,kl  ,m,1) = 0.8 * tmp2
1257                       aerosw(i,kl-1,m,1) = 0.2 * tmp2
1258                     endif
1259                   enddo    ! end do_i_block
1260                 enddo      ! end do_m_block
1261     
1262     !  ---  check print
1263     
1264     !           do i = 1, IMAX
1265     !             print *,' LEV  PRESS    TAUSAV    NEWTAU    FOR PROFILE:',&
1266     !    &                i,'  KCUTH, KCUTL =',kcuth(i),kcutl(i)
1267     !             kh = kcuth(i) + 1
1268     !             kl = kcutl(i) - 10
1269     !             do k = kh, kl, -1
1270     !               write(6,71) NLP1-k,prsl(i,k),aersav(i,k),aerosw(i,k,1,1)
1271     !             enddo
1272     !           enddo
1273               endif        ! end if_laddsw_block
1274     
1275     !  ---  lw: add volcanic aerosol optical depth to the background value
1276     
1277               if ( laddlw ) then
1278                 if ( NBDIR == 1 ) then
1279     
1280                   tmp1 = (0.55 / 11.0) ** 1.2
1281                   do i = 1, IMAX
1282                     kh = kcuth(i)
1283                     kl = kcutl(i)
1284                     do k = kl, kh
1285                       tmp2 = tmp1 * ((prsi(i,k) - prsi(i,k+1)) / delp(i))
1286                       aerolw(i,k,1,1) = aerolw(i,k,1,1) + tmp2*volcae(i)
1287                     enddo
1288                   enddo    ! end do_i_block
1289     
1290                 else
1291     
1292                   do m = 1, NBDIR
1293                     tmp1 = (0.275e-4 * (wvnlw2(m) + wvnlw1(m))) ** 1.2
1294     
1295                     do i = 1, IMAX
1296                       kh = kcuth(i)
1297                       kl = kcutl(i)
1298                       do k = kl, kh
1299                         tmp2 = tmp1 * ((prsi(i,k)-prsi(i,k+1)) / delp(i))
1300                         aerolw(i,k,m,1) = aerolw(i,k,m,1) + tmp2*volcae(i)
1301                       enddo
1302                     enddo    ! end do_i_block
1303                   enddo      ! end do_m_block
1304     
1305                 endif      ! end if_NBDIR_block
1306               endif        ! end if_laddlw_block
1307     
1308             endif                           ! end if_iflip_block
1309     
1310     !  ---  adding volcanic optical depth to stratospheric layers
1311     
1312     
1313           endif   ! end if_lvolflg_block
1314     
1315     !
1316           return
1317     !...................................
1318           end subroutine setaer
1319     !-----------------------------------
1320     
1321     
1322     
1323     !-----------------------------------
1324           subroutine clim_aerinit                                           &
1325     !...................................
1326     !  ---  inputs:
1327          &     ( NWVTOT,solfwv,soltot,NWVTIR,eirfwv,                        &
1328          &       NBDSW,NBDIR,NBDSWLW, imon, me                              &
1329     !  ---  outputs: ( none )
1330          &     )
1331     
1332     !  ==================================================================  !
1333     !                                                                      !
1334     !  subprogram : clim_aerinit                                           !
1335     !                                                                      !
1336     !    this is the initialization progrmam for climatological aerosols   !
1337     !                                                                      !
1338     !    it reads in monthly global distribution of aerosol profiles in    !
1339     !    five degree horizontal resolution. Then, it reads and maps the    !
1340     !    tabulated aerosol optical spectral data onto corresponding sw     !
1341     !    radiation spectral bands.                                         !
1342     !                                                                      !
1343     !  ====================  defination of variables  ===================  !
1344     !                                                                      !
1345     !  inputs:                                                             !
1346     !   NWVTOT           - total num of wave numbers used in sw spectrum   !
1347     !   solfwv(NWVTOT)   - solar flux for each individual wavenumber (w/m2)!
1348     !   soltot           - total solar flux for the spectrual range  (w/m2)!
1349     !   NWVTIR           - total num of wave numbers used in the ir region !
1350     !   eirfwv(NWVTIR)   - ir flux(273k) for each individual wavenum (w/m2)!
1351     !   NBDSW            - num of bands calculated for sw aeros opt prop   !
1352     !   NBDIR            - num of bands calculated for lw aeros opt prop   !
1353     !   NBDSWLW          - total num of bands calc for sw+lw aeros opt prop!
1354     !   imon             - month of the year                               !
1355     !   me               - print message control flag                      !
1356     !                                                                      !
1357     !  outputs: (to the module variables)                                  !
1358     !                                                                      !
1359     !  module variables:                                                   !
1360     !     NBDSW   - total number of sw spectral bands                      !
1361     !     wvnum1,wvnum2 (NSWSTR:NSWEND)                                    !
1362     !             - start/end wavenumbers for each of sw bands             !
1363     !     NBDLW   - total number of lw spectral bands                      !
1364     !     wvnlw1,wvnlw2 (NBDLW)                                            !
1365     !             - start/end wavenumbers for each of lw bands             !
1366     !     NBDSWLW - total number of sw+lw bands used in this version       !
1367     !     extrhi  - extinction coef for rh-indep aeros         NCM1*NBDSWLW!
1368     !     scarhi  - scattering coef for rh-indep aeros         NCM1*NBDSWLW!
1369     !     ssarhi  - single-scat-alb for rh-indep aeros         NCM1*NBDSWLW!
1370     !     asyrhi  - asymmetry factor for rh-indep aeros        NCM1*NBDSWLW!
1371     !     extrhd  - extinction coef for rh-dep aeros    NRHLEV*NCM2*NBDSWLW!
1372     !     scarhd  - scattering coef for rh-dep aeros    NRHLEV*NCM2*NBDSWLW!
1373     !     ssarhd  - single-scat-alb for rh-dep aeros    NRHLEV*NCM2*NBDSWLW!
1374     !     asyrhd  - asymmetry factor for rh-dep aeros   NRHLEV*NCM2*NBDSWLW!
1375     !                                                                      !
1376     !     kprfg   - aerosols profile index                IMXAE*JMXAE      !
1377     !     idxcg   - aerosols component index              NXC*IMXAE*JMXAE  !
1378     !     cmixg   - aerosols component mixing ratio       NXC*IMXAE*JMXAE  !
1379     !     denng   - aerosols number density               NXC*IMXAE*JMXAE  !
1380     !                                                                      !
1381     !  major local variables:                                              !
1382     !   for handling spectral band structures                              !
1383     !     iendwv   - ending wvnum (cm**-1) for each band  NAERBND          !
1384     !   for handling optical properties of rh independent species (NCM1)   !
1385     !         1. insoluble        (inso); 2. soot             (soot);      !
1386     !         3. mineral nuc mode (minm); 4. mineral acc mode (miam);      !
1387     !         5. mineral coa mode (micm); 6. mineral transport(mitr).      !
1388     !     rhidext0 - extinction coefficient             NAERBND*NCM1       !
1389     !     rhidsca0 - scattering coefficient             NAERBND*NCM1       !
1390     !     rhidssa0 - single scattering albedo           NAERBND*NCM1       !
1391     !     rhidasy0 - asymmetry parameter                NAERBND*NCM1       !
1392     !   for handling optical properties of rh ndependent species (NCM2)    !
1393     !         1. water soluble    (waso); 2. sea salt acc mode(ssam);      !
1394     !         3. sea salt coa mode(sscm); 4. sulfate droplets (suso).      !
1395     !         rh level (NRHLEV): 00%, 50%, 70%, 80%, 90%, 95%, 98%, 99%    !
1396     !     rhdpext0 - extinction coefficient             NAERBND,NRHLEV,NCM2!
1397     !     rhdpsca0 - scattering coefficient             NAERBND,NRHLEV,NCM2!
1398     !     rhdpssa0 - single scattering albedo           NAERBND,NRHLEV,NCM2!
1399     !     rhdpasy0 - asymmetry parameter                NAERBND,NRHLEV,NCM2!
1400     !   for handling optical properties of stratospheric bkgrnd aerosols   !
1401     !     straext0 - extingction coefficients             NAERBND          !
1402     !                                                                      !
1403     !  usage:    call aerinit                                              !
1404     !                                                                      !
1405     !  subprograms called:  optavg                                         !
1406     !                                                                      !
1407     !  ==================================================================  !
1408     !
1409           implicit none
1410     
1411     !  ---  inputs:
1412           integer, intent(in) :: NWVTOT,NWVTIR,NBDSW,NBDIR,NBDSWLW,imon,me
1413     
1414           real (kind=kind_phys), intent(in) :: solfwv(:),soltot, eirfwv(:)
1415     
1416     !  ---  output: ( none )
1417     
1418     !  ---  locals:
1419           real (kind=kind_io8) :: cmix(NXC), denn, tem
1420           integer              :: idxc(NXC), kprf
1421     
1422           integer, dimension(NAERBND) :: iendwv
1423     
1424           real (kind=kind_phys), dimension(NAERBND,NCM1)       ::           &
1425          &       rhidext0, rhidsca0, rhidssa0, rhidasy0
1426           real (kind=kind_phys), dimension(NAERBND,NRHLEV,NCM2)::           &
1427          &       rhdpext0, rhdpsca0, rhdpssa0, rhdpasy0
1428           real (kind=kind_phys), dimension(NAERBND)            :: straext0
1429     
1430           real (kind=kind_phys), dimension(NBDSW,NAERBND) :: solwaer
1431           real (kind=kind_phys), dimension(NBDSW)         :: solbnd
1432           real (kind=kind_phys), dimension(NBDIR,NAERBND) :: eirwaer
1433           real (kind=kind_phys), dimension(NBDIR)         :: eirbnd
1434           real (kind=kind_phys) :: sumsol, sumir
1435     
1436           integer, dimension(NBDSW) :: nv1, nv2
1437           integer, dimension(NBDIR) :: nr1, nr2
1438     
1439           integer :: i, j, k, m, mb, nc, iy, ib, ii, id, iw, iw1, iw2
1440           logical :: file_exist
1441     
1442           character :: cline*80, ctyp*3, aerosol_file*24
1443     
1444     !     data aerosol_file / 'climaeropac_global.txt  ' /
1445           data aerosol_file / 'aerosol.dat  ' /
1446     
1447     !===>  ...  begin here
1448     
1449           if ( .not. allocated(kprfg) ) then
1450             allocate ( kprfg (    IMXAE,JMXAE) )
1451             allocate ( cmixg (NXC,IMXAE,JMXAE) )
1452             allocate ( denng (NXC,IMXAE,JMXAE) )
1453             allocate ( idxcg (NXC,IMXAE,JMXAE) )
1454           endif
1455     
1456     !  --- ...  reading climatological aerosols data
1457     
1458           inquire (file=aerosol_file, exist=file_exist)
1459     
1460           if ( file_exist ) then
1461             open (unit=NIAERCM,file=aerosol_file,status='OLD',              &
1462          &        form='FORMATTED')
1463             rewind (NIAERCM)
1464     
1465             if ( me == 0 ) then
1466               print *,'   Opened aerosol data file: ',aerosol_file
1467             endif
1468           else
1469             print *,'    Requested aerosol data file "',aerosol_file,       &
1470          &          '" not found!'
1471             print *,'    *** Stopped in subroutine AERINIT !!'
1472             stop
1473           endif              ! end if_file_exist_block
1474     
1475           cmixg = f_zero
1476           denng = f_zero
1477           idxcg = 0
1478     
1479     !  --- ...  loop over 12 month global distribution
1480     
1481           Lab_do_12mon : do m = 1, 12
1482     
1483             read(NIAERCM,12) cline
1484       12    format(a80/)
1485     
1486             if ( m /= imon ) then
1487     !         if ( me == 0 ) print *,'  *** Skipped ',cline
1488     
1489               do j = 1, JMXAE
1490               do i = 1, IMXAE
1491                 read(NIAERCM,*) id
1492               enddo
1493               enddo
1494             else
1495               if ( me == 0 ) print *,'  --- Reading ',cline
1496     
1497               do j = 1, JMXAE
1498               do i = 1, IMXAE
1499                 read(NIAERCM,14) (idxc(k),cmix(k),k=1,NXC),kprf,denn,nc,ctyp
1500       14        format(5(i2,e11.4),i2,f8.2,i3,1x,a3)
1501     
1502                 kprfg(i,j)     = kprf
1503                 denng(1,i,j)   = denn       ! num density of 1st layer
1504                 if ( kprf >= 6 ) then
1505                   denng(2,i,j) = cmix(NXC)  ! num density of 2dn layer
1506                 else
1507                   denng(2,i,j) = f_zero
1508                 endif
1509     
1510                 tem = f_one
1511                 do k = 1, NXC-1
1512                   idxcg(k,i,j) = idxc(k)    ! component index
1513                   cmixg(k,i,j) = cmix(k)    ! component mixing ratio
1514                   tem          = tem - cmix(k)
1515                 enddo
1516                 idxcg(NXC,i,j) = idxc(NXC)
1517                 cmixg(NXC,i,j) = tem        ! to make sure all add to 1.
1518               enddo
1519               enddo
1520     
1521               if ( .not. lclmin ) then
1522                 close (NIAERCM)
1523                 exit  Lab_do_12mon
1524               endif
1525             endif     ! end if_m_block
1526     
1527           enddo  Lab_do_12mon
1528     
1529     !  --  check print
1530     
1531     !     print *,'  IDXCG :'
1532     !     print 16,idxcg
1533     ! 16  format(40i3)
1534     !     print *,'  CMIXG :'
1535     !     print 17,cmixg
1536     !     print *,'  DENNG :'
1537     !     print 17,denng
1538     !     print *,'  KPRFG :'
1539     !     print 17,kprfg
1540     ! 17  format(8e16.9)
1541     
1542           if ( .not. lclmin ) then
1543     
1544     !  --- ...  already done optical property interpolation, exit
1545     
1546             return
1547     
1548           else
1549     
1550     !  --- ...  aloocate and input aerosol optical data
1551     
1552             if ( .not. allocated( haer ) ) then
1553               allocate ( haer   (NDM,NAE) )
1554               allocate ( prsref (NDM,NAE) )
1555             endif
1556     
1557             if ( .not. allocated( extrhi ) ) then
1558               allocate ( extrhi (       NCM1,NBDSWLW) )
1559               allocate ( scarhi (       NCM1,NBDSWLW) )
1560               allocate ( ssarhi (       NCM1,NBDSWLW) )
1561               allocate ( asyrhi (       NCM1,NBDSWLW) )
1562               allocate ( extrhd (NRHLEV,NCM2,NBDSWLW) )
1563               allocate ( scarhd (NRHLEV,NCM2,NBDSWLW) )
1564               allocate ( ssarhd (NRHLEV,NCM2,NBDSWLW) )
1565               allocate ( asyrhd (NRHLEV,NCM2,NBDSWLW) )
1566               allocate ( extstra(            NBDSWLW) )
1567             endif
1568     
1569             read(NIAERCM,21) cline   ! ending wave num for 61 aeros spectral bands
1570       21    format(a80)
1571             read(NIAERCM,22) iendwv(:)
1572       22    format(13i6)
1573     
1574             read(NIAERCM,21) cline   ! atmos scale height for 5 domains, 7 profs
1575             read(NIAERCM,24) haer(:,:)
1576       24    format(20f4.1)
1577     
1578             read(NIAERCM,21) cline   ! reference pressure for 5 domains, 7 profs
1579             read(NIAERCM,26) prsref(:,:)
1580       26    format(10f7.2)
1581     
1582             read(NIAERCM,21) cline   ! rh indep ext coef for 61 bands, 6 species
1583             read(NIAERCM,28) rhidext0(:,:)
1584       28    format(8e10.3)
1585     
1586             read(NIAERCM,21) cline   ! rh indep sca coef for 61 bands, 6 species
1587             read(NIAERCM,28) rhidsca0(:,:)
1588     
1589             read(NIAERCM,21) cline   ! rh indep ssa coef for 61 bands, 6 species
1590             read(NIAERCM,28) rhidssa0(:,:)
1591     
1592             read(NIAERCM,21) cline   ! rh indep asy coef for 61 bands, 6 species
1593             read(NIAERCM,28) rhidasy0(:,:)
1594     
1595             read(NIAERCM,21) cline   ! rh dep ext coef for 61 bands, 8 rh lev, 4 species
1596             read(NIAERCM,28) rhdpext0(:,:,:)
1597     
1598             read(NIAERCM,21) cline   ! rh dep sca coef for 61 bands, 8 rh lev, 4 species
1599             read(NIAERCM,28) rhdpsca0(:,:,:)
1600     
1601             read(NIAERCM,21) cline   ! rh dep ssa coef for 61 bands, 8 rh lev, 4 species
1602             read(NIAERCM,28) rhdpssa0(:,:,:)
1603     
1604             read(NIAERCM,21) cline   ! rh dep asy coef for 61 bands, 8 rh lev, 4 species
1605             read(NIAERCM,28) rhdpasy0(:,:,:)
1606     
1607             read(NIAERCM,21) cline   ! stratospheric background aeros for 61 bands
1608             read(NIAERCM,28) straext0(:)
1609     
1610             lclmin = .false.
1611     
1612     !  --- ...  compute solar flux weights and interval indices for mapping
1613     !           spectral bands between sw radiation and aerosol data
1614     
1615             solbnd (:)   = f_zero
1616             solwaer(:,:) = f_zero
1617     
1618             do ib = 1, NBDSW
1619               mb = ib + NSWSTR - 1
1620               ii = 1
1621               iw1 = nint(wvnum1(mb))
1622               iw2 = nint(wvnum2(mb))
1623     
1624     !
1625     ! ---  locate the spectral band for 550nm (for aod diag)
1626     !
1627               if (10000./iw1 >= 0.55 .and.
1628          &        10000./iw2 <=0.55 )  then
1629                   nv_aod =  ib
1630               endif
1631     
1632               Lab_swdowhile : do while ( iw1 > iendwv(ii) )
1633                 if ( ii == NAERBND ) exit Lab_swdowhile
1634                 ii = ii + 1
1635               enddo  Lab_swdowhile
1636     
1637               sumsol = f_zero
1638               nv1(ib) = ii
1639     
1640               do iw = iw1, iw2
1641                 solbnd(ib) = solbnd(ib) + solfwv(iw)
1642                 sumsol = sumsol + solfwv(iw)
1643     
1644                 if ( iw == iendwv(ii) ) then
1645                   solwaer(ib,ii) = sumsol
1646     
1647                   if ( ii < NAERBND ) then
1648                     sumsol = f_zero
1649                     ii = ii + 1
1650                   endif
1651                 endif
1652               enddo
1653     
1654               if ( iw2 /= iendwv(ii) ) then
1655                 solwaer(ib,ii) = sumsol
1656               endif
1657     
1658               nv2(ib) = ii
1659     !         frcbnd(ib) = solbnd(ib) / soltot
1660             enddo     ! end do_ib_block for sw
1661     
1662     !  --- ...  compute ir flux weights and interval indices for mapping
1663     !           spectral bands between lw radiation and aerosol data
1664     
1665             eirbnd (:)   = f_zero
1666             eirwaer(:,:) = f_zero
1667     
1668             do ib = 1, NBDIR
1669               ii = 1
1670               if ( NBDIR == 1 ) then
1671     !           iw1 = 250                   ! corresponding 40 mu
1672                 iw1 = 400                   ! corresponding 25 mu
1673                 iw2 = 2500                  ! corresponding 4  mu
1674               else
1675                 iw1 = nint(wvnlw1(ib))
1676                 iw2 = nint(wvnlw2(ib))
1677               endif
1678     
1679               Lab_lwdowhile : do while ( iw1 > iendwv(ii) )
1680                 if ( ii == NAERBND ) exit Lab_lwdowhile
1681                 ii = ii + 1
1682               enddo  Lab_lwdowhile
1683     
1684               sumir = f_zero
1685               nr1(ib) = ii
1686     
1687               do iw = iw1, iw2
1688                 eirbnd(ib) = eirbnd(ib) + eirfwv(iw)
1689                 sumir  = sumir  + eirfwv(iw)
1690     
1691                 if ( iw == iendwv(ii) ) then
1692                   eirwaer(ib,ii) = sumir
1693     
1694                   if ( ii < NAERBND ) then
1695                     sumir = f_zero
1696                     ii = ii + 1
1697                   endif
1698                 endif
1699               enddo
1700     
1701               if ( iw2 /= iendwv(ii) ) then
1702                 eirwaer(ib,ii) = sumir
1703               endif
1704     
1705               nr2(ib) = ii
1706             enddo     ! end do_ib_block for lw
1707     
1708     !  ---  compute spectral band mean properties for each species
1709     
1710             call optavg
1711     !  ---  inputs:  (in scope variables)
1712     !  ---  outputs: (in scope variables)
1713     
1714     !  ---  check print
1715     
1716     !       do ib = 1, NBDSW
1717     !         print *,' After optavg, for sw band:',ib
1718     !         print *,'  extrhi:', extrhi(:,ib)
1719     !         print *,'  scarhi:', scarhi(:,ib)
1720     !         print *,'  ssarhi:', ssarhi(:,ib)
1721     !         print *,'  asyrhi:', asyrhi(:,ib)
1722     !         mb = ib + NSWSTR - 1
1723     !         print *,'  wvnum1,wvnum2 :',wvnum1(mb),wvnum2(mb)
1724     !         do i = 1, NRHLEV
1725     !           print *,'  extrhd for rhlev:',i
1726     !           print *,extrhd(i,:,ib)
1727     !           print *,'  scarhd for rhlev:',i
1728     !           print *,scarhd(i,:,ib)
1729     !           print *,'  ssarhd for rhlev:',i
1730     !           print *,ssarhd(i,:,ib)
1731     !           print *,'  asyrhd for rhlev:',i
1732     !           print *,asyrhd(i,:,ib)
1733     !         enddo
1734     !         print *,' extstra:', extstra(ib)
1735     !       enddo
1736     !       print *,'  wvnlw1 :',wvnlw1
1737     !       print *,'  wvnlw2 :',wvnlw2
1738     !       do ib = 1, NBDIR
1739     !         ii = NBDSW + ib
1740     !         print *,' After optavg, for lw band:',ib
1741     !         print *,'  extrhi:', extrhi(:,ii)
1742     !         print *,'  scarhi:', scarhi(:,ii)
1743     !         print *,'  ssarhi:', ssarhi(:,ii)
1744     !         print *,'  asyrhi:', asyrhi(:,ii)
1745     !         do i = 1, NRHLEV
1746     !           print *,'  extrhd for rhlev:',i
1747     !           print *,extrhd(i,:,ii)
1748     !           print *,'  scarhd for rhlev:',i
1749     !           print *,scarhd(i,:,ii)
1750     !           print *,'  ssarhd for rhlev:',i
1751     !           print *,ssarhd(i,:,ii)
1752     !           print *,'  asyrhd for rhlev:',i
1753     !           print *,asyrhd(i,:,ii)
1754     !         enddo
1755     !         print *,' extstra:', extstra(ii)
1756     !       enddo
1757     
1758           endif  ! end if_lclmin_block
1759     
1760     ! =================
1761           contains
1762     ! =================
1763     
1764     !-----------------------------
1765           subroutine optavg
1766     !.............................
1767     !  ---  inputs:  (in scope variables)
1768     !  ---  outputs: (in scope variables)
1769     
1770     ! ==================================================================== !
1771     !                                                                      !
1772     ! subprogram: optavg                                                   !
1773     !                                                                      !
1774     !   compute mean aerosols optical properties over each sw radiation    !
1775     !   spectral band for each of the species components.  This program    !
1776     !   follows gfdl's approach for thick cloud opertical property in      !
1777     !   sw radiation scheme (2000).                                        !
1778     !                                                                      !
1779     !  ====================  defination of variables  ===================  !
1780     !                                                                      !
1781     ! input arguments:                                                     !
1782     !   nv1,nv2 (NBDSW)  - start/end spectral band indices of aerosol data !
1783     !                      for each sw radiation spectral band             !
1784     !   nr1,nr2 (NBDIR)  - start/end spectral band indices of aerosol data !
1785     !                      for each ir radiation spectral band             !
1786     !   solwaer (NBDSW,NAERBND)                                            !
1787     !                    - solar flux weight over each sw radiation band   !
1788     !                      vs each aerosol data spectral band              !
1789     !   eirwaer (NBDIR,NAERBND)                                            !
1790     !                    - ir flux weight over each lw radiation band      !
1791     !                      vs each aerosol data spectral band              !
1792     !   solbnd  (NBDSW)  - solar flux weight over each sw radiation band   !
1793     !   eirbnd  (NBDIR)  - ir flux weight over each lw radiation band      !
1794     !   NBDSW            - total number of sw spectral bands               !
1795     !   NBDIR            - total number of lw spectral bands               !
1796     !   NBDSWLW          - total number of sw+lw spectral bands            !
1797     !                                                                      !
1798     ! output arguments: (to module variables)                              !
1799     !                                                                      !
1800     !  ==================================================================  !
1801     !
1802           implicit none
1803     
1804     !  ---  inputs:
1805     !  ---  output:
1806     
1807     !  ---  locals:
1808           real (kind=kind_phys) :: sumk, sums, sumok, sumokg, sumreft,      &
1809          &       sp, refb, reft, rsolbd, rirbd
1810     
1811           integer :: ib, nb, ni, nh, nc
1812     !
1813     !===> ...  begin here
1814     !
1815     !  --- ...  loop for each sw radiation spectral band
1816     
1817           do nb = 1, NBDSW
1818             rsolbd = f_one / solbnd(nb)
1819     
1820     !  ---  for rh independent aerosol species
1821     
1822             do nc = 1, NCM1
1823               sumk    = f_zero
1824               sums    = f_zero
1825               sumok   = f_zero
1826               sumokg  = f_zero
1827               sumreft = f_zero
1828     
1829               do ni = nv1(nb), nv2(nb)
1830                 sp   = sqrt( (f_one - rhidssa0(ni,nc))                      &
1831          &           / (f_one - rhidssa0(ni,nc)*rhidasy0(ni,nc)) )
1832                 reft = (f_one - sp) / (f_one + sp)
1833                 sumreft = sumreft + reft*solwaer(nb,ni)
1834     
1835                 sumk    = sumk    + rhidext0(ni,nc)*solwaer(nb,ni)
1836                 sums    = sums    + rhidsca0(ni,nc)*solwaer(nb,ni)
1837                 sumok   = sumok   + rhidssa0(ni,nc)*solwaer(nb,ni)          &
1838          &              * rhidext0(ni,nc)
1839                 sumokg  = sumokg  + rhidssa0(ni,nc)*solwaer(nb,ni)          &
1840          &              * rhidext0(ni,nc)*rhidasy0(ni,nc)
1841               enddo
1842     
1843               refb = sumreft * rsolbd
1844     
1845               extrhi(nc,nb) = sumk   * rsolbd
1846               scarhi(nc,nb) = sums   * rsolbd
1847               asyrhi(nc,nb) = sumokg / (sumok + 1.0e-10)
1848               ssarhi(nc,nb) = 4.0*refb                                      &
1849          &         / ( (f_one+refb)**2 - asyrhi(nc,nb)*(f_one-refb)**2 )
1850             enddo   ! end do_nc_block for rh-ind aeros
1851     
1852     !  ---  for rh dependent aerosols species
1853     
1854             do nc = 1, NCM2
1855               do nh = 1, NRHLEV
1856                 sumk    = f_zero
1857                 sums    = f_zero
1858                 sumok   = f_zero
1859                 sumokg  = f_zero
1860                 sumreft = f_zero
1861     
1862                 do ni = nv1(nb), nv2(nb)
1863                   sp   = sqrt( (f_one - rhdpssa0(ni,nh,nc))                 &
1864          &             / (f_one - rhdpssa0(ni,nh,nc)*rhdpasy0(ni,nh,nc)) )
1865                   reft = (f_one - sp) / (f_one + sp)
1866                   sumreft = sumreft + reft*solwaer(nb,ni)
1867     
1868                   sumk    = sumk    + rhdpext0(ni,nh,nc)*solwaer(nb,ni)
1869                   sums    = sums    + rhdpsca0(ni,nh,nc)*solwaer(nb,ni)
1870                   sumok   = sumok   + rhdpssa0(ni,nh,nc)*solwaer(nb,ni)     &
1871          &                * rhdpext0(ni,nh,nc)
1872                   sumokg  = sumokg  + rhdpssa0(ni,nh,nc)*solwaer(nb,ni)     &
1873          &                * rhdpext0(ni,nh,nc)*rhdpasy0(ni,nh,nc)
1874                 enddo
1875     
1876                 refb = sumreft * rsolbd
1877     
1878                 extrhd(nh,nc,nb) = sumk   * rsolbd
1879                 scarhd(nh,nc,nb) = sums   * rsolbd
1880                 asyrhd(nh,nc,nb) = sumokg / (sumok + 1.0e-10)
1881                 ssarhd(nh,nc,nb) = 4.0*refb                                 &
1882          &         / ( (f_one+refb)**2 - asyrhd(nh,nc,nb)*(f_one-refb)**2 )
1883               enddo   ! end do_nh_block
1884             enddo   ! end do_nc_block for rh-dep aeros
1885     
1886     !  ---  for stratospheric background aerosols
1887     
1888             sumk = f_zero
1889             do ni = nv1(nb), nv2(nb)
1890               sumk = sumk + straext0(ni)*solwaer(nb,ni)
1891             enddo
1892     
1893             extstra(nb) = sumk * rsolbd
1894     
1895     !  ---  check print
1896     !       if ( nb > 6 .and. nb < 10) then
1897     !         print *,' in optavg for sw band',nb
1898     !         print *,'  nv1, nv2:',nv1(nb),nv2(nb)
1899     !         print *,'  solwaer:',solwaer(nb,nv1(nb):nv2(nb))
1900     !         print *,'  extrhi:', extrhi(:,nb)
1901     !         do i = 1, NRHLEV
1902     !           print *,'  extrhd for rhlev:',i
1903     !           print *,extrhd(i,:,nb)
1904     !         enddo
1905     !         print *,'  sumk, rsolbd, extstra:',sumk,rsolbd,extstra(nb)
1906     !       endif
1907     
1908           enddo   !  end do_nb_block for sw
1909     
1910     !  --- ...  loop for each lw radiation spectral band
1911     
1912           do nb = 1, NBDIR
1913     
1914             ib = NBDSW + nb
1915             rirbd = f_one / eirbnd(nb)
1916     
1917     !  ---  for rh independent aerosol species
1918     
1919             do nc = 1, NCM1
1920               sumk    = f_zero
1921               sums    = f_zero
1922               sumok   = f_zero
1923               sumokg  = f_zero
1924               sumreft = f_zero
1925     
1926               do ni = nr1(nb), nr2(nb)
1927                 sp   = sqrt( (f_one - rhidssa0(ni,nc))                      &
1928          &           / (f_one - rhidssa0(ni,nc)*rhidasy0(ni,nc)) )
1929                 reft = (f_one - sp) / (f_one + sp)
1930                 sumreft = sumreft + reft*eirwaer(nb,ni)
1931     
1932                 sumk    = sumk    + rhidext0(ni,nc)*eirwaer(nb,ni)
1933                 sums    = sums    + rhidsca0(ni,nc)*eirwaer(nb,ni)
1934                 sumok   = sumok   + rhidssa0(ni,nc)*eirwaer(nb,ni)          &
1935          &              * rhidext0(ni,nc)
1936                 sumokg  = sumokg  + rhidssa0(ni,nc)*eirwaer(nb,ni)          &
1937          &              * rhidext0(ni,nc)*rhidasy0(ni,nc)
1938               enddo
1939     
1940               refb = sumreft * rirbd
1941     
1942               extrhi(nc,ib) = sumk   * rirbd
1943               scarhi(nc,ib) = sums   * rirbd
1944               asyrhi(nc,ib) = sumokg / (sumok + 1.0e-10)
1945               ssarhi(nc,ib) = 4.0*refb                                         &
1946          &         / ( (f_one+refb)**2 - asyrhi(nc,ib)*(f_one-refb)**2 )
1947             enddo   ! end do_nc_block for rh-ind aeros
1948     
1949     !  ---  for rh dependent aerosols species
1950     
1951             do nc = 1, NCM2
1952               do nh = 1, NRHLEV
1953                 sumk    = f_zero
1954                 sums    = f_zero
1955                 sumok   = f_zero
1956                 sumokg  = f_zero
1957                 sumreft = f_zero
1958     
1959                 do ni = nr1(nb), nr2(nb)
1960                   sp   = sqrt( (f_one - rhdpssa0(ni,nh,nc))                 &
1961          &           / (f_one - rhdpssa0(ni,nh,nc)*rhdpasy0(ni,nh,nc)) )
1962                   reft = (f_one - sp) / (f_one + sp)
1963                   sumreft = sumreft + reft*eirwaer(nb,ni)
1964     
1965                   sumk    = sumk    + rhdpext0(ni,nh,nc)*eirwaer(nb,ni)
1966                   sums    = sums    + rhdpsca0(ni,nh,nc)*eirwaer(nb,ni)
1967                   sumok   = sumok   + rhdpssa0(ni,nh,nc)*eirwaer(nb,ni)     &
1968          &                * rhdpext0(ni,nh,nc)
1969                   sumokg  = sumokg  + rhdpssa0(ni,nh,nc)*eirwaer(nb,ni)     &
1970          &                * rhdpext0(ni,nh,nc)*rhdpasy0(ni,nh,nc)
1971                 enddo
1972     
1973                 refb = sumreft * rirbd
1974     
1975                 extrhd(nh,nc,ib) = sumk   * rirbd
1976                 scarhd(nh,nc,ib) = sums   * rirbd
1977                 asyrhd(nh,nc,ib) = sumokg / (sumok + 1.0e-10)
1978                 ssarhd(nh,nc,ib) = 4.0*refb                                    &
1979          &         / ( (f_one+refb)**2 - asyrhd(nh,nc,ib)*(f_one-refb)**2 )
1980               enddo   ! end do_nh_block
1981             enddo   ! end do_nc_block for rh-dep aeros
1982     
1983     !  ---  for stratospheric background aerosols
1984     
1985             sumk = f_zero
1986             do ni = nr1(nb), nr2(nb)
1987               sumk = sumk + straext0(ni)*eirwaer(nb,ni)
1988             enddo
1989     
1990             extstra(ib) = sumk * rirbd
1991     
1992     !  ---  check print
1993     !       if ( nb >= 1 .and. nb < 5) then
1994     !         print *,' in optavg for ir band:',nb
1995     !         print *,'  nr1, nr2:',nr1(nb),nr2(nb)
1996     !         print *,'  eirwaer:',eirwaer(nb,nr1(nb):nr2(nb))
1997     !         print *,'  extrhi:', extrhi(:,ib)
1998     !         do i = 1, NRHLEV
1999     !           print *,'  extrhd for rhlev:',i
2000     !           print *,extrhd(i,:,ib)
2001     !         enddo
2002     !         print *,'  sumk, rirbd, extstra:',sumk,rirbd,extstra(ib)
2003     !       endif
2004     
2005           enddo   !  end do_nb_block for lw
2006     
2007     !
2008           return
2009     !................................
2010           end subroutine optavg
2011     !--------------------------------
2012     !
2013     !...................................
2014           end subroutine clim_aerinit
2015     !-----------------------------------
2016     
2017     
2018     !-----------------------------------
2019           subroutine setclimaer                                             &
2020     !...................................
2021     
2022     !  ---  inputs:
2023          &     ( alon,alat,prsi,rhlay,dz,hz,NBDSWLW,                        &
2024          &       IMAX,NLAY,NLP1, iflip, lsswr,lslwr,                        &
2025     !  ---  outputs:
2026          &       aerosw,aerolw                                              &
2027          &     )
2028     
2029     !  ==================================================================  !
2030     !                                                                      !
2031     !  setaer maps the 5 degree global climatological aerosol data set     !
2032     !  onto model grids                                                    !
2033     !                                                                      !
2034     !  inputs:                                                             !
2035     !     NBDSW   - total number of sw spectral bands               1      !
2036     !     alon, alat                                             IMAX      !
2037     !             - longitude and latitude of given points in degree       !
2038     !     prsi    - pressure at interface              mb   IMAX*NLP1      !
2039     !     rhlay   - layer mean relative humidity            IMAX*NLAY      !
2040     !     dz      - layer thickness                    m    IMAX*NLAY      !
2041     !     hz      - level high                         m    IMAX*NLP1      !
2042     !     NBDSWLW - total number of sw+ir bands for aeros opt prop  1      !
2043     !     IMAX    - horizontal dimension of arrays                  1      !
2044     !     NLAY,NLP1-vertical dimensions of arrays                   1      !
2045     !     iflip   - control flag for direction of vertical index    1      !
2046     !               =0: index from toa to surface                          !
2047     !               =1: index from surface to toa                          !
2048     !     lsswr,lslwr                                                      !
2049     !             - logical flag for sw/lw radiation calls          1      !
2050     !                                                                      !
2051     !  outputs:                                                            !
2052     !     aerosw - aeros opt properties for sw      IMAX*NLAY*NBDSW*NF_AESW!
2053     !               (:,:,:,1): optical depth                               !
2054     !               (:,:,:,2): single scattering albedo                    !
2055     !               (:,:,:,3): asymmetry parameter                         !
2056     !     aerolw - aeros opt properties for lw      IMAX*NLAY*NBDLW*NF_AELW!
2057     !               (:,:,:,1): optical depth                               !
2058     !               (:,:,:,2): single scattering albedo                    !
2059     !               (:,:,:,3): asymmetry parameter                         !
2060     !                                                                      !
2061     !  module parameters and constants:                                    !
2062     !     NBDSW   - total number of sw bands for aeros opt prop     1      !
2063     !     NBDIR   - total number of ir bands for aeros opt prop     1      !
2064     !                                                                      !
2065     !  module variable: (set by subroutine clim_aerinit)                   !
2066     !     kprfg   - aerosols profile index                IMXAE*JMXAE      !
2067     !     idxcg   - aerosols component index              NXC*IMXAE*JMXAE  !
2068     !     cmixg   - aerosols component mixing ratio       NXC*IMXAE*JMXAE  !
2069     !     denng   - aerosols number density               NXC*IMXAE*JMXAE  !
2070     !                                                                      !
2071     !  usage:    call setclimaer                                           !
2072     !                                                                      !
2073     !  subprograms called:  radclimaer                                     !
2074     !                                                                      !
2075     !  ==================================================================  !
2076     !
2077           implicit none
2078     
2079     !  ---  inputs:
2080           integer, intent(in) :: IMAX,NLAY,NLP1,iflip,NBDSWLW
2081           logical, intent(in) :: lsswr, lslwr
2082     
2083           real (kind=kind_phys), dimension(:,:), intent(in) :: prsi,        &
2084          &       rhlay, dz, hz
2085           real (kind=kind_phys), dimension(:),   intent(in) :: alon, alat
2086     
2087     !  ---  outputs:
2088           real (kind=kind_phys), dimension(:,:,:,:), intent(out) ::         &
2089          &       aerosw, aerolw
2090     
2091     !  ---  locals:
2092           real (kind=kind_phys), dimension(NXC)  :: cmix, denn
2093           integer,               dimension(NXC)  :: idxc
2094     
2095           real (kind=kind_phys), dimension(NLAY) :: delz, rh1, dz1
2096           integer,               dimension(NLAY) :: idmaer
2097     
2098           real (kind=kind_phys), dimension(NLAY,NBDSWLW):: tauae,ssaae,asyae
2099     !test real (kind=kind_phys), dimension(IMAX,NLAY) :: aersav
2100     
2101           real (kind=kind_phys) :: tmp1, tmp2
2102     
2103           integer               :: i, i1, i2, j1, j2, k, m, m1, kp
2104     
2105     !  ---  conversion constants
2106           real (kind=kind_phys), parameter :: dltg = 360.0 / float(IMXAE)
2107           real (kind=kind_phys), parameter :: hdlt = 0.5 * dltg
2108     
2109     !
2110     !===>  ...  begin here
2111     !
2112     !  ---  map grid in longitude direction
2113     
2114           lab_do_IMAX : do i = 1, IMAX
2115             i2 = 1
2116             j2 = 1
2117     
2118             lab_do_IMXAE : do i1 = 1, IMXAE
2119               tmp1 = dltg * (i1 - 1) + hdlt
2120     
2121               if (abs(alon(i)-tmp1) <= hdlt) then
2122                 i2 = i1
2123                 exit lab_do_IMXAE
2124               endif
2125             enddo  lab_do_IMXAE
2126     
2127     !  ---  map grid in latitude direction
2128     
2129             lab_do_JMXAE : do j1 = 1, JMXAE
2130               tmp2 = 90.0 - dltg * (j1 - 1)
2131     
2132               if (abs(alat(i)-tmp2) <= hdlt) then
2133                 j2 = j1
2134                 exit lab_do_JMXAE
2135               endif
2136             enddo  lab_do_JMXAE
2137     
2138             do m = 1, NXC
2139               idxc(m) = idxcg(m,i2,j2)
2140               cmix(m) = cmixg(m,i2,j2)
2141               denn(m) = denng(m,i2,j2)
2142             enddo
2143     
2144             kp = kprfg(i2,j2)
2145     
2146             do k = 1, NLAY
2147               rh1(k) = rhlay(i,k)
2148               dz1(k) = dz   (i,k)
2149             enddo
2150     
2151     !  ---  compute vertical domain indices
2152     
2153             lab_if_flip : if (iflip == 1) then        ! input from sfc to toa
2154     
2155     !  ---  setup domain index array and effective layer thickness
2156     
2157               i1 = 1
2158               do k = 1, NLAY
2159                 if (prsi(i,k+1) < prsref(i1,kp)) then
2160                   i1 = i1 + 1
2161                   if (i1 == 2 .and. prsref(2,kp) == prsref(3,kp)) then
2162                     i1 = 3
2163                   endif
2164                 endif
2165                 idmaer(k) = i1
2166     
2167                 tmp1 = haer(i1,kp)
2168                 if (tmp1 > f_zero) then
2169                   tmp2 = f_one / tmp1
2170                   delz(k) = tmp1 * (exp(-hz(i,k)*tmp2)-exp(-hz(i,k+1)*tmp2))
2171                 else
2172                   delz(k) = dz1(k)
2173                 endif
2174               enddo
2175     
2176             else  lab_if_flip                         ! input from toa to sfc
2177     
2178     !  ---  setup domain index array and modified layer thickness
2179     
2180               i1 = 1
2181               do k = NLAY, 1, -1
2182                 if (prsi(i,k) < prsref(i1,kp)) then
2183                   i1 = i1 + 1
2184                   if (i1 == 2 .and. prsref(2,kp) == prsref(3,kp)) then
2185                     i1 = 3
2186                   endif
2187                 endif
2188                 idmaer(k) = i1
2189     
2190                 tmp1 = haer(i1,kp)
2191                 if (tmp1 > f_zero) then
2192                   tmp2   = f_one / tmp1
2193                   delz(k) = tmp1 * (exp(-hz(i,k+1)*tmp2)-exp(-hz(i,k)*tmp2))
2194                 else
2195                   delz(k) = dz1(k)
2196                 endif
2197               enddo
2198     
2199             endif  lab_if_flip
2200     
2201     !  ---  check print
2202     
2203     !       print *,' in setclimaer, profile:',i
2204     !       print *,'  rh   :',rh1
2205     !       print *,'  dz   :',dz1
2206     !       print *,'  delz :',delz
2207     !       print *,'  idmaer:',idmaer
2208     
2209     !  ---  calculate sw/lw aerosol optical properties for the
2210     !       corresponding frequency bands
2211     
2212             call radclimaer
2213     !  ---  inputs:  (in scope variables)
2214     !  ---  outputs: (in scope variables)
2215     
2216             if ( lsswr ) then
2217     
2218               if ( laswflg ) then
2219     
2220                 do m = 1, NBDSW
2221                   do k = 1, NLAY
2222                     aerosw(i,k,m,1) = tauae(k,m)
2223                     aerosw(i,k,m,2) = ssaae(k,m)
2224                     aerosw(i,k,m,3) = asyae(k,m)
2225                   enddo
2226                 enddo
2227     
2228               else
2229     
2230                 aerosw(:,:,:,:) = f_zero
2231     
2232               endif
2233     
2234     !  ---  testing use
2235     !         do k = 1, NLAY
2236     !           aersav(i,k) = tauae(k,1)
2237     !test     enddo
2238             endif     ! end if_lsswr_block
2239     
2240             if ( lslwr ) then
2241     
2242               if ( lalwflg ) then
2243     
2244                 if ( NBDIR == 1 ) then
2245                   m1 = NBDSW + 1
2246                   do m = 1, NBDLW
2247                     do k = 1, NLAY
2248                       aerolw(i,k,m,1) = tauae(k,m1)
2249                       aerolw(i,k,m,2) = ssaae(k,m1)
2250                       aerolw(i,k,m,3) = asyae(k,m1)
2251                     enddo
2252                   enddo
2253                 else
2254                   do m = 1, NBDLW
2255                     m1 = NBDSW + m
2256                     do k = 1, NLAY
2257                       aerolw(i,k,m,1) = tauae(k,m1)
2258                       aerolw(i,k,m,2) = ssaae(k,m1)
2259                       aerolw(i,k,m,3) = asyae(k,m1)
2260                     enddo
2261                   enddo
2262                 endif
2263     
2264               else
2265     
2266                 aerolw(:,:,:,:) = f_zero
2267     
2268               endif
2269             endif     ! end if_lslwr_block
2270     
2271           enddo  lab_do_IMAX
2272     
2273     ! =================
2274           contains
2275     ! =================
2276     
2277     !--------------------------------
2278           subroutine radclimaer
2279     !................................
2280     
2281     !  ---  inputs:  (in scope variables)
2282     !  ---  outputs: (in scope variables)
2283     
2284     !  ==================================================================  !
2285     !                                                                      !
2286     !  compute aerosols optical properties in NBDSW sw bands. there are    !
2287     !  seven different vertical profile structures. in the troposphere,    !
2288     !  aerosol distribution at each grid point is composed from up to      !
2289     !  six components out of a total of ten different substances.          !
2290     !                                                                      !
2291     !  ref: wmo report wcp-112 (1986)                                      !
2292     !                                                                      !
2293     !  input variables:                                                    !
2294     !     idxc   - indices of aerosol components         -     NXC         !
2295     !     cmix   - mixing ratioes of aerosol components  -     NXC         !
2296     !     denn   - aerosol number densities              -     NXC         !
2297     !     rh1    - relative humidity                     -     NLAY        !
2298     !     delz   - effective layer thickness             km    NLAY        !
2299     !     idmaer - aerosol domain index                  -     NLAY        !
2300     !     NXC    - number of different aerosol components-     1           !
2301     !     NLAY   - vertical dimensions                   -     1           !
2302     !     iflip  - control flag for direction of vertical index            !
2303     !               =0: index from toa to surface                          !
2304     !               =1: index from surface to toa                          !
2305     !                                                                      !
2306     !  output variables:                                                   !
2307     !     tauae  - optical depth                         -     NLAY*NBDSW  !
2308     !     ssaae  - single scattering albedo              -     NLAY*NBDSW  !
2309     !     asyae  - asymmetry parameter                   -     NLAY*NBDSW  !
2310     !                                                                      !
2311     !  ==================================================================  !
2312     !
2313           implicit none
2314     
2315     !
2316           real (kind=kind_phys) :: crt1, crt2
2317           parameter (crt1=30.0, crt2=0.03333)
2318     
2319     !  ---  inputs:
2320     !  ---  outputs:
2321     
2322     !  ---  locals:
2323           real (kind=kind_phys) :: cm, hd, hdi, sig0u, sig0l, ratio, tt0,   &
2324          &      ex00, sc00, ss00, as00, ex01, sc01, ss01, as01,     tt1,    &
2325          &      ex02, sc02, ss02, as02, ex03, sc03, ss03, as03,     tt2,    &
2326          &      ext1, sca1, ssa1, asy1, drh0, drh1, rdrh
2327     
2328           integer :: ih1, ih2, kk, idom, icmp, ib, ii, ic, ic1
2329     
2330     !
2331     !===> ... loop over vertical layers from top to surface
2332     !
2333           lab_do_layer : do kk = 1, NLAY
2334     
2335     ! --- linear interp coeffs for rh-dep species
2336     
2337             ih2 = 1
2338             do while ( rh1(kk) > rhlev(ih2) )
2339               ih2 = ih2 + 1
2340               if ( ih2 > NRHLEV ) exit
2341             enddo
2342             ih1 = max( 1, ih2-1 )
2343             ih2 = min( NRHLEV, ih2 )
2344     
2345             drh0 = rhlev(ih2) - rhlev(ih1)
2346             drh1 = rh1(kk) - rhlev(ih1)
2347             if ( ih1 == ih2 ) then
2348               rdrh = f_zero
2349             else
2350               rdrh = drh1 / drh0
2351             endif
2352     
2353     ! --- assign optical properties in each domain
2354     
2355             idom = idmaer(kk)
2356     
2357             lab_if_idom : if (idom == 5) then
2358     ! --- 5th domain - upper stratosphere assume no aerosol
2359     
2360               do ib = 1, NBDSWLW
2361                 tauae(kk,ib) = f_zero
2362                 if ( ib <= NBDSW ) then
2363                   ssaae(kk,ib) = 0.99
2364                   asyae(kk,ib) = 0.696
2365                 else
2366                   ssaae(kk,ib) = 0.5
2367                   asyae(kk,ib) = 0.3
2368                 endif
2369               enddo
2370     
2371             elseif (idom == 4) then    lab_if_idom
2372     ! --- 4th domain - stratospheric layers
2373     
2374               do ib = 1, NBDSWLW
2375                 tauae(kk,ib) = extstra(ib) * delz(kk)
2376                 if ( ib <= NBDSW ) then
2377                   ssaae(kk,ib) = 0.99
2378                   asyae(kk,ib) = 0.696
2379                 else
2380                   ssaae(kk,ib) = 0.5
2381                   asyae(kk,ib) = 0.3
2382                 endif
2383               enddo
2384     
2385             elseif (idom == 3) then    lab_if_idom
2386     ! --- 3rd domain - free tropospheric layers
2387     !   1:inso 0.17e-3; 2:soot 0.4; 7:waso 0.59983; n:730
2388     
2389               do ib = 1, NBDSWLW
2390                 ex01 = extrhi(1,ib)
2391                 sc01 = scarhi(1,ib)
2392                 ss01 = ssarhi(1,ib)
2393                 as01 = asyrhi(1,ib)
2394     
2395                 ex02 = extrhi(2,ib)
2396                 sc02 = scarhi(2,ib)
2397                 ss02 = ssarhi(2,ib)
2398                 as02 = asyrhi(2,ib)
2399     
2400                 ex03 = extrhd(ih1,1,ib)                                     &
2401          &           + rdrh * (extrhd(ih2,1,ib) - extrhd(ih1,1,ib))
2402                 sc03 = scarhd(ih1,1,ib)                                     &
2403          &           + rdrh * (scarhd(ih2,1,ib) - scarhd(ih1,1,ib))
2404                 ss03 = ssarhd(ih1,1,ib)                                     &
2405          &           + rdrh * (ssarhd(ih2,1,ib) - ssarhd(ih1,1,ib))
2406                 as03 = asyrhd(ih1,1,ib)                                     &
2407          &           + rdrh * (asyrhd(ih2,1,ib) - asyrhd(ih1,1,ib))
2408     
2409                 ext1 = 0.17e-3*ex01 + 0.4*ex02 + 0.59983*ex03
2410                 sca1 = 0.17e-3*sc01 + 0.4*sc02 + 0.59983*sc03
2411                 ssa1 = 0.17e-3*ss01*ex01 + 0.4*ss02*ex02 + 0.59983*ss03*ex03
2412                 asy1 = 0.17e-3*as01*sc01 + 0.4*as02*sc02 + 0.59983*as03*sc03
2413     
2414                 tauae(kk,ib) = ext1 * 730.0 * delz(kk)
2415                 ssaae(kk,ib) = min(f_one, ssa1/ext1)
2416                 asyae(kk,ib) = min(f_one, asy1/sca1)
2417               enddo
2418     
2419             elseif (idom == 1) then    lab_if_idom
2420     ! --- 1st domain - mixing layer
2421     
2422               lab_do_ib : do ib = 1, NBDSWLW
2423                 ext1 = f_zero
2424                 sca1 = f_zero
2425                 ssa1 = f_zero
2426                 asy1 = f_zero
2427     
2428                 lab_do_icmp : do icmp = 1, NXC
2429                   ic = idxc(icmp)
2430                   cm = cmix(icmp)
2431     
2432                   lab_if_ic : if (ic > NCM1) then
2433                     ic1 = ic - NCM1
2434     
2435                     ex00 = extrhd(ih1,ic1,ib)                               &
2436          &               + rdrh * (extrhd(ih2,ic1,ib) - extrhd(ih1,ic1,ib))
2437                     sc00 = scarhd(ih1,ic1,ib)                               &
2438          &               + rdrh * (scarhd(ih2,ic1,ib) - scarhd(ih1,ic1,ib))
2439                     ss00 = ssarhd(ih1,ic1,ib)                               &
2440          &               + rdrh * (ssarhd(ih2,ic1,ib) - ssarhd(ih1,ic1,ib))
2441                     as00 = asyrhd(ih1,ic1,ib)                               &
2442          &               + rdrh * (asyrhd(ih2,ic1,ib) - asyrhd(ih1,ic1,ib))
2443     
2444                     ext1 = ext1 + cm * ex00
2445                     sca1 = sca1 + cm * sc00
2446                     ssa1 = ssa1 + cm * ss00 * ex00
2447                     asy1 = asy1 + cm * as00 * sc00
2448                   else if (ic > 0) then     lab_if_ic
2449                     ext1 = ext1 + cm * extrhi(ic,ib)
2450                     sca1 = sca1 + cm * scarhi(ic,ib)
2451                     ssa1 = ssa1 + cm * ssarhi(ic,ib) * extrhi(ic,ib)
2452                     asy1 = asy1 + cm * asyrhi(ic,ib) * scarhi(ic,ib)
2453                   endif  lab_if_ic
2454     
2455                 enddo  lab_do_icmp
2456     
2457                 tauae(kk,ib) = ext1 * denn(1) * delz(kk)
2458                 ssaae(kk,ib) = min(f_one, ssa1/ext1)
2459                 asyae(kk,ib) = min(f_one, asy1/sca1)
2460               enddo  lab_do_ib
2461     
2462             elseif (idom == 2) then    lab_if_idom
2463     ! --- 2nd domain - mineral transport layers
2464     
2465               do ib = 1, NBDSWLW
2466                 tauae(kk,ib) = extrhi(6,ib) * denn(2) * delz(kk)
2467                 ssaae(kk,ib) = ssarhi(6,ib)
2468                 asyae(kk,ib) = asyrhi(6,ib)
2469               enddo
2470     
2471             else  lab_if_idom
2472     ! --- domain index out off range, assume no aerosol
2473     
2474               do ib = 1, NBDSWLW
2475                 tauae(kk,ib) = f_zero
2476                 ssaae(kk,ib) = f_one
2477                 asyae(kk,ib) = f_zero
2478               enddo
2479     
2480     !         write(6,19) kk,idom
2481     ! 19      format(/'  ***  ERROR in sub AEROS: domain index out'         &
2482     !    &,            ' of range!  K, IDOM =',3i5,' ***')
2483     !         stop 19
2484     
2485             endif  lab_if_idom
2486     
2487           enddo  lab_do_layer
2488     !
2489     !===> ... smooth profile at domain boundaries
2490     !
2491           if ( iflip == 0 ) then      ! input from toa to sfc
2492     
2493             do ib = 1, NBDSWLW
2494             do kk = 2, NLAY
2495               if ( tauae(kk,ib) > f_zero ) then
2496                 ratio = tauae(kk-1,ib) / tauae(kk,ib)
2497               else
2498                 ratio = f_one
2499               endif
2500     
2501               tt0 = tauae(kk,ib) + tauae(kk-1,ib)
2502               tt1 = 0.2 * tt0
2503               tt2 = tt0 - tt1
2504     
2505               if ( ratio > crt1 ) then
2506                 tauae(kk,ib)   = tt1
2507                 tauae(kk-1,ib) = tt2
2508               endif
2509     
2510               if ( ratio < crt2 ) then
2511                 tauae(kk,ib)   = tt2
2512                 tauae(kk-1,ib) = tt1
2513               endif
2514             enddo   ! do_kk_loop
2515             enddo   ! do_ib_loop
2516     
2517           else                      ! input from sfc to toa
2518     
2519             do ib = 1, NBDSWLW
2520             do kk = NLAY-1, 1, -1
2521               if ( tauae(kk,ib) > f_zero ) then
2522                 ratio = tauae(kk+1,ib) / tauae(kk,ib)
2523               else
2524                 ratio = f_one
2525               endif
2526     
2527               tt0 = tauae(kk,ib) + tauae(kk+1,ib)
2528               tt1 = 0.2 * tt0
2529               tt2 = tt0 - tt1
2530     
2531               if ( ratio > crt1 ) then
2532                 tauae(kk,ib)   = tt1
2533                 tauae(kk+1,ib) = tt2
2534               endif
2535     
2536               if ( ratio < crt2 ) then
2537                 tauae(kk,ib)   = tt2
2538                 tauae(kk+1,ib) = tt1
2539               endif
2540             enddo   ! do_kk_loop
2541             enddo   ! do_ib_loop
2542     
2543           endif
2544     
2545     !
2546           return
2547     !................................
2548           end subroutine radclimaer
2549     !--------------------------------
2550     !
2551     !...................................
2552           end subroutine setclimaer
2553     !-----------------------------------
2554     
2555     
2556     ! =======================================================================
2557     ! GOCART code modification starts here (Sarah lu)  ---------------------!
2558     !!
2559     !! gocart_init : set_aerspc, rd_gocart_clim, rd_gocart_luts, optavg_grt
2560     !! setgocartaer: aeropt_grt, map_aermr
2561     
2562     !-----------------------------------
2563           subroutine gocart_init                                            &
2564     !...................................
2565     !  ---  inputs:
2566          &     ( NWVTOT,solfwv,soltot,NWVTIR,eirfwv,                        &
2567          &       NBDSW,NBDIR,NBDSWLW,imon,me,raddt,fdaer                    &
2568     !  ---  outputs: ( none )
2569          &     )
2570     
2571     !  ==================================================================  !
2572     !                                                                      !
2573     !  subprogram : gocart_init                                            !
2574     !                                                                      !
2575     !    this is the initialization program for gocart aerosols            !
2576     !                                                                      !
2577     !    - determine weight and index for aerosol composition/luts         !
2578     !    - read in monthly global distribution of gocart aerosols          !
2579     !    - read and map the tabulated aerosol optical spectral data        !
2580     !        onto corresponding sw/lw radiation spectral bands.            !
2581     !                                                                      !
2582     !  ====================  defination of variables  ===================  !
2583     !                                                                      !
2584     !  inputs:                                                             !
2585     !   NWVTOT           - total num of wave numbers used in sw spectrum   !
2586     !   solfwv(NWVTOT)   - solar flux for each individual wavenumber (w/m2)!
2587     !   soltot           - total solar flux for the spectrual range  (w/m2)!
2588     !   NWVTIR           - total num of wave numbers used in the ir region !
2589     !   eirfwv(NWVTIR)   - ir flux(273k) for each individual wavenum (w/m2)!
2590     !   NBDSW            - num of bands calculated for sw aeros opt prop   !
2591     !   NBDIR            - num of bands calculated for lw aeros opt prop   !
2592     !   NBDSWLW          - total num of bands calc for sw+lw aeros opt prop!
2593     !   imon             - month of the year                               !
2594     !   me               - print message control flag                      !
2595     !                                                                      !
2596     !  outputs: (to the module variables)                                  !
2597     !                                                                      !
2598     !  module variables:                                                   !
2599     !     NBDSW   - total number of sw spectral bands                      !
2600     !     wvnum1,wvnum2 (NSWSTR:NSWEND)                                    !
2601     !             - start/end wavenumbers for each of sw bands             !
2602     !     NBDLW   - total number of lw spectral bands                      !
2603     !     wvnlw1,wvnlw2 (NBDLW)                                            !
2604     !             - start/end wavenumbers for each of lw bands             !
2605     !     NBDSWLW - total number of sw+lw bands used in this version       !
2606     !     extrhi_grt  - extinction coef for rh-indep aeros  KCM1*NBDSWLW   !
2607     !     ssarhi_grt  - single-scat-alb for rh-indep aeros  KCM1*NBDSWLW   !
2608     !     asyrhi_grt  - asymmetry factor for rh-indep aeros KCM1*NBDSWLW   !
2609     !     extrhd_grt  - extinction coef for rh-dep aeros    KRHLEV*KCM2*NBDSWLW!
2610     !     ssarhd_grt  - single-scat-alb for rh-dep aeros    KRHLEV*KCM2*NBDSWLW!
2611     !     asyrhd_grt  - asymmetry factor for rh-dep aeros   KRHLEV*KCM2*NBDSWLW!
2612     !     ctaer       - merging coefficients for fcst/clim fields          !
2613     !     get_fcst    - option to get fcst aerosol fields                  !
2614     !     get_clim    - option to get clim aerosol fields                  !
2615     !     dm_indx  - index for aer spec to be included in aeropt calculations  !
2616     !     dmfcs_indx  - index for prognostic aerosol fields                !
2617     !     psclmg      - geos3/4-gocart pressure            IMXG*JMXG*KMXG  !
2618     !     dmclmg      - geos3-gocart aerosol dry mass   IMXG*JMXG*KMXG*NMXG!
2619     !                   or geos4-gocart aerosol mixing ratio               !
2620     !                                                                      !
2621     !  usage:    call gocart_init                                          !
2622     !                                                                      !
2623     !  subprograms called:  set_aerspc, rd_gocart_clim,                    !
2624     !                       rd_gocart_luts, optavg_grt                     !
2625     !                                                                      !
2626     !  ==================================================================  !
2627     
2628           implicit none
2629     
2630     !  ---  inputs:
2631           integer, intent(in) :: NWVTOT,NWVTIR,NBDSW,NBDIR,NBDSWLW,imon,me
2632     
2633           real (kind=kind_phys), intent(in) :: raddt, fdaer
2634     
2635           real (kind=kind_phys), intent(in) :: solfwv(:),soltot, eirfwv(:)
2636     
2637     !  ---  output: ( none )
2638     
2639     !  ---  locals:
2640     
2641           real (kind=kind_phys), dimension(NBDSW,KAERBND) :: solwaer
2642           real (kind=kind_phys), dimension(NBDSW)         :: solbnd
2643           real (kind=kind_phys), dimension(NBDIR,KAERBND) :: eirwaer
2644           real (kind=kind_phys), dimension(NBDIR)         :: eirbnd
2645           real (kind=kind_phys) :: sumsol, sumir
2646     
2647           integer, dimension(NBDSW) :: nv1, nv2
2648           integer, dimension(NBDIR) :: nr1, nr2
2649     
2650           integer :: i, mb, ib, ii, iw, iw1, iw2
2651     
2652     !===>  ...  begin here
2653     
2654     !--------------------------------------------------------------------------
2655     !  (1) determine aerosol specification index and merging coefficients
2656     !--------------------------------------------------------------------------
2657     
2658           if ( .not. lgrtint ) then
2659     
2660     !  --- ...  already done aerspc initialization, continue
2661     
2662             continue
2663     
2664           else
2665     
2666     !  --- ...  set aerosol specification index and merging coefficients
2667     
2668             call set_aerspc(raddt,fdaer)
2669     !  ---  inputs:  (in scope variables)
2670     !  ---  outputs: (in scope variables)
2671     
2672           endif  ! end if_lgrtinit_block
2673     
2674     !
2675     !--------------------------------------------------------------------------
2676     !  (2) read gocart climatological data
2677     !--------------------------------------------------------------------------
2678     
2679     !  --- ...  read gocart climatological data, if needed
2680     
2681           if ( get_clim ) then
2682     
2683             call rd_gocart_clim
2684     !  ---  inputs:  (in scope variables)
2685     !  ---  outputs: (in scope variables)
2686     
2687           endif
2688     
2689     !
2690     !--------------------------------------------------------------------------
2691     !  (3) read and map the tabulated aerosol optical spectral data
2692     !           onto corresponding radiation spectral bands
2693     !--------------------------------------------------------------------------
2694     
2695           if ( .not. lgrtint ) then
2696     
2697     !  --- ...  already done optical property interpolation, exit
2698     
2699             return
2700     
2701           else
2702     
2703     !  --- ...  reset lgrtint
2704     
2705             lgrtint = .false.
2706     
2707     !  --- ...  read tabulated aerosol optical input data
2708             call rd_gocart_luts
2709     !  ---  inputs:  (in scope variables)
2710     !  ---  outputs: (in scope variables)
2711     
2712     !  --- ...  compute solar flux weights and interval indices for mapping
2713     !           spectral bands between sw radiation and aerosol data
2714     
2715             solbnd (:)   = f_zero
2716             solwaer(:,:) = f_zero
2717     
2718             nv_aod = 1
2719     
2720             do ib = 1, NBDSW
2721               mb = ib + NSWSTR - 1
2722               ii = 1
2723               iw1 = nint(wvnum1(mb))
2724               iw2 = nint(wvnum2(mb))
2725     !
2726     ! ---  locate the spectral band for 550nm (for aod diag)
2727     !
2728               if (10000./iw1 >= 0.55 .and.
2729          &        10000./iw2 <=0.55 )  then
2730                   nv_aod =  ib
2731               endif
2732     
2733               Lab_swdowhile : do while ( iw1 > iendwv_grt(ii) )
2734                 if ( ii == KAERBND ) exit Lab_swdowhile
2735                 ii = ii + 1
2736               enddo  Lab_swdowhile
2737     
2738               sumsol = f_zero
2739               nv1(ib) = ii
2740     
2741               do iw = iw1, iw2
2742                 solbnd(ib) = solbnd(ib) + solfwv(iw)
2743                 sumsol = sumsol + solfwv(iw)
2744     
2745                 if ( iw == iendwv_grt(ii) ) then
2746                   solwaer(ib,ii) = sumsol
2747     
2748                   if ( ii < KAERBND ) then
2749                     sumsol = f_zero
2750                     ii = ii + 1
2751                   endif
2752                 endif
2753               enddo
2754     
2755               if ( iw2 /= iendwv_grt(ii) ) then
2756                 solwaer(ib,ii) = sumsol
2757               endif
2758     
2759               nv2(ib) = ii
2760     
2761               if((me==0) .and. lckprnt) print *,'RAD-nv1,nv2:',
2762          &        ib,iw1,iw2,nv1(ib),iendwv_grt(nv1(ib)),
2763          &        nv2(ib),iendwv_grt(nv2(ib)),
2764          &        10000./iw1, 10000./iw2
2765             enddo     ! end do_ib_block for sw
2766     
2767     ! --- check the spectral range for the nv_550 band
2768             if((me==0) .and. lckprnt) then
2769               mb = nv_aod + NSWSTR - 1
2770               iw1 = nint(wvnum1(mb))
2771               iw2 = nint(wvnum2(mb))
2772                print *,'RAD-nv_aod:',
2773          &       nv_aod, iw1, iw2, 10000./iw1, 10000./iw2
2774             endif
2775     !
2776     !  --- ...  compute ir flux weights and interval indices for mapping
2777     !           spectral bands between lw radiation and aerosol data
2778     
2779             eirbnd (:)   = f_zero
2780             eirwaer(:,:) = f_zero
2781     
2782             do ib = 1, NBDIR
2783               ii = 1
2784               if ( NBDIR == 1 ) then
2785                 iw1 = 400                   ! corresponding 25 mu
2786                 iw2 = 2500                  ! corresponding 4  mu
2787               else
2788                 iw1 = nint(wvnlw1(ib))
2789                 iw2 = nint(wvnlw2(ib))
2790               endif
2791     
2792               Lab_lwdowhile : do while ( iw1 > iendwv_grt(ii) )
2793                 if ( ii == KAERBND ) exit Lab_lwdowhile
2794                 ii = ii + 1
2795               enddo  Lab_lwdowhile
2796     
2797               sumir = f_zero
2798               nr1(ib) = ii
2799     
2800               do iw = iw1, iw2
2801                 eirbnd(ib) = eirbnd(ib) + eirfwv(iw)
2802                 sumir  = sumir  + eirfwv(iw)
2803     
2804                 if ( iw == iendwv_grt(ii) ) then
2805                   eirwaer(ib,ii) = sumir
2806     
2807                   if ( ii < KAERBND ) then
2808                     sumir = f_zero
2809                     ii = ii + 1
2810                   endif
2811                 endif
2812               enddo
2813     
2814               if ( iw2 /= iendwv_grt(ii) ) then
2815                 eirwaer(ib,ii) = sumir
2816               endif
2817     
2818               nr2(ib) = ii
2819     
2820               if(me==0 .and. lckprnt) print *,'RAD-nr1,nr2:',
2821          &        ib,iw1,iw2,nr1(ib),iendwv_grt(nr1(ib)),
2822          &        nr2(ib),iendwv_grt(nr2(ib)),
2823          &        10000./iw1, 10000./iw2
2824             enddo     ! end do_ib_block for lw
2825     
2826     !  ---  compute spectral band mean properties for each species
2827     
2828             call optavg_grt
2829     !  ---  inputs:  (in scope variables)
2830     !  ---  outputs: (in scope variables)
2831     
2832             if(me==0 .and. lckprnt) then
2833               print *, 'RAD -After optavg_grt, sw band info'
2834               do ib = 1, NBDSW
2835                mb = ib + NSWSTR - 1
2836                print *,'RAD -wvnum1,wvnum2: ',ib,wvnum1(mb),wvnum2(mb)
2837                print *,'RAD -lamda1,lamda2: ',ib,10000./wvnum1(mb),
2838          &                                   10000./wvnum2(mb)
2839                print *,'RAD -extrhi_grt:', extrhi_grt(:,ib)
2840     !          do i = 1, KRHLEV
2841                do i = 1, KRHLEV, 10
2842                  print *, 'RAD -extrhd_grt:',i,rhlev_grt(i),
2843          &                                extrhd_grt(i,:,ib)
2844                enddo
2845               enddo
2846               print *, 'RAD -After optavg_grt, lw band info'
2847               do ib = 1, NBDIR
2848                ii = NBDSW + ib
2849                print *,'RAD -wvnlw1,wvnlw2: ',ib,wvnlw1(ib),wvnlw2(ib)
2850                print *,'RAD -lamda1,lamda2: ',ib,10000./wvnlw1(ib),
2851          &                                   10000./wvnlw2(ib)
2852                print *,'RAD -extrhi_grt:', extrhi_grt(:,ii)
2853     !          do i = 1, KRHLEV
2854                do i = 1, KRHLEV, 10
2855                  print *, 'RAD -extrhd_grt:',i,rhlev_grt(i),
2856          &                                extrhd_grt(i,:,ii)
2857                enddo
2858               enddo
2859             endif
2860     
2861     !  --- ...  dealoocate input data arrays no longer needed
2862             deallocate ( iendwv_grt   )
2863             if ( allocated(rhidext0_grt) ) then
2864               deallocate ( rhidext0_grt )
2865               deallocate ( rhidssa0_grt )
2866               deallocate ( rhidasy0_grt )
2867             endif
2868             if ( allocated(rhdpext0_grt) ) then
2869               deallocate ( rhdpext0_grt )
2870               deallocate ( rhdpssa0_grt )
2871               deallocate ( rhdpasy0_grt )
2872             endif
2873     
2874           endif  ! end if_lgrtinit_block
2875     
2876     ! =================
2877           contains
2878     ! =================
2879     
2880     !-----------------------------
2881           subroutine set_aerspc(raddt,fdaer)
2882     !.............................
2883     !  ---  inputs:  (in scope variables)
2884     !  ---  outputs: (in scope variables)
2885     
2886     ! ==================================================================== !
2887     !                                                                      !
2888     ! subprogram: set_aerspc                                               !
2889     !                                                                      !
2890     ! determine merging coefficients ctaer;                                !
2891     ! set up aerosol specification: num_gridcomp, gridcomp, dm_indx,       !
2892     !                       dmfcs_indx, isoot, iwaso, isuso, issam, isscm  !
2893     !                                                                      !
2894     ! Aerosol optical properties (ext, ssa, asy) are determined from       !
2895     ! NMGX (<=12) aerosol species                                          !
2896     ! ==> DU: dust1 (4 sub-micron bins), dust2, dust3, dust4, dust5        !
2897     !     BC: soot_phobic, soot_philic                                     !
2898     !     OC: waso_phobic, waso_philic                                     !
2899     !     SU: suso (=so4)                                                  !
2900     !     SS: ssam (accumulation mode), sscm (coarse mode)                 !
2901     !                                                                      !
2902     ! The current version only supports prognostic aerosols (from GOCART   !
2903     ! in-line calculations) and climo aerosols (from GEOS-GOCART runs)    !
2904     !                                                                      !
2905     !  ==================================================================  !
2906     !
2907           implicit none
2908     
2909     !  ---  inputs:
2910           real (kind=kind_phys), intent(in) :: raddt, fdaer
2911     !  ---  output:
2912     
2913     !  ---  local:
2914     !     real (kind=kind_phys)     :: raddt
2915           integer                   :: i, indxr
2916           character*2               :: tp, gridcomp_tmp(max_num_gridcomp)
2917     
2918     !! ===> determine ctaer (user specified weight for fcst fields)
2919     !     raddt = min(fhswr,fhlwr) / 24.
2920           if( fdaer >= 99999. ) ctaer = f_one
2921           if((fdaer>0.).and.(fdaer<99999.)) ctaer=exp(-raddt/fdaer)
2922     
2923           if(me==0 .and. lckprnt) then
2924             print *, 'RAD -raddt, fdaer,ctaer: ', raddt, fdaer, ctaer
2925             if (ctaer == f_one ) then
2926               print *, 'LU -aerosol fields determined from fcst'
2927             elseif (ctaer == f_zero) then
2928               print *, 'LU -aerosol fields determined from clim'
2929             else
2930               print *, 'LU -aerosol fields determined from fcst/clim'
2931             endif
2932           endif
2933     
2934     !! ===> determine get_fcst and get_clim
2935     !!    if fcst is chosen (ctaer == f_one ), set get_clim to F
2936     !!    if clim is chosen (ctaer == f_zero), set get_fcst to F
2937           if ( ctaer == f_one  )  get_clim = .false.
2938           if ( ctaer == f_zero )  get_fcst = .false.
2939     
2940     !! ===> determine aerosol species to be included in the calculations
2941     !!      of aerosol optical properties (ext, ssa, asy)
2942     
2943     !*  If climo option is chosen, the aerosol composition is hardwired
2944     !*  to full package. If not, the composition is determined from
2945     !*  tracer_config on-the-fly (full package or subset)
2946           lab_if_fcst : if ( get_fcst ) then
2947     
2948     !!      use tracer_config to determine num_gridcomp and gridcomp
2949             if ( gfs_phy_tracer%doing_GOCART )  then
2950              if ( gfs_phy_tracer%doing_DU )  then
2951                 num_gridcomp  =  num_gridcomp  + 1
2952                 gridcomp_tmp(num_gridcomp) = 'DU'
2953              endif
2954              if ( gfs_phy_tracer%doing_SU ) then
2955                 num_gridcomp  =  num_gridcomp  + 1
2956                 gridcomp_tmp(num_gridcomp) = 'SU'
2957              endif
2958              if ( gfs_phy_tracer%doing_SS ) then
2959                 num_gridcomp  =  num_gridcomp  + 1
2960                 gridcomp_tmp(num_gridcomp) = 'SS'
2961              endif
2962              if ( gfs_phy_tracer%doing_OC ) then
2963                 num_gridcomp  =  num_gridcomp  + 1
2964                 gridcomp_tmp(num_gridcomp) = 'OC'
2965              endif
2966              if ( gfs_phy_tracer%doing_BC ) then
2967                 num_gridcomp  =  num_gridcomp  + 1
2968                 gridcomp_tmp(num_gridcomp) = 'BC'
2969              endif
2970     !
2971              if ( num_gridcomp > 0 ) then
2972                allocate ( gridcomp(num_gridcomp) )
2973                gridcomp(1:num_gridcomp) = gridcomp_tmp(1:num_gridcomp)
2974              else
2975                print *,'ERROR: prognostic aerosols not found,abort',me
2976                stop 1000
2977              endif
2978     
2979             else      ! gfs_phy_tracer%doing_GOCART=F
2980     
2981              print *,'ERROR: prognostic aerosols option off, abort',me
2982             stop 1001
2983     
2984             endif     ! end_if_gfs_phy_tracer%doing_GOCART_if_
2985     
2986           else lab_if_fcst
2987     
2988     !!      set to full package (max_num_gridcomp and max_gridcomp)
2989             num_gridcomp = max_num_gridcomp
2990             allocate ( gridcomp(num_gridcomp) )
2991             gridcomp(1:num_gridcomp) = max_gridcomp(1:num_gridcomp)
2992     
2993           endif lab_if_fcst
2994     
2995     !!
2996     !! Aerosol specification is determined as such:
2997     !! A. For radiation-aerosol feedback, the specification is based on the aeropt
2998     !!    routine from Mian Chin and Hongbin Yu (hydrophobic and hydrophilic for
2999     !!    OC/BC; submicron and supermicron for SS, 8-bins (with 4 subgroups for the
3000     !!    the submicron bin) for DU, and SO4 for SU)
3001     !! B. For transport, the specification is determined from GOCART in-line module
3002     !! C. For LUTS, (waso, soot, ssam, sscm, suso, dust) is used, based on the
3003     !!    the OPAC climo aerosol scheme (implemented by Yu-Tai Hou)
3004     
3005     !!=== <A> determine dm_indx and NMXG
3006           indxr = 0
3007           dm_indx%waso_phobic = -999         ! OC
3008           dm_indx%soot_phobic = -999         ! BC
3009           dm_indx%ssam = -999                ! SS
3010           dm_indx%suso = -999                ! SU
3011           dm_indx%dust1 = -999               ! DU
3012           do i = 1, num_gridcomp
3013              tp = gridcomp(i)
3014              select case ( tp )
3015              case ( 'OC')    ! consider hydrophobic and hydrophilic
3016                dm_indx%waso_phobic = indxr + 1
3017                dm_indx%waso_philic = indxr + 2
3018                indxr = indxr + 2
3019              case ( 'BC')    ! consider hydrophobic and hydrophilic
3020                dm_indx%soot_phobic = indxr + 1
3021                dm_indx%soot_philic = indxr + 2
3022                indxr = indxr + 2
3023              case ( 'SS')    ! consider submicron and supermicron
3024                dm_indx%ssam = indxr + 1
3025                dm_indx%sscm = indxr + 2
3026                indxr = indxr + 2
3027              case ( 'SU')    ! consider SO4 only
3028                dm_indx%suso = indxr + 1
3029                indxr = indxr + 1
3030              case ( 'DU')    ! consider all 5 bins
3031                dm_indx%dust1 = indxr + 1
3032                dm_indx%dust2 = indxr + 2
3033                dm_indx%dust3 = indxr + 3
3034                dm_indx%dust4 = indxr + 4
3035                dm_indx%dust5 = indxr + 5
3036                indxr = indxr + 5
3037              case default
3038                print *,'ERROR: aerosol species not supported, abort',me
3039                stop 1002
3040              end select
3041           enddo
3042     !!
3043           NMXG       = indxr      ! num of gocart aer spec for opt cal
3044     !!
3045     
3046     !!=== <B> determine dmfcs_indx
3047     !!    SS: 5-bins are considered for transport while only two groups
3048     !!        (accumulation/coarse modes) are considered for radiation
3049     !!    DU: 5-bins are considered for transport while 8 bins (with the
3050     !!        submicorn bin exptended to 4 bins) are considered for radiation
3051     !!    SU: DMS, SO2, and MSA are not considered for radiation
3052     
3053           if ( get_fcst ) then
3054              if ( gfs_phy_tracer%doing_OC )  then
3055                 dmfcs_indx%ocphobic = trcindx ('ocphobic', gfs_phy_tracer)
3056                 dmfcs_indx%ocphilic = trcindx ('ocphilic', gfs_phy_tracer)
3057              endif
3058              if ( gfs_phy_tracer%doing_BC )  then
3059                 dmfcs_indx%bcphobic = trcindx ('bcphobic', gfs_phy_tracer)
3060                 dmfcs_indx%bcphilic = trcindx ('bcphilic', gfs_phy_tracer)
3061              endif
3062              if ( gfs_phy_tracer%doing_SS )  then
3063                 dmfcs_indx%ss001 = trcindx ('ss001', gfs_phy_tracer)
3064                 dmfcs_indx%ss002 = trcindx ('ss002', gfs_phy_tracer)
3065                 dmfcs_indx%ss003 = trcindx ('ss003', gfs_phy_tracer)
3066                 dmfcs_indx%ss004 = trcindx ('ss004', gfs_phy_tracer)
3067                 dmfcs_indx%ss005 = trcindx ('ss005', gfs_phy_tracer)
3068              endif
3069              if ( gfs_phy_tracer%doing_SU )  then
3070                 dmfcs_indx%so4 = trcindx ('so4', gfs_phy_tracer)
3071              endif
3072              if ( gfs_phy_tracer%doing_DU )  then
3073                 dmfcs_indx%du001 = trcindx ('du001', gfs_phy_tracer)
3074                 dmfcs_indx%du002 = trcindx ('du002', gfs_phy_tracer)
3075                 dmfcs_indx%du003 = trcindx ('du003', gfs_phy_tracer)
3076                 dmfcs_indx%du004 = trcindx ('du004', gfs_phy_tracer)
3077                 dmfcs_indx%du005 = trcindx ('du005', gfs_phy_tracer)
3078              endif
3079           endif
3080     
3081     !!
3082     !!=== <C> determin KCM, KCM1, KCM2
3083     !!    DU: submicron bin (dust1) contains 4 sub-groups (e.g., hardwire
3084     !!        8 bins for aerosol optical properties luts)
3085     !!    OC/BC: while hydrophobic aerosols are rh-independent, the luts
3086     !!        for hydrophilic aerosols are used (e.g., use the coeff
3087     !!        corresponding to rh=0)
3088     !!
3089           indxr = 1
3090           isoot = -999
3091           iwaso = -999
3092           isuso = -999
3093           issam = -999
3094           isscm = -999
3095           do i = 1, num_gridcomp
3096              tp = gridcomp(i)
3097              if ( tp /= 'DU' ) then  !<--- non-dust aerosols
3098                select case ( tp )
3099                case ( 'OC ')
3100                  iwaso = indxr
3101                case ( 'BC ')
3102                  isoot = indxr
3103                case ( 'SU ')
3104                  isuso = indxr
3105                case ( 'SS ')
3106                  issam = indxr
3107                  isscm = indxr + 1
3108                end select
3109                if ( tp /= 'SS' ) then
3110                  indxr = indxr + 1
3111                else
3112                  indxr = indxr + 2
3113                endif
3114              else                   !<--- dust aerosols
3115                KCM1 =  8            ! num of rh independent aer species
3116              endif
3117           enddo
3118           KCM2 = indxr - 1          ! num of rh dependent aer species
3119           KCM  = KCM1 + KCM2        ! total num of aer species
3120     
3121     !!
3122     !! check print starts here
3123           if( me == 0 .and. lckprnt) then
3124            print *, 'RAD -num_gridcomp:', num_gridcomp
3125            print *, 'RAD -gridcomp    :', gridcomp(:)
3126            print *, 'RAD -NMXG:',  NMXG
3127            print *, 'RAD -dm_indx ===> '
3128            print *, 'RAD -aerspc: dust1=', dm_indx%dust1
3129            print *, 'RAD -aerspc: dust2=', dm_indx%dust2
3130            print *, 'RAD -aerspc: dust3=', dm_indx%dust3
3131            print *, 'RAD -aerspc: dust4=', dm_indx%dust4
3132            print *, 'RAD -aerspc: dust5=', dm_indx%dust5
3133            print *, 'RAD -aerspc: ssam=',  dm_indx%ssam
3134            print *, 'RAD -aerspc: sscm=',  dm_indx%sscm
3135            print *, 'RAD -aerspc: suso=',  dm_indx%suso
3136            print *, 'RAD -aerspc: waso_phobic=',dm_indx%waso_phobic
3137            print *, 'RAD -aerspc: waso_philic=',dm_indx%waso_philic
3138            print *, 'RAD -aerspc: soot_phobic=',dm_indx%soot_phobic
3139            print *, 'RAD -aerspc: soot_philic=',dm_indx%soot_philic
3140     
3141            print *, 'RAD -KCM1 =', KCM1
3142            print *, 'RAD -KCM2 =', KCM2
3143            print *, 'RAD -KCM  =', KCM
3144            if ( KCM2 > 0 ) then
3145              print *, 'RAD -aerspc: issam=', issam
3146              print *, 'RAD -aerspc: isscm=', isscm
3147              print *, 'RAD -aerspc: isuso=', isuso
3148              print *, 'RAD -aerspc: iwaso=', iwaso
3149              print *, 'RAD -aerspc: isoot=', isoot
3150            endif
3151     
3152            if ( get_fcst ) then
3153              print *, 'RAD -dmfcs_indx ===> '
3154              print *, 'RAD -trc_du001=',dmfcs_indx%du001
3155              print *, 'RAD -trc_du002=',dmfcs_indx%du002
3156              print *, 'RAD -trc_du003=',dmfcs_indx%du003
3157              print *, 'RAD -trc_du004=',dmfcs_indx%du004
3158              print *, 'RAD -trc_du005=',dmfcs_indx%du005
3159              print *, 'RAD -trc_so4  =',dmfcs_indx%so4
3160              print *, 'RAD -trc_ocphobic=',dmfcs_indx%ocphobic
3161              print *, 'RAD -trc_ocphilic=',dmfcs_indx%ocphilic
3162              print *, 'RAD -trc_bcphobic=',dmfcs_indx%bcphobic
3163              print *, 'RAD -trc_bcphilic=',dmfcs_indx%bcphilic
3164              print *, 'RAD -trc_ss001=',dmfcs_indx%ss001
3165              print *, 'RAD -trc_ss002=',dmfcs_indx%ss002
3166              print *, 'RAD -trc_ss003=',dmfcs_indx%ss003
3167              print *, 'RAD -trc_ss004=',dmfcs_indx%ss004
3168              print *, 'RAD -trc_ss005=',dmfcs_indx%ss005
3169            endif
3170           endif
3171     !! check print ends here
3172     
3173           return
3174     !                                                                      !
3175           end subroutine set_aerspc
3176     
3177     !-----------------------------------
3178     !-----------------------------
3179           subroutine rd_gocart_luts
3180     !.............................
3181     !  ---  inputs:  (in scope variables)
3182     !  ---  outputs: (in scope variables)
3183     
3184     ! ==================================================================== !
3185     ! subprogram: rd_gocart_luts                                           !
3186     !   read input gocart aerosol optical data from Mie code calculations  !
3187     !                                                                      !
3188     ! Remarks (Quanhua (Mark) Liu, JCSDA, June 2008)                       !
3189     !  The LUT is for NCEP selected 61 wave numbers and 6 aerosols         !
3190     !  (dust, soot, suso, waso, ssam, and sscm) and 36 aerosol effective   !
3191     !  size in microns.                                                    !
3192     !                                                                      !
3193     !  The LUT is computed using Mie code with a logorithm size            !
3194     !  distribution for each of 36 effective sizes. The standard deviation !
3195     !  sigma of the size, and min/max size follows Chin et al. 2000        !
3196     !  For each effective size, it corresponds a relative humidity value.  !
3197     !                                                                      !
3198     !  The LUT contains the density, sigma, relative humidity, mean mode   !
3199     !  radius, effective size, mass extinction coefficient, single         !
3200     !  scattering albedo, asymmetry factor, and phase function             !
3201     !                                                                      !
3202     !  ==================================================================  !
3203     !
3204           implicit none
3205     
3206     !  ---  inputs:
3207     !  ---  output:
3208     
3209     !  ---  locals:
3210           INTEGER, PARAMETER :: NP = 100, NP2 = 2*NP, nWave=100,
3211          &                      nAero=6, n_p=36
3212           INTEGER :: NW, NS, nH, n_bin
3213           real (kind=kind_io8), Dimension( NP2 ) :: Angle, Cos_Angle,
3214          &                                          Cos_Weight
3215           real (kind=kind_io8), Dimension(n_p,nAero) :: RH, rm, reff
3216           real (kind=kind_io8), Dimension(nWave,n_p,nAero) ::
3217          &                      ext0, sca0, asy0
3218           real (kind=kind_io8), Dimension(NP2,n_p,nWave,nAero) :: ph0
3219           real (kind=kind_io8) :: wavelength(nWave), density(nAero),
3220          &                        sigma(nAero), wave,n_fac,PI,t1,s1,g1
3221           CHARACTER(len=80) :: AerosolName(nAero)
3222           INTEGER    :: i, j, k, l, ij
3223     
3224           character  :: aerosol_file*30
3225           logical    :: file_exist
3226           integer    :: indx_dust(8)          ! map 36 dust bins to gocart size bins
3227     
3228           data aerosol_file  /"NCEP_AEROSOL.bin"/
3229           data AerosolName/ ' Dust ', ' Soot ', ' SUSO ', ' WASO ',
3230          &                  ' SSAM ', ' SSCM '/
3231     
3232     !! 8 dust bins
3233     !!  1       2       3       4       5       6       7       8
3234     !! .1-.18, .18-.3, .3-.6, 0.6-1.0, 1.0-1.8, 1.8-3, 3-6,  6-10  <-- def
3235     !!  0.1399  0.2399  0.4499 0.8000 1.3994  2.3964 4.4964  7.9887 <-- reff
3236           data indx_dust/4, 8, 12, 18, 21, 24, 30, 36/
3237     
3238           PI = acos(-1.d0)
3239     
3240     ! -- allocate aerosol optical data
3241           if ( .not. allocated( iendwv_grt ) ) then
3242             allocate ( iendwv_grt (KAERBND) )
3243           endif
3244           if (.not. allocated(rhidext0_grt) .and. KCM1 > 0 ) then
3245             allocate ( rhidext0_grt(KAERBND,KCM1))
3246             allocate ( rhidssa0_grt(KAERBND,KCM1))
3247             allocate ( rhidasy0_grt(KAERBND,KCM1))
3248           endif
3249           if (.not. allocated(rhdpext0_grt) .and. KCM2 > 0 ) then
3250             allocate ( rhdpext0_grt(KAERBND,KRHLEV,KCM2))
3251             allocate ( rhdpssa0_grt(KAERBND,KRHLEV,KCM2))
3252             allocate ( rhdpasy0_grt(KAERBND,KRHLEV,KCM2))
3253           endif
3254     
3255     ! -- read luts
3256           inquire (file = aerosol_file, exist = file_exist)
3257     
3258           if ( file_exist ) then
3259             if(me==0 .and. lckprnt) print *,'RAD -open :',aerosol_file
3260             close (NIAERCM)
3261             open (unit=NIAERCM,file=aerosol_file,status='OLD',              &
3262          &        form='UNFORMATTED')
3263           else
3264             print *,'    Requested aerosol data file "',aerosol_file,       &
3265          &          '" not found!', me
3266             print *,'    *** Stopped in subroutine RD_GOCART_LUTS !!'
3267             stop 1003
3268           endif              ! end if_file_exist_block
3269     
3270           READ(NIAERCM) (Cos_Angle(i),i=1,NP)
3271           READ(NIAERCM) (Cos_Weight(i),i=1,NP)
3272           READ(NIAERCM)
3273           READ(NIAERCM)
3274           READ(NIAERCM) NW,NS
3275           READ(NIAERCM)
3276           READ(NIAERCM) (wavelength(i),i=1,NW)
3277     
3278     ! --- check nAero and NW
3279           if (NW /= KAERBND) then
3280             print *, "Incorrect spectral band, abort ", NW
3281             stop 1004
3282           endif
3283     
3284     ! --- convert wavelength to wavenumber
3285           do i = 1, KAERBND
3286            iendwv_grt(i) = 10000. / wavelength(i)
3287            if(me==0 .and. lckprnt) print *,'RAD -wn,lamda:',
3288          &           i,iendwv_grt(i),wavelength(i)
3289           enddo
3290     
3291           DO j = 1, nAero
3292             if(me==0 .and. lckprnt) print *,'RAD -read LUTs:',
3293          &            j,AerosolName(j)
3294             READ(NIAERCM)
3295             READ(NIAERCM)
3296             READ(NIAERCM) n_bin, density(j), sigma(j)
3297             READ(NIAERCM)
3298             READ(NIAERCM) (RH(i,j),i=1, n_bin)
3299             READ(NIAERCM)
3300             READ(NIAERCM) (rm(i,j),i=1, n_bin)
3301             READ(NIAERCM)
3302             READ(NIAERCM) (reff(i,j),i=1, n_bin)
3303     
3304     ! --- check n_bin
3305             if (n_bin /= KRHLEV ) then
3306               print *, "Incorrect rh levels, abort ", n_bin
3307               stop 1005
3308             endif
3309     
3310     ! --- read luts
3311             DO k = 1, NW
3312               READ(NIAERCM) wave,(ext0(k,L,j),L=1,n_bin)
3313               READ(NIAERCM) (sca0(k,L,j),L=1,n_bin)
3314               READ(NIAERCM) (asy0(k,L,j),L=1,n_bin)
3315               READ(NIAERCM) (ph0(1:NP2,L,k,j),L=1,n_bin)
3316             END DO
3317     
3318     ! --- map luts input to module variables
3319             if (AerosolName(j) == ' Dust ' ) then
3320              if ( KCM1 > 0) then    !<-- only if rh independent aerosols are needed
3321               do i = 1, KCM1
3322                rhidext0_grt(1:KAERBND,i)=ext0(1:KAERBND,indx_dust(i),j)
3323                rhidssa0_grt(1:KAERBND,i)=sca0(1:KAERBND,indx_dust(i),j)
3324                rhidasy0_grt(1:KAERBND,i)=asy0(1:KAERBND,indx_dust(i),j)
3325               enddo
3326              endif
3327             else
3328              if ( KCM2 > 0) then    !<-- only if rh dependent aerosols are needed
3329               if (AerosolName(j) == ' Soot ') ij = isoot
3330               if (AerosolName(j) == ' SUSO ') ij = isuso
3331               if (AerosolName(j) == ' WASO ') ij = iwaso
3332               if (AerosolName(j) == ' SSAM ') ij = issam
3333               if (AerosolName(j) == ' SSCM ') ij = isscm
3334               if ( ij .ne. -999 ) then
3335                  rhdpext0_grt(1:KAERBND,1:KRHLEV,ij) =
3336          &               ext0(1:KAERBND,1:KRHLEV,j)
3337                  rhdpssa0_grt(1:KAERBND,1:KRHLEV,ij) =
3338          &               sca0(1:KAERBND,1:KRHLEV,j)
3339                  rhdpasy0_grt(1:KAERBND,1:KRHLEV,ij) =
3340          &               asy0(1:KAERBND,1:KRHLEV,j)
3341               endif   ! if_ij
3342              endif    ! if_KCM2
3343             endif
3344           END DO
3345     
3346           return
3347     !...................................
3348           end subroutine rd_gocart_luts
3349     !-----------------------------------
3350     !                                                                      !
3351     !-----------------------------
3352           subroutine optavg_grt
3353     !.............................
3354     !  ---  inputs:  (in scope variables)
3355     !  ---  outputs: (in scope variables)
3356     
3357     ! ==================================================================== !
3358     !                                                                      !
3359     ! subprogram: optavg_grt                                               !
3360     !                                                                      !
3361     !   compute mean aerosols optical properties over each sw/lw radiation !
3362     !   spectral band for each of the species components.  This program    !
3363     !   follows gfdl's approach for thick cloud opertical property in      !
3364     !   sw radiation scheme (2000).                                        !
3365     !                                                                      !
3366     !  ====================  defination of variables  ===================  !
3367     !                                                                      !
3368     ! input arguments:                                                     !
3369     !   nv1,nv2 (NBDSW)  - start/end spectral band indices of aerosol data !
3370     !                      for each sw radiation spectral band             !
3371     !   nr1,nr2 (NBDIR)  - start/end spectral band indices of aerosol data !
3372     !                      for each ir radiation spectral band             !
3373     !   solwaer (NBDSW,KAERBND)                                            !
3374     !                    - solar flux weight over each sw radiation band   !
3375     !                      vs each aerosol data spectral band              !
3376     !   eirwaer (NBDIR,KAERBND)                                            !
3377     !                    - ir flux weight over each lw radiation band      !
3378     !                      vs each aerosol data spectral band              !
3379     !   solbnd  (NBDSW)  - solar flux weight over each sw radiation band   !
3380     !   eirbnd  (NBDIR)  - ir flux weight over each lw radiation band      !
3381     !   NBDSW            - total number of sw spectral bands               !
3382     !   NBDIR            - total number of lw spectral bands               !
3383     !   NBDSWLW          - total number of sw+lw spectral bands            !
3384     !                                                                      !
3385     ! output arguments: (to module variables)                              !
3386     !                                                                      !
3387     !  ==================================================================  !
3388     !
3389           implicit none
3390     
3391     !  ---  inputs:
3392     !  ---  output:
3393     
3394     !  ---  locals:
3395           real (kind=kind_phys) :: sumk, sumok, sumokg, sumreft,            &
3396          &       sp, refb, reft, rsolbd, rirbd
3397     
3398           integer :: ib, nb, ni, nh, nc
3399     !
3400     !===> ...  begin here
3401     
3402     !  --- ...  allocate aerosol optical data
3403           if (.not. allocated(extrhd_grt) .and. KCM2 > 0 ) then
3404             allocate ( extrhd_grt(KRHLEV,KCM2,NBDSWLW) )
3405             allocate ( ssarhd_grt(KRHLEV,KCM2,NBDSWLW) )
3406             allocate ( asyrhd_grt(KRHLEV,KCM2,NBDSWLW) )
3407           endif
3408           if (.not. allocated(extrhi_grt) .and. KCM1 > 0 ) then
3409             allocate ( extrhi_grt(KCM1,NBDSWLW) )
3410             allocate ( ssarhi_grt(KCM1,NBDSWLW) )
3411             allocate ( asyrhi_grt(KCM1,NBDSWLW) )
3412           endif
3413     !
3414     !  --- ...  loop for each sw radiation spectral band
3415     
3416           do nb = 1, NBDSW
3417             rsolbd = f_one / solbnd(nb)
3418     
3419     !  ---  for rh independent aerosol species
3420     
3421             lab_rhi: if (KCM1 >  0 ) then
3422             do nc = 1, KCM1
3423               sumk    = f_zero
3424               sumok   = f_zero
3425               sumokg  = f_zero
3426               sumreft = f_zero
3427     
3428               do ni = nv1(nb), nv2(nb)
3429                 sp   = sqrt( (f_one - rhidssa0_grt(ni,nc))                  &
3430          &           / (f_one - rhidssa0_grt(ni,nc)*rhidasy0_grt(ni,nc)) )
3431                 reft = (f_one - sp) / (f_one + sp)
3432                 sumreft = sumreft + reft*solwaer(nb,ni)
3433     
3434                 sumk    = sumk    + rhidext0_grt(ni,nc)*solwaer(nb,ni)
3435                 sumok   = sumok   + rhidssa0_grt(ni,nc)*solwaer(nb,ni)      &
3436          &              * rhidext0_grt(ni,nc)
3437                 sumokg  = sumokg  + rhidssa0_grt(ni,nc)*solwaer(nb,ni)      &
3438          &              * rhidext0_grt(ni,nc)*rhidasy0_grt(ni,nc)
3439               enddo
3440     
3441               refb = sumreft * rsolbd
3442     
3443               extrhi_grt(nc,nb) = sumk   * rsolbd
3444               asyrhi_grt(nc,nb) = sumokg / (sumok + 1.0e-10)
3445               ssarhi_grt(nc,nb) = 4.0*refb                                  &
3446          &      / ( (f_one+refb)**2 - asyrhi_grt(nc,nb)*(f_one-refb)**2 )
3447     
3448             enddo   ! end do_nc_block for rh-ind aeros
3449             endif lab_rhi
3450     
3451     !  ---  for rh dependent aerosols species
3452     
3453             lab_rhd: if (KCM2 > 0 ) then
3454             do nc = 1, KCM2
3455               do nh = 1, KRHLEV
3456                 sumk    = f_zero
3457                 sumok   = f_zero
3458                 sumokg  = f_zero
3459                 sumreft = f_zero
3460     
3461                 do ni = nv1(nb), nv2(nb)
3462                   sp   = sqrt( (f_one - rhdpssa0_grt(ni,nh,nc))             &
3463          &        /(f_one-rhdpssa0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)))
3464                   reft = (f_one - sp) / (f_one + sp)
3465                   sumreft = sumreft + reft*solwaer(nb,ni)
3466     
3467                   sumk    = sumk   + rhdpext0_grt(ni,nh,nc)*solwaer(nb,ni)
3468                   sumok   = sumok  + rhdpssa0_grt(ni,nh,nc)*solwaer(nb,ni)  &
3469          &                * rhdpext0_grt(ni,nh,nc)
3470                   sumokg  = sumokg + rhdpssa0_grt(ni,nh,nc)*solwaer(nb,ni)  &
3471          &                * rhdpext0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)
3472                 enddo
3473     
3474                 refb = sumreft * rsolbd
3475     
3476                 extrhd_grt(nh,nc,nb) = sumk   * rsolbd
3477                 asyrhd_grt(nh,nc,nb) = sumokg / (sumok + 1.0e-10)
3478                 ssarhd_grt(nh,nc,nb) = 4.0*refb                             &
3479          &      /((f_one+refb)**2 - asyrhd_grt(nh,nc,nb)*(f_one-refb)**2)
3480               enddo   ! end do_nh_block
3481             enddo   ! end do_nc_block for rh-dep aeros
3482             endif lab_rhd
3483     
3484           enddo   !  end do_nb_block for sw
3485     
3486     !  --- ...  loop for each lw radiation spectral band
3487     
3488           do nb = 1, NBDIR
3489     
3490             ib = NBDSW + nb
3491             rirbd = f_one / eirbnd(nb)
3492     
3493     !  ---  for rh independent aerosol species
3494     
3495             lab_rhi_lw: if (KCM1 > 0 ) then
3496             do nc = 1, KCM1
3497               sumk    = f_zero
3498               sumok   = f_zero
3499               sumokg  = f_zero
3500               sumreft = f_zero
3501     
3502               do ni = nr1(nb), nr2(nb)
3503                 sp   = sqrt( (f_one - rhidssa0_grt(ni,nc))                  &
3504          &      / (f_one - rhidssa0_grt(ni,nc)*rhidasy0_grt(ni,nc)) )
3505                 reft = (f_one - sp) / (f_one + sp)
3506                 sumreft = sumreft + reft*eirwaer(nb,ni)
3507     
3508                 sumk    = sumk    + rhidext0_grt(ni,nc)*eirwaer(nb,ni)
3509                 sumok   = sumok   + rhidssa0_grt(ni,nc)*eirwaer(nb,ni)      &
3510          &              * rhidext0_grt(ni,nc)
3511                 sumokg  = sumokg  + rhidssa0_grt(ni,nc)*eirwaer(nb,ni)      &
3512          &              * rhidext0_grt(ni,nc)*rhidasy0_grt(ni,nc)
3513               enddo
3514     
3515               refb = sumreft * rirbd
3516     
3517               extrhi_grt(nc,ib) = sumk   * rirbd
3518               asyrhi_grt(nc,ib) = sumokg / (sumok + 1.0e-10)
3519               ssarhi_grt(nc,ib) = 4.0*refb                                  &
3520          &    / ( (f_one+refb)**2 - asyrhi_grt(nc,ib)*(f_one-refb)**2 )
3521             enddo   ! end do_nc_block for rh-ind aeros
3522             endif lab_rhi_lw
3523     
3524     !  ---  for rh dependent aerosols species
3525     
3526             lab_rhd_lw: if (KCM2 > 0 ) then
3527             do nc = 1, KCM2
3528               do nh = 1, KRHLEV
3529                 sumk    = f_zero
3530                 sumok   = f_zero
3531                 sumokg  = f_zero
3532                 sumreft = f_zero
3533     
3534                 do ni = nr1(nb), nr2(nb)
3535                   sp   = sqrt( (f_one - rhdpssa0_grt(ni,nh,nc))             &
3536          &        /(f_one - rhdpssa0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)) )
3537                   reft = (f_one - sp) / (f_one + sp)
3538                   sumreft = sumreft + reft*eirwaer(nb,ni)
3539     
3540                   sumk    = sumk  + rhdpext0_grt(ni,nh,nc)*eirwaer(nb,ni)
3541                   sumok   = sumok + rhdpssa0_grt(ni,nh,nc)*eirwaer(nb,ni)   &
3542          &                * rhdpext0_grt(ni,nh,nc)
3543                   sumokg  = sumokg+ rhdpssa0_grt(ni,nh,nc)*eirwaer(nb,ni)   &
3544          &                * rhdpext0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)
3545                 enddo
3546     
3547                 refb = sumreft * rirbd
3548     
3549                 extrhd_grt(nh,nc,ib) = sumk   * rirbd
3550                 asyrhd_grt(nh,nc,ib) = sumokg / (sumok + 1.0e-10)
3551                 ssarhd_grt(nh,nc,ib) = 4.0*refb                             &
3552          &      /((f_one+refb)**2 - asyrhd_grt(nh,nc,ib)*(f_one-refb)**2 )
3553               enddo   ! end do_nh_block
3554             enddo   ! end do_nc_block for rh-dep aeros
3555             endif lab_rhd_lw
3556     
3557           enddo   !  end do_nb_block for lw
3558     
3559     !
3560           return
3561     !................................
3562           end subroutine optavg_grt
3563     !--------------------------------
3564     !
3565     !-----------------------------------
3566           subroutine rd_gocart_clim                                         &
3567     !...................................
3568     !  ---  inputs:  (in scope variables)
3569     !  ---  outputs: (in scope variables)
3570     
3571     !  ==================================================================  !
3572     !                                                                      !
3573     ! subprogram: rd_gocart_clim                                           !
3574     !                                                                      !
3575     !   1. read in aerosol dry mass and surface pressure from GEOS3-GOCART !
3576     !      C3.1 2000 monthly data set                                      !
3577     !      or aerosol mixing ratio and surface pressure from GEOS4-GOCART  !
3578     !      2000-2007 averaged monthly data set                             !
3579     !   2. compute goes lat/lon array (for horizontal mapping)             !
3580     !                                                                      !
3581     !  ====================  defination of variables  ===================  !
3582     !                                                                      !
3583     ! inputs arguments:                                                    !
3584     !     imon    - month of the year                                      !
3585     !     me      - print message control flag                             !
3586     !                                                                      !
3587     ! outputs arguments: (to the module variables)                         !
3588     !     psclmg   - pressure (sfc to toa)    cb   IMXG*JMXG*KMXG          !
3589     !     dmclmg   - aerosol dry mass/mixing ratio     IMXG*JMXG*KMXG*NMXG !
3590     !     geos_rlon - goes longitude          deg      IMXG                !
3591     !     geos_rlat - goes latitude           deg      JMXG                !
3592     !                                                                      !
3593     !  usage:    call rd_gocart_clim                                       !
3594     !                                                                      !
3595     !  program history:                                                    !
3596     !    05/18/2010  ---  Lu    Add the option to read GEOS4-GOCART climo  !
3597     !  ==================================================================  !
3598     !
3599           implicit none
3600     
3601     !  ---  inputs:
3602     !  ---  output:
3603     
3604     !  ---  locals:
3605           integer, parameter :: MAXSPC = 5
3606           real (kind=kind_io4), parameter  :: PINT = 0.01
3607           real (kind=kind_io4), parameter  :: EPSQ = 0.0
3608     
3609           integer         :: i, j, k, numspci, ii
3610           integer         :: icmp, nrecl, nt1, nt2, nn(MAXSPC)
3611           character       :: ymd*6, yr*4, mn*2, tp*2,
3612          &                   fname*30, fin*30, aerosol_file*40
3613           logical         :: file_exist
3614     
3615           real (kind=kind_io4), dimension(KMXG)             :: sig
3616           real (kind=kind_io4), dimension(IMXG,JMXG)        :: ps
3617           real (kind=kind_io4), dimension(IMXG,JMXG,KMXG)   :: temp
3618           real (kind=kind_io4), dimension(IMXG,JMXG,KMXG,MAXSPC):: buff
3619           real (kind=kind_phys)   :: pstmp
3620     
3621     !     Add the following variables for GEOS4-GOCART
3622           real (kind=kind_io4), dimension(KMXG):: hyam, hybm
3623           real (kind=kind_io4)                 :: p0
3624     
3625           data yr /'2000'/           !!<=== use 2000 as the climo proxy
3626     
3627     !* sigma_coordinate for GEOS3-GOCART
3628     !* P(i,j,k) = PINT + SIG(k) * (PS(i,j) - PINT)
3629           data SIG  /
3630          &     9.98547E-01,9.94147E-01,9.86350E-01,9.74300E-01,9.56950E-01,
3631          &     9.33150E-01,9.01750E-01,8.61500E-01,8.11000E-01,7.50600E-01,
3632          &     6.82900E-01,6.10850E-01,5.37050E-01,4.63900E-01,3.93650E-01,
3633          &     3.28275E-01,2.69500E-01,2.18295E-01,1.74820E-01,1.38840E-01,
3634          &     1.09790E-01,8.66900E-02,6.84150E-02,5.39800E-02,4.25750E-02,
3635          &     3.35700E-02,2.39900E-02,1.36775E-02,5.01750E-03,5.30000E-04 /
3636     
3637     !* hybrid_sigma_pressure_coordinate for GEOS4-GOCART
3638     !* p(i,j,k) = a(k)*p0 + b(k)*ps(i,j)
3639           data hyam/
3640          &   0, 0.0062694, 0.02377049, 0.05011813, 0.08278809, 0.1186361,
3641          &   0.1540329, 0.1836373, 0.2043698, 0.2167788, 0.221193,
3642          &   0.217729, 0.2062951, 0.1865887, 0.1615213, 0.1372958,
3643          &   0.1167039, 0.09920014, 0.08432171, 0.06656809, 0.04765031,
3644          &   0.03382346, 0.0237648, 0.01435208, 0.00659734, 0.002826232,
3645          &   0.001118959, 0.0004086494, 0.0001368611, 3.750308e-05/
3646     
3647           data hybm /
3648          &   0.992555, 0.9642, 0.90556, 0.816375, 0.703815, 0.576585,
3649          &   0.44445, 0.324385, 0.226815, 0.149165, 0.089375,
3650          &   0.045865, 0.017485, 0.00348, 0, 0, 0, 0, 0,
3651          &   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 /
3652     
3653           data p0 /1013.25 /
3654     
3655     !===>  ...  begin here
3656     
3657     ! --- allocate and initialize gocart climatological data
3658           if ( .not. allocated (dmclmg) ) then
3659             allocate ( dmclmg(IMXG,JMXG,KMXG,NMXG) )
3660             allocate ( psclmg(IMXG,JMXG,KMXG) )
3661             allocate ( molwgt(NMXG) )
3662           endif
3663     
3664           dmclmg(:,:,:,:)  = f_zero
3665           psclmg(:,:,:)    = f_zero
3666           molwgt(:)        = f_zero
3667     
3668     ! --- allocate and initialize geos lat and lon arrays
3669           if ( .not. allocated ( geos_rlon  )) then
3670               allocate (geos_rlon(IMXG))
3671               allocate (geos_rlat(JMXG))
3672           endif
3673     
3674           geos_rlon(:) = f_zero
3675           geos_rlat(:) = f_zero
3676     
3677     ! --- compute geos lat and lon arrays
3678           do i = 1, IMXG
3679              geos_rlon(i)     = -180. + (i-1)* dltx
3680           end do
3681           do j = 2, JMXG-1
3682              geos_rlat(j)     = -90. + (j-1)* dlty
3683           end do
3684           geos_rlat(1)      = -89.5
3685           geos_rlat(JMXG)   =  89.5
3686     
3687     ! --- determine whether GEOS3 or GEOS4 data set is provided
3688           if ( gocart_climo == 'xxxx' ) then
3689             gocart_climo='0000'
3690     ! check geos3-gocart climo
3691             aerosol_file = '200001.PS.avg'
3692             inquire (file = aerosol_file, exist = file_exist)
3693             if ( file_exist ) gocart_climo='ver3'
3694     ! check geos4-gocart climo
3695             aerosol_file = 'gocart_climo_2000x2007_ps_01.bin'
3696             inquire (file = aerosol_file, exist = file_exist)
3697             if ( file_exist ) gocart_climo='ver4'
3698           endif
3699     !
3700     !
3701     ! --- read ps (sfc pressure) and compute 3d pressure field (psclmg)
3702     !
3703           write(mn,'(i2.2)') imon
3704           ymd = yr//mn
3705           aerosol_file = 'null'
3706           if ( gocart_climo == 'ver3' ) then
3707             aerosol_file = ymd//'.PS.avg'
3708           elseif ( gocart_climo == 'ver4' ) then
3709             aerosol_file = 'gocart_climo_2000x2007_ps_'//mn//'.bin'
3710           endif
3711     !
3712           inquire (file = aerosol_file, exist = file_exist)
3713           lab_if_ps : if ( file_exist ) then
3714     
3715            close(NIAERCM)
3716            if ( gocart_climo == 'ver3' ) then
3717             nrecl = 4 * (IMXG * JMXG)
3718             open(NIAERCM, file=trim(aerosol_file),
3719          &       access='direct',recl=nrecl)
3720             read(NIAERCM, rec=1) ps
3721             do j = 1, JMXG
3722               do i = 1, IMXG
3723                 do k = 1, KMXG
3724                   pstmp = pint + sig(k) * (ps(i,j) - pint)
3725                   psclmg(i,j,k) = 0.1 * pstmp       ! convert mb to cb
3726                 enddo
3727               enddo
3728             enddo
3729     
3730            elseif ( gocart_climo == 'ver4' ) then
3731              open(NIAERCM, file=trim(aerosol_file),
3732          &        status='old', form='unformatted')
3733              read(NIAERCM) ps(:,:)
3734              do j = 1, JMXG
3735                do i = 1, IMXG
3736                  do k = 1, KMXG
3737                   pstmp = hyam(k)*p0 + hybm(k)*ps(i,j)
3738                   psclmg(i,j,k) = 0.1 * pstmp       ! convert mb to cb
3739                 enddo
3740               enddo
3741              enddo
3742     
3743            endif     !  ---- end if_gocart_climo
3744     
3745           else lab_if_ps
3746     
3747             print *,' *** Requested aerosol data file "',
3748          &          trim(aerosol_file),  '" not found!'
3749             print *,' *** Stopped in RD_GOCART_CLIM ! ', me
3750             stop 1006
3751           endif   lab_if_ps
3752     !
3753     ! --- read aerosol dry mass (g/m3) or mixing ratios (mol/mol,kg/kg)
3754     !
3755           lab_do_icmp : do icmp = 1, num_gridcomp
3756     
3757              tp = gridcomp(icmp)
3758     
3759     !        determine aerosol_file
3760              aerosol_file = 'null'
3761              if ( gocart_climo == 'ver3' ) then
3762                if(tp == 'DU')   fname='.DU.STD.tv20.g.avg'
3763                if(tp == 'SS')   fname='.SS.STD.tv17.g.avg'
3764                if(tp == 'SU')   fname='.SU.STD.tv15.g.avg'
3765                if(tp == 'OC')   fname='.CC.STD.tv15.g.avg'
3766                if(tp == 'BC')   fname='.CC.STD.tv15.g.avg'
3767                aerosol_file=ymd//trim(fname)
3768              elseif ( gocart_climo == 'ver4' ) then
3769                fin = 'gocart_climo_2000x2007_'
3770                if(tp == 'DU')   fname=trim(fin)//'du_'
3771                if(tp == 'SS')   fname=trim(fin)//'ss_'
3772                if(tp == 'SU')   fname=trim(fin)//'su_'
3773                if(tp == 'OC')   fname=trim(fin)//'cc_'
3774                if(tp == 'BC')   fname=trim(fin)//'cc_'
3775                aerosol_file=trim(fname)//mn//'.bin'
3776              endif
3777     
3778              numspci = 4
3779              if(tp == 'DU')   numspci = 5
3780              inquire (file=trim(aerosol_file), exist = file_exist)
3781              lab_if_aer: if ( file_exist ) then
3782     !
3783               close(NIAERCM)
3784               if ( gocart_climo == 'ver3' ) then
3785               nrecl = 4 * numspci * (IMXG * JMXG * KMXG + 3)
3786               open (NIAERCM, file=trim(aerosol_file),
3787          &         access='direct', recl=nrecl)
3788               read(NIAERCM,rec=1)(nt1,nt2,nn(i),buff(:,:,:,i),i=1,numspci)
3789     
3790               elseif ( gocart_climo == 'ver4' ) then
3791                open (NIAERCM, file=trim(aerosol_file),
3792          &           status='old', form='unformatted')
3793                do i = 1, numspci
3794                  do k = 1, KMXG
3795                   read(NIAERCM) temp(:,:,k)
3796                   buff(:,:,k,i) =  temp(:,:,k)
3797                  enddo
3798                enddo
3799               endif
3800     
3801     !!===> fill dmclmg with working array buff
3802               select case ( tp )
3803     
3804     ! fill in DU from DU: du1, du2, du3, du4, du5
3805               case ('DU' )
3806                if ( dm_indx%dust1 /= -999) then
3807                 do ii = 1, 5
3808                  dmclmg(:,:,:,dm_indx%dust1+ii-1) = buff(:,:,:,ii)
3809                 enddo
3810                else
3811                 print *, 'ERROR: invalid DU index, abort! ',me
3812                 stop 1007
3813                endif
3814     
3815     ! fill in BC from CC: bc_phobic, oc_phobic, bc_philic, oc_philic
3816               case ('BC' )
3817                if ( dm_indx%soot_phobic /= -999) then
3818                 dmclmg(:,:,:,dm_indx%soot_phobic)=buff(:,:,:,1)
3819                 dmclmg(:,:,:,dm_indx%soot_philic)=buff(:,:,:,3)
3820                 molwgt(dm_indx%soot_phobic) = 12.
3821                 molwgt(dm_indx%soot_philic) = 12.
3822                else
3823                 print *, 'ERROR: invalid BC index, abort! ',me
3824                 stop 1008
3825                endif
3826     
3827     ! fill in SU from SU: dms, so2, so4, msa
3828               case ('SU' )
3829                if ( dm_indx%suso /= -999) then
3830                 dmclmg(:,:,:,dm_indx%suso) = buff(:,:,:,3)
3831                 molwgt(dm_indx%suso) = 96.
3832                else
3833                 print *, 'ERROR: invalid SU index, abort! ',me
3834                 stop 1009
3835                endif
3836     
3837     ! fill in OC from CC: bc_phobic, oc_phobic, bc_philic, oc_philic
3838               case ('OC' )
3839                if ( dm_indx%waso_phobic /= -999) then
3840                 dmclmg(:,:,:,dm_indx%waso_phobic) =
3841          &                           1.4*buff(:,:,:,2)
3842                 dmclmg(:,:,:,dm_indx%waso_philic) =
3843          &                           1.4*buff(:,:,:,4)
3844                 molwgt(dm_indx%waso_phobic) = 12.
3845                 molwgt(dm_indx%waso_philic) = 12.
3846                else
3847                 print *, 'ERROR: invalid OC index, abort! ',me
3848                 stop 1010
3849                endif
3850     
3851     ! fill in SS from SS: ss1, ss2, ss3, ss4
3852               case ('SS' )
3853                if ( dm_indx%ssam /= -999) then
3854                 dmclmg(:,:,:,dm_indx%ssam) = buff(:,:,:,1)
3855                 dmclmg(:,:,:,dm_indx%sscm) = buff(:,:,:,2) +
3856          +                     buff(:,:,:,3)+buff(:,:,:,4)
3857                else
3858                 print *, 'ERROR: invalid SS index, abort! ',me
3859                 stop  1011
3860                endif
3861     
3862               case default
3863     
3864                 print *, 'ERROR: invalid aerosol species, abort ',tp
3865                 stop  1012
3866     
3867               end select
3868     
3869              else   lab_if_aer
3870               print *,' *** Requested aerosol data file "',aerosol_file,
3871          &            '" not found!'
3872               print *,' *** Stopped in RD_GOCART_CLIM ! ', me
3873               stop 1013
3874              endif  lab_if_aer
3875     
3876            enddo lab_do_icmp
3877     
3878           return
3879     !...................................
3880           end subroutine rd_gocart_clim
3881     !-----------------------------------
3882     !
3883     !...................................
3884           end subroutine gocart_init
3885     !-----------------------------------
3886     
3887     !-----------------------------------
3888           subroutine setgocartaer                                           &
3889     !...................................
3890     
3891     !  ---  inputs:
3892          &     ( alon,alat,prslk,rhlay,dz,hz,NBDSWLW,                       &
3893          &       prsl,tlay,qlay,ozlay,                                      &
3894          &       IMAX,NLAY,NLP1, iflip, lsswr,lslwr,                        &
3895     !  ---  outputs:
3896          &       aerosw,aerolw                                              &
3897          &,      tau_gocart                                                 &
3898          &     )
3899     
3900     
3901     !  ==================================================================  !
3902     !                                                                      !
3903     !  setgocartaer computes sw + lw aerosol optical properties for gocart !
3904     !  aerosol species (merged from fcst and clim fields)                  !
3905     !                                                                      !
3906     !  inputs:                                                             !
3907     !     alon, alat                                             IMAX      !
3908     !             - longitude and latitude of given points in degree       !
3909     !     prslk   - pressure                           cb   IMAX*NLAY      !
3910     !     rhlay   - layer mean relative humidity            IMAX*NLAY      !
3911     !     dz      - layer thickness                    m    IMAX*NLAY      !
3912     !     hz      - level high                         m    IMAX*NLP1      !
3913     !     NBDSWLW - total number of sw+ir bands for aeros opt prop  1      !
3914     !     prsl    - layer mean pressure                mb   IMAX*NLAY      !
3915     !     tlay    - layer mean temperature             k    IMAX*NLAY      !
3916     !     qlay    - layer mean specific humidity       g/g  IMAX*NLAY      !
3917     !     ozlay   - layer mean specific tracer         g/g  IMAX*NLAY*NTRAC!
3918     !     IMAX    - horizontal dimension of arrays                  1      !
3919     !     NLAY,NLP1-vertical dimensions of arrays                   1      !
3920     !     iflip   - control flag for direction of vertical index    1      !
3921     !               =0: index from toa to surface                          !
3922     !               =1: index from surface to toa                          !
3923     !     lsswr,lslwr                                                      !
3924     !             - logical flag for sw/lw radiation calls          1      !
3925     !                                                                      !
3926     !  outputs:                                                            !
3927     !     aerosw - aeros opt properties for sw      IMAX*NLAY*NBDSW*NF_AESW!
3928     !               (:,:,:,1): optical depth                               !
3929     !               (:,:,:,2): single scattering albedo                    !
3930     !               (:,:,:,3): asymmetry parameter                         !
3931     !     aerolw - aeros opt properties for lw      IMAX*NLAY*NBDLW*NF_AELW!
3932     !               (:,:,:,1): optical depth                               !
3933     !               (:,:,:,2): single scattering albedo                    !
3934     !               (:,:,:,3): asymmetry parameter                         !
3935     !     tau_gocart - 550nm aeros opt depth     IMAX*NLAY*MAX_NUM_GRIDCOMP!
3936     !                                                                      !
3937     !  module parameters and constants:                                    !
3938     !     NBDSW   - total number of sw bands for aeros opt prop     1      !
3939     !     NBDIR   - total number of ir bands for aeros opt prop     1      !
3940     !                                                                      !
3941     !  module variable: (set by subroutine gocart_init)                    !
3942     !     dmclmg  - aerosols dry mass/mixing ratios   IMXG*JMXG*KMXG*NMXG  !
3943     !     psclmg  - pressure                cb        IMXG*JMXG*KMXG       !
3944     !                                                                      !
3945     !  usage:    call setgocartaer                                         !
3946     !                                                                      !
3947     !  subprograms called:  map_aermr, aeropt_grt                          !
3948     !                                                                      !
3949     !  ==================================================================  !
3950     !
3951           implicit none
3952     
3953     !  ---  inputs:
3954           integer, intent(in) :: IMAX,NLAY,NLP1,iflip,NBDSWLW
3955           logical, intent(in) :: lsswr, lslwr
3956     
3957           real (kind=kind_phys), dimension(:,:), intent(in) :: prslk,       &
3958          &       prsl, rhlay, tlay, qlay, dz, hz
3959           real (kind=kind_phys), dimension(:),   intent(in) :: alon, alat
3960           real (kind=kind_phys), dimension(:,:,:), intent(in) :: ozlay
3961     
3962     !  ---  outputs:
3963           real (kind=kind_phys), dimension(:,:,:,:), intent(out) ::         &
3964          &       aerosw, aerolw
3965           real (kind=kind_phys), dimension(:,:,:), intent(out) :: tau_gocart
3966     
3967     !  ---  locals:
3968           real (kind=kind_phys), dimension(NLAY) :: rh1, dz1
3969           real (kind=kind_phys), dimension(NLAY,NBDSWLW)::tauae,ssaae,asyae
3970           real (kind=kind_phys), dimension(NLAY,max_num_gridcomp) ::
3971          &                       tauae_gocart
3972     
3973           real (kind=kind_phys) :: tmp1, tmp2
3974     
3975           integer               :: i, i1, i2, j1, j2, k, m, m1, kp
3976     
3977     !     prognostic aerosols on gfs grids
3978           real (kind=kind_phys), dimension(:,:,:),allocatable:: aermr,dmfcs
3979     
3980     ! aerosol (dry mass) on gfs grids/levels
3981           real (kind=kind_phys), dimension(:,:), allocatable :: dmanl,dmclm,
3982          &                                                      dmclmx
3983           real (kind=kind_phys), dimension(KMXG)     :: pstmp, pkstr
3984           real (kind=kind_phys) :: ptop, psfc, tem, plv, tv, rho
3985     
3986     !  ---  conversion constants
3987           real (kind=kind_phys), parameter :: hdltx = 0.5 * dltx
3988           real (kind=kind_phys), parameter :: hdlty = 0.5 * dlty
3989     
3990     !===>  ...  begin here
3991     !
3992           if ( .not. allocated(dmanl) ) then
3993             allocate ( dmclmx(KMXG,NMXG) )
3994             allocate ( dmanl(NLAY,NMXG) )
3995             allocate ( dmclm(NLAY,NMXG) )
3996     
3997             allocate ( aermr(IMAX,NLAY,NMXG) )
3998             allocate ( dmfcs(IMAX,NLAY,NMXG) )
3999           endif
4000     !
4001     ! --  map input tracer array (ozlay) to local tracer array (aermr)
4002     !
4003           dmfcs(:,:,:) = f_zero
4004           lab_if_fcst : if ( get_fcst ) then
4005     
4006             call map_aermr
4007     !  ---  inputs:  (in scope variables)
4008     !  ---  outputs: (in scope variables)
4009     
4010           endif       lab_if_fcst
4011     !
4012     ! --  map geos-gocart climo (dmclmg) to gfs grids (dmclm)
4013     !
4014           lab_do_IMAX : do i = 1, IMAX
4015     
4016             dmclm(:,:) = f_zero
4017     
4018             lab_if_clim : if ( get_clim ) then
4019     !  ---  map grid in longitude direction
4020               i2 = 1
4021               j2 = 1
4022               tmp1 = alon(i)
4023               if (tmp1 > 180.) tmp1 = tmp1 - 360.0
4024               lab_do_IMXG : do i1 = 1, IMXG
4025                 tmp2 = geos_rlon(i1)
4026                 if (tmp2 > 180.) tmp2 = tmp2 - 360.0
4027                 if (abs(tmp1-tmp2) <= hdltx) then
4028                   i2 = i1
4029                   exit lab_do_IMXG
4030                 endif
4031               enddo  lab_do_IMXG
4032     
4033     !  ---  map grid in latitude direction
4034               lab_do_JMXG : do j1 = 1, JMXG
4035                 if (abs(alat(i)-geos_rlat(j1)) <= hdlty) then
4036                   j2 = j1
4037                   exit lab_do_JMXG
4038                 endif
4039               enddo  lab_do_JMXG
4040     
4041     !  ---  update local arrays pstmp and dmclmx
4042               pstmp(:)= psclmg(i2,j2,:)*1000.0      ! cb to Pa
4043               dmclmx(:,:) = dmclmg(i2,j2,:,:)
4044     
4045     !  ---  map geos-gocart climo (dmclmx) to gfs level (dmclm)
4046               pkstr(:)=fpkap(pstmp(:))
4047               psfc = pkstr(1)                       ! pressure at sfc
4048               ptop = pkstr(KMXG)                    ! pressure at toa
4049     
4050     !  ---  map grid in verical direction (follow how ozone is mapped
4051     !       in radiation_gases routine)
4052               do k = 1, NLAY
4053                kp = k                              ! from sfc to toa
4054                if(iflip==0) kp = NLAY - k + 1      ! from toa to sfc
4055                tmp1 = prslk(i,kp)
4056     
4057                do m1 = 1, KMXG - 1                 ! from sfc to toa
4058                  if(tmp1 > pkstr(m1+1) .and. tmp1 <= pkstr(m1)) then
4059                    tmp2 = f_one / (pkstr(m1)-pkstr(m1+1))
4060                    tem = (pkstr(m1) - tmp1) * tmp2
4061                    dmclm(kp,:) = tem * dmclmx(m1+1,:)+
4062          &                   (f_one-tem) * dmclmx(m1,:)
4063                  endif
4064                enddo
4065     
4066     !*         if(tmp1 > psfc) dmclm(kp,:) = dmclmx(1,:)
4067     !*         if(tmp1 < ptop) dmclm(kp,:) = dmclmx(KMXG,:)
4068     
4069               enddo
4070             endif    lab_if_clim
4071     !
4072     !  ---  compute fcst/clim merged aerosol loading (dmanl) and the
4073     !       radiation optical properties (aerosw, aerolw)
4074     !
4075             do k = 1, NLAY
4076     
4077     !  ---  map global to local arrays (rh1 and dz1)
4078               rh1(k) = rhlay(i,k)
4079               dz1(k) = dz   (i,k)
4080     
4081     !  ---  convert from mixing ratio to dry mass (g/m3)
4082               plv = 100. * prsl(i,k)       ! convert pressure from mb to Pa
4083               tv  = tlay(i,k) * (f_one+con_fvirt*qlay(i,k))  ! virtual temp in K
4084               rho = plv / (con_rd * tv)    ! air density in kg/m3
4085               if ( get_fcst ) then
4086                 do m = 1,  NMXG            ! mixing ratio (g/g)
4087                 dmfcs(i,k,m) = max(1000.*(rho*aermr(i,k,m)),f_zero)
4088                 enddo     ! m_do_loop
4089               endif
4090               if ( get_clim .and. (gocart_climo == 'ver4') ) then
4091                 do m = 1,  NMXG
4092                  dmclm(k,m)=1000.*dmclm(k,m)*rho !mixing ratio (g/g)
4093                  if ( molwgt(m) /= 0. ) then     !mixing ratio (mol/mol)
4094                   dmclm(k,m)=dmclm(k,m) * (molwgt(m)/con_amd)
4095                  endif
4096                 enddo     ! m_do_loop
4097               endif
4098     
4099     
4100     !  ---  determine dmanl from dmclm and dmfcs
4101               do m = 1, NMXG
4102                  dmanl(k,m)= ctaer*dmfcs(i,k,m) +
4103          &                 ( f_one-ctaer)*dmclm(k,m)
4104               enddo
4105             enddo
4106     
4107     !  ---  calculate sw/lw aerosol optical properties for the
4108     !       corresponding frequency bands
4109     
4110             call aeropt_grt                                                 &
4111     !  ---  inputs:  (in scope variables)
4112     !  ---  outputs: (in scope variables)
4113     
4114             if ( lsswr ) then
4115     
4116               if ( laswflg ) then
4117     
4118                 do m = 1, NBDSW
4119                   do k = 1, NLAY
4120                     aerosw(i,k,m,1) = tauae(k,m)
4121                     aerosw(i,k,m,2) = ssaae(k,m)
4122                     aerosw(i,k,m,3) = asyae(k,m)
4123                   enddo
4124                 enddo
4125     !
4126                 do k = 1, NLAY
4127                   do m = 1, max_num_gridcomp
4128                    tau_gocart(i,k,m) = tauae_gocart(k,m)
4129                   enddo
4130                 enddo
4131     
4132               else
4133     
4134                 aerosw(:,:,:,:) = f_zero
4135     
4136               endif
4137     
4138             endif     ! end if_lsswr_block
4139     
4140             if ( lslwr ) then
4141     
4142               if ( lalwflg ) then
4143     
4144                 if ( NBDIR == 1 ) then
4145                   m1 = NBDSW + 1
4146                   do m = 1, NBDLW
4147                     do k = 1, NLAY
4148                       aerolw(i,k,m,1) = tauae(k,m1)
4149                       aerolw(i,k,m,2) = ssaae(k,m1)
4150                       aerolw(i,k,m,3) = asyae(k,m1)
4151                     enddo
4152                   enddo
4153                 else
4154                   do m = 1, NBDLW
4155                     m1 = NBDSW + m
4156                     do k = 1, NLAY
4157                       aerolw(i,k,m,1) = tauae(k,m1)
4158                       aerolw(i,k,m,2) = ssaae(k,m1)
4159                       aerolw(i,k,m,3) = asyae(k,m1)
4160                     enddo
4161                   enddo
4162                 endif
4163     
4164               else
4165     
4166                 aerolw(:,:,:,:) = f_zero
4167     
4168               endif
4169             endif     ! end if_lslwr_block
4170     
4171           enddo  lab_do_IMAX
4172     
4173     ! =================
4174           contains
4175     ! =================
4176     
4177     !-----------------------------
4178           subroutine map_aermr
4179     !.............................
4180     !  ---  inputs:  (in scope variables)
4181     !  ---  outputs: (in scope variables)
4182     
4183     ! ==================================================================== !
4184     !                                                                      !
4185     ! subprogram: map_aermr                                                !
4186     !                                                                      !
4187     !   map input tracer fields (ozlay) to local tracer array (aermr)      !
4188     !                                                                      !
4189     !  ====================  defination of variables  ===================  !
4190     !                                                                      !
4191     ! input arguments:                                                     !
4192     !     IMAX    - horizontal dimension of arrays                  1      !
4193     !     NLAY    - vertical dimensions of arrays                   1      !
4194     !     ozlay   - layer tracer mass mixing ratio     g/g  IMAX*NLAY*NTRAC!
4195     ! output arguments: (to module variables)                              !
4196     !     aermr   - layer aerosol mass mixing ratio    g/g  IMAX*NLAY*NMXG !
4197     !                                                                      !
4198     ! note:                                                                !
4199     !  NTRAC is the number of tracers excluding water vapor                !
4200     !  NMXG is the number of prognostic aerosol species                    !
4201     !  ==================================================================  !
4202     !
4203           implicit none
4204     
4205     !  ---  inputs:
4206     !  ---  output:
4207     
4208     !  ---  local:
4209           integer    :: i, indx, ii
4210           character  :: tp*2
4211     
4212     ! initialize
4213           aermr(:,:,:) = f_zero
4214           ii = 1        !! <---- ozlay does not contain q
4215     
4216     ! ==>  DU: du1 (submicron bins), du2, du3, du4, du5
4217            if( gfs_phy_tracer%doing_DU ) then
4218              aermr(:,:,dm_indx%dust1) = ozlay(:,:,dmfcs_indx%du001-ii)
4219              aermr(:,:,dm_indx%dust2) = ozlay(:,:,dmfcs_indx%du002-ii)
4220              aermr(:,:,dm_indx%dust3) = ozlay(:,:,dmfcs_indx%du003-ii)
4221              aermr(:,:,dm_indx%dust4) = ozlay(:,:,dmfcs_indx%du004-ii)
4222              aermr(:,:,dm_indx%dust5) = ozlay(:,:,dmfcs_indx%du005-ii)
4223            endif
4224     
4225     ! ==>  OC: oc_phobic, oc_philic
4226            if( gfs_phy_tracer%doing_OC ) then
4227              aermr(:,:,dm_indx%waso_phobic) =
4228          &                     ozlay(:,:,dmfcs_indx%ocphobic-ii)
4229              aermr(:,:,dm_indx%waso_philic) =
4230          &                     ozlay(:,:,dmfcs_indx%ocphilic-ii)
4231            endif
4232     
4233     ! ==>  BC: bc_phobic, bc_philic
4234            if( gfs_phy_tracer%doing_BC ) then
4235              aermr(:,:,dm_indx%soot_phobic) =
4236          &                     ozlay(:,:,dmfcs_indx%bcphobic-ii)
4237              aermr(:,:,dm_indx%soot_philic) =
4238          &                     ozlay(:,:,dmfcs_indx%bcphilic-ii)
4239            endif
4240     
4241     ! ==>  SS: ss1, ss2 (submicron bins), ss3, ss4, ss5
4242            if( gfs_phy_tracer%doing_SS ) then
4243               aermr(:,:,dm_indx%ssam)=ozlay(:,:,dmfcs_indx%ss001-ii)
4244          &                          + ozlay(:,:,dmfcs_indx%ss002-ii)
4245               aermr(:,:,dm_indx%sscm)=ozlay(:,:,dmfcs_indx%ss003-ii)
4246          &                          + ozlay(:,:,dmfcs_indx%ss004-ii)
4247          &                          + ozlay(:,:,dmfcs_indx%ss005-ii)
4248            endif
4249     
4250     ! ==>  SU: so4
4251            if( gfs_phy_tracer%doing_SU ) then
4252               aermr(:,:,dm_indx%suso) = ozlay(:,:,dmfcs_indx%so4-ii)
4253            endif
4254     
4255           return
4256     !...................................
4257           end subroutine map_aermr
4258     !-----------------------------------
4259     
4260     
4261     !-----------------------------------
4262           subroutine aeropt_grt                                             &
4263     !...................................
4264     !  ---  inputs:  (in scope variables)
4265     !  ---  outputs: (in scope variables)
4266     
4267     !  ==================================================================  !
4268     !                                                                      !
4269     !  subprogram: aeropt_grt                                              !
4270     !                                                                      !
4271     !  compute aerosols optical properties in NBDSWLW sw/lw bands.         !
4272     !  Aerosol distribution at each grid point is composed from up to      !
4273     !  NMXG aerosol species (from NUM_GRIDCOMP components).                !
4274     !                                                                      !
4275     !  input variables:                                                    !
4276     !     dmanl  - aerosol dry mass                     g/m3   NLAY*NMXG   !
4277     !     rh1    - relative humidity                     %     NLAY        !
4278     !     dz1    - layer thickness                       km    NLAY        !
4279     !     NLAY   - vertical dimensions                   -     1           !
4280     !     iflip  - control flag for direction of vertical index            !
4281     !               =0: index from toa to surface                          !
4282     !               =1: index from surface to toa                          !
4283     !                                                                      !
4284     !  output variables:                                                   !
4285     !     tauae  - aerosol optical depth                 -   NLAY*NBDSWLW  !
4286     !     ssaae  - aerosol single scattering albedo      -   NLAY*NBDSWLW  !
4287     !     asyae  - aerosol asymmetry parameter           -   NLAY*NBDSWLW  !
4288     !                                                                      !
4289     !  ==================================================================  !
4290     !
4291           implicit none
4292     
4293     !  ---  inputs:
4294     !  ---  outputs:
4295     
4296     !  ---  locals:
4297           real (kind=kind_phys) :: aerdm
4298           real (kind=kind_phys) :: ext1, ssa1, asy1, ex00, ss00, as00,
4299          &                         ex01, ss01, as01, exint
4300           real (kind=kind_phys) :: tau, ssa, asy,
4301          &                         sum_tau, sum_ssa, sum_asy
4302     
4303     !  ---  subgroups for sub-micron dust
4304     !  ---  corresponds to 0.1-0.18, 0.18-0.3, 0.3-0.6, 0.6-1.0 micron
4305     
4306           real (kind=kind_phys) ::   fd(4)
4307           data fd  / 0.01053,0.08421,0.25263,0.65263 /
4308     
4309           character    :: tp*2
4310           integer      :: icmp, n, kk, ib, ih2, ih1, ii, ij, ijk
4311           real (kind=kind_phys) :: drh0, drh1, rdrh
4312     
4313           real (kind=kind_phys) :: qmin  !<--lower bound for opt calc
4314           data qmin  / 1.e-20 /
4315     
4316     !===>  ...  begin here
4317     
4318     ! --- initialize (assume no aerosols)
4319           tauae = f_zero
4320           ssaae = f_one
4321           asyae = f_zero
4322     
4323           tauae_gocart = f_zero
4324     
4325     !===> ... loop over vertical layers
4326     !
4327           lab_do_layer : do kk = 1, NLAY
4328     
4329     ! --- linear interp coeffs for rh-dep species
4330     
4331             ih2 = 1
4332             do while ( rh1(kk) > rhlev_grt(ih2) )
4333               ih2 = ih2 + 1
4334               if ( ih2 > KRHLEV ) exit
4335             enddo
4336             ih1 = max( 1, ih2-1 )
4337             ih2 = min( KRHLEV, ih2 )
4338     
4339             drh0 = rhlev_grt(ih2) - rhlev_grt(ih1)
4340             drh1 = rh1(kk) - rhlev_grt(ih1)
4341             if ( ih1 == ih2 ) then
4342               rdrh = f_zero
4343             else
4344               rdrh = drh1 / drh0
4345             endif
4346     
4347     ! --- loop through sw/lw spectral bands
4348     
4349             lab_do_ib : do ib = 1, NBDSWLW
4350               sum_tau = f_zero
4351               sum_ssa = f_zero
4352               sum_asy = f_zero
4353     
4354     ! --- loop through aerosol grid components
4355               lab_do_icmp : do icmp = 1, NUM_GRIDCOMP
4356                 ext1 = f_zero
4357                 ssa1 = f_zero
4358                 asy1 = f_zero
4359     
4360                 tp = gridcomp(icmp)
4361     
4362                 select case ( tp )
4363     
4364     ! -- dust aerosols: no humidification effect
4365                 case ( 'DU')
4366                   do n = 1, KCM1
4367     
4368                     if (n <= 4) then
4369                       aerdm = dmanl(kk,dm_indx%dust1) * fd(n)
4370                     else
4371                       aerdm = dmanl(kk,dm_indx%dust1+n-4 )
4372                     endif
4373     
4374                     if (aerdm < qmin) aerdm = f_zero
4375                     ex00 = extrhi_grt(n,ib)*(1000.*dz1(kk))*aerdm
4376                     ss00 = ssarhi_grt(n,ib)
4377                     as00 = asyrhi_grt(n,ib)
4378                     ext1 = ext1 + ex00
4379                     ssa1 = ssa1 + ex00 * ss00
4380                     asy1 = asy1 + ex00 * ss00 * as00
4381     
4382                   enddo
4383     
4384     ! -- suso aerosols: with humidification effect
4385                 case ( 'SU')
4386                   ij = isuso
4387                   exint = extrhd_grt(ih1,ij,ib)
4388          &          + rdrh*(extrhd_grt(ih2,ij,ib) - extrhd_grt(ih1,ij,ib))
4389                   ss00 = ssarhd_grt(ih1,ij,ib)
4390          &          + rdrh*(ssarhd_grt(ih2,ij,ib) - ssarhd_grt(ih1,ij,ib))
4391                   as00 = asyrhd_grt(ih1,ij,ib)
4392          &          + rdrh*(asyrhd_grt(ih2,ij,ib) - asyrhd_grt(ih1,ij,ib))
4393     
4394                   aerdm = dmanl(kk, dm_indx%suso)
4395                   if (aerdm < qmin) aerdm = f_zero
4396                   ex00 = exint*(1000.*dz1(kk))*aerdm
4397                   ext1 = ex00
4398                   ssa1 = ex00 * ss00
4399                   asy1 = ex00 * ss00 * as00
4400     
4401     ! -- seasalt aerosols: with humidification effect
4402                 case ( 'SS')
4403                   do n = 1, 2                  !<---- ssam, sscm
4404                     ij = issam + (n-1)
4405                     exint = extrhd_grt(ih1,ij,ib)
4406          &          + rdrh*(extrhd_grt(ih2,ij,ib) - extrhd_grt(ih1,ij,ib))
4407                     ss00 = ssarhd_grt(ih1,ij,ib)
4408          &          + rdrh*(ssarhd_grt(ih2,ij,ib) - ssarhd_grt(ih1,ij,ib))
4409                     as00 = asyrhd_grt(ih1,ij,ib)
4410          &          + rdrh*(asyrhd_grt(ih2,ij,ib) - asyrhd_grt(ih1,ij,ib))
4411     
4412                     aerdm = dmanl(kk, dm_indx%ssam+n-1)
4413                     if (aerdm < qmin) aerdm = f_zero
4414                     ex00 = exint*(1000.*dz1(kk))*aerdm
4415                     ext1 = ext1 + ex00
4416                     ssa1 = ssa1 + ex00 * ss00
4417                     asy1 = asy1 + ex00 * ss00 * as00
4418     
4419                   enddo
4420     
4421     ! -- organic carbon/black carbon:
4422     !    using 'waso' and 'soot' for hydrophilic OC and BC
4423     !    using 'waso' and 'soot' at RH=0 for hydrophobic OC and BC
4424                 case ( 'OC', 'BC')
4425                   if(tp == 'OC') then
4426                      ii = dm_indx%waso_phobic
4427                      ij = iwaso
4428                   else
4429                      ii = dm_indx%soot_phobic
4430                      ij = isoot
4431                   endif
4432     
4433     ! --- hydrophobic
4434                   aerdm = dmanl(kk, ii)
4435                   if (aerdm < qmin) aerdm = f_zero
4436                   ex00 = extrhd_grt(1,ij,ib)*(1000.*dz1(kk))*aerdm
4437                   ss00 = ssarhd_grt(1,ij,ib)
4438                   as00 = asyrhd_grt(1,ij,ib)
4439     ! --- hydrophilic
4440                   aerdm = dmanl(kk, ii+1)
4441                   if (aerdm < qmin) aerdm = f_zero
4442                   exint = extrhd_grt(ih1,ij,ib)
4443          &         + rdrh*(extrhd_grt(ih2,ij,ib) - extrhd_grt(ih1,ij,ib))
4444                   ex01 = exint*(1000.*dz1(kk))*aerdm
4445                   ss01 = ssarhd_grt(ih1,ij,ib)
4446          &          + rdrh*(ssarhd_grt(ih2,ij,ib) - ssarhd_grt(ih1,ij,ib))
4447                   as01 = asyrhd_grt(ih1,ij,ib)
4448          &          + rdrh*(asyrhd_grt(ih2,ij,ib) - asyrhd_grt(ih1,ij,ib))
4449     
4450                   ext1 = ex00 + ex01
4451                   ssa1 = (ex00 * ss00) + (ex01 * ss01)
4452                   asy1 = (ex00 * ss00 * as00) + (ex01 * ss01 * as01)
4453     
4454                 end select
4455     
4456     ! --- determine tau, ssa, asy for each grid component
4457                 tau = ext1
4458                 if (ext1 > f_zero) ssa=min(f_one,ssa1/ext1)
4459                 if (ssa1 > f_zero) asy=min(f_one,asy1/ssa1)
4460     
4461     ! --- save tau at 550 nm for each grid component
4462                 if ( ib == nv_aod ) then
4463                   do ijk = 1, max_num_gridcomp
4464                     if ( tp == max_gridcomp(ijk) )  then
4465                        tauae_gocart(kk,ijk) = tau
4466                     endif
4467                   enddo
4468                 endif
4469     
4470     ! --- update sum_tau, sum_ssa, sum_asy
4471                 sum_tau = sum_tau + tau
4472                 sum_ssa = sum_ssa + tau * ssa
4473                 sum_asy = sum_asy + tau * ssa * asy
4474     
4475               enddo lab_do_icmp
4476     
4477     
4478     ! --- determine total tau, ssa, asy for aerosol mixture
4479               tauae(kk,ib) = sum_tau
4480               if (sum_tau > f_zero) ssaae(kk,ib) = sum_ssa / sum_tau
4481               if (sum_ssa > f_zero) asyae(kk,ib) = sum_asy / sum_ssa
4482     
4483             enddo   lab_do_ib
4484     
4485           enddo lab_do_layer
4486     
4487     !
4488           return
4489     !...................................
4490           end subroutine aeropt_grt
4491     !--------------------------------
4492     
4493     !................................
4494           end subroutine setgocartaer
4495     !--------------------------------
4496     !
4497     ! GOCART code modification end here (Sarah Lu)  ------------------------!
4498     ! =======================================================================
4499     
4500     
4501     !..........................................!
4502           end module module_radiation_aerosols !
4503     !==========================================!
4504