Changeset 3113

Timestamp:

Nov 3, 2023, 2:36:01 PM (14 months ago)

Author:

llange

Message:

Mars PCM
Introducing qsoil to model H2O adsorption/desorption in the subsurface. For now, I've fixed the number of tracers in the subsurface to three (H2O vap, H2O ice, H2O ads).
The model of adsorption/desorption will follow later.
LL

Location:

trunk/LMDZ.MARS

Files:

• changelog.txt (modified) (1 diff)
• libf/phymars/comsoil_h.F90 (modified) (5 diffs)
• libf/phymars/conf_phys.F (modified) (2 diffs)
• libf/phymars/dyn1d/init_testphys1d_mod.F90 (modified) (2 diffs)
• libf/phymars/dyn1d/testphys1d.F90 (modified) (2 diffs)
• libf/phymars/phyetat0_mod.F90 (modified) (4 diffs)
• libf/phymars/phyredem.F90 (modified) (5 diffs)
• libf/phymars/physiq_mod.F (modified) (3 diffs)
• libf/phymars/soil_settings.F (modified) (5 diffs)

trunk/LMDZ.MARS/changelog.txt

@@ -4306 +4306 @@
 == 03/11/2023 == JBC
 Small fix on "field_def_physics_mars.xml".
+
+== 03/11/2023 == LL
+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).
+

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

@@ -20 +20 @@
   real,save,allocatable :: mthermdiff(:,:,:)  ! (FC) mid-layer thermal diffusivity
   real,save,allocatable :: thermdiff(:,:,:)   ! (FC) inter-layer thermal diffusivity
-  real,save,allocatable :: coefq(:)         ! (FC) q_{k+1/2} coefficients
+  real,save,allocatable :: coefq(:)           ! (FC) q_{k+1/2} coefficients
   real,save,allocatable :: coefd(:,:,:)       ! (FC) d_k coefficients
   real,save,allocatable :: alph(:,:,:)        ! (FC) alpha_k coefficients
   real,save,allocatable :: beta(:,:,:)        ! beta_k coefficients
   real,save :: mu
-  real,save,allocatable :: flux_geo(:,:)       ! Geothermal Flux (W/m^2)
+  real,save,allocatable :: flux_geo(:,:)      ! Geothermal Flux (W/m^2)
 
 !$OMP THREADPRIVATE(tsoil,mthermdiff,thermdiff,coefq,coefd,alph,beta,mu,flux_geo)
+
+
+! Subsurface tracers:
+  logical,save :: adsorption_soil             ! boolean to call adosrption (or not)
+  integer, parameter :: nqsoil = 3            ! number of subsurface tracers, only three when working with water
+  real,save,allocatable :: qsoil(:,:,:,:)     ! subsurface tracers (kg/m^3 of regol)
+  integer, parameter :: igcm_h2o_vap_soil = 1
+  integer, parameter :: igcm_h2o_ice_soil = 2
+  integer, parameter :: igcm_h2o_vap_ads  = 3
+!$OMP THREADPRIVATE(adsorption_soil,qsoil)
+
 
 contains
@@ -48 +59 @@
     allocate(beta(ngrid,nsoilmx-1,nslope))
     allocate(flux_geo(ngrid,nslope))
+    allocate(qsoil(ngrid,nsoilmx,nqsoil,nslope))
  
   end subroutine ini_comsoil_h
@@ -68 +80 @@
     if (allocated(beta)) deallocate(beta)
     if (allocated(flux_geo)) deallocate(flux_geo)
+    if (allocated(qsoil))  deallocate(qsoil)
   end subroutine end_comsoil_h
 
@@ -84 +97 @@
     allocate(beta(ngrid,nsoilmx-1,nslope))
     allocate(flux_geo(ngrid,nslope))
+    allocate(qsoil(ngrid,nsoilmx,nqsoil,nslope))
  
   end subroutine ini_comsoil_h_slope_var
@@ -100 +114 @@
     if (allocated(beta)) deallocate(beta)
     if (allocated(flux_geo)) deallocate(flux_geo)
+    if (allocated(qsoil))  deallocate(qsoil)
 
   end subroutine end_comsoil_h_slope_var

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

@@ -40 +40 @@
      &                           lag_layer 
       use microphys_h, only: mteta
+      use comsoil_h, only: adsorption_soil
 
       IMPLICIT NONE
