Line data Source code
1 : !! This module is a coupler between the CAM model and the Community Aerosol
2 : !! and Radiation Model for Atmospheres (CARMA) microphysics model. It adds the
3 : !! capabilities of CARMA to CAM, allowing for binned microphysics studies
4 : !! within the CAM framework. This module supports the CAM physics interface, and
5 : !! uses the CARMA and CARMASTATE objects to perform the microphysical
6 : !! calculations.
7 : !!
8 : !! @author Chuck Bardeen
9 : !! @version July 2009
10 : module carma_intr
11 :
12 : use carma_precision_mod, only: f
13 : use carma_enums_mod, only: I_OPTICS_FIXED, I_OPTICS_MIXED_CORESHELL, I_OPTICS_MIXED_VOLUME, &
14 : I_OPTICS_MIXED_MAXWELL, I_OPTICS_SULFATE, I_CNSTTYPE_PROGNOSTIC, I_HYBRID, RC_OK, RC_ERROR, &
15 : I_WTPCT_H2SO4, I_PETTERS
16 : use carma_constants_mod, only : GRAV, REARTH, WTMOL_AIR, WTMOL_H2O, R_AIR, CP, RKAPPA, NWAVE, &
17 : CARMA_NAME_LEN, CARMA_SHORT_NAME_LEN, PI, CAM_FILL, RGAS, RM2CGS, RAD2DEG, CLDFRC_INCLOUD, MAXCLDAERDIAG
18 : use carma_types_mod, only : carma_type, carmastate_type
19 : use carma_flags_mod, only : carma_flag, carma_do_fixedinit, carma_model, carma_do_wetdep, carma_do_emission, &
20 : carma_do_pheat, carma_do_substep, carma_do_thermo, carma_do_cldice, carma_diags_file, &
21 : carma_do_grow, carma_ndebugpkgs, carma_conmax, carma_cstick, carma_tstick, carma_vf_const, carma_sulfnuc_method, &
22 : carma_rhcrit, carma_rad_feedback, carma_minsubsteps, carma_maxsubsteps, carma_gstickl, carma_gsticki, &
23 : carma_maxretries, carma_dt_threshold, carma_ds_threshold, carma_do_vtran, carma_do_vdiff, carma_do_pheatatm, &
24 : carma_do_partialinit, carma_do_optics, carma_do_incloud, carma_do_explised, carma_do_drydep, carma_do_detrain, &
25 : carma_do_coremasscheck, carma_do_coag, carma_do_clearsky, carma_do_cldliq, carma_do_aerosol, carma_dgc_threshold, &
26 : carma_diags_packages, carma_ndiagpkgs
27 : use carma_model_mod, only : NGAS, NBIN, NELEM, NGROUP, NMIE_WTP, NREFIDX, MIE_RH, NMIE_RH, NSOLUTE
28 : use carma_model_mod, only : mie_rh, mie_wtp, is_convtran1, CARMAMODEL_DiagnoseBulk, CARMAMODEL_DiagnoseBins, &
29 : CARMAMODEL_Detrain, CARMAMODEL_OutputDiagnostics, CARMAMODEL_CreateOpticsFile, CARMAMODEL_WetDeposition, &
30 : CARMAMODEL_EmitParticle, CARMAMODEL_InitializeParticle, CARMAMODEL_DefineModel, CARMAMODEL_InitializeModel, &
31 : CARMAMODEL_OutputBudgetDiagnostics, CARMAMODEL_OutputCloudborneDiagnostics, CARMAMODEL_CalculateCloudborneDiagnostics
32 : use carmaelement_mod, only : CARMAELEMENT_Get
33 : use carmagas_mod, only : CARMAGAS_Get
34 : use carmagroup_mod, only : CARMAGROUP_Get
35 : use carmastate_mod, only : CARMASTATE_CreateFromReference, CARMASTATE_SetGas, CARMASTATE_Step, CARMASTATE_GetBin, &
36 : CARMASTATE_GetGas, CARMASTATE_GetState, CARMASTATE_Get, CARMASTATE_Create, CARMASTATE_SetBin, CARMASTATE_Destroy
37 : use carma_mod, only : CARMA_Get, CARMA_Create, CARMA_Initialize, CARMA_Destroy
38 :
39 : use shr_kind_mod, only: r8 => shr_kind_r8
40 : use spmd_utils, only: masterproc, mpicom
41 : use ppgrid, only: pcols, pver, pverp
42 : use ref_pres, only: pref_mid, pref_edge, pref_mid_norm, psurf_ref
43 : use physics_types, only: physics_state, physics_ptend, physics_ptend_init, &
44 : set_dry_to_wet, physics_state_copy
45 : use physconst, only: cpair
46 : use constituents, only: pcnst, cnst_add, cnst_get_ind, &
47 : cnst_name, cnst_longname
48 : use cam_abortutils, only: endrun
49 : use physics_buffer, only: physics_buffer_desc, pbuf_add_field, pbuf_old_tim_idx, &
50 : pbuf_get_index, pbuf_get_field, dtype_r8, pbuf_set_field
51 : use pio, only: var_desc_t
52 : use radconstants, only: nlwbands, nswbands
53 : use wv_sat_methods, only: wv_sat_qsat_water
54 :
55 : implicit none
56 :
57 : private
58 : save
59 :
60 : ! Public interfaces
61 :
62 : ! CAM Physics Interface
63 : public carma_register ! register consituents
64 : public carma_is_active ! retrns true if this package is active (microphysics = .true.)
65 : public carma_implements_cnst ! returns true if consituent is implemented by this package
66 : public carma_init_cnst ! initialize constituent mixing ratios, if not read from initial file
67 : public carma_init ! initialize timestep independent variables
68 : public carma_final ! finalize the CARMA module
69 : public carma_timestep_init ! initialize timestep dependent variables
70 : public carma_timestep_tend ! interface to tendency computation
71 : public carma_accumulate_stats ! collect stats from all MPI tasks
72 :
73 :
74 : ! Other Microphysics
75 : public carma_emission_tend ! calculate tendency from emission source function
76 : public carma_calculate_cloudborne_diagnostics ! calculate model specific budget diagnostics for cloudborne aerosols
77 : public carma_output_cloudborne_diagnostics ! output model specific budget diagnostics for cloudborne aerosols
78 : public carma_output_budget_diagnostics ! calculate and output model specific aerosol budget terms
79 : public carma_wetdep_tend ! calculate tendency from wet deposition
80 : public :: carma_restart_init
81 : public :: carma_restart_write
82 : public :: carma_restart_read
83 :
84 : ! Microphysics info from CAM state
85 : !
86 : ! NOTE: These calls can be used in CAM when the CAM state is available, but the CARMASTATE
87 : ! is not available. These will return the instantaneous values instead of relying on
88 : ! pbuf fields that might be from the previous timestep.
89 : public carma_get_bin
90 : public carma_get_bin_cld
91 : public carma_get_dry_radius
92 : public carma_get_elem_for_group
93 : public carma_get_group_by_name
94 : public carma_get_kappa
95 : public carma_get_number
96 : public carma_get_number_cld
97 : public carma_get_total_mmr
98 : public carma_get_total_mmr_cld
99 : public carma_get_wet_radius
100 : public carma_get_bin_rmass
101 : public carma_set_bin
102 : public carma_get_sad
103 : public :: carma_get_wght_pct
104 : public :: carma_effecitive_radius
105 :
106 : ! NOTE: This is required by physpkg.F90, since the carma_intr.F90 stub in physics/cam
107 : ! does not have access to carma_constant.F90, but needs to also provide a defintion
108 : ! for MAXCLDAERDIAG. Thus the definition of this variable needs to come from
109 : ! carma_intr.F90.
110 : public :: MAXCLDAERDIAG
111 :
112 : ! Private data
113 :
114 : ! Particle Group Statistics
115 :
116 : ! Gridbox average
117 : integer, parameter :: NGPDIAGS = 13 ! Number of particle diagnostics ...
118 : integer, parameter :: GPDIAGS_ND = 1 ! Number density
119 : integer, parameter :: GPDIAGS_AD = 2 ! Surface area density
120 : integer, parameter :: GPDIAGS_MD = 3 ! Mass density
121 : integer, parameter :: GPDIAGS_RE = 4 ! Effective Radius
122 : integer, parameter :: GPDIAGS_RM = 5 ! Mitchell [2002] Effective Radius
123 : integer, parameter :: GPDIAGS_JN = 6 ! Nucleation Rate
124 : integer, parameter :: GPDIAGS_MR = 7 ! Mass Mixing Ratio
125 : integer, parameter :: GPDIAGS_EX = 8 ! Extinction
126 : integer, parameter :: GPDIAGS_OD = 9 ! Optical Depth
127 : integer, parameter :: GPDIAGS_VM = 10 ! Mass Weighted Fall Velocity
128 : integer, parameter :: GPDIAGS_PA = 11 ! Projected Area
129 : integer, parameter :: GPDIAGS_AR = 12 ! Area Ratio
130 : integer, parameter :: GPDIAGS_VR = 13 ! Volatile Mixing Ratio
131 :
132 : ! Particle Bin (Element) Statistics
133 : integer, parameter :: NBNDIAGS = 5 ! Number of bin surface diagnostics ...
134 : integer, parameter :: BNDIAGS_TP = 1 ! Delta Particle Temperature [K]
135 : integer, parameter :: BNDIAGS_WETR = 2 ! wet radius
136 : integer, parameter :: BNDIAGS_ND = 3 ! Number density
137 : integer, parameter :: BNDIAGS_RO = 4 ! particle density
138 : integer, parameter :: BNDIAGS_VR = 5 ! Volatile Mixing Ratio
139 :
140 : ! Surface
141 : integer, parameter :: NSBDIAGS = 2 ! Number of bin surface diagnostics ...
142 : integer, parameter :: SBDIAGS_DD = 1 ! Dry deposition flux [kg/m2/s]
143 : integer, parameter :: SBDIAGS_VD = 2 ! Dry deposition velocity [cm/s]
144 :
145 :
146 : ! Gas Statistics
147 : integer, parameter :: NGSDIAGS = 5 ! Number of gas diagnostics ...
148 : integer, parameter :: GSDIAGS_SI = 1 ! saturation wrt ice
149 : integer, parameter :: GSDIAGS_SL = 2 ! saturation wrt water
150 : integer, parameter :: GSDIAGS_EI = 3 ! equilibrium vp wrt ice
151 : integer, parameter :: GSDIAGS_EL = 4 ! equilibrium vp wrt water
152 : integer, parameter :: GSDIAGS_WT = 5 ! weight percent composition for aerosols
153 :
154 : ! Step Statistics
155 : integer, parameter :: NSPDIAGS = 2 ! Number of step diagnostics ...
156 : integer, parameter :: SPDIAGS_NSTEP = 1 ! number of substeps
157 : integer, parameter :: SPDIAGS_LNSTEP = 2 ! ln(number of substeps)
158 :
159 : ! Defaults not in the namelist
160 : character(len=10), parameter :: carma_mixtype = 'wet' ! mixing ratio type for CARMA constituents
161 : integer :: LUNOPRT = 6 ! lun for output
162 :
163 : ! Constituent Mappings
164 : integer :: icnst4elem(NELEM, NBIN) ! constituent index for a carma element
165 : integer :: icnst4gas(NGAS) ! constituent index for a carma gas
166 :
167 : character(len=16) :: btndname(NGROUP, NBIN) ! names of group per bin tendencies
168 : character(len=16) :: etndname(NELEM, NBIN) ! names of element tendencies
169 : character(len=16) :: gtndname(NGAS) ! names of gas tendencies
170 :
171 : ! Flags to indicate whether each constituent could have a CARMA tendency.
172 : logical :: lq_carma(pcnst)
173 :
174 : ! The CARMA object stores the configuration inforamtion about CARMA, only one is
175 : ! is needed per MPI task. In the future, this could potentially be turned into one
176 : ! per model to allow multiple models with different numbers of bins, ... to be
177 : ! run simultaneously. However, it is more complicated than that, since some of the
178 : ! globals here would need to be put into the carma or another object and some sort
179 : ! of callback mechanism is needed to call the correct model implementations for
180 : ! the model specific methods.
181 : type(carma_type), target :: carma ! the carma object
182 :
183 :
184 : ! Physics Buffer Indicies
185 : integer :: ipbuf4gas(NGAS)=-1 ! physics buffer index for a carma gas
186 : integer :: ipbuf4t=-1 ! physics buffer index for a carma temperature
187 : integer :: ipbuf4sati(NGAS)=-1 ! physics buffer index for a carma saturation over ice
188 : integer :: ipbuf4satl(NGAS)=-1 ! physics buffer index for a carma saturation over liquid
189 :
190 : ! Globals used for a reference atmosphere.
191 : real(kind=f) :: carma_t_ref(pver) = -huge(1._f) ! midpoint temperature (Pa)
192 : real(kind=f) :: carma_h2o_ref(pver) = -huge(1._f) ! h2o mmmr (kg/kg)
193 : real(kind=f) :: carma_h2so4_ref(pver) = -huge(1._f) ! h2so4 mmr (kg/kg)
194 :
195 : type(var_desc_t) :: t_ref_desc
196 : type(var_desc_t) :: h2o_ref_desc
197 : type(var_desc_t) :: h2so4_ref_desc
198 :
199 : ! Globals used for total statistics
200 : real(kind=f) :: glob_max_nsubstep = 0._f
201 : real(kind=f) :: glob_max_nretry = 0._f
202 : real(kind=f) :: glob_nstep = 0._f
203 : real(kind=f) :: glob_nsubstep = 0._f
204 : real(kind=f) :: glob_nretry = 0._f
205 :
206 : real(kind=f) :: step_max_nsubstep = 0._f
207 : real(kind=f) :: step_max_nretry = 0._f
208 : real(kind=f) :: step_nstep = 0._f
209 : real(kind=f) :: step_nsubstep = 0._f
210 : real(kind=f) :: step_nretry = 0._f
211 :
212 : contains
213 :
214 :
215 : !! Read the names of the constituents from CARMA and automatically create
216 : !! a list of names based on the constituents and the number of size
217 : !! bins. A naming convention is used to map from CARMA constiuent & bin to
218 : !! CAM constituent, with the smallest bin being <shortname>01, then next
219 : !! shortname<02>, ...
220 : !!
221 : !! A check is done to see if the CARMA gases are already present. If so,
222 : !! they gases are linked; otherwise, the gas is added to CAM. The shortname
223 : !! of the gas is used as the constituent name.
224 : !!
225 : !! NOTE: This call is part of the CAM Physics Interface
226 : !!
227 : !! @author Chuck Bardeen
228 : !! @version May-2009
229 3072 : subroutine carma_register
230 : use radconstants, only : get_sw_spectral_boundaries, get_lw_spectral_boundaries
231 : use cam_logfile, only : iulog
232 : use cam_control_mod, only : initial_run
233 : use physconst, only: gravit, p_rearth=>rearth, mwdry, mwh2o
234 : use phys_control, only: phys_getopts
235 :
236 : implicit none
237 :
238 : integer :: ielem ! element index
239 : integer :: ibin ! bin index
240 : integer :: igas ! gas index
241 : integer :: igroup ! group index
242 : integer :: rc ! CARMA return code
243 : character(len=8) :: c_name ! constituent name
244 : character(len=50) :: c_longname ! constituent long name
245 : real(r8) :: wave(NWAVE) ! wavelength band centers (cm)
246 : real(r8) :: dwave(NWAVE) ! wavelength band width (cm)
247 : logical :: do_wave_emit(NWAVE) ! do emission in band?
248 : real(r8) :: r(NBIN) ! particle radius (cm)
249 : real(r8) :: rmass(NBIN) ! particle mass (g)
250 : character(len=8) :: shortname ! short (CAM) name
251 : character(len=8) :: grp_short ! group short (CAM) name
252 : character(len=50) :: name ! long (CARMA) name
253 : real(r8) :: wtmol ! gas molecular weight
254 : integer :: cnsttype ! constituent type
255 : integer :: maxbin ! last prognostic bin
256 : logical :: ndropmixed ! tracer is vertically mixed in ndrop
257 :
258 : character(len=16) :: radiation_scheme ! CAM's radiation package.
259 :
260 : ! Initialize the return code.
261 1536 : rc = 0
262 :
263 : ! Some constants are set on the fly in CAM, so initialize them and any derived "constants" here.
264 : ! Some of them are needed in CARMA_DefineModel and CARMA_Initialize.
265 1536 : GRAV = gravit * RM2CGS
266 1536 : REARTH = p_rearth * RM2CGS
267 1536 : WTMOL_AIR = mwdry
268 1536 : WTMOL_H2O = mwh2o
269 1536 : R_AIR = RGAS / WTMOL_AIR
270 1536 : CP = cpair * 1.e7_r8 / 1000._r8
271 1536 : RKAPPA = R_AIR / CP
272 :
273 : ! Setup the lun for output.
274 1536 : LUNOPRT = iulog
275 :
276 : ! Find out which radiation scheme is active.
277 1536 : call phys_getopts(radiation_scheme_out = radiation_scheme)
278 :
279 : ! Get the wavelength centers for the CAM longwave and shortwave bands
280 : ! from the radiation code.
281 :
282 : ! Can do this 'in place'; set wave to lower boundaries of all bands,
283 : ! dwave to upper boundaries.
284 1536 : call get_lw_spectral_boundaries( wave(:nlwbands ), dwave(:nlwbands ), 'cm')
285 1536 : call get_sw_spectral_boundaries( wave( nlwbands+1:), dwave( nlwbands+1:), 'cm')
286 :
287 : ! Now make dwave the difference and wave the average
288 47616 : dwave = dwave - wave
289 47616 : wave = wave + (dwave / 2._r8)
290 :
291 : ! NOTE: RRTMG does not calculate emission in the shortwave bands and the first and
292 : ! last shortwave bands overlap the longwave bands. At least the first and last bands
293 : ! needs to be excluded and perhaps all of the shortwave bands need to be excluded or
294 : ! double counting will happen for particle emission. The details of this should
295 : ! probably be moved into the radiation code.
296 26112 : do_wave_emit(:nlwbands) = .TRUE.
297 : ! do_wave_emit(nlwbands+1:) = .FALSE.
298 23040 : do_wave_emit(nlwbands+1:) = .TRUE.
299 1536 : do_wave_emit(nlwbands+1) = .FALSE.
300 1536 : do_wave_emit(NWAVE) = .FALSE.
301 :
302 : ! Create the CARMA object that will contain all the information about the
303 : ! how CARMA is configured.
304 : call CARMA_Create(carma, NBIN, NELEM, NGROUP, NSOLUTE, NGAS, NWAVE, rc, &
305 1536 : LUNOPRT=LUNOPRT, wave=wave, dwave=dwave, do_wave_emit=do_wave_emit, NREFIDX=NREFIDX)
306 1536 : if (rc < 0) call endrun('carma_register::CARMA_Create failed.')
307 :
308 : ! Define the microphysical model.
309 1536 : call CARMAMODEL_DefineModel(carma, rc)
310 1536 : if (rc < 0) call endrun('carma_register::CARMA_DefineModel failed.')
311 :
312 1536 : if (masterproc) then
313 2 : write(LUNOPRT,*) ''
314 2 : write(LUNOPRT,*) 'CARMA general settings for ', trim(carma_model), ' model : '
315 2 : write(LUNOPRT,*) ' carma_do_aerosol = ', carma_do_aerosol
316 2 : write(LUNOPRT,*) ' carma_do_coremasscheck = ',carma_do_coremasscheck
317 2 : write(LUNOPRT,*) ' carma_do_cldice = ', carma_do_cldice
318 2 : write(LUNOPRT,*) ' carma_do_cldliq = ', carma_do_cldliq
319 2 : write(LUNOPRT,*) ' carma_do_clearsky = ', carma_do_clearsky
320 2 : write(LUNOPRT,*) ' carma_do_coag = ', carma_do_coag
321 2 : write(LUNOPRT,*) ' carma_do_detrain = ', carma_do_detrain
322 2 : write(LUNOPRT,*) ' carma_do_drydep = ', carma_do_drydep
323 2 : write(LUNOPRT,*) ' carma_do_emission = ', carma_do_emission
324 2 : write(LUNOPRT,*) ' carma_do_fixedinit = ', carma_do_fixedinit
325 2 : write(LUNOPRT,*) ' carma_do_grow = ', carma_do_grow
326 2 : write(LUNOPRT,*) ' carma_do_explised = ', carma_do_explised
327 2 : write(LUNOPRT,*) ' carma_do_incloud = ', carma_do_incloud
328 2 : write(LUNOPRT,*) ' carma_do_partialinit= ', carma_do_partialinit
329 2 : write(LUNOPRT,*) ' carma_do_pheat = ', carma_do_pheat
330 2 : write(LUNOPRT,*) ' carma_do_pheatatm = ', carma_do_pheatatm
331 2 : write(LUNOPRT,*) ' carma_do_substep = ', carma_do_substep
332 2 : write(LUNOPRT,*) ' carma_do_thermo = ', carma_do_thermo
333 2 : write(LUNOPRT,*) ' carma_do_vdiff = ', carma_do_vdiff
334 2 : write(LUNOPRT,*) ' carma_do_vtran = ', carma_do_vtran
335 2 : write(LUNOPRT,*) ' carma_do_wetdep = ', carma_do_wetdep
336 2 : write(LUNOPRT,*) ' carma_dgc_threshold = ', carma_dgc_threshold
337 2 : write(LUNOPRT,*) ' carma_ds_threshold = ', carma_ds_threshold
338 2 : write(LUNOPRT,*) ' carma_dt_threshold = ', carma_dt_threshold
339 2 : write(LUNOPRT,*) ' carma_cstick = ', carma_cstick
340 2 : write(LUNOPRT,*) ' carma_gsticki = ', carma_gsticki
341 2 : write(LUNOPRT,*) ' carma_gstickl = ', carma_gstickl
342 2 : write(LUNOPRT,*) ' carma_tstick = ', carma_tstick
343 2 : write(LUNOPRT,*) ' carma_rhcrit = ', carma_rhcrit
344 2 : write(LUNOPRT,*) ' carma_conmax = ', carma_conmax
345 2 : write(LUNOPRT,*) ' carma_minsubsteps = ', carma_minsubsteps
346 2 : write(LUNOPRT,*) ' carma_maxsubsteps = ', carma_maxsubsteps
347 2 : write(LUNOPRT,*) ' carma_maxretries = ', carma_maxretries
348 2 : write(LUNOPRT,*) ' carma_vf_const = ', carma_vf_const
349 2 : write(LUNOPRT,*) ' cldfrc_incloud = ', CLDFRC_INCLOUD
350 2 : write(LUNOPRT,*) ' carma_rad_feedback = ', carma_rad_feedback
351 2 : write(LUNOPRT,*) ' carma_sulfnuc_method = ', carma_sulfnuc_method
352 2 : write(LUNOPRT,*) ''
353 : endif
354 :
355 : ! Intialize the model based upon the namelist configuration.
356 : !
357 : ! NOTE: When used with CAM, the latents heats (of melting and evaporation)
358 : ! need to be constant for energy balance to work. This allows them to match the
359 : ! assumptions made in the CAM energy checking and microphysics code.
360 : call CARMA_Initialize(carma, &
361 : rc, &
362 : sulfnucl_method = carma_sulfnuc_method, &
363 : do_coremasscheck = carma_do_coremasscheck, &
364 : do_clearsky = carma_do_clearsky, &
365 : do_cnst_rlh = .true., &
366 : do_coag = carma_do_coag, &
367 : do_detrain = carma_do_detrain, &
368 : do_drydep = carma_do_drydep, &
369 : do_fixedinit = carma_do_fixedinit, &
370 : do_grow = carma_do_grow, &
371 : do_explised = carma_do_explised, &
372 : do_incloud = carma_do_incloud, &
373 : do_partialinit= carma_do_partialinit, &
374 : do_pheat = carma_do_pheat, &
375 : do_pheatatm = carma_do_pheatatm, &
376 : do_print_init = masterproc, &
377 : do_substep = carma_do_substep, &
378 : do_thermo = carma_do_thermo, &
379 : do_vdiff = carma_do_vdiff, &
380 : do_vtran = carma_do_vtran, &
381 : minsubsteps = carma_minsubsteps, &
382 : maxsubsteps = carma_maxsubsteps, &
383 : maxretries = carma_maxretries, &
384 : vf_const = carma_vf_const, &
385 : conmax = carma_conmax, &
386 : cstick = carma_cstick, &
387 : dt_threshold = carma_dt_threshold, &
388 : gsticki = carma_gsticki, &
389 : gstickl = carma_gstickl, &
390 1536 : tstick = carma_tstick)
391 1536 : if (rc < 0) call endrun('carma_register::CARMA_Initialize failed.')
392 :
393 1536 : ndropmixed = carma_model(:10)=='trop_strat'
394 :
395 : ! The elements and gases from CARMA need to be added as constituents in
396 : ! CAM (if they don't already exist). For the elements, each radius bin
397 : ! needs to be its own constiuent in CAM.
398 : !
399 : ! Some rules about the constituents:
400 : ! 1) The shortname must be 8 characters or less, so the CARMA short name
401 : ! is limited to 6 characters and 2 characters for the bin number.
402 : ! 2) The molecular weight is in kg/kmol.
403 : ! 3) The specific heat at constant pressure is in J/kg/K.
404 : ! 4) The consituents are added sequentially.
405 :
406 : ! Add a CAM constituents for each bin of each element.
407 12288 : do ielem = 1, NELEM
408 :
409 10752 : call CARMAELEMENT_Get(carma, ielem, rc, igroup=igroup, shortname=shortname, name=name)
410 10752 : if (rc < 0) call endrun('carma_register::CARMAELEMENT_Get failed.')
411 :
412 10752 : call CARMAGROUP_Get(carma, igroup, rc, cnsttype=cnsttype, r=r, rmass=rmass, maxbin=maxbin, shortname=grp_short)
413 10752 : if (rc < 0) call endrun('carma_register::CARMAGROUP_Get failed.')
414 :
415 : ! For prognostic groups, all of the bins need to be represented as actual CAM
416 : ! constituents. Diagnostic groups are determined from state information that
417 : ! is already present in CAM, and thus their bins only exist in CARMA.
418 248832 : do ibin = 1, NBIN
419 215040 : write(btndname(igroup, ibin), '(A, I2.2)') trim(grp_short), ibin
420 :
421 225792 : if (cnsttype == I_CNSTTYPE_PROGNOSTIC) then
422 : ! Bins past maxbin are treated as diagnostic even if the group
423 : ! is prognostic and thus are not advected in the paerent model.
424 215040 : if (ibin <= maxbin) then
425 :
426 215040 : write(c_name, '(A, I2.2)') trim(shortname), ibin
427 215040 : write(c_longname, '(A, e11.4, A)') trim(name) // ', ', r(ibin)*1.e4_r8, ' um'
428 :
429 : ! The molecular weight seems to be used for molecular diffusion, which
430 : ! doesn't make sense for particles. The CAM solvers are unstable if the
431 : ! mass provided is large.
432 430080 : call cnst_add(c_name, WTMOL_AIR, cpair, 0._r8, icnst4elem(ielem, ibin), &
433 215040 : longname=c_longname, mixtype=carma_mixtype, is_convtran1=is_convtran1(igroup), &
434 860160 : ndropmixed=ndropmixed )
435 : end if
436 : end if
437 : end do
438 : end do
439 :
440 : ! Find the constituent for the gas or add it if not found.
441 4608 : do igas = 1, NGAS
442 :
443 3072 : call CARMAGAS_Get(carma, igas, rc, shortname=shortname, name=name, wtmol=wtmol)
444 3072 : if (rc < 0) call endrun('carma_register::CARMAGAS_Get failed.')
445 :
446 : ! Is the gas already defined?
447 3072 : call cnst_get_ind(shortname, icnst4gas(igas))
448 :
449 : ! For substepping, we need to store the last mmr values for the gas.
450 3072 : call pbuf_add_field('CG' // shortname, 'global',dtype_r8, (/pcols, pver/), ipbuf4gas(igas))
451 :
452 : ! For substepping, we need to store the last supersaturations.
453 3072 : call pbuf_add_field('CI' // shortname, 'global',dtype_r8, (/pcols, pver/), ipbuf4sati(igas))
454 10752 : call pbuf_add_field('CL' // shortname, 'global',dtype_r8, (/pcols, pver/), ipbuf4satl(igas))
455 : end do
456 :
457 :
458 : ! For substepping, we need to store the temperature.
459 1536 : call pbuf_add_field('CT', 'global',dtype_r8, (/pcols, pver/), ipbuf4t)
460 :
461 :
462 : ! Create the optical properties files needed for RRTMG radiative transfer
463 : ! calculations.
464 : !
465 : ! NOTE: This only needs to be done once at the start of the run and does not need
466 : ! to be done for restarts.
467 : !
468 : ! NOTE: We only want to do this with RRTMG. If CAM_RT is being used, then skip this.
469 1536 : if ((masterproc) .and. (initial_run) .and. (radiation_scheme == "rrtmg") .and. (carma_do_optics)) then
470 0 : call CARMA_CreateOpticsFile(carma, rc)
471 0 : if (rc < 0) call endrun('carma_register::carma_CreateOpticsFiles failed.')
472 : end if
473 :
474 1536 : return
475 : end subroutine carma_register
476 :
477 :
478 : !! Returns true if the CARMA package is active
479 : !!
480 : !! NOTE: This call is part of the CAM Physics Interface
481 : !!
482 : !! @author Chuck Bardeen
483 : !! @version May 2009
484 0 : function carma_is_active()
485 : implicit none
486 :
487 : logical :: carma_is_active
488 :
489 0 : carma_is_active = carma_flag
490 :
491 : return
492 : end function carma_is_active
493 :
494 :
495 : !! Returns true if specified constituent is implemented by CARMA
496 : !!
497 : !! NOTE: This call is part of the CAM Physics Interface
498 : !!
499 : !! @author Chuck Bardeen
500 : !! @version May 2009
501 92160 : function carma_implements_cnst(name)
502 : implicit none
503 :
504 : character(len=*), intent(in) :: name !! constituent name
505 : logical :: carma_implements_cnst ! return value
506 :
507 : integer :: igroup
508 : integer :: ielem
509 : integer :: ibin
510 : integer :: igas
511 : integer :: rc
512 :
513 : integer :: cnsttype ! constituent type
514 : integer :: maxbin ! last prognostic bin
515 :
516 92160 : rc = 0
517 :
518 92160 : carma_implements_cnst = .false.
519 :
520 : ! Check each bin to see if it this constituent.
521 138240 : do ielem = 1, NELEM
522 1935360 : do ibin = 1, NBIN
523 1889280 : call CARMAELEMENT_Get(carma, ielem, rc, igroup=igroup)
524 1889280 : if (rc < 0) call endrun('carma_implements_cnst::CARMAELEMENT_Get failed.')
525 :
526 1889280 : call CARMAGROUP_Get(carma, igroup, rc, cnsttype=cnsttype, maxbin=maxbin)
527 1889280 : if (rc < 0) call endrun('carma_implements_cnst::CARMAGROUP_Get failed.')
528 :
529 5713920 : if (cnsttype == I_CNSTTYPE_PROGNOSTIC) then
530 :
531 : ! Bins past maxbin are treated as diagnostic even if the group
532 : ! is prognostic and thus are not advected in the parent model.
533 1889280 : if (ibin <= maxbin) then
534 :
535 1889280 : if (name == cnst_name(icnst4elem(ielem, ibin))) then
536 92160 : carma_implements_cnst = .true.
537 : return
538 : end if
539 : end if
540 : end if
541 : end do
542 : end do
543 :
544 : ! Check each gas to see if it this constituent.
545 0 : do igas = 1, NGAS
546 0 : if (name == cnst_name(icnst4gas(igas))) then
547 92160 : carma_implements_cnst = .true.
548 : return
549 : end if
550 : end do
551 :
552 : return
553 : end function carma_implements_cnst
554 :
555 :
556 : !! Initialize items in CARMA that only need to be initialized once. This
557 : !! routine is called after carma_register has been called.
558 : !!
559 : !! NOTE: This call is part of the CAM Physics Interface
560 : !!
561 : !! @author Chuck Bardeen
562 : !! @version May 2009
563 1536 : subroutine carma_init(pbuf2d)
564 : use cam_history, only: addfld, add_default, horiz_only
565 : use wrap_nf
566 : use time_manager, only: is_first_step
567 : use phys_control, only: phys_getopts
568 :
569 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
570 :
571 : integer :: ielem ! element index
572 : integer :: ibin ! bin index
573 : integer :: igas ! gas index
574 : integer :: igroup ! group index
575 : integer :: icnst ! constituent index
576 : integer :: rc ! CARMA return code
577 : character(len=8) :: sname ! short (CAM) name
578 : integer :: cnsttype ! constituent type
579 : integer :: maxbin ! last prognostic bin
580 : logical :: is_cloud ! is the group a cloud?
581 : logical :: do_drydep ! is dry deposition enabled?
582 : integer :: ncore ! number of core elements in the group
583 :
584 : logical :: history_carma
585 : logical :: history_carma_srf_flx
586 : integer :: astat
587 :
588 : 1 format(a6,4x,a11,4x,a11,4x,a11)
589 : 2 format(i6,4x,3(1pe11.3,4x))
590 :
591 : ! Initialize the return code.
592 1536 : rc = 0
593 :
594 1536 : call phys_getopts(history_carma_out=history_carma, history_carma_srf_flx_out=history_carma_srf_flx)
595 :
596 : ! Set names of constituent sources and declare them as history variables; howver,
597 : ! only prognostic variables have.
598 1536 : lq_carma(:) = .false.
599 :
600 12288 : do ielem = 1, NELEM
601 227328 : do ibin = 1, NBIN
602 215040 : call CARMAELEMENT_Get(carma, ielem, rc, igroup=igroup)
603 215040 : if (rc < 0) call endrun('carma_init::CARMAELEMENT_Get failed.')
604 :
605 215040 : call CARMAGROUP_Get(carma, igroup, rc, cnsttype=cnsttype, maxbin=maxbin, do_drydep=do_drydep)
606 215040 : if (rc < 0) call endrun('carma_init::CARMAGROUP_Get failed.')
607 :
608 655872 : if (cnsttype == I_CNSTTYPE_PROGNOSTIC) then
609 :
610 : ! Bins past maxbin are treated as diagnostic even if the group
611 : ! is prognostic and thus are not advected in the parent model.
612 215040 : if (ibin <= maxbin) then
613 :
614 215040 : icnst = icnst4elem(ielem, ibin)
615 :
616 : ! Indicate that this is a constituent whose tendency could be changed by
617 : ! CARMA.
618 215040 : lq_carma(icnst) = .true.
619 :
620 215040 : etndname(ielem, ibin) = trim(cnst_name(icnst))
621 :
622 430080 : call addfld(cnst_name(icnst), (/ 'lev' /), 'A', 'kg/kg', cnst_longname(icnst))
623 215040 : if (history_carma) then
624 215040 : call add_default(cnst_name(icnst), 1, ' ')
625 : end if
626 :
627 215040 : call addfld(trim(etndname(ielem, ibin))//'TC', (/ 'lev' /), 'A', 'kg/kg/s', &
628 645120 : trim(cnst_name(icnst)) // ' tendency')
629 215040 : call addfld(trim(etndname(ielem, ibin))//'SF', horiz_only, 'A', 'kg/m2/s', &
630 430080 : trim(cnst_name(icnst)) // ' surface emission')
631 215040 : call addfld(trim(etndname(ielem, ibin))//'EM', (/ 'lev' /), 'A', 'kg/kg/s', &
632 645120 : trim(cnst_name(icnst)) // ' emission')
633 215040 : call addfld(trim(etndname(ielem, ibin))//'WD', (/ 'lev' /), 'A', 'kg/kg/s', &
634 645120 : trim(cnst_name(icnst)) // ' wet deposition')
635 215040 : call addfld(trim(etndname(ielem, ibin))//'SW', horiz_only, 'A', 'kg/m2/s', &
636 430080 : trim(cnst_name(icnst)) // ' wet deposition flux at surface')
637 :
638 215040 : if (history_carma_srf_flx) then
639 215040 : call add_default(trim(etndname(ielem, ibin))//'EM', 1, ' ')
640 215040 : call add_default(trim(etndname(ielem, ibin))//'SF', 1, ' ')
641 215040 : call add_default(trim(etndname(ielem, ibin))//'SW', 1, ' ')
642 : end if
643 :
644 215040 : if (do_drydep) then
645 215040 : call addfld(trim(etndname(ielem, ibin))//'DD', horiz_only, 'A', 'kg/m2/s ', &
646 430080 : trim(cnst_name(icnst)) // ' dry deposition')
647 215040 : if (history_carma_srf_flx) then
648 215040 : call add_default(trim(etndname(ielem, ibin))//'DD', 1, ' ')
649 : end if
650 : end if
651 :
652 215040 : if (carma_do_pheat) then
653 0 : call addfld(trim(etndname(ielem, ibin))//'TP', (/ 'lev' /), 'A', 'K ', &
654 0 : trim(cnst_name(icnst)) // ' delta particle temperature')
655 : end if
656 : end if
657 : end if
658 : end do
659 : end do
660 :
661 4608 : do igroup = 1, NGROUP
662 3072 : call CARMAGROUP_Get(carma, igroup, rc, shortname=sname, is_cloud=is_cloud, do_drydep=do_drydep, ncore=ncore)
663 3072 : if (rc < 0) call endrun('carma_init::CARMAGROUP_GetGroup failed.')
664 :
665 : ! Gridbox average
666 : !
667 : ! NOTE: Would like use flag_xf_fill for the reffective radius fields, but cam_history
668 : ! currently only supports fill values in the entire column.
669 6144 : call addfld(trim(sname)//'ND', (/ 'lev' /), 'A', '#/cm3', trim(sname) // ' number density')
670 6144 : call addfld(trim(sname)//'AD', (/ 'lev' /), 'A', 'um2/cm3', trim(sname) // ' surface area density')
671 6144 : call addfld(trim(sname)//'MD', (/ 'lev' /), 'A', 'g/cm3', trim(sname) // ' mass density')
672 6144 : call addfld(trim(sname)//'RE', (/ 'lev' /), 'A', 'um', trim(sname) // ' effective radius')
673 6144 : call addfld(trim(sname)//'RM', (/ 'lev' /), 'A', 'um', trim(sname) // ' Mitchell effective radius')
674 6144 : call addfld(trim(sname)//'JN', (/ 'lev' /), 'A', '#/cm3/s', trim(sname) // ' nucleation rate')
675 6144 : call addfld(trim(sname)//'MR', (/ 'lev' /), 'A', 'kg/kg', trim(sname) // ' mass mixing ratio')
676 6144 : call addfld(trim(sname)//'EX', (/ 'lev' /), 'A', 'km-1', trim(sname) // ' extinction')
677 6144 : call addfld(trim(sname)//'OD', (/ 'lev' /), 'A', ' ', trim(sname) // ' optical depth')
678 6144 : call addfld(trim(sname)//'PA', (/ 'lev' /), 'A', 'cm2', trim(sname) // ' projected area')
679 6144 : call addfld(trim(sname)//'AR', (/ 'lev' /), 'A', ' ', trim(sname) // ' area ratio')
680 6144 : call addfld(trim(sname)//'VM', (/ 'lev' /), 'A', 'm/s', trim(sname) // ' fall velocity')
681 6144 : call addfld(trim(sname)//'VR', (/ 'lev' /), 'A', 'kg/kg', trim(sname) // ' volatile mass mixing ratio')
682 :
683 3072 : if (history_carma) then
684 3072 : call add_default(trim(sname)//'ND', 1, ' ')
685 3072 : call add_default(trim(sname)//'AD', 1, ' ')
686 3072 : call add_default(trim(sname)//'MD', 1, ' ')
687 3072 : call add_default(trim(sname)//'RE', 1, ' ')
688 3072 : call add_default(trim(sname)//'RM', 1, ' ')
689 3072 : call add_default(trim(sname)//'MR', 1, ' ')
690 3072 : call add_default(trim(sname)//'EX', 1, ' ')
691 3072 : call add_default(trim(sname)//'OD', 1, ' ')
692 3072 : call add_default(trim(sname)//'PA', 1, ' ')
693 3072 : call add_default(trim(sname)//'AR', 1, ' ')
694 3072 : call add_default(trim(sname)//'VM', 1, ' ')
695 3072 : call add_default(trim(sname)//'VR', 1, ' ')
696 :
697 3072 : if (carma_do_grow) then
698 3072 : call add_default(trim(sname)//'JN', 1, ' ')
699 : end if
700 : end if
701 :
702 : ! Per bin stats ..
703 3072 : if (do_drydep) then
704 64512 : do ibin = 1, NBIN
705 122880 : call addfld(trim(btndname(igroup, ibin))//'VD', horiz_only, 'A', 'm/s', &
706 187392 : trim(btndname(igroup, ibin)) // ' dry deposition velocity')
707 : end do
708 : end if
709 :
710 69120 : do ibin = 1, NBIN
711 122880 : call addfld(trim(btndname(igroup, ibin))//'ND', (/ 'lev' /), 'A', '#/cm3', &
712 245760 : trim(btndname(igroup, ibin)) // ' number density')
713 61440 : call addfld(trim(btndname(igroup, ibin))//'WR', (/ 'lev' /), 'A', 'um', &
714 184320 : trim(btndname(igroup, ibin)) // ' wet radius')
715 61440 : call addfld(trim(btndname(igroup, ibin))//'RO', (/ 'lev' /), 'A', 'g/cm3', &
716 184320 : trim(btndname(igroup, ibin)) // ' wet particle density')
717 61440 : call addfld(trim(btndname(igroup, ibin))//'VR', (/ 'lev' /), 'A', 'um', &
718 184320 : trim(btndname(igroup, ibin)) // ' volatile mixing ratio')
719 :
720 :
721 64512 : if ((carma_ndebugpkgs > 0) .and. (ncore > 0)) then
722 0 : call addfld(trim(btndname(igroup, ibin))//'LCFM', horiz_only, 'A', 'kg/m2', trim(btndname(igroup, ibin)) // ' CARMA local mass fixer fail mass ')
723 0 : call addfld(trim(btndname(igroup, ibin))//'LCFP', horiz_only, 'A', 'probability', trim(btndname(igroup, ibin)) // ' CARMA mass local fail PDF')
724 0 : call addfld(trim(btndname(igroup, ibin))//'LCR', (/ 'lev' /), 'A', 'kg/kg', trim(btndname(igroup, ibin)) // ' CARMA local mass fix MMR')
725 0 : call addfld(trim(btndname(igroup, ibin))//'LCP', (/ 'lev' /), 'A', 'probability', trim(btndname(igroup, ibin)) // ' CARMA local fix PDF')
726 :
727 0 : if (carma_diags_file > 0) then
728 0 : call add_default(trim(btndname(igroup, ibin))//'LCFM', carma_diags_file, ' ')
729 0 : call add_default(trim(btndname(igroup, ibin))//'LCFP', carma_diags_file, ' ')
730 0 : call add_default(trim(btndname(igroup, ibin))//'LCR', carma_diags_file, ' ')
731 0 : call add_default(trim(btndname(igroup, ibin))//'LCP', carma_diags_file, ' ')
732 : end if
733 : end if
734 :
735 : end do
736 :
737 : end do
738 :
739 4608 : do igas = 1, NGAS
740 3072 : icnst = icnst4gas(igas)
741 :
742 : ! Indicate that this is a constituent whose tendency could be changed by
743 : ! CARMA.
744 3072 : lq_carma(icnst) = .true.
745 3072 : gtndname(igas) = trim(cnst_name(icnst)) // 'TC'
746 :
747 6144 : call addfld(gtndname(igas), (/ 'lev' /), 'A', 'kg/kg/s', trim(cnst_name(icnst)) // ' CARMA tendency')
748 :
749 : call addfld(trim(cnst_name(icnst))//'SI', (/ 'lev' /), 'A', 'ratio', &
750 6144 : trim(cnst_name(icnst)) // ' saturation wrt ice')
751 : call addfld(trim(cnst_name(icnst))//'SL', (/ 'lev' /), 'A', 'ratio', &
752 6144 : trim(cnst_name(icnst)) // ' saturation wrt liquid')
753 : call addfld(trim(cnst_name(icnst))//'EI', (/ 'lev' /), 'A', 'mol/mol', &
754 6144 : trim(cnst_name(icnst)) // ' equilibrium vmr wrt ice')
755 : call addfld(trim(cnst_name(icnst))//'EL', (/ 'lev' /), 'A', 'mol/mol', &
756 6144 : trim(cnst_name(icnst)) // ' equilibrium vmr wrt liquid')
757 : call addfld(trim(cnst_name(icnst))//'WT', (/ 'lev' /), 'A', '%', &
758 6144 : trim(cnst_name(icnst)) // ' weight percent aerosol composition')
759 :
760 4608 : if (history_carma) then
761 3072 : call add_default(trim(cnst_name(icnst))//'SI', 1, ' ')
762 3072 : call add_default(trim(cnst_name(icnst))//'SL', 1, ' ')
763 : end if
764 : end do
765 :
766 1536 : if (carma_do_thermo) then
767 0 : call addfld('CRTT', (/ 'lev' /), 'A', 'K/s', ' CARMA temperature tendency')
768 : end if
769 :
770 : ! Add fields for diagnostic fields, and make them defaults on the first tape.
771 1536 : if (carma_do_substep) then
772 3072 : call addfld('CRNSTEP', (/ 'lev' /), 'A', ' ', 'number of carma substeps')
773 3072 : call addfld('CRLNSTEP', (/ 'lev' /), 'A', ' ', 'ln(number of carma substeps)')
774 :
775 1536 : if (history_carma) then
776 1536 : call add_default('CRNSTEP', 1, ' ')
777 1536 : call add_default('CRLNSTEP', 1, ' ')
778 : end if
779 : end if
780 :
781 :
782 : ! Set up the reference atmosphere that can be used for fixed initialization. This is
783 : ! an approximate atmospheric used to define average fall velocities, coagulation
784 : ! kernels, and growth parameters.
785 1536 : if (carma_do_fixedinit) then
786 :
787 : ! NOTE: Reading the initial condtion file using the supplied routines must
788 : ! be done outside of masterproc, so does this in all threads before deciding
789 : ! if it will be used. The initial condition file is only opened on an initial run.
790 0 : if (is_first_step()) then
791 0 : call carma_getT(carma_t_ref)
792 0 : if (carma%f_igash2o /= 0) call carma_getH2O(carma_h2o_ref)
793 0 : if (carma%f_igash2So4 /= 0) call carma_getH2SO4(carma_h2so4_ref)
794 : end if
795 : end if
796 :
797 1536 : if (is_first_step()) then
798 : ! initialize physics buffer fields
799 2304 : do igas = 1, NGAS
800 1536 : call pbuf_set_field(pbuf2d, ipbuf4gas(igas), 0.0_r8)
801 1536 : call pbuf_set_field(pbuf2d, ipbuf4sati(igas), 0.0_r8)
802 2304 : call pbuf_set_field(pbuf2d, ipbuf4satl(igas), 0.0_r8)
803 : end do
804 768 : call pbuf_set_field(pbuf2d, ipbuf4t, 0.0_r8)
805 : endif
806 :
807 : ! Do a model specific initialization.
808 1536 : call CARMAMODEL_InitializeModel(carma, lq_carma, pbuf2d, rc)
809 1536 : if (rc < 0) call endrun('carma_init::CARMA_InitializeModel failed.')
810 :
811 1536 : return
812 1536 : end subroutine carma_init
813 :
814 :
815 : !! Finalize (cleanup allocations) in the CARMA model.
816 : !!
817 : !! NOTE: This call is part of the CAM Physics Interface
818 : !!
819 : !! @author Chuck Bardeen
820 : !! @version October 2009
821 3072 : subroutine carma_final
822 : implicit none
823 :
824 : integer :: rc ! CARMA return code
825 : integer :: LUNOPRT ! logical unit number for output
826 : logical :: do_print ! do print output?
827 :
828 : 2 format(' carma_final: overall substepping statistics',/,&
829 : ' max nsubstep=',1F9.0,/,' avg nsubstep=',1F9.2,/,&
830 : ' max nretry=',1F9.0,/,' avg nretry=',1F10.4)
831 :
832 : ! Initialize the return code.
833 1536 : rc = 0
834 :
835 : ! Output the end of run statistics for CARMA
836 1536 : if (carma_do_substep) then
837 1536 : if (masterproc) then
838 2 : call CARMA_Get(carma, rc, do_print=do_print, LUNOPRT=LUNOPRT)
839 2 : if (rc < 0) call endrun('carma_final::CARMA_Get failed.')
840 :
841 2 : if (glob_nstep > 0) then
842 4 : if (do_print) write(LUNOPRT,2) glob_max_nsubstep, &
843 2 : glob_nsubstep / glob_nstep, &
844 2 : glob_max_nretry, &
845 4 : glob_nretry / glob_nstep
846 : else
847 0 : if (do_print) write(LUNOPRT,2) glob_max_nsubstep, &
848 0 : 0._r8, &
849 0 : glob_max_nretry, &
850 0 : 0._r8
851 : end if
852 : end if
853 : end if
854 :
855 :
856 : ! Do a model specific initialization.
857 1536 : call CARMA_Destroy(carma, rc)
858 1536 : if (rc < 0) call endrun('carma_final::CARMA_Destroy failed.')
859 :
860 1536 : return
861 1536 : end subroutine carma_final
862 :
863 :
864 : !! Initialization that needs to be done prior to each timestep.
865 : !!
866 : !! NOTE: This call is part of the CAM Physics Interface
867 : !!
868 : !! @author Chuck Bardeen
869 : !! @version May-2009
870 16128 : subroutine carma_timestep_init
871 : implicit none
872 :
873 16128 : if (.not. carma_flag) return
874 :
875 : ! Reset the stats, so that they are per timestep values.
876 16128 : step_max_nsubstep = 0._f
877 16128 : step_max_nretry = 0._f
878 16128 : step_nstep = 0._f
879 16128 : step_nsubstep = 0._f
880 16128 : step_nretry = 0._f
881 :
882 16128 : return
883 : end subroutine carma_timestep_init
884 :
885 :
886 : !! Calculates the tendencies for all of the constituents handled by CARMA.
887 : !! To do this:
888 : !!
889 : !! - a CARMASTATE object is created
890 : !! - it is set to the current CAM state
891 : !! - a new state is determined by CARMA
892 : !! - the difference between these states is used to determine the tendencies
893 : !! - statistics arecollected and reported
894 : !!
895 : !! NOTE: This call is part of the CAM Physics Interface
896 : !!
897 : !! NOTE: Need to add code for getting/putting last fields into the physics
898 : !! buffer from substeping.
899 : !!
900 : !! @author Chuck Bardeen
901 : !! @version May-2009
902 15978240 : subroutine carma_timestep_tend(state, cam_in, cam_out, ptend, dt, pbuf, dlf, rliq, prec_str, snow_str, &
903 : prec_sed, snow_sed, ustar, obklen)
904 : use time_manager, only: get_nstep, is_first_step
905 : use camsrfexch, only: cam_in_t, cam_out_t
906 : use planck, only: planckIntensity
907 :
908 : implicit none
909 :
910 : type(physics_state), intent(in) :: state !! physics state variables
911 : type(cam_in_t), intent(in) :: cam_in !! surface inputs
912 : type(cam_out_t), intent(inout) :: cam_out !! cam output to surface models
913 : type(physics_ptend), intent(out) :: ptend !! constituent tendencies
914 : real(r8), intent(in) :: dt !! time step (s)
915 : type(physics_buffer_desc), pointer :: pbuf(:) !! physics buffer
916 : real(r8), intent(in), optional :: dlf(pcols,pver) !! Detraining cld H20 from convection (kg/kg/s)
917 : real(r8), intent(inout), optional :: rliq(pcols) !! vertical integral of liquid not yet in q(ixcldliq)
918 : real(r8), intent(out), optional :: prec_str(pcols) !! [Total] sfc flux of precip from stratiform (m/s)
919 : real(r8), intent(out), optional :: snow_str(pcols) !! [Total] sfc flux of snow from stratiform (m/s)
920 : real(r8), intent(out), optional :: prec_sed(pcols) !! total precip from cloud sedimentation (m/s)
921 : real(r8), intent(out), optional :: snow_sed(pcols) !! snow from cloud ice sedimentation (m/s)
922 : real(r8), intent(in), optional :: ustar(pcols) !! friction velocity (m/s)
923 : real(r8), intent(in), optional :: obklen(pcols) !! Obukhov length [ m ]
924 :
925 : ! Local variables
926 72960 : type(physics_state) :: state_loc ! local physics state using wet mmr
927 : type(carma_type), pointer :: carma_ptr ! the carma state object
928 72960 : type(carmastate_type) :: cstate ! the carma state object
929 : integer :: igroup ! group index
930 : integer :: ielem ! element index
931 : integer :: ibin ! bin index
932 : integer :: igas ! gas index
933 : integer :: icol ! column index
934 : integer :: icnst ! constituent index
935 : integer :: icnst_q ! H2O constituent index
936 : integer :: rc ! CARMA return code
937 : integer :: cnsttype ! constituent type
938 : integer :: maxbin ! last prognostic bin
939 : real(r8) :: spdiags(pcols, pver, NSPDIAGS) ! CARMA step diagnostic output
940 : real(r8) :: gsdiags(pcols, pver, NGAS, NGSDIAGS) ! CARMA gas diagnostic output
941 : real(r8) :: gpdiags(pcols, pver, NGROUP, NGPDIAGS) ! CARMA group diagnostic output
942 : real(r8) :: sbdiags(pcols, NBIN, NELEM, NSBDIAGS) ! CARMA surface bin diagnostic output
943 : real(r8) :: bndiags(pcols, pver, NBIN, NELEM, NBNDIAGS) ! CARMA bin diagnostic output
944 : real(r8) :: newstate(pver) ! next state for a physics state field
945 : real(r8) :: dz(pver) ! z width
946 : real(r8) :: satice(pver) ! saturation wrt ice
947 : real(r8) :: satliq(pver) ! saturation wrt liquid
948 : real(r8) :: eqice(pver) ! equil vp wrt ice
949 : real(r8) :: eqliq(pver) ! equil vp wrt liquid
950 : real(r8) :: wtpct(pver) ! weight percent aerosol composition
951 : real(r8) :: time ! the total elapsed time (s)
952 : real(r8) :: r(NBIN) ! particle radius (cm)
953 : real(r8) :: rmass(NBIN) ! particle mass (g)
954 : real(r8) :: rrat(NBIN) ! particle maximum radius ratio ()
955 : real(r8) :: arat(NBIN) ! particle area ration ()
956 : real(r8) :: nd(pver) ! number density (cm-3)
957 : real(r8) :: ad(pver) ! area density (um2/cm3)
958 : real(r8) :: md(pver) ! mass density (g cm-3)
959 : real(r8) :: mr(pver) ! mass mixing ratio (kg/kg)
960 : real(r8) :: re(pver) ! effective radius (um)
961 : real(r8) :: rm(pver) ! Mitchell effective radius (um)
962 : real(r8) :: ex(pver) ! extinction (km-1)
963 : real(r8) :: od(pver) ! optical depth
964 : real(r8) :: re2(pver) ! N(r)*r^2 (cm2)
965 : real(r8) :: re3(pver) ! N(r)*r^3 (cm3)
966 : real(r8) :: pa(pver) ! Projected Area (cm2)
967 : real(r8) :: ar(pver) ! Area Ratio
968 : real(r8) :: vm(pver) ! Massweighted fall velocity (cm2)
969 : real(r8) :: jn(pver) ! nucleation (cm-3)
970 : real(r8) :: totalmmr(pver) ! total particle mmr (kg/kg)
971 : real(r8) :: numberDensity(pver) ! number density (cm-3)
972 : real(r8) :: nucleationRate(pver) ! nucleation rate (cm-3 s-1)
973 : real(r8) :: extinctionCoefficient(pver) ! extinction coefficient (cm2)
974 : real(r8) :: r_wet(pver) ! wet radius (um)
975 : real(r8) :: rhop_wet(pver) ! wet particle density (g/cm3)
976 : real(r8) :: dd ! dry deposition (kg/m2)
977 : real(r8) :: vd ! dry deposition velocity (cm/s)
978 : real(r8) :: vf(pverp) ! fall velocity (cm/s)
979 : real(r8) :: dtpart(pver) ! delta particle temperature (K)
980 72960 : real(r8), pointer, dimension(:, :) :: t_ptr ! last temperature
981 72960 : real(r8), pointer, dimension(:, :) :: gc_ptr ! last gas mmr
982 72960 : real(r8), pointer, dimension(:, :) :: sati_ptr ! last saturation wrt ice
983 72960 : real(r8), pointer, dimension(:, :) :: satl_ptr ! last saturation wrt liquid
984 72960 : real(r8), pointer, dimension(:, :, :) :: su_ptr ! shortwave flux up (W/m2)
985 72960 : real(r8), pointer, dimension(:, :, :) :: sd_ptr ! shortwave flux down (W/m2)
986 72960 : real(r8), pointer, dimension(:, :, :) :: lu_ptr ! longwave flux up (W/m2)
987 72960 : real(r8), pointer, dimension(:, :, :) :: ld_ptr ! longwave flux down (W/m2)
988 72960 : real(r8), pointer, dimension(:,:) :: tnd_qsnow ! external tendency on snow mass (kg/kg/s)
989 72960 : real(r8), pointer, dimension(:,:) :: tnd_nsnow ! external tendency on snow number(#/kg/s)
990 72960 : real(r8), pointer, dimension(:,:) :: re_ice ! ice effective radius (m)
991 : integer :: iz
992 : real(r8) :: cldfrc(pver) ! cloud fraction [fraction]
993 : real(r8) :: rhcrit(pver) ! relative humidity for onset of liquid clouds [fraction]
994 : real(r8) :: lndram ! land aerodynamical resistance [s/m]
995 : real(r8) :: lndfv ! surface friction velocity from land [m/s]
996 : real(r8) :: ocnram ! ocean aerodynamical resistance [s/m]
997 : real(r8) :: ocnfv ! surface friction velocity from ocean [m/s]
998 : real(r8) :: iceram ! sea ice aerodynamical resistance [s/m]
999 : real(r8) :: icefv ! surface friction velocity from sea ice [m/s]
1000 : real(r8) :: radint(pver,NWAVE) ! radiative intensity (W/m2/sr/cm)
1001 : real(kind=f) :: dwave(NWAVE) ! the wavelengths widths (cm)
1002 : real(kind=f) :: wave(NWAVE) ! the center wavelengths (cm)
1003 :
1004 : integer :: max_nsubstep
1005 : real(kind=f) :: max_nretry
1006 : real(kind=f) :: nstep
1007 : integer :: nsubstep
1008 : real(kind=f) :: nretry
1009 : real(kind=f) :: zsubsteps(pver)
1010 : logical :: is_cloud ! is the group a cloud?
1011 : logical :: is_ice ! is the group ice?
1012 : logical :: grp_do_drydep ! is dry depostion enabled for group?
1013 : logical :: do_drydep ! is dry depostion enabled?
1014 : logical :: do_detrain ! do convective detrainment?
1015 : integer :: iwvl
1016 : real(r8), parameter :: zzocen = 0.0001_r8 ! Ocean aerodynamic roughness length [m]
1017 : real(r8), parameter :: zzsice = 0.0400_r8 ! Sea ice aerodynamic roughness length [m]
1018 :
1019 :
1020 : ! Initialize the return code.
1021 72960 : rc = 0
1022 :
1023 : ! Initialize the output tendency structure.
1024 72960 : call physics_ptend_init(ptend,state%psetcols,'CARMA', ls=carma_do_thermo, lq=lq_carma)
1025 :
1026 72960 : if (present(prec_sed)) prec_sed(:) = 0._f
1027 72960 : if (present(snow_sed)) snow_sed(:) = 0._f
1028 72960 : if (present(prec_str)) prec_str(:) = 0._f
1029 72960 : if (present(snow_str)) snow_str(:) = 0._f
1030 :
1031 72960 : if (.not. carma_flag) return
1032 :
1033 : ! Determine the current time in seconds.
1034 72960 : time = dt * get_nstep() - 1
1035 :
1036 : ! The CARMA interface assumes that mass mixing ratios are relative to a
1037 : ! wet atmosphere, so convert any dry mass mixing ratios to wet.
1038 72960 : call physics_state_copy(state, state_loc)
1039 72960 : call set_dry_to_wet(state_loc, convert_cnst_type='dry')
1040 :
1041 72960 : spdiags(:, :, :) = 0.0_r8
1042 72960 : gpdiags(:, :, :, :) = 0.0_r8
1043 72960 : gsdiags(:, :, :, :) = 0.0_r8
1044 72960 : sbdiags(:, :, :, :) = 0.0_r8
1045 72960 : bndiags(:, :, :, :, :) = 0.0_r8
1046 :
1047 : ! Find the constituent index for water vapor.
1048 72960 : call cnst_get_ind('Q', icnst_q)
1049 :
1050 : ! Get pointers into pbuf ...
1051 72960 : call pbuf_get_field(pbuf, ipbuf4t, t_ptr)
1052 :
1053 : ! If doing particle heating, then get pointers to the spectral flux data provided
1054 : ! by the radiation code in the physics buffer.
1055 : !
1056 : ! NOTE: RRTMG can now be done in a subset of all levels, rather than the full
1057 : ! model height. Any implications for this code have not been worked out.
1058 72960 : if (carma_do_pheat) then
1059 0 : call pbuf_get_field(pbuf, pbuf_get_index("SU"), su_ptr)
1060 0 : call pbuf_get_field(pbuf, pbuf_get_index("SD"), sd_ptr)
1061 0 : call pbuf_get_field(pbuf, pbuf_get_index("LU"), lu_ptr)
1062 0 : call pbuf_get_field(pbuf, pbuf_get_index("LD"), ld_ptr)
1063 : end if
1064 :
1065 : ! Cloud ice pbuf fields
1066 72960 : if (carma_do_cldice) then
1067 0 : call pbuf_get_field(pbuf, pbuf_get_index("TND_QSNOW"), tnd_qsnow)
1068 0 : call pbuf_get_field(pbuf, pbuf_get_index("TND_NSNOW"), tnd_nsnow)
1069 0 : call pbuf_get_field(pbuf, pbuf_get_index("RE_ICE"), re_ice)
1070 : end if
1071 :
1072 :
1073 : ! Create a CARMASTATE object which contains state information about one
1074 : ! column of the atmosphere.
1075 72960 : carma_ptr => carma
1076 :
1077 :
1078 : ! If initializing CARMASTATE from a reference state, do it before entering the main
1079 : ! loop.
1080 : !
1081 72960 : call CARMA_Get(carma, rc, do_drydep=do_drydep)
1082 72960 : if (rc < 0) call endrun('carma_timestep_tend::CARMA_Get failed.')
1083 :
1084 72960 : if (carma_do_fixedinit) then
1085 : call CARMASTATE_CreateFromReference(cstate, &
1086 : carma_ptr, &
1087 : time, &
1088 : dt, &
1089 : pver, &
1090 : I_HYBRID, &
1091 : 40._r8, &
1092 : 255._r8, &
1093 : pref_mid_norm, &
1094 : pref_edge/psurf_ref, &
1095 : pref_mid(:), &
1096 : pref_edge(:), &
1097 : carma_t_ref(:), &
1098 : rc, &
1099 : qh2o=carma_h2o_ref, &
1100 0 : qh2so4=carma_h2so4_ref)
1101 0 : if (rc < 0) call endrun('carma_timestep_tend::CARMASTATE_CreateFromReference failed.')
1102 : end if
1103 :
1104 :
1105 : ! Process each column.
1106 1123584 : do icol = 1, state_loc%ncol
1107 :
1108 1050624 : if (is_first_step()) then
1109 3926016 : t_ptr(icol,:) = state_loc%t(icol,:)
1110 : end if
1111 :
1112 : ! For particle heating, need to get the incoming radiative intensity from
1113 : ! the radiation code.
1114 : !
1115 : ! The radiation code can optionally provide the flux up and down per band in W/m2,
1116 : ! when the compute_spectral_flux namelist variable is provided to the radiation. This
1117 : ! data needs to be scaled to a radiative intensity by assuming it is isotrotropic.
1118 1050624 : radint(:,:) = 0._f
1119 :
1120 1050624 : if (carma_do_pheat) then
1121 0 : call CARMA_Get(carma, rc, dwave=dwave, wave=wave)
1122 0 : if (rc < 0) call endrun('carma_timestep_tend::CARMA_Get failed.')
1123 :
1124 : ! CARMA may run before the radiation code for the very first time step.
1125 : ! In that case, the lu, ld, su and sd values are NaN. NaN will crash
1126 : ! the model, so instead substitute an approximation that is roughly a
1127 : ! nighttime (su=sd=0) with a black body temperature of the grid point
1128 : ! temperature (lu=ld=B(T)).
1129 : !
1130 : ! NOTE: planckIntensity is in erg/cm2/s/sr/cm and lu is in W/m2,
1131 : ! so some conversion factors are needed.
1132 0 : if (is_first_step()) then
1133 0 : su_ptr(icol, :, :) = 0._r8
1134 0 : sd_ptr(icol, :, :) = 0._r8
1135 :
1136 0 : do iwvl = 1, nlwbands
1137 0 : do iz = 1, pver
1138 0 : lu_ptr(icol, iz, iwvl) = planckIntensity(wave(iwvl), state_loc%t(icol, iz)) / 1e7_f * 1e4_f * dwave(iwvl) * PI
1139 : end do
1140 0 : lu_ptr(icol, pverp, iwvl) = lu_ptr(icol, pver, iwvl)
1141 :
1142 0 : ld_ptr(icol, 2:pverp, iwvl) = lu_ptr(icol, 1:pver, iwvl)
1143 0 : ld_ptr(icol, 1, iwvl) = lu_ptr(icol, 2, iwvl)
1144 : end do
1145 : end if
1146 :
1147 :
1148 0 : do iwvl = 1, nlwbands
1149 0 : radint(:, iwvl) = (lu_ptr(icol, 2:, iwvl) + ld_ptr(icol, :pver, iwvl)) / 2._r8 / PI / dwave(iwvl)
1150 : end do
1151 :
1152 0 : do iwvl = 1, nswbands
1153 0 : radint(:, nlwbands+iwvl) = (su_ptr(icol, 2:, iwvl) + sd_ptr(icol, :pver, iwvl)) / 2._r8 / PI / dwave(nlwbands+iwvl)
1154 : end do
1155 : end if
1156 :
1157 : call CARMASTATE_Create(cstate, &
1158 : carma_ptr, &
1159 : time, &
1160 : dt, &
1161 : pver, &
1162 : I_HYBRID, &
1163 0 : state_loc%lat(icol) * RAD2DEG, &
1164 0 : state_loc%lon(icol) * RAD2DEG, &
1165 : pref_mid_norm, &
1166 : pref_edge/psurf_ref, &
1167 0 : state_loc%pmid(icol, :), &
1168 0 : state_loc%pint(icol, :), &
1169 0 : state_loc%t(icol, :), &
1170 : rc, &
1171 0 : qh2o=state_loc%q(icol, :, icnst_q), &
1172 : told=t_ptr(icol, :), &
1173 370870272 : radint=radint)
1174 1050624 : if (rc < 0) call endrun('carma_timestep_tend::CARMASTATE_Create failed.')
1175 :
1176 :
1177 : ! Store information about the CARMA particles.
1178 :
1179 : ! For prognostic groups, the mass of the particles for each bin is stored as
1180 : ! a unique constituent within CAM.
1181 8404992 : do ielem = 1, NELEM
1182 7354368 : call CARMAELEMENT_Get(carma, ielem, rc, igroup=igroup)
1183 7354368 : if (rc < 0) call endrun('carma_timestep_tend::CARMAELEMENT_Get failed.')
1184 :
1185 7354368 : call CARMAGROUP_Get(carma, igroup, rc, cnsttype=cnsttype, maxbin=maxbin)
1186 7354368 : if (rc < 0) call endrun('carma_timestep_tend::CARMAGROUP_Get failed.')
1187 :
1188 23113728 : if (cnsttype == I_CNSTTYPE_PROGNOSTIC) then
1189 :
1190 : ! For prognostic groups, set the bin from the corresponding constituent.
1191 154441728 : do ibin = 1, NBIN
1192 :
1193 : ! Bins past maxbin are treated as diagnostic even if the group
1194 : ! is prognostic and thus are not advected in the parent model.
1195 154441728 : if (ibin <= maxbin) then
1196 10443202560 : call CARMASTATE_SetBin(cstate, ielem, ibin, state_loc%q(icol, :, icnst4elem(ielem, ibin)), rc)
1197 147087360 : if (rc < 0) call endrun('carma_timestep_tend::CARMASTATE_SetBin failed.')
1198 : else
1199 0 : newstate(:) = 0._f
1200 :
1201 0 : call CARMASTATE_SetBin(cstate, ielem, ibin, newstate, rc)
1202 0 : if (rc < 0) call endrun('carma_timestep_tend::CARMASTATE_SetBin failed.')
1203 : end if
1204 : end do
1205 : end if
1206 : end do
1207 :
1208 : ! Store information about CARMA gases.
1209 3151872 : do igas = 1, NGAS
1210 2101248 : call pbuf_get_field(pbuf, ipbuf4gas(igas), gc_ptr)
1211 2101248 : call pbuf_get_field(pbuf, ipbuf4sati(igas), sati_ptr)
1212 2101248 : call pbuf_get_field(pbuf, ipbuf4satl(igas), satl_ptr)
1213 :
1214 : ! Handle the initial case where we don't have last values.
1215 2101248 : if (is_first_step()) then
1216 7852032 : gc_ptr(icol,:) = state_loc%q(icol, :, icnst4gas(igas))
1217 7852032 : sati_ptr(icol,:) = -1._f
1218 7852032 : satl_ptr(icol,:) = -1._f
1219 : end if
1220 :
1221 2101248 : call CARMASTATE_SetGas(cstate, igas, state_loc%q(icol, :, icnst4gas(igas)), rc, &
1222 151289856 : mmr_old=gc_ptr(icol,:), satice_old=sati_ptr(icol,:), satliq_old=satl_ptr(icol,:))
1223 9455616 : if (rc < 0) call endrun('carma_timestep_tend::CARMASTATE_SetGas failed.')
1224 : end do
1225 :
1226 :
1227 1050624 : call CARMAMODEL_DiagnoseBins(carma, cstate, state_loc, pbuf, icol, dt, rc, rliq=rliq, prec_str=prec_str, snow_str=snow_str)
1228 1050624 : if (rc < 0) call endrun('carma_timestep_tend::CARMA_DiagnoseBins failed.')
1229 :
1230 :
1231 : ! If the model supports detraining of condensed water from convection, then pass
1232 : ! along the condensed H2O.
1233 1050624 : call CARMA_Get(carma, rc, do_detrain=do_detrain)
1234 1050624 : if (rc < 0) call endrun('CARMA_Detrain::CARMA_Get failed.')
1235 :
1236 1050624 : if (do_detrain) then
1237 : call CARMAMODEL_Detrain(carma, cstate, cam_in, dlf, state_loc, icol, dt, rc, rliq=rliq, prec_str=prec_str, &
1238 0 : snow_str=snow_str, tnd_qsnow=tnd_qsnow, tnd_nsnow=tnd_nsnow)
1239 0 : if (rc < 0) call endrun('carma_timestep_tend::CARMA_Detrain failed.')
1240 : end if
1241 :
1242 :
1243 : ! Now that detrainment has happened, determine the cloud fractions.
1244 : ! These will be used to scale the cloud amount to go from gridbox average to in-cloud
1245 : ! values and back.
1246 : !
1247 : ! For the cirrus model, assume the cloud fraction is just the ice cloud fraction.
1248 1050624 : call CARMA_CloudFraction(carma, cstate, cam_in, state_loc, icol, cldfrc, rhcrit, rc)
1249 1050624 : if (rc < 0) call endrun('carma_timestep_tend::carma_CloudFraction failed.')
1250 :
1251 : ! A fixed value for rhcrit can be specified in the namelist rather than using the
1252 : ! one from the cloud fraction.
1253 1050624 : if (carma_rhcrit /= 0._f) then
1254 74594304 : rhcrit(:) = carma_rhcrit
1255 : end if
1256 :
1257 :
1258 : ! For dry deposition, provide a surface friction velocity and an aerodynamic
1259 : ! resistance for each of the land surface types. The values for the land come
1260 : ! from the land model, but those for ocean and sea ice need to be calculated.
1261 1050624 : if (do_drydep) then
1262 :
1263 : ! Land
1264 1050624 : lndfv = cam_in%fv(icol)
1265 1050624 : lndram = cam_in%ram1(icol)
1266 :
1267 : ! Ocean
1268 1050624 : ocnfv = ustar(icol)
1269 1050624 : ocnram = 0._r8
1270 1050624 : if (cam_in%ocnfrac(icol) > 0._r8) then
1271 : call CARMA_calcram(ocnfv, &
1272 : zzocen, &
1273 0 : state_loc%pdel(icol, pver), &
1274 0 : state_loc%pmid(icol, pver), &
1275 0 : state_loc%t(icol, pver), &
1276 : obklen(icol), &
1277 704075 : ocnram)
1278 : end if
1279 :
1280 : ! Sea Ice
1281 1050624 : icefv = ustar(icol)
1282 1050624 : iceram = 0._r8
1283 1050624 : if (cam_in%icefrac(icol) > 0._r8) then
1284 : call CARMA_calcram(ocnfv, &
1285 : zzocen, &
1286 0 : state_loc%pdel(icol, pver), &
1287 0 : state_loc%pmid(icol, pver), &
1288 0 : state_loc%t(icol, pver), &
1289 : obklen(icol), &
1290 148388 : iceram)
1291 : end if
1292 : end if
1293 :
1294 :
1295 : ! Advance the microphysics one timestep.
1296 : call CARMASTATE_Step(cstate, rc, cldfrc=cldfrc, rhcrit=rhcrit, &
1297 : lndfv=lndfv, ocnfv=ocnfv, icefv=icefv, lndram=lndram, &
1298 1050624 : ocnram=ocnram, iceram=iceram, lndfrac=cam_in%landfrac(icol), &
1299 1050624 : ocnfrac=cam_in%ocnfrac(icol), icefrac=cam_in%icefrac(icol))
1300 1050624 : if (rc < 0) call endrun('carma_timestep_tend::CARMA_Step failed.')
1301 :
1302 :
1303 : ! Get the results for the CARMA particles.
1304 :
1305 : ! For diagnostic groups, a special routine needs to be called to determine how
1306 : ! bins affect the bulk state, since there is not an individual constituent for
1307 : ! each bin.
1308 : !
1309 : ! NOTE: To work around an XL Fortran compiler bug, the optional arguments can only
1310 : ! be passed when defined.
1311 1050624 : if (present(rliq)) then
1312 : call CARMAMODEL_DiagnoseBulk(carma, cstate, cam_out, state_loc, pbuf, ptend, icol, dt, rc, &
1313 : rliq=rliq, prec_str=prec_str, snow_str=snow_str, prec_sed=prec_sed, &
1314 0 : snow_sed=snow_sed, tnd_qsnow=tnd_qsnow, tnd_nsnow=tnd_nsnow, re_ice=re_ice)
1315 : else
1316 1050624 : call CARMAMODEL_DiagnoseBulk(carma, cstate, cam_out, state_loc, pbuf, ptend, icol, dt, rc)
1317 : end if
1318 1050624 : if (rc < 0) call endrun('carma_timestep_tend::CARMASTATE_DiagnoseBulk failed.')
1319 :
1320 :
1321 : ! Calculate the group statistics for all elements.
1322 74594304 : dz(:) = state_loc%zi(icol, 1:pver) - state_loc%zi(icol, 2:pverp)
1323 :
1324 8404992 : do ielem = 1, NELEM
1325 :
1326 7354368 : call CARMAELEMENT_Get(carma, ielem, rc, igroup=igroup)
1327 7354368 : if (rc < 0) call endrun('carma_timestep_tend::CARMAELEMENT_Get failed.')
1328 :
1329 : call CARMAGROUP_Get(carma, igroup, rc, cnsttype=cnsttype, r=r, rmass=rmass, maxbin=maxbin, &
1330 7354368 : is_cloud=is_cloud, is_ice=is_ice, do_drydep=grp_do_drydep, rrat=rrat, arat=arat)
1331 7354368 : if (rc < 0) call endrun('carma_timestep_tend::CARMAGROUP_Get failed.')
1332 :
1333 : ! Intialize the group totals
1334 7354368 : nd(:) = 0.0_r8
1335 7354368 : ad(:) = 0.0_r8
1336 7354368 : md(:) = 0.0_r8
1337 7354368 : mr(:) = 0.0_r8
1338 7354368 : re(:) = 0.0_r8
1339 7354368 : rm(:) = 0.0_r8
1340 7354368 : ex(:) = 0.0_r8
1341 7354368 : od(:) = 0.0_r8
1342 7354368 : re2(:) = 0.0_r8
1343 7354368 : re3(:) = 0.0_r8
1344 7354368 : jn(:) = 0.0_r8
1345 7354368 : pa(:) = 0.0_r8
1346 7354368 : vm(:) = 0.0_r8
1347 7354368 : ar(:) = 0.0_r8
1348 :
1349 154441728 : do ibin = 1, NBIN
1350 : call CARMASTATE_GetBin(cstate, ielem, ibin, newstate(:), rc, &
1351 : numberDensity=numberDensity, nucleationRate=nucleationRate, r_wet=r_wet, &
1352 147087360 : rhop_wet=rhop_wet, sedimentationflux=dd, vd=vd, vf=vf, dtpart=dtpart, totalmmr=totalmmr)
1353 147087360 : if (rc < 0) call endrun('carma_timestep_tend::CARMASTATE_GetBin failed.')
1354 :
1355 : ! For prognostic groups, set the tendency from the corresponding constituents.
1356 147087360 : if (cnsttype == I_CNSTTYPE_PROGNOSTIC) then
1357 :
1358 : ! Bins past maxbin are treated as diagnostic even if the group
1359 : ! is prognostic and thus are not advected in the paerent model.
1360 147087360 : if (ibin <= maxbin) then
1361 :
1362 147087360 : icnst = icnst4elem(ielem, ibin)
1363 :
1364 : ! Update the consituent tendency.
1365 10443202560 : ptend%q(icol, :, icnst) = (newstate(:) - state_loc%q(icol, :, icnst)) / dt
1366 :
1367 147087360 : if (grp_do_drydep) then
1368 147087360 : sbdiags(icol, ibin, ielem, SBDIAGS_DD) = dd
1369 147087360 : sbdiags(icol, ibin, ielem, SBDIAGS_VD) = - vd / 100._r8
1370 : end if
1371 : end if
1372 : end if
1373 :
1374 : ! Calculate the total densities.
1375 : !
1376 : ! NOTE: Convert AD to um2/cm3.
1377 147087360 : if (numberDensity(1) /= CAM_FILL) then
1378 2983772160 : nd(:) = nd(:) + numberDensity(:)
1379 2983772160 : re2(:) = re2(:) + numberDensity(:) * ((r(ibin)*rrat(ibin))**2)
1380 2983772160 : re3(:) = re3(:) + numberDensity(:) * ((r(ibin)*rrat(ibin))**3)
1381 2983772160 : ad(:) = ad(:) + numberDensity(:) * 4.0_r8 * PI * (r(ibin)**2) * 1.0e8_r8
1382 2983772160 : md(:) = md(:) + numberDensity(:) * rmass(ibin)
1383 2983772160 : mr(:) = mr(:) + totalmmr(:)
1384 2983772160 : pa(:) = pa(:) + numberDensity(:) * PI * ((r(ibin) * rrat(ibin))**2) * arat(ibin)
1385 2983772160 : vm(:) = vm(:) + numberDensity(:) * rmass(ibin) * vf(2:) / 100._f
1386 :
1387 : ! Calculate the optical depth and extinction.
1388 : !
1389 : ! NOTE: Assume Qext = 2 for optical depth. This can be pulled out of CARMA
1390 : ! mie claculations later.
1391 : !
1392 : ! Convert extinction coefficient to km-1.
1393 2983772160 : extinctionCoefficient(:) = 2.0_r8 * PI * (r(ibin)**2)
1394 2983772160 : ex(:) = ex(:) + numberDensity(:) * extinctionCoefficient(:) * 1e5_r8
1395 2983772160 : od(:) = od(:) + numberDensity(:) * extinctionCoefficient(:) * dz(:) * 100._r8
1396 : end if
1397 :
1398 10443202560 : bndiags(icol,:,ibin,ielem,BNDIAGS_VR) = bndiags(icol,:,ibin,ielem,BNDIAGS_VR) + totalmmr(:)
1399 10443202560 : gpdiags(icol, :, igroup, GPDIAGS_VR) = gpdiags(icol, :, igroup, GPDIAGS_VR) + totalmmr(:)
1400 :
1401 : ! Particle temperatures from particle heating.
1402 147087360 : if (carma_do_pheat) then
1403 0 : bndiags(icol, :, ibin, ielem, BNDIAGS_TP) = dtpart(:)
1404 : end if
1405 :
1406 147087360 : if (nucleationRate(1) /= CAM_FILL) then
1407 2983772160 : jn(:) = jn(:) + nucleationRate(:)
1408 : end if
1409 :
1410 : ! Output nd and wet radius for each bin.
1411 10443202560 : r_wet = r_wet * 1e4_r8 ! cm to um
1412 10443202560 : bndiags(icol,:,ibin,ielem,BNDIAGS_WETR) = r_wet(:)
1413 10443202560 : bndiags(icol,:,ibin,ielem,BNDIAGS_ND) = numberDensity(:)
1414 10597644288 : bndiags(icol,:,ibin,ielem,BNDIAGS_RO) = rhop_wet(:)
1415 : end do
1416 :
1417 : ! If this is the number element for the group, then write out the
1418 : ! statistics.
1419 23113728 : if (numberDensity(1) /= CAM_FILL) then
1420 :
1421 : ! Calculate the effective radius (total volume / total area). Places
1422 : ! with no surface area will cause NaN values.
1423 : !
1424 : ! NOTE: Convert RE to um.
1425 590450688 : where (re2(:) > 0.0_r8)
1426 : re(:) = (re3(:) / re2(:)) * 1e4_r8
1427 : rm(:) = (3._r8 / 4._r8) * (md(:) / (0.917_r8 * pa(:))) * 1e4_r8
1428 : ar(:) = pa(:) / PI / re2(:)
1429 : end where
1430 :
1431 149188608 : where (md(:) > 0.0_r8)
1432 : vm(:) = vm(:) / md(:)
1433 : end where
1434 :
1435 : ! Store the statistics.
1436 :
1437 : ! Gridbox average
1438 149188608 : gpdiags(icol, :, igroup, GPDIAGS_ND) = nd
1439 149188608 : gpdiags(icol, :, igroup, GPDIAGS_AD) = ad
1440 149188608 : gpdiags(icol, :, igroup, GPDIAGS_MD) = md
1441 149188608 : gpdiags(icol, :, igroup, GPDIAGS_RE) = re
1442 149188608 : gpdiags(icol, :, igroup, GPDIAGS_RM) = rm
1443 149188608 : gpdiags(icol, :, igroup, GPDIAGS_MR) = mr
1444 149188608 : gpdiags(icol, :, igroup, GPDIAGS_EX) = ex
1445 149188608 : gpdiags(icol, :, igroup, GPDIAGS_OD) = od
1446 149188608 : gpdiags(icol, :, igroup, GPDIAGS_VM) = vm
1447 149188608 : gpdiags(icol, :, igroup, GPDIAGS_PA) = pa
1448 149188608 : gpdiags(icol, :, igroup, GPDIAGS_AR) = ar
1449 :
1450 2101248 : if (nucleationRate(1) /= CAM_FILL) then
1451 149188608 : gpdiags(icol, :, igroup, GPDIAGS_JN) = jn
1452 : end if
1453 : end if
1454 : end do
1455 :
1456 :
1457 : ! Get the results for the CARMA gases.
1458 3151872 : do igas = 1, NGAS
1459 2101248 : call pbuf_get_field(pbuf, ipbuf4gas(igas), gc_ptr)
1460 2101248 : call pbuf_get_field(pbuf, ipbuf4sati(igas), sati_ptr)
1461 2101248 : call pbuf_get_field(pbuf, ipbuf4satl(igas), satl_ptr)
1462 :
1463 : call CARMASTATE_GetGas(cstate, igas, newstate(:), rc, satice=satice, satliq=satliq, &
1464 2101248 : eqice=eqice, eqliq=eqliq, wtpct=wtpct)
1465 2101248 : if (rc < 0) call endrun('carma_timestep_tend::CARMASTATE_GetGas failed.')
1466 :
1467 2101248 : icnst = icnst4gas(igas)
1468 :
1469 149188608 : ptend%q(icol, :, icnst) = (newstate(:) - state_loc%q(icol, :, icnst)) / dt
1470 :
1471 149188608 : gsdiags(icol, :, igas, GSDIAGS_SI) = satice(:)
1472 149188608 : gsdiags(icol, :, igas, GSDIAGS_SL) = satliq(:)
1473 149188608 : gsdiags(icol, :, igas, GSDIAGS_EI) = eqice(:)
1474 149188608 : gsdiags(icol, :, igas, GSDIAGS_EL) = eqliq(:)
1475 149188608 : gsdiags(icol, :, igas, GSDIAGS_WT) = wtpct(:)
1476 :
1477 : ! Store the values needed for substepping in the physics buffer.
1478 149188608 : gc_ptr(icol,:) = newstate(:)
1479 149188608 : sati_ptr(icol, :) = satice(:)
1480 156542976 : satl_ptr(icol, :) = satliq(:)
1481 : end do
1482 :
1483 :
1484 : ! Get the results for temperature.
1485 1050624 : call CARMASTATE_GetState(cstate, rc, t=newstate(:))
1486 1050624 : if (rc < 0) call endrun('carma_timestep_tend::CARMASTATE_GetState failed.')
1487 :
1488 : ! Store the values needed for substepping in the physics buffer.
1489 74594304 : t_ptr(icol,:) = newstate(:)
1490 :
1491 1050624 : if (carma_do_thermo) then
1492 0 : ptend%s(icol, :) = (newstate(:) - state_loc%t(icol, :)) * cpair / dt
1493 : endif
1494 :
1495 :
1496 : ! Get the substepping statistics
1497 5326080 : if (carma_do_substep) then
1498 1050624 : call CARMASTATE_Get(cstate, rc, zsubsteps=zsubsteps)
1499 1050624 : if (rc < 0) call endrun('carma_timestep_tend::CARMASTATE_Get failed.')
1500 :
1501 74594304 : spdiags(icol, :, SPDIAGS_NSTEP) = zsubsteps(:)
1502 74594304 : where (zsubsteps/=0.0_r8)
1503 : spdiags(icol, :, SPDIAGS_LNSTEP) = log(zsubsteps(:))
1504 : end where
1505 : end if
1506 : end do
1507 :
1508 :
1509 : ! Report substep diagnostics
1510 72960 : if (carma_do_substep) then
1511 : call CARMASTATE_Get(cstate, rc, max_nsubstep=max_nsubstep, max_nretry=max_nretry, &
1512 72960 : nstep=nstep, nsubstep=nsubstep, nretry=nretry)
1513 72960 : if (rc < 0) call endrun('carma_timestep_tend::CARMASTATE_Get failed.')
1514 :
1515 : !$OMP CRITICAL
1516 72960 : step_max_nsubstep = max(step_max_nsubstep, real(max_nsubstep, f))
1517 72960 : step_max_nretry = max(step_max_nretry, max_nretry)
1518 :
1519 72960 : step_nstep = step_nstep + nstep
1520 72960 : step_nsubstep = step_nsubstep + real(nsubstep, f)
1521 72960 : step_nretry = step_nretry + nretry
1522 : !$OMP END CRITICAL
1523 : end if
1524 :
1525 : ! The CARMASTATE object is no longer needed.
1526 72960 : call CARMASTATE_Destroy(cstate, rc)
1527 72960 : if (rc < 0) call endrun('carma_timestep_tend::CARMASTATE_Destroy failed.')
1528 :
1529 : ! Output diagnostic fields.
1530 72960 : call carma_output_diagnostics(state_loc, ptend, pbuf, cam_in, gpdiags, sbdiags, gsdiags, spdiags, bndiags)
1531 :
1532 145920 : end subroutine carma_timestep_tend
1533 :
1534 :
1535 : !! Get the index for the constituents array for the specified bin
1536 : !! of the specified element.
1537 : !!
1538 : !! @author Yunqian Zhu, Francis Vitt
1539 : !! @version September-2022
1540 : subroutine carma_getcnstforbin(ielem, ibin, icnst)
1541 : implicit none
1542 :
1543 : integer, intent(in) :: ielem, ibin
1544 : integer, intent(out) :: icnst
1545 :
1546 : icnst = icnst4elem(ielem,ibin)
1547 : return
1548 72960 : end subroutine carma_getcnstforbin
1549 :
1550 : !! Collect CARMA substep statistics from all MPI tasks.
1551 : !!
1552 : !! @author Chuck Bardeen
1553 : !! @version May-2009
1554 14592 : subroutine carma_accumulate_stats()
1555 : #if ( defined SPMD )
1556 : use mpishorthand
1557 : #endif
1558 : implicit none
1559 :
1560 : integer :: istat
1561 : integer :: rc
1562 : real(kind=f) :: wrk
1563 : integer :: LUNOPRT ! logical unit number for output
1564 : logical :: do_print ! do print output?
1565 :
1566 : ! Define formats
1567 : 1 format(' carma: max nsubstep=',1F9.0,3x,'avg nsubstep=',1F9.2,3x,'max nretry=',1F9.0,3x,'avg nretry=',1F10.4)
1568 :
1569 14592 : if (carma_do_substep) then
1570 :
1571 14592 : call CARMA_Get(carma, rc, do_print=do_print, LUNOPRT=LUNOPRT)
1572 14592 : if (rc < 0) call endrun('carma_init::CARMA_Get failed.')
1573 :
1574 : #ifdef SPMD
1575 14592 : call mpi_allreduce(step_max_nsubstep, wrk, 1, mpir8, mpi_max, mpicom, istat)
1576 14592 : if( istat /= MPI_SUCCESS ) then
1577 0 : if (do_print) write(LUNOPRT,*) 'carma_timestep_tend: MPI_ALLREDUCE for max_nsubstep failed; error = ',istat
1578 0 : call endrun
1579 : end if
1580 14592 : step_max_nsubstep = wrk
1581 14592 : glob_max_nsubstep = max(glob_max_nsubstep, wrk)
1582 :
1583 14592 : call mpi_allreduce(step_max_nretry, wrk, 1, mpir8, mpi_max, mpicom, istat)
1584 14592 : if( istat /= MPI_SUCCESS ) then
1585 0 : if (do_print) write(LUNOPRT,*) 'carma_timestep_tend: MPI_ALLREDUCE for max_nsubstep failed; error = ',istat
1586 0 : call endrun
1587 : end if
1588 14592 : step_max_nretry = wrk
1589 14592 : glob_max_nretry = max(glob_max_nretry, wrk)
1590 :
1591 14592 : call mpi_allreduce(step_nstep, wrk, 1, mpir8, mpi_sum, mpicom, istat)
1592 14592 : if( istat /= MPI_SUCCESS ) then
1593 0 : if (do_print) write(LUNOPRT,*) 'carma_timestep_tend: MPI_ALLREDUCE for nstep failed; error = ',istat
1594 0 : call endrun
1595 : end if
1596 14592 : step_nstep = wrk
1597 14592 : glob_nstep = glob_nstep + wrk
1598 :
1599 14592 : call mpi_allreduce(step_nsubstep, wrk, 1, mpir8, mpi_sum, mpicom, istat)
1600 14592 : if( istat /= MPI_SUCCESS ) then
1601 0 : if (do_print) write(LUNOPRT,*) 'carma_timestep_tend: MPI_ALLREDUCE for nsubstep failed; error = ',istat
1602 0 : call endrun
1603 : end if
1604 14592 : step_nsubstep = wrk
1605 14592 : glob_nsubstep = glob_nsubstep + wrk
1606 :
1607 14592 : call mpi_allreduce(step_nretry, wrk, 1, mpir8, mpi_sum, mpicom, istat)
1608 14592 : if( istat /= MPI_SUCCESS ) then
1609 0 : if (do_print) write(LUNOPRT,*) 'carma_timestep_tend: MPI_ALLREDUCE for nretry failed; error = ',istat
1610 0 : call endrun
1611 : end if
1612 14592 : step_nretry = wrk
1613 14592 : glob_nretry = glob_nretry + wrk
1614 : #else
1615 :
1616 : ! For single CPU or OMP, just set the globals directly.
1617 : glob_max_nsubstep = max(glob_max_nsubstep, step_max_nsubstep)
1618 : glob_max_nretry = max(glob_max_nretry, step_max_nretry)
1619 : glob_nstep = glob_nstep + step_nstep
1620 : glob_nsubstep = glob_nsubstep + step_nsubstep
1621 : glob_nretry = glob_nretry + step_nretry
1622 :
1623 : #endif
1624 :
1625 14592 : if (masterproc) then
1626 19 : if (step_nstep > 0) then
1627 38 : if (do_print) write(LUNOPRT,1) step_max_nsubstep, &
1628 19 : step_nsubstep / step_nstep, &
1629 19 : step_max_nretry, &
1630 38 : step_nretry / step_nstep
1631 : else
1632 0 : if (do_print) write(LUNOPRT,1) step_max_nsubstep, &
1633 0 : 0._r8, &
1634 0 : step_max_nretry, &
1635 0 : 0._r8
1636 : end if
1637 : end if
1638 : end if
1639 :
1640 14592 : end subroutine carma_accumulate_stats
1641 :
1642 :
1643 : !! Set initial mass mixing ratios of constituents, if nothing is specifed
1644 : !! in the initial conditions file.
1645 : !!
1646 : !! NOTE: This call is part of the CAM Physics Interface
1647 : !!
1648 : !! @author Chuck Bardeen
1649 : !! @version May-2009
1650 92160 : subroutine carma_init_cnst(name, latvals, lonvals, mask, q)
1651 : implicit none
1652 :
1653 : character(len=*), intent(in) :: name !! constituent name
1654 : real(r8), intent(in) :: latvals(:) !! lat in degrees (ncol)
1655 : real(r8), intent(in) :: lonvals(:) !! lon in degrees (ncol)
1656 : logical, intent(in) :: mask(:) !! Only initialize where .true.
1657 : real(r8), intent(out) :: q(:,:) !! mass mixing ratio (gcol, lev)
1658 :
1659 : integer :: igroup ! group index
1660 : integer :: ielem ! element index
1661 : integer :: ilev ! level index
1662 : integer :: ibin ! bin index
1663 : integer :: icnst ! constituent index
1664 : integer :: rc ! CARMA return code
1665 : integer :: cnsttype ! constituent type
1666 : integer :: maxbin ! last prognostic bin
1667 :
1668 : ! Initialize the return code.
1669 92160 : rc = 0
1670 :
1671 : ! Determine the element an bin for the particle
1672 737280 : do ielem = 1, NELEM
1673 13639680 : do ibin = 1, NBIN
1674 :
1675 12902400 : call CARMAELEMENT_Get(carma, ielem, rc, igroup=igroup)
1676 12902400 : if (rc < 0) call endrun('carma_timestep_tend::CARMAELEMENT_Get failed.')
1677 :
1678 12902400 : call CARMAGROUP_Get(carma, igroup, rc, cnsttype=cnsttype, maxbin=maxbin)
1679 12902400 : if (rc < 0) call endrun('carma_timestep_tend::CARMAGROUP_Get failed.')
1680 :
1681 39352320 : if (cnsttype == I_CNSTTYPE_PROGNOSTIC) then
1682 :
1683 : ! Bins past maxbin are treated as diagnostic even if the group
1684 : ! is prognostic and thus are not advected in the paerent model.
1685 12902400 : if (ibin <= maxbin) then
1686 :
1687 12902400 : icnst = icnst4elem(ielem, ibin)
1688 :
1689 12902400 : if (cnst_name(icnst) == name) then
1690 :
1691 : ! By default, initialize all constituents to 0.
1692 6543360 : do ilev = 1, size(q, 2)
1693 167823360 : where(mask)
1694 6451200 : q(:, ilev) = 0.0_r8
1695 : end where
1696 : end do
1697 :
1698 92160 : call CARMAMODEL_InitializeParticle(carma, ielem, ibin, latvals, lonvals, mask, q, rc)
1699 92160 : if (rc < 0) call endrun('carma_init_cnst::CARMA_InitializeParticle failed.')
1700 : end if
1701 : end if
1702 : end if
1703 : end do
1704 : end do
1705 :
1706 : ! NOTE: There is currently no initialization for gases, but it could be
1707 : ! added here.
1708 :
1709 92160 : return
1710 : end subroutine carma_init_cnst
1711 :
1712 : !! Calculate amounts of cloudborne aerosols for use in budget diagnostics. This should
1713 : !! be called before the timestep, and the results passed to CARMA_output_cloudborne_diagnostics()
1714 : !! after the timestep to calculate the tendencies and write them out the the history files.
1715 : !!
1716 : !! NOTE: The exact fields that are calculated are determined by the particular CARMA model.
1717 : !!
1718 : !! @author Chuck Bardeen
1719 : !! @version January-2023
1720 0 : subroutine carma_calculate_cloudborne_diagnostics(state, pbuf, aerclddiag)
1721 :
1722 : implicit none
1723 :
1724 : type(physics_state), intent(in) :: state !! Physics state variables - before CARMA
1725 : type(physics_buffer_desc), pointer :: pbuf(:) !! physics buffer
1726 : real(r8), intent(out) :: aerclddiag(pcols, MAXCLDAERDIAG) !! previous cloudborne diagnostics
1727 :
1728 : integer :: rc
1729 :
1730 0 : call CARMAMODEL_CalculateCloudborneDiagnostics(carma, state, pbuf, aerclddiag, rc)
1731 :
1732 0 : return
1733 : end subroutine carma_calculate_cloudborne_diagnostics
1734 :
1735 :
1736 : !! Output cloudborne aerosol budget tendencies to the history files for physics packages
1737 : !! other than CARMA that may be affecting the CARMA aerosols. Since cloudborne aerosols
1738 : !! are not in the physics_state, you must call CARMA_calculate_cloudborne_diagnostics()
1739 : !! before the timestep tend to capture the prior state. This call will calculate the
1740 : !! final state and output the difference as a tendency. This may be useful for
1741 : !! debugging and for calculating aerosol budgets.
1742 : !!
1743 : !! @author Chuck Bardeen
1744 : !! @version January-2023
1745 0 : subroutine carma_output_cloudborne_diagnostics(state, pbuf, pname, dt, oldaerclddiag)
1746 :
1747 : implicit none
1748 :
1749 : type(physics_state), intent(in) :: state !! Physics state variables - before CARMA
1750 : type(physics_buffer_desc), pointer :: pbuf(:) !! physics buffer
1751 : character(*), intent(in) :: pname !! short name of the physics package
1752 : real(r8), intent(in) :: dt !! timestep (s)
1753 : real(r8), intent(in) :: oldaerclddiag(pcols, MAXCLDAERDIAG) !! previous cloudborne diagnostics
1754 :
1755 : integer :: rc
1756 : integer :: i
1757 :
1758 : ! Check to make sure the the package is in the packages list.
1759 0 : do i = 1, carma_ndiagpkgs
1760 0 : if (trim(carma_diags_packages(i)) .eq. trim(pname)) then
1761 :
1762 : ! Allow models to output their own diagnostics related to aerosol
1763 : ! budgets related to physics packages other than CARMA
1764 0 : call CARMAMODEL_OutputCloudborneDiagnostics(carma, state, pbuf, dt, pname, oldaerclddiag, rc)
1765 0 : exit
1766 : end if
1767 : end do
1768 :
1769 0 : return
1770 : end subroutine carma_output_cloudborne_diagnostics
1771 :
1772 :
1773 :
1774 : !! Output budget tendencies to the history files for physics packages
1775 : !! other than CARMA that may be affecting the CARMA aerosols. This can be
1776 : !! called for any physics package that is using ptend to modify the CARMA
1777 : !! aerosol, and may be useful for debugging and for calculating aerosol budgets.
1778 : !!
1779 : !! All the columns in the chunk should be output at the same time.
1780 : !!
1781 : !! @author Chuck Bardeen
1782 : !! @version January-2023
1783 0 : subroutine carma_output_budget_diagnostics(state, ptend, old_cflux, cflux, dt, pname)
1784 :
1785 : implicit none
1786 :
1787 : type(physics_state), intent(in) :: state !! Physics state variables - before CARMA
1788 : type(physics_ptend), intent(in) :: ptend !! indivdual parameterization tendencies
1789 : real(r8) :: old_cflux(pcols,pcnst) !! cam_in%clfux from before the timestep_tend
1790 : real(r8) :: cflux(pcols,pcnst) !! cam_in%clfux from after the timestep_tend
1791 : real(r8), intent(in) :: dt !! timestep (s)
1792 : character(*), intent(in) :: pname !! short name of the physics package
1793 :
1794 : integer :: rc
1795 :
1796 : integer :: i
1797 :
1798 : ! Check to make sure the the package is in the packages list.
1799 0 : do i = 1, carma_ndiagpkgs
1800 0 : if (trim(carma_diags_packages(i)) .eq. trim(pname)) then
1801 :
1802 : ! Allow models to output their own diagnostics related to aerosol
1803 : ! budgets related to physics packages other than CARMA
1804 0 : call CARMAMODEL_OutputBudgetDiagnostics(carma, icnst4elem, icnst4gas, state, ptend, old_cflux, cflux, dt, pname, rc)
1805 0 : exit
1806 : end if
1807 : end do
1808 :
1809 0 : return
1810 : end subroutine carma_output_budget_diagnostics
1811 :
1812 : !! Outputs tracer tendencies and diagnositc fields to the history files.
1813 : !! All the columns in the chunk should be output at the same time.
1814 : !!
1815 : !! @author Chuck Bardeen
1816 : !! @version May-2009
1817 72960 : subroutine carma_output_diagnostics(state, ptend, pbuf, cam_in, gpdiags, sbdiags, gsdiags, spdiags, bndiags)
1818 : use cam_history, only: outfld
1819 : use camsrfexch, only: cam_in_t
1820 :
1821 : implicit none
1822 :
1823 : type(physics_state), intent(in) :: state !! Physics state variables - before CARMA
1824 : type(physics_ptend), intent(in) :: ptend !! indivdual parameterization tendencies
1825 : type(physics_buffer_desc), pointer :: pbuf(:) !! physics buffer
1826 : type(cam_in_t), intent(in) :: cam_in !! surface inputs
1827 : real(r8), intent(in), dimension(pcols, pver, NGROUP, NGPDIAGS) :: gpdiags !! CARMA group diagnostic output
1828 : real(r8), intent(in), dimension(pcols, NBIN, NELEM, NSBDIAGS) :: sbdiags !! CARMA surface bin diagnostic output
1829 : real(r8), intent(in), dimension(pcols, pver, NGAS, NGSDIAGS) :: gsdiags !! CARMA gas diagnostic output
1830 : real(r8), intent(in), dimension(pcols, pver, NSPDIAGS) :: spdiags !! CARMA step diagnostic output
1831 : real(r8), intent(in), dimension(pcols, pver, NBIN, NELEM, NBNDIAGS) :: bndiags !! CARMA bin diagnostic output
1832 :
1833 : ! Local variables
1834 : integer :: igroup ! group index
1835 : integer :: ielem ! element index
1836 : integer :: ibin ! bin index
1837 : integer :: igas ! gas index
1838 : integer :: ienconc ! element index for group's concentration element
1839 : integer :: icnst ! constituent index
1840 : integer :: lchnk ! chunk identifier
1841 : integer :: rc ! CARMA return code
1842 : character(len=8) :: sname ! short (CAM) name
1843 : integer :: cnsttype ! constituent type
1844 : integer :: maxbin ! last prognostic bin
1845 : logical :: is_cloud ! is the group a cloud?
1846 : logical :: do_drydep ! is dry deposition enabled?
1847 :
1848 : character(len=*), parameter :: subname = 'carma_output_diagnostics'
1849 :
1850 : ! Initialize the return code.
1851 72960 : rc = 0
1852 :
1853 : ! Check each column int the chunk.
1854 72960 : lchnk = state%lchnk
1855 :
1856 : ! Output step diagnostics.
1857 72960 : if (carma_do_substep) then
1858 72960 : call outfld('CRNSTEP', spdiags(:, :, SPDIAGS_NSTEP), pcols, lchnk)
1859 72960 : call outfld('CRLNSTEP', spdiags(:, :, SPDIAGS_LNSTEP), pcols, lchnk)
1860 : end if
1861 :
1862 : ! Output the particle tendencies.
1863 583680 : do ielem = 1, NELEM
1864 10798080 : do ibin = 1, NBIN
1865 :
1866 10214400 : call CARMAELEMENT_Get(carma, ielem, rc, igroup=igroup)
1867 10214400 : if (rc < 0) call endrun(subname//'::CARMAELEMENT_Get failed.')
1868 :
1869 10214400 : call CARMAGROUP_Get(carma, igroup, rc, cnsttype=cnsttype, maxbin=maxbin, do_drydep=do_drydep)
1870 10214400 : if (rc < 0) call endrun(subname//'::CARMAGROUP_Get failed.')
1871 :
1872 31153920 : if (cnsttype == I_CNSTTYPE_PROGNOSTIC) then
1873 :
1874 : ! Bins past maxbin are treated as diagnostic even if the group
1875 : ! is prognostic and thus are not advected in the paerent model.
1876 10214400 : if (ibin <= maxbin) then
1877 :
1878 10214400 : icnst = icnst4elem(ielem, ibin)
1879 :
1880 10214400 : call outfld(trim(etndname(ielem, ibin))//'TC', ptend%q(:, :, icnst), pcols, lchnk)
1881 :
1882 10214400 : if (do_drydep) then
1883 10214400 : call outfld(trim(etndname(ielem, ibin))//'DD', sbdiags(:, ibin, ielem, SBDIAGS_DD), pcols, lchnk)
1884 : end if
1885 :
1886 10214400 : if (carma_do_pheat) then
1887 :
1888 : ! Only specified for the number density element of the group.
1889 0 : if (bndiags(1, 1, ibin, ielem, BNDIAGS_TP) /= CAM_FILL) then
1890 0 : call outfld(trim(etndname(ielem, ibin))//'TP', bndiags(:, :, ibin, ielem, BNDIAGS_TP), pcols, lchnk)
1891 : end if
1892 : end if
1893 : end if
1894 : end if
1895 : end do
1896 : end do
1897 :
1898 : ! Output the particle diagnostics.
1899 218880 : do igroup = 1, NGROUP
1900 145920 : call CARMAGROUP_Get(carma, igroup, rc, shortname=sname, is_cloud=is_cloud, do_drydep=do_drydep, ienconc=ienconc)
1901 145920 : if (rc < 0) call endrun(subname//'::CARMAGROUP_Get failed.')
1902 :
1903 : ! Gridbox average
1904 145920 : call outfld(trim(sname)//'ND', gpdiags(:, :, igroup, GPDIAGS_ND), pcols, lchnk)
1905 145920 : call outfld(trim(sname)//'AD', gpdiags(:, :, igroup, GPDIAGS_AD), pcols, lchnk)
1906 145920 : call outfld(trim(sname)//'MD', gpdiags(:, :, igroup, GPDIAGS_MD), pcols, lchnk)
1907 145920 : call outfld(trim(sname)//'RE', gpdiags(:, :, igroup, GPDIAGS_RE), pcols, lchnk)
1908 145920 : call outfld(trim(sname)//'RM', gpdiags(:, :, igroup, GPDIAGS_RM), pcols, lchnk)
1909 145920 : call outfld(trim(sname)//'JN', gpdiags(:, :, igroup, GPDIAGS_JN), pcols, lchnk)
1910 145920 : call outfld(trim(sname)//'MR', gpdiags(:, :, igroup, GPDIAGS_MR), pcols, lchnk)
1911 145920 : call outfld(trim(sname)//'EX', gpdiags(:, :, igroup, GPDIAGS_EX), pcols, lchnk)
1912 145920 : call outfld(trim(sname)//'OD', gpdiags(:, :, igroup, GPDIAGS_OD), pcols, lchnk)
1913 145920 : call outfld(trim(sname)//'PA', gpdiags(:, :, igroup, GPDIAGS_PA), pcols, lchnk)
1914 145920 : call outfld(trim(sname)//'AR', gpdiags(:, :, igroup, GPDIAGS_AR), pcols, lchnk)
1915 145920 : call outfld(trim(sname)//'VM', gpdiags(:, :, igroup, GPDIAGS_VM), pcols, lchnk)
1916 145920 : call outfld(trim(sname)//'VR', gpdiags(:, :, igroup, GPDIAGS_VR), pcols, lchnk)
1917 :
1918 145920 : if (do_drydep) then
1919 3064320 : do ibin = 1, NBIN
1920 3064320 : call outfld(trim(btndname(igroup, ibin))//'VD', sbdiags(:, ibin, ienconc, SBDIAGS_VD), pcols, lchnk)
1921 : end do
1922 : end if
1923 :
1924 3283200 : do ibin = 1,NBIN
1925 2918400 : call outfld(trim(btndname(igroup, ibin))//'ND',bndiags(:, :, ibin, ienconc, BNDIAGS_ND), pcols, lchnk)
1926 2918400 : call outfld(trim(btndname(igroup, ibin))//'WR',bndiags(:, :, ibin, ienconc, BNDIAGS_WETR), pcols, lchnk)
1927 2918400 : call outfld(trim(btndname(igroup, ibin))//'RO',bndiags(:, :, ibin, ienconc, BNDIAGS_RO), pcols, lchnk)
1928 3064320 : call outfld(trim(btndname(igroup, ibin))//'VR',bndiags(:, :, ibin, ienconc, BNDIAGS_VR), pcols, lchnk)
1929 : end do
1930 : end do
1931 :
1932 : ! Output the gas tendencies.
1933 218880 : do igas = 1, NGAS
1934 145920 : icnst = icnst4gas(igas)
1935 :
1936 145920 : call outfld(gtndname(igas), ptend%q(:, :, icnst), pcols, lchnk)
1937 :
1938 : ! Output the supersaturations.
1939 145920 : call outfld(trim(cnst_name(icnst))//'SI', gsdiags(:, :, igas, GSDIAGS_SI), pcols, lchnk)
1940 145920 : call outfld(trim(cnst_name(icnst))//'SL', gsdiags(:, :, igas, GSDIAGS_SL), pcols, lchnk)
1941 145920 : call outfld(trim(cnst_name(icnst))//'EI', gsdiags(:, :, igas, GSDIAGS_EI), pcols, lchnk)
1942 145920 : call outfld(trim(cnst_name(icnst))//'EL', gsdiags(:, :, igas, GSDIAGS_EL), pcols, lchnk)
1943 218880 : call outfld(trim(cnst_name(icnst))//'WT', gsdiags(:, :, igas, GSDIAGS_WT), pcols, lchnk)
1944 : end do
1945 :
1946 : ! Output the temperature tendency.
1947 72960 : if (carma_do_thermo) then
1948 0 : call outfld('CRTT', ptend%s(:, :) / cpair, pcols, lchnk)
1949 : end if
1950 :
1951 : ! Allow models to output their own diagnostics
1952 72960 : call CARMAMODEL_OutputDiagnostics(carma, icnst4elem, state, ptend, pbuf, cam_in, rc)
1953 :
1954 72960 : return
1955 72960 : end subroutine carma_output_diagnostics
1956 :
1957 : !! Calculate the emissions for CARMA aerosols. This is taken from
1958 : !! the routine aerosol_emis_intr in aerosol_intr.F90 and dust_emis_intr in
1959 : !! dust_intr.F90 by Phil Rasch.
1960 : !!
1961 : !! @author Chuck Bardeen
1962 : !! @version May-2009
1963 72960 : subroutine carma_emission_tend (state, ptend, cam_in, dt, pbuf)
1964 72960 : use cam_history, only: outfld
1965 : use camsrfexch, only: cam_in_t
1966 :
1967 : implicit none
1968 :
1969 : type(physics_state), intent(in ) :: state !! physics state
1970 : type(physics_ptend), intent(inout) :: ptend !! physics state tendencies
1971 : type(cam_in_t), intent(inout) :: cam_in !! surface inputs
1972 : real(r8), intent(in) :: dt !! time step (s)
1973 : type(physics_buffer_desc), pointer :: pbuf(:) !! physics buffer
1974 :
1975 : integer :: lchnk ! chunk identifier
1976 : integer :: ncol ! number of columns in chunk
1977 : integer :: igroup ! group index
1978 : integer :: ielem ! element index
1979 : integer :: ibin ! bin index
1980 : integer :: icnst ! consituent index
1981 : real(r8) :: tendency(pcols, pver) ! constituent tendency (kg/kg/s)
1982 : real(r8) :: surfaceFlux(pcols) ! constituent surface flux (kg/m^2/s)
1983 : integer :: cnsttype ! constituent type
1984 : integer :: maxbin ! last prognostic bin
1985 : integer :: rc ! CARMA return code
1986 :
1987 : ! Initialize the return code.
1988 72960 : rc = 0
1989 :
1990 : ! Initialize the output tendency structure.
1991 72960 : call physics_ptend_init(ptend,state%psetcols, 'CARMA (emission)', lq=lq_carma)
1992 :
1993 72960 : if (.not. carma_flag) return
1994 72960 : if (.not. carma_do_emission) return
1995 :
1996 72960 : ncol = state%ncol
1997 72960 : lchnk = state%lchnk
1998 :
1999 : ! Provide emissions rates for particles.
2000 : !
2001 : ! NOTE: This can only be done for prognostic groups.
2002 583680 : do ielem = 1, NELEM
2003 510720 : call CARMAELEMENT_Get(carma, ielem, rc, igroup=igroup)
2004 510720 : if (rc < 0) call endrun('carma_drydep_tend::CARMAELEMENT_Get failed.')
2005 :
2006 510720 : call CARMAGROUP_Get(carma, igroup, rc, cnsttype=cnsttype, maxbin=maxbin)
2007 510720 : if (rc < 0) call endrun('carma_drydep_tend::CARMAGROUP_Get failed.')
2008 :
2009 1605120 : if (cnsttype == I_CNSTTYPE_PROGNOSTIC) then
2010 :
2011 10725120 : do ibin = 1, NBIN
2012 :
2013 : ! Bins past maxbin are treated as diagnostic even if the group
2014 : ! is prognostic and thus are not advected in the paerent model.
2015 10725120 : if (ibin <= maxbin) then
2016 :
2017 10214400 : icnst = icnst4elem(ielem, ibin)
2018 :
2019 10214400 : call CARMAMODEL_EmitParticle(carma, ielem, ibin, icnst, dt, state, cam_in, tendency, surfaceFlux, pbuf, rc)
2020 10214400 : if (rc < 0) call endrun('carma_emission_tend::CARMA_EmitParticle failed.')
2021 :
2022 : ! Add any surface flux here.
2023 157301760 : cam_in%cflx(:ncol, icnst) = surfaceFlux(:ncol)
2024 10214400 : call outfld(trim(cnst_name(icnst))//'SF', cam_in%cflx(:ncol, icnst), ncol, lchnk)
2025 :
2026 : ! For emissions into the atmosphere, put the emission here.
2027 11021337600 : ptend%q(:ncol, :pver, icnst) = tendency(:ncol, :pver)
2028 11021337600 : call outfld(trim(cnst_name(icnst))//'EM', ptend%q(:ncol, :, icnst), ncol, lchnk)
2029 : end if
2030 : enddo
2031 : end if
2032 : enddo
2033 :
2034 : ! No emissions rate is set up for gases, but it could be added here.
2035 :
2036 : return
2037 72960 : end subroutine carma_emission_tend
2038 :
2039 :
2040 : !! Calculate the wet deposition for the CARMA aerosols. This is taken from
2041 : !! the routine aerosol_wet_int in aerosol_intr.F90 and dust_wet_intr in
2042 : !! dust_intr.F90 by Phil Rasch.
2043 : !!
2044 : !! Method:
2045 : !! Use a modified version of the scavenging parameterization described in
2046 : !! Barth et al, 2000, JGR (sulfur cycle paper)
2047 : !! Rasch et al, 2001, JGR (INDOEX paper)
2048 : !!
2049 : !! @author Chuck Bardeen
2050 : !! @version May-2009
2051 0 : subroutine carma_wetdep_tend(state, ptend, dt, pbuf, dlf, cam_out)
2052 72960 : use cam_history, only: outfld
2053 : use phys_control, only: cam_physpkg_is
2054 : use wetdep, only: clddiag, wetdepa_v1, wetdepa_v2
2055 : use camsrfexch, only: cam_out_t
2056 : use physconst, only: gravit
2057 :
2058 : implicit none
2059 :
2060 : real(r8), intent(in) :: dt !! time step (s)
2061 : type(physics_state), intent(in ) :: state !! physics state
2062 : type(physics_ptend), intent(inout) :: ptend !! physics state tendencies
2063 : type(physics_buffer_desc), pointer :: pbuf(:) !! physics buffer
2064 : real(r8), intent(in) :: dlf(pcols,pver) !! Detrainment of convective condensate
2065 : type(cam_out_t), intent(inout) :: cam_out !! cam output to surface models
2066 :
2067 : ! local vars
2068 : real(r8) :: rainmr(pcols,pver) ! mixing ratio of rain within cloud volume
2069 : real(r8) :: cldv(pcols,pver) ! cloudy volume undergoing wet chem and scavenging
2070 : real(r8) :: cldvcu(pcols,pver) ! Convective precipitation area, top interface of current layer
2071 : real(r8) :: cldvst(pcols,pver) ! Stratiform precipitation area, top interface of current layer
2072 : integer :: ielem ! element index
2073 : integer :: igroup ! group index
2074 : integer :: ibin ! bin index
2075 : integer :: icnst ! constituent index
2076 : real(r8) :: conicw(pcols,pver) ! convective in-cloud water
2077 : real(r8) :: cmfdqr(pcols,pver) ! convective production of rain
2078 : real(r8) :: cldc(pcols,pver) ! convective cloud fraction, currently empty
2079 : real(r8) :: clds(pcols,pver) ! Stratiform cloud fraction
2080 : real(r8) :: evapc(pcols,pver) ! Evaporation rate of convective precipitation
2081 : real(r8) :: iscavt(pcols, pver)
2082 : real(r8) :: scavt(pcols, pver)
2083 : integer :: ixcldliq
2084 : integer :: ixcldice
2085 : real(r8) :: totcond(pcols, pver) ! total condensate
2086 : real(r8) :: solfac(pcols, pver) ! solubility factor
2087 : real(r8) :: solfac_in ! solubility factor
2088 : real(r8) :: scavcoef ! scavenging Coefficient
2089 : logical :: do_wetdep
2090 : integer :: ncol ! number of columns
2091 : integer :: lchnk ! chunk identifier
2092 : integer :: rc ! CARMA return code
2093 : real(r8) :: z_scavcoef(pcols,pver) ! Dana and Hales coefficient (/mm)
2094 : integer :: cnsttype ! constituent type
2095 : integer :: k
2096 : real(r8) :: sflx(pcols) ! Surface Flux (kg/m2/s)
2097 : integer :: maxbin
2098 :
2099 : ! physics buffer
2100 : integer itim_old
2101 0 : real(r8), pointer, dimension(:,:) :: cldn ! cloud fraction
2102 0 : real(r8), pointer, dimension(:,:) :: cme
2103 0 : real(r8), pointer, dimension(:,:) :: prain
2104 0 : real(r8), pointer, dimension(:,:) :: evapr
2105 0 : real(r8), pointer, dimension(:,:) :: icwmrdp ! in cloud water mixing ratio, deep convection
2106 0 : real(r8), pointer, dimension(:,:) :: rprddp ! rain production, deep convection
2107 0 : real(r8), pointer, dimension(:,:) :: icwmrsh ! in cloud water mixing ratio, deep convection
2108 0 : real(r8), pointer, dimension(:,:) :: rprdsh ! rain production, deep convection
2109 0 : real(r8), pointer, dimension(:,:,:) :: fracis ! fraction of transported species that are insoluble
2110 0 : real(r8), pointer, dimension(:,:) :: sh_frac ! Shallow convective cloud fraction
2111 0 : real(r8), pointer, dimension(:,:) :: dp_frac ! Deep convective cloud fraction
2112 0 : real(r8), pointer, dimension(:,:) :: evapcsh ! Evaporation rate of shallow convective precipitation >=0.
2113 0 : real(r8), pointer, dimension(:,:) :: evapcdp ! Evaporation rate of deep convective precipitation >=0.
2114 :
2115 : ! Initialize the return code.
2116 0 : rc = 0
2117 :
2118 : ! Initialize the output tendency structure.
2119 0 : call physics_ptend_init(ptend,state%psetcols, 'CARMA (wetdep)', lq=lq_carma)
2120 :
2121 0 : if (.not. carma_flag) return
2122 0 : if (.not. carma_do_wetdep) return
2123 :
2124 0 : ncol = state%ncol
2125 0 : lchnk = state%lchnk
2126 :
2127 : ! Associate pointers with physics buffer fields
2128 0 : itim_old = pbuf_old_tim_idx()
2129 :
2130 0 : call pbuf_get_field(pbuf, pbuf_get_index('CLD'), cldn, (/1,1,itim_old/),(/pcols,pver,1/))
2131 0 : call pbuf_get_field(pbuf, pbuf_get_index('QME'), cme )
2132 0 : call pbuf_get_field(pbuf, pbuf_get_index('PRAIN'), prain )
2133 0 : call pbuf_get_field(pbuf, pbuf_get_index('NEVAPR'), evapr )
2134 0 : call pbuf_get_field(pbuf, pbuf_get_index('FRACIS'), fracis )
2135 0 : call pbuf_get_field(pbuf, pbuf_get_index('ICWMRDP'), icwmrdp )
2136 0 : call pbuf_get_field(pbuf, pbuf_get_index('RPRDDP'), rprddp )
2137 0 : call pbuf_get_field(pbuf, pbuf_get_index('ICWMRSH'), icwmrsh )
2138 0 : call pbuf_get_field(pbuf, pbuf_get_index('RPRDSH'), rprdsh )
2139 :
2140 : ! sum deep and shallow convection contributions
2141 0 : conicw(:ncol,:) = icwmrdp(:ncol,:) + icwmrsh(:ncol,:)
2142 0 : cmfdqr(:ncol,:) = rprddp(:ncol,:) + rprdsh(:ncol,:)
2143 :
2144 0 : call pbuf_get_field(pbuf, pbuf_get_index('SH_FRAC'), sh_frac )
2145 0 : call pbuf_get_field(pbuf, pbuf_get_index('DP_FRAC'), dp_frac )
2146 0 : call pbuf_get_field(pbuf, pbuf_get_index('NEVAPR_SHCU'), evapcsh )
2147 0 : call pbuf_get_field(pbuf, pbuf_get_index('NEVAPR_DPCU'), evapcdp )
2148 :
2149 0 : cldc(:ncol,:) = dp_frac(:ncol,:) + sh_frac(:ncol,:) ! Sungsu included this.
2150 0 : evapc(:ncol,:) = evapcsh(:ncol,:) + evapcdp(:ncol,:) ! Sungsu included this.
2151 0 : clds(:ncol,:) = cldn(:ncol,:) - cldc(:ncol,:) ! Stratiform cloud fraction
2152 :
2153 :
2154 0 : cmfdqr(:ncol,:) = rprddp(:ncol,:) + rprdsh(:ncol,:)
2155 :
2156 : ! fields needed for wet scavenging
2157 : call clddiag( state%t, state%pmid, state%pdel, cmfdqr, evapc, cldn, cldc, clds, cme, evapr, prain, &
2158 0 : cldv, cldvcu, cldvst, rainmr, ncol )
2159 :
2160 0 : call cnst_get_ind('CLDICE', ixcldice)
2161 0 : call cnst_get_ind('CLDLIQ', ixcldliq)
2162 0 : totcond(:ncol,:) = state%q(:ncol,:,ixcldliq) + &
2163 0 : state%q(:ncol,:,ixcldice)
2164 :
2165 : ! Iterate over each particle and calculate a tendency from wet
2166 : ! scavenging for it.
2167 0 : do ielem = 1, NELEM
2168 :
2169 : ! NOTE: This can only be done for prognistic groups.
2170 :
2171 0 : call CARMAELEMENT_Get(carma, ielem, rc, igroup=igroup)
2172 0 : if (rc < 0) call endrun('carma_wetdep_tend::CARMAELEMENT_Get failed.')
2173 :
2174 : call CARMAGROUP_Get(carma, igroup, rc, cnsttype=cnsttype, do_wetdep=do_wetdep, &
2175 0 : solfac=solfac_in, scavcoef=scavcoef, maxbin=maxbin)
2176 0 : if (rc < 0) call endrun('carma_wetdep_tend::CARMAGROUP_Get failed.')
2177 :
2178 0 : solfac(:,:) = solfac_in
2179 :
2180 0 : if ((do_wetdep) .and. (cnsttype == I_CNSTTYPE_PROGNOSTIC)) then
2181 :
2182 0 : do ibin = 1, NBIN
2183 :
2184 : ! Bins past maxbin are treated as diagnostic even if the group
2185 : ! is prognostic and thus are not advected in the parent model.
2186 0 : if (ibin <= maxbin) then
2187 :
2188 0 : icnst = icnst4elem(ielem, ibin)
2189 :
2190 0 : scavt = 0._r8
2191 :
2192 : ! The scavenging coefficient might be calculated as a function of
2193 : ! the aerosol bin at each grid point. However, for now, we will just
2194 : ! use a constant value for each group.
2195 0 : z_scavcoef(:, :) = scavcoef
2196 :
2197 0 : if (cam_physpkg_is('cam5') .or. cam_physpkg_is('cam6')) then
2198 :
2199 : call wetdepa_v2( &
2200 : state%pmid, &
2201 : state%q, &
2202 : state%pdel, &
2203 : cldn, &
2204 : cldc, &
2205 : cmfdqr, &
2206 : evapc, &
2207 : conicw, &
2208 : prain, &
2209 : cme, &
2210 : evapr, &
2211 : totcond, &
2212 : state%q(:, :, icnst), &
2213 : dt, &
2214 : scavt, &
2215 : iscavt, &
2216 : cldvcu, &
2217 : cldvst, &
2218 : dlf, &
2219 : fracis(:, :, icnst), &
2220 : solfac, &
2221 : ncol, &
2222 0 : z_scavcoef)
2223 :
2224 0 : else if (cam_physpkg_is('cam4')) then
2225 :
2226 : call wetdepa_v1(state%t, &
2227 : state%pmid, &
2228 : state%q, &
2229 : state%pdel, &
2230 : cldn, &
2231 : cldc, &
2232 : cmfdqr, &
2233 : conicw, &
2234 : prain, &
2235 : cme, &
2236 : evapr, &
2237 : totcond, &
2238 : state%q(:, :, icnst), &
2239 : dt, &
2240 : scavt, &
2241 : iscavt, &
2242 : cldv, &
2243 : fracis(:, :, icnst), &
2244 : solfac_in, &
2245 : ncol, &
2246 0 : z_scavcoef)
2247 : else
2248 :
2249 0 : call endrun('carma_wetdep_tend:: No wet deposition routine is available for this configuration.')
2250 : end if
2251 :
2252 0 : ptend%q(:, :, icnst) = scavt
2253 0 : call outfld(trim(cnst_name(icnst))//'WD', ptend%q(:, :, icnst), pcols, lchnk)
2254 :
2255 : !
2256 : ! ptend%q(kg/kg air/s) * pdel(Pa) / gravit (m/s2) => (kg/m2/s)
2257 : ! note: 1Pa = 1 kg air * (m/s2) / m2
2258 0 : sflx(:) = 0._r8
2259 :
2260 0 : do k = 1,pver
2261 0 : sflx(:ncol) = sflx(:ncol) - ptend%q(:ncol, k, icnst) * state%pdel(:ncol,k) / gravit
2262 : enddo
2263 :
2264 0 : call outfld(trim(cnst_name(icnst))//'SW', sflx, pcols, lchnk)
2265 :
2266 : ! Add this to the surface amount of the constituent
2267 0 : call CARMAMODEL_WetDeposition(carma, ielem, ibin, sflx, cam_out, state, rc)
2268 :
2269 : end if
2270 : end do
2271 : end if
2272 : end do
2273 :
2274 : return
2275 0 : end subroutine carma_wetdep_tend
2276 :
2277 :
2278 : !! This routine creates files containing optical properties for each radiatively
2279 : !! active particle type. These optical properties are used by the RRTMG radiation
2280 : !! code to include the impact of CARMA particles in the radiative transfer
2281 : !! calculation.
2282 : !!
2283 : !! The opticsType that is specified for the group determines how the optical
2284 : !! properties will be generated for that group. Each group can use a different
2285 : !! optics method if needed. Refractive indices need for these calculation are
2286 : !! are specified in the group's elements rather than at the group level. This
2287 : !! allows various mixing approaches to be used to determine the refractive index
2288 : !! for the particle as a whole.
2289 : !!
2290 : !! The I_OPTICS_MIXED_YU2105 and I_OPTICS_SULFATE_YU2015 optics methods are
2291 : !! designed to trop_strat models as define in the Yu et al. (2015) paper. The
2292 : !! other optics types can be applied more generically to a number of different
2293 : !! aerosol/cloud models.
2294 : !!
2295 : !! NOTE: The format of this file is determined by the needs of the radiative tranfer
2296 : !! code, so ideally a routine would exist in that module that could create a file
2297 : !! with the proper format. Since that doesn't exist, we do it all here.
2298 0 : subroutine CARMA_CreateOpticsFile(carma, rc)
2299 :
2300 : implicit none
2301 :
2302 : type(carma_type), intent(inout) :: carma !! the carma object
2303 : integer, intent(out) :: rc !! return code, negative indicates failure
2304 :
2305 : ! Local variables
2306 : integer :: igroup
2307 : logical :: do_mie
2308 : integer :: cnsttype ! constituent type
2309 : integer :: opticsType
2310 :
2311 : ! Assume success.
2312 0 : rc = 0
2313 :
2314 : ! Process each group that is defined in the model.
2315 0 : do igroup = 1, NGROUP
2316 :
2317 : ! Get the necessary group properties.
2318 0 : call CARMAGROUP_Get(carma, igroup, rc, do_mie=do_mie, cnsttype=cnsttype, iopticstype=opticsType)
2319 0 : if (rc < 0) call endrun('carma_CreateOpticsFile::CARMAGROUP_Get failed.')
2320 :
2321 : ! Are we supposed to do the mie calculation for this group?
2322 0 : if ((do_mie) .and. (cnsttype == I_CNSTTYPE_PROGNOSTIC)) then
2323 :
2324 : ! What type of calculation is needed for this group?
2325 : !
2326 : ! NOTE: Some of these calculations generate optical properties as single mass
2327 : ! coefficients, while others are lookup tables designed around multiple
2328 : ! dimensions.
2329 0 : select case (opticsType)
2330 :
2331 : ! This is for fixed composition, but the particle may swell in response
2332 : ! to changes in RH. Only one refractive index specified at the group level.
2333 : !
2334 : ! NOTE: This is what was used by the first CARMA models that were radiatively
2335 : ! active.
2336 : case (I_OPTICS_FIXED)
2337 0 : call CARMA_CreateOpticsFile_Fixed(carma, igroup, rc)
2338 0 : if (rc < 0) call endrun('carma_CreateOpticsFile::CreateOpticsFile_Fixed failed.')
2339 :
2340 : ! This is similar to Yu (2015) in that handles mixed particles treated as
2341 : ! core shell particles; however the dimensions of the lookup table are the
2342 : ! the radii and the refractive indicies, so it can be used with various
2343 : ! aerosol configurations (not just as in the Yu(2015)).
2344 : case(I_OPTICS_MIXED_CORESHELL)
2345 0 : call endrun('carma_CreateOpticsFile mixed_coreshell has not been implemented.')
2346 :
2347 : ! This is similar to MAM4, in that a volume mixing approach is used to
2348 : ! mixed both the core and the shell together and thus only one radius and
2349 : ! one refractive index are needed in the lookup table.
2350 : case(I_OPTICS_MIXED_VOLUME)
2351 0 : call endrun('carma_CreateOpticsFile mixed_volume has not been implemented.')
2352 :
2353 : ! This is similar to "mixed_volume", except that Maxwell-Garnett mixing
2354 : ! is used instead of volume mixing.
2355 : case(I_OPTICS_MIXED_MAXWELL)
2356 0 : call endrun('carma_CreateOpticsFile mixed_maxwell has not been implemented.')
2357 :
2358 : ! This is for a pure sulfate group where the table is based upon weight
2359 : ! percent; however, unlike sulfate_Yu, the refractive index of the sulfate
2360 : ! changes with the weight percent of H2SO4.
2361 : case(I_OPTICS_SULFATE)
2362 0 : call CARMA_CreateOpticsFile_Sulfate(carma, igroup, rc)
2363 0 : if (rc < 0) call endrun('carma_CreateOpticsFile::CreateOpticsFile_Sulfate failed.')
2364 :
2365 : ! Other types are not generically useful are are particular to the
2366 : ! specific model, so thos are handled by model specific code. These
2367 : ! include:
2368 : ! I_OPTICS_MIXED_YU2015
2369 : ! I_OPTICS_MIXED_YU_H2O
2370 : ! I_OPTICS_SULFATE_YU2015
2371 : case default
2372 0 : call CARMAMODEL_CreateOpticsFile(carma, igroup, opticsType, rc)
2373 :
2374 : end select
2375 : end if
2376 : end do
2377 :
2378 0 : return
2379 0 : end subroutine CARMA_CreateOpticsFile
2380 :
2381 :
2382 : !! This routine creates files containing optical properties for each radiatively
2383 : !! active particle type. These optical properties are used by the RRTMG radiation
2384 : !! code to include the impact of CARMA particles in the radiative transfer
2385 : !! calculation.
2386 : !!
2387 : !! NOTE: The format of this file is determined by the needs of the radiative tranfer
2388 : !! code, so ideally a routine would exist in that module that could create a file
2389 : !! with the proper format. Since that doesn't exist, we do it all here.
2390 0 : subroutine CARMA_CreateOpticsFile_Fixed(carma, igroup, rc)
2391 : use wrap_nf
2392 : use wetr, only : getwetr
2393 :
2394 : implicit none
2395 :
2396 : type(carma_type), intent(inout) :: carma !! the carma object
2397 : integer, intent(in) :: igroup !! group index
2398 : integer, intent(out) :: rc !! return code, negative indicates failure
2399 :
2400 : ! Local variables
2401 : integer :: ibin, iwave, irh
2402 : integer :: irhswell
2403 : integer :: ienconc
2404 : real(kind=f) :: rho(NBIN), rhopwet
2405 : real(kind=f) :: r(NBIN), rmass(NBIN), rlow(NBIN), rup(NBIN)
2406 : real(kind=f) :: wave(NWAVE)
2407 : complex(kind=f) :: refidx(NWAVE, NREFIDX)
2408 : character(len=CARMA_NAME_LEN) :: name
2409 : character(len=CARMA_SHORT_NAME_LEN) :: shortname
2410 : logical :: do_mie
2411 : integer :: fid
2412 : integer :: rhdim, lwdim, swdim
2413 : integer :: rhvar, lwvar, swvar
2414 : integer :: abs_lw_var
2415 : integer :: ext_sw_var, ssa_sw_var, asm_sw_var
2416 : integer :: omdim, andim, namedim
2417 : integer :: omvar, anvar, namevar
2418 : integer :: dimids(2)
2419 : integer :: denvar, slogvar, dryrvar, rminvar, rmaxvar, hygrovar, ntmvar
2420 : real(kind=f) :: abs_lw(NMIE_RH, nlwbands)
2421 : real(kind=f) :: ext_sw(NMIE_RH, nswbands)
2422 : real(kind=f) :: ssa_sw(NMIE_RH, nswbands)
2423 : real(kind=f) :: asm_sw(NMIE_RH, nswbands)
2424 : character(len=8) :: c_name ! constituent name
2425 : character(len=32) :: aer_name ! long enough for both aername and name
2426 : character(len=255) :: filepath
2427 : real(kind=f) :: rwet
2428 : real(kind=f) :: Qext
2429 : real(kind=f) :: Qsca
2430 : real(kind=f) :: asym
2431 : integer :: start_text(2), count_text(2)
2432 : integer :: sw_r_refidx_var, sw_i_refidx_var, lw_r_refidx_var, lw_i_refidx_var
2433 : integer :: nrh
2434 : integer :: cnsttype ! constituent type
2435 : integer :: maxbin ! last prognostic bin
2436 : integer :: LUNOPRT ! logical unit number for output
2437 : logical :: do_print ! do print output?
2438 : integer :: ret
2439 :
2440 :
2441 : ! Assume success.
2442 0 : rc = 0
2443 :
2444 : ! Get the wavelength structure.
2445 0 : call CARMA_GET(carma, rc, wave=wave, do_print=do_print, LUNOPRT=LUNOPRT)
2446 0 : if (rc < 0) call endrun('carma_CreateOpticsFile::CARMA_Get failed.')
2447 :
2448 : ! Get the necessary group properties.
2449 : call CARMAGROUP_Get(carma, igroup, rc, do_mie=do_mie, name=name, shortname=shortname, r=r, &
2450 : rlow=rlow, rup=rup, rmass=rmass, irhswell=irhswell, &
2451 0 : ienconc=ienconc, cnsttype=cnsttype, maxbin=maxbin)
2452 0 : if (rc < 0) call endrun('carma_CreateOpticsFile::CARMAGROUP_Get failed.')
2453 :
2454 0 : call CARMAELEMENT_Get(carma, ienconc, rc, rho=rho, refidx=refidx)
2455 0 : if (rc < 0) call endrun('carma_CreateOpticsFile::CARMAELEMENT_Get failed.')
2456 :
2457 : ! A file needs to be created for each bin.
2458 0 : do ibin = 1, NBIN
2459 :
2460 : ! Bins past maxbin are treated as diagnostic even if the group
2461 : ! is prognostic and thus are not advected in the paerent model.
2462 0 : if (ibin <= maxbin) then
2463 :
2464 0 : write(c_name, '(A, I2.2)') trim(shortname), ibin
2465 :
2466 : ! Construct the path to the file. Each model will have its own subdirectory
2467 : ! where the optical property files are stored.
2468 0 : filepath = trim(carma_model) // '_' // trim(c_name) // '_rrtmg.nc'
2469 :
2470 0 : if (do_print) write(LUNOPRT,*) 'Creating CARMA optics file ... ', trim(filepath)
2471 :
2472 : ! Create the file.
2473 0 : call wrap_create(filepath, NF90_CLOBBER, fid)
2474 :
2475 : ! For non-hygroscopic, only use 1 RH value.
2476 0 : if (irhswell /= 0) then
2477 0 : nrh = NMIE_RH
2478 : else
2479 0 : nrh = min(NMIE_RH, 1)
2480 : end if
2481 :
2482 : ! Define the dimensions: rh, lwbands, swbands
2483 0 : call wrap_def_dim(fid, 'rh_idx', nrh, rhdim)
2484 0 : call wrap_def_dim(fid, 'lw_band', nlwbands, lwdim)
2485 0 : call wrap_def_dim(fid, 'sw_band', nswbands, swdim)
2486 :
2487 0 : write(LUNOPRT,*) "Defined rh_idx, lw_band, and sw_band dims."
2488 :
2489 0 : dimids(1) = rhdim
2490 0 : call wrap_def_var(fid, 'rh', NF90_DOUBLE, 1, dimids(1:1), rhvar)
2491 :
2492 0 : dimids(1) = lwdim
2493 0 : call wrap_def_var(fid, 'lw_band', NF90_DOUBLE, 1, dimids(1:1), lwvar)
2494 :
2495 0 : dimids(1) = swdim
2496 0 : call wrap_def_var(fid, 'sw_band', NF90_DOUBLE, 1, dimids(1:1), swvar)
2497 :
2498 0 : write(LUNOPRT,*) "Defined rh_idx, lw_band, and sw_band vars."
2499 :
2500 0 : call wrap_put_att_text(fid, rhvar, 'units', 'fraction')
2501 0 : call wrap_put_att_text(fid, lwvar, 'units', 'm')
2502 0 : call wrap_put_att_text(fid, swvar, 'units', 'm')
2503 :
2504 0 : call wrap_put_att_text(fid, rhvar, 'long_name', 'relative humidity')
2505 0 : call wrap_put_att_text(fid, lwvar, 'long_name', 'longwave bands')
2506 0 : call wrap_put_att_text(fid, swvar, 'long_name', 'shortwave bands')
2507 :
2508 : ! Define the variables: abs_lw, ext_sw, ssa_sw, asm_sw
2509 0 : dimids(1) = rhdim
2510 0 : dimids(2) = lwdim
2511 0 : call wrap_def_var(fid, 'abs_lw', NF90_DOUBLE, 2, dimids, abs_lw_var)
2512 :
2513 0 : write(LUNOPRT,*) "Defined abs_lw."
2514 :
2515 0 : call wrap_put_att_text(fid, abs_lw_var, 'units', 'meter^2 kilogram^-1')
2516 :
2517 0 : dimids(1) = rhdim
2518 0 : dimids(2) = swdim
2519 0 : call wrap_def_var(fid, 'ext_sw', NF90_DOUBLE, 2, dimids, ext_sw_var)
2520 0 : call wrap_def_var(fid, 'ssa_sw', NF90_DOUBLE, 2, dimids, ssa_sw_var)
2521 0 : call wrap_def_var(fid, 'asm_sw', NF90_DOUBLE, 2, dimids, asm_sw_var)
2522 :
2523 0 : write(LUNOPRT,*) "Defined ext_sw, ssa_sw, and asm_sw."
2524 :
2525 0 : call wrap_put_att_text(fid, ssa_sw_var, 'units', 'fraction')
2526 0 : call wrap_put_att_text(fid, ext_sw_var, 'units', 'meter^2 kilogram^-1')
2527 0 : call wrap_put_att_text(fid, asm_sw_var, 'units', '-')
2528 :
2529 : ! Define the variables for the refractive indicies.
2530 0 : dimids(1) = swdim
2531 0 : call wrap_def_var(fid, 'refindex_real_aer_sw', NF90_DOUBLE, 1, dimids(1:1), sw_r_refidx_var)
2532 0 : call wrap_def_var(fid, 'refindex_im_aer_sw', NF90_DOUBLE, 1, dimids(1:1), sw_i_refidx_var)
2533 :
2534 0 : write(LUNOPRT,*) "Defined lw refindex."
2535 :
2536 0 : dimids(1) = lwdim
2537 0 : call wrap_def_var(fid, 'refindex_real_aer_lw', NF90_DOUBLE, 1, dimids(1:1), lw_r_refidx_var)
2538 0 : call wrap_def_var(fid, 'refindex_im_aer_lw', NF90_DOUBLE, 1, dimids(1:1), lw_i_refidx_var)
2539 :
2540 0 : write(LUNOPRT,*) "Defined sw refindex."
2541 :
2542 0 : call wrap_put_att_text(fid, sw_r_refidx_var, 'units', '-')
2543 0 : call wrap_put_att_text(fid, sw_i_refidx_var, 'units', '-')
2544 0 : call wrap_put_att_text(fid, lw_r_refidx_var, 'units', '-')
2545 0 : call wrap_put_att_text(fid, lw_i_refidx_var, 'units', '-')
2546 :
2547 0 : call wrap_put_att_text(fid, sw_r_refidx_var, 'long_name', 'real refractive index of aerosol - shortwave')
2548 0 : call wrap_put_att_text(fid, sw_i_refidx_var, 'long_name', 'imaginary refractive index of aerosol - shortwave')
2549 0 : call wrap_put_att_text(fid, lw_r_refidx_var, 'long_name', 'real refractive index of aerosol - longwave')
2550 0 : call wrap_put_att_text(fid, lw_i_refidx_var, 'long_name', 'imaginary refractive index of aerosol - longwave')
2551 :
2552 :
2553 : ! Define fields that define the aerosol properties.
2554 0 : call wrap_def_dim(fid, 'opticsmethod_len', 32, omdim)
2555 0 : dimids(1) = omdim
2556 0 : call wrap_def_var(fid, 'opticsmethod', NF90_CHAR, 1, dimids(1:1), omvar)
2557 :
2558 0 : write(LUNOPRT,*) "Defined omdim."
2559 :
2560 0 : call wrap_def_dim(fid, 'namelength', 20, andim)
2561 0 : dimids(1) = andim
2562 0 : call wrap_def_var(fid, 'aername', NF90_CHAR, 1, dimids(1:1), anvar)
2563 :
2564 0 : write(LUNOPRT,*) "Defined aername."
2565 :
2566 0 : call wrap_def_dim(fid, 'name_len', 32, namedim)
2567 0 : dimids(1) = namedim
2568 0 : call wrap_def_var(fid, 'name', NF90_CHAR, 1, dimids(1:1), namevar)
2569 :
2570 0 : write(LUNOPRT,*) "Defined name."
2571 :
2572 0 : call wrap_def_var(fid, 'density', NF90_DOUBLE, 0, dimids(1:0), denvar)
2573 0 : call wrap_def_var(fid, 'sigma_logr', NF90_DOUBLE, 0, dimids(1:0), slogvar)
2574 0 : call wrap_def_var(fid, 'dryrad', NF90_DOUBLE, 0, dimids(1:0), dryrvar)
2575 0 : call wrap_def_var(fid, 'radmin_aer', NF90_DOUBLE, 0, dimids(1:0), rminvar)
2576 0 : call wrap_def_var(fid, 'radmax_aer', NF90_DOUBLE, 0, dimids(1:0), rmaxvar)
2577 0 : call wrap_def_var(fid, 'hygroscopicity', NF90_DOUBLE, 0, dimids(1:0), hygrovar)
2578 0 : call wrap_def_var(fid, 'num_to_mass_ratio', NF90_DOUBLE, 0, dimids(1:0), ntmvar)
2579 :
2580 0 : call wrap_put_att_text(fid, denvar, 'units', 'kg m^-3')
2581 0 : call wrap_put_att_text(fid, slogvar, 'units', '-')
2582 0 : call wrap_put_att_text(fid, dryrvar, 'units', 'm')
2583 0 : call wrap_put_att_text(fid, rminvar, 'units', 'm')
2584 0 : call wrap_put_att_text(fid, rmaxvar, 'units', 'm')
2585 0 : call wrap_put_att_text(fid, hygrovar, 'units', '-')
2586 0 : call wrap_put_att_text(fid, ntmvar, 'units', 'kg^-1')
2587 :
2588 0 : call wrap_put_att_text(fid, denvar, 'long_name', 'aerosol material density')
2589 0 : call wrap_put_att_text(fid, slogvar, 'long_name', 'geometric standard deviation of aerosol')
2590 0 : call wrap_put_att_text(fid, dryrvar, 'long_name', 'dry number mode radius of aerosol')
2591 0 : call wrap_put_att_text(fid, rminvar, 'long_name', 'minimum dry radius of aerosol for bin')
2592 0 : call wrap_put_att_text(fid, rmaxvar, 'long_name', 'maximum dry radius of aerosol for bin')
2593 0 : call wrap_put_att_text(fid, hygrovar, 'long_name', 'hygroscopicity of aerosol')
2594 0 : call wrap_put_att_text(fid, ntmvar, 'long_name', 'ratio of number to mass of aerosol')
2595 :
2596 :
2597 0 : write(LUNOPRT,*) "Defined all variables."
2598 :
2599 : ! End the defintion phase of the netcdf file.
2600 0 : call wrap_enddef(fid)
2601 :
2602 : ! Write out the dimensions.
2603 0 : call wrap_put_var_realx(fid, rhvar, mie_rh(:nrh))
2604 0 : call wrap_put_var_realx(fid, lwvar, wave(:nlwbands) * 1e-2_f)
2605 0 : call wrap_put_var_realx(fid, swvar, wave(nlwbands+1:) * 1e-2_f)
2606 :
2607 : ! Write out the refractive indicies.
2608 0 : call wrap_put_var_realx(fid, sw_r_refidx_var, real(refidx(nlwbands+1:, 1)))
2609 0 : call wrap_put_var_realx(fid, sw_i_refidx_var, aimag(refidx(nlwbands+1:, 1)))
2610 0 : call wrap_put_var_realx(fid, lw_r_refidx_var, real(refidx(:nlwbands, 1)))
2611 0 : call wrap_put_var_realx(fid, lw_i_refidx_var, aimag(refidx(:nlwbands, 1)))
2612 :
2613 :
2614 : ! Pad the names out with spaces.
2615 0 : aer_name = ' '
2616 0 : aer_name(1:len(trim(c_name))) = c_name
2617 :
2618 0 : start_text(1) = 1
2619 0 : count_text(1) = 32
2620 0 : call wrap_put_vara_text(fid, namevar, start_text, count_text, (/ aer_name /))
2621 0 : count_text(1) = 20
2622 0 : call wrap_put_vara_text(fid, anvar, start_text, count_text, (/ aer_name /))
2623 :
2624 : ! These fields control whether the particle is treated as a CCN. For now,
2625 : ! set these so that CARMA particles are not considered as CCN by the
2626 : ! CAM microphysics.
2627 0 : if (irhswell /= 0) then
2628 0 : count_text(1) = len('hygroscopic ')
2629 0 : call wrap_put_vara_text(fid, omvar, start_text, count_text, (/ 'hygroscopic ' /))
2630 : else
2631 0 : count_text(1) = len('insoluble ')
2632 0 : call wrap_put_vara_text(fid, omvar, start_text, count_text, (/ 'insoluble ' /))
2633 : end if
2634 :
2635 0 : call wrap_put_var_realx(fid, denvar, (/ rho(ibin) * 1e-3_f / 1e-6_f /))
2636 0 : call wrap_put_var_realx(fid, slogvar, (/ 0._f /))
2637 0 : call wrap_put_var_realx(fid, dryrvar, (/ r(ibin) * 1e-2_f /))
2638 0 : call wrap_put_var_realx(fid, rminvar, (/ rlow(ibin) * 1e-2_f /))
2639 0 : call wrap_put_var_realx(fid, rmaxvar, (/ rup(ibin) * 1e-2_f /))
2640 0 : call wrap_put_var_realx(fid, hygrovar, (/ 0._f /))
2641 0 : call wrap_put_var_realx(fid, ntmvar, (/ 1._f / rmass(ibin) / 1e-3_f /))
2642 :
2643 : ! Iterate over a range of relative humidities, since the particle may swell
2644 : ! with relative humidity which will change its optical properties.
2645 0 : do irh = 1, nrh
2646 :
2647 : ! Determine the wet radius.
2648 0 : call getwetr(carma, igroup, mie_rh(irh), r(ibin), rwet, rho(ibin), rhopwet, rc)
2649 0 : if (rc < 0) call endrun('carma_CreateOpticsFile::wetr failed.')
2650 :
2651 : ! Calculate at each wavelength.
2652 0 : do iwave = 1, NWAVE
2653 :
2654 : ! Using Mie code, calculate the optical properties: extinction coefficient,
2655 : ! single scattering albedo and asymmetry factor.
2656 : ! Assume the particle is homogeneous (no core).
2657 : !
2658 : ! NOTE: nmon, df, rmon and falpha are only used for fractal particles.
2659 : call mie(carma, &
2660 0 : carma%f_group(igroup)%f_imiertn, &
2661 : rwet, &
2662 0 : carma%f_wave(iwave), &
2663 0 : real(carma%f_group(igroup)%f_nmon(ibin),kind=f), &
2664 0 : carma%f_group(igroup)%f_df(ibin), &
2665 : carma%f_group(igroup)%f_rmon, &
2666 : carma%f_group(igroup)%f_falpha, &
2667 0 : refidx(iwave, 1), &
2668 : 0.0_f, &
2669 : refidx(iwave, 1), &
2670 : Qext, &
2671 : Qsca, &
2672 : asym, &
2673 0 : rc)
2674 0 : if (rc < 0) call endrun('carma_CreateOpticsFile::mie failed.')
2675 :
2676 : ! Calculate the shortwave and longwave properties?
2677 : !
2678 : ! NOTE: miess is in cgs units, but the optics file needs to be in mks
2679 : ! units, so perform the necessary conversions.
2680 0 : if (iwave <= nlwbands) then
2681 :
2682 : ! Longwave just needs absorption: abs_lw.
2683 0 : abs_lw(irh, iwave) = (Qext - Qsca) * PI * (rwet * 1e-2_f)**2 / (rmass(ibin) * 1e-3_f)
2684 : else
2685 :
2686 : ! Shortwave needs extinction, single scattering albedo and asymmetry factor:
2687 : ! ext_sw, ssa_sw and asm_sw.
2688 0 : ext_sw(irh, iwave - nlwbands) = Qext * PI * (rwet * 1e-2_f)**2 / (rmass(ibin) * 1e-3_f)
2689 0 : ssa_sw(irh, iwave - nlwbands) = Qsca / Qext
2690 0 : asm_sw(irh, iwave - nlwbands) = asym
2691 : end if
2692 : end do
2693 : end do
2694 :
2695 : ! Write out the longwave fields.
2696 0 : ret = nf90_put_var (fid, abs_lw_var, abs_lw(:nrh, :))
2697 0 : if (ret/=NF90_NOERR) then
2698 0 : write(iulog,*)'CARMA_CreateOpticsFile: error writing varid =', abs_lw_var
2699 0 : call handle_error (ret)
2700 : end if
2701 :
2702 : ! Write out the shortwave fields.
2703 0 : ret = nf90_put_var (fid, ext_sw_var, ext_sw(:nrh, :))
2704 0 : if (ret/=NF90_NOERR) then
2705 0 : write(iulog,*)'CARMA_CreateOpticsFile: error writing varid =', ext_sw_var
2706 0 : call handle_error (ret)
2707 : end if
2708 0 : ret = nf90_put_var (fid, ssa_sw_var, ssa_sw(:nrh, :))
2709 0 : if (ret/=NF90_NOERR) then
2710 0 : write(iulog,*)'CARMA_CreateOpticsFile: error writing varid =', ssa_sw_var
2711 0 : call handle_error (ret)
2712 : end if
2713 0 : ret = nf90_put_var (fid, asm_sw_var, asm_sw(:nrh, :))
2714 0 : if (ret/=NF90_NOERR) then
2715 0 : write(iulog,*)'CARMA_CreateOpticsFile: error writing varid =', asm_sw_var
2716 0 : call handle_error (ret)
2717 : end if
2718 :
2719 : ! Close the file.
2720 0 : call wrap_close(fid)
2721 : end if
2722 : end do
2723 :
2724 0 : return
2725 : end subroutine CARMA_CreateOpticsFile_Fixed
2726 :
2727 :
2728 :
2729 : !! This routine creates files containing optical properties for the pure sulfate group
2730 : !! following Yu et al. (2015). These optical properties are used by the RRTMG radiation
2731 : !! code to include the impact of CARMA particles in the radiative transfer
2732 : !! calculation.
2733 0 : subroutine CARMA_CreateOpticsFile_Sulfate(carma, igroup, rc)
2734 : use wrap_nf
2735 : use wetr, only : getwetr
2736 :
2737 : implicit none
2738 :
2739 : type(carma_type), intent(inout) :: carma !! the carma object
2740 : integer, intent(in) :: igroup !! group index
2741 : integer, intent(out) :: rc !! return code, negative indicates failure
2742 :
2743 : ! Local variables
2744 : integer :: ibin, iwave, iwtp
2745 : integer :: irhswell
2746 : integer :: imiertn
2747 : integer :: ienconc
2748 : real(kind=f) :: rho(NBIN), rhopwet
2749 : real(kind=f) :: r(NBIN), rmass(NBIN), rlow(NBIN), rup(NBIN)
2750 : real(kind=f) :: wave(NWAVE)
2751 : complex(kind=f) :: refidx(NWAVE)
2752 : complex(kind=f) :: refidxS(NWAVE, NREFIDX)
2753 : complex(kind=f) :: refidxW(NWAVE)
2754 : character(len=CARMA_NAME_LEN) :: name
2755 : character(len=CARMA_SHORT_NAME_LEN) :: shortname
2756 : integer :: fid
2757 : integer :: rhdim, lwdim, swdim, wtpdim
2758 : integer :: rhvar, lwvar, swvar, wtp_var
2759 : integer :: rwetvar
2760 : integer :: abs_lw_wtp_var, qabs_lw_wtp_var
2761 : integer :: ext_sw_wtp_var, ssa_sw_wtp_var, asm_sw_wtp_var, qext_sw_wtp_var
2762 : integer :: omdim, andim, namedim
2763 : integer :: omvar, anvar, namevar
2764 : integer :: dimids(2)
2765 : integer :: denvar, slogvar, dryrvar, rminvar, rmaxvar, hygrovar, ntmvar
2766 : real(kind=f) :: abs_lw_wtp(NMIE_WTP, nlwbands)
2767 : real(kind=f) :: qabs_lw_wtp(NMIE_WTP, nlwbands)
2768 : real(kind=f) :: ext_sw_wtp(NMIE_WTP, nswbands)
2769 : real(kind=f) :: qext_sw_wtp(NMIE_WTP, nswbands)
2770 : real(kind=f) :: ssa_sw_wtp(NMIE_WTP, nswbands)
2771 : real(kind=f) :: asm_sw_wtp(NMIE_WTP, nswbands)
2772 : character(len=8) :: c_name ! constituent name
2773 : character(len=32) :: aer_name ! long enough for both aername and name
2774 : character(len=255) :: filepath
2775 : real(kind=f) :: rwet
2776 : real(kind=f) :: Qext
2777 : real(kind=f) :: Qsca
2778 : real(kind=f) :: asym
2779 : integer :: start_text(2), count_text(2)
2780 : integer :: sw_r_refidx_var, sw_i_refidx_var, lw_r_refidx_var, lw_i_refidx_var
2781 : integer :: cnsttype ! constituent type
2782 : integer :: maxbin ! last prognostic bin
2783 : integer :: LUNOPRT ! logical unit number for output
2784 : logical :: do_print ! do print output?
2785 : integer :: ret
2786 : real(kind=f) :: volwater
2787 : real(kind=f) :: volsulfate
2788 : real(kind=f) :: volshell
2789 : integer :: igash2o
2790 :
2791 :
2792 : ! Assume success.
2793 0 : rc = 0
2794 :
2795 : ! Get the wavelength structure.
2796 0 : call CARMA_GET(carma, rc, wave=wave, do_print=do_print, LUNOPRT=LUNOPRT, igash2o=igash2o)
2797 0 : if (rc < 0) call endrun('carma_CreateOpticsFile::CARMA_Get failed.')
2798 :
2799 : ! Get the necessary group properties.
2800 : call CARMAGROUP_Get(carma, igroup, rc, name=name, shortname=shortname, r=r, &
2801 : rlow=rlow, rup=rup, rmass=rmass, irhswell=irhswell, &
2802 0 : ienconc=ienconc, cnsttype=cnsttype, maxbin=maxbin, imiertn=imiertn)
2803 0 : if (rc < 0) call endrun('carma_CreateOpticsFile::CARMAGROUP_Get failed.')
2804 :
2805 : ! Get the necessary element properties.
2806 0 : call CARMAELEMENT_Get(carma, ienconc, rc, rho=rho, refidx=refidxS)
2807 0 : if (rc < 0) call endrun('carma_CreateOpticsFile::CARMAELEMENT_Get failed.')
2808 :
2809 : ! Get the refractive index for water.
2810 0 : call CARMAGAS_Get(carma, igash2o, rc, refidx=refidxW)
2811 0 : if (rc < 0) call endrun('carma_CreateOpticsFile::CARMAGAS_Get failed.')
2812 :
2813 : ! A file needs to be created for each bin.
2814 0 : do ibin = 1, NBIN
2815 :
2816 : ! Bins past maxbin are treated as diagnostic even if the group
2817 : ! is prognostic and thus are not advected in the paerent model.
2818 0 : if (ibin <= maxbin) then
2819 :
2820 0 : write(c_name, '(A, I2.2)') trim(shortname), ibin
2821 :
2822 : ! Construct the path to the file. Each model will have its own subdirectory
2823 : ! where the optical property files are stored.
2824 0 : filepath = trim(carma_model) // '_' // trim(c_name) // '_rrtmg.nc'
2825 :
2826 0 : if (do_print) write(LUNOPRT,*) 'Creating CARMA optics file ... ', trim(filepath)
2827 :
2828 : ! Create the file.
2829 0 : call wrap_create(filepath, NF90_CLOBBER, fid)
2830 :
2831 : ! Define the dimensions: rh, lwbands, swbands
2832 0 : call wrap_def_dim(fid, 'rh_idx', NMIE_RH, rhdim)
2833 0 : call wrap_def_dim(fid, 'lw_band', nlwbands, lwdim)
2834 0 : call wrap_def_dim(fid, 'sw_band', nswbands, swdim)
2835 :
2836 0 : call wrap_def_dim(fid, 'wgtpct', NMIE_WTP, wtpdim)
2837 :
2838 0 : dimids(1) = rhdim
2839 0 : call wrap_def_var(fid, 'rh', NF90_DOUBLE, 1, dimids(1), rhvar)
2840 0 : call wrap_def_var(fid, 'rwet',NF90_DOUBLE, 1, dimids(1), rwetvar)
2841 :
2842 0 : dimids(1) = lwdim
2843 0 : call wrap_def_var(fid, 'lw_band', NF90_DOUBLE, 1, dimids(1), lwvar)
2844 :
2845 0 : dimids(1) = swdim
2846 0 : call wrap_def_var(fid, 'sw_band', NF90_DOUBLE, 1, dimids(1), swvar)
2847 :
2848 0 : dimids(1) = wtpdim
2849 0 : call wrap_def_var(fid, 'wgtpct', NF90_DOUBLE, 1, dimids(1), wtp_var)
2850 :
2851 0 : call wrap_put_att_text(fid, rhvar, 'units', 'fraction')
2852 0 : call wrap_put_att_text(fid, rwetvar, 'units', 'cm')
2853 0 : call wrap_put_att_text(fid, lwvar, 'units', 'm')
2854 0 : call wrap_put_att_text(fid, swvar, 'units', 'm')
2855 :
2856 0 : call wrap_put_att_text(fid, wtp_var,'units', 'unitless')
2857 0 : call wrap_put_att_text(fid, wtp_var,'long_name', 'weight percent')
2858 :
2859 0 : call wrap_put_att_text(fid, rhvar, 'long_name', 'relative humidity')
2860 0 : call wrap_put_att_text(fid, rwetvar, 'long_name', 'wet radius')
2861 0 : call wrap_put_att_text(fid, lwvar, 'long_name', 'longwave bands')
2862 0 : call wrap_put_att_text(fid, swvar, 'long_name', 'shortwave bands')
2863 :
2864 : ! Define the variables: abs_lw, ext_sw, ssa_sw, asm_sw
2865 : ! Define 2-dimension (:nrh,:nswbands) LW optics properties: abs_lw, qabs_lw
2866 0 : dimids(1) = wtpdim
2867 0 : dimids(2) = lwdim
2868 0 : call wrap_def_var(fid, 'abs_lw_wtp', NF90_DOUBLE, 2, dimids(1:2), abs_lw_wtp_var)
2869 0 : call wrap_def_var(fid, 'qabs_lw_wtp',NF90_DOUBLE, 2, dimids(1:2), qabs_lw_wtp_var)
2870 :
2871 0 : call wrap_put_att_text(fid, abs_lw_wtp_var, 'units', 'meter^2 kilogram^-1')
2872 0 : call wrap_put_att_text(fid, qabs_lw_wtp_var,'units', '-')
2873 :
2874 : ! Define 2-dimension (:nrh,:nswbands) optics properties: ext_sw, qext_sw, ssa_sw, asm_sw
2875 0 : dimids(1) = wtpdim
2876 0 : dimids(2) = swdim
2877 0 : call wrap_def_var(fid, 'ext_sw_wtp', NF90_DOUBLE, 2, dimids(1:2), ext_sw_wtp_var)
2878 0 : call wrap_def_var(fid, 'qext_sw_wtp',NF90_DOUBLE, 2, dimids(1:2), qext_sw_wtp_var)
2879 0 : call wrap_def_var(fid, 'ssa_sw_wtp', NF90_DOUBLE, 2, dimids(1:2), ssa_sw_wtp_var)
2880 0 : call wrap_def_var(fid, 'asm_sw_wtp', NF90_DOUBLE, 2, dimids(1:2), asm_sw_wtp_var)
2881 :
2882 0 : call wrap_put_att_text(fid, ssa_sw_wtp_var, 'units', 'fraction')
2883 0 : call wrap_put_att_text(fid, qext_sw_wtp_var,'units', '-')
2884 0 : call wrap_put_att_text(fid, ext_sw_wtp_var, 'units', 'meter^2 kilogram^-1')
2885 0 : call wrap_put_att_text(fid, asm_sw_wtp_var, 'units', '-')
2886 :
2887 : ! Define the variables for the refractive indicies.
2888 0 : dimids(1) = swdim
2889 0 : call wrap_def_var(fid, 'refindex_real_aer_sw', NF90_DOUBLE, 1, dimids(1), sw_r_refidx_var)
2890 0 : call wrap_def_var(fid, 'refindex_im_aer_sw', NF90_DOUBLE, 1, dimids(1), sw_i_refidx_var)
2891 :
2892 0 : dimids(1) = lwdim
2893 0 : call wrap_def_var(fid, 'refindex_real_aer_lw', NF90_DOUBLE, 1, dimids(1), lw_r_refidx_var)
2894 0 : call wrap_def_var(fid, 'refindex_im_aer_lw', NF90_DOUBLE, 1, dimids(1), lw_i_refidx_var)
2895 :
2896 0 : call wrap_put_att_text(fid, sw_r_refidx_var, 'units', '-')
2897 0 : call wrap_put_att_text(fid, sw_i_refidx_var, 'units', '-')
2898 0 : call wrap_put_att_text(fid, lw_r_refidx_var, 'units', '-')
2899 0 : call wrap_put_att_text(fid, lw_i_refidx_var, 'units', '-')
2900 :
2901 0 : call wrap_put_att_text(fid, sw_r_refidx_var, 'long_name', 'real refractive index of aerosol - shortwave')
2902 0 : call wrap_put_att_text(fid, sw_i_refidx_var, 'long_name', 'imaginary refractive index of aerosol - shortwave')
2903 0 : call wrap_put_att_text(fid, lw_r_refidx_var, 'long_name', 'real refractive index of aerosol - longwave')
2904 0 : call wrap_put_att_text(fid, lw_i_refidx_var, 'long_name', 'imaginary refractive index of aerosol - longwave')
2905 :
2906 :
2907 : ! Define fields that define the aerosol properties.
2908 0 : call wrap_def_dim(fid, 'opticsmethod_len', 32, omdim)
2909 0 : dimids(1) = omdim
2910 0 : call wrap_def_var(fid, 'opticsmethod', NF90_CHAR, 1, dimids(1), omvar)
2911 :
2912 0 : call wrap_def_dim(fid, 'namelength', 20, andim)
2913 0 : dimids(1) = andim
2914 0 : call wrap_def_var(fid, 'aername', NF90_CHAR, 1, dimids(1), anvar)
2915 :
2916 0 : call wrap_def_dim(fid, 'name_len', 32, namedim)
2917 0 : dimids(1) = namedim
2918 0 : call wrap_def_var(fid, 'name', NF90_CHAR, 1, dimids, namevar)
2919 :
2920 0 : call wrap_def_var(fid, 'density', NF90_DOUBLE, 0, dimids(1), denvar)
2921 0 : call wrap_def_var(fid, 'sigma_logr', NF90_DOUBLE, 0, dimids(1), slogvar)
2922 0 : call wrap_def_var(fid, 'dryrad', NF90_DOUBLE, 0, dimids(1), dryrvar)
2923 0 : call wrap_def_var(fid, 'radmin_aer', NF90_DOUBLE, 0, dimids(1), rminvar)
2924 0 : call wrap_def_var(fid, 'radmax_aer', NF90_DOUBLE, 0, dimids(1), rmaxvar)
2925 0 : call wrap_def_var(fid, 'hygroscopicity', NF90_DOUBLE, 0, dimids(1), hygrovar)
2926 0 : call wrap_def_var(fid, 'num_to_mass_ratio', NF90_DOUBLE, 0, dimids(1), ntmvar)
2927 :
2928 0 : call wrap_put_att_text(fid, denvar, 'units', 'kg m^-3')
2929 0 : call wrap_put_att_text(fid, slogvar, 'units', '-')
2930 0 : call wrap_put_att_text(fid, dryrvar, 'units', 'm')
2931 0 : call wrap_put_att_text(fid, rminvar, 'units', 'm')
2932 0 : call wrap_put_att_text(fid, rmaxvar, 'units', 'm')
2933 0 : call wrap_put_att_text(fid, hygrovar, 'units', '-')
2934 0 : call wrap_put_att_text(fid, ntmvar, 'units', 'kg^-1')
2935 :
2936 0 : call wrap_put_att_text(fid, denvar, 'long_name', 'aerosol material density')
2937 0 : call wrap_put_att_text(fid, slogvar, 'long_name', 'geometric standard deviation of aerosol')
2938 0 : call wrap_put_att_text(fid, dryrvar, 'long_name', 'dry number mode radius of aerosol')
2939 0 : call wrap_put_att_text(fid, rminvar, 'long_name', 'minimum dry radius of aerosol for bin')
2940 0 : call wrap_put_att_text(fid, rmaxvar, 'long_name', 'maximum dry radius of aerosol for bin')
2941 0 : call wrap_put_att_text(fid, hygrovar, 'long_name', 'hygroscopicity of aerosol')
2942 0 : call wrap_put_att_text(fid, ntmvar, 'long_name', 'ratio of number to mass of aerosol')
2943 :
2944 : ! End the defintion phase of the netcdf file.
2945 0 : call wrap_enddef(fid)
2946 :
2947 : ! Write out the dimensions.
2948 0 : call wrap_put_var_realx(fid, rhvar, mie_rh(:))
2949 0 : call wrap_put_var_realx(fid, lwvar, wave(:nlwbands) * 1e-2_f)
2950 0 : call wrap_put_var_realx(fid, swvar, wave(nlwbands+1:) * 1e-2_f)
2951 :
2952 0 : call wrap_put_var_realx(fid, wtp_var, mie_wtp(:)*100._f)
2953 :
2954 : ! Write out the refractive indicies.
2955 0 : call wrap_put_var_realx(fid, sw_r_refidx_var, real(refidxS(nlwbands+1:, 1)))
2956 0 : call wrap_put_var_realx(fid, sw_i_refidx_var, aimag(refidxS(nlwbands+1:, 1)))
2957 0 : call wrap_put_var_realx(fid, lw_r_refidx_var, real(refidxS(:nlwbands, 1)))
2958 0 : call wrap_put_var_realx(fid, lw_i_refidx_var, aimag(refidxS(:nlwbands, 1)))
2959 :
2960 : ! Pad the names out with spaces.
2961 0 : aer_name = ' '
2962 0 : aer_name(1:len(trim(c_name))) = c_name
2963 :
2964 0 : start_text(1) = 1
2965 0 : count_text(1) = 32
2966 0 : call wrap_put_vara_text(fid, namevar, start_text, count_text, (/ aer_name /))
2967 0 : count_text(1) = 20
2968 0 : call wrap_put_vara_text(fid, anvar, start_text, count_text, (/ aer_name /))
2969 :
2970 0 : count_text(1) = len('hygroscopic_wtp ')
2971 0 : call wrap_put_vara_text(fid, omvar, start_text, count_text, (/ 'hygroscopic_wtp ' /))
2972 :
2973 0 : call wrap_put_var_realx(fid, denvar, (/ rho(ibin) * 1e-3_f / 1e-6_f /))
2974 0 : call wrap_put_var_realx(fid, slogvar, (/ 0._f /))
2975 0 : call wrap_put_var_realx(fid, dryrvar, (/ r(ibin) * 1e-2_f /))
2976 0 : call wrap_put_var_realx(fid, rminvar, (/ rlow(ibin) * 1e-2_f /))
2977 0 : call wrap_put_var_realx(fid, rmaxvar, (/ rup(ibin) * 1e-2_f /))
2978 0 : call wrap_put_var_realx(fid, hygrovar, (/ 0.6_f /))
2979 0 : call wrap_put_var_realx(fid, ntmvar, (/ 1._f / rmass(ibin) / 1e-3_f /))
2980 :
2981 : ! For now, ext_sw(:nrh, :nswbands) and ext_sw_coreshell(:nrh, :nswbands, :ncoreshellratio) both are calculated
2982 : ! Since other aerosols in CAM may use ext_sw rather than ext_sw_coreshell
2983 : ! Modified by Pengfei Yu
2984 : ! April.1, 2012
2985 :
2986 : ! calculate qext and ext for pure sulfate dependent on weight percent
2987 : ! ideally qext is based on (wgt,temp,wave), however Beyer et al. (1996) Figure 5
2988 : ! shows sulfate density is roughly 0.006 g/cm3/k, I negelet temp dimension, assuming temp = 270 K
2989 : ! In code, sulfate density is precisely calculated to determine wet raidus
2990 0 : do iwtp = 1, NMIE_WTP
2991 :
2992 : ! NOTE: Weight percent is normal a result of the getwetr calculation. To build the
2993 : ! table based upon weight percent, we need to pass in the desired value and a
2994 : ! reference temperature. In that case, the RH is ignored.
2995 0 : call getwetr(carma, igroup, mie_rh(1), r(ibin), rwet, rho(ibin), rhopwet, rc, wgtpct=mie_wtp(iwtp)*100._f, temp=270._f)
2996 0 : if (rc < 0) call endrun('carma_CreateOpticsFile::wetr failed.')
2997 :
2998 : ! This is not in Yu (2015), but rather than using the refractive
2999 : ! index of H2SO4 for the shell, do a volume mix of water and H2SO4
3000 : ! for the refractive index of the shell.
3001 0 : volwater = rwet**3._f - r(ibin)**3._f
3002 0 : volsulfate = r(ibin)**3._f
3003 0 : volshell = volwater + volsulfate
3004 0 : if (volshell > 0._f) then
3005 0 : refidx(:) = (volwater / volshell) * refidxW(:) + (volsulfate / volshell) * refidxS(:, 1)
3006 : else
3007 0 : refidx(:) = refidxS(:, 1)
3008 : end if
3009 :
3010 : ! Calculate at each wavelength.
3011 0 : do iwave = 1, NWAVE
3012 :
3013 : ! Using Mie code, calculate the optical properties: extinction coefficient,
3014 : ! single scattering albedo and asymmetry factor.
3015 : ! Assume the particle is homogeneous (no core).
3016 : !
3017 : ! NOTE: The refractive index for sulfate changes with RH/weight percent, which
3018 : ! is not reflected in this code.
3019 : call mie(carma, &
3020 : imiertn, &
3021 : rwet, &
3022 0 : wave(iwave), &
3023 : 0._f, &
3024 : 3.0_f, &
3025 : 0.0_f, &
3026 : 1.0_f, &
3027 : refidx(iwave), &
3028 : 0.0_f, &
3029 : refidx(iwave), &
3030 : Qext, &
3031 : Qsca, &
3032 : asym, &
3033 0 : rc)
3034 0 : if (rc < 0) call endrun('carma_CreateOpticsFile::mie failed.')
3035 :
3036 : ! Calculate the shortwave and longwave properties?
3037 : !
3038 : ! NOTE: miess is in cgs units, but the optics file needs to be in mks
3039 : ! units, so perform the necessary conversions.
3040 0 : if (iwave <= nlwbands) then
3041 :
3042 : ! Longwave just needs absorption: abs_lw.
3043 0 : qabs_lw_wtp(iwtp, iwave) = (Qext - Qsca) ! absorption per particle
3044 0 : abs_lw_wtp (iwtp, iwave) = (Qext - Qsca) * PI * (rwet * 1e-2_f)**2 / (rmass(ibin) * 1e-3_f)
3045 : else
3046 :
3047 : ! Shortwave needs extinction, single scattering albedo and asymmetry factor:
3048 : ! ext_sw, ssa_sw and asm_sw.
3049 0 : qext_sw_wtp(iwtp, iwave - nlwbands) = Qext ! extinction per particle
3050 0 : ext_sw_wtp (iwtp, iwave - nlwbands) = Qext * PI * (rwet * 1e-2_f)**2 / (rmass(ibin) * 1e-3_f)
3051 0 : ssa_sw_wtp (iwtp, iwave - nlwbands) = Qsca / Qext
3052 0 : asm_sw_wtp (iwtp, iwave - nlwbands) = asym
3053 : end if
3054 : end do ! iwave
3055 : end do ! iwtp
3056 :
3057 : ! Write out the longwave fields.
3058 0 : ret = nf90_put_var(fid, abs_lw_wtp_var, abs_lw_wtp (:, :))
3059 0 : if (ret /= NF90_NOERR) then
3060 0 : write(iulog,*)'CARMA_CreateOpticsFile_SulfateYu: error writing varid =', fid, abs_lw_wtp_var
3061 0 : call handle_error(ret)
3062 : end if
3063 :
3064 0 : ret = nf90_put_var(fid, qabs_lw_wtp_var, qabs_lw_wtp(:, :))
3065 0 : if (ret /= NF90_NOERR) then
3066 0 : write(iulog,*)'CARMA_CreateOpticsFile_SulfateYu: error writing varid =', qabs_lw_wtp_var
3067 0 : call handle_error(ret)
3068 : end if
3069 :
3070 : ! Write out the shortwave fields.
3071 0 : ret = nf90_put_var(fid, ext_sw_wtp_var, ext_sw_wtp (:, :))
3072 0 : if (ret /= NF90_NOERR) then
3073 0 : write(iulog,*)'CARMA_CreateOpticsFile_SulfateYu: error writing varid =', ext_sw_wtp_var
3074 0 : call handle_error(ret)
3075 : end if
3076 :
3077 0 : ret = nf90_put_var(fid, qext_sw_wtp_var,qext_sw_wtp(:, :))
3078 0 : if (ret /= NF90_NOERR) then
3079 0 : write(iulog,*)'CARMA_CreateOpticsFile_SulfateYu: error writing varid =', qext_sw_wtp_var
3080 0 : call handle_error(ret)
3081 : end if
3082 :
3083 0 : ret = nf90_put_var(fid, ssa_sw_wtp_var, ssa_sw_wtp (:, :))
3084 0 : if (ret /= NF90_NOERR) then
3085 0 : write(iulog,*)'CARMA_CreateOpticsFile_SulfateYu: error writing varid =', ssa_sw_wtp_var
3086 0 : call handle_error(ret)
3087 : end if
3088 :
3089 0 : ret = nf90_put_var(fid, asm_sw_wtp_var, asm_sw_wtp (:, :))
3090 0 : if (ret /= NF90_NOERR) then
3091 0 : write(iulog,*)'CARMA_CreateOpticsFile_SulfateYu: error writing varid =', asm_sw_wtp_var
3092 0 : call handle_error(ret)
3093 : end if
3094 :
3095 : ! Close the file.
3096 0 : call wrap_close(fid)
3097 : end if
3098 : end do
3099 :
3100 0 : return
3101 : end subroutine CARMA_CreateOpticsFile_Sulfate
3102 :
3103 :
3104 : !! Calculate the aerodynamic resistance for dry deposition.
3105 : !!
3106 : !! This is based upon Seinfeld and Pandis (1998) page 963, and
3107 : !! is similar to the calcram routine in drydep_mod.F90;
3108 : !! however, it enables independent determination of the aerodynamic
3109 : !! resistance each surface type (land, ocean and ice).
3110 : !!
3111 : !! @author Tianyi Fan
3112 : !! @version Aug 2011
3113 852463 : subroutine CARMA_calcram(ustar, z0, pdel, pmid, tmid, obklen, ram)
3114 : use shr_const_mod, only: shr_const_karman
3115 : use physconst, only: rair, gravit
3116 :
3117 : implicit none
3118 :
3119 : ! input and output argument
3120 : real(r8), intent(in) :: ustar ! friction velocity
3121 : real(r8), intent(in) :: z0 ! roughness length
3122 : real(r8), intent(in) :: pdel ! layer depth in pressure [Pa]
3123 : real(r8), intent(in) :: pmid ! layer mid-point pressure [Pa]
3124 : real(r8), intent(in) :: tmid ! layer mid-point temperature [K]
3125 : real(r8), intent(in) :: obklen ! Monin-Obukhov length
3126 : real(r8), intent(out) :: ram ! aerodynamic resistance
3127 :
3128 : ! local varibles
3129 : real(r8) :: z ! half the layer height
3130 : real(r8) :: psi ! stability parameter for z
3131 : real(r8) :: psi0 ! stability parameter for z0
3132 : real(r8) :: nu ! temparory variable
3133 : real(r8) :: nu0 ! temparory variable
3134 : real(r8), parameter :: xkar = shr_const_karman
3135 :
3136 :
3137 : ! Use half the layer height like Ganzefeld and Lelieveld, 1995
3138 852463 : z = pdel * rair * tmid / pmid / gravit / 2._r8
3139 :
3140 852463 : if (obklen .eq. 0._r8) then
3141 : psi = 0._r8
3142 : psi0 = 0._r8
3143 : else
3144 852463 : psi = min(max(z / obklen, -1._r8), 1._r8)
3145 852463 : psi0 = min(max(z0 / obklen, -1._r8), 1._r8)
3146 : endif
3147 :
3148 : ! Stable
3149 852463 : if (psi > 0._r8) then
3150 172628 : ram = 1._r8 / xkar / ustar * (log(z / z0) + 4.7_r8 * (psi - psi0))
3151 :
3152 : ! Unstable
3153 679835 : else if (psi < 0._r8) then
3154 679835 : nu = (1._r8 - 15._r8 *psi)**(.25_r8)
3155 679835 : nu0 = (1._r8 - 15._r8 *psi0)**(.25_r8)
3156 :
3157 679835 : if (ustar /= 0._r8) then
3158 : ram = 1._r8 / xkar / ustar * (log(z / z0) + &
3159 : log(((nu0**2 + 1._r8) * (nu0 + 1._r8)**2) / &
3160 : ((nu**2 + 1._r8) * (nu + 1._r8)**2)) + &
3161 679835 : 2._r8 * (atan(nu) - atan(nu0)))
3162 : else
3163 0 : ram = 0._r8
3164 : end if
3165 :
3166 : ! Neutral
3167 : else
3168 0 : ram = 1._r8 / xkar / ustar * log(z / z0)
3169 : end if
3170 :
3171 852463 : return
3172 : end subroutine CARMA_calcram
3173 :
3174 : !---------------------------------------------------------------------------
3175 : ! define fields for reference profiles in cam restart file
3176 : !---------------------------------------------------------------------------
3177 1536 : subroutine CARMA_restart_init( File )
3178 : use cam_pio_utils, only: cam_pio_def_dim
3179 : use pio, only: file_desc_t, pio_def_var, pio_double
3180 :
3181 : ! arguments
3182 : type(file_desc_t),intent(inout) :: File ! pio File pointer
3183 :
3184 : ! local variables
3185 : integer :: levid, ierr
3186 :
3187 1536 : if (carma_do_fixedinit) then
3188 0 : call cam_pio_def_dim(File, 'lev', pver, levid, existOK=.true.)
3189 0 : ierr = pio_def_var(File, 'CARMA_REF_T', pio_double, (/ levid /), t_ref_desc)
3190 0 : ierr = pio_def_var(File, 'CARMA_REF_H2O', pio_double, (/ levid /), h2o_ref_desc)
3191 0 : ierr = pio_def_var(File, 'CARMA_REF_H2SO4', pio_double, (/ levid /), h2so4_ref_desc)
3192 : endif
3193 :
3194 1536 : end subroutine CARMA_restart_init
3195 :
3196 : !---------------------------------------------------------------------------
3197 : ! write reference profiles to restart file
3198 : !---------------------------------------------------------------------------
3199 1536 : subroutine CARMA_restart_write(File)
3200 1536 : use pio, only: file_desc_t, pio_put_var
3201 :
3202 : ! arguments
3203 : type(file_desc_t), intent(inout) :: File
3204 :
3205 : ! local variables
3206 : integer ::ierr
3207 :
3208 1536 : if (carma_do_fixedinit) then
3209 0 : ierr = pio_put_var(File, t_ref_desc, carma_t_ref)
3210 0 : if (carma%f_igash2o /= 0) then
3211 0 : ierr = pio_put_var(File, h2o_ref_desc, carma_h2o_ref)
3212 : endif
3213 0 : if (carma%f_igash2So4 /= 0) then
3214 0 : ierr = pio_put_var(File, h2so4_ref_desc, carma_h2so4_ref)
3215 : endif
3216 : endif
3217 :
3218 1536 : end subroutine CARMA_restart_write
3219 :
3220 : !---------------------------------------------------------------------------
3221 : ! read reference profiles from restart file
3222 : !---------------------------------------------------------------------------
3223 768 : subroutine CARMA_restart_read(File)
3224 : use pio, only: file_desc_t, pio_inq_varid, pio_get_var
3225 :
3226 : ! arguments
3227 : type(file_desc_t),intent(inout) :: File ! pio File pointer
3228 :
3229 : ! local variables
3230 : integer :: ierr, varid
3231 : character(len=*), parameter :: subname = 'CARMA_restart_read: '
3232 :
3233 768 : if (carma_do_fixedinit) then
3234 0 : ierr = pio_inq_varid(File, 'CARMA_REF_T', varid)
3235 0 : if (varid>0) then
3236 0 : ierr = pio_get_var(File, varid, carma_t_ref)
3237 : else
3238 0 : call endrun(subname//'restart file must include CARMA_REF_T')
3239 : endif
3240 0 : ierr = pio_inq_varid(File, 'CARMA_REF_H2O', varid)
3241 0 : if (varid>0) then
3242 0 : ierr = pio_get_var(File, varid, carma_h2o_ref)
3243 0 : else if (carma%f_igash2o /= 0) then
3244 0 : call endrun(subname//'restart file must include CARMA_REF_H2O')
3245 : endif
3246 0 : ierr = pio_inq_varid(File, 'CARMA_REF_H2SO4', varid)
3247 0 : if (varid>0) then
3248 0 : ierr = pio_get_var(File, varid, carma_h2so4_ref)
3249 0 : else if (carma%f_igash2So4 /= 0) then
3250 0 : call endrun(subname//'restart file must include CARMA_REF_H2SO4')
3251 : endif
3252 : endif
3253 :
3254 768 : end subroutine CARMA_restart_read
3255 :
3256 : !! Get the mixing ratio for the specified element and bin.
3257 : !!
3258 : !! @author Chuck Bardeen
3259 : !! @version Aug 2023
3260 95591986194 : subroutine carma_get_bin(state, ielem, ibin, mmr, rc)
3261 : type(physics_state), intent(in) :: state !! physics state variables
3262 : integer, intent(in) :: ielem !! element index
3263 : integer, intent(in) :: ibin !! bin index
3264 : real(r8), intent(out) :: mmr(pcols,pver) !! mass mixing ratio (kg/kg)
3265 : integer, intent(out) :: rc !! return code
3266 :
3267 : integer :: ncol
3268 :
3269 : ! default return code
3270 95591986194 : rc = RC_OK
3271 :
3272 95591986194 : ncol = state%ncol
3273 :
3274 : ! Check the group and bin ranges
3275 95591986194 : if ((ielem < 1) .or. (ielem .gt. NELEM)) then
3276 0 : write(LUNOPRT, *) 'carma_get_bin:: ERROR - Invalid element id, ', ielem
3277 0 : rc = RC_ERROR
3278 0 : return
3279 : end if
3280 :
3281 95591986194 : if ((ibin < 1) .or. (ibin .gt. NBIN)) then
3282 0 : write(LUNOPRT, *) 'carma_get_bin:: ERROR - Invalid bin id, ', ibin
3283 0 : rc = RC_ERROR
3284 0 : return
3285 : end if
3286 :
3287 : ! Get the element from the physics state
3288 >1032*10^11 : mmr(:ncol, :) = state%q(:ncol, :, icnst4elem(ielem, ibin))
3289 :
3290 : return
3291 : end subroutine
3292 :
3293 : !! Get the mixing ratio for the specified element and bin.
3294 3095857800 : subroutine carma_get_bin_cld(pbuf, ielem, ibin, ncol, nlev, mmr, rc)
3295 : type(physics_buffer_desc), pointer :: pbuf(:) !! physics buffer
3296 : integer, intent(in) :: ielem !! element index
3297 : integer, intent(in) :: ibin !! bin index
3298 : integer, intent(in) :: ncol,nlev !! dimensions
3299 : real(r8), intent(out) :: mmr(:,:) !! mass mixing ratio (kg/kg)
3300 : integer, intent(out) :: rc !! return code
3301 :
3302 3095857800 : real(r8), pointer :: mmr_ptr(:,:)
3303 : character(len=8) :: shortname ! short (CAM) name
3304 : character(len=16) :: c_name
3305 : integer :: idx
3306 :
3307 : ! default return code
3308 3095857800 : rc = RC_OK
3309 :
3310 3095857800 : call CARMAELEMENT_Get(carma, ielem, rc, shortname=shortname)
3311 :
3312 3095857800 : write(c_name, '(A, I2.2)') trim(shortname), ibin
3313 :
3314 3095857800 : idx = pbuf_get_index('CLD'//trim(c_name))
3315 3095857800 : call pbuf_get_field(pbuf, idx, mmr_ptr)
3316 :
3317 >33431*10^8 : mmr(:ncol,:nlev) = mmr_ptr(:ncol,:nlev)
3318 :
3319 3095857800 : end subroutine carma_get_bin_cld
3320 :
3321 : !! Determine the dry radius and dry density for the particular bin.
3322 : !!
3323 : !! @author Chuck Bardeen
3324 : !! @version Aug 2023
3325 160204800 : subroutine carma_get_dry_radius(state, igroup, ibin, rdry, rhopdry, rc)
3326 :
3327 : implicit none
3328 :
3329 : type(physics_state), intent(in) :: state !! physics state variables
3330 : integer, intent(in) :: igroup !! group index
3331 : integer, intent(in) :: ibin !! bin index
3332 : real(r8), intent(out) :: rdry(:,:) !! dry radius (m)
3333 : real(r8), intent(out) :: rhopdry(:,:) !! dry density (kg/m3)
3334 : integer, intent(out) :: rc !! return code
3335 :
3336 : real(r8) :: rhoelem(NBIN) ! element density (g/cm3)
3337 : real(r8) :: totvol(pcols,pver) ! total volume (m3/kg)
3338 : real(r8) :: totmmr(pcols,pver) ! total mmr (kg/kg)
3339 : real(r8) :: mmr(pcols, pver) ! mass mixing ratio (kg/kg)
3340 : real(r8) :: nmr(pcols, pver) ! number mixing ratio (#/kg)
3341 : integer :: nelems ! number of elements in group
3342 : integer :: ielems(NELEM) ! element indexes for group
3343 : integer :: ncol
3344 : integer :: i
3345 : integer :: ielem
3346 :
3347 : ! default return code
3348 160204800 : rc = RC_OK
3349 :
3350 160204800 : ncol = state%ncol
3351 :
3352 : ! Check the group and bin ranges
3353 160204800 : if ((igroup < 1) .or. (igroup .gt. NGROUP)) then
3354 0 : write(LUNOPRT, *) 'carma_get_dry_radius:: ERROR - Invalid group id, ', igroup
3355 0 : rc = RC_ERROR
3356 0 : return
3357 : end if
3358 :
3359 160204800 : if ((ibin < 1) .or. (ibin .gt. NBIN)) then
3360 0 : write(LUNOPRT, *) 'carma_get_dry_radius:: ERROR - Invalid bin id, ', ibin
3361 0 : rc = RC_ERROR
3362 0 : return
3363 : end if
3364 :
3365 : ! Iterate over all of the composition and determine the dry volume and dry radius.
3366 160204800 : call carma_get_elem_for_group(igroup, nelems, ielems, rc)
3367 160204800 : if (rc < 0) return
3368 :
3369 >17286*10^7 : totvol(:ncol, :) = 0._r8
3370 >17286*10^7 : totmmr(:ncol, :) = 0._r8
3371 >17286*10^7 : rhopdry(:ncol, :)= 0._r8
3372 >17286*10^7 : rdry(:ncol, :) = 0._r8
3373 :
3374 962841600 : do i = 1, nelems
3375 802636800 : ielem = ielems(i)
3376 :
3377 802636800 : call CARMAELEMENT_Get(carma, ielem, rc, rho=rhoelem)
3378 802636800 : if (rc < 0) return
3379 :
3380 802636800 : call carma_get_bin(state, ielem, ibin, mmr, rc)
3381 802636800 : if (rc < 0) return
3382 :
3383 >86604*10^7 : totmmr(:ncol, :) = totmmr(:ncol, :) + mmr(:ncol, :)
3384 >86620*10^7 : totvol(:ncol, :) = totvol(:ncol, :) + mmr(:ncol, :) / (rhoelem(ibin) / 1.e3_r8 * 1.e6_r8)
3385 : end do
3386 :
3387 : ! Add checks for totvol = 0 and nmr = 0
3388 >17286*10^7 : where(totvol(:ncol, :)>0._r8)
3389 160204800 : rhopdry(:ncol, :) = totmmr(:ncol, :) / totvol(:ncol, :)
3390 : end where
3391 :
3392 160204800 : call carma_get_number(state, igroup, ibin, nmr, rc)
3393 160204800 : if (rc < 0) return
3394 :
3395 >17286*10^7 : where(nmr(:ncol, :)>0._r8)
3396 160204800 : rdry(:ncol, :) = ((3._r8 * totvol(:ncol, :) / nmr(:ncol, :)) / (4._r8 * PI)) ** (1._r8 / 3._r8)
3397 : !rdry(:ncol, :) = ((three_o_fourpi* totvol(:ncol, :) / nmr(:ncol, :))) ** onethird
3398 : end where
3399 :
3400 : return
3401 : end subroutine carma_get_dry_radius
3402 :
3403 :
3404 : !! Get the number of elements and list of element ids for a group. This includes
3405 : !! the concentration elements and the core masses.
3406 : !!
3407 : !! @author Chuck Bardeen
3408 : !! @version Aug 2023
3409 34477941998 : subroutine carma_get_elem_for_group(igroup, nelems, ielems, rc)
3410 : integer, intent(in) :: igroup !! group index
3411 : integer, intent(out) :: nelems !! number of elements in group
3412 : integer, intent(out) :: ielems(NELEM) !! indexes of elements in group
3413 : integer, intent(out) :: rc !! return code
3414 :
3415 : integer :: ienconc
3416 : integer :: ncore
3417 : integer :: icorelem(NELEM)
3418 :
3419 : ! default return code
3420 17238970999 : rc = RC_OK
3421 :
3422 : ! Check the group range.
3423 17238970999 : if ((igroup < 1) .or. (igroup .gt. NGROUP)) then
3424 0 : write(LUNOPRT, *) 'carma_get_elem_for_group:: ERROR - Invalid group id, ', igroup
3425 0 : rc = RC_ERROR
3426 0 : return
3427 : end if
3428 :
3429 17238970999 : call CARMAGROUP_Get(carma, igroup, rc, ienconc=ienconc, ncore=ncore, icorelem=icorelem)
3430 :
3431 17238970999 : nelems = ncore + 1
3432 17238970999 : ielems(1) = ienconc
3433 :
3434 17238970999 : if (ncore .gt. 0) then
3435 97738647594 : ielems(2:ncore+1) = icorelem(1:ncore)
3436 : end if
3437 :
3438 : return
3439 : end subroutine
3440 :
3441 :
3442 : !! Get the CARMA group id a group name.
3443 : !!
3444 : !! @author Chuck Bardeen
3445 : !! @version Aug 2023
3446 2298266873 : subroutine carma_get_group_by_name(shortname, igroup, rc)
3447 : character(len=*), intent(in) :: shortname !! the group short name
3448 : integer, intent(out) :: igroup !! group index
3449 : integer, intent(out) :: rc !! return code
3450 :
3451 : integer :: i
3452 : character(len=32) :: name
3453 :
3454 : ! default return code
3455 2298266873 : rc = RC_OK
3456 :
3457 2298266873 : igroup = -1
3458 :
3459 : ! Check the short names of each group for one that matches
3460 3535465627 : do i = 1, NGROUP
3461 3535465627 : call CARMAGROUP_Get(carma, i, rc, shortname=name)
3462 :
3463 3535465627 : if (trim(shortname) .eq. trim(name)) then
3464 2298266873 : igroup = i
3465 2298266873 : exit
3466 : end if
3467 : end do
3468 :
3469 2298266873 : if (igroup .eq. -1) then
3470 0 : write(LUNOPRT, *) 'carma_get_group_by_name:: ERROR - group not found, ', shortname
3471 0 : rc = RC_ERROR
3472 0 : return
3473 : end if
3474 :
3475 : return
3476 : end subroutine
3477 :
3478 :
3479 : !! Get the CARMA group id and bin id from a compound name xxxxxxnn, where xxxxxx is the
3480 : !! name of the group and nn is the two digit bin number.
3481 : !!
3482 : !! @author Chuck Bardeen
3483 : !! @version Aug 2023
3484 : subroutine carma_get_group_and_bin_by_name(shortname, igroup, ibin, rc)
3485 : character(len=*), intent(out) :: shortname !! the group short name
3486 : integer, intent(out) :: igroup !! group index
3487 : integer, intent(out) :: ibin !! bin index
3488 : integer, intent(out) :: rc !! return code
3489 :
3490 : integer :: i
3491 : character(len=32) :: name
3492 : character(len=32) :: groupname
3493 : character(len=32) :: binname
3494 :
3495 : ! default return code
3496 : rc = RC_OK
3497 :
3498 : igroup = -1
3499 : ibin = -1
3500 :
3501 : if (len(shortname) <= 2) then
3502 : write(LUNOPRT, *) 'carma_get_group_and_bin_by_name:: ERROR - Illegal shortname, ' // shortname
3503 : rc = RC_ERROR
3504 : return
3505 : end if
3506 :
3507 : ! Check the short names of each group for one that matches
3508 : groupname = shortname(:len(shortname)-2)
3509 : binname = shortname(len(shortname)-2:)
3510 :
3511 : call carma_get_group_by_name(groupname, igroup, rc)
3512 : if (rc < 0) return
3513 :
3514 : read(binname, *) ibin
3515 :
3516 : return
3517 : end subroutine
3518 :
3519 :
3520 : !! Determine a mass weighted kappa for the entire particle.
3521 : !!
3522 : !! @author Chuck Bardeen
3523 : !! @version Aug 2023
3524 15097198999 : subroutine carma_get_kappa(state, igroup, ibin, kappa, rc)
3525 :
3526 : implicit none
3527 :
3528 : type(physics_state), intent(in) :: state !! physics state variables
3529 : integer, intent(in) :: igroup !! group index
3530 : integer, intent(in) :: ibin !! bin index
3531 : real(r8), intent(out) :: kappa(:,:) !! kappa value for the entire particle
3532 : integer, intent(out) :: rc !! return code
3533 :
3534 : real(r8) :: totmmr(pcols,pver) ! total mmr (kg/kg)
3535 : real(r8) :: mmr(pcols,pver) ! element mmr (kg/kg)
3536 : real(r8) :: kappaelem ! element kappa
3537 : integer :: ncol
3538 : integer :: nelems
3539 : integer :: ielems(NELEM)
3540 : integer :: i
3541 : integer :: ielem
3542 :
3543 : ! default return code
3544 15097198999 : rc = RC_OK
3545 :
3546 15097198999 : ncol = state%ncol
3547 :
3548 : ! Check the group and bin ranges
3549 15097198999 : if ((igroup < 1) .or. (igroup .gt. NGROUP)) then
3550 0 : write(LUNOPRT, *) 'carma_get_kappa:: ERROR - Invalid group id, ', igroup
3551 0 : rc = RC_ERROR
3552 0 : return
3553 : end if
3554 :
3555 15097198999 : if ((ibin < 1) .or. (igroup .gt. NBIN)) then
3556 0 : write(LUNOPRT, *) 'carma_get_kappa:: ERROR - Invalid bin id, ', ibin
3557 0 : rc = RC_ERROR
3558 0 : return
3559 : end if
3560 :
3561 : ! Iterate over all of the composition and determine the total mass.
3562 15097198999 : call carma_get_elem_for_group(igroup, nelems, ielems, rc)
3563 15097198999 : if (rc < 0) return
3564 :
3565 >16307*10^9 : totmmr(:ncol, :) = 0._r8
3566 >16307*10^9 : kappa(:ncol, :) = 0._r8
3567 :
3568 >10568*10^7 : do i = 1, nelems
3569 90583193994 : ielem = ielems(i)
3570 :
3571 90583193994 : call carma_get_bin(state, ielem, ibin, mmr, rc)
3572 90583193994 : if (rc < 0) return
3573 :
3574 90583193994 : call CARMAELEMENT_Get(carma, ielem, rc, kappa=kappaelem)
3575 :
3576 >97844*10^9 : kappa(:ncol, :) = kappa(:ncol, :) + mmr(:ncol, :) * kappaelem
3577 >97950*10^9 : totmmr(:ncol, :) = totmmr(:ncol, :) + mmr(:ncol, :)
3578 : end do
3579 :
3580 : ! Figure out the average kappa.q
3581 >16307*10^9 : where (totmmr(:ncol,:) .gt. 0._r8)
3582 15097198999 : kappa(:ncol,:) = kappa(:ncol,:) / totmmr(:ncol,:)
3583 : end where
3584 :
3585 : return
3586 : end subroutine
3587 :
3588 :
3589 : !! Get the number mixing ratio for the group. This is the number of particles per
3590 : !! density of air.
3591 : !!
3592 : !! @author Chuck Bardeen
3593 : !! @version Aug 2023
3594 1103948400 : subroutine carma_get_number(state, igroup, ibin, nmr, rc)
3595 :
3596 : implicit none
3597 :
3598 : type(physics_state), intent(in) :: state !! physics state variables
3599 : integer, intent(in) :: igroup !! group index
3600 : integer, intent(in) :: ibin !! bin index
3601 : real(r8), intent(out) :: nmr(pcols,pver) !! number mixing ratio (#/kg)
3602 : integer, intent(out) :: rc !! return code
3603 :
3604 2207896800 : real(r8) :: rmass(carma%f_NBIN) ! the bin mass (g)
3605 : real(r8) :: totmmr(pcols,pver) ! total mmr (kg/kg)
3606 : integer :: ncol
3607 :
3608 : ! default return code
3609 1103948400 : rc = RC_OK
3610 :
3611 1103948400 : ncol = state%ncol
3612 :
3613 : ! Check the group and bin ranges
3614 1103948400 : if ((igroup < 1) .or. (igroup .gt. NGROUP)) then
3615 0 : write(LUNOPRT, *) 'carma_get_number:: ERROR - Invalid group id, ', igroup
3616 0 : rc = RC_ERROR
3617 0 : return
3618 : end if
3619 :
3620 1103948400 : if ((ibin < 1) .or. (igroup .gt. NBIN)) then
3621 0 : write(LUNOPRT, *) 'carma_get_number:: ERROR - Invalid bin id, ', ibin
3622 0 : rc = RC_ERROR
3623 0 : return
3624 : end if
3625 :
3626 : ! Get the mass in each bin
3627 1103948400 : call CARMAGROUP_Get(carma, igroup, rc, rmass=rmass)
3628 1103948400 : if (rc < 0) return
3629 :
3630 : ! Get the total mmr in the bin
3631 1103948400 : call carma_get_total_mmr(state, igroup, ibin, totmmr, rc)
3632 1103948400 : if (rc < 0) return
3633 :
3634 : ! Get the mmr is the total mass divided by rmass, but need to convert rmass
3635 : ! to kg.
3636 >11919*10^8 : nmr(:ncol, :) = totmmr(:ncol, :) / (rmass(ibin) / 1.e3_r8)
3637 :
3638 : return
3639 : end subroutine carma_get_number
3640 :
3641 877618800 : subroutine carma_get_number_cld(pbuf, igroup, ibin, ncol, nlev, nmr, rc)
3642 :
3643 : implicit none
3644 :
3645 : type(physics_buffer_desc),pointer :: pbuf(:) !! physics buffer
3646 : integer, intent(in) :: igroup !! group index
3647 : integer, intent(in) :: ibin !! bin index
3648 : integer, intent(in) :: ncol,nlev !! dimensions
3649 : real(r8), intent(out) :: nmr(pcols,pver) !! number mixing ratio (#/kg)
3650 : integer, intent(out) :: rc !! return code
3651 :
3652 1755237600 : real(r8) :: rmass(carma%f_NBIN) ! the bin mass (g)
3653 : real(r8) :: totmmr(pcols,pver) ! total mmr (kg/kg)
3654 :
3655 : ! default return code
3656 877618800 : rc = RC_OK
3657 :
3658 : ! Check the group and bin ranges
3659 877618800 : if ((igroup < 1) .or. (igroup .gt. NGROUP)) then
3660 0 : write(LUNOPRT, *) 'carma_get_number:: ERROR - Invalid group id, ', igroup
3661 0 : rc = RC_ERROR
3662 0 : return
3663 : end if
3664 :
3665 877618800 : if ((ibin < 1) .or. (igroup .gt. NBIN)) then
3666 0 : write(LUNOPRT, *) 'carma_get_number:: ERROR - Invalid bin id, ', ibin
3667 0 : rc = RC_ERROR
3668 0 : return
3669 : end if
3670 :
3671 : ! Get the mass in each bin
3672 877618800 : call CARMAGROUP_Get(carma, igroup, rc, rmass=rmass)
3673 877618800 : if (rc < 0) return
3674 :
3675 : ! Get the total mmr in the bin
3676 877618800 : call carma_get_total_mmr_cld(pbuf, igroup, ibin, ncol, nlev, totmmr, rc)
3677 877618800 : if (rc < 0) return
3678 :
3679 : ! Get the mmr is the total mass divided by rmass, but need to convert rmass
3680 : ! to kg.
3681 >94771*10^7 : nmr(:ncol, :) = totmmr(:ncol, :) / (rmass(ibin) / 1.e3_r8)
3682 :
3683 : return
3684 : end subroutine carma_get_number_cld
3685 :
3686 :
3687 : !! Get the mixing ratio for the group. This is the total of all the elements that
3688 : !! make up the group.
3689 : !!
3690 : !! @author Chuck Bardeen
3691 : !! @version Aug 2023
3692 2207896800 : subroutine carma_get_total_mmr(state, igroup, ibin, totmmr, rc)
3693 :
3694 : implicit none
3695 :
3696 : type(physics_state), intent(in) :: state !! physics state variables
3697 : integer, intent(in) :: igroup !! group index
3698 : integer, intent(in) :: ibin !! bin index
3699 : real(r8), intent(out) :: totmmr(pcols,pver) !! total mmr (kg/kg)
3700 : integer, intent(out) :: rc !! return code
3701 :
3702 : real(r8) :: mmr(pcols, pver) ! mmr (kg/kg)
3703 : integer :: i
3704 : integer :: nelems
3705 : integer :: ielems(NELEM)
3706 : integer :: ielem
3707 : integer :: ncol
3708 :
3709 : ! default return code
3710 1103948400 : rc = RC_OK
3711 :
3712 1103948400 : ncol = state%ncol
3713 :
3714 : ! Check the group and bin ranges
3715 1103948400 : if ((igroup < 1) .or. (igroup .gt. NGROUP)) then
3716 0 : write(LUNOPRT, *) 'carma_get_total_mmr:: ERROR - Invalid group id, ', igroup
3717 0 : rc = RC_ERROR
3718 0 : return
3719 : end if
3720 :
3721 1103948400 : if ((ibin < 1) .or. (ibin .gt. NBIN)) then
3722 0 : write(LUNOPRT, *) 'carma_get_total_mmr:: ERROR - Invalid bin id, ', ibin
3723 0 : rc = RC_ERROR
3724 0 : return
3725 : end if
3726 :
3727 : ! Iterate over all of the composition and determine the total mass.
3728 1103948400 : call carma_get_elem_for_group(igroup, nelems, ielems, rc)
3729 1103948400 : if (rc < 0) return
3730 :
3731 >11919*10^8 : totmmr(:ncol, :) = 0._r8
3732 :
3733 5310103800 : do i = 1, nelems
3734 4206155400 : ielem = ielems(i)
3735 :
3736 4206155400 : call carma_get_bin(state, ielem, ibin, mmr, rc)
3737 4206155400 : if (rc < 0) return
3738 :
3739 >45422*10^8 : totmmr(:ncol, :) = totmmr(:ncol, :) + mmr(:ncol, :)
3740 : end do
3741 :
3742 : return
3743 : end subroutine carma_get_total_mmr
3744 :
3745 1755237600 : subroutine carma_get_total_mmr_cld(pbuf, igroup, ibin, ncol, nlev, totmmr, rc)
3746 :
3747 : type(physics_buffer_desc),pointer :: pbuf(:) !! physics buffer
3748 : integer, intent(in) :: igroup !! group index
3749 : integer, intent(in) :: ibin !! bin index
3750 : integer, intent(in) :: ncol,nlev !! dimensions
3751 : real(r8), intent(out) :: totmmr(pcols,pver) !! total mmr (kg/kg)
3752 : integer, intent(out) :: rc !! return code
3753 :
3754 : real(r8) :: mmr(pcols, pver) ! mmr (kg/kg)
3755 : integer :: i
3756 : integer :: nelems
3757 : integer :: ielems(NELEM)
3758 : integer :: ielem
3759 :
3760 : ! default return code
3761 877618800 : rc = RC_OK
3762 :
3763 :
3764 : ! Check the group and bin ranges
3765 877618800 : if ((igroup < 1) .or. (igroup .gt. NGROUP)) then
3766 0 : write(LUNOPRT, *) 'carma_get_total_mmr:: ERROR - Invalid group id, ', igroup
3767 0 : rc = RC_ERROR
3768 0 : return
3769 : end if
3770 :
3771 877618800 : if ((ibin < 1) .or. (ibin .gt. NBIN)) then
3772 0 : write(LUNOPRT, *) 'carma_get_total_mmr:: ERROR - Invalid bin id, ', ibin
3773 0 : rc = RC_ERROR
3774 0 : return
3775 : end if
3776 :
3777 : ! Iterate over all of the composition and determine the total mass.
3778 877618800 : call carma_get_elem_for_group(igroup, nelems, ielems, rc)
3779 877618800 : if (rc < 0) return
3780 :
3781 >94771*10^7 : totmmr(:ncol, :) = 0._r8
3782 :
3783 3973476600 : do i = 1, nelems
3784 3095857800 : ielem = ielems(i)
3785 :
3786 3095857800 : call carma_get_bin_cld(pbuf, ielem, ibin, ncol, nlev, mmr, rc)
3787 3095857800 : if (rc < 0) return
3788 :
3789 >33439*10^8 : totmmr(:ncol, :) = totmmr(:ncol, :) + mmr(:ncol, :)
3790 : end do
3791 :
3792 : end subroutine carma_get_total_mmr_cld
3793 :
3794 5836800 : subroutine carma_get_sad(state, igroup, ibin, sad, rc)
3795 : type(physics_state), intent(in) :: state !! physics state variables
3796 : integer, intent(in) :: igroup !! group index
3797 : integer, intent(in) :: ibin !! bin index
3798 : real(r8), intent(out) :: sad(pcols,pver) !! surface area dens (cm2/cm3)
3799 : integer, intent(out) :: rc !! return code
3800 :
3801 : real(r8) :: nmr(pcols,pver) !! number mixing ratio (#/kg)
3802 : real(r8) :: rwet(pcols,pver) !! wet radius (m)
3803 : real(r8) :: rhopwet(pcols,pver) !! wet density (kg/m3)
3804 : real(r8) :: rhoa(pcols,pver) !! air density (kg/m3)
3805 : real(r8) :: ndens(pcols,pver) !! number density (#/m3)
3806 :
3807 : integer :: ncol
3808 :
3809 5836800 : rc = RC_OK
3810 :
3811 5836800 : call carma_get_wet_radius(state, igroup, ibin, rwet, rhopwet, rc)
3812 5836800 : call carma_get_number(state, igroup, ibin, nmr, rc)
3813 :
3814 5836800 : ncol = state%ncol
3815 :
3816 6297907200 : rhoa(:ncol,:) = (state%pmid(:ncol,:) * 10._r8) / (R_AIR * state%t(:ncol,:)) / 1.e3_r8 * 1.e6_r8 ! air density (kg/m3)
3817 :
3818 6297907200 : ndens(:ncol,:) = nmr(:ncol,:) * rhoa(:ncol,:) ! #/m3
3819 :
3820 6297907200 : sad(:ncol,:) = 4.0_r8 * PI * ndens(:ncol,:) * (rwet(:ncol,:)**2) * 1.e-2_r8 ! cm2/cm3
3821 :
3822 5836800 : end subroutine carma_get_sad
3823 :
3824 :
3825 : !! Find the wet radius and wet density for the group and bin specified.
3826 : !!
3827 : !! NOTE: Groups can be configured with different methods to determine the wet
3828 : !! radius, so multiple methods need to be supported and code from rhopart and
3829 : !! wetr need to be included in this routine.
3830 : !!
3831 : !! @author Chuck Bardeen
3832 : !! @version Aug 2023
3833 89856000 : subroutine carma_get_wet_radius(state, igroup, ibin, rwet, rhopwet, rc)
3834 : use wetr, only: getwetr
3835 :
3836 : type(physics_state), intent(in) :: state !! physics state variables
3837 : integer, intent(in) :: igroup !! group index
3838 : integer, intent(in) :: ibin !! bin index
3839 : real(r8), intent(out) :: rwet(pcols,pver) !! wet radius (m)
3840 : real(r8), intent(out) :: rhopwet(pcols,pver) !! wet density (kg/m3)
3841 : integer, intent(inout) :: rc !! return code
3842 :
3843 : real(r8) :: rdry(pcols,pver) !! dry radius (m)
3844 : real(r8) :: rhopdry(pcols,pver) !! dry density (kg/m3)
3845 : real(r8) :: rhoa(pcols,pver) !! air density (kg/m3)
3846 : real(r8) :: kappa(pcols,pver) !! dry radius (m)
3847 : real(r8) :: es !! saturation vapor pressure
3848 : real(r8) :: qs !! saturation specific humidity
3849 : real(r8) :: relhum !! relative humidity
3850 : real(r8) :: wvpres !! water eq. vaper pressure (dynes/cm2)
3851 : real(r8) :: watcon !! water concentration (g/cm3)
3852 : real(r8) :: dryden !! dry density (g/cm3)
3853 : real(r8) :: dryrad !! dry radius (cm)
3854 : integer :: icol
3855 : integer :: iz
3856 : integer :: ncol
3857 : integer :: iq
3858 : integer :: irhswell
3859 :
3860 : ! default return code
3861 29952000 : rc = RC_OK
3862 :
3863 29952000 : ncol = state%ncol
3864 :
3865 : ! Check the group and bin ranges
3866 29952000 : if ((igroup < 1) .or. (igroup .gt. NGROUP)) then
3867 0 : write(LUNOPRT, *) 'carma_get_total_mmr:: ERROR - Invalid group id, ', igroup
3868 0 : rc = RC_ERROR
3869 : !return
3870 : end if
3871 29952000 : if (rc/=RC_OK) then
3872 0 : call endrun('carma_get_wet_radius ERROR1: rc = ',rc)
3873 : end if
3874 :
3875 29952000 : if ((ibin < 1) .or. (ibin .gt. NBIN)) then
3876 0 : write(LUNOPRT, *) 'carma_get_total_mmr:: ERROR - Invalid bin id, ', ibin
3877 0 : rc = RC_ERROR
3878 : !return
3879 : end if
3880 29952000 : if (rc/=RC_OK) then
3881 0 : call endrun('carma_get_wet_radius ERROR2: rc = ',rc)
3882 : end if
3883 :
3884 : ! Get the constiuent index for water vapor (Q)
3885 29952000 : call cnst_get_ind("Q", iq)
3886 :
3887 : ! The wet radius can be configured differently for each group, so we
3888 : ! need to use getwetr to handle those differences. This requires repeating
3889 : ! some code that is in rhopart to use getwetr properly. There may be a
3890 : ! better way to do this, but for now we will duplicate the code.
3891 29952000 : call carma_get_dry_radius(state, igroup, ibin, rdry, rhopdry, rc)
3892 : !if (rc < 0) return
3893 29952000 : if (rc/=RC_OK) then
3894 0 : call endrun('carma_get_wet_radius ERROR3: rc = ',rc)
3895 : end if
3896 :
3897 : ! Calculate the dry air density at each level, using the ideal gas law.
3898 32318208000 : rhoa(:ncol, :) = (state%pmid(:ncol, :) * 10._r8) / (R_AIR * state%t(:ncol, :)) / 1.e3_r8 * 1.e6_r8
3899 :
3900 29952000 : call CARMAGROUP_Get(carma, igroup, rc, irhswell=irhswell)
3901 : !if (rc < 0) return
3902 29952000 : if (rc/=RC_OK) then
3903 0 : call endrun('carma_get_wet_radius ERROR4: rc = ',rc)
3904 : end if
3905 :
3906 461260800 : do icol = 1, ncol
3907 30652876800 : do iz = 1, pver
3908 30622924800 : if (rdry(icol, iz)>0._r8) then
3909 : ! Get relative humidity and vapor pressure
3910 29310079993 : call wv_sat_qsat_water(state%t(icol,iz), state%pmid(icol,iz), es, qs)
3911 :
3912 : ! NOTE: getwetr is in cgs units, so some conversions are needed from the
3913 : ! mks values
3914 29310079993 : wvpres = es * 10._r8 ! dynes/cm2
3915 29310079993 : relhum = state%q(icol,iz,iq) / qs
3916 29310079993 : watcon = state%q(icol,iz,iq) * rhoa(icol, iz) * 1.e-3_r8 ! g/cm3
3917 29310079993 : dryden = rhopdry(icol,iz) * 1.e-3_r8 ! g/cm3
3918 29310079993 : dryrad = rdry(icol,iz) * 1.e2_r8 ! cm
3919 :
3920 : ! If humidity affects the particle, then determine the equilbirium
3921 : ! radius and density based upon the relative humidity.
3922 : !
3923 29310079993 : if (irhswell == I_WTPCT_H2SO4) then
3924 :
3925 : call getwetr(carma, igroup, relhum, dryrad, rwet(icol, iz), dryden, rhopwet(icol,iz), rc, &
3926 : h2o_mass=watcon, h2o_vp=wvpres, temp=state%t(icol,iz))
3927 14214416994 : if (rc/=RC_OK) then
3928 0 : call endrun('carma_get_wet_radius ERROR5: rc = ',rc) ! <======
3929 : end if
3930 :
3931 30191325998 : else if (irhswell == I_PETTERS) then
3932 :
3933 15095662999 : call carma_get_kappa(state, igroup, ibin, kappa, rc)
3934 15095662999 : if (rc/=RC_OK) then
3935 0 : call endrun('carma_get_wet_radius carma_get_kappa ERROR: rc = ',rc)
3936 : end if
3937 :
3938 : call getwetr(carma, igroup, relhum, dryrad, rwet(icol, iz), dryden, rhopwet(icol,iz), rc, &
3939 : h2o_mass=watcon, h2o_vp=wvpres, temp=state%t(icol,iz), kappa=kappa(icol,iz))
3940 15095662999 : if (rc/=RC_OK) then
3941 0 : call endrun('carma_get_wet_radius ERROR6: rc = ',rc)
3942 : end if
3943 :
3944 : else ! I_GERBER and I_FITZGERALD
3945 :
3946 : call getwetr(carma, igroup, relhum, dryrad, rwet(icol, iz), dryden, rhopwet(icol,iz), rc )
3947 0 : if (rc/=RC_OK) then
3948 0 : call endrun('carma_get_wet_radius ERROR7: rc = ',rc)
3949 : end if
3950 :
3951 : end if
3952 : else
3953 881536007 : rhopwet(icol,iz) = 0._r8
3954 881536007 : rwet(icol, iz) = 0._r8
3955 : end if
3956 : end do
3957 : end do
3958 :
3959 : ! Convert rwet and rhopwet to mks units
3960 32318208000 : rwet(:ncol,:) = rwet(:ncol,:) * 1.e-2_r8 ! cm --> m
3961 32318208000 : rhopwet(:ncol,:) = rhopwet(:ncol,:) * 1.e3_r8 ! g/cm3 --> kg/m3
3962 :
3963 29952000 : if (rc/=RC_OK) then
3964 0 : call endrun('carma_get_wet_radius ERROR8: rc = ',rc)
3965 : end if
3966 :
3967 29952000 : return
3968 : end subroutine
3969 :
3970 :
3971 : !! Provides the tendency (in kg/kg/s) required to change the element and bin from
3972 : !! the current state to the desired mmr.
3973 : !!
3974 : !! NOTE: The caller needs to make sure that the lq flags are set in ptend for the
3975 : !! particular tracer. Perhaps we need a routine that will set lq to true for all
3976 : !! the fields that could be set by CARMA to be used by the caller of this routine.
3977 : !!
3978 : !! @author Chuck Bardeen
3979 : !! @version Aug 2023
3980 0 : subroutine carma_set_bin(state, ielem, ibin, mmr, dt, ptend, rc)
3981 :
3982 : implicit none
3983 :
3984 : type(physics_state), intent(in) :: state !! physics state variables
3985 : integer, intent(in) :: ielem !! element index
3986 : integer, intent(in) :: ibin !! bin index
3987 : real(r8), intent(in) :: mmr(pcols,pver) !! mass mixing ratio (kg/kg)
3988 : integer :: dt !! timestep size (sec)
3989 : type(physics_ptend), intent(inout) :: ptend !! constituent tendencies
3990 : integer, intent(out) :: rc !! return code
3991 :
3992 : integer :: ncol
3993 : integer :: icnst
3994 :
3995 : ! default return code
3996 0 : rc = RC_OK
3997 :
3998 0 : ncol = state%ncol
3999 :
4000 : ! Check the element and bin ranges
4001 0 : if ((ielem < 1) .or. (ielem .gt. NELEM)) then
4002 0 : write(LUNOPRT, *) 'carma_set_bin:: ERROR - Invalid element id, ', ielem
4003 0 : rc = RC_ERROR
4004 0 : return
4005 : end if
4006 :
4007 0 : if ((ibin < 1) .or. (ibin .gt. NBIN)) then
4008 0 : write(LUNOPRT, *) 'carma_set_binr:: ERROR - Invalid bin id, ', ibin
4009 0 : rc = RC_ERROR
4010 0 : return
4011 : end if
4012 :
4013 : ! Determine the tendency needed to make state into mmr for this tracer.
4014 0 : icnst = icnst4elem(ielem, ibin)
4015 0 : ptend%q(:ncol, :, icnst) = (mmr(:ncol, :) - state%q(:ncol, :, icnst)) / dt
4016 :
4017 : return
4018 : end subroutine
4019 :
4020 346037273 : subroutine carma_get_bin_rmass(igroup, ibin, mass, rc)
4021 :
4022 : integer, intent(in) :: igroup !! group index
4023 : integer, intent(in) :: ibin !! bin index
4024 : real(r8),intent(out) :: mass ! grams ???
4025 : integer, intent(out) :: rc !! return code
4026 :
4027 0 : real(r8) :: rmass(carma%f_NBIN) ! the bin mass (g)
4028 :
4029 : ! default return code
4030 346037273 : rc = RC_OK
4031 7266782733 : rmass = rmass
4032 :
4033 346037273 : call CARMAGROUP_Get(carma, igroup, rc, rmass=rmass) ! rmass in g
4034 346037273 : if (rc /= RC_OK) return
4035 :
4036 346037273 : mass = rmass(ibin)*1.e-03_r8 ! convert to kg
4037 :
4038 346037273 : end subroutine carma_get_bin_rmass
4039 :
4040 76800 : function carma_get_wght_pct(ncol,nlev,state) result(wtpct)
4041 : use sulfate_utils, only: wtpct_tabaz
4042 :
4043 : integer, intent(in) :: ncol,nlev
4044 : type(physics_state), intent(in) :: state !! Physics state variables - before CARMA
4045 :
4046 : real(r8) :: wtpct(ncol,nlev)
4047 :
4048 : integer :: rc !! return code
4049 : real(r8) :: pvapl, es, qs, gc_cgs, rhoa
4050 : integer :: icol, ilev
4051 :
4052 76800 : rc = RC_OK
4053 :
4054 5452800 : do ilev = 1,nlev
4055 82867200 : do icol = 1,ncol
4056 : ! Get relative humidity and vapor pressure
4057 :
4058 77414400 : call wv_sat_qsat_water(state%t(icol,ilev), state%pmid(icol,ilev), es, qs) ! es = Saturation vapor pressure in Pa
4059 :
4060 77414400 : pvapl = es * 10._r8 ! Pa -> dynes/cm2
4061 :
4062 77414400 : rhoa = (state%pmid(icol,ilev) * 10._r8) / (R_AIR * state%t(icol,ilev)) ! grams/cm3
4063 :
4064 77414400 : gc_cgs = state%q(icol,ilev,icnst4gas(carma%f_igash2o)) * rhoa ! h2o grams/cm3
4065 :
4066 77414400 : wtpct(icol,ilev) = wtpct_tabaz(carma, state%t(icol,ilev), gc_cgs, pvapl, rc)
4067 :
4068 160204800 : if (rc/=RC_OK) then
4069 0 : call endrun('carma_get_wght_pct: rc = ',rc)
4070 : end if
4071 : end do
4072 : end do
4073 :
4074 153600 : end function carma_get_wght_pct
4075 :
4076 145920 : function carma_effecitive_radius(state) result(rad)
4077 :
4078 : type(physics_state), intent(in) :: state !! physics state variables
4079 : real(r8) :: rad(pcols,pver) ! effective radius (cm)
4080 :
4081 : integer :: igroup, ibin, rc, ncol
4082 : real(r8) :: rwet(pcols,pver) !! wet radius (m)
4083 : real(r8) :: rho(pcols,pver) !! density (kg/m3)
4084 : real(r8) :: nmr(pcols,pver) !! num/kg
4085 : real(r8) :: rtmp3(pcols,pver)
4086 : real(r8) :: rtmp2(pcols,pver)
4087 :
4088 145920 : rc = RC_OK
4089 :
4090 145920 : rtmp2(:,:) = 0.0_r8
4091 145920 : rtmp3(:,:) = 0.0_r8
4092 :
4093 145920 : ncol = state%ncol
4094 :
4095 437760 : do igroup = 1, NGROUP
4096 6274560 : do ibin = 1, NBIN
4097 :
4098 5836800 : call carma_get_number(state, igroup, ibin, nmr, rc)
4099 5836800 : call carma_get_wet_radius(state, igroup, ibin, rwet, rho, rc)
4100 5836800 : if (rc/=RC_OK) then
4101 0 : call endrun('carma_effecitive_radius -- carma_get_wet_radius ERROR: rc = ',rc)
4102 : end if
4103 :
4104 6297907200 : rtmp3(:ncol,:) = rtmp3(:ncol,:) + nmr(:ncol,:)*(rwet(:ncol,:)**3)
4105 6304035840 : rtmp2(:ncol,:) = rtmp2(:ncol,:) + nmr(:ncol,:)*(rwet(:ncol,:)**2)
4106 :
4107 : end do
4108 : end do
4109 :
4110 157447680 : rad(:ncol,:) = (rtmp3(:ncol,:)/rtmp2(:ncol,:))*100._r8 ! cm
4111 :
4112 222720 : end function carma_effecitive_radius
4113 :
4114 : end module carma_intr
|
