Index: /trunk/LMDZ.MARS/changelog.txt
===================================================================
--- /trunk/LMDZ.MARS/changelog.txt	(revision 3097)
+++ /trunk/LMDZ.MARS/changelog.txt	(revision 3098)
@@ -4273,4 +4273,20 @@
 Prettyfied solarlong.F and made it a module. Likewise for conf_phys.F
 
+== 23/10/2023 == EV
+We added the to module vdifc the possibily of subsurface intercation, mostly to have the option of buried glaicer that can lose ice and create polar layers.
+if no ice is on the surface the ice the atmosphere can directly interact with the subsurface ice.
+if ice is on the surface it can interact with the subsurface ice.
+**For now it is not working with more then one subslope. I will add this feature soon**
+In order to use it there are two flags that needs to be set to true.
+paleoclimate=.true.
+lag_layer=.true.
+in addition one need to give the following parametrs in startfi or run.def:
+d_coef=X for the diffusion coeficent (4e-4 is the defualt)
+<h2o_ice_depth=X> if its zero or lower there is no subsurface ice (default -1).
+outputs: zdq_ssi, momentarily zdqsdif after the interaction with the SSIzdq_subtimestep,  momentarily flux after the interaction with the SSIzdq_ssi_frost momentarily zdqsdif after the interaction of the frost with the SSI
+q1, atmospheric first layer water vapor quatity
+qeq, SSI water vapor quantity
+
+
 Some code tidying:
 Made pi in module comcstfi_h a parameter (and not a variable recomputed at
Index: /trunk/LMDZ.MARS/libf/phymars/conf_phys.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 3097)
+++ /trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 3098)
@@ -37,5 +37,6 @@
       use aeropacity_mod, only: iddist, topdustref
       USE mod_phys_lmdz_transfert_para, ONLY: bcast
-      USE paleoclimate_mod,ONLY: paleoclimate,albedo_perenialco2
+      USE paleoclimate_mod,ONLY: paleoclimate,albedo_perenialco2,
+     &                           lag_layer 
       use microphys_h, only: mteta
 
@@ -338,8 +339,14 @@
 ! PALEOCLIMATE
 
-         write(*,*)"Is it paleoclimate run?"
+         write(*,*)"Using lag layer??"
+         lag_layer=.false.
+         call getin_p("lag_layer",lag_layer)
+         write(*,*) " lag_layer = ", lag_layer
+
+        write(*,*)"Is it paleoclimate run?"
          paleoclimate=.false. ! default value
          call getin_p("paleoclimate",paleoclimate)
          write(*,*)" paleoclimate = ",paleoclimate
+
 
          write(*,*)"Albedo for perenial CO2 ice?"
Index: /trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90	(revision 3097)
+++ /trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90	(revision 3098)
@@ -477,4 +477,5 @@
     call getin("subsurface_ice_depth",ice_depth)
 
+
     z0(1) = z0_default ! default value for roughness
     write(*,*) 'Surface roughness length z0 (m)?'
@@ -716,4 +717,8 @@
         write(*,*) 'Prescribed atmospheric water vapor profile'
         write(*,*) 'Unless it reaches saturation (maximal value)'
+    else if (atm_wat_profile .eq. 2) then
+        write(*,*) 'Prescribed atmospheric water vapor profile'
+        write(*,*) 'like present day northern midlatitude'
+        write(*,*) 'Unless it reaches saturation (maximal value)'
     else
         error stop 'Water vapor profile value not correct!'
Index: /trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90	(revision 3097)
+++ /trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90	(revision 3098)
@@ -179,5 +179,17 @@
             endif
         endif
-    endif
+    !EV profile
+!            IF(atm_wat_profile.eq.2) THEN
+!             DO ilayer=1,16
+!              q(1,ilayer,igcm_h2o_vap)=min(zqsat(ilayer),6e-5)
+!             ENDDO! ilayer=1,16
+!             DO ilayer=17,22
+!              q(1,ilayer,igcm_h2o_vap)=min(zqsat(ilayer),3e-5)
+!             ENDDO! ilayer=17,22
+!             DO ilayer=23,nlayer
+!              q(1,ilayer,igcm_h2o_vap)=0
+!             ENDDO! ilayer=23,nlayer
+!             endif
+     endif
 
     ! Call physics
Index: /trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90	(revision 3097)
+++ /trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90	(revision 3098)
@@ -15,8 +15,10 @@
 !$OMP THREADPRIVATE(paleoclimate)
 
