source: TOOLS/MOSAIX/ComputeNemoCoast.py

### ===========================================================================
###
### Compute a mask for coastal point of NEMO domain
###
### ===========================================================================
##
##  MOSAIX is under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt"
##  file for an english version of the licence and
##  "Licence_CeCILL_V2-fr.txt" for a french version.
##
##  Permission is hereby granted, free of charge, to any person or
##  organization obtaining a copy of the software and accompanying
##  documentation covered by this license (the "Software") to use,
##  reproduce, display, distribute, execute, and transmit the
##  Software, and to prepare derivative works of the Software, and to
##  permit third-parties to whom the Software is furnished to do so,
##  all subject to the following:
##
##  Warning, to install, configure, run, use any of MOSAIX 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).  Authors 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
##
##
## SVN information
__Author__   = "$Author$"
__Date__     = "$Date$"
__Revision__ = "$Revision$"
__Id__       = "$Id$"
__HeadURL    = "$HeadURL$"
##
import netCDF4
import nemo
import numpy as np
import getopt, sys

def usage () :
    texte = """Compute a mask for coastal point of NEMO domain

    %(prog)s usage :
python %(prog)s [-h|--help] [-d|--debug] [-i <orca grid file>] [-n <perio>]
 -d         | --debug         : debug
 -i <file>  | --input=<file>  : input file  (default: none)
 -n <perio> | --perio=<perio> : periodicity type (default: try to guess)
    #   1, 4, 6 : Cyclic on i dimension (generaly longitudes)
    #   2       : Obsolete (was symmetric condition at southern boundary ?)
    #   3, 4    : North fold T-point pivot (legacy ORCA2)
    #   5, 6    : North fold F-point pivot (ORCA1, ORCA025, ORCA2 with new grid for paleo)
    If <perio> is no specified,  %(prog)s will try to guess it from the grid dimension 
example :
    python %(prog)s -n 4 -i ORCA2.3_coordinates_mask.nc
    """
    print ( texte % { 'prog':sys.argv[0] } )
    
## Default input parameters
nperio      = None
GridFile    = None
Debug       = False

## Command line options
try:
    myopts, myargs = getopt.getopt ( sys.argv[1:], 'i:n:h', [ 'input=', 'nperio=', 'debug=', 'help' ] )
except getopt.GetoptError as cmdle :
    print ( "Command line error : "+str(cmdle)+"\n" )
    usage ()
    sys.exit(1)

for myopt, myval in myopts :
    if myopt in [ '-h', '--help'   ] : usage () ; sys.exit (0) ;
    if myopt in [ '-i', '--input'  ] : GridFile = myval ;
    if myopt in [ '-n', '--nperio' ] : nperio   = int(myval) ;
    if myopt in [ '-d', '--debug'  ] : Debug    = True ;
##
if GridFile == None :
    print ("Input grid file not specified\n")
    usage ( )
    sys.exit(-1)

print ("Input file : " + GridFile)   

## Open grid file
GridFile = netCDF4.Dataset ( GridFile, "r+" )

if nperio == None : # Try to get periodicity for the grid
    print ("Trying to guess nperio parameter")
    jpoi = GridFile.dimensions["x_grid_T"].size
    jpoj = GridFile.dimensions["y_grid_T"].size
    print ("Grid dimensions: ("+str(jpoj)+", "+str(jpoi)+")")
    if (jpoj, jpoi) == (149, 182):
        print ("ORCA 2 grid found: nperio may vary for this configuration")
    if (jpoj, jpoi) == (332, 292):
        nperio = 6
        print ("ORCA1 grid found" )
    if (jpoj, jpoi) == (332, 362):
        nperio = 6
        print ("eORCA1 grid found" )
####    
if nperio == None :
    print ("Periodicity not specified and not found\n")
    usage ( )
    sys.exit(-1)

###########################################################
    
print ("Periodicity : "+ str(nperio) )

if nperio in  [1, 4, 6] : print ("Cyclic on i dimension (generaly longitudes)")
if nperio in  [2,     ] : print ("Obsolete (was symmetric condition at southern boundary ?)")
if nperio in  [3, 4   ] : print ("North fold T-point pivot (legacy ORCA2)")
if nperio in  [5, 6   ] : print ("North fold F-point pivot (ORCA1, ORCA025, ORCA2 with new grid for paleo)")

###    
OceMask = GridFile.variables["maskutil_T"]
OceMask = np.array ( OceMask)

OceMask = nemo.lbc ( np.int_(OceMask[:,:]), nperio=nperio, cd_type='T' )

Temp = OceMask*0
Temp[1:-1, 1:-1] =   \
       + OceMask[0:-2, 1:-1] + OceMask[2:  , 1:-1] \
       + OceMask[1:-1, 0:-2] + OceMask[1:-1, 2:  ] \
       + OceMask[0:-2, 0:-2] + OceMask[1:-1, 0:-2] \
       + OceMask[1:-1, 1:-1] + OceMask[0:-2, 1:-1] 

CoastCrit = Temp.max()
print ("Maximum number of neighbors : "+str(CoastCrit) )

Temp = nemo.lbc ( Temp,  nperio=nperio, cd_type='T' )

OceCoastal = np.where (OceMask == 1, True, False) * np.where (Temp < CoastCrit, True, False)
OceCoastal = np.where (OceCoastal, 1, 0)

varOceCoastal = GridFile.createVariable  ("OceCoastal", "i4", ( "y_grid_T", "x_grid_T") )
varOceCoastal.cell_measures = "area: area_grid_T"
varOceCoastal.coordinates   = "nav_lat_grid_T nav_lon_grid_T"

varOceCoastal[:,:] = OceCoastal[:,:]


GridFile.sync()
GridFile.close()