Changeset 4151


Ignore:
Timestamp:
11/22/18 11:45:17 (3 years ago)
Author:
omamce
Message:

O.M. : adding comments to cotes_etal

File:
1 edited

Legend:

Unmodified
Added
Removed
  • TOOLS/MOSAIX/cotes_etal.py

    r4148 r4151  
    11# -*- Mode: python -*- 
    2 import netCDF4 
    3 import numpy as np 
    4 import nemo 
    5 from scipy import ndimage 
    6 import sys, os, platform, argparse, textwrap, time 
    7  
    8 zero = np.dtype('float64').type(0.0) 
    9 zone = np.dtype('float64').type(1.0) 
    10 epsfrac = np.dtype('float64').type(1.0E-10) 
    11 pi  = np.pi 
    12 rad = pi/np.dtype('float64').type(180.0) 
    13 ra = np.dtype('float64').type(6400000.0) # Earth radius 
    14  
    15 def geodist (plon1, plat1, plon2, plat2) : 
    16       """Distance between two points (on sphere)""" 
    17       zs = np.sin (rad*plat1) * np.sin (rad*plat2) +  np.cos (rad*plat1) * np.cos (rad*plat2) * np.cos(rad*(plon2-plon1)) 
    18       zs = np.maximum (-zone, np.minimum (zone, zs)) 
    19       geodist =  np.arccos (zs) 
    20       return geodist 
    21  
     2### =========================================================================== 
    223### 
    23 ###   A partir d'un fichier de poids, complete avec les rivieres et 
    24 ###   les poids pour le run off 
    25 ###   On met a zero les poids sur les points non cotiers, et on renormalise 
     4### Compute runoff weights. 
     5### For LMDZ only. Not suitable for DYNAMICO 
    266### 
    27  
    28     
    29 ### Lecture de la ligne de comande 
     7### =========================================================================== 
     8## 
     9##  Warning, to install, configure, run, use any of Olivier Marti's 
     10##  software or to read the associated documentation you'll need at least 
     11##  one (1) brain in a reasonably working order. Lack of this implement 
     12##  will void any warranties (either express or implied). 
     13##  O. Marti assumes no responsability for errors, omissions, 
     14##  data loss, or any other consequences caused directly or indirectly by 
     15##  the usage of his software by incorrectly or partially configured 
     16##  personal. 
     17## 
    3018# SVN information 
    3119__Author__   = "$Author$" 
     
    3523__HeadURL__  = "$HeadURL$" 
    3624__SVN_Date__ = "$SVN_Date: $" 
    37  
    38 ### 
    39  
    40 ### ===== Handling command line parameters ================================================== 
     25## 
     26 
     27import netCDF4 
     28import numpy as np 
     29import nemo 
     30from scipy import ndimage 
     31import sys, os, platform, argparse, textwrap, time 
     32 
     33## Userful constants 
     34zero = np.dtype('float64').type(0.0) 
     35zone = np.dtype('float64').type(1.0) 
     36epsfrac = np.dtype('float64').type(1.0E-10) 
     37pi  = np.pi 
     38rad = pi/np.dtype('float64').type(180.0) 
     39ra = np.dtype('float64').type(6400000.0) # Earth radius 
     40 
     41## Functions 
     42def geodist (plon1, plat1, plon2, plat2) : 
     43      """Distance between two points (on sphere)""" 
     44      zs = np.sin (rad*plat1) * np.sin (rad*plat2) +  np.cos (rad*plat1) * np.cos (rad*plat2) * np.cos(rad*(plon2-plon1)) 
     45      zs = np.maximum (-zone, np.minimum (zone, zs)) 
     46      geodist =  np.arccos (zs) 
     47      return geodist 
     48 
     49 
     50    
     51 
     52### ===== Reading command line parameters ================================================== 
    4153# Creating a parser 
    4254parser = argparse.ArgumentParser ( 
    4355    description = """Compute calving weights""", 
    44     epilog='-------- This is the end of the help message --------') 
     56    epilog='-------- End of the help message --------') 
    4557 
    4658# Adding arguments 
     
    4860parser.add_argument ('--atm'   , help='atm model name (LMD*)', type=str, default='LMD9695'    ) 
    4961parser.add_argument ('--atmCoastWidth', help='width of the coastal band in atmosphere (in grid points)', default=1 ) 
    50 parser.add_argument ('--oceCoastWidth', help='width of the coastal band in ocean (in grid points)', default=1 ) 
     62parser.add_argument ('--oceCoastWidth', help='width of the coastal band in ocean (in grid points)'     , default=2 ) 
    5163parser.add_argument ('--grids', help='grids file name', default='grids.nc' ) 
    5264parser.add_argument ('--areas', help='masks file name', default='areas.nc' ) 
     
    5567parser.add_argument ('--output', help='output rmp file name', default='rmp_tlmd_to_torc_runoff.nc' ) 
    5668 
    57  
    5869# Parse command line 
    5970myargs = parser.parse_args() 
     
    6879atm_Name = myargs.atm 
    6980oce_Name = myargs.oce 
    70 # Largeur de la bande cotiere cote atm 
     81# Width of the coastal band (land points) in the atmopshere 
    7182atmCoastWidth = myargs.atmCoastWidth 
    72 # Largeur de la bande cotiere cote oce 
     83# Width of the coastal band (ocean points) in the ocean  
    7384oceCoastWidth = myargs.oceCoastWidth 
    7485 
     
    94105atm_grid_dims       = atm_grid_area.shape 
    95106(atm_nvertex, atm_jpj, atm_jpi) = atm_grid_corner_lat.shape 
    96 if atm_jpi == 1 : atm_perio = 0 # Icosahedre 
    97 else            : atm_perio = 1 # 
     107atm_perio = 0 
    98108atm_grid_pmask = nemo.lbc_mask (atm_grid_imask, 'T', atm_perio) 
    99109atm_address = np.reshape ( np.arange(atm_jpj*atm_jpi), (atm_jpj, atm_jpi) ) 
     
    115125oce_grid_size = oce_jpj*oce_jpi 
    116126 
    117 print ('Calculs points terre/oce/cote sur grille atmosphere' ) 
    118  
    119 # Fill closed sea with image processing library 
     127## Fill closed sea with image processing library 
    120128oce_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' ) 
    121129 
     130##  
    122131print ("Determination d'une bande cotiere ocean") 
    123132 
    124 oceLand  = np.where (oce_grid_imask[:] < 0.5, True, False) 
    125 oceOcean = np.where (oce_grid_imask[:] > 0.5, True, False) 
    126  
    127 oceCoast = np.where (ndimage.uniform_filter(oceOcean.astype(float), size=1+2*oceCoastWidth)<1.0,True,False) & oceOcean 
     133oceLand  = np.where (oce_grid_pmask[:] < 0.5, True, False) 
     134oceOcean = np.where (oce_grid_pmask[:] > 0.5, True, False) 
     135 
     136NNocean = 1+2*oceCoastWidth 
     137oceOceanFiltered = ndimage.uniform_filter(oceOcean.astype(float), size=NNocean) 
     138oceCoast = np.where (oceOceanFiltered<(1.0-0.5/(NNocean**2)),True,False) & oceOcean 
    128139oceCoast = nemo.lbc_mask (oceCoast, oce_perio, 'T') 
    129140 
    130 print ('Nombre de point oceLand  : ' + str(oceLand.sum()) ) 
    131 print ('Nombre de point oceOcean : ' + str(oceOcean.sum()) ) 
    132 print ('Nombre de point oceCoast : ' + str(oceCoast.sum()) ) 
    133  
    134 # Faire une liste de ces points 
     141print ('Number of points in oceLand  : ' + str(oceLand.sum()) ) 
     142print ('Number of points in oceOcean : ' + str(oceOcean.sum()) ) 
     143print ('Number of points in oceCoast : ' + str(oceCoast.sum()) ) 
     144 
     145# Arrays with coastal points only 
    135146oceCoast_grid_center_lon = oce_grid_center_lon[oceCoast] 
    136147oceCoast_grid_center_lat = oce_grid_center_lat[oceCoast] 
     
    142153 
    143154print ("Determination d'une bande cotiere atmosphere " ) 
    144 atmLand  = np.where (o2aFrac[:] < epsfrac, True, False) 
    145 atmFrac  = np.where (o2aFrac[:] < epsfrac, True, False) & np.where (o2aFrac[:] < (zone-epsfrac), True, False) 
    146 atmOcean = np.where (o2aFrac[:] < (zone-epsfrac), True, False) 
    147  
    148 atmCoast = np.where (ndimage.uniform_filter(atmLand.astype(float), size=1+2*atmCoastWidth)<1.0,True,False) & atmLand 
     155atmLand      = np.where (o2aFrac[:] < epsfrac       , True, False) 
     156atmLandFrac  = np.where (o2aFrac[:] < zone-epsfrac  , True, False) 
     157atmFrac      = np.where (o2aFrac[:] > epsfrac       , True, False) & np.where (o2aFrac[:] < (zone-epsfrac), True, False) 
     158atmOcean     = np.where (o2aFrac[:] < (zone-epsfrac), True, False) 
     159atmOceanFrac = np.where (o2aFrac[:] > epsfrac       , True, False) 
     160 
     161NNatm = 1+2*atmCoastWidth 
     162atmLandFiltered = ndimage.uniform_filter(atmLand.astype(float), size=NNatm) 
     163atmCoast = np.where (atmLandFiltered<(1.0-0.5/(NNatm**2)),True,False) & atmLandFrac 
    149164atmCoast = nemo.lbc_mask (atmCoast, 1, 'T') 
    150165 
    151 print ('Nombre de point atmLand  : ' + str(atmLand.sum()) ) 
    152 print ('Nombre de point atmOcean : ' + str(atmOcean.sum()) ) 
    153 print ('Nombre de point atmCoast : ' + str(atmCoast.sum()) ) 
    154  
    155 # Faire une liste de ces points 
     166print ('Number of points in atmLand  : ' + str(atmLand.sum()) ) 
     167print ('Number of points in atmOcean : ' + str(atmOcean.sum()) ) 
     168print ('Number of points in atmCoast : ' + str(atmCoast.sum()) ) 
     169 
     170# Arrays with coastal points only 
    156171atmCoast_grid_center_lon = atm_grid_center_lon[atmCoast] 
    157172atmCoast_grid_center_lat = atm_grid_center_lat[atmCoast] 
     
    165180oce_address  = np.empty ( shape=(0), dtype=np.int32   ) 
    166181 
     182## Loop on atmosphere coastal points 
    167183for ja in np.arange(len(atmCoast_grid_pmask)) : 
    168184    z_dist = geodist ( atmCoast_grid_center_lon[ja], atmCoast_grid_center_lat[ja], oceCoast_grid_center_lon, oceCoast_grid_center_lat) 
     
    185201         
    186202 
    187 print ('fini') 
     203print ('End of loop') 
    188204 
    189205num_links = remap_matrix.shape[0] 
    190206 
    191 # Output file 
     207### Output file 
    192208runoff = myargs.output 
    193209f_runoff = netCDF4.Dataset ( runoff, 'w', format='NETCDF3_64BIT' ) 
     
    203219f_runoff.associatedFiles = grids + " " + areas + " " + masks 
    204220f_runoff.directory       = os.getcwd () 
    205 f_runoff.description     = "Generated with ???" 
     221f_runoff.description     = "Generated with cotes_etal.py" 
    206222f_runoff.title           = runoff 
    207223f_runoff.Program         = "Generated by " + sys.argv[0] + " with flags " + str(sys.argv[1:]) 
     
    242258 
    243259v_remap_matrix[:] = remap_matrix 
    244 v_atm_address [:] = atm_address 
    245 v_oce_address [:] = oce_address 
     260v_atm_address [:] = atm_address + 1 # OASIS uses Fortran style indexing 
     261v_oce_address [:] = oce_address + 1 
    246262 
    247263v_atm_grid_dims       = f_runoff.createVariable ( 'src_grid_dims'      , 'i4', ('src_grid_rank',) ) 
     
    256272v_atm_grid_area       = f_runoff.createVariable ( 'src_grid_area'      , 'f8', ('src_grid_size',)  ) 
    257273v_atm_grid_area.long_name="Grid area" ; v_atm_grid_area.standard_name="cell_area" ; v_atm_grid_area.units="m2" 
    258 v_atm_grid_imask      = f_runoff.createVariable ( 'src_grid_imask'     , 'f8', ('src_grid_size',)  ) 
     274v_atm_grid_imask      = f_runoff.createVariable ( 'src_grid_imask'     , 'i4', ('src_grid_size',)  ) 
    259275v_atm_grid_imask.long_name="Land-sea mask" ; v_atm_grid_imask.units="Land:1, Ocean:0" 
     276v_atm_grid_pmask      = f_runoff.createVariable ( 'src_grid_pmask'     , 'i4', ('src_grid_size',)  ) 
     277v_atm_grid_pmask.long_name="Land-sea mask (periodicity removed)" ; v_atm_grid_pmask.units="Land:1, Ocean:0" 
    260278 
    261279v_atm_grid_dims      [:] = atm_grid_dims 
     
    266284v_atm_grid_area      [:] = atm_grid_area[:].ravel() 
    267285v_atm_grid_imask     [:] = atm_grid_imask[:].ravel() 
     286v_atm_grid_pmask     [:] = atm_grid_pmask[:].ravel() 
    268287 
    269288# -- 
     
    274293v_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" 
    275294v_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" 
    276 v_oce_grid_corner_lon = f_runoff.createVariable ( 'oce_grid_corner_lon', 'f8', ('dst_grid_size', 'dst_grid_corners')  ) 
     295v_oce_grid_corner_lon = f_runoff.createVariable ( 'dst_grid_corner_lon', 'f8', ('dst_grid_size', 'dst_grid_corners')  ) 
    277296v_oce_grid_corner_lat = f_runoff.createVariable ( 'dst_grid_corner_lat', 'f8', ('dst_grid_size', 'dst_grid_corners')  ) 
    278297v_oce_grid_corner_lon.units="degrees_east" 
    279298v_oce_grid_corner_lat.units="degrees_north" 
    280 v_oce_grid_area       = f_runoff.createVariable ( 'dst_grid_area'      , 'f8', ('dst_grid_size',) ) 
     299v_oce_grid_area       = f_runoff.createVariable ( 'dst_grid_area'  , 'f8', ('dst_grid_size',) ) 
    281300v_oce_grid_area.long_name="Grid area" ; v_oce_grid_area.standard_name="cell_area" ; v_oce_grid_area.units="m2" 
    282 v_oce_grid_imask      = f_runoff.createVariable ( 'dst_grid_imask'     , 'f8', ('dst_grid_size',)  ) 
     301v_oce_grid_imask      = f_runoff.createVariable ( 'dst_grid_imask'     , 'i4', ('dst_grid_size',)  ) 
    283302v_oce_grid_imask.long_name="Land-sea mask" ; v_oce_grid_imask.units="Land:1, Ocean:0" 
     303v_oce_grid_pmask      = f_runoff.createVariable ( 'dst_grid_pmask'     , 'i4', ('dst_grid_size',)  ) 
     304v_oce_grid_pmask.long_name="Land-sea mask (periodicity removed)" ; v_oce_grid_pmask.units="Land:1, Ocean:0" 
    284305 
    285306v_oce_grid_dims      [:] = oce_grid_dims 
     
    290311v_oce_grid_area      [:] = oce_grid_area[:].ravel() 
    291312v_oce_grid_imask     [:] = oce_grid_imask[:].ravel() 
    292  
    293 # For diags 
    294 v_atm_lon_addressed = f_runoff.createVariable ( 'src_lon_addressed', 'f8', ('num_links',) ) 
    295 v_atm_lat_addressed = f_runoff.createVariable ( 'src_lat_addressed', 'f8', ('num_links',) ) 
    296 v_oce_lon_addressed = f_runoff.createVariable ( 'dst_lon_addressed', 'f8', ('num_links',) ) 
    297 v_oce_lat_addressed = f_runoff.createVariable ( 'dst_lat_addressed', 'f8', ('num_links',) ) 
     313v_oce_grid_pmask     [:] = oce_grid_pmask[:].ravel() 
     314 
     315## For diag, not used by OASIS 
     316v_atm_lon_addressed   = f_runoff.createVariable ( 'src_lon_addressed'  , 'f8', ('num_links',) ) 
     317v_atm_lat_addressed   = f_runoff.createVariable ( 'src_lat_addressed'  , 'f8', ('num_links',) ) 
     318v_atm_area_addressed  = f_runoff.createVariable ( 'src_area_addressed' , 'f8', ('num_links',) ) 
     319v_atm_imask_addressed = f_runoff.createVariable ( 'src_imask_addressed', 'i4', ('num_links',) ) 
     320v_atm_pmask_addressed = f_runoff.createVariable ( 'src_pmask_addressed', 'i4', ('num_links',) ) 
     321 
     322v_oce_lon_addressed   = f_runoff.createVariable ( 'dst_lon_addressed'  , 'f8', ('num_links',) ) 
     323v_oce_lat_addressed   = f_runoff.createVariable ( 'dst_lat_addressed'  , 'f8', ('num_links',) ) 
     324v_oce_area_addressed  = f_runoff.createVariable ( 'dst_area_addressed' , 'f8', ('num_links',) ) 
     325v_oce_imask_addressed = f_runoff.createVariable ( 'dst_imask_addressed', 'i4', ('num_links',) ) 
     326v_oce_pmask_addressed = f_runoff.createVariable ( 'dst_pmask_addressed', 'i4', ('num_links',) ) 
    298327 
    299328v_atm_lon_addressed.long_name="Longitude" ; v_atm_lon_addressed.standard_name="longitude" ; v_atm_lon_addressed.units="degrees_east" 
    300 v_atm_lat_addressed.long_name="Latitude"  ; v_atm_lat_addressed.standard_name="latitude"  ; v_atm_lat_addressed.units="degrees_east" 
    301 v_atm_lon_addressed [:] = atm_grid_center_lon.ravel()[atm_address].ravel() 
    302 v_atm_lat_addressed [:] = atm_grid_center_lat.ravel()[atm_address].ravel() 
     329v_atm_lat_addressed.long_name="Latitude"  ; v_atm_lat_addressed.standard_name="latitude"  ; v_atm_lat_addressed.units="degrees_north" 
     330v_atm_lon_addressed  [:] = atm_grid_center_lon.ravel()[atm_address].ravel() 
     331v_atm_lat_addressed  [:] = atm_grid_center_lat.ravel()[atm_address].ravel() 
     332v_atm_area_addressed [:] = atm_grid_area.ravel()[atm_address].ravel() 
     333v_atm_imask_addressed[:] = 1-atm_grid_imask.ravel()[atm_address].ravel() 
     334v_atm_pmask_addressed[:] = 1-atm_grid_pmask.ravel()[atm_address].ravel() 
    303335 
    304336v_oce_lon_addressed.long_name="Longitude" ; v_oce_lon_addressed.standard_name="longitude" ; v_oce_lon_addressed.units="degrees_east" 
    305 v_oce_lat_addressed.long_name="Latitude"  ; v_oce_lat_addressed.standard_name="latitude"  ; v_oce_lat_addressed.units="degrees_east" 
    306 v_oce_lon_addressed [:] = oce_grid_center_lon.ravel()[oce_address].ravel() 
    307 v_oce_lat_addressed [:] = oce_grid_center_lat.ravel()[oce_address].ravel() 
     337v_oce_lat_addressed.long_name="Latitude"  ; v_oce_lat_addressed.standard_name="latitude"  ; v_oce_lat_addressed.units="degrees_north" 
     338v_oce_lon_addressed  [:] = oce_grid_center_lon.ravel()[oce_address].ravel() 
     339v_oce_lat_addressed  [:] = oce_grid_center_lat.ravel()[oce_address].ravel() 
     340v_oce_area_addressed [:] = oce_grid_area.ravel()[oce_address].ravel() 
     341v_oce_imask_addressed[:] = 1-oce_grid_imask.ravel()[oce_address].ravel() 
     342v_oce_pmask_addressed[:] = 1-oce_grid_pmask.ravel()[oce_address].ravel() 
     343 
     344v_atmLand         = f_runoff.createVariable ( 'atmLand'        , 'i4', ('src_grid_size',) ) 
     345v_atmLandFiltered = f_runoff.createVariable ( 'atmLandFiltered', 'f4', ('src_grid_size',) ) 
     346v_atmLandFrac     = f_runoff.createVariable ( 'atmLandFrac'    , 'i4', ('src_grid_size',) ) 
     347v_atmFrac         = f_runoff.createVariable ( 'atmFrac'        , 'i4', ('src_grid_size',) ) 
     348v_atmOcean        = f_runoff.createVariable ( 'atmOcean'       , 'i4', ('src_grid_size',) ) 
     349v_atmOceanFrac    = f_runoff.createVariable ( 'atmOceanFrac'   , 'i4', ('src_grid_size',) ) 
     350v_atmCoast        = f_runoff.createVariable ( 'atmCoast'       , 'i4', ('src_grid_size',) )  
     351 
     352v_atmLand[:]         = atmLand.ravel() 
     353v_atmLandFrac[:]     = atmLandFrac.ravel() 
     354v_atmLandFiltered[:] = atmLandFiltered.ravel() 
     355v_atmFrac[:]         = atmFrac.ravel() 
     356v_atmOcean[:]        = atmOcean.ravel() 
     357v_atmOceanFrac[:]    = atmOceanFrac.ravel() 
     358v_atmCoast[:]        = atmCoast.ravel() 
     359 
     360v_oceLand          = f_runoff.createVariable ( 'oceLand'         , 'i4', ('dst_grid_size',) ) 
     361v_oceOcean         = f_runoff.createVariable ( 'oceOcean'        , 'i4', ('dst_grid_size',) ) 
     362v_oceOceanFiltered = f_runoff.createVariable ( 'oceOceanFiltered', 'f4', ('dst_grid_size',) ) 
     363v_oceCoast         = f_runoff.createVariable ( 'oceCoast'        , 'i4', ('dst_grid_size',) ) 
     364 
     365v_oceLand[:]      = oceLand.ravel() 
     366v_oceOcean[:]     = oceOcean.ravel() 
     367v_oceOceanFiltered[:]     = oceOceanFiltered.ravel() 
     368v_oceCoast[:]     = oceCoast.ravel() 
    308369 
    309370f_runoff.close () 
Note: See TracChangeset for help on using the changeset viewer.