Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

interp_fcn.F @ 2754

Last change on this file since 2754 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 235.2 KB

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

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

Download in other formats:

Original Format