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 | |
---|
26 | import netCDF4 |
---|
27 | import numpy as np |
---|
28 | import glob |
---|
29 | import sys, argparse, textwrap |
---|
30 | import nemo |
---|
31 | import os |
---|
32 | |
---|
33 | # Creating a parser |
---|
34 | parser = 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 |
---|
39 | parser.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' ) |
---|
40 | parser.add_argument ( '--grids' , help='grids file', type=str, default='grids.nc' ) |
---|
41 | parser.add_argument ( '--areas' , help='areas file', type=str, default='areas.nc' ) |
---|
42 | parser.add_argument ( '--masks' , help='masks file', type=str, default='masks.nc' ) |
---|
43 | parser.add_argument ( '--fracs' , help='fracs file', type=str, default=None ) |
---|
44 | parser.add_argument ( '--guess' , help='If set, build grids, areas and masks file name from the dia file name', action="store_true", default=False ) |
---|
45 | parser.add_argument ( '--src_int' , help='Source field is integrated' , action="store_true" ) |
---|
46 | parser.add_argument ( '--dst_int' , help='Destination field is integrated' , action="store_true" ) |
---|
47 | parser.add_argument ( '--src_nofrac', help='Do not apply frac on source' , action="store_true" ) |
---|
48 | parser.add_argument ( '--dst_nofrac', help='Do not apply frac on destination', action="store_true" ) |
---|
49 | |
---|
50 | # Parse command line |
---|
51 | myargs = parser.parse_args() |
---|
52 | |
---|
53 | # |
---|
54 | name_input = myargs.file |
---|
55 | |
---|
56 | # Get dir and file |
---|
57 | InDir = os.path.dirname (os.path.abspath(myargs.file)) |
---|
58 | InFile = os.path.basename (os.path.abspath(myargs.file)) |
---|
59 | |
---|
60 | print ('InDir : ' + InDir ) |
---|
61 | print ('InFile : ' + InFile ) |
---|
62 | |
---|
63 | Elements = InFile.split('_') |
---|
64 | |
---|
65 | # Get models |
---|
66 | Name_src = Elements[1][1:] ; Name_dst = Elements[3][1:] |
---|
67 | grid_src = Elements[1][0] ; grid_dst = Elements[3][0] |
---|
68 | |
---|
69 | Grid_src = grid_src.upper() ; Grid_dst = grid_dst.upper() |
---|
70 | if Grid_src in ( 'A', 'C') : Grid_src = 'T' |
---|
71 | if Grid_dst in ( 'A', 'C') : Grid_dst = 'T' |
---|
72 | |
---|
73 | # src side |
---|
74 | if Name_src.count('ico') != 0 or Name_src.count('lmd') != 0 : name_src = Name_src ; ShortNames = True |
---|
75 | else : |
---|
76 | if Name_src.count('ICO') : name_src = 'ico' ; ShortNames = False |
---|
77 | if Name_src.count('LMD') : name_src = 'lmd' ; ShortNames = False |
---|
78 | |
---|
79 | if Name_dst.count('orc') != 0 : name_src = Name_src ; ShortNames = True |
---|
80 | else : |
---|
81 | if Name_src.count('ORC')!= 0 : name_src = 'orc' ; ShortNames = False |
---|
82 | |
---|
83 | # dst side |
---|
84 | if Name_dst.count('ico')!= 0 or Name_dst.count('lmd') != 0 : name_dst = Name_dst ; ShortNames = True |
---|
85 | else : |
---|
86 | if Name_dst.count('ICO') : name_dst = 'ico' ; ShortNames = False |
---|
87 | if Name_dst.count('LMD') : name_dst = 'lmd' ; ShortNames = False |
---|
88 | |
---|
89 | if Name_dst.count('orc') != 0 : name_dst = Name_dst ; ShortNames = True |
---|
90 | else : |
---|
91 | if Name_dst.count('ORC') != 0 : name_dst = 'orc' ; ShortNames = False |
---|
92 | |
---|
93 | if Name_src.count('ORC') : CplModel = Name_src + 'x' + Name_dst |
---|
94 | if Name_dst.count('ORC') : CplModel = Name_dst + 'x' + Name_src |
---|
95 | |
---|
96 | print ('CplModel : ' + CplModel ) |
---|
97 | print ('Input file : ' + name_input ) |
---|
98 | |
---|
99 | # Coordonnees et masques |
---|
100 | if 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' |
---|
103 | else: |
---|
104 | name_frc = myargs.fracs |
---|
105 | |
---|
106 | if myargs.guess : |
---|
107 | n_grids = InDir + '/grids_' + CplModel + '.nc' |
---|
108 | n_masks = InDir + '/masks_' + CplModel + '.nc' |
---|
109 | n_areas = InDir + '/areas_' + CplModel + '.nc' |
---|
110 | else : |
---|
111 | n_grids = myargs.grids |
---|
112 | n_masks = myargs.masks |
---|
113 | n_areas = myargs.areas |
---|
114 | |
---|
115 | print ('Opening grids file: ' + n_grids ) |
---|
116 | f_grids = netCDF4.Dataset ( n_grids ) |
---|
117 | print ('Opening masks file: ' + n_masks ) |
---|
118 | f_masks = netCDF4.Dataset ( n_masks ) |
---|
119 | print ('Opening areas file: ' + n_areas ) |
---|
120 | f_areas = netCDF4.Dataset ( n_areas ) |
---|
121 | print ('Opening fracs file: ' + name_frc ) |
---|
122 | f_frac = netCDF4.Dataset ( name_frc ) |
---|
123 | |
---|
124 | # |
---|
125 | msk_src = np.float64(1.0) - f_masks.variables[ grid_src + name_src + '.msk'][:].squeeze() |
---|
126 | lon_src = f_grids.variables [grid_src + name_src + '.lon'][:].squeeze() |
---|
127 | lat_src = f_grids.variables [grid_src + name_src + '.lat'][:].squeeze() |
---|
128 | if myargs.src_int : srf_src = np.ones ( (msk_src.shape), dtype=np.float64 )[:] |
---|
129 | else : srf_src = f_areas.variables [grid_src + name_src + '.srf'][:].squeeze() |
---|
130 | |
---|
131 | # |
---|
132 | msk_dst = np.float64(1.0) - f_masks.variables[grid_dst + name_dst + '.msk'][:] |
---|
133 | lon_dst = f_grids.variables [grid_dst + name_dst + '.lon'][:].squeeze() |
---|
134 | lat_dst = f_grids.variables [grid_dst + name_dst + '.lat'][:].squeeze() |
---|
135 | if myargs.dst_int : srf_dst = np.ones ( (msk_dst.shape), dtype=np.float64)[:] |
---|
136 | else : srf_dst = f_areas.variables [grid_dst + name_dst + '.srf'][:].squeeze() |
---|
137 | |
---|
138 | # |
---|
139 | if 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 | |
---|
143 | if name_src == 'orc' : frc_src = np.ones ( (msk_src.shape), dtype=np.float64) |
---|
144 | |
---|
145 | if 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 | |
---|
149 | if name_dst == 'orc' : frc_dst = np.ones ( (msk_dst.shape), dtype=np.float64) |
---|
150 | |
---|
151 | # Periodicity |
---|
152 | nperio_src = 0 ; nperio_dst = 0 |
---|
153 | |
---|
154 | if Name_src == 'ORCA2.3' : nperio_src = 4 |
---|
155 | if Name_src == 'eORCA1.2' : nperio_src = 6 |
---|
156 | if Name_src == 'eORCA025' : nperio_src = 6 |
---|
157 | |
---|
158 | if Name_dst == 'ORCA2.3' : nperio_dst = 4 |
---|
159 | if Name_dst == 'eORCA1.2' : nperio_dst = 6 |
---|
160 | if Name_dst == 'eORCA025' : nperio_dst = 6 |
---|
161 | |
---|
162 | # NEMO periodicity. |
---|
163 | if Name_src in ('ORCA2.3', 'eORCA1.2', 'eORCA025') : msk_src = nemo.lbc_mask ( msk_src, nperio=nperio_src, cd_type=Grid_src ) |
---|
164 | if 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 |
---|
167 | area_src = np.sum ( srf_src * msk_src * frc_src ) |
---|
168 | area_dst = np.sum ( srf_dst * msk_dst * frc_dst ) |
---|
169 | |
---|
170 | print ( Name_src ) |
---|
171 | print ( 'mask: {:12.3} {:12.3} {:12.3}'.format( np.min(msk_src), np.max(msk_src), np.sum(msk_src) ) ) |
---|
172 | print ( 'frac: {:12.3} {:12.3} {:12.3}'.format( np.min(frc_src), np.max(frc_src), np.sum(frc_src) ) ) |
---|
173 | print ( 'area: {:12.3} {:12.3} {:12.3}'.format( np.min(srf_src), np.max(srf_src), np.sum(srf_src) ) ) |
---|
174 | |
---|
175 | print ( Name_dst ) |
---|
176 | print ( 'mask: {:12.3} {:12.3} {:12.3}'.format( np.min(msk_dst), np.max(msk_dst), np.sum(msk_dst) ) ) |
---|
177 | print ( 'frac: {:12.3} {:12.3} {:12.3}'.format( np.min(frc_dst), np.max(frc_dst), np.sum(frc_dst) ) ) |
---|
178 | print ( 'area: {:12.3} {:12.3} {:12.3}'.format( np.min(srf_dst), np.max(srf_dst), np.sum(srf_dst) ) ) |
---|
179 | print (' ') |
---|
180 | print ( "Surfaces : %14.6e %14.6e %11.3e"%( area_src, area_dst, (area_src-area_dst)/(area_src+area_dst)*0.5 ) ) |
---|
181 | print (' ') |
---|
182 | |
---|
183 | ## |
---|
184 | f_input = netCDF4.Dataset( name_input ) |
---|
185 | # Loop over fields |
---|
186 | for 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 | #======================================= |
---|