Index: trunk/LMDZ.MARS/changelog.txt
===================================================================
--- trunk/LMDZ.MARS/changelog.txt	(revision 3120)
+++ trunk/LMDZ.MARS/changelog.txt	(revision 3121)
@@ -4325,2 +4325,9 @@
 == 09/11/2023 == JBC + EV
 Correction of a bug introduced in r3059: in case of 'watercaptag=.true.', 'watercap' was prevented to evolve as it should be because of 'zqsurf' in "vdifc_mod.F90" which was forced to be 0 if negative.
+
+== 10/11/2023 == LL
+* Fixing index bug in soilwater.F The adsorption/desorption model now run
+with the subtimestep (but still not in 3D yet).
+* Small changes for the boundary conditions when considering
+adsorption/desorption
+
Index: trunk/LMDZ.MARS/libf/phymars/soilwater.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/soilwater.F90	(revision 3120)
+++ trunk/LMDZ.MARS/libf/phymars/soilwater.F90	(revision 3121)
@@ -80,5 +80,5 @@
 real, intent(in) :: pc(ngrid, nlayer)           ! Coefficients zc
 real, intent(in) :: pd(ngrid, nlayer)           ! Coefficients zd
-real, intent(in) :: pdqsdifpot(ngrid, nq)       ! Potential pdqsdif (without subsurface - atmosphere exchange)
+real, intent(in) :: pdqsdifpot(ngrid)       ! Potential pdqsdif (without subsurface - atmosphere exchange)
 
 real, intent(in) :: pplev(ngrid, nlayer+1)      ! Atmospheric pressure levels
@@ -90,5 +90,5 @@
 ! ----------------------------------------------------------------------
 real, intent(inout) :: pqsoil(ngrid, nsoil, nqsoil) ! Subsurface tracers (water vapor and ice) 
-real, intent(in) :: pqsurf(ngrid, nq)           ! Water ice on the surface
+real, intent(in) :: pqsurf(ngrid)           ! Water ice on the surface
                                                       ! (in kg.m - 3 : m - 3 of pore air for water vapor and m - 3 of regolith for water ice and adsorbed water)
 ! Outputs :
@@ -128,5 +128,4 @@
 real*8 :: alpha0(ngrid) 
 real*8 :: beta0(ngrid) 
-real    :: co2ice(ngrid)                        ! CO2 ice on the surface
 
 ! Adsorbtion/Desorption variables and parameters
@@ -247,4 +246,5 @@
 integer :: sublim_n                       ! number of sublimation loop runs
 integer :: satur_n                        ! number of saturation loop runs
+
 
 
@@ -286,5 +286,4 @@
 ! ================= 
 
-co2ice =  pqsurf(:,igcm_co2)
 
 if (firstcall_soil) then
@@ -648,5 +647,5 @@
             alpha0(ig) = porosity_ice_free(ig, 1) * pb(ig, 1) / delta0(ig)
             
-            beta0(ig) = ( pb(ig, 2) * pc(ig, 2) + pq(ig, 1, igcm_h2o_vap) * pa(ig, 1) + pqsurf(ig, igcm_h2o_ice) ) &
+            beta0(ig) = ( pb(ig, 2) * pc(ig, 2) + pq(ig, 1, igcm_h2o_vap) * pa(ig, 1) + pqsurf(ig) ) &
                   / delta0(ig)
                   
@@ -818,5 +817,5 @@
                         flux(ig, 0) = porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep &
                               * ( znsoil(ig, 1) / rho(ig) - beta0(ig) - alpha0(ig) * znsoil(ig, 1) / rho(ig) )   
-                  elseif (pqsurf(ig, igcm_h2o_ice).gt.0.) then
+                  elseif (pqsurf(ig).gt.0.) then
                         ! assume that the surface is covered by waterice (if it is co2ice it should not call this subroutine at all)
                         flux(ig, 0) = D_mid(ig, 1) / midlayer_dz(ig, 0) * ( znsoil(ig, 1) - qsat_surf(ig) * rhosurf(ig) ) 
@@ -824,5 +823,5 @@
                   
                   ! check if the ground would take up water from the surface but there is none
