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