CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
module_sf_noahmp_glacier.F90
1#define CCPP
2
4
6module noahmp_glacier_globals
7
8 use machine , only : kind_phys
9 use sfc_diff, only : stability
10 use module_sf_noahmplsm, only : sfcdif4
11
12 implicit none
13
14! ==================================================================================================
15!------------------------------------------------------------------------------------------!
16! physical constants: !
17!------------------------------------------------------------------------------------------!
18
19 real (kind=kind_phys), parameter :: grav = 9.80616
20 real (kind=kind_phys), parameter :: sb = 5.67e-08
21 real (kind=kind_phys), parameter :: vkc = 0.40
22 real (kind=kind_phys), parameter :: tfrz = 273.16
23 real (kind=kind_phys), parameter :: hsub = 2.8440e06
24 real (kind=kind_phys), parameter :: hvap = 2.5104e06
25 real (kind=kind_phys), parameter :: hfus = 0.3336e06
26 real (kind=kind_phys), parameter :: cwat = 4.188e06
27 real (kind=kind_phys), parameter :: cice = 2.094e06
28 real (kind=kind_phys), parameter :: cpair = 1004.64
29 real (kind=kind_phys), parameter :: tkwat = 0.6
30 real (kind=kind_phys), parameter :: tkice = 2.2
31 real (kind=kind_phys), parameter :: tkair = 0.023
32 real (kind=kind_phys), parameter :: rair = 287.04
33 real (kind=kind_phys), parameter :: rw = 461.269
34 real (kind=kind_phys), parameter :: denh2o = 1000.
35 real (kind=kind_phys), parameter :: denice = 917.
36
37! =====================================options for different schemes================================
38
41
42 INTEGER :: OPT_ALB != 2 !(suggested 2)
43
46
47 INTEGER :: OPT_SNF != 1 !(suggested 1)
48
52
53 INTEGER :: OPT_TBOT != 2 !(suggested 2)
54
57
58 INTEGER :: OPT_STC != 1 !(suggested 1)
59
62
63 INTEGER :: OPT_GLA != 1 !(suggested 1)
64
65 INTEGER :: OPT_SFC != 1 !(suggested 1)
66 INTEGER :: OPT_TRS != 1 !(suggested 2)
67
68! adjustable parameters for snow processes
69
70 REAL, PARAMETER :: Z0SNO = 0.002
71 REAL, PARAMETER :: SSI = 0.03
72 REAL, PARAMETER :: SWEMX = 1.00
74
75!------------------------------------------------------------------------------------------!
76end module noahmp_glacier_globals
77!------------------------------------------------------------------------------------------!
78
81 use noahmp_glacier_globals
82#ifndef CCPP
83 use module_wrf_utl
84#endif
85 implicit none
86
88 public :: noahmp_glacier
89
90 private :: atm_glacier
91 private :: energy_glacier
92 private :: thermoprop_glacier
93 private :: csnow_glacier
94 private :: radiation_glacier
95 private :: snow_age_glacier
96 private :: snowalb_bats_glacier
97 private :: snowalb_class_glacier
98 private :: glacier_flux
99 private :: sfcdif1_glacier
100 private :: tsnosoi_glacier
101 private :: hrt_glacier
102 private :: hstep_glacier
103 private :: rosr12_glacier
104 private :: phasechange_glacier
105
106 private :: water_glacier
107 private :: snowwater_glacier
108 private :: snowfall_glacier
109 private :: combine_glacier
110 private :: divide_glacier
111 private :: combo_glacier
112 private :: compact_glacier
113 private :: snowh2o_glacier
114
115 private :: error_glacier
116
117contains
118!
119! ==================================================================================================
120
122 subroutine noahmp_glacier (&
123 iloc ,jloc ,cosz ,nsnow ,nsoil ,dt , & ! in : time/space/model-related
124 sfctmp ,sfcprs ,uu ,vv ,q2 ,soldn , & ! in : forcing
125 prcp ,lwdn ,tbot ,zlvl ,ficeold ,zsoil , & ! in : forcing
126 thsfc_loc ,prslkix ,prsik1x ,prslk1x , &
127 psfc ,pblhx ,iz0tlnd ,itime , &
128 sigmaf1 ,garea1 ,psi_opt , & ! in :
129 ep_1 ,ep_2 ,epsm1 ,cp , &
130 qsnow ,sneqvo ,albold ,cm ,ch ,isnow , & ! in/out :
131 sneqv ,smc ,zsnso ,snowh ,snice ,snliq , & ! in/out :
132 tg ,stc ,sh2o ,tauss ,qsfc , & ! in/out :
133 fsa ,fsr ,fira ,fsh ,fgev ,ssoil , & ! out :
134 trad ,edir ,runsrf ,runsub ,sag ,albedo , & ! out :
135 qsnbot ,ponding ,ponding1 ,ponding2 ,t2m,q2e ,z0h_total , & ! out :
136#ifdef CCPP
137 emissi ,fpice ,ch2b , esnow , albsnd , albsni , &
138 errmsg ,errflg)
139#else
140 emissi ,fpice ,ch2b , esnow. , albsnd , albsni)
141#endif
142
143
144! --------------------------------------------------------------------------------------------------
145! initial code: guo-yue niu, oct. 2007
146! modified to glacier: michael barlage, june 2012
147! --------------------------------------------------------------------------------------------------
148 implicit none
149! --------------------------------------------------------------------------------------------------
150! input
151 integer , intent(in) :: iloc
152 integer , intent(in) :: jloc
153 real (kind=kind_phys) , intent(in) :: cosz
154 integer , intent(in) :: nsnow
155 integer , intent(in) :: nsoil
156 integer , intent(in) :: psi_opt
157
158 real (kind=kind_phys) , intent(in) :: dt
159 real (kind=kind_phys) , intent(in) :: sfctmp
160 real (kind=kind_phys) , intent(in) :: sfcprs
161 real (kind=kind_phys) , intent(in) :: uu
162 real (kind=kind_phys) , intent(in) :: vv
163 real (kind=kind_phys) , intent(in) :: q2
164 real (kind=kind_phys) , intent(in) :: soldn
165 real (kind=kind_phys) , intent(in) :: prcp
166 real (kind=kind_phys) , intent(in) :: lwdn
167 real (kind=kind_phys) , intent(in) :: tbot
168 real (kind=kind_phys) , intent(in) :: zlvl
169 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(in) :: ficeold
170 real (kind=kind_phys), dimension( 1:nsoil), intent(in) :: zsoil
171 logical , intent(in) :: thsfc_loc
172 real (kind=kind_phys) , intent(in) :: prslkix
173 real (kind=kind_phys) , intent(in) :: prsik1x
174 real (kind=kind_phys) , intent(in) :: prslk1x
175
176 real (kind=kind_phys) , intent(in) :: psfc ! surface pressure
177 real (kind=kind_phys) , intent(in) :: pblhx ! pbl height
178 real (kind=kind_phys) , intent(in) :: ep_1
179 real (kind=kind_phys) , intent(in) :: ep_2
180 real (kind=kind_phys) , intent(in) :: epsm1
181 real (kind=kind_phys) , intent(in) :: cp
182 integer , intent(in) :: iz0tlnd !
183 integer , intent(in) :: itime
184
185 real (kind=kind_phys) , intent(in) :: sigmaf1
186 real (kind=kind_phys) , intent(in) :: garea1
187
188
189
190! input/output : need arbitary intial values
191 real (kind=kind_phys) , intent(inout) :: qsnow
192 real (kind=kind_phys) , intent(inout) :: sneqvo
193 real (kind=kind_phys) , intent(inout) :: albold
194 real (kind=kind_phys) , intent(inout) :: cm
195 real (kind=kind_phys) , intent(inout) :: ch
196
197! prognostic variables
198 integer , intent(inout) :: isnow
199 real (kind=kind_phys) , intent(inout) :: sneqv
200 real (kind=kind_phys), dimension( 1:nsoil), intent(inout) :: smc
201 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: zsnso
202 real (kind=kind_phys) , intent(inout) :: snowh
203 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(inout) :: snice
204 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(inout) :: snliq
205 real (kind=kind_phys) , intent(inout) :: tg
206 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: stc
207 real (kind=kind_phys), dimension( 1:nsoil), intent(inout) :: sh2o
208 real (kind=kind_phys) , intent(inout) :: tauss
209 real (kind=kind_phys) , intent(inout) :: qsfc
210
211! output
212 real (kind=kind_phys) , intent(out) :: fsa
213 real (kind=kind_phys) , intent(out) :: fsr
214 real (kind=kind_phys) , intent(out) :: fira
215 real (kind=kind_phys) , intent(out) :: fsh
216 real (kind=kind_phys) , intent(out) :: fgev
217 real (kind=kind_phys) , intent(out) :: ssoil
218 real (kind=kind_phys) , intent(out) :: trad
219 real (kind=kind_phys) , intent(out) :: edir
220 real (kind=kind_phys) , intent(out) :: runsrf
221 real (kind=kind_phys) , intent(out) :: runsub
222 real (kind=kind_phys) , intent(out) :: sag
223 real (kind=kind_phys) , intent(out) :: albedo
224 real (kind=kind_phys) , intent(out) :: qsnbot
225 real (kind=kind_phys) , intent(out) :: ponding
226 real (kind=kind_phys) , intent(out) :: ponding1
227 real (kind=kind_phys) , intent(out) :: ponding2
228 real (kind=kind_phys) , intent(out) :: t2m
229 real (kind=kind_phys) , intent(out) :: q2e
230 real (kind=kind_phys) , intent(out) :: z0h_total
231 real (kind=kind_phys) , intent(out) :: emissi
232 real (kind=kind_phys) , intent(out) :: fpice
233 real (kind=kind_phys) , intent(out) :: ch2b
234 real (kind=kind_phys) , intent(out) :: esnow
235 real (kind=kind_phys), dimension(1:2) , intent(out) :: albsnd
236 real (kind=kind_phys), dimension(1:2) , intent(out) :: albsni
237
238
239#ifdef CCPP
240 character(len=*), intent(inout) :: errmsg
241 integer, intent(inout) :: errflg
242#endif
243
244! local
245 integer :: iz
246 integer, dimension(-nsnow+1:nsoil) :: imelt
247 real (kind=kind_phys) :: rhoair
248 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: dzsnso
249 real (kind=kind_phys) :: thair
250 real (kind=kind_phys) :: qair
251 real (kind=kind_phys) :: eair
252 real (kind=kind_phys), dimension( 1: 2) :: solad
253 real (kind=kind_phys), dimension( 1: 2) :: solai
254 real (kind=kind_phys), dimension( 1:nsoil) :: sice
255 real (kind=kind_phys), dimension(-nsnow+1: 0) :: snicev
256 real (kind=kind_phys), dimension(-nsnow+1: 0) :: snliqv
257 real (kind=kind_phys), dimension(-nsnow+1: 0) :: epore
258 real (kind=kind_phys) :: qdew
259 real (kind=kind_phys) :: qvap
260 real (kind=kind_phys) :: lathea
261 real (kind=kind_phys) :: qmelt
262 real (kind=kind_phys) :: swdown
263 real (kind=kind_phys) :: beg_wb
264 real (kind=kind_phys) :: zbot = -8.0
265
266 character*256 message
267
268! --------------------------------------------------------------------------------------------------
269! re-process atmospheric forcing
270
271 call atm_glacier (ep_2, epsm1,sfcprs ,sfctmp ,q2 ,soldn ,cosz ,thair , &
272 qair ,eair ,rhoair ,solad ,solai ,swdown )
273
274 beg_wb = sneqv
275
276! snow/soil layer thickness (m); interface depth: zsnso < 0; layer thickness dzsnso > 0
277
278 do iz = isnow+1, nsoil
279 if(iz == isnow+1) then
280 dzsnso(iz) = - zsnso(iz)
281 else
282 dzsnso(iz) = zsnso(iz-1) - zsnso(iz)
283 end if
284 end do
285
286! compute energy budget (momentum & energy fluxes and phase changes)
287
288 call energy_glacier (nsnow ,nsoil ,isnow ,dt ,qsnow ,rhoair , & !in
289 eair ,sfcprs ,qair ,sfctmp ,lwdn ,uu , & !in
290 vv ,solad ,solai ,cosz ,zlvl , & !in
291 tbot ,zbot ,zsnso ,dzsnso ,sigmaf1 ,garea1 , & !in
292 thsfc_loc ,prslkix ,prsik1x ,prslk1x , & !in
293 psfc ,pblhx ,iz0tlnd ,itime ,psi_opt , &
294 ep_1, ep_2, epsm1, cp, &
295 tg ,stc ,snowh ,sneqv ,sneqvo ,sh2o , & !inout
296 smc ,snice ,snliq ,albold ,cm ,ch , & !inout
297#ifdef CCPP
298 tauss ,qsfc ,errmsg ,errflg , & !inout
299#else
300 tauss ,qsfc , & !inout
301#endif
302 imelt ,snicev ,snliqv ,epore ,qmelt ,ponding , & !out
303 sag ,fsa ,fsr ,fira ,fsh ,fgev , & !out
304 trad ,t2m ,ssoil ,lathea ,q2e ,emissi , & !out
305 ch2b ,albsnd ,albsni ,z0h_total) !out
306
307#ifdef CCPP
308 if (errflg /= 0) return
309#endif
310
311 sice = max(0.0, smc - sh2o)
312 sneqvo = sneqv
313
314 qvap = max( fgev/lathea, 0.) ! positive part of fgev [mm/s] > 0
315 qdew = abs( min(fgev/lathea, 0.)) ! negative part of fgev [mm/s] > 0
316 edir = qvap - qdew
317
318! compute water budgets (water storages, et components, and runoff)
319
320 call water_glacier (nsnow ,nsoil ,imelt ,dt ,prcp ,sfctmp , & !in
321 qvap ,qdew ,ficeold ,zsoil , & !in
322 isnow ,snowh ,sneqv ,snice ,snliq ,stc , & !inout
323 dzsnso ,sh2o ,sice ,ponding ,zsnso ,fsh , & !inout
324 runsrf ,runsub ,qsnow ,ponding1 ,ponding2 ,qsnbot , & !out
325 fpice ,esnow) !out
326
327 if(opt_gla == 2) then
328 edir = qvap - qdew
329 fgev = edir * lathea
330 end if
331
332! if(maxval(sice) < 0.0001) then
333! write(message,*) "glacier has melted at:",iloc,jloc," are you sure this should be a glacier point?"
334! call wrf_debug(10,trim(message))
335! end if
336
337! water and energy balance check
338
339 call error_glacier (iloc ,jloc ,swdown ,fsa ,fsr ,fira , &
340 fsh ,fgev ,ssoil ,sag ,prcp ,edir , &
341#ifdef CCPP
342 runsrf ,runsub ,sneqv ,dt ,beg_wb ,errmsg , errflg )
343#else
344 runsrf ,runsub ,sneqv ,dt ,beg_wb )
345#endif
346
347#ifdef CCPP
348 if (errflg /= 0) return
349#endif
350
351 if(snowh <= 1.e-6 .or. sneqv <= 1.e-3) then
352 snowh = 0.0
353 sneqv = 0.0
354 end if
355
356 if(swdown.ne.0.) then
357 albedo = fsr / swdown
358 else
359 albedo = -999.9
360 end if
361
362
363 end subroutine noahmp_glacier
364! ==================================================================================================
367 subroutine atm_glacier (ep_2, epsm1, sfcprs ,sfctmp ,q2 ,soldn ,cosz ,thair , &
368 qair ,eair ,rhoair ,solad ,solai , &
369 swdown )
370! --------------------------------------------------------------------------------------------------
371! re-process atmospheric forcing
372! --------------------------------------------------------------------------------------------------
373 implicit none
374! --------------------------------------------------------------------------------------------------
375! inputs
376
377 real (kind=kind_phys) , intent(in) :: ep_2
378 real (kind=kind_phys) , intent(in) :: epsm1
379 real (kind=kind_phys) , intent(in) :: sfcprs
380 real (kind=kind_phys) , intent(in) :: sfctmp
381 real (kind=kind_phys) , intent(in) :: q2
382 real (kind=kind_phys) , intent(in) :: soldn
383 real (kind=kind_phys) , intent(in) :: cosz
384
385! outputs
386
387 real (kind=kind_phys) , intent(out) :: thair
388 real (kind=kind_phys) , intent(out) :: qair
389 real (kind=kind_phys) , intent(out) :: eair
390 real (kind=kind_phys), dimension( 1: 2), intent(out) :: solad
391 real (kind=kind_phys), dimension( 1: 2), intent(out) :: solai
392 real (kind=kind_phys) , intent(out) :: rhoair
393 real (kind=kind_phys) , intent(out) :: swdown
394
395!locals
396
397 real (kind=kind_phys) :: pair
398! --------------------------------------------------------------------------------------------------
399
400 pair = sfcprs ! atm bottom level pressure (pa)
401 thair = sfctmp * (sfcprs/pair)**(rair/cpair)
402! qair = q2 / (1.0+q2) ! mixing ratio to specific humidity [kg/kg]
403 qair = q2 ! in wrf, driver converts to specific humidity
404
405 eair = qair*sfcprs / (ep_2-epsm1*qair)
406 rhoair = (sfcprs+epsm1*eair) / (rair*sfctmp)
407
408 if(cosz <= 0.) then
409 swdown = 0.
410 else
411 swdown = soldn
412 end if
413
414 solad(1) = swdown*0.7*0.5 ! direct vis
415 solad(2) = swdown*0.7*0.5 ! direct nir
416 solai(1) = swdown*0.3*0.5 ! diffuse vis
417 solai(2) = swdown*0.3*0.5 ! diffuse nir
418
419 end subroutine atm_glacier
420! ==================================================================================================
421! --------------------------------------------------------------------------------------------------
424 subroutine energy_glacier (nsnow ,nsoil ,isnow ,dt ,qsnow ,rhoair , & !in
425 eair ,sfcprs ,qair ,sfctmp ,lwdn ,uu , & !in
426 vv ,solad ,solai ,cosz ,zref , & !in
427 tbot ,zbot ,zsnso ,dzsnso ,sigmaf1 ,garea1 , & !in
428 thsfc_loc ,prslkix ,prsik1x ,prslk1x , & !in
429 psfc ,pblhx ,iz0tlnd ,itime ,psi_opt , &
430 ep_1, ep_2, epsm1, cp, &
431 tg ,stc ,snowh ,sneqv ,sneqvo ,sh2o , & !inout
432 smc ,snice ,snliq ,albold ,cm ,ch , & !inout
433#ifdef CCPP
434 tauss ,qsfc ,errmsg ,errflg , & !inout
435#else
436 tauss ,qsfc , & !inout
437#endif
438 imelt ,snicev ,snliqv ,epore ,qmelt ,ponding , & !out
439 sag ,fsa ,fsr ,fira ,fsh ,fgev , & !out
440 trad ,t2m ,ssoil ,lathea ,q2e ,emissi , & !out
441 ch2b ,albsnd ,albsni ,z0h_total) !out
442
443! --------------------------------------------------------------------------------------------------
444! --------------------------------------------------------------------------------------------------
445! use noahmp_veg_parameters
446! use noahmp_rad_parameters
447! --------------------------------------------------------------------------------------------------
448 implicit none
449! --------------------------------------------------------------------------------------------------
450! inputs
451 integer , intent(in) :: nsnow
452 integer , intent(in) :: nsoil
453 integer , intent(in) :: psi_opt
454
455 integer , intent(in) :: isnow
456 real (kind=kind_phys) , intent(in) :: dt
457 real (kind=kind_phys) , intent(in) :: qsnow
458 real (kind=kind_phys) , intent(in) :: rhoair
459 real (kind=kind_phys) , intent(in) :: eair
460 real (kind=kind_phys) , intent(in) :: sfcprs
461 real (kind=kind_phys) , intent(in) :: qair
462 real (kind=kind_phys) , intent(in) :: sfctmp
463 real (kind=kind_phys) , intent(in) :: lwdn
464 real (kind=kind_phys) , intent(in) :: uu
465 real (kind=kind_phys) , intent(in) :: vv
466 real (kind=kind_phys) , dimension( 1: 2), intent(in) :: solad
467 real (kind=kind_phys) , dimension( 1: 2), intent(in) :: solai
468 real (kind=kind_phys) , intent(in) :: cosz
469 real (kind=kind_phys) , intent(in) :: zref
470 real (kind=kind_phys) , intent(in) :: tbot
471 real (kind=kind_phys) , intent(in) :: zbot
472 real (kind=kind_phys) , dimension(-nsnow+1:nsoil), intent(in) :: zsnso
473 real (kind=kind_phys) , dimension(-nsnow+1:nsoil), intent(in) :: dzsnso
474
475 logical , intent(in) :: thsfc_loc
476 real (kind=kind_phys) , intent(in) :: prslkix ! in exner function
477 real (kind=kind_phys) , intent(in) :: prsik1x ! in exner function
478 real (kind=kind_phys) , intent(in) :: prslk1x ! in exner function
479
480 real (kind=kind_phys) , intent(in) :: pblhx
481 real (kind=kind_phys) , intent(in) :: psfc
482 real (kind=kind_phys) , intent(in) :: ep_1
483 real (kind=kind_phys) , intent(in) :: ep_2
484 real (kind=kind_phys) , intent(in) :: epsm1
485 real (kind=kind_phys) , intent(in) :: cp
486 integer , intent(in) :: iz0tlnd
487 integer , intent(in) :: itime
488
489 real (kind=kind_phys) , intent(in) :: sigmaf1
490 real (kind=kind_phys) , intent(in) :: garea1
491
492! input & output
493 real (kind=kind_phys) , intent(inout) :: tg
494 real (kind=kind_phys) , dimension(-nsnow+1:nsoil), intent(inout) :: stc
495 real (kind=kind_phys) , intent(inout) :: snowh
496 real (kind=kind_phys) , intent(inout) :: sneqv
497 real (kind=kind_phys) , intent(inout) :: sneqvo
498 real (kind=kind_phys) , dimension( 1:nsoil), intent(inout) :: sh2o
499 real (kind=kind_phys) , dimension( 1:nsoil), intent(inout) :: smc
500 real (kind=kind_phys) , dimension(-nsnow+1: 0), intent(inout) :: snice
501 real (kind=kind_phys) , dimension(-nsnow+1: 0), intent(inout) :: snliq
502 real (kind=kind_phys) , intent(inout) :: albold
503 real (kind=kind_phys) , intent(inout) :: cm
504 real (kind=kind_phys) , intent(inout) :: ch
505 real (kind=kind_phys) , intent(inout) :: tauss
506 real (kind=kind_phys) , intent(inout) :: qsfc
507
508#ifdef CCPP
509 character(len=*) , intent(inout) :: errmsg
510 integer , intent(inout) :: errflg
511#endif
512
513! outputs
514 integer, dimension(-nsnow+1:nsoil) , intent(out) :: imelt
515 real (kind=kind_phys) , dimension(-nsnow+1: 0), intent(out) :: snicev
516 real (kind=kind_phys) , dimension(-nsnow+1: 0), intent(out) :: snliqv
517 real (kind=kind_phys) , dimension(-nsnow+1: 0), intent(out) :: epore
518 real (kind=kind_phys) , intent(out) :: qmelt
519 real (kind=kind_phys) , intent(out) :: ponding
520 real (kind=kind_phys) , intent(out) :: sag
521 real (kind=kind_phys) , intent(out) :: fsa
522 real (kind=kind_phys) , intent(out) :: fsr
523 real (kind=kind_phys) , intent(out) :: fira
524 real (kind=kind_phys) , intent(out) :: fsh
525 real (kind=kind_phys) , intent(out) :: fgev
526 real (kind=kind_phys) , intent(out) :: trad
527 real (kind=kind_phys) , intent(out) :: t2m
528 real (kind=kind_phys) , intent(out) :: ssoil
529 real (kind=kind_phys) , intent(out) :: lathea
530 real (kind=kind_phys) , intent(out) :: q2e
531 real (kind=kind_phys) , intent(out) :: emissi
532 real (kind=kind_phys) , intent(out) :: ch2b
533 real (kind=kind_phys), dimension(1:2) , intent(out) :: albsnd
534 real (kind=kind_phys), dimension(1:2) , intent(out) :: albsni
535 real (kind=kind_phys) , intent(out) :: z0h_total
536
537
538! local
539 real (kind=kind_phys) :: ur
540 real (kind=kind_phys) :: zlvl
541 real (kind=kind_phys) :: rsurf
542 real (kind=kind_phys) :: zpd
543 real (kind=kind_phys) :: z0mg
544 real (kind=kind_phys) :: emg
545 real (kind=kind_phys) :: fire
546 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: fact
547 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: df
548 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: hcpct
549 real (kind=kind_phys) :: gamma
550 real (kind=kind_phys) :: rhsur
551
552! ---------------------------------------------------------------------------------------------------
553
554! wind speed at reference height: ur >= 1
555
556 ur = max( sqrt(uu**2.+vv**2.), 1. )
557
558! roughness length and displacement height
559
560 z0mg = z0sno
561 zpd = snowh
562
563 zlvl = zpd + zref
564
565! thermal properties of soil, snow, lake, and frozen soil
566
567 call thermoprop_glacier (nsoil ,nsnow ,isnow ,dzsnso , & !in
568 dt ,snowh ,snice ,snliq , & !in
569 df ,hcpct ,snicev ,snliqv ,epore , & !out
570 fact ) !out
571
572! solar radiation: absorbed & reflected by the ground
573
574 call radiation_glacier (dt ,tg ,sneqvo ,sneqv ,cosz , & !in
575 qsnow ,solad ,solai , & !in
576 albold ,tauss , & !inout
577 sag ,fsr ,fsa ,albsnd ,albsni) !out
578
579! vegetation and ground emissivity
580
581 emg = 0.98
582
583! soil surface resistance for ground evap.
584
585 rhsur = 1.0
586 rsurf = 1.0
587
588! set psychrometric constant
589
590 lathea = hsub
591 gamma = cpair*sfcprs/(ep_2*lathea)
592
593! surface temperatures of the ground and energy fluxes
594
595 call glacier_flux (nsoil ,nsnow ,emg ,isnow ,df ,dzsnso ,z0mg , & !in
596 zlvl ,zpd ,qair ,sfctmp ,rhoair ,sfcprs , & !in
597 ur ,gamma ,rsurf ,lwdn ,rhsur ,smc , & !in
598 eair ,stc ,sag ,snowh ,lathea ,sh2o , & !in
599 thsfc_loc ,prslkix ,prsik1x ,prslk1x , &
600 psfc ,pblhx ,iz0tlnd ,itime ,uu ,vv , &
601 sigmaf1 ,garea1 ,psi_opt ,ep_1, ep_2, epsm1, cp, & !in
602#ifdef CCPP
603 cm ,ch ,tg ,qsfc ,errmsg ,errflg , & !inout
604#else
605 cm ,ch ,tg ,qsfc , & !inout
606#endif
607 fira ,fsh ,fgev ,ssoil , & !out
608 t2m ,q2e ,ch2b ,z0h_total) !out
609
610!energy balance at surface: sag=(irb+shb+evb+ghb)
611
612 fire = lwdn + fira
613
614 if(fire <=0.) then
615#ifdef CCPP
616 errflg = 1
617 errmsg = "stop in noah-mp: emitted longwave <0"
618 return
619#else
620 call wrf_error_fatal("stop in noah-mp: emitted longwave <0")
621#endif
622 end if
623
624 ! compute a net emissivity
625 emissi = emg
626
627 ! when we're computing a trad, subtract from the emitted ir the
628 ! reflected portion of the incoming lwdn, so we're just
629 ! considering the ir originating in the canopy/ground system.
630
631 trad = ( ( fire - (1-emissi)*lwdn ) / (emissi*sb) ) ** 0.25
632
633! 3l snow & 4l soil temperatures
634
635 call tsnosoi_glacier (nsoil ,nsnow ,isnow ,dt ,tbot , & !in
636 ssoil ,snowh ,zbot ,zsnso ,df , & !in
637 hcpct , & !in
638 stc ) !inout
639
640! adjusting snow surface temperature
641 if(opt_stc == 2) then
642 if (snowh > 0.05 .and. tg > tfrz) tg = tfrz
643 end if
644
645! energy released or consumed by snow & ice
646
647 call phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & !in
648 dzsnso , & !in
649 stc ,snice ,snliq ,sneqv ,snowh , & !inout
650 smc ,sh2o , & !inout
651 qmelt ,imelt ,ponding ) !out
652
653
654 end subroutine energy_glacier
655! ==================================================================================================
658 subroutine thermoprop_glacier (nsoil ,nsnow ,isnow ,dzsnso , & !in
659 dt ,snowh ,snice ,snliq , & !in
660 df ,hcpct ,snicev ,snliqv ,epore , & !out
661 fact ) !out
662! -------------------------------------------------------------------------------------------------
663! -------------------------------------------------------------------------------------------------
664 implicit none
665! --------------------------------------------------------------------------------------------------
666! inputs
667 integer , intent(in) :: nsoil
668 integer , intent(in) :: nsnow
669 integer , intent(in) :: isnow
670 real (kind=kind_phys) , intent(in) :: dt
671 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(in) :: snice
672 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(in) :: snliq
673 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: dzsnso
674 real (kind=kind_phys) , intent(in) :: snowh
675
676! outputs
677 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(out) :: df
678 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(out) :: hcpct
679 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(out) :: snicev
680 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(out) :: snliqv
681 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(out) :: epore
682 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(out) :: fact
683! --------------------------------------------------------------------------------------------------
684! locals
685
686 integer :: iz, iz2
687 real (kind=kind_phys), dimension(-nsnow+1: 0) :: cvsno
688 real (kind=kind_phys), dimension(-nsnow+1: 0) :: tksno
689 real (kind=kind_phys) :: zmid
690! --------------------------------------------------------------------------------------------------
691
692! compute snow thermal conductivity and heat capacity
693
694 call csnow_glacier (isnow ,nsnow ,nsoil ,snice ,snliq ,dzsnso , & !in
695 tksno ,cvsno ,snicev ,snliqv ,epore ) !out
696
697 do iz = isnow+1, 0
698 df(iz) = tksno(iz)
699 hcpct(iz) = cvsno(iz)
700 end do
701
702! compute soil thermal properties (using noah glacial ice approximations)
703
704 do iz = 1, nsoil
705 zmid = 0.5 * (dzsnso(iz))
706 do iz2 = 1, iz-1
707 zmid = zmid + dzsnso(iz2)
708 end do
709 hcpct(iz) = 1.e6 * ( 0.8194 + 0.1309*zmid )
710 df(iz) = 0.32333 + ( 0.10073 * zmid )
711 end do
712
713! combine a temporary variable used for melting/freezing of snow and frozen soil
714
715 do iz = isnow+1,nsoil
716 fact(iz) = dt/(hcpct(iz)*dzsnso(iz))
717 end do
718
719! snow/soil interface
720
721 if(isnow == 0) then
722 df(1) = (df(1)*dzsnso(1)+0.35*snowh) / (snowh +dzsnso(1))
723 else
724 df(1) = (df(1)*dzsnso(1)+df(0)*dzsnso(0)) / (dzsnso(0)+dzsnso(1))
725 end if
726
727
728 end subroutine thermoprop_glacier
729! ==================================================================================================
730! --------------------------------------------------------------------------------------------------
733 subroutine csnow_glacier (isnow ,nsnow ,nsoil ,snice ,snliq ,dzsnso , & !in
734 tksno ,cvsno ,snicev ,snliqv ,epore ) !out
735! --------------------------------------------------------------------------------------------------
736! snow bulk density,volumetric capacity, and thermal conductivity
737!---------------------------------------------------------------------------------------------------
738 implicit none
739!---------------------------------------------------------------------------------------------------
740! inputs
741
742 integer, intent(in) :: isnow
743 integer , intent(in) :: nsnow
744 integer , intent(in) :: nsoil
745 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(in) :: snice
746 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(in) :: snliq
747 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: dzsnso
748
749! outputs
750
751 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(out) :: cvsno
752 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(out) :: tksno
753 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(out) :: snicev
754 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(out) :: snliqv
755 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(out) :: epore
756
757! locals
758
759 integer :: iz
760 real (kind=kind_phys), dimension(-nsnow+1: 0) :: bdsnoi
761
762!---------------------------------------------------------------------------------------------------
763! thermal capacity of snow
764
765 do iz = isnow+1, 0
766 snicev(iz) = min(1., snice(iz)/(dzsnso(iz)*denice) )
767 epore(iz) = 1. - snicev(iz)
768 snliqv(iz) = min(epore(iz),snliq(iz)/(dzsnso(iz)*denh2o))
769 enddo
770
771 do iz = isnow+1, 0
772 bdsnoi(iz) = (snice(iz)+snliq(iz))/dzsnso(iz)
773 cvsno(iz) = cice*snicev(iz)+cwat*snliqv(iz)
774! cvsno(iz) = 0.525e06 ! constant
775 enddo
776
777! thermal conductivity of snow
778
779 do iz = isnow+1, 0
780! tksno(iz) = 3.2217e-6*bdsnoi(iz)**2. ! stieglitz(yen,1965)
781! tksno(iz) = 2e-2+2.5e-6*bdsnoi(iz)*bdsnoi(iz) ! anderson, 1976
782! tksno(iz) = 0.35 ! constant
783 tksno(iz) = 2.576e-6*bdsnoi(iz)**2. + 0.074 ! verseghy (1991)
784! tksno(iz) = 2.22*(bdsnoi(iz)/1000.)**1.88 ! douvill(yen, 1981)
785 enddo
786
787 end subroutine csnow_glacier
788!===================================================================================================
791 subroutine radiation_glacier (dt ,tg ,sneqvo ,sneqv ,cosz , & !in
792 qsnow ,solad ,solai , & !in
793 albold ,tauss , & !inout
794 sag ,fsr ,fsa ,albsnd ,albsni) !out
795! --------------------------------------------------------------------------------------------------
796 implicit none
797! --------------------------------------------------------------------------------------------------
798! input
799 real (kind=kind_phys), intent(in) :: dt
800 real (kind=kind_phys), intent(in) :: tg
801 real (kind=kind_phys), intent(in) :: sneqvo
802 real (kind=kind_phys), intent(in) :: sneqv
803 real (kind=kind_phys), intent(in) :: cosz
804 real (kind=kind_phys), intent(in) :: qsnow
805 real (kind=kind_phys), dimension(1:2) , intent(in) :: solad
806 real (kind=kind_phys), dimension(1:2) , intent(in) :: solai
807
808! inout
809 real (kind=kind_phys), intent(inout) :: albold
810 real (kind=kind_phys), intent(inout) :: tauss
811 real (kind=kind_phys), dimension(1:2) :: albsnd
812 real (kind=kind_phys), dimension(1:2) :: albsni
813
814! output
815 real (kind=kind_phys), intent(out) :: sag
816 real (kind=kind_phys), intent(out) :: fsr
817 real (kind=kind_phys), intent(out) :: fsa
818
819! local
820 integer :: ib
821 integer :: nband
822 real (kind=kind_phys) :: fage
823 real (kind=kind_phys) :: alb
824 real (kind=kind_phys) :: abs
825 real (kind=kind_phys) :: ref
826 real (kind=kind_phys) :: fsno
827 real (kind=kind_phys), dimension(1:2) :: albice
828
829 real (kind=kind_phys),parameter :: mpe = 1.e-6
830
831! --------------------------------------------------------------------------------------------------
832
833 nband = 2
834 albsnd = 0.0
835 albsni = 0.0
836 albice(1) = 0.80 !albedo land ice: 1=vis, 2=nir
837 albice(2) = 0.55
838
839! snow age
840
841 call snow_age_glacier (dt,tg,sneqvo,sneqv,tauss,fage)
842
843! snow albedos: age even when sun is not present
844
845 if(opt_alb == 1) &
846 call snowalb_bats_glacier (nband,cosz,fage,albsnd,albsni)
847 if(opt_alb == 2) then
848 call snowalb_class_glacier(nband,qsnow,dt,alb,albold,albsnd,albsni)
849 albold = alb
850 end if
851
852! zero summed solar fluxes
853
854 sag = 0.
855 fsa = 0.
856 fsr = 0.
857
858 fsno = 0.0
859 if(sneqv > 0.0) fsno = 1.0
860
861! loop over nband wavebands
862
863 do ib = 1, nband
864
865 albsnd(ib) = albice(ib)*(1.-fsno) + albsnd(ib)*fsno
866 albsni(ib) = albice(ib)*(1.-fsno) + albsni(ib)*fsno
867
868! solar radiation absorbed by ground surface
869
870 abs = solad(ib)*(1.-albsnd(ib)) + solai(ib)*(1.-albsni(ib))
871 sag = sag + abs
872 fsa = fsa + abs
873
874 ref = solad(ib)*albsnd(ib) + solai(ib)*albsni(ib)
875 fsr = fsr + ref
876
877 end do
878
879 end subroutine radiation_glacier
880! ==================================================================================================
882 subroutine snow_age_glacier (dt,tg,sneqvo,sneqv,tauss,fage)
883! --------------------------------------------------------------------------------------------------
884 implicit none
885! ------------------------ code history ------------------------------------------------------------
886! from bats
887! ------------------------ input/output variables --------------------------------------------------
888!input
889 real (kind=kind_phys), intent(in) :: dt
890 real (kind=kind_phys), intent(in) :: tg
891 real (kind=kind_phys), intent(in) :: sneqvo
892 real (kind=kind_phys), intent(in) :: sneqv
893
894! inout
895 real (kind=kind_phys), intent(inout) :: tauss
896
897!output
898 real (kind=kind_phys), intent(out) :: fage
899
900!local
901 real (kind=kind_phys) :: tage
902 real (kind=kind_phys) :: age1
903 real (kind=kind_phys) :: age2
904 real (kind=kind_phys) :: age3
905 real (kind=kind_phys) :: dela
906 real (kind=kind_phys) :: sge
907 real (kind=kind_phys) :: dels
908 real (kind=kind_phys) :: dela0
909 real (kind=kind_phys) :: arg
910! see yang et al. (1997) j.of climate for detail.
911!---------------------------------------------------------------------------------------------------
912
913 if(sneqv.le.0.0) then
914 tauss = 0.
915 else if (sneqv.gt.800.) then
916 tauss = 0.
917 else
918! tauss = 0.
919 dela0 = 1.e-6*dt
920 arg = 5.e3*(1./tfrz-1./tg)
921 age1 = exp(arg)
922 age2 = exp(amin1(0.,10.*arg))
923 age3 = 0.3
924 tage = age1+age2+age3
925 dela = dela0*tage
926 dels = amax1(0.0,sneqv-sneqvo) / swemx
927 sge = (tauss+dela)*(1.0-dels)
928 tauss = amax1(0.,sge)
929 endif
930
931 fage= tauss/(tauss+1.)
932
933 end subroutine snow_age_glacier
934! ==================================================================================================
935! --------------------------------------------------------------------------------------------------
937 subroutine snowalb_bats_glacier (nband,cosz,fage,albsnd,albsni)
938! --------------------------------------------------------------------------------------------------
939 implicit none
940! --------------------------------------------------------------------------------------------------
941! input
942
943 integer,intent(in) :: nband
944
945 real (kind=kind_phys),intent(in) :: cosz
946 real (kind=kind_phys),intent(in) :: fage
947
948! output
949
950 real (kind=kind_phys), dimension(1:2),intent(out) :: albsnd
951 real (kind=kind_phys), dimension(1:2),intent(out) :: albsni
952! ---------------------------------------------------------------------------------------------
953
954 real (kind=kind_phys) :: fzen
955 real (kind=kind_phys) :: cf1
956 real (kind=kind_phys) :: sl2
957 real (kind=kind_phys) :: sl1
958 real (kind=kind_phys) :: sl
959 real (kind=kind_phys), parameter :: c1 = 0.2
960 real (kind=kind_phys), parameter :: c2 = 0.5
961! real (kind=kind_phys), parameter :: c1 = 0.2 * 2. !< double the default to match sleepers river's
962! real (kind=kind_phys), parameter :: c2 = 0.5 * 2. !< snow surface albedo (double aging effects)
963! ---------------------------------------------------------------------------------------------
964! zero albedos for all points
965
966 albsnd(1: nband) = 0.
967 albsni(1: nband) = 0.
968
969! when cosz > 0
970
971 sl=2.0
972 sl1=1./sl
973 sl2=2.*sl
974 cf1=((1.+sl1)/(1.+sl2*cosz)-sl1)
975 fzen=amax1(cf1,0.)
976
977 albsni(1)=0.95*(1.-c1*fage)
978 albsni(2)=0.65*(1.-c2*fage)
979
980 albsnd(1)=albsni(1)+0.4*fzen*(1.-albsni(1)) ! vis direct
981 albsnd(2)=albsni(2)+0.4*fzen*(1.-albsni(2)) ! nir direct
982
983 end subroutine snowalb_bats_glacier
984! ==================================================================================================
985! --------------------------------------------------------------------------------------------------
987 subroutine snowalb_class_glacier (nband,qsnow,dt,alb,albold,albsnd,albsni)
988! --------------------------------------------------------------------------------------------------
989 implicit none
990! --------------------------------------------------------------------------------------------------
991! input
992
993 integer,intent(in) :: nband
994
995 real (kind=kind_phys),intent(in) :: qsnow
996 real (kind=kind_phys),intent(in) :: dt
997 real (kind=kind_phys),intent(in) :: albold
998
999! in & out
1000
1001 real (kind=kind_phys), intent(inout) :: alb !
1002! output
1003
1004 real (kind=kind_phys), dimension(1:2),intent(out) :: albsnd
1005 real (kind=kind_phys), dimension(1:2),intent(out) :: albsni
1006! ---------------------------------------------------------------------------------------------
1007
1008! ---------------------------------------------------------------------------------------------
1009! zero albedos for all points
1010
1011 albsnd(1: nband) = 0.
1012 albsni(1: nband) = 0.
1013
1014! when cosz > 0
1015
1016 alb = 0.55 + (albold-0.55) * exp(-0.01*dt/3600.)
1017
1018! 1 mm fresh snow(swe) -- 10mm snow depth, assumed the fresh snow density 100kg/m3
1019! here assume 1cm snow depth will fully cover the old snow
1020
1021 if (qsnow > 0.) then
1022 alb = alb + min(qsnow*dt,swemx) * (0.84-alb)/(swemx)
1023 endif
1024
1025 albsni(1)= alb ! vis diffuse
1026 albsni(2)= alb ! nir diffuse
1027 albsnd(1)= alb ! vis direct
1028 albsnd(2)= alb ! nir direct
1029
1030 end subroutine snowalb_class_glacier
1031! ==================================================================================================
1035 subroutine glacier_flux (nsoil ,nsnow ,emg ,isnow ,df ,dzsnso ,z0m , & !in
1036 zlvl ,zpd ,qair ,sfctmp ,rhoair ,sfcprs , & !in
1037 ur ,gamma ,rsurf ,lwdn ,rhsur ,smc , & !in
1038 eair ,stc ,sag ,snowh ,lathea ,sh2o , & !in
1039 thsfc_loc ,prslkix ,prsik1x ,prslk1x , &
1040 psfc ,pblhx ,iz0tlnd ,itime ,uu ,vv , &
1041 sigmaf1 ,garea1 ,psi_opt ,ep_1, ep_2, epsm1, cp, & !in
1042#ifdef CCPP
1043 cm ,ch ,tgb ,qsfc ,errmsg ,errflg , & !inout
1044#else
1045 cm ,ch ,tgb ,qsfc , & !inout
1046#endif
1047 irb ,shb ,evb ,ghb , & !out
1048 t2mb ,q2b ,ehb2 ,z0h_total) !out
1049
1050! --------------------------------------------------------------------------------------------------
1051! use newton-raphson iteration to solve ground (tg) temperature
1052! that balances the surface energy budgets for glacier.
1053
1054! bare soil:
1055! -sab + irb[tg] + shb[tg] + evb[tg] + ghb[tg] = 0
1056! ----------------------------------------------------------------------
1057! use module_model_constants
1058! ----------------------------------------------------------------------
1059 implicit none
1060! ----------------------------------------------------------------------
1061! input
1062 integer, intent(in) :: nsnow
1063 integer, intent(in) :: nsoil
1064 integer, intent(in) :: psi_opt
1065
1066 real (kind=kind_phys), intent(in) :: emg
1067 integer, intent(in) :: isnow
1068 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: df
1069 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: dzsnso
1070 real (kind=kind_phys), intent(in) :: z0m
1071 real (kind=kind_phys), intent(in) :: zlvl
1072 real (kind=kind_phys), intent(in) :: zpd
1073 real (kind=kind_phys), intent(in) :: qair
1074 real (kind=kind_phys), intent(in) :: sfctmp
1075 real (kind=kind_phys), intent(in) :: rhoair
1076 real (kind=kind_phys), intent(in) :: sfcprs
1077 real (kind=kind_phys), intent(in) :: ur
1078 real (kind=kind_phys), intent(in) :: gamma
1079 real (kind=kind_phys), intent(in) :: rsurf
1080 real (kind=kind_phys), intent(in) :: lwdn
1081 real (kind=kind_phys), intent(in) :: rhsur
1082 real (kind=kind_phys), intent(in) :: eair
1083 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: stc
1084 real (kind=kind_phys), dimension( 1:nsoil), intent(in) :: smc
1085 real (kind=kind_phys), dimension( 1:nsoil), intent(in) :: sh2o
1086 real (kind=kind_phys), intent(in) :: sag
1087 real (kind=kind_phys), intent(in) :: snowh
1088 real (kind=kind_phys), intent(in) :: lathea
1089
1090 logical , intent(in) :: thsfc_loc !way to th tmp
1091 real (kind=kind_phys), intent(in) :: prslkix ! in exner function
1092 real (kind=kind_phys), intent(in) :: prsik1x ! in exner function
1093 real (kind=kind_phys), intent(in) :: prslk1x ! in exner function
1094
1095 real (kind=kind_phys) , intent(in) :: pblhx
1096 real (kind=kind_phys) , intent(in) :: psfc
1097 real (kind=kind_phys) , intent(in) :: ep_1
1098 real (kind=kind_phys) , intent(in) :: ep_2
1099 real (kind=kind_phys) , intent(in) :: epsm1
1100 real (kind=kind_phys) , intent(in) :: cp
1101 integer , intent(in) :: iz0tlnd
1102 integer , intent(in) :: itime
1103 real (kind=kind_phys) , intent(in) :: uu
1104 real (kind=kind_phys) , intent(in) :: vv
1105
1106 real (kind=kind_phys), intent(in) :: sigmaf1 !
1107 real (kind=kind_phys), intent(in) :: garea1 !
1108
1109! input/output
1110 real (kind=kind_phys), intent(inout) :: cm
1111 real (kind=kind_phys), intent(inout) :: ch
1112 real (kind=kind_phys), intent(inout) :: tgb
1113 real (kind=kind_phys), intent(inout) :: qsfc
1114
1115#ifdef CCPP
1116 character(len=*), intent(inout) :: errmsg
1117 integer, intent(inout) :: errflg
1118#endif
1119
1120! output
1121! -sab + irb[tg] + shb[tg] + evb[tg] + ghb[tg] = 0
1122 real (kind=kind_phys), intent(out) :: irb
1123 real (kind=kind_phys), intent(out) :: shb
1124 real (kind=kind_phys), intent(out) :: evb
1125 real (kind=kind_phys), intent(out) :: ghb
1126 real (kind=kind_phys), intent(out) :: t2mb
1127 real (kind=kind_phys), intent(out) :: q2b
1128 real (kind=kind_phys), intent(out) :: ehb2
1129 real (kind=kind_phys), intent(out) :: z0h_total
1130
1131
1132! local variables
1133 integer :: niterb
1134 integer :: niter
1135
1136 real (kind=kind_phys) :: mpe
1137 real (kind=kind_phys) :: dtg
1138 integer :: mozsgn
1139 real (kind=kind_phys) :: mozold
1140 real (kind=kind_phys) :: fm2
1141 real (kind=kind_phys) :: fh2
1142 real (kind=kind_phys) :: ch2
1143 real (kind=kind_phys) :: h
1144 real (kind=kind_phys) :: fv
1145 real (kind=kind_phys) :: cir
1146 real (kind=kind_phys) :: cgh
1147 real (kind=kind_phys) :: csh
1148 real (kind=kind_phys) :: cev
1149 real (kind=kind_phys) :: cq2b
1150 integer :: iter
1151 real (kind=kind_phys) :: z0h
1152
1153 real (kind=kind_phys) :: qfx
1154 real (kind=kind_phys) :: cq2
1155
1156
1157 real(kind=kind_phys) :: rb1i ! bulk richardson #
1158 real(kind=kind_phys) :: fm10i ! fm10 over land ice
1159
1160 real(kind=kind_phys) :: stress1i! wind stress m2 S-2
1161
1162 real(kind=kind_phys) :: wspd1i
1163 real(kind=kind_phys) :: flhc1i
1164 real(kind=kind_phys) :: flqc1i
1165
1166 real(kind=kind_phys) :: tv1i ! virtual potential temp @ ref level
1167
1168 real(kind=kind_phys) :: thv1i ! virtual potential temp @ ref level
1169 real(kind=kind_phys) :: tvsi ! surface virtual temp
1170 real(kind=kind_phys) :: zlvli ! ref. level
1171
1172 real(kind=kind_phys) :: snwd ! snow depth in mm
1173
1174 real(kind=kind_phys) :: reyni ! roughness Reynolds #
1175 real(kind=kind_phys) :: virtfaci! virutal factor
1176
1177 real(kind=kind_phys) :: tem1,tem2,zvfun1,gdx
1178 real(kind=kind_phys), parameter :: z0lo=0.1, z0up=1.0
1179
1180 real (kind=kind_phys) :: moz
1181 real (kind=kind_phys) :: fm
1182 real (kind=kind_phys) :: fh
1183 real (kind=kind_phys) :: ramb
1184 real (kind=kind_phys) :: rahb
1185 real (kind=kind_phys) :: rawb
1186 real (kind=kind_phys) :: estg
1187 real (kind=kind_phys) :: destg
1188 real (kind=kind_phys) :: esatw
1189 real (kind=kind_phys) :: esati
1190 real (kind=kind_phys) :: dsatw
1191 real (kind=kind_phys) :: dsati
1192 real (kind=kind_phys) :: a
1193 real (kind=kind_phys) :: b
1194 real (kind=kind_phys) :: t, tdc
1195 real (kind=kind_phys), dimension( 1:nsoil) :: sice
1196 real (kind=kind_phys) :: czil
1197
1198 tdc(t) = min( 50., max(-50.,(t-tfrz)) )
1199 czil=0.1
1200
1201! -----------------------------------------------------------------
1202! initialization variables that do not depend on stability iteration
1203! -----------------------------------------------------------------
1204 niterb = 5
1205 niter = 1
1206
1207 mpe = 1e-6
1208 dtg = 0.
1209 mozsgn = 0
1210 mozold = 0.
1211 moz = 0.
1212
1213 h = 0.
1214
1215 fh2 = 0.
1216 qfx = 0.
1217
1218
1219! the following only applies to opt_sfc =3, opt_sfc = 1 still done its old way
1220
1221 snwd = snowh*1000.0
1222 zlvli = zlvl - zpd
1223
1224! fv = ustarx ! the input maybe too high for glacial
1225 fv = ur*vkc/log(zlvli/z0m)
1226 reyni = fv*z0m/(1.5e-05) !introduction of fv dependent z0h for the iter
1227
1228 if (opt_trs == 1) then
1229 z0h = z0m
1230 elseif (opt_trs == 2) then
1231 z0h = z0m*exp(-czil*0.4*258.2*sqrt(fv*z0m))
1232 elseif (opt_trs == 3) then
1233 z0h = z0m*0.1
1234 elseif (opt_trs == 4) then
1235 if (reyni .gt. 2.0) then
1236 z0h = z0m/exp(2.46*(reyni)**0.25 - log(7.4)) !Brutsaert 1982
1237 else
1238 z0h = z0m/exp(-log(0.397)) !Brusaert 1982, table 4
1239 endif
1240 endif
1241
1242 z0h_total = z0h
1243
1244 virtfaci = 1.0 + 0.61 * max(qair, 1.e-8)
1245 tv1i = sfctmp * virtfaci ! virt tmp @ middle
1246
1247 if(thsfc_loc) then ! Use local potential temperature
1248 thv1i = sfctmp * prslkix * virtfaci
1249 else ! Use potential temperature reference to 1000 hPa
1250 thv1i = sfctmp / prslk1x * virtfaci
1251 endif
1252
1253 if ( ur < 2.0) niter = 2
1254
1255 cir = emg*sb
1256 cgh = 2.*df(isnow+1)/dzsnso(isnow+1)
1257
1258! -----------------------------------------------------------------
1259 tem1 = (z0m - z0lo) / (z0up - z0lo)
1260 tem1 = min(max(tem1, 0.0_kind_phys), 1.0_kind_phys)
1261 tem2 = max(sigmaf1, 0.1_kind_phys)
1262 zvfun1= sqrt(tem1 * tem2)
1263 gdx=sqrt(garea1)
1264
1265 if(opt_sfc == 1 .or. opt_sfc == 2 .or. opt_sfc == 4) then !Add option for sfc scheme,use '1' for both '1'/'2'
1266 loop3: do iter = 1, niterb ! begin stability iteration
1267 if(opt_sfc == 1 .or. opt_sfc == 2) then
1268
1269! for now, only allow sfcdif1 until others can be fixed
1270
1271 call sfcdif1_glacier(iter ,zlvl ,zpd ,z0h ,z0m , & !in
1272 qair ,sfctmp ,h ,rhoair ,mpe ,ur , & !in
1273#ifdef CCPP
1274 & moz ,mozsgn ,fm ,fh ,fm2 ,fh2,fv, errmsg, errflg, & !inout
1275#else
1276 & moz ,mozsgn ,fm ,fh ,fm2 ,fh2,fv , & !inout
1277#endif
1278 & cm ,ch ,ch2) !out
1279
1280#ifdef CCPP
1281 if (errflg /= 0) return
1282#endif
1283 endif
1284
1285 if(opt_sfc == 4) then
1286
1287 call sfcdif4(1 ,1 ,uu ,vv ,sfctmp , & !allow location for use in the driver
1288 sfcprs ,psfc ,pblhx ,gdx ,z0m , &
1289 ep_1, ep_2, cp, &
1290 itime ,snwd ,1 ,psi_opt, &
1291 tgb ,qair ,zlvl ,iz0tlnd,qsfc , & ! use zlvli?
1292 h ,qfx ,cm ,ch ,ch2 , & ! ch2 = cq2 most of times
1293 cq2 ,moz ,fv ,rb1i, fm, fh, &
1294 stress1i,fm10i ,fh2 , wspd1i ,flhc1i ,flqc1i) ! some are for use in the driver call
1295
1296
1297 ! Undo the multiplication by windspeed that SFCDIF4
1298 ! applies to exchange coefficients CH and CM:
1299
1300 ch = ch / wspd1i
1301 cm = cm / wspd1i
1302 ch2 = ch2 / wspd1i
1303 cq2 = cq2 / wspd1i
1304
1305 if(snwd > 0.) then
1306 cm = min(0.01,cm)
1307 ch = min(0.01,ch)
1308 ch2 = min(0.01,ch2)
1309 cq2 = min(0.01,cq2)
1310 end if
1311
1312 endif ! 4
1313
1314 if(opt_sfc == 1 .or. opt_sfc == 2 .or. opt_sfc == 3) then
1315 ramb = max(1.,1./(cm*ur))
1316 rahb = max(1.,1./(ch*ur))
1317 elseif(opt_sfc == 4) then
1318 ramb = max(1.,1./(cm*wspd1i) )
1319 rahb = max(1.,1./(ch*wspd1i) )
1320 endif
1321
1322 rawb = rahb
1323
1324! es and d(es)/dt evaluated at tg
1325
1326 t = tdc(tgb)
1327 call esat(t, esatw, esati, dsatw, dsati)
1328 if (t .gt. 0.) then
1329 estg = esatw
1330 destg = dsatw
1331 else
1332 estg = esati
1333 destg = dsati
1334 end if
1335
1336 csh = rhoair*cpair/rahb
1337 if(snowh > 0.0 .or. opt_gla == 1) then
1338 cev = rhoair*cpair/gamma/(rsurf+rawb)
1339 else
1340 cev = 0.0 ! don't allow any sublimation of glacier in opt_gla=2
1341 end if
1342
1343! surface fluxes and dtg
1344
1345 irb = cir * tgb**4 - emg*lwdn
1346 shb = csh * (tgb - sfctmp )
1347 evb = cev * (estg*rhsur - eair )
1348 ghb = cgh * (tgb - stc(isnow+1))
1349
1350 b = sag-irb-shb-evb-ghb
1351 a = 4.*cir*tgb**3 + csh + cev*destg + cgh
1352 dtg = b/a
1353
1354 irb = irb + 4.*cir*tgb**3*dtg
1355 shb = shb + csh*dtg
1356 evb = evb + cev*destg*dtg
1357 ghb = ghb + cgh*dtg
1358
1359! update ground surface temperature
1360 tgb = tgb + dtg
1361
1362! for m-o length
1363 h = csh * (tgb - sfctmp)
1364
1365 t = tdc(tgb)
1366 call esat(t, esatw, esati, dsatw, dsati)
1367 if (t .gt. 0.) then
1368 estg = esatw
1369 else
1370 estg = esati
1371 end if
1372 qsfc = ep_2*(estg*rhsur)/(sfcprs+epsm1*(estg*rhsur))
1373 qfx = (qsfc-qair)*cev*gamma/cpair
1374
1375 end do loop3 ! end stability iteration
1376 end if
1377
1378 if (opt_sfc == 3) then
1379
1380 do iter = 1, niter
1381
1382 if(thsfc_loc) then ! Use local potential temperature
1383 tvsi = tgb * virtfaci
1384 else ! Use potential temperature referenced to 1000 hPa
1385 tvsi = tgb/prsik1x * virtfaci
1386 endif
1387
1388 call stability &
1389 (zlvli, zvfun1, gdx,tv1i,thv1i, ur, z0m, z0h, tvsi, grav,thsfc_loc, &
1390 rb1i, fm,fh,fm10i,fh2,cm,ch,stress1i,fv)
1391
1392! maybe need to add some sorts of err handling if CCPP
1393
1394 ramb = max(1.,1./(cm*ur))
1395 rahb = max(1.,1./(ch*ur))
1396 rawb = rahb
1397
1398! es and d(es)/dt evaluated at tg
1399
1400 t = tdc(tgb)
1401 call esat(t, esatw, esati, dsatw, dsati)
1402 if (t .gt. 0.) then
1403 estg = esatw
1404 destg = dsatw
1405 else
1406 estg = esati
1407 destg = dsati
1408 end if
1409
1410 csh = rhoair*cpair/rahb
1411
1412 if(snowh > 0.0 .or. opt_gla == 1) then
1413 cev = rhoair*cpair/gamma/(rsurf+rawb)
1414 else
1415 cev = 0.0 ! don't allow any sublimation of glacier in opt_gla=2
1416 end if
1417
1418! surface fluxes and dtg
1419
1420 irb = cir * tgb**4 - emg*lwdn
1421 shb = csh * (tgb - sfctmp )
1422 evb = cev * (estg*rhsur - eair )
1423 ghb = cgh * (tgb - stc(isnow+1))
1424
1425 b = sag-irb-shb-evb-ghb
1426 a = 4.*cir*tgb**3 + csh + cev*destg + cgh
1427 dtg = b/a
1428
1429 irb = irb + 4.*cir*tgb**3*dtg
1430 shb = shb + csh*dtg
1431 evb = evb + cev*destg*dtg
1432 ghb = ghb + cgh*dtg
1433
1434! update ground surface temperature to update cm/ch
1435
1436 tgb = tgb + dtg
1437
1438 t = tdc(tgb)
1439 call esat(t, esatw, esati, dsatw, dsati)
1440 if (t .gt. 0.) then
1441 estg = esatw
1442 else
1443 estg = esati
1444 end if
1445 qsfc = ep_2*(estg*rhsur)/(sfcprs+epsm1*(estg*rhsur))
1446
1447 end do !sfc_diff3 iter
1448 end if !sfc_diff3
1449
1450! -----------------------------------------------------------------
1451
1452! if snow on ground and tg > tfrz: reset tg = tfrz. reevaluate ground fluxes.
1453
1454 sice = smc - sh2o
1455 if(opt_stc == 1 .or. opt_stc ==3) then
1456 if ((maxval(sice) > 0.0 .or. snowh > 0.0) .and. tgb > tfrz .and. opt_gla == 1) then
1457 tgb = tfrz
1458 t = tdc(tgb) ! mb: recalculate estg
1459 call esat(t, esatw, esati, dsatw, dsati)
1460 estg = esati
1461 qsfc = ep_2*(estg*rhsur)/(sfcprs+epsm1*(estg*rhsur))
1462 irb = cir * tgb**4 - emg*lwdn
1463 shb = csh * (tgb - sfctmp)
1464 evb = cev * (estg*rhsur - eair ) !estg reevaluate ?
1465 ghb = sag - (irb+shb+evb)
1466 end if
1467 end if
1468
1469! 2m air temperature
1470 ehb2 = fv*vkc/(log((2.+z0h)/z0h)-fh2)
1471 cq2b = ehb2
1472! for opt_sfc 3
1473 if (opt_sfc ==3) then
1474 ehb2 = fv*vkc/fh2
1475 cq2b = ehb2
1476 endif
1477
1478 if (opt_sfc == 4) then
1479 ehb2 = ch2 * wspd1i ! need conductance,z0h from sfcdif4
1480 cq2b = cq2 * wspd1i ! conductance
1481 endif
1482
1483 if (ehb2.lt.1.e-5 ) then
1484 t2mb = tgb
1485 q2b = qsfc
1486 else
1487 t2mb = tgb - shb/(rhoair*cpair) * 1./ehb2
1488 q2b = qsfc - evb/(lathea*rhoair)*(1./cq2b + rsurf)
1489 endif
1490
1491! update ch
1492 ch = 1./rahb
1493
1494 end subroutine glacier_flux
1495! ==================================================================================================
1497 subroutine esat(t, esw, esi, desw, desi)
1498!---------------------------------------------------------------------------------------------------
1501 implicit none
1502!---------------------------------------------------------------------------------------------------
1503! in
1504
1505 real (kind=kind_phys), intent(in) :: t
1506
1507!out
1508
1509 real (kind=kind_phys), intent(out) :: esw
1510 real (kind=kind_phys), intent(out) :: esi
1511 real (kind=kind_phys), intent(out) :: desw
1512 real (kind=kind_phys), intent(out) :: desi
1513
1514! local
1515
1516 real (kind=kind_phys) :: a0,a1,a2,a3,a4,a5,a6
1517 real (kind=kind_phys) :: b0,b1,b2,b3,b4,b5,b6
1518 real (kind=kind_phys) :: c0,c1,c2,c3,c4,c5,c6
1519 real (kind=kind_phys) :: d0,d1,d2,d3,d4,d5,d6
1520
1521 parameter(a0=6.107799961 , a1=4.436518521e-01, &
1522 a2=1.428945805e-02, a3=2.650648471e-04, &
1523 a4=3.031240396e-06, a5=2.034080948e-08, &
1524 a6=6.136820929e-11)
1525
1526 parameter(b0=6.109177956 , b1=5.034698970e-01, &
1527 b2=1.886013408e-02, b3=4.176223716e-04, &
1528 b4=5.824720280e-06, b5=4.838803174e-08, &
1529 b6=1.838826904e-10)
1530
1531 parameter(c0= 4.438099984e-01, c1=2.857002636e-02, &
1532 c2= 7.938054040e-04, c3=1.215215065e-05, &
1533 c4= 1.036561403e-07, c5=3.532421810e-10, &
1534 c6=-7.090244804e-13)
1535
1536 parameter(d0=5.030305237e-01, d1=3.773255020e-02, &
1537 d2=1.267995369e-03, d3=2.477563108e-05, &
1538 d4=3.005693132e-07, d5=2.158542548e-09, &
1539 d6=7.131097725e-12)
1540
1541 esw = 100.*(a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*a6))))))
1542 esi = 100.*(b0+t*(b1+t*(b2+t*(b3+t*(b4+t*(b5+t*b6))))))
1543 desw = 100.*(c0+t*(c1+t*(c2+t*(c3+t*(c4+t*(c5+t*c6))))))
1544 desi = 100.*(d0+t*(d1+t*(d2+t*(d3+t*(d4+t*(d5+t*d6))))))
1545
1546 end subroutine esat
1547! ==================================================================================================
1550 subroutine sfcdif1_glacier(iter ,zlvl ,zpd ,z0h ,z0m , & !in
1551 qair ,sfctmp ,h ,rhoair ,mpe ,ur , & !in
1552#ifdef CCPP
1553 & moz ,mozsgn ,fm ,fh ,fm2 ,fh2 , fv, & !inout
1554 & errmsg ,errflg , & !inout
1555#else
1556 & moz ,mozsgn ,fm ,fh ,fm2 ,fh2 , fv, & !inout
1557#endif
1558 & cm ,ch ,ch2 ) !out
1559! -------------------------------------------------------------------------------------------------
1560! computing surface drag coefficient cm for momentum and ch for heat
1561! -------------------------------------------------------------------------------------------------
1562 implicit none
1563! -------------------------------------------------------------------------------------------------
1564! inputs
1565 integer, intent(in) :: iter
1566 real (kind=kind_phys), intent(in) :: zlvl
1567 real (kind=kind_phys), intent(in) :: zpd
1568 real (kind=kind_phys), intent(in) :: z0h
1569 real (kind=kind_phys), intent(in) :: z0m
1570 real (kind=kind_phys), intent(in) :: qair
1571 real (kind=kind_phys), intent(in) :: sfctmp
1572 real (kind=kind_phys), intent(in) :: h
1573 real (kind=kind_phys), intent(in) :: rhoair
1574 real (kind=kind_phys), intent(in) :: mpe
1575 real (kind=kind_phys), intent(in) :: ur
1576
1577! in & out
1578 real (kind=kind_phys), intent(inout) :: moz
1579 integer, intent(inout) :: mozsgn
1580 real (kind=kind_phys), intent(inout) :: fm
1581 real (kind=kind_phys), intent(inout) :: fh
1582 real (kind=kind_phys), intent(inout) :: fm2
1583 real (kind=kind_phys), intent(inout) :: fh2
1584 real (kind=kind_phys), intent(inout) :: fv
1585
1586#ifdef CCPP
1587 character(len=*), intent(inout) :: errmsg
1588 integer, intent(inout) :: errflg
1589#endif
1590
1591! outputs
1592 real (kind=kind_phys), intent(out) :: cm
1593 real (kind=kind_phys), intent(out) :: ch
1594 real (kind=kind_phys), intent(out) :: ch2
1595
1596! locals
1597 real (kind=kind_phys) :: mozold
1598 real (kind=kind_phys) :: tmpcm
1599 real (kind=kind_phys) :: tmpch
1600 real (kind=kind_phys) :: mol
1601 real (kind=kind_phys) :: tvir
1602 real (kind=kind_phys) :: tmp1,tmp2,tmp3
1603 real (kind=kind_phys) :: fmnew
1604 real (kind=kind_phys) :: fhnew
1605 real (kind=kind_phys) :: moz2
1606 real (kind=kind_phys) :: tmpcm2
1607 real (kind=kind_phys) :: tmpch2
1608 real (kind=kind_phys) :: fm2new
1609 real (kind=kind_phys) :: fh2new
1610 real (kind=kind_phys) :: tmp12,tmp22,tmp32
1611
1612 real (kind=kind_phys) :: cmfm, chfh, cm2fm2, ch2fh2
1613
1614
1615! -------------------------------------------------------------------------------------------------
1616! monin-obukhov stability parameter moz for next iteration
1617
1618 mozold = moz
1619
1620 if(zlvl <= zpd) then
1621 write(*,*) 'critical glacier problem: zlvl <= zpd; model stops', zlvl, zpd
1622#ifdef CCPP
1623 errflg = 1
1624 errmsg = "stop in noah-mp glacier"
1625 return
1626#else
1627 call wrf_error_fatal("stop in noah-mp glacier")
1628#endif
1629 endif
1630
1631 tmpcm = log((zlvl-zpd) / z0m)
1632 tmpch = log((zlvl-zpd) / z0h)
1633 tmpcm2 = log((2.0 + z0m) / z0m)
1634 tmpch2 = log((2.0 + z0h) / z0h)
1635
1636 if(iter == 1) then
1637 fv = 0.0
1638 moz = 0.0
1639 mol = 0.0
1640 moz2 = 0.0
1641 else
1642 tvir = (1. + 0.61*qair) * sfctmp
1643 tmp1 = vkc * (grav/tvir) * h/(rhoair*cpair)
1644 if (abs(tmp1) .le. mpe) tmp1 = mpe
1645 mol = -1. * fv**3 / tmp1
1646 moz = min( (zlvl-zpd)/mol, 1.)
1647 moz2 = min( (2.0 + z0h)/mol, 1.)
1648 endif
1649
1650! accumulate number of times moz changes sign.
1651
1652 if (mozold*moz .lt. 0.) mozsgn = mozsgn+1
1653 if (mozsgn .ge. 2) then
1654 moz = 0.
1655 fm = 0.
1656 fh = 0.
1657 moz2 = 0.
1658 fm2 = 0.
1659 fh2 = 0.
1660 endif
1661
1662! evaluate stability-dependent variables using moz from prior iteration
1663 if (moz .lt. 0.) then
1664 tmp1 = (1. - 16.*moz)**0.25
1665 tmp2 = log((1.+tmp1*tmp1)/2.)
1666 tmp3 = log((1.+tmp1)/2.)
1667 fmnew = 2.*tmp3 + tmp2 - 2.*atan(tmp1) + 1.5707963
1668 fhnew = 2*tmp2
1669
1670! 2-meter
1671 tmp12 = (1. - 16.*moz2)**0.25
1672 tmp22 = log((1.+tmp12*tmp12)/2.)
1673 tmp32 = log((1.+tmp12)/2.)
1674 fm2new = 2.*tmp32 + tmp22 - 2.*atan(tmp12) + 1.5707963
1675 fh2new = 2*tmp22
1676 else
1677 fmnew = -5.*moz
1678 fhnew = fmnew
1679 fm2new = -5.*moz2
1680 fh2new = fm2new
1681 endif
1682
1683! except for first iteration, weight stability factors for previous
1684! iteration to help avoid flip-flops from one iteration to the next
1685
1686 if (iter == 1) then
1687 fm = fmnew
1688 fh = fhnew
1689 fm2 = fm2new
1690 fh2 = fh2new
1691 else
1692 fm = 0.5 * (fm+fmnew)
1693 fh = 0.5 * (fh+fhnew)
1694 fm2 = 0.5 * (fm2+fm2new)
1695 fh2 = 0.5 * (fh2+fh2new)
1696 endif
1697
1698! exchange coefficients
1699
1700 fh = min(fh,0.9*tmpch)
1701 fm = min(fm,0.9*tmpcm)
1702 fh2 = min(fh2,0.9*tmpch2)
1703 fm2 = min(fm2,0.9*tmpcm2)
1704
1705 cmfm = tmpcm-fm
1706 chfh = tmpch-fh
1707 cm2fm2 = tmpcm2-fm2
1708 ch2fh2 = tmpch2-fh2
1709 if(abs(cmfm) <= mpe) cmfm = mpe
1710 if(abs(chfh) <= mpe) chfh = mpe
1711 if(abs(cm2fm2) <= mpe) cm2fm2 = mpe
1712 if(abs(ch2fh2) <= mpe) ch2fh2 = mpe
1713 cm = vkc*vkc/(cmfm*cmfm)
1714 ch = vkc*vkc/(cmfm*chfh)
1715 ch2 = vkc*vkc/(cm2fm2*ch2fh2)
1716
1717! friction velocity
1718
1719 fv = ur * sqrt(cm)
1720 ch2 = vkc*fv/ch2fh2
1721
1722 end subroutine sfcdif1_glacier
1723! ==================================================================================================
1725 subroutine tsnosoi_glacier (nsoil ,nsnow ,isnow ,dt ,tbot , & !in
1726 ssoil ,snowh ,zbot ,zsnso ,df , & !in
1727 hcpct , & !in
1728 stc ) !inout
1729! --------------------------------------------------------------------------------------------------
1733! --------------------------------------------------------------------------------------------------
1734 implicit none
1735! --------------------------------------------------------------------------------------------------
1736!input
1737
1738 integer, intent(in) :: nsoil
1739 integer, intent(in) :: nsnow
1740 integer, intent(in) :: isnow
1741
1742 real (kind=kind_phys), intent(in) :: dt
1743 real (kind=kind_phys), intent(in) :: tbot
1744 real (kind=kind_phys), intent(in) :: ssoil
1745 real (kind=kind_phys), intent(in) :: snowh
1746 real (kind=kind_phys), intent(in) :: zbot
1747 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: zsnso
1748 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: df
1749 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: hcpct
1750
1751!input and output
1752
1753 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: stc
1754
1755!local
1756
1757 integer :: iz
1758 real (kind=kind_phys) :: zbotsno
1759 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: ai, bi, ci, rhsts
1760 real (kind=kind_phys) :: eflxb
1761 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: phi
1762
1763! ----------------------------------------------------------------------
1764
1765! prescribe solar penetration into ice/snow
1766
1767 phi(isnow+1:nsoil) = 0.
1768
1769! adjust zbot from soil surface to zbotsno from snow surface
1770
1771 zbotsno = zbot - snowh !from snow surface
1772
1773! compute ice temperatures
1774
1775 call hrt_glacier (nsnow ,nsoil ,isnow ,zsnso , &
1776 stc ,tbot ,zbotsno ,df , &
1777 hcpct ,ssoil ,phi , &
1778 ai ,bi ,ci ,rhsts , &
1779 eflxb )
1780
1781 call hstep_glacier (nsnow ,nsoil ,isnow ,dt , &
1782 ai ,bi ,ci ,rhsts , &
1783 stc )
1784
1785 end subroutine tsnosoi_glacier
1786! ==================================================================================================
1787! ----------------------------------------------------------------------
1789 subroutine hrt_glacier (nsnow ,nsoil ,isnow ,zsnso , & !in
1790 stc ,tbot ,zbot ,df , & !in
1791 hcpct ,ssoil ,phi , & !in
1792 ai ,bi ,ci ,rhsts , & !out
1793 botflx ) !out
1794! ----------------------------------------------------------------------
1795! ----------------------------------------------------------------------
1799! ----------------------------------------------------------------------
1800 implicit none
1801! ----------------------------------------------------------------------
1802! input
1803
1804 integer, intent(in) :: nsoil
1805 integer, intent(in) :: nsnow
1806 integer, intent(in) :: isnow
1807 real (kind=kind_phys), intent(in) :: tbot
1808 real (kind=kind_phys), intent(in) :: zbot
1810 real (kind=kind_phys), intent(in) :: ssoil
1811 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: zsnso
1812 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: stc
1813 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: df
1814 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: hcpct
1815 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: phi
1816
1817! output
1818
1819 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(out) :: rhsts
1820 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(out) :: ai
1821 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(out) :: bi
1822 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(out) :: ci
1823 real (kind=kind_phys), intent(out) :: botflx
1824
1825! local
1826
1827 integer :: k
1828 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: ddz
1829 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: denom
1830 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: dtsdz
1831 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: eflux
1832 real (kind=kind_phys) :: temp1
1833! ----------------------------------------------------------------------
1834
1835 do k = isnow+1, nsoil
1836 if (k == isnow+1) then
1837 denom(k) = - zsnso(k) * hcpct(k)
1838 temp1 = - zsnso(k+1)
1839 ddz(k) = 2.0 / temp1
1840 dtsdz(k) = 2.0 * (stc(k) - stc(k+1)) / temp1
1841 eflux(k) = df(k) * dtsdz(k) - ssoil - phi(k)
1842 else if (k < nsoil) then
1843 denom(k) = (zsnso(k-1) - zsnso(k)) * hcpct(k)
1844 temp1 = zsnso(k-1) - zsnso(k+1)
1845 ddz(k) = 2.0 / temp1
1846 dtsdz(k) = 2.0 * (stc(k) - stc(k+1)) / temp1
1847 eflux(k) = (df(k)*dtsdz(k) - df(k-1)*dtsdz(k-1)) - phi(k)
1848 else if (k == nsoil) then
1849 denom(k) = (zsnso(k-1) - zsnso(k)) * hcpct(k)
1850 temp1 = zsnso(k-1) - zsnso(k)
1851 if(opt_tbot == 1) then
1852 botflx = 0.
1853 end if
1854 if(opt_tbot == 2) then
1855 dtsdz(k) = (stc(k) - tbot) / ( 0.5*(zsnso(k-1)+zsnso(k)) - zbot)
1856 botflx = -df(k) * dtsdz(k)
1857 end if
1858 eflux(k) = (-botflx - df(k-1)*dtsdz(k-1) ) - phi(k)
1859 end if
1860 end do
1861
1862 do k = isnow+1, nsoil
1863 if (k == isnow+1) then
1864 ai(k) = 0.0
1865 ci(k) = - df(k) * ddz(k) / denom(k)
1866 if (opt_stc == 1 .or. opt_stc == 3) then
1867 bi(k) = - ci(k)
1868 end if
1869 if (opt_stc == 2) then
1870 bi(k) = - ci(k) + df(k)/(0.5*zsnso(k)*zsnso(k)*hcpct(k))
1871 end if
1872 else if (k < nsoil) then
1873 ai(k) = - df(k-1) * ddz(k-1) / denom(k)
1874 ci(k) = - df(k ) * ddz(k ) / denom(k)
1875 bi(k) = - (ai(k) + ci(k))
1876 else if (k == nsoil) then
1877 ai(k) = - df(k-1) * ddz(k-1) / denom(k)
1878 ci(k) = 0.0
1879 bi(k) = - (ai(k) + ci(k))
1880 end if
1881 rhsts(k) = eflux(k)/ (-denom(k))
1882 end do
1883
1884 end subroutine hrt_glacier
1885! ==================================================================================================
1886! ----------------------------------------------------------------------
1888 subroutine hstep_glacier (nsnow ,nsoil ,isnow ,dt , & !in
1889 ai ,bi ,ci ,rhsts , & !inout
1890 stc ) !inout
1891! ----------------------------------------------------------------------
1893! ----------------------------------------------------------------------
1894 implicit none
1895! ----------------------------------------------------------------------
1896! input
1897
1898 integer, intent(in) :: nsoil
1899 integer, intent(in) :: nsnow
1900 integer, intent(in) :: isnow
1901 real (kind=kind_phys), intent(in) :: dt
1902
1903! output & input
1904 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: ai
1905 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: bi
1906 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: ci
1907 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: stc
1908 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: rhsts
1909
1910! local
1911 integer :: k
1912 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: rhstsin
1913 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: ciin
1914! ----------------------------------------------------------------------
1915
1916 do k = isnow+1,nsoil
1917 rhsts(k) = rhsts(k) * dt
1918 ai(k) = ai(k) * dt
1919 bi(k) = 1. + bi(k) * dt
1920 ci(k) = ci(k) * dt
1921 end do
1922
1923! copy values for input variables before call to rosr12
1924
1925 do k = isnow+1,nsoil
1926 rhstsin(k) = rhsts(k)
1927 ciin(k) = ci(k)
1928 end do
1929
1930! solve the tri-diagonal matrix equation
1931
1932 call rosr12_glacier (ci,ai,bi,ciin,rhstsin,rhsts,isnow+1,nsoil,nsnow)
1933
1934! update snow & soil temperature
1935
1936 do k = isnow+1,nsoil
1937 stc(k) = stc(k) + ci(k)
1938 end do
1939
1940 end subroutine hstep_glacier
1941! ==================================================================================================
1943 subroutine rosr12_glacier (p,a,b,c,d,delta,ntop,nsoil,nsnow)
1944! ----------------------------------------------------------------------
1945! subroutine rosr12
1946! ----------------------------------------------------------------------
1947! invert (solve) the tri-diagonal matrix problem shown below:
1948! ### ### ### ### ### ###
1949! #b(1), c(1), 0 , 0 , 0 , . . . , 0 # # # # #
1950! #a(2), b(2), c(2), 0 , 0 , . . . , 0 # # # # #
1951! # 0 , a(3), b(3), c(3), 0 , . . . , 0 # # # # d(3) #
1952! # 0 , 0 , a(4), b(4), c(4), . . . , 0 # # p(4) # # d(4) #
1953! # 0 , 0 , 0 , a(5), b(5), . . . , 0 # # p(5) # # d(5) #
1954! # . . # # . # = # . #
1955! # . . # # . # # . #
1956! # . . # # . # # . #
1957! # 0 , . . . , 0 , a(m-2), b(m-2), c(m-2), 0 # #p(m-2)# #d(m-2)#
1958! # 0 , . . . , 0 , 0 , a(m-1), b(m-1), c(m-1)# #p(m-1)# #d(m-1)#
1959! # 0 , . . . , 0 , 0 , 0 , a(m) , b(m) # # p(m) # # d(m) #
1960! ### ### ### ### ### ###
1961! ----------------------------------------------------------------------
1962 implicit none
1963
1964 integer, intent(in) :: ntop
1965 integer, intent(in) :: nsoil,nsnow
1966 integer :: k, kk
1967
1968 real (kind=kind_phys), dimension(-nsnow+1:nsoil),intent(in):: a, b, d
1969 real (kind=kind_phys), dimension(-nsnow+1:nsoil),intent(inout):: c,p,delta
1970
1971! ----------------------------------------------------------------------
1972! initialize eqn coef c for the lowest soil layer
1973! ----------------------------------------------------------------------
1974 c(nsoil) = 0.0
1975 p(ntop) = - c(ntop) / b(ntop)
1976! ----------------------------------------------------------------------
1977! solve the coefs for the 1st soil layer
1978! ----------------------------------------------------------------------
1979 delta(ntop) = d(ntop) / b(ntop)
1980! ----------------------------------------------------------------------
1981! solve the coefs for soil layers 2 thru nsoil
1982! ----------------------------------------------------------------------
1983 do k = ntop+1,nsoil
1984 p(k) = - c(k) * ( 1.0 / (b(k) + a(k) * p(k -1)) )
1985 delta(k) = (d(k) - a(k)* delta(k -1))* (1.0/ (b(k) + a(k)&
1986 * p(k -1)))
1987 end do
1988! ----------------------------------------------------------------------
1989! set p to delta for lowest soil layer
1990! ----------------------------------------------------------------------
1991 p(nsoil) = delta(nsoil)
1992! ----------------------------------------------------------------------
1993! adjust p for soil layers 2 thru nsoil
1994! ----------------------------------------------------------------------
1995 do k = ntop+1,nsoil
1996 kk = nsoil - k + (ntop-1) + 1
1997 p(kk) = p(kk) * p(kk +1) + delta(kk)
1998 end do
1999! ----------------------------------------------------------------------
2000 end subroutine rosr12_glacier
2001! ----------------------------------------------------------------------
2002! ==================================================================================================
2004 subroutine phasechange_glacier (nsnow ,nsoil ,isnow ,dt ,fact , & !in
2005 dzsnso , & !in
2006 stc ,snice ,snliq ,sneqv ,snowh , & !inout
2007 smc ,sh2o , & !inout
2008 qmelt ,imelt ,ponding ) !out
2009! ----------------------------------------------------------------------
2011! ----------------------------------------------------------------------
2012 implicit none
2013! ----------------------------------------------------------------------
2014! inputs
2015
2016 integer, intent(in) :: nsnow
2017 integer, intent(in) :: nsoil
2018 integer, intent(in) :: isnow
2019 real (kind=kind_phys), intent(in) :: dt
2020 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: fact
2021 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: dzsnso
2022
2023! inputs/outputs
2024
2025 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: stc
2026 real (kind=kind_phys), dimension(-nsnow+1:0) , intent(inout) :: snice
2027 real (kind=kind_phys), dimension(-nsnow+1:0) , intent(inout) :: snliq
2028 real (kind=kind_phys), intent(inout) :: sneqv
2029 real (kind=kind_phys), intent(inout) :: snowh
2030 real (kind=kind_phys), dimension( 1:nsoil), intent(inout) :: sh2o
2031 real (kind=kind_phys), dimension( 1:nsoil), intent(inout) :: smc
2032
2033! outputs
2034 real (kind=kind_phys), intent(out) :: qmelt
2035 integer, dimension(-nsnow+1:nsoil), intent(out) :: imelt
2036 real (kind=kind_phys), intent(out) :: ponding
2037
2038! local
2039
2040 integer :: j,k
2041 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: hm
2042 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: xm
2043 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: wmass0
2044 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: wice0
2045 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: wliq0
2046 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: mice
2047 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: mliq
2048 real (kind=kind_phys), dimension(-nsnow+1:nsoil) :: heatr
2049 real (kind=kind_phys) :: temp1
2050 real (kind=kind_phys) :: propor
2051 real (kind=kind_phys) :: xmf
2052
2053! ----------------------------------------------------------------------
2054! initialization
2055
2056 qmelt = 0.
2057 ponding = 0.
2058 xmf = 0.
2059
2060 do j = isnow+1,0 ! all snow layers
2061 mice(j) = snice(j)
2062 mliq(j) = snliq(j)
2063 end do
2064
2065 do j = isnow+1,0 ! all snow layers; do ice later
2066 imelt(j) = 0
2067 hm(j) = 0.
2068 xm(j) = 0.
2069 wice0(j) = mice(j)
2070 wliq0(j) = mliq(j)
2071 wmass0(j) = mice(j) + mliq(j)
2072 enddo
2073
2074 do j = isnow+1,0
2075 if (mice(j) > 0. .and. stc(j) >= tfrz) then ! melting
2076 imelt(j) = 1
2077 endif
2078 if (mliq(j) > 0. .and. stc(j) < tfrz) then ! freezing
2079 imelt(j) = 2
2080 endif
2081
2082 enddo
2083
2084! calculate the energy surplus and loss for melting and freezing
2085
2086 do j = isnow+1,0
2087 if (imelt(j) > 0) then
2088 hm(j) = (stc(j)-tfrz)/fact(j)
2089 stc(j) = tfrz
2090 endif
2091
2092 if (imelt(j) == 1 .and. hm(j) < 0.) then
2093 hm(j) = 0.
2094 imelt(j) = 0
2095 endif
2096 if (imelt(j) == 2 .and. hm(j) > 0.) then
2097 hm(j) = 0.
2098 imelt(j) = 0
2099 endif
2100 xm(j) = hm(j)*dt/hfus
2101 enddo
2102
2103! the rate of melting and freezing for snow without a layer, opt_gla==1 treated below
2104
2105if (opt_gla == 2) then
2106
2107 if (isnow == 0 .and. sneqv > 0. .and. stc(1) >= tfrz) then
2108 hm(1) = (stc(1)-tfrz)/fact(1) ! available heat
2109 stc(1) = tfrz ! set t to freezing
2110 xm(1) = hm(1)*dt/hfus ! total snow melt possible
2111
2112 temp1 = sneqv
2113 sneqv = max(0.,temp1-xm(1)) ! snow remaining
2114 propor = sneqv/temp1 ! fraction melted
2115 snowh = max(0.,propor * snowh) ! new snow height
2116 heatr(1) = hm(1) - hfus*(temp1-sneqv)/dt ! excess heat
2117 if (heatr(1) > 0.) then
2118 xm(1) = heatr(1)*dt/hfus
2119 stc(1) = stc(1) + fact(1)*heatr(1) ! re-heat ice
2120 else
2121 xm(1) = 0. ! heat used up
2122 hm(1) = 0.
2123 endif
2124 qmelt = max(0.,(temp1-sneqv))/dt ! melted snow rate
2125 xmf = hfus*qmelt ! melted snow energy
2126 ponding = temp1-sneqv ! melt water
2127 endif
2128
2129end if ! opt_gla == 2
2130
2131! the rate of melting and freezing for snow
2132
2133 do j = isnow+1,0
2134 if (imelt(j) > 0 .and. abs(hm(j)) > 0.) then
2135
2136 heatr(j) = 0.
2137 if (xm(j) > 0.) then
2138 mice(j) = max(0., wice0(j)-xm(j))
2139 heatr(j) = hm(j) - hfus*(wice0(j)-mice(j))/dt
2140 else if (xm(j) < 0.) then
2141 mice(j) = min(wmass0(j), wice0(j)-xm(j))
2142 heatr(j) = hm(j) - hfus*(wice0(j)-mice(j))/dt
2143 endif
2144
2145 mliq(j) = max(0.,wmass0(j)-mice(j))
2146
2147 if (abs(heatr(j)) > 0.) then
2148 stc(j) = stc(j) + fact(j)*heatr(j)
2149 if (mliq(j)*mice(j)>0.) stc(j) = tfrz
2150 endif
2151
2152 qmelt = qmelt + max(0.,(wice0(j)-mice(j)))/dt
2153
2154 endif
2155 enddo
2156
2157if (opt_gla == 1) then ! operate on the ice layers
2158
2159 do j = 1, nsoil ! all soil layers
2160 mliq(j) = sh2o(j) * dzsnso(j) * 1000.
2161 mice(j) = (smc(j) - sh2o(j)) * dzsnso(j) * 1000.
2162 end do
2163
2164 do j = 1,nsoil ! all layers
2165 imelt(j) = 0
2166 hm(j) = 0.
2167 xm(j) = 0.
2168 wice0(j) = mice(j)
2169 wliq0(j) = mliq(j)
2170 wmass0(j) = mice(j) + mliq(j)
2171 enddo
2172
2173 do j = 1,nsoil
2174 if (mice(j) > 0. .and. stc(j) >= tfrz) then ! melting
2175 imelt(j) = 1
2176 endif
2177 if (mliq(j) > 0. .and. stc(j) < tfrz) then ! freezing
2178 imelt(j) = 2
2179 endif
2180
2181 ! if snow exists, but its thickness is not enough to create a layer
2182 if (isnow == 0 .and. sneqv > 0. .and. j == 1) then
2183 if (stc(j) >= tfrz) then
2184 imelt(j) = 1
2185 endif
2186 endif
2187 enddo
2188
2189! calculate the energy surplus and loss for melting and freezing
2190
2191 do j = 1,nsoil
2192 if (imelt(j) > 0) then
2193 hm(j) = (stc(j)-tfrz)/fact(j)
2194 stc(j) = tfrz
2195 endif
2196
2197 if (imelt(j) == 1 .and. hm(j) < 0.) then
2198 hm(j) = 0.
2199 imelt(j) = 0
2200 endif
2201 if (imelt(j) == 2 .and. hm(j) > 0.) then
2202 hm(j) = 0.
2203 imelt(j) = 0
2204 endif
2205 xm(j) = hm(j)*dt/hfus
2206 enddo
2207
2208! the rate of melting and freezing for snow without a layer, needs more work.
2209
2210 if (isnow == 0 .and. sneqv > 0. .and. xm(1) > 0.) then
2211 temp1 = sneqv
2212 sneqv = max(0.,temp1-xm(1))
2213 propor = sneqv/temp1
2214 snowh = max(0.,propor * snowh)
2215 heatr(1) = hm(1) - hfus*(temp1-sneqv)/dt
2216 if (heatr(1) > 0.) then
2217 xm(1) = heatr(1)*dt/hfus
2218 hm(1) = heatr(1)
2219 imelt(1) = 1
2220 else
2221 xm(1) = 0.
2222 hm(1) = 0.
2223 imelt(1) = 0
2224 endif
2225 qmelt = max(0.,(temp1-sneqv))/dt
2226 xmf = hfus*qmelt
2227 ponding = temp1-sneqv
2228 endif
2229
2230! the rate of melting and freezing for soil
2231
2232 do j = 1,nsoil
2233 if (imelt(j) > 0 .and. abs(hm(j)) > 0.) then
2234
2235 heatr(j) = 0.
2236 if (xm(j) > 0.) then
2237 mice(j) = max(0., wice0(j)-xm(j))
2238 heatr(j) = hm(j) - hfus*(wice0(j)-mice(j))/dt
2239 else if (xm(j) < 0.) then
2240 mice(j) = min(wmass0(j), wice0(j)-xm(j))
2241 heatr(j) = hm(j) - hfus*(wice0(j)-mice(j))/dt
2242 endif
2243
2244 mliq(j) = max(0.,wmass0(j)-mice(j))
2245
2246 if (abs(heatr(j)) > 0.) then
2247 stc(j) = stc(j) + fact(j)*heatr(j)
2248 if (j <= 0) then ! snow
2249 if (mliq(j)*mice(j)>0.) stc(j) = tfrz
2250 end if
2251 endif
2252
2253 if (j > 0) xmf = xmf + hfus * (wice0(j)-mice(j))/dt
2254
2255 if (j < 1) then
2256 qmelt = qmelt + max(0.,(wice0(j)-mice(j)))/dt
2257 endif
2258 endif
2259 enddo
2260 heatr = 0.0
2261 xm = 0.0
2262
2263! deal with residuals in ice/soil
2264
2265! first remove excess heat by reducing temperature of layers
2266
2267 if (any(stc(1:4) > tfrz) .and. any(stc(1:4) < tfrz)) then
2268 do j = 1,nsoil
2269 if ( stc(j) > tfrz ) then
2270 heatr(j) = (stc(j)-tfrz)/fact(j)
2271 do k = 1,nsoil
2272 if (j .ne. k .and. stc(k) < tfrz .and. heatr(j) > 0.1) then
2273 heatr(k) = (stc(k)-tfrz)/fact(k)
2274 if (abs(heatr(k)) > heatr(j)) then ! layer absorbs all
2275 heatr(k) = heatr(k) + heatr(j)
2276 stc(k) = tfrz + heatr(k)*fact(k)
2277 heatr(j) = 0.0
2278 else
2279 heatr(j) = heatr(j) + heatr(k)
2280 heatr(k) = 0.0
2281 stc(k) = tfrz
2282 end if
2283 end if
2284 end do
2285 stc(j) = tfrz + heatr(j)*fact(j)
2286 end if
2287 end do
2288 end if
2289
2290! now remove excess cold by increasing temperature of layers (may not be necessary with above loop)
2291
2292 if (any(stc(1:4) > tfrz) .and. any(stc(1:4) < tfrz)) then
2293 do j = 1,nsoil
2294 if ( stc(j) < tfrz ) then
2295 heatr(j) = (stc(j)-tfrz)/fact(j)
2296 do k = 1,nsoil
2297 if (j .ne. k .and. stc(k) > tfrz .and. heatr(j) < -0.1) then
2298 heatr(k) = (stc(k)-tfrz)/fact(k)
2299 if (heatr(k) > abs(heatr(j))) then ! layer absorbs all
2300 heatr(k) = heatr(k) + heatr(j)
2301 stc(k) = tfrz + heatr(k)*fact(k)
2302 heatr(j) = 0.0
2303 else
2304 heatr(j) = heatr(j) + heatr(k)
2305 heatr(k) = 0.0
2306 stc(k) = tfrz
2307 end if
2308 end if
2309 end do
2310 stc(j) = tfrz + heatr(j)*fact(j)
2311 end if
2312 end do
2313 end if
2314
2315! now remove excess heat by melting ice
2316
2317 if (any(stc(1:4) > tfrz) .and. any(mice(1:4) > 0.)) then
2318 do j = 1,nsoil
2319 if ( stc(j) > tfrz ) then
2320 heatr(j) = (stc(j)-tfrz)/fact(j)
2321 xm(j) = heatr(j)*dt/hfus
2322 do k = 1,nsoil
2323 if (j .ne. k .and. mice(k) > 0. .and. xm(j) > 0.1) then
2324 if (mice(k) > xm(j)) then ! layer absorbs all
2325 mice(k) = mice(k) - xm(j)
2326 xmf = xmf + hfus * xm(j)/dt
2327 stc(k) = tfrz
2328 xm(j) = 0.0
2329 else
2330 xm(j) = xm(j) - mice(k)
2331 xmf = xmf + hfus * mice(k)/dt
2332 mice(k) = 0.0
2333 stc(k) = tfrz
2334 end if
2335 mliq(k) = max(0.,wmass0(k)-mice(k))
2336 end if
2337 end do
2338 heatr(j) = xm(j)*hfus/dt
2339 stc(j) = tfrz + heatr(j)*fact(j)
2340 end if
2341 end do
2342 end if
2343
2344! now remove excess cold by freezing liquid of layers (may not be necessary with above loop)
2345
2346 if (any(stc(1:4) < tfrz) .and. any(mliq(1:4) > 0.)) then
2347 do j = 1,nsoil
2348 if ( stc(j) < tfrz ) then
2349 heatr(j) = (stc(j)-tfrz)/fact(j)
2350 xm(j) = heatr(j)*dt/hfus
2351 do k = 1,nsoil
2352 if (j .ne. k .and. mliq(k) > 0. .and. xm(j) < -0.1) then
2353 if (mliq(k) > abs(xm(j))) then ! layer absorbs all
2354 mice(k) = mice(k) - xm(j)
2355 xmf = xmf + hfus * xm(j)/dt
2356 stc(k) = tfrz
2357 xm(j) = 0.0
2358 else
2359 xm(j) = xm(j) + mliq(k)
2360 xmf = xmf - hfus * mliq(k)/dt
2361 mice(k) = wmass0(k)
2362 stc(k) = tfrz
2363 end if
2364 mliq(k) = max(0.,wmass0(k)-mice(k))
2365 end if
2366 end do
2367 heatr(j) = xm(j)*hfus/dt
2368 stc(j) = tfrz + heatr(j)*fact(j)
2369 end if
2370 end do
2371 end if
2372
2373end if ! opt_gla == 1
2374
2375 do j = isnow+1,0 ! snow
2376 snliq(j) = mliq(j)
2377 snice(j) = mice(j)
2378 end do
2379
2380 do j = 1, nsoil ! soil
2381 if(opt_gla == 1) then
2382 sh2o(j) = mliq(j) / (1000. * dzsnso(j))
2383 sh2o(j) = max(0.0,min(1.0,sh2o(j)))
2384! smc(j) = (mliq(j) + mice(j)) / (1000. * dzsnso(j))
2385 elseif(opt_gla == 2) then
2386 sh2o(j) = 0.0 ! ice, assume all frozen...forever
2387 end if
2388 smc(j) = 1.0
2389 end do
2390
2391 end subroutine phasechange_glacier
2392! ==================================================================================================
2394 subroutine water_glacier (nsnow ,nsoil ,imelt ,dt ,prcp ,sfctmp , & !in
2395 qvap ,qdew ,ficeold ,zsoil , & !in
2396 isnow ,snowh ,sneqv ,snice ,snliq ,stc , & !inout
2397 dzsnso ,sh2o ,sice ,ponding ,zsnso ,fsh , & !inout
2398 runsrf ,runsub ,qsnow ,ponding1 ,ponding2 ,qsnbot , & !out
2399 fpice ,esnow) !out
2400! ----------------------------------------------------------------------
2401! code history:
2402! initial code: guo-yue niu, oct. 2007
2403! ----------------------------------------------------------------------
2404 implicit none
2405! ----------------------------------------------------------------------
2406! input
2407 integer, intent(in) :: nsnow
2408 integer, intent(in) :: nsoil
2409 integer, dimension(-nsnow+1:0) , intent(in) :: imelt
2410 real (kind=kind_phys), intent(in) :: dt
2411 real (kind=kind_phys), intent(in) :: prcp
2412 real (kind=kind_phys), intent(in) :: sfctmp
2413 real (kind=kind_phys), intent(inout) :: qvap
2414 real (kind=kind_phys), intent(inout) :: qdew
2415 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(in) :: ficeold
2416 real (kind=kind_phys), dimension( 1:nsoil), intent(in) :: zsoil
2417
2418! input/output
2419 integer, intent(inout) :: isnow
2420 real (kind=kind_phys), intent(inout) :: snowh
2421 real (kind=kind_phys), intent(inout) :: sneqv
2422 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(inout) :: snice
2423 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(inout) :: snliq
2424 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: stc
2425 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: dzsnso
2426 real (kind=kind_phys), dimension( 1:nsoil), intent(inout) :: sh2o
2427 real (kind=kind_phys), dimension( 1:nsoil), intent(inout) :: sice
2428 real (kind=kind_phys) , intent(inout) :: ponding
2429 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: zsnso
2430 real (kind=kind_phys) , intent(inout) :: fsh
2431
2432! output
2433 real (kind=kind_phys), intent(out) :: runsrf
2434 real (kind=kind_phys), intent(out) :: runsub
2435 real (kind=kind_phys), intent(out) :: qsnow
2436 real (kind=kind_phys), intent(out) :: ponding1
2437 real (kind=kind_phys), intent(out) :: ponding2
2438 real (kind=kind_phys), intent(out) :: qsnbot
2439 real (kind=kind_phys), intent(out) :: fpice
2440 real (kind=kind_phys), intent(out) :: esnow
2441
2442! local
2443 real (kind=kind_phys) :: qrain
2444 real (kind=kind_phys) :: qseva
2445 real (kind=kind_phys) :: qsdew
2446 real (kind=kind_phys) :: qsnfro
2447 real (kind=kind_phys) :: qsnsub
2448 real (kind=kind_phys) :: snowhin
2449 real (kind=kind_phys) :: snoflow
2450 real (kind=kind_phys) :: bdfall
2451 real (kind=kind_phys) :: replace
2452 real (kind=kind_phys), dimension( 1:nsoil) :: sice_save
2453 real (kind=kind_phys), dimension( 1:nsoil) :: sh2o_save
2454 integer :: ilev
2455
2456
2457! ----------------------------------------------------------------------
2458! initialize
2459
2460 snoflow = 0.
2461 runsub = 0.
2462 runsrf = 0.
2463 sice_save = sice
2464 sh2o_save = sh2o
2465
2466! --------------------------------------------------------------------
2467! partition precipitation into rain and snow (from canwater)
2468
2469! jordan (1991)
2470
2471 if(opt_snf == 1 .or. opt_snf == 4) then
2472 if(sfctmp > tfrz+2.5)then
2473 fpice = 0.
2474 else
2475 if(sfctmp <= tfrz+0.5)then
2476 fpice = 1.0
2477 else if(sfctmp <= tfrz+2.)then
2478 fpice = 1.-(-54.632 + 0.2*sfctmp)
2479 else
2480 fpice = 0.6
2481 endif
2482 endif
2483 endif
2484
2485 if(opt_snf == 2) then
2486 if(sfctmp >= tfrz+2.2) then
2487 fpice = 0.
2488 else
2489 fpice = 1.0
2490 endif
2491 endif
2492
2493 if(opt_snf == 3) then
2494 if(sfctmp >= tfrz) then
2495 fpice = 0.
2496 else
2497 fpice = 1.0
2498 endif
2499 endif
2500! print*, 'fpice: ',fpice
2501
2502! hedstrom nr and jw pomeroy (1998), hydrol. processes, 12, 1611-1625
2503! fresh snow density
2504
2505 bdfall = min(120.,67.92+51.25*exp((sfctmp-tfrz)/2.59)) !mb: change to min v3.7
2506
2507 qrain = prcp * (1.-fpice)
2508 qsnow = prcp * fpice
2509 snowhin = qsnow/bdfall
2510! print *, 'qrain, qsnow',qrain,qsnow,qrain*dt,qsnow*dt
2511
2512! sublimation, frost, evaporation, and dew
2513
2514 qsnsub = qvap ! send total sublimation/frost to snowwater and deal with it there
2515 qsnfro = qdew
2516 esnow = qsnsub*2.83e+6
2517
2518 call snowwater_glacier (nsnow ,nsoil ,imelt ,dt ,sfctmp , & !in
2519 snowhin ,qsnow ,qsnfro ,qsnsub ,qrain , & !in
2520 ficeold ,zsoil , & !in
2521 isnow ,snowh ,sneqv ,snice ,snliq , & !inout
2522 sh2o ,sice ,stc ,dzsnso ,zsnso , & !inout
2523 fsh , & !inout
2524 qsnbot ,snoflow ,ponding1 ,ponding2) !out
2525
2526 !ponding: melting water from snow when there is no layer
2527
2528 runsrf = (ponding+ponding1+ponding2)/dt
2529
2530 if(isnow == 0) then
2531 runsrf = runsrf + qsnbot + qrain
2532 else
2533 runsrf = runsrf + qsnbot
2534 endif
2535
2536
2537 if(opt_gla == 1) then
2538 replace = 0.0
2539 do ilev = 1,nsoil
2540 replace = replace + dzsnso(ilev)*(sice(ilev) - sice_save(ilev) + sh2o(ilev) - sh2o_save(ilev))
2541 end do
2542 replace = replace * 1000.0 / dt ! convert to [mm/s]
2543
2544 sice = min(1.0,sice_save)
2545 elseif(opt_gla == 2) then
2546 sice = 1.0
2547 end if
2548 sh2o = 1.0 - sice
2549
2550 ! use runsub as a water balancer, snoflow is snow that disappears, replace is
2551 ! water from below that replaces glacier loss
2552
2553 if(opt_gla == 1) then
2554 runsub = snoflow + replace
2555 elseif(opt_gla == 2) then
2556 runsub = snoflow
2557 qvap = qsnsub
2558 qdew = qsnfro
2559 end if
2560
2561 end subroutine water_glacier
2562! ==================================================================================================
2563! ----------------------------------------------------------------------
2565 subroutine snowwater_glacier (nsnow ,nsoil ,imelt ,dt ,sfctmp , & !in
2566 snowhin ,qsnow ,qsnfro ,qsnsub ,qrain , & !in
2567 ficeold ,zsoil , & !in
2568 isnow ,snowh ,sneqv ,snice ,snliq , & !inout
2569 sh2o ,sice ,stc ,dzsnso ,zsnso , & !inout
2570 fsh , & !inout
2571 qsnbot ,snoflow ,ponding1 ,ponding2) !out
2572! ----------------------------------------------------------------------
2573 implicit none
2574! ----------------------------------------------------------------------
2575! input
2576 integer, intent(in) :: nsnow
2577 integer, intent(in) :: nsoil
2578 integer, dimension(-nsnow+1:0) , intent(in) :: imelt
2579 real (kind=kind_phys), intent(in) :: dt
2580 real (kind=kind_phys), intent(in) :: sfctmp
2581 real (kind=kind_phys), intent(in) :: snowhin
2582 real (kind=kind_phys), intent(in) :: qsnow
2583 real (kind=kind_phys), intent(inout) :: qsnfro
2584 real (kind=kind_phys), intent(inout) :: qsnsub
2585 real (kind=kind_phys), intent(in) :: qrain
2586 real (kind=kind_phys), dimension(-nsnow+1:0) , intent(in) :: ficeold
2587 real (kind=kind_phys), dimension( 1:nsoil), intent(in) :: zsoil
2588
2589! input & output
2590 integer, intent(inout) :: isnow
2591 real (kind=kind_phys), intent(inout) :: snowh
2592 real (kind=kind_phys), intent(inout) :: sneqv
2593 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(inout) :: snice
2594 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(inout) :: snliq
2595 real (kind=kind_phys), dimension( 1:nsoil), intent(inout) :: sh2o
2596 real (kind=kind_phys), dimension( 1:nsoil), intent(inout) :: sice
2597 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: stc
2598 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: dzsnso
2599 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: zsnso
2600 real (kind=kind_phys), intent(inout) :: fsh
2601
2602! output
2603 real (kind=kind_phys), intent(out) :: qsnbot
2604 real (kind=kind_phys), intent(out) :: snoflow
2605 real (kind=kind_phys), intent(out) :: ponding1
2606 real (kind=kind_phys), intent(out) :: ponding2
2607
2608! local
2609 integer :: iz
2610 real (kind=kind_phys) :: bdsnow
2611 real (kind=kind_phys),parameter :: mwd = 100.
2612! ----------------------------------------------------------------------
2613 snoflow = 0.0
2614 ponding1 = 0.0
2615 ponding2 = 0.0
2616
2617 call snowfall_glacier (nsoil ,nsnow ,dt ,qsnow ,snowhin, & !in
2618 sfctmp , & !in
2619 isnow ,snowh ,dzsnso ,stc ,snice , & !inout
2620 snliq ,sneqv ) !inout
2621
2622 if(isnow < 0) then !when more than one layer
2623 call compact_glacier (nsnow ,nsoil ,dt ,stc ,snice , & !in
2624 snliq ,imelt ,ficeold, & !in
2625 isnow ,dzsnso ) !inout
2626
2627 call combine_glacier (nsnow ,nsoil , & !in
2628 isnow ,sh2o ,stc ,snice ,snliq , & !inout
2629 dzsnso ,sice ,snowh ,sneqv , & !inout
2630 ponding1 ,ponding2) !out
2631
2632 call divide_glacier (nsnow ,nsoil , & !in
2633 isnow ,stc ,snice ,snliq ,dzsnso ) !inout
2634 end if
2635
2636!set empty snow layers to zero
2637
2638 do iz = -nsnow+1, isnow
2639 snice(iz) = 0.
2640 snliq(iz) = 0.
2641 stc(iz) = 0.
2642 dzsnso(iz)= 0.
2643 zsnso(iz) = 0.
2644 enddo
2645
2646 call snowh2o_glacier (nsnow ,nsoil ,dt ,qsnfro ,qsnsub , & !in
2647 qrain , & !in
2648 isnow ,dzsnso ,snowh ,sneqv ,snice , & !inout
2649 snliq ,sh2o ,sice ,stc , & !inout
2650 ponding1 ,ponding2 ,fsh , & !inout
2651 qsnbot ) !out
2652
2653!to obtain equilibrium state of snow in glacier region
2654
2655 if(sneqv > mwd .and. isnow /= 0) then ! 100 mm -> maximum water depth
2656 bdsnow = snice(0) / dzsnso(0)
2657 snoflow = (sneqv - mwd)
2658 snice(0) = snice(0) - snoflow
2659 dzsnso(0) = dzsnso(0) - snoflow/bdsnow
2660 snoflow = snoflow / dt
2661 end if
2662
2663! sum up snow mass for layered snow
2664
2665 if(isnow /= 0) then
2666 sneqv = 0.
2667 snowh = 0.
2668 do iz = isnow+1,0
2669 sneqv = sneqv + snice(iz) + snliq(iz)
2670 snowh = snowh + dzsnso(iz)
2671 enddo
2672 end if
2673
2674! reset zsnso and layer thinkness dzsnso
2675
2676 do iz = isnow+1, 0
2677 dzsnso(iz) = -dzsnso(iz)
2678 end do
2679
2680 dzsnso(1) = zsoil(1)
2681 do iz = 2,nsoil
2682 dzsnso(iz) = (zsoil(iz) - zsoil(iz-1))
2683 end do
2684
2685 zsnso(isnow+1) = dzsnso(isnow+1)
2686 do iz = isnow+2 ,nsoil
2687 zsnso(iz) = zsnso(iz-1) + dzsnso(iz)
2688 enddo
2689
2690 do iz = isnow+1 ,nsoil
2691 dzsnso(iz) = -dzsnso(iz)
2692 end do
2693
2694 end subroutine snowwater_glacier
2695! ==================================================================================================
2697 subroutine snowfall_glacier (nsoil ,nsnow ,dt ,qsnow ,snowhin , & !in
2698 sfctmp , & !in
2699 isnow ,snowh ,dzsnso ,stc ,snice , & !inout
2700 snliq ,sneqv ) !inout
2701! ----------------------------------------------------------------------
2704! ----------------------------------------------------------------------
2705 implicit none
2706! ----------------------------------------------------------------------
2707! input
2708
2709 integer, intent(in) :: nsoil
2710 integer, intent(in) :: nsnow
2711 real (kind=kind_phys), intent(in) :: dt
2712 real (kind=kind_phys), intent(in) :: qsnow
2713 real (kind=kind_phys), intent(in) :: snowhin
2714 real (kind=kind_phys), intent(in) :: sfctmp
2715
2716! input and output
2717
2718 integer, intent(inout) :: isnow
2719 real (kind=kind_phys), intent(inout) :: snowh
2720 real (kind=kind_phys), intent(inout) :: sneqv
2721 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: dzsnso
2722 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: stc
2723 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(inout) :: snice
2724 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(inout) :: snliq
2725
2726! local
2727
2728 integer :: newnode
2729! ----------------------------------------------------------------------
2730 newnode = 0
2731
2732! shallow snow / no layer
2733
2734 if(isnow == 0 .and. qsnow > 0.) then
2735 snowh = snowh + snowhin * dt
2736 sneqv = sneqv + qsnow * dt
2737 end if
2738
2739! creating a new layer
2740
2741 if(isnow == 0 .and. qsnow>0. .and. snowh >= 0.05) then
2742 isnow = -1
2743 newnode = 1
2744 dzsnso(0)= snowh
2745 snowh = 0.
2746 stc(0) = min(273.16, sfctmp) ! temporary setup
2747 snice(0) = sneqv
2748 snliq(0) = 0.
2749 end if
2750
2751! snow with layers
2752
2753 if(isnow < 0 .and. newnode == 0 .and. qsnow > 0.) then
2754 snice(isnow+1) = snice(isnow+1) + qsnow * dt
2755 dzsnso(isnow+1) = dzsnso(isnow+1) + snowhin * dt
2756 endif
2757
2758! ----------------------------------------------------------------------
2759 end subroutine snowfall_glacier
2760! ==================================================================================================
2761! ----------------------------------------------------------------------
2763 subroutine compact_glacier (nsnow ,nsoil ,dt ,stc ,snice , & !in
2764 snliq ,imelt ,ficeold, & !in
2765 isnow ,dzsnso ) !inout
2766! ----------------------------------------------------------------------
2767! ----------------------------------------------------------------------
2768 implicit none
2769! ----------------------------------------------------------------------
2770! input
2771 integer, intent(in) :: nsoil
2772 integer, intent(in) :: nsnow
2773 integer, dimension(-nsnow+1:0) , intent(in) :: imelt
2774 real (kind=kind_phys), intent(in) :: dt
2775 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(in) :: stc
2776 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(in) :: snice
2777 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(in) :: snliq
2778 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(in) :: ficeold
2779
2780! input and output
2781 integer, intent(inout) :: isnow
2782 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: dzsnso
2783
2784! local
2785 real (kind=kind_phys), parameter :: c2 = 21.e-3
2786 real (kind=kind_phys), parameter :: c3 = 2.5e-6
2787 real (kind=kind_phys), parameter :: c4 = 0.04
2788 real (kind=kind_phys), parameter :: c5 = 2.0
2789 real (kind=kind_phys), parameter :: dm = 100.0
2790 real (kind=kind_phys), parameter :: eta0 = 1.8e+6
2791 !according to anderson, it is between 0.52e6~1.38e6
2792 real (kind=kind_phys) :: burden
2793 real (kind=kind_phys) :: ddz1
2794 real (kind=kind_phys) :: ddz2
2795 real (kind=kind_phys) :: ddz3
2796 real (kind=kind_phys) :: dexpf
2797 real (kind=kind_phys) :: td
2798 real (kind=kind_phys) :: pdzdtc
2799 real (kind=kind_phys) :: void
2800 real (kind=kind_phys) :: wx
2801 real (kind=kind_phys) :: bi
2802 real (kind=kind_phys), dimension(-nsnow+1:0) :: fice
2803
2804 integer :: j
2805
2806! ----------------------------------------------------------------------
2807 burden = 0.0
2808
2809 do j = isnow+1, 0
2810
2811 wx = snice(j) + snliq(j)
2812 fice(j) = snice(j) / wx
2813 void = 1. - (snice(j)/denice + snliq(j)/denh2o) / dzsnso(j)
2814
2815 ! allow compaction only for non-saturated node and higher ice lens node.
2816 if (void > 0.001 .and. snice(j) > 0.1) then
2817 bi = snice(j) / dzsnso(j)
2818 td = max(0.,tfrz-stc(j))
2819 dexpf = exp(-c4*td)
2820
2821 ! settling as a result of destructive metamorphism
2822
2823 ddz1 = -c3*dexpf
2824
2825 if (bi > dm) ddz1 = ddz1*exp(-46.0e-3*(bi-dm))
2826
2827 ! liquid water term
2828
2829 if (snliq(j) > 0.01*dzsnso(j)) ddz1=ddz1*c5
2830
2831 ! compaction due to overburden
2832
2833 ddz2 = -(burden+0.5*wx)*exp(-0.08*td-c2*bi)/eta0 ! 0.5*wx -> self-burden
2834
2835 ! compaction occurring during melt
2836
2837 if (imelt(j) == 1) then
2838 ddz3 = max(0.,(ficeold(j) - fice(j))/max(1.e-6,ficeold(j)))
2839 ddz3 = - ddz3/dt ! sometimes too large
2840 else
2841 ddz3 = 0.
2842 end if
2843
2844 ! time rate of fractional change in dz (units of s-1)
2845
2846 pdzdtc = (ddz1 + ddz2 + ddz3)*dt
2847 pdzdtc = max(-0.5,pdzdtc)
2848
2849 ! the change in dz due to compaction
2850
2851 dzsnso(j) = dzsnso(j)*(1.+pdzdtc)
2852 dzsnso(j) = min(max(dzsnso(j),(snliq(j)+snice(j))/500.0),(snliq(j)+snice(j))/50.0) ! limit adjustment to a reasonable density
2853 end if
2854
2855 ! pressure of overlying snow
2856
2857 burden = burden + wx
2858
2859 end do
2860
2861 end subroutine compact_glacier
2862! ==================================================================================================
2864 subroutine combine_glacier (nsnow ,nsoil , & !in
2865 isnow ,sh2o ,stc ,snice ,snliq , & !inout
2866 dzsnso ,sice ,snowh ,sneqv , & !inout
2867 ponding1 ,ponding2) !inout
2868! ----------------------------------------------------------------------
2869 implicit none
2870! ----------------------------------------------------------------------
2871! input
2872
2873 integer, intent(in) :: nsnow
2874 integer, intent(in) :: nsoil
2875
2876! input and output
2877
2878 integer, intent(inout) :: isnow
2879 real (kind=kind_phys), dimension( 1:nsoil), intent(inout) :: sh2o
2880 real (kind=kind_phys), dimension( 1:nsoil), intent(inout) :: sice
2881 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: stc
2882 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(inout) :: snice
2883 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(inout) :: snliq
2884 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: dzsnso
2885 real (kind=kind_phys), intent(inout) :: sneqv
2886 real (kind=kind_phys), intent(inout) :: snowh
2887 real (kind=kind_phys), intent(inout) :: ponding1
2888 real (kind=kind_phys), intent(inout) :: ponding2
2889
2890! local variables:
2891
2892 integer :: i,j,k,l
2893 integer :: isnow_old
2894 integer :: mssi
2895 integer :: neibor
2896 real (kind=kind_phys) :: zwice
2897 real (kind=kind_phys) :: zwliq
2898 real (kind=kind_phys) :: dzmin(3)
2899 data dzmin /0.045, 0.05, 0.2/
2900! data dzmin /0.025, 0.025, 0.1/ ! mb: change limit
2901!-----------------------------------------------------------------------
2902
2903 isnow_old = isnow
2904
2905 do j = isnow_old+1,0
2906 if (snice(j) <= .1) then
2907 if(j /= 0) then
2908 snliq(j+1) = snliq(j+1) + snliq(j)
2909 snice(j+1) = snice(j+1) + snice(j)
2910 else
2911 if (isnow_old < -1) then
2912 snliq(j-1) = snliq(j-1) + snliq(j)
2913 snice(j-1) = snice(j-1) + snice(j)
2914 else
2915 ponding1 = ponding1 + snliq(j) ! isnow will get set to zero below
2916 sneqv = snice(j) ! ponding will get added to ponding from
2917 snowh = dzsnso(j) ! phasechange which should be zero here
2918 snliq(j) = 0.0 ! because there it was only calculated
2919 snice(j) = 0.0 ! for thin snow
2920 dzsnso(j) = 0.0
2921 endif
2922! sh2o(1) = sh2o(1)+snliq(j)/(dzsnso(1)*1000.)
2923! sice(1) = sice(1)+snice(j)/(dzsnso(1)*1000.)
2924 endif
2925
2926 ! shift all elements above this down by one.
2927 if (j > isnow+1 .and. isnow < -1) then
2928 do i = j, isnow+2, -1
2929 stc(i) = stc(i-1)
2930 snliq(i) = snliq(i-1)
2931 snice(i) = snice(i-1)
2932 dzsnso(i)= dzsnso(i-1)
2933 end do
2934 end if
2935 isnow = isnow + 1
2936 end if
2937 end do
2938
2939! to conserve water in case of too large surface sublimation
2940
2941 if(sice(1) < 0.) then
2942 sh2o(1) = sh2o(1) + sice(1)
2943 sice(1) = 0.
2944 end if
2945
2946 if(isnow ==0) return ! mb: get out if no longer multi-layer
2947
2948 sneqv = 0.
2949 snowh = 0.
2950 zwice = 0.
2951 zwliq = 0.
2952
2953 do j = isnow+1,0
2954 sneqv = sneqv + snice(j) + snliq(j)
2955 snowh = snowh + dzsnso(j)
2956 zwice = zwice + snice(j)
2957 zwliq = zwliq + snliq(j)
2958 end do
2959
2960! check the snow depth - all snow gone
2961! the liquid water assumes ponding on soil surface.
2962
2963! if (snowh < 0.025 .and. isnow < 0 ) then ! mb: change limit
2964 if (snowh < 0.05 .and. isnow < 0 ) then
2965 isnow = 0
2966 sneqv = zwice
2967 ponding2 = ponding2 + zwliq ! limit of isnow < 0 means input ponding
2968 if(sneqv <= 0.) snowh = 0. ! should be zero; see above
2969 end if
2970
2971! if (snowh < 0.05 ) then
2972! isnow = 0
2973! sneqv = zwice
2974! sh2o(1) = sh2o(1) + zwliq / (dzsnso(1) * 1000.)
2975! if(sneqv <= 0.) snowh = 0.
2976! end if
2977
2978! check the snow depth - snow layers combined
2979
2980 if (isnow < -1) then
2981
2982 isnow_old = isnow
2983 mssi = 1
2984
2985 do i = isnow_old+1,0
2986 if (dzsnso(i) < dzmin(mssi)) then
2987
2988 if (i == isnow+1) then
2989 neibor = i + 1
2990 else if (i == 0) then
2991 neibor = i - 1
2992 else
2993 neibor = i + 1
2994 if ((dzsnso(i-1)+dzsnso(i)) < (dzsnso(i+1)+dzsnso(i))) neibor = i-1
2995 end if
2996
2997 ! node l and j are combined and stored as node j.
2998 if (neibor > i) then
2999 j = neibor
3000 l = i
3001 else
3002 j = i
3003 l = neibor
3004 end if
3005
3006 call combo_glacier (dzsnso(j), snliq(j), snice(j), &
3007 stc(j), dzsnso(l), snliq(l), snice(l), stc(l) )
3008
3009 ! now shift all elements above this down one.
3010 if (j-1 > isnow+1) then
3011 do k = j-1, isnow+2, -1
3012 stc(k) = stc(k-1)
3013 snice(k) = snice(k-1)
3014 snliq(k) = snliq(k-1)
3015 dzsnso(k) = dzsnso(k-1)
3016 end do
3017 end if
3018
3019 ! decrease the number of snow layers
3020 isnow = isnow + 1
3021 if (isnow >= -1) exit
3022 else
3023
3024 ! the layer thickness is greater than the prescribed minimum value
3025 mssi = mssi + 1
3026
3027 end if
3028 end do
3029
3030 end if
3031
3032 end subroutine combine_glacier
3033! ==================================================================================================
3034
3035! ----------------------------------------------------------------------
3037 subroutine combo_glacier(dz, wliq, wice, t, dz2, wliq2, wice2, t2)
3038! ----------------------------------------------------------------------
3039 implicit none
3040! ----------------------------------------------------------------------
3041
3042! ----------------------------------------------------------------------s
3043! input
3044
3045 real (kind=kind_phys), intent(in) :: dz2
3046 real (kind=kind_phys), intent(in) :: wliq2
3047 real (kind=kind_phys), intent(in) :: wice2
3048 real (kind=kind_phys), intent(in) :: t2
3049 real (kind=kind_phys), intent(inout) :: dz
3050 real (kind=kind_phys), intent(inout) :: wliq
3051 real (kind=kind_phys), intent(inout) :: wice
3052 real (kind=kind_phys), intent(inout) :: t
3053
3054! local
3055
3056 real (kind=kind_phys) :: dzc
3057 real (kind=kind_phys) :: wliqc
3058 real (kind=kind_phys) :: wicec
3059 real (kind=kind_phys) :: tc
3060 real (kind=kind_phys) :: h
3061 real (kind=kind_phys) :: h2
3062 real (kind=kind_phys) :: hc
3063
3064!-----------------------------------------------------------------------
3065
3066 dzc = dz+dz2
3067 wicec = (wice+wice2)
3068 wliqc = (wliq+wliq2)
3069 h = (cice*wice+cwat*wliq) * (t-tfrz)+hfus*wliq
3070 h2= (cice*wice2+cwat*wliq2) * (t2-tfrz)+hfus*wliq2
3071
3072 hc = h + h2
3073 if(hc < 0.)then
3074 tc = tfrz + hc/(cice*wicec + cwat*wliqc)
3075 else if (hc.le.hfus*wliqc) then
3076 tc = tfrz
3077 else
3078 tc = tfrz + (hc - hfus*wliqc) / (cice*wicec + cwat*wliqc)
3079 end if
3080
3081 dz = dzc
3082 wice = wicec
3083 wliq = wliqc
3084 t = tc
3085
3086 end subroutine combo_glacier
3087! ==================================================================================================
3089 subroutine divide_glacier (nsnow ,nsoil , & !in
3090 isnow ,stc ,snice ,snliq ,dzsnso ) !inout
3091! ----------------------------------------------------------------------
3092 implicit none
3093! ----------------------------------------------------------------------
3094! input
3095
3096 integer, intent(in) :: nsnow
3097 integer, intent(in) :: nsoil
3098
3099! input and output
3100
3101 integer , intent(inout) :: isnow
3102 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: stc
3103 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(inout) :: snice
3104 real (kind=kind_phys), dimension(-nsnow+1: 0), intent(inout) :: snliq
3105 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: dzsnso
3106
3107! local variables:
3108
3109 integer :: j
3110 integer :: msno
3111 real (kind=kind_phys) :: drr
3112 real (kind=kind_phys), dimension( 1:nsnow) :: dz
3113 real (kind=kind_phys), dimension( 1:nsnow) :: swice
3114 real (kind=kind_phys), dimension( 1:nsnow) :: swliq
3115 real (kind=kind_phys), dimension( 1:nsnow) :: tsno
3116 real (kind=kind_phys) :: zwice
3117 real (kind=kind_phys) :: zwliq
3118 real (kind=kind_phys) :: propor
3119 real (kind=kind_phys) :: dtdz
3120! ----------------------------------------------------------------------
3121
3122 do j = 1,nsnow
3123 if (j <= abs(isnow)) then
3124 dz(j) = dzsnso(j+isnow)
3125 swice(j) = snice(j+isnow)
3126 swliq(j) = snliq(j+isnow)
3127 tsno(j) = stc(j+isnow)
3128 end if
3129 end do
3130
3131 msno = abs(isnow)
3132
3133 if (msno == 1) then
3134 ! specify a new snow layer
3135 if (dz(1) > 0.05) then
3136 msno = 2
3137 dz(1) = dz(1)/2.
3138 swice(1) = swice(1)/2.
3139 swliq(1) = swliq(1)/2.
3140 dz(2) = dz(1)
3141 swice(2) = swice(1)
3142 swliq(2) = swliq(1)
3143 tsno(2) = tsno(1)
3144 end if
3145 end if
3146
3147 if (msno > 1) then
3148 if (dz(1) > 0.05) then
3149 drr = dz(1) - 0.05
3150 propor = drr/dz(1)
3151 zwice = propor*swice(1)
3152 zwliq = propor*swliq(1)
3153 propor = 0.05/dz(1)
3154 swice(1) = propor*swice(1)
3155 swliq(1) = propor*swliq(1)
3156 dz(1) = 0.05
3157
3158 call combo_glacier (dz(2), swliq(2), swice(2), tsno(2), drr, &
3159 zwliq, zwice, tsno(1))
3160
3161 ! subdivide a new layer
3162! if (msno <= 2 .and. dz(2) > 0.20) then ! mb: change limit
3163 if (msno <= 2 .and. dz(2) > 0.10) then
3164 msno = 3
3165 dtdz = (tsno(1) - tsno(2))/((dz(1)+dz(2))/2.)
3166 dz(2) = dz(2)/2.
3167 swice(2) = swice(2)/2.
3168 swliq(2) = swliq(2)/2.
3169 dz(3) = dz(2)
3170 swice(3) = swice(2)
3171 swliq(3) = swliq(2)
3172 tsno(3) = tsno(2) - dtdz*dz(2)/2.
3173 if (tsno(3) >= tfrz) then
3174 tsno(3) = tsno(2)
3175 else
3176 tsno(2) = tsno(2) + dtdz*dz(2)/2.
3177 endif
3178
3179 end if
3180 end if
3181 end if
3182
3183 if (msno > 2) then
3184 if (dz(2) > 0.2) then
3185 drr = dz(2) - 0.2
3186 propor = drr/dz(2)
3187 zwice = propor*swice(2)
3188 zwliq = propor*swliq(2)
3189 propor = 0.2/dz(2)
3190 swice(2) = propor*swice(2)
3191 swliq(2) = propor*swliq(2)
3192 dz(2) = 0.2
3193 call combo_glacier (dz(3), swliq(3), swice(3), tsno(3), drr, &
3194 zwliq, zwice, tsno(2))
3195 end if
3196 end if
3197
3198 isnow = -msno
3199
3200 do j = isnow+1,0
3201 dzsnso(j) = dz(j-isnow)
3202 snice(j) = swice(j-isnow)
3203 snliq(j) = swliq(j-isnow)
3204 stc(j) = tsno(j-isnow)
3205 end do
3206
3207
3208! do j = isnow+1,nsoil
3209! write(*,'(i5,7f10.3)') j, dzsnso(j), snice(j), snliq(j),stc(j)
3210! end do
3211
3212 end subroutine divide_glacier
3213! ==================================================================================================
3215 subroutine snowh2o_glacier (nsnow ,nsoil ,dt ,qsnfro ,qsnsub , & !in
3216 qrain , & !in
3217 isnow ,dzsnso ,snowh ,sneqv ,snice , & !inout
3218 snliq ,sh2o ,sice ,stc , & !inout
3219 ponding1 ,ponding2 ,fsh , & !inout
3220 qsnbot ) !out
3221! ----------------------------------------------------------------------
3224! ----------------------------------------------------------------------
3225 implicit none
3226! ----------------------------------------------------------------------
3227! input
3228
3229 integer, intent(in) :: nsnow
3230 integer, intent(in) :: nsoil
3231 real (kind=kind_phys), intent(in) :: dt
3232 real (kind=kind_phys), intent(inout) :: qsnfro
3233 real (kind=kind_phys), intent(inout) :: qsnsub
3234 real (kind=kind_phys), intent(in) :: qrain
3235
3236! output
3237
3238 real (kind=kind_phys), intent(out) :: qsnbot
3239
3240! input and output
3241
3242 integer, intent(inout) :: isnow
3243 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: dzsnso
3244 real (kind=kind_phys), intent(inout) :: snowh
3245 real (kind=kind_phys), intent(inout) :: sneqv
3246 real (kind=kind_phys), dimension(-nsnow+1:0), intent(inout) :: snice
3247 real (kind=kind_phys), dimension(-nsnow+1:0), intent(inout) :: snliq
3248 real (kind=kind_phys), dimension( 1:nsoil), intent(inout) :: sh2o
3249 real (kind=kind_phys), dimension( 1:nsoil), intent(inout) :: sice
3250 real (kind=kind_phys), dimension(-nsnow+1:nsoil), intent(inout) :: stc
3251 real (kind=kind_phys), intent(inout) :: ponding1
3252 real (kind=kind_phys), intent(inout) :: ponding2
3253 real (kind=kind_phys), intent(inout) :: fsh
3254
3255! local variables:
3256
3257 integer :: j
3258 real (kind=kind_phys) :: qin
3259 real (kind=kind_phys) :: qout
3260 real (kind=kind_phys) :: wgdif
3261 real (kind=kind_phys), dimension(-nsnow+1:0) :: vol_liq
3262 real (kind=kind_phys), dimension(-nsnow+1:0) :: vol_ice
3263 real (kind=kind_phys), dimension(-nsnow+1:0) :: epore
3264 real (kind=kind_phys) :: propor, temp
3265! ----------------------------------------------------------------------
3266
3267!for the case when sneqv becomes '0' after 'combine'
3268
3269 if(sneqv == 0.) then
3270 if(opt_gla == 1) then
3271 sice(1) = sice(1) + (qsnfro-qsnsub)*dt/(dzsnso(1)*1000.)
3272 elseif(opt_gla == 2) then
3273 fsh = fsh - (qsnfro-qsnsub)*hsub
3274 qsnfro = 0.0
3275 qsnsub = 0.0
3276 end if
3277 end if
3278
3279! for shallow snow without a layer
3280! snow surface sublimation may be larger than existing snow mass. to conserve water,
3281! excessive sublimation is used to reduce soil water. smaller time steps would tend
3282! to aviod this problem.
3283
3284 if(isnow == 0 .and. sneqv > 0.) then
3285 if(opt_gla == 1) then
3286 temp = sneqv
3287 sneqv = sneqv - qsnsub*dt + qsnfro*dt
3288 propor = sneqv/temp
3289 snowh = max(0.,propor * snowh)
3290 elseif(opt_gla == 2) then
3291 fsh = fsh - (qsnfro-qsnsub)*hsub
3292 qsnfro = 0.0
3293 qsnsub = 0.0
3294 end if
3295
3296 if(sneqv < 0.) then
3297 sice(1) = sice(1) + sneqv/(dzsnso(1)*1000.)
3298 sneqv = 0.
3299 snowh = 0.
3300 end if
3301 if(sice(1) < 0.) then
3302 sh2o(1) = sh2o(1) + sice(1)
3303 sice(1) = 0.
3304 end if
3305 end if
3306
3307 if(snowh <= 1.e-8 .or. sneqv <= 1.e-6) then
3308 snowh = 0.0
3309 sneqv = 0.0
3310 end if
3311
3312! for deep snow
3313
3314 if ( isnow < 0 ) then !kwm added this if statement to prevent out-of-bounds array references
3315
3316 wgdif = snice(isnow+1) - qsnsub*dt + qsnfro*dt
3317 snice(isnow+1) = wgdif
3318 if (wgdif < 1.e-6 .and. isnow <0) then
3319 call combine_glacier (nsnow ,nsoil , & !in
3320 isnow ,sh2o ,stc ,snice ,snliq , & !inout
3321 dzsnso ,sice ,snowh ,sneqv , & !inout
3322 ponding1, ponding2 ) !inout
3323 endif
3324 !kwm: subroutine combine can change isnow to make it 0 again?
3325 if ( isnow < 0 ) then !kwm added this if statement to prevent out-of-bounds array references
3326 snliq(isnow+1) = snliq(isnow+1) + qrain * dt
3327 snliq(isnow+1) = max(0., snliq(isnow+1))
3328 endif
3329
3330 endif !kwm -- can the endif be moved toward the end of the subroutine (just set qsnbot=0)?
3331
3332! porosity and partial volume
3333
3334 !kwm looks to me like loop index / if test can be simplified.
3335
3336 do j = -nsnow+1, 0
3337 if (j >= isnow+1) then
3338 vol_ice(j) = min(1., snice(j)/(dzsnso(j)*denice))
3339 epore(j) = 1. - vol_ice(j)
3340 vol_liq(j) = min(epore(j),snliq(j)/(dzsnso(j)*denh2o))
3341 end if
3342 end do
3343
3344 qin = 0.
3345 qout = 0.
3346
3347 !kwm looks to me like loop index / if test can be simplified.
3348
3349 do j = -nsnow+1, 0
3350 if (j >= isnow+1) then
3351 snliq(j) = snliq(j) + qin
3352 if (j <= -1) then
3353 if (epore(j) < 0.05 .or. epore(j+1) < 0.05) then
3354 qout = 0.
3355 else
3356 qout = max(0.,(vol_liq(j)-ssi*epore(j))*dzsnso(j))
3357 qout = min(qout,(1.-vol_ice(j+1)-vol_liq(j+1))*dzsnso(j+1))
3358 end if
3359 else
3360 qout = max(0.,(vol_liq(j) - ssi*epore(j))*dzsnso(j))
3361 end if
3362 qout = qout*1000.
3363 snliq(j) = snliq(j) - qout
3364 qin = qout
3365 end if
3366 end do
3367
3368! liquid water from snow bottom to soil
3369
3370 qsnbot = qout / dt ! mm/s
3371
3372 end subroutine snowh2o_glacier
3373! ********************* end of water subroutines ******************************************
3374! ==================================================================================================
3376 subroutine error_glacier (iloc ,jloc ,swdown ,fsa ,fsr ,fira , &
3377 fsh ,fgev ,ssoil ,sag ,prcp ,edir , &
3378#ifdef CCPP
3379 runsrf ,runsub ,sneqv ,dt ,beg_wb, errmsg , errflg )
3380#else
3381 runsrf ,runsub ,sneqv ,dt ,beg_wb )
3382#endif
3383! --------------------------------------------------------------------------------------------------
3385! --------------------------------------------------------------------------------------------------
3386 implicit none
3387! --------------------------------------------------------------------------------------------------
3388! inputs
3389 integer , intent(in) :: iloc
3390 integer , intent(in) :: jloc
3391 real (kind=kind_phys) , intent(in) :: swdown
3392 real (kind=kind_phys) , intent(in) :: fsa
3393 real (kind=kind_phys) , intent(in) :: fsr
3394 real (kind=kind_phys) , intent(in) :: fira
3395 real (kind=kind_phys) , intent(in) :: fsh
3396 real (kind=kind_phys) , intent(in) :: fgev
3397 real (kind=kind_phys) , intent(in) :: ssoil
3398 real (kind=kind_phys) , intent(in) :: sag
3399
3400 real (kind=kind_phys) , intent(in) :: prcp
3401 real (kind=kind_phys) , intent(in) :: edir
3402 real (kind=kind_phys) , intent(in) :: runsrf
3403 real (kind=kind_phys) , intent(in) :: runsub
3404 real (kind=kind_phys) , intent(in) :: sneqv
3405 real (kind=kind_phys) , intent(in) :: dt
3406 real (kind=kind_phys) , intent(in) :: beg_wb
3407
3408#ifdef CCPP
3409 character(len=*) , intent(inout) :: errmsg
3410 integer , intent(inout) :: errflg
3411#endif
3412
3413 real (kind=kind_phys) :: end_wb
3414 real (kind=kind_phys) :: errwat
3415 real (kind=kind_phys) :: erreng
3416 real (kind=kind_phys) :: errsw
3417 character(len=256) :: message
3418! --------------------------------------------------------------------------------------------------
3419 errsw = swdown - (fsa + fsr)
3420 if (errsw > 0.01) then ! w/m2
3421 write(*,*) "sag =",sag
3422 write(*,*) "fsa =",fsa
3423 write(*,*) "fsr =",fsr
3424 write(message,*) 'errsw =',errsw
3425#ifdef CCPP
3426 errflg = 1
3427 errmsg = trim(message)//new_line('A')//"radiation budget problem in noahmp glacier"
3428 return
3429#else
3430 call wrf_message(trim(message))
3431 call wrf_error_fatal("radiation budget problem in noahmp glacier")
3432#endif
3433 end if
3434
3435 erreng = sag-(fira+fsh+fgev+ssoil)
3436 if(erreng > 0.01) then
3437 write(message,*) 'erreng =',erreng
3438#ifdef CCPP
3439 errmsg = trim(message)
3440#else
3441 call wrf_message(trim(message))
3442#endif
3443 write(message,'(i6,1x,i6,1x,5f10.4)')iloc,jloc,sag,fira,fsh,fgev,ssoil
3444#ifdef CCPP
3445 errflg = 1
3446 errmsg = trim(errmsg)//new_line('A')//"energy budget problem in noahmp glacier"
3447 return
3448#else
3449 call wrf_message(trim(message))
3450 call wrf_error_fatal("energy budget problem in noahmp glacier")
3451#endif
3452 end if
3453
3454 end_wb = sneqv
3455 errwat = end_wb-beg_wb-(prcp-edir-runsrf-runsub)*dt
3456
3457
3458 end subroutine error_glacier
3459! ==================================================================================================
3460
3463 subroutine noahmp_options_glacier(iopt_alb ,iopt_snf ,iopt_tbot, iopt_stc, iopt_gla,&
3464 iopt_sfc, iopt_trs)
3465
3466 implicit none
3467
3468 integer, intent(in) :: iopt_alb
3469 integer, intent(in) :: iopt_snf
3470 integer, intent(in) :: iopt_tbot
3471 integer, intent(in) :: iopt_stc
3473 integer, intent(in) :: iopt_gla
3474 integer, intent(in) :: iopt_sfc
3475 integer, intent(in) :: iopt_trs
3476
3477! -------------------------------------------------------------------------------------------------
3478
3479 opt_alb = iopt_alb
3480 opt_snf = iopt_snf
3481 opt_tbot = iopt_tbot
3482 opt_stc = iopt_stc
3483 opt_gla = iopt_gla
3484 opt_sfc = iopt_sfc
3485 opt_trs = iopt_trs
3486
3487 end subroutine noahmp_options_glacier
3488
3489end module noahmp_glacier_routines
3490! ==================================================================================================
3491
3493
3495 use noahmp_glacier_globals
3496
3497end module module_sf_noahmp_glacier
3498
subroutine, public stability(z1, zvfun, gdx, tv1, thv1, wind, z0max, ztmax, tvs, grav, thsfc_loc, rb, fm, fh, fm10, fh2, cm, ch, stress, ustar)
Definition sfc_diff.f:481
subroutine, private atm_glacier(ep_2, epsm1, sfcprs, sfctmp, q2, soldn, cosz, thair, qair, eair, rhoair, solad, solai, swdown)
re-process atmospheric forcing
subroutine esat(t, esw, esi, desw, desi)
subroutine, private thermoprop_glacier(nsoil, nsnow, isnow, dzsnso, dt, snowh, snice, snliq, df, hcpct, snicev, snliqv, epore, fact)
calculate thermal properties of soil, snow, lake, and frozen soil.
subroutine, private combo_glacier(dz, wliq, wice, t, dz2, wliq2, wice2, t2)
subroutine, private sfcdif1_glacier(iter,zlvl,zpd,z0h,z0m, qair,sfctmp,h,rhoair,mpe,ur, ifdef ccpp
compute surface drag coefficient cm for momentum and ch for heat
subroutine, private rosr12_glacier(p, a, b, c, d, delta, ntop, nsoil, nsnow)
subroutine, private snowh2o_glacier(nsnow, nsoil, dt, qsnfro, qsnsub, qrain, isnow, dzsnso, snowh, sneqv, snice, snliq, sh2o, sice, stc, ponding1, ponding2, fsh, qsnbot)
subroutine, private hrt_glacier(nsnow, nsoil, isnow, zsnso, stc, tbot, zbot, df, hcpct, ssoil, phi, ai, bi, ci, rhsts, botflx)
subroutine, private tsnosoi_glacier(nsoil, nsnow, isnow, dt, tbot, ssoil, snowh, zbot, zsnso, df, hcpct, stc)
subroutine sfcdif4(iloc, jloc, ux, vx, t1d, p1d, psfcpa, pblhx, dx, znt, ep_1, ep_2, cp, itime, snwh, isice, psi_opt, tsk, qx, zlvl, iz0tlnd, qsfc, hfx, qfx, cm, chs, chs2, cqs2, rmolx, ust, rbx, fmx, fhx, stressx, fm10x, fh2x, wspdx, flhcx, flqcx)
subroutine, private compact_glacier(nsnow, nsoil, dt, stc, snice, snliq, imelt, ficeold, isnow, dzsnso)
subroutine, private snowalb_bats_glacier(nband, cosz, fage, albsnd, albsni)
subroutine, private error_glacier(iloc,jloc,swdown,fsa,fsr,fira, fsh,fgev,ssoil,sag,prcp,edir, ifdef ccpp
subroutine, private combine_glacier(nsnow, nsoil, isnow, sh2o, stc, snice, snliq, dzsnso, sice, snowh, sneqv, ponding1, ponding2)
subroutine, private radiation_glacier(dt, tg, sneqvo, sneqv, cosz, qsnow, solad, solai, albold, tauss, sag, fsr, fsa, albsnd, albsni)
Compute solar radiation: absorbed & reflected by the ground.
subroutine, private divide_glacier(nsnow, nsoil, isnow, stc, snice, snliq, dzsnso)
subroutine, private snowwater_glacier(nsnow, nsoil, imelt, dt, sfctmp, snowhin, qsnow, qsnfro, qsnsub, qrain, ficeold, zsoil, isnow, snowh, sneqv, snice, snliq, sh2o, sice, stc, dzsnso, zsnso, fsh, qsnbot, snoflow, ponding1, ponding2)
subroutine, private water_glacier(nsnow, nsoil, imelt, dt, prcp, sfctmp, qvap, qdew, ficeold, zsoil, isnow, snowh, sneqv, snice, snliq, stc, dzsnso, sh2o, sice, ponding, zsnso, fsh, runsrf, runsub, qsnow, ponding1, ponding2, qsnbot, fpice, esnow)
subroutine, private snowfall_glacier(nsoil, nsnow, dt, qsnow, snowhin, sfctmp, isnow, snowh, dzsnso, stc, snice, snliq, sneqv)
subroutine, private phasechange_glacier(nsnow, nsoil, isnow, dt, fact, dzsnso, stc, snice, snliq, sneqv, snowh, smc, sh2o, qmelt, imelt, ponding)
subroutine, public noahmp_glacier( iloc,jloc,cosz,nsnow,nsoil,dt, sfctmp,sfcprs,uu,vv,q2,soldn, prcp,lwdn,tbot,zlvl,ficeold,zsoil, thsfc_loc,prslkix,prsik1x,prslk1x, psfc,pblhx,iz0tlnd,itime, sigmaf1,garea1,psi_opt, ep_1,ep_2,epsm1,cp, qsnow,sneqvo,albold,cm,ch,isnow, sneqv,smc,zsnso,snowh,snice,snliq, tg,stc,sh2o,tauss,qsfc, fsa,fsr,fira,fsh,fgev,ssoil, trad,edir,runsrf,runsub,sag,albedo, qsnbot,ponding,ponding1,ponding2,t2m, q2e,z0h_total, ifdef ccpp
subroutine, private csnow_glacier(isnow, nsnow, nsoil, snice, snliq, dzsnso, tksno, cvsno, snicev, snliqv, epore)
snow bulk density, volumetric capacity, and thermal conductivity
subroutine, private energy_glacier(nsnow,nsoil,isnow,dt,qsnow,rhoair, eair,sfcprs,qair,sfctmp,lwdn,uu, vv,solad,solai,cosz,zref, tbot,zbot,zsnso,dzsnso,sigmaf1,garea1, thsfc_loc,prslkix,prsik1x,prslk1x, psfc,pblhx,iz0tlnd,itime,psi_opt, ep_1, ep_2, epsm1, cp, tg,stc,snowh,sneqv,sneqvo,sh2o, smc,snice,snliq,albold,cm,ch, ifdef ccpp
Compute energy budget (momentum & energy fluxes and phase changes).
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,...
subroutine, private snowalb_class_glacier(nband, qsnow, dt, alb, albold, albsnd, albsni)
subroutine, private snow_age_glacier(dt, tg, sneqvo, sneqv, tauss, fage)
subroutine, private hstep_glacier(nsnow, nsoil, isnow, dt, ai, bi, ci, rhsts, stc)
subroutine, private glacier_flux(nsoil,nsnow,emg,isnow,df,dzsnso,z0m, zlvl,zpd,qair,sfctmp,rhoair,sfcprs, ur,gamma,rsurf,lwdn,rhsur,smc, eair,stc,sag,snowh,lathea,sh2o, thsfc_loc,prslkix,prsik1x,prslk1x, psfc,pblhx,iz0tlnd,itime,uu,vv, sigmaf1,garea1,psi_opt,ep_1, ep_2, epsm1, cp, ifdef ccpp
use newton-raphson iteration to solve ground (tg) temperature that balances the surface energy budget...
subroutine, public noahmp_options_glacier(iopt_alb, iopt_snf, iopt_tbot, iopt_stc, iopt_gla, iopt_sfc, iopt_trs)
This module contains the CCPP-compliant GFS surface layer scheme.
Definition sfc_diff.f:7