source: TOOLS/MOSAIX/ComputeNemoCoast.py @ 3725

### ===========================================================================
###
### Compute a mask for coastal point of NEMO domain
###
### ===========================================================================
##
##  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$"
__Date__     = "$Date$"
__Revision__ = "$Revision$"
__Id__       = "$Id$"
__HeadURL    = "$HeadURL$"
##
import netCDF4
import nemo
import numpy as np
import getopt, sys

def usage () :
    texte = """%(prog)s usage :
python %(prog)s [-d] [-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=4, 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()