Index: trunk/WRF.COMMON/INTERFACES_V4/dynphy_wrf_generic_lmd/callphysiq_mod.F
===================================================================
--- trunk/WRF.COMMON/INTERFACES_V4/dynphy_wrf_generic_lmd/callphysiq_mod.F	(revision 2866)
+++ trunk/WRF.COMMON/INTERFACES_V4/dynphy_wrf_generic_lmd/callphysiq_mod.F	(revision 2868)
@@ -43,10 +43,8 @@
 
   call allocate_comm_wrf(klon,llm)
-
 ! Call physics package with required inputs/outputs
   CALL physiq(klon,           & ! ngrid
               llm,            & ! nlayer
               nqtot,          & ! nq
-!              noms,          & ! nametrac
               debut_split,    & ! firstcall
               lafin_split,    & ! lastcall
Index: trunk/WRF.COMMON/INTERFACES_V4/dynphy_wrf_generic_lmd/update_inputs_physiq_mod.F
===================================================================
--- trunk/WRF.COMMON/INTERFACES_V4/dynphy_wrf_generic_lmd/update_inputs_physiq_mod.F	(revision 2866)
+++ trunk/WRF.COMMON/INTERFACES_V4/dynphy_wrf_generic_lmd/update_inputs_physiq_mod.F	(revision 2868)
@@ -27,6 +27,6 @@
   REAL :: sec,nsec
   
-  IF (JULYR .le. 8999) THEN
-    if (tlocked .eqv. .false.) THEN
+  !IF (JULYR .le. 8999) THEN
+  !  if (tlocked .eqv. .false.) THEN
       JH_cur_split = (GMT + elaps/3600.) !! universal time (0<JH_cur_split<1): JH_cur_split=0.5 at 12:00 UT 
       JH_cur_split = MODULO(JH_cur_split,24.)   !! the two arguments of MODULO must be of the same type 
@@ -36,32 +36,33 @@
       MY = (JULYR-2000) + (86400*(JULDAY - 1)+3600.0*GMT+elaps)/31968000
       MY = INT(MY)
-    ELSE
-      JH_cur_split = (GMT)! + elaps/420000.) !! universal time (0<JH_cur_split<1): JH_cur_split=0.5 at 12:00 UT 
-      JH_cur_split = MODULO(JH_cur_split,24.)   !! the two arguments of MODULO must be of the same type 
-      JH_cur_split = JH_cur_split / 24.
-      JD_cur = (JULDAY - 1 + INT((3600*GMT)/86400))!+elaps)/1.008e7))
-      JD_cur = MODULO(int(JD_cur),2)
-      MY = (JULYR-2000) + (86400*(JULDAY - 1)+3600*GMT+elaps)/31968000
-      MY = INT(MY)
-    ENDIF
-  ELSE
-    if (tlocked .eqv. .false.) THEN
-      JH_cur_split = lct_input - lon_input / 15. + elaps/1500.0
-      JD_cur =  INT((sec*(lct_input - lon_input / 15.) + elaps)/36000)
-      !ptime = lct_input - lon_input / 15. + elaps/3600.
-      !pday =  INT((3600*(lct_input - lon_input / 15.) + elaps)/86400)
-    ELSE
-      JH_cur_split = lct_input - lon_input / 15. !+ elaps/1500.0
-      !pday =  INT((sec*(lct_input - lon_input / 15.)+ elaps)/36000)
-      JD_cur =  INT((sec*(lct_input - lon_input / 15.))/3600)
-      !print*,'ptime',ptime
-      !print*,'pday',pday 
-    !pday =  INT((3600*(lct_input - lon_input / 15.) + elaps)/86400)
-      JH_cur_split = MODULO(ptime,24.)
-      JH_cur_split = JH_cur_split / 24.
-      JD_cur = MODULO(int(pday),365)
-      MY = 2024
-    ENDIF
-  ENDIF
+  !  ELSE
+  !    JH_cur_split = (GMT)! + elaps/420000.) !! universal time (0<JH_cur_split<1): JH_cur_split=0.5 at 12:00 UT 
+  !    JH_cur_split = MODULO(JH_cur_split,24.)   !! the two arguments of MODULO must be of the same type 
+  !    JH_cur_split = JH_cur_split / 24.
+  !    JD_cur = (JULDAY - 1 + INT((3600*GMT)/86400))!+elaps)/1.008e7))
+  !    JD_cur = MODULO(int(JD_cur),2)
+  !    MY = (JULYR-2000) + (86400*(JULDAY - 1)+3600*GMT+elaps)/31968000
+  !    MY = INT(MY)
+  !  ENDIF
+  !ELSE
+  !  if (tlocked .eqv. .false.) THEN
+  !    JH_cur_split = lct_input - lon_input / 15. + elaps/1500.0
+  !    JD_cur =  INT((sec*(lct_input - lon_input / 15.) + elaps)/36000)
+  !    !ptime = lct_input - lon_input / 15. + elaps/3600.
+  !    !pday =  INT((3600*(lct_input - lon_input / 15.) + elaps)/86400)
+  !  ELSE
+  !    JH_cur_split = lct_input - lon_input / 15. !+ elaps/1500.0
+  !    !pday =  INT((sec*(lct_input - lon_input / 15.)+ elaps)/36000)
+  !    JD_cur =  INT((sec*(lct_input - lon_input / 15.))/3600)
+  !    !print*,'ptime',ptime
+  !    !print*,'pday',pday 
+  !  !pday =  INT((3600*(lct_input - lon_input / 15.) + elaps)/86400)
+  !    JH_cur_split = MODULO(ptime,24.)
+  !    JH_cur_split = JH_cur_split / 24.
+  !    JD_cur = MODULO(int(pday),365)
+  !    MY = 2024
+  !  ENDIF
+  !ENDIF
+
 
 END SUBROUTINE update_inputs_physiq_time
