Index: LMDZ6/branches/Ocean_skin/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/physiq_mod.F90	(revision 3428)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/physiq_mod.F90	(revision 3605)
@@ -25,5 +25,5 @@
     USE dimphy
     USE infotrac_phy, ONLY: nqtot, nbtr, nqo, type_trac
-    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo
+    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, grid1dTo2d_glo, grid_type, unstructured
     USE mod_phys_lmdz_para
     USE iophy
@@ -38,4 +38,7 @@
 #ifdef CPP_Dust
     USE phytracr_spl_mod, ONLY: phytracr_spl
+#endif
+#ifdef CPP_StratAer
+    USE strataer_mod, ONLY: strataer_init
 #endif
     USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
@@ -117,4 +120,6 @@
        zustar, zu10m, zv10m, rh2m, qsat2m, &
        zq2m, zt2m, weak_inversion, &
+       zq2m_cor,zt2m_cor,zu10m_cor,zv10m_cor, & ! pour corriger d'un bug
+       zrh2m_cor,zqsat2m_cor, &
        zt2m_min_mon, zt2m_max_mon,   &         ! pour calcul_divers.h
        t2m_min_mon, t2m_max_mon,  &            ! pour calcul_divers.h
@@ -170,5 +175,5 @@
        !    Deep convective variables used in phytrac
        pmflxr, pmflxs,  &
-       wdtrainA, wdtrainM,  &
+       wdtrainA, wdtrainS, wdtrainM,  &
        upwd, dnwd, &
        ep,  &
@@ -180,4 +185,5 @@
        ev, &
        elij, &
+       qtaa, &
        clw, &
        epmlmMm, eplaMm, &
@@ -243,6 +249,6 @@
 #endif
     USE indice_sol_mod
-    USE phytrac_mod, ONLY : phytrac
-    USE carbon_cycle_mod, ONLY : infocfields_init
+    USE phytrac_mod, ONLY : phytrac_init, phytrac
+    USE carbon_cycle_mod, ONLY : infocfields_init, RCO2_glo, carbon_cycle_rad
 
 #ifdef CPP_RRTM
@@ -265,4 +271,12 @@
     USE ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
     USE VERTICAL_LAYERS_MOD, ONLY: aps,bps
+    USE etat0_limit_unstruct_mod
+#ifdef CPP_XIOS
+    USE xios, ONLY: xios_update_calendar, xios_context_finalize
+#endif
+    USE limit_read_mod, ONLY : init_limit_read
+    USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
+    USE readaerosol_mod, ONLY : init_aero_fromfile
+    USE readaerosolstrato_m, ONLY : init_readaerosolstrato
 
     IMPLICIT NONE
@@ -323,10 +337,11 @@
     include "dimpft.h"
     !======================================================================
+    LOGICAL, SAVE :: ok_volcan ! pour activer les diagnostics volcaniques
+    !$OMP THREADPRIVATE(ok_volcan)
     LOGICAL ok_cvl  ! pour activer le nouveau driver pour convection KE
     PARAMETER (ok_cvl=.TRUE.)
     LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
     PARAMETER (ok_gust=.FALSE.)
-    integer iflag_radia     ! active ou non le rayonnement (MPL)
-    save iflag_radia
+    INTEGER, SAVE :: iflag_radia     ! active ou non le rayonnement (MPL)
     !$OMP THREADPRIVATE(iflag_radia)
     !======================================================================
@@ -363,21 +378,21 @@
     !======================================================================
     LOGICAL ok_journe ! sortir le fichier journalier
-    save ok_journe
+    SAVE ok_journe
     !$OMP THREADPRIVATE(ok_journe)
     !
     LOGICAL ok_mensuel ! sortir le fichier mensuel
-    save ok_mensuel
+    SAVE ok_mensuel
     !$OMP THREADPRIVATE(ok_mensuel)
     !
     LOGICAL ok_instan ! sortir le fichier instantane
-    save ok_instan
+    SAVE ok_instan
     !$OMP THREADPRIVATE(ok_instan)
     !
     LOGICAL ok_LES ! sortir le fichier LES 
-    save ok_LES                            
+    SAVE ok_LES                            
     !$OMP THREADPRIVATE(ok_LES)                  
     !
     LOGICAL callstats ! sortir le fichier stats 
-    save callstats                            
+    SAVE callstats                            
     !$OMP THREADPRIVATE(callstats)                  
     !
@@ -385,9 +400,9 @@
     PARAMETER (ok_region=.FALSE.)
     !======================================================================
-    real seuil_inversion
-    save seuil_inversion
+    REAL seuil_inversion
+    SAVE seuil_inversion
     !$OMP THREADPRIVATE(seuil_inversion)
-    integer iflag_ratqs
-    save iflag_ratqs
+    INTEGER iflag_ratqs
+    SAVE iflag_ratqs
     !$OMP THREADPRIVATE(iflag_ratqs)
     real facteur
@@ -396,9 +411,9 @@
     REAL tau_overturning_th(klon)
 
-    integer lmax_th(klon)
-    integer limbas(klon)
-    real ratqscth(klon,klev)
-    real ratqsdiff(klon,klev)
-    real zqsatth(klon,klev)
+    INTEGER lmax_th(klon)
+    INTEGER limbas(klon)
+    REAL ratqscth(klon,klev)
+    REAL ratqsdiff(klon,klev)
+    REAL zqsatth(klon,klev)
 
     !======================================================================
@@ -497,16 +512,15 @@
     CHARACTER*3 region
     PARAMETER(region='3d')
-    logical ok_hf
-    !
-    save ok_hf
+    LOGICAL ok_hf
+    !
+    SAVE ok_hf
     !$OMP THREADPRIVATE(ok_hf)
 
-    INTEGER,PARAMETER :: longcles=20
-    REAL,SAVE :: clesphy0(longcles)
+    INTEGER, PARAMETER :: longcles=20
+    REAL, SAVE :: clesphy0(longcles)
     !$OMP THREADPRIVATE(clesphy0)
     !
     ! Variables propres a la physique
-    INTEGER itap
-    SAVE itap                   ! compteur pour la physique
+    INTEGER, SAVE :: itap         ! compteur pour la physique
     !$OMP THREADPRIVATE(itap)
 
@@ -514,5 +528,5 @@
     !$OMP THREADPRIVATE(abortphy)
     !
-    REAL,save ::  solarlong0
+    REAL,SAVE ::  solarlong0
     !$OMP THREADPRIVATE(solarlong0)
 
@@ -531,7 +545,5 @@
     ! Variables liees a la convection de K. Emanuel (sb):
     !
-    REAL bas, top             ! cloud base and top levels
-    SAVE bas
-    SAVE top
+    REAL, SAVE :: bas, top             ! cloud base and top levels
     !$OMP THREADPRIVATE(bas, top)
     !------------------------------------------------------------------
@@ -551,5 +563,6 @@
     ! Variables li\'ees \`a la poche froide (jyg)
 
-    REAL mip(klon,klev)  ! mass flux shed by the adiab ascent at each level
+!!    REAL mipsh(klon,klev)  ! mass flux shed by the adiab ascent at each level
+!!      Moved to phys_state_var_mod
     !
     REAL wape_prescr, fip_prescr
@@ -568,14 +581,11 @@
 !!    REAL, DIMENSION(klon,klev)     :: dql_sat
 
-    real, save :: alp_bl_prescr=0.
-    real, save :: ale_bl_prescr=0.
-
-    real, save :: wake_s_min_lsp=0.1
-
+    REAL, SAVE :: alp_bl_prescr=0.
+    REAL, SAVE :: ale_bl_prescr=0.
+    REAL, SAVE :: wake_s_min_lsp=0.1
     !$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
     !$OMP THREADPRIVATE(wake_s_min_lsp)
 
-
-    real ok_wk_lsp(klon)
+    REAL ok_wk_lsp(klon)
 
     !RC
@@ -590,4 +600,10 @@
                                                         ! gust-front in the grid cell.
     !$OMP THREADPRIVATE(iflag_alp_wk_cond)
