CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
radiation_aerosols.f
1
4
5! ========================================================== !!!!!
6! 'module_radiation_aerosols' description !!!!!
7! ========================================================== !!!!!
8! !
9! this module contains climatological atmospheric aerosol schemes for!
10! radiation computations. !
11! !
12! in the module, the externally callable subroutines are : !
13! !
14! 'aer_init' -- initialization !
15! inputs: !
16! ( NLAY, me ) !
17! outputs: !
18! ( errflg, errmsg ) !
19! !
20! 'aer_update' -- updating aerosol data !
21! inputs: !
22! ( iyear, imon, me ) !
23! outputs: !
24! ( errflg, errmsg ) !
25! !
26! 'setaer' -- mapping aeros profile, compute aeros opticals !
27! inputs: !
28! (prsi,prsl,prslk,tvly,rhlay,slmsk,tracer,aerfld,xlon,xlat, !
29! IMAX,NLAY,NLP1, lsswr,lslwr, !
30! outputs: !
31! (aerosw,aerolw,aerodp,errmsg,errflg) !
32! !
33! !
34! external modules referenced: !
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 asymmetry 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! oct 23, 2010 --- Hsin-mu lin modified subr setclimaer to !
68! interpolate the 5 degree aerosol data to small domain based on!
69! the nearby 4 points instead of previous nearby assignment by !
70! using the 5 degree data. this process will eliminate the dsw !
71! jagged edges in the east conus where aerosol effect are lagre.!
72! dec 2010 --- y.-t. hou modified and optimized bi-linear!
73! horizontal interpolation in subr setclimaer. added safe guard !
74! measures in lat/lon indexing and added sea/land mask variable !
75! slmsk as input field to help aerosol profile selection. !
76! jan 2011 --- y.-t. hou divided the program into two !
77! separated interchangeable modules: a climatology aerosol !
78! module, and a gocart aerosol scheme module. the stratospheric !
79! volcanic aerosol part is still within the two driver modules, !
80! and may also become a separate one in the further development.!
81! unified in/out argument list for both clim and gocart types of!
82! schemes and added vertically integrated aer-opt-dep, aerodp, !
83! to replace tau_gocart as optional output for various species. !
84! aug 2012 --- y.-t. hou changed the initialization subr !
85! 'aerinit' into two parts: 'aer_init' is called at the start !
86! of run to set up module parameters; and 'aer_update' is !
87! called within the time loop to check and update data sets. !
88! nov 2012 --- y.-t. hou modified control parameters thru!
89! module 'physparam'. !
90! jan 2013 --- sarah lu and y.-t. hou reintegrate both !
91! opac-clim and gocart schemes into one module to make the !
92! program best utilize common components. added aerosol model !
93! scheme selection control variable iaer_mdl to the namelist. !
94! Aug 2013 --- s. moorthi - merge sarah's gocart changes with!
95! yutai's changes !
96! 13Feb2014 --- Sarah lu - compute aod at 550nm !
97! jun 2018 --- h-m lin and y-t hou updated spectral band !
98! mapping method for aerosol optical properties. controled by !
99! internal variable lmap_new through namelist variable iaer. !
100! may 2019 --- sarah lu, restore the gocart option, allowing !
101! aerosol ext, ssa, asy determined from MERRA2 monthly climo !
102! with new spectral band mapping method !
103! !
104! references for opac climatological aerosols: !
105! hou et al. 2002 (ncep office note 441) !
106! hess et al. 1998 - bams v79 831-844 !
107! !
108! references for gocart interactive aerosols: !
109! chin et al., 2000 - jgr, v105, 24671-24687 !
110! colarco et al., 2010 - jgr, v115, D14207 !
111! !
112! references for merra2 aerosol reanalysis: !
113! randles et al., 2017 - jclim, v30, 6823-6850 !
114! buchard et al., 2017 - jclim, v30, 6851-6871 !
115! !
116! references for stratosperic volcanical aerosols: !
117! sato et al. 1993 - jgr, v98, d12, 22987-22994 !
118! !
119! !
120!!!!! ========================================================== !!!!!
121!!!!! end descriptions !!!!!
122!!!!! ========================================================== !!!!!
123
124
127 module module_radiation_aerosols
128!
129 use machine, only : kind_phys, kind_io4, kind_io8
130 use module_iounitdef, only : niaercm
131 use module_radsw_parameters, only : nbdsw, wvnsw1=>wvnum1, &
132 & nswstr, wvnsw2=>wvnum2
133 use module_radlw_parameters, only : nbdlw, wvnlw1, wvnlw2
134!
135 use funcphys, only : fpkap
136 use aerclm_def, only : ntrcaerm
137
138!
139 implicit none
140!
141 private
142
143! --- version tag and last revision date
144 character(40), parameter :: &
145 & VTAGAER='NCEP-Radiation_aerosols v5.2 Jan 2013 '
146! & VTAGAER='NCEP-Radiation_aerosols v5.1 Nov 2012 '
147! & VTAGAER='NCEP-Radiation_aerosols v5.0 Aug 2012 '
148
149! --- general use parameter constants:
151 integer, parameter, public :: nf_aesw = 3
153 integer, parameter, public :: nf_aelw = 3
155 integer, parameter, public :: nlwstr = 1
157 integer, parameter, public :: nspc = 5
159 integer, parameter, public :: nspc1 = nspc + 1
160
161 real (kind=kind_phys), parameter :: f_zero = 0.0
162 real (kind=kind_phys), parameter :: f_one = 1.0
163
164! --- module control parameters set in subroutine "aer_init"
167 integer, save :: nswbnd = nbdsw
170 integer, save :: nlwbnd = nbdlw
172 integer, save :: nswlwbd = nbdsw+nbdlw
173! LW aerosols effect control flag
174! =.true.:aerosol effect is included in LW radiation
175! =.false.:aerosol effect is not included in LW radiation
176 logical, save :: lalwflg = .true.
177! SW aerosols effect control flag
178! =.true.:aerosol effect is included in SW radiation
179! =.false.:aerosol effect is not included in SW radiation
180 logical, save :: laswflg = .true.
181! stratospheric volcanic aerosol effect flag
182! =.true.:historical events of stratosphere volcanic aerosol effect
183! is included radiation (limited by data availability)
184! =.false.:volcanic aerosol effect is not included in radiation
185 logical, save :: lavoflg = .true.
186
187 logical, save :: lmap_new = .true. ! use new mapping method (set in aer_init)
188
189! --------------------------------------------------------------------- !
190! section-1 : module variables for spectral band interpolation !
191! similar to gfdl-sw treatment (2000 version) !
192! --------------------------------------------------------------------- !
193
194! --- parameter constants:
196 integer, parameter, public :: nwvsol = 151
197
199 integer, parameter, public :: nwvtot = 57600
201 integer, parameter, public :: nwvtir = 4000
202
204 integer, dimension(NWVSOL), save :: nwvns0
205
206 data nwvns0 / 100, 11, 14, 18, 24, 33, 50, 83, 12, 12, &
207 & 13, 15, 15, 17, 18, 20, 21, 24, 26, 30, 32, 37, 42, &
208 & 47, 55, 64, 76, 91, 111, 139, 179, 238, 333, 41, 42, 45, &
209 & 46, 48, 51, 53, 55, 58, 61, 64, 68, 71, 75, 79, 84, &
210 & 89, 95, 101, 107, 115, 123, 133, 142, 154, 167, 181, 197, 217, &
211 & 238, 263, 293, 326, 368, 417, 476, 549, 641, 758, 909, 101, 103, &
212 & 105, 108, 109, 112, 115, 117, 119, 122, 125, 128, 130, 134, 137, &
213 & 140, 143, 147, 151, 154, 158, 163, 166, 171, 175, 181, 185, 190, &
214 & 196, 201, 207, 213, 219, 227, 233, 240, 248, 256, 264, 274, 282, &
215 & 292, 303, 313, 325, 337, 349, 363, 377, 392, 408, 425, 444, 462, &
216 & 483, 505, 529, 554, 580, 610, 641, 675, 711, 751, 793, 841, 891, &
217 & 947,1008,1075,1150,1231,1323,1425,1538,1667,1633,14300 /
218
220 real (kind=kind_phys), dimension(NWVSOL), save :: s0intv
221
222 data s0intv( 1: 50) / &
223 & 1.60000e-6, 2.88000e-5, 3.60000e-5, 4.59200e-5, 6.13200e-5, &
224 & 8.55000e-5, 1.28600e-4, 2.16000e-4, 2.90580e-4, 3.10184e-4, &
225 & 3.34152e-4, 3.58722e-4, 3.88050e-4, 4.20000e-4, 4.57056e-4, &
226 & 4.96892e-4, 5.45160e-4, 6.00600e-4, 6.53600e-4, 7.25040e-4, &
227 & 7.98660e-4, 9.11200e-4, 1.03680e-3, 1.18440e-3, 1.36682e-3, &
228 & 1.57560e-3, 1.87440e-3, 2.25500e-3, 2.74500e-3, 3.39840e-3, &
229 & 4.34000e-3, 5.75400e-3, 7.74000e-3, 9.53050e-3, 9.90192e-3, &
230 & 1.02874e-2, 1.06803e-2, 1.11366e-2, 1.15830e-2, 1.21088e-2, &
231 & 1.26420e-2, 1.32250e-2, 1.38088e-2, 1.44612e-2, 1.51164e-2, &
232 & 1.58878e-2, 1.66500e-2, 1.75140e-2, 1.84450e-2, 1.94106e-2 /
233 data s0intv( 51:100) / &
234 & 2.04864e-2, 2.17248e-2, 2.30640e-2, 2.44470e-2, 2.59840e-2, &
235 & 2.75940e-2, 2.94138e-2, 3.13950e-2, 3.34800e-2, 3.57696e-2, &
236 & 3.84054e-2, 4.13490e-2, 4.46880e-2, 4.82220e-2, 5.22918e-2, &
237 & 5.70078e-2, 6.19888e-2, 6.54720e-2, 6.69060e-2, 6.81226e-2, &
238 & 6.97788e-2, 7.12668e-2, 7.27100e-2, 7.31610e-2, 7.33471e-2, &
239 & 7.34814e-2, 7.34717e-2, 7.35072e-2, 7.34939e-2, 7.35202e-2, &
240 & 7.33249e-2, 7.31713e-2, 7.35462e-2, 7.36920e-2, 7.23677e-2, &
241 & 7.25023e-2, 7.24258e-2, 7.20766e-2, 7.18284e-2, 7.32757e-2, &
242 & 7.31645e-2, 7.33277e-2, 7.36128e-2, 7.33752e-2, 7.28965e-2, &
243 & 7.24924e-2, 7.23307e-2, 7.21050e-2, 7.12620e-2, 7.10903e-2 /
244 data s0intv(101:151) / 7.12714e-2, &
245 & 7.08012e-2, 7.03752e-2, 7.00350e-2, 6.98639e-2, 6.90690e-2, &
246 & 6.87621e-2, 6.52080e-2, 6.65184e-2, 6.60038e-2, 6.47615e-2, &
247 & 6.44831e-2, 6.37206e-2, 6.24102e-2, 6.18698e-2, 6.06320e-2, &
248 & 5.83498e-2, 5.67028e-2, 5.51232e-2, 5.48645e-2, 5.12340e-2, &
249 & 4.85581e-2, 4.85010e-2, 4.79220e-2, 4.44058e-2, 4.48718e-2, &
250 & 4.29373e-2, 4.15242e-2, 3.81744e-2, 3.16342e-2, 2.99615e-2, &
251 & 2.92740e-2, 2.67484e-2, 1.76904e-2, 1.40049e-2, 1.46224e-2, &
252 & 1.39993e-2, 1.19574e-2, 1.06386e-2, 1.00980e-2, 8.63808e-3, &
253 & 6.52736e-3, 4.99410e-3, 4.39350e-3, 2.21676e-3, 1.33812e-3, &
254 & 1.12320e-3, 5.59000e-4, 3.60000e-4, 2.98080e-4, 7.46294e-5 /
255
256 real (kind=kind_phys), dimension(NBDSW), save :: wvn_sw1, wvn_sw2
257 real (kind=kind_phys), dimension(NBDLW), save :: wvn_lw1, wvn_lw2
258! --------------------------------------------------------------------- !
259! section-2 : module variables for stratospheric volcanic aerosols !
260! from historical data (sato et al. 1993) !
261! --------------------------------------------------------------------- !
262
263! --- parameter constants:
265 integer, parameter :: minvyr = 1850
267 integer, parameter :: maxvyr = 1999
268
270 integer, allocatable, save :: ivolae(:,:,:)
271
272! --- static control variables:
274 integer :: kyrstr
276 integer :: kyrend
278 integer :: kyrsav
280 integer :: kmonsav
281
282! --------------------------------------------------------------------- !
283! section-3 : module variables for opac climatological aerosols !
284! optical properties (hess et al. 1989) !
285! --------------------------------------------------------------------- !
286
287! --- parameters and constants:
289 integer, parameter :: nxc = 5
291 integer, parameter :: nae = 7
293 integer, parameter :: ndm = 5
295 integer, parameter :: imxae = 72
297 integer, parameter :: jmxae = 37
299 integer, parameter :: naerbnd=61
301 integer, parameter :: nrhlev =8
303 integer, parameter :: ncm1 = 6
305 integer, parameter :: ncm2 = 4
306 integer, parameter :: ncm = ncm1+ncm2
307
309 real (kind=kind_phys), dimension(NRHLEV), save :: rhlev
310 data rhlev(:) / 0.0, 0.5, 0.7, 0.8, 0.9, 0.95, 0.98, 0.99 /
311
312! --- the following arrays are for climatological data that are
313! allocated and read in subroutine 'clim_aerinit'.
314! - global aerosol distribution:
315! haer (NDM,NAE) - scale height of aerosols (km)
316! prsref(NDM,NAE) - ref pressure lev (sfc to toa) in mb (100Pa)
317! sigref(NDM,NAE) - ref sigma lev (sfc to toa)
318
320 real (kind=kind_phys), save, dimension(NDM,NAE) :: haer
322 real (kind=kind_phys), save, dimension(NDM,NAE) :: prsref
324 real (kind=kind_phys), save, dimension(NDM,NAE) :: sigref
325
326! --- the following arrays are allocate and setup in subr 'clim_aerinit'
327! - for relative humidity independent aerosol optical properties:
328! species : insoluble (inso); soot (soot);
329! mineral nuc mode (minm); mineral acc mode (miam);
330! mineral coa mode (micm); mineral transport(mitr).
331! extrhi(NCM1,NSWLWBD) - extinction coefficient for sw+lw spectral band
332! scarhi(NCM1,NSWLWBD) - scattering coefficient for sw+lw spectral band
333! ssarhi(NCM1,NSWLWBD) - single scattering albedo for sw+lw spectral band
334! asyrhi(NCM1,NSWLWBD) - asymmetry parameter for sw+lw spectral band
335! - for relative humidity dependent aerosol optical properties:
336! species : water soluble (waso); sea salt acc mode(ssam);
337! sea salt coa mode(sscm); sulfate droplets (suso).
338! rh level: 00%, 50%, 70%, 80%, 90%, 95%, 98%, 99%
339! extrhd(NRHLEV,NCM2,NSWLWBD) - extinction coefficient for sw+lw band
340! scarhd(NRHLEV,NCM2,NSWLWBD) - scattering coefficient for sw+lw band
341! ssarhd(NRHLEV,NCM2,NSWLWBD) - single scattering albedo for sw+lw band
342! asyrhd(NRHLEV,NCM2,NSWLWBD) - asymmetry parameter for sw+lw band
343! - for stratospheric aerosols optical properties:
344! extstra(NSWLWBD) - extinction coefficient for sw+lw band
345
346 real (kind=kind_phys), allocatable, save, dimension(:,:) :: &
347 & extrhi, scarhi, ssarhi, asyrhi
348 real (kind=kind_phys), allocatable, save, dimension(:,:,:) :: &
349 & extrhd, scarhd, ssarhd, asyrhd
350 real (kind=kind_phys), allocatable, save, dimension(:) :: &
351 & extstra
352
353! --- the following arrays are calculated in subr 'clim_aerinit'
354! - for topospheric aerosol profile distibution:
355! kprfg ( IMXAE*JMXAE) - aeros profile index
356! idxcg (NXC*IMXAE*JMXAE) - aeros component index
357! cmixg (NXC*IMXAE*JMXAE) - aeros component mixing ratio
358! denng ( 2 *IMXAE*JMXAE) - aerosols number density
359
361
363 real (kind=kind_phys), dimension(NXC,IMXAE,JMXAE), save :: cmixg
365 real (kind=kind_phys), dimension( 2 ,IMXAE,JMXAE), save :: denng
367 integer, dimension(NXC,IMXAE,JMXAE), save :: idxcg
369 integer, dimension( IMXAE,JMXAE), save :: kprfg
370
371! --------------------------------------------------------------------- !
372! section-4 : module variables for gocart aerosol optical properties !
373! --------------------------------------------------------------------- !
375
376! --- parameters and constants:
378 integer, parameter :: kaerbndd=61
379 integer, parameter :: kaerbndi=56
381 integer, parameter :: krhlev =36
383 integer, parameter :: kcm1 = 5
385 integer, parameter :: kcm2 = 10
387 integer, parameter :: kcm = kcm1 + kcm2
388
389 real (kind=kind_phys), dimension(KRHLEV) :: rhlev_grt &
390 data rhlev_grt(:)/ .00, .05, .10, .15, .20, .25, .30, .35, &
391 & .40, .45, .50, .55, .60, .65, .70, .75, .80, .81, .82, &
392 & .83, .84, .85, .86, .87, .88, .89, .90, .91, .92, .93, &
393 & .94, .95, .96, .97, .98, .99 /
394
397! extrhi_grt(KCM1,NSWLWBD) - extinction coefficient for sw+lw band
398! scarhi_grt(KCM1,NSWLWBD) - scattering coefficient for sw+lw band
399! ssarhi_grt(KCM1,NSWLWBD) - single scattering albedo for sw+lw band
400! asyrhi_grt(KCM1,NSWLWBD) - asymmetry parameter for sw+lw band
401 real (kind=kind_phys),allocatable,save,dimension(:,:) :: &
402 & extrhi_grt, extrhi_grt_550, scarhi_grt, ssarhi_grt, asyrhi_grt
403!
407! extrhd_grt(KRHLEV,KCM2,NSWLWBD) - extinction coefficient for sw+lw band
408! scarhd_grt(KRHLEV,KCM2,NSWLWBD) - scattering coefficient for sw+lw band
409! ssarhd_grt(KRHLEV,KCM2,NSWLWBD) - single scattering albedo for sw+lw band
410! asyrhd_grt(KRHLEV,KCM2,NSWLWBD) - asymmetry parameter for sw+lw band
411
413 real (kind=kind_phys),allocatable,save,dimension(:,:,:) :: &
414 & extrhd_grt, extrhd_grt_550, scarhd_grt, ssarhd_grt, asyrhd_grt
415
417 integer, parameter :: num_gc = 5
418 character*2 :: gridcomp(num_gc)
419 integer, dimension (num_gc):: num_radius, radius_lower
420 integer, dimension (num_gc):: trc_to_aod
421
422 data gridcomp /'DU', 'SS', 'SU', 'BC', 'OC'/
423 data num_radius /5, 5, 1, 2, 2 /
424 data radius_lower /1, 6, 11, 12, 14 /
425 data trc_to_aod /1, 5, 4, 2, 3/ ! dust, soot, waso, suso, ssam
426
427! =======================================================================
428! --------------------------------------------------------------------- !
429! section-5 : module variables for aod diagnostic !
430! --------------------------------------------------------------------- !
431!! --- the following are for diagnostic purpose to output aerosol optical depth
432! aod from 10 components are grouped into 5 major different species:
433! 1:dust (inso,minm,miam,micm,mitr); 2:black carbon (soot)
434! 3:water soluble (waso); 4:sulfate (suso); 5:sea salt (ssam,sscm)
435!
436! idxspc (NCM) - index conversion array
437! lspcaod - logical flag for aod from individual species
438!
440 integer, dimension(NCM) :: idxspc
441 data idxspc / 1, 2, 1, 1, 1, 1, 3, 5, 5, 4 /
442!
443! - wvn550 is the wavenumber (1/cm) of wavelenth 550nm for diagnostic aod output
444! nv_aod is the sw spectral band covering wvn550 (comp in aer_init)
445!
447 real (kind=kind_phys), parameter :: wvn550 = 1.0e4/0.55
449 integer, save :: nv_aod = 1
450 integer :: i550,id550
451
452! --- public interfaces
453
454 public aer_init, aer_update, setaer
455
456! =================
457 contains
458! =================
459
494!-----------------------------------
495 subroutine aer_init &
496 & ( nlay, me, iaermdl, iaerflg, lalw1bd, aeros_file, con_pi, &
497 & con_t0c, con_c, con_boltz, con_plnk, errflg, errmsg)
498
499! ================================================================== !
500! !
501! aer_init is the initialization program to set up necessary !
502! parameters and working arrays. !
503! !
504! inputs: !
505! NLAY - number of model vertical layers (not used) !
506! me - print message control flag !
507! iaermdl - tropospheric aerosol model scheme flag !
508! =0 opac-clim; =1 gocart-clim, =2 gocart-prognostic !
509! =5 opac-clim new spectral mapping !
510! lalw1bd = logical lw aeros propty 1 band vs multi-band cntl flag !
511! =t use 1 broad band optical property !
512! =f use multi bands optical property !
513! !
514! outputs: (CCPP error handling) !
515! errmsg - CCPP error message !
516! errflg - CCPP error flag !
517! !
518! internal module variables: !
519! lalwflg - logical lw aerosols effect control flag !
520! =t compute lw aerosol optical prop !
521! laswflg - logical sw aerosols effect control flag !
522! =t compute sw aerosol optical prop !
523! lavoflg - logical stratosphere volcanic aerosol control flag !
524! =t include volcanic aerosol effect !
525! !
526! internal module constants: !
527! NWVSOL - num of wvnum regions where solar flux is constant !
528! NWVTOT - total num of wave numbers used in sw spectrum !
529! NWVTIR - total num of wave numbers used in the ir region !
530! NSWBND - total number of sw spectral bands !
531! NLWBND - total number of lw spectral bands !
532! !
533! usage: call aer_init !
534! !
535! subprograms called: clim_aerinit, gocart_aerinit, !
536! wrt_aerlog, set_volcaer, set_spectrum, !
537! !
538! ================================================================== !
539
540! --- inputs:
541 integer, intent(in) :: nlay, me, iaermdl, iaerflg
542 logical, intent(in) :: lalw1bd
543 character(len=26),intent(in) :: aeros_file
544 real(kind_phys), intent(in) :: con_pi,con_t0c, con_c, con_boltz, &
545 & con_plnk
546! --- output:
547 integer, intent(out) :: errflg
548 character(len=*), intent(out) :: errmsg
549
550! --- locals:
551 real (kind=kind_phys), dimension(NWVTOT) :: solfwv ! one wvn sol flux
552 real (kind=kind_phys), dimension(NWVTIR) :: eirfwv ! one wvn ir flux
553!
554!===> ... begin here
555!
556
557! Initialize CCPP error handling variables
558 errmsg = ''
559 errflg = 0
560
561 kyrstr = 1
562 kyrend = 1
563 kyrsav = 1
564 kmonsav = 1
565
566 laswflg= (mod(iaerflg,10) > 0) ! control flag for sw tropospheric aerosol
567 lalwflg= (mod(iaerflg/10,10) > 0) ! control flag for lw tropospheric aerosol
568 lavoflg= (mod(iaerflg/100,10) >0) ! control flag for stratospheric volcanic aeros
569
571
572 if ( me == 0 ) then
573
574 call wrt_aerlog(iaermdl, iaerflg, lalw1bd, errflg, errmsg) ! write aerosol param info to log file
575! --- inputs: (in scope variables)
576! --- outputs: (CCPP error handling)
577
578 endif
579
580 if ( iaerflg == 0 ) return ! return without any aerosol calculations
581
582! --- ... in sw, aerosols optical properties are computed for each radiation
583! spectral band; while in lw, optical properties can be calculated
584! for either only one broad band or for each of the lw radiation bands
585
586 if ( laswflg ) then
587 nswbnd = nbdsw
588 else
589 nswbnd = 0
590 endif
591
592 if ( lalwflg ) then
593 if ( lalw1bd ) then
594 nlwbnd = 1
595 else
596 nlwbnd = nbdlw
597 endif
598 else
599 nlwbnd = 0
600 endif
601
602 nswlwbd = nswbnd + nlwbnd
603
604 wvn_sw1(:) = wvnsw1(:)
605 wvn_sw2(:) = wvnsw2(:)
606 wvn_lw1(:) = wvnlw1(:)
607 wvn_lw2(:) = wvnlw2(:)
608
609! note: for result consistency, the defalt opac-clim aeros setting still use
610! old spectral band mapping. use iaermdl=5 to use new mapping method
611
612 if ( iaermdl == 0 ) then ! opac-climatology scheme
613 lmap_new = .false.
614
615 wvn_sw1(2:nbdsw-1) = wvn_sw1(2:nbdsw-1) + 1
616 wvn_lw1(2:nbdlw) = wvn_lw1(2:nbdlw) + 1
617 else
618 lmap_new = .true.
619 endif
620
621 if ( iaerflg /= 100 ) then
622
625
626 call set_spectrum(con_pi, con_t0c, con_c, con_boltz, con_plnk, &
627 & errflg, errmsg)
628! --- inputs: (module constants)
629! --- outputs: (ccpp error handling)
630
632
633 if ( iaermdl==0 .or. iaermdl==5 ) then ! opac-climatology scheme
634 call clim_aerinit &
635! --- inputs:
636 & ( solfwv, eirfwv, me, aeros_file, &
637! --- outputs:
638 & errflg, errmsg)
639
640 elseif ( iaermdl==1 .or. iaermdl==2 ) then ! gocart clim/prog scheme
641
642 call gocart_aerinit &
643! --- inputs:
644 & ( solfwv, eirfwv, me, &
645! --- outputs:
646 & errflg, errmsg)
647
648 else
649 if ( me == 0 ) then
650 print *,' !!! ERROR in aerosol model scheme selection', &
651 & ' iaermdl =',iaermdl
652 errflg = 1
653 errmsg = 'ERROR(aer_init): aerosol model scheme selected'// &
654 & 'is invalid'
655 return
656 endif
657 endif
658
659 endif ! end if_iaerflg_block
660
663
664 if ( lavoflg ) then
665
666 call set_volcaer(errflg, errmsg)
667! --- inputs: (module variables)
668! --- outputs: (module variables: ccpp error handling)
669
670 endif ! end if_lavoflg_block
671
672
673! =================
674 contains
675! =================
676
678!--------------------------------
679 subroutine wrt_aerlog(iaermdl, iaerflg, lalw1bd, errflg, errmsg)
680! ================================================================== !
681! !
682! subprogram : wrt_aerlog !
683! !
684! write aerosol parameter configuration to run log file. !
685! !
686! ==================== defination of variables =================== !
687! !
688! internal module variables: !
689! lalwflg - toposphere lw aerosol effect: =f:no; =t:yes !
690! laswflg - toposphere sw aerosol effect: =f:no; =t:yes !
691! lavoflg - stratosphere volcanic aeros effect: =f:no; =t:yes !
692! !
693! inputs: !
694! iaerflg - aerosol effect control flag: 3-digits (volc,lw,sw) !
695! iaermdl - tropospheric aerosol model scheme flag !
696! !
697! outputs: !
698! errmsg - CCPP error message !
699! errflg - CCPP error flag !
700! !
701! subroutines called: none !
702! !
703! usage: call wrt_aerlog !
704! !
705! ================================================================== !
706
707! --- inputs: ()
708 integer, intent(in) :: iaermdl, iaerflg
709 logical, intent(in) :: lalw1bd
710! --- output: (CCPP error handling)
711 integer, intent(out) :: errflg
712 character(len=*), intent(out) :: errmsg
713! --- locals:
714
715!
716!===> ... begin here
717!
718
719! Initialize CCPP error handling variables
720 errmsg = ''
721 errflg = 0
722
723 print *, vtagaer ! print out version tag
724
725 if ( iaermdl==0 .or. iaermdl==5 ) then
726 print *,' - Using OPAC-seasonal climatology for tropospheric', &
727 & ' aerosol effect'
728 elseif ( iaermdl == 1 ) then
729 print *,' - Using GOCART-climatology for tropospheric', &
730 & ' aerosol effect'
731 elseif ( iaermdl == 2 ) then
732 print *,' - Using GOCART-prognostic aerosols for tropospheric', &
733 & ' aerosol effect'
734 else
735 print *,' !!! ERROR in selection of aerosol model scheme', &
736 & ' IAER_MDL =',iaermdl
737 errflg = 1
738 errmsg = 'ERROR(wrt_aerlog): Selected aerosol model scheme is'//&
739 & 'is invalid'
740 return
741 endif ! end_if_iaermdl_block
742
743 print *,' IAER=',iaerflg,' LW-trop-aer=',lalwflg, &
744 & ' SW-trop-aer=',laswflg,' Volc-aer=',lavoflg
745
746 if ( iaerflg <= 0 ) then ! turn off all aerosol effects
747 print *,' - No tropospheric/volcanic aerosol effect included'
748 print *,' Input values of aerosol optical properties to' &
749 & ,' both SW and LW radiations are set to zeros'
750 else
751 if ( iaerflg >= 100 ) then ! incl stratospheric volcanic aerosols
752 print *,' - Include stratospheric volcanic aerosol effect'
753 else ! no stratospheric volcanic aerosols
754 print *,' - No stratospheric volcanic aerosol effect'
755 endif
756
757 if ( laswflg ) then ! chcek for sw effect
758 print *,' - Compute multi-band aerosol optical' &
759 & ,' properties for SW input parameters'
760 else
761 print *,' - No SW radiation aerosol effect, values of' &
762 & ,' aerosol properties to SW input are set to zeros'
763 endif
764
765 if ( lalwflg ) then ! check for lw effect
766 if ( lalw1bd ) then
767 print *,' - Compute 1 broad-band aerosol optical' &
768 & ,' properties for LW input parameters'
769 else
770 print *,' - Compute multi-band aerosol optical' &
771 & ,' properties for LW input parameters'
772 endif
773 else
774 print *,' - No LW radiation aerosol effect, values of' &
775 & ,' aerosol properties to LW input are set to zeros'
776 endif
777 endif ! end if_iaerflg_block
778!
779 return
780!................................
781 end subroutine wrt_aerlog
782!--------------------------------
783
787 subroutine set_spectrum(con_pi, con_t0c, con_c, con_boltz, &
788 & con_plnk, errflg, errmsg)
789
790! ================================================================== !
791! !
792! subprogram : set_spectrum !
793! !
794! define the one wavenumber solar fluxes based on toa solar spectral!
795! distrobution, and define the one wavenumber ir fluxes based on !
796! black-body emission distribution at a predefined temperature. !
797! !
798! ==================== defination of variables =================== !
799! !
820! !
821! subroutines called: none !
822! !
823! usage: call set_spectrum !
824! !
825! ================================================================== !
826
827! --- inputs: (module constants)
828! integer :: NWVTOT, NWVTIR
829! --- inputs: (CCPP Interstitials)
830 real(kind_phys),intent(in) :: con_pi, con_t0c, con_c, con_boltz, &
831 & con_plnk
832
833! --- output: (in-scope variables)
834! real (kind=kind_phys), dimension(NWVTOT) :: solfwv ! one wvn sol flux
835! real (kind=kind_phys), dimension(NWVTIR) :: eirfwv ! one wvn ir flux
836! --- output: (CCPP error-handling)
837 integer, intent(out) :: errflg
838 character(len=*), intent(out) :: errmsg
839! --- locals:
840 real (kind=kind_phys) :: soltot, tmp1, tmp2, tmp3
841
842 integer :: nb, nw, nw1, nw2, nmax, nmin
843
844! Initialize CCPP error handling variables
845 errmsg = ''
846 errflg = 0
847!
848!===> ... begin here
849!
850! nmax = min( NWVTOT, nint( maxval(wvnsw2) ))
851! nmin = max( 1, nint( minval(wvnsw1) ))
852
853! --- check print
854! print *,' MINWVN, MAXWVN = ',nmin, nmax
855! --- ... define the one wavenumber solar fluxes based on toa solar
856! spectral distribution
857
858! soltot1 = f_zero
859! soltot = f_zero
860 do nb = 1, nwvsol
861 if ( nb == 1 ) then
862 nw1 = 1
863 else
864 nw1 = nw1 + nwvns0(nb-1)
865 endif
866
867 nw2 = nw1 + nwvns0(nb) - 1
868
869 do nw = nw1, nw2
870 solfwv(nw) = s0intv(nb)
871! soltot1 = soltot1 + s0intv(nb)
872! if ( nw >= nmin .and. nw <= nmax ) then
873! soltot = soltot + s0intv(nb)
874! endif
875 enddo
876 enddo
877
878! --- ... define the one wavenumber ir fluxes based on black-body
879! emission distribution at a predefined temperature
880
881 tmp1 = (con_pi + con_pi) * con_plnk * con_c* con_c
882 tmp2 = con_plnk * con_c / (con_boltz * con_t0c)
883
884!$omp parallel do private(nw,tmp3)
885 do nw = 1, nwvtir
886 tmp3 = 100.0 * nw
887 eirfwv(nw) = (tmp1 * tmp3**3) / (exp(tmp2*tmp3) - 1.0)
888 enddo
889!
890 return
891!................................
892 end subroutine set_spectrum
893!--------------------------------
894
895
897!-----------------------------
898 subroutine set_volcaer(errflg, errmsg)
899!.............................
900! --- inputs: ( none ) !
901! outputs: (CCPP error handling) !
902! errflg - CCPP error flag !
903! errmsg - CCPP error message !
904! ================================================================== !
905! !
906! subprogram : set_volcaer !
907! !
908! this is the initialization progrmam for stratospheric volcanic !
909! aerosols. !
910! !
911! subroutines called: none !
912! !
913! usage: call set_volcaer !
914! !
915! ================================================================== !
916
917! --- inputs: (none)
918
919! --- output: (CCPP error handling)
920! integer :: ivolae(:,:,:)
921 integer, intent(out) :: errflg
922 character(len=*), intent(out) :: errmsg
923! --- locals:
924!
925!===> ... begin here
926!
927
928! Initialize CCPP error handling variables
929 errmsg = ''
930 errflg = 0
931
932! --- allocate data space
933
934 if ( .not. allocated(ivolae) ) then
935 allocate ( ivolae(12,4,10) ) ! for 12-mon,4-lat_zone,10-year
936 endif
937!
938 return
939!................................
940 end subroutine set_volcaer
941!--------------------------------
942!
943!...................................
944 end subroutine aer_init
945!-----------------------------------
946
947
957 subroutine clim_aerinit &
958 & ( solfwv, eirfwv, me, aeros_file, & ! --- inputs
959 & errflg, errmsg) ! --- outputs
960
961! ================================================================== !
962! !
963! clim_aerinit is the opac-climatology aerosol initialization program !
964! to set up necessary parameters and working arrays. !
965! !
966! inputs: !
967! solfwv(NWVTOT) - solar flux for each individual wavenumber (w/m2)!
968! eirfwv(NWVTIR) - ir flux(273k) for each individual wavenum (w/m2)!
969! me - print message control flag !
970! aeros_file - external aerosol data file name !
971! !
972! outputs: (CCPP error handling) !
973! errflg - CCPP error flag !
974! errmsg - CCPP error message !
975! !
976! internal module variables: !
977! lalwflg - logical lw aerosols effect control flag !
978! =t compute lw aerosol optical prop !
979! laswflg - logical sw aerosols effect control flag !
980! =t compute sw aerosol optical prop !
981! !
982! module constants: !
983! NWVSOL - num of wvnum regions where solar flux is constant !
984! NWVTOT - total num of wave numbers used in sw spectrum !
985! NWVTIR - total num of wave numbers used in the ir region !
986! NSWBND - total number of sw spectral bands !
987! NLWBND - total number of lw spectral bands !
988! NAERBND - number of bands for climatology aerosol data !
989! NCM1 - number of rh independent aeros species !
990! NCM2 - number of rh dependent aeros species !
991! !
992! usage: call clim_aerinit !
993! !
994! subprograms called: set_aercoef, optavg !
995! !
996! ================================================================== !
997
998! --- inputs:
999 real (kind=kind_phys), dimension(:) :: solfwv ! one wvn sol flux
1000 real (kind=kind_phys), dimension(:) :: eirfwv ! one wvn ir flux
1001 integer, intent(in) :: me
1002 character(len=26), intent(in) :: aeros_file
1003! --- output: (CCPP error handling)
1004 integer, intent(out) :: errflg
1005 character(len=*), intent(out) :: errmsg
1006
1007! --- locals:
1008 real (kind=kind_phys), dimension(NAERBND,NCM1) :: &
1009 & rhidext0, rhidsca0, rhidssa0, rhidasy0
1010 real (kind=kind_phys), dimension(NAERBND,NRHLEV,NCM2):: &
1011 & rhdpext0, rhdpsca0, rhdpssa0, rhdpasy0
1012 real (kind=kind_phys), dimension(NAERBND) :: straext0
1013
1014 real (kind=kind_phys), dimension(NSWBND,NAERBND) :: solwaer
1015 real (kind=kind_phys), dimension(NSWBND) :: solbnd
1016 real (kind=kind_phys), dimension(NLWBND,NAERBND) :: eirwaer
1017 real (kind=kind_phys), dimension(NLWBND) :: eirbnd
1018
1019 integer, dimension(NSWBND) :: nv1, nv2
1020 integer, dimension(NLWBND) :: nr1, nr2
1021!
1022!===> ... begin here
1023!
1024! Initialize CCPP error handling variables
1025 errmsg = ''
1026 errflg = 0
1027
1028! --- ... invoke tropospheric aerosol initialization
1029
1031 call set_aercoef(aeros_file, errflg, errmsg)
1032! --- inputs: (in-scope variables, module constants)
1033! --- outputs: (module variables)
1034
1035
1036! =================
1037 contains
1038! =================
1039
1044!--------------------------------
1045 subroutine set_aercoef(aeros_file,errflg, errmsg)
1046!................................
1047! --- inputs: (in-scope variables, module constants)
1048! --- outputs: (CCPP error handling)
1049
1050! ================================================================== !
1051! !
1052! subprogram : set_aercoef !
1053! !
1054! this is the initialization progrmam for climatological aerosols !
1055! !
1056! the program reads and maps the pre-tabulated aerosol optical !
1057! spectral data onto corresponding sw radiation spectral bands. !
1058! !
1059! ==================== defination of variables =================== !
1060! !
1061! inputs: (in-scope variables, module constants) !
1062! solfwv(:) - real, solar flux for individual wavenumber (w/m2) !
1063! eirfwv(:) - real, lw flux(273k) for individual wavenum (w/m2) !
1064! me - integer, select cpu number as print control flag !
1065! !
1066! outputs: (to the module variables) !
1067! outputs: (CCPP error handling) !
1068! errflg - CCPP error flag !
1069! errmsg - CCPP error message !
1070! !
1071! external module variables: !
1072! lalwflg - module control flag for lw trop-aer: =f:no; =t:yes !
1073! laswflg - module control flag for sw trop-aer: =f:no; =t:yes !
1074! aeros_file- external aerosol data file name !
1075! !
1076! internal module variables: !
1077! IMXAE - number of longitude points in global aeros data set !
1078! JMXAE - number of latitude points in global aeros data set !
1079! wvnsw1,wvnsw2 (NSWSTR:NSWEND) !
1080! - start/end wavenumbers for each of sw bands !
1081! wvnlw1,wvnlw2 ( 1:NBDLW) !
1082! - start/end wavenumbers for each of lw bands !
1083! NSWLWBD - total num of bands (sw+lw) for aeros optical properties!
1084! NSWBND - number of sw spectral bands actually invloved !
1085! NLWBND - number of lw spectral bands actually invloved !
1086! NIAERCM - unit number for reading input data set !
1087! extrhi - extinction coef for rh-indep aeros NCM1*NSWLWBD!
1088! scarhi - scattering coef for rh-indep aeros NCM1*NSWLWBD!
1089! ssarhi - single-scat-alb for rh-indep aeros NCM1*NSWLWBD!
1090! asyrhi - asymmetry factor for rh-indep aeros NCM1*NSWLWBD!
1091! extrhd - extinction coef for rh-dep aeros NRHLEV*NCM2*NSWLWBD!
1092! scarhd - scattering coef for rh-dep aeros NRHLEV*NCM2*NSWLWBD!
1093! ssarhd - single-scat-alb for rh-dep aeros NRHLEV*NCM2*NSWLWBD!
1094! asyrhd - asymmetry factor for rh-dep aeros NRHLEV*NCM2*NSWLWBD!
1095! !
1096! major local variables: !
1097! for handling spectral band structures !
1098! iendwv - ending wvnum (cm**-1) for each band NAERBND !
1099! for handling optical properties of rh independent species (NCM1) !
1100! 1. insoluble (inso); 2. soot (soot); !
1101! 3. mineral nuc mode (minm); 4. mineral acc mode (miam); !
1102! 5. mineral coa mode (micm); 6. mineral transport(mitr). !
1103! rhidext0 - extinction coefficient NAERBND*NCM1 !
1104! rhidsca0 - scattering coefficient NAERBND*NCM1 !
1105! rhidssa0 - single scattering albedo NAERBND*NCM1 !
1106! rhidasy0 - asymmetry parameter NAERBND*NCM1 !
1107! for handling optical properties of rh ndependent species (NCM2) !
1108! 1. water soluble (waso); 2. sea salt acc mode(ssam); !
1109! 3. sea salt coa mode(sscm); 4. sulfate droplets (suso). !
1110! rh level (NRHLEV): 00%, 50%, 70%, 80%, 90%, 95%, 98%, 99% !
1111! rhdpext0 - extinction coefficient NAERBND,NRHLEV,NCM2!
1112! rhdpsca0 - scattering coefficient NAERBND,NRHLEV,NCM2!
1113! rhdpssa0 - single scattering albedo NAERBND,NRHLEV,NCM2!
1114! rhdpasy0 - asymmetry parameter NAERBND,NRHLEV,NCM2!
1115! for handling optical properties of stratospheric bkgrnd aerosols !
1116! straext0 - extingction coefficients NAERBND !
1117! !
1118! usage: call set_aercoef !
1119! !
1120! subprograms called: optavg !
1121! !
1122! ================================================================== !
1123!
1124! --- inputs: ( none )
1125 character(len=26),intent(in) :: aeros_file
1126! --- output: (CCPP error handling)
1127 integer, intent(out) :: errflg
1128 character(len=*), intent(out) :: errmsg
1129
1130! --- locals:
1131 integer, dimension(NAERBND) :: iendwv
1132
1133 integer :: i, j, k, m, mb, ib, ii, id, iw, iw1, iw2, ik, ibs, ibe
1134
1135 real (kind=kind_phys) :: sumsol, sumir, fac, tmp, wvs, wve
1136
1137 logical :: file_exist
1138 character :: cline*80
1139!
1140!===> ... begin here
1141!
1142
1143! Initialize CCPP error handling variables
1144 errmsg = ''
1145 errflg = 0
1146
1149
1150 inquire (file=aeros_file, exist=file_exist)
1151
1152 if ( file_exist ) then
1153 close (niaercm)
1154 open (unit=niaercm,file=aeros_file,status='OLD', &
1155 & action='read',form='FORMATTED')
1156 rewind(niaercm)
1157 else
1158 print *,' Requested aerosol data file "',aeros_file, &
1159 & '" not found!'
1160 print *,' *** Stopped in subroutine aero_init !!'
1161 errflg = 1
1162 errmsg = 'ERROR(set_aercoef): Requested aerosol data file '// &
1163 & aeros_file//' not found'
1164 return
1165 endif ! end if_file_exist_block
1166
1167! --- ... skip monthly global distribution
1168
1169 do m = 1, 12
1170 read (niaercm,12) cline
1171 12 format(a80/)
1172
1173 do j = 1, jmxae
1174 do i = 1, imxae
1175 read(niaercm,*) id
1176 enddo
1177 enddo
1178 enddo ! end do_m_block
1179
1180! --- ... aloocate and input aerosol optical data
1181
1182 if ( .not. allocated( extrhi ) ) then
1183 allocate ( extrhi( ncm1,nswlwbd) )
1184 allocate ( scarhi( ncm1,nswlwbd) )
1185 allocate ( ssarhi( ncm1,nswlwbd) )
1186 allocate ( asyrhi( ncm1,nswlwbd) )
1187 allocate ( extrhd(nrhlev,ncm2,nswlwbd) )
1188 allocate ( scarhd(nrhlev,ncm2,nswlwbd) )
1189 allocate ( ssarhd(nrhlev,ncm2,nswlwbd) )
1190 allocate ( asyrhd(nrhlev,ncm2,nswlwbd) )
1191 allocate ( extstra( nswlwbd) )
1192 endif
1193
1195 read(niaercm,21) cline
1196 21 format(a80)
1197 read(niaercm,22) iendwv(:)
1198 22 format(13i6)
1199
1201 read(niaercm,21) cline
1202 read(niaercm,24) haer(:,:)
1203 24 format(20f4.1)
1204
1206 read(niaercm,21) cline
1207 read(niaercm,26) prsref(:,:)
1208 26 format(10f7.2)
1209
1211 read(niaercm,21) cline
1212 read(niaercm,28) rhidext0(:,:)
1213 28 format(8e10.3)
1214
1216 read(niaercm,21) cline
1217 read(niaercm,28) rhidsca0(:,:)
1218
1220 read(niaercm,21) cline
1221 read(niaercm,28) rhidssa0(:,:)
1222
1224 read(niaercm,21) cline
1225 read(niaercm,28) rhidasy0(:,:)
1226
1228 read(niaercm,21) cline
1229 read(niaercm,28) rhdpext0(:,:,:)
1230
1232 read(niaercm,21) cline
1233 read(niaercm,28) rhdpsca0(:,:,:)
1234
1236 read(niaercm,21) cline
1237 read(niaercm,28) rhdpssa0(:,:,:)
1238
1240 read(niaercm,21) cline
1241 read(niaercm,28) rhdpasy0(:,:,:)
1242
1244 read(niaercm,21) cline
1245 read(niaercm,28) straext0(:)
1246
1247 close (niaercm)
1248
1251
1252 sigref(:,:) = 0.001 * prsref(:,:)
1253
1256
1257 if ( laswflg ) then
1258 solbnd(:) = f_zero
1259!$omp parallel do private(i,j)
1260 do j=1,naerbnd
1261 do i=1,nswbnd
1262 solwaer(i,j) = f_zero
1263 enddo
1264 enddo
1265
1266 ibs = 1
1267 ibe = 1
1268 wvs = wvn_sw1(1)
1269 wve = wvn_sw1(1)
1270 nv_aod = 1
1271 do ib = 2, nswbnd
1272 mb = ib + nswstr - 1
1273 if ( wvn_sw2(mb) >= wvn550 .and. wvn550 >= wvn_sw1(mb) ) then
1274 nv_aod = ib ! sw band number covering 550nm wavelenth
1275 endif
1276
1277 if ( wvn_sw1(mb) < wvs ) then
1278 wvs = wvn_sw1(mb)
1279 ibs = ib
1280 endif
1281 if ( wvn_sw1(mb) > wve ) then
1282 wve = wvn_sw1(mb)
1283 ibe = ib
1284 endif
1285 enddo
1286
1287!$omp parallel do private(ib,mb,ii,iw1,iw2,iw,sumsol,fac,tmp,ibs,ibe)
1288 do ib = 1, nswbnd
1289 mb = ib + nswstr - 1
1290 ii = 1
1291 iw1 = nint(wvn_sw1(mb))
1292 iw2 = nint(wvn_sw2(mb))
1293
1294 lab_swdowhile : do while ( iw1 > iendwv(ii) )
1295 if ( ii == naerbnd ) exit lab_swdowhile
1296 ii = ii + 1
1297 enddo lab_swdowhile
1298
1299 if ( lmap_new ) then
1300 if (ib == ibs) then
1301 sumsol = f_zero
1302 else
1303 sumsol = -0.5 * solfwv(iw1)
1304 endif
1305 if (ib == ibe) then
1306 fac = f_zero
1307 else
1308 fac = -0.5
1309 endif
1310 solbnd(ib) = sumsol
1311 else
1312 sumsol = f_zero
1313 endif
1314 nv1(ib) = ii
1315
1316 do iw = iw1, iw2
1317 solbnd(ib) = solbnd(ib) + solfwv(iw)
1318 sumsol = sumsol + solfwv(iw)
1319
1320 if ( iw == iendwv(ii) ) then
1321 solwaer(ib,ii) = sumsol
1322
1323 if ( ii < naerbnd ) then
1324 sumsol = f_zero
1325 ii = ii + 1
1326 endif
1327 endif
1328 enddo
1329
1330 if ( iw2 /= iendwv(ii) ) then
1331 solwaer(ib,ii) = sumsol
1332 endif
1333
1334 if ( lmap_new ) then
1335 tmp = fac * solfwv(iw2)
1336 solwaer(ib,ii) = solwaer(ib,ii) + tmp
1337 solbnd(ib) = solbnd(ib) + tmp
1338 endif
1339
1340 nv2(ib) = ii
1341! frcbnd(ib) = solbnd(ib) / soltot
1342 enddo ! end do_ib_block for sw
1343 endif ! end if_laswflg_block
1344
1347
1348 if ( lalwflg ) then
1349 eirbnd(:) = f_zero
1350!$omp parallel do private(i,j)
1351 do j=1,naerbnd
1352 do i=1,nlwbnd
1353 eirwaer(i,j) = f_zero
1354 enddo
1355 enddo
1356
1357 ibs = 1
1358 ibe = 1
1359 if (nlwbnd > 1 ) then
1360 wvs = wvn_lw1(1)
1361 wve = wvn_lw1(1)
1362 do ib = 2, nlwbnd
1363 mb = ib + nlwstr - 1
1364 if ( wvn_lw1(mb) < wvs ) then
1365 wvs = wvn_lw1(mb)
1366 ibs = ib
1367 endif
1368 if ( wvn_lw1(mb) > wve ) then
1369 wve = wvn_lw1(mb)
1370 ibe = ib
1371 endif
1372 enddo
1373 endif
1374
1375!$omp parallel do private(ib,ii,iw1,iw2,iw,mb,sumir,fac,tmp,ibs,ibe)
1376 do ib = 1, nlwbnd
1377 ii = 1
1378 if ( nlwbnd == 1 ) then
1379! iw1 = 250 ! corresponding 40 mu
1380 iw1 = 400 ! corresponding 25 mu
1381 iw2 = 2500 ! corresponding 4 mu
1382 else
1383 mb = ib + nlwstr - 1
1384 iw1 = nint(wvn_lw1(mb))
1385 iw2 = nint(wvn_lw2(mb))
1386 endif
1387
1388 lab_lwdowhile : do while ( iw1 > iendwv(ii) )
1389 if ( ii == naerbnd ) exit lab_lwdowhile
1390 ii = ii + 1
1391 enddo lab_lwdowhile
1392
1393 if ( lmap_new ) then
1394 if (ib == ibs) then
1395 sumir = f_zero
1396 else
1397 sumir = -0.5 * eirfwv(iw1)
1398 endif
1399 if (ib == ibe) then
1400 fac = f_zero
1401 else
1402 fac = -0.5
1403 endif
1404 eirbnd(ib) = sumir
1405 else
1406 sumir = f_zero
1407 endif
1408 nr1(ib) = ii
1409
1410 do iw = iw1, iw2
1411 eirbnd(ib) = eirbnd(ib) + eirfwv(iw)
1412 sumir = sumir + eirfwv(iw)
1413
1414 if ( iw == iendwv(ii) ) then
1415 eirwaer(ib,ii) = sumir
1416
1417 if ( ii < naerbnd ) then
1418 sumir = f_zero
1419 ii = ii + 1
1420 endif
1421 endif
1422 enddo
1423
1424 if ( iw2 /= iendwv(ii) ) then
1425 eirwaer(ib,ii) = sumir
1426 endif
1427
1428 if ( lmap_new ) then
1429 tmp = fac * eirfwv(iw2)
1430 eirwaer(ib,ii) = eirwaer(ib,ii) + tmp
1431 eirbnd(ib) = eirbnd(ib) + tmp
1432 endif
1433
1434 nr2(ib) = ii
1435 enddo ! end do_ib_block for lw
1436 endif ! end if_lalwflg_block
1437
1440
1441 call optavg
1442! --- inputs: (in-scope variables, module variables)
1443! --- outputs: (module variables)
1444
1445! --- check print
1446! do ib = 1, NSWBND
1447! print *,' After optavg, for sw band:',ib
1448! print *,' extrhi:', extrhi(:,ib)
1449! print *,' scarhi:', scarhi(:,ib)
1450! print *,' ssarhi:', ssarhi(:,ib)
1451! print *,' asyrhi:', asyrhi(:,ib)
1452! mb = ib + NSWSTR - 1
1453! print *,' wvnsw1,wvnsw2 :',wvnsw1(mb),wvnsw2(mb)
1454! do i = 1, NRHLEV
1455! print *,' extrhd for rhlev:',i
1456! print *,extrhd(i,:,ib)
1457! print *,' scarhd for rhlev:',i
1458! print *,scarhd(i,:,ib)
1459! print *,' ssarhd for rhlev:',i
1460! print *,ssarhd(i,:,ib)
1461! print *,' asyrhd for rhlev:',i
1462! print *,asyrhd(i,:,ib)
1463! enddo
1464! print *,' extstra:', extstra(ib)
1465! enddo
1466! print *,' wvnlw1 :',wvnlw1
1467! print *,' wvnlw2 :',wvnlw2
1468! do ib = 1, NLWBND
1469! ii = NSWBND + ib
1470! print *,' After optavg, for lw band:',ib
1471! print *,' extrhi:', extrhi(:,ii)
1472! print *,' scarhi:', scarhi(:,ii)
1473! print *,' ssarhi:', ssarhi(:,ii)
1474! print *,' asyrhi:', asyrhi(:,ii)
1475! do i = 1, NRHLEV
1476! print *,' extrhd for rhlev:',i
1477! print *,extrhd(i,:,ii)
1478! print *,' scarhd for rhlev:',i
1479! print *,scarhd(i,:,ii)
1480! print *,' ssarhd for rhlev:',i
1481! print *,ssarhd(i,:,ii)
1482! print *,' asyrhd for rhlev:',i
1483! print *,asyrhd(i,:,ii)
1484! enddo
1485! print *,' extstra:', extstra(ii)
1486! enddo
1487!
1488 return
1489!................................
1490 end subroutine set_aercoef
1491!--------------------------------
1492
1497!--------------------------------
1498 subroutine optavg
1499!................................
1500! --- inputs: (in-scope variables, module variables
1501! --- outputs: (module variables)
1502
1503! ==================================================================== !
1504! !
1505! subprogram: optavg !
1506! !
1507! compute mean aerosols optical properties over each sw radiation !
1508! spectral band for each of the species components. This program !
1509! follows gfdl's approach for thick cloud opertical property in !
1510! sw radiation scheme (2000). !
1511! !
1512! ==================== defination of variables =================== !
1513! !
1514! major input variables: !
1515! nv1,nv2 (NSWBND) - start/end spectral band indices of aerosol data !
1516! for each sw radiation spectral band !
1517! nr1,nr2 (NLWBND) - start/end spectral band indices of aerosol data !
1518! for each ir radiation spectral band !
1519! solwaer (NSWBND,NAERBND) !
1520! - solar flux weight over each sw radiation band !
1521! vs each aerosol data spectral band !
1522! eirwaer (NLWBND,NAERBND) !
1523! - ir flux weight over each lw radiation band !
1524! vs each aerosol data spectral band !
1525! solbnd (NSWBND) - solar flux weight over each sw radiation band !
1526! eirbnd (NLWBND) - ir flux weight over each lw radiation band !
1527! NSWBND - total number of sw spectral bands !
1528! NLWBND - total number of lw spectral bands !
1529! !
1530! external module variables: !
1531! laswflg - control flag for sw spectral region !
1532! lalwflg - control flag for lw spectral region !
1533! !
1534! output variables: (to module variables) !
1535! !
1536! ================================================================== !
1537
1538! --- inputs:
1539! --- output:
1540
1541! --- locals:
1542 real (kind=kind_phys) :: sumk, sums, sumok, sumokg, sumreft, &
1543 & sp, refb, reft, rsolbd, rirbd
1544
1545 integer :: ib, nb, ni, nh, nc
1546!
1547!===> ... begin here
1548!
1549! --- ... loop for each sw radiation spectral band
1550
1551 if ( laswflg ) then
1552
1553!$omp parallel do private(nb,nc,sumk,sums,sumok,sumokg,sumreft)
1554!$omp+private(ni,nh,sp,reft,refb,rsolbd)
1555 do nb = 1, nswbnd
1556 rsolbd = f_one / solbnd(nb)
1557
1558! --- for rh independent aerosol species
1559
1560 do nc = 1, ncm1 ! --- for rh independent aerosol species
1561 sumk = f_zero
1562 sums = f_zero
1563 sumok = f_zero
1564 sumokg = f_zero
1565 sumreft = f_zero
1566
1567 do ni = nv1(nb), nv2(nb)
1568 sp = sqrt( (f_one - rhidssa0(ni,nc)) &
1569 & / (f_one - rhidssa0(ni,nc)*rhidasy0(ni,nc)) )
1570 reft = (f_one - sp) / (f_one + sp)
1571 sumreft = sumreft + reft*solwaer(nb,ni)
1572
1573 sumk = sumk + rhidext0(ni,nc)*solwaer(nb,ni)
1574 sums = sums + rhidsca0(ni,nc)*solwaer(nb,ni)
1575 sumok = sumok + rhidssa0(ni,nc)*solwaer(nb,ni) &
1576 & * rhidext0(ni,nc)
1577 sumokg = sumokg + rhidssa0(ni,nc)*solwaer(nb,ni) &
1578 & * rhidext0(ni,nc)*rhidasy0(ni,nc)
1579 enddo
1580
1581 refb = sumreft * rsolbd
1582
1583 extrhi(nc,nb) = sumk * rsolbd
1584 scarhi(nc,nb) = sums * rsolbd
1585 asyrhi(nc,nb) = sumokg / (sumok + 1.0e-10)
1586 ssarhi(nc,nb) = 4.0*refb &
1587 & / ( (f_one+refb)**2 - asyrhi(nc,nb)*(f_one-refb)**2 )
1588 enddo ! end do_nc_block for rh-ind aeros
1589
1590
1591 do nc = 1, ncm2 ! --- for rh dependent aerosols species
1592 do nh = 1, nrhlev
1593 sumk = f_zero
1594 sums = f_zero
1595 sumok = f_zero
1596 sumokg = f_zero
1597 sumreft = f_zero
1598
1599 do ni = nv1(nb), nv2(nb)
1600 sp = sqrt( (f_one - rhdpssa0(ni,nh,nc)) &
1601 & / (f_one - rhdpssa0(ni,nh,nc)*rhdpasy0(ni,nh,nc)) )
1602 reft = (f_one - sp) / (f_one + sp)
1603 sumreft = sumreft + reft*solwaer(nb,ni)
1604
1605 sumk = sumk + rhdpext0(ni,nh,nc)*solwaer(nb,ni)
1606 sums = sums + rhdpsca0(ni,nh,nc)*solwaer(nb,ni)
1607 sumok = sumok + rhdpssa0(ni,nh,nc)*solwaer(nb,ni) &
1608 & * rhdpext0(ni,nh,nc)
1609 sumokg = sumokg + rhdpssa0(ni,nh,nc)*solwaer(nb,ni) &
1610 & * rhdpext0(ni,nh,nc)*rhdpasy0(ni,nh,nc)
1611 enddo
1612
1613 refb = sumreft * rsolbd
1614
1615 extrhd(nh,nc,nb) = sumk * rsolbd
1616 scarhd(nh,nc,nb) = sums * rsolbd
1617 asyrhd(nh,nc,nb) = sumokg / (sumok + 1.0e-10)
1618 ssarhd(nh,nc,nb) = 4.0*refb &
1619 & / ( (f_one+refb)**2 - asyrhd(nh,nc,nb)*(f_one-refb)**2 )
1620 enddo ! end do_nh_block
1621 enddo ! end do_nc_block for rh-dep aeros
1622
1623! --- for stratospheric background aerosols
1624
1625 sumk = f_zero
1626 do ni = nv1(nb), nv2(nb)
1627 sumk = sumk + straext0(ni)*solwaer(nb,ni)
1628 enddo
1629
1630 extstra(nb) = sumk * rsolbd
1631
1632! --- check print
1633! if ( nb > 6 .and. nb < 10) then
1634! print *,' in optavg for sw band',nb
1635! print *,' nv1, nv2:',nv1(nb),nv2(nb)
1636! print *,' solwaer:',solwaer(nb,nv1(nb):nv2(nb))
1637! print *,' extrhi:', extrhi(:,nb)
1638! do i = 1, NRHLEV
1639! print *,' extrhd for rhlev:',i
1640! print *,extrhd(i,:,nb)
1641! enddo
1642! print *,' sumk, rsolbd, extstra:',sumk,rsolbd,extstra(nb)
1643! endif
1644
1645 enddo ! end do_nb_block for sw
1646 endif ! end if_laswflg_block
1647
1648! --- ... loop for each lw radiation spectral band
1649
1650 if ( lalwflg ) then
1651
1652!$omp parallel do private(nb,ib,nc,rirbd,sumk,sums,sumok,sumokg,sumreft)
1653!$omp+private(ni,nh,sp,reft,refb,rsolbd)
1654 do nb = 1, nlwbnd
1655
1656 ib = nswbnd + nb
1657 rirbd = f_one / eirbnd(nb)
1658
1659 do nc = 1, ncm1 ! --- for rh independent aerosol species
1660 sumk = f_zero
1661 sums = f_zero
1662 sumok = f_zero
1663 sumokg = f_zero
1664 sumreft = f_zero
1665
1666 do ni = nr1(nb), nr2(nb)
1667 sp = sqrt( (f_one - rhidssa0(ni,nc)) &
1668 & / (f_one - rhidssa0(ni,nc)*rhidasy0(ni,nc)) )
1669 reft = (f_one - sp) / (f_one + sp)
1670 sumreft = sumreft + reft*eirwaer(nb,ni)
1671
1672 sumk = sumk + rhidext0(ni,nc)*eirwaer(nb,ni)
1673 sums = sums + rhidsca0(ni,nc)*eirwaer(nb,ni)
1674 sumok = sumok + rhidssa0(ni,nc)*eirwaer(nb,ni) &
1675 & * rhidext0(ni,nc)
1676 sumokg = sumokg + rhidssa0(ni,nc)*eirwaer(nb,ni) &
1677 & * rhidext0(ni,nc)*rhidasy0(ni,nc)
1678 enddo
1679
1680 refb = sumreft * rirbd
1681
1682 extrhi(nc,ib) = sumk * rirbd
1683 scarhi(nc,ib) = sums * rirbd
1684 asyrhi(nc,ib) = sumokg / (sumok + 1.0e-10)
1685 ssarhi(nc,ib) = 4.0*refb &
1686 & / ( (f_one+refb)**2 - asyrhi(nc,ib)*(f_one-refb)**2 )
1687 enddo ! end do_nc_block for rh-ind aeros
1688
1689 do nc = 1, ncm2 ! --- for rh dependent aerosols species
1690 do nh = 1, nrhlev
1691 sumk = f_zero
1692 sums = f_zero
1693 sumok = f_zero
1694 sumokg = f_zero
1695 sumreft = f_zero
1696
1697 do ni = nr1(nb), nr2(nb)
1698 sp = sqrt( (f_one - rhdpssa0(ni,nh,nc)) &
1699 & / (f_one - rhdpssa0(ni,nh,nc)*rhdpasy0(ni,nh,nc)) )
1700 reft = (f_one - sp) / (f_one + sp)
1701 sumreft = sumreft + reft*eirwaer(nb,ni)
1702
1703 sumk = sumk + rhdpext0(ni,nh,nc)*eirwaer(nb,ni)
1704 sums = sums + rhdpsca0(ni,nh,nc)*eirwaer(nb,ni)
1705 sumok = sumok + rhdpssa0(ni,nh,nc)*eirwaer(nb,ni) &
1706 & * rhdpext0(ni,nh,nc)
1707 sumokg = sumokg + rhdpssa0(ni,nh,nc)*eirwaer(nb,ni) &
1708 & * rhdpext0(ni,nh,nc)*rhdpasy0(ni,nh,nc)
1709 enddo
1710
1711 refb = sumreft * rirbd
1712
1713 extrhd(nh,nc,ib) = sumk * rirbd
1714 scarhd(nh,nc,ib) = sums * rirbd
1715 asyrhd(nh,nc,ib) = sumokg / (sumok + 1.0e-10)
1716 ssarhd(nh,nc,ib) = 4.0*refb &
1717 & / ( (f_one+refb)**2 - asyrhd(nh,nc,ib)*(f_one-refb)**2 )
1718 enddo ! end do_nh_block
1719 enddo ! end do_nc_block for rh-dep aeros
1720
1721! --- for stratospheric background aerosols
1722
1723 sumk = f_zero
1724 do ni = nr1(nb), nr2(nb)
1725 sumk = sumk + straext0(ni)*eirwaer(nb,ni)
1726 enddo
1727
1728 extstra(ib) = sumk * rirbd
1729
1730! --- check print
1731! if ( nb >= 1 .and. nb < 5) then
1732! print *,' in optavg for ir band:',nb
1733! print *,' nr1, nr2:',nr1(nb),nr2(nb)
1734! print *,' eirwaer:',eirwaer(nb,nr1(nb):nr2(nb))
1735! print *,' extrhi:', extrhi(:,ib)
1736! do i = 1, NRHLEV
1737! print *,' extrhd for rhlev:',i
1738! print *,extrhd(i,:,ib)
1739! enddo
1740! print *,' sumk, rirbd, extstra:',sumk,rirbd,extstra(ib)
1741! endif
1742
1743 enddo ! end do_nb_block for lw
1744 endif ! end if_lalwflg_block
1745!
1746 return
1747!................................
1748 end subroutine optavg
1749!--------------------------------
1750!
1751!...................................
1752 end subroutine clim_aerinit
1753!-----------------------------------
1754
1755
1763!-----------------------------------
1764 subroutine aer_update &
1765 & ( iyear, imon, me, iaermdl, aeros_file, errflg, errmsg )
1766
1767! ================================================================== !
1768! !
1769! aer_update checks and update time varying climatology aerosol !
1770! data sets. !
1771! !
1772! inputs: size !
1773! iyear - 4-digit calender year 1 !
1774! imon - month of the year 1 !
1775! me - print message control flag 1 !
1776! iaermdl - tropospheric aerosol model scheme flag 1 !
1777! aeros_file - external aerosol data file name len=26 !
1778! !
1779! outputs: (CCPP error handling) len=* !
1780! errmsg - CCPP error message 1 !
1781! errflg - CCPP error flag !
1782! !
1783! internal module variables: !
1784! lalwflg - control flag for tropospheric lw aerosol !
1785! laswflg - control flag for tropospheric sw aerosol !
1786! lavoflg - control flag for stratospheric volcanic aerosol !
1787! !
1788! usage: call aero_update !
1789! !
1790! subprograms called: trop_update, volc_update !
1791! !
1792! ================================================================== !
1793
1794! --- inputs:
1795 integer, intent(in) :: iyear, imon, me, iaermdl
1796 character(len=26),intent(in) :: aeros_file
1797! --- output: (CCPP error-handling)
1798 integer, intent(out) :: errflg
1799 character(len=*), intent(out) :: errmsg
1800! --- locals: ( none )
1801!
1802!===> ... begin here
1803!
1804
1805! Initialize CCPP error handling variables
1806 errmsg = ''
1807 errflg = 0
1808
1809 if ( imon < 1 .or. imon > 12 ) then
1810 print *,' ***** ERROR in specifying requested month !!! ', &
1811 & 'imon=', imon
1812 print *,' ***** STOPPED in subroutinte aer_update !!!'
1813 errflg = 1
1814 errmsg = 'ERROR(aer_update): Requested month not valid'
1815 return
1816 endif
1817
1819 if ( lalwflg .or. laswflg ) then
1820
1821 if ( iaermdl == 0 .or. iaermdl==5 ) then ! opac-climatology scheme
1822 call trop_update(aeros_file, errflg, errmsg)
1823 endif
1824
1825 endif
1826
1828 if ( lavoflg ) then
1829 call volc_update(errflg, errmsg)
1830 endif
1831
1832
1833! =================
1834 contains
1835! =================
1836
1839!--------------------------------
1840 subroutine trop_update(aeros_file, errflg, errmsg)
1841
1842! ================================================================== !
1843! !
1844! subprogram : trop_update !
1845! !
1846! updates the monthly global distribution of aerosol profiles in !
1847! five degree horizontal resolution. !
1848! !
1849! ==================== defination of variables =================== !
1850! !
1851! inputs: (in-scope variables, module constants) !
1852! imon - integer, month of the year !
1853! me - integer, print message control flag !
1854! inputs: (CCPP Interstitials) !
1855! aeros_file - external aerosol data file name !
1856! !
1857! outputs: (module variables) !
1858!
1859! outputs: (CCPP error-handling) !
1860! errmsg - Error message !
1861! errflg - Error flag !
1862! !
1863! internal module variables: !
1864! kprfg ( IMXAE*JMXAE) - aeros profile index !
1865! idxcg (NXC*IMXAE*JMXAE) - aeros component index !
1866! cmixg (NXC*IMXAE*JMXAE) - aeros component mixing ratio !
1867! denng ( 2 *IMXAE*JMXAE) - aerosols number density !
1868! !
1869! NIAERCM - unit number for input data set !
1870! !
1871! subroutines called: none !
1872! !
1873! usage: call trop_update !
1874! !
1875! ================================================================== !
1876
1877! --- inputs: (CCPP Interstitials)
1878 character(len=26),intent(in) :: aeros_file
1879! --- output: (CCPP error handling)
1880 integer, intent(out) :: errflg
1881 character(len=*), intent(out) :: errmsg
1882
1883! --- locals:
1884! real (kind=kind_io8) :: cmix(NXC), denn, tem
1885 real (kind=kind_phys) :: cmix(nxc), denn, tem
1886 integer :: idxc(nxc), kprf
1887
1888 integer :: i, id, j, k, m, nc
1889 logical :: file_exist
1890
1891 character :: cline*80, ctyp*3
1892!
1893!===> ... begin here
1894!
1895
1896! Initialize CCPP error handling variables
1897 errmsg = ''
1898 errflg = 0
1899
1900! --- ... reading climatological aerosols data
1901
1902 inquire (file=aeros_file, exist=file_exist)
1903
1904 if ( file_exist ) then
1905 close(niaercm)
1906 open (unit=niaercm,file=aeros_file,status='OLD', &
1907 & action='read',form='FORMATTED')
1908 rewind(niaercm)
1909
1910 if ( me == 0 ) then
1911 print *,' Opened aerosol data file: ',aeros_file
1912 endif
1913 else
1914 print *,' Requested aerosol data file "',aeros_file, &
1915 & '" not found!'
1916 print *,' *** Stopped in subroutine trop_update !!'
1917 errflg = 1
1918 errmsg = 'ERROR(trop_update):Requested aerosol data file '// &
1919 & aeros_file // ' not found.'
1920 return
1921 endif ! end if_file_exist_block
1922
1923!$omp parallel do private(i,j,m)
1924 do j = 1, jmxae
1925 do i = 1, imxae
1926 do m = 1, nxc
1927 idxcg(m,i,j) = 0
1928 cmixg(m,i,j) = f_zero
1929 enddo
1930 enddo
1931 enddo
1932
1933!$omp parallel do private(i,j)
1934 do j = 1, jmxae
1935 do i = 1, imxae
1936 denng(1,i,j) = f_zero
1937 denng(2,i,j) = f_zero
1938 enddo
1939 enddo
1940
1941! --- ... loop over 12 month global distribution
1942
1943 lab_do_12mon : do m = 1, 12
1944
1945 read(niaercm,12) cline
1946 12 format(a80/)
1947
1948 if ( m /= imon ) then
1949! if ( me == 0 ) print *,' *** Skipped ',cline
1950
1951 do j = 1, jmxae
1952 do i = 1, imxae
1953 read(niaercm,*) id
1954 enddo
1955 enddo
1956 else
1957 if ( me == 0 ) print *,' --- Reading ',cline
1958
1959 do j = 1, jmxae
1960 do i = 1, imxae
1961 read(niaercm,14) (idxc(k),cmix(k),k=1,nxc),kprf,denn,nc, &
1962 & ctyp
1963 14 format(5(i2,e11.4),i2,f8.2,i3,1x,a3)
1964
1965 kprfg(i,j) = kprf
1966 denng(1,i,j) = denn ! num density of 1st layer
1967 if ( kprf >= 6 ) then
1968 denng(2,i,j) = cmix(nxc) ! num density of 2dn layer
1969 else
1970 denng(2,i,j) = f_zero
1971 endif
1972
1973 tem = f_one
1974 do k = 1, nxc-1
1975 idxcg(k,i,j) = idxc(k) ! component index
1976 cmixg(k,i,j) = cmix(k) ! component mixing ratio
1977 tem = tem - cmix(k)
1978 enddo
1979 idxcg(nxc,i,j) = idxc(nxc)
1980 cmixg(nxc,i,j) = tem ! to make sure all add to 1.
1981 enddo
1982 enddo
1983
1984 close (niaercm)
1985 exit lab_do_12mon
1986 endif ! end if_m_block
1987
1988 enddo lab_do_12mon
1989
1990! -- check print
1991
1992! print *,' IDXCG :'
1993! print 16,idxcg
1994! 16 format(40i3)
1995! print *,' CMIXG :'
1996! print 17,cmixg
1997! print *,' DENNG :'
1998! print 17,denng
1999! print *,' KPRFG :'
2000! print 17,kprfg
2001! 17 format(8e16.9)
2002!
2003 return
2004!................................
2005 end subroutine trop_update
2006!--------------------------------
2007
2008
2011!--------------------------------
2012 subroutine volc_update(errflg, errmsg)
2013!................................
2014! --- inputs: (in scope variables, module variables)
2015! --- outputs: (CCPP error handling)
2016
2017! ================================================================== !
2018! !
2019! subprogram : volc_update !
2020! !
2021! searches historical volcanic data sets to find and read in !
2022! monthly 45-degree lat-zone band data of optical depth. !
2023! !
2024! ==================== defination of variables =================== !
2025! !
2026! inputs: (in-scope variables, module constants) !
2027! iyear - integer, 4-digit calender year 1 !
2028! imon - integer, month of the year 1 !
2029! me - integer, print message control flag 1 !
2030! NIAERCM - integer, unit number for input data set 1 !
2031! !
2032! outputs: (module variables) !
2033! ivolae - integer, monthly, 45-deg lat-zone volc odp 12*4*10 !
2034! kyrstr - integer, starting year of data in the input file !
2035! kyrend - integer, ending year of data in the input file !
2036! kyrsav - integer, the year of data in use in the input file !
2037! kmonsav - integer, the month of data in use in the input file !
2038! !
2039! outputs: (CCPP error-handling) !
2040! errmsg - Error message !
2041! errflg - Error flag !
2042! !
2043! subroutines called: none !
2044! !
2045! usage: call volc_aerinit !
2046! !
2047! ================================================================== !
2048
2049! --- inputs: (in-scope variables, module constants)
2050! integer :: iyear, imon, me, NIAERCM
2051
2052! --- output: (module variables)
2053! integer :: ivolae(:,:,:), kyrstr, kyrend, kyrsav, kmonsav
2054! --- output: (CCPP error-handling)
2055 integer, intent(out) :: errflg
2056 character(len=*), intent(out) :: errmsg
2057
2058! --- locals:
2059 integer :: i, j, k
2060 logical :: file_exist
2061
2062 character :: cline*80, volcano_file*32
2063 data volcano_file / 'volcanic_aerosols_1850-1859.txt ' /
2064!
2065!===> ... begin here
2066!
2067
2068! Initialize CCPP error handling variables
2069 errmsg = ''
2070 errflg = 0
2071
2072 kmonsav = imon
2073
2074 if ( kyrstr<=iyear .and. iyear<=kyrend ) then ! use previously input data
2075 kyrsav = iyear
2076 return
2077 else ! need to input new data
2078 kyrsav = iyear
2079 kyrstr = iyear - mod(iyear,10)
2080 kyrend = kyrstr + 9
2081
2082! --- check print
2083! print *,' kyrstr, kyrend, kyrsav, kmonsav =', &
2084! & kyrstr,kyrend,kyrsav,kmonsav
2085
2086 if ( iyear < minvyr .or. iyear > maxvyr ) then
2087! if ( .not. allocated(ivolae) ) then
2088! allocate ( ivolae(12,4,10) ) ! for 12-mon,4-lat_zone,10-year
2089! endif
2090 ivolae(:,:,:) = 1 ! set as lowest value
2091 if ( me == 0 ) then
2092 print *,' Request volcanic date out of range,', &
2093 & ' optical depth set to lowest value'
2094 endif
2095 else
2096 write(volcano_file(19:27),60) kyrstr,kyrend
2097 60 format(i4.4,'-',i4.4)
2098
2099 inquire (file=volcano_file, exist=file_exist)
2100 if ( file_exist ) then
2101 close(niaercm)
2102 open (unit=niaercm,file=volcano_file,status='OLD', &
2103 & action='read',form='FORMATTED')
2104
2105 read(niaercm,62) cline
2106 62 format(a80)
2107
2108! --- check print
2109 if ( me == 0 ) then
2110 print *,' Opened volcanic data file: ',volcano_file
2111 print *, cline
2112 endif
2113
2114 do k = 1, 10
2115 do j = 1, 4
2116 read(niaercm,64) (ivolae(i,j,k),i=1,12)
2117 64 format(12i5)
2118 enddo
2119 enddo
2120
2121 close (niaercm)
2122 else
2123 print *,' Requested volcanic data file "', &
2124 & volcano_file,'" not found!'
2125 print *,' *** Stopped in subroutine VOLC_AERINIT !!'
2126 errflg = 1
2127 errmsg = 'ERROR(volc_update): Requested volcanic data '// &
2128 & 'file '//volcano_file//' not found!'
2129 return
2130 endif ! end if_file_exist_block
2131
2132 endif ! end if_iyear_block
2133 endif ! end if_kyrstr_block
2134
2135! --- check print
2136 if ( me == 0 ) then
2137 k = mod(kyrsav,10) + 1
2138 print *,' CHECK: Sample Volcanic data used for month, year:', &
2139 & imon, iyear
2140 print *, ivolae(kmonsav,:,k)
2141 endif
2142!
2143 return
2144!................................
2145 end subroutine volc_update
2146!--------------------------------
2147!
2148!...................................
2149 end subroutine aer_update
2150!-----------------------------------
2151
2177!-----------------------------------
2178 subroutine setaer &
2179 & ( prsi,prsl,prslk,tvly,rhlay,slmsk,tracer,aerfld,xlon,xlat, & ! --- inputs
2180 & imax,nlay,nlp1, lsswr,lslwr,iaermdl,iaerflg,top_at_1, &
2181 & con_pi,con_rd,con_g,aerosw,aerolw, & ! --- outputs
2182 & aerodp, ext550, errflg, errmsg &
2183 & )
2184
2185! ================================================================== !
2186! !
2187! setaer computes aerosols optical properties !
2188! !
2189! inputs: size !
2190! prsi - pressure at interface mb IMAX*NLP1 !
2191! prsl - layer mean pressure mb IMAX*NLAY !
2192! prslk - exner function = (p/p0)**rocp IMAX*NLAY !
2193! tvly - layer virtual temperature k IMAX*NLAY !
2194! rhlay - layer mean relative humidity IMAX*NLAY !
2195! slmsk - sea/land mask (sea:0,land:1,sea-ice:2) IMAX !
2196! tracer - aerosol tracer concentration IMAX*NLAY*NTRAC !
2197! aerfld - prescribed aerosol mixing rat IMAX*NLAY*NTRCAER!
2198! xlon - longitude of given points in radiance IMAX !
2199! ok for both 0->2pi or -pi->+pi ranges !
2200! xlat - latitude of given points in radiance IMAX !
2201! default to pi/2 -> -pi/2, otherwise see in-line comment!
2202! IMAX - horizontal dimension of arrays 1 !
2203! NLAY,NLP1-vertical dimensions of arrays 1 !
2204! lsswr,lslwr !
2205! - logical flags for sw/lw radiation calls 1 !
2206! con_pi - Physical constant (pi) !
2207! con_t0c - Physical constant (temperature kelvin at zero celcius) !
2208! con_c - Physical constant (speed of light) !
2209! iaermdl - tropospheric aerosol model scheme flag !
2210! iaerflg - aerosol effect control flag !
2211! top_at_1 - Vertical ordering convection flag !
2212! !
2213! outputs: !
2214! aerosw - aeros opt properties for sw IMAX*NLAY*NBDSW*NF_AESW!
2215! (:,:,:,1): optical depth !
2216! (:,:,:,2): single scattering albedo !
2217! (:,:,:,3): asymmetry parameter !
2218! aerolw - aeros opt properties for lw IMAX*NLAY*NBDLW*NF_AELW!
2219! (:,:,:,1): optical depth !
2220! (:,:,:,2): single scattering albedo !
2221! (:,:,:,3): asymmetry parameter !
2222! tau_gocart - 550nm aeros opt depth IMAX*NLAY*MAX_NUM_GRIDCOMP!
2223!! aerodp - vertically integrated optical depth IMAX*NSPC1 !
2224! !
2225! errflg - CCPP error flag !
2226! errmsg - CCPP error message !
2227! !
2228! internal module variable: !
2229! laswflg - tropospheric aerosol control flag for sw radiation !
2230! =f: no sw aeros calc. =t: do sw aeros calc. !
2231! lalwflg - tropospheric aerosol control flag for lw radiation !
2232! =f: no lw aeros calc. =t: do lw aeros calc. !
2233! lavoflg - control flag for stratospheric vocanic aerosols !
2234! =t: add volcanic aerosols to the background aerosols !
2235! internal module variable: (set by subroutine aer_init) !
2236! ivolae - stratosphere volcanic aerosol optical depth (fac 1.e4) !
2237! 12*4*10 !
2238! usage: call setaer !
2239! !
2240! subprograms called: aer_property !
2241! !
2242! ================================================================== !
2243
2244! --- inputs:
2245 integer, intent(in) :: imax, nlay, nlp1, iaermdl, iaerflg
2246 real (kind=kind_phys), intent(in) :: con_pi, con_rd, con_g
2247 real (kind=kind_phys), dimension(:,:), intent(in) :: prsi, prsl, &
2248 & prslk, tvly, rhlay
2249 real (kind=kind_phys), dimension(:), intent(in) :: xlon, xlat, &
2250 & slmsk
2251 real (kind=kind_phys), dimension(:,:,:),intent(in):: tracer
2252 real (kind=kind_phys), dimension(:,:,:),intent(in):: aerfld
2253
2254 logical, intent(in) :: lsswr, lslwr, top_at_1
2255
2256
2257! --- outputs:
2258 real (kind=kind_phys), dimension(:,:,:,:), intent(out) :: &
2259 & aerosw, aerolw
2260
2261 real (kind=kind_phys), dimension(:,:) , intent(out) :: aerodp
2262 real (kind=kind_phys), dimension(:,:) , intent(out) :: ext550
2263 integer, intent(out) :: errflg
2264 character(len=*), intent(out) :: errmsg
2265
2266! --- locals:
2267 real (kind=kind_phys), parameter :: psrfh = 5.0 ! ref press (mb) for upper bound
2268
2269 real (kind=kind_phys), dimension(IMAX) :: alon,alat,volcae,rdelp
2270! real (kind=kind_phys), dimension(IMAX) :: sumodp
2271 real (kind=kind_phys) :: prsln(nlp1),hz(imax,nlp1),dz(imax,nlay)
2272 real (kind=kind_phys) :: tmp1, tmp2, psrfl
2273
2274 integer :: kcutl(imax), kcuth(imax)
2275 integer :: i, i1, j, k, m, mb, kh, kl
2276
2277 logical :: laddsw=.false., laersw=.false.
2278 logical :: laddlw=.false., laerlw=.false.
2279
2280! --- conversion constants
2281 real (kind=kind_phys) :: rdg
2282 real (kind=kind_phys) :: rovg
2283
2284!===> ... begin here
2285 rdg = 180._kind_phys / con_pi
2286 rovg = 0.001_kind_phys * con_rd / con_g
2287
2288! Initialize CCPP error handling variables
2289 errmsg = ''
2290 errflg = 0
2291
2292 do m = 1, nf_aesw
2293 do j = 1, nbdsw
2294 do k = 1, nlay
2295 do i = 1, imax
2296 aerosw(i,k,j,m) = f_zero
2297 enddo
2298 enddo
2299 enddo
2300 enddo
2301
2302 do m = 1, nf_aelw
2303 do j = 1, nbdlw
2304 do k = 1, nlay
2305 do i = 1, imax
2306 aerolw(i,k,j,m) = f_zero
2307 enddo
2308 enddo
2309 enddo
2310 enddo
2311
2312! sumodp = f_zero
2313 do i = 1, imax
2314 do k = 1, nspc1
2315 aerodp(i,k) = f_zero
2316 enddo
2317 enddo
2318 ext550(:,:) = f_zero
2319
2320 if ( .not. (lsswr .or. lslwr) ) then
2321 return
2322 endif
2323
2324 if ( iaerflg <= 0 ) then
2325 return
2326 endif
2327
2328 laersw = lsswr .and. laswflg
2329 laerlw = lslwr .and. lalwflg
2330
2332
2333 do i = 1, imax
2334 alon(i) = xlon(i) * rdg
2335 if (alon(i) < f_zero) alon(i) = alon(i) + 360.0
2336 alat(i) = xlat(i) * rdg ! if xlat in pi/2 -> -pi/2 range
2337! alat(i) = 90.0 - xlat(i)*rdg ! if xlat in 0 -> pi range
2338 enddo
2339
2341
2342 if ( laswflg .or. lalwflg ) then
2343
2344 lab_do_imax : do i = 1, imax
2345
2346 lab_if_flip : if (.not. top_at_1) then ! input from sfc to toa
2347
2348 do k = 1, nlay
2349 prsln(k) = log(prsi(i,k))
2350 enddo
2351 prsln(nlp1)= log(prsl(i,nlay))
2352
2353 do k = nlay, 1, -1
2354 dz(i,k) = rovg * (prsln(k) - prsln(k+1)) * tvly(i,k)
2355 enddo
2356 dz(i,nlay) = 2.0 * dz(i,nlay)
2357
2358 hz(i,1) = f_zero
2359 do k = 1, nlay
2360 hz(i,k+1) = hz(i,k) + dz(i,k)
2361 enddo
2362
2363 else lab_if_flip ! input from toa to sfc
2364
2365 prsln(1) = log(prsl(i,1))
2366 do k = 2, nlp1
2367 prsln(k) = log(prsi(i,k))
2368 enddo
2369
2370 do k = 1, nlay
2371 dz(i,k) = rovg * (prsln(k+1) - prsln(k)) * tvly(i,k)
2372 enddo
2373 dz(i,1) = 2.0 * dz(i,1)
2374
2375 hz(i,nlp1) = f_zero
2376 do k = nlay, 1, -1
2377 hz(i,k) = hz(i,k+1) + dz(i,k)
2378 enddo
2379
2380 endif lab_if_flip
2381
2382 enddo lab_do_imax
2383
2384
2394
2395 if ( iaermdl==0 .or. iaermdl==5 ) then ! use opac aerosol climatology
2396
2397 call aer_property &
2398! --- inputs:
2399 & ( prsi,prsl,prslk,tvly,rhlay,dz,hz,tracer, &
2400 & alon,alat,slmsk, laersw,laerlw, &
2401 & imax,nlay,nlp1,top_at_1, &
2402! & IMAX,NLAY,NLP1,NSPC1, &
2403! --- outputs:
2404 & aerosw,aerolw,aerodp,errflg,errmsg &
2405 & )
2406
2407!
2408 elseif ( iaermdl==1 .or. iaermdl==2) then ! use gocart aerosols
2409
2410 call aer_property_gocart &
2411! --- inputs:
2412 & ( prsi,prsl,prslk,tvly,rhlay,dz,hz,tracer,aerfld, &
2413 & alon,alat,slmsk,laersw,laerlw,con_rd, &
2414 & imax,nlay,nlp1, &
2415! --- outputs:
2416 & aerosw,aerolw,aerodp,ext550,errflg,errmsg &
2417 & )
2418 endif ! end if_iaerflg_block
2419
2420
2421! --- check print
2422! do m = 1, NBDSW
2423! print *,' *** CHECK AEROSOLS PROPERTIES FOR SW BAND =',m, &
2424! & ' ***'
2425! do k = 1, 10
2426! print *,' LEVEL :',k
2427! print *,' TAUAER:',aerosw(:,k,m,1)
2428! print *,' SSAAER:',aerosw(:,k,m,2)
2429! print *,' ASYAER:',aerosw(:,k,m,3)
2430! enddo
2431! enddo
2432! print *,' *** CHECK AEROSOLS OPTICAL DEPTH FOR 550nm REGION'
2433! print *, aerodp(:,1)
2434! if ( laod_out ) then
2435! do m = 1, NSPC1
2436! print *,' *** CHECK AEROSOLS OPTICAL DEPTH FOR SPECIES:', &
2437! & m
2438! print *, aerodp(:,m)
2439! sumodp(:) = sumodp(:) + aerodp(:,m)
2440! enddo
2441
2442!
2443! print *,' *** CHECK AEROSOLS OPTICAL DEPTH FOR ALL SPECIES:'
2444! print *, sumodp(:)
2445! endif
2446! do m = 1, NBDLW
2447! print *,' *** CHECK AEROSOLS PROPERTIES FOR LW BAND =',m, &
2448! & ' ***'
2449! do k = 1, 10
2450! print *,' LEVEL :',k
2451! print *,' TAUAER:',aerolw(:,k,m,1)
2452! print *,' SSAAER:',aerolw(:,k,m,2)
2453! print *,' ASYAER:',aerolw(:,k,m,3)
2454! enddo
2455! enddo
2456
2457 endif ! end if_laswflg_or_lalwflg_block
2458
2466! --- ... stratosphere volcanic forcing
2467
2468 if ( lavoflg ) then
2469
2470 if ( iaerflg == 100 ) then
2471 laddsw = lsswr
2472 laddlw = lslwr
2473 else
2474 laddsw = lsswr .and. laswflg
2475 laddlw = lslwr .and. lalwflg
2476 endif
2477
2478 i1 = mod(kyrsav, 10) + 1
2479
2480! --- select data in 4 lat bands, interpolation at the boundaires
2481
2482 do i = 1, imax
2483 if ( alat(i) > 46.0 ) then
2484 volcae(i) = 1.0e-4 * ivolae(kmonsav,1,i1)
2485 else if ( alat(i) > 44.0 ) then
2486 volcae(i) = 5.0e-5 &
2487 & * (ivolae(kmonsav,1,i1) + ivolae(kmonsav,2,i1))
2488 else if ( alat(i) > 1.0 ) then
2489 volcae(i) = 1.0e-4 * ivolae(kmonsav,2,i1)
2490 else if ( alat(i) > -1.0 ) then
2491 volcae(i) = 5.0e-5 &
2492 & * (ivolae(kmonsav,2,i1) + ivolae(kmonsav,3,i1))
2493 else if ( alat(i) >-44.0 ) then
2494 volcae(i) = 1.0e-4 * ivolae(kmonsav,3,i1)
2495 else if ( alat(i) >-46.0 ) then
2496 volcae(i) = 5.0e-5 &
2497 & * (ivolae(kmonsav,3,i1) + ivolae(kmonsav,4,i1))
2498 else
2499 volcae(i) = 1.0e-4 * ivolae(kmonsav,4,i1)
2500 endif
2501 enddo
2502
2503 if (top_at_1) then ! input data from toa to sfc
2504
2505! --- find lower boundary of stratosphere
2506
2507 do i = 1, imax
2508
2509 tmp1 = abs( alat(i) )
2510 if ( tmp1 > 70.0 ) then ! polar, fixed at 25000pa (250mb)
2511 psrfl = 250.0
2512 elseif ( tmp1 < 20.0 ) then ! tropic, fixed at 15000pa (150mb)
2513 psrfl = 150.0
2514 else ! mid-lat, interpolation
2515 psrfl = 110.0 + 2.0*tmp1
2516 endif
2517
2518 kcuth(i) = nlay - 1
2519 kcutl(i) = 2
2520 rdelp(i) = f_one / prsi(i,2)
2521
2522 lab_do_kcuth0 : do k = 2, nlay-2
2523 if ( prsi(i,k) >= psrfh ) then
2524 kcuth(i) = k - 1
2525 exit lab_do_kcuth0
2526 endif
2527 enddo lab_do_kcuth0
2528
2529 lab_do_kcutl0 : do k = 2, nlay-2
2530 if ( prsi(i,k) >= psrfl ) then
2531 kcutl(i) = k - 1
2532 rdelp(i) = f_one / (prsi(i,k) - prsi(i,kcuth(i)))
2533 exit lab_do_kcutl0
2534 endif
2535 enddo lab_do_kcutl0
2536 enddo
2537
2538! --- sw: add volcanic aerosol optical depth to the background value
2539
2540 if ( laddsw ) then
2541 do m = 1, nbdsw
2542 mb = nswstr + m - 1
2543
2544 if ( wvn_sw1(mb) > 20000 ) then ! range of wvlth < 0.5mu
2545 tmp2 = 0.74
2546 elseif ( wvn_sw2(mb) < 20000 ) then ! range of wvlth > 0.5mu
2547 tmp2 = 1.14
2548 else ! range of wvlth in btwn
2549 tmp2 = 0.94
2550 endif
2551 tmp1 = (0.275e-4 * (wvn_sw2(mb)+wvn_sw1(mb))) ** tmp2
2552
2553 do i = 1, imax
2554 kh = kcuth(i)
2555 kl = kcutl(i)
2556 do k = kh, kl
2557 tmp2 = tmp1 * ((prsi(i,k+1) - prsi(i,k)) * rdelp(i))
2558 aerosw(i,k,m,1) = aerosw(i,k,m,1) + tmp2*volcae(i)
2559 enddo
2560
2561! --- smoothing profile at boundary if needed
2562
2563 if ( aerosw(i,kl,m,1) > 10.*aerosw(i,kl+1,m,1) ) then
2564 tmp2 = aerosw(i,kl,m,1) + aerosw(i,kl+1,m,1)
2565 aerosw(i,kl ,m,1) = 0.8 * tmp2
2566 aerosw(i,kl+1,m,1) = 0.2 * tmp2
2567 endif
2568 enddo ! end do_i_block
2569 enddo ! end do_m_block
2570
2571! --- check print
2572! do i = 1, IMAX
2573! print *,' LEV PRESS TAU FOR PROFILE:',i, &
2574! & ' KCUTH, KCUTL =',kcuth(i),kcutl(i)
2575! kh = kcuth(i) - 1
2576! kl = kcutl(i) + 10
2577! do k = kh, kl
2578! write(6,71) k, prsl(i,k), aerosw(i,k,1,1)
2579! 71 format(i3,2e11.4)
2580! enddo
2581! enddo
2582
2583 endif ! end if_laddsw_block
2584
2585! --- lw: add volcanic aerosol optical depth to the background value
2586
2587 if ( laddlw ) then
2588 if ( nlwbnd == 1 ) then
2589
2590 tmp1 = (0.55 / 11.0) ** 1.2
2591 do i = 1, imax
2592 kh = kcuth(i)
2593 kl = kcutl(i)
2594 do k = kh, kl
2595 tmp2 = tmp1 * ((prsi(i,k+1) - prsi(i,k)) * rdelp(i)) &
2596 & * volcae(i)
2597 do m = 1, nbdlw
2598 aerolw(i,k,m,1) = aerolw(i,k,m,1) + tmp2
2599 enddo
2600 enddo
2601 enddo ! end do_i_block
2602
2603 else
2604
2605 do m = 1, nbdlw
2606 tmp1 = (0.275e-4 * (wvn_lw2(m) + wvn_lw1(m))) ** 1.2
2607
2608 do i = 1, imax
2609 kh = kcuth(i)
2610 kl = kcutl(i)
2611 do k = kh, kl
2612 tmp2 = tmp1 * ((prsi(i,k+1)-prsi(i,k)) * rdelp(i))
2613 aerolw(i,k,m,1) = aerolw(i,k,m,1) + tmp2*volcae(i)
2614 enddo
2615 enddo ! end do_i_block
2616 enddo ! end do_m_block
2617
2618 endif ! end if_NLWBND_block
2619 endif ! end if_laddlw_block
2620
2621 else ! input data from sfc to toa
2622
2623! --- find lower boundary of stratosphere
2624
2625 do i = 1, imax
2626
2627 tmp1 = abs( alat(i) )
2628 if ( tmp1 > 70.0 ) then ! polar, fixed at 25000pa (250mb)
2629 psrfl = 250.0
2630 elseif ( tmp1 < 20.0 ) then ! tropic, fixed at 15000pa (150mb)
2631 psrfl = 150.0
2632 else ! mid-lat, interpolation
2633 psrfl = 110.0 + 2.0*tmp1
2634 endif
2635
2636 kcuth(i) = 2
2637 kcutl(i) = nlay - 1
2638 rdelp(i) = f_one / prsi(i,nlay-1)
2639
2640 lab_do_kcuth1 : do k = nlay-1, 2, -1
2641 if ( prsi(i,k) >= psrfh ) then
2642 kcuth(i) = k
2643 exit lab_do_kcuth1
2644 endif
2645 enddo lab_do_kcuth1
2646
2647 lab_do_kcutl1 : do k = nlay, 2, -1
2648 if ( prsi(i,k) >= psrfl ) then
2649 kcutl(i) = k
2650 rdelp(i) = f_one / (prsi(i,k) - prsi(i,kcuth(i)+1))
2651 exit lab_do_kcutl1
2652 endif
2653 enddo lab_do_kcutl1
2654 enddo
2655
2656! --- sw: add volcanic aerosol optical depth to the background value
2657
2658 if ( laddsw ) then
2659 do m = 1, nbdsw
2660 mb = nswstr + m - 1
2661
2662 if ( wvn_sw1(mb) > 20000 ) then ! range of wvlth < 0.5mu
2663 tmp2 = 0.74
2664 elseif ( wvn_sw2(mb) < 20000 ) then ! range of wvlth > 0.5mu
2665 tmp2 = 1.14
2666 else ! range of wvlth in btwn
2667 tmp2 = 0.94
2668 endif
2669 tmp1 = (0.275e-4 * (wvn_sw2(mb)+wvn_sw1(mb))) ** tmp2
2670
2671 do i = 1, imax
2672 kh = kcuth(i)
2673 kl = kcutl(i)
2674 do k = kl, kh
2675 tmp2 = tmp1 * ((prsi(i,k) - prsi(i,k+1)) * rdelp(i))
2676 aerosw(i,k,m,1) = aerosw(i,k,m,1) + tmp2*volcae(i)
2677 enddo
2678
2679! --- smoothing profile at boundary if needed
2680
2681 if ( aerosw(i,kl,m,1) > 10.*aerosw(i,kl-1,m,1) ) then
2682 tmp2 = aerosw(i,kl,m,1) + aerosw(i,kl-1,m,1)
2683 aerosw(i,kl ,m,1) = 0.8 * tmp2
2684 aerosw(i,kl-1,m,1) = 0.2 * tmp2
2685 endif
2686 enddo ! end do_i_block
2687 enddo ! end do_m_block
2688
2689! --- check print
2690! do i = 1, IMAX
2691! print *,' LEV PRESS TAU FOR PROFILE:',i, &
2692! & ' KCUTH, KCUTL =',kcuth(i),kcutl(i)
2693! kh = kcuth(i) + 1
2694! kl = kcutl(i) - 10
2695! do k = kh, kl, -1
2696! write(6,71) NLP1-k,prsl(i,k),aerosw(i,k,1,1)
2697! enddo
2698! enddo
2699
2700 endif ! end if_laddsw_block
2701
2702! --- lw: add volcanic aerosol optical depth to the background value
2703
2704 if ( laddlw ) then
2705 if ( nlwbnd == 1 ) then
2706
2707 tmp1 = (0.55 / 11.0) ** 1.2
2708 do i = 1, imax
2709 kh = kcuth(i)
2710 kl = kcutl(i)
2711 do k = kl, kh
2712 tmp2 = tmp1 * ((prsi(i,k) - prsi(i,k+1)) * rdelp(i)) &
2713 & * volcae(i)
2714 do m = 1, nbdlw
2715 aerolw(i,k,m,1) = aerolw(i,k,m,1) + tmp2
2716 enddo
2717 enddo
2718 enddo ! end do_i_block
2719
2720 else
2721
2722 do m = 1, nbdlw
2723 tmp1 = (0.275e-4 * (wvn_lw2(m) + wvn_lw1(m))) ** 1.2
2724
2725 do i = 1, imax
2726 kh = kcuth(i)
2727 kl = kcutl(i)
2728 do k = kl, kh
2729 tmp2 = tmp1 * ((prsi(i,k)-prsi(i,k+1)) * rdelp(i))
2730 aerolw(i,k,m,1) = aerolw(i,k,m,1) + tmp2*volcae(i)
2731 enddo
2732 enddo ! end do_i_block
2733 enddo ! end do_m_block
2734
2735 endif ! end if_NLWBND_block
2736 endif ! end if_laddlw_block
2737
2738 endif ! end if_top_at_1_block
2739
2740 endif ! end if_lavoflg_block
2741!
2742 return
2743!...................................
2744 end subroutine setaer
2745!-----------------------------------
2746
2774!-----------------------------------
2775 subroutine aer_property &
2776 & ( prsi,prsl,prslk,tvly,rhlay,dz,hz,tracer, & ! --- inputs:
2777 & alon,alat,slmsk, laersw,laerlw, &
2778 & imax,nlay,nlp1,top_at_1, &
2779 & aerosw,aerolw,aerodp,errflg,errmsg & ! --- outputs:
2780 & )
2781
2782! ================================================================== !
2783! !
2784! aer_property maps the 5 degree global climatological aerosol data !
2785! set onto model grids, and compute aerosol optical properties for sw !
2786! and lw radiations. !
2787! !
2788! inputs: !
2789! prsi - pressure at interface mb IMAX*NLP1 !
2790! prsl - layer mean pressure (not used) IMAX*NLAY !
2791! prslk - exner function=(p/p0)**rocp (not used) IMAX*NLAY !
2792! tvly - layer virtual temperature (not used) IMAX*NLAY !
2793! rhlay - layer mean relative humidity IMAX*NLAY !
2794! dz - layer thickness m IMAX*NLAY !
2795! hz - level high m IMAX*NLP1 !
2796! tracer - aer tracer concentrations (not used) IMAX*NLAY*NTRAC!
2797! alon, alat IMAX !
2798! - longitude and latitude of given points in degree !
2799! slmsk - sea/land mask (sea:0,land:1,sea-ice:2) IMAX !
2800! laersw,laerlw 1 !
2801! - logical flag for sw/lw aerosol calculations !
2802! IMAX - horizontal dimension of arrays 1 !
2803! NLAY,NLP1-vertical dimensions of arrays 1 !
2804!! NSPC - num of species for optional aod output fields 1 !
2805! top_at_1 - vertical ordering flag !
2806! !
2807! outputs: !
2808! aerosw - aeros opt properties for sw IMAX*NLAY*NBDSW*NF_AESW!
2809! (:,:,:,1): optical depth !
2810! (:,:,:,2): single scattering albedo !
2811! (:,:,:,3): asymmetry parameter !
2812! aerolw - aeros opt properties for lw IMAX*NLAY*NBDLW*NF_AELW!
2813! (:,:,:,1): optical depth !
2814! (:,:,:,2): single scattering albedo !
2815! (:,:,:,3): asymmetry parameter !
2816!! aerodp - vertically integrated aer-opt-depth IMAX*NSPC+1 !
2817! !
2818! errflg - CCPP error flag !
2819! errmsg - CCPP error message !
2820! !
2821! module parameters and constants: !
2822! NSWBND - total number of actual sw spectral bands computed !
2823! NLWBND - total number of actual lw spectral bands computed !
2824! NSWLWBD - total number of sw+lw bands computed !
2825! !
2826! module variable: (set by subroutine aer_init) !
2827! kprfg - aerosols profile index IMXAE*JMXAE !
2828! 1:ant 2:arc 3:cnt 4:mar 5:des 6:marme 7:cntme !
2829! idxcg - aerosols component index NXC*IMXAE*JMXAE !
2830! 1:inso 2:soot 3:minm 4:miam 5:micm !
2831! 6:mitr 7:waso 8:ssam 9:sscm 10:suso !
2832! cmixg - aerosols component mixing ratio NXC*IMXAE*JMXAE !
2833! denng - aerosols number density 2 *IMXAE*JMXAE !
2834! 1:for domain-1 2:domain-2 (prof marme/cntme only) !
2835! !
2836! usage: call aer_property !
2837! !
2838! subprograms called: radclimaer !
2839! !
2840! ================================================================== !
2841
2842! --- inputs:
2843 integer, intent(in) :: IMAX, NLAY, NLP1
2844! integer, intent(in) :: IMAX, NLAY, NLP1, NSPC
2845 logical, intent(in) :: laersw, laerlw, top_at_1
2846
2847 real (kind=kind_phys), dimension(:,:), intent(in) :: prsi, prsl, &
2848 & prslk, tvly, rhlay, dz, hz
2849 real (kind=kind_phys), dimension(:), intent(in) :: alon, alat, &
2850 & slmsk
2851 real (kind=kind_phys), dimension(:,:,:),intent(in):: tracer
2852
2853! --- outputs:
2854 real (kind=kind_phys), dimension(:,:,:,:), intent(out) :: &
2855 & aerosw, aerolw
2856 real (kind=kind_phys), dimension(:,:) , intent(out) :: aerodp
2857 integer, intent(out) :: errflg
2858 character(len=*), intent(out) :: errmsg
2859
2860! --- locals:
2861 real (kind=kind_phys), dimension(NCM) :: cmix
2862 real (kind=kind_phys), dimension( 2) :: denn
2863 real (kind=kind_phys), dimension(NSPC) :: spcodp
2864
2865 real (kind=kind_phys), dimension(NLAY) :: delz, rh1, dz1
2866 integer, dimension(NLAY) :: idmaer
2867
2868 real (kind=kind_phys), dimension(NLAY,NSWLWBD):: tauae,ssaae,asyae
2869!test real (kind=kind_phys), dimension(IMAX,NLAY) :: aersav
2870
2871 real (kind=kind_phys) :: tmp1, tmp2, rps, dtmp, h1
2872 real (kind=kind_phys) :: wi, wj, w11, w12, w21, w22
2873
2874 integer :: i, ii, i1, i2, i3, j1, j2, j3, k, m, m1, &
2875 & kp, kpa, kpi, kpj
2876
2877! --- conversion constants
2878 real (kind=kind_phys), parameter :: dltg = 360.0 / float(imxae)
2879 real (kind=kind_phys), parameter :: hdlt = 0.5 * dltg
2880 real (kind=kind_phys), parameter :: rdlt = 1.0 / dltg
2881
2882!
2883!===> ... begin here
2884!
2885
2886! Initialize CCPP error handling variables
2887 errmsg = ''
2888 errflg = 0
2889
2893
2894 i1 = 1
2895 i2 = 2
2896 j1 = 1
2897 j2 = 2
2898
2899 lab_do_imax : do i = 1, imax
2900
2901! --- map grid in longitude direction, lon from 0 to 355 deg resolution
2902
2903! print *,' Seeking lon index for point i =',i
2904 i3 = i1
2905 lab_do_imxae : do while ( i3 <= imxae )
2906 tmp1 = dltg * (i3 - 1)
2907 dtmp = alon(i) - tmp1
2908! print *,' alon, i3, tlon, dlon =',alon(i),i3,tmp1,dtmp
2909
2910 if ( dtmp > dltg ) then
2911 i3 = i3 + 1
2912 if ( i3 > imxae ) then
2913 print *,' ERROR! In setclimaer alon>360. ipt =',i, &
2914 & ', dltg,alon,tlon,dlon =',dltg,alon(i),tmp1,dtmp
2915 errflg = 1
2916 errmsg = 'ERROR(aer_property)'
2917 return
2918 endif
2919 elseif ( dtmp >= f_zero ) then
2920 i1 = i3
2921 i2 = mod(i3,imxae) + 1
2922 wi = dtmp * rdlt
2923 if ( dtmp <= hdlt ) then
2924 kpi = i3
2925 else
2926 kpi = i2
2927 endif
2928! print *,' found i1, i2, wi =',i1,i2,wi
2929 exit lab_do_imxae
2930 else
2931 i3 = i3 - 1
2932 if ( i3 < 1 ) then
2933 print *,' ERROR! In setclimaer alon< 0. ipt =',i, &
2934 & ', dltg,alon,tlon,dlon =',dltg,alon(i),tmp1,dtmp
2935 errflg = 1
2936 errmsg = 'ERROR(aer_property)'
2937 return
2938 endif
2939 endif
2940 enddo lab_do_imxae
2941
2942! --- map grid in latitude direction, lat from 90n to 90s in 5 deg resolution
2943
2944! print *,' Seeking lat index for point i =',i
2945 j3 = j1
2946 lab_do_jmxae : do while ( j3 <= jmxae )
2947 tmp2 = 90.0 - dltg * (j3 - 1)
2948 dtmp = tmp2 - alat(i)
2949! print *,' alat, j3, tlat, dlat =',alat(i),j3,tmp2,dtmp
2950
2951 if ( dtmp > dltg ) then
2952 j3 = j3 + 1
2953 if ( j3 >= jmxae ) then
2954 print *,' ERROR! In setclimaer alat<-90. ipt =',i, &
2955 & ', dltg,alat,tlat,dlat =',dltg,alat(i),tmp2,dtmp
2956 errflg = 1
2957 errmsg = 'ERROR(aer_property)'
2958 return
2959 endif
2960 elseif ( dtmp >= f_zero ) then
2961 j1 = j3
2962 j2 = j3 + 1
2963 wj = dtmp * rdlt
2964 if ( dtmp <= hdlt ) then
2965 kpj = j3
2966 else
2967 kpj = j2
2968 endif
2969! print *,' found j1, j2, wj =',j1,j2,wj
2970 exit lab_do_jmxae
2971 else
2972 j3 = j3 - 1
2973 if ( j3 < 1 ) then
2974 print *,' ERROR! In setclimaer alat>90. ipt =',i, &
2975 & ', dltg,alat,tlat,dlat =',dltg,alat(i),tmp2,dtmp
2976 errflg = 1
2977 errmsg = 'ERROR(aer_property)'
2978 return
2979 endif
2980 endif
2981 enddo lab_do_jmxae
2982
2985
2986 kp = kprfg(kpi,kpj) ! nearest typical aeros profile as default
2987 kpa = max( kprfg(i1,j1),kprfg(i1,j2),kprfg(i2,j1),kprfg(i2,j2) )
2988 h1 = haer(1,kp)
2989 denn(2) = f_zero
2990 ii = 1
2991
2992 if ( kp /= kpa ) then
2993 if ( kpa == 6 ) then ! if ocean prof with mineral aeros overlay
2994 ii = 2 ! need 2 types of densities
2995 if ( slmsk(i) > f_zero ) then ! but actually a land/sea-ice point
2996 kp = 7 ! reset prof index to land
2997 h1 = 0.5*(haer(1,6) + haer(1,7)) ! use a transition scale hight
2998 else
2999 kp = kpa
3000 h1 = haer(1,6)
3001 endif
3002 elseif ( kpa == 7 ) then ! if land prof with mineral aeros overlay
3003 ii = 2 ! need 2 types of densities
3004 if ( slmsk(i) <= f_zero ) then ! but actually an ocean point
3005 kp = 6 ! reset prof index to ocean
3006 h1 = 0.5*(haer(1,6) + haer(1,7)) ! use a transition scale hight
3007 else
3008 kp = kpa
3009 h1 = haer(1,7)
3010 endif
3011 else ! lower atmos without mineral aeros overlay
3012! h1 = 0.5*(haer(1,kp) + haer(1,kpa)) ! use a transition scale hight
3013 h1 = haer(1,kpa)
3014 kp = kpa
3015 endif
3016 endif
3017
3019
3020 w11 = (f_one-wi) * (f_one-wj)
3021 w12 = (f_one-wi) * wj
3022 w21 = wi * (f_one-wj)
3023 w22 = wi * wj
3024
3025! --- check print
3026! print *,' Grid pt', i,', alon, alat =',alon(i),alat(i), &
3027! & ', tlon, tlat =',tmp1,tmp2
3028! print *,' lon grid index i1, i2 =',i1,i2,', weight wi =',wi
3029! print *,' lat grid index j1, j2 =',j1,j2,', weight wj =',wj
3030! print *,' bi-linear weights w11,w21,w12,w22 =',w11,w21,w12,w22
3031! print *,' kp,kpa,slmsk,h1 =',kp,m1,slmsk(i),h1
3032
3035
3036 do m = 1, ii ! ii=1 for domain 1; =2 for domain 2.
3037 denn(m) = w11*denng(m,i1,j1) + w12*denng(m,i1,j2) &
3038 & + w21*denng(m,i2,j1) + w22*denng(m,i2,j2)
3039 enddo ! end_do_m_loop
3040
3042
3043 cmix(:) = f_zero
3044 do m = 1, nxc
3045 ii = idxcg(m,i1,j1)
3046 if ( ii > 0 ) then
3047 cmix(ii) = cmix(ii) + w11*cmixg(m,i1,j1)
3048 endif
3049 ii = idxcg(m,i1,j2)
3050 if ( ii > 0 ) then
3051 cmix(ii) = cmix(ii) + w12*cmixg(m,i1,j2)
3052 endif
3053 ii = idxcg(m,i2,j1)
3054 if ( ii > 0 ) then
3055 cmix(ii) = cmix(ii) + w21*cmixg(m,i2,j1)
3056 endif
3057 ii = idxcg(m,i2,j2)
3058 if ( ii > 0 ) then
3059 cmix(ii) = cmix(ii) + w22*cmixg(m,i2,j2)
3060 endif
3061 enddo ! end_do_m_loop
3062
3063! --- check print
3064! print *,' denn =',denn(:)
3065! print *,' cmix =',cmix(:)
3066
3069
3070 do k = 1, nlay
3071 rh1(k) = rhlay(i,k)
3072 dz1(k) = dz(i,k)
3073 enddo
3074
3075 lab_if_flip : if (.not. top_at_1) then ! input from sfc to toa
3076
3077 if ( prsi(i,1) > 100.0 ) then
3078 rps = f_one / prsi(i,1)
3079 else
3080 print *,' !!! (1) Error in subr radiation_aerosols:', &
3081 & ' unrealistic surface pressure =', i,prsi(i,1)
3082 errflg = 1
3083 errmsg = 'ERROR(aer_property): Unrealistic surface pressure'
3084 return
3085 endif
3086
3087 ii = 1
3088 do k = 1, nlay
3089 if (prsi(i,k+1)*rps < sigref(ii,kp)) then
3090 ii = ii + 1
3091 if (ii == 2 .and. prsref(2,kp) == prsref(3,kp)) then
3092 ii = 3
3093 endif
3094 endif
3095 idmaer(k) = ii
3096
3097 if ( ii > 1 ) then
3098 tmp1 = haer(ii,kp)
3099 else
3100 tmp1 = h1
3101 endif
3102
3103 if (tmp1 > f_zero) then
3104 tmp2 = f_one / tmp1
3105 delz(k) = tmp1 * (exp(-hz(i,k)*tmp2)-exp(-hz(i,k+1)*tmp2))
3106 else
3107 delz(k) = dz1(k)
3108 endif
3109 enddo
3110
3111 else lab_if_flip ! input from toa to sfc
3112
3113 if ( prsi(i,nlp1) > 100.0 ) then
3114 rps = 1.0 / prsi(i,nlp1)
3115 else
3116 print *,' !!! (2) Error in subr radiation_aerosols:', &
3117 & ' unrealistic surface pressure =', i,prsi(i,nlp1)
3118 endif
3119
3120 ii = 1
3121 do k = nlay, 1, -1
3122 if (prsi(i,k)*rps < sigref(ii,kp)) then
3123 ii = ii + 1
3124 if (ii == 2 .and. prsref(2,kp) == prsref(3,kp)) then
3125 ii = 3
3126 endif
3127 endif
3128 idmaer(k) = ii
3129
3130 if ( ii > 1 ) then
3131 tmp1 = haer(ii,kp)
3132 else
3133 tmp1 = h1
3134 endif
3135
3136 if (tmp1 > f_zero) then
3137 tmp2 = f_one / tmp1
3138 delz(k) = tmp1 * (exp(-hz(i,k+1)*tmp2)-exp(-hz(i,k)*tmp2))
3139 else
3140 delz(k) = dz1(k)
3141 endif
3142 enddo
3143
3144 endif lab_if_flip
3145
3146! --- check print
3147
3148! print *,' in setclimaer, profile:',i
3149! print *,' rh :',rh1
3150! print *,' dz :',dz1
3151! print *,' delz :',delz
3152! print *,' idmaer:',idmaer
3153
3156
3157 call radclimaer(top_at_1)
3158! --- inputs: (in-scope variables)
3159! --- outputs: (in-scope variables)
3160
3161 if ( laersw ) then
3162
3163 do m = 1, nbdsw
3164 do k = 1, nlay
3165 aerosw(i,k,m,1) = tauae(k,m)
3166 aerosw(i,k,m,2) = ssaae(k,m)
3167 aerosw(i,k,m,3) = asyae(k,m)
3168 enddo
3169 enddo
3170
3171! --- total aod (optional)
3172 do k = 1, nlay
3173 aerodp(i,1) = aerodp(i,1) + tauae(k,nv_aod)
3174 enddo
3175
3176! --- for diagnostic output (optional)
3177 do m = 1, nspc
3178 aerodp(i,m+1) = spcodp(m)
3179 enddo
3180
3181 endif ! end if_larsw_block
3182
3183 if ( laerlw ) then
3184
3185 if ( nlwbnd == 1 ) then
3186 m1 = nswbnd + 1
3187 do m = 1, nbdlw
3188 do k = 1, nlay
3189 aerolw(i,k,m,1) = tauae(k,m1)
3190 aerolw(i,k,m,2) = ssaae(k,m1)
3191 aerolw(i,k,m,3) = asyae(k,m1)
3192 enddo
3193 enddo
3194 else
3195 do m = 1, nbdlw
3196 m1 = nswbnd + m
3197 do k = 1, nlay
3198 aerolw(i,k,m,1) = tauae(k,m1)
3199 aerolw(i,k,m,2) = ssaae(k,m1)
3200 aerolw(i,k,m,3) = asyae(k,m1)
3201 enddo
3202 enddo
3203 endif
3204
3205 endif ! end if_laerlw_block
3206
3207 enddo lab_do_imax
3208
3209! =================
3210 contains
3211! =================
3212
3217!--------------------------------
3218 subroutine radclimaer(top_at_1)
3219!................................
3220
3221! --- inputs: (in scope variables)
3222! --- outputs: (in scope variables)
3223
3224! ================================================================== !
3225! !
3226! compute aerosols optical properties in NSWLWBD bands. there are !
3227! seven different vertical profile structures. in the troposphere, !
3228! aerosol distribution at each grid point is composed from up to !
3229! six components out of a total of ten different substances. !
3230! !
3231! ref: wmo report wcp-112 (1986) !
3232! !
3233! input variables: !
3234! cmix - mixing ratioes of aerosol components - NCM !
3235! denn - aerosol number densities - 2 !
3236! rh1 - relative humidity - NLAY !
3237! delz - effective layer thickness km NLAY !
3238! idmaer - aerosol domain index - NLAY !
3239! NXC - number of different aerosol components- 1 !
3240! NLAY - vertical dimensions - 1 !
3241! !
3242! output variables: !
3243! tauae - optical depth - NLAY*NSWLWBD!
3244! ssaae - single scattering albedo - NLAY*NSWLWBD!
3245! asyae - asymmetry parameter - NLAY*NSWLWBD!
3246!! aerodp - vertically integrated aer-opt-depth - IMAX*NSPC+1 !
3247! !
3248! ================================================================== !
3249!
3250 real (kind=kind_phys) :: crt1, crt2
3251 parameter(crt1=30.0, crt2=0.03333)
3252
3253! --- inputs:
3254 logical, intent(in) :: top_at_1
3255! --- outputs:
3256
3257! --- locals:
3258 real (kind=kind_phys) :: cm, hd, hdi, sig0u, sig0l, ratio, tt0, &
3259 & ex00, sc00, ss00, as00, ex01, sc01, ss01, as01, tt1, &
3260 & ex02, sc02, ss02, as02, ex03, sc03, ss03, as03, tt2, &
3261 & ext1, sca1, ssa1, asy1, drh0, drh1, rdrh
3262
3263 integer :: ih1, ih2, kk, idom, icmp, ib, ii, ic, ic1
3264 integer :: idx
3265
3266!===> ... begin here
3267
3268 spcodp = f_zero
3269
3270!===> ... loop over vertical layers from top to surface
3271
3272 lab_do_layer : do kk = 1, nlay
3273
3274! --- linear interp coeffs for rh-dep species
3275
3276 ih2 = 1
3277 do while ( rh1(kk) > rhlev(ih2) )
3278 ih2 = ih2 + 1
3279 if ( ih2 > nrhlev ) exit
3280 enddo
3281 ih1 = max( 1, ih2-1 )
3282 ih2 = min( nrhlev, ih2 )
3283
3284 drh0 = rhlev(ih2) - rhlev(ih1)
3285 drh1 = rh1(kk) - rhlev(ih1)
3286 if ( ih1 == ih2 ) then
3287 rdrh = f_zero
3288 else
3289 rdrh = drh1 / drh0
3290 endif
3291
3292! --- assign optical properties in each domain
3293
3294 idom = idmaer(kk)
3295
3296 lab_if_idom : if (idom == 5) then
3297! --- 5th domain - upper stratosphere assume no aerosol
3298
3299 do ib = 1, nswlwbd
3300 tauae(kk,ib) = f_zero
3301 if ( ib <= nswbnd ) then
3302 ssaae(kk,ib) = 0.99
3303 asyae(kk,ib) = 0.696
3304 else
3305 ssaae(kk,ib) = 0.5
3306 asyae(kk,ib) = 0.3
3307 endif
3308 enddo
3309
3310 elseif (idom == 4) then lab_if_idom
3311! --- 4th domain - stratospheric layers
3312
3313 do ib = 1, nswlwbd
3314 tauae(kk,ib) = extstra(ib) * delz(kk)
3315 if ( ib <= nswbnd ) then
3316 ssaae(kk,ib) = 0.99
3317 asyae(kk,ib) = 0.696
3318 else
3319 ssaae(kk,ib) = 0.5
3320 asyae(kk,ib) = 0.3
3321 endif
3322 enddo
3323
3324! --- compute aod from individual species' contribution (optional)
3325 idx = idxspc(10) ! for sulfate
3326 spcodp(idx) = spcodp(idx) + tauae(kk,nv_aod)
3327
3328 elseif (idom == 3) then lab_if_idom
3329! --- 3rd domain - free tropospheric layers
3330! 1:inso 0.17e-3; 2:soot 0.4; 7:waso 0.59983; n:730
3331
3332 do ib = 1, nswlwbd
3333 ex01 = extrhi(1,ib)
3334 sc01 = scarhi(1,ib)
3335 ss01 = ssarhi(1,ib)
3336 as01 = asyrhi(1,ib)
3337
3338 ex02 = extrhi(2,ib)
3339 sc02 = scarhi(2,ib)
3340 ss02 = ssarhi(2,ib)
3341 as02 = asyrhi(2,ib)
3342
3343 ex03 = extrhd(ih1,1,ib) &
3344 & + rdrh * (extrhd(ih2,1,ib) - extrhd(ih1,1,ib))
3345 sc03 = scarhd(ih1,1,ib) &
3346 & + rdrh * (scarhd(ih2,1,ib) - scarhd(ih1,1,ib))
3347 ss03 = ssarhd(ih1,1,ib) &
3348 & + rdrh * (ssarhd(ih2,1,ib) - ssarhd(ih1,1,ib))
3349 as03 = asyrhd(ih1,1,ib) &
3350 & + rdrh * (asyrhd(ih2,1,ib) - asyrhd(ih1,1,ib))
3351
3352 ext1 = 0.17e-3*ex01 + 0.4*ex02 + 0.59983*ex03
3353 sca1 = 0.17e-3*sc01 + 0.4*sc02 + 0.59983*sc03
3354 ssa1 = 0.17e-3*ss01*ex01 + 0.4*ss02*ex02 + 0.59983*ss03*ex03
3355 asy1 = 0.17e-3*as01*sc01 + 0.4*as02*sc02 + 0.59983*as03*sc03
3356
3357 tauae(kk,ib) = ext1 * 730.0 * delz(kk)
3358 ssaae(kk,ib) = min(f_one, ssa1/ext1)
3359 asyae(kk,ib) = min(f_one, asy1/sca1)
3360
3361! --- compute aod from individual species' contribution (optional)
3362 if ( ib==nv_aod ) then
3363 spcodp(1) = spcodp(1) + 0.17e-3*ex01*730.0*delz(kk) ! dust (inso) #1
3364 spcodp(2) = spcodp(2) + 0.4 *ex02*730.0*delz(kk) ! black carbon #2
3365 spcodp(3) = spcodp(3) + 0.59983*ex03*730.0*delz(kk) ! water soluble #7
3366 endif
3367
3368 enddo
3369
3370 elseif (idom == 1) then lab_if_idom
3371! --- 1st domain - mixing layer
3372
3373 lab_do_ib : do ib = 1, nswlwbd
3374 ext1 = f_zero
3375 sca1 = f_zero
3376 ssa1 = f_zero
3377 asy1 = f_zero
3378
3379 lab_do_icmp : do icmp = 1, ncm
3380 ic = icmp
3381 idx = idxspc(icmp)
3382
3383 cm = cmix(icmp)
3384 lab_if_cm : if ( cm > f_zero ) then
3385
3386 lab_if_ic : if ( ic <= ncm1 ) then ! component withour rh dep
3387 tt0 = cm * extrhi(ic,ib)
3388 ext1 = ext1 + tt0
3389 sca1 = sca1 + cm * scarhi(ic,ib)
3390 ssa1 = ssa1 + cm * ssarhi(ic,ib) * extrhi(ic,ib)
3391 asy1 = asy1 + cm * asyrhi(ic,ib) * scarhi(ic,ib)
3392 else lab_if_ic ! component with rh dep
3393 ic1 = ic - ncm1
3394
3395 ex00 = extrhd(ih1,ic1,ib) &
3396 & + rdrh * (extrhd(ih2,ic1,ib) - extrhd(ih1,ic1,ib))
3397 sc00 = scarhd(ih1,ic1,ib) &
3398 & + rdrh * (scarhd(ih2,ic1,ib) - scarhd(ih1,ic1,ib))
3399 ss00 = ssarhd(ih1,ic1,ib) &
3400 & + rdrh * (ssarhd(ih2,ic1,ib) - ssarhd(ih1,ic1,ib))
3401 as00 = asyrhd(ih1,ic1,ib) &
3402 & + rdrh * (asyrhd(ih2,ic1,ib) - asyrhd(ih1,ic1,ib))
3403
3404 tt0 = cm * ex00
3405 ext1 = ext1 + tt0
3406 sca1 = sca1 + cm * sc00
3407 ssa1 = ssa1 + cm * ss00 * ex00
3408 asy1 = asy1 + cm * as00 * sc00
3409 endif lab_if_ic
3410
3411! --- compute aod from individual species' contribution (optional)
3412 if ( ib==nv_aod ) then
3413 spcodp(idx) = spcodp(idx) + tt0*denn(1)*delz(kk) ! idx for dif species
3414 endif
3415
3416 endif lab_if_cm
3417 enddo lab_do_icmp
3418
3419 tauae(kk,ib) = ext1 * denn(1) * delz(kk)
3420 ssaae(kk,ib) = min(f_one, ssa1/ext1)
3421 asyae(kk,ib) = min(f_one, asy1/sca1)
3422 enddo lab_do_ib
3423
3424 elseif (idom == 2) then lab_if_idom
3425! --- 2nd domain - mineral transport layers
3426
3427 do ib = 1, nswlwbd
3428 tauae(kk,ib) = extrhi(6,ib) * denn(2) * delz(kk)
3429 ssaae(kk,ib) = ssarhi(6,ib)
3430 asyae(kk,ib) = asyrhi(6,ib)
3431 enddo
3432
3433! --- compute aod from individual species' contribution (optional)
3434 spcodp(1) = spcodp(1) + tauae(kk,nv_aod) ! dust
3435
3436 else lab_if_idom
3437! --- domain index out off range, assume no aerosol
3438
3439 do ib = 1, nswlwbd
3440 tauae(kk,ib) = f_zero
3441 ssaae(kk,ib) = f_one
3442 asyae(kk,ib) = f_zero
3443 enddo
3444
3445! write(6,19) kk,idom
3446! 19 format(/' *** ERROR in sub AEROS: domain index out' &
3447! &, ' of range! K, IDOM =',3i5,' ***')
3448! stop 19
3449
3450 endif lab_if_idom
3451
3452 enddo lab_do_layer
3453
3454!
3455!===> ... smooth profile at domain boundaries
3456!
3457 if (top_at_1) then ! input from toa to sfc
3458
3459 do ib = 1, nswlwbd
3460 do kk = 2, nlay
3461 if ( tauae(kk,ib) > f_zero ) then
3462 ratio = tauae(kk-1,ib) / tauae(kk,ib)
3463 else
3464 ratio = f_one
3465 endif
3466
3467 tt0 = tauae(kk,ib) + tauae(kk-1,ib)
3468 tt1 = 0.2 * tt0
3469 tt2 = tt0 - tt1
3470
3471 if ( ratio > crt1 ) then
3472 tauae(kk,ib) = tt1
3473 tauae(kk-1,ib) = tt2
3474 endif
3475
3476 if ( ratio < crt2 ) then
3477 tauae(kk,ib) = tt2
3478 tauae(kk-1,ib) = tt1
3479 endif
3480 enddo ! do_kk_loop
3481 enddo ! do_ib_loop
3482
3483 else ! input from sfc to toa
3484
3485 do ib = 1, nswlwbd
3486 do kk = nlay-1, 1, -1
3487 if ( tauae(kk,ib) > f_zero ) then
3488 ratio = tauae(kk+1,ib) / tauae(kk,ib)
3489 else
3490 ratio = f_one
3491 endif
3492
3493 tt0 = tauae(kk,ib) + tauae(kk+1,ib)
3494 tt1 = 0.2 * tt0
3495 tt2 = tt0 - tt1
3496
3497 if ( ratio > crt1 ) then
3498 tauae(kk,ib) = tt1
3499 tauae(kk+1,ib) = tt2
3500 endif
3501
3502 if ( ratio < crt2 ) then
3503 tauae(kk,ib) = tt2
3504 tauae(kk+1,ib) = tt1
3505 endif
3506 enddo ! do_kk_loop
3507 enddo ! do_ib_loop
3508
3509 endif
3510
3511!
3512 return
3513!................................
3514 end subroutine radclimaer
3515!--------------------------------
3516!
3517!...................................
3518 end subroutine aer_property
3519!-----------------------------------
3520
3530!-----------------------------------
3531 subroutine gocart_aerinit &
3532 & ( solfwv, eirfwv, me, &
3533 & errflg, errmsg)
3534
3535! ================================================================== !
3536! !
3537! subprogram : gocart_aerinit !
3538! !
3539! gocart_aerinit is the gocart aerosol initialization program !
3540! to set up necessary parameters and working arrays. !
3541! !
3542! inputs: !
3543! solfwv(NWVTOT) - solar flux for each individual wavenumber (w/m2)!
3544! eirfwv(NWVTIR) - ir flux(273k) for each individual wavenum (w/m2)!
3545! me - print message control flag !
3546! !
3547! outputs: (CCPP error handling) !
3548! errflg - CCPP error flag !
3549! errmsg - CCPP error message !
3550! !
3551! module variables: !
3552! NWVSOL - num of wvnum regions where solar flux is constant !
3553! NWVTOT - total num of wave numbers used in sw spectrum !
3554! NWVTIR - total num of wave numbers used in the ir region !
3555! NSWBND - total number of sw spectral bands !
3556! NLWBND - total number of lw spectral bands !
3557! NAERBND - number of bands for climatology aerosol data !
3558! KCM1 - number of rh independent aeros species !
3559! KCM2 - number of rh dependent aeros species !
3560! !
3561! usage: call gocart_init !
3562! !
3563! subprograms called: rd_gocart_luts, optavg_gocart !
3564! !
3565! ================================================================== !
3566
3567 implicit none
3568
3569! --- inputs:
3570 real (kind=kind_phys), dimension(:) :: solfwv ! one wvn sol flux
3571 real (kind=kind_phys), dimension(:) :: eirfwv ! one wvn ir flux
3572
3573 integer, intent(in) :: me
3574
3575! --- output: (CCPP error handling)
3576 integer, intent(out) :: errflg
3577 character(len=*), intent(out) :: errmsg
3578
3579! --- locals:
3580 real (kind=kind_phys), dimension(kaerbndi,kcm1) :: &
3581 & rhidext0_grt, rhidsca0_grt, rhidssa0_grt, rhidasy0_grt
3582 real (kind=kind_phys), dimension(kaerbndd,krhlev,kcm2):: &
3583 & rhdpext0_grt, rhdpsca0_grt, rhdpssa0_grt, rhdpasy0_grt
3584
3585 real (kind=kind_phys), dimension(nswbnd,kaerbndd) :: solwaer
3586 real (kind=kind_phys), dimension(nswbnd) :: solbnd
3587 real (kind=kind_phys), dimension(nlwbnd,kaerbndd) :: eirwaer
3588 real (kind=kind_phys), dimension(nlwbnd) :: eirbnd
3589
3590 real (kind=kind_phys), dimension(nswbnd,kaerbndi) :: solwaer_du
3591 real (kind=kind_phys), dimension(nswbnd) :: solbnd_du
3592 real (kind=kind_phys), dimension(nlwbnd,kaerbndi) :: eirwaer_du
3593 real (kind=kind_phys), dimension(nlwbnd) :: eirbnd_du
3594
3595 integer, dimension(nswbnd) :: nv1, nv2, nv1_du, nv2_du
3596 integer, dimension(nlwbnd) :: nr1, nr2, nr1_du, nr2_du
3597
3598 integer, dimension(kaerbndd) :: iendwv
3599 integer, dimension(kaerbndi) :: iendwv_du
3600 real (kind=kind_phys), dimension(kaerbndd) :: wavelength
3601 real (kind=kind_phys), dimension(kaerbndi) :: wavelength_du
3602 real (kind=kind_phys) :: sumsol, sumir, sumsol_du, sumir_du
3603
3604 integer :: i, j, k, mb, ib, ii, iix, iw, iw1, iw2
3605
3606!
3607!===> ... begin here
3608
3609! Initialize CCPP error handling variables
3610 errmsg = ''
3611 errflg = 0
3612
3613!
3614! --- ... invoke gocart aerosol initialization
3615
3616
3617 if (kcm /= ntrcaerm ) then
3618 print *, 'ERROR in # of gocart aer species',kcm
3619 errflg = 1
3620 errmsg = 'ERROR(gocart_init): Incorrect # of species'
3621 return
3622 endif
3623
3624! --- ... aloocate and input aerosol optical data
3625
3626 if ( .not. allocated( extrhi_grt ) ) then
3627 allocate ( extrhi_grt( kcm1,nswlwbd) )
3628 allocate ( extrhi_grt_550( kcm1,1) )
3629 allocate ( scarhi_grt( kcm1,nswlwbd) )
3630 allocate ( ssarhi_grt( kcm1,nswlwbd) )
3631 allocate ( asyrhi_grt( kcm1,nswlwbd) )
3632 allocate ( extrhd_grt(krhlev,kcm2,nswlwbd) )
3633 allocate ( extrhd_grt_550(krhlev,kcm2,1) )
3634 allocate ( scarhd_grt(krhlev,kcm2,nswlwbd) )
3635 allocate ( ssarhd_grt(krhlev,kcm2,nswlwbd) )
3636 allocate ( asyrhd_grt(krhlev,kcm2,nswlwbd) )
3637 endif
3638
3639! --- ... read tabulated GOCART aerosols optical data
3640
3641 call rd_gocart_luts
3642! --- inputs: (in scope variables, module variables)
3643! --- outputs: (in scope variables)
3644
3645! --- ... convert wavelength to wavenumber
3646! wavelength and wavelength_du are read-in by rd_gocart_luts
3647
3648 do i = 1, kaerbndd
3649 iendwv(i) = int(10000. / wavelength(i))
3650 enddo
3651
3652 do i = 1, kaerbndi
3653 iendwv_du(i) = int(10000. / wavelength_du(i))
3654 enddo
3655
3656! --- ... compute solar flux weights and interval indices for mapping
3657! spectral bands between sw radiation and aerosol data
3658
3659 if ( laswflg ) then
3660 solbnd(:) = f_zero
3661 solbnd_du(:)= f_zero
3662 do i=1,nswbnd
3663 do j=1,kaerbndd
3664 solwaer(i,j) = f_zero
3665 enddo
3666 do j=1,kaerbndi
3667 solwaer_du(i,j) = f_zero
3668 enddo
3669 enddo
3670
3671 do ib = 1, nswbnd
3672 mb = ib + nswstr - 1
3673 ii = 1
3674 iix = 1
3675 iw1 = nint(wvnsw1(mb))
3676 iw2 = nint(wvnsw2(mb))
3677
3678 if ( wvnsw2(mb)>=wvn550 .and. wvn550>=wvnsw1(mb) ) then
3679 nv_aod = ib ! sw band number covering 550nm wavelenth
3680 endif
3681
3682! -- for rd-dependent
3683 do while ( iw1 > iendwv(ii) )
3684 if ( ii == kaerbndd ) exit
3685 ii = ii + 1
3686 enddo
3687 sumsol = f_zero
3688 nv1(ib) = ii
3689
3690! -- for rd-independent
3691 do while ( iw1 > iendwv_du(iix) )
3692 if ( iix == kaerbndi ) exit
3693 iix = iix + 1
3694 enddo
3695 sumsol_du = f_zero
3696 nv1_du(ib) = iix
3697
3698 do iw = iw1, iw2
3699! -- for rd-dependent
3700 solbnd(ib) = solbnd(ib) + solfwv(iw)
3701 sumsol = sumsol + solfwv(iw)
3702
3703 if ( iw == iendwv(ii) ) then
3704 solwaer(ib,ii) = sumsol
3705 if ( ii < kaerbndd ) then
3706 sumsol = f_zero
3707 ii = ii + 1
3708 endif
3709 endif
3710
3711! -- for rd-independent
3712 solbnd_du(ib) = solbnd_du(ib) + solfwv(iw)
3713 sumsol_du = sumsol_du + solfwv(iw)
3714
3715 if ( iw == iendwv_du(iix) ) then
3716 solwaer_du(ib,iix) = sumsol_du
3717 if ( iix < kaerbndi ) then
3718 sumsol_du = f_zero
3719 iix = iix + 1
3720 endif
3721 endif
3722 enddo
3723
3724 if ( iw2 /= iendwv(ii) ) then
3725 solwaer(ib,ii) = sumsol
3726 endif
3727 if ( iw2 /= iendwv_du(iix) ) then
3728 solwaer_du(ib,iix) = sumsol_du
3729 endif
3730
3731 nv2(ib) = ii
3732 nv2_du(ib) = iix
3733 enddo ! end do_ib_block for sw
3734 endif ! end if_laswflg_block
3735
3736! --- ... compute lw flux weights and interval indices for mapping
3737! spectral bands between lw radiation and aerosol data
3738
3739 if ( lalwflg ) then
3740 eirbnd(:) = f_zero
3741 eirbnd_du(:) = f_zero
3742 do i=1,nlwbnd
3743 do j=1,kaerbndd
3744 eirwaer(i,j) = f_zero
3745 enddo
3746 do j=1,kaerbndi
3747 eirwaer_du(i,j) = f_zero
3748 enddo
3749 enddo
3750
3751 do ib = 1, nlwbnd
3752 ii = 1
3753 iix = 1
3754 if ( nlwbnd == 1 ) then
3755 iw1 = 400 ! corresponding 25 mu
3756 iw2 = 2500 ! corresponding 4 mu
3757 else
3758 mb = ib + nlwstr - 1
3759 iw1 = nint(wvnlw1(mb))
3760 iw2 = nint(wvnlw2(mb))
3761 endif
3762
3763! -- for rd-dependent
3764 do while ( iw1 > iendwv(ii) )
3765 if ( ii == kaerbndd ) exit
3766 ii = ii + 1
3767 enddo
3768 sumir = f_zero
3769 nr1(ib) = ii
3770
3771! -- for rd-independent
3772 do while ( iw1 > iendwv_du(iix) )
3773 if ( iix == kaerbndi ) exit
3774 iix = iix + 1
3775 enddo
3776 sumir_du = f_zero
3777 nr1_du(ib) = iix
3778
3779 do iw = iw1, iw2
3780! -- for rd-dependent
3781 eirbnd(ib) = eirbnd(ib) + eirfwv(iw)
3782 sumir = sumir + eirfwv(iw)
3783
3784 if ( iw == iendwv(ii) ) then
3785 eirwaer(ib,ii) = sumir
3786
3787 if ( ii < kaerbndd ) then
3788 sumir = f_zero
3789 ii = ii + 1
3790 endif
3791 endif
3792
3793! -- for rd-independent
3794 eirbnd_du(ib) = eirbnd_du(ib) + eirfwv(iw)
3795 sumir_du = sumir_du + eirfwv(iw)
3796
3797 if ( iw == iendwv_du(iix) ) then
3798 eirwaer_du(ib,iix) = sumir_du
3799
3800 if ( iix < kaerbndi ) then
3801 sumir_du = f_zero
3802 iix = iix + 1
3803 endif
3804 endif
3805 enddo
3806
3807 if ( iw2 /= iendwv(ii) ) then
3808 eirwaer(ib,ii) = sumir
3809 endif
3810 if ( iw2 /= iendwv_du(iix) ) then
3811 eirwaer_du(ib,iix) = sumir_du
3812 endif
3813
3814 nr2(ib) = ii
3815 nr2_du(ib) = iix
3816 enddo ! end do_ib_block for lw
3817 endif ! end if_lalwflg_block
3818
3819! --- compute spectral band mean properties for each species
3820
3821 call optavg_gocart
3822! --- inputs: (in-scope variables, module variables)
3823! --- outputs: (module variables)
3824
3825
3826! --- check print
3827! if (me == 0) then
3828! do ib = 1, NSWBND
3829! mb = ib + NSWSTR - 1
3830! print *, ' wvnsw1,wvnsw2 :',wvnsw1(mb),wvnsw2(mb)
3831! print *, ' After optavg_gocart, for sw band:',ib
3832! print *, ' extrhi:', extrhi_grt(:,ib)
3833! print *, ' scarhi:', scarhi_grt(:,ib)
3834! print *, ' ssarhi:', ssarhi_grt(:,ib)
3835! print *, ' asyrhi:', asyrhi_grt(:,ib)
3836! do i = 1, KRHLEV
3837! print *, ' extrhd for rhlev:',i
3838! print *, extrhd_grt(i,:,ib)
3839! print *, ' scarhd for rhlev:',i
3840! print *, scarhd_grt(i,:,ib)
3841! print *, ' ssarhd for rhlev:',i
3842! print *, ssarhd_grt(i,:,ib)
3843! print *, ' asyrhd for rhlev:',i
3844! print *, asyrhd_grt(i,:,ib)
3845! enddo
3846! enddo
3847! print *, ' wvnlw1 :',wvnlw1
3848! print *, ' wvnlw2 :',wvnlw2
3849! do ib = 1, NLWBND
3850! ii = NSWBND + ib
3851! print *,' After optavg_gocart, for lw band:',ib
3852! print *,' extrhi_grt:', extrhi_grt(:,ii)
3853! print *,' scarhi_grt:', scarhi_grt(:,ii)
3854! print *,' ssarhi_grt:', ssarhi_grt(:,ii)
3855! print *,' asyrhi_grt:', asyrhi_grt(:,ii)
3856! do i = 1, KRHLEV
3857! print *,' extrhd for rhlev:',i
3858! print *, extrhd_grt(i,:,ib)
3859! print *,' scarhd for rhlev:',i
3860! print *, scarhd_grt(i,:,ib)
3861! print *,' ssarhd for rhlev:',i
3862! print *, ssarhd_grt(i,:,ib)
3863! print *,' asyrhd for rhlev:',i
3864! print *, asyrhd_grt(i,:,ib)
3865! enddo
3866! enddo
3867! endif
3868
3869! =================
3870 contains
3871! =================
3872
3873!-----------------------------
3876 subroutine rd_gocart_luts
3877!.............................
3878! --- inputs: (in scope variables, module variables)
3879! --- outputs: (in scope variables)
3880
3881! ==================================================================== !
3882! !
3883! subprogram: rd_gocart_luts !
3884! read GMAO pre-tabultaed aerosol optical data for dust, seasalt, !
3885! sulfate, black carbon, and organic carbon aerosols !
3886! !
3887! major local variables: !
3888! for handling spectral band structures !
3889! iendwv - ending wvnum (cm**-1) for each band kaerbndd !
3890! iendwv_du - ending wvnum (cm**-1) for each band kaerbndi !
3891! for handling optical properties of rh independent species (kcm1) !
3892! 1=du001, 2=du002, 3=du003, 4=du004, 5=du005 !
3893! rhidext0_grt - extinction coefficient kaerbndi*kcm1 !
3894! rhidsca0_grt - scattering coefficient kaerbndi*kcm1 !
3895! rhidssa0_grt - single scattering albedo kaerbndi*kcm1 !
3896! rhidasy0_grt - asymmetry parameter kaerbndi*kcm1 !
3897! for handling optical properties of rh ndependent species (kcm2) !
3898! 1=ss001, 2=ss002, 3=ss003, 4=ss004, 5=ss005, 6=so4, !
3899! 7=bcphobic, 8=bcphilic, 9=ocphobic, 10=ocphilic !
3900! rhdpext0_grt - extinction coefficient kaerbndd*krhlev*kcm2!
3901! rhdpsca0_grt - scattering coefficient kaerbndd*krhlev*kcm2!
3902! rhdpssa0_grt - single scattering albedo kaerbndd*krhlev*kcm2!
3903! rhdpasy0_grt - asymmetry parameter kaerbndd*krhlev*kcm2!
3904! !
3905! usage: call rd_gocart_luts !
3906! !
3907! ================================================================== !
3908!
3909 implicit none
3910
3911! --- inputs: (none)
3912! --- output: (none)
3913
3914! --- locals:
3915 integer :: iradius, ik, ibeg
3916 integer, parameter :: numspc = 5 ! # of aerosol species
3917
3918! - input tabulated aerosol optical spectral data from GSFC
3919 real, dimension(kaerbndd) :: lambda ! wavelength (m) for non-dust
3920 real, dimension(kaerbndi) :: lambda_du ! wavelength (m) for dust
3921 real, dimension(krhlev) :: rh ! relative humidity (fraction)
3922 real, dimension(kaerbndd,krhlev,numspc) :: bext! extinction efficiency (m2/kg)
3923 real, dimension(kaerbndd,krhlev,numspc) :: bsca! scattering efficiency (m2/kg)
3924 real, dimension(kaerbndd,krhlev,numspc) :: g ! asymmetry factor (dimensionless)
3925 real, dimension(kaerbndi,krhlev,numspc) :: bext_du! extinction efficiency (m2/kg)
3926 real, dimension(kaerbndi,krhlev,numspc) :: bsca_du! scattering efficiency (m2/kg)
3927 real, dimension(kaerbndi,krhlev,numspc) :: g_du ! asymmetry factor (dimensionless)
3928!
3929 logical :: file_exist
3930 character*50 :: fin, dummy
3931
3932! --- read LUTs for dust aerosols
3933 fin='optics_'//gridcomp(1)//'.dat'
3934 inquire (file=trim(fin), exist=file_exist)
3935 if ( file_exist ) then
3936 close(niaercm)
3937 open (unit=niaercm, file=fin, status='OLD')
3938 rewind(niaercm)
3939 else
3940 print *,' Requested luts file ',trim(fin),' not found'
3941 print *,' ** Stopped in rd_gocart_luts ** '
3942 errflg = 1
3943 errmsg = 'Requested luts file '//trim(fin)//' not found'
3944 return
3945 endif ! end if_file_exist_block
3946
3947 iradius = 5
3948! read lambda and compute mpwavelength (m)
3949 read(niaercm,'(a40)') dummy
3950 read(niaercm,*) (lambda_du(i), i=1, kaerbndi)
3951! read rh, relative humidity (fraction)
3952 read(niaercm,'(a40)') dummy
3953 read(niaercm,*) (rh(i), i=1, krhlev)
3954! read bext (m2 (kg dry mass)-1)
3955 do k = 1, iradius
3956 read(niaercm,'(a40)') dummy
3957 do j=1, krhlev
3958 read(niaercm,*) (bext_du(i,j,k), i=1,kaerbndi)
3959 enddo
3960 enddo
3961! read bsca (m2 (kg dry mass)-1)
3962 do k = 1, iradius
3963 read(niaercm,'(a40)') dummy
3964 do j=1, krhlev
3965 read(niaercm,*) (bsca_du(i,j,k), i=1, kaerbndi)
3966 enddo
3967 enddo
3968! read g (dimensionless)
3969 do k = 1, iradius
3970 read(niaercm,'(a40)') dummy
3971 do j=1, krhlev
3972 read(niaercm,*) (g_du(i,j,k), i=1, kaerbndi)
3973 enddo
3974 enddo
3975
3976! fill rhidext0 local arrays for dust aerosols (flip i-index)
3977 do i = 1, kaerbndi ! convert from m to micron
3978 j = kaerbndi -i + 1 ! flip i-index
3979 wavelength_du(j) = 1.e6 * lambda_du(i)
3980 if (int(wavelength_du(j)*100) == 55) then
3981 id550=j
3982 endif
3983 enddo
3984 do k = 1, iradius
3985 do i = 1, kaerbndi
3986 ii = kaerbndi -i + 1
3987 rhidext0_grt(ii,k) = bext_du(i,1,k)
3988 rhidsca0_grt(ii,k) = bsca_du(i,1,k)
3989 if ( bext_du(i,1,k) /= f_zero) then
3990 rhidssa0_grt(ii,k) = bsca_du(i,1,k)/bext_du(i,1,k)
3991 else
3992 rhidssa0_grt(ii,k) = f_one
3993 endif
3994 rhidasy0_grt(ii,k) = g_du(i,1,k)
3995 enddo
3996 enddo
3997
3998! --- read LUTs for non-dust aerosols
3999 do ib = 2, num_gc ! loop thru SS, SU, BC, OC
4000 fin='optics_'//gridcomp(ib)//'.dat'
4001 inquire (file=trim(fin), exist=file_exist)
4002 if ( file_exist ) then
4003 close(niaercm)
4004 open (unit=niaercm, file=fin, status='OLD')
4005 rewind(niaercm)
4006 else
4007 print *,' Requested luts file ',trim(fin),' not found'
4008 print *,' ** Stopped in rd_gocart_luts ** '
4009 errflg = 1
4010 errmsg = 'Requested luts file '//trim(fin)//' not found'
4011 return
4012 endif ! end if_file_exist_block
4013
4014 ibeg = radius_lower(ib) - kcm1
4015 iradius = num_radius(ib)
4016
4017! read lambda and compute mpwavelength (m)
4018 read(niaercm,'(a40)') dummy
4019 read(niaercm,*) (lambda(i), i=1, kaerbndd)
4020! read rh, relative humidity (fraction)
4021 read(niaercm,'(a40)') dummy
4022 read(niaercm,*) (rh(i), i=1, krhlev)
4023! read bext
4024 do k = 1, iradius
4025 read(niaercm,'(a40)') dummy
4026 do j=1, krhlev
4027 read(niaercm,*) (bext(i,j,k), i=1,kaerbndd)
4028 enddo
4029 enddo
4030! read bsca
4031 do k = 1, iradius
4032 read(niaercm,'(a40)') dummy
4033 do j=1, krhlev
4034 read(niaercm,*) (bsca(i,j,k), i=1, kaerbndd)
4035 enddo
4036 enddo
4037! read g
4038 do k = 1, iradius
4039 read(niaercm,'(a40)') dummy
4040 do j=1, krhlev
4041 read(niaercm,*) (g(i,j,k), i=1, kaerbndd)
4042 enddo
4043 enddo
4044
4045! fill rhdpext0 local arrays for non-dust aerosols (flip i-index)
4046 do i = 1, kaerbndd ! convert from m to micron
4047 j = kaerbndd -i + 1 ! flip i-index
4048 wavelength(j) = 1.e6 * lambda(i)
4049 if (int(wavelength(j)*100) == 55) then
4050 i550=j
4051 endif
4052 enddo
4053 do k = 1, iradius
4054 ik = ibeg + k - 1
4055 do i = 1, kaerbndd
4056 ii = kaerbndd -i + 1
4057 do j = 1, krhlev
4058 rhdpext0_grt(ii,j,ik) = bext(i,j,k)
4059 rhdpsca0_grt(ii,j,ik) = bsca(i,j,k)
4060 if ( bext(i,j,k) /= f_zero) then
4061 rhdpssa0_grt(ii,j,ik) = bsca(i,j,k)/bext(i,j,k)
4062 else
4063 rhdpssa0_grt(ii,j,ik) = f_one
4064 endif
4065 rhdpasy0_grt(ii,j,ik) = g(i,j,k)
4066 enddo
4067 enddo
4068 enddo
4069
4070 enddo !! ib-loop
4071
4072 return
4073!...................................
4074 end subroutine rd_gocart_luts
4075!-----------------------------------
4076
4077!--------------------------------
4082 subroutine optavg_gocart
4083!................................
4084! --- inputs: (in-scope variables, module variables)
4085! --- outputs: (module variables)
4086
4087! ==================================================================== !
4088! !
4089! subprogram: optavg_gocart !
4090! !
4091! compute mean aerosol optical properties over each sw radiation !
4092! spectral band for each of the species components. This program !
4093! follows optavg routine (in turn follows gfdl's approach for thick !
4094! cloud opertical property in sw radiation scheme (2000). !
4095! !
4096! ==================== defination of variables =================== !
4097! !
4098! major input variables: !
4099! nv1,nv2 (nswbnd) - start/end spectral band indices of aerosol data !
4100! for each sw radiation spectral band !
4101! nr1,nr2 (nlwbnd) - start/end spectral band indices of aerosol data !
4102! for each ir radiation spectral band !
4103! nv1_du,nv2_du(nswbnd) - start/end spectral band indices of aer data!
4104! for each sw radiation spectral band !
4105! nr1_du,nr2_du(nlwbnd) - start/end spectral band indices of aer data!
4106! for each ir radiation spectral band !
4107! solwaer (nswbnd,kaerbndd) !
4108! - solar flux weight over each sw radiation band !
4109! vs each aerosol data spectral band !
4110! eirwaer (nlwbnd,kaerbndd) !
4111! - ir flux weight over each lw radiation band !
4112! vs each aerosol data spectral band !
4113! solwaer_du (nswbnd,kaerbndi) !
4114! - solar flux weight over each sw radiation band !
4115! vs each aerosol data spectral band !
4116! eirwaer_du (nlwbnd,kaerbndi) !
4117! - ir flux weight over each lw radiation band !
4118! vs each aerosol data spectral band !
4119! solbnd (nswbnd) - solar flux weight over each sw radiation band !
4120! eirbnd (nlwbnd) - ir flux weight over each lw radiation band !
4121! solbnd_du(nswbnd) - solar flux weight over each sw radiation band !
4122! eirbnd_du(nlwbnd) - ir flux weight over each lw radiation band !
4123! nswbnd - total number of sw spectral bands !
4124! nlwbnd - total number of lw spectral bands !
4125! !
4126! external module variables: !
4127! laswflg - control flag for sw spectral region !
4128! lalwflg - control flag for lw spectral region !
4129! !
4130! output variables: (to module variables) !
4131! !
4132! ================================================================== !
4133
4134! --- inputs:
4135! --- output:
4136
4137! --- locals:
4138 real (kind=kind_phys) :: sumk, sums, sumok, sumokg, sumreft, &
4139 & sp, refb, reft, rsolbd, rirbd
4140
4141 integer :: ib, nb, ni, nh, nc
4142!
4143!===> ... begin here
4144!
4145! --- ... loop for each sw radiation spectral band
4146
4147 if ( laswflg ) then
4148 do nb = 1, nswbnd
4149 rsolbd = f_one / solbnd_du(nb)
4150 do nc = 1, kcm1 ! --- for rh independent aerosol species
4151 sumk = f_zero
4152 sums = f_zero
4153 sumok = f_zero
4154 sumokg = f_zero
4155 sumreft = f_zero
4156
4157 do ni = nv1_du(nb), nv2_du(nb)
4158 sp = sqrt( (f_one - rhidssa0_grt(ni,nc)) &
4159 & / (f_one - rhidssa0_grt(ni,nc)*rhidasy0_grt(ni,nc)) )
4160 reft = (f_one - sp) / (f_one + sp)
4161 sumreft = sumreft + reft*solwaer_du(nb,ni)
4162
4163 sumk = sumk + rhidext0_grt(ni,nc)*solwaer_du(nb,ni)
4164 sums = sums + rhidsca0_grt(ni,nc)*solwaer_du(nb,ni)
4165 sumok = sumok + rhidssa0_grt(ni,nc)*solwaer_du(nb,ni) &
4166 & * rhidext0_grt(ni,nc)
4167 sumokg = sumokg + rhidssa0_grt(ni,nc)*solwaer_du(nb,ni) &
4168 & * rhidext0_grt(ni,nc)*rhidasy0_grt(ni,nc)
4169 enddo
4170
4171 refb = sumreft * rsolbd
4172
4173 extrhi_grt(nc,nb) = sumk * rsolbd
4174 if (nb==nv_aod) then
4175 extrhi_grt_550(nc,1) = rhidext0_grt(id550,nc)
4176 endif
4177 scarhi_grt(nc,nb) = sums * rsolbd
4178 asyrhi_grt(nc,nb) = sumokg / (sumok + 1.0e-10)
4179 ssarhi_grt(nc,nb) = 4.0*refb &
4180 & / ( (f_one+refb)**2 - asyrhi_grt(nc,nb)*(f_one-refb)**2 )
4181 enddo ! end do_nc_block for rh-ind aeros
4182
4183 rsolbd = f_one / solbnd(nb)
4184 do nc = 1, kcm2 ! --- for rh dependent aerosol species
4185 do nh = 1, krhlev
4186 sumk = f_zero
4187 sums = f_zero
4188 sumok = f_zero
4189 sumokg = f_zero
4190 sumreft = f_zero
4191
4192 do ni = nv1(nb), nv2(nb)
4193 sp = sqrt( (f_one - rhdpssa0_grt(ni,nh,nc)) &
4194 & /(f_one-rhdpssa0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)))
4195 reft = (f_one - sp) / (f_one + sp)
4196 sumreft = sumreft + reft*solwaer(nb,ni)
4197
4198 sumk = sumk + rhdpext0_grt(ni,nh,nc)*solwaer(nb,ni)
4199 sums = sums + rhdpsca0_grt(ni,nh,nc)*solwaer(nb,ni)
4200 sumok = sumok + rhdpssa0_grt(ni,nh,nc)*solwaer(nb,ni) &
4201 & * rhdpext0_grt(ni,nh,nc)
4202 sumokg = sumokg + rhdpssa0_grt(ni,nh,nc)*solwaer(nb,ni)&
4203 & * rhdpext0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)
4204 enddo
4205
4206 refb = sumreft * rsolbd
4207
4208 extrhd_grt(nh,nc,nb) = sumk * rsolbd
4209 if (nb==nv_aod) then
4210 extrhd_grt_550(nh,nc,1) = rhdpext0_grt(i550,nh,nc)
4211 endif
4212 scarhd_grt(nh,nc,nb) = sums * rsolbd
4213 asyrhd_grt(nh,nc,nb) = sumokg / (sumok + 1.0e-10)
4214 ssarhd_grt(nh,nc,nb) = 4.0*refb &
4215 & /((f_one+refb)**2 - asyrhd_grt(nh,nc,nb)*(f_one-refb)**2)
4216
4217 enddo ! end do_nh_block
4218 enddo ! end do_nc_block for rh-dep aeros
4219
4220 enddo ! end do_nb_block for sw
4221 endif ! end if_laswflg_block
4222
4223! --- ... loop for each lw radiation spectral band
4224
4225 if ( lalwflg ) then
4226
4227 do nb = 1, nlwbnd
4228
4229 ib = nswbnd + nb
4230
4231 rirbd = f_one / eirbnd_du(nb)
4232 do nc = 1, kcm1 ! --- for rh independent aerosol species
4233 sumk = f_zero
4234 sums = f_zero
4235 sumok = f_zero
4236 sumokg = f_zero
4237 sumreft = f_zero
4238
4239 do ni = nr1_du(nb), nr2_du(nb)
4240 sp = sqrt( (f_one - rhidssa0_grt(ni,nc)) &
4241 & / (f_one - rhidssa0_grt(ni,nc)*rhidasy0_grt(ni,nc)) )
4242 reft = (f_one - sp) / (f_one + sp)
4243 sumreft = sumreft + reft*eirwaer_du(nb,ni)
4244
4245 sumk = sumk + rhidext0_grt(ni,nc)*eirwaer_du(nb,ni)
4246 sums = sums + rhidsca0_grt(ni,nc)*eirwaer_du(nb,ni)
4247 sumok = sumok + rhidssa0_grt(ni,nc)*eirwaer_du(nb,ni) &
4248 & * rhidext0_grt(ni,nc)
4249 sumokg = sumokg + rhidssa0_grt(ni,nc)*eirwaer_du(nb,ni) &
4250 & * rhidext0_grt(ni,nc)*rhidasy0_grt(ni,nc)
4251 enddo
4252
4253 refb = sumreft * rirbd
4254
4255 extrhi_grt(nc,ib) = sumk * rirbd
4256 scarhi_grt(nc,ib) = sums * rirbd
4257 asyrhi_grt(nc,ib) = sumokg / (sumok + 1.0e-10)
4258 ssarhi_grt(nc,ib) = 4.0*refb &
4259 & / ( (f_one+refb)**2 - asyrhi_grt(nc,ib)*(f_one-refb)**2 )
4260
4261 enddo ! end do_nc_block for rh-ind aeros
4262
4263 rirbd = f_one / eirbnd(nb)
4264 do nc = 1, kcm2 ! --- for rh dependent aerosol species
4265 do nh = 1, krhlev
4266 sumk = f_zero
4267 sums = f_zero
4268 sumok = f_zero
4269 sumokg = f_zero
4270 sumreft = f_zero
4271
4272 do ni = nr1(nb), nr2(nb)
4273 sp = sqrt( (f_one - rhdpssa0_grt(ni,nh,nc)) &
4274 & /(f_one-rhdpssa0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)))
4275 reft = (f_one - sp) / (f_one + sp)
4276 sumreft = sumreft + reft*eirwaer(nb,ni)
4277
4278 sumk = sumk + rhdpext0_grt(ni,nh,nc)*eirwaer(nb,ni)
4279 sums = sums + rhdpsca0_grt(ni,nh,nc)*eirwaer(nb,ni)
4280 sumok = sumok + rhdpssa0_grt(ni,nh,nc)*eirwaer(nb,ni) &
4281 & * rhdpext0_grt(ni,nh,nc)
4282 sumokg = sumokg+ rhdpssa0_grt(ni,nh,nc)*eirwaer(nb,ni) &
4283 & * rhdpext0_grt(ni,nh,nc)*rhdpasy0_grt(ni,nh,nc)
4284 enddo
4285
4286 refb = sumreft * rirbd
4287
4288 extrhd_grt(nh,nc,ib) = sumk * rirbd
4289 scarhd_grt(nh,nc,ib) = sums * rirbd
4290 asyrhd_grt(nh,nc,ib) = sumokg / (sumok + 1.0e-10)
4291 ssarhd_grt(nh,nc,ib) = 4.0*refb &
4292 & /((f_one+refb)**2 - asyrhd_grt(nh,nc,ib)*(f_one-refb)**2)
4293 enddo ! end do_nh_block
4294 enddo ! end do_nc_block for rh-dep aeros
4295
4296 enddo ! end do_nb_block for lw
4297 endif ! end if_lalwflg_block
4298!
4299 return
4300 return
4301!...................................
4302 end subroutine optavg_gocart
4303!-----------------------------------
4304
4305!...................................
4306 end subroutine gocart_aerinit
4307!-----------------------------------
4308
4336!-----------------------------------
4337 subroutine aer_property_gocart &
4338!...................................
4339
4340! --- inputs:
4341 & ( prsi,prsl,prslk,tvly,rhlay,dz,hz,tracer,aerfld, &
4342 & alon,alat,slmsk, laersw,laerlw,con_rd, &
4343 & imax,nlay,nlp1, &
4344! --- outputs:
4345 & aerosw,aerolw,aerodp,ext550,errflg,errmsg &
4346 & )
4347
4348! ================================================================== !
4349! !
4350! aer_property_gocart maps prescribed gocart aerosol data set onto !
4351! model grids, and compute aerosol optical properties for sw and !
4352! lw radiations. !
4353! !
4354! inputs: !
4355! prsi - pressure at interface mb IMAX*NLP1 !
4356! prsl - layer mean pressure (not used) IMAX*NLAY !
4357! prslk - exner function=(p/p0)**rocp (not used) IMAX*NLAY !
4358! tvly - layer virtual temperature (not used) IMAX*NLAY !
4359! rhlay - layer mean relative humidity IMAX*NLAY !
4360! dz - layer thickness m IMAX*NLAY !
4361! hz - level high m IMAX*NLP1 !
4362! tracer - aer tracer concentrations (not used) IMAX*NLAY*NTRAC!
4363! aerfld - prescribed aer tracer mixing ratios IMAX*NLAY*NTRCAER!
4364! alon, alat IMAX !
4365! - longitude and latitude of given points in degree !
4366! slmsk - sea/land mask (sea:0,land:1,sea-ice:2) IMAX !
4367! laersw,laerlw 1 !
4368! - logical flag for sw/lw aerosol calculations !
4369! IMAX - horizontal dimension of arrays 1 !
4370! NLAY,NLP1-vertical dimensions of arrays 1 !
4371! con_rd - Physical constant (gas constant for dry air) !
4372! !
4373! outputs: !
4374! aerosw - aeros opt properties for sw IMAX*NLAY*NBDSW*NF_AESW!
4375! (:,:,:,1): optical depth !
4376! (:,:,:,2): single scattering albedo !
4377! (:,:,:,3): asymmetry parameter !
4378! aerolw - aeros opt properties for lw IMAX*NLAY*NBDLW*NF_AELW!
4379! (:,:,:,1): optical depth !
4380! (:,:,:,2): single scattering albedo !
4381! (:,:,:,3): asymmetry parameter !
4382! aerodp - vertically integrated aer-opt-depth IMAX*NSPC+1 !
4383! errflg - CCPP error flag !
4384! errmsg - CCPP error message !
4385! !
4386! module parameters and constants: !
4387! NSWBND - total number of actual sw spectral bands computed !
4388! NLWBND - total number of actual lw spectral bands computed !
4389! NSWLWBD - total number of sw+lw bands computed !
4390! !
4391! module variable: (set by subroutine aer_init) !
4392! !
4393! usage: call aer_property_gocart !
4394! !
4395! ================================================================== !
4396
4397! --- inputs:
4398 integer, intent(in) :: IMAX, NLAY, NLP1
4399 logical, intent(in) :: laersw, laerlw
4400 real (kind=kind_phys), intent(in) :: con_rd
4401 real (kind=kind_phys), dimension(:,:), intent(in) :: prsi, prsl, &
4402 & prslk, tvly, rhlay, dz, hz
4403 real (kind=kind_phys), dimension(:), intent(in) :: alon, alat, &
4404 & slmsk
4405 real (kind=kind_phys), dimension(:,:,:),intent(in):: tracer
4406 real (kind=kind_phys), dimension(:,:,:),intent(in):: aerfld
4407
4408! --- outputs:
4409 real (kind=kind_phys), dimension(:,:,:,:), intent(out) :: &
4410 & aerosw, aerolw
4411 real (kind=kind_phys), dimension(:,:) , intent(out) :: aerodp
4412 real (kind=kind_phys), dimension(:,:) , intent(out) :: ext550
4413 integer, intent(out) :: errflg
4414 character(len=*), intent(out) :: errmsg
4415
4416! --- locals:
4417 real (kind=kind_phys), dimension(nlay,nswlwbd):: tauae,ssaae,asyae
4418 real (kind=kind_phys), dimension(nlay,1):: tauae_550
4419 real (kind=kind_phys), dimension(nlay,nspc) :: spcodp
4420
4421 real (kind=kind_phys),dimension(nlay,kcm) :: aerms
4422 real (kind=kind_phys),dimension(nlay) :: dz1, rh1
4423 real (kind=kind_phys) :: plv, tv, rho
4424 integer :: i, m, m1, k
4425
4426!
4427!===> ... begin here
4428!
4429
4430! Initialize CCPP error handling variables
4431 errmsg = ''
4432 errflg = 0
4433
4434 lab_do_imaxg : do i = 1, imax
4435
4436! --- initialize tauae, ssaae, asyae
4437 do m = 1, nswlwbd
4438 do k = 1, nlay
4439 tauae(k,m) = f_zero
4440 if (m==nv_aod) then
4441 tauae_550(k,1) = f_zero
4442 endif
4443 ssaae(k,m) = f_one
4444 asyae(k,m) = f_zero
4445 enddo
4446 enddo
4447
4448! --- set floor value for aerms (kg/m3)
4449 do k = 1, nlay
4450 do m = 1, kcm
4451 aerms(k,m) = 1.e-15
4452 enddo
4453 enddo
4454
4455 do k = 1, nlay
4456 do m = 1, nspc
4457 spcodp(k,m) = f_zero
4458 enddo
4459 enddo
4460
4461 do k = 1, nlay
4462 rh1(k) = rhlay(i,k) !
4463 dz1(k) = 1000.*dz(i,k) ! thickness converted from km to m
4464 plv = 100.*prsl(i,k) ! convert pressure from mb to Pa
4465 tv = tvly(i,k) ! virtual temp in K
4466 rho = plv / ( con_rd * tv) ! air density in kg/m3
4467
4468 do m = 1, kcm
4469 aerms(k,m) = aerfld(i,k,m)*rho ! dry mass (kg/m3)
4470 enddo
4471!
4472! --- calculate sw/lw aerosol optical properties for the
4473! corresponding frequency bands
4474
4475 call aeropt
4476! --- inputs: (in-scope variables)
4477! --- outputs: (in-scope variables)
4478
4479 enddo ! end_do_k_loop
4480
4481! ----------------------------------------------------------------------
4482
4483! --- update aerosw and aerolw arrays
4484 if ( laersw ) then
4485
4486 do m = 1, nbdsw
4487 do k = 1, nlay
4488 aerosw(i,k,m,1) = tauae(k,m)
4489 aerosw(i,k,m,2) = ssaae(k,m)
4490 aerosw(i,k,m,3) = asyae(k,m)
4491 enddo
4492 enddo
4493
4494! --- update diagnostic aod arrays
4495 do k = 1, nlay
4496 aerodp(i,1) = aerodp(i,1) + tauae_550(k,1)
4497 ext550(i,k) = tauae_550(k,1)
4498 do m = 1, nspc
4499 aerodp(i,m+1) = aerodp(i,m+1)+spcodp(k,m)
4500 enddo
4501 enddo
4502
4503 endif ! end if_larsw_block
4504
4505 if ( laerlw ) then
4506
4507 if ( nlwbnd == 1 ) then
4508 m1 = nswbnd + 1
4509 do m = 1, nbdlw
4510 do k = 1, nlay
4511 aerolw(i,k,m,1) = tauae(k,m1)
4512 aerolw(i,k,m,2) = ssaae(k,m1)
4513 aerolw(i,k,m,3) = asyae(k,m1)
4514 enddo
4515 enddo
4516 else
4517 do m = 1, nbdlw
4518 m1 = nswbnd + m
4519 do k = 1, nlay
4520 aerolw(i,k,m,1) = tauae(k,m1)
4521 aerolw(i,k,m,2) = ssaae(k,m1)
4522 aerolw(i,k,m,3) = asyae(k,m1)
4523 enddo
4524 enddo
4525 endif
4526
4527 endif ! end if_laerlw_block
4528
4529 enddo lab_do_imaxg
4530
4531! =================
4532 contains
4533! =================
4534
4535!--------------------------------
4538 subroutine aeropt
4539!................................
4540
4541! --- inputs: (in scope variables)
4542! --- outputs: (in scope variables)
4543
4544! ================================================================== !
4545! !
4546! compute aerosols optical properties in NSWLWBD bands for gocart !
4547! aerosol species !
4548! !
4549! input variables: !
4550! rh1 - relative humidity % NLAY !
4551! dz1 - layer thickness m NLAY !
4552! aerms - aerosol mass concentration kg/m3 NLAY*KCM !
4553! NLAY - vertical dimensions - 1 !
4554! !
4555! output variables: !
4556! tauae - optical depth - NLAY*NSWLWBD!
4557! ssaae - single scattering albedo - NLAY*NSWLWBD!
4558! asyae - asymmetry parameter - NLAY*NSWLWBD!
4559! aerodp - vertically integrated aer-opt-depth - IMAX*NSPC+1 !
4560! !
4561! ================================================================== !
4562
4563! --- inputs:
4564! --- outputs:
4565
4566! --- locals:
4567 real (kind=kind_phys) :: drh0, drh1, rdrh
4568 real (kind=kind_phys) :: cm, ext01, ext01_550, sca01,asy01,ssa01
4569 real (kind=kind_phys) :: ext1, ext1_550, asy1, ssa1,sca1,tau_550
4570 real (kind=kind_phys) :: sum_tau,sum_asy,sum_ssa,tau,asy,ssa
4571 real (kind=kind_phys) :: sum_tau_550
4572 integer :: ih1, ih2, nbin, ib, ntrc, ktrc
4573
4574! --- linear interp coeffs for rh-dep species
4575 ih2 = 1
4576 do while ( rh1(k) > rhlev_grt(ih2) )
4577 ih2 = ih2 + 1
4578 if ( ih2 > krhlev ) exit
4579 enddo
4580 ih1 = max( 1, ih2-1 )
4581 ih2 = min( krhlev, ih2 )
4582
4583 drh0 = rhlev_grt(ih2) - rhlev_grt(ih1)
4584 drh1 = rh1(k) - rhlev_grt(ih1)
4585 if ( ih1 == ih2 ) then
4586 rdrh = f_zero
4587 else
4588 rdrh = drh1 / drh0
4589 endif
4590
4591! --- compute optical properties for each spectral bands
4592 do ib = 1, nswlwbd
4593
4594 sum_tau = f_zero
4595 if (ib == nv_aod ) then
4596 sum_tau_550 = f_zero
4597 ext1_550 = f_zero
4598 endif
4599 sum_ssa = f_zero
4600 sum_asy = f_zero
4601
4602! --- determine tau, ssa, asy for dust aerosols
4603 ext1 = f_zero
4604 asy1 = f_zero
4605 sca1 = f_zero
4606 ssa1 = f_zero
4607 asy = f_zero
4608 ssa = f_zero
4609 do m = 1, kcm1
4610 cm = max(aerms(k,m),0.0) * dz1(k)
4611 ext1 = ext1 + cm*extrhi_grt(m,ib)
4612 if (ib == nv_aod) then
4613 ext1_550 = ext1_550 + cm*extrhi_grt_550(m,1)
4614 endif
4615 sca1 = sca1 + cm*scarhi_grt(m,ib)
4616 ssa1 = ssa1 + cm*extrhi_grt(m,ib) * ssarhi_grt(m,ib)
4617 asy1 = asy1 + cm*scarhi_grt(m,ib) * asyrhi_grt(m,ib)
4618 enddo ! m-loop
4619 tau = ext1
4620 if (ext1 > f_zero) ssa=min(f_one, ssa1/ext1)
4621 if (sca1 > f_zero) asy=min(f_one, asy1/sca1)
4622
4623! --- update aod from individual species
4624 if ( ib==nv_aod ) then
4625 tau_550 = ext1_550
4626 spcodp(k,1) = tau_550
4627 sum_tau_550 = sum_tau_550 + tau_550
4628 endif
4629! --- update sum_tau, sum_ssa, sum_asy
4630 sum_tau = sum_tau + tau
4631 sum_ssa = sum_ssa + tau * ssa
4632 sum_asy = sum_asy + tau * ssa * asy
4633
4634! --- determine tau, ssa, asy for non-dust aerosols
4635 do ntrc = 2, nspc
4636 ext1 = f_zero
4637 if ( ib==nv_aod ) then
4638 ext1_550 = f_zero
4639 endif
4640 asy1 = f_zero
4641 sca1 = f_zero
4642 ssa1 = f_zero
4643 ktrc = trc_to_aod(ntrc)
4644 do nbin = 1, num_radius(ntrc)
4645 m1 = radius_lower(ntrc) + nbin - 1
4646 m = m1 - num_radius(1) ! exclude dust aerosols
4647 cm = max(aerms(k,m1),0.0) * dz1(k)
4648 ext01 = extrhd_grt(ih1,m,ib) + &
4649 & rdrh * (extrhd_grt(ih2,m,ib)-extrhd_grt(ih1,m,ib))
4650 if ( ib==nv_aod ) then
4651 ext01_550 = extrhd_grt_550(ih1,m,1) + &
4652 & rdrh * (extrhd_grt_550(ih2,m,1)-extrhd_grt_550(ih1,m,1))
4653 endif
4654 sca01 = scarhd_grt(ih1,m,ib) + &
4655 & rdrh * (scarhd_grt(ih2,m,ib)-scarhd_grt(ih1,m,ib))
4656 ssa01 = ssarhd_grt(ih1,m,ib) + &
4657 & rdrh * (ssarhd_grt(ih2,m,ib)-ssarhd_grt(ih1,m,ib))
4658 asy01 = asyrhd_grt(ih1,m,ib) + &
4659 & rdrh * (asyrhd_grt(ih2,m,ib)-asyrhd_grt(ih1,m,ib))
4660 ext1 = ext1 + cm*ext01
4661 if ( ib==nv_aod ) then
4662 ext1_550 = ext1_550 + cm*ext01_550
4663 endif
4664 sca1 = sca1 + cm*sca01
4665 ssa1 = ssa1 + cm*ext01 * ssa01
4666 asy1 = asy1 + cm*sca01 * asy01
4667 enddo ! end_do_nbin_loop
4668 tau = ext1
4669 if (ext1 > f_zero) ssa=min(f_one, ssa1/ext1)
4670 if (sca1 > f_zero) asy=min(f_one, asy1/sca1)
4671! --- update aod from individual species
4672 if ( ib==nv_aod ) then
4673 tau_550 = ext1_550
4674 spcodp(k,ktrc) = tau_550
4675 sum_tau_550 = sum_tau_550 + tau_550
4676 endif
4677! --- update sum_tau, sum_ssa, sum_asy
4678 sum_tau = sum_tau + tau
4679 sum_ssa = sum_ssa + tau * ssa
4680 sum_asy = sum_asy + tau * ssa * asy
4681 enddo ! end_do_ntrc_loop
4682
4683! --- determine total tau, ssa, asy for aerosol mixture
4684 tauae(k,ib) = sum_tau
4685 if ( ib==nv_aod ) then
4686 tauae_550(k,1) = sum_tau_550
4687 endif
4688 if (sum_tau > f_zero) ssaae(k,ib) = sum_ssa / sum_tau
4689 if (sum_ssa > f_zero) asyae(k,ib) = sum_asy / sum_ssa
4690
4691 enddo ! end_do_ib_loop
4692!
4693 return
4694!................................
4695 end subroutine aeropt
4696!--------------------------------
4697
4698!...................................
4699 end subroutine aer_property_gocart
4700!-----------------------------------
4701!
4702! =======================================================================
4703!..........................................!
4704 end module module_radiation_aerosols !
4705!==========================================!
module_iounitdef
this module defines fortran unit numbers for input/output data files for the ncep gfs model.
Definition iounitdef.f:53
module_radiation_aerosols
This module contains climatological atmospheric aerosol schemes for radiation computations.
Definition radiation_aerosols.f:127
module_radlw_parameters
This module contains LW band parameters set up.
Definition radlw_param.f:61
module_radsw_parameters
This module is for specifying the band structures and program parameters used by the RRTMG-SW scheme.
Definition radsw_param.f:62