source: TOOLS/BUG_OASIS_BILINEAR/correct_OasisBilinear.py @ 3325

Last change on this file since 3325 was 3325, checked in by omamce, 5 years ago

O.M. : premier enregistrement.
Utilitaire pour corriger les poids du bilineaire

File size: 5.7 KB
Line 
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
25import cdms2, MV2, numpy as np
26import os, time, shutil
27from 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:
31cdms2.setNetcdfShuffleFlag     (0)
32cdms2.setNetcdfDeflateFlag     (0)
33cdms2.setNetcdfDeflateLevelFlag(0)
34
35#
36rad = np.rad2deg (1.0) 
37
38## Copy initial rmp file with bilinear weights
39shutil.copyfile ( 'rmp_tlmd_to_torc_BILINEAR_ComputedByOasis.nc' , 'rmp_tlmd_to_torc_BILINEAR_Corrected.nc' )
40
41## Read the new file
42rmp_file             = cdms2.open   ( 'rmp_tlmd_to_torc_BILINEAR_Corrected.nc' , 'r+')
43
44src_grid_center_lat  = rmp_file ( 'src_grid_center_lat' )
45src_grid_center_lon  = rmp_file ( 'src_grid_center_lon' )
46dst_grid_center_lat  = rmp_file ( 'dst_grid_center_lat' )
47dst_grid_center_lon  = rmp_file ( 'dst_grid_center_lon' )
48src_address          = rmp_file ( 'src_address'         )
49dst_address          = rmp_file ( 'dst_address'         )
50remap_matrix         = rmp_file ( 'remap_matrix'        )
51src_grid_dims        = rmp_file ( 'src_grid_dims'       )
52dst_grid_dims        = rmp_file ( 'dst_grid_dims'       )
53
54src_grid_corner_lat  = rmp_file ( 'src_grid_corner_lat' )
55src_grid_corner_lon  = rmp_file ( 'src_grid_corner_lon' )
56src_grid_area        = rmp_file ( 'src_grid_area'       )
57src_grid_imask       = rmp_file ( 'src_grid_imask'      )
58
59dst_grid_corner_lat  = rmp_file ( 'dst_grid_corner_lat' )
60dst_grid_corner_lon  = rmp_file ( 'dst_grid_corner_lon' )
61dst_grid_area        = rmp_file ( 'dst_grid_area'       )
62dst_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
71print 'Dimensions '
72print 'src_grid_size : ', src_grid_size, '  src_grid_dims : ', src_grid_dims
73print 'dst_grid_size : ', dst_grid_size, '  dst_grid_dims : ', dst_grid_dims
74print ' '
75
76## Get list of wrong weights
77list_ind = np.where ( np.abs(remap_matrix) > 1.0 )[0]
78print 'List of indexes of buggy points : ', list_ind
79print ' '
80
81for 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
125rmp_file.write ( remap_matrix )
126rmp_file.write ( src_address  )
127 
128rmp_file.close ()
129
130## That's all folks !!!
Note: See TracBrowser for help on using the repository browser.