+
+    INTEGER,  SAVE               :: iflag_bug_t2m_ipslcm61=1 !
+    !$OMP THREADPRIVATE(iflag_bug_t2m_ipslcm61)
+    INTEGER,  SAVE               :: iflag_bug_t2m_stab_ipslcm61=-1 !
+    !$OMP THREADPRIVATE(iflag_bug_t2m_stab_ipslcm61)
+
     REAL t_w(klon,klev),q_w(klon,klev) ! temperature and moisture profiles in the wake region
     REAL t_x(klon,klev),q_x(klon,klev) ! temperature and moisture profiles in the off-wake region
@@ -727,5 +743,5 @@
     REAL :: jD_eq
 
-    LOGICAL, parameter :: new_orbit = .true.
+    LOGICAL, parameter :: new_orbit = .TRUE.
 
     !
@@ -913,6 +929,6 @@
     INTEGER kcbot(klon), kctop(klon), kdtop(klon)
     !
-    real ratqsbas,ratqshaut,tau_ratqs
-    save ratqsbas,ratqshaut,tau_ratqs
+    REAL ratqsbas,ratqshaut,tau_ratqs
+    SAVE ratqsbas,ratqshaut,tau_ratqs
     !$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
     REAL, SAVE :: ratqsp0=50000., ratqsdp=20000.
@@ -920,19 +936,15 @@
 
     ! Parametres lies au nouveau schema de nuages (SB, PDF)
-    real fact_cldcon
-    real facttemps
-    logical ok_newmicro
-    save ok_newmicro
+    REAL, SAVE :: fact_cldcon
+    REAL, SAVE :: facttemps
+    !$OMP THREADPRIVATE(fact_cldcon,facttemps)
+    LOGICAL, SAVE :: ok_newmicro
     !$OMP THREADPRIVATE(ok_newmicro)
-    !real ref_liq_pi(klon,klev), ref_ice_pi(klon,klev)
-    save fact_cldcon,facttemps
-    !$OMP THREADPRIVATE(fact_cldcon,facttemps)
-
-    integer iflag_cld_th
-    save iflag_cld_th
+
+    INTEGER, SAVE :: iflag_cld_th
     !$OMP THREADPRIVATE(iflag_cld_th)
 !IM logical ptconv(klon,klev)  !passe dans phys_local_var_mod
     !IM cf. AM 081204 BEG
-    logical ptconvth(klon,klev)
+    LOGICAL ptconvth(klon,klev)
     !IM cf. AM 081204 END
     !
@@ -941,5 +953,4 @@
     !======================================================================
     !
-
     !
 !JLD    integer itau_w   ! pas de temps ecriture = itap + itau_phy
@@ -1007,16 +1018,15 @@
 !JLD    REAL zstophy, zout
 
-    character*20 modname
-    character*80 abort_message
-    logical, save ::  ok_sync, ok_sync_omp
+    CHARACTER*20 modname
+    CHARACTER*80 abort_message
+    LOGICAL, SAVE ::  ok_sync, ok_sync_omp
     !$OMP THREADPRIVATE(ok_sync)
-    real date0
+    REAL date0
 
     ! essai writephys
-    integer fid_day, fid_mth, fid_ins
-    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)
+    INTEGER fid_day, fid_mth, fid_ins
+    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)
     REAL ztsol(klon)
     REAL q2m(klon,nbsrf)  ! humidite a 2m
@@ -1070,5 +1080,5 @@
     ! Declaration des constantes et des fonctions thermodynamiques
     !
-    LOGICAL,SAVE :: first=.true.
+    LOGICAL,SAVE :: first=.TRUE.
     !$OMP THREADPRIVATE(first)
 
@@ -1106,5 +1116,5 @@
     ! Declarations pour Simulateur COSP
     !============================================================
-    real :: mr_ozone(klon,klev)
+    real :: mr_ozone(klon,klev), phicosp(klon,klev)
 
     !IM stations CFMIP
@@ -1164,10 +1174,11 @@
     REAL zzz
     !albedo SB >>>
-    real,dimension(6),save :: SFRWL
+    REAL,DIMENSION(6), SAVE :: SFRWL
+!$OMP THREADPRIVATE(SFRWL)
     !albedo SB <<<
 
     !--OB variables for mass fixer (hard coded for now)
-    logical, parameter :: mass_fixer=.false.
-    real qql1(klon),qql2(klon),corrqql
+    LOGICAL, PARAMETER :: mass_fixer=.FALSE.
+    REAL qql1(klon),qql2(klon),corrqql
 
     REAL pi
@@ -1183,4 +1194,8 @@
     pdtphys=pdtphys_
     CALL update_time(pdtphys)
+    phys_tstep=NINT(pdtphys)
+#ifdef CPP_XIOS
+    IF (.NOT. debut .AND. is_omp_master) CALL xios_update_calendar(itap+1)
+#endif
 
     !======================================================================
@@ -1211,15 +1226,17 @@
 
     ! Quick check on pressure levels:
-    call assert(paprs(:, nbp_lev + 1) < paprs(:, nbp_lev), &
+    CALL assert(paprs(:, nbp_lev + 1) < paprs(:, nbp_lev), &
             "physiq_mod paprs bad order")
 
     IF (first) THEN 
+       CALL init_etat0_limit_unstruct
+       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
        !CR:nvelles variables convection/poches froides
 
-       print*, '================================================='
-       print*, 'Allocation des variables locales et sauvegardees'
+       WRITE(lunout,*) '================================================='
+       WRITE(lunout,*) 'Allocation des variables locales et sauvegardees'
+       WRITE(lunout,*) '================================================='
        CALL phys_local_var_init
        !
-       pasphys=pdtphys
        !     appel a la lecture du run.def physique
        CALL conf_phys(ok_journe, ok_mensuel, &
@@ -1230,5 +1247,6 @@
             fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
             iflag_cld_th,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
-            ok_ade, ok_aie, ok_alw, ok_cdnc, aerosol_couple,  chemistry_couple, &
+            ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, aerosol_couple, &
+            chemistry_couple, &
             flag_aerosol, flag_aerosol_strat, flag_aer_feedback, new_aod, &
             flag_bc_internal_mixture, bl95_b0, bl95_b1, &
@@ -1239,4 +1257,11 @@
        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 CPP_StratAer
+       CALL strataer_init
+#endif
+
        print*, '================================================='
        !
@@ -1245,5 +1270,6 @@
           WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
                '(H2Ov, H2Ol, H2Oi) but nqo=', nqo, '. Might as well stop here.'
-          STOP
+          abort_message='see above'
+          CALL abort_physic(modname,abort_message,1)
        ENDIF
 
@@ -1258,5 +1284,5 @@
 
        itau_con=0
-       first=.false.
+       first=.FALSE.
 
     ENDIF  ! first
@@ -1287,5 +1313,11 @@
 ! secondes
        tau_gl=86400.*tau_gl
-       print*,'debut physiq_mod tau_gl=',tau_gl
+       WRITE(lunout,*) 'debut physiq_mod tau_gl=',tau_gl
+
+       iflag_bug_t2m_ipslcm61 = 1
+       CALL getin_p('iflag_bug_t2m_ipslcm61', iflag_bug_t2m_ipslcm61)
+       iflag_bug_t2m_stab_ipslcm61 = -1
+       CALL getin_p('iflag_bug_t2m_stab_ipslcm61', iflag_bug_t2m_stab_ipslcm61)
+
        CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond) 
        CALL getin_p('random_notrig_max',random_notrig_max)
@@ -1318,4 +1350,21 @@
        CALL getin_p('NVM',nvm_lmdz)
 
+       WRITE(lunout,*) 'iflag_alp_wk_cond=',  iflag_alp_wk_cond
+       WRITE(lunout,*) 'random_ntrig_max=',   random_notrig_max
+       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,*) 'ratqsp0=',            ratqsp0
+       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_bug_split_th=',    ok_bug_split_th
+       WRITE(lunout,*) 'fl_ebil=',            fl_ebil
+       WRITE(lunout,*) 'fl_cor_ebil=',        fl_cor_ebil
+       WRITE(lunout,*) 'iflag_phytrac=',      iflag_phytrac
+       WRITE(lunout,*) 'NVM=',                nvm_lmdz
+
        !--PC: defining fields to be exchanged between LMDz, ORCHIDEE and NEMO
        WRITE(lunout,*) 'Call to infocfields from physiq'
