Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/add_phys_tend_mod.F90
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/add_phys_tend_mod.F90	(revision 3055)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/add_phys_tend_mod.F90	(revision 3065)
@@ -1,2 +1,5 @@
+!
+! $Id$
+!
 !
 MODULE add_phys_tend_mod
@@ -82,5 +85,5 @@
 END SUBROUTINE add_pbl_tend
 !
-! $Id: add_phys_tend.F90 2611 2016-08-03 15:41:26Z jyg $
+! $Id$
 !
 SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,zdqi,paprs,text, &
@@ -98,5 +101,5 @@
 
 USE dimphy, ONLY: klon, klev
-USE phys_state_var_mod, ONLY : dtime
+USE phys_state_var_mod, ONLY : phys_tstep
 USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, q_seri, t_seri
 USE phys_state_var_mod, ONLY: ftsol
@@ -451,17 +454,17 @@
     ! ------------------------------------------------
 
-    d_qw_col(:) = (zqw_col(:,2)-zqw_col(:,1))/dtime
-    d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/dtime
-    d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/dtime
+    d_qw_col(:) = (zqw_col(:,2)-zqw_col(:,1))/phys_tstep
+    d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/phys_tstep
+    d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/phys_tstep
     d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:)
 
-    d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/dtime
-
-    d_h_dair_col(:) = (zh_dair_col(:,2)-zh_dair_col(:,1))/dtime
-    d_h_qw_col(:) = (zh_qw_col(:,2)-zh_qw_col(:,1))/dtime
-    d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/dtime
-    d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/dtime
-
-    d_h_col = (zh_col(:,2)-zh_col(:,1))/dtime 
+    d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/phys_tstep
+
+    d_h_dair_col(:) = (zh_dair_col(:,2)-zh_dair_col(:,1))/phys_tstep
+    d_h_qw_col(:) = (zh_qw_col(:,2)-zh_qw_col(:,1))/phys_tstep
+    d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/phys_tstep
+    d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/phys_tstep
+
+    d_h_col = (zh_col(:,2)-zh_col(:,1))/phys_tstep 
 
   end if ! end if (fl_ebil .GT. 0)
@@ -494,5 +497,5 @@
 !======================================================================
 
-USE phys_state_var_mod, ONLY : dtime, ftsol
+USE phys_state_var_mod, ONLY : phys_tstep, ftsol
 USE geometry_mod, ONLY: longitude_deg, latitude_deg
 USE print_control_mod, ONLY: prt_level
@@ -621,10 +624,10 @@
     ! ------------------------------------------------
 
-    d_qw_col(:) = (zqw_col(:,2)-zqw_col(:,1))/dtime
-    d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/dtime
-    d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/dtime
+    d_qw_col(:) = (zqw_col(:,2)-zqw_col(:,1))/phys_tstep
+    d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/phys_tstep
+    d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/phys_tstep
     d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:)
 
-    d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/dtime
+    d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/phys_tstep
 
    print *,'zdu ', zdu
@@ -632,10 +635,10 @@
    print *,'d_ek_col, zek_col(2), zek_col(1) ',d_ek_col(1), zek_col(1,2), zek_col(1,1)
 
-    d_h_dair_col(:) = (zh_dair_col(:,2)-zh_dair_col(:,1))/dtime
-    d_h_qw_col(:) = (zh_qw_col(:,2)-zh_qw_col(:,1))/dtime
-    d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/dtime
-    d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/dtime
-
-    d_h_col = (zh_col(:,2)-zh_col(:,1))/dtime 
+    d_h_dair_col(:) = (zh_dair_col(:,2)-zh_dair_col(:,1))/phys_tstep
+    d_h_qw_col(:) = (zh_qw_col(:,2)-zh_qw_col(:,1))/phys_tstep
+    d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/phys_tstep
+    d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/phys_tstep
+
+    d_h_col = (zh_col(:,2)-zh_col(:,1))/phys_tstep 
 
   end if ! end if (fl_ebil .GT. 0)
@@ -716,5 +719,5 @@
 
 USE dimphy, ONLY: klon, klev
-USE phys_state_var_mod, ONLY : dtime
+USE phys_state_var_mod, ONLY : phys_tstep
 USE phys_state_var_mod, ONLY : topsw, toplw, solsw, sollw, rain_con, snow_con
 USE geometry_mod, ONLY: longitude_deg, latitude_deg
Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/calcul_divers.h
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/calcul_divers.h	(revision 3055)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/calcul_divers.h	(revision 3065)
@@ -1,2 +1,5 @@
+!
+! $Id$
+!
 !
 ! $Header$
@@ -14,5 +17,5 @@
 
 ! Calcul fin de journee : total_rain, nday_rain
-      IF(MOD(itap,NINT(un_jour/dtime)).EQ.0) THEN
+      IF(MOD(itap,NINT(un_jour/phys_tstep)).EQ.0) THEN
 !        print*,'calcul nday_rain itap ',itap
          DO i = 1, klon
@@ -23,6 +26,6 @@
 
 ! Initialisation fin de mois
-      IF(MOD(itap-itapm1,NINT(mth_len*un_jour/dtime)).EQ.0) THEN
-        itapm1=itapm1+NINT(mth_len*un_jour/dtime)
+      IF(MOD(itap-itapm1,NINT(mth_len*un_jour/phys_tstep)).EQ.0) THEN
+        itapm1=itapm1+NINT(mth_len*un_jour/phys_tstep)
 !       print*,'initialisation itapm1 ',itapm1
       ENDIF
@@ -35,5 +38,5 @@
      t2m_max_mon=0.
   ENDIF
