Changeset 2885

Timestamp:

Jan 31, 2023, 9:52:30 PM (22 months ago)

Author:

romain.vande

Message:

Mars PCM:
Move a endif misplaced
RV

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• comsoil_h_PEM.F90 (modified) (3 diffs)
• criterion_ice_stop_mod_slope.F90 (modified) (4 diffs)
• pem.F90 (modified) (11 diffs)
• pemetat0.F90 (modified) (1 diff)
• pemredem.F90 (modified) (3 diffs)
• read_data_GCM.F90 (modified) (5 diffs)

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

@@ -22 +22 @@
   real, save, allocatable :: h2o_adsorbded_phys(:,:,:)  ! h2o that is in the regolith [kg/m^2]
   LOGICAL soil_pem  ! True by default, to run with the subsurface physic. Read in pem.def
+  real,save,allocatable,dimension(:) :: water_reservoir      ! Large reserve of water   [kg/m^2]
 
 contains
@@ -44 +45 @@
     allocate(co2_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
     allocate(h2o_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
+    allocate(water_reservoir(ngrid))
   end subroutine ini_comsoil_h_PEM
 
@@ -64 +66 @@
     if (allocated(co2_adsorbded_phys)) deallocate(co2_adsorbded_phys)
     if (allocated(h2o_adsorbded_phys)) deallocate(h2o_adsorbded_phys)
+    if (allocated(water_reservoir)) deallocate(water_reservoir)
   end subroutine end_comsoil_h_PEM
 

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

@@ -2 +2 @@
 ! $Id $
 !
-SUBROUTINE criterion_ice_stop_slope(cell_area,ini_surf,qsurf,STOPPING,ngrid,initial_h2o_ice,global_ave_press_GCM,global_ave_press_new,nslope)
+SUBROUTINE criterion_ice_stop_slope(cell_area,ini_surf,qsurf,STOPPING,STOPPING_ps,ngrid,initial_h2o_ice,global_ave_press_GCM,global_ave_press_new,nslope)
 
   USE temps_mod_evol, ONLY: alpha_criterion
@@ -28 +28 @@
 
 !   OUTPUT
-  LOGICAL, intent(out) :: STOPPING              ! Logical : is the criterion reached?
+  LOGICAL, intent(out) :: STOPPING,STOPPING_ps              ! Logical : is the criterion reached?
 
 !   local:
@@ -39 +39 @@
 !   initialisation to false
     STOPPING=.FALSE.
+    STOPPING_ps=.FALSE.
 
 !   computation of the actual surface
@@ -67 +68 @@
   if(global_ave_press_new.LT.global_ave_press_GCM*(0.9) .OR. &
      global_ave_press_new.GT.global_ave_press_GCM*(1.1)) then
-    STOPPING=.TRUE. 
+    STOPPING_ps=.TRUE. 
     print *, "Reason of stopping : The global pressure reach the threshold:"
     print *, "Current global pressure=", global_ave_press_new

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

@@ -37 +37 @@
                            albedodat, zmea, zstd, zsig, zgam, zthe, &
                            hmons, summit, base,albedo_h2o_frost, &
-                           frost_albedo_threshold,emissiv
+                           frost_albedo_threshold,emissiv,watercaptag
       use dimradmars_mod, only: totcloudfrac, albedo
       use dust_param_mod, only: tauscaling
@@ -90 +90 @@
                               TI_PEM,inertiedat_PEM,alph_PEM, beta_PEM,         & ! soil thermal inertia          
                               tsoil_PEM, mlayer_PEM,layer_PEM,                  & !number of subsurface layers, soil mid layer depths
-                              fluxgeo, co2_adsorbded_phys, h2o_adsorbded_phys     ! geothermal flux, mass of co2 & h2o in the regolith
+                              fluxgeo, co2_adsorbded_phys, h2o_adsorbded_phys,  & ! geothermal flux, mass of co2 & h2o in the regolith
+                              water_reservoir
       use adsorption_mod, only : regolith_adsorption
 
@@ -187 +188 @@
       LOGICAL :: STOPPING_co2     ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
       LOGICAL :: STOPPING_1_co2   ! Logical : is there still water ice to sublimate?
+      LOGICAL :: STOPPING_pressure
+      INTEGER :: criterion_stop            !Which criterion is reached : 1=surf h20, 2=surf_co2, 3=ps, 4=orb.param
+
 
       REAL, dimension(:,:,:),allocatable  :: q_co2_GCM ! Lon x Lat x Time : mass mixing ratio of co2 in the first layer [kg/kg]
@@ -275 +279 @@
     INTEGER :: year_iter_max             ! maximum number of iterations before stopping
     REAL :: timestep                     ! timestep [s] 
+    REAL WC_sum
+
 #ifdef CPP_STD
 !     INTEGER :: nsplit_phys=1
@@ -418 +424 @@
        enddo
      enddo
+
+     call surfini(ngrid,qsurf)
 
 #else
@@ -793 +801 @@
       watersurf_density_phys_ave,watersoil_density_phys_PEM_ave, &
       co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded)
+
+      do ig=1,ngrid
+         do islope=1,nslope
+           qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+watercap_slope(ig,islope)+water_reservoir(ig)
+         enddo
+      enddo
 
     if(soil_pem) then
@@ -1139 +1153 @@
       call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
 
-      call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,ngrid,initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
+      call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,STOPPING_pressure,ngrid,initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
 
       year_iter=year_iter+dt_pem
@@ -1146 +1160 @@
       if (STOPPING_water) then
         print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water 
+        criterion_stop=1
+      endif
+      if (STOPPING_1_water) then
+        print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
+        criterion_stop=1
+      endif
+      if (STOPPING_co2) then
+        print *, "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2
+        criterion_stop=2
+      endif
+      if (STOPPING_1_co2) then
+        print *, "STOPPING because tendencies on co2 ice=0, see message above", STOPPING_1_co2
+        criterion_stop=2
+      endif
+      if (STOPPING_pressure) then
+        print *, "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure
+        criterion_stop=3
       endif
       if (year_iter.ge.year_iter_max) then
         print *, "STOPPING because maximum number of iterations reached"
+        criterion_stop=4
       endif
-      if (STOPPING_1_water) then
-        print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
-      endif
-      if (STOPPING_co2) then
-        print *, "STOPPING because surface of co2 ice sublimating is too low or global pressure changed too much, see message above", STOPPING_co2
-      endif
-      if (STOPPING_1_co2) then
-        print *, "STOPPING because tendencies on co2 ice=0, see message above", STOPPING_1_co2
-      endif
-
-
-
-      if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2)  then
+
+
+
+      if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2 .or. STOPPING_pressure)  then
         exit
       else
