Index: DZ6/trunk/libf/phylmd/calcul_divers.h
===================================================================
--- /LMDZ6/trunk/libf/phylmd/calcul_divers.h	(revision 5292)
+++ 	(revision )
@@ -1,89 +1,0 @@
-!
-! $Id$
-!
-! itap : nombre de pas de temps de la physique
-! itapm1 : somme du nombre de pas de temps du/des mois precedent/s.
-! (itap - itapm1) : nombre de pas de temps du mois courant
-!
-! Ne pas modifier les IFs ci-dessous impliquant itapm1,
-! autrement les resultats seront faux !!
-!
-! Ici on utilise MOD(itap - itapm1,NINT(mth_len*un_jour/phys_tstep)).EQ.1)
-! pour detecter le debut de chaque mois lorsque  l on tourne par an.
-!
-! IM, 26.05.2023
-!
-!
-! Initialisations itapm1 du premier mois
-      IF(itap.EQ.1) THEN
-         itapm1=0
-!        print*,'initialisation itap=1 itapm1 ',itapm1
-      ENDIF
-
-!
-! Initialisation debut de mois
-      IF(itap.EQ.itapm1+1) THEN
-         ndayrain_mth(:)=0.
-!        print*,'Initialisation ndayrain_mth ',itap
-      ENDIF
-!
-! Initialisation debut de chaque jour
-      IF(MOD(itap,NINT(un_jour/phys_tstep)).EQ.1) THEN
-        nday_rain(:)=0.
-!       print*,'initialisation mois suivants day_rain itap',itap
-      ENDIF
-!
-! Calcul a chaque pas de temps de la physique
-         DO i = 1, klon
-            total_rain(i)=rain_fall(i)+snow_fall(i)  
-            IF(total_rain(i).GT.0.) nday_rain(i)=1.
-         ENDDO
-!
-! Cumul en fin de journee
-      IF(MOD(itap,NINT(un_jour/phys_tstep)).EQ.0) THEN
-         DO i = 1, klon
-            ndayrain_mth(i)=ndayrain_mth(i)+nday_rain(i)
-         ENDDO
-      ENDIF
-!
-! Initialisation fin de mois
-! Ne pas changer le IF ci-dessous, car le compteur itapm1 est augmente
-! apres, dans la boucle !!!
-! IM, 260523
-      IF(MOD(itap-itapm1,NINT(mth_len*un_jour/phys_tstep)).EQ.0) THEN
-        itapm1=itapm1+NINT(mth_len*un_jour/phys_tstep)
-!       print*,'fmois i i1 mjt',itap,itapm1,mth_len,un_jour,phys_tstep
-      ENDIF
-!
-! calcul temperatures minimale et maximale moyennees sur le mois 
-!
-!initialisation debut de mois pour les fichiers mensuels annuels
-  IF(itap.EQ.itapm1+1) THEN
-     t2m_min_mon=0.
-     t2m_max_mon=0.
-  ENDIF
-!
-!initialisation debut de journee pour les fichiers mensuels annuels
-  IF(MOD(itap,NINT(un_jour/phys_tstep)).EQ.1) THEN
-     zt2m_min_mon=zt2m
-     zt2m_max_mon=zt2m
-  ENDIF
-!
-!calcul sur tous les pas de temps pour les fichiers mensuels annuels
-     DO i = 1, klon
-        zt2m_min_mon(i)=MIN(zt2m(i),zt2m_min_mon(i))
-        zt2m_max_mon(i)=MAX(zt2m(i),zt2m_max_mon(i))
-     ENDDO
-!
-!fin journee
-  IF(MOD(itap,NINT(un_jour/phys_tstep)).EQ.0) THEN
-   t2m_min_mon=t2m_min_mon+zt2m_min_mon 
-   t2m_max_mon=t2m_max_mon+zt2m_max_mon 
-  ENDIF
-!
-!fin mois
-  IF(itap==itapm1) THEN
-   t2m_min_mon=t2m_min_mon/mth_len
-   t2m_max_mon=t2m_max_mon/mth_len
-  ENDIF
-!
Index: /LMDZ6/trunk/libf/phylmd/calcul_divers_mod_h.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/calcul_divers_mod_h.f90	(revision 5293)
+++ /LMDZ6/trunk/libf/phylmd/calcul_divers_mod_h.f90	(revision 5293)
@@ -0,0 +1,105 @@
+MODULE calcul_divers_mod_h
+  IMPLICIT NONE; PRIVATE
+  PUBLIC calcul_divers
+
+CONTAINS
+
+  SUBROUTINE calcul_divers(itap, itapm1, un_jour)
+    ! IM, 26.05.2023
+    ! Ne pas modifier les IFs ci-dessous impliquant itapm1,
+    ! autrement les resultats seront faux !!
+
+    ! On utilise MOD(itap - itapm1,NINT(mth_len*un_jour/phys_tstep)).EQ.1)
+    ! pour detecter le debut de chaque mois lorsque l'on tourne par an.
+
+    ! itap : nombre de pas de temps de la physique
+    ! itapm1 : somme du nombre de pas de temps du/des mois precedent/s.
+    ! (itap - itapm1) : nombre de pas de temps du mois courant
+    USE dimphy, ONLY: klon
+    USE phys_state_var_mod, ONLY: phys_tstep, ndayrain_mth, nday_rain, total_rain, rain_fall, snow_fall
+    USE phys_local_var_mod, ONLY: t2m_min_mon, t2m_max_mon, zt2m_min_mon, zt2m_max_mon, zt2m
+    USE phys_cal_mod, ONLY: mth_len
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(IN) :: itap
+    INTEGER, INTENT(INOUT) :: itapm1
+    REAL, INTENT(IN) :: un_jour
+
+    INTEGER :: i
+
+    ! Initialisations itapm1 du premier mois
+    IF(itap==1) THEN
+      itapm1 = 0
+      !        PRINT*,'initialisation itap=1 itapm1 ',itapm1
+    ENDIF
+
+    ! Initialisation debut de mois
+    IF(itap==itapm1 + 1) THEN
+      ndayrain_mth(:) = 0.
+      !        PRINT*,'Initialisation ndayrain_mth ',itap
+    ENDIF
+
+    ! Initialisation debut de chaque jour
+    IF(MOD(itap, NINT(un_jour / phys_tstep))==1) THEN
+      nday_rain(:) = 0.
+      !       PRINT*,'initialisation mois suivants day_rain itap',itap
+    ENDIF
+
+    ! Calcul a chaque pas de temps de la physique
+    DO i = 1, klon
+      total_rain(i) = rain_fall(i) + snow_fall(i)
+      IF(total_rain(i)>0.) nday_rain(i) = 1.
+    ENDDO
+
+    ! Cumul en fin de journee
+    IF(MOD(itap, NINT(un_jour / phys_tstep))==0) THEN
+      DO i = 1, klon
+        ndayrain_mth(i) = ndayrain_mth(i) + nday_rain(i)
+      ENDDO
+    ENDIF
+
+    ! Initialisation fin de mois
+    ! Ne pas changer le IF ci-dessous, car le compteur itapm1 est augmente
+    ! apres, dans la boucle !!!
+    ! IM, 260523
+    IF(MOD(itap - itapm1, NINT(mth_len * un_jour / phys_tstep))==0) THEN
+      itapm1 = itapm1 + NINT(mth_len * un_jour / phys_tstep)
+      !       PRINT*,'fmois i i1 mjt',itap,itapm1,mth_len,un_jour,phys_tstep
+    ENDIF
+
+    ! calcul temperatures minimale et maximale moyennees sur le mois
+
+    !initialisation debut de mois pour les fichiers mensuels annuels
+    IF(itap==itapm1 + 1) THEN
+      t2m_min_mon = 0.
+      t2m_max_mon = 0.
+    ENDIF
+
+    !initialisation debut de journee pour les fichiers mensuels annuels
+    IF(MOD(itap, NINT(un_jour / phys_tstep))==1) THEN
+      zt2m_min_mon = zt2m
+      zt2m_max_mon = zt2m
+    ENDIF
+
+    !calcul sur tous les pas de temps pour les fichiers mensuels annuels
+    DO i = 1, klon
+      zt2m_min_mon(i) = MIN(zt2m(i), zt2m_min_mon(i))
+      zt2m_max_mon(i) = MAX(zt2m(i), zt2m_max_mon(i))
+    ENDDO
+
+    !fin journee
+    IF(MOD(itap, NINT(un_jour / phys_tstep))==0) THEN
+      t2m_min_mon = t2m_min_mon + zt2m_min_mon
+      t2m_max_mon = t2m_max_mon + zt2m_max_mon
+    ENDIF
+
+    !fin mois
+    IF(itap==itapm1) THEN
+      t2m_min_mon = t2m_min_mon / mth_len
+      t2m_max_mon = t2m_max_mon / mth_len
+    ENDIF
+
+  END SUBROUTINE calcul_divers
+
+END MODULE calcul_divers_mod_h
Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 5292)
+++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 5293)
@@ -18,5 +18,5 @@
 ! For clarity, the "USE" section is now arranged in alphabetical order,
 ! with a separate section for CPP keys
