Changeset 4982 for LMDZ6/trunk/libf/phylmdiso

Timestamp:

Jun 15, 2024, 5:17:08 PM (6 months ago)

Author:

crisi

Message:

suppress isotope_params.def + update physiq_mod + proof of concept of 3rd dimension with reevap routine

Location:

LMDZ6/trunk/libf/phylmdiso

Files:

• add_wake_tend.F90 (modified) (4 diffs)
• isotopes_mod.F90 (modified) (11 diffs)
• isotopes_routines_mod.F90 (modified) (7 diffs)
• isotopes_verif_mod.F90 (modified) (1 diff)
• isotrac_routines_mod.F90 (modified) (1 diff)
• output_physiqex_mod.F90 (added)
• phyetat0_mod.F90 (modified) (11 diffs)
• phys_local_var_mod.F90 (modified) (32 diffs)
• physiq_mod.F90 (modified) (74 diffs)
• physiqex_mod.F90 (added)
• reevap.F90 (modified) (4 diffs)

LMDZ6/trunk/libf/phylmdiso/add_wake_tend.F90

@@ -1 @@
-SUBROUTINE add_wake_tend(zddeltat, zddeltaq, zds, zddensaw, zddensw, zoccur, text, abortphy &
+SUBROUTINE add_wake_tend(zddeltat, zddeltaq, zds, zdas, zddensw, zddensaw, zoccur, text, abortphy &
 #ifdef ISO
         , zddeltaxt &
@@ -13 +13 @@
 
 USE dimphy, ONLY: klon, klev
-USE phys_state_var_mod, ONLY: wake_deltat, wake_deltaq, wake_s,  &
-                              awake_dens, wake_dens
+USE phys_state_var_mod, ONLY: wake_deltat, wake_deltaq, wake_s, awake_s,  &
+                              wake_dens, awake_dens
 
 USE print_control_mod, ONLY: prt_level
@@ -26 +26 @@
 !------------
   REAL, DIMENSION(klon, klev),   INTENT (IN)         :: zddeltat, zddeltaq
-  REAL, DIMENSION(klon),         INTENT (IN)         :: zds, zddensaw, zddensw
+  REAL, DIMENSION(klon),         INTENT (IN)         :: zds, zdas, zddensw, zddensaw
   INTEGER, DIMENSION(klon),      INTENT (IN)         :: zoccur
   CHARACTER*(*),                 INTENT (IN)         :: text
@@ -79 +79 @@
            IF (zoccur(i) .GE. 1) THEN
              wake_s(i)     = wake_s(i)    + zds(i)
+             awake_s(i)    = awake_s(i)    + zdas(i)
+             wake_dens(i)  = wake_dens(i) + zddensw(i)
              awake_dens(i) = awake_dens(i) + zddensaw(i)
-             wake_dens(i)  = wake_dens(i) + zddensw(i)
            ELSE
              wake_s(i)     = 0.
+             awake_s(i)    = 0.
+             wake_dens(i)  = 0.
              awake_dens(i) = 0.
-             wake_dens(i)  = 0.
            ENDIF   ! (zoccur(i) .GE. 1)
          END DO

LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90

@@ -20 +20 @@
 
    !--- Variables not depending on isotopes
-   REAL,    SAVE :: pxtmelt, pxtice, pxtmin, pxtmax
-!$OMP THREADPRIVATE(pxtmelt, pxtice, pxtmin, pxtmax)
-   REAL,    SAVE :: tdifexp, tv0cin, thumxt1
-!$OMP THREADPRIVATE(tdifexp, tv0cin, thumxt1)
+   REAL,    SAVE :: thumxt1
+   !$OMP THREADPRIVATE(thumxt1)
    INTEGER, SAVE :: ntot 
 !$OMP THREADPRIVATE(ntot)
@@ -30 +28 @@
    REAL,    SAVE :: P_veg
 !$OMP THREADPRIVATE(P_veg)
-   REAL,    SAVE :: musi, lambda_sursat
-!$OMP THREADPRIVATE(musi, lambda_sursat)
-   REAL,    SAVE :: Kd
-!$OMP THREADPRIVATE(Kd)
-   REAL,    SAVE :: rh_cste_surf_cond, T_cste_surf_cond
-!$OMP THREADPRIVATE(rh_cste_surf_cond, T_cste_surf_cond)
+   REAL,    SAVE :: lambda_sursat
+!$OMP THREADPRIVATE(lambda_sursat)
    LOGICAL, SAVE :: bidouille_anti_divergence    ! T: regularly, xteau <- q to avoid slow drifts
 !$OMP THREADPRIVATE(bidouille_anti_divergence)
@@ -54 +48 @@
    REAL,    SAVE :: sstlatcrit, dsstlatcrit
 !$OMP THREADPRIVATE(sstlatcrit, dsstlatcrit)
-   REAL,    SAVE :: deltaO18_oce
-!$OMP THREADPRIVATE(deltaO18_oce)
    INTEGER, SAVE :: albedo_prescrit              ! 0: default ; 1: constant albedo
 !$OMP THREADPRIVATE(albedo_prescrit)
@@ -88 +80 @@
    REAL,    SAVE :: fac_modif_evaoce
 !$OMP THREADPRIVATE(fac_modif_evaoce)
+   REAL,    SAVE :: deltaO18_oce
+!$OMP THREADPRIVATE(deltaO18_oce)
    INTEGER, SAVE :: ok_bidouille_wake
 !$OMP THREADPRIVATE(ok_bidouille_wake)
@@ -106 +100 @@
                     alpha_liq_sol, Rdefault, Rmethox
 !$OMP THREADPRIVATE(alpha_liq_sol, Rdefault, Rmethox)
-   REAL, SAVE ::    fac_coeff_eq17_liq, fac_coeff_eq17_ice
-!$OMP THREADPRIVATE(fac_coeff_eq17_liq, fac_coeff_eq17_ice)
+!   REAL, SAVE ::    fac_coeff_eq17_liq, fac_coeff_eq17_ice
+!!$OMP THREADPRIVATE(fac_coeff_eq17_liq, fac_coeff_eq17_ice)
+
+   !--- H2[18]O reference
+   REAL, PARAMETER :: fac_enrichoce18=0.0005
+   REAL, PARAMETER :: alpha_liq_sol_O18=1.00291
+   REAL, PARAMETER :: talph1_O18=1137.
+   REAL, PARAMETER :: talph2_O18=-0.4156
+   REAL, PARAMETER :: talph3_O18=-2.0667E-3
+   REAL, PARAMETER :: talps1_O18=11.839
+   REAL, PARAMETER :: talps2_O18=-0.028244
+   REAL, PARAMETER :: tdifrel_O18=1./0.9723
+   REAL, PARAMETER :: tkcin0_O18=0.006
+   REAL, PARAMETER :: tkcin1_O18=0.000285
+   REAL, PARAMETER :: tkcin2_O18=0.00082
+   REAL, PARAMETER :: fac_coeff_eq17_liq=0.529
+   REAL, PARAMETER :: fac_coeff_eq17_ice=0.529
+
+   !---- Parameters that do not depend on the nature of water isotopes:
+   REAL, PARAMETER :: pxtmelt = 273.15 ! temperature at which ice formation starts
+   REAL, PARAMETER :: pxtice  = 273.15-10.0 ! -- temperature at which all condensate is ice:
+   REAL, PARAMETER :: pxtmin = 273.15 - 120.0   ! On ne calcule qu'au dessus de -120°C
+   REAL, PARAMETER :: pxtmax = 273.15 +  60.0   ! On ne calcule qu'au dessus de +60°C
+   REAL, PARAMETER :: tdifexp = 0.58 ! -- a constant for alpha_eff for equilibrium below cloud base:
+   REAL, PARAMETER :: tv0cin  = 7.0 ! wind threshold (m/s) for smooth/rough regime (Merlivat and Jouzel 1979)
+   REAL, PARAMETER :: musi=1.0  ! facteurs lambda et mu dans Si=musi-lambda*T
+   REAL, PARAMETER :: Kd=2.5e-9 ! m2/s ! diffusion dans le sol
+   REAL, PARAMETER :: rh_cste_surf_cond = 0.6 ! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir
+   REAL, PARAMETER :: T_cste_surf_cond = 288.0
+
 
    !--- Negligible lower thresholds: no need to check for absurd values under these lower limits
@@ -140 +162 @@
    INTEGER :: ixt
 
-   !--- H2[18]O reference
-   REAL :: fac_enrichoce18, alpha_liq_sol_O18, &
-           talph1_O18, talph2_O18, talph3_O18, talps1_O18, talps2_O18, &
-           tkcin0_O18, tkcin1_O18, tkcin2_O18, tdifrel_O18  
-
+ 
    !--- For H2[17]O
    REAL    :: fac_kcin, pente_MWL
@@ -152 +170 @@
    LOGICAL, PARAMETER ::   ok_nocinsat = .FALSE. ! if T: no sursaturation effect for ice
    LOGICAL, PARAMETER :: Rdefault_smow = .FALSE. ! if T: Rdefault=smow; if F: nul
+   LOGICAL, PARAMETER :: tnat1 = .TRUE. ! If T: all tnats are 1.
 
    !--- For [3]H
@@ -157 +176 @@
 
    CHARACTER(LEN=maxlen) :: modname, sxt
-   REAL, ALLOCATABLE :: tmp(:)
 
    modname = 'iso_init'
@@ -214 +232 @@
       CALL get_in('lat_max_albedo', lat_max_albedo,  100.)
    END IF
-   deltaO18_oce=0.0
    CALL get_in('deltaP_BL',           deltaP_BL,     10.0)
    CALL get_in('ruissellement_pluie', ruissellement_pluie, 0)
@@ -249 +266 @@
    CALL get_in('ok_prod_nucl_tritium', ok_prod_nucl_tritium, .FALSE., .FALSE.)
 
-   !--------------------------------------------------------------
-   ! Parameters that do not depend on the nature of water isotopes:
-   !--------------------------------------------------------------
-   ! -- temperature at which ice condensate starts to form (valeur ECHAM?):
-   pxtmelt = 273.15
-
-   ! -- temperature at which all condensate is ice:
-   pxtice  = 273.15-10.0
-
-   !- -- test PHASE
-!   pxtmelt = 273.15 - 10.0
-!   pxtice  = 273.15 - 30.0
-
-   ! -- minimum temperature to calculate fractionation coeff
-   pxtmin = 273.15 - 120.0   ! On ne calcule qu'au dessus de -120°C
-   pxtmax = 273.15 +  60.0   ! On ne calcule qu'au dessus de +60°C
-   !    Remarque: les coeffs ont ete mesures seulement jusq'à -40!
-
-   ! -- a constant for alpha_eff for equilibrium below cloud base:
-   tdifexp = 0.58
-   tv0cin  = 7.0
-
-   ! facteurs lambda et mu dans Si=musi-lambda*T
-   musi=1.0
-   if (ok_nocinsat) lambda_sursat = 0.0          ! no sursaturation effect
-
-   ! diffusion dans le sol
-   Kd=2.5e-9 ! m2/s   
-
-   ! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir
-   rh_cste_surf_cond = 0.6
-    T_cste_surf_cond = 288.0
+   ! Ocean composition
+   CALL get_in('deltaO18_oce',  deltaO18_oce, 0.0)
    
    CALL msg('iso_O18, iso_HDO, iso_eau = '//TRIM(strStack(int2str([iso_O18, iso_HDO, iso_eau]))), modname)
 
    !--------------------------------------------------------------
-   ! Parameters that depend on the nature of water isotopes:
+   ! Isotope fractionation factors and a few isotopic constants
    !--------------------------------------------------------------
-   IF(getKey('tnat',    tnat,    isoName)) CALL abort_physic(modname, 'can''t get tnat',    1)
-   IF(getKey('toce',    toce,    isoName)) CALL abort_physic(modname, 'can''t get toce',    1)
-   IF(getKey('tcorr',   tcorr,   isoName)) CALL abort_physic(modname, 'can''t get tcorr',   1)
-   IF(getKey('talph1',  talph1,  isoName)) CALL abort_physic(modname, 'can''t get talph1',  1)
-   IF(getKey('talph2',  talph2,  isoName)) CALL abort_physic(modname, 'can''t get talph2',  1)
-   IF(getKey('talph3',  talph3,  isoName)) CALL abort_physic(modname, 'can''t get talph3',  1)
-   IF(getKey('talps1',  talps1,  isoName)) CALL abort_physic(modname, 'can''t get talps1',  1)
-   IF(getKey('talps2',  talps2,  isoName)) CALL abort_physic(modname, 'can''t get talps2',  1)
-   IF(getKey('tkcin0',  tkcin0,  isoName)) CALL abort_physic(modname, 'can''t get tkcin0',  1)
-   IF(getKey('tkcin1',  tkcin1,  isoName)) CALL abort_physic(modname, 'can''t get tkcin1',  1)
-   IF(getKey('tkcin2',  tkcin2,  isoName)) CALL abort_physic(modname, 'can''t get tkcin2',  1)
-   IF(getKey('tdifrel', tdifrel, isoName)) CALL abort_physic(modname, 'can''t get tdifrel', 1)
-   IF(getKey('alpha_liq_sol', alpha_liq_sol, isoName)) CALL abort_physic(modname, 'can''t get alpha_liq_sol',  1)
-   IF(getKey('Rdefault',Rdefault,isoName)) CALL abort_physic(modname, 'can''t get Rdefault',1)
-   IF(getKey('Rmethox', Rmethox, isoName)) CALL abort_physic(modname, 'can''t get Rmethox', 1)
+   ALLOCATE(tkcin0(niso))
+   ALLOCATE(tkcin1(niso))
+   ALLOCATE(tkcin2(niso))
+   ALLOCATE(tnat(niso))
+   ALLOCATE(tdifrel(niso))
+   ALLOCATE(toce(niso))
+   ALLOCATE(tcorr(niso))
+   ALLOCATE(talph1(niso))
+   ALLOCATE(talph2(niso))
+   ALLOCATE(talph3(niso))
+   ALLOCATE(talps1(niso))
+   ALLOCATE(talps2(niso))
+   ALLOCATE(alpha_liq_sol(niso))
+   ALLOCATE(Rdefault(niso))
+   ALLOCATE(Rmethox(niso))
+
+   do ixt=1,niso
+     if (ixt.eq.iso_HTO) then  ! Tritium
+       tkcin0(ixt) = 0.01056
+       tkcin1(ixt) = 0.0005016
+       tkcin2(ixt) = 0.0014432
+       if (tnat1) then 
+               tnat(ixt)=1 
+       else 
+               tnat(ixt)=0.
+       endif
+       toce(ixt)=4.0E-19 ! rapport T/H = 0.2 TU Dreisigacker and Roether 1978
+       tcorr(ixt)=1.
+       tdifrel(ixt)=1./0.968
+       talph1(ixt)=46480.
+       talph2(ixt)=-103.87
+       talph3(ixt)=0.
+       talps1(ixt)=46480.
+       talps2(ixt)=-103.87
+       alpha_liq_sol(ixt)=1. 
+       Rmethox(ixt)=0.0
+     endif
+     if (ixt.eq.iso_O17) then  ! O17
+       pente_MWL=0.528 
+       tdifrel(ixt)=1./0.98555 ! valeur utilisée en 1D et dans modèle de LdG ! tdifrel(ixt)=1./0.985452 ! donné par Amaelle
+       fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0) ! fac_kcin=0.5145 ! donné par Amaelle
+       tkcin0(ixt) = tkcin0_O18*fac_kcin
+       tkcin1(ixt) = tkcin1_O18*fac_kcin
+       tkcin2(ixt) = tkcin2_O18*fac_kcin
+       if (tnat1) then 
+               tnat(ixt)=1 
+       else 
+               tnat(ixt)=0.004/100. ! O17 représente 0.004% de l'oxygène
+       endif
+       toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)**pente_MWL
+       tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL ! donné par Amaelle            
+       talph1(ixt)=talph1_O18
+       talph2(ixt)=talph2_O18
+       talph3(ixt)=talph3_O18
+       talps1(ixt)=talps1_O18
+       talps2(ixt)=talps2_O18     
+       alpha_liq_sol(ixt)=(alpha_liq_sol_O18)**fac_coeff_eq17_liq 
+       Rdefault(ixt)=tnat(ixt)*(-3.15/1000.0+1.0)
+       Rmethox(ixt)=(230./1000.+1.)*tnat(ixt) !Zahn et al 2006
+     endif
+     if (ixt.eq.iso_O18) then  ! Oxygene18
+       tkcin0(ixt) = tkcin0_O18
+       tkcin1(ixt) = tkcin1_O18
+       tkcin2(ixt) = tkcin2_O18
+       if (tnat1) then
+               tnat(ixt)=1 
+       else 
+               tnat(ixt)=2005.2E-6
+       endif
+       toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)
+       tcorr(ixt)=1.0+fac_enrichoce18
+       tdifrel(ixt)=tdifrel_O18
+       talph1(ixt)=talph1_O18
+       talph2(ixt)=talph2_O18
+       talph3(ixt)=talph3_O18
+       talps1(ixt)=talps1_O18
+       talps2(ixt)=talps2_O18
+       alpha_liq_sol(ixt)=alpha_liq_sol_O18   
+       Rdefault(ixt)=tnat(ixt)*(-6.0/1000.0+1.0)
+       Rmethox(ixt)=(130./1000.+1.)*tnat(ixt) !Zahn et al 2006  
+     endif
+     if (ixt.eq.iso_HDO) then ! Deuterium
+       pente_MWL=8.0
+       tdifrel(ixt)=1./0.9755 !          fac_kcin=0.88
+       fac_kcin= (tdifrel(ixt)-1)/(tdifrel_O18-1)
+       tkcin0(ixt) = tkcin0_O18*fac_kcin
+       tkcin1(ixt) = tkcin1_O18*fac_kcin
+       tkcin2(ixt) = tkcin2_O18*fac_kcin
+       if (tnat1) then 
+               tnat(ixt)=1
+       else
+               tnat(ixt)=155.76E-6
+       endif
+       toce(ixt)=tnat(ixt)*(1.0+pente_MWL*deltaO18_oce/1000.0)
+       tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL          
+       talph1(ixt)=24844.
+       talph2(ixt)=-76.248
+       talph3(ixt)=52.612E-3
+       talps1(ixt)=16288.
+       talps2(ixt)=-0.0934
+       !ZXalpha_liq_sol=1.0192 ! Weston, Ralph, 1955
+       alpha_liq_sol(ixt)=1.0212 
+       ! valeur de Lehmann & Siegenthaler, 1991, Journal of
+       ! Glaciology, vol 37, p 23
+       Rdefault(ixt)=tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0)
+       Rmethox(ixt)=tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006
+     endif
+     if (ixt.eq.iso_eau) then ! Oxygene16
+       tkcin0(ixt) = 0.0
+       tkcin1(ixt) = 0.0
+       tkcin2(ixt) = 0.0
+       tnat(ixt)=1.
+       toce(ixt)=tnat(ixt)
+       tcorr(ixt)=1.0
+       tdifrel(ixt)=1.
+       talph1(ixt)=0.
+       talph2(ixt)=0.
+       talph3(ixt)=0.
+       talps1(ixt)=0.
+       talph3(ixt)=0.
+       alpha_liq_sol(ixt)=1.
+       Rdefault(ixt)=tnat(ixt)*1.0
+       Rmethox(ixt)=1.0
+     endif
+   enddo ! ixt=1,niso
 
    IF(.NOT.Rdefault_smow) then