@@ -80,5 +81,4 @@
       
                  !! tableau dans tracer_mod.F90
-  nqtot=nq
   if ((TRACER_MODE .eq. 1) .or. (TRACER_MODE .ge. 42)) THEN
     nqtot=2
@@ -107,9 +107,4 @@
 SUBROUTINE update_inputs_physiq_constants
 
-   !USE module_model_constants
-   !use comcstfi_h, only: omeg,mugaz
-   !use planete_h,  only: year_day,periheli,aphelie, &
-   !                       peri_day,obliquit,emin_turb, &
-   !                       lmixmin
    use planete_mod, only: year_day, periastr, apoastr, peri_day,&
                        obliquit, z0, lmixmin, emin_turb
@@ -117,67 +112,8 @@
                          emisice,dtemisice
    !                       z0_default
-   !use comsoil_h, only: volcapa
    use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, cpp, r
-   !! comcstfi_h
    use phys_state_var_mod, only :cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tsea_ice
-   !use iniorbit
-   !use time_phylmdz_mod, only: dtphys, daysec,day_ini
-  
-  
-     !open(17,file='controle.txt',form='formatted',status='old')
-     !rewind(17)
-     !read(17,*)
-     !read(17,*) 
-     !read(17,*) day_ini !(tab0+3)
-     !read(17,*) 
-     !read(17,*) !tab0+5)
-     !read(17,*) omeg !(tab0+6)
-     !read(17,*) !(tab0+7)
-     !read(17,*)  !(tab0+8)
-     !read(17,*)  !(tab0+9)
-     !read(17,*) daysec
-     !read(17,*) dtphys !tab0+11)
-     !read(17,*)
-     !read(17,*)
-     !read(17,*) year_day !(tab0+14)
-     !read(17,*) periastr !tab0+15)
-     !read(17,*) apoastr !tab0+16)
-     !read(17,*) peri_day !tab0+17)
-     !read(17,*) obliquit !tab0+18)
-     !read(17,*) z0
-     !read(17,*)
-     !read(17,*)
-     !read(17,*)
-     !read(17,*)
-     !read(17,*) emisice(1)
-     !read(17,*) emisice(2)
-     !read(17,*) emissiv
-     !read(17,*)
-     !read(17,*)
-     !read(17,*)
-     !read(17,*)
-     !read(17,*) iceradius(1)
-     !read(17,*) iceradius(2)
-     !read(17,*) dtemisice(1)
-     !read(17,*) dtemisice(2)
-     !close(17) 
-     !cpp=(8.314511/(mugaz/1000.0))/rcp
-     !print*,'cpp',cpp
-     !print*,'g',g
-
-     !emissiv(:)=EMIS
-     !cloudfrac(:,:)=0.5
-     !totcloudfrac(:)=0.5
-     !hice(:)=0.
-     !rnat(:)=0.
-     !pctsrf_sic(:)=0.
-     !tsea_ice(:)=0.
-     !qsurf(:,:) = 0. 
-     !print*,'iceradius',iceradius,'dtemisice',dtemisice
-     !print*,'apoastr,periastr,year_day,peri_day,obliq',apoastr,periastr,year_day,peri_day,obliquit
-     !print*,'emissiv',emissiv
- 
-     ! SUBROUTINE iniorbit(apoastr,periastr,year_day,peri_day,obliq)
-
+   !JL22 this routine does not do anything for the generic interface
+   ! The various use abave can surely be removed.
 END SUBROUTINE update_inputs_physiq_constants
 
