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