Changeset 6229 for TOOLS

Timestamp:

08/09/22 09:31:27 (21 months ago)

Author:

omamce

Message:

O.M: CPLRESTART. Corrections on create_flxat.py

File:

• TOOLS/CPLRESTART/create_flxat.py (modified) (3 diffs)

TOOLS/CPLRESTART/create_flxat.py

@@ -54 +54 @@
 d_In = xr.open_dataset (myargs.input)
 
-lon       = d_In.lon.astype(np.float64)
-lat       = d_In.lat.astype(np.float64)
-evap_oce  = d_In.evap_oce.astype(np.float64)
-evap_sic  = d_In.evap_sic.astype(np.float64)
-fract_oce = d_In.fract_oce.astype(np.float64)
-fract_sic = d_In.fract_sic.astype(np.float64)
-precip    = d_In.precip.astype(np.float64)
-snow      = d_In.snow.astype(np.float64)
-soll      = d_In.soll.astype(np.float64)
-sols      = d_In.sols.astype(np.float64)
-taux_oce  = d_In.taux_oce.astype(np.float64)
-taux_sic  = d_In.taux_sic.astype(np.float64)
-tauy_oce  = d_In.tauy_oce.astype(np.float64)
-tauy_sic  = d_In.tauy_sic.astype(np.float64)
-wind10m   = d_In.wind10m.astype(np.float64)
+lon       = d_In.lon.astype(np.float64).values
+lat       = d_In.lat.astype(np.float64).values
+evap_oce  = d_In.evap_oce.astype(np.float64).values
+evap_sic  = d_In.evap_sic.astype(np.float64).values
+fract_oce = d_In.fract_oce.astype(np.float64).values
+fract_sic = d_In.fract_sic.astype(np.float64).values
+precip    = d_In.precip.astype(np.float64).values
+snow      = d_In.snow.astype(np.float64).values
+soll      = d_In.soll.astype(np.float64).values
+sols      = d_In.sols.astype(np.float64).values
+taux_oce  = d_In.taux_oce.astype(np.float64).values
+taux_sic  = d_In.taux_sic.astype(np.float64).values
+tauy_oce  = d_In.tauy_oce.astype(np.float64).values
+tauy_sic  = d_In.tauy_sic.astype(np.float64).values
+wind10m   = d_In.wind10m.astype(np.float64).values
+
+##
+yxshape = lat.shape
+ny, nx = yxshape
 
 ## Computations
-
+np.seterr (divide='ignore', invalid='ignore')
 fract_oce_plus_sic  = (fract_oce + fract_sic) ; ## ocean fraction
