Changeset 3297

Timestamp:

Apr 8, 2024, 3:49:41 PM (8 months ago)

Author:

jbclement

Message:

PEM:
Integration of the module "layering_mod.F90" with the rest of the PEM:

• The linked list data structure representative of layered deposits is converted into an array which can be outputed in the "restratpem.nc" files. This array has dimensions (ngrid,nslope,nb_str_max,6) where 'nb_str_max' is the maximum number of 'stratum' through the layerings and '6' is the number of properties of 'stratum';
• this structure can also be read from "startpem.nc" files to initialize PEM runs;
• The layering algorithm is now used in the main PEM loop to make the layerings evolve.

JBC

Location:

trunk/LMDZ.COMMON

Files:

• arch/arch-gfortran.fcm (modified) (2 diffs)
• arch/arch-gfortran.path (modified) (1 diff)
• libf/evolution/changelog.txt (modified) (1 diff)
• libf/evolution/info_PEM_mod.F90 (modified) (1 diff)
• libf/evolution/iostart_PEM.F90 (modified) (56 diffs)
• libf/evolution/layering_mod.F90 (modified) (27 diffs)
• libf/evolution/pem.F90 (modified) (9 diffs)
• libf/evolution/pemetat0.F90 (modified) (7 diffs)
• libf/evolution/pemredem.F90 (modified) (4 diffs)

trunk/LMDZ.COMMON/arch/arch-gfortran.fcm

@@ -1 @@
-%COMPILER            gfortran
-%LINK                gfortran
+%COMPILER            mpif90
+%LINK                mpif90
 %AR                  ar
 %MAKE                make
@@ -7 +7 @@
 %CPP_DEF             LAPACK
 %FPP_DEF             NC_DOUBLE
-%BASE_FFLAGS         -c -fdefault-real-8 -fdefault-double-8 -ffree-line-length-none -fno-align-commons -fallow-argument-mismatch
-%PROD_FFLAGS         -O3
-%DEV_FFLAGS          -O
-%DEBUG_FFLAGS        -ffpe-trap=invalid,zero,overflow -fbounds-check -g3 -O0 -fstack-protector-all -finit-real=snan -fbacktrace
+%BASE_FFLAGS         -pg -c -fdefault-real-8 -fdefault-double-8 -ffree-line-length-none -fno-align-commons -fallow-argument-mismatch
+%PROD_FFLAGS         -O2
+%DEV_FFLAGS          -O -Wall -fbounds-check
+%DEBUG_FFLAGS        -fcheck=all -ffree-line-length-0 -Wall -ffpe-trap=invalid,zero,overflow -fbounds-check -g3 -O0 -fstack-protector-all -finit-real=nan -fbacktrace
 %C_COMPILER          gcc
-%C_OPTIM             -O3
-%MPI_FFLAGS
+%C_OPTIM             -O2
+%MPI_FFLAGS          -I/usr/lib/x86_64-linux-gnu/openmpi/include
 %OMP_FFLAGS         
-%BASE_LD     
-%MPI_LD
+%BASE_LD             -pg
+%MPI_LD              -lmpi
 %OMP_LD              

trunk/LMDZ.COMMON/arch/arch-gfortran.path

@@ -1 +1 @@
 ROOT=$PWD
+NETCDF="$HOME/netcdf4hdf5"
 
 NETCDF_LIBDIR="-L$NETCDF/lib"

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

@@ -265 +265 @@
 == 08/04/2024 == JBC
 Correction of "launch_pem.sh" in the deftank: the number of years to be simulated was not respected because of extra PCM runs at the end of the simulation due to the order of PCM/PEM runs inside the loop.
+
+== 08/04/2024 == JBC
+Integration of the module "layering_mod.F90" with the rest of the PEM:
+    - The linked list data structure representative of layered deposits is converted into an array which can be outputed in the "restratpem.nc" files. This array has dimensions (ngrid,nslope,nb_str_max,6) where 'nb_str_max' is the maximum number of 'stratum' through the layerings  and '6' is the number of properties of 'stratum';
+    - this structure can also be read from "startpem.nc" files to initialize PEM runs;
+    - The layering algorithm is now used in the main PEM loop to make the layerings evolve.

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

