Ignore:
Timestamp:
Dec 6, 2023, 4:02:06 PM (12 months ago)
Author:
jbclement
Message:

PEM:

  • Simplification of the algorithm managing the stopping criteria;
  • Complete rework of the ice management in the PEM (H2O & CO2);

    Subroutines to evolve the H2O and CO2 ice are now in the same module "evol_ice_mod.F90".
    Tendencies are computed from the variation of "ice + frost" between the 2 PCM runs.
    Evolving ice in the PEM is now called 'h2o_ice' or 'co2_ice' (not anymore in 'qsurf' and free of 'water_reservoir').
    Default value 'ini_h2o_bigreservoir' (= 10 m) initializes the H2O ice of the first PEM run where there is 'watercap'. For the next PEM runs, initialization is done with the value kept in "startpem.nc". CO2 ice is taken from 'perennial_co2ice' of the PCM (paleoclimate flag must be true).
    Simplification of the condition to compute the surface ice cover needed for the stopping criteria.
    Frost ('qsurf') is not evolved by the PEM and given back to the PCM.
    New default threshold value 'inf_h2oice_threshold' (= 2 m) to decide at the end of the PEM run if the H2O ice should be 'watercap' or not for the next PCM runs. If H2O ice cannot be 'watercap', then the remaining H2O ice is transferred to the frost ('qsurf').

  • Renaming of variables/subroutines for clarity;
  • Some cleanings throughout the code;
  • Small updates in files of the deftank.

