Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3082)
+++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3088)
@@ -101,2 +101,9 @@
     - Correction of a bug in the initialization of constants. The correct modules are now used: 'comcstfi_h' (and no longer 'comconst_mod'!) in the general case and 'comcstfi_mod' in the case of generic model;
     - Addition of the variable 'ecritpem' in "run_PEM.def" to set the frequency of outputs in the "diagfi.nc". By default, 'ecritpem = 1' which means there is one output at each PEM year.
+
+== 16/10/2023 == JBC
+    - Addition of the bash script "concat_diagpem.sh" in deftank to concatenate all the "diagpem.nc" outputs into one file;
+    - Addition of "writediagpem_mod.F90" inspired by the subroutine "writediagfi.F" in phymars/;
+    - Renaming of PEM outputs ('startfi_PEM' -> 'startpem' and 'diagfi' -> 'diagpem');
+    - Gathering of all "diagpem.nc" outputs at the same place at the end of the main PEM loop;
+    - Update of "launch_pem.sh", "README" and "modify_startfi_orbit.sh" in the deftank.
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3082)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3088)
@@ -8,5 +8,5 @@
 !    I_f Compute tendencies & Save initial situation
 !    I_g Save initial PCM situation
-!    I_h Read the PEMstart
+!    I_h Read the startpem.nc
 !    I_i Compute orbit criterion
 
@@ -16,11 +16,12 @@
 !    II_c CO2 & H2O glaciers flows
 !    II_d Update surface and soil temperatures
-!    II_e Update the tendencies
-!    II_f Checking the stopping criterion
+!    II_e Outputs
+!    II_f Update the tendencies
+!    II_g Checking the stopping criterion
 
 ! III Output
 !    III_a Update surface value for the PCM start files
 !    III_b Write restart_evol.nc and restartfi_evol.nc
-!    III_c Write restartfi_PEM.nc
+!    III_c Write restartpem.nc
 !------------------------
 
@@ -69,4 +70,5 @@
 use recomp_tend_co2_slope_mod,    only: recomp_tend_co2_slope
 use soil_pem_compute_mod,         only: soil_pem_compute
+use writediagpem_mod,             only: writediagpem
 
 #ifndef CPP_STD
@@ -598,7 +600,7 @@
 !------------------------
 ! I   Initialization
-!    I_h Read the PEMstart
-!------------------------
-call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
+!    I_h Read the startpem.nc
+!------------------------
+call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
               porefillingice_thickness,tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,          &
               tsoil_phys_PEM_timeseries,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:),                      &
@@ -665,5 +667,4 @@
         enddo
     enddo
-    call writediagfi(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new)
 
     if (adsorption_pem) then
@@ -741,12 +742,4 @@
     call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
 
