Changeset 3527 for trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90

Timestamp:

Nov 20, 2024, 3:53:19 PM (25 hours ago)

Author:

jbclement

Message:

PEM:

• Removing redundant Norbert Schorghofer's subroutines/parameters (follow-up of r3526);
• Making all Norbert Schorghofer's subroutines with modern explicit interface via modules;
• Cleaning of "glaciers_mod.F90";
• Optimization for the computation of ice density according to temperature by using a function.

JBC

File:

• trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90 (modified) (4 diffs)

trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90

@@ +1 @@
+MODULE dyn_ss_ice_m_mod
+
+implicit none
+
+CONTAINS
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -15 +21 @@
 ! orbital elements remain constant
 !***********************************************************************
-  use miscparameters, only : pi, d2r, NMAX, marsDay, solsperyear 
-  !use allinterfaces
+  use constants_marspem_mod, only: sec_per_sol
+  use fast_subs_mars, only: psv, icelayer_mars, NMAX
+#ifndef CPP_STD
+    use comcstfi_h,   only: pi
+#else
+    use comcstfi_mod, only: pi
+#endif
   implicit none
   integer, parameter :: NP=1   ! # of sites
@@ -105 +116 @@
  !    print *,'  latitude (deg)',latitude(k),' rho*c (J/m^3/K)',rhoc(k),' thIn=',thIn(k)
  !    print *,'  total pressure=',p0(k),'partial pressure=',pfrost(k)
-     delta = thIn(k)/rhoc(k)*sqrt(marsDay/pi)
+     delta = thIn(k)/rhoc(k)*sqrt(sec_per_sol/pi)
  !    print *,'  skin depths (m)',delta,delta*sqrt(solsperyear)
      call soilthprop(porosity,1.d0,rhoc(k),thIn(k),1,newrhoc,newti,icefrac)
@@ -170 +181 @@
   
 end subroutine dyn_ss_ice_m
+
+!=======================================================================
+
+subroutine soilthprop(porosity,fill,rhocobs,tiobs,layertype, &
+     &     newrhoc,newti,icefrac)
+!***********************************************************************
+! soilthprop: assign thermal properties of icy soil or dirty ice
+!
+!     porositiy = void space / total volume
+!     rhof = density of free ice in space not occupied by regolith [kg/m^3]
+!     fill = rhof/icedensity <=1 (only relevant for layertype 1)
+!     rhocobs = heat capacity per volume of dry regolith [J/m^3]
+!     tiobs = thermal inertia of dry regolith [SI-units]
+!     layertype: 1=interstitial ice, 2=pure ice or ice with dirt
+!                3=pure ice, 4=ice-cemented soil, 5=custom values
+!     icefrac: fraction of ice in icelayer (only relevant for layertype 2)
+!     output are newti and newrhoc
+!***********************************************************************
+  implicit none
+  integer, intent(IN) :: layertype
+  real(8), intent(IN) :: porosity, fill, rhocobs, tiobs
+  real(8), intent(OUT) :: newti, newrhoc
+  real(8), intent(IN) :: icefrac
+  real(8) kobs, cice, icedensity, kice
+  !parameter (cice=2000.d0, icedensity=926.d0, kice=2.4d0) ! unaffected by scaling
+  parameter (cice=1540.d0, icedensity=927.d0, kice=3.2d0) ! at 198 Kelvin
+  real(8) fA, ki0, ki, k
+  real(8), parameter :: kw=3. ! Mellon et al., JGR 102, 19357 (1997)
+
+  kobs = tiobs**2/rhocobs
+  ! k, rhoc, and ti are defined in between grid points
+  ! rhof and T are defined on grid points
+
+  newrhoc = -9999.
+  newti  = -9999.
+
+  select case (layertype)
+  case (1) ! interstitial ice
+     newrhoc = rhocobs + porosity*fill*icedensity*cice
+     if (fill>0.) then
+        !--linear addition (option A)
+        k = porosity*fill*kice + kobs
+        !--Mellon et al. 1997 (option B)
+        ki0 = porosity/(1/kobs-(1-porosity)/kw)
+        fA = sqrt(fill)
+        ki = (1-fA)*ki0 + fA*kice
+        !k = kw*ki/((1-porosity)*ki+porosity*kw)
+     else
+        k = kobs
+     endif
+     newti = sqrt(newrhoc*k)
+     
+  case (2)  ! massive ice (pure or dirty ice)
+     newrhoc = rhocobs*(1.-icefrac)/(1.-porosity) + icefrac*icedensity*cice
+     k = icefrac*kice + (1.-icefrac)*kw
+     newti = sqrt(newrhoc*k)
+  
+  case (3)  ! all ice, special case of layertype 2, which doesn't use porosity
+     newrhoc = icedensity*cice
+     k = kice 
+     newti = sqrt(newrhoc*k)
+  
+  case (4)  ! pores completely filled with ice, special case of layertype 1
+     newrhoc = rhocobs + porosity*icedensity*cice
+     k = porosity*kice + kobs ! option A, end-member case of type 1, option A 
+     !k = kw*kice/((1-porosity)*kice+porosity*kw) ! option B, harmonic average
+     newti = sqrt(newrhoc*k)
+
+  case (5)  ! custom values
+     ! values from Mellon et al. (2004) for ice-cemented soil
+     newrhoc = 2018.*1040.
+     k = 2.5
+     newti = sqrt(newrhoc*k)
+
+  case default
+     error stop 'invalid layer type'
+     
+  end select
+  
+end subroutine soilthprop
+
+
+!=======================================================================
+      
+      real*8 function frostpoint(p)
+!     inverse of psv
+!     input is partial pressure [Pascal]
+!     output is temperature [Kelvin]
+      implicit none
+      real*8 p
+      
+!-----inverse of parametrization 1
+!      real*8 DHmelt,DHvap,DHsub,R,pt,Tt
+!      parameter (DHmelt=6008.,DHvap=45050.)
+!      parameter (DHsub=DHmelt+DHvap)
+!      parameter (R=8.314,pt=6.11e2,Tt=273.16)
+!      frostpoint = 1./(1./Tt-R/DHsub*log(p/pt))
+      
+!-----inverse of parametrization 2
+!     inverse of eq. (2) in Murphy & Koop (2005)
+      real*8 A,B
+      parameter (A=-6143.7, B=28.9074)
+      frostpoint = A / (log(p) - B)
+
+!-----approximate inverse of parametrization 3
+!     eq. (8) in Murphy & Koop (2005)
+!      frostpoint = (1.814625*log(p) + 6190.134)/(29.120 - log(p))
+      
+      end function frostpoint
+
+END MODULE dyn_ss_ice_m_mod