@@ -308 +405 @@
    ENDIF
    write(*,*) 'Rdefault=',Rdefault
+   write(*,*) 'toce=',toce
 
    !--- Sensitivity test: no kinetic effect in sfc evaporation

LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90

@@ -1525 +1525 @@
 #endif       
        pxtfra=max(min(pxtfra,alpha_max),0.0)
-
 
       end subroutine fractcalk_liq
@@ -15922 +15921 @@
 
       ! verif
-!      text="phyisoetat0 67"
-!      write(*,*) 'snow(8,1)=',snow(8,1)
-!      write(*,*) 'xtsnow(4,8,1)=',xtsnow(4,8,1)
 #ifdef ISOVERIF
       do i=1,klon
@@ -15934 +15930 @@
          enddo !do ixt=1,niso
       enddo !do i=1,klon
-#endif      
-#ifdef ISOVERIF
       do i=1,klon
          if (iso_eau.gt.0) then
@@ -16021 +16015 @@
          endif 
        enddo !do i=1,klon
-
 #endif
       !end verif
@@ -16128 +16121 @@
           deltaD_run_off_lic_0(ixt)=deltaD_sol(ixt)
           deltaD_land_ice(ixt)=deltaD_snow(ixt)
-          call fractcalk_liq(ixt, 283.0, alpha(ixt))            
+          call fractcalk_liq(ixt, 283.0, alpha(ixt))    
         enddo !do ixt=1,niso
         call calcul_kcin(2.0,kcin) 
@@ -18830 +18823 @@
         if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO)) then
             if (q.gt.ridicule) then
+                    write(*,*) 'xt,q=',xt,q
+                    write(*,*) 'alpha=',alpha
+                    write(*,*) 'toce,kcin,h0=',toce,kcin,h0
+                    write(*,*) 'RMerlivat=',RMerlivat
                 call iso_verif_aberrant_encadre( xt/q, 'isotopes_routines_mod 18930b: iso_init_ideal')
             endif
@@ -18902 +18899 @@
 end subroutine appel_stewart_debug
 
+
+subroutine dispatch(klon,klev,qx,q_seri,xt_seri,ql_seri,xtl_seri,qs_seri,xts_seri)
+
+use infotrac_phy, ONLY: nqtot,nqo,ivap,iliq,isol,iqIsoPha,ntraciso=>ntiso
+implicit none
+
+! inputs
+integer, intent(in) :: klon,klev
+real,dimension(klon,klev,nqtot), intent(in) ::qx
+
+! outputs
+real,dimension(klon,klev), intent(out) ::q_seri,ql_seri,qs_seri
+real,dimension(ntraciso,klon,klev), intent(out) :: xt_seri,xtl_seri,xts_seri
+
+! locals
+integer :: i,k,ixt
+
+do k=1,klev
+do i=1,klon
+    q_seri(i,k)  = qx(i,k,ivap)
+    ql_seri(i,k) = qx(i,k,iliq)
+    IF (nqo.EQ.2) THEN             !--vapour and liquid only
+             qs_seri(i,k) = 0.
+    ELSE IF (nqo.ge.3) THEN        !--vapour, liquid and ice
+             qs_seri(i,k) = qx(i,k,isol)
+    ENDIF
+    do ixt=1,ntraciso
+          xt_seri(ixt,i,k)  = qx(i,k,iqIsoPha(ixt,ivap))
+          xtl_seri(ixt,i,k) = qx(i,k,iqIsoPha(ixt,iliq))
+          if (nqo.eq.2) then
+             xts_seri(ixt,i,k) = 0.
+          else if (nqo.eq.3) then
+             xts_seri(ixt,i,k) = qx(i,k,iqIsoPha(ixt,isol))
+          endif
+    enddo !do ixt=1,niso
+
+enddo
+enddo
+
+end subroutine dispatch
+
+subroutine together(klon,klev,qx,q_seri,xt_seri,ql_seri,xtl_seri,qs_seri,xts_seri)
+
+use infotrac_phy, ONLY: nqtot,nqo,ivap,iliq,isol,iqIsoPha,ntraciso=>ntiso
+implicit none
+
+! inputs
+integer, intent(in) :: klon,klev
+real,dimension(klon,klev), intent(in) ::q_seri,ql_seri,qs_seri
+real,dimension(ntraciso,klon,klev), intent(in) :: xt_seri,xtl_seri,xts_seri
+
+! inputs
+real,dimension(klon,klev,nqtot), intent(out) ::qx
+
+! locals
+integer :: i,k,ixt
+
+do k=1,klev
+do i=1,klon  
+    qx(i,k,ivap)  = q_seri(i,k)
+    qx(i,k,iliq) = ql_seri(i,k)
+    IF (nqo.ge.3) THEN        !--vapour, liquid and ice
+             qx(i,k,isol) = qs_seri(i,k)
+    ENDIF
+    do ixt=1,ntraciso
+          qx(i,k,iqIsoPha(ixt,ivap)) = xt_seri(ixt,i,k)  
+          qx(i,k,iqIsoPha(ixt,iliq)) = xtl_seri(ixt,i,k)
+          if (nqo.ge.3) then
+             qx(i,k,iqIsoPha(ixt,isol)) = xts_seri(ixt,i,k)
+          endif
+    enddo !do ixt=1,niso
+
+enddo
+enddo
+
+end subroutine together
+
+
 END MODULE isotopes_routines_mod
 #endif

LMDZ6/trunk/libf/phylmdiso/isotopes_verif_mod.F90

@@ -1042 +1042 @@
             write(*,*) 'deltaD=',deltaD
             write(*,*) 'Dexcess=',dexcess
+            write(*,*) 'tnat=',tnat
 !            stop
             iso_verif_O18_aberrant_nostop=1

LMDZ6/trunk/libf/phylmdiso/isotrac_routines_mod.F90

@@ -681 +681 @@
               Eqi_prime_cas(il)=Eqi_prime(cas(il)) &
      &           *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
-              Eqi_cas(il)=Eqi(il) &
+              Eqi_cas(il)=Eqi(cas(il)) & ! corr bug Camille 15 juin 2024
      &           *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
             else

LMDZ6/trunk/libf/phylmdiso/phyetat0_mod.F90

@@ -29 +29 @@
        solsw, solswfdiff, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
        wake_deltat, wake_delta_pbl_TKE, delta_tsurf, beta_aridity, wake_fip, wake_pe, &
-       wake_s, wake_dens, awake_dens, cv_gen, zgam, zmax0, zmea, zpic, zsig, &
+       wake_s, awake_s, wake_dens, awake_dens, cv_gen, zgam, zmax0, zmea, zpic, zsig, &
 #ifdef ISO
        fxtevap, xtsol, xt_ancien, xtl_ancien, xts_ancien, wake_deltaxt, &
@@ -48 +48 @@
   USE time_phylmdz_mod, ONLY: init_iteration, pdtphys, itau_phy
   USE wxios, ONLY: missing_val_xios => missing_val, using_xios