-fract_oce_norm = xr.where (fract_oce_plus_sic >  0.0, (fract_oce)/fract_oce_plus_sic, 0.0) # free ocean vs. total ocen fraction
-fract_sic_norm = xr.where (fract_oce_plus_sic >  0.0, (fract_sic)/fract_oce_plus_sic, 0.0) # sea ice vs. total ocean fraction
-##
-COTOTRAI = precip-snow    ##  Liquid precip
-COTOTSNO = snow           ##  Solid precipitation
-COTOTEVA = evap_oce*fract_oce_norm + evap_sic*fract_sic_norm
-COICEVAP = evap_sic       ## Evaporation on sea ice
-COQSRMIX = sols           ## Heat flux short wave
-COQNSMIX = soll           ## Heat flux minus short wave
-COSHFICE = sols           ## Heat flux short wave over sea ice
-CONSFICE = soll           ## Heat flux minus short wave over sea ice
-CODFLXDT = -20.0          ## W/m2 - dQ/dt
-COCALVIN =   0.0          ## Calving of icebergs, solid
-COLIQRUN =   0.0          ## River run-off , liquid
-COWINDSP = wind10m        ## Wind speed at 10m high
-COTAUMOD =   0.0          ## Wind stress modulus
+fract_oce_norm = np.where (fract_oce_plus_sic >  0.0, fract_oce/fract_oce_plus_sic, 0.0) # free ocean vs. total ocen fraction
+fract_sic_norm = np.where (fract_oce_plus_sic >  0.0, fract_sic/fract_oce_plus_sic, 0.0) # sea ice vs. total ocean fraction
+##
+COTOTRAI = xr.DataArray (precip-snow, dims=('y', 'x'))
+COTOTRAI.attrs['long_name'] = 'Liquid precip'
+COTOTRAI.attrs['coordinates'] = "lat lon"
+
+COTOTSNO = xr.DataArray (snow       , dims=('y', 'x'))
+COTOTSNO.attrs['long_name'] ='Solid precipitation'
+COTOTSNO.attrs['coordinates'] = "lat lon"
+
+COTOTEVA = xr.DataArray (evap_oce*fract_oce_norm + evap_sic*fract_sic_norm, dims=('y', 'x'))  
+COTOTEVA.attrs['coordinates'] = "lat lon"
+
+COICEVAP = xr.DataArray (evap_sic   , dims=('y', 'x'))
+COICEVAP.attrs['long_name'] = 'Evaporation on sea ice'
+COICEVAP.attrs['coordinates'] = "lat lon"
+
+COQSRMIX = xr.DataArray (sols       , dims=('y', 'x'))
+COQSRMIX.attrs['long_name'] = 'Heat flux short wave'
+COQSRMIX.attrs['units'] = 'W/m2'
+COQSRMIX.attrs['coordinates'] = "lat lon"
+
+COQNSMIX = xr.DataArray (soll       , dims=('y', 'x'))
+COQNSMIX.attrs['long_name'] = 'Heat flux minus short wave'
+COQNSMIX.attrs['units'] = 'W/m2'
+COQNSMIX.attrs['coordinates'] = "lat lon"
+
+COSHFICE = xr.DataArray (sols       , dims=('y', 'x'))
+COSHFICE.attrs['long_name'] = 'Heat flux short wave over sea ice'
+COSHFICE.attrs['units'] = 'W/m2'
+COSHFICE.attrs['coordinates'] = "lat lon"
+
+CONSFICE = xr.DataArray (soll       , dims=('y', 'x'))
+CONSFICE.attrs['long_name'] = 'Heat flux minus short wave over sea ice'
+CONSFICE.attrs['units'] = 'W/m2'
+CONSFICE.attrs['coordinates'] = "lat lon"
+
+COWINDSP = xr.DataArray (wind10m    , dims=('y', 'x'))
+COWINDSP.attrs['long_name'] = 'Wind speed at 10m high'
+COWINDSP.attrs['units'] = 'm/s'
+COWINDSP.attrs['coordinates'] = "lat lon"
+
+CODFLXDT = xr.DataArray (-20.0*np.ones(yxshape), dims=('y', 'x'))
+CODFLXDT.attrs['long_name'] = 'dQ/dT - Derivative over temperatue of turbulent heat fluxes'
+CODFLXDT.attrs['units'] = 'W/m2/K'     
+CODFLXDT.attrs['coordinates'] = "lat lon"
+
+COCALVIN = xr.DataArray (  0.0*np.ones(yxshape), dims=('y', 'x'))
+COCALVIN.attrs['long_name'] = 'Calving of icebergs, solid'
+COCALVIN.attrs['coordinates'] = "lat lon"
+
+COLIQRUN = xr.DataArray (  0.0*np.ones(yxshape), dims=('y', 'x'))
+COLIQRUN.attrs['long_name'] = 'River run-off, liquid'
+COLIQRUN.attrs['coordinates'] = "lat lon"
 
 ## Wind stress
-tau_x = (taux_oce*fract_oce_norm + taux_sic*fract_sic_norm) ; 
-tau_y = (tauy_oce*fract_oce_norm + tauy_sic*fract_sic_norm) ;
-COTAUMOD = np.sqrt ( (tau_x*tau_x + tau_y*tau_y) )  ;
+tau_x = (taux_oce*fract_oce_norm + taux_sic*fract_sic_norm)
+tau_y = (tauy_oce*fract_oce_norm + tauy_sic*fract_sic_norm)
+COTAUMOD = xr.DataArray (np.sqrt ( (tau_x*tau_x + tau_y*tau_y) ) , dims=('y', 'x'))
+COTAUMOD.attrs['long_name'] = 'Wind stress modulus'
+COTAUMOD.attrs['units'] = 'Pa'
+COTAUMOD.attrs['coordinates'] = "lat lon"
 
 ## Wind stress, from east/north components to geocentric
-rad = np.deg2rad(1.0)
-COTAUXXU = -tau_x * np.sin(rad * lon) - tau_y * np.sin(rad * lat) * np.cos(rad * lon) 
-COTAUYYU =  tau_x * np.cos(rad * lon) - tau_y * np.sin(rad * lat) * np.sin(rad * lon)
-COTAUZZU =  tau_y * np.cos(rad * lat)
+rad = np.deg2rad (1.0)
+COTAUXXU = xr.DataArray (-tau_x * np.sin(rad * lon) - tau_y * np.sin(rad * lat) * np.cos(rad * lon) , dims=('y', 'x'))
+COTAUYYU = xr.DataArray ( tau_x * np.cos(rad * lon) - tau_y * np.sin(rad * lat) * np.sin(rad * lon) , dims=('y', 'x'))
+COTAUZZU = xr.DataArray ( tau_y * np.cos(rad * lat)                                                 , dims=('y', 'x'))
+
 
 ## Value at North Pole