-  IF(MOD(itap,NINT(un_jour/dtime)).EQ.1) THEN
+  IF(MOD(itap,NINT(un_jour/phys_tstep)).EQ.1) THEN
      zt2m_min_mon=zt2m
      zt2m_max_mon=zt2m
@@ -45,5 +48,5 @@
      ENDDO
 !fin de journee
-  IF(MOD(itap,NINT(un_jour/dtime)).EQ.0) THEN
+  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 
Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/clesphys.h
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/clesphys.h	(revision 3055)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/clesphys.h	(revision 3065)
@@ -65,6 +65,4 @@
 !IM freq_outNMC : frequences de sortie fichiers niveaux de pression (histmthNMC, histdayNMC, histhfNMC)
 !IM freq_calNMC : frequences de calcul fis. hist*NMC.nc
-!IM pasphys : pas de temps de physique (secondes)
-       REAL pasphys
        LOGICAL ok_histNMC(3)
        INTEGER levout_histNMC(3)
@@ -106,5 +104,5 @@
      &     , min_wind_speed,f_gust_wk,f_gust_bl,f_qsat_oce,f_z0qh_oce   &
      &     , z0m_seaice,z0h_seaice                                      &
-     &     , pasphys            , freq_outNMC, freq_calNMC              &
+     &     , freq_outNMC, freq_calNMC                                   &
      &     , lonmin_ins, lonmax_ins, latmin_ins, latmax_ins             &
      &     , freq_ISCCP, ecrit_ISCCP, freq_COSP, freq_AIRS              &
Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/climoz_mod.F90
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/climoz_mod.F90	(revision 3065)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/climoz_mod.F90	(revision 3065)
@@ -0,0 +1,80 @@
+!
+! $Id$
+!
+! This module may contains all subroutine and function used
+! for ozone climatology
+! code re-design is needed
+
+MODULE climoz_mod
+
+
+
+
+
+CONTAINS
+
+
+  SUBROUTINE create_climoz(read_climoz)
+  USE mod_grid_phy_lmdz, ONLY :  grid_type, unstructured
+  USE regr_lat_time_climoz_m
+  USE mod_phys_lmdz_para
+#ifdef CPP_XIOS
+  USE XIOS
+#endif
+  IMPLICIT NONE
+  INTEGER, INTENT(IN)  :: read_climoz
+
+    !--- OZONE CLIMATOLOGY
+    IF (grid_type==unstructured) THEN
+#ifdef CPP_XIOS    
+      IF(read_climoz>=1) THEN 
+        IF (is_master) CALL regr_lat_time_climoz(read_climoz,.FALSE.)
+        IF (is_omp_master) THEN
+          CALL xios_set_field_attr( "tro3_reg",  enabled=.TRUE.)
+          CALL xios_set_field_attr( "tro3_read", enabled=.TRUE.)
+          IF (read_climoz==2) THEN
+            CALL xios_set_field_attr( "tro3_daylight_reg",  enabled=.TRUE.)
+            CALL xios_set_field_attr( "tro3_daylight_read", enabled=.TRUE.)
+          ENDIF
+        ENDIF
+      ENDIF
+#endif    
+    ENDIF
+
+  END SUBROUTINE create_climoz
+
+
+  
+  SUBROUTINE get_ozone_var(name,press_in_edg,paprs,v3)
+  USE dimphy
+  USE mod_phys_lmdz_para
+#ifdef CPP_XIOS
+  USE xios
+#endif
+  USE regr1_step_av_m, only: regr1_step_av
+  IMPLICIT NONE
+    CHARACTER(LEN=*), INTENT(in):: name(:) ! of the NetCDF variables      
+    REAL, INTENT(IN):: press_in_edg(:)    ! edges of pressure intervals for input data, in Pa, in strictly ascending order
+    REAL, INTENT(IN):: paprs(:, :) ! (klon, klev + 1)
+    REAL, INTENT(OUT):: v3(:, :, :) ! (klon, klev, size(name))
+    
+    REAL  :: v1_mpi(klon_mpi, size(press_in_edg) - 1, size(v3,3))
+    REAL  :: v1(klon, size(press_in_edg) - 1, size(v3,3))
+    INTEGER :: m,i
+
+#ifdef CPP_XIOS    
+    DO m=1,size(name)
+      IF (is_omp_master) CALL xios_recv_field(name(m),v1_mpi(:,:,m))
+      CALL scatter_omp(v1_mpi,v1)
+      ! Regrid in pressure at each horizontal position:
+      DO i = 1, klon
+         v3(i, klev:1:-1, m) = regr1_step_av(v1(i,:,m), press_in_edg, &
+              paprs(i, klev+1:1:-1))
+         ! (invert order of indices because "paprs" is in descending order)
+      END DO
+    ENDDO
+#endif
+  END SUBROUTINE get_ozone_var  
+  
+
+END MODULE climoz_mod
Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/conf_phys_m.F90
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/conf_phys_m.F90	(revision 3055)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/conf_phys_m.F90	(revision 3065)
@@ -1,4 +1,4 @@
 !
-! $Id: conf_phys.F90 1668 2012-10-12 10:47:37Z idelkadi $
+! $Id$
 !
 !
@@ -29,4 +29,5 @@
     USE mod_grid_phy_lmdz, ONLY: klon_glo
     USE print_control_mod, ONLY: lunout
+    USE phys_state_var_mod, ONLY: phys_tstep
 
     include "conema3.h"
@@ -1692,10 +1693,10 @@
     !Config Desc = freq_calNMC(2) = frequence de calcul fichiers histdayNMC
     !Config Desc = freq_calNMC(3) = frequence de calcul fichiers histhfNMC