@@ -1369,4 +1418,8 @@
        ENDIF
 
+       tau_aero(:,:,:,:) = 1.e-15
+       piz_aero(:,:,:,:) = 1.
+       cg_aero(:,:,:,:)  = 0.
+
        IF (aerosol_couple .AND. (config_inca /= "aero" &
             .AND. config_inca /= "aeNP ")) THEN
@@ -1376,6 +1429,4 @@
           CALL abort_physic (modname,abort_message,1)
        ENDIF
-
-
 
        rnebcon0(:,:) = 0.0
@@ -1417,14 +1468,12 @@
        ! pour obtenir le meme resultat.
 !jyg for fh<
-!!       dtime=pdtphys
-       dtime=NINT(pdtphys)
-       WRITE(lunout,*) 'Pas de temps dtime pdtphys ',dtime,pdtphys
-       IF (abs(dtime-pdtphys)>1.e-10) THEN
+       WRITE(lunout,*) 'Pas de temps phys_tstep pdtphys ',phys_tstep,pdtphys
+       IF (abs(phys_tstep-pdtphys)>1.e-10) THEN
           abort_message='pas de temps doit etre entier en seconde pour orchidee et XIOS'
           CALL abort_physic(modname,abort_message,1)
        ENDIF
 !>jyg
-       IF (MOD(NINT(86400./dtime),nbapp_rad).EQ.0) THEN
-          radpas = NINT( 86400./dtime)/nbapp_rad
+       IF (MOD(NINT(86400./phys_tstep),nbapp_rad).EQ.0) THEN
+          radpas = NINT( 86400./phys_tstep)/nbapp_rad
        ELSE 
           WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
@@ -1436,9 +1485,9 @@
           CALL abort_physic(modname,abort_message,1)
        ENDIF
-       IF (nbapp_cv .EQ. 0) nbapp_cv=86400./dtime
-       IF (nbapp_wk .EQ. 0) nbapp_wk=86400./dtime
+       IF (nbapp_cv .EQ. 0) nbapp_cv=86400./phys_tstep
+       IF (nbapp_wk .EQ. 0) nbapp_wk=86400./phys_tstep
        print *,'physiq, nbapp_cv, nbapp_wk ',nbapp_cv,nbapp_wk
-       IF (MOD(NINT(86400./dtime),nbapp_cv).EQ.0) THEN
-          cvpas_0 = NINT( 86400./dtime)/nbapp_cv
+       IF (MOD(NINT(86400./phys_tstep),nbapp_cv).EQ.0) THEN
+          cvpas_0 = NINT( 86400./phys_tstep)/nbapp_cv
           cvpas = cvpas_0
        print *,'physiq, cvpas ',cvpas
@@ -1450,9 +1499,9 @@
           abort_message='nbre de pas de temps physique n est pas multiple ' &
                // 'de nbapp_cv'
-          call abort_physic(modname,abort_message,1)
-       ENDIF
-       IF (MOD(NINT(86400./dtime),nbapp_wk).EQ.0) THEN
-          wkpas = NINT( 86400./dtime)/nbapp_wk
-       print *,'physiq, wkpas ',wkpas
+          CALL abort_physic(modname,abort_message,1)
+       ENDIF
+       IF (MOD(NINT(86400./phys_tstep),nbapp_wk).EQ.0) THEN
+          wkpas = NINT( 86400./phys_tstep)/nbapp_wk
+!       print *,'physiq, wkpas ',wkpas
        ELSE 
           WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
@@ -1462,22 +1511,208 @@
           abort_message='nbre de pas de temps physique n est pas multiple ' &
                // 'de nbapp_wk'
-          call abort_physic(modname,abort_message,1)
+          CALL abort_physic(modname,abort_message,1)
        ENDIF
        ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+       CALL init_iophy_new(latitude_deg,longitude_deg)
+
+          !===================================================================
+          !IM stations CFMIP
+          nCFMIP=npCFMIP
+          OPEN(98,file='npCFMIP_param.data',status='old', &
+               form='formatted',iostat=iostat)
+          IF (iostat == 0) THEN
+             READ(98,*,end=998) nCFMIP
+998          CONTINUE
+             CLOSE(98)
+             CONTINUE
+             IF(nCFMIP.GT.npCFMIP) THEN
+                print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
+                CALL abort_physic("physiq", "", 1)
+             ELSE
+                print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
+             ENDIF
+
+             !
+             ALLOCATE(tabCFMIP(nCFMIP))
+             ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
+             ALLOCATE(tabijGCM(nCFMIP))
+             ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
+             ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
+             !
+             ! lecture des nCFMIP stations CFMIP, de leur numero 
+             ! et des coordonnees geographiques lonCFMIP, latCFMIP
+             !
+             CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP,  &
+                  lonCFMIP, latCFMIP)
+             !
+             ! identification des
+             ! 1) coordonnees lonGCM, latGCM des points CFMIP dans la
+             ! grille de LMDZ
+             ! 2) indices points tabijGCM de la grille physique 1d sur
+             ! klon points
+             ! 3) indices iGCM, jGCM de la grille physique 2d
+             !
+             CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
+                  tabijGCM, lonGCM, latGCM, iGCM, jGCM)
+             !
+          ELSE
+             ALLOCATE(tabijGCM(0))
+             ALLOCATE(lonGCM(0), latGCM(0))
+             ALLOCATE(iGCM(0), jGCM(0))
+          ENDIF
+
+#ifdef CPP_IOIPSL
+
+       !$OMP MASTER
+       ! FH : if ok_sync=.true. , the time axis is written at each time step
+       ! in the output files. Only at the end in the opposite case
+       ok_sync_omp=.FALSE.
+       CALL getin('ok_sync',ok_sync_omp)
+       CALL phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
+            iGCM,jGCM,lonGCM,latGCM, &
+            jjmp1,nlevSTD,clevSTD,rlevSTD, phys_tstep,ok_veget, &
+            type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
+            ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
+            read_climoz, phys_out_filestations, &
+            new_aod, aerosol_couple, &
+            flag_aerosol_strat, pdtphys, paprs, pphis,  &
+            pplay, lmax_th, ptconv, ptconvth, ivap,  &
+            d_u, d_t, qx, d_qx, zmasse, ok_sync_omp)
+       !$OMP END MASTER
+       !$OMP BARRIER
+       ok_sync=ok_sync_omp
+
+       freq_outNMC(1) = ecrit_files(7)
+       freq_outNMC(2) = ecrit_files(8)
+       freq_outNMC(3) = ecrit_files(9)
+       WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
+       WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
+       WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
+
+#ifndef CPP_XIOS
+       CALL ini_paramLMDZ_phy(phys_tstep,nid_ctesGCM)
+#endif
+
+#endif
+       ecrit_reg = ecrit_reg * un_jour
+       ecrit_tra = ecrit_tra * un_jour
+
+       !XXXPB Positionner date0 pour initialisation de ORCHIDEE
+       date0 = jD_ref 
+       WRITE(*,*) 'physiq date0 : ',date0
+       !
+
+!       CALL create_climoz(read_climoz)
+      IF (.NOT. create_etat0_limit) CALL init_aero_fromfile(flag_aerosol)  !! initialise aero from file for XIOS interpolation (unstructured_grid)
+      IF (.NOT. create_etat0_limit) CALL init_readaerosolstrato(flag_aerosol_strat)  !! initialise aero strato from file for XIOS interpolation (unstructured_grid)
+
+#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
+        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,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_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, &
+               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, &
+               phicosp,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
+
+       !
+       !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+       ! Nouvelle initialisation pour le rayonnement RRTM
+       !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+       CALL iniradia(klon,klev,paprs(1,1:klev+1))
+       ! Initialisation des champs dans phytrac qui sont utilisés par phys_output_write
+       IF (iflag_phytrac == 1 ) THEN
+          CALL phytrac_init()
+        ENDIF
+
+       CALL phys_output_write(itap, pdtphys, paprs, pphis,                    &
+                              pplay, lmax_th, aerosol_couple,                 &
+                              ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, new_aod, ok_sync,&
+                              ptconv, read_climoz, clevSTD,                   &
+                              ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
+                              flag_aerosol, flag_aerosol_strat, ok_cdnc)
+
+#ifdef CPP_XIOS
+       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
        CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