-
 if myargs.IsUnstructured == 'yes' :
     ## Value at North Pole for DYNAMICO grid 
@@ -119 +165 @@
 
 ##
+COTAUXXU.attrs['long_name'] = 'Wind stress in geocentric referential - x-component'
+COTAUYYU.attrs['long_name'] = 'Wind stress in geocentric referential - y-component'
+COTAUZZU.attrs['long_name'] = 'Wind stress in geocentric referential - z-component'
+
+COTAUXXU.attrs['units'] = 'Pa'
+COTAUYYU.attrs['units'] = 'Pa'
+COTAUZZU.attrs['units'] = 'Pa'
+
+COTAUXXU.attrs['coordinates'] = "lat lon"
+COTAUYYU.attrs['coordinates'] = "lat lon"
+COTAUZZU.attrs['coordinates'] = "lat lon"
+
 COTAUXXV = COTAUXXU ; COTAUYYV = COTAUYYU ; COTAUZZV = COTAUZZU
-##
+
 
 ## Creates final Dataset
-
-d_Out = xr.Dataset ( {
+lon = xr.DataArray (lon, dims=('y', 'x'))
+lon.attrs['name']      = 'longitude'
+lon.attrs['long_name'] = 'Longitude'
+lon.attrs['units']     = 'degrees_east'
+lon.attrs['bounds']    = 'bounds_lon'
+
+lat = xr.DataArray (lat, dims=('y', 'x'))
+lat.attrs['name']      = 'latitude'
+lat.attrs['long_name'] = 'Latitude'
+lat.attrs['units']     = 'degrees_north'
+lat.attrs['bounds']    = 'bounds_lat'
+
+bounds_lon = d_In.bounds_lon.values.astype (np.float64)
+bounds_lat = d_In.bounds_lat.values.astype (np.float64)
+nvertex = bounds_lon.shape[-1]
+
+bounds_lon = xr.DataArray ( np.reshape (bounds_lon, (ny,nx,nvertex)), dims=('y', 'x', 'nvertex'))
+bounds_lat = xr.DataArray ( np.reshape (bounds_lat, (ny,nx,nvertex)), dims=('y', 'x', 'nvertex'))
+
+bounds_lon.attrs['units'] = 'degrees_east'
+bounds_lat.attrs['units'] = 'degrees_north'
+
+d_Out = xr.Dataset ( {'lat':lat, 'lon':lon, 'bounds_lon':bounds_lon, 'bounds_lat':bounds_lat, 
     'COTOTRAI':COTOTRAI,
     'COTOTSNO':COTOTSNO,
@@ -145 +224 @@
     'COTAUZZV':COTAUZZV } )
 
+d_Out.attrs['Associated files'] = myargs.input
+d_Out.attrs['Conventions']     = "CF-1.6"
+d_Out.attrs['source']          = "IPSL Earth system model"
+d_Out.attrs['group']           = "ICMC IPSL Climate Modelling Center"
+d_Out.attrs['Institution']     = "IPSL https.//www.ipsl.fr"
+try    : d_Out.attrs['directory'] = os.getcwd ()
+except : pass
+try    : d_Out.attrs['HOSTNAME'] = platform.node ()
+except : pass
+try    : d_Out.attrs['LOGNAME']  = os.getlogin ()
+except : pass
+try    : d_Out.attrs['Python']   = "Python version " +  platform.python_version ()
+except : pass
+try    : d_Out.attrs['OS']       = platform.system ()
+except : pass
+try    : d_Out.attrs['release']  = platform.release ()
+except : pass
+try    : d_Out.attrs['hardware'] = platform.machine ()
+except : pass
+d_Out.attrs['conventions']     = "SCRIP"
+d_Out.attrs['dest_grid']       = "curvilinear"
+d_Out.attrs['Model']           = "IPSL CM6"
+d_Out.attrs['SVN_Author']      = "$Author: omamce $"
+d_Out.attrs['SVN_Date']        = "$Date: 2022-07-06 11:06:07 +0200 (Wed, 06 Jul 2022) $"
+d_Out.attrs['SVN_Revision']    = "$Revision: 6190 $"
+d_Out.attrs['SVN_Id']          = "$Id: CalvingWeights.py 6190 2022-07-06 09:06:07Z omamce $"
+d_Out.attrs['SVN_HeadURL']     = "$HeadURL: svn+ssh://omamce@forge.ipsl.jussieu.fr/ipsl/forge/projets/igcmg/svn/TOOLS/MOSAIX/CalvingWeights.py $"
+
+
 d_Out.to_netcdf ( myargs.output, format=NcFormat )
-
+d_Out.close()
 ## ===========================================================================
 ##