Changeset 965 for trunk/LMDZ.GENERIC/libf

Timestamp:

May 22, 2013, 9:10:28 AM (12 years ago)

Author:

emillour

Message:

Common dynamics and generic/universal GCM:

• LMDZ.COMMON: minor bug fix on the computation of physics mesh area in gcm.F
• LMDZ.UNIVERSAL: missing clean initialization of tab_cntrl(:) array in phyredem.F90
• LMDZ.GENERIC: minor bug fix in hydrol.F90, only output runoff if it is used. Update output routines so that all outputs files (stats, diagfi.nc, diagsoil.nc, diagspecIR.nc and diagspecVI.nc) can be generated when running LMDZ.UNIVERSAL in MPI mode.

EM

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• hydrol.F90 (modified) (2 diffs)
• inistats.F (modified) (4 diffs)
• iniwrite_specIR.F (modified) (2 diffs)
• iniwrite_specVI.F (modified) (2 diffs)
• iniwritesoil.F90 (modified) (7 diffs)
• mkstat.F90 (modified) (3 diffs)
• physiq.F90 (modified) (1 diff)
• writediagfi.F (modified) (15 diffs)
• writediagsoil.F90 (modified) (12 diffs)
• writediagspecIR.F (modified) (15 diffs)
• writediagspecVI.F (modified) (13 diffs)
• wstats.F90 (modified) (11 diffs)

trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90

@@ -110 +110 @@
          
          
-         ALLOCATE(runoff(ngrid))
+         if (activerunoff) ALLOCATE(runoff(ngrid))
 
          ivap=igcm_h2o_vap
@@ -360 +360 @@
       enddo
 
-      call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)
+      if (activerunoff) then
+        call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)
+      endif
 
       return

trunk/LMDZ.GENERIC/libf/phystd/inistats.F

@@ -1 +1 @@
       subroutine inistats(ierr)
 
+#ifdef CPP_PARA
+      use mod_phys_lmdz_para, only : is_master
+#endif
       implicit none
 
@@ -11 +14 @@
 #include "netcdf.inc"
 
+#ifndef CPP_PARA
+      logical,parameter :: is_master=.true.
+#endif
       integer,intent(out) :: ierr
       integer :: nid
@@ -47 +53 @@
          pseudoalt(l)=-10.*log(presnivs(l)/preff)   
       enddo
+
+      if (is_master) then
+      ! only the master needs do this
 
       ierr = NF_CREATE("stats.nc",NF_CLOBBER,nid)
@@ -115 +124 @@
       ierr=NF_CLOSE(nid)
 
+      endif ! of if (is_master)
       end

trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F

@@ -1 @@
-      SUBROUTINE iniwrite_specIR(nid,idayref,phis)
+      SUBROUTINE iniwrite_specIR(nid,idayref)
 
       use radinc_h, only: L_NSPECTI
       use radcommon_h, only: WNOI,DWNI
-      use comsoil_h
+!      use comsoil_h
 
       implicit none
@@ -40 +40 @@
 c   ----------
 
-      integer nid        ! NetCDF file ID
-      INTEGER*4 idayref  ! date (initial date for this run)
-      REAL phis(ip1jmp1) ! surface geopotential
+      integer,intent(in) :: nid        ! NetCDF file ID
+      INTEGER*4,intent(in) :: idayref  ! date (initial date for this run)
 
 c   Local:

trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F

@@ -1 @@
-      SUBROUTINE iniwrite_specVI(nid,idayref,phis)
+      SUBROUTINE iniwrite_specVI(nid,idayref)
 
       use radinc_h, only: L_NSPECTV
@@ -42 +42 @@
       integer nid        ! NetCDF file ID
       INTEGER*4 idayref  ! date (initial date for this run)
-      REAL phis(ip1jmp1) ! surface geopotential
 
 c   Local:

trunk/LMDZ.GENERIC/libf/phystd/iniwritesoil.F90

@@ -1 +1 @@
 subroutine iniwritesoil(nid,ngrid)
 
-use comsoil_h
+use comsoil_h, only : inertiedat, mlayer
+#ifdef CPP_PARA
+use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
+use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
+#endif
 
 ! initialization routine for 'writediagoil'. Here we create/define
@@ -34 +38 @@
 integer :: i,j,l,ig0
 
+#ifdef CPP_PARA
+! Added to work in parallel mode
+real dx3_glop(klon_glo,nsoilmx)
+real dx3_glo(iim,jjp1,nsoilmx) ! to store a global 3D data set
+#else
+logical,parameter :: is_master=.true.
+logical,parameter :: is_mpi_root=.true.
+#endif
+
 ! 1. Define the dimensions
+if (is_master) then
 ! Switch to NetCDF define mode
 ierr=NF_REDEF(nid)
@@ -183 +197 @@
 ! Note no need to write time variable here; it is done in writediagsoil.
 
+endif ! of if (is_master) 
+
 ! 3. Other variables to be included
 
 ! 3.1 mesh area surrounding each horizontal point
+if (is_master) then
 ierr=NF_REDEF(nid) ! switch to NetCDF define mode
 
@@ -219 +236 @@
 endif
 
+endif ! of if (is_master)
+
+
 ! 3.2 Thermal inertia
+if (is_master) then
 ierr=NF_REDEF(nid) ! switch to NetCDF define mode
 
@@ -241 +262 @@
 ierr=NF_PUT_ATT_TEXT(nid,varid,"units",len_trim(text),text)
 
+endif !of if (is_master)
+
 ! Recast data along 'dynamics' grid