@@ -38 +38 @@
     call execute_command_line('sed -i "1s/.*/'//trim(ich1)//' '//trim(ich2)//' '//trim(fch)//'/" info_PEM.txt',cmdstat = cstat)
     if (cstat > 0) then
-        error stop 'info_PEM: command exection failed!'
+        error stop 'info_PEM: command execution failed!'
     else if (cstat < 0) then
         error stop 'info_PEM: command execution not supported!'

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

@@ -2 +2 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
-!!! Purpose: Read and write specific netcdf for the PEM 
-!!! 
-!!! 
+!!! Purpose: Read and write specific netcdf for the PEM
+!!!
+!!!
 !!! Author: LL, inspired by iostart from the PCM
 !!!
@@ -13 +13 @@
     INTEGER,SAVE :: nid_start ! NetCDF file identifier for startfi.nc file
     INTEGER,SAVE :: nid_restart ! NetCDF file identifier for restartfi.nc file
-    
+
     ! restartfi.nc file dimension identifiers: (see open_restartphy())
-    INTEGER,SAVE :: idim1 ! "index" dimension
-    INTEGER,SAVE :: idim2 ! "physical_points" dimension
-    INTEGER,SAVE :: idim3 ! "subsurface_layers" dimension
-    INTEGER,SAVE :: idim4 ! "nlayer_plus_1" dimension
-    INTEGER,SAVE :: idim5 ! "number_of_advected_fields" dimension
-    INTEGER,SAVE :: idim6 ! "nlayer" dimension
-    INTEGER,SAVE :: idim7 ! "Time" dimension
-    INTEGER,SAVE :: idim8 ! "nslope" dimension
-    INTEGER,SAVE :: idim9 ! "inter slope" dimension
+    INTEGER,SAVE :: idim1  ! "index" dimension
+    INTEGER,SAVE :: idim2  ! "physical_points" dimension
+    INTEGER,SAVE :: idim3  ! "subsurface_layers" dimension
+    INTEGER,SAVE :: idim4  ! "nlayer_plus_1" dimension
+    INTEGER,SAVE :: idim5  ! "number_of_advected_fields" dimension
+    INTEGER,SAVE :: idim6  ! "nlayer" dimension
+    INTEGER,SAVE :: idim7  ! "Time" dimension
+    INTEGER,SAVE :: idim8  ! "nslope" dimension
+    INTEGER,SAVE :: idim9  ! "inter slope" dimension
+    INTEGER,SAVE :: idim10 ! "nb_str_max" dimension
     INTEGER,SAVE :: timeindex ! current time index (for time-dependent fields)
     INTEGER,PARAMETER :: length=100 ! size of tab_cntrl array
-    
+
     INTERFACE get_field
       MODULE PROCEDURE Get_field_r1,Get_field_r2,Get_field_r3
     END INTERFACE get_field
-    
+
     INTERFACE get_var
       MODULE PROCEDURE get_var_r0,Get_var_r1,Get_var_r2,Get_var_r3
@@ -48 +49 @@
     PUBLIC inquire_field, inquire_field_ndims
     PUBLIC open_startphy,close_startphy,open_restartphy,close_restartphy
-    
+
 CONTAINS
 
@@ -66 +67 @@
       ENDIF
     ENDIF
-    
+
     CALL bcast(nid_start) ! tell all procs about nid_start
-  
+
   END SUBROUTINE open_startphy
 
@@ -93 +94 @@
     INTEGER :: varid
     INTEGER :: ierr
-    
+
     IF (is_master) THEN
       ierr=NF90_INQ_VARID(nid_start,Field_name,varid)
@@ -109 +110 @@
 
   FUNCTION inquire_field_ndims(Field_name)
-  ! give the number of dimensions of "Field_name" stored in the input file 
+  ! give the number of dimensions of "Field_name" stored in the input file
   USE netcdf, only: nf90_inq_varid, nf90_inquire_variable, &
                     NF90_NOERR, nf90_strerror
@@ -118 +119 @@
     INTEGER :: varid
     INTEGER :: ierr
-    
+
     IF (is_master) THEN
       ierr=nf90_inq_varid(nid_start,Field_name,varid)
@@ -145 +146 @@
     INTEGER :: varid
     INTEGER :: ierr
-    
+
     IF (is_master) THEN
       ierr=NF90_INQ_DIMID(nid_start,Field_name,varid)
@@ -169 +170 @@
     INTEGER :: varid
     INTEGER :: ierr
-    
+
     IF (is_master) THEN
       ierr=nf90_inq_dimid(nid_start,Field_name,varid)
@@ -194 +195 @@
     CHARACTER(LEN=*),INTENT(IN)    :: Field_name
     REAL,INTENT(INOUT)               :: Field(:)
-    LOGICAL,INTENT(OUT),OPTIONAL   :: found 
+    LOGICAL,INTENT(OUT),OPTIONAL   :: found
     INTEGER,INTENT(IN),OPTIONAL    :: timeindex ! time index of sought data
 
@@ -211 +212 @@
       CALL Get_field_rgen(field_name,field,1,corners,edges)
     ENDIF
-      
+
   END SUBROUTINE Get_Field_r1
-  
+
   SUBROUTINE Get_Field_r2(field_name,field,found,timeindex)
   ! For a "3D" horizontal-vertical field
@@ -220 +221 @@
     CHARACTER(LEN=*),INTENT(IN)    :: Field_name
     REAL,INTENT(INOUT)               :: Field(:,:)
-    LOGICAL,INTENT(OUT),OPTIONAL   :: found 
+    LOGICAL,INTENT(OUT),OPTIONAL   :: found
     INTEGER,INTENT(IN),OPTIONAL    :: timeindex ! time index of sought data
 
@@ -232 +233 @@
       corners(3)=timeindex
     endif
-    
+
     IF (PRESENT(found)) THEN
       CALL Get_field_rgen(field_name,field,size(field,2),&
@@ -241 +242 @@
     ENDIF
 
-      
+
   END SUBROUTINE Get_Field_r2
-  
+
   SUBROUTINE Get_Field_r3(field_name,field,found,timeindex)
   ! for a "4D" field surf/alt/??
@@ -250 +251 @@
     CHARACTER(LEN=*),INTENT(IN)    :: Field_name
     REAL,INTENT(INOUT)               :: Field(:,:,:)
-    LOGICAL,INTENT(OUT),OPTIONAL   :: found 
+    LOGICAL,INTENT(OUT),OPTIONAL   :: found
     INTEGER,INTENT(IN),OPTIONAL    :: timeindex ! time index of sought data
 
@@ -263 +264 @@
       corners(4)=timeindex
     endif
-    
+
     IF (PRESENT(found)) THEN
       CALL Get_field_rgen(field_name,field,size(field,2)*size(field,3),&
@@ -271 +272 @@
                           corners,edges)
     ENDIF
-      
+
   END SUBROUTINE Get_Field_r3
-  
+
   SUBROUTINE Get_field_rgen(field_name,field,field_size, &
                             corners,edges,found)
@@ -287 +288 @@
     INTEGER,INTENT(IN) :: edges(4)
     LOGICAL,OPTIONAL :: found
-    
+
     REAL    :: field_glo(klon_glo,field_size)
     LOGICAL :: tmp_found
     INTEGER :: varid
     INTEGER :: ierr
-    
+
     IF (is_master) THEN
-  
+
       ierr=NF90_INQ_VARID(nid_start,Field_name,varid)
-      
+
       IF (ierr==NF90_NOERR) THEN
         CALL body(field_glo)
@@ -303 +304 @@
         tmp_found=.FALSE.
       ENDIF
-    
-    ENDIF
-    
+
+    ENDIF
+
     CALL bcast(tmp_found)
 
@@ -311 +312 @@
       CALL scatter(field_glo,field)
     ENDIF
-    
+
     IF (PRESENT(found)) THEN
       found=tmp_found
@@ -320 +321 @@
       ENDIF
     ENDIF
- 
-    
+
+
     CONTAINS
-     
+
      SUBROUTINE body(field_glo)
        REAL :: field_glo(klon_glo*field_size)
          ierr=NF90_GET_VAR(nid_start,varid,field_glo,corners,edges)
          IF (ierr/=NF90_NOERR) THEN
-           ! La variable exist dans le fichier mais la lecture a echouee. 
+           ! La variable exist dans le fichier mais la lecture a echouee.
            write(*,*) 'get_field_rgen: Failed reading <'//field_name//'>'
 
@@ -353 +354 @@
   SUBROUTINE get_var_r0(var_name,var,found)
   ! Get a scalar from input file
-  IMPLICIT NONE  
+  IMPLICIT NONE
     CHARACTER(LEN=*),INTENT(IN)  :: var_name
     REAL,INTENT(INOUT)             :: var
@@ -359 +360 @@
 
     REAL                         :: varout(1)
-    
+
     IF (PRESENT(found)) THEN
       CALL Get_var_rgen(var_name,varout,size(varout),found)
@@ -366 +367 @@
     ENDIF
     var=varout(1)
- 
+
   END SUBROUTINE get_var_r0
 
   SUBROUTINE get_var_r1(var_name,var,found)
   ! Get a vector from input file
-  IMPLICIT NONE  
+  IMPLICIT NONE
     CHARACTER(LEN=*),INTENT(IN)  :: var_name
     REAL,INTENT(INOUT)             :: var(:)
     LOGICAL,OPTIONAL,INTENT(OUT) :: found
-    
+
     IF (PRESENT(found)) THEN
       CALL Get_var_rgen(var_name,var,size(var),found)
@@ -381 +382 @@
       CALL Get_var_rgen(var_name,var,size(var))
     ENDIF
-  
+
   END SUBROUTINE get_var_r1
 
   SUBROUTINE get_var_r2(var_name,var,found)
   ! Get a 2D field from input file
-  IMPLICIT NONE  
+  IMPLICIT NONE
     CHARACTER(LEN=*),INTENT(IN)  :: var_name
     REAL,INTENT(OUT)             :: var(:,:)
     LOGICAL,OPTIONAL,INTENT(OUT) :: found
-    
+
     IF (PRESENT(found)) THEN
       CALL Get_var_rgen(var_name,var,size(var),found)
@@ -396 +397 @@
       CALL Get_var_rgen(var_name,var,size(var))
     ENDIF
-  
+
   END SUBROUTINE get_var_r2
 
   SUBROUTINE get_var_r3(var_name,var,found)
   ! Get a 3D field frominput file
-  IMPLICIT NONE  
+  IMPLICIT NONE
     CHARACTER(LEN=*),INTENT(IN)  :: var_name
     REAL,INTENT(INOUT)             :: var(:,:,:)
     LOGICAL,OPTIONAL,INTENT(OUT) :: found
-    
+
     IF (PRESENT(found)) THEN
       CALL Get_var_rgen(var_name,var,size(var),found)
@@ -411 +412 @@
       CALL Get_var_rgen(var_name,var,size(var))
     ENDIF
-  
+
   END SUBROUTINE get_var_r3
 
@@ -424 +425 @@
     REAL             :: var(var_size)
     LOGICAL,OPTIONAL :: found
-    
+
     LOGICAL :: tmp_found
     INTEGER :: varid
     INTEGER :: ierr
-    
+
     IF (is_mpi_root .AND. is_omp_root) THEN
-  
+
       ierr=NF90_INQ_VARID(nid_start,var_name,varid)
-      
+
       IF (ierr==NF90_NOERR) THEN
         ierr=NF90_GET_VAR(nid_start,varid,var)
@@ -443 +444 @@
         tmp_found=.FALSE.
       ENDIF
-    
-    ENDIF
-    
+
+    ENDIF
+
     CALL bcast(tmp_found)
 
@@ -451 +452 @@
       CALL bcast(var)
     ENDIF
-    
+
     IF (PRESENT(found)) THEN
       found=tmp_found
@@ -481 +482 @@
   USE comsoil_h_PEM, only:  nsoilmx_PEM
   USE comslope_mod, only: nslope
+  use layering_mod, only: nb_str_max
   IMPLICIT NONE
     CHARACTER(LEN=*),INTENT(IN) :: filename
@@ -490 +492 @@
     nqmx=nqtot
 #endif
-    
+
     IF (is_master) THEN
       if (.not.already_created) then
@@ -526 +528 @@
         CALL abort_physic("open_restartphy","Failed defining index",1)
       ENDIF
-      
+
       ierr=NF90_DEF_DIM(nid_restart,"physical_points",klon_glo,idim2)
       IF (ierr/=NF90_NOERR) THEN
@@ -533 +535 @@
         CALL abort_physic("open_restartphy","Failed defining physical_points",1)
       ENDIF
-      
+
       ierr=NF90_DEF_DIM(nid_restart,"subsurface_layers",nsoilmx_PEM,idim3)
       IF (ierr/=NF90_NOERR) THEN
@@ -540 +542 @@
         CALL abort_physic("open_restartphy","Failed defining subsurface_layers",1)
       ENDIF
-      
+
       ierr=NF90_DEF_DIM(nid_restart,"nlayer_plus_1",klevp1,idim4)
       IF (ierr/=NF90_NOERR) THEN
@@ -547 +549 @@
         CALL abort_physic("open_restartphy","Failed defining nlayer_plus_1",1)
       ENDIF
-      
+
       if (nqmx.ne.0) then
         ierr=NF90_DEF_DIM(nid_restart,"number_of_advected_fields",nqmx,idim5)
@@ -566 +568 @@
         CALL abort_physic("open_restartphy","Failed defining nlayer",1)
       ENDIF
-      
+
       ierr=NF90_DEF_DIM(nid_restart,"Time",NF90_UNLIMITED,idim7)
       IF (ierr/=NF90_NOERR) THEN
@@ -588 +590 @@
       ENDIF
 
+      ierr=NF90_DEF_DIM(nid_restart,"nb_str_max",nb_str_max,idim10)
+      IF (ierr/=NF90_NOERR) THEN
+        write(*,*)'phyredem: problem defining nb_str_max dimension'
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL ABORT
+      ENDIF
 
       ierr=NF90_ENDDEF(nid_restart)
@@ -608 +616 @@
       ierr = NF90_CLOSE (nid_restart)
     ENDIF
- 
+
   END SUBROUTINE close_restartphy
 
@@ -618 +626 @@
   REAL,INTENT(IN)                :: field(:)
   REAL,OPTIONAL,INTENT(IN)       :: time
-  
+
   IF (present(time)) THEN
     ! if timeindex is present, it is a time-dependent variable
@@ -625 +633 @@
     CALL put_field_rgen(field_name,title,field,1)
   ENDIF
-  
+
   END SUBROUTINE put_field_r1
 
@@ -635 +643 @@
   REAL,INTENT(IN)                :: field(:,:)
   REAL,OPTIONAL,INTENT(IN)       :: time
-  
+
   IF (present(time)) THEN
     ! if timeindex is present, it is a time-dependent variable
@@ -642 +650 @@
     CALL put_field_rgen(field_name,title,field,size(field,2))
   ENDIF
-  
+
   END SUBROUTINE put_field_r2
 
@@ -652 +660 @@
   REAL,INTENT(IN)                :: field(:,:,:)
   REAL,OPTIONAL,INTENT(IN)       :: time
-  
+
   IF (present(time)) THEN
     ! if timeindex is present, it is a time-dependent variable
     CALL put_field_rgen(field_name,title,field,size(field,2)*size(field,3),&
                         time)
-  ELSE  
+  ELSE
     CALL put_field_rgen(field_name,title,field,size(field,2)*size(field,3))
   ENDIF
-  
+
   END SUBROUTINE put_field_r3
-  
+
   SUBROUTINE put_field_rgen(field_name,title,field,field_size,time)
   USE netcdf
@@ -668 +676 @@
   USE comsoil_h_PEM, only:  nsoilmx_PEM
   USE comslope_mod, ONLY: nslope
+  use layering_mod, only: nb_str_max
   USE mod_grid_phy_lmdz
   USE mod_phys_lmdz_para
@@ -676 +685 @@
   REAL,INTENT(IN)                :: field(klon,field_size)
   REAL,OPTIONAL,INTENT(IN)       :: time
-  
+
   REAL                           :: field_glo(klon_glo,field_size)
   INTEGER                        :: ierr
   INTEGER                        :: nvarid
-   
+
     CALL gather(field,field_glo)
-    
+
     IF (is_master) THEN
 
@@ -898 +907 @@
         endif ! of if (.not.present(time))
 
+      ELSE IF (field_size == nb_str_max) THEN
+        ! input is a 2D "stratification" array
+        if (.not.present(time)) then ! for a time-independent field
+          ierr = NF90_REDEF(nid_restart)
+#ifdef NC_DOUBLE
+          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                            (/idim2,idim8/),nvarid)
+#else
+          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
+                            (/idim2,idim8/),nvarid)
+#endif
+          if (ierr.ne.NF90_NOERR) then
+            write(*,*)"put_field_rgen error: failed to define"//trim(field_name)
+            write(*,*)trim(nf90_strerror(ierr))
+          endif
+          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
+          ierr = NF90_ENDDEF(nid_restart)
+          ierr = NF90_PUT_VAR(nid_restart,nvarid,field_glo)
+        else
+          ! check if the variable has already been defined:
+          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
+          if (ierr/=NF90_NOERR) then ! variable not found, define it
+            ierr=NF90_REDEF(nid_restart)
+#ifdef NC_DOUBLE
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                              (/idim2,idim10,idim7/),nvarid)
+#else
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
+                              (/idim2,idim10,idim7/),nvarid)
+#endif
+           if (ierr.ne.NF90_NOERR) then
+              write(*,*)"put_field_rgen error: failed to define"//trim(field_name)
+              write(*,*)trim(nf90_strerror(ierr))
+            endif
+            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
+            ierr=NF90_ENDDEF(nid_restart)
+          endif
+          ! Write the variable
+          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
+                            start=(/1,1,timeindex/))
+
+        endif ! of if (.not.present(time))
+
       ELSE
         write(*,*) "Error phyredem(put_field_rgen) : wrong dimension for ",trim(field_name)
