Changeset 2500 for LMDZ5/trunk

Timestamp:

Apr 28, 2016, 4:58:11 PM (9 years ago)

Author:

jbmadeleine

Message:

Ajouts de commentaires dans fisrtilp pour clarifier la lecture.
Ajout de flags P1> P2> P3> pour se situer dans le code plus facilement.
Aucune incidence sur les resultats du modele.

File:

• LMDZ5/trunk/libf/phylmd/fisrtilp.F90 (modified) (36 diffs)

LMDZ5/trunk/libf/phylmd/fisrtilp.F90

@@ -21 +21 @@
   ! Objet: condensation et precipitation stratiforme.
   !        schema de nuage
+  ! Fusion de fisrt (physique sursaturation, P. LeVan K. Laval)
+  !             et ilp (il pleut, L. Li)
+  ! Principales parties:
+  ! P1> Evaporation de la precipitation (qui vient du niveau k+1)
+  ! P2> Formation du nuage (en k) 
+  ! P3> Formation de la precipitation (en k)
   !======================================================================
   !======================================================================
@@ -28 +34 @@
 
   !
-  ! Arguments:
+  ! Principaux inputs:
   !
   REAL dtime ! intervalle du temps (s)
@@ -35 +41 @@
   REAL t(klon,klev) ! temperature (K)
   REAL q(klon,klev) ! humidite specifique (kg/kg)
+  !
+  ! Principaux outputs:
+  !
   REAL d_t(klon,klev) ! incrementation de la temperature (K)
   REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
@@ -46 +55 @@
   REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
   REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s) 
+  !
+  ! Autres arguments
+  ! 
   REAL ztv(klon,klev)
   REAL zqta(klon,klev),fraca(klon,klev) 
@@ -155 +167 @@
   !     Pour la conversion eau-neige
   REAL zlh_solid(klon), zm_solid
-  !IM 
-  !ym      INTEGER klevm1
   !---------------------------------------------------------------
   !
   ! Fonctions en ligne:
   !
-  REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
+  REAL fallvs,fallvc ! Vitesse de chute pour cristaux de glace
+                     ! (Heymsfield & Donner, 1990)
   REAL zzz
   include "YOETHF.h"
@@ -278 +289 @@
   zfrac_lessi = 0.
 
-  !AA----------------------------------------------------------
+  !AA==================================================================
   !
   ncoreczq=0
-  ! Boucle verticale (du haut vers le bas)
-  !
-  !IM : klevm1
-  !ym      klevm1=klev-1
+  ! BOUCLE VERTICALE (DU HAUT VERS LE BAS)
+  !
   DO k = klev, 1, -1
      !
-     !AA----------------------------------------------------------
-     !
+     !AA===============================================================
+     !
+     ! Initialisation temperature et vapeur
      DO i = 1, klon
         zt(i)=t(i,k)
@@ -312 +322 @@
         ENDDO
      ENDIF
-     !
-     !
-     ! Calculer l'evaporation de la precipitation
-     !
-
-
-     ! Calculer l'evaporation de la precipitation
-     !
-
-
+     ! ----------------------------------------------------------------
+     ! P1> Debut evaporation de la precipitation
+     ! ----------------------------------------------------------------
      IF (evap_prec) THEN
         DO i = 1, klon
@@ -328 +331 @@
            IF (zrfl(i)+zifl(i).GT.0.) THEN
 !>AJ
+              ! Calcul du qsat
               IF (thermcep) THEN
                  zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
@@ -350 +354 @@
            IF (zrfl(i)+zifl(i).GT.0.) THEN
 !>AJ
+                ! Evap max pour ne pas saturer la fraction sous le nuage
                 zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
+                ! Calcul de l'evaporation du flux de precip herite
+                !   d'au-dessus
                 zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
                      * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
                 zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
                      * RG*dtime/(paprs(i,k)-paprs(i,k+1))
+                ! Seuil pour ne pas saturer la fraction sous le nuage
                 zqev = MIN (zqev, zqevt)
+                ! Nouveau flux de precip
                 zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
                      /RG/dtime
-       
-                ! pour la glace, on ré-évapore toute la précip dans la
-                ! couche du dessous
-                ! la glace venant de la couche du dessus est simplement
-                ! dans la couche du dessous.
-       
+                ! Aucun flux liquide pour T < t_coup
                 IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
-       
+                ! Nouvelle vapeur 
                 zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
                      * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+                ! Nouvelle temperature (chaleur latente)
                 zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
                      * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
@@ -384 +389 @@
 !>AJ
      !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-     ! Modification de l'évaporation avec la glace
-     ! Différentiation entre précipitation liquide et solide 
-     ! On suppose que coef_evai=2*coef_eva
+     ! Modification de l'evaporation avec la glace
+     ! Differentiation entre precipitation liquide et solide 
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
      