-                  if ((.not.exchange(ig)).and.(pqsurf(ig, igcm_h2o_ice).eq.0.).and.(flux(ig, 0).lt.0.)) then
+                  if ((.not.exchange(ig)).and.(pqsurf(ig).eq.0.).and.(flux(ig, 0).lt.0.)) then
                         flux(ig, 0) = 0.
                   endif
@@ -926,5 +925,5 @@
                   ! calculate the temporty mixing ratio above the surface
                   zq1temp2(ig) = beta0(ig) + alpha0(ig) * znsoil(ig, 1) / rho(ig)
-                  
+
                   ! calculate the flux from the subsurface
                   zdqsdifrego(ig) = porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep * (zq1temp2(ig) - znsoil(ig, 1) / rho(ig) )
@@ -933,5 +932,5 @@
                   ! calculate the flux from the subsurface
                   zdqsdifrego(ig) = -D_mid(ig, 1) / midlayer_dz(ig, 0) * (qsat_surf(ig) * rhosurf(ig) - znsoil(ig, 1)) ! faut - il faire intervenir la porosite ?
-                  if ((pqsurf(ig, igcm_h2o_ice).eq.0.).and.(zdqsdifrego(ig).lt.0.)) then
+                  if ((pqsurf(ig).eq.0.).and.(zdqsdifrego(ig).lt.0.)) then
                         zdqsdifrego(ig) = 0.
                   endif
@@ -946,9 +945,9 @@
             if ( (.not.exchange(ig)) &
             .and. ( -(zdqsdifrego(ig) * ptimestep) &
-            .gt.( pqsurf(ig, igcm_h2o_ice) + pdqsdifpot(ig, igcm_h2o_ice) * ptimestep) ) &
-            .and.( (pqsurf(ig, igcm_h2o_ice) + pdqsdifpot(ig, igcm_h2o_ice) * ptimestep).gt.0. ) ) then
+            .gt.( pqsurf(ig)  + pdqsdifpot(ig) * ptimestep) ) &
+            .and.( (pqsurf(ig) + pdqsdifpot(ig) * ptimestep).gt.0. ) ) then
 
                   ! calculate a new flux from the subsurface
-                  zdqsdifrego(ig) = -( pqsurf(ig, igcm_h2o_ice) + pdqsdifpot(ig, igcm_h2o_ice) * ptimestep ) / ptimestep
+                  zdqsdifrego(ig) = -( pqsurf(ig) + pdqsdifpot(ig) * ptimestep ) / ptimestep
                   
 !                  ! check case where there is CO2 ice on the surface but qsurf = 0
Index: trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 3120)
+++ trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 3121)
@@ -955,11 +955,7 @@
 c          make_tsub : sous pas de temps adaptatif
 c          la subroutine est a la fin du fichier
-           if (adsorption_soil) then  
-              nsubtimestep(:) = 1
-           else 
               call make_tsub(ngrid,pdtsrf(:,islope),zqsurf,
      &                    ptimestep,dtmax,watercaptag,
      &                    nsubtimestep)
-           endif
 c           Calculation for turbulent exchange with the surface (for ice)
 c           initialization of ztsrf, which is surface temperature in
@@ -1030,5 +1026,4 @@
 
              if (adsorption_soil) then  
-
                  call soilwater(1,nlay,nq,nsoil, nqsoil,
      &                     ztsrf(ig),ptsoil(ig,:,islope),subtimestep,
@@ -1047,5 +1042,4 @@
                          zdqsdif_tot(ig)=
      &                        -zqsurf(ig)/subtimestep
-                        zdqsdif_tot(ig) = zdqsdif_surf(ig) 
                      else
                          zdqsdif_tot(ig) = zdqsdif_surf(ig) +
@@ -1055,5 +1049,5 @@
               endif ! adsorption
 
-             if(.not.watercaptag(ig)) then
+             if(.not.watercaptag(ig).and.(.not.adsorption_soil)) then
                if ((-zdqsdif_tot(ig)*subtimestep)
      &            .gt.(zqsurf(ig))) then
@@ -1185,5 +1179,5 @@
 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
+c             Take into account the H2O latent heat impact on the surface temperature
               if (latentheat_surfwater) then
                 lh=(2834.3-0.28*(ztsrf(ig)-To)-
@@ -1198,5 +1192,4 @@
               ENDDO
 c             Subtimestep water budget :
-
               ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig,islope) 
      &                 + zdtsrf(ig,islope))*subtimestep