@@ -272,9 +208,5 @@
    latitude_deg(:) = plat(:)/DEGRAD
    cell_area(:) = parea(:)
-   !call planetwide_sumval(parea,totarea_planet)
-   !print*,'parea',parea(1)
-   !totarea=SSUM(ngrid,parea,1)
    totarea=ngrid*parea(1)
-   !totarea_planet=SSUM(ngrid,parea,1)
    totarea_planet=ngrid*parea(1)
 
@@ -294,5 +226,4 @@
    DO k=1, nlayer
    read(12,*) znw(k)
-   !write(6,*) 'read level ', k,grid%znw(k)
    ENDDO
    close(12)
@@ -319,5 +250,5 @@
             JULYR,TRACER_MODE,&
             M_ALBEDO,CST_AL,&
-            M_TSURF,M_EMISS,M_CO2ICE,&
+            P_TSURF,M_EMISS,M_CO2ICE,&
             M_GW,M_Z0,CST_Z0,&
             M_H2OICE,&
@@ -335,9 +266,7 @@
    REAL, INTENT(IN  ) :: CST_AL, phisfi_val, CST_Z0
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)  :: &
-     M_ALBEDO,M_TSURF,M_EMISS,M_CO2ICE,M_H2OICE,M_Z0
+     M_ALBEDO,P_TSURF,M_EMISS,M_CO2ICE,M_H2OICE,M_Z0
    REAL, DIMENSION( ims:ime, 5, jms:jme ), INTENT(IN   )  :: M_GW  
 
-   !print*,'ALLOCATED(phisfi)',ALLOCATED(phisfi)
-   !print*,'size phisfi',size(phisfi)
    DO j = jps,jpe
    DO i = ips,ipe
@@ -352,6 +281,4 @@
      !---------------------!
      phisfi(subs) = phisfi_val
-!     print*,'size phisfi',size(phisfi)
-     !print*,'phisfi',phisfi(subs)
      !---------------!
      ! Ground albedo !
@@ -383,32 +310,9 @@
      !----------------------------!
      z0 = CST_Z0
-     !IF (JULYR .le. 8999) THEN
-     ! IF (CST_Z0 == 0) THEN
-     !  z0(subs) = M_Z0(i,j)
-     ! ELSE
-     !  z0(subs) = CST_Z0
-     !  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** SET CONSTANT SURF ROUGHNESS (m) ',CST_Z0
-     ! ENDIF
-     !ELSE
-     ! IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION z0 (m) ', CST_Z0
-     ! z0(subs)=CST_Z0
-     !ENDIF
-     !!!!! ADDITIONAL SECURITY. THIS MIGHT HAPPEN WITH OLD INIT FILES.
-     !IF (z0(subs) == 0.) THEN
-     !  IF ( (i == ips) .AND. (j == jps) ) PRINT *, 'WELL, z0 is 0, this is no good. Setting to old defaults value 0.01 m'
-     !  z0(subs) = 0.01
-     !ENDIF
-     !!!!! ADDITIONAL SECURITY. INTERP+SMOOTH IN GEOGRID MIGHT YIELD NEGATIVE Z0 !!!
-     !IF (z0(subs) < 0.) THEN
-     !  PRINT *, 'WELL, z0 is NEGATIVE, this is impossible. better stop here.'
-     !  PRINT *, 'advice --> correct interpolation / smoothing of z0 in WPS'
-     !  PRINT *, '           -- or check the constant value set in namelist.input'
-     !  STOP
-     !ENDIF
-
+     
      !-----------------------------------------------!
      ! Ground temperature, emissivity, CO2 ice cover !
      !-----------------------------------------------!
