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