@@ -912 +964 @@
 
     ENDIF ! of IF (is_master)
-    
-  END SUBROUTINE put_field_rgen  
-  
+
+  END SUBROUTINE put_field_rgen
+
   SUBROUTINE put_var_r0(var_name,title,var)
   ! Put a scalar in file
@@ -922 +974 @@
      REAL,INTENT(IN)             :: var
      REAL                        :: varin(1)
-     
+
      varin(1)=var
-     
+
      CALL put_var_rgen(var_name,title,varin,size(varin))
 
@@ -936 +988 @@
      CHARACTER(LEN=*),INTENT(IN) :: title
      REAL,INTENT(IN)             :: var(:)
-     
+
      CALL put_var_rgen(var_name,title,var,size(var))
 
   END SUBROUTINE put_var_r1
- 
+
   SUBROUTINE put_var_r2(var_name,title,var)
   ! Put a 2D field in file
@@ -947 +999 @@
      CHARACTER(LEN=*),INTENT(IN) :: title
      REAL,INTENT(IN)             :: var(:,:)
-     
+
      CALL put_var_rgen(var_name,title,var,size(var))
 
-  END SUBROUTINE put_var_r2     
-  
+  END SUBROUTINE put_var_r2
+
   SUBROUTINE put_var_r3(var_name,title,var)
   ! Put a 3D field in file
