source: CPL/oasis3-mct/branches/OASIS3-MCT_2.0_branch/examples/test_rmp_esmf/interp_esmf_unstructured_to_curv_one_step.ncl @ 4775

Last change on this file since 4775 was 4775, checked in by aclsce, 4 years ago
  • Imported oasis3-mct from Cerfacs svn server (not suppotred anymore).

The version has been extracted from ​https://oasis3mct.cerfacs.fr/svn/branches/OASIS3-MCT_2.0_branch/oasis3-mct@1818
File size: 16.4 KB
Line 
1;----------------------------------------------------------------------
2; This script regrids dummy data from a curvilinear grid to an
3; unstructured grid. Both grids are on one NetCDF file.
4;
5; There's an issue with this script if you try to apply SrcMask2D,
6; because it causes the destination grid to be reduced in the number
7; of points (24,572 to 24,067). The fix seems to be in not setting
8; SrcMask2D, but I would like to figure out why this is the case.
9;----------------------------------------------------------------------
10;--------------------------------------------
11;---- Predefined ESMF modules
12;--------------------------------------------
13load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
14load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
15load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
16;
17; Loading a personal ESMF regridding script to provide
18; a needed bug fix for curvilinear to unstructured.
19;
20;load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
21load "./ESMF_regridding.ncl"
22
23begin
24;--------------------------------------------
25;--- PATHS
26;--------------------------------------------
27  RUNDIR="wkdir_esmf_"+str_lower(srcgrid)+"_"+str_lower(tgtgrid)+"_"+str_lower(interp_method)
28  PATH1="/home/evian/coquart/ESMF/ESMF_with_NCL/toy_calcul_weights_with_esmf/RUNDIR/"
29  PATH2="/home/evian/coquart/ESMF/ESMF_with_NCL/toy_calcul_weights_with_esmf/RUNDIR/"
30  PATH1="./"
31  PATH2="./"
32
33; WARNING : DOES NOT WORK YET FOR CONSERVATIVE REMAPPING
34  if (str_lower(interp_method) .eq. str_lower("conserve")) then
35     print("Does not work yet for conservative remapping")
36     exit
37  end if
38;--------------------------------------------
39;--- Source and Target grid Files used with oasis3-mct for SCRIP
40;--------------------------------------------
41  gridsfile = PATH1+"grids.nc"    ; File with grids at OASIS3-MCT description
42  masksfile = PATH1+"masks.nc"    ; File with masks at OASIS3-MCT description
43
44;--------------------------------------------
45;--- Files
46;-------------------------------------------- 
47  funcini     = PATH2+"fldin.nc"     ; File with the initial field
48  funcinterp  = PATH2+"fldou.nc"     ; File with the remapping field, error
49  srcGridName = PATH2+"source_grid_esmf.nc"
50  tgtGridName = PATH2+"target_grid_esmf.nc"
51  wgtFile     = PATH2+"rmp_"+srcgrid+"to"+tgtgrid+"_" + interp_method + "_withfillValue.nc"
52  wgtFile1    = PATH2+"rmp_"+srcgrid+"to"+tgtgrid+"_neareststod.nc"
53  wgtFilennei = PATH2+"rmp_"+srcgrid+"to"+tgtgrid+"_" + interp_method + ".nc"
54  oa3wgtFile  = PATH2+"rmp_esmf_for_oa3mct_"+srcgrid+"to"+tgtgrid+"_" + interp_method + ".nc"
55
56;--------------------------------------------
57;--- Constants
58;--------------------------------------------
59  msk_field_value = -1
60  msk_err_value = -10000
61  weight_fill_val = -50000
62  pi = 3.14159265359
63  dp_length = 1.2*pi
64  dp_conv   = pi/180.
65
66;--------------------------------------------
67;--- Cleaning: remove any pre-existing file
68;--------------------------------------------
69  system("/bin/rm -f "+funcini)
70  system("/bin/rm -f "+funcinterp)
71  system("/bin/rm -f "+srcGridName)
72  system("/bin/rm -f "+tgtGridName)
73  system("/bin/rm -f "+wgtFile) 
74  system("/bin/rm -f "+wgtFile1) 
75  system("/bin/rm -f "+wgtFilennei)
76  system("/bin/rm -f "+oa3wgtFile)
77  system("/bin/rm -f PET0.RegridWeightGen.Log")
78
79;--------------------------------------------
80;--- Coordinate names
81;--------------------------------------------
82  srcgrdlon = srcgrid+"_lon"
83  srcgrdlat = srcgrid+"_lat"
84  srcgrdclo = srcgrid+"_clo"
85  srcgrdcla = srcgrid+"_cla"
86  srcgrdmsk = srcgrid+"_msk"
87  tgtgrdlon = tgtgrid+"_lon"
88  tgtgrdlat = tgtgrid+"_lat"
89  tgtgrdclo = tgtgrid+"_clo"
90  tgtgrdcla = tgtgrid+"_cla"
91  tgtgrdmsk = tgtgrid+"_msk"
92
93;--------------------------------------------
94;--- Read grids, masks, areas arrays
95;--------------------------------------------
96  gfile   = addfile(gridsfile,"r")
97  mfile   = addfile(masksfile,"r")
98 
99  nt=0
100  src_lon  = gfile->$srcgrdlon$(nt,:)
101  src_lat  = gfile->$srcgrdlat$(nt,:)
102  src_clo  = gfile->$srcgrdclo$(:,nt,:)
103  src_cla  = gfile->$srcgrdcla$(:,nt,:)
104  src_msk  = mfile->$srcgrdmsk$(nt,:)
105  tgt_lon  = gfile->$tgtgrdlon$
106  tgt_lat  = gfile->$tgtgrdlat$
107  tgt_clo  = gfile->$tgtgrdclo$
108  tgt_cla  = gfile->$tgtgrdcla$
109  tgt_msk  = mfile->$tgtgrdmsk$
110
111  nlatsrc=dimsizes(src_lat)
112  nlonsrc=dimsizes(src_lon)
113  nlattgt=dimsizes(tgt_lat(:,0))
114  nlontgt=dimsizes(tgt_lon(0,:))
115
116;--------------------------------------------
117;--- Calculate the function to regrid
118;--------------------------------------------
119  func_ana_src = 2. - cos(pi*(acos(cos(src_lat*dp_conv)*cos(src_lon*dp_conv))/dp_length))
120
121;--------------------------------------------
122;--- Calculate the function on the destination grid for the error
123;--------------------------------------------
124  func_ana_tgt = 2. - cos(pi*(acos(cos(tgt_lat*dp_conv)*cos(tgt_lon*dp_conv))/dp_length))
125
126;----------------------------------------------------------------------
127;--- Regridding all in one step
128;----------------------------------------------------------------------
129  Opt                   = True
130  Opt@SrcGridName       = srcGridName
131  Opt@SrcGridLat        = src_lat
132  Opt@SrcGridLon        = src_lon
133  Opt@SrcMask2D         = 1-src_msk     ; Don't set this, otherwise you may get error.
134  Opt@DstGridName       = tgtGridName
135  Opt@DstGridLat        = tgt_lat
136  Opt@DstGridLon        = tgt_lon   
137  Opt@DstMask2D         = 1-tgt_msk
138
139  Opt@InterpMethod      = str_lower(interp_method)
140  Opt@WgtFileName       = wgtFile
141
142  if(Opt@InterpMethod.eq.str_lower("conserve")) then
143    Opt@SrcGridCornerLat  = ndtooned(src_cla)
144    Opt@SrcGridCornerLon  = ndtooned(src_clo)
145  end if
146
147  if ( interp_method .eq. str_lower("conserve") ) then
148    dims = getfilevardims(gfile,tgtgrdcla)
149;---Make sure grid is nlat x nlon x 4
150    Opt@DstGridCornerLat  = gfile->$tgtgrdcla$($dims(1)$|:,$dims(2)$|:,$dims(0)$|:)     
151    Opt@DstGridCornerLon  = gfile->$tgtgrdclo$($dims(1)$|:,$dims(2)$|:,$dims(0)$|:)
152    printMinMax(Opt@DstGridCornerLat,0)
153    printMinMax(Opt@DstGridCornerLon,0)
154  end if
155
156
157
158  Opt@ForceOverwrite    = True
159  Opt@PrintTimings      = True
160  Opt@Debug             = True
161;
162;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
163; Points that receive no values are at _FillValue=9.96920996838687e+36 (blank points)
164  func_regrid = ESMF_regrid(func_ana_src,Opt)
165  printVarSummary(func_regrid)
166  undef_val=func_regrid@_FillValue
167  printVarSummary(undef_val)
168  delete (Opt)
169  func_regrid_msk_fillvalue  = where((1-tgt_msk) .eq. 1, func_regrid, msk_field_value)
170;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
171;--------------------------------------------------------------------------
172;   Do the regridding "neareststod" all in one step to fill blank points
173;--------------------------------------------------------------------------
174  Opt=True
175  Opt@SrcGridName    = srcGridName
176  Opt@DstGridName    = tgtGridName
177  Opt@WgtFileName    = wgtFile1
178  Opt@SrcGridLat     = src_lat
179  Opt@SrcGridLon     = src_lon
180  Opt@DstGridLat     = tgt_lat
181  Opt@DstGridLon     = tgt_lon
182  Opt@DstMask2D      = 1-tgt_msk
183  Opt@SrcMask2D      = 1-src_msk
184
185  Opt@InterpMethod   = "neareststod"
186;  Opt@Pole           = "none"
187  Opt@ForceOverwrite = True
188  Opt@PrintTimings   = True
189  Opt@Debug          = True
190
191;---You can choose to skip any of these three steps.
192;  Opt@SkipSrcGrid   = True    ; Will assume source grid file has been generated
193;  Opt@SkipDstGrid   = True    ; Will assume target grid file has been generated
194;  Opt@SkipWgtGen    = True    ; Will assume weights file has been generated
195
196  func_neareststod = ESMF_regrid(func_ana_src,Opt)
197  func_neareststod = where((1-tgt_msk) .eq. 1, func_neareststod, msk_field_value)
198;----------------------------------------------------------------------
199;--- Combine interp_method + neareststod
200;----------------------------------------------------------------------
201  wfile          = addfile(wgtFile,"r") 
202;
203  S_regrid       = wfile->$"S"$
204  dst_lat_w      = wfile->$"yc_b"$
205  dst_lon_w      = wfile->$"xc_b"$
206  src_lat_w      = wfile->$"yc_a"$
207  src_lon_w      = wfile->$"xc_a"$
208  dst_address    = wfile->$"row"$
209  src_address    = wfile->$"col"$
210  frac_dst       = wfile->$"frac_b"$
211  numlinks      = dimsizes(S_regrid)
212  dst_size      = dimsizes(dst_lat_w)
213  src_size      = dimsizes(src_lat_w)
214  print("dst_size_w")
215  print(dst_size)
216;
217  wfile1            = addfile(wgtFile1,"r")
218  S_neareststod     = wfile1->$"S"$
219  dst_lat_n         = wfile1->$"yc_b"$
220  dst_lon_n         = wfile1->$"xc_b"$
221  src_lat_n         = wfile1->$"yc_a"$
222  dst_neareststod   = wfile1->$"row"$
223  src_neareststod   = wfile1->$"col"$
224  links_neareststod = dimsizes(S_neareststod)
225  dst_size_n    = dimsizes(dst_lat_n)
226  print("dst_size_n")
227  print(dst_size_n)
228;
229  if (interp_method .eq. str_lower("neareststod")) then
230;--- Just write the new weights into the remapping file
231     print("Nothing to do")
232  else
233      func_oned                = ndtooned(func_regrid_msk_fillvalue)
234      dst_missing              = new((/dst_size/),"integer")
235      dst_missing@_FillValue   = doubletointeger(undef_val)
236
237      count = -1
238      do i=0,dst_size-1
239         if (ismissing(func_oned(i))) then
240            count = count+1
241            dst_missing(count) = i+1
242         end if
243      end do
244; count > 0 : there are some points that need nearest neighbours
245      if (count .ge. 0) then
246      dst_missing_good    = new((/count+1/),"integer")
247      do i=0,count
248         dst_missing_good(i) = dst_missing(i)
249      end do
250      delete(dst_missing)
251      print("count")
252      print(count)
253      print("dst_missing_good")
254;      print(dst_missing_good)
255;
256;--- Replace weights where function = _FillValue by 1. (neareststod)
257      dst_address_new    = new((/numlinks+count+1/),"integer")
258      src_address_new    = new((/numlinks+count+1/),"integer")
259      S_regrid_new       = new((/numlinks+count+1/),double)
260      print("New size numlinks+count")
261      print(numlinks+count+1)
262      dst_address_new(0:numlinks-1)  = dst_address
263      src_address_new(0:numlinks-1)  = src_address
264      S_regrid_new(0:numlinks-1)     = S_regrid
265      do j=0,count
266         do i=0,links_neareststod-1
267            if ( dst_neareststod(i) .eq. dst_missing_good(j) ) then
268                dst_address_new(numlinks+j) = dst_missing_good(j)
269                src_address_new(numlinks+j) = src_neareststod(i)
270                S_regrid_new(numlinks+j)    = S_neareststod(i)
271                break
272            end if
273         end do
274     end do
275     delete(dst_missing_good)
276; count = 0
277     else
278     dst_address_new    = new((/numlinks/),"integer")
279     src_address_new    = new((/numlinks/),"integer")
280     S_regrid_new       = new((/numlinks/),double)
281     dst_address_new    = dst_address
282     src_address_new    = src_address
283     S_regrid_new       = S_regrid
284     end if
285
286     delete(S_regrid)
287     delete(dst_address)
288     delete(src_address)
289
290     igood = ind (.not. ismissing(S_regrid_new))
291     S_regrid = S_regrid_new(igood)
292     src_address = src_address_new(igood)
293     dst_address = dst_address_new(igood)
294
295; endif neareststod
296     delete(S_regrid_new)
297     delete(dst_address_new)
298     delete(src_address_new)
299     end if
300
301;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
302; Recalculate the interpolated function with the new weights with ESMF
303; Points that receive no values are at _FillValue=9.96920996838687e+36 (blank points)
304  if (interp_method .eq. str_lower("neareststod")) then
305;--- Just write the new weights into the remapping file
306     print("Nothing to do")
307  else
308  filennei = addfile(wgtFilennei,"c")
309  globalAtt               = True
310  globalAtt@title         = "ESMF Offline Regridding Weight Generator"
311  globalAtt@normalization = "destarea"
312  if (interp_method .eq. str_lower("conserve")) then
313     globalAtt@map_method = "Conservative remapping"
314     globalAtt@ESMF_regrid_method = "First-order Conservative"
315  end if
316  if (interp_method .eq. str_lower("bilinear")) then
317     globalAtt@map_method = "Bilinear remapping"
318     globalAtt@ESMF_regrid_method = "Bilinear"
319  end if
320  if (interp_method .eq. str_lower("patch")) then
321     globalAtt@map_method = "Bilinear remapping"
322     globalAtt@ESMF_regrid_method = "Higher-order Patch"
323  end if
324  if (interp_method .eq. str_lower("neareststod")) then
325     globalAtt@map_method = "Bilinear remapping"
326     globalAtt@ESMF_regrid_method = "Nearest source to destination"
327  end if
328  globalAtt@conventions   = "NCAR-CSM"
329  globalAtt@domain_a      = "source_grid_file.nc"
330  globalAtt@domain_b      = "destination_grid_file.nc"
331  globalAtt@grid_file_src = "source_grid_file.nc"
332  globalAtt@grid_file_dst = "destination_grid_file.nc"
333  globalAtt@CVS_revision  = "6.3.0r"
334;
335  fileattdef( filennei, globalAtt )
336  vars = getfilevarnames (wfile)
337  nv=dimsizes(vars)
338  print(nv)
339  S_nnei=S_regrid
340  num_links= dimsizes(S_regrid)
341  var_temp=new((/1/),integer)
342  var_temp(0)=1
343  src_address!0="n_s"
344  dst_address!0="n_s"
345  S_nnei!0="n_s"
346  var_temp!0="num_wgts"
347  delete(src_address@_FillValue)
348  delete(dst_address@_FillValue)
349  delete(S_nnei@_FillValue)
350  do i = 0,nv-4
351     filennei->$vars(i)$=wfile->$vars(i)$
352  end do
353  filennei->$"col"$=src_address
354  filennei->$"row"$=dst_address
355  filennei->$"S"$=S_nnei
356  filennei->$"var_temp"$=var_temp
357
358  Opt                = True
359  Opt@Debug          = True
360  Opt@PrintTimings   = True
361  func_regrid_new = ESMF_regrid_with_weights(func_ana_src,wgtFilennei,Opt)
362;  func_regrid_new = ESMF_regrid_with_weights(func_ana_src,wgtFile,Opt)
363  printVarSummary(func_regrid_new)
364
365  func_regrid_msk_new  = where((1-tgt_msk) .eq. 1, func_regrid_new, msk_field_value)
366
367  delete(var_temp)
368  end if
369;--------------------------------------------
370;--- Calculate the error on the target grid
371;--------------------------------------------
372  func_ana_tgt       = where((1-tgt_msk) .eq. 1, func_ana_tgt, msk_field_value)
373  if ( (interp_method .eq. "conserve") .or. (interp_method .eq. "bilinear") .or. (interp_method .eq. "patch") ) then
374       error = where((1-tgt_msk) .eq. 1, abs(((func_ana_tgt - func_regrid_msk_new)/func_regrid_msk_new))*100, msk_err_value)
375  else
376       error = where((1-tgt_msk) .eq. 1, abs(((func_ana_tgt - func_regrid_msk_fillvalue)/func_regrid_msk_fillvalue))*100, msk_err_value)
377  end if
378  print("Max of the error (equal -10000 on land)")
379  printMinMax(error,True)
380  error_min = where(error .eq. msk_err_value, - msk_err_value, error)
381  print("Min of the error (equal 10000 on land)")
382  printMinMax(error_min,True)
383
384;--------------------------------------------
385;--- Output files: to plot ESMF results and create the rmp file for OASIS3-MCT
386;--------------------------------------------
387  ifile = addfile(funcini,"c")
388  ifile->$srcgrdlon$=src_lon
389  ifile->$srcgrdlat$=src_lat
390  ifile->$"field_in"$=(func_ana_src)
391
392  ofile = addfile(funcinterp,"c")
393  ofile->$tgtgrdlon$=tgt_lon
394  ofile->$tgtgrdlat$=tgt_lat
395  ofile->$"field_ana"$=(func_ana_tgt)
396  ofile->$"field_undef"$=(func_regrid)
397  ofile->$"field_msk_undef"$=(func_regrid_msk_fillvalue)
398  if ( (interp_method .eq. "conserve") .or. (interp_method .eq. "bilinear") .or. (interp_method .eq. "patch") ) then
399     ofile->$"field_ou"$=(func_regrid_msk_new)
400  end if
401  ofile->$"error"$=(error)
402
403;---- Write the rmp file for OASIS3-MCT
404; if bilinear interpolation : num_wgts=1
405; if conservative interpolation : num_wgts=3
406; if nearest neighbour : num_wgts=1
407
408;--- Normalize to obtain Fracarea option for OASIS3-MCT for conservative interpolation
409  if (interp_method .eq. str_lower("conserve")) then
410     do i=0,num_links-1
411        if (frac_dst(dst_address(i)-1) .ne. 0.0) then
412           S_regrid(i)=S_regrid(i)/frac_dst(dst_address(i)-1)
413        end if
414     end do
415  end if
416
417    oa3file = addfile(oa3wgtFile,"c")
418    dst_lat_w!0="dst_grid_size"
419    src_lat_w!0="src_grid_size"
420    src_address!0="num_links"
421    dst_address!0="num_links"
422    oa3file->$"src_grid_center_lat"$=src_lat_w
423    oa3file->$"dst_grid_center_lat"$=dst_lat_w
424    oa3file->$"src_address"$=src_address
425    oa3file->$"dst_address"$=dst_address
426    if (interp_method .eq. "conserve") then
427       var_temp=S_regrid
428       num_links=dimsizes(var_temp)
429       num_wgts=3
430       var = new ((/num_links,num_wgts/),double)
431       var(:,0) = var_temp
432       var(:,1) = 0.
433       var(:,2) = 0.
434       var!0="num_links"
435       var!1="num_wgts"
436       oa3file->$"remap_matrix"$=var
437  else
438       var_temp=S_regrid
439       num_links=dimsizes(var_temp)
440       num_wgts=1
441       var = new ((/num_links,num_wgts/),double)
442       var(:,0) = var_temp
443       var!0="num_links"
444       var!1="num_wgts"
445       oa3file->$"remap_matrix"$=var
446  end if
447
448end
Note: See TracBrowser for help on using the repository browser.