source: TOOLS/CPLRESTART/create_flxat.py @ 6764

Last change on this file since 6764 was 6669, checked in by omamce, 7 months ago

O.M. : CPLRESTART

  • More comments
  • Python cleanup
  • More use of nemo.py functionalities
File size: 14.5 KB
Line 
1### ===========================================================================
2###
3### Part of CPL Restart IPSL packages
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###
17### Documentation : https://forge.ipsl.jussieu.fr/igcmg/wiki/IPSLCM6/MOSAIX
18###
19## SVN information
20#  $Author: omamce $
21#  $Date: 2020-09-18 17:02:09 +0200 (Fri, 18 Sep 2020) $
22#  $Revision: 5157 $
23#  $Id: CreateRestartAtm4Oasis.bash 5157 2020-09-18 15:02:09Z omamce $
24#  $HeadURL: svn+ssh://omamce@forge.ipsl.jussieu.fr/ipsl/forge/projets/igcmg/svn/TOOLS/CPLRESTART/CreateRestartAtm4Oasis.bash $
25#
26'''Create atmospheric restart
27
28## SVN information
29Author   = "$Author: omamce $"
30Date     = "$Date: 2023-10-25 17:11:15 +0200 (Wed, 25 Oct 2023) $"
31Revision = "$Revision: 6666 $"
32Id       = "$Id: RunoffWeights.py 6666 2023-10-25 15:11:15Z omamce $"
33HeadURL  = "$HeadURL: svn+ssh://omamce@forge.ipsl.jussieu.fr/ipsl/forge/projets/igcmg/svn/TOOLS/MOSAIX/RunoffWeights.py $"
34'''
35
36# SVN information
37__SVN__ = ({
38    'Author'   : '$Author: omamce $',
39    'Date'     : '$Date: 2023-10-25 17:11:15 +0200 (Wed, 25 Oct 2023) $',
40    'Revision' : '$Revision: 6666 $',
41    'Id'       : '$Id: RunoffWeights.py 6666 2023-10-25 15:11:15Z omamce $',
42    'HeadURL'  : '$HeadURL: svn+ssh://omamce@forge.ipsl.jussieu.fr/ipsl/forge/projets/igcmg/svn/TOOLS/MOSAIX/RunoffWeights.py $',
43    'SVN_Date' : '$SVN_Date: $',
44    })
45##
46
47import argparse
48import textwrap
49import time
50import os
51import sys
52import platform
53import numpy as np
54import xarray as xr
55
56## Reading the command line arguments
57#
58# Creating a parser
59# The first step in using the argparse is creating an ArgumentParser object:
60parser = argparse.ArgumentParser (description = textwrap.dedent ("""
61   create_flxat
62   """), epilog='-------- This is the end of the help message --------')
63
64# Adding arguments
65parser.add_argument ('-v', '--verbose', help='verbosity', action="store_true",
66                         default=False )
67parser.add_argument ('-l', '--level', type=int, default=0,
68                         choices=[0, 1, 2], help='verbosity level')
69parser.add_argument ('-i', '--input'  , metavar='inputFile'  , nargs='?',
70                         default="flxat_fields_notime.nc",
71                         help="input file")
72parser.add_argument ('-o', '--output' , metavar='outputFile', nargs='?',
73                         default="tmp_flxat.nc",
74                         help="output file" )
75parser.add_argument ('-f', '--format' , metavar='NetCDF_Format', nargs='?',
76                         default="NETCDF4",
77                         choices=["NETCDF4","NETCDF4_CLASSIC", "NETCDF3_64BIT",
78                                   "NETCDF3_CLASSIC", "64bits"],
79                         help="NetCDF format")
80parser.add_argument ('--IsUnstructured',  choices=['True', 'False'],
81                         required=True )
82
83# Parse command line
84myargs = parser.parse_args ()
85
86IsUnstructured = myargs.IsUnstructured in [ True, 'True', 'true', 'TRUE' ]
87NcFormat = myargs.format
88if NcFormat == '64bit' : NcFormat = 'NETCDF4'
89
90## Read Data - promote them to 64 bits
91d_In = xr.open_dataset (myargs.input)
92
93lon = d_In.lon.astype(float).values
94lat = d_In.lat.astype(float).values
95
96if IsUnstructured :
97    dims = lon.shape
98else :
99    dims = ( lat.shape[0], lon.shape[0] )
100
101# Try to read variables
102# Set them to zero if not possible
103if 'evap_oce' in d_In.variables  :
104    evap_oce  = d_In.evap_oce[0].squeeze().astype(float).values
105else                             :
106    evap_oce  = np.zeros ( dims )
107if 'evap_sic' in d_In.variables  :
108    evap_sic  = d_In.evap_sic[0].squeeze().astype(float).values
109else :
110    evap_sic  = np.zeros ( dims )
111if 'fract_oce' in d_In.variables :
112    fract_oce = d_In.fract_oce[0].squeeze().astype(float).values
113else :
114    fract_oce = np.ones ( dims )
115if 'fract_sic' in d_In.variables :
116    fract_sic = d_In.fract_sic[0].squeeze().astype(float).values
117else :
118    fract_sic = np.zeros ( dims )
119if 'precip' in d_In.variables    :
120    precip    = d_In.precip[0].squeeze().astype(float).values
121else :
122    evap_oce  = np.zeros ( dims )
123if 'snow' in d_In.variables      :
124    snow      = d_In.snow[0].squeeze().astype(float).values
125else :
126    snow      = np.zeros ( dims )
127if 'soll' in d_In.variables      :
128    soll      = d_In.soll[0].squeeze().astype(float).values
129else :
130    soll      = np.zeros ( dims )
131if 'sols' in d_In.variables      :
132    sols      = d_In.sols[0].squeeze().astype(float).values
133else :
134    sols      = np.zeros ( dims )
135if 'taux_oce' in d_In.variables  :
136    taux_oce  = d_In.taux_oce[0].squeeze().astype(float).values
137else :
138    taux_oce  = np.zeros ( dims )
139if 'taux_sic' in d_In.variables  :
140    taux_sic  = d_In.taux_sic[0].squeeze().astype(float).values
141else :
142    taux_sic  = np.zeros ( dims )
143if 'tauy_oce' in d_In.variables  :
144    tauy_oce  = d_In.tauy_oce[0].squeeze().astype(float).values
145else :
146    tauy_oce  = np.zeros ( dims )
147if 'tauy_sic' in d_In.variables  :
148    tauy_sic  = d_In.tauy_sic[0].squeeze().astype(float).values
149else :
150    tauy_sic  = np.zeros ( dims )
151if 'wind10m' in d_In.variables   :
152    wind10m   = d_In.wind10m[0].squeeze().astype(float).values
153else :
154    wind10m   = np.zeros ( dims )
155
156if IsUnstructured :
157    lon       = np.expand_dims ( lon       , axis=1 )
158    lat       = np.expand_dims ( lat       , axis=1 )
159    evap_oce  = np.expand_dims ( evap_oce  , axis=1 )
160    evap_sic  = np.expand_dims ( evap_sic  , axis=1 )
161    fract_oce = np.expand_dims ( fract_oce , axis=1 )
162    fract_sic = np.expand_dims ( fract_sic , axis=1 )
163    precip    = np.expand_dims ( precip    , axis=1 )
164    snow      = np.expand_dims ( snow      , axis=1 )
165    soll      = np.expand_dims ( soll      , axis=1 )
166    sols      = np.expand_dims ( sols      , axis=1 )
167    taux_oce  = np.expand_dims ( taux_oce  , axis=1 )
168    taux_sic  = np.expand_dims ( taux_sic  , axis=1 )
169    tauy_oce  = np.expand_dims ( tauy_oce  , axis=1 )
170    tauy_sic  = np.expand_dims ( tauy_sic  , axis=1 )
171    wind10m   = np.expand_dims ( wind10m   , axis=1 )
172else :
173    lon2 = lat [:, np.newaxis]*0 + lon [np.newaxis, :]
174    lat2 = lat [:, np.newaxis]   + lon [np.newaxis, :]*0
175    lon = lon2
176    lat = lat2
177
178##
179yxshape = lat.shape
180ny, nx = yxshape
181
182## Computations
183np.seterr (divide='ignore', invalid='ignore')
184
185## ocean fraction
186fract_oce_plus_sic  = fract_oce + fract_sic
187# free ocean vs. total ocen fraction
188fract_oce_norm = np.where (fract_oce_plus_sic >  0.0,
189                               fract_oce/fract_oce_plus_sic, 0.0)
190 # sea ice vs. total ocean fraction
191fract_sic_norm = np.where (fract_oce_plus_sic >  0.0,
192                               fract_sic/fract_oce_plus_sic, 0.0)
193##
194COTOTRAI = xr.DataArray (precip-snow, dims=('y', 'x'))
195COTOTRAI.attrs['long_name'] = 'Liquid precip'
196COTOTRAI.attrs['coordinates'] = "lat lon"
197
198COTOTSNO = xr.DataArray (snow       , dims=('y', 'x'))
199COTOTSNO.attrs['long_name'] ='Solid precipitation'
200COTOTSNO.attrs['coordinates'] = "lat lon"
201
202COTOTEVA = xr.DataArray (evap_oce*fract_oce_norm + evap_sic*fract_sic_norm,
203                             dims=('y', 'x'))
204COTOTEVA.attrs['coordinates'] = "lat lon"
205
206COICEVAP = xr.DataArray (evap_sic   , dims=('y', 'x'))
207COICEVAP.attrs['long_name'] = 'Evaporation on sea ice'
208COICEVAP.attrs['coordinates'] = "lat lon"
209
210COQSRMIX = xr.DataArray (sols       , dims=('y', 'x'))
211COQSRMIX.attrs['long_name'] = 'Heat flux short wave'
212COQSRMIX.attrs['units'] = 'W/m2'
213COQSRMIX.attrs['coordinates'] = "lat lon"
214
215COQNSMIX = xr.DataArray (soll       , dims=('y', 'x'))
216COQNSMIX.attrs['long_name'] = 'Heat flux minus short wave'
217COQNSMIX.attrs['units'] = 'W/m2'
218COQNSMIX.attrs['coordinates'] = "lat lon"
219
220COSHFICE = xr.DataArray (sols       , dims=('y', 'x'))
221COSHFICE.attrs['long_name'] = 'Heat flux short wave over sea ice'
222COSHFICE.attrs['units'] = 'W/m2'
223COSHFICE.attrs['coordinates'] = "lat lon"
224
225CONSFICE = xr.DataArray (soll       , dims=('y', 'x'))
226CONSFICE.attrs['long_name'] = 'Heat flux minus short wave over sea ice'
227CONSFICE.attrs['units'] = 'W/m2'
228CONSFICE.attrs['coordinates'] = "lat lon"
229
230COWINDSP = xr.DataArray (wind10m    , dims=('y', 'x'))
231COWINDSP.attrs['long_name'] = 'Wind speed at 10m high'
232COWINDSP.attrs['units'] = 'm/s'
233COWINDSP.attrs['coordinates'] = "lat lon"
234
235CODFLXDT = xr.DataArray (-20.0*np.ones(yxshape), dims=('y', 'x'))
236CODFLXDT.attrs['long_name'] = 'dQ/dT - Derivative over temperature of turbulent heat fluxes'
237CODFLXDT.attrs['units'] = 'W/m2/K'
238CODFLXDT.attrs['coordinates'] = "lat lon"
239
240COCALVIN = xr.DataArray (  0.0*np.ones(yxshape), dims=('y', 'x'))
241COCALVIN.attrs['long_name'] = 'Calving of icebergs, solid'
242COCALVIN.attrs['coordinates'] = "lat lon"
243
244COLIQRUN = xr.DataArray (  0.0*np.ones(yxshape), dims=('y', 'x'))
245COLIQRUN.attrs['long_name'] = 'River run-off, liquid'
246COLIQRUN.attrs['coordinates'] = "lat lon"
247
248## Wind stress
249tau_x = taux_oce*fract_oce_norm + taux_sic*fract_sic_norm
250tau_y = tauy_oce*fract_oce_norm + tauy_sic*fract_sic_norm
251COTAUMOD = xr.DataArray (np.sqrt ( (tau_x*tau_x + tau_y*tau_y) ) , dims=('y', 'x'))
252COTAUMOD.attrs.update ( {'long_name':'Wind stress modulus', 'units':'Pa',
253                             'coordinates':'lat lon' })
254
255## Wind stress, from east/north components to geocentric
256rad = np.deg2rad (1.0)
257COTAUXXU = xr.DataArray (-tau_x * np.sin(rad * lon)
258                        - tau_y * np.sin(rad * lat) * np.cos(rad * lon) ,
259                             dims=('y', 'x'))
260COTAUYYU = xr.DataArray ( tau_x * np.cos(rad * lon)
261                        - tau_y * np.sin(rad * lat) * np.sin(rad * lon) ,
262                              dims=('y', 'x'))
263COTAUZZU = xr.DataArray ( tau_y * np.cos(rad * lat)                                                 , dims=('y', 'x'))
264
265## Value at North Pole
266if IsUnstructured :
267    ## Value at North Pole for DYNAMICO grid
268    COTAUXXU = xr.where ( lat >= 89.999, -tau_y, COTAUXXU)
269    COTAUYYU = xr.where ( lat >= 89.999,  tau_x, COTAUYYU)
270    ## Value at South Pole for DYNAMICO grid ?
271
272else :
273    ## Value at North Pole for LMDZ lon/lat grid
274    COTAUXXU[0,:] = ( -tau_x [0, 0] )
275    COTAUYYU[0,:] = ( -tau_y [0, 0] )
276    COTAUZZU[0,:] =  0.0
277    ## Value at south Pole for LMDZ lon/lat grid
278    COTAUXXU[-1,:] = ( -tau_x [-1, 0] )
279    COTAUYYU[-1,:] = ( -tau_y [-1, 0] )
280
281##
282COTAUXXU.attrs.update ({'long_name':'Wind stress in geocentric referential - x-component',
283                            'units':'Pa', 'coordinates':'lat lon', 'Grid':'U' })
284COTAUYYU.attrs.update ({'long_name':'Wind stress in geocentric referential - y-component',
285                            'units':'Pa', 'coordinates':'lat lon', 'Grid':'U' })
286COTAUZZU.attrs.update ({'long_name':'Wind stress in geocentric referential - z-component',
287                            'units':'Pa', 'coordinates':'lat lon', 'Grid':'U' })
288
289COTAUXXV = COTAUXXU
290COTAUYYV = COTAUYYU
291COTAUZZV = COTAUZZU
292
293COTAUXXV.attrs.update ( {'Grid':'V'} )
294COTAUYYV.attrs.update ( {'Grid':'V'} )
295COTAUZZV.attrs.update ( {'Grid':'V'} )
296
297## check if bounds for lon lat are present and add them to dataset or drop them
298ifbnds = 'bounds_lon' in d_In.data_vars and 'bounds_lat' in d_In.data_vars
299
300## Creates final Dataset
301lon = xr.DataArray (lon, dims=('y', 'x'))
302lon.attrs.update ({'name':'longitude', 'long_name':'Longitude',
303                       'units':'degrees_east', 'standard_name':'longitude' })
304
305lat = xr.DataArray (lat, dims=('y', 'x'))
306lat.attrs.update ({'name':'latitude', 'long_name':'Latitude',
307                       'units':'degrees_north', 'standard_name':'latitude' })
308
309if ifbnds :
310    lon.attrs['bounds'] = 'bounds_lon'
311    lat.attrs['bounds'] = 'bounds_lat'
312
313    bounds_lon = d_In.bounds_lon.values.astype (float)
314    bounds_lat = d_In.bounds_lat.values.astype (float)
315    nvertex = bounds_lon.shape[-1]
316
317    bounds_lon = xr.DataArray ( np.reshape (bounds_lon, (ny, nx, nvertex)),
318                                    dims=('y', 'x', 'nvertex'))
319    bounds_lat = xr.DataArray ( np.reshape (bounds_lat, (ny, nx, nvertex)),
320                                    dims=('y', 'x', 'nvertex'))
321
322    bounds_lon.attrs['units'] = 'degrees_east'
323    bounds_lat.attrs['units'] = 'degrees_north'
324
325# Prepares dictionnary to export dataset to netcdf file with or without bounds
326dictdataset = {'lat':lat, 'lon':lon }
327if ifbnds:
328    dictdataset.update ( {'bounds_lon':bounds_lon,
329                          'bounds_lat':bounds_lat} )
330dictdataset.update ( {
331    'COTOTRAI':COTOTRAI, 'COTOTSNO':COTOTSNO, 'COTOTEVA':COTOTEVA,
332    'COICEVAP':COICEVAP, 'COQSRMIX':COQSRMIX, 'COQNSMIX':COQNSMIX,
333    'COSHFICE':COSHFICE, 'CONSFICE':CONSFICE, 'CODFLXDT':CODFLXDT,
334    'COCALVIN':COCALVIN, 'COLIQRUN':COLIQRUN, 'COWINDSP':COWINDSP,
335    'COTAUMOD':COTAUMOD, 'COTAUXXU':COTAUXXU, 'COTAUYYU':COTAUYYU,
336    'COTAUZZU':COTAUZZU, 'COTAUXXV':COTAUXXV, 'COTAUYYV':COTAUYYV,
337    'COTAUZZV':COTAUZZV } )
338
339d_Out = xr.Dataset (dictdataset)
340
341d_Out.attrs.update ( {
342    'AssociatedFiles' : myargs.input,
343    'Conventions'     : "CF-1.6",
344    'source'          : "IPSL Earth system model",
345    'group'           : "ICMC IPSL Climate Modelling Center",
346    'Institution'     : "IPSL https://www.ipsl.fr",
347    'Model'           : "IPSL CM",
348    'description'     : "Fields needed by OASIS-MCT",
349    'Program'         : f'Generated by {sys.argv[0]} with flags ' + ' '.join (sys.argv[1:]),
350    'timeStamp'       : time.asctime (),
351    'directory'       : os.getcwd (),
352    'HOSTNAME'        : platform.node (),
353    'LOGNAME'         : os.getlogin (),
354    'Python'          : "Python version " +  platform.python_version (),
355    'OS'              : platform.system (),
356    'release'         : platform.release (),
357    'hardware'        : platform.machine (),
358    'SVN_Author'      : "$Author: omamce $",
359    'SVN_Date'        : "$Date: 2022-07-06 11:06:07 +0200 (Wed, 06 Jul 2022) $",
360    'SVN_Revision'    : "$Revision: 6190 $",
361    'SVN_Id'          : "$Id: CalvingWeights.py 6190 2022-07-06 09:06:07Z omamce $",
362    'SVN_HeadURL'     : "$HeadURL: svn+ssh://omamce@forge.ipsl.jussieu.fr/ipsl/forge/projets/igcmg/svn/TOOLS/MOSAIX/CalvingWeights.py $",
363    } )
364
365d_Out.to_netcdf ( myargs.output, format=NcFormat )
366d_Out.close ()
367## ===========================================================================
368##
369##                               That's all folk's !!!
370##
371## ===========================================================================
Note: See TracBrowser for help on using the repository browser.