+#ifdef CPP_PARA
+  ! gather field on a "global" (without redundant longitude) array
+  call Gather(inertiedat,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
+    data3(1:iim,:,:)=dx3_glo(1:iim,:,:)
+    data3(iip1,:,:)=data3(1,:,:)
+  endif
+!$OMP END MASTER
+!$OMP BARRIER
+#else
 ! Note: inertiedat is known from comsoil.h
-
 do l=1,nsoilmx
   ! handle the poles
@@ -250 +285 @@
     data3(i,jjp1,l)=inertiedat(ngrid,l)
   enddo
-  !!! THIS WILL NOT WORK IN PARALLEL !!!!
   ! rest of the grid
   do j=2,jjm
@@ -260 +294 @@
   enddo
 enddo ! of do l=1,nsoilmx
-
-! Write data2 to file
-ierr=NF_ENDDEF(nid) ! switch out of NetCDF define mode
-! Write
-#ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(nid,varid,data3)
-#else
-ierr=NF_PUT_VAR_REAL(nid,varid,data3)
-#endif
-if (ierr.ne.NF_NOERR) then
+#endif
+
+! Write data3 to file
+if (is_master) then
+ ierr=NF_ENDDEF(nid) ! switch out of NetCDF define mode
+ ! Write
+#ifdef NC_DOUBLE
+ ierr=NF_PUT_VAR_DOUBLE(nid,varid,data3)
+#else
+ ierr=NF_PUT_VAR_REAL(nid,varid,data3)
+#endif
+ if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not write th_inertia variable"
-endif
+ endif
+endif ! of if (is_master)
 
 end subroutine iniwritesoil

trunk/LMDZ.GENERIC/libf/phystd/mkstat.F90

@@ -9 +9 @@
 !  Yann W. july 2003
 
+#ifdef CPP_PARA
+use mod_phys_lmdz_para, only : is_master
+#endif
 
 implicit none
@@ -32 +35 @@
 integer :: meanid,sdid
 !integer, dimension(4) :: dimout
+#ifndef CPP_PARA
+logical,parameter :: is_master=.true.
+#endif
 
 ! Incrementation of count for the last step, which is not done in wstats
 count(istime)=count(istime)+1
+
+if (is_master) then
+! only the master needs do this
 
 ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
@@ -161 +170 @@
 ierr= NF_CLOSE(nid)
 
+endif ! of if (is_master)
+
 end

trunk/LMDZ.GENERIC/libf/phystd/physiq.F90

@@ -1686 +1686 @@
         call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
 
+!     Subsurface temperatures
+        call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
+        call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
+
 !     Total energy balance diagnostics
         if(callrad.and.(.not.newtonian))then

trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F

@@ -41 +41 @@
 
       USE surfdat_h
- 
+#ifdef CPP_PARA
+      USE mod_phys_lmdz_para, only : is_parallel, is_mpi_root,
+     &                               is_master, gather
+      USE mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
+#endif
       implicit none
 
@@ -70 +74 @@
       real*4,save :: date
 
-      REAL phis(ip1jmp1)
+      REAL phis(ip1jmp1) ! surface geopotential on extended lonxlat grid
 
       integer irythme
@@ -99 +103 @@
       integer dimvert
 
+#ifdef CPP_PARA
+! Added to work in parallel mode
+      real dx3_glop(klon_glo,llm)
+      real dx3_glo(iim,jjp1,llm) ! to store a global 3D data set
+      real dx2_glop(klon_glo)
+      real dx2_glo(iim,jjp1)     ! to store a global 2D (surface) data set
+      real px2(ngrid)
+!      real dx1_glo(llm)          ! to store a 1D (column) data set
+!      real dx0_glo
+      real phisfi_glo(klon_glo) ! surface geopotential on global physics grid
+#else
+      logical,parameter :: is_parallel=.false.
+      logical,parameter :: is_master=.true.
+      logical,parameter :: is_mpi_root=.true.
+      real phisfi_glo(ngrid) ! surface geopotential on global physics grid
+#endif
 !***************************************************************
 !Sortie des variables au rythme voulu
@@ -170 +190 @@
 
 #ifdef CPP_PARA
+          ! Gather phisfi() geopotential on physics grid
+          call Gather(phisfi,phisfi_glo)
+#else
+         phisfi_glo(:)=phisfi(:)
+#endif
+
          !! parallel: we cannot use the usual writediagfi method
-         call iophys_ini
-#else
+!!         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)
@@ -193 +221 @@
 
          ! write "header" of file (longitudes, latitudes, geopotential, ...)
-         call gr_fi_dyn(1,ngrid,iip1,jjp1,phisfi,phis)
+         call gr_fi_dyn(1,size(phisfi_glo),iip1,jjp1,phisfi_glo,phis)
          call iniwrite(nid,day_ini,phis)
 
+         endif ! of if (is_master)
+
       else
-         ! Open the NetCDF file
-         ierr = NF_OPEN(fichnom,NF_WRITE,nid)
-#endif
+
+         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')
 
@@ -223 +258 @@
       if ( MOD(zitau+1,irythme) .eq.0.) then
 
-#ifdef CPP_PARA
-         !! parallel: we cannot use the usual writediagfi method
-         if (dim .eq. 2) then
-             dimvert = 1
-         else if (dim == 3) then
-             dimvert = llm
-         endif
-         call iophys_ecrit(nom,dimvert,titre,unite,px)
-#else
+!#ifdef CPP_PARA
+!         !! parallel: we cannot use the usual writediagfi method
+!         if (dim .eq. 2) then
+!             dimvert = 1
+!         else if (dim == 3) then
+!             dimvert = llm
+!         endif
+!         call iophys_ecrit(nom,dimvert,titre,unite,px)
+!#else
 
 ! Compute/write/extend 'Time' coordinate (date given in days)
@@ -238 +273 @@
 !--------------------------------------------------------
 
+        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)
@@ -261 +298 @@
         end if ! of if (nom.eq.firstnom)
 
+        endif ! of if (is_master)
+
 !Case of a 3D variable
 !---------------------
         if (dim.eq.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:iim,:,:)=dx3_glo(1:iim,:,:)
