Changeset 2416

Timestamp:

Dec 24, 2015, 12:58:33 PM (9 years ago)

Author:

jyg

Message:

Accounting for situations where the adiabatic
ascent does not go above pLCL - 50 hPa in
cv3_cine.F90.

File:

• LMDZ5/trunk/libf/phylmd/cv3_cine.F90 (modified) (4 diffs)

LMDZ5/trunk/libf/phylmd/cv3_cine.F90

@@ -34 +34 @@
   INTEGER itop(nloc), ineg(nloc), ilow(nloc)
   INTEGER ifst(nloc), isublcl(nloc)
-  LOGICAL lswitch(nloc), lswitch1(nloc), lswitch2(nloc)
+  LOGICAL lswitch(nloc), lswitch1(nloc), lswitch2(nloc), lswitch3(nloc)
   LOGICAL exist_lfc(nloc)
   REAL dpmax
@@ -161 +161 @@
   END DO
 
+  ! 1.2.1 Recompute itop (=1st layer with positive buoyancy above ineg)
+  ! -------------------------------------------------------------------
+
+  DO il = 1, ncum
+    IF (lswitch(il)) THEN
+      itop(il) = nl - 1
+    END IF
+  END DO
+
+  DO k = nl, 1, -1
+    DO il = 1, ncum
+      IF (lswitch(il)) THEN
+        IF (k>=ineg(il) .AND. buoy(il,k)>0) THEN
+          itop(il) = k
+        END IF
+      END IF
+    END DO
+  END DO
+
+  ! If there is no layer with positive buoyancy above ineg, set Plfc, 
+  ! Cina and Cinb to arbitrary extreme values.
+  DO il = 1, ncum
+    IF (lswitch(il) .AND. itop(il) == nl - 1) THEN
+      plfc(il) = 1.121
+      cinb(il) = -1121.
+      cina(il) = -1122.
+    END IF
+  END DO
+
+  DO il = 1, ncum
+    lswitch3(il) = itop(il) < nl -1
+    lswitch(il) = lswitch1(il) .AND. lswitch2(il) .AND. lswitch3(il)
+  END DO
+
   DO il = 1, ncum
     IF (lswitch(il)) THEN
       cinb(il) = 0.
 
-      ! 1.2.1  Calcul de la pression du niveau de flot. nulle juste au-dessus
+      ! 1.2.2  Calcul de la pression du niveau de flot. nulle juste au-dessus
       ! de LCL
       ! ---------------------------------------------------------------------------
@@ -171 +205 @@
         ! In order to get P0, one may interpolate linearly buoyancies
         ! between P(ineg) and P(ineg-1).
-        p0(il) = (buoy(il,ineg(il))*p(il,ineg(il)-1)-buoy(il,ineg( &
-          il)-1)*p(il,ineg(il)))/(buoy(il,ineg(il))-buoy(il,ineg(il)-1))
+        p0(il) = (buoy(il,ineg(il))*p(il,ineg(il)-1)-buoy(il,ineg(il)-1)*p(il,ineg(il)))/ &
+          (buoy(il,ineg(il))-buoy(il,ineg(il)-1))
       ELSE
         ! In order to get P0, one has to interpolate between P(ineg) and
@@ -180 +214 @@
       END IF
     END IF
-  END DO
-
-  ! 1.2.2 Recompute itop (=1st layer with positive buoyancy above ineg)
-  ! -------------------------------------------------------------------
-  DO il = 1, ncum
-    IF (lswitch(il)) THEN
-      itop(il) = nl - 1
-    END IF
-  END DO
-
-  DO k = nl, 1, -1
-    DO il = 1, ncum
-      IF (lswitch(il)) THEN
-        IF (k>=ineg(il) .AND. buoy(il,k)>0) THEN
-          itop(il) = k
-        END IF
-      END IF
-    END DO
   END DO