A first version of the code. Not tested yet !!

.gitignore

+# Mac OS X
+*.DS_Store
+.DS_Store?
+._*
+.Spotlight-V100
+.Trashes
+ehthumbs.db
+Thumbs.db
+
+# Xcode
+*.pbxuser
+*.mode1v3
+*.mode2v3
+*.perspectivev3
+*.xcuserstate
+project.xcworkspace/
+xcuserdata/
+
+# Generated files
+*.[oa]
+*.pyc
+*.exe
+*.mod
+*.pdone
+*.bak
+
+# Backup files
+*~.nib
+*.*~
+*~
+
+# Specific files
+*.o*
+*.png
+Weights/
+__pycache__/
+*.so
+*.log
\ No newline at end of file

DomainDecomp.py

+#
+#
+#
+import numpy as np
+import os, sys
+from netCDF4 import Dataset
+from mpi4py import MPI
+comm = MPI.COMM_WORLD
+#
+sys.path.append('../')
+#
+# Get the information from the configuration file.
+#
+import configparser
+config=configparser.ConfigParser({'SaveWeights':'true', "DiagLon":"0.0, 0.0", "DiagLat":"0.0, 0.0"})
+config.read("../run.def")
+#
+import getargs
+log_master, log_world = getargs.getLogger(__name__)
+INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
+INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
+#
+import ModelGrid as MG
+import Partition as PA
+import RPPtools as RPP
+#
+import matplotlib.pylab as plt
+import savefile
+#
+gg = MG.GlobalGrid()
+#
+szhalo=1
+nbc=comm.Get_size()
+rank=comm.Get_rank()
+#
+# Open write unit
+#
+if not os.path.exists('Output'):
+    os.makedirs('Output')
+wunit = open("Output/TestDomain_out_"+str(rank)+".txt", "w")
+#
+wunit.write("START :: i = 0 :"+str(gg.ni-1)+" >> j = 0 :"+str(gg.nj-1)+" SZ="+str(gg.ni*gg.nj)+" nbland="+str(gg.nbland)+'\n')
+wunit.write("======================================"+'\n')
+#
+part = PA.partition(gg.ni, gg.nj, gg.land, comm, nbc, szhalo, rank, wunit=wunit)
+
+wunit.write("Domain on proc : "+str(part.ni)+str(part.nj)+" Offset from global domain : "+str(part.istart)+str(part.jstart)+'\n')
+wunit.write("Domain with halo on proc : "+str(part.nih)+str(part.njh)+" Offset with halo : "+str(part.ihstart)+str(part.jhstart)+'\n')
+wunit.write(" i inner = "+str(part.istart)+":"+str(part.istart+part.ni-1)+" i w.halo = "+str(part.ihstart)+":"+str(part.ihstart+part.nih-1))
+wunit.write(" j inner = "+str(part.jstart)+":"+str(part.jstart+part.nj-1)+" j w.halo = "+str(part.jhstart)+":"+str(part.jhstart+part.njh-1)+'\n')
+wunit.write("=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+======================"+'\n')
+
+x=np.zeros((part.njh,part.nih))
+ihoff=part.istart-part.ihstart
+jhoff=part.jstart-part.jhstart
+x[jhoff:jhoff+part.nj,ihoff:ihoff+part.ni]=99
+
+wunit.write("====== TEST FIELD START ======================"+'\n')
+for j in range(part.njh) :
+    wunit.write(str(x[j,:].astype(int))+'\n')
+
+part.sendtohalo(x)
+
+wunit.write("====== TEST FIELD AFTER HALO EXCHANGE ======================"+'\n')
+for j in range(part.njh) :
+    wunit.write(str(x[j,:].astype(int))+'\n')
+
+x=np.zeros((part.njh,part.nih))
+ihoff=part.istart-part.ihstart
+jhoff=part.jstart-part.jhstart
+x[jhoff:jhoff+part.nj,ihoff:ihoff+part.ni]=rank
+
+wunit.write("====== GATHER FIELD START ======================"+'\n')
+for j in range(part.njh) :
+    wunit.write(str(x[j,:].astype(int))+'\n')
+
+xgather = part.gather(x)
+
+if rank == 0 :
+    print("Type of gathered variable : ", xgather.astype(int))
+
+modelgrid=MG.ModelGrid(part.ihstart+gg.igstart, part.nih, part.jhstart+gg.jgstart, part.njh)
+INFO("Longitude interval on proc "+str(rank)+" = "+str(modelgrid.box_land[0][0])+" : "+str(modelgrid.box_land[0][1]))
+INFO("Latitude interval on proc "+str(rank)+" = "+str(modelgrid.box_land[1][0])+" : "+str(modelgrid.box_land[1][1]))
+
+rland=np.zeros((modelgrid.nbland))
+rland[:] = rank
+
+xland = part.landscatgat(modelgrid, rland)
+
+if rank == 0 :
+    xland[~np.isnan(xland)] = 1
+    xland[np.isnan(xland)] = 0
+    print("Difference between global and and the gathered version :",np.sum(xland-gg.land))
+
+#
+# Test the landsendtohalo function
+#
+x=np.zeros((part.njh,part.nih))
+x[:,:]=-1
+ihoff=part.istart-part.ihstart
+jhoff=part.jstart-part.jhstart
+x[jhoff:jhoff+part.nj,ihoff:ihoff+part.ni]=rank
+
+lx = modelgrid.landgather(x)
+part.landsendtohalo(lx)
+xland = part.landscatgat(modelgrid, lx)
+ 
+if rank == 0 :
+    xland[~np.isnan(xland)] = 1
+    xland[np.isnan(xland)] = 0
+    print("Difference between global and gathered version of land only :",np.sum(xland-gg.land))
+
+print("len areas : ", len(modelgrid.area))
+areas = part.landscatgat(modelgrid, np.array(modelgrid.area))
+if rank == 0 :
+    print(areas/1000000000.)
+wunit.close()
+
+y=np.zeros((nbc,modelgrid.nbland))
+print("Test : ", y.shape, modelgrid.nbland)
+y[rank,:]=rank
+
+savefile.dumpnetcdf("TestDomainDecomp.nc", gg, modelgrid, part, rland, x, y)

