1 | SUBROUTINE scriprmp (dst_array, src_array, src_size, dst_size, |
---|
2 | $ src_mask, dst_mask, |
---|
3 | $ src_lon, src_lat, nlon_src, nlat_src, |
---|
4 | $ dst_lon, dst_lat, nlon_dst, nlat_dst, |
---|
5 | $ map_method, cdgrdtyp, id_nosper, |
---|
6 | $ src_name, dst_name, |
---|
7 | $ normalize_opt, order, rst_type, n_srch_bins, |
---|
8 | $ lextrapdone, rl_varmul, id_scripvoi) |
---|
9 | C***** |
---|
10 | C ***************************** |
---|
11 | C * OASIS ROUTINE - LEVEL 3 * |
---|
12 | C * ------------- ------- * |
---|
13 | C ***************************** |
---|
14 | C |
---|
15 | C**** *scriprmp* - SCRIP remapping |
---|
16 | C |
---|
17 | C Purpose: |
---|
18 | C ------- |
---|
19 | C Main routine of SCRIP remapping |
---|
20 | C - finds out, whether remapping matrix exists |
---|
21 | C - drives calculation of missing remapping matrices |
---|
22 | C - performs matrix multiplication |
---|
23 | C |
---|
24 | C Interface: |
---|
25 | C --------- |
---|
26 | C *CALL* *scriprmp (dst_array, src_array, src_size, dst_size, |
---|
27 | C src_mask, dst_mask, |
---|
28 | C src_lon, src_lat, nlon_src, nlat_src, |
---|
29 | C dst_lon, dst_lat, nlon_dst, nlat_dst, |
---|
30 | C map_method, cdgrdtyp, id_nosper, |
---|
31 | C src_name, dst_name, |
---|
32 | C normalize_opt, order, rst_type, n_srch_bins, |
---|
33 | C lextrapdone, rl_varmul) |
---|
34 | C |
---|
35 | C Called from: |
---|
36 | C ----------- |
---|
37 | C interp |
---|
38 | C |
---|
39 | C Input: |
---|
40 | C ----- |
---|
41 | C src_array : field on source grid(real 1D) |
---|
42 | C src_size : source grid size (integer) |
---|
43 | C dst_size : target grid size (integer) |
---|
44 | C src_mask : source grid mask (INTEGER) |
---|
45 | C dst_mask : target grid mask (INTEGER) |
---|
46 | C src_lon : source grid longitudes (real 1D) |
---|
47 | C src_lat : source grid latitudes (real 1D) |
---|
48 | C nlon_src : number of source grid longitudes (integer) |
---|
49 | C nlat_src : number of source grid latitudes (integer) |
---|
50 | C dst_lon : target grid longitudes (real 1D) |
---|
51 | C dst_lat : target grid latitudes (real 1D) |
---|
52 | C nlon_dst : number of destination grid longitudes (integer) |
---|
53 | C nlat_dst : number of destination grid latitudes (integer) |
---|
54 | C map_method: remapping method (character*8) |
---|
55 | C cdgrdtyp : source grid type (character*8) |
---|
56 | C id_nosper : number of overlapping for source grid |
---|
57 | C src_name : source grid name (character*8) |
---|
58 | C dst_name : target grid name (character*8) |
---|
59 | C normalize_opt: option for normalization (character*8) |
---|
60 | C order : order of conservative remapping (character*8) |
---|
61 | C rst_type : type of scrip search restriction (character*8) |
---|
62 | C n_srch_bins : number of seach bins (integer) |
---|
63 | C lextrapdone : logical, true if EXTRAP done on field |
---|
64 | C rl_varmul : Gaussian variance (for GAUSWGT) |
---|
65 | C id_scripvoi : number of neighbour for DISTWGT and GAUSWGT |
---|
66 | C |
---|
67 | C Output: |
---|
68 | C ------ |
---|
69 | C dst_array : field on target grid (real 1D) |
---|
70 | C |
---|
71 | C Externals: |
---|
72 | C --------- |
---|
73 | C corners, scrip, gradient, gradient_bilin |
---|
74 | C |
---|
75 | C History: |
---|
76 | C ------- |
---|
77 | C Version Programmer Date Description |
---|
78 | C ------- ---------- ---- ----------- |
---|
79 | C 2.0 V.Gayler 2001/11/09 created |
---|
80 | C 2.5 D.Declat 2002/07/08 completed |
---|
81 | C 2.5 D.Declat 2002/08/01 the mask of the tgt grid |
---|
82 | C taken into account |
---|
83 | C 2.5 D.Declat 2002/08/09 'NONE' and 'CONSERV' for conserv |
---|
84 | C check pole |
---|
85 | C |
---|
86 | C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
87 | C* ---------------------------- Modules used ---------------------------- |
---|
88 | C |
---|
89 | USE grids |
---|
90 | USE remap_vars |
---|
91 | USE mod_unit |
---|
92 | USE mod_printing |
---|
93 | C |
---|
94 | C* ---------------------------- Implicit -------------------------------- |
---|
95 | C |
---|
96 | IMPLICIT NONE |
---|
97 | C |
---|
98 | C* ---------------------------- Include files --------------------------- |
---|
99 | C |
---|
100 | INCLUDE 'netcdf.inc' |
---|
101 | |
---|
102 | C |
---|
103 | C* ---------------------------- Intent In ------------------------------- |
---|
104 | C |
---|
105 | INTEGER (kind=int_kind) :: |
---|
106 | $ nlon_src, ! number of source grid longitudes |
---|
107 | $ nlat_src, ! number of source grid latitudes |
---|
108 | $ nlon_dst, ! number of destination grid longitudes |
---|
109 | $ nlat_dst, ! number of destination grid latitudes |
---|
110 | $ src_size, ! number of source grid cells |
---|
111 | $ dst_size, ! number of destination grid cells |
---|
112 | $ n_srch_bins, ! number of search bins fos SCRIP |
---|
113 | $ src_mask(src_size), ! source grid mask |
---|
114 | $ dst_mask(dst_size), ! target grid mask |
---|
115 | $ id_nosper ! number of overlapping points for source grid |
---|
116 | C |
---|
117 | REAL (kind=real_kind) :: |
---|
118 | $ src_array(src_size), ! source grid array |
---|
119 | $ src_lon(src_size), ! source grid longitudes |
---|
120 | $ src_lat(src_size), ! source grid latitudes |
---|
121 | $ dst_lon(dst_size), ! target grid longitudes |
---|
122 | $ dst_lat(dst_size) ! target grid latitudes |
---|
123 | C |
---|
124 | CHARACTER*8 :: |
---|
125 | $ map_method, ! remapping method |
---|
126 | $ cdgrdtyp, ! source grid type |
---|
127 | $ src_name, ! source grid name |
---|
128 | $ dst_name, ! target grid name |
---|
129 | $ normalize_opt, ! option for normalization |
---|
130 | $ order, ! order of conservative remapping |
---|
131 | $ rst_type ! type of scrip search restriction |
---|
132 | C |
---|
133 | LOGICAL :: |
---|
134 | $ lextrapdone ! logical, true if EXTRAP done on field |
---|
135 | C |
---|
136 | REAL (kind=real_kind) :: |
---|
137 | $ rl_varmul ! Gaussian variance (for GAUSWGT) |
---|
138 | C |
---|
139 | INTEGER (kind=int_kind) :: |
---|
140 | $ id_scripvoi ! number of neighbour for DISTWGT and GAUSWGT |
---|
141 | C |
---|
142 | C* ---------------------------- Intent Out ------------------------------- |
---|
143 | C |
---|
144 | REAL (kind=real_kind):: |
---|
145 | $ dst_array(dst_size) ! array on destination grid |
---|
146 | C |
---|
147 | C* ---------------------------- Local declarations ---------------------- |
---|
148 | C |
---|
149 | CHARACTER*12 :: |
---|
150 | $ cweight ! string for weights |
---|
151 | C |
---|
152 | CHARACTER*11 :: |
---|
153 | $ csrcadd, ! string for source grid addresses |
---|
154 | $ cdstadd ! string for destination grid addresses |
---|
155 | C |
---|
156 | CHARACTER*13 :: |
---|
157 | $ cdstare, ! string for destination grid area |
---|
158 | $ cdstfra ! string for destination grid frac |
---|
159 | C |
---|
160 | CHARACTER (char_len) :: |
---|
161 | $ crmpfile, ! name of the SCRIP matrix file |
---|
162 | $ cmapping ! mapping name |
---|
163 | C |
---|
164 | INTEGER (kind=int_kind) :: |
---|
165 | $ n, ! looping indicee |
---|
166 | $ num_links, ! number of intersections |
---|
167 | $ ncrn_src, ncrn_dst, ! number of grid cell corners |
---|
168 | $ src_rank, dst_rank, ! source / target grid rank |
---|
169 | $ num_wgts, ! number of weights |
---|
170 | $ src_dims(2), ! source grid dimensions |
---|
171 | $ dst_dims(2), ! target grid dimensions |
---|
172 | $ sou_mask(src_size), ! source grid mask |
---|
173 | $ tgt_mask(dst_size) ! target grid mask |
---|
174 | C |
---|
175 | REAL (kind=real_kind) :: |
---|
176 | $ dst_area(dst_size), ! target grid area |
---|
177 | $ dst_frac(dst_size), ! target grid frac |
---|
178 | $ dst_err(dst_size) ! target grid error |
---|
179 | C |
---|
180 | C* netCDF-declarations |
---|
181 | INTEGER (kind=int_kind) :: |
---|
182 | $ stat, ! netCDF error status |
---|
183 | $ nc_scpid, ! file id |
---|
184 | $ dimid, ! dimension id |
---|
185 | $ varid ! variable id |
---|
186 | C |
---|
187 | LOGICAL :: |
---|
188 | $ lcalc ! calculate matrix? |
---|
189 | C |
---|
190 | INTEGER (kind=int_kind), DIMENSION(:), ALLOCATABLE :: |
---|
191 | $ src_addr, ! addresses of source grid |
---|
192 | $ dst_addr ! addresses of destination grid |
---|
193 | C |
---|
194 | REAL (kind=real_kind), DIMENSION(:), ALLOCATABLE :: |
---|
195 | $ gradient_lat, ! latitudinal gradient (conservative rmp.) |
---|
196 | $ gradient_lon, ! longitudinal gradient (conservative rmp.) |
---|
197 | $ gradient_i, ! gradient in i direction (bilinear rmp.) |
---|
198 | $ gradient_j, ! gradient in j direction (bilinear rmp.) |
---|
199 | $ gradient_ij ! cross gradient (bilinear rmp.) |
---|
200 | C |
---|
201 | REAL (kind=real_kind), DIMENSION (:,:), ALLOCATABLE :: |
---|
202 | $ weights ! Remapping weights |
---|
203 | C |
---|
204 | REAL (kind=real_kind), DIMENSION (:,:), ALLOCATABLE :: |
---|
205 | $ src_corner_lon, ! longitudes of source grid corners |
---|
206 | $ src_corner_lat, ! latitudes of source grid corners |
---|
207 | $ dst_corner_lon, ! longitudes of destination grid corners |
---|
208 | $ dst_corner_lat ! latitudes of destination grid corners |
---|
209 | C |
---|
210 | C* pole contribution |
---|
211 | REAL (kind=real_kind) :: |
---|
212 | $ moy_tmp_S, moy_tmp_N, ! field average at the pole |
---|
213 | $ moy_err_S, moy_err_N, ! error average at the pole |
---|
214 | $ latpol_N, latpol_S |
---|
215 | C |
---|
216 | INTEGER (kind=int_kind) :: |
---|
217 | $ compt_S, compt_N ! number of cells representating a pole |
---|
218 | C |
---|
219 | C* ---------------------------- Poema verses ---------------------------- |
---|
220 | C |
---|
221 | C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
222 | C |
---|
223 | C* 1. Initialization |
---|
224 | C -------------- |
---|
225 | C |
---|
226 | IF (nlogprt .GE. 2) THEN |
---|
227 | WRITE (UNIT = nulou,FMT = *) ' ' |
---|
228 | WRITE (UNIT = nulou,FMT = *) ' ' |
---|
229 | WRITE (UNIT = nulou,FMT = *) |
---|
230 | $ ' ROUTINE scriprmp - Level ?' |
---|
231 | WRITE (UNIT = nulou,FMT = *) |
---|
232 | $ ' **************** *******' |
---|
233 | WRITE (UNIT = nulou,FMT = *) ' ' |
---|
234 | WRITE (UNIT = nulou,FMT = *) ' SCRIP remapping' |
---|
235 | WRITE (UNIT = nulou,FMT = *) ' ' |
---|
236 | ENDIF |
---|
237 | |
---|
238 | C |
---|
239 | C* -- get the name of the file containig the remapping matrix |
---|
240 | C |
---|
241 | SELECT CASE (map_method) |
---|
242 | CASE ('CONSERV') ! conservative remapping |
---|
243 | cmapping = |
---|
244 | $ src_name(1:4)//' to '//dst_name(1:4)//' '//map_method// |
---|
245 | $ ' '//normalize_opt(1:4)//' remapping' |
---|
246 | |
---|
247 | crmpfile = |
---|
248 | $ 'rmp_'//src_name(1:4)//'_to_'//dst_name(1:4)//'_'// |
---|
249 | $ map_method(1:7)//'_'//normalize_opt(1:8)//'.nc' |
---|
250 | CASE DEFAULT |
---|
251 | cmapping = |
---|
252 | $ src_name(1:4)//' to '//dst_name(1:4)//' '//map_method// |
---|
253 | $ ' remapping' |
---|
254 | |
---|
255 | crmpfile = |
---|
256 | $ 'rmp_'//src_name(1:4)//'_to_'//dst_name(1:4)//'_'// |
---|
257 | $ map_method(1:7)//'.nc' |
---|
258 | END SELECT |
---|
259 | C |
---|
260 | IF (nlogprt .GE. 2) THEN |
---|
261 | WRITE (UNIT = nulou,FMT = *) |
---|
262 | $ ' SCRIP filename : ', crmpfile |
---|
263 | WRITE (UNIT = nulou,FMT = *) ' ' |
---|
264 | ENDIF |
---|
265 | |
---|
266 | C**** |
---|
267 | C |
---|
268 | C* -- character strings of weights and addresses |
---|
269 | C |
---|
270 | csrcadd = 'src_address' |
---|
271 | cdstadd = 'dst_address' |
---|
272 | cweight = 'remap_matrix' |
---|
273 | cdstare = 'dst_grid_area' |
---|
274 | cdstfra = 'dst_grid_frac' |
---|
275 | C |
---|
276 | C* -- logical to calculate matrix |
---|
277 | C |
---|
278 | lcalc = .TRUE. |
---|
279 | C |
---|
280 | C* -- find out whether remapping file exists |
---|
281 | C |
---|
282 | stat = NF_OPEN(crmpfile, NF_NOWRITE, nc_scpid) |
---|
283 | IF (stat == NF_NOERR) THEN |
---|
284 | lcalc = .FALSE. |
---|
285 | C |
---|
286 | IF (nlogprt .GE. 2) THEN |
---|
287 | WRITE (UNIT = nulou,FMT = *) |
---|
288 | $ ' SCRIP file opened - no matrix calculation needed' |
---|
289 | WRITE (UNIT = nulou,FMT = *) ' ' |
---|
290 | ENDIF |
---|
291 | ENDIF |
---|
292 | C |
---|
293 | C* -- setting of the mask for the source AND the target grid |
---|
294 | C |
---|
295 | WHERE (src_mask .eq. 1) |
---|
296 | sou_mask = 0 |
---|
297 | END WHERE |
---|
298 | WHERE (src_mask .eq. 0) |
---|
299 | sou_mask = 1 |
---|
300 | END WHERE |
---|
301 | C |
---|
302 | WHERE (dst_mask .eq. 1) |
---|
303 | tgt_mask = 0 |
---|
304 | END WHERE |
---|
305 | WHERE (dst_mask .eq. 0) |
---|
306 | tgt_mask = 1 |
---|
307 | END WHERE |
---|
308 | |
---|
309 | IF (lcalc) THEN |
---|
310 | C |
---|
311 | C* 2. Calculate SCRIP remapping matrix |
---|
312 | C -------------------------------- |
---|
313 | C |
---|
314 | SELECT CASE (normalize_opt) |
---|
315 | CASE ('FRACNNEI') |
---|
316 | normalize_opt = 'FRACAREA' |
---|
317 | lfracnnei = .true. |
---|
318 | END SELECT |
---|
319 | C |
---|
320 | IF (nlogprt .GE. 2) THEN |
---|
321 | WRITE (UNIT = nulou,FMT = *) |
---|
322 | $ ' Calculation of SCRIP remapping matrix: method = ', |
---|
323 | $ map_method |
---|
324 | WRITE (UNIT = nulou,FMT = *) ' ' |
---|
325 | ENDIF |
---|
326 | |
---|
327 | C* -- get grid cell corners for conservative remapping |
---|
328 | C |
---|
329 | ncrn_src = 4. |
---|
330 | ncrn_dst = 4. |
---|
331 | C |
---|
332 | IF (map_method == 'CONSERV') THEN |
---|
333 | C |
---|
334 | ALLOCATE(src_corner_lon(ncrn_src,src_size), |
---|
335 | $ src_corner_lat(ncrn_src,src_size), |
---|
336 | $ dst_corner_lon(ncrn_dst,dst_size), |
---|
337 | $ dst_corner_lat(ncrn_dst,dst_size)) |
---|
338 | C |
---|
339 | CALL corners(nlon_src, nlat_src, ncrn_src, src_lon, src_lat, |
---|
340 | $ src_name, cdgrdtyp, src_corner_lon, src_corner_lat) |
---|
341 | CALL corners(nlon_dst, nlat_dst, ncrn_dst, dst_lon, dst_lat, |
---|
342 | $ dst_name, cdgrdtyp, dst_corner_lon, dst_corner_lat) |
---|
343 | ENDIF |
---|
344 | C |
---|
345 | C* -- initialization of grid arrays for SCRIP |
---|
346 | C |
---|
347 | src_dims(1) = nlon_src |
---|
348 | src_dims(2) = nlat_src |
---|
349 | dst_dims(1) = nlon_dst |
---|
350 | dst_dims(2) = nlat_dst |
---|
351 | src_rank = 2 |
---|
352 | dst_rank = 2 |
---|
353 | ! Modifier car src_rank et dst_rank n'est pas toujours =2. |
---|
354 | C |
---|
355 | C |
---|
356 | CALL grid_init(map_method, rst_type, n_srch_bins, |
---|
357 | $ src_size, dst_size, src_dims, dst_dims, |
---|
358 | $ src_rank, dst_rank, ncrn_src, ncrn_dst, |
---|
359 | $ sou_mask, tgt_mask, src_name, dst_name, |
---|
360 | $ src_lat, src_lon, dst_lat, dst_lon, |
---|
361 | $ src_corner_lat, src_corner_lon, |
---|
362 | $ dst_corner_lat, dst_corner_lon) |
---|
363 | |
---|
364 | C |
---|
365 | C* -- calculation of weights and addresses using SCRIP-library |
---|
366 | C |
---|
367 | CALL scrip(crmpfile, cmapping, map_method, normalize_opt, |
---|
368 | $ lextrapdone, rl_varmul, id_scripvoi) |
---|
369 | C |
---|
370 | IF (map_method == 'CONSERV') THEN |
---|
371 | DEALLOCATE(src_corner_lon, src_corner_lat, |
---|
372 | $ dst_corner_lon, dst_corner_lat) |
---|
373 | ENDIF |
---|
374 | C |
---|
375 | C* -- open just created SCRIP matrix file |
---|
376 | C |
---|
377 | CALL hdlerr(NF_OPEN |
---|
378 | $ (crmpfile, NF_NOWRITE, nc_scpid), 'scriprmp') |
---|
379 | C |
---|
380 | IF (nlogprt .GE. 2) THEN |
---|
381 | WRITE (UNIT = nulou,FMT = *) |
---|
382 | $ ' SCRIP file created and opened' |
---|
383 | WRITE (UNIT = nulou,FMT = *) ' ' |
---|
384 | CALL FLUSH(nulou) |
---|
385 | ENDIF |
---|
386 | ENDIF |
---|
387 | C |
---|
388 | C* 3. Read weights and addresses |
---|
389 | C -------------------------- |
---|
390 | C |
---|
391 | C* -- get matrix size |
---|
392 | C |
---|
393 | CALL hdlerr(NF_INQ_DIMID |
---|
394 | $ (nc_scpid, 'num_links', dimid), 'scriprmp') |
---|
395 | CALL hdlerr(NF_INQ_DIMLEN |
---|
396 | $ (nc_scpid, dimid, num_links), 'scriprmp') |
---|
397 | C |
---|
398 | C* -- get number of weights |
---|
399 | C |
---|
400 | SELECT CASE (map_method) |
---|
401 | CASE ('CONSERV') ! conservative remapping |
---|
402 | num_wgts = 3. |
---|
403 | CASE ('BILINEAR') ! bilinear remapping |
---|
404 | num_wgts = 1. |
---|
405 | CASE ('BICUBIC') ! bicubic remapping |
---|
406 | IF (cdgrdtyp .eq. 'D') THEN |
---|
407 | num_wgts = 1. |
---|
408 | ELSE |
---|
409 | num_wgts = 4. |
---|
410 | ENDIF |
---|
411 | CASE ('DISTWGT') ! distance weighted averaging |
---|
412 | num_wgts = 1. |
---|
413 | CASE ('GAUSWGT') ! distance gaussian weighted averaging |
---|
414 | num_wgts = 1. |
---|
415 | END SELECT |
---|
416 | C |
---|
417 | C* -- array allocation |
---|
418 | C |
---|
419 | ALLOCATE (src_addr(num_links), dst_addr(num_links), |
---|
420 | $ weights(num_wgts,num_links)) |
---|
421 | C |
---|
422 | C* -- read source grid addresses and weights |
---|
423 | C |
---|
424 | CALL hdlerr(NF_INQ_VARID |
---|
425 | $ (nc_scpid, csrcadd, varid), 'scriprmp') |
---|
426 | CALL hdlerr(NF_GET_VAR_INT |
---|
427 | $ (nc_scpid, varid, src_addr), 'scriprmp') |
---|
428 | CALL hdlerr(NF_INQ_VARID |
---|
429 | $ (nc_scpid, cdstadd, varid), 'scriprmp') |
---|
430 | CALL hdlerr(NF_GET_VAR_INT |
---|
431 | $ (nc_scpid, varid, dst_addr), 'scriprmp') |
---|
432 | CALL hdlerr(NF_INQ_VARID |
---|
433 | $ (nc_scpid, cweight, varid), 'scriprmp') |
---|
434 | CALL hdlerr(NF_GET_VAR_DOUBLE |
---|
435 | $ (nc_scpid, varid, weights), 'scriprmp') |
---|
436 | CALL hdlerr(NF_INQ_VARID |
---|
437 | $ (nc_scpid, cdstare, varid), 'scriprmp') |
---|
438 | CALL hdlerr(NF_GET_VAR_DOUBLE |
---|
439 | $ (nc_scpid, varid, dst_area), 'scriprmp') |
---|
440 | CALL hdlerr(NF_INQ_VARID |
---|
441 | $ (nc_scpid, cdstfra, varid), 'scriprmp') |
---|
442 | CALL hdlerr(NF_GET_VAR_DOUBLE |
---|
443 | $ (nc_scpid, varid, dst_frac), 'scriprmp') |
---|
444 | C |
---|
445 | C* 4. Do the matrix multiplication |
---|
446 | C ---------------------------- |
---|
447 | C |
---|
448 | SELECT CASE (map_method) |
---|
449 | C |
---|
450 | CASE ('CONSERV') ! conservative remapping |
---|
451 | C |
---|
452 | SELECT CASE (order) |
---|
453 | C |
---|
454 | CASE ('FIRST') ! first order remapping |
---|
455 | C |
---|
456 | DO n = 1, num_links |
---|
457 | IF (src_addr(n) .ne. 0) |
---|
458 | $ dst_array(dst_addr(n)) = dst_array(dst_addr(n)) |
---|
459 | $ + weights(1,n) * src_array(src_addr(n)) |
---|
460 | ENDDO |
---|
461 | C |
---|
462 | CASE ('SECOND') ! second order remapping |
---|
463 | C ! (including gradients) |
---|
464 | IF (cdgrdtyp .ne. 'LR') THEN |
---|
465 | WRITE (UNIT = nulou,FMT = *) |
---|
466 | $ 'Field gradient cannot be calculated' |
---|
467 | WRITE (UNIT = nulou,FMT = *) |
---|
468 | $ 'by Oasis as grid is not logically rectangular' |
---|
469 | CALL HALTE('STOP in scriprmp (CONSERV)') |
---|
470 | ENDIF |
---|
471 | ALLOCATE(gradient_lat(src_size), gradient_lon(src_size)) |
---|
472 | C |
---|
473 | call gradient(nlon_src, nlat_src, src_array, sou_mask, |
---|
474 | $ src_lat, src_lon, id_nosper, |
---|
475 | $ gradient_lat, gradient_lon) |
---|
476 | C |
---|
477 | DO n = 1, num_links |
---|
478 | IF (src_addr(n) .ne. 0) |
---|
479 | $ dst_array(dst_addr(n)) = dst_array(dst_addr(n)) |
---|
480 | $ + weights(1,n) * src_array(src_addr(n)) |
---|
481 | $ + weights(2,n) * gradient_lat(src_addr(n)) |
---|
482 | $ + weights(3,n) * gradient_lon(src_addr(n)) |
---|
483 | ENDDO |
---|
484 | DEALLOCATE(gradient_lat, gradient_lon) |
---|
485 | C |
---|
486 | END SELECT ! order |
---|
487 | C |
---|
488 | CASE ('BILINEAR') ! bilinear remapping |
---|
489 | C |
---|
490 | DO n = 1, num_links |
---|
491 | IF (src_addr(n) .ne. 0) |
---|
492 | $ dst_array(dst_addr(n)) = dst_array(dst_addr(n)) |
---|
493 | $ + weights(1,n) * src_array(src_addr(n)) |
---|
494 | ENDDO |
---|
495 | C |
---|
496 | CASE ('BICUBIC') ! bicubic remapping |
---|
497 | C |
---|
498 | SELECT CASE (cdgrdtyp) ! |
---|
499 | CASE ('LR') ! logically rectangular |
---|
500 | |
---|
501 | ALLOCATE(gradient_i(src_size), gradient_j(src_size), |
---|
502 | $ gradient_ij(src_size)) |
---|
503 | C |
---|
504 | CALL gradient_bicubic(nlon_src,nlat_src,src_array, |
---|
505 | $ sou_mask, src_lat, src_lon, id_nosper, |
---|
506 | $ gradient_i, gradient_j, gradient_ij) |
---|
507 | C |
---|
508 | DO n = 1, num_links |
---|
509 | IF (src_addr(n) .ne. 0) |
---|
510 | $ dst_array(dst_addr(n)) = dst_array(dst_addr(n)) |
---|
511 | $ + weights(1,n) * src_array(src_addr(n)) |
---|
512 | $ + weights(2,n) * gradient_i(src_addr(n)) |
---|
513 | $ + weights(3,n) * gradient_j(src_addr(n)) |
---|
514 | $ + weights(4,n) * gradient_ij(src_addr(n)) |
---|
515 | ENDDO |
---|
516 | C |
---|
517 | DEALLOCATE(gradient_i, gradient_j, gradient_ij) |
---|
518 | C |
---|
519 | CASE ('D') !reduced |
---|
520 | C |
---|
521 | DO n = 1, num_links |
---|
522 | IF (src_addr(n) .ne. 0) |
---|
523 | $ dst_array(dst_addr(n)) = dst_array(dst_addr(n)) |
---|
524 | $ + weights(1,n) * src_array(src_addr(n)) |
---|
525 | ENDDO |
---|
526 | C |
---|
527 | END SELECT |
---|
528 | C |
---|
529 | CASE ('DISTWGT') ! distance weighted average |
---|
530 | C |
---|
531 | DO n = 1, num_links |
---|
532 | IF (src_addr(n) .ne. 0) |
---|
533 | $ dst_array(dst_addr(n)) = dst_array(dst_addr(n)) |
---|
534 | $ + weights(1,n) * src_array(src_addr(n)) |
---|
535 | ENDDO |
---|
536 | C |
---|
537 | CASE ('GAUSWGT') ! distance gaussian weighted average |
---|
538 | C |
---|
539 | DO n = 1, num_links |
---|
540 | IF (src_addr(n) .ne. 0) |
---|
541 | $ dst_array(dst_addr(n)) = dst_array(dst_addr(n)) |
---|
542 | $ + weights(1,n) * src_array(src_addr(n)) |
---|
543 | ENDDO |
---|
544 | C |
---|
545 | END SELECT ! remapping method |
---|
546 | C |
---|
547 | C* 5. Check the cells on the poles |
---|
548 | C ---------------------------- |
---|
549 | C |
---|
550 | C* -- For the north pole _N and the south pole _S |
---|
551 | C |
---|
552 | latpol_N = pi*half |
---|
553 | latpol_S = -pi*half |
---|
554 | compt_N = 0 |
---|
555 | moy_tmp_N = 0d0 |
---|
556 | moy_ERR_N = 0d0 |
---|
557 | compt_S = 0 |
---|
558 | moy_tmp_S = 0d0 |
---|
559 | moy_ERR_S = 0d0 |
---|
560 | C |
---|
561 | DO n = 1, dst_size |
---|
562 | IF (dst_lat(n) == latpol_N) THEN |
---|
563 | moy_tmp_N = moy_tmp_N + dst_array(n) |
---|
564 | moy_err_N = moy_err_N + dst_err(n) |
---|
565 | compt_N = compt_N + 1 |
---|
566 | ELSE IF (dst_lat(n) == latpol_S) THEN |
---|
567 | moy_tmp_S = moy_tmp_S + dst_array(n) |
---|
568 | moy_err_S = moy_err_S + dst_err(n) |
---|
569 | compt_S = compt_S + 1 |
---|
570 | END IF |
---|
571 | END DO |
---|
572 | C |
---|
573 | IF (compt_N/=0) THEN |
---|
574 | moy_tmp_N = moy_tmp_N/compt_N |
---|
575 | moy_err_N = moy_err_N/compt_N |
---|
576 | DO n = 1, num_links |
---|
577 | IF (dst_lat(dst_addr(n)) == latpol_N) THEN |
---|
578 | dst_array(dst_addr(n)) = moy_tmp_N |
---|
579 | dst_err(dst_addr(n)) = moy_err_N |
---|
580 | END IF |
---|
581 | END DO |
---|
582 | PRINT*, 'Points at the pole N' |
---|
583 | PRINT*, 'Average value at the pole : ', moy_tmp_N |
---|
584 | PRINT*, 'Average error at the pole : ', moy_err_N |
---|
585 | ELSE |
---|
586 | PRINT*, 'No point at the pole N' |
---|
587 | END IF |
---|
588 | C |
---|
589 | IF (compt_S/=0) THEN |
---|
590 | moy_tmp_S = moy_tmp_S/compt_S |
---|
591 | moy_err_S = moy_err_S/compt_S |
---|
592 | DO n = 1, num_links |
---|
593 | IF (dst_lat(dst_addr(n)) == latpol_S) THEN |
---|
594 | dst_array(dst_addr(n)) = moy_tmp_S |
---|
595 | dst_err(dst_addr(n)) = moy_err_S |
---|
596 | END IF |
---|
597 | END DO |
---|
598 | PRINT*, 'Points at the pole S' |
---|
599 | PRINT*, 'Average value at the pole : ', moy_tmp_S |
---|
600 | PRINT*, 'Average error at the pole : ', moy_err_S |
---|
601 | ELSE |
---|
602 | PRINT*, 'No point at the pole S' |
---|
603 | END IF |
---|
604 | C |
---|
605 | DEALLOCATE (src_addr, dst_addr, weights) |
---|
606 | C |
---|
607 | C* 6. Close remapping file |
---|
608 | C -------------------- |
---|
609 | C |
---|
610 | CALL hdlerr(NF_CLOSE(nc_scpid), 'scriprmp') |
---|
611 | C |
---|
612 | C* 7. End of routine |
---|
613 | C -------------- |
---|
614 | C |
---|
615 | IF (nlogprt .GE. 2) THEN |
---|
616 | WRITE (UNIT = nulou,FMT = *) ' ' |
---|
617 | WRITE (UNIT = nulou,FMT = *) |
---|
618 | $ ' --------- End of routine scriprmp ---------' |
---|
619 | CALL FLUSH (nulou) |
---|
620 | ENDIF |
---|
621 | RETURN |
---|
622 | END |
---|