-    real, save, allocatable :: lag_h2o_ice(:,:)  ! Thickness of the lag before H2O ice [m]
+    real, save, allocatable :: h2o_ice_depth(:,:)  ! Thickness of the lag before H2O ice [m]
     real, save, allocatable :: lag_co2_ice(:,:)  ! Thickness of the lag before CO2 ice [m]
     real, save :: albedo_perenialco2             ! Albedo for perenial co2 ice [1]
-!$OMP THREADPRIVATE(lag_h2o_ice,lag_co2_ice,albedo_perenialco2)
+    real, save, allocatable :: d_coef(:,:)  ! Diffusion coeficent
+    LOGICAL,SAVE :: lag_layer ! does lag layer is present?
+!$OMP THREADPRIVATE(h2o_ice_depth,lag_co2_ice,albedo_perenialco2)
 
     CONTAINS
@@ -29,6 +31,7 @@
   integer,intent(in) :: nslope ! number of slope within a mesh 
 
-    allocate(lag_h2o_ice(ngrid,nslope))
+    allocate(h2o_ice_depth(ngrid,nslope))
     allocate(lag_co2_ice(ngrid,nslope))
+    allocate(d_coef(ngrid,nslope))
   end subroutine ini_paleoclimate_h
 
@@ -36,5 +39,6 @@
 
   implicit none
-    if (allocated(lag_h2o_ice)) deallocate(lag_h2o_ice)
+    if (allocated(d_coef)) deallocate(d_coef)
+    if (allocated(h2o_ice_depth)) deallocate(h2o_ice_depth)
     if (allocated(lag_co2_ice)) deallocate(lag_co2_ice)
   end subroutine end_paleoclimate_h
Index: /trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 3097)
+++ /trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 3098)
@@ -27,5 +27,5 @@
   use comsoil_h, only: flux_geo
   USE comslope_mod, ONLY: nslope, major_slope
-  USE paleoclimate_mod, ONLY: paleoclimate, lag_h2o_ice, lag_co2_ice
+  USE paleoclimate_mod, ONLY: paleoclimate, h2o_ice_depth, lag_co2_ice, d_coef
   USE comcstfi_h, only: pi
   use geometry_mod, only: latitude
@@ -744,10 +744,19 @@
   if (startphy_file) then
 ! Depth of H2O lag
-   call get_field("lag_h2o_ice",lag_h2o_ice,found,indextime)
-   if (.not.found) then
-     write(*,*) "phyetat0: Failed loading <lag_h2o_ice> : ", &
-                          "<lag_h2o_ice> is set as -1 (no subsurface ice)"
-     lag_h2o_ice(:,:) = -1.
-   endif
+   call get_field("h2o_ice_depth",h2o_ice_depth,found,indextime)
+   if (.not.found) then
+     write(*,*) "phyetat0: Failed loading <h2o_ice_depth> : ", &
+                          "<h2o_ice_depth> is set as -1 (no subsurface ice)"
+     h2o_ice_depth(:,:) = -1.
+   endif
+
+! Diffuesion coeficent 
+   call get_field("d_coef",d_coef,found,indextime)
+   if (.not.found) then
+     write(*,*) "phyetat0: Failed loading <d_coef> : ", &
+                          "<d_coef> is set as 4e-4 (defualt value)"
+     d_coef(:,:) = 4e-4
+   endif
+
 
 !  Depth of CO2 lag
@@ -774,7 +783,8 @@
    endif ! notfound
   else ! no startfiphyle
-     lag_h2o_ice(:,:) = -1.
+     h2o_ice_depth(:,:) = -1.
      lag_co2_ice(:,:) = -1.
      perenial_co2ice(:,:) = 0.
+     d_coef(:,:)= 4.e-4
      if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then
        DO islope = 1,nslope
@@ -786,7 +796,8 @@
   endif !startphy_file
 else
-  lag_h2o_ice(:,:) = -1.
+  h2o_ice_depth(:,:) = -1.
   lag_co2_ice(:,:) = -1.
   perenial_co2ice(:,:) = 0.
+  d_coef(:,:)= 4.e-4
 endif !paleoclimate
 
Index: /trunk/LMDZ.MARS/libf/phymars/phyredem.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 3097)
+++ /trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 3098)
@@ -25,5 +25,5 @@
   use time_phylmdz_mod, only: daysec
   use comslope_mod, ONLY: nslope