+
 !jyg<
-       IF (klon_glo==1) THEN
-          IF (iflag_pbl > 1) THEN          
-              pbl_tke(:,:,is_ave) = 0.
-              DO nsrf=1,nbsrf
-                DO k = 1,klev+1
-                     pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
-                         +pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
-                ENDDO
-              ENDDO
-          ELSE   ! (iflag_pbl > 1)
-              pbl_tke(:,:,:) = 0.
-          ENDIF  ! (iflag_pbl > 1)
+       IF (iflag_pbl<=1) THEN
+          ! No TKE for Standard Physics
+          pbl_tke(:,:,:)=0.
+
+       ELSE IF (klon_glo==1) THEN
+          pbl_tke(:,:,is_ave) = 0.
+          DO nsrf=1,nbsrf
+            DO k = 1,klev+1
+                 pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
+                     +pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
+            ENDDO
+          ENDDO
+        ELSE
+          pbl_tke(:,:,is_ave) = 0. !ym missing init : maybe must be initialized in the same way that for klon_glo==1 ??
 !>jyg
        ENDIF
@@ -1496,14 +1731,11 @@
        ENDIF
 
-       CALL printflag( tabcntr0,radpas,ok_journe, &
-            ok_instan, ok_region )
-       !
-       IF (ABS(dtime-pdtphys).GT.0.001) THEN
-          WRITE(lunout,*) 'Pas physique n est pas correct',dtime, &
-               pdtphys
-          abort_message='Pas physique n est pas correct '
-          !           call abort_physic(modname,abort_message,1)
-          dtime=pdtphys
-       ENDIF
+!       IF (ABS(phys_tstep-pdtphys).GT.0.001) THEN
+!          WRITE(lunout,*) 'Pas physique n est pas correct',phys_tstep, &
+!               pdtphys
+!          abort_message='Pas physique n est pas correct '
+!          !           call abort_physic(modname,abort_message,1)
+!          phys_tstep=pdtphys
+!       ENDIF
        IF (nlon .NE. klon) THEN
           WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,  &
@@ -1519,5 +1751,5 @@
        ENDIF
        !
-       IF (dtime*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"
@@ -1581,53 +1813,8 @@
           !        enddo
 
-          !===================================================================
-          !IM stations CFMIP
-          nCFMIP=npCFMIP
-          OPEN(98,file='npCFMIP_param.data',status='old', &
-               form='formatted',iostat=iostat)
-          IF (iostat == 0) THEN
-             READ(98,*,end=998) nCFMIP
-998          CONTINUE
-             CLOSE(98)
-             CONTINUE
-             IF(nCFMIP.GT.npCFMIP) THEN
-                print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
-                CALL abort_physic("physiq", "", 1)
-             ELSE
-                print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
-             ENDIF
-
-             !
-             ALLOCATE(tabCFMIP(nCFMIP))
-             ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
-             ALLOCATE(tabijGCM(nCFMIP))
-             ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
-             ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
-             !
-             ! lecture des nCFMIP stations CFMIP, de leur numero 
-             ! et des coordonnees geographiques lonCFMIP, latCFMIP
-             !
-             CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP,  &
-                  lonCFMIP, latCFMIP)
-             !
-             ! identification des
-             ! 1) coordonnees lonGCM, latGCM des points CFMIP dans la
-             ! grille de LMDZ
-             ! 2) indices points tabijGCM de la grille physique 1d sur
-             ! klon points
-             ! 3) indices iGCM, jGCM de la grille physique 2d
-             !
-             CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
-                  tabijGCM, lonGCM, latGCM, iGCM, jGCM)
-             !
-          ELSE
-             ALLOCATE(tabijGCM(0))
-             ALLOCATE(lonGCM(0), latGCM(0))
-             ALLOCATE(iGCM(0), jGCM(0))
-          ENDIF
-       ELSE
-          ALLOCATE(tabijGCM(0))
-          ALLOCATE(lonGCM(0), latGCM(0))
-          ALLOCATE(iGCM(0), jGCM(0))
+       !ELSE
+       !   ALLOCATE(tabijGCM(0))
+       !   ALLOCATE(lonGCM(0), latGCM(0))
+       !   ALLOCATE(iGCM(0), jGCM(0))
        ENDIF
 
@@ -1665,5 +1852,5 @@
        !
        !
-       lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
+       lmt_pas = NINT(86400./phys_tstep * 1.0)   ! tous les jours
        WRITE(lunout,*)'La frequence de lecture surface est de ',  &
             lmt_pas
@@ -1681,4 +1868,13 @@
        !   Initialisation des sorties
        !=============================================================
+
+#ifdef CPP_XIOS
+       ! Get "missing_val" value from XML files (from temperature variable)
+       !$OMP MASTER
+       CALL xios_get_field_attr("temp",default_value=missing_val_omp)
+       !$OMP END MASTER
+       !$OMP BARRIER
+       missing_val=missing_val_omp
+#endif
 
 #ifdef CPP_XIOS 
@@ -1693,43 +1889,7 @@
 #endif 
 
-#ifdef CPP_IOIPSL
-
-       !$OMP MASTER
-       ! FH : if ok_sync=.true. , the time axis is written at each time step
-       ! in the output files. Only at the end in the opposite case
-       ok_sync_omp=.false.
-       CALL getin('ok_sync',ok_sync_omp)
-       CALL phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
-            iGCM,jGCM,lonGCM,latGCM, &
-            jjmp1,nlevSTD,clevSTD,rlevSTD, dtime,ok_veget, &
-            type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
-            ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
-            read_climoz, phys_out_filestations, &
-            new_aod, aerosol_couple, &
-            flag_aerosol_strat, pdtphys, paprs, pphis,  &
-            pplay, lmax_th, ptconv, ptconvth, ivap,  &
-            d_u, d_t, qx, d_qx, zmasse, ok_sync_omp)
-       !$OMP END MASTER
-       !$OMP BARRIER
-       ok_sync=ok_sync_omp
-
-       freq_outNMC(1) = ecrit_files(7)
-       freq_outNMC(2) = ecrit_files(8)
-       freq_outNMC(3) = ecrit_files(9)
-       WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
-       WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
-       WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
-
-#ifndef CPP_XIOS
-       CALL ini_paramLMDZ_phy(dtime,nid_ctesGCM)
-#endif
-
-#endif
-       ecrit_reg = ecrit_reg * un_jour
-       ecrit_tra = ecrit_tra * un_jour
-
-       !XXXPB Positionner date0 pour initialisation de ORCHIDEE
-       date0 = jD_ref 
-       WRITE(*,*) 'physiq date0 : ',date0
+
+       CALL printflag( tabcntr0,radpas,ok_journe, &
+            ok_instan, ok_region )
        !
        !
@@ -1792,10 +1952,4 @@
 #endif
        ENDIF
-       !
-       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-       ! Nouvelle initialisation pour le rayonnement RRTM
-       ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-       CALL iniradia(klon,klev,paprs(1,1:klev+1))
 
        !$omp single
@@ -1815,14 +1969,14 @@
 
        !albedo SB >>>
-       select case(nsw)
-       case(2)
+       SELECT CASE(nsw)
+       CASE(2)
           SFRWL(1)=0.45538747
           SFRWL(2)=0.54461211
-       case(4)
+       CASE(4)
           SFRWL(1)=0.45538747
           SFRWL(2)=0.32870591
           SFRWL(3)=0.18568763
           SFRWL(4)=3.02191470E-02