-    !Config Def  = pasphys
-    !Config Help =
-    !
-    freq_calNMC_omp(1) = pasphys
-    freq_calNMC_omp(2) = pasphys
-    freq_calNMC_omp(3) = pasphys
+    !Config Def  = phys_tstep
+    !Config Help =
+    !
+    freq_calNMC_omp(1) = phys_tstep
+    freq_calNMC_omp(2) = phys_tstep
+    freq_calNMC_omp(3) = phys_tstep
     CALL getin('freq_calNMC',freq_calNMC_omp)
     !
Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/grid_noro_m.F90
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/grid_noro_m.F90	(revision 3055)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/grid_noro_m.F90	(revision 3065)
@@ -1,2 +1,5 @@
+!
+! $Id$
+!
 MODULE grid_noro_m
 !
@@ -334,6 +337,4 @@
   imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
   jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
-! IF(imar/=iim)   CALL abort_gcm(TRIM(modname),'imar/=iim'  ,1)
-! IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
   iext=imdp/10
   xpi = ACOS(-1.)
Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/ini_undefSTD.F90
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/ini_undefSTD.F90	(revision 3055)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/ini_undefSTD.F90	(revision 3065)
@@ -55,5 +55,5 @@
 
     IF (n==1 .AND. itap-itapm1==1 .OR. n>1 .AND. mod(itap,nint( &
-        freq_outnmc(n)/dtime))==1) THEN
+        freq_outnmc(n)/phys_tstep))==1) THEN
       ! print*,'ini_undefSTD n itap',n,itap
       DO k = 1, nlevstd
Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/moy_undefSTD.F90
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/moy_undefSTD.F90	(revision 3055)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/moy_undefSTD.F90	(revision 3065)
@@ -75,5 +75,5 @@
 
     IF (n==1 .AND. itap==itapm1 .OR. n>1 .AND. mod(itap,nint(freq_outnmc(n)/ & 
-        dtime))==0) THEN
+        phys_tstep))==0) THEN
 
       ! print*,'moy_undefSTD n itap itapm1',n,itap,itapm1
@@ -140,5 +140,5 @@
         END DO !i
       END DO !k
-    END IF !MOD(itap,NINT(freq_outNMC(n)/dtime)).EQ.0
+    END IF !MOD(itap,NINT(freq_outNMC(n)/phys_tstep)).EQ.0
 
   END DO !n
Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/phyetat0.F90
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/phyetat0.F90	(revision 3055)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/phyetat0.F90	(revision 3065)
@@ -9,5 +9,5 @@
   USE pbl_surface_mod,  ONLY : pbl_surface_init
   USE surface_data,     ONLY : type_ocean, version_ocean
-  USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, dtime, &
+  USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, phys_tstep, &
        qsol, fevap, z0m, z0h, agesno, &
        du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
@@ -466,5 +466,5 @@
 
   IF ( type_ocean == 'slab' ) THEN
-      CALL ocean_slab_init(dtime, pctsrf)
+      CALL ocean_slab_init(phys_tstep, pctsrf)
       IF (nslay.EQ.1) THEN
         found=phyetat0_get(1,tslab,"tslab01","tslab",0.)
@@ -515,5 +515,5 @@
   ! Initialize module ocean_cpl_mod for the case of coupled ocean
   IF ( type_ocean == 'couple' ) THEN
-     CALL ocean_cpl_init(dtime, longitude_deg, latitude_deg)
+     CALL ocean_cpl_init(phys_tstep, longitude_deg, latitude_deg)
   ENDIF
 
Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/phys_state_var_mod.F90	(revision 3055)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/phys_state_var_mod.F90	(revision 3065)
@@ -23,6 +23,6 @@
 !$OMP THREADPRIVATE(cvpas)
 !$OMP THREADPRIVATE(wkpas)
-      REAL, SAVE :: dtime, solaire_etat0
-!$OMP THREADPRIVATE(dtime, solaire_etat0)
+      REAL, SAVE :: phys_tstep, solaire_etat0
+!$OMP THREADPRIVATE(phys_tstep, solaire_etat0)
 
       REAL, ALLOCATABLE, SAVE :: pctsrf(:,:)
Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/physiq_mod.F90	(revision 3055)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/physiq_mod.F90	(revision 3065)
@@ -1113,4 +1113,5 @@
     !albedo SB >>>
     real,dimension(6),save :: SFRWL
+!$OMP THREADPRIVATE(SFRWL)
     !albedo SB <<<
 
@@ -1165,5 +1166,4 @@
        CALL phys_local_var_init
        !
