Changeset 3230 for trunk/LMDZ.MARS

Timestamp:

Feb 20, 2024, 10:14:22 PM (9 months ago)

Author:

llange

Message:

Mars PCM

• Move soil_tifeedback into a module waterice_tifeedback_mod.F90
• The flag to call tifeedback has been changed from "tifeedback" to surfaceice_tifeedback for clarity
• Add the possibility to change the thermal inertia while pore ices is forming in the soil: to do so, use the flag poreice_tifeedback. The computation is done in waterice_tifeedback_mod.F90.

For now, surfaceice_tifeedback and poreice_tifeedback can not be called together, we might think about how to merge later.
LL

Location:

trunk/LMDZ.MARS

Files:

• changelog.txt (modified) (1 diff)
• libf/phymars/callkeys.h (modified) (3 diffs)
• libf/phymars/comsoil_h.F90 (modified) (1 diff)
• libf/phymars/conf_phys.F (modified) (2 diffs)
• libf/phymars/physiq_mod.F (modified) (8 diffs)
• libf/phymars/soil.F (modified) (1 diff)
• libf/phymars/soil_tifeedback.F (deleted)
• libf/phymars/soilwater.F90 (modified) (4 diffs)
• libf/phymars/vdifc_mod.F (modified) (4 diffs)
• libf/phymars/waterice_tifeedback_mod.F90 (added)

trunk/LMDZ.MARS/changelog.txt

@@ -4511 +4511 @@
 Remove interactive checking with XIOS whether a field should be sent to it;
 some yet unresolved issues arise when using this in mixed MPI-OpenMP mode...
+
+== 20/02/2024 == LL
+* Move soil_tifeedback into a module waterice_tifeedback_mod.F90
+* The flag to call tifeedback has been changed from "tifeedback" to surfaceice_tifeedback for clarity
+* Add the possibility to change the thermal inertia while pore ices is forming in the soil: to do so, use the flag poreice_tifeedback. The computation is done i
+n waterice_tifeedback_mod.F90.
+For now, surfaceice_tifeedback and poreice_tifeedback can not be called together
+, we might think about how to merge later.

trunk/LMDZ.MARS/libf/phymars/callkeys.h

@@ -12 +12 @@
      &   ,callg2d,linear,rayleigh                                       &
      &   ,scavenging,sedimentation                                      &
-     &   ,activice,water,tifeedback,microphys,supersat,caps,photochem   &
+     &   ,activice,water,microphys,supersat,caps,photochem              &
      &   ,calltherm,callrichsl,callslope,tituscap,callyamada4,co2clouds &
      &   ,co2useh2o,meteo_flux,activeco2ice,CLFvaryingCO2,spantCO2      &
@@ -18 +18 @@
      &   ,latentheat_surfwater,gwd_convective_source,startphy_file      &
      &   ,hdo,hdofrac,cst_cap_albedo,temp_dependent_m,refill_watercap   &
-     &   ,cloud_adapt_ts, callatke
+     &   ,cloud_adapt_ts,callatke,surfaceice_tifeedback                 &
+     &   , poreice_tifeedback
 !$OMP THREADPRIVATE(/callkeys_l/)
 
@@ -79 +80 @@
       logical cloud_adapt_ts ! adaptative timestep for cloud microphysics
       logical sedimentation
-      logical activice,tifeedback,supersat,caps
+      logical activice,surfaceice_tifeedback,supersat,caps
+      logical poreice_tifeedback
       logical co2clouds,co2useh2o,meteo_flux,CLFvaryingCO2,satindexco2
       logical activeco2ice

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

@@ -30 +30 @@
 
 ! Subsurface tracers:
-logical, save                                  :: adsorption_soil ! boolean to call adosrption (or not)
-real,    save                                  :: choice_ads      ! Choice for adsorption isotherm (3 means no adsorption, see soilwater.F90)
-real,    save, allocatable, dimension(:,:,:,:) :: qsoil           ! subsurface tracers (kg/m^3 of regol)
-integer, parameter :: nqsoil = 3                                  ! number of subsurface tracers, only three when working with water
+logical, save                                  :: adsorption_soil      ! boolean to call adosrption (or not)
+real,    save                                  :: choice_ads           ! Choice for adsorption isotherm (3 means no adsorption, see soilwater.F90)
+real,    save, allocatable, dimension(:,:,:,:) :: qsoil                ! subsurface tracers (kg/m^3 of regol)
+integer, parameter :: nqsoil = 3                                       ! number of subsurface tracers, only three when working with water
 integer, parameter :: igcm_h2o_vap_soil = 1
 integer, parameter :: igcm_h2o_ice_soil = 2
 integer, parameter :: igcm_h2o_vap_ads  = 3