+     ! Evap max pour ne pas saturer la fraction sous le nuage
          zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
      !    zqev0 = MAX (0.0, zqs(i)-zq(i) ) 
 
      !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-     ! On différencie qsat pour l'eau et la glace
+     ! On differencie qsat pour l'eau et la glace
      ! Si zdelta=1. --> glace
      ! Si zdelta=0. --> eau liquide
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-         
+        
+         ! Calcul du qsat par rapport a l'eau liquide
          qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
          qsl= MIN(0.5,qsl)
@@ -403 +409 @@
          qsl= qsl*zcor
          
+         ! Calcul de l'evaporation du flux de precip herite
+         !   d'au-dessus
+         ! Formulation en racine du flux de precip
+         ! (Klemp & Wilhelmson, 1978; Sundqvist, 1988)
          zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
               *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
@@ -408 +418 @@
               *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 
 
+         
+         ! Calcul du qsat par rapport a la glace
          qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
          qsi= MIN(0.5,qsi)
@@ -413 +425 @@
          qsi= qsi*zcor
 
+         ! Calcul de la sublimation du flux de precip solide herite
+         !   d'au-dessus
          zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
               *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
@@ -420 +434 @@
               
      !JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-     ! Vérification sur l'évaporation
+     ! Verification sur l'evaporation
+     ! On s'assure qu'on ne sature pas 
+     ! la fraction sous le nuage sinon on
+     ! repartit zqev0 en gardant la proportion
+     ! liquide / glace
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      
@@ -436 +454 @@
              ENDIF
          ENDIF
-     
+         ! Nouveaux flux de precip liquide et solide 
          zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
                                  /RG/dtime)
          zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
                                  /RG/dtime)
-         
-     ! Pour la glace, on révapore toute la précip dans la couche du dessous
-     ! la glace venant de la couche du dessus est simplement dans la couche
-     ! du dessous.
-     
-     !    IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
-!         print*,zrfl(i),zrfln(i),zqevt,zqevti,RLMLT,'fluxdeprecip'
+
+         ! Mise a jour de la vapeur, temperature et flux de precip
          zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
                   * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
@@ -466 +479 @@
            zmelt = ((zt(i)-273.15)/(ztfondue-273.15))             ! jyg
            zmelt = MIN(MAX(zmelt,0.),1.)
+           ! Fusion de la glace
            zrfl(i)=zrfl(i)+zmelt*zifl(i)
            zifl(i)=zifl(i)*(1.-zmelt)
 !           print*,zt(i),'octavio1'
+           ! Chaleur latente de fusion
            zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
                       *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
@@ -481 +496 @@
       ENDIF ! (.NOT. ice_thermo) 
      
+     ! ----------------------------------------------------------------
+     ! Fin evaporation de la precipitation
+     ! ----------------------------------------------------------------
      ENDIF ! (evap_prec)
      !
@@ -518 +536 @@
 !       endif
 
+     ! ----------------------------------------------------------------
+     ! P2> Formation du nuage
+     ! ----------------------------------------------------------------
      IF (cpartiel) THEN
 
@@ -534 +555 @@
         if (iflag_pdf.eq.0) then
 
+           ! version creneau de (Li, 1998)
            do i=1,klon
               zdelq = min(ratqs(i,k),0.99) * zq(i)
@@ -575 +597 @@
            end if
 
-!CR: variation de qsat avec T en présence de glace ou non
+!CR: variation de qsat avec T en presence de glace ou non
 !initialisations
            do i=1,klon
@@ -587 +609 @@
 !Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
            if (iflag_fisrtilp_qsat.ge.0) then
+             ! Iteration pour condensation avec variation de qsat(T) 
+             ! -----------------------------------------------------
              do iter=1,iflag_fisrtilp_qsat+1
                 
@@ -603 +627 @@
                     endif
                  endif
+                 ! Calcul des PDF lognormales
                  zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
                  zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
@@ -647 +672 @@
 
                  else
-
-!calcul de la fraction de glace
-!CR: on utilise la nouvelle fonction de JBM pour l ancien calcul
-!                 zfice(i) = icefrac_lsc(Tbef(i), t_glace_min, &
-!                                     t_glace_max, exposant_glace)
-!                 zfice(i) = 1.0 - (Tbef(i)-ztglace) / (RTT-ztglace)
-!                 zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
-!                 zfice(i) = zfice(i)**nexpo
+                 ! Iteration pour convergence avec qsat(T)
                  if (iflag_t_glace.eq.1) then
                  CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
@@ -699 +717 @@
   
           
-             enddo
+             enddo ! iter=1,iflag_fisrtilp_qsat+1
+             ! Fin d'iteration pour condensation avec variation de qsat(T) 
+             ! -----------------------------------------------------------
            endif
 
