Changeset 5224 for LMDZ6/branches/Amaury_dev/libf/phylmd/physiq_mod.F90

Timestamp:

Sep 24, 2024, 10:47:17 AM (4 weeks ago)

Author:

abarral

Message:

Merge r5204 r5205
Light lint
Correct missing IOIPSL includes

File:

• LMDZ6/branches/Amaury_dev/libf/phylmd/physiq_mod.F90 (modified) (20 diffs)

LMDZ6/branches/Amaury_dev/libf/phylmd/physiq_mod.F90

@@ -1 @@
-! $Id$
-
-!#define IO_DEBUG
 MODULE physiq_mod
   IMPLICIT NONE
@@ -70 +67 @@
     USE tracinca_mod, ONLY: config_inca
     USE tropopause_m, ONLY: dyn_tropopause
-    USE ice_sursat_mod, ONLY: flight_init, airplane
     USE lmdz_vampir
     USE lmdz_writefield_phy
@@ -139 +135 @@
             ! [Variables internes non sauvegardees de la physique]
             ! Variables locales pour effectuer les appels en serie
-            t_seri, q_seri, ql_seri, qs_seri, qbs_seri, u_seri, v_seri, tr_seri, rneb_seri, &
+            t_seri,q_seri,ql_seri,qs_seri,qbs_seri, &
+            u_seri,v_seri,cf_seri,rvc_seri,tr_seri, &
             rhcl, &
             ! Dynamic tendencies (diagnostics)
-            d_t_dyn, d_q_dyn, d_ql_dyn, d_qs_dyn, d_qbs_dyn, d_u_dyn, d_v_dyn, d_tr_dyn, d_rneb_dyn, &
+            d_t_dyn, d_q_dyn, d_ql_dyn, d_qs_dyn, d_qbs_dyn, &
+            d_u_dyn, d_v_dyn,d_cf_dyn,d_rvc_dyn,d_tr_dyn, &
             d_q_dyn2d, d_ql_dyn2d, d_qs_dyn2d, d_qbs_dyn2d, &
             ! Physic tendencies
@@ -314 +312 @@
             pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb, &
             distcltop, temp_cltop, &
-            zqsatl, zqsats, &
-            qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
+            !-- LSCP - condensation and ice supersaturation variables
+            qsub, qissr, qcld, subfra, issrfra, gamma_cond, ratio_qi_qtot, &
+            dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
+            dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
+            !-- LSCP - aviation and contrails variables
             Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
+            dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
+            !
             cldemi, &
             cldfra, cldtau, fiwc, &
@@ -496 +499 @@
 
     ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional), blowing snow (optional)
-    INTEGER, SAVE :: ivap, iliq, isol, irneb, ibs
-    !$OMP THREADPRIVATE(ivap, iliq, isol, irneb, ibs)
+    INTEGER, SAVE :: ivap, iliq, isol, ibs, icf, irvc
+    !$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, irvc)
 
 
@@ -1317 +1320 @@
       iliq = strIdx(tracers(:)%name, addPhase('H2O', 'l'))
       isol = strIdx(tracers(:)%name, addPhase('H2O', 's'))
-      irneb = strIdx(tracers(:)%name, addPhase('H2O', 'r'))
       ibs = strIdx(tracers(:)%name, addPhase('H2O', 'b'))
+      icf = strIdx(tracers(:)%name, addPhase('H2O', 'f'))
+      irvc = strIdx(tracers(:)%name, addPhase('H2O', 'c'))
       !       CALL init_etat0_limit_unstruct
       !       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
@@ -1360 +1364 @@
       ENDIF
 
-      IF (ok_ice_sursat.AND.(iflag_ice_thermo==0)) THEN
-        WRITE (lunout, *) ' ok_ice_sursat=y requires iflag_ice_thermo=1 as well'
+      IF (ok_ice_supersat.AND.(iflag_ice_thermo==0)) THEN
+        WRITE (lunout, *) ' ok_ice_supersat=y requires iflag_ice_thermo=1 as well'
         abort_message = 'see above'
         CALL abort_physic(modname, abort_message, 1)
       ENDIF
 
-      IF (ok_ice_sursat.AND.(nqo<4)) THEN
-        WRITE (lunout, *) ' ok_ice_sursat=y requires 4 H2O tracers ', &
-                '(H2O_g, H2O_l, H2O_s, H2O_r) but nqo=', nqo, '. Might as well stop here.'
+      IF (ok_ice_supersat.AND.(nqo<5)) THEN
+        WRITE (lunout, *) ' ok_ice_supersat=y requires 5 H2O tracers ', &
+                '(H2O_g, H2O_l, H2O_s, H2O_f, H2O_c) but nqo=', nqo, '. Might as well stop here.'
         abort_message = 'see above'
         CALL abort_physic(modname, abort_message, 1)
       ENDIF
 
