Line data Source code
1 : !! The CARMA state module contains the atmospheric data for use with the CARMA
2 : !! module. This implementation has been customized to work within other model
3 : !! frameworks. CARMA adds a lot of extra state information (atmospheric
4 : !! properties, fall velocities, coagulation kernels, growth kernels, ...) and
5 : !! thus has a large memory footprint. Because only one column will be operated
6 : !! upon at a time per thread, only one cstate object needs to be instantiated
7 : !! at a time and each cstate object only represents one column. This keeps
8 : !! the memory requirements of CARMA to a minimum.
9 : !!
10 : !! @version Feb-2009
11 : !! @author Chuck Bardeen, Pete Colarco, Jamie Smith
12 : !
13 : ! NOTE: Documentation for this code can be generated automatically using f90doc,
14 : ! which is freely available from:
15 : ! http://erikdemaine.org/software/f90doc/
16 : ! Comment lines with double comment characters are processed by f90doc, and there are
17 : ! some special characters added to the comments to control the documentation process.
18 : ! In addition to the special characters mentioned in the f990doc documentation, html
19 : ! formatting tags (e.g. <i></i>, <sup></sup>, ...) can also be added to the f90doc
20 : ! comments.
21 : module carmastate_mod
22 :
23 : ! This module maps the parents models constants into the constants need by CARMA.
24 : ! NOTE: CARMA constants are in CGS units, while the parent models are typically in
25 : ! MKS units.
26 : use carma_precision_mod
27 : use carma_enums_mod
28 : use carma_constants_mod
29 : use carma_types_mod
30 :
31 : ! cstate explicitly declares all variables.
32 : implicit none
33 :
34 : ! All cstate variables and procedures are private except those explicitly
35 : ! declared to be public.
36 : private
37 :
38 : ! Declare the public methods.
39 : public CARMASTATE_Create
40 : public CARMASTATE_CreateFromReference
41 : public CARMASTATE_Destroy
42 : public CARMASTATE_Get
43 : public CARMASTATE_GetBin
44 : public CARMASTATE_GetDetrain
45 : public CARMASTATE_GetGas
46 : public CARMASTATE_GetState
47 : public CARMASTATE_SetBin
48 : public CARMASTATE_SetDetrain
49 : public CARMASTATE_SetGas
50 : public CARMASTATE_SetState
51 : public CARMASTATE_Step
52 :
53 : contains
54 :
55 : ! These are the methods that provide the interface between the parent model and
56 : ! the atmospheric state data of the CARMA microphysical model. There are many other
57 : ! methods that are not in this file that are used to implement the microphysical
58 : ! calculations needed by the CARMA model. These other methods are in effect private
59 : ! methods of the CARMA module, but are in individual files since that is the way that
60 : ! CARMA has traditionally been structured and where users may want to extend or
61 : ! replace code to affect the microphysics.
62 :
63 : !! Create the CARMASTATE object, which contains information about the
64 : !! atmospheric state. Internally, CARMA uses CGS units, but this interface uses
65 : !! MKS units which are more commonly used in parent models. The units and grid
66 : !! orientation depend on the grid type:
67 : !!
68 : !! - igridv
69 : !! - I_CART : Cartesian coordinates, units in [m], bottom at NZ=1
70 : !! - I_SIG : Sigma coordinates, unitless [P/P0], top at NZ=1
71 : !! - I_HYBRID : Hybrid coordinates, unitless [~P/P0], top at NZ=1
72 : !!
73 : !! NOTE: The supplied CARMA object should already have been created, configured,
74 : !! and initialized.
75 : !!
76 : !! NOTE: The relative humidity is optional, but needs to be supplied if particles
77 : !! are subject to swelling based upon relative humidity. The specific humdity can
78 : !! can be specified instead. If both are specified, then the realtive humidity is
79 : !! used.
80 : !!
81 : !! @author Chuck Bardeen
82 : !! @version Feb-2009
83 : !! @see CARMA_Create
84 : !! @see CARMA_Initialize
85 : !! @see CARMASTATE_Destroy
86 1050624 : subroutine CARMASTATE_Create(cstate, carma_ptr, time, dtime, NZ, igridv, &
87 3151872 : xc, yc, zc, zl, p, pl, t, rc, qh2o, relhum, told, radint)
88 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
89 : type(carma_type), pointer, intent(in) :: carma_ptr !! (in) the carma object
90 : real(kind=f), intent(in) :: time !! the model time [s]
91 : real(kind=f), intent(in) :: dtime !! the timestep size [s]
92 : integer, intent(in) :: NZ !! the number of vertical grid points
93 : integer, intent(in) :: igridv !! vertical grid type
94 : real(kind=f), intent(in) :: xc !! x at center
95 : real(kind=f), intent(in) :: yc !! y at center
96 : real(kind=f), intent(in) :: zc(NZ) !! z at center
97 : real(kind=f), intent(in) :: zl(NZ+1) !! z at edge
98 : real(kind=f), intent(in) :: p(NZ) !! pressure at center [Pa]
99 : real(kind=f), intent(in) :: pl(NZ+1) !! presssure at edge [Pa]
100 : real(kind=f), intent(in) :: t(NZ) !! temperature at center [K]
101 : integer, intent(out) :: rc !! return code, negative indicates failure
102 : real(kind=f), intent(in) , optional :: qh2o(NZ) !! specific humidity at center [mmr]
103 : real(kind=f), intent(in) , optional :: relhum(NZ) !! relative humidity at center [fraction]
104 : real(kind=f), intent(in) , optional :: told(NZ) !! previous temperature at center [K]
105 : real(kind=f), intent(in) , optional :: radint(NZ,carma_ptr%f_NWAVE) !! radiative intensity [W/m2/sr/cm]
106 :
107 : integer :: iz
108 : real(kind=f) :: rvap
109 : real(kind=f) :: pvap_liq
110 : real(kind=f) :: pvap_ice
111 : real(kind=f) :: gc_cgs
112 :
113 : ! Assume success.
114 1050624 : rc = RC_OK
115 :
116 : ! Save the defintion of the number of comonents involved in the microphysics.
117 1050624 : cstate%f_carma => carma_ptr
118 :
119 : ! Save the model timing.
120 1050624 : cstate%f_time = time
121 1050624 : cstate%f_dtime_orig = dtime
122 1050624 : cstate%f_dtime = dtime
123 1050624 : cstate%f_nretries = 0
124 :
125 : ! Save the grid dimensions.
126 1050624 : cstate%f_NZ = NZ
127 1050624 : cstate%f_NZP1 = NZ+1
128 :
129 : ! Save the grid definition.
130 1050624 : cstate%f_igridv = igridv
131 :
132 : ! Store away the grid location information.
133 1050624 : cstate%f_yc = yc
134 1050624 : cstate%f_xc = yc
135 :
136 : ! Allocate all the dynamic variables related to state.
137 1050624 : call CARMASTATE_Allocate(cstate, rc)
138 1050624 : if (rc < 0) return
139 :
140 34670592 : cstate%f_zc(:) = zc(:)
141 35721216 : cstate%f_zl(:) = zl(:)
142 :
143 : ! Store away the grid state, doing any necessary unit conversions from MKS to CGS.
144 34670592 : cstate%f_p(:) = p(:) * RPA2CGS
145 35721216 : cstate%f_pl(:) = pl(:) * RPA2CGS
146 34670592 : cstate%f_t(:) = t(:)
147 :
148 7640137728 : cstate%f_pcd(:,:,:) = 0._f
149 :
150 1050624 : if (carma_ptr%f_do_substep) then
151 1050624 : if (present(told)) then
152 34670592 : cstate%f_told(:) = told
153 : else
154 0 : if (carma_ptr%f_do_print) write(carma_ptr%f_LUNOPRT,*) "CARMASTATE_Create: Error - Need to specify told when substepping."
155 0 : rc = RC_ERROR
156 :
157 0 : return
158 : end if
159 : end if
160 :
161 : ! Calculate the metrics, ...
162 : ! if Cartesian coordinates were specifed, then the units need to be converted
163 : ! from MKS to CGS.
164 1050624 : if (cstate%f_igridv == I_CART) then
165 0 : cstate%f_zc = cstate%f_zc * RM2CGS
166 0 : cstate%f_zl = cstate%f_zl * RM2CGS
167 : end if
168 :
169 : ! Initialize the state of the atmosphere.
170 1050624 : call setupatm(carma_ptr, cstate, carma_ptr%f_do_fixedinit, rc)
171 1050624 : if (rc < 0) return
172 :
173 : ! Set the realtive humidity. If necessary, it will be calculated from
174 : ! the specific humidity.
175 1050624 : if (present(relhum)) then
176 0 : cstate%f_relhum(:) = relhum(:)
177 1050624 : else if (present(qh2o)) then
178 :
179 : ! Define gas constant for this gas
180 1050624 : rvap = RGAS/WTMOL_H2O
181 :
182 : ! Calculate relative humidity
183 34670592 : do iz = 1, NZ
184 33619968 : call vaporp_h2o_murphy2005(carma_ptr, cstate, iz, rc, pvap_liq, pvap_ice)
185 33619968 : if (rc < 0) return
186 :
187 33619968 : gc_cgs = qh2o(iz)*cstate%f_rhoa_wet(iz) / cstate%f_zmet(iz)
188 34670592 : cstate%f_relhum(iz) = ( gc_cgs * rvap * t(iz)) / pvap_liq
189 : enddo
190 : end if
191 :
192 : ! Need for vertical transport.
193 : !
194 : ! NOTE: How should these be set? Optional parameters?
195 1050624 : if (carma_ptr%f_do_vtran) then
196 243744768 : cstate%f_ftoppart(:,:) = 0._f
197 243744768 : cstate%f_fbotpart(:,:) = 0._f
198 243744768 : cstate%f_pc_topbnd(:,:) = 0._f
199 243744768 : cstate%f_pc_botbnd(:,:) = 0._f
200 : end if
201 :
202 : ! Radiative intensity for particle heating.
203 : !
204 : ! W/m2/sr/cm -> erg/s/cm2/sr/cm
205 1050624 : if (carma_ptr%f_do_grow) then
206 1041168384 : if (present(radint)) cstate%f_radint(:,:) = radint(:,:) * 1e7_f / 1e4_f
207 : end if
208 :
209 : return
210 3151872 : end subroutine CARMASTATE_Create
211 :
212 :
213 : !! Create the CARMASTATE object, which contains information about the
214 : !! atmospheric state.
215 : !!
216 : !! This call is similar to CARMASTATE_Create, but differs in that all the
217 : !! initialization happens here based on the the fixed state information provided rather
218 : !! than occurring in CARMASTATE_Step.
219 : !!
220 : !! This call should be done before CARMASTATE_Create when do_fixedinit has been
221 : !! specified. The temperatures and pressures specified here should be the reference
222 : !! state used for all columns, not an actual column from the model.
223 : !!
224 : !! A water vapor profile is optional, but is used whenever either qh2o (preferred)
225 : !! or relhum have been provided. If this is not provided, then initialization will
226 : !! be done on a dry profile. If particle swelling occurs, initialization will be
227 : !! done on the wet radius; however, most of the initialized values will not get
228 : !! recalculated as the wet radius changes.
229 : !!
230 : !! CARMASTATE_Create should still be called again after this call with the actual
231 : !! column of state information from the model. The initialization will be done once
232 : !! from the reference state, but the microphysical calculations will be done on the
233 : !! model state. Multiple CARMASTATE_Create ... CARMASTATE_Step calls can be done
234 : !! before a CARMASTATE_Destroy. This reduces the amount of memory allocations and
235 : !! when used with do_fixedinit, reduces the amount of time spent initializing.
236 : !!
237 : !! @author Chuck Bardeen
238 : !! @version June-2010
239 : !! @see CARMA_Create
240 : !! @see CARMA_Initialize
241 : !! @see CARMASTATE_Destroy
242 0 : subroutine CARMASTATE_CreateFromReference(cstate, carma_ptr, time, dtime, NZ, igridv, &
243 0 : xc, yc, zc, zl, p, pl, t, rc, qh2o, relhum, qh2so4)
244 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
245 : type(carma_type), pointer, intent(in) :: carma_ptr !! (in) the carma object
246 : real(kind=f), intent(in) :: time !! the model time [s]
247 : real(kind=f), intent(in) :: dtime !! the timestep size [s]
248 : integer, intent(in) :: NZ !! the number of vertical grid points
249 : integer, intent(in) :: igridv !! vertical grid type
250 : real(kind=f), intent(in) :: xc !! x at center
251 : real(kind=f), intent(in) :: yc !! y at center
252 : real(kind=f), intent(in) :: zc(NZ) !! z at center
253 : real(kind=f), intent(in) :: zl(NZ+1) !! z at edge
254 : real(kind=f), intent(in) :: p(NZ) !! pressure at center [Pa]
255 : real(kind=f), intent(in) :: pl(NZ+1) !! presssure at edge [Pa]
256 : real(kind=f), intent(in) :: t(NZ) !! temperature at center [K]
257 : integer, intent(out) :: rc !! return code, negative indicates failure
258 : real(kind=f), intent(in) , optional :: qh2o(NZ) !! specific humidity at center [mmr]
259 : real(kind=f), intent(in) , optional :: relhum(NZ) !! relative humidity at center [fraction]
260 : real(kind=f), intent(in) , optional :: qh2so4(NZ) !! H2SO4 mass mixing ratio at center [mmr]
261 :
262 : integer :: iz
263 : integer :: igas
264 : real(kind=f) :: rvap
265 : real(kind=f) :: pvap_liq
266 : real(kind=f) :: pvap_ice
267 : real(kind=f) :: gc_cgs
268 :
269 : ! Assume success.
270 0 : rc = RC_OK
271 :
272 : ! Save the defintion of the number of comonents involved in the microphysics.
273 0 : cstate%f_carma => carma_ptr
274 :
275 : ! Save the model timing.
276 0 : cstate%f_time = time
277 0 : cstate%f_dtime_orig = dtime
278 0 : cstate%f_dtime = dtime
279 0 : cstate%f_nretries = 0
280 :
281 : ! Save the grid dimensions.
282 0 : cstate%f_NZ = NZ
283 0 : cstate%f_NZP1 = NZ+1
284 :
285 : ! Save the grid definition.
286 0 : cstate%f_igridv = igridv
287 :
288 : ! Store away the grid location information.
289 0 : cstate%f_xc = xc
290 0 : cstate%f_yc = yc
291 :
292 : ! Allocate all the dynamic variables related to state.
293 0 : call CARMASTATE_Allocate(cstate, rc)
294 0 : if (rc < 0) return
295 :
296 0 : cstate%f_zc(:) = zc(:)
297 0 : cstate%f_zl(:) = zl(:)
298 :
299 : ! Store away the grid state, doing any necessary unit conversions from MKS to CGS.
300 0 : cstate%f_p(:) = p(:) * RPA2CGS
301 0 : cstate%f_pl(:) = pl(:) * RPA2CGS
302 0 : cstate%f_t(:) = t(:)
303 :
304 0 : cstate%f_pcd(:,:,:) = 0._f
305 :
306 : ! Calculate the metrics, ...
307 : ! if Cartesian coordinates were specifed, then the units need to be converted
308 : ! from MKS to CGS.
309 0 : if (cstate%f_igridv == I_CART) then
310 0 : cstate%f_zc = cstate%f_zc * RM2CGS
311 0 : cstate%f_zl = cstate%f_zl * RM2CGS
312 : end if
313 :
314 : ! Initialize the state of the atmosphere.
315 0 : call setupatm(carma_ptr, cstate, .false., rc)
316 0 : if (rc < 0) return
317 :
318 : ! If the model uses a gas, then set the relative and
319 : ! specific humidities.
320 0 : if (carma_ptr%f_igash2o /= 0) then
321 :
322 0 : if (present(qh2o)) then
323 0 : cstate%f_gc(:, carma_ptr%f_igash2o) = qh2o(:) * cstate%f_rhoa_wet(:)
324 :
325 : ! Define gas constant for this gas
326 0 : rvap = RGAS/WTMOL_H2O
327 :
328 : ! Calculate relative humidity
329 0 : do iz = 1, NZ
330 0 : call vaporp_h2o_murphy2005(carma_ptr, cstate, iz, rc, pvap_liq, pvap_ice)
331 0 : if (rc < 0) return
332 :
333 0 : gc_cgs = qh2o(iz) * cstate%f_rhoa_wet(iz) / cstate%f_zmet(iz)
334 0 : cstate%f_relhum(iz) = (gc_cgs * rvap * t(iz)) / pvap_liq
335 : enddo
336 :
337 0 : else if (present(relhum)) then
338 0 : cstate%f_relhum(:) = relhum
339 :
340 : ! Define gas constant for this gas
341 0 : rvap = RGAS/WTMOL_H2O
342 :
343 : ! Calculate specific humidity
344 0 : do iz = 1, NZ
345 0 : call vaporp_h2o_murphy2005(carma_ptr, cstate, iz, rc, pvap_liq, pvap_ice)
346 0 : if (rc < 0) return
347 :
348 0 : gc_cgs = (rvap * t(iz)) / (pvap_liq * relhum(iz))
349 0 : cstate%f_gc(iz, carma_ptr%f_igash2o) = gc_cgs * &
350 0 : cstate%f_zmet(iz) / &
351 0 : cstate%f_rhoa_wet(iz)
352 : enddo
353 : end if
354 : end if
355 :
356 : ! If the model uses sulfuric acid, then set that gas concentration.
357 0 : if (carma_ptr%f_igash2so4 /= 0) then
358 0 : if (present(qh2so4)) then
359 0 : cstate%f_gc(:, carma_ptr%f_igash2so4) = qh2so4(:) * cstate%f_rhoa_wet(:)
360 : end if
361 : end if
362 :
363 : ! Determine the gas supersaturations.
364 0 : do iz = 1, cstate%f_NZ
365 0 : do igas = 1, cstate%f_carma%f_NGAS
366 0 : call supersat(cstate%f_carma, cstate, iz, igas, rc)
367 0 : if (rc < 0) return
368 : end do
369 : end do
370 :
371 : ! Need for vertical transport.
372 : !
373 : ! NOTE: How should these be set? Optional parameters?
374 0 : if (carma_ptr%f_do_vtran) then
375 0 : cstate%f_ftoppart(:,:) = 0._f
376 0 : cstate%f_fbotpart(:,:) = 0._f
377 0 : cstate%f_pc_topbnd(:,:) = 0._f
378 0 : cstate%f_pc_botbnd(:,:) = 0._f
379 : end if
380 :
381 :
382 : ! Now do the initialization that is normally done in CARMASTATE_Step. However
383 : ! here it is done using the reference atmosphere.
384 :
385 : ! Determine the particle densities.
386 0 : call rhopart(cstate%f_carma, cstate, rc)
387 0 : if (rc < 0) return
388 :
389 : ! Save off the wet radius and wet density as reference values to be used
390 : ! later to scale process rates based upon changes to the wet radius and
391 : ! wet density when particle swelling is used.
392 0 : cstate%f_r_ref(:,:,:) = cstate%f_r_wet(:,:,:)
393 0 : cstate%f_rhop_ref(:,:,:) = cstate%f_rhop_wet(:,:,:)
394 :
395 : ! If configured for fixed initialization, then we will lose some accuracy
396 : ! in the calculation of the fall velocities, growth kernels, ... and in return
397 : ! will gain a significant performance by not having to initialize as often.
398 :
399 : ! Initialize the vertical transport.
400 0 : if (cstate%f_carma%f_do_vtran .or. cstate%f_carma%f_do_coag .or. cstate%f_carma%f_do_grow) then
401 0 : call setupvf(cstate%f_carma, cstate, rc)
402 :
403 0 : if (cstate%f_carma%f_do_vdiff) then
404 0 : call setupbdif(cstate%f_carma, cstate, rc)
405 : end if
406 : end if
407 :
408 : ! Intialize the nucleation, growth and evaporation.
409 0 : if (cstate%f_carma%f_do_grow) then
410 0 : call setupgrow(cstate%f_carma, cstate, rc)
411 0 : if (rc < 0) return
412 :
413 0 : call setupgkern(cstate%f_carma, cstate, rc)
414 0 : if (rc < 0) return
415 :
416 0 : call setupnuc(cstate%f_carma, cstate, rc)
417 0 : if (rc < 0) return
418 : end if
419 :
420 : ! Initialize the coagulation.
421 0 : if (cstate%f_carma%f_do_coag) then
422 0 : call setupckern(cstate%f_carma, cstate, rc)
423 0 : if (rc < 0) return
424 : end if
425 :
426 : return
427 0 : end subroutine CARMASTATE_CreateFromReference
428 :
429 :
430 1050624 : subroutine CARMASTATE_Allocate(cstate, rc)
431 : type(carmastate_type), intent(inout) :: cstate
432 : integer, intent(out) :: rc
433 :
434 : ! Local Variables
435 : integer :: ier
436 : integer :: NZ
437 : integer :: NZP1
438 : integer :: NGROUP
439 : integer :: NELEM
440 : integer :: NBIN
441 : integer :: NGAS
442 : integer :: NWAVE
443 :
444 : ! Assume success.
445 1050624 : rc = RC_OK
446 :
447 : ! Check to see if the arrays are already allocated. If so, just reuse the
448 : ! existing allocations.
449 :
450 : ! Allocate the variables needed for setupatm.
451 1050624 : if (.not. (allocated(cstate%f_zmet))) then
452 :
453 72960 : NZ = cstate%f_NZ
454 72960 : NZP1 = cstate%f_NZP1
455 72960 : NGROUP = cstate%f_carma%f_NGROUP
456 72960 : NELEM = cstate%f_carma%f_NELEM
457 72960 : NBIN = cstate%f_carma%f_NBIN
458 72960 : NGAS = cstate%f_carma%f_NGAS
459 72960 : NWAVE = cstate%f_carma%f_NWAVE
460 :
461 : allocate( &
462 : cstate%f_zmet(NZ), &
463 : cstate%f_zmetl(NZP1), &
464 : cstate%f_zc(NZ), &
465 : cstate%f_dz(NZ), &
466 : cstate%f_zl(NZP1), &
467 : cstate%f_pc(NZ,NBIN,NELEM), &
468 : cstate%f_pcd(NZ,NBIN,NELEM), &
469 : cstate%f_pc_surf(NBIN,NELEM), &
470 : cstate%f_sedimentationflux(NBIN,NELEM), &
471 : cstate%f_gc(NZ,NGAS), &
472 : cstate%f_cldfrc(NZ), &
473 : cstate%f_rhcrit(NZ), &
474 : cstate%f_rhop(NZ,NBIN,NGROUP), &
475 : cstate%f_r_wet(NZ,NBIN,NGROUP), &
476 : cstate%f_rlow_wet(NZ,NBIN,NGROUP), &
477 : cstate%f_rup_wet(NZ,NBIN,NGROUP), &
478 : cstate%f_rhop_wet(NZ,NBIN,NGROUP), &
479 : cstate%f_r_ref(NZ,NBIN,NGROUP), &
480 : cstate%f_rhop_ref(NZ,NBIN,NGROUP), &
481 : cstate%f_rhoa(NZ), &
482 : cstate%f_rhoa_wet(NZ), &
483 : cstate%f_t(NZ), &
484 : cstate%f_p(NZ), &
485 : cstate%f_pl(NZP1), &
486 : cstate%f_relhum(NZ), &
487 : cstate%f_wtpct(NZ), &
488 : cstate%f_rmu(NZ), &
489 : cstate%f_thcond(NZ), &
490 : cstate%f_thcondnc(NZ,NBIN,NGROUP), &
491 : cstate%f_dpc_sed(NBIN,NELEM), &
492 : cstate%f_pconmax(NZ,NGROUP), &
493 : cstate%f_pcl(NZ,NBIN,NELEM), &
494 : cstate%f_kappahygro(NZ,NBIN,NGROUP), &
495 5107200 : stat=ier)
496 72960 : if (ier /= 0) then
497 0 : if (cstate%f_carma%f_do_print) then
498 : write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_Allocate::&
499 0 : &ERROR allocating atmosphere arrays, status=", ier
500 : end if
501 0 : rc = RC_ERROR
502 0 : return
503 : end if
504 :
505 2407680 : cstate%f_relhum(:) = 0._f
506 530565120 : cstate%f_pc(:,:,:) = 0._f
507 530565120 : cstate%f_pcd(:,:,:) = 0._f
508 16926720 : cstate%f_pc_surf(:,:) = 0._f
509 16926720 : cstate%f_sedimentationflux(:,:) = 0._f
510 2407680 : cstate%f_cldfrc(:) = 1._f
511 2407680 : cstate%f_rhcrit(:) = 1._f
512 :
513 : ! Allocate the last fields if they are needed for substepping.
514 72960 : if (cstate%f_carma%f_do_substep) then
515 : allocate( &
516 : cstate%f_gcl(NZ,NGAS), &
517 : cstate%f_d_gc(NZ,NGAS), &
518 : cstate%f_told(NZ), &
519 : cstate%f_d_t(NZ), &
520 : cstate%f_zsubsteps(NZ), &
521 656640 : stat=ier)
522 72960 : if (ier /= 0) then
523 0 : if (cstate%f_carma%f_do_print) then
524 : write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_Allocate::&
525 0 : &ERROR allocating stepping arrays, status=", ier
526 : end if
527 0 : rc = RC_ERROR
528 0 : return
529 : endif
530 :
531 : ! Initialize
532 4888320 : cstate%f_gcl(:,:) = 0._f
533 4888320 : cstate%f_d_gc(:,:) = 0._f
534 2407680 : cstate%f_told(:) = 0._f
535 2407680 : cstate%f_d_t(:) = 0._f
536 2407680 : cstate%f_zsubsteps(:) = 0._f
537 :
538 : ! When substepping is enabled, we want to initialize these statistics once for
539 : ! the life of the object.
540 72960 : cstate%f_max_nsubstep = 0
541 72960 : cstate%f_max_nretry = 0._f
542 72960 : cstate%f_nstep = 0._f
543 72960 : cstate%f_nsubstep = 0
544 72960 : cstate%f_nretry = 0._f
545 : endif
546 :
547 :
548 : ! Allocate the variables needed for setupvf.
549 : !
550 : ! NOTE: Coagulation and dry deposition also need bpm, vf and re.
551 : if (cstate%f_carma%f_do_vtran .or. cstate%f_carma%f_do_coag .or. &
552 72960 : cstate%f_carma%f_do_grow .or. cstate%f_carma%f_do_drydep) then
553 : allocate( &
554 : cstate%f_bpm(NZ,NBIN,NGROUP), &
555 : cstate%f_vf(NZP1,NBIN,NGROUP), &
556 : cstate%f_re(NZ,NBIN,NGROUP), &
557 : cstate%f_dkz(NZP1,NBIN,NGROUP), &
558 : cstate%f_ftoppart(NBIN,NELEM), &
559 : cstate%f_fbotpart(NBIN,NELEM), &
560 : cstate%f_pc_topbnd(NBIN,NELEM), &
561 : cstate%f_pc_botbnd(NBIN,NELEM), &
562 : cstate%f_vd(NBIN, NGROUP), &
563 1969920 : stat=ier)
564 72960 : if (ier /= 0) then
565 0 : if (cstate%f_carma%f_do_print) then
566 : write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_Allocate::&
567 0 : &ERROR allocating vertical transport arrays, status=", ier
568 : end if
569 0 : rc = RC_ERROR
570 0 : return
571 : endif
572 :
573 : ! Initialize
574 96526080 : cstate%f_bpm(:,:,:) = 0._f
575 99444480 : cstate%f_vf(:,:,:) = 0._f
576 96526080 : cstate%f_re(:,:,:) = 0._f
577 99444480 : cstate%f_dkz(:,:,:) = 0._f
578 16926720 : cstate%f_ftoppart(:,:) = 0._f
579 16926720 : cstate%f_fbotpart(:,:) = 0._f
580 16926720 : cstate%f_pc_topbnd(:,:) = 0._f
581 16926720 : cstate%f_pc_botbnd(:,:) = 0._f
582 3137280 : cstate%f_vd(:, :) = 0._f
583 : end if
584 :
585 :
586 :
587 72960 : if (cstate%f_carma%f_NGAS > 0) then
588 : allocate( &
589 : cstate%f_pvapl(NZ,NGAS), &
590 : cstate%f_pvapi(NZ,NGAS), &
591 : cstate%f_supsatl(NZ,NGAS), &
592 : cstate%f_supsati(NZ,NGAS), &
593 : cstate%f_supsatlold(NZ,NGAS), &
594 : cstate%f_supsatiold(NZ,NGAS), &
595 1021440 : stat=ier)
596 72960 : if (ier /= 0) then
597 0 : if (cstate%f_carma%f_do_print) then
598 : write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_Allocate:: "&
599 0 : //"ERROR allocating gas arrays, status=", ier
600 : end if
601 0 : rc = RC_ERROR
602 0 : return
603 : endif
604 : end if
605 :
606 :
607 72960 : if (cstate%f_carma%f_do_grow) then
608 : allocate( &
609 : cstate%f_diffus(NZ,NGAS), &
610 : cstate%f_rlhe(NZ,NGAS), &
611 : cstate%f_rlhm(NZ,NGAS), &
612 : cstate%f_surfctwa(NZ), &
613 : cstate%f_surfctiw(NZ), &
614 : cstate%f_surfctia(NZ), &
615 : cstate%f_akelvin(NZ,NGAS), &
616 : cstate%f_akelvini(NZ,NGAS), &
617 : cstate%f_ft(NZ,NBIN,NGROUP), &
618 : cstate%f_gro(NZ,NBIN,NGROUP), &
619 : cstate%f_gro1(NZ,NBIN,NGROUP), &
620 : cstate%f_gro2(NZ,NGROUP), &
621 : cstate%f_scrit(NZ,NBIN,NGROUP), &
622 : cstate%f_rnuclg(NBIN,NGROUP,NGROUP),&
623 : cstate%f_rhompe(NBIN,NELEM), &
624 : cstate%f_rnucpe(NBIN,NELEM), &
625 : cstate%f_pc_nucl(NZ,NBIN,NELEM), &
626 : cstate%f_growpe(NBIN,NELEM), &
627 : cstate%f_evappe(NBIN,NELEM), &
628 : cstate%f_evcore(NELEM), &
629 : cstate%f_growlg(NBIN,NGROUP), &
630 : cstate%f_evaplg(NBIN,NGROUP), &
631 : cstate%f_gasprod(NGAS), &
632 : cstate%f_rlheat(NZ), &
633 : cstate%f_radint(NZ,NWAVE), &
634 : cstate%f_partheat(NZ), &
635 : cstate%f_dtpart(NZ,NBIN,NGROUP), &
636 : cstate%f_cmf(NBIN,NGROUP), &
637 : cstate%f_totevap(NBIN,NGROUP), &
638 5107200 : stat=ier)
639 72960 : if (ier /= 0) then
640 0 : if (cstate%f_carma%f_do_print) then
641 : write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_Allocate::&
642 0 : &ERROR allocating growth arrays, status=", ier
643 : end if
644 0 : rc = RC_ERROR
645 0 : return
646 : endif
647 :
648 72303360 : cstate%f_radint(:,:) = 0._f
649 : end if
650 :
651 72960 : if (cstate%f_carma%f_do_coag) then
652 : allocate( &
653 : cstate%f_coaglg(NZ,NBIN,NGROUP), &
654 : cstate%f_coagpe(NZ,NBIN,NELEM), &
655 : cstate%f_ckernel(NZ,NBIN,NBIN,NGROUP,NGROUP), &
656 1094400 : stat = ier)
657 72960 : if (ier /= 0) then
658 0 : if (cstate%f_carma%f_do_print) then
659 : write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_Allocate::&
660 0 : &ERROR allocating coag arrays, status=", ier
661 : end if
662 0 : rc = RC_ERROR
663 0 : return
664 : end if
665 :
666 : ! Initialize
667 96526080 : cstate%f_coaglg(:,:,:) = 0._f
668 530565120 : cstate%f_coagpe(:,:,:) = 0._f
669 3858635520 : cstate%f_ckernel(:,:,:,:,:) = 0._f
670 : end if
671 : end if
672 :
673 : return
674 : end subroutine CARMASTATE_Allocate
675 :
676 :
677 : !! The routine should be called when the carma state object is no longer needed.
678 : !! It deallocates any memory allocations made by CARMA during CARMASTATE_Create(),
679 : !! and failure to call this routine could result in memory leaks.
680 : !!
681 : !! @author Chuck Bardeen
682 : !! @version Feb-2009
683 : !! @see CARMASTATE_Create
684 72960 : subroutine CARMASTATE_Destroy(cstate, rc)
685 : type(carmastate_type), intent(inout) :: cstate
686 : integer, intent(out) :: rc
687 :
688 : ! Local variables
689 : integer :: ier
690 :
691 : ! Assume success.
692 72960 : rc = RC_OK
693 :
694 : ! Check to see if the arrays are already allocated. If so, deallocate them.
695 :
696 : ! Allocate the variables needed for setupatm.
697 72960 : if (allocated(cstate%f_zmet)) then
698 :
699 : deallocate( &
700 : cstate%f_zmet, &
701 : cstate%f_zmetl, &
702 : cstate%f_zc, &
703 : cstate%f_dz, &
704 : cstate%f_zl, &
705 : cstate%f_pc, &
706 : cstate%f_pcd, &
707 : cstate%f_pc_surf, &
708 : cstate%f_sedimentationflux, &
709 : cstate%f_gc, &
710 : cstate%f_cldfrc, &
711 : cstate%f_rhcrit, &
712 : cstate%f_rhop, &
713 : cstate%f_r_wet, &
714 : cstate%f_rlow_wet, &
715 : cstate%f_rup_wet, &
716 : cstate%f_rhop_wet, &
717 : cstate%f_r_ref, &
718 : cstate%f_rhop_ref, &
719 : cstate%f_rhoa, &
720 : cstate%f_rhoa_wet, &
721 : cstate%f_t, &
722 : cstate%f_p, &
723 : cstate%f_pl, &
724 : cstate%f_relhum, &
725 : cstate%f_wtpct, &
726 : cstate%f_rmu, &
727 : cstate%f_thcond, &
728 : cstate%f_thcondnc, &
729 : cstate%f_dpc_sed, &
730 : cstate%f_pconmax, &
731 : cstate%f_pcl, &
732 : cstate%f_kappahygro, &
733 72960 : stat=ier)
734 72960 : if (ier /= 0) then
735 0 : if (cstate%f_carma%f_do_print) then
736 : write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_Destroy::&
737 0 : &ERROR deallocating atmosphere arrays, status=", ier
738 : end if
739 0 : rc = RC_ERROR
740 0 : return
741 : end if
742 :
743 : ! Allocate the last fields if they are needed for substepping stepping.
744 72960 : if (allocated(cstate%f_gcl)) then
745 : deallocate( &
746 : cstate%f_gcl, &
747 : cstate%f_d_gc, &
748 : cstate%f_told, &
749 : cstate%f_d_t, &
750 : cstate%f_zsubsteps, &
751 72960 : stat=ier)
752 72960 : if (ier /= 0) then
753 0 : if (cstate%f_carma%f_do_print) then
754 : write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_Destroy::&
755 0 : &ERROR deallocating stepping arrays, status=", ier
756 : end if
757 0 : rc = RC_ERROR
758 0 : return
759 : endif
760 : endif
761 :
762 : ! Allocate the variables needed for setupvf.
763 : !
764 : ! NOTE: Coagulation also needs bpm, vf and re.
765 72960 : if (allocated(cstate%f_bpm)) then
766 : deallocate( &
767 : cstate%f_bpm, &
768 : cstate%f_vf, &
769 : cstate%f_re, &
770 : cstate%f_dkz, &
771 : cstate%f_ftoppart, &
772 : cstate%f_fbotpart, &
773 : cstate%f_pc_topbnd, &
774 : cstate%f_pc_botbnd, &
775 : cstate%f_vd, &
776 72960 : stat=ier)
777 72960 : if (ier /= 0) then
778 0 : if (cstate%f_carma%f_do_print) then
779 : write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_Destroy::&
780 0 : &ERROR deallocating vertical transport arrays, status=", ier
781 : end if
782 0 : rc = RC_ERROR
783 0 : return
784 : endif
785 : end if
786 :
787 72960 : if (allocated(cstate%f_diffus)) then
788 : deallocate( &
789 : cstate%f_diffus, &
790 : cstate%f_rlhe, &
791 : cstate%f_rlhm, &
792 : cstate%f_surfctwa, &
793 : cstate%f_surfctiw, &
794 : cstate%f_surfctia, &
795 : cstate%f_akelvin, &
796 : cstate%f_akelvini, &
797 : cstate%f_ft, &
798 : cstate%f_gro, &
799 : cstate%f_gro1, &
800 : cstate%f_gro2, &
801 : cstate%f_scrit, &
802 : cstate%f_rnuclg,&
803 : cstate%f_rnucpe, &
804 : cstate%f_rhompe, &
805 : cstate%f_pc_nucl, &
806 : cstate%f_growpe, &
807 : cstate%f_evappe, &
808 : cstate%f_evcore, &
809 : cstate%f_growlg, &
810 : cstate%f_evaplg, &
811 : cstate%f_gasprod, &
812 : cstate%f_rlheat, &
813 : cstate%f_radint, &
814 : cstate%f_partheat, &
815 : cstate%f_dtpart, &
816 : cstate%f_cmf, &
817 : cstate%f_totevap, &
818 72960 : stat=ier)
819 72960 : if (ier /= 0) then
820 0 : if (cstate%f_carma%f_do_print) then
821 : write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_Destroy::&
822 0 : &ERROR deallocating growth arrays, status=", ier
823 : end if
824 0 : rc = RC_ERROR
825 0 : return
826 : endif
827 : end if
828 :
829 72960 : if (allocated(cstate%f_pvapl)) then
830 : deallocate( &
831 : cstate%f_pvapl, &
832 : cstate%f_pvapi, &
833 : cstate%f_supsatl, &
834 : cstate%f_supsati, &
835 : cstate%f_supsatlold, &
836 : cstate%f_supsatiold, &
837 72960 : stat=ier)
838 72960 : if (ier /= 0) then
839 0 : if (cstate%f_carma%f_do_print) then
840 : write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_Destroy::&
841 0 : &ERROR deallocating gas arrays, status=", ier
842 : end if
843 0 : rc = RC_ERROR
844 0 : return
845 : endif
846 : end if
847 :
848 72960 : if (allocated(cstate%f_coaglg)) then
849 : deallocate( &
850 : cstate%f_coaglg, &
851 : cstate%f_coagpe, &
852 : cstate%f_ckernel, &
853 72960 : stat = ier)
854 72960 : if (ier /= 0) then
855 0 : if (cstate%f_carma%f_do_print) then
856 : write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_Destroy::&
857 0 : &ERROR deallocating coag arrays, status=", ier
858 : end if
859 0 : rc = RC_ERROR
860 0 : return
861 : end if
862 : end if
863 : end if
864 :
865 : return
866 : end subroutine CARMASTATE_Destroy
867 :
868 :
869 : !! The routine performs the main CARMA processing for one timestep of
870 : !! the parent model. The state variables should have all been set before
871 : !! calling CARMASTATE_Step(). When this routine returns, the state will
872 : !! have been updated to reflect the changes from the CARMA microphysics.
873 : !! If tendencies are desired, then the difference between the final and
874 : !! initial state will need to be computed by the caller.
875 : !!
876 : !! NIOTE: xxxfv, xxxram and xxxfrac need to be specified for dry deposition.
877 : !!
878 : !! @author Chuck Bardeen
879 : !! @version Feb-2009
880 2101248 : subroutine CARMASTATE_Step(cstate, rc, cldfrc, rhcrit, lndfv, ocnfv, icefv, lndram, ocnram, iceram, lndfrac, ocnfrac, icefrac)
881 : type(carmastate_type), intent(inout) :: cstate
882 : integer, intent(out) :: rc
883 : real(kind=f), intent(in), optional :: cldfrc(cstate%f_NZ) !! cloud fraction [fraction]
884 : real(kind=f), intent(in), optional :: rhcrit(cstate%f_NZ) !! relative humidity for onset of liquid clouds [fraction]
885 : real(kind=f), intent(in), optional :: lndfv !! the surface friction velocity over land [m/s]
886 : real(kind=f), intent(in), optional :: ocnfv !! the surface friction velocity over ocean [m/s]
887 : real(kind=f), intent(in), optional :: icefv !! the surface friction velocity over ice [m/s]
888 : real(kind=f), intent(in), optional :: lndram !! the aerodynamic resistance over land [s/m]
889 : real(kind=f), intent(in), optional :: ocnram !! the aerodynamic resistance over ocean [s/m]
890 : real(kind=f), intent(in), optional :: iceram !! the aerodynamic resistance over ice [s/m]
891 : real(kind=f), intent(in), optional :: lndfrac !! land fraction
892 : real(kind=f), intent(in), optional :: ocnfrac !! ocn fraction
893 : real(kind=f), intent(in), optional :: icefrac !! ice fraction
894 :
895 : integer :: iz ! vertical index
896 : integer :: igas ! gas index
897 :
898 : ! Assume success.
899 1050624 : rc = RC_OK
900 :
901 : ! Store the cloud fraction if specified
902 34670592 : cstate%f_cldfrc(:) = 1._f
903 34670592 : cstate%f_rhcrit(:) = 1._f
904 :
905 34670592 : if (present(cldfrc)) cstate%f_cldfrc(:) = cldfrc(:)
906 34670592 : if (present(rhcrit)) cstate%f_rhcrit(:) = rhcrit(:)
907 :
908 : ! Determine the gas supersaturations.
909 34670592 : do iz = 1, cstate%f_NZ
910 101910528 : do igas = 1, cstate%f_carma%f_NGAS
911 67239936 : call supersat(cstate%f_carma, cstate, iz, igas, rc)
912 100859904 : if (rc < 0) return
913 : end do
914 : end do
915 :
916 : ! Determine the particle densities.
917 1050624 : call rhopart(cstate%f_carma, cstate, rc)
918 1050624 : if (rc < 0) return
919 :
920 : ! We have to hold off initialization until now, because the particle density
921 : ! (rhop) can not be determined until the particle masses are known (i.e. after
922 : ! CARMASTATE_SetBin), because rhop is used to determine the fall velocity.
923 : !
924 : ! NOTE: If configured for fixed initialization, then we will lose some accuracy
925 : ! in the calculation of the fall velocities, growth kernels, ... and in return
926 : ! will gain a significant performance by not having to initialize as often.
927 : !
928 : ! NOTE: If configured for partial initialized in conjunction with fixed
929 : ! initialization, then do the fall velocity (and growth) initialization which
930 : ! is relatively quick, but skip the recalculation of the coagulation kernels
931 : ! which is relatively expensive. This could be useful for particles that have
932 : ! a wet radius that is different from the dry radius or where there are large
933 : ! changes from the average conditions (temperature, water vapor, ...) used in
934 : ! the fixed initialization.
935 1050624 : if ((.not. cstate%f_carma%f_do_fixedinit) .or. &
936 : (cstate%f_carma%f_do_partialinit)) then
937 :
938 : ! Initialize the vertical transport.
939 1050624 : if (cstate%f_carma%f_do_vtran .or. cstate%f_carma%f_do_coag .or. cstate%f_carma%f_do_grow) then
940 1050624 : call setupvf(cstate%f_carma, cstate, rc)
941 :
942 1050624 : if (cstate%f_carma%f_do_vdiff) then
943 1050624 : call setupbdif(cstate%f_carma, cstate, rc)
944 : end if
945 : end if
946 :
947 : ! Initialize the nucleation, growth and evaporation.
948 1050624 : if (cstate%f_carma%f_do_grow) then
949 1050624 : call setupgrow(cstate%f_carma, cstate, rc)
950 1050624 : if (rc < RC_OK) return
951 :
952 1050624 : call setupgkern(cstate%f_carma, cstate, rc)
953 1050624 : if (rc < RC_OK) return
954 :
955 1050624 : call setupnuc(cstate%f_carma, cstate, rc)
956 1050624 : if (rc < RC_OK) return
957 : end if
958 :
959 : ! Initialize the coagulation.
960 1050624 : if (cstate%f_carma%f_do_coag .and. &
961 : (.not. cstate%f_carma%f_do_fixedinit)) then
962 1050624 : call setupckern(cstate%f_carma, cstate, rc)
963 1050624 : if (rc < RC_OK) return
964 : end if
965 : end if
966 :
967 : ! Initialize the dry deposition
968 : !
969 : ! NOTE: This is tied to the surface fields that vary from column to column,
970 : ! so it needs to get calculated here whether using fixed or full initialization.
971 1050624 : if (cstate%f_carma%f_do_drydep) then
972 : if (present(lndfv) .and. present(lndram) .and. present(lndfrac) .and. &
973 : present(ocnfv) .and. present(ocnram) .and. present(ocnfrac) .and. &
974 1050624 : present(icefv) .and. present(iceram) .and. present(icefrac)) then
975 :
976 : ! NOTE: Need to convert surfric and ram from mks to cgs units.
977 : call setupvdry(cstate%f_carma, cstate, &
978 : lndfv * 100._f, ocnfv * 100._f, icefv * 100._f, &
979 : lndram / 100._f, ocnram / 100._f, iceram / 100._f, &
980 1050624 : lndfrac, ocnfrac, icefrac, rc)
981 1050624 : if (rc < RC_OK) return
982 : else
983 : write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_Step: &
984 : &do_drydep requires that the optional inputs xxxfv, xxxram &
985 0 : &and xxxfrac be provided."
986 0 : rc = RC_ERROR
987 0 : return
988 : end if
989 : end if
990 :
991 1050624 : call fixcorecol(cstate%f_carma, cstate, rc)
992 :
993 : ! Calculate the impact of microphysics upon the state.
994 1050624 : call step(cstate%f_carma, cstate, rc)
995 :
996 1050624 : call fixcorecol(cstate%f_carma, cstate, rc)
997 :
998 1050624 : return
999 2101248 : end subroutine CARMASTATE_Step
1000 :
1001 :
1002 : ! Query, Control and State I/O
1003 :
1004 : !! Gets the mass mixing ratio for the gas (igas). After a call to CARMA_Step(),
1005 : !! the new mass mixing ratio of the gas can be retrieved.
1006 : !!
1007 : !! @author Chuck Bardeen
1008 : !! @version Feb-2009
1009 : !! @see CARMA_AddGas
1010 : !! @see CARMA_GetGas
1011 : !! @see CARMA_Step
1012 : !! @see CARMASTATE_SetGas
1013 1050624 : subroutine CARMASTATE_Get(cstate, rc, max_nsubstep, max_nretry, nstep, nsubstep, nretry, zsubsteps, xc, yc)
1014 : type(carmastate_type), intent(in) :: cstate !! the carma state object
1015 : integer, intent(out) :: rc !! return code, negative indicates failure
1016 : integer, optional, intent(out) :: max_nsubstep !! maximum number of substeps in a step
1017 : real(kind=f), optional, intent(out) :: max_nretry !! maximum number of retries in a step
1018 : real(kind=f), optional, intent(out) :: nstep !! total number of steps taken
1019 : integer, optional, intent(out) :: nsubstep !! total number of substeps taken
1020 : real(kind=f), optional, intent(out) :: nretry !! total number of retries taken
1021 : real(kind=f), optional, intent(out) :: zsubsteps(cstate%f_NZ) !! number of substeps taken per vertical grid point
1022 : real(kind=f), optional, intent(out) :: xc !! x location at center
1023 : real(kind=f), optional, intent(out) :: yc !! y location at center
1024 :
1025 : ! Assume success.
1026 1123584 : rc = RC_OK
1027 :
1028 1123584 : if (present(max_nsubstep)) max_nsubstep = cstate%f_max_nsubstep
1029 1123584 : if (present(max_nretry)) max_nretry = cstate%f_max_nretry
1030 1123584 : if (present(nstep)) nstep = cstate%f_nstep
1031 1123584 : if (present(nsubstep)) nsubstep = cstate%f_nsubstep
1032 1123584 : if (present(nretry)) nretry = cstate%f_nretry
1033 34743552 : if (present(zsubsteps)) zsubsteps = cstate%f_zsubsteps
1034 1123584 : if (present(xc)) xc = cstate%f_xc
1035 1123584 : if (present(yc)) yc = cstate%f_yc
1036 :
1037 1123584 : return
1038 1123584 : end subroutine CARMASTATE_Get
1039 :
1040 :
1041 : !! Gets the mass of the bins (ibin) for each particle element (ielem). After the
1042 : !! CARMA_Step() call, new particle concentrations are determined. The number density
1043 : !! and the nucleation rate are only calculated if the element is the number density
1044 : !! element for the group.
1045 : !!
1046 : !! @author Chuck Bardeen
1047 : !! @version Feb-2009
1048 : !! @see CARMA_AddElement
1049 : !! @see CARMA_AddGroup
1050 : !! @see CARMA_Step
1051 : !! @see CARMASTATE_SetBin
1052 420249600 : subroutine CARMASTATE_GetBin(cstate, ielem, ibin, mmr, rc, &
1053 924549120 : nmr, numberDensity, nucleationRate, r_wet, rhop_wet, rhop_dry, &
1054 693411840 : surface, sedimentationflux, vf, vd, dtpart, kappa, totalmmr)
1055 : type(carmastate_type), intent(in) :: cstate !! the carma state object
1056 : integer, intent(in) :: ielem !! the element index
1057 : integer, intent(in) :: ibin !! the bin index
1058 : real(kind=f), intent(out) :: mmr(cstate%f_NZ) !! the bin mass mixing ratio [kg/kg]
1059 : integer, intent(out) :: rc !! return code negative indicates failure
1060 : real(kind=f), optional, intent(out) :: nmr(cstate%f_NZ) !! number mixing ratio [#/kg]
1061 : real(kind=f), optional, intent(out) :: numberDensity(cstate%f_NZ) !! number density [#/cm3]
1062 : real(kind=f), optional, intent(out) :: nucleationRate(cstate%f_NZ) !! nucleation rate [1/cm3/s]
1063 : real(kind=f), optional, intent(out) :: r_wet(cstate%f_NZ) !! wet particle radius [cm]
1064 : real(kind=f), optional, intent(out) :: rhop_wet(cstate%f_NZ) !! wet particle density [g/cm3]
1065 : real(kind=f), optional, intent(out) :: rhop_dry(cstate%f_NZ) !! dry particle density [g/cm3]
1066 : real(kind=f), optional, intent(out) :: surface !! particle mass on the surface [kg/m2]
1067 : real(kind=f), optional, intent(out) :: sedimentationflux !! particle sedimentation mass flux to surface [kg/m2/s]
1068 : real(kind=f), optional, intent(out) :: vf(cstate%f_NZ+1) !! fall velocity [cm/s]
1069 : real(kind=f), optional, intent(out) :: vd !! deposition velocity [cm/s]
1070 : real(kind=f), optional, intent(out) :: dtpart(cstate%f_NZ) !! delta particle temperature [K]
1071 : real(kind=f), optional, intent(out) :: kappa(cstate%f_NZ) !! hygroscopicity parameter
1072 : real(kind=f), optional, intent(out) :: totalmmr(cstate%f_NZ) !! mmr of the entire particle (kg/m3)
1073 :
1074 : integer :: ienconc ! index of element that is the particle concentration for the group
1075 : integer :: igroup ! Group containing this bin
1076 : integer :: icore ! coremass index for group
1077 :
1078 : ! Assume success.
1079 420249600 : rc = RC_OK
1080 :
1081 : ! Determine the particle group for the bin.
1082 420249600 : igroup = cstate%f_carma%f_element(ielem)%f_igroup
1083 :
1084 : ! Make sure there are enough elements allocated.
1085 420249600 : if (ielem > cstate%f_carma%f_NELEM) then
1086 0 : if (cstate%f_carma%f_do_print) write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_SetBin:: ERROR - The specifed element (", &
1087 0 : ielem, ") is larger than the number of elements (", cstate%f_carma%f_NELEM, ")."
1088 0 : rc = RC_ERROR
1089 0 : return
1090 : end if
1091 :
1092 : ! Make sure there are enough bins allocated.
1093 420249600 : if (ibin > cstate%f_carma%f_NBIN) then
1094 0 : if (cstate%f_carma%f_do_print) write(cstate%f_carma%f_LUNOPRT, *) "CARMA_SetBin:: ERROR - The specifed bin (", &
1095 0 : ibin, ") is larger than the number of bins (", cstate%f_carma%f_NBIN, ")."
1096 0 : rc = RC_ERROR
1097 0 : return
1098 : end if
1099 :
1100 : ! Use the specified mass mixing ratio and the air density to determine the mass
1101 : ! of the particles in g/x/y/z.
1102 13868236800 : mmr(:) = cstate%f_pc(:, ibin, ielem) / cstate%f_rhoa_wet(:)
1103 :
1104 : ! Handle the special cases for different types of elements ...
1105 420249600 : if ((cstate%f_carma%f_element(ielem)%f_itype == I_INVOLATILE) .or. &
1106 : (cstate%f_carma%f_element(ielem)%f_itype == I_VOLATILE)) then
1107 1386823680 : mmr(:) = mmr(:) * cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1108 378224640 : else if (cstate%f_carma%f_element(ielem)%f_itype == I_CORE2MOM) then
1109 0 : mmr(:) = mmr(:) / cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1110 : end if
1111 :
1112 : ! NOTE: The concentration element will be handled different, we want to return
1113 : ! the mmr of the element NOT the mmr of the total particle, so you need to subtract
1114 : ! the sum of the core masses.
1115 420249600 : ienconc = cstate%f_carma%f_group(igroup)%f_ienconc
1116 :
1117 420249600 : if (ienconc == ielem) then
1118 :
1119 : ! Subtract the core massed from the total mass
1120 231137280 : do icore = 1, cstate%f_carma%f_group(igroup)%f_ncore
1121 6282731520 : mmr(:) = mmr(:) - cstate%f_pc(:, ibin, cstate%f_carma%f_group(igroup)%f_icorelem(icore)) / cstate%f_rhoa_wet(:)
1122 : end do
1123 : end if
1124 :
1125 : ! NOTE: Could add a check here to make sure that the conc_md is greater than or
1126 : ! equal to 0.
1127 13868236800 : where (cstate%f_pc(:, ibin, ienconc) <= 0.0_f) mmr(:) = 0.0_f
1128 :
1129 : ! Do they also want the mass concentration of particles at the surface?
1130 420249600 : if (present(surface)) then
1131 :
1132 : ! Convert from g/cm2 to kg/m2
1133 0 : surface = cstate%f_pc_surf(ibin, ielem) * 1e4_f / 1e3_f
1134 :
1135 : ! Handle the special cases for different types of elements ...
1136 0 : if ((cstate%f_carma%f_element(ielem)%f_itype == I_INVOLATILE) .or. &
1137 : (cstate%f_carma%f_element(ielem)%f_itype == I_VOLATILE)) then
1138 0 : surface = surface * cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1139 0 : else if (cstate%f_carma%f_element(ielem)%f_itype == I_CORE2MOM) then
1140 0 : surface = surface / cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1141 : end if
1142 :
1143 : ! NOTE: The concentration element will be handled different, we want to return
1144 : ! the mass concentration of the element NOT the mass concentration of the total
1145 : ! particle, so you need to subtract the sum of the core masses.
1146 0 : if (ienconc == ielem) then
1147 :
1148 : ! Subtract the core massed from the total mass
1149 0 : do icore = 1, cstate%f_carma%f_group(igroup)%f_ncore
1150 0 : surface = surface - cstate%f_pc_surf(ibin, cstate%f_carma%f_group(igroup)%f_icorelem(icore)) * 1e4_f / 1e3_f
1151 : end do
1152 : end if
1153 : end if
1154 :
1155 : ! Do they also want the mass flux of particles that sediment to the surface?
1156 420249600 : if (present(sedimentationflux)) then
1157 :
1158 : ! Convert from g/cm2 to kg/m2
1159 231137280 : sedimentationflux = cstate%f_sedimentationflux(ibin, ielem) * 1e4_f / 1e3_f
1160 :
1161 : ! Handle the special cases for different types of elements ...
1162 231137280 : if ((cstate%f_carma%f_element(ielem)%f_itype == I_INVOLATILE) .or. &
1163 : (cstate%f_carma%f_element(ielem)%f_itype == I_VOLATILE)) then
1164 42024960 : sedimentationflux = sedimentationflux * cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1165 189112320 : else if (cstate%f_carma%f_element(ielem)%f_itype == I_CORE2MOM) then
1166 0 : sedimentationflux = sedimentationflux / cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1167 : end if
1168 :
1169 : ! NOTE: The concentration element will be handled different, we want to return
1170 : ! the sedimentation flux of the element NOT the sedimentation flux of the total
1171 : ! particle, so you need to subtract the sum of the core masses.
1172 231137280 : if (ienconc == ielem) then
1173 :
1174 : ! Subtract the core massed from the total mass
1175 231137280 : do icore = 1, cstate%f_carma%f_group(igroup)%f_ncore
1176 231137280 : sedimentationflux = sedimentationflux - cstate%f_sedimentationflux(ibin, cstate%f_carma%f_group(igroup)%f_icorelem(icore)) * 1e4_f / 1e3_f
1177 : end do
1178 : end if
1179 : end if
1180 :
1181 : ! Is the hygroscopicity parameter requested?
1182 420249600 : if (present(kappa)) then
1183 0 : kappa = cstate%f_kappahygro(:, ibin, igroup)
1184 : end if
1185 :
1186 :
1187 : ! If this is the partcile # element, then determine some other statistics.
1188 420249600 : if (ienconc == ielem) then
1189 1386823680 : if (present(totalmmr)) totalmmr(:) = (cstate%f_pc(:, ibin, ielem) * cstate%f_carma%f_group(igroup)%f_rmass(ibin)) / cstate%f_rhoa_wet(:)
1190 42024960 : if (present(nmr)) nmr(:) = (cstate%f_pc(:, ibin, ielem) / cstate%f_rhoa_wet(:)) * 1000._f
1191 1386823680 : if (present(numberDensity)) numberDensity(:) = cstate%f_pc(:, ibin, ielem) / cstate%f_zmet(:)
1192 1386823680 : if (present(r_wet)) r_wet(:) = cstate%f_r_wet(:, ibin, igroup)
1193 1386823680 : if (present(rhop_wet)) rhop_wet(:) = cstate%f_rhop_wet(:, ibin, igroup)
1194 42024960 : if (present(rhop_dry)) rhop_dry(:) = cstate%f_rhop(:, ibin, igroup)
1195 :
1196 42024960 : if (cstate%f_carma%f_do_vtran) then
1197 1428848640 : if (present(vf)) vf(:) = cstate%f_vf(:, ibin, igroup) * cstate%f_zmetl(:)
1198 : else
1199 0 : if (present(vf)) vf(:) = CAM_FILL
1200 : end if
1201 :
1202 42024960 : if (cstate%f_carma%f_do_drydep) then
1203 42024960 : if (present(vd)) then
1204 42024960 : if (cstate%f_igridv .eq. I_CART) then
1205 0 : vd = cstate%f_vd(ibin, igroup) * cstate%f_zmetl(1)
1206 : else
1207 42024960 : vd = cstate%f_vd(ibin, igroup) * cstate%f_zmetl(cstate%f_NZP1)
1208 : end if
1209 : end if
1210 : else
1211 0 : if (present(vd)) vd = CAM_FILL
1212 : end if
1213 :
1214 42024960 : if (cstate%f_carma%f_do_grow) then
1215 42024960 : if (present(nucleationRate)) nucleationRate(:) = cstate%f_pc_nucl(:, ibin, ielem) / &
1216 1386823680 : cstate%f_zmet(:) / cstate%f_dtime
1217 : else
1218 0 : if (present(nucleationRate)) nucleationRate(:) = CAM_FILL
1219 : end if
1220 :
1221 42024960 : if (cstate%f_carma%f_do_pheat) then
1222 0 : if (present(dtpart)) dtpart(:) = cstate%f_dtpart(:, ibin, igroup)
1223 : else
1224 1386823680 : if (present(dtpart)) dtpart(:) = CAM_FILL
1225 : end if
1226 : else
1227 6429818880 : if (present(totalmmr)) totalmmr(:) = CAM_FILL
1228 378224640 : if (present(nmr)) nmr(:) = CAM_FILL
1229 6429818880 : if (present(numberDensity)) numberDensity(:) = CAM_FILL
1230 6429818880 : if (present(nucleationRate)) nucleationRate(:) = CAM_FILL
1231 6429818880 : if (present(r_wet)) r_wet(:) = CAM_FILL
1232 6429818880 : if (present(rhop_wet)) rhop_wet(:) = CAM_FILL
1233 378224640 : if (present(rhop_dry)) rhop_dry(:) = CAM_FILL
1234 6429818880 : if (present(dtpart)) dtpart(:) = CAM_FILL
1235 6618931200 : if (present(vf)) vf(:) = CAM_FILL
1236 378224640 : if (present(vd)) vd = CAM_FILL
1237 : end if
1238 :
1239 : return
1240 1807073280 : end subroutine CARMASTATE_GetBin
1241 :
1242 :
1243 : !! Gets the mass of the detrained condensate for the bins (ibin) for each particle
1244 : !! element (ielem) in the grid.
1245 : !!
1246 : !!
1247 : !! @author Chuck Bardeen
1248 : !! @version Feb-2009
1249 : !! @see CARMA_AddElement
1250 : !! @see CARMA_AddGroup
1251 : !! @see CARMA_Step
1252 : !! @see CARMASTATE_SetDetrain
1253 0 : subroutine CARMASTATE_GetDetrain(cstate, ielem, ibin, mmr, rc, nmr, numberDensity, r_wet, rhop_wet)
1254 : type(carmastate_type), intent(in) :: cstate !! the carma state object
1255 : integer, intent(in) :: ielem !! the element index
1256 : integer, intent(in) :: ibin !! the bin index
1257 : real(kind=f), intent(out) :: mmr(cstate%f_NZ) !! the bin mass mixing ratio [kg/kg]
1258 : integer, intent(out) :: rc !! return code negative indicates failure
1259 : real(kind=f), optional, intent(out) :: nmr(cstate%f_NZ) !! number mixing ratio [#/kg]
1260 : real(kind=f), optional, intent(out) :: numberDensity(cstate%f_NZ) !! number density [#/cm3]
1261 : real(kind=f), optional, intent(out) :: r_wet(cstate%f_NZ) !! wet particle radius [cm]
1262 : real(kind=f), optional, intent(out) :: rhop_wet(cstate%f_NZ) !! wet particle density [g/cm3]
1263 :
1264 : integer :: ienconc !! index of element that is the particle concentration for the group
1265 : integer :: igroup ! Group containing this bin
1266 :
1267 : ! Assume success.
1268 0 : rc = RC_OK
1269 :
1270 : ! Determine the particle group for the bin.
1271 0 : igroup = cstate%f_carma%f_element(ielem)%f_igroup
1272 :
1273 : ! Make sure there are enough elements allocated.
1274 0 : if (ielem > cstate%f_carma%f_NELEM) then
1275 0 : if (cstate%f_carma%f_do_print) write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_SetDetrain:: ERROR - The specifed element (", &
1276 0 : ielem, ") is larger than the number of elements (", cstate%f_carma%f_NELEM, ")."
1277 0 : rc = RC_ERROR
1278 0 : return
1279 : end if
1280 :
1281 : ! Make sure there are enough bins allocated.
1282 0 : if (ibin > cstate%f_carma%f_NBIN) then
1283 0 : if (cstate%f_carma%f_do_print) write(cstate%f_carma%f_LUNOPRT, *) "CARMA_SetDetrainin:: ERROR - The specifed bin (", &
1284 0 : ibin, ") is larger than the number of bins (", cstate%f_carma%f_NBIN, ")."
1285 0 : rc = RC_ERROR
1286 0 : return
1287 : end if
1288 :
1289 :
1290 : ! Use the specified mass mixing ratio and the air density to determine the mass
1291 : ! of the particles in g/x/y/z.
1292 0 : mmr(:) = cstate%f_pcd(:, ibin, ielem) / cstate%f_rhoa_wet(:)
1293 :
1294 :
1295 : ! Handle the special cases for different types of elements ...
1296 0 : if ((cstate%f_carma%f_element(ielem)%f_itype == I_INVOLATILE) .or. &
1297 : (cstate%f_carma%f_element(ielem)%f_itype == I_VOLATILE)) then
1298 0 : mmr(:) = mmr(:) * cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1299 0 : else if (cstate%f_carma%f_element(ielem)%f_itype == I_CORE2MOM) then
1300 0 : mmr(:) = mmr(:) / cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1301 : end if
1302 :
1303 : ! If this is the partcile # element, then determine some other statistics.
1304 0 : ienconc = cstate%f_carma%f_group(igroup)%f_ienconc
1305 0 : if (ienconc == ielem) then
1306 0 : if (present(nmr)) nmr(:) = (cstate%f_pcd(:, ibin, ielem) / cstate%f_rhoa_wet(:)) * 1000._f
1307 0 : if (present(numberDensity)) numberDensity(:) = cstate%f_pcd(:, ibin, ielem) / cstate%f_zmet(:)
1308 0 : if (present(r_wet)) r_wet(:) = cstate%f_r_wet(:, ibin, igroup)
1309 0 : if (present(rhop_wet)) rhop_wet(:) = cstate%f_rhop_wet(:, ibin, igroup)
1310 : else
1311 0 : if (present(nmr)) nmr(:) = CAM_FILL
1312 0 : if (present(numberDensity)) numberDensity(:) = CAM_FILL
1313 : end if
1314 :
1315 : return
1316 0 : end subroutine CARMASTATE_GetDetrain
1317 :
1318 :
1319 : !! Gets the mass mixing ratio for the gas (igas). After a call to CARMA_Step(),
1320 : !! the new mass mixing ratio of the gas can be retrieved.
1321 : !!
1322 : !! @author Chuck Bardeen
1323 : !! @version Feb-2009
1324 : !! @see CARMA_AddGas
1325 : !! @see CARMA_GetGas
1326 : !! @see CARMA_Step
1327 : !! @see CARMASTATE_SetGas
1328 10506240 : subroutine CARMASTATE_GetGas(cstate, igas, mmr, rc, satice, satliq, eqice, eqliq, wtpct)
1329 : type(carmastate_type), intent(in) :: cstate !! the carma state object
1330 : integer, intent(in) :: igas !! the gas index
1331 : real(kind=f), intent(out) :: mmr(cstate%f_NZ) !! the gas mass mixing ratio [kg/kg]
1332 : integer, intent(out) :: rc !! return code, negative indicates failure
1333 : real(kind=f), optional, intent(out) :: satice(cstate%f_NZ) !! the gas saturation wrt ice
1334 : real(kind=f), optional, intent(out) :: satliq(cstate%f_NZ) !! the gas saturation wrt liquid
1335 : real(kind=f), optional, intent(out) :: eqice(cstate%f_NZ) !! the gas vapor pressure wrt ice
1336 : real(kind=f), optional, intent(out) :: eqliq(cstate%f_NZ) !! the gas vapor pressure wrt liquid
1337 : real(kind=f), optional, intent(out) :: wtpct(cstate%f_NZ) !! weight percent aerosol composition
1338 :
1339 : ! Assume success.
1340 2101248 : rc = RC_OK
1341 :
1342 : ! Make sure there are enough gases allocated.
1343 2101248 : if (igas > cstate%f_carma%f_NGAS) then
1344 0 : if (cstate%f_carma%f_do_print) write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_GetGas:: ERROR - The specifed gas (", &
1345 0 : igas, ") is larger than the number of gases (", cstate%f_carma%f_NGAS, ")."
1346 0 : rc = RC_ERROR
1347 0 : return
1348 : end if
1349 :
1350 : ! Use the specified mass mixing ratio and the air density to determine the mass
1351 : ! of the gas in g/x/y/z.
1352 69341184 : mmr(:) = cstate%f_gc(:, igas) / cstate%f_rhoa_wet(:)
1353 :
1354 69341184 : if (present(satice)) satice(:) = cstate%f_supsati(:, igas) + 1._f
1355 69341184 : if (present(satliq)) satliq(:) = cstate%f_supsatl(:, igas) + 1._f
1356 69341184 : if (present(eqice)) eqice(:) = cstate%f_pvapi(:, igas) / cstate%f_p(:)
1357 69341184 : if (present(eqliq)) eqliq(:) = cstate%f_pvapl(:, igas) / cstate%f_p(:)
1358 69341184 : if (present(wtpct)) wtpct(:) = cstate%f_wtpct(:)
1359 :
1360 : return
1361 10506240 : end subroutine CARMASTATE_GetGas
1362 :
1363 :
1364 : !! Gets information about the state of the atmosphere. After the CARMA_Step() call,
1365 : !! a new atmospheric state is determined.
1366 : !!
1367 : !! @author Chuck Bardeen
1368 : !! @version Feb-2009
1369 : !! @see CARMA_Step
1370 : !! @see CARMASTATE_Create
1371 1050624 : subroutine CARMASTATE_GetState(cstate, rc, t, p, rhoa_wet, rlheat)
1372 : type(carmastate_type), intent(in) :: cstate !! the carma state object
1373 : integer, intent(out) :: rc !! return code, negative indicates failure
1374 : real(kind=f), optional, intent(out) :: t(cstate%f_NZ) !! the air temperature [K]
1375 : real(kind=f), optional, intent(out) :: p(cstate%f_NZ) !! the air pressure [Pa]
1376 : real(kind=f), optional, intent(out) :: rhoa_wet(cstate%f_NZ) !! air density [kg m-3]
1377 : real(kind=f), optional, intent(out) :: rlheat(cstate%f_NZ) !! latent heat [K/s]
1378 :
1379 : ! Assume success.
1380 1050624 : rc = RC_OK
1381 :
1382 : ! Return the temperature, pressure, and/or density.
1383 34670592 : if (present(t)) t(:) = cstate%f_t(:)
1384 :
1385 : ! DYNE -> Pa
1386 1050624 : if (present(p)) p(:) = cstate%f_p(:) / RPA2CGS
1387 :
1388 : ! Convert rhoa from the scaled units to mks.
1389 1050624 : if (present(rhoa_wet)) rhoa_wet(:) = (cstate%f_rhoa_wet(:) / cstate%f_zmet(:)) * 1e6_f / 1e3_f
1390 :
1391 1050624 : if (present(rlheat)) rlheat(:) = cstate%f_rlheat(:)
1392 :
1393 1050624 : return
1394 2101248 : end subroutine CARMASTATE_GetState
1395 :
1396 :
1397 : !! Sets the mass of the bins (ibin) for each particle element (ielem) in the grid.
1398 : !! This call should be made after CARMASTATE_Create() and before CARMA_Step().
1399 : !!
1400 : !! @author Chuck Bardeen
1401 : !! @version Feb-2009
1402 : !! @see CARMA_AddBin
1403 : !! @see CARMA_Step
1404 : !! @see CARMASTATE_GetBin
1405 336199680 : subroutine CARMASTATE_SetBin(cstate, ielem, ibin, mmr, rc, surface)
1406 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
1407 : integer, intent(in) :: ielem !! the element index
1408 : integer, intent(in) :: ibin !! the bin index
1409 : real(kind=f), intent(in) :: mmr(cstate%f_NZ) !! the bin mass mixing ratio [kg/kg]
1410 : integer, intent(out) :: rc !! return code, negative indicates failure
1411 : real(kind=f), optional, intent(in) :: surface !! element mass on the surface [kg/m2]
1412 :
1413 : integer :: igroup ! Group containing this bin
1414 672399360 : real(kind=f) :: conc_md(cstate%f_NZ) ! mass density of concentration element (g/x/y/z)
1415 : real(kind=f) :: conc_mc ! mass concentration of concentration element (kg/m2)
1416 : integer :: ienconc
1417 : integer :: icore
1418 :
1419 : ! Assume success.
1420 336199680 : rc = RC_OK
1421 :
1422 : ! Determine the particle group for the bin.
1423 336199680 : igroup = cstate%f_carma%f_element(ielem)%f_igroup
1424 :
1425 : ! Make sure there are enough elements allocated.
1426 336199680 : if (ielem > cstate%f_carma%f_NELEM) then
1427 0 : if (cstate%f_carma%f_do_print) write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_SetBin:: ERROR - The specifed element (", &
1428 0 : ielem, ") is larger than the number of elements (", cstate%f_carma%f_NELEM, ")."
1429 0 : rc = RC_ERROR
1430 0 : return
1431 : end if
1432 :
1433 : ! Make sure there are enough bins allocated.
1434 336199680 : if (ibin > cstate%f_carma%f_NBIN) then
1435 0 : if (cstate%f_carma%f_do_print) write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_SetBin:: ERROR - The specifed bin (", &
1436 0 : ibin, ") is larger than the number of bins (", cstate%f_carma%f_NBIN, ")."
1437 0 : rc = RC_ERROR
1438 0 : return
1439 : end if
1440 :
1441 336199680 : ienconc = cstate%f_carma%f_group(igroup)%f_ienconc
1442 :
1443 : ! Determine the current mass density of the species represented by the concentration element
1444 : ! (concentration_element*rmass - sum(coremass))
1445 11094589440 : conc_md(:) = cstate%f_pc(:, ibin, ienconc) * cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1446 :
1447 3172884480 : do icore = 1, cstate%f_carma%f_group(igroup)%f_ncore
1448 93946798080 : conc_md(:) = conc_md(:) - cstate%f_pc(:, ibin, cstate%f_carma%f_group(igroup)%f_icorelem(icore))
1449 : end do
1450 :
1451 : ! Use the specified mass mixing ratio and the air density to determine the mass
1452 : ! of the particles in g/x/y/z.
1453 11094589440 : cstate%f_pc(:, ibin, ielem) = mmr(:) * cstate%f_rhoa_wet(:)
1454 :
1455 : ! Handle the special cases for different types of elements ...
1456 336199680 : if ((cstate%f_carma%f_element(ielem)%f_itype == I_INVOLATILE) .or. &
1457 : (cstate%f_carma%f_element(ielem)%f_itype == I_VOLATILE)) then
1458 1386823680 : cstate%f_pc(:, ibin, ielem) = cstate%f_pc(:, ibin, ielem) / cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1459 294174720 : else if (cstate%f_carma%f_element(ielem)%f_itype == I_CORE2MOM) then
1460 0 : cstate%f_pc(:, ibin, ielem) = cstate%f_pc(:, ibin, ielem) * cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1461 : end if
1462 :
1463 : ! Update the concentration element to be the concentration of the particle (i.e. mass density
1464 : ! of the concentration element species plus the sum of the coremasses).
1465 : !
1466 : ! NOTE: Recalculating this from scratch makes sure that the concentration element * rmass is
1467 : ! greater than or equal to the sum of the cores masses and you don't end up with a
1468 : ! negative concentration species mass.
1469 336199680 : if (ielem /= ienconc) then
1470 9707765760 : cstate%f_pc(:, ibin, ienconc) = conc_md(:) / cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1471 : end if
1472 :
1473 3172884480 : do icore = 1, cstate%f_carma%f_group(igroup)%f_ncore
1474 93946798080 : cstate%f_pc(:, ibin, ienconc) = cstate%f_pc(:, ibin, ienconc) + cstate%f_pc(:, ibin, cstate%f_carma%f_group(igroup)%f_icorelem(icore)) / cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1475 : end do
1476 :
1477 : ! If they specified an initial mass of particles on the surface, then use that
1478 : ! value.
1479 336199680 : if (present(surface)) then
1480 :
1481 : ! Determine the current mass density of the species represented by the concentration element
1482 : ! (concentration_element*rmass - sum(coremass))
1483 0 : conc_mc = cstate%f_pc_surf(ibin, ienconc) * cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1484 :
1485 0 : do icore = 1, cstate%f_carma%f_group(igroup)%f_ncore
1486 0 : conc_mc = conc_mc - cstate%f_pc_surf(ibin, cstate%f_carma%f_group(igroup)%f_icorelem(icore))
1487 : end do
1488 :
1489 : ! Convert from g/cm2 to kg/m2
1490 0 : cstate%f_pc_surf(ibin, ielem) = surface / 1e4_f * 1e3_f
1491 :
1492 : ! Handle the special cases for different types of elements ...
1493 0 : if ((cstate%f_carma%f_element(ielem)%f_itype == I_INVOLATILE) .or. &
1494 : (cstate%f_carma%f_element(ielem)%f_itype == I_VOLATILE)) then
1495 0 : cstate%f_pc_surf(ibin, ielem) = cstate%f_pc_surf(ibin, ielem) / cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1496 0 : else if (cstate%f_carma%f_element(ielem)%f_itype == I_CORE2MOM) then
1497 0 : cstate%f_pc_surf(ibin, ielem) = cstate%f_pc_surf(ibin, ielem) * cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1498 : end if
1499 :
1500 : ! Update the concentration element to be the concentration of the particle (i.e. mass density
1501 : ! of the concentration element species plus the sum of the coremasses).
1502 : !
1503 : ! NOTE: Recalculating this from scratch makes sure that the concentration element * rmass is
1504 : ! greater than or equal to the sum of the cores masses and you don't end up with a
1505 : ! negative concentration species mass.
1506 0 : if (ielem /= ienconc) then
1507 0 : cstate%f_pc_surf(ibin, ienconc) = conc_mc / cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1508 : end if
1509 :
1510 0 : do icore = 1, cstate%f_carma%f_group(igroup)%f_ncore
1511 0 : cstate%f_pc_surf(ibin, ienconc) = cstate%f_pc_surf(ibin, ienconc) + cstate%f_pc_surf(ibin, cstate%f_carma%f_group(igroup)%f_icorelem(icore)) / cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1512 : end do
1513 : else
1514 336199680 : cstate%f_pc_surf(ibin, ielem) = 0.0_f
1515 : end if
1516 :
1517 : return
1518 : end subroutine CARMASTATE_SetBin
1519 :
1520 :
1521 : !! Sets the mass of the detrained condensate for the bins (ibin) for each particle
1522 : !! element (ielem) in the grid. This call should be made after CARMASTATE_Create()
1523 : !! and before CARMA_Step().
1524 : !!
1525 : !! @author Chuck Bardeen
1526 : !! @version May-2010
1527 : !! @see CARMA_AddBin
1528 : !! @see CARMA_Step
1529 : !! @see CARMASTATE_GetDetrain
1530 0 : subroutine CARMASTATE_SetDetrain(cstate, ielem, ibin, mmr, rc)
1531 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
1532 : integer, intent(in) :: ielem !! the element index
1533 : integer, intent(in) :: ibin !! the bin index
1534 : real(kind=f), intent(in) :: mmr(cstate%f_NZ) !! the bin mass mixing ratio [kg/kg]
1535 : integer, intent(out) :: rc !! return code, negative indicates failure
1536 :
1537 : integer :: igroup ! Group containing this bin
1538 :
1539 : ! Assume success.
1540 0 : rc = RC_OK
1541 :
1542 : ! Determine the particle group for the bin.
1543 0 : igroup = cstate%f_carma%f_element(ielem)%f_igroup
1544 :
1545 : ! Make sure there are enough elements allocated.
1546 0 : if (ielem > cstate%f_carma%f_NELEM) then
1547 0 : if (cstate%f_carma%f_do_print) write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_SetDetrain:: ERROR - The specifed element (", &
1548 0 : ielem, ") is larger than the number of elements (", cstate%f_carma%f_NELEM, ")."
1549 0 : rc = RC_ERROR
1550 0 : return
1551 : end if
1552 :
1553 : ! Make sure there are enough bins allocated.
1554 0 : if (ibin > cstate%f_carma%f_NBIN) then
1555 0 : if (cstate%f_carma%f_do_print) write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_SetDetrain:: ERROR - The specifed bin (", &
1556 0 : ibin, ") is larger than the number of bins (", cstate%f_carma%f_NBIN, ")."
1557 0 : rc = RC_ERROR
1558 0 : return
1559 : end if
1560 :
1561 : ! Use the specified mass mixing ratio and the air density to determine the mass
1562 : ! of the particles in g/x/y/z.
1563 0 : cstate%f_pcd(:, ibin, ielem) = mmr(:) * cstate%f_rhoa_wet(:)
1564 :
1565 : ! Handle the special cases for different types of elements ...
1566 0 : if ((cstate%f_carma%f_element(ielem)%f_itype == I_INVOLATILE) .or. &
1567 : (cstate%f_carma%f_element(ielem)%f_itype == I_VOLATILE)) then
1568 0 : cstate%f_pcd(:, ibin, ielem) = cstate%f_pcd(:, ibin, ielem) / cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1569 0 : else if (cstate%f_carma%f_element(ielem)%f_itype == I_CORE2MOM) then
1570 0 : cstate%f_pcd(:, ibin, ielem) = cstate%f_pcd(:, ibin, ielem) * cstate%f_carma%f_group(igroup)%f_rmass(ibin)
1571 : end if
1572 :
1573 : return
1574 : end subroutine CARMASTATE_SetDetrain
1575 :
1576 :
1577 :
1578 : !! Sets the mass of the gas (igas) in the grid. This call should be made after
1579 : !! CARMASTATE_Create() and before CARMA_Step().
1580 : !!
1581 : !! @author Chuck Bardeen
1582 : !! @version Feb-2009
1583 : !! @see CARMA_AddGas
1584 : !! @see CARMA_GetGas
1585 : !! @see CARMA_InitializeStep
1586 : !! @see CARMA_Step
1587 6303744 : subroutine CARMASTATE_SetGas(cstate, igas, mmr, rc, mmr_old, satice_old, satliq_old)
1588 : type(carmastate_type), intent(inout) :: cstate !! the carma object
1589 : integer, intent(in) :: igas !! the gas index
1590 : real(kind=f), intent(in) :: mmr(cstate%f_NZ) !! the gas mass mixing ratio [kg/kg]
1591 : integer, intent(out) :: rc !! return code, negative indicates failure
1592 : real(kind=f), intent(in), optional :: mmr_old(cstate%f_NZ) !! the previous gas mass mixing ratio [kg/kg]
1593 : real(kind=f), intent(inout), optional :: satice_old(cstate%f_NZ) !! the previous gas saturation wrt ice, calculates if -1
1594 : real(kind=f), intent(inout), optional :: satliq_old(cstate%f_NZ) !! the previous gas saturation wrt liquid, calculates if -1
1595 :
1596 4202496 : real(kind=f) :: tnew(cstate%f_NZ)
1597 : integer :: iz
1598 : logical :: calculateOld
1599 :
1600 : ! Assume success.
1601 2101248 : rc = RC_OK
1602 :
1603 : ! Make sure there are enough gases allocated.
1604 2101248 : if (igas > cstate%f_carma%f_NGAS) then
1605 0 : if (cstate%f_carma%f_do_print) write(cstate%f_carma%f_LUNOPRT, *) "CARMASTATE_SetGas:: ERROR - The specifed gas (", &
1606 0 : igas, ") is larger than the number of gases (", cstate%f_carma%f_NGAS, ")."
1607 0 : rc = RC_ERROR
1608 0 : return
1609 : end if
1610 :
1611 2101248 : if (cstate%f_carma%f_do_substep) then
1612 2101248 : if (.not. present(mmr_old)) then
1613 0 : if (cstate%f_carma%f_do_print) then
1614 : write(cstate%f_carma%f_LUNOPRT,*) "CARMASTATE_SetGas: &
1615 0 : &Error - Need to specify mmr_old, satic_old, satliq_old when substepping."
1616 : end if
1617 0 : rc = RC_ERROR
1618 :
1619 0 : return
1620 :
1621 : else
1622 69341184 : cstate%f_gcl(:, igas) = mmr_old(:) * cstate%f_rhoa_wet(:) * cstate%f_t(:) / cstate%f_told(:)
1623 :
1624 : ! A value of -1 for the saturation ratio means that it needs to be calculated from the old temperature
1625 : ! and the old gc.
1626 : !
1627 : ! NOTE: This is typically just a problem for the first step, so we just need to get close.
1628 2101248 : calculateOld = .false.
1629 2101248 : if (present(satice_old) .and. present(satliq_old)) then
1630 131604480 : if (any(satice_old(:) == -1._f) .or. any(satliq_old(:) == -1._f)) calculateOld = .true.
1631 : else
1632 : calculateOld = .true.
1633 : end if
1634 :
1635 : if (calculateOld) then
1636 :
1637 : ! This is a bit of a hack, because of the way CARMA has the vapor pressure and saturation
1638 : ! routines implemented.
1639 :
1640 : ! Temporarily set the temperature and gc of to the old state
1641 :
1642 3649536 : tnew(:) = cstate%f_t(:)
1643 3649536 : cstate%f_t(:) = cstate%f_told(:)
1644 :
1645 3649536 : cstate%f_gc(:, igas) = mmr_old(:) * cstate%f_rhoa_wet(:)
1646 :
1647 3649536 : do iz = 1, cstate%f_NZ
1648 3538944 : call supersat(cstate%f_carma, cstate, iz, igas, rc)
1649 3538944 : if (rc < RC_OK) return
1650 :
1651 3538944 : if (present(satice_old)) then
1652 3538944 : if (satice_old(iz) == -1._f) then
1653 3538944 : cstate%f_supsatiold(iz, igas) = cstate%f_supsati(iz, igas)
1654 : else
1655 0 : cstate%f_supsatiold(iz, igas) = satice_old(iz) - 1._f
1656 : endif
1657 : else
1658 0 : cstate%f_supsatiold(iz, igas) = cstate%f_supsati(iz, igas)
1659 : end if
1660 :
1661 3649536 : if (present(satliq_old)) then
1662 3538944 : if (satliq_old(iz) == -1._f) then
1663 3538944 : cstate%f_supsatlold(iz, igas) = cstate%f_supsatl(iz, igas)
1664 : else
1665 0 : cstate%f_supsatlold(iz, igas) = satliq_old(iz) - 1._f
1666 : endif
1667 : else
1668 0 : cstate%f_supsatlold(iz, igas) = cstate%f_supsatl(iz, igas)
1669 : end if
1670 : end do
1671 :
1672 3649536 : cstate%f_t(:) = tnew(:)
1673 :
1674 : else
1675 67682304 : cstate%f_supsatiold(:, igas) = satice_old(:) - 1._f
1676 65691648 : cstate%f_supsatlold(:, igas) = satliq_old(:) - 1._f
1677 : end if
1678 : end if
1679 : end if
1680 :
1681 : ! Use the specified mass mixing ratio and the air density to determine the mass
1682 : ! of the gas in g/x/y/z.
1683 69341184 : cstate%f_gc(:, igas) = mmr(:) * cstate%f_rhoa_wet(:)
1684 :
1685 : return
1686 6303744 : end subroutine CARMASTATE_SetGas
1687 :
1688 :
1689 : !! Sets information about the state of the atmosphere.
1690 : !!
1691 : !! @author Chuck Bardeen
1692 : !! @version Feb-2009
1693 : !! @see CARMA_Step
1694 : !! @see CARMASTATE_Create
1695 0 : subroutine CARMASTATE_SetState(cstate, rc, t, rhoa_wet)
1696 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
1697 : integer, intent(out) :: rc !! return code, negative indicates failure
1698 : real(kind=f), optional, intent(in) :: t(cstate%f_NZ) !! the air temperature [K]
1699 : real(kind=f), optional, intent(in) :: rhoa_wet(cstate%f_NZ) !! air density [kg m-3]
1700 :
1701 : ! Assume success.
1702 0 : rc = RC_OK
1703 :
1704 : ! Return the temperature or density.
1705 0 : if (present(t)) cstate%f_t(:) = t(:)
1706 :
1707 : ! Convert rhoa from mks to the scaled units.
1708 0 : if (present(rhoa_wet)) cstate%f_rhoa_wet(:) = (rhoa_wet(:) * cstate%f_zmet(:)) / 1e6_f * 1e3_f
1709 :
1710 0 : return
1711 0 : end subroutine CARMASTATE_SetState
1712 :
1713 : end module carmastate_mod
|