Environment

+#!/bin/bash
+#
+#########################################
+#
+# Script to build and maintain en environment to run
+# the python code in this directory. We verify that
+# all needed packages are available.
+# This is based on the anaconda python system : https://www.anaconda.com/distribution/
+#
+# Usage : source Environment
+#
+#########################################
+#
+# Functions
+#
+function hasanaconda (){
+    condahost="camelot climserv ciclad"
+    h=$(hostname)
+    res=1
+    for ch in ${condahost} ; do
+	if [ $(echo $h | grep ${ch} | wc -c) != "0" ] ; then
+	    res=0
+	fi
+    done
+    return $res
+}
+function checkenv (){
+    rm -f checkenv.py
+    cat <<EOF > checkenv.py
+import mpi4py
+import spherical_geometry
+import netCDF4
+import matplotlib
+import mpl_toolkits
+import configparser
+import numba
+print ("OK")
+EOF
+    python checkenv.py > /dev/null
+    if [ $? -gt 0 ] ; then
+	echo "Environment not OK. Some packages must be missing."
+	echo "This was tested with the checkenv.py code."
+    else
+	rm -f checkenv.py
+    fi
+}
+#
+# Main
+#
+if ( hasanaconda ) ; then
+    #
+    # Make sure the right version of Python
+    #
+    module unload python
+    module load python/3.6-anaconda50
+    #
+    # Test if the MPI environment is installed.
+    #
+    if [ ! -e ${HOME}/.conda/envs/MPI ] ; then
+	#
+	# Create the MPI environment as it does not exist.
+	#
+	conda create -y -n MPI mpi4py numba astropy netcdf4 matplotlib basemap ConfigParser
+	source activate MPI
+	conda install -y --channel astropy spherical-geometry
+    fi
+    #
+    # Activate the MPI environment
+    #
+    source activate MPI
+fi
+#
+# Verify we can load all the needed modules.
+#
+checkenv