-      IF (ok_plane_h2o.AND..NOT.ok_ice_sursat) THEN
-        WRITE (lunout, *) ' ok_plane_h2o=y requires ok_ice_sursat=y '
+      IF (ok_plane_h2o.AND..NOT.ok_ice_supersat) THEN
+        WRITE (lunout, *) ' ok_plane_h2o=y requires ok_ice_supersat=y '
         abort_message = 'see above'
         CALL abort_physic(modname, abort_message, 1)
       ENDIF
 
-      IF (ok_plane_contrail.AND..NOT.ok_ice_sursat) THEN
-        WRITE (lunout, *) ' ok_plane_contrail=y requires ok_ice_sursat=y '
+      IF (ok_plane_contrail.AND..NOT.ok_ice_supersat) THEN
+        WRITE (lunout, *) ' ok_plane_contrail=y requires ok_ice_supersat=y '
         abort_message = 'see above'
         CALL abort_physic(modname, abort_message, 1)
@@ -1386 +1390 @@
 
       IF (ok_bs) THEN
-        IF ((ok_ice_sursat.AND.nqo <5).OR.(.NOT.ok_ice_sursat.AND.nqo<4)) THEN
+        IF ((ok_ice_supersat.AND.nqo <6).OR.(.NOT.ok_ice_supersat.AND.nqo<4)) THEN
           WRITE (lunout, *) 'activation of blowing snow needs a specific H2O tracer', &
                   'but nqo=', nqo
@@ -1795 +1799 @@
               RG, RD, RCPD, RKAPPA, RLVTT, RETV)
       CALL ratqs_ini(klon, klev, iflag_thermals, lunout, nbsrf, is_lic, is_ter, RG, RV, RD, RCPD, RLSTT, RLVTT, RTT)
-      CALL lscp_ini(pdtphys, lunout, prt_level, ok_ice_sursat, iflag_ratqs, fl_cor_ebil, RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI)
+      CALL lscp_ini(pdtphys, lunout, prt_level, ok_ice_supersat, iflag_ratqs, fl_cor_ebil, &
+              RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RV, RG, RPI, EPS_W)
       CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, &
               RVTMP2, RTT, RD, RG, RV, RPI)
@@ -2265 +2270 @@
               sollwdown(:))
 
+      !--Init for LSCP - condensation
+      ratio_qi_qtot(:,:) = 0.
+
     ENDIF
 
@@ -2361 +2369 @@
         ql_seri(i, k) = qx(i, k, iliq)
         qbs_seri(i, k) = 0.
+        cf_seri(i, k) = 0.
+        rvc_seri(i, k) = 0.
         !CR: ATTENTION, on rajoute la variable glace
         IF (nqo==2) THEN             !--vapour and liquid only
           qs_seri(i, k) = 0.
-          rneb_seri(i, k) = 0.
         ELSE IF (nqo==3) THEN        !--vapour, liquid and ice
           qs_seri(i, k) = qx(i, k, isol)
-          rneb_seri(i, k) = 0.
-        ELSE IF (nqo>=4) THEN        !--vapour, liquid, ice and rneb and blowing snow
+        ELSE IF (nqo>=4) THEN        !--vapour, liquid, ice, blowing snow, cloud fraction and cloudy water vapor to total water vapor ratio
           qs_seri(i, k) = qx(i, k, isol)
-          IF (ok_ice_sursat) THEN
-            rneb_seri(i, k) = qx(i, k, irneb)
-          ENDIF
+          IF (ok_ice_supersat) THEN
+            cf_seri(i, k) = qx(i, k, icf)
+            rvc_seri(i, k) = qx(i, k, irvc)
+          END IF
           IF (ok_bs) THEN
             qbs_seri(i, k) = qx(i, k, ibs)
-          ENDIF
-        ENDIF
-      ENDDO
-    ENDDO
+          END IF
+        END IF
+      END DO
+    END DO
     ! Lea Raillard qs_ini for cloud phase param.
     qs_ini(:, :) = qs_seri(:, :)
@@ -2450 +2459 @@
       d_qs_dyn(:, :) = (qs_seri(:, :) - qs_ancien(:, :)) / phys_tstep
       d_qbs_dyn(:, :) = (qbs_seri(:, :) - qbs_ancien(:, :)) / phys_tstep
+      d_cf_dyn(:, :) = (cf_seri(:, :) - cf_ancien(:, :)) / phys_tstep
+      d_rvc_dyn(:, :) = (rvc_seri(:, :) - rvc_ancien(:, :))/phys_tstep
       CALL water_int(klon, klev, q_seri, zmasse, zx_tmp_fi2d)
       d_q_dyn2d(:) = (zx_tmp_fi2d(:) - prw_ancien(:)) / phys_tstep
@@ -2461 +2472 @@
       IF (nqtot > nqo) d_tr_dyn(:, :, :) = (tr_seri(:, :, :) - tr_ancien(:, :, :)) / phys_tstep
       ! !! RomP <<<