@@ -958 +1010 @@
      CHARACTER(LEN=*),INTENT(IN) :: title
      REAL,INTENT(IN)             :: var(:,:,:)
-     
+
      CALL put_var_rgen(var_name,title,var,size(var))
 
@@ -976 +1028 @@
      INTEGER,INTENT(IN)          :: var_size
      REAL,INTENT(IN)             :: var(var_size)
-     
+
      INTEGER :: ierr
      INTEGER :: nvarid
      INTEGER :: idim1d
      logical,save :: firsttime=.true.
-         
+
     IF (is_master) THEN
 
@@ -998 +1050 @@
           IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
           ierr=NF90_ENDDEF(nid_restart)
-          
+
           firsttime=.false.
         endif
@@ -1019 +1071 @@
         idim1d=idim1
       ELSEIF (var_size==nsoilmx_PEM) THEN
-        ! We know it is an  "mlayer" kind of 1D array
+        ! We know it is an "mlayer" kind of 1D array
         idim1d=idim3
       ELSEIF (var_size==nslope+1) THEN
-        ! We know it is an  "inter slope" kind of 1D array
+        ! We know it is an "inter slope" kind of 1D array
         idim1d=idim9
-      ELSE 
+      ELSEIF (var_name == "nb_str_max") THEN
+        ! We know it is a kind of stratification array
+        idim1d = idim10
+      ELSE
         write(*,*) "put_var_rgen error : wrong dimension"
         write(*,*) "  var_size =",var_size
@@ -1051 +1106 @@
       ENDIF
     ENDIF ! of IF (is_master)
-    
-  END SUBROUTINE put_var_rgen     
+
+  END SUBROUTINE put_var_rgen
 
 END MODULE iostart_PEM

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

@@ -14 +14 @@
 ! Numerical parameter
 real, parameter :: eps = epsilon(1.), tol = 1.e2*eps
+integer         :: nb_str_max
 
 ! Physical parameters
@@ -29 +30 @@
 type :: stratum
     real               :: thickness          ! Layer thickness [m]
-    real               :: elevation          ! Layer elevation (top height from the surface) [m]
+    real               :: top_elevation      ! Layer top_elevation (top height from the surface) [m]
     real               :: co2ice_volfrac     ! CO2 ice volumic fraction
     real               :: h2oice_volfrac     ! H2O ice volumic fraction
@@ -45 +46 @@
 end type layering
 
+! Array of pointers to "stratum"
+type ptrarray
+    type(stratum), pointer :: p
+end type ptrarray
+
 !=======================================================================
 contains
+! Procedures to manipulate the data structure:
+!     > print_stratum
+!     > add_stratum
+!     > insert_stratum
+!     > rm_stratum
+!     > get_stratum
+!     > ini_layering
+!     > del_layering
+!     > print_layering
+!     > get_nb_str_max
+!     > stratif2array
+!     > array2stratif
+!     > print_stratif
+! Procedures for the algorithm to evolve the layering:
+!     > make_layering
+!     > subl_co2ice_layering
+!     > subl_h2oice_layering
+!     > lift_dust_lags
 !=======================================================================
 ! To display a stratum
@@ -58 +82 @@
 !---- Code
 write(*,'(a20,es13.6)') 'Thickness      : ', this%thickness
-write(*,'(a20,es13.6)') 'Elevation      : ', this%elevation
-write(*,'(a20,es13.6)') 'CO2 ice content: ', this%co2ice_volfrac
-write(*,'(a20,es13.6)') 'H2O ice content: ', this%h2oice_volfrac
-write(*,'(a20,es13.6)') 'Dust content   : ', this%dust_volfrac
-write(*,'(a20,es13.6)') 'Air content    : ', this%air_volfrac
+write(*,'(a20,es13.6)') 'Top elevation  : ', this%top_elevation
+write(*,'(a20,es13.6)') 'CO2 ice (m3/m3): ', this%co2ice_volfrac
+write(*,'(a20,es13.6)') 'H2O ice (m3/m3): ', this%h2oice_volfrac
+write(*,'(a20,es13.6)') 'Dust (m3/m3)   : ', this%dust_volfrac
+write(*,'(a20,es13.6)') 'Air (m3/m3)    : ', this%air_volfrac
 
 END SUBROUTINE print_stratum
 
 !=======================================================================
-! To display a layering
-SUBROUTINE print_lay(this)
-
-implicit none
-
-!---- Arguments
-type(layering), intent(in) :: this
-
-!---- Local variables
-type(stratum), pointer :: current
-integer                :: i
-
-!---- Code
-current => this%top
-
-i = this%nb_str
-do while (associated(current))
-    write(*,'(a8,i4)') 'Stratum ', i
-    call print_stratum(current)
-    write(*,*)
-    current => current%down
-    i = i - 1
-enddo
-
-END SUBROUTINE print_lay
-
-!=======================================================================
 ! To add a stratum to the top of a layering
-SUBROUTINE add_stratum(this,thickness,elevation,co2ice,h2oice,dust,air)
+SUBROUTINE add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,air)
 
 implicit none
@@ -101 +98 @@
 !---- Arguments
 type(layering), intent(inout)  :: this
-real, intent(in), optional :: thickness, elevation, co2ice, h2oice, dust, air
+real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, air
 
 !---- Local variables
@@ -111 +108 @@
 nullify(str%up,str%down)
 str%thickness = 1.
-str%elevation = 1.
+str%top_elevation = 1.
 str%co2ice_volfrac = 0.
 str%h2oice_volfrac = 0.
@@ -117 +114 @@
 str%air_volfrac = 1.
 if (present(thickness)) str%thickness = thickness
-if (present(elevation)) str%elevation = elevation
+if (present(top_elevation)) str%top_elevation = top_elevation
 if (present(co2ice)) str%co2ice_volfrac = co2ice
 if (present(h2oice)) str%h2oice_volfrac = h2oice
@@ -144 +141 @@
 !=======================================================================
 ! To insert a stratum after another one in a layering
-SUBROUTINE insert_stratum(this,str_prev,thickness,elevation,co2ice,h2oice,dust,air)
+SUBROUTINE insert_stratum(this,str_prev,thickness,top_elevation,co2ice,h2oice,dust,air)
 
 implicit none
@@ -151 +148 @@
 type(layering),         intent(inout) :: this
 type(stratum), pointer, intent(inout) :: str_prev
-real, intent(in), optional :: thickness, elevation, co2ice, h2oice, dust, air
+real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, air
 
 !---- Local variables
@@ -158 +155 @@
 !---- Code
 if (.not. associated(str_prev%up)) then ! If str_prev is the top stratum of the layering
-    call add_stratum(this,thickness,elevation,co2ice,h2oice,dust,air)
+    call add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,air)
 else
     ! Creation of the stratum
@@ -164 +161 @@
     nullify(str%up,str%down)
     str%thickness = 1.
-    str%elevation = 1.
+    str%top_elevation = 1.
     str%co2ice_volfrac = 0.
     str%h2oice_volfrac = 0.
@@ -170 +167 @@
     str%air_volfrac = 1.
     if (present(thickness)) str%thickness = thickness
-    if (present(elevation)) str%elevation = elevation
+    if (present(top_elevation)) str%top_elevation = top_elevation
     if (present(co2ice)) str%co2ice_volfrac = co2ice
     if (present(h2oice)) str%h2oice_volfrac = h2oice
@@ -189 +186 @@
     str%up%down => str
 
-    ! Correct the value of elevation for the upper strata
+    ! Correct the value of top_elevation for the upper strata
     current => str%up
     do while (associated(current))
-        current%elevation = current%down%elevation + current%thickness
+        current%top_elevation = current%down%top_elevation + current%thickness
         current => current%up
     enddo