-     tsurf(subs) = M_TSURF(i,j)
+     tsurf(subs) = P_TSURF(i,j)
      emis(subs) = M_EMISS(i,j)     
      !do i=1,noceanmx
@@ -466,5 +370,5 @@
             M_TI,CST_TI,&
             M_ISOIL,M_DSOIL,&
-            M_TSOIL,M_TSURF)
+            M_TSOIL,P_TSURF)
 
    use comsoil_h, only: inertiedat,mlayer,layer,volcapa
@@ -476,5 +380,5 @@
    REAL, INTENT(IN  ) :: CST_TI 
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)  :: &
-     M_TI, M_TSURF
+     M_TI, P_TSURF
    REAL, DIMENSION( ims:ime, nsoil, jms:jme ), INTENT(IN)  :: &
      M_TSOIL, M_ISOIL, M_DSOIL
@@ -541,9 +445,5 @@
        IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** no tsoil. set it to tsurf.'
        do k=1,nsoil
-        !print*,'M_TSURF(i,j)',M_TSURF(1,:)
-        !print*,'size M_TSURF',size(M_TSURF)
-        !print*,'size tsoil',size(tsoil)
-        tsoil(subs,k) = M_TSURF(i,j)
-        !print*,'tsoil(subs,k)',tsoil(subs,k)
+        tsoil(subs,k) = P_TSURF(i,j)
        enddo
      ENDIF
Index: trunk/WRF.COMMON/INTERFACES_V4/dynphy_wrf_generic_lmd/update_outputs_physiq_mod.F
===================================================================
--- trunk/WRF.COMMON/INTERFACES_V4/dynphy_wrf_generic_lmd/update_outputs_physiq_mod.F	(revision 2866)
+++ trunk/WRF.COMMON/INTERFACES_V4/dynphy_wrf_generic_lmd/update_outputs_physiq_mod.F	(revision 2868)
@@ -9,5 +9,5 @@
             ips,ipe,jps,jpe,&
             TRACER_MODE,&
-            M_TSURF,M_CO2ICE,&
+            P_TSURF,M_CO2ICE,&
             M_H2OICE)
 
@@ -20,18 +20,18 @@
    INTEGER :: i,j,subs
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)  :: &
-     M_TSURF,M_CO2ICE,M_H2OICE
-
-   DO j = jps,jpe
-   DO i = ips,ipe
-
-     !-----------------------------------!
-     ! 1D subscript for physics "cursor" !
-     !-----------------------------------!
-     subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
-
-     !-------------------------------------------------------!
-     ! Save key variables for restart and output and nesting !  
-     !-------------------------------------------------------!
-     M_TSURF(i,j) = tsurf(subs)
+     P_TSURF,M_CO2ICE,M_H2OICE
+
+   DO j = jps,jpe
+   DO i = ips,ipe
+
+     !-----------------------------------!
+     ! 1D subscript for physics "cursor" !
+     !-----------------------------------!
+     subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
+
+     !-------------------------------------------------------!
+     ! Save key variables for restart and output and nesting !  
+     !-------------------------------------------------------!
+     P_TSURF(i,j) = tsurf(subs)
 
    ENDDO
