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

Timestamp:

Apr 26, 2024, 5:03:35 PM (7 months ago)

Author:

jbclement

Message:

PEM:
Small correction about the dimension of an array.
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• conf_pem.F90 (modified) (3 diffs)
• deftank/launch_pem.sh (modified) (1 diff)
• deftank/run_PEM.def (modified) (2 diffs)
• evol_ice_mod.F90 (modified) (2 diffs)
• layering_mod.F90 (modified) (1 diff)
• pem.F90 (modified) (7 diffs)
• pemetat0.F90 (modified) (4 diffs)
• pemredem.F90 (modified) (2 diffs)

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

@@ -29 +29 @@
 use glaciers_mod,          only: h2oice_flow, co2ice_flow, inf_h2oice_threshold, metam_co2ice_threshold, metam_h2oice_threshold, metam_h2oice, metam_co2ice
 use ice_table_mod,         only: icetable_equilibrium, icetable_dynamic
+use layering_mod,          only: layering_algo
 
 implicit none
@@ -103 +104 @@
 write(*,*)  'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp
 
-!#---------- Layering parameters ----------#
 fluxgeo = 0.
 call getin('fluxgeo',fluxgeo)
@@ -174 +174 @@
 call getin('co2ice_flow',co2ice_flow)
 
+!#---------- Layering parameters ----------#
+layering_algo = .false.
+call getin('layering',layering_algo)
+
 END SUBROUTINE conf_pem
 

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

@@ -18 +18 @@
 fi
 
+# Save the current value of LC_NUMERIC and set it to a locale that uses a dot as the decimal separator
+OLD_LC_NUMERIC=$LC_NUMERIC
+LC_NUMERIC=en_US.UTF-8
+
+
+#####################################################################
 # A few parameters that might need be changed depending on your setup:
 # Path to the arch.env to source:
 source ../trunk/LMDZ.COMMON/arch.env
 
-# Save the current value of LC_NUMERIC and set it to a locale that uses a dot as the decimal separator
-OLD_LC_NUMERIC=$LC_NUMERIC
-LC_NUMERIC=en_US.UTF-8
-
-
-#####################################################################
 #---------- Modify here the number of years to be simulated ------------
 ## set the number of years, either Martian or Earth years:

trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def

@@ -52 +52 @@
 # reg_thprop_dependp=.false.
 
-#---------- Layering parameters ----------#
 # Value of the geothermal flux? Default = 0.
 # fluxgeo=0.
@@ -93 +92 @@
 # co2ice_flow=.true.
 
+!#---------- Layering parameters ----------#
+# Do you want to run with the layering algorithm? Default = .false.
+# layering=.true.
+
 # Some definitions for the physics, in file 'callphys.def'
 INCLUDEDEF=callphys.def

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

@@ -106 +106 @@
     enddo
 
+    new_tendencies = 0.
     if (abs(pos_tend) < 1.e-10 .or. abs(neg_tend) < 1.e-10) then
         write(*,*) "Reason of stopping: there is no sublimating or condensing h2o ice!"
@@ -112 +113 @@
         write(*,*) "This can be due to the absence of h2o ice in the PCM run!"
         stopPEM = 2
-        new_tendencies = 0.
     else
         ! We adapt the tendencies to conserve h2o and do only exchange between grid points

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

@@ -15 +15 @@
 real, parameter :: eps = epsilon(1.), tol = 1.e2*eps
 integer         :: nb_str_max
+logical         :: layering_algo
 
 ! Physical parameters

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

@@ -70 +70 @@
 use writediagpem_mod,           only: writediagpem, writediagsoilpem
 use co2condens_mod,             only: CO2cond_ps
-use layering_mod,               only: tend_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max
+use layering_mod,               only: tend_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo
 
 #ifndef CPP_STD
@@ -351 +351 @@
     status = nf90_close(ncid)
 