@@ -214 +211 @@
 !---- Code
 ! Protect the 3 sub-surface strata from removing
-if (str%elevation <= 0.) return
+if (str%top_elevation <= 0.) return
 
 ! Decrement the number of strata
 this%nb_str = this%nb_str - 1
 
-! Remove the stratum from the stratification
+! Remove the stratum from the layering
 str%down%up => str%up
 if (associated(str%up)) then ! If it is not the last one at the top of the layering
@@ -235 +232 @@
 
 !=======================================================================
+! To get a specific stratum in a layering
+FUNCTION get_stratum(lay,i) RESULT(str)
+
+implicit none
+
+!---- Arguments
+type(layering), intent(in) :: lay
+integer,        intent(in) :: i
+type(stratum), pointer     :: str
+
+!---- Local variables
+integer :: k
+
+!---- Code
+if (i < 1 .or. i > lay%nb_str) error stop 'get_stratum: requested number is beyond the number of strata of the layering!'
+k = 1
+str => lay%bottom
+do while (k /= i)
+    str => str%up
+    k = k + 1
+enddo
+
+END FUNCTION get_stratum
+
+!=======================================================================
 ! To initialize the layering
 SUBROUTINE ini_layering(this)
@@ -279 +301 @@
 
 !=======================================================================
-! To get a specific stratum in a layering
-FUNCTION get_stratum(lay,i) RESULT(str)
-
-implicit none
-
-!---- Arguments
-type(layering), intent(in) :: lay
-integer,        intent(in) :: i
-type(stratum), pointer     :: str
-
-!---- Local variables
-integer :: k
-
-!---- Code
-if (i < 1 .or. i > lay%nb_str) error stop 'get_stratum: requested number is beyond the number of strata of the layering!'
-k = 1
-str => lay%bottom
-do while (k /= i)
-    str => str%up
-    k = k + 1
-enddo
-
-END FUNCTION get_stratum
-
-!=======================================================================
-! To sublimate CO2 ice in the layering
-SUBROUTINE subl_co2ice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag)
-
-implicit none
-
-!---- Arguments
-type(layering),         intent(inout) :: this
-type(stratum), pointer, intent(inout) :: str
-logical,                intent(inout) :: new_lag
-real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old
-
-!---- Local variables
-real               :: r_subl, hsubl_dust, h_lag
+! To display a layering
+SUBROUTINE print_layering(this)
+
+implicit none
+
+!---- Arguments
+type(layering), intent(in) :: this
+
+!---- Local variables
 type(stratum), pointer :: current
-
-!---- Code
-! How much dust does CO2 ice sublimation involve to create the lag?
-r_subl = h2subl/(1. - co2ice_porosity) &
-        /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity)))
-hsubl_dust = r_subl*str%dust_volfrac*str%thickness
-! Update of properties in the sublimating stratum
-str%thickness = str%thickness - h2subl/(1. - co2ice_porosity) - hsubl_dust
-str%elevation = str%elevation - h2subl/(1. - co2ice_porosity) - hsubl_dust
-str%co2ice_volfrac = (h_co2ice_old - h2subl)/str%thickness
-str%h2oice_volfrac = h_h2oice_old/str%thickness
-str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
-str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
-! Correct the value of elevation for the upper strata
-current => str%up
+integer                :: i
+
+!---- Code
+current => this%top
+
+i = this%nb_str
 do while (associated(current))
-    current%elevation = current%down%elevation + current%thickness
-    current => current%up
-enddo
-! Remove the sublimating stratum if there is no ice anymore
-if (str%thickness < tol) call rm_stratum(this,str)
-! New stratum = dust lag
-h_lag = hsubl_dust/(1. - dry_lag_porosity)
-if (h_lag > 0.) then ! Only if the dust lag is building up
-    if (new_lag) then
-        call insert_stratum(this,str,h_lag,str%elevation + h_lag,0.,0.,1. - dry_lag_porosity,dry_lag_porosity)
-        new_lag = .false.
-    else
-        str%up%thickness = str%up%thickness + h_lag
-        str%up%elevation = str%up%elevation + h_lag
-    endif
-endif
-
-END SUBROUTINE subl_co2ice_layering
-
-!=======================================================================
-! To sublimate H2O ice in the layering
-SUBROUTINE subl_h2oice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag)
-
-implicit none
-
-!---- Arguments
-type(layering),         intent(inout) :: this
-type(stratum), pointer, intent(inout) :: str
-logical,                intent(inout) :: new_lag
-real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old
-
-!---- Local variables
-real               :: r_subl, hsubl_dust, h_lag
+    write(*,'(a8,i4)') 'Stratum ', i
+    call print_stratum(current)
+    write(*,*)
+    current => current%down
+    i = i - 1
+enddo
+
+END SUBROUTINE print_layering
+
+!=======================================================================
+! To get the maximum number of "stratum" in the stratification (layerings)
+FUNCTION get_nb_str_max(stratif,ngrid,nslope) RESULT(nb_str_max)
+
+implicit none
+
+!---- Arguments
+type(layering), dimension(ngrid,nslope), intent(in) :: stratif
+integer,                                 intent(in) :: ngrid, nslope
+integer :: nb_str_max
+
+!---- Local variables
+integer :: ig, islope
+
+!---- Code
+nb_str_max = 0
+do islope = 1,nslope
+    do ig = 1,ngrid
+        nb_str_max = max(stratif(ig,islope)%nb_str,nb_str_max)
+    enddo
+enddo
+
+END FUNCTION get_nb_str_max
+
+!=======================================================================
+! To convert the stratification data structure (layerings) into an array able to be outputted
+SUBROUTINE stratif2array(stratif,ngrid,nslope,stratif_array)
+
+implicit none
+
+!---- Arguments
+integer,                                 intent(in) :: ngrid, nslope
+type(layering), dimension(ngrid,nslope), intent(in) :: stratif
+real, dimension(:,:,:,:), allocatable, intent(inout) :: stratif_array
+
+!---- Local variables
 type(stratum), pointer :: current