-  USE paleoclimate_mod, ONLY: paleoclimate, lag_h2o_ice, lag_co2_ice  
+  USE paleoclimate_mod, ONLY: paleoclimate, h2o_ice_depth, lag_co2_ice  
   implicit none
  
@@ -156,5 +156,5 @@
   ! Paleoclimate outputs
   if(paleoclimate) then
-    call put_field("lag_h2o_ice","Depth of the H2O lags",lag_h2o_ice)
+    call put_field("h2o_ice_depth","Depth to the fisrt H2O ice",h2o_ice_depth)
     call put_field("lag_co2_ice","Depth of the CO2 lags",lag_co2_ice)
   endif
Index: /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 3097)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 3098)
@@ -1514,5 +1514,6 @@
      &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,
      &        zcondicea_co2microp,sensibFlux,
-     &        dustliftday,local_time,watercap,dwatercap_dif)
+     &        dustliftday,local_time,watercap,dwatercap_dif,
+     &        tsoil,nsoilmx)
 
           DO ig=1,ngrid
Index: /trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 3097)
+++ /trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 3098)
@@ -13,6 +13,6 @@
      $                pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,
      $                hfmax,pcondicea_co2microp,sensibFlux,
-     $                dustliftday,local_time,watercap, dwatercap_dif)
-
+     $                dustliftday,local_time,watercap, dwatercap_dif,
+     $                tsoil,nsoil)
       use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
      &                      igcm_dust_submicron, igcm_h2o_vap,
@@ -32,4 +32,6 @@
      &                      subslope_dist,major_slope,iflat
       use microphys_h, only: To
+      use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer
+      use comsoil_h, only: layer, mlayer
 
       IMPLICIT NONE
@@ -61,5 +63,5 @@
 c   ----------
 
-      INTEGER,INTENT(IN) :: ngrid,nlay
+      INTEGER,INTENT(IN) :: ngrid,nlay,nsoil
       REAL,INTENT(IN) :: ptimestep
       REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
@@ -75,4 +77,5 @@
       REAL,INTENT(IN) :: pcapcal(ngrid,nslope)
       REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
+      REAL,INTENT(IN) :: tsoil(ngrid,nsoil,nslope)
 
 c    Argument added for condensation:
@@ -102,5 +105,5 @@
       REAL :: pt(ngrid,nlay)
  
-      INTEGER ilev,ig,ilay,nlev,islope
+      INTEGER ilev,ig,ilay,nlev,islope,ik,lice
 
       REAL z4st,zdplanck(ngrid)
@@ -118,4 +121,11 @@
       REAL zcst1
       REAL zu2(ngrid)
+      REAL Tice(ngrid,nslope) ! subsurface temperature where ice is located.
+      REAL qeq(ngrid,nslope) ! saturation water vapor in the subsurface
+      REAL dist_up(ngrid,nslope) !distance from ice to layer above
+      REAL dist_down(ngrid,nslope) !distance from ice to layer down
+      REAL dist_sum(ngrid,nslope) ! sum of distance
+      REAL zdqsdif_ssi(ngrid,nslope) !SSI flux
+      REAL zdqsdif_ssi_frost(ngrid,nslope) !SSI-frost interaction
 
       EXTERNAL SSUM,SCOPY
@@ -155,4 +165,6 @@
       REAL rho(ngrid) ! near surface air density
       REAL qsat(ngrid)
+      REAL qsat2(ngrid,nslope)
+      REAL resist(ngrid,nslope) !subsurface ice flux reduction coef
 
       REAL hdoflux(ngrid,nslope)       ! value of vapour flux of HDO
@@ -272,5 +284,5 @@
       zdqsdif(1:ngrid)=0
       dwatercap_dif(1:ngrid,1:nslope)=0
-
+!      h2o_ice_depth(1:ngrid,1:nslope)=5
 c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
 c       avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
@@ -914,5 +926,5 @@
 c      Utilisation d'un sous pas de temps afin
 c      de decrire le flux de chaleur latente 
-
+         
         do iq=1,nq  
           if ((water).and.(iq.eq.igcm_h2o_vap)) then
@@ -920,4 +932,7 @@
           DO islope = 1,nslope
            DO ig=1,ngrid
+
+             
+
              zqsurf(ig)=pqsurf(ig,igcm_h2o_ice,islope)/
      &             cos(pi*def_slope_mean(islope)/180.)
@@ -943,5 +958,4 @@
             zq_tmp_vap(ig,:,:) =zq(ig,:,:)
 c           Debut du sous pas de temps
