Changeset 181 for LMDZ.3.3/branches/rel-LF/libf/phylmd

Timestamp:

Mar 27, 2001, 12:31:52 PM (24 years ago)

Author:

lmdzadmin

Message:

Traitement de la fonte de la neige dans le cas du modele de sol simplifie
est sorti de calcul_flux
LF

File:

• LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90 (modified) (8 diffs)

LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90

@@ -329 +329 @@
      &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
 
+       CALL fonte_neige( klon, knon, nisurf, dtime, &
+     &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &
+     &   precip_rain, precip_snow, snow, qsol,  &
+     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
+     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
+     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+
     else
         CALL albsno(klon,agesno,alb_neig_grid)  
@@ -476 +483 @@
      &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
 
+    if (ocean /= 'couple') then
+      call fonte_neige( klon, knon, nisurf, dtime, &
+     &   tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
+     &   precip_rain, precip_snow, snow, qsol,  &
+     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
+     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
+     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+    endif
 !
 ! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean
@@ -526 +541 @@
 
     call calcul_fluxs( klon, knon, nisurf, dtime, &
+     &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &
+     &   precip_rain, precip_snow, snow, qsol,  &
+     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
+     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
+     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+
+    call fonte_neige( klon, knon, nisurf, dtime, &
      &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &
      &   precip_rain, precip_snow, snow, qsol,  &
@@ -1921 +1943 @@
   real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
   real, dimension(klon) :: zx_sl, zx_k1
-  real, dimension(klon) :: zx_h_ts, zx_q_0 , d_ts
+  real, dimension(klon) :: zx_q_0 , d_ts
   real                  :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
   real                  :: bilan_f, fq_fonte
@@ -2023 +2045 @@
     zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
 
-
 ! Tsurface
     tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
@@ -2035 +2056 @@
 ! Y'a-t-il fonte de neige?
 !
-    fonte_neige = (nisurf /= is_oce) .AND. &
-     & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
-     & .AND. (tsurf_new(i) >= RTT)
-    if (fonte_neige) tsurf_new(i) = RTT  
-    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
+!    fonte_neige = (nisurf /= is_oce) .AND. &
+!     & (snow(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_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(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)
@@ -2054 +2075 @@
 ! 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 (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)
@@ -2196 +2217 @@
 !#########################################################################
 !
+
+  SUBROUTINE fonte_neige( klon, knon, nisurf, dtime, &
+     & tsurf, p1lay, cal, beta, coef1lay, ps, &
+     & precip_rain, precip_snow, snow, qsol, &
+     & radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
+     & petAcoef, peqAcoef, petBcoef, peqBcoef, &
+     & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+
+! Routine de traitement de la fonte de la neige dans le cas du traitement
+! de sol simplifié
+!
+! LF 03/2001
+! 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
+!   coef1lay     coefficient d'echange
+!   ps           pression au sol
+!   precip_rain  precipitations liquides
+!   precip_snow  precipitations solides
+!   snow         champs hauteur de neige
+!   qsol         humidite du sol
+!   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
+!   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
+!
+
+#include "YOETHF.inc"
+#include "FCTTRE.inc"
+
+! Parametres d'entree
+  integer, intent(IN) :: knon, nisurf, klon
+  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, coef1lay
+  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
+  real, dimension(klon), intent(INOUT) :: snow, qsol
+
+! Parametres sorties
+  real, dimension(klon), intent(INOUT):: tsurf_new, evap, fluxsens, fluxlat
+  real, dimension(klon), intent(INOUT):: dflux_s, dflux_l
+
+! 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, zx_coef
+  real, dimension(klon) :: zx_sl, zx_k1
+  real, dimension(klon) :: zx_q_0 , d_ts
+  real                  :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
+  real                  :: bilan_f, fq_fonte
+  REAL                  :: subli, fsno
+  real, parameter :: t_grnd = 271.35, t_coup = 273.15
+!! PB temporaire en attendant mieux pour le modele de neige
+  REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)
+!
+  logical, save         :: check = .true.
+  character (len = 20)  :: modname = 'fonte_neige'
+  logical, save         :: neige_fond = .false.
+  real, save            :: max_eau_sol = 150.0
+  character (len = 80) :: abort_message 
+  logical,save         :: first = .true.,second=.false.
+
+  if (check) write(*,*)'Entree ', modname,' surface = ',nisurf
+
+! Initialisations
+  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).LT.t_coup) THEN
+        zx_qs = qsats(tsurf(i)) / ps(i)
+        zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
+     &                    / zx_pkh(i)
+      ELSE
+        zx_qs = qsatl(tsurf(i)) / ps(i)
+        zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
+     &               / zx_pkh(i)
+      ENDIF
+    ENDIF
+    zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
+    zx_qsat(i) = zx_qs
+    zx_coef(i) = coef1lay(i) &
+     & * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &
+     & * p1lay(i)/(RD*t1lay(i))
+  ENDDO
+! === Calcul de la temperature de surface ===
+! 
+! zx_sl = chaleur latente d'evaporation ou de sublimation
+!
+  do i = 1, knon
+    zx_sl(i) = RLVTT
+    if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
+    zx_k1(i) = zx_coef(i)
+  enddo
+
+
+  do i = 1, knon
+! Q
+    zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
+    zx_mq(i) = beta(i) * zx_k1(i) * &
+     &             (peqAcoef(i) - zx_qsat(i) &
+     &                          + zx_dq_s_dt(i) * tsurf(i)) &
+     &             / zx_oq(i)
+    zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
+     &                              / zx_oq(i)
+
+! H
+    zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
+    zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
+    zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
+  enddo
+
+!
+! Y'a-t-il fonte de neige?
+!
+  do i = 1, knon
+    neige_fond = ((snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
+     & .AND. tsurf_new(i) >= RTT)
+    if (neige_fond) then
+      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))
+      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)
+      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) 
+    endif
+  enddo
+
+  END SUBROUTINE fonte_neige
+!
+!#########################################################################
+!
   END MODULE interface_surf