Changeset 3331 for trunk/LMDZ.COMMON

Timestamp:

May 16, 2024, 6:40:20 PM (6 months ago)

Author:

jbclement

Message:

PEM:
Update of the layering algorithm + corrections of wrong lines commited in r3330.
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• layering_mod.F90 (modified) (8 diffs)
• pem.F90 (modified) (2 diffs)

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

@@ -315 +315 @@
 == 16/05/2024 == JBC
 Improvement of a flag-like variable (more robust as a logical) to know where co2 ice was at the beginning.
+
+== 16/05/2024 == JBC
+Update of the layering algorithm + corrections of wrong lines commited in r3330.

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

@@ -474 +474 @@
 htend_dust = tend_dust/rho_dust
 
-if (htend_dust < 0.) then ! Dust lifting
+if (htend_dust < 0.) then ! Dust lifting only
     if (abs(htend_co2ice) > eps .or. abs(htend_h2oice) > eps) error stop 'Situation not managed: dust lifting + CO2/H2O ice condensation/sublimation!'
     write(*,'(a)') ' Stratification -> Dust lifting'
@@ -506 +506 @@
 else
 
-!------ Dust sedimentation
+!------ Dust sedimentation only
     if (abs(htend_co2ice) < eps .and. abs(htend_h2oice) < eps) then
         write(*,'(a)') ' Stratification -> Dust sedimentation'
         ! New stratum at the layering top by sedimentation of dust
         thickness = htend_dust/(1. - dry_regolith_porosity)
-        if (thickness > 0.) then ! Only if the stratum is builiding up
+        if (thickness > 0.) then ! Only if the stratum is building up
              if (new_str) then
                  call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,0.,htend_dust/thickness,dry_regolith_porosity)
@@ -526 +526 @@
         ! New stratum at the layering top by condensation of CO2 and H2O ice
         thickness = htend_co2ice/(1. - co2ice_porosity) + htend_h2oice/(1. - h2oice_porosity) + htend_dust
-        if (thickness > 0.) then ! Only if the stratum is builiding up
+        if (thickness > 0.) then ! Only if the stratum is building up
              if (new_str) then
                  call add_stratum(this,thickness,this%top%top_elevation + thickness,htend_co2ice/thickness,htend_h2oice/thickness,htend_dust/thickness,1. - (htend_co2ice/thickness + htend_h2oice/thickness + htend_dust/thickness))
@@ -550 +550 @@
                 h2subl_tot = 0.
                 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
-            else if (h_co2ice_old < eps) then ! There is nothing to sublimate so we move to the underlying stratum
-                current1 => current1%down
-                new_lag1 = .true.
+            else if (h_co2ice_old < eps) then ! There is nothing to sublimate so we stop here
+                !current1 => current1%down ! Move to the underlying stratum
+                !new_lag1 = .true.
+                exit
             else ! Only a fraction can sublimate and so we move to the underlying stratum
                 h2subl = h_co2ice_old
@@ -567 +568 @@
         ! New stratum at the layering top by condensation of H2O ice
         thickness = htend_h2oice/(1. - h2oice_porosity) + htend_dust
-        if (thickness > 0.) then ! Only if the stratum is builiding up
+        if (thickness > 0.) then ! Only if the stratum is building up
              if (new_str) then
                  call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,htend_h2oice/thickness,htend_dust/thickness,1. - (htend_h2oice/thickness + htend_dust/thickness))
@@ -591 +592 @@
                 h2subl_tot = 0.
                 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
-            else if (h_co2ice_old < eps) then ! There is nothing to sublimate so we move to the underlying stratum
-                current1 => current1%down
-                new_lag1 = .true.
+            else if (h_co2ice_old < eps) then ! There is nothing to sublimate so we stop here
+                !current1 => current1%down ! Move to the underlying stratum
+                !new_lag1 = .true.
+                exit
             else ! Only a fraction can sublimate and so we move to the underlying stratum
                 h2subl = h_co2ice_old
@@ -617 +619 @@
                 h2subl_tot = 0.
                 call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag2)
-            else if (h_h2oice_old < eps) then ! There is nothing to sublimate so we move to the underlying stratum
-                current2 => current2%down
-                new_lag2 = .true.
+            else if (h_h2oice_old < eps) then ! There is nothing to sublimate so we stop here
+                !current2 => current1%down ! Move to the underlying stratum
+                !new_lag2 = .true.
+                exit
             else ! Only a fraction can sublimate and so we move to the underlying stratum
                 h2subl = h_h2oice_old
@@ -682 +685 @@
 endif
 
-
 ! New stratum = dust lag
 h_lag = hsubl_dust/(1. - dry_lag_porosity)

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

@@ -792 +792 @@
         do islope = 1,nslope
             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.
+                co2ice_disappeared(ig,islope) = .true.
                 if (latitude_deg(ig) > 0) then
                     do ig_loop = ig,ngrid
@@ -816 +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