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

Last change on this file since 4775 was 4775, checked in by aclsce, 5 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.5 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$
101  src_lat  = gfile->$srcgrdlat$
102  src_clo  = gfile->$srcgrdclo$
103  src_cla  = gfile->$srcgrdcla$
104  src_msk  = mfile->$srcgrdmsk$
105  tgt_lon  = gfile->$tgtgrdlon$(nt,:)
106  tgt_lat  = gfile->$tgtgrdlat$(nt,:)
107  tgt_clo  = gfile->$tgtgrdclo$(:,nt,:)
108  tgt_cla  = gfile->$tgtgrdcla$(:,nt,:)
109  tgt_msk  = mfile->$tgtgrdmsk$(nt,:)
110
111  nlatsrc=dimsizes(src_lat(:,0))
112  nlonsrc=dimsizes(src_lon(0,:))
113  nlattgt=dimsizes(tgt_lat)
114  nlontgt=dimsizes(tgt_lon)
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  Opt@DstGridType       = "unstructured"
139
140  Opt@InterpMethod      = str_lower(interp_method)
141  Opt@WgtFileName       = wgtFile
142
143  if ( interp_method .eq. str_lower("conserve") ) then
144    dims = getfilevardims(gfile,srcgrdcla)
145;---Make sure grid is nlat x nlon x 4
146    Opt@SrcGridCornerLat  = gfile->$srcgrdcla$($dims(1)$|:,$dims(2)$|:,$dims(0)$|:)     
147    Opt@SrcGridCornerLon  = gfile->$srcgrdclo$($dims(1)$|:,$dims(2)$|:,$dims(0)$|:)
148    printMinMax(Opt@SrcGridCornerLat,0)
149    printMinMax(Opt@SrcGridCornerLon,0)
150  end if
151
152  if(Opt@InterpMethod.eq.str_lower("conserve")) then
153    Opt@DstGridCornerLat  = ndtooned(tgt_cla)
154    Opt@DstGridCornerLon  = ndtooned(tgt_clo)
155  end if
156
157  Opt@ForceOverwrite    = True
158  Opt@PrintTimings      = True
159  Opt@Debug             = True
160;
161;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
162; Points that receive no values are at _FillValue=9.96920996838687e+36 (blank points)
163  func_regrid = ESMF_regrid(func_ana_src,Opt)
164  printVarSummary(func_regrid)
165  undef_val=func_regrid@_FillValue
166  printVarSummary(undef_val)
167  delete (Opt)
168  func_regrid_msk_fillvalue  = where((1-tgt_msk) .eq. 1, func_regrid, msk_field_value)
169;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
170;--------------------------------------------------------------------------
171;   Do the regridding "neareststod" all in one step to fill blank points
172;--------------------------------------------------------------------------
173  Opt=True
174  Opt@SrcGridName    = srcGridName
175  Opt@DstGridName    = tgtGridName
176  Opt@WgtFileName    = wgtFile1
177  Opt@SrcGridLat     = src_lat
178  Opt@SrcGridLon     = src_lon
179  Opt@DstGridLat     = tgt_lat
180  Opt@DstGridLon     = tgt_lon
181  Opt@DstMask2D      = 1-tgt_msk
182  Opt@SrcMask2D      = 1-src_msk
183
184  Opt@DstGridType       = "unstructured"
185
186  Opt@InterpMethod   = "neareststod"
187;  Opt@Pole           = "none"
188  Opt@ForceOverwrite = True
189  Opt@PrintTimings   = True
190  Opt@Debug          = True
191
192;---You can choose to skip any of these three steps.
193;  Opt@SkipSrcGrid   = True    ; Will assume source grid file has been generated
194;  Opt@SkipDstGrid   = True    ; Will assume target grid file has been generated
195;  Opt@SkipWgtGen    = True    ; Will assume weights file has been generated
196
197  func_neareststod = ESMF_regrid(func_ana_src,Opt)
198  func_neareststod = where((1-tgt_msk) .eq. 1, func_neareststod, msk_field_value)
199;----------------------------------------------------------------------
200;--- Combine interp_method + neareststod
201;----------------------------------------------------------------------
202  wfile          = addfile(wgtFile,"r") 
203;
204  S_regrid       = wfile->$"S"$
205  dst_lat_w      = wfile->$"yc_b"$
206  dst_lon_w      = wfile->$"xc_b"$
207  src_lat_w      = wfile->$"yc_a"$
208  src_lon_w      = wfile->$"xc_a"$
209  dst_address    = wfile->$"row"$
210  src_address    = wfile->$"col"$
211  frac_dst       = wfile->$"frac_b"$
212  numlinks      = dimsizes(S_regrid)
213  dst_size      = dimsizes(dst_lat_w)
214  src_size      = dimsizes(src_lat_w)
215  print("dst_size_w")
216  print(dst_size)
217;
218  wfile1            = addfile(wgtFile1,"r")
219  S_neareststod     = wfile1->$"S"$
220  dst_lat_n         = wfile1->$"yc_b"$
221  dst_lon_n         = wfile1->$"xc_b"$
222  src_lat_n         = wfile1->$"yc_a"$
223  dst_neareststod   = wfile1->$"row"$
224  src_neareststod   = wfile1->$"col"$
225  links_neareststod = dimsizes(S_neareststod)
226  dst_size_n    = dimsizes(dst_lat_n)
227  print("dst_size_n")
228  print(dst_size_n)
229;
230  if (interp_method .eq. str_lower("neareststod")) then
231;--- Just write the new weights into the remapping file
232     print("Nothing to do")
233  else
234      func_oned                = ndtooned(func_regrid_msk_fillvalue)
235      dst_missing              = new((/dst_size/),"integer")
236      dst_missing@_FillValue   = doubletointeger(undef_val)
237
238      count = -1
239      do i=0,dst_size-1
240         if (ismissing(func_oned(i))) then
241            count = count+1
242            dst_missing(count) = i+1
243         end if
244      end do
245; count > 0 : there are some points that need nearest neighbours
246      if (count .ge. 0) then
247      dst_missing_good    = new((/count+1/),"integer")
248      do i=0,count
249         dst_missing_good(i) = dst_missing(i)
250      end do
251      delete(dst_missing)
252      print("count")
253      print(count)
254      print("dst_missing_good")
255;      print(dst_missing_good)
256;
257;--- Replace weights where function = _FillValue by 1. (neareststod)
258      dst_address_new    = new((/numlinks+count+1/),"integer")
259      src_address_new    = new((/numlinks+count+1/),"integer")
260      S_regrid_new       = new((/numlinks+count+1/),double)
261      print("New size numlinks+count")
262      print(numlinks+count+1)
263      dst_address_new(0:numlinks-1)  = dst_address
264      src_address_new(0:numlinks-1)  = src_address
265      S_regrid_new(0:numlinks-1)     = S_regrid
266      do j=0,count
267         do i=0,links_neareststod-1
268            if ( dst_neareststod(i) .eq. dst_missing_good(j) ) then
269                dst_address_new(numlinks+j) = dst_missing_good(j)
270                src_address_new(numlinks+j) = src_neareststod(i)
271                S_regrid_new(numlinks+j)    = S_neareststod(i)
272                break
273            end if
274         end do
275     end do
276     delete(dst_missing_good)
277; count = 0
278     else
279     dst_address_new    = new((/numlinks/),"integer")
280     src_address_new    = new((/numlinks/),"integer")
281     S_regrid_new       = new((/numlinks/),double)
282     dst_address_new    = dst_address
283     src_address_new    = src_address
284     S_regrid_new       = S_regrid
285     end if
286
287     delete(S_regrid)
288     delete(dst_address)
289     delete(src_address)
290
291     igood = ind (.not. ismissing(S_regrid_new))
292     S_regrid = S_regrid_new(igood)
293     src_address = src_address_new(igood)
294     dst_address = dst_address_new(igood)
295
296; endif neareststod
297     delete(S_regrid_new)
298     delete(dst_address_new)
299     delete(src_address_new)
300     end if
301
302;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
303; Recalculate the interpolated function with the new weights with ESMF
304; Points that receive no values are at _FillValue=9.96920996838687e+36 (blank points)
305  if (interp_method .eq. str_lower("neareststod")) then
306;--- Just write the new weights into the remapping file
307     print("Nothing to do")
308  else
309  filennei = addfile(wgtFilennei,"c")
310  globalAtt               = True
311  globalAtt@title         = "ESMF Offline Regridding Weight Generator"
312  globalAtt@normalization = "destarea"
313  if (interp_method .eq. str_lower("conserve")) then
314     globalAtt@map_method = "Conservative remapping"
315     globalAtt@ESMF_regrid_method = "First-order Conservative"
316  end if
317  if (interp_method .eq. str_lower("bilinear")) then
318     globalAtt@map_method = "Bilinear remapping"
319     globalAtt@ESMF_regrid_method = "Bilinear"
320  end if
321  if (interp_method .eq. str_lower("patch")) then
322     globalAtt@map_method = "Bilinear remapping"
323     globalAtt@ESMF_regrid_method = "Higher-order Patch"
324  end if
325  if (interp_method .eq. str_lower("neareststod")) then
326     globalAtt@map_method = "Bilinear remapping"
327     globalAtt@ESMF_regrid_method = "Nearest source to destination"
328  end if
329  globalAtt@conventions   = "NCAR-CSM"
330  globalAtt@domain_a      = "source_grid_file.nc"
331  globalAtt@domain_b      = "destination_grid_file.nc"
332  globalAtt@grid_file_src = "source_grid_file.nc"
333  globalAtt@grid_file_dst = "destination_grid_file.nc"
334  globalAtt@CVS_revision  = "6.3.0r"
335;
336  fileattdef( filennei, globalAtt )
337  vars = getfilevarnames (wfile)
338  nv=dimsizes(vars)
339  print(nv)
340  S_nnei=S_regrid
341  num_links= dimsizes(S_regrid)
342  var_temp=new((/1/),integer)
343  var_temp(0)=1
344  src_address!0="n_s"
345  dst_address!0="n_s"
346  S_nnei!0="n_s"
347  var_temp!0="num_wgts"
348  delete(src_address@_FillValue)
349  delete(dst_address@_FillValue)
350  delete(S_nnei@_FillValue)
351  do i = 0,nv-4
352     filennei->$vars(i)$=wfile->$vars(i)$
353  end do
354  filennei->$"col"$=src_address
355  filennei->$"row"$=dst_address
356  filennei->$"S"$=S_nnei
357  filennei->$"var_temp"$=var_temp
358
359  Opt                = True
360  Opt@Debug          = True
361  Opt@PrintTimings   = True
362  func_regrid_new = ESMF_regrid_with_weights(func_ana_src,wgtFilennei,Opt)
363;  func_regrid_new = ESMF_regrid_with_weights(func_ana_src,wgtFile,Opt)
364  printVarSummary(func_regrid_new)
365
366  func_regrid_msk_new  = where((1-tgt_msk) .eq. 1, func_regrid_new, msk_field_value)
367
368  delete(var_temp)
369  end if
370;--------------------------------------------
371;--- Calculate the error on the target grid
372;--------------------------------------------
373  func_ana_tgt       = where((1-tgt_msk) .eq. 1, func_ana_tgt, msk_field_value)
374  if ( (interp_method .eq. "conserve") .or. (interp_method .eq. "bilinear") .or. (interp_method .eq. "patch") ) then
375       error = where((1-tgt_msk) .eq. 1, abs(((func_ana_tgt - func_regrid_msk_new)/func_regrid_msk_new))*100, msk_err_value)
376  else
377       error = where((1-tgt_msk) .eq. 1, abs(((func_ana_tgt - func_regrid_msk_fillvalue)/func_regrid_msk_fillvalue))*100, msk_err_value)
378  end if
379  print("Max of the error (equal -10000 on land)")
380  printMinMax(error,True)
381  error_min = where(error .eq. msk_err_value, - msk_err_value, error)
382  print("Min of the error (equal 10000 on land)")
383  printMinMax(error_min,True)
384
385;--------------------------------------------
386;--- Output files: to plot ESMF results and create the rmp file for OASIS3-MCT
387;--------------------------------------------
388  ifile = addfile(funcini,"c")
389  ifile->$srcgrdlon$=src_lon
390  ifile->$srcgrdlat$=src_lat
391  ifile->$"field_in"$=(func_ana_src)
392
393  ofile = addfile(funcinterp,"c")
394  ofile->$tgtgrdlon$=tgt_lon
395  ofile->$tgtgrdlat$=tgt_lat
396  ofile->$"field_ana"$=(func_ana_tgt)
397  ofile->$"field_undef"$=(func_regrid)
398  ofile->$"field_msk_undef"$=(func_regrid_msk_fillvalue)
399  if ( (interp_method .eq. "conserve") .or. (interp_method .eq. "bilinear") .or. (interp_method .eq. "patch") ) then
400     ofile->$"field_ou"$=(func_regrid_msk_new)
401  end if
402  ofile->$"error"$=(error)
403
404;---- Write the rmp file for OASIS3-MCT
405; if bilinear interpolation : num_wgts=1
406; if conservative interpolation : num_wgts=3
407; if nearest neighbour : num_wgts=1
408
409;--- Normalize to obtain Fracarea option for OASIS3-MCT for conservative interpolation
410  if (interp_method .eq. str_lower("conserve")) then
411     do i=0,num_links-1
412        if (frac_dst(dst_address(i)-1) .ne. 0.0) then
413           S_regrid(i)=S_regrid(i)/frac_dst(dst_address(i)-1)
414        end if
415     end do
416  end if
417
418    oa3file = addfile(oa3wgtFile,"c")
419    dst_lat_w!0="dst_grid_size"
420    src_lat_w!0="src_grid_size"
421    src_address!0="num_links"
422    dst_address!0="num_links"
423    oa3file->$"src_grid_center_lat"$=src_lat_w
424    oa3file->$"dst_grid_center_lat"$=dst_lat_w
425    oa3file->$"src_address"$=src_address
426    oa3file->$"dst_address"$=dst_address
427    if (interp_method .eq. "conserve") then
428       var_temp=S_regrid
429       num_links=dimsizes(var_temp)
430       num_wgts=3
431       var = new ((/num_links,num_wgts/),double)
432       var(:,0) = var_temp
433       var(:,1) = 0.
434       var(:,2) = 0.
435       var!0="num_links"
436       var!1="num_wgts"
437       oa3file->$"remap_matrix"$=var
438  else
439       var_temp=S_regrid
440       num_links=dimsizes(var_temp)
441       num_wgts=1
442       var = new ((/num_links,num_wgts/),double)
443       var(:,0) = var_temp
444       var!0="num_links"
445       var!1="num_wgts"
446       oa3file->$"remap_matrix"$=var
447  end if
448
449end
Note: See TracBrowser for help on using the repository browser.