-    do islope = 1,nslope
-        write(str2(1:2),'(i2.2)') islope
-        call writediagfi(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
-        call writediagfi(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
-        call writediagfi(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
-        call writediagfi(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
-    enddo
-
 !------------------------
 ! II  Run
@@ -754,18 +747,9 @@
 !------------------------
     write(*,*) "CO2 glacier flows"
-
     if (co2glaciersflow) call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries, &
                                                global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
 
     write(*,*) "H2O glacier flows"
-
     if (h2oglaciersflow) call h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,qsurf(:,igcm_h2o_ice,:),flag_h2oflow,flag_h2oflow_mesh)
-
-    do islope = 1,nslope
-        write(str2(1:2),'(i2.2)') islope
-        call writediagfi(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
-        call writediagfi(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
-        call writediagfi(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
-    enddo
 
 !------------------------
@@ -845,9 +829,4 @@
     enddo
 
-    do islope = 1,nslope
-        write(str2(1:2),'(i2.2)') islope
-        call writediagfi(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
-    enddo
-
     if (soil_pem) then
 
@@ -919,5 +898,20 @@
 !------------------------
 ! II  Run
-!    II_e Update the tendencies
+!    II_e Outputs
+!------------------------
+    call writediagpem(ngrid,'ps_ave','Global average pressure','Pa',0,(/global_ave_press_new/))
+    do islope = 1,nslope
+        write(str2(1:2),'(i2.2)') islope
+        call writediagpem(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
+        call writediagpem(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
+        call writediagpem(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
+        call writediagpem(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
+        call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
+        call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
+    enddo
+
+!------------------------
+! II  Run
+!    II_f Update the tendencies
 !------------------------
     write(*,*) "Adaptation of the new co2 tendencies to the current pressure"
@@ -927,5 +921,5 @@
 !------------------------
 ! II  Run
-!    II_f Checking the stopping criterion
+!    II_g Checking the stopping criterion
 !------------------------
     call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice)
@@ -1136,12 +1130,12 @@
 !------------------------
 ! III Output
-!    III_c Write restartfi_PEM.nc
-!------------------------
-call pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
+!    III_c Write restartpem.nc
+!------------------------
+call pemdem0("restartpem.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
              float(day_ini),0.,nslope,def_slope,subslope_dist)
-call pemdem1("restartfi_PEM.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
+call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
              TI_PEM, porefillingice_depth,porefillingice_thickness,         &
              co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
-write(*,*) "restartfi_PEM.nc has been written"
+write(*,*) "restartpem.nc has been written"
 
 call info_run_PEM(year_iter,criterion_stop,i_myear,n_myear)
Index: trunk/LMDZ.COMMON/libf/evolution/writediagpem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/writediagpem.F90	(revision 3088)
+++ trunk/LMDZ.COMMON/libf/evolution/writediagpem.F90	(revision 3088)
@@ -0,0 +1,583 @@
+MODULE writediagpem_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE writediagpem(ngrid,nom,titre,unite,dim,px)
+
+!  Ecriture de variables diagnostiques au choix dans la physique
+!  dans un fichier NetCDF nomme  'diagfi'. Ces variables peuvent etre
+!  3d (ex : temperature), 2d (ex : temperature de surface), ou
+!  0d (pour un scalaire qui ne depend que du temps : ex : la longitude
+!  solaire)
+!  (ou encore 1d, dans le cas de testphys1d, pour sortir une colonne)
+!  La periode d'ecriture est donnee par
+!  "ecritphy " regle dans le fichier de controle de run :  run.def
+!
+!    writediagpem peut etre appele de n'importe quelle subroutine
+!    de la physique, plusieurs fois. L'initialisation et la creation du
+!    fichier se fait au tout premier appel.
+!
+! WARNING : les variables dynamique (u,v,t,q,ps)
+!  sauvees par writediagpem avec une
+! date donnee sont legerement differentes que dans le fichier histoire car
+! on ne leur a pas encore ajoute de la dissipation et de la physique !!!
+! IL est  RECOMMANDE d'ajouter les tendance physique a ces variables
+! avant l'ecriture dans diagfi (cf. physiq.F)
+!
+! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4
+!         Oct 2011 Francois: enable having a 'diagfi.def' file to select
+!                            at runtime, which variables to put in file
+!         Oct 2023 JB: conversion into Fortran 90 with module for the PEM
+!
+!  parametres (input) :
+!  ----------
+!      ngrid : nombres de point ou est calcule la physique
+!                (ngrid = 2+(jjm-1)*iim - 1/jjm)
+!                (= nlon ou klon dans la physique terrestre)
+!      unit : unite logique du fichier de sortie (toujours la meme)
+!      nom  : nom de la variable a sortir (chaine de caracteres)
+!      titre: titre de la variable (chaine de caracteres)
+!      unite : unite de la variable (chaine de caracteres)
+!      px : variable a sortir (real 0, 1, 2, ou 3d)
+!      dim : dimension de px : 0, 1, 2, ou 3 dimensions
+!
+!=================================================================
+use surfdat_h,          only: phisfi
+use geometry_mod,       only: cell_area
+use time_phylmdz_mod,   only: ecritphy, day_step, iphysiq, day_ini
+use mod_phys_lmdz_para, only: is_parallel, is_mpi_root, is_master, gather
+use mod_grid_phy_lmdz,  only: klon_glo, Grid1Dto2D_glo, nbp_lon, nbp_lat, nbp_lev, grid_type, unstructured
+
+implicit none
+
+! Commons
+include "netcdf.inc"
+
+! Arguments on input:
+integer,                        intent(in) :: ngrid
+character(len=*),               intent(in) :: nom, titre, unite
+integer,                        intent(in) :: dim
+real, dimension(ngrid,nbp_lev), intent(in) :: px
+
+! Local variables:
+real*4, dimension(nbp_lon + 1,nbp_lat,nbp_lev) :: dx3    ! to store a 3D data set
+real*4, dimension(nbp_lon + 1,nbp_lat)         :: dx2    ! to store a 2D (surface) data set
+real*4, dimension(nbp_lev)                     :: dx1    ! to store a 1D (column) data set
+real*4                                         :: dx0
+real*4, dimension(1,nbp_lev)                   :: dx3_1d ! to store a profile with 1D model
+real*4                                         :: dx2_1d ! to store a surface value with 1D model
+real*4, save                                   :: date
+
+!$OMP THREADPRIVATE(date)
+real, dimension((nbp_lon + 1),nbp_lat) :: phis
+real, dimension((nbp_lon + 1),nbp_lat) :: area
+integer                                :: irythme, ierr, ierr2, i, j, l, ig0
+integer,                          save :: zitau = 0
+character(27),                    save :: firstnom='1234567890'
+!$OMP THREADPRIVATE(zitau,firstnom)
+
+! Ajouts
+integer,           save :: ntime = 0
+!$OMP THREADPRIVATE(ntime)
+integer                 :: idim, varid
+integer                 :: nid
+character(*), parameter :: fichnom = "diagfi.nc"
+integer, dimension(4)   :: id
+integer, dimension(4)   :: edges, corner
+
+! Added to use diagfi.def to select output variable
+logical,                                  save :: diagfi_def
+logical                                        :: getout
+integer,                                  save :: n_nom_def
+integer                                        :: n
+integer, parameter                             :: n_nom_def_max = 199
+character(120), dimension(n_nom_def_max), save :: nom_def
+logical,                                  save :: firstcall = .true.
+!$OMP THREADPRIVATE(firstcall) !diagfi_def,n_nom_def,nom_def read in diagfi.def
+
+#ifdef CPP_PARA
+! Added to work in parallel mode
+    real, dimension(klon_glo,nbp_lev)        :: dx3_glop
+    real, dimension(nbp_lon,nbp_lat,nbp_lev) :: dx3_glo  ! to store a global 3D data set
+    real, dimension(klon_glo)                :: dx2_glop
+    real, dimension(nbp_lon,nbp_lat)         :: dx2_glo  ! to store a global 2D (surface) data set
+    real, dimension(ngrid)                   :: px2
+!    real, dimension(nbp_lev)                :: dx1_glo   ! to store a 1D (column) data set
+!    real :: dx0_glo
+    real, dimension(klon_glo) :: phisfi_glo ! surface geopotential on global physics grid
+    real, dimension(klon_glo) :: areafi_glo ! mesh area on global physics grid
+#else
+    real, dimension(ngrid) :: phisfi_glo ! surface geopotential on global physics grid
+    real, dimension(ngrid) :: areafi_glo ! mesh area on global physics grid
+#endif
+
+if (grid_type == unstructured) return
+
+!***************************************************************
+!Sortie des variables au rythme voulu
+irythme = int(ecritphy) ! output rate set by ecritphy
+
+!***************************************************************
+! At very first call, check if there is a "diagfi.def" to use and read it
+! -----------------------------------------------------------------------
+IF (firstcall) THEN
+    firstcall=.false.
+
+!$OMP MASTER
+    ! Open diagfi.def definition file if there is one:
+        open(99,file="diagfi.def",status='old',form='formatted',iostat=ierr2)
+
+        if (ierr2 == 0) then
+            diagfi_def=.true.
+            write(*,*) "******************"
+            write(*,*) "Reading diagfi.def"
+            write(*,*) "******************"
+            do n=1,n_nom_def_max
+                read(99,fmt='(a)',end=88) nom_def(n)
+                write(*,*) 'Output in diagfi: ', trim(nom_def(n))
+            enddo
+ 88         continue
+            if (n.ge.n_nom_def_max) then
+                write(*,*)"n_nom_def_max too small in writediagpem.F:",n
+                call abort_physic("writediagpem","n_nom_def_max too small",1)
+            end if
+            n_nom_def=n-1
+            close(99)
+        else
+            diagfi_def=.false.
+        endif
+!$OMP END MASTER
+!$OMP BARRIER
+ENDIF ! of IF (firstcall)
+
+! Get out of write_diagfi if there is diagfi.def AND variable not listed
+!  ---------------------------------------------------------------------
+if (diagfi_def) then
+    getout=.true.
+    do n=1,n_nom_def
+        if (trim(nom_def(n)).eq.nom) getout=.false.
+    enddo
+    if (getout) return
+endif
+
+! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file
+! ------------------------------------------------------------------------
+! (at very first call to the subroutine, in accordance with diagfi.def)
+if (firstnom.eq.'1234567890') then ! .true. for the very first valid
+    ! call to this subroutine; now set 'firstnom'
+    firstnom = nom
+    ! just to be sure, check that firstnom is large enough to hold nom
+    if (len_trim(firstnom).lt.len_trim(nom)) then
+        write(*,*) "writediagpem: Error !!!"
+        write(*,*) "   firstnom string not long enough!!"
+        write(*,*) "   increase its size to at least ",len_trim(nom)
+        call abort_physic("writediagpem","firstnom too short",1)
+    endif
+
+    zitau = -1 ! initialize zitau
+
+#ifdef CPP_PARA
+        ! Gather phisfi() geopotential on physics grid
+        call Gather(phisfi,phisfi_glo)
+        ! Gather cell_area() mesh area on physics grid
+        call Gather(cell_area,areafi_glo)
+#else
+        phisfi_glo(:)=phisfi(:)
+        areafi_glo(:)=cell_area(:)
+#endif
+
+    !! parallel: we cannot use the usual writediagpem method
+    !! call iophys_ini
+    if (is_master) then
+        ! only the master is required to do this
+        ! Create the NetCDF file
+        ierr = NF_CREATE(fichnom,NF_CLOBBER,nid)
+        ! Define the 'Time' dimension
+        ierr = nf_def_dim(nid,"Time",NF_UNLIMITED,idim)
+        ! Define the 'Time' variable
+    !#ifdef NC_DOUBLE
+    !    ierr = NF_DEF_VAR (nid, "Time", NF_DOUBLE, 1, idim,varid)
+    !#else
+        ierr = NF_DEF_VAR (nid, "Time", NF_FLOAT, 1, idim,varid)
+    !#endif
+        ! Add a long_name attribute
+        ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name",4,"Time")
+        ! Add a units attribute
+        ierr = NF_PUT_ATT_TEXT(nid, varid,'units',29,"days since 0000-00-0 00:00:00")
+        ! Switch out of NetCDF Define mode
+        ierr = NF_ENDDEF(nid)
+
+        ! Build phis() and area()
+        IF (klon_glo>1) THEN
+            do i=1,nbp_lon+1 ! poles
+            phis(i,1)=phisfi_glo(1)
+            phis(i,nbp_lat)=phisfi_glo(klon_glo)
+            ! for area, divide at the poles by nbp_lon
+            area(i,1)=areafi_glo(1)/nbp_lon
+            area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
+            enddo
+            do j=2,nbp_lat-1
+            ig0= 1+(j-2)*nbp_lon
+            do i=1,nbp_lon
+                phis(i,j)=phisfi_glo(ig0+i)
+                area(i,j)=areafi_glo(ig0+i)
+            enddo
+            ! handle redundant point in longitude
+            phis(nbp_lon+1,j)=phis(1,j)
+            area(nbp_lon+1,j)=area(1,j)
+            enddo
+        ENDIF
+
+        ! write "header" of file (longitudes, latitudes, geopotential, ...)
+        IF (klon_glo>1) THEN ! general 3D case
+            call iniwrite(nid,day_ini,phis,area,nbp_lon+1,nbp_lat)
+        ELSE ! 1D model
+            call iniwrite(nid,day_ini,phisfi_glo(1),areafi_glo(1),1,1)
+        ENDIF
+
+    endif ! of if (is_master)
+
+    else
+
+    if (is_master) then
+        ! only the master is required to do this
+        ! Open the NetCDF file
+        ierr = NF_OPEN(fichnom,NF_WRITE,nid)
+    endif ! of if (is_master)
+
+endif ! if (firstnom.eq.'1234567890')
+
+! Increment time index 'zitau' if it is the "fist call" (at given time level)
+! to writediagpem
+!------------------------------------------------------------------------
+if (nom == firstnom) zitau = zitau + iphysiq
+
+!--------------------------------------------------------
+! Write the variables to output file if it's time to do so
+!--------------------------------------------------------
+if (MOD(zitau+1,irythme) == 0.) then
+! Compute/write/extend 'Time' coordinate (date given in days)
+! (done every "first call" (at given time level) to writediagpem)
+! Note: date is incremented as 1 step ahead of physics time
+!--------------------------------------------------------
+    if (is_master) then
+        ! only the master is required to do this
+        if (nom.eq.firstnom) then
+        ! We have identified a "first call" (at given date)
+            ntime=ntime+1 ! increment # of stored time steps
+            ! compute corresponding date (in days and fractions thereof)
+            date=(zitau +1.)/day_step
+            ! Get NetCDF ID of 'Time' variable
+            ierr= NF_INQ_VARID(nid,"Time",varid)
+            ! Write (append) the new date to the 'Time' array
+!#ifdef NC_DOUBLE
+!    ierr= NF_PUT_VARA_DOUBLE(nid,varid,[ntime],[1],[date])
+!#else
+            ierr= NF_PUT_VARA_REAL(nid,varid,[ntime],[1],[date])
+!#endif
+            if (ierr.ne.NF_NOERR) then
+                write(*,*) "***** PUT_VAR matter in writediagpem_nc"
+                write(*,*) "***** with time"
+                write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
+!                call abort
+            endif
+
+            write(6,*)'WRITEDIAGFI: date= ', date
+        endif ! of if (nom.eq.firstnom)
+    endif ! of if (is_master)
+
+!Case of a 3D variable
+!---------------------
+    if (dim == 3) then
+#ifdef CPP_PARA
+        ! Gather field on a "global" (without redundant longitude) array
+        call Gather(px,dx3_glop)
+!$OMP MASTER
+        if (is_mpi_root) then
+            call Grid1Dto2D_glo(dx3_glop,dx3_glo)
+            ! copy dx3_glo() to dx3(:) and add redundant longitude
+            dx3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
+            dx3(nbp_lon+1,:,:)=dx3(1,:,:)
+        endif
+!$OMP END MASTER
+!$OMP BARRIER
+#else
+!         Passage variable physique -->  variable dynamique
+!         recast (copy) variable from physics grid to dynamics grid
+        IF (klon_glo>1) THEN ! General case
+            DO l=1,nbp_lev
+                DO i=1,nbp_lon+1
+                    dx3(i,1,l)=px(1,l)
+                    dx3(i,nbp_lat,l)=px(ngrid,l)
+                ENDDO
+                DO j=2,nbp_lat-1
+                    ig0= 1+(j-2)*nbp_lon
+                    DO i=1,nbp_lon
+                        dx3(i,j,l)=px(ig0+i,l)
+                    ENDDO
+                    dx3(nbp_lon+1,j,l)=dx3(1,j,l)
+                ENDDO
+            ENDDO
+        ELSE ! 1D model case
+            dx3_1d(1,1:nbp_lev)=px(1,1:nbp_lev)
+        ENDIF
+#endif
+!       Ecriture du champs
+        if (is_master) then
+! only the master writes to output
+! name of the variable
+            ierr= NF_INQ_VARID(nid,nom,varid)
+            if (ierr /= NF_NOERR) then
+! corresponding dimensions
+                ierr = NF_INQ_DIMID(nid,"longitude",id(1))
+                ierr = NF_INQ_DIMID(nid,"latitude",id(2))
+                ierr = NF_INQ_DIMID(nid,"altitude",id(3))
+                ierr = NF_INQ_DIMID(nid,"Time",id(4))
+
+! Create the variable if it doesn't exist yet
+                write (*,*) "=========================="
+                write (*,*) "DIAGFI: creating variable ",trim(nom)
+                call def_var(nid,nom,titre,unite,4,id,varid,ierr)
+
+            else
+                if (ntime==0) then
+                    write(*,*) "DIAGFI Error: failed creating variable ",trim(nom)
+                    write(*,*) "it seems it already exists!"
+                    call abort_physic("writediagpem",trim(nom)//" already exists",1)
+                endif
+            endif
+
+            corner(1)=1
+            corner(2)=1
+            corner(3)=1
+            corner(4)=ntime
+
+            IF (klon_glo==1) THEN
+                edges(1)=1
+            ELSE
+                edges(1)=nbp_lon+1
+            ENDIF
+            edges(2)=nbp_lat
+            edges(3)=nbp_lev
+            edges(4)=1
+!#ifdef NC_DOUBLE
+!            ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
+!#else
+!            write(*,*)"test:  nid=",nid," varid=",varid
+!            write(*,*)"       corner()=",corner
+!            write(*,*)"       edges()=",edges
+!            write(*,*)"       dx3()=",dx3
+            IF (klon_glo>1) THEN ! General case
+                ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
+            ELSE
+                ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3_1d)
+            ENDIF
+!#endif
+
+            if (ierr.ne.NF_NOERR) then
+                write(*,*) "***** PUT_VAR problem in writediagpem"
+                write(*,*) "***** with dx3: ",trim(nom)
+                write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
+                call abort_physic("writediagpem","failed writing "//trim(nom),1)
+            endif
+
+        endif !of if (is_master)
+
+!Case of a 2D variable
+!---------------------
+    else if (dim == 2) then
+
+#ifdef CPP_PARA
+        ! Gather field on a "global" (without redundant longitude) array
+        px2(:)=px(:,1)
+        call Gather(px2,dx2_glop)
+!$OMP MASTER
+        if (is_mpi_root) then
+            call Grid1Dto2D_glo(dx2_glop,dx2_glo)
+            ! copy dx2_glo() to dx2(:) and add redundant longitude
+            dx2(1:nbp_lon,:)=dx2_glo(1:nbp_lon,:)
+            dx2(nbp_lon+1,:)=dx2(1,:)
+        endif
+!$OMP END MASTER
+!$OMP BARRIER
+#else
+!         Passage variable physique -->  physique dynamique
+!         recast (copy) variable from physics grid to dynamics grid
+        IF (klon_glo>1) THEN ! General case
+            DO i=1,nbp_lon+1
+                dx2(i,1)=px(1,1)
+                dx2(i,nbp_lat)=px(ngrid,1)
+            ENDDO
+            DO j=2,nbp_lat-1
+                ig0= 1+(j-2)*nbp_lon
+                DO i=1,nbp_lon
+                    dx2(i,j)=px(ig0+i,1)
+                ENDDO
+                dx2(nbp_lon+1,j)=dx2(1,j)
+            ENDDO
+        ELSE ! 1D model case
+            dx2_1d=px(1,1)
+        ENDIF
+#endif
+
+        if (is_master) then
+        ! only the master writes to output
+!            write (*,*) 'In  writediagpem, on sauve:  ' , nom
+!            write (*,*) 'In  writediagpem. Estimated date = ' ,date
+            ierr= NF_INQ_VARID(nid,nom,varid)
+            if (ierr /= NF_NOERR) then
+! corresponding dimensions
+                ierr= NF_INQ_DIMID(nid,"longitude",id(1))
+                ierr= NF_INQ_DIMID(nid,"latitude",id(2))
+                ierr= NF_INQ_DIMID(nid,"Time",id(3))
+
+! Create the variable if it doesn't exist yet
+                write (*,*) "=========================="
+                write (*,*) "DIAGFI: creating variable ",trim(nom)
+
+                call def_var(nid,nom,titre,unite,3,id,varid,ierr)
+
+            else
+                if (ntime==0) then
+                    write(*,*) "DIAGFI Error: failed creating variable ",trim(nom)
+                    write(*,*) "it seems it already exists!"
+                    call abort_physic("writediagpem",trim(nom)//" already exists",1)
+                endif
+            endif
+
+            corner(1)=1
+            corner(2)=1
+            corner(3)=ntime
+            IF (klon_glo==1) THEN
+                edges(1)=1
+            ELSE
+                edges(1)=nbp_lon+1
+            ENDIF
+            edges(2)=nbp_lat
+            edges(3)=1
+
+
+!#ifdef NC_DOUBLE
+!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx2)
+!#else
+            IF (klon_glo>1) THEN ! General case
+                ierr = NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2)
+            ELSE
+                ierr = NF_PUT_VARA_REAL(nid,varid,corner,edges,[dx2_1d])
+            ENDIF
+!#endif
+
+            if (ierr.ne.NF_NOERR) then
+                write(*,*) "***** PUT_VAR matter in writediagpem"
+                write(*,*) "***** with dx2: ",trim(nom)
+                write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
+                call abort_physic("writediagpem","failed writing "//trim(nom),1)
+            endif
+
+        endif !of if (is_master)
+
+!Case of a 1D variable (ie: a column)
+!---------------------------------------------------
+    else if (dim.eq.1) then
+        if (is_parallel) then
+            write(*,*) "writediagpem error: dim=1 not implemented ","in parallel mode. Problem for ",trim(nom)
+                call abort_physic("writediagpem","failed writing "//trim(nom),1)
+        endif
+!       Passage variable physique -->  physique dynamique
+!       recast (copy) variable from physics grid to dynamics grid
+        do l=1,nbp_lev
+            dx1(l)=px(1,l)
+        enddo
+
+        ierr= NF_INQ_VARID(nid,nom,varid)
+        if (ierr /= NF_NOERR) then
+! corresponding dimensions
+            ierr= NF_INQ_DIMID(nid,"altitude",id(1))
+            ierr= NF_INQ_DIMID(nid,"Time",id(2))
+
+! Create the variable if it doesn't exist yet
+
+            write (*,*) "=========================="
+            write (*,*) "DIAGFI: creating variable ",trim(nom)
+
+            call def_var(nid,nom,titre,unite,2,id,varid,ierr)
+
+        else
+            if (ntime==0) then
+                write(*,*) "DIAGFI Error: failed creating variable ",trim(nom)
+                write(*,*) "it seems it already exists!"
+                call abort_physic("writediagpem",trim(nom)//" already exists",1)
+            endif
+        endif
+
+        corner(1)=1
+        corner(2)=ntime
+        edges(1)=nbp_lev
+        edges(2)=1
+!#ifdef NC_DOUBLE
+!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx1)
+!#else
+        ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx1)
+!#endif
+
+        if (ierr /= NF_NOERR) then
+            write(*,*) "***** PUT_VAR problem in writediagpem"
+            write(*,*) "***** with dx1: ",trim(nom)
+            write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
+            call abort_physic("writediagpem","failed writing "//trim(nom),1)
+        endif
+
+!Case of a 0D variable (ie: a time-dependent scalar)
+!---------------------------------------------------
+    else if (dim == 0) then
+
+        dx0 = px (1,1)
+
+        if (is_master) then
+            ! only the master writes to output
+            ierr= NF_INQ_VARID(nid,nom,varid)
+            if (ierr /= NF_NOERR) then
+! corresponding dimensions
+                ierr= NF_INQ_DIMID(nid,"Time",id(1))
+
+! Create the variable if it doesn't exist yet
+                write (*,*) "=========================="
+                write (*,*) "DIAGFI: creating variable ",trim(nom)
+                call def_var(nid,nom,titre,unite,1,id,varid,ierr)
+
+            else
+                if (ntime==0) then
+                    write(*,*) "DIAGFI Error: failed creating variable ",trim(nom)
+                    write(*,*) "it seems it already exists!"
+                    call abort_physic("writediagpem",trim(nom)//" already exists",1)
+                endif
+            endif
+
+            corner(1)=ntime
+            edges(1)=1
+
+!#ifdef NC_DOUBLE
+!            ierr = NF_PUT_VARA_DOUBLE (nid,varid,[corner(1)],[edges(1)],[dx0])
+!#else
+            ierr= NF_PUT_VARA_REAL(nid,varid,[corner(1)],[edges(1)],[dx0])
+!#endif
+            if (ierr.ne.NF_NOERR) then
+                write(*,*) "***** PUT_VAR matter in writediagpem"
+                write(*,*) "***** with dx0: ",trim(nom)
+                write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
+                call abort_physic("writediagpem","failed writing "//trim(nom),1)
+            endif
+        endif !of if (is_master)
+    endif ! of if (dim.eq.3) elseif(dim.eq.2)...
+endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
+
+if (is_master) ierr= NF_CLOSE(nid)
+
+END SUBROUTINE writediagpem
+
+END MODULE writediagpem_mod
+
