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

Timestamp:

Nov 20, 2024, 3:53:19 PM (3 days 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_fast_subs_mars.F90 (modified) (10 diffs)

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

@@ -1 @@
-module thermalmodelparam_mars
+MODULE fast_subs_mars
+
+implicit none
+
   ! parameters for thermal model
   ! they are only used in the subroutines below
@@ -8 +11 @@
   real(8), parameter :: emiss0 = 1.  ! emissivity of dry surface
   integer, parameter :: EQUILTIME =15 ! [Mars years]
-end module thermalmodelparam_mars
-
-
+  
+  integer, parameter :: NMAX = 1000
+
+CONTAINS
 !*****************************************************
 ! Subroutines for fast method
@@ -25 +29 @@
 ! latitude  [degree]
 !*************************************************************************
-  use miscparameters, only : d2r, NMAX, icedensity
-  use allinterfaces, except_this_one => icelayer_mars
+  use ice_table_mod, only: rho_ice
+  use fast_subs_univ, only: icechanges
   !use omp_lib
   implicit none
@@ -53 +57 @@
      Diff = 4e-4*600./p0(k)
      fracIR = 0.04*p0(k)/600.; fracDust = 0.02*p0(k)/600.
-     B = Diff*bigstep*86400.*365.24/(porosity*icedensity)
+     B = Diff*bigstep*86400.*365.24/(porosity*927.)
+     !B = Diff*bigstep*86400.*365.24/(porosity*rho_ice(Tb(),'h2o'))
      
      typeT = -9
@@ -140 +145 @@
 !  Tb>0 initializes temperatures with Tb 
 !***********************************************************************
-  use miscparameters
-  use thermalmodelparam_mars
-  use allinterfaces, except_this_one => ajsub_mars
+  use fast_subs_univ, only: depths_avmeth, equildepth
+  use constants_marspem_mod, only: sols_per_my, sec_per_sol
   implicit none
   integer, intent(IN) :: nz, typeT
@@ -166 +170 @@
   real(8) Tsurf, Tco2frost, albedo, Fsurf, m, dE, emiss, T(NMAX)
   real(8) Told(nz), Fsurfold, Tsurfold, Tmean0, avrho2
-  real(8) rhosatav0, rhosatav(nz), rlow 
-  real(8), external :: psv, tfrostco2
+  real(8) rhosatav0, rhosatav(nz), rlow
   
-  tmax = EQUILTIME*solsperyear
+  tmax = EQUILTIME*sols_per_my
   nsteps = int(tmax/dt)     ! calculate total number of timesteps
 
@@ -207 +210 @@
   do n=0,nsteps-1
      time = (n+1)*dt         !   time at n+1 
-     tdays = time*(marsDay/earthDay) ! parenthesis may improve roundoff
+   !  tdays = time*(sec_per_sol/earthDay) ! parenthesis may improve roundoff
    !  call generalorbit(tdays,a,ecc,omega,eps,marsLs,marsDec,marsR)
    !  HA = 2.*pi*mod(time,1.d0)  ! hour angle
@@ -217 +220 @@
      Told(1:nz) = T(1:nz)
      if (m<=0. .or. Tsurf>Tco2frost) then
-       ! call conductionQ(nz,z,dt*marsDay,Qn,Qnp1,T,ti,rhocv,emiss, &
+       ! call conductionQ(nz,z,dt*sec_per_sol,Qn,Qnp1,T,ti,rhocv,emiss, &
        !      &           Tsurf,Fgeotherm,Fsurf)
      endif
      if (Tsurf<Tco2frost .or. m>0.) then ! CO2 condensation
         T(1:nz) = Told(1:nz)
-        !call conductionT(nz,z,dt*marsDay,T,Tsurfold,Tco2frost,ti, &
+        !call conductionT(nz,z,dt*sec_per_sol,T,Tsurfold,Tco2frost,ti, &
              !&              rhocv,Fgeotherm,Fsurf) 
         Tsurf = Tco2frost
     !    dE = (- Qn - Qnp1 + Fsurfold + Fsurf + &
     !         &      emiss*sigSB*(Tsurfold**4+Tsurf**4)/2.
-        m = m + dt*marsDay*dE/Lco2frost
+        m = m + dt*sec_per_sol*dE/Lco2frost
      endif
      if (Tsurf>Tco2frost .or. m<=0.) then
@@ -238 +241 @@
      !Qn=Qnp1
      
-     if (time>=tmax-solsperyear) then
+     if (time>=tmax-sols_per_my) then
         Tmean1 = Tmean1 + Tsurf
         Tmean3 = Tmean3 + T(nz)
@@ -301 +304 @@
 tfrostco2 = 3182.48/(23.3494+log(100./p))
 end function
+
+!=======================================================================
+
+      real*8 function psv(T)
+!     saturation vapor pressure of H2O ice [Pascal]
+!     input is temperature [Kelvin]
+      implicit none
+      real*8 T
+
+!-----parametrization 1
+!      real*8 DHmelt,DHvap,DHsub,R,pt,Tt,C
+!      parameter (DHmelt=6008.,DHvap=45050.)
+!      parameter (DHsub=DHmelt+DHvap) ! sublimation enthalpy [J/mol]
+!      parameter (R=8.314,pt=6.11e2,Tt=273.16)
+!      C = (DHsub/R)*(1./T - 1./Tt)
+!      psv = pt*exp(-C)
+
+!-----parametrization 2
+!     eq. (2) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
+!     differs from parametrization 1 by only 0.1%
+      real*8 A,B
+      parameter (A=-6143.7, B=28.9074)
+      psv = exp(A/T+B)  ! Clapeyron
+
+!-----parametrization 3      
+!     eq. (7) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
+!     psv = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T)
+      
+      end function psv
+
+END MODULE fast_subs_mars