F90subroutines/Makefile

+#
+#
+obj = defprec.o ioipsl.o constantes_var.o haversine.o grid.o routing_reg.o routing_tools.o
+#
+FC = gfortran
+FFLAG = -fPIC
+FDBG = -g -fbounds-check
+#
+FPY = f2py
+PDBG =  --debug
+#
+all : routing_interface.so
+
+clean :
+	rm -f compilation.log *.o *.so *.mod routing_interface.so *~ 
+
+routing_interface.so : $(obj) routing_interface.f90
+	rm -f compilation.log
+	$(FPY) -c routing_interface.f90 $(PDBG) --fcompiler=gfortran -m routing_interface $(obj) > compilation.log
+
+defprec.o : defprec.f90
+	$(FC) $(FFLAG) $(FDBG) -c defprec.f90
+
+ioipsl.o : ioipsl.f90
+	$(FC) $(FFLAG) $(FDBG) -c ioipsl.f90
+
+constantes_var.o : constantes_var.f90
+	$(FC) $(FFLAG) $(FDBG) -c constantes_var.f90
+
+haversine.o : haversine.f90
+	$(FC) $(FFLAG) $(FDBG) -c haversine.f90
+
+grid.o : defprec.f90 haversine.o constantes_var.o grid.f90
+	$(FC) $(FFLAG) $(FDBG) -c grid.f90
+
+routing_tools.o : defprec.o grid.o constantes_var.o routing_tools.f90
+	$(FC) $(FFLAG) $(FDBG) -c routing_tools.f90
+
+routing_reg.o : defprec.f90 grid.o constantes_var.o routing_tools.o routing_reg.f90
+	$(FC) $(FFLAG) $(FDBG) -c routing_reg.f90

F90subroutines/constantes_var.f90

+MODULE constantes_var
+  !
+  USE defprec
+  !
+  IMPLICIT NONE
+  !
+  REAL(r_std), PARAMETER            :: pi = 3.141592653589793238
+  REAL(r_std), PARAMETER            :: R_Earth = 6378000.
+  REAL(r_std), PARAMETER            :: mincos  = 0.0001   !! Minimum cosine value used for interpolation (unitless)
+  !
+  REAL(r_std), PARAMETER            :: zero = 0.0
+  REAL(r_std), PARAMETER            :: un = 1.0
+  INTEGER(i_std), PARAMETER         :: undef_int = 999999999
+  REAL(r_std), PARAMETER            :: undef_sechiba = 1.E+20
+  REAL(r_std), PARAMETER            :: val_exp = 999999.
+  !
+  INTEGER(i_std), PARAMETER         :: mpi_rank = 1
+  !
+END MODULE constantes_var

F90subroutines/defprec.f90

+MODULE defprec 
+!-
+! $Id: def.prec 386 2008-09-04 08:38:48Z bellier $
+!-
+! This software is governed by the CeCILL license
+! See IOIPSL/IOIPSL_License_CeCILL.txt
+!!--------------------------------------------------------------------
+!! The module "defprec" set default precision for computation
+!!
+!! This module should be used by every modules
+!! to keep the right precision for every variable
+!!--------------------------------------------------------------------
+!?INTEGERS of KIND 1 are not supported on all computers
+!?INTEGER,PARAMETER :: i_1=SELECTED_INT_KIND(2)
+  INTEGER,PARAMETER :: i_2=SELECTED_INT_KIND(4)
+  INTEGER,PARAMETER :: i_4=SELECTED_INT_KIND(9)
+  INTEGER,PARAMETER :: i_8=SELECTED_INT_KIND(13)
+  INTEGER,PARAMETER :: r_4=SELECTED_REAL_KIND(6,37)
+  INTEGER,PARAMETER :: r_8=SELECTED_REAL_KIND(15,307)
+  INTEGER,PARAMETER :: i_std=i_4, r_std=r_4
+!-----------------
+END MODULE defprec

