Context Navigation

interp_fcn.F

Last change on this file was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 197.1 KB

Line
1	!WRF:MEDIATION_LAYER:INTERPOLATIONFUNCTION
2	!
3
4	#if (DA_CORE != 1)
5	#define MM5_SINT
6	#endif
7	!#define DUMBCOPY
8
9	! Note, NMM-specific routines moved to end. 20080612. JM
10
11	SUBROUTINE interp_fcn ( cfld, & ! CD field
12	cids, cide, ckds, ckde, cjds, cjde, &
13	cims, cime, ckms, ckme, cjms, cjme, &
14	cits, cite, ckts, ckte, cjts, cjte, &
15	nfld, & ! ND field
16	nids, nide, nkds, nkde, njds, njde, &
17	nims, nime, nkms, nkme, njms, njme, &
18	nits, nite, nkts, nkte, njts, njte, &
19	shw, & ! stencil half width for interp
20	imask, & ! interpolation mask
21	xstag, ystag, & ! staggering of field
22	ipos, jpos, & ! Position of lower left of nest in CD
23	nri, nrj ) ! nest ratios
24	USE module_timing
25	USE module_configure
26	IMPLICIT NONE
27
28
29	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
30	cims, cime, ckms, ckme, cjms, cjme, &
31	cits, cite, ckts, ckte, cjts, cjte, &
32	nids, nide, nkds, nkde, njds, njde, &
33	nims, nime, nkms, nkme, njms, njme, &
34	nits, nite, nkts, nkte, njts, njte, &
35	shw, &
36	ipos, jpos, &
37	nri, nrj
38	LOGICAL, INTENT(IN) :: xstag, ystag
39
40	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
41	REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
42	INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
43
44	! Local
45
46	!logical first
47
48	INTEGER ci, cj, ck, ni, nj, nk, ip, jp, ioff, joff, nioff, njoff
49	#ifdef MM5_SINT
50	INTEGER nfx, ior
51	PARAMETER (ior=2)
52	INTEGER nf
53	REAL psca(cims:cime,cjms:cjme,nri*nrj)
54	LOGICAL icmask( cims:cime, cjms:cjme )
55	INTEGER i,j,k
56	#endif
57
58	! Iterate over the ND tile and compute the values
59	! from the CD tile.
60
61	#ifdef MM5_SINT
62
63	ioff = 0 ; joff = 0
64	nioff = 0 ; njoff = 0
65	IF ( xstag ) THEN
66	ioff = (nri-1)/2
67	nioff = nri
68	ENDIF
69	IF ( ystag ) THEN
70	joff = (nrj-1)/2
71	njoff = nrj
72	ENDIF
73
74	nfx = nri * nrj
75	!$OMP PARALLEL DO &
76	!$OMP PRIVATE ( i,j,k,ni,nj,ci,cj,ip,jp,nk,ck,nf,icmask,psca )
77	DO k = ckts, ckte
78	icmask = .FALSE.
79	DO nf = 1,nfx
80	DO j = cjms,cjme
81	nj = (j-jpos) * nrj + ( nrj / 2 + 1 ) ! j point on nest
82	DO i = cims,cime
83	ni = (i-ipos) * nri + ( nri / 2 + 1 ) ! i point on nest
84	if ( ni .ge. nits-nioff-1 .and. ni .le. nite+nioff+1 .and. nj .ge. njts-njoff-1 .and. nj .le. njte+njoff+1 ) then
85	! if ( imask(ni,nj) .eq. 1 .or. imask(ni-nioff,nj-njoff) .eq. 1 ) then
86	! icmask( i, j ) = .TRUE.
87	! endif
88	if ( imask(ni,nj) .eq. 1 ) then
89	icmask( i, j ) = .TRUE.
90	endif
91	if (ni .ge. nits-nioff .and. nj .ge. njts-njoff ) then
92	if ( imask(ni-nioff,nj-njoff) .eq. 1) then
93	icmask( i, j ) = .TRUE.
94	endif
95	endif
96	endif
97	psca(i,j,nf) = cfld(i,k,j)
98	ENDDO
99	ENDDO
100	ENDDO
101
102	! tile dims in this call to sint are 1-over to account for the fact
103	! that the number of cells on the nest local subdomain is not
104	! necessarily a multiple of the nest ratio in a given dim.
105	! this could be a little less ham-handed.
106
107	!call start_timing
108
109	CALL sint( psca, &
110	cims, cime, cjms, cjme, icmask, &
111	cits-1, cite+1, cjts-1, cjte+1, nrj*nri, xstag, ystag )
112
113	!call end_timing( ' sint ' )
114
115	DO nj = njts, njte+joff
116	cj = jpos + (nj-1) / nrj ! j coord of CD point
117	jp = mod ( nj-1 , nrj ) ! coord of ND w/i CD point
118	nk = k
119	ck = nk
120	DO ni = nits, nite+ioff
121	ci = ipos + (ni-1) / nri ! i coord of CD point
122	ip = mod ( ni-1 , nri ) ! coord of ND w/i CD point
123	if ( imask ( ni, nj ) .eq. 1 .or. imask ( ni-ioff, nj-joff ) .eq. 1 ) then
124	nfld( ni-ioff, nk, nj-joff ) = psca( ci , cj, ip+1 + (jp)*nri )
125	endif
126	ENDDO
127	ENDDO
128	ENDDO
129	!$OMP END PARALLEL DO
130	#endif
131
132	#ifdef DUMBCOPY
133	!write(0,'(") cims:cime, ckms:ckme, cjms:cjme ",6i4)')cims,cime, ckms,ckme, cjms,cjme
134	!write(0,'(") nims:nime, nkms:nkme, njms:njme ",6i4)')nims,nime, nkms,nkme, njms,njme
135	!write(0,'(") cits:cite, ckts:ckte, cjts:cjte ",6i4)')cits,cite, ckts,ckte, cjts,cjte
136	!write(0,'(") nits:nite, nkts:nkte, njts:njte ",6i4)')nits,nite, nkts,nkte, njts,njte
137
138	DO nj = njts, njte
139	cj = jpos + (nj-1) / nrj ! j coord of CD point
140	jp = mod ( nj , nrj ) ! coord of ND w/i CD point
141	DO nk = nkts, nkte
142	ck = nk
143	DO ni = nits, nite
144	ci = ipos + (ni-1) / nri ! j coord of CD point
145	ip = mod ( ni , nri ) ! coord of ND w/i CD point
146	! This is a trivial implementation of the interp_fcn; just copies
147	! the values from the CD into the ND
148	if ( imask ( ni, nj ) .eq. 1 ) then
149	nfld( ni, nk, nj ) = cfld( ci , ck , cj )
150	endif
151	ENDDO
152	ENDDO
153	ENDDO
154	#endif
155
156	RETURN
157
158	END SUBROUTINE interp_fcn
159
160	!==================================
161	! this is the default function used in feedback.
162
163	SUBROUTINE copy_fcn ( cfld, & ! CD field
164	cids, cide, ckds, ckde, cjds, cjde, &
165	cims, cime, ckms, ckme, cjms, cjme, &
166	cits, cite, ckts, ckte, cjts, cjte, &
167	nfld, & ! ND field
168	nids, nide, nkds, nkde, njds, njde, &
169	nims, nime, nkms, nkme, njms, njme, &
170	nits, nite, nkts, nkte, njts, njte, &
171	shw, & ! stencil half width for interp
172	imask, & ! interpolation mask
173	xstag, ystag, & ! staggering of field
174	ipos, jpos, & ! Position of lower left of nest in CD
175	nri, nrj ) ! nest ratios
176	USE module_configure
177	IMPLICIT NONE
178
179
180	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
181	cims, cime, ckms, ckme, cjms, cjme, &
182	cits, cite, ckts, ckte, cjts, cjte, &
183	nids, nide, nkds, nkde, njds, njde, &
184	nims, nime, nkms, nkme, njms, njme, &
185	nits, nite, nkts, nkte, njts, njte, &
186	shw, &
187	ipos, jpos, &
188	nri, nrj
189	LOGICAL, INTENT(IN) :: xstag, ystag
190
191	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ), INTENT(OUT) :: cfld
192	REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ),INTENT(IN) :: nfld
193	INTEGER, DIMENSION ( nims:nime, njms:njme ),INTENT(IN) :: imask
194
195	! Local
196
197	INTEGER ci, cj, ck, ni, nj, nk, ip, jp, ioff, joff, ioffa, joffa
198	INTEGER :: icmin,icmax,jcmin,jcmax
199	INTEGER :: istag,jstag, ipoints,jpoints,ijpoints
200	INTEGER , PARAMETER :: passes = 2
201	INTEGER spec_zone
202
203	! Loop over the coarse grid in the area of the fine mesh. Do not
204	! process the coarse grid values that are along the lateral BC
205	! provided to the fine grid. Since that is in the specified zone
206	! for the fine grid, it should not be used in any feedback to the
207	! coarse grid as it should not have changed.
208
209	! Due to peculiarities of staggering, it is simpler to handle the feedback
210	! for the staggerings based upon whether it is a even ratio (2::1, 4::1, etc.) or
211	! an odd staggering ratio (3::1, 5::1, etc.).
212
213	! Though there are separate grid ratios for the i and j directions, this code
214	! is not general enough to handle aspect ratios .NE. 1 for the fine grid cell.
215
216	! These are local integer increments in the looping. Basically, istag=1 means
217	! that we will assume one less point in the i direction. Note that ci and cj
218	! have a maximum value that is decreased by istag and jstag, respectively.
219
220	! Horizontal momentum feedback is along the face, not within the cell. For a
221	! 3::1 ratio, temperature would use 9 pts for feedback, while u and v use
222	! only 3 points for feedback from the nest to the parent.
223
224	CALL nl_get_spec_zone( 1 , spec_zone )
225	istag = 1 ; jstag = 1
226	IF ( xstag ) istag = 0
227	IF ( ystag ) jstag = 0
228
229	IF( MOD(nrj,2) .NE. 0) THEN ! odd refinement ratio
230
231	IF ( ( .NOT. xstag ) .AND. ( .NOT. ystag ) ) THEN
232	DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
233	nj = (cj-jpos)*nrj + jstag + 1
234	DO ck = ckts, ckte
235	nk = ck
236	DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
237	ni = (ci-ipos)*nri + istag + 1
238	cfld( ci, ck, cj ) = 0.
239	DO ijpoints = 1 , nri * nrj
240	ipoints = MOD((ijpoints-1),nri) + 1 - nri/2 - 1
241	jpoints = (ijpoints-1)/nri + 1 - nrj/2 - 1
242	cfld( ci, ck, cj ) = cfld( ci, ck, cj ) + &
243	1./REAL(nrinrj) nfld( ni+ipoints , nk , nj+jpoints )
244	END DO
245	! cfld( ci, ck, cj ) = 1./9. * &
246	! ( nfld( ni-1, nk , nj-1) + &
247	! nfld( ni , nk , nj-1) + &
248	! nfld( ni+1, nk , nj-1) + &
249	! nfld( ni-1, nk , nj ) + &
250	! nfld( ni , nk , nj ) + &
251	! nfld( ni+1, nk , nj ) + &
252	! nfld( ni-1, nk , nj+1) + &
253	! nfld( ni , nk , nj+1) + &
254	! nfld( ni+1, nk , nj+1) )
255	ENDDO
256	ENDDO
257	ENDDO
258
259	ELSE IF ( ( xstag ) .AND. ( .NOT. ystag ) ) THEN
260	DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
261	nj = (cj-jpos)*nrj + jstag + 1
262	DO ck = ckts, ckte
263	nk = ck
264	DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
265	ni = (ci-ipos)*nri + istag + 1
266	cfld( ci, ck, cj ) = 0.
267	DO ijpoints = (nri+1)/2 , (nri+1)/2 + nri*(nri-1) , nri
268	ipoints = MOD((ijpoints-1),nri) + 1 - nri/2 - 1
269	jpoints = (ijpoints-1)/nri + 1 - nrj/2 - 1
270	cfld( ci, ck, cj ) = cfld( ci, ck, cj ) + &
271	1./REAL(nri ) * nfld( ni+ipoints , nk , nj+jpoints )
272	END DO
273	! cfld( ci, ck, cj ) = 1./3. * &
274	! ( nfld( ni , nk , nj-1) + &
275	! nfld( ni , nk , nj ) + &
276	! nfld( ni , nk , nj+1) )
277	ENDDO
278	ENDDO
279	ENDDO
280
281	ELSE IF ( ( .NOT. xstag ) .AND. ( ystag ) ) THEN
282	DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
283	nj = (cj-jpos)*nrj + jstag + 1
284	DO ck = ckts, ckte
285	nk = ck
286	DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
287	ni = (ci-ipos)*nri + istag + 1
288	cfld( ci, ck, cj ) = 0.
289	DO ijpoints = ( nrjnrj +1 )/2 - nrj/2 , ( nrjnrj +1 )/2 - nrj/2 + nrj-1
290	ipoints = MOD((ijpoints-1),nri) + 1 - nri/2 - 1
291	jpoints = (ijpoints-1)/nri + 1 - nrj/2 - 1
292	cfld( ci, ck, cj ) = cfld( ci, ck, cj ) + &
293	1./REAL( nrj) * nfld( ni+ipoints , nk , nj+jpoints )
294	END DO
295	! cfld( ci, ck, cj ) = 1./3. * &
296	! ( nfld( ni-1, nk , nj ) + &
297	! nfld( ni , nk , nj ) + &
298	! nfld( ni+1, nk , nj ) )
299	ENDDO
300	ENDDO
301	ENDDO
302
303	END IF
304
305	! Even refinement ratio
306
307	ELSE IF ( MOD(nrj,2) .EQ. 0) THEN
308	IF ( ( .NOT. xstag ) .AND. ( .NOT. ystag ) ) THEN
309
310	! This is a simple schematic of the feedback indexing used in the even
311	! ratio nests. For simplicity, a 2::1 ratio is depicted. Only the
312	! mass variable staggering is shown.
313	! Each of
314	! the boxes with a "T" and four small "t" represents a coarse grid (CG)
315	! cell, that is composed of four (2::1 ratio) fine grid (FG) cells.
316
317	! Shown below is the area of the CG that is in the area of the FG. The
318	! first grid point of the depicted CG is the starting location of the nest
319	! in the parent domain (ipos,jpos - i_parent_start and j_parent_start from
320	! the namelist).
321
322	! For each of the CG points, the feedback loop is over each of the FG points
323	! within the CG cell. For a 2::1 ratio, there are four total points (this is
324	! the ijpoints loop). The feedback value to the CG is the arithmetic mean of
325	! all of the FG values within each CG cell.
326
327	! \|-------------\|\|-------------\| \|-------------\|\|-------------\|
328	! \| t t \|\| t t \| \| t t \|\| t t \|
329	! jpos+ \| \|\| \| \| \|\| \|
330	! (njde-njds)- \| T \|\| T \| \| T \|\| T \|
331	! jstag \| \|\| \| \| \|\| \|
332	! \| t t \|\| t t \| \| t t \|\| t t \|
333	! \|-------------\|\|-------------\| \|-------------\|\|-------------\|
334	! \|-------------\|\|-------------\| \|-------------\|\|-------------\|
335	! \| t t \|\| t t \| \| t t \|\| t t \|
336	! \| \|\| \| \| \|\| \|
337	! \| T \|\| T \| \| T \|\| T \|
338	! \| \|\| \| \| \|\| \|
339	! \| t t \|\| t t \| \| t t \|\| t t \|
340	! \|-------------\|\|-------------\| \|-------------\|\|-------------\|
341	!
342	! ...
343	! ...
344	! ...
345	! ...
346	! ...
347
348	! \|-------------\|\|-------------\| \|-------------\|\|-------------\|
349	! jpoints = 1 \| t t \|\| t t \| \| t t \|\| t t \|
350	! \| \|\| \| \| \|\| \|
351	! \| T \|\| T \| \| T \|\| T \|
352	! \| \|\| \| \| \|\| \|
353	! jpoints = 0, \| t t \|\| t t \| \| t t \|\| t t \|
354	! nj=3 \|-------------\|\|-------------\| \|-------------\|\|-------------\|
355	! \|-------------\|\|-------------\| \|-------------\|\|-------------\|
356	! jpoints = 1 \| t t \|\| t t \| \| t t \|\| t t \|
357	! \| \|\| \| \| \|\| \|
358	! jpos \| T \|\| T \| ... \| T \|\| T \|
359	! \| \|\| \| ... \| \|\| \|
360	! jpoints = 0, \| t t \|\| t t \| ... \| t t \|\| t t \|
361	! nj=1 \|-------------\|\|-------------\| \|-------------\|\|-------------\|
362	! ^ ^
363	! \| \|
364	! \| \|
365	! ipos ipos+
366	! ni = 1 3 (nide-nids)/nri
367	! ipoints= 0 1 0 1 -istag
368	!
369
370	! For performance benefits, users can comment out the inner most loop (and cfld=0) and
371	! hardcode the loop feedback. For example, it is set up to run a 2::1 ratio
372	! if uncommented. This lacks generality, but is likely to gain timing benefits
373	! with compilers unable to unroll inner loops that do not have parameterized sizes.
374
375	! The extra +1 ---------/ and the extra -1 ----\ (both for ci and cj)
376	! / \ keeps the feedback out of the
377	! / \ outer row/col, since that CG data
378	! / \ specified the nest boundary originally
379	! / \ This
380	! / \ is just
381	! / \ a sentence to not end a line
382	! / \ with a stupid backslash
383	DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
384	nj = (cj-jpos)*nrj + jstag
385	DO ck = ckts, ckte
386	nk = ck
387	DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
388	ni = (ci-ipos)*nri + istag
389	cfld( ci, ck, cj ) = 0.
390	DO ijpoints = 1 , nri * nrj
391	ipoints = MOD((ijpoints-1),nri)
392	jpoints = (ijpoints-1)/nri
393	cfld( ci, ck, cj ) = cfld( ci, ck, cj ) + &
394	1./REAL(nrinrj) nfld( ni+ipoints , nk , nj+jpoints )
395	END DO
396	! cfld( ci, ck, cj ) = 1./4. * &
397	! ( nfld( ni , nk , nj ) + &
398	! nfld( ni+1, nk , nj ) + &
399	! nfld( ni , nk , nj+1) + &
400	! nfld( ni+1, nk , nj+1) )
401	END DO
402	END DO
403	END DO
404
405	! U
406
407	ELSE IF ( ( xstag ) .AND. ( .NOT. ystag ) ) THEN
408	! \|---------------\|
409	! \| \|
410	! jpoints = 1 u u \|
411	! \| \|
412	! U \|
413	! \| \|
414	! jpoints = 0, u u \|
415	! nj=3 \| \|
416	! \|---------------\|
417	! \|---------------\|
418	! \| \|
419	! jpoints = 1 u u \|
420	! \| \|
421	! jpos U \|
422	! \| \|
423	! jpoints = 0, u u \|
424	! nj=1 \| \|
425	! \|---------------\|
426	!
427	! ^
428	! \|
429	! \|
430	! ipos
431	! ni = 1 3
432	! ipoints= 0 1 0
433	!
434
435	DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
436	nj = (cj-jpos)*nrj + 1
437	DO ck = ckts, ckte
438	nk = ck
439	DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
440	ni = (ci-ipos)*nri + 1
441	cfld( ci, ck, cj ) = 0.
442	DO ijpoints = 1 , nri*nrj , nri
443	ipoints = MOD((ijpoints-1),nri)
444	jpoints = (ijpoints-1)/nri
445	cfld( ci, ck, cj ) = cfld( ci, ck, cj ) + &
446	1./REAL(nri ) * nfld( ni+ipoints , nk , nj+jpoints )
447	END DO
448	! cfld( ci, ck, cj ) = 1./2. * &
449	! ( nfld( ni , nk , nj ) + &
450	! nfld( ni , nk , nj+1) )
451	ENDDO
452	ENDDO
453	ENDDO
454
455	! V
456
457	ELSE IF ( ( .NOT. xstag ) .AND. ( ystag ) ) THEN
458	DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
459	nj = (cj-jpos)*nrj + 1
460	DO ck = ckts, ckte
461	nk = ck
462	DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
463	ni = (ci-ipos)*nri + 1
464	cfld( ci, ck, cj ) = 0.
465	DO ijpoints = 1 , nri
466	ipoints = MOD((ijpoints-1),nri)
467	jpoints = (ijpoints-1)/nri
468	cfld( ci, ck, cj ) = cfld( ci, ck, cj ) + &
469	1./REAL(nri ) * nfld( ni+ipoints , nk , nj+jpoints )
470	END DO
471	! cfld( ci, ck, cj ) = 1./2. * &
472	! ( nfld( ni , nk , nj ) + &
473	! nfld( ni+1, nk , nj ) )
474	ENDDO
475	ENDDO
476	ENDDO
477	END IF
478	END IF
479
480	RETURN
481
482	END SUBROUTINE copy_fcn
483
484	!==================================
485	! this is the 1pt function used in feedback.
486
487	SUBROUTINE copy_fcnm ( cfld, & ! CD field
488	cids, cide, ckds, ckde, cjds, cjde, &
489	cims, cime, ckms, ckme, cjms, cjme, &
490	cits, cite, ckts, ckte, cjts, cjte, &
491	nfld, & ! ND field
492	nids, nide, nkds, nkde, njds, njde, &
493	nims, nime, nkms, nkme, njms, njme, &
494	nits, nite, nkts, nkte, njts, njte, &
495	shw, & ! stencil half width for interp
496	imask, & ! interpolation mask
497	xstag, ystag, & ! staggering of field
498	ipos, jpos, & ! Position of lower left of nest in CD
499	nri, nrj ) ! nest ratios
500	USE module_configure
501	USE module_wrf_error
502	IMPLICIT NONE
503
504
505	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
506	cims, cime, ckms, ckme, cjms, cjme, &
507	cits, cite, ckts, ckte, cjts, cjte, &
508	nids, nide, nkds, nkde, njds, njde, &
509	nims, nime, nkms, nkme, njms, njme, &
510	nits, nite, nkts, nkte, njts, njte, &
511	shw, &
512	ipos, jpos, &
513	nri, nrj
514	LOGICAL, INTENT(IN) :: xstag, ystag
515
516	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ), INTENT(OUT) :: cfld
517	REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ), INTENT(IN) :: nfld
518	INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: imask
519
520	! Local
521
522	INTEGER ci, cj, ck, ni, nj, nk, ip, jp, ioff, joff, ioffa, joffa
523	INTEGER :: icmin,icmax,jcmin,jcmax
524	INTEGER :: istag,jstag, ipoints,jpoints,ijpoints
525	INTEGER , PARAMETER :: passes = 2
526	INTEGER spec_zone
527
528	CALL nl_get_spec_zone( 1, spec_zone )
529	istag = 1 ; jstag = 1
530	IF ( xstag ) istag = 0
531	IF ( ystag ) jstag = 0
532
533	IF( MOD(nrj,2) .NE. 0) THEN ! odd refinement ratio
534
535	DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
536	nj = (cj-jpos)*nrj + jstag + 1
537	DO ck = ckts, ckte
538	nk = ck
539	DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
540	ni = (ci-ipos)*nri + istag + 1
541	cfld( ci, ck, cj ) = nfld( ni , nk , nj )
542	ENDDO
543	ENDDO
544	ENDDO
545
546	ELSE ! even refinement ratio, pick nearest neighbor on SW corner
547	DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
548	nj = (cj-jpos)*nrj + 1
549	DO ck = ckts, ckte
550	nk = ck
551	DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
552	ni = (ci-ipos)*nri + 1
553	ipoints = nri/2 -1
554	jpoints = nrj/2 -1
555	cfld( ci, ck, cj ) = nfld( ni+ipoints , nk , nj+jpoints )
556	END DO
557	END DO
558	END DO
559
560	END IF
561
562	RETURN
563
564	END SUBROUTINE copy_fcnm
565
566	!==================================
567	! this is the 1pt function used in feedback for integers
568
569	SUBROUTINE copy_fcni ( cfld, & ! CD field
570	cids, cide, ckds, ckde, cjds, cjde, &
571	cims, cime, ckms, ckme, cjms, cjme, &
572	cits, cite, ckts, ckte, cjts, cjte, &
573	nfld, & ! ND field
574	nids, nide, nkds, nkde, njds, njde, &
575	nims, nime, nkms, nkme, njms, njme, &
576	nits, nite, nkts, nkte, njts, njte, &
577	shw, & ! stencil half width for interp
578	imask, & ! interpolation mask
579	xstag, ystag, & ! staggering of field
580	ipos, jpos, & ! Position of lower left of nest in CD
581	nri, nrj ) ! nest ratios
582	USE module_configure
583	USE module_wrf_error
584	IMPLICIT NONE
585
586
587	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
588	cims, cime, ckms, ckme, cjms, cjme, &
589	cits, cite, ckts, ckte, cjts, cjte, &
590	nids, nide, nkds, nkde, njds, njde, &
591	nims, nime, nkms, nkme, njms, njme, &
592	nits, nite, nkts, nkte, njts, njte, &
593	shw, &
594	ipos, jpos, &
595	nri, nrj
596	LOGICAL, INTENT(IN) :: xstag, ystag
597
598	INTEGER, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ), INTENT(OUT) :: cfld
599	INTEGER, DIMENSION ( nims:nime, nkms:nkme, njms:njme ), INTENT(IN) :: nfld
600	INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: imask
601
602	! Local
603
604	INTEGER ci, cj, ck, ni, nj, nk, ip, jp, ioff, joff, ioffa, joffa
605	INTEGER :: icmin,icmax,jcmin,jcmax
606	INTEGER :: istag,jstag, ipoints,jpoints,ijpoints
607	INTEGER , PARAMETER :: passes = 2
608	INTEGER spec_zone
609
610	CALL nl_get_spec_zone( 1, spec_zone )
611	istag = 1 ; jstag = 1
612	IF ( xstag ) istag = 0
613	IF ( ystag ) jstag = 0
614
615	IF( MOD(nrj,2) .NE. 0) THEN ! odd refinement ratio
616
617	DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
618	nj = (cj-jpos)*nrj + jstag + 1
619	DO ck = ckts, ckte
620	nk = ck
621	DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
622	ni = (ci-ipos)*nri + istag + 1
623	cfld( ci, ck, cj ) = nfld( ni , nk , nj )
624	ENDDO
625	ENDDO
626	ENDDO
627
628	ELSE ! even refinement ratio
629	DO cj = MAX(jpos+spec_zone,cjts),MIN(jpos+(njde-njds)/nrj-jstag-spec_zone,cjte)
630	nj = (cj-jpos)*nrj + 1
631	DO ck = ckts, ckte
632	nk = ck
633	DO ci = MAX(ipos+spec_zone,cits),MIN(ipos+(nide-nids)/nri-istag-spec_zone,cite)
634	ni = (ci-ipos)*nri + 1
635	ipoints = nri/2 -1
636	jpoints = nrj/2 -1
637	cfld( ci, ck, cj ) = nfld( ni+ipoints , nk , nj+jpoints )
638	END DO
639	END DO
640	END DO
641
642	END IF
643
644	RETURN
645
646	END SUBROUTINE copy_fcni
647
648	!==================================
649
650	SUBROUTINE bdy_interp ( cfld, & ! CD field
651	cids, cide, ckds, ckde, cjds, cjde, &
652	cims, cime, ckms, ckme, cjms, cjme, &
653	cits, cite, ckts, ckte, cjts, cjte, &
654	nfld, & ! ND field
655	nids, nide, nkds, nkde, njds, njde, &
656	nims, nime, nkms, nkme, njms, njme, &
657	nits, nite, nkts, nkte, njts, njte, &
658	shw, & ! stencil half width
659	imask, & ! interpolation mask
660	xstag, ystag, & ! staggering of field
661	ipos, jpos, & ! Position of lower left of nest in CD
662	nri, nrj, & ! nest ratios
663	cbdy_xs, nbdy_xs, &
664	cbdy_xe, nbdy_xe, &
665	cbdy_ys, nbdy_ys, &
666	cbdy_ye, nbdy_ye, &
667	cbdy_txs, nbdy_txs, &
668	cbdy_txe, nbdy_txe, &
669	cbdy_tys, nbdy_tys, &
670	cbdy_tye, nbdy_tye, &
671	cdt, ndt &
672	) ! boundary arrays
673	USE module_configure
674	IMPLICIT NONE
675
676	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
677	cims, cime, ckms, ckme, cjms, cjme, &
678	cits, cite, ckts, ckte, cjts, cjte, &
679	nids, nide, nkds, nkde, njds, njde, &
680	nims, nime, nkms, nkme, njms, njme, &
681	nits, nite, nkts, nkte, njts, njte, &
682	shw, &
683	ipos, jpos, &
684	nri, nrj
685
686	LOGICAL, INTENT(IN) :: xstag, ystag
687
688	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
689	REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
690	INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
691	REAL, DIMENSION( * ), INTENT(INOUT) :: cbdy_xs, cbdy_txs, nbdy_xs, nbdy_txs
692	REAL, DIMENSION( * ), INTENT(INOUT) :: cbdy_xe, cbdy_txe, nbdy_xe, nbdy_txe
693	REAL, DIMENSION( * ), INTENT(INOUT) :: cbdy_ys, cbdy_tys, nbdy_ys, nbdy_tys
694	REAL, DIMENSION( * ), INTENT(INOUT) :: cbdy_ye, cbdy_tye, nbdy_ye, nbdy_tye
695	REAL cdt, ndt
696
697	! Local
698
699	INTEGER nijds, nijde, spec_bdy_width
700
701	nijds = min(nids, njds)
702	nijde = max(nide, njde)
703	CALL nl_get_spec_bdy_width( 1, spec_bdy_width )
704
705	CALL bdy_interp1( cfld, & ! CD field
706	cids, cide, ckds, ckde, cjds, cjde, &
707	cims, cime, ckms, ckme, cjms, cjme, &
708	cits, cite, ckts, ckte, cjts, cjte, &
709	nfld, & ! ND field
710	nijds, nijde , spec_bdy_width , &
711	nids, nide, nkds, nkde, njds, njde, &
712	nims, nime, nkms, nkme, njms, njme, &
713	nits, nite, nkts, nkte, njts, njte, &
714	shw, imask, &
715	xstag, ystag, & ! staggering of field
716	ipos, jpos, & ! Position of lower left of nest in CD
717	nri, nrj, &
718	cbdy_xs, nbdy_xs, &
719	cbdy_xe, nbdy_xe, &
720	cbdy_ys, nbdy_ys, &
721	cbdy_ye, nbdy_ye, &
722	cbdy_txs, nbdy_txs, &
723	cbdy_txe, nbdy_txe, &
724	cbdy_tys, nbdy_tys, &
725	cbdy_tye, nbdy_tye, &
726	cdt, ndt &
727	)
728
729	RETURN
730
731	END SUBROUTINE bdy_interp
732
733	SUBROUTINE bdy_interp1( cfld, & ! CD field
734	cids, cide, ckds, ckde, cjds, cjde, &
735	cims, cime, ckms, ckme, cjms, cjme, &
736	cits, cite, ckts, ckte, cjts, cjte, &
737	nfld, & ! ND field
738	nijds, nijde, spec_bdy_width , &
739	nids, nide, nkds, nkde, njds, njde, &
740	nims, nime, nkms, nkme, njms, njme, &
741	nits, nite, nkts, nkte, njts, njte, &
742	shw1, &
743	imask, & ! interpolation mask
744	xstag, ystag, & ! staggering of field
745	ipos, jpos, & ! Position of lower left of nest in CD
746	nri, nrj, &
747	cbdy_xs, bdy_xs, &
748	cbdy_xe, bdy_xe, &
749	cbdy_ys, bdy_ys, &
750	cbdy_ye, bdy_ye, &
751	cbdy_txs, bdy_txs, &
752	cbdy_txe, bdy_txe, &
753	cbdy_tys, bdy_tys, &
754	cbdy_tye, bdy_tye, &
755	cdt, ndt &
756	)
757
758	USE module_configure
759	use module_state_description
760	IMPLICIT NONE
761
762	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
763	cims, cime, ckms, ckme, cjms, cjme, &
764	cits, cite, ckts, ckte, cjts, cjte, &
765	nids, nide, nkds, nkde, njds, njde, &
766	nims, nime, nkms, nkme, njms, njme, &
767	nits, nite, nkts, nkte, njts, njte, &
768	shw1, & ! ignore
769	ipos, jpos, &
770	nri, nrj
771	INTEGER, INTENT(IN) :: nijds, nijde, spec_bdy_width
772	LOGICAL, INTENT(IN) :: xstag, ystag
773
774	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ), INTENT(INOUT) :: cfld
775	REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ), INTENT(INOUT) :: nfld
776	INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
777	REAL, DIMENSION ( * ), INTENT(INOUT) :: cbdy_xs, cbdy_txs ! not used
778	REAL, DIMENSION ( * ), INTENT(INOUT) :: cbdy_xe, cbdy_txe ! not used
779	REAL, DIMENSION ( * ), INTENT(INOUT) :: cbdy_ys, cbdy_tys ! not used
780	REAL, DIMENSION ( * ), INTENT(INOUT) :: cbdy_ye, cbdy_tye ! not used
781	REAL :: cdt, ndt
782	REAL, DIMENSION ( njms:njme, nkms:nkme, spec_bdy_width ), INTENT(INOUT) :: bdy_xs, bdy_txs
783	REAL, DIMENSION ( njms:njme, nkms:nkme, spec_bdy_width ), INTENT(INOUT) :: bdy_xe, bdy_txe
784	REAL, DIMENSION ( nims:nime, nkms:nkme, spec_bdy_width ), INTENT(INOUT) :: bdy_ys, bdy_tys
785	REAL, DIMENSION ( nims:nime, nkms:nkme, spec_bdy_width ), INTENT(INOUT) :: bdy_ye, bdy_tye
786
787	! Local
788
789	REAL*8 rdt
790	INTEGER ci, cj, ck, ni, nj, nk, ni1, nj1, nk1, ip, jp, ioff, joff
791	#ifdef MM5_SINT
792	INTEGER nfx, ior
793	PARAMETER (ior=2)
794	INTEGER nf
795	REAL psca1(cims:cime,cjms:cjme,nri*nrj)
796	REAL psca(cims:cime,cjms:cjme,nri*nrj)
797	LOGICAL icmask( cims:cime, cjms:cjme )
798	INTEGER i,j,k
799	#endif
800	INTEGER shw
801	INTEGER spec_zone
802	INTEGER relax_zone
803	INTEGER sz
804	INTEGER n2ci,n
805	INTEGER n2cj
806
807	! statement functions for converting a nest index to coarse
808	n2ci(n) = (n+ipos*nri-1)/nri
809	n2cj(n) = (n+jpos*nrj-1)/nrj
810
811	rdt = 1.D0/cdt
812
813	shw = 0
814
815	ioff = 0 ; joff = 0
816	IF ( xstag ) ioff = (nri-1)/2
817	IF ( ystag ) joff = (nrj-1)/2
818
819	! Iterate over the ND tile and compute the values
820	! from the CD tile.
821
822	#ifdef MM5_SINT
823	CALL nl_get_spec_zone( 1, spec_zone )
824	CALL nl_get_relax_zone( 1, relax_zone )
825	sz = MIN(MAX( spec_zone, relax_zone + 1 ),spec_bdy_width)
826
827	nfx = nri * nrj
828
829	!$OMP PARALLEL DO &
830	!$OMP PRIVATE ( i,j,k,ni,nj,ni1,nj1,ci,cj,ip,jp,nk,ck,nf,icmask,psca,psca1 )
831	DO k = ckts, ckte
832
833	DO nf = 1,nfx
834	DO j = cjms,cjme
835	nj = (j-jpos) * nrj + ( nrj / 2 + 1 ) ! j point on nest
836	DO i = cims,cime
837	ni = (i-ipos) * nri + ( nri / 2 + 1 ) ! i point on nest
838	psca1(i,j,nf) = cfld(i,k,j)
839	ENDDO
840	ENDDO
841	ENDDO
842	! hopefully less ham handed but still correct and more efficient
843	! sintb ignores icmask so it does not matter that icmask is not set
844	!
845	! SOUTH BDY
846	IF ( njts .ge. njds .and. njts .le. njds + sz + joff ) THEN
847	CALL sintb( psca1, psca, &
848	cims, cime, cjms, cjme, icmask, &
849	n2ci(nits)-1, n2ci(nite)+1, n2cj(MAX(njts,njds)), n2cj(MIN(njte,njds+sz+joff)), nrj*nri, xstag, ystag )
850	ENDIF
851	! NORTH BDY
852	IF ( njte .le. njde .and. njte .ge. njde - sz - joff ) THEN
853	CALL sintb( psca1, psca, &
854	cims, cime, cjms, cjme, icmask, &
855	n2ci(nits)-1, n2ci(nite)+1, n2cj(MAX(njts,njde-sz-joff)), n2cj(MIN(njte,njde-1+joff)), nrj*nri, xstag, ystag )
856	ENDIF
857	! WEST BDY
858	IF ( nits .ge. nids .and. nits .le. nids + sz + ioff ) THEN
859	CALL sintb( psca1, psca, &
860	cims, cime, cjms, cjme, icmask, &
861	n2ci(MAX(nits,nids)), n2ci(MIN(nite,nids+sz+ioff)), n2cj(njts)-1, n2cj(njte)+1, nrj*nri, xstag, ystag )
862	ENDIF
863	! EAST BDY
864	IF ( nite .le. nide .and. nite .ge. nide - sz - ioff ) THEN
865	CALL sintb( psca1, psca, &
866	cims, cime, cjms, cjme, icmask, &
867	n2ci(MAX(nits,nide-sz-ioff)), n2ci(MIN(nite,nide-1+ioff)), n2cj(njts)-1, n2cj(njte)+1, nrj*nri, xstag, ystag )
868	ENDIF
869
870	DO nj1 = MAX(njds,njts-1), MIN(njde+joff,njte+joff+1)
871	cj = jpos + (nj1-1) / nrj ! j coord of CD point
872	jp = mod ( nj1-1 , nrj ) ! coord of ND w/i CD point
873	nk = k
874	ck = nk
875	DO ni1 = MAX(nids,nits-1), MIN(nide+ioff,nite+ioff+1)
876	ci = ipos + (ni1-1) / nri ! j coord of CD point
877	ip = mod ( ni1-1 , nri ) ! coord of ND w/i CD point
878
879	ni = ni1-ioff
880	nj = nj1-joff
881
882	IF ( ( ni.LT.nids) .OR. (nj.LT.njds) ) THEN
883	CYCLE
884	END IF
885
886	!bdy contains the value at t-dt. psca contains the value at t
887	!compute dv/dt and store in bdy_t
888	!afterwards store the new value of v at t into bdy
889	! WEST
890	IF ( ni .ge. nids .and. ni .lt. nids + sz ) THEN
891	bdy_txs( nj,k,ni ) = rdt(psca(ci+shw,cj+shw,ip+1+(jp)nri)-nfld(ni,k,nj))
892	bdy_xs( nj,k,ni ) = psca(ci+shw,cj+shw,ip+1+(jp)*nri )
893	ENDIF
894
895	! SOUTH
896	IF ( nj .ge. njds .and. nj .lt. njds + sz ) THEN
897	bdy_tys( ni,k,nj ) = rdt(psca(ci+shw,cj+shw,ip+1+(jp)nri)-nfld(ni,k,nj))
898	bdy_ys( ni,k,nj ) = psca(ci+shw,cj+shw,ip+1+(jp)*nri )
899	ENDIF
900
901	! EAST
902	IF ( xstag ) THEN
903	IF ( ni .ge. nide - sz + 1 .AND. ni .le. nide ) THEN
904	bdy_txe( nj,k,nide-ni+1 ) = rdt(psca(ci+shw,cj+shw,ip+1+(jp)nri)-nfld(ni,k,nj))
905	bdy_xe( nj,k,nide-ni+1 ) = psca(ci+shw,cj+shw,ip+1+(jp)*nri )
906	ENDIF
907	ELSE
908	IF ( ni .ge. nide - sz .AND. ni .le. nide-1 ) THEN
909	bdy_txe( nj,k,nide-ni ) = rdt(psca(ci+shw,cj+shw,ip+1+(jp)nri)-nfld(ni,k,nj))
910	bdy_xe( nj,k,nide-ni ) = psca(ci+shw,cj+shw,ip+1+(jp)*nri )
911	ENDIF
912	ENDIF
913
914	! NORTH
915	IF ( ystag ) THEN
916	IF ( nj .ge. njde - sz + 1 .AND. nj .le. njde ) THEN
917	bdy_tye( ni,k,njde-nj+1 ) = rdt(psca(ci+shw,cj+shw,ip+1+(jp)nri)-nfld(ni,k,nj))
918	bdy_ye( ni,k,njde-nj+1 ) = psca(ci+shw,cj+shw,ip+1+(jp)*nri )
919	ENDIF
920	ELSE
921	IF ( nj .ge. njde - sz .AND. nj .le. njde-1 ) THEN
922	bdy_tye(ni,k,njde-nj ) = rdt(psca(ci+shw,cj+shw,ip+1+(jp)nri)-nfld(ni,k,nj))
923	bdy_ye( ni,k,njde-nj ) = psca(ci+shw,cj+shw,ip+1+(jp)*nri )
924	ENDIF
925	ENDIF
926
927	ENDDO
928	ENDDO
929	ENDDO
930	!$OMP END PARALLEL DO
931	#endif
932
933	RETURN
934
935	END SUBROUTINE bdy_interp1
936
937
938
939	SUBROUTINE interp_fcni( cfld, & ! CD field
940	cids, cide, ckds, ckde, cjds, cjde, &
941	cims, cime, ckms, ckme, cjms, cjme, &
942	cits, cite, ckts, ckte, cjts, cjte, &
943	nfld, & ! ND field
944	nids, nide, nkds, nkde, njds, njde, &
945	nims, nime, nkms, nkme, njms, njme, &
946	nits, nite, nkts, nkte, njts, njte, &
947	shw, & ! stencil half width
948	imask, & ! interpolation mask
949	xstag, ystag, & ! staggering of field
950	ipos, jpos, & ! Position of lower left of nest in CD
951	nri, nrj ) ! nest ratios
952	USE module_configure
953	IMPLICIT NONE
954
955
956	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
957	cims, cime, ckms, ckme, cjms, cjme, &
958	cits, cite, ckts, ckte, cjts, cjte, &
959	nids, nide, nkds, nkde, njds, njde, &
960	nims, nime, nkms, nkme, njms, njme, &
961	nits, nite, nkts, nkte, njts, njte, &
962	shw, &
963	ipos, jpos, &
964	nri, nrj
965	LOGICAL, INTENT(IN) :: xstag, ystag
966
967	INTEGER, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
968	INTEGER, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
969	INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
970
971	! Local
972
973	INTEGER ci, cj, ck, ni, nj, nk, ip, jp
974
975	! Iterate over the ND tile and compute the values
976	! from the CD tile.
977
978	!write(0,'("cits:cite, ckts:ckte, cjts:cjte ",6i4)')cits,cite, ckts,ckte, cjts,cjte
979	!write(0,'("nits:nite, nkts:nkte, njts:njte ",6i4)')nits,nite, nkts,nkte, njts,njte
980
981	DO nj = njts, njte
982	cj = jpos + (nj-1) / nrj ! j coord of CD point
983	jp = mod ( nj , nrj ) ! coord of ND w/i CD point
984	DO nk = nkts, nkte
985	ck = nk
986	DO ni = nits, nite
987	ci = ipos + (ni-1) / nri ! j coord of CD point
988	ip = mod ( ni , nri ) ! coord of ND w/i CD point
989	! This is a trivial implementation of the interp_fcn; just copies
990	! the values from the CD into the ND
991	nfld( ni, nk, nj ) = cfld( ci , ck , cj )
992	ENDDO
993	ENDDO
994	ENDDO
995
996	RETURN
997
998	END SUBROUTINE interp_fcni
999
1000	SUBROUTINE interp_fcnm( cfld, & ! CD field
1001	cids, cide, ckds, ckde, cjds, cjde, &
1002	cims, cime, ckms, ckme, cjms, cjme, &
1003	cits, cite, ckts, ckte, cjts, cjte, &
1004	nfld, & ! ND field
1005	nids, nide, nkds, nkde, njds, njde, &
1006	nims, nime, nkms, nkme, njms, njme, &
1007	nits, nite, nkts, nkte, njts, njte, &
1008	shw, & ! stencil half width
1009	imask, & ! interpolation mask
1010	xstag, ystag, & ! staggering of field
1011	ipos, jpos, & ! Position of lower left of nest in CD
1012	nri, nrj ) ! nest ratios
1013	USE module_configure
1014	IMPLICIT NONE
1015
1016
1017	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
1018	cims, cime, ckms, ckme, cjms, cjme, &
1019	cits, cite, ckts, ckte, cjts, cjte, &
1020	nids, nide, nkds, nkde, njds, njde, &
1021	nims, nime, nkms, nkme, njms, njme, &
1022	nits, nite, nkts, nkte, njts, njte, &
1023	shw, &
1024	ipos, jpos, &
1025	nri, nrj
1026	LOGICAL, INTENT(IN) :: xstag, ystag
1027
1028	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
1029	REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
1030	INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
1031
1032	! Local
1033
1034	INTEGER ci, cj, ck, ni, nj, nk, ip, jp
1035
1036	! Iterate over the ND tile and compute the values
1037	! from the CD tile.
1038
1039	!write(0,'("mask cits:cite, ckts:ckte, cjts:cjte ",6i4)')cits,cite, ckts,ckte, cjts,cjte
1040	!write(0,'("mask nits:nite, nkts:nkte, njts:njte ",6i4)')nits,nite, nkts,nkte, njts,njte
1041
1042	DO nj = njts, njte
1043	cj = jpos + (nj-1) / nrj ! j coord of CD point
1044	jp = mod ( nj , nrj ) ! coord of ND w/i CD point
1045	DO nk = nkts, nkte
1046	ck = nk
1047	DO ni = nits, nite
1048	ci = ipos + (ni-1) / nri ! j coord of CD point
1049	ip = mod ( ni , nri ) ! coord of ND w/i CD point
1050	! This is a trivial implementation of the interp_fcn; just copies
1051	! the values from the CD into the ND
1052	nfld( ni, nk, nj ) = cfld( ci , ck , cj )
1053	ENDDO
1054	ENDDO
1055	ENDDO
1056
1057	RETURN
1058
1059	END SUBROUTINE interp_fcnm
1060
1061	SUBROUTINE interp_mask_land_field ( enable, & ! says whether to allow interpolation or just the bcasts
1062	cfld, & ! CD field
1063	cids, cide, ckds, ckde, cjds, cjde, &
1064	cims, cime, ckms, ckme, cjms, cjme, &
1065	cits, cite, ckts, ckte, cjts, cjte, &
1066	nfld, & ! ND field
1067	nids, nide, nkds, nkde, njds, njde, &
1068	nims, nime, nkms, nkme, njms, njme, &
1069	nits, nite, nkts, nkte, njts, njte, &
1070	shw, & ! stencil half width
1071	imask, & ! interpolation mask
1072	xstag, ystag, & ! staggering of field
1073	ipos, jpos, & ! Position of lower left of nest in CD
1074	nri, nrj, & ! nest ratios
1075	clu, nlu )
1076
1077	USE module_configure
1078	USE module_wrf_error
1079
1080	IMPLICIT NONE
1081
1082
1083	LOGICAL, INTENT(IN) :: enable
1084	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
1085	cims, cime, ckms, ckme, cjms, cjme, &
1086	cits, cite, ckts, ckte, cjts, cjte, &
1087	nids, nide, nkds, nkde, njds, njde, &
1088	nims, nime, nkms, nkme, njms, njme, &
1089	nits, nite, nkts, nkte, njts, njte, &
1090	shw, &
1091	ipos, jpos, &
1092	nri, nrj
1093	LOGICAL, INTENT(IN) :: xstag, ystag
1094
1095	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
1096	REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
1097	INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
1098
1099	REAL, DIMENSION ( cims:cime, cjms:cjme ) :: clu
1100	REAL, DIMENSION ( nims:nime, njms:njme ) :: nlu
1101
1102	! Local
1103
1104	INTEGER ci, cj, ck, ni, nj, nk, ip, jp
1105	INTEGER :: icount , ii , jj , ist , ien , jst , jen , iswater
1106	REAL :: avg , sum , dx , dy
1107	INTEGER , PARAMETER :: max_search = 5
1108	CHARACTER*120 message
1109
1110	! Find out what the water value is.
1111
1112	CALL nl_get_iswater(1,iswater)
1113
1114	! Right now, only mass point locations permitted.
1115
1116	IF ( ( .NOT. xstag) .AND. ( .NOT. ystag ) ) THEN
1117
1118	! Loop over each i,k,j in the nested domain.
1119
1120	IF ( enable ) THEN
1121
1122	DO nj = njts, njte
1123	IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
1124	cj = ( nj + (nrj/2)-1 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
1125	ELSE
1126	cj = ( nj + (nrj-1)/2 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
1127	END IF
1128	DO nk = nkts, nkte
1129	ck = nk
1130	DO ni = nits, nite
1131	IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
1132	ci = ( ni + (nri/2)-1 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
1133	ELSE
1134	ci = ( ni + (nri-1)/2 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
1135	END IF
1136
1137
1138
1139
1140	!
1141	! (ci,cj+1) (ci+1,cj+1)
1142	! - -------------
1143	! 1-dy \| \| \|
1144	! \| \| \|
1145	! - \| * \|
1146	! dy \| \| (ni,nj) \|
1147	! \| \| \|
1148	! - -------------
1149	! (ci,cj) (ci+1,cj)
1150	!
1151	! \|--\|--------\|
1152	! dx 1-dx
1153
1154
1155	! For odd ratios, at (nri+1)/2, we are on the coarse grid point, so dx = 0
1156
1157	IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
1158	dx = ( REAL ( MOD ( ni+(nri-1)/2 , nri ) ) + 0.5 ) / REAL ( nri )
1159	ELSE
1160	dx = REAL ( MOD ( ni+(nri-1)/2 , nri ) ) / REAL ( nri )
1161	END IF
1162	IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
1163	dy = ( REAL ( MOD ( nj+(nrj-1)/2 , nrj ) ) + 0.5 ) / REAL ( nrj )
1164	ELSE
1165	dy = REAL ( MOD ( nj+(nrj-1)/2 , nrj ) ) / REAL ( nrj )
1166	END IF
1167
1168	! This is a "land only" field. If this is a water point, no operations required.
1169
1170	IF ( ( NINT(nlu(ni ,nj )) .EQ. iswater ) ) THEN
1171	! noop
1172	! nfld(ni,nk,nj) = 1.e20
1173	nfld(ni,nk,nj) = -1
1174
1175	! If this is a nested land point, and the surrounding coarse values are all land points,
1176	! then this is a simple 4-pt interpolation.
1177
1178	ELSE IF ( ( NINT(nlu(ni ,nj )) .NE. iswater ) .AND. &
1179	( NINT(clu(ci ,cj )) .NE. iswater ) .AND. &
1180	( NINT(clu(ci+1,cj )) .NE. iswater ) .AND. &
1181	( NINT(clu(ci ,cj+1)) .NE. iswater ) .AND. &
1182	( NINT(clu(ci+1,cj+1)) .NE. iswater ) ) THEN
1183	nfld(ni,nk,nj) = ( 1. - dx ) * ( ( 1. - dy ) * cfld(ci ,ck,cj ) + &
1184	dy * cfld(ci ,ck,cj+1) ) + &
1185	dx * ( ( 1. - dy ) * cfld(ci+1,ck,cj ) + &
1186	dy * cfld(ci+1,ck,cj+1) )
1187
1188	! If this is a nested land point and there are NO coarse land values surrounding,
1189	! we temporarily punt.
1190
1191	ELSE IF ( ( NINT(nlu(ni ,nj )) .NE. iswater ) .AND. &
1192	( NINT(clu(ci ,cj )) .EQ. iswater ) .AND. &
1193	( NINT(clu(ci+1,cj )) .EQ. iswater ) .AND. &
1194	( NINT(clu(ci ,cj+1)) .EQ. iswater ) .AND. &
1195	( NINT(clu(ci+1,cj+1)) .EQ. iswater ) ) THEN
1196	! nfld(ni,nk,nj) = -1.e20
1197	nfld(ni,nk,nj) = -1
1198
1199	! If there are some water points and some land points, take an average.
1200
1201	ELSE IF ( NINT(nlu(ni ,nj )) .NE. iswater ) THEN
1202	icount = 0
1203	sum = 0
1204	IF ( NINT(clu(ci ,cj )) .NE. iswater ) THEN
1205	icount = icount + 1
1206	sum = sum + cfld(ci ,ck,cj )
1207	END IF
1208	IF ( NINT(clu(ci+1,cj )) .NE. iswater ) THEN
1209	icount = icount + 1
1210	sum = sum + cfld(ci+1,ck,cj )
1211	END IF
1212	IF ( NINT(clu(ci ,cj+1)) .NE. iswater ) THEN
1213	icount = icount + 1
1214	sum = sum + cfld(ci ,ck,cj+1)
1215	END IF
1216	IF ( NINT(clu(ci+1,cj+1)) .NE. iswater ) THEN
1217	icount = icount + 1
1218	sum = sum + cfld(ci+1,ck,cj+1)
1219	END IF
1220	nfld(ni,nk,nj) = sum / REAL ( icount )
1221	END IF
1222	END DO
1223	END DO
1224	END DO
1225
1226
1227	! Get an average of the whole domain for problem locations.
1228
1229	sum = 0
1230	icount = 0
1231	DO nj = njts, njte
1232	DO nk = nkts, nkte
1233	DO ni = nits, nite
1234	IF ( ( nfld(ni,nk,nj) .GT. -1.e19 ) .AND. ( nfld(ni,nk,nj) .LT. 1.e19 ) ) THEN
1235	icount = icount + 1
1236	sum = sum + nfld(ni,nk,nj)
1237	END IF
1238	END DO
1239	END DO
1240	END DO
1241	ELSE
1242	sum = 0.
1243	icount = 0
1244	ENDIF
1245	CALL wrf_dm_bcast_real( sum, 1 )
1246	CALL wrf_dm_bcast_integer( icount, 1 )
1247	IF ( enable ) THEN
1248	IF ( icount .GT. 0 ) THEN
1249	avg = sum / REAL ( icount )
1250
1251	! OK, if there were any of those island situations, we try to search a bit broader
1252	! of an area in the coarse grid.
1253
1254	DO nj = njts, njte
1255	DO nk = nkts, nkte
1256	DO ni = nits, nite
1257	IF ( nfld(ni,nk,nj) .LT. -1.e19 ) THEN
1258	IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
1259	cj = ( nj + (nrj/2)-1 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
1260	ELSE
1261	cj = ( nj + (nrj-1)/2 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
1262	END IF
1263	IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
1264	ci = ( ni + (nri/2)-1 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
1265	ELSE
1266	ci = ( ni + (nri-1)/2 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
1267	END IF
1268	ist = MAX (ci-max_search,cits)
1269	ien = MIN (ci+max_search,cite,cide-1)
1270	jst = MAX (cj-max_search,cjts)
1271	jen = MIN (cj+max_search,cjte,cjde-1)
1272	icount = 0
1273	sum = 0
1274	DO jj = jst,jen
1275	DO ii = ist,ien
1276	IF ( NINT(clu(ii,jj)) .NE. iswater ) THEN
1277	icount = icount + 1
1278	sum = sum + cfld(ii,nk,jj)
1279	END IF
1280	END DO
1281	END DO
1282	IF ( icount .GT. 0 ) THEN
1283	nfld(ni,nk,nj) = sum / REAL ( icount )
1284	ELSE
1285	! CALL wrf_error_fatal ( "horizontal interp error - island" )
1286	write(message,*) 'horizontal interp error - island, using average ', avg
1287	CALL wrf_message ( message )
1288	nfld(ni,nk,nj) = avg
1289	END IF
1290	END IF
1291	END DO
1292	END DO
1293	END DO
1294	ENDIF
1295	ENDIF
1296	ELSE
1297	CALL wrf_error_fatal ( "only unstaggered fields right now" )
1298	END IF
1299
1300	END SUBROUTINE interp_mask_land_field
1301
1302	SUBROUTINE interp_mask_water_field ( enable, & ! says whether to allow interpolation or just the bcasts
1303	cfld, & ! CD field
1304	cids, cide, ckds, ckde, cjds, cjde, &
1305	cims, cime, ckms, ckme, cjms, cjme, &
1306	cits, cite, ckts, ckte, cjts, cjte, &
1307	nfld, & ! ND field
1308	nids, nide, nkds, nkde, njds, njde, &
1309	nims, nime, nkms, nkme, njms, njme, &
1310	nits, nite, nkts, nkte, njts, njte, &
1311	shw, & ! stencil half width
1312	imask, & ! interpolation mask
1313	xstag, ystag, & ! staggering of field
1314	ipos, jpos, & ! Position of lower left of nest in CD
1315	nri, nrj, & ! nest ratios
1316	clu, nlu )
1317
1318	USE module_configure
1319	USE module_wrf_error
1320
1321	IMPLICIT NONE
1322
1323
1324	LOGICAL, INTENT(IN) :: enable
1325	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
1326	cims, cime, ckms, ckme, cjms, cjme, &
1327	cits, cite, ckts, ckte, cjts, cjte, &
1328	nids, nide, nkds, nkde, njds, njde, &
1329	nims, nime, nkms, nkme, njms, njme, &
1330	nits, nite, nkts, nkte, njts, njte, &
1331	shw, &
1332	ipos, jpos, &
1333	nri, nrj
1334	LOGICAL, INTENT(IN) :: xstag, ystag
1335
1336	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
1337	REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
1338	INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
1339
1340	REAL, DIMENSION ( cims:cime, cjms:cjme ) :: clu
1341	REAL, DIMENSION ( nims:nime, njms:njme ) :: nlu
1342
1343	! Local
1344
1345	INTEGER ci, cj, ck, ni, nj, nk, ip, jp
1346	INTEGER :: icount , ii , jj , ist , ien , jst , jen , iswater
1347	REAL :: avg , sum , dx , dy
1348	INTEGER , PARAMETER :: max_search = 5
1349
1350	! Find out what the water value is.
1351
1352	CALL nl_get_iswater(1,iswater)
1353
1354	! Right now, only mass point locations permitted.
1355
1356	IF ( ( .NOT. xstag) .AND. ( .NOT. ystag ) ) THEN
1357
1358	IF ( enable ) THEN
1359	! Loop over each i,k,j in the nested domain.
1360
1361	DO nj = njts, njte
1362	IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
1363	cj = ( nj + (nrj/2)-1 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
1364	ELSE
1365	cj = ( nj + (nrj-1)/2 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
1366	END IF
1367	DO nk = nkts, nkte
1368	ck = nk
1369	DO ni = nits, nite
1370	IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
1371	ci = ( ni + (nri/2)-1 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
1372	ELSE
1373	ci = ( ni + (nri-1)/2 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
1374	END IF
1375
1376
1377
1378
1379	!
1380	! (ci,cj+1) (ci+1,cj+1)
1381	! - -------------
1382	! 1-dy \| \| \|
1383	! \| \| \|
1384	! - \| * \|
1385	! dy \| \| (ni,nj) \|
1386	! \| \| \|
1387	! - -------------
1388	! (ci,cj) (ci+1,cj)
1389	!
1390	! \|--\|--------\|
1391	! dx 1-dx
1392
1393
1394	! At ni=2, we are on the coarse grid point, so dx = 0
1395
1396	IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
1397	dx = ( REAL ( MOD ( ni+(nri-1)/2 , nri ) ) + 0.5 ) / REAL ( nri )
1398	ELSE
1399	dx = REAL ( MOD ( ni+(nri-1)/2 , nri ) ) / REAL ( nri )
1400	END IF
1401	IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
1402	dy = ( REAL ( MOD ( nj+(nrj-1)/2 , nrj ) ) + 0.5 ) / REAL ( nrj )
1403	ELSE
1404	dy = REAL ( MOD ( nj+(nrj-1)/2 , nrj ) ) / REAL ( nrj )
1405	END IF
1406
1407	! This is a "water only" field. If this is a land point, no operations required.
1408
1409	IF ( ( NINT(nlu(ni ,nj )) .NE. iswater ) ) THEN
1410	! noop
1411	! nfld(ni,nk,nj) = 1.e20
1412	nfld(ni,nk,nj) = -1
1413
1414	! If this is a nested water point, and the surrounding coarse values are all water points,
1415	! then this is a simple 4-pt interpolation.
1416
1417	ELSE IF ( ( NINT(nlu(ni ,nj )) .EQ. iswater ) .AND. &
1418	( NINT(clu(ci ,cj )) .EQ. iswater ) .AND. &
1419	( NINT(clu(ci+1,cj )) .EQ. iswater ) .AND. &
1420	( NINT(clu(ci ,cj+1)) .EQ. iswater ) .AND. &
1421	( NINT(clu(ci+1,cj+1)) .EQ. iswater ) ) THEN
1422	nfld(ni,nk,nj) = ( 1. - dx ) * ( ( 1. - dy ) * cfld(ci ,ck,cj ) + &
1423	dy * cfld(ci ,ck,cj+1) ) + &
1424	dx * ( ( 1. - dy ) * cfld(ci+1,ck,cj ) + &
1425	dy * cfld(ci+1,ck,cj+1) )
1426
1427	! If this is a nested water point and there are NO coarse water values surrounding,
1428	! we temporarily punt.
1429
1430	ELSE IF ( ( NINT(nlu(ni ,nj )) .EQ. iswater ) .AND. &
1431	( NINT(clu(ci ,cj )) .NE. iswater ) .AND. &
1432	( NINT(clu(ci+1,cj )) .NE. iswater ) .AND. &
1433	( NINT(clu(ci ,cj+1)) .NE. iswater ) .AND. &
1434	( NINT(clu(ci+1,cj+1)) .NE. iswater ) ) THEN
1435	! nfld(ni,nk,nj) = -1.e20
1436	nfld(ni,nk,nj) = -1
1437
1438	! If there are some land points and some water points, take an average.
1439
1440	ELSE IF ( NINT(nlu(ni ,nj )) .EQ. iswater ) THEN
1441	icount = 0
1442	sum = 0
1443	IF ( NINT(clu(ci ,cj )) .EQ. iswater ) THEN
1444	icount = icount + 1
1445	sum = sum + cfld(ci ,ck,cj )
1446	END IF
1447	IF ( NINT(clu(ci+1,cj )) .EQ. iswater ) THEN
1448	icount = icount + 1
1449	sum = sum + cfld(ci+1,ck,cj )
1450	END IF
1451	IF ( NINT(clu(ci ,cj+1)) .EQ. iswater ) THEN
1452	icount = icount + 1
1453	sum = sum + cfld(ci ,ck,cj+1)
1454	END IF
1455	IF ( NINT(clu(ci+1,cj+1)) .EQ. iswater ) THEN
1456	icount = icount + 1
1457	sum = sum + cfld(ci+1,ck,cj+1)
1458	END IF
1459	nfld(ni,nk,nj) = sum / REAL ( icount )
1460	END IF
1461	END DO
1462	END DO
1463	END DO
1464
1465	! Get an average of the whole domain for problem locations.
1466
1467	sum = 0
1468	icount = 0
1469	DO nj = njts, njte
1470	DO nk = nkts, nkte
1471	DO ni = nits, nite
1472	IF ( ( nfld(ni,nk,nj) .GT. -1.e19 ) .AND. ( nfld(ni,nk,nj) .LT. 1.e19 ) ) THEN
1473	icount = icount + 1
1474	sum = sum + nfld(ni,nk,nj)
1475	END IF
1476	END DO
1477	END DO
1478	END DO
1479	ELSE
1480	sum = 0.
1481	icount = 0
1482	ENDIF
1483	CALL wrf_dm_bcast_real( sum, 1 )
1484	CALL wrf_dm_bcast_integer( icount, 1 )
1485	IF ( enable ) THEN
1486	IF ( icount .NE. 0 ) THEN
1487	avg = sum / REAL ( icount )
1488
1489
1490	! OK, if there were any of those lake situations, we try to search a bit broader
1491	! of an area in the coarse grid.
1492
1493	DO nj = njts, njte
1494	DO nk = nkts, nkte
1495	DO ni = nits, nite
1496	IF ( nfld(ni,nk,nj) .LT. -1.e19 ) THEN
1497	IF ( MOD ( nrj , 2 ) .EQ. 0 ) THEN
1498	cj = ( nj + (nrj/2)-1 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
1499	ELSE
1500	cj = ( nj + (nrj-1)/2 ) / nrj + jpos -1 ! first coarse position equal to or below nest point
1501	END IF
1502	IF ( MOD ( nri , 2 ) .EQ. 0 ) THEN
1503	ci = ( ni + (nri/2)-1 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
1504	ELSE
1505	ci = ( ni + (nri-1)/2 ) / nri + ipos -1 ! first coarse position equal to or to the left of nest point
1506	END IF
1507	ist = MAX (ci-max_search,cits)
1508	ien = MIN (ci+max_search,cite,cide-1)
1509	jst = MAX (cj-max_search,cjts)
1510	jen = MIN (cj+max_search,cjte,cjde-1)
1511	icount = 0
1512	sum = 0
1513	DO jj = jst,jen
1514	DO ii = ist,ien
1515	IF ( NINT(clu(ii,jj)) .EQ. iswater ) THEN
1516	icount = icount + 1
1517	sum = sum + cfld(ii,nk,jj)
1518	END IF
1519	END DO
1520	END DO
1521	IF ( icount .GT. 0 ) THEN
1522	nfld(ni,nk,nj) = sum / REAL ( icount )
1523	ELSE
1524	! CALL wrf_error_fatal ( "horizontal interp error - lake" )
1525	print *,'horizontal interp error - lake, using average ',avg
1526	nfld(ni,nk,nj) = avg
1527	END IF
1528	END IF
1529	END DO
1530	END DO
1531	END DO
1532	ENDIF
1533	ENDIF
1534	ELSE
1535	CALL wrf_error_fatal ( "only unstaggered fields right now" )
1536	END IF
1537
1538	END SUBROUTINE interp_mask_water_field
1539
1540	SUBROUTINE none
1541	END SUBROUTINE none
1542
1543	SUBROUTINE smoother ( cfld , &
1544	cids, cide, ckds, ckde, cjds, cjde, &
1545	cims, cime, ckms, ckme, cjms, cjme, &
1546	cits, cite, ckts, ckte, cjts, cjte, &
1547	nids, nide, nkds, nkde, njds, njde, &
1548	nims, nime, nkms, nkme, njms, njme, &
1549	nits, nite, nkts, nkte, njts, njte, &
1550	xstag, ystag, & ! staggering of field
1551	ipos, jpos, & ! Position of lower left of nest in
1552	nri, nrj &
1553	)
1554
1555	USE module_configure
1556	IMPLICIT NONE
1557
1558	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
1559	cims, cime, ckms, ckme, cjms, cjme, &
1560	cits, cite, ckts, ckte, cjts, cjte, &
1561	nids, nide, nkds, nkde, njds, njde, &
1562	nims, nime, nkms, nkme, njms, njme, &
1563	nits, nite, nkts, nkte, njts, njte, &
1564	nri, nrj, &
1565	ipos, jpos
1566	LOGICAL, INTENT(IN) :: xstag, ystag
1567	INTEGER :: smooth_option, feedback , spec_zone
1568
1569	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
1570
1571	! If there is no feedback, there can be no smoothing.
1572
1573	CALL nl_get_feedback ( 1, feedback )
1574	IF ( feedback == 0 ) RETURN
1575	CALL nl_get_spec_zone ( 1, spec_zone )
1576
1577	! These are the 2d smoothers used on the fedback data. These filters
1578	! are run on the coarse grid data (after the nested info has been
1579	! fedback). Only the area of the nest in the coarse grid is filtered.
1580
1581	CALL nl_get_smooth_option ( 1, smooth_option )
1582
1583	IF ( smooth_option == 0 ) THEN
1584	! no op
1585	ELSE IF ( smooth_option == 1 ) THEN
1586	CALL sm121 ( cfld , &
1587	cids, cide, ckds, ckde, cjds, cjde, &
1588	cims, cime, ckms, ckme, cjms, cjme, &
1589	cits, cite, ckts, ckte, cjts, cjte, &
1590	xstag, ystag, & ! staggering of field
1591	nids, nide, nkds, nkde, njds, njde, &
1592	nims, nime, nkms, nkme, njms, njme, &
1593	nits, nite, nkts, nkte, njts, njte, &
1594	nri, nrj, &
1595	ipos, jpos & ! Position of lower left of nest in
1596	)
1597	ELSE IF ( smooth_option == 2 ) THEN
1598	CALL smdsm ( cfld , &
1599	cids, cide, ckds, ckde, cjds, cjde, &
1600	cims, cime, ckms, ckme, cjms, cjme, &
1601	cits, cite, ckts, ckte, cjts, cjte, &
1602	xstag, ystag, & ! staggering of field
1603	nids, nide, nkds, nkde, njds, njde, &
1604	nims, nime, nkms, nkme, njms, njme, &
1605	nits, nite, nkts, nkte, njts, njte, &
1606	nri, nrj, &
1607	ipos, jpos & ! Position of lower left of nest in
1608	)
1609	END IF
1610
1611	END SUBROUTINE smoother
1612
1613	SUBROUTINE sm121 ( cfld , &
1614	cids, cide, ckds, ckde, cjds, cjde, &
1615	cims, cime, ckms, ckme, cjms, cjme, &
1616	cits, cite, ckts, ckte, cjts, cjte, &
1617	xstag, ystag, & ! staggering of field
1618	nids, nide, nkds, nkde, njds, njde, &
1619	nims, nime, nkms, nkme, njms, njme, &
1620	nits, nite, nkts, nkte, njts, njte, &
1621	nri, nrj, &
1622	ipos, jpos & ! Position of lower left of nest in
1623	)
1624
1625	USE module_configure
1626	IMPLICIT NONE
1627
1628	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
1629	cims, cime, ckms, ckme, cjms, cjme, &
1630	cits, cite, ckts, ckte, cjts, cjte, &
1631	nids, nide, nkds, nkde, njds, njde, &
1632	nims, nime, nkms, nkme, njms, njme, &
1633	nits, nite, nkts, nkte, njts, njte, &
1634	nri, nrj, &
1635	ipos, jpos
1636	LOGICAL, INTENT(IN) :: xstag, ystag
1637
1638	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
1639	REAL, DIMENSION ( cims:cime, cjms:cjme ) :: cfldnew
1640
1641	INTEGER :: i , j , k , loop
1642	INTEGER :: istag,jstag
1643
1644	INTEGER, PARAMETER :: smooth_passes = 1 ! More passes requires a larger stencil (currently 48 pt)
1645
1646	istag = 1 ; jstag = 1
1647	IF ( xstag ) istag = 0
1648	IF ( ystag ) jstag = 0
1649
1650	! Simple 1-2-1 smoother.
1651
1652	smoothing_passes : DO loop = 1 , smooth_passes
1653
1654	DO k = ckts , ckte
1655
1656	! Initialize dummy cfldnew
1657
1658	DO i = MAX(ipos,cits-3) , MIN(ipos+(nide-nids)/nri,cite+3)
1659	DO j = MAX(jpos,cjts-3) , MIN(jpos+(njde-njds)/nrj,cjte+3)
1660	cfldnew(i,j) = cfld(i,k,j)
1661	END DO
1662	END DO
1663
1664	! 1-2-1 smoothing in the j direction first,
1665
1666	DO i = MAX(ipos+2,cits-2) , MIN(ipos+(nide-nids)/nri-2-istag,cite+2)
1667	DO j = MAX(jpos+2,cjts-2) , MIN(jpos+(njde-njds)/nrj-2-jstag,cjte+2)
1668	cfldnew(i,j) = 0.25 * ( cfld(i,k,j+1) + 2.*cfld(i,k,j) + cfld(i,k,j-1) )
1669	END DO
1670	END DO
1671
1672	! then 1-2-1 smoothing in the i direction last
1673
1674	DO j = MAX(jpos+2,cjts-2) , MIN(jpos+(njde-njds)/nrj-2-jstag,cjte+2)
1675	DO i = MAX(ipos+2,cits-2) , MIN(ipos+(nide-nids)/nri-2-istag,cite+2)
1676	cfld(i,k,j) = 0.25 * ( cfldnew(i+1,j) + 2.*cfldnew(i,j) + cfldnew(i-1,j) )
1677	END DO
1678	END DO
1679
1680	END DO
1681
1682	END DO smoothing_passes
1683
1684	END SUBROUTINE sm121
1685
1686	SUBROUTINE smdsm ( cfld , &
1687	cids, cide, ckds, ckde, cjds, cjde, &
1688	cims, cime, ckms, ckme, cjms, cjme, &
1689	cits, cite, ckts, ckte, cjts, cjte, &
1690	xstag, ystag, & ! staggering of field
1691	nids, nide, nkds, nkde, njds, njde, &
1692	nims, nime, nkms, nkme, njms, njme, &
1693	nits, nite, nkts, nkte, njts, njte, &
1694	nri, nrj, &
1695	ipos, jpos & ! Position of lower left of nest in
1696	)
1697
1698	USE module_configure
1699	IMPLICIT NONE
1700
1701	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
1702	cims, cime, ckms, ckme, cjms, cjme, &
1703	cits, cite, ckts, ckte, cjts, cjte, &
1704	nids, nide, nkds, nkde, njds, njde, &
1705	nims, nime, nkms, nkme, njms, njme, &
1706	nits, nite, nkts, nkte, njts, njte, &
1707	nri, nrj, &
1708	ipos, jpos
1709	LOGICAL, INTENT(IN) :: xstag, ystag
1710
1711	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
1712	REAL, DIMENSION ( cims:cime, cjms:cjme ) :: cfldnew
1713
1714	REAL , DIMENSION ( 2 ) :: xnu
1715	INTEGER :: i , j , k , loop , n
1716	INTEGER :: istag,jstag
1717
1718	INTEGER, PARAMETER :: smooth_passes = 1 ! More passes requires a larger stencil (currently 48 pt)
1719
1720	xnu = (/ 0.50 , -0.52 /)
1721
1722	istag = 1 ; jstag = 1
1723	IF ( xstag ) istag = 0
1724	IF ( ystag ) jstag = 0
1725
1726	! The odd number passes of this are the "smoother", the even
1727	! number passes are the "de-smoother" (note the different signs on xnu).
1728
1729	smoothing_passes : DO loop = 1 , smooth_passes * 2
1730
1731	n = 2 - MOD ( loop , 2 )
1732
1733	DO k = ckts , ckte
1734
1735	DO i = MAX(ipos+2,cits-2) , MIN(ipos+(nide-nids)/nri-2-istag,cite+2)
1736	DO j = MAX(jpos+2,cjts-2) , MIN(jpos+(njde-njds)/nrj-2-jstag,cjte+2)
1737	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))
1738	END DO
1739	END DO
1740
1741	DO i = MAX(ipos+2,cits-2) , MIN(ipos+(nide-nids)/nri-2-istag,cite+2)
1742	DO j = MAX(jpos+2,cjts-2) , MIN(jpos+(njde-njds)/nrj-2-jstag,cjte+2)
1743	cfld(i,k,j) = cfldnew(i,j)
1744	END DO
1745	END DO
1746
1747	DO j = MAX(jpos+2,cjts-2) , MIN(jpos+(njde-njds)/nrj-2-jstag,cjte+2)
1748	DO i = MAX(ipos+2,cits-2) , MIN(ipos+(nide-nids)/nri-2-istag,cite+2)
1749	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))
1750	END DO
1751	END DO
1752
1753	DO j = MAX(jpos+2,cjts-2) , MIN(jpos+(njde-njds)/nrj-2-jstag,cjte+2)
1754	DO i = MAX(ipos+2,cits-2) , MIN(ipos+(nide-nids)/nri-2-istag,cite+2)
1755	cfld(i,k,j) = cfldnew(i,j)
1756	END DO
1757	END DO
1758
1759	END DO
1760
1761	END DO smoothing_passes
1762
1763	END SUBROUTINE smdsm
1764
1765	!==================================
1766	! this is used to modify a field over the nest so we can see where the nest is
1767
1768	SUBROUTINE mark_domain ( cfld, & ! CD field
1769	cids, cide, ckds, ckde, cjds, cjde, &
1770	cims, cime, ckms, ckme, cjms, cjme, &
1771	cits, cite, ckts, ckte, cjts, cjte, &
1772	nfld, & ! ND field
1773	nids, nide, nkds, nkde, njds, njde, &
1774	nims, nime, nkms, nkme, njms, njme, &
1775	nits, nite, nkts, nkte, njts, njte, &
1776	shw, & ! stencil half width for interp
1777	imask, & ! interpolation mask
1778	xstag, ystag, & ! staggering of field
1779	ipos, jpos, & ! Position of lower left of nest in CD
1780	nri, nrj ) ! nest ratios
1781	USE module_configure
1782	USE module_wrf_error
1783	IMPLICIT NONE
1784
1785
1786	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
1787	cims, cime, ckms, ckme, cjms, cjme, &
1788	cits, cite, ckts, ckte, cjts, cjte, &
1789	nids, nide, nkds, nkde, njds, njde, &
1790	nims, nime, nkms, nkme, njms, njme, &
1791	nits, nite, nkts, nkte, njts, njte, &
1792	shw, &
1793	ipos, jpos, &
1794	nri, nrj
1795	LOGICAL, INTENT(IN) :: xstag, ystag
1796
1797	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ), INTENT(OUT) :: cfld
1798	REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ), INTENT(IN) :: nfld
1799	INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: imask
1800
1801	! Local
1802
1803	INTEGER ci, cj, ck, ni, nj, nk, ip, jp, ioff, joff, ioffa, joffa
1804	INTEGER :: icmin,icmax,jcmin,jcmax
1805	INTEGER :: istag,jstag, ipoints,jpoints,ijpoints
1806
1807	istag = 1 ; jstag = 1
1808	IF ( xstag ) istag = 0
1809	IF ( ystag ) jstag = 0
1810
1811	DO cj = MAX(jpos+1,cjts),MIN(jpos+(njde-njds)/nrj-jstag-1,cjte)
1812	nj = (cj-jpos)*nrj + jstag + 1
1813	DO ck = ckts, ckte
1814	nk = ck
1815	DO ci = MAX(ipos+1,cits),MIN(ipos+(nide-nids)/nri-istag-1,cite)
1816	ni = (ci-ipos)*nri + istag + 1
1817	cfld( ci, ck, cj ) = 9021000. !magic number: Beverly Hills * 100.
1818	ENDDO
1819	ENDDO
1820	ENDDO
1821
1822	END SUBROUTINE mark_domain
1823
1824	#if ( NMM_CORE == 1 )
1825	!=======================================================================================
1826	! E grid interpolation for mass with addition of terrain adjustments. First routine
1827	! pertains to initial conditions and the next one corresponds to boundary conditions
1828	! This is gopal's doing
1829	!=======================================================================================
1830
1831	SUBROUTINE interp_mass_nmm (cfld, & ! CD field
1832	cids, cide, ckds, ckde, cjds, cjde, &
1833	cims, cime, ckms, ckme, cjms, cjme, &
1834	cits, cite, ckts, ckte, cjts, cjte, &
1835	nfld, & ! ND field
1836	nids, nide, nkds, nkde, njds, njde, &
1837	nims, nime, nkms, nkme, njms, njme, &
1838	nits, nite, nkts, nkte, njts, njte, &
1839	shw, & ! stencil half width for interp
1840	imask, & ! interpolation mask
1841	xstag, ystag, & ! staggering of field
1842	ipos, jpos, & ! Position of lower left of nest in CD
1843	nri, nrj, & ! nest ratios
1844	CII, IIH, CJJ, JJH, CBWGT1, HBWGT1, & ! south-western grid locs and weights
1845	CBWGT2, HBWGT2, CBWGT3, HBWGT3, & ! note that "C"ourse grid ones are
1846	CBWGT4, HBWGT4, & ! dummys for weights
1847	CZ3d, Z3d, & ! Z3d interpolated from CZ3d
1848	CFIS,FIS, & ! CFIS dummy on fine domain
1849	CSM,SM, & ! CSM is dummy
1850	CPDTOP,PDTOP, &
1851	CPTOP,PTOP, &
1852	CPSTD,PSTD, &
1853	CKZMAX,KZMAX )
1854
1855	USE MODULE_MODEL_CONSTANTS
1856	USE module_timing
1857	IMPLICIT NONE
1858
1859	LOGICAL,INTENT(IN) :: xstag, ystag
1860	INTEGER,INTENT(IN) :: ckzmax,kzmax
1861	INTEGER,INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
1862	cims, cime, ckms, ckme, cjms, cjme, &
1863	cits, cite, ckts, ckte, cjts, cjte, &
1864	nids, nide, nkds, nkde, njds, njde, &
1865	nims, nime, nkms, nkme, njms, njme, &
1866	nits, nite, nkts, nkte, njts, njte, &
1867	shw,ipos,jpos,nri,nrj
1868
1869	INTEGER,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: IMASK
1870
1871	! parent domain
1872
1873	INTEGER,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CII,CJJ ! dummy
1874	REAL,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CBWGT1,CBWGT2,CBWGT3
1875	REAL,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CBWGT4,CFIS,CSM
1876	REAL,DIMENSION(cims:cime,cjms:cjme,ckms:ckme), INTENT(IN) :: CFLD
1877	REAL,DIMENSION(cims:cime,cjms:cjme,1:KZMAX), INTENT(IN) :: CZ3d
1878	REAL,DIMENSION(1:KZMAX), INTENT(IN) :: CPSTD
1879	REAL,INTENT(IN) :: CPDTOP,CPTOP
1880
1881	! nested domain
1882
1883	INTEGER,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: IIH,JJH
1884	REAL,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: HBWGT1,HBWGT2,HBWGT3
1885	REAL,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: HBWGT4
1886	REAL,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: FIS,SM
1887	REAL,DIMENSION(nims:nime,njms:njme,nkms:nkme), INTENT(INOUT) :: NFLD
1888	REAL,DIMENSION(1:KZMAX), INTENT(IN) :: PSTD
1889	REAL,DIMENSION(nims:nime,njms:njme,1:KZMAX), INTENT(OUT) :: Z3d
1890	REAL,INTENT(IN) :: PDTOP,PTOP
1891
1892	! local
1893
1894	INTEGER,PARAMETER :: JTB=134
1895	REAL, PARAMETER :: LAPSR=6.5E-3,GI=1./G, D608=0.608
1896	REAL, PARAMETER :: COEF3=R_DGILAPSR
1897	INTEGER :: I,J,K,IDUM
1898	REAL :: dlnpdz,tvout,pmo
1899	REAL,DIMENSION(nims:nime,njms:njme) :: ZS,DUM2d
1900	REAL,DIMENSION(JTB) :: PIN,ZIN,Y2,PIO,ZOUT,DUM1,DUM2
1901	!-----------------------------------------------------------------------------------------------------
1902	!
1903	!*** CHECK DOMAIN BOUNDS BEFORE INTERPOLATION
1904	!
1905	DO J=NJTS,MIN(NJTE,NJDE-1)
1906	DO I=NITS,MIN(NITE,NIDE-1)
1907	IF(IIH(i,j).LT.(CIDS-shw) .OR. IIH(i,j).GT.(CIDE+shw)) &
1908	CALL wrf_error_fatal ('mass points:check domain bounds along x' )
1909	IF(JJH(i,j).LT.(CJDS-shw) .OR. JJH(i,j).GT.(CJDE+shw)) &
1910	CALL wrf_error_fatal ('mass points:check domain bounds along y' )
1911	ENDDO
1912	ENDDO
1913
1914	IF(KZMAX .GT. (JTB-10)) &
1915	CALL wrf_error_fatal ('mass points: increase JTB in interp_mass_nmm')
1916
1917	! WRITE(21,*)'------------- MED NEST INITIAL 1 ----------------'
1918	! DO J=NJTS,MIN(NJTE,NJDE-1)
1919	! DO I=NITS,MIN(NITE,NIDE-1)
1920	! WRITE(21,*)I,J,IMASK(I,J),NFLD(I,1,J)
1921	! ENDDO
1922	! ENDDO
1923	! WRITE(21,*)
1924
1925	!
1926	!*** DEFINE LOCAL TOPOGRAPHY ON THE BASIS OF FIS. ALSO CHECK IF SM IS LAND (SM=0) OVER TOPO
1927	!*** YOU DON'T WANT MOUNTAINS INSIDE WATER BODIES!
1928	!
1929
1930	DO J=NJTS,MIN(NJTE,NJDE-1)
1931	DO I=NITS,MIN(NITE,NIDE-1)
1932	ZS(I,J)=FIS(I,J)/G
1933	ENDDO
1934	ENDDO
1935
1936	!
1937	!*** Interpolate GPMs DERIVED FROM STANDARD ATMOSPHERIC LAPSE RATE FROM THE PARENT TO
1938	!*** THE NESTED DOMAIN
1939	!
1940	!*** INDEX CONVENTIONS
1941	!*** HBWGT4
1942	!*** 4
1943	!***
1944	!***
1945	!***
1946	!*** h
1947	!*** 1 2
1948	!*** HBWGT1 HBWGT2
1949	!***
1950	!***
1951	!*** 3
1952	!*** HBWGT3
1953
1954	Z3d=0.0
1955	DO K=NKTS,KZMAX ! Please note that we are still in isobaric surfaces
1956	DO J=NJTS,MIN(NJTE,NJDE-1)
1957	DO I=NITS,MIN(NITE,NIDE-1)
1958	!
1959	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7
1960	Z3d(I,J,K) = HBWGT1(I,J)*CZ3d(IIH(I,J), JJH(I,J) ,K) &
1961	+ HBWGT2(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J) ,K) &
1962	+ HBWGT3(I,J)*CZ3d(IIH(I,J), JJH(I,J)-1,K) &
1963	+ HBWGT4(I,J)*CZ3d(IIH(I,J), JJH(I,J)+1,K)
1964	ELSE
1965	Z3d(I,J,K) = HBWGT1(I,J)*CZ3d(IIH(I,J), JJH(I,J) ,K) &
1966	+ HBWGT2(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J) ,K) &
1967	+ HBWGT3(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J)-1,K) &
1968	+ HBWGT4(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J)+1,K)
1969
1970	ENDIF
1971	!
1972	ENDDO
1973	ENDDO
1974	ENDDO
1975
1976	! RECONSTRUCT PDs ON THE BASIS OF TOPOGRAPHY AND THE INTERPOLATED HEIGHTS
1977
1978	DO J=NJTS,MIN(NJTE,NJDE-1)
1979	DO I=NITS,MIN(NITE,NIDE-1)
1980	!
1981	IF (ZS(I,J) .LT. Z3d(I,J,1)) THEN
1982	dlnpdz = (log(PSTD(1))-log(PSTD(2)) )/(Z3d(i,j,1)-Z3d(i,j,2))
1983	dum2d(i,j) = exp(log(PSTD(1)) + dlnpdz*(ZS(I,J) - Z3d(i,j,1)))
1984	dum2d(i,j) = dum2d(i,j) - PDTOP -PTOP
1985	ELSE ! target level bounded by input levels
1986	DO K =NKTS,KZMAX-1 ! still in the isobaric surfaces
1987	IF(ZS(I,J) .GE. Z3d(I,J,K) .AND. ZS(I,J) .LT. Z3d(I,J,K+1))THEN
1988	dlnpdz = (log(PSTD(K))-log(PSTD(K+1)) ) /(Z3d(I,J,K)-Z3d(I,J,K+1))
1989	dum2d(i,j) = exp(log(PSTD(K)) + dlnpdz*(ZS(I,J)- Z3d(I,J,K)))
1990	dum2d(i,j) = dum2d(i,j) - PDTOP -PTOP
1991	ENDIF
1992	ENDDO
1993	ENDIF
1994	IF(ZS(I,J) .GE. Z3d(I,J,KZMAX))THEN
1995	WRITE(0,*)'I=',I,'J=',J,'K=',KZMAX,'TERRAIN HEIGHT',ZS(I,J),'Z3d',Z3d(I,J,KZMAX)
1996	CALL wrf_error_fatal3 ( "interp_fcn.b" , 176 , "MOUNTAIN TOO HIGH TO FIT THE MODEL DEPTH")
1997	ENDIF
1998	!
1999	ENDDO
2000	ENDDO
2001
2002	DO K=NKDS,NKDE ! NKTE is 1, nevertheless let us pretend religious
2003	DO J=NJTS,MIN(NJTE,NJDE-1)
2004	DO I=NITS,MIN(NITE,NIDE-1)
2005	IF(IMASK(I,J) .NE. 1)THEN
2006	NFLD(I,J,K)= dum2d(i,j) ! PD defined in the nested domain
2007	ENDIF
2008	ENDDO
2009	ENDDO
2010	ENDDO
2011
2012	!
2013	END SUBROUTINE interp_mass_nmm
2014	!
2015	!--------------------------------------------------------------------------------------
2016
2017	SUBROUTINE nmm_bdymass_hinterp ( cfld, & ! CD field
2018	cids, cide, ckds, ckde, cjds, cjde, &
2019	cims, cime, ckms, ckme, cjms, cjme, &
2020	cits, cite, ckts, ckte, cjts, cjte, &
2021	nfld, & ! ND field
2022	nids, nide, nkds, nkde, njds, njde, &
2023	nims, nime, nkms, nkme, njms, njme, &
2024	nits, nite, nkts, nkte, njts, njte, &
2025	shw, & ! stencil half width
2026	imask, & ! interpolation mask
2027	xstag, ystag, & ! staggering of field
2028	ipos, jpos, & ! Position of lower left of nest in CD
2029	nri, nrj, & ! nest ratios
2030	c_bxs,n_bxs, &
2031	c_bxe,n_bxe, &
2032	c_bys,n_bys, &
2033	c_bye,n_bye, &
2034	c_btxs,n_btxs, &
2035	c_btxe,n_btxe, &
2036	c_btys,n_btys, &
2037	c_btye,n_btye, &
2038	CTEMP_B,NTEMP_B, & ! These temp arrays should be removed
2039	CTEMP_BT,NTEMP_BT, & ! later on
2040	CII, IIH, CJJ, JJH, CBWGT1, HBWGT1, & ! south-western grid locs and weights
2041	CBWGT2, HBWGT2, CBWGT3, HBWGT3, & ! note that "C"ourse grid ones are
2042	CBWGT4, HBWGT4, & ! dummys
2043	CZ3d, Z3d, & ! Z3d dummy on nested domain
2044	CFIS,FIS, & ! CFIS dummy on fine domain
2045	CSM,SM, & ! CSM is dummy
2046	CPDTOP,PDTOP, &
2047	CPTOP,PTOP, &
2048	CPSTD,PSTD, &
2049	CKZMAX,KZMAX )
2050
2051
2052	USE MODULE_MODEL_CONSTANTS
2053	USE module_configure
2054	USE module_wrf_error
2055
2056	IMPLICIT NONE
2057
2058
2059	INTEGER, INTENT(IN) :: ckzmax,kzmax
2060	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
2061	cims, cime, ckms, ckme, cjms, cjme, &
2062	cits, cite, ckts, ckte, cjts, cjte, &
2063	nids, nide, nkds, nkde, njds, njde, &
2064	nims, nime, nkms, nkme, njms, njme, &
2065	nits, nite, nkts, nkte, njts, njte, &
2066	shw, &
2067	ipos, jpos, &
2068	nri, nrj
2069
2070
2071	REAL, DIMENSION ( cims:cime, cjms:cjme, ckms:ckme ), INTENT(OUT) :: ctemp_b,ctemp_bt
2072	REAL, DIMENSION ( nims:nime, njms:njme, nkms:nkme ), INTENT(OUT) :: ntemp_b,ntemp_bt
2073	LOGICAL, INTENT(IN) :: xstag, ystag
2074	REAL, DIMENSION( * ), INTENT(INOUT) :: c_bxs,n_bxs,c_bxe,n_bxe,c_bys,n_bys,c_bye,n_bye
2075	REAL, DIMENSION( * ), INTENT(INOUT) :: c_btxs,n_btxs,c_btxe,n_btxe,c_btys,n_btys,c_btye,n_btye
2076
2077	! parent domain
2078
2079	INTEGER,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: IMASK
2080	INTEGER,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CII,CJJ ! dummy
2081	REAL,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CBWGT1,CBWGT2,CBWGT3
2082	REAL,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CBWGT4,CFIS,CSM
2083	REAL,DIMENSION(cims:cime,cjms:cjme,ckms:ckme), INTENT(IN) :: CFLD
2084	REAL,DIMENSION(cims:cime,cjms:cjme,1:KZMAX), INTENT(IN) :: CZ3d
2085	REAL,DIMENSION(1:KZMAX), INTENT(IN) :: CPSTD
2086	REAL,INTENT(IN) :: CPDTOP,CPTOP
2087
2088	! nested domain
2089
2090	INTEGER,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: IIH,JJH
2091	REAL,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: HBWGT1,HBWGT2,HBWGT3
2092	REAL,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: HBWGT4
2093	REAL,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: FIS,SM
2094	REAL,DIMENSION(nims:nime,njms:njme,nkms:nkme), INTENT(INOUT) :: NFLD
2095	REAL,DIMENSION(1:KZMAX), INTENT(IN) :: PSTD
2096	REAL,DIMENSION(nims:nime,njms:njme,1:KZMAX), INTENT(OUT) :: Z3d
2097	REAL,INTENT(IN) :: PDTOP,PTOP
2098
2099	! Local
2100
2101	INTEGER :: nijds, nijde, spec_bdy_width,i,j,k
2102	REAL :: dlnpdz,dum2d
2103	REAL,DIMENSION(nims:nime,njms:njme) :: zs
2104
2105	INTEGER,PARAMETER :: JTB=134
2106	INTEGER :: ii,jj
2107	REAL, DIMENSION (nims:nime,njms:njme) :: CWK1,CWK2,CWK3,CWK4
2108
2109	nijds = min(nids, njds)
2110	nijde = max(nide, njde)
2111	CALL nl_get_spec_bdy_width( 1, spec_bdy_width )
2112
2113	!
2114	!*** DEFINE LOCAL TOPOGRAPHY ON THE BASIS OF FIS. ASLO CHECK IF SM IS LAND (SM=0) OVER TOPO
2115	!*** YOU DON'T WANT MOUNTAINS INSIDE WATER BODIES!
2116	!
2117	DO J=NJTS,MIN(NJTE,NJDE-1)
2118	DO I=NITS,MIN(NITE,NIDE-1)
2119	ZS(I,J)=FIS(I,J)/G
2120	ENDDO
2121	ENDDO
2122
2123	! X start boundary
2124
2125	NMM_XS: IF(NITS .EQ. NIDS)THEN
2126	! WRITE(0,*)'ENTERING X1 START BOUNDARY AT MASS POINTS',NJTS,MIN(NJTE,NJDE-1)
2127	I = NIDS
2128
2129	DO K=NKTS,KZMAX
2130	DO J = NJTS,MIN(NJTE,NJDE-1)
2131	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7 of the parent domain
2132	Z3d(I,J,K) = HBWGT1(I,J)*CZ3d(IIH(I,J), JJH(I,J) ,K) &
2133	+ HBWGT2(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J) ,K) &
2134	+ HBWGT3(I,J)*CZ3d(IIH(I,J), JJH(I,J)-1,K) &
2135	+ HBWGT4(I,J)*CZ3d(IIH(I,J), JJH(I,J)+1,K)
2136	ELSE
2137	Z3d(I,J,K) = HBWGT1(I,J)*CZ3d(IIH(I,J), JJH(I,J) ,K) &
2138	+ HBWGT2(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J) ,K) &
2139	+ HBWGT3(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J)-1,K) &
2140	+ HBWGT4(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J)+1,K)
2141	ENDIF
2142	END DO
2143	END DO
2144
2145	DO J = NJTS,MIN(NJTE,NJDE-1)
2146	IF(MOD(J,2) .NE. 0)THEN
2147	IF (ZS(I,J) .LT. Z3d(I,J,2)) THEN ! level 2 has to be changed
2148	dlnpdz = (log(PSTD(1))-log(PSTD(2)) )/(Z3d(I,J,1)-Z3d(I,J,2))
2149	dum2d = exp(log(PSTD(1)) + dlnpdz*(ZS(I,J) - Z3d(I,J,1)))
2150	CWK1(I,J) = dum2d -PDTOP -PTOP
2151	ELSE ! target level bounded by input levels
2152	DO K =NKTS,KZMAX-1
2153	IF(ZS(I,J) .GE. Z3d(I,J,K) .AND. ZS(I,J) .LT. Z3d(I,J,K+1))THEN
2154	dlnpdz = (log(PSTD(K))-log(PSTD(K+1)) ) /(Z3d(I,J,K)-Z3d(I,J,K+1))
2155	dum2d = exp(log(PSTD(K)) + dlnpdz*(ZS(I,J)- Z3d(I,J,K)))
2156	CWK1(I,J) = dum2d -PDTOP -PTOP
2157	ENDIF
2158	ENDDO
2159	ENDIF
2160	IF(ZS(I,J) .GE. Z3d(I,J,KZMAX))THEN
2161	WRITE(0,*)'I=',I,'J=',J,'K=',K,'TERRAIN HEIGHT',ZS(I,J),'Z3d',Z3d(I,J,KZMAX)
2162	CALL wrf_error_fatal("BC:MOUNTAIN TOO HIGH TO FIT THE MODEL DEPTH")
2163	ENDIF
2164	ELSE
2165	CWK1(I,J)=0.
2166	ENDIF
2167	ENDDO
2168
2169	DO J = NJTS,MIN(NJTE,NJDE-1)
2170	DO K = NKDS,NKDE
2171	ntemp_b(i,j,k) = CWK1(I,J)
2172	ntemp_bt(i,j,k) = 0.0
2173	END DO
2174	END DO
2175	ENDIF NMM_XS
2176
2177	! X end boundary
2178
2179	NMM_XE: IF(NITE-1 .EQ. NIDE-1)THEN
2180	! WRITE(0,*)'ENTERING X END BOUNDARY AT MASS POINTS',NJTS,MIN(NJTE,NJDE-1)
2181	I = NIDE-1
2182	II = NIDE - I
2183
2184	DO K=NKTS,KZMAX
2185	DO J=NJTS,MIN(NJTE,NJDE-1)
2186	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7
2187	Z3d(I,J,K) = HBWGT1(I,J)*CZ3d(IIH(I,J), JJH(I,J) ,K) &
2188	+ HBWGT2(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J) ,K) &
2189	+ HBWGT3(I,J)*CZ3d(IIH(I,J), JJH(I,J)-1,K) &
2190	+ HBWGT4(I,J)*CZ3d(IIH(I,J), JJH(I,J)+1,K)
2191	ELSE
2192	Z3d(I,J,K) = HBWGT1(I,J)*CZ3d(IIH(I,J), JJH(I,J) ,K) &
2193	+ HBWGT2(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J) ,K) &
2194	+ HBWGT3(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J)-1,K) &
2195	+ HBWGT4(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J)+1,K)
2196	ENDIF
2197	ENDDO
2198	ENDDO
2199
2200	DO J = NJTS,MIN(NJTE,NJDE-1)
2201	IF(MOD(J,2) .NE.0)THEN ! 1,3,5,7 of nested domain
2202	IF (ZS(I,J) .LT. Z3d(I,J,2)) THEN ! level 2 has to be changed
2203	dlnpdz = (log(PSTD(1))-log(PSTD(2)) )/(Z3d(I,J,1)-Z3d(I,J,2))
2204	dum2d = exp(log(PSTD(1)) + dlnpdz*(ZS(I,J) - Z3d(I,J,1)))
2205	CWK2(I,J) = dum2d -PDTOP -PTOP
2206	ELSE ! target level bounded by input levels
2207	DO K =NKTS,KZMAX-1
2208	IF(ZS(I,J) .GE. Z3d(I,J,K) .AND. ZS(I,J) .LT. Z3d(I,J,K+1))THEN
2209	dlnpdz = (log(PSTD(K))-log(PSTD(K+1)) ) /(Z3d(I,J,K)-Z3d(I,J,K+1))
2210	dum2d = exp(log(PSTD(K)) + dlnpdz*(ZS(I,J)- Z3d(I,J,K)))
2211	CWK2(I,J) = dum2d -PDTOP -PTOP
2212	ENDIF
2213	ENDDO
2214	ENDIF
2215	IF(ZS(I,J) .GE. Z3d(I,J,KZMAX))THEN
2216	WRITE(0,*)'I=',I,'J=',J,'K=',K,'TERRAIN HEIGHT',ZS(I,J),'Z3d',Z3d(I,J,KZMAX)
2217	CALL wrf_error_fatal("BC:MOUNTAIN TOO HIGH TO FIT THE MODEL DEPTH")
2218	ENDIF
2219	ELSE
2220	CWK2(I,J) = 0.0
2221	ENDIF
2222	ENDDO
2223
2224	DO J = NJTS,MIN(NJTE,NJDE-1)
2225	DO K = NKDS,NKDE
2226	ntemp_b(i,j,k) = CWK2(I,J)
2227	ntemp_bt(i,j,k) = 0.0
2228	END DO
2229	END DO
2230	ENDIF NMM_XE
2231
2232	! Y start boundary
2233
2234	NMM_YS: IF(NJTS .EQ. NJDS)THEN
2235	! WRITE(20,*)'ENTERING Y START BOUNDARY AT MASS POINTS',NITS,MIN(NITE,NIDE-1)
2236	J = NJDS
2237	DO K=NKTS,KZMAX
2238	DO I = NITS,MIN(NITE,NIDE-1)
2239	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7
2240	Z3d(I,J,K) = HBWGT1(I,J)*CZ3d(IIH(I,J), JJH(I,J) ,K) &
2241	+ HBWGT2(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J) ,K) &
2242	+ HBWGT3(I,J)*CZ3d(IIH(I,J), JJH(I,J)-1,K) &
2243	+ HBWGT4(I,J)*CZ3d(IIH(I,J), JJH(I,J)+1,K)
2244	ELSE
2245	Z3d(I,J,K) = HBWGT1(I,J)*CZ3d(IIH(I,J), JJH(I,J) ,K) &
2246	+ HBWGT2(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J) ,K) &
2247	+ HBWGT3(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J)-1,K) &
2248	+ HBWGT4(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J)+1,K)
2249	ENDIF
2250	END DO
2251	END DO
2252
2253	DO I = NITS,MIN(NITE,NIDE-1)
2254	IF (ZS(I,J) .LT. Z3d(I,J,2)) THEN ! level 2 has to be changed
2255	dlnpdz = (log(PSTD(1))-log(PSTD(2)) )/(Z3d(I,J,1)-Z3d(I,J,2))
2256	dum2d = exp(log(PSTD(1)) + dlnpdz*(ZS(I,J) - Z3d(I,J,1)))
2257	CWK3(I,J) = dum2d -PDTOP -PTOP
2258	ELSE ! target level bounded by input levels
2259	DO K =NKTS,KZMAX-1
2260	IF(ZS(I,J) .GE. Z3d(I,J,K) .AND. ZS(I,J) .LT. Z3d(I,J,K+1))THEN
2261	dlnpdz = (log(PSTD(K))-log(PSTD(K+1)) ) /(Z3d(I,J,K)-Z3d(I,J,K+1))
2262	dum2d = exp(log(PSTD(K)) + dlnpdz*(ZS(I,J)- Z3d(I,J,K)))
2263	CWK3(I,J) = dum2d -PDTOP -PTOP
2264	ENDIF
2265	ENDDO
2266	ENDIF
2267	IF(ZS(I,J) .GE. Z3d(I,J,KZMAX))THEN
2268	WRITE(0,*)'I=',I,'J=',J,'K=',K,'TERRAIN HEIGHT',ZS(I,J),'Z3d',Z3d(I,J,KZMAX)
2269	CALL wrf_error_fatal("BC:MOUNTAIN TOO HIGH TO FIT THE MODEL DEPTH")
2270	ENDIF
2271	ENDDO
2272
2273	DO K = NKDS, NKDE
2274	DO I = NITS,MIN(NITE,NIDE-1)
2275	ntemp_b(i,j,k) = CWK3(I,J)
2276	ntemp_bt(i,j,k) = 0.0
2277	END DO
2278	END DO
2279	END IF NMM_YS
2280
2281	! Y end boundary
2282
2283	NMM_YE: IF(NJTE-1 .EQ. NJDE-1)THEN
2284	! WRITE(20,*)'ENTERING Y END BOUNDARY AT MASS POINTS',NITS,MIN(NITE,NIDE-1)
2285	J = NJDE-1
2286	JJ = NJDE - J
2287	DO K=NKTS,KZMAX
2288	DO I = NITS,MIN(NITE,NIDE-1)
2289	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7
2290	Z3d(I,J,K) = HBWGT1(I,J)*CZ3d(IIH(I,J), JJH(I,J) ,K) &
2291	+ HBWGT2(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J) ,K) &
2292	+ HBWGT3(I,J)*CZ3d(IIH(I,J), JJH(I,J)-1,K) &
2293	+ HBWGT4(I,J)*CZ3d(IIH(I,J), JJH(I,J)+1,K)
2294	ELSE
2295	Z3d(I,J,K) = HBWGT1(I,J)*CZ3d(IIH(I,J), JJH(I,J) ,K) &
2296	+ HBWGT2(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J) ,K) &
2297	+ HBWGT3(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J)-1,K) &
2298	+ HBWGT4(I,J)*CZ3d(IIH(I,J)+1,JJH(I,J)+1,K)
2299	ENDIF
2300	END DO
2301	END DO
2302
2303	DO I = NITS,MIN(NITE,NIDE-1)
2304	IF (ZS(I,J) .LT. Z3d(I,J,2)) THEN ! level 2 has to be changed
2305	dlnpdz = (log(PSTD(1))-log(PSTD(2)) )/(Z3d(I,J,1)-Z3d(I,J,2))
2306	dum2d = exp(log(PSTD(1)) + dlnpdz*(ZS(I,J) - Z3d(I,J,1)))
2307	CWK4(I,J) = dum2d -PDTOP -PTOP
2308	ELSE ! target level bounded by input levels
2309	DO K =NKTS,KZMAX-1
2310	IF(ZS(I,J) .GE. Z3d(I,J,K) .AND. ZS(I,J) .LT. Z3d(I,J,K+1))THEN
2311	dlnpdz = (log(PSTD(K))-log(PSTD(K+1)) ) /(Z3d(I,J,K)-Z3d(I,J,K+1))
2312	dum2d = exp(log(PSTD(K)) + dlnpdz*(ZS(I,J)- Z3d(I,J,K)))
2313	CWK4(I,J) = dum2d -PDTOP -PTOP
2314	ENDIF
2315	ENDDO
2316	ENDIF
2317	IF(ZS(I,J) .GE. Z3d(I,J,KZMAX))THEN
2318	WRITE(0,*)'I=',I,'J=',J,'K=',K,'TERRAIN HEIGHT',ZS(I,J),'Z3d',Z3d(I,J,KZMAX)
2319	CALL wrf_error_fatal("BC:MOUNTAIN TOO HIGH TO FIT THE MODEL DEPTH")
2320	ENDIF
2321	ENDDO
2322
2323	DO K = NKDS,NKDE
2324	DO I = NITS,MIN(NITE,NIDE-1)
2325	ntemp_b(i,j,k) = CWK4(I,J)
2326	ntemp_bt(i,j,k) = 0.0
2327	END DO
2328	END DO
2329	END IF NMM_YE
2330
2331	RETURN
2332
2333	END SUBROUTINE nmm_bdymass_hinterp
2334	!
2335	!=======================================================================================
2336	!
2337	! ADDED FOR INCLUDING MOISTURE AND THERMODYNAMIC ENERGY BALANCE
2338	!
2339	!=======================================================================================
2340
2341	SUBROUTINE interp_scalar_nmm (cfld, & ! CD field
2342	cids,cide,ckds,ckde,cjds,cjde, &
2343	cims,cime,ckms,ckme,cjms,cjme, &
2344	cits,cite,ckts,ckte,cjts,cjte, &
2345	nfld, & ! ND field
2346	nids,nide,nkds,nkde,njds,njde, &
2347	nims,nime,nkms,nkme,njms,njme, &
2348	nits,nite,nkts,nkte,njts,njte, &
2349	shw, & ! stencil half width for interp
2350	imask, & ! interpolation mask
2351	xstag,ystag, & ! staggering of field
2352	ipos,jpos, & ! Position of lower left of nest in CD
2353	nri,nrj, & ! nest ratios
2354	CII, IIH, CJJ, JJH, CBWGT1, HBWGT1, & ! south-western grid locs and weights
2355	CBWGT2, HBWGT2, CBWGT3, HBWGT3, & ! note that "C"ourse grid ones are
2356	CBWGT4, HBWGT4, & ! dummys for weights
2357	CC3d,C3d, &
2358	CPD,PD, &
2359	CPSTD,PSTD, &
2360	CPDTOP,PDTOP, &
2361	CPTOP,PTOP, &
2362	CETA1,ETA1,CETA2,ETA2 )
2363
2364	USE MODULE_MODEL_CONSTANTS
2365	USE module_timing
2366	IMPLICIT NONE
2367
2368	LOGICAL,INTENT(IN) :: xstag, ystag
2369	INTEGER,INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
2370	cims, cime, ckms, ckme, cjms, cjme, &
2371	cits, cite, ckts, ckte, cjts, cjte, &
2372	nids, nide, nkds, nkde, njds, njde, &
2373	nims, nime, nkms, nkme, njms, njme, &
2374	nits, nite, nkts, nkte, njts, njte, &
2375	shw,ipos,jpos,nri,nrj
2376
2377	INTEGER,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: IMASK
2378
2379	! parent domain
2380
2381	INTEGER,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CII,CJJ ! dummy
2382	REAL,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CBWGT1,CBWGT2
2383	REAL,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CBWGT3,CBWGT4
2384
2385	REAL,DIMENSION(cims:cime,cjms:cjme,ckms:ckme), INTENT(IN) :: CFLD
2386	REAL,DIMENSION(cims:cime,cjms:cjme,ckms:ckme), INTENT(IN) :: CC3d ! scalar input on constant pressure levels
2387	REAL,DIMENSION(ckms:ckme), INTENT(IN) :: CPSTD
2388	REAL,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CPD
2389	REAL,DIMENSION(ckms:ckme), INTENT(IN) :: CETA1,CETA2
2390	REAL, INTENT(IN) :: CPDTOP,CPTOP
2391
2392	! nested domain
2393
2394	INTEGER,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: IIH,JJH
2395	REAL,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: HBWGT1,HBWGT2
2396	REAL,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: HBWGT3,HBWGT4
2397
2398	REAL,DIMENSION(nims:nime,njms:njme,nkms:nkme), INTENT(OUT):: NFLD ! This is scalar on hybrid levels
2399	REAL,DIMENSION(nims:nime,njms:njme,nkms:nkme), INTENT(OUT):: C3d ! Scalar on constant pressure levels
2400	REAL,DIMENSION(nkms:nkme), INTENT(IN) :: PSTD
2401	REAL,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: PD
2402	REAL,DIMENSION(nkms:nkme), INTENT(IN) :: ETA1,ETA2
2403	REAL,INTENT(IN) :: PDTOP,PTOP
2404
2405	! local
2406
2407	INTEGER,PARAMETER :: JTB=134
2408	INTEGER :: I,J,K
2409	REAL,DIMENSION(JTB) :: PIN,CIN,Y2,PIO,PTMP,COUT,DUM1,DUM2
2410
2411	!-----------------------------------------------------------------------------------------------------
2412	!
2413	!
2414	! *** CHECK VERTICAL BOUNDS BEFORE USING SPLINE OR LINEAR INTERPOLATION
2415	!
2416	IF(nkme .GT. (JTB-10) .OR. NKDE .GT. (JTB-10)) &
2417	CALL wrf_error_fatal ('mass points: increase JTB in interp_mass_nmm')
2418
2419	!
2420	! FIRST, HORIZONTALLY INTERPOLATE MOISTURE NOW AVAILABLE ON CONSTANT PRESSURE SURFACE (LEVELS) FROM THE
2421	! PARENT TO THE NESTED DOMAIN
2422	!
2423	!*** INDEX CONVENTIONS
2424	!*** HBWGT4
2425	!*** 4
2426	!***
2427	!***
2428	!***
2429	!*** h
2430	!*** 1 2
2431	!*** HBWGT1 HBWGT2
2432	!***
2433	!***
2434	!*** 3
2435	!*** HBWGT3
2436
2437	C3d=0.0
2438	DO K=NKDS,NKDE-1 ! Please note that we are still in isobaric surfaces
2439	DO J=NJTS,MIN(NJTE,NJDE-1)
2440	DO I=NITS,MIN(NITE,NIDE-1)
2441	IF(IMASK(I,J) .NE. 1)THEN
2442	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7
2443	C3d(I,J,K) = HBWGT1(I,J)*CC3d(IIH(I,J), JJH(I,J) ,K) &
2444	+ HBWGT2(I,J)*CC3d(IIH(I,J)+1,JJH(I,J) ,K) &
2445	+ HBWGT3(I,J)*CC3d(IIH(I,J), JJH(I,J)-1,K) &
2446	+ HBWGT4(I,J)*CC3d(IIH(I,J), JJH(I,J)+1,K)
2447
2448	ELSE
2449	C3d(I,J,K) = HBWGT1(I,J)*CC3d(IIH(I,J), JJH(I,J) ,K) &
2450	+ HBWGT2(I,J)*CC3d(IIH(I,J)+1,JJH(I,J) ,K) &
2451	+ HBWGT3(I,J)*CC3d(IIH(I,J)+1,JJH(I,J)-1,K) &
2452	+ HBWGT4(I,J)*CC3d(IIH(I,J)+1,JJH(I,J)+1,K)
2453
2454	ENDIF
2455	ENDIF
2456	ENDDO
2457	ENDDO
2458	ENDDO
2459
2460	!
2461	! RECOVER THE SCALARS FROM CONSTANT PRESSURE SURFACES (LEVELS) ON TO HYBRID SURFACES
2462	!
2463	DO J=NJTS,MIN(NJTE,NJDE-1)
2464	DO I=NITS,MIN(NITE,NIDE-1)
2465	IF(IMASK(I,J) .NE. 1)THEN
2466
2467	! clean local array before use of spline or linear interpolation
2468
2469	CIN=0.;PIN=0.;Y2=0;PIO=0.;PTMP=0.;COUT=0.;DUM1=0.;DUM2=0.
2470
2471	DO K=NKDS+1,NKDE ! inputs at standard levels
2472	PIN(K-1) = EXP((ALOG(PSTD(NKDE-K+1))+ALOG(PSTD(NKDE-K+2)))*0.5)
2473	CIN(K-1) = C3d(I,J,NKDE-K+1)
2474	ENDDO
2475
2476	Y2(1 )=0.
2477	Y2(NKDE-1)=0.
2478
2479	DO K=NKDS,NKDE ! target points in model interface levels (pint)
2480	PTMP(K) = ETA1(K)PDTOP + ETA2(K)PD(I,J) + PTOP
2481	ENDDO
2482
2483	DO K=NKDS,NKDE-1 ! target points in model levels
2484	PIO(K) = EXP((ALOG(PTMP(K))+ALOG(PTMP(K+1)))*0.5)
2485	ENDDO
2486
2487	IF(PTMP(1) .GE. PSTD(1))THEN ! if lower boundary is higher than PMSL(1) re-set lower boundary
2488	PIN(NKDE-1) = PIO(1) ! be consistent with target. This may not happen at all
2489	WRITE(0,*)'WARNING: NESTED DOMAIN PRESSURE AT LOWEST LEVEL HIGHER THAN PSTD'
2490	WRITE(0,*)'I,J,PIO(1),PSTD(1)',I,J,PIO(1),PSTD(1)
2491	ENDIF
2492
2493	CALL SPLINE2(I,J,JTB,NKDE-1,PIN,CIN,Y2,NKDE-1,PIO,COUT,DUM1,DUM2) ! interpolate
2494
2495	DO K=1,NKDE-1
2496	NFLD(I,J,K)= COUT(K) ! scalar in the nested domain
2497	ENDDO
2498
2499	ENDIF
2500	ENDDO
2501	ENDDO
2502
2503	END SUBROUTINE interp_scalar_nmm
2504	!
2505	!===========================================================================================
2506	!
2507	SUBROUTINE nmm_bdy_scalar (cfld, & ! CD field
2508	cids,cide,ckds,ckde,cjds,cjde, &
2509	cims,cime,ckms,ckme,cjms,cjme, &
2510	cits,cite,ckts,ckte,cjts,cjte, &
2511	nfld, & ! ND field
2512	nids,nide,nkds,nkde,njds,njde, &
2513	nims,nime,nkms,nkme,njms,njme, &
2514	nits,nite,nkts,nkte,njts,njte, &
2515	shw, & ! stencil half width for interp
2516	imask, & ! interpolation mask
2517	xstag,ystag, & ! staggering of field
2518	ipos,jpos, & ! Position of lower left of nest in CD
2519	nri,nrj, & ! nest ratios
2520	c_bxs,n_bxs, &
2521	c_bxe,n_bxe, &
2522	c_bys,n_bys, &
2523	c_bye,n_bye, &
2524	c_btxs,n_btxs, &
2525	c_btxe,n_btxe, &
2526	c_btys,n_btys, &
2527	c_btye,n_btye, &
2528	cdt, ndt, &
2529	CTEMP_B,NTEMP_B, & ! to be removed
2530	CTEMP_BT,NTEMP_BT, &
2531	CII, IIH, CJJ, JJH, CBWGT1, HBWGT1, & ! south-western grid locs and weights
2532	CBWGT2, HBWGT2, CBWGT3, HBWGT3, & ! note that "C"ourse grid ones are
2533	CBWGT4, HBWGT4, & ! dummys for weights
2534	CC3d,C3d, &
2535	CPD,PD, &
2536	CPSTD,PSTD, &
2537	CPDTOP,PDTOP, &
2538	CPTOP,PTOP, &
2539	CETA1,ETA1,CETA2,ETA2 )
2540	USE MODULE_MODEL_CONSTANTS
2541	USE module_timing
2542	IMPLICIT NONE
2543
2544	LOGICAL,INTENT(IN) :: xstag, ystag
2545	REAL, INTENT(INOUT) :: cdt, ndt
2546	INTEGER,INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
2547	cims, cime, ckms, ckme, cjms, cjme, &
2548	cits, cite, ckts, ckte, cjts, cjte, &
2549	nids, nide, nkds, nkde, njds, njde, &
2550	nims, nime, nkms, nkme, njms, njme, &
2551	nits, nite, nkts, nkte, njts, njte, &
2552	shw,ipos,jpos,nri,nrj
2553	REAL, DIMENSION ( cims:cime, cjms:cjme, ckms:ckme ), INTENT(OUT) :: ctemp_b,ctemp_bt
2554	REAL, DIMENSION ( nims:nime, njms:njme, nkms:nkme ), INTENT(OUT) :: ntemp_b,ntemp_bt
2555	REAL, DIMENSION( * ), INTENT(INOUT) :: c_bxs,n_bxs,c_bxe,n_bxe,c_bys,n_bys,c_bye,n_bye
2556	REAL, DIMENSION( * ), INTENT(INOUT) :: c_btxs,n_btxs,c_btxe,n_btxe,c_btys,n_btys,c_btye,n_btye
2557
2558
2559	INTEGER,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: IMASK
2560
2561	! parent domain
2562
2563	INTEGER,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CII,CJJ ! dummy
2564	REAL,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CBWGT1,CBWGT2
2565	REAL,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CBWGT3,CBWGT4
2566	REAL,DIMENSION(cims:cime,cjms:cjme,ckms:ckme), INTENT(IN) :: CFLD
2567	REAL,DIMENSION(cims:cime,cjms:cjme,ckms:ckme), INTENT(IN) :: CC3d ! scalar input on constant pressure levels
2568	REAL,DIMENSION(ckms:ckme), INTENT(IN) :: CPSTD
2569	REAL,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CPD
2570	REAL,DIMENSION(ckms:ckme), INTENT(IN) :: CETA1,CETA2
2571	REAL, INTENT(IN) :: CPDTOP,CPTOP
2572
2573	! nested domain
2574
2575	INTEGER,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: IIH,JJH
2576	REAL,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: HBWGT1,HBWGT2
2577	REAL,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: HBWGT3,HBWGT4
2578	REAL,DIMENSION(nims:nime,njms:njme,nkms:nkme), INTENT(OUT):: NFLD
2579	REAL,DIMENSION(nims:nime,njms:njme,nkms:nkme), INTENT(OUT):: C3d !Scalar on constant pressure levels
2580	REAL,DIMENSION(nkms:nkme), INTENT(IN) :: PSTD
2581	REAL,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: PD
2582	REAL,DIMENSION(nkms:nkme), INTENT(IN) :: ETA1,ETA2
2583	REAL,INTENT(IN) :: PDTOP,PTOP
2584
2585	! local
2586
2587	INTEGER,PARAMETER :: JTB=134
2588	INTEGER :: I,J,K,II,JJ
2589	REAL,DIMENSION(JTB) :: PIN,CIN,Y2,PIO,PTMP,COUT,DUM1,DUM2
2590	REAL, DIMENSION (nims:nime,njms:njme,nkms:nkme) :: CWK1,CWK2,CWK3,CWK4
2591	!-----------------------------------------------------------------------------------------------------
2592	!
2593	!
2594	! *** CHECK VERTICAL BOUNDS BEFORE USING SPLINE INTERPOLATION
2595	!
2596	IF(nkme .GT. (JTB-10) .OR. NKDE .GT. (JTB-10)) &
2597	CALL wrf_error_fatal ('mass points: increase JTB in interp_mass_nmm')
2598
2599	! X start boundary
2600
2601	NMM_XS: IF(NITS .EQ. NIDS)THEN
2602	! WRITE(0,*)'ENTERING X1 START BOUNDARY AT T POINTS',NJTS,MIN(NJTE,NJDE-1)
2603	I = NIDS
2604	DO K=NKDS,NKDE-1 ! Please note that we are still in isobaric surfaces
2605	DO J = NJTS,MIN(NJTE,NJDE-1)
2606	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7 of the parent domain
2607	C3d(I,J,K) = HBWGT1(I,J)*CC3d(IIH(I,J), JJH(I,J) ,K) &
2608	+ HBWGT2(I,J)*CC3d(IIH(I,J)+1,JJH(I,J) ,K) &
2609	+ HBWGT3(I,J)*CC3d(IIH(I,J), JJH(I,J)-1,K) &
2610	+ HBWGT4(I,J)*CC3d(IIH(I,J), JJH(I,J)+1,K)
2611	ELSE
2612	C3d(I,J,K) = HBWGT1(I,J)*CC3d(IIH(I,J), JJH(I,J) ,K) &
2613	+ HBWGT2(I,J)*CC3d(IIH(I,J)+1,JJH(I,J) ,K) &
2614	+ HBWGT3(I,J)*CC3d(IIH(I,J)+1,JJH(I,J)-1,K) &
2615	+ HBWGT4(I,J)*CC3d(IIH(I,J)+1,JJH(I,J)+1,K)
2616	ENDIF
2617	ENDDO
2618	ENDDO
2619	!
2620	DO J=NJTS,MIN(NJTE,NJDE-1)
2621	IF(MOD(J,2) .NE. 0)THEN
2622	CIN=0.;PIN=0.;Y2=0;PIO=0.;PTMP=0.;COUT=0.;DUM1=0.;DUM2=0. ! clean up local array
2623	DO K=NKDS+1,NKDE ! inputs at standard levels
2624	PIN(K-1) = EXP((ALOG(PSTD(NKDE-K+1))+ALOG(PSTD(NKDE-K+2)))*0.5)
2625	CIN(K-1) = C3d(I,J,NKDE-K+1)
2626	ENDDO
2627	Y2(1 )=0.
2628	Y2(NKDE-1)=0.
2629	DO K=NKDS,NKDE ! target points in model interface levels (pint)
2630	PTMP(K) = ETA1(K)PDTOP + ETA2(K)PD(I,J) + PTOP
2631	ENDDO
2632	DO K=NKDS,NKDE-1 ! target points in model levels
2633	PIO(K) = EXP((ALOG(PTMP(K))+ALOG(PTMP(K+1)))*0.5)
2634	ENDDO
2635	IF(PTMP(1) .GE. PSTD(1))THEN ! if lower boundary is higher than PMSL(1) re-set lower boundary
2636	PIN(NKDE-1) = PIO(1) ! be consistent with target. This may not happen at all
2637	WRITE(0,*)'WARNING: NESTED DOMAIN PRESSURE AT LOWEST LEVEL HIGHER THAN PSTD'
2638	WRITE(0,*)'I,J,PIO(1),PSTD(1)',I,J,PIO(1),PSTD(1)
2639	ENDIF
2640
2641	CALL SPLINE2(I,J,JTB,NKDE-1,PIN,CIN,Y2,NKDE-1,PIO,COUT,DUM1,DUM2) ! interpolate
2642
2643	DO K=1,NKDE-1
2644	CWK1(I,J,K)= COUT(K) ! scalar in the nested domain
2645	ENDDO
2646	ELSE
2647	DO K=NKDS,NKDE-1
2648	CWK1(I,J,K)=0.0
2649	ENDDO
2650	ENDIF
2651	ENDDO
2652
2653	DO J = NJTS,MIN(NJTE,NJDE-1)
2654	DO K = NKDS,NKDE-1
2655	ntemp_b(i,j,k) = CWK1(I,J,K)
2656	ntemp_bt(i,j,k) = 0.0
2657	END DO
2658	END DO
2659
2660	ENDIF NMM_XS
2661
2662
2663	! X end boundary
2664
2665	NMM_XE: IF(NITE-1 .EQ. NIDE-1)THEN
2666	! WRITE(0,*)'ENTERING X END BOUNDARY AT T POINTS',NJTS,MIN(NJTE,NJDE-1)
2667	I = NIDE-1
2668	DO K=NKDS,NKDE-1 ! Please note that we are still in isobaric surfaces
2669	DO J = NJTS,MIN(NJTE,NJDE-1)
2670	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7 of the parent domain
2671	C3d(I,J,K) = HBWGT1(I,J)*CC3d(IIH(I,J), JJH(I,J) ,K) &
2672	+ HBWGT2(I,J)*CC3d(IIH(I,J)+1,JJH(I,J) ,K) &
2673	+ HBWGT3(I,J)*CC3d(IIH(I,J), JJH(I,J)-1,K) &
2674	+ HBWGT4(I,J)*CC3d(IIH(I,J), JJH(I,J)+1,K)
2675	ELSE
2676	C3d(I,J,K) = HBWGT1(I,J)*CC3d(IIH(I,J), JJH(I,J) ,K) &
2677	+ HBWGT2(I,J)*CC3d(IIH(I,J)+1,JJH(I,J) ,K) &
2678	+ HBWGT3(I,J)*CC3d(IIH(I,J)+1,JJH(I,J)-1,K) &
2679	+ HBWGT4(I,J)*CC3d(IIH(I,J)+1,JJH(I,J)+1,K)
2680	ENDIF
2681	ENDDO
2682	ENDDO
2683
2684	DO J=NJTS,MIN(NJTE,NJDE-1)
2685	IF(MOD(J,2) .NE. 0)THEN
2686	CIN=0.;PIN=0.;Y2=0;PIO=0.;PTMP=0.;COUT=0.;DUM1=0.;DUM2=0. ! clean up local array
2687	DO K=NKDS+1,NKDE ! inputs at standard levels
2688	PIN(K-1) = EXP((ALOG(PSTD(NKDE-K+1))+ALOG(PSTD(NKDE-K+2)))*0.5)
2689	CIN(K-1) = C3d(I,J,NKDE-K+1)
2690	ENDDO
2691	Y2(1 )=0.
2692	Y2(NKDE-1)=0.
2693	DO K=NKDS,NKDE ! target points in model interface levels (pint)
2694	PTMP(K) = ETA1(K)PDTOP + ETA2(K)PD(I,J) + PTOP
2695	ENDDO
2696	DO K=NKDS,NKDE-1 ! target points in model levels
2697	PIO(K) = EXP((ALOG(PTMP(K))+ALOG(PTMP(K+1)))*0.5)
2698	ENDDO
2699	IF(PTMP(1) .GE. PSTD(1))THEN ! if lower boundary is higher than PMSL(1) re-set lower boundary
2700	PIN(NKDE-1) = PIO(1) ! be consistent with target. This may not happen at all
2701	WRITE(0,*)'WARNING: NESTED DOMAIN PRESSURE AT LOWEST LEVEL HIGHER THAN PSTD'
2702	WRITE(0,*)'I,J,PIO(1),PSTD(1)',I,J,PIO(1),PSTD(1)
2703	ENDIF
2704
2705	CALL SPLINE2(I,J,JTB,NKDE-1,PIN,CIN,Y2,NKDE-1,PIO,COUT,DUM1,DUM2) ! interpolate
2706
2707	DO K=1,NKDE-1
2708	CWK2(I,J,K)= COUT(K) ! scalar in the nested domain
2709	ENDDO
2710	ELSE
2711	DO K=NKDS,NKDE-1
2712	CWK2(I,J,K)=0.0
2713	ENDDO
2714	ENDIF
2715	ENDDO
2716
2717	DO J = NJTS,MIN(NJTE,NJDE-1)
2718	DO K = NKDS,MIN(NKTE,NKDE-1)
2719	ntemp_b(i,j,k) = CWK2(I,J,K)
2720	ntemp_bt(i,j,k) = 0.0
2721	END DO
2722	END DO
2723
2724	ENDIF NMM_XE
2725
2726	! Y start boundary
2727
2728	NMM_YS: IF(NJTS .EQ. NJDS)THEN
2729	! WRITE(0,*)'ENTERING Y START BOUNDARY AT T POINTS',NITS,MIN(NITE,NIDE-1)
2730	J = NJDS
2731	DO K=NKDS,NKDE-1
2732	DO I = NITS,MIN(NITE,NIDE-1)
2733	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7 of the parent domain
2734	C3d(I,J,K) = HBWGT1(I,J)*CC3d(IIH(I,J), JJH(I,J) ,K) &
2735	+ HBWGT2(I,J)*CC3d(IIH(I,J)+1,JJH(I,J) ,K) &
2736	+ HBWGT3(I,J)*CC3d(IIH(I,J), JJH(I,J)-1,K) &
2737	+ HBWGT4(I,J)*CC3d(IIH(I,J), JJH(I,J)+1,K)
2738	ELSE
2739	C3d(I,J,K) = HBWGT1(I,J)*CC3d(IIH(I,J), JJH(I,J) ,K) &
2740	+ HBWGT2(I,J)*CC3d(IIH(I,J)+1,JJH(I,J) ,K) &
2741	+ HBWGT3(I,J)*CC3d(IIH(I,J)+1,JJH(I,J)-1,K) &
2742	+ HBWGT4(I,J)*CC3d(IIH(I,J)+1,JJH(I,J)+1,K)
2743
2744	ENDIF
2745	ENDDO
2746	ENDDO
2747	!
2748	DO I=NITS,MIN(NITE,NIDE-1)
2749	CIN=0.;PIN=0.;Y2=0;PIO=0.;PTMP=0.;COUT=0.;DUM1=0.;DUM2=0. ! clean up local array
2750	DO K=NKDS+1,NKDE ! inputs at standard levels
2751	PIN(K-1) = EXP((ALOG(PSTD(NKDE-K+1))+ALOG(PSTD(NKDE-K+2)))*0.5)
2752	CIN(K-1) = C3d(I,J,NKDE-K+1)
2753	ENDDO
2754	Y2(1 )=0.
2755	Y2(NKDE-1)=0.
2756	DO K=NKDS,NKDE ! target points in model interface levels (pint)
2757	PTMP(K) = ETA1(K)PDTOP + ETA2(K)PD(I,J) + PTOP
2758	ENDDO
2759	DO K=NKDS,NKDE-1 ! target points in model levels
2760	PIO(K) = EXP((ALOG(PTMP(K))+ALOG(PTMP(K+1)))*0.5)
2761	ENDDO
2762	IF(PTMP(1) .GE. PSTD(1))THEN ! if lower boundary is higher than PMSL(1) re-set lower boundary
2763	PIN(NKDE-1) = PIO(1) ! be consistent with target. This may not happen at all
2764	WRITE(0,*)'WARNING: NESTED DOMAIN PRESSURE AT LOWEST LEVEL HIGHER THAN PSTD'
2765	WRITE(0,*)'I,J,PIO(1),PSTD(1)',I,J,PIO(1),PSTD(1)
2766	ENDIF
2767
2768	CALL SPLINE2(I,J,JTB,NKDE-1,PIN,CIN,Y2,NKDE-1,PIO,COUT,DUM1,DUM2) ! interpolate
2769
2770	DO K=1,NKDE-1
2771	CWK3(I,J,K)= COUT(K) ! scalar in the nested domain
2772	ENDDO
2773	ENDDO
2774
2775	DO K = NKDS,NKDE-1
2776	DO I = NITS,MIN(NITE,NIDE-1)
2777	ntemp_b(i,J,K) = CWK3(I,J,K)
2778	ntemp_bt(i,J,K) = 0.0
2779	ENDDO
2780	ENDDO
2781
2782	ENDIF NMM_YS
2783
2784	! Y end boundary
2785
2786	NMM_YE: IF(NJTE-1 .EQ. NJDE-1)THEN
2787	! WRITE(0,*)'ENTERING Y END BOUNDARY AT T POINTS',NITS,MIN(NITE,NIDE-1)
2788	J = NJDE-1
2789	DO K=NKDS,NKDE-1
2790	DO I = NITS,MIN(NITE,NIDE-1)
2791	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7 of the parent domain
2792	C3d(I,J,K) = HBWGT1(I,J)*CC3d(IIH(I,J), JJH(I,J) ,K) &
2793	+ HBWGT2(I,J)*CC3d(IIH(I,J)+1,JJH(I,J) ,K) &
2794	+ HBWGT3(I,J)*CC3d(IIH(I,J), JJH(I,J)-1,K) &
2795	+ HBWGT4(I,J)*CC3d(IIH(I,J), JJH(I,J)+1,K)
2796	ELSE
2797	C3d(I,J,K) = HBWGT1(I,J)*CC3d(IIH(I,J), JJH(I,J) ,K) &
2798	+ HBWGT2(I,J)*CC3d(IIH(I,J)+1,JJH(I,J) ,K) &
2799	+ HBWGT3(I,J)*CC3d(IIH(I,J)+1,JJH(I,J)-1,K) &
2800	+ HBWGT4(I,J)*CC3d(IIH(I,J)+1,JJH(I,J)+1,K)
2801
2802	ENDIF
2803	ENDDO
2804	ENDDO
2805
2806	DO I=NITS,MIN(NITE,NIDE-1)
2807	CIN=0.;PIN=0.;Y2=0;PIO=0.;PTMP=0.;COUT=0.;DUM1=0.;DUM2=0. ! clean up local array
2808	DO K=NKDS+1,NKDE ! inputs at standard levels
2809	PIN(K-1) = EXP((ALOG(PSTD(NKDE-K+1))+ALOG(PSTD(NKDE-K+2)))*0.5)
2810	CIN(K-1) = C3d(I,J,NKDE-K+1)
2811	ENDDO
2812	Y2(1 )=0.
2813	Y2(NKDE-1)=0.
2814	DO K=NKDS,NKDE ! target points in model interface levels (pint)
2815	PTMP(K) = ETA1(K)PDTOP + ETA2(K)PD(I,J) + PTOP
2816	ENDDO
2817	DO K=NKDS,NKDE-1 ! target points in model levels
2818	PIO(K) = EXP((ALOG(PTMP(K))+ALOG(PTMP(K+1)))*0.5)
2819	ENDDO
2820	IF(PTMP(1) .GE. PSTD(1))THEN ! if lower boundary is higher than PMSL(1) re-set lower boundary
2821	PIN(NKDE-1) = PIO(1) ! be consistent with target. This may not happen at all
2822	WRITE(0,*)'WARNING: NESTED DOMAIN PRESSURE AT LOWEST LEVEL HIGHER THAN PSTD'
2823	WRITE(0,*)'I,J,PIO(1),PSTD(1)',I,J,PIO(1),PSTD(1)
2824	ENDIF
2825
2826	CALL SPLINE2(I,J,JTB,NKDE-1,PIN,CIN,Y2,NKDE-1,PIO,COUT,DUM1,DUM2) ! interpolate
2827
2828	DO K=1,NKDE-1
2829	CWK4(I,J,K)= COUT(K) ! scalar in the nested domain
2830	ENDDO
2831	ENDDO
2832
2833	DO K = NKDS,NKDE-1
2834	DO I = NITS,MIN(NITE,NIDE-1)
2835	ntemp_b(i,J,K) = CWK4(I,J,K)
2836	ntemp_bt(i,J,K) = 0.0
2837	END DO
2838	END DO
2839
2840	ENDIF NMM_YE
2841
2842	END SUBROUTINE nmm_bdy_scalar
2843	!
2844	!
2845	!=======================================================================================
2846	SUBROUTINE SPLINE2(I,J,JTBX,NOLD,XOLD,YOLD,Y2,NNEW,XNEW,YNEW,P,Q)
2847	!
2848	! ******************************************************************
2849	! * *
2850	! * THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE *
2851	! * PROGRAMED FOR A SMALL SCALAR MACHINE. *
2852	! * *
2853	! * PROGRAMER Z. JANJIC *
2854	! * *
2855	! * NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION. MUST BE GE 3. *
2856	! * XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
2857	! * FUNCTION ARE GIVEN. MUST BE IN ASCENDING ORDER. *
2858	! * YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD. *
2859	! * Y2 - THE SECOND DERIVATIVES AT THE POINTS XOLD. IF NATURAL *
2860	! * SPLINE IS FITTED Y2(1)=0. AND Y2(NOLD)=0. MUST BE *
2861	! * SPECIFIED. *
2862	! * NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED. *
2863	! * XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
2864	! * FUNCTION ARE CALCULATED. XNEW(K) MUST BE GE XOLD(1) *
2865	! * AND LE XOLD(NOLD). *
2866	! * YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED. *
2867	! * P, Q - AUXILIARY VECTORS OF THE LENGTH NOLD-2. *
2868	! * *
2869	! ******************************************************************
2870	!---------------------------------------------------------------------
2871	IMPLICIT NONE
2872	!---------------------------------------------------------------------
2873	INTEGER,INTENT(IN) :: I,J,JTBX,NNEW,NOLD
2874	REAL,DIMENSION(JTBX),INTENT(IN) :: XNEW,XOLD,YOLD
2875	REAL,DIMENSION(JTBX),INTENT(INOUT) :: P,Q,Y2
2876	REAL,DIMENSION(JTBX),INTENT(OUT) :: YNEW
2877	!
2878	INTEGER :: II,JJ,K,K1,K2,KOLD,NOLDM1
2879	REAL :: AK,BK,CK,DEN,DX,DXC,DXL,DXR,DYDXL,DYDXR &
2880	,RDX,RTDXC,X,XK,XSQ,Y2K,Y2KP1
2881	!---------------------------------------------------------------------
2882
2883	! debug
2884
2885	II=9999
2886	JJ=9999
2887	IF(I.eq.II.and.J.eq.JJ)THEN
2888	WRITE(0,*)'DEBUG in SPLINE2: I,J',I,J
2889	WRITE(0,*)'DEBUG in SPLINE2:HSO= ',xnew(1:nold)
2890	DO K=1,NOLD
2891	WRITE(0,*)'DEBUG in SPLINE2:L,ZETAI,PINTI= ' &
2892	,K,YOLD(K),XOLD(K)
2893	ENDDO
2894	ENDIF
2895	!
2896	NOLDM1=NOLD-1
2897	!
2898	DXL=XOLD(2)-XOLD(1)
2899	DXR=XOLD(3)-XOLD(2)
2900	DYDXL=(YOLD(2)-YOLD(1))/DXL
2901	DYDXR=(YOLD(3)-YOLD(2))/DXR
2902	RTDXC=0.5/(DXL+DXR)
2903	!
2904	P(1)= RTDXC(6.(DYDXR-DYDXL)-DXL*Y2(1))
2905	Q(1)=-RTDXC*DXR
2906	!
2907	IF(NOLD.EQ.3)GO TO 150
2908	!---------------------------------------------------------------------
2909	K=3
2910	!
2911	100 DXL=DXR
2912	DYDXL=DYDXR
2913	DXR=XOLD(K+1)-XOLD(K)
2914	DYDXR=(YOLD(K+1)-YOLD(K))/DXR
2915	DXC=DXL+DXR
2916	DEN=1./(DXL*Q(K-2)+DXC+DXC)
2917	!
2918	P(K-1)= DEN(6.(DYDXR-DYDXL)-DXL*P(K-2))
2919	Q(K-1)=-DEN*DXR
2920	!
2921	K=K+1
2922	IF(K.LT.NOLD)GO TO 100
2923	!-----------------------------------------------------------------------
2924	150 K=NOLDM1
2925	!
2926	200 Y2(K)=P(K-1)+Q(K-1)*Y2(K+1)
2927	!
2928	K=K-1
2929	IF(K.GT.1)GO TO 200
2930	!-----------------------------------------------------------------------
2931	K1=1
2932	!
2933	300 XK=XNEW(K1)
2934	!
2935	DO 400 K2=2,NOLD
2936	!
2937	IF(XOLD(K2).GT.XK)THEN
2938	KOLD=K2-1
2939	GO TO 450
2940	ENDIF
2941	!
2942	400 CONTINUE
2943	!
2944	YNEW(K1)=YOLD(NOLD)
2945	GO TO 600
2946	!
2947	450 IF(K1.EQ.1)GO TO 500
2948	IF(K.EQ.KOLD)GO TO 550
2949	!
2950	500 K=KOLD
2951	!
2952	Y2K=Y2(K)
2953	Y2KP1=Y2(K+1)
2954	DX=XOLD(K+1)-XOLD(K)
2955	RDX=1./DX
2956	!
2957	AK=.1666667RDX(Y2KP1-Y2K)
2958	BK=0.5*Y2K
2959	CK=RDX(YOLD(K+1)-YOLD(K))-.1666667DX*(Y2KP1+Y2K+Y2K)
2960	!
2961	550 X=XK-XOLD(K)
2962	XSQ=X*X
2963	!
2964	YNEW(K1)=AKXSQX+BKXSQ+CKX+YOLD(K)
2965
2966	! debug
2967
2968	IF(I.eq.II.and.J.eq.JJ)THEN
2969	WRITE(0,*) 'DEBUG:: k1,xnew(k1),ynew(k1): ', K1,XNEW(k1),YNEW(k1)
2970	ENDIF
2971
2972	!
2973	600 K1=K1+1
2974	IF(K1.LE.NNEW)GO TO 300
2975
2976	RETURN
2977
2978	END SUBROUTINE SPLINE2
2979
2980	!=======================================================================================
2981	! E grid interpolation for H and V points
2982	!=======================================================================================
2983
2984	SUBROUTINE interp_h_nmm (cfld, & ! CD field
2985	cids, cide, ckds, ckde, cjds, cjde, &
2986	cims, cime, ckms, ckme, cjms, cjme, &
2987	cits, cite, ckts, ckte, cjts, cjte, &
2988	nfld, & ! ND field
2989	nids, nide, nkds, nkde, njds, njde, &
2990	nims, nime, nkms, nkme, njms, njme, &
2991	nits, nite, nkts, nkte, njts, njte, &
2992	shw, & ! stencil half width for interp
2993	imask, & ! interpolation mask
2994	xstag, ystag, & ! staggering of field
2995	ipos, jpos, & ! Position of lower left of nest in CD
2996	nri, nrj, & ! nest ratios
2997	CII, IIH, CJJ, JJH, CBWGT1, HBWGT1, & ! south-western grid locs and weights
2998	CBWGT2, HBWGT2, CBWGT3, HBWGT3, & ! note that "C"ourse grid ones are
2999	CBWGT4, HBWGT4 ) ! dummys for weights
3000	USE module_timing
3001	IMPLICIT NONE
3002
3003	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3004	cims, cime, ckms, ckme, cjms, cjme, &
3005	cits, cite, ckts, ckte, cjts, cjte, &
3006	nids, nide, nkds, nkde, njds, njde, &
3007	nims, nime, nkms, nkme, njms, njme, &
3008	nits, nite, nkts, nkte, njts, njte, &
3009	shw, &
3010	ipos, jpos, &
3011	nri, nrj
3012	LOGICAL, INTENT(IN) :: xstag, ystag
3013
3014	REAL, DIMENSION ( cims:cime, cjms:cjme, ckms:ckme ) :: cfld
3015	REAL, DIMENSION ( nims:nime, njms:njme, nkms:nkme ) :: nfld
3016	REAL, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CBWGT1,CBWGT2,CBWGT3,CBWGT4 ! dummy
3017	REAL, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
3018	INTEGER, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CII,CJJ ! dummy
3019	INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: IIH,JJH
3020	INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
3021
3022	! local
3023	INTEGER i,j,k
3024	!
3025	!*** CHECK DOMAIN BOUNDS BEFORE INTERPOLATION
3026	!
3027	DO J=NJTS,MIN(NJTE,NJDE-1)
3028	DO I=NITS,MIN(NITE,NIDE-1)
3029	IF(IIH(i,j).LT.(CIDS-shw) .OR. IIH(i,j).GT.(CIDE+shw)) &
3030	CALL wrf_error_fatal ('hpoints:check domain bounds along x' )
3031	IF(JJH(i,j).LT.(CJDS-shw) .OR. JJH(i,j).GT.(CJDE+shw)) &
3032	CALL wrf_error_fatal ('hpoints:check domain bounds along y' )
3033	ENDDO
3034	ENDDO
3035	!
3036	!*** INDEX CONVENTIONS
3037	!*** HBWGT4
3038	!*** 4
3039	!***
3040	!***
3041	!***
3042	!*** h
3043	!*** 1 2
3044	!*** HBWGT1 HBWGT2
3045	!***
3046	!***
3047	!*** 3
3048	!*** HBWGT3
3049
3050	DO K=NKDS,NKDE
3051	DO J=NJTS,MIN(NJTE,NJDE-1)
3052	DO I=NITS,MIN(NITE,NIDE-1)
3053	IF(IMASK(I,J) .NE. 1)THEN
3054	!
3055	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7
3056	NFLD(I,J,K) = HBWGT1(I,J)*CFLD(IIH(I,J), JJH(I,J) ,K) &
3057	+ HBWGT2(I,J)*CFLD(IIH(I,J)+1,JJH(I,J) ,K) &
3058	+ HBWGT3(I,J)*CFLD(IIH(I,J), JJH(I,J)-1,K) &
3059	+ HBWGT4(I,J)*CFLD(IIH(I,J), JJH(I,J)+1,K)
3060	ELSE
3061	NFLD(I,J,K) = HBWGT1(I,J)*CFLD(IIH(I,J), JJH(I,J) ,K) &
3062	+ HBWGT2(I,J)*CFLD(IIH(I,J)+1,JJH(I,J) ,K) &
3063	+ HBWGT3(I,J)*CFLD(IIH(I,J)+1,JJH(I,J)-1,K) &
3064	+ HBWGT4(I,J)*CFLD(IIH(I,J)+1,JJH(I,J)+1,K)
3065	ENDIF
3066	!
3067	ENDIF
3068	ENDDO
3069	ENDDO
3070	ENDDO
3071
3072	END SUBROUTINE interp_h_nmm
3073	!
3074	SUBROUTINE interp_v_nmm (cfld, & ! CD field
3075	cids, cide, ckds, ckde, cjds, cjde, &
3076	cims, cime, ckms, ckme, cjms, cjme, &
3077	cits, cite, ckts, ckte, cjts, cjte, &
3078	nfld, & ! ND field
3079	nids, nide, nkds, nkde, njds, njde, &
3080	nims, nime, nkms, nkme, njms, njme, &
3081	nits, nite, nkts, nkte, njts, njte, &
3082	shw, & ! stencil half width for interp
3083	imask, & ! interpolation mask
3084	xstag, ystag, & ! staggering of field
3085	ipos, jpos, & ! Position of lower left of nest in CD
3086	nri, nrj, & ! nest ratios
3087	CII, IIV, CJJ, JJV, CBWGT1, VBWGT1, & ! south-western grid locs and weights
3088	CBWGT2, VBWGT2, CBWGT3, VBWGT3, & ! note that "C"ourse grid ones are
3089	CBWGT4, VBWGT4 ) ! dummys
3090	USE module_timing
3091	IMPLICIT NONE
3092
3093	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3094	cims, cime, ckms, ckme, cjms, cjme, &
3095	cits, cite, ckts, ckte, cjts, cjte, &
3096	nids, nide, nkds, nkde, njds, njde, &
3097	nims, nime, nkms, nkme, njms, njme, &
3098	nits, nite, nkts, nkte, njts, njte, &
3099	shw, &
3100	ipos, jpos, &
3101	nri, nrj
3102	LOGICAL, INTENT(IN) :: xstag, ystag
3103
3104	REAL, DIMENSION ( cims:cime, cjms:cjme, ckms:ckme ) :: cfld
3105	REAL, DIMENSION ( nims:nime, njms:njme, nkms:nkme ) :: nfld
3106	REAL, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CBWGT1,CBWGT2,CBWGT3,CBWGT4 ! dummy
3107	REAL, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
3108	INTEGER, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CII,CJJ ! dummy
3109	INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: IIV,JJV
3110	INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
3111
3112	! local
3113	INTEGER i,j,k
3114
3115
3116	!
3117	!*** CHECK DOMAIN BOUNDS BEFORE INTERPOLATION
3118	!
3119	DO J=NJTS,MIN(NJTE,NJDE-1)
3120	DO I=NITS,MIN(NITE,NIDE-1)
3121	IF(IIV(i,j).LT.(CIDS-shw) .OR. IIV(i,j).GT.(CIDE+shw)) &
3122	CALL wrf_error_fatal ('vpoints:check domain bounds along x' )
3123	IF(JJV(i,j).LT.(CJDS-shw) .OR. JJV(i,j).GT.(CJDE+shw)) &
3124	CALL wrf_error_fatal ('vpoints:check domain bounds along y' )
3125	ENDDO
3126	ENDDO
3127	!
3128	!*** INDEX CONVENTIONS
3129	!*** VBWGT4
3130	!*** 4
3131	!***
3132	!***
3133	!***
3134	!*** h
3135	!*** 1 2
3136	!*** VBWGT1 VBWGT2
3137	!***
3138	!***
3139	!*** 3
3140	!*** VBWGT3
3141
3142	DO K=NKDS,NKDE
3143	DO J=NJTS,MIN(NJTE,NJDE-1)
3144	DO I=NITS,MIN(NITE,NIDE-1)
3145	IF(IMASK(I,J) .NE. 1)THEN
3146
3147	IF(MOD(JJV(I,J),2) .NE. 0)THEN ! 1,3,5,7
3148	NFLD(I,J,K) = VBWGT1(I,J)*CFLD(IIV(I,J), JJV(I,J) ,K) &
3149	+ VBWGT2(I,J)*CFLD(IIV(I,J)+1,JJV(I,J) ,K) &
3150	+ VBWGT3(I,J)*CFLD(IIV(I,J)+1,JJV(I,J)-1,K) &
3151	+ VBWGT4(I,J)*CFLD(IIV(I,J)+1,JJV(I,J)+1,K)
3152	ELSE
3153	NFLD(I,J,K) = VBWGT1(I,J)*CFLD(IIV(I,J), JJV(I,J) ,K) &
3154	+ VBWGT2(I,J)*CFLD(IIV(I,J)+1,JJV(I,J) ,K) &
3155	+ VBWGT3(I,J)*CFLD(IIV(I,J), JJV(I,J)-1,K) &
3156	+ VBWGT4(I,J)*CFLD(IIV(I,J), JJV(I,J)+1,K)
3157	ENDIF
3158
3159	ENDIF
3160	ENDDO
3161	ENDDO
3162	ENDDO
3163
3164	END SUBROUTINE interp_v_nmm
3165	!
3166	!=======================================================================================
3167	! E grid nearest neighbour interpolation for H points.
3168	! This routine assumes cfld and nfld are in IJK
3169	!=======================================================================================
3170	!
3171	SUBROUTINE interp_hnear_nmm (cfld, & ! CD field
3172	cids, cide, ckds, ckde, cjds, cjde, &
3173	cims, cime, ckms, ckme, cjms, cjme, &
3174	cits, cite, ckts, ckte, cjts, cjte, &
3175	nfld, & ! ND field
3176	nids, nide, nkds, nkde, njds, njde, &
3177	nims, nime, nkms, nkme, njms, njme, &
3178	nits, nite, nkts, nkte, njts, njte, &
3179	shw, & ! stencil half width for interp
3180	imask, & ! interpolation mask
3181	xstag, ystag, & ! staggering of field
3182	ipos, jpos, & ! Position of lower left of nest in CD
3183	nri, nrj, & ! nest ratios
3184	CII, IIH, CJJ, JJH, CBWGT1, HBWGT1, & ! south-western grid locs and weights
3185	CBWGT2, HBWGT2, CBWGT3, HBWGT3, & ! note that "C"ourse grid ones are
3186	CBWGT4, HBWGT4 ) ! just dummys
3187	USE module_timing
3188	IMPLICIT NONE
3189
3190	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3191	cims, cime, ckms, ckme, cjms, cjme, &
3192	cits, cite, ckts, ckte, cjts, cjte, &
3193	nids, nide, nkds, nkde, njds, njde, &
3194	nims, nime, nkms, nkme, njms, njme, &
3195	nits, nite, nkts, nkte, njts, njte, &
3196	shw, &
3197	ipos, jpos, &
3198	nri, nrj
3199	LOGICAL, INTENT(IN) :: xstag, ystag
3200
3201	REAL, DIMENSION ( cims:cime, cjms:cjme, ckms:ckme ) :: cfld
3202	REAL, DIMENSION ( nims:nime, njms:njme, nkms:nkme ) :: nfld
3203	REAL, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CBWGT1,CBWGT2,CBWGT3,CBWGT4 ! dummy
3204	REAL, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
3205	INTEGER, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CII,CJJ ! dummy
3206	INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: IIH,JJH
3207	INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
3208
3209	! local
3210
3211	LOGICAL FLIP
3212	INTEGER i,j,k,n
3213	REAL SUM,AMAXVAL
3214	REAL, DIMENSION (4, nims:nime, njms:njme ) :: NBWGT
3215
3216	!
3217	!*** INDEX CONVENTIONS
3218	!*** NBWGT4=0
3219	!*** 4
3220	!***
3221	!***
3222	!***
3223	!*** h
3224	!*** 1 2
3225	!*** NBWGT1=1 NBWGT2=0
3226	!***
3227	!***
3228	!*** 3
3229	!*** NBWGT3=0
3230
3231	DO J=NJTS,MIN(NJTE,NJDE-1)
3232	DO I=NITS,MIN(NITE,NIDE-1)
3233	IF(IMASK(I,J) .NE. 1)THEN
3234	NBWGT(1,I,J)=HBWGT1(I,J)
3235	NBWGT(2,I,J)=HBWGT2(I,J)
3236	NBWGT(3,I,J)=HBWGT3(I,J)
3237	NBWGT(4,I,J)=HBWGT4(I,J)
3238	ENDIF
3239	ENDDO
3240	ENDDO
3241
3242	DO J=NJTS,MIN(NJTE,NJDE-1)
3243	DO I=NITS,MIN(NITE,NIDE-1)
3244	IF(IMASK(I,J) .NE. 1)THEN
3245
3246	AMAXVAL=0.
3247	DO N=1,4
3248	AMAXVAL=amax1(NBWGT(N,I,J),AMAXVAL)
3249	ENDDO
3250
3251	FLIP=.TRUE.
3252	SUM=0.0
3253	DO N=1,4
3254	IF(AMAXVAL .EQ. NBWGT(N,I,J) .AND. FLIP)THEN
3255	NBWGT(N,I,J)=1.0
3256	FLIP=.FALSE.
3257	ELSE
3258	NBWGT(N,I,J)=0.0
3259	ENDIF
3260	SUM=SUM+NBWGT(N,I,J)
3261	IF(SUM .GT. 1.0)CALL wrf_error_fatal ( "horizontal interp error - interp_hnear_nmm" )
3262	ENDDO
3263
3264	ENDIF
3265	ENDDO
3266	ENDDO
3267
3268	DO K=NKDS,NKDE
3269	DO J=NJTS,MIN(NJTE,NJDE-1)
3270	DO I=NITS,MIN(NITE,NIDE-1)
3271	IF(IMASK(I,J) .NE. 1)THEN
3272	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7
3273	NFLD(I,J,K) = NBWGT(1,I,J)*CFLD(IIH(I,J), JJH(I,J) ,K) &
3274	+ NBWGT(2,I,J)*CFLD(IIH(I,J)+1,JJH(I,J) ,K) &
3275	+ NBWGT(3,I,J)*CFLD(IIH(I,J), JJH(I,J)-1,K) &
3276	+ NBWGT(4,I,J)*CFLD(IIH(I,J), JJH(I,J)+1,K)
3277	ELSE
3278	NFLD(I,J,K) = NBWGT(1,I,J)*CFLD(IIH(I,J), JJH(I,J) ,K) &
3279	+ NBWGT(2,I,J)*CFLD(IIH(I,J)+1,JJH(I,J) ,K) &
3280	+ NBWGT(3,I,J)*CFLD(IIH(I,J)+1,JJH(I,J)-1,K) &
3281	+ NBWGT(4,I,J)*CFLD(IIH(I,J)+1,JJH(I,J)+1,K)
3282	ENDIF
3283	ENDIF
3284	ENDDO
3285	ENDDO
3286	ENDDO
3287
3288	END SUBROUTINE interp_hnear_nmm
3289	!
3290	!=======================================================================================
3291	! E grid nearest neighbour interpolation for H points.
3292	! This routine assumes cfld and nfld are in IKJ or ILJ
3293	!=======================================================================================
3294	!
3295	SUBROUTINE interp_hnear_ikj_nmm (cfld, & ! CD field
3296	cids, cide, ckds, ckde, cjds, cjde, &
3297	cims, cime, ckms, ckme, cjms, cjme, &
3298	cits, cite, ckts, ckte, cjts, cjte, &
3299	nfld, & ! ND field
3300	nids, nide, nkds, nkde, njds, njde, &
3301	nims, nime, nkms, nkme, njms, njme, &
3302	nits, nite, nkts, nkte, njts, njte, &
3303	shw, & ! stencil half width for interp
3304	imask, & ! interpolation mask
3305	xstag, ystag, & ! staggering of field
3306	ipos, jpos, & ! Position of lower left of nest in CD
3307	nri, nrj, & ! nest ratios
3308	CII, IIH, CJJ, JJH, CBWGT1, HBWGT1, & ! south-western grid locs and weights
3309	CBWGT2, HBWGT2, CBWGT3, HBWGT3, & ! note that "C"ourse grid ones are
3310	CBWGT4, HBWGT4 ) ! just dummys
3311	USE module_timing
3312	IMPLICIT NONE
3313
3314	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3315	cims, cime, ckms, ckme, cjms, cjme, &
3316	cits, cite, ckts, ckte, cjts, cjte, &
3317	nids, nide, nkds, nkde, njds, njde, &
3318	nims, nime, nkms, nkme, njms, njme, &
3319	nits, nite, nkts, nkte, njts, njte, &
3320	shw, &
3321	ipos, jpos, &
3322	nri, nrj
3323	LOGICAL, INTENT(IN) :: xstag, ystag
3324
3325	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfld
3326	REAL, DIMENSION ( nims:nime, nkms:nkme, njms:njme ) :: nfld
3327	REAL, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CBWGT1,CBWGT2,CBWGT3,CBWGT4 ! dummy
3328	REAL, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
3329	INTEGER, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CII,CJJ ! dummy
3330	INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: IIH,JJH
3331	INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
3332
3333	! local
3334
3335	LOGICAL FLIP
3336	INTEGER i,j,k,n
3337	REAL SUM,AMAXVAL
3338	REAL, DIMENSION (4, nims:nime, njms:njme ) :: NBWGT
3339
3340	!
3341	!*** INDEX CONVENTIONS
3342	!*** NBWGT4=0
3343	!*** 4
3344	!***
3345	!***
3346	!***
3347	!*** h
3348	!*** 1 2
3349	!*** NBWGT1=1 NBWGT2=0
3350	!***
3351	!***
3352	!*** 3
3353	!*** NBWGT3=0
3354
3355	DO J=NJTS,MIN(NJTE,NJDE-1)
3356	DO I=NITS,MIN(NITE,NIDE-1)
3357	IF(IMASK(I,J) .NE. 1)THEN
3358	NBWGT(1,I,J)=HBWGT1(I,J)
3359	NBWGT(2,I,J)=HBWGT2(I,J)
3360	NBWGT(3,I,J)=HBWGT3(I,J)
3361	NBWGT(4,I,J)=HBWGT4(I,J)
3362	ENDIF
3363	ENDDO
3364	ENDDO
3365
3366	DO J=NJTS,MIN(NJTE,NJDE-1)
3367	DO I=NITS,MIN(NITE,NIDE-1)
3368	IF(IMASK(I,J) .NE. 1)THEN
3369
3370	AMAXVAL=0.
3371	DO N=1,4
3372	AMAXVAL=amax1(NBWGT(N,I,J),AMAXVAL)
3373	ENDDO
3374
3375	FLIP=.TRUE.
3376	SUM=0.0
3377	DO N=1,4
3378	IF(AMAXVAL .EQ. NBWGT(N,I,J) .AND. FLIP)THEN
3379	NBWGT(N,I,J)=1.0
3380	FLIP=.FALSE.
3381	ELSE
3382	NBWGT(N,I,J)=0.0
3383	ENDIF
3384	SUM=SUM+NBWGT(N,I,J)
3385	IF(SUM .GT. 1.0)CALL wrf_error_fatal ( "horizontal interp error - interp_hnear_nmm" )
3386	ENDDO
3387
3388	ENDIF
3389	ENDDO
3390	ENDDO
3391
3392	DO J=NJTS,MIN(NJTE,NJDE-1)
3393	DO K=NKDS,NKDE
3394	DO I=NITS,MIN(NITE,NIDE-1)
3395	IF(IMASK(I,J) .NE. 1)THEN
3396	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7
3397	NFLD(I,K,J) = NBWGT(1,I,J)*CFLD(IIH(I,J), K,JJH(I,J) ) &
3398	+ NBWGT(2,I,J)*CFLD(IIH(I,J)+1,K,JJH(I,J) ) &
3399	+ NBWGT(3,I,J)*CFLD(IIH(I,J), K,JJH(I,J)-1) &
3400	+ NBWGT(4,I,J)*CFLD(IIH(I,J), K,JJH(I,J)+1)
3401	ELSE
3402	NFLD(I,K,J) = NBWGT(1,I,J)*CFLD(IIH(I,J), K,JJH(I,J) ) &
3403	+ NBWGT(2,I,J)*CFLD(IIH(I,J)+1,K,JJH(I,J) ) &
3404	+ NBWGT(3,I,J)*CFLD(IIH(I,J)+1,K,JJH(I,J)-1) &
3405	+ NBWGT(4,I,J)*CFLD(IIH(I,J)+1,K,JJH(I,J)+1)
3406	ENDIF
3407	ENDIF
3408	ENDDO
3409	ENDDO
3410	ENDDO
3411
3412	END SUBROUTINE interp_hnear_ikj_nmm
3413	!
3414	!=======================================================================================
3415	! E grid nearest neighbour interpolation for integer H points
3416	!=======================================================================================
3417	!
3418	SUBROUTINE interp_int_hnear_nmm (cfld, & ! CD field; integers
3419	cids, cide, ckds, ckde, cjds, cjde, &
3420	cims, cime, ckms, ckme, cjms, cjme, &
3421	cits, cite, ckts, ckte, cjts, cjte, &
3422	nfld, & ! ND field; integers
3423	nids, nide, nkds, nkde, njds, njde, &
3424	nims, nime, nkms, nkme, njms, njme, &
3425	nits, nite, nkts, nkte, njts, njte, &
3426	shw, & ! stencil half width for interp
3427	imask, & ! interpolation mask
3428	xstag, ystag, & ! staggering of field
3429	ipos, jpos, & ! lower left of nest in CD
3430	nri, nrj, & ! nest ratios
3431	CII, IIH, CJJ, JJH, CBWGT1, HBWGT1, & ! s-w grid locs and weights
3432	CBWGT2, HBWGT2, CBWGT3, HBWGT3, & ! note that "C"ourse grid ones are
3433	CBWGT4, HBWGT4 ) ! just dummys
3434	USE module_timing
3435	IMPLICIT NONE
3436
3437	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3438	cims, cime, ckms, ckme, cjms, cjme, &
3439	cits, cite, ckts, ckte, cjts, cjte, &
3440	nids, nide, nkds, nkde, njds, njde, &
3441	nims, nime, nkms, nkme, njms, njme, &
3442	nits, nite, nkts, nkte, njts, njte, &
3443	shw, &
3444	ipos, jpos, &
3445	nri, nrj
3446	LOGICAL, INTENT(IN) :: xstag, ystag
3447
3448	INTEGER, DIMENSION ( cims:cime, cjms:cjme, ckms:ckme ) :: cfld
3449	INTEGER, DIMENSION ( nims:nime, njms:njme, nkms:nkme ) :: nfld
3450	REAL, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CBWGT1,CBWGT2,CBWGT3,CBWGT4 ! dummy
3451	REAL, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
3452	INTEGER, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CII,CJJ ! dummy
3453	INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: IIH,JJH
3454	INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
3455
3456	! local
3457
3458	LOGICAL FLIP
3459	INTEGER i,j,k,n
3460	REAL SUM,AMAXVAL
3461	REAL, DIMENSION (4, nims:nime, njms:njme ) :: NBWGT
3462
3463	!
3464	!*** INDEX CONVENTIONS
3465	!*** NBWGT4=0
3466	!*** 4
3467	!***
3468	!***
3469	!***
3470	!*** h
3471	!*** 1 2
3472	!*** NBWGT1=1 NBWGT2=0
3473	!***
3474	!***
3475	!*** 3
3476	!*** NBWGT3=0
3477
3478	DO J=NJTS,MIN(NJTE,NJDE-1)
3479	DO I=NITS,MIN(NITE,NIDE-1)
3480	IF(IMASK(I,J) .NE. 1)THEN
3481	NBWGT(1,I,J)=HBWGT1(I,J)
3482	NBWGT(2,I,J)=HBWGT2(I,J)
3483	NBWGT(3,I,J)=HBWGT3(I,J)
3484	NBWGT(4,I,J)=HBWGT4(I,J)
3485	ENDIF
3486	ENDDO
3487	ENDDO
3488
3489	DO J=NJTS,MIN(NJTE,NJDE-1)
3490	DO I=NITS,MIN(NITE,NIDE-1)
3491	IF(IMASK(I,J) .NE. 1)THEN
3492
3493	AMAXVAL=0.
3494	DO N=1,4
3495	AMAXVAL=amax1(NBWGT(N,I,J),AMAXVAL)
3496	ENDDO
3497
3498	FLIP=.TRUE.
3499	SUM=0.0
3500	DO N=1,4
3501	IF(AMAXVAL .EQ. NBWGT(N,I,J) .AND. FLIP)THEN
3502	NBWGT(N,I,J)=1.0
3503	FLIP=.FALSE.
3504	ELSE
3505	NBWGT(N,I,J)=0.0
3506	ENDIF
3507	SUM=SUM+NBWGT(N,I,J)
3508	IF(SUM .GT. 1.0)CALL wrf_error_fatal ( "horizontal interp error - interp_hnear_nmm" )
3509	ENDDO
3510	!
3511	ENDIF
3512	ENDDO
3513	ENDDO
3514
3515	DO J=NJTS,MIN(NJTE,NJDE-1)
3516	DO K=NKTS,NKTS
3517	DO I=NITS,MIN(NITE,NIDE-1)
3518	IF(IMASK(I,J) .NE. 1)THEN
3519	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7
3520	NFLD(I,J,K) = NBWGT(1,I,J)*CFLD(IIH(I,J), JJH(I,J) ,K) &
3521	+ NBWGT(2,I,J)*CFLD(IIH(I,J)+1,JJH(I,J) ,K) &
3522	+ NBWGT(3,I,J)*CFLD(IIH(I,J), JJH(I,J)-1,K) &
3523	+ NBWGT(4,I,J)*CFLD(IIH(I,J), JJH(I,J)+1,K)
3524	ELSE
3525	NFLD(I,J,K) = NBWGT(1,I,J)*CFLD(IIH(I,J), JJH(I,J) ,K) &
3526	+ NBWGT(2,I,J)*CFLD(IIH(I,J)+1,JJH(I,J) ,K) &
3527	+ NBWGT(3,I,J)*CFLD(IIH(I,J)+1,JJH(I,J)-1,K) &
3528	+ NBWGT(4,I,J)*CFLD(IIH(I,J)+1,JJH(I,J)+1,K)
3529	ENDIF
3530	ENDIF
3531	ENDDO
3532	ENDDO
3533	ENDDO
3534
3535	END SUBROUTINE interp_int_hnear_nmm
3536	!
3537	!--------------------------------------------------------------------------------------
3538	!
3539	SUBROUTINE nmm_bdy_hinterp (cfld, & ! CD field
3540	cids, cide, ckds, ckde, cjds, cjde, &
3541	cims, cime, ckms, ckme, cjms, cjme, &
3542	cits, cite, ckts, ckte, cjts, cjte, &
3543	nfld, & ! ND field
3544	nids, nide, nkds, nkde, njds, njde, &
3545	nims, nime, nkms, nkme, njms, njme, &
3546	nits, nite, nkts, nkte, njts, njte, &
3547	shw, & ! stencil half width
3548	imask, & ! interpolation mask
3549	xstag, ystag, & ! staggering of field
3550	ipos, jpos, & ! Position of lower left of nest in CD
3551	nri, nrj, & ! nest ratios
3552	c_bxs,n_bxs, &
3553	c_bxe,n_bxe, &
3554	c_bys,n_bys, &
3555	c_bye,n_bye, &
3556	c_btxs,n_btxs, &
3557	c_btxe,n_btxe, &
3558	c_btys,n_btys, &
3559	c_btye,n_btye, &
3560	CTEMP_B,NTEMP_B, & ! These temp arrays should be removed
3561	CTEMP_BT,NTEMP_BT, & ! later on
3562	CII, IIH, CJJ, JJH, CBWGT1, HBWGT1, & ! south-western grid locs and weights
3563	CBWGT2, HBWGT2, CBWGT3, HBWGT3, & ! note that "C"ourse grid ones are
3564	CBWGT4, HBWGT4 ) ! dummys
3565
3566	! use module_state_description
3567	USE module_configure
3568	USE module_wrf_error
3569
3570	IMPLICIT NONE
3571
3572
3573	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3574	cims, cime, ckms, ckme, cjms, cjme, &
3575	cits, cite, ckts, ckte, cjts, cjte, &
3576	nids, nide, nkds, nkde, njds, njde, &
3577	nims, nime, nkms, nkme, njms, njme, &
3578	nits, nite, nkts, nkte, njts, njte, &
3579	shw, &
3580	ipos, jpos, &
3581	nri, nrj
3582
3583	LOGICAL, INTENT(IN) :: xstag, ystag
3584
3585	REAL, DIMENSION ( cims:cime, cjms:cjme, ckms:ckme ) :: cfld
3586	REAL, DIMENSION ( nims:nime, njms:njme, nkms:nkme ) :: nfld
3587	!
3588	REAL, DIMENSION ( cims:cime, cjms:cjme, ckms:ckme ) :: ctemp_b,ctemp_bt
3589	REAL, DIMENSION ( nims:nime, njms:njme, nkms:nkme ) :: ntemp_b,ntemp_bt
3590	!
3591	INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
3592	REAL, DIMENSION( * ), INTENT(INOUT) :: c_bxs,n_bxs,c_bxe,n_bxe,c_bys,n_bys,c_bye,n_bye
3593	REAL, DIMENSION( * ), INTENT(INOUT) :: c_btxs,n_btxs,c_btxe,n_btxe,c_btys,n_btys,c_btye,n_btye
3594	REAL, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CBWGT1,CBWGT2,CBWGT3,CBWGT4 ! dummy
3595	REAL, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
3596	INTEGER, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CII,CJJ ! dummy
3597	INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: IIH,JJH
3598
3599	! Local
3600
3601	INTEGER :: i,j,k
3602	REAL, DIMENSION ( nims:nime, njms:njme, nkms:nkme ) :: cwk1,cwk2,cwk3,cwk4
3603
3604	! X start boundary
3605
3606	NMM_XS: IF(NITS .EQ. NIDS)THEN
3607	! WRITE(0,*)'ENTERING X1 START BOUNDARY AT MASS POINTS',NJTS,MIN(NJTE,NJDE-1)
3608	I = NIDS
3609	DO K = NKDS,NKDE
3610	DO J = NJTS,MIN(NJTE,NJDE-1)
3611	IF(MOD(J,2) .NE.0)THEN ! 1,3,5,7 of nested domain
3612	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7 of the parent domain
3613	CWK1(I,J,K) = HBWGT1(I,J)*CFLD(IIH(I,J), JJH(I,J) ,K) &
3614	+ HBWGT2(I,J)*CFLD(IIH(I,J)+1,JJH(I,J) ,K) &
3615	+ HBWGT3(I,J)*CFLD(IIH(I,J), JJH(I,J)-1,K) &
3616	+ HBWGT4(I,J)*CFLD(IIH(I,J), JJH(I,J)+1,K)
3617
3618
3619	ELSE
3620	CWK1(I,J,K) = HBWGT1(I,J)*CFLD(IIH(I,J), JJH(I,J) ,K) &
3621	+ HBWGT2(I,J)*CFLD(IIH(I,J)+1,JJH(I,J) ,K) &
3622	+ HBWGT3(I,J)*CFLD(IIH(I,J)+1,JJH(I,J)-1,K) &
3623	+ HBWGT4(I,J)*CFLD(IIH(I,J)+1,JJH(I,J)+1,K)
3624	ENDIF
3625	ELSE
3626	CWK1(I,J,K) = 0.0 ! even rows at mass points of the nested domain
3627	ENDIF
3628	ntemp_b(i,J,K) = CWK1(I,J,K)
3629	ntemp_bt(i,J,K) = 0.0
3630	END DO
3631	END DO
3632	ENDIF NMM_XS
3633
3634	! X end boundary
3635
3636	NMM_XE: IF(NITE-1 .EQ. NIDE-1)THEN
3637	! WRITE(0,*)'ENTERING X END BOUNDARY AT MASS POINTS',NJTS,MIN(NJTE,NJDE-1)
3638	I = NIDE-1
3639	DO K = NKDS,NKDE
3640	DO J = NJTS,MIN(NJTE,NJDE-1)
3641	IF(MOD(J,2) .NE.0)THEN ! 1,3,5,7 of the nested domain
3642	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7 of the parent domain
3643	CWK2(I,J,K) = HBWGT1(I,J)*CFLD(IIH(I,J), JJH(I,J) ,K) &
3644	+ HBWGT2(I,J)*CFLD(IIH(I,J)+1,JJH(I,J) ,K) &
3645	+ HBWGT3(I,J)*CFLD(IIH(I,J), JJH(I,J)-1,K) &
3646	+ HBWGT4(I,J)*CFLD(IIH(I,J), JJH(I,J)+1,K)
3647	ELSE
3648	CWK2(I,J,K) = HBWGT1(I,J)*CFLD(IIH(I,J), JJH(I,J) ,K) &
3649	+ HBWGT2(I,J)*CFLD(IIH(I,J)+1,JJH(I,J) ,K) &
3650	+ HBWGT3(I,J)*CFLD(IIH(I,J)+1,JJH(I,J)-1,K) &
3651	+ HBWGT4(I,J)*CFLD(IIH(I,J)+1,JJH(I,J)+1,K)
3652	ENDIF
3653	ELSE
3654	CWK2(I,J,K) = 0.0 ! even rows at mass points
3655	ENDIF
3656	ntemp_b(i,J,K) = CWK2(I,J,K)
3657	ntemp_bt(i,J,K) = 0.0
3658	END DO
3659	END DO
3660	ENDIF NMM_XE
3661
3662	! Y start boundary
3663
3664	NMM_YS: IF(NJTS .EQ. NJDS)THEN
3665	! WRITE(0,*)'ENTERING Y START BOUNDARY AT MASS POINTS',NITS,MIN(NITE,NIDE-1)
3666	J = NJDS
3667	DO K = NKDS, NKDE
3668	DO I = NITS,MIN(NITE,NIDE-1)
3669	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7
3670	CWK3(I,J,K) = HBWGT1(I,J)*CFLD(IIH(I,J), JJH(I,J) ,K) &
3671	+ HBWGT2(I,J)*CFLD(IIH(I,J)+1,JJH(I,J) ,K) &
3672	+ HBWGT3(I,J)*CFLD(IIH(I,J), JJH(I,J)-1,K) &
3673	+ HBWGT4(I,J)*CFLD(IIH(I,J), JJH(I,J)+1,K)
3674	ELSE
3675	CWK3(I,J,K) = HBWGT1(I,J)*CFLD(IIH(I,J), JJH(I,J) ,K) &
3676	+ HBWGT2(I,J)*CFLD(IIH(I,J)+1,JJH(I,J) ,K) &
3677	+ HBWGT3(I,J)*CFLD(IIH(I,J)+1,JJH(I,J)-1,K) &
3678	+ HBWGT4(I,J)*CFLD(IIH(I,J)+1,JJH(I,J)+1,K)
3679	ENDIF
3680	ntemp_b(i,J,K) = CWK3(I,J,K)
3681	ntemp_bt(i,J,K) = 0.0
3682	END DO
3683	END DO
3684	END IF NMM_YS
3685
3686	! Y end boundary
3687
3688	NMM_YE: IF(NJTE-1 .EQ. NJDE-1)THEN
3689	! WRITE(0,*)'ENTERING Y END BOUNDARY AT MASS POINTS',NITS,MIN(NITE,NIDE-1)
3690	J = NJDE-1
3691	DO K = NKDS,NKDE
3692	DO I = NITS,MIN(NITE,NIDE-1)
3693	IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7
3694	CWK4(I,J,K) = HBWGT1(I,J)*CFLD(IIH(I,J), JJH(I,J) ,K) &
3695	+ HBWGT2(I,J)*CFLD(IIH(I,J)+1,JJH(I,J) ,K) &
3696	+ HBWGT3(I,J)*CFLD(IIH(I,J), JJH(I,J)-1,K) &
3697	+ HBWGT4(I,J)*CFLD(IIH(I,J), JJH(I,J)+1,K)
3698	ELSE
3699	CWK4(I,J,K) = HBWGT1(I,J)*CFLD(IIH(I,J), JJH(I,J) ,K) &
3700	+ HBWGT2(I,J)*CFLD(IIH(I,J)+1,JJH(I,J) ,K) &
3701	+ HBWGT3(I,J)*CFLD(IIH(I,J)+1,JJH(I,J)-1,K) &
3702	+ HBWGT4(I,J)*CFLD(IIH(I,J)+1,JJH(I,J)+1,K)
3703
3704	ENDIF
3705	ntemp_b(i,J,K) = CWK4(I,J,K)
3706	ntemp_bt(i,J,K) = 0.0
3707	END DO
3708	END DO
3709	END IF NMM_YE
3710
3711	RETURN
3712
3713	END SUBROUTINE nmm_bdy_hinterp
3714
3715	!--------------------------------------------------------------------------------------
3716
3717	SUBROUTINE nmm_bdy_vinterp ( cfld, & ! CD field
3718	cids, cide, ckds, ckde, cjds, cjde, &
3719	cims, cime, ckms, ckme, cjms, cjme, &
3720	cits, cite, ckts, ckte, cjts, cjte, &
3721	nfld, & ! ND field
3722	nids, nide, nkds, nkde, njds, njde, &
3723	nims, nime, nkms, nkme, njms, njme, &
3724	nits, nite, nkts, nkte, njts, njte, &
3725	shw, & ! stencil half width
3726	imask, & ! interpolation mask
3727	xstag, ystag, & ! staggering of field
3728	ipos, jpos, & ! Position of lower left of nest in CD
3729	nri, nrj, & ! nest ratios
3730	c_bxs,n_bxs, &
3731	c_bxe,n_bxe, &
3732	c_bys,n_bys, &
3733	c_bye,n_bye, &
3734	c_btxs,n_btxs, &
3735	c_btxe,n_btxe, &
3736	c_btys,n_btys, &
3737	c_btye,n_btye, &
3738	CTEMP_B,NTEMP_B, & ! These temp arrays should be removed
3739	CTEMP_BT,NTEMP_BT, & ! later on
3740	CII, IIV, CJJ, JJV, CBWGT1, VBWGT1, & ! south-western grid locs and weights
3741	CBWGT2, VBWGT2, CBWGT3, VBWGT3, & ! note that "C"ourse grid ones are
3742	CBWGT4, VBWGT4 ) ! dummys
3743
3744	! use module_state_description
3745	USE module_configure
3746	USE module_wrf_error
3747
3748	IMPLICIT NONE
3749
3750
3751	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3752	cims, cime, ckms, ckme, cjms, cjme, &
3753	cits, cite, ckts, ckte, cjts, cjte, &
3754	nids, nide, nkds, nkde, njds, njde, &
3755	nims, nime, nkms, nkme, njms, njme, &
3756	nits, nite, nkts, nkte, njts, njte, &
3757	shw, &
3758	ipos, jpos, &
3759	nri, nrj
3760
3761	LOGICAL, INTENT(IN) :: xstag, ystag
3762
3763	REAL, DIMENSION ( cims:cime, cjms:cjme, ckms:ckme ) :: cfld
3764	REAL, DIMENSION ( nims:nime, njms:njme, nkms:nkme ) :: nfld
3765	!
3766	REAL, DIMENSION ( cims:cime, cjms:cjme, ckms:ckme ) :: ctemp_b,ctemp_bt
3767	REAL, DIMENSION ( nims:nime, njms:njme, nkms:nkme ) :: ntemp_b,ntemp_bt
3768	!
3769	INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
3770	REAL, DIMENSION( * ), INTENT(INOUT) :: c_bxs,n_bxs,c_bxe,n_bxe,c_bys,n_bys,c_bye,n_bye
3771	REAL, DIMENSION( * ), INTENT(INOUT) :: c_btxs,n_btxs,c_btxe,n_btxe,c_btys,n_btys,c_btye,n_btye
3772	REAL, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CBWGT1,CBWGT2,CBWGT3,CBWGT4 ! dummy
3773	REAL, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
3774	INTEGER, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CII,CJJ ! dummy
3775	INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: IIV,JJV
3776
3777	! Local
3778
3779	INTEGER :: i,j,k
3780	REAL, DIMENSION ( nims:nime, njms:njme, nkms:nkme ) :: cwk1,cwk2,cwk3,cwk4
3781
3782	! X start boundary
3783
3784	NMM_XS: IF(NITS .EQ. NIDS)THEN
3785	! WRITE(0,*)'ENTERING X START BOUNDARY AT VELOCITY POINTS',NITS,NIDS,NJTS,MIN(NJTE,NJDE-1)
3786	I = NIDS
3787	DO K = NKDS,NKDE
3788	DO J = NJTS,MIN(NJTE,NJDE-1)
3789	IF(MOD(J,2) .EQ.0)THEN ! 1,3,5,7 of nested domain
3790	IF(MOD(JJV(I,J),2) .NE. 0)THEN ! 1,3,5,7 of the parent domain
3791	CWK1(I,J,K) = VBWGT1(I,J)*CFLD(IIV(I,J), JJV(I,J) ,K) &
3792	+ VBWGT2(I,J)*CFLD(IIV(I,J)+1,JJV(I,J) ,K) &
3793	+ VBWGT3(I,J)*CFLD(IIV(I,J)+1,JJV(I,J)-1,K) &
3794	+ VBWGT4(I,J)*CFLD(IIV(I,J)+1,JJV(I,J)+1,K)
3795	ELSE
3796	CWK1(I,J,K) = VBWGT1(I,J)*CFLD(IIV(I,J), JJV(I,J) ,K) &
3797	+ VBWGT2(I,J)*CFLD(IIV(I,J)+1,JJV(I,J) ,K) &
3798	+ VBWGT3(I,J)*CFLD(IIV(I,J), JJV(I,J)-1,K) &
3799	+ VBWGT4(I,J)*CFLD(IIV(I,J), JJV(I,J)+1,K)
3800	ENDIF
3801	ELSE
3802	CWK1(I,J,K) = 0.0 ! odd rows along J, at mass points have zero velocity
3803	ENDIF
3804	ntemp_b(i,J,K) = CWK1(I,J,K)
3805	ntemp_bt(i,J,K) = 0.0
3806	END DO
3807	END DO
3808	ENDIF NMM_XS
3809
3810	! X end boundary
3811
3812	NMM_XE: IF(NITE-1 .EQ. NIDE-1)THEN
3813	! WRITE(0,*)'ENTERING X END BOUNDARY AT VELOCITY POINTS',NITE-1,NIDE-1,NJTS,MIN(NJTE,NJDE-1)
3814	I = NIDE-1
3815	DO K = NKDS,NKDE
3816	DO J = NJTS,MIN(NJTE,NJDE-1)
3817	IF(MOD(J,2) .EQ.0)THEN ! 1,3,5,7 of the nested domain
3818	IF(MOD(JJV(I,J),2) .NE. 0)THEN ! 1,3,5,7 of the parent domain
3819	CWK2(I,J,K) = VBWGT1(I,J)*CFLD(IIV(I,J), JJV(I,J) ,K) &
3820	+ VBWGT2(I,J)*CFLD(IIV(I,J)+1,JJV(I,J) ,K) &
3821	+ VBWGT3(I,J)*CFLD(IIV(I,J)+1,JJV(I,J)-1,K) &
3822	+ VBWGT4(I,J)*CFLD(IIV(I,J)+1,JJV(I,J)+1,K)
3823	ELSE
3824	CWK2(I,J,K) = VBWGT1(I,J)*CFLD(IIV(I,J), JJV(I,J) ,K) &
3825	+ VBWGT2(I,J)*CFLD(IIV(I,J)+1,JJV(I,J) ,K) &
3826	+ VBWGT3(I,J)*CFLD(IIV(I,J), JJV(I,J)-1,K) &
3827	+ VBWGT4(I,J)*CFLD(IIV(I,J), JJV(I,J)+1,K)
3828	ENDIF
3829	ELSE
3830	CWK2(I,J,K) = 0.0 ! odd rows at mass points
3831	ENDIF
3832	ntemp_b(i,J,K) = CWK2(I,J,K)
3833	ntemp_bt(i,J,K) = 0.0
3834	END DO
3835	END DO
3836	ENDIF NMM_XE
3837
3838	! Y start boundary
3839
3840	NMM_YS: IF(NJTS .EQ. NJDS)THEN
3841	! WRITE(0,*)'ENTERING Y START BOUNDARY AT VELOCITY POINTS',NJTS,NJDS,NITS,MIN(NITE,NIDE-1)
3842	J = NJDS
3843	DO K = NKDS, NKDE
3844	DO I = NITS,MIN(NITE,NIDE-2) ! NIDE-1 SHOULD NOT MATTER IF WE FILL UP PHANTOM CELL
3845	IF(MOD(JJV(I,J),2) .NE. 0)THEN ! 1,3,5,7
3846	CWK3(I,J,K) = VBWGT1(I,J)*CFLD(IIV(I,J), JJV(I,J) ,K) &
3847	+ VBWGT2(I,J)*CFLD(IIV(I,J)+1,JJV(I,J) ,K) &
3848	+ VBWGT3(I,J)*CFLD(IIV(I,J)+1,JJV(I,J)-1,K) &
3849	+ VBWGT4(I,J)*CFLD(IIV(I,J)+1,JJV(I,J)+1,K)
3850	ELSE
3851	CWK3(I,J,K) = VBWGT1(I,J)*CFLD(IIV(I,J), JJV(I,J) ,K) &
3852	+ VBWGT2(I,J)*CFLD(IIV(I,J)+1,JJV(I,J) ,K) &
3853	+ VBWGT3(I,J)*CFLD(IIV(I,J), JJV(I,J)-1,K) &
3854	+ VBWGT4(I,J)*CFLD(IIV(I,J), JJV(I,J)+1,K)
3855	ENDIF
3856	ntemp_b(i,J,K) = CWK3(I,J,K)
3857	ntemp_bt(i,J,K) = 0.0
3858	END DO
3859	END DO
3860	END IF NMM_YS
3861
3862	! Y end boundary
3863
3864	NMM_YE: IF(NJTE-1 .EQ. NJDE-1)THEN
3865	! WRITE(0,*)'ENTERING Y END BOUNDARY AT VELOCITY POINTS',NJTE-1,NJDE-1,NITS,MIN(NITE,NIDE-1)
3866	J = NJDE-1
3867	DO K = NKDS,NKDE
3868	DO I = NITS,MIN(NITE,NIDE-2) ! NIDE-1 SHOULD NOT MATTER IF WE FILL UP PHANTOM CELL
3869	IF(MOD(JJV(I,J),2) .NE. 0)THEN ! 1,3,5,7
3870	CWK4(I,J,K) = VBWGT1(I,J)*CFLD(IIV(I,J), JJV(I,J) ,K) &
3871	+ VBWGT2(I,J)*CFLD(IIV(I,J)+1,JJV(I,J) ,K) &
3872	+ VBWGT3(I,J)*CFLD(IIV(I,J)+1,JJV(I,J)-1,K) &
3873	+ VBWGT4(I,J)*CFLD(IIV(I,J)+1,JJV(I,J)+1,K)
3874	ELSE
3875	CWK4(I,J,K) = VBWGT1(I,J)*CFLD(IIV(I,J), JJV(I,J) ,K) &
3876	+ VBWGT2(I,J)*CFLD(IIV(I,J)+1,JJV(I,J) ,K) &
3877	+ VBWGT3(I,J)*CFLD(IIV(I,J), JJV(I,J)-1,K) &
3878	+ VBWGT4(I,J)*CFLD(IIV(I,J), JJV(I,J)+1,K)
3879	ENDIF
3880	ntemp_b(i,J,K) = CWK4(I,J,K)
3881	ntemp_bt(i,J,K) = 0.0
3882	END DO
3883	END DO
3884	END IF NMM_YE
3885
3886	RETURN
3887
3888	END SUBROUTINE nmm_bdy_vinterp
3889
3890	!
3891	!=======================================================================================
3892	! E grid interpolation: simple copy from parent to mother domain
3893	!=======================================================================================
3894	!
3895
3896	SUBROUTINE nmm_copy ( cfld, & ! CD field
3897	cids, cide, ckds, ckde, cjds, cjde, &
3898	cims, cime, ckms, ckme, cjms, cjme, &
3899	cits, cite, ckts, ckte, cjts, cjte, &
3900	nfld, & ! ND field
3901	nids, nide, nkds, nkde, njds, njde, &
3902	nims, nime, nkms, nkme, njms, njme, &
3903	nits, nite, nkts, nkte, njts, njte, &
3904	shw, & ! stencil half width
3905	imask, & ! interpolation mask
3906	xstag, ystag, & ! staggering of field
3907	ipos, jpos, & ! Position of lower left of nest in CD
3908	nri, nrj, & ! nest ratios
3909	CII, IIH, CJJ, JJH )
3910
3911	USE module_timing
3912	IMPLICIT NONE
3913
3914	LOGICAL, INTENT(IN) :: xstag, ystag
3915	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3916	cims, cime, ckms, ckme, cjms, cjme, &
3917	cits, cite, ckts, ckte, cjts, cjte, &
3918	nids, nide, nkds, nkde, njds, njde, &
3919	nims, nime, nkms, nkme, njms, njme, &
3920	nits, nite, nkts, nkte, njts, njte, &
3921	shw, &
3922	ipos, jpos, &
3923	nri, nrj
3924	REAL, DIMENSION ( cims:cime, cjms:cjme, ckms:ckme ), INTENT(IN) :: cfld
3925	REAL, DIMENSION ( nims:nime, njms:njme, nkms:nkme ), INTENT(INOUT) :: nfld
3926	INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: imask
3927	INTEGER, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CII,CJJ ! dummy
3928	INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: IIH,JJH
3929
3930	! local
3931	INTEGER i,j,k
3932
3933	DO J=NJTS,MIN(NJTE,NJDE-1)
3934	DO K=NKTS,NKTE
3935	DO I=NITS,MIN(NITE,NIDE-1)
3936	NFLD(I,J,K) = CFLD(IIH(I,J),JJH(I,J),K)
3937	ENDDO
3938	ENDDO
3939	ENDDO
3940
3941	RETURN
3942
3943	END SUBROUTINE nmm_copy
3944	!
3945	!=======================================================================================
3946	! E grid test for mass point coincidence
3947	!=======================================================================================
3948	!
3949	SUBROUTINE test_nmm (cfld, & ! CD field
3950	cids, cide, ckds, ckde, cjds, cjde, &
3951	cims, cime, ckms, ckme, cjms, cjme, &
3952	cits, cite, ckts, ckte, cjts, cjte, &
3953	nfld, & ! ND field
3954	nids, nide, nkds, nkde, njds, njde, &
3955	nims, nime, nkms, nkme, njms, njme, &
3956	nits, nite, nkts, nkte, njts, njte, &
3957	shw, & ! stencil half width for interp
3958	imask, & ! interpolation mask
3959	xstag, ystag, & ! staggering of field
3960	ipos, jpos, & ! Position of lower left of nest in CD
3961	nri, nrj, & ! nest ratios
3962	CII, IIH, CJJ, JJH, CBWGT1, HBWGT1, & ! south-western grid locs and weights
3963	CBWGT2, HBWGT2, CBWGT3, HBWGT3, & ! note that "C"ourse grid ones are
3964	CBWGT4, HBWGT4 ) ! dummys for weights
3965	USE module_timing
3966	IMPLICIT NONE
3967
3968	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
3969	cims, cime, ckms, ckme, cjms, cjme, &
3970	cits, cite, ckts, ckte, cjts, cjte, &
3971	nids, nide, nkds, nkde, njds, njde, &
3972	nims, nime, nkms, nkme, njms, njme, &
3973	nits, nite, nkts, nkte, njts, njte, &
3974	shw, &
3975	ipos, jpos, &
3976	nri, nrj
3977	LOGICAL, INTENT(IN) :: xstag, ystag
3978
3979	REAL, DIMENSION ( cims:cime, cjms:cjme, ckms:ckme ) :: cfld
3980	REAL, DIMENSION ( nims:nime, njms:njme, nkms:nkme ) :: nfld
3981	REAL, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CBWGT1,CBWGT2,CBWGT3,CBWGT4 ! dummy
3982	REAL, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
3983	INTEGER, DIMENSION ( cims:cime, cjms:cjme ), INTENT(IN) :: CII,CJJ ! dummy
3984	INTEGER, DIMENSION ( nims:nime, njms:njme ), INTENT(IN) :: IIH,JJH
3985	INTEGER, DIMENSION ( nims:nime, njms:njme ) :: imask
3986
3987	! local
3988	INTEGER i,j,k
3989	REAL,PARAMETER :: error=0.0001,error1=1.0
3990	REAL :: diff
3991	!
3992	!*** CHECK DOMAIN BOUNDS BEFORE INTERPOLATION
3993	!
3994	DO J=NJTS,MIN(NJTE,NJDE-1)
3995	DO I=NITS,MIN(NITE,NIDE-1)
3996	IF(IIH(i,j).LT.(CIDS-shw) .OR. IIH(i,j).GT.(CIDE+shw)) &
3997	CALL wrf_error_fatal ('hpoints:check domain bounds along x' )
3998	IF(JJH(i,j).LT.(CJDS-shw) .OR. JJH(i,j).GT.(CJDE+shw)) &
3999	CALL wrf_error_fatal ('hpoints:check domain bounds along y' )
4000	ENDDO
4001	ENDDO
4002
4003	!
4004	!*** INDEX CONVENTIONS
4005	!*** HBWGT4
4006	!*** 4
4007	!***
4008	!***
4009	!***
4010	!*** h
4011	!*** 1 2
4012	!*** HBWGT1 HBWGT2
4013	!***
4014	!***
4015	!*** 3
4016	!*** HBWGT3
4017
4018
4019	! WRITE(0,*)NITS,MIN(NITE,NIDE-1),CITS,CITE
4020	DO J=NJTS,MIN(NJTE,NJDE-1)
4021	DO K=NKDS,NKDE
4022	DO I=NITS,MIN(NITE,NIDE-1)
4023	IF(ABS(1.0-HBWGT1(I,J)) .LE. ERROR)THEN
4024	DIFF=ABS(NFLD(I,J,K)-CFLD(IIH(I,J),JJH(I,J),K))
4025	IF(DIFF .GT. ERROR)THEN
4026	CALL wrf_debug(1,"dyn_nmm: NON-COINCIDENT, NESTED MASS POINT")
4027	WRITE(0,*)I,IIH(I,J),J,JJH(I,J),HBWGT1(I,J),NFLD(I,J,K),CFLD(IIH(I,J),JJH(I,J),K),DIFF
4028	ENDIF
4029	IF(DIFF .GT. ERROR1)THEN
4030	WRITE(0,*)I,IIH(I,J),J,JJH(I,J),HBWGT1(I,J),NFLD(I,J,K),CFLD(IIH(I,J),JJH(I,J),K),DIFF
4031	CALL wrf_error_fatal ('dyn_nmm: NON-COINCIDENT, NESTED MASS POINT')
4032	ENDIF
4033	ENDIF
4034	ENDDO
4035	ENDDO
4036	ENDDO
4037
4038	END SUBROUTINE test_nmm
4039
4040	!==================================
4041	! this is the default function used in nmm feedback at mass points.
4042
4043	SUBROUTINE nmm_feedback ( cfld, & ! CD field
4044	cids, cide, ckds, ckde, cjds, cjde, &
4045	cims, cime, ckms, ckme, cjms, cjme, &
4046	cits, cite, ckts, ckte, cjts, cjte, &
4047	nfld, & ! ND field
4048	nids, nide, nkds, nkde, njds, njde, &
4049	nims, nime, nkms, nkme, njms, njme, &
4050	nits, nite, nkts, nkte, njts, njte, &
4051	shw, & ! stencil half width for interp
4052	imask, & ! interpolation mask
4053	xstag, ystag, & ! staggering of field
4054	ipos, jpos, & ! Position of lower left of nest in CD
4055	nri, nrj, & ! nest ratios
4056	CII, IIH, CJJ, JJH, &
4057	CBWGT1, HBWGT1, CBWGT2, HBWGT2, &
4058	CBWGT3, HBWGT3, CBWGT4, HBWGT4 )
4059	USE module_configure
4060	IMPLICIT NONE
4061
4062
4063	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
4064	cims, cime, ckms, ckme, cjms, cjme, &
4065	cits, cite, ckts, ckte, cjts, cjte, &
4066	nids, nide, nkds, nkde, njds, njde, &
4067	nims, nime, nkms, nkme, njms, njme, &
4068	nits, nite, nkts, nkte, njts, njte, &
4069	shw, &
4070	ipos, jpos, &
4071	nri, nrj
4072	INTEGER,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CII,CJJ ! dummy
4073	INTEGER,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: IIH,JJH
4074	REAL,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CBWGT1,CBWGT2,CBWGT3,CBWGT4
4075	REAL,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
4076	LOGICAL, INTENT(IN) :: xstag, ystag
4077
4078	REAL, DIMENSION ( cims:cime, cjms:cjme, ckms:ckme ), INTENT(OUT) :: cfld
4079	REAL, DIMENSION ( nims:nime, njms:njme, nkms:nkme ), INTENT(IN) :: nfld
4080	INTEGER, DIMENSION ( nims:nime, njms:njme ),INTENT(IN) :: imask
4081
4082	! Local
4083
4084	INTEGER ci, cj, ck, ni, nj, nk, ip, jp, ioff, joff, ioffa, joffa
4085	INTEGER :: icmin,icmax,jcmin,jcmax
4086	INTEGER :: is, ipoints,jpoints,ijpoints
4087	INTEGER , PARAMETER :: passes = 2
4088	REAL :: AVGH
4089
4090	!=====================================================================================
4091	!
4092
4093	IF(nri .ne. 3 .OR. nrj .ne. 3) &
4094	CALL wrf_error_fatal ('Feedback works for only 1:3 ratios, currently. Modify the namelist' )
4095
4096	! WRITE(0,*)'SIMPLE FEED BACK IS SWITCHED ON FOR MASS'
4097
4098	CFLD = 9999.0
4099
4100	DO ck = ckts, ckte
4101	nk = ck
4102	DO cj = MAX(jpos+1,cjts),MIN(jpos+(njde-njds)/nrj-1,cjte) ! exclude top and bottom BCs
4103	nj = (cj-jpos)*nrj + 1
4104	if(mod(cj,2) .eq. 0)THEN
4105	is=0 ! even rows for mass points (2,4,6,8)
4106	else
4107	is=1 ! odd rows for mass points (1,3,5,7)
4108	endif
4109	DO ci = MAX(ipos+is,cits),MIN(ipos+(nide-nids)/nri-1,cite) ! excludes LBCs
4110	ni = (ci-ipos)*nri + 2 -is
4111	IF(IS==0)THEN ! (2,4,6,8)
4112	! AVGH = NFLD(NI,NJ+1,NK) + NFLD(NI,NJ-1,NK) + NFLD(NI+1,NJ+1,NK)+ NFLD(NI+1,NJ-1,NK) &
4113	! + NFLD(NI+1,NJ,NK) + NFLD(NI-1,NJ,NK) + NFLD(NI,NJ+2,NK) + NFLD(NI,NJ-2,NK) &
4114	! + NFLD(NI+1,NJ+2,NK)+ NFLD(NI-1,NJ+2,NK)+ NFLD(NI+1,NJ-2,NK)+ NFLD(NI-1,NJ-2,NK)
4115
4116	AVGH = NFLD(NI,NJ+2,NK) &
4117	+ NFLD(NI ,NJ+1,NK) + NFLD(NI+1,NJ+1,NK) &
4118	+ NFLD(NI-1,NJ ,NK) + NFLD(NI,NJ ,NK) + NFLD(NI+1,NJ ,NK) &
4119	+ NFLD(NI ,NJ-1,NK) + NFLD(NI+1,NJ-1,NK) &
4120	+ NFLD(NI,NJ-2,NK)
4121
4122	ELSE
4123	! AVGH = NFLD(NI,NJ+1,NK) + NFLD(NI,NJ-1,NK) + NFLD(NI-1,NJ+1,NK)+ NFLD(NI-1,NJ-1,NK) &
4124	! + NFLD(NI+1,NJ,NK) + NFLD(NI-1,NJ,NK) + NFLD(NI,NJ+2,NK) + NFLD(NI,NJ-2,NK) &
4125	! + NFLD(NI+1,NJ+2,NK)+ NFLD(NI-1,NJ+2,NK)+ NFLD(NI+1,NJ-2,NK)+ NFLD(NI-1,NJ-2,NK)
4126
4127	AVGH = NFLD(NI,NJ+2,NK) &
4128	+ NFLD(NI-1,NJ+1,NK) + NFLD(NI,NJ+1,NK) &
4129	+ NFLD(NI-1,NJ ,NK) + NFLD(NI,NJ ,NK) + NFLD(NI+1,NJ ,NK) &
4130	+ NFLD(NI-1,NJ-1,NK) + NFLD(NI,NJ-1,NK) &
4131	+ NFLD(NI,NJ-2,NK)
4132
4133	ENDIF
4134	!dusan CFLD(CI,CK,CJ) = 0.5CFLD(CI,CK,CJ) + 0.5(NFLD(NI,NK,NJ)+AVGH)/13.0
4135	! CFLD(CI,CJ,CK) = (NFLD(NI,NJ,NK)+AVGH)/13.0
4136	CFLD(CI,CJ,CK) = AVGH/9.0
4137	ENDDO
4138	ENDDO
4139	ENDDO
4140
4141	END SUBROUTINE nmm_feedback
4142
4143	!===========================================================================================
4144
4145	SUBROUTINE nmm_vfeedback ( cfld, & ! CD field
4146	cids, cide, ckds, ckde, cjds, cjde, &
4147	cims, cime, ckms, ckme, cjms, cjme, &
4148	cits, cite, ckts, ckte, cjts, cjte, &
4149	nfld, & ! ND field
4150	nids, nide, nkds, nkde, njds, njde, &
4151	nims, nime, nkms, nkme, njms, njme, &
4152	nits, nite, nkts, nkte, njts, njte, &
4153	shw, & ! stencil half width for interp
4154	imask, & ! interpolation mask
4155	xstag, ystag, & ! staggering of field
4156	ipos, jpos, & ! Position of lower left of nest in CD
4157	nri, nrj, & ! nest ratios
4158	CII, IIV, CJJ, JJV, &
4159	CBWGT1, VBWGT1, CBWGT2, VBWGT2, &
4160	CBWGT3, VBWGT3, CBWGT4, VBWGT4 )
4161	USE module_configure
4162	IMPLICIT NONE
4163
4164
4165	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
4166	cims, cime, ckms, ckme, cjms, cjme, &
4167	cits, cite, ckts, ckte, cjts, cjte, &
4168	nids, nide, nkds, nkde, njds, njde, &
4169	nims, nime, nkms, nkme, njms, njme, &
4170	nits, nite, nkts, nkte, njts, njte, &
4171	shw, &
4172	ipos, jpos, &
4173	nri, nrj
4174	INTEGER,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CII,CJJ ! dummy
4175	INTEGER,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: IIV,JJV
4176	REAL,DIMENSION(cims:cime,cjms:cjme), INTENT(IN) :: CBWGT1,CBWGT2,CBWGT3,CBWGT4
4177	REAL,DIMENSION(nims:nime,njms:njme), INTENT(IN) :: VBWGT1,VBWGT2,VBWGT3,VBWGT4
4178	LOGICAL, INTENT(IN) :: xstag, ystag
4179
4180	REAL, DIMENSION ( cims:cime, cjms:cjme, ckms:ckme ), INTENT(OUT) :: cfld
4181	REAL, DIMENSION ( nims:nime, njms:njme, nkms:nkme ), INTENT(IN) :: nfld
4182	INTEGER, DIMENSION ( nims:nime, njms:njme ),INTENT(IN) :: imask
4183
4184	! Local
4185
4186	INTEGER ci, cj, ck, ni, nj, nk, ip, jp, ioff, joff, ioffa, joffa
4187	INTEGER :: icmin,icmax,jcmin,jcmax
4188	INTEGER :: is, ipoints,jpoints,ijpoints
4189	INTEGER , PARAMETER :: passes = 2
4190	REAL :: AVGV
4191
4192	!=====================================================================================
4193	!
4194
4195	IF(nri .ne. 3 .OR. nrj .ne. 3) &
4196	CALL wrf_error_fatal ('Feedback works for only 1:3 ratios, currently. Modify the namelist')
4197
4198	! WRITE(0,*)'SIMPLE FEED BACK IS SWITCHED ON FOR VELOCITY'
4199
4200	CFLD = 9999.0
4201
4202	DO ck = ckts, ckte
4203	nk = ck
4204	DO cj = MAX(jpos+1,cjts),MIN(jpos+(njde-njds)/nrj-1,cjte) ! exclude top and bottom BCs
4205	nj = (cj-jpos)*nrj + 1
4206	if(mod(cj,2) .eq. 0)THEN
4207	is=1 ! even rows for velocity points (2,4,6,8)
4208	else
4209	is=0 ! odd rows for velocity points (1,3,5,7)
4210	endif
4211	DO ci = MAX(ipos+is,cits),MIN(ipos+(nide-nids)/nri-1,cite) ! excludes LBCs
4212	ni = (ci-ipos)*nri + 2 -is
4213	IF(IS==0)THEN ! (1,3,5,7)
4214	! AVGV = NFLD(NI,NK,NJ+1) + NFLD(NI,NK,NJ-1) + NFLD(NI+1,NK,NJ+1)+ NFLD(NI+1,NK,NJ-1) &
4215	! + NFLD(NI+1,NK,NJ) + NFLD(NI-1,NK,NJ) + NFLD(NI,NK,NJ+2) + NFLD(NI,NK,NJ-2) &
4216	! + NFLD(NI+1,NK,NJ+2)+ NFLD(NI-1,NK,NJ+2)+ NFLD(NI+1,NK,NJ-2)+ NFLD(NI-1,NK,NJ-2)
4217
4218	AVGV = NFLD(NI,NJ+2,NK) &
4219	+ NFLD(NI ,NJ+1,NK) + NFLD(NI+1,NJ+1,NK) &
4220	+ NFLD(NI-1,NJ ,NK) + NFLD(NI,NJ ,NK) + NFLD(NI+1,NJ ,NK) &
4221	+ NFLD(NI ,NJ-1,NK) + NFLD(NI+1,NJ-1,NK) &
4222	+ NFLD(NI,NJ-2,NK)
4223
4224	ELSE
4225	! AVGV = NFLD(NI,NK,NJ+1) + NFLD(NI,NK,NJ-1) + NFLD(NI-1,NK,NJ+1)+ NFLD(NI-1,NK,NJ-1) &
4226	! + NFLD(NI+1,NK,NJ) + NFLD(NI-1,NK,NJ) + NFLD(NI,NK,NJ+2) + NFLD(NI,NK,NJ-2) &
4227	! + NFLD(NI+1,NK,NJ+2)+ NFLD(NI-1,NK,NJ+2)+ NFLD(NI+1,NK,NJ-2)+ NFLD(NI-1,NK,NJ-2)
4228
4229	AVGV = NFLD(NI,NJ+2,NK) &
4230	+ NFLD(NI-1,NJ+1,NK) + NFLD(NI,NJ+1,NK) &
4231	+ NFLD(NI-1,NJ ,NK) + NFLD(NI,NJ ,NK) + NFLD(NI+1,NJ ,NK) &
4232	+ NFLD(NI-1,NJ-1,NK) + NFLD(NI,NJ-1,NK) &
4233	+ NFLD(NI,NJ-2,NK)
4234
4235	ENDIF
4236	!dusan CFLD(CI,CK,CJ) = 0.5CFLD(CI,CK,CJ) + 0.5(NFLD(NI,NK,NJ)+AVGV)/13.0
4237	! CFLD(CI,CK,CJ) = (NFLD(NI,NK,NJ)+AVGV)/13.0
4238	CFLD(CI,CJ,CK) = AVGV/9.0
4239	ENDDO
4240	ENDDO
4241	ENDDO
4242
4243	END SUBROUTINE nmm_vfeedback
4244
4245
4246	SUBROUTINE nmm_smoother ( cfld , &
4247	cids, cide, ckds, ckde, cjds, cjde, &
4248	cims, cime, ckms, ckme, cjms, cjme, &
4249	cits, cite, ckts, ckte, cjts, cjte, &
4250	nids, nide, nkds, nkde, njds, njde, &
4251	nims, nime, nkms, nkme, njms, njme, &
4252	nits, nite, nkts, nkte, njts, njte, &
4253	xstag, ystag, &
4254	ipos, jpos, &
4255	nri, nrj &
4256	)
4257
4258	USE module_configure
4259	IMPLICIT NONE
4260
4261	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
4262	cims, cime, ckms, ckme, cjms, cjme, &
4263	cits, cite, ckts, ckte, cjts, cjte, &
4264	nids, nide, nkds, nkde, njds, njde, &
4265	nims, nime, nkms, nkme, njms, njme, &
4266	nits, nite, nkts, nkte, njts, njte, &
4267	nri, nrj, &
4268	ipos, jpos
4269	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ), INTENT(INOUT) :: cfld
4270	LOGICAL, INTENT(IN) :: xstag, ystag
4271
4272
4273	! Local
4274
4275	INTEGER :: feedback
4276	INTEGER, PARAMETER :: smooth_passes = 5
4277
4278	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfldnew
4279	INTEGER :: ci, cj, ck
4280	INTEGER :: is, npass
4281	REAL :: AVGH
4282
4283	RETURN
4284	! If there is no feedback, there can be no smoothing.
4285
4286	CALL nl_get_feedback ( 1, feedback )
4287	IF ( feedback == 0 ) RETURN
4288
4289	WRITE(0,*)'SIMPLE SMOOTHER IS SWITCHED ON FOR HEIGHT'
4290
4291	DO npass = 1, smooth_passes
4292
4293	DO cj = MAX(jpos+1,cjts),MIN(jpos+(njde-njds)/nrj-1,cjte) ! exclude top and bottom BCs
4294	if(mod(cj,2) .eq. 0)THEN
4295	is=0 ! even rows for mass points (2,4,6,8)
4296	else
4297	is=1 ! odd rows for mass points (1,3,5,7)
4298	endif
4299	DO ck = ckts, ckte
4300	DO ci = MAX(ipos+is,cits),MIN(ipos+(nide-nids)/nri-1,cite) ! excludes LBCs
4301	IF(IS==0)THEN ! (2,4,6,8)
4302	AVGH = CFLD(CI,CK,CJ+1) + CFLD(CI,CK,CJ-1) + CFLD(CI+1,CK,CJ+1) + CFLD(CI+1,CK,CJ-1)
4303	ELSE
4304	AVGH = CFLD(CI,CK,CJ+1) + CFLD(CI,CK,CJ-1) + CFLD(CI-1,CK,CJ+1) + CFLD(CI-1,CK,CJ-1)
4305	ENDIF
4306	CFLDNEW(CI,CK,CJ) = (AVGH + 4*CFLD(CI,CK,CJ)) / 8.0
4307	ENDDO
4308	ENDDO
4309	ENDDO
4310
4311	DO cj = MAX(jpos+1,cjts),MIN(jpos+(njde-njds)/nrj-1,cjte) ! exclude top and bottom BCs
4312	if(mod(cj,2) .eq. 0)THEN
4313	is=0 ! even rows for mass points (2,4,6,8)
4314	else
4315	is=1 ! odd rows for mass points (1,3,5,7)
4316	endif
4317	DO ck = ckts, ckte
4318	DO ci = MAX(ipos+is,cits),MIN(ipos+(nide-nids)/nri-1,cite) ! excludes LBCs
4319	CFLD(CI,CK,CJ) = CFLDNEW(CI,CK,CJ)
4320	ENDDO
4321	ENDDO
4322	ENDDO
4323
4324	ENDDO ! do npass
4325
4326	END SUBROUTINE nmm_smoother
4327
4328
4329	SUBROUTINE nmm_vsmoother ( cfld , &
4330	cids, cide, ckds, ckde, cjds, cjde, &
4331	cims, cime, ckms, ckme, cjms, cjme, &
4332	cits, cite, ckts, ckte, cjts, cjte, &
4333	nids, nide, nkds, nkde, njds, njde, &
4334	nims, nime, nkms, nkme, njms, njme, &
4335	nits, nite, nkts, nkte, njts, njte, &
4336	xstag, ystag, &
4337	ipos, jpos, &
4338	nri, nrj &
4339	)
4340
4341	USE module_configure
4342	IMPLICIT NONE
4343
4344	INTEGER, INTENT(IN) :: cids, cide, ckds, ckde, cjds, cjde, &
4345	cims, cime, ckms, ckme, cjms, cjme, &
4346	cits, cite, ckts, ckte, cjts, cjte, &
4347	nids, nide, nkds, nkde, njds, njde, &
4348	nims, nime, nkms, nkme, njms, njme, &
4349	nits, nite, nkts, nkte, njts, njte, &
4350	nri, nrj, &
4351	ipos, jpos
4352	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ), INTENT(INOUT) :: cfld
4353	LOGICAL, INTENT(IN) :: xstag, ystag
4354
4355
4356	! Local
4357
4358	INTEGER :: feedback
4359	INTEGER, PARAMETER :: smooth_passes = 5
4360
4361	REAL, DIMENSION ( cims:cime, ckms:ckme, cjms:cjme ) :: cfldnew
4362	INTEGER :: ci, cj, ck
4363	INTEGER :: is, npass
4364	REAL :: AVGV
4365
4366	RETURN
4367	! If there is no feedback, there can be no smoothing.
4368
4369	CALL nl_get_feedback ( 1, feedback )
4370	IF ( feedback == 0 ) RETURN
4371
4372	WRITE(0,*)'SIMPLE SMOOTHER IS SWITCHED ON FOR VELOCITY'
4373
4374	DO npass = 1, smooth_passes
4375
4376	DO cj = MAX(jpos+1,cjts),MIN(jpos+(njde-njds)/nrj-1,cjte) ! exclude top and bottom BCs
4377	if(mod(cj,2) .eq. 0)THEN
4378	is=1 ! even rows for mass points (2,4,6,8)
4379	else
4380	is=0 ! odd rows for mass points (1,3,5,7)
4381	endif
4382	DO ck = ckts, ckte
4383	DO ci = MAX(ipos+is,cits),MIN(ipos+(nide-nids)/nri-1,cite) ! excludes LBCs
4384	IF(IS==0)THEN ! (2,4,6,8)
4385	AVGV = CFLD(CI,CK,CJ+1) + CFLD(CI,CK,CJ-1) + CFLD(CI+1,CK,CJ+1) + CFLD(CI+1,CK,CJ-1)
4386	ELSE
4387	AVGV = CFLD(CI,CK,CJ+1) + CFLD(CI,CK,CJ-1) + CFLD(CI-1,CK,CJ+1) + CFLD(CI-1,CK,CJ-1)
4388	ENDIF
4389	CFLDNEW(CI,CK,CJ) = (AVGV + 4*CFLD(CI,CK,CJ)) / 8.0
4390	ENDDO
4391	ENDDO
4392	ENDDO
4393
4394	DO cj = MAX(jpos+1,cjts),MIN(jpos+(njde-njds)/nrj-1,cjte) ! exclude top and bottom BCs
4395	if(mod(cj,2) .eq. 0)THEN
4396	is=1 ! even rows for mass points (2,4,6,8)
4397	else
4398	is=0 ! odd rows for mass points (1,3,5,7)
4399	endif
4400	DO ck = ckts, ckte
4401	DO ci = MAX(ipos+is,cits),MIN(ipos+(nide-nids)/nri-1,cite) ! excludes LBCs
4402	CFLD(CI,CK,CJ) = CFLDNEW(CI,CK,CJ)
4403	ENDDO
4404	ENDDO
4405	ENDDO
4406
4407	ENDDO
4408
4409	END SUBROUTINE nmm_vsmoother
4410	!======================================================================================
4411	! End of gopal's doing
4412	!======================================================================================
4413	#endif

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format