-! PLEASE try to follow this rule 
+! PLEASE try to follow this rule
 
     USE ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
@@ -84,5 +84,5 @@
     USE lmdz_thermcell_ini, ONLY : thermcell_ini, iflag_thermals_tenv
     USE lmdz_thermcell_dtke, ONLY : thermcell_dtke
-    USE lmdz_blowing_snow_ini, ONLY : blowing_snow_ini , qbst_bs 
+    USE lmdz_blowing_snow_ini, ONLY : blowing_snow_ini , qbst_bs
     USE lmdz_lscp_ini, ONLY : lscp_ini
     USE lmdz_ratqs_main, ONLY : ratqs_main
@@ -143,5 +143,5 @@
        t_seri,q_seri,ql_seri,qs_seri,qbs_seri, &
        u_seri,v_seri,cf_seri,rvc_seri,tr_seri, &
-       rhcl, &        
+       rhcl, &
        ! Dynamic tendencies (diagnostics)
        d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_qbs_dyn, &
@@ -219,5 +219,5 @@
        topswai_aerop, solswai_aerop,   &
        topswad0_aerop, solswad0_aerop, &
-       topsw_aerop, topsw0_aerop,      & 
+       topsw_aerop, topsw0_aerop,      &
        solsw_aerop, solsw0_aerop,      &
        topswcf_aerop, solswcf_aerop,   &
@@ -250,6 +250,4 @@
        zustar, zu10m, zv10m, rh2m, qsat2m, &
        zq2m, zt2m, zn2mout, weak_inversion, &
-       zt2m_min_mon, zt2m_max_mon,   &         ! pour calcul_divers.h
-       t2m_min_mon, t2m_max_mon,  &            ! pour calcul_divers.h
        !
        s_pblh_x, s_pblh_w, &
@@ -275,5 +273,5 @@
        !
        wake_k, &
-       alp_wake, & 
+       alp_wake, &
        wake_h, wake_omg, &
                        ! tendencies of delta T and delta q:
@@ -287,5 +285,5 @@
 !!!       d_s_vdf, d_dens_a_vdf, d_dens_vdf, & ! due to vertical diffusion
 !!!       d_s_the, d_dens_a_the, d_dens_the, & ! due to thermals
-       !                                  
+       !
        ptconv, ratqsc, &
        wbeff, convoccur, zmax_th, &
@@ -338,5 +336,5 @@
        t2m, fluxlat,  &
        fsollw, evap_pot,  &
-       fsolsw, wfbils, wfevap, &  
+       fsolsw, wfbils, wfevap, &
        prfl, psfl,bsfl, fraca, Vprecip,  &
        zw2,  &
@@ -348,9 +346,9 @@
        wwriteSTD, phiwriteSTD, &              !pour calcul_STDlev.h
        qwriteSTD, twriteSTD, rhwriteSTD, &    !pour calcul_STDlev.h
-       ! 
+       !
        beta_prec,  &
        rneb,  &
        zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic, &
-       zxfluxt,zxfluxq 
+       zxfluxt,zxfluxq
        !
        USE phys_local_var_mod, ONLY: zfice, dNovrN, ptconv
@@ -363,4 +361,5 @@
        USE alpale_mod
        USE yoethf_mod_h
+       USE calcul_divers_mod_h, ONLY: calcul_divers
 
     IMPLICIT NONE
@@ -373,5 +372,5 @@
     !!AA                  -  uniformisation des parametrisations ds phytrac
     !!AA                  -  stockage des moyennes des champs necessaires
-    !!AA                     en mode traceur off-line 
+    !!AA                     en mode traceur off-line
     !!======================================================================
     !!   CLEFS CPP POUR LES IO
@@ -439,5 +438,5 @@
     ! Clef iflag_cycle_diurne controlant l'activation du cycle diurne:
     ! en attente du codage des cles par Fred
-    ! iflag_cycle_diurne est initialise par conf_phys et se trouve 
+    ! iflag_cycle_diurne est initialise par conf_phys et se trouve
     ! dans clesphys.h (IM)
     !======================================================================
@@ -472,11 +471,11 @@
     !$OMP THREADPRIVATE(ok_instan)
     !
-    LOGICAL ok_LES ! sortir le fichier LES 
-    SAVE ok_LES                            
-    !$OMP THREADPRIVATE(ok_LES)                  
-    !
-    LOGICAL callstats ! sortir le fichier stats 
-    SAVE callstats                            
-    !$OMP THREADPRIVATE(callstats)                  
+    LOGICAL ok_LES ! sortir le fichier LES
+    SAVE ok_LES
+    !$OMP THREADPRIVATE(ok_LES)
+    !
+    LOGICAL callstats ! sortir le fichier stats
+    SAVE callstats
+    !$OMP THREADPRIVATE(callstats)
     !
     LOGICAL ok_region ! sortir le fichier regional
@@ -486,7 +485,7 @@
     SAVE seuil_inversion
     !$OMP THREADPRIVATE(seuil_inversion)
-    
-    
-    
+
+
+
     real facteur
 
@@ -551,5 +550,5 @@
 !!    real wght_cvfd(klon,klev)
 !!    ! Variables pour le lessivage convectif