-       case(6)
+       CASE(6)
           SFRWL(1)=1.28432794E-03
           SFRWL(2)=0.12304168
@@ -1831,5 +1985,5 @@
           SFRWL(5)=0.18568763
           SFRWL(6)=3.02191470E-02
-       end select
+       END SELECT
 
 
@@ -1870,4 +2024,6 @@
       sollwdown(:)= sollwdown(:) + betalwoff *(-1.*ZFLDN0(:,1) - &
                     sollwdown(:))
+
+
     ENDIF
     !
@@ -1891,5 +2047,5 @@
     ! on the surface fraction.
     !
-    CALL change_srf_frac(itap, dtime, days_elapsed+1,  &
+    CALL change_srf_frac(itap, phys_tstep, days_elapsed+1,  &
          pctsrf, fevap, z0m, z0h, agesno,              &
          falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
@@ -1903,5 +2059,4 @@
 #endif
     ENDIF
-
 
     ! Tendances bidons pour les processus qui n'affectent pas certaines
@@ -2005,4 +2160,13 @@
        ENDDO
     ENDIF
+!
+! Temporary solutions adressing ticket #104 and the non initialisation of tr_ancien
+! LF
+    IF (debut) THEN
+      WRITE(lunout,*)' WARNING: tr_ancien initialised to tr_seri'
+      DO iq = nqo+1, nqtot
+           tr_ancien(:,:,iq-nqo)=tr_seri(:,:,iq-nqo)
+      ENDDO
+    ENDIF
     !
     DO i = 1, klon
@@ -2021,20 +2185,20 @@
     IF (ancien_ok) THEN
     !
-       d_u_dyn(:,:)  = (u_seri(:,:)-u_ancien(:,:))/dtime
-       d_v_dyn(:,:)  = (v_seri(:,:)-v_ancien(:,:))/dtime
-       d_t_dyn(:,:)  = (t_seri(:,:)-t_ancien(:,:))/dtime
-       d_q_dyn(:,:)  = (q_seri(:,:)-q_ancien(:,:))/dtime
-       d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/dtime
-       d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/dtime
+       d_u_dyn(:,:)  = (u_seri(:,:)-u_ancien(:,:))/phys_tstep
+       d_v_dyn(:,:)  = (v_seri(:,:)-v_ancien(:,:))/phys_tstep
+       d_t_dyn(:,:)  = (t_seri(:,:)-t_ancien(:,:))/phys_tstep
+       d_q_dyn(:,:)  = (q_seri(:,:)-q_ancien(:,:))/phys_tstep
+       d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
+       d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
        CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
-       d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/dtime
+       d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
        CALL water_int(klon,klev,ql_seri,zmasse,zx_tmp_fi2d)
-       d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/dtime
+       d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/phys_tstep
        CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
-       d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/dtime
+       d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
        ! !! RomP >>>   td dyn traceur
        IF (nqtot.GT.nqo) THEN     ! jyg
           DO iq = nqo+1, nqtot      ! jyg
-              d_tr_dyn(:,:,iq-nqo)=(tr_seri(:,:,iq-nqo)-tr_ancien(:,:,iq-nqo))/dtime ! jyg
+              d_tr_dyn(:,:,iq-nqo)=(tr_seri(:,:,iq-nqo)-tr_ancien(:,:,iq-nqo))/phys_tstep ! jyg
           ENDDO
        ENDIF
@@ -2138,5 +2302,5 @@
                       ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1),    &
                       time_climoz )
-          END IF
+          ENDIF
           ! Convert from mole fraction of ozone to column density of ozone in a
           ! cell, in kDU:
@@ -2158,5 +2322,5 @@
             (du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
                'eva',abortphy,flag_inhib_tend,itap,0)
-    call prt_enerbil('eva',itap)
+    CALL prt_enerbil('eva',itap)
 
     !=========================================================================
@@ -2213,5 +2377,5 @@
           !  bit comparable a l ancienne formulation cycle_diurne=true
           !  on integre entre gmtime et gmtime+radpas
-          zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
+          zdtime=phys_tstep*REAL(radpas) ! pas de temps du rayonnement (s)
           CALL zenang(zlongi,jH_cur,0.0,zdtime, &
                latitude_deg,longitude_deg,rmu0,fract)
@@ -2230,6 +2394,6 @@
           !  premier pas de temps de la physique pendant lequel
           !  itaprad=0
-          zdtime1=dtime*REAL(-MOD(itaprad,radpas)-1)      
-          zdtime2=dtime*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)
@@ -2237,5 +2401,5 @@
           ! Calcul des poids
           !
-          zdtime1=-dtime !--on corrige le rayonnement pour representer le
+          zdtime1=-phys_tstep !--on corrige le rayonnement pour representer le
           zdtime2=0.0    !--pas de temps de la physique qui se termine
           CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
@@ -2295,4 +2459,7 @@
        !
        !-------gustiness calculation-------!
+       !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
@@ -2312,8 +2479,6 @@
        ENDIF
 
-
-
        CALL pbl_surface(  &
-            dtime,     date0,     itap,    days_elapsed+1, &
+            phys_tstep,     date0,     itap,    days_elapsed+1, &
             debut,     lafin, &
             longitude_deg, latitude_deg, rugoro,  zrmu0,      &
@@ -2382,7 +2547,20 @@
        ENDIF
 
-
-
-
+!add limitation for t,q at and wind at 10m
+        if ( iflag_bug_t2m_ipslcm61 == 0 ) THEN
+          CALL borne_var_surf( klon,klev,nbsrf,                 &
+            iflag_bug_t2m_stab_ipslcm61,                        &
+            t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),    &
+            ftsol,zxqsurf,pctsrf,paprs,                         &
+            t2m, q2m, u10m, v10m,                               &
+            zt2m_cor, zq2m_cor, zu10m_cor, zv10m_cor,           &
+            zrh2m_cor, zqsat2m_cor)
+        ELSE
+          zt2m_cor(:)=zt2m(:)
+          zq2m_cor(:)=zq2m(:)
+          zu10m_cor(:)=zu10m(:)
+          zv10m_cor(:)=zv10m(:)
+          zqsat2m_cor=999.999
+        ENDIF
 
        !---------------------------------------------------------------------
@@ -2397,5 +2575,5 @@
                'vdf',abortphy,flag_inhib_tend,itap,0)
        ENDIF
-       call prt_enerbil('vdf',itap)
+       CALL prt_enerbil('vdf',itap)
        !--------------------------------------------------------------------
 
@@ -2482,7 +2660,7 @@
        DO i = 1, klon
           conv_q(i,k) = d_q_dyn(i,k)  &
-               + d_q_vdf(i,k)/dtime
+               + d_q_vdf(i,k)/phys_tstep
           conv_t(i,k) = d_t_dyn(i,k)  &
-               + d_t_vdf(i,k)/dtime
+               + d_t_vdf(i,k)/phys_tstep
        ENDDO
     ENDDO
@@ -2528,4 +2706,5 @@
     pmflxs(:,:) = 0.
     wdtrainA(:,:) = 0.
+    wdtrainS(:,:) = 0.
     wdtrainM(:,:) = 0.
     upwd(:,:) = 0.
@@ -2543,4 +2722,5 @@
     elij(:,:,:)=0.
     ev(:,:)=0.
+    qtaa(:,:)=0.
     clw(:,:)=0.
     sij(:,:,:)=0.
@@ -2549,9 +2729,9 @@
        abort_message ='reactiver le call conlmd dans physiq.F'
        CALL abort_physic (modname,abort_message,1)
-       !     CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
+       !     CALL conlmd (phys_tstep, paprs, pplay, t_seri, q_seri, conv_q,
        !    .             d_t_con, d_q_con,
        !    .             rain_con, snow_con, ibas_con, itop_con)
     ELSE IF (iflag_con.EQ.2) THEN
