Changeset 4151

Timestamp:

11/22/18 11:45:17 (5 years ago)

Author:

omamce

Message:

O.M. : adding comments to cotes_etal

File:

• TOOLS/MOSAIX/cotes_etal.py (modified) (16 diffs)

TOOLS/MOSAIX/cotes_etal.py

@@ -1 +1 @@
 # -*- Mode: python -*-
-import netCDF4
-import numpy as np
-import nemo
-from scipy import ndimage
-import sys, os, platform, argparse, textwrap, time
-
-zero = np.dtype('float64').type(0.0)
-zone = np.dtype('float64').type(1.0)
-epsfrac = np.dtype('float64').type(1.0E-10)
-pi  = np.pi
-rad = pi/np.dtype('float64').type(180.0)
-ra = np.dtype('float64').type(6400000.0) # Earth radius
-
-def geodist (plon1, plat1, plon2, plat2) :
-      """Distance between two points (on sphere)"""
-      zs = np.sin (rad*plat1) * np.sin (rad*plat2) +  np.cos (rad*plat1) * np.cos (rad*plat2) * np.cos(rad*(plon2-plon1))
-      zs = np.maximum (-zone, np.minimum (zone, zs))
-      geodist =  np.arccos (zs)
-      return geodist
-
+### ===========================================================================
 ###
-###   A partir d'un fichier de poids, complete avec les rivieres et
-###   les poids pour le run off
-###   On met a zero les poids sur les points non cotiers, et on renormalise
+### Compute runoff weights.
+### For LMDZ only. Not suitable for DYNAMICO
 ###
-
-   
-### Lecture de la ligne de comande
+### ===========================================================================
+##
+##  Warning, to install, configure, run, use any of Olivier Marti's
+##  software or to read the associated documentation you'll need at least
+##  one (1) brain in a reasonably working order. Lack of this implement
+##  will void any warranties (either express or implied).
+##  O. Marti assumes no responsability for errors, omissions,
+##  data loss, or any other consequences caused directly or indirectly by
+##  the usage of his software by incorrectly or partially configured
+##  personal.
+##
 # SVN information
 __Author__   = "$Author$"
@@ -35 +23 @@
 __HeadURL__  = "$HeadURL$"
 __SVN_Date__ = "$SVN_Date: $"
-
-###
-
-### ===== Handling command line parameters ==================================================
+##
+
+import netCDF4
+import numpy as np
+import nemo
+from scipy import ndimage
+import sys, os, platform, argparse, textwrap, time
+
+## Userful constants
+zero = np.dtype('float64').type(0.0)
+zone = np.dtype('float64').type(1.0)
+epsfrac = np.dtype('float64').type(1.0E-10)
+pi  = np.pi
+rad = pi/np.dtype('float64').type(180.0)
+ra = np.dtype('float64').type(6400000.0) # Earth radius
+
+## Functions
+def geodist (plon1, plat1, plon2, plat2) :
+      """Distance between two points (on sphere)"""
+      zs = np.sin (rad*plat1) * np.sin (rad*plat2) +  np.cos (rad*plat1) * np.cos (rad*plat2) * np.cos(rad*(plon2-plon1))
+      zs = np.maximum (-zone, np.minimum (zone, zs))
+      geodist =  np.arccos (zs)
+      return geodist
+
+
+   
+
+### ===== Reading command line parameters ==================================================
 # Creating a parser
 parser = argparse.ArgumentParser (
     description = """Compute calving weights""",
-    epilog='-------- This is the end of the help message --------')
+    epilog='-------- End of the help message --------')
 
 # Adding arguments
@@ -48 +60 @@
 parser.add_argument ('--atm'   , help='atm model name (LMD*)', type=str, default='LMD9695'    )
 parser.add_argument ('--atmCoastWidth', help='width of the coastal band in atmosphere (in grid points)', default=1 )
-parser.add_argument ('--oceCoastWidth', help='width of the coastal band in ocean (in grid points)', default=1 )
+parser.add_argument ('--oceCoastWidth', help='width of the coastal band in ocean (in grid points)'     , default=2 )
 parser.add_argument ('--grids', help='grids file name', default='grids.nc' )
 parser.add_argument ('--areas', help='masks file name', default='areas.nc' )
@@ -55 +67 @@
 parser.add_argument ('--output', help='output rmp file name', default='rmp_tlmd_to_torc_runoff.nc' )
 