-  use netcdf, only:  nf90_fill_real
+  use netcdf, only: missing_val_netcdf => nf90_fill_real
   use config_ocean_skin_m, only: activate_ocean_skin
 #ifdef ISO
@@ -111 +111 @@
   REAL Rland_ice(niso,klon)
 #endif
+
+  IF (using_xios) THEN
+    missing_val=missing_val_xios
+  ELSE
+    missing_val=missing_val_netcdf
+  ENDIF
+
   ! FH1D
   !     real iolat(jjm+1)
@@ -116 +123 @@
 
   ! Ouvrir le fichier contenant l'etat initial:
-  IF (using_xios) THEN
-    missing_val = missing_val_xios
-  ELSE
-    missing_val =  nf90_fill_real
-  ENDIF
 
   CALL open_startphy(fichnom)
@@ -323 +325 @@
 
 !===================================================================
+! Lecture dans le cas iflag_pbl_surface =1
+!===================================================================
+
+   if ( iflag_physiq <= 1 ) then
+!===================================================================
   ! Lecture des temperatures du sol profond:
 !===================================================================
@@ -350 +357 @@
   found=phyetat0_get(snow_fall,"snow_f","snow fall",0.)
   found=phyetat0_get(rain_fall,"rain_f","rain fall",0.)
-
   IF (ok_bs) THEN
      found=phyetat0_get(bs_fall,"bs_f","blowing snow fall",0.)
@@ -404 +410 @@
   ENDIF
 
+  endif ! iflag_physiq <= 1
+
   ! Lecture de l'age de la neige:
   found=phyetat0_srf(agesno,"AGESNO","SNOW AGE",0.001)
@@ -427 +435 @@
      prbsw_ancien(:)=0.
   ENDIF
-
 
   ! Ehouarn: addtional tests to check if t_ancien, q_ancien contain
@@ -450 +457 @@
   ENDIF
 
-
   found=phyetat0_get(clwcon,"CLWCON","CLWCON",0.)
   found=phyetat0_get(rnebcon,"RNEBCON","RNEBCON",0.)
@@ -485 +491 @@
   found=phyetat0_get(wake_deltaq,"WAKE_DELTAQ","Delta hum. wake/env",0.)
   found=phyetat0_get(wake_s,"WAKE_S","Wake frac. area",0.)
+  found=phyetat0_get(awake_s,"AWAKE_S","Active Wake frac. area",0.)
 !jyg<
 !  Set wake_dens to -1000. when there is no restart so that the actual
@@ -664 +671 @@
 !        write(*,*) 'xtsnow(:,994,2)=',xtsnow(:,994,2)
 !#endif
-
+  if ( iflag_physiq <= 1 ) then
   CALL pbl_surface_init(fder, snow, qsurf, tsoil)
 #ifdef ISO
   CALL pbl_surface_init_iso(xtsnow,Rland_ice)
 #endif
+  endif
 
   ! Initialize module ocean_cpl_mod for the case of coupled ocean

LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90

@@ -14 +14 @@
       REAL, SAVE, ALLOCATABLE :: ql_seri(:,:),qs_seri(:,:)
       !$OMP THREADPRIVATE(ql_seri,qs_seri)
+      REAL, SAVE, ALLOCATABLE :: qx_seri(:,:,:)
+      !$OMP THREADPRIVATE(qx_seri)
       REAL, SAVE, ALLOCATABLE :: qbs_seri(:,:)
       !$OMP THREADPRIVATE(qbs_seri)
@@ -64 +66 @@
       REAL, SAVE, ALLOCATABLE :: d_t_eva(:,:),d_q_eva(:,:),d_ql_eva(:,:),d_qi_eva(:,:)
       !$OMP THREADPRIVATE(d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
+      REAL, SAVE, ALLOCATABLE :: d_qx_eva(:,:,:)
+      !$OMP THREADPRIVATE(d_qx_eva)
       REAL, SAVE, ALLOCATABLE :: d_t_lscst(:,:),d_q_lscst(:,:)
       !$OMP THREADPRIVATE(d_t_lscst,d_q_lscst)
@@ -124 +128 @@
       REAL, SAVE, ALLOCATABLE :: xts_seri(:,:,:)
       !$OMP THREADPRIVATE( xts_seri)
+      REAL, SAVE, ALLOCATABLE :: xtbs_seri(:,:,:)
+      !$OMP THREADPRIVATE( xtbs_seri)
       REAL, SAVE, ALLOCATABLE :: d_xt_eva(:,:,:)
       !$OMP THREADPRIVATE( d_xt_eva)
@@ -134 +140 @@
       REAL, SAVE, ALLOCATABLE :: d_xt_dyn(:,:,:)
       !$OMP THREADPRIVATE( d_xt_dyn)
-      REAL, SAVE, ALLOCATABLE :: d_xtl_dyn(:,:,:), d_xts_dyn(:,:,:)
-      !$OMP THREADPRIVATE(d_xtl_dyn, d_xts_dyn)
+      REAL, SAVE, ALLOCATABLE :: d_xtl_dyn(:,:,:), d_xts_dyn(:,:,:), d_xtbs_dyn(:,:,:)
+      !$OMP THREADPRIVATE(d_xtl_dyn, d_xts_dyn, d_xtbs_dyn)
       REAL, SAVE, ALLOCATABLE :: d_xt_con(:,:,:)
       !$OMP THREADPRIVATE( d_xt_con)
@@ -292 +298 @@
 !$OMP THREADPRIVATE(toplwad0_aerop, sollwad0_aerop)
 
+!AI 08 2023 ajout pour Ecrad
+      REAL,ALLOCATABLE,SAVE :: topswad_aero_s2(:), solswad_aero_s2(:)
+!$OMP THREADPRIVATE(topswad_aero_s2, solswad_aero_s2)
+      REAL,ALLOCATABLE,SAVE :: topswai_aero_s2(:), solswai_aero_s2(:)
+!$OMP THREADPRIVATE(topswai_aero_s2, solswai_aero_s2)
+      REAL,ALLOCATABLE,SAVE :: topswad0_aero_s2(:), solswad0_aero_s2(:)
+!$OMP THREADPRIVATE(topswad0_aero_s2, solswad0_aero_s2)
+      REAL,ALLOCATABLE,SAVE :: topsw_aero_s2(:,:), topsw0_aero_s2(:,:)
+!$OMP THREADPRIVATE(topsw_aero_s2, topsw0_aero_s2)
+      REAL,ALLOCATABLE,SAVE :: solsw_aero_s2(:,:), solsw0_aero_s2(:,:)
+!$OMP THREADPRIVATE(solsw_aero_s2, solsw0_aero_s2)
+      REAL,ALLOCATABLE,SAVE :: topswcf_aero_s2(:,:), solswcf_aero_s2(:,:)
+!$OMP THREADPRIVATE(topswcf_aero_s2, solswcf_aero_s2)
+! additional LW variables CK
+      REAL,ALLOCATABLE,SAVE :: toplwad_aero_s2(:), sollwad_aero_s2(:)
+!$OMP THREADPRIVATE(toplwad_aero_s2, sollwad_aero_s2)
+      REAL,ALLOCATABLE,SAVE :: toplwai_aero_s2(:), sollwai_aero_s2(:)
+!$OMP THREADPRIVATE(toplwai_aero_s2, sollwai_aero_s2)
+      REAL,ALLOCATABLE,SAVE :: toplwad0_aero_s2(:), sollwad0_aero_s2(:)
+!$OMP THREADPRIVATE(toplwad0_aero_s2, sollwad0_aero_s2)
+
 !Ajout de celles n??cessaires au phys_output_write_mod
       REAL, SAVE, ALLOCATABLE :: tal1(:), pal1(:), pab1(:), pab2(:)
@@ -300 +327 @@
 !$OMP THREADPRIVATE(sens, flwp, fiwp)
 !!
-!FC
+!FC 
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: zxfluxt, zxfluxq
 !$OMP THREADPRIVATE(zxfluxt, zxfluxq)
@@ -348 +375 @@
 !$OMP THREADPRIVATE(cdragm, cdragh)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: cldh, cldl, cldm, cldq, cldt, qsat2m
-!$OMP THREADPRIVATE(cldh, cldl, cldm, cldq, cldt, qsat2m )
+!$OMP THREADPRIVATE(cldh, cldl, cldm, cldq, cldt, qsat2m)
 !AS: cldhjn, cldljn, cldmjn,cldtjn pas utilisés en tant que variables, juste noms de diagnostics
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: JrNt
@@ -573 +600 @@
 !$OMP THREADPRIVATE(phiwriteSTD, qwriteSTD, twriteSTD, rhwriteSTD)
 
+
 ! ug et d'autres encore:
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: beta_prec
@@ -580 +608 @@
       REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: pfraclr,pfracld
 !$OMP THREADPRIVATE(pfraclr,pfracld)
-      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: distcltop
-!$OMP THREADPRIVATE(distcltop)
-      REAL, SAVE, ALLOCATABLE :: temp_cltop(:,:) 
-!$OMP THREADPRIVATE(temp_cltop)
 
 ! variables de sorties MM
@@ -645 +669 @@
       REAL, SAVE, ALLOCATABLE :: fcontrP(:,:)
       !$OMP THREADPRIVATE(fcontrP)
+      REAL, SAVE, ALLOCATABLE :: distcltop(:,:)
+      !$OMP THREADPRIVATE(distcltop)
+      REAL, SAVE, ALLOCATABLE :: temp_cltop(:,:)
+      !$OMP THREADPRIVATE(temp_cltop)
+
 
 !--POPRECIP variables
@@ -675 +704 @@
 
 
-
+      
 
 
@@ -681 +710 @@
 !
 ! variables for stratospheric aerosol
+      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: d_q_emiss
+!$OMP THREADPRIVATE(d_q_emiss)
       REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: R2SO4
 !$OMP THREADPRIVATE(R2SO4)
@@ -695 +726 @@
       REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: SO2_lifetime
 !$OMP THREADPRIVATE(SO2_lifetime)
+      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: H2SO4_lifetime
+!$OMP THREADPRIVATE(H2SO4_lifetime)
+      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: O3_clim
+!$OMP THREADPRIVATE(O3_clim)
       REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: alpha_bin
 !$OMP THREADPRIVATE(alpha_bin)
@@ -759 +794 @@
 SUBROUTINE phys_local_var_init
 USE dimphy
-USE infotrac_phy, ONLY : nbtr
+USE infotrac_phy, ONLY : nbtr,nqtot
 #ifdef ISO
 USE infotrac_phy, ONLY : ntraciso=>ntiso,niso
@@ -770 +805 @@
 IMPLICIT NONE
       ALLOCATE(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev), qbs_seri(klon,klev))
+      ALLOCATE(qx_seri(klon,klev,nqtot))
       ALLOCATE(u_seri(klon,klev),v_seri(klon,klev))
       ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
+      ALLOCATE(pbl_eps(klon,klev+1,nbsrf+1))
+      pbl_eps(:,:,:)=0.
       l_mix(:,:,:)=0.;l_mixmin(:,:,:)=0.;wprime(:,:,:)=0. ! doit etre initialse car pas toujours remplis
-      ALLOCATE(pbl_eps(klon,klev+1,nbsrf+1))      
+      ALLOCATE(rhcl(klon,klev))
       ALLOCATE(tr_seri(klon,klev,nbtr))
       ALLOCATE(d_t_dyn(klon,klev),d_q_dyn(klon,klev))
@@ -795 +833 @@
       ALLOCATE(d_u_ajs(klon,klev),d_v_ajs(klon,klev))
       ALLOCATE(d_t_eva(klon,klev),d_q_eva(klon,klev))
+      ALLOCATE(d_qx_eva(klon,klev,nqtot))
       ALLOCATE(d_ql_eva(klon,klev),d_qi_eva(klon,klev))
       ALLOCATE(d_t_lscst(klon,klev),d_q_lscst(klon,klev))
@@ -808 +847 @@
       allocate(xtl_seri(ntraciso,klon,klev))
       allocate(xts_seri(ntraciso,klon,klev))
+      allocate(xtbs_seri(ntraciso,klon,klev))
       allocate(d_xt_dyn(ntraciso,klon,klev))
       allocate(d_xtl_dyn(ntraciso,klon,klev))
       allocate(d_xts_dyn(ntraciso,klon,klev))
+      allocate(d_xtbs_dyn(ntraciso,klon,klev))
       allocate(d_xt_con(ntraciso,klon,klev))
       allocate(d_xt_wake(ntraciso,klon,klev))
@@ -841 +882 @@
       ALLOCATE(d_u_lif(klon,klev),d_v_lif(klon,klev))
       ALLOCATE(d_ts(klon,nbsrf), d_tr(klon,klev,nbtr))
+
+! aerosols
+      ALLOCATE(m_allaer(klon,klev,naero_tot))
 ! Special RRTM
       ALLOCATE(ZLWFT0_i(klon,klev+1),ZSWFT0_i(klon,klev+1),ZFLDN0(klon,klev+1))
@@ -919 +963 @@
       ALLOCATE(toplwad0_aerop(klon), sollwad0_aerop(klon))
 
