[4084] | 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 | ## |
---|
[4153] | 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 | # |
---|
[4086] | 19 | ## SVN information |
---|
| 20 | __Author__ = "$Author$" |
---|
| 21 | __Date__ = "$Date$" |
---|
| 22 | __Revision__ = "$Revision$" |
---|
| 23 | __Id__ = "$Id: $" |
---|
| 24 | __HeadURL = "$HeadURL$" |
---|
| 25 | |
---|
[4153] | 26 | import netCDF4 |
---|
[4084] | 27 | import numpy as np |
---|
| 28 | import glob |
---|
| 29 | import sys, argparse, textwrap |
---|
| 30 | import nemo |
---|
[4153] | 31 | import os |
---|
[4084] | 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 |
---|
[4089] | 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" ) |
---|
[4084] | 49 | |
---|
| 50 | # Parse command line |
---|
| 51 | myargs = parser.parse_args() |
---|
| 52 | |
---|
| 53 | # |
---|
| 54 | name_input = myargs.file |
---|
| 55 | |
---|
[4153] | 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)) |
---|
[4084] | 59 | |
---|
[4153] | 60 | print ('InDir : ' + InDir ) |
---|
| 61 | print ('InFile : ' + InFile ) |
---|
| 62 | |
---|
| 63 | Elements = InFile.split('_') |
---|
| 64 | |
---|
| 65 | # Get models |
---|
[4084] | 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 | |
---|
[4153] | 96 | print ('CplModel : ' + CplModel ) |
---|
| 97 | print ('Input file : ' + name_input ) |
---|
[4084] | 98 | |
---|
| 99 | # Coordonnees et masques |
---|
| 100 | if myargs.fracs == None : |
---|
[4153] | 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' |
---|
[4084] | 103 | else: |
---|
| 104 | name_frc = myargs.fracs |
---|
| 105 | |
---|
| 106 | if myargs.guess : |
---|
[4153] | 107 | n_grids = InDir + '/grids_' + CplModel + '.nc' |
---|
| 108 | n_masks = InDir + '/masks_' + CplModel + '.nc' |
---|
| 109 | n_areas = InDir + '/areas_' + CplModel + '.nc' |
---|
[4084] | 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 ) |
---|
[4153] | 116 | f_grids = netCDF4.Dataset ( n_grids ) |
---|
[4084] | 117 | print ('Opening masks file: ' + n_masks ) |
---|
[4153] | 118 | f_masks = netCDF4.Dataset ( n_masks ) |
---|
[4084] | 119 | print ('Opening areas file: ' + n_areas ) |
---|
[4153] | 120 | f_areas = netCDF4.Dataset ( n_areas ) |
---|
[4084] | 121 | print ('Opening fracs file: ' + name_frc ) |
---|
[4153] | 122 | f_frac = netCDF4.Dataset ( name_frc ) |
---|
[4084] | 123 | |
---|
| 124 | # |
---|
[4153] | 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() |
---|
[4084] | 130 | |
---|
| 131 | # |
---|
[4153] | 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() |
---|
[4084] | 137 | |
---|
| 138 | # |
---|
| 139 | if name_src in ( 'ico', 'lmd' ) : |
---|
| 140 | if myargs.src_nofrac : frc_src = np.ones ( (msk_src.shape), dtype=np.float64) |
---|
[4153] | 141 | else : frc_src = f_frac.variables ['OceFrac'][:].squeeze() |
---|
| 142 | |
---|
[4084] | 143 | if name_src == 'orc' : frc_src = np.ones ( (msk_src.shape), dtype=np.float64) |
---|
[4153] | 144 | |
---|
[4084] | 145 | if name_dst in ( 'ico', 'lmd' ) : |
---|
| 146 | if myargs.dst_nofrac : frc_dst = np.ones ( (msk_dst.shape), dtype=np.float64) |
---|
[4153] | 147 | else : frc_dst = f_frac.variables['OceFrac'][:].squeeze() |
---|
| 148 | |
---|
[4084] | 149 | if name_dst == 'orc' : frc_dst = np.ones ( (msk_dst.shape), dtype=np.float64) |
---|
| 150 | |
---|
| 151 | # Periodicity |
---|
[4089] | 152 | nperio_src = 0 ; nperio_dst = 0 |
---|
[4084] | 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 ) |
---|
[4089] | 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 ) |
---|
[4084] | 169 | |
---|
| 170 | print ( Name_src ) |
---|
[4153] | 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) ) ) |
---|
[4084] | 174 | |
---|
| 175 | print ( Name_dst ) |
---|
[4153] | 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) ) ) |
---|
[4084] | 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 | ## |
---|
[4153] | 184 | f_input = netCDF4.Dataset( name_input ) |
---|
[4084] | 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 |
---|
[4153] | 192 | v_src = f_input [ name_src ][:] |
---|
| 193 | v_dst = f_input [ name_dst ][:] |
---|
[4084] | 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 | #======================================= |
---|