-       CALL conflx(dtime, paprs, pplay, t_seri, q_seri, &
+       CALL conflx(phys_tstep, paprs, pplay, t_seri, q_seri, &
             conv_t, conv_q, -evap, omega, &
             d_t_con, d_q_con, rain_con, snow_con, &
@@ -2629,5 +2809,5 @@
 
 !jyg<
-       CALL alpale( debut, itap, dtime, paprs, omega, t_seri,   &
+       CALL alpale( debut, itap, phys_tstep, paprs, omega, t_seri,   &
                     alp_offset, it_wape_prescr,  wape_prescr, fip_prescr, &
                     ale_bl_prescr, alp_bl_prescr, &
@@ -2671,5 +2851,5 @@
           !c          CALL concvl (iflag_con,iflag_clos,
           CALL concvl (iflag_clos, &
-               dtime, paprs, pplay, k_upper_cv, t_x,q_x, &
+               phys_tstep, paprs, pplay, k_upper_cv, t_x,q_x, &
                t_w,q_w,wake_s, &
                u_seri,v_seri,tr_seri,nbtr_tmp, &
@@ -2679,13 +2859,13 @@
                rain_con, snow_con, ibas_con, itop_con, sigd, &
                ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
-               Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
+               Ma,mipsh,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
                pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
                                 ! RomP >>>
                                 !!     .        pmflxr,pmflxs,da,phi,mp,
                                 !!     .        ftd,fqd,lalim_conv,wght_th)
-               pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
+               pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,qtaa,clw,elij, &
                ftd,fqd,lalim_conv,wght_th, &
                ev, ep,epmlmMm,eplaMm, &
-               wdtrainA,wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
+               wdtrainA, wdtrainS, wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
                tau_cld_cv,coefw_cld_cv,epmax_diag)
 
@@ -2741,6 +2921,6 @@
               DO k=1,klev
                  DO i=1,klon
-                    ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/dtime
-                    fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/dtime
+                    ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/phys_tstep
+                    fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/phys_tstep
                     d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
                     d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
@@ -2754,5 +2934,5 @@
 
           ! MAF conema3 ne contient pas les traceurs
-          CALL conema3 (dtime, &
+          CALL conema3 (phys_tstep, &
                paprs,pplay,t_seri,q_seri, &
                u_seri,v_seri,tr_seri,ntra, &
@@ -2865,5 +3045,5 @@
     CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
          'convection',abortphy,flag_inhib_tend,itap,0)
-    call prt_enerbil('convection',itap)
+    CALL prt_enerbil('convection',itap)
 
     !-------------------------------------------------------------------------
@@ -2886,5 +3066,5 @@
                snow_con(i))*cell_area(i)/REAL(klon)
        ENDDO
-       zx_t = zx_t/za*dtime
+       zx_t = zx_t/za*phys_tstep
        WRITE(lunout,*)"Precip=", zx_t
     ENDIF
@@ -2900,5 +3080,5 @@
        ENDDO
        DO i = 1, klon
-          z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
+          z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*phys_tstep) &
                /z_apres(i)
        ENDDO
@@ -2937,6 +3117,6 @@
                 M_dwn(i,k)   = dnwd0(i,k)
                 M_up(i,k)    = upwd(i,k)
-                dt_a(i,k)    = d_t_con(i,k)/dtime - ftd(i,k) 
-                dq_a(i,k)    = d_q_con(i,k)/dtime - fqd(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
@@ -2946,7 +3126,7 @@
              DO k = 1,klev
                 dt_dwn(:,k)= dt_dwn(:,k)+ &
-                     ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
+                     ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/phys_tstep
                 dq_dwn(:,k)= dq_dwn(:,k)+ &
-                     ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
+                     ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/phys_tstep
              ENDDO
           ELSEIF (iflag_wake==3) THEN
@@ -2959,7 +3139,7 @@
                       ! l'eau se reevapore).
                       dt_dwn(i,k)= dt_dwn(i,k)+ &
-                           ok_wk_lsp(i)*d_t_lsc(i,k)/dtime
+                           ok_wk_lsp(i)*d_t_lsc(i,k)/phys_tstep
                       dq_dwn(i,k)= dq_dwn(i,k)+ &
-                           ok_wk_lsp(i)*d_q_lsc(i,k)/dtime
+                           ok_wk_lsp(i)*d_q_lsc(i,k)/phys_tstep
                    ENDIF
                 ENDDO
@@ -2969,5 +3149,5 @@
           !
           !calcul caracteristiques de la poche froide
-          CALL calWAKE (iflag_wake_tend, paprs, pplay, dtime, &
+          CALL calWAKE (iflag_wake_tend, paprs, pplay, phys_tstep, &
                t_seri, q_seri, omega,  &
                dt_dwn, dq_dwn, M_dwn, M_up,  &
@@ -2996,5 +3176,5 @@
        CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
             abortphy,flag_inhib_tend,itap,0)
-       call prt_enerbil('wake',itap)
+       CALL prt_enerbil('wake',itap)
        !------------------------------------------------------------------------
 
@@ -3005,5 +3185,5 @@
             (d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk, wake_k, &
              'wake', abortphy) 
-          call prt_enerbil('wake',itap)
+          CALL prt_enerbil('wake',itap)
        ENDIF   ! (iflag_wake_tend .GT. 0.)
        !
@@ -3016,5 +3196,5 @@
        IF (iflag_alp_wk_cond .GT. 0.) THEN
 
-         CALL alpale_wk(dtime, cell_area, wake_k, wake_s, wake_dens, wake_fip_0, &
+         CALL alpale_wk(phys_tstep, cell_area, wake_k, wake_s, wake_dens, wake_fip_0, &
                         wake_fip)
        ELSE
@@ -3147,5 +3327,5 @@
                    (d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wake_k, 'the', abortphy) 
              ENDIF
-             call prt_enerbil('the',itap)
+             CALL prt_enerbil('the',itap)
           !
           ENDIF  ! (mod(iflag_pbl_split/10,10) .GE. 1)
@@ -3153,8 +3333,8 @@
           CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,  &
                              dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
-          call prt_enerbil('thermals',itap)
+          CALL prt_enerbil('thermals',itap)
           !
 !
-          CALL alpale_th( dtime, lmax_th, t_seri, cell_area,  &
+          CALL alpale_th( phys_tstep, lmax_th, t_seri, cell_area,  &
                           cin, s2, n2,  &
                           ale_bl_trig, ale_bl_stat, ale_bl,  &
@@ -3216,5 +3396,5 @@
           CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
                'ajsb',abortphy,flag_inhib_tend,itap,0)
-          call prt_enerbil('ajsb',itap)
+          CALL prt_enerbil('ajsb',itap)
           d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
           d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
@@ -3246,5 +3426,5 @@
     ENDIF
     !
-    CALL fisrtilp(dtime,paprs,pplay, &
+    CALL fisrtilp(phys_tstep,paprs,pplay, &
          t_seri, q_seri,ptconv,ratqs, &
          d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
@@ -3267,5 +3447,5 @@
     CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
          'lsc',abortphy,flag_inhib_tend,itap,0)
-    call prt_enerbil('lsc',itap)
+    CALL prt_enerbil('lsc',itap)
     rain_num(:)=0.
     DO k = 1, klev
@@ -3306,5 +3486,5 @@
                + snow_lsc(i))*cell_area(i)/REAL(klon)
        ENDDO
-       zx_t = zx_t/za*dtime
+       zx_t = zx_t/za*phys_tstep
        WRITE(lunout,*)"Precip=", zx_t
     ENDIF
@@ -3526,11 +3706,9 @@
        calday = REAL(days_elapsed + 1) + jH_cur
 
-       CALL chemtime(itap+itau_phy-1, date0, dtime, itap)
-       IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
-          CALL AEROSOL_METEO_CALC( &
-               calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
-               prfl,psfl,pctsrf,cell_area, &
-               latitude_deg,longitude_deg,u10m,v10m)
-       ENDIF
+       CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap)
+       CALL AEROSOL_METEO_CALC( &
+            calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
+            prfl,psfl,pctsrf,cell_area, &
+            latitude_deg,longitude_deg,u10m,v10m)
 
        zxsnow_dummy(:) = 0.0
@@ -3615,5 +3793,5 @@
 #else
                    !--climatologies or INCA aerosols
-                   CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, &
+                   CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, &
                         new_aod, flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
                         pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
@@ -3705,5 +3883,5 @@
              CALL readaerosolstrato1_rrtm(debut)
             ELSEIF (flag_aerosol_strat.EQ.2) THEN
-             CALL readaerosolstrato2_rrtm(debut)
+             CALL readaerosolstrato2_rrtm(debut, ok_volcan)
             ELSE
              abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
@@ -3717,4 +3895,6 @@
 #endif
           ENDIF
+       ELSE
+          tausum_aero(:,:,id_STRAT_phy) = 0. 
        ENDIF
 !
@@ -3895,9 +4075,13 @@
           RCFC11 = RCFC11_act
           RCFC12 = RCFC12_act
+          ! 
+          !--interactive CO2 in ppm from carbon cycle
+          IF (carbon_cycle_rad.AND..NOT.debut) THEN
+            RCO2=RCO2_glo
+          ENDIF
           !
           IF (prt_level .GE.10) THEN
              print *,' ->radlwsw, number 1 '
           ENDIF
-
           !
           CALL radlwsw &
@@ -3909,5 +4093,6 @@
                t_seri,q_seri,wo, &
                cldfrarad, cldemirad, cldtaurad, &
-               ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
+               ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie,  ok_volcan, &
+               flag_aerosol, &
                flag_aerosol_strat, flag_aer_feedback, &
                tau_aero, piz_aero, cg_aero, &
@@ -3920,4 +4105,5 @@
                ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
                heat,heat0,cool,cool0,albpla, &
+               heat_volc,cool_volc, &
                topsw,toplw,solsw,sollw, &
                sollwdown, &
@@ -3994,5 +4180,6 @@
                      t_seri,q_seri,wo, &
                      cldfrarad, cldemirad, cldtaurad, &
-                     ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
+                     ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie,  ok_volcan, &
+                     flag_aerosol, &
                      flag_aerosol_strat, flag_aer_feedback, &
                      tau_aero, piz_aero, cg_aero, &
@@ -4005,4 +4192,5 @@
                      ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
                      heatp,heat0p,coolp,cool0p,albplap, &
+                     heat_volc,cool_volc, &
                      topswp,toplwp,solswp,sollwp, &
                      sollwdownp, &
@@ -4072,14 +4260,14 @@
 
     DO k=1, klev
-       d_t_swr(:,k)=swradcorr(:)*heat(:,k)*dtime/RDAY
-       d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)*dtime/RDAY
-       d_t_lwr(:,k)=-cool(:,k)*dtime/RDAY
-       d_t_lw0(:,k)=-cool0(:,k)*dtime/RDAY
+       d_t_swr(:,k)=swradcorr(:)*heat(:,k)*phys_tstep/RDAY
+       d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)*phys_tstep/RDAY
+       d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY
+       d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY
     ENDDO
 
     CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
-    call prt_enerbil('SW',itap)
+    CALL prt_enerbil('SW',itap)
     CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
-    call prt_enerbil('LW',itap)
+    CALL prt_enerbil('LW',itap)
 
     !
@@ -4131,5 +4319,5 @@
        IF (ok_strato) THEN
 
-          CALL drag_noro_strato(0,klon,klev,dtime,paprs,pplay, &
+          CALL drag_noro_strato(0,klon,klev,phys_tstep,paprs,pplay, &
                zmea,zstd, zsig, zgam, zthe,zpic,zval, &
                igwd,idx,itest, &
@@ -4139,5 +4327,5 @@
 
        ELSE
-          CALL drag_noro(klon,klev,dtime,paprs,pplay, &
+          CALL drag_noro(klon,klev,phys_tstep,paprs,pplay, &
                zmea,zstd, zsig, zgam, zthe,zpic,zval, &
                igwd,idx,itest, &
@@ -4152,5 +4340,5 @@
        CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
             abortphy,flag_inhib_tend,itap,0)
-       call prt_enerbil('oro',itap)
+       CALL prt_enerbil('oro',itap)
        !----------------------------------------------------------------------
        !
@@ -4180,5 +4368,5 @@
        IF (ok_strato) THEN
 
-          CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
+          CALL lift_noro_strato(klon,klev,phys_tstep,paprs,pplay, &
                latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
                igwd,idx,itest, &
@@ -4188,5 +4376,5 @@
 
        ELSE
-          CALL lift_noro(klon,klev,dtime,paprs,pplay, &
+          CALL lift_noro(klon,klev,phys_tstep,paprs,pplay, &
                latitude_deg,zmea,zstd,zpic, &
                itest, &
@@ -4199,5 +4387,5 @@
        CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
             'lif', abortphy,flag_inhib_tend,itap,0)
-       call prt_enerbil('lif',itap)
+       CALL prt_enerbil('lif',itap)
     ENDIF ! fin de test sur ok_orolf
 
@@ -4208,5 +4396,5 @@
        du_gwd_hines=0.
        dv_gwd_hines=0.
-       CALL hines_gwd(klon, klev, dtime, paprs, pplay, latitude_deg, t_seri, &
+       CALL hines_gwd(klon, klev, phys_tstep, paprs, pplay, latitude_deg, t_seri, &
             u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
             du_gwd_hines, dv_gwd_hines)
@@ -4214,7 +4402,7 @@
        zvstr_gwd_hines=0.
        DO k = 1, klev
-          zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/dtime &
+          zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/phys_tstep &
                * (paprs(:, k)-paprs(:, k+1))/rg
-          zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/dtime &
+          zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/phys_tstep &
                * (paprs(:, k)-paprs(:, k+1))/rg
        ENDDO
@@ -4223,9 +4411,10 @@
        CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
             dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
-       call prt_enerbil('hin',itap)
+       CALL prt_enerbil('hin',itap)
     ENDIF
 
     IF (.not. ok_hines .and. ok_gwd_rando) then
-       CALL acama_GWD_rando(DTIME, pplay, latitude_deg, t_seri, u_seri, &
+       ! ym missing init for east_gwstress & west_gwstress -> added in phys_local_var_mod
+       CALL acama_GWD_rando(PHYS_TSTEP, pplay, latitude_deg, t_seri, u_seri, &
             v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
             dv_gwd_front, east_gwstress, west_gwstress)
@@ -4233,7 +4422,7 @@
        zvstr_gwd_front=0.
        DO k = 1, klev
-          zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/dtime &
+          zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/phys_tstep &
                * (paprs(:, k)-paprs(:, k+1))/rg
-          zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/dtime &
+          zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/phys_tstep &
                * (paprs(:, k)-paprs(:, k+1))/rg
        ENDDO
@@ -4241,20 +4430,20 @@
        CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
             paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
-       call prt_enerbil('front_gwd_rando',itap)
+       CALL prt_enerbil('front_gwd_rando',itap)
     ENDIF
 
     IF (ok_gwd_rando) THEN
-       CALL FLOTT_GWD_rando(DTIME, pplay, t_seri, u_seri, v_seri, &
+       CALL FLOTT_GWD_rando(PHYS_TSTEP, pplay, t_seri, u_seri, v_seri, &
             rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
             du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
        CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
             paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
-       call prt_enerbil('flott_gwd_rando',itap)
+       CALL prt_enerbil('flott_gwd_rando',itap)
        zustr_gwd_rando=0.
        zvstr_gwd_rando=0.
        DO k = 1, klev
-          zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/dtime &
+          zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/phys_tstep &
                * (paprs(:, k)-paprs(:, k+1))/rg
-          zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/dtime &
+          zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/phys_tstep &
                * (paprs(:, k)-paprs(:, k+1))/rg
        ENDDO
@@ -4276,7 +4465,7 @@
     DO k = 1, klev
        DO i = 1, klon
-          zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
+          zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/phys_tstep* &
                (paprs(i,k)-paprs(i,k+1))/rg
-          zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
+          zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/phys_tstep* &
                (paprs(i,k)-paprs(i,k+1))/rg
        ENDDO
@@ -4296,11 +4485,11 @@
     !IM cf. FLott END
     !DC Calcul de la tendance due au methane
-    IF(ok_qch4) THEN
+    IF (ok_qch4) THEN
        CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
        ! ajout de la tendance d'humidite due au methane
-       d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*dtime
+       d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*phys_tstep
        CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
             'q_ch4', abortphy,flag_inhib_tend,itap,0)
-       d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/dtime
+       d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/phys_tstep
     ENDIF
     !
@@ -4313,6 +4502,4 @@
 ! Inititialization
 !------------------
-
-    
 
        addtkeoro=0   
@@ -4326,13 +4513,11 @@
        alphatkeoro=min(max(0.,alphatkeoro),1.)
 
-       smallscales_tkeoro=.false.   
+       smallscales_tkeoro=.FALSE.   
        CALL getin_p('smallscales_tkeoro',smallscales_tkeoro) 
 
 
-        dtadd(:,:)=0.
-        duadd(:,:)=0.
-        dvadd(:,:)=0.
-
-
+       dtadd(:,:)=0.
+       duadd(:,:)=0.
+       dvadd(:,:)=0.
 
 ! Choices for addtkeoro:
@@ -4349,5 +4534,4 @@
 
 
-
   IF (addtkeoro .EQ. 1 ) THEN
 
@@ -4357,7 +4541,5 @@
   ELSE IF (addtkeoro .EQ. 2) THEN
 
-
-
-       IF (smallscales_tkeoro) THEN
+     IF (smallscales_tkeoro) THEN
        igwd=0
        DO i=1,klon
@@ -4382,13 +4564,10 @@
              igwd=igwd+1
              idx(igwd)=i
-          ENDIF
-       ENDDO
-
-       END IF
-
-
-
-
-       CALL drag_noro_strato(addtkeoro,klon,klev,dtime,paprs,pplay, &
+        ENDIF
+       ENDDO
+
+     ENDIF
+
+     CALL drag_noro_strato(addtkeoro,klon,klev,phys_tstep,paprs,pplay, &
                zmea,zstd, zsig, zgam, zthe,zpic,zval, &
                igwd,idx,itest, &
@@ -4397,13 +4576,12 @@
                d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
 
-            zustrdr(:)=0.
-            zvstrdr(:)=0.
-            zulow(:)=0.
-            zvlow(:)=0.
-
-            duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
-            dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
- END IF
-    
+     zustrdr(:)=0.
+     zvstrdr(:)=0.
+     zulow(:)=0.
+     zvlow(:)=0.
+
+     duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
+     dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
+  ENDIF
 
 
@@ -4416,9 +4594,7 @@
 
 
-
        ENDIF
 !      -----
 !===============================================================
-
 
 
@@ -4431,5 +4607,5 @@
        ! adeclarer 
 #ifdef CPP_COSP
-       IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
+       IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
 
           IF (prt_level .GE.10) THEN
@@ -4439,5 +4615,5 @@
           !       print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
           !     s        ref_liq,ref_ice
-          CALL phys_cosp(itap,dtime,freq_cosp, &
+          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, &
@@ -4462,5 +4638,5 @@
 
 #ifdef CPP_COSP2
-       IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
+       IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
 
           IF (prt_level .GE.10) THEN
@@ -4470,5 +4646,5 @@
                  print*,'Dans physiq.F avant appel '
           !     s        ref_liq,ref_ice
-          CALL phys_cosp2(itap,dtime,freq_cosp, &
+          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, &
@@ -4485,4 +4661,28 @@
 #endif
 
+#ifdef CPP_COSPV2
+       IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
+
+          IF (prt_level .GE.10) THEN
+             print*,'freq_cosp',freq_cosp
+          ENDIF
+          mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
+                 print*,'Dans physiq.F avant appel '
+          !     s        ref_liq,ref_ice
+          CALL lmdz_cosp_interface(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
+
     ENDIF  !ok_cosp
 
@@ -4492,5 +4692,5 @@
   IF (ok_airs) then
 
-  IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/dtime)).EQ.0) THEN
+  IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/phys_tstep)).EQ.0) THEN
      write(*,*) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
      CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
@@ -4547,5 +4747,5 @@
     CALL phytrac ( &
          itap,     days_elapsed+1,    jH_cur,   debut, &
-         lafin,    dtime,     u, v,     t, &
+         lafin,    phys_tstep,     u, v,     t, &
          paprs,    pplay,     pmfu,     pmfd, &
          pen_u,    pde_u,     pen_d,    pde_d, &
@@ -4582,5 +4782,5 @@
             cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
             frac_impa, frac_nucl, &
-            pphis,cell_area,dtime,itap, &
+            pphis,cell_area,phys_tstep,itap, &
             qx(:,:,ivap),da,phi,mp,upwd,dnwd)
 
@@ -4653,5 +4853,5 @@
 
        CALL chemhook_end ( &
-            dtime, &
+            phys_tstep, &
             pplay, &
             t_seri, &
@@ -4688,12 +4888,12 @@
     DO k = 1, klev
        DO i = 1, klon
-          d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
-          d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
-          d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
-          d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
-          d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
+          d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / phys_tstep
+          d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / phys_tstep
+          d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / phys_tstep
+          d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / phys_tstep
+          d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
           !CR: on ajoute le contenu en glace
           IF (nqo.eq.3) THEN
-             d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / dtime
+             d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
           ENDIF
        ENDDO
@@ -4707,6 +4907,6 @@
           DO  k = 1, klev
              DO  i = 1, klon
-                ! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
-                d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / dtime
+                ! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / phys_tstep
+                d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / phys_tstep
              ENDDO
           ENDDO
@@ -4877,5 +5077,5 @@
     CALL phys_output_write(itap, pdtphys, paprs, pphis,  &
          pplay, lmax_th, aerosol_couple,                 &
-         ok_ade, ok_aie, ivap, iliq, isol, new_aod,      & 
+         ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, new_aod,      & 
          ok_sync, ptconv, read_climoz, clevSTD,          &
          ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
@@ -4890,5 +5090,5 @@
 
 ! On remet des variables a .false. apres un premier appel
-    if (debut) then
+    IF (debut) THEN
 #ifdef CPP_XIOS
       swaero_diag=.FALSE.
@@ -4898,5 +5098,5 @@
 !      write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
 
-      IF (is_master) then
+      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. & 
@@ -4929,13 +5129,13 @@
            xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
            ok_4xCO2atm=.TRUE. 
-      endif
+      ENDIF
       !$OMP BARRIER
-      call bcast(swaero_diag)
-      call bcast(swaerofree_diag)
-      call bcast(dryaod_diag)
-      call bcast(ok_4xCO2atm)
+      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
-    endif
+    ENDIF
 
     !====================================================================
@@ -4962,19 +5162,22 @@
        !         write(97) u_seri,v_seri,t_seri,q_seri
        !         close(97)
-       !$OMP MASTER
-       IF (read_climoz >= 1) THEN
-          IF (is_mpi_root) THEN
-             CALL nf95_close(ncid_climoz)
-          ENDIF
-          DEALLOCATE(press_edg_climoz) ! pointer
-          DEALLOCATE(press_cen_climoz) ! pointer
-       ENDIF
-       !$OMP END MASTER
-       print *,' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
+      
+       IF (is_omp_master) THEN
+       
+         IF (read_climoz >= 1) THEN
+           IF (is_mpi_root) CALL nf95_close(ncid_climoz)
+            DEALLOCATE(press_edg_climoz) ! pointer
+            DEALLOCATE(press_cen_climoz) ! pointer
+         ENDIF
+       
+       ENDIF
+#ifdef CPP_XIOS
+       IF (is_omp_master) CALL xios_context_finalize
+#endif
+       WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
     ENDIF
 
     !      first=.false.
 
-
   END SUBROUTINE physiq
 