@@ -1197 +1220 @@
         qsurf(ig,igcm_co2_ice)=co2ice(ig)
 #endif
-      ENDDO ! of DO ig=1,ngrid
-! H2o ice
-      DO ig = 1,ngrid
-        qsurf(ig,igcm_h2o_ice) = 0.
-        DO islope = 1,nslope
-          qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
-                                   * subslope_dist(ig,islope) /          &
-                                  cos(pi*def_slope_mean(islope)/180.)
-        ENDDO
       ENDDO ! of DO ig=1,ngrid
 
@@ -1319 +1333 @@
        enddo
      enddo
+
+     DO ig=1,ngrid
+       if(watercaptag(ig)) then
+         print *, "INNN"
+         WC_sum=0.
+         DO islope=1,nslope
+           if(qsurf_slope(ig,igcm_h2o_ice,islope).GT. (watercap_slope(ig,islope)+water_reservoir(ig))) then
+             qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-(watercap_slope(ig,islope)+water_reservoir(ig))
+           else
+             watercap_slope(ig,islope)=watercap_slope(ig,islope)+qsurf_slope(ig,igcm_h2o_ice,islope)-(watercap_slope(ig,islope)+water_reservoir(ig)/nslope)
+             qsurf_slope(ig,igcm_h2o_ice,islope)=0.
+           endif
+           WC_sum=WC_sum+watercap_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
+           watercap_slope(ig,islope)=0.
+         enddo  
+           water_reservoir(ig)=water_reservoir(ig)+WC_sum
+       endif
+     enddo
+
+! H2o ice
+      DO ig = 1,ngrid
+        qsurf(ig,igcm_h2o_ice) = 0.
+        DO islope = 1,nslope
+          qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
+                                   * subslope_dist(ig,islope) /          &
+                                  cos(pi*def_slope_mean(islope)/180.)
+        ENDDO
+      ENDDO ! of DO ig=1,ngrid
+
+     DO ig=1,ngrid
+!       DO islope=1,nslope
+         if(qsurf(ig,igcm_h2o_ice).GT.500) then
+            watercaptag(ig)=.true.
+       DO islope=1,nslope
+            qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-250
+            water_reservoir(ig)=water_reservoir(ig)+250
+       ENDDO            
+         endif
+!       enddo
+     enddo
+
+     DO ig=1,ngrid
+!       DO islope=1,nslope
+         if(water_reservoir(ig).LT. 10) then
+            watercaptag(ig)=.false.
+            qsurf(ig,igcm_h2o_ice)=qsurf(ig,igcm_h2o_ice)+water_reservoir(ig)
+       DO islope=1,nslope
+            qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+water_reservoir(ig)
+       ENDDO
+         endif
+!       enddo
+     enddo
+
+      DO ig = 1,ngrid
+        watercap(ig) = 0.
+        DO islope = 1,nslope
+          watercap(ig) = watercap(ig) + watercap_slope(ig,islope) &
+                                   * subslope_dist(ig,islope) /          &
+                                  cos(pi*def_slope_mean(islope)/180.)
+        ENDDO
+      ENDDO ! of DO ig=1,ngrid
+
      
 !------------------------