@@ -1000 +1001 @@
 c        ----------------------------------------------------------
 
+! Adsorption
+         adsorption_soil = .false.
+         call getin_p("adsorption_soil",adsorption_soil) 
+         if (adsorption_soil .and. (.not. water)) then
+              write(*,*)"Adsorption can be run only if water = True"
+              call abort_physic(modname,
+     &        "Adsorption must be used with water = true",1)
+         endif  
+c        ----------------------------------------------------------
+
 ! THERMOSPHERE
 

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

@@ -20 +20 @@
 use iostart,                  only: open_startphy, get_var, close_startphy
 use physics_distribution_mod, only: init_physics_distribution
-use comsoil_h,                only: volcapa, nsoilmx, inertiesoil, inertiedat, layer, mlayer, flux_geo, tsoil
+use comsoil_h,                only: volcapa, nsoilmx, inertiesoil, inertiedat, layer, mlayer, flux_geo, tsoil, qsoil
 use comvert_mod,              only: ap, bp, aps, bps, pa, preff, presnivs, pseudoalt, scaleheight
 use dimradmars_mod,           only: tauvis, totcloudfrac, albedo
@@ -661 +661 @@
 flux_geo(:,:) = flux_geo_tmp
 
+! Initialize soil content
+! -----------------
+if (.not. startfiles_1D) then
+    qsoil(:,:,:,:) = 0.
+endif
+
 ! Initialize depths
 ! -----------------

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

@@ -1 +1 @@
 PROGRAM testphys1d
 
-use comsoil_h,           only: inertiedat, inertiesoil, nsoilmx, tsoil
+use comsoil_h,           only: inertiedat, inertiesoil, nsoilmx, tsoil, nqsoil, qsoil
 use surfdat_h,           only: albedodat, perenial_co2ice, watercap, tsurf, emis, qsurf
 use comslope_mod,        only: def_slope, subslope_dist
@@ -127 +127 @@
                   llm,nq,dttestphys,float(day0),0.,cell_area,    &
                   albedodat,inertiedat,def_slope,subslope_dist)
-    call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,dttestphys,time,       &
-                  tsurf,tsoil,inertiesoil,albedo,emis,q2,qsurf,tauscaling, &
+    call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,nqsoil,dttestphys,time,       &
+                  tsurf,tsoil,inertiesoil,albedo,emis,q2,qsurf,qsoil,tauscaling, &
                   totcloudfrac,wstar,watercap,perenial_co2ice)
 endif !(.not. therestartfi)

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

@@ -8 +8 @@
 contains
 
