| 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 | #======================================= |
|---|