-       pasphys=pdtphys
        !     appel a la lecture du run.def physique
        CALL conf_phys(ok_journe, ok_mensuel, &
@@ -1354,14 +1354,13 @@
        ! 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
+       phys_tstep=NINT(pdtphys)
+       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 ', &
@@ -1373,9 +1372,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 = NINT( 86400./dtime)/nbapp_cv
+       IF (MOD(NINT(86400./phys_tstep),nbapp_cv).EQ.0) THEN
+          cvpas = NINT( 86400./phys_tstep)/nbapp_cv
        print *,'physiq, cvpas ',cvpas
        ELSE 
@@ -1388,6 +1387,6 @@
           call abort_physic(modname,abort_message,1)
        ENDIF
-       IF (MOD(NINT(86400./dtime),nbapp_wk).EQ.0) THEN
-          wkpas = NINT( 86400./dtime)/nbapp_wk
+       IF (MOD(NINT(86400./phys_tstep),nbapp_wk).EQ.0) THEN
+          wkpas = NINT( 86400./phys_tstep)/nbapp_wk
        print *,'physiq, wkpas ',wkpas
        ELSE 
@@ -1431,11 +1430,11 @@
             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,  &
@@ -1451,5 +1450,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"
@@ -1597,5 +1596,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
@@ -1649,5 +1648,5 @@
        CALL phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
             iGCM,jGCM,lonGCM,latGCM, &
-            jjmp1,nlevSTD,clevSTD,rlevSTD, dtime,ok_veget, &
+            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, &
@@ -1669,5 +1668,5 @@
 
 #ifndef CPP_XIOS
-       CALL ini_paramLMDZ_phy(dtime,nid_ctesGCM)
+       CALL ini_paramLMDZ_phy(phys_tstep,nid_ctesGCM)
 #endif
 
@@ -1810,5 +1809,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)
@@ -1935,20 +1934,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
@@ -2126,5 +2125,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)
@@ -2143,6 +2142,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)
@@ -2150,5 +2149,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, &
@@ -2225,5 +2224,5 @@
 
        CALL pbl_surface(  &
-            dtime,     date0,     itap,    days_elapsed+1, &
+            phys_tstep,     date0,     itap,    days_elapsed+1, &
             debut,     lafin, &
             longitude_deg, latitude_deg, rugoro,  zrmu0,      &
@@ -2392,7 +2391,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
@@ -2432,9 +2431,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, &
@@ -2512,5 +2511,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, &
@@ -2549,5 +2548,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, &
@@ -2588,6 +2587,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)
@@ -2601,5 +2600,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, &
@@ -2725,5 +2724,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
@@ -2739,5 +2738,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
@@ -2776,6 +2775,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
@@ -2785,7 +2784,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
@@ -2798,7 +2797,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
@@ -2808,5 +2807,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,  &
@@ -2849,5 +2848,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
@@ -2984,5 +2983,5 @@
           !
 !
-          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,  &
@@ -3074,5 +3073,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, &
@@ -3134,5 +3133,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
@@ -3354,5 +3353,5 @@
        calday = REAL(days_elapsed + 1) + jH_cur
 
-       CALL chemtime(itap+itau_phy-1, date0, dtime, itap)
+       CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap)
        IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
           CALL AEROSOL_METEO_CALC( &
@@ -3877,8 +3876,8 @@
 
     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
 
@@ -3936,5 +3935,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, &
@@ -3944,5 +3943,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, &
@@ -3985,5 +3984,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, &
@@ -3993,5 +3992,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, &
@@ -4013,5 +4012,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)
@@ -4019,7 +4018,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
@@ -4032,5 +4031,5 @@
 
     IF (.not. ok_hines .and. ok_gwd_rando) then
-       CALL acama_GWD_rando(DTIME, pplay, latitude_deg, t_seri, u_seri, &
+       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)
@@ -4038,7 +4037,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
@@ -4050,5 +4049,5 @@
 
     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)
@@ -4059,7 +4058,7 @@
        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
@@ -4081,7 +4080,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
@@ -4104,8 +4103,8 @@
        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
     !