-!!    ! RomP >>> 
+!!    ! RomP >>>
 !!    real phi2(klon,klev,klev)
 !!    real d1a(klon,klev),dam(klon,klev)
@@ -572,8 +571,8 @@
     INTEGER n
     !ym      INTEGER npoints
-    !ym      PARAMETER(npoints=klon) 
+    !ym      PARAMETER(npoints=klon)
     !
     INTEGER nregISCtot
-    PARAMETER(nregISCtot=1) 
+    PARAMETER(nregISCtot=1)
     !
     ! imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties
@@ -584,5 +583,5 @@
     ! direction j (latitude)
 !JLD    INTEGER imin_debut, nbpti
-!JLD    INTEGER jmin_debut, nbptj 
+!JLD    INTEGER jmin_debut, nbptj
     !IM: region='3d' <==> sorties en global
     CHARACTER*3 region
@@ -673,5 +672,5 @@
     INTEGER,  SAVE               :: iflag_alp_wk_cond=0 ! wake: if =0, then Alp_wk is the average lifting
                                                         ! power provided by the wakes; else, Alp_wk is the
-                                                        ! lifting power conditionned on the presence of a 
+                                                        ! lifting power conditionned on the presence of a
                                                         ! gust-front in the grid cell.
     !$OMP THREADPRIVATE(iflag_alp_wk_cond)
@@ -696,8 +695,8 @@
     !
 !!!    INTEGER, SAVE, DIMENSION(klon)   :: wake_k
-!!!    !$OMP THREADPRIVATE(wake_k) 
+!!!    !$OMP THREADPRIVATE(wake_k)
     !
     !jyg<
-    !cc      REAL wake_pe(klon)              ! Wake potential energy - WAPE 
+    !cc      REAL wake_pe(klon)              ! Wake potential energy - WAPE
     !>jyg
 
@@ -716,9 +715,9 @@
     REAL d_q_adjwk(klon,klev)                !jyg
     LOGICAL,SAVE :: ok_adjwk=.FALSE.
-    !$OMP THREADPRIVATE(ok_adjwk) 
+    !$OMP THREADPRIVATE(ok_adjwk)
     INTEGER,SAVE :: iflag_adjwk=0            !jyg
     !$OMP THREADPRIVATE(iflag_adjwk)         !jyg
     REAL,SAVE :: oliqmax=999.,oicemax=999.
-    !$OMP THREADPRIVATE(oliqmax,oicemax) 
+    !$OMP THREADPRIVATE(oliqmax,oicemax)
     REAL, SAVE :: alp_offset
     !$OMP THREADPRIVATE(alp_offset)
@@ -728,5 +727,5 @@
     !$OMP THREADPRIVATE(dqcon_multistep_max)
 
-  
+
     !
     !RR:fin declarations poches froides
@@ -735,5 +734,5 @@
     REAL ztv(klon,klev),ztva(klon,klev)
     REAL zpspsk(klon,klev)
-    REAL ztla(klon,klev),zqla(klon,klev) 
+    REAL ztla(klon,klev),zqla(klon,klev)
     REAL zthl(klon,klev)
 
@@ -741,5 +740,5 @@
 
     !--------Stochastic Boundary Layer Triggering: ALE_BL--------
-    !---Propri\'et\'es du thermiques au LCL 
+    !---Propri\'et\'es du thermiques au LCL
 !    real zlcl_th(klon)          ! Altitude du LCL calcul\'e
     ! continument (pcon dans
@@ -749,6 +748,6 @@
     real w_conv(klon)           ! Vitesse verticale de grande \'echelle au LCL
     real tke0(klon,klev+1)      ! TKE au d\'ebut du pas de temps
-    real therm_tke_max0(klon)   ! TKE dans les thermiques au LCL 
-    real env_tke_max0(klon)     ! TKE dans l'environnement au LCL 
+    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)
@@ -758,5 +757,5 @@
 
     REAL,SAVE :: random_notrig_max=1.
-    !$OMP THREADPRIVATE(random_notrig_max) 
+    !$OMP THREADPRIVATE(random_notrig_max)
 
     !--------Statistical Boundary Layer Closure: ALP_BL--------
@@ -767,9 +766,9 @@
     !-------Activer les tendances de TKE due a l'orograp??ie---------
      INTEGER, SAVE :: addtkeoro
-    !$OMP THREADPRIVATE(addtkeoro) 
+    !$OMP THREADPRIVATE(addtkeoro)
      REAL, SAVE :: alphatkeoro
-    !$OMP THREADPRIVATE(alphatkeoro) 
+    !$OMP THREADPRIVATE(alphatkeoro)
      LOGICAL, SAVE :: smallscales_tkeoro
-    !$OMP THREADPRIVATE(smallscales_tkeoro) 
+    !$OMP THREADPRIVATE(smallscales_tkeoro)
 
 
@@ -786,5 +785,5 @@
     !
     !AA
-    !AA  Pour phytrac 
+    !AA  Pour phytrac
     REAL u1(klon)             ! vents dans la premiere couche U
     REAL v1(klon)             ! vents dans la premiere couche V
@@ -795,5 +794,5 @@
     REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
     REAL frac_nucl(klon,klev) ! idem (nucleation)
-    ! RomP >>> 
+    ! RomP >>>
     REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
     ! RomP <<<
@@ -820,6 +819,6 @@
     INTEGER lmt_pas
     SAVE lmt_pas                ! frequence de mise a jour
-    !$OMP THREADPRIVATE(lmt_pas) 
-    real zmasse(klon, nbp_lev),exner(klon, nbp_lev) 
+    !$OMP THREADPRIVATE(lmt_pas)
+    real zmasse(klon, nbp_lev),exner(klon, nbp_lev)
     !     (column-density of mass of air in a cell, in kg m-2)
     real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
@@ -866,5 +865,5 @@
     REAL radocond(klon,klev)  ! eau condensee nuageuse
     !
-    !XXX PB 
+    !XXX PB
     REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
     REAL fluxqbs(klon,klev, nbsrf)   ! flux turbulent de neige soufflee
@@ -919,7 +918,7 @@
     ! La convection n'est pas calculee tous les pas, il faut donc
     !                      sauvegarder les sorties de la convection
-    !ym      SAVE  
-    !ym      SAVE  
-    !ym      SAVE  
+    !ym      SAVE
+    !ym      SAVE
+    !ym      SAVE
     !
     INTEGER itapcv, itapwk
@@ -961,8 +960,8 @@
     REAL, DIMENSION(klon,klev)  :: d_q_ch4_dtime
     !
-    ! Flag pour pouvoir ne pas ajouter les tendances. 
+    ! Flag pour pouvoir ne pas ajouter les tendances.
     ! Par defaut, les tendances doivente etre ajoutees et
     ! flag_inhib_tend = 0
-    ! flag_inhib_tend > 0 : tendances non ajoutees, avec un nombre 
+    ! flag_inhib_tend > 0 : tendances non ajoutees, avec un nombre
     ! croissant de print quand la valeur du flag augmente
     !!! attention, ce flag doit etre change avec prudence !!!
