Changeset 2420 for LMDZ5/trunk/libf/phylmd/cv3p2_closure.F90

Timestamp:

Jan 5, 2016, 4:37:49 PM (8 years ago)

Author:

crio

Message:

Nouvelle option d'epluchage de l'ascendance adiabatique dans le schema d'Emanuel: epluchage fonction de B/w2 au lieu de w. S'active avec iflag_mix_adiab=1 (valeur par defaut iflag_mix_adiab=0). Fonctionne avec iflag_mix=0 et iflag_mix=1.
Correction de bugs dans le schema de convection pour le calcul de inb, cape et buoy (sous le meme flag pour l'instant).
New option for the erosion of the adiabatic ascent in the Emanuel scheme: erosion function of B/w2 instead of w. Activated by iflag_mix_adiab=1 (standard value iflag_mix_adiab=0). Should work with iflag_mix=0 and iflag_mix=1.
Various bug corrections in the convection scheme for the computation of inb, cape, buoy (protected by the same flag for now).

File:

• LMDZ5/trunk/libf/phylmd/cv3p2_closure.F90 (modified) (5 diffs)

LMDZ5/trunk/libf/phylmd/cv3p2_closure.F90

@@ -59 +59 @@
   REAL                                               :: deltap, fac, w, amu
   REAL, DIMENSION (nloc, nd)                         :: rhodp               ! Factor such that m=rhodp*sig*w
+  REAL                                               :: dz
   REAL                                               :: pbmxup
   REAL, DIMENSION (nloc, nd)                         :: dtmin, sigold
@@ -85 +86 @@
   REAL, DIMENSION (nloc)                             :: alp2                  ! Alp with offset
 
+!CR: variables for new erosion of adiabiatic ascent
+  REAL, DIMENSION (nloc, nd)                         :: mad, me, betalim, beta_coef
+  REAL, DIMENSION (nloc, nd)                         :: med, md
+  REAL                                               :: coef_peel
+  PARAMETER (coef_peel=0.25)
+
   REAL                                               :: sigmax
   PARAMETER (sigmax=0.1)
@@ -120 +127 @@
     END DO
   END DO
+
+!CR: initializations for erosion of adiabatic ascent 
+  DO k = 1,nl
+    DO il = 1, ncum
+        mad(il,k)=0.
+        me(il,k)=0.
+        betalim(il,k)=1.
+        wlim(il,k)=0.
+    ENDDO
+  ENDDO
 
   ! -------------------------------------------------------
@@ -496 +513 @@
       IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
 
+        IF (iflag_mix_adiab.eq.1) THEN
+!CR:computation of cape from LCL: keep flag or to modify in all cases?
+        deltap = min(plcl(il), ph(il,k-1)) - min(plcl(il), ph(il,k))
+        ELSE
         deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
+        ENDIF
         cape(il) = cape(il) + rrd*buoy(il, k-1)*deltap/p(il, k-1)
         cape(il) = amax1(0.0, cape(il))
@@ -681 +703 @@
                          (k,w0(igout,k),sig(igout,k), k=icb(igout),inb(igout))
 
+!CR: new erosion of adiabatic ascent: modification of m
+!computation of the sum of ascending fluxes 
+  IF (iflag_mix_adiab.eq.1) THEN
+
+!Verification sum(me)=sum(m)
+  DO k = 1,nl+1
+    DO il = 1, ncum
+       md(il,k)=0.
+       med(il,k)=0.
+    ENDDO
+  ENDDO
+
+  DO k = nl,1,-1
+    DO il = 1, ncum
+           md(il,k)=md(il,k+1)+m(il,k+1)
+    ENDDO
+  ENDDO
+
+  DO k = nl,1,-1
+    DO il = 1, ncum
+        IF ((k>=(icb(il))) .AND. (k<=inb(il))) THEN
+           mad(il,k)=mad(il,k+1)+m(il,k+1)
+        ENDIF
+!        print*,"mad",il,k,mad(il,k)
+    ENDDO
+  ENDDO
+
+!CR: erosion of each adiabatic ascent during its ascent
+
+!Computation of erosion coefficient beta_coef
+  DO k = 1, nl
+    DO il = 1, ncum
+       IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k).gt.0.)) THEN     
+!          print*,"beta_coef",il,k,icb(il),inb(il),buoy(il,k),tv(il,k),wlim(il,k),wlim(il,k+1)
+          beta_coef(il,k)=RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2
+       ELSE
+          beta_coef(il,k)=0.
+       ENDIF
+    ENDDO
+  ENDDO
+
+!  print*,"apres beta_coef"
+
+  DO k = 1, nl
+    DO il = 1, ncum
+
+      IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
+
+!        print*,"dz",il,k,tv(il, k-1)
+        dz = (ph(il,k-1)-ph(il,k))/(p(il, k-1)/(rrd*tv(il, k-1))*RG)
+        betalim(il,k)=betalim(il,k-1)*exp(-1.*beta_coef(il,k-1)*dz)
+!        betalim(il,k)=betalim(il,k-1)*exp(-RG*coef_peel*buoy(il,k-1)/tv(il,k-1)/5.**2*dz)
+!        print*,"me",il,k,mlim(il,k),buoy(il,k),wlim(il,k),mad(il,k)
+        dz = (ph(il,k)-ph(il,k+1))/(p(il, k)/(rrd*tv(il, k))*RG)
+!        me(il,k)=betalim(il,k)*(m(il,k)+RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz*mad(il,k))
+        me(il,k)=betalim(il,k)*(m(il,k)+beta_coef(il,k)*dz*mad(il,k))
+!        print*,"B/w2",il,k,RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz    
+      
+      END IF
+        
+!Modification of m
+      m(il,k)=me(il,k) 
+    END DO
+  END DO
+ 
+!  DO il = 1, ncum
+!     dz = (ph(il,icb(il))-ph(il,icb(il)+1))/(p(il, icb(il))/(rrd*tv(il, icb(il)))*RG)
+!     m(il,icb(il))=m(il,icb(il))+RG*coef_peel*buoy(il,icb(il))/tv(il,icb(il)) &
+!                  /((wlim(il,icb(il))+wlim(il,icb(il)+1))/2.)**2*dz*mad(il,icb(il))
+!     print*,"wlim(icb)",icb(il),wlim(il,icb(il)),m(il,icb(il))
+!  ENDDO
+
+!Verification sum(me)=sum(m)
+  DO k = nl,1,-1
+    DO il = 1, ncum
+           med(il,k)=med(il,k+1)+m(il,k+1)
+!           print*,"somme(me),somme(m)",il,k,icb(il),med(il,k),md(il,k),me(il,k),m(il,k),wlim(il,k)
+    ENDDO
+  ENDDO
+
+
+  ENDIF !(iflag_mix_adiab)
+!RC
+
   ! c 3. Compute final cloud base mass flux;
   ! c    set iflag to 3 if cloud base mass flux is exceedingly small and is