-    call iniphysiq('startfi.nc',iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys,rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
+    call iniphysiq(iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys,rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
 #else
     ps_start_PCM(1) = ps(1)
@@ -589 +589 @@
 
 allocate(stratif(ngrid,nslope))
-do islope = 1,nslope
-    do i = 1,ngrid
-        call ini_layering(stratif(i,islope))
-    enddo
-enddo
+if (layering_algo) then
+    do islope = 1,nslope
+        do i = 1,ngrid
+            call ini_layering(stratif(i,islope))
+        enddo
+    enddo
+endif
 
 call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
@@ -655 +657 @@
 year_iter = 0
 stopPEM = 0
-allocate(new_str(ngrid,nslope),new_lag1(ngrid,nslope),new_lag2(ngrid,nslope),current1(ngrid,nslope),current2(ngrid,nslope))
-new_str = .true.
-new_lag1 = .true.
-new_lag2 = .true.
-do islope = 1,nslope
-    do ig = 1,ngrid
-        current1(ig,islope)%p => stratif(ig,islope)%top
-        current2(ig,islope)%p => stratif(ig,islope)%top
-    enddo
-enddo
+if (layering_algo) then
+    allocate(new_str(ngrid,nslope),new_lag1(ngrid,nslope),new_lag2(ngrid,nslope),current1(ngrid,nslope),current2(ngrid,nslope))
+    new_str = .true.
+    new_lag1 = .true.
+    new_lag2 = .true.
+    do islope = 1,nslope
+        do ig = 1,ngrid
+            current1(ig,islope)%p => stratif(ig,islope)%top
+            current2(ig,islope)%p => stratif(ig,islope)%top
+        enddo
+    enddo
+endif
 
 do while (year_iter < 10 .and. i_myear < n_myear)
@@ -748 +752 @@
     call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,tend_h2o_ice,stopPEM)
     call evol_co2_ice(ngrid,nslope,co2_ice,tend_co2_ice)
-    do islope = 1,nslope
-        do ig = 1,ngrid
-            call make_layering(stratif(ig,islope),tend_co2_ice(ig,islope),tend_h2o_ice(ig,islope),tend_dust,new_str(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p)
-        enddo
-    enddo
+    if (layering_algo) then
+        do islope = 1,nslope
+            do ig = 1,ngrid
+                call make_layering(stratif(ig,islope),tend_co2_ice(ig,islope),tend_h2o_ice(ig,islope),tend_dust,new_str(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p)
+            enddo
+        enddo
+    endif
 
 !------------------------
@@ -1122 +1128 @@
 !    III_c Write the "restartpem.nc"
 !------------------------
-nb_str_max = get_nb_str_max(stratif,ngrid,nslope) ! Get the maximum number of "stratum" in the stratification (layerings)
+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, &
@@ -1136 +1142 @@
 write(*,*) "LL & RV & JBC: so far, so good!"
 
-do islope = 1,nslope
-    do i = 1,ngrid
-        call del_layering(stratif(i,islope))
-    enddo
-enddo
-deallocate(stratif,new_str,new_lag1,new_lag2,current1,current2)
+if (layering_algo) then
+    do islope = 1,nslope
+        do i = 1,ngrid
+            call del_layering(stratif(i,islope))
+        enddo
+    enddo
+    deallocate(new_str,new_lag1,new_lag2,current1,current2)
+endif
 deallocate(vmr_co2_PCM,ps_timeseries,tsurf_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
 deallocate(co2_ice_PCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfavg_before_saved)
 deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
 deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
-deallocate(co2_ice_ini)
+deallocate(co2_ice_ini,stratif)
 !----------------------------- END OUTPUT ------------------------------
 

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

@@ -23 +23 @@
 use compute_soiltemp_mod,       only: ini_tsoil_pem, compute_tsoil_pem
 use comslope_mod,               only: def_slope_mean, subslope_dist
-use layering_mod,               only: layering, array2stratif
+use layering_mod,               only: layering, array2stratif, nb_str_max, layering_algo
 
 #ifndef CPP_STD
@@ -78 +78 @@
 logical                                 :: startpem_file                   ! boolean to check if we read the startfile or not
 real, dimension(:,:,:,:), allocatable   :: stratif_array                   ! Array for stratification (layerings)
-real                                    :: nb_str_max
 
 #ifdef CPP_STD
@@ -118 +117 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     ! Stratification (layerings)
-    found = inquire_dimension("nb_str_max")
-    if (.not. found) then
-        write(*,*) 'Pemetat0: failed loading <nb_str_max>!'
-        write(*,*) '''nb_str_max'' is set to 10!'
-        nb_str_max = 10.
-    else
-        nb_str_max = inquire_dimension_length('nb_str_max')
+    nb_str_max = 3
+    if (layering_algo) then
+        found = inquire_dimension("nb_str_max")
+        if (.not. found) then
+            write(*,*) 'Pemetat0: failed loading <nb_str_max>!'
+            write(*,*) '''nb_str_max'' is set to 3!'
+        else
+            nb_str_max = int(inquire_dimension_length('nb_str_max'))
+        endif
+        allocate(stratif_array(ngrid,nslope,nb_str_max,6))
+        stratif_array = 0.
+        do islope = 1,nslope
+            write(num,'(i2.2)') islope
+            call get_field('stratif_slope'//num//'_thickness',stratif_array(:,islope,:,1),found)
+            call get_field('stratif_slope'//num//'_top_elevation',stratif_array(:,islope,:,2),found2)
+            found = found .or. found2
+            call get_field('stratif_slope'//num//'_co2ice_volfrac',stratif_array(:,islope,:,3),found2)
+            found = found .or. found2
+            call get_field('stratif_slope'//num//'_h2oice_volfrac',stratif_array(:,islope,:,4),found2)
+            found = found .or. found2
+            call get_field('stratif_slope'//num//'_dust_volfrac',stratif_array(:,islope,:,5),found2)
+            found = found .or. found2
+            call get_field('stratif_slope'//num//'_air_volfrac',stratif_array(:,islope,:,6),found2)
+            found = found .or. found2
+            if (.not. found) then
+                write(*,*) 'Pemetat0: failed loading at least one of the properties of <stratif_slope'//num//'>'
+            endif ! not found
+        enddo ! islope
+        call array2stratif(stratif_array,ngrid,nslope,stratif)
+        deallocate(stratif_array)
     endif
-    allocate(stratif_array(ngrid,nslope,int(nb_str_max),6))
-    stratif_array = 0.
-    do islope = 1,nslope
-        write(num,'(i2.2)') islope
-        call get_field('stratif_slope'//num//'_thickness',stratif_array(:,islope,:,1),found)
-        call get_field('stratif_slope'//num//'_top_elevation',stratif_array(:,islope,:,2),found2)
-        found = found .or. found2
-        call get_field('stratif_slope'//num//'_co2ice_volfrac',stratif_array(:,islope,:,3),found2)
-        found = found .or. found2
-        call get_field('stratif_slope'//num//'_h2oice_volfrac',stratif_array(:,islope,:,4),found2)
-        found = found .or. found2
-        call get_field('stratif_slope'//num//'_dust_volfrac',stratif_array(:,islope,:,5),found2)
-        found = found .or. found2
-        call get_field('stratif_slope'//num//'_air_volfrac',stratif_array(:,islope,:,6),found2)
-        found = found .or. found2
-        if (.not. found) then
-            write(*,*) 'Pemetat0: failed loading at least one of the properties of <stratif_slope'//num//'>'
-        endif ! not found
-    enddo ! islope
-    call array2stratif(stratif_array,ngrid,nslope,stratif)
-    deallocate(stratif_array)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -358 +359 @@
 
     ! Stratification (layerings)
-    write(*,*)'So the stratification (layerings) is initialized with only the 3 sub-surface strata.'
     nb_str_max = 3
+    if (layering_algo) write(*,*)'So the stratification (layerings) is initialized with only the 3 sub-surface strata.'
 
     ! h2o ice

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

@@ -64 +64 @@
 use comsoil_h_PEM, only: inertiedat_PEM, soil_pem
 use time_evol_mod, only: year_bp_ini, convert_years
-use layering_mod,  only: layering, nb_str_max, stratif2array, print_layering
+use layering_mod,  only: layering, nb_str_max, stratif2array, print_layering, layering_algo
 
 implicit none
@@ -96 +96 @@
 call put_field('h2o_ice','h2o_ice',h2o_ice,Year)
 
-allocate(stratif_array(ngrid,nslope,nb_str_max,6))
-call stratif2array(stratif,ngrid,nslope,stratif_array)
-do islope = 1,nslope
-    write(num,fmt='(i2.2)') islope
-    call put_field('stratif_slope'//num//'_thickness','Layering thickness',stratif_array(:,islope,:,1),Year)
-    call put_field('stratif_slope'//num//'_top_elevation','Layering top elevation',stratif_array(:,islope,:,2),Year)
-    call put_field('stratif_slope'//num//'_co2ice_volfrac','Layering CO2 ice volume fraction',stratif_array(:,islope,:,3),Year)
-    call put_field('stratif_slope'//num//'_h2oice_volfrac','Layering H2O ice volume fraction',stratif_array(:,islope,:,4),Year)
-    call put_field('stratif_slope'//num//'_dust_volfrac','Layering dust volume fraction',stratif_array(:,islope,:,5),Year)
-    call put_field('stratif_slope'//num//'_air_volfrac','Layering air volume fraction',stratif_array(:,islope,:,6),Year)
-enddo
-deallocate(stratif_array)
+if (layering_algo) then
+    allocate(stratif_array(ngrid,nslope,nb_str_max,6))
+    call stratif2array(stratif,ngrid,nslope,stratif_array)
+    do islope = 1,nslope
+        write(num,fmt='(i2.2)') islope
+        call put_field('stratif_slope'//num//'_thickness','Layering thickness',stratif_array(:,islope,:,1),Year)
+        call put_field('stratif_slope'//num//'_top_elevation','Layering top elevation',stratif_array(:,islope,:,2),Year)
+        call put_field('stratif_slope'//num//'_co2ice_volfrac','Layering CO2 ice volume fraction',stratif_array(:,islope,:,3),Year)
+        call put_field('stratif_slope'//num//'_h2oice_volfrac','Layering H2O ice volume fraction',stratif_array(:,islope,:,4),Year)
+        call put_field('stratif_slope'//num//'_dust_volfrac','Layering dust volume fraction',stratif_array(:,islope,:,5),Year)
+        call put_field('stratif_slope'//num//'_air_volfrac','Layering air volume fraction',stratif_array(:,islope,:,6),Year)
+    enddo
+    deallocate(stratif_array)
+endif
 
 if (soil_pem) then