Changeset 4887 for LMDZ6/trunk/libf/phylmd/physiq_mod.F90

Timestamp:

Apr 2, 2024, 1:14:33 PM (7 weeks ago)

Author:

oboucher

Message:

Correction of another small bug on the water mass fixer (not activated)
Introduction of a water_budget diagnostic that quantifies non-conservation from the LMDZ physics

File:

• LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (modified) (4 diffs)

LMDZ6/trunk/libf/phylmd/physiq_mod.F90

@@ -259 +259 @@
        JrNt,                             &
        dthmin, evap, snowerosion,fder, plcl, plfc,   &
-       prw, prlw, prsw, prbsw,                  &
+       prw, prlw, prsw, prbsw, water_budget,         &
        s_lcl, s_pblh, s_pblt, s_therm,   &
        cdragm, cdragh,                   &
@@ -2464 +2464 @@
     qql1(:)=0.0
     DO k = 1, klev
-      qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k)+qbs_seri(:,k))*zmasse(:,k)
+      qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k))*zmasse(:,k)
+      IF (nqo >= 3) THEN
+        qql1(:)=qql1(:)+qs_seri(:,k)*zmasse(:,k)
+      ENDIF
+      IF (ok_bs) THEN
+        qql1(:)=qql1(:)+qbs_seri(:,k)*zmasse(:,k)
+      ENDIF
     ENDDO
     ENDIF
@@ -5465 +5471 @@
     t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
 
-    !=======================================================================
-    !   SORTIES
-    !=======================================================================
-    !
-    !IM initialisation + calculs divers diag AMIP2
-    !
-    include "calcul_divers.h"
-    !
-    !IM Interpolation sur les niveaux de pression du NMC
-    !   -------------------------------------------------
-    !
-    include "calcul_STDlev.h"
-    !
-    ! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
-    CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
-    !
+    !==================================================================
+    !--OB water mass fixer for the physics
+    !--water profiles are corrected to force mass conservation of water
+    !--currently flag is turned off
+    !==================================================================
+    IF (mass_fixer) THEN
+    qql2(:)=0.0
+    DO k = 1, klev
+      qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k))*zmasse(:,k)
+      IF (nqo >= 3) THEN
+        qql2(:)=qql2(:)+qs_seri(:,k)*zmasse(:,k)
+      ENDIF
+      IF (ok_bs) THEN
+        qql2(:)=qql2(:)+qbs_seri(:,k)*zmasse(:,k)
+      ENDIF
+    ENDDO
+
+#ifdef CPP_StratAer
+    IF (ok_qemiss) THEN
+       DO k = 1, klev
+          qql1(:) = qql1(:)+d_q_emiss(:,k)*zmasse(:,k)
+       ENDDO
+    ENDIF
+#endif
+    IF (ok_qch4) THEN
+       DO k = 1, klev
+          qql1(:) = qql1(:)+d_q_ch4_dtime(:,k)*zmasse(:,k)
+       ENDDO
+    ENDIF
+    
+    DO i = 1, klon
+      !--compute ratio of what q+ql should be with conservation to what it is
+      IF (ok_bs) THEN
+        corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i)-bs_fall(i))*pdtphys)/qql2(i)
+      ELSE
+        corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
+      ENDIF
+      DO k = 1, klev
+        q_seri(i,k) =q_seri(i,k)*corrqql
+        ql_seri(i,k)=ql_seri(i,k)*corrqql
+        IF (nqo >= 3) THEN
+          qs_seri(i,k)=qs_seri(i,k)*corrqql
+        ENDIF
+        IF (ok_bs) THEN
+          qbs_seri(i,k)=qbs_seri(i,k)*corrqql
+        ENDIF
+      ENDDO
+    ENDDO
+    ENDIF
+    !--fin mass fixer
+
     !cc prw  = eau precipitable
     !   prlw = colonne eau liquide 
     !   prlw = colonne eau solide
     !   prbsw = colonne neige soufflee
+    !   water_budget = non-conservation residual from the LMDZ physics
+    !                  (should be equal to machine precision if mass fixer is activated)
     prw(:) = 0.
     prlw(:) = 0.
     prsw(:) = 0.
     prbsw(:) = 0.
+    water_budget(:) = 0.0
     DO k = 1, klev
        prw(:)  = prw(:)  + q_seri(:,k)*zmasse(:,k)
        prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
-       prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
-       prbsw(:)= prbsw(:) + qbs_seri(:,k)*zmasse(:,k)
+       water_budget(:) = water_budget(:) + (q_seri(:,k)-qx(:,k,ivap)+ql_seri(:,k)-qx(:,k,iliq))*zmasse(:,k)
+       IF (nqo >= 3) THEN
+         prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
+         water_budget(:) = water_budget(:) + (qs_seri(:,k)-qx(:,k,isol))*zmasse(:,k)
+       ENDIF
+       IF (nqo >= 4 .AND. ok_bs) THEN
+         prbsw(:)= prbsw(:) + qbs_seri(:,k)*zmasse(:,k)
+         water_budget(:) = water_budget(:) + (qbs_seri(:,k)-qx(:,k,ibs))*zmasse(:,k)
+       ENDIF
     ENDDO
+    water_budget(:)=water_budget(:)+(rain_fall(:)+snow_fall(:)-evap(:))*pdtphys
+    IF (ok_bs) THEN
+      water_budget(:)=water_budget(:)+bs_fall(:)*pdtphys
+    ENDIF
+
+    !=======================================================================
+    !   SORTIES
+    !=======================================================================
+    !
+    !IM initialisation + calculs divers diag AMIP2
+    !
+    include "calcul_divers.h"
+    !
+    !IM Interpolation sur les niveaux de pression du NMC
+    !   -------------------------------------------------
+    !
+    include "calcul_STDlev.h"
+    !
+    ! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
+    CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
+    !
     !
     IF (ANY(type_trac == ['inca','inco'])) THEN
@@ -5579 +5651 @@
     !IM global posePB      include "write_bilKP_ave.h"
     !
-
-    !--OB mass fixer 
-    !--profile is corrected to force mass conservation of water
-    IF (mass_fixer) THEN
-    qql2(:)=0.0
-    DO k = 1, klev
-      qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k)+qbs_seri(:,k))*zmasse(:,k)
-    ENDDO
-
-#ifdef CPP_StratAer
-    IF (ok_qemiss) THEN
-       DO k = 1, klev
-          qql1(:) = qql1(:)+d_q_emiss(:,k)*zmasse(:,k)
-       ENDDO
-    ENDIF
-#endif
-    IF (ok_qch4) THEN
-       DO k = 1, klev
-          qql1(:) = qql1(:)+d_q_ch4_dtime(:,k)*zmasse(:,k)
-       ENDDO
-    ENDIF
-    
-    DO i = 1, klon
-      !--compute ratio of what q+ql should be with conservation to what it is
-      corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
-      DO k = 1, klev
-        q_seri(i,k) =q_seri(i,k)*corrqql
-        ql_seri(i,k)=ql_seri(i,k)*corrqql
-        qs_seri(i,k)=qs_seri(i,k)*corrqql
-        qbs_seri(i,k)=qbs_seri(i,k)*corrqql
-      ENDDO
-    ENDDO
-    ENDIF
-    !--fin mass fixer
-
     ! Sauvegarder les valeurs de t et q a la fin de la physique:
     !