Line data Source code
1 : ! Include shortname defintions, so that the F77 code does not have to be modified to
2 : ! reference the CARMA structure.
3 : #include "carma_globaer.h"
4 :
5 : !! This routine manages the calculations that update state variables
6 : !! of the model with new values at the current simulation time.
7 : !!
8 : !! @author Bardeen
9 : !! @version Jan 2012
10 1050624 : subroutine newstate(carma, cstate, rc)
11 :
12 : ! types
13 : use carma_precision_mod
14 : use carma_enums_mod
15 : use carma_constants_mod
16 : use carma_types_mod
17 : use carmastate_mod
18 : use carma_mod
19 :
20 : implicit none
21 :
22 : type(carma_type), intent(in) :: carma !! the carma object
23 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
24 : integer, intent(inout) :: rc !! return code, negative indicates failure
25 :
26 2101248 : real(kind=f) :: pc_orig(NZ,NBIN,NELEM)
27 2101248 : real(kind=f) :: gc_orig(NZ,NGAS)
28 2101248 : real(kind=f) :: t_orig(NZ)
29 2101248 : real(kind=f) :: cldfrc_orig(NZ)
30 2101248 : real(kind=f) :: scale_cldfrc(NZ)
31 2101248 : real(kind=f) :: pc_cloudy(NZ,NBIN,NELEM)
32 2101248 : real(kind=f) :: gc_cloudy(NZ,NGAS)
33 2101248 : real(kind=f) :: t_cloudy(NZ)
34 2101248 : real(kind=f) :: rlheat_cloudy(NZ)
35 2101248 : real(kind=f) :: partheat_cloudy(NZ)
36 2101248 : real(kind=f) :: zsubsteps_cloudy(NZ)
37 2101248 : real(kind=f) :: pc_clear(NZ,NBIN,NELEM)
38 2101248 : real(kind=f) :: gc_clear(NZ,NGAS)
39 2101248 : real(kind=f) :: t_clear(NZ)
40 2101248 : real(kind=f) :: rlheat_clear(NZ)
41 2101248 : real(kind=f) :: partheat_clear(NZ)
42 2101248 : real(kind=f) :: zsubsteps_clear(NZ)
43 2101248 : real(kind=f) :: scale_threshold(NZ)
44 : integer :: igroup
45 : integer :: igas
46 : integer :: ielem
47 : integer :: ibin
48 : integer :: iz
49 :
50 : ! Calculate changes due to vertical transport
51 1050624 : if (do_vtran) then
52 :
53 1050624 : call vertical(carma, cstate, rc)
54 1050624 : if (rc < RC_OK) return
55 : endif
56 :
57 1050624 : call fixcorecol(carma, cstate, rc)
58 :
59 : ! There can be two phases to the microphysics: in-cloud and clear sky. Particles
60 : ! that are tagged as "In-cloud" will only be processed in the in-cloud loop, and their
61 : ! concentrations will be scaled by the cloud fraction since it is assumed to be all
62 : ! in-cloud. Other particle types will be process in-cloud and out of cloud; however,
63 : ! their mass is assumed to be a gridbox average.
64 :
65 : ! If doing doing in-cloud processing, then scale the parameters for in-cloud concentrations.
66 : !
67 : ! NOTE: Don't want to do this before sedimentation, since sedimentation doesn't take into
68 : ! account the varying cloud fractions, and thus a particle scaled at one level and cloud
69 : ! fraction would be scaled inappropriately at another level and cloud fraction.
70 : !
71 : ! NOTE: All detrainment also happens only in the in-cloud portion.
72 1050624 : if (do_incloud) then
73 :
74 : ! First do the in-cloud processing.
75 :
76 : ! Convert "cloud" particles to in-cloud values.
77 : !
78 : ! NOTE: If a particle is a "cloud" particle, it means that the entire mass of the
79 : ! particle is in the incloud portion of the grid box. Particle that are not "cloud
80 : ! particles" have their mass spread throughout the grid box.
81 0 : pc_orig(:,:,:) = pc(:,:,:)
82 0 : gc_orig(:,:) = gc(:,:)
83 0 : t_orig(:) = t(:)
84 :
85 : ! If the cloud fraction gets too small it causes the microphysics to require a
86 : ! lot of substeps. Enforce a minimum cloud fraction for the purposes of scaling
87 : ! to incloud values.
88 0 : scale_cldfrc(:) = max(CLDFRC_MIN, cldfrc(:))
89 0 : scale_cldfrc(:) = min(1._f - CLDFRC_MIN, scale_cldfrc(:))
90 :
91 0 : do ielem = 1, NELEM
92 0 : igroup = igelem(ielem)
93 :
94 0 : if (is_grp_cloud(igroup)) then
95 0 : do ibin = 1, NBIN
96 0 : pc(:, ibin, ielem) = pc(:, ibin, ielem) / scale_cldfrc(:)
97 0 : pcd(:, ibin, ielem) = pcd(:, ibin, ielem) / scale_cldfrc(:)
98 : end do
99 : end if
100 : end do
101 :
102 0 : call newstate_calc(carma, cstate, scale_cldfrc(:), rc)
103 0 : if (rc < RC_OK) return
104 :
105 0 : call fixcorecol(carma, cstate, rc)
106 :
107 : ! Save the new in-cloud values for the gas, particle and temperature fields.
108 0 : pc_cloudy(:,:,:) = pc(:,:,:)
109 0 : gc_cloudy(:,:) = gc(:,:)
110 0 : t_cloudy(:) = t(:)
111 0 : rlheat_cloudy(:) = rlheat(:)
112 0 : partheat_cloudy(:) = partheat(:)
113 :
114 0 : if (do_substep) zsubsteps_cloudy(:) = zsubsteps(:)
115 :
116 : ! Now do the clear sky portion, using the original gridbox average concentrations.
117 : ! This is optional. If clear sky is not selected then all of the microphysics is
118 : ! done in-cloud.
119 0 : pc(:,:,:) = pc_orig(:,:,:)
120 0 : gc(:,:) = gc_orig(:,:)
121 0 : t(:) = t_orig(:)
122 :
123 0 : if (do_clearsky) then
124 :
125 : ! Convert "cloud" particles to clear sky values.
126 : !
127 : ! NOTE: If a particle is a "cloud" particle, it means that the entire mass of the
128 : ! particle is in the in-cloud portion of the grid box. They have no mass in the
129 : ! clear sky portion.
130 0 : do ielem = 1, NELEM
131 0 : igroup = igelem(ielem)
132 :
133 0 : if (is_grp_cloud(igroup)) then
134 0 : pc(:, :, ielem) = 0._f
135 0 : pcd(:, :, ielem) = 0._f
136 : end if
137 : end do
138 :
139 : ! Don't let the supersaturation be scaled by setting the cloud fraction used
140 : ! by the saturation code to 1.0. Any clouds formed in-situ in the clear sky
141 : ! are assumed to fill the grid box.
142 0 : cldfrc_orig(:) = cldfrc(:)
143 0 : cldfrc(:) = 1._f
144 :
145 : ! Recalculate supersaturation.
146 0 : do igas = 1, NGAS
147 0 : do iz = 1, NZ
148 0 : call supersat(carma, cstate, iz, igas, rc)
149 0 : if (rc < RC_OK) return
150 : end do
151 : end do
152 :
153 0 : call newstate_calc(carma, cstate, (1._f - scale_cldfrc(:)), rc)
154 0 : if (rc < RC_OK) return
155 :
156 0 : call fixcorecol(carma, cstate, rc)
157 :
158 : ! Restore the cloud fraction
159 0 : cldfrc(:) = cldfrc_orig(:)
160 :
161 : ! Save the new clear sky values for the gas, particle and temperature fields.
162 0 : pc_clear(:,:,:) = pc(:,:,:)
163 0 : gc_clear(:,:) = gc(:,:)
164 0 : t_clear(:) = t(:)
165 0 : rlheat_clear(:) = rlheat(:)
166 0 : partheat_clear(:) = partheat(:)
167 :
168 0 : if (do_substep) zsubsteps_clear(:) = zsubsteps(:)
169 :
170 : ! If not doing a clear sky calculation, then the clear sky portion reamins
171 : ! the same except for any contribution from advection.
172 : else
173 :
174 : ! NOTE: If a particle is a "cloud" particle, it means that the entire mass of the
175 : ! particle is in the in-cloud portion of the grid box. They have no mass in the
176 : ! clear sky portion.
177 0 : do ielem = 1, NELEM
178 0 : igroup = igelem(ielem)
179 :
180 0 : if (is_grp_cloud(igroup)) then
181 0 : pc_clear(:, :, ielem) = 0._f
182 : else
183 0 : pc_clear(:, :, ielem) = pc(:, :, ielem)
184 : end if
185 : end do
186 :
187 0 : do igas = 1, NGAS
188 0 : gc_clear(:,:) = gc(:,:)
189 : end do
190 :
191 0 : t_clear(:) = t(:)
192 0 : rlheat_clear(:) = 0._f
193 0 : partheat_clear(:) = 0._f
194 :
195 : ! If substepping, then add the advected part that is being doled out over
196 : ! the substeps.
197 0 : if (do_substep) then
198 0 : do igas = 1, NGAS
199 0 : gc_clear(:, igas) = gc_clear(:, igas) + d_gc(:, igas)
200 : end do
201 0 : t_clear(:) = t_clear(:) + d_t(:)
202 :
203 0 : zsubsteps_clear(:) = 0._f
204 : end if
205 : end if
206 :
207 :
208 : ! Add up the changes to the particle from the cloudy and clear sky components.
209 0 : do ielem = 1, NELEM
210 0 : igroup = igelem(ielem)
211 :
212 0 : do ibin = 1, NBIN
213 0 : pc(:, ibin, ielem) = (1._f - scale_cldfrc(:)) * pc_clear(:, ibin, ielem) + scale_cldfrc(:) * pc_cloudy(:, ibin, ielem)
214 : end do
215 : end do
216 :
217 0 : t(:) = (1._f - scale_cldfrc(:)) * t_clear(:) + scale_cldfrc(:) * t_cloudy(:)
218 :
219 0 : if (do_grow) then
220 0 : rlheat(:) = (1._f - scale_cldfrc(:)) * rlheat_clear(:) + scale_cldfrc(:) * rlheat_cloudy(:)
221 0 : partheat(:) = (1._f - scale_cldfrc(:)) * partheat_clear(:) + scale_cldfrc(:) * partheat_cloudy(:)
222 : end if
223 :
224 0 : do igas = 1, NGAS
225 0 : gc(:, igas) = (1._f - scale_cldfrc(:)) * gc_clear(:, igas) + scale_cldfrc(:) * gc_cloudy(:, igas)
226 :
227 : ! Recalculate gridbox average supersaturation.
228 0 : do iz = 1, NZ
229 0 : call supersat(carma, cstate, iz, igas, rc)
230 0 : if (rc < RC_OK) return
231 : end do
232 : end do
233 :
234 0 : if (do_substep) zsubsteps(:) = zsubsteps_clear(:) + zsubsteps_cloudy(:)
235 :
236 :
237 :
238 : ! No special in-cloud/clear sky processing, everything is gridbox average.
239 : else
240 34670592 : scale_threshold(:) = 1._f
241 1050624 : call newstate_calc(carma, cstate, scale_threshold, rc)
242 1050624 : if (rc < RC_OK) return
243 :
244 1050624 : call fixcorecol(carma, cstate, rc)
245 : end if
246 :
247 : ! Return to caller with new state computed
248 : return
249 1050624 : end
|
