Changeset 3170 for trunk/LMDZ.COMMON/libf/evolution

Timestamp:

Jan 4, 2024, 11:37:33 AM (11 months ago)

Author:

llange

Message:

Mars PEM
Fixing a small bug: the subroutine compute_icetable was always called, even if tthe option 'icetable_equilibrium' was set to false in the run_PEM.def. It is now fixed by adding a flag before the call.
LL

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• pem.F90 (modified) (2 diffs)
• pemetat0.F90 (modified) (5 diffs)

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

@@ -186 +186 @@
 - 'ini_h2o_bigreservoir' is renamed into 'ini_huge_h2oice'.
 - Some cleanings, in particular for the main program "pem.F90".
+
+== 04/01/2024 == LL
+Fixing a small bug: the subroutine compute_icetable was always called, even if tthe option 'icetable_equilibrium' was set to false in the run_PEM.def. It is now fixed by adding a flag before the call.

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

@@ -55 +55 @@
 use recomp_orb_param_mod,       only: recomp_orb_param
 use ice_table_mod,              only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
-                                      ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
+                                      ini_ice_table_porefilling, icetable_equilibrium, computeice_table_equilibrium,compute_massh2o_exchange_ssi
 use soil_thermalproperties_mod, only: update_soil_thermalproperties
 use time_phylmdz_mod,           only: daysec, dtphys, day_end
@@ -838 +838 @@
 
 ! II_d.3 Update the ice table
-        write(*,*) "Compute ice table"
-        porefillingice_thickness_prev_iter = porefillingice_thickness
-        call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
-        call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
-
+        if (icetable_equilibrium) then
+            write(*,*) "Compute ice table"
+            porefillingice_thickness_prev_iter = porefillingice_thickness
+            call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
+            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
+        endif
 ! II_d.4 Update the soil thermal properties
         call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)

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

@@ -16 +16 @@
 use comsoil_h,                  only: volcapa, inertiedat
 use adsorption_mod,             only: regolith_adsorption, adsorption_pem
-use ice_table_mod,              only: computeice_table_equilibrium
+use ice_table_mod,              only: computeice_table_equilibrium, icetable_equilibrium
 use constants_marspem_mod,      only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
 use soil_thermalproperties_mod, only: update_soil_thermalproperties
@@ -102 +102 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-!0. Check if the start_PEM exist.
+!0.1 Check if the start_PEM exist.
 
 inquire(file = filename,exist = startpem_file)
@@ -108 +108 @@
 write(*,*)'Is start PEM?',startpem_file
 
+!0.2 Set to default values
+ice_table = -1.  ! by default, no ice table
+ice_table_thickness = -1.
 !1. Run
 if (startpem_file) then
@@ -272 +275 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !3. Ice Table
-        call get_field("ice_table",ice_table,found)
-        if (.not. found) then
-            write(*,*)'PEM settings: failed loading <ice_table>'
-            write(*,*)'will reconstruct the values of the ice table given the current state'
-            ice_table = -1  ! by default, no ice table
-            ice_table_thickness = -1  ! by default, no ice table
-            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table,ice_table_thickness)
-            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
-            do islope = 1,nslope
-                call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
-            enddo
+        if(icetable_equilibrium) then
+            call get_field("ice_table",ice_table,found)
+            if (.not. found) then
+                write(*,*)'PEM settings: failed loading <ice_table>'
+                write(*,*)'will reconstruct the values of the ice table given the current state'
+                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table,ice_table_thickness)
+                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
+                do islope = 1,nslope
+                    call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
+                enddo
+            endif
+            write(*,*) 'PEMETAT0: ICE TABLE done'
         endif
-
-        write(*,*) 'PEMETAT0: ICE TABLE done'
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -442 +444 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !c) Ice table
-        ice_table = -1.  ! by default, no ice table
-        ice_table_thickness = -1.
-        call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table,ice_table_thickness)
-        call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
-        do islope = 1,nslope
-            call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
-        enddo
-        write(*,*) 'PEMETAT0: Ice table done'
-
+        if(icetable_equilibrium) then
+            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table,ice_table_thickness)
+            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
+            do islope = 1,nslope
+                call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
+            enddo
+            write(*,*) 'PEMETAT0: Ice table done'
+        endif
+        
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !d) Regolith adsorbed