Changeset 2887

Timestamp:

Feb 1, 2023, 9:34:37 AM (22 months ago)

Author:

llange

Message:

Mars PCM
Include geothermal flux as boundary conditions when solving the
diffusion equation for the soil temperature
By default, the geothermal flux is set to 0. To have geothermal flux >
0., it must be specified in the startfi. (I suggest to use a value of
30e-3 in this case)
Will be use with the PEM for studies of "paleosubsurfaces"
LL

Location:

trunk/LMDZ.MARS/libf/phymars

Files:

• comsoil_h.F90 (modified) (3 diffs)
• phyetat0_mod.F90 (modified) (2 diffs)
• phyredem.F90 (modified) (2 diffs)
• physiq_mod.F (modified) (3 diffs)
• soil.F (modified) (2 diffs)
• soil_settings.F (modified) (2 diffs)

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

@@ -25 +25 @@
   real,save,allocatable :: beta(:,:)        ! beta_k coefficients
   real,save :: mu
+  real,save,allocatable :: flux_geo(:)       ! Geothermal Flux (W/m^2)
 
-!$OMP THREADPRIVATE(tsoil,mthermdiff,thermdiff,coefq,coefd,alph,beta,mu )
+!$OMP THREADPRIVATE(tsoil,mthermdiff,thermdiff,coefq,coefd,alph,beta,mu,flux_geo)
 
 contains
@@ -47 +48 @@
     allocate(alph(ngrid,nsoilmx-1))
     allocate(beta(ngrid,nsoilmx-1))
-
+    allocate(flux_geo(ngrid))
  
   end subroutine ini_comsoil_h
@@ -66 +67 @@
     if (allocated(alph)) deallocate(alph)
     if (allocated(beta)) deallocate(beta)
-
+    if (allocated(flux_geo)) deallocate(flux_geo)
   end subroutine end_comsoil_h
 

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

@@ -24 +24 @@
   use dust_param_mod, only: dustscaling_mode
   USE ioipsl_getin_p_mod, ONLY : getin_p
-
+  use comsoil_h, only: flux_geo
   implicit none
   
@@ -644 +644 @@
   ! as well as thermal inertia and volumetric heat capacity
   call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
+else
+    flux_geo(:) = 0.
 endif ! of if (startphy_file)
 

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

@@ -162 +162 @@
   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
   implicit none
   
@@ -314 +314 @@
   endif
 
+
+  ! Geothermal Flux
+     call put_field('flux_geo','Geothermal flux',flux_geo,time)
   ! Close file
   call close_restartphy

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

@@ -3598 +3598 @@
 
          ! Write soil temperature
-!        call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
-!    &                     3,tsoil)
+        call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
+     &                     3,tsoil)
          ! Write surface temperature
 !        call writediagsoil(ngrid,"tsurf","Surface temperature","K",
@@ -3628 +3628 @@
 c        CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
 c        CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
+        call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
+     &                     3,tsoil)
 
 ! THERMALS STUFF 1D
@@ -3977 +3979 @@
      &                   log(zplay(1,nlayer)/zplay(1,nlayer-1))*
      &                   rnew(1,nlayer)*tmean/g
+
+
 #endif
 

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

@@ -6 +6 @@
       use comsoil_h, only: layer, mlayer, volcapa,
      &                     mthermdiff, thermdiff, coefq,
-     &                     coefd, alph, beta, mu
+     &                     coefd, alph, beta, mu,flux_geo
       use surfdat_h, only: watercaptag, inert_h2o_ice
 
@@ -141 +141 @@
       do ig=1,ngrid
         beta(ig,nsoil-1)=coefq(nsoil-1)*tsoil(ig,nsoil)
-     &                   /(coefq(nsoil-1)+coefd(ig,nsoil-1))
+     &            /(coefq(nsoil-1)+coefd(ig,nsoil-1))
+     &            +flux_geo(ig)/(coefq(nsoil-1) + coefd(ig,nsoil-1))
       enddo
+   
 ! Other layers
       do ik=nsoil-2,1,-1

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

@@ -2 +2 @@
 
 !      use netcdf
-      use comsoil_h, only: layer, mlayer, inertiedat, volcapa
+      use comsoil_h, only: layer, mlayer, inertiedat, volcapa,flux_geo
       use iostart, only: inquire_field_ndims, get_var, get_field,
      &                   inquire_field, inquire_dimension_length
@@ -435 +435 @@
         deallocate(oldtsoil)
       endif ! of if (interpol)
-      
-! 6. Report min and max values of soil temperatures and thermal inertias
+
+
+      
+! 6. Load Geothermal Flux
+! ----------------------------------------------------------------------
+
+
+      call get_field("flux_geo",flux_geo,found,timeindex=indextime)
+      if (.not.found) then
+           write(*,*) "soil_settings: Failed loading <flux_geo>"
+           write(*,*) "soil_settings: Will put  <flux_geo> to 0."
+           flux_geo(:) = 0.
+      endif
+
+
+      
+! 7. Report min and max values of soil temperatures and thermal inertias
 ! ----------------------------------------------------------------------