-
 # Parse command line
 myargs = parser.parse_args()
@@ -68 +79 @@
 atm_Name = myargs.atm
 oce_Name = myargs.oce
-# Largeur de la bande cotiere cote atm
+# Width of the coastal band (land points) in the atmopshere
 atmCoastWidth = myargs.atmCoastWidth
-# Largeur de la bande cotiere cote oce
+# Width of the coastal band (ocean points) in the ocean 
 oceCoastWidth = myargs.oceCoastWidth
 
@@ -94 +105 @@
 atm_grid_dims       = atm_grid_area.shape
 (atm_nvertex, atm_jpj, atm_jpi) = atm_grid_corner_lat.shape
-if atm_jpi == 1 : atm_perio = 0 # Icosahedre
-else            : atm_perio = 1 #
+atm_perio = 0
 atm_grid_pmask = nemo.lbc_mask (atm_grid_imask, 'T', atm_perio)
 atm_address = np.reshape ( np.arange(atm_jpj*atm_jpi), (atm_jpj, atm_jpi) )
@@ -115 +125 @@
 oce_grid_size = oce_jpj*oce_jpi
 
-print ('Calculs points terre/oce/cote sur grille atmosphere' )
-
-# Fill closed sea with image processing library
+## Fill closed sea with image processing library
 oce_grid_imask = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(oce_grid_imask, nperio=oce_perio, cd_type='T')), nperio=oce_perio, cd_type='T' )
 
+## 
 print ("Determination d'une bande cotiere ocean")
 
-oceLand  = np.where (oce_grid_imask[:] < 0.5, True, False)
-oceOcean = np.where (oce_grid_imask[:] > 0.5, True, False)
-
-oceCoast = np.where (ndimage.uniform_filter(oceOcean.astype(float), size=1+2*oceCoastWidth)<1.0,True,False) & oceOcean
+oceLand  = np.where (oce_grid_pmask[:] < 0.5, True, False)
+oceOcean = np.where (oce_grid_pmask[:] > 0.5, True, False)
+
+NNocean = 1+2*oceCoastWidth
+oceOceanFiltered = ndimage.uniform_filter(oceOcean.astype(float), size=NNocean)
+oceCoast = np.where (oceOceanFiltered<(1.0-0.5/(NNocean**2)),True,False) & oceOcean
 oceCoast = nemo.lbc_mask (oceCoast, oce_perio, 'T')
 
-print ('Nombre de point oceLand  : ' + str(oceLand.sum()) )
-print ('Nombre de point oceOcean : ' + str(oceOcean.sum()) )
-print ('Nombre de point oceCoast : ' + str(oceCoast.sum()) )
-
-# Faire une liste de ces points
+print ('Number of points in oceLand  : ' + str(oceLand.sum()) )
+print ('Number of points in oceOcean : ' + str(oceOcean.sum()) )
+print ('Number of points in oceCoast : ' + str(oceCoast.sum()) )
+
+# Arrays with coastal points only
 oceCoast_grid_center_lon = oce_grid_center_lon[oceCoast]
 oceCoast_grid_center_lat = oce_grid_center_lat[oceCoast]
@@ -142 +153 @@
 
 print ("Determination d'une bande cotiere atmosphere " )
-atmLand  = np.where (o2aFrac[:] < epsfrac, True, False)
-atmFrac  = np.where (o2aFrac[:] < epsfrac, True, False) & np.where (o2aFrac[:] < (zone-epsfrac), True, False)
-atmOcean = np.where (o2aFrac[:] < (zone-epsfrac), True, False)
-
-atmCoast = np.where (ndimage.uniform_filter(atmLand.astype(float), size=1+2*atmCoastWidth)<1.0,True,False) & atmLand
+atmLand      = np.where (o2aFrac[:] < epsfrac       , True, False)
+atmLandFrac  = np.where (o2aFrac[:] < zone-epsfrac  , True, False)
+atmFrac      = np.where (o2aFrac[:] > epsfrac       , True, False) & np.where (o2aFrac[:] < (zone-epsfrac), True, False)
+atmOcean     = np.where (o2aFrac[:] < (zone-epsfrac), True, False)
+atmOceanFrac = np.where (o2aFrac[:] > epsfrac       , True, False)
+
+NNatm = 1+2*atmCoastWidth
+atmLandFiltered = ndimage.uniform_filter(atmLand.astype(float), size=NNatm)
+atmCoast = np.where (atmLandFiltered<(1.0-0.5/(NNatm**2)),True,False) & atmLandFrac
 atmCoast = nemo.lbc_mask (atmCoast, 1, 'T')
 
