6module noahmp_glacier_globals
8 use machine ,
only : kind_phys
10 use module_sf_noahmplsm,
only :
sfcdif4
19 real (kind=kind_phys),
parameter :: grav = 9.80616
20 real (kind=kind_phys),
parameter :: sb = 5.67e-08
21 real (kind=kind_phys),
parameter :: vkc = 0.40
22 real (kind=kind_phys),
parameter :: tfrz = 273.16
23 real (kind=kind_phys),
parameter :: hsub = 2.8440e06
24 real (kind=kind_phys),
parameter :: hvap = 2.5104e06
25 real (kind=kind_phys),
parameter :: hfus = 0.3336e06
26 real (kind=kind_phys),
parameter :: cwat = 4.188e06
27 real (kind=kind_phys),
parameter :: cice = 2.094e06
28 real (kind=kind_phys),
parameter :: cpair = 1004.64
29 real (kind=kind_phys),
parameter :: tkwat = 0.6
30 real (kind=kind_phys),
parameter :: tkice = 2.2
31 real (kind=kind_phys),
parameter :: tkair = 0.023
32 real (kind=kind_phys),
parameter :: rair = 287.04
33 real (kind=kind_phys),
parameter :: rw = 461.269
34 real (kind=kind_phys),
parameter :: denh2o = 1000.
35 real (kind=kind_phys),
parameter :: denice = 917.
70 REAL,
PARAMETER :: Z0SNO = 0.002
71 REAL,
PARAMETER :: SSI = 0.03
72 REAL,
PARAMETER :: SWEMX = 1.00
76end module noahmp_glacier_globals
81 use noahmp_glacier_globals
123 iloc ,jloc ,cosz ,nsnow ,nsoil ,dt , & ! in : time/space/model-related
124 sfctmp ,sfcprs ,uu ,vv ,q2 ,soldn , & ! in : forcing
125 prcp ,lwdn ,tbot ,zlvl ,ficeold ,zsoil , & ! in : forcing
126 thsfc_loc ,prslkix ,prsik1x ,prslk1x , &
127 psfc ,pblhx ,iz0tlnd ,itime , &
128 sigmaf1 ,garea1 ,psi_opt , & ! in :
129 ep_1 ,ep_2 ,epsm1 ,cp , &
130 qsnow ,sneqvo ,albold ,cm ,ch ,isnow , & ! in/out :
131 sneqv ,smc ,zsnso ,snowh ,snice ,snliq , & ! in/out :
132 tg ,stc ,sh2o ,tauss ,qsfc , & ! in/out :
133 fsa ,fsr ,fira ,fsh ,fgev ,ssoil , & ! out :
134 trad ,edir ,runsrf ,runsub ,sag ,albedo , & ! out :
135 qsnbot ,ponding ,ponding1 ,ponding2 ,t2m,q2e ,z0h_total , & ! out :
137 emissi ,fpice ,ch2b , esnow , albsnd , albsni , &
140 emissi ,fpice ,ch2b , esnow. , albsnd , albsni)
151 integer ,
intent(in) :: iloc
152 integer ,
intent(in) :: jloc
153 real (kind=kind_phys) ,
intent(in) :: cosz
154 integer ,
intent(in) :: nsnow
155 integer ,
intent(in) :: nsoil
156 integer ,
intent(in) :: psi_opt
158 real (kind=kind_phys) ,
intent(in) :: dt
159 real (kind=kind_phys) ,
intent(in) :: sfctmp
160 real (kind=kind_phys) ,
intent(in) :: sfcprs
161 real (kind=kind_phys) ,
intent(in) :: uu
162 real (kind=kind_phys) ,
intent(in) :: vv
163 real (kind=kind_phys) ,
intent(in) :: q2
164 real (kind=kind_phys) ,
intent(in) :: soldn
165 real (kind=kind_phys) ,
intent(in) :: prcp
166 real (kind=kind_phys) ,
intent(in) :: lwdn
167 real (kind=kind_phys) ,
intent(in) :: tbot
168 real (kind=kind_phys) ,
intent(in) :: zlvl
169 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(in) :: ficeold
170 real (kind=kind_phys),
dimension( 1:nsoil),
intent(in) :: zsoil
171 logical ,
intent(in) :: thsfc_loc
172 real (kind=kind_phys) ,
intent(in) :: prslkix
173 real (kind=kind_phys) ,
intent(in) :: prsik1x
174 real (kind=kind_phys) ,
intent(in) :: prslk1x
176 real (kind=kind_phys) ,
intent(in) :: psfc
177 real (kind=kind_phys) ,
intent(in) :: pblhx
178 real (kind=kind_phys) ,
intent(in) :: ep_1
179 real (kind=kind_phys) ,
intent(in) :: ep_2
180 real (kind=kind_phys) ,
intent(in) :: epsm1
181 real (kind=kind_phys) ,
intent(in) :: cp
182 integer ,
intent(in) :: iz0tlnd
183 integer ,
intent(in) :: itime
185 real (kind=kind_phys) ,
intent(in) :: sigmaf1
186 real (kind=kind_phys) ,
intent(in) :: garea1
191 real (kind=kind_phys) ,
intent(inout) :: qsnow
192 real (kind=kind_phys) ,
intent(inout) :: sneqvo
193 real (kind=kind_phys) ,
intent(inout) :: albold
194 real (kind=kind_phys) ,
intent(inout) :: cm
195 real (kind=kind_phys) ,
intent(inout) :: ch
198 integer ,
intent(inout) :: isnow
199 real (kind=kind_phys) ,
intent(inout) :: sneqv
200 real (kind=kind_phys),
dimension( 1:nsoil),
intent(inout) :: smc
201 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: zsnso
202 real (kind=kind_phys) ,
intent(inout) :: snowh
203 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(inout) :: snice
204 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(inout) :: snliq
205 real (kind=kind_phys) ,
intent(inout) :: tg
206 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: stc
207 real (kind=kind_phys),
dimension( 1:nsoil),
intent(inout) :: sh2o
208 real (kind=kind_phys) ,
intent(inout) :: tauss
209 real (kind=kind_phys) ,
intent(inout) :: qsfc
212 real (kind=kind_phys) ,
intent(out) :: fsa
213 real (kind=kind_phys) ,
intent(out) :: fsr
214 real (kind=kind_phys) ,
intent(out) :: fira
215 real (kind=kind_phys) ,
intent(out) :: fsh
216 real (kind=kind_phys) ,
intent(out) :: fgev
217 real (kind=kind_phys) ,
intent(out) :: ssoil
218 real (kind=kind_phys) ,
intent(out) :: trad
219 real (kind=kind_phys) ,
intent(out) :: edir
220 real (kind=kind_phys) ,
intent(out) :: runsrf
221 real (kind=kind_phys) ,
intent(out) :: runsub
222 real (kind=kind_phys) ,
intent(out) :: sag
223 real (kind=kind_phys) ,
intent(out) ::
albedo
224 real (kind=kind_phys) ,
intent(out) :: qsnbot
225 real (kind=kind_phys) ,
intent(out) :: ponding
226 real (kind=kind_phys) ,
intent(out) :: ponding1
227 real (kind=kind_phys) ,
intent(out) :: ponding2
228 real (kind=kind_phys) ,
intent(out) :: t2m
229 real (kind=kind_phys) ,
intent(out) :: q2e
230 real (kind=kind_phys) ,
intent(out) :: z0h_total
231 real (kind=kind_phys) ,
intent(out) :: emissi
232 real (kind=kind_phys) ,
intent(out) :: fpice
233 real (kind=kind_phys) ,
intent(out) :: ch2b
234 real (kind=kind_phys) ,
intent(out) :: esnow
235 real (kind=kind_phys),
dimension(1:2) ,
intent(out) :: albsnd
236 real (kind=kind_phys),
dimension(1:2) ,
intent(out) :: albsni
240 character(len=*),
intent(inout) :: errmsg
241 integer,
intent(inout) :: errflg
246 integer,
dimension(-nsnow+1:nsoil) :: imelt
247 real (kind=kind_phys) :: rhoair
248 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: dzsnso
249 real (kind=kind_phys) :: thair
250 real (kind=kind_phys) :: qair
251 real (kind=kind_phys) :: eair
252 real (kind=kind_phys),
dimension( 1: 2) :: solad
253 real (kind=kind_phys),
dimension( 1: 2) :: solai
254 real (kind=kind_phys),
dimension( 1:nsoil) :: sice
255 real (kind=kind_phys),
dimension(-nsnow+1: 0) :: snicev
256 real (kind=kind_phys),
dimension(-nsnow+1: 0) :: snliqv
257 real (kind=kind_phys),
dimension(-nsnow+1: 0) :: epore
258 real (kind=kind_phys) :: qdew
259 real (kind=kind_phys) :: qvap
260 real (kind=kind_phys) :: lathea
261 real (kind=kind_phys) :: qmelt
262 real (kind=kind_phys) :: swdown
263 real (kind=kind_phys) :: beg_wb
264 real (kind=kind_phys) :: zbot = -8.0
266 character*256 message
271 call atm_glacier (ep_2, epsm1,sfcprs ,sfctmp ,q2 ,soldn ,cosz ,thair , &
272 qair ,eair ,rhoair ,solad ,solai ,swdown )
278 do iz = isnow+1, nsoil
279 if(iz == isnow+1)
then
280 dzsnso(iz) = - zsnso(iz)
282 dzsnso(iz) = zsnso(iz-1) - zsnso(iz)
289 eair ,sfcprs ,qair ,sfctmp ,lwdn ,uu , &
290 vv ,solad ,solai ,cosz ,zlvl , &
291 tbot ,zbot ,zsnso ,dzsnso ,sigmaf1 ,garea1 , &
292 thsfc_loc ,prslkix ,prsik1x ,prslk1x , &
293 psfc ,pblhx ,iz0tlnd ,itime ,psi_opt , &
294 ep_1, ep_2, epsm1, cp, &
295 tg ,stc ,snowh ,sneqv ,sneqvo ,sh2o , &
296 smc ,snice ,snliq ,albold ,cm ,ch , &
298 tauss ,qsfc ,errmsg ,errflg , &
302 imelt ,snicev ,snliqv ,epore ,qmelt ,ponding , &
303 sag ,fsa ,fsr ,fira ,fsh ,fgev , &
304 trad ,t2m ,ssoil ,lathea ,q2e ,emissi , &
305 ch2b ,albsnd ,albsni ,z0h_total)
308 if (errflg /= 0)
return
311 sice = max(0.0, smc - sh2o)
314 qvap = max( fgev/lathea, 0.)
315 qdew = abs( min(fgev/lathea, 0.))
320 call water_glacier (nsnow ,nsoil ,imelt ,dt ,prcp ,sfctmp , &
321 qvap ,qdew ,ficeold ,zsoil , &
322 isnow ,snowh ,sneqv ,snice ,snliq ,stc , &
323 dzsnso ,sh2o ,sice ,ponding ,zsnso ,fsh , &
324 runsrf ,runsub ,qsnow ,ponding1 ,ponding2 ,qsnbot , &
327 if(opt_gla == 2)
then
340 fsh ,fgev ,ssoil ,sag ,prcp ,edir , &
342 runsrf ,runsub ,sneqv ,dt ,beg_wb ,errmsg , errflg )
344 runsrf ,runsub ,sneqv ,dt ,beg_wb )
348 if (errflg /= 0)
return
351 if(snowh <= 1.e-6 .or. sneqv <= 1.e-3)
then
356 if(swdown.ne.0.)
then
367 subroutine atm_glacier (ep_2, epsm1, sfcprs ,sfctmp ,q2 ,soldn ,cosz ,thair , &
368 qair ,eair ,rhoair ,solad ,solai , &
377 real (kind=kind_phys) ,
intent(in) :: ep_2
378 real (kind=kind_phys) ,
intent(in) :: epsm1
379 real (kind=kind_phys) ,
intent(in) :: sfcprs
380 real (kind=kind_phys) ,
intent(in) :: sfctmp
381 real (kind=kind_phys) ,
intent(in) :: q2
382 real (kind=kind_phys) ,
intent(in) :: soldn
383 real (kind=kind_phys) ,
intent(in) :: cosz
387 real (kind=kind_phys) ,
intent(out) :: thair
388 real (kind=kind_phys) ,
intent(out) :: qair
389 real (kind=kind_phys) ,
intent(out) :: eair
390 real (kind=kind_phys),
dimension( 1: 2),
intent(out) :: solad
391 real (kind=kind_phys),
dimension( 1: 2),
intent(out) :: solai
392 real (kind=kind_phys) ,
intent(out) :: rhoair
393 real (kind=kind_phys) ,
intent(out) :: swdown
397 real (kind=kind_phys) :: pair
401 thair = sfctmp * (sfcprs/pair)**(rair/cpair)
405 eair = qair*sfcprs / (ep_2-epsm1*qair)
406 rhoair = (sfcprs+epsm1*eair) / (rair*sfctmp)
414 solad(1) = swdown*0.7*0.5
415 solad(2) = swdown*0.7*0.5
416 solai(1) = swdown*0.3*0.5
417 solai(2) = swdown*0.3*0.5
425 eair ,sfcprs ,qair ,sfctmp ,lwdn ,uu , & !in
426 vv ,solad ,solai ,cosz ,zref , & !in
427 tbot ,zbot ,zsnso ,dzsnso ,sigmaf1 ,garea1 , & !in
428 thsfc_loc ,prslkix ,prsik1x ,prslk1x , & !in
429 psfc ,pblhx ,iz0tlnd ,itime ,psi_opt , &
430 ep_1, ep_2, epsm1, cp, &
431 tg ,stc ,snowh ,sneqv ,sneqvo ,sh2o , & !inout
432 smc ,snice ,snliq ,albold ,cm ,ch , & !inout
434 tauss ,qsfc ,errmsg ,errflg , & !inout
436 tauss ,qsfc , & !inout
438 imelt ,snicev ,snliqv ,epore ,qmelt ,ponding , & !out
439 sag ,fsa ,fsr ,fira ,fsh ,fgev , & !out
440 trad ,t2m ,ssoil ,lathea ,q2e ,emissi , & !out
441 ch2b ,albsnd ,albsni ,z0h_total)
451 integer ,
intent(in) :: nsnow
452 integer ,
intent(in) :: nsoil
453 integer ,
intent(in) :: psi_opt
455 integer ,
intent(in) :: isnow
456 real (kind=kind_phys) ,
intent(in) :: dt
457 real (kind=kind_phys) ,
intent(in) :: qsnow
458 real (kind=kind_phys) ,
intent(in) :: rhoair
459 real (kind=kind_phys) ,
intent(in) :: eair
460 real (kind=kind_phys) ,
intent(in) :: sfcprs
461 real (kind=kind_phys) ,
intent(in) :: qair
462 real (kind=kind_phys) ,
intent(in) :: sfctmp
463 real (kind=kind_phys) ,
intent(in) :: lwdn
464 real (kind=kind_phys) ,
intent(in) :: uu
465 real (kind=kind_phys) ,
intent(in) :: vv
466 real (kind=kind_phys) ,
dimension( 1: 2),
intent(in) :: solad
467 real (kind=kind_phys) ,
dimension( 1: 2),
intent(in) :: solai
468 real (kind=kind_phys) ,
intent(in) :: cosz
469 real (kind=kind_phys) ,
intent(in) :: zref
470 real (kind=kind_phys) ,
intent(in) :: tbot
471 real (kind=kind_phys) ,
intent(in) :: zbot
472 real (kind=kind_phys) ,
dimension(-nsnow+1:nsoil),
intent(in) :: zsnso
473 real (kind=kind_phys) ,
dimension(-nsnow+1:nsoil),
intent(in) :: dzsnso
475 logical ,
intent(in) :: thsfc_loc
476 real (kind=kind_phys) ,
intent(in) :: prslkix
477 real (kind=kind_phys) ,
intent(in) :: prsik1x
478 real (kind=kind_phys) ,
intent(in) :: prslk1x
480 real (kind=kind_phys) ,
intent(in) :: pblhx
481 real (kind=kind_phys) ,
intent(in) :: psfc
482 real (kind=kind_phys) ,
intent(in) :: ep_1
483 real (kind=kind_phys) ,
intent(in) :: ep_2
484 real (kind=kind_phys) ,
intent(in) :: epsm1
485 real (kind=kind_phys) ,
intent(in) :: cp
486 integer ,
intent(in) :: iz0tlnd
487 integer ,
intent(in) :: itime
489 real (kind=kind_phys) ,
intent(in) :: sigmaf1
490 real (kind=kind_phys) ,
intent(in) :: garea1
493 real (kind=kind_phys) ,
intent(inout) :: tg
494 real (kind=kind_phys) ,
dimension(-nsnow+1:nsoil),
intent(inout) :: stc
495 real (kind=kind_phys) ,
intent(inout) :: snowh
496 real (kind=kind_phys) ,
intent(inout) :: sneqv
497 real (kind=kind_phys) ,
intent(inout) :: sneqvo
498 real (kind=kind_phys) ,
dimension( 1:nsoil),
intent(inout) :: sh2o
499 real (kind=kind_phys) ,
dimension( 1:nsoil),
intent(inout) :: smc
500 real (kind=kind_phys) ,
dimension(-nsnow+1: 0),
intent(inout) :: snice
501 real (kind=kind_phys) ,
dimension(-nsnow+1: 0),
intent(inout) :: snliq
502 real (kind=kind_phys) ,
intent(inout) :: albold
503 real (kind=kind_phys) ,
intent(inout) :: cm
504 real (kind=kind_phys) ,
intent(inout) :: ch
505 real (kind=kind_phys) ,
intent(inout) :: tauss
506 real (kind=kind_phys) ,
intent(inout) :: qsfc
509 character(len=*) ,
intent(inout) :: errmsg
510 integer ,
intent(inout) :: errflg
514 integer,
dimension(-nsnow+1:nsoil) ,
intent(out) :: imelt
515 real (kind=kind_phys) ,
dimension(-nsnow+1: 0),
intent(out) :: snicev
516 real (kind=kind_phys) ,
dimension(-nsnow+1: 0),
intent(out) :: snliqv
517 real (kind=kind_phys) ,
dimension(-nsnow+1: 0),
intent(out) :: epore
518 real (kind=kind_phys) ,
intent(out) :: qmelt
519 real (kind=kind_phys) ,
intent(out) :: ponding
520 real (kind=kind_phys) ,
intent(out) :: sag
521 real (kind=kind_phys) ,
intent(out) :: fsa
522 real (kind=kind_phys) ,
intent(out) :: fsr
523 real (kind=kind_phys) ,
intent(out) :: fira
524 real (kind=kind_phys) ,
intent(out) :: fsh
525 real (kind=kind_phys) ,
intent(out) :: fgev
526 real (kind=kind_phys) ,
intent(out) :: trad
527 real (kind=kind_phys) ,
intent(out) :: t2m
528 real (kind=kind_phys) ,
intent(out) :: ssoil
529 real (kind=kind_phys) ,
intent(out) :: lathea
530 real (kind=kind_phys) ,
intent(out) :: q2e
531 real (kind=kind_phys) ,
intent(out) :: emissi
532 real (kind=kind_phys) ,
intent(out) :: ch2b
533 real (kind=kind_phys),
dimension(1:2) ,
intent(out) :: albsnd
534 real (kind=kind_phys),
dimension(1:2) ,
intent(out) :: albsni
535 real (kind=kind_phys) ,
intent(out) :: z0h_total
539 real (kind=kind_phys) :: ur
540 real (kind=kind_phys) :: zlvl
541 real (kind=kind_phys) :: rsurf
542 real (kind=kind_phys) :: zpd
543 real (kind=kind_phys) :: z0mg
544 real (kind=kind_phys) :: emg
545 real (kind=kind_phys) :: fire
546 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: fact
547 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: df
548 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: hcpct
549 real (kind=kind_phys) :: gamma
550 real (kind=kind_phys) :: rhsur
556 ur = max( sqrt(uu**2.+vv**2.), 1. )
568 dt ,snowh ,snice ,snliq , &
569 df ,hcpct ,snicev ,snliqv ,epore , &
575 qsnow ,solad ,solai , &
577 sag ,fsr ,fsa ,albsnd ,albsni)
591 gamma = cpair*sfcprs/(ep_2*lathea)
595 call glacier_flux (nsoil ,nsnow ,emg ,isnow ,df ,dzsnso ,z0mg , &
596 zlvl ,zpd ,qair ,sfctmp ,rhoair ,sfcprs , &
597 ur ,gamma ,rsurf ,lwdn ,rhsur ,smc , &
598 eair ,stc ,sag ,snowh ,lathea ,sh2o , &
599 thsfc_loc ,prslkix ,prsik1x ,prslk1x , &
600 psfc ,pblhx ,iz0tlnd ,itime ,uu ,vv , &
601 sigmaf1 ,garea1 ,psi_opt ,ep_1, ep_2, epsm1, cp, &
603 cm ,ch ,tg ,qsfc ,errmsg ,errflg , &
607 fira ,fsh ,fgev ,ssoil , &
608 t2m ,q2e ,ch2b ,z0h_total)
617 errmsg =
"stop in noah-mp: emitted longwave <0"
620 call wrf_error_fatal(
"stop in noah-mp: emitted longwave <0")
631 trad = ( ( fire - (1-emissi)*lwdn ) / (emissi*sb) ) ** 0.25
636 ssoil ,snowh ,zbot ,zsnso ,df , &
641 if(opt_stc == 2)
then
642 if (snowh > 0.05 .and. tg > tfrz) tg = tfrz
649 stc ,snice ,snliq ,sneqv ,snowh , &
651 qmelt ,imelt ,ponding )
659 dt ,snowh ,snice ,snliq , & !in
660 df ,hcpct ,snicev ,snliqv ,epore , & !out
667 integer ,
intent(in) :: nsoil
668 integer ,
intent(in) :: nsnow
669 integer ,
intent(in) :: isnow
670 real (kind=kind_phys) ,
intent(in) :: dt
671 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(in) :: snice
672 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(in) :: snliq
673 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: dzsnso
674 real (kind=kind_phys) ,
intent(in) :: snowh
677 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(out) :: df
678 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(out) :: hcpct
679 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(out) :: snicev
680 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(out) :: snliqv
681 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(out) :: epore
682 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(out) :: fact
687 real (kind=kind_phys),
dimension(-nsnow+1: 0) :: cvsno
688 real (kind=kind_phys),
dimension(-nsnow+1: 0) :: tksno
689 real (kind=kind_phys) :: zmid
694 call csnow_glacier (isnow ,nsnow ,nsoil ,snice ,snliq ,dzsnso , &
695 tksno ,cvsno ,snicev ,snliqv ,epore )
699 hcpct(iz) = cvsno(iz)
705 zmid = 0.5 * (dzsnso(iz))
707 zmid = zmid + dzsnso(iz2)
709 hcpct(iz) = 1.e6 * ( 0.8194 + 0.1309*zmid )
710 df(iz) = 0.32333 + ( 0.10073 * zmid )
715 do iz = isnow+1,nsoil
716 fact(iz) = dt/(hcpct(iz)*dzsnso(iz))
722 df(1) = (df(1)*dzsnso(1)+0.35*snowh) / (snowh +dzsnso(1))
724 df(1) = (df(1)*dzsnso(1)+df(0)*dzsnso(0)) / (dzsnso(0)+dzsnso(1))
734 tksno ,cvsno ,snicev ,snliqv ,epore )
742 integer,
intent(in) :: isnow
743 integer ,
intent(in) :: nsnow
744 integer ,
intent(in) :: nsoil
745 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(in) :: snice
746 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(in) :: snliq
747 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: dzsnso
751 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(out) :: cvsno
752 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(out) :: tksno
753 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(out) :: snicev
754 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(out) :: snliqv
755 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(out) :: epore
760 real (kind=kind_phys),
dimension(-nsnow+1: 0) :: bdsnoi
766 snicev(iz) = min(1., snice(iz)/(dzsnso(iz)*denice) )
767 epore(iz) = 1. - snicev(iz)
768 snliqv(iz) = min(epore(iz),snliq(iz)/(dzsnso(iz)*denh2o))
772 bdsnoi(iz) = (snice(iz)+snliq(iz))/dzsnso(iz)
773 cvsno(iz) = cice*snicev(iz)+cwat*snliqv(iz)
783 tksno(iz) = 2.576e-6*bdsnoi(iz)**2. + 0.074
792 qsnow ,solad ,solai , & !in
793 albold ,tauss , & !inout
794 sag ,fsr ,fsa ,albsnd ,albsni)
799 real (kind=kind_phys),
intent(in) :: dt
800 real (kind=kind_phys),
intent(in) :: tg
801 real (kind=kind_phys),
intent(in) :: sneqvo
802 real (kind=kind_phys),
intent(in) :: sneqv
803 real (kind=kind_phys),
intent(in) :: cosz
804 real (kind=kind_phys),
intent(in) :: qsnow
805 real (kind=kind_phys),
dimension(1:2) ,
intent(in) :: solad
806 real (kind=kind_phys),
dimension(1:2) ,
intent(in) :: solai
809 real (kind=kind_phys),
intent(inout) :: albold
810 real (kind=kind_phys),
intent(inout) :: tauss
811 real (kind=kind_phys),
dimension(1:2) :: albsnd
812 real (kind=kind_phys),
dimension(1:2) :: albsni
815 real (kind=kind_phys),
intent(out) :: sag
816 real (kind=kind_phys),
intent(out) :: fsr
817 real (kind=kind_phys),
intent(out) :: fsa
822 real (kind=kind_phys) :: fage
823 real (kind=kind_phys) :: alb
824 real (kind=kind_phys) :: abs
825 real (kind=kind_phys) :: ref
826 real (kind=kind_phys) :: fsno
827 real (kind=kind_phys),
dimension(1:2) :: albice
829 real (kind=kind_phys),
parameter :: mpe = 1.e-6
847 if(opt_alb == 2)
then
859 if(sneqv > 0.0) fsno = 1.0
865 albsnd(ib) = albice(ib)*(1.-fsno) + albsnd(ib)*fsno
866 albsni(ib) = albice(ib)*(1.-fsno) + albsni(ib)*fsno
870 abs = solad(ib)*(1.-albsnd(ib)) + solai(ib)*(1.-albsni(ib))
874 ref = solad(ib)*albsnd(ib) + solai(ib)*albsni(ib)
889 real (kind=kind_phys),
intent(in) :: dt
890 real (kind=kind_phys),
intent(in) :: tg
891 real (kind=kind_phys),
intent(in) :: sneqvo
892 real (kind=kind_phys),
intent(in) :: sneqv
895 real (kind=kind_phys),
intent(inout) :: tauss
898 real (kind=kind_phys),
intent(out) :: fage
901 real (kind=kind_phys) :: tage
902 real (kind=kind_phys) :: age1
903 real (kind=kind_phys) :: age2
904 real (kind=kind_phys) :: age3
905 real (kind=kind_phys) :: dela
906 real (kind=kind_phys) :: sge
907 real (kind=kind_phys) :: dels
908 real (kind=kind_phys) :: dela0
909 real (kind=kind_phys) :: arg
913 if(sneqv.le.0.0)
then
915 else if (sneqv.gt.800.)
then
920 arg = 5.e3*(1./tfrz-1./tg)
922 age2 = exp(amin1(0.,10.*arg))
924 tage = age1+age2+age3
926 dels = amax1(0.0,sneqv-sneqvo) / swemx
927 sge = (tauss+dela)*(1.0-dels)
928 tauss = amax1(0.,sge)
931 fage= tauss/(tauss+1.)
943 integer,
intent(in) :: nband
945 real (kind=kind_phys),
intent(in) :: cosz
946 real (kind=kind_phys),
intent(in) :: fage
950 real (kind=kind_phys),
dimension(1:2),
intent(out) :: albsnd
951 real (kind=kind_phys),
dimension(1:2),
intent(out) :: albsni
954 real (kind=kind_phys) :: fzen
955 real (kind=kind_phys) :: cf1
956 real (kind=kind_phys) :: sl2
957 real (kind=kind_phys) :: sl1
958 real (kind=kind_phys) :: sl
959 real (kind=kind_phys),
parameter :: c1 = 0.2
960 real (kind=kind_phys),
parameter :: c2 = 0.5
966 albsnd(1: nband) = 0.
967 albsni(1: nband) = 0.
974 cf1=((1.+sl1)/(1.+sl2*cosz)-sl1)
977 albsni(1)=0.95*(1.-c1*fage)
978 albsni(2)=0.65*(1.-c2*fage)
980 albsnd(1)=albsni(1)+0.4*fzen*(1.-albsni(1))
981 albsnd(2)=albsni(2)+0.4*fzen*(1.-albsni(2))
993 integer,
intent(in) :: nband
995 real (kind=kind_phys),
intent(in) :: qsnow
996 real (kind=kind_phys),
intent(in) :: dt
997 real (kind=kind_phys),
intent(in) :: albold
1001 real (kind=kind_phys),
intent(inout) :: alb
1004 real (kind=kind_phys),
dimension(1:2),
intent(out) :: albsnd
1005 real (kind=kind_phys),
dimension(1:2),
intent(out) :: albsni
1011 albsnd(1: nband) = 0.
1012 albsni(1: nband) = 0.
1016 alb = 0.55 + (albold-0.55) * exp(-0.01*dt/3600.)
1021 if (qsnow > 0.)
then
1022 alb = alb + min(qsnow*dt,swemx) * (0.84-alb)/(swemx)
1036 zlvl ,zpd ,qair ,sfctmp ,rhoair ,sfcprs , & !in
1037 ur ,gamma ,rsurf ,lwdn ,rhsur ,smc , & !in
1038 eair ,stc ,sag ,snowh ,lathea ,sh2o , & !in
1039 thsfc_loc ,prslkix ,prsik1x ,prslk1x , &
1040 psfc ,pblhx ,iz0tlnd ,itime ,uu ,vv , &
1041 sigmaf1 ,garea1 ,psi_opt ,ep_1, ep_2, epsm1, cp, & !in
1043 cm ,ch ,tgb ,qsfc ,errmsg ,errflg , & !inout
1045 cm ,ch ,tgb ,qsfc , & !inout
1047 irb ,shb ,evb ,ghb , & !out
1048 t2mb ,q2b ,ehb2 ,z0h_total)
1062 integer,
intent(in) :: nsnow
1063 integer,
intent(in) :: nsoil
1064 integer,
intent(in) :: psi_opt
1066 real (kind=kind_phys),
intent(in) :: emg
1067 integer,
intent(in) :: isnow
1068 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: df
1069 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: dzsnso
1070 real (kind=kind_phys),
intent(in) :: z0m
1071 real (kind=kind_phys),
intent(in) :: zlvl
1072 real (kind=kind_phys),
intent(in) :: zpd
1073 real (kind=kind_phys),
intent(in) :: qair
1074 real (kind=kind_phys),
intent(in) :: sfctmp
1075 real (kind=kind_phys),
intent(in) :: rhoair
1076 real (kind=kind_phys),
intent(in) :: sfcprs
1077 real (kind=kind_phys),
intent(in) :: ur
1078 real (kind=kind_phys),
intent(in) :: gamma
1079 real (kind=kind_phys),
intent(in) :: rsurf
1080 real (kind=kind_phys),
intent(in) :: lwdn
1081 real (kind=kind_phys),
intent(in) :: rhsur
1082 real (kind=kind_phys),
intent(in) :: eair
1083 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: stc
1084 real (kind=kind_phys),
dimension( 1:nsoil),
intent(in) :: smc
1085 real (kind=kind_phys),
dimension( 1:nsoil),
intent(in) :: sh2o
1086 real (kind=kind_phys),
intent(in) :: sag
1087 real (kind=kind_phys),
intent(in) :: snowh
1088 real (kind=kind_phys),
intent(in) :: lathea
1090 logical ,
intent(in) :: thsfc_loc
1091 real (kind=kind_phys),
intent(in) :: prslkix
1092 real (kind=kind_phys),
intent(in) :: prsik1x
1093 real (kind=kind_phys),
intent(in) :: prslk1x
1095 real (kind=kind_phys) ,
intent(in) :: pblhx
1096 real (kind=kind_phys) ,
intent(in) :: psfc
1097 real (kind=kind_phys) ,
intent(in) :: ep_1
1098 real (kind=kind_phys) ,
intent(in) :: ep_2
1099 real (kind=kind_phys) ,
intent(in) :: epsm1
1100 real (kind=kind_phys) ,
intent(in) :: cp
1101 integer ,
intent(in) :: iz0tlnd
1102 integer ,
intent(in) :: itime
1103 real (kind=kind_phys) ,
intent(in) :: uu
1104 real (kind=kind_phys) ,
intent(in) :: vv
1106 real (kind=kind_phys),
intent(in) :: sigmaf1
1107 real (kind=kind_phys),
intent(in) :: garea1
1110 real (kind=kind_phys),
intent(inout) :: cm
1111 real (kind=kind_phys),
intent(inout) :: ch
1112 real (kind=kind_phys),
intent(inout) :: tgb
1113 real (kind=kind_phys),
intent(inout) :: qsfc
1116 character(len=*),
intent(inout) :: errmsg
1117 integer,
intent(inout) :: errflg
1122 real (kind=kind_phys),
intent(out) :: irb
1123 real (kind=kind_phys),
intent(out) :: shb
1124 real (kind=kind_phys),
intent(out) :: evb
1125 real (kind=kind_phys),
intent(out) :: ghb
1126 real (kind=kind_phys),
intent(out) :: t2mb
1127 real (kind=kind_phys),
intent(out) :: q2b
1128 real (kind=kind_phys),
intent(out) :: ehb2
1129 real (kind=kind_phys),
intent(out) :: z0h_total
1136 real (kind=kind_phys) :: mpe
1137 real (kind=kind_phys) :: dtg
1139 real (kind=kind_phys) :: mozold
1140 real (kind=kind_phys) :: fm2
1141 real (kind=kind_phys) :: fh2
1142 real (kind=kind_phys) :: ch2
1143 real (kind=kind_phys) :: h
1144 real (kind=kind_phys) :: fv
1145 real (kind=kind_phys) :: cir
1146 real (kind=kind_phys) :: cgh
1147 real (kind=kind_phys) :: csh
1148 real (kind=kind_phys) :: cev
1149 real (kind=kind_phys) :: cq2b
1151 real (kind=kind_phys) :: z0h
1153 real (kind=kind_phys) :: qfx
1154 real (kind=kind_phys) :: cq2
1157 real(kind=kind_phys) :: rb1i
1158 real(kind=kind_phys) :: fm10i
1160 real(kind=kind_phys) :: stress1i
1162 real(kind=kind_phys) :: wspd1i
1163 real(kind=kind_phys) :: flhc1i
1164 real(kind=kind_phys) :: flqc1i
1166 real(kind=kind_phys) :: tv1i
1168 real(kind=kind_phys) :: thv1i
1169 real(kind=kind_phys) :: tvsi
1170 real(kind=kind_phys) :: zlvli
1172 real(kind=kind_phys) :: snwd
1174 real(kind=kind_phys) :: reyni
1175 real(kind=kind_phys) :: virtfaci
1177 real(kind=kind_phys) :: tem1,tem2,zvfun1,gdx
1178 real(kind=kind_phys),
parameter :: z0lo=0.1, z0up=1.0
1180 real (kind=kind_phys) :: moz
1181 real (kind=kind_phys) :: fm
1182 real (kind=kind_phys) :: fh
1183 real (kind=kind_phys) :: ramb
1184 real (kind=kind_phys) :: rahb
1185 real (kind=kind_phys) :: rawb
1186 real (kind=kind_phys) :: estg
1187 real (kind=kind_phys) :: destg
1188 real (kind=kind_phys) :: esatw
1189 real (kind=kind_phys) :: esati
1190 real (kind=kind_phys) :: dsatw
1191 real (kind=kind_phys) :: dsati
1192 real (kind=kind_phys) :: a
1193 real (kind=kind_phys) :: b
1194 real (kind=kind_phys) :: t, tdc
1195 real (kind=kind_phys),
dimension( 1:nsoil) :: sice
1196 real (kind=kind_phys) :: czil
1198 tdc(t) = min( 50., max(-50.,(t-tfrz)) )
1225 fv = ur*vkc/log(zlvli/z0m)
1226 reyni = fv*z0m/(1.5e-05)
1228 if (opt_trs == 1)
then
1230 elseif (opt_trs == 2)
then
1231 z0h = z0m*exp(-czil*0.4*258.2*sqrt(fv*z0m))
1232 elseif (opt_trs == 3)
then
1234 elseif (opt_trs == 4)
then
1235 if (reyni .gt. 2.0)
then
1236 z0h = z0m/exp(2.46*(reyni)**0.25 - log(7.4))
1238 z0h = z0m/exp(-log(0.397))
1244 virtfaci = 1.0 + 0.61 * max(qair, 1.e-8)
1245 tv1i = sfctmp * virtfaci
1248 thv1i = sfctmp * prslkix * virtfaci
1250 thv1i = sfctmp / prslk1x * virtfaci
1253 if ( ur < 2.0) niter = 2
1256 cgh = 2.*df(isnow+1)/dzsnso(isnow+1)
1259 tem1 = (z0m - z0lo) / (z0up - z0lo)
1260 tem1 = min(max(tem1, 0.0_kind_phys), 1.0_kind_phys)
1261 tem2 = max(sigmaf1, 0.1_kind_phys)
1262 zvfun1= sqrt(tem1 * tem2)
1265 if(opt_sfc == 1 .or. opt_sfc == 2 .or. opt_sfc == 4)
then
1266 loop3:
do iter = 1, niterb
1267 if(opt_sfc == 1 .or. opt_sfc == 2)
then
1272 qair ,sfctmp ,h ,rhoair ,mpe ,ur , &
1274 & moz ,mozsgn ,fm ,fh ,fm2 ,fh2,fv, errmsg, errflg, &
1276 & moz ,mozsgn ,fm ,fh ,fm2 ,fh2,fv , &
1281 if (errflg /= 0)
return
1285 if(opt_sfc == 4)
then
1287 call sfcdif4(1 ,1 ,uu ,vv ,sfctmp , &
1288 sfcprs ,psfc ,pblhx ,gdx ,z0m , &
1290 itime ,snwd ,1 ,psi_opt, &
1291 tgb ,qair ,zlvl ,iz0tlnd,qsfc , &
1292 h ,qfx ,cm ,ch ,ch2 , &
1293 cq2 ,moz ,fv ,rb1i, fm, fh, &
1294 stress1i,fm10i ,fh2 , wspd1i ,flhc1i ,flqc1i)
1314 if(opt_sfc == 1 .or. opt_sfc == 2 .or. opt_sfc == 3)
then
1315 ramb = max(1.,1./(cm*ur))
1316 rahb = max(1.,1./(ch*ur))
1317 elseif(opt_sfc == 4)
then
1318 ramb = max(1.,1./(cm*wspd1i) )
1319 rahb = max(1.,1./(ch*wspd1i) )
1327 call esat(t, esatw, esati, dsatw, dsati)
1336 csh = rhoair*cpair/rahb
1337 if(snowh > 0.0 .or. opt_gla == 1)
then
1338 cev = rhoair*cpair/gamma/(rsurf+rawb)
1345 irb = cir * tgb**4 - emg*lwdn
1346 shb = csh * (tgb - sfctmp )
1347 evb = cev * (estg*rhsur - eair )
1348 ghb = cgh * (tgb - stc(isnow+1))
1350 b = sag-irb-shb-evb-ghb
1351 a = 4.*cir*tgb**3 + csh + cev*destg + cgh
1354 irb = irb + 4.*cir*tgb**3*dtg
1356 evb = evb + cev*destg*dtg
1363 h = csh * (tgb - sfctmp)
1366 call esat(t, esatw, esati, dsatw, dsati)
1372 qsfc = ep_2*(estg*rhsur)/(sfcprs+epsm1*(estg*rhsur))
1373 qfx = (qsfc-qair)*cev*gamma/cpair
1378 if (opt_sfc == 3)
then
1383 tvsi = tgb * virtfaci
1385 tvsi = tgb/prsik1x * virtfaci
1389 (zlvli, zvfun1, gdx,tv1i,thv1i, ur, z0m, z0h, tvsi, grav,thsfc_loc, &
1390 rb1i, fm,fh,fm10i,fh2,cm,ch,stress1i,fv)
1394 ramb = max(1.,1./(cm*ur))
1395 rahb = max(1.,1./(ch*ur))
1401 call esat(t, esatw, esati, dsatw, dsati)
1410 csh = rhoair*cpair/rahb
1412 if(snowh > 0.0 .or. opt_gla == 1)
then
1413 cev = rhoair*cpair/gamma/(rsurf+rawb)
1420 irb = cir * tgb**4 - emg*lwdn
1421 shb = csh * (tgb - sfctmp )
1422 evb = cev * (estg*rhsur - eair )
1423 ghb = cgh * (tgb - stc(isnow+1))
1425 b = sag-irb-shb-evb-ghb
1426 a = 4.*cir*tgb**3 + csh + cev*destg + cgh
1429 irb = irb + 4.*cir*tgb**3*dtg
1431 evb = evb + cev*destg*dtg
1439 call esat(t, esatw, esati, dsatw, dsati)
1445 qsfc = ep_2*(estg*rhsur)/(sfcprs+epsm1*(estg*rhsur))
1455 if(opt_stc == 1 .or. opt_stc ==3)
then
1456 if ((maxval(sice) > 0.0 .or. snowh > 0.0) .and. tgb > tfrz .and. opt_gla == 1)
then
1459 call esat(t, esatw, esati, dsatw, dsati)
1461 qsfc = ep_2*(estg*rhsur)/(sfcprs+epsm1*(estg*rhsur))
1462 irb = cir * tgb**4 - emg*lwdn
1463 shb = csh * (tgb - sfctmp)
1464 evb = cev * (estg*rhsur - eair )
1465 ghb = sag - (irb+shb+evb)
1470 ehb2 = fv*vkc/(log((2.+z0h)/z0h)-fh2)
1473 if (opt_sfc ==3)
then
1478 if (opt_sfc == 4)
then
1483 if (ehb2.lt.1.e-5 )
then
1487 t2mb = tgb - shb/(rhoair*cpair) * 1./ehb2
1488 q2b = qsfc - evb/(lathea*rhoair)*(1./cq2b + rsurf)
1497 subroutine esat(t, esw, esi, desw, desi)
1505 real (kind=kind_phys),
intent(in) :: t
1509 real (kind=kind_phys),
intent(out) :: esw
1510 real (kind=kind_phys),
intent(out) :: esi
1511 real (kind=kind_phys),
intent(out) :: desw
1512 real (kind=kind_phys),
intent(out) :: desi
1516 real (kind=kind_phys) :: a0,a1,a2,a3,a4,a5,a6
1517 real (kind=kind_phys) :: b0,b1,b2,b3,b4,b5,b6
1518 real (kind=kind_phys) :: c0,c1,c2,c3,c4,c5,c6
1519 real (kind=kind_phys) :: d0,d1,d2,d3,d4,d5,d6
1521 parameter(a0=6.107799961 , a1=4.436518521e-01, &
1522 a2=1.428945805e-02, a3=2.650648471e-04, &
1523 a4=3.031240396e-06, a5=2.034080948e-08, &
1526 parameter(b0=6.109177956 , b1=5.034698970e-01, &
1527 b2=1.886013408e-02, b3=4.176223716e-04, &
1528 b4=5.824720280e-06, b5=4.838803174e-08, &
1531 parameter(c0= 4.438099984e-01, c1=2.857002636e-02, &
1532 c2= 7.938054040e-04, c3=1.215215065e-05, &
1533 c4= 1.036561403e-07, c5=3.532421810e-10, &
1534 c6=-7.090244804e-13)
1536 parameter(d0=5.030305237e-01, d1=3.773255020e-02, &
1537 d2=1.267995369e-03, d3=2.477563108e-05, &
1538 d4=3.005693132e-07, d5=2.158542548e-09, &
1541 esw = 100.*(a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*a6))))))
1542 esi = 100.*(b0+t*(b1+t*(b2+t*(b3+t*(b4+t*(b5+t*b6))))))
1543 desw = 100.*(c0+t*(c1+t*(c2+t*(c3+t*(c4+t*(c5+t*c6))))))
1544 desi = 100.*(d0+t*(d1+t*(d2+t*(d3+t*(d4+t*(d5+t*d6))))))
1551 qair ,sfctmp ,h ,rhoair ,mpe ,ur , & !in
1553 & moz ,mozsgn ,fm ,fh ,fm2 ,fh2 , fv, & !inout
1554 & errmsg ,errflg , & !inout
1556 & moz ,mozsgn ,fm ,fh ,fm2 ,fh2 , fv, & !inout
1565 integer,
intent(in) :: iter
1566 real (kind=kind_phys),
intent(in) :: zlvl
1567 real (kind=kind_phys),
intent(in) :: zpd
1568 real (kind=kind_phys),
intent(in) :: z0h
1569 real (kind=kind_phys),
intent(in) :: z0m
1570 real (kind=kind_phys),
intent(in) :: qair
1571 real (kind=kind_phys),
intent(in) :: sfctmp
1572 real (kind=kind_phys),
intent(in) :: h
1573 real (kind=kind_phys),
intent(in) :: rhoair
1574 real (kind=kind_phys),
intent(in) :: mpe
1575 real (kind=kind_phys),
intent(in) :: ur
1578 real (kind=kind_phys),
intent(inout) :: moz
1579 integer,
intent(inout) :: mozsgn
1580 real (kind=kind_phys),
intent(inout) :: fm
1581 real (kind=kind_phys),
intent(inout) :: fh
1582 real (kind=kind_phys),
intent(inout) :: fm2
1583 real (kind=kind_phys),
intent(inout) :: fh2
1584 real (kind=kind_phys),
intent(inout) :: fv
1587 character(len=*),
intent(inout) :: errmsg
1588 integer,
intent(inout) :: errflg
1592 real (kind=kind_phys),
intent(out) :: cm
1593 real (kind=kind_phys),
intent(out) :: ch
1594 real (kind=kind_phys),
intent(out) :: ch2
1597 real (kind=kind_phys) :: mozold
1598 real (kind=kind_phys) :: tmpcm
1599 real (kind=kind_phys) :: tmpch
1600 real (kind=kind_phys) :: mol
1601 real (kind=kind_phys) :: tvir
1602 real (kind=kind_phys) :: tmp1,tmp2,tmp3
1603 real (kind=kind_phys) :: fmnew
1604 real (kind=kind_phys) :: fhnew
1605 real (kind=kind_phys) :: moz2
1606 real (kind=kind_phys) :: tmpcm2
1607 real (kind=kind_phys) :: tmpch2
1608 real (kind=kind_phys) :: fm2new
1609 real (kind=kind_phys) :: fh2new
1610 real (kind=kind_phys) :: tmp12,tmp22,tmp32
1612 real (kind=kind_phys) :: cmfm, chfh, cm2fm2, ch2fh2
1620 if(zlvl <= zpd)
then
1621 write(*,*)
'critical glacier problem: zlvl <= zpd; model stops', zlvl, zpd
1624 errmsg =
"stop in noah-mp glacier"
1627 call wrf_error_fatal(
"stop in noah-mp glacier")
1631 tmpcm = log((zlvl-zpd) / z0m)
1632 tmpch = log((zlvl-zpd) / z0h)
1633 tmpcm2 = log((2.0 + z0m) / z0m)
1634 tmpch2 = log((2.0 + z0h) / z0h)
1642 tvir = (1. + 0.61*qair) * sfctmp
1643 tmp1 = vkc * (grav/tvir) * h/(rhoair*cpair)
1644 if (abs(tmp1) .le. mpe) tmp1 = mpe
1645 mol = -1. * fv**3 / tmp1
1646 moz = min( (zlvl-zpd)/mol, 1.)
1647 moz2 = min( (2.0 + z0h)/mol, 1.)
1652 if (mozold*moz .lt. 0.) mozsgn = mozsgn+1
1653 if (mozsgn .ge. 2)
then
1663 if (moz .lt. 0.)
then
1664 tmp1 = (1. - 16.*moz)**0.25
1665 tmp2 = log((1.+tmp1*tmp1)/2.)
1666 tmp3 = log((1.+tmp1)/2.)
1667 fmnew = 2.*tmp3 + tmp2 - 2.*atan(tmp1) + 1.5707963
1671 tmp12 = (1. - 16.*moz2)**0.25
1672 tmp22 = log((1.+tmp12*tmp12)/2.)
1673 tmp32 = log((1.+tmp12)/2.)
1674 fm2new = 2.*tmp32 + tmp22 - 2.*atan(tmp12) + 1.5707963
1692 fm = 0.5 * (fm+fmnew)
1693 fh = 0.5 * (fh+fhnew)
1694 fm2 = 0.5 * (fm2+fm2new)
1695 fh2 = 0.5 * (fh2+fh2new)
1700 fh = min(fh,0.9*tmpch)
1701 fm = min(fm,0.9*tmpcm)
1702 fh2 = min(fh2,0.9*tmpch2)
1703 fm2 = min(fm2,0.9*tmpcm2)
1709 if(abs(cmfm) <= mpe) cmfm = mpe
1710 if(abs(chfh) <= mpe) chfh = mpe
1711 if(abs(cm2fm2) <= mpe) cm2fm2 = mpe
1712 if(abs(ch2fh2) <= mpe) ch2fh2 = mpe
1713 cm = vkc*vkc/(cmfm*cmfm)
1714 ch = vkc*vkc/(cmfm*chfh)
1715 ch2 = vkc*vkc/(cm2fm2*ch2fh2)
1726 ssoil ,snowh ,zbot ,zsnso ,df , & !in
1738 integer,
intent(in) :: nsoil
1739 integer,
intent(in) :: nsnow
1740 integer,
intent(in) :: isnow
1742 real (kind=kind_phys),
intent(in) :: dt
1743 real (kind=kind_phys),
intent(in) :: tbot
1744 real (kind=kind_phys),
intent(in) :: ssoil
1745 real (kind=kind_phys),
intent(in) :: snowh
1746 real (kind=kind_phys),
intent(in) :: zbot
1747 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: zsnso
1748 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: df
1749 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: hcpct
1753 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: stc
1758 real (kind=kind_phys) :: zbotsno
1759 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: ai, bi, ci, rhsts
1760 real (kind=kind_phys) :: eflxb
1761 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: phi
1767 phi(isnow+1:nsoil) = 0.
1771 zbotsno = zbot - snowh
1776 stc ,tbot ,zbotsno ,df , &
1777 hcpct ,ssoil ,phi , &
1778 ai ,bi ,ci ,rhsts , &
1782 ai ,bi ,ci ,rhsts , &
1790 stc ,tbot ,zbot ,df , & !in
1791 hcpct ,ssoil ,phi , & !in
1792 ai ,bi ,ci ,rhsts , & !out
1804 integer,
intent(in) :: nsoil
1805 integer,
intent(in) :: nsnow
1806 integer,
intent(in) :: isnow
1807 real (kind=kind_phys),
intent(in) :: tbot
1808 real (kind=kind_phys),
intent(in) :: zbot
1810 real (kind=kind_phys),
intent(in) :: ssoil
1811 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: zsnso
1812 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: stc
1813 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: df
1814 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: hcpct
1815 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: phi
1819 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(out) :: rhsts
1820 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(out) :: ai
1821 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(out) :: bi
1822 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(out) :: ci
1823 real (kind=kind_phys),
intent(out) :: botflx
1828 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: ddz
1829 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: denom
1830 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: dtsdz
1831 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: eflux
1832 real (kind=kind_phys) :: temp1
1835 do k = isnow+1, nsoil
1836 if (k == isnow+1)
then
1837 denom(k) = - zsnso(k) * hcpct(k)
1838 temp1 = - zsnso(k+1)
1839 ddz(k) = 2.0 / temp1
1840 dtsdz(k) = 2.0 * (stc(k) - stc(k+1)) / temp1
1841 eflux(k) = df(k) * dtsdz(k) - ssoil - phi(k)
1842 else if (k < nsoil)
then
1843 denom(k) = (zsnso(k-1) - zsnso(k)) * hcpct(k)
1844 temp1 = zsnso(k-1) - zsnso(k+1)
1845 ddz(k) = 2.0 / temp1
1846 dtsdz(k) = 2.0 * (stc(k) - stc(k+1)) / temp1
1847 eflux(k) = (df(k)*dtsdz(k) - df(k-1)*dtsdz(k-1)) - phi(k)
1848 else if (k == nsoil)
then
1849 denom(k) = (zsnso(k-1) - zsnso(k)) * hcpct(k)
1850 temp1 = zsnso(k-1) - zsnso(k)
1851 if(opt_tbot == 1)
then
1854 if(opt_tbot == 2)
then
1855 dtsdz(k) = (stc(k) - tbot) / ( 0.5*(zsnso(k-1)+zsnso(k)) - zbot)
1856 botflx = -df(k) * dtsdz(k)
1858 eflux(k) = (-botflx - df(k-1)*dtsdz(k-1) ) - phi(k)
1862 do k = isnow+1, nsoil
1863 if (k == isnow+1)
then
1865 ci(k) = - df(k) * ddz(k) / denom(k)
1866 if (opt_stc == 1 .or. opt_stc == 3)
then
1869 if (opt_stc == 2)
then
1870 bi(k) = - ci(k) + df(k)/(0.5*zsnso(k)*zsnso(k)*hcpct(k))
1872 else if (k < nsoil)
then
1873 ai(k) = - df(k-1) * ddz(k-1) / denom(k)
1874 ci(k) = - df(k ) * ddz(k ) / denom(k)
1875 bi(k) = - (ai(k) + ci(k))
1876 else if (k == nsoil)
then
1877 ai(k) = - df(k-1) * ddz(k-1) / denom(k)
1879 bi(k) = - (ai(k) + ci(k))
1881 rhsts(k) = eflux(k)/ (-denom(k))
1889 ai ,bi ,ci ,rhsts , & !inout
1898 integer,
intent(in) :: nsoil
1899 integer,
intent(in) :: nsnow
1900 integer,
intent(in) :: isnow
1901 real (kind=kind_phys),
intent(in) :: dt
1904 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: ai
1905 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: bi
1906 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: ci
1907 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: stc
1908 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: rhsts
1912 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: rhstsin
1913 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: ciin
1916 do k = isnow+1,nsoil
1917 rhsts(k) = rhsts(k) * dt
1919 bi(k) = 1. + bi(k) * dt
1925 do k = isnow+1,nsoil
1926 rhstsin(k) = rhsts(k)
1932 call rosr12_glacier (ci,ai,bi,ciin,rhstsin,rhsts,isnow+1,nsoil,nsnow)
1936 do k = isnow+1,nsoil
1937 stc(k) = stc(k) + ci(k)
1964 integer,
intent(in) :: ntop
1965 integer,
intent(in) :: nsoil,nsnow
1968 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in):: a, b, d
1969 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout):: c,p,delta
1975 p(ntop) = - c(ntop) / b(ntop)
1979 delta(ntop) = d(ntop) / b(ntop)
1984 p(k) = - c(k) * ( 1.0 / (b(k) + a(k) * p(k -1)) )
1985 delta(k) = (d(k) - a(k)* delta(k -1))* (1.0/ (b(k) + a(k)&
1991 p(nsoil) = delta(nsoil)
1996 kk = nsoil - k + (ntop-1) + 1
1997 p(kk) = p(kk) * p(kk +1) + delta(kk)
2006 stc ,snice ,snliq ,sneqv ,snowh , & !inout
2007 smc ,sh2o , & !inout
2008 qmelt ,imelt ,ponding )
2016 integer,
intent(in) :: nsnow
2017 integer,
intent(in) :: nsoil
2018 integer,
intent(in) :: isnow
2019 real (kind=kind_phys),
intent(in) :: dt
2020 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: fact
2021 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: dzsnso
2025 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: stc
2026 real (kind=kind_phys),
dimension(-nsnow+1:0) ,
intent(inout) :: snice
2027 real (kind=kind_phys),
dimension(-nsnow+1:0) ,
intent(inout) :: snliq
2028 real (kind=kind_phys),
intent(inout) :: sneqv
2029 real (kind=kind_phys),
intent(inout) :: snowh
2030 real (kind=kind_phys),
dimension( 1:nsoil),
intent(inout) :: sh2o
2031 real (kind=kind_phys),
dimension( 1:nsoil),
intent(inout) :: smc
2034 real (kind=kind_phys),
intent(out) :: qmelt
2035 integer,
dimension(-nsnow+1:nsoil),
intent(out) :: imelt
2036 real (kind=kind_phys),
intent(out) :: ponding
2041 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: hm
2042 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: xm
2043 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: wmass0
2044 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: wice0
2045 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: wliq0
2046 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: mice
2047 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: mliq
2048 real (kind=kind_phys),
dimension(-nsnow+1:nsoil) :: heatr
2049 real (kind=kind_phys) :: temp1
2050 real (kind=kind_phys) :: propor
2051 real (kind=kind_phys) :: xmf
2071 wmass0(j) = mice(j) + mliq(j)
2075 if (mice(j) > 0. .and. stc(j) >= tfrz)
then
2078 if (mliq(j) > 0. .and. stc(j) < tfrz)
then
2087 if (imelt(j) > 0)
then
2088 hm(j) = (stc(j)-tfrz)/fact(j)
2092 if (imelt(j) == 1 .and. hm(j) < 0.)
then
2096 if (imelt(j) == 2 .and. hm(j) > 0.)
then
2100 xm(j) = hm(j)*dt/hfus
2105if (opt_gla == 2)
then
2107 if (isnow == 0 .and. sneqv > 0. .and. stc(1) >= tfrz)
then
2108 hm(1) = (stc(1)-tfrz)/fact(1)
2110 xm(1) = hm(1)*dt/hfus
2113 sneqv = max(0.,temp1-xm(1))
2114 propor = sneqv/temp1
2115 snowh = max(0.,propor * snowh)
2116 heatr(1) = hm(1) - hfus*(temp1-sneqv)/dt
2117 if (heatr(1) > 0.)
then
2118 xm(1) = heatr(1)*dt/hfus
2119 stc(1) = stc(1) + fact(1)*heatr(1)
2124 qmelt = max(0.,(temp1-sneqv))/dt
2126 ponding = temp1-sneqv
2134 if (imelt(j) > 0 .and. abs(hm(j)) > 0.)
then
2137 if (xm(j) > 0.)
then
2138 mice(j) = max(0., wice0(j)-xm(j))
2139 heatr(j) = hm(j) - hfus*(wice0(j)-mice(j))/dt
2140 else if (xm(j) < 0.)
then
2141 mice(j) = min(wmass0(j), wice0(j)-xm(j))
2142 heatr(j) = hm(j) - hfus*(wice0(j)-mice(j))/dt
2145 mliq(j) = max(0.,wmass0(j)-mice(j))
2147 if (abs(heatr(j)) > 0.)
then
2148 stc(j) = stc(j) + fact(j)*heatr(j)
2149 if (mliq(j)*mice(j)>0.) stc(j) = tfrz
2152 qmelt = qmelt + max(0.,(wice0(j)-mice(j)))/dt
2157if (opt_gla == 1)
then
2160 mliq(j) = sh2o(j) * dzsnso(j) * 1000.
2161 mice(j) = (smc(j) - sh2o(j)) * dzsnso(j) * 1000.
2170 wmass0(j) = mice(j) + mliq(j)
2174 if (mice(j) > 0. .and. stc(j) >= tfrz)
then
2177 if (mliq(j) > 0. .and. stc(j) < tfrz)
then
2182 if (isnow == 0 .and. sneqv > 0. .and. j == 1)
then
2183 if (stc(j) >= tfrz)
then
2192 if (imelt(j) > 0)
then
2193 hm(j) = (stc(j)-tfrz)/fact(j)
2197 if (imelt(j) == 1 .and. hm(j) < 0.)
then
2201 if (imelt(j) == 2 .and. hm(j) > 0.)
then
2205 xm(j) = hm(j)*dt/hfus
2210 if (isnow == 0 .and. sneqv > 0. .and. xm(1) > 0.)
then
2212 sneqv = max(0.,temp1-xm(1))
2213 propor = sneqv/temp1
2214 snowh = max(0.,propor * snowh)
2215 heatr(1) = hm(1) - hfus*(temp1-sneqv)/dt
2216 if (heatr(1) > 0.)
then
2217 xm(1) = heatr(1)*dt/hfus
2225 qmelt = max(0.,(temp1-sneqv))/dt
2227 ponding = temp1-sneqv
2233 if (imelt(j) > 0 .and. abs(hm(j)) > 0.)
then
2236 if (xm(j) > 0.)
then
2237 mice(j) = max(0., wice0(j)-xm(j))
2238 heatr(j) = hm(j) - hfus*(wice0(j)-mice(j))/dt
2239 else if (xm(j) < 0.)
then
2240 mice(j) = min(wmass0(j), wice0(j)-xm(j))
2241 heatr(j) = hm(j) - hfus*(wice0(j)-mice(j))/dt
2244 mliq(j) = max(0.,wmass0(j)-mice(j))
2246 if (abs(heatr(j)) > 0.)
then
2247 stc(j) = stc(j) + fact(j)*heatr(j)
2249 if (mliq(j)*mice(j)>0.) stc(j) = tfrz
2253 if (j > 0) xmf = xmf + hfus * (wice0(j)-mice(j))/dt
2256 qmelt = qmelt + max(0.,(wice0(j)-mice(j)))/dt
2267 if (any(stc(1:4) > tfrz) .and. any(stc(1:4) < tfrz))
then
2269 if ( stc(j) > tfrz )
then
2270 heatr(j) = (stc(j)-tfrz)/fact(j)
2272 if (j .ne. k .and. stc(k) < tfrz .and. heatr(j) > 0.1)
then
2273 heatr(k) = (stc(k)-tfrz)/fact(k)
2274 if (abs(heatr(k)) > heatr(j))
then
2275 heatr(k) = heatr(k) + heatr(j)
2276 stc(k) = tfrz + heatr(k)*fact(k)
2279 heatr(j) = heatr(j) + heatr(k)
2285 stc(j) = tfrz + heatr(j)*fact(j)
2292 if (any(stc(1:4) > tfrz) .and. any(stc(1:4) < tfrz))
then
2294 if ( stc(j) < tfrz )
then
2295 heatr(j) = (stc(j)-tfrz)/fact(j)
2297 if (j .ne. k .and. stc(k) > tfrz .and. heatr(j) < -0.1)
then
2298 heatr(k) = (stc(k)-tfrz)/fact(k)
2299 if (heatr(k) > abs(heatr(j)))
then
2300 heatr(k) = heatr(k) + heatr(j)
2301 stc(k) = tfrz + heatr(k)*fact(k)
2304 heatr(j) = heatr(j) + heatr(k)
2310 stc(j) = tfrz + heatr(j)*fact(j)
2317 if (any(stc(1:4) > tfrz) .and. any(mice(1:4) > 0.))
then
2319 if ( stc(j) > tfrz )
then
2320 heatr(j) = (stc(j)-tfrz)/fact(j)
2321 xm(j) = heatr(j)*dt/hfus
2323 if (j .ne. k .and. mice(k) > 0. .and. xm(j) > 0.1)
then
2324 if (mice(k) > xm(j))
then
2325 mice(k) = mice(k) - xm(j)
2326 xmf = xmf + hfus * xm(j)/dt
2330 xm(j) = xm(j) - mice(k)
2331 xmf = xmf + hfus * mice(k)/dt
2335 mliq(k) = max(0.,wmass0(k)-mice(k))
2338 heatr(j) = xm(j)*hfus/dt
2339 stc(j) = tfrz + heatr(j)*fact(j)
2346 if (any(stc(1:4) < tfrz) .and. any(mliq(1:4) > 0.))
then
2348 if ( stc(j) < tfrz )
then
2349 heatr(j) = (stc(j)-tfrz)/fact(j)
2350 xm(j) = heatr(j)*dt/hfus
2352 if (j .ne. k .and. mliq(k) > 0. .and. xm(j) < -0.1)
then
2353 if (mliq(k) > abs(xm(j)))
then
2354 mice(k) = mice(k) - xm(j)
2355 xmf = xmf + hfus * xm(j)/dt
2359 xm(j) = xm(j) + mliq(k)
2360 xmf = xmf - hfus * mliq(k)/dt
2364 mliq(k) = max(0.,wmass0(k)-mice(k))
2367 heatr(j) = xm(j)*hfus/dt
2368 stc(j) = tfrz + heatr(j)*fact(j)
2381 if(opt_gla == 1)
then
2382 sh2o(j) = mliq(j) / (1000. * dzsnso(j))
2383 sh2o(j) = max(0.0,min(1.0,sh2o(j)))
2385 elseif(opt_gla == 2)
then
2395 qvap ,qdew ,ficeold ,zsoil , & !in
2396 isnow ,snowh ,sneqv ,snice ,snliq ,stc , & !inout
2397 dzsnso ,sh2o ,sice ,ponding ,zsnso ,fsh , & !inout
2398 runsrf ,runsub ,qsnow ,ponding1 ,ponding2 ,qsnbot , & !out
2407 integer,
intent(in) :: nsnow
2408 integer,
intent(in) :: nsoil
2409 integer,
dimension(-nsnow+1:0) ,
intent(in) :: imelt
2410 real (kind=kind_phys),
intent(in) :: dt
2411 real (kind=kind_phys),
intent(in) :: prcp
2412 real (kind=kind_phys),
intent(in) :: sfctmp
2413 real (kind=kind_phys),
intent(inout) :: qvap
2414 real (kind=kind_phys),
intent(inout) :: qdew
2415 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(in) :: ficeold
2416 real (kind=kind_phys),
dimension( 1:nsoil),
intent(in) :: zsoil
2419 integer,
intent(inout) :: isnow
2420 real (kind=kind_phys),
intent(inout) :: snowh
2421 real (kind=kind_phys),
intent(inout) :: sneqv
2422 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(inout) :: snice
2423 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(inout) :: snliq
2424 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: stc
2425 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: dzsnso
2426 real (kind=kind_phys),
dimension( 1:nsoil),
intent(inout) :: sh2o
2427 real (kind=kind_phys),
dimension( 1:nsoil),
intent(inout) :: sice
2428 real (kind=kind_phys) ,
intent(inout) :: ponding
2429 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: zsnso
2430 real (kind=kind_phys) ,
intent(inout) :: fsh
2433 real (kind=kind_phys),
intent(out) :: runsrf
2434 real (kind=kind_phys),
intent(out) :: runsub
2435 real (kind=kind_phys),
intent(out) :: qsnow
2436 real (kind=kind_phys),
intent(out) :: ponding1
2437 real (kind=kind_phys),
intent(out) :: ponding2
2438 real (kind=kind_phys),
intent(out) :: qsnbot
2439 real (kind=kind_phys),
intent(out) :: fpice
2440 real (kind=kind_phys),
intent(out) :: esnow
2443 real (kind=kind_phys) :: qrain
2444 real (kind=kind_phys) :: qseva
2445 real (kind=kind_phys) :: qsdew
2446 real (kind=kind_phys) :: qsnfro
2447 real (kind=kind_phys) :: qsnsub
2448 real (kind=kind_phys) :: snowhin
2449 real (kind=kind_phys) :: snoflow
2450 real (kind=kind_phys) :: bdfall
2451 real (kind=kind_phys) :: replace
2452 real (kind=kind_phys),
dimension( 1:nsoil) :: sice_save
2453 real (kind=kind_phys),
dimension( 1:nsoil) :: sh2o_save
2471 if(opt_snf == 1 .or. opt_snf == 4)
then
2472 if(sfctmp > tfrz+2.5)
then
2475 if(sfctmp <= tfrz+0.5)
then
2477 else if(sfctmp <= tfrz+2.)
then
2478 fpice = 1.-(-54.632 + 0.2*sfctmp)
2485 if(opt_snf == 2)
then
2486 if(sfctmp >= tfrz+2.2)
then
2493 if(opt_snf == 3)
then
2494 if(sfctmp >= tfrz)
then
2505 bdfall = min(120.,67.92+51.25*exp((sfctmp-tfrz)/2.59))
2507 qrain = prcp * (1.-fpice)
2508 qsnow = prcp * fpice
2509 snowhin = qsnow/bdfall
2516 esnow = qsnsub*2.83e+6
2519 snowhin ,qsnow ,qsnfro ,qsnsub ,qrain , &
2521 isnow ,snowh ,sneqv ,snice ,snliq , &
2522 sh2o ,sice ,stc ,dzsnso ,zsnso , &
2524 qsnbot ,snoflow ,ponding1 ,ponding2)
2528 runsrf = (ponding+ponding1+ponding2)/dt
2531 runsrf = runsrf + qsnbot + qrain
2533 runsrf = runsrf + qsnbot
2537 if(opt_gla == 1)
then
2540 replace = replace + dzsnso(ilev)*(sice(ilev) - sice_save(ilev) + sh2o(ilev) - sh2o_save(ilev))
2542 replace = replace * 1000.0 / dt
2544 sice = min(1.0,sice_save)
2545 elseif(opt_gla == 2)
then
2553 if(opt_gla == 1)
then
2554 runsub = snoflow + replace
2555 elseif(opt_gla == 2)
then
2566 snowhin ,qsnow ,qsnfro ,qsnsub ,qrain , & !in
2567 ficeold ,zsoil , & !in
2568 isnow ,snowh ,sneqv ,snice ,snliq , & !inout
2569 sh2o ,sice ,stc ,dzsnso ,zsnso , & !inout
2571 qsnbot ,snoflow ,ponding1 ,ponding2)
2576 integer,
intent(in) :: nsnow
2577 integer,
intent(in) :: nsoil
2578 integer,
dimension(-nsnow+1:0) ,
intent(in) :: imelt
2579 real (kind=kind_phys),
intent(in) :: dt
2580 real (kind=kind_phys),
intent(in) :: sfctmp
2581 real (kind=kind_phys),
intent(in) :: snowhin
2582 real (kind=kind_phys),
intent(in) :: qsnow
2583 real (kind=kind_phys),
intent(inout) :: qsnfro
2584 real (kind=kind_phys),
intent(inout) :: qsnsub
2585 real (kind=kind_phys),
intent(in) :: qrain
2586 real (kind=kind_phys),
dimension(-nsnow+1:0) ,
intent(in) :: ficeold
2587 real (kind=kind_phys),
dimension( 1:nsoil),
intent(in) :: zsoil
2590 integer,
intent(inout) :: isnow
2591 real (kind=kind_phys),
intent(inout) :: snowh
2592 real (kind=kind_phys),
intent(inout) :: sneqv
2593 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(inout) :: snice
2594 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(inout) :: snliq
2595 real (kind=kind_phys),
dimension( 1:nsoil),
intent(inout) :: sh2o
2596 real (kind=kind_phys),
dimension( 1:nsoil),
intent(inout) :: sice
2597 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: stc
2598 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: dzsnso
2599 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: zsnso
2600 real (kind=kind_phys),
intent(inout) :: fsh
2603 real (kind=kind_phys),
intent(out) :: qsnbot
2604 real (kind=kind_phys),
intent(out) :: snoflow
2605 real (kind=kind_phys),
intent(out) :: ponding1
2606 real (kind=kind_phys),
intent(out) :: ponding2
2610 real (kind=kind_phys) :: bdsnow
2611 real (kind=kind_phys),
parameter :: mwd = 100.
2619 isnow ,snowh ,dzsnso ,stc ,snice , &
2624 snliq ,imelt ,ficeold, &
2628 isnow ,sh2o ,stc ,snice ,snliq , &
2629 dzsnso ,sice ,snowh ,sneqv , &
2633 isnow ,stc ,snice ,snliq ,dzsnso )
2638 do iz = -nsnow+1, isnow
2648 isnow ,dzsnso ,snowh ,sneqv ,snice , &
2649 snliq ,sh2o ,sice ,stc , &
2650 ponding1 ,ponding2 ,fsh , &
2655 if(sneqv > mwd .and. isnow /= 0)
then
2656 bdsnow = snice(0) / dzsnso(0)
2657 snoflow = (sneqv - mwd)
2658 snice(0) = snice(0) - snoflow
2659 dzsnso(0) = dzsnso(0) - snoflow/bdsnow
2660 snoflow = snoflow / dt
2669 sneqv = sneqv + snice(iz) + snliq(iz)
2670 snowh = snowh + dzsnso(iz)
2677 dzsnso(iz) = -dzsnso(iz)
2680 dzsnso(1) = zsoil(1)
2682 dzsnso(iz) = (zsoil(iz) - zsoil(iz-1))
2685 zsnso(isnow+1) = dzsnso(isnow+1)
2686 do iz = isnow+2 ,nsoil
2687 zsnso(iz) = zsnso(iz-1) + dzsnso(iz)
2690 do iz = isnow+1 ,nsoil
2691 dzsnso(iz) = -dzsnso(iz)
2699 isnow ,snowh ,dzsnso ,stc ,snice , & !inout
2709 integer,
intent(in) :: nsoil
2710 integer,
intent(in) :: nsnow
2711 real (kind=kind_phys),
intent(in) :: dt
2712 real (kind=kind_phys),
intent(in) :: qsnow
2713 real (kind=kind_phys),
intent(in) :: snowhin
2714 real (kind=kind_phys),
intent(in) :: sfctmp
2718 integer,
intent(inout) :: isnow
2719 real (kind=kind_phys),
intent(inout) :: snowh
2720 real (kind=kind_phys),
intent(inout) :: sneqv
2721 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: dzsnso
2722 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: stc
2723 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(inout) :: snice
2724 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(inout) :: snliq
2734 if(isnow == 0 .and. qsnow > 0.)
then
2735 snowh = snowh + snowhin * dt
2736 sneqv = sneqv + qsnow * dt
2741 if(isnow == 0 .and. qsnow>0. .and. snowh >= 0.05)
then
2746 stc(0) = min(273.16, sfctmp)
2753 if(isnow < 0 .and. newnode == 0 .and. qsnow > 0.)
then
2754 snice(isnow+1) = snice(isnow+1) + qsnow * dt
2755 dzsnso(isnow+1) = dzsnso(isnow+1) + snowhin * dt
2764 snliq ,imelt ,ficeold, & !in
2771 integer,
intent(in) :: nsoil
2772 integer,
intent(in) :: nsnow
2773 integer,
dimension(-nsnow+1:0) ,
intent(in) :: imelt
2774 real (kind=kind_phys),
intent(in) :: dt
2775 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(in) :: stc
2776 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(in) :: snice
2777 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(in) :: snliq
2778 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(in) :: ficeold
2781 integer,
intent(inout) :: isnow
2782 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: dzsnso
2785 real (kind=kind_phys),
parameter :: c2 = 21.e-3
2786 real (kind=kind_phys),
parameter :: c3 = 2.5e-6
2787 real (kind=kind_phys),
parameter :: c4 = 0.04
2788 real (kind=kind_phys),
parameter :: c5 = 2.0
2789 real (kind=kind_phys),
parameter :: dm = 100.0
2790 real (kind=kind_phys),
parameter :: eta0 = 1.8e+6
2792 real (kind=kind_phys) :: burden
2793 real (kind=kind_phys) :: ddz1
2794 real (kind=kind_phys) :: ddz2
2795 real (kind=kind_phys) :: ddz3
2796 real (kind=kind_phys) :: dexpf
2797 real (kind=kind_phys) :: td
2798 real (kind=kind_phys) :: pdzdtc
2799 real (kind=kind_phys) :: void
2800 real (kind=kind_phys) :: wx
2801 real (kind=kind_phys) :: bi
2802 real (kind=kind_phys),
dimension(-nsnow+1:0) :: fice
2811 wx = snice(j) + snliq(j)
2812 fice(j) = snice(j) / wx
2813 void = 1. - (snice(j)/denice + snliq(j)/denh2o) / dzsnso(j)
2816 if (void > 0.001 .and. snice(j) > 0.1)
then
2817 bi = snice(j) / dzsnso(j)
2818 td = max(0.,tfrz-stc(j))
2825 if (bi > dm) ddz1 = ddz1*exp(-46.0e-3*(bi-dm))
2829 if (snliq(j) > 0.01*dzsnso(j)) ddz1=ddz1*c5
2833 ddz2 = -(burden+0.5*wx)*exp(-0.08*td-c2*bi)/eta0
2837 if (imelt(j) == 1)
then
2838 ddz3 = max(0.,(ficeold(j) - fice(j))/max(1.e-6,ficeold(j)))
2846 pdzdtc = (ddz1 + ddz2 + ddz3)*dt
2847 pdzdtc = max(-0.5,pdzdtc)
2851 dzsnso(j) = dzsnso(j)*(1.+pdzdtc)
2852 dzsnso(j) = min(max(dzsnso(j),(snliq(j)+snice(j))/500.0),(snliq(j)+snice(j))/50.0)
2857 burden = burden + wx
2865 isnow ,sh2o ,stc ,snice ,snliq , & !inout
2866 dzsnso ,sice ,snowh ,sneqv , & !inout
2873 integer,
intent(in) :: nsnow
2874 integer,
intent(in) :: nsoil
2878 integer,
intent(inout) :: isnow
2879 real (kind=kind_phys),
dimension( 1:nsoil),
intent(inout) :: sh2o
2880 real (kind=kind_phys),
dimension( 1:nsoil),
intent(inout) :: sice
2881 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: stc
2882 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(inout) :: snice
2883 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(inout) :: snliq
2884 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: dzsnso
2885 real (kind=kind_phys),
intent(inout) :: sneqv
2886 real (kind=kind_phys),
intent(inout) :: snowh
2887 real (kind=kind_phys),
intent(inout) :: ponding1
2888 real (kind=kind_phys),
intent(inout) :: ponding2
2893 integer :: isnow_old
2896 real (kind=kind_phys) :: zwice
2897 real (kind=kind_phys) :: zwliq
2898 real (kind=kind_phys) :: dzmin(3)
2899 data dzmin /0.045, 0.05, 0.2/
2905 do j = isnow_old+1,0
2906 if (snice(j) <= .1)
then
2908 snliq(j+1) = snliq(j+1) + snliq(j)
2909 snice(j+1) = snice(j+1) + snice(j)
2911 if (isnow_old < -1)
then
2912 snliq(j-1) = snliq(j-1) + snliq(j)
2913 snice(j-1) = snice(j-1) + snice(j)
2915 ponding1 = ponding1 + snliq(j)
2927 if (j > isnow+1 .and. isnow < -1)
then
2928 do i = j, isnow+2, -1
2930 snliq(i) = snliq(i-1)
2931 snice(i) = snice(i-1)
2932 dzsnso(i)= dzsnso(i-1)
2941 if(sice(1) < 0.)
then
2942 sh2o(1) = sh2o(1) + sice(1)
2946 if(isnow ==0)
return
2954 sneqv = sneqv + snice(j) + snliq(j)
2955 snowh = snowh + dzsnso(j)
2956 zwice = zwice + snice(j)
2957 zwliq = zwliq + snliq(j)
2964 if (snowh < 0.05 .and. isnow < 0 )
then
2967 ponding2 = ponding2 + zwliq
2968 if(sneqv <= 0.) snowh = 0.
2980 if (isnow < -1)
then
2985 do i = isnow_old+1,0
2986 if (dzsnso(i) < dzmin(mssi))
then
2988 if (i == isnow+1)
then
2990 else if (i == 0)
then
2994 if ((dzsnso(i-1)+dzsnso(i)) < (dzsnso(i+1)+dzsnso(i))) neibor = i-1
2998 if (neibor > i)
then
3007 stc(j), dzsnso(l), snliq(l), snice(l), stc(l) )
3010 if (j-1 > isnow+1)
then
3011 do k = j-1, isnow+2, -1
3013 snice(k) = snice(k-1)
3014 snliq(k) = snliq(k-1)
3015 dzsnso(k) = dzsnso(k-1)
3021 if (isnow >= -1)
exit
3045 real (kind=kind_phys),
intent(in) :: dz2
3046 real (kind=kind_phys),
intent(in) :: wliq2
3047 real (kind=kind_phys),
intent(in) :: wice2
3048 real (kind=kind_phys),
intent(in) :: t2
3049 real (kind=kind_phys),
intent(inout) :: dz
3050 real (kind=kind_phys),
intent(inout) :: wliq
3051 real (kind=kind_phys),
intent(inout) :: wice
3052 real (kind=kind_phys),
intent(inout) :: t
3056 real (kind=kind_phys) :: dzc
3057 real (kind=kind_phys) :: wliqc
3058 real (kind=kind_phys) :: wicec
3059 real (kind=kind_phys) :: tc
3060 real (kind=kind_phys) :: h
3061 real (kind=kind_phys) :: h2
3062 real (kind=kind_phys) :: hc
3067 wicec = (wice+wice2)
3068 wliqc = (wliq+wliq2)
3069 h = (cice*wice+cwat*wliq) * (t-tfrz)+hfus*wliq
3070 h2= (cice*wice2+cwat*wliq2) * (t2-tfrz)+hfus*wliq2
3074 tc = tfrz + hc/(cice*wicec + cwat*wliqc)
3075 else if (hc.le.hfus*wliqc)
then
3078 tc = tfrz + (hc - hfus*wliqc) / (cice*wicec + cwat*wliqc)
3090 isnow ,stc ,snice ,snliq ,dzsnso )
3096 integer,
intent(in) :: nsnow
3097 integer,
intent(in) :: nsoil
3101 integer ,
intent(inout) :: isnow
3102 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: stc
3103 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(inout) :: snice
3104 real (kind=kind_phys),
dimension(-nsnow+1: 0),
intent(inout) :: snliq
3105 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: dzsnso
3111 real (kind=kind_phys) :: drr
3112 real (kind=kind_phys),
dimension( 1:nsnow) :: dz
3113 real (kind=kind_phys),
dimension( 1:nsnow) :: swice
3114 real (kind=kind_phys),
dimension( 1:nsnow) :: swliq
3115 real (kind=kind_phys),
dimension( 1:nsnow) :: tsno
3116 real (kind=kind_phys) :: zwice
3117 real (kind=kind_phys) :: zwliq
3118 real (kind=kind_phys) :: propor
3119 real (kind=kind_phys) :: dtdz
3123 if (j <= abs(isnow))
then
3124 dz(j) = dzsnso(j+isnow)
3125 swice(j) = snice(j+isnow)
3126 swliq(j) = snliq(j+isnow)
3127 tsno(j) = stc(j+isnow)
3135 if (dz(1) > 0.05)
then
3138 swice(1) = swice(1)/2.
3139 swliq(1) = swliq(1)/2.
3148 if (dz(1) > 0.05)
then
3151 zwice = propor*swice(1)
3152 zwliq = propor*swliq(1)
3154 swice(1) = propor*swice(1)
3155 swliq(1) = propor*swliq(1)
3158 call combo_glacier (dz(2), swliq(2), swice(2), tsno(2), drr, &
3159 zwliq, zwice, tsno(1))
3163 if (msno <= 2 .and. dz(2) > 0.10)
then
3165 dtdz = (tsno(1) - tsno(2))/((dz(1)+dz(2))/2.)
3167 swice(2) = swice(2)/2.
3168 swliq(2) = swliq(2)/2.
3172 tsno(3) = tsno(2) - dtdz*dz(2)/2.
3173 if (tsno(3) >= tfrz)
then
3176 tsno(2) = tsno(2) + dtdz*dz(2)/2.
3184 if (dz(2) > 0.2)
then
3187 zwice = propor*swice(2)
3188 zwliq = propor*swliq(2)
3190 swice(2) = propor*swice(2)
3191 swliq(2) = propor*swliq(2)
3193 call combo_glacier (dz(3), swliq(3), swice(3), tsno(3), drr, &
3194 zwliq, zwice, tsno(2))
3201 dzsnso(j) = dz(j-isnow)
3202 snice(j) = swice(j-isnow)
3203 snliq(j) = swliq(j-isnow)
3204 stc(j) = tsno(j-isnow)
3217 isnow ,dzsnso ,snowh ,sneqv ,snice , & !inout
3218 snliq ,sh2o ,sice ,stc , & !inout
3219 ponding1 ,ponding2 ,fsh , & !inout
3229 integer,
intent(in) :: nsnow
3230 integer,
intent(in) :: nsoil
3231 real (kind=kind_phys),
intent(in) :: dt
3232 real (kind=kind_phys),
intent(inout) :: qsnfro
3233 real (kind=kind_phys),
intent(inout) :: qsnsub
3234 real (kind=kind_phys),
intent(in) :: qrain
3238 real (kind=kind_phys),
intent(out) :: qsnbot
3242 integer,
intent(inout) :: isnow
3243 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: dzsnso
3244 real (kind=kind_phys),
intent(inout) :: snowh
3245 real (kind=kind_phys),
intent(inout) :: sneqv
3246 real (kind=kind_phys),
dimension(-nsnow+1:0),
intent(inout) :: snice
3247 real (kind=kind_phys),
dimension(-nsnow+1:0),
intent(inout) :: snliq
3248 real (kind=kind_phys),
dimension( 1:nsoil),
intent(inout) :: sh2o
3249 real (kind=kind_phys),
dimension( 1:nsoil),
intent(inout) :: sice
3250 real (kind=kind_phys),
dimension(-nsnow+1:nsoil),
intent(inout) :: stc
3251 real (kind=kind_phys),
intent(inout) :: ponding1
3252 real (kind=kind_phys),
intent(inout) :: ponding2
3253 real (kind=kind_phys),
intent(inout) :: fsh
3258 real (kind=kind_phys) :: qin
3259 real (kind=kind_phys) :: qout
3260 real (kind=kind_phys) :: wgdif
3261 real (kind=kind_phys),
dimension(-nsnow+1:0) :: vol_liq
3262 real (kind=kind_phys),
dimension(-nsnow+1:0) :: vol_ice
3263 real (kind=kind_phys),
dimension(-nsnow+1:0) :: epore
3264 real (kind=kind_phys) :: propor, temp
3269 if(sneqv == 0.)
then
3270 if(opt_gla == 1)
then
3271 sice(1) = sice(1) + (qsnfro-qsnsub)*dt/(dzsnso(1)*1000.)
3272 elseif(opt_gla == 2)
then
3273 fsh = fsh - (qsnfro-qsnsub)*hsub
3284 if(isnow == 0 .and. sneqv > 0.)
then
3285 if(opt_gla == 1)
then
3287 sneqv = sneqv - qsnsub*dt + qsnfro*dt
3289 snowh = max(0.,propor * snowh)
3290 elseif(opt_gla == 2)
then
3291 fsh = fsh - (qsnfro-qsnsub)*hsub
3297 sice(1) = sice(1) + sneqv/(dzsnso(1)*1000.)
3301 if(sice(1) < 0.)
then
3302 sh2o(1) = sh2o(1) + sice(1)
3307 if(snowh <= 1.e-8 .or. sneqv <= 1.e-6)
then
3314 if ( isnow < 0 )
then
3316 wgdif = snice(isnow+1) - qsnsub*dt + qsnfro*dt
3317 snice(isnow+1) = wgdif
3318 if (wgdif < 1.e-6 .and. isnow <0)
then
3320 isnow ,sh2o ,stc ,snice ,snliq , &
3321 dzsnso ,sice ,snowh ,sneqv , &
3322 ponding1, ponding2 )
3325 if ( isnow < 0 )
then
3326 snliq(isnow+1) = snliq(isnow+1) + qrain * dt
3327 snliq(isnow+1) = max(0., snliq(isnow+1))
3337 if (j >= isnow+1)
then
3338 vol_ice(j) = min(1., snice(j)/(dzsnso(j)*denice))
3339 epore(j) = 1. - vol_ice(j)
3340 vol_liq(j) = min(epore(j),snliq(j)/(dzsnso(j)*denh2o))
3350 if (j >= isnow+1)
then
3351 snliq(j) = snliq(j) + qin
3353 if (epore(j) < 0.05 .or. epore(j+1) < 0.05)
then
3356 qout = max(0.,(vol_liq(j)-ssi*epore(j))*dzsnso(j))
3357 qout = min(qout,(1.-vol_ice(j+1)-vol_liq(j+1))*dzsnso(j+1))
3360 qout = max(0.,(vol_liq(j) - ssi*epore(j))*dzsnso(j))
3363 snliq(j) = snliq(j) - qout
3377 fsh ,fgev ,ssoil ,sag ,prcp ,edir , &
3379 runsrf ,runsub ,sneqv ,dt ,beg_wb, errmsg , errflg )
3381 runsrf ,runsub ,sneqv ,dt ,beg_wb )
3389 integer ,
intent(in) :: iloc
3390 integer ,
intent(in) :: jloc
3391 real (kind=kind_phys) ,
intent(in) :: swdown
3392 real (kind=kind_phys) ,
intent(in) :: fsa
3393 real (kind=kind_phys) ,
intent(in) :: fsr
3394 real (kind=kind_phys) ,
intent(in) :: fira
3395 real (kind=kind_phys) ,
intent(in) :: fsh
3396 real (kind=kind_phys) ,
intent(in) :: fgev
3397 real (kind=kind_phys) ,
intent(in) :: ssoil
3398 real (kind=kind_phys) ,
intent(in) :: sag
3400 real (kind=kind_phys) ,
intent(in) :: prcp
3401 real (kind=kind_phys) ,
intent(in) :: edir
3402 real (kind=kind_phys) ,
intent(in) :: runsrf
3403 real (kind=kind_phys) ,
intent(in) :: runsub
3404 real (kind=kind_phys) ,
intent(in) :: sneqv
3405 real (kind=kind_phys) ,
intent(in) :: dt
3406 real (kind=kind_phys) ,
intent(in) :: beg_wb
3409 character(len=*) ,
intent(inout) :: errmsg
3410 integer ,
intent(inout) :: errflg
3413 real (kind=kind_phys) :: end_wb
3414 real (kind=kind_phys) :: errwat
3415 real (kind=kind_phys) :: erreng
3416 real (kind=kind_phys) :: errsw
3417 character(len=256) :: message
3419 errsw = swdown - (fsa + fsr)
3420 if (errsw > 0.01)
then
3421 write(*,*)
"sag =",sag
3422 write(*,*)
"fsa =",fsa
3423 write(*,*)
"fsr =",fsr
3424 write(message,*)
'errsw =',errsw
3427 errmsg = trim(message)//new_line(
'A')//
"radiation budget problem in noahmp glacier"
3430 call wrf_message(trim(message))
3431 call wrf_error_fatal(
"radiation budget problem in noahmp glacier")
3435 erreng = sag-(fira+fsh+fgev+ssoil)
3436 if(erreng > 0.01)
then
3437 write(message,*)
'erreng =',erreng
3439 errmsg = trim(message)
3441 call wrf_message(trim(message))
3443 write(message,
'(i6,1x,i6,1x,5f10.4)')iloc,jloc,sag,fira,fsh,fgev,ssoil
3446 errmsg = trim(errmsg)//new_line(
'A')//
"energy budget problem in noahmp glacier"
3449 call wrf_message(trim(message))
3450 call wrf_error_fatal(
"energy budget problem in noahmp glacier")
3455 errwat = end_wb-beg_wb-(prcp-edir-runsrf-runsub)*dt
3468 integer,
intent(in) :: iopt_alb
3469 integer,
intent(in) :: iopt_snf
3470 integer,
intent(in) :: iopt_tbot
3471 integer,
intent(in) :: iopt_stc
3473 integer,
intent(in) :: iopt_gla
3474 integer,
intent(in) :: iopt_sfc
3475 integer,
intent(in) :: iopt_trs
3481 opt_tbot = iopt_tbot
3495 use noahmp_glacier_globals
subroutine, public stability(z1, zvfun, gdx, tv1, thv1, wind, z0max, ztmax, tvs, grav, thsfc_loc, rb, fm, fh, fm10, fh2, cm, ch, stress, ustar)
subroutine, private atm_glacier(ep_2, epsm1, sfcprs, sfctmp, q2, soldn, cosz, thair, qair, eair, rhoair, solad, solai, swdown)
re-process atmospheric forcing
subroutine esat(t, esw, esi, desw, desi)
subroutine, private thermoprop_glacier(nsoil, nsnow, isnow, dzsnso, dt, snowh, snice, snliq, df, hcpct, snicev, snliqv, epore, fact)
calculate thermal properties of soil, snow, lake, and frozen soil.
subroutine, private combo_glacier(dz, wliq, wice, t, dz2, wliq2, wice2, t2)
subroutine, private sfcdif1_glacier(iter,zlvl,zpd,z0h,z0m, qair,sfctmp,h,rhoair,mpe,ur, ifdef ccpp
compute surface drag coefficient cm for momentum and ch for heat
subroutine, private rosr12_glacier(p, a, b, c, d, delta, ntop, nsoil, nsnow)
subroutine, private snowh2o_glacier(nsnow, nsoil, dt, qsnfro, qsnsub, qrain, isnow, dzsnso, snowh, sneqv, snice, snliq, sh2o, sice, stc, ponding1, ponding2, fsh, qsnbot)
subroutine, private hrt_glacier(nsnow, nsoil, isnow, zsnso, stc, tbot, zbot, df, hcpct, ssoil, phi, ai, bi, ci, rhsts, botflx)
subroutine, private tsnosoi_glacier(nsoil, nsnow, isnow, dt, tbot, ssoil, snowh, zbot, zsnso, df, hcpct, stc)
subroutine sfcdif4(iloc, jloc, ux, vx, t1d, p1d, psfcpa, pblhx, dx, znt, ep_1, ep_2, cp, itime, snwh, isice, psi_opt, tsk, qx, zlvl, iz0tlnd, qsfc, hfx, qfx, cm, chs, chs2, cqs2, rmolx, ust, rbx, fmx, fhx, stressx, fm10x, fh2x, wspdx, flhcx, flqcx)
subroutine, private compact_glacier(nsnow, nsoil, dt, stc, snice, snliq, imelt, ficeold, isnow, dzsnso)
subroutine, private snowalb_bats_glacier(nband, cosz, fage, albsnd, albsni)
subroutine, private error_glacier(iloc,jloc,swdown,fsa,fsr,fira, fsh,fgev,ssoil,sag,prcp,edir, ifdef ccpp
subroutine, private combine_glacier(nsnow, nsoil, isnow, sh2o, stc, snice, snliq, dzsnso, sice, snowh, sneqv, ponding1, ponding2)
subroutine, private radiation_glacier(dt, tg, sneqvo, sneqv, cosz, qsnow, solad, solai, albold, tauss, sag, fsr, fsa, albsnd, albsni)
Compute solar radiation: absorbed & reflected by the ground.
subroutine, private divide_glacier(nsnow, nsoil, isnow, stc, snice, snliq, dzsnso)
subroutine, private snowwater_glacier(nsnow, nsoil, imelt, dt, sfctmp, snowhin, qsnow, qsnfro, qsnsub, qrain, ficeold, zsoil, isnow, snowh, sneqv, snice, snliq, sh2o, sice, stc, dzsnso, zsnso, fsh, qsnbot, snoflow, ponding1, ponding2)
subroutine, private water_glacier(nsnow, nsoil, imelt, dt, prcp, sfctmp, qvap, qdew, ficeold, zsoil, isnow, snowh, sneqv, snice, snliq, stc, dzsnso, sh2o, sice, ponding, zsnso, fsh, runsrf, runsub, qsnow, ponding1, ponding2, qsnbot, fpice, esnow)
subroutine, private snowfall_glacier(nsoil, nsnow, dt, qsnow, snowhin, sfctmp, isnow, snowh, dzsnso, stc, snice, snliq, sneqv)
subroutine, private phasechange_glacier(nsnow, nsoil, isnow, dt, fact, dzsnso, stc, snice, snliq, sneqv, snowh, smc, sh2o, qmelt, imelt, ponding)
subroutine, public noahmp_glacier( iloc,jloc,cosz,nsnow,nsoil,dt, sfctmp,sfcprs,uu,vv,q2,soldn, prcp,lwdn,tbot,zlvl,ficeold,zsoil, thsfc_loc,prslkix,prsik1x,prslk1x, psfc,pblhx,iz0tlnd,itime, sigmaf1,garea1,psi_opt, ep_1,ep_2,epsm1,cp, qsnow,sneqvo,albold,cm,ch,isnow, sneqv,smc,zsnso,snowh,snice,snliq, tg,stc,sh2o,tauss,qsfc, fsa,fsr,fira,fsh,fgev,ssoil, trad,edir,runsrf,runsub,sag,albedo, qsnbot,ponding,ponding1,ponding2,t2m, q2e,z0h_total, ifdef ccpp
subroutine, private csnow_glacier(isnow, nsnow, nsoil, snice, snliq, dzsnso, tksno, cvsno, snicev, snliqv, epore)
snow bulk density, volumetric capacity, and thermal conductivity
subroutine, private energy_glacier(nsnow,nsoil,isnow,dt,qsnow,rhoair, eair,sfcprs,qair,sfctmp,lwdn,uu, vv,solad,solai,cosz,zref, tbot,zbot,zsnso,dzsnso,sigmaf1,garea1, thsfc_loc,prslkix,prsik1x,prslk1x, psfc,pblhx,iz0tlnd,itime,psi_opt, ep_1, ep_2, epsm1, cp, tg,stc,snowh,sneqv,sneqvo,sh2o, smc,snice,snliq,albold,cm,ch, ifdef ccpp
Compute energy budget (momentum & energy fluxes and phase changes).
subroutine albedo(parameters, vegtyp, ist, ice, nsoil, dt, cosz, fage, elai, esai, tg, tv, snowh, fsno, fwet, smc, sneqvo, sneqv, qsnow, fveg, iloc, jloc, albold, tauss, albgrd, albgri, albd, albi, fabd, fabi, ftdd, ftid, ftii, fsun, frevi, frevd, fregd, fregi, bgap, wgap, albsnd, albsni)
surface albedos. also fluxes (per unit incoming direct and diffuse radiation) reflected,...
subroutine, private snowalb_class_glacier(nband, qsnow, dt, alb, albold, albsnd, albsni)
subroutine, private snow_age_glacier(dt, tg, sneqvo, sneqv, tauss, fage)
subroutine, private hstep_glacier(nsnow, nsoil, isnow, dt, ai, bi, ci, rhsts, stc)
subroutine, private glacier_flux(nsoil,nsnow,emg,isnow,df,dzsnso,z0m, zlvl,zpd,qair,sfctmp,rhoair,sfcprs, ur,gamma,rsurf,lwdn,rhsur,smc, eair,stc,sag,snowh,lathea,sh2o, thsfc_loc,prslkix,prsik1x,prslk1x, psfc,pblhx,iz0tlnd,itime,uu,vv, sigmaf1,garea1,psi_opt,ep_1, ep_2, epsm1, cp, ifdef ccpp
use newton-raphson iteration to solve ground (tg) temperature that balances the surface energy budget...
subroutine, public noahmp_options_glacier(iopt_alb, iopt_snf, iopt_tbot, iopt_stc, iopt_gla, iopt_sfc, iopt_trs)
This module contains the CCPP-compliant GFS surface layer scheme.