@@ -4195,5 +4194,5 @@
 
 
-       CALL drag_noro_strato(addtkeoro,klon,klev,dtime,paprs,pplay, &
+       CALL drag_noro_strato(addtkeoro,klon,klev,phys_tstep,paprs,pplay, &
                zmea,zstd, zsig, zgam, zthe,zpic,zval, &
                igwd,idx,itest, &
@@ -4236,5 +4235,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
@@ -4244,5 +4243,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, &
@@ -4273,5 +4272,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,&
@@ -4328,5 +4327,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, &
@@ -4363,5 +4362,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)
 
@@ -4434,5 +4433,5 @@
 
        CALL chemhook_end ( &
-            dtime, &
+            phys_tstep, &
             pplay, &
             t_seri, &
@@ -4469,12 +4468,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
@@ -4488,6 +4487,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
Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/regr_lat_time_climoz_m.F90
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/regr_lat_time_climoz_m.F90	(revision 3065)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/regr_lat_time_climoz_m.F90	(revision 3065)
@@ -0,0 +1,496 @@
+! $Id$
+module regr_lat_time_climoz_m
+
+  ! Author: Lionel GUEZ
+
+  implicit none
+
+  private
+  public regr_lat_time_climoz
+
+contains
+
+  subroutine regr_lat_time_climoz(read_climoz, regr_lat)
+
+    ! "regr_lat_time_climoz" stands for "regrid latitude time
+    ! climatology ozone".
+
+    ! This procedure reads a climatology of ozone from a NetCDF file,
+    ! regrids it in latitude and time, and writes the regridded field
+    ! to a new NetCDF file.
+
+    ! The input field depends on time, pressure level and latitude.
+
+    ! If the input field has missing values, they must be signaled by
+    ! the "missing_value" attribute.
+
+    ! We assume that the input field is a step function of latitude
+    ! and that the input latitude coordinate gives the centers of steps.
+    ! Regridding in latitude is made by averaging, with a cosine of
+    ! latitude factor.
+    ! The target LMDZ latitude grid is the "scalar" grid: "rlatu".
+    ! The values of "rlatu" are taken to be the centers of intervals.
+
+    ! We assume that in the input file:
+
+    ! -- Latitude is in degrees.
+
+    ! -- Latitude and pressure are strictly monotonic (as all NetCDF
+    ! coordinate variables should be).
+
+    ! -- The time coordinate is in ascending order (even though we do
+    ! not use its values).
+    ! The input file may contain either values for 12 months or values
+    ! for 14 months.
+    ! If there are 14 months then we assume that we have (in that order):
+    ! December, January, February, ..., November, December, January
+
+    ! -- Missing values are contiguous, at the bottom of
+    ! the vertical domain and at the latitudinal boundaries.
+
+    ! If values are all missing at a given latitude and date, then we
+    ! replace those missing values by values at the closest latitude,
+    ! equatorward, with valid values.
+    ! Then, at each latitude and each date, the missing values are replaced
+    ! by the lowest valid value above missing values.
+
+    ! Regridding in time is by linear interpolation.
+    ! Monthly values are processed to get daily values, on the basis
+    ! of a 360-day calendar.
+    ! If there are 14 months, we use the first December value to
+    ! interpolate values between January 1st and mid-January.
+    ! We use the last January value to interpolate values between
+    ! mid-December and end of December.
+    ! If there are only 12 months in the input file then we assume
+    ! periodicity for interpolation at the beginning and at the end of the
+    ! year.
+    use mod_grid_phy_lmdz, ONLY : nbp_lat, grid_type, unstructured
+    use regr1_step_av_m, only: regr1_step_av
+    use regr3_lint_m, only: regr3_lint
+    use netcdf95, only: handle_err, nf95_close, nf95_get_att, nf95_gw_var, &
+         nf95_inq_dimid, nf95_inq_varid, nf95_inquire_dimension, nf95_open, &
+         nf95_put_var
+    use netcdf, only: nf90_get_att, nf90_get_var, nf90_noerr, nf90_nowrite
+    use assert_m, only: assert
+    use regular_lonlat_mod, only : boundslat_reg, south
+    use nrtype, only: pi
+    use regular_lonlat_mod, only : lat_reg
+    implicit none
+    integer, intent(in):: read_climoz ! read ozone climatology
+    ! Allowed values are 1 and 2
+    ! 1: read a single ozone climatology that will be used day and night
+    ! 2: read two ozone climatologies, the average day and night
+    ! climatology and the daylight climatology
+
+    ! Variables local to the procedure:
+    LOGICAL, intent(in), OPTIONAL :: regr_lat
+
+    integer n_plev ! number of pressure levels in the input data
+    integer n_lat ! number of latitudes in the input data
+    integer n_month ! number of months in the input data
+    integer n_lat_out ! number of latitudes in the output data
+
+    real, pointer:: latitude(:)
+    ! (of input data, converted to rad, sorted in strictly ascending order)
+
+    real, allocatable:: lat_in_edg(:)
+    ! (edges of latitude intervals for input data, in rad, in strictly
+    ! ascending order)
+
+    real, pointer:: plev(:)
+    ! pressure levels of input data, sorted in strictly ascending
+    ! order, converted to hPa
+
+    logical desc_lat ! latitude in descending order in the input file
+    logical desc_plev ! pressure levels in descending order in the input file
+
+    real, allocatable:: o3_in(:, :, :, :)
+    ! (n_lat, n_plev, n_month, read_climoz)
+    ! ozone climatologies from the input file
+    ! "o3_in(j, k, :, :)" is at latitude "latitude(j)" and pressure
+    ! level "plev(k)".
+    ! Third dimension is month index, first value may be December or January.
+    ! "o3_in(:, :, :, 1)" is for the day- night average, "o3_in(:, :, :, 2)"
+    ! is for daylight.
+
+    real missing_value
+
+    real, allocatable:: o3_regr_lat(:, :, :, :)
+    ! (nbp_lat, n_plev, 0:13, read_climoz)
+    ! mean of "o3_in" over a latitude interval of LMDZ
+    ! First dimension is latitude interval.
+    ! The latitude interval for "o3_regr_lat(j,:, :, :)" contains "rlatu(j)".
+    ! If "j" is between 2 and "nbp_lat-1" then the interval is:
+    ! [rlatv(j), rlatv(j-1)]
+    ! If "j" is 1 or "nbp_lat" then the interval is:
+    ! [rlatv(1), pi / 2]
+    ! or:
+    ! [- pi / 2, rlatv(nbp_lat-1)]
+    ! respectively.
+    ! "o3_regr_lat(:, k, :, :)" is for pressure level "plev(k)".
+    ! Third dimension is month number, 1 for January.
+    ! "o3_regr_lat(:, :, :, 1)" is average day and night,
+    ! "o3_regr_lat(:, :, :, 2)" is for daylight.
+
+    real, allocatable:: o3_out(:, :, :, :)
+    ! (nbp_lat, n_plev, 360, read_climoz)
+    ! regridded ozone climatology
+    ! "o3_out(j, k, l, :)" is at latitude "rlatu(j)", pressure
+    ! level "plev(k)" and date "January 1st 0h" + "tmidday(l)", in a
+    ! 360-day calendar.
+    ! "o3_out(:, :, :, 1)" is average day and night,
+    ! "o3_out(:, :, :, 2)" is for daylight.
+
+    integer j, k, l,m
+
+    ! For NetCDF:
+    integer ncid_in, ncid_out ! IDs for input and output files
+    integer varid_plev, varid_time, varid, ncerr, dimid
+    character(len=80) press_unit ! pressure unit
+
+    integer varid_in(read_climoz), varid_out(read_climoz)
+    ! index 1 is for average ozone day and night, index 2 is for
+    ! daylight ozone.
+
+    real, parameter:: tmidmonth(0:13) = (/(-15. + 30. * l, l = 0, 13)/)
+    ! (time to middle of month, in days since January 1st 0h, in a
+    ! 360-day calendar)
+    ! (We add values -15 and 375 so that, for example, day 3 of the year is
+    ! interpolated between the December and the January value.)
+
+    real, parameter:: tmidday(360) = (/(l + 0.5, l = 0, 359)/)
+    ! (time to middle of day, in days since January 1st 0h, in a
+    ! 360-day calendar)
+    logical :: regr_lat_
+    REAL,ALLOCATABLE :: lat_out(:)
+    
+    !---------------------------------
+    IF (PRESENT(regr_lat)) THEN
+      regr_lat_=regr_lat
+    ELSE
+      regr_lat_=.TRUE.
+    ENDIF
+
+    print *, "Call sequence information: regr_lat_time_climoz"
+    call assert(read_climoz == 1 .or. read_climoz == 2, "regr_lat_time_climoz")
+
+    call nf95_open("climoz.nc", nf90_nowrite, ncid_in)
+
+    ! Get coordinates from the input file:
+
+    call nf95_inq_varid(ncid_in, "latitude", varid)
+    call nf95_gw_var(ncid_in, varid, latitude)
+    ! Convert from degrees to rad, because we will take the sine of latitude:
+    latitude = latitude / 180. * pi
+    n_lat = size(latitude)
+    ! We need to supply the latitudes to "regr1_step_av" in
+    ! ascending order, so invert order if necessary:
+    desc_lat = latitude(1) > latitude(n_lat)
+    if (desc_lat) latitude = latitude(n_lat:1:-1)
+
+    ! Compute edges of latitude intervals:
+    allocate(lat_in_edg(n_lat + 1))
+    lat_in_edg(1) = - pi / 2
+    forall (j = 2:n_lat) lat_in_edg(j) = (latitude(j - 1) + latitude(j)) / 2
+    lat_in_edg(n_lat + 1) = pi / 2
+
+    call nf95_inq_varid(ncid_in, "plev", varid)
+    call nf95_gw_var(ncid_in, varid, plev)
+    n_plev = size(plev)
+    ! We only need the pressure coordinate to copy it to the output file.
+    ! The program "gcm" will assume that pressure levels are in
+    ! ascending order in the regridded climatology so invert order if
+    ! necessary:
+    desc_plev = plev(1) > plev(n_plev)
+    if (desc_plev) plev = plev(n_plev:1:-1)
+    call nf95_get_att(ncid_in, varid, "units", press_unit)
+    if (press_unit == "Pa") then
+       ! Convert to hPa:
+       plev = plev / 100.
+    elseif (press_unit /= "hPa") then
+       print *, "regr_lat_time_climoz: the only recognized units are Pa " &
+            // "and hPa."
+       stop 1
+    end if
+
+    IF (regr_lat_) THEN 
+      n_lat_out=nbp_lat
+      ALLOCATE(lat_out(n_lat_out))
+      lat_out=lat_reg
+    ELSE 
+      n_lat_out=n_lat
+      ALLOCATE(lat_out(n_lat_out))
+      lat_out=latitude
+    ENDIF
+
+    ! Create the output file and get the variable IDs:
+    call prepare_out(ncid_in, n_lat_out, lat_out, n_plev, ncid_out, varid_out, varid_plev, &
+         varid_time)
+
+    deallocate(latitude) ! pointer
+
+    ! Write remaining coordinate variables:
+    call nf95_put_var(ncid_out, varid_plev, plev)
+    call nf95_put_var(ncid_out, varid_time, tmidday)
+
+    deallocate(plev) ! pointer
+
+    ! Get the  number of months:
+    call nf95_inq_dimid(ncid_in, "time", dimid)
+    call nf95_inquire_dimension(ncid_in, dimid, nclen=n_month)
+
+    allocate(o3_in(n_lat, n_plev, n_month, read_climoz))
+
+    call nf95_inq_varid(ncid_in, "tro3", varid_in(1))
+    ncerr = nf90_get_var(ncid_in, varid_in(1), o3_in(:, :, :, 1))
+    call handle_err("regr_lat_time_climoz nf90_get_var tro3", ncerr, ncid_in)
+
+    if (read_climoz == 2) then
+       call nf95_inq_varid(ncid_in, "tro3_daylight", varid_in(2))
+       ncerr = nf90_get_var(ncid_in, varid_in(2), o3_in(:, :, :, 2))
+       call handle_err("regr_lat_time_climoz nf90_get_var tro3_daylight", &
+            ncerr, ncid_in, varid_in(2))
+    end if
+
+    if (desc_lat) o3_in = o3_in(n_lat:1:-1, :, :, :)
+    if (desc_plev) o3_in = o3_in(:, n_plev:1:-1, :, :)
+
+    do m = 1, read_climoz
+       ncerr = nf90_get_att(ncid_in, varid_in(m), "missing_value", &
+            missing_value)
+       if (ncerr == nf90_noerr) then
+          do l = 1, n_month
+             ! Take care of latitudes where values are all missing:
+
+             ! Next to the south pole:
+             j = 1
+             do while (o3_in(j, 1, l, m) == missing_value)
+                j = j + 1
+             end do
+             if (j > 1) o3_in(:j-1, :, l, m) = &
+                  spread(o3_in(j, :, l, m), dim=1, ncopies=j-1)
+             
+             ! Next to the north pole:
+             j = n_lat
+             do while (o3_in(j, 1, l, m) == missing_value)
+                j = j - 1
+             end do
+             if (j < n_lat) o3_in(j+1:, :, l, m) = &
+                  spread(o3_in(j, :, l, m), dim=1, ncopies=n_lat-j)
+
+             ! Take care of missing values at high pressure:
+             do j = 1, n_lat
+                ! Find missing values, starting from top of atmosphere
+                ! and going down.
+                ! We have already taken care of latitudes full of
+                ! missing values so the highest level has a valid value.
+                k = 2
+                do while  (o3_in(j, k, l, m) /= missing_value .and. k < n_plev)
+                   k = k + 1
+                end do
+                ! Replace missing values with the valid value at the
+                ! lowest level above missing values:
+                if (o3_in(j, k, l, m) == missing_value) &
+                     o3_in(j, k:n_plev, l, m) = o3_in(j, k-1, l, m)
+             end do
+          end do
+       else
+          print *, "regr_lat_time_climoz: field ", m, &
+               ", no missing value attribute"
+       end if
+    end do
+
+    call nf95_close(ncid_in)
+
+    allocate(o3_out(n_lat_out, n_plev, 360, read_climoz))
+    IF (regr_lat_) THEN
+      allocate(o3_regr_lat(nbp_lat, n_plev, 0:13, read_climoz))
+
+      ! Regrid in latitude:
+      ! We average with respect to sine of latitude, which is
+      ! equivalent to weighting by cosine of latitude:
+      if (n_month == 12) then
+         print *, &
+              "Found 12 months in ozone climatologies, assuming periodicity..."
+         o3_regr_lat(nbp_lat:1:-1, :, 1:12, :) = regr1_step_av(o3_in, &
+              xs=sin(lat_in_edg), xt=sin((/- pi / 2, boundslat_reg(nbp_lat-1:1:-1,south), pi / 2/)))
+         ! (invert order of indices in "o3_regr_lat" because "rlatu" is
+         ! in descending order)
+
+         ! Duplicate January and December values, in preparation of time
+         ! interpolation:
+         o3_regr_lat(:, :, 0, :) = o3_regr_lat(:, :, 12, :)
+         o3_regr_lat(:, :, 13, :) = o3_regr_lat(:, :, 1, :)
+      else
+         print *, "Using 14 months in ozone climatologies..."
+         o3_regr_lat(nbp_lat:1:-1, :, :, :) = regr1_step_av(o3_in, &
+              xs=sin(lat_in_edg), xt=sin((/- pi / 2, boundslat_reg(nbp_lat-1:1:-1,south), pi / 2/)))
+         ! (invert order of indices in "o3_regr_lat" because "rlatu" is
+         ! in descending order)
+      end if
+      ! Regrid in time by linear interpolation:
+      o3_out = regr3_lint(o3_in, tmidmonth, tmidday)
+    ELSE
+      ! Regrid in time by linear interpolation:
+      o3_out = regr3_lint(o3_in, tmidmonth, tmidday)
+    ENDIF
+
+    ! Write to file:
+
+    do m = 1, read_climoz
+       IF (grid_type==unstructured) THEN
+         ! Doing spatial interpolation from XIOS need to have some point in longitude
+         ! waiting zonal mean operation from XIOS
+         call nf95_put_var(ncid_out, varid_out(m), SPREAD(o3_out(n_lat_out:1:-1, :, :, m),1,4))
+       ELSE
+         call nf95_put_var(ncid_out, varid_out(m), o3_out(n_lat_out:1:-1, :, :, m))
+       ENDIF
+       
+       ! (The order of "rlatu" is inverted in the output file)
+    end do
+
+    call nf95_close(ncid_out)
+
+  end subroutine regr_lat_time_climoz
+
+  !********************************************
+
+  subroutine prepare_out(ncid_in, n_lat, lat, n_plev, ncid_out, varid_out, varid_plev, &
+       varid_time)
+
+    ! This subroutine creates the NetCDF output file, defines
+    ! dimensions and variables, and writes one of the coordinate variables.
+
+    use netcdf95, only: nf95_create, nf95_def_dim, nf95_def_var, &
+         nf95_put_att, nf95_enddef, nf95_copy_att, nf95_put_var
+    use netcdf, only: nf90_clobber,nf90_netcdf4, nf90_float, nf90_global, nf90_unlimited
+    use nrtype, only: pi
+    use mod_grid_phy_lmdz, ONLY :  grid_type, unstructured
+
+    integer, intent(in):: ncid_in, n_lat, n_plev
+    real, intent(in):: lat(:)
+    integer, intent(out):: ncid_out, varid_plev, varid_time
+
+    integer, intent(out):: varid_out(:) ! dim(1 or 2)
+    ! "varid_out(1)" is for average ozone day and night,
+    ! "varid_out(2)" is for daylight ozone.
+
+    ! Variables local to the procedure:
+
+    integer ncerr
+    integer dimid_rlatu, dimid_rlonv, dimid_plev, dimid_time
+    integer varid_rlonv, varid_rlatu
+
+    !---------------------------
+
+    print *, "Call sequence information: prepare_out"
+
+    call nf95_create("climoz_LMDZ.nc", NF90_CLOBBER + NF90_NETCDF4, ncid_out)
+
+    ! Dimensions:
+    call nf95_def_dim(ncid_out, "time",  nf90_unlimited, dimid_time)
+    call nf95_def_dim(ncid_out, "plev", n_plev, dimid_plev)
+    call nf95_def_dim(ncid_out, "rlatu", n_lat, dimid_rlatu)
+    if (grid_type==unstructured) call nf95_def_dim(ncid_out, "rlonv", 4, dimid_rlonv)
+
+    ! Define coordinate variables:
+
+    call nf95_def_var(ncid_out, "time", nf90_float, dimid_time, varid_time)
+    call nf95_put_att(ncid_out, varid_time, "units", "days since 2000-1-1")
+    call nf95_put_att(ncid_out, varid_time, "calendar", "360_day")
+    call nf95_put_att(ncid_out, varid_time, "standard_name", "time")
+
+    call nf95_def_var(ncid_out, "plev", nf90_float, dimid_plev, varid_plev)
+    call nf95_put_att(ncid_out, varid_plev, "units", "millibar")
+    call nf95_put_att(ncid_out, varid_plev, "standard_name", "air_pressure")
+    call nf95_put_att(ncid_out, varid_plev, "long_name", "air pressure")
+
+    call nf95_def_var(ncid_out, "rlatu", nf90_float, dimid_rlatu, varid_rlatu)
+    call nf95_put_att(ncid_out, varid_rlatu, "units", "degrees_north")
+    call nf95_put_att(ncid_out, varid_rlatu, "standard_name", "latitude")
+
+    if (grid_type==unstructured) then
+      call nf95_def_var(ncid_out, "rlonv", nf90_float, dimid_rlonv, varid_rlonv)
+      call nf95_put_att(ncid_out, varid_rlonv, "units", "degrees_east")
+      call nf95_put_att(ncid_out, varid_rlonv, "standard_name", "longitude")
+    endif
+    
+    ! Define the primary variables:
+
+    if (grid_type==unstructured) then
+        call nf95_def_var(ncid_out, "tro3", nf90_float, &
+         (/dimid_rlonv, dimid_rlatu, dimid_plev, dimid_time/), varid_out(1))
+    else
+        call nf95_def_var(ncid_out, "tro3", nf90_float, &
+         (/dimid_rlatu, dimid_plev, dimid_time/), varid_out(1))
+    endif  
+    
+    call nf95_put_att(ncid_out, varid_out(1), "long_name", &
+         "ozone mole fraction")
+    call nf95_put_att(ncid_out, varid_out(1), "standard_name", &
+         "mole_fraction_of_ozone_in_air")
+
+    if (size(varid_out) == 2) then
+      if (grid_type==unstructured) then
+         call nf95_def_var(ncid_out, "tro3_daylight", nf90_float, &
+             (/dimid_rlonv, dimid_rlatu, dimid_plev, dimid_time/), varid_out(2))
+      else
+         call nf95_def_var(ncid_out, "tro3_daylight", nf90_float, &
+             (/ dimid_rlatu, dimid_plev, dimid_time/), varid_out(2))
+      endif
+       
+       call nf95_put_att(ncid_out, varid_out(2), "long_name", &
+            "ozone mole fraction in daylight")
+    end if
+
+    ! Global attributes:
+
+    ! The following commands, copying attributes, may fail.
+    ! That is OK.
+    ! It should just mean that the attribute is not defined in the input file.
+
+    call nf95_copy_att(ncid_in, nf90_global, "Conventions", ncid_out, &
+         nf90_global, ncerr)
+    call handle_err_copy_att("Conventions")
+
+    call nf95_copy_att(ncid_in, nf90_global, "title", ncid_out, nf90_global, &
+         ncerr)
+    call handle_err_copy_att("title")
+
+    call nf95_copy_att(ncid_in, nf90_global, "institution", ncid_out, &
+         nf90_global, ncerr)
+    call handle_err_copy_att("institution")
+
+    call nf95_copy_att(ncid_in, nf90_global, "source", ncid_out, nf90_global, &
+         ncerr)
+    call handle_err_copy_att("source")
+
+    call nf95_put_att(ncid_out, nf90_global, "comment", "Regridded for LMDZ")
+
+    call nf95_enddef(ncid_out)
+
+    ! Write one of the coordinate variables:
+    call nf95_put_var(ncid_out, varid_rlatu, lat(n_lat:1:-1) / pi * 180.)
+    ! (convert from rad to degrees and sort in ascending order)
+    if (grid_type==unstructured)  call nf95_put_var(ncid_out, varid_rlonv, (/ 0., 90., 180.,270. /))
+
+  contains
+
+    subroutine handle_err_copy_att(att_name)
+
+      use netcdf, only: nf90_noerr, nf90_strerror
+
+      character(len=*), intent(in):: att_name
+
+      !----------------------------------------
+
+      if (ncerr /= nf90_noerr) then
+         print *, "regr_lat_time_climoz_m prepare_out nf95_copy_att " &
+              // att_name // " -- " // trim(nf90_strerror(ncerr))
+      end if
+
+    end subroutine handle_err_copy_att
+
+  end subroutine prepare_out
+
+end module regr_lat_time_climoz_m
Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/undefSTD.F90
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/undefSTD.F90	(revision 3055)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/undefSTD.F90	(revision 3065)
@@ -43,5 +43,4 @@
   ! PARAMETER(klevSTD=17)
   INTEGER itap
-  ! REAL dtime
 
   ! variables locales
@@ -67,8 +66,8 @@
 
 
-    ! calcul variables tous les freq_calNMC(n)/dtime pas de temps
+    ! calcul variables tous les freq_calNMC(n)/phys_tstep pas de temps
     ! de la physique
 
-    IF (mod(itap,nint(freq_calnmc(n)/dtime))==0) THEN
+    IF (mod(itap,nint(freq_calnmc(n)/phys_tstep))==0) THEN
       DO k = 1, nlevstd
         DO i = 1, klon
@@ -103,5 +102,5 @@
       END DO !k
 
-    END IF !MOD(itap,NINT(freq_calNMC(n)/dtime)).EQ.0
+    END IF !MOD(itap,NINT(freq_calNMC(n)/phys_tstep)).EQ.0
 
   END DO !n
