Changeset 2993


Ignore:
Timestamp:
Jul 11, 2023, 7:25:48 PM (18 months ago)
Author:
llange
Message:

RS PEM
Update CO2 glacier flow law after discussions with I.Smith. The model
used before from Nye et al. was false as it overestimates CO2 ice
thickness required to start the flow
The algorithm now used an easier approach with outputs from Smith et al., JGR Planets 2022. See the full description in the note I'll send by
mail.
LL

~

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2985 r2993  
    1       module co2glaciers_mod
    2         implicit none
    3         LOGICAL  co2glaciersflow ! True by default, to compute co2 ice flow. Read in  pem.def
     1module co2glaciers_mod
     2       
     3 implicit none
     4 LOGICAL  co2glaciersflow ! True by default, to compute co2 ice flow. Read in  pem.def
    45
    56!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    3334
    3435! Inputs:
    35       INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat ! # number of time sample, physical points, subslopes, index of the flat subslope
    36       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
     36      INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat ! number of time sample, physical points, subslopes, index of the flat subslope
     37      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
    3738      REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
    3839      REAL,INTENT(IN) :: ps_GCM(ngrid,timelen)      ! Physical x Time field: surface pressure given by the GCM [Pa]
    … …  
    5556      call compute_hmaxglaciers_co2(ngrid,nslope,iflat,Tcond,def_slope_mean,hmax)
    5657
    57       call  transfer_co2ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,co2ice_slope,flag_co2flow,flag_co2flow_mesh)
     58      call transfer_co2ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,co2ice_slope,flag_co2flow,flag_co2flow_mesh)
    5859   RETURN
    5960end subroutine
    … …  
    7071!!!          before initating flow
    7172!!!         
    72 !!! Author: LL, based on theoretical work by A.Grau Galofre (LPG)
     73!!! Author: LL,based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
    7374!!!
    7475!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    8788      REAL,INTENT(OUT) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum co2 thickness before flaw [m]
    8889! Local
    89       INTEGER,PARAMETER :: n = 7 ! flow law exponent Nye et al., 2000
    90       REAL,PARAMETER :: Rg = 8.3145 ! gas constant [J/K/mol]
    91       REAL,PARAMETER :: Q = 59000. ! Activation Energy [J/mol], Nye et al., 2000
    92       DOUBLE PRECISION,PARAMETER :: C = 1.8138e-21 ! Nye et al., 2000 [s/m^n]
    93       DOUBLE PRECISION :: Ad = 1.202e11 ! Softness prefactor [MPa^-n] Nye et al., 2000
    94       REAL :: Ro,Ho, S,ratio ! gemoetry from Nye et al., 2000
    95       DOUBLE PRECISION :: A,Ao ! softness parameter [s/m^n]
    96       DOUBLE PRECISION :: C1 ! intermediate variable
    97       DOUBLE PRECISION :: t_0 ! relaxation time (assuming radial symetry) [s]
    98       DOUBLE PRECISION :: u ! characteristic horizontal deformation rate [m/s]
    99       DOUBLE PRECISION :: tau_d ! characteristic basal drag, understood as the strest that an ice CO2 mass flowing under its weight balanced by viscosity
     90      DOUBLE PRECISION :: tau_d = 5.e3 ! characteristic basal drag, understood as the stress that an ice CO2 mass flowing under its weight balanced by viscosity. Value obtained from I.Smith
    10091      REAL :: rho_co2(ngrid) ! co2 ice density [kg/m^3]
    10192      INTEGER :: ig,islope ! loop variables
    10293      REAL :: slo_angle
    10394
    104 ! 0. Geometry parameters
    105       Ro = 200e3
    106       Ho = 3000.
    107       ratio = 2./3.
    108       S = Ho/Ro*1/((2+1./n)*(1+1./n))*(1-ratio**(1+1./n))**(1./(2+1./n)-1)*ratio**(1+1./n-1)
    109 ! 1. Flow parameters
    110      Ao = 3**(1./(2*n+2))*Ad
    111      do ig = 1,n
    112      Ao = Ao*1e-6
    113      enddo
    114 ! 2. Compute rho_co2
     95! 1. Compute rho_co2
    11596    DO ig = 1,ngrid
    11697      rho_co2(ig) = (1.72391 - 2.53e-4*Tcond(ig)-2.87*1e-7*Tcond(ig)**2)*1e3  ! Mangan et al. 2017
    11798    ENDDO
     99
    118100! 3. Compute max thickness
    119101    DO ig = 1,ngrid
    120        A = Ao*exp(-Q/(Rg*Tcond(ig)))
    121        C1 = A*(rho_co2(ig)*g)**(float(n))/float(n+2)
    122        t_0 = Ro/(C1*(5*n+3))*(float(2*n+1)/float(n+1))**n
    123        u = Ro/t_0
    124        tau_d = (u*Ho*(rho_co2(ig)*g)**2*S**2/(2*A))**(1./(n+2))     
    125102       DO islope = 1,nslope
    126103          if(islope.eq.iflat) then
    … …  
    256233        ENDDO
    257234        Tcond(ig) = ave/timelen
    258        ENDDO
     235      ENDDO
    259236RETURN
    260237
    261238
    262       end subroutine
    263       end module
     239end subroutine
     240end module
Note: See TracChangeset for help on using the changeset viewer.