-print ('Nombre de point atmLand  : ' + str(atmLand.sum()) )
-print ('Nombre de point atmOcean : ' + str(atmOcean.sum()) )
-print ('Nombre de point atmCoast : ' + str(atmCoast.sum()) )
-
-# Faire une liste de ces points
+print ('Number of points in atmLand  : ' + str(atmLand.sum()) )
+print ('Number of points in atmOcean : ' + str(atmOcean.sum()) )
+print ('Number of points in atmCoast : ' + str(atmCoast.sum()) )
+
+# Arrays with coastal points only
 atmCoast_grid_center_lon = atm_grid_center_lon[atmCoast]
 atmCoast_grid_center_lat = atm_grid_center_lat[atmCoast]
@@ -165 +180 @@
 oce_address  = np.empty ( shape=(0), dtype=np.int32   )
 
+## Loop on atmosphere coastal points
 for ja in np.arange(len(atmCoast_grid_pmask)) :
     z_dist = geodist ( atmCoast_grid_center_lon[ja], atmCoast_grid_center_lat[ja], oceCoast_grid_center_lon, oceCoast_grid_center_lat)
@@ -185 +201 @@
         
 
-print ('fini')
+print ('End of loop')
 
 num_links = remap_matrix.shape[0]
 
-# Output file
+### Output file
 runoff = myargs.output
 f_runoff = netCDF4.Dataset ( runoff, 'w', format='NETCDF3_64BIT' )
@@ -203 +219 @@
 f_runoff.associatedFiles = grids + " " + areas + " " + masks
 f_runoff.directory       = os.getcwd ()
-f_runoff.description     = "Generated with ???"
+f_runoff.description     = "Generated with cotes_etal.py"
 f_runoff.title           = runoff
 f_runoff.Program         = "Generated by " + sys.argv[0] + " with flags " + str(sys.argv[1:])
@@ -242 +258 @@
 
 v_remap_matrix[:] = remap_matrix
-v_atm_address [:] = atm_address
-v_oce_address [:] = oce_address
+v_atm_address [:] = atm_address + 1 # OASIS uses Fortran style indexing
+v_oce_address [:] = oce_address + 1
 
 v_atm_grid_dims       = f_runoff.createVariable ( 'src_grid_dims'      , 'i4', ('src_grid_rank',) )
@@ -256 +272 @@
 v_atm_grid_area       = f_runoff.createVariable ( 'src_grid_area'      , 'f8', ('src_grid_size',)  )
 v_atm_grid_area.long_name="Grid area" ; v_atm_grid_area.standard_name="cell_area" ; v_atm_grid_area.units="m2"
-v_atm_grid_imask      = f_runoff.createVariable ( 'src_grid_imask'     , 'f8', ('src_grid_size',)  )
+v_atm_grid_imask      = f_runoff.createVariable ( 'src_grid_imask'     , 'i4', ('src_grid_size',)  )
 v_atm_grid_imask.long_name="Land-sea mask" ; v_atm_grid_imask.units="Land:1, Ocean:0"
+v_atm_grid_pmask      = f_runoff.createVariable ( 'src_grid_pmask'     , 'i4', ('src_grid_size',)  )
+v_atm_grid_pmask.long_name="Land-sea mask (periodicity removed)" ; v_atm_grid_pmask.units="Land:1, Ocean:0"
 
 v_atm_grid_dims      [:] = atm_grid_dims
@@ -266 +284 @@
 v_atm_grid_area      [:] = atm_grid_area[:].ravel()
 v_atm_grid_imask     [:] = atm_grid_imask[:].ravel()
+v_atm_grid_pmask     [:] = atm_grid_pmask[:].ravel()
 
 # --
@@ -274 +293 @@
 v_oce_grid_center_lon.units='degrees_east'  ; v_oce_grid_center_lon.long_name='Longitude' ; v_oce_grid_center_lon.long_name='longitude' ; v_oce_grid_center_lon.bounds="dst_grid_corner_lon"
 v_oce_grid_center_lat.units='degrees_north' ; v_oce_grid_center_lat.long_name='Latitude'  ; v_oce_grid_center_lat.long_name='latitude'  ; v_oce_grid_center_lat.bounds="dst_grid_corner_lat"
