Changeset 5544

Timestamp:

Feb 17, 2025, 1:29:29 PM (3 days ago)

Author:

jyg

Message:

Bug fix in cv3_routines.f90 (in cv3_unsat and cv3_yield) concerning water conservation.
And small improvement in the coding of ice fraction in convective precipitation.

File:

• LMDZ6/trunk/libf/phylmd/cv3_routines.f90 (modified) (12 diffs)

LMDZ6/trunk/libf/phylmd/cv3_routines.f90

@@ -2907 +2907 @@
     DO il = 1, ncum
       IF (i<=inb(il) .AND. lwork(il)) THEN
-        wdtrain(il) = grav*ep(il, i)*m(il, i)*clw(il, i)
-        wdtrainS(il, i) = wdtrain(il)/grav                                            !   Ps   jyg
-!!        wdtrainA(il, i) = wdtrain(il)/grav                                          !   Ps   RomP
+        wdtrainS(il, i) = ep(il, i)*m(il, i)*clw(il, i)                               ! jyg
+!!jyg        wdtrain(il) = grav*ep(il, i)*m(il, i)*clw(il, i)
+!!jyg        wdtrainS(il, i) = wdtrain(il)/grav                                            !   Ps   jyg
+!!jyg        wdtrainA(il, i) = wdtrain(il)/grav                                          !   Ps   RomP
       END IF
     END DO
@@ -2919 +2920 @@
             awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
             awat = max(awat, 0.0)
-            wdtrain(il) = wdtrain(il) + grav*awat*ment(il, j, i)
-            wdtrainM(il, i) = wdtrain(il)/grav - wdtrainS(il, i)    !   Pm  jyg
-!!            wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i)  !   Pm  RomP
+            wdtrainM(il, i) = wdtrainM(il, i) + awat*ment(il, j, i)                   ! jyg
+!!jyg            wdtrain(il) = wdtrain(il) + grav*awat*ment(il, j, i)
+!!jyg            wdtrainM(il, i) = wdtrain(il)/grav - wdtrainS(il, i)    !   Pm  jyg
+!!jyg            wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i)  !   Pm  RomP
           END IF
         END DO
@@ -2938 +2940 @@
               IF (i<=inb(il) .AND. lwork(il)) THEN
                 wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i))/(1. - qta(il, i-1))    !   Pa   jygprl
-                wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i)
+!!jyg                wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i)
               END IF
             END DO
@@ -2949 +2951 @@
               IF (i<=inb(il) .AND. lwork(il)) THEN
                 wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i))                        !   Pa   jygprl
-                wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i)
+!!jyg                wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i)
               END IF
             END DO
@@ -2958 +2960 @@
     ENDIF  ! (cvflag_prec_eject)
 
+    IF ( i > 1) THEN
+      DO il = 1, ncum
+        IF (i<=inb(il) .AND. lwork(il)) THEN
+          wdtrain(il) = grav*(wdtrainS(il,i) + wdtrainM(il,i) + wdtrainA(il,i))
+        END IF
+      END DO
+    ENDIF  ! ( i > 1)
 
 ! ***    find rain water and evaporation using provisional   ***
@@ -3139 +3148 @@
           ice(il, i) = ice(il, i) - fondue(il, i)
 
-          IF (water(il,i)+ice(il,i)<1.E-30) THEN
-            faci(il, i) = 0.
-          ELSE
-            faci(il, i) = ice(il, i)/(water(il,i)+ice(il,i))
-          END IF
+!!          IF (water(il,i)+ice(il,i)<1.E-30) THEN
+!!            faci(il, i) = 0.
+!!          ELSE
+!!            faci(il, i) = ice(il, i)/(water(il,i)+ice(il,i))
+!!          END IF
+
+            faci(il,i) = ice(il, i)/max((water(il,i)+ice(il,i)), smallestreal)
 
 !           water(il,i)=water(il,i+1)+(1.-fraci(il,i))*e6+(1.-faci(il,i))*f6
@@ -3423 +3434 @@
 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
-
   RETURN
 
@@ -3531 +3541 @@
       REAL, DIMENSION (nloc, nd)                         :: sigment, qtment             ! cld
       REAL, DIMENSION (nloc, nd, nd)                     :: qdet
-      REAL sumdq !jyg
+!!      REAL sumdq !jyg
 ! 
 ! -------------------------------------------------------------
@@ -4007 +4017 @@
 ! ***                      through each level                          ***
 
-
 !jyg<
 !!  DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
@@ -4026 +4035 @@
 DO il = 1, ncum
    amp1(il) = upwd(il,i+1)
-   ad(il) = dnwd(il,i)
+!!jyg20250215   ad(il) = dnwd(il,i)
+   ad(il) = -dnwd(il,i)
 ENDDO
 !-----------------------------------------------------------
@@ -4360 +4370 @@
 500 END DO
 
-!JYG<
-!Conservation de l'eau
-!   sumdq = 0.
-!   DO k = 1, nl
-!     sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
-!   END DO
-!   PRINT *, 'cv3_yield, apres 500, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
-!JYG>
+!!!JYG<
+!!!Conservation de l'eau
+!!   sumdq = 0.
+!!   DO k = 1, nl
+!!     sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
+!!   END DO
+!!   PRINT *, 'cv3_yield, apres 500, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
+!!!JYG>
 ! ***   move the detrainment at level inb down to level inb-1   ***
 ! ***        in such a way as to preserve the vertically        ***
@@ -4402 +4412 @@
   END DO
 
-!JYG<
-!Conservation de l'eau
-!   sumdq = 0.
-!   DO k = 1, nl
-!     sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
-!   END DO
-!   PRINT *, 'cv3_yield, apres 503, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
-!JYG>
+!!!JYG<
+!!!Conservation de l'eau
+!!   sumdq = 0.
+!!   DO k = 1, nl
+!!     sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
+!!   END DO
+!!   PRINT *, 'cv3_yield, apres 503, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
+!!!JYG>
 
 !AC!      do j=1,ntra