Changeset 3330 for trunk/LMDZ.COMMON

Timestamp:

May 16, 2024, 11:38:21 AM (6 months ago)

Author:

jbclement

Message:

PEM:
Improvement of a flag-like variable (more robust as a logical) to know where co2 ice was at the beginning.
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• pem.F90 (modified) (12 diffs)

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

@@ -313 +313 @@
 - Cleaning the previous cleaning commit.
 
+== 16/05/2024 == JBC
+Improvement of a flag-like variable (more robust as a logical) to know where co2 ice was at the beginning.

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

@@ -7 +7 @@
 !    I_e Initialization of the PEM variable and soil
 !    I_f Compute tendencies
-!    I_g Save initial PCM situation
-!    I_h Read the "startpem.nc"
+!    I_g Read the "startpem.nc"
+!    I_h Save the initial situation
 !    I_i Compute orbit criterion
 
@@ -45 +45 @@
 use evol_ice_mod,               only: evol_co2_ice, evol_h2o_ice
 use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
-                                      TI_PEM,                           & ! soil thermal inertia
-                                      tsoil_PEM, layer_PEM, &             ! Soil temp, number of subsurface layers, soil mid layer depths
-                                      fluxgeo                             ! Geothermal flux for the PEM and PCM
+                                      TI_PEM,               & ! soil thermal inertia
+                                      tsoil_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
+                                      fluxgeo                 ! Geothermal flux for the PEM and PCM
 use adsorption_mod,             only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
                                       ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
@@ -178 +178 @@
 
 ! Variables for co2_ice evolution
-real,  dimension(:,:),   allocatable :: co2_ice           ! co2 ice in the PEM
-real,  dimension(:,:),   allocatable :: co2_ice_ini       ! Initial amount of co2 ice in the PEM
+real,    dimension(:,:), allocatable :: co2_ice           ! co2 ice in the PEM
+logical, dimension(:,:), allocatable :: is_co2ice_ini     ! Was there co2 ice initially in the PEM?
 real,  dimension(:,:,:), allocatable :: min_co2_ice       ! Minimum of co2 ice at each point for the first year [kg/m^2]
 real                                 :: co2ice_ini_surf   ! Initial surface of sublimating co2 ice
 logical, dimension(:,:), allocatable :: ini_co2ice_sublim ! Logical array to know if co2 ice is sublimating
-real,   dimension(:,:),  allocatable :: vmr_co2_PCM       ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
-real,   dimension(:,:),  allocatable :: vmr_co2_PEM_phys  ! Physics x Times co2 volume mixing ratio used in the PEM
-real,   dimension(:,:),  allocatable :: q_co2_PEM_phys    ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
+real,    dimension(:,:), allocatable :: vmr_co2_PCM       ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
+real,    dimension(:,:), allocatable :: vmr_co2_PEM_phys  ! Physics x Times co2 volume mixing ratio used in the PEM
+real,    dimension(:,:), allocatable :: q_co2_PEM_phys    ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
 
 ! Variables for stratification (layering) evolution
@@ -223 +223 @@
 real                                  :: totmassh2o_adsorbded               ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
 logical                               :: bool_sublim                        ! logical to check if there is sublimation or not
-logical, dimension(:,:),  allocatable :: co2ice_has_disappeared             ! logical to check if a co2ice reservoir has already disappeared at a previous timestep
+logical, dimension(:,:),  allocatable :: co2ice_disappeared                 ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
 real, dimension(:,:),     allocatable :: porefillingice_thickness_prev_iter ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
 real, dimension(:),       allocatable :: delta_h2o_icetablesublim           ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
@@ -497 +497 @@
 allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
 allocate(delta_co2_adsorbded(ngrid))
-allocate(co2ice_has_disappeared(ngrid,nslope))
+allocate(co2ice_disappeared(ngrid,nslope))
 allocate(porefillingice_thickness_prev_iter(ngrid,nslope))
 allocate(delta_h2o_icetablesublim(ngrid))
@@ -556 +556 @@
 !------------------------
 ! I   Initialization
-!    I_g Save initial PCM situation
-!------------------------
-! We save the places where h2o ice is sublimating
-! We compute the surface of h2o ice sublimating
-allocate(ini_co2ice_sublim(ngrid,nslope),ini_h2oice_sublim(ngrid,nslope))
-co2ice_ini_surf = 0.
-h2oice_ini_surf = 0.
-Total_surface = 0.
-ini_co2ice_sublim = .false.
-ini_h2oice_sublim = .false.
-co2ice_has_disappeared(:,:) = .false.
-do i = 1,ngrid
-    Total_surface = Total_surface + cell_area(i)
-    do islope = 1,nslope
-        if (tend_co2_ice(i,islope) < 0.) then
-            ini_co2ice_sublim(i,islope) = .true.
-            co2ice_ini_surf = co2ice_ini_surf + cell_area(i)*subslope_dist(i,islope)
-        endif
-        if (tend_h2o_ice(i,islope) < 0.) then
-            ini_h2oice_sublim(i,islope) = .true.
-            h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
-        endif
-    enddo
-enddo
-
-write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2ice_ini_surf
-write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2oice_ini_surf
-write(*,*) "Total surface of the planet =", Total_surface
-allocate(zplev_PCM(ngrid,nlayer + 1))
-
-do ig = 1,ngrid
-    zplev_PCM(ig,:) = ap + bp*ps_start_PCM(ig)
-enddo
-
-global_avg_press_old = sum(cell_area*ps_start_PCM)/Total_surface
-global_avg_press_PCM = global_avg_press_old
-global_avg_press_new = global_avg_press_old
-write(*,*) "Initial global average pressure =", global_avg_press_PCM
-
-!------------------------
-! I   Initialization
-!    I_h Read the "startpem.nc"
+!    I_g Read the "startpem.nc"
 !------------------------
 allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope))
