Changeset 5452 for LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90

Timestamp:

Dec 23, 2024, 6:12:29 PM (40 hours ago)

Author:

aborella

Message:

First implementation of the contrails parameterisation
Lacks the emission of H2O + the impact on radiative transfer

File:

• LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90 (modified) (12 diffs)

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

@@ -327 +327 @@
        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, &
+       rcont_seri, d_rcont_dyn, flight_dist, flight_h2o, &
+       Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
+       dcf_avi, dqi_avi, dqvc_avi, &
        !
        cldemi,  &
@@ -511 +512 @@
     !
     ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional), blowing snow (optional)
-    INTEGER,SAVE :: ivap, iliq, isol, ibs, icf, irvc
-!$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, irvc)
+    INTEGER,SAVE :: ivap, iliq, isol, ibs, icf, irvc, ircont
+!$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, irvc, ircont)
     !
     !
@@ -1356 +1357 @@
        icf  = strIdx(tracers(:)%name, addPhase('H2O', 'f'))
        irvc = strIdx(tracers(:)%name, addPhase('H2O', 'c'))
+       ircont = strIdx(tracers(:)%name, addPhase('H2O', 'a'))
 !       CALL init_etat0_limit_unstruct
 !       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
@@ -1418 +1420 @@
        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)
+       ENDIF
+
+       IF (ok_plane_contrail.AND.(nqo.LT.6)) THEN
+          WRITE (lunout, *) ' ok_plane_contrail=y requires 6 H2O tracers ', &
+              '(H2O_g, H2O_l, H2O_s, H2O_f, H2O_c, H2O_a) but nqo=', nqo, '. Might as well stop here.'
           abort_message='see above'
           CALL abort_physic(modname,abort_message,1)
@@ -1843 +1852 @@
    &    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_supersat,iflag_ratqs,fl_cor_ebil, &
+       CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_supersat,&
+                     ok_plane_contrail,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, &
@@ -2433 +2443 @@
           cf_seri(i,k) = 0.
           rvc_seri(i,k)= 0.
+          rcont_seri(i,k)= 0.
           !CR: ATTENTION, on rajoute la variable glace
           IF (nqo.EQ.2) THEN             !--vapour and liquid only
@@ -2443 +2454 @@
                cf_seri(i,k) = qx(i,k,icf)
                rvc_seri(i,k) = qx(i,k,irvc)
+             ENDIF
+             IF (ok_plane_contrail) THEN
+               rcont_seri(i,k) = qx(i,k,ircont)
              ENDIF
              IF (ok_bs) THEN
@@ -2521 +2535 @@
        d_cf_dyn(:,:) = (cf_seri(:,:)-cf_ancien(:,:))/phys_tstep
        d_rvc_dyn(:,:)= (rvc_seri(:,:)-rvc_ancien(:,:))/phys_tstep
+       d_rcont_dyn(:,:)= (rcont_seri(:,:)-rcont_ancien(:,:))/phys_tstep
        CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
        d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
@@ -2542 +2557 @@
        d_cf_dyn(:,:) = 0.0
        d_rvc_dyn(:,:)= 0.0
+       d_rcont_dyn(:,:)= 0.0
        d_q_dyn2d(:)  = 0.0
        d_ql_dyn2d(:) = 0.0
@@ -3889 +3905 @@
          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, &
+         rcont_seri, flight_dist, flight_h2o, &
+         Tcritcont, qcritcont, potcontfraP, potcontfraNP, dcf_avi, dqi_avi, dqvc_avi, &
          cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
          qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
@@ -5606 +5622 @@
              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
+             IF (nqo.ge.6 .and. ok_plane_contrail) THEN
+               d_qx(i,k,ircont) = ( rcont_seri(i,k) - qx(i,k,ircont) ) / phys_tstep
+             ENDIF
           ENDIF
 
@@ -5642 +5661 @@
     cf_ancien(:,:) = cf_seri(:,:)
     rvc_ancien(:,:)= rvc_seri(:,:)
+    rcont_ancien(:,:)= rcont_seri(:,:)
     CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
     CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)