@@ -970,5 +969,5 @@
 !!    INTEGER :: flag_inhib_tend = 2
     !
-    ! Logical switch to a bug : reseting to 0 convective variables at the 
+    ! Logical switch to a bug : reseting to 0 convective variables at the
     ! begining of physiq.
     LOGICAL, SAVE :: ok_bug_cv_trac = .TRUE.
@@ -980,5 +979,5 @@
     !$OMP THREADPRIVATE(ok_bug_split_th)
 
-    ! Logical switch to a bug : modifying directly wake_deltat  by adding 
+    ! 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.
@@ -1061,5 +1060,5 @@
     !
     REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
-    REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D 
+    REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
 !JLD    REAL zx_tmp_2d(nbp_lon,nbp_lat)
 !JLD    REAL zx_lon(nbp_lon,nbp_lat)
@@ -1096,5 +1095,5 @@
     ! essai writephys
     INTEGER fid_day, fid_mth, fid_ins
-    PARAMETER (fid_ins = 1, fid_day = 2, fid_mth = 3) 
+    PARAMETER (fid_ins = 1, fid_day = 2, fid_mth = 3)
     INTEGER prof2d_on, prof3d_on, prof2d_av, prof3d_av
     PARAMETER (prof2d_on = 1, prof3d_on = 2, prof2d_av = 3, prof3d_av = 4)
@@ -1111,5 +1110,5 @@
     INTEGER :: naero
     ! Aerosol optical properties
-    CHARACTER*4, DIMENSION(naero_grp) :: rfname 
+    CHARACTER*4, DIMENSION(naero_grp) :: rfname
     REAL, DIMENSION(klon,klev)     :: mass_solu_aero ! total mass
     ! concentration
@@ -1129,10 +1128,10 @@
     LOGICAL, SAVE :: aerosol_couple ! true  : calcul des aerosols dans INCA
     ! false : lecture des aerosol dans un fichier
-    !$OMP THREADPRIVATE(aerosol_couple)    
-    LOGICAL, SAVE :: chemistry_couple ! true  : use INCA chemistry O3 
+    !$OMP THREADPRIVATE(aerosol_couple)
+    LOGICAL, SAVE :: chemistry_couple ! true  : use INCA chemistry O3
     ! false : use offline chemistry O3
-    !$OMP THREADPRIVATE(chemistry_couple)    
-    INTEGER, SAVE :: flag_aerosol 
-    !$OMP THREADPRIVATE(flag_aerosol) 
+    !$OMP THREADPRIVATE(chemistry_couple)
+    INTEGER, SAVE :: flag_aerosol
+    !$OMP THREADPRIVATE(flag_aerosol)
     LOGICAL, SAVE :: flag_bc_internal_mixture
     !$OMP THREADPRIVATE(flag_bc_internal_mixture)
@@ -1191,5 +1190,5 @@
     INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
     REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
-    !$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP) 
+    !$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP)
     INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
     REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
@@ -1230,5 +1229,5 @@
 
     !lwoff=y : offset LW CRE for radiation code and other schemes
-    REAL, SAVE :: betalwoff 
+    REAL, SAVE :: betalwoff
     !OMP THREADPRIVATE(betalwoff)
 !
@@ -1245,5 +1244,5 @@
     !albedo SB <<<
 
-    !--Lea Raillard qs_ini 
+    !--Lea Raillard qs_ini
     REAL, dimension(klon,klev) :: qs_ini
 
@@ -1291,5 +1290,5 @@
     ! set-up call to alerte function
     call_alert = (alert_first_call .AND. is_master)
-    
+
     ! Ehouarn: set value of jjmp1 since it is no longer a "fixed parameter"
     jjmp1=nbp_lat
@@ -1304,5 +1303,5 @@
 
     IF (using_xios) THEN
-      ! switch to XIOS LMDZ physics context 
+      ! switch to XIOS LMDZ physics context
       IF (.NOT. debut .AND. is_omp_master) THEN
         CALL wxios_set_context()
@@ -1324,5 +1323,5 @@
             'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
        write(lunout,*) &
-            nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys 
+            nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
 
        write(lunout,*) 'paprs, play, phi, u, v, t'
@@ -1341,5 +1340,5 @@
             "physiq_mod paprs bad order")
 
-    IF (first) THEN 
+    IF (first) THEN
        ivap = strIdx(tracers(:)%name, addPhase('H2O', 'g'))
        iliq = strIdx(tracers(:)%name, addPhase('H2O', 'l'))
@@ -1376,5 +1375,5 @@
        CALL phys_state_var_init(read_climoz)
        CALL phys_output_var_init
-       IF (read_climoz>=1 .AND. create_etat0_limit .AND. grid_type==unstructured) & 
+       IF (read_climoz>=1 .AND. create_etat0_limit .AND. grid_type==unstructured) &
           CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
 
@@ -1437,5 +1436,5 @@
     ENDIF  ! first
 
-    !ym => necessaire pour iflag_con != 2    
+    !ym => necessaire pour iflag_con != 2
     pmfd(:,:) = 0.
     pen_u(:,:) = 0.
@@ -1465,7 +1464,7 @@
        call getin_p('iflag_thermcell_tke', iflag_thermcell_tke)                          !
 
-       CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond) 
+       CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond)
        CALL getin_p('random_notrig_max',random_notrig_max)
-       CALL getin_p('ok_adjwk',ok_adjwk) 
+       CALL getin_p('ok_adjwk',ok_adjwk)
        IF (ok_adjwk) iflag_adjwk=2  ! for compatibility with older versions
        ! iflag_adjwk: ! 0 = Default: no convective adjustment of w-region
@@ -1481,5 +1480,5 @@
        iflag_wake_tend = 0
        CALL getin_p('iflag_wake_tend',iflag_wake_tend)
-       ok_bad_ecmwf_thermo=.TRUE. ! By default thermodynamical constants are set 
+       ok_bad_ecmwf_thermo=.TRUE. ! By default thermodynamical constants are set
                                   ! in rrtm/suphec.F90 (and rvtmp2 is set to 0).
        CALL getin_p('ok_bad_ecmwf_thermo',ok_bad_ecmwf_thermo)
@@ -1509,11 +1508,11 @@
        WRITE(lunout,*) 'ok_adjwk=',           ok_adjwk
        WRITE(lunout,*) 'iflag_adjwk=',        iflag_adjwk
-       WRITE(lunout,*) 'qtcon_multistep_max=',dtcon_multistep_max 
-       WRITE(lunout,*) 'qdcon_multistep_max=',dqcon_multistep_max 
+       WRITE(lunout,*) 'qtcon_multistep_max=',dtcon_multistep_max
+       WRITE(lunout,*) 'qdcon_multistep_max=',dqcon_multistep_max
        WRITE(lunout,*) 'ratqsp0=',            ratqsp0
-       WRITE(lunout,*) 'ratqsdp=',            ratqsdp 
+       WRITE(lunout,*) 'ratqsdp=',            ratqsdp
        WRITE(lunout,*) 'iflag_wake_tend=',    iflag_wake_tend