+!AI Ajout Ecrad (3Deffect)
+      ALLOCATE(topswad_aero_s2(klon), solswad_aero_s2(klon))
+      ALLOCATE(topswai_aero_s2(klon), solswai_aero_s2(klon))
+      ALLOCATE(topswad0_aero_s2(klon), solswad0_aero_s2(klon))
+      ALLOCATE(topsw_aero_s2(klon,naero_grp), topsw0_aero_s2(klon,naero_grp))
+      ALLOCATE(solsw_aero_s2(klon,naero_grp), solsw0_aero_s2(klon,naero_grp))
+      ALLOCATE(topswcf_aero_s2(klon,naero_grp), solswcf_aero_s2(klon,naero_grp))
+! additional LW variables CK
+      ALLOCATE(toplwad_aero_s2(klon), sollwad_aero_s2(klon))
+      ALLOCATE(toplwai_aero_s2(klon), sollwai_aero_s2(klon))
+      ALLOCATE(toplwad0_aero_s2(klon), sollwad0_aero_s2(klon))
+
 ! FH Ajout de celles necessaires au phys_output_write_mod
 
@@ -1038 +1094 @@
       ALLOCATE(wfevap(klon, nbsrf))
       ALLOCATE(evap_pot(klon, nbsrf))
-! FC
+! FC 
       ALLOCATE(zxfluxq(klon,klev),zxfluxt(klon,klev))
 !
@@ -1045 +1101 @@
       ALLOCATE(pmflxr(klon, klev+1), pmflxs(klon, klev+1))
       ALLOCATE(wdtrainA(klon,klev),wdtrainS(klon,klev),wdtrainM(klon,klev))
-      ALLOCATE(dnwd(klon, klev), upwd(klon, klev) )
+      ALLOCATE(dnwd(klon, klev), upwd(klon, klev))
       ALLOCATE(ep(klon,klev))                          ! epmax_cape
-      ALLOCATE(da(klon,klev), mp(klon,klev) )
-      ALLOCATE(phi(klon,klev,klev) )
-      ALLOCATE(wght_cvfd(klon,klev) )
-      ALLOCATE(phi2(klon,klev,klev) )
+      ALLOCATE(da(klon,klev), mp(klon,klev))
+      ALLOCATE(phi(klon,klev,klev))
+      ALLOCATE(wght_cvfd(klon,klev))
+      ALLOCATE(phi2(klon,klev,klev))
       ALLOCATE(d1a(klon,klev), dam(klon,klev))
-      ALLOCATE(ev(klon,klev) )
-      ALLOCATE(elij(klon,klev,klev) )
-      ALLOCATE(qtaa(klon,klev) )
-      ALLOCATE(clw(klon,klev) )
-      ALLOCATE(epmlmMm(klon,klev,klev), eplaMm(klon,klev) )
-      ALLOCATE(sij(klon,klev,klev) )
+      ALLOCATE(ev(klon,klev))
+      ALLOCATE(elij(klon,klev,klev))
+      ALLOCATE(qtaa(klon,klev))
+      ALLOCATE(clw(klon,klev))
+      ALLOCATE(epmlmMm(klon,klev,klev), eplaMm(klon,klev))
+      ALLOCATE(sij(klon,klev,klev))
 #ifdef ISO
       ALLOCATE(xtwdtrainA(ntraciso,klon,klev))
@@ -1128 +1184 @@
 
 !--POPRECIP variables
-      ALLOCATE(dqreva(klon,klev),dqssub(klon,klev))
-      ALLOCATE(qraindiag(klon,klev), qsnowdiag(klon,klev)) 
+      ALLOCATE(qraindiag(klon,klev), qsnowdiag(klon,klev))
+      ALLOCATE(dqreva(klon,klev), dqssub(klon,klev))
       ALLOCATE(dqrauto(klon,klev), dqrcol(klon,klev), dqrmelt(klon,klev), dqrfreez(klon,klev))
       ALLOCATE(dqsauto(klon,klev), dqsagg(klon,klev), dqsrim(klon,klev), dqsmelt(klon,klev), dqsfreez(klon,klev))
 
 #ifdef CPP_StratAer
+      ALLOCATE (d_q_emiss(klon,klev))
       ALLOCATE (R2SO4(klon,klev))
       ALLOCATE (DENSO4(klon,klev))
@@ -1147 +1204 @@
       ALLOCATE (OCS_lifetime(klon,klev))
       ALLOCATE (SO2_lifetime(klon,klev))
+      ALLOCATE (H2SO4_lifetime(klon,klev))
+      ALLOCATE (O3_clim(klon,klev))
       ALLOCATE (alpha_bin(nbands_sw_rrtm+nbands_lw_rrtm+nwave,nbtr))
       ALLOCATE (piz_bin(nbands_sw_rrtm+nbands_lw_rrtm+nwave,nbtr))
@@ -1180 +1239 @@
 USE indice_sol_mod
 IMPLICIT NONE
-      DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri, qbs_seri)
+      DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri, qbs_seri,qx_seri)
       DEALLOCATE(u_seri,v_seri)
       DEALLOCATE(l_mixmin,l_mix,wprime)
       DEALLOCATE(pbl_eps)
-
+      DEALLOCATE(rhcl)
       DEALLOCATE(tr_seri)
       DEALLOCATE(d_t_dyn,d_q_dyn)
@@ -1206 +1265 @@
       DEALLOCATE(d_u_ajs,d_v_ajs)
       DEALLOCATE(d_t_eva,d_q_eva)
+      DEALLOCATE(d_qx_eva)
       DEALLOCATE(d_ql_eva,d_qi_eva)
       DEALLOCATE(d_t_lscst,d_q_lscst)
@@ -1214 +1274 @@
       DEALLOCATE(d_t_bs,d_q_bs,d_qbs_bs)
 #ifdef ISO
-      deallocate(xt_seri,xtl_seri,xts_seri)
+      deallocate(xt_seri,xtl_seri,xts_seri,xtbs_seri)
       DEALLOCATE(d_xtl_eva,d_xti_eva)
-      deallocate(d_xt_dyn,d_xtl_dyn,d_xts_dyn)
+      deallocate(d_xt_dyn,d_xtl_dyn,d_xts_dyn,d_xtbs_dyn)
       deallocate(d_xt_con)
       deallocate(d_xt_wake)
@@ -1306 +1366 @@
       DEALLOCATE(solsw_aerop, solsw0_aerop)
       DEALLOCATE(topswcf_aerop, solswcf_aerop)
-
+!AI Aerosols
+      DEALLOCATE(m_allaer)
 !CK LW diagnostics
       DEALLOCATE(toplwad_aerop, sollwad_aerop)
       DEALLOCATE(toplwai_aerop, sollwai_aerop)
       DEALLOCATE(toplwad0_aerop, sollwad0_aerop)
+
+!AI Ajout pour Ecrad (3Deffect)
+      DEALLOCATE(topswad_aero_s2, solswad_aero_s2)
+      DEALLOCATE(topswai_aero_s2, solswai_aero_s2)
+      DEALLOCATE(topswad0_aero_s2, solswad0_aero_s2)
+      DEALLOCATE(topsw_aero_s2, topsw0_aero_s2)
+      DEALLOCATE(solsw_aero_s2, solsw0_aero_s2)
+      DEALLOCATE(topswcf_aero_s2, solswcf_aero_s2)
+!CK LW diagnostics
+      DEALLOCATE(toplwad_aero_s2, sollwad_aero_s2)
+      DEALLOCATE(toplwai_aero_s2, sollwai_aero_s2)
+      DEALLOCATE(toplwad0_aero_s2, sollwad0_aero_s2)      
 
 ! FH Ajout de celles necessaires au phys_output_write_mod
@@ -1352 +1425 @@
       DEALLOCATE(zxfqcalving, zxfluxlat)
       DEALLOCATE(zxrunofflic)
-      DEALLOCATE(zxustartlic, zxrhoslic)
+      DEALLOCATE(zxustartlic, zxrhoslic, zxqsaltlic)
       DEALLOCATE(zxtsol, snow_lsc, zxfqfonte, zxqsurf)
       DEALLOCATE(rain_lsc)
@@ -1420 +1493 @@
       DEALLOCATE(upwd, dnwd)
       DEALLOCATE(ep)
-      DEALLOCATE(da, mp )
-      DEALLOCATE(phi )
-      DEALLOCATE(wght_cvfd )
-      DEALLOCATE(phi2 )
+      DEALLOCATE(da, mp)
+      DEALLOCATE(phi)
+      DEALLOCATE(wght_cvfd)
+      DEALLOCATE(phi2)
       DEALLOCATE(d1a, dam)
-      DEALLOCATE(ev )
-      DEALLOCATE(elij )
-      DEALLOCATE(qtaa )
-      DEALLOCATE(clw )
-      DEALLOCATE(epmlmMm, eplaMm )
-      DEALLOCATE(sij )
+      DEALLOCATE(ev)
+      DEALLOCATE(elij)
+      DEALLOCATE(qtaa)
+      DEALLOCATE(clw)
+      DEALLOCATE(epmlmMm, eplaMm)
+      DEALLOCATE(sij)
 #ifdef ISO
       DEALLOCATE(xtwdtrainA)
@@ -1470 +1543 @@
       DEALLOCATE(rneb)
       DEALLOCATE(pfraclr,pfracld)
+      DEALLOCATE (zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic)
       DEALLOCATE(distcltop)
       DEALLOCATE(temp_cltop)
-      DEALLOCATE (zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic)
 #ifdef ISO
       DEALLOCATE (zxxtsnow,xtVprecip,xtVprecipi,pxtrfl,pxtsfl)
@@ -1493 +1566 @@
 
 !--POPRECIP variables
