[3325] | 1 | # -*- Mode: python -*- |
---|
| 2 | ### - |
---|
| 3 | ### - Correct rmp file computed by OASIS MCT |
---|
| 4 | ## - |
---|
| 5 | ## - |
---|
| 6 | ## - Olivier Marti (olivier.marti@lsce.ipsl.fr) |
---|
| 7 | ## - Institut Pierre Simon Laplace |
---|
| 8 | ## - Laboratoire des Sciences du Climat et de l'Environnment |
---|
| 9 | ## - |
---|
| 10 | # |
---|
| 11 | # This program overcomes an interpolation bug in Scrip, uses by OASIS-MCT |
---|
| 12 | # |
---|
| 13 | # In IPSLCM6 at resolution ORCA1 x LMD 144x144, the bilinear interpolation of the wind stress is bugged. |
---|
| 14 | # Two NEMO points around 0N/180W takes there wind stress values from points at 0N/0E. |
---|
| 15 | # |
---|
| 16 | # The bilinear interpolation compute weights and adresses that are written in a file : |
---|
| 17 | # rmp_tlmd_to_torc_BILINEAR_ComputedByOasis.nc. |
---|
| 18 | # The program copies the file to rmp_tlmd_to_torc_BILINEAR_Corrected.nc and makes a |
---|
| 19 | # few corrections of adresses and weigths. |
---|
| 20 | # |
---|
| 21 | ## See : https://forge.ipsl.jussieu.fr/igcmg/wiki/IPSLCM6/Coupling/BugBilinear |
---|
| 22 | # |
---|
| 23 | |
---|
| 24 | ## - Needed Python modules |
---|
| 25 | import cdms2, MV2, numpy as np |
---|
| 26 | import os, time, shutil |
---|
| 27 | from cdms2.coord import TransientVirtualAxis |
---|
| 28 | |
---|
| 29 | ## - NetCDF parameters |
---|
| 30 | # - If you want to turn that off or set different values of compression use the functions: |
---|
| 31 | cdms2.setNetcdfShuffleFlag (0) |
---|
| 32 | cdms2.setNetcdfDeflateFlag (0) |
---|
| 33 | cdms2.setNetcdfDeflateLevelFlag(0) |
---|
| 34 | |
---|
| 35 | # |
---|
| 36 | rad = np.rad2deg (1.0) |
---|
| 37 | |
---|
| 38 | ## Copy initial rmp file with bilinear weights |
---|
| 39 | shutil.copyfile ( 'rmp_tlmd_to_torc_BILINEAR_ComputedByOasis.nc' , 'rmp_tlmd_to_torc_BILINEAR_Corrected.nc' ) |
---|
| 40 | |
---|
| 41 | ## Read the new file |
---|
| 42 | rmp_file = cdms2.open ( 'rmp_tlmd_to_torc_BILINEAR_Corrected.nc' , 'r+') |
---|
| 43 | |
---|
| 44 | src_grid_center_lat = rmp_file ( 'src_grid_center_lat' ) |
---|
| 45 | src_grid_center_lon = rmp_file ( 'src_grid_center_lon' ) |
---|
| 46 | dst_grid_center_lat = rmp_file ( 'dst_grid_center_lat' ) |
---|
| 47 | dst_grid_center_lon = rmp_file ( 'dst_grid_center_lon' ) |
---|
| 48 | src_address = rmp_file ( 'src_address' ) |
---|
| 49 | dst_address = rmp_file ( 'dst_address' ) |
---|
| 50 | remap_matrix = rmp_file ( 'remap_matrix' ) |
---|
| 51 | src_grid_dims = rmp_file ( 'src_grid_dims' ) |
---|
| 52 | dst_grid_dims = rmp_file ( 'dst_grid_dims' ) |
---|
| 53 | |
---|
| 54 | src_grid_corner_lat = rmp_file ( 'src_grid_corner_lat' ) |
---|
| 55 | src_grid_corner_lon = rmp_file ( 'src_grid_corner_lon' ) |
---|
| 56 | src_grid_area = rmp_file ( 'src_grid_area' ) |
---|
| 57 | src_grid_imask = rmp_file ( 'src_grid_imask' ) |
---|
| 58 | |
---|
| 59 | dst_grid_corner_lat = rmp_file ( 'dst_grid_corner_lat' ) |
---|
| 60 | dst_grid_corner_lon = rmp_file ( 'dst_grid_corner_lon' ) |
---|
| 61 | dst_grid_area = rmp_file ( 'dst_grid_area' ) |
---|
| 62 | dst_grid_imask = rmp_file ( 'dst_grid_imask' ) |
---|
| 63 | |
---|
| 64 | # Get dimensions |
---|
| 65 | ( src_grid_size ) = src_grid_center_lat.shape |
---|
| 66 | ( dst_grid_size ) = dst_grid_center_lat.shape |
---|
| 67 | ( num_links, num_wgts) = remap_matrix.shape |
---|
| 68 | ( src_grid_rank ) = src_grid_dims.shape |
---|
| 69 | ( dst_grid_rank ) = dst_grid_dims.shape |
---|
| 70 | |
---|
| 71 | print 'Dimensions ' |
---|
| 72 | print 'src_grid_size : ', src_grid_size, ' src_grid_dims : ', src_grid_dims |
---|
| 73 | print 'dst_grid_size : ', dst_grid_size, ' dst_grid_dims : ', dst_grid_dims |
---|
| 74 | print ' ' |
---|
| 75 | |
---|
| 76 | ## Get list of wrong weights |
---|
| 77 | list_ind = np.where ( np.abs(remap_matrix) > 1.0 )[0] |
---|
| 78 | print 'List of indexes of buggy points : ', list_ind |
---|
| 79 | print ' ' |
---|
| 80 | |
---|
| 81 | for ind in list_ind[:]: |
---|
| 82 | print 'Working on point : ', ind, remap_matrix[ind] |
---|
| 83 | ## Compute 2D indexes |
---|
| 84 | ## warning : ind is a 'C' index |
---|
| 85 | src_grid_ind_2 = (src_address[ind]+1)/src_grid_dims[0] + 1 |
---|
| 86 | src_grid_ind_1 = (src_address[ind]+1) - ( src_grid_ind_2 - 1 ) * src_grid_dims[0] |
---|
| 87 | |
---|
| 88 | dst_grid_ind_2 = dst_address[ind]/dst_grid_dims[0] + 1 |
---|
| 89 | dst_grid_ind_1 = dst_address[ind] - ( dst_grid_ind_2 - 1 ) * dst_grid_dims[0] |
---|
| 90 | |
---|
| 91 | print src_address[ind], src_grid_ind_1, src_grid_ind_2, src_grid_center_lon[src_address[ind]]*rad, src_grid_center_lat[src_address[ind]]*rad |
---|
| 92 | print dst_address[ind], dst_grid_ind_1, dst_grid_ind_2, dst_grid_center_lon[dst_address[ind]]*rad, dst_grid_center_lat[dst_address[ind]]*rad |
---|
| 93 | print dst_grid_center_lon[dst_address[ind]]*rad - src_grid_center_lon[src_address[ind]]*rad, |
---|
| 94 | print dst_grid_center_lat[dst_address[ind]]*rad - src_grid_center_lat[src_address[ind]]*rad |
---|
| 95 | |
---|
| 96 | # |
---|
| 97 | new_src_grid_ind_1 = np.mod ( (src_grid_ind_1 - src_grid_dims[0]/2)-1, src_grid_dims[0]) + 1 # Decalage de 180 degres ... |
---|
| 98 | new_src_grid_ind = ( src_grid_ind_2 - 1 ) * src_grid_dims[0] + new_src_grid_ind_1 - 1 |
---|
| 99 | |
---|
| 100 | print src_address[ind], new_src_grid_ind, new_src_grid_ind_1, src_grid_ind_2 |
---|
| 101 | src_address [ind] = new_src_grid_ind |
---|
| 102 | |
---|
| 103 | dist = dst_grid_center_lon[dst_address[ind]-1]*rad - src_grid_center_lon[src_address[ind]-1]*rad |
---|
| 104 | |
---|
| 105 | if dist == 0.0: |
---|
| 106 | remap_matrix [ind] = 1.0 |
---|
| 107 | other_src_grid_ind = 0 ; delta_x = 0.0 |
---|
| 108 | else: |
---|
| 109 | if dist > 0.0: |
---|
| 110 | other_src_grid_ind_1 = np.mod ( (new_src_grid_ind_1 - 1 )-1, src_grid_dims[0]) + 1 |
---|
| 111 | if dist < 0.0: |
---|
| 112 | other_src_grid_ind_1 = np.mod ( (new_src_grid_ind_1 + 1 )-1, src_grid_dims[0]) + 1 |
---|
| 113 | |
---|
| 114 | other_src_grid_ind = ( src_grid_ind_2 - 1 ) * src_grid_dims[0] + other_src_grid_ind_1 - 1 |
---|
| 115 | delta_x = src_grid_center_lon[new_src_grid_ind]*rad- src_grid_center_lon[other_src_grid_ind]*rad |
---|
| 116 | remap_matrix [ind] = 1.0 - dist / delta_x |
---|
| 117 | |
---|
| 118 | print 'New ', new_src_grid_ind, other_src_grid_ind, dist, delta_x, remap_matrix[ind] |
---|
| 119 | print src_address[ind], new_src_grid_ind_1, src_grid_ind_2, src_grid_center_lon[src_address[ind]-1]*rad, src_grid_center_lat[src_address[ind]-1]*rad |
---|
| 120 | print dst_address[ind], dst_grid_ind_1, dst_grid_ind_2, dst_grid_center_lon[dst_address[ind]-1]*rad, dst_grid_center_lat[dst_address[ind]-1]*rad |
---|
| 121 | print dst_grid_center_lon[dst_address[ind]-1]*rad - src_grid_center_lon[src_address[ind]-1]*rad, |
---|
| 122 | print dst_grid_center_lat[dst_address[ind]-1]*rad - src_grid_center_lat[src_address[ind]-1]*rad |
---|
| 123 | |
---|
| 124 | ## Write corrected fields |
---|
| 125 | rmp_file.write ( remap_matrix ) |
---|
| 126 | rmp_file.write ( src_address ) |
---|
| 127 | |
---|
| 128 | rmp_file.close () |
---|
| 129 | |
---|
| 130 | ## That's all folks !!! |
---|