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