1 | ### =========================================================================== |
---|
2 | ### |
---|
3 | ### Checks integrals |
---|
4 | ### |
---|
5 | ### =========================================================================== |
---|
6 | ## |
---|
7 | ## MOSAIX is under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt" |
---|
8 | ## file for an english version of the licence and |
---|
9 | ## "Licence_CeCILL_V2-fr.txt" for a french version. |
---|
10 | ## |
---|
11 | ## Permission is hereby granted, free of charge, to any person or |
---|
12 | ## organization obtaining a copy of the software and accompanying |
---|
13 | ## documentation covered by this license (the "Software") to use, |
---|
14 | ## reproduce, display, distribute, execute, and transmit the |
---|
15 | ## Software, and to prepare derivative works of the Software, and to |
---|
16 | ## permit third-parties to whom the Software is furnished to do so, |
---|
17 | ## all subject to the following: |
---|
18 | ## |
---|
19 | ## Warning, to install, configure, run, use any of MOSAIX software or |
---|
20 | ## to read the associated documentation you'll need at least one (1) |
---|
21 | ## brain in a reasonably working order. Lack of this implement will |
---|
22 | ## void any warranties (either express or implied). Authors assumes |
---|
23 | ## no responsability for errors, omissions, data loss, or any other |
---|
24 | ## consequences caused directly or indirectly by the usage of his |
---|
25 | ## software by incorrectly or partially configured |
---|
26 | ## |
---|
27 | # |
---|
28 | # 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 |
---|
29 | # |
---|
30 | ## SVN information |
---|
31 | __Author__ = "$Author$" |
---|
32 | __Date__ = "$Date$" |
---|
33 | __Revision__ = "$Revision$" |
---|
34 | __Id__ = "$Id: $" |
---|
35 | __HeadURL = "$HeadURL$" |
---|
36 | |
---|
37 | import netCDF4 |
---|
38 | import numpy as np |
---|
39 | import glob |
---|
40 | import sys, argparse, textwrap |
---|
41 | import nemo |
---|
42 | import os |
---|
43 | |
---|
44 | # Creating a parser |
---|
45 | parser = argparse.ArgumentParser ( |
---|
46 | description = """Compute integral of matching coupling fields on atmosphere and ocean side.""", |
---|
47 | epilog='-------- This is the end of the help message --------') |
---|
48 | |
---|
49 | # Adding arguments |
---|
50 | 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' ) |
---|
51 | parser.add_argument ( '--grids' , help='grids file', type=str, default='grids.nc' ) |
---|
52 | parser.add_argument ( '--areas' , help='areas file', type=str, default='areas.nc' ) |
---|
53 | parser.add_argument ( '--masks' , help='masks file', type=str, default='masks.nc' ) |
---|
54 | parser.add_argument ( '--fracs' , help='fracs file', type=str, default=None ) |
---|
55 | parser.add_argument ( '--guess' , help='If set, build grids, areas and masks file name from the dia file name', action="store_true", default=False ) |
---|
56 | parser.add_argument ( '--src_int' , help='Source field is integrated' , action="store_true" ) |
---|
57 | parser.add_argument ( '--dst_int' , help='Destination field is integrated' , action="store_true" ) |
---|
58 | parser.add_argument ( '--src_nofrac', help='Do not apply frac on source' , action="store_true" ) |
---|
59 | parser.add_argument ( '--dst_nofrac', help='Do not apply frac on destination', action="store_true" ) |
---|
60 | |
---|
61 | # Parse command line |
---|
62 | myargs = parser.parse_args() |
---|
63 | |
---|
64 | # |
---|
65 | name_input = myargs.file |
---|
66 | |
---|
67 | # Get dir and file |
---|
68 | InDir = os.path.dirname (os.path.abspath(myargs.file)) |
---|
69 | InFile = os.path.basename (os.path.abspath(myargs.file)) |
---|
70 | |
---|
71 | print ('InDir : ' + InDir ) |
---|
72 | print ('InFile : ' + InFile ) |
---|
73 | |
---|
74 | Elements = InFile.split('_') |
---|
75 | |
---|
76 | # Get models |
---|
77 | Name_src = Elements[1][1:] ; Name_dst = Elements[3][1:] |
---|
78 | grid_src = Elements[1][0] ; grid_dst = Elements[3][0] |
---|
79 | |
---|
80 | Grid_src = grid_src.upper() ; Grid_dst = grid_dst.upper() |
---|
81 | if Grid_src in ( 'A', 'C') : Grid_src = 'T' |
---|
82 | if Grid_dst in ( 'A', 'C') : Grid_dst = 'T' |
---|
83 | |
---|
84 | # src side |
---|
85 | if Name_src.count('ico') != 0 or Name_src.count('lmd') != 0 : name_src = Name_src ; ShortNames = True |
---|
86 | else : |
---|
87 | if Name_src.count('ICO') : name_src = 'ico' ; ShortNames = False |
---|
88 | if Name_src.count('LMD') : name_src = 'lmd' ; ShortNames = False |
---|
89 | |
---|
90 | if Name_dst.count('orc') != 0 : name_src = Name_src ; ShortNames = True |
---|
91 | else : |
---|
92 | if Name_src.count('ORC')!= 0 : name_src = 'orc' ; ShortNames = False |
---|
93 | |
---|
94 | # dst side |
---|
95 | if Name_dst.count('ico')!= 0 or Name_dst.count('lmd') != 0 : name_dst = Name_dst ; ShortNames = True |
---|
96 | else : |
---|
97 | if Name_dst.count('ICO') : name_dst = 'ico' ; ShortNames = False |
---|
98 | if Name_dst.count('LMD') : name_dst = 'lmd' ; ShortNames = False |
---|
99 | |
---|
100 | if Name_dst.count('orc') != 0 : name_dst = Name_dst ; ShortNames = True |
---|
101 | else : |
---|
102 | if Name_dst.count('ORC') != 0 : name_dst = 'orc' ; ShortNames = False |
---|
103 | |
---|
104 | if Name_src.count('ORC') : CplModel = Name_src + 'x' + Name_dst |
---|
105 | if Name_dst.count('ORC') : CplModel = Name_dst + 'x' + Name_src |
---|
106 | |
---|
107 | print ('CplModel : ' + CplModel ) |
---|
108 | print ('Input file : ' + name_input ) |
---|
109 | |
---|
110 | # Coordonnees et masques |
---|
111 | if myargs.fracs == None : |
---|
112 | if Name_src.count('ORC') : name_frc = InDir + '/' + Name_dst + '_grid_maskFrom_' + Name_src + '.nc' |
---|
113 | if Name_dst.count('ORC') : name_frc = InDir + '/' + Name_src + '_grid_maskFrom_' + Name_dst + '.nc' |
---|
114 | else: |
---|
115 | name_frc = myargs.fracs |
---|
116 | |
---|
117 | if myargs.guess : |
---|
118 | n_grids = InDir + '/grids_' + CplModel + '.nc' |
---|
119 | n_masks = InDir + '/masks_' + CplModel + '.nc' |
---|
120 | n_areas = InDir + '/areas_' + CplModel + '.nc' |
---|
121 | else : |
---|
122 | n_grids = myargs.grids |
---|
123 | n_masks = myargs.masks |
---|
124 | n_areas = myargs.areas |
---|
125 | |
---|
126 | print ('Opening grids file: ' + n_grids ) |
---|
127 | f_grids = netCDF4.Dataset ( n_grids ) |
---|
128 | print ('Opening masks file: ' + n_masks ) |
---|
129 | f_masks = netCDF4.Dataset ( n_masks ) |
---|
130 | print ('Opening areas file: ' + n_areas ) |
---|
131 | f_areas = netCDF4.Dataset ( n_areas ) |
---|
132 | print ('Opening fracs file: ' + name_frc ) |
---|
133 | f_frac = netCDF4.Dataset ( name_frc ) |
---|
134 | |
---|
135 | # |
---|
136 | msk_src = np.float64(1.0) - f_masks.variables[ grid_src + name_src + '.msk'][:].squeeze() |
---|
137 | lon_src = f_grids.variables [grid_src + name_src + '.lon'][:].squeeze() |
---|
138 | lat_src = f_grids.variables [grid_src + name_src + '.lat'][:].squeeze() |
---|
139 | if myargs.src_int : srf_src = np.ones ( (msk_src.shape), dtype=np.float64 )[:] |
---|
140 | else : srf_src = f_areas.variables [grid_src + name_src + '.srf'][:].squeeze() |
---|
141 | |
---|
142 | # |
---|
143 | msk_dst = np.float64(1.0) - f_masks.variables[grid_dst + name_dst + '.msk'][:] |
---|
144 | lon_dst = f_grids.variables [grid_dst + name_dst + '.lon'][:].squeeze() |
---|
145 | lat_dst = f_grids.variables [grid_dst + name_dst + '.lat'][:].squeeze() |
---|
146 | if myargs.dst_int : srf_dst = np.ones ( (msk_dst.shape), dtype=np.float64)[:] |
---|
147 | else : srf_dst = f_areas.variables [grid_dst + name_dst + '.srf'][:].squeeze() |
---|
148 | |
---|
149 | # |
---|
150 | if name_src in ( 'ico', 'lmd' ) : |
---|
151 | if myargs.src_nofrac : frc_src = np.ones ( (msk_src.shape), dtype=np.float64) |
---|
152 | else : frc_src = f_frac.variables ['OceFrac'][:].squeeze() |
---|
153 | |
---|
154 | if name_src == 'orc' : frc_src = np.ones ( (msk_src.shape), dtype=np.float64) |
---|
155 | |
---|
156 | if name_dst in ( 'ico', 'lmd' ) : |
---|
157 | if myargs.dst_nofrac : frc_dst = np.ones ( (msk_dst.shape), dtype=np.float64) |
---|
158 | else : frc_dst = f_frac.variables['OceFrac'][:].squeeze() |
---|
159 | |
---|
160 | if name_dst == 'orc' : frc_dst = np.ones ( (msk_dst.shape), dtype=np.float64) |
---|
161 | |
---|
162 | # Periodicity |
---|
163 | nperio_src = 0 ; nperio_dst = 0 |
---|
164 | |
---|
165 | if Name_src == 'ORCA2.3' : nperio_src = 4 |
---|
166 | if Name_src == 'eORCA1.2' : nperio_src = 6 |
---|
167 | if Name_src == 'eORCA025' : nperio_src = 6 |
---|
168 | |
---|
169 | if Name_dst == 'ORCA2.3' : nperio_dst = 4 |
---|
170 | if Name_dst == 'eORCA1.2' : nperio_dst = 6 |
---|
171 | if Name_dst == 'eORCA025' : nperio_dst = 6 |
---|
172 | |
---|
173 | # NEMO periodicity. |
---|
174 | if Name_src in ('ORCA2.3', 'eORCA1.2', 'eORCA025') : msk_src = nemo.lbc_mask ( msk_src, nperio=nperio_src, cd_type=Grid_src ) |
---|
175 | if Name_dst in ('ORCA2.3', 'eORCA1.2', 'eORCA025') : msk_dst = nemo.lbc_mask ( msk_dst, nperio=nperio_dst, cd_type=Grid_dst ) |
---|
176 | |
---|
177 | # Surfaces |
---|
178 | area_src = np.sum ( srf_src * msk_src * frc_src ) |
---|
179 | area_dst = np.sum ( srf_dst * msk_dst * frc_dst ) |
---|
180 | |
---|
181 | print ( Name_src ) |
---|
182 | print ( 'mask: {:12.3} {:12.3} {:12.3}'.format( np.min(msk_src), np.max(msk_src), np.sum(msk_src) ) ) |
---|
183 | print ( 'frac: {:12.3} {:12.3} {:12.3}'.format( np.min(frc_src), np.max(frc_src), np.sum(frc_src) ) ) |
---|
184 | print ( 'area: {:12.3} {:12.3} {:12.3}'.format( np.min(srf_src), np.max(srf_src), np.sum(srf_src) ) ) |
---|
185 | |
---|
186 | print ( Name_dst ) |
---|
187 | print ( 'mask: {:12.3} {:12.3} {:12.3}'.format( np.min(msk_dst), np.max(msk_dst), np.sum(msk_dst) ) ) |
---|
188 | print ( 'frac: {:12.3} {:12.3} {:12.3}'.format( np.min(frc_dst), np.max(frc_dst), np.sum(frc_dst) ) ) |
---|
189 | print ( 'area: {:12.3} {:12.3} {:12.3}'.format( np.min(srf_dst), np.max(srf_dst), np.sum(srf_dst) ) ) |
---|
190 | print (' ') |
---|
191 | print ( "Surfaces : %14.6e %14.6e %11.3e"%( area_src, area_dst, (area_src-area_dst)/(area_src+area_dst)*0.5 ) ) |
---|
192 | print (' ') |
---|
193 | |
---|
194 | ## |
---|
195 | f_input = netCDF4.Dataset( name_input ) |
---|
196 | # Loop over fields |
---|
197 | for num in np.arange (1, 7) : |
---|
198 | name_src = 'field{:02}_src'.format(num) |
---|
199 | name_dst = 'field{:02}_dst'.format(num) |
---|
200 | print ( '{:02}'.format(num) + ':' + name_src + ':' + name_dst) |
---|
201 | |
---|
202 | # Reading field |
---|
203 | v_src = f_input [ name_src ][:] |
---|
204 | v_dst = f_input [ name_dst ][:] |
---|
205 | |
---|
206 | # Compute integrals |
---|
207 | sum_src = np.sum (v_src * srf_src * msk_src * frc_src) |
---|
208 | sum_dst = np.sum (v_dst * srf_dst * msk_dst * frc_dst) |
---|
209 | |
---|
210 | # Compute average |
---|
211 | mean_src = sum_src / area_src |
---|
212 | mean_dst = sum_dst / area_dst |
---|
213 | |
---|
214 | # |
---|
215 | print ( "Integrals : %14.6e %14.6e %11.3e"%( sum_src , sum_dst , (sum_src - sum_dst)/(sum_src + sum_dst)*0.5 ) ) |
---|
216 | print ( "Averages : %14.6e %14.6e %11.3e"%( mean_src, mean_dst, (mean_src-mean_dst)/(mean_src+mean_dst)*0.5 ) ) |
---|
217 | print ( " " ) |
---|
218 | #======================================= |
---|