@@ -1386 +1462 @@
 
          call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
-                      tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys,h2o_adsorbded_phys)
+                      tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
+
+     call info_run_PEM(year_iter, criterion_stop)
 
       print *, "restartfi_PEM.nc has been written"

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

@@ -4 +4 @@
    
    use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
-   use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM
+   use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM, water_reservoir
    use comsoil_h, only:  volcapa,inertiedat 
    use adsorption_mod, only : regolith_adsorption
    USE temps_mod_evol, ONLY: year_PEM
+#ifndef CPP_STD   
+   use surfdat_h, only: watercaptag
+#endif
 
 

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

@@ -73 +73 @@
 subroutine pemdem1(filename,year_iter,nsoil_PEM,ngrid,nslope, &
                     tsoil_slope_PEM,inertiesoil_slope_PEM,ice_table, &
-                    m_co2_regolith,m_h2o_regolith)
+                    m_co2_regolith,m_h2o_regolith,water_reservoir)
   ! write time-dependent variable to restart file
   use iostart_PEM, only : open_restartphy, close_restartphy, & 
@@ -96 +96 @@
   real,intent(in) :: m_co2_regolith(ngrid,nsoil_PEM,nslope)
   real,intent(in) :: m_h2o_regolith(ngrid,nsoil_PEM,nslope)
+  real,intent(IN) :: water_reservoir(ngrid)
   integer :: islope
   CHARACTER*2 :: num  
@@ -115 +116 @@
   ! set time counter in file
   call put_var("Time","Year of simulation",year_tot)
+
+  call put_field("water_reservoir","water_reservoir", &
+                 water_reservoir,year_tot)
 
     if(soil_pem) then

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

@@ -69 +69 @@
   REAL, ALLOCATABLE ::  co2_ice_s(:,:,:)                               ! co2 ice, mesh averaged, of the concatenated file [kg/m^2]
   REAL, ALLOCATABLE ::  h2o_ice_s_slope(:,:,:,:)                       ! h2o ice per slope of the concatenated file [kg/m^2]
+  REAL, ALLOCATABLE :: watercap_slope(:,:,:,:)
 
 !-----------------------------------------------------------------------
@@ -80 +81 @@
       allocate(co2_ice_s(iim+1,jjm+1,timelen))
       allocate(h2o_ice_s_slope(iim+1,jjm+1,nslope,timelen))
+      allocate(watercap_slope(iim_input+1,jjm_input+1,nslope,timelen))
       allocate(h2o_ice_s(iim+1,jjm+1,timelen))
 
@@ -133 +135 @@
 
      print *, "Downloading data for h2o_ice_s_slope done"
+     print *, "Downloading data for watercap_slope ..."
+
+DO islope=1,nslope
+  write(num,fmt='(i2.2)') islope
+  call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
+ENDDO
+
+     print *, "Downloading data for watercap_slope done"
      print *, "Downloading data for tsurf_slope ..."
 
@@ -198 +208 @@
 
   print *, "Computing the min of h2o_ice_slope"
-  min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope,4)
+  min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope+watercap_slope,4)
   print *, "Computing the min of co2_ice_slope"
   min_co2_ice_slope(:,:,:)=minval(co2_ice_slope,4)
@@ -225 +235 @@
             min_co2_ice_slope(i,j,islope)  = 0.
           endif
-          if (min_h2o_ice_slope(i,j,islope).LT.0) then
-            min_h2o_ice_slope(i,j,islope)  = 0.
-          endif
+!          if (min_h2o_ice_slope(i,j,islope).LT.0) then
+!            min_h2o_ice_slope(i,j,islope)  = 0.
+!          endif
        ENDDO
     ENDDO