source: TOOLS/MOSAIX/ComputeNemoCoast.py @ 3740

Last change on this file since 3740 was 3725, checked in by omamce, 6 years ago

O.M. : first step to implement run-off

  • Property svn:keywords set to Date Revision HeadURL Author Id
File size: 4.5 KB
Line 
1### ===========================================================================
2###
3### Compute a mask for coastal point of NEMO domain
4###
5### ===========================================================================
6##
7##  Warning, to install, configure, run, use any of Olivier Marti's
8##  software or to read the associated documentation you'll need at least
9##  one (1) brain in a reasonably working order. Lack of this implement
10##  will void any warranties (either express or implied).
11##  O. Marti assumes no responsability for errors, omissions,
12##  data loss, or any other consequences caused directly or indirectly by
13##  the usage of his software by incorrectly or partially configured
14##  personal.
15##
16## SVN information
17__Author__   = "$Author$"
18__Date__     = "$Date$"
19__Revision__ = "$Revision$"
20__Id__       = "$Id$"
21__HeadURL    = "$HeadURL$"
22##
23import netCDF4
24import nemo
25import numpy as np
26import getopt, sys
27
28def usage () :
29    texte = """%(prog)s usage :
30python %(prog)s [-d] [-i <orca grid file>] [-n <perio>]
31 -d         | --debug         : debug
32 -i <file>  | --input=<file>  : input file  (default: none)
33 -n <perio> | --perio=<perio> : periodicity type (default: try to guess)
34    #   1, 4, 6 : Cyclic on i dimension (generaly longitudes)
35    #   2       : Obsolete (was symmetric condition at southern boundary ?)
36    #   3, 4    : North fold T-point pivot (legacy ORCA2)
37    #   5, 6    : North fold F-point pivot (ORCA1, ORCA025, ORCA2 with new grid for paleo)
38    If <perio> is no specified,  %(prog)s will try to guess it from the grid dimension
39example :
40    python %(prog)s -n 4 -i ORCA2.3_coordinates_mask.nc
41    """
42    print ( texte % { 'prog':sys.argv[0] } )
43   
44## Default input parameters
45nperio      = None
46GridFile    = None
47Debug       = False
48
49## Command line options
50try:
51    myopts, myargs = getopt.getopt ( sys.argv[1:], 'i:n:h', [ 'input=', 'nperio=', 'debug=', '--help' ] )
52except getopt.GetoptError as cmdle :
53    print ( "Command line error : "+str(cmdle)+"\n" )
54    usage ()
55    sys.exit(1)
56
57for myopt, myval in myopts :
58    if myopt in [ '-h', '--help'   ] : usage () ; sys.exit (0) ;
59    if myopt in [ '-i', '--input'  ] : GridFile = myval ;
60    if myopt in [ '-n', '--nperio' ] : nperio   = int(myval) ;
61    if myopt in [ '-d', '--debug'  ] : Debug    = True ;
62##
63if GridFile == None :
64    print ("Input grid file not specified\n")
65    usage ( )
66    sys.exit(-1)
67
68print ("Input file  :" + GridFile)   
69
70## Open grid file
71GridFile = netCDF4.Dataset ( GridFile, "r+" )
72
73if nperio == None : # Try to get periodicity for the grid
74    print ("Trying to guess nperio parameter")
75    jpoi = GridFile.dimensions["x_grid_T"].size
76    jpoj = GridFile.dimensions["y_grid_T"].size
77    print ("Grid dimensions: ("+str(jpoj)+", "+str(jpoi)+")")
78    if (jpoj, jpoi) == (149, 182):
79        print ("ORCA 2 grid found: nperio may vary for this configuration")
80    if (jpoj, jpoi) == (332, 292):
81        nperio = 6
82        print ("ORCA1 grid found" )
83    if (jpoj, jpoi) == (332, 362):
84        nperio = 6
85        print ("eORCA1 grid found" )
86####   
87if nperio == None :
88    print ("Periodicity not specified and not found\n")
89    usage ( )
90    sys.exit(-1)
91
92###########################################################
93   
94print ("Periodicity : "+ str(nperio) )
95
96if nperio in  [1, 4, 6] : print ("Cyclic on i dimension (generaly longitudes)")
97if nperio in  [2,     ] : print ("Obsolete (was symmetric condition at southern boundary ?)")
98if nperio in  [3, 4   ] : print ("North fold T-point pivot (legacy ORCA2)")
99if nperio in  [5, 6   ] : print ("North fold F-point pivot (ORCA1, ORCA025, ORCA2 with new grid for paleo)")
100
101###   
102OceMask = GridFile.variables["maskutil_T"]
103OceMask = np.array ( OceMask)
104
105OceMask = nemo.lbc ( np.int_(OceMask[:,:]), nperio=nperio, cd_type='T' )
106
107Temp = OceMask*0
108Temp[1:-1, 1:-1] =   \
109       + OceMask[0:-2, 1:-1] + OceMask[2:  , 1:-1] \
110       + OceMask[1:-1, 0:-2] + OceMask[1:-1, 2:  ] \
111       + OceMask[0:-2, 0:-2] + OceMask[1:-1, 0:-2] \
112       + OceMask[1:-1, 1:-1] + OceMask[0:-2, 1:-1] 
113
114CoastCrit = Temp.max()
115print ("Maximum number of neighbors : "+str(CoastCrit) )
116
117
118Temp = nemo.lbc ( Temp,  nperio=4, cd_type='T' )
119
120OceCoastal = np.where (OceMask == 1, True, False) * np.where (Temp < CoastCrit, True, False)
121OceCoastal = np.where (OceCoastal, 1, 0)
122
123varOceCoastal = GridFile.createVariable  ("OceCoastal", "i4", ( "y_grid_T", "x_grid_T") )
124varOceCoastal.cell_measures = "area: area_grid_T"
125varOceCoastal.coordinates   = "nav_lat_grid_T nav_lon_grid_T"
126
127varOceCoastal[:,:] = OceCoastal[:,:]
128
129
130GridFile.sync()
131GridFile.close()
Note: See TracBrowser for help on using the repository browser.