Changeset 3527 for trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90

Timestamp:

Nov 20, 2024, 3:53:19 PM (26 hours ago)

Author:

jbclement

Message:

PEM:

• Removing redundant Norbert Schorghofer's subroutines/parameters (follow-up of r3526);
• Making all Norbert Schorghofer's subroutines with modern explicit interface via modules;
• Cleaning of "glaciers_mod.F90";
• Optimization for the computation of ice density according to temperature by using a function.

JBC

File:

• trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90 (modified) (9 diffs)

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

@@ -184 +184 @@
         do islope = 1,nslope
             call compute_Tice_pem(nsoil,tsoil(ig,:,islope),tsurf(ig,islope),icetable_depth(ig,islope),Tice)
-            rho = -3.5353e-4*Tice**2 + 0.0351*Tice + 933.5030 ! Rottgers, 2012
-            delta_m_h2o(ig) = delta_m_h2o(ig) + porosity*rho*(icetable_thickness(ig,islope) - icetable_thickness_old(ig,islope)) & ! convention > 0. <=> it condenses
+            delta_m_h2o(ig) = delta_m_h2o(ig) + porosity*rho_ice(Tice,'h2o')*(icetable_thickness(ig,islope) - icetable_thickness_old(ig,islope)) & ! convention > 0. <=> it condenses
                               *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
         enddo
@@ -195 +194 @@
             do isoil = 1,nsoil
                 call compute_Tice_pem(nsoil,tsoil(ig,:,islope),tsurf(ig,islope),mlayer_PEM(isoil - 1),Tice)
-                rho = -3.5353e-4*Tice**2 + 0.0351*Tice + 933.5030 ! Rottgers, 2012
-                delta_m_h2o(ig) = delta_m_h2o(ig) + rho*(ice_porefilling(ig,isoil,islope) - ice_porefilling_old(ig,isoil,islope)) & ! convention > 0. <=> it condenses
+                delta_m_h2o(ig) = delta_m_h2o(ig) + rho_ice(Tice,'h2o')*(ice_porefilling(ig,isoil,islope) - ice_porefilling_old(ig,isoil,islope)) & ! convention > 0. <=> it condenses
                                   *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
             enddo
@@ -258 +256 @@
     index_tmp = nsoilmx_PEM
     do ilay = 1,nsoilmx_PEM
-        call constriction(porefill(ilay),eta(ilay))
+        eta(ilay) = constriction(porefill(ilay))
     enddo
 else
     index_tmp = index_IS
     do ilay = 1,index_IS - 1
-        call constriction(porefill(ilay),eta(ilay))
+        eta(ilay) = constriction(porefill(ilay))
     enddo
     do ilay = index_IS,nsoilmx_PEM
@@ -350 +348 @@
 
 !-----------------------------------------------------------------------
-SUBROUTINE constriction(porefill,eta)
+FUNCTION constriction(porefill) RESULT(eta)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -363 +361 @@
 
 real, intent(in) :: porefill ! pore filling fraction
-real, intent(out) :: eta ! constriction
+real :: eta ! constriction
 
 if (porefill <= 0.) then
@@ -373 +371 @@
 endif
 
-END SUBROUTINE constriction
+END FUNCTION constriction
 
 !-----------------------------------------------------------------------
@@ -386 +384 @@
 
 use comsoil_h_PEM, only: layer_PEM, mlayer_PEM
+use abort_pem_mod, only: abort_pem
 
 implicit none
@@ -422 +421 @@
     endif
     if (indexice < 0) then
-        call abort_physic("compute_Tice_pem","subsurface ice is below the last soil layer",1)
+        call abort_pem("compute_Tice_pem","Subsurface ice is below the last soil layer!",1)
     else
         if(indexice >= 1) then ! Linear inteprolation between soil temperature
@@ -434 +433 @@
 END SUBROUTINE compute_Tice_pem
 
+!-----------------------------------------------------------------------
+FUNCTION rho_ice(T,ice) RESULT(rho)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Compute subsurface ice density
+!!!
+!!! Author: JBC
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+use abort_pem_mod, only: abort_pem
+
+implicit none
+
+real,         intent(in) :: T
+character(3), intent(in) :: ice
+real :: rho
+
+select case (trim(adjustl(ice)))
+    case('h2o')
+        rho = -3.5353e-4*T**2 + 0.0351*T + 933.5030 ! Rottgers 2012
+    case('co2')
+        rho = (1.72391 - 2.53e-4*T-2.87*1e-7*T**2)*1.e3 ! Mangan et al. 2017
+    case default
+        call abort_pem("rho_ice","Type of ice not known!",1)
+end select
+
+END FUNCTION rho_ice
+
 END MODULE ice_table_mod