Index: trunk/WRF.COMMON/INTERFACES_V4/module_lmd_driver.F
===================================================================
--- trunk/WRF.COMMON/INTERFACES_V4/module_lmd_driver.F	(revision 2866)
+++ trunk/WRF.COMMON/INTERFACES_V4/module_lmd_driver.F	(revision 2868)
@@ -26,4 +26,5 @@
         P8W,DZ8W,T8W,Z,HT,MUT,DNW, &
         U,V,TH,T,P,EXNER,RHO, &
+        P_HYD, P_HYD_W, &
         PTOP, &
         RADT, &
@@ -35,5 +36,5 @@
         planet_type, &
         M_ALBEDO,M_TI,M_CO2ICE,M_EMISS, &
-        M_H2OICE,M_TSOIL,M_Q2,M_TSURF, &
+        M_H2OICE,M_TSOIL,M_Q2,P_TSURF, &
         M_FLUXRAD,M_WSTAR,M_ISOIL,M_DSOIL,&
         M_Z0, CST_Z0, M_GW, &
@@ -41,6 +42,5 @@
         CST_AL, CST_TI, &
         isfflx, diff_opt, km_opt, &
-        HISTORY_INTERVAL, &
-        HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,DT_AJS,&
+        HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,&
         CLOUDFRAC,TOTCLOUDFRAC,RH, &
         DQICE,DQVAP,REEVAP,SURFRAIN,ALBEQ,FLUXTOP_DN,FLUXABS_SW,FLUXTOP_LW,FLUXSURF_SW,&
@@ -94,5 +94,4 @@
 REAL, INTENT(IN  ) :: CST_AL, CST_TI 
 REAL, INTENT(IN  ) :: PTOP
-INTEGER, INTENT(IN   ) :: HISTORY_INTERVAL
 ! 2D arrays          
 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)  :: &
@@ -103,5 +102,5 @@
      SLPX,SLPY, &
      M_CO2ICE,M_H2OICE, &
-     M_TSURF, M_Z0, &
+     P_TSURF, M_Z0, &
      M_FLUXRAD,M_WSTAR, &
      PSFC,TSK
@@ -116,8 +115,8 @@
 ! 3D arrays 
 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
-     dz8w,p8w,p,exner,t,t8w,rho,u,v,z,th
+     dz8w,p8w,p,exner,t,t8w,rho,u,v,z,th,p_hyd,p_hyd_w
 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
      RTHBLTEN,RUBLTEN,RVBLTEN, &
-     HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,DT_AJS,RDUST,VMR_ICE,RICE,&
+     HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,RDUST,VMR_ICE,RICE,&
      CLOUDFRAC,RH,DQICE,DQVAP,DTLSC,DTRAIN,DT_MOIST,H2OICE_REFF
 REAL, DIMENSION( ims:ime, kms:kme+1, jms:jme ), INTENT(INOUT ) :: &
@@ -156,5 +155,4 @@
    ! <------ outputs:
    !     physical tendencies
-   REAL*8,DIMENSION(:,:),ALLOCATABLE :: pdtheta
    ! ... intermediate arrays
    REAL, DIMENSION(:), ALLOCATABLE  :: &
@@ -171,26 +169,7 @@
    LOGICAL, SAVE :: firstcall = .true.
    INTEGER, SAVE :: previous_id = 0
-!**************************************************
-! IMPORTANT: pour travailler avec grid nesting
-! --- deux comportements distincts du save
-! ... ex1: ferme planeto avec PGF+MPI: standard
-! ... ex2: iDataPlex avec IFORT+MPI: SPECIAL_NEST_SAVE 
-!**************************************************
-#ifdef SPECIAL_NEST_SAVE
-      !  une dimension supplementaire liee au nest
-      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
-             dp_save
-      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
-             du_save, dv_save, dt_save,dtheta_save
-      REAL, DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: &
-             dq_save      
-#else
-      REAL, DIMENSION(:), ALLOCATABLE, SAVE :: &
-             dp_save
-      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
-             du_save, dv_save, dt_save,dtheta_save
-      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
-             dq_save      
-#endif
+   REAL, DIMENSION(:), ALLOCATABLE, SAVE :: dp_save
+   REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: du_save, dv_save, dt_save
+   REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: dq_save      
 
 !!!IDEALIZED IDEALIZED
