Changeset 3115 for trunk

Timestamp:

Nov 3, 2023, 6:02:56 PM (13 months ago)

Author:

llange

Message:

Mars PCM

• Introducing the possibily to compute water adsorption / desorption -routine soilwater.F90) -. FOR NOW IT WORKS ONLY IN 1D, DON'T TEST IN 3D (by default, adsorption is not called, if not specified in the callphys.def). By default, when using the adsorption, the number of subtimestep used in the water subimation scheme is 1 (otherwise, it crashes)

• New grid for the soil layers (better resolution) as it is needed to solve the diffusion equations. It does not increase the CPU time.

PYM & LL

Location:

trunk/LMDZ.MARS

Files:

• changelog.txt (modified) (1 diff)
• deftank/field_def_physics_mars.xml (modified) (2 diffs)
• libf/phymars/comsoil_h.F90 (modified) (1 diff)
• libf/phymars/dyn1d/init_testphys1d_mod.F90 (modified) (1 diff)
• libf/phymars/physiq_mod.F (modified) (1 diff)
• libf/phymars/soil_settings.F (modified) (1 diff)
• libf/phymars/soilwater.F90 (added)
• libf/phymars/vdifc_mod.F (modified) (22 diffs)

trunk/LMDZ.MARS/changelog.txt

@@ -4310 +4310 @@
 Introducing qsoil to model H2O adsorption/desorption in the subsurface. For now, I've fixed the number of tracers in the subsurface to t
 hree (H2O vap, H2O ice, H2O ads).
-
+Changing the soil grid (better resolution) as given by Pierre Yves Meslin
+Adding the possibility to compute H2O adsorption desorption and exchange with
+the near-surface atmosphere, only works in 1D for now. By default, subtimestep
+for water sublimation is 1 when adsorption is activated (otherwise it crashes)
+
+

trunk/LMDZ.MARS/deftank/field_def_physics_mars.xml

@@ -695 +695 @@
                    unit="m/s" />
 
+            <!-- Subsurface tracers  (adsorption) -->
+
+            <field id="flux_rego"
+                   long_name="flux of water from the regolith"
+                   unit="kg/m^2" />
+            <field id="mass_h2o_soil"
+                   long_name="Mass of subsurface water column at each point"
+                   unit="kg m-2" />
+            <field id="mass_ice_soil"
+                   long_name="Mass of subsurface ice at each point"
+                   unit="kg m-2" />
+            <field id="nsurf"
+                   long_name="surface water vapor density"
+                   unit="kg/m^3" />
+
+
         </field_group>
 
@@ -1227 +1243 @@
                    long_name="Flux after all contributions"
                    unit="kg.m-2.s-1" />
+            <field id="flux_soillayer"
+                   long_name="flux of water between the soil layers"
+                   unit="kg.m-2.s-1" />
+            <field id="ice_saturation_soil"
+                   long_name="Water ice saturation in the soil layers"
+                   unit="Percent" />
+            <field id="znsoil"
+                   long_name="Water vapor soil concentration"
+                   unit="kg m-3 of pore air" />
+            <field id="nsatsoil"
+                   long_name="subsurface water vapor saturation density"
+                   unit="kg/m^3" />
+            <field id="adswater"
+                   long_name="subsurface adsorbed water"
+                   unit="kg/m^3" />
+            <field id="coeff_diffusion_soil"
+                   long_name="interlayer diffusion coefficient"
+                   unit="m^2/s" />
+
 
         </field_group>

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

@@ -3 +3 @@
 implicit none
 ! nsoilmx : number of subterranean layers
-  integer, parameter :: nsoilmx = 18 
+  integer, parameter :: nsoilmx = 57 
 
   real,save,allocatable,dimension(:) :: layer      ! soil layer depths

trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90

@@ -670 +670 @@
 ! -----------------
 do isoil = 0,nsoil - 1