-v_oce_grid_corner_lon = f_runoff.createVariable ( 'oce_grid_corner_lon', 'f8', ('dst_grid_size', 'dst_grid_corners')  )
+v_oce_grid_corner_lon = f_runoff.createVariable ( 'dst_grid_corner_lon', 'f8', ('dst_grid_size', 'dst_grid_corners')  )
 v_oce_grid_corner_lat = f_runoff.createVariable ( 'dst_grid_corner_lat', 'f8', ('dst_grid_size', 'dst_grid_corners')  )
 v_oce_grid_corner_lon.units="degrees_east"
 v_oce_grid_corner_lat.units="degrees_north"
-v_oce_grid_area       = f_runoff.createVariable ( 'dst_grid_area'      , 'f8', ('dst_grid_size',)  )
+v_oce_grid_area       = f_runoff.createVariable ( 'dst_grid_area'  , 'f8', ('dst_grid_size',) )
 v_oce_grid_area.long_name="Grid area" ; v_oce_grid_area.standard_name="cell_area" ; v_oce_grid_area.units="m2"
-v_oce_grid_imask      = f_runoff.createVariable ( 'dst_grid_imask'     , 'f8', ('dst_grid_size',)  )
+v_oce_grid_imask      = f_runoff.createVariable ( 'dst_grid_imask'     , 'i4', ('dst_grid_size',)  )
 v_oce_grid_imask.long_name="Land-sea mask" ; v_oce_grid_imask.units="Land:1, Ocean:0"
+v_oce_grid_pmask      = f_runoff.createVariable ( 'dst_grid_pmask'     , 'i4', ('dst_grid_size',)  )
+v_oce_grid_pmask.long_name="Land-sea mask (periodicity removed)" ; v_oce_grid_pmask.units="Land:1, Ocean:0"
 
 v_oce_grid_dims      [:] = oce_grid_dims
@@ -290 +311 @@
 v_oce_grid_area      [:] = oce_grid_area[:].ravel()
 v_oce_grid_imask     [:] = oce_grid_imask[:].ravel()
-
-# For diags
-v_atm_lon_addressed = f_runoff.createVariable ( 'src_lon_addressed', 'f8', ('num_links',) )
-v_atm_lat_addressed = f_runoff.createVariable ( 'src_lat_addressed', 'f8', ('num_links',) )
-v_oce_lon_addressed = f_runoff.createVariable ( 'dst_lon_addressed', 'f8', ('num_links',) )
-v_oce_lat_addressed = f_runoff.createVariable ( 'dst_lat_addressed', 'f8', ('num_links',) )
+v_oce_grid_pmask     [:] = oce_grid_pmask[:].ravel()
+
+## For diag, not used by OASIS
+v_atm_lon_addressed   = f_runoff.createVariable ( 'src_lon_addressed'  , 'f8', ('num_links',) )
+v_atm_lat_addressed   = f_runoff.createVariable ( 'src_lat_addressed'  , 'f8', ('num_links',) )
+v_atm_area_addressed  = f_runoff.createVariable ( 'src_area_addressed' , 'f8', ('num_links',) )
+v_atm_imask_addressed = f_runoff.createVariable ( 'src_imask_addressed', 'i4', ('num_links',) )
+v_atm_pmask_addressed = f_runoff.createVariable ( 'src_pmask_addressed', 'i4', ('num_links',) )
+
+v_oce_lon_addressed   = f_runoff.createVariable ( 'dst_lon_addressed'  , 'f8', ('num_links',) )
+v_oce_lat_addressed   = f_runoff.createVariable ( 'dst_lat_addressed'  , 'f8', ('num_links',) )
+v_oce_area_addressed  = f_runoff.createVariable ( 'dst_area_addressed' , 'f8', ('num_links',) )
+v_oce_imask_addressed = f_runoff.createVariable ( 'dst_imask_addressed', 'i4', ('num_links',) )
+v_oce_pmask_addressed = f_runoff.createVariable ( 'dst_pmask_addressed', 'i4', ('num_links',) )
 
 v_atm_lon_addressed.long_name="Longitude" ; v_atm_lon_addressed.standard_name="longitude" ; v_atm_lon_addressed.units="degrees_east"