@@ -270,5 +249,4 @@
 ! DIMENSIONS FOR THE PHYSICS !
 !----------------------------!
-!day_ini = JULDAY - 1      !! GCM convention  !! pday at firstcall is day_ini
 wappel_phys = RADT
 zdt_split = dt*wappel_phys            ! physical timestep (s)
@@ -282,7 +260,4 @@
 lastcall = .false.
 ! **** needed but hardcoded
-
-          !PRINT *, ips, ipe, jps, jpe
-          !PRINT *, ngrid
 
 elaps = (float(itimestep)-1.)*dt  ! elapsed seconds of simulation
@@ -300,28 +275,4 @@
   !firstcall=.false. !! just in case you'd want to get rid of the physics 
 test=0
-#ifdef SPECIAL_NEST_SAVE
-PRINT *, 'ALLOCATED dp_save ???', ALLOCATED( dp_save ), id
-IF( .NOT. ALLOCATED( dp_save  ) ) THEN
-   PRINT *, '**** check **** OK I ALLOCATE one save ARRAY for all NESTS ', max_dom, id
-   !! here are the arrays that would be useful to save physics tendencies
-   ALLOCATE(dp_save(ngrid,max_dom))
-   ALLOCATE(du_save(ngrid,nlayer,max_dom))
-   ALLOCATE(dv_save(ngrid,nlayer,max_dom))
-   ALLOCATE(dt_save(ngrid,nlayer,max_dom))
-   ALLOCATE(dtheta_save(ngrid,nlayer,max_dom))
-   ALLOCATE(dq_save(ngrid,nlayer,nq,max_dom))
-   dp_save(:,:)=0.    !! initialize these arrays ...
-   du_save(:,:,:)=0.
-   dv_save(:,:,:)=0.
-   dt_save(:,:,:)=0.
-   dtheta_save(:,:,:)=0.
-   dq_save(:,:,:,:)=0.
-ENDIF
-IF (id .lt. max_dom) THEN
-   flag_first_restart=.true.
-ELSE
-   flag_first_restart=.false.
-ENDIF
-#else
 IF( .NOT. ALLOCATED( dp_save  ) ) THEN
 ALLOCATE(dp_save(ngrid)) !! here are the arrays that would be useful to save physics tendencies
@@ -329,5 +280,4 @@
 ALLOCATE(dv_save(ngrid,nlayer))
 ALLOCATE(dt_save(ngrid,nlayer))
-ALLOCATE(dtheta_save(ngrid,nlayer))
 ALLOCATE(dq_save(ngrid,nlayer,nq))
 ENDIF
@@ -336,9 +286,8 @@
 dv_save(:,:)=0.
 dt_save(:,:)=0.
-dtheta_save(:,:)=0.
 dq_save(:,:,:)=0.
 flag_first_restart=.false.
-#endif
-ELSE
+
+ELSE ! is_first_step
 !-------------------------------------------------!
 ! what is done for the other steps of simulation  !
@@ -396,5 +345,4 @@
 !----------!
 CALL allocate_interface(ngrid,nlayer,nq)
-ALLOCATE(pdtheta(ngrid,nlayer))
 !!!
 !!! BIG LOOP : 1. no call for physics, used saved values
@@ -402,12 +350,4 @@
 IF (test.NE.0) THEN
 print *,'** ',planet_type,'** NO CALL FOR PHYSICS, go to next step...',test
-#ifdef SPECIAL_NEST_SAVE
-zdpsrf_omp(:)=dp_save(:,id)
-zdufi_omp(:,:)=du_save(:,:,id)
-zdvfi_omp(:,:)=dv_save(:,:,id)
-zdtfi_omp(:,:)=dt_save(:,:,id)
-pdtheta(:,:)=dtheta_save(:,:,id)
-zdqfi_omp(:,:,:)=dq_save(:,:,:,id)
-#else
 print*,'else'
 zdpsrf_omp(:)=dp_save(:)
@@ -415,11 +355,9 @@
 zdvfi_omp(:,:)=dv_save(:,:)
 zdtfi_omp(:,:)=dt_save(:,:)