@@ -616 +575 @@
               watersurf_density_ave,watersoil_density_PEM_ave,co2_adsorbded_phys,delta_co2_adsorbded,                &
               h2o_adsorbded_phys,delta_h2o_adsorbded,stratif)
-
-allocate(co2_ice_ini(ngrid,nslope))
-co2_ice_ini = co2_ice
 
 delta_h2o_icetablesublim = 0.
@@ -649 +605 @@
 !------------------------
 ! I   Initialization
+!    I_h Save the initial situation
+!------------------------
+! We save the places where h2o ice is sublimating
+! We compute the surface of h2o ice sublimating
+allocate(ini_co2ice_sublim(ngrid,nslope),ini_h2oice_sublim(ngrid,nslope),is_co2ice_ini(ngrid,nslope))
+co2ice_ini_surf = 0.
+h2oice_ini_surf = 0.
+Total_surface = 0.
+ini_co2ice_sublim = .false.
+ini_h2oice_sublim = .false.
+is_co2ice_ini = .false.
+co2ice_disappeared = .false.
+do i = 1,ngrid
+    Total_surface = Total_surface + cell_area(i)
+    do islope = 1,nslope
+        if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true.
+        if (tend_co2_ice(i,islope) < 0.) then
+            ini_co2ice_sublim(i,islope) = .true.
+            co2ice_ini_surf = co2ice_ini_surf + cell_area(i)*subslope_dist(i,islope)
+        endif
+        if (tend_h2o_ice(i,islope) < 0.) then
+            ini_h2oice_sublim(i,islope) = .true.
+            h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
+        endif
+    enddo
+enddo
+
+write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2ice_ini_surf
+write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2oice_ini_surf
+write(*,*) "Total surface of the planet =", Total_surface
+allocate(zplev_PCM(ngrid,nlayer + 1))
+
+do ig = 1,ngrid
+    zplev_PCM(ig,:) = ap + bp*ps_start_PCM(ig)
+enddo
+
+global_avg_press_old = sum(cell_area*ps_start_PCM)/Total_surface
+global_avg_press_PCM = global_avg_press_old
+global_avg_press_new = global_avg_press_old
+write(*,*) "Initial global average pressure =", global_avg_press_PCM
+
+!------------------------
+! I   Initialization
 !    I_i Compute orbit criterion
 !------------------------
@@ -792 +791 @@
     do ig = 1,ngrid
         do islope = 1,nslope
-            if (co2_ice_ini(ig,islope) > 0. .and. co2_ice(ig,islope) < 1.e-10 .and. (.not.(co2ice_has_disappeared(ig,islope)))) then ! co2 ice disappeared, look for closest point without co2ice
-                co2ice_has_disappeared(ig,islope) = .true.
+            if (is_co2ice_ini(ig,islope) .and. co2_ice(ig,islope) < 1.e-10 .and. .not. co2ice_disappeared(ig,islope)) then ! co2 ice disappeared, look for closest point without co2ice
+                co2ice_disappeared = .true.
                 if (latitude_deg(ig) > 0) then
                     do ig_loop = ig,ngrid
                         do islope_loop = islope,iflat,-1
-                            if (co2_ice_ini(ig_loop,islope_loop) < 1.e-10 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
+                            if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
                                 tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
                                 bool_sublim = .true.
@@ -808 +807 @@
                     do ig_loop = ig,1,-1
                         do islope_loop = islope,iflat
-                            if (co2_ice_ini(ig_loop,islope_loop) < 1.e-10 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
+                            if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
                                 tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
                                 bool_sublim = .true.
@@ -817 +816 @@
                     enddo
                 endif
+                is_co2ice_ini(ig,islope) = .false.
                 if ((co2_ice(ig,islope) < 1.e-10) .and. (h2o_ice(ig,islope) > frost_albedo_threshold)) then
                     albedo(ig,1,islope) = albedo_h2o_frost
@@ -1159 +1159 @@
 deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
 deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
-deallocate(co2_ice_ini,ini_co2ice_sublim,ini_h2oice_sublim,stratif)
+deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif)
 !----------------------------- END OUTPUT ------------------------------