F90subroutines/grid.f90

+MODULE grid
+  !
+  USE defprec
+  USE haversine
+  USE constantes_var
+  !
+  IMPLICIT NONE
+  !
+  INTEGER(i_std), SAVE                                     :: NbSegments  !! Number of segments in 
+  !                                                                  the polygone defining the grid box
+  INTEGER(i_std), SAVE                                     :: NbNeighb    !! Number of neighbours
+  !
+  REAL(r_std), ALLOCATABLE, DIMENSION (:), SAVE       :: contfrac
+  INTEGER(i_std), ALLOCATABLE, DIMENSION (:,:), SAVE  :: neighbours
+  REAL(r_std), ALLOCATABLE, DIMENSION(:,:), SAVE      :: headings
+  REAL(r_std), ALLOCATABLE, DIMENSION(:,:), SAVE      :: seglength
+  REAL(r_std), ALLOCATABLE, DIMENSION(:), SAVE        :: area
+  !
+  REAL(r_std), ALLOCATABLE, DIMENSION (:,:), SAVE     :: lalo
+  !
+  INTEGER(i_std), ALLOCATABLE, DIMENSION (:,:), SAVE  :: neighbours_g
+  REAL(r_std), ALLOCATABLE, DIMENSION (:,:), SAVE     :: lalo_g
+  REAL(r_std), ALLOCATABLE, DIMENSION(:), SAVE        :: area_g
+  REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:), SAVE    :: corners_g
+  REAL(r_std), ALLOCATABLE, DIMENSION(:,:), SAVE      :: headings_g
+  
+CONTAINS
+
+  SUBROUTINE grid_init(npts, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_in)
+
+    !! INPUT VARIABLES
+    INTEGER(i_std), INTENT(in) :: npts   !! Atmospheric Domain size (unitless)
+    INTEGER(i_std), INTENT(in) :: nbseg  !! Number of segments in the polygone
+    REAL(r_std), INTENT(in)    :: polygons_in(npts,2*nbseg+1,2)
+    REAL(r_std), INTENT(in)    :: centre_in(npts,2)
+    REAL(r_std), INTENT(in)    :: contfrac_in(npts)
+    REAL(r_std), INTENT(in)    :: area_in(npts)
+    INTEGER, INTENT(in) :: neighbours_in(npts,nbseg*2)
+    !
+    !
+    !! Local variables
+    INTEGER(i_std) :: i, is, iss
+    REAL(r_std), DIMENSION(npts,2*nbseg) :: lonpoly, latpoly
+    LOGICAL :: debug = .FALSE.
+    !
+    ! Define dimension of polygones
+    !
+    NbSegments = nbseg
+    NbNeighb = 2*NbSegments
+    !
+    ! Allocate memory and attribute values
+    !
+    IF ( .NOT. ALLOCATED(contfrac) ) THEN
+       ALLOCATE(contfrac(npts))
+       contfrac(:) = contfrac_in(:)
+    ENDIF
+    IF ( (.NOT.ALLOCATED(neighbours))) THEN
+       ALLOCATE(neighbours(npts,NbNeighb))
+       neighbours(:,:) = neighbours_in
+    ENDIF
+    IF ( (.NOT.ALLOCATED(neighbours_g))) THEN
+       ALLOCATE(neighbours_g(npts,NbNeighb))
+       neighbours_g(:,:) = neighbours_in
+    ENDIF
+    IF ( (.NOT.ALLOCATED(headings))) THEN
+       ALLOCATE(headings(npts,NbNeighb))
+       headings(:,:) = val_exp
+    ENDIF
+    IF ( (.NOT.ALLOCATED(headings_g))) THEN
+       ALLOCATE(headings_g(npts,NbNeighb))
+       headings_g(:,:) = val_exp
+    ENDIF
+    IF ( (.NOT.ALLOCATED(seglength))) THEN
+       ALLOCATE(seglength(npts,NbSegments))
+       seglength(:,:) = val_exp
+    ENDIF
+    IF ( (.NOT.ALLOCATED(corners_g))) THEN
+       ALLOCATE(corners_g(npts,NbSegments,2))
+       corners_g(:,:,:) = val_exp
+    ENDIF
+    IF ( (.NOT.ALLOCATED(area))) THEN
+       ALLOCATE(area(npts))
+       area(:) = area_in(:)
+    ENDIF
+    IF ( (.NOT.ALLOCATED(area_g))) THEN
+       ALLOCATE(area_g(npts))
+       area_g(:) = area_in(:)
+    ENDIF
+    IF ( (.NOT.ALLOCATED(lalo))) THEN
+       ALLOCATE(lalo(npts,2))
+       lalo(:,:) = val_exp
+    ENDIF
+    IF ( (.NOT.ALLOCATED(lalo_g))) THEN
+       ALLOCATE(lalo_g(npts,2))
+       lalo_g(:,:) = val_exp
+    ENDIF
+    !
+    ! Compute other variables we might need
+    !
+    !
+    ! Save the longitude and latitudes of the nbseg corners_g (=vertices)
+    !
+    DO i=1,npts
+       lonpoly(i,:) = polygons_in(i,1:2*nbseg,1)
+       latpoly(i,:) = polygons_in(i,1:2*nbseg,2)
+       !
+       lalo(i,1) = centre_in(i,2)
+       lalo(i,2) = centre_in(i,1)
+       lalo_g(i,1) = centre_in(i,2)
+       lalo_g(i,2) = centre_in(i,1)
+    ENDDO
+    !
+    ! Get the heading normal to the 4 segments and through the 4 corners_g.
+    !
+    CALL haversine_polyheadings(npts, nbseg, lonpoly, latpoly, centre_in, headings)
+    !
+    ! Order the points of the polygone in clockwise order Starting with the northern most
+    !
+    CALL haversine_polysort(npts, nbseg, lonpoly, latpoly, headings, neighbours)
+    !
+    DO i=1,npts
+       DO is=1,nbseg
+          iss=(is-1)*2+1
+          corners_g(i,is,1) = lonpoly(i,iss)
+          corners_g(i,is,2) = latpoly(i,iss)
+       ENDDO
+       headings_g(i,:) =  headings(i,:)
+    ENDDO
+    !
+    IF ( debug ) THEN
+       DO is=1,npts
+          CALL grid_diag(is, npts, centre_in)
+       ENDDO
+    ENDIF
+    !
+  END SUBROUTINE grid_init
+
+  SUBROUTINE grid_diag(is, npts, center)
+    !
+    !
+    INTEGER(i_std), INTENT(in) :: is, npts
+    REAL(r_std), INTENT(in)    :: center(npts,2)
+    !
+    CHARACTER(LEN=25) :: NW, NN, NE, EE, SE, SS, SW, WW, blk
+    !
+    blk = "              --       "
+    IF (neighbours(is,8) > 0) THEN
+       WRITE(NW,'(2F12.3)') center(neighbours(is,8),:)
+    ELSE
+       NW=blk
+    ENDIF
+    IF (neighbours(is,1) > 0) THEN
+       WRITE(NN,'(2F12.3)') center(neighbours(is,1),:)
+    ELSE
+       NN=blk
+    ENDIF
+    IF (neighbours(is,2) > 0) THEN
+       WRITE(NE,'(2F12.3)') center(neighbours(is,2),:)
+    ELSE
+       NE=blk
+    ENDIF
+    IF (neighbours(is,3) > 0) THEN
+       WRITE(EE,'(2F12.3)') center(neighbours(is,3),:)
+    ELSE
+       EE=blk
+    ENDIF    
+    IF (neighbours(is,4) > 0) THEN
+       WRITE(SE,'(2F12.3)') center(neighbours(is,4),:)
+    ELSE
+       SE=blk
+    ENDIF
+    IF (neighbours(is,5) > 0) THEN
+       WRITE(SS,'(2F12.3)') center(neighbours(is,5),:)
+    ELSE
+       SS=blk
+    ENDIF
+    IF (neighbours(is,6) > 0) THEN
+       WRITE(SW,'(2F12.3)') center(neighbours(is,6),:)
+    ELSE
+       SW=blk
+    ENDIF
+    IF (neighbours(is,7) > 0) THEN
+       WRITE(WW,'(2F12.3)') center(neighbours(is,7),:)
+    ELSE
+       WW=blk
+    ENDIF
+    WRITE(*,*) "AROUND : ", center(is,:)
+    WRITE(*,*) "---------------------------------------------------------------"
+    WRITE(*,*) NW," <--> ",       NN, " <--> ",            NE
+    WRITE(*,*) WW,"         XXXXXXXXXXXXX               ", EE
+    WRITE(*,*) SW," <--> ",       SS, " <--> ",            SE
+    WRITE(*,*) "---------------------------------------------------------------"
+    WRITE(*,*) headings(is,8)," <--> ",headings(is,1)," <--> ",  headings(is,2)
+    WRITE(*,*) headings(is,7)," <--> 000000000000000000 <--> ",  headings(is,3)
+    WRITE(*,*) headings(is,6)," <--> ",headings(is,5)," <--> ",  headings(is,4)
+    WRITE(*,*) "---------------------------------------------------------------"
+    
+  END SUBROUTINE grid_diag
+  
+END MODULE grid

