| blob | caf0ba70a4e6e105b58d5f4f026d94edbf61eefa |
1 MODULE module_sf_lake
3 ! The lake scheme was retrieved from the Community Land Model version 4.5
4 ! (Oleson et al. 2013) with some modifications by Gu et al. (2013). It is a
5 ! one-dimensional mass and energy balance scheme with 20-25 model layers,
6 ! including up to 5 snow layers on the lake ice, 10 water layers, and 10 soil
7 ! layers on the lake bottom. The lake scheme is used with actual lake points and
8 ! lake depth derived from the WPS, and it also can be used with user defined
9 ! lake points and lake depth in WRF (lake_min_elev and lakedepth_default).
10 ! The lake scheme is independent of a land surface scheme and therefore
11 ! can be used with any land surface scheme embedded in WRF. The lake scheme
12 ! developments and evaluations were included in Subin et al. (2012) and Gu et al. (2013)
13 !
14 ! Subin et al. 2012: Improved lake model for climate simulations, J. Adv. Model.
15 ! Earth Syst., 4, M02001. DOI:10.1029/2011MS000072;
16 ! Gu et al. 2013: Calibration and validation of lake surface temperature simulations
17 ! with the coupled WRF-Lake model. Climatic Change, 1-13, 10.1007/s10584-013-0978-y.
19 USE module_wrf_error
20 USE module_model_constants, ONLY : rcp
22 implicit none
23 integer, parameter :: r8 = selected_real_kind(12)
25 integer, parameter :: nlevsoil = 10 ! number of soil layers
26 integer, parameter :: nlevlake = 10 ! number of lake layers
27 integer, parameter :: nlevsnow = 5 ! maximum number of snow layers
29 integer,parameter :: lbp = 1 ! pft-index bounds
30 integer,parameter :: ubp = 1
31 integer,parameter :: lbc = 1 ! column-index bounds
32 integer,parameter :: ubc = 1
33 integer,parameter :: num_shlakec = 1 ! number of columns in lake filter
34 integer,parameter :: filter_shlakec(1) = 1 ! lake filter (columns)
35 integer,parameter :: num_shlakep = 1 ! number of pfts in lake filter
36 integer,parameter :: filter_shlakep(1) = 1 ! lake filter (pfts)
37 integer,parameter :: pcolumn(1) = 1
38 integer,parameter :: pgridcell(1) = 1
39 integer,parameter :: cgridcell(1) = 1 ! gridcell index of column
40 integer,parameter :: clandunit(1) = 1 ! landunit index of column
42 integer,parameter :: begg = 1
43 integer,parameter :: endg = 1
44 integer,parameter :: begl = 1
45 integer,parameter :: endl = 1
46 integer,parameter :: begc = 1
47 integer,parameter :: endc = 1
48 integer,parameter :: begp = 1
49 integer,parameter :: endp = 1
51 integer,parameter :: column =1
52 logical,parameter :: lakpoi(1) = .true.
57 !Initialize physical constants:
58 real(r8), parameter :: vkc = 0.4_r8 !von Karman constant [-]
59 real(r8), parameter :: pie = 3.141592653589793_r8 ! pi
60 real(r8), parameter :: grav = 9.80616_r8 !gravity constant [m/s2]
61 real(r8), parameter :: sb = 5.67e-8_r8 !stefan-boltzmann constant [W/m2/K4]
62 real(r8), parameter :: tfrz = 273.16_r8 !freezing temperature [K]
63 real(r8), parameter :: denh2o = 1.000e3_r8 !density of liquid water [kg/m3]
64 real(r8), parameter :: denice = 0.917e3_r8 !density of ice [kg/m3]
65 real(r8), parameter :: cpice = 2.11727e3_r8 !Specific heat of ice [J/kg-K]
66 real(r8), parameter :: cpliq = 4.188e3_r8 !Specific heat of water [J/kg-K]
67 real(r8), parameter :: hfus = 3.337e5_r8 !Latent heat of fusion for ice [J/kg]
68 real(r8), parameter :: hvap = 2.501e6_r8 !Latent heat of evap for water [J/kg]
69 real(r8), parameter :: hsub = 2.501e6_r8+3.337e5_r8 !Latent heat of sublimation [J/kg]
70 real(r8), parameter :: rair = 287.0423_r8 !gas constant for dry air [J/kg/K]
71 real(r8), parameter :: cpair = 1.00464e3_r8 !specific heat of dry air [J/kg/K]
72 real(r8), parameter :: tcrit = 2.5 !critical temperature to determine rain or snow
73 real(r8), parameter :: tkwat = 0.6 !thermal conductivity of water [W/m/k]
74 real(r8), parameter :: tkice = 2.290 !thermal conductivity of ice [W/m/k]
75 real(r8), parameter :: tkairc = 0.023 !thermal conductivity of air [W/m/k]
76 real(r8), parameter :: bdsno = 250. !bulk density snow (kg/m**3)
78 real(r8), public, parameter :: spval = 1.e36 !special value for missing data (ocean)
80 real, parameter :: depth_c = 50. ! below the level t_lake3d will be 277.0 !mchen
83 ! These are tunable constants
84 real(r8), parameter :: wimp = 0.05 !Water impremeable if porosity less than wimp
85 real(r8), parameter :: ssi = 0.033 !Irreducible water saturation of snow
86 real(r8), parameter :: cnfac = 0.5 !Crank Nicholson factor between 0 and 1
89 ! Initialize water type constants
90 integer,parameter :: istsoil = 1 !soil "water" type
91 integer, private :: i ! loop index
92 real(r8) :: dtime ! land model time step (sec)
94 real(r8) :: zlak(1:nlevlake) !lake z (layers)
95 real(r8) :: dzlak(1:nlevlake) !lake dz (thickness)
96 real(r8) :: zsoi(1:nlevsoil) !soil z (layers)
97 real(r8) :: dzsoi(1:nlevsoil) !soil dz (thickness)
98 real(r8) :: zisoi(0:nlevsoil) !soil zi (interfaces)
101 real(r8) :: sand(19) ! percent sand
102 real(r8) :: clay(19) ! percent clay
104 data(sand(i), i=1,19)/92.,80.,66.,20.,5.,43.,60.,&
105 10.,32.,51., 6.,22.,39.7,0.,100.,54.,17.,100.,92./
107 data(clay(i), i=1,19)/ 3., 5.,10.,15.,5.,18.,27.,&
108 33.,33.,41.,47.,58.,14.7,0., 0., 8.5,54., 0., 3./
111 ! real(r8) :: dtime ! land model time step (sec)
112 real(r8) :: watsat(1,nlevsoil) ! volumetric soil water at saturation (porosity)
113 real(r8) :: tksatu(1,nlevsoil) ! thermal conductivity, saturated soil [W/m-K]
114 real(r8) :: tkmg(1,nlevsoil) ! thermal conductivity, soil minerals [W/m-K]
115 real(r8) :: tkdry(1,nlevsoil) ! thermal conductivity, dry soil (W/m/Kelvin)
116 real(r8) :: csol(1,nlevsoil) ! heat capacity, soil solids (J/m**3/Kelvin)
117 CONTAINS
120 SUBROUTINE Lake( t_phy ,p8w ,dz8w ,qvcurr ,& !i
121 u_phy ,v_phy , glw ,emiss ,&
122 rainbl ,dtbl ,swdown ,albedo ,&
123 xlat_urb2d ,z_lake3d ,dz_lake3d ,lakedepth2d ,&
124 watsat3d ,csol3d ,tkmg3d ,tkdry3d ,&
125 tksatu3d ,ivgtyp ,ht ,xland ,&
126 iswater, xice, xice_threshold, lake_min_elev ,&
127 ids ,ide ,jds ,jde ,&
128 kds ,kde ,ims ,ime ,&
129 jms ,jme ,kms ,kme ,&
130 its ,ite ,jts ,jte ,&
131 kts ,kte ,&
132 h2osno2d ,snowdp2d ,snl2d ,z3d ,& !h
133 dz3d ,zi3d ,h2osoi_vol3d ,h2osoi_liq3d ,&
134 h2osoi_ice3d ,t_grnd2d ,t_soisno3d ,t_lake3d ,&
135 savedtke12d ,lake_icefrac3d ,&
136 #if (EM_CORE==1)
137 ! lakemask ,lakeflag ,&
138 lakemask ,&
139 #endif
140 hfx ,lh ,grdflx ,tsk ,& !o
141 qfx ,t2 ,th2 ,q2 )
143 !==============================================================================
144 ! This subroutine was first edited by Hongping Gu and Jiming Jin for coupling
145 ! 07/20/2010
146 !==============================================================================
147 IMPLICIT NONE
149 !in:
151 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
152 ims,ime, jms,jme, kms,kme, &
153 its,ite, jts,jte, kts,kte
154 INTEGER , INTENT (IN) :: iswater
155 REAL, INTENT(IN) :: xice_threshold
156 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: XICE
157 #if (EM_CORE==1)
158 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: LAKEMASK
159 ! INTEGER, INTENT(IN):: LAKEFLAG
160 #endif
162 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: t_phy
163 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: p8w
164 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: dz8w
165 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: qvcurr
166 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: U_PHY
167 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: V_PHY
168 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: glw
169 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: emiss
170 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: rainbl
171 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: swdown
172 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(INOUT) :: albedo
173 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: XLAND
174 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: XLAT_URB2D
175 INTEGER, DIMENSION( ims:ime, jms:jme ) ,INTENT(INOUT) :: IVGTYP
176 REAL, INTENT(IN) :: dtbl
178 REAL, DIMENSION( ims:ime,1:nlevlake,jms:jme ),INTENT(IN) :: z_lake3d
179 REAL, DIMENSION( ims:ime,1:nlevlake,jms:jme ),INTENT(IN) :: dz_lake3d
180 REAL, DIMENSION( ims:ime,1:nlevsoil,jms:jme ),INTENT(IN) :: watsat3d
181 REAL, DIMENSION( ims:ime,1:nlevsoil,jms:jme ),INTENT(IN) :: csol3d
182 REAL, DIMENSION( ims:ime,1:nlevsoil,jms:jme ),INTENT(IN) :: tkmg3d
183 REAL, DIMENSION( ims:ime,1:nlevsoil,jms:jme ),INTENT(IN) :: tkdry3d
184 REAL, DIMENSION( ims:ime,1:nlevsoil,jms:jme ),INTENT(IN) :: tksatu3d
185 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: lakedepth2d
186 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(IN) :: ht
187 REAL ,INTENT(IN) :: lake_min_elev
189 !out:
190 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: HFX
191 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: LH
192 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: GRDFLX
193 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: TSK
194 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: QFX
195 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: T2
196 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: TH2
197 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(OUT) :: Q2
199 !in&out:
201 real, dimension(ims:ime,jms:jme ) ,intent(inout) :: savedtke12d
202 real, dimension(ims:ime,jms:jme ) ,intent(inout) :: snowdp2d, &
203 h2osno2d, &
204 snl2d, &
205 t_grnd2d
207 real, dimension( ims:ime,1:nlevlake, jms:jme ) ,INTENT(inout) :: t_lake3d, &
208 lake_icefrac3d
209 real, dimension( ims:ime,-nlevsnow+1:nlevsoil, jms:jme ) ,INTENT(inout) :: t_soisno3d, &
210 h2osoi_ice3d, &
211 h2osoi_liq3d, &
212 h2osoi_vol3d, &
213 z3d, &
214 dz3d
215 real, dimension( ims:ime,-nlevsnow+0:nlevsoil, jms:jme ) ,INTENT(inout) :: zi3d
218 !local variable:
220 REAL :: SFCTMP,PBOT,PSFC,ZLVL,Q2K,EMISSI,LWDN,PRCP,SOLDN,SOLNET
221 INTEGER :: C,i,j,k
224 !tempory varibles in:
225 real(r8) :: forc_t(1) ! atmospheric temperature (Kelvin)
226 real(r8) :: forc_pbot(1) ! atm bottom level pressure (Pa)
227 real(r8) :: forc_psrf(1) ! atmospheric surface pressure (Pa)
228 real(r8) :: forc_hgt(1) ! atmospheric reference height (m)
229 real(r8) :: forc_hgt_q(1) ! observational height of humidity [m]
230 real(r8) :: forc_hgt_t(1) ! observational height of temperature [m]
231 real(r8) :: forc_hgt_u(1) ! observational height of wind [m]
232 real(r8) :: forc_q(1) ! atmospheric specific humidity (kg/kg)
233 real(r8) :: forc_u(1) ! atmospheric wind speed in east direction (m/s)
234 real(r8) :: forc_v(1) ! atmospheric wind speed in north direction (m/s)
235 ! real(r8) :: forc_rho(1) ! density (kg/m**3)
236 real(r8) :: forc_lwrad(1) ! downward infrared (longwave) radiation (W/m**2)
237 real(r8) :: prec(1) ! snow or rain rate [mm/s]
238 real(r8) :: sabg(1) ! solar radiation absorbed by ground (W/m**2)
239 real(r8) :: lat(1) ! latitude (radians)
240 real(r8) :: z_lake(1,nlevlake) ! layer depth for lake (m)
241 real(r8) :: dz_lake(1,nlevlake) ! layer thickness for lake (m)
243 real(r8) :: lakedepth(1) ! column lake depth (m)
244 logical :: do_capsnow(1) ! true => do snow capping
246 !in&out
247 real(r8) :: h2osoi_vol(1,-nlevsnow+1:nlevsoil) ! volumetric soil water (0<=h2osoi_vol<=watsat)[m3/m3]
248 real(r8) :: t_grnd(1) ! ground temperature (Kelvin)
249 real(r8) :: h2osno(1) ! snow water (mm H2O)
250 real(r8) :: snowdp(1) ! snow height (m)
251 real(r8) :: z(1,-nlevsnow+1:nlevsoil) ! layer depth for snow & soil (m)
252 real(r8) :: dz(1,-nlevsnow+1:nlevsoil) ! layer thickness for soil or snow (m)
253 real(r8) :: t_soisno(1,-nlevsnow+1:nlevsoil) ! soil (or snow) temperature (Kelvin)
254 real(r8) :: t_lake(1,nlevlake) ! lake temperature (Kelvin)
255 integer :: snl(1) ! number of snow layers
256 real(r8) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) ! liquid water (kg/m2)
257 real(r8) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) ! ice lens (kg/m2)
258 real(r8) :: savedtke1(1) ! top level eddy conductivity from previous timestep (W/m.K)
259 real(r8) :: zi(1,-nlevsnow+0:nlevsoil) ! interface level below a "z" level (m)
260 real(r8) :: lake_icefrac(1,nlevlake) ! mass fraction of lake layer that is frozen
263 !out:
264 real(r8) :: eflx_gnet(1) !net heat flux into ground (W/m**2)
265 real(r8) :: eflx_lwrad_net(1) ! net infrared (longwave) rad (W/m**2) [+ = to atm]
266 real(r8) :: eflx_sh_tot(1) ! total sensible heat flux (W/m**2) [+ to atm]
267 real(r8) :: eflx_lh_tot(1) ! total latent heat flux (W/m8*2) [+ to atm]
268 real(r8) :: t_ref2m(1) ! 2 m height surface air temperature (Kelvin)
269 real(r8) :: q_ref2m(1) ! 2 m height surface specific humidity (kg/kg)
270 real(r8) :: taux(1) ! wind (shear) stress: e-w (kg/m/s**2)
271 real(r8) :: tauy(1) ! wind (shear) stress: n-s (kg/m/s**2)
272 real(r8) :: ram1(1) ! aerodynamical resistance (s/m)
273 ! for calculation of decay of eddy diffusivity with depth
274 ! Change the type variable to pass back to WRF.
275 real(r8) :: z0mg(1) ! roughness length over ground, momentum (m(
278 dtime = dtbl
280 DO J = jts,jte
281 DO I = its,ite
283 SFCTMP = t_phy(i,1,j)
284 PBOT = p8w(i,2,j)
285 PSFC = P8w(i,1,j)
286 ZLVL = 0.5 * dz8w(i,1,j)
287 Q2K = qvcurr(i,1,j)/(1.0 + qvcurr(i,1,j))
288 EMISSI = EMISS(I,J)
289 LWDN = GLW(I,J)*EMISSI
290 PRCP = RAINBL(i,j)/dtbl
291 SOLDN = SWDOWN(I,J) ! SOLDN is total incoming solar
292 SOLNET = SOLDN*(1.-ALBEDO(I,J)) ! use mid-day albedo to determine net downward solar
293 ! (no solar zenith angle correction)
294 ! IF (XLAND(I,J).GT.1.5) THEN
296 ! if ( xice(i,j).gt.xice_threshold) then
297 ! ivgtyp(i,j) = iswater
298 ! xland(i,j) = 2.
299 ! lake_icefrac3d(i,1,j) = xice(i,j)
300 ! endif
302 #if (EM_CORE==1)
303 if (lakemask(i,j).eq.1) THEN
304 #else
305 if (ivgtyp(i,j)==iswater.and.ht(i,j)>= lake_min_elev ) THEN
306 #endif
308 do c = 1,column
310 forc_t(c) = SFCTMP ! [K]
311 forc_pbot(c) = PBOT
312 forc_psrf(c) = PSFC
313 forc_hgt(c) = ZLVL ! [m]
314 forc_hgt_q(c) = ZLVL ! [m]
315 forc_hgt_t(c) = ZLVL ! [m]
316 forc_hgt_u(c) = ZLVL ! [m]
317 forc_q(c) = Q2K ! [kg/kg]
318 forc_u(c) = U_PHY(I,1,J)
319 forc_v(c) = V_PHY(I,1,J)
320 ! forc_rho(c) = SFCPRS / (287.04 * SFCTMP * (1.0+ 0.61 * Q2K)) ![kg/m/m/m]
321 forc_lwrad(c) = LWDN ! [W/m/m]
322 prec(c) = PRCP ! [mm/s]
323 sabg(c) = SOLNET
324 lat(c) = XLAT_URB2D(I,J)*pie/180 ! [radian]
325 do_capsnow(c) = .false.
327 lakedepth(c) = lakedepth2d(i,j)
328 savedtke1(c) = savedtke12d(i,j)
329 snowdp(c) = snowdp2d(i,j)
330 h2osno(c) = h2osno2d(i,j)
331 snl(c) = snl2d(i,j)
332 t_grnd(c) = t_grnd2d(i,j)
333 do k = 1,nlevlake
334 t_lake(c,k) = t_lake3d(i,k,j)
335 lake_icefrac(c,k) = lake_icefrac3d(i,k,j)
336 z_lake(c,k) = z_lake3d(i,k,j)
337 dz_lake(c,k) = dz_lake3d(i,k,j)
338 enddo
339 do k = -nlevsnow+1,nlevsoil
340 t_soisno(c,k) = t_soisno3d(i,k,j)
341 h2osoi_ice(c,k) = h2osoi_ice3d(i,k,j)
342 h2osoi_liq(c,k) = h2osoi_liq3d(i,k,j)
343 h2osoi_vol(c,k) = h2osoi_vol3d(i,k,j)
344 z(c,k) = z3d(i,k,j)
345 dz(c,k) = dz3d(i,k,j)
346 enddo
347 do k = -nlevsnow+0,nlevsoil
348 zi(c,k) = zi3d(i,k,j)
349 enddo
350 do k = 1,nlevsoil
351 watsat(c,k) = watsat3d(i,k,j)
352 csol(c,k) = csol3d(i,k,j)
353 tkmg(c,k) = tkmg3d(i,k,j)
354 tkdry(c,k) = tkdry3d(i,k,j)
355 tksatu(c,k) = tksatu3d(i,k,j)
356 enddo
358 enddo
359 CALL LakeMain(forc_t,forc_pbot,forc_psrf,forc_hgt,forc_hgt_q, & !I
360 forc_hgt_t,forc_hgt_u,forc_q, forc_u, &
361 forc_v,forc_lwrad,prec, sabg,lat, &
362 z_lake,dz_lake,lakedepth,do_capsnow, &
363 h2osno,snowdp,snl,z,dz,zi, & !H
364 h2osoi_vol,h2osoi_liq,h2osoi_ice, &
365 t_grnd,t_soisno,t_lake, &
366 savedtke1,lake_icefrac, &
367 eflx_lwrad_net,eflx_gnet, & !O
368 eflx_sh_tot,eflx_lh_tot, &
369 t_ref2m,q_ref2m, &
370 taux,tauy,ram1,z0mg)
373 do c = 1,column
374 HFX(I,J) = eflx_sh_tot(c) ![W/m/m]
375 LH(I,J) = eflx_lh_tot(c) !W/m/m]
376 GRDFLX(I,J) = eflx_gnet(c) !W/m/m]
377 TSK(I,J) = t_grnd(c) ![K]
378 T2(I,J) = t_ref2m(c)
379 TH2(I,J) = T2(I,J)*(1.E5/PSFC)**RCP
380 Q2(I,J) = q_ref2m(c)
381 albedo(i,j) = ( 0.6 * lake_icefrac(c,1) ) + ( (1.0-lake_icefrac(c,1)) * 0.08)
383 if( tsk(i,j) >= tfrz ) then
384 qfx(i,j) = eflx_lh_tot(c)/hvap
385 else
386 qfx(i,j) = eflx_lh_tot(c)/hsub ! heat flux (W/m^2)=>mass flux(kg/(sm^2))
387 endif
388 enddo
390 ! Renew Lake State Varialbes:(14)
391 do c = 1,column
393 savedtke12d(i,j) = savedtke1(c)
394 snowdp2d(i,j) = snowdp(c)
395 h2osno2d(i,j) = h2osno(c)
396 snl2d(i,j) = snl(c)
397 t_grnd2d(i,j) = t_grnd(c)
398 do k = 1,nlevlake
399 t_lake3d(i,k,j) = t_lake(c,k)
400 lake_icefrac3d(i,k,j) = lake_icefrac(c,k)
401 enddo
402 do k = -nlevsnow+1,nlevsoil
403 z3d(i,k,j) = z(c,k)
404 dz3d(i,k,j) = dz(c,k)
405 t_soisno3d(i,k,j) = t_soisno(c,k)
406 h2osoi_liq3d(i,k,j) = h2osoi_liq(c,k)
407 h2osoi_ice3d(i,k,j) = h2osoi_ice(c,k)
408 h2osoi_vol3d(i,k,j) = h2osoi_vol(c,k)
409 enddo
410 do k = -nlevsnow+0,nlevsoil
411 zi3d(i,k,j) = zi(c,k)
412 enddo
414 enddo
416 endif
417 ! ENDIF ! if xland = 2
418 ENDDO
419 ENDDO
421 END SUBROUTINE Lake
424 SUBROUTINE LakeMain(forc_t,forc_pbot,forc_psrf,forc_hgt,forc_hgt_q, & !I
425 forc_hgt_t,forc_hgt_u,forc_q, forc_u, &
426 forc_v,forc_lwrad,prec, sabg,lat, &
427 z_lake,dz_lake,lakedepth,do_capsnow, &
428 h2osno,snowdp,snl,z,dz,zi, & !H
429 h2osoi_vol,h2osoi_liq,h2osoi_ice, &
430 t_grnd,t_soisno,t_lake, &
431 savedtke1,lake_icefrac, &
432 eflx_lwrad_net,eflx_gnet, & !O
433 eflx_sh_tot,eflx_lh_tot, &
434 t_ref2m,q_ref2m, &
435 taux,tauy,ram1,z0mg)
436 implicit none
437 !in:
439 real(r8),intent(in) :: forc_t(1) ! atmospheric temperature (Kelvin)
440 real(r8),intent(in) :: forc_pbot(1) ! atm bottom level pressure (Pa)
441 real(r8),intent(in) :: forc_psrf(1) ! atmospheric surface pressure (Pa)
442 real(r8),intent(in) :: forc_hgt(1) ! atmospheric reference height (m)
443 real(r8),intent(in) :: forc_hgt_q(1) ! observational height of humidity [m]
444 real(r8),intent(in) :: forc_hgt_t(1) ! observational height of temperature [m]
445 real(r8),intent(in) :: forc_hgt_u(1) ! observational height of wind [m]
446 real(r8),intent(in) :: forc_q(1) ! atmospheric specific humidity (kg/kg)
447 real(r8),intent(in) :: forc_u(1) ! atmospheric wind speed in east direction (m/s)
448 real(r8),intent(in) :: forc_v(1) ! atmospheric wind speed in north direction (m/s)
449 ! real(r8),intent(in) :: forc_rho(1) ! density (kg/m**3)
450 real(r8),intent(in) :: forc_lwrad(1) ! downward infrared (longwave) radiation (W/m**2)
451 real(r8),intent(in) :: prec(1) ! snow or rain rate [mm/s]
452 real(r8),intent(in) :: sabg(1) ! solar radiation absorbed by ground (W/m**2)
453 real(r8),intent(in) :: lat(1) ! latitude (radians)
454 real(r8),intent(in) :: z_lake(1,nlevlake) ! layer depth for lake (m)
455 real(r8),intent(in) :: dz_lake(1,nlevlake) ! layer thickness for lake (m)
457 real(r8), intent(in) :: lakedepth(1) ! column lake depth (m)
458 !!!!!!!!!!!!!!!!tep(in),hydro(in)
459 ! real(r8), intent(in) :: watsat(1,1:nlevsoil) ! volumetric soil water at saturation (porosity)
460 !!!!!!!!!!!!!!!!hydro
461 logical , intent(in) :: do_capsnow(1) ! true => do snow capping
465 !in&out
466 real(r8),intent(inout) :: h2osoi_vol(1,-nlevsnow+1:nlevsoil) ! volumetric soil water (0<=h2osoi_vol<=watsat)[m3/m3]
467 real(r8),intent(inout) :: t_grnd(1) ! ground temperature (Kelvin)
468 real(r8),intent(inout) :: h2osno(1) ! snow water (mm H2O)
469 real(r8),intent(inout) :: snowdp(1) ! snow height (m)
470 real(r8),intent(inout) :: z(1,-nlevsnow+1:nlevsoil) ! layer depth for snow & soil (m)
471 real(r8),intent(inout) :: dz(1,-nlevsnow+1:nlevsoil) ! layer thickness for soil or snow (m)
472 real(r8),intent(inout) :: t_soisno(1,-nlevsnow+1:nlevsoil) ! soil (or snow) temperature (Kelvin)
473 real(r8),intent(inout) :: t_lake(1,nlevlake) ! lake temperature (Kelvin)
474 integer ,intent(inout) :: snl(1) ! number of snow layers
475 real(r8),intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) ! liquid water (kg/m2)
476 real(r8),intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) ! ice lens (kg/m2)
477 real(r8),intent(inout) :: savedtke1(1) ! top level eddy conductivity from previous timestep (W/m.K)
478 real(r8),intent(inout) :: zi(1,-nlevsnow+0:nlevsoil) ! interface level below a "z" level (m)
479 real(r8),intent(inout) :: lake_icefrac(1,nlevlake) ! mass fraction of lake layer that is frozen
482 !out:
483 real(r8),intent(out) :: eflx_gnet(1) !net heat flux into ground (W/m**2)
484 real(r8),intent(out) :: eflx_lwrad_net(1) ! net infrared (longwave) rad (W/m**2) [+ = to atm]
485 real(r8),intent(out) :: eflx_sh_tot(1) ! total sensible heat flux (W/m**2) [+ to atm]
486 real(r8),intent(out) :: eflx_lh_tot(1) ! total latent heat flux (W/m8*2) [+ to atm]
487 real(r8),intent(out) :: t_ref2m(1) ! 2 m height surface air temperature (Kelvin)
488 real(r8),intent(out) :: q_ref2m(1) ! 2 m height surface specific humidity (kg/kg)
489 real(r8),intent(out) :: taux(1) ! wind (shear) stress: e-w (kg/m/s**2)
490 real(r8),intent(out) :: tauy(1) ! wind (shear) stress: n-s (kg/m/s**2)
491 real(r8),intent(out) :: ram1(1) ! aerodynamical resistance (s/m)
492 ! for calculation of decay of eddy diffusivity with depth
493 ! Change the type variable to pass back to WRF.
494 real(r8),intent(out) :: z0mg(1) ! roughness length over ground, momentum (m(
497 !local output
499 real(r8) :: begwb(1) ! water mass begining of the time step
500 real(r8) :: t_veg(1) ! vegetation temperature (Kelvin)
501 real(r8) :: eflx_soil_grnd(1) ! soil heat flux (W/m**2) [+ = into soil]
502 real(r8) :: eflx_lh_grnd(1) ! ground evaporation heat flux (W/m**2) [+ to atm]
503 real(r8) :: eflx_sh_grnd(1) ! sensible heat flux from ground (W/m**2) [+ to atm]
504 real(r8) :: eflx_lwrad_out(1) ! emitted infrared (longwave) radiation (W/m**2)
505 real(r8) :: qflx_evap_tot(1) ! qflx_evap_soi + qflx_evap_veg + qflx_tran_veg
506 real(r8) :: qflx_evap_soi(1) ! soil evaporation (mm H2O/s) (+ = to atm)
507 real(r8) :: qflx_prec_grnd(1) ! water onto ground including canopy runoff [kg/(m2 s)]
508 real(r8) :: forc_snow(1) ! snow rate [mm/s]
509 real(r8) :: forc_rain(1) ! rain rate [mm/s]
510 real(r8) :: ws(1) ! surface friction velocity (m/s)
511 real(r8) :: ks(1) ! coefficient passed to ShalLakeTemperature
512 real(r8) :: qflx_snomelt(1) !snow melt (mm H2O /s) tem(out),snowwater(in)
513 integer :: imelt(1,-nlevsnow+1:nlevsoil) !flag for melting (=1), freezing (=2), Not=0 (new)
514 real(r8) :: endwb(1) ! water mass end of the time step
515 real(r8) :: snowage(1) ! non dimensional snow age [-]
516 real(r8) :: snowice(1) ! average snow ice lens
517 real(r8) :: snowliq(1) ! average snow liquid water
518 real(r8) :: t_snow(1) ! vertically averaged snow temperature
519 real(r8) :: qflx_drain(1) ! sub-surface runoff (mm H2O /s)
520 real(r8) :: qflx_surf(1) ! surface runoff (mm H2O /s)
521 real(r8) :: qflx_infl(1) ! infiltration (mm H2O /s)
522 real(r8) :: qflx_qrgwl(1) ! qflx_surf at glaciers, wetlands, lakes
523 real(r8) :: qcharge(1) ! aquifer recharge rate (mm/s)
524 real(r8) :: qflx_snowcap(1) ! excess precipitation due to snow capping (mm H2O /s) [+]
525 real(r8) :: qflx_snowcap_col(1) ! excess precipitation due to snow capping (mm H2O /s) [+]
526 real(r8) :: qflx_snow_grnd_pft(1) ! snow on ground after interception (mm H2O/s) [+]
527 real(r8) :: qflx_snow_grnd_col(1) ! snow on ground after interception (mm H2O/s) [+]
528 real(r8) :: qflx_rain_grnd(1) ! rain on ground after interception (mm H2O/s) [+]
529 real(r8) :: frac_iceold(1,-nlevsnow+1:nlevsoil) ! fraction of ice relative to the tot water
530 real(r8) :: qflx_evap_tot_col(1) !pft quantity averaged to the column (assuming one pft)
531 real(r8) :: soilalpha(1) !factor that reduces ground saturated specific humidity (-)
532 real(r8) :: zwt(1) !water table depth
533 real(r8) :: fcov(1) !fractional area with water table at surface
534 real(r8) :: rootr_column(1,1:nlevsoil) !effective fraction of roots in each soil layer
535 real(r8) :: qflx_evap_grnd(1) ! ground surface evaporation rate (mm H2O/s) [+]
536 real(r8) :: qflx_sub_snow(1) ! sublimation rate from snow pack (mm H2O /s) [+]
537 real(r8) :: qflx_dew_snow(1) ! surface dew added to snow pack (mm H2O /s) [+]
538 real(r8) :: qflx_dew_grnd(1) ! ground surface dew formation (mm H2O /s) [+]
539 real(r8) :: qflx_rain_grnd_col(1) !rain on ground after interception (mm H2O/s) [+]
542 ! lat = lat*pie/180 ! [radian]
544 if (prec(1)> 0.) then
545 if ( forc_t(1) > (tfrz + tcrit)) then
546 forc_rain(1) = prec(1)
547 forc_snow(1) = 0.
548 ! flfall(1) = 1.
549 else
550 forc_rain(1) = 0.
551 forc_snow(1) = prec(1)
553 ! if ( forc_t(1) <= tfrz) then
554 ! flfall(1) = 0.
555 ! else if ( forc_t(1) <= tfrz+2.) then
556 ! flfall(1) = -54.632 + 0.2 * forc_t(1)
557 ! else
558 ! flfall(1) = 0.4
559 endif
560 else
561 forc_rain(1) = 0.
562 forc_snow(1) = 0.
563 ! flfall(1) = 1.
564 endif
566 CALL ShalLakeFluxes(forc_t,forc_pbot,forc_psrf,forc_hgt,forc_hgt_q, & !i
567 forc_hgt_t,forc_hgt_u,forc_q, &
568 forc_u,forc_v,forc_lwrad,forc_snow, &
569 forc_rain,t_grnd,h2osno,snowdp,sabg,lat, &
570 dz,dz_lake,t_soisno,t_lake,snl,h2osoi_liq, &
571 h2osoi_ice,savedtke1, &
572 qflx_prec_grnd,qflx_evap_soi,qflx_evap_tot, & !o
573 eflx_sh_grnd,eflx_lwrad_out,eflx_lwrad_net, &
574 eflx_soil_grnd,eflx_sh_tot,eflx_lh_tot, &
575 eflx_lh_grnd,t_veg,t_ref2m,q_ref2m,taux,tauy, &
576 ram1,ws,ks,eflx_gnet,z0mg)
579 CALL ShalLakeTemperature(t_grnd,h2osno,sabg,dz,dz_lake,z,zi, & !i
580 z_lake,ws,ks,snl,eflx_gnet,lakedepth, &
581 lake_icefrac,snowdp, & !i&o
582 eflx_sh_grnd,eflx_sh_tot,eflx_soil_grnd, & !o
583 t_lake,t_soisno,h2osoi_liq, &
584 h2osoi_ice,savedtke1, &
585 frac_iceold,qflx_snomelt,imelt)
589 CALL ShalLakeHydrology(dz_lake,forc_rain,forc_snow, & !i
590 begwb,qflx_evap_tot,forc_t,do_capsnow, &
591 t_grnd,qflx_evap_soi, &
592 qflx_snomelt,imelt,frac_iceold, & !i add by guhp
593 z,dz,zi,snl,h2osno,snowdp,lake_icefrac,t_lake, & !i&o
594 endwb,snowage,snowice,snowliq,t_snow, & !o
595 t_soisno,h2osoi_ice,h2osoi_liq,h2osoi_vol, &
596 qflx_drain,qflx_surf,qflx_infl,qflx_qrgwl, &
597 qcharge,qflx_prec_grnd,qflx_snowcap, &
598 qflx_snowcap_col,qflx_snow_grnd_pft, &
599 qflx_snow_grnd_col,qflx_rain_grnd, &
600 qflx_evap_tot_col,soilalpha,zwt,fcov, &
601 rootr_column,qflx_evap_grnd,qflx_sub_snow, &
602 qflx_dew_snow,qflx_dew_grnd,qflx_rain_grnd_col)
604 !==================================================================================
605 ! !DESCRIPTION:
606 ! Calculation of Shallow Lake Hydrology. Full hydrology of snow layers is
607 ! done. However, there is no infiltration, and the water budget is balanced with
609 END SUBROUTINE LakeMain
612 SUBROUTINE ShalLakeFluxes(forc_t,forc_pbot,forc_psrf,forc_hgt,forc_hgt_q, & !i
613 forc_hgt_t,forc_hgt_u,forc_q, &
614 forc_u,forc_v,forc_lwrad,forc_snow, &
615 forc_rain,t_grnd,h2osno,snowdp,sabg,lat, &
616 dz,dz_lake,t_soisno,t_lake,snl,h2osoi_liq, &
617 h2osoi_ice,savedtke1, &
618 qflx_prec_grnd,qflx_evap_soi,qflx_evap_tot, & !o
619 eflx_sh_grnd,eflx_lwrad_out,eflx_lwrad_net, &
620 eflx_soil_grnd,eflx_sh_tot,eflx_lh_tot, &
621 eflx_lh_grnd,t_veg,t_ref2m,q_ref2m,taux,tauy, &
622 ram1,ws,ks,eflx_gnet,z0mg)
623 !==============================================================================
624 ! DESCRIPTION:
625 ! Calculates lake temperatures and surface fluxes for shallow lakes.
626 !
627 ! Shallow lakes have variable depth, possible snow layers above, freezing & thawing of lake water,
628 ! and soil layers with active temperature and gas diffusion below.
629 !
630 ! WARNING: This subroutine assumes lake columns have one and only one pft.
631 !
632 ! REVISION HISTORY:
633 ! Created by Zack Subin, 2009
634 ! Reedited by Hongping Gu, 2010
635 !==============================================================================
637 ! implicit none
639 implicit none
641 !in:
643 real(r8),intent(in) :: forc_t(1) ! atmospheric temperature (Kelvin)
644 real(r8),intent(in) :: forc_pbot(1) ! atmospheric pressure (Pa)
645 real(r8),intent(in) :: forc_psrf(1) ! atmospheric surface pressure (Pa)
646 real(r8),intent(in) :: forc_hgt(1) ! atmospheric reference height (m)
647 real(r8),intent(in) :: forc_hgt_q(1) ! observational height of humidity [m]
648 real(r8),intent(in) :: forc_hgt_t(1) ! observational height of temperature [m]
649 real(r8),intent(in) :: forc_hgt_u(1) ! observational height of wind [m]
650 real(r8),intent(in) :: forc_q(1) ! atmospheric specific humidity (kg/kg)
651 real(r8),intent(in) :: forc_u(1) ! atmospheric wind speed in east direction (m/s)
652 real(r8),intent(in) :: forc_v(1) ! atmospheric wind speed in north direction (m/s)
653 real(r8),intent(in) :: forc_lwrad(1) ! downward infrared (longwave) radiation (W/m**2)
654 ! real(r8),intent(in) :: forc_rho(1) ! density (kg/m**3)
655 real(r8),intent(in) :: forc_snow(1) ! snow rate [mm/s]
656 real(r8),intent(in) :: forc_rain(1) ! rain rate [mm/s]
657 real(r8),intent(in) :: h2osno(1) ! snow water (mm H2O)
658 real(r8),intent(in) :: snowdp(1) ! snow height (m)
659 real(r8),intent(in) :: sabg(1) ! solar radiation absorbed by ground (W/m**2)
660 real(r8),intent(in) :: lat(1) ! latitude (radians)
661 real(r8),intent(in) :: dz(1,-nlevsnow+1:nlevsoil) ! layer thickness for soil or snow (m)
662 real(r8),intent(in) :: dz_lake(1,nlevlake) ! layer thickness for lake (m)
663 real(r8),intent(in) :: t_soisno(1,-nlevsnow+1:nlevsoil) ! soil (or snow) temperature (Kelvin)
664 real(r8),intent(in) :: t_lake(1,nlevlake) ! lake temperature (Kelvin)
665 integer ,intent(in) :: snl(1) ! number of snow layers
666 real(r8),intent(in) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) ! liquid water (kg/m2)
667 real(r8),intent(in) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) ! ice lens (kg/m2)
668 real(r8),intent(in) :: savedtke1(1) ! top level eddy conductivity from previous timestep (W/m.K)
670 !inout:
671 real(r8),intent(inout) :: t_grnd(1) ! ground temperature (Kelvin)
672 !out:
673 real(r8),intent(out):: qflx_prec_grnd(1) ! water onto ground including canopy runoff [kg/(m2 s)]
674 real(r8),intent(out):: qflx_evap_soi(1) ! soil evaporation (mm H2O/s) (+ = to atm)
675 real(r8),intent(out):: qflx_evap_tot(1) ! qflx_evap_soi + qflx_evap_veg + qflx_tran_veg
676 real(r8),intent(out):: eflx_sh_grnd(1) ! sensible heat flux from ground (W/m**2) [+ to atm]
677 real(r8),intent(out):: eflx_lwrad_out(1) ! emitted infrared (longwave) radiation (W/m**2)
678 real(r8),intent(out):: eflx_lwrad_net(1) ! net infrared (longwave) rad (W/m**2) [+ = to atm]
679 real(r8),intent(out):: eflx_soil_grnd(1) ! soil heat flux (W/m**2) [+ = into soil]
680 real(r8),intent(out):: eflx_sh_tot(1) ! total sensible heat flux (W/m**2) [+ to atm]
681 real(r8),intent(out):: eflx_lh_tot(1) ! total latent heat flux (W/m8*2) [+ to atm]
682 real(r8),intent(out):: eflx_lh_grnd(1) ! ground evaporation heat flux (W/m**2) [+ to atm]
683 real(r8),intent(out):: t_veg(1) ! vegetation temperature (Kelvin)
684 real(r8),intent(out):: t_ref2m(1) ! 2 m height surface air temperature (Kelvin)
685 real(r8),intent(out):: q_ref2m(1) ! 2 m height surface specific humidity (kg/kg)
686 real(r8),intent(out):: taux(1) ! wind (shear) stress: e-w (kg/m/s**2)
687 real(r8),intent(out):: tauy(1) ! wind (shear) stress: n-s (kg/m/s**2)
688 real(r8),intent(out):: ram1(1) ! aerodynamical resistance (s/m)
689 real(r8),intent(out):: ws(1) ! surface friction velocity (m/s)
690 real(r8),intent(out):: ks(1) ! coefficient passed to ShalLakeTemperature
691 ! for calculation of decay of eddy diffusivity with depth
692 real(r8),intent(out):: eflx_gnet(1) !net heat flux into ground (W/m**2)
693 ! Change the type variable to pass back to WRF.
694 real(r8),intent(out):: z0mg(1) ! roughness length over ground, momentum (m(
698 !OTHER LOCAL VARIABLES:
700 integer , parameter :: islak = 2 ! index of lake, 1 = deep lake, 2 = shallow lake
701 integer , parameter :: niters = 3 ! maximum number of iterations for surface temperature
702 real(r8), parameter :: beta1 = 1._r8 ! coefficient of convective velocity (in computing W_*) [-]
703 real(r8), parameter :: emg = 0.97_r8 ! ground emissivity (0.97 for snow)
704 real(r8), parameter :: zii = 1000._r8! convective boundary height [m]
705 real(r8), parameter :: tdmax = 277._r8 ! temperature of maximum water density
706 real(r8) :: forc_th(1) ! atmospheric potential temperature (Kelvin)
707 real(r8) :: forc_vp(1) !atmospheric vapor pressure (Pa)
708 real(r8) :: forc_rho(1) ! density (kg/m**3)
709 integer :: i,fc,fp,g,c,p ! do loop or array index
710 integer :: fncopy ! number of values in pft filter copy
711 integer :: fnold ! previous number of pft filter values
712 integer :: fpcopy(num_shlakep) ! pft filter copy for iteration loop
713 integer :: iter ! iteration index
714 integer :: nmozsgn(lbp:ubp) ! number of times moz changes sign
715 integer :: jtop(lbc:ubc) ! top level for each column (no longer all 1)
716 ! real(r8) :: dtime ! land model time step (sec)
717 real(r8) :: ax ! used in iteration loop for calculating t_grnd (numerator of NR solution)
718 real(r8) :: bx ! used in iteration loop for calculating t_grnd (denomin. of NR solution)
719 real(r8) :: degdT ! d(eg)/dT
720 real(r8) :: dqh(lbp:ubp) ! diff of humidity between ref. height and surface
721 real(r8) :: dth(lbp:ubp) ! diff of virtual temp. between ref. height and surface
722 real(r8) :: dthv ! diff of vir. poten. temp. between ref. height and surface
723 real(r8) :: dzsur(lbc:ubc) ! 1/2 the top layer thickness (m)
724 real(r8) :: eg ! water vapor pressure at temperature T [pa]
725 real(r8) :: htvp(lbc:ubc) ! latent heat of vapor of water (or sublimation) [j/kg]
726 real(r8) :: obu(lbp:ubp) ! monin-obukhov length (m)
727 real(r8) :: obuold(lbp:ubp) ! monin-obukhov length of previous iteration
728 real(r8) :: qsatg(lbc:ubc) ! saturated humidity [kg/kg]
729 real(r8) :: qsatgdT(lbc:ubc) ! d(qsatg)/dT
730 real(r8) :: qstar ! moisture scaling parameter
731 real(r8) :: ram(lbp:ubp) ! aerodynamical resistance [s/m]
732 real(r8) :: rah(lbp:ubp) ! thermal resistance [s/m]
733 real(r8) :: raw(lbp:ubp) ! moisture resistance [s/m]
734 real(r8) :: stftg3(lbp:ubp) ! derivative of fluxes w.r.t ground temperature
735 real(r8) :: temp1(lbp:ubp) ! relation for potential temperature profile
736 real(r8) :: temp12m(lbp:ubp) ! relation for potential temperature profile applied at 2-m
737 real(r8) :: temp2(lbp:ubp) ! relation for specific humidity profile
738 real(r8) :: temp22m(lbp:ubp) ! relation for specific humidity profile applied at 2-m
739 real(r8) :: tgbef(lbc:ubc) ! initial ground temperature
740 real(r8) :: thm(lbc:ubc) ! intermediate variable (forc_t+0.0098*forc_hgt_t)
741 real(r8) :: thv(lbc:ubc) ! virtual potential temperature (kelvin)
742 real(r8) :: thvstar ! virtual potential temperature scaling parameter
743 real(r8) :: tksur ! thermal conductivity of snow/soil (w/m/kelvin)
744 real(r8) :: tsur ! top layer temperature
745 real(r8) :: tstar ! temperature scaling parameter
746 real(r8) :: um(lbp:ubp) ! wind speed including the stablity effect [m/s]
747 real(r8) :: ur(lbp:ubp) ! wind speed at reference height [m/s]
748 real(r8) :: ustar(lbp:ubp) ! friction velocity [m/s]
749 real(r8) :: wc ! convective velocity [m/s]
750 real(r8) :: zeta ! dimensionless height used in Monin-Obukhov theory
751 real(r8) :: zldis(lbp:ubp) ! reference height "minus" zero displacement height [m]
752 real(r8) :: displa(lbp:ubp) ! displacement (always zero) [m]
753 ! real(r8) :: z0mg(lbp:ubp) ! roughness length over ground, momentum [m]
754 real(r8) :: z0hg(lbp:ubp) ! roughness length over ground, sensible heat [m]
755 real(r8) :: z0qg(lbp:ubp) ! roughness length over ground, latent heat [m]
756 real(r8) :: beta(2) ! fraction solar rad absorbed at surface: depends on lake type
757 real(r8) :: u2m ! 2 m wind speed (m/s)
758 real(r8) :: u10(1) ! 10-m wind (m/s) (for dust model)
759 real(r8) :: fv(1) ! friction velocity (m/s) (for dust model)
761 real(r8) :: fm(lbp:ubp) ! needed for BGC only to diagnose 10m wind speed
762 real(r8) :: bw ! partial density of water (ice + liquid)
763 real(r8) :: t_grnd_temp ! Used in surface flux correction over frozen ground
764 real(r8) :: betaprime(lbc:ubc) ! Effective beta: 1 for snow layers, beta(islak) otherwise
765 character*256 :: message
766 ! This assumes all radiation is absorbed in the top snow layer and will need
767 ! to be changed for CLM 4.
768 !
769 ! Constants for lake temperature model
770 !
771 data beta/0.4_r8, 0.4_r8/ ! (deep lake, shallow lake)
772 ! This is the energy absorbed at the lake surface if no snow.
773 ! data za /0.6_r8, 0.5_r8/
774 ! data eta /0.1_r8, 0.5_r8/
775 !-----------------------------------------------------------------------
778 ! dtime = get_step_size()
780 ! Begin calculations
782 !dir$ concurrent
783 !cdir nodep
784 forc_th(1) = forc_t(1) * (forc_psrf(1)/ forc_pbot(1))**(rair/cpair)
785 forc_vp(1) = forc_q(1) * forc_pbot(1)/ (0.622 + 0.378 * forc_q(1))
786 forc_rho(1) = (forc_pbot(1) - 0.378 * forc_vp(1)) / (rair * forc_t(1))
788 do fc = 1, num_shlakec
789 c = filter_shlakec(fc)
790 g = cgridcell(c)
792 ! Surface temperature and fluxes
794 ! Find top layer
795 if (snl(c) > 0 .or. snl(c) < -5) then
796 WRITE(message,*) 'snl is not defined in ShalLakeFluxesMod'
797 CALL wrf_message(message)
798 CALL wrf_error_fatal("snl: out of range value")
799 end if
800 ! if (snl(c) /= 0) then
801 ! write(6,*)'snl is not equal to zero in ShalLakeFluxesMod'
802 ! call endrun()
803 ! end if
804 jtop(c) = snl(c) + 1
807 if (snl(c) < 0) then
808 betaprime(c) = 1._r8 !Assume all solar rad. absorbed at the surface of the top snow layer.
809 dzsur(c) = dz(c,jtop(c))/2._r8
810 else
811 betaprime(c) = beta(islak)
812 dzsur(c) = dz_lake(c,1)/2._r8
813 end if
814 ! Originally this was 1*dz, but shouldn't it be 1/2?
816 ! Saturated vapor pressure, specific humidity and their derivatives
817 ! at lake surface
819 call QSat(t_grnd(c), forc_pbot(g), eg, degdT, qsatg(c), qsatgdT(c))
821 ! Potential, virtual potential temperature, and wind speed at the
822 ! reference height
824 thm(c) = forc_t(g) + 0.0098_r8*forc_hgt_t(g) ! intermediate variable
825 thv(c) = forc_th(g)*(1._r8+0.61_r8*forc_q(g)) ! virtual potential T
826 end do
828 !dir$ concurrent
829 !cdir nodep
830 do fp = 1, num_shlakep
831 p = filter_shlakep(fp)
832 c = pcolumn(p)
833 g = pgridcell(p)
835 nmozsgn(p) = 0
836 obuold(p) = 0._r8
837 displa(p) = 0._r8
839 ! Roughness lengths
842 ! changed by Hongping Gu
843 ! if (t_grnd(c) >= tfrz) then ! for unfrozen lake
844 ! z0mg(p) = 0.01_r8
845 ! else ! for frozen lake
846 ! ! Is this okay even if it is snow covered? What is the roughness over
847 ! non-veg. snow?
848 ! z0mg(p) = 0.04_r8
849 ! end if
851 if (t_grnd(c) >= tfrz) then ! for unfrozen lake
852 z0mg(p) = 0.001_r8 !original 0.01
853 else if(snl(c) == 0 ) then ! for frozen lake
854 ! Is this okay even if it is snow covered? What is the roughness over
855 ! non-veg. snow?
856 z0mg(p) = 0.005_r8 !original 0.04, now for frozen lake without snow
857 else ! for frozen lake with snow
858 z0mg(p) = 0.0024_r8
859 end if
864 z0hg(p) = z0mg(p)
865 z0qg(p) = z0mg(p)
867 ! Latent heat
869 #if (defined PERGRO)
870 htvp(c) = hvap
871 #else
872 if (t_grnd(c) > tfrz) then
873 htvp(c) = hvap
874 else
875 htvp(c) = hsub
876 end if
877 #endif
878 ! Zack Subin, 3/26/09: Shouldn't this be the ground temperature rather than the air temperature above?
879 ! I'll change it for now.
881 ! Initialize stability variables
883 ur(p) = max(1.0_r8,sqrt(forc_u(g)*forc_u(g)+forc_v(g)*forc_v(g)))
884 dth(p) = thm(c)-t_grnd(c)
885 dqh(p) = forc_q(g)-qsatg(c)
886 dthv = dth(p)*(1._r8+0.61_r8*forc_q(g))+0.61_r8*forc_th(g)*dqh(p)
887 zldis(p) = forc_hgt_u(g) - 0._r8
889 ! Initialize Monin-Obukhov length and wind speed
891 call MoninObukIni(ur(p), thv(c), dthv, zldis(p), z0mg(p), um(p), obu(p))
893 end do
895 iter = 1
896 fncopy = num_shlakep
897 fpcopy(1:num_shlakep) = filter_shlakep(1:num_shlakep)
899 ! Begin stability iteration
901 ITERATION : do while (iter <= niters .and. fncopy > 0)
903 ! Determine friction velocity, and potential temperature and humidity
904 ! profiles of the surface boundary layer
906 call FrictionVelocity(pgridcell,forc_hgt,forc_hgt_u, & !i
907 forc_hgt_t,forc_hgt_q, & !i
908 lbp, ubp, fncopy, fpcopy, & !i
909 displa, z0mg, z0hg, z0qg, & !i
910 obu, iter, ur, um, & !i
911 ustar,temp1, temp2, temp12m, temp22m, & !o
912 u10,fv, & !o
913 fm) !i&o
915 !dir$ concurrent
916 !cdir nodep
917 do fp = 1, fncopy
918 p = fpcopy(fp)
919 c = pcolumn(p)
920 g = pgridcell(p)
922 tgbef(c) = t_grnd(c)
923 if (t_grnd(c) > tfrz .and. t_lake(c,1) > tfrz .and. snl(c) == 0) then
924 tksur = savedtke1(c)
925 ! Set this to the eddy conductivity from the last
926 ! timestep, as the molecular conductivity will be orders of magnitude too small.
927 ! Will have to deal with first timestep.
928 tsur = t_lake(c,1)
929 else if (snl(c) == 0) then !frozen but no snow layers
930 tksur = tkice
931 tsur = t_lake(c,1)
932 else
933 !Need to calculate thermal conductivity of the top snow layer
934 bw = (h2osoi_ice(c,jtop(c))+h2osoi_liq(c,jtop(c)))/dz(c,jtop(c))
935 tksur = tkairc + (7.75e-5_r8 *bw + 1.105e-6_r8*bw*bw)*(tkice-tkairc)
936 tsur = t_soisno(c,jtop(c))
937 end if
939 ! Determine aerodynamic resistances
941 ram(p) = 1._r8/(ustar(p)*ustar(p)/um(p))
942 rah(p) = 1._r8/(temp1(p)*ustar(p))
943 raw(p) = 1._r8/(temp2(p)*ustar(p))
944 ram1(p) = ram(p) !pass value to global variable
946 ! Get derivative of fluxes with respect to ground temperature
948 stftg3(p) = emg*sb*tgbef(c)*tgbef(c)*tgbef(c)
950 ! Changed surface temperature from t_lake(c,1) to tsur.
951 ! Also adjusted so that if there are snow layers present, all radiation is absorbed in the top layer.
952 ax = betaprime(c)*sabg(p) + emg*forc_lwrad(g) + 3._r8*stftg3(p)*tgbef(c) &
953 + forc_rho(g)*cpair/rah(p)*thm(c) &
954 - htvp(c)*forc_rho(g)/raw(p)*(qsatg(c)-qsatgdT(c)*tgbef(c) - forc_q(g)) &
955 + tksur*tsur/dzsur(c)
956 !Changed sabg(p) and to betaprime(c)*sabg(p).
957 bx = 4._r8*stftg3(p) + forc_rho(g)*cpair/rah(p) &
958 + htvp(c)*forc_rho(g)/raw(p)*qsatgdT(c) + tksur/dzsur(c)
960 t_grnd(c) = ax/bx
962 ! Update htvp
963 #ifndef PERGRO
964 if (t_grnd(c) > tfrz) then
965 htvp(c) = hvap
966 else
967 htvp(c) = hsub
968 end if
969 #endif
971 ! Surface fluxes of momentum, sensible and latent heat
972 ! using ground temperatures from previous time step
974 eflx_sh_grnd(p) = forc_rho(g)*cpair*(t_grnd(c)-thm(c))/rah(p)
975 qflx_evap_soi(p) = forc_rho(g)*(qsatg(c)+qsatgdT(c)*(t_grnd(c)-tgbef(c))-forc_q(g))/raw(p)
977 ! Re-calculate saturated vapor pressure, specific humidity and their
978 ! derivatives at lake surface
980 call QSat(t_grnd(c), forc_pbot(g), eg, degdT, qsatg(c), qsatgdT(c))
982 dth(p)=thm(c)-t_grnd(c)
983 dqh(p)=forc_q(g)-qsatg(c)
985 tstar = temp1(p)*dth(p)
986 qstar = temp2(p)*dqh(p)
988 thvstar=tstar*(1._r8+0.61_r8*forc_q(g)) + 0.61_r8*forc_th(g)*qstar
989 zeta=zldis(p)*vkc * grav*thvstar/(ustar(p)**2*thv(c))
991 if (zeta >= 0._r8) then !stable
992 zeta = min(2._r8,max(zeta,0.01_r8))
993 um(p) = max(ur(p),0.1_r8)
994 else !unstable
995 zeta = max(-100._r8,min(zeta,-0.01_r8))
996 wc = beta1*(-grav*ustar(p)*thvstar*zii/thv(c))**0.333_r8
997 um(p) = sqrt(ur(p)*ur(p)+wc*wc)
998 end if
999 obu(p) = zldis(p)/zeta
1001 if (obuold(p)*obu(p) < 0._r8) nmozsgn(p) = nmozsgn(p)+1
1003 obuold(p) = obu(p)
1005 end do ! end of filtered pft loop
1007 iter = iter + 1
1008 if (iter <= niters ) then
1009 ! Rebuild copy of pft filter for next pass through the ITERATION loop
1011 fnold = fncopy
1012 fncopy = 0
1013 do fp = 1, fnold
1014 p = fpcopy(fp)
1015 if (nmozsgn(p) < 3) then
1016 fncopy = fncopy + 1
1017 fpcopy(fncopy) = p
1018 end if
1019 end do ! end of filtered pft loop
1020 end if
1022 end do ITERATION ! end of stability iteration
1024 !dir$ concurrent
1025 !cdir nodep
1026 do fp = 1, num_shlakep
1027 p = filter_shlakep(fp)
1028 c = pcolumn(p)
1029 g = pgridcell(p)
1031 ! If there is snow on the ground and t_grnd > tfrz: reset t_grnd = tfrz.
1032 ! Re-evaluate ground fluxes.
1033 ! h2osno > 0.5 prevents spurious fluxes.
1034 ! note that qsatg and qsatgdT should be f(tgbef) (PET: not sure what this
1035 ! comment means)
1036 ! Zack Subin, 3/27/09: Since they are now a function of whatever t_grnd was before cooling
1037 ! to freezing temperature, then this value should be used in the derivative correction term.
1038 ! Should this happen if the lake temperature is below freezing, too? I'll assume that for now.
1039 ! Also, allow convection if ground temp is colder than lake but warmer than 4C, or warmer than
1040 ! lake which is warmer than freezing but less than 4C.
1041 !#ifndef SHLAKETEST
1042 if ( (h2osno(c) > 0.5_r8 .or. t_lake(c,1) <= tfrz) .and. t_grnd(c) > tfrz) then
1043 !#else
1044 ! if ( t_lake(c,1) <= tfrz .and. t_grnd(c) > tfrz) then
1045 !#endif
1046 t_grnd_temp = t_grnd(c)
1047 t_grnd(c) = tfrz
1048 eflx_sh_grnd(p) = forc_rho(g)*cpair*(t_grnd(c)-thm(c))/rah(p)
1049 qflx_evap_soi(p) = forc_rho(g)*(qsatg(c)+qsatgdT(c)*(t_grnd(c)-t_grnd_temp) - forc_q(g))/raw(p)
1050 else if ( (t_lake(c,1) > t_grnd(c) .and. t_grnd(c) > tdmax) .or. &
1051 (t_lake(c,1) < t_grnd(c) .and. t_lake(c,1) > tfrz .and. t_grnd(c) < tdmax) ) then
1052 ! Convective mixing will occur at surface
1053 t_grnd_temp = t_grnd(c)
1054 t_grnd(c) = t_lake(c,1)
1055 eflx_sh_grnd(p) = forc_rho(g)*cpair*(t_grnd(c)-thm(c))/rah(p)
1056 qflx_evap_soi(p) = forc_rho(g)*(qsatg(c)+qsatgdT(c)*(t_grnd(c)-t_grnd_temp) - forc_q(g))/raw(p)
1057 end if
1059 ! Update htvp
1060 #ifndef PERGRO
1061 if (t_grnd(c) > tfrz) then
1062 htvp(c) = hvap
1063 else
1064 htvp(c) = hsub
1065 end if
1066 #endif
1068 ! Net longwave from ground to atmosphere
1070 ! eflx_lwrad_out(p) = (1._r8-emg)*forc_lwrad(g) + stftg3(p)*(-3._r8*tgbef(c)+4._r8*t_grnd(c))
1071 ! What is tgbef doing in this equation? Can't it be exact now? --Zack Subin, 4/14/09
1072 eflx_lwrad_out(p) = (1._r8-emg)*forc_lwrad(g) + emg*sb*t_grnd(c)**4
1074 ! Ground heat flux
1076 eflx_soil_grnd(p) = sabg(p) + forc_lwrad(g) - eflx_lwrad_out(p) - &
1077 eflx_sh_grnd(p) - htvp(c)*qflx_evap_soi(p)
1078 !Why is this sabg(p) and not beta*sabg(p)??
1079 !I've kept this as the incorrect sabg so that the energy balance check will be correct.
1080 !This is the effective energy flux into the ground including the lake solar absorption
1081 !below the surface. The variable eflx_gnet will be used to pass the actual heat flux
1082 !from the ground interface into the lake.
1084 taux(p) = -forc_rho(g)*forc_u(g)/ram(p)
1085 tauy(p) = -forc_rho(g)*forc_v(g)/ram(p)
1087 eflx_sh_tot(p) = eflx_sh_grnd(p)
1088 qflx_evap_tot(p) = qflx_evap_soi(p)
1089 eflx_lh_tot(p) = htvp(c)*qflx_evap_soi(p)
1090 eflx_lh_grnd(p) = htvp(c)*qflx_evap_soi(p)
1091 #if (defined LAKEDEBUG)
1092 write(message,*) 'c, sensible heat = ', c, eflx_sh_tot(p), 'latent heat = ', eflx_lh_tot(p) &
1093 , 'ground temp = ', t_grnd(c), 'h2osno = ', h2osno(c)
1094 CALL wrf_message(message)
1095 if (abs(eflx_sh_tot(p)) > 1500 .or. abs(eflx_lh_tot(p)) > 1500) then
1096 write(message,*)'WARNING: SH, LH = ', eflx_sh_tot(p), eflx_lh_tot(p)
1097 CALL wrf_message(message)
1098 end if
1099 if (abs(eflx_sh_tot(p)) > 10000 .or. abs(eflx_lh_tot(p)) > 10000 &
1100 .or. abs(t_grnd(c)-288)>200 ) CALL wrf_error_fatal ( 't_grnd is out of range' )
1101 #endif
1102 ! 2 m height air temperature
1103 t_ref2m(p) = thm(c) + temp1(p)*dth(p)*(1._r8/temp12m(p) - 1._r8/temp1(p))
1105 ! 2 m height specific humidity
1106 q_ref2m(p) = forc_q(g) + temp2(p)*dqh(p)*(1._r8/temp22m(p) - 1._r8/temp2(p))
1108 ! Energy residual used for melting snow
1109 ! Effectively moved to ShalLakeTemp
1111 ! Prepare for lake layer temperature calculations below
1112 ! fin(c) = betaprime * sabg(p) + forc_lwrad(g) - (eflx_lwrad_out(p) + &
1113 ! eflx_sh_tot(p) + eflx_lh_tot(p))
1114 ! NOW this is just the net ground heat flux calculated below.
1116 eflx_gnet(p) = betaprime(c) * sabg(p) + forc_lwrad(g) - (eflx_lwrad_out(p) + &
1117 eflx_sh_tot(p) + eflx_lh_tot(p))
1118 ! This is the actual heat flux from the ground interface into the lake, not including
1119 ! the light that penetrates the surface.
1121 ! u2m = max(1.0_r8,ustar(p)/vkc*log(2._r8/z0mg(p)))
1122 ! u2 often goes below 1 m/s; it seems like the only reason for this minimum is to
1123 ! keep it from being zero in the ks equation below; 0.1 m/s is a better limit for
1124 ! stable conditions --ZS
1125 u2m = max(0.1_r8,ustar(p)/vkc*log(2._r8/z0mg(p)))
1127 ws(c) = 1.2e-03_r8 * u2m
1128 ks(c) = 6.6_r8*sqrt(abs(sin(lat(g))))*(u2m**(-1.84_r8))
1130 end do
1132 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1133 ! End of surface flux relevant code in original BiogeophysicsLakeMod until history loop.
1135 ! The following are needed for global average on history tape.
1137 !dir$ concurrent
1138 !cdir nodep
1139 do fp = 1, num_shlakep
1140 p = filter_shlakep(fp)
1141 c = pcolumn(p)
1142 g = pgridcell(p)
1143 ! t_veg(p) = forc_t(g)
1144 !This is an odd choice, since elsewhere t_veg = t_grnd for bare ground.
1145 !Zack Subin, 4/09
1146 t_veg(p) = t_grnd(c)
1147 eflx_lwrad_net(p) = eflx_lwrad_out(p) - forc_lwrad(g)
1148 qflx_prec_grnd(p) = forc_rain(g) + forc_snow(g)
1149 end do
1151 END SUBROUTINE ShalLakeFluxes
1153 SUBROUTINE ShalLakeTemperature(t_grnd,h2osno,sabg,dz,dz_lake,z,zi, & !i
1154 z_lake,ws,ks,snl,eflx_gnet,lakedepth, &
1155 lake_icefrac,snowdp, & !i&o
1156 eflx_sh_grnd,eflx_sh_tot,eflx_soil_grnd, & !o
1157 t_lake,t_soisno,h2osoi_liq, &
1158 h2osoi_ice,savedtke1, &
1159 frac_iceold,qflx_snomelt,imelt)
1160 !=======================================================================================================
1161 ! !DESCRIPTION:
1162 ! Calculates temperatures in the 20-25 layer column of (possible) snow,
1163 ! lake water, and soil beneath lake.
1164 ! Snow and soil temperatures are determined as in SoilTemperature, except
1165 ! for appropriate boundary conditions at the top of the snow (the flux is fixed
1166 ! to be the ground heat flux calculated in ShalLakeFluxes), the bottom of the snow
1167 ! (adjacent to top lake layer), and the top of the soil (adjacent to the bottom
1168 ! lake layer). Also, the soil is assumed to be always fully saturated (ShalLakeHydrology
1169 ! will have to insure this). The whole column is solved simultaneously as one tridiagonal matrix.
1170 ! Lake temperatures are determined from the Hostetler model as before, except now:
1171 ! i) Lake water layers can freeze by any fraction and release latent heat; thermal
1172 ! and mechanical properties are adjusted for ice fraction.
1173 ! ii) Convective mixing (though not eddy diffusion) still occurs for frozen lakes.
1174 ! iii) No sunlight is absorbed in the lake if there are snow layers.
1175 ! iv) Light is allowed to reach the top soil layer (where it is assumed to be completely absorbed).
1176 ! v) Lakes have variable depth, set ultimately in surface data set but now in initShalLakeMod.
1177 !
1178 ! Eddy + molecular diffusion:
1179 ! d ts d d ts 1 ds
1180 ! ---- = -- [(km + ke) ----] + -- --
1181 ! dt dz dz cw dz
1182 !
1183 ! where: ts = temperature (kelvin)
1184 ! t = time (s)
1185 ! z = depth (m)
1186 ! km = molecular diffusion coefficient (m**2/s)
1187 ! ke = eddy diffusion coefficient (m**2/s)
1188 ! cw = heat capacity (j/m**3/kelvin)
1189 ! s = heat source term (w/m**2)
1190 !
1191 ! Shallow lakes are allowed to have variable depth, set in _____.
1192 !
1193 ! For shallow lakes: ke > 0 if unfrozen,
1194 ! and convective mixing occurs WHETHER OR NOT frozen. (See e.g. Martynov...)
1195 !
1196 ! Use the Crank-Nicholson method to set up tridiagonal system of equations to
1197 ! solve for ts at time n+1, where the temperature equation for layer i is
1198 ! r_i = a_i [ts_i-1] n+1 + b_i [ts_i] n+1 + c_i [ts_i+1] n+1
1199 !
1200 ! The solution conserves energy as:
1201 !
1202 ! [For lake layers]
1203 ! cw*([ts( 1)] n+1 - [ts( 1)] n)*dz( 1)/dt + ... +
1204 ! cw*([ts(nlevlake)] n+1 - [ts(nlevlake)] n)*dz(nlevlake)/dt = fin
1205 ! But now there is phase change, so cv is not constant and there is
1206 ! latent heat.
1207 !
1208 ! where:
1209 ! [ts] n = old temperature (kelvin)
1210 ! [ts] n+1 = new temperature (kelvin)
1211 ! fin = heat flux into lake (w/m**2)
1212 ! = betaprime*sabg + forc_lwrad - eflx_lwrad_out - eflx_sh_tot - eflx_lh_tot
1213 ! (This is now the same as the ground heat flux.)
1214 ! + phi(1) + ... + phi(nlevlake) + phi(top soil level)
1215 ! betaprime = beta(islak) for no snow layers, and 1 for snow layers.
1216 ! This assumes all radiation is absorbed in the top snow layer and will need
1217 ! to be changed for CLM 4.
1218 !
1219 ! WARNING: This subroutine assumes lake columns have one and only one pft.
1220 !
1221 ! Outline:
1222 ! 1!) Initialization
1223 ! 2!) Lake density
1224 ! 3!) Diffusivity
1225 ! 4!) Heat source term from solar radiation penetrating lake
1226 ! 5!) Set thermal props and find initial energy content
1227 ! 6!) Set up vectors for tridiagonal matrix solution
1228 ! 7!) Solve tridiagonal and back-substitute
1229 ! 8!) (Optional) Do first energy check using temperature change at constant heat capacity.
1230 ! 9!) Phase change
1231 ! 9.5!) (Optional) Do second energy check using temperature change and latent heat, considering changed heat capacity.
1232 ! Also do soil water balance check.
1233 !10!) Convective mixing
1234 !11!) Do final energy check to detect small numerical errors (especially from convection)
1235 ! and dump small imbalance into sensible heat, or pass large errors to BalanceCheckMod for abort.
1236 !
1237 ! REVISION HISTORY:
1238 ! Created by Zack Subin, 2009.
1239 ! Reedited by Hongping Gu, 2010.
1240 !=========================================================================================================
1243 ! use TridiagonalMod , only : Tridiagonal
1245 implicit none
1247 !in:
1248 real(r8), intent(in) :: t_grnd(1) ! ground temperature (Kelvin)
1249 real(r8), intent(inout) :: h2osno(1) ! snow water (mm H2O)
1250 real(r8), intent(in) :: sabg(1) ! solar radiation absorbed by ground (W/m**2)
1251 real(r8), intent(in) :: dz(1,-nlevsnow + 1:nlevsoil) ! layer thickness for snow & soil (m)
1252 real(r8), intent(in) :: dz_lake(1,nlevlake) ! layer thickness for lake (m)
1253 real(r8), intent(in) :: z(1,-nlevsnow+1:nlevsoil) ! layer depth for snow & soil (m)
1254 real(r8), intent(in) :: zi(1,-nlevsnow+0:nlevsoil) ! interface level below a "z" level (m)
1255 ! the other z and dz variables
1256 real(r8), intent(in) :: z_lake(1,nlevlake) ! layer depth for lake (m)
1257 real(r8), intent(in) :: ws(1) ! surface friction velocity (m/s)
1258 real(r8), intent(in) :: ks(1) ! coefficient passed to ShalLakeTemperature
1259 ! for calculation of decay of eddy diffusivity with depth
1260 integer , intent(in) :: snl(1) ! negative of number of snow layers
1261 real(r8), intent(inout) :: eflx_gnet(1) ! net heat flux into ground (W/m**2) at the surface interface
1262 real(r8), intent(in) :: lakedepth(1) ! column lake depth (m)
1264 ! real(r8), intent(in) :: watsat(1,nlevsoil) ! volumetric soil water at saturation (porosity)
1265 real(r8), intent(inout) :: snowdp(1) !snow height (m)
1266 !out:
1268 real(r8), intent(out) :: eflx_sh_grnd(1) ! sensible heat flux from ground (W/m**2) [+ to atm]
1269 real(r8), intent(out) :: eflx_sh_tot(1) ! total sensible heat flux (W/m**2) [+ to atm]
1270 real(r8), intent(out) :: eflx_soil_grnd(1) ! heat flux into snow / lake (W/m**2) [+ = into soil]
1271 ! Here this includes the whole lake radiation absorbed.
1272 #if (defined SHLAKETEST)
1273 real(r8), intent(out) :: qmelt(1) ! snow melt [mm/s] [temporary]
1274 #endif
1275 real(r8), intent(inout) :: t_lake(1,nlevlake) ! lake temperature (Kelvin)
1276 real(r8), intent(inout) :: t_soisno(1,-nlevsnow+1:nlevsoil) ! soil (or snow) temperature (Kelvin)
1277 real(r8), intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) ! liquid water (kg/m2) [for snow & soil layers]
1278 real(r8), intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) ! ice lens (kg/m2) [for snow & soil layers]
1279 real(r8), intent(inout) :: lake_icefrac(1,nlevlake) ! mass fraction of lake layer that is frozen
1280 real(r8), intent(out) :: savedtke1(1) ! top level thermal conductivity (W/mK)
1281 real(r8), intent(out) :: frac_iceold(1,-nlevsnow+1:nlevsoil) ! fraction of ice relative to the tot water
1282 real(r8), intent(out) :: qflx_snomelt(1) !snow melt (mm H2O /s)
1283 integer, intent(out) :: imelt(1,-nlevsnow+1:nlevsoil) !flag for melting (=1), freezing (=2), Not=0 (new)
1286 ! OTHER LOCAL VARIABLES:
1288 integer , parameter :: islak = 2 ! index of lake, 1 = deep lake, 2 = shallow lake
1289 real(r8), parameter :: p0 = 1._r8 ! neutral value of turbulent prandtl number
1290 integer :: i,j,fc,fp,g,c,p ! do loop or array index
1291 ! real(r8) :: dtime ! land model time step (sec)
1292 real(r8) :: beta(2) ! fraction solar rad absorbed at surface: depends on lake type
1293 real(r8) :: za(2) ! base of surface absorption layer (m): depends on lake type
1294 real(r8) :: eta(2) ! light extinction coefficient (/m): depends on lake type
1295 real(r8) :: cwat ! specific heat capacity of water (j/m**3/kelvin)
1296 real(r8) :: cice_eff ! effective heat capacity of ice (using density of
1297 ! water because layer depth is not adjusted when freezing
1298 real(r8) :: cfus ! effective heat of fusion per unit volume
1299 ! using water density as above
1300 real(r8) :: km ! molecular diffusion coefficient (m**2/s)
1301 real(r8) :: tkice_eff ! effective conductivity since layer depth is constant
1302 real(r8) :: a(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) ! "a" vector for tridiagonal matrix
1303 real(r8) :: b(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) ! "b" vector for tridiagonal matrix
1304 real(r8) :: c1(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) ! "c" vector for tridiagonal matrix
1305 real(r8) :: r(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) ! "r" vector for tridiagonal solution
1306 real(r8) :: rhow(lbc:ubc,nlevlake) ! density of water (kg/m**3)
1307 real(r8) :: phi(lbc:ubc,nlevlake) ! solar radiation absorbed by layer (w/m**2)
1308 real(r8) :: kme(lbc:ubc,nlevlake) ! molecular + eddy diffusion coefficient (m**2/s)
1309 real(r8) :: rsfin ! relative flux of solar radiation into layer
1310 real(r8) :: rsfout ! relative flux of solar radiation out of layer
1311 real(r8) :: phi_soil(lbc:ubc) ! solar radiation into top soil layer (W/m**2)
1312 real(r8) :: ri ! richardson number
1313 real(r8) :: fin(lbc:ubc) ! net heat flux into lake at ground interface (w/m**2)
1314 real(r8) :: ocvts(lbc:ubc) ! (cwat*(t_lake[n ])*dz
1315 real(r8) :: ncvts(lbc:ubc) ! (cwat*(t_lake[n+1])*dz
1316 real(r8) :: ke ! eddy diffusion coefficient (m**2/s)
1317 real(r8) :: zin ! depth at top of layer (m)
1318 real(r8) :: zout ! depth at bottom of layer (m)
1319 real(r8) :: drhodz ! d [rhow] /dz (kg/m**4)
1320 real(r8) :: n2 ! brunt-vaisala frequency (/s**2)
1321 real(r8) :: num ! used in calculating ri
1322 real(r8) :: den ! used in calculating ri
1323 real(r8) :: tav_froz(lbc:ubc) ! used in aver temp for convectively mixed layers (C)
1324 real(r8) :: tav_unfr(lbc:ubc) ! "
1325 real(r8) :: nav(lbc:ubc) ! used in aver temp for convectively mixed layers
1326 real(r8) :: phidum ! temporary value of phi
1327 real(r8) :: iceav(lbc:ubc) ! used in calc aver ice for convectively mixed layers
1328 real(r8) :: qav(lbc:ubc) ! used in calc aver heat content for conv. mixed layers
1329 integer :: jtop(lbc:ubc) ! top level for each column (no longer all 1)
1330 real(r8) :: cv (lbc:ubc,-nlevsnow+1:nlevsoil) !heat capacity of soil/snow [J/(m2 K)]
1331 real(r8) :: tk (lbc:ubc,-nlevsnow+1:nlevsoil) !thermal conductivity of soil/snow [W/(m K)]
1332 !(at interface below, except for j=0)
1333 real(r8) :: cv_lake (lbc:ubc,1:nlevlake) !heat capacity [J/(m2 K)]
1334 real(r8) :: tk_lake (lbc:ubc,1:nlevlake) !thermal conductivity at layer node [W/(m K)]
1335 real(r8) :: cvx (lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !heat capacity for whole column [J/(m2 K)]
1336 real(r8) :: tkix(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !thermal conductivity at layer interfaces
1337 !for whole column [W/(m K)]
1338 real(r8) :: tx(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) ! temperature of whole column [K]
1339 real(r8) :: tktopsoillay(lbc:ubc) ! thermal conductivity [W/(m K)]
1340 real(r8) :: fnx(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !heat diffusion through the layer interface below [W/m2]
1341 real(r8) :: phix(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !solar source term for whole column [W/m**2]
1342 real(r8) :: zx(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !interface depth (+ below surface) for whole column [m]
1343 real(r8) :: dzm !used in computing tridiagonal matrix [m]
1344 real(r8) :: dzp !used in computing tridiagonal matrix [m]
1345 integer :: jprime ! j - nlevlake
1346 real(r8) :: factx(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !coefficient used in computing tridiagonal matrix
1347 real(r8) :: t_lake_bef(lbc:ubc,1:nlevlake) !beginning lake temp for energy conservation check [K]
1348 real(r8) :: t_soisno_bef(lbc:ubc,-nlevsnow+1:nlevsoil) !beginning soil temp for E cons. check [K]
1349 real(r8) :: lhabs(lbc:ubc) ! total per-column latent heat abs. from phase change (J/m^2)
1350 real(r8) :: esum1(lbc:ubc) ! temp for checking energy (J/m^2)
1351 real(r8) :: esum2(lbc:ubc) ! ""
1352 real(r8) :: zsum(lbc:ubc) ! temp for putting ice at the top during convection (m)
1353 real(r8) :: wsum(lbc:ubc) ! temp for checking water (kg/m^2)
1354 real(r8) :: wsum_end(lbc:ubc) ! temp for checking water (kg/m^2)
1355 real(r8) :: errsoi(1) ! soil/lake energy conservation error (W/m**2)
1356 real(r8) :: eflx_snomelt(1) !snow melt heat flux (W/m**2)
1357 CHARACTER*256 :: message
1358 !
1359 ! Constants for lake temperature model
1360 !
1361 data beta/0.4_r8, 0.4_r8/ ! (deep lake, shallow lake)
1362 data za /0.6_r8, 0.6_r8/
1363 ! For now, keep beta and za for shallow lake the same as deep lake, until better data is found.
1364 ! It looks like eta is key and that larger values give better results for shallow lakes. Use
1365 ! empirical expression from Hakanson (below). This is still a very unconstrained parameter
1366 ! that deserves more attention.
1367 ! Some radiation will be allowed to reach the soil.
1368 !-----------------------------------------------------------------------
1371 ! 1!) Initialization
1372 ! Determine step size
1374 ! dtime = get_step_size()
1376 ! Initialize constants
1377 cwat = cpliq*denh2o ! water heat capacity per unit volume
1378 cice_eff = cpice*denh2o !use water density because layer depth is not adjusted
1379 !for freezing
1380 cfus = hfus*denh2o ! latent heat per unit volume
1381 tkice_eff = tkice * denice/denh2o !effective conductivity since layer depth is constant
1382 km = tkwat/cwat ! a constant (molecular diffusivity)
1384 ! Begin calculations
1386 !dir$ concurrent
1387 !cdir nodep
1388 do fc = 1, num_shlakec
1389 c = filter_shlakec(fc)
1391 ! Initialize Ebal quantities computed below
1393 ocvts(c) = 0._r8
1394 ncvts(c) = 0._r8
1395 esum1(c) = 0._r8
1396 esum2(c) = 0._r8
1398 end do
1400 ! Initialize set of previous time-step variables as in DriverInit,
1401 ! which is currently not called over lakes. This has to be done
1402 ! here because phase change will occur in this routine.
1403 ! Ice fraction of snow at previous time step
1405 do j = -nlevsnow+1,0
1406 !dir$ concurrent
1407 !cdir nodep
1408 do fc = 1, num_shlakec
1409 c = filter_shlakec(fc)
1410 if (j >= snl(c) + 1) then
1411 frac_iceold(c,j) = h2osoi_ice(c,j)/(h2osoi_liq(c,j)+h2osoi_ice(c,j))
1412 end if
1413 end do
1414 end do
1416 ! Sum soil water.
1417 do j = 1, nlevsoil
1418 !dir$ concurrent
1419 !cdir nodep
1420 do fc = 1, num_shlakec
1421 c = filter_shlakec(fc)
1422 if (j == 1) wsum(c) = 0._r8
1423 wsum(c) = wsum(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
1424 end do
1425 end do
1427 !dir$ concurrent
1428 !cdir nodep
1429 do fp = 1, num_shlakep
1430 p = filter_shlakep(fp)
1431 c = pcolumn(p)
1434 ! Prepare for lake layer temperature calculations below
1436 ! fin(c) = betaprime * sabg(p) + forc_lwrad(g) - (eflx_lwrad_out(p) + &
1437 ! eflx_sh_tot(p) + eflx_lh_tot(p))
1438 ! fin(c) now passed from ShalLakeFluxes as eflx_gnet
1439 fin(c) = eflx_gnet(p)
1441 end do
1443 ! 2!) Lake density
1445 do j = 1, nlevlake
1446 !dir$ concurrent
1447 !cdir nodep
1448 do fc = 1, num_shlakec
1449 c = filter_shlakec(fc)
1450 rhow(c,j) = (1._r8 - lake_icefrac(c,j)) * &
1451 1000._r8*( 1.0_r8 - 1.9549e-05_r8*(abs(t_lake(c,j)-277._r8))**1.68_r8 ) &
1452 + lake_icefrac(c,j)*denice
1453 ! Allow for ice fraction; assume constant ice density.
1454 ! Is this the right weighted average?
1455 ! Using this average will make sure that surface ice is treated properly during
1456 ! convective mixing.
1457 end do
1458 end do
1460 ! 3!) Diffusivity and implied thermal "conductivity" = diffusivity * cwat
1461 do j = 1, nlevlake-1
1462 !dir$ prefervector
1463 !dir$ concurrent
1464 !cdir nodep
1465 do fc = 1, num_shlakec
1466 c = filter_shlakec(fc)
1467 drhodz = (rhow(c,j+1)-rhow(c,j)) / (z_lake(c,j+1)-z_lake(c,j))
1468 n2 = grav / rhow(c,j) * drhodz
1469 ! Fixed sign error here: our z goes up going down into the lake, so no negative
1470 ! sign is needed to make this positive unlike in Hostetler. --ZS
1471 num = 40._r8 * n2 * (vkc*z_lake(c,j))**2
1472 den = max( (ws(c)**2) * exp(-2._r8*ks(c)*z_lake(c,j)), 1.e-10_r8 )
1473 ri = ( -1._r8 + sqrt( max(1._r8+num/den, 0._r8) ) ) / 20._r8
1474 if (t_grnd(c) > tfrz .and. t_lake(c,1) > tfrz .and. snl(c) == 0) then
1475 ! ke = vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._r8+37._r8*ri*ri)
1477 if( t_lake(c,1) > 277.15_r8 ) then
1478 if (lakedepth(c) > 15.0 ) then
1479 ke = 1.e+2_r8*vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._r8+37._r8*ri*ri)
1480 else
1481 ke = vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._r8+37._r8*ri*ri)
1482 endif
1483 else
1484 if (lakedepth(c) > 15.0 ) then
1485 if (lakedepth(c) > 150.0 ) then
1486 ke = 1.e+5_r8*vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._r8+37._r8*ri*ri)
1487 else
1488 ke =1.e+4_r8*vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._r8+37._r8*ri*ri)
1489 end if
1490 else
1491 ke = vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._r8+37._r8*ri*ri)
1492 endif
1493 end if
1495 kme(c,j) = km + ke
1496 tk_lake(c,j) = kme(c,j)*cwat
1497 ! If there is some ice in this layer (this should rarely happen because the surface
1498 ! is unfrozen and it will be unstable), still use the cwat to get out the tk b/c the eddy
1499 ! diffusivity equation assumes water.
1500 else
1501 kme(c,j) = km
1502 tk_lake(c,j) = tkwat*tkice_eff / ( (1._r8-lake_icefrac(c,j))*tkice_eff &
1503 + tkwat*lake_icefrac(c,j) )
1504 ! Assume the resistances add as for the calculation of conductivities at layer interfaces.
1505 end if
1506 end do
1507 end do
1509 !dir$ concurrent
1510 !cdir nodep
1511 do fc = 1, num_shlakec
1512 c = filter_shlakec(fc)
1514 j = nlevlake
1515 kme(c,nlevlake) = kme(c,nlevlake-1)
1517 if (t_grnd(c) > tfrz .and. t_lake(c,1) > tfrz .and. snl(c) == 0) then
1518 tk_lake(c,j) = tk_lake(c,j-1)
1519 else
1520 tk_lake(c,j) = tkwat*tkice_eff / ( (1._r8-lake_icefrac(c,j))*tkice_eff &
1521 + tkwat*lake_icefrac(c,j) )
1522 end if
1524 ! Use in surface flux calculation for next timestep.
1525 savedtke1(c) = kme(c,1)*cwat ! Will only be used if unfrozen
1526 ! set number of column levels for use by Tridiagonal below
1527 jtop(c) = snl(c) + 1
1528 end do
1530 ! 4!) Heat source term: unfrozen lakes only
1531 do j = 1, nlevlake
1532 !dir$ concurrent
1533 !cdir nodep
1534 do fp = 1, num_shlakep
1535 p = filter_shlakep(fp)
1536 c = pcolumn(p)
1538 ! Set eta(:), the extinction coefficient, according to L Hakanson, Aquatic Sciences, 1995
1539 ! (regression of Secchi Depth with lake depth for small glacial basin lakes), and the
1540 ! Poole & Atkins expression for extinction coeffient of 1.7 / Secchi Depth (m).
1541 #ifndef ETALAKE
1542 eta(:) = 1.1925_r8*lakedepth(c)**(-0.424)
1543 #else
1544 eta(:) = ETALAKE
1545 #endif
1547 zin = z_lake(c,j) - 0.5_r8*dz_lake(c,j)
1548 zout = z_lake(c,j) + 0.5_r8*dz_lake(c,j)
1549 rsfin = exp( -eta(islak)*max( zin-za(islak),0._r8 ) )
1550 rsfout = exp( -eta(islak)*max( zout-za(islak),0._r8 ) )
1552 ! Let rsfout for bottom layer go into soil.
1553 ! This looks like it should be robust even for pathological cases,
1554 ! like lakes thinner than za.
1555 if (t_grnd(c) > tfrz .and. t_lake(c,1) > tfrz .and. snl(c) == 0) then
1556 phidum = (rsfin-rsfout) * sabg(p) * (1._r8-beta(islak))
1557 if (j == nlevlake) then
1558 phi_soil(c) = rsfout * sabg(p) * (1._r8-beta(islak))
1559 end if
1560 else if (j == 1 .and. snl(c) == 0) then !if frozen but no snow layers
1561 phidum = sabg(p) * (1._r8-beta(islak))
1562 else !radiation absorbed at surface
1563 phidum = 0._r8
1564 if (j == nlevlake) phi_soil(c) = 0._r8
1565 end if
1566 phi(c,j) = phidum
1568 end do
1569 end do
1571 ! 5!) Set thermal properties and check initial energy content.
1573 ! For lake
1574 do j = 1, nlevlake
1575 !dir$ concurrent
1576 !cdir nodep
1577 do fc = 1, num_shlakec
1578 c = filter_shlakec(fc)
1580 cv_lake(c,j) = dz_lake(c,j) * (cwat*(1._r8-lake_icefrac(c,j)) + cice_eff*lake_icefrac(c,j))
1581 end do
1582 end do
1584 ! For snow / soil
1585 call SoilThermProp_Lake (snl,dz,zi,z,t_soisno,h2osoi_liq,h2osoi_ice, &
1586 tk, cv, tktopsoillay)
1588 ! Sum cv*t_lake for energy check
1589 ! Include latent heat term, and correction for changing heat capacity with phase change.
1591 ! This will need to be over all soil / lake / snow layers. Lake is below.
1592 do j = 1, nlevlake
1593 !dir$ concurrent
1594 !cdir nodep
1595 do fc = 1, num_shlakec
1596 c = filter_shlakec(fc)
1598 ! ocvts(c) = ocvts(c) + cv_lake(c,j)*t_lake(c,j) &
1599 ocvts(c) = ocvts(c) + cv_lake(c,j)*(t_lake(c,j)-tfrz) &
1600 + cfus*dz_lake(c,j)*(1._r8-lake_icefrac(c,j)) !&
1601 ! + (cwat-cice_eff)*lake_icefrac(c)*tfrz*dz_lake(c,j) !enthalpy reconciliation term
1602 t_lake_bef(c,j) = t_lake(c,j)
1603 end do
1604 end do
1606 ! Now do for soil / snow layers
1607 do j = -nlevsnow + 1, nlevsoil
1608 !dir$ concurrent
1609 !cdir nodep
1610 do fc = 1, num_shlakec
1611 c = filter_shlakec(fc)
1613 if (j >= jtop(c)) then
1614 ! ocvts(c) = ocvts(c) + cv(c,j)*t_soisno(c,j) &
1615 ocvts(c) = ocvts(c) + cv(c,j)*(t_soisno(c,j)-tfrz) &
1616 + hfus*h2osoi_liq(c,j) !&
1617 ! + (cpliq-cpice)*h2osoi_ice(c,j)*tfrz !enthalpy reconciliation term
1618 if (j == 1 .and. h2osno(c) > 0._r8 .and. j == jtop(c)) then
1619 ocvts(c) = ocvts(c) - h2osno(c)*hfus
1620 end if
1621 t_soisno_bef(c,j) = t_soisno(c,j)
1622 if(abs(t_soisno(c,j)-288) > 150) then
1623 WRITE( message,* ) 'WARNING: Extreme t_soisno at c, level',c, j
1624 CALL wrf_error_fatal ( message )
1625 endif
1626 end if
1627 end do
1628 end do
1630 !!!!!!!!!!!!!!!!!!!
1631 ! 6!) Set up vector r and vectors a, b, c1 that define tridiagonal matrix
1633 ! Heat capacity and resistance of snow without snow layers (<1cm) is ignored during diffusion,
1634 ! but its capacity to absorb latent heat may be used during phase change.
1636 ! Set up interface depths, zx, heat capacities, cvx, solar source terms, phix, and temperatures, tx.
1637 do j = -nlevsnow+1, nlevlake+nlevsoil
1638 !dir$ prefervector
1639 !dir$ concurrent
1640 !cdir nodep
1641 do fc = 1,num_shlakec
1642 c = filter_shlakec(fc)
1644 jprime = j - nlevlake
1646 if (j >= jtop(c)) then
1647 if (j < 1) then !snow layer
1648 zx(c,j) = z(c,j)
1649 cvx(c,j) = cv(c,j)
1650 phix(c,j) = 0._r8
1651 tx(c,j) = t_soisno(c,j)
1652 else if (j <= nlevlake) then !lake layer
1653 zx(c,j) = z_lake(c,j)
1654 cvx(c,j) = cv_lake(c,j)
1655 phix(c,j) = phi(c,j)
1656 tx(c,j) = t_lake(c,j)
1657 else !soil layer
1658 zx(c,j) = zx(c,nlevlake) + dz_lake(c,nlevlake)/2._r8 + z(c,jprime)
1659 cvx(c,j) = cv(c,jprime)
1660 if (j == nlevlake + 1) then !top soil layer
1661 phix(c,j) = phi_soil(c)
1662 else !middle or bottom soil layer
1663 phix(c,j) = 0._r8
1664 end if
1665 tx(c,j) = t_soisno(c,jprime)
1666 end if
1667 end if
1669 end do
1670 end do
1672 ! Determine interface thermal conductivities, tkix
1674 do j = -nlevsnow+1, nlevlake+nlevsoil
1675 !dir$ prefervector
1676 !dir$ concurrent
1677 !cdir nodep
1678 do fc = 1,num_shlakec
1679 c = filter_shlakec(fc)
1681 jprime = j - nlevlake
1683 if (j >= jtop(c)) then
1684 if (j < 0) then !non-bottom snow layer
1685 tkix(c,j) = tk(c,j)
1686 else if (j == 0) then !bottom snow layer
1687 dzp = zx(c,j+1) - zx(c,j)
1688 tkix(c,j) = tk_lake(c,1)*tk(c,j)*dzp / &
1689 (tk(c,j)*z_lake(c,1) + tk_lake(c,1)*(-z(c,j)) )
1690 ! tk(c,0) is the conductivity at the middle of that layer, as defined in SoilThermProp_Lake
1691 else if (j < nlevlake) then !non-bottom lake layer
1692 tkix(c,j) = ( tk_lake(c,j)*tk_lake(c,j+1) * (dz_lake(c,j+1)+dz_lake(c,j)) ) &
1693 / ( tk_lake(c,j)*dz_lake(c,j+1) + tk_lake(c,j+1)*dz_lake(c,j) )
1694 else if (j == nlevlake) then !bottom lake layer
1695 dzp = zx(c,j+1) - zx(c,j)
1696 tkix(c,j) = (tktopsoillay(c)*tk_lake(c,j)*dzp / &
1697 (tktopsoillay(c)*dz_lake(c,j)/2._r8 + tk_lake(c,j)*z(c,1) ) )
1698 ! tktopsoillay is the conductivity at the middle of that layer, as defined in SoilThermProp_Lake
1699 else !soil layer
1700 tkix(c,j) = tk(c,jprime)
1701 end if
1702 end if
1704 end do
1705 end do
1708 ! Determine heat diffusion through the layer interface and factor used in computing
1709 ! tridiagonal matrix and set up vector r and vectors a, b, c1 that define tridiagonal
1710 ! matrix and solve system
1712 do j = -nlevsnow+1, nlevlake+nlevsoil
1713 !dir$ prefervector
1714 !dir$ concurrent
1715 !cdir nodep
1716 do fc = 1,num_shlakec
1717 c = filter_shlakec(fc)
1718 if (j >= jtop(c)) then
1719 if (j < nlevlake+nlevsoil) then !top or interior layer
1720 factx(c,j) = dtime/cvx(c,j)
1721 fnx(c,j) = tkix(c,j)*(tx(c,j+1)-tx(c,j))/(zx(c,j+1)-zx(c,j))
1722 else !bottom soil layer
1723 factx(c,j) = dtime/cvx(c,j)
1724 fnx(c,j) = 0._r8 !not used
1725 end if
1726 end if
1727 enddo
1728 end do
1730 do j = -nlevsnow+1,nlevlake+nlevsoil
1731 !dir$ prefervector
1732 !dir$ concurrent
1733 !cdir nodep
1734 do fc = 1,num_shlakec
1735 c = filter_shlakec(fc)
1736 if (j >= jtop(c)) then
1737 if (j == jtop(c)) then !top layer
1738 dzp = zx(c,j+1)-zx(c,j)
1739 a(c,j) = 0._r8
1740 b(c,j) = 1+(1._r8-cnfac)*factx(c,j)*tkix(c,j)/dzp
1741 c1(c,j) = -(1._r8-cnfac)*factx(c,j)*tkix(c,j)/dzp
1742 r(c,j) = tx(c,j) + factx(c,j)*( fin(c) + phix(c,j) + cnfac*fnx(c,j) )
1743 else if (j < nlevlake+nlevsoil) then !middle layer
1744 dzm = (zx(c,j)-zx(c,j-1))
1745 dzp = (zx(c,j+1)-zx(c,j))
1746 a(c,j) = - (1._r8-cnfac)*factx(c,j)* tkix(c,j-1)/dzm
1747 b(c,j) = 1._r8+ (1._r8-cnfac)*factx(c,j)*(tkix(c,j)/dzp + tkix(c,j-1)/dzm)
1748 c1(c,j) = - (1._r8-cnfac)*factx(c,j)* tkix(c,j)/dzp
1749 r(c,j) = tx(c,j) + cnfac*factx(c,j)*( fnx(c,j) - fnx(c,j-1) ) + factx(c,j)*phix(c,j)
1750 else !bottom soil layer
1751 dzm = (zx(c,j)-zx(c,j-1))
1752 a(c,j) = - (1._r8-cnfac)*factx(c,j)*tkix(c,j-1)/dzm
1753 b(c,j) = 1._r8+ (1._r8-cnfac)*factx(c,j)*tkix(c,j-1)/dzm
1754 c1(c,j) = 0._r8
1755 r(c,j) = tx(c,j) - cnfac*factx(c,j)*fnx(c,j-1)
1756 end if
1757 end if
1758 enddo
1759 end do
1760 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1763 ! 7!) Solve for tdsolution
1765 call Tridiagonal(lbc, ubc, -nlevsnow + 1, nlevlake + nlevsoil, jtop, num_shlakec, filter_shlakec, &
1766 a, b, c1, r, tx)
1768 ! Set t_soisno and t_lake
1769 do j = -nlevsnow+1, nlevlake + nlevsoil
1770 !dir$ concurrent
1771 !cdir nodep
1772 do fc = 1, num_shlakec
1773 c = filter_shlakec(fc)
1775 jprime = j - nlevlake
1777 ! Don't do anything with invalid snow layers.
1778 if (j >= jtop(c)) then
1779 if (j < 1) then !snow layer
1780 t_soisno(c,j) = tx(c,j)
1781 else if (j <= nlevlake) then !lake layer
1782 t_lake(c,j) = tx(c,j)
1783 else !soil layer
1784 t_soisno(c,jprime) = tx(c,j)
1785 end if
1786 end if
1787 end do
1788 end do
1790 !!!!!!!!!!!!!!!!!!!!!!!
1792 ! 8!) Sum energy content and total energy into lake for energy check. Any errors will be from the
1793 ! Tridiagonal solution.
1795 #if (defined LAKEDEBUG)
1796 do j = 1, nlevlake
1797 !dir$ concurrent
1798 !cdir nodep
1799 do fc = 1, num_shlakec
1800 c = filter_shlakec(fc)
1802 esum1(c) = esum1(c) + (t_lake(c,j)-t_lake_bef(c,j))*cv_lake(c,j)
1803 esum2(c) = esum2(c) + (t_lake(c,j)-tfrz)*cv_lake(c,j)
1804 end do
1805 end do
1807 do j = -nlevsnow+1, nlevsoil
1808 !dir$ concurrent
1809 !cdir nodep
1810 do fc = 1, num_shlakec
1811 c = filter_shlakec(fc)
1813 if (j >= jtop(c)) then
1814 esum1(c) = esum1(c) + (t_soisno(c,j)-t_soisno_bef(c,j))*cv(c,j)
1815 esum2(c) = esum2(c) + (t_soisno(c,j)-tfrz)*cv(c,j)
1816 end if
1817 end do
1818 end do
1820 !dir$ concurrent
1821 !cdir nodep
1822 do fp = 1, num_shlakep
1823 p = filter_shlakep(fp)
1824 c = pcolumn(p)
1825 ! Again assuming only one pft per column
1826 ! esum1(c) = esum1(c) + lhabs(c)
1827 errsoi(c) = esum1(c)/dtime - eflx_soil_grnd(p)
1828 ! eflx_soil_grnd includes all the solar radiation absorbed in the lake,
1829 ! unlike eflx_gnet
1830 if(abs(errsoi(c)) > 1.e-5_r8) then
1831 WRITE( message,* )'Primary soil energy conservation error in shlake &
1832 column during Tridiagonal Solution,', 'error (W/m^2):', c, errsoi(c)
1833 CALL wrf_error_fatal ( message )
1834 end if
1835 end do
1836 ! This has to be done before convective mixing because the heat capacities for each layer
1837 ! will get scrambled.
1839 #endif
1841 !!!!!!!!!!!!!!!!!!!!!!!
1843 ! 9!) Phase change
1844 call PhaseChange_Lake (snl,h2osno,dz,dz_lake, & !i
1845 t_soisno,h2osoi_liq,h2osoi_ice, & !i&o
1846 lake_icefrac,t_lake, snowdp, & !i&o
1847 qflx_snomelt,eflx_snomelt,imelt, & !o
1848 cv, cv_lake, & !i&o
1849 lhabs) !o
1851 !!!!!!!!!!!!!!!!!!!!!!!
1853 ! 9.5!) Second energy check and water check. Now check energy balance before and after phase
1854 ! change, considering the possibility of changed heat capacity during phase change, by
1855 ! using initial heat capacity in the first step, final heat capacity in the second step,
1856 ! and differences from tfrz only to avoid enthalpy correction for (cpliq-cpice)*melt*tfrz.
1857 ! Also check soil water sum.
1859 #if (defined LAKEDEBUG)
1860 do j = 1, nlevlake
1861 !dir$ concurrent
1862 !cdir nodep
1863 do fc = 1, num_shlakec
1864 c = filter_shlakec(fc)
1866 esum2(c) = esum2(c) - (t_lake(c,j)-tfrz)*cv_lake(c,j)
1867 end do
1868 end do
1870 do j = -nlevsnow+1, nlevsoil
1871 !dir$ concurrent
1872 !cdir nodep
1873 do fc = 1, num_shlakec
1874 c = filter_shlakec(fc)
1876 if (j >= jtop(c)) then
1877 esum2(c) = esum2(c) - (t_soisno(c,j)-tfrz)*cv(c,j)
1878 end if
1879 end do
1880 end do
1882 !dir$ concurrent
1883 !cdir nodep
1884 do fp = 1, num_shlakep
1885 p = filter_shlakep(fp)
1886 c = pcolumn(p)
1887 ! Again assuming only one pft per column
1888 esum2(c) = esum2(c) - lhabs(c)
1889 errsoi(c) = esum2(c)/dtime
1890 if(abs(errsoi(c)) > 1.e-5_r8) then
1891 write(message,*)'Primary soil energy conservation error in shlake column during Phase Change, error (W/m^2):', &
1892 c, errsoi(c)
1893 CALL wrf_error_fatal ( message )
1894 end if
1895 end do
1897 ! Check soil water
1898 ! Sum soil water.
1899 do j = 1, nlevsoil
1900 !dir$ concurrent
1901 !cdir nodep
1902 do fc = 1, num_shlakec
1903 c = filter_shlakec(fc)
1904 if (j == 1) wsum_end(c) = 0._r8
1905 wsum_end(c) = wsum_end(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
1906 if (j == nlevsoil) then
1907 if (abs(wsum(c)-wsum_end(c))>1.e-7_r8) then
1908 write(message,*)'Soil water balance error during phase change in ShalLakeTemperature.', &
1909 'column, error (kg/m^2):', c, wsum_end(c)-wsum(c)
1910 CALL wrf_error_fatal ( message )
1911 end if
1912 end if
1913 end do
1914 end do
1916 #endif
1918 !!!!!!!!!!!!!!!!!!!!!!!!!!
1919 ! 10!) Convective mixing: make sure fracice*dz is conserved, heat content c*dz*T is conserved, and
1920 ! all ice ends up at the top. Done over all lakes even if frozen.
1921 ! Either an unstable density profile or ice in a layer below an incompletely frozen layer will trigger.
1923 !Recalculate density
1924 do j = 1, nlevlake
1925 !dir$ concurrent
1926 !cdir nodep
1927 do fc = 1, num_shlakec
1928 c = filter_shlakec(fc)
1929 rhow(c,j) = (1._r8 - lake_icefrac(c,j)) * &
1930 1000._r8*( 1.0_r8 - 1.9549e-05_r8*(abs(t_lake(c,j)-277._r8))**1.68_r8 ) &
1931 + lake_icefrac(c,j)*denice
1932 end do
1933 end do
1935 do j = 1, nlevlake-1
1936 !dir$ concurrent
1937 !cdir nodep
1938 do fc = 1, num_shlakec
1939 c = filter_shlakec(fc)
1940 qav(c) = 0._r8
1941 nav(c) = 0._r8
1942 iceav(c) = 0._r8
1943 end do
1945 do i = 1, j+1
1946 !dir$ concurrent
1947 !cdir nodep
1948 do fc = 1, num_shlakec
1949 c = filter_shlakec(fc)
1950 if (rhow(c,j) > rhow(c,j+1) .or. &
1951 (lake_icefrac(c,j) < 1._r8 .and. lake_icefrac(c,j+1) > 0._r8) ) then
1952 #if (defined LAKEDEBUG)
1953 if (i==1) then
1954 write(message,*), 'Convective Mixing in column ', c, '.'
1955 CALL wrf_message(message)
1956 endif
1957 #endif
1958 qav(c) = qav(c) + dz_lake(c,i)*(t_lake(c,i)-tfrz) * &
1959 ((1._r8 - lake_icefrac(c,i))*cwat + lake_icefrac(c,i)*cice_eff)
1960 ! tav(c) = tav(c) + t_lake(c,i)*dz_lake(c,i)
1961 iceav(c) = iceav(c) + lake_icefrac(c,i)*dz_lake(c,i)
1962 nav(c) = nav(c) + dz_lake(c,i)
1963 end if
1964 end do
1965 end do
1967 !dir$ concurrent
1968 !cdir nodep
1969 do fc = 1, num_shlakec
1970 c = filter_shlakec(fc)
1971 if (rhow(c,j) > rhow(c,j+1) .or. &
1972 (lake_icefrac(c,j) < 1._r8 .and. lake_icefrac(c,j+1) > 0._r8) ) then
1973 qav(c) = qav(c)/nav(c)
1974 iceav(c) = iceav(c)/nav(c)
1975 !If the average temperature is above freezing, put the extra energy into the water.
1976 !If it is below freezing, take it away from the ice.
1977 if (qav(c) > 0._r8) then
1978 tav_froz(c) = 0._r8 !Celsius
1979 tav_unfr(c) = qav(c) / ((1._r8 - iceav(c))*cwat)
1980 else if (qav(c) < 0._r8) then
1981 tav_froz(c) = qav(c) / (iceav(c)*cice_eff)
1982 tav_unfr(c) = 0._r8 !Celsius
1983 else
1984 tav_froz(c) = 0._r8
1985 tav_unfr(c) = 0._r8
1986 end if
1987 end if
1988 end do
1990 do i = 1, j+1
1991 !dir$ concurrent
1992 !cdir nodep
1993 do fc = 1, num_shlakec
1994 c = filter_shlakec(fc)
1995 if (nav(c) > 0._r8) then
1996 ! if(0==1) then
1998 !Put all the ice at the top.!
1999 !If the average temperature is above freezing, put the extra energy into the water.
2000 !If it is below freezing, take it away from the ice.
2001 !For the layer with both ice & water, be careful to use the average temperature
2002 !that preserves the correct total heat content given what the heat capacity of that
2003 !layer will actually be.
2004 if (i == 1) zsum(c) = 0._r8
2005 if ((zsum(c)+dz_lake(c,i))/nav(c) <= iceav(c)) then
2006 lake_icefrac(c,i) = 1._r8
2007 t_lake(c,i) = tav_froz(c) + tfrz
2008 else if (zsum(c)/nav(c) < iceav(c)) then
2009 lake_icefrac(c,i) = (iceav(c)*nav(c) - zsum(c)) / dz_lake(c,i)
2010 ! Find average value that preserves correct heat content.
2011 t_lake(c,i) = ( lake_icefrac(c,i)*tav_froz(c)*cice_eff &
2012 + (1._r8 - lake_icefrac(c,i))*tav_unfr(c)*cwat ) &
2013 / ( lake_icefrac(c,i)*cice_eff + (1-lake_icefrac(c,i))*cwat ) + tfrz
2014 else
2015 lake_icefrac(c,i) = 0._r8
2016 t_lake(c,i) = tav_unfr(c) + tfrz
2017 end if
2018 zsum(c) = zsum(c) + dz_lake(c,i)
2020 rhow(c,i) = (1._r8 - lake_icefrac(c,i)) * &
2021 1000._r8*( 1.0_r8 - 1.9549e-05_r8*(abs(t_lake(c,i)-277._r8))**1.68_r8 ) &
2022 + lake_icefrac(c,i)*denice
2023 end if
2024 end do
2025 end do
2026 end do
2028 !!!!!!!!!!!!!!!!!!!!!!!
2029 ! 11!) Re-evaluate thermal properties and sum energy content.
2030 ! For lake
2031 do j = 1, nlevlake
2032 !dir$ concurrent
2033 !cdir nodep
2034 do fc = 1, num_shlakec
2035 c = filter_shlakec(fc)
2037 cv_lake(c,j) = dz_lake(c,j) * (cwat*(1._r8-lake_icefrac(c,j)) + cice_eff*lake_icefrac(c,j))
2038 #if (defined LAKEDEBUG)
2039 write(message,*)'Lake Ice Fraction, c, level:', c, j, lake_icefrac(c,j)
2040 CALL wrf_message(message)
2041 #endif
2042 end do
2043 end do
2044 ! For snow / soil
2045 ! call SoilThermProp_Lake(lbc, ubc, num_shlakec, filter_shlakec, tk, cv, tktopsoillay)
2046 call SoilThermProp_Lake (snl,dz,zi,z,t_soisno,h2osoi_liq,h2osoi_ice, &
2047 tk, cv, tktopsoillay)
2050 ! Do as above to sum energy content
2051 do j = 1, nlevlake
2052 !dir$ concurrent
2053 !cdir nodep
2054 do fc = 1, num_shlakec
2055 c = filter_shlakec(fc)
2057 ! ncvts(c) = ncvts(c) + cv_lake(c,j)*t_lake(c,j) &
2058 ncvts(c) = ncvts(c) + cv_lake(c,j)*(t_lake(c,j)-tfrz) &
2059 + cfus*dz_lake(c,j)*(1._r8-lake_icefrac(c,j)) !&
2060 ! + (cwat-cice_eff)*lake_icefrac(c)*tfrz*dz_lake(c,j) !enthalpy reconciliation term
2061 fin(c) = fin(c) + phi(c,j)
2062 end do
2063 end do
2065 do j = -nlevsnow + 1, nlevsoil
2066 !dir$ concurrent
2067 !cdir nodep
2068 do fc = 1, num_shlakec
2069 c = filter_shlakec(fc)
2071 if (j >= jtop(c)) then
2072 ! ncvts(c) = ncvts(c) + cv(c,j)*t_soisno(c,j) &
2073 ncvts(c) = ncvts(c) + cv(c,j)*(t_soisno(c,j)-tfrz) &
2074 + hfus*h2osoi_liq(c,j) !&
2075 ! + (cpliq-cpice)*h2osoi_ice(c,j)*tfrz !enthalpy reconciliation term
2076 if (j == 1 .and. h2osno(c) > 0._r8 .and. j == jtop(c)) then
2077 ncvts(c) = ncvts(c) - h2osno(c)*hfus
2078 end if
2079 end if
2080 if (j == 1) fin(c) = fin(c) + phi_soil(c)
2081 end do
2082 end do
2085 ! Check energy conservation.
2087 do fp = 1, num_shlakep
2088 p = filter_shlakep(fp)
2089 c = pcolumn(p)
2090 errsoi(c) = (ncvts(c)-ocvts(c)) / dtime - fin(c)
2091 #ifndef LAKEDEBUG
2092 ! if (abs(errsoi(c)) < 0.10_r8) then ! else send to Balance Check and abort
2093 if (abs(errsoi(c)) < 10._r8) then ! else send to Balance Check and abort
2094 #else
2095 if (abs(errsoi(c)) < 1._r8) then
2096 #endif
2097 eflx_sh_tot(p) = eflx_sh_tot(p) - errsoi(c)
2098 eflx_sh_grnd(p) = eflx_sh_grnd(p) - errsoi(c)
2099 eflx_soil_grnd(p) = eflx_soil_grnd(p) + errsoi(c)
2100 eflx_gnet(p) = eflx_gnet(p) + errsoi(c)
2101 ! if (abs(errsoi(c)) > 1.e-3_r8) then
2102 if (abs(errsoi(c)) > 1.e-1_r8) then
2103 write(message,*)'errsoi incorporated into sensible heat in ShalLakeTemperature: c, (W/m^2):', c, errsoi(c)
2104 CALL wrf_message(message)
2105 end if
2106 errsoi(c) = 0._r8
2107 #if (defined LAKEDEBUG)
2108 else
2109 write(message,*)'Soil Energy Balance Error at column, ', c, 'G, fintotal, column E tendency = ', &
2110 eflx_gnet(p), fin(c), (ncvts(c)-ocvts(c)) / dtime
2111 CALL wrf_message(message)
2112 #endif
2113 end if
2114 end do
2115 ! This loop assumes only one point per column.
2117 end subroutine ShalLakeTemperature
2119 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2120 !-----------------------------------------------------------------------
2121 !BOP
2122 !
2123 ! ROUTINE: SoilThermProp_Lake
2124 !
2125 ! !INTERFACE:
2126 subroutine SoilThermProp_Lake (snl,dz,zi,z,t_soisno,h2osoi_liq,h2osoi_ice, &
2127 tk, cv, tktopsoillay)
2129 !
2130 ! !DESCRIPTION:
2131 ! Calculation of thermal conductivities and heat capacities of
2132 ! snow/soil layers
2133 ! (1) The volumetric heat capacity is calculated as a linear combination
2134 ! in terms of the volumetric fraction of the constituent phases.
2135 !
2136 ! (2) The thermal conductivity of soil is computed from the algorithm of
2137 ! Johansen (as reported by Farouki 1981), and of snow is from the
2138 ! formulation used in SNTHERM (Jordan 1991).
2139 ! The thermal conductivities at the interfaces between two neighboring
2140 ! layers (j, j+1) are derived from an assumption that the flux across
2141 ! the interface is equal to that from the node j to the interface and the
2142 ! flux from the interface to the node j+1.
2143 !
2144 ! For lakes, the proper soil layers (not snow) should always be saturated.
2145 !
2146 ! !USES:
2148 implicit none
2149 !in
2151 integer , intent(in) :: snl(1) ! number of snow layers
2152 ! real(r8), intent(in) :: h2osno(1) ! snow water (mm H2O)
2153 ! real(r8), intent(in) :: watsat(1,nlevsoil) ! volumetric soil water at saturation (porosity)
2154 ! real(r8), intent(in) :: tksatu(1,nlevsoil) ! thermal conductivity, saturated soil [W/m-K]
2155 ! real(r8), intent(in) :: tkmg(1,nlevsoil) ! thermal conductivity, soil minerals [W/m-K]
2156 ! real(r8), intent(in) :: tkdry(1,nlevsoil) ! thermal conductivity, dry soil (W/m/Kelvin)
2157 ! real(r8), intent(in) :: csol(1,nlevsoil) ! heat capacity, soil solids (J/m**3/Kelvin)
2158 real(r8), intent(in) :: dz(1,-nlevsnow+1:nlevsoil) ! layer thickness (m)
2159 real(r8), intent(in) :: zi(1,-nlevsnow+0:nlevsoil) ! interface level below a "z" level (m)
2160 real(r8), intent(in) :: z(1,-nlevsnow+1:nlevsoil) ! layer depth (m)
2161 real(r8), intent(in) :: t_soisno(1,-nlevsnow+1:nlevsoil) ! soil temperature (Kelvin)
2162 real(r8), intent(in) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) ! liquid water (kg/m2)
2163 real(r8), intent(in) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) ! ice lens (kg/m2)
2165 !out
2166 real(r8), intent(out) :: cv(lbc:ubc,-nlevsnow+1:nlevsoil) ! heat capacity [J/(m2 K)]
2167 real(r8), intent(out) :: tk(lbc:ubc,-nlevsnow+1:nlevsoil) ! thermal conductivity [W/(m K)]
2168 real(r8), intent(out) :: tktopsoillay(lbc:ubc) ! thermal conductivity [W/(m K)]
2169 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2170 ! !CALLED FROM:
2171 ! subroutine ShalLakeTemperature in this module.
2172 !
2173 ! !REVISION HISTORY:
2174 ! 15 September 1999: Yongjiu Dai; Initial code
2175 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
2176 ! 2/13/02, Peter Thornton: migrated to new data structures
2177 ! 7/01/03, Mariana Vertenstein: migrated to vector code
2178 ! 4/09, Zack Subin, adjustment for ShalLake code.
2179 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2180 ! !LOCAL VARIABLES:
2181 !
2182 ! local pointers to original implicit in scalars
2183 !
2184 ! integer , pointer :: clandunit(:) ! column's landunit
2185 ! integer , pointer :: ityplun(:) ! landunit type
2186 !
2187 !EOP
2190 ! OTHER LOCAL VARIABLES:
2192 integer :: l,c,j ! indices
2193 integer :: fc ! lake filtered column indices
2194 real(r8) :: bw ! partial density of water (ice + liquid)
2195 real(r8) :: dksat ! thermal conductivity for saturated soil (j/(k s m))
2196 real(r8) :: dke ! kersten number
2197 real(r8) :: fl ! fraction of liquid or unfrozen water to total water
2198 real(r8) :: satw ! relative total water content of soil.
2199 real(r8) :: thk(lbc:ubc,-nlevsnow+1:nlevsoil) ! thermal conductivity of layer
2200 character*256 :: message
2202 ! Thermal conductivity of soil from Farouki (1981)
2204 do j = -nlevsnow+1,nlevsoil
2205 !dir$ concurrent
2206 !cdir nodep
2207 do fc = 1, num_shlakec
2208 c = filter_shlakec(fc)
2210 ! Only examine levels from 1->nlevsoil
2211 if (j >= 1) then
2212 ! l = clandunit(c)
2213 ! if (ityplun(l) /= istwet .AND. ityplun(l) /= istice) then
2214 ! This could be altered later for allowing this to be over glaciers.
2216 ! Soil should be saturated.
2217 #if (defined LAKEDEBUG)
2218 satw = (h2osoi_liq(c,j)/denh2o + h2osoi_ice(c,j)/denice)/(dz(c,j)*watsat(c,j))
2219 ! satw = min(1._r8, satw)
2220 if (satw < 0.999_r8) then
2221 write(message,*)'WARNING: soil layer unsaturated in SoilThermProp_Lake, satw, j = ', satw, j
2222 CALL wrf_error_fatal ( message )
2223 end if
2224 ! Could use denice because if it starts out frozen, the volume of water will go below sat.,
2225 ! since we're not yet doing excess ice.
2226 ! But take care of this in HydrologyLake.
2227 #endif
2228 satw = 1._r8
2229 fl = h2osoi_liq(c,j)/(h2osoi_ice(c,j)+h2osoi_liq(c,j))
2230 if (t_soisno(c,j) >= tfrz) then ! Unfrozen soil
2231 dke = max(0._r8, log10(satw) + 1.0_r8)
2232 dksat = tksatu(c,j)
2233 else ! Frozen soil
2234 dke = satw
2235 dksat = tkmg(c,j)*0.249_r8**(fl*watsat(c,j))*2.29_r8**watsat(c,j)
2236 endif
2237 thk(c,j) = dke*dksat + (1._r8-dke)*tkdry(c,j)
2238 ! else
2239 ! thk(c,j) = tkwat
2240 ! if (t_soisno(c,j) < tfrz) thk(c,j) = tkice
2241 ! endif
2242 endif
2244 ! Thermal conductivity of snow, which from Jordan (1991) pp. 18
2245 ! Only examine levels from snl(c)+1 -> 0 where snl(c) < 1
2246 if (snl(c)+1 < 1 .AND. (j >= snl(c)+1) .AND. (j <= 0)) then
2247 bw = (h2osoi_ice(c,j)+h2osoi_liq(c,j))/dz(c,j)
2248 thk(c,j) = tkairc + (7.75e-5_r8 *bw + 1.105e-6_r8*bw*bw)*(tkice-tkairc)
2249 end if
2251 end do
2252 end do
2254 ! Thermal conductivity at the layer interface
2256 ! Have to correct for the fact that bottom snow layer and top soil layer border lake.
2257 ! For the first case, the snow layer conductivity for the middle of the layer will be returned.
2258 ! Because the interfaces are below the soil layers, the conductivity for the top soil layer
2259 ! will have to be returned separately.
2260 do j = -nlevsnow+1,nlevsoil
2261 !dir$ concurrent
2262 !cdir nodep
2263 do fc = 1,num_shlakec
2264 c = filter_shlakec(fc)
2265 if (j >= snl(c)+1 .AND. j <= nlevsoil-1 .AND. j /= 0) then
2266 tk(c,j) = thk(c,j)*thk(c,j+1)*(z(c,j+1)-z(c,j)) &
2267 /(thk(c,j)*(z(c,j+1)-zi(c,j))+thk(c,j+1)*(zi(c,j)-z(c,j)))
2268 else if (j == 0) then
2269 tk(c,j) = thk(c,j)
2270 else if (j == nlevsoil) then
2271 tk(c,j) = 0._r8
2272 end if
2273 ! For top soil layer.
2274 if (j == 1) tktopsoillay(c) = thk(c,j)
2275 end do
2276 end do
2278 ! Soil heat capacity, from de Vires (1963)
2280 do j = 1, nlevsoil
2281 !dir$ concurrent
2282 !cdir nodep
2283 do fc = 1,num_shlakec
2284 c = filter_shlakec(fc)
2285 ! l = clandunit(c)
2286 ! if (ityplun(l) /= istwet .AND. ityplun(l) /= istice) then
2287 cv(c,j) = csol(c,j)*(1-watsat(c,j))*dz(c,j) + &
2288 (h2osoi_ice(c,j)*cpice + h2osoi_liq(c,j)*cpliq)
2289 ! else
2290 ! cv(c,j) = (h2osoi_ice(c,j)*cpice + h2osoi_liq(c,j)*cpliq)
2291 ! endif
2292 ! if (j == 1) then
2293 ! if (snl(c)+1 == 1 .AND. h2osno(c) > 0._r8) then
2294 ! cv(c,j) = cv(c,j) + cpice*h2osno(c)
2295 ! end if
2296 ! end if
2297 ! Won't worry about heat capacity for thin snow on lake with no snow layers.
2298 enddo
2299 end do
2301 ! Snow heat capacity
2303 do j = -nlevsnow+1,0
2304 !dir$ concurrent
2305 !cdir nodep
2306 do fc = 1,num_shlakec
2307 c = filter_shlakec(fc)
2308 if (snl(c)+1 < 1 .and. j >= snl(c)+1) then
2309 cv(c,j) = cpliq*h2osoi_liq(c,j) + cpice*h2osoi_ice(c,j)
2310 end if
2311 end do
2312 end do
2314 end subroutine SoilThermProp_Lake
2317 !-----------------------------------------------------------------------
2318 !BOP
2319 !
2320 ! ROUTINE: PhaseChange_Lake
2321 !
2322 ! !INTERFACE:
2323 subroutine PhaseChange_Lake (snl,h2osno,dz,dz_lake, & !i
2324 t_soisno,h2osoi_liq,h2osoi_ice, & !i&o
2325 lake_icefrac,t_lake, snowdp, & !i&o
2326 qflx_snomelt,eflx_snomelt,imelt, & !o
2327 cv, cv_lake, & !i&o
2328 lhabs) !o
2329 !=============================================================================================
2330 ! !DESCRIPTION:
2331 ! Calculation of the phase change within snow, soil, & lake layers:
2332 ! (1) Check the conditions for which the phase change may take place,
2333 ! i.e., the layer temperature is great than the freezing point
2334 ! and the ice mass is not equal to zero (i.e. melting),
2335 ! or the layer temperature is less than the freezing point
2336 ! and the liquid water mass is greater than the allowable supercooled
2337 ! (i.e. freezing).
2338 ! (2) Assess the amount of phase change from the energy excess (or deficit)
2339 ! after setting the layer temperature to freezing point, depending on
2340 ! how much water or ice is available.
2341 ! (3) Re-adjust the ice and liquid mass, and the layer temperature: either to
2342 ! the freezing point if enough water or ice is available to fully compensate,
2343 ! or to a remaining temperature.
2344 ! The specific heats are assumed constant. Potential cycling errors resulting from
2345 ! this assumption will be trapped at the end of ShalLakeTemperature.
2346 ! !CALLED FROM:
2347 ! subroutine ShalLakeTemperature in this module
2348 !
2349 ! !REVISION HISTORY:
2350 ! 04/2009 Zack Subin: Initial code
2351 !==============================================================================================
2352 ! !USES:
2353 !
2354 ! !ARGUMENTS:
2355 implicit none
2356 !in:
2358 integer , intent(in) :: snl(1) !number of snow layers
2359 real(r8), intent(inout) :: h2osno(1) !snow water (mm H2O)
2360 real(r8), intent(in) :: dz(1,-nlevsnow+1:nlevsoil) !layer thickness (m)
2361 real(r8), intent(in) :: dz_lake(1,nlevlake) !lake layer thickness (m)
2362 ! Needed in case snow height is less than critical value.
2364 !inout:
2366 real(r8), intent(inout) :: snowdp(1) !snow height (m)
2367 real(r8), intent(inout) :: t_soisno(1,-nlevsnow+1:nlevsoil) !soil temperature (Kelvin)
2368 real(r8), intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) !liquid water (kg/m2)
2369 real(r8), intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) !ice lens (kg/m2)
2370 real(r8), intent(inout) :: lake_icefrac(1,nlevlake) ! mass fraction of lake layer that is frozen
2371 real(r8), intent(inout) :: t_lake(1,nlevlake) ! lake temperature (Kelvin)
2372 !out:
2374 real(r8), intent(out) :: qflx_snomelt(1) !snow melt (mm H2O /s)
2375 real(r8), intent(out) :: eflx_snomelt(1) !snow melt heat flux (W/m**2)
2376 integer, intent(out) :: imelt(1,-nlevsnow+1:nlevsoil) !flag for melting (=1), freezing (=2), Not=0 (new)
2377 !What's the sign of this? Is it just output?
2378 real(r8), intent(inout) :: cv(lbc:ubc,-nlevsnow+1:nlevsoil) ! heat capacity [J/(m2 K)]
2379 real(r8), intent(inout) :: cv_lake (lbc:ubc,1:nlevlake) ! heat capacity [J/(m2 K)]
2380 real(r8), intent(out):: lhabs(lbc:ubc) ! total per-column latent heat abs. (J/m^2)
2383 ! OTHER LOCAL VARIABLES:
2385 integer :: j,c,g !do loop index
2386 integer :: fc !lake filtered column indices
2387 ! real(r8) :: dtime !land model time step (sec)
2388 real(r8) :: heatavail !available energy for melting or freezing (J/m^2)
2389 real(r8) :: heatrem !energy residual or loss after melting or freezing
2390 real(r8) :: melt !actual melting (+) or freezing (-) [kg/m2]
2391 real(r8), parameter :: smallnumber = 1.e-7_r8 !to prevent tiny residuals from rounding error
2392 logical :: dophasechangeflag
2393 !-----------------------------------------------------------------------
2395 ! dtime = get_step_size()
2397 ! Initialization
2399 !dir$ concurrent
2400 !cdir nodep
2401 do fc = 1,num_shlakec
2402 c = filter_shlakec(fc)
2404 qflx_snomelt(c) = 0._r8
2405 eflx_snomelt(c) = 0._r8
2406 lhabs(c) = 0._r8
2407 end do
2409 do j = -nlevsnow+1,0
2410 !dir$ concurrent
2411 !cdir nodep
2412 do fc = 1,num_shlakec
2413 c = filter_shlakec(fc)
2415 if (j >= snl(c) + 1) imelt(c,j) = 0
2416 end do
2417 end do
2419 ! Check for case of snow without snow layers and top lake layer temp above freezing.
2421 !dir$ concurrent
2422 !cdir nodep
2423 do fc = 1,num_shlakec
2424 c = filter_shlakec(fc)
2426 if (snl(c) == 0 .and. h2osno(c) > 0._r8 .and. t_lake(c,1) > tfrz) then
2427 heatavail = (t_lake(c,1) - tfrz) * cv_lake(c,1)
2428 melt = min(h2osno(c), heatavail/hfus)
2429 heatrem = max(heatavail - melt*hfus, 0._r8)
2430 !catch small negative value to keep t at tfrz
2431 t_lake(c,1) = tfrz + heatrem/(cv_lake(c,1))
2432 snowdp(c) = snowdp(c)*(1._r8 - melt/h2osno(c))
2433 h2osno(c) = h2osno(c) - melt
2434 lhabs(c) = lhabs(c) + melt*hfus
2435 qflx_snomelt(c) = qflx_snomelt(c) + melt
2436 ! Prevent tiny residuals
2437 if (h2osno(c) < smallnumber) h2osno(c) = 0._r8
2438 if (snowdp(c) < smallnumber) snowdp(c) = 0._r8
2439 end if
2440 end do
2442 ! Lake phase change
2444 do j = 1,nlevlake
2445 !dir$ concurrent
2446 !cdir nodep
2447 do fc = 1,num_shlakec
2448 c = filter_shlakec(fc)
2450 dophasechangeflag = .false.
2451 if (t_lake(c,j) > tfrz .and. lake_icefrac(c,j) > 0._r8) then ! melting
2452 dophasechangeflag = .true.
2453 heatavail = (t_lake(c,j) - tfrz) * cv_lake(c,j)
2454 melt = min(lake_icefrac(c,j)*denh2o*dz_lake(c,j), heatavail/hfus)
2455 !denh2o is used because layer thickness is not adjusted for freezing
2456 heatrem = max(heatavail - melt*hfus, 0._r8)
2457 !catch small negative value to keep t at tfrz
2458 else if (t_lake(c,j) < tfrz .and. lake_icefrac(c,j) < 1._r8) then !freezing
2459 dophasechangeflag = .true.
2460 heatavail = (t_lake(c,j) - tfrz) * cv_lake(c,j)
2461 melt = max(-(1._r8-lake_icefrac(c,j))*denh2o*dz_lake(c,j), heatavail/hfus)
2462 !denh2o is used because layer thickness is not adjusted for freezing
2463 heatrem = min(heatavail - melt*hfus, 0._r8)
2464 !catch small positive value to keep t at tfrz
2465 end if
2466 ! Update temperature and ice fraction.
2467 if (dophasechangeflag) then
2468 lake_icefrac(c,j) = lake_icefrac(c,j) - melt/(denh2o*dz_lake(c,j))
2469 lhabs(c) = lhabs(c) + melt*hfus
2470 ! Update heat capacity
2471 cv_lake(c,j) = cv_lake(c,j) + melt*(cpliq-cpice)
2472 t_lake(c,j) = tfrz + heatrem/cv_lake(c,j)
2473 ! Prevent tiny residuals
2474 if (lake_icefrac(c,j) > 1._r8 - smallnumber) lake_icefrac(c,j) = 1._r8
2475 if (lake_icefrac(c,j) < smallnumber) lake_icefrac(c,j) = 0._r8
2476 end if
2477 end do
2478 end do
2480 ! Snow & soil phase change
2482 do j = -nlevsnow+1,nlevsoil
2483 !dir$ concurrent
2484 !cdir nodep
2485 do fc = 1,num_shlakec
2486 c = filter_shlakec(fc)
2487 dophasechangeflag = .false.
2489 if (j >= snl(c) + 1) then
2491 if (t_soisno(c,j) > tfrz .and. h2osoi_ice(c,j) > 0._r8) then ! melting
2492 dophasechangeflag = .true.
2493 heatavail = (t_soisno(c,j) - tfrz) * cv(c,j)
2494 melt = min(h2osoi_ice(c,j), heatavail/hfus)
2495 heatrem = max(heatavail - melt*hfus, 0._r8)
2496 !catch small negative value to keep t at tfrz
2497 if (j <= 0) then !snow
2498 imelt(c,j) = 1
2499 qflx_snomelt(c) = qflx_snomelt(c) + melt
2500 end if
2501 else if (t_soisno(c,j) < tfrz .and. h2osoi_liq(c,j) > 0._r8) then !freezing
2502 dophasechangeflag = .true.
2503 heatavail = (t_soisno(c,j) - tfrz) * cv(c,j)
2504 melt = max(-h2osoi_liq(c,j), heatavail/hfus)
2505 heatrem = min(heatavail - melt*hfus, 0._r8)
2506 !catch small positive value to keep t at tfrz
2507 if (j <= 0) then !snow
2508 imelt(c,j) = 2
2509 qflx_snomelt(c) = qflx_snomelt(c) + melt
2510 ! Does this works for both signs of melt in SnowHydrology? I think
2511 ! qflx_snomelt(c) is just output.
2512 end if
2513 end if
2515 ! Update temperature and soil components.
2516 if (dophasechangeflag) then
2517 h2osoi_ice(c,j) = h2osoi_ice(c,j) - melt
2518 h2osoi_liq(c,j) = h2osoi_liq(c,j) + melt
2519 lhabs(c) = lhabs(c) + melt*hfus
2520 ! Update heat capacity
2521 cv(c,j) = cv(c,j) + melt*(cpliq-cpice)
2522 t_soisno(c,j) = tfrz + heatrem/cv(c,j)
2523 ! Prevent tiny residuals
2524 if (h2osoi_ice(c,j) < smallnumber) h2osoi_ice(c,j) = 0._r8
2525 if (h2osoi_liq(c,j) < smallnumber) h2osoi_liq(c,j) = 0._r8
2526 end if
2528 end if
2529 end do
2530 end do
2532 ! Update eflx_snomelt(c)
2533 !dir$ concurrent
2534 !cdir nodep
2535 do fc = 1,num_shlakec
2536 c = filter_shlakec(fc)
2537 eflx_snomelt(c) = qflx_snomelt(c)*hfus
2538 end do
2539 !!!
2541 end subroutine PhaseChange_Lake
2544 subroutine ShalLakeHydrology(dz_lake,forc_rain,forc_snow, & !i
2545 begwb,qflx_evap_tot,forc_t,do_capsnow, &
2546 t_grnd,qflx_evap_soi, &
2547 qflx_snomelt,imelt,frac_iceold, & !i add by guhp
2548 z,dz,zi,snl,h2osno,snowdp,lake_icefrac,t_lake, & !i&o
2549 endwb,snowage,snowice,snowliq,t_snow, & !o
2550 t_soisno,h2osoi_ice,h2osoi_liq,h2osoi_vol, &
2551 qflx_drain,qflx_surf,qflx_infl,qflx_qrgwl, &
2552 qcharge,qflx_prec_grnd,qflx_snowcap, &
2553 qflx_snowcap_col,qflx_snow_grnd_pft, &
2554 qflx_snow_grnd_col,qflx_rain_grnd, &
2555 qflx_evap_tot_col,soilalpha,zwt,fcov, &
2556 rootr_column,qflx_evap_grnd,qflx_sub_snow, &
2557 qflx_dew_snow,qflx_dew_grnd,qflx_rain_grnd_col)
2559 !==================================================================================
2560 ! !DESCRIPTION:
2561 ! Calculation of Shallow Lake Hydrology. Full hydrology of snow layers is
2562 ! done. However, there is no infiltration, and the water budget is balanced with
2563 ! qflx_qrgwl. Lake water mass is kept constant. The soil is simply maintained at
2564 ! volumetric saturation if ice melting frees up pore space. Likewise, if the water
2565 ! portion alone at some point exceeds pore capacity, it is reduced. This is consistent
2566 ! with the possibility of initializing the soil layer with excess ice. The only
2567 ! real error with that is that the thermal conductivity will ignore the excess ice
2568 ! (and accompanying thickness change).
2569 !
2570 ! If snow layers are present over an unfrozen lake, and the top layer of the lake
2571 ! is capable of absorbing the latent heat without going below freezing,
2572 ! the snow-water is runoff and the latent heat is subtracted from the lake.
2573 !
2574 ! WARNING: This subroutine assumes lake columns have one and only one pft.
2575 !
2576 ! Sequence is:
2577 ! ShalLakeHydrology:
2578 ! Do needed tasks from Hydrology1, Biogeophysics2, & top of Hydrology2.
2579 ! -> SnowWater: change of snow mass and snow water onto soil
2580 ! -> SnowCompaction: compaction of snow layers
2581 ! -> CombineSnowLayers: combine snow layers that are thinner than minimum
2582 ! -> DivideSnowLayers: subdivide snow layers that are thicker than maximum
2583 ! Add water to soil if melting has left it with open pore space.
2584 ! Cleanup and do water balance.
2585 ! If snow layers are found above a lake with unfrozen top layer, whose top
2586 ! layer has enough heat to melt all the snow ice without freezing, do so
2587 ! and eliminate the snow layers.
2588 !
2589 ! !REVISION HISTORY:
2590 ! Created by Zack Subin, 2009
2591 !
2592 !============================================================================================
2594 ! USES:
2595 !
2596 implicit none
2598 ! in:
2600 ! integer , intent(in) :: clandunit(1) ! column's landunit
2601 ! integer , intent(in) :: ityplun(1) ! landunit type
2602 ! real(r8), intent(in) :: watsat(1,1:nlevsoil) ! volumetric soil water at saturation (porosity)
2603 real(r8), intent(in) :: dz_lake(1,nlevlake) ! layer thickness for lake (m)
2604 real(r8), intent(in) :: forc_rain(1) ! rain rate [mm/s]
2605 real(r8), intent(in) :: forc_snow(1) ! snow rate [mm/s]
2606 real(r8), intent(in) :: qflx_evap_tot(1) ! qflx_evap_soi + qflx_evap_veg + qflx_tran_veg
2607 real(r8), intent(in) :: forc_t(1) ! atmospheric temperature (Kelvin)
2608 #if (defined OFFLINE)
2609 real(r8), intent(in) :: flfall(1) ! fraction of liquid water within falling precipitation
2610 #endif
2611 logical , intent(in) :: do_capsnow(1) ! true => do snow capping
2612 real(r8), intent(in) :: t_grnd(1) ! ground temperature (Kelvin)
2613 real(r8), intent(in) :: qflx_evap_soi(1) ! soil evaporation (mm H2O/s) (+ = to atm)
2614 real(r8), intent(in) :: qflx_snomelt(1) !snow melt (mm H2O /s)
2615 integer, intent(in) :: imelt(1,-nlevsnow+1:nlevsoil) !flag for melting (=1), freezing (=2), Not=0
2617 !inout:
2619 real(r8), intent(inout) :: begwb(1) ! water mass begining of the time step
2621 ! inout:
2624 real(r8), intent(inout) :: z(1,-nlevsnow+1:nlevsoil) ! layer depth (m)
2625 real(r8), intent(inout) :: dz(1,-nlevsnow+1:nlevsoil) ! layer thickness depth (m)
2626 real(r8), intent(inout) :: zi(1,-nlevsnow+0:nlevsoil) ! interface depth (m)
2627 integer , intent(inout) :: snl(1) ! number of snow layers
2628 real(r8), intent(inout) :: h2osno(1) ! snow water (mm H2O)
2629 real(r8), intent(inout) :: snowdp(1) ! snow height (m)
2630 real(r8), intent(inout) :: lake_icefrac(1,nlevlake) ! mass fraction of lake layer that is frozen
2631 real(r8), intent(inout) :: t_lake(1,nlevlake) ! lake temperature (Kelvin)
2633 real(r8), intent(inout) :: frac_iceold(1,-nlevsnow+1:nlevsoil) ! fraction of ice relative to the tot water
2634 ! out:
2637 real(r8), intent(out) :: endwb(1) ! water mass end of the time step
2638 real(r8), intent(out) :: snowage(1) ! non dimensional snow age [-]
2639 real(r8), intent(out) :: snowice(1) ! average snow ice lens
2640 real(r8), intent(out) :: snowliq(1) ! average snow liquid water
2641 real(r8), intent(out) :: t_snow(1) ! vertically averaged snow temperature
2642 real(r8), intent(out) :: t_soisno(1,-nlevsnow+1:nlevsoil) ! snow temperature (Kelvin)
2643 real(r8), intent(out) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) ! ice lens (kg/m2)
2644 real(r8), intent(out) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) ! liquid water (kg/m2)
2645 real(r8), intent(out) :: h2osoi_vol(1,-nlevsnow+1:nlevsoil) ! volumetric soil water (0<=h2osoi_vol<=watsat)[m3/m3]
2646 real(r8), intent(out) :: qflx_drain(1) ! sub-surface runoff (mm H2O /s)
2647 real(r8), intent(out) :: qflx_surf(1) ! surface runoff (mm H2O /s)
2648 real(r8), intent(out) :: qflx_infl(1) ! infiltration (mm H2O /s)
2649 real(r8), intent(out) :: qflx_qrgwl(1) ! qflx_surf at glaciers, wetlands, lakes
2650 real(r8), intent(out) :: qcharge(1) ! aquifer recharge rate (mm/s)
2651 real(r8), intent(out) :: qflx_prec_grnd(1) ! water onto ground including canopy runoff [kg/(m2 s)]
2652 real(r8), intent(out) :: qflx_snowcap(1) ! excess precipitation due to snow capping (mm H2O /s) [+]
2653 real(r8), intent(out) :: qflx_snowcap_col(1) ! excess precipitation due to snow capping (mm H2O /s) [+]
2654 real(r8), intent(out) :: qflx_snow_grnd_pft(1) ! snow on ground after interception (mm H2O/s) [+]
2655 real(r8), intent(out) :: qflx_snow_grnd_col(1) ! snow on ground after interception (mm H2O/s) [+]
2656 real(r8), intent(out) :: qflx_rain_grnd(1) ! rain on ground after interception (mm H2O/s) [+]
2657 real(r8), intent(out) :: qflx_evap_tot_col(1) !pft quantity averaged to the column (assuming one pft)
2658 real(r8) ,intent(out) :: soilalpha(1) !factor that reduces ground saturated specific humidity (-)
2659 real(r8), intent(out) :: zwt(1) !water table depth
2660 real(r8), intent(out) :: fcov(1) !fractional area with water table at surface
2661 real(r8), intent(out) :: rootr_column(1,1:nlevsoil) !effective fraction of roots in each soil layer
2662 real(r8), intent(out) :: qflx_evap_grnd(1) ! ground surface evaporation rate (mm H2O/s) [+]
2663 real(r8), intent(out) :: qflx_sub_snow(1) ! sublimation rate from snow pack (mm H2O /s) [+]
2664 real(r8), intent(out) :: qflx_dew_snow(1) ! surface dew added to snow pack (mm H2O /s) [+]
2665 real(r8), intent(out) :: qflx_dew_grnd(1) ! ground surface dew formation (mm H2O /s) [+]
2666 real(r8), intent(out) :: qflx_rain_grnd_col(1) !rain on ground after interception (mm H2O/s) [+]
2668 ! Block of biogeochem currently not used.
2669 #ifndef SHLAKE
2670 real(r8), pointer :: sucsat(:,:) ! minimum soil suction (mm)
2671 real(r8), pointer :: bsw(:,:) ! Clapp and Hornberger "b"
2672 real(r8), pointer :: bsw2(:,:) ! Clapp and Hornberger "b" for CN code
2673 real(r8), pointer :: psisat(:,:) ! soil water potential at saturation for CN code (MPa)
2674 real(r8), pointer :: vwcsat(:,:) ! volumetric water content at saturation for CN code (m3/m3)
2675 real(r8), pointer :: wf(:) ! soil water as frac. of whc for top 0.5 m
2676 real(r8), pointer :: soilpsi(:,:) ! soil water potential in each soil layer (MPa)
2677 real(r8) :: psi,vwc,fsat ! temporary variables for soilpsi calculation
2678 #if (defined DGVM) || (defined CN)
2679 real(r8) :: watdry ! temporary
2680 real(r8) :: rwat(lbc:ubc) ! soil water wgted by depth to maximum depth of 0.5 m
2681 real(r8) :: swat(lbc:ubc) ! same as rwat but at saturation
2682 real(r8) :: rz(lbc:ubc) ! thickness of soil layers contributing to rwat (m)
2683 real(r8) :: tsw ! volumetric soil water to 0.5 m
2684 real(r8) :: stsw ! volumetric soil water to 0.5 m at saturation
2685 #endif
2686 #endif
2689 ! OTHER LOCAL VARIABLES:
2691 integer :: p,fp,g,l,c,j,fc,jtop ! indices
2692 integer :: num_shlakesnowc ! number of column snow points
2693 integer :: filter_shlakesnowc(ubc-lbc+1) ! column filter for snow points
2694 integer :: num_shlakenosnowc ! number of column non-snow points
2695 integer :: filter_shlakenosnowc(ubc-lbc+1) ! column filter for non-snow points
2696 ! real(r8) :: dtime ! land model time step (sec)
2697 integer :: newnode ! flag when new snow node is set, (1=yes, 0=no)
2698 real(r8) :: dz_snowf ! layer thickness rate change due to precipitation [mm/s]
2699 real(r8) :: bifall ! bulk density of newly fallen dry snow [kg/m3]
2700 real(r8) :: fracsnow(lbp:ubp) ! frac of precipitation that is snow
2701 real(r8) :: fracrain(lbp:ubp) ! frac of precipitation that is rain
2702 real(r8) :: qflx_prec_grnd_snow(lbp:ubp) ! snow precipitation incident on ground [mm/s]
2703 real(r8) :: qflx_prec_grnd_rain(lbp:ubp) ! rain precipitation incident on ground [mm/s]
2704 real(r8) :: qflx_evap_soi_lim ! temporary evap_soi limited by top snow layer content [mm/s]
2705 real(r8) :: h2osno_temp ! temporary h2osno [kg/m^2]
2706 real(r8), parameter :: snow_bd = 250._r8 !constant snow bulk density (only used in special case here) [kg/m^3]
2707 real(r8) :: sumsnowice(lbc:ubc) ! sum of snow ice if snow layers found above unfrozen lake [kg/m&2]
2708 logical :: unfrozen(lbc:ubc) ! true if top lake layer is unfrozen with snow layers above
2709 real(r8) :: heatrem ! used in case above [J/m^2]
2710 real(r8) :: heatsum(lbc:ubc) ! used in case above [J/m^2]
2711 real(r8) :: qflx_top_soil(1) !net water input into soil from top (mm/s)
2712 character*256 :: message
2714 #if (defined LAKEDEBUG)
2715 real(r8) :: snow_water(lbc:ubc) ! temporary sum of snow water for Bal Check [kg/m^2]
2716 #endif
2717 !-----------------------------------------------------------------------
2720 ! Determine step size
2722 ! dtime = get_step_size()
2724 ! Add soil water to water balance.
2725 do j = 1, nlevsoil
2726 !dir$ concurrent
2727 !cdir nodep
2728 do fc = 1, num_shlakec
2729 c = filter_shlakec(fc)
2730 begwb(c) = begwb(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
2731 end do
2732 end do
2734 !!!!!!!!!!!!!!!!!!!!!!!!!!!
2736 ! Do precipitation onto ground, etc., from Hydrology1.
2738 !dir$ concurrent
2739 !cdir nodep
2740 do fp = 1, num_shlakep
2741 p = filter_shlakep(fp)
2742 g = pgridcell(p)
2743 ! l = plandunit(p)
2744 c = pcolumn(p)
2746 ! Precipitation onto ground (kg/(m2 s))
2747 ! ! PET, 1/18/2005: Added new terms for mass balance correction
2748 ! ! due to dynamic pft weight shifting (column-level h2ocan_loss)
2749 ! ! Because the fractionation between rain and snow is indeterminate if
2750 ! ! rain + snow = 0, I am adding this very small flux only to the rain
2751 ! ! components.
2752 ! Not relevant unless PFTs are added to lake later.
2753 ! if (frac_veg_nosno(p) == 0) then
2754 qflx_prec_grnd_snow(p) = forc_snow(g)
2755 qflx_prec_grnd_rain(p) = forc_rain(g) !+ h2ocan_loss(c)
2756 ! else
2757 ! qflx_prec_grnd_snow(p) = qflx_through_snow(p) + (qflx_candrip(p) * fracsnow(p))
2758 ! qflx_prec_grnd_rain(p) = qflx_through_rain(p) + (qflx_candrip(p) * fracrain(p)) + h2ocan_loss(c)
2759 ! end if
2760 qflx_prec_grnd(p) = qflx_prec_grnd_snow(p) + qflx_prec_grnd_rain(p)
2762 if (do_capsnow(c)) then
2763 qflx_snowcap(p) = qflx_prec_grnd_snow(p) + qflx_prec_grnd_rain(p)
2764 qflx_snow_grnd_pft(p) = 0._r8
2765 qflx_rain_grnd(p) = 0._r8
2766 else
2767 qflx_snowcap(p) = 0._r8
2768 #if (defined OFFLINE)
2769 qflx_snow_grnd_pft(p) = qflx_prec_grnd(p)*(1._r8-flfall(g)) ! ice onto ground (mm/s)
2770 qflx_rain_grnd(p) = qflx_prec_grnd(p)*flfall(g) ! liquid water onto ground (mm/s)
2771 #else
2772 qflx_snow_grnd_pft(p) = qflx_prec_grnd_snow(p) ! ice onto ground (mm/s)
2773 qflx_rain_grnd(p) = qflx_prec_grnd_rain(p) ! liquid water onto ground (mm/s)
2774 #endif
2775 end if
2776 ! Assuming one PFT; needed for below
2777 qflx_snow_grnd_col(c) = qflx_snow_grnd_pft(p)
2778 qflx_rain_grnd_col(c) = qflx_rain_grnd(p)
2780 end do ! (end pft loop)
2782 ! Determine snow height and snow water
2784 !dir$ concurrent
2785 !cdir nodep
2786 do fc = 1, num_shlakec
2787 c = filter_shlakec(fc)
2788 ! l = clandunit(c)
2789 g = cgridcell(c)
2791 ! Use Alta relationship, Anderson(1976); LaChapelle(1961),
2792 ! U.S.Department of Agriculture Forest Service, Project F,
2793 ! Progress Rep. 1, Alta Avalanche Study Center:Snow Layer Densification.
2795 if (do_capsnow(c)) then
2796 dz_snowf = 0._r8
2797 else
2798 if (forc_t(g) > tfrz + 2._r8) then
2799 bifall=50._r8 + 1.7_r8*(17.0_r8)**1.5_r8
2800 else if (forc_t(g) > tfrz - 15._r8) then
2801 bifall=50._r8 + 1.7_r8*(forc_t(g) - tfrz + 15._r8)**1.5_r8
2802 else
2803 bifall=50._r8
2804 end if
2805 dz_snowf = qflx_snow_grnd_col(c)/bifall
2806 snowdp(c) = snowdp(c) + dz_snowf*dtime
2807 h2osno(c) = h2osno(c) + qflx_snow_grnd_col(c)*dtime ! snow water equivalent (mm)
2808 end if
2810 ! if (itype(l)==istwet .and. t_grnd(c)>tfrz) then
2811 ! h2osno(c)=0._r8
2812 ! snowdp(c)=0._r8
2813 ! snowage(c)=0._r8
2814 ! end if
2815 ! Take care of this later in function.
2817 ! When the snow accumulation exceeds 10 mm, initialize snow layer
2818 ! Currently, the water temperature for the precipitation is simply set
2819 ! as the surface air temperature
2821 newnode = 0 ! flag for when snow node will be initialized
2822 if (snl(c) == 0 .and. qflx_snow_grnd_col(c) > 0.0_r8 .and. snowdp(c) >= 0.01_r8) then
2823 newnode = 1
2824 snl(c) = -1
2825 dz(c,0) = snowdp(c) ! meter
2826 z(c,0) = -0.5_r8*dz(c,0)
2827 zi(c,-1) = -dz(c,0)
2828 snowage(c) = 0._r8 ! snow age
2829 t_soisno(c,0) = min(tfrz, forc_t(g)) ! K
2830 h2osoi_ice(c,0) = h2osno(c) ! kg/m2
2831 h2osoi_liq(c,0) = 0._r8 ! kg/m2
2832 frac_iceold(c,0) = 1._r8
2833 end if
2835 ! The change of ice partial density of surface node due to precipitation.
2836 ! Only ice part of snowfall is added here, the liquid part will be added
2837 ! later.
2839 if (snl(c) < 0 .and. newnode == 0) then
2840 h2osoi_ice(c,snl(c)+1) = h2osoi_ice(c,snl(c)+1)+dtime*qflx_snow_grnd_col(c)
2841 dz(c,snl(c)+1) = dz(c,snl(c)+1)+dz_snowf*dtime
2842 end if
2844 end do
2846 ! Calculate sublimation and dew, adapted from HydrologyLake and Biogeophysics2.
2848 !dir$ concurrent
2849 !cdir nodep
2850 do fp = 1,num_shlakep
2851 p = filter_shlakep(fp)
2852 c = pcolumn(p)
2853 jtop = snl(c)+1
2855 ! Use column variables here
2856 qflx_evap_grnd(c) = 0._r8
2857 qflx_sub_snow(c) = 0._r8
2858 qflx_dew_snow(c) = 0._r8
2859 qflx_dew_grnd(c) = 0._r8
2861 if (jtop <= 0) then ! snow layers
2862 j = jtop
2863 ! Assign ground evaporation to sublimation from soil ice or to dew
2864 ! on snow or ground
2866 if (qflx_evap_soi(p) >= 0._r8) then
2867 ! for evaporation partitioning between liquid evap and ice sublimation,
2868 ! use the ratio of liquid to (liquid+ice) in the top layer to determine split
2869 ! Since we're not limiting evap over lakes, but still can't remove more from top
2870 ! snow layer than there is there, create temp. limited evap_soi.
2871 qflx_evap_soi_lim = min(qflx_evap_soi(p), (h2osoi_liq(c,j)+h2osoi_ice(c,j))/dtime)
2872 if ((h2osoi_liq(c,j)+h2osoi_ice(c,j)) > 0._r8) then
2873 qflx_evap_grnd(c) = max(qflx_evap_soi_lim*(h2osoi_liq(c,j)/(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._r8)
2874 else
2875 qflx_evap_grnd(c) = 0._r8
2876 end if
2877 qflx_sub_snow(c) = qflx_evap_soi_lim - qflx_evap_grnd(c)
2878 else
2879 if (t_grnd(c) < tfrz) then
2880 qflx_dew_snow(c) = abs(qflx_evap_soi(p))
2881 else
2882 qflx_dew_grnd(c) = abs(qflx_evap_soi(p))
2883 end if
2884 end if
2885 ! Update the pft-level qflx_snowcap
2886 ! This was moved in from Hydrology2 to keep all pft-level
2887 ! calculations out of Hydrology2
2888 if (do_capsnow(c)) qflx_snowcap(p) = qflx_snowcap(p) + qflx_dew_snow(c) + qflx_dew_grnd(c)
2890 else ! No snow layers: do as in HydrologyLake but with actual clmtype variables
2891 if (qflx_evap_soi(p) >= 0._r8) then
2892 ! Sublimation: do not allow for more sublimation than there is snow
2893 ! after melt. Remaining surface evaporation used for infiltration.
2894 qflx_sub_snow(c) = min(qflx_evap_soi(p), h2osno(c)/dtime)
2895 qflx_evap_grnd(c) = qflx_evap_soi(p) - qflx_sub_snow(c)
2896 else
2897 if (t_grnd(c) < tfrz-0.1_r8) then
2898 qflx_dew_snow(c) = abs(qflx_evap_soi(p))
2899 else
2900 qflx_dew_grnd(c) = abs(qflx_evap_soi(p))
2901 end if
2902 end if
2904 ! Update snow pack for dew & sub.
2905 h2osno_temp = h2osno(c)
2906 if (do_capsnow(c)) then
2907 h2osno(c) = h2osno(c) - qflx_sub_snow(c)*dtime
2908 qflx_snowcap(p) = qflx_snowcap(p) + qflx_dew_snow(c) + qflx_dew_grnd(c)
2909 else
2910 h2osno(c) = h2osno(c) + (-qflx_sub_snow(c)+qflx_dew_snow(c))*dtime
2911 end if
2912 if (h2osno_temp > 0._r8) then
2913 snowdp(c) = snowdp(c) * h2osno(c) / h2osno_temp
2914 else
2915 snowdp(c) = h2osno(c)/snow_bd !Assume a constant snow bulk density = 250.
2916 end if
2918 #if (defined PERGRO)
2919 if (abs(h2osno(c)) < 1.e-10_r8) h2osno(c) = 0._r8
2920 #else
2921 h2osno(c) = max(h2osno(c), 0._r8)
2922 #endif
2924 end if
2926 qflx_snowcap_col(c) = qflx_snowcap(p)
2928 end do
2931 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2932 ! Determine initial snow/no-snow filters (will be modified possibly by
2933 ! routines CombineSnowLayers and DivideSnowLayers below
2935 call BuildSnowFilter(lbc, ubc, num_shlakec, filter_shlakec,snl, & !i
2936 num_shlakesnowc, filter_shlakesnowc, num_shlakenosnowc, filter_shlakenosnowc) !o
2938 ! Determine the change of snow mass and the snow water onto soil
2940 call SnowWater(lbc, ubc, num_shlakesnowc, filter_shlakesnowc, & !i
2941 num_shlakenosnowc, filter_shlakenosnowc, & !i
2942 snl,do_capsnow,qflx_snomelt,qflx_rain_grnd, & !i
2943 qflx_sub_snow,qflx_evap_grnd, & !i
2944 qflx_dew_snow,qflx_dew_grnd,dz, & !i
2945 h2osoi_ice,h2osoi_liq, & !i&o
2946 qflx_top_soil) !o
2949 ! Determine soil hydrology
2950 ! Here this consists only of making sure that soil is saturated even as it melts and 10%
2951 ! of pore space opens up. Conversely, if excess ice is melting and the liquid water exceeds the
2952 ! saturation value, then remove water.
2954 do j = 1,nlevsoil
2955 !dir$ concurrent
2956 !cdir nodep
2957 do fc = 1, num_shlakec
2958 c = filter_shlakec(fc)
2960 if (h2osoi_vol(c,j) < watsat(c,j)) then
2961 h2osoi_liq(c,j) = (watsat(c,j)*dz(c,j) - h2osoi_ice(c,j)/denice)*denh2o
2962 ! h2osoi_vol will be updated below, and this water addition will come from qflx_qrgwl
2963 else if (h2osoi_liq(c,j) > watsat(c,j)*denh2o*dz(c,j)) then
2964 h2osoi_liq(c,j) = watsat(c,j)*denh2o*dz(c,j)
2965 end if
2967 end do
2968 end do
2969 !!!!!!!!!!
2971 ! if (.not. is_perpetual()) then
2972 if (1==1) then
2974 ! Natural compaction and metamorphosis.
2976 call SnowCompaction(lbc, ubc, num_shlakesnowc, filter_shlakesnowc, &!i
2977 snl,imelt,frac_iceold,t_soisno, &!i
2978 h2osoi_ice,h2osoi_liq, &!i
2979 dz) !&o
2981 ! Combine thin snow elements
2983 call CombineSnowLayers(lbc, ubc, & !i
2984 num_shlakesnowc, filter_shlakesnowc, & !i&o
2985 snl,h2osno,snowdp,dz,zi, & !i&o
2986 t_soisno,h2osoi_ice,h2osoi_liq, & !i&o
2987 z) !o
2990 ! Divide thick snow elements
2992 call DivideSnowLayers(lbc, ubc, & !i
2993 num_shlakesnowc, filter_shlakesnowc, & !i&o
2994 snl,dz,zi,t_soisno, & !i&o
2995 h2osoi_ice,h2osoi_liq, & !i&o
2996 z) !o
2999 else
3001 do fc = 1, num_shlakesnowc
3002 c = filter_shlakesnowc(fc)
3003 h2osno(c) = 0._r8
3004 end do
3005 do j = -nlevsnow+1,0
3006 do fc = 1, num_shlakesnowc
3007 c = filter_shlakesnowc(fc)
3008 if (j >= snl(c)+1) then
3009 h2osno(c) = h2osno(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
3010 end if
3011 end do
3012 end do
3014 end if
3016 ! Check for snow layers above lake with unfrozen top layer. Mechanically,
3017 ! the snow will fall into the lake and melt or turn to ice. If the top layer has
3018 ! sufficient heat to melt the snow without freezing, then that will be done.
3019 ! Otherwise, the top layer will undergo freezing, but only if the top layer will
3020 ! not freeze completely. Otherwise, let the snow layers persist and melt by diffusion.
3021 !dir$ concurrent
3022 !cdir nodep
3023 do fc = 1, num_shlakec
3024 c = filter_shlakec(fc)
3026 if (t_lake(c,1) > tfrz .and. lake_icefrac(c,1) == 0._r8 .and. snl(c) < 0) then
3027 unfrozen(c) = .true.
3028 else
3029 unfrozen(c) = .false.
3030 end if
3031 end do
3033 do j = -nlevsnow+1,0
3034 !dir$ concurrent
3035 !cdir nodep
3036 do fc = 1, num_shlakec
3037 c = filter_shlakec(fc)
3039 if (unfrozen(c)) then
3040 if (j == -nlevsnow+1) then
3041 sumsnowice(c) = 0._r8
3042 heatsum(c) = 0._r8
3043 end if
3044 if (j >= snl(c)+1) then
3045 sumsnowice(c) = sumsnowice(c) + h2osoi_ice(c,j)
3046 heatsum(c) = heatsum(c) + h2osoi_ice(c,j)*cpice*(tfrz - t_soisno(c,j)) &
3047 + h2osoi_liq(c,j)*cpliq*(tfrz - t_soisno(c,j))
3048 end if
3049 end if
3050 end do
3051 end do
3053 !dir$ concurrent
3054 !cdir nodep
3055 do fc = 1, num_shlakec
3056 c = filter_shlakec(fc)
3058 if (unfrozen(c)) then
3059 heatsum(c) = heatsum(c) + sumsnowice(c)*hfus
3060 heatrem = (t_lake(c,1) - tfrz)*cpliq*denh2o*dz_lake(c,1) - heatsum(c)
3062 if (heatrem + denh2o*dz_lake(c,1)*hfus > 0._r8) then
3063 ! Remove snow and subtract the latent heat from the top layer.
3064 h2osno(c) = 0._r8
3065 snl(c) = 0
3066 ! The rest of the bookkeeping for the removed snow will be done below.
3067 #if (defined LAKEDEBUG)
3068 write(message,*)'Snow layers removed above unfrozen lake for column, snowice:', &
3069 c, sumsnowice(c)
3070 CALL wrf_message(message)
3071 #endif
3072 if (heatrem > 0._r8) then ! simply subtract the heat from the layer
3073 t_lake(c,1) = t_lake(c,1) - heatrem/(cpliq*denh2o*dz_lake(c,1))
3074 else !freeze part of the layer
3075 t_lake(c,1) = tfrz
3076 lake_icefrac(c,1) = -heatrem/(denh2o*dz_lake(c,1)*hfus)
3077 end if
3078 end if
3079 end if
3080 end do
3081 !!!!!!!!!!!!
3083 ! Set snow age to zero if no snow
3085 !dir$ concurrent
3086 !cdir nodep
3087 do fc = 1, num_shlakesnowc
3088 c = filter_shlakesnowc(fc)
3089 if (snl(c) == 0) then
3090 snowage(c) = 0._r8
3091 end if
3092 end do
3094 ! Set empty snow layers to zero
3096 do j = -nlevsnow+1,0
3097 !dir$ concurrent
3098 !cdir nodep
3099 do fc = 1, num_shlakesnowc
3100 c = filter_shlakesnowc(fc)
3101 if (j <= snl(c) .and. snl(c) > -nlevsnow) then
3102 h2osoi_ice(c,j) = 0._r8
3103 h2osoi_liq(c,j) = 0._r8
3104 t_soisno(c,j) = 0._r8
3105 dz(c,j) = 0._r8
3106 z(c,j) = 0._r8
3107 zi(c,j-1) = 0._r8
3108 end if
3109 end do
3110 end do
3112 ! Build new snow filter
3114 call BuildSnowFilter(lbc, ubc, num_shlakec, filter_shlakec, snl,& !i
3115 num_shlakesnowc, filter_shlakesnowc, num_shlakenosnowc, filter_shlakenosnowc) !o
3117 ! Vertically average t_soisno and sum of h2osoi_liq and h2osoi_ice
3118 ! over all snow layers for history output
3120 !dir$ concurrent
3121 !cdir nodep
3122 do fc = 1, num_shlakesnowc
3123 c = filter_shlakesnowc(fc)
3124 t_snow(c) = 0._r8
3125 snowice(c) = 0._r8
3126 snowliq(c) = 0._r8
3127 end do
3128 !dir$ concurrent
3129 !cdir nodep
3130 do fc = 1, num_shlakenosnowc
3131 c = filter_shlakenosnowc(fc)
3132 t_snow(c) = spval
3133 snowice(c) = spval
3134 snowliq(c) = spval
3135 end do
3137 do j = -nlevsnow+1, 0
3138 !dir$ concurrent
3139 !cdir nodep
3140 do fc = 1, num_shlakesnowc
3141 c = filter_shlakesnowc(fc)
3142 if (j >= snl(c)+1) then
3143 t_snow(c) = t_snow(c) + t_soisno(c,j)
3144 snowice(c) = snowice(c) + h2osoi_ice(c,j)
3145 snowliq(c) = snowliq(c) + h2osoi_liq(c,j)
3146 end if
3147 end do
3148 end do
3150 ! Determine ending water balance and volumetric soil water
3152 !dir$ concurrent
3153 !cdir nodep
3154 do fc = 1, num_shlakec
3156 c = filter_shlakec(fc)
3157 if (snl(c) < 0) t_snow(c) = t_snow(c)/abs(snl(c))
3158 endwb(c) = h2osno(c)
3159 end do
3161 do j = 1, nlevsoil
3162 !dir$ concurrent
3163 !cdir nodep
3164 do fc = 1, num_shlakec
3165 c = filter_shlakec(fc)
3166 endwb(c) = endwb(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
3167 h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice)
3168 end do
3169 end do
3171 #if (defined LAKEDEBUG)
3172 ! Check to make sure snow water adds up correctly.
3173 do j = -nlevsnow+1,0
3174 !dir$ concurrent
3175 !cdir nodep
3176 do fc = 1, num_shlakec
3177 c = filter_shlakec(fc)
3179 jtop = snl(c)+1
3180 if(j == jtop) snow_water(c) = 0._r8
3181 if(j >= jtop) then
3182 snow_water(c) = snow_water(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
3183 if(j == 0 .and. abs(snow_water(c)-h2osno(c))>1.e-7_r8) then
3184 write(message,*)'h2osno does not equal sum of snow layers in ShalLakeHydrology:', &
3185 'column, h2osno, sum of snow layers =', c, h2osno(c), snow_water(c)
3186 CALL wrf_error_fatal ( message )
3187 end if
3188 end if
3189 end do
3190 end do
3191 #endif
3193 !!!!!!!!!!!!!
3194 ! Do history variables and set special landunit runoff (adapted from end of HydrologyLake)
3195 !dir$ concurrent
3196 !cdir nodep
3197 do fp = 1,num_shlakep
3198 p = filter_shlakep(fp)
3199 c = pcolumn(p)
3200 g = pgridcell(p)
3202 qflx_infl(c) = 0._r8
3203 qflx_surf(c) = 0._r8
3204 qflx_drain(c) = 0._r8
3205 rootr_column(c,:) = spval
3206 soilalpha(c) = spval
3207 zwt(c) = spval
3208 fcov(c) = spval
3209 qcharge(c) = spval
3210 ! h2osoi_vol(c,:) = spval
3212 ! Insure water balance using qflx_qrgwl
3213 qflx_qrgwl(c) = forc_rain(g) + forc_snow(g) - qflx_evap_tot(p) - (endwb(c)-begwb(c))/dtime
3214 #if (defined LAKEDEBUG)
3215 write(message,*)'c, rain, snow, evap, endwb, begwb, qflx_qrgwl:', &
3216 c, forc_rain(g), forc_snow(g), qflx_evap_tot(p), endwb(c), begwb(c), qflx_qrgwl(c)
3217 CALL wrf_message(message)
3218 #endif
3220 ! The pft average must be done here for output to history tape
3221 qflx_evap_tot_col(c) = qflx_evap_tot(p)
3222 end do
3224 !!!!!!!!!!!!!
3225 !For now, bracket off the remaining biogeochem code. May need to bring it back
3226 !to do soil carbon and methane beneath lakes.
3227 #if (defined CN)
3228 #ifndef SHLAKE
3229 do j = 1, nlevsoil
3230 !dir$ concurrent
3231 !cdir nodep
3232 do fc = 1, num_soilc
3233 c = filter_soilc(fc)
3235 if (h2osoi_liq(c,j) > 0._r8) then
3236 vwc = h2osoi_liq(c,j)/(dz(c,j)*denh2o)
3238 ! the following limit set to catch very small values of
3239 ! fractional saturation that can crash the calculation of psi
3241 fsat = max(vwc/vwcsat(c,j), 0.001_r8)
3242 psi = psisat(c,j) * (fsat)**bsw2(c,j)
3243 soilpsi(c,j) = min(max(psi,-15.0_r8),0._r8)
3244 else
3245 soilpsi(c,j) = -15.0_r8
3246 end if
3247 end do
3248 end do
3249 #endif
3250 #endif
3252 #if (defined DGVM) || (defined CN)
3253 #ifndef SHLAKE
3254 ! Available soil water up to a depth of 0.5 m.
3255 ! Potentially available soil water (=whc) up to a depth of 0.5 m.
3256 ! Water content as fraction of whc up to a depth of 0.5 m.
3258 !dir$ concurrent
3259 !cdir nodep
3260 do c = lbc,ubc
3261 l = clandunit(c)
3262 if (ityplun(l) == istsoil) then
3263 rwat(c) = 0._r8
3264 swat(c) = 0._r8
3265 rz(c) = 0._r8
3266 end if
3267 end do
3269 do j = 1, nlevsoil
3270 !dir$ concurrent
3271 !cdir nodep
3272 do c = lbc,ubc
3273 l = clandunit(c)
3274 if (ityplun(l) == istsoil) then
3275 if (z(c,j)+0.5_r8*dz(c,j) <= 0.5_r8) then
3276 watdry = watsat(c,j) * (316230._r8/sucsat(c,j)) ** (-1._r8/bsw(c,j))
3277 rwat(c) = rwat(c) + (h2osoi_vol(c,j)-watdry) * dz(c,j)
3278 swat(c) = swat(c) + (watsat(c,j) -watdry) * dz(c,j)
3279 rz(c) = rz(c) + dz(c,j)
3280 end if
3281 end if
3282 end do
3283 end do
3285 !dir$ concurrent
3286 !cdir nodep
3287 do c = lbc,ubc
3288 l = clandunit(c)
3289 if (ityplun(l) == istsoil) then
3290 if (rz(c) /= 0._r8) then
3291 tsw = rwat(c)/rz(c)
3292 stsw = swat(c)/rz(c)
3293 else
3294 watdry = watsat(c,1) * (316230._r8/sucsat(c,1)) ** (-1._r8/bsw(c,1))
3295 tsw = h2osoi_vol(c,1) - watdry
3296 stsw = watsat(c,1) - watdry
3297 end if
3298 wf(c) = tsw/stsw
3299 else
3300 wf(c) = 1.0_r8
3301 end if
3302 end do
3304 #endif
3305 #endif
3307 end subroutine ShalLakeHydrology
3309 subroutine QSat (T, p, es, esdT, qs, qsdT)
3310 !
3311 ! !DESCRIPTION:
3312 ! Computes saturation mixing ratio and the change in saturation
3313 ! mixing ratio with respect to temperature.
3314 ! Reference: Polynomial approximations from:
3315 ! Piotr J. Flatau, et al.,1992: Polynomial fits to saturation
3316 ! vapor pressure. Journal of Applied Meteorology, 31, 1507-1513.
3317 !
3318 ! !USES:
3319 !
3320 ! !ARGUMENTS:
3321 implicit none
3322 real(r8), intent(in) :: T ! temperature (K)
3323 real(r8), intent(in) :: p ! surface atmospheric pressure (pa)
3324 real(r8), intent(out) :: es ! vapor pressure (pa)
3325 real(r8), intent(out) :: esdT ! d(es)/d(T)
3326 real(r8), intent(out) :: qs ! humidity (kg/kg)
3327 real(r8), intent(out) :: qsdT ! d(qs)/d(T)
3328 !
3329 ! !CALLED FROM:
3330 ! subroutine Biogeophysics1 in module Biogeophysics1Mod
3331 ! subroutine BiogeophysicsLake in module BiogeophysicsLakeMod
3332 ! subroutine CanopyFluxesMod CanopyFluxesMod
3333 !
3334 ! !REVISION HISTORY:
3335 ! 15 September 1999: Yongjiu Dai; Initial code
3336 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
3337 !
3338 !EOP
3339 !
3340 ! !LOCAL VARIABLES:
3341 !
3342 real(r8) :: T_limit
3343 real(r8) :: td,vp,vp1,vp2
3344 !
3345 ! For water vapor (temperature range 0C-100C)
3346 !
3347 real(r8), parameter :: a0 = 6.11213476
3348 real(r8), parameter :: a1 = 0.444007856
3349 real(r8), parameter :: a2 = 0.143064234e-01
3350 real(r8), parameter :: a3 = 0.264461437e-03
3351 real(r8), parameter :: a4 = 0.305903558e-05
3352 real(r8), parameter :: a5 = 0.196237241e-07
3353 real(r8), parameter :: a6 = 0.892344772e-10
3354 real(r8), parameter :: a7 = -0.373208410e-12
3355 real(r8), parameter :: a8 = 0.209339997e-15
3356 !
3357 ! For derivative:water vapor
3358 !
3359 real(r8), parameter :: b0 = 0.444017302
3360 real(r8), parameter :: b1 = 0.286064092e-01
3361 real(r8), parameter :: b2 = 0.794683137e-03
3362 real(r8), parameter :: b3 = 0.121211669e-04
3363 real(r8), parameter :: b4 = 0.103354611e-06
3364 real(r8), parameter :: b5 = 0.404125005e-09
3365 real(r8), parameter :: b6 = -0.788037859e-12
3366 real(r8), parameter :: b7 = -0.114596802e-13
3367 real(r8), parameter :: b8 = 0.381294516e-16
3368 !
3369 ! For ice (temperature range -75C-0C)
3370 !
3371 real(r8), parameter :: c0 = 6.11123516
3372 real(r8), parameter :: c1 = 0.503109514
3373 real(r8), parameter :: c2 = 0.188369801e-01
3374 real(r8), parameter :: c3 = 0.420547422e-03
3375 real(r8), parameter :: c4 = 0.614396778e-05
3376 real(r8), parameter :: c5 = 0.602780717e-07
3377 real(r8), parameter :: c6 = 0.387940929e-09
3378 real(r8), parameter :: c7 = 0.149436277e-11
3379 real(r8), parameter :: c8 = 0.262655803e-14
3380 !
3381 ! For derivative:ice
3382 !
3383 real(r8), parameter :: d0 = 0.503277922
3384 real(r8), parameter :: d1 = 0.377289173e-01
3385 real(r8), parameter :: d2 = 0.126801703e-02
3386 real(r8), parameter :: d3 = 0.249468427e-04
3387 real(r8), parameter :: d4 = 0.313703411e-06
3388 real(r8), parameter :: d5 = 0.257180651e-08
3389 real(r8), parameter :: d6 = 0.133268878e-10
3390 real(r8), parameter :: d7 = 0.394116744e-13
3391 real(r8), parameter :: d8 = 0.498070196e-16
3392 !-----------------------------------------------------------------------
3394 T_limit = T - tfrz
3395 if (T_limit > 100.0) T_limit=100.0
3396 if (T_limit < -75.0) T_limit=-75.0
3398 td = T_limit
3399 if (td >= 0.0) then
3400 es = a0 + td*(a1 + td*(a2 + td*(a3 + td*(a4 &
3401 + td*(a5 + td*(a6 + td*(a7 + td*a8)))))))
3402 esdT = b0 + td*(b1 + td*(b2 + td*(b3 + td*(b4 &
3403 + td*(b5 + td*(b6 + td*(b7 + td*b8)))))))
3404 else
3405 es = c0 + td*(c1 + td*(c2 + td*(c3 + td*(c4 &
3406 + td*(c5 + td*(c6 + td*(c7 + td*c8)))))))
3407 esdT = d0 + td*(d1 + td*(d2 + td*(d3 + td*(d4 &
3408 + td*(d5 + td*(d6 + td*(d7 + td*d8)))))))
3409 endif
3411 es = es * 100. ! pa
3412 esdT = esdT * 100. ! pa/K
3414 vp = 1.0 / (p - 0.378*es)
3415 vp1 = 0.622 * vp
3416 vp2 = vp1 * vp
3418 qs = es * vp1 ! kg/kg
3419 qsdT = esdT * vp2 * p ! 1 / K
3421 end subroutine QSat
3424 subroutine Tridiagonal (lbc, ubc, lbj, ubj, jtop, numf, filter, &
3425 a, b, c, r, u)
3426 !
3427 ! !DESCRIPTION:
3428 ! Tridiagonal matrix solution
3429 !
3430 ! !USES:
3431 ! use shr_kind_mod, only: r8 => shr_kind_r8
3432 !
3433 ! !ARGUMENTS:
3434 implicit none
3435 integer , intent(in) :: lbc, ubc ! lbinning and ubing column indices
3436 integer , intent(in) :: lbj, ubj ! lbinning and ubing level indices
3437 integer , intent(in) :: jtop(lbc:ubc) ! top level for each column
3438 integer , intent(in) :: numf ! filter dimension
3439 integer , intent(in) :: filter(1:numf) ! filter
3440 real(r8), intent(in) :: a(lbc:ubc, lbj:ubj) ! "a" left off diagonal of tridiagonal matrix
3441 real(r8), intent(in) :: b(lbc:ubc, lbj:ubj) ! "b" diagonal column for tridiagonal matrix
3442 real(r8), intent(in) :: c(lbc:ubc, lbj:ubj) ! "c" right off diagonal tridiagonal matrix
3443 real(r8), intent(in) :: r(lbc:ubc, lbj:ubj) ! "r" forcing term of tridiagonal matrix
3444 real(r8), intent(inout) :: u(lbc:ubc, lbj:ubj) ! solution
3445 !
3446 ! !CALLED FROM:
3447 ! subroutine BiogeophysicsLake in module BiogeophysicsLakeMod
3448 ! subroutine SoilTemperature in module SoilTemperatureMod
3449 ! subroutine SoilWater in module HydrologyMod
3450 !
3451 ! !REVISION HISTORY:
3452 ! 15 September 1999: Yongjiu Dai; Initial code
3453 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
3454 ! 1 July 2003: Mariana Vertenstein; modified for vectorization
3455 !
3456 !EOP
3457 !
3458 ! !OTHER LOCAL VARIABLES:
3459 !
3460 integer :: j,ci,fc !indices
3461 real(r8) :: gam(lbc:ubc,lbj:ubj) !temporary
3462 real(r8) :: bet(lbc:ubc) !temporary
3463 !-----------------------------------------------------------------------
3465 ! Solve the matrix
3467 !dir$ concurrent
3468 !cdir nodep
3469 do fc = 1,numf
3470 ci = filter(fc)
3471 bet(ci) = b(ci,jtop(ci))
3472 end do
3474 do j = lbj, ubj
3475 !dir$ prefervector
3476 !dir$ concurrent
3477 !cdir nodep
3478 do fc = 1,numf
3479 ci = filter(fc)
3480 if (j >= jtop(ci)) then
3481 if (j == jtop(ci)) then
3482 u(ci,j) = r(ci,j) / bet(ci)
3483 else
3484 gam(ci,j) = c(ci,j-1) / bet(ci)
3485 bet(ci) = b(ci,j) - a(ci,j) * gam(ci,j)
3486 u(ci,j) = (r(ci,j) - a(ci,j)*u(ci,j-1)) / bet(ci)
3487 end if
3488 end if
3489 end do
3490 end do
3492 !Cray X1 unroll directive used here as work-around for compiler issue 2003/10/20
3493 !dir$ unroll 0
3494 do j = ubj-1,lbj,-1
3495 !dir$ prefervector
3496 !dir$ concurrent
3497 !cdir nodep
3498 do fc = 1,numf
3499 ci = filter(fc)
3500 if (j >= jtop(ci)) then
3501 u(ci,j) = u(ci,j) - gam(ci,j+1) * u(ci,j+1)
3502 end if
3503 end do
3504 end do
3506 end subroutine Tridiagonal
3509 subroutine SnowWater(lbc, ubc, num_snowc, filter_snowc, & !i
3510 num_nosnowc, filter_nosnowc, & !i
3511 snl,do_capsnow,qflx_snomelt,qflx_rain_grnd, & !i
3512 qflx_sub_snow,qflx_evap_grnd, & !i
3513 qflx_dew_snow,qflx_dew_grnd,dz, & !i
3514 h2osoi_ice,h2osoi_liq, & !i&o
3515 qflx_top_soil) !o
3516 !===============================================================================
3517 ! !DESCRIPTION:
3518 ! Evaluate the change of snow mass and the snow water onto soil.
3519 ! Water flow within snow is computed by an explicit and non-physical
3520 ! based scheme, which permits a part of liquid water over the holding
3521 ! capacity (a tentative value is used, i.e. equal to 0.033*porosity) to
3522 ! percolate into the underlying layer. Except for cases where the
3523 ! porosity of one of the two neighboring layers is less than 0.05, zero
3524 ! flow is assumed. The water flow out of the bottom of the snow pack will
3525 ! participate as the input of the soil water and runoff. This subroutine
3526 ! uses a filter for columns containing snow which must be constructed prior
3527 ! to being called.
3528 !
3529 ! !REVISION HISTORY:
3530 ! 15 September 1999: Yongjiu Dai; Initial code
3531 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
3532 ! 15 November 2000: Mariana Vertenstein
3533 ! 2/26/02, Peter Thornton: Migrated to new data structures.
3534 !=============================================================================
3535 ! !USES:
3536 ! use clmtype
3538 implicit none
3540 !in:
3541 integer, intent(in) :: lbc, ubc ! column bounds
3542 integer, intent(in) :: num_snowc ! number of snow points in column filter
3543 integer, intent(in) :: filter_snowc(ubc-lbc+1) ! column filter for snow points
3544 integer, intent(in) :: num_nosnowc ! number of non-snow points in column filter
3545 integer, intent(in) :: filter_nosnowc(ubc-lbc+1) ! column filter for non-snow points
3547 integer , intent(in) :: snl(1) !number of snow layers
3548 logical , intent(in) :: do_capsnow(1) !true => do snow capping
3549 real(r8), intent(in) :: qflx_snomelt(1) !snow melt (mm H2O /s)
3550 real(r8), intent(in) :: qflx_rain_grnd(1) !rain on ground after interception (mm H2O/s) [+]
3551 real(r8), intent(in) :: qflx_sub_snow(1) !sublimation rate from snow pack (mm H2O /s) [+]
3552 real(r8), intent(in) :: qflx_evap_grnd(1) !ground surface evaporation rate (mm H2O/s) [+]
3553 real(r8), intent(in) :: qflx_dew_snow(1) !surface dew added to snow pack (mm H2O /s) [+]
3554 real(r8), intent(in) :: qflx_dew_grnd(1) !ground surface dew formation (mm H2O /s) [+]
3555 real(r8), intent(in) :: dz(1,-nlevsnow+1:nlevsoil) !layer depth (m)
3558 !inout:
3560 real(r8), intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) !ice lens (kg/m2)
3561 real(r8), intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) !liquid water (kg/m2)
3563 !out:
3565 real(r8), intent(out) :: qflx_top_soil(1) !net water input into soil from top (mm/s)
3568 ! OTHER LOCAL VARIABLES:
3570 integer :: c, j, fc !do loop/array indices
3571 real(r8) :: qin(lbc:ubc) !water flow into the elmement (mm/s)
3572 real(r8) :: qout(lbc:ubc) !water flow out of the elmement (mm/s)
3573 real(r8) :: wgdif !ice mass after minus sublimation
3574 real(r8) :: vol_liq(lbc:ubc,-nlevsnow+1:0) !partial volume of liquid water in layer
3575 real(r8) :: vol_ice(lbc:ubc,-nlevsnow+1:0) !partial volume of ice lens in layer
3576 real(r8) :: eff_porosity(lbc:ubc,-nlevsnow+1:0) !effective porosity = porosity - vol_ice
3577 !-----------------------------------------------------------------------
3578 ! Renew the mass of ice lens (h2osoi_ice) and liquid (h2osoi_liq) in the
3579 ! surface snow layer resulting from sublimation (frost) / evaporation (condense)
3581 !dir$ concurrent
3582 !cdir nodep
3583 do fc = 1,num_snowc
3584 c = filter_snowc(fc)
3585 if (do_capsnow(c)) then
3586 wgdif = h2osoi_ice(c,snl(c)+1) - qflx_sub_snow(c)*dtime
3587 h2osoi_ice(c,snl(c)+1) = wgdif
3588 if (wgdif < 0.) then
3589 h2osoi_ice(c,snl(c)+1) = 0.
3590 h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) + wgdif
3591 end if
3592 h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) - qflx_evap_grnd(c) * dtime
3593 else
3594 wgdif = h2osoi_ice(c,snl(c)+1) + (qflx_dew_snow(c) - qflx_sub_snow(c)) * dtime
3595 h2osoi_ice(c,snl(c)+1) = wgdif
3596 if (wgdif < 0.) then
3597 h2osoi_ice(c,snl(c)+1) = 0.
3598 h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) + wgdif
3599 end if
3600 h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) + &
3601 (qflx_rain_grnd(c) + qflx_dew_grnd(c) - qflx_evap_grnd(c)) * dtime
3602 end if
3603 h2osoi_liq(c,snl(c)+1) = max(0._r8, h2osoi_liq(c,snl(c)+1))
3604 end do
3606 ! Porosity and partial volume
3608 do j = -nlevsnow+1, 0
3609 !dir$ concurrent
3610 !cdir nodep
3611 do fc = 1, num_snowc
3612 c = filter_snowc(fc)
3613 if (j >= snl(c)+1) then
3614 vol_ice(c,j) = min(1._r8, h2osoi_ice(c,j)/(dz(c,j)*denice))
3615 eff_porosity(c,j) = 1. - vol_ice(c,j)
3616 vol_liq(c,j) = min(eff_porosity(c,j),h2osoi_liq(c,j)/(dz(c,j)*denh2o))
3617 end if
3618 end do
3619 end do
3621 ! Capillary forces within snow are usually two or more orders of magnitude
3622 ! less than those of gravity. Only gravity terms are considered.
3623 ! the genernal expression for water flow is "K * ss**3", however,
3624 ! no effective parameterization for "K". Thus, a very simple consideration
3625 ! (not physically based) is introduced:
3626 ! when the liquid water of layer exceeds the layer's holding
3627 ! capacity, the excess meltwater adds to the underlying neighbor layer.
3629 qin(:) = 0._r8
3631 do j = -nlevsnow+1, 0
3632 !dir$ concurrent
3633 !cdir nodep
3634 do fc = 1, num_snowc
3635 c = filter_snowc(fc)
3636 if (j >= snl(c)+1) then
3637 h2osoi_liq(c,j) = h2osoi_liq(c,j) + qin(c)
3638 if (j <= -1) then
3639 ! No runoff over snow surface, just ponding on surface
3640 if (eff_porosity(c,j) < wimp .OR. eff_porosity(c,j+1) < wimp) then
3641 qout(c) = 0._r8
3642 else
3643 qout(c) = max(0._r8,(vol_liq(c,j)-ssi*eff_porosity(c,j))*dz(c,j))
3644 qout(c) = min(qout(c),(1.-vol_ice(c,j+1)-vol_liq(c,j+1))*dz(c,j+1))
3645 end if
3646 else
3647 qout(c) = max(0._r8,(vol_liq(c,j) - ssi*eff_porosity(c,j))*dz(c,j))
3648 end if
3649 qout(c) = qout(c)*1000.
3650 h2osoi_liq(c,j) = h2osoi_liq(c,j) - qout(c)
3651 qin(c) = qout(c)
3652 end if
3653 end do
3654 end do
3656 !dir$ concurrent
3657 !cdir nodep
3658 do fc = 1, num_snowc
3659 c = filter_snowc(fc)
3660 ! Qout from snow bottom
3661 qflx_top_soil(c) = qout(c) / dtime
3662 end do
3664 !dir$ concurrent
3665 !cdir nodep
3666 do fc = 1, num_nosnowc
3667 c = filter_nosnowc(fc)
3668 qflx_top_soil(c) = qflx_rain_grnd(c) + qflx_snomelt(c)
3669 end do
3671 end subroutine SnowWater
3673 subroutine SnowCompaction(lbc, ubc, num_snowc, filter_snowc, &!i
3674 snl,imelt,frac_iceold,t_soisno, &!i
3675 h2osoi_ice,h2osoi_liq, &!i
3676 dz) !i&o
3679 !================================================================================
3680 ! !DESCRIPTION:
3681 ! Determine the change in snow layer thickness due to compaction and
3682 ! settling.
3683 ! Three metamorphisms of changing snow characteristics are implemented,
3684 ! i.e., destructive, overburden, and melt. The treatments of the former
3685 ! two are from SNTHERM.89 and SNTHERM.99 (1991, 1999). The contribution
3686 ! due to melt metamorphism is simply taken as a ratio of snow ice
3687 ! fraction after the melting versus before the melting.
3688 !
3689 ! CALLED FROM:
3690 ! subroutine Hydrology2 in module Hydrology2Mod
3691 !
3692 ! REVISION HISTORY:
3693 ! 15 September 1999: Yongjiu Dai; Initial code
3694 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
3695 ! 2/28/02, Peter Thornton: Migrated to new data structures
3696 !==============================================================================
3697 ! USES:
3698 ! use clmtype
3699 !
3700 ! !ARGUMENTS:
3701 implicit none
3703 !in:
3704 integer, intent(in) :: lbc, ubc ! column bounds
3705 integer, intent(in) :: num_snowc ! number of column snow points in column filter
3706 integer, intent(in) :: filter_snowc(ubc-lbc+1) ! column filter for snow points
3707 integer, intent(in) :: snl(1) !number of snow layers
3708 integer, intent(in) :: imelt(1,-nlevsnow+1:nlevsoil) !flag for melting (=1), freezing (=2), Not=0
3709 real(r8), intent(in) :: frac_iceold(1,-nlevsnow+1:nlevsoil) !fraction of ice relative to the tot water
3710 real(r8), intent(in) :: t_soisno(1,-nlevsnow+1:nlevsoil) !soil temperature (Kelvin)
3711 real(r8), intent(in) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) !ice lens (kg/m2)
3712 real(r8), intent(in) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) !liquid water (kg/m2)
3714 !inout:
3716 real(r8), intent(inout) :: dz(1,-nlevsnow+1:nlevsoil) !layer depth (m)
3718 ! OTHER LOCAL VARIABLES:
3720 integer :: j, c, fc ! indices
3721 real(r8), parameter :: c2 = 23.e-3 ! [m3/kg]
3722 real(r8), parameter :: c3 = 2.777e-6 ! [1/s]
3723 real(r8), parameter :: c4 = 0.04 ! [1/K]
3724 real(r8), parameter :: c5 = 2.0 !
3725 real(r8), parameter :: dm = 100.0 ! Upper Limit on Destructive Metamorphism Compaction [kg/m3]
3726 real(r8), parameter :: eta0 = 9.e+5 ! The Viscosity Coefficient Eta0 [kg-s/m2]
3727 real(r8) :: burden(lbc:ubc) ! pressure of overlying snow [kg/m2]
3728 real(r8) :: ddz1 ! Rate of settling of snowpack due to destructive metamorphism.
3729 real(r8) :: ddz2 ! Rate of compaction of snowpack due to overburden.
3730 real(r8) :: ddz3 ! Rate of compaction of snowpack due to melt [1/s]
3731 real(r8) :: dexpf ! expf=exp(-c4*(273.15-t_soisno)).
3732 real(r8) :: fi ! Fraction of ice relative to the total water content at current time step
3733 real(r8) :: td ! t_soisno - tfrz [K]
3734 real(r8) :: pdzdtc ! Nodal rate of change in fractional-thickness due to compaction [fraction/s]
3735 real(r8) :: void ! void (1 - vol_ice - vol_liq)
3736 real(r8) :: wx ! water mass (ice+liquid) [kg/m2]
3737 real(r8) :: bi ! partial density of ice [kg/m3]
3739 !-----------------------------------------------------------------------
3742 ! Begin calculation - note that the following column loops are only invoked if snl(c) < 0
3744 burden(:) = 0._r8
3746 do j = -nlevsnow+1, 0
3747 !dir$ concurrent
3748 !cdir nodep
3749 do fc = 1, num_snowc
3750 c = filter_snowc(fc)
3751 if (j >= snl(c)+1) then
3753 wx = h2osoi_ice(c,j) + h2osoi_liq(c,j)
3754 void = 1. - (h2osoi_ice(c,j)/denice + h2osoi_liq(c,j)/denh2o) / dz(c,j)
3756 ! Allow compaction only for non-saturated node and higher ice lens node.
3757 if (void > 0.001 .and. h2osoi_ice(c,j) > .1) then
3758 bi = h2osoi_ice(c,j) / dz(c,j)
3759 fi = h2osoi_ice(c,j) / wx
3760 td = tfrz-t_soisno(c,j)
3761 dexpf = exp(-c4*td)
3763 ! Settling as a result of destructive metamorphism
3765 ddz1 = -c3*dexpf
3766 if (bi > dm) ddz1 = ddz1*exp(-46.0e-3*(bi-dm))
3768 ! Liquid water term
3770 if (h2osoi_liq(c,j) > 0.01*dz(c,j)) ddz1=ddz1*c5
3772 ! Compaction due to overburden
3774 ddz2 = -burden(c)*exp(-0.08*td - c2*bi)/eta0
3776 ! Compaction occurring during melt
3778 if (imelt(c,j) == 1) then
3779 ddz3 = - 1./dtime * max(0._r8,(frac_iceold(c,j) - fi)/frac_iceold(c,j))
3780 else
3781 ddz3 = 0._r8
3782 end if
3784 ! Time rate of fractional change in dz (units of s-1)
3786 pdzdtc = ddz1 + ddz2 + ddz3
3788 ! The change in dz due to compaction
3790 dz(c,j) = dz(c,j) * (1.+pdzdtc*dtime)
3791 end if
3793 ! Pressure of overlying snow
3795 burden(c) = burden(c) + wx
3797 end if
3798 end do
3799 end do
3801 end subroutine SnowCompaction
3803 subroutine CombineSnowLayers(lbc, ubc, & !i
3804 num_snowc, filter_snowc, & !i&o
3805 snl,h2osno,snowdp,dz,zi, & !i&o
3806 t_soisno,h2osoi_ice,h2osoi_liq, & !i&o
3807 z) !o
3808 !==========================================================================
3809 ! !DESCRIPTION:
3810 ! Combine snow layers that are less than a minimum thickness or mass
3811 ! If the snow element thickness or mass is less than a prescribed minimum,
3812 ! then it is combined with a neighboring element. The subroutine
3813 ! clm\_combo.f90 then executes the combination of mass and energy.
3814 ! !CALLED FROM:
3815 ! subroutine Hydrology2 in module Hydrology2Mod
3816 !
3817 ! !REVISION HISTORY:
3818 ! 15 September 1999: Yongjiu Dai; Initial code
3819 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
3820 ! 2/28/02, Peter Thornton: Migrated to new data structures.
3821 !=========================================================================
3822 ! !USES:
3823 ! use clmtype
3824 !
3825 ! !ARGUMENTS:
3826 implicit none
3827 !in:
3828 integer, intent(in) :: lbc, ubc ! column bounds
3829 ! integer, intent(in) :: clandunit(1) !landunit index for each column
3830 ! integer, intent(in) :: ityplun(1) !landunit type
3832 !inout:
3833 integer, intent(inout) :: num_snowc ! number of column snow points in column filter
3834 integer, intent(inout) :: filter_snowc(ubc-lbc+1) ! column filter for snow points
3835 integer , intent(inout) :: snl(1) !number of snow layers
3836 real(r8), intent(inout) :: h2osno(1) !snow water (mm H2O)
3837 real(r8), intent(inout) :: snowdp(1) !snow height (m)
3838 real(r8), intent(inout) :: dz(1,-nlevsnow+1:nlevsoil) !layer depth (m)
3839 real(r8), intent(inout) :: zi(1,-nlevsnow+0:nlevsoil) !interface level below a "z" level (m)
3840 real(r8), intent(inout) :: t_soisno(1,-nlevsnow+1:nlevsoil) !soil temperature (Kelvin)
3841 real(r8), intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) !ice lens (kg/m2)
3842 real(r8), intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) !liquid water (kg/m2)
3844 !out:
3846 real(r8), intent(out) :: z(1,-nlevsnow+1:nlevsoil) !layer thickness (m)
3847 !
3848 !EOP
3849 !
3850 ! !OTHER LOCAL VARIABLES:
3851 !
3852 integer :: c, fc ! column indices
3853 integer :: i,k ! loop indices
3854 integer :: j,l ! node indices
3855 integer :: msn_old(lbc:ubc) ! number of top snow layer
3856 integer :: mssi(lbc:ubc) ! node index
3857 integer :: neibor ! adjacent node selected for combination
3858 real(r8):: zwice(lbc:ubc) ! total ice mass in snow
3859 real(r8):: zwliq (lbc:ubc) ! total liquid water in snow
3860 real(r8):: dzmin(5) ! minimum of top snow layer
3862 data dzmin /0.010, 0.015, 0.025, 0.055, 0.115/
3863 !-----------------------------------------------------------------------
3865 ! Check the mass of ice lens of snow, when the total is less than a small value,
3866 ! combine it with the underlying neighbor.
3868 !dir$ concurrent
3869 !cdir nodep
3870 do fc = 1, num_snowc
3871 c = filter_snowc(fc)
3872 msn_old(c) = snl(c)
3873 end do
3875 ! The following loop is NOT VECTORIZED
3877 do fc = 1, num_snowc
3878 c = filter_snowc(fc)
3879 ! l = clandunit(c)
3880 do j = msn_old(c)+1,0
3881 if (h2osoi_ice(c,j) <= .1) then
3882 ! if (ityplun(l) == istsoil) then
3883 ! h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + h2osoi_liq(c,j)
3884 ! h2osoi_ice(c,j+1) = h2osoi_ice(c,j+1) + h2osoi_ice(c,j)
3885 ! else if (ityplun(l) /= istsoil .and. j /= 0) then
3886 h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + h2osoi_liq(c,j)
3887 h2osoi_ice(c,j+1) = h2osoi_ice(c,j+1) + h2osoi_ice(c,j)
3888 ! end if
3890 ! shift all elements above this down one.
3891 if (j > snl(c)+1 .and. snl(c) < -1) then
3892 do i = j, snl(c)+2, -1
3893 t_soisno(c,i) = t_soisno(c,i-1)
3894 h2osoi_liq(c,i) = h2osoi_liq(c,i-1)
3895 h2osoi_ice(c,i) = h2osoi_ice(c,i-1)
3896 dz(c,i) = dz(c,i-1)
3897 end do
3898 end if
3899 snl(c) = snl(c) + 1
3900 end if
3901 end do
3902 end do
3904 !dir$ concurrent
3905 !cdir nodep
3906 do fc = 1, num_snowc
3907 c = filter_snowc(fc)
3908 h2osno(c) = 0._r8
3909 snowdp(c) = 0._r8
3910 zwice(c) = 0._r8
3911 zwliq(c) = 0._r8
3912 end do
3914 do j = -nlevsnow+1,0
3915 !dir$ concurrent
3916 !cdir nodep
3917 do fc = 1, num_snowc
3918 c = filter_snowc(fc)
3919 if (j >= snl(c)+1) then
3920 h2osno(c) = h2osno(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
3921 snowdp(c) = snowdp(c) + dz(c,j)
3922 zwice(c) = zwice(c) + h2osoi_ice(c,j)
3923 zwliq(c) = zwliq(c) + h2osoi_liq(c,j)
3924 end if
3925 end do
3926 end do
3928 ! Check the snow depth - all snow gone
3929 ! The liquid water assumes ponding on soil surface.
3931 !dir$ concurrent
3932 !cdir nodep
3933 do fc = 1, num_snowc
3934 c = filter_snowc(fc)
3935 ! l = clandunit(c)
3936 if (snowdp(c) < 0.01 .and. snowdp(c) > 0.) then
3937 snl(c) = 0
3938 h2osno(c) = zwice(c)
3939 if (h2osno(c) <= 0.) snowdp(c) = 0._r8
3940 ! if (ityplun(l) == istsoil) h2osoi_liq(c,1) = h2osoi_liq(c,1) + zwliq(c) !change by guhp
3941 end if
3942 end do
3944 ! Check the snow depth - snow layers combined
3945 ! The following loop IS NOT VECTORIZED
3947 do fc = 1, num_snowc
3948 c = filter_snowc(fc)
3950 ! Two or more layers
3952 if (snl(c) < -1) then
3954 msn_old(c) = snl(c)
3955 mssi(c) = 1
3957 do i = msn_old(c)+1,0
3958 if (dz(c,i) < dzmin(mssi(c))) then
3960 if (i == snl(c)+1) then
3961 ! If top node is removed, combine with bottom neighbor.
3962 neibor = i + 1
3963 else if (i == 0) then
3964 ! If the bottom neighbor is not snow, combine with the top neighbor.
3965 neibor = i - 1
3966 else
3967 ! If none of the above special cases apply, combine with the thinnest neighbor
3968 neibor = i + 1
3969 if ((dz(c,i-1)+dz(c,i)) < (dz(c,i+1)+dz(c,i))) neibor = i-1
3970 end if
3972 ! Node l and j are combined and stored as node j.
3973 if (neibor > i) then
3974 j = neibor
3975 l = i
3976 else
3977 j = i
3978 l = neibor
3979 end if
3981 call Combo (dz(c,j), h2osoi_liq(c,j), h2osoi_ice(c,j), &
3982 t_soisno(c,j), dz(c,l), h2osoi_liq(c,l), h2osoi_ice(c,l), t_soisno(c,l) )
3984 ! Now shift all elements above this down one.
3985 if (j-1 > snl(c)+1) then
3986 do k = j-1, snl(c)+2, -1
3987 t_soisno(c,k) = t_soisno(c,k-1)
3988 h2osoi_ice(c,k) = h2osoi_ice(c,k-1)
3989 h2osoi_liq(c,k) = h2osoi_liq(c,k-1)
3990 dz(c,k) = dz(c,k-1)
3991 end do
3992 end if
3994 ! Decrease the number of snow layers
3995 snl(c) = snl(c) + 1
3996 if (snl(c) >= -1) EXIT
3998 else
4000 ! The layer thickness is greater than the prescribed minimum value
4001 mssi(c) = mssi(c) + 1
4003 end if
4004 end do
4006 end if
4008 end do
4010 ! Reset the node depth and the depth of layer interface
4012 do j = 0, -nlevsnow+1, -1
4013 !dir$ concurrent
4014 !cdir nodep
4015 do fc = 1, num_snowc
4016 c = filter_snowc(fc)
4017 if (j >= snl(c) + 1) then
4018 z(c,j) = zi(c,j) - 0.5*dz(c,j)
4019 zi(c,j-1) = zi(c,j) - dz(c,j)
4020 end if
4021 end do
4022 end do
4024 end subroutine CombineSnowLayers
4026 subroutine DivideSnowLayers(lbc, ubc, & !i
4027 num_snowc, filter_snowc, & !i&o
4028 snl,dz,zi,t_soisno, & !i&o
4029 h2osoi_ice,h2osoi_liq, & !i&o
4030 z) !o
4033 !============================================================================
4034 ! !DESCRIPTION:
4035 ! Subdivides snow layers if they exceed their prescribed maximum thickness.
4036 ! !CALLED FROM:
4037 ! subroutine Hydrology2 in module Hydrology2Mod
4038 !
4039 ! !REVISION HISTORY:
4040 ! 15 September 1999: Yongjiu Dai; Initial code
4041 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
4042 ! 2/28/02, Peter Thornton: Migrated to new data structures.
4043 !============================================================================
4044 ! !USES:
4045 ! use clmtype
4046 !
4047 ! !ARGUMENTS:
4048 implicit none
4050 !in:
4051 integer, intent(in) :: lbc, ubc ! column bounds
4053 !inout:
4055 integer, intent(inout) :: num_snowc ! number of column snow points in column filter
4056 integer, intent(inout) :: filter_snowc(ubc-lbc+1) ! column filter for snow points
4057 integer , intent(inout) :: snl(1) !number of snow layers
4058 real(r8), intent(inout) :: dz(1,-nlevsnow+1:nlevsoil) !layer depth (m)
4059 real(r8), intent(inout) :: zi(1,-nlevsnow+0:nlevsoil) !interface level below a "z" level (m)
4060 real(r8), intent(inout) :: t_soisno(1,-nlevsnow+1:nlevsoil) !soil temperature (Kelvin)
4061 real(r8), intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) !ice lens (kg/m2)
4062 real(r8), intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) !liquid water (kg/m2)
4064 !out:
4066 real(r8), intent(out) :: z(1,-nlevsnow+1:nlevsoil) !layer thickness (m)
4070 ! OTHER LOCAL VARIABLES:
4072 integer :: j, c, fc ! indices
4073 real(r8) :: drr ! thickness of the combined [m]
4074 integer :: msno ! number of snow layer 1 (top) to msno (bottom)
4075 real(r8) :: dzsno(lbc:ubc,nlevsnow) ! Snow layer thickness [m]
4076 real(r8) :: swice(lbc:ubc,nlevsnow) ! Partial volume of ice [m3/m3]
4077 real(r8) :: swliq(lbc:ubc,nlevsnow) ! Partial volume of liquid water [m3/m3]
4078 real(r8) :: tsno(lbc:ubc ,nlevsnow) ! Nodel temperature [K]
4079 real(r8) :: zwice ! temporary
4080 real(r8) :: zwliq ! temporary
4081 real(r8) :: propor ! temporary
4082 !-----------------------------------------------------------------------
4084 ! Begin calculation - note that the following column loops are only invoked
4085 ! for snow-covered columns
4087 do j = 1,nlevsnow
4088 !dir$ concurrent
4089 !cdir nodep
4090 do fc = 1, num_snowc
4091 c = filter_snowc(fc)
4092 if (j <= abs(snl(c))) then
4093 dzsno(c,j) = dz(c,j+snl(c))
4094 swice(c,j) = h2osoi_ice(c,j+snl(c))
4095 swliq(c,j) = h2osoi_liq(c,j+snl(c))
4096 tsno(c,j) = t_soisno(c,j+snl(c))
4097 end if
4098 end do
4099 end do
4101 !dir$ concurrent
4102 !cdir nodep
4103 do fc = 1, num_snowc
4104 c = filter_snowc(fc)
4106 msno = abs(snl(c))
4108 if (msno == 1) then
4109 ! Specify a new snow layer
4110 if (dzsno(c,1) > 0.03) then
4111 msno = 2
4112 dzsno(c,1) = dzsno(c,1)/2.
4113 swice(c,1) = swice(c,1)/2.
4114 swliq(c,1) = swliq(c,1)/2.
4115 dzsno(c,2) = dzsno(c,1)
4116 swice(c,2) = swice(c,1)
4117 swliq(c,2) = swliq(c,1)
4118 tsno(c,2) = tsno(c,1)
4119 end if
4120 end if
4122 if (msno > 1) then
4123 if (dzsno(c,1) > 0.02) then
4124 drr = dzsno(c,1) - 0.02
4125 propor = drr/dzsno(c,1)
4126 zwice = propor*swice(c,1)
4127 zwliq = propor*swliq(c,1)
4128 propor = 0.02/dzsno(c,1)
4129 swice(c,1) = propor*swice(c,1)
4130 swliq(c,1) = propor*swliq(c,1)
4131 dzsno(c,1) = 0.02
4133 call Combo (dzsno(c,2), swliq(c,2), swice(c,2), tsno(c,2), drr, &
4134 zwliq, zwice, tsno(c,1))
4136 ! Subdivide a new layer
4137 if (msno <= 2 .and. dzsno(c,2) > 0.07) then
4138 msno = 3
4139 dzsno(c,2) = dzsno(c,2)/2.
4140 swice(c,2) = swice(c,2)/2.
4141 swliq(c,2) = swliq(c,2)/2.
4142 dzsno(c,3) = dzsno(c,2)
4143 swice(c,3) = swice(c,2)
4144 swliq(c,3) = swliq(c,2)
4145 tsno(c,3) = tsno(c,2)
4146 end if
4147 end if
4148 end if
4150 if (msno > 2) then
4151 if (dzsno(c,2) > 0.05) then
4152 drr = dzsno(c,2) - 0.05
4153 propor = drr/dzsno(c,2)
4154 zwice = propor*swice(c,2)
4155 zwliq = propor*swliq(c,2)
4156 propor = 0.05/dzsno(c,2)
4157 swice(c,2) = propor*swice(c,2)
4158 swliq(c,2) = propor*swliq(c,2)
4159 dzsno(c,2) = 0.05
4161 call Combo (dzsno(c,3), swliq(c,3), swice(c,3), tsno(c,3), drr, &
4162 zwliq, zwice, tsno(c,2))
4164 ! Subdivided a new layer
4165 if (msno <= 3 .and. dzsno(c,3) > 0.18) then
4166 msno = 4
4167 dzsno(c,3) = dzsno(c,3)/2.
4168 swice(c,3) = swice(c,3)/2.
4169 swliq(c,3) = swliq(c,3)/2.
4170 dzsno(c,4) = dzsno(c,3)
4171 swice(c,4) = swice(c,3)
4172 swliq(c,4) = swliq(c,3)
4173 tsno(c,4) = tsno(c,3)
4174 end if
4175 end if
4176 end if
4178 if (msno > 3) then
4179 if (dzsno(c,3) > 0.11) then
4180 drr = dzsno(c,3) - 0.11
4181 propor = drr/dzsno(c,3)
4182 zwice = propor*swice(c,3)
4183 zwliq = propor*swliq(c,3)
4184 propor = 0.11/dzsno(c,3)
4185 swice(c,3) = propor*swice(c,3)
4186 swliq(c,3) = propor*swliq(c,3)
4187 dzsno(c,3) = 0.11
4189 call Combo (dzsno(c,4), swliq(c,4), swice(c,4), tsno(c,4), drr, &
4190 zwliq, zwice, tsno(c,3))
4192 ! Subdivided a new layer
4193 if (msno <= 4 .and. dzsno(c,4) > 0.41) then
4194 msno = 5
4195 dzsno(c,4) = dzsno(c,4)/2.
4196 swice(c,4) = swice(c,4)/2.
4197 swliq(c,4) = swliq(c,4)/2.
4198 dzsno(c,5) = dzsno(c,4)
4199 swice(c,5) = swice(c,4)
4200 swliq(c,5) = swliq(c,4)
4201 tsno(c,5) = tsno(c,4)
4202 end if
4203 end if
4204 end if
4206 if (msno > 4) then
4207 if (dzsno(c,4) > 0.23) then
4208 drr = dzsno(c,4) - 0.23
4209 propor = drr/dzsno(c,4)
4210 zwice = propor*swice(c,4)
4211 zwliq = propor*swliq(c,4)
4212 propor = 0.23/dzsno(c,4)
4213 swice(c,4) = propor*swice(c,4)
4214 swliq(c,4) = propor*swliq(c,4)
4215 dzsno(c,4) = 0.23
4217 call Combo (dzsno(c,5), swliq(c,5), swice(c,5), tsno(c,5), drr, &
4218 zwliq, zwice, tsno(c,4))
4219 end if
4220 end if
4222 snl(c) = -msno
4224 end do
4226 do j = -nlevsnow+1,0
4227 !dir$ concurrent
4228 !cdir nodep
4229 do fc = 1, num_snowc
4230 c = filter_snowc(fc)
4231 if (j >= snl(c)+1) then
4232 dz(c,j) = dzsno(c,j-snl(c))
4233 h2osoi_ice(c,j) = swice(c,j-snl(c))
4234 h2osoi_liq(c,j) = swliq(c,j-snl(c))
4235 t_soisno(c,j) = tsno(c,j-snl(c))
4236 end if
4237 end do
4238 end do
4240 do j = 0, -nlevsnow+1, -1
4241 !dir$ concurrent
4242 !cdir nodep
4243 do fc = 1, num_snowc
4244 c = filter_snowc(fc)
4245 if (j >= snl(c)+1) then
4246 z(c,j) = zi(c,j) - 0.5*dz(c,j)
4247 zi(c,j-1) = zi(c,j) - dz(c,j)
4248 end if
4249 end do
4250 end do
4252 end subroutine DivideSnowLayers
4254 subroutine Combo(dz, wliq, wice, t, dz2, wliq2, wice2, t2)
4255 !
4256 ! !DESCRIPTION:
4257 ! Combines two elements and returns the following combined
4258 ! variables: dz, t, wliq, wice.
4259 ! The combined temperature is based on the equation:
4260 ! the sum of the enthalpies of the two elements =
4261 ! that of the combined element.
4262 !
4263 ! !USES:
4264 !
4265 ! !ARGUMENTS:
4266 implicit none
4267 real(r8), intent(in) :: dz2 ! nodal thickness of 2 elements being combined [m]
4268 real(r8), intent(in) :: wliq2 ! liquid water of element 2 [kg/m2]
4269 real(r8), intent(in) :: wice2 ! ice of element 2 [kg/m2]
4270 real(r8), intent(in) :: t2 ! nodal temperature of element 2 [K]
4271 real(r8), intent(inout) :: dz ! nodal thickness of 1 elements being combined [m]
4272 real(r8), intent(inout) :: wliq ! liquid water of element 1
4273 real(r8), intent(inout) :: wice ! ice of element 1 [kg/m2]
4274 real(r8), intent(inout) :: t ! nodel temperature of elment 1 [K]
4275 !
4276 ! !CALLED FROM:
4277 ! subroutine CombineSnowLayers in this module
4278 ! subroutine DivideSnowLayers in this module
4279 !
4280 ! !REVISION HISTORY:
4281 ! 15 September 1999: Yongjiu Dai; Initial code
4282 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
4283 !
4284 !EOP
4285 !
4286 ! !LOCAL VARIABLES:
4287 !
4288 real(r8) :: dzc ! Total thickness of nodes 1 and 2 (dzc=dz+dz2).
4289 real(r8) :: wliqc ! Combined liquid water [kg/m2]
4290 real(r8) :: wicec ! Combined ice [kg/m2]
4291 real(r8) :: tc ! Combined node temperature [K]
4292 real(r8) :: h ! enthalpy of element 1 [J/m2]
4293 real(r8) :: h2 ! enthalpy of element 2 [J/m2]
4294 real(r8) :: hc ! temporary
4295 !-----------------------------------------------------------------------
4297 dzc = dz+dz2
4298 wicec = (wice+wice2)
4299 wliqc = (wliq+wliq2)
4300 h = (cpice*wice+cpliq*wliq) * (t-tfrz)+hfus*wliq
4301 h2= (cpice*wice2+cpliq*wliq2) * (t2-tfrz)+hfus*wliq2
4303 hc = h + h2
4304 if(hc < 0.)then
4305 tc = tfrz + hc/(cpice*wicec + cpliq*wliqc)
4306 else if (hc.le.hfus*wliqc) then
4307 tc = tfrz
4308 else
4309 tc = tfrz + (hc - hfus*wliqc) / (cpice*wicec + cpliq*wliqc)
4310 end if
4312 dz = dzc
4313 wice = wicec
4314 wliq = wliqc
4315 t = tc
4317 end subroutine Combo
4319 subroutine BuildSnowFilter(lbc, ubc, num_nolakec, filter_nolakec,snl, & !i
4320 num_snowc, filter_snowc, & !o
4321 num_nosnowc, filter_nosnowc) !o
4322 !
4323 ! !DESCRIPTION:
4324 ! Constructs snow filter for use in vectorized loops for snow hydrology.
4325 !
4326 ! !USES:
4327 ! use clmtype
4328 !
4329 ! !ARGUMENTS:
4330 implicit none
4331 integer, intent(in) :: lbc, ubc ! column bounds
4332 integer, intent(in) :: num_nolakec ! number of column non-lake points in column filter
4333 integer, intent(in) :: filter_nolakec(ubc-lbc+1) ! column filter for non-lake points
4334 integer, intent(in) :: snl(1) ! number of snow layers
4335 integer, intent(out) :: num_snowc ! number of column snow points in column filter
4336 integer, intent(out) :: filter_snowc(ubc-lbc+1) ! column filter for snow points
4337 integer, intent(out) :: num_nosnowc ! number of column non-snow points in column filter
4338 integer, intent(out) :: filter_nosnowc(ubc-lbc+1) ! column filter for non-snow points
4339 !
4340 ! !CALLED FROM:
4341 ! subroutine Hydrology2 in Hydrology2Mod
4342 ! subroutine CombineSnowLayers in this module
4343 !
4344 ! !REVISION HISTORY:
4345 ! 2003 July 31: Forrest Hoffman
4346 !
4347 ! !LOCAL VARIABLES:
4348 ! local pointers to implicit in arguments
4349 !
4350 !EOP
4351 !
4352 ! !OTHER LOCAL VARIABLES:
4353 integer :: fc, c
4354 !-----------------------------------------------------------------------
4357 ! Build snow/no-snow filters for other subroutines
4359 num_snowc = 0
4360 num_nosnowc = 0
4361 do fc = 1, num_nolakec
4362 c = filter_nolakec(fc)
4363 if (snl(c) < 0) then
4364 num_snowc = num_snowc + 1
4365 filter_snowc(num_snowc) = c
4366 else
4367 num_nosnowc = num_nosnowc + 1
4368 filter_nosnowc(num_nosnowc) = c
4369 end if
4370 end do
4372 end subroutine BuildSnowFilter
4376 subroutine FrictionVelocity(pgridcell,forc_hgt,forc_hgt_u, & !i
4377 forc_hgt_t,forc_hgt_q, & !i
4378 lbp, ubp, fn, filterp, & !i
4379 displa, z0m, z0h, z0q, & !i
4380 obu, iter, ur, um, & !i
4381 ustar,temp1, temp2, temp12m, temp22m, & !o
4382 u10,fv, & !o
4383 fm) !i&o
4385 !=============================================================================
4386 ! !DESCRIPTION:
4387 ! Calculation of the friction velocity, relation for potential
4388 ! temperature and humidity profiles of surface boundary layer.
4389 ! The scheme is based on the work of Zeng et al. (1998):
4390 ! Intercomparison of bulk aerodynamic algorithms for the computation
4391 ! of sea surface fluxes using TOGA CORE and TAO data. J. Climate,
4392 ! Vol. 11, 2628-2644.
4393 !
4394 ! !REVISION HISTORY:
4395 ! 15 September 1999: Yongjiu Dai; Initial code
4396 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
4397 ! 12/19/01, Peter Thornton
4398 ! Added arguments to eliminate passing clm derived type into this function.
4399 ! Created by Mariana Vertenstein
4400 !============================================================================
4401 ! !USES:
4402 ! use clmtype
4403 !!use clm_atmlnd, only : clm_a2l
4404 !
4405 ! !ARGUMENTS:
4406 implicit none
4408 !in:
4410 integer , intent(in) :: pgridcell(1) ! pft's gridcell index
4411 real(r8), intent(in) :: forc_hgt(1) ! atmospheric reference height (m)
4412 real(r8), intent(in) :: forc_hgt_u(1) ! observational height of wind [m]
4413 real(r8), intent(in) :: forc_hgt_t(1) ! observational height of temperature [m]
4414 real(r8), intent(in) :: forc_hgt_q(1) ! observational height of humidity [m]
4415 integer , intent(in) :: lbp, ubp ! pft array bounds
4416 integer , intent(in) :: fn ! number of filtered pft elements
4417 integer , intent(in) :: filterp(fn) ! pft filter
4418 real(r8), intent(in) :: displa(lbp:ubp) ! displacement height (m)
4419 real(r8), intent(in) :: z0m(lbp:ubp) ! roughness length over vegetation, momentum [m]
4420 real(r8), intent(in) :: z0h(lbp:ubp) ! roughness length over vegetation, sensible heat [m]
4421 real(r8), intent(in) :: z0q(lbp:ubp) ! roughness length over vegetation, latent heat [m]
4422 real(r8), intent(in) :: obu(lbp:ubp) ! monin-obukhov length (m)
4423 integer, intent(in) :: iter ! iteration number
4424 real(r8), intent(in) :: ur(lbp:ubp) ! wind speed at reference height [m/s]
4425 real(r8), intent(in) :: um(lbp:ubp) ! wind speed including the stablity effect [m/s]
4427 !out:
4429 real(r8), intent(out) :: ustar(lbp:ubp) ! friction velocity [m/s]
4430 real(r8), intent(out) :: temp1(lbp:ubp) ! relation for potential temperature profile
4431 real(r8), intent(out) :: temp12m(lbp:ubp) ! relation for potential temperature profile applied at 2-m
4432 real(r8), intent(out) :: temp2(lbp:ubp) ! relation for specific humidity profile
4433 real(r8), intent(out) :: temp22m(lbp:ubp) ! relation for specific humidity profile applied at 2-m
4434 real(r8), intent(out) :: u10(1) ! 10-m wind (m/s) (for dust model)
4435 real(r8), intent(out) :: fv(1) ! friction velocity (m/s) (for dust model)
4437 !inout:
4438 real(r8), intent(inout) :: fm(lbp:ubp) ! needed for DGVM only to diagnose 10m wind
4440 ! OTHER LOCAL VARIABLES:
4442 real(r8), parameter :: zetam = 1.574_r8 ! transition point of flux-gradient relation (wind profile)
4443 real(r8), parameter :: zetat = 0.465_r8 ! transition point of flux-gradient relation (temp. profile)
4444 integer :: f ! pft-filter index
4445 integer :: p ! pft index
4446 integer :: g ! gridcell index
4447 real(r8):: zldis(lbp:ubp) ! reference height "minus" zero displacement heght [m]
4448 real(r8):: zeta(lbp:ubp) ! dimensionless height used in Monin-Obukhov theory
4449 #if (defined DGVM) || (defined DUST)
4450 real(r8) :: tmp1,tmp2,tmp3,tmp4 ! Used to diagnose the 10 meter wind
4451 real(r8) :: fmnew ! Used to diagnose the 10 meter wind
4452 real(r8) :: fm10 ! Used to diagnose the 10 meter wind
4453 real(r8) :: zeta10 ! Used to diagnose the 10 meter wind
4454 #endif
4455 !------------------------------------------------------------------------------
4458 ! Adjustment factors for unstable (moz < 0) or stable (moz > 0) conditions.
4460 #if (!defined PERGRO)
4462 !dir$ concurrent
4463 !cdir nodep
4464 do f = 1, fn
4465 p = filterp(f)
4466 g = pgridcell(p)
4468 ! Wind profile
4470 zldis(p) = forc_hgt_u(g)-displa(p)
4471 zeta(p) = zldis(p)/obu(p)
4472 if (zeta(p) < -zetam) then
4473 ustar(p) = vkc*um(p)/(log(-zetam*obu(p)/z0m(p))&
4474 - StabilityFunc1(-zetam) &
4475 + StabilityFunc1(z0m(p)/obu(p)) &
4476 + 1.14_r8*((-zeta(p))**0.333_r8-(zetam)**0.333_r8))
4477 else if (zeta(p) < 0._r8) then
4478 ustar(p) = vkc*um(p)/(log(zldis(p)/z0m(p))&
4479 - StabilityFunc1(zeta(p))&
4480 + StabilityFunc1(z0m(p)/obu(p)))
4481 else if (zeta(p) <= 1._r8) then
4482 ustar(p) = vkc*um(p)/(log(zldis(p)/z0m(p)) + 5._r8*zeta(p) -5._r8*z0m(p)/obu(p))
4483 else
4484 ustar(p) = vkc*um(p)/(log(obu(p)/z0m(p))+5._r8-5._r8*z0m(p)/obu(p) &
4485 +(5._r8*log(zeta(p))+zeta(p)-1._r8))
4486 end if
4488 ! Temperature profile
4490 zldis(p) = forc_hgt_t(g)-displa(p)
4491 zeta(p) = zldis(p)/obu(p)
4492 if (zeta(p) < -zetat) then
4493 temp1(p) = vkc/(log(-zetat*obu(p)/z0h(p))&
4494 - StabilityFunc2(-zetat) &
4495 + StabilityFunc2(z0h(p)/obu(p)) &
4496 + 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(p))**(-0.333_r8)))
4497 else if (zeta(p) < 0._r8) then
4498 temp1(p) = vkc/(log(zldis(p)/z0h(p)) &
4499 - StabilityFunc2(zeta(p)) &
4500 + StabilityFunc2(z0h(p)/obu(p)))
4501 else if (zeta(p) <= 1._r8) then
4502 temp1(p) = vkc/(log(zldis(p)/z0h(p)) + 5._r8*zeta(p) - 5._r8*z0h(p)/obu(p))
4503 else
4504 temp1(p) = vkc/(log(obu(p)/z0h(p)) + 5._r8 - 5._r8*z0h(p)/obu(p) &
4505 + (5._r8*log(zeta(p))+zeta(p)-1._r8))
4506 end if
4508 ! Humidity profile
4510 if (forc_hgt_q(g) == forc_hgt_t(g) .and. z0q(p) == z0h(p)) then
4511 temp2(p) = temp1(p)
4512 else
4513 zldis(p) = forc_hgt_q(g)-displa(p)
4514 zeta(p) = zldis(p)/obu(p)
4515 if (zeta(p) < -zetat) then
4516 temp2(p) = vkc/(log(-zetat*obu(p)/z0q(p)) &
4517 - StabilityFunc2(-zetat) &
4518 + StabilityFunc2(z0q(p)/obu(p)) &
4519 + 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(p))**(-0.333_r8)))
4520 else if (zeta(p) < 0._r8) then
4521 temp2(p) = vkc/(log(zldis(p)/z0q(p)) &
4522 - StabilityFunc2(zeta(p)) &
4523 + StabilityFunc2(z0q(p)/obu(p)))
4524 else if (zeta(p) <= 1._r8) then
4525 temp2(p) = vkc/(log(zldis(p)/z0q(p)) + 5._r8*zeta(p)-5._r8*z0q(p)/obu(p))
4526 else
4527 temp2(p) = vkc/(log(obu(p)/z0q(p)) + 5._r8 - 5._r8*z0q(p)/obu(p) &
4528 + (5._r8*log(zeta(p))+zeta(p)-1._r8))
4529 end if
4530 endif
4532 ! Temperature profile applied at 2-m
4534 zldis(p) = 2.0_r8 + z0h(p)
4535 zeta(p) = zldis(p)/obu(p)
4536 if (zeta(p) < -zetat) then
4537 temp12m(p) = vkc/(log(-zetat*obu(p)/z0h(p))&
4538 - StabilityFunc2(-zetat) &
4539 + StabilityFunc2(z0h(p)/obu(p)) &
4540 + 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(p))**(-0.333_r8)))
4541 else if (zeta(p) < 0._r8) then
4542 temp12m(p) = vkc/(log(zldis(p)/z0h(p)) &
4543 - StabilityFunc2(zeta(p)) &
4544 + StabilityFunc2(z0h(p)/obu(p)))
4545 else if (zeta(p) <= 1._r8) then
4546 temp12m(p) = vkc/(log(zldis(p)/z0h(p)) + 5._r8*zeta(p) - 5._r8*z0h(p)/obu(p))
4547 else
4548 temp12m(p) = vkc/(log(obu(p)/z0h(p)) + 5._r8 - 5._r8*z0h(p)/obu(p) &
4549 + (5._r8*log(zeta(p))+zeta(p)-1._r8))
4550 end if
4552 ! Humidity profile applied at 2-m
4554 if (z0q(p) == z0h(p)) then
4555 temp22m(p) = temp12m(p)
4556 else
4557 zldis(p) = 2.0_r8 + z0q(p)
4558 zeta(p) = zldis(p)/obu(p)
4559 if (zeta(p) < -zetat) then
4560 temp22m(p) = vkc/(log(-zetat*obu(p)/z0q(p)) - &
4561 StabilityFunc2(-zetat) + StabilityFunc2(z0q(p)/obu(p)) &
4562 + 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(p))**(-0.333_r8)))
4563 else if (zeta(p) < 0._r8) then
4564 temp22m(p) = vkc/(log(zldis(p)/z0q(p)) - &
4565 StabilityFunc2(zeta(p))+StabilityFunc2(z0q(p)/obu(p)))
4566 else if (zeta(p) <= 1._r8) then
4567 temp22m(p) = vkc/(log(zldis(p)/z0q(p)) + 5._r8*zeta(p)-5._r8*z0q(p)/obu(p))
4568 else
4569 temp22m(p) = vkc/(log(obu(p)/z0q(p)) + 5._r8 - 5._r8*z0q(p)/obu(p) &
4570 + (5._r8*log(zeta(p))+zeta(p)-1._r8))
4571 end if
4572 end if
4574 #if (defined DGVM) || (defined DUST)
4575 ! diagnose 10-m wind for dust model (dstmbl.F)
4576 ! Notes from C. Zender's dst.F:
4577 ! According to Bon96 p. 62, the displacement height d (here displa) is
4578 ! 0.0 <= d <= 0.34 m in dust source regions (i.e., regions w/o trees).
4579 ! Therefore d <= 0.034*z1 and may safely be neglected.
4580 ! Code from LSM routine SurfaceTemperature was used to obtain u10
4582 zldis(p) = forc_hgt_u(g)-displa(p)
4583 zeta(p) = zldis(p)/obu(p)
4584 if (min(zeta(p), 1._r8) < 0._r8) then
4585 tmp1 = (1._r8 - 16._r8*min(zeta(p),1._r8))**0.25_r8
4586 tmp2 = log((1._r8+tmp1*tmp1)/2._r8)
4587 tmp3 = log((1._r8+tmp1)/2._r8)
4588 fmnew = 2._r8*tmp3 + tmp2 - 2._r8*atan(tmp1) + 1.5707963_r8
4589 else
4590 fmnew = -5._r8*min(zeta(p),1._r8)
4591 endif
4592 if (iter == 1) then
4593 fm(p) = fmnew
4594 else
4595 fm(p) = 0.5_r8 * (fm(p)+fmnew)
4596 end if
4597 zeta10 = min(10._r8/obu(p), 1._r8)
4598 if (zeta(p) == 0._r8) zeta10 = 0._r8
4599 if (zeta10 < 0._r8) then
4600 tmp1 = (1.0_r8 - 16.0_r8 * zeta10)**0.25_r8
4601 tmp2 = log((1.0_r8 + tmp1*tmp1)/2.0_r8)
4602 tmp3 = log((1.0_r8 + tmp1)/2.0_r8)
4603 fm10 = 2.0_r8*tmp3 + tmp2 - 2.0_r8*atan(tmp1) + 1.5707963_r8
4604 else ! not stable
4605 fm10 = -5.0_r8 * zeta10
4606 end if
4607 tmp4 = log(forc_hgt(g) / 10._r8)
4608 u10(p) = ur(p) - ustar(p)/vkc * (tmp4 - fm(p) + fm10)
4609 fv(p) = ustar(p)
4610 #endif
4612 end do
4613 #endif
4616 #if (defined PERGRO)
4618 !===============================================================================
4619 ! The following only applies when PERGRO is defined
4620 !===============================================================================
4622 !dir$ concurrent
4623 !cdir nodep
4624 do f = 1, fn
4625 p = filterp(f)
4626 g = pgridcell(p)
4628 zldis(p) = forc_hgt_u(g)-displa(p)
4629 zeta(p) = zldis(p)/obu(p)
4630 if (zeta(p) < -zetam) then ! zeta < -1
4631 ustar(p) = vkc * um(p) / log(-zetam*obu(p)/z0m(p))
4632 else if (zeta(p) < 0._r8) then ! -1 <= zeta < 0
4633 ustar(p) = vkc * um(p) / log(zldis(p)/z0m(p))
4634 else if (zeta(p) <= 1._r8) then ! 0 <= ztea <= 1
4635 ustar(p)=vkc * um(p)/log(zldis(p)/z0m(p))
4636 else ! 1 < zeta, phi=5+zeta
4637 ustar(p)=vkc * um(p)/log(obu(p)/z0m(p))
4638 endif
4640 zldis(p) = forc_hgt_t(g)-displa(p)
4641 zeta(p) = zldis(p)/obu(p)
4642 if (zeta(p) < -zetat) then
4643 temp1(p)=vkc/log(-zetat*obu(p)/z0h(p))
4644 else if (zeta(p) < 0._r8) then
4645 temp1(p)=vkc/log(zldis(p)/z0h(p))
4646 else if (zeta(p) <= 1._r8) then
4647 temp1(p)=vkc/log(zldis(p)/z0h(p))
4648 else
4649 temp1(p)=vkc/log(obu(p)/z0h(p))
4650 end if
4652 zldis(p) = forc_hgt_q(g)-displa(p)
4653 zeta(p) = zldis(p)/obu(p)
4654 if (zeta(p) < -zetat) then
4655 temp2(p)=vkc/log(-zetat*obu(p)/z0q(p))
4656 else if (zeta(p) < 0._r8) then
4657 temp2(p)=vkc/log(zldis(p)/z0q(p))
4658 else if (zeta(p) <= 1._r8) then
4659 temp2(p)=vkc/log(zldis(p)/z0q(p))
4660 else
4661 temp2(p)=vkc/log(obu(p)/z0q(p))
4662 end if
4664 zldis(p) = 2.0_r8 + z0h(p)
4665 zeta(p) = zldis(p)/obu(p)
4666 if (zeta(p) < -zetat) then
4667 temp12m(p)=vkc/log(-zetat*obu(p)/z0h(p))
4668 else if (zeta(p) < 0._r8) then
4669 temp12m(p)=vkc/log(zldis(p)/z0h(p))
4670 else if (zeta(p) <= 1._r8) then
4671 temp12m(p)=vkc/log(zldis(p)/z0h(p))
4672 else
4673 temp12m(p)=vkc/log(obu(p)/z0h(p))
4674 end if
4676 zldis(p) = 2.0_r8 + z0q(p)
4677 zeta(p) = zldis(p)/obu(p)
4678 if (zeta(p) < -zetat) then
4679 temp22m(p)=vkc/log(-zetat*obu(p)/z0q(p))
4680 else if (zeta(p) < 0._r8) then
4681 temp22m(p)=vkc/log(zldis(p)/z0q(p))
4682 else if (zeta(p) <= 1._r8) then
4683 temp22m(p)=vkc/log(zldis(p)/z0q(p))
4684 else
4685 temp22m(p)=vkc/log(obu(p)/z0q(p))
4686 end if
4687 #if (defined DGVM) || (defined DUST)
4688 ! diagnose 10-m wind for dust model (dstmbl.F)
4689 ! Notes from C. Zender's dst.F:
4690 ! According to Bon96 p. 62, the displacement height d (here displa) is
4691 ! 0.0 <= d <= 0.34 m in dust source regions (i.e., regions w/o trees).
4692 ! Therefore d <= 0.034*z1 and may safely be neglected.
4693 ! Code from LSM routine SurfaceTemperature was used to obtain u10
4695 zldis(p) = forc_hgt_u(g)-displa(p)
4696 zeta(p) = zldis(p)/obu(p)
4697 if (min(zeta(p), 1._r8) < 0._r8) then
4698 tmp1 = (1._r8 - 16._r8*min(zeta(p),1._r8))**0.25_r8
4699 tmp2 = log((1._r8+tmp1*tmp1)/2._r8)
4700 tmp3 = log((1._r8+tmp1)/2._r8)
4701 fmnew = 2._r8*tmp3 + tmp2 - 2._r8*atan(tmp1) + 1.5707963_r8
4702 else
4703 fmnew = -5._r8*min(zeta(p),1._r8)
4704 endif
4705 if (iter == 1) then
4706 fm(p) = fmnew
4707 else
4708 fm(p) = 0.5_r8 * (fm(p)+fmnew)
4709 end if
4710 zeta10 = min(10._r8/obu(p), 1._r8)
4711 if (zeta(p) == 0._r8) zeta10 = 0._r8
4712 if (zeta10 < 0._r8) then
4713 tmp1 = (1.0_r8 - 16.0_r8 * zeta10)**0.25_r8
4714 tmp2 = log((1.0_r8 + tmp1*tmp1)/2.0_r8)
4715 tmp3 = log((1.0_r8 + tmp1)/2.0_r8)
4716 fm10 = 2.0_r8*tmp3 + tmp2 - 2.0_r8*atan(tmp1) + 1.5707963_r8
4717 else ! not stable
4718 fm10 = -5.0_r8 * zeta10
4719 end if
4720 tmp4 = log(forc_hgt(g) / 10._r8)
4721 u10(p) = ur(p) - ustar(p)/vkc * (tmp4 - fm(p) + fm10)
4722 fv(p) = ustar(p)
4723 #endif
4724 end do
4726 #endif
4728 end subroutine FrictionVelocity
4730 ! !IROUTINE: StabilityFunc
4731 !
4732 ! !INTERFACE:
4733 real(r8) function StabilityFunc1(zeta)
4734 !
4735 ! !DESCRIPTION:
4736 ! Stability function for rib < 0.
4737 !
4738 ! !USES:
4739 ! use shr_const_mod, only: SHR_CONST_PI
4740 !Zack Subin, 7/8/08
4741 !
4742 ! !ARGUMENTS:
4743 implicit none
4744 real(r8), intent(in) :: zeta ! dimensionless height used in Monin-Obukhov theory
4745 !
4746 ! !CALLED FROM:
4747 ! subroutine FrictionVelocity in this module
4748 !
4749 ! !REVISION HISTORY:
4750 ! 15 September 1999: Yongjiu Dai; Initial code
4751 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
4752 !
4753 !EOP
4754 !
4755 ! !LOCAL VARIABLES:
4756 real(r8) :: chik, chik2
4757 !------------------------------------------------------------------------------
4759 chik2 = sqrt(1._r8-16._r8*zeta)
4760 chik = sqrt(chik2)
4761 StabilityFunc1 = 2._r8*log((1._r8+chik)*0.5_r8) &
4762 !Changed to pie, Zack Subin, 7/9/08
4763 + log((1._r8+chik2)*0.5_r8)-2._r8*atan(chik)+pie*0.5_r8
4765 end function StabilityFunc1
4767 !------------------------------------------------------------------------------
4768 !BOP
4769 !
4770 ! !IROUTINE: StabilityFunc2
4771 !
4772 ! !INTERFACE:
4773 real(r8) function StabilityFunc2(zeta)
4774 !
4775 ! !DESCRIPTION:
4776 ! Stability function for rib < 0.
4777 !
4778 ! !USES:
4779 !Removed by Zack Subin, 7/9/08
4780 ! use shr_const_mod, only: SHR_CONST_PI
4781 !
4782 ! !ARGUMENTS:
4783 implicit none
4784 real(r8), intent(in) :: zeta ! dimensionless height used in Monin-Obukhov theory
4785 !
4786 ! !CALLED FROM:
4787 ! subroutine FrictionVelocity in this module
4788 !
4789 ! !REVISION HISTORY:
4790 ! 15 September 1999: Yongjiu Dai; Initial code
4791 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
4792 !
4793 !EOP
4794 !
4795 ! !LOCAL VARIABLES:
4796 real(r8) :: chik2
4797 !------------------------------------------------------------------------------
4799 chik2 = sqrt(1._r8-16._r8*zeta)
4800 StabilityFunc2 = 2._r8*log((1._r8+chik2)*0.5_r8)
4802 end function StabilityFunc2
4804 !-----------------------------------------------------------------------
4805 !BOP
4806 !
4807 ! !IROUTINE: MoninObukIni
4808 !
4809 ! !INTERFACE:
4810 subroutine MoninObukIni (ur, thv, dthv, zldis, z0m, um, obu)
4811 !
4812 ! !DESCRIPTION:
4813 ! Initialization of the Monin-Obukhov length.
4814 ! The scheme is based on the work of Zeng et al. (1998):
4815 ! Intercomparison of bulk aerodynamic algorithms for the computation
4816 ! of sea surface fluxes using TOGA CORE and TAO data. J. Climate,
4817 ! Vol. 11, 2628-2644.
4818 !
4819 ! !USES:
4820 !
4821 ! !ARGUMENTS:
4822 implicit none
4823 real(r8), intent(in) :: ur ! wind speed at reference height [m/s]
4824 real(r8), intent(in) :: thv ! virtual potential temperature (kelvin)
4825 real(r8), intent(in) :: dthv ! diff of vir. poten. temp. between ref. height and surface
4826 real(r8), intent(in) :: zldis ! reference height "minus" zero displacement heght [m]
4827 real(r8), intent(in) :: z0m ! roughness length, momentum [m]
4828 real(r8), intent(out) :: um ! wind speed including the stability effect [m/s]
4829 real(r8), intent(out) :: obu ! monin-obukhov length (m)
4830 !
4831 ! !CALLED FROM:
4832 ! subroutine BareGroundFluxes in module BareGroundFluxesMod.F90
4833 ! subroutine BiogeophysicsLake in module BiogeophysicsLakeMod.F90
4834 ! subroutine CanopyFluxes in module CanopyFluxesMod.F90
4835 !
4836 ! !REVISION HISTORY:
4837 ! 15 September 1999: Yongjiu Dai; Initial code
4838 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
4839 !
4840 !EOP
4841 !
4842 ! !LOCAL VARIABLES:
4843 !
4844 real(r8) :: wc ! convective velocity [m/s]
4845 real(r8) :: rib ! bulk Richardson number
4846 real(r8) :: zeta ! dimensionless height used in Monin-Obukhov theory
4847 real(r8) :: ustar ! friction velocity [m/s]
4848 !-----------------------------------------------------------------------
4850 ! Initial values of u* and convective velocity
4852 ustar=0.06_r8
4853 wc=0.5_r8
4854 if (dthv >= 0._r8) then
4855 um=max(ur,0.1_r8)
4856 else
4857 um=sqrt(ur*ur+wc*wc)
4858 endif
4860 rib=grav*zldis*dthv/(thv*um*um)
4861 #if (defined PERGRO)
4862 rib = 0._r8
4863 #endif
4865 if (rib >= 0._r8) then ! neutral or stable
4866 zeta = rib*log(zldis/z0m)/(1._r8-5._r8*min(rib,0.19_r8))
4867 zeta = min(2._r8,max(zeta,0.01_r8 ))
4868 else ! unstable
4869 zeta=rib*log(zldis/z0m)
4870 zeta = max(-100._r8,min(zeta,-0.01_r8 ))
4871 endif
4873 obu=zldis/zeta
4875 end subroutine MoninObukIni
4877 subroutine LakeDebug( str )
4879 IMPLICIT NONE
4880 CHARACTER*(*), str
4882 CALL wrf_debug( 0 , TRIM(str) )
4884 end subroutine LakeDebug
4886 SUBROUTINE lakeini(IVGTYP, ISLTYP, HT, SNOW, & !i
4887 lake_min_elev, restart, lakedepth_default, lake_depth, &
4888 lakedepth2d, savedtke12d, snowdp2d, h2osno2d, & !o
4889 snl2d, t_grnd2d, t_lake3d, lake_icefrac3d, &
4890 z_lake3d, dz_lake3d, t_soisno3d, h2osoi_ice3d, &
4891 h2osoi_liq3d, h2osoi_vol3d, z3d, dz3d, &
4892 zi3d, watsat3d, csol3d, tkmg3d, &
4893 iswater, xice, xice_threshold, xland, tsk, &
4894 #if (EM_CORE == 1)
4895 lakemask, lakeflag, &
4896 #endif
4897 lake_depth_flag, use_lakedepth, &
4898 tkdry3d, tksatu3d, lake, its, ite, jts, jte, &
4899 ims,ime, jms,jme)
4901 !==============================================================================
4902 ! This subroutine was first edited by Hongping Gu for coupling
4903 ! 07/20/2010
4904 !==============================================================================
4906 USE module_wrf_error
4907 implicit none
4909 INTEGER , INTENT (IN) :: iswater
4910 REAL, INTENT(IN) :: xice_threshold
4911 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: XICE
4912 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN):: TSK
4913 REAL, DIMENSION( ims:ime, jms:jme ) ,INTENT(INOUT) :: XLAND
4915 #if (EM_CORE == 1)
4916 REAL, DIMENSION( ims:ime , jms:jme ) :: LAKEMASK
4917 INTEGER , INTENT (IN) :: lakeflag
4918 #endif
4919 INTEGER , INTENT (INOUT) :: lake_depth_flag
4920 INTEGER , INTENT (IN) :: use_lakedepth
4922 LOGICAL , INTENT(IN) :: restart
4923 INTEGER, INTENT(IN ) :: ims,ime, jms,jme
4924 INTEGER, INTENT(IN ) :: its,ite, jts,jte
4925 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: IVGTYP, &
4926 ISLTYP
4927 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: HT
4928 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SNOW
4929 real, intent(in) :: lakedepth_default,lake_min_elev
4931 real, dimension(ims:ime,jms:jme ),intent(out) :: lakedepth2d, &
4932 savedtke12d
4933 real, dimension(ims:ime,jms:jme ),intent(out) :: snowdp2d, &
4934 h2osno2d, &
4935 snl2d, &
4936 t_grnd2d
4938 real, dimension( ims:ime,1:nlevlake, jms:jme ),INTENT(out) :: t_lake3d, &
4939 lake_icefrac3d, &
4940 z_lake3d, &
4941 dz_lake3d
4942 real, dimension( ims:ime,-nlevsnow+1:nlevsoil, jms:jme ),INTENT(out) :: t_soisno3d, &
4943 h2osoi_ice3d, &
4944 h2osoi_liq3d, &
4945 h2osoi_vol3d, &
4946 z3d, &
4947 dz3d
4948 real, dimension( ims:ime,1:nlevsoil, jms:jme ),INTENT(out) :: watsat3d, &
4949 csol3d, &
4950 tkmg3d, &
4951 tkdry3d, &
4952 tksatu3d
4953 real, dimension( ims:ime,-nlevsnow+0:nlevsoil, jms:jme ),INTENT(out) :: zi3d
4955 LOGICAL, DIMENSION( ims:ime, jms:jme ),intent(out) :: lake
4956 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: lake_depth
4958 real, dimension( ims:ime,1:nlevsoil, jms:jme ) :: bsw3d, &
4959 bsw23d, &
4960 psisat3d, &
4961 vwcsat3d, &
4962 watdry3d, &
4963 watopt3d, &
4964 hksat3d, &
4965 sucsat3d, &
4966 clay3d, &
4967 sand3d
4968 integer :: n,i,j,k,ib,lev,bottom ! indices
4969 real(r8),dimension(ims:ime,jms:jme ) :: bd2d ! bulk density of dry soil material [kg/m^3]
4970 real(r8),dimension(ims:ime,jms:jme ) :: tkm2d ! mineral conductivity
4971 real(r8),dimension(ims:ime,jms:jme ) :: xksat2d ! maximum hydraulic conductivity of soil [mm/s]
4972 real(r8),dimension(ims:ime,jms:jme ) :: depthratio2d ! ratio of lake depth to standard deep lake depth
4973 real(r8),dimension(ims:ime,jms:jme ) :: clay2d ! temporary
4974 real(r8),dimension(ims:ime,jms:jme ) :: sand2d ! temporary
4976 real(r8) :: scalez = 0.025_r8 ! Soil layer thickness discretization (m)
4977 logical,parameter :: arbinit = .true.
4978 real,parameter :: defval = -999.0
4979 integer :: isl
4980 integer :: numb_lak ! for debug
4981 character*256 :: message
4983 IF ( RESTART ) RETURN
4985 DO j = jts,jte
4986 DO i = its,ite
4987 snowdp2d(i,j) = snow(i,j)*0.005 ! SNOW in kg/m^2 and snowdp in m
4988 h2osno2d(i,j) = snow(i,j) ! mm
4989 ENDDO
4990 ENDDO
4992 ! initialize all the grid with default value
4993 DO j = jts,jte
4994 DO i = its,ite
4996 lakedepth2d(i,j) = defval
4997 snl2d(i,j) = defval
4998 do k = -nlevsnow+1,nlevsoil
4999 h2osoi_liq3d(i,k,j) = defval
5000 h2osoi_ice3d(i,k,j) = defval
5001 t_soisno3d(i,k,j) = defval
5002 z3d(i,k,j) = defval
5003 dz3d(i,k,j) = defval
5004 enddo
5005 do k = 1,nlevlake
5006 t_lake3d(i,k,j) = defval
5007 lake_icefrac3d(i,k,j) = defval
5008 z_lake3d(i,k,j) = defval
5009 dz_lake3d(i,k,j) = defval
5010 enddo
5012 ENDDO
5013 ENDDO
5015 ! judge whether the grid is lake grid
5016 numb_lak = 0
5017 do i=its,ite
5018 do j=jts,jte
5019 #if (EM_CORE==1)
5020 IF (lakeflag.eq.0) THEN
5021 if(ht(i,j)>=lake_min_elev) then
5022 if ( xice(i,j).gt.xice_threshold) then !mchen
5023 ivgtyp(i,j) = iswater
5024 xland(i,j) = 2.
5025 lake_icefrac3d(i,1,j) = xice(i,j)
5026 xice(i,j)=0.0
5027 endif
5028 endif
5030 if(ivgtyp(i,j)==iswater.and.ht(i,j)>=lake_min_elev) then
5031 lake(i,j) = .true.
5032 lakemask(i,j) = 1
5033 numb_lak = numb_lak + 1
5034 else
5035 lake(i,j) = .false.
5036 lakemask(i,j) = 0
5037 end if
5038 ELSE
5039 if(lakemask(i,j).eq.1) then
5040 lake(i,j) = .true.
5041 numb_lak = numb_lak + 1
5042 if ( xice(i,j).gt.xice_threshold) then !mchen
5043 ivgtyp(i,j) = iswater
5044 xland(i,j) = 2.
5045 lake_icefrac3d(i,1,j) = xice(i,j)
5046 xice(i,j)=0.0
5047 endif
5048 else
5049 lake(i,j) = .false.
5050 endif
5051 ENDIF ! end if lakeflag=0
5052 #else
5053 if(ht(i,j)>=lake_min_elev) then
5054 if ( xice(i,j).gt.xice_threshold) then !mchen
5055 ivgtyp(i,j) = iswater
5056 xland(i,j) = 2.
5057 lake_icefrac3d(i,1,j) = xice(i,j)
5058 xice(i,j)=0.0
5059 endif
5060 endif
5061 if(ivgtyp(i,j)==iswater.and.ht(i,j)>=lake_min_elev) then
5062 lake(i,j) = .true.
5063 numb_lak = numb_lak + 1
5064 else
5065 lake(i,j) = .false.
5066 end if
5068 #endif
5069 end do
5070 end do
5071 write(message,*) "the total number of lake grid is :", numb_lak
5072 CALL wrf_message(message)
5073 ! CALL LakeDebug(msg)
5074 ! initialize lake grid
5076 DO j = jts,jte
5077 DO i = its,ite
5079 if ( lake(i,j) ) then
5081 ! t_soisno3d(i,:,j) = tsk(i,j)
5082 ! t_lake3d(i,:,j) = tsk(i,j)
5083 ! t_grnd2d(i,j) = tsk(i,j)
5085 z3d(i,:,j) = 0.0
5086 dz3d(i,:,j) = 0.0
5087 zi3d(i,:,j) = 0.0
5088 h2osoi_liq3d(i,:,j) = 0.0
5089 h2osoi_ice3d(i,:,j) = 0.0
5090 lake_icefrac3d(i,:,j) = 0.0
5091 h2osoi_vol3d(i,:,j) = 0.0
5092 snl2d(i,j) = 0.0
5093 if ( use_lakedepth.eq.1 .and.lake_depth_flag.eq.0 ) then !mchen
5094 call wrf_error_fatal ( 'STOP: You need lake-depth information. Rerun WPS or set use_lakedepth = 0')
5095 end if
5096 if ( use_lakedepth.eq.0 .and.lake_depth_flag.eq.1 ) then !mchen
5097 lake_depth_flag = 0
5098 end if
5099 if ( lake_depth_flag.eq.1 ) then
5101 if (lake_depth(i,j) > 0.0) then
5102 lakedepth2d(i,j) = lake_depth(i,j)
5103 else
5104 if ( lakedepth_default > 0.0 ) then
5105 lakedepth2d(i,j) = lakedepth_default
5106 else
5107 lakedepth2d(i,j) = spval
5108 endif
5109 endif
5111 else
5112 if ( lakedepth_default > 0.0 ) then
5113 lakedepth2d(i,j) = lakedepth_default
5114 else
5115 lakedepth2d(i,j) = spval
5116 endif
5117 endif
5118 endif
5120 ENDDO
5121 ENDDO
5124 #ifndef EXTRALAKELAYERS
5125 ! dzlak(1) = 0.1_r8
5126 ! dzlak(2) = 1._r8
5127 ! dzlak(3) = 2._r8
5128 ! dzlak(4) = 3._r8
5129 ! dzlak(5) = 4._r8
5130 ! dzlak(6) = 5._r8
5131 ! dzlak(7) = 7._r8
5132 ! dzlak(8) = 7._r8
5133 ! dzlak(9) = 10.45_r8
5134 ! dzlak(10)= 10.45_r8
5135 !
5136 ! zlak(1) = 0.05_r8
5137 ! zlak(2) = 0.6_r8
5138 ! zlak(3) = 2.1_r8
5139 ! zlak(4) = 4.6_r8
5140 ! zlak(5) = 8.1_r8
5141 ! zlak(6) = 12.6_r8
5142 ! zlak(7) = 18.6_r8
5143 ! zlak(8) = 25.6_r8
5144 ! zlak(9) = 34.325_r8
5145 ! zlak(10)= 44.775_r8
5146 dzlak(1) = 0.1_r8
5147 dzlak(2) = 0.1_r8
5148 dzlak(3) = 0.1_r8
5149 dzlak(4) = 0.1_r8
5150 dzlak(5) = 0.1_r8
5151 dzlak(6) = 0.1_r8
5152 dzlak(7) = 0.1_r8
5153 dzlak(8) = 0.1_r8
5154 dzlak(9) = 0.1_r8
5155 dzlak(10)= 0.1_r8
5157 zlak(1) = 0.05_r8
5158 zlak(2) = 0.15_r8
5159 zlak(3) = 0.25_r8
5160 zlak(4) = 0.35_r8
5161 zlak(5) = 0.45_r8
5162 zlak(6) = 0.55_r8
5163 zlak(7) = 0.65_r8
5164 zlak(8) = 0.75_r8
5165 zlak(9) = 0.85_r8
5166 zlak(10)= 0.95_r8
5167 #else
5168 dzlak(1) =0.1_r8
5169 dzlak(2) =0.25_r8
5170 dzlak(3) =0.25_r8
5171 dzlak(4) =0.25_r8
5172 dzlak(5) =0.25_r8
5173 dzlak(6) =0.5_r8
5174 dzlak(7) =0.5_r8
5175 dzlak(8) =0.5_r8
5176 dzlak(9) =0.5_r8
5177 dzlak(10) =0.75_r8
5178 dzlak(11) =0.75_r8
5179 dzlak(12) =0.75_r8
5180 dzlak(13) =0.75_r8
5181 dzlak(14) =2_r8
5182 dzlak(15) =2_r8
5183 dzlak(16) =2.5_r8
5184 dzlak(17) =2.5_r8
5185 dzlak(18) =3.5_r8
5186 dzlak(19) =3.5_r8
5187 dzlak(20) =3.5_r8
5188 dzlak(21) =3.5_r8
5189 dzlak(22) =5.225_r8
5190 dzlak(23) =5.225_r8
5191 dzlak(24) =5.225_r8
5192 dzlak(25) =5.225_r8
5194 zlak(1) = dzlak(1)/2._r8
5195 do k = 2,nlevlake
5196 zlak(k) = zlak(k-1) + (dzlak(k-1)+dzlak(k))/2._r8
5197 end do
5198 #endif
5200 ! "0" refers to soil surface and "nlevsoil" refers to the bottom of model soil
5202 do j = 1, nlevsoil
5203 zsoi(j) = scalez*(exp(0.5_r8*(j-0.5_r8))-1._r8) !node depths
5204 enddo
5206 dzsoi(1) = 0.5_r8*(zsoi(1)+zsoi(2)) !thickness b/n two interfaces
5207 do j = 2,nlevsoil-1
5208 dzsoi(j)= 0.5_r8*(zsoi(j+1)-zsoi(j-1))
5209 enddo
5210 dzsoi(nlevsoil) = zsoi(nlevsoil)-zsoi(nlevsoil-1)
5212 zisoi(0) = 0._r8
5213 do j = 1, nlevsoil-1
5214 zisoi(j) = 0.5_r8*(zsoi(j)+zsoi(j+1)) !interface depths
5215 enddo
5216 zisoi(nlevsoil) = zsoi(nlevsoil) + 0.5_r8*dzsoi(nlevsoil)
5219 !!!!!!!!!!!!!!!!!!begin to initialize lake variables!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5221 DO j = jts,jte
5222 DO i = its,ite
5224 if ( lake(i,j) ) then
5226 ! Soil hydraulic and thermal properties
5227 isl = ISLTYP(i,j)
5228 if (isl == 14 ) isl = isl + 1
5229 do k = 1,nlevsoil
5230 sand3d(i,k,j) = sand(isl)
5231 clay3d(i,k,j) = clay(isl)
5232 enddo
5234 do k = 1,nlevsoil
5235 clay2d(i,j) = clay3d(i,k,j)
5236 sand2d(i,j) = sand3d(i,k,j)
5237 watsat3d(i,k,j) = 0.489_r8 - 0.00126_r8*sand2d(i,j)
5238 bd2d(i,j) = (1._r8-watsat3d(i,k,j))*2.7e3_r8
5239 xksat2d(i,j) = 0.0070556_r8 *( 10._r8**(-0.884_r8+0.0153_r8*sand2d(i,j)) ) ! mm/s
5240 tkm2d(i,j) = (8.80_r8*sand2d(i,j)+2.92_r8*clay2d(i,j))/(sand2d(i,j)+clay2d(i,j)) ! W/(m K)
5242 bsw3d(i,k,j) = 2.91_r8 + 0.159_r8*clay2d(i,j)
5243 bsw23d(i,k,j) = -(3.10_r8 + 0.157_r8*clay2d(i,j) - 0.003_r8*sand2d(i,j))
5244 psisat3d(i,k,j) = -(exp((1.54_r8 - 0.0095_r8*sand2d(i,j) + 0.0063_r8*(100.0_r8-sand2d(i,j) &
5245 -clay2d(i,j)))*log(10.0_r8))*9.8e-5_r8)
5246 vwcsat3d(i,k,j) = (50.5_r8 - 0.142_r8*sand2d(i,j) - 0.037_r8*clay2d(i,j))/100.0_r8
5247 hksat3d(i,k,j) = xksat2d(i,j)
5248 sucsat3d(i,k,j) = 10._r8 * ( 10._r8**(1.88_r8-0.0131_r8*sand2d(i,j)) )
5249 tkmg3d(i,k,j) = tkm2d(i,j) ** (1._r8- watsat3d(i,k,j))
5250 tksatu3d(i,k,j) = tkmg3d(i,k,j)*0.57_r8**watsat3d(i,k,j)
5251 tkdry3d(i,k,j) = (0.135_r8*bd2d(i,j) + 64.7_r8) / (2.7e3_r8 - 0.947_r8*bd2d(i,j))
5252 csol3d(i,k,j) = (2.128_r8*sand2d(i,j)+2.385_r8*clay2d(i,j)) / (sand2d(i,j)+clay2d(i,j))*1.e6_r8 ! J/(m3 K)
5253 watdry3d(i,k,j) = watsat3d(i,k,j) * (316230._r8/sucsat3d(i,k,j)) ** (-1._r8/bsw3d(i,k,j))
5254 watopt3d(i,k,j) = watsat3d(i,k,j) * (158490._r8/sucsat3d(i,k,j)) ** (-1._r8/bsw3d(i,k,j))
5255 end do
5256 if (lakedepth2d(i,j) == spval) then
5257 lakedepth2d(i,j) = zlak(nlevlake) + 0.5_r8*dzlak(nlevlake)
5258 z_lake3d(i,1:nlevlake,j) = zlak(1:nlevlake)
5259 dz_lake3d(i,1:nlevlake,j) = dzlak(1:nlevlake)
5260 else
5261 depthratio2d(i,j) = lakedepth2d(i,j) / (zlak(nlevlake) + 0.5_r8*dzlak(nlevlake))
5262 z_lake3d(i,1,j) = zlak(1)
5263 dz_lake3d(i,1,j) = dzlak(1)
5264 dz_lake3d(i,2:nlevlake,j) = dzlak(2:nlevlake)*depthratio2d(i,j)
5265 z_lake3d(i,2:nlevlake,j) = zlak(2:nlevlake)*depthratio2d(i,j) + dz_lake3d(i,1,j)*(1._r8 - depthratio2d(i,j))
5266 end if
5267 ! initial t_lake3d here
5268 t_soisno3d(i,1,j) = tsk(i,j)
5269 t_lake3d(i,1,j) = tsk(i,j)
5270 t_grnd2d(i,j) = 277.0
5271 do k = 2, nlevlake
5272 if(z_lake3d(i,k,j).le.depth_c) then
5273 t_soisno3d(i,k,j)=tsk(i,j)+(277.0-tsk(i,j))/depth_c*z_lake3d(i,k,j)
5274 t_lake3d(i,k,j)=tsk(i,j)+(277.0-tsk(i,j))/depth_c*z_lake3d(i,k,j)
5275 else
5276 t_soisno3d(i,k,j) = 277.0
5277 t_lake3d(i,k,j) = 277.0
5278 end if
5279 enddo
5280 !end initial t_lake3d here
5281 z3d(i,1:nlevsoil,j) = zsoi(1:nlevsoil)
5282 zi3d(i,0:nlevsoil,j) = zisoi(0:nlevsoil)
5283 dz3d(i,1:nlevsoil,j) = dzsoi(1:nlevsoil)
5284 savedtke12d(i,j) = tkwat ! Initialize for first timestep.
5287 if (snowdp2d(i,j) < 0.01_r8) then
5288 snl2d(i,j) = 0
5289 dz3d(i,-nlevsnow+1:0,j) = 0._r8
5290 z3d (i,-nlevsnow+1:0,j) = 0._r8
5291 zi3d(i,-nlevsnow+0:0,j) = 0._r8
5292 else
5293 if ((snowdp2d(i,j) >= 0.01_r8) .and. (snowdp2d(i,j) <= 0.03_r8)) then
5294 snl2d(i,j) = -1
5295 dz3d(i,0,j) = snowdp2d(i,j)
5296 else if ((snowdp2d(i,j) > 0.03_r8) .and. (snowdp2d(i,j) <= 0.04_r8)) then
5297 snl2d(i,j) = -2
5298 dz3d(i,-1,j) = snowdp2d(i,j)/2._r8
5299 dz3d(i, 0,j) = dz3d(i,-1,j)
5300 else if ((snowdp2d(i,j) > 0.04_r8) .and. (snowdp2d(i,j) <= 0.07_r8)) then
5301 snl2d(i,j) = -2
5302 dz3d(i,-1,j) = 0.02_r8
5303 dz3d(i, 0,j) = snowdp2d(i,j) - dz3d(i,-1,j)
5304 else if ((snowdp2d(i,j) > 0.07_r8) .and. (snowdp2d(i,j) <= 0.12_r8)) then
5305 snl2d(i,j) = -3
5306 dz3d(i,-2,j) = 0.02_r8
5307 dz3d(i,-1,j) = (snowdp2d(i,j) - 0.02_r8)/2._r8
5308 dz3d(i, 0,j) = dz3d(i,-1,j)
5309 else if ((snowdp2d(i,j) > 0.12_r8) .and. (snowdp2d(i,j) <= 0.18_r8)) then
5310 snl2d(i,j) = -3
5311 dz3d(i,-2,j) = 0.02_r8
5312 dz3d(i,-1,j) = 0.05_r8
5313 dz3d(i, 0,j) = snowdp2d(i,j) - dz3d(i,-2,j) - dz3d(i,-1,j)
5314 else if ((snowdp2d(i,j) > 0.18_r8) .and. (snowdp2d(i,j) <= 0.29_r8)) then
5315 snl2d(i,j) = -4
5316 dz3d(i,-3,j) = 0.02_r8
5317 dz3d(i,-2,j) = 0.05_r8
5318 dz3d(i,-1,j) = (snowdp2d(i,j) - dz3d(i,-3,j) - dz3d(i,-2,j))/2._r8
5319 dz3d(i, 0,j) = dz3d(i,-1,j)
5320 else if ((snowdp2d(i,j) > 0.29_r8) .and. (snowdp2d(i,j) <= 0.41_r8)) then
5321 snl2d(i,j) = -4
5322 dz3d(i,-3,j) = 0.02_r8
5323 dz3d(i,-2,j) = 0.05_r8
5324 dz3d(i,-1,j) = 0.11_r8
5325 dz3d(i, 0,j) = snowdp2d(i,j) - dz3d(i,-3,j) - dz3d(i,-2,j) - dz3d(i,-1,j)
5326 else if ((snowdp2d(i,j) > 0.41_r8) .and. (snowdp2d(i,j) <= 0.64_r8)) then
5327 snl2d(i,j) = -5
5328 dz3d(i,-4,j) = 0.02_r8
5329 dz3d(i,-3,j) = 0.05_r8
5330 dz3d(i,-2,j) = 0.11_r8
5331 dz3d(i,-1,j) = (snowdp2d(i,j) - dz3d(i,-4,j) - dz3d(i,-3,j) - dz3d(i,-2,j))/2._r8
5332 dz3d(i, 0,j) = dz3d(i,-1,j)
5333 else if (snowdp2d(i,j) > 0.64_r8) then
5334 snl2d(i,j) = -5
5335 dz3d(i,-4,j) = 0.02_r8
5336 dz3d(i,-3,j) = 0.05_r8
5337 dz3d(i,-2,j) = 0.11_r8
5338 dz3d(i,-1,j) = 0.23_r8
5339 dz3d(i, 0,j)=snowdp2d(i,j)-dz3d(i,-4,j)-dz3d(i,-3,j)-dz3d(i,-2,j)-dz3d(i,-1,j)
5340 endif
5341 end if
5343 do k = 0, snl2d(i,j)+1, -1
5344 z3d(i,k,j) = zi3d(i,k,j) - 0.5_r8*dz3d(i,k,j)
5345 zi3d(i,k-1,j) = zi3d(i,k,j) - dz3d(i,k,j)
5346 end do
5348 ! 3:subroutine makearbinit
5350 if (snl2d(i,j) < 0) then
5351 do k = snl2d(i,j)+1, 0
5352 ! Be careful because there may be new snow layers with bad temperatures like 0 even if
5353 ! coming from init. con. file.
5354 if(arbinit .or. t_soisno3d(i,k,j) > 300 .or. t_soisno3d(i,k,j) < 200) t_soisno3d(i,k,j) = 250._r8
5355 enddo
5356 end if
5358 do k = 1, nlevsoil
5359 if(arbinit .or. t_soisno3d(i,k,j) > 1000 .or. t_soisno3d(i,k,j) < 0) t_soisno3d(i,k,j) = t_lake3d(i,nlevlake,j)
5360 end do
5362 do k = 1, nlevlake
5363 if(arbinit .or. lake_icefrac3d(i,k,j) > 1._r8 .or. lake_icefrac3d(i,k,j) < 0._r8) then
5364 if(t_lake3d(i,k,j) >= tfrz) then
5365 lake_icefrac3d(i,k,j) = 0._r8
5366 else
5367 lake_icefrac3d(i,k,j) = 1._r8
5368 end if
5369 end if
5370 end do
5372 do k = 1,nlevsoil
5373 if (arbinit .or. h2osoi_vol3d(i,k,j) > 10._r8 .or. h2osoi_vol3d(i,k,j) < 0._r8) h2osoi_vol3d(i,k,j) = 1.0_r8
5374 h2osoi_vol3d(i,k,j) = min(h2osoi_vol3d(i,k,j),watsat3d(i,k,j))
5376 ! soil layers
5377 if (t_soisno3d(i,k,j) <= tfrz) then
5378 h2osoi_ice3d(i,k,j) = dz3d(i,k,j)*denice*h2osoi_vol3d(i,k,j)
5379 h2osoi_liq3d(i,k,j) = 0._r8
5380 else
5381 h2osoi_ice3d(i,k,j) = 0._r8
5382 h2osoi_liq3d(i,k,j) = dz3d(i,k,j)*denh2o*h2osoi_vol3d(i,k,j)
5383 endif
5384 enddo
5386 do k = -nlevsnow+1, 0
5387 if (k > snl2d(i,j)) then
5388 h2osoi_ice3d(i,k,j) = dz3d(i,k,j)*bdsno
5389 h2osoi_liq3d(i,k,j) = 0._r8
5390 end if
5391 end do
5393 end if !lake(i,j)
5394 ENDDO
5395 ENDDO
5397 END SUBROUTINE lakeini
5399 END MODULE module_sf_lake