-pdtheta(:,:)=dtheta_save(:,:)
 zdqfi_omp(:,:,:)=dq_save(:,:,:)
-#endif
 !!!
 !!! BIG LOOP : 2. calculate physical tendencies
 !!!
-ELSE
+ELSE ! if (test.EQ.0)
 !----------!
 ! ALLOCATE !
@@ -489,5 +427,7 @@
 !--------------------------------------!
 dz8w_prof(:) = dz8w(i,kps:kpe,j)   ! dz between full levels (m)   
-!p_prof(:) = p(i,kps:kpe,j)         ! pressure half level (Pa) >> zplay_omp
+p_prof(:) = p_hyd(i,kps:kpe,j)         ! hydrostatic pressure at layers >> zplay_omp
+p8w_prof(:) = p_hyd_w(i,kps:kpe+1,j)         ! hydrostatic pressure at levels
+!JL22 using hydrostatic pressures to conserve mass
 t_prof(:) = t(i,kps:kpe,j)         ! temperature half level (K) >> pt
 t8w_prof(:) = t8w(i,kps:kpe,j)     ! temperature full level (K)
@@ -495,22 +435,4 @@
 v_prof(:) = v(i,kps:kpe,j)         ! meridional wind (A-grid: unstaggered) half level (m/s) >> pv
 z_prof(:) = z(i,kps:kpe,j)         ! geopotential height half level (m) >> zphi_omp/g
-
-
-!------------------------------------------------------------!
-! reconstructing hydrostatic level pressure to conserve mass !
-!------------------------------------------------------------!
-p8w_prof(kpe+1) = ptop
-if (TRACER_MODE .ge. 42) then
-  do k=kpe,kps,-1
-    p8w_prof(k) = p8w_prof(k+1) - dnw(k) * mut(i,j) * (1.d0 + moist(i,k,j,P_QV))
-    p_prof(k) = 0.5d0*(p8w_prof(k+1)+p8w_prof(k))         ! pressure half level (Pa) >> zplay_omp
-  enddo
-else
-  do k=kpe,kps,-1
-    p8w_prof(k) = p8w_prof(k+1) - dnw(k) * mut(i,j)
-    p_prof(k) = 0.5d0*(p8w_prof(k+1)+p8w_prof(k))         ! pressure half level (Pa) >> zplay_omp
-  enddo
-endif
-!p8w_prof(:) = p8w(i,kps:kpe,j)     ! pressure full level (Pa) >> zplev_omp
 
 
@@ -644,5 +566,5 @@
             M_TI,CST_TI,&
             M_ISOIL,M_DSOIL,&
