Changeset 1901 for LMDZ5/trunk

Timestamp:

Nov 19, 2013, 10:29:06 AM (11 years ago)

Author:

idelkadi

Message:

Prise en compte de la variation de l'humidite a saturation (qsat) avec la
variation de temperature associee a la condensation dans le calcul de
l'eau nuageuse dans fisrtilp. Introduction d'une boucle de convergence
avec un nombre maximum d'iterations fixe a iflag_fisrtilp_qsat.
Si iflag_fisrtilp_qsat=0: calcul sans variation de qsat avec T
Si iflag_fisrtilp_qsat=-1: calcul bogue

File:

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

LMDZ5/trunk/libf/phylmd/fisrtilp.F90

@@ -100 +100 @@
   !
   INTEGER i, k, n, kk
-  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5   
+  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5  
+  REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
+  LOGICAL convergence(klon)
+  REAL DDT0
+  PARAMETER (DDT0=.01)
+  INTEGER n_i, iter
+ 
   REAL zrfl(klon), zrfln(klon), zqev, zqevt 
   REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
@@ -467 +473 @@
      ! de l'eau condensee:
      !
+!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
+       if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
+         write(*,*) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
+        " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
+         stop
+       endif
 
      IF (cpartiel) THEN
@@ -526 +538 @@
 
            do i=1,klon
+              Tbef(i)=zt(i)
+              qlbef(i)=0. 
               if (lognormale(i)) then
                  zpdf_sig(i)=ratqs(i,k)*zq(i)
@@ -538 +552 @@
                  zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
                  zpdf_e2(i)=1.-erf(zpdf_e2(i))
-              endif
-           enddo
-
-           do i=1,klon
-              if (lognormale(i)) then
+             
                  if (zpdf_e1(i).lt.1.e-10) then
                     rneb(i,k)=0.
@@ -550 +560 @@
                     zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
                  endif
+
+                 qlbef(i)=max(0.,zqn(i)-zqs(i))
+                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
+                 denom=1.+rneb(i,k)*zdqs(i)
+                 DT(i)=num/denom
               endif
-
            enddo
+  
+           n_i=1
+!               do while ((abs(DT(i)).gt.DDT0).and.((zqn(i)-zqs(i)).gt.0.))
+! iflag_fisrtilp_qsat=0: qsat ne varie pas avec T
+! iflag_fisrtilp_qsat > 1 : nombre d iterations max pour converger sur le calcul de qsat(T)
+!               do while (n_i.le.iflag_fisrtilp_qsat)
+           if (iflag_fisrtilp_qsat.ge.1) then
+           do iter=1,iflag_fisrtilp_qsat
+              do i=1,klon
+                 convergence(i)=abs(DT(i)).gt.DDT0
+                 if (convergence(i).and.lognormale(i)) then
+                 Tbef(i)=Tbef(i)+DT(i)                 
+                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(i)))
+                 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
+                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
+                 zqs(i)= R2ES * FOEEW(Tbef(i),zdelta)/pplay(i,k)
+                 zqs(i)=MIN(0.5,zqs(i))
+                 zcor=1./(1.-retv*zqs(i))
+                 zqs(i)=zqs(i)*zcor
+                 zdqs(i)=FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
+
+                 zpdf_sig(i)=ratqs(i,k)*zq(i)
+                 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
+                 zpdf_delta(i)=log(zq(i)/zqs(i))
+                 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
+                 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
+                 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
+                 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
+                 zpdf_e1(i)=1.-erf(zpdf_e1(i))
+                 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
+                 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
+                 zpdf_e2(i)=1.-erf(zpdf_e2(i))
+              
+                 if (zpdf_e1(i).lt.1.e-10) then
+                    rneb(i,k)=0.
+                    zqn(i)=zqs(i)
+                 else
+                    rneb(i,k)=0.5*zpdf_e1(i)
+                    zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
+                 endif
+
+                 qlbef(i)=max(0.,zqn(i)-zqs(i))    
+                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
+                 denom=1.+rneb(i,k)*zdqs(i)
+                 DT(i)=num/denom
+                 n_i=n_i+1
+                 endif
+              enddo
+              enddo
+
+           endif
 
 
         endif ! iflag_pdf
 
+        if (iflag_fisrtilp_qsat.eq.-1) then 
+!CR: ATTENTION option fausse mais a existe
         DO i=1,klon
            IF (rneb(i,k) .LE. 0.0) THEN
@@ -565 +632 @@
            ELSE IF (rneb(i,k) .GE. 1.0) THEN
               zqn(i) = zq(i)
-              rneb(i,k) = 1.0                  
-              zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1.+iflag_fisrtilp_qsat*zdqs(i))
+              rneb(i,k) = 1.0                 
+              zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
               rhcl(i,k)=1.0
            ELSE
-              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+iflag_fisrtilp_qsat*zdqs(i))
+              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
               rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
            ENDIF
         ENDDO
+
+        ELSE
+
+        DO i=1,klon
+           IF (rneb(i,k) .LE. 0.0) THEN
+              zqn(i) = 0.0
+              rneb(i,k) = 0.0
+              zcond(i) = 0.0
+              rhcl(i,k)=zq(i)/zqs(i)
+           ELSE IF (rneb(i,k) .GE. 1.0) THEN
+              zqn(i) = zq(i)
+              rneb(i,k) = 1.0
+              zcond(i) = MAX(0.0,zqn(i)-zqs(i))
+              rhcl(i,k)=1.0
+           ELSE
+              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
+              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
+           ENDIF
+        ENDDO
+
+
+        ENDIF
+
         !         do i=1,klon
         !            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
@@ -606 +696 @@
 !AJ<
      IF (.NOT. ice_thermo) THEN
-       DO i = 1, klon
-         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
-       ENDDO
+        if (iflag_fisrtilp_qsat.lt.1) then
+           DO i = 1, klon
+             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
+           ENDDO
+        else if (iflag_fisrtilp_qsat.gt.0) then
+           DO i= 1, klon
+             if (lognormale(i)) then
+             zt(i)=Tbef(i)
+             else
+             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
+             endif
+           ENDDO
+        endif
      ELSE
-       DO i = 1, klon
-         zfice(i) = 1.0 - (zt(i)-ztglace) / (273.15-ztglace)
-         zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
-         zfice(i) = zfice(i)**nexpo
-         zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
+         if (iflag_fisrtilp_qsat.lt.1) then
+           DO i = 1, klon
+              zfice(i) = 1.0 - (zt(i)-ztglace) / (273.15-ztglace)
+              zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
+              zfice(i) = zfice(i)**nexpo
+              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
                        +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
+           ENDDO
+         else
+           DO i=1, klon
+              zfice(i) = 1.0 - (zt(i)-ztglace) / (273.15-ztglace)
+              zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
+              zfice(i) = zfice(i)**nexpo
+!CR: ATTENTION zt different de Tbef: à corriger
+              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
+                       +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
+           ENDDO
+         endif
 !         print*,zt(i),zrfl(i),zifl(i),'temp1'
-       ENDDO
      ENDIF
 !>AJ