83 USE yomhook
,ONLY : lhook, dr_hook
84 USE parkind1
,ONLY : jprb
100 c_MO_q_stab = 5.0 , &
101 c_MO_u_conv = 15.0 , &
102 c_MO_t_conv = 15.0 , &
103 c_MO_q_conv = 15.0 , &
104 c_MO_u_exp = 0.25 , &
107 z0u_ice_rough = 1.0E-03 , &
108 c_z0u_smooth = 0.1 , &
109 c_z0u_rough = 1.23E-02 , &
110 c_z0u_rough_L = 1.00E-01 , &
111 c_z0u_ftch_f = 0.70 , &
112 c_z0u_ftch_ex = 0.3333333 , &
113 c_z0t_rough_1 = 4.0 , &
114 c_z0t_rough_2 = 3.2 , &
115 c_z0t_rough_3 = 0.5 , &
116 c_z0q_rough_1 = 4.0 , &
117 c_z0q_rough_2 = 4.2 , &
118 c_z0q_rough_3 = 0.5 , &
119 c_z0t_ice_b0s = 1.250 , &
120 c_z0t_ice_b0t = 0.149 , &
121 c_z0t_ice_b1t = -0.550 , &
122 c_z0t_ice_b0r = 0.317 , &
123 c_z0t_ice_b1r = -0.565 , &
124 c_z0t_ice_b2r = -0.183 , &
125 c_z0q_ice_b0s = 1.610 , &
126 c_z0q_ice_b0t = 0.351 , &
127 c_z0q_ice_b1t = -0.628 , &
128 c_z0q_ice_b0r = 0.396 , &
129 c_z0q_ice_b1r = -0.512 , &
130 c_z0q_ice_b2r = -0.180 , &
131 Re_z0s_ice_t = 2.5 , &
138 REAL ,
PARAMETER :: &
142 REAL ,
PARAMETER :: &
146 REAL ,
PARAMETER :: &
147 tpsf_C_StefBoltz = 5.67E-08 , &
148 tpsf_R_dryair = 2.8705E+02 , &
149 tpsf_R_watvap = 4.6151E+02 , &
150 tpsf_c_a_p = 1.005E+03 , &
151 tpsf_L_evap = 2.501E+06 , &
152 tpsf_nu_u_a = 1.50E-05 , &
153 tpsf_kappa_t_a = 2.20E-05 , &
154 tpsf_kappa_q_a = 2.40E-05
157 REAL ,
PARAMETER :: &
158 tpsf_Rd_o_Rv = tpsf_R_dryair/tpsf_R_watvap , &
159 tpsf_alpha_q = (1.-tpsf_Rd_o_Rv)/tpsf_Rd_o_Rv
162 REAL ,
PARAMETER :: &
177 REAL,
PARAMETER,
PRIVATE :: Z_=-HUGE(0.0)
186 REAL ,
PARAMETER :: &
187 u_wind_min_sf = 1.0E-02 , &
188 u_star_min_sf = 1.0E-04 , &
189 z0t_min_sf = 1.0E-11 , &
190 c_accur_sf = 1.0E-07 , &
194 REAL ,
PARAMETER :: &
264 REAL ,
INTENT(IN) :: &
279 REAL ,
PARAMETER :: &
280 c_lmmgo_1 = 43.057924 , &
283 INTEGER ,
PARAMETER :: &
285 REAL ,
PARAMETER,
DIMENSION (nband_coef) :: &
286 corr_cl_tot = (/0.70, 0.45, 0.32, &
287 0.23, 0.18, 0.13/) , &
288 corr_cl_low = (/0.76, 0.49, 0.35, &
289 0.26, 0.20, 0.15/) , &
290 corr_cl_midhigh = (/0.46, 0.30, 0.21, &
292 REAL ,
PARAMETER :: &
298 REAL ,
PARAMETER :: &
299 c_watvap_corr_min = 0.6100 , &
300 c_watvap_corr_max = 0.7320 , &
301 c_watvap_corr_e = 0.0050
311 c_cl_midhigh_corr , &
315 REAL(KIND=JPRB) :: zhook_handle
322 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_LWRADATM',0,zhook_handle)
323 f_wvpres_corr = c_watvap_corr_min + c_watvap_corr_e*sqrt(e_a)
324 f_wvpres_corr = min(f_wvpres_corr, c_watvap_corr_max)
327 IF(t_a.LT.t_low)
THEN
328 c_cl_tot_corr = corr_cl_tot(1)
329 c_cl_low_corr = corr_cl_low(1)
330 c_cl_midhigh_corr = corr_cl_midhigh(1)
331 ELSE IF(t_a.GE.t_low+(nband_coef-1)*del_t)
THEN
332 c_cl_tot_corr = corr_cl_tot(nband_coef)
333 c_cl_low_corr = corr_cl_low(nband_coef)
334 c_cl_midhigh_corr = corr_cl_midhigh(nband_coef)
338 IF(t_a.GE.t_corr.AND.t_a.LT.t_corr+del_t)
THEN
339 c_cl_tot_corr = (t_a-t_corr)/del_t
340 c_cl_low_corr = corr_cl_low(i) + (corr_cl_low(i+1)-corr_cl_low(i))*c_cl_tot_corr
341 c_cl_midhigh_corr = corr_cl_midhigh(i) + (corr_cl_midhigh(i+1)-corr_cl_midhigh(i))*c_cl_tot_corr
342 c_cl_tot_corr = corr_cl_tot(i) + (corr_cl_tot(i+1)-corr_cl_tot(i))*c_cl_tot_corr
344 t_corr = t_corr + del_t
348 IF(cl_low.LT.0.)
THEN
349 f_cloud_corr = 1. + c_cl_tot_corr*cl_tot*cl_tot
351 f_cloud_corr = (1. + c_cl_low_corr*cl_low*cl_low) &
352 * (1. + c_cl_midhigh_corr*(cl_tot*cl_tot-cl_low*cl_low))
366 * f_wvpres_corr*f_cloud_corr
367 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_LWRADATM',1,zhook_handle)
430 REAL ,
INTENT(IN) :: &
437 REAL(KIND=JPRB) :: zhook_handle
445 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_LWRADWSFC',0,zhook_handle)
447 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_LWRADWSFC',1,zhook_handle)
463 u_a, t_a, q_a, t_s, p_a, h_ice, &
464 q_momentum, q_sensible, q_latent, q_watvap,&
465 ri, z0u_ini, z0t_ini, qsat_out, &
466 q_latenti, q_sublim )
526 REAL ,
INTENT(IN) :: &
544 REAL ,
INTENT(INOUT) :: &
546 REAL ,
INTENT(OUT) :: &
557 INTEGER ,
PARAMETER :: &
609 REAL(KIND=JPRB) :: zhook_handle
622 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_MOMSENLAT',0,zhook_handle)
630 q_mom_mol = -tpsf_nu_u_a*u_a/height_u
631 q_sen_mol = -tpsf_kappa_t_a*(t_a-t_s)/height_tq
632 q_lat_mol = -tpsf_kappa_q_a*(q_a-q_s)/height_tq
638 par_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)
639 IF(par_conv_visc.GT.0.)
THEN
641 par_conv_visc = (par_conv_visc*tpl_grav/tpsf_nu_u_a)**num_1o3_sf
642 q_sen_con = c_free_conv*sqrt(tpsf_kappa_t_a)*par_conv_visc
643 q_sen_con = q_sen_con*(t_s-t_a)
644 q_lat_con = c_free_conv*sqrt(tpsf_kappa_q_a)*par_conv_visc
645 q_lat_con = q_lat_con*(q_s-q_a)
647 l_conv_visc = .false.
657 r_z = height_tq/height_u
658 ri_cr = c_mo_t_stab/c_mo_u_stab**2*r_z
659 ri = tpl_grav*((t_a-t_s)/t_s+tpsf_alpha_q*(q_a-q_s))/max(u_a,u_wind_min_sf)**2
660 ri = min(xrimax,ri*height_u/pr_neutral)
662 turb_fluxes:
IF(u_a.LT.u_wind_min_sf.OR.ri.GT.ri_cr-c_small_sf)
THEN
675 zol = sqrt(1.-4.*(c_mo_u_stab-r_z*c_mo_t_stab)*ri)
676 zol = zol - 1. + 2.*c_mo_u_stab*ri
677 zol = zol/2./c_mo_u_stab/c_mo_u_stab/(ri_cr-ri)
681 u_star_previter = ri*max(1., sqrt(r_z*c_mo_t_conv/c_mo_u_conv))
682 DO WHILE (delta.GT.c_accur_sf.AND.n_iter.LT.n_iter_max)
683 fun = u_star_previter**2*(c_mo_u_conv*u_star_previter-1.) &
684 + ri**2*(1.-r_z*c_mo_t_conv*u_star_previter)
685 fun_prime = 3.*c_mo_u_conv*u_star_previter**2 &
686 - 2.*u_star_previter - r_z*c_mo_t_conv*ri**2
687 zol = u_star_previter - fun/fun_prime
688 delta = abs(zol-u_star_previter)/max(c_accur_sf, abs(zol+u_star_previter))
689 u_star_previter = zol
699 CALL
sfcflx_roughness(fetch, u_a, u_star_min_sf, h_ice, c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
702 u_star_st = u_star_thresh
703 CALL
sfcflx_roughness(fetch, u_a, u_star_st, h_ice, c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
705 psi_u = c_mo_u_stab*zol*(1.-min(z0u_sf/height_u, 1.))
707 psi_t = (1.-c_mo_u_conv*zol)**c_mo_u_exp
708 psi_q = (1.-c_mo_u_conv*zol*min(z0u_sf/height_u, 1.))**c_mo_u_exp
709 psi_u = 2.*(atan(psi_t)-atan(psi_q)) &
710 + 2.*log((1.+psi_q)/(1.+psi_t)) &
711 + log((1.+psi_q*psi_q)/(1.+psi_t*psi_t))
713 u_a_thresh = u_star_thresh/c_karman*(log(height_u/z0u_sf)+psi_u)
720 u_star_previter = max( sqrt( - q_momentum / rho_a ) , u_star_thresh )
723 IF(u_a.LE.u_a_thresh)
THEN
724 DO WHILE (delta.GT.c_accur_sf.AND.n_iter.LT.n_iter_max)
725 CALL
sfcflx_roughness(fetch, u_a, min(u_star_thresh, u_star_previter), h_ice, &
726 c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
728 psi_u = c_mo_u_stab*zol*(1.-min(z0u_sf/height_u, 1.))
729 fun = log(height_u/z0u_sf) + psi_u
730 fun_prime = (fun + 1. + c_mo_u_stab*zol*min(z0u_sf/height_u, 1.))/c_karman
731 fun = fun*u_star_previter/c_karman - u_a
733 psi_t = (1.-c_mo_u_conv*zol)**c_mo_u_exp
734 psi_q = (1.-c_mo_u_conv*zol*min(z0u_sf/height_u, 1.))**c_mo_u_exp
735 psi_u = 2.*(atan(psi_t)-atan(psi_q)) &
736 + 2.*log((1.+psi_q)/(1.+psi_t)) &
737 + log((1.+psi_q*psi_q)/(1.+psi_t*psi_t))
738 fun = log(height_u/z0u_sf) + psi_u
739 fun_prime = (fun + 1./psi_q)/c_karman
740 fun = fun*u_star_previter/c_karman - u_a
742 u_star_st = u_star_previter - fun/fun_prime
743 delta = abs((u_star_st-u_star_previter)/(u_star_st+u_star_previter))
744 u_star_previter = u_star_st
748 DO WHILE (delta.GT.c_accur_sf.AND.n_iter.LT.n_iter_max.AND.z0t_sf>z0t_min_sf)
749 CALL
sfcflx_roughness(fetch, u_a, max(u_star_thresh, u_star_previter), h_ice, &
750 c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
752 psi_u = c_mo_u_stab*zol*(1.-min(z0u_sf/height_u, 1.))
753 fun = log(height_u/z0u_sf) + psi_u
754 fun_prime = (fun - 2. - 2.*c_mo_u_stab*zol*min(z0u_sf/height_u, 1.))/c_karman
755 fun = fun*u_star_previter/c_karman - u_a
757 psi_t = (1.-c_mo_u_conv*zol)**c_mo_u_exp
758 psi_q = (1.-c_mo_u_conv*zol*min(z0u_sf/height_u, 1.))**c_mo_u_exp
759 psi_u = 2.*(atan(psi_t)-atan(psi_q)) &
760 + 2.*log((1.+psi_q)/(1.+psi_t)) &
761 + log((1.+psi_q*psi_q)/(1.+psi_t*psi_t))
762 fun = log(height_u/z0u_sf) + psi_u
763 fun_prime = (fun - 2./psi_q)/c_karman
764 fun = fun*u_star_previter/c_karman - u_a
766 IF(h_ice.GE.h_ice_min_flk)
THEN
767 u_star_st = c_karman*u_a/max(c_small_sf, log(height_u/z0u_sf)+psi_u)
768 u_star_previter = u_star_st
770 u_star_st = u_star_previter - fun/fun_prime
772 delta = abs((u_star_st-u_star_previter)/(u_star_st+u_star_previter))
773 u_star_previter = u_star_st
783 IF (n_iter == n_iter_max .OR. z0t_sf<=z0t_min_sf)
THEN
784 u_star_st = max( sqrt( - q_momentum / rho_a ) , u_star_min_sf )
786 c_z0u_fetch, u_star_previter, z0u_sf, z0t_sf, z0q_sf)
788 u_star_st = c_karman*u_a/max(c_small_sf, log(height_u/z0u_sf))
803 q_mom_tur = -u_star_st*u_star_st
806 CALL
sfcflx_roughness(fetch, u_a, u_star_st, h_ice, c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
809 psi_t = c_mo_t_stab*r_z*zol*(1.-min(z0t_sf/height_tq, 1.))
810 psi_q = c_mo_q_stab*r_z*zol*(1.-min(z0q_sf/height_tq, 1.))
815 psi_u = (1.-c_mo_t_conv*r_z*zol)**c_mo_t_exp
816 psi_t = (1.-c_mo_t_conv*r_z*zol*min(z0t_sf/height_tq, 1.))**c_mo_t_exp
817 psi_t = 2.*log((1.+psi_t)/(1.+psi_u))
818 psi_u = (1.-c_mo_q_conv*r_z*zol)**c_mo_q_exp
819 psi_q = (1.-c_mo_q_conv*r_z*zol*min(z0q_sf/height_tq, 1.))**c_mo_q_exp
820 psi_q = 2.*log((1.+psi_q)/(1.+psi_u))
826 q_sen_tur = -(t_a-t_s)*u_star_st*c_karman/pr_neutral &
827 / max(c_small_sf, log(height_tq/z0t_sf)+psi_t)
828 q_lat_tur = -(q_a-q_s)*u_star_st*c_karman/sc_neutral &
829 / max(c_small_sf, log(height_tq/z0q_sf)+psi_q)
837 q_momentum = min(q_mom_tur, q_mom_mol, q_mom_con)
839 IF(abs(q_sen_tur).GE.abs(q_sen_con))
THEN
840 q_sensible = q_sen_tur
842 q_sensible = q_sen_con
844 IF(abs(q_sensible).LT.abs(q_sen_mol))
THEN
845 q_sensible = q_sen_mol
847 IF(abs(q_lat_tur).GE.abs(q_lat_con))
THEN
852 IF(abs(q_latent).LT.abs(q_lat_mol))
THEN
856 IF(abs(q_sen_tur).GE.abs(q_sen_mol))
THEN
857 q_sensible = q_sen_tur
859 q_sensible = q_sen_mol
861 IF(abs(q_lat_tur).GE.abs(q_lat_mol))
THEN
872 q_momentum = q_momentum*rho_a
873 q_sensible = q_sensible*rho_a*tpsf_c_a_p
874 q_watvap = q_latent*rho_a
876 IF(h_ice.GE.h_ice_min_flk)
THEN
877 q_latent = q_watvap * (tpsf_l_evap + tpl_l_f)
878 q_latenti = q_watvap * (tpsf_l_evap + tpl_l_f)
881 q_latent = q_watvap * tpsf_l_evap
888 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_MOMSENLAT',1,zhook_handle)
953 REAL ,
INTENT(IN) :: &
961 REAL(KIND=JPRB) :: zhook_handle
969 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_RHOAIR',0,zhook_handle)
970 sfcflx_rhoair = p/tpsf_r_dryair/t/(1.+(1./tpsf_rd_o_rv-1.)*q)
971 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_RHOAIR',1,zhook_handle)
987 c_z0u_fetch, u_star_thresh, z0u, z0t, z0q)
1047 REAL ,
INTENT(IN) :: &
1054 REAL ,
INTENT(OUT) :: &
1065 REAL(KIND=JPRB) :: zhook_handle
1071 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_ROUGHNESS',0,zhook_handle)
1072 water_or_ice:
IF(h_ice.LT.h_ice_min_flk)
THEN
1075 c_z0u_fetch = max(u_a, u_wind_min_sf)**2/tpl_grav/fetch
1076 c_z0u_fetch = c_z0u_rough + c_z0u_ftch_f*c_z0u_fetch**c_z0u_ftch_ex
1077 c_z0u_fetch = min(c_z0u_fetch, c_z0u_rough_l)
1080 u_star_thresh = (c_z0u_smooth/c_z0u_fetch*tpl_grav*tpsf_nu_u_a)**num_1o3_sf
1083 re_s = u_star**3/tpsf_nu_u_a/tpl_grav
1084 re_s_thresh = c_z0u_smooth/c_z0u_fetch
1087 IF(re_s.LE.re_s_thresh)
THEN
1088 z0u = c_z0u_smooth*tpsf_nu_u_a/u_star
1090 z0u = c_z0u_fetch*u_star*u_star/tpl_grav
1093 z0q = c_z0u_fetch*max(re_s, re_s_thresh)
1094 z0t = c_z0t_rough_1*z0q**c_z0t_rough_3 - c_z0t_rough_2
1095 z0q = c_z0q_rough_1*z0q**c_z0q_rough_3 - c_z0q_rough_2
1096 z0t = z0u*exp(-c_karman/pr_neutral*z0t)
1097 z0q = z0u*exp(-c_karman/sc_neutral*z0q)
1102 c_z0u_fetch = c_z0u_rough
1105 u_star_thresh = c_z0u_smooth*tpsf_nu_u_a/z0u_ice_rough
1108 z0u = max(z0u_ice_rough, c_z0u_smooth*tpsf_nu_u_a/u_star)
1111 re_s = max(u_star*z0u/tpsf_nu_u_a, c_accur_sf)
1114 IF(re_s.LE.re_z0s_ice_t)
THEN
1115 z0t = c_z0t_ice_b0t + c_z0t_ice_b1t*log(re_s)
1116 z0t = min(z0t, c_z0t_ice_b0s)
1117 z0q = c_z0q_ice_b0t + c_z0q_ice_b1t*log(re_s)
1118 z0q = min(z0q, c_z0q_ice_b0s)
1120 z0t = c_z0t_ice_b0r + c_z0t_ice_b1r*log(re_s) + c_z0t_ice_b2r*log(re_s)**2
1121 z0q = c_z0q_ice_b0r + c_z0q_ice_b1r*log(re_s) + c_z0q_ice_b2r*log(re_s)**2
1127 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_ROUGHNESS',1,zhook_handle)
1194 REAL ,
INTENT(IN) :: &
1203 REAL ,
PARAMETER :: &
1206 b2w_vap = 17.2693882 , &
1207 b2i_vap = 21.8745584 , &
1210 REAL(KIND=JPRB) :: zhook_handle
1218 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_SATWVPRES',0,zhook_handle)
1219 IF(h_ice.LT.h_ice_min_flk)
THEN
1224 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_SATWVPRES',1,zhook_handle)
1287 REAL ,
INTENT(IN) :: &
1294 REAL(KIND=JPRB) :: zhook_handle
1302 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_SPECHUM',0,zhook_handle)
1303 sfcflx_spechum = tpsf_rd_o_rv*wvpres/(p-(1.-tpsf_rd_o_rv)*wvpres)
1304 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_SPECHUM',1,zhook_handle)
1368 REAL ,
INTENT(IN) :: &
1377 REAL(KIND=JPRB) :: zhook_handle
1385 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_WVPRESWETBULB',0,zhook_handle)
1387 - tpsf_c_a_p*p/tpsf_l_evap/tpsf_rd_o_rv*(t_dry-t_wetbulb)
1388 IF (lhook) CALL dr_hook(
'SFCFLX:SFCFLX_WVPRESWETBULB',1,zhook_handle)
real function sfcflx_lwradwsfc(emis, pts)
subroutine sfcflx_momsenlat(height_u, height_tq, fetch, U_a, T_a, q_a, T_s, P_a, h_ice, Q_momentum, Q_sensible, Q_latent, Q_watvap, Ri, z0u_ini, z0t_ini, Qsat_out, Q_latenti, Q_sublim)
real function sfcflx_spechum(wvpres, P)
real function sfcflx_lwradatm(T_a, e_a, cl_tot, cl_low)
subroutine sfcflx_roughness(fetch, U_a, u_star, h_ice, c_z0u_fetch, u_star_thresh, z0u, z0t, z0q)
real function sfcflx_rhoair(T, q, P)
real function sfcflx_wvpreswetbulb(T_dry, T_wetbulb, satwvpres_bulb, P)
real function sfcflx_satwvpres(T, h_ice)