+            dx3(iip1,:,:)=dx3(1,:,:)
+          endif
+!$OMP END MASTER
+!$OMP BARRIER
+#else
 !         Passage variable physique -->  variable dynamique
 !         recast (copy) variable from physics grid to dynamics grid
@@ -280 +332 @@
              ENDDO
            ENDDO
-
+#endif
 !         Ecriture du champs
 
-!         write (*,*) 'In  writediagfi, on sauve:  ' , nom
-!         write (*,*) 'In  writediagfi. Estimated date = ' ,date
+          if (is_master) then
+           ! only the master writes to output
 ! name of the variable
            ierr= NF_INQ_VARID(nid,nom,varid)
@@ -314 +366 @@
 !           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
            ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
 !#endif
@@ -320 +376 @@
               write(*,*) "***** PUT_VAR problem in writediagfi"
               write(*,*) "***** with ",nom
-              write(*,*) 'ierr=', ierr
+              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
 c             call abort
            endif 
 
+          endif !of if (is_master)
+
 !Case of a 2D variable
 !---------------------
 
         else if (dim.eq.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:iim,:)=dx2_glo(1:iim,:)
+            dx2(iip1,:)=dx2(1,:)
+          endif
+!$OMP END MASTER
+!$OMP BARRIER
+#else
 
 !         Passage variable physique -->  physique dynamique
@@ -343 +416 @@
                 dx2(iip1,j)=dx2(1,j)
              ENDDO
-
+#endif
+
+          if (is_master) then
+           ! only the master writes to output
 !         write (*,*) 'In  writediagfi, on sauve:  ' , nom
 !         write (*,*) 'In  writediagfi. Estimated date = ' ,date
@@ -379 +455 @@
               write(*,*) "***** PUT_VAR matter in writediagfi"
               write(*,*) "***** with ",nom
-              write(*,*) 'ierr=', ierr
+              write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
 c             call abort
            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(*,*) "writediagfi error: dim=1 not implemented ",
+     &                 "in parallel mode"
+           stop
+         endif
 !         Passage variable physique -->  physique dynamique
 !         recast (copy) variable from physics grid to dynamics grid
@@ -430 +513 @@
 
         else if (dim.eq.0) then
+         if (is_parallel) then
+           write(*,*) "writediagfi error: dim=0 not implemented ",
+     &                 "in parallel mode"
+           stop
+         endif
+
            dx0 = px (1,1)
 
@@ -462 +551 @@
 
         endif ! of if (dim.eq.3) elseif(dim.eq.2)...
-#endif
 
       endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
 
-#ifndef CPP_PARA
-      ierr= NF_CLOSE(nid)
-#endif
+      if (is_master) then
+        ierr= NF_CLOSE(nid)
+      endif
 
       end
+

trunk/LMDZ.GENERIC/libf/phystd/writediagsoil.F90

@@ -10 +10 @@
 ! (yielding the sought time series of the variable)
 
-use comsoil_h
+#ifdef CPP_PARA
+use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
+use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
+#endif
 
 implicit none
@@ -31 +34 @@
 
 ! Local variables:
-real,dimension(iip1,jjp1,nsoilmx) :: data3 ! to store 3D data
+real,dimension(iip1,jjp1,nsoilmx) :: data3 ! to store 3D data on extended lonxlat grid
 ! Note iip1,jjp1 known from paramet.h; nsoilmx known from dimphys.h
-real,dimension(iip1,jjp1) :: data2 ! to store 2D data
+real,dimension(iip1,jjp1) :: data2 ! to store 2D data on extended lonxlat grid
 real :: data0 ! to store 0D data
 integer :: i,j,l ! for loops
@@ -53 +56 @@
 integer,dimension(4) :: edges,corners
 
+#ifdef CPP_PARA
+! Added to work in parallel mode
+real dx3_glop(klon_glo,nsoilmx)
+real dx3_glo(iim,jjp1,nsoilmx) ! to store a global 3D data set
+real dx2_glop(klon_glo)
+real dx2_glo(iim,jjp1)     ! to store a global 2D (surface) data set
+real px2(ngrid)
+#else
+logical,parameter :: is_master=.true.
+logical,parameter :: is_mpi_root=.true.
+#endif
+
 ! 1. Initialization step
 if (firstname.eq."1234567890") then
@@ -72 +87 @@
   
   ! Create output NetCDF file
-  ierr=NF_CREATE(filename,NF_CLOBBER,nid)
-  if (ierr.ne.NF_NOERR) then
+  if (is_master) then
+   ierr=NF_CREATE(filename,NF_CLOBBER,nid)
+   if (ierr.ne.NF_NOERR) then
     write(*,*)'writediagsoil: Error, failed creating file '//trim(filename)
     stop
-  endif
-  
+   endif
+  endif ! of if (is_master)
+
   ! Define dimensions and axis attributes
   call iniwritesoil(nid,ngrid)
@@ -86 +103 @@
 else
   ! If not an initialization call, simply open the NetCDF file
-  ierr=NF_OPEN(filename,NF_WRITE,nid)
+  if (is_master) then
+   ierr=NF_OPEN(filename,NF_WRITE,nid)
+  endif
 endif ! of if (firstname.eq."1234567890")
 
@@ -105 +124 @@
     ! Note: day_step is known from control.h
     
-    ! Get NetCDF ID for "time"
-    ierr=NF_INQ_VARID(nid,"time",varid)
-    ! Add the current value of date to the "time" array
+    if (is_master) then
+     ! Get NetCDF ID for "time"
+     ierr=NF_INQ_VARID(nid,"time",varid)
+     ! Add the current value of 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
+     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(*,*)"writediagsoil: Failed writing date to time variable"
       stop 
-    endif
+     endif
+    endif ! of if (is_master)
   endif ! of if (name.eq.firstname)
 
@@ -122 +143 @@
 if (dimpx.eq.3) then ! Case of a 3D variable
   ! A. Recast data along 'dynamics' grid
