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

Timestamp:

Sep 18, 2024, 4:09:02 PM (2 months ago)

Author:

jbclement

Message:

PEM:

• Correction of the way stopping criteria interact with each other because of a situation where the algorithm does not stop even though a stopping criterion was met.
• Update and corrections for the layering algorithm: there is now only one subsurface stratum.
• Improvement of the launching script in the case of the 1D PEM: it ends automatically after PCM/PEM run crash.

JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• deftank/lib_launchPEM.sh (modified) (5 diffs)
• layering_mod.F90 (modified) (10 diffs)
• pem.F90 (modified) (2 diffs)
• pemetat0.F90 (modified) (7 diffs)
• stopping_crit_mod.F90 (modified) (2 diffs)

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

@@ -422 +422 @@
 == 16/09/2024 == JBC
 Corrections and improvements for the launching script.
+
+== 18/09/2024 == JBC
+- Correction of the way stopping criteria interact with each other because of a situation where the algorithm does not stop even though a stopping criterion was met.
+- Update and corrections for the layering algorithm: there is now only one subsurface stratum.
+- Improvement of the launching script in the case of the 1D PEM: it ends automatically after PCM/PEM run crash.

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

@@ -175 +175 @@
             sed -i "s/^k=[0-9]\+$/k=$(echo "3 - $nPCM_ini" | bc -l)/" PCMrun.job
             ./PCMrun.job
+            if [ $? -ne 0 ]; then
+                errlaunch
+            fi
         else # 3D model
             cp PCMrun.job PCMrun${iPCM}.job
@@ -197 +200 @@
                 sed -i "s/^k=[0-9]\+$/k=$(echo "$i + 2 - $nPCM_ini" | bc -l)/" PCMrun.job
                 ./PCMrun.job
+                if [ $? -ne 0 ]; then
+                    errlaunch
+                fi
             else # 3D model
                 cp PCMrun.job PCMrun${iPCM}.job
@@ -219 +225 @@
         if [ $1 -eq 1 ]; then # 1D model
             ./PEMrun.job
+            if [ $? -ne 0 ]; then
+                errlaunch
+            fi
         else # 3D model
             sed -i -E "s/($name_job[^0-9]*[0-9]*[^0-9]*)[0-9]+$/\1${iPEM}/" PEMrun.job
@@ -245 +254 @@
         if [ $1 -eq 1 ]; then # 1D model
             ./PEMrun.job
+            if [ $? -ne 0 ]; then
+                errlaunch
+            fi
         else # 3D model
             sed -i -E "s/($name_job[^0-9]*[0-9]*[^0-9]*)[0-9]+$/\1${iPEM}/" PEMrun.job
@@ -351 +363 @@
 # arg1: model dimension
 relaunchPEM() {
-    iPEM=$irelaunch
+    iPEM=$(($irelaunch + 1))
     iPCM=$(($nPCM_ini + $nPCM*($irelaunch - 1) + 1))
     i_myear=$(awk "NR==$(($iPEM + 1)) {print \$1}" "info_PEM.txt")

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

@@ -273 +273 @@
 ! Creation of three strata at the bottom of the layering
 ! 1) Dry porous deep regolith
-call add_stratum(this,18.,-12.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity)
+!call add_stratum(this,18.,-12.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity)
 ! 2) Ice-cemented regolith (ice table)
-call add_stratum(this,10.,-2.,0.,1. - dry_regolith_porosity,dry_regolith_porosity,0.)
+!call add_stratum(this,10.,-2.,0.,1. - dry_regolith_porosity,dry_regolith_porosity,0.)
 ! 3) Dry porous regolith
 call add_stratum(this,2.,0.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity)
@@ -404 +404 @@
 
 !---- Local variables
-type(stratum), pointer :: current
-integer                :: ig, islope, k
+integer :: ig, islope, k
 
 !---- Code
 do islope = 1,nslope
     do ig = 1,ngrid
