Changeset 3121

Timestamp:

Nov 10, 2023, 2:12:31 PM (13 months ago)

Author:

llange

Message:

MARS PCM

• 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

LL

Location:

trunk/LMDZ.MARS

Files:

• changelog.txt (modified) (1 diff)
• libf/phymars/soilwater.F90 (modified) (11 diffs)
• libf/phymars/vdifc_mod.F (modified) (6 diffs)

trunk/LMDZ.MARS/changelog.txt

@@ -4325 +4325 @@
 == 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
+

trunk/LMDZ.MARS/libf/phymars/soilwater.F90

@@ -80 +80 @@
 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 +90 @@
 ! ----------------------------------------------------------------------
 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 +128 @@
 real*8 :: alpha0(ngrid) 
 real*8 :: beta0(ngrid) 
-real    :: co2ice(ngrid)                        ! CO2 ice on the surface
 
 ! Adsorbtion/Desorption variables and parameters
@@ -247 +246 @@
 integer :: sublim_n                       ! number of sublimation loop runs
 integer :: satur_n                        ! number of saturation loop runs
+
 
 
@@ -286 +286 @@
 ! ================= 
 
-co2ice =  pqsurf(:,igcm_co2)
 
 if (firstcall_soil) then
@@ -648 +647 @@
             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 +817 @@
                         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 +823 @@
                   
                   ! 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 +925 @@
                   ! 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 +932 @@
                   ! 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 +945 @@
             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

trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F

@@ -955 +955 @@
 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 +1026 @@
 
              if (adsorption_soil) then  
-
                  call soilwater(1,nlay,nq,nsoil, nqsoil,
      &                     ztsrf(ig),ptsoil(ig,:,islope),subtimestep,
@@ -1047 +1042 @@
                          zdqsdif_tot(ig)=
      &                        -zqsurf(ig)/subtimestep
-                        zdqsdif_tot(ig) = zdqsdif_surf(ig) 
                      else
                          zdqsdif_tot(ig) = zdqsdif_surf(ig) +
@@ -1055 +1049 @@
               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 +1179 @@
 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 +1192 @@
               ENDDO
 c             Subtimestep water budget :
-
               ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig,islope) 
      &                 + zdtsrf(ig,islope))*subtimestep