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