-            M_TSOIL,M_TSURF)
+            M_TSOIL,P_TSURF)
 !!!
 CALL update_inputs_physiq_surf( &
@@ -651,5 +573,5 @@
             JULYR,TRACER_MODE,&
             M_ALBEDO,CST_AL,&
-            M_TSURF,M_EMISS,M_CO2ICE,&
+            P_TSURF,M_EMISS,M_CO2ICE,&
             M_GW,M_Z0,CST_Z0,&
             M_H2OICE,&
@@ -678,5 +600,4 @@
 !!!!!!!!!!!!!!!!!!!!!!
 
-call_physics : IF (wappel_phys .ne. 0.) THEN
 !!! initialize tendencies (security)
 zdpsrf_omp(:)=0.
@@ -684,6 +605,7 @@
 zdvfi_omp(:,:)=0.
 zdtfi_omp(:,:)=0.
-pdtheta(:,:)=0.
 zdqfi_omp(:,:,:)=0.
+
+call_physics : IF (wappel_phys .ne. 0.) THEN
 !print *, '** ',planet_type,'** CALL TO LMD PHYSICS'
 !!!
@@ -699,44 +621,12 @@
 !!!
 
-!! specific scenario. necessary to add the right amount of dust.
-#ifdef DUSTSTORM
-IF (firstcall .EQV. .true.) THEN
-  zdqfi_omp(:,:,:) = zdqfi_omp(:,:,:) / dt
-ENDIF
-#endif
-
-IF (planet_type .eq. "venus" ) THEN
-  DO j=jps,jpe
-  DO i=ips,ipe
-    do k=kps,kpe
-      subs=(j-jps)*(ipe-ips+1)+(i-ips+1)
-      tk1=(ztfi_omp(subs,k)**nu + nu*TT00**nu*log((p1000mb/zplay_omp(subs,k))**rcp))**(1/nu)
-      tk2=((ztfi_omp(subs,k) + zdtfi_omp(subs,k))**nu + nu*TT00**nu*log((p1000mb/zplay_omp(subs,k))**rcp))**(1/nu)
-      pdtheta(subs,k)=tk2-tk1
-    enddo
-  ENDDO
-  ENDDO
-ENDIF
-
-!print *, '** ',planet_type,'** CALL TO LMD PHYSICS DONE'
-
 !---------------------------------------------------------------------------------!
 ! PHYSIQ TENDENCIES ARE SAVED TO BE SPLIT WITHIN INTERMEDIATE DYNAMICAL TIMESTEPS !
 !---------------------------------------------------------------------------------!
-#ifdef SPECIAL_NEST_SAVE
-dp_save(:,id)=zdpsrf_omp(:)
-du_save(:,:,id)=zdufi_omp(:,:)
-dv_save(:,:,id)=zdvfi_omp(:,:)
-dt_save(:,:,id)=zdtfi_omp(:,:)
-dtheta_save(:,:,id)=pdtheta(:,:)
-dq_save(:,:,:,id)=zdqfi_omp(:,:,:)
-#else
 dp_save(:)=zdpsrf_omp(:)
 du_save(:,:)=zdufi_omp(:,:)
 dv_save(:,:)=zdvfi_omp(:,:)
 dt_save(:,:)=zdtfi_omp(:,:)
-dtheta_save(:,:)=pdtheta(:,:)
 dq_save(:,:,:)=zdqfi_omp(:,:,:)
-#endif
 
 !! OUTPUT OUTPUT OUTPUT
@@ -749,5 +639,5 @@
             ips,ipe,jps,jpe,&
             TRACER_MODE,&
-            M_TSURF,M_CO2ICE,&
+            P_TSURF,M_CO2ICE,&
             M_H2OICE)
 !!!
@@ -812,14 +702,10 @@
 
     ! zonal wind
-  RUBLTEN(i,kps:kpe,j) = zdufi_omp(subs,kps:kpe) 
+  RUBLTEN(i,kps:kpe,j) = zdufi_omp(subs,kps:kpe)
     ! meridional wind
   RVBLTEN(i,kps:kpe,j) = zdvfi_omp(subs,kps:kpe) 
     ! potential temperature
     ! (dT = dtheta * exner for isobaric coordinates or if pressure variations are negligible)
-  IF (planet_type .eq. "venus" ) THEN
-    RTHBLTEN(i,kps:kpe,j) = pdtheta(subs,kps:kpe)
-  ELSE
-    RTHBLTEN(i,kps:kpe,j) = zdtfi_omp(subs,kps:kpe) / exner(i,kps:kpe,j)
-  ENDIF
+  RTHBLTEN(i,kps:kpe,j) = zdtfi_omp(subs,kps:kpe) / exner(i,kps:kpe,j)
     ! update surface pressure (cf CO2 cycle in physics)
     ! here dt is needed
@@ -841,6 +727,6 @@
       scalar(i,kps:kpe,j,P_QH2O_ICE)=scalar(i,kps:kpe,j,P_QH2O_ICE) &
            +zdqfi_omp(subs,kps:kpe,2)*dt * (1.d0+moist(i,kps:kpe,j,P_QV))
-      ! if you want to use this mode, RTHBLTEN should be corrected as below. 
-      ! we keep it like that for the moment for testing.
+       ! if you want to use this mode, RTHBLTEN should be corrected as below. 
+       ! we keep it like that for the moment for testing.
     CASE(43)
       moist(i,kps:kpe,j,P_QV)=moist(i,kps:kpe,j,P_QV) &
@@ -864,5 +750,4 @@
 ENDDO
 CALL deallocate_interface
-DEALLOCATE(pdtheta)
 !!*****!!
 !! END !!