Changeset 3537 for trunk/LMDZ.COMMON/libf

Timestamp:

Dec 6, 2024, 5:47:19 PM (2 weeks ago)

Author:

jbclement

Message:

PEM:
Small fixes to initialize and output ice table-related variables and in the launching script.
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• deftank/launchPEM.sh (modified) (1 diff)
• deftank/lib_launchPEM.sh (modified) (5 diffs)
• pem.F90 (modified) (4 diffs)
• pemetat0.F90 (modified) (8 diffs)
• pemredem.F90 (modified) (3 diffs)

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

@@ -507 +507 @@
 == 04/12/2024 == JBC
 Fixing the writing format for "info_PEM.txt".
+
+== 06/12/2024 == JBC
+Small fixes to initialize and output ice table-related variables and in the launching script.

trunk/LMDZ.COMMON/libf/evolution/deftank/launchPEM.sh

@@ -114 +114 @@
             echo "The number of PCM runs between each PEM run has been modified from $nPCM_old to $nPCM."
         fi
-        if [ $n_myear -ne $n_myear_old ]; then
+        if [ "$(echo "$n_myear != $n_myear_old" | bc)" -eq 1 ]; then
             echo "The number of initial PCM runs has been modified from $n_myear_old to $n_myear."
         fi

trunk/LMDZ.COMMON/libf/evolution/deftank/lib_launchPEM.sh

@@ -295 +295 @@
     iPCM=$(($irelaunch + 1))
     cleanfiles diags/diagfi .nc $irelaunch
+    cleanfiles diags/diagsoil .nc $irelaunch
     cleanfiles diags/data2reshape .nc $irelaunch
     cleanfiles "out_PCM/run" "" $irelaunch
@@ -310 +311 @@
         iPEM=1
         cleanfiles diags/diagpem .nc $iPEM
+        cleanfiles diags/diagsoilpem .nc $iPEM
         cleanfiles "out_PEM/run" "" $iPEM
         cleanfiles starts/restart1D_postPEM .txt $iPEM
@@ -333 +335 @@
         il=$(echo "($irelaunch - $nPCM_ini + 1)%$nPCM + 1" | bc)
         cleanfiles diags/diagpem .nc $iPEM
+        cleanfiles diags/diagsoilpem .nc $iPEM
         cleanfiles "out_PEM/run" "" $iPEM
         cleanfiles starts/restart1D_postPEM .txt $iPEM
@@ -366 +369 @@
     sed -i "1s/.*/$i_myear $n_myear $convert_years $iPCM $iPEM $nPCM $nPCM_ini/" info_PEM.txt
     cleanfiles diags/diagfi .nc $(($iPCM - 1))
+    cleanfiles diags/diagsoil .nc $(($iPCM - 1))
     cleanfiles "out_PCM/run" "" $(($iPCM - 1))
     cleanfiles starts/restart1D .txt $(($iPCM - 1))
@@ -372 +376 @@
     cleanfiles diags/data2reshape .nc $(($iPCM - 1))
     cleanfiles diags/diagpem .nc $irelaunch
+    cleanfiles diags/diagsoilpem .nc $irelaunch
     cleanfiles "out_PEM/run" "" $irelaunch
     cleanfiles starts/restart1D_postPEM .txt $irelaunch

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

@@ -995 +995 @@
         if (icetable_equilibrium) then
             call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
-            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope))
+            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table thickness','m',2,icetable_thickness(:,islope))
         else if (icetable_dynamic) then
             call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
-            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope))
         endif
 
@@ -1004 +1003 @@
             call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil_PEM','K',3,tsoil_PEM(:,:,islope))
             call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI_PEM','K',3,TI_PEM(:,:,islope))
+            if (icetable_dynamic) call writediagpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope))
             if (adsorption_pem) then
                 call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbded_phys(:,:,islope))
@@ -1070 +1070 @@
 !------------------------
 ! III_a.1 Ice update (for startfi)
-
 watercap = 0.
 perennial_co2ice = co2_ice
@@ -1210 +1209 @@
 if (layering_algo) nb_str_max = get_nb_str_max(stratif,ngrid,nslope) ! Get the maximum number of "stratum" in the stratification (layerings)
 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,   &
+call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
              co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif)
 write(*,*) "restartpem.nc has been written"

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

@@ -7 +7 @@
 !=======================================================================
 
-SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table_depth,ice_table_thickness, &
+SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness, &
                     ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,d_h2oice,d_co2ice,co2_ice,h2o_ice,    &
                     global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,                    &
@@ -16 +16 @@
 use comsoil_h,                  only: volcapa, inertiedat
 use adsorption_mod,             only: regolith_adsorption, adsorption_pem
-use ice_table_mod,              only: computeice_table_equilibrium, icetable_depth, icetable_thickness, icetable_equilibrium, icetable_dynamic
+use ice_table_mod,              only: computeice_table_equilibrium, icetable_equilibrium, icetable_dynamic
 use constants_marspem_mod,      only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
 use soil_thermalproperties_mod, only: update_soil_thermalproperties