-v_atm_lat_addressed.long_name="Latitude"  ; v_atm_lat_addressed.standard_name="latitude"  ; v_atm_lat_addressed.units="degrees_east"
-v_atm_lon_addressed [:] = atm_grid_center_lon.ravel()[atm_address].ravel()
-v_atm_lat_addressed [:] = atm_grid_center_lat.ravel()[atm_address].ravel()
+v_atm_lat_addressed.long_name="Latitude"  ; v_atm_lat_addressed.standard_name="latitude"  ; v_atm_lat_addressed.units="degrees_north"
+v_atm_lon_addressed  [:] = atm_grid_center_lon.ravel()[atm_address].ravel()
+v_atm_lat_addressed  [:] = atm_grid_center_lat.ravel()[atm_address].ravel()
+v_atm_area_addressed [:] = atm_grid_area.ravel()[atm_address].ravel()
+v_atm_imask_addressed[:] = 1-atm_grid_imask.ravel()[atm_address].ravel()
+v_atm_pmask_addressed[:] = 1-atm_grid_pmask.ravel()[atm_address].ravel()
 
 v_oce_lon_addressed.long_name="Longitude" ; v_oce_lon_addressed.standard_name="longitude" ; v_oce_lon_addressed.units="degrees_east"
-v_oce_lat_addressed.long_name="Latitude"  ; v_oce_lat_addressed.standard_name="latitude"  ; v_oce_lat_addressed.units="degrees_east"
-v_oce_lon_addressed [:] = oce_grid_center_lon.ravel()[oce_address].ravel()
-v_oce_lat_addressed [:] = oce_grid_center_lat.ravel()[oce_address].ravel()
+v_oce_lat_addressed.long_name="Latitude"  ; v_oce_lat_addressed.standard_name="latitude"  ; v_oce_lat_addressed.units="degrees_north"
+v_oce_lon_addressed  [:] = oce_grid_center_lon.ravel()[oce_address].ravel()
+v_oce_lat_addressed  [:] = oce_grid_center_lat.ravel()[oce_address].ravel()
+v_oce_area_addressed [:] = oce_grid_area.ravel()[oce_address].ravel()
+v_oce_imask_addressed[:] = 1-oce_grid_imask.ravel()[oce_address].ravel()
+v_oce_pmask_addressed[:] = 1-oce_grid_pmask.ravel()[oce_address].ravel()
+
+v_atmLand         = f_runoff.createVariable ( 'atmLand'        , 'i4', ('src_grid_size',) )
+v_atmLandFiltered = f_runoff.createVariable ( 'atmLandFiltered', 'f4', ('src_grid_size',) )
+v_atmLandFrac     = f_runoff.createVariable ( 'atmLandFrac'    , 'i4', ('src_grid_size',) )
+v_atmFrac         = f_runoff.createVariable ( 'atmFrac'        , 'i4', ('src_grid_size',) )
+v_atmOcean        = f_runoff.createVariable ( 'atmOcean'       , 'i4', ('src_grid_size',) )
+v_atmOceanFrac    = f_runoff.createVariable ( 'atmOceanFrac'   , 'i4', ('src_grid_size',) )
+v_atmCoast        = f_runoff.createVariable ( 'atmCoast'       , 'i4', ('src_grid_size',) ) 
+
+v_atmLand[:]         = atmLand.ravel()
+v_atmLandFrac[:]     = atmLandFrac.ravel()
+v_atmLandFiltered[:] = atmLandFiltered.ravel()
+v_atmFrac[:]         = atmFrac.ravel()
+v_atmOcean[:]        = atmOcean.ravel()
+v_atmOceanFrac[:]    = atmOceanFrac.ravel()
+v_atmCoast[:]        = atmCoast.ravel()
+
+v_oceLand          = f_runoff.createVariable ( 'oceLand'         , 'i4', ('dst_grid_size',) )
+v_oceOcean         = f_runoff.createVariable ( 'oceOcean'        , 'i4', ('dst_grid_size',) )
+v_oceOceanFiltered = f_runoff.createVariable ( 'oceOceanFiltered', 'f4', ('dst_grid_size',) )
+v_oceCoast         = f_runoff.createVariable ( 'oceCoast'        , 'i4', ('dst_grid_size',) )
+
+v_oceLand[:]      = oceLand.ravel()
+v_oceOcean[:]     = oceOcean.ravel()
+v_oceOceanFiltered[:]     = oceOceanFiltered.ravel()
+v_oceCoast[:]     = oceCoast.ravel()
 
 f_runoff.close ()