-subroutine phyetat0 (fichnom,tab0,Lmodif,nsoil,ngrid,nlay,nq, &
-                     day_ini,time0,tsurf,tsoil,albedo,emis,q2,qsurf, &
+subroutine phyetat0 (fichnom,tab0,Lmodif,nsoil,ngrid,nlay,nq,nqsoil, &
+                     day_ini,time0,tsurf,tsoil,albedo,emis,q2,qsurf,qsoil, &
                      tauscaling,totcloudfrac,wstar,watercap,perenial_co2ice, &
                      def_slope,def_slope_mean,subslope_dist)
@@ -59 +59 @@
   integer,intent(in) :: nlay ! # of atmospheric layers
   integer,intent(in) :: nq
+  integer,intent(in) :: nqsoil ! # of tracers in the soil
   integer :: day_ini
   real :: time0
@@ -69 +70 @@
   real,intent(out) :: q2(ngrid,nlay+1) ! 
   real,intent(out) :: qsurf(ngrid,nq,nslope) ! tracers on surface
+  real,intent(out) :: qsoil(ngrid,nsoil,nqsoil,nslope) ! tracers in the subsurface
   real,intent(out) :: tauscaling(ngrid) ! dust conversion factor
   real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
@@ -714 +716 @@
   ! Call to soil_settings, in order to read soil temperatures,
   ! as well as thermal inertia and volumetric heat capacity
-  call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
+  call soil_settings(nid_start,ngrid,nsoil,nqsoil,tsurf,tsoil,qsoil,indextime)
 else
     flux_geo(:,:) = 0.

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

@@ -164 +164 @@
 end subroutine physdem0
 
-subroutine physdem1(filename,nsoil,ngrid,nlay,nq, &
+subroutine physdem1(filename,nsoil,ngrid,nlay,nq,nqsoil, &
                     phystep,time,tsurf,tsoil,inertiesoil, &
-                    albedo,emis,q2,qsurf,&
+                    albedo,emis,q2,qsurf,qsoil,&
                     tauscaling,totcloudfrac,wstar, &
                     watercap,perenial_co2ice)
@@ -177 +177 @@
   use dust_rad_adjust_mod, only: dust_rad_adjust_prev,dust_rad_adjust_next
   use dust_param_mod, only: dustscaling_mode
-  use comsoil_h,only: flux_geo
+  use comsoil_h,only: flux_geo,adsorption_soil,igcm_h2o_vap_soil, &
+                      igcm_h2o_ice_soil,igcm_h2o_vap_ads
   use comslope_mod, only: nslope
   use paleoclimate_mod, only: paleoclimate
@@ -189 +190 @@
   integer,intent(in) :: nlay
   integer,intent(in) :: nq
+  integer,intent(in) :: nqsoil
   real,intent(in) :: phystep
   real,intent(in) :: time
@@ -198 +200 @@
   real,intent(in) :: q2(ngrid,nlay+1)
   real,intent(in) :: qsurf(ngrid,nq,nslope)
+  real,intent(in) :: qsoil(ngrid,nsoil,nqsoil,nslope)
   real,intent(in) :: tauscaling(ngrid)
   real,intent(in) :: totcloudfrac(ngrid)
@@ -339 +342 @@
   endif
 
-
   ! Geothermal Flux
      call put_field('flux_geo','Geothermal flux',flux_geo,time)
+
+  ! Adsorption
+  if (adsorption_soil) then
+     call put_field("h2o_vap_soil","subsurface water vapour", &
+                        qsoil(:,:,igcm_h2o_vap_soil,:), time)
+     call put_field("h2o_ice_soil","subsurface water ice", &
+                        qsoil(:,:,igcm_h2o_ice_soil,:), time)
+     call put_field("h2o_vap_ads", "adsorbed water", &
+                        qsoil(:,:,igcm_h2o_vap_ads,:), time)
+  endif
+
   ! Close file
   call close_restartphy

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

@@ -48 +48 @@
       use comsoil_h, only: inertiedat, inertiesoil,! dat: soil thermal inertia for present climate, inertiesoil is the TI read in the start
      &                     tsoil, nsoilmx,!number of subsurface layers
-     &                     mlayer,layer ! soil mid layer depths
+     &                     mlayer,layer, ! soil mid layer depths
+     &                     nqsoil,qsoil  ! adsorption
       use geometry_mod, only: longitude, latitude, cell_area,
      &                        cell_area_for_lonlat_outputs,longitude_deg
@@ -599 +600 @@
 ! GCM. Read netcdf initial physical parameters.
          CALL phyetat0 ("startfi.nc",0,0,
-     &         nsoilmx,ngrid,nlayer,nq,
+     &         nsoilmx,ngrid,nlayer,nq,nqsoil,
      &         day_ini,time_phys,
      &         tsurf,tsoil,albedo,emis,
-     &         q2,qsurf,tauscaling,totcloudfrac,wstar,
+     &         q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar,
      &         watercap,perenial_co2ice,
      &         def_slope,def_slope_mean,subslope_dist)
@@ -2621 +2622 @@
      .          icount,' date=',ztime_fin
             
-          call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
+          call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,
      .                ptimestep,ztime_fin,
      .                tsurf,tsoil,inertiesoil,albedo,
-     .                emis,q2,qsurf,tauscaling,totcloudfrac,wstar,
-     .                watercap,perenial_co2ice)
+     .                emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,
+     .                wstar,watercap,perenial_co2ice)
           ENDIF ! of IF (write_restart)
 

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

@@ -1 @@
-      subroutine soil_settings(nid,ngrid,nsoil,tsurf,tsoil,indextime)
+      subroutine soil_settings(nid,ngrid,nsoil,nqsoil,tsurf,tsoil,
+     &                         qsoil,indextime)
 
 !      use netcdf
       use comsoil_h, only: layer, mlayer, inertiedat, inertiesoil,
-     &                     volcapa,flux_geo
+     &                     volcapa,flux_geo,adsorption_soil,
+     &                     igcm_h2o_vap_soil,igcm_h2o_ice_soil,
+     &                     igcm_h2o_vap_ads
       use iostart, only: inquire_field_ndims, get_var, get_field,
      &                   inquire_field, inquire_dimension_length