+#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
+    data3(1:iim,:,:)=dx3_glo(1:iim,:,:)
+    data3(iip1,:,:)=data3(1,:,:)
+  endif
+!$OMP END MASTER
+!$OMP BARRIER
+#else
   do l=1,nsoilmx
     ! handle the poles
@@ -137 +171 @@
     enddo
   enddo
+#endif
   
   ! B. Write (append) the variable to the NetCDF file
+  if (is_master) then
   ! B.1. Get the ID of the variable
   ierr=NF_INQ_VARID(nid,name,varid)
@@ -176 +212 @@
                " to file "//trim(filename)//" at time",date
   endif
+  endif ! of if (is_master)
 
 elseif (dimpx.eq.2) then ! Case of a 2D variable
+
   ! A. Recast data along 'dynamics' grid
+#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 dx3_glo() to dx3(:) and add redundant longitude
+    data2(1:iim,:)=dx2_glo(1:iim,:)
+    data2(iip1,:)=data2(1,:)
+  endif
+!$OMP END MASTER
+!$OMP BARRIER
+#else
   ! handle the poles
   do i=1,iip1
@@ -192 +244 @@
     data2(iip1,j)=data2(1,j) ! extra (modulo) longitude
   enddo
+#endif
 
   ! B. Write (append) the variable to the NetCDF file
+  if (is_master) then
   ! B.1. Get the ID of the variable
   ierr=NF_INQ_VARID(nid,name,varid)
@@ -228 +282 @@
                " to file "//trim(filename)//" at time",date
   endif
+  endif ! of if (is_master)
 
 elseif (dimpx.eq.0) then ! Case of a 0D variable
+#ifdef CPP_PARA
+  write(*,*) "writediagsoil: dimps==0 case not implemented in // mode!!"
+  stop
+#endif
   ! A. Copy data value
   data0=px(1,1)
@@ -267 +326 @@
 
 ! 4. Close the NetCDF file
-ierr=NF_CLOSE(nid)
+if (is_master) then
+  ierr=NF_CLOSE(nid)
+endif
 
 end subroutine writediagsoil

trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F

@@ -1 @@
-      subroutine writediagspecIR(ngrid,nom,titre,unite,dim,px)
+      subroutine writediagspecIR(ngrid,nom,titre,unite,dimpx,px)
 
 !  Ecriture de variables diagnostiques au choix dans la physique 
@@ -32 +32 @@
 !      unite : unite de la variable (chaine de caracteres)
 !      px : variable a sortir (real 0, 2, ou 3d)
-!      dim : dimension de px : 0, 2, ou 3 dimensions
+!      dimpx : dimension de px : 0, 2, ou 3 dimensions
 !
 !=================================================================
@@ -43 +43 @@
 ! Addition by RW (2010) to allow OLR to be saved in .nc format
       use radinc_h, only : L_NSPECTI
-      USE surfdat_h
+!      USE surfdat_h, only : phisfi
+#ifdef CPP_PARA
+      use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
+      use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
+#endif
 
       implicit none
@@ -60 +64 @@
 
 ! Arguments on input:
-      integer ngrid
-      character (len=*) :: nom,titre,unite
-      integer dim
-      real px(ngrid,L_NSPECTI)
+      integer,intent(in) :: ngrid
+      character (len=*),intent(in) :: nom,titre,unite
+      integer,intent(in) :: dimpx
+      real,intent(in) :: px(ngrid,L_NSPECTI)
 
 ! Local variables:
@@ -73 +77 @@
       real date
 
-      REAL phis(ip1jmp1)
+!      REAL phis(ip1jmp1)
 
       integer irythme
@@ -99 +103 @@
        real dx3(iip1,jjp1,L_NSPECTI) ! to store the data set
 
+#ifdef CPP_PARA
+! Added to work in parallel mode
+      real dx3_glop(klon_glo,L_NSPECTI)
+      real dx3_glo(iim,jjp1,L_NSPECTI) ! to store a global 3D data set
+#else
+      logical,parameter :: is_master=.true.
+      logical,parameter :: is_mpi_root=.true.
+#endif
 
 !***************************************************************
@@ -109 +121 @@
 
 !***************************************************************
-
-! The following test is here to enforce that writediagfi is not used with the
-! 1D version of the GCM
-!not anymore (JL12)
-      if (ngrid.eq.-1) return
-      
-c     nom=trim((nom))
-c     unite=trim((unite))
-c     titre=trim((titre))
 
 ! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file
@@ -130 +133 @@
          ! just to be sure, check that firstnom is large enough to hold nom
          if (len_trim(firstnom).lt.len_trim(nom)) then
-           write(*,*) "writediagfi: Error !!!"
+           write(*,*) "writediagspecIR: Error !!!"
            write(*,*) "   firstnom string not long enough!!"
            write(*,*) "   increase its size to at least ",len_trim(nom)
@@ -137 +140 @@
 
          ! Create the NetCDF file
+         if (is_master) then
          ierr = NF_CREATE(fichnom, NF_CLOBBER, nid)
          ! Define the 'Time' dimension
@@ -155 +159 @@
          ierr = NF_ENDDEF(nid)
 
-         ! write "header" of file (longitudes, latitudes, geopotential, ...)
-         call gr_fi_dyn(1,ngrid,iip1,jjp1,phisfi,phis)
-!         call iniwrite(nid,day_ini,phis)
-         call iniwrite_specIR(nid,day_ini,phis)
+!         call gr_fi_dyn(1,size(phisfi_glo),iip1,jjp1,phisfi_glo,phis)
+         ! write "header" of file (longitudes, latitudes, area, ...)
+         call iniwrite_specIR(nid,day_ini)
+         endif ! of if (is_master)
 
          zitau = -1 ! initialize zitau
       else