F90subroutines/ioipsl.f90

+MODULE ioipsl
+  
+  IMPLICIT NONE
+
+  PUBLIC ::  ipslerr, ipslerr_p
+  
+  INTEGER, PARAMETER :: numout=6
+  INTEGER :: ipslout=6, ilv_cur=0, ilv_max=0
+  LOGICAL :: ioipsl_debug=.FALSE., lact_mode=.TRUE.
+  
+CONTAINS
+  !
+  SUBROUTINE ipslerr (plev,pcname,pstr1,pstr2,pstr3)
+!---------------------------------------------------------------------
+!! The "ipslerr" routine
+!! allows to handle the messages to the user.
+!!
+!! INPUT
+!!
+!! plev   : Category of message to be reported to the user
+!!          1 = Note to the user
+!!          2 = Warning to the user
+!!          3 = Fatal error
+!! pcname : Name of subroutine which has called ipslerr
+!! pstr1   
+!! pstr2  : Strings containing the explanations to the user
+!! pstr3
+!---------------------------------------------------------------------
+    IMPLICIT NONE
+!-
+    INTEGER :: plev
+    CHARACTER(LEN=*) :: pcname,pstr1,pstr2,pstr3
+!-
+    CHARACTER(LEN=30),DIMENSION(3) :: pemsg = &
+         &  (/ "NOTE TO THE USER FROM ROUTINE ", &
+         &     "WARNING FROM ROUTINE          ", &
+         &     "FATAL ERROR FROM ROUTINE      " /)
+!---------------------------------------------------------------------
+    IF ( (plev >= 1).AND.(plev <= 3) ) THEN
+       ilv_cur = plev
+       ilv_max = MAX(ilv_max,plev)
+       WRITE(ipslout,'(/,A," ",A)') TRIM(pemsg(plev)),TRIM(pcname)
+       WRITE(ipslout,'(3(" --> ",A,/))') TRIM(pstr1),TRIM(pstr2),TRIM(pstr3)
+    ENDIF
+    IF ( (plev == 3).AND.lact_mode) THEN
+       WRITE(ipslout,'("Fatal error from IOIPSL. STOP in ipslerr with code")')
+       STOP 1
+    ENDIF
+!---------------------
+  END SUBROUTINE ipslerr
+!
+  SUBROUTINE ipslerr_p (plev,pcname,pstr1,pstr2,pstr3)
+!---------------------------------------------------------------------
+!! The "ipslerr" routine
+!! allows to handle the messages to the user.
+!!
+!! INPUT
+!!
+!! plev   : Category of message to be reported to the user
+!!          1 = Note to the user
+!!          2 = Warning to the user
+!!          3 = Fatal error
+!! pcname : Name of subroutine which has called ipslerr
+!! pstr1