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