Changeset 3512

Timestamp:

Nov 12, 2024, 6:35:33 PM (6 weeks ago)

Author:

jbclement

Message:

PEM:
Few corrections related to r3498 (time step from integer to real) and r3493 (Norbert Schorghofer's subroutines for dynamic ice table) in order to make the code work properly.
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• NS_dyn_ss_ice_m.F90 (modified) (2 diffs)
• NS_fast_subs_mars.F90 (modified) (1 diff)
• changelog.txt (modified) (1 diff)
• ice_table_mod.F90 (modified) (1 diff)
• info_PEM_mod.F90 (modified) (2 diffs)
• pem.F90 (modified) (2 diffs)
• pemetat0.F90 (modified) (3 diffs)
• writediagpem.F90 (modified) (3 diffs)

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

@@ -9 +9 @@
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
 
 SUBROUTINE dyn_ss_ice_m(ssi_depth_in,T1,Tb,nz,thIn,p0,pfrost,porefill_in,porefill,ssi_depth)
@@ -69 +67 @@
   !call setgrid(nz,z,zmax,zfac)
   l1=2.e-4
-  do iloop=0,nz
-    z(iloop) = l1*(1+iloop**2.9*(1-exp(-real(iloop)/20.)))
+  do iloop=0,nz - 1
+    z(iloop + 1) = l1*(1+iloop**2.9*(1-exp(-real(iloop)/20.)))
   enddo
   

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

@@ -116 +116 @@
      call icechanges(nz,z(:),typeF,avdrho(k),avdrhoP(k),ypp(:), &
           & Diff,porosity,icefrac,bigstep,zdepthT(k),porefill(:,k),typeG)
-
+     if (typeP>2) then
      if (mode2 .and. porefill(typeP,k)>=1. .and. porefill(typeP-1,k)==0.) then  ! nothing changed
         porefill(typeP-1,k)=1.  ! paste a layer
    !     write(34,*) '# mode 2 growth occurred',typeP,typeF,typeT
+     endif
      endif
 

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

@@ -477 +477 @@
 - Renaming of the tendencies in the PEM with the prefix 'd_' instead of 'tend_';
 - Modification of the PEM time step type from integer to real. As a consequence, all time variables are now of real type. This change adds the possibility to consider fractions of year as time step.
+
+== 12/11/2024 == JBC
+Few corrections related to r3498 (time step from integer to real) and r3493 (Norbert Schorghofer's subroutines for dynamic ice table) in order to make the code work properly.

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

@@ -28 +28 @@
 
 allocate(icetable_depth(ngrid,nslope))
-if (icetable_equilibrium) then
-    allocate(icetable_thickness(ngrid,nslope))
-else if (icetable_dynamic) then
-    allocate(ice_porefilling(ngrid,nsoil,nslope))
-endif
+allocate(icetable_thickness(ngrid,nslope))
+allocate(ice_porefilling(ngrid,nsoil,nslope))
 
 END SUBROUTINE ini_ice_table

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

@@ -31 +31 @@
 logical       :: ok
 integer       :: cstat
-character(10) :: ich1, ich2, ich3, ich4, ich5, ich6, fch
+character(10) :: frmt
+character(20) :: ich1, ich2, ich3, ich4, fch1, fch2, fch3
 
 !----- Code
 inquire(file = 'info_PEM.txt',exist = ok)
 if (ok) then
-    write(ich1,'(f0.4)') i_myear
-    write(ich2,'(f0.4)') n_myear
-    write(fch,'(f0.4)') convert_years ! 4 digits to the right of the decimal point to respect the precision of Martian year in "launch_pem.sh"
-    write(ich3,'(i0)') iPCM
-    write(ich4,'(i0)') iPEM + 1
-    write(ich5,'(i0)') nPCM
-    write(ich6,'(i0)') nPCM_ini
-    call execute_command_line('sed -i "1s/.*/'//trim(ich1)//' '//trim(ich2)//' '//trim(fch)//' '//trim(ich3)//' '//trim(ich4)//' '//trim(ich5)//' '//trim(ich6)//'/" info_PEM.txt',cmdstat = cstat)
+    write(fch1,'(f'//int2str(nb_digits(i_myear) + 5)//'.4)') i_myear
+    write(fch2,'(f'//int2str(nb_digits(n_myear) + 5)//'.4)') n_myear
+    write(fch3,'(f6.4)') convert_years ! 4 digits to the right of the decimal point to respect the precision of Martian year in "launch_pem.sh"
+    write(ich1,'(i0)') iPCM
+    write(ich2,'(i0)') iPEM + 1
+    write(ich3,'(i0)') nPCM
+    write(ich4,'(i0)') nPCM_ini
+    call execute_command_line('sed -i "1s/.*/'//trim(fch1)//' '//trim(fch2)//' '//trim(fch3)//' '//trim(ich1)//' '//trim(ich2)//' '//trim(ich3)//' '//trim(ich4)//'/" info_PEM.txt',cmdstat = cstat)
     if (cstat > 0) then
         error stop 'info_PEM: command execution failed!'
@@ -60 +61 @@
 END SUBROUTINE info_PEM
 
+!=======================================================================
+
+FUNCTION int2str(i) RESULT(str)
+! Function to convert an integer into a string
+
+integer, intent(in) :: i
+character(20)       :: str
+
+if (nb_digits(real(i)) > len(str)) error stop 'int2str [info_PEM_mod]: invalid integer for conversion!'
+write(str,'(i0)') i
+str = trim(adjustl(str))
+
+END FUNCTION int2str
+
+!=======================================================================
+
+FUNCTION nb_digits(x) RESULT(idigits)
+! Function to give the number of digits for the integer part of a real number
+
+real, intent(in) :: x
+integer          :: idigits
+
+idigits = 1
+! If x /= 0 then:
+if (abs(x) > 1.e-10) idigits = int(log10(abs(x))) + 1
+
+END FUNCTION nb_digits
+
 END MODULE info_PEM_mod

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

@@ -941 +941 @@
             enddo
             deallocate(porefill)
-            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
+            !call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
         endif
 
@@ -949 +949 @@
 ! II_d.5 Update the mass of the regolith adsorbed
         if (adsorption_pem) then
-            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,                   &
+            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,                           &
                                      h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
                                      h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)

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

@@ -326 +326 @@
                 write(*,*)'PEM settings: failed loading <ice_table_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)
+                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,TI_PEM)
                 do islope = 1,nslope
@@ -335 +335 @@
             if (.not. found) then
                 write(*,*)'PEM settings: failed loading <ice_porefilling>'
-                ice_porefilling = 0.
+                ice_porefilling = 1.
             endif
             write(*,*) 'PEMETAT0: ICE TABLE done'
@@ -510 +510 @@
             write(*,*) 'PEMETAT0: Ice table done'
         else if (icetable_dynamic) then
-            ice_porefilling = 0.
-            ice_table_depth = 0.
+            ice_porefilling = 1.
+            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,TI_PEM)
             do islope = 1,nslope

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

@@ -51 +51 @@
 use mod_phys_lmdz_para, only: is_parallel, is_mpi_root, is_master, gather
 use mod_grid_phy_lmdz,  only: klon_glo, Grid1Dto2D_glo, nbp_lon, nbp_lat, nbp_lev, grid_type, unstructured
-use time_evol_mod,      only: ecritpem, dt
+use time_evol_mod,      only: ecritpem
 
 implicit none
@@ -600 +600 @@
 use mod_grid_phy_lmdz,  only: klon_glo, Grid1Dto2D_glo, nbp_lon, nbp_lat
 use mod_grid_phy_lmdz,  only: grid_type, unstructured
-use time_evol_mod,      only: ecritpem, dt
+use time_evol_mod,      only: ecritpem
 use iniwritesoil_mod,   only: iniwritesoil
 
@@ -747 +747 @@
 if (name.eq.firstname) then
   ! if we run across 'firstname', then it is a new time step
-  zitau = zitau + dt
+  zitau = zitau + 1
 endif