-
-!---- Code
-! How much dust does CO2 ice sublimation involve to create the lag?
-r_subl = h2subl/(1. - h2oice_porosity) &
-        /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity)))
-hsubl_dust = r_subl*str%dust_volfrac*str%thickness
-! Update of properties in the sublimating stratum
-str%thickness = str%thickness - h2subl/(1. - h2oice_porosity) - hsubl_dust
-str%elevation = str%elevation - h2subl/(1. - h2oice_porosity) - hsubl_dust
-str%co2ice_volfrac = h_co2ice_old/str%thickness
-str%h2oice_volfrac = (h_h2oice_old - h2subl)/str%thickness
-str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
-str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
-! Correct the value of elevation for the upper strata
-current => str%up
-do while (associated(current))
-    current%elevation = current%down%elevation + current%thickness
-    current => current%up
-enddo
-! Remove the sublimating stratum if there is no ice anymore
-if (str%thickness < tol) call rm_stratum(this,str)
-! New stratum = dust lag
-h_lag = hsubl_dust/(1. - dry_lag_porosity)
-if (h_lag > 0.) then ! Only if the dust lag is building up
-    if (new_lag) then
-        call insert_stratum(this,str,h_lag,str%elevation + h_lag,0.,0.,1. - dry_lag_porosity,dry_lag_porosity)
-        new_lag = .false.
-    else
-        str%up%thickness = str%up%thickness + h_lag
-        str%up%elevation = str%up%elevation + h_lag
-    endif
-endif
-
-END SUBROUTINE subl_h2oice_layering
-
-!=======================================================================
-! To lift dust in the layering
-SUBROUTINE lift_dust_lags(this,str,h2lift)
-
-implicit none
-
-!---- Arguments
-type(layering),         intent(inout) :: this
-type(stratum), pointer, intent(inout) :: str
-real, intent(in) :: h2lift
-
-!---- Code
-! Update of properties in the eroding dust lag
-str%thickness = str%thickness - h2lift/(1. - str%air_volfrac)
-str%elevation = str%elevation - h2lift/(1. - str%air_volfrac)
-! Remove the eroding dust lag if there is no dust anymore
-if (str%thickness < tol) call rm_stratum(this,str)
-
-END SUBROUTINE lift_dust_lags
+integer                :: ig, islope, k
+
+!---- Code
+stratif_array = 0.
+do islope = 1,nslope
+    do ig = 1,ngrid
+        current => stratif(ig,islope)%bottom
+        k = 1
+        do while (associated(current))
+            stratif_array(ig,islope,k,1) = current%thickness
+            stratif_array(ig,islope,k,2) = current%top_elevation
+            stratif_array(ig,islope,k,3) = current%co2ice_volfrac
+            stratif_array(ig,islope,k,4) = current%h2oice_volfrac
+            stratif_array(ig,islope,k,5) = current%dust_volfrac
+            stratif_array(ig,islope,k,6) = current%air_volfrac
+            current => current%up
+            k = k + 1
+        enddo
+    enddo
+enddo
+
+END SUBROUTINE stratif2array
+
+!=======================================================================
+! To convert the stratification array into the stratification data structure (layerings)
+SUBROUTINE array2stratif(stratif_array,ngrid,nslope,stratif)
+
+implicit none
+
+!---- Arguments
+integer,                               intent(in) :: ngrid, nslope
+real, dimension(:,:,:,:), allocatable, intent(in) :: stratif_array
+type(layering), dimension(ngrid,nslope), intent(inout) :: stratif
+
+!---- Local variables
+type(stratum), pointer :: current
+integer                :: ig, islope, k
+
+!---- Code
+do islope = 1,nslope
+    do ig = 1,ngrid
+        current => stratif(ig,islope)%bottom
+        do k = 1,3
+            current%thickness      = stratif_array(ig,islope,k,1)
+            current%top_elevation  = stratif_array(ig,islope,k,2)
+            current%co2ice_volfrac = stratif_array(ig,islope,k,3)
+            current%h2oice_volfrac = stratif_array(ig,islope,k,4)
+            current%dust_volfrac   = stratif_array(ig,islope,k,5)
+            current%air_volfrac    = stratif_array(ig,islope,k,6)
+            current => current%up
+        enddo
+        do k = 4,size(stratif_array,3)
+            if (all(abs(stratif_array(ig,islope,k,:)) < eps)) exit
+            call add_stratum(stratif(ig,islope),stratif_array(ig,islope,k,1),stratif_array(ig,islope,k,2),stratif_array(ig,islope,k,3),stratif_array(ig,islope,k,4),stratif_array(ig,islope,k,5),stratif_array(ig,islope,k,6))
+        enddo
+    enddo
+enddo
+
+END SUBROUTINE array2stratif
+
+!=======================================================================
+! To display a stratification (layerings)
+SUBROUTINE print_stratif(this,ngrid,nslope)
+
+implicit none
+
+!---- Arguments
+integer,                                 intent(in) :: ngrid, nslope
+type(layering), dimension(ngrid,nslope), intent(in) :: this
+
+!---- Local variables
+integer :: ig, islope
+
+!---- Code
+do ig = 1,ngrid
+    write(*,'(a10,i4)') 'Grid cell ', ig
+    do islope = 1,nslope
+        write(*,'(a13,i1)') 'Slope number ', islope
+        call print_layering(this(ig,islope))
+        write(*,*) '~~~~~~~~~~'
+    enddo
+    write(*,*) '--------------------'
+enddo
+
+END SUBROUTINE print_stratif
 
 !=======================================================================
 ! Layering algorithm
-SUBROUTINE make_layering(this,k,htend_co2ice,htend_h2oice,htend_dust,new_str,new_lag1,new_lag2,stopPEM,current1,current2)
+SUBROUTINE make_layering(this,tend_co2ice,tend_h2oice,tend_dust,new_str,new_lag1,new_lag2,stopPEM,current1,current2)
 
 implicit none
@@ -432 +459 @@
 type(layering),         intent(inout) :: this
 type(stratum), pointer, intent(inout) :: current1, current2
-logical,                intent(inout) :: new_str, new_lag1, new_lag2, stopPEM
-real, intent(in) :: htend_co2ice, htend_h2oice, htend_dust
-integer,  intent(in) :: k
-
-!---- Local variables
-real :: h_co2ice_old, h_h2oice_old, h_dust_old, thickness, h2subl, h2subl_tot, h2lift, h2lift_tot
-
-!---- Code
+logical,                intent(inout) :: new_str, new_lag1, new_lag2
+integer,                intent(inout) :: stopPEM
+real,     intent(in) :: tend_co2ice, tend_h2oice, tend_dust
+
+!---- Local variables
+real :: htend_co2ice, htend_h2oice, htend_dust
+real :: h_co2ice_old, h_h2oice_old, h_dust_old
+real :: thickness, h2subl, h2subl_tot, h2lift, h2lift_tot
+
+!---- Code
+htend_co2ice = tend_co2ice/rho_co2ice
+htend_h2oice = tend_h2oice/rho_h2oice
+htend_dust = tend_dust/rho_dust
+
 if (htend_dust < 0.) then ! Dust lifting
     if (abs(htend_co2ice) > eps .or. abs(htend_h2oice) > eps) error stop 'Situation not managed: dust lifting + CO2/H2O ice condensation/sublimation!'
-    write(*,'(a5,i3,a)') 'Year ', k, ' -> Dust lifting'
+    write(*,'(a)') ' Stratification -> Dust lifting'
     h2lift_tot = abs(htend_dust)
     do while (h2lift_tot > 0. .and. associated(current1))
@@ -460 +493 @@
             endif
         else ! No, we stop
-            write(*,'(a5,i3,a)') 'Year ', k, ' -> Dust lifting: no dust lag!'
-            stopPEM = .true.
+            write(*,'(a)') ' Stratification -> Dust lifting: no dust lag!'
+            stopPEM = 7
             exit
         endif
     enddo
     if (h2lift_tot > 0.) then
-        write(*,'(a5,i3,a)') 'Year ', k, ' -> Not enough dust in the lag layers to complete the lifting!'
-        stopPEM = .true.
+        write(*,'(a)') ' Stratification -> Not enough dust in the lag layers to complete the lifting!'
+        stopPEM = 7
     endif
 
@@ -474 +507 @@
 !------ Dust sedimentation
     if (abs(htend_co2ice) < eps .and. abs(htend_h2oice) < eps) then
-        write(*,'(a5,i3,a)') 'Year ', k, ' -> Dust sedimentation'
+        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 (new_str) then
-                 call add_stratum(this,thickness,this%top%elevation + thickness,0.,0.,htend_dust/thickness,dry_regolith_porosity)
+                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,0.,htend_dust/thickness,dry_regolith_porosity)
                  new_str = .false.
              else
                  this%top%thickness = this%top%thickness + thickness
-                 this%top%elevation = this%top%elevation + thickness
+                 this%top%top_elevation = this%top%top_elevation + thickness
              endif
         endif
@@ -489 +522 @@
 !------ Condensation of CO2 ice + H2O ice
     else if ((htend_co2ice >= 0. .and. htend_h2oice > 0.) .or. (htend_co2ice > 0. .and. htend_h2oice >= 0.)) then
-        write(*,'(a5,i3,a)') 'Year ', k, ' -> CO2 and H2O ice condensation'
+        write(*,'(a)') ' Stratification -> CO2 and H2O ice condensation'
         ! 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 (new_str) then
-                 call add_stratum(this,thickness,this%top%elevation + thickness,htend_co2ice/thickness,htend_h2oice/thickness,htend_dust/thickness,1. - (htend_co2ice/thickness + htend_h2oice/thickness + htend_dust/thickness))
+                 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))
                  new_str = .false.
              else
                  this%top%thickness = this%top%thickness + thickness
