46 real(
rp) :: a0, a1, a2, c1, c2, c3
47 real(
rp) :: esat, qsat, deltaqsat, cp
50 integer,
parameter :: data_length_max=10000
59 real(
rp) :: albedomin = 0.4_rp
60 real(
rp) :: albedomax = 0.85_rp
61 real(
rp) :: zmin = 0.01_rp
64 real(
rp),
parameter :: epsilon = 0.97_rp
65 real(
rp),
parameter :: sigma = 5.67e-08_rp
67 real(
rp),
parameter :: ch = 0.002_rp
68 real(
rp),
parameter :: ce = 0.0021_rp
69 real(
rp),
parameter :: lv = 2.5e6_rp
70 real(
rp),
parameter :: lf = 3.34e5_rp
72 logical :: albedo_const = .true.
73 logical :: debug = .false.
75 integer :: zn_flag, ts_flag, sflag
85 real(
rp) :: snow_conductivity = 0.42_rp
86 real(
rp) :: water_content = 0.1_rp
87 real(
rp) :: snow_heat_capacityrho = 8.4e+05_rp
88 real(
rp) :: snow_rho = 400.0_rp
89 real(
rp) :: snowdepth_initial = 0.0_rp
90 real(
rp) :: albedo_value = 0.686_rp
92 namelist / param_land_phy_snow_ky90 / &
96 snow_heat_capacityrho, &
106 log_info(
"LAND_PHY_SNOW_KY90_setup",*)
'Setup'
110 read(
io_fid_conf,nml=param_land_phy_snow_ky90,iostat=ierr)
112 log_info(
"LAND_PHY_SNOW_KY90_setup",*)
'Not found namelist. Default used.'
113 elseif( ierr > 0 )
then
114 log_error(
"LAND_PHY_SNOW_KY90_setup",*)
'Not appropriate names in namelist PARAM_LAND_PHY_SNOW_KY90. Check!'
117 log_nml(param_land_phy_snow_ky90)
119 lambdas = snow_conductivity
121 csrhos = snow_heat_capacityrho
125 albedo = albedo_value
134 LIA, LIS, LIE, LJA, LJS, LJE, &
135 SFLX_water, SFLX_ENGI, & ! [IN]
136 PRSA, TA, QA, & ! [IN]
139 SFLX_RAD_dn, & ! [IN]
140 exists_land, dt, & ! [IN]
141 TSNOW, SWE, & ! [INOUT]
142 SDepth, SDzero, & ! [INOUT]
143 nosnowsec, & ! [INOUT]
146 SFLX_LH, SFLX_QV, & ! [OUT]
147 SFLX_QV_ENGI, & ! [OUT]
148 SFLX_GH, SNOW_LAND_GH, & ! [OUT]
149 SNOW_LAND_Water, & ! [OUT]
156 qsatf => atmos_saturation_psat_all
164 integer,
intent(in) :: lia, lis, lie
165 integer,
intent(in) :: lja, ljs, lje
168 real(
rp),
intent(in) :: sflx_water(lia,lja)
169 real(
rp),
intent(in) :: sflx_engi (lia,lja)
170 real(
rp),
intent(in) :: prsa (lia,lja)
171 real(
rp),
intent(in) :: ta (lia,lja)
172 real(
rp),
intent(in) :: wa (lia,lja)
173 real(
rp),
intent(in) :: ua (lia,lja)
174 real(
rp),
intent(in) :: va (lia,lja)
175 real(
rp),
intent(in) :: qa (lia,lja)
176 real(
rp),
intent(in) :: dens (lia,lja)
178 real(
dp),
intent(in) :: dt
179 logical,
intent(in) :: exists_land(lia,lja)
182 real(
rp),
intent(inout) :: tsnow (lia,lja)
183 real(
rp),
intent(inout) :: swe (lia,lja)
184 real(
rp),
intent(inout) :: sdepth (lia,lja)
185 real(
rp),
intent(inout) :: sdzero (lia,lja)
186 real(
rp),
intent(inout) :: nosnowsec (lia,lja)
189 real(
rp),
intent(out) :: salbedo (lia,lja,2)
190 real(
rp),
intent(out) :: sflx_sh (lia,lja)
191 real(
rp),
intent(out) :: sflx_lh (lia,lja)
192 real(
rp),
intent(out) :: sflx_qv (lia,lja)
193 real(
rp),
intent(out) :: sflx_qv_engi (lia,lja)
194 real(
rp),
intent(out) :: sflx_gh (lia,lja)
195 real(
rp),
intent(out) :: snow_land_gh (lia,lja)
196 real(
rp),
intent(out) :: snow_land_water(lia,lja)
197 real(
rp),
intent(out) :: snow_frac (lia,lja)
199 real(
rp) :: qcc (lia,lja)
200 real(
rp) :: qfusion (lia,lja)
201 real(
rp) :: melt (lia,lja)
202 real(
rp) :: swemelt (lia,lja)
210 real(
rp),
parameter :: uabs_min = 0.1_rp
214 real(
rp) :: qdry, psat
215 real(
rp) :: sflx_sw_dn, sflx_lw_dn
221 log_progress(*)
'Land / physics / snow / KY90'
226 if( ( exists_land(i,j) ) .and. &
227 ( swe(i,j)>0. .or. sflx_water(i,j)>0. ) )
then
229 uabs = sqrt( wa(i,j)**2 + ua(i,j)**2 + va(i,j)**2 )
234 call qsatf( ta(i,j), psat )
235 qasat = epsvap * psat / ( prsa(i,j) - ( 1.0_rp-epsvap ) * psat )
238 log_info(
"LAND_PHY_SNOW_KY90",*)
"RH, ",rh,
" DENS (1.289 org) ",dens(i,j)
250 znsnow1 = sdzero(i,j)
280 snow_land_gh(i,j) = snow_land_gh(i,j) / dt
281 snow_land_water(i,j) = swemelt(i,j) / dt
283 if ( swe1 <= 0. .and. swe(i,j) <= 0. )
then
284 snow_frac(i,j) = 0.0_rp
286 snow_frac(i,j) = 1.0_rp
292 sdzero(i,j) = znsnow1
295 log_info(
"LAND_PHY_SNOW_KY90",*)
"SNOW_frac, SWE, TSNOW", snow_frac(i,j), swe(i,j), tsnow(i,j)
300 snow_land_water(i,j) = sflx_water(i,j)
301 snow_frac(i,j) = 0.0_rp
307 salbedo(i,j,:) = 0.0_rp
308 sflx_sh(i,j) = 0.0_rp
309 sflx_lh(i,j) = 0.0_rp
310 sflx_gh(i,j) = 0.0_rp
311 sflx_qv(i,j) = 0.0_rp
313 qfusion(i,j) = 0.0_rp
315 swemelt(i,j) = 0.0_rp
316 snow_land_gh(i,j) = 0.0_rp
320 call file_history_in( qcc(:,:),
'LAND_SNOW_QCC',
'Heat used for changing temperature profile',
'J/m2', dim_type=
'XY' )
321 call file_history_in( qfusion(:,:),
'LAND_SNOW_QFUSION',
'Heat used for phase change of snow',
'J/m2', dim_type=
'XY' )
322 call file_history_in( melt(:,:),
'LAND_SNOW_MELT',
'Heat used for snow melt',
'J/m2', dim_type=
'XY' )
323 call file_history_in( swemelt(:,:),
'LAND_SNOW_SWEMELT',
'Equivalent water of melt snow',
'kg/m2', dim_type=
'XY' )
339 nosnowsec, & ! [INOUT]
340 ALBEDO_out, & ! [OUT]
343 LATENTFLUX, & ! [OUT]
351 Gflux2land, & ! [OUT]
369 real(RP),
intent(inout) :: TSNOW
370 real(RP),
intent(inout) :: SWE
371 real(RP),
intent(inout) :: DEPTH
372 real(RP),
intent(inout) :: ZNSNOW
374 real(RP),
intent(inout) :: nosnowsec
375 real(RP),
intent(out) :: ALBEDO_out
376 real(RP),
intent(out) :: Emiss
379 real(RP),
intent(out) :: HFLUX
380 real(RP),
intent(out) :: LATENTFLUX
381 real(RP),
intent(out) :: GFLUX
382 real(RP),
intent(out) :: EvapFLX
383 real(RP),
intent(out) :: Evap_ENGI
385 real(RP),
intent(out) :: QCC
386 real(RP),
intent(out) :: QFUSION
387 real(RP),
intent(out) :: MELT
388 real(RP),
intent(out) :: SWEMELT
389 real(RP),
intent(out) :: Gflux2land
393 real(RP),
intent(in) :: SFLX_SNOW
394 real(RP),
intent(in) :: SFLX_ENGI
395 real(RP),
intent(in) :: TA
396 real(RP),
intent(in) :: UA
397 real(RP),
intent(in) :: RH
398 real(RP),
intent(in) :: DENS
399 real(RP),
intent(in) :: SW
400 real(RP),
intent(in) :: LW
401 real(DP),
intent(in) :: time
415 real(RP) :: DELTADEPTH
434 snow = sflx_snow * time
444 depth0 = depth0 + (snow /
rhosnow)
445 znsnow0 = znsnow0 + (snow /
rhosnow)
448 log_info(
"SNOW_ky90_main",*)
"UA, SNOW,SFLX_SNOW,time : ", ua, snow, sflx_snow, time
449 log_info(
"SNOW_ky90_main",*)
"SWE , TSNOW, and TA : ", swe0, tsnow0, ta
450 log_info(
"SNOW_ky90_main",*)
"DEPTH is: ", depth0
451 log_info(
"SNOW_ky90_main",*)
"ZN beginning: ", znsnow0
456 if (snow > 0.0_rp)
then
459 nosnowsec = nosnowsec + time
461 if ( .not. albedo_const )
then
462 albedo = albedomin + (albedomax-albedomin)*exp(-1.0_rp*nosnowsec/real(4.0_rp*24.0_rp*3600.0_rp))
466 emiss = 1.0_rp - epsilon
468 log_info(
"SNOW_ky90_main",*)
"Albedo ",albedo
473 call groundflux (tsnow0, ta, ua, rh, dens, albedo, sw, lw, &
474 gflux, rflux, sflux, linflux, loutflux, hflux, latentflux)
477 evapflx = latentflux / lv
479 evap_engi = evapflx * (
cv_ice * tsnow0 -
lhf )
480 gflux = gflux + sflx_engi - evap_engi
482 swe0 = swe0 - evapflx * time
483 deltadepth = (evapflx * time) /
rhosnow
484 depth0 = depth0 - deltadepth
485 znsnow0 = znsnow0 - deltadepth
490 call check_allsnowmelt (gflux, tsnow0, znsnow0, depth0, sflag, time)
492 log_info(
"SNOW_ky90_main",*)
'LAND/snow: All snow melt'
493 qcc = 0.5_rp*csrhos*znsnow0*(t0-tsnow0)
502 gflux2land = gflux*time - (qcc + qfusion + melt)
506 call cal_param (znsnow0, tsnow0, gflux, ta, ua, rh, dens, lw, time)
509 call check_applicability (gflux, tsnow0, znsnow0, ta, ua, rh, dens, lw, gres, beta, time)
510 if ((gres >= 0.0_rp).and.(beta >= 0.0_rp))
then
511 log_info(
"SNOW_ky90_main",*)
'LAND/snow model is not appropriate',gres,beta
512 qcc = 0.5_rp*csrhos*znsnow0*(t0-tsnow0)
523 call snowdepth ( gflux, znsnow0, znsnow, time)
526 log_info(
"SNOW_ky90_main",*)
"ZN is: ", znsnow
528 IF(znsnow < zmin)
THEN
531 log_info(
"SNOW_ky90_main",*)
"ZN is replaced to: ", znsnow ,
" to ", zmin
534 ELSE IF(znsnow > depth0)
THEN
537 log_info(
"SNOW_ky90_main",*)
"ZN is replaced to: ", znsnow ,
" to ", depth0
543 call equation415(lambdas, c2, znsnow, rh, qsat, tsnow0, znsnow0, gflux, ta, ua, dens, lw, tsnow, time)
546 log_info(
"SNOW_ky90_main",*)
'TSNOW is: ', tsnow
553 call check_res(znsnow0, znsnow, tsnow0, tsnow, gflux, ta, ua, rh, dens, lw,
"1", time)
554 IF (znsnow < zmin)
THEN
557 log_info(
"SNOW_ky90_main",*)
"ZN is updated/replaced to: ", znsnow ,
" to ", zmin
560 ELSE IF(znsnow > depth0)
THEN
563 log_info(
"SNOW_ky90_main",*)
"ZN is updated/replaced to: ", znsnow ,
" to ", depth0
571 log_info(
"SNOW_ky90_main",*)
'TSNOW is updated: ', tsnow
572 log_info(
"SNOW_ky90_main",*)
'ZNSNOW is updated: ', znsnow
575 call check_res( znsnow0, znsnow, tsnow0, tsnow, gflux, ta, ua, rh, dens, lw,
"0", time)
578 IF ( (znsnow-zmin)/zmin < 0.00001 )
THEN
580 gflux, csrhos, znsnow0, tsnow0, znsnow, tsnow, &
581 melt, qcc, qfusion, time)
584 gflux, csrhos, znsnow0, tsnow0, znsnow, tsnow, &
585 melt, qcc, qfusion, time)
590 gflux2land = gflux*time - (qcc + qfusion + melt)
592 log_info(
"SNOW_ky90_main",*)
"### ZN_flag = ", zn_flag,
"TS_flag = ", ts_flag
593 log_info(
"SNOW_ky90_main",*)
'### Heat flux from snowpack to land surface: ', gflux2land
596 deltadepth = melt / ((1.0_rp-
w0)*lf*
rhosnow)
599 depth = depth0 - deltadepth
600 if (depth < znsnow)
then
603 log_info(
"SNOW_ky90_main",*)
"replace ZNSNOW <= DEPTH"
609 log_info(
"SNOW_ky90_main",*)
'MELT in water equivalent is: ', swemelt
610 log_info(
"SNOW_ky90_main",*)
'SWE0 is: ', swe
611 log_info(
"SNOW_ky90_main",*)
'DELTADEPTH is: ', deltadepth
612 log_info(
"SNOW_ky90_main",*)
'DEPTH0 is: ', depth
613 log_info(
"SNOW_ky90_main",*)
'ZNSNOW0 is: ', znsnow
623 subroutine groundflux (TS, TA, UA, RH, rhoair, ALPHA, SW, LW, &
624 GFLUX, RFLUX, SFLUX, LINFLUX, LOUTFLUX, HFLUX, LATENTFLUX)
627 real(RP),
intent(in) :: TS, TA, UA, RH, rhoair, ALPHA, SW, LW
628 real(RP),
intent(out) :: GFLUX, RFLUX, SFLUX, LINFLUX, LOUTFLUX, HFLUX, LATENTFLUX
630 esat = 0.6112_rp * exp( (17.67_rp * (ta-273.15_rp)) / (ta-29.66_rp) )
631 qsat = 0.622_rp * esat / 101.325_rp
632 deltaqsat = 0.622_rp * lv * qsat /(287.04_rp * (ta**2))
633 cp = 1004.67_rp * (1.0_rp + (0.84_rp * qsat))
636 sflux = (1.0_rp-alpha) * sw
637 rflux = epsilon * (lw-(sigma*(ts**4)))
638 linflux = epsilon * lw
639 loutflux = -1.0_rp * (epsilon * sigma * (ts**4))
640 hflux = cp * rhoair * ch * ua * (ts-ta)
641 latentflux = lv * rhoair * ce * ua * ((1-rh)*qsat + (deltaqsat*(ts-ta)))
643 gflux = (sflux + rflux - hflux - latentflux)
646 log_info(
"LAND_PHY_SNOW_KY90_groundflux",*)
"-------------- groundflux --------------"
647 log_info_cont(*)
"GFLUX is: ", gflux
648 log_info_cont(*)
"SFLUX: ", sflux
649 log_info_cont(*)
"RFLUX: ", rflux, linflux+loutflux
650 log_info_cont(*)
" (LONG in: ", linflux,
")"
651 log_info_cont(*)
" (LONG out: ", loutflux,
")"
652 log_info_cont(*)
"HFLUX is: ", hflux
653 log_info_cont(*)
"LATENT FLUX: ", latentflux
660 subroutine check_allsnowmelt (GFLUX, TS1, ZN1, D, sflag, time)
665 real(RP),
intent(in) :: GFLUX, TS1, ZN1, D
666 real(DP),
intent(in) ::time
667 integer,
intent(out) :: sflag
668 real(RP) :: energy_in, energy_use
669 real(RP) :: energy_use_ripe, energy_use_melt
671 energy_in = gflux * time
673 energy_use_ripe = 0.5_rp*csrhos*zn1*(t0-ts1) +
w0*
rhosnow*lf*zn1
674 energy_use_melt = (1.0_rp-
w0)*
rhosnow*lf*d
675 energy_use = energy_use_ripe + energy_use_melt
677 if(energy_in >= energy_use)
then
684 log_info(
"LAND_PHY_SNOW_KY90_check_allSnowMelt",*)
"Energy in =",energy_in
685 log_info(
"LAND_PHY_SNOW_KY90_check_allSnowMelt",*)
"Energy use =",energy_use
692 subroutine cal_param (ZN1, TS1, GFLUX, TA, UA, RH, rhoair, LW, time)
697 real(RP),
intent(in) :: ZN1, TS1, TA, UA, RH, rhoair, LW
698 real(RP),
intent(in) :: GFLUX
699 real(DP),
intent(in) :: time
702 c2 = (4.0_rp*epsilon*sigma*(ta**3)) + (cp*rhoair*ch*ua) + (lv*rhoair*ce*ua*deltaqsat)
705 a0 = lambdas*(c3*zn1 + (c1*zn1*(t0-ts1)- gflux*time))
706 a1 = c2*(c3*zn1 + c1*zn1*(t0-ts1) - gflux*time) - c3*lambdas
708 ( epsilon*(lw-(sigma*(ta**4))) + (c2*(ta-t0)) - (lv*rhoair*ce*ua*(1.0_rp-rh)*qsat) ) &
712 log_info(
"LAND_PHY_SNOW_KY90_cal_param",*)
"-------------- snowdepth --------------"
713 log_info_cont(*)
"C1",c1
714 log_info_cont(*)
"C2",c2,(4.0_rp*epsilon*sigma*(ta**3)),(cp*rhoair*ch*ua), (lv*rhoair*ce*ua*deltaqsat)
715 log_info_cont(*)
"C3",c3
716 log_info_cont(*)
"A0",a0
717 log_info_cont(*)
"A1",a1
718 log_info_cont(*)
"A2",a2
724 subroutine check_applicability (GFLUX, TS1, ZN1, TA, UA, RH, rhoair, LW, GFLUX_res, beta, time)
729 real(RP),
intent(in) :: GFLUX, TS1, ZN1
730 real(RP),
intent(in) :: TA, UA, RH, rhoair, LW
731 real(DP),
intent(in) :: time
732 real(RP),
intent(out) :: GFLUX_res, beta
733 real(RP) :: energy_in, energy_use_max
735 energy_in = gflux * time
736 energy_use_max = 0.5_rp*csrhos*zn1*(t0-ts1) +
w0*
rhosnow*lf*zn1
738 gflux_res = energy_in - energy_use_max
741 beta = (lambdas/c2)* &
742 ( t0 - ( ta*c2 - lv*rhoair*ce*ua*(1.0_rp-rh)*qsat + epsilon*(lw-sigma*(ta**4)))/c2 )
752 real(RP),
intent(in) :: GFLUX
753 real(RP),
intent(in) :: ZN1
754 real(DP),
intent(in) :: time
755 real(RP),
intent(out) :: ZN2
764 zn2 = ((-1.0_rp*a1) - ((a1**2.0_rp - 4.0_rp*a2*a0)**0.5_rp)) / (2.0_rp*a2)
767 log_info(
"LAND_PHY_SNOW_KY90_snowdepth",*)
"ZN old = ",zn1,
"ZN new = ",zn2
773 subroutine equation415(LAMBDAS, C2, ZN2, RH, QSAT, TS1, ZN1, GFLUX, TA, UA, rhoair, LW, TS2, time)
778 real(RP),
intent(in) :: TA, UA, rhoair, LW, ZN2, GFLUX, TS1,ZN1
779 real(RP),
intent(in) :: C2, LAMBDAS, RH, QSAT
780 real(DP),
intent(in) :: time
781 real(RP),
intent(out) :: TS2
784 ts2 = ((lambdas*t0) + (ta*c2*zn2) - (zn2*lv*rhoair*ce*ua*(1.0_rp-rh)*qsat) &
785 + zn2*epsilon*(lw - (sigma*(ta**4))))/ (lambdas + (c2*zn2))
788 log_info(
"LAND_PHY_SNOW_KY90_equation415",*)
"-------------- equation415 --------------"
790 ts_check = gflux*time/(0.5*csrhos*zn2) + t0 - zn1*(t0-ts1)/zn2 -
w0*
rhosnow*lf*(zn1-zn2)/(0.5*csrhos*zn2)
792 log_info_cont(*)
"compare ",ts2, ts_check, ts2-ts_check
797 end subroutine equation415
804 real(RP),
intent(in) :: ZN1
805 real(RP),
intent(in) :: TS, GFLUX
806 real(DP),
intent(in) :: time
807 real(RP),
intent(out) :: ZN2
809 zn2 = zn1 + ((c1*zn1*(t0-ts) - gflux*time) / c3)
812 log_info(
"LAND_PHY_SNOW_KY90_recalculateZ",*)
"-------------- recalculate Z --------------"
820 MELT, QCC, QFUSION, time)
826 real(RP),
intent(in) :: GFLUX, CSRHOS, ZN1, TS1, TS2, ZN2
827 real(DP),
intent(in) :: time
828 real(RP),
intent(out) :: MELT, QCC, QFUSION
830 qcc = 0.5_rp*csrhos*(zn1*(t0-ts1)-zn2*(t0-ts2))
832 melt = ( gflux*time - qcc - qfusion )
835 log_info(
"LAND_PHY_SNOW_KY90_calculationMO",*)
"--------------------------MELT----------------"
836 log_info_cont(*)
"GFLUX*time is: ", gflux*time
837 log_info_cont(*)
"QCC is : ", qcc
838 log_info_cont(*)
"QFUSION is : ", qfusion
839 log_info_cont(*)
"QMELT is : ", melt
841 log_info_cont(*) qcc+qfusion+melt
842 log_info_cont(*)
"diff= ", qcc + qfusion + melt - (gflux*time)
845 if ( abs(qcc+qfusion+melt - (gflux*time)) > 10.)
then
846 log_error(
"LAND_PHY_SNOW_KY90_calculationMO",*)
"Calculation is fault. Model would include bugs. Please check! Melt"
856 MELT, QCC, QFUSION, time)
860 real(RP),
intent(in) :: GFLUX, CSRHOS, ZN1, TS1, TS2, ZN2
861 real(DP),
intent(in) :: time
862 real(RP),
intent(out) :: MELT, QCC, QFUSION
864 qcc = 0.5_rp*csrhos*(zn1*(t0-ts1)-zn2*(t0-ts2))
869 log_info(
"LAND_PHY_SNOW_KY90_calculationNoMO",*)
"--------------------------NOMELT----------------"
870 log_info_cont(*)
"GFLUX*time is: ", gflux*time
871 log_info_cont(*)
"QCC is : ", qcc
872 log_info_cont(*)
"QFUSION is : ", qfusion
873 log_info_cont(*)
"QMELT is : ", melt
875 log_info_cont(*) qcc+qfusion+melt
876 log_info_cont(*)
"diff= ", qcc +qfusion - (gflux*time)
888 subroutine check_res(ZN1, ZN2, TS1, TS2, GFLUX, TA, UA, RH, rhoair, LW, flag, time)
892 real(RP),
intent(in) :: ZN1, ZN2, TS1, TS2, GFLUX
893 real(RP),
intent(in) :: TA, UA, RH, rhoair, LW
894 real(DP),
intent(in) :: time
895 character(len=*) :: flag
901 r1 = c1 * (zn1*(t0-ts1) - zn2*(t0-ts2)) + c3*(zn1-zn2)-gflux*time
902 r2 = zn2*( epsilon*lw-epsilon*sigma*(ta**4) - (c2*(ts2-ta)) - (lv*rhoair*ce*ua*(1.0_rp-rh)*qsat) ) &
904 r3 = ((lambdas*t0) + (ta*c2*zn2) - (zn2*lv*rhoair*ce*ua*(1.0_rp-rh)*qsat) &
905 + zn2*epsilon*( lw - (sigma*(ta**4))))/ (lambdas + (c2*zn2))-ts2
909 log_info(
"LAND_PHY_SNOW_KY90_check_res",*)
"R1 is : ", r1,
"flag = ", flag
910 log_info(
"LAND_PHY_SNOW_KY90_check_res",*)
"R2 is : ", r2, r3,
"flag = ", flag
912 if(abs(r1)>10000)
then
913 log_info(
"LAND_PHY_SNOW_KY90_check_res",*) c1 * (zn1*(t0-ts1) - zn2*(t0-ts2)), c3*(zn1-zn2), -gflux*time
921 subroutine cal_r1r2(ZN1, TS1, GFLUX, TA, UA, RH, rhoair, LW, time)
925 real(RP),
intent(in) :: ZN1, TS1
926 real(RP),
intent(in) :: GFLUX, TA, UA, RH, rhoair, LW
927 real(DP),
intent(in) :: time
935 real(RP) :: ts_r1,ts_r2,zn_r1,zn_r2
936 character(len=3) :: ttt =
""
938 real(RP) :: a,b,c,d,e,f,g
940 character(len=H_LONG) :: fname
951 call io_get_fname(fname,
'check_R1-R2_zn-base'//ttt//
'.dat')
952 open(fid, file=fname, status=
'unknown')
955 zn2 = zn0 + 9.0_rp*0.0001_rp*real((j-1),kind=rp)
958 ts_r1 = t0 - (c1*zn1*(t0-ts1) + c3*(zn1-zn2) - gflux*time)/(c1*zn2)
960 ts_r2 = ((lambdas*t0) + (ta*c2*zn2) - (zn2*lv*rhoair*ce*ua*(1.0_rp-rh)*qsat) &
961 + zn2*epsilon*( lw - (sigma*(ta**4))))/ (lambdas + (c2*zn2))
963 if ( debug )
write(fid,
'(3f15.5)') zn2, ts_r1, ts_r2
966 zn2 = -1.0_rp + 2.0_rp*0.00001_rp*real((j-1),kind=rp)
969 ts_r1 = t0 - (c1*zn1*(t0-ts1) + c3*(zn1-zn2) - gflux*time)/(c1*zn2)
971 ts_r2 = ((lambdas*t0) + (ta*c2*zn2) - (zn2*lv*rhoair*ce*ua*(1.0_rp-rh)*qsat) &
972 + zn2*epsilon*( lw - (sigma*(ta**4))))/ (lambdas + (c2*zn2))
974 if ( debug )
write(fid,
'(3f15.5)') zn2, ts_r1, ts_r2
977 zn2 = 1.0_rp + 9.0_rp*0.0001_rp*real((j-1),kind=rp)
980 ts_r1 = t0 - (c1*zn1*(t0-ts1) + c3*(zn1-zn2) - gflux*time)/(c1*zn2)
982 ts_r2 = ((lambdas*t0) + (ta*c2*zn2) - (zn2*lv*rhoair*ce*ua*(1.0_rp-rh)*qsat) &
983 + zn2*epsilon*( lw - (sigma*(ta**4))))/ (lambdas + (c2*zn2))
985 if ( debug )
write(70,
'(3f15.5)') zn2, ts_r1, ts_r2
992 b=-1.0_rp * (c1*zn1*(t0-ts1) + c3*zn1 - gflux*time)
996 e=( ta*c2 - (lv*rhoair*ce*ua*(1.0_rp-rh)*qsat) + epsilon*( lw - (sigma*(ta**4))))
997 f=(lambdas + (c2*zn2))
1001 call io_get_fname(fname,
'check_R1-R2-grad'//ttt//
'.dat')
1002 open(fid, file=fname, status=
'unknown')
1003 write(fid,
'(5f15.5)') b/c, d,e,lambdas,c2
1006 call io_get_fname(fname,
'check_R1-R2_ts-base'//ttt//
'.dat')
1007 open(fid, file=fname, status=
'unknown')
1009 ts2 = ts0 + 150.0_rp*0.00001_rp*real((i-1),kind=rp)
1011 zn_r1 = (c1*zn1*(t0-ts1) + c3*zn1 - gflux*time)/(c1*(t0-ts2)+c3)
1012 zn_r2 = -1.0_rp*lambdas*(t0-ts2)/(lw*epsilon-epsilon*sigma*(ta**4) - (c2*(ts2-ta)) - (lv*rhoair*ce*ua*(1.0_rp-rh)*qsat))
1014 write(fid,
'(3f15.5)') ts2, zn_r1, zn_r2
module atmosphere / hydrometeor
real(rp), public lhf
latent heat of fusion for use [J/kg]
real(rp), public cv_ice
CV for ice [J/kg/K].
module atmosphere / saturation
real(rp), public const_eps
small number
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
integer, public const_i_lw
long-wave radiation index
real(rp), public const_epsvap
Rdry / Rvap.
integer, public const_i_sw
short-wave radiation index
module coupler / surface-atmospehre
integer, parameter, public i_r_ir
integer, parameter, public i_r_direct
integer, parameter, public i_r_nir
integer, parameter, public n_rad_rgn
integer, parameter, public i_r_diffuse
integer, parameter, public n_rad_dir
integer, parameter, public i_r_vis
subroutine, public io_get_fname(outstr, instr, rank, ext, len)
generate process specific filename
integer, public io_fid_conf
Config file ID.
integer function, public io_get_available_fid()
search & get available file ID
module land / physics / snow / ky90
subroutine cal_r1r2(ZN1, TS1, GFLUX, TA, UA, RH, rhoair, LW, time)
subroutine, public land_phy_snow_ky90_setup
Setup.
subroutine cal_param(ZN1, TS1, GFLUX, TA, UA, RH, rhoair, LW, time)
subroutine snow_ky90_main(TSNOW, SWE, DEPTH, ZNSNOW, nosnowsec, ALBEDO_out, Emiss, HFLUX, LATENTFLUX, GFLUX, EvapFLX, Evap_ENGI, QCC, QFUSION, MELT, SWEMELT, Gflux2land, SFLX_SNOW, SFLX_ENGI, TA, UA, RH, DENS, SW, LW, time)
snow model main routine
subroutine snowdepth(GFLUX, ZN1, ZN2, time)
subroutine recalculatez(ZN1, TS, GFLUX, ZN2, time)
subroutine check_res(ZN1, ZN2, TS1, TS2, GFLUX, TA, UA, RH, rhoair, LW, flag, time)
subroutine check_applicability(GFLUX, TS1, ZN1, TA, UA, RH, rhoair, LW, GFLUX_res, beta, time)
subroutine, public land_phy_snow_ky90(LIA, LIS, LIE, LJA, LJS, LJE, SFLX_water, SFLX_ENGI, PRSA, TA, QA, WA, UA, VA, DENS, SFLX_RAD_dn, exists_land, dt, TSNOW, SWE, SDepth, SDzero, nosnowsec, Salbedo, SFLX_SH, SFLX_LH, SFLX_QV, SFLX_QV_ENGI, SFLX_GH, SNOW_LAND_GH, SNOW_LAND_Water, SNOW_frac)
Main routine for land submodel.
subroutine groundflux(TS, TA, UA, RH, rhoair, ALPHA, SW, LW, GFLUX, RFLUX, SFLUX, LINFLUX, LOUTFLUX, HFLUX, LATENTFLUX)
subroutine calculationmo(GFLUX, CSRHOS, ZN1, TS1, ZN2, TS2, MELT, QCC, QFUSION, time)
subroutine calculationnomo(GFLUX, CSRHOS, ZN1, TS1, ZN2, TS2, MELT, QCC, QFUSION, time)
subroutine, public prc_abort
Abort Process.
integer, parameter, public dp
integer, parameter, public rp