-
+REAL, parameter :: porosity_reg  = 0.45
 !$OMP THREADPRIVATE(adsorption_soil,qsoil,choice_ads)
 

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

@@ -683 +683 @@
 ! thermal inertia feedback
          write(*,*) "Activate the thermal inertia feedback ?"
-         tifeedback=.false. ! default value
-         call getin_p("tifeedback",tifeedback)
-         write(*,*) " tifeedback = ",tifeedback
+         surfaceice_tifeedback=.false. ! default value
+         call getin_p("surfaceice_tifeedback",surfaceice_tifeedback)
+         write(*,*) " surfaceice_tifeedback = ",surfaceice_tifeedback
 
 ! Test of incompatibility:
 
-         if (tifeedback.and..not.water) then
-           print*,'if tifeedback is used,'
+         if (surfaceice_tifeedback.and..not.water) then
+           print*,'if surfaceice_tifeedback is used,'
            print*,'water should be used too'
            call abort_physic(modname,
-     &          "tifeedback requires water",1)
-         endif
-
-         if (tifeedback.and..not.callsoil) then
-           print*,'if tifeedback is used,'
+     &          "surfaceice_tifeedback requires water",1)
+         endif
+
+         if (surfaceice_tifeedback.and..not.callsoil) then
+           print*,'if surfaceice_tifeedback is used,'
            print*,'callsoil should be used too'
            call abort_physic(modname,
-     &          "tifeedback requires callsoil",1)
+     &          "surfaceice_tifeedback requires callsoil",1)
          endif
 
@@ -1041 +1041 @@
          if (adsorption_soil) choice_ads = 1
          call getin_p("choice_ads",choice_ads) 