-                 this%top%elevation = this%top%elevation + thickness
+                 this%top%top_elevation = this%top%top_elevation + thickness
              endif
         endif
@@ -504 +537 @@
 !------ Sublimation of CO2 ice + Condensation of H2O ice
     else if (htend_co2ice < 0. .and. htend_h2oice >= 0. ) then
-        write(*,'(a5,i3,a)') 'Year ', k, ' -> Sublimation of CO2 ice + Condensation of H2O ice'
+        write(*,'(a)') ' Stratification -> Sublimation of CO2 ice + Condensation of H2O ice'
         ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
         h2subl_tot = abs(htend_co2ice)
@@ -528 +561 @@
         enddo
         if (h2subl_tot > 0.) then
-            write(*,'(a5,i3,a)') 'Year ', k, ' -> Not enough CO2 ice in the layering to complete the sublimation!'
-            stopPEM = .true.
+            write(*,'(a)') ' Stratification -> Not enough CO2 ice in the layering to complete the sublimation!'
+            stopPEM = 7
         endif
         ! New stratum at the layering top by condensation of H2O ice
@@ -535 +568 @@
         if (thickness > 0.) then ! Only if the stratum is builiding up
              if (new_str) then
-                 call add_stratum(this,thickness,this%top%elevation + thickness,0.,htend_h2oice/thickness,htend_dust/thickness,1. - (htend_h2oice/thickness + htend_dust/thickness))
+                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,htend_h2oice/thickness,htend_dust/thickness,1. - (htend_h2oice/thickness + htend_dust/thickness))
                  new_str = .false.
              else
                  this%top%thickness = this%top%thickness + thickness
-                 this%top%elevation = this%top%elevation + thickness
+                 this%top%top_elevation = this%top%top_elevation + thickness
              endif
         endif
@@ -545 +578 @@
 !------ Sublimation of CO2 ice + H2O ice
     else if (htend_co2ice < 0. .and. htend_h2oice < 0.) then
-        write(*,'(a5,i3,a)') 'Year ', k, ' -> Sublimation of CO2 and H2O ice'
+        write(*,'(a)') ' Stratification -> Sublimation of CO2 and H2O ice'
         ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
         h2subl_tot = abs(htend_co2ice)
@@ -569 +602 @@
         enddo
         if (h2subl_tot > 0.) then
-            write(*,'(a5,i3,a)') 'Year ', k, ' -> Not enough CO2 ice in the layering to complete the sublimation!'
-            stopPEM = .true.
+            write(*,'(a)') ' Stratification -> Not enough CO2 ice in the layering to complete the sublimation!'
+            stopPEM = 7
         endif
         ! H2O ice sublimation in the considered stratum + New stratum for dust lag
@@ -595 +628 @@
         enddo
         if (h2subl_tot > 0.) then
-            write(*,'(a5,i3,a)') 'Year ', k, ' -> Not enough H2O ice in the layering to complete the sublimation!'
-            stopPEM = .true.
+            write(*,'(a)') ' Stratification -> Not enough H2O ice in the layering to complete the sublimation!'
+            stopPEM = 7
         endif
 
@@ -607 +640 @@
 END SUBROUTINE make_layering
 
+!=======================================================================
+! To sublimate CO2 ice in the layering
+SUBROUTINE subl_co2ice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag)
+
+implicit none
+
+!---- Arguments
+type(layering),         intent(inout) :: this
+type(stratum), pointer, intent(inout) :: str
+logical,                intent(inout) :: new_lag
+real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old
+
+!---- Local variables
+real                   :: r_subl, hsubl_dust, h_lag
+type(stratum), pointer :: current
+
+!---- Code
+! How much dust does CO2 ice sublimation involve to create the lag?
+r_subl = h2subl/(1. - co2ice_porosity) &
+        /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity)))
+hsubl_dust = r_subl*str%dust_volfrac*str%thickness
+! Update of properties in the sublimating stratum
+str%thickness = str%thickness - h2subl/(1. - co2ice_porosity) - hsubl_dust
+str%top_elevation = str%top_elevation - h2subl/(1. - co2ice_porosity) - hsubl_dust
+str%co2ice_volfrac = (h_co2ice_old - h2subl)/str%thickness
+str%h2oice_volfrac = h_h2oice_old/str%thickness
+str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
+str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
+! Correct the value of top_elevation for the upper strata
+current => str%up
+do while (associated(current))
+    current%top_elevation = current%down%top_elevation + current%thickness
+    current => current%up
+enddo
+! Remove the sublimating stratum if there is no ice anymore
+if (str%thickness < tol) call rm_stratum(this,str)
+! New stratum = dust lag
+h_lag = hsubl_dust/(1. - dry_lag_porosity)
+if (h_lag > 0.) then ! Only if the dust lag is building up
+    if (new_lag) then
+        call insert_stratum(this,str,h_lag,str%top_elevation + h_lag,0.,0.,1. - dry_lag_porosity,dry_lag_porosity)
+        new_lag = .false.
+    else
+        str%up%thickness = str%up%thickness + h_lag
+        str%up%top_elevation = str%up%top_elevation + h_lag
+    endif
+endif
+
+END SUBROUTINE subl_co2ice_layering
+
+!=======================================================================
+! To sublimate H2O ice in the layering
+SUBROUTINE subl_h2oice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag)
+
+implicit none
+
+!---- Arguments
+type(layering),         intent(inout) :: this
+type(stratum), pointer, intent(inout) :: str
+logical,                intent(inout) :: new_lag
+real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old
+
+!---- Local variables
+real                   :: r_subl, hsubl_dust, h_lag
+type(stratum), pointer :: current
+
+!---- Code
+! How much dust does CO2 ice sublimation involve to create the lag?
+r_subl = h2subl/(1. - h2oice_porosity) &
+        /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity)))
+hsubl_dust = r_subl*str%dust_volfrac*str%thickness
+! Update of properties in the sublimating stratum
+str%thickness = str%thickness - h2subl/(1. - h2oice_porosity) - hsubl_dust
+str%top_elevation = str%top_elevation - h2subl/(1. - h2oice_porosity) - hsubl_dust
+str%co2ice_volfrac = h_co2ice_old/str%thickness
+str%h2oice_volfrac = (h_h2oice_old - h2subl)/str%thickness
+str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
+str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
+! Correct the value of top_elevation for the upper strata
+current => str%up
+do while (associated(current))
+    current%top_elevation = current%down%top_elevation + current%thickness
+    current => current%up
+enddo
+! Remove the sublimating stratum if there is no ice anymore
+if (str%thickness < tol) call rm_stratum(this,str)
+! New stratum = dust lag
+h_lag = hsubl_dust/(1. - dry_lag_porosity)
+if (h_lag > 0.) then ! Only if the dust lag is building up
+    if (new_lag) then
+        call insert_stratum(this,str,h_lag,str%top_elevation + h_lag,0.,0.,1. - dry_lag_porosity,dry_lag_porosity)
+        new_lag = .false.
+    else
+        str%up%thickness = str%up%thickness + h_lag
+        str%up%top_elevation = str%up%top_elevation + h_lag
+    endif
+endif
+
+END SUBROUTINE subl_h2oice_layering
+
+!=======================================================================
+! To lift dust in the layering
+SUBROUTINE lift_dust_lags(this,str,h2lift)
+
+implicit none
+
+!---- Arguments
+type(layering),         intent(inout) :: this
+type(stratum), pointer, intent(inout) :: str
+real, intent(in) :: h2lift
+
+!---- Code
+! Update of properties in the eroding dust lag
+str%thickness = str%thickness - h2lift/(1. - str%air_volfrac)
+str%top_elevation = str%top_elevation - h2lift/(1. - str%air_volfrac)
+! Remove the eroding dust lag if there is no dust anymore
+if (str%thickness < tol) call rm_stratum(this,str)
+
+END SUBROUTINE lift_dust_lags
+
 END MODULE layering_mod

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