-      DEALLOCATE(qraindiag, qsnowdiag) 
-      DEALLOCATE(dqreva,dqssub)
-      DEALLOCATE(dqrauto,dqrcol,dqrmelt,dqrfreez)
-      DEALLOCATE(dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez)
+      DEALLOCATE(qraindiag, qsnowdiag)
+      DEALLOCATE(dqreva, dqssub)
+      DEALLOCATE(dqrauto, dqrcol, dqrmelt, dqrfreez)
+      DEALLOCATE(dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
 
 #ifdef CPP_StratAer
 ! variables for strat. aerosol CK
+      DEALLOCATE (d_q_emiss)
       DEALLOCATE (R2SO4)
       DEALLOCATE (DENSO4)
@@ -1507 +1581 @@
       DEALLOCATE (SO2_lifetime)
       DEALLOCATE (OCS_lifetime)
+      DEALLOCATE (H2SO4_lifetime)
+      DEALLOCATE (O3_clim)
       DEALLOCATE (alpha_bin)
       DEALLOCATE (piz_bin)

LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90

@@ -39 +39 @@
     USE ioipsl_getin_p_mod, ONLY : getin_p
     USE indice_sol_mod
-    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac
+    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac,ivap,iliq,isol
     USE readTracFiles_mod, ONLY: addPhase
     USE strings_mod,  ONLY: strIdx
@@ -94 +94 @@
     USE phys_output_var_mod, ONLY : cloud_cover_sw, cloud_cover_sw_s2
 
+
     !USE cmp_seri_mod
 !    USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
@@ -118 +119 @@
     USE chem_rep, ONLY: Init_chem_rep_xjour, d_q_rep, d_ql_rep, d_qi_rep, &
                         ptrop, ttrop, ztrop, gravit, itroprep, Z1, Z2, fac, B
+    USE strataer_local_var_mod
+    USE strataer_emiss_mod, ONLY: strataer_emiss_init
 #endif
 #if defined INCA || defined REPROBUS
@@ -132 +135 @@
 
 #ifdef CPP_StratAer
+    USE phys_local_var_mod, ONLY: d_q_emiss
     USE strataer_local_var_mod
     USE strataer_nuc_mod, ONLY: strataer_nuc_init
     USE strataer_emiss_mod, ONLY: strataer_emiss_init
 #endif
-
 
     USE lmdz_xios, ONLY: xios_update_calendar, xios_context_finalize
@@ -154 +157 @@
         & modif_ratqs,essai_convergence,iso_init,ridicule_rain,tnat, &
         & ridicule,ridicule_snow
-    USE isotopes_routines_mod, ONLY: iso_tritium
+    USE isotopes_routines_mod, ONLY: iso_tritium,dispatch,together
 #ifdef ISOVERIF
     USE isotopes_verif_mod, ONLY: errmax,errmaxrel, &
@@ -189 +192 @@
 !!!!!!!!!!!!!!!!!!  END "USE" for CPP keys !!!!!!!!!!!!!!!!!!!!!!
 
+USE physiqex_mod, ONLY : physiqex
 USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
        ! [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, &
+       qx_seri, & ! CR
+       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, &
@@ -207 +213 @@
        !
        d_t_eva,d_q_eva,d_ql_eva,d_qi_eva, &
+       d_qx_eva, &
        d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc, &
        d_t_lscst,d_q_lscst, &
@@ -245 +252 @@
        toplwai_aero,sollwai_aero,   &
        toplwad0_aero,sollwad0_aero, &
+       !pour Ecrad
+       topswad_aero_s2, solswad_aero_s2,   &
+       topswai_aero_s2, solswai_aero_s2,   &
+       topswad0_aero_s2, solswad0_aero_s2, &
+       topsw_aero_s2, topsw0_aero_s2,      &
+       solsw_aero_s2, solsw0_aero_s2,      &
+       topswcf_aero_s2, solswcf_aero_s2,   &
+       !LW diagnostics
+       toplwad_aero_s2, sollwad_aero_s2,   &
+       toplwai_aero_s2, sollwai_aero_s2,   &
+       toplwad0_aero_s2, sollwad0_aero_s2, &
        !
        topsw_aero,solsw_aero,       &
@@ -263 +281 @@
        toplwai_aerop, sollwai_aerop,   &
        toplwad0_aerop, sollwad0_aerop, &
+       !pour Ecrad
+       topswad_aero_s2, solswad_aero_s2,   &
+       topswai_aero_s2, solswai_aero_s2,   &
+       topswad0_aero_s2, solswad0_aero_s2, &
+       topsw_aero_s2, topsw0_aero_s2,      &
+       solsw_aero_s2, solsw0_aero_s2,      &
+       topswcf_aero_s2, solswcf_aero_s2,   &
+       !LW diagnostics
+       toplwad_aero_s2, sollwad_aero_s2,   &
+       toplwai_aero_s2, sollwai_aero_s2,   &
+       toplwad0_aero_s2, sollwad0_aero_s2, &
        !
        ptstar, pt0, slp, &
@@ -271 +300 @@
        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,                   &
@@ -359 +388 @@
        t2m, fluxlat,  &
        fsollw, evap_pot,  &
-       fsolsw, wfbils, wfevap,  &
+       fsolsw, wfbils, wfevap, &  
        prfl, psfl,bsfl, fraca, Vprecip,  &
        zw2,  &
@@ -373 +402 @@
        rneb,  &
        zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic, &
-       zxfluxt,zxfluxq
-
-
-#ifdef ISO
-       USE phys_local_var_mod, ONLY: xt_seri,xtl_seri,xts_seri, &
+       zxfluxt,zxfluxq 
+
+
+#ifdef ISO
+       USE phys_local_var_mod, ONLY: xt_seri,xtl_seri,xts_seri,xtbs_seri, &
        d_xt_eva,d_xtl_eva,d_xti_eva,           &
-       d_xt_dyn,d_xtl_dyn,d_xts_dyn,           &
+       d_xt_dyn,d_xtl_dyn,d_xts_dyn,d_xtbs_dyn, &
        d_xt_con, d_xt_wake,                    &
        d_xt_ajsb, d_xt_ajs,                    &
@@ -405 +434 @@
        USE phys_output_var_mod, ONLY: scdnc, cldncl, reffclwtop, lcc, reffclws, &
        reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra
+       USE output_physiqex_mod, ONLY: output_physiqex
 
 
@@ -549 +579 @@
     !
     ! 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, irneb, ibs
+!!$OMP THREADPRIVATE(ivap, iliq, isol, irneb, ibs)
+! Camille Risi 25 juillet 2023: ivap,iliq,isol deja definis dans infotrac_phy. 
+! Et ils sont utiles ailleurs que dans physiq_mod (ex:
+! reevap -> je commente les 2 lignes au dessus et je laisse la definition
+! plutot dans infotrac_phy
+    INTEGER,SAVE :: irneb, ibs
+!$OMP THREADPRIVATE(irneb, ibs)
+!
     !
     ! Variables argument:
@@ -805 +841 @@
     real therm_tke_max0(klon)   ! TKE dans les thermiques au LCL 
     real env_tke_max0(klon)     ! TKE dans l'environnement au LCL 
+    INTEGER, SAVE :: iflag_thermcell_tke ! transtport TKE by thermals
+    !$OMP THREADPRIVATE(iflag_thermcell_tke)
 
 !JLD    !---D\'eclenchement stochastique
@@ -893 +931 @@
     EXTERNAL ajsec     ! ajustement sec
     EXTERNAL conlmd    ! convection (schema LMD)
-    !KE43
     EXTERNAL conema3  ! convect4.3
-    !AA 
-    ! JBM (3/14) fisrtilp_tr not loaded
-    ! EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
-    !                          ! stockage des coefficients necessaires au
-    !                          ! lessivage OFF-LINE et ON-LINE
     EXTERNAL hgardfou  ! verifier les temperatures
     EXTERNAL nuage     ! calculer les proprietes radiatives
@@ -921 +953 @@
     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !
-    REAL rhcl(klon,klev)    ! humiditi relative ciel clair
+!    REAL rhcl(klon,klev)    ! humiditi relative ciel clair
     REAL dialiq(klon,klev)  ! eau liquide nuageuse
     REAL diafra(klon,klev)  ! fraction nuageuse
@@ -953 +985 @@
     REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
     !
-#ifdef INCA
-    REAL zxsnow_dummy(klon)
-#endif
     REAL zsav_tsol(klon)
     !
@@ -1061 +1090 @@
     !$OMP THREADPRIVATE(ok_bug_split_th)
 
+    ! Logical switch to a bug : modifying directly wake_deltat  by adding 
+    ! the (w) dry adjustment tendency to wake_deltat
+    LOGICAL, SAVE :: ok_bug_ajs_cv = .TRUE.
+    !$OMP THREADPRIVATE(ok_bug_ajs_cv)
     !
     !********************************************************
@@ -1265 +1298 @@
     ! Declarations pour Simulateur COSP
     !============================================================
+    ! AI 10-22
+#ifdef CPP_COSP
+    include "ini_COSP.h"
+#endif
+#ifdef CPP_COSPV2
+    include "ini_COSP.h"
+#endif
     real :: mr_ozone(klon,klev), phicosp(klon,klev)
 
@@ -1340 +1380 @@
 
     REAL, dimension(klon,klev) :: t_env,q_env
+#ifdef ISO
+    real, dimension(ntraciso,klon,klev) ::  xt_env
+#endif
 
     REAL, dimension(klon) :: pr_et
@@ -1350 +1393 @@
     !AI namelist pour gerer le double appel de Ecrad
     CHARACTER(len=512) :: namelist_ecrad_file
+
+    !======================================================================!
+    ! Bifurcation vers un nouveau moniteur physique pour experimenter      !
+    ! des solutions et préparer le couplage avec la physique de MesoNH     !
+    ! 14 mai 2023                                                          !
+    !======================================================================!
+    if (debut) then                                                        !
+       iflag_physiq=0
+       call getin_p('iflag_physiq', iflag_physiq)                          !
+    endif                                                                  !
+    if ( iflag_physiq == 2 ) then
+#ifdef ISO
+       abort_message='isotopes pas encore dans physiqex'
+       CALL abort_physic(modname,abort_message,1)
+#endif
+       call physiqex (nlon,nlev, &                                         !
+       debut,lafin,pdtphys_, &                                             !
+       paprs,pplay,pphi,pphis,presnivs, &                                  !
+       u,v,rot,t,qx, &                                                     !
+       flxmass_w, &                                                        !
+       d_u, d_v, d_t, d_qx, d_ps)                                          !
+       return                                                              !
+    endif                                                                  !
+    !======================================================================!
+
 
     pi = 4. * ATAN(1.)
@@ -1366 +1434 @@
     phys_tstep=NINT(pdtphys)
     IF (.NOT. using_xios) missing_val=nf90_fill_real
-#ifdef CPP_XIOS
-! switch to XIOS LMDZ physics context 
-    IF (.NOT. debut .AND. is_omp_master) THEN
-       CALL wxios_set_context()
-       CALL xios_update_calendar(itap+1)
+
+    IF (using_xios) THEN
+      ! switch to XIOS LMDZ physics context 
+      IF (.NOT. debut .AND. is_omp_master) THEN
+        CALL wxios_set_context()
+        CALL xios_update_calendar(itap+1)
+      ENDIF
     ENDIF
-#endif
 
     !======================================================================
@@ -1378 +1447 @@
     ! Utilise notamment en 1D mais peut etre active egalement en 3D
     ! en imposant la valeur de igout.
-    !======================================================================d
+    !======================================================================
     IF (prt_level.ge.1) THEN
        igout=klon/2+1/klon
@@ -1410 +1479 @@
        irneb= strIdx(tracers(:)%name, addPhase('H2O', 'r'))
        ibs  = strIdx(tracers(:)%name, addPhase('H2O', 'b'))
-       CALL init_etat0_limit_unstruct
-       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
+!       CALL init_etat0_limit_unstruct
+!       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
        !CR:nvelles variables convection/poches froides
 
@@ -1434 +1503 @@
             read_climoz, &
             alp_offset)
+       CALL init_etat0_limit_unstruct
+       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
        CALL phys_state_var_init(read_climoz)
        CALL phys_output_var_init
        IF (read_climoz>=1 .AND. create_etat0_limit .AND. grid_type==unstructured) & 
           CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
+
+#ifdef REPROBUS
+       CALL strataer_init
+       CALL strataer_emiss_init
+#endif
 
 #ifdef CPP_StratAer
@@ -1481 +1557 @@
 
         IF (ok_bs) THEN
+#ifdef ISO
           abort_message='blowing snow cannot be activated with water isotopes yet'
           CALL abort_physic(modname,abort_message, 1)
+#endif
          IF ((ok_ice_sursat.AND.nqo .LT.5).OR.(.NOT.ok_ice_sursat.AND.nqo.LT.4)) THEN
              WRITE (lunout, *) 'activation of blowing snow needs a specific H2O tracer', &
@@ -1490 +1568 @@
          ENDIF
         ENDIF
+
        Ncvpaseq1 = 0
        dnwd0=0.0
@@ -1531 +1610 @@
        tau_gl=86400.*tau_gl
        WRITE(lunout,*) 'debut physiq_mod tau_gl=',tau_gl
+       iflag_thermcell_tke=0
+       call getin_p('iflag_thermcell_tke', iflag_thermcell_tke)                          !
 
        CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond) 
@@ -1553 +1634 @@
        CALL getin_p('ok_bug_cv_trac',ok_bug_cv_trac)
        CALL getin_p('ok_bug_split_th',ok_bug_split_th)
+       CALL getin_p('ok_bug_ajs_cv',ok_bug_ajs_cv)
        fl_ebil = 0 ! by default, conservation diagnostics are desactivated
        CALL getin_p('fl_ebil',fl_ebil)
@@ -1589 +1671 @@
        CALL infocfields_init
 
+       !AI 08 2023
 #ifdef CPP_ECRAD
        ok_3Deffect=.false.
@@ -1868 +1951 @@
       IF (.NOT. create_etat0_limit) CALL init_readaerosolstrato(flag_aerosol_strat)  !! initialise aero strato from file for XIOS interpolation (unstructured_grid)
 
+      if (ok_cosp) then
 #ifdef CPP_COSP
-      IF (ok_cosp) THEN
-!           DO k = 1, klev
-!             DO i = 1, klon
-!               phicosp(i,k) = pphi(i,k) + pphis(i)
-!             ENDDO
-!           ENDDO
+        ! A.I : Initialisations pour le 1er passage a Cosp
+        CALL ini_COSP(ref_liq_cosp0,ref_ice_cosp0,pctsrf_cosp0,zu10m_cosp0,zv10m_cosp0, &
+               zxtsol_cosp0,zx_rh_cosp0,cldfra_cosp0,rnebcon_cosp0,flwc_cosp0, &
+               fiwc_cosp0,prfl_cosp0,psfl_cosp0,pmflxr_cosp0,pmflxs_cosp0, &
+               mr_ozone_cosp0,cldtau_cosp0,cldemi_cosp0,JrNt_cosp0)
+
         CALL phys_cosp(itap,phys_tstep,freq_cosp, &
+               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
+               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
+               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
+               JrNt_cosp0,ref_liq_cosp0,ref_ice_cosp0, &
+               pctsrf_cosp0, &
+               zu10m_cosp0,zv10m_cosp0,pphis, &
+               pphi,paprs(:,1:klev),pplay,zxtsol_cosp0,t, &
+               qx(:,:,ivap),zx_rh_cosp0,cldfra_cosp0,rnebcon_cosp0,flwc_cosp0,fiwc_cosp0, &
+               prfl_cosp0(:,1:klev),psfl_cosp0(:,1:klev), &
+               pmflxr_cosp0(:,1:klev),pmflxs_cosp0(:,1:klev), &
+               mr_ozone_cosp0,cldtau_cosp0, cldemi_cosp0)
+#endif
+
+#ifdef CPP_COSP2
+        CALL ini_COSP(ref_liq_cosp0,ref_ice_cosp0,pctsrf_cosp0,zu10m_cosp0,zv10m_cosp0, &
+               zxtsol_cosp0,zx_rh_cosp0,cldfra_cosp0,rnebcon_cosp0,flwc_cosp0, &
+               fiwc_cosp0,prfl_cosp0,psfl_cosp0,pmflxr_cosp0,pmflxs_cosp0, &
+               mr_ozone_cosp0,cldtau_cosp0,cldemi_cosp0,JrNt_cosp0)
+     
+        CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
                ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
                ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
@@ -1887 +1991 @@
                pmflxr(:,1:klev),pmflxs(:,1:klev), &
                mr_ozone,cldtau, cldemi)