-         ! Open the NetCDF file
-         ierr = NF_OPEN(fichnom,NF_WRITE,nid)
+         if (is_master) then
+           ! Open the NetCDF file
+           ierr = NF_OPEN(fichnom,NF_WRITE,nid)
+         endif
       endif ! if (firstnom.eq.'1234567890')
 
@@ -191 +197 @@
            ! compute corresponding date (in days and fractions thereof)
            date= float (zitau +1)/float (day_step)
-           ! Get NetCDF ID of 'Time' variable
-           ierr= NF_INQ_VARID(nid,"Time",varid)
-
-           !print*,'in writediagfi_spec.F, time=',varid
-
-           ! Write (append) the new date to the 'Time' array
+
+           if (is_master) then
+             ! 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
+             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 writediagspec_nc"
               write(*,*) "***** with time"
               write(*,*) 'ierr=', ierr   
 c             call abort
-           endif
-
-           write(6,*)'WRITEDIAGSPEC: date= ', date
+             endif
+
+             write(6,*)'WRITEDIAGSPEC: date= ', date
+           endif ! of if (is_master)
         end if ! of if (nom.eq.firstnom)
 
@@ -216 +223 @@
 !Case of a 3D variable
 !---------------------
-        if (dim.eq.3) then
-
-!         recast (copy) variable from physics grid to dynamics grid
+        if (dimpx.eq.3) then
+
+!         A. Recast (copy) variable from physics grid to dynamics grid
+#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:iim,:,:)=dx3_glo(1:iim,:,:)
+            dx3(iip1,:,:)=dx3(1,:,:)
+          endif
+!$OMP END MASTER
+!$OMP BARRIER
+#else
            DO l=1,L_NSPECTI
              DO i=1,iip1
@@ -232 +252 @@
              ENDDO
            ENDDO
+#endif
+
+!         B. Write (append) the variable to the NetCDF file
+          if (is_master) then
 
 ! name of the variable
@@ -245 +269 @@
 
               write (*,*) "=========================="
-              write (*,*) "DIAGSPEC: creating variable ",nom
+              write (*,*) "DIAGSPECIR: creating variable ",nom
               call def_var(nid,nom,titre,unite,4,id,varid,ierr)
 
@@ -272 +296 @@
            endif 
 
-        endif ! of if (dim.eq.3) elseif(dim.eq.2)...
+          endif ! of if (is_master)
+
+        endif ! of if (dimpx.eq.3)
 
       endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
 
-      ierr= NF_CLOSE(nid)
+      ! Close the NetCDF file
+      if (is_master) then
+        ierr= NF_CLOSE(nid)
+      endif
 
       end

trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F

@@ -1 @@
-      subroutine writediagspecVI(ngrid,nom,titre,unite,dim,px)
+      subroutine writediagspecVI(ngrid,nom,titre,unite,dimpx,px)
 
 !  Ecriture de variables diagnostiques au choix dans la physique 
@@ -32 +32 @@
 !      unite : unite de la variable (chaine de caracteres)
 !      px : variable a sortir (real 0, 2, ou 3d)
-!      dim : dimension de px : 0, 2, ou 3 dimensions
+!      dimpx : dimension de px : 0, 2, ou 3 dimensions
 !
 !=================================================================
@@ -43 +43 @@
 ! Addition by RW (2010) to allow OSR to be saved in .nc format
       use radinc_h, only : L_NSPECTV
-      USE surfdat_h
+!      USE surfdat_h
+#ifdef CPP_PARA
+      use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
+      use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
+#endif
 
       implicit none
@@ -62 +66 @@
       integer ngrid
       character (len=*) :: nom,titre,unite
-      integer dim
+      integer dimpx
       real px(ngrid,L_NSPECTV)
 
@@ -73 +77 @@
       real date
 
-      REAL phis(ip1jmp1)
+!      REAL phis(ip1jmp1)
 
       integer irythme
@@ -99 +103 @@
        real dx3(iip1,jjp1,L_NSPECTV) ! to store the data set
 
+#ifdef CPP_PARA
+! Added to work in parallel mode
+      real dx3_glop(klon_glo,L_NSPECTV)
+      real dx3_glo(iim,jjp1,L_NSPECTV) ! to store a global 3D data set
+#else
+      logical,parameter :: is_master=.true.
+      logical,parameter :: is_mpi_root=.true.
+#endif
 
 !***************************************************************
@@ -108 +120 @@
 
 !***************************************************************
-
-! The following test is here to enforce that writediagfi is not used with the
-! 1D version of the GCM
-!not anymore (JL12)
-      if (ngrid.eq.-1) return
-      
-c     nom=trim((nom))
-c     unite=trim((unite))
-c     titre=trim((titre))
 
 ! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file
@@ -136 +139 @@
 
          ! Create the NetCDF file
+         if (is_master) then
          ierr = NF_CREATE(fichnom, NF_CLOBBER, nid)
          ! Define the 'Time' dimension
@@ -155 +159 @@
 
          ! write "header" of file (longitudes, latitudes, geopotential, ...)
-         call gr_fi_dyn(1,ngrid,iip1,jjp1,phisfi,phis)
+!         call gr_fi_dyn(1,ngrid,iip1,jjp1,phisfi,phis)
 !         call iniwrite(nid,day_ini,phis)
-         call iniwrite_specVI(nid,day_ini,phis)
+         call iniwrite_specVI(nid,day_ini)
+         endif ! of if (is_master)
 
          zitau = -1 ! initialize zitau
       else
-         ! Open the NetCDF file
-         ierr = NF_OPEN(fichnom,NF_WRITE,nid)
+         if (is_master) then
+           ! Open the NetCDF file
+           ierr = NF_OPEN(fichnom,NF_WRITE,nid)
+         endif
       endif ! if (firstnom.eq.'1234567890')
 
@@ -190 +197 @@
            ! compute corresponding date (in days and fractions thereof)
            date= float (zitau +1)/float (day_step)