-    mlayer(isoil) = 2.e-4*(2.**(isoil - 0.5)) ! mid-layer depth
-    layer(isoil + 1) = 2.e-4*(2.**isoil)      ! layer depth
+    mlayer(isoil) = 2.e-4*(1+isoil**2.9*(1-exp(-real(isoil)/20.))) ! mid layer depth
 enddo
+
+do isoil = 1,nsoil - 1
+    layer(isoil)=(mlayer(isoil)+mlayer(isoil-1))/2
+enddo
+layer(nsoil)=2*mlayer(nsoil-1) - mlayer(nsoil-2)
 
 ! Initialize traceurs

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

@@ -1507 +1507 @@
          zcdv(:) = 0.
 
-         CALL vdifc(ngrid,nlayer,nq,zpopsk,
+         CALL vdifc(ngrid,nlayer,nsoilmx,nq,nqsoil,zpopsk,
      $        ptimestep,capcal,lwrite,
      $        zplay,zplev,zzlay,zzlev,z0,
-     $        pu,pv,zh,pq,tsurf,emis,qsurf,
+     $        pu,pv,zh,pq,tsurf,tsoil,emis,qsurf,qsoil,
      $        zdum1,zdum2,zdh,pdq,zflubid,
      $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
      &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,
      &        zcondicea_co2microp,sensibFlux,
-     &        dustliftday,local_time,watercap,dwatercap_dif,
-     &        tsoil,nsoilmx)
+     &        dustliftday,local_time,watercap,dwatercap_dif)
 
           DO ig=1,ngrid

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

@@ -148 +148 @@
 ! 1.4 Build mlayer(), if necessary
       if (interpol) then
-      ! default mlayer distribution, following a power law:
-      !  mlayer(k)=lay1*alpha**(k-1/2)
-         lay1=2.e-4
-         alpha=2
-         do iloop=0,nsoil-1
-           mlayer(iloop)=lay1*(alpha**(iloop-0.5))
-         enddo
+      !  mlayer(k)=lay1*(1+k^(2.9)*(1-exp(-k/20)), PYM grid
+        lay1=2.e-4
+        do iloop=0,nsoil-1
+          mlayer(iloop) = lay1*(1+iloop**2.9*(1-exp(-real(iloop)/20.)))
+        enddo
       endif
 
-! 1.5 Build layer(); following the same law as mlayer()
-      ! Assuming layer distribution follows mid-layer law:
-      ! layer(k)=lay1*alpha**(k-1)
-      lay1=sqrt(mlayer(0)*mlayer(1))
-      alpha=mlayer(1)/mlayer(0)
-      do iloop=1,nsoil
-        layer(iloop)=lay1*(alpha**(iloop-1))
+! 1.5 Build layer(); 
+      do iloop=1,nsoil-1
+        layer(iloop)=(mlayer(iloop)+mlayer(iloop-1))/2
       enddo
-
+      layer(nsoil)=2*mlayer(nsoil-1) - mlayer(nsoil-2)
 
 ! 2. Volumetric heat capacity (note: it is declared in comsoil_h)

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

@@ -5 +5 @@
       CONTAINS
       
-      SUBROUTINE vdifc(ngrid,nlay,nq,ppopsk,
+      SUBROUTINE vdifc(ngrid,nlay,nsoil,nq,nqsoil,ppopsk,
      $                ptimestep,pcapcal,lecrit,
      $                pplay,pplev,pzlay,pzlev,pz0,
-     $                pu,pv,ph,pq,ptsrf,pemis,pqsurf,
+     $                pu,pv,ph,pq,ptsrf,ptsoil,pemis,pqsurf,qsoil,
      $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
      $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
      $                pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,
      $                hfmax,pcondicea_co2microp,sensibFlux,
-     $                dustliftday,local_time,watercap, dwatercap_dif,
-     $                tsoil,nsoil)
+     $                dustliftday,local_time,watercap, dwatercap_dif)
+
       use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
      &                      igcm_dust_submicron, igcm_h2o_vap,
@@ -34 +34 @@
       use microphys_h, only: To
       use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer
-      use comsoil_h, only: layer, mlayer
+      use comsoil_h, only: layer, mlayer,adsorption_soil
 
       IMPLICIT NONE
@@ -64 +64 @@
 c   ----------
 
-      INTEGER,INTENT(IN) :: ngrid,nlay,nsoil
+      INTEGER,INTENT(IN) :: ngrid,nlay,nsoil,nqsoil
       REAL,INTENT(IN) :: ptimestep
       REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1)