@@ -43 +46 @@
       integer,intent(in) :: ngrid       ! # of horizontal grid points
       integer,intent(in) :: nsoil       ! # of soil layers
-      real,intent(in) :: tsurf(ngrid,nslope)   ! surface temperature
+      integer,intent(in) :: nqsoil      ! # of tracers in the soil
+      real,intent(in) :: tsurf(ngrid,nslope)   ! surface temperature [K]
       integer,intent(in) :: indextime   ! position on time axis
 !  output:
-      real,intent(out) :: tsoil(ngrid,nsoil,nslope)     ! soil temperature
+      real,intent(out) :: tsoil(ngrid,nsoil,nslope)     ! soil temperature [K]
+      real,intent(out) :: qsoil(ngrid,nsoil,nqsoil,nslope) ! Tracers in the subsurface [kg/kg]
 
 !======================================================================
@@ -81 +86 @@
 
       logical :: found,ok
+
+      character (len=30):: txt 
       
 !======================================================================
@@ -158 +165 @@
         layer(iloop)=lay1*(alpha**(iloop-1))
       enddo
-
 
 
@@ -439 +445 @@
 
 
+! 7. Adsorption
+! ----------------------------------------------------------------------
+
+
+      if (adsorption_soil) then
+! Subsurface water vapor
+             txt="h2o_vap_soil"
+             write(*,*) 'phyetat0: loading subsurface tracer', 
+     &                               ' h2o_vap_soil'
+             call get_field(txt,qsoil(:,:,igcm_h2o_vap_soil,:),found,
+     &                                                 indextime)
+             write(*,*) 'found',found
+             write(*,*) 'igcm_',igcm_h2o_vap_soil
+             write(*,*) 'q=',qsoil(:,:,:,:)
+             if (.not.found) then
+               write(*,*) "phyetat0: Failed loading <",trim(txt),">"
+               write(*,*) "         ",trim(txt)," is set to zero"
+               qsoil(:,:,igcm_h2o_vap_soil,:)= 0. 
+             else
+               write(*,*) "phyetat0: suburface tracer <",trim(txt),
+     &              "> range:", minval(qsoil(:,:,igcm_h2o_vap_soil,:)), 
+     &                          maxval(qsoil(:,:,igcm_h2o_vap_soil,:))
+             endif
+! Subsurface ice
+               txt="h2o_ice_soil"
+               write(*,*) 'phyetat0: loading subsurface tracer', 
+     &                                ' h2o_ice_soil'
+             call get_field(txt,qsoil(:,:,igcm_h2o_ice_soil,:),found,
+     &                                                  indextime)
+             if (.not.found) then
+               write(*,*) "phyetat0: Failed loading <",trim(txt),">"
+               write(*,*) "         ",trim(txt)," is set to zero"
+               qsoil(:,:,igcm_h2o_ice_soil,:)= 0. 
+             else
+               write(*,*) "phyetat0: suburface tracer <",trim(txt),
+     &                     "> range:", 
+     &                     minval(qsoil(:,:,igcm_h2o_ice_soil,:)),
+     &                     maxval(qsoil(:,:,igcm_h2o_ice_soil,:))
+             endif
+! Adsorbed water
+             txt="h2o_vap_ads"
+               write(*,*) 'phyetat0: loading subsurface tracer', 
+     &                                ' h2o_vap_ads'
+             call get_field(txt,qsoil(:,:,igcm_h2o_vap_ads,:),found,
+     &                                                indextime)
+             if (.not.found) then
+               write(*,*) "phyetat0: Failed loading <",trim(txt),">"
+               write(*,*) "         ",trim(txt)," is set to zero"
+               qsoil(:,:,igcm_h2o_vap_ads,:)= 0. 
+             else
+               write(*,*) "phyetat0: suburface tracer <",trim(txt),">
+     &                     range:", 
+     &                     minval(qsoil(:,:,igcm_h2o_vap_ads,:)), 
+     &                     maxval(qsoil(:,:,igcm_h2o_vap_ads,:))
+             endif
+
+      endif ! of adsorption_soil
+
+
       
-! 7. Report min and max values of soil temperatures and thermal inertias
+! 8. Report min and max values of soil temperatures and thermal inertias
 ! ----------------------------------------------------------------------