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