CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
clm_lake.f90
1
3
22MODULE clm_lake
23
24 use machine, only: kind_phys, kind_dbl_prec
25
26 implicit none
27
28 private
29
30 public :: clm_lake_run, clm_lake_init, lakedebug
31
32 ! In WRF, the CLM Lake Model was hard-coded to use double precision, regardless of
33 ! precision of physics. For that reason, we retain double precision here. However,
34 ! we're not yet certain that all of CLM Lake needs to be double precision, so we
35 ! maintain a "kind_lake" to allow future experimentation in datatypes.
36 integer, parameter, public :: kind_lake = kind_dbl_prec
37
38 logical :: lakedebug = .false. ! Enable lots of checks and debug prints and errors
39 logical :: debug_print = .false. ! Enable lots of checks and debug prints and errors
40
41 logical, parameter :: pergro = .false.
42
43 logical, parameter :: use_etalake = .false.
44 real(kind_lake), parameter :: etalake = 1.1925*50**(-0.424)
45
46 ! Level counts must be consistent with model (GFS_typedefs.F90)
47 integer, parameter :: nlevsoil = 10
48 integer, parameter :: nlevlake = 10
49 integer, parameter :: nlevsnow = 5
50 real(kind_lake), parameter :: scalez = 0.025_kind_lake
51
52 integer,parameter :: lbp = 1
53 integer,parameter :: ubp = 1
54 integer,parameter :: lbc = 1
55 integer,parameter :: ubc = 1
56 integer,parameter :: num_shlakec = 1
57 integer,parameter :: filter_shlakec(1) = 1
58 integer,parameter :: num_shlakep = 1
59 integer,parameter :: filter_shlakep(1) = 1
60 integer,parameter :: pcolumn(1) = 1
61 integer,parameter :: pgridcell(1) = 1
62 integer,parameter :: cgridcell(1) = 1
63 integer,parameter :: clandunit(1) = 1
64
65 integer,parameter :: begg = 1
66 integer,parameter :: endg = 1
67 integer,parameter :: begl = 1
68 integer,parameter :: endl = 1
69 integer,parameter :: begc = 1
70 integer,parameter :: endc = 1
71 integer,parameter :: begp = 1
72 integer,parameter :: endp = 1
73
74 integer,parameter :: column =1
75 logical,parameter :: lakpoi(1) = .true.
76
77 !Initialize physical constants not available from model:
78 real(kind_lake), parameter :: tcrit = 2.5
79 real(kind_lake), parameter :: tkwat = 0.6
80 real(kind_lake), parameter :: tkice = 2.290
81 real(kind_lake), parameter :: tkairc = 0.023
82 real(kind_lake), parameter :: snow_bd = 250
83
84 ! Constants that are copied from model values by clm_lake_init:
85 real(kind_lake) :: pi
86 real(kind_lake) :: vkc
87 real(kind_lake) :: grav
88 real(kind_lake) :: sb
89 real(kind_lake) :: tfrz
90 real(kind_lake) :: denh2o
91 real(kind_lake) :: denice
92 real(kind_lake) :: cpice
93 real(kind_lake) :: cpliq
94 real(kind_lake) :: hfus
95 real(kind_lake) :: hvap
96 real(kind_lake) :: hsub
97 real(kind_lake) :: invhvap
98 real(kind_lake) :: invhsub
99 real(kind_lake) :: rair
100 real(kind_lake) :: cpair
101 real(kind_lake) :: con_eps
102 real(kind_lake) :: one_minus_con_eps
103 real(kind_lake) :: con_fvirt
104
105 real(kind_lake), public, parameter :: spval = 1.e36
106 real(kind_lake), parameter :: depth_c = 50.
107 real(kind_lake), parameter :: zero_h2o = 1e-12
108
109 ! These are tunable constants
110 real(kind_lake), parameter :: wimp = 0.05
111 real(kind_lake), parameter :: ssi = 0.033
112 real(kind_lake), parameter :: cnfac = 0.5
113
114 ! Initialize water type constants
115 integer,parameter :: istsoil = 1!soil "water" type
116
117 ! percent sand
118 real(kind_lake), parameter :: sand(19) = &
119 (/92.,80.,66.,20.,5.,43.,60.,10.,32.,51., 6.,22.,39.7,0.,100.,54.,17.,100.,92./)
120
121 ! percent clay
122 real(kind_lake), parameter :: clay(19) = &
123 (/ 3., 5.,10.,15.,5.,18.,27.,33.,33.,41.,47.,58.,14.7,0., 0., 8.5,54., 0., 3./)
124
125 ! These are initialized in clm_lake_init and are not modified elsewhere
126 real(kind_lake) :: zlak(1:nlevlake)
127 real(kind_lake) :: dzlak(1:nlevlake)
128 real(kind_lake) :: zsoi(1:nlevsoil)
129 real(kind_lake) :: dzsoi(1:nlevsoil)
130 real(kind_lake) :: zisoi(0:nlevsoil)
131
132 real, parameter :: saltlk_t(1:25) = (/ 0.5, 0.,-0.5, 3., 4., 7., 8., 12., 13., 16., 19., 21., &
133 23.5, 25., 26.,24.,23.,20.5,18., 15., 11.5, 8., 4., 1., 0.5/)
134 real, parameter :: month_length(12) = (/ 31, 29, 31, 30, 31, 30, 31, 30, 30, 31, 30, 31 /)
135 logical, parameter :: include_all_salty_locations = .false.
136
137 CONTAINS
138
139 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
140
141 logical function limit_temperature_by_climatology(xlat_d,xlon_positive)
142 implicit none
143 real(kind_phys), intent(in) :: xlat_d, xlon_positive
144 real(kind_phys) :: xlon_d
145
146 xlon_d = xlon_positive
147 if(xlon_d>180) xlon_d = xlon_d - 360
148
149 limit_temperature_by_climatology=.false.
150
160 if ((xlon_d.gt.-117.7 .and. xlon_d.lt.-111.5) .and. &
161 ! excludes Willard Bay
162 .not. (xlon_d.gt.-112.104 .and. xlon_d.lt.-112.100))then
163
164 if(xlat_d.gt.39.5 .and. xlat_d.lt.41.22) then
165 if(debug_print) then
166 print *,'The Great Salt Lake south of 41.22 N, lat,lon',xlat_d,xlon_d
167 endif
168 limit_temperature_by_climatology = .true.
169
170 elseif(( xlat_d.ge.41.22 .and. xlat_d.lt.42.) .and. .not. &
171 ! excludes Willard Bay
172 (xlat_d.gt.41.352 .and. xlat_d.lt.41.354)) then
173 if(debug_print) then
174 print *,'The Great Salt Lake north of 41.22 N xlat_d,xlon_d ',xlat_d,xlon_d
175 endif
176 !print *,'Ice fraction on the GSL ', i,j,lake_icefrac3d(i,:,j)
177 limit_temperature_by_climatology = .true.
178
179 endif ! xlat_d
180
181 endif ! xlon_d
182
183 !if(i==495.and.j==668) print *,'Willard Bay salty=',i,j,limit_temperature_by_climatology,xlat_d,xlon_d
184
185 end function limit_temperature_by_climatology
186
187 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
188
189 subroutine is_salty(xlat_d,xlon_positive, cannot_freeze, salty)
190 implicit none
191 real(kind_phys), intent(in) :: xlat_d, xlon_positive
192 logical, intent(inout) :: cannot_freeze, salty
193 real(kind_phys) :: xlon_d
194
195 xlon_d = xlon_positive
196 if(xlon_d>180) xlon_d = xlon_d - 360
197
198 ! for the Great Salt Lake
199 cannot_freeze = limit_temperature_by_climatology(xlat_d,xlon_d)
200 salty = cannot_freeze
201
202 ! --- The Mono Lake in California, salinity is 75 ppt with freezing point at
203 ! --- -4.2 C (Stan). The Mono Lake lat/long (37.9-38.2, -119.3 - 118.8)
204 if (xlon_d.gt.-119.3.and. xlon_d.lt.-118.8) then
205 if(xlat_d.gt.37.9 .and. xlat_d.lt.38.2) then
206 salty = .true.
207 if(debug_print) then
208 print *,'Salty Mono Lake, i,j',xlat_d,xlon_d
209 endif
210 endif ! xlat_d
211 endif ! xlon_d
212
213 other_locations: if(include_all_salty_locations) then
214 ! --- Caspian Sea and Dead Sea are salty too (Sam, Tanya)
215 if ( xlat_d>36.5_kind_phys .and. xlat_d<47.1_kind_phys .and. xlon_d>46.8_kind_phys .and. xlon_d<55.0_kind_phys ) then
216 if(debug_print) then
217 print *,'Salty Caspian Sea ',xlat_d,xlon_d
218 endif
219 salty = .true.
220 end if
221 if ( xlon_d>35.3 .and. xlon_d<35.6 .and. xlat_d>31.3 .and. xlat_d<31.8) then
222 if(debug_print) then
223 print *,'Salty Dead Sea ',xlat_d,xlon_d
224 endif
225 salty = .true.
226 endif
227 endif other_locations
228 !tgs --- end of special treatment for salty lakes
229 end subroutine is_salty
230
231 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
232
233 subroutine calculate_z_dz_lake(i,input_lakedepth,clm_lakedepth,z_lake,dz_lake)
234 implicit none
235 integer, intent(in) :: i
236 real(kind_phys), intent(inout) :: clm_lakedepth(:) ! lake depth used by clm
237 real(kind_phys), intent(in) :: input_lakedepth(:) ! lake depth before correction (m)
238 real(kind_lake) :: z_lake(nlevlake) ! layer depth for lake (m)
239 real(kind_lake) :: dz_lake(nlevlake) ! layer thickness for lake (m)
240 real(kind_lake) :: depthratio
241
242 if (input_lakedepth(i) == spval .or. input_lakedepth(i) < 0.1) then
243 ! This is a safeguard against:
244 ! 1. missing in the lakedepth database (== spval)
245 ! 2. errors in model cycling or unexpected changes in the orography database (< 0.1)
246 clm_lakedepth(i) = zlak(nlevlake) + 0.5_kind_lake*dzlak(nlevlake)
247 z_lake(1:nlevlake) = zlak(1:nlevlake)
248 dz_lake(1:nlevlake) = dzlak(1:nlevlake)
249 else
250 depthratio = input_lakedepth(i) / (zlak(nlevlake) + 0.5_kind_lake*dzlak(nlevlake))
251 z_lake(1) = zlak(1)
252 dz_lake(1) = dzlak(1)
253 dz_lake(2:nlevlake) = dzlak(2:nlevlake)*depthratio
254 z_lake(2:nlevlake) = zlak(2:nlevlake)*depthratio + dz_lake(1)*(1._kind_lake - depthratio)
255 end if
256
257 end subroutine calculate_z_dz_lake
258
259 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
260
264 SUBROUTINE clm_lake_run( &
265 ! Model time and metadata:
266 flag_restart, im, km, me, master, fhour, IDATE, kdt, &
267
268 ! Configuration and initialization:
269 iopt_lake, iopt_lake_clm, min_lakeice, lakedepth_default, use_lakedepth, &
270 dtp, use_lake_model, clm_lake_initialized, frac_grid, frac_ice, lkm, &
271
272 ! Atmospheric model state inputs:
273 tg3, pgr, zlvl, gt0, prsi, phii, qvcurr, gu0, gv0, xlat_d, xlon_d, &
274 ch, cm, dlwsfci, dswsfci, oro_lakedepth, wind, tsfc, &
275 flag_iter, flag_lakefreeze, ISLTYP, rainncprv, raincprv, &
276
277 ! Feedback to atmosphere:
278 evap_wat, evap_ice, hflx_wat, hflx_ice, gflx_wat, gflx_ice, &
279 ep1d_water, ep1d_ice, tsurf_water, tsurf_ice, tsfc_wat, tisfc, &
280 weasdi, snodi, hice, qss_water, qss_ice, &
281 cmm_water, cmm_ice, chh_water, chh_ice, &
282 uustar_water, uustar_ice, lake_t_snow, albedo, zorlw, &
283 zorli, lake_t2m, lake_q2m, weasd, snowd, fice, &
284 icy, &
285
286 ! Lake model internal state stored by caller:
287
288 salty, savedtke12d, snowdp2d, h2osno2d, snl2d, t_grnd2d, t_lake3d, &
289 lake_icefrac3d, t_soisno3d, h2osoi_ice3d, h2osoi_liq3d, h2osoi_vol3d, &
290 z3d, dz3d, zi3d, t1, qv1, prsl1, &
291 input_lakedepth, clm_lakedepth, cannot_freeze, &
292
293 ! Error reporting:
294 errflg, errmsg)
295
296 !==============================================================================
297 ! This subroutine was first edited by Hongping Gu and Jiming Jin for coupling
298 ! 07/20/2010
299 ! Long after, in June 2022, Sam Trahan updated it for CCPP
300 !==============================================================================
301
302 IMPLICIT NONE
303
304 !
305 ! Model time and metadata:
306 !
307 LOGICAL , INTENT (IN) :: flag_restart
308 INTEGER , INTENT (IN) :: im,km,me,master
309 INTEGER, INTENT(IN) :: idate(4), kdt
310 REAL(kind_phys), INTENT(IN) :: fhour
311 INTEGER, INTENT(IN) :: lkm
312
313 !
314 ! Configuration and initialization:
315 !
316 INTEGER, INTENT(IN) :: iopt_lake, iopt_lake_clm
317 REAL(kind_phys), INTENT(IN) :: min_lakeice, lakedepth_default, dtp
318 LOGICAL, INTENT(IN) :: use_lakedepth
319 INTEGER, DIMENSION(:), INTENT(IN) :: use_lake_model
320 REAL(kind_phys), INTENT(INOUT), OPTIONAL :: clm_lake_initialized(:)
321 LOGICAL, INTENT(IN) :: frac_grid, frac_ice
322
323 !
324 ! Atmospheric model state inputs:
325 !
326 REAL(kind_phys), DIMENSION(:), INTENT(IN):: &
327 tg3, pgr, zlvl, qvcurr, xlat_d, xlon_d, ch, cm, &
328 dlwsfci, dswsfci, oro_lakedepth, wind, &
329 t1, qv1, prsl1
330 REAL(kind_phys), DIMENSION(:), INTENT(IN), OPTIONAL :: &
331 rainncprv, raincprv
332 REAL(kind_phys), DIMENSION(:,:), INTENT(in) :: gu0, gv0, prsi, gt0, phii
333 LOGICAL, DIMENSION(:), INTENT(IN) :: flag_iter
334 LOGICAL, DIMENSION(:), INTENT(INOUT) :: flag_lakefreeze
335
336 INTEGER, DIMENSION(:), INTENT(IN) :: isltyp
337
338 !
339 ! Feedback to atmosphere:
340 !
341 REAL(kind_phys), DIMENSION(:), INTENT(INOUT) :: &
342 evap_wat, evap_ice, hflx_wat, hflx_ice, gflx_wat, gflx_ice, &
343 ep1d_water, ep1d_ice, tsurf_water, tsurf_ice, tsfc_wat, tisfc, tsfc, &
344 weasdi, snodi, hice, qss_water, qss_ice, &
345 cmm_water, cmm_ice, chh_water, chh_ice, &
346 uustar_water, uustar_ice, zorlw, zorli, weasd, snowd, fice
347 REAL(kind_phys), DIMENSION(:), INTENT(INOUT) , OPTIONAL :: &
348 lake_t_snow, albedo, lake_t2m, lake_q2m
349 LOGICAL, INTENT(INOUT) :: icy(:)
350
351 !
352 ! Lake model internal state stored by caller:
353 !
354 INTEGER, DIMENSION( : ), INTENT(INOUT), OPTIONAL :: salty
355 INTEGER, DIMENSION( : ), INTENT(INOUT), OPTIONAL :: cannot_freeze
356
357 real(kind_phys), dimension(: ), OPTIONAL ,intent(inout) :: savedtke12d, &
358 snowdp2d, &
359 h2osno2d, &
360 snl2d, &
361 t_grnd2d
362
363 real(kind_phys), dimension( :,: ), OPTIONAL, INTENT(inout) :: t_lake3d, &
364 lake_icefrac3d
365 real(kind_phys), dimension( :,-nlevsnow+1: ) ,INTENT(inout), OPTIONAL :: t_soisno3d, &
366 h2osoi_ice3d, &
367 h2osoi_liq3d, &
368 h2osoi_vol3d, &
369 z3d, &
370 dz3d
371 real(kind_phys), dimension( :,-nlevsnow+0: ) ,INTENT(inout), OPTIONAL :: zi3d
372
373 REAL(kind_phys), DIMENSION( : ) ,INTENT(INOUT), OPTIONAL :: clm_lakedepth
374 REAL(kind_phys), DIMENSION( : ) ,INTENT(INOUT), OPTIONAL :: input_lakedepth
375
376 !
377 ! Error reporting:
378 !
379 INTEGER, INTENT(OUT) :: errflg
380 CHARACTER(*), INTENT(OUT) :: errmsg
381
382
383
384 !
385 !local variables:
386 !
387
388 REAL(kind_lake) :: sfctmp,pbot,psfc,q2k,lwdn,prcp,soldn,solnet,dtime
389 INTEGER :: c,i,j,k
390
391
392 !temporary varibles in:
393 real(kind_lake) :: forc_t(1) ! atmospheric temperature (Kelvin)
394 real(kind_lake) :: forc_pbot(1) ! atm bottom level pressure (Pa)
395 real(kind_lake) :: forc_psrf(1) ! atmospheric surface pressure (Pa)
396 real(kind_lake) :: forc_hgt(1) ! atmospheric reference height (m)
397 real(kind_lake) :: forc_hgt_q(1) ! observational height of humidity [m]
398 real(kind_lake) :: forc_hgt_t(1) ! observational height of temperature [m]
399 real(kind_lake) :: forc_hgt_u(1) ! observational height of wind [m]
400 real(kind_lake) :: forc_q(1) ! atmospheric specific humidity (kg/kg)
401 real(kind_lake) :: forc_u(1) ! atmospheric wind speed in east direction (m/s)
402 real(kind_lake) :: forc_v(1) ! atmospheric wind speed in north direction (m/s)
403 real(kind_lake) :: forc_lwrad(1) ! downward infrared (longwave) radiation (W/m**2)
404 real(kind_lake) :: prec(1) ! snow or rain rate [mm/s]
405 real(kind_lake) :: sabg(1) ! solar radiation absorbed by ground (W/m**2)
406 real(kind_lake) :: lat(1) ! latitude (radians)
407 real(kind_lake) :: z_lake(1,nlevlake) ! layer depth for lake (m)
408 real(kind_lake) :: dz_lake(1,nlevlake) ! layer thickness for lake (m)
409
410 real(kind_lake) :: lakedepth(1) ! column lake depth (m)
411 logical :: do_capsnow(1) ! true => do snow capping
412
413 !in&out
414 real(kind_lake) :: h2osoi_vol(1,-nlevsnow+1:nlevsoil) ! volumetric soil water (0<=h2osoi_vol<=watsat)[m3/m3]
415 real(kind_lake) :: t_grnd(1) ! ground temperature (Kelvin)
416 real(kind_lake) :: h2osno(1) ! snow water (mm H2O)
417 real(kind_lake) :: snowdp(1) ! snow height (m)
418 real(kind_lake) :: z(1,-nlevsnow+1:nlevsoil) ! layer depth for snow & soil (m)
419 real(kind_lake) :: dz(1,-nlevsnow+1:nlevsoil) ! layer thickness for soil or snow (m)
420 real(kind_lake) :: t_soisno(1,-nlevsnow+1:nlevsoil) ! soil (or snow) temperature (Kelvin)
421 real(kind_lake) :: t_lake(1,nlevlake) ! lake temperature (Kelvin)
422 integer :: snl(1) ! number of snow layers
423 real(kind_lake) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil) ! liquid water (kg/m2)
424 real(kind_lake) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil) ! ice lens (kg/m2)
425 real(kind_lake) :: savedtke1(1) ! top level eddy conductivity from previous timestep (W/m.K)
426 real(kind_lake) :: zi(1,-nlevsnow+0:nlevsoil) ! interface level below a "z" level (m)
427 real(kind_lake) :: lake_icefrac(1,nlevlake) ! mass fraction of lake layer that is frozen
428
429
430 !out:
431 real(kind_lake) :: eflx_gnet(1) !net heat flux into ground (W/m**2)
432 real(kind_lake) :: eflx_lwrad_net(1) ! net infrared (longwave) rad (W/m**2) [+ = to atm]
433 real(kind_lake) :: eflx_sh_tot(1) ! total sensible heat flux (W/m**2) [+ to atm]
434 real(kind_lake) :: eflx_lh_tot(1) ! total latent heat flux (W/m8*2) [+ to atm]
435 real(kind_lake) :: t_ref2m(1) ! 2 m height surface air temperature (Kelvin)
436 real(kind_lake) :: q_ref2m(1) ! 2 m height surface specific humidity (kg/kg)
437 real(kind_lake) :: taux(1) ! wind (shear) stress: e-w (kg/m/s**2)
438 real(kind_lake) :: tauy(1) ! wind (shear) stress: n-s (kg/m/s**2)
439 real(kind_lake) :: ram1(1) ! aerodynamical resistance (s/m)
440 ! for calculation of decay of eddy diffusivity with depth
441 ! Change the type variable to pass back to WRF.
442 real(kind_lake) :: z0mg(1) ! roughness length over ground, momentum (m(
443 real(kind_lake) :: qfx ! mass flux, old WRF qfx(:) variable, (kg/(sm^2))
444
445 real(kind_lake) :: ustar_out(1) ! friction velocity (temporary) [m/s]
446
447 real(kind_lake) :: discard1, discard2, discard3 ! for unused temporary data
448
449 real(kind_lake) :: watsat(1,nlevsoil) ! volumetric soil water at saturation (porosity)
450 real(kind_lake) :: tksatu(1,nlevsoil) ! thermal conductivity, saturated soil [W/m-K]
451 real(kind_lake) :: tkmg(1,nlevsoil) ! thermal conductivity, soil minerals [W/m-K]
452 real(kind_lake) :: tkdry(1,nlevsoil) ! thermal conductivity, dry soil (W/m/Kelvin)
453 real(kind_lake) :: csol(1,nlevsoil) ! heat capacity, soil solids (J/m**3/Kelvin)
454
455! real(kind_lake) :: emiss ! surface emissivity
456
457 integer :: lake_points, snow_points, ice_points
458 character*255 :: message
459 logical, parameter :: feedback_to_atmosphere = .true. ! FIXME: REMOVE
460
461 real(kind_lake) :: to_radians, lat_d, lon_d, qss, tkm, bd
462 real(kind_lake) :: rho0 ! lowest model level air density
463
464 integer :: month,num1,num2,day_of_month,isl
465 real(kind_lake) :: wght1,wght2,tclim,depthratio
466
467 logical salty_flag, cannot_freeze_flag
468
469 errmsg = ' '
470 errflg = 0
471
472 if(iopt_lake/=iopt_lake_clm .or. lkm==0 .or. all(use_lake_model==0)) then
473 return ! nothing to do
474 endif
475
476 dtime=dtp
477
478 ! Initialize any uninitialized lake points.
479 call lakeini(kdt=kdt, isltyp=isltyp, gt0=gt0, snowd=snowd, weasd=weasd, &
480 lakedepth_default=lakedepth_default, fhour=fhour, &
481 oro_lakedepth=oro_lakedepth, savedtke12d=savedtke12d, snowdp2d=snowdp2d, &
482 h2osno2d=h2osno2d, snl2d=snl2d, t_grnd2d=t_grnd2d, t_lake3d=t_lake3d, &
483 lake_icefrac3d=lake_icefrac3d, &
484 t_soisno3d=t_soisno3d, h2osoi_ice3d=h2osoi_ice3d, h2osoi_liq3d=h2osoi_liq3d, &
485 h2osoi_vol3d=h2osoi_vol3d, z3d=z3d, dz3d=dz3d, zi3d=zi3d, &
486 fice=fice, hice=hice, min_lakeice=min_lakeice, &
487 tsfc=tsfc, &
488 use_lake_model=use_lake_model, use_lakedepth=use_lakedepth, &
489 im=im, prsi=prsi, xlat_d=xlat_d, xlon_d=xlon_d, &
490 clm_lake_initialized=clm_lake_initialized, input_lakedepth=input_lakedepth, &
491 tg3=tg3, clm_lakedepth=clm_lakedepth, km=km, me=me, master=master, &
492 errmsg=errmsg, errflg=errflg)
493 if(errflg/=0) then
494 return
495 endif
496
497 lake_points=0
498 snow_points=0
499 ice_points=0
500
501 to_radians = pi/180
502
503 month = idate(2)
504 day_of_month = idate(3)
505
506 num1 = month*2-1
507 if(day_of_month>15) then
508 num1 = num1 + 1
509 endif
510 num2 = num1+1
511
512 wght2 = day_of_month/month_length(month)
513 if(wght2<0 .or. wght2>1) then
514 if(debug_print) then
515 write(0,*) 'Warning: wght2 is not 0..1: ',wght2
516 endif
517 wght2 = max(0.0_kind_lake,min(1.0_kind_lake,wght2))
518 endif
519 wght1 = 1.0_kind_lake - wght2
520
521 if(debug_print ) then
522 print *,'month,num1,num2,wght1,wght2',month,num1,num2,wght1,wght2
523 endif
524
525 lake_top_loop: DO i = 1,im
526
527 if_lake_is_here: if (flag_iter(i) .and. use_lake_model(i)/=0) THEN
528
529 call is_salty(xlat_d(i),xlon_d(i),salty_flag,cannot_freeze_flag)
530
531 if(salty_flag) then
532 salty(i) = 1 ! The Great Salt Lake and Mono Lake
533 else
534 salty(i) = 0
535 endif
536
537 if(cannot_freeze_flag) then
538 cannot_freeze(i) = 1 ! only the Great Salt Lake
539 else
540 cannot_freeze(i) = 0
541 endif
542
543 sfctmp = gt0(i,1)
544 pbot = prsi(i,1)
545 psfc = pgr(i)
546 q2k = qvcurr(i)
547! EMISS = 0.99 * lake_icefrac3d(i,1) + emg * (1.0-lake_icefrac3d(i,1)) ! emg=0.97, parameter, needs to be moved to the top
548 lwdn = dlwsfci(i) ! LWDN is downward LW flux, do not use EMISS here.
549! LWDN = DLWSFCI(I)*EMISS(I)
550 ! FIXME: Should multiply PRCP by 1000
551 prcp = (raincprv(i)+rainncprv(i))/dtime ! [mm/s] use physics timestep since PRCP comes from non-surface schemes
552 soldn = dswsfci(i) ! SOLDN is total incoming solar
553 albedo(i) = ( 0.6 * lake_icefrac3d(i,1) ) + &
554 ( (1.0-lake_icefrac3d(i,1)) * 0.08)
555 solnet = soldn*(1.-albedo(i)) ! use mid-day albedo to determine net downward solar
556 ! (no solar zenith angle correction)
557
558 lake_points = lake_points+1
559
560 call calculate_z_dz_lake(i,input_lakedepth,clm_lakedepth,z_lake(1,:),dz_lake(1,:))
561
562 do c = 2,column
563 z_lake(c,:) = z_lake(1,:)
564 dz_lake(c,:) = dz_lake(1,:)
565 enddo
566
567 ! Soil hydraulic and thermal properties
568 isl = isltyp(i)
569 if (isl == 0 ) isl = 14
570 if (isl == 14 ) isl = isl + 1
571
572 watsat = 0.489_kind_lake - 0.00126_kind_lake*sand(isl)
573 csol = (2.128_kind_lake*sand(isl)+2.385_kind_lake*clay(isl)) / (sand(isl)+clay(isl))*1.e6_kind_lake ! J/(m3 K)
574 tkm = (8.80_kind_lake*sand(isl)+2.92_kind_lake*clay(isl))/(sand(isl)+clay(isl)) ! W/(m K)
575 bd = (1._kind_lake-watsat(1,1))*2.7e3_kind_lake
576 tkmg = tkm ** (1._kind_lake- watsat(1,1))
577 tkdry = (0.135_kind_lake*bd + 64.7_kind_lake) / (2.7e3_kind_lake - 0.947_kind_lake*bd)
578 tksatu = tkmg(1,1)*0.57_kind_lake**watsat(1,1)
579
580 do c = 1,column
581
582 forc_t(c) = sfctmp ! [K]
583 forc_pbot(c) = pbot ! [Pa]
584 forc_psrf(c) = psfc ! [Pa]
585 forc_hgt(c) = zlvl(i) ! [m]
586 forc_hgt_q(c) = zlvl(i) ! [m]
587 forc_hgt_t(c) = zlvl(i) ! [m]
588 forc_hgt_u(c) = zlvl(i) ! [m]
589 forc_q(c) = q2k ! [kg/kg]
590 forc_u(c) = gu0(i,1) ! [m/s]
591 forc_v(c) = gv0(i,1) ! [m/s]
592 forc_lwrad(c) = lwdn ! [W/m/m]
593 prec(c) = prcp ! [mm/s]
594 sabg(c) = solnet
595 lat(c) = xlat_d(i)*to_radians ! [radian]
596 do_capsnow(c) = .false.
597
598 lakedepth(c) = clm_lakedepth(i)
599 savedtke1(c) = savedtke12d(i)
600 snowdp(c) = snowdp2d(i)
601 h2osno(c) = h2osno2d(i)
602 snl(c) = snl2d(i)
603 t_grnd(c) = t_grnd2d(i)
604 do k = 1,nlevlake
605 t_lake(c,k) = t_lake3d(i,k)
606 lake_icefrac(c,k) = lake_icefrac3d(i,k)
607 enddo
608 do k = -nlevsnow+1,nlevsoil
609 t_soisno(c,k) = t_soisno3d(i,k)
610 h2osoi_ice(c,k) = h2osoi_ice3d(i,k)
611 h2osoi_liq(c,k) = h2osoi_liq3d(i,k)
612 h2osoi_vol(c,k) = h2osoi_vol3d(i,k)
613 z(c,k) = z3d(i,k)
614 dz(c,k) = dz3d(i,k)
615 enddo
616 do k = -nlevsnow+0,nlevsoil
617 zi(c,k) = zi3d(i,k)
618 enddo
619 enddo
620
621 eflx_lwrad_net = -9999
622 eflx_gnet = -9999
623 eflx_sh_tot = -9999
624 eflx_lh_tot = -9999
625 t_ref2m = -9999
626 q_ref2m = -9999
627 taux = -9999
628 tauy = -9999
629 ram1 = -9999
630 z0mg = -9999
631 ustar_out = -9999
632 lat_d = xlat_d(i)
633 lon_d = xlon_d(i)
634
635 CALL lakemain(forc_t,forc_pbot,forc_psrf,forc_hgt,forc_hgt_q, & !I
636 forc_hgt_t,forc_hgt_u,forc_q, forc_u, &
637 forc_v,forc_lwrad,prec, sabg,lat, &
638 z_lake,dz_lake,lakedepth,do_capsnow, &
639 h2osno,snowdp,snl,z,dz,zi, & !H
640 h2osoi_vol,h2osoi_liq,h2osoi_ice, &
641 t_grnd,t_soisno,t_lake, &
642 savedtke1,lake_icefrac, &
643 eflx_lwrad_net,eflx_gnet, & !O
644 eflx_sh_tot,eflx_lh_tot, &
645 t_ref2m,q_ref2m, dtime, &
646 watsat, tksatu, tkmg, tkdry, csol, &
647 taux,tauy,ram1,z0mg,ustar_out,errmsg,errflg, &
648 lat_d,lon_d)
649 ! Renew Lake State Variables:(14)
650 do c = 1,column
651
652 if(cannot_freeze(i) == 1) then
653 ! The Great Salt Lake
654 do k = 1,nlevlake
655 lake_icefrac(c,k) = 0._kind_lake
656 enddo
657 ! bound lake temperture with the climatology
658 tclim = tfrz + wght1*saltlk_t(num1) &
659 + wght2*saltlk_t(num2)
660 if(debug_print) print *,'GSL - Tclim,tsfc,t_lake3d',i,tclim,t_grnd(c),t_lake(c,:),t_soisno(c,:)
661 t_grnd(c) = min(tclim+3.0_kind_lake,(max(t_grnd(c),tclim-3.0_kind_lake)))
662 do k = 1,nlevlake
663 t_lake(c,k) = min(tclim+3.0_kind_lake,(max(t_lake(c,k),tclim-3.0_kind_lake)))
664 enddo
665 t_soisno(c,1) = min(tclim+3.0_kind_lake,(max(t_soisno(c,1),tclim-3.0_kind_lake)))
666 if(debug_print) print *,'GSL - after Tclim,tsfc,t_lake3d',i,tclim,t_grnd(c),t_lake(c,:),t_soisno(c,:)
667 elseif(salty(i) == 1) then
668 ! Mono Lake never freezes, its temperature is above freezing point = -4.2 C
669 t_grnd(c) = max(tfrz-4.2_kind_lake,t_grnd(c))
670 do k = 1,nlevlake
671 lake_icefrac(c,k) = 0._kind_lake
672 t_lake(c,k) = max(tfrz-4.2_kind_lake,t_lake(c,k))
673 enddo
674 t_soisno(c,1) = max(tfrz-4.2_kind_lake,t_soisno(c,1))
675 if(debug_print) print *,'Mono - tsfc,t_lake3d',i,t_grnd(c),t_lake(c,:),t_soisno(c,:)
676 endif
677
678 savedtke12d(i) = savedtke1(c)
679 snowdp2d(i) = snowdp(c)
680 h2osno2d(i) = h2osno(c)
681 snl2d(i) = snl(c)
682 t_grnd2d(i) = t_grnd(c)
683 do k = 1,nlevlake
684 t_lake3d(i,k) = t_lake(c,k)
685 lake_icefrac3d(i,k) = lake_icefrac(c,k)
686 enddo
687 do k = -nlevsnow+1,nlevsoil
688 z3d(i,k) = z(c,k)
689 dz3d(i,k) = dz(c,k)
690 t_soisno3d(i,k) = t_soisno(c,k)
691 h2osoi_liq3d(i,k) = h2osoi_liq(c,k)
692 h2osoi_ice3d(i,k) = h2osoi_ice(c,k)
693 h2osoi_vol3d(i,k) = h2osoi_vol(c,k)
694 enddo
695 do k = -nlevsnow+0,nlevsoil
696 zi3d(i,k) = zi(c,k)
697 enddo
698
699 enddo
700
701 feedback: if(feedback_to_atmosphere) then
702 c = 1
703
704 !-- The CLM output is combined for fractional ice and water
705 if( t_grnd(c) >= tfrz ) then
706 qfx = eflx_lh_tot(c)*invhvap
707 else
708 qfx = eflx_lh_tot(c)*invhsub ! heat flux (W/m^2)=>mass flux(kg/(sm^2))
709 endif
710 rho0 = prsl1(i) / (rair*t1(i)*(1.0 + con_fvirt*qv1(i)))
711 evap_wat(i) = qfx/rho0 ! kinematic_surface_upward_latent_heat_flux_over_water
712 hflx_wat(i) = eflx_sh_tot(c)/(rho0*cpair) ! kinematic_surface_upward_sensible_heat_flux_over_water
713 gflx_wat(i) = eflx_gnet(c) ![W/m/m] upward_heat_flux_in_soil_over_water
714 ep1d_water(i) = eflx_lh_tot(c) ![W/m/m] surface_upward_potential_latent_heat_flux_over_water
715 tsurf_water(i) = t_grnd(c) ![K] surface skin temperature after iteration over water
716 tsurf_ice(i) = t_grnd(c) ! surface_skin_temperature_after_iteration_over_ice
717 tsfc_wat(i) = t_grnd(c) ![K] surface skin temperature over water
718 tisfc(i) = t_grnd(c)
719 tsfc(i) = t_grnd(c)
720 lake_t2m(i) = t_ref2m(c) ![K] temperature_at_2m_from_clm_lake
721 lake_q2m(i) = q_ref2m(c) ! [frac] specific_humidity_at_2m_from_clm_lake
722 albedo(i) = ( 0.6 * lake_icefrac3d(i,1) ) + & ! mid_day_surface_albedo_over_lake
723 ( (1.0-lake_icefrac3d(i,1)) * 0.08)
724 fice(i) = lake_icefrac3d(i,1) ! sea_ice_area_fraction_of_sea_area_fraction
725 !uustar_water(i) = ustar_out(c) ! surface_friction_velocity_over_water
726 zorlw(i) = z0mg(c) ! surface_roughness_length_over_water
727
728 ! WRF variables with no equivalent in CCPP:
729 ! LH(I) = eflx_lh_tot(c)/rho1(i) ![kg*m/(kg*s)]
730 !TH2(I) = T2(I)*(1.E5/PSFC)**RCP ! potential temperature
731
732 ! Calculate qsfc from t_grnd: ! surface_specific_humidity_over_water
733 psfc = prsi(i,1)
734 discard1 = -9999
735 discard2 = -9999
736 discard3 = -9999
737 qss = qss_water(i)
738 call qsat(t_grnd(c),psfc,discard1,discard2,qss,discard3)
739 qss_water(i) = qss
740
741 ! Combined water-ice chh and cmm calculations come from Flake model:
742 chh_water(i) = ch(i)*wind(i)*1.225 ! surface_drag_mass_flux_for_heat_and_moisture_in_air_over_water
743 cmm_water(i) = cm(i)*wind(i) ! surface_drag_wind_speed_for_momentum_in_air_over_water
744
745 ice_point: if(fice(i)>=min_lakeice) then
746 ! Icy lake
747 ! Most ice variables are identical to water variables.
748 if(debug_print) print *,'Icy xlat_d(i),xlon_d(i),frac_ice,frac_grid ',xlat_d(i),xlon_d(i),frac_ice,frac_grid
749 if(frac_ice .or. frac_grid) then
750 evap_ice(i) = evap_wat(i) ! kinematic_surface_upward_latent_heat_flux_over_ice
751 hflx_ice(i) = hflx_wat(i) ! kinematic_surface_upward_sensible_heat_flux_over_ice
752 gflx_ice(i) = gflx_wat(i) ! upward_heat_flux_in_soil_over_ice
753 ep1d_ice(i) = ep1d_water(i) ! surface_upward_potential_latent_heat_flux_over_ice
754 chh_ice(i) = chh_water(i) ! surface_drag_mass_flux_for_heat_and_moisture_in_air_over_ice
755 cmm_ice(i) = cmm_water(i) ! surface_drag_wind_speed_for_momentum_in_air_over_ice
756 qss_ice(i) = qss_water(i) ! surface_specific_humidity_over_ice
757! uustar_ice(i) = uustar_water(i) ! surface_friction_velocity_over_ice
758 endif
759
760 tsurf_ice(i) = t_grnd(c) ! surface_skin_temperature_after_iteration_over_ice
761 tisfc(i) = t_grnd(c) ! surface_skin_temperature_over_ice
762 tsfc(i) = t_grnd(c) ! surface_skin_temperature_over_ice
763 weasdi(i) = h2osno(c) ! water_equivalent_accumulated_snow_depth_over_ice
764 snodi(i) = snowdp(c)*1.e3 ! surface_snow_thickness_water_equivalent_over_ice
765 weasd(i) = weasdi(i)
766 snowd(i) = snodi(c) ! surface_snow_thickness_water_equivalent_over_ice
767
768
769 if (.not. icy(i)) then
770 flag_lakefreeze(i)=.true.
771 end if
772
773 ! Ice points are icy:
774 icy(i)=.true. ! flag_nonzero_sea_ice_surface_fraction
775 ice_points = ice_points+1
776
777 zorli(i) = z0mg(c) ! surface_roughness_length_over_ice
778
779 ! Assume that, if a layer has ice, the entire layer thickness is ice.
780 hice(i) = 0 ! sea_ice_thickness
781 do k=1,nlevlake
782 if(lake_icefrac3d(i,k)>0) then
783 hice(i) = hice(i) + dz_lake(c,k)
784 endif
785 end do
786 else ! Not an ice point
787 ! On non-icy lake points, set variables relevant to
788 ! lake ice to reasonable defaults. Let LSM fill in
789 ! other variables.
790 icy(i)=.false.
791 weasdi(i) = 0
792 snodi(i) = 0
793 weasd(i) = 0
794 snowd(i) = 0
795 tisfc(i) = t_grnd(c)
796 tsurf_ice(i) = tisfc(i)
797 tsfc(i) = t_grnd(c)
798 hice(i) = 0
799 fice(i) = 0
800 endif ice_point
801
802 if(snl2d(i)<0) then
803 ! If there is snow, ice surface temperature should be snow temperature.
804 lake_t_snow(i) = t_grnd(c) ! surface_skin_temperature_over_ice
805 tisfc(i) = lake_t_snow(i) ! temperature_of_snow_on_lake
806 snow_points = snow_points+1
807 else
808 lake_t_snow(i) = -9999
809 endif
810
811 endif feedback
812
813 endif if_lake_is_here
814 ENDDO lake_top_loop
815
816 if(debug_print .and. lake_points>0 .and. (kdt<3 .or. mod(kdt,30)==3)) then
8173082 format('lake points processed in timestep ',i0,' by rank ',i0,' = ',i0,' snow=',i0,' ice=',i0)
818 print 3082,kdt,me,lake_points,snow_points,ice_points
819 endif
820
821 END SUBROUTINE clm_lake_run
822
823
824 SUBROUTINE lakemain(forc_t,forc_pbot,forc_psrf,forc_hgt,forc_hgt_q, & !I
825 forc_hgt_t,forc_hgt_u,forc_q, forc_u, &
826 forc_v,forc_lwrad,prec, sabg,lat, &
827 z_lake,dz_lake,lakedepth,do_capsnow, &
828 h2osno,snowdp,snl,z,dz,zi, & !H
829 h2osoi_vol,h2osoi_liq,h2osoi_ice, &
830 t_grnd,t_soisno,t_lake, &
831 savedtke1,lake_icefrac, &
832 eflx_lwrad_net,eflx_gnet, & !O
833 eflx_sh_tot,eflx_lh_tot, &
834 t_ref2m,q_ref2m, dtime, &
835 watsat, tksatu, tkmg, tkdry, csol, &
836 taux,tauy,ram1,z0mg,ustar_out,errmsg,errflg, xlat_d,xlon_d)
837 implicit none
838 !in:
839
840 integer, intent(inout) :: errflg
841 character(*), intent(inout) :: errmsg
842 real(kind_lake),intent(in) :: dtime
843 real(kind_lake),intent(in) :: xlat_d, xlon_d
844 real(kind_lake),intent(in) :: forc_t(1)
845 real(kind_lake),intent(in) :: forc_pbot(1)
846 real(kind_lake),intent(in) :: forc_psrf(1)
847 real(kind_lake),intent(in) :: forc_hgt(1)
848 real(kind_lake),intent(in) :: forc_hgt_q(1)
849 real(kind_lake),intent(in) :: forc_hgt_t(1)
850 real(kind_lake),intent(in) :: forc_hgt_u(1)
851 real(kind_lake),intent(in) :: forc_q(1)
852 real(kind_lake),intent(in) :: forc_u(1)
853 real(kind_lake),intent(in) :: forc_v(1)
854 ! real(kind_lake),intent(in) :: forc_rho(1) ! density (kg/m**3)
855 real(kind_lake),intent(in) :: forc_lwrad(1)
856 real(kind_lake),intent(in) :: prec(1)
857 real(kind_lake),intent(in) :: sabg(1)
858 real(kind_lake),intent(in) :: lat(1)
859 real(kind_lake),intent(in) :: z_lake(1,nlevlake)
860 real(kind_lake),intent(in) :: dz_lake(1,nlevlake)
861 real(kind_lake),intent(out) :: ustar_out(1)
862 real(kind_lake), intent(in) :: lakedepth(1)!!!!!!!!!!!!!!tep(in),hydro(in)
864 ! real(kind_lake), intent(in) :: watsat(1,1:nlevsoil) ! volumetric soil water at saturation (porosity)
865 !!!!!!!!!!!!!!!!hydro
866 logical , intent(in) :: do_capsnow(1)
867 real(kind_lake), intent(in) :: watsat(1,nlevsoil)
868 real(kind_lake), intent(in) :: tksatu(1,nlevsoil)
869 real(kind_lake), intent(in) :: tkmg(1,nlevsoil)
870 real(kind_lake), intent(in) :: tkdry(1,nlevsoil)
871 real(kind_lake), intent(in) :: csol(1,nlevsoil)
872
873
874
875 !in&out
876 real(kind_lake),intent(inout) :: h2osoi_vol(1,-nlevsnow+1:nlevsoil)
877 real(kind_lake),intent(inout) :: t_grnd(1)
878 real(kind_lake),intent(inout) :: h2osno(1)
879 real(kind_lake),intent(inout) :: snowdp(1)
880 real(kind_lake),intent(inout) :: z(1,-nlevsnow+1:nlevsoil)
881 real(kind_lake),intent(inout) :: dz(1,-nlevsnow+1:nlevsoil)
882 real(kind_lake),intent(inout) :: t_soisno(1,-nlevsnow+1:nlevsoil)
883 real(kind_lake),intent(inout) :: t_lake(1,nlevlake)
884 integer ,intent(inout) :: snl(1)
885 real(kind_lake),intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil)
886 real(kind_lake),intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil)
887 real(kind_lake),intent(inout) :: savedtke1(1)
888 real(kind_lake),intent(inout) :: zi(1,-nlevsnow+0:nlevsoil)
889 real(kind_lake),intent(inout) :: lake_icefrac(1,nlevlake)
890
891
892 !out:
893 real(kind_lake),intent(out) :: eflx_gnet(1)
894 real(kind_lake),intent(out) :: eflx_lwrad_net(1)
895 real(kind_lake),intent(out) :: eflx_sh_tot(1)
896 real(kind_lake),intent(out) :: eflx_lh_tot(1)
897 real(kind_lake),intent(out) :: t_ref2m(1)
898 real(kind_lake),intent(out) :: q_ref2m(1)
899 real(kind_lake),intent(out) :: taux(1)
900 real(kind_lake),intent(out) :: tauy(1)
901 real(kind_lake),intent(out) :: ram1(1)
904 real(kind_lake),intent(out) :: z0mg(1)
905
906
907 !local output
908
909 real(kind_lake) :: begwb(1) ! water mass begining of the time step
910 real(kind_lake) :: t_veg(1) ! vegetation temperature (Kelvin)
911 real(kind_lake) :: eflx_soil_grnd(1) ! soil heat flux (W/m**2) [+ = into soil]
912 real(kind_lake) :: eflx_lh_grnd(1) ! ground evaporation heat flux (W/m**2) [+ to atm]
913 real(kind_lake) :: eflx_sh_grnd(1) ! sensible heat flux from ground (W/m**2) [+ to atm]
914 real(kind_lake) :: eflx_lwrad_out(1) ! emitted infrared (longwave) radiation (W/m**2)
915 real(kind_lake) :: qflx_evap_tot(1) ! qflx_evap_soi + qflx_evap_veg + qflx_tran_veg
916 real(kind_lake) :: qflx_evap_soi(1) ! soil evaporation (mm H2O/s) (+ = to atm)
917 real(kind_lake) :: qflx_prec_grnd(1) ! water onto ground including canopy runoff [kg/(m2 s)]
918 real(kind_lake) :: forc_snow(1) ! snow rate [mm/s]
919 real(kind_lake) :: forc_rain(1) ! rain rate [mm/s]
920 real(kind_lake) :: ws(1) ! surface friction velocity (m/s)
921 real(kind_lake) :: ks(1) ! coefficient passed to ShalLakeTemperature
922 real(kind_lake) :: qflx_snomelt(1) !snow melt (mm H2O /s) tem(out),snowwater(in)
923 integer :: imelt(1,-nlevsnow+1:nlevsoil) !flag for melting (=1), freezing (=2), Not=0 (new)
924 real(kind_lake) :: endwb(1) ! water mass end of the time step
925 real(kind_lake) :: snowage(1) ! non dimensional snow age [-]
926 real(kind_lake) :: snowice(1) ! average snow ice lens
927 real(kind_lake) :: snowliq(1) ! average snow liquid water
928 real(kind_lake) :: t_snow(1) ! vertically averaged snow temperature
929 real(kind_lake) :: qflx_drain(1) ! sub-surface runoff (mm H2O /s)
930 real(kind_lake) :: qflx_surf(1) ! surface runoff (mm H2O /s)
931 real(kind_lake) :: qflx_infl(1) ! infiltration (mm H2O /s)
932 real(kind_lake) :: qflx_qrgwl(1) ! qflx_surf at glaciers, wetlands, lakes
933 real(kind_lake) :: qcharge(1) ! aquifer recharge rate (mm/s)
934 real(kind_lake) :: qflx_snowcap(1) ! excess precipitation due to snow capping (mm H2O /s) [+]
935 real(kind_lake) :: qflx_snowcap_col(1) ! excess precipitation due to snow capping (mm H2O /s) [+]
936 real(kind_lake) :: qflx_snow_grnd_pft(1) ! snow on ground after interception (mm H2O/s) [+]
937 real(kind_lake) :: qflx_snow_grnd_col(1) ! snow on ground after interception (mm H2O/s) [+]
938 real(kind_lake) :: qflx_rain_grnd(1) ! rain on ground after interception (mm H2O/s) [+]
939 real(kind_lake) :: frac_iceold(1,-nlevsnow+1:nlevsoil) ! fraction of ice relative to the tot water
940 real(kind_lake) :: qflx_evap_tot_col(1) !pft quantity averaged to the column (assuming one pft)
941 real(kind_lake) :: soilalpha(1) !factor that reduces ground saturated specific humidity (-)
942 real(kind_lake) :: zwt(1) !water table depth
943 real(kind_lake) :: fcov(1) !fractional area with water table at surface
944 real(kind_lake) :: rootr_column(1,1:nlevsoil) !effective fraction of roots in each soil layer
945 real(kind_lake) :: qflx_evap_grnd(1) ! ground surface evaporation rate (mm H2O/s) [+]
946 real(kind_lake) :: qflx_sub_snow(1) ! sublimation rate from snow pack (mm H2O /s) [+]
947 real(kind_lake) :: qflx_dew_snow(1) ! surface dew added to snow pack (mm H2O /s) [+]
948 real(kind_lake) :: qflx_dew_grnd(1) ! ground surface dew formation (mm H2O /s) [+]
949 real(kind_lake) :: qflx_rain_grnd_col(1) !rain on ground after interception (mm H2O/s) [+]
950 begwb = 0
951
952 ! lat = lat*pi/180 ! [radian]
953
954 if (prec(1)> 0.) then
955 if ( forc_t(1) > (tfrz + tcrit)) then
956 forc_rain(1) = prec(1)
957 forc_snow(1) = 0.
958 ! flfall(1) = 1.
959 else
960 forc_rain(1) = 0.
961 forc_snow(1) = prec(1)
962
963 ! if ( forc_t(1) <= tfrz) then
964 ! flfall(1) = 0.
965 ! else if ( forc_t(1) <= tfrz+2.) then
966 ! flfall(1) = -54.632 + 0.2 * forc_t(1)
967 ! else
968 ! flfall(1) = 0.4
969 endif
970 else
971 forc_rain(1) = 0.
972 forc_snow(1) = 0.
973 ! flfall(1) = 1.
974 endif
975
976 CALL shallakefluxes(forc_t,forc_pbot,forc_psrf,forc_hgt,forc_hgt_q, & !i
977 forc_hgt_t,forc_hgt_u,forc_q, &
978 forc_u,forc_v,forc_lwrad,forc_snow, &
979 forc_rain,t_grnd,h2osno,snowdp,sabg,lat, &
980 dz,dz_lake,t_soisno,t_lake,snl,h2osoi_liq, &
981 h2osoi_ice,savedtke1, &
982 qflx_prec_grnd,qflx_evap_soi,qflx_evap_tot, & !o
983 eflx_sh_grnd,eflx_lwrad_out,eflx_lwrad_net, &
984 eflx_soil_grnd,eflx_sh_tot,eflx_lh_tot, &
985 eflx_lh_grnd,t_veg,t_ref2m,q_ref2m,taux,tauy, &
986 ram1,ws,ks,eflx_gnet,z0mg,ustar_out,errmsg,errflg,xlat_d,xlon_d)
987 if(errflg/=0) then
988 return ! State is invalid now, so pass error to caller.
989 endif
990
991 CALL shallaketemperature(t_grnd,h2osno,sabg,dz,dz_lake,z,zi, & !i
992 z_lake,ws,ks,snl,eflx_gnet,lakedepth, &
993 lake_icefrac,snowdp, & !i&o
994 eflx_sh_grnd,eflx_sh_tot,eflx_soil_grnd, & !o
995 t_lake,t_soisno,h2osoi_liq, &
996 h2osoi_ice,savedtke1, &
997 watsat, tksatu, tkmg, tkdry, csol, dtime, &
998 frac_iceold,qflx_snomelt,imelt,errmsg,errflg,xlat_d,xlon_d)
999 if(errflg/=0) then
1000 return ! State is invalid now, so pass error to caller.
1001 endif
1002
1003 CALL shallakehydrology(dz_lake,forc_rain,forc_snow, & !i
1004 begwb,qflx_evap_tot,forc_t,do_capsnow, &
1005 t_grnd,qflx_evap_soi, &
1006 qflx_snomelt,imelt,frac_iceold, & !i add by guhp
1007 z,dz,zi,snl,h2osno,snowdp,lake_icefrac,t_lake, & !i&o
1008 endwb,snowage,snowice,snowliq,t_snow, & !o
1009 t_soisno,h2osoi_ice,h2osoi_liq,h2osoi_vol, &
1010 qflx_drain,qflx_surf,qflx_infl,qflx_qrgwl, &
1011 qcharge,qflx_prec_grnd,qflx_snowcap, &
1012 qflx_snowcap_col,qflx_snow_grnd_pft, &
1013 qflx_snow_grnd_col,qflx_rain_grnd, &
1014 qflx_evap_tot_col,soilalpha,zwt,fcov, &
1015 rootr_column,qflx_evap_grnd,qflx_sub_snow, &
1016 qflx_dew_snow,qflx_dew_grnd,qflx_rain_grnd_col, &
1017 watsat, tksatu, tkmg, tkdry, csol, &
1018 dtime,errmsg,errflg)
1019 if(errflg/=0) then
1020 return ! State is invalid now, so pass error to caller.
1021 endif
1022
1023 !==================================================================================
1024 ! !DESCRIPTION:
1025 ! Calculation of Shallow Lake Hydrology. Full hydrology of snow layers is
1026 ! done. However, there is no infiltration, and the water budget is balanced with
1027
1028 END SUBROUTINE lakemain
1029
1030
1031 ! DESCRIPTION:
1038SUBROUTINE shallakefluxes(forc_t,forc_pbot,forc_psrf,forc_hgt,forc_hgt_q, & !i
1039 forc_hgt_t,forc_hgt_u,forc_q, &
1040 forc_u,forc_v,forc_lwrad,forc_snow, &
1041 forc_rain,t_grnd,h2osno,snowdp,sabg,lat, &
1042 dz,dz_lake,t_soisno,t_lake,snl,h2osoi_liq, &
1043 h2osoi_ice,savedtke1, &
1044 qflx_prec_grnd,qflx_evap_soi,qflx_evap_tot, & !o
1045 eflx_sh_grnd,eflx_lwrad_out,eflx_lwrad_net, &
1046 eflx_soil_grnd,eflx_sh_tot,eflx_lh_tot, &
1047 eflx_lh_grnd,t_veg,t_ref2m,q_ref2m,taux,tauy, &
1048 ram1,ws,ks,eflx_gnet,z0mg,ustar_out,errmsg,errflg,xlat_d,xlon_d)
1049 !==============================================================================
1050 ! REVISION HISTORY:
1051 ! - Created by Zack Subin, 2009
1052 ! - Reedited by Hongping Gu, 2010
1053 ! - Updated for CCPP by Sam Trahan, 2022
1054 !==============================================================================
1055
1056 ! implicit none
1057
1058 implicit none
1059
1060 !in:
1061
1062 integer, intent(inout) :: errflg
1063 character(len=*), intent(inout) :: errmsg
1064 real(kind_lake),intent(in) :: xlat_d,xlon_d
1065 real(kind_lake),intent(in) :: forc_t(1)
1066 real(kind_lake),intent(in) :: forc_pbot(1)
1067 real(kind_lake),intent(in) :: forc_psrf(1)
1068 real(kind_lake),intent(in) :: forc_hgt(1)
1069 real(kind_lake),intent(in) :: forc_hgt_q(1)
1070 real(kind_lake),intent(in) :: forc_hgt_t(1)
1071 real(kind_lake),intent(in) :: forc_hgt_u(1)
1072 real(kind_lake),intent(in) :: forc_q(1)
1073 real(kind_lake),intent(in) :: forc_u(1)
1074 real(kind_lake),intent(in) :: forc_v(1)
1075 real(kind_lake),intent(in) :: forc_lwrad(1)
1076 ! real(kind_lake),intent(in) :: forc_rho(1) ! density (kg/m**3)
1077 real(kind_lake),intent(in) :: forc_snow(1)
1078 real(kind_lake),intent(in) :: forc_rain(1)
1079 real(kind_lake),intent(in) :: h2osno(1)
1080 real(kind_lake),intent(in) :: snowdp(1)
1081 real(kind_lake),intent(in) :: sabg(1)
1082 real(kind_lake),intent(in) :: lat(1)
1083 real(kind_lake),intent(in) :: dz(1,-nlevsnow+1:nlevsoil)
1084 real(kind_lake),intent(in) :: dz_lake(1,nlevlake)
1085 real(kind_lake),intent(in) :: t_soisno(1,-nlevsnow+1:nlevsoil)
1086 real(kind_lake),intent(in) :: t_lake(1,nlevlake)
1087 integer ,intent(in) :: snl(1)
1088 real(kind_lake),intent(in) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil)
1089 real(kind_lake),intent(in) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil)
1090 real(kind_lake),intent(in) :: savedtke1(1)
1091
1092 !inout:
1093 real(kind_lake),intent(inout) :: t_grnd(1)
1094 !out:
1095 real(kind_lake),intent(out):: ustar_out(1)
1096 real(kind_lake),intent(out):: qflx_prec_grnd(1)
1097 real(kind_lake),intent(out):: qflx_evap_soi(1)
1098 real(kind_lake),intent(out):: qflx_evap_tot(1)
1099 real(kind_lake),intent(out):: eflx_sh_grnd(1)
1100 real(kind_lake),intent(out):: eflx_lwrad_out(1)
1101 real(kind_lake),intent(out):: eflx_lwrad_net(1)
1102 real(kind_lake),intent(out):: eflx_soil_grnd(1)
1103 real(kind_lake),intent(out):: eflx_sh_tot(1)
1104 real(kind_lake),intent(out):: eflx_lh_tot(1)
1105 real(kind_lake),intent(out):: eflx_lh_grnd(1)
1106 real(kind_lake),intent(out):: t_veg(1)
1107 real(kind_lake),intent(out):: t_ref2m(1)
1108 real(kind_lake),intent(out):: q_ref2m(1)
1109 real(kind_lake),intent(out):: taux(1)
1110 real(kind_lake),intent(out):: tauy(1)
1111 real(kind_lake),intent(out):: ram1(1)
1112 real(kind_lake),intent(out):: ws(1)
1113 real(kind_lake),intent(out):: ks(1)
1115 real(kind_lake),intent(out):: eflx_gnet(1)
1117 real(kind_lake),intent(out):: z0mg(1)
1118
1119
1120
1121 !OTHER LOCAL VARIABLES:
1122
1123 integer , parameter :: islak = 2 ! index of lake, 1 = deep lake, 2 = shallow lake
1124 integer , parameter :: niters = 3 ! maximum number of iterations for surface temperature
1125 real(kind_lake), parameter :: beta1 = 1._kind_lake ! coefficient of convective velocity (in computing W_*) [-]
1126 real(kind_lake), parameter :: emg = 0.97_kind_lake ! ground emissivity (0.97 for water)
1127 real(kind_lake), parameter :: zii = 1000._kind_lake! convective boundary height [m]
1128 real(kind_lake), parameter :: tdmax = 277._kind_lake ! temperature of maximum water density
1129 real(kind_lake) :: forc_th(1) ! atmospheric potential temperature (Kelvin)
1130 real(kind_lake) :: forc_vp(1) !atmospheric vapor pressure (Pa)
1131 real(kind_lake) :: forc_rho(1) ! density (kg/m**3)
1132 integer :: i,fc,fp,g,c,p ! do loop or array index
1133 integer :: fncopy ! number of values in pft filter copy
1134 integer :: fnold ! previous number of pft filter values
1135 integer :: fpcopy(num_shlakep) ! pft filter copy for iteration loop
1136 integer :: iter ! iteration index
1137 integer :: nmozsgn(lbp:ubp) ! number of times moz changes sign
1138 integer :: jtop(lbc:ubc) ! top level for each column (no longer all 1)
1139 real(kind_lake) :: ax ! used in iteration loop for calculating t_grnd (numerator of NR solution)
1140 real(kind_lake) :: bx ! used in iteration loop for calculating t_grnd (denomin. of NR solution)
1141 real(kind_lake) :: degdT ! d(eg)/dT
1142 real(kind_lake) :: dqh(lbp:ubp) ! diff of humidity between ref. height and surface
1143 real(kind_lake) :: dth(lbp:ubp) ! diff of virtual temp. between ref. height and surface
1144 real(kind_lake) :: dthv ! diff of vir. poten. temp. between ref. height and surface
1145 real(kind_lake) :: dzsur(lbc:ubc) ! 1/2 the top layer thickness (m)
1146 real(kind_lake) :: eg ! water vapor pressure at temperature T [pa]
1147 real(kind_lake) :: htvp(lbc:ubc) ! latent heat of vapor of water (or sublimation) [j/kg]
1148 real(kind_lake) :: obu(lbp:ubp) ! monin-obukhov length (m)
1149 real(kind_lake) :: obuold(lbp:ubp) ! monin-obukhov length of previous iteration
1150 real(kind_lake) :: qsatg(lbc:ubc) ! saturated humidity [kg/kg]
1151 real(kind_lake) :: qsatgdT(lbc:ubc) ! d(qsatg)/dT
1152 real(kind_lake) :: qstar ! moisture scaling parameter
1153 real(kind_lake) :: ram(lbp:ubp) ! aerodynamical resistance [s/m]
1154 real(kind_lake) :: rah(lbp:ubp) ! thermal resistance [s/m]
1155 real(kind_lake) :: raw(lbp:ubp) ! moisture resistance [s/m]
1156 real(kind_lake) :: stftg3(lbp:ubp) ! derivative of fluxes w.r.t ground temperature
1157 real(kind_lake) :: temp1(lbp:ubp) ! relation for potential temperature profile
1158 real(kind_lake) :: temp12m(lbp:ubp) ! relation for potential temperature profile applied at 2-m
1159 real(kind_lake) :: temp2(lbp:ubp) ! relation for specific humidity profile
1160 real(kind_lake) :: temp22m(lbp:ubp) ! relation for specific humidity profile applied at 2-m
1161 real(kind_lake) :: tgbef(lbc:ubc) ! initial ground temperature
1162 real(kind_lake) :: thm(lbc:ubc) ! intermediate variable (forc_t+0.0098*forc_hgt_t)
1163 real(kind_lake) :: thv(lbc:ubc) ! virtual potential temperature (kelvin)
1164 real(kind_lake) :: thvstar ! virtual potential temperature scaling parameter
1165 real(kind_lake) :: tksur ! thermal conductivity of snow/soil (w/m/kelvin)
1166 real(kind_lake) :: tsur ! top layer temperature
1167 real(kind_lake) :: tstar ! temperature scaling parameter
1168 real(kind_lake) :: um(lbp:ubp) ! wind speed including the stablity effect [m/s]
1169 real(kind_lake) :: ur(lbp:ubp) ! wind speed at reference height [m/s]
1170 real(kind_lake) :: ustar(lbp:ubp) ! friction velocity [m/s]
1171 real(kind_lake) :: wc ! convective velocity [m/s]
1172 real(kind_lake) :: zeta ! dimensionless height used in Monin-Obukhov theory
1173 real(kind_lake) :: zldis(lbp:ubp) ! reference height "minus" zero displacement height [m]
1174 real(kind_lake) :: displa(lbp:ubp) ! displacement (always zero) [m]
1175 ! real(kind_lake) :: z0mg(lbp:ubp) ! roughness length over ground, momentum [m]
1176 real(kind_lake) :: z0hg(lbp:ubp) ! roughness length over ground, sensible heat [m]
1177 real(kind_lake) :: z0qg(lbp:ubp) ! roughness length over ground, latent heat [m]
1178 real(kind_lake) :: u2m ! 2 m wind speed (m/s)
1179 real(kind_lake) :: u10(1) ! 10-m wind (m/s) (for dust model)
1180 real(kind_lake) :: fv(1) ! friction velocity (m/s) (for dust model)
1181
1182 real(kind_lake) :: fm(lbp:ubp) ! needed for BGC only to diagnose 10m wind speed
1183 real(kind_lake) :: bw ! partial density of water (ice + liquid)
1184 real(kind_lake) :: t_grnd_temp ! Used in surface flux correction over frozen ground
1185 real(kind_lake) :: betaprime(lbc:ubc) ! Effective beta: 1 for snow layers, beta(islak) otherwise
1186 character*256 :: message
1187 ! tgs COARE
1188 real(kind_lake) :: tc, visc, ren
1189
1190 ! This assumes all radiation is absorbed in the top snow layer and will need
1191 ! to be changed for CLM 4.
1192 !
1193 ! Constants for lake temperature model
1194 !
1195 real(kind_lake), parameter :: beta(2) = & ! fraction solar rad absorbed at surface: depends on lake type
1196 (/0.4_kind_lake, 0.4_kind_lake/) ! (deep lake, shallow lake)
1197 ! This is the energy absorbed at the lake surface if no snow.
1198 ! data za /0.6_kind_lake, 0.5_kind_lake/
1199 ! data eta /0.1_kind_lake, 0.5_kind_lake/
1200 !-----------------------------------------------------------------------
1201
1202 ! Begin calculations
1203
1204 !dir$ concurrent
1205 !cdir nodep
1206 forc_th(1) = forc_t(1) * (forc_psrf(1)/ forc_pbot(1))**(rair/cpair)
1207 forc_vp(1) = forc_q(1) * forc_pbot(1)/ (con_eps + one_minus_con_eps * forc_q(1))
1208 forc_rho(1) = (forc_pbot(1) - one_minus_con_eps * forc_vp(1)) / (rair * forc_t(1))
1209
1210 do fc = 1, num_shlakec
1211 c = filter_shlakec(fc)
1212 g = cgridcell(c)
1213
1214 ! Surface temperature and fluxes
1215
1216 ! Find top layer
1217 if (snl(c) > 0 .or. snl(c) < -5) then
1218 errmsg='snl is not defined in ShalLakeFluxesMod; snl: out of range value'
1219 errflg=1
1220 return ! Cannot continue
1221 end if
1222 ! if (snl(c) /= 0) then
1223 ! write(6,*)'snl is not equal to zero in ShalLakeFluxesMod'
1224 ! call endrun()
1225 ! end if
1226 jtop(c) = snl(c) + 1
1227
1228
1229 if (snl(c) < 0) then
1230 betaprime(c) = 1._kind_lake !Assume all solar rad. absorbed at the surface of the top snow layer.
1231 dzsur(c) = dz(c,jtop(c))*0.5_kind_lake
1232 else
1233 betaprime(c) = beta(islak)
1234 dzsur(c) = dz_lake(c,1)*0.5_kind_lake
1235 end if
1236 ! Originally this was 1*dz, but shouldn't it be 1/2?
1237
1238 ! Saturated vapor pressure, specific humidity and their derivatives
1239 ! at lake surface
1240
1241 call qsat(t_grnd(c), forc_pbot(g), eg, degdt, qsatg(c), qsatgdt(c))
1242
1243 ! Potential, virtual potential temperature, and wind speed at the
1244 ! reference height
1245
1246 thm(c) = forc_t(g) + 0.0098_kind_lake*forc_hgt_t(g) ! intermediate variable
1247 thv(c) = forc_th(g)*(1._kind_lake+con_fvirt*forc_q(g)) ! virtual potential T
1248 end do
1249
1250 !dir$ concurrent
1251 !cdir nodep
1252 do fp = 1, num_shlakep
1253 p = filter_shlakep(fp)
1254 c = pcolumn(p)
1255 g = pgridcell(p)
1256
1257 nmozsgn(p) = 0
1258 obuold(p) = 0._kind_lake
1259 displa(p) = 0._kind_lake
1260
1261 ! Roughness lengths
1262
1263
1264 ! changed by Hongping Gu
1265 ! if (t_grnd(c) >= tfrz) then ! for unfrozen lake
1266 ! z0mg(p) = 0.01_kind_lake
1267 ! else ! for frozen lake
1268 ! ! Is this okay even if it is snow covered? What is the roughness over
1269 ! non-veg. snow?
1270 ! z0mg(p) = 0.04_kind_lake
1271 ! end if
1272
1273 if (t_grnd(c) >= tfrz) then ! for unfrozen lake
1274 z0mg(p) = 0.001_kind_lake !original 0.01
1275 else if(snl(c) == 0 ) then ! for frozen lake
1276 ! Is this okay even if it is snow covered? What is the roughness over
1277 ! non-veg. snow?
1278 z0mg(p) = 0.005_kind_lake !original 0.04, now for frozen lake without snow
1279 else ! for frozen lake with snow
1280 z0mg(p) = 0.0024_kind_lake
1281 end if
1282
1283 if(.false.) then
1284 ! This can't work since it uses ustar before ustar is initialized
1285 !- tgs - use COARE formulation for z0hg and z0qg.
1286 !-- suggestion from Ayumi Manome (GLERL), Aug. 2018
1287 !-- Charusombat et al., 2018, https://doi.org/10.5194/hess-2017-725
1288 tc=forc_t(g)-273.15_kind_lake
1289 visc=1.326e-5_kind_lake*(1._kind_lake + 6.542e-3_kind_lake*tc + 8.301e-6_kind_lake*tc*tc &
1290 - 4.84e-9_kind_lake*tc*tc*tc)
1291 visc=max(1e-7_kind_lake, visc)
1292
1293 ren = max(ustar(p)*z0mg(p)/visc, 0.1_kind_lake)
1294 z0hg(p) = (5.5e-5_kind_lake)*(ren**(-0.60_kind_lake))
1295
1296 z0hg(p) = min(z0hg(p),1.0e-4_kind_lake)
1297 z0hg(p) = max(z0hg(p),2.0e-9_kind_lake)
1298
1299 z0qg(p) = z0hg(p)
1300
1301 ! end COARE
1302 endif
1303 z0hg(p) = z0mg(p)
1304 z0qg(p) = z0mg(p)
1305
1306 ! Latent heat
1307
1308 if(pergro) then
1309 htvp(c) = hvap
1310 else
1311 if (t_grnd(c) > tfrz) then
1312 htvp(c) = hvap
1313 else
1314 htvp(c) = hsub
1315 end if
1316 endif
1317
1318 ! Zack Subin, 3/26/09: Shouldn't this be the ground temperature rather than the air temperature above?
1319 ! I'll change it for now.
1320
1321 ! Initialize stability variables
1322
1323 ur(p) = max(1.0_kind_lake,sqrt(forc_u(g)*forc_u(g)+forc_v(g)*forc_v(g)))
1324 dth(p) = thm(c)-t_grnd(c)
1325 dqh(p) = forc_q(g)-qsatg(c)
1326 dthv = dth(p)*(1._kind_lake+con_fvirt*forc_q(g))+con_fvirt*forc_th(g)*dqh(p)
1327 zldis(p) = forc_hgt_u(g) - 0._kind_lake
1328
1329 ! Initialize Monin-Obukhov length and wind speed
1330
1331 call moninobukini(ur(p), thv(c), dthv, zldis(p), z0mg(p), um(p), obu(p))
1332
1333 end do
1334
1335 iter = 1
1336 fncopy = num_shlakep
1337 fpcopy(1:num_shlakep) = filter_shlakep(1:num_shlakep)
1338
1339 ! Begin stability iteration
1340
1341 iteration : do while (iter <= niters .and. fncopy > 0)
1342
1343 ! Determine friction velocity, and potential temperature and humidity
1344 ! profiles of the surface boundary layer
1345
1346 call frictionvelocity(pgridcell,forc_hgt,forc_hgt_u, & !i
1347 forc_hgt_t,forc_hgt_q, & !i
1348 lbp, ubp, fncopy, fpcopy, & !i
1349 displa, z0mg, z0hg, z0qg, & !i
1350 obu, iter, ur, um, & !i
1351 ustar,temp1, temp2, temp12m, temp22m, & !o
1352 u10,fv, & !o
1353 fm) !i&o
1354
1355 !dir$ concurrent
1356 !cdir nodep
1357 do fp = 1, fncopy
1358 p = fpcopy(fp)
1359 c = pcolumn(p)
1360 g = pgridcell(p)
1361
1362 tgbef(c) = t_grnd(c)
1363 if (t_grnd(c) > tfrz .and. t_lake(c,1) > tfrz .and. snl(c) == 0) then
1364 tksur = savedtke1(c)
1365 ! Set this to the eddy conductivity from the last
1366 ! timestep, as the molecular conductivity will be orders of magnitude too small.
1367 ! Will have to deal with first timestep.
1368 tsur = t_lake(c,1)
1369 else if (snl(c) == 0) then !frozen but no snow layers
1370 tksur = tkice
1371 tsur = t_lake(c,1)
1372 else
1373 !Need to calculate thermal conductivity of the top snow layer
1374 bw = (h2osoi_ice(c,jtop(c))+h2osoi_liq(c,jtop(c)))/dz(c,jtop(c))
1375 tksur = tkairc + (7.75e-5_kind_lake *bw + 1.105e-6_kind_lake*bw*bw)*(tkice-tkairc)
1376 tsur = t_soisno(c,jtop(c))
1377 end if
1378
1379 ! Determine aerodynamic resistances
1380
1381 ram(p) = 1._kind_lake/(ustar(p)*ustar(p)/um(p))
1382 rah(p) = 1._kind_lake/(temp1(p)*ustar(p))
1383 raw(p) = 1._kind_lake/(temp2(p)*ustar(p))
1384 ram1(p) = ram(p) !pass value to global variable
1385
1386 ! Get derivative of fluxes with respect to ground temperature
1387
1388 stftg3(p) = emg*sb*tgbef(c)*tgbef(c)*tgbef(c)
1389
1390 ! Changed surface temperature from t_lake(c,1) to tsur.
1391 ! Also adjusted so that if there are snow layers present, all radiation is absorbed in the top layer.
1392 ax = betaprime(c)*sabg(p) + emg*forc_lwrad(g) + 3._kind_lake*stftg3(p)*tgbef(c) &
1393 + forc_rho(g)*cpair/rah(p)*thm(c) &
1394 - htvp(c)*forc_rho(g)/raw(p)*(qsatg(c)-qsatgdt(c)*tgbef(c) - forc_q(g)) &
1395 + tksur*tsur/dzsur(c)
1396 !Changed sabg(p) and to betaprime(c)*sabg(p).
1397 bx = 4._kind_lake*stftg3(p) + forc_rho(g)*cpair/rah(p) &
1398 + htvp(c)*forc_rho(g)/raw(p)*qsatgdt(c) + tksur/dzsur(c)
1399
1400 t_grnd(c) = ax/bx
1401
1402 ! Update htvp
1403 if(.not.pergro) then
1404 if (t_grnd(c) > tfrz) then
1405 htvp(c) = hvap
1406 else
1407 htvp(c) = hsub
1408 end if
1409 endif
1410
1411 ! Surface fluxes of momentum, sensible and latent heat
1412 ! using ground temperatures from previous time step
1413
1414 eflx_sh_grnd(p) = forc_rho(g)*cpair*(t_grnd(c)-thm(c))/rah(p)
1415 qflx_evap_soi(p) = forc_rho(g)*(qsatg(c)+qsatgdt(c)*(t_grnd(c)-tgbef(c))-forc_q(g))/raw(p)
1416
1417 ! Re-calculate saturated vapor pressure, specific humidity and their
1418 ! derivatives at lake surface
1419
1420 call qsat(t_grnd(c), forc_pbot(g), eg, degdt, qsatg(c), qsatgdt(c))
1421
1422 dth(p)=thm(c)-t_grnd(c)
1423 dqh(p)=forc_q(g)-qsatg(c)
1424
1425 tstar = temp1(p)*dth(p)
1426 qstar = temp2(p)*dqh(p)
1427
1428 thvstar=tstar*(1._kind_lake+con_fvirt*forc_q(g)) + con_fvirt*forc_th(g)*qstar
1429 zeta=zldis(p)*vkc * grav*thvstar/(ustar(p)**2*thv(c))
1430
1431 if (zeta >= 0._kind_lake) then !stable
1432 zeta = min(2._kind_lake,max(zeta,0.01_kind_lake))
1433 um(p) = max(ur(p),0.1_kind_lake)
1434 else !unstable
1435 zeta = max(-100._kind_lake,min(zeta,-0.01_kind_lake))
1436 wc = beta1*(-grav*ustar(p)*thvstar*zii/thv(c))**0.333_kind_lake
1437 um(p) = sqrt(ur(p)*ur(p)+wc*wc)
1438 end if
1439 obu(p) = zldis(p)/zeta
1440
1441 if (obuold(p)*obu(p) < 0._kind_lake) nmozsgn(p) = nmozsgn(p)+1
1442
1443 obuold(p) = obu(p)
1444
1445 end do ! end of filtered pft loop
1446
1447 iter = iter + 1
1448 if (iter <= niters ) then
1449 ! Rebuild copy of pft filter for next pass through the ITERATION loop
1450
1451 fnold = fncopy
1452 fncopy = 0
1453 do fp = 1, fnold
1454 p = fpcopy(fp)
1455 if (nmozsgn(p) < 3) then
1456 fncopy = fncopy + 1
1457 fpcopy(fncopy) = p
1458 end if
1459 end do ! end of filtered pft loop
1460 end if
1461
1462 end do iteration ! end of stability iteration
1463
1464 !dir$ concurrent
1465 !cdir nodep
1466 do fp = 1, num_shlakep
1467 p = filter_shlakep(fp)
1468 c = pcolumn(p)
1469 g = pgridcell(p)
1470
1471 ! If there is snow on the ground and t_grnd > tfrz: reset t_grnd = tfrz.
1472 ! Re-evaluate ground fluxes.
1473 ! h2osno > 0.5 prevents spurious fluxes.
1474 ! note that qsatg and qsatgdT should be f(tgbef) (PET: not sure what this
1475 ! comment means)
1476 ! Zack Subin, 3/27/09: Since they are now a function of whatever t_grnd was before cooling
1477 ! to freezing temperature, then this value should be used in the derivative correction term.
1478 ! Should this happen if the lake temperature is below freezing, too? I'll assume that for now.
1479 ! Also, allow convection if ground temp is colder than lake but warmer than 4C, or warmer than
1480 ! lake which is warmer than freezing but less than 4C.
1481 if ( (h2osno(c) > 0.5_kind_lake .or. t_lake(c,1) <= tfrz) .and. t_grnd(c) > tfrz) then
1482 t_grnd_temp = t_grnd(c)
1483 t_grnd(c) = tfrz
1484 eflx_sh_grnd(p) = forc_rho(g)*cpair*(t_grnd(c)-thm(c))/rah(p)
1485 qflx_evap_soi(p) = forc_rho(g)*(qsatg(c)+qsatgdt(c)*(t_grnd(c)-t_grnd_temp) - forc_q(g))/raw(p)
1486 else if ( (t_lake(c,1) > t_grnd(c) .and. t_grnd(c) > tdmax) .or. &
1487 (t_lake(c,1) < t_grnd(c) .and. t_lake(c,1) > tfrz .and. t_grnd(c) < tdmax) ) then
1488 ! Convective mixing will occur at surface
1489 t_grnd_temp = t_grnd(c)
1490 t_grnd(c) = t_lake(c,1)
1491 eflx_sh_grnd(p) = forc_rho(g)*cpair*(t_grnd(c)-thm(c))/rah(p)
1492 qflx_evap_soi(p) = forc_rho(g)*(qsatg(c)+qsatgdt(c)*(t_grnd(c)-t_grnd_temp) - forc_q(g))/raw(p)
1493 end if
1494
1495 ! Update htvp
1496 if(.not.pergro) then
1497 if (t_grnd(c) > tfrz) then
1498 htvp(c) = hvap
1499 else
1500 htvp(c) = hsub
1501 end if
1502 endif
1503
1504 ! Net longwave from ground to atmosphere
1505
1506 ! eflx_lwrad_out(p) = (1._kind_lake-emg)*forc_lwrad(g) + stftg3(p)*(-3._kind_lake*tgbef(c)+4._kind_lake*t_grnd(c))
1507 ! What is tgbef doing in this equation? Can't it be exact now? --Zack Subin, 4/14/09
1508 eflx_lwrad_out(p) = (1._kind_lake-emg)*forc_lwrad(g) + emg*sb*t_grnd(c)**4
1509
1510 ! Ground heat flux
1511
1512 eflx_soil_grnd(p) = sabg(p) + forc_lwrad(g) - eflx_lwrad_out(p) - &
1513 eflx_sh_grnd(p) - htvp(c)*qflx_evap_soi(p)
1514 !Why is this sabg(p) and not beta*sabg(p)??
1515 !I've kept this as the incorrect sabg so that the energy balance check will be correct.
1516 !This is the effective energy flux into the ground including the lake solar absorption
1517 !below the surface. The variable eflx_gnet will be used to pass the actual heat flux
1518 !from the ground interface into the lake.
1519
1520 taux(p) = -forc_rho(g)*forc_u(g)/ram(p)
1521 tauy(p) = -forc_rho(g)*forc_v(g)/ram(p)
1522
1523 eflx_sh_tot(p) = eflx_sh_grnd(p)
1524 qflx_evap_tot(p) = qflx_evap_soi(p)
1525 eflx_lh_tot(p) = htvp(c)*qflx_evap_soi(p)
1526 eflx_lh_grnd(p) = htvp(c)*qflx_evap_soi(p)
1527 if(debug_print) then
15281604 format('CLM_Lake ShalLakeFluxes: c=',i0,' sensible heat = ',f12.4,' latent heat =',f12.4, &
1529 ' ground temp = ', f12.4, ' h2osno = ', f12.4, ' at xlat_d=',f10.3,' xlon_d=',f10.3)
1530 print 1604, c, eflx_sh_tot(p), eflx_lh_tot(p), t_grnd(c), h2osno(c),xlat_d,xlon_d
1531 if (abs(eflx_sh_tot(p)) > 1500 .or. abs(eflx_lh_tot(p)) > 1500) then
15323018 format('CLM_Lake ShalLakeFluxes: WARNING: SH=',f12.4,' LH=',f12.4,' at xlat_d=',f10.3,' xlon_d=',f10.3)
1533 print 3018,eflx_sh_tot(p), eflx_lh_tot(p),xlat_d,xlon_d
1534 end if
1535 if (abs(eflx_sh_tot(p)) > 10000 .or. abs(eflx_lh_tot(p)) > 10000 &
1536 .or. abs(t_grnd(c)-288)>200 ) then
1537840 format('CLM_Lake ShalLakeFluxes: t_grnd is out of range: eflx_sh_tot(p)=',g20.12,' eflx_lh_tot(p)=',g20.12,' t_grnd(c)=',g20.12,' at p=',i0,' c=',i0,' xlat_d=',f10.3,' xlon_d=',f10.3)
1538 write(message,840) eflx_sh_tot(p),eflx_lh_tot(p),t_grnd(c),p,c,xlat_d,xlon_d
1539 ! errmsg=message
1540 ! errflg=1
1541 write(0,'(A)') trim(message)
1542 endif
1543 endif
1544 ! 2 m height air temperature
1545 t_ref2m(p) = thm(c) + temp1(p)*dth(p)*(1._kind_lake/temp12m(p) - 1._kind_lake/temp1(p))
1546
1547 ! 2 m height specific humidity
1548 q_ref2m(p) = forc_q(g) + temp2(p)*dqh(p)*(1._kind_lake/temp22m(p) - 1._kind_lake/temp2(p))
1549
1550 ! Energy residual used for melting snow
1551 ! Effectively moved to ShalLakeTemp
1552
1553 ! Prepare for lake layer temperature calculations below
1554 ! fin(c) = betaprime * sabg(p) + forc_lwrad(g) - (eflx_lwrad_out(p) + &
1555 ! eflx_sh_tot(p) + eflx_lh_tot(p))
1556 ! NOW this is just the net ground heat flux calculated below.
1557
1558 eflx_gnet(p) = betaprime(c) * sabg(p) + forc_lwrad(g) - (eflx_lwrad_out(p) + &
1559 eflx_sh_tot(p) + eflx_lh_tot(p))
1560 ! This is the actual heat flux from the ground interface into the lake, not including
1561 ! the light that penetrates the surface.
1562
1563 ! u2m = max(1.0_kind_lake,ustar(p)/vkc*log(2._kind_lake/z0mg(p)))
1564 ! u2 often goes below 1 m/s; it seems like the only reason for this minimum is to
1565 ! keep it from being zero in the ks equation below; 0.1 m/s is a better limit for
1566 ! stable conditions --ZS
1567 u2m = max(0.1_kind_lake,ustar(p)/vkc*log(2._kind_lake/z0mg(p)))
1568
1569 ws(c) = 1.2e-03_kind_lake * u2m
1570 ks(c) = 6.6_kind_lake*sqrt(abs(sin(lat(g))))*(u2m**(-1.84_kind_lake))
1571
1572 end do
1573
1574 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1575 ! End of surface flux relevant code in original BiogeophysicsLakeMod until history loop.
1576
1577 ! The following are needed for global average on history tape.
1578
1579 !dir$ concurrent
1580 !cdir nodep
1581 do fp = 1, num_shlakep
1582 p = filter_shlakep(fp)
1583 c = pcolumn(p)
1584 g = pgridcell(p)
1585 ! t_veg(p) = forc_t(g)
1586 !This is an odd choice, since elsewhere t_veg = t_grnd for bare ground.
1587 !Zack Subin, 4/09
1588 t_veg(p) = t_grnd(c)
1589 eflx_lwrad_net(p) = eflx_lwrad_out(p) - forc_lwrad(g)
1590 qflx_prec_grnd(p) = forc_rain(g) + forc_snow(g)
1591 end do
1592
1593 ustar_out(1) = ustar(1)
1594
1595
1596END SUBROUTINE shallakefluxes
1597
1598SUBROUTINE shallaketemperature(t_grnd,h2osno,sabg,dz,dz_lake,z,zi, & !i
1599 z_lake,ws,ks,snl,eflx_gnet,lakedepth, &
1600 lake_icefrac,snowdp, & !i&o
1601 eflx_sh_grnd,eflx_sh_tot,eflx_soil_grnd, & !o
1602 t_lake,t_soisno,h2osoi_liq, &
1603 h2osoi_ice,savedtke1, &
1604 watsat, tksatu, tkmg, tkdry, csol, dtime, &
1605 frac_iceold,qflx_snomelt,imelt,errmsg,errflg,xlat_d,xlon_d)
1606 !=======================================================================================================
1607 ! DESCRIPTION:
1608
1623 !
1624 ! Eddy + molecular diffusion:
1625 ! d ts d d ts 1 ds
1626 ! ---- = -- [(km + ke) ----] + -- --
1627 ! dt dz dz cw dz
1628 !
1629 ! where: ts = temperature (kelvin)
1630 ! t = time (s)
1631 ! z = depth (m)
1632 ! km = molecular diffusion coefficient (m**2/s)
1633 ! ke = eddy diffusion coefficient (m**2/s)
1634 ! cw = heat capacity (j/m**3/kelvin)
1635 ! s = heat source term (w/m**2)
1636 !
1637 ! Shallow lakes are allowed to have variable depth, set in _____.
1638 !
1639 ! For shallow lakes: ke > 0 if unfrozen,
1640 ! and convective mixing occurs WHETHER OR NOT frozen. (See e.g. Martynov...)
1641 !
1642 ! Use the Crank-Nicholson method to set up tridiagonal system of equations to
1643 ! solve for ts at time n+1, where the temperature equation for layer i is
1644 ! r_i = a_i [ts_i-1] n+1 + b_i [ts_i] n+1 + c_i [ts_i+1] n+1
1645 !
1646 ! The solution conserves energy as:
1647 !
1648 ! [For lake layers]
1649 ! cw*([ts( 1)] n+1 - [ts( 1)] n)*dz( 1)/dt + ... +
1650 ! cw*([ts(nlevlake)] n+1 - [ts(nlevlake)] n)*dz(nlevlake)/dt = fin
1651 ! But now there is phase change, so cv is not constant and there is
1652 ! latent heat.
1653 !
1654 ! where:
1655 ! [ts] n = old temperature (kelvin)
1656 ! [ts] n+1 = new temperature (kelvin)
1657 ! fin = heat flux into lake (w/m**2)
1658 ! = betaprime*sabg + forc_lwrad - eflx_lwrad_out - eflx_sh_tot - eflx_lh_tot
1659 ! (This is now the same as the ground heat flux.)
1660 ! + phi(1) + ... + phi(nlevlake) + phi(top soil level)
1661 ! betaprime = beta(islak) for no snow layers, and 1 for snow layers.
1662 ! This assumes all radiation is absorbed in the top snow layer and will need
1663 ! to be changed for CLM 4.
1664 !
1665 ! WARNING: This subroutine assumes lake columns have one and only one pft.
1666 !
1667 ! Outline:
1668 ! 1!) Initialization
1669 ! 2!) Lake density
1670 ! 3!) Diffusivity
1671 ! 4!) Heat source term from solar radiation penetrating lake
1672 ! 5!) Set thermal props and find initial energy content
1673 ! 6!) Set up vectors for tridiagonal matrix solution
1674 ! 7!) Solve tridiagonal and back-substitute
1675 ! 8!) (Optional) Do first energy check using temperature change at constant heat capacity.
1676 ! 9!) Phase change
1677 ! 9.5!) (Optional) Do second energy check using temperature change and latent heat, considering changed heat capacity.
1678 ! Also do soil water balance check.
1679 !10!) Convective mixing
1680 !11!) Do final energy check to detect small numerical errors (especially from convection)
1681 ! and dump small imbalance into sensible heat, or pass large errors to BalanceCheckMod for abort.
1682 !
1683 ! REVISION HISTORY:
1684 ! Created by Zack Subin, 2009.
1685 ! Reedited by Hongping Gu, 2010.
1686 ! Updated for CCPP by Sam Trahan, 2022.
1687 !=========================================================================================================
1688
1689
1690 implicit none
1691
1692 !in:
1693 real(kind_lake), intent(in) :: xlat_d
1694 real(kind_lake), intent(in) :: xlon_d
1695 integer, intent(inout) :: errflg
1696 real(kind_lake), intent(in) :: watsat(1,nlevsoil)
1697 real(kind_lake), intent(in) :: tksatu(1,nlevsoil)
1698 real(kind_lake), intent(in) :: tkmg(1,nlevsoil)
1699 real(kind_lake), intent(in) :: tkdry(1,nlevsoil)
1700 real(kind_lake), intent(in) :: csol(1,nlevsoil)
1701 character(*), intent(inout) :: errmsg
1702 real(kind_lake), intent(in) :: t_grnd(1)
1703 real(kind_lake), intent(inout) :: h2osno(1)
1704 real(kind_lake), intent(in) :: sabg(1)
1705 real(kind_lake), intent(in) :: dz(1,-nlevsnow + 1:nlevsoil)
1706 real(kind_lake), intent(in) :: dz_lake(1,nlevlake)
1707 real(kind_lake), intent(in) :: z(1,-nlevsnow+1:nlevsoil)
1708 real(kind_lake), intent(in) :: zi(1,-nlevsnow+0:nlevsoil)
1710 real(kind_lake), intent(in) :: z_lake(1,nlevlake)
1711 real(kind_lake), intent(in) :: ws(1)
1712 real(kind_lake), intent(in) :: ks(1)
1714 integer , intent(in) :: snl(1)
1715 real(kind_lake), intent(inout) :: eflx_gnet(1)
1716 real(kind_lake), intent(in) :: lakedepth(1)
1717
1718 ! real(kind_lake), intent(in) :: watsat(1,nlevsoil) ! volumetric soil water at saturation (porosity)
1719 real(kind_lake), intent(inout) :: snowdp(1)
1720 real(kind_lake), intent(in) :: dtime
1721 !out:
1722
1723 real(kind_lake), intent(out) :: eflx_sh_grnd(1)
1724 real(kind_lake), intent(out) :: eflx_sh_tot(1)
1725 real(kind_lake), intent(out) :: eflx_soil_grnd(1)
1727 !real(kind_lake), intent(out) :: qmelt(1) !< snow melt [mm/s] [temporary]
1728
1729 real(kind_lake), intent(inout) :: t_lake(1,nlevlake)
1730 real(kind_lake), intent(inout) :: t_soisno(1,-nlevsnow+1:nlevsoil)
1731 real(kind_lake), intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil)
1732 real(kind_lake), intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil)
1733 real(kind_lake), intent(inout) :: lake_icefrac(1,nlevlake)
1734 real(kind_lake), intent(out) :: savedtke1(1)
1735 real(kind_lake), intent(out) :: frac_iceold(1,-nlevsnow+1:nlevsoil)
1736 real(kind_lake), intent(out) :: qflx_snomelt(1)
1737 integer, intent(out) :: imelt(1,-nlevsnow+1:nlevsoil)
1738
1739
1740 ! OTHER LOCAL VARIABLES:
1741
1742 integer , parameter :: islak = 2 ! index of lake, 1 = deep lake, 2 = shallow lake
1743 real(kind_lake), parameter :: p0 = 1._kind_lake ! neutral value of turbulent prandtl number
1744 integer :: i,j,fc,fp,g,c,p ! do loop or array index
1745 real(kind_lake) :: eta(2) ! light extinction coefficient (/m): depends on lake type
1746 real(kind_lake) :: cwat ! specific heat capacity of water (j/m**3/kelvin)
1747 real(kind_lake) :: cice_eff ! effective heat capacity of ice (using density of
1748 ! water because layer depth is not adjusted when freezing
1749 real(kind_lake) :: cfus ! effective heat of fusion per unit volume
1750 ! using water density as above
1751 real(kind_lake) :: km ! molecular diffusion coefficient (m**2/s)
1752 real(kind_lake) :: tkice_eff ! effective conductivity since layer depth is constant
1753 real(kind_lake) :: a(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) ! "a" vector for tridiagonal matrix
1754 real(kind_lake) :: b(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) ! "b" vector for tridiagonal matrix
1755 real(kind_lake) :: c1(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) ! "c" vector for tridiagonal matrix
1756 real(kind_lake) :: r(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) ! "r" vector for tridiagonal solution
1757 real(kind_lake) :: rhow(lbc:ubc,nlevlake) ! density of water (kg/m**3)
1758 real(kind_lake) :: phi(lbc:ubc,nlevlake) ! solar radiation absorbed by layer (w/m**2)
1759 real(kind_lake) :: kme(lbc:ubc,nlevlake) ! molecular + eddy diffusion coefficient (m**2/s)
1760 real(kind_lake) :: rsfin ! relative flux of solar radiation into layer
1761 real(kind_lake) :: rsfout ! relative flux of solar radiation out of layer
1762 real(kind_lake) :: phi_soil(lbc:ubc) ! solar radiation into top soil layer (W/m**2)
1763 real(kind_lake) :: ri ! richardson number
1764 real(kind_lake) :: fin(lbc:ubc) ! net heat flux into lake at ground interface (w/m**2)
1765 real(kind_lake) :: ocvts(lbc:ubc) ! (cwat*(t_lake[n ])*dz
1766 real(kind_lake) :: ncvts(lbc:ubc) ! (cwat*(t_lake[n+1])*dz
1767 real(kind_lake) :: ke ! eddy diffusion coefficient (m**2/s)
1768 real(kind_lake) :: zin ! depth at top of layer (m)
1769 real(kind_lake) :: zout ! depth at bottom of layer (m)
1770 real(kind_lake) :: drhodz ! d [rhow] /dz (kg/m**4)
1771 real(kind_lake) :: n2 ! brunt-vaisala frequency (/s**2)
1772 real(kind_lake) :: num ! used in calculating ri
1773 real(kind_lake) :: den ! used in calculating ri
1774 real(kind_lake) :: tav_froz(lbc:ubc) ! used in aver temp for convectively mixed layers (C)
1775 real(kind_lake) :: tav_unfr(lbc:ubc) ! "
1776 real(kind_lake) :: nav(lbc:ubc) ! used in aver temp for convectively mixed layers
1777 real(kind_lake) :: phidum ! temporary value of phi
1778 real(kind_lake) :: iceav(lbc:ubc) ! used in calc aver ice for convectively mixed layers
1779 real(kind_lake) :: qav(lbc:ubc) ! used in calc aver heat content for conv. mixed layers
1780 integer :: jtop(lbc:ubc) ! top level for each column (no longer all 1)
1781 real(kind_lake) :: cv (lbc:ubc,-nlevsnow+1:nlevsoil) !heat capacity of soil/snow [J/(m2 K)]
1782 real(kind_lake) :: tk (lbc:ubc,-nlevsnow+1:nlevsoil) !thermal conductivity of soil/snow [W/(m K)]
1783 !(at interface below, except for j=0)
1784 real(kind_lake) :: cv_lake (lbc:ubc,1:nlevlake) !heat capacity [J/(m2 K)]
1785 real(kind_lake) :: tk_lake (lbc:ubc,1:nlevlake) !thermal conductivity at layer node [W/(m K)]
1786 real(kind_lake) :: cvx (lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !heat capacity for whole column [J/(m2 K)]
1787 real(kind_lake) :: tkix(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !thermal conductivity at layer interfaces
1788 !for whole column [W/(m K)]
1789 real(kind_lake) :: tx(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) ! temperature of whole column [K]
1790 real(kind_lake) :: tktopsoillay(lbc:ubc) ! thermal conductivity [W/(m K)]
1791 real(kind_lake) :: fnx(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !heat diffusion through the layer interface below [W/m2]
1792 real(kind_lake) :: phix(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !solar source term for whole column [W/m**2]
1793 real(kind_lake) :: zx(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !interface depth (+ below surface) for whole column [m]
1794 real(kind_lake) :: dzm !used in computing tridiagonal matrix [m]
1795 real(kind_lake) :: dzp !used in computing tridiagonal matrix [m]
1796 integer :: jprime ! j - nlevlake
1797 real(kind_lake) :: factx(lbc:ubc,-nlevsnow+1:nlevlake+nlevsoil) !coefficient used in computing tridiagonal matrix
1798 real(kind_lake) :: t_lake_bef(lbc:ubc,1:nlevlake) !beginning lake temp for energy conservation check [K]
1799 real(kind_lake) :: t_soisno_bef(lbc:ubc,-nlevsnow+1:nlevsoil) !beginning soil temp for E cons. check [K]
1800 real(kind_lake) :: lhabs(lbc:ubc) ! total per-column latent heat abs. from phase change (J/m^2)
1801 real(kind_lake) :: esum1(lbc:ubc) ! temp for checking energy (J/m^2)
1802 real(kind_lake) :: esum2(lbc:ubc) ! ""
1803 real(kind_lake) :: zsum(lbc:ubc) ! temp for putting ice at the top during convection (m)
1804 real(kind_lake) :: wsum(lbc:ubc) ! temp for checking water (kg/m^2)
1805 real(kind_lake) :: wsum_end(lbc:ubc) ! temp for checking water (kg/m^2)
1806 real(kind_lake) :: errsoi(1) ! soil/lake energy conservation error (W/m**2)
1807 real(kind_lake) :: eflx_snomelt(1) !snow melt heat flux (W/m**2)
1808 CHARACTER*256 :: message
1809 !
1810 ! Constants for lake temperature model
1811 !
1812 real(kind_lake), parameter :: beta(2) = & ! fraction solar rad absorbed at surface: depends on lake type
1813 (/0.4_kind_lake, 0.4_kind_lake/) ! (deep lake, shallow lake)
1814 real(kind_lake), parameter :: za(2) = & ! base of surface absorption layer (m): depends on lake type
1815 (/0.6_kind_lake, 0.6_kind_lake/)
1816 ! For now, keep beta and za for shallow lake the same as deep lake, until better data is found.
1817 ! It looks like eta is key and that larger values give better results for shallow lakes. Use
1818 ! empirical expression from Hakanson (below). This is still a very unconstrained parameter
1819 ! that deserves more attention.
1820 ! Some radiation will be allowed to reach the soil.
1821 !-----------------------------------------------------------------------
1822
1823
1824 ! 1!) Initialization
1825 ! Determine step size
1826
1827 ! Initialize constants
1828 cwat = cpliq*denh2o ! water heat capacity per unit volume
1829 cice_eff = cpice*denh2o !use water density because layer depth is not adjusted
1830 !for freezing
1831 cfus = hfus*denh2o ! latent heat per unit volume
1832 tkice_eff = tkice * denice/denh2o !effective conductivity since layer depth is constant
1833 km = tkwat/cwat ! a constant (molecular diffusivity)
1834
1835 ! Begin calculations
1836
1837 !dir$ concurrent
1838 !cdir nodep
1839 do fc = 1, num_shlakec
1840 c = filter_shlakec(fc)
1841
1842 ! Initialize Ebal quantities computed below
1843
1844 ocvts(c) = 0._kind_lake
1845 ncvts(c) = 0._kind_lake
1846 esum1(c) = 0._kind_lake
1847 esum2(c) = 0._kind_lake
1848
1849 end do
1850
1851 ! Initialize set of previous time-step variables as in DriverInit,
1852 ! which is currently not called over lakes. This has to be done
1853 ! here because phase change will occur in this routine.
1854 ! Ice fraction of snow at previous time step
1855
1856 do j = -nlevsnow+1,0
1857 !dir$ concurrent
1858 !cdir nodep
1859 do fc = 1, num_shlakec
1860 c = filter_shlakec(fc)
1861 if (j >= snl(c) + 1) then
1862 frac_iceold(c,j) = h2osoi_ice(c,j)/(h2osoi_liq(c,j)+h2osoi_ice(c,j))
1863 end if
1864 end do
1865 end do
1866
1867 ! Sum soil water.
1868 do j = 1, nlevsoil
1869 !dir$ concurrent
1870 !cdir nodep
1871 do fc = 1, num_shlakec
1872 c = filter_shlakec(fc)
1873 if (j == 1) wsum(c) = 0._kind_lake
1874 wsum(c) = wsum(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
1875 end do
1876 end do
1877
1878 !dir$ concurrent
1879 !cdir nodep
1880 do fp = 1, num_shlakep
1881 p = filter_shlakep(fp)
1882 c = pcolumn(p)
1883
1884
1885 ! Prepare for lake layer temperature calculations below
1886
1887 ! fin(c) = betaprime * sabg(p) + forc_lwrad(g) - (eflx_lwrad_out(p) + &
1888 ! eflx_sh_tot(p) + eflx_lh_tot(p))
1889 ! fin(c) now passed from ShalLakeFluxes as eflx_gnet
1890 fin(c) = eflx_gnet(p)
1891
1892 end do
1893
1894 ! 2!) Lake density
1895
1896 do j = 1, nlevlake
1897 !dir$ concurrent
1898 !cdir nodep
1899 do fc = 1, num_shlakec
1900 c = filter_shlakec(fc)
1901 rhow(c,j) = (1._kind_lake - lake_icefrac(c,j)) * &
1902 1000._kind_lake*( 1.0_kind_lake - 1.9549e-05_kind_lake*(abs(t_lake(c,j)-277._kind_lake))**1.68_kind_lake ) &
1903 + lake_icefrac(c,j)*denice
1904 ! Allow for ice fraction; assume constant ice density.
1905 ! Is this the right weighted average?
1906 ! Using this average will make sure that surface ice is treated properly during
1907 ! convective mixing.
1908 end do
1909 end do
1910
1911 ! 3!) Diffusivity and implied thermal "conductivity" = diffusivity * cwat
1912 do j = 1, nlevlake-1
1913 !dir$ prefervector
1914 !dir$ concurrent
1915 !cdir nodep
1916 do fc = 1, num_shlakec
1917 c = filter_shlakec(fc)
1918 drhodz = (rhow(c,j+1)-rhow(c,j)) / (z_lake(c,j+1)-z_lake(c,j))
1919 n2 = grav / rhow(c,j) * drhodz
1920 ! Fixed sign error here: our z goes up going down into the lake, so no negative
1921 ! sign is needed to make this positive unlike in Hostetler. --ZS
1922 num = 40._kind_lake * n2 * (vkc*z_lake(c,j))**2
1923 den = max( (ws(c)**2) * exp(-2._kind_lake*ks(c)*z_lake(c,j)), 1.e-10_kind_lake )
1924 ri = ( -1._kind_lake + sqrt( max(1._kind_lake+num/den, 0._kind_lake) ) ) / 20._kind_lake
1925 if (t_grnd(c) > tfrz .and. t_lake(c,1) > tfrz .and. snl(c) == 0) then
1926 ! ke = vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._kind_lake+37._kind_lake*ri*ri)
1927
1928 if( t_lake(c,1) > 277.15_kind_lake ) then
1929 if (lakedepth(c) > 15.0 ) then
1930 ke = 1.e+2_kind_lake*vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._kind_lake+37._kind_lake*ri*ri)
1931 else
1932 ke = vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._kind_lake+37._kind_lake*ri*ri)
1933 endif
1934 else
1935 if (lakedepth(c) > 15.0 ) then
1936 if (lakedepth(c) > 150.0 ) then
1937 ke = 1.e+5_kind_lake*vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._kind_lake+37._kind_lake*ri*ri)
1938 else
1939 ke =1.e+4_kind_lake*vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._kind_lake+37._kind_lake*ri*ri)
1940 end if
1941 else
1942 ke = vkc*ws(c)*z_lake(c,j)/p0 * exp(-ks(c)*z_lake(c,j)) / (1._kind_lake+37._kind_lake*ri*ri)
1943 endif
1944 end if
1945
1946 kme(c,j) = km + ke
1947 tk_lake(c,j) = kme(c,j)*cwat
1948 ! If there is some ice in this layer (this should rarely happen because the surface
1949 ! is unfrozen and it will be unstable), still use the cwat to get out the tk b/c the eddy
1950 ! diffusivity equation assumes water.
1951 else
1952 kme(c,j) = km
1953 tk_lake(c,j) = tkwat*tkice_eff / ( (1._kind_lake-lake_icefrac(c,j))*tkice_eff &
1954 + tkwat*lake_icefrac(c,j) )
1955 ! Assume the resistances add as for the calculation of conductivities at layer interfaces.
1956 end if
1957 end do
1958 end do
1959
1960 !dir$ concurrent
1961 !cdir nodep
1962 do fc = 1, num_shlakec
1963 c = filter_shlakec(fc)
1964
1965 j = nlevlake
1966 kme(c,nlevlake) = kme(c,nlevlake-1)
1967
1968 if (t_grnd(c) > tfrz .and. t_lake(c,1) > tfrz .and. snl(c) == 0) then
1969 tk_lake(c,j) = tk_lake(c,j-1)
1970 else
1971 tk_lake(c,j) = tkwat*tkice_eff / ( (1._kind_lake-lake_icefrac(c,j))*tkice_eff &
1972 + tkwat*lake_icefrac(c,j) )
1973 end if
1974
1975 ! Use in surface flux calculation for next timestep.
1976 savedtke1(c) = kme(c,1)*cwat ! Will only be used if unfrozen
1977 ! set number of column levels for use by Tridiagonal below
1978 jtop(c) = snl(c) + 1
1979 end do
1980
1981 ! 4!) Heat source term: unfrozen lakes only
1982 do j = 1, nlevlake
1983 !dir$ concurrent
1984 !cdir nodep
1985 do fp = 1, num_shlakep
1986 p = filter_shlakep(fp)
1987 c = pcolumn(p)
1988
1989 ! Set eta(:), the extinction coefficient, according to L Hakanson, Aquatic Sciences, 1995
1990 ! (regression of Secchi Depth with lake depth for small glacial basin lakes), and the
1991 ! Poole & Atkins expression for extinction coeffient of 1.7 / Secchi Depth (m).
1992 if(.not.use_etalake) then
1993 eta(:) = 1.1925_kind_lake*lakedepth(c)**(-0.424)
1994 else
1995 eta(:) = etalake
1996 endif
1997
1998 zin = z_lake(c,j) - 0.5_kind_lake*dz_lake(c,j)
1999 zout = z_lake(c,j) + 0.5_kind_lake*dz_lake(c,j)
2000 rsfin = exp( -eta(islak)*max( zin-za(islak),0._kind_lake ) )
2001 rsfout = exp( -eta(islak)*max( zout-za(islak),0._kind_lake ) )
2002
2003 ! Let rsfout for bottom layer go into soil.
2004 ! This looks like it should be robust even for pathological cases,
2005 ! like lakes thinner than za.
2006 if (t_grnd(c) > tfrz .and. t_lake(c,1) > tfrz .and. snl(c) == 0) then
2007 phidum = (rsfin-rsfout) * sabg(p) * (1._kind_lake-beta(islak))
2008 if (j == nlevlake) then
2009 phi_soil(c) = rsfout * sabg(p) * (1._kind_lake-beta(islak))
2010 end if
2011 else if (j == 1 .and. snl(c) == 0) then !if frozen but no snow layers
2012 phidum = sabg(p) * (1._kind_lake-beta(islak))
2013 else !radiation absorbed at surface
2014 phidum = 0._kind_lake
2015 if (j == nlevlake) phi_soil(c) = 0._kind_lake
2016 end if
2017 phi(c,j) = phidum
2018
2019 end do
2020 end do
2021
2022 ! 5!) Set thermal properties and check initial energy content.
2023
2024 ! For lake
2025 do j = 1, nlevlake
2026 !dir$ concurrent
2027 !cdir nodep
2028 do fc = 1, num_shlakec
2029 c = filter_shlakec(fc)
2030
2031 cv_lake(c,j) = dz_lake(c,j) * (cwat*(1._kind_lake-lake_icefrac(c,j)) + cice_eff*lake_icefrac(c,j))
2032 end do
2033 end do
2034
2035 ! For snow / soil
2036 call soilthermprop_lake (snl,dz,zi,z,t_soisno,h2osoi_liq,h2osoi_ice, &
2037 watsat, tksatu, tkmg, tkdry, csol, tk, cv, tktopsoillay,errmsg,errflg)
2038 if(errflg/=0) then
2039 ! State is no longer valid, so return error to caller
2040 ! FIXME: PUT THIS BACK return
2041 endif
2042
2043 ! Sum cv*t_lake for energy check
2044 ! Include latent heat term, and correction for changing heat capacity with phase change.
2045
2046 ! This will need to be over all soil / lake / snow layers. Lake is below.
2047 do j = 1, nlevlake
2048 !dir$ concurrent
2049 !cdir nodep
2050 do fc = 1, num_shlakec
2051 c = filter_shlakec(fc)
2052
2053 ! ocvts(c) = ocvts(c) + cv_lake(c,j)*t_lake(c,j) &
2054 ocvts(c) = ocvts(c) + cv_lake(c,j)*(t_lake(c,j)-tfrz) &
2055 + cfus*dz_lake(c,j)*(1._kind_lake-lake_icefrac(c,j)) !&
2056 ! + (cwat-cice_eff)*lake_icefrac(c)*tfrz*dz_lake(c,j) !enthalpy reconciliation term
2057 t_lake_bef(c,j) = t_lake(c,j)
2058 if(debug_print) then
2059 if (abs(xlat_d-52.1152).lt.0.1 .and. &
2060 abs(xlon_d-260.405).lt.0.1)then
2061 print *,' ocvts(c) at xlat_d,xlon_d',xlat_d,xlon_d
2062 print *,'j,dz_lake(c,j) ', j,dz_lake(c,j)
2063 print*,'cv_lake(c,j),lake_icefrac(c,j),t_lake(c,j),cfus,ocvts(c)', &
2064 cv_lake(c,j),lake_icefrac(c,j),t_lake(c,j),cfus,ocvts(c)
2065 endif
2066 endif
2067 end do
2068 end do
2069
2070 ! Now do for soil / snow layers
2071 do j = -nlevsnow + 1, nlevsoil
2072 !dir$ concurrent
2073 !cdir nodep
2074 do fc = 1, num_shlakec
2075 c = filter_shlakec(fc)
2076
2077 if (j >= jtop(c)) then
2078 ! ocvts(c) = ocvts(c) + cv(c,j)*t_soisno(c,j) &
2079 ocvts(c) = ocvts(c) + cv(c,j)*(t_soisno(c,j)-tfrz) &
2080 + hfus*h2osoi_liq(c,j) !&
2081 ! + (cpliq-cpice)*h2osoi_ice(c,j)*tfrz !enthalpy reconciliation term
2082 if(debug_print) then
2083 if (abs(xlat_d-52.1152).lt.0.1 .and. &
2084 abs(xlon_d-260.405).lt.0.1)then
2085 print *,' ocvts(c) at xlat_d,xlon_d',xlat_d,xlon_d
2086 print *,' j,jtop(c)',j,jtop(c),'h2osoi_liq(c,j) ',h2osoi_liq(c,j),'h2osoi_ice(c,j)',h2osoi_ice(c,j)
2087 print *,' cv(c,j),t_soisno(c,j),hfus,ocvts(c)',c,j,cv(c,j),t_soisno(c,j),hfus,ocvts(c)
2088 print *,' h2osno(c)',h2osno(c)
2089 endif
2090 endif
2091 if (j == 1 .and. h2osno(c) > 0._kind_lake .and. j == jtop(c)) then
2092 ocvts(c) = ocvts(c) - h2osno(c)*hfus
2093 end if
2094 t_soisno_bef(c,j) = t_soisno(c,j)
2095 if(abs(t_soisno(c,j)-288) > 150) then
209648 format('WARNING: At c=',i0,' level=',i0,' extreme t_soisno = ',f15.10)
2097 WRITE(message,48) c,j,t_soisno(c,j)
2098 ! errmsg=trim(message)
2099 ! errflg=1
2100 write(0,'(A)') trim(message)
2101 endif
2102 end if
2103 end do
2104 end do
2105
2106 !!!!!!!!!!!!!!!!!!!
2107 ! 6!) Set up vector r and vectors a, b, c1 that define tridiagonal matrix
2108
2109 ! Heat capacity and resistance of snow without snow layers (<1cm) is ignored during diffusion,
2110 ! but its capacity to absorb latent heat may be used during phase change.
2111
2112 ! Set up interface depths, zx, heat capacities, cvx, solar source terms, phix, and temperatures, tx.
2113 do j = -nlevsnow+1, nlevlake+nlevsoil
2114 !dir$ prefervector
2115 !dir$ concurrent
2116 !cdir nodep
2117 do fc = 1,num_shlakec
2118 c = filter_shlakec(fc)
2119
2120 jprime = j - nlevlake
2121
2122 if (j >= jtop(c)) then
2123 if (j < 1) then !snow layer
2124 zx(c,j) = z(c,j)
2125 cvx(c,j) = cv(c,j)
2126 phix(c,j) = 0._kind_lake
2127 tx(c,j) = t_soisno(c,j)
2128 else if (j <= nlevlake) then !lake layer
2129 zx(c,j) = z_lake(c,j)
2130 cvx(c,j) = cv_lake(c,j)
2131 phix(c,j) = phi(c,j)
2132 tx(c,j) = t_lake(c,j)
2133 else !soil layer
2134 zx(c,j) = zx(c,nlevlake) + dz_lake(c,nlevlake)*0.5_kind_lake + z(c,jprime)
2135 cvx(c,j) = cv(c,jprime)
2136 if (j == nlevlake + 1) then !top soil layer
2137 phix(c,j) = phi_soil(c)
2138 else !middle or bottom soil layer
2139 phix(c,j) = 0._kind_lake
2140 end if
2141 tx(c,j) = t_soisno(c,jprime)
2142 end if
2143 end if
2144
2145 end do
2146 end do
2147
2148 ! Determine interface thermal conductivities, tkix
2149
2150 do j = -nlevsnow+1, nlevlake+nlevsoil
2151 !dir$ prefervector
2152 !dir$ concurrent
2153 !cdir nodep
2154 do fc = 1,num_shlakec
2155 c = filter_shlakec(fc)
2156
2157 jprime = j - nlevlake
2158
2159 if (j >= jtop(c)) then
2160 if (j < 0) then !non-bottom snow layer
2161 tkix(c,j) = tk(c,j)
2162 else if (j == 0) then !bottom snow layer
2163 dzp = zx(c,j+1) - zx(c,j)
2164 tkix(c,j) = tk_lake(c,1)*tk(c,j)*dzp / &
2165 (tk(c,j)*z_lake(c,1) + tk_lake(c,1)*(-z(c,j)) )
2166 ! tk(c,0) is the conductivity at the middle of that layer, as defined in SoilThermProp_Lake
2167 else if (j < nlevlake) then !non-bottom lake layer
2168 tkix(c,j) = ( tk_lake(c,j)*tk_lake(c,j+1) * (dz_lake(c,j+1)+dz_lake(c,j)) ) &
2169 / ( tk_lake(c,j)*dz_lake(c,j+1) + tk_lake(c,j+1)*dz_lake(c,j) )
2170 else if (j == nlevlake) then !bottom lake layer
2171 dzp = zx(c,j+1) - zx(c,j)
2172 tkix(c,j) = (tktopsoillay(c)*tk_lake(c,j)*dzp / &
2173 (tktopsoillay(c)*dz_lake(c,j)*0.5_kind_lake + tk_lake(c,j)*z(c,1) ) )
2174 ! tktopsoillay is the conductivity at the middle of that layer, as defined in SoilThermProp_Lake
2175 else !soil layer
2176 tkix(c,j) = tk(c,jprime)
2177 end if
2178 end if
2179
2180 end do
2181 end do
2182
2183
2184 ! Determine heat diffusion through the layer interface and factor used in computing
2185 ! tridiagonal matrix and set up vector r and vectors a, b, c1 that define tridiagonal
2186 ! matrix and solve system
2187
2188 do j = -nlevsnow+1, nlevlake+nlevsoil
2189 !dir$ prefervector
2190 !dir$ concurrent
2191 !cdir nodep
2192 do fc = 1,num_shlakec
2193 c = filter_shlakec(fc)
2194 if (j >= jtop(c)) then
2195 if (j < nlevlake+nlevsoil) then !top or interior layer
2196 factx(c,j) = dtime/cvx(c,j)
2197 fnx(c,j) = tkix(c,j)*(tx(c,j+1)-tx(c,j))/(zx(c,j+1)-zx(c,j))
2198 else !bottom soil layer
2199 factx(c,j) = dtime/cvx(c,j)
2200 fnx(c,j) = 0._kind_lake !not used
2201 end if
2202 end if
2203 enddo
2204 end do
2205
2206 do j = -nlevsnow+1,nlevlake+nlevsoil
2207 !dir$ prefervector
2208 !dir$ concurrent
2209 !cdir nodep
2210 do fc = 1,num_shlakec
2211 c = filter_shlakec(fc)
2212 if (j >= jtop(c)) then
2213 if (j == jtop(c)) then !top layer
2214 dzp = zx(c,j+1)-zx(c,j)
2215 a(c,j) = 0._kind_lake
2216 b(c,j) = 1+(1._kind_lake-cnfac)*factx(c,j)*tkix(c,j)/dzp
2217 c1(c,j) = -(1._kind_lake-cnfac)*factx(c,j)*tkix(c,j)/dzp
2218 r(c,j) = tx(c,j) + factx(c,j)*( fin(c) + phix(c,j) + cnfac*fnx(c,j) )
2219 else if (j < nlevlake+nlevsoil) then !middle layer
2220 dzm = (zx(c,j)-zx(c,j-1))
2221 dzp = (zx(c,j+1)-zx(c,j))
2222 a(c,j) = - (1._kind_lake-cnfac)*factx(c,j)* tkix(c,j-1)/dzm
2223 b(c,j) = 1._kind_lake+ (1._kind_lake-cnfac)*factx(c,j)*(tkix(c,j)/dzp + tkix(c,j-1)/dzm)
2224 c1(c,j) = - (1._kind_lake-cnfac)*factx(c,j)* tkix(c,j)/dzp
2225 r(c,j) = tx(c,j) + cnfac*factx(c,j)*( fnx(c,j) - fnx(c,j-1) ) + factx(c,j)*phix(c,j)
2226 else !bottom soil layer
2227 dzm = (zx(c,j)-zx(c,j-1))
2228 a(c,j) = - (1._kind_lake-cnfac)*factx(c,j)*tkix(c,j-1)/dzm
2229 b(c,j) = 1._kind_lake+ (1._kind_lake-cnfac)*factx(c,j)*tkix(c,j-1)/dzm
2230 c1(c,j) = 0._kind_lake
2231 r(c,j) = tx(c,j) - cnfac*factx(c,j)*fnx(c,j-1)
2232 end if
2233 end if
2234 enddo
2235 end do
2236 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2237
2238
2239 ! 7!) Solve for tdsolution
2240
2241 call tridiagonal(lbc, ubc, -nlevsnow + 1, nlevlake + nlevsoil, jtop, num_shlakec, filter_shlakec, &
2242 a, b, c1, r, tx)
2243
2244 ! Set t_soisno and t_lake
2245 do j = -nlevsnow+1, nlevlake + nlevsoil
2246 !dir$ concurrent
2247 !cdir nodep
2248 do fc = 1, num_shlakec
2249 c = filter_shlakec(fc)
2250
2251 jprime = j - nlevlake
2252
2253 ! Don't do anything with invalid snow layers.
2254 if (j >= jtop(c)) then
2255 if (j < 1) then !snow layer
2256 t_soisno(c,j) = tx(c,j)
2257 else if (j <= nlevlake) then !lake layer
2258 t_lake(c,j) = tx(c,j)
2259 else !soil layer
2260 t_soisno(c,jprime) = tx(c,j)
2261 end if
2262 end if
2263 end do
2264 end do
2265
2266 !!!!!!!!!!!!!!!!!!!!!!!
2267
2268 ! 8!) Sum energy content and total energy into lake for energy check. Any errors will be from the
2269 ! Tridiagonal solution.
2270
2271 if_debug_energy: if (lakedebug) then
2272 do j = 1, nlevlake
2273 !dir$ concurrent
2274 !cdir nodep
2275 do fc = 1, num_shlakec
2276 c = filter_shlakec(fc)
2277
2278 esum1(c) = esum1(c) + (t_lake(c,j)-t_lake_bef(c,j))*cv_lake(c,j)
2279 esum2(c) = esum2(c) + (t_lake(c,j)-tfrz)*cv_lake(c,j)
2280 end do
2281 end do
2282
2283 do j = -nlevsnow+1, nlevsoil
2284 !dir$ concurrent
2285 !cdir nodep
2286 do fc = 1, num_shlakec
2287 c = filter_shlakec(fc)
2288
2289 if (j >= jtop(c)) then
2290 esum1(c) = esum1(c) + (t_soisno(c,j)-t_soisno_bef(c,j))*cv(c,j)
2291 esum2(c) = esum2(c) + (t_soisno(c,j)-tfrz)*cv(c,j)
2292 end if
2293 end do
2294 end do
2295
2296 !dir$ concurrent
2297 !cdir nodep
2298 do fp = 1, num_shlakep
2299 p = filter_shlakep(fp)
2300 c = pcolumn(p)
2301 ! Again assuming only one pft per column
2302 ! esum1(c) = esum1(c) + lhabs(c)
2303 errsoi(c) = esum1(c)/dtime - eflx_soil_grnd(p)
2304 ! eflx_soil_grnd includes all the solar radiation absorbed in the lake,
2305 ! unlike eflx_gnet
2306 if(abs(errsoi(c)) > .001_kind_lake) then ! 1.e-5_kind_lake) then
2307 WRITE( message,* )'Primary soil energy conservation error in shlake &
2308 column during Tridiagonal Solution,', 'error (W/m^2):', c, errsoi(c)
2309 errmsg=trim(message)
2310 errflg=1
2311 return
2312 end if
2313 end do
2314 ! This has to be done before convective mixing because the heat capacities for each layer
2315 ! will get scrambled.
2316
2317 end if if_debug_energy
2318
2319 !!!!!!!!!!!!!!!!!!!!!!!
2320
2321 ! 9!) Phase change
2322 call phasechange_lake (snl,h2osno,dz,dz_lake, & !i
2323 t_soisno,h2osoi_liq,h2osoi_ice, & !i&o
2324 lake_icefrac,t_lake, snowdp, & !i&o
2325 qflx_snomelt,eflx_snomelt,imelt, & !o
2326 cv, cv_lake, & !i&o
2327 lhabs) !o
2328
2329 !!!!!!!!!!!!!!!!!!!!!!!
2330
2331 ! 9.5!) Second energy check and water check. Now check energy balance before and after phase
2332 ! change, considering the possibility of changed heat capacity during phase change, by
2333 ! using initial heat capacity in the first step, final heat capacity in the second step,
2334 ! and differences from tfrz only to avoid enthalpy correction for (cpliq-cpice)*melt*tfrz.
2335 ! Also check soil water sum.
2336
2337 if_debug_balance: if (lakedebug) then
2338 do j = 1, nlevlake
2339 !dir$ concurrent
2340 !cdir nodep
2341 do fc = 1, num_shlakec
2342 c = filter_shlakec(fc)
2343
2344 esum2(c) = esum2(c) - (t_lake(c,j)-tfrz)*cv_lake(c,j)
2345 end do
2346 end do
2347
2348 do j = -nlevsnow+1, nlevsoil
2349 !dir$ concurrent
2350 !cdir nodep
2351 do fc = 1, num_shlakec
2352 c = filter_shlakec(fc)
2353
2354 if (j >= jtop(c)) then
2355 esum2(c) = esum2(c) - (t_soisno(c,j)-tfrz)*cv(c,j)
2356 end if
2357 end do
2358 end do
2359
2360 !dir$ concurrent
2361 !cdir nodep
2362 do fp = 1, num_shlakep
2363 p = filter_shlakep(fp)
2364 c = pcolumn(p)
2365 ! Again assuming only one pft per column
2366 esum2(c) = esum2(c) - lhabs(c)
2367 errsoi(c) = esum2(c)/dtime
2368 if(abs(errsoi(c)) > 1.e-5_kind_lake) then
2369 write(message,*)'Primary soil energy conservation error in shlake column during Phase Change, error (W/m^2):', &
2370 c, errsoi(c)
2371 errmsg=trim(message)
2372 errflg=1
2373 return
2374 end if
2375 end do
2376
2377 ! Check soil water
2378 ! Sum soil water.
2379 do j = 1, nlevsoil
2380 !dir$ concurrent
2381 !cdir nodep
2382 do fc = 1, num_shlakec
2383 c = filter_shlakec(fc)
2384 if (j == 1) wsum_end(c) = 0._kind_lake
2385 wsum_end(c) = wsum_end(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
2386 if (j == nlevsoil) then
2387 if (abs(wsum(c)-wsum_end(c))>1.e-7_kind_lake) then
2388 write(message,*)'Soil water balance error during phase change in ShalLakeTemperature.', &
2389 'column, error (kg/m^2):', c, wsum_end(c)-wsum(c)
2390 errmsg=trim(message)
2391 errflg=1
2392 return
2393 end if
2394 end if
2395 end do
2396 end do
2397
2398 endif if_debug_balance
2399
2400 !!!!!!!!!!!!!!!!!!!!!!!!!!
2401 ! 10!) Convective mixing: make sure fracice*dz is conserved, heat content c*dz*T is conserved, and
2402 ! all ice ends up at the top. Done over all lakes even if frozen.
2403 ! Either an unstable density profile or ice in a layer below an incompletely frozen layer will trigger.
2404
2405 !Recalculate density
2406 do j = 1, nlevlake
2407 !dir$ concurrent
2408 !cdir nodep
2409 do fc = 1, num_shlakec
2410 c = filter_shlakec(fc)
2411 rhow(c,j) = (1._kind_lake - lake_icefrac(c,j)) * &
2412 1000._kind_lake*( 1.0_kind_lake - 1.9549e-05_kind_lake*(abs(t_lake(c,j)-277._kind_lake))**1.68_kind_lake ) &
2413 + lake_icefrac(c,j)*denice
2414 end do
2415 end do
2416
2417 do j = 1, nlevlake-1
2418 !dir$ concurrent
2419 !cdir nodep
2420 do fc = 1, num_shlakec
2421 c = filter_shlakec(fc)
2422 qav(c) = 0._kind_lake
2423 nav(c) = 0._kind_lake
2424 iceav(c) = 0._kind_lake
2425 end do
2426
2427 do i = 1, j+1
2428 !dir$ concurrent
2429 !cdir nodep
2430 do fc = 1, num_shlakec
2431 c = filter_shlakec(fc)
2432 if (rhow(c,j) > rhow(c,j+1) .or. &
2433 (lake_icefrac(c,j) < 1._kind_lake .and. lake_icefrac(c,j+1) > 0._kind_lake) ) then
2434 if(debug_print) then
2435 if (i==1) then
2436 print *, 'Convective Ice Mixing in column ', c, 'lake_icefrac(c,j) ',lake_icefrac(c,j),lake_icefrac(c,j+1)
2437 endif
2438 endif
2439 qav(c) = qav(c) + dz_lake(c,i)*(t_lake(c,i)-tfrz) * &
2440 ((1._kind_lake - lake_icefrac(c,i))*cwat + lake_icefrac(c,i)*cice_eff)
2441 ! tav(c) = tav(c) + t_lake(c,i)*dz_lake(c,i)
2442 iceav(c) = iceav(c) + lake_icefrac(c,i)*dz_lake(c,i)
2443 nav(c) = nav(c) + dz_lake(c,i)
2444 end if
2445 end do
2446 end do
2447
2448 !dir$ concurrent
2449 !cdir nodep
2450 do fc = 1, num_shlakec
2451 c = filter_shlakec(fc)
2452 if (rhow(c,j) > rhow(c,j+1) .or. &
2453 (lake_icefrac(c,j) < 1._kind_lake .and. lake_icefrac(c,j+1) > 0._kind_lake) ) then
2454 qav(c) = qav(c)/nav(c)
2455 iceav(c) = iceav(c)/nav(c)
2456 !If the average temperature is above freezing, put the extra energy into the water.
2457 !If it is below freezing, take it away from the ice.
2458 if (qav(c) > 0._kind_lake) then
2459 tav_froz(c) = 0._kind_lake !Celsius
2460 tav_unfr(c) = qav(c) / ((1._kind_lake - iceav(c))*cwat)
2461 else if (qav(c) < 0._kind_lake) then
2462 tav_froz(c) = qav(c) / (iceav(c)*cice_eff)
2463 tav_unfr(c) = 0._kind_lake !Celsius
2464 else
2465 tav_froz(c) = 0._kind_lake
2466 tav_unfr(c) = 0._kind_lake
2467 end if
2468 end if
2469 end do
2470
2471 do i = 1, j+1
2472 !dir$ concurrent
2473 !cdir nodep
2474 do fc = 1, num_shlakec
2475 c = filter_shlakec(fc)
2476 if (nav(c) > 0._kind_lake) then
2477 ! if(0==1) then
2478
2479 !Put all the ice at the top.!
2480 !If the average temperature is above freezing, put the extra energy into the water.
2481 !If it is below freezing, take it away from the ice.
2482 !For the layer with both ice & water, be careful to use the average temperature
2483 !that preserves the correct total heat content given what the heat capacity of that
2484 !layer will actually be.
2485 if (i == 1) zsum(c) = 0._kind_lake
2486 if ((zsum(c)+dz_lake(c,i))/nav(c) <= iceav(c)) then
2487 t_lake(c,i) = tav_froz(c) + tfrz
2488 !tgs - 30jul19 - the next line is a bug and should be commented
2489 !out. This bug prevents lake ice form completely melting.
2490 ! lake_icefrac(c,i) = 1._kind_lake
2491 else if (zsum(c)/nav(c) < iceav(c)) then
2492 !tgs - change ice fraction
2493 lake_icefrac(c,i) = (iceav(c)*nav(c) - zsum(c)) / dz_lake(c,i)
2494 ! Find average value that preserves correct heat content.
2495 t_lake(c,i) = ( lake_icefrac(c,i)*tav_froz(c)*cice_eff &
2496 + (1._kind_lake - lake_icefrac(c,i))*tav_unfr(c)*cwat ) &
2497 / ( lake_icefrac(c,i)*cice_eff + (1-lake_icefrac(c,i))*cwat ) + tfrz
2498 else
2499 !tgs - remove ice
2500 lake_icefrac(c,i) = 0._kind_lake
2501 t_lake(c,i) = tav_unfr(c) + tfrz
2502 end if
2503 zsum(c) = zsum(c) + dz_lake(c,i)
2504
2505 rhow(c,i) = (1._kind_lake - lake_icefrac(c,i)) * &
2506 1000._kind_lake*( 1.0_kind_lake - 1.9549e-05_kind_lake*(abs(t_lake(c,i)-277._kind_lake))**1.68_kind_lake ) &
2507 + lake_icefrac(c,i)*denice
2508 if (debug_print .and. lake_icefrac(c,j) > 0.)print *,'rhow(c,i),lake_icefrac(c,i),t_lake(c,i)', &
2509 i,rhow(c,i),lake_icefrac(c,i),t_lake(c,i),denice
2510 end if
2511 end do
2512 end do
2513 end do
2514
2515 !!!!!!!!!!!!!!!!!!!!!!!
2516 ! 11!) Re-evaluate thermal properties and sum energy content.
2517 ! For lake
2518 do j = 1, nlevlake
2519 !dir$ concurrent
2520 !cdir nodep
2521 do fc = 1, num_shlakec
2522 c = filter_shlakec(fc)
2523
2524 cv_lake(c,j) = dz_lake(c,j) * (cwat*(1._kind_lake-lake_icefrac(c,j)) + cice_eff*lake_icefrac(c,j))
2525 if (debug_print .and. lake_icefrac(c,j) > 0.) then
2526 print *,'Lake Ice Fraction, c, level:', c, j, lake_icefrac(c,j)
2527 endif
2528 end do
2529 end do
2530 ! For snow / soil
2531 ! call SoilThermProp_Lake(lbc, ubc, num_shlakec, filter_shlakec, tk, cv, tktopsoillay)
2532 call soilthermprop_lake (snl,dz,zi,z,t_soisno,h2osoi_liq,h2osoi_ice, &
2533 watsat, tksatu, tkmg, tkdry, csol, tk, cv, tktopsoillay,errmsg,errflg)
2534
2535
2536 ! Do as above to sum energy content
2537 do j = 1, nlevlake
2538 !dir$ concurrent
2539 !cdir nodep
2540 do fc = 1, num_shlakec
2541 c = filter_shlakec(fc)
2542
2543 ! ncvts(c) = ncvts(c) + cv_lake(c,j)*t_lake(c,j) &
2544 ncvts(c) = ncvts(c) + cv_lake(c,j)*(t_lake(c,j)-tfrz) &
2545 + cfus*dz_lake(c,j)*(1._kind_lake-lake_icefrac(c,j)) !&
2546 ! + (cwat-cice_eff)*lake_icefrac(c)*tfrz*dz_lake(c,j) !enthalpy reconciliation term
2547 fin(c) = fin(c) + phi(c,j)
2548 if(debug_print) then
2549 if (abs(xlat_d-52.1152).lt.0.1 .and. &
2550 abs(xlon_d-260.405).lt.0.1)then
2551 print *,' ncvts(c) at xlat_d,xlon_d',xlat_d,xlon_d
2552 print *,' new cv_lake(c,j),t_lake(c,j),cfus,lake_icefrac(c,j),ncvts(c),fin(c)', &
2553 j,cv_lake(c,j),t_lake(c,j),cfus,lake_icefrac(c,j),ncvts(c),fin(c)
2554 print *,' new dz_lake(c,j),fin(c),phi(c,j)',c,dz_lake(c,j),fin(c),phi(c,j)
2555 endif
2556 endif
2557 end do
2558 end do
2559
2560 do j = -nlevsnow + 1, nlevsoil
2561 !dir$ concurrent
2562 !cdir nodep
2563 do fc = 1, num_shlakec
2564 c = filter_shlakec(fc)
2565
2566 if (j >= jtop(c)) then
2567 ! ncvts(c) = ncvts(c) + cv(c,j)*t_soisno(c,j) &
2568 ncvts(c) = ncvts(c) + cv(c,j)*(t_soisno(c,j)-tfrz) &
2569 + hfus*h2osoi_liq(c,j) !&
2570 ! + (cpliq-cpice)*h2osoi_ice(c,j)*tfrz !enthalpy reconciliation term
2571 if(debug_print) then
2572 if (abs(xlat_d-52.1152).lt.0.1 .and. &
2573 abs(xlon_d-260.405).lt.0.1)then
2574 print *,' ncvts(c) at xlat_d,xlon_d',xlat_d,xlon_d
2575 print *,'new j,jtop(c)',j,jtop(c),'h2osoi_liq(c,j) ',h2osoi_liq(c,j),'h2osoi_ice(c,j)',h2osoi_ice(c,j)
2576 print *,'new cv(c,j),t_soisno(c,j),hfus,ncvts(c)',c,j,cv(c,j),t_soisno(c,j),hfus,ncvts(c)
2577 print *,'new h2osno(c)',h2osno(c)
2578 endif
2579 endif
2580 if (j == 1 .and. h2osno(c) > 0._kind_lake .and. j == jtop(c)) then
2581 ncvts(c) = ncvts(c) - h2osno(c)*hfus
2582 end if
2583 end if
2584 if (j == 1) fin(c) = fin(c) + phi_soil(c)
2585 end do
2586 end do
2587
2588
2589 ! Check energy conservation.
2590
2591 do fp = 1, num_shlakep
2592 p = filter_shlakep(fp)
2593 c = pcolumn(p)
2594 errsoi(c) = (ncvts(c)-ocvts(c)) / dtime - fin(c)
2595 if(debug_print) then
2596 if (abs(xlat_d-52.1152).lt.0.1 .and. &
2597 abs(xlon_d-260.405).lt.0.1)then
2598 print *,'xlat_d,xlon_d',xlat_d,xlon_d
2599 print *,'errsoi(c),fin(c),ncvts(c),ocvts(c),dtime,lake_icefrac(c,:),h2osno(c)', &
2600 errsoi(c),fin(c),ncvts(c),ocvts(c),dtime,lake_icefrac(c,:),h2osno(c)
2601 endif
2602 endif
2603 if( .not.lakedebug ) then
2604 if (abs(errsoi(c)) < 10._kind_lake) then
2605 eflx_sh_tot(p) = eflx_sh_tot(p) - errsoi(c)
2606 eflx_sh_grnd(p) = eflx_sh_grnd(p) - errsoi(c)
2607 eflx_soil_grnd(p) = eflx_soil_grnd(p) + errsoi(c)
2608 eflx_gnet(p) = eflx_gnet(p) + errsoi(c)
2609 if(debug_print) then
2610 if (abs(errsoi(c)) > 1.e-1_kind_lake) then
2611 print *,'errsoi incorporated at xlat_d,xlon_d',xlat_d,xlon_d
2612 print *,'errsoi incorporated into sensible heat in ShalLakeTemperature: c, (W/m^2):', c, errsoi(c)
2613 end if
2614 endif
2615 errsoi(c) = 0._kind_lake
2616 endif
2617 elseif ( lakedebug) then
2618 if (abs(errsoi(c)) < 1._kind_lake) then
2619 eflx_sh_tot(p) = eflx_sh_tot(p) - errsoi(c)
2620 eflx_sh_grnd(p) = eflx_sh_grnd(p) - errsoi(c)
2621 eflx_soil_grnd(p) = eflx_soil_grnd(p) + errsoi(c)
2622 eflx_gnet(p) = eflx_gnet(p) + errsoi(c)
2623 if (abs(errsoi(c)) > 1.e-1_kind_lake) then
2624 print *,'errsoi incorporated into sensible heat in ShalLakeTemperature: c, (W/m^2):', c, errsoi(c)
2625 end if
2626 errsoi(c) = 0._kind_lake
2627 else
2628 print *,'Soil Energy Balance Error at column, ', c, 'G, fintotal, column E tendency = ', &
2629 eflx_gnet(p), fin(c), (ncvts(c)-ocvts(c)) / dtime,'xlat_d,xlon_d',xlat_d,xlon_d
2630 print *,'errsoi(c),ncvts(c),ocvts(c)',errsoi(c),ncvts(c),ocvts(c),'xlat_d,xlon_d',xlat_d,xlon_d
2631 print *,'lake_icefrac(c,:),h2osno(c)', lake_icefrac(c,:),h2osno(c)
2632 print *,'t_lake(c,:),t_soisno(c,:)',t_lake(c,:),t_soisno(c,:)
2633 end if
2634 end if ! LAKEDEBUG
2635 end do
2636 ! This loop assumes only one point per column.
2637
2638 end subroutine shallaketemperature
2639
2640 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2641 !-----------------------------------------------------------------------
2642 !BOP
2643 !
2644 ! ROUTINE: SoilThermProp_Lake
2645 !
2646 ! INTERFACE:
2647 ! DESCRIPTION:
2662 subroutine soilthermprop_lake (snl,dz,zi,z,t_soisno,h2osoi_liq,h2osoi_ice, &
2663 watsat, tksatu, tkmg, tkdry, csol, tk, cv, tktopsoillay,errmsg,errflg)
2664
2665 implicit none
2666 !in
2667
2668 integer, intent(inout) :: errflg
2669 character(*), intent(inout) :: errmsg
2670 integer , intent(in) :: snl(1)
2671 ! real(kind_lake), intent(in) :: h2osno(1) !< snow water (mm H2O)
2672 real(kind_lake), intent(in) :: watsat(1,nlevsoil)
2673 real(kind_lake), intent(in) :: tksatu(1,nlevsoil)
2674 real(kind_lake), intent(in) :: tkmg(1,nlevsoil)
2675 real(kind_lake), intent(in) :: tkdry(1,nlevsoil)
2676 real(kind_lake), intent(in) :: csol(1,nlevsoil)
2677 real(kind_lake), intent(in) :: dz(1,-nlevsnow+1:nlevsoil)
2678 real(kind_lake), intent(in) :: zi(1,-nlevsnow+0:nlevsoil)
2679 real(kind_lake), intent(in) :: z(1,-nlevsnow+1:nlevsoil)
2680 real(kind_lake), intent(in) :: t_soisno(1,-nlevsnow+1:nlevsoil)
2681 real(kind_lake), intent(in) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil)
2682 real(kind_lake), intent(in) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil)
2683
2684 !out
2685 real(kind_lake), intent(out) :: cv(lbc:ubc,-nlevsnow+1:nlevsoil)
2686 real(kind_lake), intent(out) :: tk(lbc:ubc,-nlevsnow+1:nlevsoil)
2687 real(kind_lake), intent(out) :: tktopsoillay(lbc:ubc)!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2689 ! !CALLED FROM:
2690 ! subroutine ShalLakeTemperature in this module.
2691 !
2692 ! !REVISION HISTORY:
2693 ! 15 September 1999: Yongjiu Dai; Initial code
2694 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
2695 ! 2/13/02, Peter Thornton: migrated to new data structures
2696 ! 7/01/03, Mariana Vertenstein: migrated to vector code
2697 ! 4/09, Zack Subin, adjustment for ShalLake code.
2698 ! June 2022, Sam Trahan updated for CCPP
2699 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2700 ! !LOCAL VARIABLES:
2701 !
2702 ! local pointers to original implicit in scalars
2703 !
2704 ! integer , pointer :: clandunit(:) ! column's landunit
2705 ! integer , pointer :: ityplun(:) ! landunit type
2706 !
2707 !EOP
2708
2709
2710 ! OTHER LOCAL VARIABLES:
2711
2712 integer :: l,c,j ! indices
2713 integer :: fc ! lake filtered column indices
2714 real(kind_lake) :: bw ! partial density of water (ice + liquid)
2715 real(kind_lake) :: dksat ! thermal conductivity for saturated soil (j/(k s m))
2716 real(kind_lake) :: dke ! kersten number
2717 real(kind_lake) :: fl ! fraction of liquid or unfrozen water to total water
2718 real(kind_lake) :: satw ! relative total water content of soil.
2719 real(kind_lake) :: thk(lbc:ubc,-nlevsnow+1:nlevsoil) ! thermal conductivity of layer
2720 character*256 :: message
2721
2722 real(kind_lake) :: denom
2723
2724 ! Thermal conductivity of soil from Farouki (1981)
2725
2726 do j = -nlevsnow+1,nlevsoil
2727 !dir$ concurrent
2728 !cdir nodep
2729 do fc = 1, num_shlakec
2730 c = filter_shlakec(fc)
2731
2732 ! Only examine levels from 1->nlevsoil
2733 if (j >= 1) then
2734 ! l = clandunit(c)
2735 ! if (ityplun(l) /= istwet .AND. ityplun(l) /= istice) then
2736 ! This could be altered later for allowing this to be over glaciers.
2737
2738 ! Soil should be saturated.
2739 if (lakedebug) then
2740 satw = (h2osoi_liq(c,j)/denh2o + h2osoi_ice(c,j)/denice)/(dz(c,j)*watsat(c,j))
2741 ! satw = min(1._kind_lake, satw)
2742 if (satw < 0.999_kind_lake) then
2743 write(message,*)'WARNING: soil layer unsaturated in SoilThermProp_Lake, satw, j = ', satw, j
2744 ! errmsg=trim(message)
2745 ! errflg=1
2746 write(0,'(A)') trim(message)
2747 end if
2748 ! Could use denice because if it starts out frozen, the volume of water will go below sat.,
2749 ! since we're not yet doing excess ice.
2750 ! But take care of this in HydrologyLake.
2751 endif
2752 satw = 1._kind_lake
2753 denom = (h2osoi_ice(c,j)+h2osoi_liq(c,j))
2754 if(denom>zero_h2o) then
2755 fl = h2osoi_liq(c,j)/denom
2756 else
2757 write(message,'(A,I0)') 'WARNING: zero h2osoi_ice+h2osoi_liq at j = ', j
2758 ! errmsg=trim(message)
2759 ! errflg=1
2760 fl = 0
2761 write(0,'(A)') trim(message)
2762 endif
2763 if (t_soisno(c,j) >= tfrz) then ! Unfrozen soil
2764 dke = max(0._kind_lake, log10(satw) + 1.0_kind_lake)
2765 dksat = tksatu(c,j)
2766 else ! Frozen soil
2767 dke = satw
2768 dksat = tkmg(c,j)*0.249_kind_lake**(fl*watsat(c,j))*2.29_kind_lake**watsat(c,j)
2769 endif
2770 thk(c,j) = dke*dksat + (1._kind_lake-dke)*tkdry(c,j)
2771 ! else
2772 ! thk(c,j) = tkwat
2773 ! if (t_soisno(c,j) < tfrz) thk(c,j) = tkice
2774 ! endif
2775 endif
2776
2777 ! Thermal conductivity of snow, which from Jordan (1991) pp. 18
2778 ! Only examine levels from snl(c)+1 -> 0 where snl(c) < 1
2779 if (snl(c)+1 < 1 .AND. (j >= snl(c)+1) .AND. (j <= 0)) then
2780 bw = (h2osoi_ice(c,j)+h2osoi_liq(c,j))/dz(c,j)
2781 thk(c,j) = tkairc + (7.75e-5_kind_lake *bw + 1.105e-6_kind_lake*bw*bw)*(tkice-tkairc)
2782 end if
2783
2784 end do
2785 end do
2786
2787 ! Thermal conductivity at the layer interface
2788
2789 ! Have to correct for the fact that bottom snow layer and top soil layer border lake.
2790 ! For the first case, the snow layer conductivity for the middle of the layer will be returned.
2791 ! Because the interfaces are below the soil layers, the conductivity for the top soil layer
2792 ! will have to be returned separately.
2793 do j = -nlevsnow+1,nlevsoil
2794 !dir$ concurrent
2795 !cdir nodep
2796 do fc = 1,num_shlakec
2797 c = filter_shlakec(fc)
2798 if (j >= snl(c)+1 .AND. j <= nlevsoil-1 .AND. j /= 0) then
2799 tk(c,j) = thk(c,j)*thk(c,j+1)*(z(c,j+1)-z(c,j)) &
2800 /(thk(c,j)*(z(c,j+1)-zi(c,j))+thk(c,j+1)*(zi(c,j)-z(c,j)))
2801 else if (j == 0) then
2802 tk(c,j) = thk(c,j)
2803 else if (j == nlevsoil) then
2804 tk(c,j) = 0._kind_lake
2805 end if
2806 ! For top soil layer.
2807 if (j == 1) tktopsoillay(c) = thk(c,j)
2808 end do
2809 end do
2810
2811 ! Soil heat capacity, from de Vires (1963)
2812
2813 do j = 1, nlevsoil
2814 !dir$ concurrent
2815 !cdir nodep
2816 do fc = 1,num_shlakec
2817 c = filter_shlakec(fc)
2818 ! l = clandunit(c)
2819 ! if (ityplun(l) /= istwet .AND. ityplun(l) /= istice) then
2820 cv(c,j) = csol(c,j)*(1-watsat(c,j))*dz(c,j) + &
2821 (h2osoi_ice(c,j)*cpice + h2osoi_liq(c,j)*cpliq)
2822 ! else
2823 ! cv(c,j) = (h2osoi_ice(c,j)*cpice + h2osoi_liq(c,j)*cpliq)
2824 ! endif
2825 ! if (j == 1) then
2826 ! if (snl(c)+1 == 1 .AND. h2osno(c) > 0._kind_lake) then
2827 ! cv(c,j) = cv(c,j) + cpice*h2osno(c)
2828 ! end if
2829 ! end if
2830 ! Won't worry about heat capacity for thin snow on lake with no snow layers.
2831 enddo
2832 end do
2833
2834 ! Snow heat capacity
2835
2836 do j = -nlevsnow+1,0
2837 !dir$ concurrent
2838 !cdir nodep
2839 do fc = 1,num_shlakec
2840 c = filter_shlakec(fc)
2841 if (snl(c)+1 < 1 .and. j >= snl(c)+1) then
2842 cv(c,j) = cpliq*h2osoi_liq(c,j) + cpice*h2osoi_ice(c,j)
2843 end if
2844 end do
2845 end do
2846
2847 end subroutine soilthermprop_lake
2848
2849
2850 !-----------------------------------------------------------------------
2851 !BOP
2852 !
2853 ! ROUTINE: PhaseChange_Lake
2854 !
2855 ! !INTERFACE:
2856
2857 ! DESCRIPTION:
2874 subroutine phasechange_lake (snl,h2osno,dz,dz_lake, & !i
2875 t_soisno,h2osoi_liq,h2osoi_ice, & !i&o
2876 lake_icefrac,t_lake, snowdp, & !i&o
2877 qflx_snomelt,eflx_snomelt,imelt, & !o
2878 cv, cv_lake, & !i&o
2879 lhabs) !o
2880 !=============================================================================================
2881 !CALLED FROM:
2882 ! subroutine ShalLakeTemperature in this module
2883 !
2884 !REVISION HISTORY:
2885 ! - 04/2009 Zack Subin: Initial code
2886 ! - June 2022 Sam Trahan: Modified for CCPP
2887 !==============================================================================================
2888 ! !USES:
2889 !
2890 ! !ARGUMENTS:
2891 implicit none
2892 !in:
2893
2894 integer , intent(in) :: snl(1)
2895 real(kind_lake), intent(inout) :: h2osno(1)
2896 real(kind_lake), intent(in) :: dz(1,-nlevsnow+1:nlevsoil)
2897 real(kind_lake), intent(in) :: dz_lake(1,nlevlake)
2898 ! Needed in case snow height is less than critical value.
2899
2900 !inout:
2901
2902 real(kind_lake), intent(inout) :: snowdp(1)
2903 real(kind_lake), intent(inout) :: t_soisno(1,-nlevsnow+1:nlevsoil)
2904 real(kind_lake), intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil)
2905 real(kind_lake), intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil)
2906 real(kind_lake), intent(inout) :: lake_icefrac(1,nlevlake)
2907 real(kind_lake), intent(inout) :: t_lake(1,nlevlake)
2908 !out:
2909
2910 real(kind_lake), intent(out) :: qflx_snomelt(1)
2911 real(kind_lake), intent(out) :: eflx_snomelt(1)
2912 integer, intent(out) :: imelt(1,-nlevsnow+1:nlevsoil)
2913 !What's the sign of this? Is it just output?
2914 real(kind_lake), intent(inout) :: cv(lbc:ubc,-nlevsnow+1:nlevsoil)
2915 real(kind_lake), intent(inout) :: cv_lake (lbc:ubc,1:nlevlake)
2916 real(kind_lake), intent(out):: lhabs(lbc:ubc)
2917
2918
2919 ! OTHER LOCAL VARIABLES:
2920
2921 integer :: j,c,g !do loop index
2922 integer :: fc !lake filtered column indices
2923 real(kind_lake) :: heatavail !available energy for melting or freezing (J/m^2)
2924 real(kind_lake) :: heatrem !energy residual or loss after melting or freezing
2925 real(kind_lake) :: melt !actual melting (+) or freezing (-) [kg/m2]
2926 real(kind_lake), parameter :: smallnumber = 1.e-7_kind_lake !to prevent tiny residuals from rounding error
2927 logical :: dophasechangeflag
2928 !-----------------------------------------------------------------------
2929
2930 ! Initialization
2931
2932 !dir$ concurrent
2933 !cdir nodep
2934 do fc = 1,num_shlakec
2935 c = filter_shlakec(fc)
2936
2937 qflx_snomelt(c) = 0._kind_lake
2938 eflx_snomelt(c) = 0._kind_lake
2939 lhabs(c) = 0._kind_lake
2940 end do
2941
2942 do j = -nlevsnow+1,0
2943 !dir$ concurrent
2944 !cdir nodep
2945 do fc = 1,num_shlakec
2946 c = filter_shlakec(fc)
2947
2948 if (j >= snl(c) + 1) imelt(c,j) = 0
2949 end do
2950 end do
2951
2952 ! Check for case of snow without snow layers and top lake layer temp above freezing.
2953
2954 !dir$ concurrent
2955 !cdir nodep
2956 do fc = 1,num_shlakec
2957 c = filter_shlakec(fc)
2958
2959 if (snl(c) == 0 .and. h2osno(c) > 0._kind_lake .and. t_lake(c,1) > tfrz) then
2960 heatavail = (t_lake(c,1) - tfrz) * cv_lake(c,1)
2961 melt = min(h2osno(c), heatavail/hfus)
2962 heatrem = max(heatavail - melt*hfus, 0._kind_lake)
2963 !catch small negative value to keep t at tfrz
2964 t_lake(c,1) = tfrz + heatrem/(cv_lake(c,1))
2965 snowdp(c) = snowdp(c)*(1._kind_lake - melt/h2osno(c))
2966 h2osno(c) = h2osno(c) - melt
2967 lhabs(c) = lhabs(c) + melt*hfus
2968 qflx_snomelt(c) = qflx_snomelt(c) + melt
2969 ! Prevent tiny residuals
2970 if (h2osno(c) < smallnumber) h2osno(c) = 0._kind_lake
2971 if (snowdp(c) < smallnumber) snowdp(c) = 0._kind_lake
2972 end if
2973 end do
2974
2975 ! Lake phase change
2976
2977 do j = 1,nlevlake
2978 !dir$ concurrent
2979 !cdir nodep
2980 do fc = 1,num_shlakec
2981 c = filter_shlakec(fc)
2982
2983 dophasechangeflag = .false.
2984 if (t_lake(c,j) > tfrz .and. lake_icefrac(c,j) > 0._kind_lake) then ! melting
2985 dophasechangeflag = .true.
2986 heatavail = (t_lake(c,j) - tfrz) * cv_lake(c,j)
2987 melt = min(lake_icefrac(c,j)*denh2o*dz_lake(c,j), heatavail/hfus)
2988 !denh2o is used because layer thickness is not adjusted for freezing
2989 heatrem = max(heatavail - melt*hfus, 0._kind_lake)
2990 !catch small negative value to keep t at tfrz
2991 else if (t_lake(c,j) < tfrz .and. lake_icefrac(c,j) < 1._kind_lake) then !freezing
2992 dophasechangeflag = .true.
2993 heatavail = (t_lake(c,j) - tfrz) * cv_lake(c,j)
2994 melt = max(-(1._kind_lake-lake_icefrac(c,j))*denh2o*dz_lake(c,j), heatavail/hfus)
2995 !denh2o is used because layer thickness is not adjusted for freezing
2996 heatrem = min(heatavail - melt*hfus, 0._kind_lake)
2997 !catch small positive value to keep t at tfrz
2998 end if
2999 ! Update temperature and ice fraction.
3000 if (dophasechangeflag) then
3001 lake_icefrac(c,j) = lake_icefrac(c,j) - melt/(denh2o*dz_lake(c,j))
3002 lhabs(c) = lhabs(c) + melt*hfus
3003 ! Update heat capacity
3004 cv_lake(c,j) = cv_lake(c,j) + melt*(cpliq-cpice)
3005 t_lake(c,j) = tfrz + heatrem/cv_lake(c,j)
3006 ! Prevent tiny residuals
3007 if (lake_icefrac(c,j) > 1._kind_lake - smallnumber) lake_icefrac(c,j) = 1._kind_lake
3008 if (lake_icefrac(c,j) < smallnumber) lake_icefrac(c,j) = 0._kind_lake
3009 end if
3010 end do
3011 end do
3012
3013 ! Snow & soil phase change
3014
3015 do j = -nlevsnow+1,nlevsoil
3016 !dir$ concurrent
3017 !cdir nodep
3018 do fc = 1,num_shlakec
3019 c = filter_shlakec(fc)
3020 dophasechangeflag = .false.
3021
3022 if (j >= snl(c) + 1) then
3023
3024 if (t_soisno(c,j) > tfrz .and. h2osoi_ice(c,j) > 0._kind_lake) then ! melting
3025 dophasechangeflag = .true.
3026 heatavail = (t_soisno(c,j) - tfrz) * cv(c,j)
3027 melt = min(h2osoi_ice(c,j), heatavail/hfus)
3028 heatrem = max(heatavail - melt*hfus, 0._kind_lake)
3029 !catch small negative value to keep t at tfrz
3030 if (j <= 0) then !snow
3031 imelt(c,j) = 1
3032 qflx_snomelt(c) = qflx_snomelt(c) + melt
3033 end if
3034 else if (t_soisno(c,j) < tfrz .and. h2osoi_liq(c,j) > 0._kind_lake) then !freezing
3035 dophasechangeflag = .true.
3036 heatavail = (t_soisno(c,j) - tfrz) * cv(c,j)
3037 melt = max(-h2osoi_liq(c,j), heatavail/hfus)
3038 heatrem = min(heatavail - melt*hfus, 0._kind_lake)
3039 !catch small positive value to keep t at tfrz
3040 if (j <= 0) then !snow
3041 imelt(c,j) = 2
3042 qflx_snomelt(c) = qflx_snomelt(c) + melt
3043 ! Does this works for both signs of melt in SnowHydrology? I think
3044 ! qflx_snomelt(c) is just output.
3045 end if
3046 end if
3047
3048 ! Update temperature and soil components.
3049 if (dophasechangeflag) then
3050 h2osoi_ice(c,j) = h2osoi_ice(c,j) - melt
3051 h2osoi_liq(c,j) = h2osoi_liq(c,j) + melt
3052 lhabs(c) = lhabs(c) + melt*hfus
3053 ! Update heat capacity
3054 cv(c,j) = cv(c,j) + melt*(cpliq-cpice)
3055 t_soisno(c,j) = tfrz + heatrem/cv(c,j)
3056 ! Prevent tiny residuals
3057 if (h2osoi_ice(c,j) < smallnumber) h2osoi_ice(c,j) = 0._kind_lake
3058 if (h2osoi_liq(c,j) < smallnumber) h2osoi_liq(c,j) = 0._kind_lake
3059 end if
3060
3061 end if
3062 end do
3063 end do
3064
3065 ! Update eflx_snomelt(c)
3066 !dir$ concurrent
3067 !cdir nodep
3068 do fc = 1,num_shlakec
3069 c = filter_shlakec(fc)
3070 eflx_snomelt(c) = qflx_snomelt(c)*hfus
3071 end do
3072 !!!
3073
3074 end subroutine phasechange_lake
3075
3076 ! DESCRIPTION:
3089 subroutine shallakehydrology(dz_lake,forc_rain,forc_snow, & !i
3090 begwb,qflx_evap_tot,forc_t,do_capsnow, &
3091 t_grnd,qflx_evap_soi, &
3092 qflx_snomelt,imelt,frac_iceold, & !i add by guhp
3093 z,dz,zi,snl,h2osno,snowdp,lake_icefrac,t_lake, & !i&o
3094 endwb,snowage,snowice,snowliq,t_snow, & !o
3095 t_soisno,h2osoi_ice,h2osoi_liq,h2osoi_vol, &
3096 qflx_drain,qflx_surf,qflx_infl,qflx_qrgwl, &
3097 qcharge,qflx_prec_grnd,qflx_snowcap, &
3098 qflx_snowcap_col,qflx_snow_grnd_pft, &
3099 qflx_snow_grnd_col,qflx_rain_grnd, &
3100 qflx_evap_tot_col,soilalpha,zwt,fcov, &
3101 rootr_column,qflx_evap_grnd,qflx_sub_snow, &
3102 qflx_dew_snow,qflx_dew_grnd,qflx_rain_grnd_col, &
3103 watsat, tksatu, tkmg, tkdry, csol, &
3104 dtime,errmsg,errflg)
3105
3106 !==================================================================================
3107 !
3108 ! WARNING: This subroutine assumes lake columns have one and only one pft.
3109 !
3110 ! Sequence is:
3111 ! ShalLakeHydrology:
3112 ! Do needed tasks from Hydrology1, Biogeophysics2, & top of Hydrology2.
3113 ! -> SnowWater: change of snow mass and snow water onto soil
3114 ! -> SnowCompaction: compaction of snow layers
3115 ! -> CombineSnowLayers: combine snow layers that are thinner than minimum
3116 ! -> DivideSnowLayers: subdivide snow layers that are thicker than maximum
3117 ! Add water to soil if melting has left it with open pore space.
3118 ! Cleanup and do water balance.
3119 ! If snow layers are found above a lake with unfrozen top layer, whose top
3120 ! layer has enough heat to melt all the snow ice without freezing, do so
3121 ! and eliminate the snow layers.
3122 !
3123 ! !REVISION HISTORY:
3124 ! Created by Zack Subin, 2009
3125 !
3126 !============================================================================================
3127
3128 ! USES:
3129 !
3130 implicit none
3131
3132 ! in:
3133
3134 integer, intent(inout) :: errflg
3135 character(*), intent(inout) :: errmsg
3136
3137 real(kind_lake) :: watsat(1,nlevsoil)
3138 real(kind_lake) :: tksatu(1,nlevsoil)
3139 real(kind_lake) :: tkmg(1,nlevsoil)
3140 real(kind_lake) :: tkdry(1,nlevsoil)
3141 real(kind_lake) :: csol(1,nlevsoil)
3142
3143 ! integer , intent(in) :: clandunit(1) !< column's landunit
3144 ! integer , intent(in) :: ityplun(1) !< landunit type
3145 real(kind_lake), intent(in) :: dtime
3146 real(kind_lake), intent(in) :: dz_lake(1,nlevlake)
3147 real(kind_lake), intent(in) :: forc_rain(1)
3148 real(kind_lake), intent(in) :: forc_snow(1)
3149 real(kind_lake), intent(in) :: qflx_evap_tot(1)
3150 real(kind_lake), intent(in) :: forc_t(1)
3151
3152 !real(kind_lake), intent(in),optional :: flfall(1) !< fraction of liquid water within falling precipitation (unused)
3153
3154 logical , intent(in) :: do_capsnow(1)
3155 real(kind_lake), intent(in) :: t_grnd(1)
3156 real(kind_lake), intent(in) :: qflx_evap_soi(1)
3157 real(kind_lake), intent(in) :: qflx_snomelt(1)
3158 integer, intent(in) :: imelt(1,-nlevsnow+1:nlevsoil)
3159
3160 !inout:
3161
3162 real(kind_lake), intent(inout) :: begwb(1)
3163
3164 ! inout:
3165
3166
3167 real(kind_lake), intent(inout) :: z(1,-nlevsnow+1:nlevsoil)
3168 real(kind_lake), intent(inout) :: dz(1,-nlevsnow+1:nlevsoil)
3169 real(kind_lake), intent(inout) :: zi(1,-nlevsnow+0:nlevsoil)
3170 integer , intent(inout) :: snl(1)
3171 real(kind_lake), intent(inout) :: h2osno(1)
3172 real(kind_lake), intent(inout) :: snowdp(1)
3173 real(kind_lake), intent(inout) :: lake_icefrac(1,nlevlake)
3174 real(kind_lake), intent(inout) :: t_lake(1,nlevlake)
3175
3176 real(kind_lake), intent(inout) :: frac_iceold(1,-nlevsnow+1:nlevsoil)
3177 ! out:
3178
3179
3180 real(kind_lake), intent(out) :: endwb(1)
3181 real(kind_lake), intent(out) :: snowage(1)
3182 real(kind_lake), intent(out) :: snowice(1)
3183 real(kind_lake), intent(out) :: snowliq(1)
3184 real(kind_lake), intent(out) :: t_snow(1)
3185 real(kind_lake), intent(out) :: t_soisno(1,-nlevsnow+1:nlevsoil)
3186 real(kind_lake), intent(out) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil)
3187 real(kind_lake), intent(out) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil)
3188 real(kind_lake), intent(out) :: h2osoi_vol(1,-nlevsnow+1:nlevsoil)
3189 real(kind_lake), intent(out) :: qflx_drain(1)
3190 real(kind_lake), intent(out) :: qflx_surf(1)
3191 real(kind_lake), intent(out) :: qflx_infl(1)
3192 real(kind_lake), intent(out) :: qflx_qrgwl(1)
3193 real(kind_lake), intent(out) :: qcharge(1)
3194 real(kind_lake), intent(out) :: qflx_prec_grnd(1)
3195 real(kind_lake), intent(out) :: qflx_snowcap(1)
3196 real(kind_lake), intent(out) :: qflx_snowcap_col(1)
3197 real(kind_lake), intent(out) :: qflx_snow_grnd_pft(1)
3198 real(kind_lake), intent(out) :: qflx_snow_grnd_col(1)
3199 real(kind_lake), intent(out) :: qflx_rain_grnd(1)
3200 real(kind_lake), intent(out) :: qflx_evap_tot_col(1)
3201 real(kind_lake) ,intent(out) :: soilalpha(1)
3202 real(kind_lake), intent(out) :: zwt(1)
3203 real(kind_lake), intent(out) :: fcov(1)
3204 real(kind_lake), intent(out) :: rootr_column(1,1:nlevsoil)
3205 real(kind_lake), intent(out) :: qflx_evap_grnd(1)
3206 real(kind_lake), intent(out) :: qflx_sub_snow(1)
3207 real(kind_lake), intent(out) :: qflx_dew_snow(1)
3208 real(kind_lake), intent(out) :: qflx_dew_grnd(1)
3209 real(kind_lake), intent(out) :: qflx_rain_grnd_col(1)
3210
3211 ! Block of biogeochem currently not used.
3212 real(kind_lake), pointer :: sucsat(:,:) ! minimum soil suction (mm)
3213 real(kind_lake), pointer :: bsw(:,:) ! Clapp and Hornberger "b"
3214 real(kind_lake), pointer :: bsw2(:,:) ! Clapp and Hornberger "b" for CN code
3215 real(kind_lake), pointer :: psisat(:,:) ! soil water potential at saturation for CN code (MPa)
3216 real(kind_lake), pointer :: vwcsat(:,:) ! volumetric water content at saturation for CN code (m3/m3)
3217 real(kind_lake), pointer :: wf(:) ! soil water as frac. of whc for top 0.5 m
3218 real(kind_lake), pointer :: soilpsi(:,:) ! soil water potential in each soil layer (MPa)
3219
3220 ! OTHER LOCAL VARIABLES:
3221
3222 integer :: p,fp,g,l,c,j,fc,jtop ! indices
3223 integer :: num_shlakesnowc ! number of column snow points
3224 integer :: filter_shlakesnowc(ubc-lbc+1) ! column filter for snow points
3225 integer :: num_shlakenosnowc ! number of column non-snow points
3226 integer :: filter_shlakenosnowc(ubc-lbc+1) ! column filter for non-snow points
3227 integer :: newnode ! flag when new snow node is set, (1=yes, 0=no)
3228 real(kind_lake) :: dz_snowf ! layer thickness rate change due to precipitation [mm/s]
3229 real(kind_lake) :: bifall ! bulk density of newly fallen dry snow [kg/m3]
3230 real(kind_lake) :: fracsnow(lbp:ubp) ! frac of precipitation that is snow
3231 real(kind_lake) :: fracrain(lbp:ubp) ! frac of precipitation that is rain
3232 real(kind_lake) :: qflx_prec_grnd_snow(lbp:ubp) ! snow precipitation incident on ground [mm/s]
3233 real(kind_lake) :: qflx_prec_grnd_rain(lbp:ubp) ! rain precipitation incident on ground [mm/s]
3234 real(kind_lake) :: qflx_evap_soi_lim ! temporary evap_soi limited by top snow layer content [mm/s]
3235 real(kind_lake) :: h2osno_temp ! temporary h2osno [kg/m^2]
3236 real(kind_lake) :: sumsnowice(lbc:ubc) ! sum of snow ice if snow layers found above unfrozen lake [kg/m&2]
3237 logical :: unfrozen(lbc:ubc) ! true if top lake layer is unfrozen with snow layers above
3238 real(kind_lake) :: heatrem ! used in case above [J/m^2]
3239 real(kind_lake) :: heatsum(lbc:ubc) ! used in case above [J/m^2]
3240 real(kind_lake) :: qflx_top_soil(1) !net water input into soil from top (mm/s)
3241 character*256 :: message
3242
3243 real(kind_lake),allocatable :: snow_water(:) ! temporary sum of snow water for Bal Check [kg/m^2]
3244 !-----------------------------------------------------------------------
3245
3246 ! Determine step size
3247
3248 ! Add soil water to water balance.
3249 do j = 1, nlevsoil
3250 !dir$ concurrent
3251 !cdir nodep
3252 do fc = 1, num_shlakec
3253 c = filter_shlakec(fc)
3254 begwb(c) = begwb(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
3255 end do
3256 end do
3257
3258 !!!!!!!!!!!!!!!!!!!!!!!!!!!
3259
3260 ! Do precipitation onto ground, etc., from Hydrology1.
3261
3262 !dir$ concurrent
3263 !cdir nodep
3264 do fp = 1, num_shlakep
3265 p = filter_shlakep(fp)
3266 g = pgridcell(p)
3267 ! l = plandunit(p)
3268 c = pcolumn(p)
3269
3270 ! Precipitation onto ground (kg/(m2 s))
3271 ! ! PET, 1/18/2005: Added new terms for mass balance correction
3272 ! ! due to dynamic pft weight shifting (column-level h2ocan_loss)
3273 ! ! Because the fractionation between rain and snow is indeterminate if
3274 ! ! rain + snow = 0, I am adding this very small flux only to the rain
3275 ! ! components.
3276 ! Not relevant unless PFTs are added to lake later.
3277 ! if (frac_veg_nosno(p) == 0) then
3278 qflx_prec_grnd_snow(p) = forc_snow(g)
3279 qflx_prec_grnd_rain(p) = forc_rain(g) !+ h2ocan_loss(c)
3280 ! else
3281 ! qflx_prec_grnd_snow(p) = qflx_through_snow(p) + (qflx_candrip(p) * fracsnow(p))
3282 ! qflx_prec_grnd_rain(p) = qflx_through_rain(p) + (qflx_candrip(p) * fracrain(p)) + h2ocan_loss(c)
3283 ! end if
3284 qflx_prec_grnd(p) = qflx_prec_grnd_snow(p) + qflx_prec_grnd_rain(p)
3285
3286 if (do_capsnow(c)) then
3287 qflx_snowcap(p) = qflx_prec_grnd_snow(p) + qflx_prec_grnd_rain(p)
3288 qflx_snow_grnd_pft(p) = 0._kind_lake
3289 qflx_rain_grnd(p) = 0._kind_lake
3290 else
3291 qflx_snowcap(p) = 0._kind_lake
3292 qflx_snow_grnd_pft(p) = qflx_prec_grnd_snow(p) ! ice onto ground (mm/s)
3293 qflx_rain_grnd(p) = qflx_prec_grnd_rain(p) ! liquid water onto ground (mm/s)
3294 end if
3295 ! Assuming one PFT; needed for below
3296 qflx_snow_grnd_col(c) = qflx_snow_grnd_pft(p)
3297 qflx_rain_grnd_col(c) = qflx_rain_grnd(p)
3298
3299 end do ! (end pft loop)
3300
3301 ! Determine snow height and snow water
3302
3303 !dir$ concurrent
3304 !cdir nodep
3305 do fc = 1, num_shlakec
3306 c = filter_shlakec(fc)
3307 ! l = clandunit(c)
3308 g = cgridcell(c)
3309
3310 ! Use Alta relationship, Anderson(1976); LaChapelle(1961),
3311 ! U.S.Department of Agriculture Forest Service, Project F,
3312 ! Progress Rep. 1, Alta Avalanche Study Center:Snow Layer Densification.
3313
3314 if (do_capsnow(c)) then
3315 dz_snowf = 0._kind_lake
3316 else
3317 if (forc_t(g) > tfrz + 2._kind_lake) then
3318 bifall=50._kind_lake + 1.7_kind_lake*(17.0_kind_lake)**1.5_kind_lake
3319 else if (forc_t(g) > tfrz - 15._kind_lake) then
3320 bifall=50._kind_lake + 1.7_kind_lake*(forc_t(g) - tfrz + 15._kind_lake)**1.5_kind_lake
3321 else
3322 bifall=50._kind_lake
3323 end if
3324 dz_snowf = qflx_snow_grnd_col(c)/bifall
3325 snowdp(c) = snowdp(c) + dz_snowf*dtime
3326 h2osno(c) = h2osno(c) + qflx_snow_grnd_col(c)*dtime ! snow water equivalent (mm)
3327 end if
3328
3329 ! if (itype(l)==istwet .and. t_grnd(c)>tfrz) then
3330 ! h2osno(c)=0._kind_lake
3331 ! snowdp(c)=0._kind_lake
3332 ! snowage(c)=0._kind_lake
3333 ! end if
3334 ! Take care of this later in function.
3335
3336 ! When the snow accumulation exceeds 10 mm, initialize snow layer
3337 ! Currently, the water temperature for the precipitation is simply set
3338 ! as the surface air temperature
3339
3340 newnode = 0 ! flag for when snow node will be initialized
3341 if (snl(c) == 0 .and. qflx_snow_grnd_col(c) > 0.0_kind_lake .and. snowdp(c) >= 0.01_kind_lake) then
3342 newnode = 1
3343 snl(c) = -1
3344 dz(c,0) = snowdp(c) ! meter
3345 z(c,0) = -0.5_kind_lake*dz(c,0)
3346 zi(c,-1) = -dz(c,0)
3347 snowage(c) = 0._kind_lake ! snow age
3348 t_soisno(c,0) = min(tfrz, forc_t(g)) ! K
3349 h2osoi_ice(c,0) = h2osno(c) ! kg/m2
3350 h2osoi_liq(c,0) = 0._kind_lake ! kg/m2
3351 frac_iceold(c,0) = 1._kind_lake
3352 end if
3353
3354 ! The change of ice partial density of surface node due to precipitation.
3355 ! Only ice part of snowfall is added here, the liquid part will be added
3356 ! later.
3357
3358 if (snl(c) < 0 .and. newnode == 0) then
3359 h2osoi_ice(c,snl(c)+1) = h2osoi_ice(c,snl(c)+1)+dtime*qflx_snow_grnd_col(c)
3360 dz(c,snl(c)+1) = dz(c,snl(c)+1)+dz_snowf*dtime
3361 end if
3362
3363 end do
3364
3365 ! Calculate sublimation and dew, adapted from HydrologyLake and Biogeophysics2.
3366
3367 !dir$ concurrent
3368 !cdir nodep
3369 do fp = 1,num_shlakep
3370 p = filter_shlakep(fp)
3371 c = pcolumn(p)
3372 jtop = snl(c)+1
3373
3374 ! Use column variables here
3375 qflx_evap_grnd(c) = 0._kind_lake
3376 qflx_sub_snow(c) = 0._kind_lake
3377 qflx_dew_snow(c) = 0._kind_lake
3378 qflx_dew_grnd(c) = 0._kind_lake
3379
3380 if (jtop <= 0) then ! snow layers
3381 j = jtop
3382 ! Assign ground evaporation to sublimation from soil ice or to dew
3383 ! on snow or ground
3384
3385 if (qflx_evap_soi(p) >= 0._kind_lake) then
3386 ! for evaporation partitioning between liquid evap and ice sublimation,
3387 ! use the ratio of liquid to (liquid+ice) in the top layer to determine split
3388 ! Since we're not limiting evap over lakes, but still can't remove more from top
3389 ! snow layer than there is there, create temp. limited evap_soi.
3390 qflx_evap_soi_lim = min(qflx_evap_soi(p), (h2osoi_liq(c,j)+h2osoi_ice(c,j))/dtime)
3391 if ((h2osoi_liq(c,j)+h2osoi_ice(c,j)) > 0._kind_lake) then
3392 qflx_evap_grnd(c) = max(qflx_evap_soi_lim*(h2osoi_liq(c,j)/(h2osoi_liq(c,j)+h2osoi_ice(c,j))), 0._kind_lake)
3393 else
3394 qflx_evap_grnd(c) = 0._kind_lake
3395 end if
3396 qflx_sub_snow(c) = qflx_evap_soi_lim - qflx_evap_grnd(c)
3397 else
3398 if (t_grnd(c) < tfrz) then
3399 qflx_dew_snow(c) = abs(qflx_evap_soi(p))
3400 else
3401 qflx_dew_grnd(c) = abs(qflx_evap_soi(p))
3402 end if
3403 end if
3404 ! Update the pft-level qflx_snowcap
3405 ! This was moved in from Hydrology2 to keep all pft-level
3406 ! calculations out of Hydrology2
3407 if (do_capsnow(c)) qflx_snowcap(p) = qflx_snowcap(p) + qflx_dew_snow(c) + qflx_dew_grnd(c)
3408
3409 else ! No snow layers: do as in HydrologyLake but with actual clmtype variables
3410 if (qflx_evap_soi(p) >= 0._kind_lake) then
3411 ! Sublimation: do not allow for more sublimation than there is snow
3412 ! after melt. Remaining surface evaporation used for infiltration.
3413 qflx_sub_snow(c) = min(qflx_evap_soi(p), h2osno(c)/dtime)
3414 qflx_evap_grnd(c) = qflx_evap_soi(p) - qflx_sub_snow(c)
3415 else
3416 if (t_grnd(c) < tfrz-0.1_kind_lake) then
3417 qflx_dew_snow(c) = abs(qflx_evap_soi(p))
3418 else
3419 qflx_dew_grnd(c) = abs(qflx_evap_soi(p))
3420 end if
3421 end if
3422
3423 ! Update snow pack for dew & sub.
3424 h2osno_temp = h2osno(c)
3425 if (do_capsnow(c)) then
3426 h2osno(c) = h2osno(c) - qflx_sub_snow(c)*dtime
3427 qflx_snowcap(p) = qflx_snowcap(p) + qflx_dew_snow(c) + qflx_dew_grnd(c)
3428 else
3429 h2osno(c) = h2osno(c) + (-qflx_sub_snow(c)+qflx_dew_snow(c))*dtime
3430 end if
3431 if (h2osno_temp > 0._kind_lake) then
3432 snowdp(c) = snowdp(c) * h2osno(c) / h2osno_temp
3433 else
3434 snowdp(c) = h2osno(c)/snow_bd !Assume a constant snow bulk density = 250.
3435 end if
3436
3437 if (pergro) then
3438 if (abs(h2osno(c)) < 1.e-10_kind_lake) h2osno(c) = 0._kind_lake
3439 else
3440 h2osno(c) = max(h2osno(c), 0._kind_lake)
3441 endif
3442
3443 end if
3444
3445 qflx_snowcap_col(c) = qflx_snowcap(p)
3446
3447 end do
3448
3449
3450 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3451 ! Determine initial snow/no-snow filters (will be modified possibly by
3452 ! routines CombineSnowLayers and DivideSnowLayers below
3453
3454 call buildsnowfilter(lbc, ubc, num_shlakec, filter_shlakec,snl, & !i
3455 num_shlakesnowc, filter_shlakesnowc, num_shlakenosnowc, filter_shlakenosnowc) !o
3456
3457 ! Determine the change of snow mass and the snow water onto soil
3458
3459 call snowwater(lbc, ubc, num_shlakesnowc, filter_shlakesnowc, & !i
3460 num_shlakenosnowc, filter_shlakenosnowc, & !i
3461 snl,do_capsnow,qflx_snomelt,qflx_rain_grnd, & !i
3462 qflx_sub_snow,qflx_evap_grnd, & !i
3463 qflx_dew_snow,qflx_dew_grnd,dz,dtime, & !i
3464 h2osoi_ice,h2osoi_liq, & !i&o
3465 qflx_top_soil) !o
3466
3467
3468 ! Determine soil hydrology
3469 ! Here this consists only of making sure that soil is saturated even as it melts and 10%
3470 ! of pore space opens up. Conversely, if excess ice is melting and the liquid water exceeds the
3471 ! saturation value, then remove water.
3472
3473 do j = 1,nlevsoil
3474 !dir$ concurrent
3475 !cdir nodep
3476 do fc = 1, num_shlakec
3477 c = filter_shlakec(fc)
3478
3479 if (h2osoi_vol(c,j) < watsat(c,j)) then
3480 h2osoi_liq(c,j) = (watsat(c,j)*dz(c,j) - h2osoi_ice(c,j)/denice)*denh2o
3481 ! h2osoi_vol will be updated below, and this water addition will come from qflx_qrgwl
3482 else if (h2osoi_liq(c,j) > watsat(c,j)*denh2o*dz(c,j)) then
3483 h2osoi_liq(c,j) = watsat(c,j)*denh2o*dz(c,j)
3484 end if
3485
3486 end do
3487 end do
3488 !!!!!!!!!!
3489
3490 ! if (.not. is_perpetual()) then
3491 if (1==1) then
3492
3493 ! Natural compaction and metamorphosis.
3494
3495 call snowcompaction(lbc, ubc, num_shlakesnowc, filter_shlakesnowc, &!i
3496 snl,imelt,frac_iceold,t_soisno, &!i
3497 h2osoi_ice,h2osoi_liq,dtime, &!i
3498 dz) !&o
3499
3500 ! Combine thin snow elements
3501
3502 call combinesnowlayers(lbc, ubc, & !i
3503 num_shlakesnowc, filter_shlakesnowc, & !i&o
3504 snl,h2osno,snowdp,dz,zi, & !i&o
3505 t_soisno,h2osoi_ice,h2osoi_liq, & !i&o
3506 z) !o
3507
3508
3509 ! Divide thick snow elements
3510
3511 call dividesnowlayers(lbc, ubc, & !i
3512 num_shlakesnowc, filter_shlakesnowc, & !i&o
3513 snl,dz,zi,t_soisno, & !i&o
3514 h2osoi_ice,h2osoi_liq, & !i&o
3515 z) !o
3516
3517
3518 else
3519
3520 do fc = 1, num_shlakesnowc
3521 c = filter_shlakesnowc(fc)
3522 h2osno(c) = 0._kind_lake
3523 end do
3524 do j = -nlevsnow+1,0
3525 do fc = 1, num_shlakesnowc
3526 c = filter_shlakesnowc(fc)
3527 if (j >= snl(c)+1) then
3528 h2osno(c) = h2osno(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
3529 end if
3530 end do
3531 end do
3532
3533 end if
3534
3535 ! Check for snow layers above lake with unfrozen top layer. Mechanically,
3536 ! the snow will fall into the lake and melt or turn to ice. If the top layer has
3537 ! sufficient heat to melt the snow without freezing, then that will be done.
3538 ! Otherwise, the top layer will undergo freezing, but only if the top layer will
3539 ! not freeze completely. Otherwise, let the snow layers persist and melt by diffusion.
3540 !dir$ concurrent
3541 !cdir nodep
3542 do fc = 1, num_shlakec
3543 c = filter_shlakec(fc)
3544
3545 if (t_lake(c,1) > tfrz .and. lake_icefrac(c,1) == 0._kind_lake .and. snl(c) < 0) then
3546 unfrozen(c) = .true.
3547 else
3548 unfrozen(c) = .false.
3549 end if
3550 end do
3551
3552 do j = -nlevsnow+1,0
3553 !dir$ concurrent
3554 !cdir nodep
3555 do fc = 1, num_shlakec
3556 c = filter_shlakec(fc)
3557
3558 if (unfrozen(c)) then
3559 if (j == -nlevsnow+1) then
3560 sumsnowice(c) = 0._kind_lake
3561 heatsum(c) = 0._kind_lake
3562 end if
3563 if (j >= snl(c)+1) then
3564 sumsnowice(c) = sumsnowice(c) + h2osoi_ice(c,j)
3565 heatsum(c) = heatsum(c) + h2osoi_ice(c,j)*cpice*(tfrz - t_soisno(c,j)) &
3566 + h2osoi_liq(c,j)*cpliq*(tfrz - t_soisno(c,j))
3567 end if
3568 end if
3569 end do
3570 end do
3571
3572 !dir$ concurrent
3573 !cdir nodep
3574 do fc = 1, num_shlakec
3575 c = filter_shlakec(fc)
3576
3577 if (unfrozen(c)) then
3578 heatsum(c) = heatsum(c) + sumsnowice(c)*hfus
3579 heatrem = (t_lake(c,1) - tfrz)*cpliq*denh2o*dz_lake(c,1) - heatsum(c)
3580
3581 if (heatrem + denh2o*dz_lake(c,1)*hfus > 0._kind_lake) then
3582 ! Remove snow and subtract the latent heat from the top layer.
3583 h2osno(c) = 0._kind_lake
3584 snl(c) = 0
3585 ! The rest of the bookkeeping for the removed snow will be done below.
3586 if (debug_print) then
3587 print *,'Snow layers removed above unfrozen lake for column, snowice:', &
3588 c, sumsnowice(c)
3589 endif
3590 if (heatrem > 0._kind_lake) then ! simply subtract the heat from the layer
3591 t_lake(c,1) = t_lake(c,1) - heatrem/(cpliq*denh2o*dz_lake(c,1))
3592 else !freeze part of the layer
3593 t_lake(c,1) = tfrz
3594 lake_icefrac(c,1) = -heatrem/(denh2o*dz_lake(c,1)*hfus)
3595 end if
3596 end if
3597 end if
3598 end do
3599 !!!!!!!!!!!!
3600
3601 ! Set snow age to zero if no snow
3602
3603 !dir$ concurrent
3604 !cdir nodep
3605 do fc = 1, num_shlakesnowc
3606 c = filter_shlakesnowc(fc)
3607 if (snl(c) == 0) then
3608 snowage(c) = 0._kind_lake
3609 end if
3610 end do
3611
3612 ! Set empty snow layers to zero
3613
3614 do j = -nlevsnow+1,0
3615 !dir$ concurrent
3616 !cdir nodep
3617 do fc = 1, num_shlakesnowc
3618 c = filter_shlakesnowc(fc)
3619 if (j <= snl(c) .and. snl(c) > -nlevsnow) then
3620 h2osoi_ice(c,j) = 0._kind_lake
3621 h2osoi_liq(c,j) = 0._kind_lake
3622 t_soisno(c,j) = 0._kind_lake
3623 dz(c,j) = 0._kind_lake
3624 z(c,j) = 0._kind_lake
3625 zi(c,j-1) = 0._kind_lake
3626 end if
3627 end do
3628 end do
3629
3630 ! Build new snow filter
3631
3632 call buildsnowfilter(lbc, ubc, num_shlakec, filter_shlakec, snl,& !i
3633 num_shlakesnowc, filter_shlakesnowc, num_shlakenosnowc, filter_shlakenosnowc) !o
3634
3635 ! Vertically average t_soisno and sum of h2osoi_liq and h2osoi_ice
3636 ! over all snow layers for history output
3637
3638 !dir$ concurrent
3639 !cdir nodep
3640 do fc = 1, num_shlakesnowc
3641 c = filter_shlakesnowc(fc)
3642 t_snow(c) = 0._kind_lake
3643 snowice(c) = 0._kind_lake
3644 snowliq(c) = 0._kind_lake
3645 end do
3646 !dir$ concurrent
3647 !cdir nodep
3648 do fc = 1, num_shlakenosnowc
3649 c = filter_shlakenosnowc(fc)
3650 t_snow(c) = spval
3651 snowice(c) = spval
3652 snowliq(c) = spval
3653 end do
3654
3655 do j = -nlevsnow+1, 0
3656 !dir$ concurrent
3657 !cdir nodep
3658 do fc = 1, num_shlakesnowc
3659 c = filter_shlakesnowc(fc)
3660 if (j >= snl(c)+1) then
3661 t_snow(c) = t_snow(c) + t_soisno(c,j)
3662 snowice(c) = snowice(c) + h2osoi_ice(c,j)
3663 snowliq(c) = snowliq(c) + h2osoi_liq(c,j)
3664 end if
3665 end do
3666 end do
3667
3668 ! Determine ending water balance and volumetric soil water
3669
3670 !dir$ concurrent
3671 !cdir nodep
3672 do fc = 1, num_shlakec
3673
3674 c = filter_shlakec(fc)
3675 if (snl(c) < 0) t_snow(c) = t_snow(c)/abs(snl(c))
3676 endwb(c) = h2osno(c)
3677 end do
3678
3679 do j = 1, nlevsoil
3680 !dir$ concurrent
3681 !cdir nodep
3682 do fc = 1, num_shlakec
3683 c = filter_shlakec(fc)
3684 endwb(c) = endwb(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
3685 h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice)
3686 end do
3687 end do
3688
3689 check_add_snow_water: if(lakedebug) then
3690 allocate(snow_water(lbc:ubc))
3691 ! Check to make sure snow water adds up correctly.
3692 do j = -nlevsnow+1,0
3693 !dir$ concurrent
3694 !cdir nodep
3695 do fc = 1, num_shlakec
3696 c = filter_shlakec(fc)
3697
3698 jtop = snl(c)+1
3699 if(j == jtop) snow_water(c) = 0._kind_lake
3700 if(j >= jtop) then
3701 snow_water(c) = snow_water(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
3702 if(j == 0 .and. abs(snow_water(c)-h2osno(c))>1.e-7_kind_lake) then
3703 write(message,*)'h2osno does not equal sum of snow layers in ShalLakeHydrology:', &
3704 'column, h2osno, sum of snow layers =', c, h2osno(c), snow_water(c)
3705 ! errmsg=trim(message)
3706 ! errflg=1
3707 write(0,'(A)') trim(message)
3708 end if
3709 end if
3710 end do
3711 end do
3712 deallocate(snow_water)
3713 end if check_add_snow_water
3714
3715 !!!!!!!!!!!!!
3716 ! Do history variables and set special landunit runoff (adapted from end of HydrologyLake)
3717 !dir$ concurrent
3718 !cdir nodep
3719 do fp = 1,num_shlakep
3720 p = filter_shlakep(fp)
3721 c = pcolumn(p)
3722 g = pgridcell(p)
3723
3724 qflx_infl(c) = 0._kind_lake
3725 qflx_surf(c) = 0._kind_lake
3726 qflx_drain(c) = 0._kind_lake
3727 rootr_column(c,:) = spval
3728 soilalpha(c) = spval
3729 zwt(c) = spval
3730 fcov(c) = spval
3731 qcharge(c) = spval
3732 ! h2osoi_vol(c,:) = spval
3733
3734 ! Insure water balance using qflx_qrgwl
3735 qflx_qrgwl(c) = forc_rain(g) + forc_snow(g) - qflx_evap_tot(p) - (endwb(c)-begwb(c))/dtime
3736 if (debug_print) then
3737 print *,'c, rain, snow, evap, endwb, begwb, qflx_qrgwl:', &
3738 c, forc_rain(g), forc_snow(g), qflx_evap_tot(p), endwb(c), begwb(c), qflx_qrgwl(c)
3739 endif
3740
3741 ! The pft average must be done here for output to history tape
3742 qflx_evap_tot_col(c) = qflx_evap_tot(p)
3743 end do
3744
3745 end subroutine shallakehydrology
3746
3747! DESCRIPTION:
3750 subroutine qsat (T, p, es, esdT, qs, qsdT)
3751 !
3752 ! Reference: Polynomial approximations from:
3753 ! Piotr J. Flatau, et al.,1992: Polynomial fits to saturation
3754 ! vapor pressure. Journal of Applied Meteorology, 31, 1507-1513.
3755 !
3756 ! USES:
3757 !
3758 ! ARGUMENTS:
3759 implicit none
3760 real(kind_lake), intent(in) :: T
3761 real(kind_lake), intent(in) :: p
3762 real(kind_lake), intent(out) :: es
3763 real(kind_lake), intent(out) :: esdT
3764 real(kind_lake), intent(out) :: qs
3765 real(kind_lake), intent(out) :: qsdT
3766 !
3767 ! !CALLED FROM:
3768 ! subroutine Biogeophysics1 in module Biogeophysics1Mod
3769 ! subroutine BiogeophysicsLake in module BiogeophysicsLakeMod
3770 ! subroutine CanopyFluxesMod CanopyFluxesMod
3771 !
3772 ! !REVISION HISTORY:
3773 ! 15 September 1999: Yongjiu Dai; Initial code
3774 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
3775 !
3776 !EOP
3777 !
3778 ! !LOCAL VARIABLES:
3779 !
3780 real(kind_lake) :: T_limit
3781 real(kind_lake) :: td,vp,vp1,vp2
3782 !
3783 ! For water vapor (temperature range 0C-100C)
3784 !
3785 real(kind_lake), parameter :: a0 = 6.11213476
3786 real(kind_lake), parameter :: a1 = 0.444007856
3787 real(kind_lake), parameter :: a2 = 0.143064234e-01
3788 real(kind_lake), parameter :: a3 = 0.264461437e-03
3789 real(kind_lake), parameter :: a4 = 0.305903558e-05
3790 real(kind_lake), parameter :: a5 = 0.196237241e-07
3791 real(kind_lake), parameter :: a6 = 0.892344772e-10
3792 real(kind_lake), parameter :: a7 = -0.373208410e-12
3793 real(kind_lake), parameter :: a8 = 0.209339997e-15
3794 !
3795 ! For derivative:water vapor
3796 !
3797 real(kind_lake), parameter :: b0 = 0.444017302
3798 real(kind_lake), parameter :: b1 = 0.286064092e-01
3799 real(kind_lake), parameter :: b2 = 0.794683137e-03
3800 real(kind_lake), parameter :: b3 = 0.121211669e-04
3801 real(kind_lake), parameter :: b4 = 0.103354611e-06
3802 real(kind_lake), parameter :: b5 = 0.404125005e-09
3803 real(kind_lake), parameter :: b6 = -0.788037859e-12
3804 real(kind_lake), parameter :: b7 = -0.114596802e-13
3805 real(kind_lake), parameter :: b8 = 0.381294516e-16
3806 !
3807 ! For ice (temperature range -75C-0C)
3808 !
3809 real(kind_lake), parameter :: c0 = 6.11123516
3810 real(kind_lake), parameter :: c1 = 0.503109514
3811 real(kind_lake), parameter :: c2 = 0.188369801e-01
3812 real(kind_lake), parameter :: c3 = 0.420547422e-03
3813 real(kind_lake), parameter :: c4 = 0.614396778e-05
3814 real(kind_lake), parameter :: c5 = 0.602780717e-07
3815 real(kind_lake), parameter :: c6 = 0.387940929e-09
3816 real(kind_lake), parameter :: c7 = 0.149436277e-11
3817 real(kind_lake), parameter :: c8 = 0.262655803e-14
3818 !
3819 ! For derivative:ice
3820 !
3821 real(kind_lake), parameter :: d0 = 0.503277922
3822 real(kind_lake), parameter :: d1 = 0.377289173e-01
3823 real(kind_lake), parameter :: d2 = 0.126801703e-02
3824 real(kind_lake), parameter :: d3 = 0.249468427e-04
3825 real(kind_lake), parameter :: d4 = 0.313703411e-06
3826 real(kind_lake), parameter :: d5 = 0.257180651e-08
3827 real(kind_lake), parameter :: d6 = 0.133268878e-10
3828 real(kind_lake), parameter :: d7 = 0.394116744e-13
3829 real(kind_lake), parameter :: d8 = 0.498070196e-16
3830 !-----------------------------------------------------------------------
3831
3832 t_limit = t - tfrz
3833 if (t_limit > 100.0) t_limit=100.0
3834 if (t_limit < -75.0) t_limit=-75.0
3835
3836 td = t_limit
3837 if (td >= 0.0) then
3838 es = a0 + td*(a1 + td*(a2 + td*(a3 + td*(a4 &
3839 + td*(a5 + td*(a6 + td*(a7 + td*a8)))))))
3840 esdt = b0 + td*(b1 + td*(b2 + td*(b3 + td*(b4 &
3841 + td*(b5 + td*(b6 + td*(b7 + td*b8)))))))
3842 else
3843 es = c0 + td*(c1 + td*(c2 + td*(c3 + td*(c4 &
3844 + td*(c5 + td*(c6 + td*(c7 + td*c8)))))))
3845 esdt = d0 + td*(d1 + td*(d2 + td*(d3 + td*(d4 &
3846 + td*(d5 + td*(d6 + td*(d7 + td*d8)))))))
3847 endif
3848
3849 es = es * 100. ! pa
3850 esdt = esdt * 100. ! pa/K
3851
3852 vp = 1.0 / (p - one_minus_con_eps*es)
3853 vp1 = 0.622 * vp
3854 vp2 = vp1 * vp
3855
3856 qs = es * vp1 ! kg/kg
3857 qsdt = esdt * vp2 * p ! 1 / K
3858
3859 end subroutine qsat
3860
3861
3862 subroutine tridiagonal (lbc, ubc, lbj, ubj, jtop, numf, filter, &
3863 a, b, c, r, u)
3864 !
3865 ! DESCRIPTION:
3866
3867 !
3868 ! ARGUMENTS:
3869 implicit none
3870 integer , intent(in) :: lbc, ubc
3871 integer , intent(in) :: lbj, ubj
3872 integer , intent(in) :: jtop(lbc:ubc)
3873 integer , intent(in) :: numf
3874 integer , intent(in) :: filter(1:numf)
3875 real(kind_lake), intent(in) :: a(lbc:ubc, lbj:ubj)
3876 real(kind_lake), intent(in) :: b(lbc:ubc, lbj:ubj)
3877 real(kind_lake), intent(in) :: c(lbc:ubc, lbj:ubj)
3878 real(kind_lake), intent(in) :: r(lbc:ubc, lbj:ubj)
3879 real(kind_lake), intent(inout) :: u(lbc:ubc, lbj:ubj)
3880 !
3881 ! !CALLED FROM:
3882 ! subroutine BiogeophysicsLake in module BiogeophysicsLakeMod
3883 ! subroutine SoilTemperature in module SoilTemperatureMod
3884 ! subroutine SoilWater in module HydrologyMod
3885 !
3886 ! !REVISION HISTORY:
3887 ! 15 September 1999: Yongjiu Dai; Initial code
3888 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
3889 ! 1 July 2003: Mariana Vertenstein; modified for vectorization
3890 !
3891 !EOP
3892 !
3893 ! !OTHER LOCAL VARIABLES:
3894 !
3895 integer :: j,ci,fc !indices
3896 real(kind_lake) :: gam(lbc:ubc,lbj:ubj) !temporary
3897 real(kind_lake) :: bet(lbc:ubc) !temporary
3898 !-----------------------------------------------------------------------
3899
3900 ! Solve the matrix
3901
3902 !dir$ concurrent
3903 !cdir nodep
3904 do fc = 1,numf
3905 ci = filter(fc)
3906 bet(ci) = b(ci,jtop(ci))
3907 end do
3908
3909 do j = lbj, ubj
3910 !dir$ prefervector
3911 !dir$ concurrent
3912 !cdir nodep
3913 do fc = 1,numf
3914 ci = filter(fc)
3915 if (j >= jtop(ci)) then
3916 if (j == jtop(ci)) then
3917 u(ci,j) = r(ci,j) / bet(ci)
3918 else
3919 gam(ci,j) = c(ci,j-1) / bet(ci)
3920 bet(ci) = b(ci,j) - a(ci,j) * gam(ci,j)
3921 u(ci,j) = (r(ci,j) - a(ci,j)*u(ci,j-1)) / bet(ci)
3922 end if
3923 end if
3924 end do
3925 end do
3926
3927 !Cray X1 unroll directive used here as work-around for compiler issue 2003/10/20
3928 !dir$ unroll 0
3929 do j = ubj-1,lbj,-1
3930 !dir$ prefervector
3931 !dir$ concurrent
3932 !cdir nodep
3933 do fc = 1,numf
3934 ci = filter(fc)
3935 if (j >= jtop(ci)) then
3936 u(ci,j) = u(ci,j) - gam(ci,j+1) * u(ci,j+1)
3937 end if
3938 end do
3939 end do
3940
3941 end subroutine tridiagonal
3942
3943
3944 ! DESCRIPTION:
3955 subroutine snowwater(lbc, ubc, num_snowc, filter_snowc, & !i
3956 num_nosnowc, filter_nosnowc, & !i
3957 snl,do_capsnow,qflx_snomelt,qflx_rain_grnd, & !i
3958 qflx_sub_snow,qflx_evap_grnd, & !i
3959 qflx_dew_snow,qflx_dew_grnd,dz,dtime, & !i
3960 h2osoi_ice,h2osoi_liq, & !i&o
3961 qflx_top_soil) !o
3962 !===============================================================================
3963 ! !REVISION HISTORY:
3964 ! 15 September 1999: Yongjiu Dai; Initial code
3965 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
3966 ! 15 November 2000: Mariana Vertenstein
3967 ! 2/26/02, Peter Thornton: Migrated to new data structures.
3968 !=============================================================================
3969 ! !USES:
3970 ! use clmtype
3971
3972 implicit none
3973
3974 !in:
3975 integer, intent(in) :: lbc, ubc
3976 integer, intent(in) :: num_snowc
3977 integer, intent(in) :: filter_snowc(ubc-lbc+1)
3978 integer, intent(in) :: num_nosnowc
3979 integer, intent(in) :: filter_nosnowc(ubc-lbc+1)
3980
3981 integer , intent(in) :: snl(1)
3982 logical , intent(in) :: do_capsnow(1)
3983 real(kind_lake), intent(in) :: dtime
3984 real(kind_lake), intent(in) :: qflx_snomelt(1)
3985 real(kind_lake), intent(in) :: qflx_rain_grnd(1)
3986 real(kind_lake), intent(in) :: qflx_sub_snow(1)
3987 real(kind_lake), intent(in) :: qflx_evap_grnd(1)
3988 real(kind_lake), intent(in) :: qflx_dew_snow(1)
3989 real(kind_lake), intent(in) :: qflx_dew_grnd(1)
3990 real(kind_lake), intent(in) :: dz(1,-nlevsnow+1:nlevsoil)
3991
3992
3993 !inout:
3994
3995 real(kind_lake), intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil)
3996 real(kind_lake), intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil)
3997
3998 !out:
3999
4000 real(kind_lake), intent(out) :: qflx_top_soil(1)
4001
4002
4003 ! OTHER LOCAL VARIABLES:
4004
4005 integer :: c, j, fc !do loop/array indices
4006 real(kind_lake) :: qin(lbc:ubc) !water flow into the elmement (mm/s)
4007 real(kind_lake) :: qout(lbc:ubc) !water flow out of the elmement (mm/s)
4008 real(kind_lake) :: wgdif !ice mass after minus sublimation
4009 real(kind_lake) :: vol_liq(lbc:ubc,-nlevsnow+1:0) !partial volume of liquid water in layer
4010 real(kind_lake) :: vol_ice(lbc:ubc,-nlevsnow+1:0) !partial volume of ice lens in layer
4011 real(kind_lake) :: eff_porosity(lbc:ubc,-nlevsnow+1:0) !effective porosity = porosity - vol_ice
4012 !-----------------------------------------------------------------------
4013 ! Renew the mass of ice lens (h2osoi_ice) and liquid (h2osoi_liq) in the
4014 ! surface snow layer resulting from sublimation (frost) / evaporation (condense)
4015
4016 !dir$ concurrent
4017 !cdir nodep
4018 do fc = 1,num_snowc
4019 c = filter_snowc(fc)
4020 if (do_capsnow(c)) then
4021 wgdif = h2osoi_ice(c,snl(c)+1) - qflx_sub_snow(c)*dtime
4022 h2osoi_ice(c,snl(c)+1) = wgdif
4023 if (wgdif < 0.) then
4024 h2osoi_ice(c,snl(c)+1) = 0.
4025 h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) + wgdif
4026 end if
4027 h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) - qflx_evap_grnd(c) * dtime
4028 else
4029 wgdif = h2osoi_ice(c,snl(c)+1) + (qflx_dew_snow(c) - qflx_sub_snow(c)) * dtime
4030 h2osoi_ice(c,snl(c)+1) = wgdif
4031 if (wgdif < 0.) then
4032 h2osoi_ice(c,snl(c)+1) = 0.
4033 h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) + wgdif
4034 end if
4035 h2osoi_liq(c,snl(c)+1) = h2osoi_liq(c,snl(c)+1) + &
4036 (qflx_rain_grnd(c) + qflx_dew_grnd(c) - qflx_evap_grnd(c)) * dtime
4037 end if
4038 h2osoi_liq(c,snl(c)+1) = max(0._kind_lake, h2osoi_liq(c,snl(c)+1))
4039 end do
4040
4041 ! Porosity and partial volume
4042
4043 do j = -nlevsnow+1, 0
4044 !dir$ concurrent
4045 !cdir nodep
4046 do fc = 1, num_snowc
4047 c = filter_snowc(fc)
4048 if (j >= snl(c)+1) then
4049 vol_ice(c,j) = min(1._kind_lake, h2osoi_ice(c,j)/(dz(c,j)*denice))
4050 eff_porosity(c,j) = 1. - vol_ice(c,j)
4051 vol_liq(c,j) = min(eff_porosity(c,j),h2osoi_liq(c,j)/(dz(c,j)*denh2o))
4052 end if
4053 end do
4054 end do
4055
4056 ! Capillary forces within snow are usually two or more orders of magnitude
4057 ! less than those of gravity. Only gravity terms are considered.
4058 ! the genernal expression for water flow is "K * ss**3", however,
4059 ! no effective parameterization for "K". Thus, a very simple consideration
4060 ! (not physically based) is introduced:
4061 ! when the liquid water of layer exceeds the layer's holding
4062 ! capacity, the excess meltwater adds to the underlying neighbor layer.
4063
4064 qin(:) = 0._kind_lake
4065
4066 do j = -nlevsnow+1, 0
4067 !dir$ concurrent
4068 !cdir nodep
4069 do fc = 1, num_snowc
4070 c = filter_snowc(fc)
4071 if (j >= snl(c)+1) then
4072 h2osoi_liq(c,j) = h2osoi_liq(c,j) + qin(c)
4073 if (j <= -1) then
4074 ! No runoff over snow surface, just ponding on surface
4075 if (eff_porosity(c,j) < wimp .OR. eff_porosity(c,j+1) < wimp) then
4076 qout(c) = 0._kind_lake
4077 else
4078 qout(c) = max(0._kind_lake,(vol_liq(c,j)-ssi*eff_porosity(c,j))*dz(c,j))
4079 qout(c) = min(qout(c),(1.-vol_ice(c,j+1)-vol_liq(c,j+1))*dz(c,j+1))
4080 end if
4081 else
4082 qout(c) = max(0._kind_lake,(vol_liq(c,j) - ssi*eff_porosity(c,j))*dz(c,j))
4083 end if
4084 qout(c) = qout(c)*1000.
4085 h2osoi_liq(c,j) = h2osoi_liq(c,j) - qout(c)
4086 qin(c) = qout(c)
4087 end if
4088 end do
4089 end do
4090
4091 !dir$ concurrent
4092 !cdir nodep
4093 do fc = 1, num_snowc
4094 c = filter_snowc(fc)
4095 ! Qout from snow bottom
4096 qflx_top_soil(c) = qout(c) / dtime
4097 end do
4098
4099 !dir$ concurrent
4100 !cdir nodep
4101 do fc = 1, num_nosnowc
4102 c = filter_nosnowc(fc)
4103 qflx_top_soil(c) = qflx_rain_grnd(c) + qflx_snomelt(c)
4104 end do
4105
4106 end subroutine snowwater
4107
4115 subroutine snowcompaction(lbc, ubc, num_snowc, filter_snowc, &!i
4116 snl,imelt,frac_iceold,t_soisno, &!i
4117 h2osoi_ice,h2osoi_liq,dtime, &!i
4118 dz) !i&o
4119
4120
4121 !================================================================================
4122 ! CALLED FROM:
4123 ! subroutine Hydrology2 in module Hydrology2Mod
4124 !
4125 ! REVISION HISTORY:
4126 ! 15 September 1999: Yongjiu Dai; Initial code
4127 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
4128 ! 2/28/02, Peter Thornton: Migrated to new data structures
4129 !==============================================================================
4130 ! USES:
4131 ! use clmtype
4132 !
4133 ! !ARGUMENTS:
4134 implicit none
4135
4136 !in:
4137 integer, intent(in) :: lbc, ubc
4138 integer, intent(in) :: num_snowc
4139 integer, intent(in) :: filter_snowc(ubc-lbc+1)
4140 integer, intent(in) :: snl(1)
4141 integer, intent(in) :: imelt(1,-nlevsnow+1:nlevsoil)
4142 real(kind_lake), intent(in) :: dtime
4143 real(kind_lake), intent(in) :: frac_iceold(1,-nlevsnow+1:nlevsoil)
4144 real(kind_lake), intent(in) :: t_soisno(1,-nlevsnow+1:nlevsoil)
4145 real(kind_lake), intent(in) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil)
4146 real(kind_lake), intent(in) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil)
4147
4148 !inout:
4149
4150 real(kind_lake), intent(inout) :: dz(1,-nlevsnow+1:nlevsoil)
4151
4152 ! OTHER LOCAL VARIABLES:
4153
4154 integer :: j, c, fc ! indices
4155 real(kind_lake), parameter :: c2 = 23.e-3 ! [m3/kg]
4156 real(kind_lake), parameter :: c3 = 2.777e-6 ! [1/s]
4157 real(kind_lake), parameter :: c4 = 0.04 ! [1/K]
4158 real(kind_lake), parameter :: c5 = 2.0 !
4159 real(kind_lake), parameter :: dm = 100.0 ! Upper Limit on Destructive Metamorphism Compaction [kg/m3]
4160 real(kind_lake), parameter :: eta0 = 9.e+5 ! The Viscosity Coefficient Eta0 [kg-s/m2]
4161 real(kind_lake) :: burden(lbc:ubc) ! pressure of overlying snow [kg/m2]
4162 real(kind_lake) :: ddz1 ! Rate of settling of snowpack due to destructive metamorphism.
4163 real(kind_lake) :: ddz2 ! Rate of compaction of snowpack due to overburden.
4164 real(kind_lake) :: ddz3 ! Rate of compaction of snowpack due to melt [1/s]
4165 real(kind_lake) :: dexpf ! expf=exp(-c4*(273.15-t_soisno)).
4166 real(kind_lake) :: fi ! Fraction of ice relative to the total water content at current time step
4167 real(kind_lake) :: td ! t_soisno - tfrz [K]
4168 real(kind_lake) :: pdzdtc ! Nodal rate of change in fractional-thickness due to compaction [fraction/s]
4169 real(kind_lake) :: void ! void (1 - vol_ice - vol_liq)
4170 real(kind_lake) :: wx ! water mass (ice+liquid) [kg/m2]
4171 real(kind_lake) :: bi ! partial density of ice [kg/m3]
4172
4173 !-----------------------------------------------------------------------
4174
4175
4176 ! Begin calculation - note that the following column loops are only invoked if snl(c) < 0
4177
4178 burden(:) = 0._kind_lake
4179
4180 do j = -nlevsnow+1, 0
4181 !dir$ concurrent
4182 !cdir nodep
4183 do fc = 1, num_snowc
4184 c = filter_snowc(fc)
4185 if (j >= snl(c)+1) then
4186
4187 wx = h2osoi_ice(c,j) + h2osoi_liq(c,j)
4188 void = 1. - (h2osoi_ice(c,j)/denice + h2osoi_liq(c,j)/denh2o) / dz(c,j)
4189
4190 ! Allow compaction only for non-saturated node and higher ice lens node.
4191 if (void > 0.001 .and. h2osoi_ice(c,j) > .1) then
4192 bi = h2osoi_ice(c,j) / dz(c,j)
4193 fi = h2osoi_ice(c,j) / wx
4194 td = tfrz-t_soisno(c,j)
4195 dexpf = exp(-c4*td)
4196
4197 ! Settling as a result of destructive metamorphism
4198
4199 ddz1 = -c3*dexpf
4200 if (bi > dm) ddz1 = ddz1*exp(-46.0e-3*(bi-dm))
4201
4202 ! Liquid water term
4203
4204 if (h2osoi_liq(c,j) > 0.01*dz(c,j)) ddz1=ddz1*c5
4205
4206 ! Compaction due to overburden
4207
4208 ddz2 = -burden(c)*exp(-0.08*td - c2*bi)/eta0
4209
4210 ! Compaction occurring during melt
4211
4212 if (imelt(c,j) == 1) then
4213 ddz3 = - 1./dtime * max(0._kind_lake,(frac_iceold(c,j) - fi)/frac_iceold(c,j))
4214 else
4215 ddz3 = 0._kind_lake
4216 end if
4217
4218 ! Time rate of fractional change in dz (units of s-1)
4219
4220 pdzdtc = ddz1 + ddz2 + ddz3
4221
4222 ! The change in dz due to compaction
4223
4224 dz(c,j) = dz(c,j) * (1.+pdzdtc*dtime)
4225 end if
4226
4227 ! Pressure of overlying snow
4228
4229 burden(c) = burden(c) + wx
4230
4231 end if
4232 end do
4233 end do
4234
4235 end subroutine snowcompaction
4236
4240 subroutine combinesnowlayers(lbc, ubc, & !i
4241 num_snowc, filter_snowc, & !i&o
4242 snl,h2osno,snowdp,dz,zi, & !i&o
4243 t_soisno,h2osoi_ice,h2osoi_liq, & !i&o
4244 z) !o
4245 !==========================================================================
4246 ! DESCRIPTION:
4247 ! Combine snow layers that are less than a minimum thickness or mass
4248 ! If the snow element thickness or mass is less than a prescribed minimum,
4249 ! then it is combined with a neighboring element. The subroutine
4250 ! clm\_combo.f90 then executes the combination of mass and energy.
4251 ! CALLED FROM:
4252 ! subroutine Hydrology2 in module Hydrology2Mod
4253 !
4254 ! REVISION HISTORY:
4255 ! 15 September 1999: Yongjiu Dai; Initial code
4256 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
4257 ! 2/28/02, Peter Thornton: Migrated to new data structures.
4258 !=========================================================================
4259 ! !USES:
4260 ! use clmtype
4261 !
4262 ! !ARGUMENTS:
4263 implicit none
4264 !in:
4265 integer, intent(in) :: lbc, ubc
4266 ! integer, intent(in) :: clandunit(1) !landunit index for each column
4267 ! integer, intent(in) :: ityplun(1) !landunit type
4268
4269 !inout:
4270 integer, intent(inout) :: num_snowc
4271 integer, intent(inout) :: filter_snowc(ubc-lbc+1)
4272 integer , intent(inout) :: snl(1)
4273 real(kind_lake), intent(inout) :: h2osno(1)
4274 real(kind_lake), intent(inout) :: snowdp(1)
4275 real(kind_lake), intent(inout) :: dz(1,-nlevsnow+1:nlevsoil)
4276 real(kind_lake), intent(inout) :: zi(1,-nlevsnow+0:nlevsoil)
4277 real(kind_lake), intent(inout) :: t_soisno(1,-nlevsnow+1:nlevsoil)
4278 real(kind_lake), intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil)
4279 real(kind_lake), intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil)
4280
4281 !out:
4282
4283 real(kind_lake), intent(out) :: z(1,-nlevsnow+1:nlevsoil)
4284 !
4285 !EOP
4286 !
4287 ! !OTHER LOCAL VARIABLES:
4288 !
4289 integer :: c, fc ! column indices
4290 integer :: i,k ! loop indices
4291 integer :: j,l ! node indices
4292 integer :: msn_old(lbc:ubc) ! number of top snow layer
4293 integer :: mssi(lbc:ubc) ! node index
4294 integer :: neibor ! adjacent node selected for combination
4295 real(kind_lake):: zwice(lbc:ubc) ! total ice mass in snow
4296 real(kind_lake):: zwliq (lbc:ubc) ! total liquid water in snow
4297 real(kind_lake), parameter :: dzmin(5) = & ! minimum of top snow layer
4298 (/0.010, 0.015, 0.025, 0.055, 0.115/)
4299 !-----------------------------------------------------------------------
4300
4301 ! Check the mass of ice lens of snow, when the total is less than a small value,
4302 ! combine it with the underlying neighbor.
4303
4304 !dir$ concurrent
4305 !cdir nodep
4306 do fc = 1, num_snowc
4307 c = filter_snowc(fc)
4308 msn_old(c) = snl(c)
4309 end do
4310
4311 ! The following loop is NOT VECTORIZED
4312
4313 do fc = 1, num_snowc
4314 c = filter_snowc(fc)
4315 ! l = clandunit(c)
4316 do j = msn_old(c)+1,0
4317 if (h2osoi_ice(c,j) <= .1) then
4318 ! if (ityplun(l) == istsoil) then
4319 ! h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + h2osoi_liq(c,j)
4320 ! h2osoi_ice(c,j+1) = h2osoi_ice(c,j+1) + h2osoi_ice(c,j)
4321 ! else if (ityplun(l) /= istsoil .and. j /= 0) then
4322 h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + h2osoi_liq(c,j)
4323 h2osoi_ice(c,j+1) = h2osoi_ice(c,j+1) + h2osoi_ice(c,j)
4324 ! end if
4325
4326 ! shift all elements above this down one.
4327 if (j > snl(c)+1 .and. snl(c) < -1) then
4328 do i = j, snl(c)+2, -1
4329 t_soisno(c,i) = t_soisno(c,i-1)
4330 h2osoi_liq(c,i) = h2osoi_liq(c,i-1)
4331 h2osoi_ice(c,i) = h2osoi_ice(c,i-1)
4332 dz(c,i) = dz(c,i-1)
4333 end do
4334 end if
4335 snl(c) = snl(c) + 1
4336 end if
4337 end do
4338 end do
4339
4340 !dir$ concurrent
4341 !cdir nodep
4342 do fc = 1, num_snowc
4343 c = filter_snowc(fc)
4344 h2osno(c) = 0._kind_lake
4345 snowdp(c) = 0._kind_lake
4346 zwice(c) = 0._kind_lake
4347 zwliq(c) = 0._kind_lake
4348 end do
4349
4350 do j = -nlevsnow+1,0
4351 !dir$ concurrent
4352 !cdir nodep
4353 do fc = 1, num_snowc
4354 c = filter_snowc(fc)
4355 if (j >= snl(c)+1) then
4356 h2osno(c) = h2osno(c) + h2osoi_ice(c,j) + h2osoi_liq(c,j)
4357 snowdp(c) = snowdp(c) + dz(c,j)
4358 zwice(c) = zwice(c) + h2osoi_ice(c,j)
4359 zwliq(c) = zwliq(c) + h2osoi_liq(c,j)
4360 end if
4361 end do
4362 end do
4363
4364 ! Check the snow depth - all snow gone
4365 ! The liquid water assumes ponding on soil surface.
4366
4367 !dir$ concurrent
4368 !cdir nodep
4369 do fc = 1, num_snowc
4370 c = filter_snowc(fc)
4371 ! l = clandunit(c)
4372 if (snowdp(c) < 0.01 .and. snowdp(c) > 0.) then
4373 snl(c) = 0
4374 h2osno(c) = zwice(c)
4375 if (h2osno(c) <= 0.) snowdp(c) = 0._kind_lake
4376 ! if (ityplun(l) == istsoil) h2osoi_liq(c,1) = h2osoi_liq(c,1) + zwliq(c) !change by guhp
4377 end if
4378 end do
4379
4380 ! Check the snow depth - snow layers combined
4381 ! The following loop IS NOT VECTORIZED
4382
4383 do fc = 1, num_snowc
4384 c = filter_snowc(fc)
4385
4386 ! Two or more layers
4387
4388 if (snl(c) < -1) then
4389
4390 msn_old(c) = snl(c)
4391 mssi(c) = 1
4392
4393 do i = msn_old(c)+1,0
4394 if (dz(c,i) < dzmin(mssi(c))) then
4395
4396 if (i == snl(c)+1) then
4397 ! If top node is removed, combine with bottom neighbor.
4398 neibor = i + 1
4399 else if (i == 0) then
4400 ! If the bottom neighbor is not snow, combine with the top neighbor.
4401 neibor = i - 1
4402 else
4403 ! If none of the above special cases apply, combine with the thinnest neighbor
4404 neibor = i + 1
4405 if ((dz(c,i-1)+dz(c,i)) < (dz(c,i+1)+dz(c,i))) neibor = i-1
4406 end if
4407
4408 ! Node l and j are combined and stored as node j.
4409 if (neibor > i) then
4410 j = neibor
4411 l = i
4412 else
4413 j = i
4414 l = neibor
4415 end if
4416
4417 call combo (dz(c,j), h2osoi_liq(c,j), h2osoi_ice(c,j), &
4418 t_soisno(c,j), dz(c,l), h2osoi_liq(c,l), h2osoi_ice(c,l), t_soisno(c,l) )
4419
4420 ! Now shift all elements above this down one.
4421 if (j-1 > snl(c)+1) then
4422 do k = j-1, snl(c)+2, -1
4423 t_soisno(c,k) = t_soisno(c,k-1)
4424 h2osoi_ice(c,k) = h2osoi_ice(c,k-1)
4425 h2osoi_liq(c,k) = h2osoi_liq(c,k-1)
4426 dz(c,k) = dz(c,k-1)
4427 end do
4428 end if
4429
4430 ! Decrease the number of snow layers
4431 snl(c) = snl(c) + 1
4432 if (snl(c) >= -1) EXIT
4433
4434 else
4435
4436 ! The layer thickness is greater than the prescribed minimum value
4437 mssi(c) = mssi(c) + 1
4438
4439 end if
4440 end do
4441
4442 end if
4443
4444 end do
4445
4446 ! Reset the node depth and the depth of layer interface
4447
4448 do j = 0, -nlevsnow+1, -1
4449 !dir$ concurrent
4450 !cdir nodep
4451 do fc = 1, num_snowc
4452 c = filter_snowc(fc)
4453 if (j >= snl(c) + 1) then
4454 z(c,j) = zi(c,j) - 0.5*dz(c,j)
4455 zi(c,j-1) = zi(c,j) - dz(c,j)
4456 end if
4457 end do
4458 end do
4459
4460 end subroutine combinesnowlayers
4461
4462! DESCRIPTION:
4464 subroutine dividesnowlayers(lbc, ubc, & !i
4465 num_snowc, filter_snowc, & !i&o
4466 snl,dz,zi,t_soisno, & !i&o
4467 h2osoi_ice,h2osoi_liq, & !i&o
4468 z) !o
4469
4470
4471 !============================================================================
4472 ! !CALLED FROM:
4473 ! subroutine Hydrology2 in module Hydrology2Mod
4474 !
4475 ! !REVISION HISTORY:
4476 ! 15 September 1999: Yongjiu Dai; Initial code
4477 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
4478 ! 2/28/02, Peter Thornton: Migrated to new data structures.
4479 !============================================================================
4480 ! !USES:
4481 ! use clmtype
4482 !
4483 ! !ARGUMENTS:
4484 implicit none
4485
4486 !in:
4487 integer, intent(in) :: lbc, ubc
4488
4489 !inout:
4490
4491 integer, intent(inout) :: num_snowc
4492 integer, intent(inout) :: filter_snowc(ubc-lbc+1)
4493 integer , intent(inout) :: snl(1)
4494 real(kind_lake), intent(inout) :: dz(1,-nlevsnow+1:nlevsoil)
4495 real(kind_lake), intent(inout) :: zi(1,-nlevsnow+0:nlevsoil)
4496 real(kind_lake), intent(inout) :: t_soisno(1,-nlevsnow+1:nlevsoil)
4497 real(kind_lake), intent(inout) :: h2osoi_ice(1,-nlevsnow+1:nlevsoil)
4498 real(kind_lake), intent(inout) :: h2osoi_liq(1,-nlevsnow+1:nlevsoil)
4499
4500 !out:
4501
4502 real(kind_lake), intent(out) :: z(1,-nlevsnow+1:nlevsoil)
4503
4504
4505
4506 ! OTHER LOCAL VARIABLES:
4507
4508 integer :: j, c, fc ! indices
4509 real(kind_lake) :: drr ! thickness of the combined [m]
4510 integer :: msno ! number of snow layer 1 (top) to msno (bottom)
4511 real(kind_lake) :: dzsno(lbc:ubc,nlevsnow) ! Snow layer thickness [m]
4512 real(kind_lake) :: swice(lbc:ubc,nlevsnow) ! Partial volume of ice [m3/m3]
4513 real(kind_lake) :: swliq(lbc:ubc,nlevsnow) ! Partial volume of liquid water [m3/m3]
4514 real(kind_lake) :: tsno(lbc:ubc ,nlevsnow) ! Nodel temperature [K]
4515 real(kind_lake) :: zwice ! temporary
4516 real(kind_lake) :: zwliq ! temporary
4517 real(kind_lake) :: propor ! temporary
4518 !-----------------------------------------------------------------------
4519
4520 ! Begin calculation - note that the following column loops are only invoked
4521 ! for snow-covered columns
4522
4523 do j = 1,nlevsnow
4524 !dir$ concurrent
4525 !cdir nodep
4526 do fc = 1, num_snowc
4527 c = filter_snowc(fc)
4528 if (j <= abs(snl(c))) then
4529 dzsno(c,j) = dz(c,j+snl(c))
4530 swice(c,j) = h2osoi_ice(c,j+snl(c))
4531 swliq(c,j) = h2osoi_liq(c,j+snl(c))
4532 tsno(c,j) = t_soisno(c,j+snl(c))
4533 end if
4534 end do
4535 end do
4536
4537 !dir$ concurrent
4538 !cdir nodep
4539 do fc = 1, num_snowc
4540 c = filter_snowc(fc)
4541
4542 msno = abs(snl(c))
4543
4544 if (msno == 1) then
4545 ! Specify a new snow layer
4546 if (dzsno(c,1) > 0.03) then
4547 msno = 2
4548 dzsno(c,1) = dzsno(c,1)*0.5
4549 swice(c,1) = swice(c,1)*0.5
4550 swliq(c,1) = swliq(c,1)*0.5
4551 dzsno(c,2) = dzsno(c,1)
4552 swice(c,2) = swice(c,1)
4553 swliq(c,2) = swliq(c,1)
4554 tsno(c,2) = tsno(c,1)
4555 end if
4556 end if
4557
4558 if (msno > 1) then
4559 if (dzsno(c,1) > 0.02) then
4560 drr = dzsno(c,1) - 0.02
4561 propor = drr/dzsno(c,1)
4562 zwice = propor*swice(c,1)
4563 zwliq = propor*swliq(c,1)
4564 propor = 0.02/dzsno(c,1)
4565 swice(c,1) = propor*swice(c,1)
4566 swliq(c,1) = propor*swliq(c,1)
4567 dzsno(c,1) = 0.02
4568
4569 call combo (dzsno(c,2), swliq(c,2), swice(c,2), tsno(c,2), drr, &
4570 zwliq, zwice, tsno(c,1))
4571
4572 ! Subdivide a new layer
4573 if (msno <= 2 .and. dzsno(c,2) > 0.07) then
4574 msno = 3
4575 dzsno(c,2) = dzsno(c,2)*0.5
4576 swice(c,2) = swice(c,2)*0.5
4577 swliq(c,2) = swliq(c,2)*0.5
4578 dzsno(c,3) = dzsno(c,2)
4579 swice(c,3) = swice(c,2)
4580 swliq(c,3) = swliq(c,2)
4581 tsno(c,3) = tsno(c,2)
4582 end if
4583 end if
4584 end if
4585
4586 if (msno > 2) then
4587 if (dzsno(c,2) > 0.05) then
4588 drr = dzsno(c,2) - 0.05
4589 propor = drr/dzsno(c,2)
4590 zwice = propor*swice(c,2)
4591 zwliq = propor*swliq(c,2)
4592 propor = 0.05/dzsno(c,2)
4593 swice(c,2) = propor*swice(c,2)
4594 swliq(c,2) = propor*swliq(c,2)
4595 dzsno(c,2) = 0.05
4596
4597 call combo (dzsno(c,3), swliq(c,3), swice(c,3), tsno(c,3), drr, &
4598 zwliq, zwice, tsno(c,2))
4599
4600 ! Subdivided a new layer
4601 if (msno <= 3 .and. dzsno(c,3) > 0.18) then
4602 msno = 4
4603 dzsno(c,3) = dzsno(c,3)*0.5
4604 swice(c,3) = swice(c,3)*0.5
4605 swliq(c,3) = swliq(c,3)*0.5
4606 dzsno(c,4) = dzsno(c,3)
4607 swice(c,4) = swice(c,3)
4608 swliq(c,4) = swliq(c,3)
4609 tsno(c,4) = tsno(c,3)
4610 end if
4611 end if
4612 end if
4613
4614 if (msno > 3) then
4615 if (dzsno(c,3) > 0.11) then
4616 drr = dzsno(c,3) - 0.11
4617 propor = drr/dzsno(c,3)
4618 zwice = propor*swice(c,3)
4619 zwliq = propor*swliq(c,3)
4620 propor = 0.11/dzsno(c,3)
4621 swice(c,3) = propor*swice(c,3)
4622 swliq(c,3) = propor*swliq(c,3)
4623 dzsno(c,3) = 0.11
4624
4625 call combo (dzsno(c,4), swliq(c,4), swice(c,4), tsno(c,4), drr, &
4626 zwliq, zwice, tsno(c,3))
4627
4628 ! Subdivided a new layer
4629 if (msno <= 4 .and. dzsno(c,4) > 0.41) then
4630 msno = 5
4631 dzsno(c,4) = dzsno(c,4)*0.5
4632 swice(c,4) = swice(c,4)*0.5
4633 swliq(c,4) = swliq(c,4)*0.5
4634 dzsno(c,5) = dzsno(c,4)
4635 swice(c,5) = swice(c,4)
4636 swliq(c,5) = swliq(c,4)
4637 tsno(c,5) = tsno(c,4)
4638 end if
4639 end if
4640 end if
4641
4642 if (msno > 4) then
4643 if (dzsno(c,4) > 0.23) then
4644 drr = dzsno(c,4) - 0.23
4645 propor = drr/dzsno(c,4)
4646 zwice = propor*swice(c,4)
4647 zwliq = propor*swliq(c,4)
4648 propor = 0.23/dzsno(c,4)
4649 swice(c,4) = propor*swice(c,4)
4650 swliq(c,4) = propor*swliq(c,4)
4651 dzsno(c,4) = 0.23
4652
4653 call combo (dzsno(c,5), swliq(c,5), swice(c,5), tsno(c,5), drr, &
4654 zwliq, zwice, tsno(c,4))
4655 end if
4656 end if
4657
4658 snl(c) = -msno
4659
4660 end do
4661
4662 do j = -nlevsnow+1,0
4663 !dir$ concurrent
4664 !cdir nodep
4665 do fc = 1, num_snowc
4666 c = filter_snowc(fc)
4667 if (j >= snl(c)+1) then
4668 dz(c,j) = dzsno(c,j-snl(c))
4669 h2osoi_ice(c,j) = swice(c,j-snl(c))
4670 h2osoi_liq(c,j) = swliq(c,j-snl(c))
4671 t_soisno(c,j) = tsno(c,j-snl(c))
4672 end if
4673 end do
4674 end do
4675
4676 do j = 0, -nlevsnow+1, -1
4677 !dir$ concurrent
4678 !cdir nodep
4679 do fc = 1, num_snowc
4680 c = filter_snowc(fc)
4681 if (j >= snl(c)+1) then
4682 z(c,j) = zi(c,j) - 0.5*dz(c,j)
4683 zi(c,j-1) = zi(c,j) - dz(c,j)
4684 end if
4685 end do
4686 end do
4687
4688 end subroutine dividesnowlayers
4689
4692 subroutine combo(dz, wliq, wice, t, dz2, wliq2, wice2, t2)
4693 !
4694 ! DESCRIPTION:
4697 ! The combined temperature is based on the equation:
4698 ! the sum of the enthalpies of the two elements =
4699 ! that of the combined element.
4700 !
4701 ! !USES:
4702 !
4703 ! !ARGUMENTS:
4704 implicit none
4705 real(kind_lake), intent(in) :: dz2
4706 real(kind_lake), intent(in) :: wliq2
4707 real(kind_lake), intent(in) :: wice2
4708 real(kind_lake), intent(in) :: t2
4709 real(kind_lake), intent(inout) :: dz
4710 real(kind_lake), intent(inout) :: wliq
4711 real(kind_lake), intent(inout) :: wice
4712 real(kind_lake), intent(inout) :: t
4713 !
4714 ! !CALLED FROM:
4715 ! subroutine CombineSnowLayers in this module
4716 ! subroutine DivideSnowLayers in this module
4717 !
4718 ! !REVISION HISTORY:
4719 ! 15 September 1999: Yongjiu Dai; Initial code
4720 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
4721 ! June 2022: Sam Trahan; modified for CCPP
4722 !
4723 !EOP
4724 !
4725 ! !LOCAL VARIABLES:
4726 !
4727 real(kind_lake) :: dzc ! Total thickness of nodes 1 and 2 (dzc=dz+dz2).
4728 real(kind_lake) :: wliqc ! Combined liquid water [kg/m2]
4729 real(kind_lake) :: wicec ! Combined ice [kg/m2]
4730 real(kind_lake) :: tc ! Combined node temperature [K]
4731 real(kind_lake) :: h ! enthalpy of element 1 [J/m2]
4732 real(kind_lake) :: h2 ! enthalpy of element 2 [J/m2]
4733 real(kind_lake) :: hc ! temporary
4734 !-----------------------------------------------------------------------
4735
4736 dzc = dz+dz2
4737 wicec = (wice+wice2)
4738 wliqc = (wliq+wliq2)
4739 h = (cpice*wice+cpliq*wliq) * (t-tfrz)+hfus*wliq
4740 h2= (cpice*wice2+cpliq*wliq2) * (t2-tfrz)+hfus*wliq2
4741
4742 hc = h + h2
4743 if(hc < 0.)then
4744 tc = tfrz + hc/(cpice*wicec + cpliq*wliqc)
4745 else if (hc.le.hfus*wliqc) then
4746 tc = tfrz
4747 else
4748 tc = tfrz + (hc - hfus*wliqc) / (cpice*wicec + cpliq*wliqc)
4749 end if
4750
4751 dz = dzc
4752 wice = wicec
4753 wliq = wliqc
4754 t = tc
4755
4756 end subroutine combo
4757
4759 subroutine buildsnowfilter(lbc, ubc, num_nolakec, filter_nolakec,snl, & !i
4760 num_snowc, filter_snowc, & !o
4761 num_nosnowc, filter_nosnowc) !o
4762 !
4763 ! DESCRIPTION:
4765 !
4766 ! !USES:
4767 ! use clmtype
4768 !
4769 ! !ARGUMENTS:
4770 implicit none
4771 integer, intent(in) :: lbc, ubc
4772 integer, intent(in) :: num_nolakec
4773 integer, intent(in) :: filter_nolakec(ubc-lbc+1)
4774 integer, intent(in) :: snl(1)
4775 integer, intent(out) :: num_snowc
4776 integer, intent(out) :: filter_snowc(ubc-lbc+1)
4777 integer, intent(out) :: num_nosnowc
4778 integer, intent(out) :: filter_nosnowc(ubc-lbc+1)
4779 !
4780 ! !CALLED FROM:
4781 ! subroutine Hydrology2 in Hydrology2Mod
4782 ! subroutine CombineSnowLayers in this module
4783 !
4784 ! !REVISION HISTORY:
4785 ! 2003 July 31: Forrest Hoffman
4786 ! 2022 June: Sam Trahan modified for CCPP
4787 !
4788 ! !LOCAL VARIABLES:
4789 ! local pointers to implicit in arguments
4790 !
4791 !EOP
4792 !
4793 ! !OTHER LOCAL VARIABLES:
4794 integer :: fc, c
4795 !-----------------------------------------------------------------------
4796
4797
4798 ! Build snow/no-snow filters for other subroutines
4799
4800 num_snowc = 0
4801 num_nosnowc = 0
4802 do fc = 1, num_nolakec
4803 c = filter_nolakec(fc)
4804 if (snl(c) < 0) then
4805 num_snowc = num_snowc + 1
4806 filter_snowc(num_snowc) = c
4807 else
4808 num_nosnowc = num_nosnowc + 1
4809 filter_nosnowc(num_nosnowc) = c
4810 end if
4811 end do
4812
4813 end subroutine buildsnowfilter
4814
4815
4816 ! DESCRIPTION:
4823subroutine frictionvelocity(pgridcell,forc_hgt,forc_hgt_u, & !i
4824 forc_hgt_t,forc_hgt_q, & !i
4825 lbp, ubp, fn, filterp, & !i
4826 displa, z0m, z0h, z0q, & !i
4827 obu, iter, ur, um, & !i
4828 ustar,temp1, temp2, temp12m, temp22m, & !o
4829 u10,fv, & !o
4830 fm) !i&o
4831
4832 !=============================================================================
4833 ! REVISION HISTORY:
4834 ! 15 September 1999: Yongjiu Dai; Initial code
4835 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
4836 ! 12/19/01, Peter Thornton
4837 ! Added arguments to eliminate passing clm derived type into this function.
4838 ! Created by Mariana Vertenstein
4839 ! June 2022: Sam Trahan modified for CCPP
4840 !============================================================================
4841 ! !USES:
4842 ! use clmtype
4843 !!use clm_atmlnd, only : clm_a2l
4844 !
4845 ! !ARGUMENTS:
4846 implicit none
4847
4848 !in:
4849
4850 integer , intent(in) :: pgridcell(1)
4851 real(kind_lake), intent(in) :: forc_hgt(1)
4852 real(kind_lake), intent(in) :: forc_hgt_u(1)
4853 real(kind_lake), intent(in) :: forc_hgt_t(1)
4854 real(kind_lake), intent(in) :: forc_hgt_q(1)
4855 integer , intent(in) :: lbp, ubp
4856 integer , intent(in) :: fn
4857 integer , intent(in) :: filterp(fn)
4858 real(kind_lake), intent(in) :: displa(lbp:ubp)
4859 real(kind_lake), intent(in) :: z0m(lbp:ubp)
4860 real(kind_lake), intent(in) :: z0h(lbp:ubp)
4861 real(kind_lake), intent(in) :: z0q(lbp:ubp)
4862 real(kind_lake), intent(in) :: obu(lbp:ubp)
4863 integer, intent(in) :: iter
4864 real(kind_lake), intent(in) :: ur(lbp:ubp)
4865 real(kind_lake), intent(in) :: um(lbp:ubp)
4866
4867 !out:
4868
4869 real(kind_lake), intent(out) :: ustar(lbp:ubp)
4870 real(kind_lake), intent(out) :: temp1(lbp:ubp)
4871 real(kind_lake), intent(out) :: temp12m(lbp:ubp)
4872 real(kind_lake), intent(out) :: temp2(lbp:ubp)
4873 real(kind_lake), intent(out) :: temp22m(lbp:ubp)
4874 real(kind_lake), intent(out) :: u10(1)
4875 real(kind_lake), intent(out) :: fv(1)
4876
4877 !inout:
4878 real(kind_lake), intent(inout) :: fm(lbp:ubp)
4879
4880 ! OTHER LOCAL VARIABLES:
4881
4882 real(kind_lake), parameter :: zetam = 1.574_kind_lake ! transition point of flux-gradient relation (wind profile)
4883 real(kind_lake), parameter :: zetat = 0.465_kind_lake ! transition point of flux-gradient relation (temp. profile)
4884 integer :: f ! pft-filter index
4885 integer :: p ! pft index
4886 integer :: g ! gridcell index
4887 real(kind_lake):: zldis(lbp:ubp) ! reference height "minus" zero displacement heght [m]
4888 real(kind_lake):: zeta(lbp:ubp) ! dimensionless height used in Monin-Obukhov theory
4889
4890 !------------------------------------------------------------------------------
4891
4892
4893 ! Adjustment factors for unstable (moz < 0) or stable (moz > 0) conditions.
4894
4895 if_not_pergro: if(.not.pergro) then
4896
4897 !dir$ concurrent
4898 !cdir nodep
4899 do f = 1, fn
4900 p = filterp(f)
4901 g = pgridcell(p)
4902
4903 ! Wind profile
4904
4905 zldis(p) = forc_hgt_u(g)-displa(p)
4906 zeta(p) = zldis(p)/obu(p)
4907 if (zeta(p) < -zetam) then
4908 ustar(p) = vkc*um(p)/(log(-zetam*obu(p)/z0m(p))&
4909 - stabilityfunc1(-zetam) &
4910 + stabilityfunc1(z0m(p)/obu(p)) &
4911 + 1.14_kind_lake*((-zeta(p))**0.333_kind_lake-(zetam)**0.333_kind_lake))
4912 else if (zeta(p) < 0._kind_lake) then
4913 ustar(p) = vkc*um(p)/(log(zldis(p)/z0m(p))&
4914 - stabilityfunc1(zeta(p))&
4915 + stabilityfunc1(z0m(p)/obu(p)))
4916 else if (zeta(p) <= 1._kind_lake) then
4917 ustar(p) = vkc*um(p)/(log(zldis(p)/z0m(p)) + 5._kind_lake*zeta(p) -5._kind_lake*z0m(p)/obu(p))
4918 else
4919 ustar(p) = vkc*um(p)/(log(obu(p)/z0m(p))+5._kind_lake-5._kind_lake*z0m(p)/obu(p) &
4920 +(5._kind_lake*log(zeta(p))+zeta(p)-1._kind_lake))
4921 end if
4922
4923 ! Temperature profile
4924
4925 zldis(p) = forc_hgt_t(g)-displa(p)
4926 zeta(p) = zldis(p)/obu(p)
4927 if (zeta(p) < -zetat) then
4928 temp1(p) = vkc/(log(-zetat*obu(p)/z0h(p))&
4929 - stabilityfunc2(-zetat) &
4930 + stabilityfunc2(z0h(p)/obu(p)) &
4931 + 0.8_kind_lake*((zetat)**(-0.333_kind_lake)-(-zeta(p))**(-0.333_kind_lake)))
4932 else if (zeta(p) < 0._kind_lake) then
4933 temp1(p) = vkc/(log(zldis(p)/z0h(p)) &
4934 - stabilityfunc2(zeta(p)) &
4935 + stabilityfunc2(z0h(p)/obu(p)))
4936 else if (zeta(p) <= 1._kind_lake) then
4937 temp1(p) = vkc/(log(zldis(p)/z0h(p)) + 5._kind_lake*zeta(p) - 5._kind_lake*z0h(p)/obu(p))
4938 else
4939 temp1(p) = vkc/(log(obu(p)/z0h(p)) + 5._kind_lake - 5._kind_lake*z0h(p)/obu(p) &
4940 + (5._kind_lake*log(zeta(p))+zeta(p)-1._kind_lake))
4941 end if
4942
4943 ! Humidity profile
4944
4945 if (forc_hgt_q(g) == forc_hgt_t(g) .and. z0q(p) == z0h(p)) then
4946 temp2(p) = temp1(p)
4947 else
4948 zldis(p) = forc_hgt_q(g)-displa(p)
4949 zeta(p) = zldis(p)/obu(p)
4950 if (zeta(p) < -zetat) then
4951 temp2(p) = vkc/(log(-zetat*obu(p)/z0q(p)) &
4952 - stabilityfunc2(-zetat) &
4953 + stabilityfunc2(z0q(p)/obu(p)) &
4954 + 0.8_kind_lake*((zetat)**(-0.333_kind_lake)-(-zeta(p))**(-0.333_kind_lake)))
4955 else if (zeta(p) < 0._kind_lake) then
4956 temp2(p) = vkc/(log(zldis(p)/z0q(p)) &
4957 - stabilityfunc2(zeta(p)) &
4958 + stabilityfunc2(z0q(p)/obu(p)))
4959 else if (zeta(p) <= 1._kind_lake) then
4960 temp2(p) = vkc/(log(zldis(p)/z0q(p)) + 5._kind_lake*zeta(p)-5._kind_lake*z0q(p)/obu(p))
4961 else
4962 temp2(p) = vkc/(log(obu(p)/z0q(p)) + 5._kind_lake - 5._kind_lake*z0q(p)/obu(p) &
4963 + (5._kind_lake*log(zeta(p))+zeta(p)-1._kind_lake))
4964 end if
4965 endif
4966
4967 ! Temperature profile applied at 2-m
4968
4969 zldis(p) = 2.0_kind_lake + z0h(p)
4970 zeta(p) = zldis(p)/obu(p)
4971 if (zeta(p) < -zetat) then
4972 temp12m(p) = vkc/(log(-zetat*obu(p)/z0h(p))&
4973 - stabilityfunc2(-zetat) &
4974 + stabilityfunc2(z0h(p)/obu(p)) &
4975 + 0.8_kind_lake*((zetat)**(-0.333_kind_lake)-(-zeta(p))**(-0.333_kind_lake)))
4976 else if (zeta(p) < 0._kind_lake) then
4977 temp12m(p) = vkc/(log(zldis(p)/z0h(p)) &
4978 - stabilityfunc2(zeta(p)) &
4979 + stabilityfunc2(z0h(p)/obu(p)))
4980 else if (zeta(p) <= 1._kind_lake) then
4981 temp12m(p) = vkc/(log(zldis(p)/z0h(p)) + 5._kind_lake*zeta(p) - 5._kind_lake*z0h(p)/obu(p))
4982 else
4983 temp12m(p) = vkc/(log(obu(p)/z0h(p)) + 5._kind_lake - 5._kind_lake*z0h(p)/obu(p) &
4984 + (5._kind_lake*log(zeta(p))+zeta(p)-1._kind_lake))
4985 end if
4986
4987 ! Humidity profile applied at 2-m
4988
4989 if (z0q(p) == z0h(p)) then
4990 temp22m(p) = temp12m(p)
4991 else
4992 zldis(p) = 2.0_kind_lake + z0q(p)
4993 zeta(p) = zldis(p)/obu(p)
4994 if (zeta(p) < -zetat) then
4995 temp22m(p) = vkc/(log(-zetat*obu(p)/z0q(p)) - &
4996 stabilityfunc2(-zetat) + stabilityfunc2(z0q(p)/obu(p)) &
4997 + 0.8_kind_lake*((zetat)**(-0.333_kind_lake)-(-zeta(p))**(-0.333_kind_lake)))
4998 else if (zeta(p) < 0._kind_lake) then
4999 temp22m(p) = vkc/(log(zldis(p)/z0q(p)) - &
5000 stabilityfunc2(zeta(p))+stabilityfunc2(z0q(p)/obu(p)))
5001 else if (zeta(p) <= 1._kind_lake) then
5002 temp22m(p) = vkc/(log(zldis(p)/z0q(p)) + 5._kind_lake*zeta(p)-5._kind_lake*z0q(p)/obu(p))
5003 else
5004 temp22m(p) = vkc/(log(obu(p)/z0q(p)) + 5._kind_lake - 5._kind_lake*z0q(p)/obu(p) &
5005 + (5._kind_lake*log(zeta(p))+zeta(p)-1._kind_lake))
5006 end if
5007 end if
5008 end do
5009 endif if_not_pergro
5010
5011
5012if_pergro: if (pergro) then
5013
5014 !===============================================================================
5015 ! The following only applies when PERGRO is defined
5016 !===============================================================================
5017
5018 !dir$ concurrent
5019 !cdir nodep
5020 do f = 1, fn
5021 p = filterp(f)
5022 g = pgridcell(p)
5023
5024 zldis(p) = forc_hgt_u(g)-displa(p)
5025 zeta(p) = zldis(p)/obu(p)
5026 if (zeta(p) < -zetam) then ! zeta < -1
5027 ustar(p) = vkc * um(p) / log(-zetam*obu(p)/z0m(p))
5028 else if (zeta(p) < 0._kind_lake) then ! -1 <= zeta < 0
5029 ustar(p) = vkc * um(p) / log(zldis(p)/z0m(p))
5030 else if (zeta(p) <= 1._kind_lake) then ! 0 <= ztea <= 1
5031 ustar(p)=vkc * um(p)/log(zldis(p)/z0m(p))
5032 else ! 1 < zeta, phi=5+zeta
5033 ustar(p)=vkc * um(p)/log(obu(p)/z0m(p))
5034 endif
5035
5036 zldis(p) = forc_hgt_t(g)-displa(p)
5037 zeta(p) = zldis(p)/obu(p)
5038 if (zeta(p) < -zetat) then
5039 temp1(p)=vkc/log(-zetat*obu(p)/z0h(p))
5040 else if (zeta(p) < 0._kind_lake) then
5041 temp1(p)=vkc/log(zldis(p)/z0h(p))
5042 else if (zeta(p) <= 1._kind_lake) then
5043 temp1(p)=vkc/log(zldis(p)/z0h(p))
5044 else
5045 temp1(p)=vkc/log(obu(p)/z0h(p))
5046 end if
5047
5048 zldis(p) = forc_hgt_q(g)-displa(p)
5049 zeta(p) = zldis(p)/obu(p)
5050 if (zeta(p) < -zetat) then
5051 temp2(p)=vkc/log(-zetat*obu(p)/z0q(p))
5052 else if (zeta(p) < 0._kind_lake) then
5053 temp2(p)=vkc/log(zldis(p)/z0q(p))
5054 else if (zeta(p) <= 1._kind_lake) then
5055 temp2(p)=vkc/log(zldis(p)/z0q(p))
5056 else
5057 temp2(p)=vkc/log(obu(p)/z0q(p))
5058 end if
5059
5060 zldis(p) = 2.0_kind_lake + z0h(p)
5061 zeta(p) = zldis(p)/obu(p)
5062 if (zeta(p) < -zetat) then
5063 temp12m(p)=vkc/log(-zetat*obu(p)/z0h(p))
5064 else if (zeta(p) < 0._kind_lake) then
5065 temp12m(p)=vkc/log(zldis(p)/z0h(p))
5066 else if (zeta(p) <= 1._kind_lake) then
5067 temp12m(p)=vkc/log(zldis(p)/z0h(p))
5068 else
5069 temp12m(p)=vkc/log(obu(p)/z0h(p))
5070 end if
5071
5072 zldis(p) = 2.0_kind_lake + z0q(p)
5073 zeta(p) = zldis(p)/obu(p)
5074 if (zeta(p) < -zetat) then
5075 temp22m(p)=vkc/log(-zetat*obu(p)/z0q(p))
5076 else if (zeta(p) < 0._kind_lake) then
5077 temp22m(p)=vkc/log(zldis(p)/z0q(p))
5078 else if (zeta(p) <= 1._kind_lake) then
5079 temp22m(p)=vkc/log(zldis(p)/z0q(p))
5080 else
5081 temp22m(p)=vkc/log(obu(p)/z0q(p))
5082 end if
5083 end do
5084
5085 endif if_pergro
5086
5087 end subroutine frictionvelocity
5088
5089 ! !IROUTINE: StabilityFunc
5090 !
5091 ! !INTERFACE:
5092 real(kind_lake) function stabilityfunc1(zeta)
5093 !
5094 ! DESCRIPTION:
5096 !
5097 ! !USES:
5098 ! use shr_const_mod, only: SHR_CONST_PI
5099 !Zack Subin, 7/8/08
5100 !
5101 ! !ARGUMENTS:
5102 implicit none
5103 real(kind_lake), intent(in) :: zeta
5104 !
5105 ! !CALLED FROM:
5106 ! subroutine FrictionVelocity in this module
5107 !
5108 ! !REVISION HISTORY:
5109 ! 15 September 1999: Yongjiu Dai; Initial code
5110 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
5111 ! June 2022: Sam Trahan; modified for CCPP
5112 !
5113 !EOP
5114 !
5115 ! !LOCAL VARIABLES:
5116 real(kind_lake) :: chik, chik2
5117 !------------------------------------------------------------------------------
5118
5119 chik2 = sqrt(1._kind_lake-16._kind_lake*zeta)
5120 chik = sqrt(chik2)
5121 stabilityfunc1 = 2._kind_lake*log((1._kind_lake+chik)*0.5_kind_lake) &
5122 !Changed to pie, Zack Subin, 7/9/08
5123 !Spelling corrected, changed to pi, Sam Trahan the Killjoy, 6/2/22
5124 + log((1._kind_lake+chik2)*0.5_kind_lake)-2._kind_lake*atan(chik)+pi*0.5_kind_lake
5125
5126 end function stabilityfunc1
5127
5128 !------------------------------------------------------------------------------
5129 !BOP
5130 !
5131 ! !IROUTINE: StabilityFunc2
5132 !
5133 ! !INTERFACE:
5134 real(kind_lake) function stabilityfunc2(zeta)
5135 !
5136 ! !DESCRIPTION:
5138 !
5139 ! !USES:
5140 !Removed by Zack Subin, 7/9/08
5141 ! use shr_const_mod, only: SHR_CONST_PI
5142 !
5143 ! !ARGUMENTS:
5144 implicit none
5145 real(kind_lake), intent(in) :: zeta
5146 !
5147 ! !CALLED FROM:
5148 ! subroutine FrictionVelocity in this module
5149 !
5150 ! !REVISION HISTORY:
5151 ! 15 September 1999: Yongjiu Dai; Initial code
5152 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
5153 ! June 2022: Sam Trahan modified for CCPP
5154 !
5155 !EOP
5156 !
5157 ! !LOCAL VARIABLES:
5158 real(kind_lake) :: chik2
5159 !------------------------------------------------------------------------------
5160
5161 chik2 = sqrt(1._kind_lake-16._kind_lake*zeta)
5162 stabilityfunc2 = 2._kind_lake*log((1._kind_lake+chik2)*0.5_kind_lake)
5163
5164 end function stabilityfunc2
5165
5166 !-----------------------------------------------------------------------
5167 !BOP
5168 !
5169 ! !IROUTINE: MoninObukIni
5170 !
5171 ! !INTERFACE:
5172 subroutine moninobukini (ur, thv, dthv, zldis, z0m, um, obu)
5173 !
5174 ! !DESCRIPTION:
5180 !
5181 ! !USES:
5182 !
5183 ! !ARGUMENTS:
5184 implicit none
5185 real(kind_lake), intent(in) :: ur
5186 real(kind_lake), intent(in) :: thv
5187 real(kind_lake), intent(in) :: dthv
5188 real(kind_lake), intent(in) :: zldis
5189 real(kind_lake), intent(in) :: z0m
5190 real(kind_lake), intent(out) :: um
5191 real(kind_lake), intent(out) :: obu
5192 !
5193 ! !CALLED FROM:
5194 ! subroutine BareGroundFluxes in module BareGroundFluxesMod.F90
5195 ! subroutine BiogeophysicsLake in module BiogeophysicsLakeMod.F90
5196 ! subroutine CanopyFluxes in module CanopyFluxesMod.F90
5197 !
5198 ! !REVISION HISTORY:
5199 ! 15 September 1999: Yongjiu Dai; Initial code
5200 ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
5201 ! June 2022: Sam Trahan modified for CCPP
5202 !
5203 !EOP
5204 !
5205 ! !LOCAL VARIABLES:
5206 !
5207 real(kind_lake) :: wc ! convective velocity [m/s]
5208 real(kind_lake) :: rib ! bulk Richardson number
5209 real(kind_lake) :: zeta ! dimensionless height used in Monin-Obukhov theory
5210 real(kind_lake) :: ustar ! friction velocity [m/s]
5211 !-----------------------------------------------------------------------
5212
5213 ! Initial values of u* and convective velocity
5214
5215 ustar=0.06_kind_lake
5216 wc=0.5_kind_lake
5217 if (dthv >= 0._kind_lake) then
5218 um=max(ur,0.1_kind_lake)
5219 else
5220 um=sqrt(ur*ur+wc*wc)
5221 endif
5222
5223 rib=grav*zldis*dthv/(thv*um*um)
5224 if (pergro) then
5225 rib = 0._kind_lake
5226 endif
5227
5228 if (rib >= 0._kind_lake) then ! neutral or stable
5229 zeta = rib*log(zldis/z0m)/(1._kind_lake-5._kind_lake*min(rib,0.19_kind_lake))
5230 zeta = min(2._kind_lake,max(zeta,0.01_kind_lake ))
5231 else ! unstable
5232 zeta=rib*log(zldis/z0m)
5233 zeta = max(-100._kind_lake,min(zeta,-0.01_kind_lake ))
5234 endif
5235
5236 obu=zldis/zeta
5237
5238 end subroutine moninobukini
5239
5243 subroutine clm_lake_init(con_pi,karman,con_g,con_sbc,con_t0c,rhowater,con_csol,con_cliq, &
5244 con_hfus,con_hvap,con_rd,con_cp,rholakeice,clm_lake_debug, &
5245 clm_debug_print,con_eps_model,con_fvirt_model,errmsg,errflg)
5246 implicit none
5247 real(kind_phys), intent(in) :: con_pi,karman,con_g,con_sbc,con_t0c, &
5248 rhowater,con_csol,con_cliq, con_hfus,con_hvap,con_rd,con_cp, &
5249 rholakeice,con_eps_model,con_fvirt_model
5250 INTEGER, INTENT(OUT) :: errflg
5251 CHARACTER(*), INTENT(OUT) :: errmsg
5252 logical, intent(in) :: clm_lake_debug,clm_debug_print
5253 integer :: i, j
5254
5255 lakedebug = clm_lake_debug
5256 debug_print = clm_debug_print
5257 if(debug_print) then
5258 write(0,*) 'clm_lake_init'
5259 endif
5260
5261 errflg=0
5262 errmsg=''
5263
5264 pi = con_pi
5265 vkc = karman
5266 grav = con_g
5267 sb = con_sbc
5268 tfrz = con_t0c
5269 denh2o = rhowater
5270 denice = rholakeice
5271 cpice = con_csol
5272 cpliq = con_cliq
5273 hfus = con_hfus
5274 hvap = con_hvap
5275 hsub = con_hfus+con_hvap
5276 invhvap = 1._kind_lake/hvap
5277 invhsub = 1._kind_lake/hsub
5278 rair = con_rd
5279 cpair = con_cp
5280 con_eps = con_eps_model
5281 con_fvirt = con_fvirt_model
5282 one_minus_con_eps = 1.0_kind_lake - con_eps
5283
5284 ! dzlak(1) = 0.1_kind_lake
5285 ! dzlak(2) = 1._kind_lake
5286 ! dzlak(3) = 2._kind_lake
5287 ! dzlak(4) = 3._kind_lake
5288 ! dzlak(5) = 4._kind_lake
5289 ! dzlak(6) = 5._kind_lake
5290 ! dzlak(7) = 7._kind_lake
5291 ! dzlak(8) = 7._kind_lake
5292 ! dzlak(9) = 10.45_kind_lake
5293 ! dzlak(10)= 10.45_kind_lake
5294 !
5295 ! zlak(1) = 0.05_kind_lake
5296 ! zlak(2) = 0.6_kind_lake
5297 ! zlak(3) = 2.1_kind_lake
5298 ! zlak(4) = 4.6_kind_lake
5299 ! zlak(5) = 8.1_kind_lake
5300 ! zlak(6) = 12.6_kind_lake
5301 ! zlak(7) = 18.6_kind_lake
5302 ! zlak(8) = 25.6_kind_lake
5303 ! zlak(9) = 34.325_kind_lake
5304 ! zlak(10)= 44.775_kind_lake
5305 dzlak(1) = 0.1_kind_lake
5306 dzlak(2) = 0.1_kind_lake
5307 dzlak(3) = 0.1_kind_lake
5308 dzlak(4) = 0.1_kind_lake
5309 dzlak(5) = 0.1_kind_lake
5310 dzlak(6) = 0.1_kind_lake
5311 dzlak(7) = 0.1_kind_lake
5312 dzlak(8) = 0.1_kind_lake
5313 dzlak(9) = 0.1_kind_lake
5314 dzlak(10)= 0.1_kind_lake
5315
5316 zlak(1) = 0.05_kind_lake
5317 zlak(2) = 0.15_kind_lake
5318 zlak(3) = 0.25_kind_lake
5319 zlak(4) = 0.35_kind_lake
5320 zlak(5) = 0.45_kind_lake
5321 zlak(6) = 0.55_kind_lake
5322 zlak(7) = 0.65_kind_lake
5323 zlak(8) = 0.75_kind_lake
5324 zlak(9) = 0.85_kind_lake
5325 zlak(10)= 0.95_kind_lake
5326
5327 ! "0" refers to soil surface and "nlevsoil" refers to the bottom of model soil
5328
5329 do j = 1, nlevsoil
5330 zsoi(j) = scalez*(exp(0.5_kind_lake*(j-0.5_kind_lake))-1._kind_lake) !node depths
5331 enddo
5332
5333 dzsoi(1) = 0.5_kind_lake*(zsoi(1)+zsoi(2)) !thickness b/n two interfaces
5334 do j = 2,nlevsoil-1
5335 dzsoi(j)= 0.5_kind_lake*(zsoi(j+1)-zsoi(j-1))
5336 enddo
5337 dzsoi(nlevsoil) = zsoi(nlevsoil)-zsoi(nlevsoil-1)
5338
5339 zisoi(0) = 0._kind_lake
5340 do j = 1, nlevsoil-1
5341 zisoi(j) = 0.5_kind_lake*(zsoi(j)+zsoi(j+1)) !interface depths
5342 enddo
5343 zisoi(nlevsoil) = zsoi(nlevsoil) + 0.5_kind_lake*dzsoi(nlevsoil)
5344
5345 end subroutine clm_lake_init
5346
5347 SUBROUTINE lakeini(kdt, ISLTYP, gt0, snowd, & !i
5348 weasd, lakedepth_default, fhour, &
5349 oro_lakedepth, savedtke12d, snowdp2d, h2osno2d, & !o
5350 snl2d, t_grnd2d, t_lake3d, lake_icefrac3d, &
5351 t_soisno3d, h2osoi_ice3d, &
5352 h2osoi_liq3d, h2osoi_vol3d, z3d, dz3d, &
5353 zi3d, &
5354 fice, hice, min_lakeice, tsfc, &
5355 use_lake_model, use_lakedepth, &
5356 im, prsi, &
5357 xlat_d, xlon_d, clm_lake_initialized, &
5358 input_lakedepth, tg3, clm_lakedepth, &
5359 km, me, master, errmsg, errflg)
5360
5367
5368 !==============================================================================
5369 ! This subroutine was first edited by Hongping Gu for coupling
5370 ! 07/20/2010
5371 ! Long after, in June 2022, Sam Trahan updated it for CCPP
5372 !==============================================================================
5373
5374 implicit none
5375
5376 INTEGER, INTENT(OUT) :: errflg
5377 CHARACTER(*), INTENT(OUT) :: errmsg
5378
5379 INTEGER , INTENT (IN) :: im, me, master, km, kdt
5380 REAL(KIND_PHYS), INTENT(IN) :: min_lakeice, fhour
5381 REAL(KIND_PHYS), DIMENSION(IM), INTENT(INOUT):: FICE, hice
5382 REAL(KIND_PHYS), DIMENSION(IM), INTENT(IN):: TG3, xlat_d, xlon_d
5383 REAL(KIND_PHYS), DIMENSION(IM), INTENT(IN):: tsfc
5384 REAL(KIND_PHYS), DIMENSION(IM) ,INTENT(INOUT) :: clm_lake_initialized
5385 integer, dimension(IM), intent(in) :: use_lake_model
5386 !INTEGER , INTENT (IN) :: lakeflag
5387 !INTEGER , INTENT (INOUT) :: lake_depth_flag
5388 LOGICAL, INTENT (IN) :: use_lakedepth
5389
5390 INTEGER, DIMENSION(IM), INTENT(IN) :: ISLTYP
5391 REAL(KIND_PHYS), DIMENSION(IM), INTENT(INOUT) :: snowd,weasd
5392 REAL(kind_phys), DIMENSION(IM,KM), INTENT(IN) :: gt0, prsi
5393 real(kind_phys), intent(in) :: lakedepth_default
5394
5395 real(kind_phys), dimension(IM),intent(inout) :: clm_lakedepth
5396 real(kind_phys), dimension(IM),intent(inout) :: input_lakedepth
5397 real(kind_phys), dimension(IM),intent(in) :: oro_lakedepth
5398 real(kind_phys), dimension(IM),intent(out) :: savedtke12d
5399 real(kind_phys), dimension(IM),intent(out) :: snowdp2d, &
5400 h2osno2d, &
5401 snl2d, &
5402 t_grnd2d
5403
5404 real(kind_phys), dimension(IM,nlevlake),INTENT(out) :: t_lake3d, &
5405 lake_icefrac3d
5406 real(kind_phys), dimension(IM,-nlevsnow+1:nlevsoil ),INTENT(out) :: t_soisno3d, &
5407 h2osoi_ice3d, &
5408 h2osoi_liq3d, &
5409 h2osoi_vol3d, &
5410 z3d, &
5411 dz3d
5412
5413 real(kind_phys), dimension( IM,-nlevsnow+0:nlevsoil ),INTENT(out) :: zi3d
5414
5415 !LOGICAL, DIMENSION( : ),intent(out) :: lake
5416 !REAL(KIND_PHYS), OPTIONAL, DIMENSION( : ), INTENT(IN) :: lake_depth ! no separate variable for this in CCPP
5417
5418 integer :: n,i,j,k,ib,lev,bottom ! indices
5419 real(kind_lake),dimension(1:im ) :: bd2d ! bulk density of dry soil material [kg/m^3]
5420 real(kind_lake),dimension(1:im ) :: tkm2d ! mineral conductivity
5421 real(kind_lake),dimension(1:im ) :: xksat2d ! maximum hydraulic conductivity of soil [mm/s]
5422 real(kind_lake),dimension(1:im ) :: depthratio2d ! ratio of lake depth to standard deep lake depth
5423
5424 logical,parameter :: arbinit = .false.
5425 real(kind_lake),parameter :: defval = -999.0
5426 integer :: isl
5427 integer :: numb_lak ! for debug
5428 character*256 :: message
5429 real(kind_lake) :: ht
5430 real(kind_lake) :: rhosn
5431 real(kind_lake) :: depth, lakedepth
5432
5433 logical :: climatology_limits
5434
5435 real(kind_lake) :: z_lake(nlevlake) ! layer depth for lake (m)
5436 real(kind_lake) :: dz_lake(nlevlake) ! layer thickness for lake (m)
5437
5438 integer, parameter :: xcheck=38
5439 integer, parameter :: ycheck=92
5440
5441 integer :: used_lakedepth_default, init_points, month, julday
5442 integer :: mon, iday, num2, num1, juld, day2, day1, wght1, wght2
5443 real(kind_lake) :: Tclim, watsat
5444
5445 used_lakedepth_default=0
5446
5447 errmsg = ''
5448 errflg = 0
5449
5450 !!!!!!!!!!!!!!!!!!begin to initialize lake variables!!!!!!!!!!!!!!!!!!
5451
5452 init_points=0
5453 do_init: DO i=1,im
5454 if(use_lake_model(i)==0 .or. clm_lake_initialized(i)>0) then
5455 cycle
5456 endif
5457
5458 ! To handle cold-start with bad lakedepth2d
5459 if ( use_lakedepth ) then
5460 if (oro_lakedepth(i) == 10.0 .or. oro_lakedepth(i) <= 0.) then
5461 !- 10.0 is the fill value for lake depth, in this case set to default value
5462 clm_lakedepth(i) = lakedepth_default
5463 used_lakedepth_default = used_lakedepth_default+1
5464 else
5465 clm_lakedepth(i) = oro_lakedepth(i)
5466 endif
5467 else
5468 !- all lakes are initialized with the default lake depth
5469 clm_lakedepth(i) = lakedepth_default
5470 used_lakedepth_default = used_lakedepth_default+1
5471 endif
5472
5473 if(clm_lake_initialized(i)>0) then
5474 cycle
5475 endif
5476
5477 input_lakedepth=clm_lakedepth
5478
5479 snl2d(i) = defval
5480 do k = -nlevsnow+1,nlevsoil
5481 h2osoi_liq3d(i,k) = defval
5482 h2osoi_ice3d(i,k) = defval
5483 h2osoi_vol3d(i,k) = defval
5484 t_soisno3d(i,k) = defval
5485 z3d(i,k) = defval
5486 dz3d(i,k) = defval
5487 enddo
5488 do k = 1,nlevlake
5489 t_lake3d(i,k) = defval
5490 lake_icefrac3d(i,k) = defval
5491 enddo
5492
5493 if (use_lake_model(i) == 1) then
5494 ! for lake points only
5495 z3d(i,:) = 0.0
5496 dz3d(i,:) = 0.0
5497 zi3d(i,:) = 0.0
5498 h2osoi_liq3d(i,:) = 0.0
5499 h2osoi_ice3d(i,:) = 0.0
5500 lake_icefrac3d(i,:) = 0.0
5501 h2osoi_vol3d(i,:) = 0.0
5502 snl2d(i) = 0.0
5503
5504 if(fice(i)>min_lakeice) then
5505 lake_icefrac3d(i,1) = fice(i)
5506 snowdp2d(i) = snowd(i)*1e-3 ! SNOW in kg/m^2 and snowdp in m
5507 h2osno2d(i) = weasd(i) ! mm
5508 else
5509 fice(i) = 0.
5510 snowd(i) = 0.
5511 weasd(i) = 0.
5512 snowdp2d(i) = 0.
5513 h2osno2d(i) = 0.
5514 endif
5515
5516 ! Soil hydraulic and thermal properties
5517 isl = isltyp(i)
5518 if (isl == 0 ) isl = 14
5519 if (isl == 14 ) isl = isl + 1
5520
5521 call calculate_z_dz_lake(i,input_lakedepth,clm_lakedepth,z_lake,dz_lake)
5522
5523 z3d(i,1:nlevsoil) = zsoi(1:nlevsoil)
5524 zi3d(i,0:nlevsoil) = zisoi(0:nlevsoil)
5525 dz3d(i,1:nlevsoil) = dzsoi(1:nlevsoil)
5526 savedtke12d(i) = tkwat ! Initialize for first timestep.
5527
5528
5529 if (snowdp2d(i) < 0.01_kind_lake) then
5530 snl2d(i) = 0
5531 dz3d(i,-nlevsnow+1:0) = 0._kind_lake
5532 z3d(i,-nlevsnow+1:0) = 0._kind_lake
5533 zi3d(i,-nlevsnow+0:0) = 0._kind_lake
5534 else
5535 if ((snowdp2d(i) >= 0.01_kind_lake) .and. (snowdp2d(i) <= 0.03_kind_lake)) then
5536 snl2d(i) = -1
5537 dz3d(i,0) = snowdp2d(i)
5538 else if ((snowdp2d(i) > 0.03_kind_lake) .and. (snowdp2d(i) <= 0.04_kind_lake)) then
5539 snl2d(i) = -2
5540 dz3d(i,-1) = snowdp2d(i)*0.5_kind_lake
5541 dz3d(i, 0) = dz3d(i,-1)
5542 else if ((snowdp2d(i) > 0.04_kind_lake) .and. (snowdp2d(i) <= 0.07_kind_lake)) then
5543 snl2d(i) = -2
5544 dz3d(i,-1) = 0.02_kind_lake
5545 dz3d(i, 0) = snowdp2d(i) - dz3d(i,-1)
5546 else if ((snowdp2d(i) > 0.07_kind_lake) .and. (snowdp2d(i) <= 0.12_kind_lake)) then
5547 snl2d(i) = -3
5548 dz3d(i,-2) = 0.02_kind_lake
5549 dz3d(i,-1) = (snowdp2d(i) - 0.02_kind_lake)*0.5_kind_lake
5550 dz3d(i, 0) = dz3d(i,-1)
5551 else if ((snowdp2d(i) > 0.12_kind_lake) .and. (snowdp2d(i) <= 0.18_kind_lake)) then
5552 snl2d(i) = -3
5553 dz3d(i,-2) = 0.02_kind_lake
5554 dz3d(i,-1) = 0.05_kind_lake
5555 dz3d(i, 0) = snowdp2d(i) - dz3d(i,-2) - dz3d(i,-1)
5556 else if ((snowdp2d(i) > 0.18_kind_lake) .and. (snowdp2d(i) <= 0.29_kind_lake)) then
5557 snl2d(i) = -4
5558 dz3d(i,-3) = 0.02_kind_lake
5559 dz3d(i,-2) = 0.05_kind_lake
5560 dz3d(i,-1) = (snowdp2d(i) - dz3d(i,-3) - dz3d(i,-2))*0.5_kind_lake
5561 dz3d(i, 0) = dz3d(i,-1)
5562 else if ((snowdp2d(i) > 0.29_kind_lake) .and. (snowdp2d(i) <= 0.41_kind_lake)) then
5563 snl2d(i) = -4
5564 dz3d(i,-3) = 0.02_kind_lake
5565 dz3d(i,-2) = 0.05_kind_lake
5566 dz3d(i,-1) = 0.11_kind_lake
5567 dz3d(i, 0) = snowdp2d(i) - dz3d(i,-3) - dz3d(i,-2) - dz3d(i,-1)
5568 else if ((snowdp2d(i) > 0.41_kind_lake) .and. (snowdp2d(i) <= 0.64_kind_lake)) then
5569 snl2d(i) = -5
5570 dz3d(i,-4) = 0.02_kind_lake
5571 dz3d(i,-3) = 0.05_kind_lake
5572 dz3d(i,-2) = 0.11_kind_lake
5573 dz3d(i,-1) = (snowdp2d(i) - dz3d(i,-4) - dz3d(i,-3) - dz3d(i,-2))*0.5_kind_lake
5574 dz3d(i, 0) = dz3d(i,-1)
5575 else if (snowdp2d(i) > 0.64_kind_lake) then
5576 snl2d(i) = -5
5577 dz3d(i,-4) = 0.02_kind_lake
5578 dz3d(i,-3) = 0.05_kind_lake
5579 dz3d(i,-2) = 0.11_kind_lake
5580 dz3d(i,-1) = 0.23_kind_lake
5581 dz3d(i, 0)=snowdp2d(i)-dz3d(i,-4)-dz3d(i,-3)-dz3d(i,-2)-dz3d(i,-1)
5582 endif
5583 end if
5584
5585 do k = 0, nint(snl2d(i)+1), -1
5586 z3d(i,k) = zi3d(i,k) - 0.5_kind_lake*dz3d(i,k)
5587 zi3d(i,k-1) = zi3d(i,k) - dz3d(i,k)
5588 end do
5589
5590 ! 3:subroutine makearbinit
5591
5592 ! initial t_lake3d here
5593 if(tsfc(i)<160) then
5594 write(errmsg,'(A,F20.12,A)') 'Invalid tsfc value ',tsfc(i),' found. Was tsfc not initialized?'
5595 write(0,'(A)') trim(errmsg)
5596 errflg=1
5597 return
5598 endif
5599
5600 if(lake_icefrac3d(i,1) > 0.) then
5601 depth = 0.
5602 do k=2,nlevlake
5603 depth = depth + dz_lake(k)
5604 if(hice(i) >= depth) then
5605 lake_icefrac3d(i,k) = max(0.,lake_icefrac3d(i,1)+(0.-lake_icefrac3d(i,1))/z_lake(nlevlake)*depth)
5606 else
5607 lake_icefrac3d(i,k) = 0.
5608 endif
5609 end do
5610 endif
5611
5612 t_lake3d(i,1) = tsfc(i)
5613 t_grnd2d(i) = tsfc(i)
5614 if (lake_icefrac3d(i,1) <= 0.) then
5615 t_lake3d(i,1) = max(tfrz,tsfc(i))
5616 t_grnd2d(i) = max(tfrz,tsfc(i))
5617 endif
5618 do k = 2, nlevlake
5619 if(z_lake(k).le.depth_c) then
5620 t_lake3d(i,k) = tsfc(i)+(277.2_kind_lake-tsfc(i))/depth_c*z_lake(k)
5621 else
5622 t_lake3d(i,k) = 277.2_kind_lake
5623 end if
5624 if (lake_icefrac3d(i,k) <= 0.) then
5625 t_lake3d(i,k) = max(tfrz,t_lake3d(i,k))
5626 endif
5627 enddo
5628
5629 ! initial t_soisno3d
5630 ! in snow
5631 if(snowdp2d(i) > 0.) then
5632 do k = snl2d(i)+1, 0
5633 t_soisno3d(i,k) =min(tfrz,tsfc(i))
5634 enddo
5635 endif
5636
5637 ! in soil
5638 t_soisno3d(i,1) = t_lake3d(i,nlevlake)
5639 t_soisno3d(i,nlevsoil) = tg3(i)
5640 do k = 2, nlevsoil-1
5641 t_soisno3d(i,k)=t_soisno3d(i,1)+(t_soisno3d(i,nlevsoil)-t_soisno3d(i,1))*dzsoi(k)
5642 enddo
5643
5644 if (snl2d(i) < 0) then
5645 do k = nint(snl2d(i)+1), 0
5646 ! Be careful because there may be new snow layers with bad temperatures like 0 even if
5647 ! coming from init. con. file.
5648 if(t_soisno3d(i,k) > 300 .or. t_soisno3d(i,k) < 200) t_soisno3d(i,k) = min(tfrz,tsfc(i))
5649 enddo
5650 end if
5651
5652 do k = 1,nlevsoil
5653 h2osoi_vol3d(i,k) = 1.0_kind_lake
5654 watsat = 0.489_kind_lake - 0.00126_kind_lake*sand(isl)
5655 h2osoi_vol3d(i,k) = min(h2osoi_vol3d(i,k),watsat)
5656
5657 ! soil layers
5658 if (t_soisno3d(i,k) <= tfrz) then
5659 h2osoi_ice3d(i,k) = dz3d(i,k)*denice*h2osoi_vol3d(i,k)
5660 h2osoi_liq3d(i,k) = 0._kind_lake
5661 else
5662 h2osoi_ice3d(i,k) = 0._kind_lake
5663 h2osoi_liq3d(i,k) = dz3d(i,k)*denh2o*h2osoi_vol3d(i,k)
5664 endif
5665 enddo
5666
5667 !tgs - in RAP and HRRR applications with cycled snow depth and snow
5668 !water equivalent, the actual snow density could be computed. This is
5669 !not used for now for consistency with the main Lake subroutine, where
5670 !constant snow density (250.) is used.
5671 if(h2osno2d(i).gt.0. .and. snowdp2d(i).gt.0.) then
5672 rhosn = h2osno2d(i)/snowdp2d(i)
5673 else
5674 rhosn = snow_bd ! bdsno=250.
5675 endif
5676
5677
5678 do k = -nlevsnow+1, 0
5679 if (k > snl2d(i)) then
5680 h2osoi_ice3d(i,k) = dz3d(i,k)*snow_bd
5681 h2osoi_liq3d(i,k) = 0._kind_lake
5682 end if
5683 end do
5684
5685 clm_lake_initialized(i) = 1
5686
5687 endif !if ( use_lakedepth ) then
5688 ENDDO do_init
5689
5690
5691 if(debug_print .and. init_points>0) then
5692 print *,'points initialized in clm_lake',init_points
5693 end if
5694
5695END SUBROUTINE lakeini
5696
5697END MODULE clm_lake
combo
subroutine combo(parameters, dz, wliq, wice, t, dz2, wliq2, wice2, t2)
Definition module_sf_noahmplsm.F90:7916
snowwater
subroutine snowwater(parameters, nsnow, nsoil, imelt, dt, zsoil, sfctmp, snowhin, qsnow, qsnfro, qsnsub, qrain, ficeold, iloc, jloc, isnow, snowh, sneqv, snice, snliq, sh2o, sice, stc, zsnso, dzsnso, qsnbot, snoflow, ponding1, ponding2)
Definition module_sf_noahmplsm.F90:7391
albedo
subroutine albedo(parameters, vegtyp, ist, ice, nsoil, dt, cosz, fage, elai, esai, tg, tv, snowh, fsno, fwet, smc, sneqvo, sneqv, qsnow, fveg, iloc, jloc, albold, tauss, albgrd, albgri, albd, albi, fabd, fabi, ftdd, ftid, ftii, fsun, frevi, frevd, fregd, fregi, bgap, wgap, albsnd, albsni)
surface albedos. also fluxes (per unit incoming direct and diffuse radiation) reflected,...
Definition module_sf_noahmplsm.F90:2899
clm_lake
This module contains the CLM Lake model.
Definition clm_lake.f90:22