Changeset 3493 for trunk


Ignore:
Timestamp:
Nov 5, 2024, 5:19:11 PM (7 weeks ago)
Author:
jbclement
Message:

PEM:

  • Renaming of Norbert Schorghofer's subroutines with the prefix 'NS_';
  • Making the extension of all NS's subroutines as '.F90';
  • Deletion of the wrapper subroutine;
  • Making the initialization, variables management and arguments of the main subroutine for the dynamic computation of ice table to be more suitable.

JBC

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
1 deleted
5 edited
15 moved

Legend:

Unmodified
Added
Removed

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

    r3492 r3493  
    1212
    1313
    14 SUBROUTINE dyn_ss_ice_m(ssi_depth_in,T1,T_in,nsoil,thIn,p0,pfrost,porefill_in,porefill,ssi_depth)
     14SUBROUTINE dyn_ss_ice_m(ssi_depth_in,T1,Tb,nz,thIn,p0,pfrost,porefill_in,porefill,ssi_depth)
    1515
    1616!***********************************************************************
    … …  
    2222  implicit none
    2323  integer, parameter :: NP=1   ! # of sites
    24   integer nz, i, k, nsoil, iloop
     24  integer nz, i, k, iloop
    2525  real(8)  zmax, delta, z(NMAX), icetime, porosity, icefrac
    2626  real(8), dimension(NP) ::  albedo, thIn, rhoc
    … …  
    3030  real(8), dimension(NP) :: zdepthF, zdepthE, zdepthT, zdepthG
    3131  real(8), dimension(NMAX,NP) :: porefill, porefill_in
    32   real(8), dimension(NP) :: Tb, T_in, Tmean1, Tmean3, avrho1
     32  real(8), dimension(nz) :: Tb
     33  real(8), dimension(NP) :: Tmean1, Tmean3, avrho1
    3334  real(8) tmax, tlast, avrho1prescribed(NP), l1
    3435  real(8), external :: smartzfac
    … …  
    4243
    4344  ! parameters that never ever change
    44   nz=nsoil
    4545  porosity = 0.4d0  ! porosity of till
    4646  !rhoc(:) = 1500.*800.  ! will be overwritten
    … …  
    6969  !call setgrid(nz,z,zmax,zfac)
    7070  l1=2.e-4
    71   do iloop=0,nsoil
     71  do iloop=0,nz
    7272    z(iloop) = l1*(1+iloop**2.9*(1-exp(-real(iloop)/20.)))
    7373  enddo
    … …  
    8989 
    9090  ! initializations
    91   Tb=T_in
     91  !Tb = -9999.
    9292  zdepthF(:) = -9999.
    9393

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

    r3492 r3493  
    3232  real(8), intent(IN) :: bigstep
    3333  real(8), intent(IN) :: thIn(NP), rhoc(NP), z(NMAX), porosity, pfrost(NP)
    34   real(8), intent(INOUT) :: Tb(NP), porefill(nz,NP), zdepthF(NP), zdepthT(NP)
     34  real(8), intent(INOUT) :: porefill(nz,NP), zdepthF(NP), zdepthT(NP)
     35  real(8), intent(INOUT) :: Tb(nz)
    3536  real(8), intent(OUT), dimension(NP) :: zdepthE, Tmean1, Tmean3, zdepthG
    3637  real(8), intent(IN), dimension(NP) ::  albedo, p0
    … …  
    149150  real(8), intent(OUT) :: avdrho, avdrhoP  ! difference in annual mean vapor density
    150151  real(8), intent(OUT) :: avrho1  ! mean vapor density on surface
    151   real(8), intent(INOUT) :: Tb, Tmean1
     152  real(8), intent(INOUT) :: Tmean1
     153  real(8), intent(INOUT) :: Tb(nz)
    152154  integer, intent(OUT) :: typeF  ! index of depth below which filling occurs
    153155  real(8), intent(OUT) :: zdepthE, zdepthF

  • trunk/LMDZ.COMMON/libf/evolution/NS_generalorbit.F90

    r3492 r3493  
    1 C=============================================================
    2 C Subroutine to return distance, longitude, and declination of
    3 C the sun in planetocentric coordinates from orbital elements
    4 C
    5 C INPUTS:
    6 C edays = time since perihelion (earth days)
    7 C a = semimajor axis (AU)
    8 C ecc = eccentricity
    9 C omega = Ls of perihelion, relative to equinox (radians)
    10 C eps = obliquity (radians)
    11 C
    12 C OUTPUTS:
    13 C Ls = areocentric longitude (radians)
    14 C dec = planetocentric solar declination (radians)
    15 C r = heliocentric distance (AU)
    16 C=============================================================
     1!=============================================================
     2! Subroutine to return distance, longitude, and declination of
     3! the sun in planetocentric coordinates from orbital elements
     4!
     5! INPUTS:
     6! edays = time since perihelion (earth days)
     7! a = semimajor axis (AU)
     8! ecc = eccentricity
     9! omega = Ls of perihelion, relative to equinox (radians)
     10! eps = obliquity (radians)
     11!
     12! OUTPUTS:
     13! Ls = areocentric longitude (radians)
     14! dec = planetocentric solar declination (radians)
     15! r = heliocentric distance (AU)
     16!=============================================================
    1717
    1818      SUBROUTINE generalorbit(edays,a,ecc,omega,eps,Ls,dec,r)
    … …  
    2525      real*8 M,E,nu,T,Eold
    2626
    27 c     T = orbital period (days)
     27!     T = orbital period (days)
    2828      T = sqrt(4*pi**2/(6.674e-11*1.989e30)*(a*149.598e9)**3)/86400.
    2929
    30 c     M = mean anomaly (radians)
     30!     M = mean anomaly (radians)
    3131      M = 2.*pi*edays/T  ! M=0 at perihelion
    3232
    33 c     E = eccentric anomaly
    34 c     solve M = E - ecc*sin(E) by newton method
     33!     E = eccentric anomaly
     34!     solve M = E - ecc*sin(E) by newton method
    3535      E = M
    3636      do j=1,10
    … …  
    4040      enddo
    4141
    42 c     nu = true anomaly
     42!     nu = true anomaly
    4343      !nu = acos(cos(E)-ecc/(1.-ecc*cos(E)))
    4444      !nu = sqrt(1-ecc^2)*sin(E)/(1.-ecc*cos(E))

  • trunk/LMDZ.COMMON/libf/evolution/NS_psv.F90

    r3492 r3493  
    11      real*8 function psv(T)
    2 C     saturation vapor pressure of H2O ice [Pascal]
    3 C     input is temperature [Kelvin]
     2!     saturation vapor pressure of H2O ice [Pascal]
     3!     input is temperature [Kelvin]
    44      implicit none
    55      real*8 T
    66
    7 C-----parametrization 1
    8 c      real*8 DHmelt,DHvap,DHsub,R,pt,Tt,C
    9 c      parameter (DHmelt=6008.,DHvap=45050.)
    10 c      parameter (DHsub=DHmelt+DHvap) ! sublimation enthalpy [J/mol]
    11 c      parameter (R=8.314,pt=6.11e2,Tt=273.16)
    12 c      C = (DHsub/R)*(1./T - 1./Tt)
    13 c      psv = pt*exp(-C)
     7!-----parametrization 1
     8!      real*8 DHmelt,DHvap,DHsub,R,pt,Tt,C
     9!      parameter (DHmelt=6008.,DHvap=45050.)
     10!      parameter (DHsub=DHmelt+DHvap) ! sublimation enthalpy [J/mol]
     11!      parameter (R=8.314,pt=6.11e2,Tt=273.16)
     12!      C = (DHsub/R)*(1./T - 1./Tt)
     13!      psv = pt*exp(-C)
    1414
    15 C-----parametrization 2
    16 C     eq. (2) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
    17 C     differs from parametrization 1 by only 0.1%
     15!-----parametrization 2
     16!     eq. (2) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
     17!     differs from parametrization 1 by only 0.1%
    1818      real*8 A,B
    1919      parameter (A=-6143.7, B=28.9074)
    2020      psv = exp(A/T+B)  ! Clapeyron
    2121
    22 C-----parametrization 3     
    23 C     eq. (7) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
    24 c     psv = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T)
     22!-----parametrization 3     
     23!     eq. (7) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
     24!     psv = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T)
    2525     
    2626      end
    … …  
    2929     
    3030      real*8 function frostpoint(p)
    31 C     inverse of psv
    32 C     input is partial pressure [Pascal]
    33 C     output is temperature [Kelvin]
     31!     inverse of psv
     32!     input is partial pressure [Pascal]
     33!     output is temperature [Kelvin]
    3434      implicit none
    3535      real*8 p
    3636     
    37 C-----inverse of parametrization 1
    38 c      real*8 DHmelt,DHvap,DHsub,R,pt,Tt
    39 c      parameter (DHmelt=6008.,DHvap=45050.)
    40 c      parameter (DHsub=DHmelt+DHvap)
    41 c      parameter (R=8.314,pt=6.11e2,Tt=273.16)
    42 c      frostpoint = 1./(1./Tt-R/DHsub*log(p/pt))
     37!-----inverse of parametrization 1
     38!      real*8 DHmelt,DHvap,DHsub,R,pt,Tt
     39!      parameter (DHmelt=6008.,DHvap=45050.)
     40!      parameter (DHsub=DHmelt+DHvap)
     41!      parameter (R=8.314,pt=6.11e2,Tt=273.16)
     42!      frostpoint = 1./(1./Tt-R/DHsub*log(p/pt))
    4343     
    44 C-----inverse of parametrization 2
    45 C     inverse of eq. (2) in Murphy & Koop (2005)
     44!-----inverse of parametrization 2
     45!     inverse of eq. (2) in Murphy & Koop (2005)
    4646      real*8 A,B
    4747      parameter (A=-6143.7, B=28.9074)
    4848      frostpoint = A / (log(p) - B)
    4949
    50 C-----approximate inverse of parametrization 3
    51 C     eq. (8) in Murphy & Koop (2005)
    52 c      frostpoint = (1.814625*log(p) + 6190.134)/(29.120 - log(p))
     50!-----approximate inverse of parametrization 3
     51!     eq. (8) in Murphy & Koop (2005)
     52!      frostpoint = (1.814625*log(p) + 6190.134)/(29.120 - log(p))
    5353     
    5454      end

  • trunk/LMDZ.COMMON/libf/evolution/NS_psvco2.F90

    r3492 r3493  
    11      real*8 function psvco2(T)
    2 C     solid-vapor transition for CO2
    3 C     returns saturation pressure in Pascal, input is temperature in Kelvin
     2!     solid-vapor transition for CO2
     3!     returns saturation pressure in Pascal, input is temperature in Kelvin
    44      implicit none
    55      real*8 T
    … …  
    1010
    1111      real*8 function tfrostco2(p)
    12 C     the inverse of function psvco2
    13 C     input is pressure in Pascal, output is temperature in Kelvin
     12!     the inverse of function psvco2
     13!     input is pressure in Pascal, output is temperature in Kelvin
    1414      implicit none
    1515      real*8 p
    … …  
    2020
    2121
    22 C------------------------------------------------------------------------------
    23 C     Antoine equation parameters from NIST, 154K-196K
    24 C     based on Giauque and Egan (1937)
    25 C     A=6.81228, B=1301.679, C=-3.494
    26 c     p = 10**(A-(B/(T+C)))*1.e5
     22!------------------------------------------------------------------------------
     23!     Antoine equation parameters from NIST, 154K-196K
     24!     based on Giauque and Egan (1937)
     25!     A=6.81228, B=1301.679, C=-3.494
     26!     p = 10**(A-(B/(T+C)))*1.e5
    2727
    28 C     Expressions from Int. Crit. Tabl. 3.207, based on many references
    29 C     mm2Pa = 133.32
    30 C     -135C to -56.7C (138K-216K)
    31 c     p = 10**(-0.05223*26179.3/T + 9.9082)*mm2Pa
    32 C     -183C to -135C (90K-138K)
    33 c     p = 10**(-1275.62/T + 0.006833*T + 8.3071)*mm2Pa
    34 C     Expressions from Int. Crit. Tabl. 3.208, based on Henning
    35 C     -110 to -80C (163K-193K)
    36 c     p = 10**(- 1279.11/T + 1.75*log10(T) - 0.0020757*T + 5.85242)*mm2Pa
    37 c     p = 10**(- 1352.6/T + 9.8318)*mm2Pa
     28!     Expressions from Int. Crit. Tabl. 3.207, based on many references
     29!     mm2Pa = 133.32
     30!     -135C to -56.7C (138K-216K)
     31!     p = 10**(-0.05223*26179.3/T + 9.9082)*mm2Pa
     32!     -183C to -135C (90K-138K)
     33!     p = 10**(-1275.62/T + 0.006833*T + 8.3071)*mm2Pa
     34!     Expressions from Int. Crit. Tabl. 3.208, based on Henning
     35!     -110 to -80C (163K-193K)
     36!     p = 10**(- 1279.11/T + 1.75*log10(T) - 0.0020757*T + 5.85242)*mm2Pa
     37!     p = 10**(- 1352.6/T + 9.8318)*mm2Pa
    3838     
    39 C     Mars book (1992), p959, fit by chapter authors
    40 c     p = exp(23.3494 - 3182.48/T)*100.   ! 120-160K
    41 c     p = exp(25.2194 - 3311.57/T - 6.71e-3*T)*100  ! 100-170K
    42 C     Mars book (1992), p960, based on Miller & Smythe (1970)
    43 c     p = exp(26.1228 - 3385.26/T - 9.4461e-3*T)*100 ! 123-173K
     39!     Mars book (1992), p959, fit by chapter authors
     40!     p = exp(23.3494 - 3182.48/T)*100.   ! 120-160K
     41!     p = exp(25.2194 - 3311.57/T - 6.71e-3*T)*100  ! 100-170K
     42!     Mars book (1992), p960, based on Miller & Smythe (1970)
     43!     p = exp(26.1228 - 3385.26/T - 9.4461e-3*T)*100 ! 123-173K
    4444
    45 C     Fray & Schmitt, PSS 57, 2053 (2009)
    46 C     A0=1.476e1, A1=-2.571e3, A2=-7.781e4, A3=4.325e6, A4=-1.207e8, A5=1.350e9
    47 c     p = exp(A0+A1/T+A2/T**2+A3/T**3+A4/T**4+A5/T**5)*1e5 ! 40K-194.7K
    48 c     A0=1.861e1, A1=-4.154e3, A2=1.041e5
    49 c     p = exp(A0 + A1/T + A2/T**2)*1e5  ! 194.7K-216.58K
    50 C------------------------------------------------------------------------------
     45!     Fray & Schmitt, PSS 57, 2053 (2009)
     46!     A0=1.476e1, A1=-2.571e3, A2=-7.781e4, A3=4.325e6, A4=-1.207e8, A5=1.350e9
     47!     p = exp(A0+A1/T+A2/T**2+A3/T**3+A4/T**4+A5/T**5)*1e5 ! 40K-194.7K
     48!     A0=1.861e1, A1=-4.154e3, A2=1.041e5
     49!     p = exp(A0 + A1/T + A2/T**2)*1e5  ! 194.7K-216.58K
     50!------------------------------------------------------------------------------

  • trunk/LMDZ.COMMON/libf/evolution/NS_tridag.F90

    r3492 r3493  
    1 C================================================
    2 C Tridiagonal solver
    3 C================================================
     1!================================================
     2! Tridiagonal solver
     3!================================================
    44      SUBROUTINE tridag(a,b,c,r,u,n)
    55      INTEGER n,NMAX
    … …  
    1111         stop 'tridag: rewrite equations'
    1212      endif
    13 c      if(n.gt.NMAX) then
    14 c         print *, 'tridag: too many points, set NMAX>',n
    15 c         stop
    16 c      endif
     13!      if(n.gt.NMAX) then
     14!         print *, 'tridag: too many points, set NMAX>',n
     15!         stop
     16!      endif
    1717      bet=b(1)
    1818      u(1)=r(1)/bet
    19       do 11 j=2,n
     19      do j=2,n
    2020        gam(j)=c(j-1)/bet
    2121        bet=b(j)-a(j)*gam(j)
    22 c        if(bet.eq.0.)pause 'tridag failed'
     22!        if(bet.eq.0.)pause 'tridag failed'
    2323        u(j)=(r(j)-a(j)*u(j-1))/bet
    24 11    continue
    25       do 12 j=n-1,1,-1
     24      enddo
     25      do j=n-1,1,-1
    2626        u(j)=u(j)-gam(j+1)*u(j+1)
    27 12    continue
    28       return
     27      enddo
     28!      return
    2929      END
    30 C  (C) Copr. 1986-92 Numerical Recipes Software 0(9p#31&#5(+.
     30!  (C) Copr. 1986-92 Numerical Recipes Software 0(9p#31&#5(+.

  • trunk/LMDZ.COMMON/libf/evolution/changelog.txt

    r3492 r3493  
    460460== 05/11/2024 == JBC
    461461Correction in the formula to update the potential temperature introduced in r3442.
     462
     463== 05/11/2024 == JBC
     464- Renaming of Norbert Schorghofer's subroutines with the prefix 'NS_';
     465- Making the extension of all NS's subroutines as '.F90';
     466- Deletion of the wrapper subroutine;
     467- Making the initialization, variables management and arguments of the main subroutine for the dynamic computation of ice table to be more suitable.

  • trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90

    r3490 r3493  
    1010!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1111
    12 logical, save                           :: icetable_equilibrium     ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium
    13 logical, save                           :: icetable_dynamic         ! Boolean to say if the PEM needs to recompute the icetable depth with the dynamic method
    14 real, allocatable, dimension(:,:), save :: porefillingice_depth     ! ngrid x nslope: Depth of the ice table [m]
    15 real, allocatable, dimension(:,:), save :: porefillingice_thickness ! ngrid x nslope: Thickness of the ice table [m]
     12logical, save                           :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium
     13logical, save                           :: icetable_dynamic     ! Boolean to say if the PEM needs to recompute the icetable depth with the dynamic method
     14real, allocatable, dimension(:,:)       :: icetable_depth       ! ngrid x nslope: Depth of the ice table [m]
     15real, allocatable, dimension(:,:)       :: icetable_thickness   ! ngrid x nslope: Thickness of the ice table [m]
     16real, allocatable, dimension(:,:,:)     :: ice_porefilling      ! the amout of porefilling in each layer in each grid [m^3/m^3]
    1617
    1718!-----------------------------------------------------------------------
    1819contains
    1920!-----------------------------------------------------------------------
    20 SUBROUTINE ini_ice_table_porefilling(ngrid,nslope)
     21SUBROUTINE ini_ice_table(ngrid,nslope,nsoil)
    2122
    2223implicit none
    … …  
    2425integer, intent(in) :: ngrid  ! number of atmospheric columns
    2526integer, intent(in) :: nslope ! number of slope within a mesh
    26 
    27 allocate(porefillingice_depth(ngrid,nslope))
    28 allocate(porefillingice_thickness(ngrid,nslope))
    29 
    30 END SUBROUTINE ini_ice_table_porefilling
    31 
    32 !-----------------------------------------------------------------------
    33 SUBROUTINE end_ice_table_porefilling
    34 
    35 implicit none
    36 
    37 if (allocated(porefillingice_depth)) deallocate(porefillingice_depth)
    38 if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness)
    39 
    40 END SUBROUTINE end_ice_table_porefilling
     27integer, intent(in) :: nsoil  ! number of soil layers
     28
     29allocate(icetable_depth(ngrid,nslope))
     30if (icetable_equilibrium) then
     31    allocate(icetable_thickness(ngrid,nslope))
     32else if (icetable_dynamic) then
     33    allocate(ice_porefilling(ngrid,nsoil,nslope))
     34endif
     35
     36END SUBROUTINE ini_ice_table
     37
     38!-----------------------------------------------------------------------
     39SUBROUTINE end_ice_table
     40
     41implicit none
     42
     43if (allocated(icetable_depth)) deallocate(icetable_depth)
     44if (allocated(icetable_thickness)) deallocate(icetable_thickness)
     45if (allocated(ice_porefilling)) deallocate(ice_porefilling)
     46
     47END SUBROUTINE end_ice_table
    4148
    4249!------------------------------------------------------------------------------------------------------

  • trunk/LMDZ.COMMON/libf/evolution/pem.F90

    r3492 r3493  
    5454use orbit_param_criterion_mod,  only: orbit_param_criterion
    5555use recomp_orb_param_mod,       only: recomp_orb_param
    56 use ice_table_mod,              only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
    57                                       ini_ice_table_porefilling, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium,compute_massh2o_exchange_ssi
     56use ice_table_mod,              only: icetable_depth, icetable_thickness, end_ice_table, ice_porefilling, &
     57                                      ini_ice_table, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium, compute_massh2o_exchange_ssi
    5858use soil_thermalproperties_mod, only: update_soil_thermalproperties
    5959use time_phylmdz_mod,           only: daysec, dtphys
    … …  
    204204
    205205! Variables for surface and soil
    206 real, dimension(:,:),     allocatable :: tsurf_avg                          ! Physic x SLOPE field: Averaged Surface Temperature [K]
    207 real, dimension(:,:,:),   allocatable :: tsoil_avg                          ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K]
    208 real, dimension(:,:,:),   allocatable :: tsurf_PCM_timeseries               ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K]
    209 real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries          ! IG x SOIL x SLOPE x TIMES field: Non averaged Soil Temperature [K]
    210 real, dimension(:,:,:,:), allocatable :: tsoil_anom                         ! IG x SOIL x SLOPE x TIMES field: Amplitude between instataneous and yearly average soil temperature at the beginning [K]
    211 real, dimension(:,:),     allocatable :: tsurf_avg_yr1                      ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K]
    212 real, dimension(:,:),     allocatable :: Tsoil_locslope                     ! Physic x Soil: Intermediate when computing Tsoil [K]
    213 real, dimension(:),       allocatable :: Tsurf_locslope                     ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
    214 real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries       ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
    215 real, dimension(:,:),     allocatable :: watersurf_density_avg              ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
    216 real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries   ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
    217 real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg          ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3]
    218 real, dimension(:,:),     allocatable :: Tsurfavg_before_saved              ! Surface temperature saved from previous time step [K]
    219 real, dimension(:),       allocatable :: delta_co2_adsorbded                ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
    220 real, dimension(:),       allocatable :: delta_h2o_adsorbded                ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
    221 real                                  :: totmassco2_adsorbded               ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
    222 real                                  :: totmassh2o_adsorbded               ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
    223 logical                               :: bool_sublim                        ! logical to check if there is sublimation or not
    224 logical, dimension(:,:),  allocatable :: co2ice_disappeared                 ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
    225 real, dimension(:,:),     allocatable :: porefillingice_thickness_prev_iter ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
    226 real, dimension(:),       allocatable :: delta_h2o_icetablesublim           ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
    227 
    228 ! Dynamic ice table
    229 real, dimension(:,:),     allocatable :: ss_ice_pf                          ! the amout of porefilling in each layer in each grid [m^3/m^3]
    230 real, dimension(:,:),     allocatable :: ss_ice_pf_new                      ! the amout of porefilling in each layer in each grid after the ss module[m^3/m^3]
    231 real, dimension(:,:),     allocatable :: porefillingice_depth_new           ! new pure ice table depth
     206real, dimension(:,:),     allocatable :: tsurf_avg                        ! Physic x SLOPE field: Averaged Surface Temperature [K]
     207real, dimension(:,:,:),   allocatable :: tsoil_avg                        ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K]
     208real, dimension(:,:,:),   allocatable :: tsurf_PCM_timeseries             ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K]
     209real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries        ! IG x SOIL x SLOPE x TIMES field: Non averaged Soil Temperature [K]
     210real, dimension(:,:,:,:), allocatable :: tsoil_anom                       ! IG x SOIL x SLOPE x TIMES field: Amplitude between instataneous and yearly average soil temperature at the beginning [K]
     211real, dimension(:,:),     allocatable :: tsurf_avg_yr1                    ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K]
     212real, dimension(:,:),     allocatable :: Tsoil_locslope                   ! Physic x Soil: Intermediate when computing Tsoil [K]
     213real, dimension(:),       allocatable :: Tsurf_locslope                   ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
     214real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries     ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
     215real, dimension(:,:),     allocatable :: watersurf_density_avg            ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
     216real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
     217real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3]
     218real, dimension(:,:),     allocatable :: Tsurfavg_before_saved            ! Surface temperature saved from previous time step [K]
     219real, dimension(:),       allocatable :: delta_co2_adsorbded              ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
     220real, dimension(:),       allocatable :: delta_h2o_adsorbded              ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
     221real                                  :: totmassco2_adsorbded             ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
     222real                                  :: totmassh2o_adsorbded             ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
     223logical                               :: bool_sublim                      ! logical to check if there is sublimation or not
     224logical, dimension(:,:),  allocatable :: co2ice_disappeared               ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
     225real, dimension(:,:),     allocatable :: icetable_thickness_prev_iter     ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
     226real, dimension(:),       allocatable :: delta_h2o_icetablesublim         ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
     227real, dimension(:),       allocatable :: porefill                         ! Pore filling (output) to compute the dynamic ice table
     228real                                  :: ssi_depth                        ! Ice table depth (output) to compute the dynamic ice table
     229
    232230! Some variables for the PEM run
    233231real, parameter :: year_step = 1   ! Timestep for the pem
    … …  
    557555allocate(delta_co2_adsorbded(ngrid))
    558556allocate(co2ice_disappeared(ngrid,nslope))
    559 allocate(porefillingice_thickness_prev_iter(ngrid,nslope))
     557allocate(icetable_thickness_prev_iter(ngrid,nslope))
    560558allocate(delta_h2o_icetablesublim(ngrid))
    561559allocate(delta_h2o_adsorbded(ngrid))
    … …  
    583581call end_adsorption_h_PEM
    584582call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
    585 call end_ice_table_porefilling
    586 call ini_ice_table_porefilling(ngrid,nslope)
     583call end_ice_table
     584call ini_ice_table(ngrid,nslope,nsoilmx_PEM)
    587585
    588586if (soil_pem) then
    … …  
    646644endif
    647645
    648 call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
    649               porefillingice_thickness,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries,          &
    650               tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,global_avg_press_PCM,              &
    651               watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys,delta_co2_adsorbded,                &
    652               h2o_adsorbded_phys,delta_h2o_adsorbded,stratif)
     646call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth, &
     647              icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,        &
     648              ps_timeseries,tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,               &
     649              global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys,          &
     650              delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,stratif)
    653651
    654652! We save the places where h2o ice is sublimating
    … …  
    928926! II_d.3 Update the ice table
    929927        if (icetable_equilibrium) then
    930             write(*,*) "Compute ice table"
    931             porefillingice_thickness_prev_iter = porefillingice_thickness
    932             call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
    933             call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
     928            write(*,*) "Compute ice table (equilibrium method)"
     929            icetable_thickness_prev_iter = icetable_thickness
     930            call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
     931            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
    934932        else if (icetable_dynamic) then
    935             write(*,*) "Compute ice table"
    936             porefillingice_thickness_prev_iter = porefillingice_thickness
    937             call dyn_ss_ice_m_wrapper(ngrid,nsoilmx,TI_PEM,ps,mmol(igcm_h2o_vap),tsoil_PEM,porefillingice_depth,ss_ice_pf,ss_ice_pf_new,porefillingice_depth_new)
    938            
    939             call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
     933            write(*,*) "Compute ice table (dynamic method)"
     934            allocate(porefill(nsoilmx_PEM))
     935            do ig = 1,ngrid
     936                do islope = 1,nslope
     937                    call dyn_ss_ice_m(icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM(ig,:,islope),nsoilmx_PEM,TI_PEM(ig,1,nslope),ps(ig),sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2),ice_porefilling(ig,:,islope),porefill,ssi_depth)
     938                    icetable_depth(ig,islope) = ssi_depth
     939                    ice_porefilling(ig,:,islope) = porefill
     940                enddo
     941            enddo
     942            deallocate(porefill)
     943            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
    940944        endif
    941945
    942946! II_d.4 Update the soil thermal properties
    943         call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
     947        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,TI_PEM)
    944948
    945949! II_d.5 Update the mass of the regolith adsorbed
    … …  
    952956            totmassh2o_adsorbded = 0.
    953957            do ig = 1,ngrid
    954                 do islope =1, nslope
     958                do islope = 1,nslope
    955959                    do l = 1,nsoilmx_PEM
    956                         if (l==1) then
     960                        if (l == 1) then
    957961                            totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &
    958962                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    … …  
    960964                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    961965                        else
    962                             totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
     966                            totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
    963967                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    964                             totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
     968                            totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
    965969                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    966970                        endif
    … …  
    987991        call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
    988992        if (icetable_equilibrium) then
    989             call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,porefillingice_depth(:,islope))
    990             call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,porefillingice_thickness(:,islope))
     993            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
     994            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope))
    991995        else if (icetable_dynamic) then
    992             call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,porefillingice_depth(:,islope))
    993             call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,porefillingice_thickness(:,islope))
     996            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
     997            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope))
    994998        endif
    995999
    … …  
    12051209call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
    12061210call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
    1207              TI_PEM,porefillingice_depth,porefillingice_thickness,       &
     1211             TI_PEM,icetable_depth,icetable_thickness,ice_porefilling,   &
    12081212             co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif)
    12091213write(*,*) "restartpem.nc has been written"
    … …  
    12271231deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved)
    12281232deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg)
    1229 deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
     1233deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,icetable_thickness_prev_iter)
    12301234deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif)
    12311235!----------------------------- END OUTPUT ------------------------------

  • trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90

    r3490 r3493  
    77!=======================================================================
    88
    9 SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table,ice_table_thickness, &
    10                     tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,      &
    11                     global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,              &
     9SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table_depth,ice_table_thickness,      &
     10                    ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice, &
     11                    global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,                         &
    1212                    m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif)
    1313
    … …  
    6060real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
    6161real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
    62 real, dimension(ngrid,nslope),                   intent(inout) :: ice_table           ! Ice table depth [m]
     62real, dimension(ngrid,nslope),                   intent(inout) :: ice_table_depth     ! Ice table depth [m]
    6363real, dimension(ngrid,nslope),                   intent(inout) :: ice_table_thickness ! Ice table thickness [m]
     64real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: ice_porefilling     ! Subsurface ice pore filling [m3/m3]
    6465real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst          ! instantaneous soil (mid-layer) temperature [K]
    6566real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
    … …  
    107108
    108109!0.2 Set to default values
    109 ice_table = -1. ! by default, no ice table
     110ice_table_depth = -1. ! by default, no ice table
    110111ice_table_thickness = -1.
    111112!1. Run
    … …  
    225226            do ig = 1,ngrid
    226227                if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
    227                     inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
    228                                                               (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
     228                    inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
     229                                                              (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
    229230                                                              ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
    230231
    … …  
    308309!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    309310!3. Ice Table
    310         if (icetable_equilibrium .or. icetable_dynamic) then
    311             call get_field("ice_table",ice_table,found)
     311        if (icetable_equilibrium) then
     312            call get_field("ice_table_depth",ice_table_depth,found)
    312313            if (.not. found) then
    313                 write(*,*)'PEM settings: failed loading <ice_table>'
     314                write(*,*)'PEM settings: failed loading <ice_table_depth>'
    314315                write(*,*)'will reconstruct the values of the ice table given the current state'
    315                 call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table,ice_table_thickness)
    316                 call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
     316                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
     317                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
    317318                do islope = 1,nslope
    318319                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
    319320                enddo
     321            endif
     322            write(*,*) 'PEMETAT0: ICE TABLE done'
     323        else if (icetable_dynamic) then
     324            call get_field("ice_table_depth",ice_table_depth,found)
     325            if (.not. found) then
     326                write(*,*)'PEM settings: failed loading <ice_table_depth>'
     327                write(*,*)'will reconstruct the values of the ice table given the current state'
     328                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
     329                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
     330                do islope = 1,nslope
     331                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
     332                enddo
     333            endif
     334            call get_field("ice_porefilling",ice_porefilling,found)
     335            if (.not. found) then
     336                write(*,*)'PEM settings: failed loading <ice_porefilling>'
     337                ice_porefilling = 0.
    320338            endif
    321339            write(*,*) 'PEMETAT0: ICE TABLE done'
    … …  
    354372    call close_startphy
    355373
    356 else !No startfi, let's build all by hand
     374else !No startpem, let's build all by hand
    357375
    358376    write(*,*)'There is no "startpem.nc"!'
    … …  
    484502!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    485503!c) Ice table
    486         if (icetable_equilibrium .or. icetable_dynamic) then
    487             call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table,ice_table_thickness)
    488             call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
     504        if (icetable_equilibrium) then
     505            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
     506            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
     507            do islope = 1,nslope
     508                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
     509            enddo
     510            write(*,*) 'PEMETAT0: Ice table done'
     511        else if (icetable_dynamic) then
     512            ice_porefilling = 0.
     513            ice_table_depth = 0.
     514            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
    489515            do islope = 1,nslope
    490516                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))

  • trunk/LMDZ.COMMON/libf/evolution/pemredem.F90

    r3319 r3493  
    5858
    5959SUBROUTINE pemdem1(filename,i_myear,nsoil_PEM,ngrid,nslope,tsoil_slope_PEM,inertiesoil_slope_PEM, &
    60                    ice_table_depth,ice_table_thickness,m_co2_regolith,m_h2o_regolith,h2o_ice,stratif)
     60                   ice_table_depth,ice_table_thickness,ice_porefilling,m_co2_regolith,m_h2o_regolith,h2o_ice,stratif)
    6161
    6262! write time-dependent variable to restart file
    … …  
    7878real, dimension(ngrid,nslope),           intent(in) :: ice_table_depth       ! under mesh bining according to slope
    7979real, dimension(ngrid,nslope),           intent(in) :: ice_table_thickness   ! under mesh bining according to slope
     80real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: ice_porefilling       ! under mesh bining according to slope
    8081real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: m_co2_regolith, m_h2o_regolith
    8182real, dimension(ngrid,nslope),           intent(in) :: h2o_ice
    … …  
    122123    call put_field("ice_table_depth","Depth of ice table",ice_table_depth,Year)
    123124    call put_field("ice_table_thickness","Depth of ice table",ice_table_thickness,Year)
     125    call put_field("ice_porefilling","Subsurface ice pore filling",ice_porefilling,Year)
    124126    call put_field("inertiedat_PEM","Thermal inertie of PEM ",inertiedat_PEM,Year)
    125127endif ! soil_pem
Note: See TracChangeset for help on using the changeset viewer.