Ignore:
Timestamp:
Nov 20, 2024, 3:53:19 PM (5 weeks 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

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
4 deleted
7 edited

Legend:

Unmodified
Added
Removed

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

    r3526 r3527  
     1MODULE dyn_ss_ice_m_mod
     2
     3implicit none
     4
     5CONTAINS
     6
    17!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    28!!!
    … …  
    1521! orbital elements remain constant
    1622!***********************************************************************
    17   use miscparameters, only : pi, d2r, NMAX, marsDay, solsperyear
    18   !use allinterfaces
     23  use constants_marspem_mod, only: sec_per_sol
     24  use fast_subs_mars, only: psv, icelayer_mars, NMAX
     25#ifndef CPP_STD
     26    use comcstfi_h,   only: pi
     27#else
     28    use comcstfi_mod, only: pi
     29#endif
    1930  implicit none
    2031  integer, parameter :: NP=1   ! # of sites
    … …  
    105116 !    print *,'  latitude (deg)',latitude(k),' rho*c (J/m^3/K)',rhoc(k),' thIn=',thIn(k)
    106117 !    print *,'  total pressure=',p0(k),'partial pressure=',pfrost(k)
    107      delta = thIn(k)/rhoc(k)*sqrt(marsDay/pi)
     118     delta = thIn(k)/rhoc(k)*sqrt(sec_per_sol/pi)
    108119 !    print *,'  skin depths (m)',delta,delta*sqrt(solsperyear)
    109120     call soilthprop(porosity,1.d0,rhoc(k),thIn(k),1,newrhoc,newti,icefrac)
    … …  
    170181 
    171182end subroutine dyn_ss_ice_m
     183
     184!=======================================================================
     185
     186subroutine soilthprop(porosity,fill,rhocobs,tiobs,layertype, &
     187     &     newrhoc,newti,icefrac)
     188!***********************************************************************
     189! soilthprop: assign thermal properties of icy soil or dirty ice
     190!
     191!     porositiy = void space / total volume
     192!     rhof = density of free ice in space not occupied by regolith [kg/m^3]
     193!     fill = rhof/icedensity <=1 (only relevant for layertype 1)
     194!     rhocobs = heat capacity per volume of dry regolith [J/m^3]
     195!     tiobs = thermal inertia of dry regolith [SI-units]
     196!     layertype: 1=interstitial ice, 2=pure ice or ice with dirt
     197!                3=pure ice, 4=ice-cemented soil, 5=custom values
     198!     icefrac: fraction of ice in icelayer (only relevant for layertype 2)
     199!     output are newti and newrhoc
     200!***********************************************************************
     201  implicit none
     202  integer, intent(IN) :: layertype
     203  real(8), intent(IN) :: porosity, fill, rhocobs, tiobs
     204  real(8), intent(OUT) :: newti, newrhoc
     205  real(8), intent(IN) :: icefrac
     206  real(8) kobs, cice, icedensity, kice
     207  !parameter (cice=2000.d0, icedensity=926.d0, kice=2.4d0) ! unaffected by scaling
     208  parameter (cice=1540.d0, icedensity=927.d0, kice=3.2d0) ! at 198 Kelvin
     209  real(8) fA, ki0, ki, k
     210  real(8), parameter :: kw=3. ! Mellon et al., JGR 102, 19357 (1997)
     211
     212  kobs = tiobs**2/rhocobs
     213  ! k, rhoc, and ti are defined in between grid points
     214  ! rhof and T are defined on grid points
     215
     216  newrhoc = -9999.
     217  newti  = -9999.
     218
     219  select case (layertype)
     220  case (1) ! interstitial ice
     221     newrhoc = rhocobs + porosity*fill*icedensity*cice
     222     if (fill>0.) then
     223        !--linear addition (option A)
     224        k = porosity*fill*kice + kobs
     225        !--Mellon et al. 1997 (option B)
     226        ki0 = porosity/(1/kobs-(1-porosity)/kw)
     227        fA = sqrt(fill)
     228        ki = (1-fA)*ki0 + fA*kice
     229        !k = kw*ki/((1-porosity)*ki+porosity*kw)
     230     else
     231        k = kobs
     232     endif
     233     newti = sqrt(newrhoc*k)
     234     
     235  case (2)  ! massive ice (pure or dirty ice)
     236     newrhoc = rhocobs*(1.-icefrac)/(1.-porosity) + icefrac*icedensity*cice
     237     k = icefrac*kice + (1.-icefrac)*kw
     238     newti = sqrt(newrhoc*k)
     239 
     240  case (3)  ! all ice, special case of layertype 2, which doesn't use porosity
     241     newrhoc = icedensity*cice
     242     k = kice
     243     newti = sqrt(newrhoc*k)
     244 
     245  case (4)  ! pores completely filled with ice, special case of layertype 1
     246     newrhoc = rhocobs + porosity*icedensity*cice
     247     k = porosity*kice + kobs ! option A, end-member case of type 1, option A
     248     !k = kw*kice/((1-porosity)*kice+porosity*kw) ! option B, harmonic average
     249     newti = sqrt(newrhoc*k)
     250
     251  case (5)  ! custom values
     252     ! values from Mellon et al. (2004) for ice-cemented soil
     253     newrhoc = 2018.*1040.
     254     k = 2.5
     255     newti = sqrt(newrhoc*k)
     256
     257  case default
     258     error stop 'invalid layer type'
     259     
     260  end select
     261 
     262end subroutine soilthprop
     263
     264
     265!=======================================================================
     266     
     267      real*8 function frostpoint(p)
     268!     inverse of psv
     269!     input is partial pressure [Pascal]
     270!     output is temperature [Kelvin]
     271      implicit none
     272      real*8 p
     273     
     274!-----inverse of parametrization 1
     275!      real*8 DHmelt,DHvap,DHsub,R,pt,Tt
     276!      parameter (DHmelt=6008.,DHvap=45050.)
     277!      parameter (DHsub=DHmelt+DHvap)
     278!      parameter (R=8.314,pt=6.11e2,Tt=273.16)
     279!      frostpoint = 1./(1./Tt-R/DHsub*log(p/pt))
     280     
     281!-----inverse of parametrization 2
     282!     inverse of eq. (2) in Murphy & Koop (2005)
     283      real*8 A,B
     284      parameter (A=-6143.7, B=28.9074)
     285      frostpoint = A / (log(p) - B)
     286
     287!-----approximate inverse of parametrization 3
     288!     eq. (8) in Murphy & Koop (2005)
     289!      frostpoint = (1.814625*log(p) + 6190.134)/(29.120 - log(p))
     290     
     291      end function frostpoint
     292
     293END MODULE dyn_ss_ice_m_mod

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

    r3526 r3527  
    1 module thermalmodelparam_mars
     1MODULE fast_subs_mars
     2
     3implicit none
     4
    25  ! parameters for thermal model
    36  ! they are only used in the subroutines below
    … …  
    811  real(8), parameter :: emiss0 = 1.  ! emissivity of dry surface
    912  integer, parameter :: EQUILTIME =15 ! [Mars years]
    10 end module thermalmodelparam_mars
    11 
    12 
     13 
     14  integer, parameter :: NMAX = 1000
     15
     16CONTAINS
    1317!*****************************************************
    1418! Subroutines for fast method
    … …  
    2529! latitude  [degree]
    2630!*************************************************************************
    27   use miscparameters, only : d2r, NMAX, icedensity
    28   use allinterfaces, except_this_one => icelayer_mars
     31  use ice_table_mod, only: rho_ice
     32  use fast_subs_univ, only: icechanges
    2933  !use omp_lib
    3034  implicit none
    … …  
    5357     Diff = 4e-4*600./p0(k)
    5458     fracIR = 0.04*p0(k)/600.; fracDust = 0.02*p0(k)/600.
    55      B = Diff*bigstep*86400.*365.24/(porosity*icedensity)
     59     B = Diff*bigstep*86400.*365.24/(porosity*927.)
     60     !B = Diff*bigstep*86400.*365.24/(porosity*rho_ice(Tb(),'h2o'))
    5661     
    5762     typeT = -9
    … …  
    140145!  Tb>0 initializes temperatures with Tb
    141146!***********************************************************************
    142   use miscparameters
    143   use thermalmodelparam_mars
    144   use allinterfaces, except_this_one => ajsub_mars
     147  use fast_subs_univ, only: depths_avmeth, equildepth
     148  use constants_marspem_mod, only: sols_per_my, sec_per_sol
    145149  implicit none
    146150  integer, intent(IN) :: nz, typeT
    … …  
    166170  real(8) Tsurf, Tco2frost, albedo, Fsurf, m, dE, emiss, T(NMAX)
    167171  real(8) Told(nz), Fsurfold, Tsurfold, Tmean0, avrho2
    168   real(8) rhosatav0, rhosatav(nz), rlow
    169   real(8), external :: psv, tfrostco2
     172  real(8) rhosatav0, rhosatav(nz), rlow
    170173 
    171   tmax = EQUILTIME*solsperyear
     174  tmax = EQUILTIME*sols_per_my
    172175  nsteps = int(tmax/dt)     ! calculate total number of timesteps
    173176
    … …  
    207210  do n=0,nsteps-1
    208211     time = (n+1)*dt         !   time at n+1
    209      tdays = time*(marsDay/earthDay) ! parenthesis may improve roundoff
     212   !  tdays = time*(sec_per_sol/earthDay) ! parenthesis may improve roundoff
    210213   !  call generalorbit(tdays,a,ecc,omega,eps,marsLs,marsDec,marsR)
    211214   !  HA = 2.*pi*mod(time,1.d0)  ! hour angle
    … …  
    217220     Told(1:nz) = T(1:nz)
    218221     if (m<=0. .or. Tsurf>Tco2frost) then
    219        ! call conductionQ(nz,z,dt*marsDay,Qn,Qnp1,T,ti,rhocv,emiss, &
     222       ! call conductionQ(nz,z,dt*sec_per_sol,Qn,Qnp1,T,ti,rhocv,emiss, &
    220223       !      &           Tsurf,Fgeotherm,Fsurf)
    221224     endif
    222225     if (Tsurf<Tco2frost .or. m>0.) then ! CO2 condensation
    223226        T(1:nz) = Told(1:nz)
    224         !call conductionT(nz,z,dt*marsDay,T,Tsurfold,Tco2frost,ti, &
     227        !call conductionT(nz,z,dt*sec_per_sol,T,Tsurfold,Tco2frost,ti, &
    225228             !&              rhocv,Fgeotherm,Fsurf)
    226229        Tsurf = Tco2frost
    227230    !    dE = (- Qn - Qnp1 + Fsurfold + Fsurf + &
    228231    !         &      emiss*sigSB*(Tsurfold**4+Tsurf**4)/2.
    229         m = m + dt*marsDay*dE/Lco2frost
     232        m = m + dt*sec_per_sol*dE/Lco2frost
    230233     endif
    231234     if (Tsurf>Tco2frost .or. m<=0.) then
    … …  
    238241     !Qn=Qnp1
    239242     
    240      if (time>=tmax-solsperyear) then
     243     if (time>=tmax-sols_per_my) then
    241244        Tmean1 = Tmean1 + Tsurf
    242245        Tmean3 = Tmean3 + T(nz)
    … …  
    301304tfrostco2 = 3182.48/(23.3494+log(100./p))
    302305end function
     306
     307!=======================================================================
     308
     309      real*8 function psv(T)
     310!     saturation vapor pressure of H2O ice [Pascal]
     311!     input is temperature [Kelvin]
     312      implicit none
     313      real*8 T
     314
     315!-----parametrization 1
     316!      real*8 DHmelt,DHvap,DHsub,R,pt,Tt,C
     317!      parameter (DHmelt=6008.,DHvap=45050.)
     318!      parameter (DHsub=DHmelt+DHvap) ! sublimation enthalpy [J/mol]
     319!      parameter (R=8.314,pt=6.11e2,Tt=273.16)
     320!      C = (DHsub/R)*(1./T - 1./Tt)
     321!      psv = pt*exp(-C)
     322
     323!-----parametrization 2
     324!     eq. (2) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
     325!     differs from parametrization 1 by only 0.1%
     326      real*8 A,B
     327      parameter (A=-6143.7, B=28.9074)
     328      psv = exp(A/T+B)  ! Clapeyron
     329
     330!-----parametrization 3     
     331!     eq. (7) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
     332!     psv = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T)
     333     
     334      end function psv
     335
     336END MODULE fast_subs_mars

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

    r3526 r3527  
     1MODULE fast_subs_univ
     2
     3implicit none
     4
     5CONTAINS
     6
    17!*****************************************************
    28! Commonly used subroutines for fast method
    … …  
    1925!  this is not the true (final) equilibrium depth
    2026!***********************************************************************
    21   use allinterfaces, except_this_one => equildepth
    2227  implicit none
    2328  integer, intent(IN) :: nz
    … …  
    5459!  B = Diff*bigstep/(porosity*icedensity)  [SI units]
    5560!***********************************************************************
    56   use allinterfaces, except_this_one => depths_avmeth
    5761  use math_mod, only: deriv2_simple, deriv1_onesided, deriv1, colint
     62  use ice_table_mod, only: constriction
    5863  implicit none
    5964  integer, intent(IN) :: nz, typeT
    … …  
    178183
    179184
    180 pure function constriction(porefill)
    181 ! specify constriction function here, 0<=eta<=1
    182   implicit none
    183   real(8), intent(IN) :: porefill
    184   real(8) eta, constriction
    185   if (porefill<=0.) eta = 1.
    186   if (porefill>0. .and. porefill<1.) then
    187      ! eta = 1.
    188      ! eta = 1-porefill
    189      eta = (1-porefill)**2  ! Hudson et al., JGR, 2009
    190   endif
    191   if (porefill>=1.) eta = 0.
    192   constriction = eta
    193 end function constriction
    194 
    195 
    196 
    197185pure subroutine icechanges(nz,z,typeF,avdrho,avdrhoP,ypp, &
    198186     & Diff,porosity,icefrac,bigstep,zdepthT,porefill,typeG)
    … …  
    202190! a crucial subroutine
    203191!***********************************************************
    204   use miscparameters, only : icedensity
    205   use allinterfaces, except_this_one => icechanges
    206192  use math_mod, only: colint
     193  use ice_table_mod, only: rho_ice
    207194  implicit none
    208195  integer, intent(IN) :: nz, typeF, typeG
    … …  
    213200  real(8) B, beta, integ
    214201
    215   B = Diff*bigstep*86400.*365.24/(porosity*icedensity)
     202  B = Diff*bigstep*86400.*365.24/(porosity*927.)
     203  !B = Diff*bigstep*86400.*365.24/(porosity*rho_ice(T,'h2o'))
    216204
    217205  ! advance ice table, avdrho>0 is retreat
    … …  
    232220     if (typeP==typeT) then   ! new 2011-09-01
    233221        beta = (1-icefrac)/(1-porosity)/icefrac
    234         beta = Diff*bigstep*beta*86400*365.24/icedensity
     222        beta = Diff*bigstep*beta*86400*365.24/927.
     223        !beta = Diff*bigstep*beta*86400*365.24/rho_ice(T,'h2o')
    235224        zdepthT = sqrt(2*beta*avdrho*18./8314. + zdepthT**2)
    236225     endif
    … …  
    283272  end if
    284273end subroutine icechanges
     274
     275END MODULE fast_subs_univ

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

    r3526 r3527  
    495495== 19/11/2024 == JBC
    496496Removing unused or redundant Norbert Schorghofer's subroutines (follow-up of r3493) + cleaning and some modifications of related subroutines.
     497
     498== 20/11/2024 == JBC
     499- Removing redundant Norbert Schorghofer's subroutines/parameters (follow-up of r3526);
     500- Making all Norbert Schorghofer's subroutines with modern explicit interface via modules;
     501- Cleaning of "glaciers_mod.F90";
     502- Optimization for the computation of ice density according to temperature by using a function.

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

    r3345 r3527  
    22
    33implicit none
    4 
    54
    65! Flags for ice management
    … …  
    3332implicit none
    3433
    35 ! arguments
    36 ! ---------
    37 
    3834! Inputs:
    39       integer,                           intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
    40       real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Physical points x Slopes: Distribution of the subgrid slopes
    41       real,    dimension(ngrid),         intent(in) :: def_slope_mean                ! Physical points: values of the sub grid slope angles
    42       real,    dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM                   ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
    43       real,    dimension(ngrid,timelen), intent(in) :: ps_PCM                        ! Physical x Time field: surface pressure given by the PCM [Pa]
    44       real,                              intent(in) :: global_avg_ps_PCM             ! Global averaged surface pressure from the PCM [Pa]
    45       real,                              intent(in) :: global_avg_ps_PEM             ! global averaged surface pressure during the PEM iteration [Pa]
    46      
     35!--------
     36integer,                           intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
     37real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Physical points x Slopes: Distribution of the subgrid slopes
     38real,    dimension(ngrid),         intent(in) :: def_slope_mean                ! Physical points: values of the sub grid slope angles
     39real,    dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM                   ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
     40real,    dimension(ngrid,timelen), intent(in) :: ps_PCM                        ! Physical x Time field: surface pressure given by the PCM [Pa]
     41real,                              intent(in) :: global_avg_ps_PCM             ! Global averaged surface pressure from the PCM [Pa]
     42real,                              intent(in) :: global_avg_ps_PEM             ! global averaged surface pressure during the PEM iteration [Pa]
     43
    4744! Ouputs:
    48       real, dimension(ngrid,nslope), intent(inout) :: co2ice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
    49       real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow      ! flag to see if there is flow on the subgrid slopes
    50       real, dimension(ngrid),        intent(inout) :: flag_co2flow_mesh ! same but within the mesh
    51 
     45!--------
     46real, dimension(ngrid,nslope), intent(inout) :: co2ice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
     47real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow      ! flag to see if there is flow on the subgrid slopes
     48real, dimension(ngrid),        intent(inout) :: flag_co2flow_mesh ! same but within the mesh
    5249! Local
    53       real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
    54       real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow
    55 
    56 !-----------------------------
     50!------
     51real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
     52real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow
     53
    5754write(*,*) "Flow of CO2 glaciers"
    58    
    5955call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond)
    6056call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
    … …  
    8379
    8480! Inputs:
    85       integer,intent(in) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
    86       real,intent(in) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles
    87       real,intent(in) :: Tice(ngrid,nslope) ! Ice Temperature [K]
     81integer,                       intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
     82real, dimension(ngrid,nslope), intent(in) :: subslope_dist  ! Physical points x Slopes : Distribution of the subgrid slopes
     83real, dimension(ngrid),        intent(in) :: def_slope_mean ! Slopes: values of the sub grid slope angles
     84real, dimension(ngrid,nslope), intent(in) :: Tice           ! Ice Temperature [K]
    8885! Ouputs:
    89       real,intent(inout) :: h2oice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
    90       real,intent(inout) :: flag_h2oflow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
    91       real,intent(inout) :: flag_h2oflow_mesh(ngrid) ! same but within the mesh
     86real, dimension(ngrid,nslope), intent(inout) :: h2oice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
     87real, dimension(ngrid,nslope), intent(inout) :: flag_h2oflow      ! flag to see if there is flow on the subgrid slopes
     88real, dimension(ngrid),        intent(inout) :: flag_h2oflow_mesh ! same but within the mesh
    9289! Local
    93       real :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow
    94 
    95 !-----------------------------
     90real, dimension(ngrid,nslope) :: hmax ! Physical x Slope field: maximum thickness for co2  glacier before flow
     91
    9692write(*,*) "Flow of H2O glaciers"
    97 
    9893call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
    9994call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh)
    … …  
    105100SUBROUTINE compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
    106101
     102use ice_table_mod, only: rho_ice
    107103use abort_pem_mod, only: abort_pem
    108104#ifndef CPP_STD
    … …  
    114110!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    115111!!!
    116 !!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle
    117 !!!          before initating flow
     112!!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle before initating flow
    118113!!!         
    119114!!! Author: LL,based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
    … …  
    123118implicit none
    124119
    125 ! arguments
    126 ! --------
    127 
    128120! Inputs
    129       integer,intent(in) :: ngrid,nslope ! # of grid points and subslopes
    130       integer,intent(in) :: iflat        ! index of the flat subslope
    131       real,intent(in) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg]
    132       real,intent(in) :: Tice(ngrid,nslope)     ! Physical field:  ice temperature [K]
    133       character(len=3), intent(in) :: name_ice ! Nature of the ice
     121! ------
     122integer,                       intent(in) :: ngrid, nslope  ! # of grid points and subslopes
     123integer,                       intent(in) :: iflat          ! index of the flat subslope
     124real, dimension(nslope),       intent(in) :: def_slope_mean ! Slope field: Values of the subgrid slope angles [deg]
     125real, dimension(ngrid,nslope), intent(in) :: Tice           ! Physical field:  ice temperature [K]
     126character(3),                  intent(in) :: name_ice       ! Nature of ice
    134127! Outputs
    135       real,intent(out) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum  thickness before flaw [m]
     128! -------
     129real, dimension(ngrid,nslope), intent(out) :: hmax ! Physical grid x Slope field: maximum  thickness before flaw [m]
    136130! Local
    137       DOUBLE PRECISION :: tau_d    ! characteristic basal drag, understood as the stress that an ice mass flowing under its weight balanced by viscosity. Value obtained from I.Smith
    138       real :: rho(ngrid,nslope) ! co2 ice density [kg/m^3]
    139       integer :: ig,islope ! loop variables
    140       real :: slo_angle
    141 
    142 ! 1. Compute rho
    143     if(name_ice.eq."co2") then
    144       DO ig = 1,ngrid
    145         DO islope = 1,nslope
    146         rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3  ! Mangan et al. 2017
     131! -----
     132real    :: tau_d      ! characteristic basal drag, understood as the stress that an ice mass flowing under its weight balanced by viscosity. Value obtained from I.Smith
     133integer :: ig, islope ! loop variables
     134real    :: slo_angle
     135
     136select case (trim(adjustl(name_ice)))
     137    case('h2o')
     138        tau_d = 1.e5
     139    case('co2')
    147140        tau_d = 5.e3
    148         ENDDO
    149       ENDDO
    150     elseif (name_ice.eq."h2o") then
    151       DO ig = 1,ngrid
    152         DO islope = 1,nslope
    153           rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
    154           tau_d = 1.e5
    155         ENDDO
    156       ENDDO
    157     else
    158       call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1)
    159     endif
    160 
    161 ! 3. Compute max thickness
    162     DO ig = 1,ngrid
    163        DO islope = 1,nslope
    164           if(islope.eq.iflat) then
     141    case default
     142        call abort_pem("compute_hmaxglaciers","Type of ice not known!",1)
     143end select
     144
     145do ig = 1,ngrid
     146    do islope = 1,nslope
     147        if (islope == iflat) then
    165148            hmax(ig,islope) = 1.e8
    166           else
     149        else
    167150            slo_angle = abs(def_slope_mean(islope)*pi/180.)
    168             hmax(ig,islope) = tau_d/(rho(ig,islope)*g*slo_angle)
    169           endif
    170        ENDDO
    171     ENDDO
     151            hmax(ig,islope) = tau_d/(rho_ice(Tice(ig,islope),name_ice)*g*slo_angle)
     152        endif
     153    enddo
     154enddo
     155
    172156END SUBROUTINE compute_hmaxglaciers
    173157
    … …  
    183167!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    184168
     169use ice_table_mod, only: rho_ice
    185170use abort_pem_mod, only: abort_pem
    186171#ifndef CPP_STD
    … …  
    192177implicit none
    193178
    194 ! arguments
    195 ! --------
    196 
    197179! Inputs
    198       integer, intent(in) :: ngrid,nslope               !# of physical points and subslope
    199       integer, intent(in) :: iflat                      ! index of the flat subslope
    200       real, intent(in) :: subslope_dist(ngrid,nslope)   ! Distribution of the subgrid slopes within the mesh
    201       real, intent(in) :: def_slope_mean(nslope)        ! values of the subgrid slopes
    202       real, intent(in) :: hmax(ngrid,nslope)            ! maximum height of the  glaciers before initiating flow [m]
    203       real, intent(in) :: Tice(ngrid,nslope)            ! Ice temperature[K]
    204       character(len=3), intent(in) :: name_ice              ! Nature of the ice
    205 
     180!-------
     181integer,                       intent(in) :: ngrid, nslope  ! # of physical points and subslope
     182integer,                       intent(in) :: iflat          ! index of the flat subslope
     183real, dimension(ngrid,nslope), intent(in) :: subslope_dist  ! Distribution of the subgrid slopes within the mesh
     184real, dimension(nslope),       intent(in) :: def_slope_mean ! values of the subgrid slopes
     185real, dimension(ngrid,nslope), intent(in) :: hmax           ! maximum height of the  glaciers before initiating flow [m]
     186real, dimension(ngrid,nslope), intent(in) :: Tice           ! Ice temperature[K]
     187character(3),                  intent(in) :: name_ice       ! Nature of the ice
    206188! Outputs
    207       real, intent(inout) :: qice(ngrid,nslope)      ! CO2 in the subslope [kg/m^2]
    208       real, intent(inout) :: flag_flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope
    209       real, intent(inout) :: flag_flowmesh(ngrid)    ! boolean to check if there is flow in the mesh
     189!--------
     190real, dimension(ngrid,nslope), intent(inout) :: qice          ! CO2 in the subslope [kg/m^2]
     191real, dimension(ngrid,nslope), intent(inout) :: flag_flow     ! boolean to check if there is flow on a subgrid slope
     192real, dimension(ngrid),        intent(inout) :: flag_flowmesh ! boolean to check if there is flow in the mesh
    210193! Local
    211       integer ig,islope ! loop
    212       real rho(ngrid,nslope) ! density of ice, temperature dependant [kg/m^3]
    213       integer iaval ! ice will be transfered here
    214 
    215 ! 0. Compute rho
    216     if(name_ice.eq."co2") then
    217       DO ig = 1,ngrid
    218         DO islope = 1,nslope
    219         rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3  ! Mangan et al. 2017
    220         ENDDO
    221       ENDDO
    222     elseif (name_ice.eq."h2o") then
    223       DO ig = 1,ngrid
    224         DO islope = 1,nslope
    225           rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
    226         ENDDO
    227       ENDDO
    228     else
    229       call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1)
    230     endif
    231 
    232 ! 1. Compute the transfer of ice
    233 
    234        DO ig = 1,ngrid
    235         DO islope = 1,nslope
    236           IF(islope.ne.iflat) THEN ! ice can be infinite on flat ground
     194!------
     195integer :: ig, islope ! loop
     196integer :: iaval      ! ice will be transfered here
     197
     198do ig = 1,ngrid
     199    do islope = 1,nslope
     200        if (islope /= iflat) then ! ice can be infinite on flat ground
    237201! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely  on flat ground
    238             IF(qice(ig,islope).ge.rho(ig,islope)*hmax(ig,islope) * &
    239                   cos(pi*def_slope_mean(islope)/180.)) THEN
     202            if (qice(ig,islope) >= rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)) then
    240203! Second: determine the flatest slopes possible:
    241                 IF(islope.gt.iflat) THEN
    242                   iaval=islope-1
    243                 ELSE
    244                  iaval=islope+1
    245                 ENDIF
    246                 do while ((iaval.ne.iflat).and.  &
    247                     (subslope_dist(ig,iaval).eq.0))
    248                   IF(iaval.gt.iflat) THEN
    249                      iaval=iaval-1
    250                   ELSE
    251                      iaval=iaval+1
    252                   ENDIF
     204                if (islope > iflat) then
     205                    iaval=islope-1
     206                else
     207                    iaval = islope + 1
     208                endif
     209                do while (iaval /= iflat .and. subslope_dist(ig,iaval) == 0)
     210                    if (iaval > iflat) then
     211                        iaval = iaval - 1
     212                    else
     213                        iaval = iaval + 1
     214                    endif
    253215                enddo
    254               qice(ig,iaval) = qice(ig,iaval) + &
    255                (qice(ig,islope) - rho(ig,islope)*hmax(ig,islope) *     &
    256                cos(pi*def_slope_mean(islope)/180.)) *             &
    257                subslope_dist(ig,islope)/subslope_dist(ig,iaval) * &
    258                cos(pi*def_slope_mean(iaval)/180.) /               &
    259                cos(pi*def_slope_mean(islope)/180.)
    260 
    261               qice(ig,islope)=rho(ig,islope)*hmax(ig,islope) *        &
    262                cos(pi*def_slope_mean(islope)/180.)
    263 
    264               flag_flow(ig,islope) = 1.
    265               flag_flowmesh(ig) = 1.
    266             ENDIF ! co2ice > hmax
    267           ENDIF ! iflat
    268         ENDDO !islope
    269        ENDDO !ig
     216                qice(ig,iaval) = qice(ig,iaval) + (qice(ig,islope) - rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)) &
     217                               *subslope_dist(ig,islope)/subslope_dist(ig,iaval)*cos(pi*def_slope_mean(iaval)/180.)/cos(pi*def_slope_mean(islope)/180.)
     218
     219                qice(ig,islope) = rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)
     220                flag_flow(ig,islope) = 1.
     221                flag_flowmesh(ig) = 1.
     222            endif ! co2ice > hmax
     223        endif ! iflat
     224    enddo !islope
     225enddo !ig
     226
    270227END SUBROUTINE
    271228
    … …  
    281238!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    282239 
    283 use constants_marspem_mod,only : alpha_clap_co2,beta_clap_co2
     240use constants_marspem_mod, only: alpha_clap_co2,beta_clap_co2
    284241
    285242implicit none
    … …  
    302259real    :: ave    ! Intermediate to compute average
    303260
    304 !!!!!!!!!!!!!!!!!!!!!!!!!!!!
    305261do ig = 1,ngrid
    306262    ave = 0

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

    r3525 r3527  
    184184        do islope = 1,nslope
    185185            call compute_Tice_pem(nsoil,tsoil(ig,:,islope),tsurf(ig,islope),icetable_depth(ig,islope),Tice)
    186             rho = -3.5353e-4*Tice**2 + 0.0351*Tice + 933.5030 ! Rottgers, 2012
    187             delta_m_h2o(ig) = delta_m_h2o(ig) + porosity*rho*(icetable_thickness(ig,islope) - icetable_thickness_old(ig,islope)) & ! convention > 0. <=> it condenses
     186            delta_m_h2o(ig) = delta_m_h2o(ig) + porosity*rho_ice(Tice,'h2o')*(icetable_thickness(ig,islope) - icetable_thickness_old(ig,islope)) & ! convention > 0. <=> it condenses
    188187                              *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
    189188        enddo
    … …  
    195194            do isoil = 1,nsoil
    196195                call compute_Tice_pem(nsoil,tsoil(ig,:,islope),tsurf(ig,islope),mlayer_PEM(isoil - 1),Tice)
    197                 rho = -3.5353e-4*Tice**2 + 0.0351*Tice + 933.5030 ! Rottgers, 2012
    198                 delta_m_h2o(ig) = delta_m_h2o(ig) + rho*(ice_porefilling(ig,isoil,islope) - ice_porefilling_old(ig,isoil,islope)) & ! convention > 0. <=> it condenses
     196                delta_m_h2o(ig) = delta_m_h2o(ig) + rho_ice(Tice,'h2o')*(ice_porefilling(ig,isoil,islope) - ice_porefilling_old(ig,isoil,islope)) & ! convention > 0. <=> it condenses
    199197                                  *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
    200198            enddo
    … …  
    258256    index_tmp = nsoilmx_PEM
    259257    do ilay = 1,nsoilmx_PEM
    260         call constriction(porefill(ilay),eta(ilay))
     258        eta(ilay) = constriction(porefill(ilay))
    261259    enddo
    262260else
    263261    index_tmp = index_IS
    264262    do ilay = 1,index_IS - 1
    265         call constriction(porefill(ilay),eta(ilay))
     263        eta(ilay) = constriction(porefill(ilay))
    266264    enddo
    267265    do ilay = index_IS,nsoilmx_PEM
    … …  
    350348
    351349!-----------------------------------------------------------------------
    352 SUBROUTINE constriction(porefill,eta)
     350FUNCTION constriction(porefill) RESULT(eta)
    353351
    354352!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    363361
    364362real, intent(in) :: porefill ! pore filling fraction
    365 real, intent(out) :: eta ! constriction
     363real :: eta ! constriction
    366364
    367365if (porefill <= 0.) then
    … …  
    373371endif
    374372
    375 END SUBROUTINE constriction
     373END FUNCTION constriction
    376374
    377375!-----------------------------------------------------------------------
    … …  
    386384
    387385use comsoil_h_PEM, only: layer_PEM, mlayer_PEM
     386use abort_pem_mod, only: abort_pem
    388387
    389388implicit none
    … …  
    422421    endif
    423422    if (indexice < 0) then
    424         call abort_physic("compute_Tice_pem","subsurface ice is below the last soil layer",1)
     423        call abort_pem("compute_Tice_pem","Subsurface ice is below the last soil layer!",1)
    425424    else
    426425        if(indexice >= 1) then ! Linear inteprolation between soil temperature
    … …  
    434433END SUBROUTINE compute_Tice_pem
    435434
     435!-----------------------------------------------------------------------
     436FUNCTION rho_ice(T,ice) RESULT(rho)
     437!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     438!!!
     439!!! Purpose: Compute subsurface ice density
     440!!!
     441!!! Author: JBC
     442!!!
     443!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     444
     445use abort_pem_mod, only: abort_pem
     446
     447implicit none
     448
     449real,         intent(in) :: T
     450character(3), intent(in) :: ice
     451real :: rho
     452
     453select case (trim(adjustl(ice)))
     454    case('h2o')
     455        rho = -3.5353e-4*T**2 + 0.0351*T + 933.5030 ! Rottgers 2012
     456    case('co2')
     457        rho = (1.72391 - 2.53e-4*T-2.87*1e-7*T**2)*1.e3 ! Mangan et al. 2017
     458    case default
     459        call abort_pem("rho_ice","Type of ice not known!",1)
     460end select
     461
     462END FUNCTION rho_ice
     463
    436464END MODULE ice_table_mod

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

    r3525 r3527  
    7171use co2condens_mod,             only: CO2cond_ps
    7272use layering_mod,               only: d_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo
     73use dyn_ss_ice_m_mod,           only: dyn_ss_ice_m
    7374
    7475#ifndef CPP_STD
Note: See TracChangeset for help on using the changeset viewer.