@@ -60 +60 @@
 real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
 real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
-real, dimension(ngrid,nslope),                   intent(inout) :: ice_table_depth     ! Ice table depth [m]
-real, dimension(ngrid,nslope),                   intent(inout) :: ice_table_thickness ! Ice table thickness [m]
+real, dimension(ngrid,nslope),                   intent(inout) :: icetable_depth      ! Ice table depth [m]
+real, dimension(ngrid,nslope),                   intent(inout) :: icetable_thickness  ! Ice table thickness [m]
 real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: ice_porefilling     ! Subsurface ice pore filling [m3/m3]
 real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst          ! instantaneous soil (mid-layer) temperature [K]
@@ -108 +108 @@
 
 !0.2 Set to default values
-ice_table_depth = -1. ! by default, no ice table
-ice_table_thickness = -1.
+icetable_depth = -1. ! by default, no ice table
+icetable_thickness = -1.
+ice_porefilling = 0.
 !1. Run
 if (startpem_file) then
@@ -305 +306 @@
 !3. Ice Table
         if (icetable_equilibrium) then
-            call get_field("ice_table_depth",ice_table_depth,found)
+            call get_field("icetable_depth",icetable_depth,found)
             if (.not. found) then
-                write(*,*)'PEM settings: failed loading <ice_table_depth>'
+                write(*,*)'PEM settings: failed loading <icetable_depth>'
                 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_depth,ice_table_thickness)
-                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
+                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
+                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
                 do islope = 1,nslope
                     call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
@@ -318 +319 @@
         else if (icetable_dynamic) then
             call get_field("ice_porefilling",ice_porefilling,found)
+            if (.not. found) write(*,*)'PEM settings: failed loading <ice_porefilling>'
+            call get_field("icetable_depth",icetable_depth,found)
             if (.not. found) then
-                write(*,*)'PEM settings: failed loading <ice_porefilling>'
-                ice_porefilling = 0.
-            endif
-            call get_field("ice_table_depth",ice_table_depth,found)
-            if (.not. found) then
-                write(*,*)'PEM settings: failed loading <ice_table_depth>'
+                write(*,*)'PEM settings: failed loading <icetable_depth>'
                 write(*,*)'will reconstruct the values of the ice table given the current state'
-                ice_table_depth = -9999.
-                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
+                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
                 do islope = 1,nslope
                     call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
@@ -498 +495 @@
 !c) Ice table
         if (icetable_equilibrium) then
-            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
-            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
+            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
+            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
             do islope = 1,nslope
                 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
@@ -505 +502 @@
             write(*,*) 'PEMETAT0: Ice table done'
         else if (icetable_dynamic) then
-            ice_porefilling = 0.
-            ice_table_depth = -9999.
-            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
+            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
             do islope = 1,nslope
                 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))

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

@@ -58 +58 @@
 
 SUBROUTINE pemdem1(filename,i_myear,nsoil_PEM,ngrid,nslope,tsoil_slope_PEM,inertiesoil_slope_PEM, &
-                   ice_table_depth,ice_table_thickness,ice_porefilling,m_co2_regolith,m_h2o_regolith,h2o_ice,stratif)
+                   icetable_depth,icetable_thickness,ice_porefilling,m_co2_regolith,m_h2o_regolith,h2o_ice,stratif)
 
 ! write time-dependent variable to restart file
@@ -77 +77 @@
 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_slope_PEM       ! under mesh bining according to slope
 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: inertiesoil_slope_PEM ! under mesh bining according to slope
-real, dimension(ngrid,nslope),           intent(in) :: ice_table_depth       ! under mesh bining according to slope
-real, dimension(ngrid,nslope),           intent(in) :: ice_table_thickness   ! under mesh bining according to slope
+real, dimension(ngrid,nslope),           intent(in) :: icetable_depth       ! under mesh bining according to slope
+real, dimension(ngrid,nslope),           intent(in) :: icetable_thickness   ! under mesh bining according to slope
 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: ice_porefilling       ! under mesh bining according to slope
 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: m_co2_regolith, m_h2o_regolith
@@ -122 +122 @@
         call put_field("mh2o_reg_ads_slope"//num, "Mass of h2o adsorbded in the regolith",m_h2o_regolith(:,:,islope),Year)
     enddo
-    call put_field("ice_table_depth","Depth of ice table",ice_table_depth,Year)
-    call put_field("ice_table_thickness","Depth of ice table",ice_table_thickness,Year)
+    call put_field("icetable_depth","Depth of ice table",icetable_depth,Year)
+    call put_field("icetable_thickness","Depth of ice table",icetable_thickness,Year)
     call put_field("ice_porefilling","Subsurface ice pore filling",ice_porefilling,Year)
     call put_field("inertiedat_PEM","Thermal inertie of PEM ",inertiedat_PEM,Year)