-
             DO tsub=1,nsubtimestep(ig)
 
@@ -980,7 +994,89 @@
                if ((-zdqsdif(ig)*subtimestep)
      &            .gt.(zqsurf(ig))) then
+                
+    
+                 !EV subsurface ice
+                IF(h2o_ice_depth(ig,islope) .gt. 0 .and. lag_layer)
+     &           then 
+                zdqsdif(ig)=
+     &                        -zqsurf(ig)/subtimestep
+                zqsurf(ig)=0
+
+               DO ik=0,nsoil-2 ! go through all the layers to find the ice locations
+                IF((mlayer(ik).le.h2o_ice_depth(ig,islope)).and.
+     &           (mlayer(ik+1).gt.h2o_ice_depth(ig,islope))) THEN
+                 lice = ik+1
+                EXIT
+                ENDIF
+               ENDDO !of subsurface loop
+                IF (lice .gt. 1) then !calculate the distance from the layers
+                 dist_up(ig,islope)=(h2o_ice_depth(ig,islope)
+     &                              -mlayer(lice-1))
+                 dist_down(ig,islope)=(mlayer(lice)
+     &                                -h2o_ice_depth(ig,islope))
+                 dist_sum(ig,islope)=dist_up(ig,islope)
+     &                               +dist_down(ig,islope)
+                 Tice(ig,islope)=(dist_up(ig,islope) ! Linear interp to calculate the temp
+     &                          *tsoil(ig,lice-1,islope)
+     &           /dist_sum(ig,islope))+
+     &           (dist_down(ig,islope)*tsoil(ig,lice,islope)
+     &           /dist_sum(ig,islope))
+                        ELSE
+                                Tice(ig,islope)=tsoil(ig,1,islope)
+                        ENDIF
+              call watersat(1,Tice(ig,1),pplev(ig,1)
+     &                     ,qsat2(ig,1))
+              qeq(ig,1)=(ztsrf(ig)/Tice(ig,1))
+     &                        *qsat2(ig,1)
+              resist(ig,1)=(1+(h2o_ice_depth(ig,islope)*zcdh(ig)
+     &                           /d_coef(ig,islope)))
+              !write(*,*)'R=',resist(ig,islope)
+              !write(*,*)'zice=',h2o_ice_depth(ig,islope)
+              !write(*,*)'D=',d_coef(ig,islope)
+              !write(*,*)'zcdh=',zcdh(ig)
+              zb(ig,1)=zb(ig,1)/resist(ig,1) ! change zb to account subsurface ice
+              !vdifc algorithem !!!!!needs to change  reseist io to the mean
+              !beacuse of the slopes!!!!!
+              z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+             zc(ig,nlay)=za(ig,nlay)*zq_tmp_vap(ig,nlay,iq)*z1(ig)
+             zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+             DO ilay=nlay-1,2,-1
+                 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
+     $            zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+                 zc(ig,ilay)=(za(ig,ilay)*zq_tmp_vap(ig,ilay,iq)+
+     $            zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+                 zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+             ENDDO
+             z1(ig)=1./(za(ig,1)+zb(ig,1)+
+     $          zb(ig,2)*(1.-zd(ig,2)))
+             zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+
+     $        zb(ig,2)*zc(ig,2)) * z1(ig)
+             zd(ig,1)=zb(ig,1)*z1(ig)
+             zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
+             zdqsdif_ssi(ig,1)=rho(ig)*dryness(ig)*zcdv(ig)    
+     &                      *(zq1temp(ig)-qeq(ig,islope))
+                      call write_output('zdq_ssi',
+     &                '','',zdqsdif_ssi(ig,1))
+                       call write_output('qeq',
+     &                '','',qeq(ig,1))
+                       call write_output('q1',
+     &                '','',zq1temp(ig))
+             !write(*,*)'qeq=',qeq(ig,1)
+             !write(*,*)'q1=',zq1temp(ig) 
+             !write(*,*)'zdqsdif_ssi=',zdqsdif_ssi(ig,1)
+             !I should some all the interactions with the SSI accoridng
+             !to the statistics and the evaluate
+             
+             !write(*,*) "zdq_ssi*t", zdqsdif_ssi(ig)*subtimestep
+             if (zdqsdif_ssi(ig,1)<0) then
+               !zdqsdif(ig)=zdqsdif_ssi(ig,1)-(zqsurf(ig)/subtimestep)
+               zdqsdif(ig)=zdqsdif(ig)+zdqsdif_ssi(ig,1)
+                call write_output('zdq_zdqssi',
+     &                '','',zdqsdif(ig)+zdqsdif_ssi(ig,1))
+             endif
+              ELSE
 c             pdqsdif > 0 : ice condensing
 c             pdqsdif < 0 : ice subliming
