Changeset 3983 for trunk/LMDZ.COMMON/libf/evolution/pem.F90

Timestamp:

Dec 8, 2025, 11:27:43 AM (7 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

File:

• trunk/LMDZ.COMMON/libf/evolution/pem.F90 (modified) (9 diffs)

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

@@ -38 +38 @@
 use conf_pem_mod,               only: conf_pem
 use pemredem,                   only: pemdem0, pemdem1
-use glaciers_mod,               only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, &
-                                      metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice, computeTcondCO2
+use glaciers_mod,               only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, computeTcondCO2, set_perice4PCM
 use stopping_crit_mod,          only: stopping_crit_h2o_ice, stopping_crit_co2, stopping_crit_h2o, stopFlags
 use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
@@ -83 +82 @@
 use dimradmars_mod,             only: totcloudfrac, albedo
 use dust_param_mod,             only: tauscaling
-use tracer_mod,                 only: noms, igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
+use tracer_mod,                 only: noms, mmol, igcm_h2o_vap ! Tracer names and molar masses
 use mod_phys_lmdz_para,         only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
 use planete_h,                  only: year_day
 use surfini_mod,                only: surfini
 use comcstfi_h,                 only: mugaz
+use metamorphism,               only: ini_frost_id, set_frost4PCM, iPCM_h2ofrost, iPCM_co2frost
 
 #ifndef CPP_1D
@@ -161 +161 @@
 real                                  :: ps_avg_global_old    ! constant: Global average pressure of previous time step
 real                                  :: ps_avg_global_new    ! constant: Global average pressure of current time step
-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]
-real,   dimension(:,:),   allocatable :: zplev_start0         ! Grid points x Atmospheric field: mass of the atmospheric layers in the start [kg/m^2]
-real,   dimension(:,:,:), allocatable :: zplev_timeseries_new ! Grid points x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
-real,   dimension(:,:,:), allocatable :: zplev_timeseries_old ! same but with the time series, for previous time step
+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]
+real, dimension(:,:),     allocatable :: zplev_start0         ! Grid points x Atmospheric field: mass of the atmospheric layers in the start [kg/m^2]
+real, dimension(:,:,:),   allocatable :: zplev_timeseries_new ! Grid points x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
+real, dimension(:,:,:),   allocatable :: zplev_timeseries_old ! same but with the time series, for previous time step
 type(stopFlags)                       :: stopCrit             ! Stopping criteria
 real                                  :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
-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]
+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]
 integer                               :: timelen              ! # time samples
 real                                  :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
-real, dimension(:,:),    allocatable  :: d_h2oice_new         ! Adjusted tendencies to keep balance between donor and recipient
+real, dimension(:,:),     allocatable :: d_h2oice_new         ! Adjusted tendencies to keep balance between donor and recipient
+real, dimension(:,:),     allocatable :: min_h2oice           ! Yearly minimum of H2O ice for the last year of the PCM
 real                                  :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to balance H2O
 
@@ -183 +184 @@
 real,    dimension(:,:), allocatable :: vmr_co2_PEM_phys       ! Grid points x Times co2 volume mixing ratio used in the PEM
 real,    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]
+real,    dimension(:,:), allocatable :: min_co2ice             ! Yearly minimum of CO2 ice for the last year of the PCM
 real(kind = 16)                      :: totmass_co2ice, totmass_atmco2         ! Current total CO2 masses
 real(kind = 16)                      :: totmass_co2ice_ini, totmass_atmco2_ini ! Initial total CO2 masses
@@ -413 +415 @@
 
 call surfini(ngrid,nslope,qsurf)
-
-do nnq = 1,nqtot  ! Why not using ini_tracer?
-    if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
+call ini_frost_id(nqtot,noms)
+do nnq = 1,nqtot
     if (noms(nnq) == "h2o_vap") then
         igcm_h2o_vap = nnq
         mmol(igcm_h2o_vap) = 18.
     endif
-    if (noms(nnq) == "co2") igcm_co2 = nnq
 enddo
 
@@ -448 +448 @@
 allocate(watersurf_density_avg(ngrid,nslope),watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
 allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope),d_co2ice_ini(ngrid,nslope))
-
-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)
+allocate(min_co2ice(ngrid,nslope),min_h2oice(ngrid,nslope))
+
+call 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, &
+                   d_h2oice,d_co2ice,ps_timeseries,q_h2o_PEM_phys,q_co2_PEM_phys,watersoil_density_timeseries,min_h2oice,min_co2ice)
+
 vmr_co2_PCM = q_co2_PEM_phys/(A*q_co2_PEM_phys + B)/m_co2
 d_co2ice_ini = d_co2ice
@@ -947 +950 @@
 write(*,*) '********* PEM finalization *********'
 ! III_a.1 Ice update for start file
-write(*,*) '> Reconstructing perennial ice and frost for the PCM'
-watercap = 0.
-perennial_co2ice = co2_ice
-do ig = 1,ngrid
-    ! H2O ice metamorphism
-    !if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then
-    !    h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold
-    !    qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold
-    !endif
-
-    ! Is H2O ice still considered as an infinite reservoir for the PCM?
-    if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then
-        ! There is enough ice to be considered as an infinite reservoir
-        watercaptag(ig) = .true.
-    else
-        ! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
-        watercaptag(ig) = .false.
-        qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
-        h2o_ice(ig,:) = 0.
-    endif
-
-    ! CO2 ice metamorphism
-    !if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then
-    !    perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold
-    !    qsurf(ig,igcm_co2,:) = metam_co2ice_threshold
-    !endif
-enddo
+call set_frost4PCM(qsurf)
+call set_perice4PCM(ngrid,nslope,qsurf,watercaptag,watercap,perennial_co2ice)
 
 ! III.a.3. Tsurf update for start file
@@ -1069 +1047 @@
         endif
         ! H2O frost
-        if (qsurf(ig,igcm_h2o_ice,islope) > 0.) then
+        if (qsurf(ig,iPCM_h2ofrost,islope) > 0.) then
             albedo(ig,:,islope) = albedo_h2o_frost
             emis(ig,islope) = 1.
         endif
         ! CO2 frost
-        if (qsurf(ig,igcm_co2,islope) > 0.) then
+        if (qsurf(ig,iPCM_co2frost,islope) > 0.) then
             albedo(ig,:,islope) = albedice(icap)
             emis(ig,islope) = emisice(icap)
@@ -1134 +1112 @@
 call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
 call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
-             co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,layerings_map)
+             co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,co2_ice,layerings_map)
 
 call info_PEM(i_myear_leg,stopCrit%stop_code(),i_myear,n_myear)