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 )
748 subroutine snowdepth (GFLUX, ZN1, ZN2, time)
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 --------------"
819 subroutine calculationmo(GFLUX, CSRHOS, ZN1, TS1, ZN2, TS2, &
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