Changeset 3983 for trunk/LMDZ.COMMON


Ignore:
Timestamp:
Dec 8, 2025, 11:27:43 AM (6 days ago)
Author:
jbclement
Message:

PEM:

  • Removing completely the ice metamorphism computed by a threshold at the end of the PCM (which was commented).
  • Addition of a module "metamorphism" to compute the PCM frost at the PEM beginning and give it back to the PCM at the PEM end. The frost is considered as the ice given by the PCM "startfi.nc" which is above the yearly minimum. Thereby, metamorphism is performed through this operation.
  • Ice reservoirs representation in the PEM is modernized.

JBC

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
9 edited

Legend:

Unmodified
Added
Removed

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

    r3981 r3983  
    812812- Removing the counting method where PCM years are taken into account.
    813813- CO2 ice is now a PEM reservoir such as H2O ice.
     814
     815== 08/12/2025 == JBC
     816- Removing completely the ice metamorphism computed by a threshold at the end of the PCM (which was commented).
     817- Addition of a module "metamorphism" to compute the PCM frost at the PEM beginning and give it back to the PCM at the PEM end. The frost is considered as the ice given by the PCM "startfi.nc" which is above the yearly minimum. Thereby, metamorphism is performed through this operation.
     818- Ice reservoirs representation in the PEM is modernized.

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

    r3980 r3983  
    2727use comsoil_h_pem,         only: soil_pem, fluxgeo, ini_huge_h2oice, depth_breccia, depth_bedrock, reg_thprop_dependp
    2828use adsorption_mod,        only: adsorption_pem
    29 use glaciers_mod,          only: h2oice_flow, co2ice_flow, inf_h2oice_threshold, metam_co2ice_threshold, metam_h2oice_threshold, metam_h2oice, metam_co2ice
     29use glaciers_mod,          only: h2oice_flow, co2ice_flow, inf_h2oice_threshold
    3030use ice_table_mod,         only: icetable_equilibrium, icetable_dynamic
    3131use layering_mod,          only: layering_algo, d_dust, impose_dust_ratio, dust2ice_ratio
    … …  
    166166call getin('inf_h2oice_threshold',inf_h2oice_threshold)
    167167
    168 metam_h2oice = .false.
    169 call getin('metam_h2oice',metam_h2oice)
    170 
    171 metam_h2oice_threshold = 460. ! kg.m-2 (= 0.5 m)
    172 call getin('metam_h2oice_threshold',metam_h2oice_threshold)
    173 
    174168h2oice_flow = .true.
    175169call getin('h2oice_flow',h2oice_flow)
    176 
    177 metam_co2ice = .false.
    178 call getin('metam_co2ice',metam_co2ice)
    179 
    180 metam_co2ice_threshold = 16500. ! kg.m-2 (= 10 m)
    181 call getin('metam_co2ice_threshold',metam_co2ice_threshold)
    182170
    183171co2ice_flow = .true.

  • trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def

    r3926 r3983  
    8383# inf_h2oice_threshold=460.
    8484
    85 # Do you want H2O frost to transform into perennial H2O ice? Default = .false.
    86 # metam_h2oice=.false.
    87 
    88 # Threshold to consider frost is becoming perennial H2O ice. Default = 460. kg.m-2 (= 0.5 m)
    89 # metam_h2oice_threshold=460.
    90 
    9185# Do you want the H2O ice to flow along subslope inside a cell? Default = .true.
    9286# h2oice_flow=.true.
    93 
    94 # Do you want CO2 frost to transform into perennial CO2 ice? Default = .false.
    95 # metam_co2ice=.false.
    96 
    97 # Threshold to consider frost is becoming perennial CO2 ice. Default = 16500. kg.m-2 (= 10 m)
    98 # metam_co2ice_threshold=16500.
    9987
    10088# Do you want the CO2 ice to flow along subslope inside a cell? Default = .true.

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

    r3961 r3983  
    33implicit none
    44
    5 ! Flags for ice management
     5! Different types of ice
     6type :: ice
     7    real :: h2o
     8    real :: co2
     9end type ice
     10
     11! Perennial ice manage by the PEM
     12type(ice), dimension(:,:), allocatable :: perice
     13
     14! Flags for ice flow
    615logical :: h2oice_flow  ! True by default, to compute H2O ice flow. Read in "run_PEM.def"
    716logical :: co2ice_flow  ! True by default, to compute CO2 ice flow. Read in "run_PEM.def"
    8 logical :: metam_h2oice ! False by default, to compute H2O ice metamorphism. Read in "run_PEM.def"
    9 logical :: metam_co2ice ! False by default, to compute CO2 ice metamorphism. Read in "run_PEM.def"
    10 
    11 ! Thresholds for ice management
     17
     18! Threshold to consider H2O ice as watercap
    1219real :: inf_h2oice_threshold   ! To consider the amount of H2O ice as an infinite reservoir [kg.m-2]
    13 real :: metam_h2oice_threshold ! To consider frost is becoming perennial H2O ice [kg.m-2]
    14 real :: metam_co2ice_threshold ! To consider frost is becoming perennial CO2 ice [kg.m-2]
    15 
     20
     21! Ice properties
    1622real, parameter :: rho_co2ice = 1650. ! Density of CO2 ice [kg.m-3]
    1723real, parameter :: rho_h2oice = 920.  ! Density of H2O ice [kg.m-3]
    … …  
    1925!=======================================================================
    2026contains
     27!=======================================================================
     28
     29SUBROUTINE set_perice4PCM(ngrid,nslope,PCMfrost,is_perh2oice,PCMh2oice,PCMco2ice)
     30
     31use metamorphism, only: iPCM_h2ofrost
     32use comslope_mod, only: subslope_dist, def_slope_mean
     33#ifndef CPP_1D
     34    use comconst_mod, only: pi
     35#else
     36    use comcstfi_h,   only: pi
     37#endif
     38
     39implicit none
     40
     41! Arguments
     42!----------
     43integer, intent(in) :: ngrid, nslope
     44real, dimension(:,:,:), intent(inout) :: PCMfrost
     45logical, dimension(ngrid),        intent(out) :: is_perh2oice
     46real,    dimension(ngrid,nslope), intent(out) :: PCMco2ice, PCMh2oice
     47
     48! Local variables
     49!----------------
     50integer :: i
     51
     52! Code
     53!-----
     54write(*,*) '> Reconstructing perennial ic for the PCM'
     55PCMco2ice(:,:) = perice(:,:)%co2
     56PCMh2oice(:,:) = 0. ! Because in the Mars PCM, only the variation of perennial H2O ice is monitored, not the absolute quantity
     57do i = 1,ngrid
     58    ! Is H2O ice still considered as an infinite reservoir for the PCM?
     59    if (sum(perice(i,:)%h2o*subslope_dist(i,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then
     60        ! There is enough ice to be considered as an infinite reservoir
     61        is_perh2oice(i) = .true.
     62    else
     63        ! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
     64        is_perh2oice(i) = .false.
     65        PCMfrost(i,iPCM_h2ofrost,:) = PCMfrost(i,iPCM_h2ofrost,:) + perice(i,:)%h2o
     66        perice(i,:)%h2o = 0.
     67    endif
     68enddo
     69
     70END SUBROUTINE set_perice4PCM
     71!=======================================================================
     72
     73SUBROUTINE ini_ice(ngrid,nslope)
     74
     75implicit none
     76
     77! Arguments
     78!----------
     79integer, intent(in) :: ngrid, nslope
     80
     81! Local variables
     82!----------------
     83
     84! Code
     85!-----
     86if (.not. allocated(perice)) allocate(perice(ngrid,nslope))
     87
     88END SUBROUTINE ini_ice
     89!=======================================================================
     90
     91SUBROUTINE end_ice()
     92
     93implicit none
     94
     95! Arguments
     96!----------
     97
     98! Local variables
     99!----------------
     100
     101! Code
     102!-----
     103if (allocated(perice)) deallocate(perice)
     104
     105END SUBROUTINE end_ice
     106
    21107!=======================================================================
    22108

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

    r3961 r3983  
    667667  USE comsoil_h_PEM, only:  nsoilmx_PEM
    668668  USE comslope_mod, ONLY: nslope
    669   use layering_mod, only: nb_str_max
    670669  USE mod_grid_phy_lmdz
    671670  USE mod_phys_lmdz_para

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

    r3979 r3983  
    3838use conf_pem_mod,               only: conf_pem
    3939use pemredem,                   only: pemdem0, pemdem1
    40 use glaciers_mod,               only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, &
    41                                       metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice, computeTcondCO2
     40use glaciers_mod,               only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, computeTcondCO2, set_perice4PCM
    4241use stopping_crit_mod,          only: stopping_crit_h2o_ice, stopping_crit_co2, stopping_crit_h2o, stopFlags
    4342use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
    … …  
    8382use dimradmars_mod,             only: totcloudfrac, albedo
    8483use dust_param_mod,             only: tauscaling
    85 use tracer_mod,                 only: noms, igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
     84use tracer_mod,                 only: noms, mmol, igcm_h2o_vap ! Tracer names and molar masses
    8685use mod_phys_lmdz_para,         only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
    8786use planete_h,                  only: year_day
    8887use surfini_mod,                only: surfini
    8988use comcstfi_h,                 only: mugaz
     89use metamorphism,               only: ini_frost_id, set_frost4PCM, iPCM_h2ofrost, iPCM_co2frost
    9090
    9191#ifndef CPP_1D
    … …  
    161161real                                  :: ps_avg_global_old    ! constant: Global average pressure of previous time step
    162162real                                  :: ps_avg_global_new    ! constant: Global average pressure of current time step
    163 real,   dimension(:,:),   allocatable :: zplev_new            ! Grid points x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2]
    164 real,   dimension(:,:),   allocatable :: zplev_start0         ! Grid points x Atmospheric field: mass of the atmospheric layers in the start [kg/m^2]
    165 real,   dimension(:,:,:), allocatable :: zplev_timeseries_new ! Grid points x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
    166 real,   dimension(:,:,:), allocatable :: zplev_timeseries_old ! same but with the time series, for previous time step
     163real, dimension(:,:),     allocatable :: zplev_new            ! Grid points x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2]
     164real, dimension(:,:),     allocatable :: zplev_start0         ! Grid points x Atmospheric field: mass of the atmospheric layers in the start [kg/m^2]
     165real, dimension(:,:,:),  allocatable :: zplev_timeseries_new ! Grid points x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
     166real, dimension(:,:,:),  allocatable :: zplev_timeseries_old ! same but with the time series, for previous time step
    167167type(stopFlags)                       :: stopCrit             ! Stopping criteria
    168168real                                  :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
    169 real,   dimension(:,:),   allocatable :: q_h2o_PEM_phys       ! Grid points x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
     169real, dimension(:,:),     allocatable :: q_h2o_PEM_phys       ! Grid points x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
    170170integer                               :: timelen              ! # time samples
    171171real                                  :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
    172 real, dimension(:,:),    allocatable  :: d_h2oice_new         ! Adjusted tendencies to keep balance between donor and recipient
     172real, dimension(:,:),     allocatable :: d_h2oice_new         ! Adjusted tendencies to keep balance between donor and recipient
     173real, dimension(:,:),     allocatable :: min_h2oice           ! Yearly minimum of H2O ice for the last year of the PCM
    173174real                                  :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to balance H2O
    174175
    … …  
    183184real,    dimension(:,:), allocatable :: vmr_co2_PEM_phys       ! Grid points x Times co2 volume mixing ratio used in the PEM
    184185real,    dimension(:,:), allocatable :: q_co2_PEM_phys         ! Grid points x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
     186real,    dimension(:,:), allocatable :: min_co2ice             ! Yearly minimum of CO2 ice for the last year of the PCM
    185187real(kind = 16)                      :: totmass_co2ice, totmass_atmco2         ! Current total CO2 masses
    186188real(kind = 16)                      :: totmass_co2ice_ini, totmass_atmco2_ini ! Initial total CO2 masses
    … …  
    413415
    414416call surfini(ngrid,nslope,qsurf)
    415 
    416 do nnq = 1,nqtot  ! Why not using ini_tracer?
    417     if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
     417call ini_frost_id(nqtot,noms)
     418do nnq = 1,nqtot
    418419    if (noms(nnq) == "h2o_vap") then
    419420        igcm_h2o_vap = nnq
    420421        mmol(igcm_h2o_vap) = 18.
    421422    endif
    422     if (noms(nnq) == "co2") igcm_co2 = nnq
    423423enddo
    424424
    … …  
    448448allocate(watersurf_density_avg(ngrid,nslope),watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
    449449allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope),d_co2ice_ini(ngrid,nslope))
    450 
    451 call read_PCM_data(ngrid,nslope,nsoilmx,timelen,ps_avg,tsurf_avg,tsurf_avg_yr1,tsoil_avg,tsoil_timeseries,watersurf_density_avg,d_h2oice,d_co2ice,ps_timeseries,q_h2o_PEM_phys,q_co2_PEM_phys,watersoil_density_timeseries)
     450allocate(min_co2ice(ngrid,nslope),min_h2oice(ngrid,nslope))
     451
     452call read_PCM_data(ngrid,nslope,nsoilmx,timelen,qsurf(:,iPCM_h2ofrost,:),qsurf(:,iPCM_co2frost,:),ps_avg,tsurf_avg,tsurf_avg_yr1,tsoil_avg,tsoil_timeseries,watersurf_density_avg, &
     453                   d_h2oice,d_co2ice,ps_timeseries,q_h2o_PEM_phys,q_co2_PEM_phys,watersoil_density_timeseries,min_h2oice,min_co2ice)
     454
    452455vmr_co2_PCM = q_co2_PEM_phys/(A*q_co2_PEM_phys + B)/m_co2
    453456d_co2ice_ini = d_co2ice
    … …  
    947950write(*,*) '********* PEM finalization *********'
    948951! III_a.1 Ice update for start file
    949 write(*,*) '> Reconstructing perennial ice and frost for the PCM'
    950 watercap = 0.
    951 perennial_co2ice = co2_ice
    952 do ig = 1,ngrid
    953     ! H2O ice metamorphism
    954     !if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then
    955     !    h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold
    956     !    qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold
    957     !endif
    958 
    959     ! Is H2O ice still considered as an infinite reservoir for the PCM?
    960     if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then
    961         ! There is enough ice to be considered as an infinite reservoir
    962         watercaptag(ig) = .true.
    963     else
    964         ! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
    965         watercaptag(ig) = .false.
    966         qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
    967         h2o_ice(ig,:) = 0.
    968     endif
    969 
    970     ! CO2 ice metamorphism
    971     !if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then
    972     !    perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold
    973     !    qsurf(ig,igcm_co2,:) = metam_co2ice_threshold
    974     !endif
    975 enddo
     952call set_frost4PCM(qsurf)
     953call set_perice4PCM(ngrid,nslope,qsurf,watercaptag,watercap,perennial_co2ice)
    976954
    977955! III.a.3. Tsurf update for start file
    … …  
    10691047        endif
    10701048        ! H2O frost
    1071         if (qsurf(ig,igcm_h2o_ice,islope) > 0.) then
     1049        if (qsurf(ig,iPCM_h2ofrost,islope) > 0.) then
    10721050            albedo(ig,:,islope) = albedo_h2o_frost
    10731051            emis(ig,islope) = 1.
    10741052        endif
    10751053        ! CO2 frost
    1076         if (qsurf(ig,igcm_co2,islope) > 0.) then
     1054        if (qsurf(ig,iPCM_co2frost,islope) > 0.) then
    10771055            albedo(ig,:,islope) = albedice(icap)
    10781056            emis(ig,islope) = emisice(icap)
    … …  
    11341112call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
    11351113call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
    1136              co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,layerings_map)
     1114             co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,co2_ice,layerings_map)
    11371115
    11381116call info_PEM(i_myear_leg,stopCrit%stop_code(),i_myear,n_myear)

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

    r3981 r3983  
    2525use glaciers_mod,               only: rho_co2ice, rho_h2oice
    2626use comcstfi_h,                 only: r, mugaz, pi
    27 use surfdat_h,                  only: watercaptag, perennial_co2ice
     27use surfdat_h,                  only: watercaptag, perennial_co2ice, qsurf
     28use metamorphism,               only: frost4PCM, iPCM_h2ofrost, iPCM_co2frost
    2829
    2930implicit none
    … …  
    8081!!!        4. Mass of CO2 & H2O adsorbed
    8182!!!
    82 !!! /!\ This order must be respected !
     83!!! /!\ This order must be respected!
    8384!!! Author: LL
    8485!!!
    … …  
    102103
    103104!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    104     ! h2o ice
    105     h2o_ice = 0.
     105    ! H2O ice
     106    h2o_ice(:,:) = 0.
    106107    call get_field("h2o_ice",h2o_ice,found)
    107108    if (.not. found) then
    … …  
    110111        write(*,*)'with default value ''ini_huge_h2oice'''
    111112        do ig = 1,ngrid
    112             if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
     113            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice + qsurf(ig,iPCM_h2ofrost,:) - frost4PCM(ig,:)%h2o
    113114        enddo
    114115    endif
    115116
    116     ! co2 ice
    117     co2_ice = 0.
     117    ! CO2 ice
     118    co2_ice(:,:) = 0.
    118119    call get_field("co2_ice",co2_ice,found)
    119120    if (.not. found) then
    120121        write(*,*)'Pemetat0: failed loading <co2_ice>'
    121122        write(*,*)'will reconstruct the values from perennial_co2ice'
    122         co2_ice = perennial_co2ice
     123        co2_ice(:,:) = perennial_co2ice(:,:) + qsurf(:,iPCM_co2frost,:) - frost4PCM(:,:)%co2
    123124    endif
    124125
    … …  
    208209                        delta = depth_breccia
    209210                        TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
    210                                                                   (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
    211                                                                   ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
    212                         do iloop=index_breccia + 2,index_bedrock
     211                                                                   (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
     212                                                                   ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
     213                        do iloop = index_breccia + 2,index_bedrock
    213214                            TI_PEM(ig,iloop,islope) = TI_breccia
    214215                        enddo
    … …  
    246247            do ig = 1,ngrid
    247248                if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
    248                     inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
    249                                                               (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
    250                                                               ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
     249                    inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
     250                                                                (((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
     251                                                                ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
    251252
    252253                    do iloop = index_breccia + 2,index_bedrock
    … …  
    254255                    enddo
    255256                else
    256                     do iloop=index_breccia + 1,index_bedrock
     257                    do iloop = index_breccia + 1,index_bedrock
    257258                        inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
    258259                    enddo
    … …  
    263264            delta = depth_bedrock
    264265            do ig = 1,ngrid
    265                 inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
    266                                                            (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
    267                                                            ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
    268             enddo
    269 
    270             do iloop = index_bedrock + 2, nsoil_PEM
     266                inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
     267                                                            (((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
     268                                                            ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
     269            enddo
     270
     271            do iloop = index_bedrock + 2,nsoil_PEM
    271272                do ig = 1,ngrid
    272273                    inertiedat_PEM(ig,iloop) = TI_bedrock
    … …  
    365366
    366367    ! h2o ice
    367     h2o_ice = 0.
     368    h2o_ice(:,:) = 0.
    368369    write(*,*)'So ''h2o_ice'' is initialized with default value ''ini_huge_h2oice'' where ''watercaptag'' is true.'
    369370    do ig = 1,ngrid
    370         if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
     371        if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice + qsurf(ig,iPCM_h2ofrost,:) - frost4PCM(ig,:)%h2o
    371372    enddo
    372373
    373374    ! co2 ice
    374375    write(*,*)'So ''co2_ice'' is initialized with ''perennial_co2ice'' found in the PCM.'
    375     co2_ice = perennial_co2ice
     376    co2_ice(:,:) = perennial_co2ice(:,:) + qsurf(:,iPCM_co2frost,:) - frost4PCM(:,:)%co2
    376377
    377378    ! Layerings

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

    r3981 r3983  
    5858
    5959SUBROUTINE pemdem1(filename,i_myear,nsoil_PEM,ngrid,nslope,tsoil_slope_PEM,inertiesoil_slope_PEM, &
    60                    icetable_depth,icetable_thickness,ice_porefilling,m_co2_regolith,m_h2o_regolith,h2o_ice,layerings_map)
     60                   icetable_depth,icetable_thickness,ice_porefilling,m_co2_regolith,m_h2o_regolith,h2o_ice,co2_ice,layerings_map)
    6161
    6262! write time-dependent variable to restart file

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

    r3977 r3983  
    1919contains
    2020!=======================================================================
    21 SUBROUTINE read_PCM_data(ngrid,nslope,nsoil_PCM,nsol,ps_avg,tsurf_avg,tsurf_avg_y1,tsoil_avg,tsoil_ts,watersurf_density_avg,d_h2oice,d_co2ice, &
    22                          ps_ts,q_h2o_ts,q_co2_ts,watersoil_density_ts)
    23 
    24 use compute_tend_mod, only: compute_tend
     21SUBROUTINE read_PCM_data(ngrid,nslope,nsoil_PCM,nsol,h2ofrost_PCM,co2frost_PCM,ps_avg,tsurf_avg,tsurf_avg_y1,tsoil_avg,tsoil_ts,watersurf_density_avg,d_h2oice,d_co2ice, &
     22                         ps_ts,q_h2o_ts,q_co2_ts,watersoil_density_ts,min_h2oice,min_co2ice)
     23
    2524use grid_conversion,  only: lonlat2vect
    2625use comsoil_h_PEM,    only: soil_pem
     26use compute_tend_mod, only: compute_tend
     27use metamorphism,     only: compute_frost
    2728
    2829implicit none
    … …  
    3031! Arguments
    3132!----------
    32 integer, intent(in) :: ngrid, nslope, nsoil_PCM, nsol
     33integer,                       intent(in) :: ngrid, nslope, nsoil_PCM, nsol
     34real, dimension(ngrid,nslope), intent(in) :: h2ofrost_PCM, co2frost_PCM
    3335real, dimension(ngrid),                       intent(out) :: ps_avg
    34 real, dimension(ngrid,nslope),                intent(out) :: tsurf_avg, tsurf_avg_y1, watersurf_density_avg, d_h2oice, d_co2ice
     36real, dimension(ngrid,nslope),                intent(out) :: tsurf_avg, tsurf_avg_y1, watersurf_density_avg, d_h2oice, d_co2ice, min_h2oice, min_co2ice
    3537real, dimension(ngrid,nsoil_PCM,nslope),      intent(out) :: tsoil_avg
    3638real, dimension(ngrid,nsol),                  intent(out) :: ps_ts, q_h2o_ts, q_co2_ts
    … …  
    4446real, dimension(:,:,:,:), allocatable :: var_read_4d
    4547character(:),             allocatable :: num ! For reading slope variables
    46 real, dimension(ngrid,nslope,2)       :: min_h2oice, min_co2ice, min_h2ofrost, min_co2frost
     48real, dimension(ngrid,nslope,2)       :: min_h2operice, min_co2perice, min_h2ofrost, min_co2frost
    4749
    4850! Code
    4951!-----
    5052! Initialization
    51 min_h2oice = 0.
    52 min_co2ice = 0.
     53min_h2operice = 0.
     54min_co2perice = 0.
    5355min_h2ofrost = 0.
    5456min_co2frost = 0.
    … …  
    8082    call get_var('co2ice'//num,var_read_2d)          ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2frost(:,islope,1))
    8183    call get_var('h2o_ice_s'//num,var_read_2d)       ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2ofrost(:,islope,1))
    82     call get_var('watercap'//num,var_read_2d)        ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2oice(:,islope,1))
    83     call get_var('perennial_co2ice'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2ice(:,islope,1))
     84    call get_var('watercap'//num,var_read_2d)        ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2operice(:,islope,1))
     85    call get_var('perennial_co2ice'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2perice(:,islope,1))
    8486    call get_var('tsurf'//num,var_read_2d)           ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,tsurf_avg_y1(:,islope))
    8587enddo
    … …  
    9597
    9698! Allocate and read the variables
    97 call get_var('ps',var_read_2d) ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,ps_avg)
     99call get_var('ps',var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,ps_avg)
    98100do islope = 1,nslope
    99101    if (nslope /= 1) then
    … …  
    105107    call get_var('co2ice'//num,var_read_2d)          ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2frost(:,islope,2))
    106108    call get_var('h2o_ice_s'//num,var_read_2d)       ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2ofrost(:,islope,2))
    107     call get_var('watercap'//num,var_read_2d)        ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2oice(:,islope,2))
    108     call get_var('perennial_co2ice'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2ice(:,islope,2))
     109    call get_var('watercap'//num,var_read_2d)        ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2operice(:,islope,2))
     110    call get_var('perennial_co2ice'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2perice(:,islope,2))
    109111    if (soil_pem) then
    110112        call get_var('soiltemp'//num,var_read_3d)
    … …  
    167169
    168170!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     171! Compute frost from yearly minima
     172call compute_frost(ngrid,nslope,h2ofrost_PCM,min_h2ofrost,co2frost_PCM,min_co2frost)
     173
     174!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    169175! Compute ice tendencies from yearly minima
    170176write(*,*) '> Computing surface ice tendencies'
    171 call compute_tend(ngrid,nslope,min_h2ofrost + min_h2oice,d_h2oice)
     177call compute_tend(ngrid,nslope,min_h2operice + min_h2ofrost,d_h2oice)
    172178write(*,*) 'H2O ice tendencies (min/max):', minval(d_h2oice), maxval(d_h2oice)
    173 call compute_tend(ngrid,nslope,min_co2frost + min_co2ice,d_co2ice)
     179call compute_tend(ngrid,nslope,min_co2perice + min_co2frost,d_co2ice)
    174180write(*,*) 'CO2 ice tendencies (min/max):', minval(d_co2ice), maxval(d_co2ice)
    175181
Note: See TracChangeset for help on using the changeset viewer.