JBC

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r3082 r3149  
    1 module glaciers_mod
    2        
    3  implicit none
    4  LOGICAL  co2glaciersflow ! True by default, to compute co2 ice flow. Read in  pem.def
    5  LOGICAL  h2oglaciersflow ! True by default, to compute co2 ice flow. Read in  pem.def
    6 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    7 !!!
    8 !!! Purpose: Compute CO2 glacier flows
    9 !!!
    10 !!! Author: LL
    11 !!!
    12 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    13 
     1MODULE glaciers_mod
     2
     3implicit none
     4
     5logical :: co2_ice_flow ! True by default, to compute co2 ice flow. Read in run_PEM.def
     6logical :: h2o_ice_flow ! True by default, to compute h2o ice flow. Read in run_PEM.def
     7
     8!=======================================================================
    149contains
    15 
    16 subroutine co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,co2ice,flag_co2flow,flag_co2flow_mesh)
     10!=======================================================================
     11
     12SUBROUTINE flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,co2ice,flag_co2flow,flag_co2flow_mesh)
    1713
    1814!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    2622!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2723
    28 IMPLICIT NONE
     24implicit none
    2925
    3026! arguments
    … …  
    3228
    3329! Inputs:
    34       INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
    35       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       REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
    37       REAL,INTENT(IN) :: ps_GCM(ngrid,timelen)      ! Physical x Time field: surface pressure given by the GCM [Pa]
    38       REAL,INTENT(IN) :: global_ave_ps_GCM          ! Global averaged surface pressure from the GCM [Pa]
    39       REAL,INTENT(IN) :: global_ave_ps_PEM          ! global averaged surface pressure during the PEM iteration [Pa]
     30      integer,                           intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
     31      real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Physical points x Slopes: Distribution of the subgrid slopes
     32      real,    dimension(ngrid),         intent(in) :: def_slope_mean                ! Physical points: values of the sub grid slope angles
     33      real,    dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM                   ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
     34      real,    dimension(ngrid,timelen), intent(in) :: ps_PCM                        ! Physical x Time field: surface pressure given by the PCM [Pa]
     35      real,                              intent(in) :: global_avg_ps_PCM             ! Global averaged surface pressure from the PCM [Pa]
     36      real,                              intent(in) :: global_avg_ps_PEM             ! global averaged surface pressure during the PEM iteration [Pa]
    4037     
    4138! Ouputs:
    42       REAL,INTENT(INOUT) :: co2ice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
    43       REAL,INTENT(INOUT) :: flag_co2flow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
    44       REAL,INTENT(INOUT) :: flag_co2flow_mesh(ngrid)  ! same but within the mesh
    45 
     39      real, dimension(ngrid,nslope), intent(inout) :: co2ice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
     40      real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow      ! flag to see if there is flow on the subgrid slopes
     41      real, dimension(ngrid),        intent(inout) :: flag_co2flow_mesh ! same but within the mesh
    4642
    4743! Local
    48       REAL :: Tcond(ngrid,nslope) ! Physical field: CO2 condensation temperature [K]
    49       REAL :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow
     44      real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
     45      real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow
    5046
    5147!-----------------------------
    52       call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond)
    53 
    54       call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
    55 
    56       call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,"co2",co2ice,flag_co2flow,flag_co2flow_mesh)
    57    RETURN   
    58 end subroutine
    59 
    60 
    61 
    62 
    63 
    64 subroutine h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh)
     48write(*,*) "Flow of CO2 glacier"
     49   
     50call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond)
     51call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
     52call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,"co2",co2ice,flag_co2flow,flag_co2flow_mesh)
     53
     54END SUBROUTINE flow_co2glaciers
     55
     56!=======================================================================
     57
     58SUBROUTINE flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh)
    6559
    6660!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    7468!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    7569
    76 
    77 IMPLICIT NONE
     70implicit none
    7871
    7972! arguments
    … …  
    8174
    8275! Inputs:
    83       INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
    84       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
    85       REAL,INTENT(IN) :: Tice(ngrid,nslope) ! Ice Temperature [K]
     76      integer,intent(in) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
     77      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
     78      real,intent(in) :: Tice(ngrid,nslope) ! Ice Temperature [K]
    8679! Ouputs:
    87       REAL,INTENT(INOUT) :: h2oice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
    88       REAL,INTENT(INOUT) :: flag_h2oflow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
    89       REAL,INTENT(INOUT) :: flag_h2oflow_mesh(ngrid)  ! same but within the mesh
     80      real,intent(inout) :: h2oice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
     81      real,intent(inout) :: flag_h2oflow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
     82      real,intent(inout) :: flag_h2oflow_mesh(ngrid)  ! same but within the mesh
    9083! Local
    91       REAL :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow
     84      real :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow
    9285
    9386!-----------------------------
    94 
    95       call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
    96       call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh)
    97 
    98    RETURN
    99 end subroutine
    100 
    101 
    102 subroutine compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
     87write(*,*) "Flow of H2O glaciers"
     88
     89call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
     90call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh)
     91
     92END SUBROUTINE flow_h2oglaciers
     93
     94!=======================================================================
     95
     96SUBROUTINE compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
    10397
    10498use abort_pem_mod, only: abort_pem
    … …  
    118112!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    119113
    120    IMPLICIT NONE
     114implicit none
    121115
    122116! arguments
    … …  
    124118
    125119! Inputs
    126       INTEGER,INTENT(IN) :: ngrid,nslope ! # of grid points and subslopes
    127       INTEGER,INTENT(IN) :: iflat        ! index of the flat subslope
    128       REAL,INTENT(IN) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg]
    129       REAL,INTENT(IN) :: Tice(ngrid,nslope)     ! Physical field:  ice temperature [K]
    130       character(len=3), INTENT(IN) :: name_ice ! Nature of the ice
     120      integer,intent(in) :: ngrid,nslope ! # of grid points and subslopes
     121      integer,intent(in) :: iflat        ! index of the flat subslope
     122      real,intent(in) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg]
     123      real,intent(in) :: Tice(ngrid,nslope)     ! Physical field:  ice temperature [K]
     124      character(len=3), intent(in) :: name_ice ! Nature of the ice
    131125! Outputs
    132       REAL,INTENT(OUT) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum  thickness before flaw [m]
     126      real,intent(out) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum  thickness before flaw [m]
    133127! Local
    134128      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
    135       REAL :: rho(ngrid,nslope) ! co2 ice density [kg/m^3]
    136       INTEGER :: ig,islope ! loop variables
    137       REAL :: slo_angle
     129      real :: rho(ngrid,nslope) ! co2 ice density [kg/m^3]
     130      integer :: ig,islope ! loop variables
     131      real :: slo_angle
    138132
    139133! 1. Compute rho
    … …  
    167161       ENDDO
    168162    ENDDO
    169 RETURN
    170 
    171 end subroutine
    172 
    173 
    174 
    175 
    176 subroutine transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,name_ice,qice,flag_flow,flag_flowmesh)
     163END SUBROUTINE compute_hmaxglaciers
     164
     165!=======================================================================
     166
     167SUBROUTINE transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,name_ice,qice,flag_flow,flag_flowmesh)
    177168!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    178169!!!
    … …  
    196187
    197188! 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
     189      integer, intent(in) :: ngrid,nslope               !# of physical points and subslope
     190      integer, intent(in) :: iflat                      ! index of the flat subslope
     191      real, intent(in) :: subslope_dist(ngrid,nslope)   ! Distribution of the subgrid slopes within the mesh
     192      real, intent(in) :: def_slope_mean(nslope)        ! values of the subgrid slopes
     193      real, intent(in) :: hmax(ngrid,nslope)            ! maximum height of the  glaciers before initiating flow [m]
     194      real, intent(in) :: Tice(ngrid,nslope)            ! Ice temperature[K]
     195      character(len=3), intent(in) :: name_ice              ! Nature of the ice
    205196
    206197! 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
     198      real, intent(inout) :: qice(ngrid,nslope)      ! CO2 in the subslope [kg/m^2]
     199      real, intent(inout) :: flag_flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope
     200      real, intent(inout) :: flag_flowmesh(ngrid)    ! boolean to check if there is flow in the mesh
    210201! 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
     202      integer ig,islope ! loop
     203      real rho(ngrid,nslope) ! density of ice, temperature dependant [kg/m^3]
     204      integer iaval ! ice will be transfered here
    214205
    215206! 0. Compute rho
    … …  
    268259        ENDDO !islope
    269260       ENDDO !ig
    270 RETURN
    271 end subroutine
    272 
    273      subroutine computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond)
     261END SUBROUTINE
     262
     263!=======================================================================
     264
     265SUBROUTINE computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond)
    274266!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    275267!!!
    … …  
    288280
    289281! INPUT
    290       INTEGER,INTENT(IN) :: timelen, ngrid,nslope ! # of timesample,physical points, subslopes
    291       REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical points x times field: VMR of CO2 in the first layer [mol/mol]
    292       REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical points x times field: surface pressure in the GCM [Pa]
    293       REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surfacepressure in the GCM [Pa]
    294       REAL, INTENT(IN) :: global_ave_ps_PEM ! Global averaged surface pressure computed during the PEM iteration
     282integer,                        intent(in) :: timelen, ngrid, nslope ! # of timesample, physical points, subslopes
     283real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM            ! Physical points x times field: VMR of CO2 in the first layer [mol/mol]
     284real, dimension(ngrid,timelen), intent(in) :: ps_PCM                 ! Physical points x times field: surface pressure in the PCM [Pa]
     285real,                           intent(in) :: global_avg_ps_PCM      ! Global averaged surfacepressure in the PCM [Pa]
     286real,                           intent(in) :: global_avg_ps_PEM      ! Global averaged surface pressure computed during the PEM iteration
     287
    295288! OUTPUT
    296       REAL,INTENT(OUT) :: Tcond(ngrid,nslope) ! Physical points : condensation temperature of CO2, yearly averaged
     289real, dimension(ngrid,nslope), intent(out) :: Tcond ! Physical points: condensation temperature of CO2, yearly averaged
    297290
    298291! LOCAL
    299 
    300       INTEGER :: ig,it,islope ! for loop
    301       REAL :: ave ! intermediate to compute average
     292integer :: ig, it ! For loop
     293real    :: ave    ! Intermediate to compute average
    302294
    303295!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    304 
    305 
    306       DO ig = 1,ngrid
    307         ave = 0
    308         DO it = 1,timelen
    309            ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(ig,it)*ps_GCM(ig,it)*global_ave_ps_GCM/global_ave_ps_PEM/100))
    310         ENDDO
    311         DO islope = 1,nslope
    312           Tcond(ig,islope) = ave/timelen
    313         ENDDO
    314       ENDDO
    315 RETURN
    316 
    317 
    318 end subroutine
    319 end module
     296do ig = 1,ngrid
     297    ave = 0
     298    do it = 1,timelen
     299        ave = ave + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM(ig,it)*ps_PCM(ig,it)*global_avg_ps_PCM/global_avg_ps_PEM/100))
     300    enddo
     301    Tcond(ig,:) = ave/timelen
     302enddo
     303
     304END SUBROUTINE computeTcondCO2
     305
     306END MODULE glaciers_mod
Note: See TracChangeset for help on using the changeset viewer.