-       WRITE(lunout,*) 'ok_bad_ecmwf_thermo=',ok_bad_ecmwf_thermo 
-       WRITE(lunout,*) 'ok_bug_cv_trac=',     ok_bug_cv_trac 
+       WRITE(lunout,*) 'ok_bad_ecmwf_thermo=',ok_bad_ecmwf_thermo
+       WRITE(lunout,*) 'ok_bug_cv_trac=',     ok_bug_cv_trac
        WRITE(lunout,*) 'ok_bug_split_th=',    ok_bug_split_th
        WRITE(lunout,*) 'fl_ebil=',            fl_ebil
@@ -1552,5 +1551,5 @@
        !des caracteristiques du thermique
        wght_th(:,:)=1.
-       lalim_conv(:)=1 
+       lalim_conv(:)=1
        !RC
        ustar(:,:)=0.
@@ -1578,5 +1577,5 @@
           CALL getin_p('config_inca',config_inca)
 
-       ELSE 
+       ELSE
           config_inca='none' ! default
        ENDIF
@@ -1623,5 +1622,5 @@
 
        IF (iflag_pbl>1) THEN
-          PRINT*, "Using method MELLOR&YAMADA" 
+          PRINT*, "Using method MELLOR&YAMADA"
        ENDIF
 
@@ -1641,5 +1640,5 @@
        IF (MOD(NINT(86400./phys_tstep),nbapp_rad).EQ.0) THEN
           radpas = NINT( 86400./phys_tstep)/nbapp_rad
-       ELSE 
+       ELSE
           WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
                'multiple de nbapp_rad'
@@ -1657,5 +1656,5 @@
           cvpas = cvpas_0
        print *,'physiq, cvpas ',cvpas
-       ELSE 
+       ELSE
           WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
                'multiple de nbapp_cv'
@@ -1669,5 +1668,5 @@
           wkpas = NINT( 86400./phys_tstep)/nbapp_wk
 !       print *,'physiq, wkpas ',wkpas
-       ELSE 
+       ELSE
           WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
                'multiple de nbapp_wk'
@@ -1705,5 +1704,5 @@
              ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
              !
-             ! lecture des nCFMIP stations CFMIP, de leur numero 
+             ! lecture des nCFMIP stations CFMIP, de leur numero
              ! et des coordonnees geographiques lonCFMIP, latCFMIP
              !
@@ -1761,5 +1760,5 @@
 
        !XXXPB Positionner date0 pour initialisation de ORCHIDEE
-       date0 = jD_ref 
+       date0 = jD_ref
        WRITE(*,*) 'physiq date0 : ',date0
        !
@@ -1811,5 +1810,5 @@
                mr_ozone_cosp0,cldtau_cosp0, cldemi_cosp0)
 END IF
-      ENDIF 
+      ENDIF
 
        !
@@ -1840,5 +1839,5 @@
              IF (ok_cdnc.AND.NRADLP.NE.3) THEN
              abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
-                  // 'pour ok_cdnc' 
+                  // 'pour ok_cdnc'
              CALL abort_physic(modname,abort_message,1)
              ENDIF
@@ -1849,5 +1848,5 @@
 #endif
           ENDIF
-       ENDIF   
+       ENDIF
        CALL cloud_optics_prop_ini(klon, prt_level, lunout, flag_aerosol, &
                                   & ok_cdnc, bl95_b0, &
@@ -1882,5 +1881,5 @@
                                 flag_aerosol, flag_aerosol_strat, ok_cdnc)
 ELSE
-       ! phys_output_write écrit des variables traceurs seulement si iflag_phytrac == 1 
+       ! phys_output_write écrit des variables traceurs seulement si iflag_phytrac == 1
        ! donc seulement dans ce cas on doit appeler phytrac_init()
        IF (iflag_phytrac == 1 ) THEN
@@ -1899,5 +1898,5 @@
          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
@@ -1959,5 +1958,5 @@
        ENDIF
        !