-           ! Get NetCDF ID of 'Time' variable
-           ierr= NF_INQ_VARID(nid,"Time",varid)
-
-           !print*,'in writediagfi_spec.F, time=',varid
-
-           ! Write (append) the new date to the 'Time' array
+
+           if (is_master) then
+             ! 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
+             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 writediagspec_nc"
               write(*,*) "***** with time"
               write(*,*) 'ierr=', ierr   
 c             call abort
-           endif
-
-           write(6,*)'WRITEDIAGSPEC: date= ', date
+             endif
+
+             write(6,*)'WRITEDIAGSPEC: date= ', date
+           endif ! of if (is_master)
         end if ! of if (nom.eq.firstnom)
 
@@ -215 +223 @@
 !Case of a 3D variable
 !---------------------
-        if (dim.eq.3) then
-
-!         recast (copy) variable from physics grid to dynamics grid
+        if (dimpx.eq.3) then
+
+!         A. Recast (copy) variable from physics grid to dynamics grid
+#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:iim,:,:)=dx3_glo(1:iim,:,:)
+            dx3(iip1,:,:)=dx3(1,:,:)
+          endif
+!$OMP END MASTER
+!$OMP BARRIER
+#else
            DO l=1,L_NSPECTV
              DO i=1,iip1
@@ -231 +252 @@
              ENDDO
            ENDDO
+#endif
+
+!         B. Write (append) the variable to the NetCDF file
+          if (is_master) then
 
 ! name of the variable
@@ -271 +296 @@
            endif 
 
-        endif ! of if (dim.eq.3) elseif(dim.eq.2)...
+          endif ! of if (is_master)
+
+        endif ! of if (dimpx.eq.3)
 
       endif ! of if ( MOD(zitau+1,irythme) .eq.0.)
 
-      ierr= NF_CLOSE(nid)
+      ! Close the NetCDF file
+      if (is_master) then
+        ierr= NF_CLOSE(nid)
+      endif
 
       end

trunk/LMDZ.GENERIC/libf/phystd/wstats.F90

@@ -1 +1 @@
 subroutine wstats(ngrid,nom,titre,unite,dim,px)
+
+#ifdef CPP_PARA
+use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather, klon_mpi_begin
+use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
+#endif
 
 implicit none
@@ -5 +10 @@
 #include "dimensions.h"
 #include "dimphys.h"
+#include "comconst.h"
 #include "statto.h"
 #include "netcdf.inc"
 
 integer,intent(in) :: ngrid
-character (len=*) :: nom,titre,unite
+character (len=*),intent(in) :: nom,titre,unite
 integer,intent(in) :: dim
 integer,parameter :: iip1=iim+1
 integer,parameter :: jjp1=jjm+1
-real, dimension(ngrid,llm) :: px
+real,intent(in) :: px(ngrid,llm)
 real, dimension(iip1,jjp1,llm) :: mean3d,sd3d,dx3
 real, dimension(iip1,jjp1) :: mean2d,sd2d,dx2
@@ -20 +26 @@
 integer :: ierr,varid,nbdim,nid
 integer :: meanid,sdid
-integer, dimension(4)  :: id,start,size
+integer, dimension(4)  :: id,start,sizes
 logical, save :: firstcall=.TRUE.
 integer :: l,i,j,ig0
-integer,save :: index
+integer,save :: indx
 
 integer, save :: step=0
 
-
+! Added to work in parallel mode
+#ifdef CPP_PARA
+real px3_glop(klon_glo,llm) ! to store a 3D data set on global physics grid
+real px3_glo(iim,jjp1,llm) ! to store a global 3D data set on global lonxlat grid
+real px2_glop(klon_glo) ! to store a 2D data set on global physics grid
+real px2_glo(iim,jjp1) ! to store a 2D data set on global lonxlat grid
+real px2(ngrid)
+real px3(ngrid,llm)
+#else
+! When not running in parallel mode:
+real px3_glop(ngrid,llm) ! to store a 3D data set on global physics grid
+real px3_glo(iim,jjp1,llm) ! to store a global 3D data set on global lonxlat grid
+real px2_glop(ngrid) ! to store a 2D data set on global physics grid
+real px2_glo(iim,jjp1) ! to store a 2D data set on global lonxlat grid
+logical,parameter :: is_master=.true.
+logical,parameter :: is_mpi_root=.true.
+integer,parameter :: klon_mpi_begin=1
+#endif
+
+! 1. Initialization (creation of stats.nc file)
 if (firstcall) then
    firstcall=.false.
@@ -34 +59 @@
 endif
 
-if (firstvar==nom) then ! If we're back to the first variable
+if (firstvar==nom) then ! If we're back to the first variable, increment time counter
       step = step + 1
 endif
 
 if (mod(step,istats).ne.0) then
+  ! if its not time to write to file, exit
    RETURN
 endif
 
