Line data Source code
1 : module zm_conv_evap
2 :
3 : use ccpp_kinds, only: kind_phys
4 :
5 : ! CACNOTE - Need to ccpp'ize cloud_fraction
6 : use cloud_fraction, only: cldfrc_fice
7 :
8 : implicit none
9 :
10 : save
11 : private ! Make default type private to the module
12 : !
13 : ! PUBLIC: interfaces
14 : !
15 : public zm_conv_evap_run ! evaporation of precip from ZM schemea
16 :
17 : contains
18 :
19 :
20 : !===============================================================================
21 : !> \section arg_table_zm_conv_evap_run Argument Table
22 : !! \htmlinclude zm_conv_evap_run.html
23 : !!
24 2990736 : subroutine zm_conv_evap_run(ncol, pver, pverp, &
25 : gravit, latice, latvap, tmelt, &
26 : cpres, ke, ke_lnd, zm_org, &
27 5981472 : t,pmid,pdel,q, &
28 2990736 : landfrac, &
29 2990736 : tend_s, tend_s_snwprd, tend_s_snwevmlt, tend_q, &
30 2990736 : prdprec, cldfrc, deltat, &
31 11962944 : prec, snow, ntprprd, ntsnprd, flxprec, flxsnow, prdsnow)
32 :
33 :
34 : !-----------------------------------------------------------------------
35 : ! Compute tendencies due to evaporation of rain from ZM scheme
36 : !--
37 : ! Compute the total precipitation and snow fluxes at the surface.
38 : ! Add in the latent heat of fusion for snow formation and melt, since it not dealt with
39 : ! in the Zhang-MacFarlane parameterization.
40 : ! Evaporate some of the precip directly into the environment using a Sundqvist type algorithm
41 : !-----------------------------------------------------------------------
42 :
43 : !CACNOTE - Not sure what to do about qsat_water
44 : use wv_saturation, only: qsat
45 :
46 : !------------------------------Arguments--------------------------------
47 : integer,intent(in) :: ncol ! number of columns
48 : integer,intent(in) :: pver, pverp
49 : real(kind_phys),intent(in) :: gravit ! gravitational acceleration (m s-2)
50 : real(kind_phys),intent(in) :: latice ! Latent heat of fusion (J kg-1)
51 : real(kind_phys),intent(in) :: latvap ! Latent heat of vaporization (J kg-1)
52 : real(kind_phys),intent(in) :: tmelt ! Freezing point of water (K)
53 : real(kind_phys), intent(in) :: cpres ! specific heat at constant pressure in j/kg-degk.
54 : real(kind_phys), intent(in) :: ke ! Tunable evaporation efficiency set from namelist input zmconv_ke
55 : real(kind_phys), intent(in) :: ke_lnd
56 : logical, intent(in) :: zm_org
57 : real(kind_phys),intent(in), dimension(:,:) :: t ! temperature (K) (ncol,pver)
58 : real(kind_phys),intent(in), dimension(:,:) :: pmid ! midpoint pressure (Pa) (ncol,pver)
59 : real(kind_phys),intent(in), dimension(:,:) :: pdel ! layer thickness (Pa) (ncol,pver)
60 : real(kind_phys),intent(in), dimension(:,:) :: q ! water vapor (kg/kg) (ncol,pver)
61 : real(kind_phys),intent(in), dimension(:) :: landfrac ! land fraction (ncol)
62 : real(kind_phys),intent(inout), dimension(:,:) :: tend_s ! heating rate (J/kg/s) (ncol,pver)
63 : real(kind_phys),intent(inout), dimension(:,:) :: tend_q ! water vapor tendency (kg/kg/s) (ncol,pver)
64 : real(kind_phys),intent(out), dimension(:,:) :: tend_s_snwprd ! Heating rate of snow production (ncol,pver)
65 : real(kind_phys),intent(out), dimension(:,:) :: tend_s_snwevmlt ! Heating rate of evap/melting of snow (ncol,pver)
66 :
67 :
68 :
69 : real(kind_phys), intent(in ) :: prdprec(:,:)! precipitation production (kg/ks/s) (ncol,pver)
70 : real(kind_phys), intent(in ) :: cldfrc(:,:) ! cloud fraction (ncol,pver)
71 : real(kind_phys), intent(in ) :: deltat ! time step
72 :
73 : real(kind_phys), intent(inout) :: prec(:) ! Convective-scale preciptn rate (ncol)
74 : real(kind_phys), intent(out) :: snow(:) ! Convective-scale snowfall rate (ncol)
75 :
76 : real(kind_phys), optional, intent(in), allocatable :: prdsnow(:,:) ! snow production (kg/ks/s)
77 :
78 : !
79 : !---------------------------Local storage-------------------------------
80 :
81 5981472 : real(kind_phys) :: es (ncol,pver) ! Saturation vapor pressure
82 5981472 : real(kind_phys) :: fice (ncol,pver) ! ice fraction in precip production
83 5981472 : real(kind_phys) :: fsnow_conv(ncol,pver) ! snow fraction in precip production
84 5981472 : real(kind_phys) :: qs (ncol,pver) ! saturation specific humidity
85 : real(kind_phys),intent(out) :: flxprec(:,:) ! Convective-scale flux of precip at interfaces (kg/m2/s) ! (ncol,pverp)
86 : real(kind_phys),intent(out) :: flxsnow(:,:) ! Convective-scale flux of snow at interfaces (kg/m2/s) ! (ncol,pverp)
87 : real(kind_phys),intent(out) :: ntprprd(:,:) ! net precip production in layer ! (ncol,pver)
88 : real(kind_phys),intent(out) :: ntsnprd(:,:) ! net snow production in layer ! (ncol,pver)
89 : real(kind_phys) :: work1 ! temp variable (pjr)
90 : real(kind_phys) :: work2 ! temp variable (pjr)
91 :
92 5981472 : real(kind_phys) :: evpvint(ncol) ! vertical integral of evaporation
93 5981472 : real(kind_phys) :: evpprec(ncol) ! evaporation of precipitation (kg/kg/s)
94 5981472 : real(kind_phys) :: evpsnow(ncol) ! evaporation of snowfall (kg/kg/s)
95 5981472 : real(kind_phys) :: snowmlt(ncol) ! snow melt tendency in layer
96 5981472 : real(kind_phys) :: flxsntm(ncol) ! flux of snow into layer, after melting
97 :
98 : real(kind_phys) :: kemask
99 : real(kind_phys) :: evplimit ! temp variable for evaporation limits
100 : real(kind_phys) :: rlat(ncol)
101 : real(kind_phys) :: dum
102 : real(kind_phys) :: omsm
103 :
104 : integer :: i,k ! longitude,level indices
105 : logical :: old_snow
106 :
107 :
108 : !-----------------------------------------------------------------------
109 :
110 : ! If prdsnow is passed in and allocated, then use it in the calculation, otherwise
111 : ! use the old snow calculation
112 2990736 : old_snow=.true.
113 2990736 : if (present(prdsnow)) then
114 0 : if (allocated(prdsnow)) then
115 0 : old_snow=.false.
116 : end if
117 : end if
118 :
119 : ! convert input precip to kg/m2/s
120 49938336 : prec(:ncol) = prec(:ncol)*1000._kind_phys
121 :
122 : ! determine saturation vapor pressure
123 80749872 : do k = 1,pver
124 80749872 : call qsat(t(1:ncol,k), pmid(1:ncol,k), es(1:ncol,k), qs(1:ncol,k), ncol)
125 : end do
126 : ! determine ice fraction in rain production (use cloud water parameterization fraction at present)
127 : !REMOVECAM - no longer need these when CAM is retired and pcols no longer exists
128 1301387472 : fice(:,:) = 0._kind_phys
129 1301387472 : fsnow_conv(:,:) = 0._kind_phys
130 : !REMOVECAM_END
131 2990736 : call cldfrc_fice(ncol, t(1:ncol,:), fice(1:ncol,:), fsnow_conv(1:ncol,:))
132 :
133 : ! zero the flux integrals on the top boundary
134 49938336 : flxprec(:ncol,1) = 0._kind_phys
135 49938336 : flxsnow(:ncol,1) = 0._kind_phys
136 49938336 : evpvint(:ncol) = 0._kind_phys
137 : omsm=0.9999_kind_phys
138 :
139 80749872 : do k = 1, pver
140 1301387472 : do i = 1, ncol
141 :
142 : ! Melt snow falling into layer, if necessary.
143 1220637600 : if( old_snow ) then
144 1220637600 : if (t(i,k) > tmelt) then
145 197813818 : flxsntm(i) = 0._kind_phys
146 197813818 : snowmlt(i) = flxsnow(i,k) * gravit/ pdel(i,k)
147 : else
148 1022823782 : flxsntm(i) = flxsnow(i,k)
149 1022823782 : snowmlt(i) = 0._kind_phys
150 : end if
151 : else
152 : ! make sure melting snow doesn't reduce temperature below threshold
153 0 : if (t(i,k) > tmelt) then
154 0 : dum = -latice/cpres*flxsnow(i,k)*gravit/pdel(i,k)*deltat
155 0 : if (t(i,k) + dum .le. tmelt) then
156 0 : dum = (t(i,k)-tmelt)*cpres/latice/deltat
157 0 : dum = dum/(flxsnow(i,k)*gravit/pdel(i,k))
158 0 : dum = max(0._kind_phys,dum)
159 0 : dum = min(1._kind_phys,dum)
160 : else
161 : dum = 1._kind_phys
162 : end if
163 0 : dum = dum*omsm
164 0 : flxsntm(i) = flxsnow(i,k)*(1.0_kind_phys-dum)
165 0 : snowmlt(i) = dum*flxsnow(i,k)*gravit/ pdel(i,k)
166 : else
167 0 : flxsntm(i) = flxsnow(i,k)
168 0 : snowmlt(i) = 0._kind_phys
169 : end if
170 : end if
171 :
172 : ! relative humidity depression must be > 0 for evaporation
173 1220637600 : evplimit = max(1._kind_phys - q(i,k)/qs(i,k), 0._kind_phys)
174 :
175 1220637600 : if (zm_org) then
176 0 : kemask = ke * (1._kind_phys - landfrac(i)) + ke_lnd * landfrac(i)
177 : else
178 1220637600 : kemask = ke
179 : endif
180 :
181 : ! total evaporation depends on flux in the top of the layer
182 : ! flux prec is the net production above layer minus evaporation into environmet
183 1220637600 : evpprec(i) = kemask * (1._kind_phys - cldfrc(i,k)) * evplimit * sqrt(flxprec(i,k))
184 :
185 : ! Don't let evaporation supersaturate layer (approx). Layer may already be saturated.
186 : ! Currently does not include heating/cooling change to qs
187 1220637600 : evplimit = max(0._kind_phys, (qs(i,k)-q(i,k)) / deltat)
188 :
189 : ! Don't evaporate more than is falling into the layer - do not evaporate rain formed
190 : ! in this layer but if precip production is negative, remove from the available precip
191 : ! Negative precip production occurs because of evaporation in downdrafts.
192 1220637600 : evplimit = min(evplimit, flxprec(i,k) * gravit / pdel(i,k))
193 :
194 : ! Total evaporation cannot exceed input precipitation
195 1220637600 : evplimit = min(evplimit, (prec(i) - evpvint(i)) * gravit / pdel(i,k))
196 :
197 1220637600 : evpprec(i) = min(evplimit, evpprec(i))
198 1220637600 : if( .not.old_snow ) then
199 0 : evpprec(i) = max(0._kind_phys, evpprec(i))
200 0 : evpprec(i) = evpprec(i)*omsm
201 : end if
202 :
203 :
204 : ! evaporation of snow depends on snow fraction of total precipitation in the top after melting
205 1220637600 : if (flxprec(i,k) > 0._kind_phys) then
206 : ! prevent roundoff problems
207 116856855 : work1 = min(max(0._kind_phys,flxsntm(i)/flxprec(i,k)),1._kind_phys)
208 116856855 : evpsnow(i) = evpprec(i) * work1
209 : else
210 1103780745 : evpsnow(i) = 0._kind_phys
211 : end if
212 :
213 : ! vertically integrated evaporation
214 1220637600 : evpvint(i) = evpvint(i) + evpprec(i) * pdel(i,k)/gravit
215 :
216 : ! net precip production is production - evaporation
217 1220637600 : ntprprd(i,k) = prdprec(i,k) - evpprec(i)
218 : ! net snow production is precip production * ice fraction - evaporation - melting
219 : ! the small amount added to flxprec in the work1 expression has been increased from
220 : ! 1e-36 to 8.64e-11 (1e-5 mm/day). This causes the temperature based partitioning
221 : ! scheme to be used for small flxprec amounts. This is to address error growth problems.
222 :
223 1220637600 : if( old_snow ) then
224 1220637600 : if (flxprec(i,k).gt.0._kind_phys) then
225 116856855 : work1 = min(max(0._kind_phys,flxsnow(i,k)/flxprec(i,k)),1._kind_phys)
226 : else
227 : work1 = 0._kind_phys
228 : endif
229 :
230 1220637600 : work2 = max(fsnow_conv(i,k), work1)
231 1220637600 : if (snowmlt(i).gt.0._kind_phys) work2 = 0._kind_phys
232 1220637600 : ntsnprd(i,k) = prdprec(i,k)*work2 - evpsnow(i) - snowmlt(i)
233 1220637600 : tend_s_snwprd (i,k) = prdprec(i,k)*work2*latice
234 1220637600 : tend_s_snwevmlt(i,k) = - ( evpsnow(i) + snowmlt(i) )*latice
235 : else
236 0 : ntsnprd(i,k) = prdsnow(i,k) - min(flxsnow(i,k)*gravit/pdel(i,k), evpsnow(i)+snowmlt(i))
237 0 : tend_s_snwprd (i,k) = prdsnow(i,k)*latice
238 0 : tend_s_snwevmlt(i,k) = -min(flxsnow(i,k)*gravit/pdel(i,k), evpsnow(i)+snowmlt(i) )*latice
239 : end if
240 :
241 : ! precipitation fluxes
242 1220637600 : flxprec(i,k+1) = flxprec(i,k) + ntprprd(i,k) * pdel(i,k)/gravit
243 1220637600 : flxsnow(i,k+1) = flxsnow(i,k) + ntsnprd(i,k) * pdel(i,k)/gravit
244 :
245 : ! protect against rounding error
246 1220637600 : flxprec(i,k+1) = max(flxprec(i,k+1), 0._kind_phys)
247 1220637600 : flxsnow(i,k+1) = max(flxsnow(i,k+1), 0._kind_phys)
248 :
249 : ! heating (cooling) and moistening due to evaporation
250 : ! - latent heat of vaporization for precip production has already been accounted for
251 : ! - snow is contained in prec
252 1220637600 : if( old_snow ) then
253 1220637600 : tend_s(i,k) =-evpprec(i)*latvap + ntsnprd(i,k)*latice
254 : else
255 0 : tend_s(i,k) =-evpprec(i)*latvap + tend_s_snwevmlt(i,k)
256 : end if
257 1298396736 : tend_q(i,k) = evpprec(i)
258 : end do
259 : end do
260 :
261 : ! set output precipitation rates (m/s)
262 49938336 : prec(:ncol) = flxprec(:ncol,pver+1) / 1000._kind_phys
263 49938336 : snow(:ncol) = flxsnow(:ncol,pver+1) / 1000._kind_phys
264 :
265 2990736 : end subroutine zm_conv_evap_run
266 :
267 :
268 : end module zm_conv_evap
|