-       IF (phys_tstep*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN 
+       IF (phys_tstep*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
           WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
           WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
@@ -2084,5 +2083,5 @@
        !=============================================================
 
-       IF (using_xios) THEN    
+       IF (using_xios) THEN
          ! Get "missing_val" value from XML files (from temperature variable)
          IF (is_omp_master) CALL xios_get_field_attr("temp",default_value=missing_val)
@@ -2090,5 +2089,5 @@
        ENDIF
 
-       IF (using_xios) THEN    
+       IF (using_xios) THEN
          ! Need to put this initialisation after phyetat0 as in the coupled model the XIOS context is only
          ! initialised at that moment
@@ -2096,32 +2095,32 @@
          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
        ! diagnostics are requested, otherwise there is no point in doing this
          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 
+           !--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 
+              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. 
+             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. & 
+          !--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. &
@@ -2129,5 +2128,5 @@
               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. 
+              ok_4xCO2atm=.TRUE.
            ENDIF
            !$OMP BARRIER
@@ -2311,5 +2310,5 @@
       ENDIF
       WRITE(*,*)'ok_lwoff=',ok_lwoff
-      ! 
+      !
       !lwoff=y to begin only sollw and sollwdown are set up to CS values
       sollw = sollw + betalwoff * (sollw0 - sollw)
@@ -2338,6 +2337,6 @@
     !
     !
-    ! Update fraction of the sub-surfaces (pctsrf) and 
-    ! initialize, where a new fraction has appeared, all variables depending 
+    ! Update fraction of the sub-surfaces (pctsrf) and
+    ! initialize, where a new fraction has appeared, all variables depending
     ! on the surface fraction.
     !
@@ -2390,5 +2389,5 @@
     beta_prec(:,:)=0.
     !
-    !   Output variables from the convective scheme should not be set to 0 
+    !   Output variables from the convective scheme should not be set to 0
     !   since convection is not always called at every time step.
     IF (ok_bug_cv_trac) THEN
@@ -2439,5 +2438,5 @@
     qs_ini(:,:)=qs_seri(:,:)
     !
-    !--OB water mass fixer 
+    !--OB water mass fixer
     IF (ok_water_mass_fixer) THEN
     !--store initial water burden
@@ -2481,5 +2480,5 @@
          IF(.NOT.tracers(iq)%isInPhysics) CYCLE
          itr = itr+1
-         tr_ancien(:,:,itr)=tr_seri(:,:,itr)       
+         tr_ancien(:,:,itr)=tr_seri(:,:,itr)
        enddo
     ENDIF
@@ -2681,5 +2680,5 @@
     day_since_equinox = (jD_cur + jH_cur) - jD_eq
     !
-    !   choix entre calcul de la longitude solaire vraie ou valeur fixee a 
+    !   choix entre calcul de la longitude solaire vraie ou valeur fixee a
     !   solarlong0
     IF (solarlong0<-999.) THEN
@@ -2693,5 +2692,5 @@
     ELSE
        zlongi=solarlong0  ! longitude solaire vraie
-       dist=1.            ! distance au soleil / moyenne 
+       dist=1.            ! distance au soleil / moyenne
     ENDIF
 
@@ -2715,5 +2714,5 @@
        ! recode par Olivier Boucher en sept 2015
        SELECT CASE (iflag_cycle_diurne)
-       CASE(0)  
+       CASE(0)
           !  Sans cycle diurne
           CALL angle(zlongi, latitude_deg, fract, rmu0)
@@ -2721,5 +2720,5 @@
           JrNt = 1.0
           zrmu0 = rmu0
-       CASE(1)  
+       CASE(1)
           !  Avec cycle diurne sans application des poids
           !  bit comparable a l ancienne formulation cycle_diurne=true
@@ -2733,5 +2732,5 @@
           JrNt = 0.0
           WHERE (fract.GT.0.0) JrNt = 1.0
-       CASE(2)  
+       CASE(2)
           !  Avec cycle diurne sans application des poids
           !  On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
@@ -2742,6 +2741,6 @@
           !  premier pas de temps de la physique pendant lequel
           !  itaprad=0
-          zdtime1=phys_tstep*REAL(-MOD(itaprad,radpas)-1)      
-          zdtime2=phys_tstep*REAL(radpas-MOD(itaprad,radpas)-1) 
+          zdtime1=phys_tstep*REAL(-MOD(itaprad,radpas)-1)
+          zdtime2=phys_tstep*REAL(radpas-MOD(itaprad,radpas)-1)
           CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
                latitude_deg,longitude_deg,rmu0,fract)
@@ -2758,5 +2757,5 @@
           ! Calcul du flag jour-nuit
           JrNt = 0.0
-          WHERE (zfract.GT.0.0) JrNt = 1.0 
+          WHERE (zfract.GT.0.0) JrNt = 1.0
        END SELECT
     ENDIF
@@ -2774,6 +2773,6 @@
     ! Cela implique tous les interactions des sous-surfaces et la
     ! partie diffusion turbulent du couche limit.
-    ! 
-    ! Certains varibales de sorties de pbl_surface sont utiliser que pour 
+    !
+    ! Certains varibales de sorties de pbl_surface sont utiliser que pour
     ! ecriture des fihiers hist_XXXX.nc, ces sont :
     !   qsol,      zq2m,      s_pblh,  s_lcl,
@@ -2787,5 +2786,5 @@
     !   wfbils,    fluxt,   fluxu, fluxv,
     !
-    ! Certains ne sont pas utiliser du tout : 
+    ! Certains ne sont pas utiliser du tout :
     !   dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
     !
@@ -2828,5 +2827,5 @@
        !ym : Warning gustiness non inialized for iflag_gusts=2 & iflag_gusts=3
        gustiness=0  !ym missing init
-       
+
        IF (iflag_gusts==0) THEN
           gustiness(1:klon)=0
@@ -2875,5 +2874,5 @@
                                 !albedo SB >>>
                                 ! albsol1,   albsol2,   sens,    evap,      &
-            albsol_dir,   albsol_dif,   sens,    evap, snowerosion, &  
+            albsol_dir,   albsol_dif,   sens,    evap, snowerosion, &
                                 !albedo SB <<<
             albsol3_lic,runoff,   snowhgt,   qsnow, to_ice, sissnow, &
@@ -2902,5 +2901,5 @@
             z0m, z0h,     agesno,    fsollw,  fsolsw, &
             d_ts,      fevap,     fluxlat, t2m, &
-            wfbils, wfevap, & 
+            wfbils, wfevap, &
             fluxt,   fluxu,  fluxv, &
             dsens,     devap,     zxsnow, &
@@ -2919,5 +2918,5 @@
           d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
           CALL add_wake_tend &
-             (d_deltat_vdf, d_deltaq_vdf, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy) 
+             (d_deltat_vdf, d_deltaq_vdf, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy)
        ELSE
           d_deltat_vdf(:,:) = 0.
@@ -3164,5 +3163,5 @@
              IF (iflag_adjwk == 2 .AND. OK_bug_ajs_cv) THEN
                CALL add_wake_tend &
-                 (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy) 
+                 (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy)
              ENDIF  ! (iflag_adjwk == 2 .AND. OK_bug_ajs_cv)
           ENDIF  ! (iflag_adjwk >= 1)
@@ -3170,5 +3169,5 @@
        !>jyg
        !
-        
+
 !!      print *,'physiq. q_w(1,k), q_x(1,k) ', &
 !!             (k, q_w(1,k), q_x(1,k),k=1,25)
@@ -3248,5 +3247,5 @@
           !
           !jyg<
-          ! If convective tendencies are too large, then call convection 
+          ! If convective tendencies are too large, then call convection
           !  every time step
           cvpas = cvpas_0
@@ -3411,5 +3410,5 @@
     ENDIF
 
-!!!jyg  Appel diagnostique a add_phys_tend pour tester la conservation de 
+!!!jyg  Appel diagnostique a add_phys_tend pour tester la conservation de
 !!!     l'energie dans les courants satures.
 !!    d_t_con_sat(:,:) = d_t_con(:,:) - ftd(:,:)*dtime
@@ -3417,5 +3416,5 @@
 !!    dql_sat(:,:) = (wdtrainA(:,:)+wdtrainM(:,:))*dtime/zmasse(:,:)
 !!    CALL add_phys_tend(d_u_con, d_v_con, d_t_con_sat, d_q_con_sat, dql_sat,   &
-!!                     dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,&  
+!!                     dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,&
 !!                     itap, 1)
 !!    call prt_enerbil('convection_sat',itap)
@@ -3443,5 +3442,5 @@
     !
     !==========================================================================
-    !RR:Evolution de la poche froide: on ne fait pas de separation wake/env 
+    !RR:Evolution de la poche froide: on ne fait pas de separation wake/env
     !pour la couche limite diffuse pour l instant
     !
@@ -3460,13 +3459,13 @@
           DO k=1,klev
              DO i=1,klon
-                dt_dwn(i,k)  = ftd(i,k) 
-                dq_dwn(i,k)  = fqd(i,k) 
+                dt_dwn(i,k)  = ftd(i,k)
+                dq_dwn(i,k)  = fqd(i,k)
                 M_dwn(i,k)   = dnwd0(i,k)
                 M_up(i,k)    = upwd(i,k)
-                dt_a(i,k)    = d_t_con(i,k)/phys_tstep - ftd(i,k) 
+                dt_a(i,k)    = d_t_con(i,k)/phys_tstep - ftd(i,k)
                 dq_a(i,k)    = d_q_con(i,k)/phys_tstep - fqd(i,k)
              ENDDO
           ENDDO
-         
+
           IF (iflag_wake==2) THEN
              ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
@@ -3493,5 +3492,5 @@
              ENDDO
           ENDIF
-         
+
           !
           !calcul caracteristiques de la poche froide
@@ -3531,5 +3530,5 @@
          CALL add_wake_tend &
             (d_deltat_wk, d_deltaq_wk, d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk, wake_k, &
-             'wake', abortphy) 
+             'wake', abortphy)
           CALL prt_enerbil('wake',itap)
        ENDIF   ! (iflag_wake_tend .GT. 0.)
@@ -3694,11 +3693,11 @@
                 DO i=1,klon
                    !
-                   d_deltat_the(i,k) = - d_t_ajs(i,k) 
-                   d_deltaq_the(i,k) = - d_q_ajs(i,k) 
+                   d_deltat_the(i,k) = - d_t_ajs(i,k)
+                   d_deltaq_the(i,k) = - d_q_ajs(i,k)
                    !