+! collect fields on a global physics grid
+#ifdef CPP_PARA
+ if (dim.eq.3) then
+  px3(1:ngrid,1:llm)=px(1:ngrid,1:llm)
+  ! Gather fieds on a "global" (without redundant longitude) array
+  call Gather(px3,px3_glop)
+!$OMP MASTER
+  if (is_mpi_root) then
+    call Grid1Dto2D_glo(px3_glop,px3_glo)
+    ! copy dx3_glo() to dx3(:) and add redundant longitude
+    dx3(1:iim,:,:)=px3_glo(1:iim,:,:)
+    dx3(iip1,:,:)=dx3(1,:,:)
+  endif
+!$OMP END MASTER
+!$OMP BARRIER
+ else ! dim.eq.2
+  ! Gather fieds on a "global" (without redundant longitude) array
+  px2(:)=px(:,1)
+  call Gather(px2,px2_glop)
+!$OMP MASTER
+          if (is_mpi_root) then
+            call Grid1Dto2D_glo(px2_glop,px2_glo)
+            ! copy px2_glo() to dx2(:) and add redundant longitude
+            dx2(1:iim,:)=px2_glo(1:iim,:)
+            dx2(iip1,:)=dx2(1,:)
+          endif
+!$OMP END MASTER
+!$OMP BARRIER
+ endif
+#else
+  if (dim.eq.3) then
+    px3_glop(:,1:llm)=px(:,1:llm)
+!  Passage variable physique -->  variable dynamique
+    DO l=1,llm
+      DO i=1,iim
+         px3_glo(i,1,l)=px(1,l)
+         px3_glo(i,jjp1,l)=px(ngrid,l)
+      ENDDO
+      DO j=2,jjm
+         ig0= 1+(j-2)*iim
+         DO i=1,iim
+            px3_glo(i,j,l)=px(ig0+i,l)
+         ENDDO
+      ENDDO
+    ENDDO
+  else ! dim.eq.2
+    px2_glop(:)=px(:,1)
+!    Passage variable physique -->  physique dynamique
+   DO i=1,iim
+     px2_glo(i,1)=px(1,1)
+     px2_glo(i,jjp1)=px(ngrid,1)
+   ENDDO
+   DO j=2,jjm
+     ig0= 1+(j-2)*iim
+     DO i=1,iim
+        px2_glo(i,j)=px(ig0+i,1)
+     ENDDO
+   ENDDO
+  endif
+#endif
+
+! 2. Write field to file
+
+if (is_master) then
+! only master needs do this
+
 ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
 
 namebis=trim(nom)
+
+! test: check if that variable already exists in file
 ierr= NF_INQ_VARID(nid,namebis,meanid)
 
 if (ierr.ne.NF_NOERR) then
-
+  ! variable not in file, create/define it
    if (firstvar==nom) then 
-      index=1
-      count=0
+      indx=1
+      count(:)=0
    endif
 
@@ -73 +167 @@
    namebis=trim(nom)
    call def_var(nid,namebis,titre,unite,nbdim,id,meanid,ierr)
-   call inivar(nid,meanid,ngrid,dim,index,px,ierr)
+   if (dim.eq.3) then
+     call inivar(nid,meanid,size(px3_glop,1),dim,indx,px3_glop,ierr)
+   else ! dim.eq.2
+     call inivar(nid,meanid,size(px2_glop,1),dim,indx,px2_glop,ierr)
+   endif
    namebis=trim(nom)//"_sd"
    call def_var(nid,namebis,trim(titre)//" total standard deviation over the season",unite,nbdim,id,sdid,ierr)
-   call inivar(nid,sdid,ngrid,dim,index,px,ierr)
+   if (dim.eq.3) then
+     call inivar(nid,sdid,size(px3_glop,1),dim,indx,px3_glop,ierr)
+   else ! dim.eq.2
+     call inivar(nid,sdid,size(px2_glop,1),dim,indx,px2_glop,ierr)
+   endif
 
    ierr= NF_CLOSE(nid)
@@ -82 +184 @@
 
 else
+   ! variable found in file
    namebis=trim(nom)//"_sd"
    ierr= NF_INQ_VARID(nid,namebis,sdid)
@@ -88 +191 @@
 
 if (firstvar==nom) then 
-   count(index)=count(int(index))+1
-   index=index+1
-   if (index>istime) then
-      index=1
-   endif
-endif
-
-if (count(index)==0) then
-   if (dim.eq.3) then
-      start=(/1,1,1,index/)
-      size=(/iip1,jjp1,llm,1/)
-      mean3d=0
-      sd3d=0
+   count(indx)=count(int(indx))+1
+   indx=indx+1
+   if (indx>istime) then
+      indx=1
+   endif
+endif
+
+if (count(indx)==0) then
+   ! very first time we write the variable to file
+   if (dim.eq.3) then
+      start=(/1,1,1,indx/)
+      sizes=(/iip1,jjp1,llm,1/)
+      mean3d(:,:,:)=0
+      sd3d(:,:,:)=0
    else if (dim.eq.2) then
-      start=(/1,1,index,0/)
-      size=(/iip1,jjp1,1,0/)
-      mean2d=0
-      sd2d=0
+      start=(/1,1,indx,0/)
+      sizes=(/iip1,jjp1,1,0/)
+      mean2d(:,:)=0
+      sd2d(:,:)=0
    endif
 else
-   if (dim.eq.3) then
-      start=(/1,1,1,index/)
-      size=(/iip1,jjp1,llm,1/)
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,size,mean3d)
-      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,size,sd3d)
-#else
-      ierr = NF_GET_VARA_REAL(nid,meanid,start,size,mean3d)
-      ierr = NF_GET_VARA_REAL(nid,sdid,start,size,sd3d)
+   ! load values from file
+   if (dim.eq.3) then
+      start=(/1,1,1,indx/)
+      sizes=(/iip1,jjp1,llm,1/)
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
+      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
+#else
+      ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean3d)
+      ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd3d)
 #endif
       if (ierr.ne.NF_NOERR) then
@@ -124 +229 @@
 
    else if (dim.eq.2) then
-      start=(/1,1,index,0/)
-      size=(/iip1,jjp1,1,0/)
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,size,mean2d)
-      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,size,sd2d)
-#else
-      ierr = NF_GET_VARA_REAL(nid,meanid,start,size,mean2d)
-      ierr = NF_GET_VARA_REAL(nid,sdid,start,size,sd2d)
+      start=(/1,1,indx,0/)
+      sizes=(/iip1,jjp1,1,0/)
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
+      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
+#else
+      ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean2d)
+      ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd2d)
 #endif
       if (ierr.ne.NF_NOERR) then
@@ -138 +243 @@
       endif
    endif
-endif
+endif ! of if (count(indx)==0)
+
+
+! 2.1. Build dx* (data on lon-lat grid, with redundant longitude)
 
 if (dim.eq.3) then