@@ -78 +78 @@
       REAL,INTENT(IN) :: pcapcal(ngrid,nslope)
       REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1)
-      REAL,INTENT(IN) :: tsoil(ngrid,nsoil,nslope)
+      REAL,INTENT(IN) :: ptsoil(ngrid,nsoil,nslope)
 
 c    Argument added for condensation:
@@ -99 +99 @@
       real,intent(out) :: pdqsdif(ngrid,nq,nslope)
       REAL,INTENT(in) :: dustliftday(ngrid)
+      REAL,INTENT(inout) :: qsoil(ngrid,nsoil,nqsoil,nslope) !subsurface tracers
       REAL,INTENT(in) :: local_time(ngrid)
       
@@ -147 +148 @@
 c     Subtimestep & implicit treatment of water vapor
 c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
-      REAL zdqsdif(ngrid) ! subtimestep pdqsdif for water ice
+      REAL zdqsdif_surf(ngrid) ! subtimestep pdqsdif for water ice
       REAL ztsrf(ngrid) ! temporary surface temperature in tsub
       REAL zdtsrf(ngrid,nslope) ! surface temperature tendancy in tsub
@@ -227 +228 @@
 
       character*2 str2
+
+!! Subsurface exchanges
+      LOGICAL :: exchange             ! boolean to check if exchange between the subsurface and the atmosphere can occurs
+      REAL :: zdqsdif_regolith(ngrid,nslope) ! Flux from subsurface (positive pointing outwards) (kg/m^2/s)
+      REAL zq1temp_regolith(ngrid)    ! Temporary atmospheric mixing ratio after exchange with subsurface (kg / kg)
+      REAL zdqsdif_tot(ngrid)         ! subtimestep pdqsdif for water ice
+      LOGICAL :: writeoutput          ! boolean to say to soilexchange.F if we are at the last iteration and thus if he  can write in the diagsoil
 
 c    ** un petit test de coherence
@@ -283 +291 @@
       pdqsdif(1:ngrid,1:nq,1:nslope)=0
       pdqsdif_tmp(1:ngrid,1:nq)=0
-      zdqsdif(1:ngrid)=0
+      zdqsdif_surf(1:ngrid)=0
       dwatercap_dif(1:ngrid,1:nslope)=0
+      zdqsdif_regolith(1:ngrid,1:nslope)=0
+      zq1temp_regolith(1:ngrid)=0 
+      zdqsdif_tot(1:ngrid)=0
 !      h2o_ice_depth(1:ngrid,1:nslope)=5
 c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
@@ -944 +955 @@
 c          make_tsub : sous pas de temps adaptatif
 c          la subroutine est a la fin du fichier