-                   d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i)) 
-                   d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i)) 
-                   d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i)) 
-                   d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i)) 
+                   d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
+                   d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
+                   d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
+                   d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
                    !
                 ENDDO
@@ -3707,8 +3706,8 @@
              IF (ok_bug_split_th) THEN
                CALL add_wake_tend &
-                   (d_deltat_the, d_deltaq_the, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy) 
+                   (d_deltat_the, d_deltaq_the, dsig0, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy)
              ELSE
                CALL add_wake_tend &
-                   (d_deltat_the, d_deltaq_the, dsig0, dsig0, ddens0, ddens0, wake_k, 'the', abortphy) 
+                   (d_deltat_the, d_deltaq_the, dsig0, dsig0, ddens0, ddens0, wake_k, 'the', abortphy)
              ENDIF
              CALL prt_enerbil('the',itap)
@@ -3739,5 +3738,5 @@
              !           zmax_th(i)=pphi(i,lmax_th(i))/rg
              !CR:04/05/12:correction calcul zmax
-             zmax_th(i)=zmax0(i) 
+             zmax_th(i)=zmax0(i)
           ENDDO
 
@@ -3762,5 +3761,5 @@
           ! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
           ! pour des test de convergence numerique.
-          ! Le nouveau ajsec est a priori mieux, meme pour le cas 
+          ! Le nouveau ajsec est a priori mieux, meme pour le cas
           ! iflag_thermals = 0 (l'ancienne version peut faire des tendances
           ! non nulles numeriquement pour des mailles non concernees.
@@ -3791,5 +3790,5 @@
     !
     !===================================================================
-    ! Computation of ratqs, the width (normalized) of the subrid scale 
+    ! Computation of ratqs, the width (normalized) of the subrid scale
     ! water distribution
 
@@ -3836,5 +3835,5 @@
     CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay, &
          t_seri, q_seri,qs_ini,ptconv,ratqs,sigma_qtherm, &
-         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, & 
+         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, &
          pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb,  &
          radocond, picefra, rain_lsc, snow_lsc, &
@@ -3879,5 +3878,5 @@
 !        & ,rain_lsc,snow_lsc
 !    write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
-!-JLD 
+!-JLD
     CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,dqbs0,paprs, &
          'lsc',abortphy,flag_inhib_tend,itap,0)
@@ -3926,5 +3925,5 @@
                 radocond(i,k)=radocond(i,k)+qbs_seri(i,k)
                 picefra(i,k)=(radocond(i,k)*picefra(i,k)+qbs_seri(i,k))/(radocond(i,k))
-                qbsfra=min(qbs_seri(i,k)/qbst_bs,1.0) 
+                qbsfra=min(qbs_seri(i,k)/qbst_bs,1.0)
                 cldfra(i,k)=max(cldfra(i,k),qbsfra)
              ENDDO
@@ -4047,5 +4046,5 @@
              ! profonde.
 
-             !IM/FH: 2011/02/23 
+             !IM/FH: 2011/02/23
              ! definition des points sur lesquels ls thermiques sont actifs
 
@@ -4131,5 +4130,5 @@
     ENDDO
 
-    !IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle 
+    !IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
     !   equivalente a 2m (tpote) pour diagnostique
     !
@@ -4267,5 +4266,5 @@
 
                    !
-                ELSE IF (NSW.EQ.2) THEN 
+                ELSE IF (NSW.EQ.2) THEN
                    !--for now we use the old aerosol properties
                    !
@@ -4307,5 +4306,5 @@
                   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, m_allaer)  
+                  tr_seri, mass_solu_aero, mass_solu_aero_pi, m_allaer)
 #else
                 abort_message='You should compile with -rad ecrad if running with iflag_rrtm=2'
@@ -4314,5 +4313,5 @@
           ENDIF
 
-       ELSE   !--flag_aerosol = 0 
+       ELSE   !--flag_aerosol = 0
           tausum_aero(:,:,:) = 0.
           drytausum_aero(:,:) = 0.
@@ -4368,5 +4367,5 @@
           ENDIF
        ELSE
-          tausum_aero(:,:,id_STRAT_phy) = 0. 
+          tausum_aero(:,:,id_STRAT_phy) = 0.
        ENDIF
 !
@@ -4380,12 +4379,12 @@
 #endif
        !--fin STRAT AEROSOL
-       !     
+       !
 
        ! Calculer les parametres optiques des nuages et quelques
        ! parametres pour diagnostiques:
        !
-       IF (aerosol_couple.AND.config_inca=='aero') THEN 
-          mass_solu_aero(:,:)    = ccm(:,:,1) 
-          mass_solu_aero_pi(:,:) = ccm(:,:,2) 
+       IF (aerosol_couple.AND.config_inca=='aero') THEN
+          mass_solu_aero(:,:)    = ccm(:,:,1)
+          mass_solu_aero_pi(:,:) = ccm(:,:,2)
        ENDIF
 
@@ -4398,5 +4397,5 @@
                cldtaupi, distcltop, temp_cltop, re, fl, ref_liq, ref_ice, &
                ref_liq_pi, ref_ice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
-               reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  & 
+               reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  &
                zfice, dNovrN, ptconv, rnebcon, clwcon)
 
@@ -4465,5 +4464,5 @@
        ENDIF
 
-       !lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek 
+       !lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
        IF (ok_chlorophyll) THEN
           print*,"-- reading chlorophyll"
@@ -4471,9 +4470,9 @@
        ENDIF
 
-!--if ok_suntime_rrtm we use ancillay data for RSUN 
+!--if ok_suntime_rrtm we use ancillay data for RSUN
 !--previous values are therefore overwritten
 !--this is needed for CMIP6 runs
 !--and only possible for new radiation scheme
-       IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN 
+       IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
 #ifdef CPP_RRTM
          CALL read_rsun_rrtm(debut)
@@ -4497,5 +4496,5 @@
        ENDIF
 
-       IF (aerosol_couple.AND.config_inca=='aero') THEN 
+       IF (aerosol_couple.AND.config_inca=='aero') THEN
 IF (CPPKEY_INCA) THEN
           CALL radlwsw_inca  &
@@ -4528,5 +4527,5 @@
           RCFC11 = RCFC11_act
           RCFC12 = RCFC12_act
-          ! 
+          !
           !--interactive CO2 in ppm from carbon cycle
           IF (carbon_cycle_rad) RCO2=RCO2_glo
@@ -4550,7 +4549,7 @@
                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_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
                ! Rajoute par OB pour RRTM
-               tau_aero_lw_rrtm, & 
+               tau_aero_lw_rrtm, &
                cldtaupirad, m_allaer, &
 !              zqsat, flwcrad, fiwcrad, &
@@ -4589,5 +4588,5 @@
                         sollwdown(:))
           cool = cool + betalwoff * (cool0 - cool)
- 
+
           IF (.NOT. using_xios) THEN
             !
@@ -4600,7 +4599,7 @@
                 RN2O_per.NE.RN2O_act.OR. &
                 RCFC11_per.NE.RCFC11_act.OR. &