-      ENDIF
-#endif
-
-#ifdef CPP_COSP2
-        IF (ok_cosp) THEN
-!           DO k = 1, klev
-!             DO i = 1, klon
-!               phicosp(i,k) = pphi(i,k) + pphis(i)
-!             ENDDO
-!           ENDDO
-          CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
-               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
-               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
-               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
-               JrNt,ref_liq,ref_ice, &
-               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
-               zu10m,zv10m,pphis, &
-               zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
-               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
-               prfl(:,1:klev),psfl(:,1:klev), &
-               pmflxr(:,1:klev),pmflxs(:,1:klev), &
-               mr_ozone,cldtau, cldemi)
-       ENDIF
 #endif
 
 #ifdef CPP_COSPV2
-        IF (ok_cosp) THEN
-           DO k = 1, klev
-             DO i = 1, klon
-               phicosp(i,k) = pphi(i,k) + pphis(i)
-             ENDDO
-           ENDDO
           CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
                ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
@@ -1931 +2006 @@
                pmflxr(:,1:klev),pmflxs(:,1:klev), &
                mr_ozone,cldtau, cldemi)
-       ENDIF
-#endif
+#endif
+      ENDIF
 
        !
@@ -2006 +2081 @@
 
 
-#ifdef CPP_XIOS
-       IF (is_omp_master) CALL xios_update_calendar(1)
-#endif
+       IF (using_xios) THEN
+         IF (is_omp_master) CALL xios_update_calendar(1)
+       ENDIF
+       
        IF(read_climoz>=1 .AND. create_etat0_limit) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
        CALL create_etat0_limit_unstruct
@@ -2211 +2287 @@
         IF (is_omp_master) CALL xios_get_field_attr("temp",default_value=missing_val)
         CALL bcast_omp(missing_val)
-       
+
        ! 
        ! Now we activate some double radiation call flags only if some
@@ -2539 +2615 @@
           u_seri(i,k)  = u(i,k)
           v_seri(i,k)  = v(i,k)
+          qx_seri(i,k,:)  = qx(i,k,:)
           q_seri(i,k)  = qx(i,k,ivap)
           ql_seri(i,k) = qx(i,k,iliq)
@@ -2552 +2629 @@
              qs_seri(i,k) = qx(i,k,isol)
              IF (ok_ice_sursat) THEN
-             rneb_seri(i,k) = qx(i,k,irneb)
+               rneb_seri(i,k) = qx(i,k,irneb)
              ENDIF
              IF (ok_bs) THEN
-             qbs_seri(i,k)= qx(i,k,ibs)
+               qbs_seri(i,k)= qx(i,k,ibs)
              ENDIF
-
           ENDIF
-
-
        ENDDO
     ENDDO
@@ -2582 +2656 @@
       DO k = 1, klev
        DO i = 1, klon
+          xtbs_seri(ixt,i,k) = 0.
           xt_seri(ixt,i,k)  = qx(i,k,iqIsoPha(ixt,ivap))
           xtl_seri(ixt,i,k) = qx(i,k,iqIsoPha(ixt,iliq))
@@ -2602 +2677 @@
     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
 #ifdef ISO
-#ifdef ISOVERIF
-        write(*,*) 'physiq tmp 1913'
-#endif
     do ixt=1,ntraciso
     xtql1(ixt,:)=0.0
     DO k = 1, klev
-      xtql1(ixt,:)=xtql1(ixt,:)+(xt_seri(ixt,:,k)+xtl_seri(ixt,:,k)+xts_seri(ixt,:,k))*zmasse(:,k)
-    ENDDO
+      xtql1(ixt,:)=xtql1(ixt,:)+(xt_seri(ixt,:,k)+xtl_seri(ixt,:,k))*zmasse(:,k)
+      IF (nqo >= 3) THEN
+        xtql1(ixt,:)=xtql1(ixt,:)+xts_seri(ixt,:,k)*zmasse(:,k)
+      ENDIF
+      IF (ok_bs) THEN
+        xtql1(ixt,:)=xtql1(ixt,:)+xtbs_seri(ixt,:,k)*zmasse(:,k)
+      ENDIF
+    ENDDO !DO k = 1, klev
     enddo !do ixt=1,ntraciso
 #endif
@@ -2625 +2709 @@
          IF(.NOT.tracers(iq)%isInPhysics) CYCLE
          itr = itr+1
-!#ifdef ISOVERIF
-!         write(*,*) 'physiq 1973: itr,iq=',itr,iq
-!         write(*,*) 'qx(1,1,iq)=',qx(1,1,iq)
-!#endif
-         DO  k = 1, klev
+          DO  k = 1, klev
              DO  i = 1, klon
                 tr_seri(i,k,itr) = qx(i,k,iq)
@@ -2744 +2824 @@
               d_xts_dyn(ixt,i,k) =  &
      &           (xts_seri(ixt,i,k)-xts_ancien(ixt,i,k))/phys_tstep
+              d_xtbs_dyn(ixt,i,k) =  &
+     &           (xtbs_seri(ixt,i,k)-xtbs_ancien(ixt,i,k))/phys_tstep
             enddo ! do ixt=1,ntraciso
          ENDDO
@@ -2757 +2839 @@
            call iso_verif_noNaN(d_xtl_dyn(ixt,i,k),'physiq 2220d')
            call iso_verif_noNaN(d_xts_dyn(ixt,i,k),'physiq 2220e')
+           call iso_verif_noNaN(d_xtbs_dyn(ixt,i,k),'physiq 2220f')
            enddo ! do ixt=1,ntraciso
          enddo
@@ -2839 +2922 @@
        d_qbs_dyn(:,:)=0.0
        ancien_ok = .TRUE.
+#ifdef ISO
+       d_xtbs_dyn(:,:,:)=0.0
+#endif
     ENDIF
     !
@@ -2962 +3048 @@
     ! Re-evaporer l'eau liquide nuageuse
     !
