source: TOOLS/MOSAIX/check_conserv.py @ 4153

Last change on this file since 4153 was 4153, checked in by omamce, 3 years ago

O.M. : include areas correction

  • Property svn:keywords set to Date Revision HeadURL Author
File size: 8.5 KB
Line 
1### ===========================================================================
2###
3### Checks integrals
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# python3 -i check_conserv.py --file dia_tORCA2.3_to_tLMD9695_1stOrder_Normalized_Surfacic_NoArea.nc --grids grids_ORCA2.3xLMD9695.nc --areas areas_ORCA2.3xLMD9695.nc --masks masks_ORCA2.3xLMD9695.nc --frac LMD9695_grid_maskFrom_ORCA2.3.nc
18#
19## SVN information
20__Author__   = "$Author$"
21__Date__     = "$Date$"
22__Revision__ = "$Revision$"
23__Id__       = "$Id:  $"
24__HeadURL    = "$HeadURL$"
25
26import netCDF4
27import numpy as np
28import glob
29import sys, argparse, textwrap
30import nemo
31import os
32
33# Creating a parser
34parser = argparse.ArgumentParser (
35    description = """Compute integral of matching coupling fields on atmosphere and ocean side.""",
36    epilog='-------- This is the end of the help message --------')
37
38# Adding arguments
39parser.add_argument ( '--file'      , help='Input file, e.g. dia_tICO40_to_teORCA1.2_1storder_false.nc', type=str, default='dia_tico_to_torc_1storder_false.nc' )
40parser.add_argument ( '--grids'     , help='grids file', type=str, default='grids.nc' )
41parser.add_argument ( '--areas'     , help='areas file', type=str, default='areas.nc' )
42parser.add_argument ( '--masks'     , help='masks file', type=str, default='masks.nc' )
43parser.add_argument ( '--fracs'     , help='fracs file', type=str, default=None )
44parser.add_argument ( '--guess'     , help='If set, build grids, areas and masks file name from the dia file name', action="store_true", default=False )
45parser.add_argument ( '--src_int'   , help='Source field is integrated'      , action="store_true" )
46parser.add_argument ( '--dst_int'   , help='Destination field is integrated' , action="store_true" )
47parser.add_argument ( '--src_nofrac', help='Do not apply frac on source'     , action="store_true" )
48parser.add_argument ( '--dst_nofrac', help='Do not apply frac on destination', action="store_true" )
49
50# Parse command line
51myargs = parser.parse_args()
52
53#
54name_input = myargs.file
55
56# Get dir and file
57InDir  =  os.path.dirname  (os.path.abspath(myargs.file))
58InFile =  os.path.basename (os.path.abspath(myargs.file))
59
60print ('InDir  : ' + InDir )
61print ('InFile : ' + InFile )
62
63Elements = InFile.split('_')
64
65# Get models
66Name_src = Elements[1][1:] ; Name_dst = Elements[3][1:]
67grid_src = Elements[1][0]  ; grid_dst = Elements[3][0]
68
69Grid_src = grid_src.upper() ; Grid_dst = grid_dst.upper()
70if Grid_src in ( 'A', 'C') : Grid_src = 'T'
71if Grid_dst in ( 'A', 'C') : Grid_dst = 'T'
72
73# src side
74if Name_src.count('ico') != 0 or Name_src.count('lmd') != 0 : name_src = Name_src ; ShortNames = True
75else :
76    if Name_src.count('ICO') : name_src = 'ico' ; ShortNames = False
77    if Name_src.count('LMD') : name_src = 'lmd' ; ShortNames = False
78
79if Name_dst.count('orc') != 0    : name_src = Name_src ; ShortNames = True
80else :
81    if Name_src.count('ORC')!= 0 : name_src = 'orc'    ; ShortNames = False
82       
83# dst side
84if Name_dst.count('ico')!= 0 or Name_dst.count('lmd') != 0 : name_dst = Name_dst ; ShortNames = True
85else :
86    if Name_dst.count('ICO') : name_dst = 'ico' ; ShortNames = False
87    if Name_dst.count('LMD') : name_dst = 'lmd' ; ShortNames = False
88
89if Name_dst.count('orc') != 0 : name_dst = Name_dst ; ShortNames = True
90else :
91    if Name_dst.count('ORC') != 0 : name_dst = 'orc' ; ShortNames = False
92   
93if Name_src.count('ORC') : CplModel = Name_src + 'x' + Name_dst
94if Name_dst.count('ORC') : CplModel = Name_dst + 'x' + Name_src
95
96print ('CplModel   : ' + CplModel )
97print ('Input file : ' + name_input )
98     
99# Coordonnees et masques
100if myargs.fracs == None :
101    if Name_src.count('ORC') : name_frc = InDir + '/' + Name_dst + '_grid_maskFrom_' + Name_src + '.nc'
102    if Name_dst.count('ORC') : name_frc = InDir + '/' + Name_src + '_grid_maskFrom_' + Name_dst + '.nc'
103else:
104    name_frc = myargs.fracs
105
106if myargs.guess :
107    n_grids =  InDir + '/grids_' + CplModel + '.nc'
108    n_masks =  InDir + '/masks_' + CplModel + '.nc'
109    n_areas =  InDir + '/areas_' + CplModel + '.nc'
110else :
111    n_grids = myargs.grids
112    n_masks = myargs.masks
113    n_areas = myargs.areas
114
115print ('Opening grids file: ' + n_grids )
116f_grids = netCDF4.Dataset ( n_grids )
117print ('Opening masks file: ' + n_masks )
118f_masks = netCDF4.Dataset  ( n_masks )
119print ('Opening areas file: ' + n_areas )
120f_areas = netCDF4.Dataset  ( n_areas )
121print ('Opening fracs file: ' + name_frc )   
122f_frac  = netCDF4.Dataset  ( name_frc     )
123
124#
125msk_src = np.float64(1.0) - f_masks.variables[ grid_src + name_src + '.msk'][:].squeeze()
126lon_src = f_grids.variables [grid_src + name_src + '.lon'][:].squeeze()
127lat_src = f_grids.variables [grid_src + name_src + '.lat'][:].squeeze()
128if myargs.src_int : srf_src = np.ones ( (msk_src.shape), dtype=np.float64 )[:]
129else              : srf_src = f_areas.variables [grid_src + name_src + '.srf'][:].squeeze()
130
131#
132msk_dst = np.float64(1.0) - f_masks.variables[grid_dst + name_dst + '.msk'][:]
133lon_dst = f_grids.variables [grid_dst + name_dst + '.lon'][:].squeeze()
134lat_dst = f_grids.variables [grid_dst + name_dst + '.lat'][:].squeeze()
135if myargs.dst_int : srf_dst = np.ones ( (msk_dst.shape), dtype=np.float64)[:]
136else              : srf_dst = f_areas.variables [grid_dst + name_dst + '.srf'][:].squeeze()
137
138#
139if name_src in ( 'ico', 'lmd' ) :
140    if myargs.src_nofrac : frc_src = np.ones ( (msk_src.shape), dtype=np.float64)
141    else                 : frc_src = f_frac.variables ['OceFrac'][:].squeeze()
142       
143if name_src == 'orc'     : frc_src = np.ones ( (msk_src.shape), dtype=np.float64)
144   
145if name_dst in ( 'ico', 'lmd' ) :
146    if myargs.dst_nofrac : frc_dst = np.ones ( (msk_dst.shape), dtype=np.float64)
147    else                 : frc_dst = f_frac.variables['OceFrac'][:].squeeze()
148       
149if name_dst == 'orc'     : frc_dst = np.ones ( (msk_dst.shape), dtype=np.float64)
150
151# Periodicity
152nperio_src = 0 ; nperio_dst = 0
153
154if Name_src == 'ORCA2.3'  : nperio_src = 4
155if Name_src == 'eORCA1.2' : nperio_src = 6
156if Name_src == 'eORCA025' : nperio_src = 6
157
158if Name_dst == 'ORCA2.3'  : nperio_dst = 4
159if Name_dst == 'eORCA1.2' : nperio_dst = 6
160if Name_dst == 'eORCA025' : nperio_dst = 6
161
162# NEMO periodicity.
163if Name_src in ('ORCA2.3', 'eORCA1.2', 'eORCA025') : msk_src = nemo.lbc_mask ( msk_src, nperio=nperio_src, cd_type=Grid_src )
164if Name_dst in ('ORCA2.3', 'eORCA1.2', 'eORCA025') : msk_dst = nemo.lbc_mask ( msk_dst, nperio=nperio_dst, cd_type=Grid_dst )
165   
166# Surfaces
167area_src = np.sum ( srf_src * msk_src * frc_src )
168area_dst = np.sum ( srf_dst * msk_dst * frc_dst )
169
170print ( Name_src )
171print ( 'mask: {:12.3} {:12.3} {:12.3}'.format( np.min(msk_src), np.max(msk_src), np.sum(msk_src) ) )
172print ( 'frac: {:12.3} {:12.3} {:12.3}'.format( np.min(frc_src), np.max(frc_src), np.sum(frc_src) ) )
173print ( 'area: {:12.3} {:12.3} {:12.3}'.format( np.min(srf_src), np.max(srf_src), np.sum(srf_src) ) )
174
175print ( Name_dst )
176print ( 'mask: {:12.3} {:12.3} {:12.3}'.format( np.min(msk_dst), np.max(msk_dst), np.sum(msk_dst) ) )
177print ( 'frac: {:12.3} {:12.3} {:12.3}'.format( np.min(frc_dst), np.max(frc_dst), np.sum(frc_dst) ) )
178print ( 'area: {:12.3} {:12.3} {:12.3}'.format( np.min(srf_dst), np.max(srf_dst), np.sum(srf_dst) ) )
179print (' ')
180print ( "Surfaces   : %14.6e %14.6e %11.3e"%( area_src, area_dst, (area_src-area_dst)/(area_src+area_dst)*0.5 ) )
181print (' ')
182
183##
184f_input = netCDF4.Dataset( name_input )
185# Loop over fields
186for num in np.arange (1, 7) :
187    name_src = 'field{:02}_src'.format(num)
188    name_dst = 'field{:02}_dst'.format(num)
189    print ( '{:02}'.format(num) + ':' + name_src + ':' + name_dst)
190
191    # Reading field
192    v_src = f_input [ name_src ][:]
193    v_dst = f_input [ name_dst ][:]
194
195    # Compute integrals
196    sum_src = np.sum (v_src * srf_src * msk_src * frc_src)
197    sum_dst = np.sum (v_dst * srf_dst * msk_dst * frc_dst)
198   
199    # Compute average
200    mean_src = sum_src / area_src
201    mean_dst = sum_dst / area_dst
202
203    #
204    print ( "Integrals : %14.6e %14.6e %11.3e"%( sum_src , sum_dst , (sum_src - sum_dst)/(sum_src + sum_dst)*0.5 ) )
205    print ( "Averages  : %14.6e %14.6e %11.3e"%( mean_src, mean_dst, (mean_src-mean_dst)/(mean_src+mean_dst)*0.5 ) )
206    print ( " " )
207#=======================================
Note: See TracBrowser for help on using the repository browser.