@@ -707 +727 @@
 
 !        if (iflag_fisrtilp_qsat.eq.-1) then 
-!CR: ATTENTION option fausse mais a existe: pour la re-activer, prendre iflag_fisrtilp_qsat=0 et activer les lignes suivantes: 
+!------------------------------------------
+!CR: ATTENTION option fausse mais a existe:
+! pour la re-activer, prendre iflag_fisrtilp_qsat=0 et
+! activer les lignes suivantes: 
        IF (1.eq.0) THEN 
        DO i=1,klon
@@ -724 +747 @@
               rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
            ENDIF
-        ENDDO
-        ENDIF
+       ENDDO
+       ENDIF
+!------------------------------------------
 
 !        ELSE
 
+        ! Calcul de l'eau in-cloud (zqn),
+        ! moyenne dans la maille (zcond),
+        ! fraction nuageuse (rneb) et
+        ! humidite relative ciel-clair (rhcl)
         DO i=1,klon
            IF (rneb(i,k) .LE. 0.0) THEN
@@ -749 +777 @@
 !        ENDIF
 
-        !         do i=1,klon
-        !            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
-        !            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
-        !            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
-        !c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
-        !c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
-        !c  la convection.
-        !c  ATTENTION !!! Il va falloir verifier tout ca.
-        !            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
-        !c           print*,'ZDQS ',zdqs(i)
-        !c--Olivier
-        !            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
-        !            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
-        !            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
-        !c--fin
-        !           ENDDO
-     ELSE
+     ELSE ! de IF (cpartiel) 
+        ! Cas "tout ou rien"
         DO i = 1, klon
            IF (zq(i).GT.zqs(i)) THEN
@@ -775 +788 @@
         ENDDO
      ENDIF
-     !
+     ! ----------------------------------------------------------------
+     ! Fin de formation du nuage
+     ! ----------------------------------------------------------------
+     !
+     ! Mise a jour vapeur d'eau
      DO i = 1, klon
         zq(i) = zq(i) - zcond(i)
@@ -781 +798 @@
      ENDDO
 !AJ<
+     ! Chaleur latente apres formation nuage
+     ! -------------------------------------
      IF (.NOT. ice_thermo) THEN
         if (iflag_fisrtilp_qsat.lt.1) then
@@ -825 +844 @@
        ENDIF
 !>AJ
+     ! ----------------------------------------------------------------
+     ! P3> Formation des precipitations
+     ! ----------------------------------------------------------------
      !
      ! Partager l'eau condensee en precipitation et eau liquide nuageuse
      !
+
+     ! Initialisation de zoliq (eau condensee moyenne dans la maille)
      DO i = 1, klon
         IF (rneb(i,k).GT.0.0) THEN
@@ -858 +882 @@
        ENDIF
      ENDIF
+
+     ! Calcul de radliq (eau condensee pour le rayonnement)
+     ! Iteration pour realiser une moyenne de l'eau nuageuse lors de la precip
+     ! Remarque: ce n'est donc pas l'eau restante en fin de precip mais une
+     ! eau moyenne restante dans le nuage sur la duree du pas de temps qui est
+     ! transmise au rayonnement;
+     ! ----------------------------------------------------------------
      DO i = 1, klon
         IF (rneb(i,k).GT.0.0) THEN
@@ -878 +909 @@
                  ztot = 0.0
               ELSE
-                 !  quantite d'eau a eliminer: zchau
-                 !  meme chose pour la glace: zfroi
+                 !  quantite d'eau a eliminer: zchau (Sundqvist, 1978)
+                 !  meme chose pour la glace: zfroi (Zender & Kiehl, 1997)
                  if (ptconv(i,k)) then
                     zcl   =cld_lc_con
@@ -916 +947 @@
         ENDDO  ! i = 1,klon
      ENDDO     ! n = 1,ninter
+     ! ----------------------------------------------------------------
      !
      IF (.NOT. ice_thermo) THEN
@@ -1028 +1060 @@
      ELSE
      ! JAM*************************************************
-     ! Revoir partie ci-dessous: à quoi servent psfl et prfl?
+     ! Revoir partie ci-dessous: a quoi servent psfl et prfl?
      ! *****************************************************
 
@@ -1040 +1072 @@
        ENDDO
      ENDIF
+     ! ----------------------------------------------------------------
+     ! Fin de formation des precipitations
+     ! ----------------------------------------------------------------
      !
      !
@@ -1110 +1145 @@
      ENDDO
      !
-     !AA----------------------------------------------------------
-     !                     FIN DE BOUCLE SUR K   
+     !AA===============================================================
+     !                     FIN DE LA BOUCLE VERTICALE  
   end DO
   !
-  !AA-----------------------------------------------------------
+  !AA==================================================================
   !
   ! Pluie ou neige au sol selon la temperature de la 1ere couche