! $Id: calcul_fluxs_mod.F90 5144 2024-07-29 21:01:04Z abarral $ MODULE calcul_fluxs_mod IMPLICIT NONE CONTAINS SUBROUTINE calcul_fluxs(knon, nisurf, dtime, & tsurf, p1lay, cal, beta, cdragh, cdragq, ps, & precip_rain, precip_snow, snow, qsurf, & radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, gustiness, & fqsat, petAcoef, peqAcoef, petBcoef, peqBcoef, & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa) USE dimphy, ONLY: klon USE indice_sol_mod USE sens_heat_rain_m, ONLY: sens_heat_rain USE lmdz_clesphys USE lmdz_yoethf USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep USE lmdz_yomcst ! Cette routine calcule les fluxs en h et q a l'interface et eventuellement ! une temperature de surface (au cas ou ok_veget = false) ! L. Fairhead 4/2000 ! input: ! knon nombre de points a traiter ! nisurf surface a traiter ! tsurf temperature de surface ! p1lay pression 1er niveau (milieu de couche) ! cal capacite calorifique du sol ! beta evap reelle ! cdragh coefficient d'echange temperature ! cdragq coefficient d'echange evaporation ! ps pression au sol ! precip_rain precipitations liquides ! precip_snow precipitations solides ! snow champs hauteur de neige ! runoff runoff en cas de trop plein ! petAcoef coeff. A de la resolution de la CL pour t ! peqAcoef coeff. A de la resolution de la CL pour q ! petBcoef coeff. B de la resolution de la CL pour t ! peqBcoef coeff. B de la resolution de la CL pour q ! radsol rayonnement net aus sol (LW + SW) ! dif_grnd coeff. diffusion vers le sol profond ! output: ! tsurf_new temperature au sol ! qsurf humidite de l'air au dessus du sol ! fluxsens flux de chaleur sensible ! fluxlat flux de chaleur latente ! dflux_s derivee du flux de chaleur sensible / Ts ! dflux_l derivee du flux de chaleur latente / Ts ! sens_prec_liq flux sensible lié aux echanges de precipitations liquides ! sens_prec_sol precipitations solides ! lat_prec_liq flux latent lié aux echanges de precipitations liquides ! lat_prec_sol precipitations solides ! Parametres d'entree !**************************************************************************************** INTEGER, INTENT(IN) :: knon, nisurf REAL, INTENT(IN) :: dtime REAL, DIMENSION(klon), INTENT(IN) :: petAcoef, peqAcoef REAL, DIMENSION(klon), INTENT(IN) :: petBcoef, peqBcoef REAL, DIMENSION(klon), INTENT(IN) :: ps, q1lay REAL, DIMENSION(klon), INTENT(IN) :: tsurf, p1lay, cal, beta, cdragh, cdragq REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow REAL, DIMENSION(klon), INTENT(IN) :: radsol, dif_grnd REAL, DIMENSION(klon), INTENT(IN) :: t1lay, u1lay, v1lay, gustiness REAL, INTENT(IN) :: fqsat ! correction factor on qsat (generally 0.98 over salty water, 1 everywhere else) REAL, INTENT(IN), optional :: rhoa(:) ! (knon) ! density of moist air (kg / m3) ! Parametres entree-sorties !**************************************************************************************** REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! snow pas utile ! Parametres sorties !**************************************************************************************** REAL, DIMENSION(klon), INTENT(OUT) :: qsurf REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new, evap, fluxsens, fluxlat REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l REAL, INTENT(OUT), OPTIONAL :: sens_prec_liq(:), sens_prec_sol(:) ! (knon) REAL, DIMENSION(klon), OPTIONAL :: lat_prec_liq, lat_prec_sol ! Variables locales !**************************************************************************************** INTEGER :: i REAL, DIMENSION(klon) :: zx_mh, zx_nh, zx_oh REAL, DIMENSION(klon) :: zx_mq, zx_nq, zx_oq REAL, DIMENSION(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat REAL, DIMENSION(klon) :: zx_sl, zx_coefh, zx_coefq, zx_wind REAL, DIMENSION(klon) :: d_ts REAL :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh REAL :: qsat_new, q1_new REAL, PARAMETER :: t_grnd = 271.35, t_coup = 273.15 REAL, PARAMETER :: max_eau_sol = 150.0 CHARACTER (len = 20) :: modname = 'calcul_fluxs' LOGICAL :: fonte_neige LOGICAL, SAVE :: check = .FALSE. !$OMP THREADPRIVATE(check) ! End definition !**************************************************************************************** IF (check) WRITE(*, *)'Entree ', modname, ' surface = ', nisurf IF (check) THEN WRITE(*, *)' radsol (min, max)', & MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon)) ENDIF ! Traitement neige et humidite du sol !**************************************************************************************** !!$ WRITE(*,*)'test calcul_flux, surface ', nisurf !!PB test !!$ if (nisurf == is_oce) THEN !!$ snow = 0. !!$ qsol = max_eau_sol !!$ else !!$ where (precip_snow > 0.) snow = snow + (precip_snow * dtime) !!$ where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime)) !!$! snow = max(0.0, snow + (precip_snow - evap) * dtime) !!$ where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime !!$ endif !!$ IF (nisurf /= is_ter) qsol = max_eau_sol ! Initialisation !**************************************************************************************** evap = 0. fluxsens = 0. fluxlat = 0. dflux_s = 0. dflux_l = 0. IF (PRESENT(lat_prec_liq)) lat_prec_liq = 0. IF (PRESENT(lat_prec_sol)) lat_prec_sol = 0. ! zx_qs = qsat en kg/kg !**************************************************************************************** DO i = 1, knon zx_pkh(i) = (ps(i) / ps(i))**RKAPPA IF (thermcep) THEN zdelta = MAX(0., SIGN(1., rtt - tsurf(i))) zcvm5 = R5LES * RLVTT * (1. - zdelta) + R5IES * RLSTT * zdelta zcvm5 = zcvm5 / RCPD / (1.0 + RVTMP2 * q1lay(i)) zx_qs = r2es * FOEEW(tsurf(i), zdelta) / ps(i) zx_qs = MIN(0.5, zx_qs) zcor = 1. / (1. - retv * zx_qs) zx_qs = zx_qs * zcor zx_dq_s_dh = FOEDE(tsurf(i), zdelta, zcvm5, zx_qs, zcor) & / RLVTT / zx_pkh(i) ELSE IF (tsurf(i) epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) & ! & .AND. (tsurf_new(i) >= RTT) ! if (fonte_neige) tsurf_new(i) = RTT d_ts(i) = tsurf_new(i) - tsurf(i) ! zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i) ! zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i) !== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas !== flux_t est le flux de cpt (energie sensible): j/(m**2 s) evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i) fluxlat(i) = - evap(i) * zx_sl(i) fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i) ! Derives des flux dF/dTs (W m-2 K-1): dflux_s(i) = zx_nh(i) dflux_l(i) = (zx_sl(i) * zx_nq(i)) ! Nouvelle valeure de l'humidite au dessus du sol qsat_new = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i) q1_new = peqAcoef(i) - peqBcoef(i) * evap(i) * dtime qsurf(i) = q1_new * (1. - beta(i)) + beta(i) * qsat_new ! en cas de fonte de neige ! if (fonte_neige) THEN ! bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - & ! & dif_grnd(i) * (tsurf_new(i) - t_grnd) - & ! & RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i)) ! bilan_f = max(0., bilan_f) ! fq_fonte = bilan_f / zx_sl(i) ! snow(i) = max(0., snow(i) - fq_fonte * dtime) ! qsol(i) = qsol(i) + (fq_fonte * dtime) ! endif !!$ if (nisurf == is_ter) & !!$ & run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0) !!$ qsol(i) = min(qsol(i), max_eau_sol) ! calcul de l'enthalpie des precipitations liquides et solides IF (PRESENT(sens_prec_liq)) sens_prec_liq(i) & = - sens_heat_rain(precip_rain(i) + precip_snow(i), t1lay(i), & q1lay(i), rhoa(i), rlvtt, tsurf_new(i), ps(i)) IF (PRESENT(sens_prec_sol)) sens_prec_sol(i) = 0. ! On calcule par rapport a T=0 !! sens_prec_liq(i) = rcw * (t1lay(i) - RTT) * precip_rain(i) !! sens_prec_sol(i) = rcs * (t1lay(i) - RTT) * precip_snow(i) IF (PRESENT(lat_prec_liq)) & lat_prec_liq(i) = precip_rain(i) * (RLVTT - RLVTT) IF (PRESENT(lat_prec_sol)) & lat_prec_sol(i) = precip_snow(i) * (RLSTT - RLVTT) ENDDO !************************************************************************** END SUBROUTINE calcul_fluxs !**************************************************************************************** SUBROUTINE calcul_flux_wind(knon, dtime, & u0, v0, u1, v1, gustiness, cdrag_m, & AcoefU, AcoefV, BcoefU, BcoefV, & p1lay, t1lay, & flux_u1, flux_v1) USE dimphy USE lmdz_clesphys USE lmdz_yomcst IMPLICIT NONE ! Input arguments !**************************************************************************************** INTEGER, INTENT(IN) :: knon REAL, INTENT(IN) :: dtime REAL, DIMENSION(klon), INTENT(IN) :: u0, v0 ! u and v at niveau 0 REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness ! u and v at niveau 1 REAL, DIMENSION(klon), INTENT(IN) :: cdrag_m ! cdrag pour momentum REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV REAL, DIMENSION(klon), INTENT(IN) :: p1lay ! pression 1er niveau (milieu de couche) REAL, DIMENSION(klon), INTENT(IN) :: t1lay ! temperature ! Output arguments !**************************************************************************************** REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1 REAL, DIMENSION(klon), INTENT(OUT) :: flux_v1 ! Local variables !**************************************************************************************** INTEGER :: i REAL :: mod_wind, buf !**************************************************************************************** ! Calculate the surface flux !**************************************************************************************** DO i = 1, knon mod_wind = min_wind_speed + SQRT(gustiness(i) + (u1(i) - u0(i))**2 + (v1(i) - v0(i))**2) buf = cdrag_m(i) * mod_wind * p1lay(i) / (RD * t1lay(i)) flux_u1(i) = (AcoefU(i) - u0(i)) / (1 / buf - BcoefU(i) * dtime) flux_v1(i) = (AcoefV(i) - v0(i)) / (1 / buf - BcoefV(i) * dtime) END DO END SUBROUTINE calcul_flux_wind !**************************************************************************************** END MODULE calcul_fluxs_mod