@@ -46 +46 @@
 use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
                                       TI_PEM,                           & ! soil thermal inertia
-                                      tsoil_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
+                                      tsoil_PEM, layer_PEM, &             ! Soil temp, number of subsurface layers, soil mid layer depths
                                       fluxgeo                             ! Geothermal flux for the PEM and PCM
 use adsorption_mod,             only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
@@ -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
 
 #ifndef CPP_STD
@@ -126 +127 @@
 real               :: time_phys    ! Same as in PCM
 real               :: ptimestep    ! Same as in PCM
-real               :: ztime_fin    ! Same as in PCM
+real               :: ztime_fin     ! Same as in PCM
 
 ! Variables to read start.nc
@@ -183 +184 @@
 real, dimension(:,:),   allocatable :: vmr_co2_PEM_phys ! Physics x Times co2 volume mixing ratio used in the PEM
 real, dimension(:,:),   allocatable :: q_co2_PEM_phys   ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
+
+! Variables for stratification (layering) evolution
+type(layering), dimension(:,:), allocatable :: stratif                     ! Layering (linked list of "stratum") for each grid point and slope
+type(ptrarray), dimension(:,:), allocatable :: current1, current2          ! Current active stratum in the layering
+logical,        dimension(:,:), allocatable :: new_str, new_lag1, new_lag2 ! Flags for the layering algorithm
 
 ! Variables for slopes
@@ -582 +588 @@
 deallocate(min_co2_ice,min_h2o_ice)
 
+allocate(stratif(ngrid,nslope))
+do islope = 1,nslope
+    do i = 1,ngrid
+        call ini_layering(stratif(i,islope))
+    enddo
+enddo
+
 call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
               porefillingice_thickness,tsurf_avg_yr1,tsurf_ave,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries,          &
               tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,global_avg_press_PCM,              &
               watersurf_density_ave,watersoil_density_PEM_ave,co2_adsorbded_phys,delta_co2_adsorbded,                &
-              h2o_adsorbded_phys,delta_h2o_adsorbded)
+              h2o_adsorbded_phys,delta_h2o_adsorbded,stratif)
 
 allocate(co2_ice_ini(ngrid,nslope))
@@ -642 +655 @@
 year_iter = 0
 stopPEM = 0
-
-do while (year_iter < year_iter_max .and. i_myear < n_myear)
+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
+
+do while (year_iter < 10 .and. i_myear < n_myear)
 ! II.a.1. Compute updated global pressure
     write(*,*) "Recomputing the new pressure..."
@@ -725 +748 @@
     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
 
 !------------------------
@@ -1094 +1122 @@
 !    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)
 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, &
-             TI_PEM, porefillingice_depth,porefillingice_thickness,      &
-             co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice)
+             TI_PEM,porefillingice_depth,porefillingice_thickness,       &
+             co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif)
 write(*,*) "restartpem.nc has been written"
 
@@ -1107 +1136 @@
 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)
 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)

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

@@ -10 +10 @@
                     tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,      &
                     global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,              &
-                    m_h2o_regolith_phys,deltam_h2o_regolith_phys)
-
-use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var
+                    m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif)
+
+use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var, inquire_dimension, inquire_dimension_length
 use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, ini_huge_h2oice, depth_breccia, depth_bedrock, index_breccia, index_bedrock
 use comsoil_h,                  only: volcapa, inertiedat
@@ -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
 
 #ifndef CPP_STD
@@ -35 +36 @@
 include "callkeys.h"
 
-character(len=*),               intent(in) :: filename            ! name of the startfi_PEM.nc
+character(*),                   intent(in) :: filename            ! name of the startfi_PEM.nc
 integer,                        intent(in) :: ngrid               ! # of physical grid points
 integer,                        intent(in) :: nsoil_PCM           ! # of vertical grid points in the PCM
@@ -56 +57 @@
 real, dimension(ngrid,nslope),                   intent(out) :: h2o_ice                  ! h2o ice amount [kg/m^2]
 real, dimension(ngrid,nslope),                   intent(out) :: co2_ice                  ! co2 ice amount [kg/m^2]
+type(layering), dimension(ngrid,nslope),         intent(inout) :: stratif             ! stratification (layerings)
 real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
 real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
@@ -65 +67 @@
 real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: watersoil_avg       ! surface water ice density, yearly averaged (kg/m^3)
 ! local
-real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM               ! soil temperature saved in the start [K]
-real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM                  ! soil thermal inertia saved in the start [SI]
-logical                                 :: found                        ! check if variables are found in the start
-logical                                 :: found2                       ! check if variables are found in the start
-integer                                 :: iloop, ig, islope, it, isoil ! index for loops
-real                                    :: kcond                        ! Thermal conductivity, intermediate variable [SI]
-real                                    :: delta                        ! Depth of the interface regolith-breccia, breccia -bedrock [m]
-character(2)                            :: num                          ! intermediate string to read PEM start sloped variables
-real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1                ! intermediate soil temperature during yr 1 [K]
-real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2                ! intermediate soil temperature during yr 2 [K]
-logical                                 :: startpem_file                ! boolean to check if we read the startfile or not
+real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM                  ! soil temperature saved in the start [K]
+real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM                     ! soil thermal inertia saved in the start [SI]
+logical                                 :: found, found2                   ! check if variables are found in the start
+integer                                 :: iloop, ig, islope, it, isoil, k ! index for loops
+real                                    :: kcond                           ! Thermal conductivity, intermediate variable [SI]
+real                                    :: delta                           ! Depth of the interface regolith-breccia, breccia -bedrock [m]
+character(2)                            :: num                             ! intermediate string to read PEM start sloped variables
+real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1                   ! intermediate soil temperature during yr 1 [K]
+real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2                   ! intermediate soil temperature during yr 2 [K]
+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
@@ -112 +115 @@
     ! open pem initial state file:
     call open_startphy(filename)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    ! 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')
+    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)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -320 +355 @@
 else !No startfi, let's build all by hand
 
+    write(*,*)'There is no "startpem.nc"!'
+
+    ! Stratification (layerings)
+    write(*,*)'So the stratification (layerings) is initialized with nothing more than the 3 sub-surface strata.'
+    nb_str_max = 3
+
     ! h2o ice
     h2o_ice = 0.
-    write(*,*)'There is no "startpem.nc"!'
     write(*,*)'So ''h2o_ice'' is initialized with default value ''ini_huge_h2oice'' where ''watercaptag'' is true.'
     do ig = 1,ngrid

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

@@ -58 +58 @@
 
 SUBROUTINE pemdem1(filename,i_myear,nsoil_PEM,ngrid,nslope,tsoil_slope_PEM,inertiesoil_slope_PEM, &
-                   ice_table_depth,ice_table_thickness,m_co2_regolith,m_h2o_regolith,h2o_ice)
+                   ice_table_depth,ice_table_thickness,m_co2_regolith,m_h2o_regolith,h2o_ice,stratif)
 
 ! write time-dependent variable to restart file
@@ -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
 
 implicit none
@@ -79 +80 @@
 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: m_co2_regolith, m_h2o_regolith
 real, dimension(ngrid,nslope),           intent(in) :: h2o_ice
+type(layering), dimension(ngrid,nslope), intent(in) :: stratif ! Stratification (layerings)
 
-integer       :: islope
-character(2)  :: num
-real          :: Year ! Year of the simulation
+integer                               :: islope
+character(2)                          :: num
+real                                  :: Year          ! Year of the simulation
+real, dimension(:,:,:,:), allocatable :: stratif_array ! Array for stratification (layerings)
 
 ! Open file
@@ -92 +95 @@
 call put_var("Time","Year of simulation",Year)
 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 (soil_pem) then