-c             write(*,*) "subliming more than available frost:  qsurf!"
+             write(*,*) "subliming more than available frost:  qsurf!"
                   zdqsdif(ig)=
      &                        -zqsurf(ig)/subtimestep
@@ -991,11 +1087,37 @@
      $           (-zdqsdif(ig)) *subtimestep) *z1(ig)
                  zq1temp(ig)=zc(ig,1)
+                 ENDIF !if h2o_ice_depth>0 and lag_layer
                endif  !if .not.watercaptag(ig) 
              endif ! if sublim more than surface
+
+              ! subsurface ice <--> frost interaction !EV
+             IF(h2o_ice_depth(ig,islope) .gt. 4e-4 .and. lag_layer
+     &        .and. zqsurf(ig) .gt. 0) then
+               zdqsdif_ssi_frost(ig,1)=(d_coef(ig,1)
+     &                         /h2o_ice_depth(ig,1))
+     &          *rho(ig)*dryness(ig)*(qsat(ig)-qeq(ig,1))
+               !needs to change to the mean of eq
+               zdqsdif(ig)=zdqsdif(ig)-zdqsdif_ssi_frost(ig,1) 
+                !!!!! zdsqdif_ssi_frosst need to be changed to an
+                !average
+              ELSEIF (h2o_ice_depth(ig,islope) .gt. 4e-4 .and. lag_layer
+     &        .and. watercaptag(ig)) then
+              zdqsdif_ssi_frost(ig,1)=(d_coef(ig,1)
+     &                         /h2o_ice_depth(ig,1))
+     &          *rho(ig)*dryness(ig)*(qsat(ig)-qeq(ig,1))
+               zdqsdif(ig)=zdqsdif(ig)-zdqsdif_ssi_frost(ig,1)
+               !needs to change to the mean of eq
+              ENDIF
+              call write_output('zdqsdif_ssi_frost',
+     &                '','',zdqsdif_ssi_frost(ig,1))
+              call write_output('subtimestep',
+     &                '','',subtimestep)
+             ! ENDDO !subsurface ice subslope
+
+
 
 c             Starting upward calculations for water :
 c             Actualisation de h2o_vap dans le premier niveau
              zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig)
-
 c    Take into account the H2O latent heat impact on the surface temperature
               if (latentheat_surfwater) then
@@ -1009,5 +1131,4 @@
      &                              *zq_tmp_vap(ig,ilay-1,iq)
               ENDDO
-
 c             Subtimestep water budget :
 
@@ -1016,5 +1137,7 @@
               zqsurf(ig)= zqsurf(ig)+(
      &                       zdqsdif(ig))*subtimestep
-
+              if (zqsurf(ig)<0) then 
+                      zqsurf(ig)=0
+              endif
 c             Monitoring instantaneous latent heat flux in W.m-2 :
               zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+
@@ -1033,5 +1156,4 @@
 c             Fin du sous pas de temps
             ENDDO ! tsub=1,nsubtimestep 
-
 c             Integration of subtimestep temp and water budget :
 c             (btw could also compute the post timestep temp and ice
@@ -1044,5 +1166,4 @@
      &                       cos(pi*def_slope_mean(islope)/180.))
      &                       /ptimestep
-
 c if subliming more than qsurf(ice) and on watercaptag, water
 c sublimates from watercap reservoir (dwatercap_dif is <0)
@@ -1063,5 +1184,4 @@
            ENDDO ! of DO ig=1,ngrid
           ENDDO ! islope
-
 c Some grid box averages: interface with the atmosphere
        DO ig = 1,ngrid
@@ -1075,5 +1195,4 @@
          ENDDO
        ENDDO
-
 ! Recompute values in kg/m^2 slopped
         DO ig = 1,ngrid
@@ -1167,5 +1286,10 @@
      &                          "Ground ice latent heat flux",
      &                               "W.m-2",surf_h2o_lh(:,iflat))
-
+                       call write_output('zdq_subtimestep',
+     &                '','',zdqsdif(:)*subtimestep)
+                       call write_output('zdq_end',
+     &                '','',zdqsdif(:))
+                       call write_output('vdifc_subtimestep',
+     &                '','',subtimestep)
 C       Diagnostic output for HDO
 !        if (hdo) then