-
-           call make_tsub(ngrid,pdtsrf(:,islope),zqsurf,
+           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
@@ -959 +972 @@
 c           Debut du sous pas de temps
             DO tsub=1,nsubtimestep(ig)
-
+              if(tsub.eq.nsubtimestep(ig)) writeoutput = .true.
 c           C'est parti ! 
              zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
@@ -992 +1005 @@
              zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
              if(old_wsublimation_scheme) then
-                zdqsdif(ig)=rho(ig)*dryness(ig)*zcdv(ig)
+                zdqsdif_surf(ig)=rho(ig)*dryness(ig)*zcdv(ig)
      &                           *(zq1temp(ig)-qsat(ig))
              else
-                zdqsdif(ig)=rho(ig)*dryness(ig)*zcdh(ig)
+                zdqsdif_surf(ig)=rho(ig)*dryness(ig)*zcdh(ig)
      &                           *(zq1temp(ig)-qsat(ig))
              endif
-c             write(*,*)'subliming more than available frost:  qsurf!'
-
+
+
+!!! Subsurface exchange
+! Check for subsurface exchanges
              if(.not.watercaptag(ig)) then
-               if ((-zdqsdif(ig)*subtimestep)
+                if (((-(zdqsdif_surf(ig))*
+     &             subtimestep).gt.zqsurf(ig)) 
+     &             .and.(pqsurf(ig,igcm_co2,islope).eq.0.)) then
+                      exchange = .true.
+                 else
+                      exchange = .false.
+                 endif
+             else
+                  exchange = .false.
+             endif
+                  zdqsdif_tot(ig) = zdqsdif_surf(ig)
+
+
+             if (adsorption_soil) then  
+
+                 call soilwater(1,nlay,nq,nsoil, nqsoil,
+     &                     ztsrf(ig),ptsoil(ig,:,islope),subtimestep,
+     &                     exchange,qsat(ig),zq_tmp_vap(ig,:,:),
+     &                     za(ig,:),zb(ig,:),zc(ig,:),zd(ig,:), 
+     &                     zdqsdif_surf(ig), zqsurf(ig), 
+     &                     qsoil(ig,:,:,islope), pplev(ig,1), rho(ig), 
+     &                     writeoutput,zdqsdif_regolith(ig,islope), 
+     &                     zq1temp_regolith(ig))
+
+
+
+                 if(.not.watercaptag(ig)) then
+                     if (exchange) then
+                         zq1temp(ig) = zq1temp_regolith(ig)
+                         zdqsdif_tot(ig)=
+     &                        -zqsurf(ig)/subtimestep
+                        zdqsdif_tot(ig) = zdqsdif_surf(ig) 
+                     else
+                         zdqsdif_tot(ig) = zdqsdif_surf(ig) +
+     &                           zdqsdif_regolith(ig,islope) ! boundary condition = qsat, but pdqsdif is calculated to update qsurf (including loss of surface ice to the subsurface)
+                     endif  ! of "if exchange = true"
+                 endif  ! of "if not.watercaptag"
+              endif ! adsorption
+
+             if(.not.watercaptag(ig)) then
+               if ((-zdqsdif_tot(ig)*subtimestep)
      &            .gt.(zqsurf(ig))) then
                 
@@ -1008 +1063 @@
                 IF(h2o_ice_depth(ig,islope) .gt. 0 .and. lag_layer)
      &           then 
-                zdqsdif(ig)=
+                zdqsdif_tot(ig)=
      &                        -zqsurf(ig)/subtimestep
                 zqsurf(ig)=0
@@ -1027 +1082 @@
      &                               +dist_down(ig,islope)
                  Tice(ig,islope)=(dist_up(ig,islope) ! Linear interp to calculate the temp
-     &                          *tsoil(ig,lice-1,islope)
+     &                          *ptsoil(ig,lice-1,islope)
      &           /dist_sum(ig,islope))+
-     &           (dist_down(ig,islope)*tsoil(ig,lice,islope)
+     &           (dist_down(ig,islope)*ptsoil(ig,lice,islope)
      &           /dist_sum(ig,islope))
                         ELSE
-                                Tice(ig,islope)=tsoil(ig,1,islope)
+                                Tice(ig,islope)=ptsoil(ig,1,islope)
                         ENDIF
               call watersat(1,Tice(ig,1),pplev(ig,1)
@@ -1079 +1134 @@
              !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)
+               !zdqsdif_surf(ig)=zdqsdif_ssi(ig,1)-(zqsurf(ig)/subtimestep)
+               zdqsdif_tot(ig)=zdqsdif_tot(ig)+zdqsdif_ssi(ig,1)
                 call write_output('zdq_zdqssi',
-     &                '','',zdqsdif(ig)+zdqsdif_ssi(ig,1))
+     &                '','',zdqsdif_tot(ig)+zdqsdif_ssi(ig,1))
              endif
               ELSE
@@ -1088 +1143 @@
 c             pdqsdif < 0 : ice subliming
 c             write(*,*) "subliming more than available frost:  qsurf!"
-                  zdqsdif(ig)=
+                  zdqsdif_tot(ig)=
      &                        -zqsurf(ig)/subtimestep
 c                write(*,*)'flux vers le sol=',pdqsdif(ig,nq)
@@ -1094 +1149 @@
                  zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,igcm_h2o_vap)+
      $           zb(ig,2)*zc(ig,2) +
-     $           (-zdqsdif(ig)) *subtimestep) *z1(ig)
+     $           (-zdqsdif_tot(ig)) *subtimestep) *z1(ig)
                  zq1temp(ig)=zc(ig,1)
                  ENDIF !if h2o_ice_depth>0 and lag_layer