-        current => stratif(ig,islope)%bottom
-        do k = 1,3
-            current%thickness      = stratif_array(ig,islope,k,1)
-            current%top_elevation  = stratif_array(ig,islope,k,2)
-            current%co2ice_volfrac = stratif_array(ig,islope,k,3)
-            current%h2oice_volfrac = stratif_array(ig,islope,k,4)
-            current%dust_volfrac   = stratif_array(ig,islope,k,5)
-            current%air_volfrac    = stratif_array(ig,islope,k,6)
-            current => current%up
-        enddo
-        do k = 4,size(stratif_array,3)
+        stratif(ig,islope)%bottom%thickness      = stratif_array(ig,islope,1,1)
+        stratif(ig,islope)%bottom%top_elevation  = stratif_array(ig,islope,1,2)
+        stratif(ig,islope)%bottom%co2ice_volfrac = stratif_array(ig,islope,1,3)
+        stratif(ig,islope)%bottom%h2oice_volfrac = stratif_array(ig,islope,1,4)
+        stratif(ig,islope)%bottom%dust_volfrac   = stratif_array(ig,islope,1,5)
+        stratif(ig,islope)%bottom%air_volfrac    = stratif_array(ig,islope,1,6)
+        do k = 2,size(stratif_array,3)
             if (all(abs(stratif_array(ig,islope,k,:)) < eps)) exit
             call add_stratum(stratif(ig,islope),stratif_array(ig,islope,k,1),stratif_array(ig,islope,k,2),stratif_array(ig,islope,k,3),stratif_array(ig,islope,k,4),stratif_array(ig,islope,k,5),stratif_array(ig,islope,k,6))
@@ -501 +496 @@
         else ! No, we stop
             write(*,'(a)') ' Stratification -> Dust lifting: no dust lag!'
-            stopPEM = 7
+            stopPEM = 8
             exit
         endif
@@ -507 +502 @@
     if (h2lift_tot > 0.) then
         write(*,'(a)') ' Stratification -> Not enough dust in the lag layers to complete the lifting!'
-        stopPEM = 7
+        stopPEM = 8
     endif
 
@@ -553 +548 @@
 
             ! How much is CO2 ice sublimation inhibited by the top dust lag?
-            currentloc => current1%up
             h_toplag = 0.
-            do while (associated(currentloc) .and. abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
-                h_toplag = h_toplag + currentloc%thickness
-                currentloc => currentloc%up
-            enddo
+            if (associated(current1%up)) then ! If there is an upper stratum
+                currentloc => current1%up
+                do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
+                    h_toplag = h_toplag + currentloc%thickness
+                    if (.not. associated(currentloc%up)) exit
+                    currentloc => currentloc%up
+                enddo
+            endif
             h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag)
 
@@ -584 +582 @@
         if (h2subl_tot > 0.) then
             write(*,'(a)') ' Stratification -> Not enough CO2 ice in the layering to complete the sublimation!'
-            stopPEM = 7
+            stopPEM = 8
         endif
         ! New stratum at the layering top by condensation of H2O ice
@@ -609 +607 @@
 
             ! How much is CO2 ice sublimation inhibited by the top dust lag?
-            currentloc => current1%up
             h_toplag = 0.
-            do while (associated(currentloc) .and. abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
-                h_toplag = h_toplag + currentloc%thickness
-                currentloc => currentloc%up
-            enddo
+            if (associated(current1%up)) then ! If there is an upper stratum
+                currentloc => current1%up
+                do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
+                    h_toplag = h_toplag + currentloc%thickness
+                    if (.not. associated(currentloc%up)) exit
+                    currentloc => currentloc%up
+                enddo
+            endif
             h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag)
 
@@ -640 +641 @@
         if (h2subl_tot > 0.) then
             write(*,'(a)') ' Stratification -> Not enough CO2 ice in the layering to complete the sublimation!'
-            stopPEM = 7
+            stopPEM = 8
         endif
         ! H2O ice sublimation in the considered stratum + New stratum for dust lag
@@ -649 +650 @@
             h_dust_old = current2%dust_volfrac*current2%thickness
 
-            ! How much is CO2 ice sublimation inhibited by the top dust lag?
-            currentloc => current2%up
+            ! How much is H2O ice sublimation inhibited by the top dust lag?
             h_toplag = 0.
-            do while (associated(currentloc) .and. abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
-                h_toplag = h_toplag + currentloc%thickness
-                currentloc => currentloc%up
-            enddo
+            if (associated(current2%up)) then ! If there is an upper stratum
+                currentloc => current2%up
+                do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
+                    h_toplag = h_toplag + currentloc%thickness
+                    if (.not. associated(currentloc%up)) exit
+                    currentloc => currentloc%up
+                enddo
+            endif
             h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag)
 
@@ -681 +685 @@
         if (h2subl_tot > 0.) then
             write(*,'(a)') ' Stratification -> Not enough H2O ice in the layering to complete the sublimation!'
-            stopPEM = 7
+            stopPEM = 8
         endif
 

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

@@ -1001 +1001 @@
     year_iter = year_iter + dt_pem
     i_myear = i_myear + dt_pem