-      !!d_rneb_dyn(:,:)=(rneb_seri(:,:)-rneb_ancien(:,:))/phys_tstep
-      d_rneb_dyn(:, :) = 0.0
     ELSE
       d_u_dyn(:, :) = 0.0
@@ -2470 +2479 @@
       d_ql_dyn(:, :) = 0.0
       d_qs_dyn(:, :) = 0.0
+      d_qbs_dyn(:, :) = 0.0
+      d_cf_dyn(:, :) = 0.0
+      d_rvc_dyn(:, :) = 0.0
       d_q_dyn2d(:) = 0.0
       d_ql_dyn2d(:) = 0.0
@@ -2477 +2489 @@
       IF (nqtot > nqo) d_tr_dyn(:, :, :) = 0.0
       ! !! RomP <<<
-      d_rneb_dyn(:, :) = 0.0
-      d_qbs_dyn(:, :) = 0.0
       ancien_ok = .TRUE.
     ENDIF
@@ -2585 +2595 @@
         ! "zmasse" changes a little.)
       ENDIF
+    ENDIF
+
+    !-- Needed for LSCP - condensation and ice supersaturation
+    IF (ok_ice_supersat) THEN
+      DO k = 1, klev
+        DO i = 1, klon
+          IF ((q_seri(i, k) + ql_seri(i, k) + qs_seri(i, k)) > 0.) THEN
+            ratio_qi_qtot(i, k) = qs_seri(i, k) / (q_seri(i, k) + ql_seri(i, k) + qs_seri(i, k))
+            rvc_seri(i, k) = rvc_seri(i, k) * q_seri(i, k) / (q_seri(i, k) + ql_seri(i, k) + qs_seri(i, k))
+          ELSE
+            ratio_qi_qtot(i, k) = 0.
+            rvc_seri(i, k) = 0.
+          ENDIF
+        ENDDO
+      ENDDO
     ENDIF
 
@@ -3738 +3763 @@
 
       !--mise à jour de flight_m et flight_h2o dans leur module
-      IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
-        CALL airplane(debut, pphis, pplay, paprs, t_seri)
-      ENDIF
+      !IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
+      !  CALL airplane(debut,pphis,pplay,paprs,t_seri)
+      !ENDIF
 
       CALL lscp(klon, klev, phys_tstep, missing_val, paprs, pplay, &
               t_seri, q_seri, qs_ini, ptconv, ratqs, &
-              d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, rneb_seri, &
+              d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, &
               pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb, &
               radocond, picefra, rain_lsc, snow_lsc, &
@@ -3750 +3775 @@
               prfl, psfl, rhcl, &
               zqasc, fraca, ztv, zpspsk, ztla, zthl, iflag_cld_th, &
-              iflag_ice_thermo, ok_ice_sursat, zqsatl, zqsats, distcltop, temp_cltop, &
-              pbl_tke(:, :, is_ave), pbl_eps(:, :, is_ave), qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
+              iflag_ice_thermo, distcltop, temp_cltop, &
+              pbl_tke(:, :, is_ave), pbl_eps(:, :, is_ave), &
+              cell_area, &
+              cf_seri, rvc_seri, u_seri, v_seri, &
+              qsub, qissr, qcld, subfra, issrfra, gamma_cond, ratio_qi_qtot, &
+              dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
+              dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
               Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
+              dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
               cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv, &
               qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
@@ -5421 +5452 @@
           d_qx(i, k, isol) = (qs_seri(i, k) - qx(i, k, isol)) / phys_tstep
         ENDIF
-        !--ice_sursat: nqo=4, on ajoute rneb
-        IF (nqo>=4 .AND. ok_ice_sursat) THEN
-          d_qx(i, k, irneb) = (rneb_seri(i, k) - qx(i, k, irneb)) / phys_tstep
+        !--ice_supersat: nqo=5, we add cloud fraction and cloudy water vapor to total water vapor ratio
+        IF (nqo>=5 .and. ok_ice_supersat) THEN
+          d_qx(i, k, icf) = (cf_seri(i, k) - qx(i, k, icf)) / phys_tstep
+          d_qx(i, k, irvc) = (rvc_seri(i, k) - qx(i, k, irvc)) / phys_tstep
         ENDIF
 
@@ -5458 +5490 @@
     qs_ancien(:, :) = qs_seri(:, :)
     qbs_ancien(:, :) = qbs_seri(:, :)
-    rneb_ancien(:, :) = rneb_seri(:, :)
+    cf_ancien(:, :) = cf_seri(:, :)
+    rvc_ancien(:, :) = rvc_seri(:, :)
     CALL water_int(klon, klev, q_ancien, zmasse, prw_ancien)
     CALL water_int(klon, klev, ql_ancien, zmasse, prlw_ancien)