interp_fcn.F

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