@@ -1107 +1162 @@
      &          *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) 
+               zdqsdif_tot(ig)=zdqsdif_tot(ig)-zdqsdif_ssi_frost(ig,1) 
                 !!!!! zdsqdif_ssi_frosst need to be changed to an
                 !average
@@ -1115 +1170 @@
      &                         /h2o_ice_depth(ig,1))
      &          *rho(ig)*dryness(ig)*(qsat(ig)-qeq(ig,1))
-               zdqsdif(ig)=zdqsdif(ig)-zdqsdif_ssi_frost(ig,1)
+               zdqsdif_tot(ig)=zdqsdif_tot(ig)-zdqsdif_ssi_frost(ig,1)
                !needs to change to the mean of eq
               ENDIF
@@ -1134 +1189 @@
                 lh=(2834.3-0.28*(ztsrf(ig)-To)-
      &              0.004*(ztsrf(ig)-To)*(ztsrf(ig)-To))*1.e+3
-                zdtsrf(ig,islope)=  zdqsdif(ig)*lh /pcapcal(ig,islope)
+                zdtsrf(ig,islope)=  zdqsdif_tot(ig)*lh 
+     &                              /pcapcal(ig,islope)
               endif ! (latentheat_surfwater) then
 
@@ -1146 +1202 @@
      &                 + zdtsrf(ig,islope))*subtimestep
               zqsurf(ig)= zqsurf(ig)+(
-     &                       zdqsdif(ig))*subtimestep
+     &                       zdqsdif_tot(ig))*subtimestep
               if (zqsurf(ig)<0) then 
                       zqsurf(ig)=0
@@ -1158 +1214 @@
 c             is still ice on the surface (oldschool)
                if(zqsurf(ig)
-     &           +zdqsdif(ig)*subtimestep
+     &           +zdqsdif_tot(ig)*subtimestep
      &           .gt.frost_albedo_threshold) then ! if there is still ice, T cannot exceed To
                  zdtsrf(ig,islope) = min(zdtsrf(ig,islope),
@@ -1297 +1353 @@
      &                     "W.m-2",surf_h2o_lh(:,iflat))
          call write_output('zdq_subtimestep',
-     &                     'Actual flux zdqsdif*subtimestep',
-     &                     'kg.m-2',zdqsdif(:)*subtimestep)
+     &                     'Actual flux zdqsdif_surf*subtimestep',
+     &                     'kg.m-2',zdqsdif_surf(:)*subtimestep)
          call write_output('zdq_end',
      &                     'Flux after all contributions',
-     &                     'kg.m-2.s-1',zdqsdif(:))
+     &                     'kg.m-2.s-1',zdqsdif_surf(:))
 C       Diagnostic output for HDO
 !        if (hdo) then