+  dx3(1:iim,:,:)=px3_glo(:,:,:)
+  dx3(iip1,:,:)=dx3(1,:,:)
+else ! dim.eq.2
+  dx2(1:iim,:)=px2_glo(:,:)
+  dx2(iip1,:)=dx2(1,:)
+endif
+
+
+! 2.2. Add current values to previously stored sums
+
+if (dim.eq.3) then
+
+   mean3d(:,:,:)=mean3d(:,:,:)+dx3(:,:,:)
+   sd3d(:,:,:)=sd3d(:,:,:)+dx3(:,:,:)**2
+
+#ifdef NC_DOUBLE
+   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
+   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
+#else
+   ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean3d)
+   ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd3d)
+#endif
+
+else if (dim.eq.2) then
+
+   mean2d(:,:)= mean2d(:,:)+dx2(:,:)
+   sd2d(:,:)=sd2d(:,:)+dx2(:,:)**2
+
+#ifdef NC_DOUBLE
+   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
+   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
+#else
+   ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean2d)
+   ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd2d)
+#endif
+
+endif ! of if (dim.eq.3) elseif (dim.eq.2)
+
+  ierr= NF_CLOSE(nid)
+endif ! of if (is_master)
+
+end
+
+!======================================================
+subroutine inivar(nid,varid,ngrid,dim,indx,px,ierr)
+
+implicit none
+
+include "dimensions.h"
+include "dimphys.h"
+include "netcdf.inc"
+
+integer, intent(in) :: nid,varid,dim,indx,ngrid
+real, dimension(ngrid,llm), intent(in) :: px
+integer, intent(out) :: ierr
+
+integer,parameter :: iip1=iim+1
+integer,parameter :: jjp1=jjm+1
+
+integer :: l,i,j,ig0
+integer, dimension(4) :: start,sizes
+real, dimension(iip1,jjp1,llm) :: dx3
+real, dimension(iip1,jjp1) :: dx2
+
+if (dim.eq.3) then
+
+   start=(/1,1,1,indx/)
+   sizes=(/iip1,jjp1,llm,1/)
 
 !  Passage variable physique -->  variable dynamique
@@ -158 +334 @@
    ENDDO
 
-   mean3d= mean3d+dx3
-   sd3d= sd3d+dx3**2
-
-#ifdef NC_DOUBLE
-   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,size,mean3d)
-   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,size,sd3d)
-#else
-   ierr = NF_PUT_VARA_REAL(nid,meanid,start,size,mean3d)
-   ierr = NF_PUT_VARA_REAL(nid,sdid,start,size,sd3d)
+#ifdef NC_DOUBLE
+   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx3)
+#else
+   ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx3)
 #endif
 
 else if (dim.eq.2) then
+
+      start=(/1,1,indx,0/)
+      sizes=(/iip1,jjp1,1,0/)
 
 !    Passage variable physique -->  physique dynamique
@@ -185 +359 @@
   ENDDO
 
-   mean2d= mean2d+dx2
-   sd2d= sd2d+dx2**2
-
-#ifdef NC_DOUBLE
-   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,size,mean2d)
-   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,size,sd2d)
-#else
-   ierr = NF_PUT_VARA_REAL(nid,meanid,start,size,mean2d)
-   ierr = NF_PUT_VARA_REAL(nid,sdid,start,size,sd2d)
-#endif
-
-endif
-
-ierr= NF_CLOSE(nid)
+#ifdef NC_DOUBLE
+   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx2)
+#else
+   ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx2)
+#endif
+
+endif
 
 end
-
-!======================================================
-subroutine inivar(nid,varid,ngrid,dim,index,px,ierr)
-
-implicit none
-
-include "dimensions.h"
-include "dimphys.h"
-include "netcdf.inc"
-
-integer, intent(in) :: nid,varid,dim,index,ngrid
-real, dimension(ngrid,llm), intent(in) :: px
-integer, intent(out) :: ierr
-
-integer,parameter :: iip1=iim+1
-integer,parameter :: jjp1=jjm+1
-
-integer :: l,i,j,ig0
-integer, dimension(4) :: start,size
-real, dimension(iip1,jjp1,llm) :: dx3
-real, dimension(iip1,jjp1) :: dx2
-
-if (dim.eq.3) then
-
-   start=(/1,1,1,index/)
-   size=(/iip1,jjp1,llm,1/)
-
-!  Passage variable physique -->  variable dynamique
-
-   DO l=1,llm
-      DO i=1,iip1
-         dx3(i,1,l)=px(1,l)
-         dx3(i,jjp1,l)=px(ngrid,l)
-      ENDDO
-      DO j=2,jjm
-         ig0= 1+(j-2)*iim
-         DO i=1,iim
-            dx3(i,j,l)=px(ig0+i,l)
-         ENDDO
-         dx3(iip1,j,l)=dx3(1,j,l)
-      ENDDO
-   ENDDO
-
-#ifdef NC_DOUBLE
-   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,dx3)
-#else
-   ierr = NF_PUT_VARA_REAL(nid,varid,start,size,dx3)
-#endif
-
-else if (dim.eq.2) then
-
-      start=(/1,1,index,0/)
-      size=(/iip1,jjp1,1,0/)
-
-!    Passage variable physique -->  physique dynamique
-
-  DO i=1,iip1
-     dx2(i,1)=px(1,1)
-     dx2(i,jjp1)=px(ngrid,1)
-  ENDDO
-  DO j=2,jjm
-     ig0= 1+(j-2)*iim
-     DO i=1,iim
-        dx2(i,j)=px(ig0+i,1)
-     ENDDO
-     dx2(iip1,j)=dx2(1,j)
-  ENDDO
-
-#ifdef NC_DOUBLE
-   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,dx2)
-#else
-   ierr = NF_PUT_VARA_REAL(nid,varid,start,size,dx2)
-#endif
-
-endif
-
-end