-     CALL reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
-   &         d_t_eva,d_q_eva,d_ql_eva,d_qi_eva &
-#ifdef ISO
-             ,xt_seri,xtl_seri,xts_seri,d_xt_eva,d_xtl_eva,d_xti_eva &
-#endif
-   &     )
+     CALL reevap (klon,klev,iflag_ice_thermo,t_seri,qx_seri, &
+   &         d_t_eva,d_qx_eva)
+
+     call dispatch(klon,klev,qx_seri,q_seri,xt_seri,ql_seri,xtl_seri,qs_seri,xts_seri)
+     call dispatch(klon,klev,d_qx_eva,d_q_eva,d_xt_eva,d_ql_eva,d_xtl_eva,d_qi_eva,d_xti_eva)
+
+
+#ifdef ISO
+#ifdef ISOVERIF
+ DO k = 1, klev 
+     DO i = 1, klon
+      do ixt=1,ntraciso
+      call iso_verif_noNaN(xt_seri(ixt,i,k), &
+     &     'reevap 2417: apres evap tot')
+      enddo
+      if (iso_eau.gt.0) then
+              call iso_verif_egalite_choix( &
+     &           xt_seri(iso_eau,i,k),q_seri(i,k), &
+     &          'reevap 1891, après réévap totale',errmax,errmaxrel)
+              call iso_verif_egalite_choix( &
+     &           xtl_seri(iso_eau,i,k),ql_seri(i,k), &
+     &          'reevap 2209, après réévap totale',errmax,errmaxrel)
+              call iso_verif_egalite_choix( &
+     &           xts_seri(iso_eau,i,k),qs_seri(i,k), &
+     &          'reevap 2209b, après réévap totale',errmax,errmaxrel)
+       endif !if (iso_eau.gt.0) then
+     
+      if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then  
+            if (q_seri(i,k).gt.ridicule) then  
+               if (iso_verif_o18_aberrant_nostop( &
+     &              xt_seri(iso_HDO,i,k)/q_seri(i,k), &
+     &              xt_seri(iso_O18,i,k)/q_seri(i,k), &
+     &              'reevap 2408: apres reevap totale').eq.1) then
+                  write(*,*) 'i,k,q_seri(i,k)=',i,k,q_seri(i,k)                       
+                  stop
+              endif !  if (iso_verif_o18_aberrant_nostop
+            endif !if (q_seri(i,k).gt.errmax) then   
+        endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then        
+#ifdef ISOTRAC      
+             call iso_verif_traceur(xt_seri(1,i,k), &
+     &           'reevap 2165c')
+             call iso_verif_traceur_pbidouille(xt_seri(1,i,k), &
+     &           'reevap 2165d')
+#endif
+       ENDDO
+    ENDDO               
+#endif                  
+#endif
+
 
      CALL add_phys_tend &
@@ -3099 +3228 @@
     ! Calcul de l'humidite de saturation au niveau du sol
 
+! Tests Fredho, instensibilite au pas de temps -------------------------------
+! A detruire en 2024 une fois les tests documentes et les choix faits        !
+! Conservation des variables avant l'appel à l a diffusion pour les tehrmic  !
+    if (iflag_thermals_tenv / 10 == 1 ) then                                 !
+        do k=1,klev                                                          !
+           do i=1,klon                                                       !
+              t_env(i,k)=t_seri(i,k)                                         !
+              q_env(i,k)=q_seri(i,k)   
+#ifdef ISO
+              do ixt=1,ntraciso
+                xt_env(ixt,i,k)=xt_seri(ixt,i,k)
+              enddo
+#endif
+           enddo                                                             !
+        enddo                                                                !
+    else if (iflag_thermals_tenv / 10 == 2 ) then                            !
+        do k=1,klev                                                          !
+           do i=1,klon                                                       !
+              t_env(i,k)=t_seri(i,k)                                         !
+           enddo                                                             !
+        enddo                                                                !
+    endif                                                                    !
+! Tests Fredho, instensibilite au pas de temps -------------------------------
 
 
@@ -3287 +3439 @@
           d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
           CALL add_wake_tend &
-             (d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy &
+             (d_deltat_vdf, d_deltaq_vdf, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy &
 #ifdef ISO
                ,d_deltaxt_vdf &
@@ -3320 +3472 @@
      &   )
        ENDIF
-#ifdef ISOVERIF
-        write(*,*) 'physiq tmp 2736'
-#endif
-
        CALL prt_enerbil('vdf',itap)
+
        !--------------------------------------------------------------------
 
@@ -3689 +3838 @@
                 ENDDO
              ENDDO
-             IF (iflag_adjwk == 2) THEN
+             IF (iflag_adjwk == 2 .AND. OK_bug_ajs_cv) THEN
                CALL add_wake_tend &
-                 (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy &
+                 (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy &
 #ifdef ISO
                       ,d_deltaxt_ajs_cv  &
 #endif
                  ) 
-             ENDIF  ! (iflag_adjwk == 2)
+             ENDIF  ! (iflag_adjwk == 2 .AND. OK_bug_ajs_cv)
           ENDIF  ! (iflag_adjwk >= 1)
        ENDIF ! (iflag_wake>=1)
@@ -4400 +4549 @@
        !  ====
        IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
+       WRITE(lunout,*) 'WARNING : running without dry convection. Somme intermediate variables are not properly defined in physiq_mod.F90'
+       ! Reprendre proprement les initialisation ci dessouds si on veut vraiment utiliser l'option (FH)
+          fraca(:,:)=0.
+          fm_therm(:,:)=0.
+          ztv(:,:)=t_seri(:,:)
+          zqasc(:,:)=q_seri(:,:)
+          ztla(:,:)=0.
+          zthl(:,:)=0.
+          zpspsk(:,:)=(pplay(:,:)/100000.)**RKAPPA
 
 
@@ -4491 +4649 @@
 
        IF (iflag_thermals>=1) THEN
+
+! Tests Fredho, instensibilite au pas de temps -------------------------------
+! A detruire en 2024 une fois les tests documentes et les choix faits        !
+          if (iflag_thermals_tenv /10 == 0 ) then                            !
+            do k=1,klev                                                      !
+               do i=1,klon                                                   !
+                  t_env(i,k)=t_seri(i,k)                                     !
+                  q_env(i,k)=q_seri(i,k)                                     !
+#ifdef ISO
+                  do ixt=1,ntraciso
+                        xt_env(ixt,i,k)=xt_seri(ixt,i,k)
+                  enddo
+#endif
+               enddo                                                         !
+            enddo                                                            !
+          else if (iflag_thermals_tenv / 10 == 2 ) then                      !
+            do k=1,klev                                                      !
+               do i=1,klon                                                   !
+                  q_env(i,k)=q_seri(i,k)                                     !
+#ifdef ISO
+                  do ixt=1,ntraciso
+                        xt_env(ixt,i,k)=xt_seri(ixt,i,k)
+                  enddo
+#endif
+               enddo                                                         !
+            enddo                                                            !
+          else if (iflag_thermals_tenv / 10 == 3 ) then                      !
+            do k=1,klev                                                      !
+               do i=1,klon                                                   !
+                  t_env(i,k)=t(i,k)                                          !
+                  q_env(i,k)=qx(i,k,1)                                       !
+#ifdef ISO
+                  do ixt=1,ntraciso
+                        xt_env(ixt,i,k)=xt_seri(ixt,i,k)
+                  enddo
+#endif
+               enddo                                                         !
+            enddo                                                            !
+          endif                                                              !
+! Tests Fredho, instensibilite au pas de temps ------------------------------
+
           !jyg<
 !!       IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
@@ -4499 +4698 @@
                    t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
                    q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
+                   t_env(i,k)   = t_env(i,k) - wake_s(i)*wake_deltat(i,k)
+                   q_env(i,k)   = q_env(i,k) - wake_s(i)*wake_deltaq(i,k)
                    u_therm(i,k) = u_seri(i,k)
                    v_therm(i,k) = v_seri(i,k)
@@ -4504 +4705 @@
                    do ixt=1,ntraciso
                      xt_therm(ixt,i,k) = xt_seri(ixt,i,k) - wake_s(i)*wake_deltaxt(ixt,i,k)
+                     xt_env(ixt,i,k) = xt_env(ixt,i,k) - wake_s(i)*wake_deltaxt(ixt,i,k)
                    enddo !do ixt=1,ntraciso
 #endif
@@ -4546 +4748 @@
                ,pplay,paprs,pphi,weak_inversion &
                         ! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
-               ,u_therm,v_therm,t_therm,q_therm,t_therm,q_therm,zqsat,debut &  !jyg
+               ,u_therm,v_therm,t_therm,q_therm,t_env,q_env,zqsat,debut &  !jyg
                ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
                ,fm_therm,entr_therm,detr_therm &
@@ -4565 +4767 @@
                ,zqla,ztva &
 #ifdef ISO         
-     &      ,xt_therm,d_xt_ajs &
+     &      ,xt_env,d_xt_ajs &
 #ifdef DIAGISO         
      &      ,q_the,xt_the &
@@ -4600 +4802 @@
              IF (ok_bug_split_th) THEN
                CALL add_wake_tend &
-                   (d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy &
+                   (d_deltat_the, d_deltaq_the, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy &
 #ifdef ISO
                    ,d_deltaxt_the &
@@ -4607 +4809 @@
              ELSE
                CALL add_wake_tend &
-                   (d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wake_k, 'the', abortphy &
+                   (d_deltat_the, d_deltaq_the, dsig0, dsig0, ddens0, ddens0, wake_k, 'the', abortphy &
 #ifdef ISO
                    ,d_deltaxt_the &
@@ -4636 +4838 @@
           ! Transport de la TKE par les panaches thermiques.
           ! FH : 2010/02/01
-          !     if (iflag_pbl.eq.10) then
-          !     call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
-          !    s           rg,paprs,pbl_tke)
-          !     endif
+               if (iflag_thermcell_tke==1) then
+               call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,rg,paprs,pbl_tke)
+               endif
           ! -------------------------------------------------------------------
 
@@ -4895 +5096 @@
     ELSE
 
+    ! Camille Risi mai 2024: on ne met pas à jour ici pour ne pas s'mbêter à modifier fisrtilp
     CALL fisrtilp(phys_tstep,paprs,pplay, &
          t_seri, q_seri,ptconv,ratqs, &
@@ -5560 +5762 @@
                 !
              ENDIF
+<<<<<<< .mine
+          ELSE IF (iflag_rrtm .EQ.2) THEN    ! ecrad RADIATION
+#ifdef CPP_ECRAD
+             !--climatologies or INCA aerosols
+             CALL readaerosol_optic_ecrad( debut, aerosol_couple, ok_alw, ok_volcan, &
+                  flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
+                  pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
+                  tr_seri, mass_solu_aero, mass_solu_aero_pi)  
+#else
+                abort_message='You should compile with -rad ecrad if running with iflag_rrtm=2'
+                CALL abort_physic(modname,abort_message,1)
+#endif
+||||||| .r4942
+=======
           ELSE IF (iflag_rrtm .EQ.2) THEN    ! ecrad RADIATION
 #ifdef CPP_ECRAD
@@ -5571 +5787 @@
                 CALL abort_physic(modname,abort_message,1)
 #endif
+>>>>>>> .r4981
           ENDIF
+
        ELSE   !--flag_aerosol = 0 
           tausum_aero(:,:,:) = 0.
@@ -5873 +6091 @@
                 ENDIF
                 !
-                ! AI namelist utilise pour l appel principal de radlwsw (ecrad)
                 namelist_ecrad_file='namelist_ecrad'
                 !
@@ -5917 +6134 @@
                      cloud_cover_sw)
           ENDIF !ok_4xCO2atm
+
+! A.I aout 2023
+! Effet 3D des nuages Ecrad
+! a passer : nom du ficher namelist et cles ok_3Deffect
+! a declarer comme iflag_rrtm et a lire dans physiq.def
+#ifdef CPP_ECRAD
+          IF (ok_3Deffect) then
+!                print*,'ok_3Deffect = ',ok_3Deffect  
+                namelist_ecrad_file='namelist_ecrad_s2'
+                CALL radlwsw &
+                     (debut, dist, rmu0, fract,  &
+                     paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
+                     t_seri,q_seri,wo, &
+                     cldfrarad, cldemirad, cldtaurad, &
+                     ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie,  ok_volcan, flag_volc_surfstrat, &
+                     flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
+                     tau_aero, piz_aero, cg_aero, &
+                     tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
+                     tau_aero_lw_rrtm, &
+                     cldtaupi, &
+                     zqsat, flwc, fiwc, &
+                     ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
+                     namelist_ecrad_file, &
+! A modifier              
+                     heat_s2,heat0_s2,cool_s2,cool0_s2,albpla_s2, &
+                     heat_volc,cool_volc, &
+                     topsw_s2,toplw_s2,solsw_s2,solswfdiff_s2,sollw_s2, &
+                     sollwdown_s2, &
+                     topsw0_s2,toplw0_s2,solsw0_s2,sollw0_s2, &
+                     lwdnc0_s2, lwdn0_s2, lwdn_s2, lwupc0_s2, lwup0_s2, lwup_s2,  &
+                     swdnc0_s2, swdn0_s2, swdn_s2, swupc0_s2, swup0_s2, swup_s2, &
+                     topswad_aero_s2, solswad_aero_s2, &
+                     topswai_aero_s2, solswai_aero_s2, &
+                     topswad0_aero_s2, solswad0_aero_s2, &
+                     topsw_aero_s2, topsw0_aero_s2, &
+                     solsw_aero_s2, solsw0_aero_s2, &
+                     topswcf_aero_s2, solswcf_aero_s2, &
+                                !-C. Kleinschmitt for LW diagnostics
+                     toplwad_aero_s2, sollwad_aero_s2,&
+                     toplwai_aero_s2, sollwai_aero_s2, &
+                     toplwad0_aero_s2, sollwad0_aero_s2,&
+                                !-end
+                     ZLWFT0_i, ZFLDN0, ZFLUP0, &
+                     ZSWFT0_i, ZFSDN0, ZFSUP0, &
+                     cloud_cover_sw_s2)
+          ENDIF ! ok_3Deffect
+#endif
+
        ENDIF ! aerosol_couple
        itaprad = 0
@@ -6140 +6405 @@
        d_t_hin(:, :)=0.
        CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
-            dqi0, dqbs0,paprs, 'hin', abortphy,flag_inhib_tend,itap,0 &
+            dqi0, dqbs0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
      &    ,dxt0,dxtl0,dxti0 &
@@ -6263 +6528 @@
        
        SELECT CASE(flag_emit)
-       CASE(1) ! emission volc H2O dans LMDZ
+       CASE(1) ! emission volc H2O in LMDZ
           DO ieru=1, nErupt
              IF (year_cur==year_emit_vol(ieru).AND.&
@@ -6271 +6536 @@
                 
                 IF(flag_verbose_strataer) print *,'IN physiq_mod: date=',year_cur,mth_cur,day_cur
-                ! initialisation tendance q emission
+                ! initialisation of q tendency emission
                 d_q_emiss(:,:)=0.
                 ! daily injection mass emission - NL
@@ -6278 +6543 @@
                 !
                 CALL STRATEMIT(pdtphys,pdtphys,latitude_deg,longitude_deg,t_seri,&
-                     pplay,paprs,tr_seri,&
-                     m_H2O_emiss_vol_daily,&
-                     xlat_min_vol(ieru),xlat_max_vol(ieru),&
-                     xlon_min_vol(ieru),xlon_max_vol(ieru),&
-                     altemiss_vol(ieru),&
-                     sigma_alt_vol(ieru),1,&
-                     1,nAerErupt+1,0)
+                    pplay,paprs,tr_seri,&
+                    m_H2O_emiss_vol_daily,&
+                    xlat_min_vol(ieru),xlat_max_vol(ieru),&
+                    xlon_min_vol(ieru),xlon_max_vol(ieru),&
+                    altemiss_vol(ieru),sigma_alt_vol(ieru),1,1.,&
+                    nAerErupt+1,0)
                 
                 IF(flag_verbose_strataer) print *,'IN physiq_mod: min max d_q_emiss=',&
@@ -6298 +6562 @@
     ENDIF
 #endif
-
 
 !===============================================================
@@ -6737 +7000 @@
     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
+#ifdef ISO
+      CALL abort_gcm('physiq 6936','isos pas prevus dans le mass fixer',1)
+      ! Camille Risi mai 2024: on attend d'avoir la 4e dimension qui rendra tout plus simple.
+#endif
+    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
-
-#ifdef ISO
-      DO i = 1, klon
-      do ixt=1,ntraciso
-       xtprw(ixt,i) = 0. 
-       DO k = 1, klev
-        xtprw(ixt,i) = xtprw(ixt,i) + &
-     &           xt_seri(ixt,i,k)*(paprs(i,k)-paprs(i,k+1))/RG
-       ENDDO !DO k = 1, klev
-      enddo !do ixt=1,ntraciso
-      enddo !DO i = 1, klon
-#endif
+    water_budget(:)=water_budget(:)+(rain_fall(:)+snow_fall(:)-evap(:))*pdtphys
+    IF (ok_bs) THEN
+      water_budget(:)=water_budget(:)+bs_fall(:)*pdtphys
+    ENDIF
+    ! Camille Risi mai 2024: pour les isotopes, on attend d'avoir la 4e dimension, ça rendra tout plus facile
+    ! ces variables sont diagnostiques, donc pas indispensables
+
+    !=======================================================================
+    !   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
@@ -6881 +7204 @@
     !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))*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
-      ENDDO
-    ENDDO
-#ifdef ISO
-    do ixt=1,ntraciso
-    xtql2(ixt,:)=0.0
-    DO k = 1, klev
-      xtql2(ixt,:)=xtql2(ixt,:)+(xt_seri(ixt,:,k)+xtl_seri(ixt,:,k)+xts_seri(ixt,:,k))*zmasse(:,k)
-    ENDDO
-    DO i = 1, klon
-      !--compute ratio of what q+ql should be with conservation to what it is
-      corrxtql(ixt)=(xtql1(ixt,i)+(xtevap(ixt,i)-xtrain_fall(ixt,i)-xtsnow_fall(ixt,i))*pdtphys)/xtql2(ixt,i)
-      DO k = 1, klev
-        xt_seri(ixt,i,k) =xt_seri(ixt,i,k)*corrxtql(ixt)
-        xtl_seri(ixt,i,k)=xtl_seri(ixt,i,k)*corrxtql(ixt)
-      ENDDO
-    ENDDO   
-    enddo !do ixt=1,ntraciso 
-#endif
-    ENDIF
-    !--fin mass fixer
-
     ! Sauvegarder les valeurs de t et q a la fin de la physique:
     !
@@ -6944 +7218 @@
     xtl_ancien(:,:,:)=xtl_seri(:,:,:)
     xts_ancien(:,:,:)=xts_seri(:,:,:)
+    xtbs_ancien(:,:,:)=xtbs_seri(:,:,:)
 #endif
     CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
@@ -7080 +7355 @@
          ok_sync, ptconv, read_climoz, clevSTD,          &
          ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
-         flag_aerosol, flag_aerosol_strat, ok_cdnc,t,u1,v1)
+         flag_aerosol, flag_aerosol_strat, ok_cdnc,t, u1, v1)
 #endif
 
 #ifndef CPP_XIOS
-    CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
-#endif
-
-#endif
-
-! Pour XIOS : On remet des variables a .false. apres un premier appel
-    IF (debut) THEN
-
-      IF (using_xios) THEN
-        swaero_diag=.FALSE.
-        swaerofree_diag=.FALSE.
-        dryaod_diag=.FALSE.
-        ok_4xCO2atm= .FALSE.
-!       write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
-
-        IF (is_master) THEN
-          !--setting up swaero_diag to TRUE in XIOS case 
-          IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. & 
-             xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. & 
-             xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR.  & 
-               (iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. & 
-                                   xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0"))))  & 
-             !!!--for now these fields are not in the XML files so they are omitted 
-             !!!  xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) & 
-             swaero_diag=.TRUE. 
-
-          !--setting up swaerofree_diag to TRUE in XIOS case 
-          IF (xios_field_is_active("SWdnSFCcleanclr").OR.xios_field_is_active("SWupSFCcleanclr").OR. &
-             xios_field_is_active("SWupTOAcleanclr").OR.xios_field_is_active("rsucsaf").OR.   &
-             xios_field_is_active("rsdcsaf") .OR. xios_field_is_active("LWdnSFCcleanclr").OR. &
-             xios_field_is_active("LWupTOAcleanclr")) &
-             swaerofree_diag=.TRUE. 
-
-          !--setting up dryaod_diag to TRUE in XIOS case 
-          DO naero = 1, naero_tot-1
-           IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE. 
-          ENDDO
-          !
-          !--setting up ok_4xCO2atm to TRUE in XIOS case 
-          IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. & 
-             xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
-             xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
-             xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
-             xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
-             xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
-             ok_4xCO2atm=.TRUE. 
-        ENDIF
-        !$OMP BARRIER
-        CALL bcast(swaero_diag)
-        CALL bcast(swaerofree_diag)
-        CALL bcast(dryaod_diag)
-        CALL bcast(ok_4xCO2atm)
-!        write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
-      ENDIF !using_xios
-    ENDIF
+      CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
+#endif
+
+#endif
+    ! Petit appelle de sorties pour accompagner le travail sur phyex
+    if ( iflag_physiq == 1 ) then
+        call output_physiqex(debut,jD_eq,pdtphys,presnivs,paprs,u,v,t,qx,cldfra,0.*t,0.*t,0.*t,pbl_tke,theta)
+    endif
 
     !====================================================================
@@ -7182 +7410 @@
     ! Disabling calls to the prt_alerte function
     alert_first_call = .FALSE.
+
     
     IF (lafin) THEN
@@ -7194 +7423 @@
          IF (read_climoz >= 1) THEN
            IF (is_mpi_root) CALL nf95_close(ncid_climoz)
-            DEALLOCATE(press_edg_climoz) ! pointer
-            DEALLOCATE(press_cen_climoz) ! pointer
+            DEALLOCATE(press_edg_climoz)
+            DEALLOCATE(press_cen_climoz)
          ENDIF
        
        ENDIF
+
+       IF (using_xios) THEN
+
+#ifdef INCA
+          IF (type_trac == 'inca') THEN 
+             IF (is_omp_master .AND. grid_type==unstructured) THEN 
+                CALL finalize_inca
+             ENDIF
+          ENDIF
+#endif
+
+          IF (is_omp_master .and. grid_type==unstructured) CALL xios_context_finalize
+       ENDIF
+
+       WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
        
-       IF (using_xios) THEN
-         IF (is_omp_master) CALL xios_context_finalize
-
-#ifdef INCA
-         if (type_trac == 'inca') then 
-            IF (is_omp_master .and. grid_type==unstructured) THEN 
-               CALL finalize_inca
-            ENDIF
-         endif
-#endif
-       ENDIF !using_xios
-       WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
     ENDIF
 

LMDZ6/trunk/libf/phylmdiso/reevap.F90

@@ -1 @@
-  SUBROUTINE reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
-   &         d_t_eva,d_q_eva,d_ql_eva,d_qs_eva &
-#ifdef ISO
-             ,xt_seri,xtl_seri,xts_seri,d_xt_eva,d_xtl_eva,d_xts_eva &
-#endif
-   &     )
+  SUBROUTINE reevap (klon,klev,iflag_ice_thermo,t_seri,qx, &
+   &         d_t_eva,d_qx_eva)
 
     ! flag to include modifications to ensure energy conservation (if flag >0)
     USE add_phys_tend_mod, only : fl_cor_ebil 
 #ifdef ISO