-                RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE. 
+                RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
           ENDIF
-   ! 
+   !
           IF (ok_4xCO2atm) THEN
                 !
@@ -4621,5 +4620,5 @@
                                 !albedo SB >>>
                                 !      paprs, pplay,zxtsol,albsol1, albsol2,  &
-                     paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, & 
+                     paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
                                 !albedo SB <<<
                      t_seri,q_seri,wo, &
@@ -4665,5 +4664,5 @@
 #ifdef CPP_ECRAD
           IF (ok_3Deffect) then
-!                print*,'ok_3Deffect = ',ok_3Deffect  
+!                print*,'ok_3Deffect = ',ok_3Deffect
                 namelist_ecrad_file='namelist_ecrad_s2'
                 CALL radlwsw &
@@ -4681,5 +4680,5 @@
                      ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
                      namelist_ecrad_file, &
-! A modifier              
+! A modifier
                      heat_s2,heat0_s2,cool_s2,cool0_s2,albpla_s2, &
                      heat_volc,cool_volc, &
@@ -4977,5 +4976,5 @@
     !IM calcul composantes axiales du moment angulaire et couple des montagnes
     !
-    IF (is_sequential .and. ok_orodr) THEN 
+    IF (is_sequential .and. ok_orodr) THEN
        CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
             ra,rg,romega, &
@@ -5057,16 +5056,16 @@
 !------------------
 
-       addtkeoro=0   
-       CALL getin_p('addtkeoro',addtkeoro) 
-      
+       addtkeoro=0
+       CALL getin_p('addtkeoro',addtkeoro)
+
        IF (prt_level.ge.5) &
             print*,'addtkeoro', addtkeoro
-            
-       alphatkeoro=1.   
+
+       alphatkeoro=1.
        CALL getin_p('alphatkeoro',alphatkeoro)
        alphatkeoro=min(max(0.,alphatkeoro),1.)
 
-       smallscales_tkeoro=.FALSE.   
-       CALL getin_p('smallscales_tkeoro',smallscales_tkeoro) 
+       smallscales_tkeoro=.FALSE.
+       CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
 
 
@@ -5076,5 +5075,5 @@
 
 ! Choices for addtkeoro:
-!      ** 0 no TKE tendency from orography    
+!      ** 0 no TKE tendency from orography
 !      ** 1 we include a fraction alphatkeoro of the whole tendency duoro
 !      ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
@@ -5101,5 +5100,5 @@
 ! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
 ! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
-! Mais attention, cela ne va pas dans le sens de la conservation de l'energie! 
+! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
           IF ((zstd(i).GT.1.0) .AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
              itest(i)=1
@@ -5109,5 +5108,5 @@
        ENDDO
 
-     ELSE 
+     ELSE
 
        igwd=0
@@ -5162,5 +5161,5 @@
 
     IF (ok_cosp) THEN
-       ! adeclarer 
+       ! adeclarer
 IF (CPPKEY_COSP) THEN
        IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
@@ -5292,7 +5291,7 @@
     ELSE
        sh_in(:,:) = qx(:,:,ivap)
-       IF (nqo >= 3) THEN 
+       IF (nqo >= 3) THEN
           ch_in(:,:) = qx(:,:,iliq) + qx(:,:,isol)
-       ELSE 
+       ELSE
           ch_in(:,:) = qx(:,:,iliq)
        ENDIF
@@ -5300,5 +5299,5 @@
 
 IF (CPPKEY_DUST) THEN
-    !  Avec SPLA, iflag_phytrac est forcé =1 
+    !  Avec SPLA, iflag_phytrac est forcé =1
     CALL       phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con,       &  ! I
                       pdtphys,ftsol,                                   &  ! I
@@ -5435,5 +5434,5 @@
        ENDDO
     ENDIF
-    
+
     DO i = 1, klon
       !--compute ratio of what q+ql should be with conservation to what it is
@@ -5458,5 +5457,5 @@
 
     !cc prw  = eau precipitable
-    !   prlw = colonne eau liquide 
+    !   prlw = colonne eau liquide
     !   prlw = colonne eau solide
     !   prbsw = colonne neige soufflee
@@ -5491,6 +5490,5 @@
     !
     !IM initialisation + calculs divers diag AMIP2
-    !
-    include "calcul_divers.h"
+    CALL calcul_divers(itap, itapm1, un_jour)
     !
     !IM Interpolation sur les niveaux de pression du NMC
@@ -5678,5 +5676,5 @@
 
     ! Recupere des varibles calcule dans differents modules
-    ! pour ecriture dans histxxx.nc 
+    ! pour ecriture dans histxxx.nc
 
     ! Get some variables from module fonte_neige_mod
@@ -5725,5 +5723,5 @@
     CALL phys_output_write(itap, pdtphys, paprs, pphis,  &
          pplay, lmax_th, aerosol_couple,                 &
-         ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ibs,   & 
+         ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ibs,   &
          ok_sync, ptconv, read_climoz, clevSTD,          &
          ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
@@ -5760,5 +5758,5 @@
     alert_first_call = .FALSE.
 
-    
+
     IF (lafin) THEN
        itau_phy = itau_phy + itap
@@ -5767,7 +5765,7 @@
        !         write(97) u_seri,v_seri,t_seri,q_seri
        !         close(97)
-      
+
        IF (is_omp_master) THEN
-       
+
          IF (read_climoz >= 1) THEN
            IF (is_mpi_root) CALL nf95_close(ncid_climoz)
@@ -5775,5 +5773,5 @@
             DEALLOCATE(press_cen_climoz)
          ENDIF
-       
+
        ENDIF
 
@@ -5793,5 +5791,5 @@
 
        WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
-       
+
     ENDIF
 
Index: DZ6/trunk/libf/phylmdiso/calcul_divers.h
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/calcul_divers.h	(revision 5292)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/calcul_divers.h
Index: /LMDZ6/trunk/libf/phylmdiso/calcul_divers_mod_h.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/calcul_divers_mod_h.f90	(revision 5293)
+++ /LMDZ6/trunk/libf/phylmdiso/calcul_divers_mod_h.f90	(revision 5293)
@@ -0,0 +1,1 @@
+link ../phylmd/calcul_divers_mod_h.f90
Index: /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 5292)
+++ /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 5293)
@@ -295,6 +295,4 @@
        zustar, zu10m, zv10m, rh2m, qsat2m, &
        zq2m, zt2m, zn2mout, weak_inversion, &
-       zt2m_min_mon, zt2m_max_mon,   &         ! pour calcul_divers.h
-       t2m_min_mon, t2m_max_mon,  &            ! pour calcul_divers.h
        !
        s_pblh_x, s_pblh_w, &
@@ -436,4 +434,5 @@
        USE alpale_mod
     USE yoethf_mod_h
+    USE calcul_divers_mod_h, ONLY: calcul_divers
 
     IMPLICIT NONE
@@ -7093,6 +7092,5 @@
     !
     !IM initialisation + calculs divers diag AMIP2
-    !
-    include "calcul_divers.h"
+    CALL calcul_divers(itap, itapm1, un_jour)
     !
     !IM Interpolation sur les niveaux de pression du NMC
