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 !!! |
---|