| blob | 89abd8eb04affc4fd178ed5d8c01c886e84bed33 |
1 MODULE module_interp_info
2 INTEGER , PARAMETER :: NOT_DEFINED_YET = 0
3 INTEGER , PARAMETER :: BILINEAR = 1
4 INTEGER , PARAMETER :: SINT = 2
5 INTEGER , PARAMETER :: NEAREST_NEIGHBOR = 3
6 INTEGER , PARAMETER :: QUADRATIC = 4
7 INTEGER , PARAMETER :: SPLINE = 5
8 INTEGER , PARAMETER :: SINT_NEW = 12
10 INTEGER :: interp_method_type = 0
11 CONTAINS
12 SUBROUTINE interp_info_init
13 #if (EM_CORE == 1)
14 CALL nl_get_interp_method_type ( 1 , interp_method_type )
15 #else
16 interp_method_type = 2
17 #endif
18 END SUBROUTINE interp_info_init
19 END MODULE module_interp_info
21 !WRF:MEDIATION_LAYER:INTERPOLATIONFUNCTION
22 !
23 !
25 !=========================================================================
27 SUBROUTINE interp_init
28 USE module_interp_info
29 CALL interp_info_init
30 END SUBROUTINE interp_init
32 !=========================================================================
34 SUBROUTINE interp_fcn ( cfld, & ! CD field
35 cids, cide, ckds, ckde, cjds, cjde, &
36 cims, cime, ckms, ckme, cjms, cjme, &
37 cits, cite, ckts, ckte, cjts, cjte, &
38 nfld, & ! ND field
39 nids, nide, nkds, nkde, njds, njde, &
40 nims, nime, nkms, nkme, njms, njme, &
41 nits, nite, nkts, nkte, njts, njte, &
42 shw, & ! stencil half width for interp
43 imask, & ! interpolation mask
44 xstag, ystag, & ! staggering of field
45 ipos, jpos, & ! Position of lower left of nest in CD
46 nri, nrj ) ! Nest ratio, i- and j-directions
48 USE module_interp_info
50 IMPLICIT NONE
52 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
53 cims, cime, ckms, ckme, cjms, cjme, &
54 cits, cite, ckts, ckte, cjts, cjte, &
55 nids, nide, nkds, nkde, njds, njde, &
56 nims, nime, nkms, nkme, njms, njme, &
57 nits, nite, nkts, nkte, njts, njte, &
58 shw, &
59 ipos, jpos, &
60 nri, nrj
61 LOGICAL, INTENT(IN) :: xstag, ystag
63 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
64 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
65 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
67 IF ( interp_method_type .EQ. NOT_DEFINED_YET ) THEN
68 interp_method_type = SINT
69 END IF
71 IF ( interp_method_type .EQ. BILINEAR ) THEN
72 CALL interp_fcn_blint ( cfld, & ! CD field
73 cids, cide, ckds, ckde, cjds, cjde, &
74 cims, cime, ckms, ckme, cjms, cjme, &
75 cits, cite, ckts, ckte, cjts, cjte, &
76 nfld, & ! ND field
77 nids, nide, nkds, nkde, njds, njde, &
78 nims, nime, nkms, nkme, njms, njme, &
79 nits, nite, nkts, nkte, njts, njte, &
80 shw, & ! stencil half width for interp
81 imask, & ! interpolation mask
82 xstag, ystag, & ! staggering of field
83 ipos, jpos, & ! Position of lower left of nest in CD
84 nri, nrj) ! Nest ratio, i- and j-directions
85 ELSE IF ( MOD(interp_method_type,10) .EQ. SINT ) THEN
86 CALL interp_fcn_sint ( cfld, & ! CD field
87 cids, cide, ckds, ckde, cjds, cjde, &
88 cims, cime, ckms, ckme, cjms, cjme, &
89 cits, cite, ckts, ckte, cjts, cjte, &
90 nfld, & ! ND field
91 nids, nide, nkds, nkde, njds, njde, &
92 nims, nime, nkms, nkme, njms, njme, &
93 nits, nite, nkts, nkte, njts, njte, &
94 shw, & ! stencil half width for interp
95 imask, & ! interpolation mask
96 xstag, ystag, & ! staggering of field
97 ipos, jpos, & ! Position of lower left of nest in CD
98 nri, nrj) ! Nest ratio, i- and j-directions
99 ELSE IF ( interp_method_type .EQ. NEAREST_NEIGHBOR ) THEN
100 CALL interp_fcn_nn ( cfld, & ! CD field
101 cids, cide, ckds, ckde, cjds, cjde, &
102 cims, cime, ckms, ckme, cjms, cjme, &
103 cits, cite, ckts, ckte, cjts, cjte, &
104 nfld, & ! ND field
105 nids, nide, nkds, nkde, njds, njde, &
106 nims, nime, nkms, nkme, njms, njme, &
107 nits, nite, nkts, nkte, njts, njte, &
108 shw, & ! stencil half width for interp
109 imask, & ! interpolation mask
110 xstag, ystag, & ! staggering of field
111 ipos, jpos, & ! Position of lower left of nest in CD
112 nri, nrj) ! Nest ratio, i- and j-directions
113 ELSE IF ( interp_method_type .EQ. QUADRATIC ) THEN
114 CALL interp_fcn_lagr ( cfld, & ! CD field
115 cids, cide, ckds, ckde, cjds, cjde, &
116 cims, cime, ckms, ckme, cjms, cjme, &
117 cits, cite, ckts, ckte, cjts, cjte, &
118 nfld, & ! ND field
119 nids, nide, nkds, nkde, njds, njde, &
120 nims, nime, nkms, nkme, njms, njme, &
121 nits, nite, nkts, nkte, njts, njte, &
122 shw, & ! stencil half width for interp
123 imask, & ! interpolation mask
124 xstag, ystag, & ! staggering of field
125 ipos, jpos, & ! Position of lower left of nest in CD
126 nri, nrj) ! Nest ratio, i- and j-directions
127 ELSE
128 CALL wrf_error_fatal ('Hold on there cowboy, we need to know which interpolation option you want')
129 END IF
131 END SUBROUTINE interp_fcn
133 !=========================================================================
135 ! Overlapping linear horizontal iterpolation for mass, u, and v staggerings.
137 SUBROUTINE interp_fcn_blint ( cfld, & ! CD field
138 cids, cide, ckds, ckde, cjds, cjde, &
139 cims, cime, ckms, ckme, cjms, cjme, &
140 cits, cite, ckts, ckte, cjts, cjte, &
141 nfld, & ! ND field
142 nids, nide, nkds, nkde, njds, njde, &
143 nims, nime, nkms, nkme, njms, njme, &
144 nits, nite, nkts, nkte, njts, njte, &
145 shw, & ! stencil half width for interp
146 imask, & ! interpolation mask
147 xstag, ystag, & ! staggering of field
148 ipos, jpos, & ! Position of lower left of nest in CD
149 nri, nrj) ! Nest ratio, i- and j-directions
151 IMPLICIT NONE
153 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
154 cims, cime, ckms, ckme, cjms, cjme, &
155 cits, cite, ckts, ckte, cjts, cjte, &
156 nids, nide, nkds, nkde, njds, njde, &
157 nims, nime, nkms, nkme, njms, njme, &
158 nits, nite, nkts, nkte, njts, njte, &
159 shw, &
160 ipos, jpos, &
161 nri, nrj
162 LOGICAL, INTENT(IN) :: xstag, ystag
164 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
165 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
166 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
168 ! Local
170 INTEGER ci, cj, ck, ni, nj, nk, istag, jstag, ioff, joff, i, j, k
171 REAL :: wx, wy, cfld_ll, cfld_lr, cfld_ul, cfld_ur
172 REAL :: cxp0, cxp1, nx, cyp0, cyp1, ny
174 ! Fortran functions. Yes, yes, I know, probably pretty slow.
176 REAL, EXTERNAL :: nest_loc_of_cg
177 INTEGER, EXTERNAL :: compute_CGLL
179 ! This stag stuff is to keep us away from the outer most row
180 ! and column for the unstaggered directions. We are going to
181 ! consider "U" an xstag variable and "V" a ystag variable. The
182 ! vertical staggering is handled in the actual arguments. The
183 ! ckte and nkte are the ending vertical dimensions for computations
184 ! for this particular variable.
186 IF ( xstag ) THEN
187 istag = 0
188 ioff = 1
189 ELSE
190 istag = 1
191 ioff = 0
192 END IF
194 IF ( ystag ) THEN
195 jstag = 0
196 joff = 1
197 ELSE
198 jstag = 1
199 joff = 0
200 END IF
202 ! Loop over each j-index on this tile for the nested domain.
204 j_loop : DO nj = njts, MIN(njde-jstag,njte)
206 ! This is the lower-left j-index of the CG.
208 ! Example is 3:1 ratio, mass-point staggering. We have listed six CG values
209 ! as an example: A, B, C, D, E, F. For a 3:1 ratio, each of these CG cells has
210 ! nine associated FG points.
211 ! |=========|=========|=========|
212 ! | - - - | - - - | - - - |
213 ! | | | |
214 ! | - D - | - E - | - F - |
215 ! | | | |
216 ! | 1 2 3 | 4 5 6 | 7 8 9 |
217 ! |=========|=========|=========|
218 ! | - - - | - - - | - - - |
219 ! | | | |
220 ! | - A - | - B - | - C - |
221 ! | | | |
222 ! | - - - | - - - | - - - |
223 ! |=========|=========|=========|
224 ! To interpolate to FG point 4, we will use CG points: A, B, D, E. It is adequate to
225 ! find the lower left point. The lower left (LL) point for "4" is "A". Below
226 ! are a few more points.
227 ! 2 => A
228 ! 3 => A
229 ! 4 => A
230 ! 5 => B
231 ! 6 => B
232 ! 7 => B
234 cj = compute_CGLL ( nj , jpos , nrj , jstag )
235 ny = REAL(nj)
236 cyp0 = nest_loc_of_cg ( cj , jpos , nrj , joff )
237 cyp1 = nest_loc_of_cg ( cj+1 , jpos , nrj , joff )
239 ! What is the weighting for this CG point to the FG point, j-weight only.
241 wy = ( cyp1 - ny ) / ( cyp1 - cyp0 )
243 ! Vertical dim of the nest domain.
245 k_loop : DO nk = nkts, nkte
247 ! Loop over each i-index on this tile for the nested domain.
249 i_loop : DO ni = nits, MIN(nide-istag,nite)
251 IF ( imask ( ni, nj ) .EQ. 1 ) THEN
253 ! The coarse grid location that is to the lower left of the FG point.
255 ci = compute_CGLL ( ni , ipos , nri , istag )
256 nx = REAL(ni)
257 cxp0 = nest_loc_of_cg ( ci , ipos , nri , ioff )
258 cxp1 = nest_loc_of_cg ( ci+1 , ipos , nri , ioff )
260 wx = ( cxp1 - nx ) / ( cxp1 - cxp0 )
262 ! The four surrounding CG values.
264 cfld_ll = cfld(ci ,nk,cj )
265 cfld_lr = cfld(ci+1,nk,cj )
266 cfld_ul = cfld(ci ,nk,cj+1)
267 cfld_ur = cfld(ci+1,nk,cj+1)
269 ! Bilinear interpolation in horizontal.
271 nfld( ni , nk , nj ) = wy * ( cfld_ll * wx + cfld_lr * (1.-wx) ) + &
272 (1.-wy) * ( cfld_ul * wx + cfld_ur * (1.-wx) )
274 END IF
275 END DO i_loop
276 END DO k_loop
277 END DO j_loop
279 END SUBROUTINE interp_fcn_blint
281 !=========================================================================
283 ! Overlapping linear horizontal iterpolation for longitude
285 SUBROUTINE interp_fcn_blint_ll ( cfld_inp, & ! CD field
286 cids, cide, ckds, ckde, cjds, cjde, &
287 cims, cime, ckms, ckme, cjms, cjme, &
288 cits, cite, ckts, ckte, cjts, cjte, &
289 nfld, & ! ND field
290 nids, nide, nkds, nkde, njds, njde, &
291 nims, nime, nkms, nkme, njms, njme, &
292 nits, nite, nkts, nkte, njts, njte, &
293 shw, & ! stencil half width for interp
294 imask, & ! interpolation mask
295 xstag, ystag, & ! staggering of field
296 ipos, jpos, & ! Position of lower left of nest in CD
297 nri, nrj, & ! Nest ratio, i- and j-directions
298 clat_in, nlat_in, & ! CG, FG latitude
299 cinput_from_file, ninput_from_file ) ! CG, FG T/F input from file
301 IMPLICIT NONE
303 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
304 cims, cime, ckms, ckme, cjms, cjme, &
305 cits, cite, ckts, ckte, cjts, cjte, &
306 nids, nide, nkds, nkde, njds, njde, &
307 nims, nime, nkms, nkme, njms, njme, &
308 nits, nite, nkts, nkte, njts, njte, &
309 shw, &
310 ipos, jpos, &
311 nri, nrj
312 LOGICAL, INTENT(IN) :: xstag, ystag
314 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld_inp, cfld
315 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
316 REAL, DIMENSION ( cims:cime, cjms:cjme ) :: clat_in
317 REAL, DIMENSION ( nims:nime, njms:njme ) :: nlat_in
318 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
319 LOGICAL :: cinput_from_file, ninput_from_file
321 ! Local
323 INTEGER ci, cj, ck, ni, nj, nk, istag, jstag, ioff, joff, i, j, k
324 REAL :: wx, wy, cfld_ll, cfld_lr, cfld_ul, cfld_ur
325 REAL :: cxp0, cxp1, nx, cyp0, cyp1, ny
326 LOGICAL :: probably_by_dateline
327 REAL :: max_lon, min_lon
328 LOGICAL :: probably_by_pole
329 REAL :: max_lat, min_lat
331 ! Fortran functions. Yes, yes, I know, probably pretty slow.
333 REAL, EXTERNAL :: nest_loc_of_cg
334 INTEGER, EXTERNAL :: compute_CGLL
336 ! This stag stuff is to keep us away from the outer most row
337 ! and column for the unstaggered directions. We are going to
338 ! consider "U" an xstag variable and "V" a ystag variable. The
339 ! vertical staggering is handled in the actual arguments. The
340 ! ckte and nkte are the ending vertical dimensions for computations
341 ! for this particular variable.
343 IF ( xstag ) THEN
344 istag = 0
345 ioff = 1
346 ELSE
347 istag = 1
348 ioff = 0
349 END IF
351 IF ( ystag ) THEN
352 jstag = 0
353 joff = 1
354 ELSE
355 jstag = 1
356 joff = 0
357 END IF
359 ! If this is a projection where the nest is over the pole, and
360 ! we are using the parent to interpolate the longitudes, then
361 ! we are going to have longitude troubles. If this is the case,
362 ! stop the model right away.
364 probably_by_pole = .FALSE.
365 max_lat = -90
366 min_lat = +90
367 DO nj = njts, MIN(njde-jstag,njte)
368 DO ni = nits, MIN(nide-istag,nite)
369 max_lat = MAX ( nlat_in(ni,nj) , max_lat )
370 min_lat = MIN ( nlat_in(ni,nj) , min_lat )
371 END DO
372 END DO
374 IF ( ( max_lat .GT. 85 ) .OR. ( ABS(min_lat) .GT. 85 ) ) THEN
375 probably_by_pole = .TRUE.
376 END IF
378 IF ( ( probably_by_pole ) .AND. ( .NOT. ninput_from_file ) ) THEN
379 CALL wrf_error_fatal ( 'Nest over the pole, single input domain, longitudes will be wrong' )
380 END IF
382 ! Initialize to NOT being by dateline.
384 probably_by_dateline = .FALSE.
385 max_lon = -180
386 min_lon = +180
387 DO nj = njts, MIN(njde-jstag,njte)
388 cj = compute_CGLL ( nj , jpos , nrj , jstag )
389 DO ni = nits, MIN(nide-istag,nite)
390 ci = compute_CGLL ( ni , ipos , nri , istag )
391 max_lon = MAX ( cfld_inp(ci,1,cj) , max_lon )
392 min_lon = MIN ( cfld_inp(ci,1,cj) , min_lon )
393 END DO
394 END DO
396 IF ( max_lon - min_lon .GT. 300 ) THEN
397 probably_by_dateline = .TRUE.
398 END IF
400 ! Load "continuous" longitude across the date line
402 DO cj = MIN(cjts-1,cjms), MAX(cjte+1,cjme)
403 DO ci = MIN(cits-1,cims), MAX(cite+1,cime)
404 IF ( ( cfld_inp(ci,ckts,cj) .LT. 0 ) .AND. ( probably_by_dateline ) ) THEN
405 cfld(ci,ckts,cj) = 360 + cfld_inp(ci,ckts,cj)
406 ELSE
407 cfld(ci,ckts,cj) = cfld_inp(ci,ckts,cj)
408 END IF
409 END DO
410 END DO
412 ! Loop over each j-index on this tile for the nested domain.
414 j_loop : DO nj = njts, MIN(njde-jstag,njte)
416 ! This is the lower-left j-index of the CG.
418 ! Example is 3:1 ratio, mass-point staggering. We have listed six CG values
419 ! as an example: A, B, C, D, E, F. For a 3:1 ratio, each of these CG cells has
420 ! nine associated FG points.
421 ! |=========|=========|=========|
422 ! | - - - | - - - | - - - |
423 ! | | | |
424 ! | - D - | - E - | - F - |
425 ! | | | |
426 ! | 1 2 3 | 4 5 6 | 7 8 9 |
427 ! |=========|=========|=========|
428 ! | - - - | - - - | - - - |
429 ! | | | |
430 ! | - A - | - B - | - C - |
431 ! | | | |
432 ! | - - - | - - - | - - - |
433 ! |=========|=========|=========|
434 ! To interpolate to FG point 4, we will use CG points: A, B, D, E. It is adequate to
435 ! find the lower left point. The lower left (LL) point for "4" is "A". Below
436 ! are a few more points.
437 ! 2 => A
438 ! 3 => A
439 ! 4 => A
440 ! 5 => B
441 ! 6 => B
442 ! 7 => B
444 cj = compute_CGLL ( nj , jpos , nrj , jstag )
445 ny = REAL(nj)
446 cyp0 = nest_loc_of_cg ( cj , jpos , nrj , joff )
447 cyp1 = nest_loc_of_cg ( cj+1 , jpos , nrj , joff )
449 ! What is the weighting for this CG point to the FG point, j-weight only.
451 wy = ( cyp1 - ny ) / ( cyp1 - cyp0 )
453 ! Vertical dim of the nest domain.
455 k_loop : DO nk = nkts, nkte
457 ! Loop over each i-index on this tile for the nested domain.
459 i_loop : DO ni = nits, MIN(nide-istag,nite)
461 IF ( imask ( ni, nj ) .EQ. 1 ) THEN
463 ! The coarse grid location that is to the lower left of the FG point.
465 ci = compute_CGLL ( ni , ipos , nri , istag )
466 nx = REAL(ni)
467 cxp0 = nest_loc_of_cg ( ci , ipos , nri , ioff )
468 cxp1 = nest_loc_of_cg ( ci+1 , ipos , nri , ioff )
470 wx = ( cxp1 - nx ) / ( cxp1 - cxp0 )
472 ! The four surrounding CG values.
474 cfld_ll = cfld(ci ,nk,cj )
475 cfld_lr = cfld(ci+1,nk,cj )
476 cfld_ul = cfld(ci ,nk,cj+1)
477 cfld_ur = cfld(ci+1,nk,cj+1)
479 ! Bilinear interpolation in horizontal.
481 nfld( ni , nk , nj ) = wy * ( cfld_ll * wx + cfld_lr * (1.-wx) ) + &
482 (1.-wy) * ( cfld_ul * wx + cfld_ur * (1.-wx) )
484 END IF
485 END DO i_loop
486 END DO k_loop
487 END DO j_loop
489 ! Put nested longitude back into the -180 to 180 range.
491 DO nj = njts, MIN(njde-jstag,njte)
492 DO ni = nits, MIN(nide-istag,nite)
493 IF ( nfld(ni,nkts,nj) .GT. 180 ) THEN
494 nfld(ni,nkts,nj) = -360 + nfld(ni,nkts,nj)
495 END IF
496 END DO
497 END DO
499 END SUBROUTINE interp_fcn_blint_ll
501 !=========================================================================
503 ! Lagrange interpolating polynomials, set up as a quadratic, with an average of
504 ! the overlap.
506 SUBROUTINE interp_fcn_lagr ( cfld, & ! CD field
507 cids, cide, ckds, ckde, cjds, cjde, &
508 cims, cime, ckms, ckme, cjms, cjme, &
509 cits, cite, ckts, ckte, cjts, cjte, &
510 nfld, & ! ND field
511 nids, nide, nkds, nkde, njds, njde, &
512 nims, nime, nkms, nkme, njms, njme, &
513 nits, nite, nkts, nkte, njts, njte, &
514 shw, & ! stencil half width for interp
515 imask, & ! interpolation mask
516 xstag, ystag, & ! staggering of field
517 ipos, jpos, & ! Position of lower left of nest in CD
518 nri, nrj) ! Nest ratio, i- and j-directions
520 IMPLICIT NONE
522 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
523 cims, cime, ckms, ckme, cjms, cjme, &
524 cits, cite, ckts, ckte, cjts, cjte, &
525 nids, nide, nkds, nkde, njds, njde, &
526 nims, nime, nkms, nkme, njms, njme, &
527 nits, nite, nkts, nkte, njts, njte, &
528 shw, &
529 ipos, jpos, &
530 nri, nrj
531 LOGICAL, INTENT(IN) :: xstag, ystag
533 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
534 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
535 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
537 ! Local
539 INTEGER ci, cj, ck, ni, nj, nk, istag, jstag, i, j, k
540 REAL :: nx, x0, x1, x2, x3, x
541 REAL :: ny, y0, y1, y2, y3
542 REAL :: cxm1, cxp0, cxp1, cxp2, nfld_m1, nfld_p0, nfld_p1, nfld_p2
543 REAL :: cym1, cyp0, cyp1, cyp2
544 INTEGER :: ioff, joff
546 ! Fortran functions.
548 REAL, EXTERNAL :: lagrange_quad_avg
549 REAL, EXTERNAL :: nest_loc_of_cg
550 INTEGER, EXTERNAL :: compute_CGLL
552 ! This stag stuff is to keep us away from the outer most row
553 ! and column for the unstaggered directions. We are going to
554 ! consider "U" an xstag variable and "V" a ystag variable. The
555 ! vertical staggering is handled in the actual arguments. The
556 ! ckte and nkte are the ending vertical dimensions for computations
557 ! for this particular variable.
559 ! The ioff and joff are offsets due to the staggering. It is a lot
560 ! simpler with ioff and joff if
561 ! u var => ioff=1
562 ! v var => joff=1
563 ! otherwise zero.
564 ! Note that is OPPOSITE of the istag, jstag vars. The stag variables are
565 ! used for the domain dimensions, the offset guys are used in the
566 ! determination of grid points between the CG and FG
568 IF ( xstag ) THEN
569 istag = 0
570 ioff = 1
571 ELSE
572 istag = 1
573 ioff = 0
574 END IF
576 IF ( ystag ) THEN
577 jstag = 0
578 joff = 1
579 ELSE
580 jstag = 1
581 joff = 0
582 END IF
584 ! Loop over each j-index on this tile for the nested domain.
586 j_loop : DO nj = njts, MIN(njde-jstag,njte)
588 ! This is the lower-left j-index of the CG.
590 ! Example is 3:1 ratio, mass-point staggering. We have listed sixteen CG values
591 ! as an example: A through P. For a 3:1 ratio, each of these CG cells has
592 ! nine associated FG points.
594 ! |=========|=========|=========|=========|
595 ! | - - - | - - - | - - - | - - - |
596 ! | | | | |
597 ! | - M - | - N d | - O - | - P - |
598 ! | | | | |
599 ! | - - - | - - - | - - - | - - - |
600 ! |=========|=========|=========|=========|
601 ! | - - - | - - - | - - - | - - - |
602 ! | | | | |
603 ! | - I - | - J c | - K - | - L - |
604 ! | | | | |
605 ! | - - - | - - - | - - - | - - - |
606 ! |=========|=========|=========|=========|
607 ! | - 1 2 | 3 4 5 | 6 7 8 | - - - |
608 ! | | | | |
609 ! | - E - | - F b | - G - | - H - |
610 ! | | | | |
611 ! | - - - | - - - | - - - | - - - |
612 ! |=========|=========|=========|=========|
613 ! | - - - | - - - | - - - | - - - |
614 ! | | | | |
615 ! | - A - | - B a | - C - | - D - |
616 ! | | | | |
617 ! | - - - | - - - | - - - | - - - |
618 ! |=========|=========|=========|=========|
620 ! To interpolate to FG point 4, 5, or 6 we will use CG points: A through P. It is
621 ! sufficient to find the lower left corner of a 4-point interpolation, and then extend
622 ! each side by one unit.
624 ! Here are the lower left hand corners of the following FG points:
625 ! 1 => E
626 ! 2 => E
627 ! 3 => E
628 ! 4 => F
629 ! 5 => F
630 ! 6 => F
631 ! 7 => G
632 ! 8 => G
634 cj = compute_CGLL ( nj , jpos , nrj , jstag )
636 ! Vertical dim of the nest domain.
638 k_loop : DO nk = nkts, nkte
640 ! Loop over each i-index on this tile for the nested domain.
642 i_loop : DO ni = nits, MIN(nide-istag,nite)
644 ! The coarse grid location that is to the lower left of the FG point.
646 ci = compute_CGLL ( ni , ipos , nri , istag )
648 ! To interpolate to point "*" (look in grid cell "F"):
649 ! 1. Use ABC to get a quadratic valid at "a"
650 ! Use BCD to get a quadratic valid at "a"
651 ! Average these to get the final value for "a"
652 ! 2. Use EFG to get a quadratic valid at "b"
653 ! Use FGH to get a quadratic valid at "b"
654 ! Average these to get the final value for "b"
655 ! 3. Use IJK to get a quadratic valid at "c"
656 ! Use JKL to get a quadratic valid at "c"
657 ! Average these to get the final value for "c"
658 ! 4. Use MNO to get a quadratic valid at "d"
659 ! Use NOP to get a quadratic valid at "d"
660 ! Average these to get the final value for "d"
661 ! 5. Use abc to get a quadratic valid at "*"
662 ! Use bcd to get a quadratic valid at "*"
663 ! Average these to get the final value for "*"
665 ! |=========|=========|=========|=========|
666 ! | - - - | - - - | - - - | - - - |
667 ! | | | | |
668 ! | - M - | - N d | - O - | - P - |
669 ! | | | | |
670 ! | - - - | - - - | - - - | - - - |
671 ! |=========|=========|=========|=========|
672 ! | - - - | - - - | - - - | - - - |
673 ! | | | | |
674 ! | - I - | - J c | - K - | - L - |
675 ! | | | | |
676 ! | - - - | - - - | - - - | - - - |
677 ! |=========|=========|=========|=========|
678 ! | - - - | - - * | - - - | - - - |
679 ! | | | | |
680 ! | - E - | - F b | - G - | - H - |
681 ! | | | | |
682 ! | - - - | - - - | - - - | - - - |
683 ! |=========|=========|=========|=========|
684 ! | - - - | - - - | - - - | - - - |
685 ! | | | | |
686 ! | - A - | - B a | - C - | - D - |
687 ! | | | | |
688 ! | - - - | - - - | - - - | - - - |
689 ! |=========|=========|=========|=========|
691 ! Overlapping quadratic interpolation.
693 IF ( imask ( ni, nj ) .EQ. 1 ) THEN
695 ! I-direction location of "*"
697 nx = REAL(ni)
699 ! I-direction location of "A", "E", "I", "M"
701 cxm1 = nest_loc_of_cg ( ci-1 , ipos , nri , ioff )
703 ! I-direction location of "B", "F", "J", "N"
705 cxp0 = nest_loc_of_cg ( ci , ipos , nri , ioff )
707 ! I-direction location of "C", "G", "K", "O"
709 cxp1 = nest_loc_of_cg ( ci+1 , ipos , nri , ioff )
711 ! I-direction location of "D", "H", "L", "P"
713 cxp2 = nest_loc_of_cg ( ci+2 , ipos , nri , ioff )
715 ! Value at "a"
717 nfld_m1 = lagrange_quad_avg ( nx, cxm1, cxp0, cxp1, cxp2, cfld(ci-1,nk,cj-1), cfld(ci+0,nk,cj-1), cfld(ci+1,nk,cj-1), cfld(ci+2,nk,cj-1) )
719 ! Value at "b"
721 nfld_p0 = lagrange_quad_avg ( nx, cxm1, cxp0, cxp1, cxp2, cfld(ci-1,nk,cj+0), cfld(ci+0,nk,cj+0), cfld(ci+1,nk,cj+0), cfld(ci+2,nk,cj+0) )
723 ! Value at "c"
725 nfld_p1 = lagrange_quad_avg ( nx, cxm1, cxp0, cxp1, cxp2, cfld(ci-1,nk,cj+1), cfld(ci+0,nk,cj+1), cfld(ci+1,nk,cj+1), cfld(ci+2,nk,cj+1) )
727 ! Value at "d"
729 nfld_p2 = lagrange_quad_avg ( nx, cxm1, cxp0, cxp1, cxp2, cfld(ci-1,nk,cj+2), cfld(ci+0,nk,cj+2), cfld(ci+1,nk,cj+2), cfld(ci+2,nk,cj+2) )
731 ! J-direction location of "*"
733 ny = REAL(nj)
735 ! J-direction location of "A", "B", "C", "D"
737 cym1 = nest_loc_of_cg ( cj-1 , jpos , nrj , joff )
739 ! J-direction location of "E", "F", "G", "H"
741 cyp0 = nest_loc_of_cg ( cj , jpos , nrj , joff )
743 ! J-direction location of "I", "J", "K", "L"
745 cyp1 = nest_loc_of_cg ( cj+1 , jpos , nrj , joff )
747 ! J-direction location of "M", "N", "O", "P"
749 cyp2 = nest_loc_of_cg ( cj+2 , jpos , nrj , joff )
751 ! Value at "*"
753 nfld(ni,nk,nj) = lagrange_quad_avg ( ny, cym1, cyp0, cyp1, &
754 cyp2, nfld_m1, nfld_p0, nfld_p1, nfld_p2 )
756 END IF
758 END DO i_loop
759 END DO k_loop
760 END DO j_loop
762 END SUBROUTINE interp_fcn_lagr
764 !=================================================================================
766 REAL FUNCTION lagrange_quad ( x , x0, x1, x2, y0, y1, y2 )
768 IMPLICIT NONE
770 REAL :: x , x0, x1, x2, y0, y1, y2
772 ! Lagrange = sum prod ( x - xj )
773 ! i=0,n ( j=0,n --------- * yi )
774 ! j<>i ( xi - xj )
776 ! For a quadratic, in the above equation, we are setting n=2. Three points
777 ! required for a quadratic, points x0, x1, x2 (hence n=2).
779 lagrange_quad = &
780 (x-x1)*(x-x2)*y0 / ( (x0-x1)*(x0-x2) ) + &
781 (x-x0)*(x-x2)*y1 / ( (x1-x0)*(x1-x2) ) + &
782 (x-x0)*(x-x1)*y2 / ( (x2-x0)*(x2-x1) )
784 END FUNCTION lagrange_quad
786 !=================================================================================
788 REAL FUNCTION lagrange_quad_avg ( x , x0, x1, x2, x3, y0, y1, y2, y3 )
790 IMPLICIT NONE
792 REAL, EXTERNAL :: lagrange_quad
793 REAL :: x , x0, x1, x2, x3, y0, y1, y2, y3
795 ! Since there are three points required for a quadratic, we compute it twice
796 ! (once with x0, x1, x2 and once with x1, x2, x3), and then average those values. This will
797 ! reduce overshoot. The "x" point is where we are interpolating TO.
799 ! x0 x1 x x2
800 ! x1 x x2 x3
802 lagrange_quad_avg = &
803 ! ( lagrange_quad ( x , x0, x1, x2, y0, y1, y2 ) + &
804 ! lagrange_quad ( x , x1, x2, x3, y1, y2, y3 ) ) / &
805 ! 2.
806 ( lagrange_quad ( x , x0, x1, x2, y0, y1, y2 ) * ( x2 - x ) + &
807 lagrange_quad ( x , x1, x2, x3, y1, y2, y3 ) * ( x - x1 ) ) / &
808 ( x2 - x1 )
810 END FUNCTION lagrange_quad_avg
812 !=================================================================================
814 REAL FUNCTION nest_loc_of_cg ( ci , ipos , nri , ioff )
816 ! I and J direction equations for mass and momentum values for even
817 ! and odd ratios: Given that the starting value of the nest in the
818 ! CG grid cell is defined as (1,1), what is the location of the CG
819 ! location in FG index units. Example, for a 2:1 ratio, the location
820 ! of the mass point T is 1.5 (3:1 ratio = 2, 4:1 ratio = 2.5, etc).
821 ! Note that for momentum points, the CG U point is defined as "1", the
822 ! same as the I-direction of the (1,1) location of the FG U point.
823 ! Same for V, but in the J-direction.
825 IMPLICIT NONE
827 INTEGER :: ci , ipos , nri , ioff
829 nest_loc_of_cg = &
830 ( ci - ipos ) * nri + ( 1 - ioff ) * REAL ( nri + 1 ) / 2. + ioff
832 END FUNCTION nest_loc_of_cg
834 !=================================================================================
836 FUNCTION compute_CGLL ( ni , ipos , nri , istag ) RESULT ( CGLL_loc )
838 IMPLICIT NONE
840 INTEGER , INTENT(IN ) :: ni , ipos , nri , istag
841 INTEGER :: CGLL_loc
843 ! Local vars
845 INTEGER :: starting_position , increments_of_CG_cells
846 INTEGER :: location_of_LL_wrt_this_CG
847 INTEGER :: ioff
848 INTEGER , PARAMETER :: MOMENTUM_STAG = 0
849 INTEGER , PARAMETER :: MASS_POINT_STAG = 1
851 starting_position = ipos
852 increments_of_CG_cells = ( ni - 1 ) / nri
853 ioff = MOD ( nri , 2 )
855 IF ( istag .EQ. MOMENTUM_STAG ) THEN
856 location_of_LL_wrt_this_CG = MOD ( ( ni - 1 ) , nri ) / ( nri + ioff ) - istag ! zero
857 ELSE IF ( istag .EQ. MASS_POINT_STAG ) THEN
858 location_of_LL_wrt_this_CG = ( MOD ( ( ni - 1 ) , nri ) + ioff ) / ( ( nri + ioff ) / 2 ) - istag
859 ELSE
860 CALL wrf_error_fatal ( 'Hold on there pard, there are only two staggerings I accept.' )
861 END IF
863 CGLL_loc = starting_position + increments_of_CG_cells + location_of_LL_wrt_this_CG
864 ! WRITE ( 6 , '(a,i4, i4, i4, i4)') 'ni ipos nri stag', ni, ipos, nri, istag
865 ! WRITE ( 6 , '(a,i4, i4, i4, i4)') 'strt inc loc CGLL', starting_position , increments_of_CG_cells , location_of_LL_wrt_this_CG , CGLL_loc
866 ! print *,' '
868 END FUNCTION compute_CGLL
870 !=================================================================================
872 ! Smolarkiewicz positive definite, monotonic transport.
874 SUBROUTINE interp_fcn_sint ( cfld, & ! CD field
875 cids, cide, ckds, ckde, cjds, cjde, &
876 cims, cime, ckms, ckme, cjms, cjme, &
877 cits, cite, ckts, ckte, cjts, cjte, &
878 nfld, & ! ND field
879 nids, nide, nkds, nkde, njds, njde, &
880 nims, nime, nkms, nkme, njms, njme, &
881 nits, nite, nkts, nkte, njts, njte, &
882 shw, & ! stencil half width for interp
883 imask, & ! interpolation mask
884 xstag, ystag, & ! staggering of field
885 ipos, jpos, & ! Position of lower left of nest in CD
886 nri, nrj ) ! nest ratios
887 USE module_configure
889 IMPLICIT NONE
891 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
892 cims, cime, ckms, ckme, cjms, cjme, &
893 cits, cite, ckts, ckte, cjts, cjte, &
894 nids, nide, nkds, nkde, njds, njde, &
895 nims, nime, nkms, nkme, njms, njme, &
896 nits, nite, nkts, nkte, njts, njte, &
897 shw, &
898 ipos, jpos, &
899 nri, nrj
900 LOGICAL, INTENT(IN) :: xstag, ystag
902 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
903 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
904 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
906 ! Local
908 INTEGER ci, cj, ck, ni, nj, nk, ip, jp, ioff, joff, nioff, njoff
909 INTEGER nfx, ior
910 PARAMETER (ior=2)
911 INTEGER nf
912 REAL psca(cims:cime,cjms:cjme,nri*nrj)
913 LOGICAL icmask( cims:cime, cjms:cjme )
914 INTEGER i,j,k
915 INTEGER nrio2, nrjo2
917 ! Iterate over the ND tile and compute the values
918 ! from the CD tile.
920 ioff = 0 ; joff = 0
921 nioff = 0 ; njoff = 0
922 IF ( xstag ) THEN
923 ioff = (nri-1)/2
924 nioff = nri
925 ENDIF
926 IF ( ystag ) THEN
927 joff = (nrj-1)/2
928 njoff = nrj
929 ENDIF
931 nrio2 = nri/2
932 nrjo2 = nrj/2
934 nfx = nri * nrj
935 !$OMP PARALLEL DO &
936 !$OMP PRIVATE ( i,j,k,ni,nj,ci,cj,ip,jp,nk,ck,nf,icmask,psca )
937 DO k = ckts, ckte
938 icmask = .FALSE.
939 DO nf = 1,nfx
940 DO j = cjms,cjme
941 nj = (j-jpos) * nrj + ( nrjo2 + 1 ) ! j point on nest
942 DO i = cims,cime
943 ni = (i-ipos) * nri + ( nrio2 + 1 ) ! i point on nest
944 if ( ni .ge. nits-nioff-nrio2 .and. &
945 ni .le. nite+nioff+nrio2 .and. &
946 nj .ge. njts-njoff-nrjo2 .and. &
947 nj .le. njte+njoff+nrjo2 ) then
948 if ( ni.ge.nims.and.ni.le.nime.and.nj.ge.njms.and.nj.le.njme) then
949 if ( imask(ni,nj) .eq. 1 ) then
950 icmask( i, j ) = .TRUE.
951 endif
952 endif
953 if ( ni-nioff.ge.nims.and.ni.le.nime.and.nj-njoff.ge.njms.and.nj.le.njme) then
954 if (ni .ge. nits-nioff .and. nj .ge. njts-njoff ) then
955 if ( imask(ni-nioff,nj-njoff) .eq. 1) then
956 icmask( i, j ) = .TRUE.
957 endif
958 endif
959 endif
960 endif
961 psca(i,j,nf) = cfld(i,k,j)
962 ENDDO ! i
963 ENDDO ! j
964 ENDDO ! nf
966 ! tile dims in this call to sint are 1-over to account for the fact
967 ! that the number of cells on the nest local subdomain is not
968 ! necessarily a multiple of the nest ratio in a given dim.
969 ! this could be a little less ham-handed.
971 CALL sint( psca, &
972 cims, cime, cjms, cjme, icmask, &
973 cits-1, cite+1, cjts-1, cjte+1, nrj*nri, xstag, ystag )
975 DO nj = njts, njte+joff
976 cj = jpos + (nj-1) / nrj ! j coord of CD point
977 jp = mod ( nj-1 , nrj ) ! coord of ND w/i CD point
978 nk = k
979 ck = nk
980 DO ni = nits, nite+ioff
981 ci = ipos + (ni-1) / nri ! i coord of CD point
982 ip = mod ( ni-1 , nri ) ! coord of ND w/i CD point
983 if ( ( ni-ioff .ge. nits ) .and. ( nj-joff .ge. njts ) ) then
984 if ( imask ( ni, nj ) .eq. 1 .or. imask ( ni-ioff, nj-joff ) .eq. 1 ) then
985 nfld( ni-ioff, nk, nj-joff ) = psca( ci , cj, ip+1 + (jp)*nri )
986 endif
987 endif
988 ENDDO
989 ENDDO
990 ENDDO
991 !$OMP END PARALLEL DO
993 END SUBROUTINE interp_fcn_sint
995 !=========================================================================
997 ! Nearest neighbor interpolation.
999 SUBROUTINE interp_fcn_nn ( cfld, & ! CD field
1000 cids, cide, ckds, ckde, cjds, cjde, &
1001 cims, cime, ckms, ckme, cjms, cjme, &
1002 cits, cite, ckts, ckte, cjts, cjte, &
1003 nfld, & ! ND field
1004 nids, nide, nkds, nkde, njds, njde, &
1005 nims, nime, nkms, nkme, njms, njme, &
1006 nits, nite, nkts, nkte, njts, njte, &
1007 shw, & ! stencil half width for interp
1008 imask, & ! interpolation mask
1009 xstag, ystag, & ! staggering of field
1010 ipos, jpos, & ! Position of lower left of nest in CD
1011 nri, nrj ) ! nest ratios
1012 IMPLICIT NONE
1014 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
1015 cims, cime, ckms, ckme, cjms, cjme, &
1016 cits, cite, ckts, ckte, cjts, cjte, &
1017 nids, nide, nkds, nkde, njds, njde, &
1018 nims, nime, nkms, nkme, njms, njme, &
1019 nits, nite, nkts, nkte, njts, njte, &
1020 shw, &
1021 ipos, jpos, &
1022 nri, nrj
1023 LOGICAL, INTENT(IN) :: xstag, ystag
1025 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
1026 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
1027 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
1029 INTEGER ci, cj, ck, ni, nj, nk
1031 ! Iterate over the ND tile and assign the values
1032 ! from the CD tile. This is a trivial implementation
1033 ! of the interp_fcn; just copies the values from the CD into the ND
1035 DO nj = njts, njte
1036 cj = jpos + (nj-1) / nrj ! j coord of CD point
1037 DO nk = nkts, nkte
1038 ck = nk
1039 DO ni = nits, nite
1040 if ( imask ( ni, nj ) .eq. 1 ) then
1041 ci = ipos + (ni-1) / nri ! i coord of CD point
1042 nfld( ni, nk, nj ) = cfld( ci , ck , cj )
1043 endif
1044 ENDDO
1045 ENDDO
1046 ENDDO
1048 END SUBROUTINE interp_fcn_nn
1050 !=========================================================================
1052 SUBROUTINE interp_fcn_bl ( cfld, & ! CD field
1053 cids, cide, ckds, ckde, cjds, cjde, &
1054 cims, cime, ckms, ckme, cjms, cjme, &
1055 cits, cite, ckts, ckte, cjts, cjte, &
1056 nfld, & ! ND field
1057 nids, nide, nkds, nkde, njds, njde, &
1058 nims, nime, nkms, nkme, njms, njme, &
1059 nits, nite, nkts, nkte, njts, njte, &
1060 shw, & ! stencil half width for interp
1061 imask, & ! interpolation mask
1062 xstag, ystag, & ! staggering of field
1063 ipos, jpos, & ! Position of lower left of nest in CD
1064 nri, nrj, & ! Nest ratio, i- and j-directions
1065 cht, nht, & ! topography for CG and FG
1066 ct_max_p,nt_max_p, & ! temperature (K) at max press, want CG value
1067 cght_max_p,nght_max_p, & ! height (m) at max press, want CG value
1068 cmax_p,nmax_p, & ! max pressure (Pa) in column, want CG value
1069 ct_min_p,nt_min_p, & ! temperature (K) at min press, want CG value
1070 cght_min_p,nght_min_p, & ! height (m) at min press, want CG value
1071 cmin_p,nmin_p, & ! min pressure (Pa) in column, want CG value
1072 zn, p_top ) ! eta levels
1073 USE module_timing
1074 ! USE module_configure
1075 USE module_model_constants , ONLY : g , r_d, cp, p1000mb, t0
1077 IMPLICIT NONE
1080 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
1081 cims, cime, ckms, ckme, cjms, cjme, &
1082 cits, cite, ckts, ckte, cjts, cjte, &
1083 nids, nide, nkds, nkde, njds, njde, &
1084 nims, nime, nkms, nkme, njms, njme, &
1085 nits, nite, nkts, nkte, njts, njte, &
1086 shw, &
1087 ipos, jpos, &
1088 nri, nrj
1089 LOGICAL, INTENT(IN) :: xstag, ystag
1091 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
1092 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
1093 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
1094 REAL, DIMENSION ( cims:cime, cjms:cjme ) :: cht, ct_max_p, cght_max_p, cmax_p, ct_min_p, cght_min_p, cmin_p
1095 REAL, DIMENSION ( nims:nime, njms:njme ) :: nht, nt_max_p, nght_max_p, nmax_p, nt_min_p, nght_min_p, nmin_p
1096 REAL, DIMENSION ( ckms:ckme ) :: zn
1097 REAL :: p_top
1098 REAL, EXTERNAL :: v_interp_col
1100 ! Local
1102 INTEGER ci, cj, ni, nj, nk, istag, jstag, i, j, k
1103 REAL :: wx, wy, nprs, cfld_ll, cfld_lr, cfld_ul, cfld_ur
1104 REAL , DIMENSION(ckms:ckme) :: cprs
1105 REAL :: p00 , t00 , a , tiso , p_surf
1107 ! Yes, memory sized to allow "outside the tile" indexing for horiz interpolation. This
1108 ! is really an intermediate domain that has quite a bit of usable real estate surrounding
1109 ! the tile dimensions.
1111 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cpb
1113 ! A bit larger than tile sized to allow horizontal interpolation on the CG.
1115 REAL, DIMENSION ( cits-2:cite+2, cjts-2:cjte+2 ) :: cfld_max_p, cfld_min_p
1117 ! The usual tile size for the FG local array.
1119 REAL, DIMENSION ( nits:nite, nkts:nkte, njts:njte ) :: npb
1121 ! Get base state constants
1123 CALL nl_get_base_pres ( 1 , p00 )
1124 CALL nl_get_base_temp ( 1 , t00 )
1125 CALL nl_get_base_lapse ( 1 , a )
1126 CALL nl_get_iso_temp ( 1 , tiso )
1128 ! This stag stuff is to keep us away from the outer most row
1129 ! and column for the unstaggered directions. We are going to
1130 ! consider "U" an xstag variable and "V" a ystag variable. The
1131 ! vertical staggering is handled in the actual arguments. The
1132 ! ckte and nkte are the ending vertical dimensions for computations
1133 ! for this particular variable.
1135 IF ( xstag ) THEN
1136 istag = 0
1137 ELSE
1138 istag = 1
1139 END IF
1141 IF ( ystag ) THEN
1142 jstag = 0
1143 ELSE
1144 jstag = 1
1145 END IF
1147 ! Compute the reference pressure for the CG, function only of constants and elevation.
1148 ! We extend the i,j range to allow us to do horizontal interpolation. We only need
1149 ! one extra grid cell surrounding the nest, and the intermediate domain has plenty of
1150 ! room with the halos set up for higher-order interpolations. For intermediate domains,
1151 ! it turns out that the "domain" size actually fits within the "tile" size. Yeppers,
1152 ! that is backwards from what usually happens. That intermediate domain size is a couple
1153 ! grid points larger than necessary, and the tile is a couple of grid cells larger still.
1154 ! For our low-order interpolation, we can use the tile size for the CG, and we will have
1155 ! plenty of data on our boundaries.
1157 DO j = cjts-2 , cjte+2
1158 DO i = cits-2 , cite+2
1159 p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*cht(i,j)/a/r_d ) **0.5 )
1160 DO k = ckts , ckte
1161 cpb(i,k,j) = zn(k)*(p_surf - p_top) + p_top
1162 END DO
1163 IF ( ckte .EQ. ckme ) THEN
1164 cfld_max_p(i,j) = cght_max_p(i,j) * g
1165 cfld_min_p(i,j) = cght_min_p(i,j) * g
1166 ELSE
1167 cfld_max_p(i,j) = ct_max_p(i,j) * (p1000mb/cmax_p(i,j))**(r_d/cp) - t0
1168 cfld_min_p(i,j) = ct_min_p(i,j) * (p1000mb/cmin_p(i,j))**(r_d/cp) - t0
1169 END IF
1170 END DO
1171 END DO
1173 ! Compute the reference pressure for the FG. This is actually the size of the entire
1174 ! domain, not some chopped down piece of intermediate domain, as in the parent
1175 ! grid. We do the traditional MAX(dom end -1,tile end), since we know a priori that the
1176 ! pressure is a mass point field (because the topo elevation is a mass point field).
1178 DO j = njts , MIN(njde-1,njte)
1179 DO i = nits , MIN(nide-1,nite)
1180 p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*nht(i,j)/a/r_d ) **0.5 )
1181 DO k = nkts , nkte
1182 npb(i,k,j) = zn(k)*(p_surf - p_top) + p_top
1183 END DO
1184 END DO
1185 END DO
1187 ! Loop over each j-index on this tile for the nested domain.
1189 j_loop : DO nj = njts, MIN(njde-jstag,njte)
1191 ! This is the lower-left j-index of the CG.
1193 ! Example is 3:1 ratio, mass-point staggering. We have listed six CG values
1194 ! as an example: A, B, C, D, E, F. For a 3:1 ratio, each of these CG cells has
1195 ! nine associated FG points.
1196 ! |=========|=========|=========|
1197 ! | - - - | - - - | - - - |
1198 ! | | | |
1199 ! | - D - | - E - | - F - |
1200 ! | | | |
1201 ! | 1 2 3 | 4 5 6 | 7 8 9 |
1202 ! |=========|=========|=========|
1203 ! | - - - | - - - | - - - |
1204 ! | | | |
1205 ! | - A - | - B - | - C - |
1206 ! | | | |
1207 ! | - - - | - - - | - - - |
1208 ! |=========|=========|=========|
1209 ! To interpolate to FG point 4, we will use CG points: A, B, D, E. It is adequate to
1210 ! find the lower left point. The lower left (LL) point for "4" is "A". Below
1211 ! are a few more points.
1212 ! 2 => A
1213 ! 3 => A
1214 ! 4 => A
1215 ! 5 => B
1216 ! 6 => B
1217 ! 7 => B
1218 ! We want an equation that returns the CG LL:
1219 ! CG LL = ipos (the starting point of the nest in the CG)
1220 ! + (ni-1)/nri (gives us the CG cell, based on the nri-groups of FG cells
1221 ! - istag (a correction term, this is either zero for u in the x-dir,
1222 ! since we are doing an "i" example, or 1 for anything else)
1223 ! + (MOD(ni-1,nri)+1 + nri/2)/nri (gives us specifically related CG point for each of the nri
1224 ! FG points, for example, we want points "1", "4", and "7" all
1225 ! to point to the CG at the left for the LL point)
1226 ! For grid points 4, 5, 6, we want the CG LL (sans the first two terms) to be -1, 0, 0 (which
1227 ! means that the CG point for "4" is to the left, and the CG LL point for "5" and "6"
1228 ! is in the current CG index.
1230 cj = jpos + (nj-1)/nrj - jstag + (MOD(nj-1,nrj)+1 + nrj/2)/nrj
1233 ! What is the weighting for this CG point to the FG point, j-weight only.
1235 IF ( ystag ) THEN
1236 wy = 1. - ( REAL(MOD(nj+(nrj-1)/2,nrj)) / REAL(nrj) + 1. / REAL (2 * nrj) )
1237 ELSE
1238 wy = 1. - ( REAL(MOD(nj+(nrj-1)/2,nrj)) / REAL(nrj) )
1239 END IF
1241 ! Vertical dim of the nest domain.
1243 k_loop : DO nk = nkts, nkte
1245 ! Loop over each i-index on this tile for the nested domain.
1247 i_loop : DO ni = nits, MIN(nide-istag,nite)
1249 ! The coarse grid location that is to the lower left of the FG point.
1251 ci = ipos + (ni-1)/nri - istag + (MOD(ni-1,nri)+1 + nri/2)/nri
1253 ! Weights in the x-direction.
1255 IF ( xstag ) THEN
1256 wx = 1. - ( REAL(MOD(ni+(nri-1)/2,nri)) / REAL(nri) + 1. / REAL (2 * nri) )
1257 ELSE
1258 wx = 1. - ( REAL(MOD(ni+(nri-1)/2,nri)) / REAL(nri) )
1259 END IF
1261 ! The pressure of the FG point.
1263 IF ( ( .NOT. xstag ) .AND. ( .NOT. ystag ) ) THEN
1264 nprs = npb( ni , nk , nj )
1265 ELSE IF ( xstag ) THEN
1266 nprs = ( npb( ni-1, nk , nj ) + npb( ni , nk , nj ) ) * 0.5
1267 ELSE IF ( ystag ) THEN
1268 nprs = ( npb( ni , nk , nj-1) + npb( ni , nk , nj ) ) * 0.5
1269 END IF
1271 ! The four surrounding CG values.
1273 IF ( ( .NOT. xstag ) .AND. ( .NOT. ystag ) ) THEN
1274 cprs(:) = cpb(ci ,:,cj )
1275 cfld_ll = v_interp_col ( cfld(ci ,:,cj ) , cprs(:) , ckms , ckme , ckte, nprs, ni, nj, nk, ci , cj , &
1276 cfld_max_p(ci ,cj ) , cmax_p(ci ,cj ) , cfld_min_p(ci ,cj ) , cmin_p(ci ,cj ) )
1277 cprs(:) = cpb(ci+1,:,cj )
1278 cfld_lr = v_interp_col ( cfld(ci+1,:,cj ) , cprs(:) , ckms , ckme , ckte, nprs, ni, nj, nk, ci+1, cj , &
1279 cfld_max_p(ci+1,cj ) , cmax_p(ci+1,cj ) , cfld_min_p(ci+1,cj ) , cmin_p(ci+1,cj ) )
1280 cprs(:) = cpb(ci ,:,cj+1)
1281 cfld_ul = v_interp_col ( cfld(ci ,:,cj+1) , cprs(:) , ckms , ckme , ckte, nprs, ni, nj, nk, ci , cj+1, &
1282 cfld_max_p(ci ,cj+1) , cmax_p(ci ,cj+1) , cfld_min_p(ci ,cj+1) , cmin_p(ci ,cj+1) )
1283 cprs(:) = cpb(ci+1,:,cj+1)
1284 cfld_ur = v_interp_col ( cfld(ci+1,:,cj+1) , cprs(:) , ckms , ckme , ckte, nprs, ni, nj, nk, ci+1, cj+1, &
1285 cfld_max_p(ci+1,cj+1) , cmax_p(ci+1,cj+1) , cfld_min_p(ci+1,cj+1) , cmin_p(ci+1,cj+1) )
1287 ELSE IF ( xstag ) THEN
1288 cprs(:) = ( cpb(ci ,:,cj ) + cpb(ci-1,:,cj ) )*0.5
1289 cfld_ll = v_interp_col ( cfld(ci ,:,cj ) , cprs(:) , ckms , ckme , ckte, nprs, ni, nj, nk, ci , cj , &
1290 cfld_max_p(ci ,cj ) , cmax_p(ci ,cj ) , cfld_min_p(ci ,cj ) , cmin_p(ci ,cj ) )
1291 cprs(:) = ( cpb(ci+1,:,cj ) + cpb(ci ,:,cj ) )*0.5
1292 cfld_lr = v_interp_col ( cfld(ci+1,:,cj ) , cprs(:) , ckms , ckme , ckte, nprs, ni, nj, nk, ci+1, cj , &
1293 cfld_max_p(ci+1,cj ) , cmax_p(ci+1,cj ) , cfld_min_p(ci+1,cj ) , cmin_p(ci+1,cj ) )
1294 cprs(:) = ( cpb(ci ,:,cj+1) + cpb(ci-1,:,cj+1) )*0.5
1295 cfld_ul = v_interp_col ( cfld(ci ,:,cj+1) , cprs(:) , ckms , ckme , ckte, nprs, ni, nj, nk, ci , cj+1, &
1296 cfld_max_p(ci ,cj+1) , cmax_p(ci ,cj+1) , cfld_min_p(ci ,cj+1) , cmin_p(ci ,cj+1) )
1297 cprs(:) = ( cpb(ci+1,:,cj+1) + cpb(ci ,:,cj+1) )*0.5
1298 cfld_ur = v_interp_col ( cfld(ci+1,:,cj+1) , cprs(:) , ckms , ckme , ckte, nprs, ni, nj, nk, ci+1, cj+1, &
1299 cfld_max_p(ci+1,cj+1) , cmax_p(ci+1,cj+1) , cfld_min_p(ci+1,cj+1) , cmin_p(ci+1,cj+1) )
1300 ELSE IF ( ystag ) THEN
1301 cprs(:) = ( cpb(ci ,:,cj ) + cpb(ci ,:,cj-1) )*0.5
1302 cfld_ll = v_interp_col ( cfld(ci ,:,cj ) , cprs(:) , ckms , ckme , ckte, nprs, ni, nj, nk, ci , cj , &
1303 cfld_max_p(ci ,cj ) , cmax_p(ci ,cj ) , cfld_min_p(ci ,cj ) , cmin_p(ci ,cj ) )
1304 cprs(:) = ( cpb(ci+1,:,cj ) + cpb(ci+1,:,cj-1) )*0.5
1305 cfld_lr = v_interp_col ( cfld(ci+1,:,cj ) , cprs(:) , ckms , ckme , ckte, nprs, ni, nj, nk, ci+1, cj , &
1306 cfld_max_p(ci+1,cj ) , cmax_p(ci+1,cj ) , cfld_min_p(ci+1,cj ) , cmin_p(ci+1,cj ) )
1307 cprs(:) = ( cpb(ci ,:,cj+1) + cpb(ci ,:,cj ) )*0.5
1308 cfld_ul = v_interp_col ( cfld(ci ,:,cj+1) , cprs(:) , ckms , ckme , ckte, nprs, ni, nj, nk, ci , cj+1, &
1309 cfld_max_p(ci ,cj+1) , cmax_p(ci ,cj+1) , cfld_min_p(ci ,cj+1) , cmin_p(ci ,cj+1) )
1310 cprs(:) = ( cpb(ci+1,:,cj+1) + cpb(ci+1,:,cj ) )*0.5
1311 cfld_ur = v_interp_col ( cfld(ci+1,:,cj+1) , cprs(:) , ckms , ckme , ckte, nprs, ni, nj, nk, ci+1, cj+1, &
1312 cfld_max_p(ci+1,cj+1) , cmax_p(ci+1,cj+1) , cfld_min_p(ci+1,cj+1) , cmin_p(ci+1,cj+1) )
1313 END IF
1315 ! Bilinear interpolation in horizontal with vertically corrected CG field values.
1317 nfld( ni , nk , nj ) = wy * ( cfld_ll * wx + cfld_lr * (1.-wx) ) + &
1318 (1.-wy) * ( cfld_ul * wx + cfld_ur * (1.-wx) )
1320 END DO i_loop
1321 END DO k_loop
1322 END DO j_loop
1324 ! If this is ph_2, make the values at k=1 all zero
1326 IF ( ckme .EQ. ckte ) THEN
1327 DO nj = njts,njte
1328 DO ni = nits, nite
1329 nfld(ni,nkts,nj) = 0.0
1330 END DO
1331 END DO
1332 END IF
1334 END SUBROUTINE interp_fcn_bl
1336 !==================================
1338 FUNCTION v_interp_col ( cfld_orig , cprs_orig , ckms , ckme , ckte , nprs, ni, nj, nk, ci, cj, &
1339 cfld_max_p , cmax_p , cfld_min_p , cmin_p ) RESULT ( cfld_interp )
1341 IMPLICIT NONE
1343 INTEGER , INTENT(IN) :: ni, nj, nk, ci, cj
1344 INTEGER , INTENT(IN) :: ckms , ckme , ckte
1345 REAL , DIMENSION(ckms:ckme) , INTENT(IN) :: cfld_orig , cprs_orig
1346 REAL , INTENT(IN) :: cfld_max_p , cmax_p , cfld_min_p , cmin_p
1347 REAL , INTENT(IN) :: nprs
1348 REAL :: cfld_interp
1350 ! Local
1352 INTEGER :: ck
1353 LOGICAL :: found
1354 CHARACTER(LEN=256) :: joe_mess
1355 REAL , DIMENSION(ckms:ckme+1+1) :: cfld , cprs
1357 ! Fill input arrays
1359 cfld(1) = cfld_max_p
1360 cprs(1) = cmax_p
1362 cfld(ckte+2) = cfld_min_p
1363 cprs(ckte+2) = cmin_p
1365 DO ck = ckms , ckte
1366 cfld(ck+1) = cfld_orig(ck)
1367 cprs(ck+1) = cprs_orig(ck)
1368 END DO
1370 found = .FALSE.
1372 IF ( cprs(ckms) .LT. nprs ) THEN
1373 cfld_interp = cfld(ckms)
1374 RETURN
1375 ELSE IF ( cprs(ckte+2) .GE. nprs ) THEN
1376 cfld_interp = cfld(ckte+2)
1377 RETURN
1378 END IF
1380 DO ck = ckms , ckte+1
1381 IF ( ( cprs(ck ) .GE. nprs ) .AND. &
1382 ( cprs(ck+1) .LT. nprs ) ) THEN
1383 cfld_interp = ( cfld(ck ) * ( nprs - cprs(ck+1) ) + &
1384 cfld(ck+1) * ( cprs(ck) - nprs ) ) / &
1385 ( cprs(ck) - cprs(ck+1) )
1386 RETURN
1387 END IF
1388 END DO
1390 CALL wrf_error_fatal ( 'ERROR -- vertical interpolation for nest interp cannot find trapping pressures' )
1392 END FUNCTION v_interp_col
1394 !==================================
1395 ! this is the default function used in feedback.
1397 SUBROUTINE copy_fcn ( cfld, & ! CD field
1398 cids, cide, ckds, ckde, cjds, cjde, &
1399 cims, cime, ckms, ckme, cjms, cjme, &
1400 cits, cite, ckts, ckte, cjts, cjte, &
1401 nfld, & ! ND field
1402 nids, nide, nkds, nkde, njds, njde, &
1403 nims, nime, nkms, nkme, njms, njme, &
1404 nits, nite, nkts, nkte, njts, njte, &
1405 shw, & ! stencil half width for interp
1406 imask, & ! interpolation mask
1407 xstag, ystag, & ! staggering of field
1408 ipos, jpos, & ! Position of lower left of nest in CD
1409 nri, nrj ) ! nest ratios
1410 USE module_configure
1411 IMPLICIT NONE
1414 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
1415 cims, cime, ckms, ckme, cjms, cjme, &
1416 cits, cite, ckts, ckte, cjts, cjte, &
1417 nids, nide, nkds, nkde, njds, njde, &
1418 nims, nime, nkms, nkme, njms, njme, &
1419 nits, nite, nkts, nkte, njts, njte, &
1420 shw, &
1421 ipos, jpos, &
1422 nri, nrj
1423 LOGICAL, INTENT(IN) :: xstag, ystag
1425 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ), INTENT(OUT) :: cfld
1426 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ),INTENT(IN) :: nfld
1427 INTEGER, DIMENSION ( nims:nime, njms:njme ),INTENT(IN) :: imask
1429 ! Local
1431 INTEGER ci, cj, ck, ni, nj, nk, ip, jp, ioff, joff, ioffa, joffa
1432 INTEGER :: icmin,icmax,jcmin,jcmax
1433 INTEGER :: istag,jstag, ipoints,jpoints,ijpoints
1434 INTEGER , PARAMETER :: passes = 2
1435 INTEGER spec_zone
1437 ! Loop over the coarse grid in the area of the fine mesh. Do not
1438 ! process the coarse grid values that are along the lateral BC
1439 ! provided to the fine grid. Since that is in the specified zone
1440 ! for the fine grid, it should not be used in any feedback to the
1441 ! coarse grid as it should not have changed.
1443 ! Due to peculiarities of staggering, it is simpler to handle the feedback
1444 ! for the staggerings based upon whether it is a even ratio (2::1, 4::1, etc.) or
1445 ! an odd staggering ratio (3::1, 5::1, etc.).
1447 ! Though there are separate grid ratios for the i and j directions, this code
1448 ! is not general enough to handle aspect ratios .NE. 1 for the fine grid cell.
1450 ! These are local integer increments in the looping. Basically, istag=1 means
1451 ! that we will assume one less point in the i direction. Note that ci and cj
1452 ! have a maximum value that is decreased by istag and jstag, respectively.
1454 ! Horizontal momentum feedback is along the face, not within the cell. For a
1455 ! 3::1 ratio, temperature would use 9 pts for feedback, while u and v use
1456 ! only 3 points for feedback from the nest to the parent.
1458 CALL nl_get_spec_zone( 1 , spec_zone )
1459 istag = 1 ; jstag = 1
1460 IF ( xstag ) istag = 0
1461 IF ( ystag ) jstag = 0
1463 IF( MOD(nrj,2) .NE. 0) THEN ! odd refinement ratio
1465 IF ( ( .NOT. xstag ) .AND. ( .NOT. ystag ) ) THEN
1466 DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
1467 nj = (cj-jpos)*nrj + nrj/2 + 1
1468 DO ck = ckts, ckte
1469 nk = ck
1470 DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
1471 ni = (ci-ipos)*nri + nri/2 + 1
1472 cfld( ci, ck, cj ) = 0.
1473 DO ijpoints = 1 , nri * nrj
1474 ipoints = MOD((ijpoints-1),nri) + 1 - nri/2 - 1
1475 jpoints = (ijpoints-1)/nri + 1 - nrj/2 - 1
1476 cfld( ci, ck, cj ) = cfld( ci, ck, cj ) + &
1477 1./REAL(nri*nrj) * nfld( ni+ipoints , nk , nj+jpoints )
1478 END DO
1479 ! cfld( ci, ck, cj ) = 1./9. * &
1480 ! ( nfld( ni-1, nk , nj-1) + &
1481 ! nfld( ni , nk , nj-1) + &
1482 ! nfld( ni+1, nk , nj-1) + &
1483 ! nfld( ni-1, nk , nj ) + &
1484 ! nfld( ni , nk , nj ) + &
1485 ! nfld( ni+1, nk , nj ) + &
1486 ! nfld( ni-1, nk , nj+1) + &
1487 ! nfld( ni , nk , nj+1) + &
1488 ! nfld( ni+1, nk , nj+1) )
1489 ! cfld( ci, ck, cj ) = 1./25. * &
1490 ! ( nfld( ni-2, nk , nj-2) + &
1491 ! nfld( ni-2, nk , nj-1) + &
1492 ! nfld( ni-2, nk , nj-0) + &
1493 ! nfld( ni-2, nk , nj+1) + &
1494 ! nfld( ni-2, nk , nj+2) + &
1495 ! nfld( ni-1, nk , nj-2) + &
1496 ! nfld( ni-1, nk , nj-1) + &
1497 ! nfld( ni-1, nk , nj-0) + &
1498 ! nfld( ni-1, nk , nj+1) + &
1499 ! nfld( ni-1, nk , nj+2) + &
1500 ! nfld( ni-0, nk , nj-2) + &
1501 ! nfld( ni-0, nk , nj-1) + &
1502 ! nfld( ni-0, nk , nj-0) + &
1503 ! nfld( ni-0, nk , nj+1) + &
1504 ! nfld( ni-0, nk , nj+2) + &
1505 ! nfld( ni+1, nk , nj-2) + &
1506 ! nfld( ni+1, nk , nj-1) + &
1507 ! nfld( ni+1, nk , nj-0) + &
1508 ! nfld( ni+1, nk , nj+1) + &
1509 ! nfld( ni+1, nk , nj+2) + &
1510 ! nfld( ni+2, nk , nj-2) + &
1511 ! nfld( ni+2, nk , nj-1) + &
1512 ! nfld( ni+2, nk , nj-0) + &
1513 ! nfld( ni+2, nk , nj+1) + &
1514 ! nfld( ni+2, nk , nj+2) )
1515 ENDDO
1516 ENDDO
1517 ENDDO
1519 ELSE IF ( ( xstag ) .AND. ( .NOT. ystag ) ) THEN
1520 DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
1521 nj = (cj-jpos)*nrj + nrj/2 + 1
1522 DO ck = ckts, ckte
1523 nk = ck
1524 DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
1525 ni = (ci-ipos)*nri + 1
1526 cfld( ci, ck, cj ) = 0.
1527 DO ijpoints = (nri+1)/2 , (nri+1)/2 + nri*(nri-1) , nri
1528 ipoints = MOD((ijpoints-1),nri) + 1 - nri/2 - 1
1529 jpoints = (ijpoints-1)/nri + 1 - nrj/2 - 1
1530 cfld( ci, ck, cj ) = cfld( ci, ck, cj ) + &
1531 1./REAL(nri ) * nfld( ni+ipoints , nk , nj+jpoints )
1532 END DO
1533 ! cfld( ci, ck, cj ) = 1./3. * &
1534 ! ( nfld( ni , nk , nj-1) + &
1535 ! nfld( ni , nk , nj ) + &
1536 ! nfld( ni , nk , nj+1) )
1537 ENDDO
1538 ENDDO
1539 ENDDO
1541 ELSE IF ( ( .NOT. xstag ) .AND. ( ystag ) ) THEN
1542 DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
1543 nj = (cj-jpos)*nrj + 1
1544 DO ck = ckts, ckte
1545 nk = ck
1546 DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
1547 ni = (ci-ipos)*nri + nri/2 + 1
1548 cfld( ci, ck, cj ) = 0.
1549 DO ijpoints = ( nrj*nrj +1 )/2 - nrj/2 , ( nrj*nrj +1 )/2 - nrj/2 + nrj-1
1550 ipoints = MOD((ijpoints-1),nri) + 1 - nri/2 - 1
1551 jpoints = (ijpoints-1)/nri + 1 - nrj/2 - 1
1552 cfld( ci, ck, cj ) = cfld( ci, ck, cj ) + &
1553 1./REAL( nrj) * nfld( ni+ipoints , nk , nj+jpoints )
1554 END DO
1555 ! cfld( ci, ck, cj ) = 1./3. * &
1556 ! ( nfld( ni-1, nk , nj ) + &
1557 ! nfld( ni , nk , nj ) + &
1558 ! nfld( ni+1, nk , nj ) )
1559 ENDDO
1560 ENDDO
1561 ENDDO
1563 END IF
1565 ! Even refinement ratio
1567 ELSE IF ( MOD(nrj,2) .EQ. 0) THEN
1568 IF ( ( .NOT. xstag ) .AND. ( .NOT. ystag ) ) THEN
1570 ! This is a simple schematic of the feedback indexing used in the even
1571 ! ratio nests. For simplicity, a 2::1 ratio is depicted. Only the
1572 ! mass variable staggering is shown.
1573 ! Each of
1574 ! the boxes with a "T" and four small "t" represents a coarse grid (CG)
1575 ! cell, that is composed of four (2::1 ratio) fine grid (FG) cells.
1577 ! Shown below is the area of the CG that is in the area of the FG. The
1578 ! first grid point of the depicted CG is the starting location of the nest
1579 ! in the parent domain (ipos,jpos = i_parent_start and j_parent_start from
1580 ! the namelist).
1582 ! For each of the CG points, the feedback loop is over each of the FG points
1583 ! within the CG cell. For a 2::1 ratio, there are four total points (this is
1584 ! the ijpoints loop). The feedback value to the CG is the arithmetic mean of
1585 ! all of the FG values within each CG cell.
1587 ! |-------------||-------------| |-------------||-------------|
1588 ! | t t || t t | | t t || t t |
1589 ! jpos+ | || | | || |
1590 ! (njde-njds)- | T || T | | T || T |
1591 ! jstag | || | | || |
1592 ! | t t || t t | | t t || t t |
1593 ! |-------------||-------------| |-------------||-------------|
1594 ! |-------------||-------------| |-------------||-------------|
1595 ! | t t || t t | | t t || t t |
1596 ! | || | | || |
1597 ! | T || T | | T || T |
1598 ! | || | | || |
1599 ! | t t || t t | | t t || t t |
1600 ! |-------------||-------------| |-------------||-------------|
1601 !
1602 ! ...
1603 ! ...
1604 ! ...
1605 ! ...
1606 ! ...
1608 ! |-------------||-------------| |-------------||-------------|
1609 ! jpoints = 1 | t t || t t | | t t || t t |
1610 ! | || | | || |
1611 ! | T || T | | T || T |
1612 ! | || | | || |
1613 ! jpoints = 0, | t t || t t | | t t || t t |
1614 ! nj=3 |-------------||-------------| |-------------||-------------|
1615 ! |-------------||-------------| |-------------||-------------|
1616 ! jpoints = 1 | t t || t t | | t t || t t |
1617 ! | || | | || |
1618 ! jpos --> | T || T | ... | T || T |
1619 ! | || | ... | || |
1620 ! jpoints = 0, | t t || t t | ... | t t || t t |
1621 ! nj=1 |-------------||-------------| |-------------||-------------|
1622 ! ^ ^
1623 ! | |
1624 ! | |
1625 ! ipos ipos+
1626 ! ni = 1 3 (nide-nids)/nri
1627 ! ipoints= 0 1 0 1 -istag
1628 !
1630 ! For performance benefits, users can comment out the inner most loop (and cfld=0) and
1631 ! hardcode the loop feedback. For example, it is set up to run a 2::1 ratio
1632 ! if uncommented. This lacks generality, but is likely to gain timing benefits
1633 ! with compilers unable to unroll inner loops that do not have parameterized sizes.
1635 ! The extra +1 ---------/ and the extra -1 ----\ (both for ci and cj)
1636 ! / \ keeps the feedback out of the
1637 ! / \ outer row/col, since that CG data
1638 ! / \ specified the nest boundary originally
1639 ! / \ This
1640 ! / \ is just
1641 ! / \ a sentence to not end a line
1642 ! / \ with a stupid backslash
1643 DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
1644 nj = (cj-jpos)*nrj + jstag
1645 DO ck = ckts, ckte
1646 nk = ck
1647 DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
1648 ni = (ci-ipos)*nri + istag
1649 cfld( ci, ck, cj ) = 0.
1650 DO ijpoints = 1 , nri * nrj
1651 ipoints = MOD((ijpoints-1),nri)
1652 jpoints = (ijpoints-1)/nri
1653 cfld( ci, ck, cj ) = cfld( ci, ck, cj ) + &
1654 1./REAL(nri*nrj) * nfld( ni+ipoints , nk , nj+jpoints )
1655 END DO
1656 ! cfld( ci, ck, cj ) = 1./4. * &
1657 ! ( nfld( ni , nk , nj ) + &
1658 ! nfld( ni+1, nk , nj ) + &
1659 ! nfld( ni , nk , nj+1) + &
1660 ! nfld( ni+1, nk , nj+1) )
1661 END DO
1662 END DO
1663 END DO
1665 ! U
1667 ELSE IF ( ( xstag ) .AND. ( .NOT. ystag ) ) THEN
1668 ! |---------------|
1669 ! | |
1670 ! jpoints = 1 u u |
1671 ! | |
1672 ! U |
1673 ! | |
1674 ! jpoints = 0, u u |
1675 ! nj=3 | |
1676 ! |---------------|
1677 ! |---------------|
1678 ! | |
1679 ! jpoints = 1 u u |
1680 ! | |
1681 ! jpos U |
1682 ! | |
1683 ! jpoints = 0, u u |
1684 ! nj=1 | |
1685 ! |---------------|
1686 !
1687 ! ^
1688 ! |
1689 ! |
1690 ! ipos
1691 ! ni = 1 3
1692 ! ipoints= 0 1 0
1693 !
1695 DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
1696 nj = (cj-jpos)*nrj + 1
1697 DO ck = ckts, ckte
1698 nk = ck
1699 DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
1700 ni = (ci-ipos)*nri + 1
1701 cfld( ci, ck, cj ) = 0.
1702 DO ijpoints = 1 , nri*nrj , nri
1703 ipoints = MOD((ijpoints-1),nri)
1704 jpoints = (ijpoints-1)/nri
1705 cfld( ci, ck, cj ) = cfld( ci, ck, cj ) + &
1706 1./REAL(nri ) * nfld( ni+ipoints , nk , nj+jpoints )
1707 END DO
1708 ! cfld( ci, ck, cj ) = 1./2. * &
1709 ! ( nfld( ni , nk , nj ) + &
1710 ! nfld( ni , nk , nj+1) )
1711 ENDDO
1712 ENDDO
1713 ENDDO
1715 ! V
1717 ELSE IF ( ( .NOT. xstag ) .AND. ( ystag ) ) THEN
1718 DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
1719 nj = (cj-jpos)*nrj + 1
1720 DO ck = ckts, ckte
1721 nk = ck
1722 DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
1723 ni = (ci-ipos)*nri + 1
1724 cfld( ci, ck, cj ) = 0.
1725 DO ijpoints = 1 , nri
1726 ipoints = MOD((ijpoints-1),nri)
1727 jpoints = (ijpoints-1)/nri
1728 cfld( ci, ck, cj ) = cfld( ci, ck, cj ) + &
1729 1./REAL(nri ) * nfld( ni+ipoints , nk , nj+jpoints )
1730 END DO
1731 ! cfld( ci, ck, cj ) = 1./2. * &
1732 ! ( nfld( ni , nk , nj ) + &
1733 ! nfld( ni+1, nk , nj ) )
1734 ENDDO
1735 ENDDO
1736 ENDDO
1737 END IF
1738 END IF
1740 RETURN
1742 END SUBROUTINE copy_fcn
1744 !==================================
1745 ! this is the 1pt function used in feedback.
1747 SUBROUTINE copy_fcnm ( cfld, & ! CD field
1748 cids, cide, ckds, ckde, cjds, cjde, &
1749 cims, cime, ckms, ckme, cjms, cjme, &
1750 cits, cite, ckts, ckte, cjts, cjte, &
1751 nfld, & ! ND field
1752 nids, nide, nkds, nkde, njds, njde, &
1753 nims, nime, nkms, nkme, njms, njme, &
1754 nits, nite, nkts, nkte, njts, njte, &
1755 shw, & ! stencil half width for interp
1756 imask, & ! interpolation mask
1757 xstag, ystag, & ! staggering of field
1758 ipos, jpos, & ! Position of lower left of nest in CD
1759 nri, nrj ) ! nest ratios
1760 USE module_configure
1761 USE module_wrf_error
1762 IMPLICIT NONE
1765 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
1766 cims, cime, ckms, ckme, cjms, cjme, &
1767 cits, cite, ckts, ckte, cjts, cjte, &
1768 nids, nide, nkds, nkde, njds, njde, &
1769 nims, nime, nkms, nkme, njms, njme, &
1770 nits, nite, nkts, nkte, njts, njte, &
1771 shw, &
1772 ipos, jpos, &
1773 nri, nrj
1774 LOGICAL, INTENT(IN) :: xstag, ystag
1776 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ), INTENT(OUT) :: cfld
1777 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ), INTENT(IN) :: nfld
1778 INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: imask
1780 ! Local
1782 INTEGER ci, cj, ck, ni, nj, nk, ip, jp, ioff, joff, ioffa, joffa
1783 INTEGER :: icmin,icmax,jcmin,jcmax
1784 INTEGER :: istag,jstag, ipoints,jpoints,ijpoints
1785 INTEGER , PARAMETER :: passes = 2
1786 INTEGER spec_zone
1788 CALL nl_get_spec_zone( 1, spec_zone )
1789 istag = 1 ; jstag = 1
1790 IF ( xstag ) istag = 0
1791 IF ( ystag ) jstag = 0
1793 IF( MOD(nrj,2) .NE. 0) THEN ! odd refinement ratio
1795 DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
1796 nj = (cj-jpos)*nrj + jstag + 1
1797 DO ck = ckts, ckte
1798 nk = ck
1799 DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
1800 ni = (ci-ipos)*nri + istag + 1
1801 cfld( ci, ck, cj ) = nfld( ni , nk , nj )
1802 ENDDO
1803 ENDDO
1804 ENDDO
1806 ELSE ! even refinement ratio, pick nearest neighbor on SW corner
1807 DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
1808 nj = (cj-jpos)*nrj + 1
1809 DO ck = ckts, ckte
1810 nk = ck
1811 DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
1812 ni = (ci-ipos)*nri + 1
1813 ipoints = nri/2 -1
1814 jpoints = nrj/2 -1
1815 cfld( ci, ck, cj ) = nfld( ni+ipoints , nk , nj+jpoints )
1816 END DO
1817 END DO
1818 END DO
1820 END IF
1822 RETURN
1824 END SUBROUTINE copy_fcnm
1826 !==================================
1827 ! this is the 1pt function used in feedback for integers
1829 SUBROUTINE copy_fcni ( cfld, & ! CD field
1830 cids, cide, ckds, ckde, cjds, cjde, &
1831 cims, cime, ckms, ckme, cjms, cjme, &
1832 cits, cite, ckts, ckte, cjts, cjte, &
1833 nfld, & ! ND field
1834 nids, nide, nkds, nkde, njds, njde, &
1835 nims, nime, nkms, nkme, njms, njme, &
1836 nits, nite, nkts, nkte, njts, njte, &
1837 shw, & ! stencil half width for interp
1838 imask, & ! interpolation mask
1839 xstag, ystag, & ! staggering of field
1840 ipos, jpos, & ! Position of lower left of nest in CD
1841 nri, nrj ) ! nest ratios
1842 USE module_configure
1843 USE module_wrf_error
1844 IMPLICIT NONE
1847 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
1848 cims, cime, ckms, ckme, cjms, cjme, &
1849 cits, cite, ckts, ckte, cjts, cjte, &
1850 nids, nide, nkds, nkde, njds, njde, &
1851 nims, nime, nkms, nkme, njms, njme, &
1852 nits, nite, nkts, nkte, njts, njte, &
1853 shw, &
1854 ipos, jpos, &
1855 nri, nrj
1856 LOGICAL, INTENT(IN) :: xstag, ystag
1858 INTEGER, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ), INTENT(OUT) :: cfld
1859 INTEGER, DIMENSION ( nims:nime, nkms:nkme, njms:njme ), INTENT(IN) :: nfld
1860 INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: imask
1862 ! Local
1864 INTEGER ci, cj, ck, ni, nj, nk, ip, jp, ioff, joff, ioffa, joffa
1865 INTEGER :: icmin,icmax,jcmin,jcmax
1866 INTEGER :: istag,jstag, ipoints,jpoints,ijpoints
1867 INTEGER , PARAMETER :: passes = 2
1868 INTEGER spec_zone
1870 CALL nl_get_spec_zone( 1, spec_zone )
1871 istag = 1 ; jstag = 1
1872 IF ( xstag ) istag = 0
1873 IF ( ystag ) jstag = 0
1875 IF( MOD(nrj,2) .NE. 0) THEN ! odd refinement ratio
1877 DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
1878 nj = (cj-jpos)*nrj + jstag + 1
1879 DO ck = ckts, ckte
1880 nk = ck
1881 DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
1882 ni = (ci-ipos)*nri + istag + 1
1883 cfld( ci, ck, cj ) = nfld( ni , nk , nj )
1884 ENDDO
1885 ENDDO
1886 ENDDO
1888 ELSE ! even refinement ratio
1889 DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
1890 nj = (cj-jpos)*nrj + 1
1891 DO ck = ckts, ckte
1892 nk = ck
1893 DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
1894 ni = (ci-ipos)*nri + 1
1895 ipoints = nri/2 -1
1896 jpoints = nrj/2 -1
1897 cfld( ci, ck, cj ) = nfld( ni+ipoints , nk , nj+jpoints )
1898 END DO
1899 END DO
1900 END DO
1902 END IF
1904 RETURN
1906 END SUBROUTINE copy_fcni
1908 !==================================
1910 SUBROUTINE vert_interp_vert_nesting ( cfld, & ! CD field
1911 ids, ide, kds, kde, jds, jde, &
1912 ims, ime, kms, kme, jms, jme, &
1913 its, ite, kts, kte, jts, jte, &
1914 pgrid_s_vert, pgrid_e_vert, & ! vertical dimensions of parent grid
1915 cf1_c, cf2_c, cf3_c, cfn_c, cfn1_c, & ! coarse grid extrapolation constants
1916 alt_u_c, alt_u_n)
1918 !KAL vertical interpolation for u, v, and mass points (w is below in a different subroutine) for vertical nesting
1920 IMPLICIT NONE
1921 REAL, DIMENSION ( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) :: cfld
1922 INTEGER, INTENT(IN) :: ids, ide, kds, kde, jds, jde, &
1923 ims, ime, kms, kme, jms, jme, &
1924 its, ite, kts, kte, jts, jte
1925 INTEGER, INTENT(IN) :: pgrid_s_vert, pgrid_e_vert ! vertical dimensions of the parent grid
1926 REAL, INTENT(IN) :: cf1_c, cf2_c, cf3_c, cfn_c, cfn1_c
1927 REAL, DIMENSION(pgrid_s_vert:pgrid_e_vert+1), INTENT(IN) :: alt_u_c
1928 REAL, DIMENSION(kde+1), INTENT(IN) :: alt_u_n
1931 !local
1933 INTEGER :: i,j,k
1934 REAL, DIMENSION(pgrid_s_vert:pgrid_e_vert+1) :: pro_u_c ! variable in 1D on the coarse grid
1935 REAL, DIMENSION(kde+1) :: pro_u_n
1937 DO j = jms,jme
1938 DO i = ims,ime
1940 ! pro_u_c is u on the 1D coarse grid
1942 do k = pgrid_s_vert,pgrid_e_vert-1
1943 pro_u_c(k+1) = cfld(i,k,j)
1944 enddo
1946 !KAL fill in the surface value and the top value using extrapolation
1948 pro_u_c(1 ) = cf1_c*cfld(i,1,j) &
1949 + cf2_c*cfld(i,2,j) &
1950 + cf3_c*cfld(i,3,j)
1952 pro_u_c(pgrid_e_vert+1) = cfn_c *cfld(i,pgrid_e_vert-1,j) &
1953 + cfn1_c*cfld(i,pgrid_e_vert-2,j)
1955 call inter_wrf_copy(pro_u_c, alt_u_c, pgrid_e_vert+1, pro_u_n, alt_u_n, kde+1)
1957 do k = 1,kde-1
1958 cfld(i,k,j) = pro_u_n(k+1)
1959 enddo
1961 ENDDO
1962 ENDDO
1965 END SUBROUTINE vert_interp_vert_nesting
1967 !==================================
1969 SUBROUTINE vert_interp_vert_nesting_w ( cfld, & ! CD field
1970 ids, ide, kds, kde, jds, jde, &
1971 ims, ime, kms, kme, jms, jme, &
1972 its, ite, kts, kte, jts, jte, &
1973 pgrid_s_vert, pgrid_e_vert, & ! vertical dimensions of parent grid
1974 alt_w_c, alt_w_n)
1976 !KAL vertical interpolation at w points for vertical nesting
1978 IMPLICIT NONE
1979 REAL, DIMENSION ( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) :: cfld
1980 INTEGER, INTENT(IN) :: ids, ide, kds, kde, jds, jde, &
1981 ims, ime, kms, kme, jms, jme, &
1982 its, ite, kts, kte, jts, jte
1983 INTEGER, INTENT(IN) :: pgrid_s_vert, pgrid_e_vert ! vertical dimensions of the parent grid
1984 REAL, DIMENSION(pgrid_s_vert:pgrid_e_vert), INTENT(IN) :: alt_w_c
1985 REAL, DIMENSION(kde), INTENT(IN) :: alt_w_n
1988 !local
1990 INTEGER :: i,j,k
1991 REAL, DIMENSION(pgrid_s_vert:pgrid_e_vert) :: pro_w_c ! variable in 1D on the coarse grid
1992 REAL, DIMENSION(kde) :: pro_w_n
1994 DO j = jms,jme
1995 DO i = ims,ime
1997 ! pro_w_c is w on the 1D coarse grid
1999 do k = pgrid_s_vert,pgrid_e_vert
2000 pro_w_c(k) = cfld(i,k,j)
2001 enddo
2003 call inter_wrf_copy(pro_w_c, alt_w_c, pgrid_e_vert, pro_w_n, alt_w_n, kde)
2005 do k = 1,kde
2006 cfld(i,k,j) = pro_w_n(k)
2007 enddo
2009 ENDDO
2010 ENDDO
2013 END SUBROUTINE vert_interp_vert_nesting_w
2015 !-----------------------------------------------------------------------------------------
2017 SUBROUTINE vert_interp_vert_nesting_1d ( cfld, & ! CD field
2018 ids, ide, kds, kde, jds, jde, &
2019 ims, ime, kms, kme, jms, jme, &
2020 its, ite, kts, kte, jts, jte, &
2021 pgrid_s_vert, pgrid_e_vert, & ! vertical dimensions of parent grid
2022 cf1_c, cf2_c, cf3_c, cfn_c, cfn1_c, & ! coarse grid extrapolation constants
2023 alt_u_c, alt_u_n)
2025 !KAL vertical interpolation for u, v, and mass points (w is below in a different subroutine) for vertical nesting
2027 IMPLICIT NONE
2028 REAL, DIMENSION (kms:kme),INTENT(INOUT) :: cfld
2029 INTEGER, INTENT(IN) :: ids, ide, kds, kde, jds, jde, &
2030 ims, ime, kms, kme, jms, jme, &
2031 its, ite, kts, kte, jts, jte
2032 INTEGER, INTENT(IN) :: pgrid_s_vert, pgrid_e_vert ! vertical dimensions of the parent grid
2033 REAL, INTENT(IN) :: cf1_c, cf2_c, cf3_c, cfn_c, cfn1_c
2034 REAL, DIMENSION(pgrid_s_vert:pgrid_e_vert+1), INTENT(IN) :: alt_u_c
2035 REAL, DIMENSION(kde+1), INTENT(IN) :: alt_u_n
2038 !local
2040 INTEGER :: i,j,k
2041 REAL, DIMENSION(pgrid_s_vert:pgrid_e_vert+1) :: pro_u_c ! variable in 1D on the coarse grid
2042 REAL, DIMENSION(kde+1) :: pro_u_n
2045 ! pro_u_c is u on the 1D coarse grid
2047 do k = pgrid_s_vert,pgrid_e_vert-1
2048 pro_u_c(k+1) = cfld(k)
2049 enddo
2051 !KAL fill in the surface value and the top value using extrapolation
2053 pro_u_c(1 ) = cf1_c*cfld(1) &
2054 + cf2_c*cfld(2) &
2055 + cf3_c*cfld(3)
2057 pro_u_c(pgrid_e_vert+1) = cfn_c *cfld(pgrid_e_vert-1) &
2058 + cfn1_c*cfld(pgrid_e_vert-2)
2060 call inter_wrf_copy(pro_u_c, alt_u_c, pgrid_e_vert+1, pro_u_n, alt_u_n, kde+1)
2062 do k = 1,kde-1
2063 cfld(k) = pro_u_n(k+1)
2064 enddo
2067 END SUBROUTINE vert_interp_vert_nesting_1d
2070 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2071 !KAL this is a direct copy of a subroutine from ndown, but a dependency on ndown will not work because it is not always compiled (for ideal cases), and most likely not compliled in the order needed.
2073 SUBROUTINE inter_wrf_copy(pro_c,alt_c,kde_c,pro_n,alt_n,kde_n)
2075 !KAL this routine has been added for vertical nesting
2077 IMPLICIT NONE
2078 INTEGER , INTENT(IN) :: kde_c,kde_n
2079 REAL , DIMENSION(kde_c) , INTENT(IN ) :: pro_c,alt_c
2080 REAL , DIMENSION(kde_n) , INTENT(IN ) :: alt_n
2081 REAL , DIMENSION(kde_n) , INTENT(OUT) :: pro_n
2083 real ,dimension(kde_c) :: a,b,c,d
2084 real :: p
2085 integer :: i,j
2088 call coeff_mon_wrf_copy(alt_c,pro_c,a,b,c,d,kde_c)
2090 do i = 1,kde_n-1
2092 do j=1,kde_c-1
2094 if ( (alt_n(i) .ge. alt_c(j)).and.(alt_n(i) .lt. alt_c(j+1)) ) then
2095 p = alt_n(i)-alt_c(j)
2096 pro_n(i) = p*( p*(a(j)*p+b(j))+c(j)) + d(j)
2097 goto 20
2098 endif
2099 enddo
2100 20 continue
2101 enddo
2103 pro_n(kde_n) = pro_c(kde_c)
2106 END SUBROUTINE inter_wrf_copy
2108 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
2109 !KAL this is a direct copy of a subroutine from ndown, but a dependency on ndown will not work because it is not always compiled (for ideal cases), and most likely not compliled in the order needed.
2111 subroutine coeff_mon_wrf_copy(x,y,a,b,c,d,n)
2113 !KAL this routine has been added for vertical nesting
2115 implicit none
2117 integer :: n
2118 real ,dimension(n) :: x,y,a,b,c,d
2119 real ,dimension(n) :: h,s,p,yp
2121 integer :: i
2124 do i=1,n-1
2125 h(i) = (x(i+1)-x(i))
2126 s(i) = (y(i+1)-y(i)) / h(i)
2127 enddo
2129 do i=2,n-1
2130 p(i) = (s(i-1)*h(i)+s(i)*h(i-1)) / (h(i-1)+h(i))
2131 enddo
2133 p(1) = s(1)
2134 p(n) = s(n-1)
2136 do i=1,n
2137 yp(i) = p(i)
2138 enddo
2139 !!!!!!!!!!!!!!!!!!!!!
2141 do i=2,n-1
2142 yp(i) = (sign(1.,s(i-1))+sign(1.,s(i)))* min( abs(s(i-1)),abs(s(i)),0.5*abs(p(i)))
2143 enddo
2145 do i = 1,n-1
2146 a(i) = (yp(i)+yp(i+1)-2.*s(i))/(h(i)*h(i))
2147 b(i) = (3.*s(i)-2.*yp(i)-yp(i+1))/h(i)
2148 c(i) = yp(i)
2149 d(i) = y(i)
2150 enddo
2152 end subroutine coeff_mon_wrf_copy
2154 !-----------------------------------------------------------------------------------------
2157 !==================================
2159 SUBROUTINE p2c ( cfld, & ! CD field
2160 cids, cide, ckds, ckde, cjds, cjde, &
2161 cims, cime, ckms, ckme, cjms, cjme, &
2162 cits, cite, ckts, ckte, cjts, cjte, &
2163 nfld, & ! ND field
2164 nids, nide, nkds, nkde, njds, njde, &
2165 nims, nime, nkms, nkme, njms, njme, &
2166 nits, nite, nkts, nkte, njts, njte, &
2167 shw, & ! stencil half width
2168 imask, & ! interpolation mask
2169 xstag, ystag, & ! staggering of field
2170 ipos, jpos, & ! Position of lower left of nest in CD
2171 nri, nrj & ! nest ratios
2172 )
2173 USE module_configure
2174 IMPLICIT NONE
2176 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
2177 cims, cime, ckms, ckme, cjms, cjme, &
2178 cits, cite, ckts, ckte, cjts, cjte, &
2179 nids, nide, nkds, nkde, njds, njde, &
2180 nims, nime, nkms, nkme, njms, njme, &
2181 nits, nite, nkts, nkte, njts, njte, &
2182 shw, &
2183 ipos, jpos, &
2184 nri, nrj
2186 LOGICAL, INTENT(IN) :: xstag, ystag
2188 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
2189 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
2190 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
2192 CALL interp_fcn (cfld, & ! CD field
2193 cids, cide, ckds, ckde, cjds, cjde, &
2194 cims, cime, ckms, ckme, cjms, cjme, &
2195 cits, cite, ckts, ckte, cjts, cjte, &
2196 nfld, & ! ND field
2197 nids, nide, nkds, nkde, njds, njde, &
2198 nims, nime, nkms, nkme, njms, njme, &
2199 nits, nite, nkts, nkte, njts, njte, &
2200 shw, & ! stencil half width for interp
2201 imask, & ! interpolation mask
2202 xstag, ystag, & ! staggering of field
2203 ipos, jpos, & ! Position of lower left of nest in CD
2204 nri, nrj ) ! nest ratios
2206 END SUBROUTINE p2c
2208 !==================================
2210 SUBROUTINE c2f_interp ( cfld, & ! CD field
2211 cids, cide, ckds, ckde, cjds, cjde, &
2212 cims, cime, ckms, ckme, cjms, cjme, &
2213 cits, cite, ckts, ckte, cjts, cjte, &
2214 nfld, & ! ND field
2215 nids, nide, nkds, nkde, njds, njde, &
2216 nims, nime, nkms, nkme, njms, njme, &
2217 nits, nite, nkts, nkte, njts, njte, &
2218 shw, & ! stencil half width
2219 imask, & ! interpolation mask
2220 xstag, ystag, & ! staggering of field
2221 ipos, jpos, & ! Position of lower left of nest in CD
2222 nri, nrj, & ! nest ratios
2223 ! cbdy_xs, nbdy_xs, &
2224 ! cbdy_xe, nbdy_xe, &
2225 ! cbdy_ys, nbdy_ys, &
2226 ! cbdy_ye, nbdy_ye, &
2227 ! cbdy_txs, nbdy_txs, &
2228 ! cbdy_txe, nbdy_txe, &
2229 ! cbdy_tys, nbdy_tys, &
2230 ! cbdy_tye, nbdy_tye, &
2231 parent_id,nest_id &!cyl
2232 ) ! boundary arrays
2233 USE module_configure
2234 IMPLICIT NONE
2236 !------------------------------------------------------------
2237 ! Subroutine c2f_interp interpolate field from coarse resolution domain
2238 ! to its nested domain. It is written by Dave Gill in NCAR for the purpose
2239 ! running phys/module_sf_oml.F-DPWP in only d01 and d02
2240 ! Chiaying Lee RSMAS/UM
2241 !------------------------------------------------------------
2243 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
2244 cims, cime, ckms, ckme, cjms, cjme, &
2245 cits, cite, ckts, ckte, cjts, cjte, &
2246 nids, nide, nkds, nkde, njds, njde, &
2247 nims, nime, nkms, nkme, njms, njme, &
2248 nits, nite, nkts, nkte, njts, njte, &
2249 shw, &
2250 ipos, jpos, &
2251 nri, nrj,parent_id,nest_id !cyl
2253 LOGICAL, INTENT(IN) :: xstag, ystag
2255 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
2256 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
2257 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
2258 ! REAL, DIMENSION( * ), INTENT(INOUT) :: cbdy_xs, cbdy_txs, nbdy_xs, nbdy_txs
2259 ! REAL, DIMENSION( * ), INTENT(INOUT) :: cbdy_xe, cbdy_txe, nbdy_xe, nbdy_txe
2260 ! REAL, DIMENSION( * ), INTENT(INOUT) :: cbdy_ys, cbdy_tys, nbdy_ys, nbdy_tys
2261 ! REAL, DIMENSION( * ), INTENT(INOUT) :: cbdy_ye, cbdy_tye, nbdy_ye, nbdy_tye
2262 REAL cdt, ndt
2264 ! Local
2266 INTEGER ci, cj, ck, ni, nj, nk, ip, jp
2268 ! Iterate over the ND tile and compute the values
2269 ! from the CD tile.
2271 !write(0,'("mask cits:cite, ckts:ckte, cjts:cjte ",6i4)')cits,cite, ckts,ckte, cjts,cjte
2272 !write(0,'("mask nits:nite, nkts:nkte, njts:njte ",6i4)')nits,nite, nkts,nkte, njts,njte
2273 ! write(0,*)'cyl parentid',parent_id
2274 ! write(0,*)'cyl nestid',nest_id
2275 ! If ( nest_id .le. 2 .and. (1.0/rdx .ge. 3000.0 .and. 1.0/rdy .ge. 3000.0) ) then ! cyl: only run it in the nest domain with dx, dy < 3 km
2276 If ( nest_id .eq. 3 ) then
2277 DO nj = njts, njte
2279 cj = jpos + (nj-1) / nrj ! j coord of CD point
2280 jp = mod ( nj , nrj ) ! coord of ND w/i CD point
2281 DO nk = nkts, nkte
2282 ck = nk
2283 DO ni = nits, nite
2284 ci = ipos + (ni-1) / nri ! j coord of CD point
2285 ip = mod ( ni , nri ) ! coord of ND w/i CD point
2286 ! This is a trivial implementation of the interp_fcn; just copies
2287 ! the values from the CD into the ND
2288 nfld( ni, nk, nj ) = cfld( ci , ck , cj )
2289 ENDDO
2290 ENDDO
2291 ENDDO
2292 ENDIF ! cyl
2293 RETURN
2295 END SUBROUTINE c2f_interp
2297 !==================================
2299 SUBROUTINE bdy_interp ( cfld, & ! CD field
2300 cids, cide, ckds, ckde, cjds, cjde, &
2301 cims, cime, ckms, ckme, cjms, cjme, &
2302 cits, cite, ckts, ckte, cjts, cjte, &
2303 nfld, & ! ND field
2304 nids, nide, nkds, nkde, njds, njde, &
2305 nims, nime, nkms, nkme, njms, njme, &
2306 nits, nite, nkts, nkte, njts, njte, &
2307 shw, & ! stencil half width
2308 imask, & ! interpolation mask
2309 xstag, ystag, & ! staggering of field
2310 ipos, jpos, & ! Position of lower left of nest in CD
2311 nri, nrj, & ! nest ratios
2312 cbdy_xs, nbdy_xs, & ! Boundary components, for the CG and FG,
2313 cbdy_xe, nbdy_xe, & ! for each of the four sides (x start, e end, y start, y end)
2314 cbdy_ys, nbdy_ys, & ! and also for the two components of the LBC:
2315 cbdy_ye, nbdy_ye, & ! the "starting" value and the tendency
2316 cbdy_txs, nbdy_txs, &
2317 cbdy_txe, nbdy_txe, & ! The CG is at a parent time step ahead of the FG
2318 cbdy_tys, nbdy_tys, &
2319 cbdy_tye, nbdy_tye, &
2320 cdt, ndt ) ! Time step size for CG and FG
2322 USE module_interp_info
2324 IMPLICIT NONE
2326 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
2327 cims, cime, ckms, ckme, cjms, cjme, &
2328 cits, cite, ckts, ckte, cjts, cjte, &
2329 nids, nide, nkds, nkde, njds, njde, &
2330 nims, nime, nkms, nkme, njms, njme, &
2331 nits, nite, nkts, nkte, njts, njte, &
2332 shw, &
2333 ipos, jpos, &
2334 nri, nrj
2336 LOGICAL, INTENT(IN) :: xstag, ystag
2338 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
2339 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
2340 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
2341 REAL, DIMENSION( * ), INTENT(INOUT) :: cbdy_xs, cbdy_txs, nbdy_xs, nbdy_txs
2342 REAL, DIMENSION( * ), INTENT(INOUT) :: cbdy_xe, cbdy_txe, nbdy_xe, nbdy_txe
2343 REAL, DIMENSION( * ), INTENT(INOUT) :: cbdy_ys, cbdy_tys, nbdy_ys, nbdy_tys
2344 REAL, DIMENSION( * ), INTENT(INOUT) :: cbdy_ye, cbdy_tye, nbdy_ye, nbdy_tye
2345 REAL cdt, ndt
2347 ! Local
2349 INTEGER nijds, nijde, spec_bdy_width
2351 nijds = min(nids, njds)
2352 nijde = max(nide, njde)
2353 CALL nl_get_spec_bdy_width( 1, spec_bdy_width )
2355 IF ( interp_method_type .EQ. NOT_DEFINED_YET ) THEN
2356 interp_method_type = SINT
2357 END IF
2359 IF ( interp_method_type .EQ. SINT ) THEN
2360 CALL bdy_interp1( cfld, & ! CD field
2361 cids, cide, ckds, ckde, cjds, cjde, &
2362 cims, cime, ckms, ckme, cjms, cjme, &
2363 cits, cite, ckts, ckte, cjts, cjte, &
2364 nfld, & ! ND field
2365 nijds, nijde , & ! start and end of nest LBC size in the LONG direction
2366 spec_bdy_width , & ! width of the LBC, the SHORT direction
2367 nids, nide, nkds, nkde, njds, njde, &
2368 nims, nime, nkms, nkme, njms, njme, &
2369 nits, nite, nkts, nkte, njts, njte, &
2370 shw, imask, &
2371 xstag, ystag, & ! staggering of field
2372 ipos, jpos, & ! Position of lower left of nest in CD
2373 nri, nrj, &
2374 cbdy_xs, nbdy_xs, & ! Boundary components, for the CG and FG,
2375 cbdy_xe, nbdy_xe, & ! for each of the four sides (x start, e end, y start, y end)
2376 cbdy_ys, nbdy_ys, & ! and also for the two components of the LBC:
2377 cbdy_ye, nbdy_ye, & ! the "starting" value and the tendency
2378 cbdy_txs, nbdy_txs, &
2379 cbdy_txe, nbdy_txe, & ! The CG is at a parent time step ahead of the FG
2380 cbdy_tys, nbdy_tys, &
2381 cbdy_tye, nbdy_tye, &
2382 cdt, ndt & ! Time step size for CG and FG
2383 )
2385 ELSE IF ( ( interp_method_type .EQ. BILINEAR ) .OR. &
2386 ( interp_method_type .EQ. NEAREST_NEIGHBOR ) .OR. &
2387 ( interp_method_type .EQ. QUADRATIC ) .OR. &
2388 ( interp_method_type .EQ. SPLINE ) .OR. &
2389 ( interp_method_type .EQ. SINT_NEW ) ) THEN
2390 CALL bdy_interp2( cfld, & ! CD field
2391 cids, cide, ckds, ckde, cjds, cjde, &
2392 cims, cime, ckms, ckme, cjms, cjme, &
2393 cits, cite, ckts, ckte, cjts, cjte, &
2394 nfld, & ! ND field
2395 nijds, nijde , & ! start and end of nest LBC size in the LONG direction
2396 spec_bdy_width , & ! width of the LBC, the SHORT direction
2397 nids, nide, nkds, nkde, njds, njde, &
2398 nims, nime, nkms, nkme, njms, njme, &
2399 nits, nite, nkts, nkte, njts, njte, &
2400 shw, imask, &
2401 xstag, ystag, & ! staggering of field
2402 ipos, jpos, & ! Position of lower left of nest in CD
2403 nri, nrj, &
2404 cbdy_xs, nbdy_xs, & ! Boundary components, for the CG and FG,
2405 cbdy_xe, nbdy_xe, & ! for each of the four sides (x start, e end, y start, y end)
2406 cbdy_ys, nbdy_ys, & ! and also for the two components of the LBC:
2407 cbdy_ye, nbdy_ye, & ! the "starting" value and the tendency
2408 cbdy_txs, nbdy_txs, &
2409 cbdy_txe, nbdy_txe, & ! The CG is at a parent time step ahead of the FG
2410 cbdy_tys, nbdy_tys, &
2411 cbdy_tye, nbdy_tye, &
2412 cdt, ndt & ! Time step size for CG and FG
2413 )
2415 ELSE
2416 CALL wrf_error_fatal ('Hold on there cowboy #2, we need to know which nested lateral boundary interpolation option you want')
2417 END IF
2419 END SUBROUTINE bdy_interp
2421 !==================================
2423 SUBROUTINE bdy_interp1( cfld, & ! CD field
2424 cids, cide, ckds, ckde, cjds, cjde, &
2425 cims, cime, ckms, ckme, cjms, cjme, &
2426 cits, cite, ckts, ckte, cjts, cjte, &
2427 nfld, & ! ND field
2428 nijds, nijde, spec_bdy_width , &
2429 nids, nide, nkds, nkde, njds, njde, &
2430 nims, nime, nkms, nkme, njms, njme, &
2431 nits, nite, nkts, nkte, njts, njte, &
2432 shw1, &
2433 imask, & ! interpolation mask
2434 xstag, ystag, & ! staggering of field
2435 ipos, jpos, & ! Position of lower left of nest in CD
2436 nri, nrj, &
2437 cbdy_xs, bdy_xs, &
2438 cbdy_xe, bdy_xe, &
2439 cbdy_ys, bdy_ys, &
2440 cbdy_ye, bdy_ye, &
2441 cbdy_txs, bdy_txs, &
2442 cbdy_txe, bdy_txe, &
2443 cbdy_tys, bdy_tys, &
2444 cbdy_tye, bdy_tye, &
2445 cdt, ndt &
2446 )
2448 ! USE module_configure , ONLY : nl_get_spec_zone, nl_get_relax_zone
2449 USE module_state_description
2451 IMPLICIT NONE
2453 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
2454 cims, cime, ckms, ckme, cjms, cjme, &
2455 cits, cite, ckts, ckte, cjts, cjte, &
2456 nids, nide, nkds, nkde, njds, njde, &
2457 nims, nime, nkms, nkme, njms, njme, &
2458 nits, nite, nkts, nkte, njts, njte, &
2459 shw1, & ! ignore
2460 ipos, jpos, &
2461 nri, nrj
2462 INTEGER, INTENT(IN) :: nijds, nijde, spec_bdy_width
2463 LOGICAL, INTENT(IN) :: xstag, ystag
2465 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ), INTENT(INOUT) :: cfld
2466 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ), INTENT(INOUT) :: nfld
2467 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
2468 REAL, DIMENSION ( * ), INTENT(INOUT) :: cbdy_xs, cbdy_txs ! not used
2469 REAL, DIMENSION ( * ), INTENT(INOUT) :: cbdy_xe, cbdy_txe ! not used
2470 REAL, DIMENSION ( * ), INTENT(INOUT) :: cbdy_ys, cbdy_tys ! not used
2471 REAL, DIMENSION ( * ), INTENT(INOUT) :: cbdy_ye, cbdy_tye ! not used
2472 REAL :: cdt, ndt
2473 REAL, DIMENSION ( njms:njme, nkms:nkme, spec_bdy_width ), INTENT(INOUT) :: bdy_xs, bdy_txs
2474 REAL, DIMENSION ( njms:njme, nkms:nkme, spec_bdy_width ), INTENT(INOUT) :: bdy_xe, bdy_txe
2475 REAL, DIMENSION ( nims:nime, nkms:nkme, spec_bdy_width ), INTENT(INOUT) :: bdy_ys, bdy_tys
2476 REAL, DIMENSION ( nims:nime, nkms:nkme, spec_bdy_width ), INTENT(INOUT) :: bdy_ye, bdy_tye
2478 ! Local
2480 REAL*8 rdt
2481 INTEGER ci, cj, ck, ni, nj, nk, ni1, nj1, nk1, ip, jp, ioff, joff
2482 INTEGER nfx, ior
2483 PARAMETER (ior=2)
2484 INTEGER nf
2485 REAL psca1(cims:cime,cjms:cjme,nri*nrj)
2486 REAL psca(cims:cime,cjms:cjme,nri*nrj)
2487 LOGICAL icmask( cims:cime, cjms:cjme )
2488 INTEGER i,j,k
2489 INTEGER shw
2490 INTEGER spec_zone
2491 INTEGER relax_zone
2492 INTEGER sz
2493 INTEGER n2ci,n
2494 INTEGER n2cj
2496 ! statement functions for converting a nest index to coarse
2497 n2ci(n) = (n+ipos*nri-1)/nri
2498 n2cj(n) = (n+jpos*nrj-1)/nrj
2500 rdt = 1.D0/cdt
2502 shw = 0
2504 ioff = 0 ; joff = 0
2505 IF ( xstag ) THEN
2506 ioff = MAX((nri-1)/2,1)
2507 ENDIF
2508 IF ( ystag ) THEN
2509 joff = MAX((nrj-1)/2,1)
2510 ENDIF
2512 ! Iterate over the ND tile and compute the values
2513 ! from the CD tile.
2515 CALL nl_get_spec_zone( 1, spec_zone )
2516 CALL nl_get_relax_zone( 1, relax_zone )
2517 sz = MIN(MAX( spec_zone, relax_zone + 1 ),spec_bdy_width)
2519 nfx = nri * nrj
2521 !$OMP PARALLEL DO &
2522 !$OMP PRIVATE ( i,j,k,ni,nj,ni1,nj1,ci,cj,ip,jp,nk,ck,nf,icmask,psca,psca1 )
2523 DO k = ckts, ckte
2525 DO nf = 1,nfx
2526 DO j = cjms,cjme
2527 nj = (j-jpos) * nrj + ( nrj / 2 + 1 ) ! j point on nest
2528 DO i = cims,cime
2529 ni = (i-ipos) * nri + ( nri / 2 + 1 ) ! i point on nest
2530 psca1(i,j,nf) = cfld(i,k,j)
2531 ENDDO
2532 ENDDO
2533 ENDDO
2534 ! hopefully less ham handed but still correct and more efficient
2535 ! sintb ignores icmask so it does not matter that icmask is not set
2536 !
2537 ! SOUTH BDY
2538 IF ( njts .ge. njds .and. njts .le. njds + sz + joff ) THEN
2539 CALL sintb( psca1, psca, &
2540 cims, cime, cjms, cjme, icmask, &
2541 n2ci(nits)-1, n2ci(nite)+1, n2cj(MAX(njts,njds)), n2cj(MIN(njte,njds+sz+joff)), nrj*nri, xstag, ystag )
2542 ENDIF
2543 ! NORTH BDY
2544 IF ( njte .le. njde .and. njte .ge. njde - sz - joff ) THEN
2545 CALL sintb( psca1, psca, &
2546 cims, cime, cjms, cjme, icmask, &
2547 n2ci(nits)-1, n2ci(nite)+1, n2cj(MAX(njts,njde-sz-joff)), n2cj(MIN(njte,njde-1+joff)), nrj*nri, xstag, ystag )
2548 ENDIF
2549 ! WEST BDY
2550 IF ( nits .ge. nids .and. nits .le. nids + sz + ioff ) THEN
2551 CALL sintb( psca1, psca, &
2552 cims, cime, cjms, cjme, icmask, &
2553 n2ci(MAX(nits,nids)), n2ci(MIN(nite,nids+sz+ioff)), n2cj(njts)-1, n2cj(njte)+1, nrj*nri, xstag, ystag )
2554 ENDIF
2555 ! EAST BDY
2556 IF ( nite .le. nide .and. nite .ge. nide - sz - ioff ) THEN
2557 CALL sintb( psca1, psca, &
2558 cims, cime, cjms, cjme, icmask, &
2559 n2ci(MAX(nits,nide-sz-ioff)), n2ci(MIN(nite,nide-1+ioff)), n2cj(njts)-1, n2cj(njte)+1, nrj*nri, xstag, ystag )
2560 ENDIF
2562 DO nj1 = MAX(njds,njts-1), MIN(njde+joff,njte+joff+1)
2563 cj = jpos + (nj1-1) / nrj ! j coord of CD point
2564 jp = mod ( nj1-1 , nrj ) ! coord of ND w/i CD point
2565 nk = k
2566 ck = nk
2567 DO ni1 = MAX(nids,nits-1), MIN(nide+ioff,nite+ioff+1)
2568 ci = ipos + (ni1-1) / nri ! j coord of CD point
2569 ip = mod ( ni1-1 , nri ) ! coord of ND w/i CD point
2571 ni = ni1-ioff
2572 nj = nj1-joff
2574 IF ( ( ni.LT.nids) .OR. (nj.LT.njds) ) THEN
2575 CYCLE
2576 END IF
2578 !bdy contains the value at t-dt. psca contains the value at t
2579 !compute dv/dt and store in bdy_t
2580 !afterwards store the new value of v at t into bdy
2581 ! WEST
2582 IF ( ni .ge. nids .and. ni .lt. nids + sz ) THEN
2583 bdy_txs( nj,k,ni ) = rdt*(psca(ci+shw,cj+shw,ip+1+(jp)*nri)-nfld(ni,k,nj))
2584 bdy_xs( nj,k,ni ) = nfld(ni,k,nj)
2585 ENDIF
2587 ! SOUTH
2588 IF ( nj .ge. njds .and. nj .lt. njds + sz ) THEN
2589 bdy_tys( ni,k,nj ) = rdt*(psca(ci+shw,cj+shw,ip+1+(jp)*nri)-nfld(ni,k,nj))
2590 bdy_ys( ni,k,nj ) = nfld(ni,k,nj)
2591 ENDIF
2593 ! EAST
2594 IF ( xstag ) THEN
2595 IF ( ni .ge. nide - sz + 1 .AND. ni .le. nide ) THEN
2596 bdy_txe( nj,k,nide-ni+1 ) = rdt*(psca(ci+shw,cj+shw,ip+1+(jp)*nri)-nfld(ni,k,nj))
2597 bdy_xe( nj,k,nide-ni+1 ) = nfld(ni,k,nj)
2598 ENDIF
2599 ELSE
2600 IF ( ni .ge. nide - sz .AND. ni .le. nide-1 ) THEN
2601 bdy_txe( nj,k,nide-ni ) = rdt*(psca(ci+shw,cj+shw,ip+1+(jp)*nri)-nfld(ni,k,nj))
2602 bdy_xe( nj,k,nide-ni ) = nfld(ni,k,nj)
2603 ENDIF
2604 ENDIF
2606 ! NORTH
2607 IF ( ystag ) THEN
2608 IF ( nj .ge. njde - sz + 1 .AND. nj .le. njde ) THEN
2609 bdy_tye( ni,k,njde-nj+1 ) = rdt*(psca(ci+shw,cj+shw,ip+1+(jp)*nri)-nfld(ni,k,nj))
2610 bdy_ye( ni,k,njde-nj+1 ) = nfld(ni,k,nj)
2611 ENDIF
2612 ELSE
2613 IF ( nj .ge. njde - sz .AND. nj .le. njde-1 ) THEN
2614 bdy_tye(ni,k,njde-nj ) = rdt*(psca(ci+shw,cj+shw,ip+1+(jp)*nri)-nfld(ni,k,nj))
2615 bdy_ye( ni,k,njde-nj ) = nfld(ni,k,nj)
2616 ENDIF
2617 ENDIF
2619 ENDDO
2620 ENDDO
2621 ENDDO
2622 !$OMP END PARALLEL DO
2624 RETURN
2626 END SUBROUTINE bdy_interp1
2628 !==================================
2630 SUBROUTINE bdy_interp2( cfld, & ! CD field
2631 cids, cide, ckds, ckde, cjds, cjde, &
2632 cims, cime, ckms, ckme, cjms, cjme, &
2633 cits, cite, ckts, ckte, cjts, cjte, &
2634 nfld, & ! ND field
2635 nijds, nijde, spec_bdy_width , &
2636 nids, nide, nkds, nkde, njds, njde, &
2637 nims, nime, nkms, nkme, njms, njme, &
2638 nits, nite, nkts, nkte, njts, njte, &
2639 shw1, &
2640 imask, & ! interpolation mask
2641 xstag, ystag, & ! staggering of field
2642 ipos, jpos, & ! Position of lower left of nest in CD
2643 nri, nrj, &
2644 cbdy_xs, bdy_xs, &
2645 cbdy_xe, bdy_xe, &
2646 cbdy_ys, bdy_ys, &
2647 cbdy_ye, bdy_ye, &
2648 cbdy_txs, bdy_txs, &
2649 cbdy_txe, bdy_txe, &
2650 cbdy_tys, bdy_tys, &
2651 cbdy_tye, bdy_tye, &
2652 cdt, ndt &
2653 )
2655 ! USE module_configure , ONLY : nl_get_spec_zone, nl_get_relax_zone
2656 ! USE module_state_description
2657 USE module_interp_info
2659 IMPLICIT NONE
2661 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
2662 cims, cime, ckms, ckme, cjms, cjme, &
2663 cits, cite, ckts, ckte, cjts, cjte, &
2664 nids, nide, nkds, nkde, njds, njde, &
2665 nims, nime, nkms, nkme, njms, njme, &
2666 nits, nite, nkts, nkte, njts, njte, &
2667 shw1, & ! ignore
2668 ipos, jpos, &
2669 nri, nrj
2670 INTEGER, INTENT(IN) :: nijds, nijde, spec_bdy_width
2671 LOGICAL, INTENT(IN) :: xstag, ystag
2673 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ), INTENT(INOUT) :: cfld
2674 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ), INTENT(INOUT) :: nfld
2675 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
2676 REAL, DIMENSION ( * ), INTENT(INOUT) :: cbdy_xs, cbdy_txs ! not used
2677 REAL, DIMENSION ( * ), INTENT(INOUT) :: cbdy_xe, cbdy_txe ! not used
2678 REAL, DIMENSION ( * ), INTENT(INOUT) :: cbdy_ys, cbdy_tys ! not used
2679 REAL, DIMENSION ( * ), INTENT(INOUT) :: cbdy_ye, cbdy_tye ! not used
2680 REAL :: cdt, ndt
2681 REAL, DIMENSION ( njms:njme, nkms:nkme, spec_bdy_width ), INTENT(INOUT) :: bdy_xs, bdy_txs
2682 REAL, DIMENSION ( njms:njme, nkms:nkme, spec_bdy_width ), INTENT(INOUT) :: bdy_xe, bdy_txe
2683 REAL, DIMENSION ( nims:nime, nkms:nkme, spec_bdy_width ), INTENT(INOUT) :: bdy_ys, bdy_tys
2684 REAL, DIMENSION ( nims:nime, nkms:nkme, spec_bdy_width ), INTENT(INOUT) :: bdy_ye, bdy_tye
2686 ! Local
2688 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld_horiz_interp ! mem dimensioned on purpose
2689 ! to allow interpolating routine
2690 ! to assume this is a mem
2691 ! sized array
2692 INTEGER ni, nj, nk, istag, jstag
2693 INTEGER shw
2694 INTEGER spec_zone
2695 INTEGER relax_zone
2696 INTEGER sz
2697 REAL*8 rdt
2699 shw = 0 ! dummy, not used, but needed for the calling interface
2701 ! Horizontally interpolate the CG to the FG, store in nfld_horiz_interp
2703 CALL interp_fcn ( cfld, & ! CD field
2704 cids, cide, ckds, ckde, cjds, cjde, &
2705 cims, cime, ckms, ckme, cjms, cjme, &
2706 cits, cite, ckts, ckte, cjts, cjte, &
2707 nfld_horiz_interp, & ! ND field
2708 nids, nide, nkds, nkde, njds, njde, &
2709 nims, nime, nkms, nkme, njms, njme, &
2710 MAX(nits-nri,nids),MIN(nite+nri,nide),&
2711 nkts, nkte, &
2712 MAX(njts-nrj,njds),MIN(njte+nrj,njde),&
2713 shw, & ! stencil half width for interp
2714 imask, & ! interpolation mask
2715 xstag, ystag, & ! staggering of field
2716 ipos, jpos, & ! Position of lower left of nest in CD
2717 ! ipos-1, jpos-1, & ! Position of lower left of nest in CD
2718 nri, nrj ) ! Nest ratio, i- and j-directions
2720 ! Staggering, to determine loop indexes
2722 IF ( xstag ) THEN
2723 istag = 0
2724 ELSE
2725 istag = 1
2726 END IF
2727 IF ( ystag ) THEN
2728 jstag = 0
2729 ELSE
2730 jstag = 1
2731 END IF
2733 ! CG time step reciprocal, for computing tendencies.
2735 rdt = 1.D0/cdt
2737 CALL nl_get_spec_zone( 1, spec_zone )
2738 CALL nl_get_relax_zone( 1, relax_zone )
2740 ! Belt and suspenders ... sz is just spec_bdy_width.
2742 sz = MIN(MAX( spec_zone, relax_zone + 1 ),spec_bdy_width)
2744 !$OMP PARALLEL DO &
2745 !$OMP PRIVATE ( ni,nj,nk )
2747 DO nj = MAX ( njts-nrj, njds ) , MIN ( njte+nrj, njde-jstag )
2748 DO nk = nkts, nkte
2749 DO ni = MAX( nits-nri, nids ) , MIN ( nite+nri, nide-istag )
2751 ! WEST boundary
2753 IF ( ni .LT. nids + sz ) THEN
2754 bdy_txs(nj,nk,ni) = rdt*(nfld_horiz_interp(ni,nk,nj)-nfld(ni,nk,nj))
2755 bdy_xs (nj,nk,ni) = nfld(ni,nk,nj)
2756 END IF
2758 ! SOUTH boundary
2760 IF ( nj .LT. njds + sz ) THEN
2761 bdy_tys(ni,nk,nj) = rdt*(nfld_horiz_interp(ni,nk,nj)-nfld(ni,nk,nj))
2762 bdy_ys (ni,nk,nj) = nfld(ni,nk,nj)
2763 END IF
2765 ! EAST boundary
2767 IF ( xstag ) THEN
2768 IF ( ( ni .GE. nide - sz + 1 ) .AND. ( ni .LE. nide ) ) THEN
2769 bdy_txe(nj,nk,nide-ni+1) = rdt*(nfld_horiz_interp(ni,nk,nj)-nfld(ni,nk,nj))
2770 bdy_xe (nj,nk,nide-ni+1) = nfld(ni,nk,nj)
2771 END IF
2772 ELSE
2773 IF ( ( ni .GE. nide - sz ) .AND. ( ni .LE. nide-1 ) ) THEN
2774 bdy_txe(nj,nk,nide-ni ) = rdt*(nfld_horiz_interp(ni,nk,nj)-nfld(ni,nk,nj))
2775 bdy_xe (nj,nk,nide-ni ) = nfld(ni,nk,nj)
2776 END IF
2777 END IF
2779 ! NORTH boundary
2781 IF ( ystag ) THEN
2782 IF ( ( nj .GE. njde - sz + 1 ) .AND. ( nj .LE. njde ) ) THEN
2783 bdy_tye(ni,nk,njde-nj+1) = rdt*(nfld_horiz_interp(ni,nk,nj)-nfld(ni,nk,nj))
2784 bdy_ye (ni,nk,njde-nj+1) = nfld(ni,nk,nj)
2785 END IF
2786 ELSE
2787 IF ( ( nj .GE. njde - sz ) .AND. ( nj .LE. njde-1 ) ) THEN
2788 bdy_tye(ni,nk,njde-nj ) = rdt*(nfld_horiz_interp(ni,nk,nj)-nfld(ni,nk,nj))
2789 bdy_ye (ni,nk,njde-nj ) = nfld(ni,nk,nj)
2790 END IF
2791 END IF
2793 END DO ! nest i
2794 END DO ! nest k
2795 END DO ! nest j
2797 !$OMP END PARALLEL DO
2799 END SUBROUTINE bdy_interp2
2801 !==================================
2803 SUBROUTINE interp_fcni( cfld, & ! CD field
2804 cids, cide, ckds, ckde, cjds, cjde, &
2805 cims, cime, ckms, ckme, cjms, cjme, &
2806 cits, cite, ckts, ckte, cjts, cjte, &
2807 nfld, & ! ND field
2808 nids, nide, nkds, nkde, njds, njde, &
2809 nims, nime, nkms, nkme, njms, njme, &
2810 nits, nite, nkts, nkte, njts, njte, &
2811 shw, & ! stencil half width
2812 imask, & ! interpolation mask
2813 xstag, ystag, & ! staggering of field
2814 ipos, jpos, & ! Position of lower left of nest in CD
2815 nri, nrj ) ! nest ratios
2816 USE module_configure
2817 IMPLICIT NONE
2820 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
2821 cims, cime, ckms, ckme, cjms, cjme, &
2822 cits, cite, ckts, ckte, cjts, cjte, &
2823 nids, nide, nkds, nkde, njds, njde, &
2824 nims, nime, nkms, nkme, njms, njme, &
2825 nits, nite, nkts, nkte, njts, njte, &
2826 shw, &
2827 ipos, jpos, &
2828 nri, nrj
2829 LOGICAL, INTENT(IN) :: xstag, ystag
2831 INTEGER, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
2832 INTEGER, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
2833 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
2835 ! Local
2837 INTEGER ci, cj, ck, ni, nj, nk, ip, jp
2839 ! Iterate over the ND tile and compute the values
2840 ! from the CD tile.
2842 !write(0,'("cits:cite, ckts:ckte, cjts:cjte ",6i4)')cits,cite, ckts,ckte, cjts,cjte
2843 !write(0,'("nits:nite, nkts:nkte, njts:njte ",6i4)')nits,nite, nkts,nkte, njts,njte
2845 DO nj = njts, njte
2846 cj = jpos + (nj-1) / nrj ! j coord of CD point
2847 jp = mod ( nj , nrj ) ! coord of ND w/i CD point
2848 DO nk = nkts, nkte
2849 ck = nk
2850 DO ni = nits, nite
2851 if ( imask(ni,nj) .NE. 1 ) cycle
2852 ci = ipos + (ni-1) / nri ! j coord of CD point
2853 ip = mod ( ni , nri ) ! coord of ND w/i CD point
2854 ! This is a trivial implementation of the interp_fcn; just copies
2855 ! the values from the CD into the ND
2856 nfld( ni, nk, nj ) = cfld( ci , ck , cj )
2857 ENDDO
2858 ENDDO
2859 ENDDO
2861 RETURN
2863 END SUBROUTINE interp_fcni
2865 SUBROUTINE interp_fcnm( cfld, & ! CD field
2866 cids, cide, ckds, ckde, cjds, cjde, &
2867 cims, cime, ckms, ckme, cjms, cjme, &
2868 cits, cite, ckts, ckte, cjts, cjte, &
2869 nfld, & ! ND field
2870 nids, nide, nkds, nkde, njds, njde, &
2871 nims, nime, nkms, nkme, njms, njme, &
2872 nits, nite, nkts, nkte, njts, njte, &
2873 shw, & ! stencil half width
2874 imask, & ! interpolation mask
2875 xstag, ystag, & ! staggering of field
2876 ipos, jpos, & ! Position of lower left of nest in CD
2877 nri, nrj ) ! nest ratios
2878 USE module_configure
2879 IMPLICIT NONE
2882 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
2883 cims, cime, ckms, ckme, cjms, cjme, &
2884 cits, cite, ckts, ckte, cjts, cjte, &
2885 nids, nide, nkds, nkde, njds, njde, &
2886 nims, nime, nkms, nkme, njms, njme, &
2887 nits, nite, nkts, nkte, njts, njte, &
2888 shw, &
2889 ipos, jpos, &
2890 nri, nrj
2891 LOGICAL, INTENT(IN) :: xstag, ystag
2893 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
2894 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
2895 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
2897 ! Local
2899 INTEGER ci, cj, ck, ni, nj, nk, ip, jp
2901 ! Iterate over the ND tile and compute the values
2902 ! from the CD tile.
2904 !write(0,'("mask cits:cite, ckts:ckte, cjts:cjte ",6i4)')cits,cite, ckts,ckte, cjts,cjte
2905 !write(0,'("mask nits:nite, nkts:nkte, njts:njte ",6i4)')nits,nite, nkts,nkte, njts,njte
2907 DO nj = njts, njte
2908 cj = jpos + (nj-1) / nrj ! j coord of CD point
2909 jp = mod ( nj , nrj ) ! coord of ND w/i CD point
2910 DO nk = nkts, nkte
2911 ck = nk
2912 DO ni = nits, nite
2913 ci = ipos + (ni-1) / nri ! j coord of CD point
2914 ip = mod ( ni , nri ) ! coord of ND w/i CD point
2915 ! This is a trivial implementation of the interp_fcn; just copies
2916 ! the values from the CD into the ND
2917 nfld( ni, nk, nj ) = cfld( ci , ck , cj )
2918 ENDDO
2919 ENDDO
2920 ENDDO
2922 RETURN
2924 END SUBROUTINE interp_fcnm
2926 SUBROUTINE interp_fcnm_lu( cfld, & ! CD field
2927 cids, cide, ckds, ckde, cjds, cjde, &
2928 cims, cime, ckms, ckme, cjms, cjme, &
2929 cits, cite, ckts, ckte, cjts, cjte, &
2930 nfld, & ! ND field
2931 nids, nide, nkds, nkde, njds, njde, &
2932 nims, nime, nkms, nkme, njms, njme, &
2933 nits, nite, nkts, nkte, njts, njte, &
2934 shw, & ! stencil half width
2935 imask, & ! interpolation mask
2936 xstag, ystag, & ! staggering of field
2937 ipos, jpos, & ! Position of lower left of nest in CD
2938 nri, nrj, & ! nest ratios
2939 cxlat, nxlat, &
2940 cxlong, nxlong, &
2941 cdx, ndx, &
2942 cid, nid )
2943 USE module_configure
2945 IMPLICIT NONE
2948 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
2949 cims, cime, ckms, ckme, cjms, cjme, &
2950 cits, cite, ckts, ckte, cjts, cjte, &
2951 nids, nide, nkds, nkde, njds, njde, &
2952 nims, nime, nkms, nkme, njms, njme, &
2953 nits, nite, nkts, nkte, njts, njte, &
2954 shw, &
2955 ipos, jpos, &
2956 nri, nrj, &
2957 cid, nid
2958 LOGICAL, INTENT(IN) :: xstag, ystag
2960 REAL, INTENT(IN) :: cdx, ndx
2962 REAL, INTENT(IN), DIMENSION ( cims:cime, cjms:cjme ) :: cxlat, cxlong
2963 REAL, INTENT(IN), DIMENSION ( nims:nime, njms:njme ) :: nxlat, nxlong
2966 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
2967 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
2968 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
2970 ! Local
2972 INTEGER i, ci, cj, ck, ni, nj, nk, ip, jp, ierr
2974 #ifdef TERRAIN_AND_LANDUSE
2975 INTEGER, DIMENSION(256) :: ipath ! array for integer coded ascii for passing path down to get_landuse
2977 REAL , ALLOCATABLE, DIMENSION(:,:) :: xlat_g, xlon_g, landuse_g
2978 CHARACTER*256 :: message
2979 CHARACTER*256 :: rsmas_data_path
2981 LOGICAL :: input_from_hires, input_from_file
2983 INTEGER, EXTERNAL :: get_landuse
2984 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
2986 CALL nl_get_input_from_hires( nid , input_from_hires)
2987 CALL nl_get_input_from_file ( nid , input_from_file )
2989 IF ( input_from_file .AND. input_from_hires ) THEN
2990 Write(message, '(a,i3,a)') &
2991 "Warning : input_from_file turned on for domain ", nid, ", input_from_hires disabled"
2992 CALL wrf_message(message)
2993 END IF
2995 IF ( .NOT. input_from_file .AND. input_from_hires ) THEN
2997 allocate(xlat_g(nids:nide,njds:njde))
2998 allocate(xlon_g(nids:nide,njds:njde))
2999 allocate(landuse_g(nids:nide,njds:njde))
3001 CALL nl_get_rsmas_data_path(1,rsmas_data_path)
3003 DO i = 1, LEN(TRIM(rsmas_data_path))
3004 ipath(i) = ICHAR(rsmas_data_path(i:i))
3005 ENDDO
3007 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
3008 CALL wrf_patch_to_global_real ( nxlat, xlat_g , nid, ' ' , 'xy' , &
3009 nids, nide-1 , njds , njde-1 , 1 , 1 , &
3010 nims, nime , njms , njme , 1 , 1 , &
3011 nits, nite , njts , njte , 1 , 1 )
3012 CALL wrf_patch_to_global_real ( nxlong, xlon_g, nid, ' ' , 'xy' , &
3013 nids, nide-1 , njds , njde-1 , 1 , 1 , &
3014 nims, nime , njms , njme , 1 , 1 , &
3015 nits, nite , njts , njte , 1 , 1 )
3016 IF ( wrf_dm_on_monitor() ) THEN
3017 ierr = get_landuse ( ndx/1000., xlat_g, xlon_g, &
3018 landuse_g, &
3019 nide-nids+1,njde-njds+1,nide-nids+1,njde-njds+1, &
3020 ipath, LEN(TRIM(rsmas_data_path)) )
3021 IF ( ierr == 1 ) THEN
3022 WRITE(message,fmt='(a)') 'get_landuse : aborted!'
3023 CALL wrf_error_fatal(TRIM(message))
3024 ENDIF
3025 ENDIF
3027 CALL wrf_global_to_patch_real ( landuse_g , nfld(:,1,:), nid, ' ' , 'xy' , &
3028 nids, nide-1 , njds , njde-1 , 1 , 1 , &
3029 nims, nime , njms , njme , 1 , 1 , &
3030 nits, nite , njts , njte , 1 , 1 )
3032 #else
3033 ierr = get_landuse ( ndx/1000., nxlat(nids:nide,njds:njde), nxlong(nids:nide,njds:njde), &
3034 nfld(nids:nide,1,njds:njde), &
3035 nide-nids+1,njde-njds+1,nide-nids+1,njde-njds+1, &
3036 ipath, LEN(TRIM(rsmas_data_path)) )
3037 #endif
3038 deallocate(xlat_g)
3039 deallocate(xlon_g)
3040 deallocate(landuse_g)
3041 ELSE
3042 #endif
3043 ! Iterate over the ND tile and compute the values
3044 ! from the CD tile.
3045 DO nj = njts, njte
3046 cj = jpos + (nj-1) / nrj ! j coord of CD point
3047 jp = mod ( nj , nrj ) ! coord of ND w/i CD point
3048 DO nk = nkts, nkte
3049 ck = nk
3050 DO ni = nits, nite
3051 ci = ipos + (ni-1) / nri ! j coord of CD point
3052 ip = mod ( ni , nri ) ! coord of ND w/i CD point
3053 ! This is a trivial implementation of the interp_fcn; just copies
3054 ! the values from the CD into the ND
3055 if ( imask(ni,nj) .eq. 1 ) then
3056 nfld( ni, nk, nj ) = cfld( ci , ck , cj )
3057 endif
3058 ENDDO
3059 ENDDO
3060 ENDDO
3061 #ifdef TERRAIN_AND_LANDUSE
3062 END IF
3063 #endif
3064 RETURN
3066 END SUBROUTINE interp_fcnm_lu
3069 SUBROUTINE interp_fcnm_imask( cfld, & ! CD field
3070 cids, cide, ckds, ckde, cjds, cjde, &
3071 cims, cime, ckms, ckme, cjms, cjme, &
3072 cits, cite, ckts, ckte, cjts, cjte, &
3073 nfld, & ! ND field
3074 nids, nide, nkds, nkde, njds, njde, &
3075 nims, nime, nkms, nkme, njms, njme, &
3076 nits, nite, nkts, nkte, njts, njte, &
3077 shw, & ! stencil half width
3078 imask, & ! interpolation mask
3079 xstag, ystag, & ! staggering of field
3080 ipos, jpos, & ! Position of lower left of nest in CD
3081 nri, nrj ) ! nest ratios
3082 USE module_configure
3083 IMPLICIT NONE
3086 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3087 cims, cime, ckms, ckme, cjms, cjme, &
3088 cits, cite, ckts, ckte, cjts, cjte, &
3089 nids, nide, nkds, nkde, njds, njde, &
3090 nims, nime, nkms, nkme, njms, njme, &
3091 nits, nite, nkts, nkte, njts, njte, &
3092 shw, &
3093 ipos, jpos, &
3094 nri, nrj
3095 LOGICAL, INTENT(IN) :: xstag, ystag
3097 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
3098 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
3099 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
3101 ! Local
3103 INTEGER ci, cj, ck, ni, nj, nk, ip, jp
3105 ! Iterate over the ND tile and compute the values
3106 ! from the CD tile.
3108 !write(0,'("mask cits:cite, ckts:ckte, cjts:cjte ",6i4)')cits,cite, ckts,ckte,cjts,cjte
3109 !write(0,'("mask nits:nite, nkts:nkte, njts:njte ",6i4)')nits,nite, nkts,nkte,njts,njte
3111 DO nj = njts, njte
3112 cj = jpos + (nj-1) / nrj ! j coord of CD point
3113 jp = mod ( nj , nrj ) ! coord of ND w/i CD point
3114 DO nk = nkts, nkte
3115 ck = nk
3116 DO ni = nits, nite
3117 ci = ipos + (ni-1) / nri ! j coord of CD point
3118 ip = mod ( ni , nri ) ! coord of ND w/i CD point
3119 ! This is a trivial implementation of the interp_fcn; just copies
3120 ! the values from the CD into the ND
3121 if ( imask(ni,nj) .eq. 1 ) then
3122 nfld( ni, nk, nj ) = cfld( ci , ck , cj )
3123 endif
3124 ENDDO
3125 ENDDO
3126 ENDDO
3128 RETURN
3130 END SUBROUTINE interp_fcnm_imask
3131 ! end of first block of ARW-only routines
3133 SUBROUTINE interp_mask_land_field ( enable, & ! says whether to allow interpolation or just the bcasts
3134 cfld, & ! CD field
3135 cids, cide, ckds, ckde, cjds, cjde, &
3136 cims, cime, ckms, ckme, cjms, cjme, &
3137 cits, cite, ckts, ckte, cjts, cjte, &
3138 nfld, & ! ND field
3139 nids, nide, nkds, nkde, njds, njde, &
3140 nims, nime, nkms, nkme, njms, njme, &
3141 nits, nite, nkts, nkte, njts, njte, &
3142 shw, & ! stencil half width
3143 imask, & ! interpolation mask
3144 xstag, ystag, & ! staggering of field
3145 ipos, jpos, & ! Position of lower left of nest in CD
3146 nri, nrj, & ! nest ratios
3147 clu, nlu )
3149 USE module_configure
3150 USE module_wrf_error
3151 USE module_dm, only : wrf_dm_sum_reals, wrf_dm_sum_integers
3153 IMPLICIT NONE
3156 LOGICAL, INTENT(IN) :: enable
3157 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3158 cims, cime, ckms, ckme, cjms, cjme, &
3159 cits, cite, ckts, ckte, cjts, cjte, &
3160 nids, nide, nkds, nkde, njds, njde, &
3161 nims, nime, nkms, nkme, njms, njme, &
3162 nits, nite, nkts, nkte, njts, njte, &
3163 shw, &
3164 ipos, jpos, &
3165 nri, nrj
3166 LOGICAL, INTENT(IN) :: xstag, ystag
3168 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
3169 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
3170 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
3172 REAL, DIMENSION ( cims:cime, cjms:cjme ) :: clu
3173 REAL, DIMENSION ( nims:nime, njms:njme ) :: nlu
3175 ! Local
3177 INTEGER ci, cj, ck, ni, nj, nk, ip, jp
3178 INTEGER :: icount , ii , jj , ist , ien , jst , jen , iswater, ierr
3179 REAL :: avg , sum , dx , dy
3180 INTEGER , PARAMETER :: max_search = 5
3181 CHARACTER(LEN=255) :: message
3182 INTEGER :: icount_n(nkts:nkte), idummy(nkts:nkte)
3183 REAL :: avg_n(nkts:nkte), sum_n(nkts:nkte), dummy(nkts:nkte)
3185 ! Find out what the water value is.
3187 CALL nl_get_iswater(1,iswater)
3189 ! Right now, only mass point locations permitted.
3191 IF ( ( .NOT. xstag) .AND. ( .NOT. ystag ) ) THEN
3193 ! Loop over each i,k,j in the nested domain.
3195 IF ( enable ) THEN
3197 DO nj = njts, njte
3198 IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
3199 cj = ( nj + (nrj/2)-1 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
3200 ELSE
3201 cj = ( nj + (nrj-1)/2 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
3202 END IF
3203 DO nk = nkts, nkte
3204 ck = nk
3205 DO ni = nits, nite
3206 IF ( imask(ni, nj) .NE. 1 ) cycle
3207 IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
3208 ci = ( ni + (nri/2)-1 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
3209 ELSE
3210 ci = ( ni + (nri-1)/2 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
3211 END IF
3216 !
3217 ! (ci,cj+1) (ci+1,cj+1)
3218 ! - -------------
3219 ! 1-dy | | |
3220 ! | | |
3221 ! - | * |
3222 ! dy | | (ni,nj) |
3223 ! | | |
3224 ! - -------------
3225 ! (ci,cj) (ci+1,cj)
3226 !
3227 ! |--|--------|
3228 ! dx 1-dx
3231 ! For odd ratios, at (nri+1)/2, we are on the coarse grid point, so dx = 0
3233 IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
3234 dx = ( REAL ( MOD ( ni+(nri-1)/2 , nri ) ) + 0.5 ) / REAL ( nri )
3235 ELSE
3236 dx = REAL ( MOD ( ni+(nri-1)/2 , nri ) ) / REAL ( nri )
3237 END IF
3238 IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
3239 dy = ( REAL ( MOD ( nj+(nrj-1)/2 , nrj ) ) + 0.5 ) / REAL ( nrj )
3240 ELSE
3241 dy = REAL ( MOD ( nj+(nrj-1)/2 , nrj ) ) / REAL ( nrj )
3242 END IF
3244 ! This is a "land only" field. If this is a water point, no operations required.
3246 IF ( ( NINT(nlu(ni ,nj )) .EQ. iswater ) ) THEN
3247 nfld(ni,nk,nj) = cfld(ci ,ck,cj )
3249 ! If this is a nested land point, and the surrounding coarse values are all land points,
3250 ! then this is a simple 4-pt interpolation.
3252 ELSE IF ( ( NINT(nlu(ni ,nj )) .NE. iswater ) .AND. &
3253 ( NINT(clu(ci ,cj )) .NE. iswater ) .AND. &
3254 ( NINT(clu(ci+1,cj )) .NE. iswater ) .AND. &
3255 ( NINT(clu(ci ,cj+1)) .NE. iswater ) .AND. &
3256 ( NINT(clu(ci+1,cj+1)) .NE. iswater ) ) THEN
3257 nfld(ni,nk,nj) = ( 1. - dx ) * ( ( 1. - dy ) * cfld(ci ,ck,cj ) + &
3258 dy * cfld(ci ,ck,cj+1) ) + &
3259 dx * ( ( 1. - dy ) * cfld(ci+1,ck,cj ) + &
3260 dy * cfld(ci+1,ck,cj+1) )
3262 ! If this is a nested land point and there are NO coarse land values surrounding,
3263 ! we temporarily punt.
3265 ELSE IF ( ( NINT(nlu(ni ,nj )) .NE. iswater ) .AND. &
3266 ( NINT(clu(ci ,cj )) .EQ. iswater ) .AND. &
3267 ( NINT(clu(ci+1,cj )) .EQ. iswater ) .AND. &
3268 ( NINT(clu(ci ,cj+1)) .EQ. iswater ) .AND. &
3269 ( NINT(clu(ci+1,cj+1)) .EQ. iswater ) ) THEN
3270 nfld(ni,nk,nj) = -1
3272 ! If there are some water points and some land points, take an average.
3274 ELSE IF ( NINT(nlu(ni ,nj )) .NE. iswater ) THEN
3275 icount = 0
3276 sum = 0
3277 IF ( NINT(clu(ci ,cj )) .NE. iswater ) THEN
3278 icount = icount + 1
3279 sum = sum + cfld(ci ,ck,cj )
3280 END IF
3281 IF ( NINT(clu(ci+1,cj )) .NE. iswater ) THEN
3282 icount = icount + 1
3283 sum = sum + cfld(ci+1,ck,cj )
3284 END IF
3285 IF ( NINT(clu(ci ,cj+1)) .NE. iswater ) THEN
3286 icount = icount + 1
3287 sum = sum + cfld(ci ,ck,cj+1)
3288 END IF
3289 IF ( NINT(clu(ci+1,cj+1)) .NE. iswater ) THEN
3290 icount = icount + 1
3291 sum = sum + cfld(ci+1,ck,cj+1)
3292 END IF
3293 nfld(ni,nk,nj) = sum / REAL ( icount )
3294 END IF
3295 END DO
3296 END DO
3297 END DO
3300 ! Get an average of the whole domain for problem locations.
3302 sum_n = 0
3303 icount_n = 0
3304 DO nj = njts, njte
3305 DO nk = nkts, nkte
3306 DO ni = nits, nite
3307 IF ( nfld(ni,nk,nj) .NE. -1 ) THEN
3308 IF ( NINT(nlu(ni,nj)) .NE. iswater ) THEN
3309 icount_n(nk) = icount_n(nk) + 1
3310 sum_n(nk) = sum_n(nk) + nfld(ni,nk,nj)
3311 END IF
3312 END IF
3313 END DO
3314 END DO
3315 END DO
3317 CALL wrf_dm_sum_reals( sum_n(nkts:nkte), dummy(nkts:nkte))
3318 sum_n = dummy
3319 CALL wrf_dm_sum_integers(icount_n(nkts:nkte), idummy(nkts:nkte))
3320 icount_n = idummy
3321 DO nk = nkts, nkte
3322 IF ( icount_n(nk) .GT. 0 ) &
3323 avg_n(nk) = sum_n(nk) / icount_n(nk)
3324 END DO
3325 ENDIF
3327 IF ( enable ) THEN
3328 IF ( ANY(nfld .EQ. -1) ) THEN
3330 ! OK, if there were any of those island situations, we try to search a bit broader
3331 ! of an area in the coarse grid.
3333 DO nj = njts, njte
3334 DO nk = nkts, nkte
3335 DO ni = nits, nite
3336 IF ( imask(ni, nj) .NE. 1 ) cycle
3337 IF ( nfld(ni,nk,nj) .EQ. -1 ) THEN
3338 IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
3339 cj = ( nj + (nrj/2)-1 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
3340 ELSE
3341 cj = ( nj + (nrj-1)/2 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
3342 END IF
3343 IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
3344 ci = ( ni + (nri/2)-1 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
3345 ELSE
3346 ci = ( ni + (nri-1)/2 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
3347 END IF
3348 ist = MAX (ci-max_search,cits)
3349 ien = MIN (ci+max_search,cite,cide-1)
3350 jst = MAX (cj-max_search,cjts)
3351 jen = MIN (cj+max_search,cjte,cjde-1)
3352 icount = 0
3353 sum = 0
3354 DO jj = jst,jen
3355 DO ii = ist,ien
3356 IF ( NINT(clu(ii,jj)) .NE. iswater ) THEN
3357 icount = icount + 1
3358 sum = sum + cfld(ii,nk,jj)
3359 END IF
3360 END DO
3361 END DO
3362 IF ( icount .GT. 0 ) THEN
3363 nfld(ni,nk,nj) = sum / REAL ( icount )
3364 ELSE
3365 Write(message,fmt='(a,i4,a,i4,a,f10.4)') &
3366 'horizontal interp error - island (', ni, ',', nj, '), using average ', avg_n(nk)
3367 CALL wrf_message ( message )
3368 nfld(ni,nk,nj) = avg_n(nk)
3369 END IF
3370 END IF
3371 END DO
3372 END DO
3373 END DO
3374 ENDIF
3375 ENDIF
3376 ELSE
3377 CALL wrf_error_fatal ( "only unstaggered fields right now" )
3378 END IF
3380 END SUBROUTINE interp_mask_land_field
3382 SUBROUTINE interp_mask_water_field ( enable, & ! says whether to allow interpolation or just the bcasts
3383 cfld, & ! CD field
3384 cids, cide, ckds, ckde, cjds, cjde, &
3385 cims, cime, ckms, ckme, cjms, cjme, &
3386 cits, cite, ckts, ckte, cjts, cjte, &
3387 nfld, & ! ND field
3388 nids, nide, nkds, nkde, njds, njde, &
3389 nims, nime, nkms, nkme, njms, njme, &
3390 nits, nite, nkts, nkte, njts, njte, &
3391 shw, & ! stencil half width
3392 imask, & ! interpolation mask
3393 xstag, ystag, & ! staggering of field
3394 ipos, jpos, & ! Position of lower left of nest in CD
3395 nri, nrj, & ! nest ratios
3396 clu, nlu, cflag, nflag )
3398 USE module_configure
3399 USE module_wrf_error
3400 USE module_dm, only : wrf_dm_sum_reals, wrf_dm_sum_integers
3402 IMPLICIT NONE
3405 LOGICAL, INTENT(IN) :: enable
3406 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3407 cims, cime, ckms, ckme, cjms, cjme, &
3408 cits, cite, ckts, ckte, cjts, cjte, &
3409 nids, nide, nkds, nkde, njds, njde, &
3410 nims, nime, nkms, nkme, njms, njme, &
3411 nits, nite, nkts, nkte, njts, njte, &
3412 shw, &
3413 ipos, jpos, &
3414 nri, nrj, cflag, nflag
3415 LOGICAL, INTENT(IN) :: xstag, ystag
3417 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
3418 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
3419 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
3421 REAL, DIMENSION ( cims:cime, cjms:cjme ) :: clu
3422 REAL, DIMENSION ( nims:nime, njms:njme ) :: nlu
3424 ! Local
3426 INTEGER ci, cj, ck, ni, nj, nk, ip, jp
3427 INTEGER :: icount , ii , jj , ist , ien , jst , jen, ierr
3428 REAL :: avg , sum , dx , dy
3429 INTEGER , PARAMETER :: max_search = 5
3430 INTEGER :: icount_n(nkts:nkte), idummy(nkts:nkte)
3431 REAL :: avg_n(nkts:nkte), sum_n(nkts:nkte), dummy(nkts:nkte)
3432 CHARACTER(LEN=255) :: message
3434 ! Right now, only mass point locations permitted.
3436 IF ( ( .NOT. xstag) .AND. ( .NOT. ystag ) ) THEN
3438 IF ( enable ) THEN
3439 ! Loop over each i,k,j in the nested domain.
3441 DO nj = njts, njte
3442 IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
3443 cj = ( nj + (nrj/2)-1 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
3444 ELSE
3445 cj = ( nj + (nrj-1)/2 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
3446 END IF
3447 DO nk = nkts, nkte
3448 ck = nk
3449 DO ni = nits, nite
3450 !dave IF ( imask(ni, nj) .NE. 1 ) cycle
3451 IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
3452 ci = ( ni + (nri/2)-1 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
3453 ELSE
3454 ci = ( ni + (nri-1)/2 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
3455 END IF
3460 !
3461 ! (ci,cj+1) (ci+1,cj+1)
3462 ! - -------------
3463 ! 1-dy | | |
3464 ! | | |
3465 ! - | * |
3466 ! dy | | (ni,nj) |
3467 ! | | |
3468 ! - -------------
3469 ! (ci,cj) (ci+1,cj)
3470 !
3471 ! |--|--------|
3472 ! dx 1-dx
3475 ! At ni=2, we are on the coarse grid point, so dx = 0
3477 IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
3478 dx = ( REAL ( MOD ( ni+(nri-1)/2 , nri ) ) + 0.5 ) / REAL ( nri )
3479 ELSE
3480 dx = REAL ( MOD ( ni+(nri-1)/2 , nri ) ) / REAL ( nri )
3481 END IF
3482 IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
3483 dy = ( REAL ( MOD ( nj+(nrj-1)/2 , nrj ) ) + 0.5 ) / REAL ( nrj )
3484 ELSE
3485 dy = REAL ( MOD ( nj+(nrj-1)/2 , nrj ) ) / REAL ( nrj )
3486 END IF
3488 ! This is a "water only" field. If this is a land point, no operations required.
3490 IF ( ( NINT(nlu(ni ,nj )) .NE. nflag ) ) THEN
3491 nfld(ni,nk,nj) = cfld(ci ,ck,cj )
3493 ! If this is a nested water point, and the surrounding coarse values are all water points,
3494 ! then this is a simple 4-pt interpolation.
3496 ELSE IF ( ( NINT(nlu(ni ,nj )) .EQ. nflag ) .AND. &
3497 ( NINT(clu(ci ,cj )) .EQ. nflag ) .AND. &
3498 ( NINT(clu(ci+1,cj )) .EQ. nflag ) .AND. &
3499 ( NINT(clu(ci ,cj+1)) .EQ. nflag ) .AND. &
3500 ( NINT(clu(ci+1,cj+1)) .EQ. nflag ) ) THEN
3501 nfld(ni,nk,nj) = ( 1. - dx ) * ( ( 1. - dy ) * cfld(ci ,ck,cj ) + &
3502 dy * cfld(ci ,ck,cj+1) ) + &
3503 dx * ( ( 1. - dy ) * cfld(ci+1,ck,cj ) + &
3504 dy * cfld(ci+1,ck,cj+1) )
3506 ! If this is a nested water point and there are NO coarse water values surrounding,
3507 ! we temporarily punt.
3509 ELSE IF ( ( NINT(nlu(ni ,nj )) .EQ. nflag ) .AND. &
3510 ( NINT(clu(ci ,cj )) .NE. nflag ) .AND. &
3511 ( NINT(clu(ci+1,cj )) .NE. nflag ) .AND. &
3512 ( NINT(clu(ci ,cj+1)) .NE. nflag ) .AND. &
3513 ( NINT(clu(ci+1,cj+1)) .NE. nflag ) ) THEN
3514 nfld(ni,nk,nj) = -4
3516 ! If there are some land points and some water points, take an average.
3518 ELSE IF ( NINT(nlu(ni ,nj )) .EQ. nflag ) THEN
3519 icount = 0
3520 sum = 0
3521 IF ( NINT(clu(ci ,cj )) .EQ. nflag ) THEN
3522 icount = icount + 1
3523 sum = sum + cfld(ci ,ck,cj )
3524 END IF
3525 IF ( NINT(clu(ci+1,cj )) .EQ. nflag ) THEN
3526 icount = icount + 1
3527 sum = sum + cfld(ci+1,ck,cj )
3528 END IF
3529 IF ( NINT(clu(ci ,cj+1)) .EQ. nflag ) THEN
3530 icount = icount + 1
3531 sum = sum + cfld(ci ,ck,cj+1)
3532 END IF
3533 IF ( NINT(clu(ci+1,cj+1)) .EQ. nflag ) THEN
3534 icount = icount + 1
3535 sum = sum + cfld(ci+1,ck,cj+1)
3536 END IF
3537 nfld(ni,nk,nj) = sum / REAL ( icount )
3538 END IF
3539 END DO
3540 END DO
3541 END DO
3543 ! Get an average of the whole domain for problem locations.
3545 sum_n = 0
3546 icount_n = 0
3547 DO nj = njts, njte
3548 DO nk = nkts, nkte
3549 DO ni = nits, nite
3550 IF ( nfld(ni,nk,nj) .NE. -1 ) THEN
3551 IF ( NINT(nlu(ni,nj)) .EQ. nflag ) THEN
3552 icount_n(nk) = icount_n(nk) + 1
3553 sum_n(nk) = sum_n(nk) + nfld(ni,nk,nj)
3554 END IF
3555 END IF
3556 END DO
3557 END DO
3558 END DO
3560 CALL wrf_dm_sum_reals( sum_n(nkts:nkte), dummy(nkts:nkte))
3561 sum_n = dummy
3562 CALL wrf_dm_sum_integers(icount_n(nkts:nkte), idummy(nkts:nkte))
3563 icount_n = idummy
3564 DO nk = nkts, nkte
3565 IF ( icount_n(nk) .GT. 0 ) &
3566 avg_n(nk) = sum_n(nk) / icount_n(nk)
3567 END DO
3568 ENDIF
3570 IF ( enable ) THEN
3571 IF ( ANY(nfld .EQ. -4) ) THEN
3573 ! OK, if there were any of those lake situations, we try to search a bit broader
3574 ! of an area in the coarse grid.
3576 DO nj = njts, njte
3577 DO nk = nkts, nkte
3578 DO ni = nits, nite
3579 !dave IF ( imask(ni, nj) .NE. 1 ) cycle
3580 IF ( nfld(ni,nk,nj) .EQ. -4 ) THEN
3581 IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
3582 cj = ( nj + (nrj/2)-1 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
3583 ELSE
3584 cj = ( nj + (nrj-1)/2 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
3585 END IF
3586 IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
3587 ci = ( ni + (nri/2)-1 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
3588 ELSE
3589 ci = ( ni + (nri-1)/2 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
3590 END IF
3591 ist = MAX (ci-max_search,cits)
3592 ien = MIN (ci+max_search,cite,cide-1)
3593 jst = MAX (cj-max_search,cjts)
3594 jen = MIN (cj+max_search,cjte,cjde-1)
3595 icount = 0
3596 sum = 0
3597 DO jj = jst,jen
3598 DO ii = ist,ien
3599 IF ( NINT(clu(ii,jj)) .EQ. nflag ) THEN
3600 icount = icount + 1
3601 sum = sum + cfld(ii,nk,jj)
3602 END IF
3603 END DO
3604 END DO
3605 IF ( icount .GT. 0 ) THEN
3606 nfld(ni,nk,nj) = sum / REAL ( icount )
3607 ELSE
3608 Write(message,fmt='(a,i4,a,i4,a,f10.4)') &
3609 'horizontal interp error - lake (', ni, ',', nj, '), using average ', avg_n(nk)
3610 CALL wrf_message ( message )
3611 nfld(ni,nk,nj) = avg_n(nk)
3612 END IF
3613 END IF
3614 END DO
3615 END DO
3616 END DO
3617 ENDIF
3618 ENDIF
3619 ELSE
3620 CALL wrf_error_fatal ( "only unstaggered fields right now" )
3621 END IF
3623 END SUBROUTINE interp_mask_water_field
3625 SUBROUTINE p2c_mask ( cfld, & ! CD field
3626 cids, cide, ckds, ckde, cjds, cjde, &
3627 cims, cime, ckms, ckme, cjms, cjme, &
3628 cits, cite, ckts, ckte, cjts, cjte, &
3629 nfld, & ! ND field
3630 nids, nide, nkds, nkde, njds, njde, &
3631 nims, nime, nkms, nkme, njms, njme, &
3632 nits, nite, nkts, nkte, njts, njte, &
3633 shw, & ! stencil half width
3634 imask, & ! interpolation mask
3635 xstag, ystag, & ! staggering of field
3636 ipos, jpos, & ! Position of lower left of nest in CD
3637 nri, nrj, & ! nest ratios
3638 clu, nlu, & ! land use categories
3639 ctslb,ntslb, & ! soil temps
3640 cnum_soil_layers,nnum_soil_layers, & ! number of soil layers for tslb
3641 ciswater, niswater ) ! iswater category
3643 USE module_configure
3644 USE module_wrf_error
3646 IMPLICIT NONE
3648 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3649 cims, cime, ckms, ckme, cjms, cjme, &
3650 cits, cite, ckts, ckte, cjts, cjte, &
3651 nids, nide, nkds, nkde, njds, njde, &
3652 nims, nime, nkms, nkme, njms, njme, &
3653 nits, nite, nkts, nkte, njts, njte, &
3654 shw, &
3655 ipos, jpos, &
3656 nri, nrj, &
3657 cnum_soil_layers, nnum_soil_layers, &
3658 ciswater, niswater
3660 LOGICAL, INTENT(IN) :: xstag, ystag
3662 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
3663 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
3664 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
3666 REAL, DIMENSION ( cims:cime, cjms:cjme ) :: clu
3667 REAL, DIMENSION ( nims:nime, njms:njme ) :: nlu
3669 REAL, DIMENSION ( cims:cime, 1:cnum_soil_layers, cjms:cjme ) :: ctslb
3670 REAL, DIMENSION ( nims:nime, 1:nnum_soil_layers, njms:njme ) :: ntslb
3672 ! Local
3674 INTEGER ci, cj, ck, ni, nj, nk
3675 INTEGER :: icount
3676 REAL :: sum , dx , dy
3678 ! Right now, only mass point locations permitted.
3680 IF ( ( .NOT. xstag) .AND. ( .NOT. ystag ) ) THEN
3682 ! Loop over each i,k,j in the nested domain.
3684 DO nj = njts, MIN(njde-1,njte)
3685 IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
3686 cj = ( nj + (nrj/2)-1 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
3687 ELSE
3688 cj = ( nj + (nrj-1)/2 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
3689 END IF
3690 DO nk = nkts, nkte
3691 ck = nk
3692 DO ni = nits, MIN(nide-1,nite)
3693 IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
3694 ci = ( ni + (nri/2)-1 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
3695 ELSE
3696 ci = ( ni + (nri-1)/2 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
3697 END IF
3699 !
3700 ! (ci,cj+1) (ci+1,cj+1)
3701 ! - -------------
3702 ! 1-dy | | |
3703 ! | | |
3704 ! - | * |
3705 ! dy | | (ni,nj) |
3706 ! | | |
3707 ! - -------------
3708 ! (ci,cj) (ci+1,cj)
3709 !
3710 ! |--|--------|
3711 ! dx 1-dx
3714 ! At ni=2, we are on the coarse grid point, so dx = 0
3716 IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
3717 dx = ( REAL ( MOD ( ni+(nri-1)/2 , nri ) ) + 0.5 ) / REAL ( nri )
3718 ELSE
3719 dx = REAL ( MOD ( ni+(nri-1)/2 , nri ) ) / REAL ( nri )
3720 END IF
3721 IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
3722 dy = ( REAL ( MOD ( nj+(nrj-1)/2 , nrj ) ) + 0.5 ) / REAL ( nrj )
3723 ELSE
3724 dy = REAL ( MOD ( nj+(nrj-1)/2 , nrj ) ) / REAL ( nrj )
3725 END IF
3727 ! This is a "water only" field. If this is a land point, no operations required.
3729 IF ( ( NINT(nlu(ni ,nj )) .NE. niswater ) ) THEN
3730 nfld(ni,nk,nj) = 273.18
3732 ! If this is a nested water point, and the surrounding coarse values are all water points,
3733 ! then this is a simple 4-pt interpolation.
3735 ELSE IF ( ( NINT(nlu(ni ,nj )) .EQ. niswater ) .AND. &
3736 ( NINT(clu(ci ,cj )) .EQ. niswater ) .AND. &
3737 ( NINT(clu(ci+1,cj )) .EQ. niswater ) .AND. &
3738 ( NINT(clu(ci ,cj+1)) .EQ. niswater ) .AND. &
3739 ( NINT(clu(ci+1,cj+1)) .EQ. niswater ) ) THEN
3740 nfld(ni,nk,nj) = ( 1. - dx ) * ( ( 1. - dy ) * cfld(ci ,ck,cj ) + &
3741 dy * cfld(ci ,ck,cj+1) ) + &
3742 dx * ( ( 1. - dy ) * cfld(ci+1,ck,cj ) + &
3743 dy * cfld(ci+1,ck,cj+1) )
3745 ! If this is a nested water point and there are NO coarse water values surrounding,
3746 ! we manufacture something from the deepest CG soil temp.
3748 ELSE IF ( ( NINT(nlu(ni ,nj )) .EQ. niswater ) .AND. &
3749 ( NINT(clu(ci ,cj )) .NE. niswater ) .AND. &
3750 ( NINT(clu(ci+1,cj )) .NE. niswater ) .AND. &
3751 ( NINT(clu(ci ,cj+1)) .NE. niswater ) .AND. &
3752 ( NINT(clu(ci+1,cj+1)) .NE. niswater ) ) THEN
3753 nfld(ni,nk,nj) = ( 1. - dx ) * ( ( 1. - dy ) * ctslb(ci ,cnum_soil_layers,cj ) + &
3754 dy * ctslb(ci ,cnum_soil_layers,cj+1) ) + &
3755 dx * ( ( 1. - dy ) * ctslb(ci+1,cnum_soil_layers,cj ) + &
3756 dy * ctslb(ci+1,cnum_soil_layers,cj+1) )
3758 ! If there are some land points and some water points, take an average of the water points.
3760 ELSE IF ( NINT(nlu(ni ,nj )) .EQ. niswater ) THEN
3761 icount = 0
3762 sum = 0
3763 IF ( NINT(clu(ci ,cj )) .EQ. niswater ) THEN
3764 icount = icount + 1
3765 sum = sum + cfld(ci ,ck,cj )
3766 END IF
3767 IF ( NINT(clu(ci+1,cj )) .EQ. niswater ) THEN
3768 icount = icount + 1
3769 sum = sum + cfld(ci+1,ck,cj )
3770 END IF
3771 IF ( NINT(clu(ci ,cj+1)) .EQ. niswater ) THEN
3772 icount = icount + 1
3773 sum = sum + cfld(ci ,ck,cj+1)
3774 END IF
3775 IF ( NINT(clu(ci+1,cj+1)) .EQ. niswater ) THEN
3776 icount = icount + 1
3777 sum = sum + cfld(ci+1,ck,cj+1)
3778 END IF
3779 nfld(ni,nk,nj) = sum / REAL ( icount )
3780 END IF
3781 END DO
3782 END DO
3783 END DO
3785 ELSE
3786 CALL wrf_error_fatal ( "only unstaggered fields right now" )
3787 END IF
3789 END SUBROUTINE p2c_mask
3791 SUBROUTINE none
3792 END SUBROUTINE none
3794 SUBROUTINE smoother ( cfld , &
3795 cids, cide, ckds, ckde, cjds, cjde, &
3796 cims, cime, ckms, ckme, cjms, cjme, &
3797 cits, cite, ckts, ckte, cjts, cjte, &
3798 nids, nide, nkds, nkde, njds, njde, &
3799 nims, nime, nkms, nkme, njms, njme, &
3800 nits, nite, nkts, nkte, njts, njte, &
3801 xstag, ystag, & ! staggering of field
3802 ipos, jpos, & ! Position of lower left of nest in
3803 nri, nrj &
3804 )
3806 USE module_configure
3807 IMPLICIT NONE
3809 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3810 cims, cime, ckms, ckme, cjms, cjme, &
3811 cits, cite, ckts, ckte, cjts, cjte, &
3812 nids, nide, nkds, nkde, njds, njde, &
3813 nims, nime, nkms, nkme, njms, njme, &
3814 nits, nite, nkts, nkte, njts, njte, &
3815 nri, nrj, &
3816 ipos, jpos
3817 LOGICAL, INTENT(IN) :: xstag, ystag
3818 INTEGER :: smooth_option, feedback , spec_zone
3820 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
3822 ! If there is no feedback, there can be no smoothing.
3824 CALL nl_get_feedback ( 1, feedback )
3825 IF ( feedback == 0 ) RETURN
3826 CALL nl_get_spec_zone ( 1, spec_zone )
3828 ! These are the 2d smoothers used on the fedback data. These filters
3829 ! are run on the coarse grid data (after the nested info has been
3830 ! fedback). Only the area of the nest in the coarse grid is filtered.
3832 CALL nl_get_smooth_option ( 1, smooth_option )
3834 IF ( smooth_option == 0 ) THEN
3835 ! no op
3836 ELSE IF ( smooth_option == 1 ) THEN
3837 CALL sm121 ( cfld , &
3838 cids, cide, ckds, ckde, cjds, cjde, &
3839 cims, cime, ckms, ckme, cjms, cjme, &
3840 cits, cite, ckts, ckte, cjts, cjte, &
3841 xstag, ystag, & ! staggering of field
3842 nids, nide, nkds, nkde, njds, njde, &
3843 nims, nime, nkms, nkme, njms, njme, &
3844 nits, nite, nkts, nkte, njts, njte, &
3845 nri, nrj, &
3846 ipos, jpos & ! Position of lower left of nest in
3847 )
3848 ELSE IF ( smooth_option == 2 ) THEN
3849 CALL smdsm ( cfld , &
3850 cids, cide, ckds, ckde, cjds, cjde, &
3851 cims, cime, ckms, ckme, cjms, cjme, &
3852 cits, cite, ckts, ckte, cjts, cjte, &
3853 xstag, ystag, & ! staggering of field
3854 nids, nide, nkds, nkde, njds, njde, &
3855 nims, nime, nkms, nkme, njms, njme, &
3856 nits, nite, nkts, nkte, njts, njte, &
3857 nri, nrj, &
3858 ipos, jpos & ! Position of lower left of nest in
3859 )
3860 END IF
3862 END SUBROUTINE smoother
3864 SUBROUTINE sm121 ( cfld , &
3865 cids, cide, ckds, ckde, cjds, cjde, &
3866 cims, cime, ckms, ckme, cjms, cjme, &
3867 cits, cite, ckts, ckte, cjts, cjte, &
3868 xstag, ystag, & ! staggering of field
3869 nids, nide, nkds, nkde, njds, njde, &
3870 nims, nime, nkms, nkme, njms, njme, &
3871 nits, nite, nkts, nkte, njts, njte, &
3872 nri, nrj, &
3873 ipos, jpos & ! Position of lower left of nest in
3874 )
3876 USE module_configure
3877 IMPLICIT NONE
3879 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3880 cims, cime, ckms, ckme, cjms, cjme, &
3881 cits, cite, ckts, ckte, cjts, cjte, &
3882 nids, nide, nkds, nkde, njds, njde, &
3883 nims, nime, nkms, nkme, njms, njme, &
3884 nits, nite, nkts, nkte, njts, njte, &
3885 nri, nrj, &
3886 ipos, jpos
3887 LOGICAL, INTENT(IN) :: xstag, ystag
3889 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
3890 REAL, DIMENSION ( cims:cime, cjms:cjme ) :: cfldnew
3892 INTEGER :: i , j , k , loop
3893 INTEGER :: istag,jstag
3895 INTEGER, PARAMETER :: smooth_passes = 1 ! More passes requires a larger stencil (currently 48 pt)
3897 istag = 1 ; jstag = 1
3898 IF ( xstag ) istag = 0
3899 IF ( ystag ) jstag = 0
3901 ! Simple 1-2-1 smoother.
3903 smoothing_passes : DO loop = 1 , smooth_passes
3905 DO k = ckts , ckte
3907 ! Initialize dummy cfldnew
3909 DO i = MAX(ipos,cits-3) , MIN(ipos+(nide-nids)/nri,cite+3)
3910 DO j = MAX(jpos,cjts-3) , MIN(jpos+(njde-njds)/nrj,cjte+3)
3911 cfldnew(i,j) = cfld(i,k,j)
3912 END DO
3913 END DO
3915 ! 1-2-1 smoothing in the j direction first,
3917 DO i = MAX(ipos+2,cits-2) , MIN(ipos+(nide-nids)/nri-2-istag,cite+2)
3918 DO j = MAX(jpos+2,cjts-2) , MIN(jpos+(njde-njds)/nrj-2-jstag,cjte+2)
3919 cfldnew(i,j) = 0.25 * ( cfld(i,k,j+1) + 2.*cfld(i,k,j) + cfld(i,k,j-1) )
3920 END DO
3921 END DO
3923 ! then 1-2-1 smoothing in the i direction last
3925 DO j = MAX(jpos+2,cjts-2) , MIN(jpos+(njde-njds)/nrj-2-jstag,cjte+2)
3926 DO i = MAX(ipos+2,cits-2) , MIN(ipos+(nide-nids)/nri-2-istag,cite+2)
3927 cfld(i,k,j) = 0.25 * ( cfldnew(i+1,j) + 2.*cfldnew(i,j) + cfldnew(i-1,j) )
3928 END DO
3929 END DO
3931 END DO
3933 END DO smoothing_passes
3935 END SUBROUTINE sm121
3937 SUBROUTINE smdsm ( cfld , &
3938 cids, cide, ckds, ckde, cjds, cjde, &
3939 cims, cime, ckms, ckme, cjms, cjme, &
3940 cits, cite, ckts, ckte, cjts, cjte, &
3941 xstag, ystag, & ! staggering of field
3942 nids, nide, nkds, nkde, njds, njde, &
3943 nims, nime, nkms, nkme, njms, njme, &
3944 nits, nite, nkts, nkte, njts, njte, &
3945 nri, nrj, &
3946 ipos, jpos & ! Position of lower left of nest in
3947 )
3949 USE module_configure
3950 IMPLICIT NONE
3952 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3953 cims, cime, ckms, ckme, cjms, cjme, &
3954 cits, cite, ckts, ckte, cjts, cjte, &
3955 nids, nide, nkds, nkde, njds, njde, &
3956 nims, nime, nkms, nkme, njms, njme, &
3957 nits, nite, nkts, nkte, njts, njte, &
3958 nri, nrj, &
3959 ipos, jpos
3960 LOGICAL, INTENT(IN) :: xstag, ystag
3962 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
3963 REAL, DIMENSION ( cims:cime, cjms:cjme ) :: cfldnew
3965 REAL , DIMENSION ( 2 ) :: xnu
3966 INTEGER :: i , j , k , loop , n
3967 INTEGER :: istag,jstag
3969 INTEGER, PARAMETER :: smooth_passes = 1 ! More passes requires a larger stencil (currently 48 pt)
3971 xnu = (/ 0.50 , -0.52 /)
3973 istag = 1 ; jstag = 1
3974 IF ( xstag ) istag = 0
3975 IF ( ystag ) jstag = 0
3977 ! The odd number passes of this are the "smoother", the even
3978 ! number passes are the "de-smoother" (note the different signs on xnu).
3980 smoothing_passes : DO loop = 1 , smooth_passes * 2
3982 n = 2 - MOD ( loop , 2 )
3984 DO k = ckts , ckte
3986 DO i = MAX(ipos+2,cits-2) , MIN(ipos+(nide-nids)/nri-2-istag,cite+2)
3987 DO j = MAX(jpos+2,cjts-2) , MIN(jpos+(njde-njds)/nrj-2-jstag,cjte+2)
3988 cfldnew(i,j) = cfld(i,k,j) + xnu(n) * ((cfld(i,k,j+1) + cfld(i,k,j-1)) * 0.5-cfld(i,k,j))
3989 END DO
3990 END DO
3992 DO i = MAX(ipos+2,cits-2) , MIN(ipos+(nide-nids)/nri-2-istag,cite+2)
3993 DO j = MAX(jpos+2,cjts-2) , MIN(jpos+(njde-njds)/nrj-2-jstag,cjte+2)
3994 cfld(i,k,j) = cfldnew(i,j)
3995 END DO
3996 END DO
3998 DO j = MAX(jpos+2,cjts-2) , MIN(jpos+(njde-njds)/nrj-2-jstag,cjte+2)
3999 DO i = MAX(ipos+2,cits-2) , MIN(ipos+(nide-nids)/nri-2-istag,cite+2)
4000 cfldnew(i,j) = cfld(i,k,j) + xnu(n) * ((cfld(i+1,k,j) + cfld(i-1,k,j)) * 0.5-cfld(i,k,j))
4001 END DO
4002 END DO
4004 DO j = MAX(jpos+2,cjts-2) , MIN(jpos+(njde-njds)/nrj-2-jstag,cjte+2)
4005 DO i = MAX(ipos+2,cits-2) , MIN(ipos+(nide-nids)/nri-2-istag,cite+2)
4006 cfld(i,k,j) = cfldnew(i,j)
4007 END DO
4008 END DO
4010 END DO
4012 END DO smoothing_passes
4014 END SUBROUTINE smdsm
4016 !==================================
4017 ! this is used to modify a field over the nest so we can see where the nest is
4019 SUBROUTINE mark_domain ( cfld, & ! CD field
4020 cids, cide, ckds, ckde, cjds, cjde, &
4021 cims, cime, ckms, ckme, cjms, cjme, &
4022 cits, cite, ckts, ckte, cjts, cjte, &
4023 nfld, & ! ND field
4024 nids, nide, nkds, nkde, njds, njde, &
4025 nims, nime, nkms, nkme, njms, njme, &
4026 nits, nite, nkts, nkte, njts, njte, &
4027 shw, & ! stencil half width for interp
4028 imask, & ! interpolation mask
4029 xstag, ystag, & ! staggering of field
4030 ipos, jpos, & ! Position of lower left of nest in CD
4031 nri, nrj ) ! nest ratios
4032 USE module_configure
4033 USE module_wrf_error
4034 IMPLICIT NONE
4037 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
4038 cims, cime, ckms, ckme, cjms, cjme, &
4039 cits, cite, ckts, ckte, cjts, cjte, &
4040 nids, nide, nkds, nkde, njds, njde, &
4041 nims, nime, nkms, nkme, njms, njme, &
4042 nits, nite, nkts, nkte, njts, njte, &
4043 shw, &
4044 ipos, jpos, &
4045 nri, nrj
4046 LOGICAL, INTENT(IN) :: xstag, ystag
4048 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ), INTENT(OUT) :: cfld
4049 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ), INTENT(IN) :: nfld
4050 INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: imask
4052 ! Local
4054 INTEGER ci, cj, ck, ni, nj, nk, ip, jp, ioff, joff, ioffa, joffa
4055 INTEGER :: icmin,icmax,jcmin,jcmax
4056 INTEGER :: istag,jstag, ipoints,jpoints,ijpoints
4058 istag = 1 ; jstag = 1
4059 IF ( xstag ) istag = 0
4060 IF ( ystag ) jstag = 0
4062 DO cj = MAX(jpos+1,cjts),MIN(jpos+(njde-njds)/nrj-jstag-1,cjte)
4063 nj = (cj-jpos)*nrj + jstag + 1
4064 DO ck = ckts, ckte
4065 nk = ck
4066 DO ci = MAX(ipos+1,cits),MIN(ipos+(nide-nids)/nri-istag-1,cite)
4067 ni = (ci-ipos)*nri + istag + 1
4068 cfld( ci, ck, cj ) = 9021000. !magic number: Beverly Hills * 100.
4069 ENDDO
4070 ENDDO
4071 ENDDO
4073 END SUBROUTINE mark_domain
4075 SUBROUTINE interp_mask_field ( enable, & ! says whether to allow interpolation or just the bcasts
4076 cfld, & ! CD field
4077 cids, cide, ckds, ckde, cjds, cjde, &
4078 cims, cime, ckms, ckme, cjms, cjme, &
4079 cits, cite, ckts, ckte, cjts, cjte, &
4080 nfld, & ! ND field
4081 nids, nide, nkds, nkde, njds, njde, &
4082 nims, nime, nkms, nkme, njms, njme, &
4083 nits, nite, nkts, nkte, njts, njte, &
4084 shw, & ! stencil half width
4085 imask, & ! interpolation mask
4086 xstag, ystag, & ! staggering of field
4087 ipos, jpos, & ! Position of lower left of nest in CD
4088 nri, nrj, & ! nest ratios
4089 clu, nlu, cflag, nflag )
4091 USE module_configure
4092 USE module_wrf_error
4093 USE module_dm , only : wrf_dm_sum_reals, wrf_dm_sum_integers
4095 IMPLICIT NONE
4098 LOGICAL, INTENT(IN) :: enable
4099 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
4100 cims, cime, ckms, ckme, cjms, cjme, &
4101 cits, cite, ckts, ckte, cjts, cjte, &
4102 nids, nide, nkds, nkde, njds, njde, &
4103 nims, nime, nkms, nkme, njms, njme, &
4104 nits, nite, nkts, nkte, njts, njte, &
4105 shw, & ! stencil half width
4106 ipos, jpos, &
4107 nri, nrj, cflag, nflag
4108 LOGICAL, INTENT(IN) :: xstag, ystag
4110 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
4111 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
4112 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
4114 REAL, DIMENSION ( cims:cime, cjms:cjme ) :: clu
4115 REAL, DIMENSION ( nims:nime, njms:njme ) :: nlu
4117 ! Local
4119 INTEGER ci, cj, ck, ni, nj, nk, ip, jp
4120 INTEGER :: icount, ii , jj , ist , ien , jst , jen , iswater, ierr
4121 REAL :: avg, sum, dx , dy
4122 INTEGER :: icount_water(nkts:nkte), icount_land(nkts:nkte), idummy(nkts:nkte)
4123 REAL :: avg_water(nkts:nkte), avg_land(nkts:nkte), sum_water(nkts:nkte), sum_land(nkts:nkte), dummy(nkts:nkte)
4124 CHARACTER (len=256) :: message
4125 CHARACTER (len=256) :: a_mess
4127 ! Find out what the water value is.
4129 !CALL nl_get_iswater(1,iswater)
4130 iswater = nflag
4132 ! Right now, only mass point locations permitted.
4134 IF ( ( .NOT. xstag) .AND. ( .NOT. ystag ) ) THEN
4136 ! Loop over each i,k,j in the nested domain.
4138 IF ( enable ) THEN
4139 DO nj = njts, njte
4140 ! first coarse position equal to or below nest point
4141 IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
4142 cj = ( nj + (nrj/2)-1 ) / nrj + jpos -1
4143 ELSE
4144 cj = ( nj + (nrj-1)/2 ) / nrj + jpos -1
4145 END IF
4146 DO nk = nkts, nkte
4147 ck = nk
4148 DO ni = nits, nite
4149 IF ( imask(ni, nj) .NE. 1 ) cycle
4150 ! first coarse position equal to or to the left of nest point
4151 IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
4152 ci = ( ni + (nri/2)-1 ) / nri + ipos -1
4153 ELSE
4154 ci = ( ni + (nri-1)/2 ) / nri + ipos -1
4155 END IF
4157 !
4158 ! (ci,cj+1) (ci+1,cj+1)
4159 ! - -------------
4160 ! 1-dy | | |
4161 ! | | |
4162 ! - | * |
4163 ! dy | | (ni,nj) |
4164 ! | | |
4165 ! - -------------
4166 ! (ci,cj) (ci+1,cj)
4167 !
4168 ! |--|--------|
4169 ! dx 1-dx
4171 ! For odd ratios, at (nri+1)/2, we are on the coarse grid point, so dx = 0
4173 IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
4174 dx = ( REAL ( MOD ( ni+(nri-1)/2 , nri ) ) + 0.5 ) / REAL ( nri )
4175 ELSE
4176 dx = REAL ( MOD ( ni+(nri-1)/2 , nri ) ) / REAL ( nri )
4177 END IF
4178 IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
4179 dy = ( REAL ( MOD ( nj+(nrj-1)/2 , nrj ) ) + 0.5 ) / REAL ( nrj )
4180 ELSE
4181 dy = REAL ( MOD ( nj+(nrj-1)/2 , nrj ) ) / REAL ( nrj )
4182 END IF
4184 ! Nested cell is a water cell.
4185 IF ( ( NINT(nlu(ni, nj)) .EQ. iswater ) ) THEN
4187 ! If the surrounding coarse values are all WATER points,
4188 ! i.e. open water, this is a simple 4-pt interpolation.
4189 ! If the surrounding coarse values are all LAND points,
4190 ! i.e. this is a 1-cell lake, we have no better way to
4191 ! come up with the value than to do a simple 4-pt interpolation.
4193 IF ( ALL( clu(ci:ci+1,cj:cj+1) == iswater ) .OR. &
4194 ALL( clu(ci:ci+1,cj:cj+1) /= iswater ) ) THEN
4196 nfld(ni,nk,nj) = ( 1. - dx ) * ( ( 1. - dy ) * cfld(ci ,ck,cj ) + &
4197 dy * cfld(ci ,ck,cj+1) ) + &
4198 dx * ( ( 1. - dy ) * cfld(ci+1,ck,cj ) + &
4199 dy * cfld(ci+1,ck,cj+1) )
4201 ! If there are some land points and some water points, take an average.
4202 ELSE
4203 icount = 0
4204 sum = 0
4205 IF ( NINT(clu(ci ,cj )) .EQ. iswater ) THEN
4206 icount = icount + 1
4207 sum = sum + cfld(ci ,ck,cj )
4208 END IF
4209 IF ( NINT(clu(ci+1,cj )) .EQ. iswater ) THEN
4210 icount = icount + 1
4211 sum = sum + cfld(ci+1,ck,cj )
4212 END IF
4213 IF ( NINT(clu(ci ,cj+1)) .EQ. iswater ) THEN
4214 icount = icount + 1
4215 sum = sum + cfld(ci ,ck,cj+1)
4216 END IF
4217 IF ( NINT(clu(ci+1,cj+1)) .EQ. iswater ) THEN
4218 icount = icount + 1
4219 sum = sum + cfld(ci+1,ck,cj+1)
4220 END IF
4221 nfld(ni,nk,nj) = sum / REAL ( icount )
4222 END IF
4224 ! Nested cell is a land cell.
4225 ELSE IF ( ( NINT(nlu(ni, nj)) .NE. iswater ) ) THEN
4227 ! If the surrounding coarse values are all LAND points,
4228 ! this is a simple 4-pt interpolation.
4229 ! If the surrounding coarse values are all WATER points,
4230 ! i.e. this is a 1-cell island, we have no better way to
4231 ! come up with the value than to do a simple 4-pt interpolation.
4233 IF ( ALL( clu(ci:ci+1,cj:cj+1) == iswater ) .OR. &
4234 ALL( clu(ci:ci+1,cj:cj+1) /= iswater ) ) THEN
4236 nfld(ni,nk,nj) = ( 1. - dx ) * ( ( 1. - dy ) * cfld(ci ,ck,cj ) + &
4237 dy * cfld(ci ,ck,cj+1) ) + &
4238 dx * ( ( 1. - dy ) * cfld(ci+1,ck,cj ) + &
4239 dy * cfld(ci+1,ck,cj+1) )
4241 ! If there are some water points and some land points, take an average.
4242 ELSE
4243 icount = 0
4244 sum = 0
4245 IF ( NINT(clu(ci ,cj )) .NE. iswater ) THEN
4246 icount = icount + 1
4247 sum = sum + cfld(ci ,ck,cj )
4248 END IF
4249 IF ( NINT(clu(ci+1,cj )) .NE. iswater ) THEN
4250 icount = icount + 1
4251 sum = sum + cfld(ci+1,ck,cj )
4252 END IF
4253 IF ( NINT(clu(ci ,cj+1)) .NE. iswater ) THEN
4254 icount = icount + 1
4255 sum = sum + cfld(ci ,ck,cj+1)
4256 END IF
4257 IF ( NINT(clu(ci+1,cj+1)) .NE. iswater ) THEN
4258 icount = icount + 1
4259 sum = sum + cfld(ci+1,ck,cj+1)
4260 END IF
4261 nfld(ni,nk,nj) = sum / REAL ( icount )
4263 END IF
4264 END IF
4266 END DO
4267 END DO
4268 END DO
4270 END IF
4271 ELSE
4272 CALL wrf_error_fatal ( "only unstaggered fields right now" )
4273 END IF
4275 END SUBROUTINE interp_mask_field
4278 SUBROUTINE interp_mask_soil ( enable, & ! says whether to allow interpolation or just the bcasts
4279 cfld, & ! CD field
4280 cids, cide, ckds, ckde, cjds, cjde, &
4281 cims, cime, ckms, ckme, cjms, cjme, &
4282 cits, cite, ckts, ckte, cjts, cjte, &
4283 nfld, & ! ND field
4284 nids, nide, nkds, nkde, njds, njde, &
4285 nims, nime, nkms, nkme, njms, njme, &
4286 nits, nite, nkts, nkte, njts, njte, &
4287 shw, & ! stencil half width
4288 imask, & ! interpolation mask
4289 xstag, ystag, & ! staggering of field
4290 ipos, jpos, & ! Position of lower left of nest in CD
4291 nri, nrj, & ! nest ratios
4292 clu, nlu )
4294 USE module_configure
4295 USE module_wrf_error
4296 USE module_dm , only : wrf_dm_sum_real, wrf_dm_sum_integer
4298 IMPLICIT NONE
4301 LOGICAL, INTENT(IN) :: enable
4302 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
4303 cims, cime, ckms, ckme, cjms, cjme, &
4304 cits, cite, ckts, ckte, cjts, cjte, &
4305 nids, nide, nkds, nkde, njds, njde, &
4306 nims, nime, nkms, nkme, njms, njme, &
4307 nits, nite, nkts, nkte, njts, njte, &
4308 shw, & ! stencil half width
4309 ipos, jpos, &
4310 nri, nrj
4311 LOGICAL, INTENT(IN) :: xstag, ystag
4313 INTEGER, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
4314 INTEGER, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
4315 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
4317 REAL, DIMENSION ( cims:cime, cjms:cjme ) :: clu
4318 REAL, DIMENSION ( nims:nime, njms:njme ) :: nlu
4320 ! Local
4322 INTEGER ci, cj, ck, ni, nj, nk, ip, jp
4323 INTEGER :: icount, ii , jj , ist , ien , jst , jen , iswater, num_soil_cat, ierr
4324 REAL :: avg, sum, dx , dy
4325 INTEGER , ALLOCATABLE :: icount_water(:,: ), icount_land(:,:)
4326 INTEGER , PARAMETER :: max_search = 5
4327 CHARACTER*120 message
4328 INTEGER, PARAMETER :: isoilwater = 14
4330 CALL nl_get_iswater(1,iswater)
4331 CALL nl_get_num_soil_cat(1,num_soil_cat)
4333 allocate (icount_water(nkms:nkme,1:num_soil_cat))
4334 allocate ( icount_land(nkms:nkme,1:num_soil_cat))
4336 ! Right now, only mass point locations permitted.
4338 IF ( ( .NOT. xstag) .AND. ( .NOT. ystag ) ) THEN
4340 ! Loop over each i,k,j in the nested domain.
4342 IF ( enable ) THEN
4344 DO nj = njts, njte
4345 cj = jpos + (nj-1) / nrj ! j coord of CD point
4346 DO nk = nkts, nkte
4347 ck = nk
4348 DO ni = nits, nite
4349 ci = ipos + (ni-1) / nri ! j coord of CD point
4351 IF ( imask(ni, nj) .NE. 1 ) cycle
4353 IF ( ( NINT(nlu(ni, nj)) .EQ. iswater ) ) then
4355 IF ( ( NINT(clu(ci ,cj )) .EQ. iswater ) ) then
4356 nfld(ni,nk,nj) = cfld(ci,ck,cj)
4357 ELSE
4358 nfld(ni,nk,nj) = -1
4359 ENDIF
4361 ELSE IF ( ( NINT(nlu(ni, nj)) .NE. iswater ) ) THEN
4363 IF ( ( NINT(clu(ci ,cj )) .NE. iswater ) ) THEN
4364 nfld(ni,nk,nj) = cfld(ci,ck,cj)
4365 ELSE
4366 nfld(ni,nk,nj) = -1
4367 ENDIF
4369 END IF
4370 END DO
4371 END DO
4372 END DO
4374 DO nj = njts, njte
4375 DO nk = nkts, nkte
4376 DO ni = nits, nite
4377 IF ( imask(ni, nj) .NE. 1 ) cycle
4378 IF ( nfld(ni,nk,nj) .EQ. -1 ) THEN
4379 IF ( NINT(nlu(ni,nj)) .EQ. iswater ) THEN
4380 nfld(ni,nk,nj) = isoilwater
4381 END IF
4382 END IF
4383 END DO
4384 END DO
4385 END DO
4386 #if 0
4387 IF ( ANY(nfld .EQ. -1) ) THEN
4389 ! Get an average of the whole domain for problem locations.
4391 sum_water = 0
4392 icount_water = 0
4393 sum_land = 0
4394 icount_land = 0
4395 avg_water = 0
4396 avg_land = 0
4398 DO nj = njts, njte
4399 cj = jpos + (nj-1) / nrj ! j coord of CD point
4400 DO nk = nkts, nkte
4401 DO ni = nits, nite
4402 ci = ipos + (ni-1) / nri ! j coord of CD point
4403 IF ( imask(ni, nj) .NE. 1 ) cycle
4404 IF ( nfld(ni,nk,nj) .EQ. -1 ) THEN
4405 ist = MAX (ci-max_search,cits)
4406 ien = MIN (ci+max_search,cite,cide-1)
4407 jst = MAX (cj-max_search,cjts)
4408 jen = MIN (cj+max_search,cjte,cjde-1)
4409 DO jj = jst,jen
4410 DO ii = ist,ien
4411 IF ( NINT(clu(ii,jj)) .NE. iswater ) THEN
4412 icount_land(nk,cfld(ii,nk,jj)) = icount_land(nk,cfld(ii,nk,jj)) +1
4413 END IF
4414 END DO
4415 END DO
4416 IF ( maxval(icount_land(nk,:)) .GT. 0 .and. maxloc(icount_land(nk,:)) .ne. isoilwater ) then
4417 nfld(ni,nk,nj) = maxloc(icount_land(nk,:))
4418 END IF
4419 END IF
4420 END DO
4421 END DO
4422 END DO
4424 END IF ! nfld = -1
4427 IF ( ANY(nfld .EQ. -1) ) THEN
4428 sum_water = 0
4429 icount_water = 0
4430 sum_land = 0
4431 icount_land = 0
4432 avg_water = 0
4433 avg_land = 0
4435 DO nj = njts, njte
4436 DO nk = nkts, nkte
4437 DO ni = nits, nite
4438 IF ( nlu(ni,nj ) .NE. iswater ) THEN
4439 icount_land(nk,nfld(ni,nk,nj)) = icount_land(nk,nfld(ni,nk,nj)) +1
4440 END IF
4441 ENDDO
4442 ENDDO
4443 ENDDO
4445 DO nj = njts, njte
4446 DO nk = nkts, nkte
4447 DO ni = nits, nite
4448 IF ( imask(ni, nj) .NE. 1 ) cycle
4449 IF ( nfld(ni,nk,nj) .EQ. -1 .and. maxloc(icount_land(nk,:)) .ne. isoilwater) THEN
4450 nfld(ni,nk,nj) = MAXLOC(icount_land(nk,:))
4451 END IF
4452 ENDDO
4453 ENDDO
4454 ENDDO
4455 END IF ! nfld = -1
4456 #endif
4458 IF ( ANY(nfld .EQ. -1) ) THEN
4459 DO nj = njts, njte
4460 DO nk = nkts, nkte
4461 DO ni = nits, nite
4462 IF ( imask(ni, nj) .NE. 1 ) cycle
4463 IF ( nfld(ni,nk,nj) .EQ. -1 ) THEN
4464 nfld(ni,nk,nj) = 8
4465 END IF
4466 ENDDO
4467 ENDDO
4468 ENDDO
4469 END IF ! nfld = -1
4471 END IF ! enable
4472 ELSE
4473 CALL wrf_error_fatal ( "only unstaggered fields right now" )
4474 END IF
4477 deallocate (icount_water)
4478 deallocate (icount_land)
4480 END SUBROUTINE interp_mask_soil
4482 !=========================================================================
4484 ! Lagrange interpolating polynomials, set up as a quadratic, with an average of
4485 ! the overlap. Specifically for longitude near the date line.
4487 SUBROUTINE interp_fcn_lagr_ll ( cfld_inp, & ! CD field
4488 cids, cide, ckds, ckde, cjds, cjde, &
4489 cims, cime, ckms, ckme, cjms, cjme, &
4490 cits, cite, ckts, ckte, cjts, cjte, &
4491 nfld, & ! ND field
4492 nids, nide, nkds, nkde, njds, njde, &
4493 nims, nime, nkms, nkme, njms, njme, &
4494 nits, nite, nkts, nkte, njts, njte, &
4495 shw, & ! stencil half width for interp
4496 imask, & ! interpolation mask
4497 xstag, ystag, & ! staggering of field
4498 ipos, jpos, & ! Position of lower left of nest in CD
4499 nri, nrj, & ! Nest ratio, i- and j-directions
4500 clat_in, nlat_in, & ! CG, FG latitude
4501 cinput_from_file, ninput_from_file ) ! CG, FG T/F input from file
4503 IMPLICIT NONE
4505 INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
4506 cims, cime, ckms, ckme, cjms, cjme, &
4507 cits, cite, ckts, ckte, cjts, cjte, &
4508 nids, nide, nkds, nkde, njds, njde, &
4509 nims, nime, nkms, nkme, njms, njme, &
4510 nits, nite, nkts, nkte, njts, njte, &
4511 shw, &
4512 ipos, jpos, &
4513 nri, nrj
4514 LOGICAL, INTENT(IN) :: xstag, ystag
4516 REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld_inp, cfld
4517 REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
4518 REAL, DIMENSION ( cims:cime, cjms:cjme ) :: clat_in
4519 REAL, DIMENSION ( nims:nime, njms:njme ) :: nlat_in
4520 INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
4521 LOGICAL :: cinput_from_file, ninput_from_file
4523 ! Local
4525 INTEGER ci, cj, ck, ni, nj, nk, istag, jstag, i, j, k
4526 REAL :: nx, x0, x1, x2, x3, x
4527 REAL :: ny, y0, y1, y2, y3
4528 REAL :: cxm1, cxp0, cxp1, cxp2, nfld_m1, nfld_p0, nfld_p1, nfld_p2
4529 REAL :: cym1, cyp0, cyp1, cyp2
4530 INTEGER :: ioff, joff
4531 LOGICAL :: probably_by_dateline
4532 REAL :: max_lon, min_lon
4533 LOGICAL :: probably_by_pole
4534 REAL :: max_lat, min_lat
4536 ! Fortran functions.
4538 REAL, EXTERNAL :: lagrange_quad_avg
4539 REAL, EXTERNAL :: nest_loc_of_cg
4540 INTEGER, EXTERNAL :: compute_CGLL
4542 ! This stag stuff is to keep us away from the outer most row
4543 ! and column for the unstaggered directions. We are going to
4544 ! consider "U" an xstag variable and "V" a ystag variable. The
4545 ! vertical staggering is handled in the actual arguments. The
4546 ! ckte and nkte are the ending vertical dimensions for computations
4547 ! for this particular variable.
4549 ! The ioff and joff are offsets due to the staggering. It is a lot
4550 ! simpler with ioff and joff if
4551 ! u var => ioff=1
4552 ! v var => joff=1
4553 ! otherwise zero.
4554 ! Note that is OPPOSITE of the istag, jstag vars. The stag variables are
4555 ! used for the domain dimensions, the offset guys are used in the
4556 ! determination of grid points between the CG and FG
4558 IF ( xstag ) THEN
4559 istag = 0
4560 ioff = 1
4561 ELSE
4562 istag = 1
4563 ioff = 0
4564 END IF
4566 IF ( ystag ) THEN
4567 jstag = 0
4568 joff = 1
4569 ELSE
4570 jstag = 1
4571 joff = 0
4572 END IF
4574 ! If this is a projection where the nest is over the pole, and
4575 ! we are using the parent to interpolate the longitudes, then
4576 ! we are going to have longitude troubles. If this is the case,
4577 ! stop the model right away.
4579 probably_by_pole = .FALSE.
4580 max_lat = -90
4581 min_lat = +90
4582 DO nj = njts, MIN(njde-jstag,njte)
4583 DO ni = nits, MIN(nide-istag,nite)
4584 max_lat = MAX ( nlat_in(ni,nj) , max_lat )
4585 min_lat = MIN ( nlat_in(ni,nj) , min_lat )
4586 END DO
4587 END DO
4589 IF ( ( max_lat .GT. 85 ) .OR. ( ABS(min_lat) .GT. 85 ) ) THEN
4590 probably_by_pole = .TRUE.
4591 END IF
4593 IF ( ( probably_by_pole ) .AND. ( .NOT. ninput_from_file ) ) THEN
4594 CALL wrf_error_fatal ( 'Nest over the pole, single input domain, longitudes will be wrong' )
4595 END IF
4597 ! Initialize to NOT being by dateline.
4599 probably_by_dateline = .FALSE.
4600 max_lon = -180
4601 min_lon = +180
4602 DO nj = njts, MIN(njde-jstag,njte)
4603 cj = compute_CGLL ( nj , jpos , nrj , jstag )
4604 DO ni = nits, MIN(nide-istag,nite)
4605 ci = compute_CGLL ( ni , ipos , nri , istag )
4606 max_lon = MAX ( cfld_inp(ci,1,cj) , max_lon )
4607 min_lon = MIN ( cfld_inp(ci,1,cj) , min_lon )
4608 END DO
4609 END DO
4611 IF ( max_lon - min_lon .GT. 300 ) THEN
4612 probably_by_dateline = .TRUE.
4613 END IF
4615 ! Load "continuous" longitude across the date line
4617 DO cj = MIN(cjts-1,cjms), MAX(cjte+1,cjme)
4618 DO ci = MIN(cits-1,cims), MAX(cite+1,cime)
4619 IF ( ( cfld_inp(ci,ckts,cj) .LT. 0 ) .AND. ( probably_by_dateline ) ) THEN
4620 cfld(ci,ckts,cj) = 360 + cfld_inp(ci,ckts,cj)
4621 ELSE
4622 cfld(ci,ckts,cj) = cfld_inp(ci,ckts,cj)
4623 END IF
4624 END DO
4625 END DO
4627 ! Loop over each j-index on this tile for the nested domain.
4629 j_loop : DO nj = njts, MIN(njde-jstag,njte)
4631 ! This is the lower-left j-index of the CG.
4633 ! Example is 3:1 ratio, mass-point staggering. We have listed sixteen CG values
4634 ! as an example: A through P. For a 3:1 ratio, each of these CG cells has
4635 ! nine associated FG points.
4637 ! |=========|=========|=========|=========|
4638 ! | - - - | - - - | - - - | - - - |
4639 ! | | | | |
4640 ! | - M - | - N d | - O - | - P - |
4641 ! | | | | |
4642 ! | - - - | - - - | - - - | - - - |
4643 ! |=========|=========|=========|=========|
4644 ! | - - - | - - - | - - - | - - - |
4645 ! | | | | |
4646 ! | - I - | - J c | - K - | - L - |
4647 ! | | | | |
4648 ! | - - - | - - - | - - - | - - - |
4649 ! |=========|=========|=========|=========|
4650 ! | - 1 2 | 3 4 5 | 6 7 8 | - - - |
4651 ! | | | | |
4652 ! | - E - | - F b | - G - | - H - |
4653 ! | | | | |
4654 ! | - - - | - - - | - - - | - - - |
4655 ! |=========|=========|=========|=========|
4656 ! | - - - | - - - | - - - | - - - |
4657 ! | | | | |
4658 ! | - A - | - B a | - C - | - D - |
4659 ! | | | | |
4660 ! | - - - | - - - | - - - | - - - |
4661 ! |=========|=========|=========|=========|
4663 ! To interpolate to FG point 4, 5, or 6 we will use CG points: A through P. It is
4664 ! sufficient to find the lower left corner of a 4-point interpolation, and then extend
4665 ! each side by one unit.
4667 ! Here are the lower left hand corners of the following FG points:
4668 ! 1 => E
4669 ! 2 => E
4670 ! 3 => E
4671 ! 4 => F
4672 ! 5 => F
4673 ! 6 => F
4674 ! 7 => G
4675 ! 8 => G
4677 cj = compute_CGLL ( nj , jpos , nrj , jstag )
4679 ! Vertical dim of the nest domain.
4681 k_loop : DO nk = nkts, nkte
4683 ! Loop over each i-index on this tile for the nested domain.
4685 i_loop : DO ni = nits, MIN(nide-istag,nite)
4687 ! The coarse grid location that is to the lower left of the FG point.
4689 ci = compute_CGLL ( ni , ipos , nri , istag )
4691 ! To interpolate to point "*" (look in grid cell "F"):
4692 ! 1. Use ABC to get a quadratic valid at "a"
4693 ! Use BCD to get a quadratic valid at "a"
4694 ! Average these to get the final value for "a"
4695 ! 2. Use EFG to get a quadratic valid at "b"
4696 ! Use FGH to get a quadratic valid at "b"
4697 ! Average these to get the final value for "b"
4698 ! 3. Use IJK to get a quadratic valid at "c"
4699 ! Use JKL to get a quadratic valid at "c"
4700 ! Average these to get the final value for "c"
4701 ! 4. Use MNO to get a quadratic valid at "d"
4702 ! Use NOP to get a quadratic valid at "d"
4703 ! Average these to get the final value for "d"
4704 ! 5. Use abc to get a quadratic valid at "*"
4705 ! Use bcd to get a quadratic valid at "*"
4706 ! Average these to get the final value for "*"
4708 ! |=========|=========|=========|=========|
4709 ! | - - - | - - - | - - - | - - - |
4710 ! | | | | |
4711 ! | - M - | - N d | - O - | - P - |
4712 ! | | | | |
4713 ! | - - - | - - - | - - - | - - - |
4714 ! |=========|=========|=========|=========|
4715 ! | - - - | - - - | - - - | - - - |
4716 ! | | | | |
4717 ! | - I - | - J c | - K - | - L - |
4718 ! | | | | |
4719 ! | - - - | - - - | - - - | - - - |
4720 ! |=========|=========|=========|=========|
4721 ! | - - - | - - * | - - - | - - - |
4722 ! | | | | |
4723 ! | - E - | - F b | - G - | - H - |
4724 ! | | | | |
4725 ! | - - - | - - - | - - - | - - - |
4726 ! |=========|=========|=========|=========|
4727 ! | - - - | - - - | - - - | - - - |
4728 ! | | | | |
4729 ! | - A - | - B a | - C - | - D - |
4730 ! | | | | |
4731 ! | - - - | - - - | - - - | - - - |
4732 ! |=========|=========|=========|=========|
4734 ! Overlapping quadratic interpolation.
4736 IF ( imask ( ni, nj ) .EQ. 1 ) THEN
4738 ! I-direction location of "*"
4740 nx = REAL(ni)
4742 ! I-direction location of "A", "E", "I", "M"
4744 cxm1 = nest_loc_of_cg ( ci-1 , ipos , nri , ioff )
4746 ! I-direction location of "B", "F", "J", "N"
4748 cxp0 = nest_loc_of_cg ( ci , ipos , nri , ioff )
4750 ! I-direction location of "C", "G", "K", "O"
4752 cxp1 = nest_loc_of_cg ( ci+1 , ipos , nri , ioff )
4754 ! I-direction location of "D", "H", "L", "P"
4756 cxp2 = nest_loc_of_cg ( ci+2 , ipos , nri , ioff )
4758 ! Value at "a"
4760 nfld_m1 = lagrange_quad_avg ( nx, cxm1, cxp0, cxp1, cxp2, &
4761 cfld(ci-1,nk,cj-1), cfld(ci+0,nk,cj-1), &
4762 cfld(ci+1,nk,cj-1), cfld(ci+2,nk,cj-1) )
4764 ! Value at "b"
4766 nfld_p0 = lagrange_quad_avg ( nx, cxm1, cxp0, cxp1, cxp2, &
4767 cfld(ci-1,nk,cj+0), cfld(ci+0,nk,cj+0), &
4768 cfld(ci+1,nk,cj+0), cfld(ci+2,nk,cj+0) )
4770 ! Value at "c"
4772 nfld_p1 = lagrange_quad_avg ( nx, cxm1, cxp0, cxp1, cxp2, &
4773 cfld(ci-1,nk,cj+1), cfld(ci+0,nk,cj+1), &
4774 cfld(ci+1,nk,cj+1), cfld(ci+2,nk,cj+1) )
4776 ! Value at "d"
4778 nfld_p2 = lagrange_quad_avg ( nx, cxm1, cxp0, cxp1, cxp2, &
4779 cfld(ci-1,nk,cj+2), cfld(ci+0,nk,cj+2), &
4780 cfld(ci+1,nk,cj+2), cfld(ci+2,nk,cj+2) )
4782 ! J-direction location of "*"
4784 ny = REAL(nj)
4786 ! J-direction location of "A", "B", "C", "D"
4788 cym1 = nest_loc_of_cg ( cj-1 , jpos , nrj , joff )
4790 ! J-direction location of "E", "F", "G", "H"
4792 cyp0 = nest_loc_of_cg ( cj , jpos , nrj , joff )
4794 ! J-direction location of "I", "J", "K", "L"
4796 cyp1 = nest_loc_of_cg ( cj+1 , jpos , nrj , joff )
4798 ! J-direction location of "M", "N", "O", "P"
4800 cyp2 = nest_loc_of_cg ( cj+2 , jpos , nrj , joff )
4802 ! Value at "*"
4804 nfld(ni,nk,nj) = lagrange_quad_avg ( ny, cym1, cyp0, cyp1, &
4805 cyp2, nfld_m1, nfld_p0, nfld_p1, nfld_p2 )
4807 END IF
4809 END DO i_loop
4810 END DO k_loop
4811 END DO j_loop
4813 ! Put nested longitude back into the -180 to 180 range.
4815 DO nj = njts, MIN(njde-jstag,njte)
4816 DO ni = nits, MIN(nide-istag,nite)
4817 IF ( nfld(ni,nkts,nj) .GT. 180 ) THEN
4818 nfld(ni,nkts,nj) = -360 + nfld(ni,nkts,nj)
4819 END IF
4820 END DO
4821 END DO
4823 END SUBROUTINE interp_fcn_lagr_ll