Corrected some stuff and ensured that the diagnostics from the F90 part are...

Corrected some stuff and ensured that the diagnostics from the F90 part are now written into files. This allows for finer diagnostics of the computation there.

.gitignore

@@ -39,6 +39,7 @@ xcuserdata/
 Weights/
 __pycache__/
 run.def
+Out*.txt
 #
 # Test configurations and directories
 #

Diags/plotgraphs.py

@@ -6,9 +6,9 @@ import numpy as np
 from netCDF4 import Dataset
 import sys
 
-NCfile="../MEDCORDEX_test_graph.nc"
+NCfile="../SPortugal_test_graph_n4.nc"
 GraphicFile = "test.png"
-title="Test"
+title="From file : "+NCfile
 #
 #####################################################################
 #
@@ -101,8 +101,8 @@ for il in range(landindex.nbland) :
     j,i = landindex.land2ij(il)
     xt,yt=m(lon_full[j,i],lat_full[j,i])
     plt.text(xt,yt,str(il))
-    #
-    for ih in range(nhtu) :
+    act_nhtu=np.count_nonzero(routetohtu[:,j,i] > 0)
+    for ih in range(act_nhtu) :
         lo = cglon[ih,j,i]
         la = cglat[ih,j,i]
         xl,yl=m(lo, la)
@@ -117,7 +117,7 @@ for il in range(landindex.nbland) :
         else :
             print("Unforseen coding of outlof HTU : ",int(routetohtu[ih,j,i]-1))
             sys.exit()
-        #
+        
         if int(routetohtu[ih,j,i]-1) < nhtu :
             ile =  int(routetogrid[ih,j,i]-1)
             je,ie = landindex.land2ij(ile)

DocumentationInterface

-initatmgrid(polygons_in,centre_in,area_in,contfrac_in,neighbours_in,[nbpt,nbseg])
+initatmgrid(rank_bn,nbcore,polygons_in,centre_in,area_in,contfrac_in,neighbours_in,[nbpt,nbseg])
 
 Wrapper for ``initatmgrid``.
 
 Parameters
 ----------
+rank_bn : input int
+nbcore : input int
 polygons_in : input rank-3 array('f') with bounds (nbpt,2 * nbseg + 1,2)
 centre_in : input rank-2 array('f') with bounds (nbpt,2)
 area_in : input rank-1 array('f') with bounds (nbpt)

F90subroutines/ioipsl.f90

@@ -4,7 +4,7 @@ MODULE ioipsl
 
   PUBLIC ::  ipslerr, ipslerr_p
   
-  INTEGER, PARAMETER :: numout=6
+  INTEGER, SAVE :: numout=6
   INTEGER :: ipslout=6, ilv_cur=0, ilv_max=0
   LOGICAL :: ioipsl_debug=.FALSE., lact_mode=.TRUE.
   

F90subroutines/routing_interface.f90

 
-SUBROUTINE initatmgrid(nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_in)
+SUBROUTINE initatmgrid(rank, nbcore, nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_in)
   !
   USE grid
+  USE ioipsl
   !
   IMPLICIT NONE
   !
   !! INPUT VARIABLES
+  INTEGER, INTENT(in) :: rank   !! rank of processor
+  INTEGER, INTENT(in) :: nbcore
   INTEGER, INTENT(in) :: nbpt   !! Atmospheric Domain size (unitless)
   INTEGER, INTENT(in) :: nbseg  !! Number of segments in the polygone
   REAL, INTENT(in) :: polygons_in(nbpt,2*nbseg+1,2)
@@ -17,6 +20,7 @@ SUBROUTINE initatmgrid(nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in
   !! Local
   INTEGER :: i
   INTEGER :: neighbours_loc(nbpt,nbseg*2)
+  CHARACTER(LEN=33) :: outfname
   !
   DO i=1,nbpt
      ! Move from C to F indexing.
@@ -25,8 +29,20 @@ SUBROUTINE initatmgrid(nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in
   !
   CALL grid_init(nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_loc)
   !
+  numout=100+rank
+  WRITE(outfname,'(A15,I4.4,A6,I4.4,A4)') "Out_RoutingReg_",rank,"_from_",nbcore,".txt"
+  OPEN(unit=numout, file=outfname)
+  !
 END SUBROUTINE initatmgrid
 
+SUBROUTINE closeinterface()
+  !
+  USE ioipsl
+  !
+  CLOSE(unit=numout)
+  !
+END SUBROUTINE closeinterface
+
 SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area,  max_basins, min_topoind, &
      & lon_rel, lat_rel, trip, basins, topoindex, fac, hierarchy, nbi, nbj, area_bx, trip_bx, basin_bx, &
      & topoind_bx, fac_bx, hierarchy_bx, lon_bx, lat_bx, lshead_bx, nwbas)

Interface.py

@@ -52,11 +52,18 @@ if gendoc.lower() == "true" :
 #
 # initatmgrid : Initialises the grid.f90 module and passes the description of the atmospheric grid.
 #
-def initatmgrid(nbpt, modelgrid) :
+def initatmgrid(rank, nbcore, nbpt, modelgrid) :
     print("INITATMGRID corners", np.array(modelgrid.polyll).shape)
     print("INITATMGRID area", np.array(modelgrid.area).shape)
     print("INITATMGRID neighbours", np.array(modelgrid.neighbours).shape)
-    routing_interface.initatmgrid( modelgrid.polyll, modelgrid.coordll, modelgrid.area, modelgrid.contfrac, modelgrid.neighbours)
+    routing_interface.initatmgrid(rank, nbcore, modelgrid.polyll, modelgrid.coordll, modelgrid.area, modelgrid.contfrac, modelgrid.neighbours)
+    return
+#
+#
+#
+def closeinterface(comm) :
+    comm.Barrier()
+    routing_interface.closeinterface()
     return
 #
 #
@@ -267,7 +274,19 @@ class HydroGraph :
         #
         # The field route_togrid is with indices on the local grid. That needs to be converted to the global grid.
         #
+        itarget=-1
+        for il in range(procgrid.nbland) :
+            lo = procgrid.lon_full[procgrid.indP[il][0],procgrid.indP[il][1]]
+            la = procgrid.lat_full[procgrid.indP[il][0],procgrid.indP[il][1]]
+            d=np.sqrt((lo-3.13)**2+(la-39.70)**2)
+            if d < 0.05 :
+                itarget = il
+                
+        if itarget >+ 0 :
+            print part.rank, itarget, " Before route_togrid = ", self.route_togrid[itarget,:]
         grgrid = part.l2glandindex(self.route_togrid[:,:])
+        if itarget >+ 0 :
+            print part.rank, itarget, " After route_togrid = ", self.route_togrid[itarget,:]
         rgrid = procgrid.landscatter(grgrid, order='F')
         rgrid = rgrid.astype(vtyp, copy=False)
         rgrid[rgrid >= RPP.IntFillValue] = NCFillValue

RoutingPreProc.py

@@ -130,7 +130,7 @@ print("nbasmax : ", nbasmax)
 #
 hydrodata = HG.HydroData(hydrogrid.ncfile, hydrogrid.box, sub_index)
 
-IF.initatmgrid(nbpt, modelgrid)
+IF.initatmgrid(rank, nbcore, nbpt, modelgrid)
 
 
 
@@ -177,3 +177,5 @@ hsuper.fetch()
 hgraph = IF.HydroGraph(nbasmax, hsuper)
 
 hgraph.dumpnetcdf(OutGraphFile, gg, modelgrid, part)
+
+IF.closeinterface(comm)