-    USE infotrac_phy, ONLY: ntiso   
+    USE infotrac_phy, ONLY: ntiso,nqtot,ivap,iliq,isol,iqWIsoPha
 #ifdef ISOVERIF
     USE isotopes_verif_mod
@@ -23 +19 @@
 
     INTEGER klon,klev,iflag_ice_thermo
-    REAL, DIMENSION(klon,klev), INTENT(in) :: t_seri,q_seri,ql_seri,qs_seri
-    REAL, DIMENSION(klon,klev), INTENT(out) :: d_t_eva,d_q_eva,d_ql_eva,d_qs_eva
+    REAL, DIMENSION(klon,klev), INTENT(in) :: t_seri
+    REAL, DIMENSION(klon,klev,nqtot), INTENT(in) ::     qx
+    REAL, DIMENSION(klon,klev), INTENT(out) :: d_t_eva
+    REAL, DIMENSION(klon,klev,nqtot), INTENT(out) ::        d_qx_eva
 
     REAL za,zb,zdelta,zlvdcp,zlsdcp
-    INTEGER i,k
+    INTEGER i,k,ixt,ivapcur,iliqcur,isolcur    
 
-#ifdef ISO
-    REAL, DIMENSION(ntiso,klon,klev), INTENT(in) :: xt_seri,xtl_seri,xts_seri
-    REAL, DIMENSION(ntiso,klon,klev), INTENT(out) :: d_xt_eva,d_xtl_eva,d_xts_eva
-    integer ixt
-#endif
 
     !--------Stochastic Boundary Layer Triggering: ALE_BL--------
@@ -42 +35 @@
     !IM 100106 BEG : pouvoir sortir les ctes de la physique
     !
+do ixt=1,1+ntiso
     ! Re-evaporer l'eau liquide nuageuse
     !
+    iliqcur= iqWIsoPha(ixt,iliq)    
+    ivapcur= iqWIsoPha(ixt,ivap)    
+    isolcur= iqWIsoPha(ixt,isol)   
 !print *,'rrevap ; fl_cor_ebil:',fl_cor_ebil,' iflag_ice_thermo:',iflag_ice_thermo,' RVTMP2',RVTMP2
     DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
        DO i = 1, klon
-        if (fl_cor_ebil .GT. 0) then
-          zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*(q_seri(i,k)+ql_seri(i,k)+qs_seri(i,k)))
-          zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*(q_seri(i,k)+ql_seri(i,k)+qs_seri(i,k)))
-        else
-          zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
-          !jyg<
-          !  Attention : Arnaud a propose des formules completement differentes
-          !                  A verifier !!!
-          zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
-        end if
-          IF (iflag_ice_thermo .EQ. 0) THEN
+
+        if (ixt.eq.1) then
+         if (fl_cor_ebil .GT. 0) then
+          !zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*(q_seri(i,k)+ql_seri(i,k)+qs_seri(i,k)))
+          !zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*(q_seri(i,k)+ql_seri(i,k)+qs_seri(i,k)))
+          zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*(qx(i,k,ivapcur)+qx(i,k,iliqcur)+qx(i,k,isolcur)))
+          zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*(qx(i,k,ivapcur)+qx(i,k,iliqcur)+qx(i,k,isolcur)))
+         else
+           zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*qx(i,k,ivapcur))
+           !jyg<
+           !  Attention : Arnaud a propose des formules completement differentes
+           !                  A verifier !!!
+           zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*qx(i,k,ivapcur))
+         end if
+         IF (iflag_ice_thermo .EQ. 0) THEN
              zlsdcp=zlvdcp
-          ENDIF
+         ENDIF
           !>jyg
-
-          IF (iflag_ice_thermo.eq.0) THEN   
-             !pas necessaire a priori
-
-             zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
-  zdelta = 0.
-             zb = MAX(0.0,ql_seri(i,k))
-             za = - MAX(0.0,ql_seri(i,k)) &
-                  * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
-             d_t_eva(i,k) = za
-             d_q_eva(i,k) = zb
-             d_ql_eva(i,k) = -ql_seri(i,k)
-             d_qs_eva(i,k) = 0.
-
-#ifdef ISO
-         do ixt=1,ntiso
-            zb = MAX(0.0,xtl_seri(ixt,i,k))
-            d_xt_eva(ixt,i,k) = zb
-            d_xtl_eva(ixt,i,k) = -xtl_seri(ixt,i,k)
-            d_xts_eva(ixt,i,k) = 0.
-         enddo
-#ifdef ISOVERIF
-      do ixt=1,ntiso
-        call iso_verif_noNaN(xt_seri(ixt,i,k), &
-     &     'reevap 2417: apres evap tot')
-      enddo
-      if (iso_eau.gt.0) then
-              call iso_verif_egalite_choix( &
-     &           xt_seri(iso_eau,i,k),q_seri(i,k), &
-     &          'reevap 1891+, après reevap totale',errmax,errmaxrel)
-              call iso_verif_egalite_choix( &
-     &           xtl_seri(iso_eau,i,k),ql_seri(i,k), &
-     &          'reevap 2209+, après reevap totale',errmax,errmaxrel)
-       endif !if (iso_eau.gt.0) then       
-      if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then  
-            if (q_seri(i,k).gt.ridicule) then  
-               if (iso_verif_o18_aberrant_nostop( &
-     &              xt_seri(iso_HDO,i,k)/q_seri(i,k), &
-     &              xt_seri(iso_O18,i,k)/q_seri(i,k), &
-     &              'reevap 2315: apres reevap totale').eq.1) then
-                  write(*,*) 'i,k,q_seri(i,k)=',i,k,q_seri(i,k)
-                  write(*,*) 'd_q_eva(i,k)=',d_q_eva(i,k)
-                  write(*,*) 'deltaD(d_q_eva(i,k))=',deltaD(d_xt_eva(iso_HDO,i,k)/d_q_eva(i,k))
-                  write(*,*) 'deltaO18(d_q_eva(i,k))=',deltaO(d_xt_eva(iso_O18,i,k)/d_q_eva(i,k))
-                  stop
-              endif !  if (iso_verif_o18_aberrant_nostop 
-            endif !if (q_seri(i,k).gt.errmax) then   
-        endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then 
-#ifdef ISOTRAC      
-             call iso_verif_traceur(xt_seri(1,i,k), &
-     &           'reevap 2165a')
-             call iso_verif_traceur_pbidouille(xt_seri(1,i,k), &
-     &           'reevap 2165b')
-#endif               
-         
-#endif           
-#endif
-
-          ELSE
-
              !CR: on r\'e-\'evapore eau liquide et glace
 
@@ -127 +68 @@
              !        za = - MAX(0.0,ql_seri(i,k)) &
              !             * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
-             zb = MAX(0.0,ql_seri(i,k)+qs_seri(i,k))
-             za = - MAX(0.0,ql_seri(i,k))*zlvdcp & 
-                  - MAX(0.0,qs_seri(i,k))*zlsdcp
+             za = - MAX(0.0,qx(i,k,iliqcur))*zlvdcp & 
+                  - MAX(0.0,qx(i,k,iliqcur))*zlsdcp
              d_t_eva(i,k) = za
-             d_q_eva(i,k) = zb
-             d_ql_eva(i,k) = -ql_seri(i,k)
-             d_qs_eva(i,k) = -qs_seri(i,k)
 
-#ifdef ISO
-         do ixt=1,ntiso
-            zb = MAX(0.0,xtl_seri(ixt,i,k)+xts_seri(ixt,i,k))
-            d_xt_eva(ixt,i,k) = zb
-            d_xtl_eva(ixt,i,k) = -xtl_seri(ixt,i,k)
-            d_xts_eva(ixt,i,k) = -xts_seri(ixt,i,k)
-         enddo
+          endif !if (ixt.eq.1) then
 
-#ifdef ISOVERIF
-      do ixt=1,ntiso
-      call iso_verif_noNaN(xt_seri(ixt,i,k), &
-     &     'reevap 2417: apres evap tot')
-      enddo
-      if (iso_eau.gt.0) then
-              call iso_verif_egalite_choix( &
-     &           xt_seri(iso_eau,i,k),q_seri(i,k), &
-     &          'reevap 1891, après réévap totale',errmax,errmaxrel)
-              call iso_verif_egalite_choix( &
-     &           xtl_seri(iso_eau,i,k),ql_seri(i,k), &
-     &          'reevap 2209, après réévap totale',errmax,errmaxrel)
-              call iso_verif_egalite_choix( &
-     &           xts_seri(iso_eau,i,k),qs_seri(i,k), &
-     &          'reevap 2209b, après réévap totale',errmax,errmaxrel)
-       endif !if (iso_eau.gt.0) then
-     
-      if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then  
-            if (q_seri(i,k).gt.ridicule) then  
-               if (iso_verif_o18_aberrant_nostop( &
-     &              xt_seri(iso_HDO,i,k)/q_seri(i,k), &
-     &              xt_seri(iso_O18,i,k)/q_seri(i,k), &
-     &              'reevap 2408: apres reevap totale').eq.1) then
-                  write(*,*) 'i,k,q_seri(i,k)=',i,k,q_seri(i,k)                       
-                  stop
-              endif !  if (iso_verif_o18_aberrant_nostop 
-            endif !if (q_seri(i,k).gt.errmax) then   
-        endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then        
-#ifdef ISOTRAC      
-             call iso_verif_traceur(xt_seri(1,i,k), &
-     &           'reevap 2165c')
-             call iso_verif_traceur_pbidouille(xt_seri(1,i,k), &
-     &           'reevap 2165d')
-#endif                  
-#endif                  
-#endif
+            !zb = MAX(0.0,ql_seri(i,k)+qs_seri(i,k))
+            !d_q_eva(i,k) = zb
+            !d_ql_eva(i,k) = -ql_seri(i,k)
+            !d_qs_eva(i,k) = -qs_seri(i,k)
 
-          ENDIF
+            zb = MAX(0.0,qx(i,k,iliqcur)+qx(i,k,isolcur))
+            d_qx_eva(i,k,ivapcur) = zb
+            d_qx_eva(i,k,iliqcur) = -qx(i,k,iliqcur)
+            d_qx_eva(i,k,isolcur) = -qx(i,k,isolcur)
+
 
        ENDDO
     ENDDO
 
+    enddo ! do ixt=1,1+niso*(nzone +1)
+    
+
 RETURN