-    if (year_iter >= year_iter_max) stopPEM = 5
-    if (i_myear >= n_myear) stopPEM = 6
+    if (stopPEM <= 0 .and. year_iter >= year_iter_max) stopPEM = 5
+    if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6
     call system_clock(c2)
-    if (timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7
-    if (stopPEM > 0) then
+    if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7
+    if (stopPEM <= 0) then
         select case (stopPEM)
             case(1)
@@ -1021 +1021 @@
             case(7)
                 write(*,*) "STOPPING because the time limit for the PEM job will be reached soon:", stopPEM
+            case(8)
+                write(*,*) "STOPPING because the layering algorithm met an hasty end:", stopPEM
             case default
                 write(*,*) "STOPPING with unknown stopping criterion code:", stopPEM

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

@@ -116 +116 @@
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+#ifndef CPP_STD
+    ! h2o ice
+    h2o_ice = 0.
+    call get_field("h2o_ice",h2o_ice,found)
+    if (.not. found) then
+        write(*,*)'Pemetat0: failed loading <h2o_ice>'
+        write(*,*)'will reconstruct the values from watercaptag'
+        write(*,*)'with default value ''ini_huge_h2oice'''
+        do ig = 1,ngrid
+            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
+        enddo
+    else
+        ! The variations of infinite reservoirs during the PCM years are taken into account
+        h2o_ice = h2o_ice + watercap
+    endif
+
+    ! co2 ice
+    co2_ice = perennial_co2ice
+#endif
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     ! Stratification (layerings)
-    nb_str_max = 3
+    nb_str_max = 1
     if (layering_algo) then
         found = inquire_dimension("nb_str_max")
@@ -149 +170 @@
     endif
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-#ifndef CPP_STD
-    ! h2o ice
-    h2o_ice = 0.
-    call get_field("h2o_ice",h2o_ice,found)
-    if (.not. found) then
-        write(*,*)'Pemetat0: failed loading <h2o_ice>'
-        write(*,*)'will reconstruct the values from watercaptag'
-        write(*,*)'with default value ''ini_huge_h2oice'''
-        do ig = 1,ngrid
-            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
-        enddo
-    else
-        ! The variations of infinite reservoirs during the PCM years are taken into account
-        h2o_ice = h2o_ice + watercap
-    endif
-
-    ! co2 ice
-    co2_ice = perennial_co2ice
-#endif
-
     if (soil_pem) then
 
@@ -186 +186 @@
                         !!! transition
                         delta = depth_breccia
-                        TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
-                                                                 (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
-                                                                 ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
-                        do iloop=index_breccia+2,index_bedrock
+                        TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
+                                                                  (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
+                                                                  ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
+                        do iloop=index_breccia + 2,index_bedrock
                             TI_PEM(ig,iloop,islope) = TI_breccia
                         enddo
@@ -224 +224 @@
             delta = depth_breccia
             do ig = 1,ngrid
-                if (inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
+                if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
                     inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
                                                               (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
                                                               ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
 
-                    do iloop = index_breccia+2,index_bedrock
+                    do iloop = index_breccia + 2,index_bedrock
                         inertiedat_PEM(ig,iloop) = TI_breccia
                     enddo
                 else
-                    do iloop=index_breccia+1,index_bedrock
+                    do iloop=index_breccia + 1,index_bedrock
                         inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
                     enddo
@@ -242 +242 @@
             delta = depth_bedrock
             do ig = 1,ngrid
-                inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
-                                                          (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
-                                                          ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
+                inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
+                                                           (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
+                                                           ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
             enddo
 
@@ -358 +358 @@
     write(*,*)'There is no "startpem.nc"!'
 
-    ! Stratification (layerings)
-    nb_str_max = 3
-    if (layering_algo) write(*,*)'So the stratification (layerings) is initialized with only the 3 sub-surface strata.'
-
     ! h2o ice
     h2o_ice = 0.
@@ -372 +368 @@
     write(*,*)'So ''co2_ice'' is initialized with ''perennial_co2ice''.'
     co2_ice = perennial_co2ice
+
+    ! Stratification (layerings)
+    nb_str_max = 1
+    if (layering_algo) write(*,*)'So the stratification (layerings) is initialized with 1 sub-surface stratum.'
 
     if (soil_pem) then

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

@@ -42 +42 @@
 
 !=======================================================================
+if (stopPEM < 0) return
+
 ! Computation of the present surface of h2o ice still sublimating
 h2oice_now_surf = 0.
@@ -104 +106 @@
 
 !=======================================================================
+if (stopPEM < 0) return
+
 ! Computation of the present surface of co2 ice still sublimating
 co2ice_now_surf = 0.