-           
+         
+         poreice_tifeedback = .false.
+         call getin_p("poreice_tifeedback",poreice_tifeedback)
+         write(*,*) 'poreice_tifeedback=',poreice_tifeedback 
+         if (poreice_tifeedback .and. (.not. adsorption_soil)) then
+              write(*,*)"Pore ice TI feedback can be run only if 
+     &                   adsorption_soil = True"
+              call abort_physic(modname,
+     &        "Pore ice TI feedback be used with adsorption_soil
+     &         = true",1)
+         endif 
 c        ----------------------------------------------------------
 

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

@@ -123 +123 @@
       use pbl_parameters_mod, only: pbl_parameters
       use lmdz_atke_turbulence_ini, only : atke_ini
+      use waterice_tifeedback_mod, only : waterice_tifeedback
       IMPLICIT NONE
 c=======================================================================
@@ -282 +283 @@
       REAL inertiesoil_tifeedback(ngrid,nsoilmx,nslope) ! Time varying subsurface
                                       ! thermal inertia (J.s-1/2.m-2.K-1)
-                                      ! (used only when tifeedback=.true.)
+                                      ! (used only when tifeedback surface or pore =.true.)
 c     Variables used by the CO2 clouds microphysical scheme:
       DOUBLE PRECISION riceco2(ngrid,nlayer)   ! co2 ice geometric mean radius (m)
@@ -540 +541 @@
 
 ! Variable for ice table
-      REAL :: rhowater_surf(ngrid,nslope)         ! Water density at the surface [kg/m^3]
-      REAL :: rhowater_surf_sat(ngrid,nslope)     ! Water density at the surface at saturation [kg/m^3]
-      REAL :: rhowater_soil(ngrid,nsoilmx,nslope) ! Water density in soil layers [kg/m^3]
-      REAL,PARAMETER :: alpha_clap_h2o = 28.9074  ! Coeff for Clapeyron law [/]
-      REAL,PARAMETER :: beta_clap_h2o = -6143.7   ! Coeff for Clapeyron law [K]
-      REAL :: pvap_surf(ngrid)                    ! Water vapor partial pressure in first layer [Pa]
-      REAL,PARAMETER :: m_co2 = 44.01E-3          ! CO2 molecular mass [kg/mol]
-      REAL,PARAMETER :: m_noco2 = 33.37E-3        ! Non condensible mol mass [kg/mol]
-      REAL :: ztmp1,ztmp2                         ! intermediate variables to compute the mean molar mass of the layer
-
+      REAL :: rhowater_surf(ngrid,nslope)            ! Water density at the surface [kg/m^3]
+      REAL :: rhowater_surf_sat(ngrid,nslope)        ! Water density at the surface at saturation [kg/m^3]
+      REAL :: rhowater_soil(ngrid,nsoilmx,nslope)    ! Water density in soil layers [kg/m^3]
+      REAL,PARAMETER :: alpha_clap_h2o = 28.9074     ! Coeff for Clapeyron law [/]
+      REAL,PARAMETER :: beta_clap_h2o = -6143.7      ! Coeff for Clapeyron law [K]
+      REAL :: pvap_surf(ngrid)                       ! Water vapor partial pressure in first layer [Pa]
+      REAL,PARAMETER :: m_co2 = 44.01E-3             ! CO2 molecular mass [kg/mol]
+      REAL,PARAMETER :: m_noco2 = 33.37E-3           ! Non condensible mol mass [kg/mol]
+      REAL :: ztmp1,ztmp2                            ! intermediate variables to compute the mean molar mass of the layer
+      REAL :: pore_icefraction(ngrid,nsoilmx,nslope) ! ice filling fraction in the pores
 ! Variable for the computation of the TKE with parameterization from ATKE working group
       REAL :: viscom                              ! kinematic molecular viscosity for momentum
@@ -694 +695 @@
          IF (callsoil) THEN
 c           Thermal inertia feedback:
-            IF (tifeedback) THEN
+            IF (surfaceice_tifeedback.or.poreice_tifeedback) THEN
              DO islope = 1,nslope
-                CALL soil_tifeedback(ngrid,nsoilmx,
-     s               qsurf(:,:,islope),
-     s               inertiesoil_tifeedback(:,:,islope))
+                CALL waterice_tifeedback(ngrid,nsoilmx,nslope,
+     s               qsurf(:,igcm_h2o_ice,:),pore_icefraction,
+     s               inertiesoil_tifeedback)
              ENDDO
                 CALL soil(ngrid,nsoilmx,firstcall,
@@ -1498 +1499 @@
 c        Calling vdif (Martian version WITH CO2 condensation)
          dwatercap_dif(:,:) = 0.
-
          CALL vdifc(ngrid,nlayer,nsoilmx,nq,nqsoil,zpopsk,
      $        ptimestep,capcal,lwrite,
      $        zplay,zplev,zzlay,zzlev,z0,
-     $        pu,pv,zh,pq,tsurf,tsoil,emis,qsurf,qsoil,
+     $        pu,pv,zh,pq,tsurf,tsoil,emis,qsurf,
+     $        qsoil,pore_icefraction,
      $        zdum1,zdum2,zdh,pdq,zflubid,
      $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
@@ -2385 +2386 @@
       IF (callsoil) THEN
 c       Thermal inertia feedback
-        IF (tifeedback) THEN
-         DO islope = 1,nslope
-           CALL soil_tifeedback(ngrid,nsoilmx,
-     s               qsurf(:,:,islope),
-     s               inertiesoil_tifeedback(:,:,islope))
-         ENDDO
+        IF (surfaceice_tifeedback.or.poreice_tifeedback) THEN
+        
+           CALL waterice_tifeedback(ngrid,nsoilmx,nslope,
+     s               qsurf(:,igcm_h2o_ice,:),pore_icefraction,
+     s               inertiesoil_tifeedback(:,:,:))
+     
          CALL soil(ngrid,nsoilmx,.false.,inertiesoil_tifeedback,
      s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
@@ -3551 +3552 @@
      &                         '',albedo(:,1,islope))
          ENDDO
-              if (tifeedback) then
+              if (surfaceice_tifeedback.or.poreice_tifeedback) then
                  call write_output("soiltemp",
      &                              "Soil temperature","K",
@@ -3558 +3559 @@
      &                       'Soil Thermal Inertia',
      &         'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,iflat))
+
          do islope=1,nslope
           write(str2(1:2),'(i2.2)') islope

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

@@ -41 +41 @@
 
 ! 0. Initialisations and preprocessing step
-      if (firstcall.or.tifeedback) then
+      if(firstcall.or.surfaceice_tifeedback.or.poreice_tifeedback) then
       ! note: firstcall is set to .true. or .false. by the caller
       !       and not changed by soil.F 

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

@@ -1 +1 @@
 subroutine soilwater(ngrid, nlayer, nq, nsoil, nqsoil, ptsrf, ptsoil, ptimestep,  & 
       exchange, qsat_surf, pq, pa, pb, pc, pd, pdqsdifpot, pqsurf,  & 
-      pqsoil, pplev, rhoatmo, writeoutput, zdqsdifrego, zq1temp2) 
-
-
-      use comsoil_h, only: igcm_h2o_vap_soil, igcm_h2o_ice_soil, igcm_h2o_vap_ads, layer, mlayer, choice_ads
+      pqsoil, pplev, rhoatmo, writeoutput, zdqsdifrego, zq1temp2, saturation_water_ice) 
+
+
+      use comsoil_h, only: igcm_h2o_vap_soil, igcm_h2o_ice_soil, igcm_h2o_vap_ads, layer, mlayer, choice_ads, porosity_reg
       use comcstfi_h
       use tracer_mod
@@ -92 +92 @@
 ! Outputs :
 ! ---------------------------------------------------------------------- 
-real, intent(out) :: zdqsdifrego(ngrid)         ! Flux from subsurface (positive pointing outwards)
-real, intent(out) :: zq1temp2(ngrid)            ! Temporary atmospheric mixing ratio after exchange with subsurface (kg / kg)
+real, intent(out) :: zdqsdifrego(ngrid)                     ! Flux from subsurface (positive pointing outwards)
+real, intent(out) :: zq1temp2(ngrid)                        ! Temporary atmospheric mixing ratio after exchange with subsurface (kg / kg)
+real*8, intent(out) :: saturation_water_ice(ngrid, nsoil)   ! Water pore ice saturation level (formerly Sw)
 
 ! Outputs for the output files
@@ -162 +163 @@
 real*8, allocatable, save :: tortuosity(:, :)         ! Tortuosity factor (formerly tortuo)
 
-real*8 :: saturation_water_ice(ngrid, nsoil)          ! Water pore ice saturation level (formerly Sw)
 real*8 :: saturation_water_ice_inter(ngrid, nsoil)          ! Water pore ice saturation level at the interlayer
 
@@ -346 +346 @@
                         
                         ! These properties are defined here in order to enable custom profiles
-                        porosity_ice_free(ig, ik) = 0.45D0
+                        porosity_ice_free(ig, ik) = porosity_reg
                         tortuosity(ig, ik) = 1.5D0
                         rho_soil(ig, ik) = 1.3D3 ! in kg/m3 of regolith (incl. porosity)

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

@@ -9 +9 @@
      $                pplay,pplev,pzlay,pzlev,pz0,
      $                pu,pv,ph,pq,ptsrf,ptsoil,pemis,pqsurf,qsoil,
-     $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
+     $                pore_icefraction,pdufi,pdvfi,pdhfi,
+     $                pdqfi,pfluxsrf,
      $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
      $                pdqdif,pdqsdif,wstar,
@@ -239 +240 @@
       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
-      
+      REAL, INTENT(out) :: pore_icefraction(ngrid,nsoil,nslope) ! ice filling fraction in the pores
 !! Water buyoncy
       LOGICAL :: virtual
@@ -307 +308 @@
       h2o_ice_depth(1:ngrid,1:nslope)=1
       virtual = .false.
-      
+      pore_icefraction(:,:,:) = 0.
 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
@@ -1084 +1085 @@
      &                     qsoil(ig,:,:,islope), pplev(ig,1), rho(ig), 
      &                     writeoutput,zdqsdif_regolith(ig,islope), 
-     &                     zq1temp_regolith(ig))
+     &                     zq1temp_regolith(ig),
+     &                     pore_icefraction(ig,:,islope))