587SUBROUTINE flake_radflux ( depth_w, albedo_water, albedo_ice, albedo_snow, &
588 opticpar_water, opticpar_ice, opticpar_snow, &
629real(kind = kind_phys),
INTENT(IN) :: &
636TYPE (opticpar_medium),
INTENT(IN) :: &
643INTEGER (KIND = iintegers) :: &
650 IF(h_ice_p_flk.GE.h_ice_min_flk)
THEN
651 IF(h_snow_p_flk.GE.h_snow_min_flk)
THEN
652 i_snow_flk = i_atm_flk*(1.0-albedo_snow)
654 DO i=1, opticpar_snow%nband_optic
655 i_bot_flk = i_bot_flk + &
656 opticpar_snow%frac_optic(i)*exp(-opticpar_snow%extincoef_optic(i)*h_snow_p_flk)
658 i_ice_flk = i_snow_flk*i_bot_flk
660 i_snow_flk = i_atm_flk
661 i_ice_flk = i_atm_flk*(1.0-albedo_ice)
664 DO i=1, opticpar_ice%nband_optic
665 i_bot_flk = i_bot_flk + &
666 opticpar_ice%frac_optic(i)*exp(-opticpar_ice%extincoef_optic(i)*h_ice_p_flk)
668 i_w_flk = i_ice_flk*i_bot_flk
670 i_snow_flk = i_atm_flk
671 i_ice_flk = i_atm_flk
672 i_w_flk = i_atm_flk*(1.0-albedo_water)
675 IF(h_ml_p_flk.GE.h_ml_min_flk)
THEN
677 DO i=1, opticpar_water%nband_optic
678 i_bot_flk = i_bot_flk + &
679 opticpar_water%frac_optic(i)*exp(-opticpar_water%extincoef_optic(i)*h_ml_p_flk)
683 i_h_flk = i_w_flk*i_bot_flk
689 DO i=1, opticpar_water%nband_optic
690 i_bot_flk = i_bot_flk + &
691 opticpar_water%frac_optic(i)*exp(-opticpar_water%extincoef_optic(i)*depth_w)
693 i_bot_flk = i_w_flk*i_bot_flk
695 IF(h_ml_p_flk.GE.h_ml_min_flk)
THEN
697 DO i=1, opticpar_water%nband_optic
698 i_intm_0_h_flk = i_intm_0_h_flk + &
699 opticpar_water%frac_optic(i)/opticpar_water%extincoef_optic(i)* &
700 (1.0 - exp(-opticpar_water%extincoef_optic(i)*h_ml_p_flk))
702 i_intm_0_h_flk = i_w_flk*i_intm_0_h_flk/h_ml_p_flk
704 i_intm_0_h_flk = i_h_flk
707 IF(h_ml_p_flk.LE.depth_w-h_ml_min_flk)
THEN
709 DO i=1, opticpar_water%nband_optic
710 i_intm_h_d_flk = i_intm_h_d_flk + &
711 opticpar_water%frac_optic(i)/opticpar_water%extincoef_optic(i)* &
712 ( exp(-opticpar_water%extincoef_optic(i)*h_ml_p_flk) &
713 - exp(-opticpar_water%extincoef_optic(i)*depth_w) )
715 i_intm_h_d_flk = i_w_flk*i_intm_h_d_flk/(depth_w-h_ml_p_flk)
717 i_intm_h_d_flk = i_h_flk
722 IF(depth_bs.GE.h_ml_min_flk) then
724 DO i=1, opticpar_water%nband_optic
725 i_intm_d_h_flk = i_intm_d_h_flk + &
726 opticpar_water%frac_optic(i)/opticpar_water%extincoef_optic(i)* &
727 ( exp(-opticpar_water%extincoef_optic(i)*depth_w) &
728 - exp(-opticpar_water%extincoef_optic(i)*(depth_w+depth_bs)) )
730 i_intm_d_h_flk = i_w_flk*i_intm_d_h_flk/depth_bs
732 i_intm_d_h_flk = i_bot_flk
737 DO i=1, opticpar_water%nband_optic
738 i_hh_flk = i_hh_flk + &
739 opticpar_water%frac_optic(i)*exp(-opticpar_water%extincoef_optic(i)*(depth_w+depth_bs))
741 i_hh_flk = i_w_flk*i_hh_flk
755SUBROUTINE flake_main ( depthw, depthbs, T_bs, par_Coriolis, &
756 extincoef_water_typ, &
757 del_time, T_sfc_p, T_sfc_n, T_bot_2_in, &
806real(kind = kind_phys),
INTENT(IN) :: &
812 extincoef_water_typ , &
818real(kind = kind_phys) :: &
824real(kind = kind_phys),
INTENT(OUT) :: &
837INTEGER (KIND = iintegers) :: &
841real(kind = kind_phys) :: &
855real(kind = kind_phys) :: &
862real(kind = kind_phys) :: &
873real(kind = kind_phys) :: ct, ctt, cq
874real(kind = kind_phys) :: lake_depth_max
910t_snow_n_flk = t_snow_p_flk
911t_ice_n_flk = t_ice_p_flk
912t_wml_n_flk = t_wml_p_flk
913t_mnw_n_flk = t_mnw_p_flk
914t_bot_n_flk = t_bot_p_flk
915t_b1_n_flk = t_b1_p_flk
916h_snow_n_flk = h_snow_p_flk
917h_ice_n_flk = h_ice_p_flk
918h_ml_n_flk = h_ml_p_flk
919h_b1_n_flk = h_b1_p_flk
940IF(h_ice_p_flk.GE.h_ice_min_flk)
THEN
941 IF(h_ml_p_flk.LE.h_ml_min_flk)
THEN
942 q_w_flk = -tpl_kappa_w*(t_bot_p_flk-t_wml_p_flk)/depth_w
943 phi_t_pr0_flk = phi_t_pr0_1*c_t_p_flk-phi_t_pr0_2
944 q_w_flk = q_w_flk*max(phi_t_pr0_flk, 1.0)
953q_star_flk = q_w_flk + i_w_flk + i_h_flk - 2.0*i_intm_0_h_flk
963IF(lflk_botsed_use)
THEN
964 q_bot_flk = -tpl_kappa_w*(t_b1_p_flk-t_bot_p_flk)/max(h_b1_p_flk, h_b1_min_flk)*phi_b1_pr0
966 q_bot_flk = -tpl_kappa_w*(t_bot_2_in-t_bot_p_flk)/max(depth_bs, h_b1_min_flk)*phi_b1_pr0
994l_ice_create = .false.
995l_ice_meltabove = .false.
997ice_exist:
IF(h_ice_p_flk.LT.h_ice_min_flk)
THEN
999 l_ice_create = t_wml_p_flk.LE.(tpl_t_f+c_small_flk).AND.q_w_flk.LT.0.0
1000 IF(l_ice_create)
THEN
1001 d_h_ice_dt = -q_w_flk/tpl_rho_i/tpl_l_f
1002 h_ice_n_flk = h_ice_p_flk + d_h_ice_dt*del_time
1003 t_ice_n_flk = tpl_t_f + h_ice_n_flk*q_w_flk/tpl_kappa_i/phi_i_pr0_lin
1004 d_h_snow_dt = dmsnowdt_flk/tpl_rho_s_min
1005 h_snow_n_flk = h_snow_p_flk + d_h_snow_dt*del_time
1006 phi_i_pr1_flk = phi_i_pr1_lin &
1007 + phi_i_ast_mr*min(1.0, h_ice_n_flk/h_ice_max)
1008 r_h_icesnow = phi_i_pr1_flk/phi_s_pr0_lin*tpl_kappa_i/flake_snowheatconduct(h_snow_n_flk) &
1009 * h_snow_n_flk/max(h_ice_n_flk, h_ice_min_flk)
1010 t_snow_n_flk = t_ice_n_flk + r_h_icesnow*(t_ice_n_flk-tpl_t_f)
1015 l_snow_exists = h_snow_p_flk.GE.h_snow_min_flk
1017 melting:
IF(t_snow_p_flk.GE.(tpl_t_f-c_small_flk))
THEN
1019 IF(l_snow_exists)
THEN
1020 flk_str_1 = q_snow_flk + i_snow_flk - i_ice_flk
1021 IF(flk_str_1.GE.0.0)
THEN
1022 l_ice_meltabove = .true.
1023 d_h_snow_dt = (-flk_str_1/tpl_l_f+dmsnowdt_flk)/flake_snowdensity(h_snow_p_flk)
1024 d_h_ice_dt = -(i_ice_flk - i_w_flk - q_w_flk)/tpl_l_f/tpl_rho_i
1027 flk_str_1 = q_ice_flk + i_ice_flk - i_w_flk - q_w_flk
1028 IF(flk_str_1.GE.0.0)
THEN
1029 l_ice_meltabove = .true.
1030 d_h_ice_dt = -flk_str_1/tpl_l_f/tpl_rho_i
1031 d_h_snow_dt = dmsnowdt_flk/tpl_rho_s_min
1034 IF(l_ice_meltabove)
THEN
1035 h_ice_n_flk = h_ice_p_flk + d_h_ice_dt *del_time
1036 h_snow_n_flk = h_snow_p_flk + d_h_snow_dt*del_time
1037 t_ice_n_flk = tpl_t_f
1038 t_snow_n_flk = tpl_t_f
1043 no_melting:
IF(.NOT.l_ice_meltabove)
THEN
1045 d_h_snow_dt = flake_snowdensity(h_snow_p_flk)
1046 IF(d_h_snow_dt.LT.tpl_rho_s_max)
THEN
1047 flk_str_1 = h_snow_p_flk*tpl_gamma_rho_s/tpl_rho_w_r
1048 flk_str_1 = flk_str_1/(1.0-flk_str_1)
1052 d_h_snow_dt = dmsnowdt_flk/d_h_snow_dt/(1.0+flk_str_1)
1053 h_snow_n_flk = h_snow_p_flk + d_h_snow_dt*del_time
1055 phi_i_pr0_flk = h_ice_p_flk/h_ice_max
1056 c_i_flk = c_i_lin - c_i_mr*(1.0+phi_i_ast_mr)*phi_i_pr0_flk
1057 phi_i_pr1_flk = phi_i_pr1_lin + phi_i_ast_mr*phi_i_pr0_flk
1058 phi_i_pr0_flk = phi_i_pr0_lin - phi_i_pr0_flk
1060 h_ice_threshold = max(1.0, 2.0*c_i_flk*tpl_c_i*(tpl_t_f-t_ice_p_flk)/tpl_l_f)
1061 h_ice_threshold = phi_i_pr0_flk/c_i_flk*tpl_kappa_i/tpl_rho_i/tpl_c_i*h_ice_threshold
1062 h_ice_threshold = sqrt(h_ice_threshold*del_time)
1063 h_ice_threshold = min(0.9*h_ice_max, max(h_ice_threshold, h_ice_min_flk))
1066 IF(h_ice_p_flk.LT.h_ice_threshold)
THEN
1068 IF(l_snow_exists)
THEN
1069 flk_str_1 = q_snow_flk + i_snow_flk - i_w_flk
1071 flk_str_1 = q_ice_flk + i_ice_flk - i_w_flk
1073 d_h_ice_dt = -(flk_str_1-q_w_flk)/tpl_l_f/tpl_rho_i
1074 h_ice_n_flk = h_ice_p_flk + d_h_ice_dt *del_time
1075 t_ice_n_flk = tpl_t_f + h_ice_n_flk*flk_str_1/tpl_kappa_i/phi_i_pr0_flk
1079 d_h_ice_dt = tpl_kappa_i*(tpl_t_f-t_ice_p_flk)/h_ice_p_flk*phi_i_pr0_flk
1080 d_h_ice_dt = (q_w_flk+d_h_ice_dt)/tpl_l_f/tpl_rho_i
1081 h_ice_n_flk = h_ice_p_flk + d_h_ice_dt*del_time
1083 r_ti_icesnow = tpl_c_i*(tpl_t_f-t_ice_p_flk)/tpl_l_f
1084 r_tstar_icesnow = 1.0 - c_i_flk
1085 IF(l_snow_exists)
THEN
1086 r_h_icesnow = phi_i_pr1_flk/phi_s_pr0_lin*tpl_kappa_i/flake_snowheatconduct(h_snow_p_flk) &
1087 * h_snow_p_flk/h_ice_p_flk
1088 r_rho_c_icesnow = flake_snowdensity(h_snow_p_flk)*tpl_c_s/tpl_rho_i/tpl_c_i
1098 r_tstar_icesnow = r_tstar_icesnow*r_ti_icesnow
1104 flk_str_2 = q_snow_flk+i_snow_flk-i_w_flk
1105 flk_str_1 = c_i_flk*h_ice_p_flk + (1.0+c_s_lin*r_h_icesnow)*r_rho_c_icesnow*h_snow_p_flk
1106 d_t_ice_dt = -(1.0-2.0*c_s_lin)*r_h_icesnow*(tpl_t_f-t_ice_p_flk) &
1107 * tpl_c_s*dmsnowdt_flk
1109 r_tstar_icesnow = r_tstar_icesnow*r_ti_icesnow
1110 flk_str_2 = q_ice_flk+i_ice_flk-i_w_flk
1111 flk_str_1 = c_i_flk*h_ice_p_flk
1114 d_t_ice_dt = d_t_ice_dt + tpl_kappa_i*(tpl_t_f-t_ice_p_flk)/h_ice_p_flk*phi_i_pr0_flk &
1115 * (1.0-r_tstar_icesnow)
1116 d_t_ice_dt = d_t_ice_dt - r_tstar_icesnow*q_w_flk
1117 d_t_ice_dt = d_t_ice_dt + flk_str_2
1118 d_t_ice_dt = d_t_ice_dt/tpl_rho_i/tpl_c_i
1119 d_t_ice_dt = d_t_ice_dt/flk_str_1
1120 t_ice_n_flk = t_ice_p_flk + d_t_ice_dt*del_time
1123 phi_i_pr1_flk = min(1.0, h_ice_n_flk/h_ice_max)
1124 phi_i_pr1_flk = phi_i_pr1_lin + phi_i_ast_mr*phi_i_pr1_flk
1125 r_h_icesnow = phi_i_pr1_flk/phi_s_pr0_lin*tpl_kappa_i/flake_snowheatconduct(h_snow_n_flk) &
1126 *h_snow_n_flk/max(h_ice_n_flk, h_ice_min_flk)
1127 t_snow_n_flk = t_ice_n_flk + r_h_icesnow*(t_ice_n_flk-tpl_t_f)
1134h_ice_n_flk = min(h_ice_n_flk, h_ice_max)
1137t_snow_n_flk = min(t_snow_n_flk, tpl_t_f)
1138t_ice_n_flk = min(t_ice_n_flk, tpl_t_f)
1150 t_snow_n_flk = max(t_snow_n_flk, 184.0)
1151 t_ice_n_flk = max(t_ice_n_flk, 184.0)
1155IF(h_ice_n_flk.LT.h_ice_min_flk)
THEN
1157 t_ice_n_flk = tpl_t_f
1159 t_snow_n_flk = tpl_t_f
1160 l_ice_create = .false.
1161ELSE IF(h_snow_n_flk.LT.h_snow_min_flk)
THEN
1163 t_snow_n_flk = t_ice_n_flk
1171IF(l_ice_create) q_w_flk = 0.0
1172d_t_mnw_dt = (q_w_flk - q_bot_flk + i_w_flk - i_bot_flk)/tpl_rho_w_r/tpl_c_w/depth_w
1183t_mnw_n_flk = t_mnw_p_flk + d_t_mnw_dt*del_time
1185t_mnw_n_flk = max(t_mnw_n_flk, tpl_t_f)
1195htc_water:
IF(h_ice_n_flk.GE.h_ice_min_flk)
THEN
1197 t_mnw_n_flk = min(t_mnw_n_flk, tpl_t_r)
1199 t_wml_n_flk = tpl_t_f
1201 IF(l_ice_create)
THEN
1202 IF(h_ml_p_flk.GE.depth_w-h_ml_min_flk)
THEN
1206 h_ml_n_flk = h_ml_p_flk
1207 c_t_n_flk = c_t_p_flk
1209 t_bot_n_flk = t_wml_n_flk - (t_wml_n_flk-t_mnw_n_flk)/c_t_n_flk/(1.0-h_ml_n_flk/depth_w)
1212 ELSE IF(t_bot_p_flk.LT.tpl_t_r)
THEN
1213 h_ml_n_flk = h_ml_p_flk
1214 c_t_n_flk = c_t_p_flk
1215 t_bot_n_flk = t_wml_n_flk - (t_wml_n_flk-t_mnw_n_flk)/c_t_n_flk/(1.0-h_ml_n_flk/depth_w)
1219 t_bot_n_flk = tpl_t_r
1220 IF(h_ml_p_flk.GE.c_small_flk)
THEN
1221 c_t_n_flk = c_t_p_flk
1222 h_ml_n_flk = depth_w*(1.0-(t_wml_n_flk-t_mnw_n_flk)/(t_wml_n_flk-t_bot_n_flk)/c_t_n_flk)
1223 h_ml_n_flk = max(h_ml_n_flk, 0.0)
1225 h_ml_n_flk = h_ml_p_flk
1226 c_t_n_flk = (t_wml_n_flk-t_mnw_n_flk)/(t_wml_n_flk-t_bot_n_flk)
1227 c_t_n_flk = min(c_t_max, max(c_t_n_flk, c_t_min))
1231 t_bot_n_flk = min(t_bot_n_flk, tpl_t_r)
1236 flk_str_1 = flake_buoypar(t_wml_p_flk)*q_star_flk/tpl_rho_w_r/tpl_c_w
1237 IF(flk_str_1.LT.0.0)
THEN
1238 w_star_sfc_flk = (-flk_str_1*h_ml_p_flk)**(1.0/3.0)
1240 w_star_sfc_flk = 0.0
1249 conv_equil_h_scale = -q_w_flk/max(i_w_flk, c_small_flk)
1250 IF(conv_equil_h_scale.GT.0.0 .AND. conv_equil_h_scale.LT.1.0 &
1251 .AND. t_wml_p_flk.GT.tpl_t_r)
THEN
1252 conv_equil_h_scale = sqrt(6.0*conv_equil_h_scale) &
1253 + 2.0*conv_equil_h_scale/(1.0-conv_equil_h_scale)
1254 conv_equil_h_scale = min(depth_w, conv_equil_h_scale/extincoef_water_typ)
1257 conv_equil_h_scale = 0.0
1261 n_t_mean = flake_buoypar(0.5*(t_wml_p_flk+t_bot_p_flk))*(t_wml_p_flk-t_bot_p_flk)
1262 IF(h_ml_p_flk.LE.depth_w-h_ml_min_flk)
THEN
1263 tmp = n_t_mean/(depth_w-h_ml_p_flk)
1264 n_t_mean = sqrt(abs(tmp))
1271 d_c_t_dt = max(w_star_sfc_flk, u_star_w_flk, u_star_min_flk)**2_iintegers
1272 d_c_t_dt = n_t_mean*(depth_w-h_ml_p_flk)**2_iintegers &
1273 / c_relax_c/d_c_t_dt
1274 d_c_t_dt = (c_t_max-c_t_min)/max(d_c_t_dt, c_small_flk)
1279 c_tt_flk = c_tt_1*c_t_p_flk-c_tt_2
1280 c_q_flk = 2.0*c_tt_flk/c_t_p_flk
1282 mixing_regime:
IF(flk_str_1.LT.0.0)
THEN
1284 c_t_n_flk = c_t_p_flk + d_c_t_dt*del_time
1285 c_t_n_flk = min(c_t_max, max(c_t_n_flk, c_t_min))
1286 d_c_t_dt = (c_t_n_flk-c_t_p_flk)/del_time
1288 IF(h_ml_p_flk.LE.depth_w-h_ml_min_flk)
THEN
1289 IF(h_ml_p_flk.LE.h_ml_min_flk)
THEN
1290 d_h_ml_dt = c_cbl_1/c_cbl_2*max(w_star_sfc_flk, c_small_flk)
1299 r_h_icesnow = depth_w/h_ml_p_flk
1300 r_rho_c_icesnow = r_h_icesnow-1.0
1301 r_ti_icesnow = c_t_p_flk/c_tt_flk
1302 r_tstar_icesnow = (r_ti_icesnow/2.0-1.0)*r_rho_c_icesnow + 1.0
1303 d_h_ml_dt = -q_star_flk*(r_tstar_icesnow*(1.0+c_cbl_1)-1.0) - q_bot_flk
1304 d_h_ml_dt = d_h_ml_dt/tpl_rho_w_r/tpl_c_w
1305 flk_str_2 = (depth_w-h_ml_p_flk)*(t_wml_p_flk-t_bot_p_flk)*c_tt_2/c_tt_flk*d_c_t_dt
1306 d_h_ml_dt = d_h_ml_dt + flk_str_2
1307 flk_str_2 = i_bot_flk + (r_ti_icesnow-1.0)*i_h_flk - r_ti_icesnow*i_intm_h_d_flk
1308 flk_str_2 = flk_str_2 + (r_ti_icesnow-2.0)*r_rho_c_icesnow*(i_h_flk-i_intm_0_h_flk)
1309 flk_str_2 = flk_str_2/tpl_rho_w_r/tpl_c_w
1310 d_h_ml_dt = d_h_ml_dt + flk_str_2
1311 flk_str_2 = -c_cbl_2*r_tstar_icesnow*q_star_flk/tpl_rho_w_r/tpl_c_w/max(w_star_sfc_flk, c_small_flk)
1312 flk_str_2 = flk_str_2 + c_t_p_flk*(t_wml_p_flk-t_bot_p_flk)
1313 d_h_ml_dt = d_h_ml_dt/flk_str_2
1350 d_h_ml_dt = max(d_h_ml_dt, c_small_flk)
1351 h_ml_n_flk = h_ml_p_flk + d_h_ml_dt*del_time
1352 h_ml_n_flk = max(h_ml_min_flk, min(h_ml_n_flk, depth_w))
1354 h_ml_n_flk = depth_w
1359 d_h_ml_dt = max(u_star_w_flk, u_star_min_flk)
1360 zm_h_scale = (abs(par_coriolis)/c_sbl_zm_n + n_t_mean/c_sbl_zm_i)*d_h_ml_dt**2_iintegers
1361 zm_h_scale = zm_h_scale + flk_str_1/c_sbl_zm_s
1362 zm_h_scale = max(zm_h_scale, c_small_flk)
1363 zm_h_scale = d_h_ml_dt**3_iintegers/zm_h_scale
1364 zm_h_scale = max(h_ml_min_flk, min(zm_h_scale, h_ml_max_flk))
1365 zm_h_scale = max(zm_h_scale, conv_equil_h_scale)
1375 d_h_ml_dt = c_relax_h*d_h_ml_dt/zm_h_scale*del_time
1376 h_ml_n_flk = zm_h_scale - (zm_h_scale-h_ml_p_flk)*exp(-d_h_ml_dt)
1377 h_ml_n_flk = max(h_ml_min_flk, min(h_ml_n_flk, depth_w))
1378 d_h_ml_dt = (h_ml_n_flk-h_ml_p_flk)/del_time
1380 IF(h_ml_n_flk.LE.h_ml_p_flk) &
1381 d_c_t_dt = -d_c_t_dt
1382 c_t_n_flk = c_t_p_flk + d_c_t_dt*del_time
1383 c_t_n_flk = min(c_t_max, max(c_t_n_flk, c_t_min))
1384 d_c_t_dt = (c_t_n_flk-c_t_p_flk)/del_time
1395 END IF mixing_regime
1401 IF(h_ml_n_flk.LE.depth_w-h_ml_min_flk)
THEN
1403 IF(h_ml_n_flk.GT.h_ml_p_flk)
THEN
1404 r_h_icesnow = h_ml_p_flk/depth_w
1405 r_rho_c_icesnow = 1.0-r_h_icesnow
1406 r_ti_icesnow = 0.5*c_t_p_flk*r_rho_c_icesnow+c_tt_flk*(2.0*r_h_icesnow-1.0)
1407 r_tstar_icesnow = (0.5+c_tt_flk-c_q_flk)/r_ti_icesnow
1408 r_ti_icesnow = (1.0-c_t_p_flk*r_rho_c_icesnow)/r_ti_icesnow
1410 d_t_bot_dt = (q_w_flk-q_bot_flk+i_w_flk-i_bot_flk)/tpl_rho_w_r/tpl_c_w
1411 d_t_bot_dt = d_t_bot_dt - c_t_p_flk*(t_wml_p_flk-t_bot_p_flk)*d_h_ml_dt
1412 d_t_bot_dt = d_t_bot_dt*r_tstar_icesnow/depth_w
1414 flk_str_2 = i_intm_h_d_flk - (1.0-c_q_flk)*i_h_flk - c_q_flk*i_bot_flk
1415 flk_str_2 = flk_str_2*r_ti_icesnow/(depth_w-h_ml_p_flk)/tpl_rho_w_r/tpl_c_w
1416 d_t_bot_dt = d_t_bot_dt + flk_str_2
1418 flk_str_2 = (1.0-c_tt_2*r_ti_icesnow)/c_t_p_flk
1419 flk_str_2 = flk_str_2*(t_wml_p_flk-t_bot_p_flk)*d_c_t_dt
1420 d_t_bot_dt = d_t_bot_dt + flk_str_2
1426 t_bot_n_flk = t_bot_p_flk + d_t_bot_dt*del_time
1427 t_bot_n_flk = max(t_bot_n_flk, tpl_t_f)
1428 flk_str_2 = (t_bot_n_flk-tpl_t_r)*flake_buoypar(t_mnw_n_flk)
1429 IF(flk_str_2.LT.0.0) t_bot_n_flk = tpl_t_r
1430 t_wml_n_flk = c_t_n_flk*(1.0-h_ml_n_flk/depth_w)
1431 t_wml_n_flk = (t_mnw_n_flk-t_bot_n_flk*t_wml_n_flk)/(1.0-t_wml_n_flk)
1432 t_wml_n_flk = max(t_wml_n_flk, tpl_t_f)
1437 h_ml_n_flk = depth_w
1438 t_wml_n_flk = t_mnw_n_flk
1439 t_bot_n_flk = t_mnw_n_flk
1453use_sediment:
IF(lflk_botsed_use)
THEN
1455 IF(h_b1_p_flk.GE.depth_bs-h_b1_min_flk)
THEN
1457 t_b1_p_flk = t_bot_p_flk
1460 flk_str_1 = 2.0*phi_b1_pr0/(1.0-c_b1)*tpl_kappa_w/tpl_rho_w_r/tpl_c_w*del_time
1461 h_ice_threshold = sqrt(flk_str_1)
1462 h_ice_threshold = min(0.9*depth_bs, h_ice_threshold)
1463 flk_str_2 = c_b2/(1.0-c_b2)*(t_bs-t_b1_p_flk)/(depth_bs-h_b1_p_flk)
1465 IF(h_b1_p_flk.LT.h_ice_threshold)
THEN
1466 h_b1_n_flk = sqrt(h_b1_p_flk**2_iintegers+flk_str_1)
1467 d_h_b1_dt = (h_b1_n_flk-h_b1_p_flk)/del_time
1469 flk_str_1 = (q_bot_flk+i_bot_flk)/h_b1_p_flk/tpl_rho_w_r/tpl_c_w
1470 flk_str_1 = flk_str_1 - (1.0-c_b1)*(t_bot_n_flk-t_bot_p_flk)/del_time
1471 d_h_b1_dt = (1.0-c_b1)*(t_bot_p_flk-t_b1_p_flk)/h_b1_p_flk + c_b1*flk_str_2
1472 d_h_b1_dt = flk_str_1/d_h_b1_dt
1473 h_b1_n_flk = h_b1_p_flk + d_h_b1_dt*del_time
1475 d_t_b1_dt = flk_str_2*d_h_b1_dt
1476 t_b1_n_flk = t_b1_p_flk + d_t_b1_dt*del_time
1501 l_snow_exists = h_b1_n_flk.GE.depth_bs-h_b1_min_flk &
1502 .OR. h_b1_n_flk.LT.h_b1_min_flk &
1503 .OR.(t_bot_n_flk-t_b1_n_flk)*(t_bs-t_b1_n_flk).LE.0.0
1504 IF(l_snow_exists)
THEN
1505 h_b1_n_flk = depth_bs
1518 if ( abs(t_bot_p_flk-t_bot_2_in) .lt. 0.01 )
then
1521 t_bot_2_out = t_bot_2_in
1525 ctt = 11.0/18.0 * ct - 7.0/45.0
1529 flk_str_1 = ( q_bot_flk * cq + i_bot_flk - i_intm_d_h_flk/depth_bs ) /tpl_rho_w_r/tpl_c_w
1530 flk_str_1 = flk_str_1 - depth_bs * (0.5-ctt) * (t_bot_n_flk-t_bot_p_flk)/del_time
1531 flk_str_1 = flk_str_1 - ctt/ct*( (q_bot_flk+i_bot_flk-i_hh_flk)/tpl_rho_w_r/tpl_c_w - &
1532 depth_bs * ( 1.0 - ct ) * (t_bot_n_flk-t_bot_p_flk)/del_time )
1533 flk_str_2 = ctt * (t_bot_p_flk-t_bot_2_in)
1534 if(abs(flk_str_2)<0.01)
then
1537 d_h_d_dt = flk_str_1/flk_str_2
1541 flk_str_1 = (q_bot_flk+i_bot_flk-i_hh_flk)/tpl_rho_w_r/tpl_c_w
1542 flk_str_1 = flk_str_1 - depth_bs*(1.0-ct)*(t_bot_n_flk-t_bot_p_flk)/del_time
1543 flk_str_1 = flk_str_1 - ct * (t_bot_p_flk-t_bot_2_in) * d_h_d_dt
1544 flk_str_2 = ct * depth_bs
1545 d_t_h_dt = flk_str_1/flk_str_2
1548 t_bot_2_out = t_bot_2_in + d_t_h_dt * del_time
1550 if ( (t_bot_2_out-tpl_t_r)*(t_bot_n_flk-tpl_t_r) .le. 0.0 )
then
1551 t_bot_2_out = tpl_t_r
1552 else if ( (t_bot_n_flk-tpl_t_r)/(t_bot_2_out-tpl_t_r) .le. 1.0 )
then
1553 t_bot_2_out = t_bot_n_flk
1557 flk_str_2 = depth_w + depth_bs
1559 depth_w = depth_w - d_h_d_dt * del_time
1561 depth_w = max(max(lake_depth_max,0.3*flk_str_2), &
1562 min(depth_w, min(0.4*flk_str_2, 80.0)))
1563 depth_bs = flk_str_2 - depth_w
1574flk_str_2 = (t_wml_n_flk-t_bot_n_flk)*flake_buoypar(t_mnw_n_flk)
1575IF(flk_str_2.LT.0.0)
THEN
1587 h_ml_n_flk = depth_w
1588 t_wml_n_flk = t_mnw_n_flk
1589 t_bot_n_flk = t_mnw_n_flk
1600IF(h_snow_n_flk.GE.h_snow_min_flk)
THEN
1601 t_sfc_n = t_snow_n_flk
1602ELSE IF(h_ice_n_flk.GE.h_ice_min_flk)
THEN
1603 t_sfc_n = t_ice_n_flk
1605 t_sfc_n = t_wml_n_flk
1607if(t_sfc_n.lt.200.0 .or. t_sfc_n.gt.350.0)
then
1608 print*,
'h_snow_n_flk=',h_snow_n_flk,
' h_Snow_min_flk=',h_snow_min_flk
1609 print*,
'h_ice_n_flk=',h_ice_n_flk,
' h_Ice_min_flk= ',h_ice_min_flk
1610 print*,
'T_wML_n_flk=',t_wml_n_flk
1611 print*,
'T_sfc_n=',t_sfc_n
2171SUBROUTINE sfcflx_momsenlat ( height_u, height_tq, fetch, &
2172 U_a, T_a, q_a, T_s, P_a, h_ice, &
2173 Q_momentum, Q_sensible, Q_latent, Q_watvap,&
2212real(kind = kind_phys),
INTENT(IN) :: &
2225real(kind = kind_phys),
INTENT(OUT) :: &
2233INTEGER (KIND = iintegers),
PARAMETER :: &
2242INTEGER (KIND = iintegers) :: &
2247real(kind = kind_phys) :: &
2253real(kind = kind_phys) :: &
2265real(kind = kind_phys) :: &
2299wvpres_s = sfcflx_satwvpres(t_s, h_ice)
2300q_s = sfcflx_spechum(wvpres_s, p_a)
2301rho_a = sfcflx_rhoair(t_s, q_s, p_a)
2308q_mom_mol = -tpsf_nu_u_a*u_a/height_u
2309q_sen_mol = -tpsf_kappa_t_a*(t_a-t_s)/height_tq
2310q_lat_mol = -tpsf_kappa_q_a*(q_a-q_s)/height_tq
2316par_conv_visc = (t_s-t_a)/t_s*sqrt(tpsf_kappa_t_a) + (q_s-q_a)*tpsf_alpha_q*sqrt(tpsf_kappa_q_a)
2317IF(par_conv_visc.GT.0.0)
THEN
2318 l_conv_visc = .true.
2319 par_conv_visc = (par_conv_visc*tpl_grav/tpsf_nu_u_a)**num_1o3_sf
2320 q_sen_con = c_free_conv*sqrt(tpsf_kappa_t_a)*par_conv_visc
2321 q_sen_con = q_sen_con*(t_s-t_a)
2322 q_lat_con = c_free_conv*sqrt(tpsf_kappa_q_a)*par_conv_visc
2323 q_lat_con = q_lat_con*(q_s-q_a)
2325 l_conv_visc = .false.
2335r_z = height_tq/height_u
2336ri_cr = c_mo_t_stab/c_mo_u_stab**2_iintegers*r_z
2337ri = tpl_grav*((t_a-t_s)/t_s+tpsf_alpha_q*(q_a-q_s))/max(u_a,u_wind_min_sf)**2_iintegers
2338ri = ri*height_u/pr_neutral
2340turb_fluxes:
IF(u_a.LT.u_wind_min_sf.OR.ri.GT.ri_cr-c_small_sf)
THEN
2351 zol = sqrt(1.0-4.0*(c_mo_u_stab-r_z*c_mo_t_stab)*ri)
2352 zol = zol - 1.0 + 2.0*c_mo_u_stab*ri
2353 zol = zol/2.0/c_mo_u_stab/c_mo_u_stab/(ri_cr-ri)
2355 n_iter = 0_iintegers
2357 u_star_previter = ri*max(1.0, sqrt(r_z*c_mo_t_conv/c_mo_u_conv))
2358 DO WHILE (delta.GT.c_accur_sf.AND.n_iter.LT.n_iter_max)
2359 fun = u_star_previter**2_iintegers*(c_mo_u_conv*u_star_previter-1.0) &
2360 + ri**2_iintegers*(1.0-r_z*c_mo_t_conv*u_star_previter)
2361 fun_prime = 3.0*c_mo_u_conv*u_star_previter**2_iintegers &
2362 - 2.0*u_star_previter - r_z*c_mo_t_conv*ri**2_iintegers
2363 zol = u_star_previter - fun/fun_prime
2364 delta = abs(zol-u_star_previter)/max(c_accur_sf, abs(zol+u_star_previter))
2365 u_star_previter = zol
2366 n_iter = n_iter + 1_iintegers
2375CALL sfcflx_roughness (fetch, u_a, u_star_min_sf, h_ice, c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
2378u_star_st = u_star_thresh
2379CALL sfcflx_roughness (fetch, u_a, u_star_st, h_ice, c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
2381 psi_u = c_mo_u_stab*zol*(1.0-min(z0u_sf/height_u, 1.0))
2383 psi_t = (1.0-c_mo_u_conv*zol)**c_mo_u_exp
2384 psi_q = (1.0-c_mo_u_conv*zol*min(z0u_sf/height_u, 1.0))**c_mo_u_exp
2385 psi_u = 2.0*(atan(psi_t)-atan(psi_q)) &
2386 + 2.0*log((1.0+psi_q)/(1.0+psi_t)) &
2387 + log((1.0+psi_q*psi_q)/(1.0+psi_t*psi_t))
2389u_a_thresh = u_star_thresh/c_karman*(log(height_u/z0u_sf)+psi_u)
2394u_star_previter = u_star_thresh
2395IF(u_a.LE.u_a_thresh)
THEN
2396 DO WHILE (delta.GT.c_accur_sf.AND.n_iter.LT.n_iter_max)
2397 CALL sfcflx_roughness (fetch, u_a, min(u_star_thresh, u_star_previter), h_ice, &
2398 c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
2400 psi_u = c_mo_u_stab*zol*(1.0-min(z0u_sf/height_u, 1.0))
2401 fun = log(height_u/z0u_sf) + psi_u
2402 fun_prime = (fun + 1.0 + c_mo_u_stab*zol*min(z0u_sf/height_u, 1.0))/c_karman
2403 fun = fun*u_star_previter/c_karman - u_a
2405 psi_t = (1.0-c_mo_u_conv*zol)**c_mo_u_exp
2406 psi_q = (1.0-c_mo_u_conv*zol*min(z0u_sf/height_u, 1.0))**c_mo_u_exp
2407 psi_u = 2.0*(atan(psi_t)-atan(psi_q)) &
2408 + 2.0*log((1.0+psi_q)/(1.0+psi_t)) &
2409 + log((1.0+psi_q*psi_q)/(1.0+psi_t*psi_t))
2410 fun = log(height_u/z0u_sf) + psi_u
2411 fun_prime = (fun + 1.0/psi_q)/c_karman
2412 fun = fun*u_star_previter/c_karman - u_a
2414 u_star_st = u_star_previter - fun/fun_prime
2415 delta = abs((u_star_st-u_star_previter)/(u_star_st+u_star_previter))
2416 u_star_previter = u_star_st
2417 n_iter = n_iter + 1_iintegers
2420 DO WHILE (delta.GT.c_accur_sf.AND.n_iter.LT.n_iter_max)
2421 CALL sfcflx_roughness (fetch, u_a, max(u_star_thresh, u_star_previter), h_ice, &
2422 c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
2424 psi_u = c_mo_u_stab*zol*(1.0-min(z0u_sf/height_u, 1.0))
2425 fun = log(height_u/z0u_sf) + psi_u
2426 fun_prime = (fun - 2.0 - 2.0*c_mo_u_stab*zol*min(z0u_sf/height_u, 1.0))/c_karman
2427 fun = fun*u_star_previter/c_karman - u_a
2429 psi_t = (1.0-c_mo_u_conv*zol)**c_mo_u_exp
2430 psi_q = (1.0-c_mo_u_conv*zol*min(z0u_sf/height_u, 1.0))**c_mo_u_exp
2431 psi_u = 2.0*(atan(psi_t)-atan(psi_q)) &
2432 + 2.0*log((1.0+psi_q)/(1.0+psi_t)) &
2433 + log((1.0+psi_q*psi_q)/(1.0+psi_t*psi_t))
2434 fun = log(height_u/z0u_sf) + psi_u
2435 fun_prime = (fun - 2.0/psi_q)/c_karman
2436 fun = fun*u_star_previter/c_karman - u_a
2438 IF(h_ice.GE.h_ice_min_flk)
THEN
2439 u_star_st = c_karman*u_a/max(c_small_sf, log(height_u/z0u_sf)+psi_u)
2440 u_star_previter = u_star_st
2442 u_star_st = u_star_previter - fun/fun_prime
2444 delta = abs((u_star_st-u_star_previter)/(u_star_st+u_star_previter))
2445 u_star_previter = u_star_st
2446 n_iter = n_iter + 1_iintegers
2462q_mom_tur = -u_star_st*u_star_st
2465CALL sfcflx_roughness (fetch, u_a, u_star_st, h_ice, c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
2467 psi_t = c_mo_t_stab*r_z*zol*(1.0-min(z0t_sf/height_tq, 1.0))
2468 psi_q = c_mo_q_stab*r_z*zol*(1.0-min(z0q_sf/height_tq, 1.0))
2473 psi_u = (1.0-c_mo_t_conv*r_z*zol)**c_mo_t_exp
2474 psi_t = (1.0-c_mo_t_conv*r_z*zol*min(z0t_sf/height_tq, 1.0))**c_mo_t_exp
2476 psi_t = abs(2.0*log((1.0+psi_t)/(1.0+psi_u)))
2477 psi_u = (1.0-c_mo_q_conv*r_z*zol)**c_mo_q_exp
2478 psi_q = (1.0-c_mo_q_conv*r_z*zol*min(z0q_sf/height_tq, 1.0))**c_mo_q_exp
2480 psi_q = abs(2.0*log((1.0+psi_q)/(1.0+psi_u)))
2486q_sen_tur = -(t_a-t_s)*u_star_st*c_karman/pr_neutral &
2487 / max(c_small_sf, log(height_tq/z0t_sf)+psi_t)
2488if(max(c_small_sf, log(height_tq/z0t_sf)+psi_t) .lt. 10e-6)
then
2489 write(0,*)
'inside flake'
2490 write(0,*) q_sen_tur, t_a, t_s, u_star_st, c_karman, pr_neutral
2491 write(0,*) c_small_sf,height_tq,z0t_sf,psi_t
2492 write(0,*)
'nominator= ', (t_a-t_s)*u_star_st*c_karman/pr_neutral
2493 write(0,*)
'denominator= ',max(c_small_sf, log(height_tq/z0t_sf)+psi_t)
2495q_lat_tur = -(q_a-q_s)*u_star_st*c_karman/sc_neutral &
2496 / max(c_small_sf, log(height_tq/z0q_sf)+psi_q)
2497if(q_lat_tur .gt. 6.0e-4)
then
2498 q_lat_tur = -(q_a-q_s)*u_star_st*c_karman/3.0 &
2499 / max(c_small_sf, log(height_tq/z0q_sf)+psi_q)
2500 write(0,*)
'Q_lat_tur= ',q_lat_tur
2501 write(0,135) q_a,q_s,u_star_st,c_karman
2502 write(0,136) max(c_small_sf,log(height_tq/z0q_sf)+psi_q),c_small_sf, log(height_tq/z0q_sf),psi_q
2504135
format(1x,4(f16.4))
2505136
format(1x,4(f16.4))
2513q_momentum = min(q_mom_tur, q_mom_mol, q_mom_con)
2515 IF(abs(q_sen_tur).GE.abs(q_sen_con))
THEN
2516 q_sensible = q_sen_tur
2518 q_sensible = q_sen_con
2520 IF(abs(q_sensible).LT.abs(q_sen_mol))
THEN
2521 q_sensible = q_sen_mol
2523 IF(abs(q_lat_tur).GE.abs(q_lat_con))
THEN
2524 q_latent = q_lat_tur
2526 q_latent = q_lat_con
2528 IF(abs(q_latent).LT.abs(q_lat_mol))
THEN
2529 q_latent = q_lat_mol
2532 IF(abs(q_sen_tur).GE.abs(q_sen_mol))
THEN
2533 q_sensible = q_sen_tur
2535 q_sensible = q_sen_mol
2537 IF(abs(q_lat_tur).GE.abs(q_lat_mol))
THEN
2538 q_latent = q_lat_tur
2540 q_latent = q_lat_mol
2548q_momentum = q_momentum*rho_a
2552q_watvap = q_latent*rho_a
2555IF(h_ice.GE.h_ice_min_flk) q_latent = q_latent + tpl_l_f
2557q_latent = q_watvap*tpsf_l_evap
2558if(q_latent .gt. 2000.00)
then
2559 write(0,145)
'final Q_watvap= ',q_watvap,
'tpsf_L_evap= ',tpsf_l_evap,
'Q_latent= ', q_latent
2562145
format(a17,e12.5,1x,a13,1x,f10.2,1x,a10,1x,e12.4)
2564u_star_a_sf = u_star_st
2565q_mom_a_sf = q_momentum
2566q_sens_a_sf = q_sensible
2567q_lat_a_sf = q_latent
2568q_watvap_a_sf = q_watvap
2571 127
format(1x, 3(f16.5,1x))
2955SUBROUTINE flake_interface ( dMsnowdt_in, I_atm_in, Q_atm_lw_in, height_u_in, height_tq_in, &
2956 U_a_in, T_a_in, q_a_in, P_a_in, &
2958 depth_w, fetch, depth_bs, T_bs, par_Coriolis, del_time, &
2959 T_snow_in, T_ice_in, T_mnw_in, T_wML_in, T_bot_in, T_B1_in, &
2960 C_T_in, h_snow_in, h_ice_in, h_ML_in, H_B1_in, T_sfc_p, &
2961 ch, cm, albedo_water, water_extinc, &
2963 T_snow_out, T_ice_out, T_mnw_out, T_wML_out, T_bot_out, &
2964 T_B1_out, C_T_out, h_snow_out, h_ice_out, h_ML_out, &
2965 H_B1_out, T_sfc_n, hflx_out, evap_out, gflx_out, lflx_out, &
2967 T_bot_2_in, T_bot_2_out,ustar, q_sfc, chh, cmm )
3022 t_snow_p_flk, t_snow_n_flk , &
3023 t_ice_p_flk, t_ice_n_flk , &
3024 t_mnw_p_flk, t_mnw_n_flk , &
3025 t_wml_p_flk, t_wml_n_flk , &
3026 t_bot_p_flk, t_bot_n_flk , &
3027 t_b1_p_flk, t_b1_n_flk , &
3028 c_t_p_flk, c_t_n_flk , &
3029 h_snow_p_flk, h_snow_n_flk , &
3030 h_ice_p_flk, h_ice_n_flk , &
3031 h_ml_p_flk, h_ml_n_flk , &
3032 h_b1_p_flk, h_b1_n_flk , &
3055 sfcflx_lwradwsfc , &
3068real(kind = kind_phys),
INTENT(IN) :: &
3083real(kind = kind_phys),
INTENT(IN) :: &
3092real(kind = kind_phys),
INTENT(IN) :: &
3110real(kind = kind_phys) :: &
3115TYPE (opticpar_medium) :: &
3122real(kind = kind_phys),
INTENT(OUT) :: &
3147real(kind = kind_phys) :: &
3156LOGICAL lflk_botsed_use
3158real(kind = kind_phys),
parameter :: tpl_rho_w_r = 1.0e+03
3159real(kind = kind_phys),
parameter :: tpl_t_f = 273.15
3160real(kind = kind_phys),
parameter :: h_snow_min_flk = 1.0e-5
3161real(kind = kind_phys),
parameter :: h_ice_min_flk = 1.0e-9
3166 lflk_botsed_use = .false.
3179albedo_ice = albedo_whiteice_ref
3181albedo_snow = albedo_drysnow_ref
3182opticpar_water%extincoef_optic(1) = water_extinc
3190opticpar_water = opticpar_water_ref
3191opticpar_ice = opticpar_ice_opaque
3192opticpar_snow = opticpar_snow_opaque
3202t_snow_p_flk = t_snow_in
3203t_ice_p_flk = t_ice_in
3204t_mnw_p_flk = t_mnw_in
3205t_wml_p_flk = t_wml_in
3206t_bot_p_flk = t_bot_in
3209h_snow_p_flk = h_snow_in
3210h_ice_p_flk = h_ice_in
3213t_bot_2_in_flk = t_bot_2_in
3216 120
format(1x,6(f12.5,1x))
3221dmsnowdt_flk = dmsnowdt_in
3228CALL flake_radflux ( depth_w, albedo_water, albedo_ice, albedo_snow, &
3229 opticpar_water, opticpar_ice, opticpar_snow, &
3236q_w_flk = q_atm_lw_in
3237q_w_flk = q_w_flk - sfcflx_lwradwsfc(t_sfc_p)
3243CALL sfcflx_momsenlat ( height_u_in, height_tq_in, fetch, &
3244 u_a_in, t_a_in, q_a_in, t_sfc_p, p_a_in, h_ice_p_flk, &
3245 q_momentum, q_sensible, q_latent, q_watvap, q_sfc, rho_a )
3248u_star_w_flk = sqrt(-q_momentum/tpl_rho_w_r)
3255q_w_flk = q_w_flk - q_sensible - q_latent
3256IF(h_ice_p_flk.GE.h_ice_min_flk)
THEN
3257 IF(h_snow_p_flk.GE.h_snow_min_flk)
THEN
3258 q_snow_flk = q_w_flk
3275CALL flake_main ( depth_w, depth_bs, t_bs, par_coriolis, &
3276 opticpar_water%extincoef_optic(1), &
3277 del_time, t_sfc_p, t_sfc_n, t_bot_2_in_flk, &
3284t_snow_out = t_snow_n_flk
3285t_ice_out = t_ice_n_flk
3286t_mnw_out = t_mnw_n_flk
3287t_wml_out = t_wml_n_flk
3288t_bot_out = t_bot_n_flk
3289t_b1_out = t_b1_n_flk
3291h_snow_out = h_snow_n_flk
3292h_ice_out = h_ice_n_flk
3293h_ml_out = h_ml_n_flk
3294h_b1_out = h_b1_n_flk
3295hflx_out = q_sensible
3300chh = ch * u_a_in * rho_a