module_bc.F

References to this file elsewhere.
1 !WRF:MODEL_LAYER:BOUNDARY
2 !
3 
4 MODULE module_bc
5 
6    USE module_configure
7    USE module_wrf_error
8    IMPLICIT NONE
9 
10 !   TYPE bcs
11 !
12 !     LOGICAL                     :: periodic_x
13 !     LOGICAL                     :: symmetric_xs
14 !     LOGICAL                     :: symmetric_xe
15 !     LOGICAL                     :: open_xs
16 !     LOGICAL                     :: open_xe
17 !     LOGICAL                     :: periodic_y
18 !     LOGICAL                     :: symmetric_ys
19 !     LOGICAL                     :: symmetric_ye
20 !     LOGICAL                     :: open_ys
21 !     LOGICAL                     :: open_ye
22 !     LOGICAL                     :: nested
23 !     LOGICAL                     :: specified
24 !     LOGICAL                     :: top_radiation
25 !
26 !   END TYPE bcs
27 
28 !  set the bdyzone.  We are hardwiring this here and we'll
29 !  decide later where it should be set and stored
30 
31    INTEGER, PARAMETER            :: bdyzone = 4
32    INTEGER, PARAMETER            :: bdyzone_x = bdyzone
33    INTEGER, PARAMETER            :: bdyzone_y = bdyzone
34 
35    INTERFACE stuff_bdy
36      MODULE PROCEDURE stuff_bdy_new , stuff_bdy_old
37    END INTERFACE
38 
39    INTERFACE stuff_bdytend
40      MODULE PROCEDURE stuff_bdytend_new , stuff_bdytend_old
41    END INTERFACE
42 
43 CONTAINS
44 
45   SUBROUTINE boundary_condition_check ( config_flags, bzone, error, gn )
46 
47 !  this routine checks the boundary condition logicals 
48 !  to make sure that the boundary conditions are not over
49 !  or under specified.  The routine also checks that the
50 !  boundary zone is sufficiently sized for the specified
51 !  boundary conditions
52 
53   IMPLICIT NONE
54 
55   TYPE( grid_config_rec_type ) config_flags
56 
57   INTEGER, INTENT(IN   ) :: bzone, gn
58   INTEGER, INTENT(INOUT) :: error
59 
60 ! local variables
61 
62   INTEGER :: xs_bc, xe_bc, ys_bc, ye_bc, bzone_min
63   INTEGER :: nprocx, nprocy
64 
65   CALL wrf_debug( 100 , ' checking boundary conditions for grid ' )
66 
67   error = 0
68   xs_bc = 0
69   xe_bc = 0
70   ys_bc = 0
71   ye_bc = 0
72 
73 !  sum the number of conditions specified for each lateral boundary.
74 !  obviously, this number should be 1
75 
76   IF( config_flags%periodic_x ) THEN
77     xs_bc = xs_bc+1
78     xe_bc = xe_bc+1
79   ENDIF
80 
81   IF( config_flags%periodic_y ) THEN
82     ys_bc = ys_bc+1
83     ye_bc = ye_bc+1
84   ENDIF
85 
86   IF( config_flags%symmetric_xs ) xs_bc = xs_bc + 1
87   IF( config_flags%symmetric_xe ) xe_bc = xe_bc + 1
88   IF( config_flags%open_xs )      xs_bc = xs_bc + 1
89   IF( config_flags%open_xe )      xe_bc = xe_bc + 1
90 
91 
92   IF( config_flags%symmetric_ys ) ys_bc = ys_bc + 1
93   IF( config_flags%symmetric_ye ) ye_bc = ye_bc + 1
94   IF( config_flags%open_ys )      ys_bc = ys_bc + 1
95   IF( config_flags%open_ye )      ye_bc = ye_bc + 1
96 
97   IF( config_flags%nested ) THEN
98      xs_bc = xs_bc + 1
99      xe_bc = xe_bc + 1
100      ys_bc = ys_bc + 1
101      ye_bc = ye_bc + 1
102    ENDIF
103 
104   IF( config_flags%specified ) THEN
105      IF( .NOT. config_flags%periodic_x)xs_bc = xs_bc + 1
106      IF( .NOT. config_flags%periodic_x)xe_bc = xe_bc + 1
107      ys_bc = ys_bc + 1
108      ye_bc = ye_bc + 1
109    ENDIF
110 
111 !  check the number of conditions for each boundary
112 
113    IF( (xs_bc /= 1) .or. &
114        (xe_bc /= 1) .or. &
115        (ys_bc /= 1) .or. &
116        (ye_bc /= 1)         ) THEN
117 
118      error = 1
119 
120      write( wrf_err_message ,*) ' *** Error in boundary condition specification '
121      CALL wrf_message ( wrf_err_message )
122      write( wrf_err_message ,*) ' boundary conditions at xs ', xs_bc
123      CALL wrf_message ( wrf_err_message )
124      write( wrf_err_message ,*) ' boundary conditions at xe ', xe_bc
125      CALL wrf_message ( wrf_err_message )
126      write( wrf_err_message ,*) ' boundary conditions at ys ', ys_bc
127      CALL wrf_message ( wrf_err_message )
128      write( wrf_err_message ,*) ' boundary conditions at ye ', ye_bc
129      CALL wrf_message ( wrf_err_message )
130      write( wrf_err_message ,*) ' boundary conditions logicals are '
131      CALL wrf_message ( wrf_err_message )
132      write( wrf_err_message ,*) ' periodic_x   ',config_flags%periodic_x
133      CALL wrf_message ( wrf_err_message )
134      write( wrf_err_message ,*) ' periodic_y   ',config_flags%periodic_y
135      CALL wrf_message ( wrf_err_message )
136      write( wrf_err_message ,*) ' symmetric_xs ',config_flags%symmetric_xs
137      CALL wrf_message ( wrf_err_message )
138      write( wrf_err_message ,*) ' symmetric_xe ',config_flags%symmetric_xe
139      CALL wrf_message ( wrf_err_message )
140      write( wrf_err_message ,*) ' symmetric_ys ',config_flags%symmetric_ys
141      CALL wrf_message ( wrf_err_message )
142      write( wrf_err_message ,*) ' symmetric_ye ',config_flags%symmetric_ye
143      CALL wrf_message ( wrf_err_message )
144      write( wrf_err_message ,*) ' open_xs      ',config_flags%open_xs
145      CALL wrf_message ( wrf_err_message )
146      write( wrf_err_message ,*) ' open_xe      ',config_flags%open_xe
147      CALL wrf_message ( wrf_err_message )
148      write( wrf_err_message ,*) ' open_ys      ',config_flags%open_ys
149      CALL wrf_message ( wrf_err_message )
150      write( wrf_err_message ,*) ' open_ye      ',config_flags%open_ye
151      CALL wrf_message ( wrf_err_message )
152      write( wrf_err_message ,*) ' nested       ',config_flags%nested
153      CALL wrf_message ( wrf_err_message )
154      write( wrf_err_message ,*) ' specified    ',config_flags%specified
155      CALL wrf_error_fatal( ' *** Error in boundary condition specification ' )
156 
157    ENDIF
158 
159 !  now check to see if boundary zone size is sufficient.
160 !  we could have the necessary boundary zone size be returned
161 !  to the calling routine.
162 
163    IF( config_flags%periodic_x   .or. &
164        config_flags%periodic_y   .or. &
165        config_flags%symmetric_xs .or. &
166        config_flags%symmetric_xe .or. &
167        config_flags%symmetric_ys .or. &
168        config_flags%symmetric_ye        )  THEN
169 
170        bzone_min = MAX( 1,                                  &
171                         (config_flags%h_mom_adv_order+1)/2, &
172                         (config_flags%h_sca_adv_order+1)/2 )
173 
174        IF( bzone < bzone_min) THEN  
175 
176          error = 2
177          WRITE ( wrf_err_message , * ) ' boundary zone not large enough '
178          CALL wrf_message ( wrf_err_message )
179          WRITE ( wrf_err_message , * ) ' boundary zone specified      ',bzone
180          CALL wrf_message ( wrf_err_message )
181          WRITE ( wrf_err_message , * ) ' minimum boundary zone needed ',bzone_min
182          CALL wrf_error_fatal ( wrf_err_message )
183 
184        ENDIF
185    ENDIF
186 
187    CALL wrf_debug ( 100 , ' boundary conditions OK for grid ' )
188 
189    END subroutine boundary_condition_check
190 
191 !--------------------------------------------------------------------------
192    SUBROUTINE set_physical_bc2d( dat, variable_in,  &
193                                  config_flags,           & 
194                                  ids,ide, jds,jde,   & ! domain dims
195                                  ims,ime, jms,jme,   & ! memory dims
196                                  ips,ipe, jps,jpe,   & ! patch  dims
197                                  its,ite, jts,jte   )      
198 
199 !  This subroutine sets the data in the boundary region, by direct
200 !  assignment if possible, for periodic and symmetric (wall)
201 !  boundary conditions.  Currently, we are only doing 1 variable
202 !  at a time - lots of overhead, so maybe this routine can be easily 
203 !  inlined later or we could pass multiple variables -
204 !  would probably want a largestep and smallstep version.
205 
206 !  15 Jan 99, Dave
207 !  Modified the incoming its,ite,jts,jte to truly be the tile size.
208 !  This required modifying the loop limits when the "istag" or "jstag"
209 !  is used, as this is only required at the end of the domain.
210 
211       IMPLICIT NONE
212 
213       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde
214       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme
215       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe
216       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte
217       CHARACTER,    INTENT(IN   )    :: variable_in
218 
219       CHARACTER                      :: variable
220 
221       REAL,  DIMENSION( ims:ime , jms:jme ) :: dat
222       TYPE( grid_config_rec_type ) config_flags
223 
224       INTEGER  :: i, j, istag, jstag, itime
225 
226       LOGICAL  :: debug, open_bc_copy
227 
228 !------------
229 
230       debug = .false.
231 
232       open_bc_copy = .false.
233 
234       variable = variable_in
235       IF ( variable_in .ge. 'A' .and. variable_in .le. 'Z' ) THEN
236         variable = CHAR( ICHAR(variable_in) - ICHAR('A') + ICHAR('a') )
237       ENDIF
238       IF ((variable == 'u') .or. (variable == 'v') .or.  &
239           (variable == 'w') .or. (variable == 't') .or.  &
240           (variable == 'x') .or. (variable == 'y') .or.  &
241           (variable == 'r') .or. (variable == 'p') ) open_bc_copy = .true.
242 
243 !  begin, first set a staggering variable
244 
245       istag = -1
246       jstag = -1
247 
248       IF ((variable == 'u') .or. (variable == 'x')) istag = 0
249       IF ((variable == 'v') .or. (variable == 'y')) jstag = 0
250 
251       if(debug) then
252         write(6,*) ' in bc2d, var is ',variable, istag, jstag
253         write(6,*) ' b.cs are ',  &
254       config_flags%periodic_x,  &
255       config_flags%periodic_y
256       end if
257       
258 
259 
260 !  periodic conditions.
261 !  note, patch must cover full range in periodic dir, or else
262 !  its intra-patch communication that is handled elsewheres.
263 !  symmetry conditions can always be handled here, because no
264 !  outside patch communication is needed
265 
266       periodicity_x:  IF( ( config_flags%periodic_x ) ) THEN 
267         IF ( ( ids == ips ) .and.  ( ide == ipe ) ) THEN  ! test if east and west both on-processor 
268           IF ( its == ids ) THEN
269 
270             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
271             DO i = 0,-(bdyzone-1),-1
272               dat(ids+i-1,j) = dat(ide+i-1,j)
273             ENDDO
274             ENDDO
275 
276           ENDIF
277 
278           IF ( ite == ide ) THEN
279 
280             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
281 !!          DO i = 1 , bdyzone
282             DO i = -istag , bdyzone
283               dat(ide+i+istag,j) = dat(ids+i+istag,j)
284             ENDDO
285             ENDDO
286 
287           ENDIF
288         ENDIF
289 
290       ELSE 
291 
292         symmetry_xs: IF( ( config_flags%symmetric_xs ) .and.  &
293                          ( its == ids )                  )  THEN
294 
295           IF ( (variable /= 'u') .and. (variable /= 'x') ) THEN
296 
297             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
298             DO i = 1, bdyzone
299               dat(ids-i,j) = dat(ids+i-1,j) !  here, dat(0) = dat(1), etc
300             ENDDO                             !  symmetry about dat(0.5) (u=0 pt)
301             ENDDO
302 
303           ELSE
304 
305             IF( variable == 'u' ) THEN
306 
307               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
308               DO i = 0, bdyzone-1
309                 dat(ids-i,j) = - dat(ids+i,j) ! here, u(0) = - u(2), etc
310               ENDDO                             !  normal b.c symmetry at u(1)
311               ENDDO
312 
313             ELSE
314 
315               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
316               DO i = 0, bdyzone-1
317                 dat(ids-i,j) =   dat(ids+i,j) ! here, phi(0) = phi(2), etc
318               ENDDO                             !  normal b.c symmetry at phi(1)
319               ENDDO
320 
321             END IF
322 
323           ENDIF
324 
325         ENDIF symmetry_xs
326 
327 
328 !  now the symmetry boundary at xe
329 
330         symmetry_xe: IF( ( config_flags%symmetric_xe ) .and.  &
331                          ( ite == ide )                  )  THEN
332 
333           IF ( (variable /= 'u') .and. (variable /= 'x') ) THEN
334 
335             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
336             DO i = 1, bdyzone
337               dat(ide+i-1,j) = dat(ide-i,j)  !  sym. about dat(ide-0.5)
338             ENDDO
339             ENDDO
340 
341           ELSE
342 
343             IF (variable == 'u' ) THEN
344 
345               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
346               DO i = 0, bdyzone-1
347                 dat(ide+i,j) = - dat(ide-i,j)  ! u(ide+1) = - u(ide-1), etc.
348               ENDDO
349               ENDDO
350 
351 
352             ELSE
353 
354               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
355               DO i = 0, bdyzone-1
356                 dat(ide+i,j) = dat(ide-i,j)  !  phi(ide+1) = phi(ide-1), etc.
357               ENDDO
358               ENDDO
359 
360             END IF
361 
362           END IF 
363 
364         END IF symmetry_xe
365 
366 
367 !  set open b.c in X copy into boundary zone here.  WCS, 19 March 2000
368 
369         open_xs: IF( ( config_flags%open_xs   .or. &
370                        config_flags%specified .or. &
371                        config_flags%nested            ) .and.  &
372                          ( its == ids ) .and. open_bc_copy  )  THEN
373 
374             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
375               dat(ids-1,j) = dat(ids,j) !  here, dat(0) = dat(1)
376               dat(ids-2,j) = dat(ids,j)
377               dat(ids-3,j) = dat(ids,j)
378             ENDDO
379 
380         ENDIF open_xs
381 
382 
383 !  now the open boundary copy at xe
384 
385         open_xe: IF( ( config_flags%open_xe   .or. &
386                        config_flags%specified .or. &
387                        config_flags%nested            ) .and.  &
388                           ( ite == ide ) .and. open_bc_copy  )  THEN
389 
390           IF ( variable /= 'u' .and. variable /= 'x') THEN
391 
392             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
393               dat(ide  ,j) = dat(ide-1,j) 
394               dat(ide+1,j) = dat(ide-1,j) 
395               dat(ide+2,j) = dat(ide-1,j) 
396             ENDDO
397 
398           ELSE
399 
400             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
401               dat(ide+1,j) = dat(ide,j)
402               dat(ide+2,j) = dat(ide,j)
403               dat(ide+3,j) = dat(ide,j)
404             ENDDO
405 
406           END IF 
407 
408         END IF open_xe
409 
410 !  end open b.c in X copy into boundary zone addition.  WCS, 19 March 2000
411 
412       END IF periodicity_x
413 
414 !  same procedure in y
415 
416       periodicity_y:  IF( ( config_flags%periodic_y ) ) THEN
417         IF ( ( jds == jps ) .and. ( jde == jpe ) )  THEN    ! test of both north and south on processor
418 
419           IF( jts == jds ) then
420 
421             DO j = 0, -(bdyzone-1), -1
422             DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
423               dat(i,jds+j-1) = dat(i,jde+j-1)
424             ENDDO
425             ENDDO
426 
427           END IF
428 
429           IF( jte == jde ) then
430 
431             DO j = -jstag, bdyzone
432             DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
433               dat(i,jde+j+jstag) = dat(i,jds+j+jstag)
434             ENDDO
435             ENDDO
436 
437           END IF
438 
439         END IF
440 
441       ELSE
442 
443         symmetry_ys: IF( ( config_flags%symmetric_ys ) .and.  &
444                          ( jts == jds)                   )  THEN
445 
446           IF ( (variable /= 'v') .and. (variable /= 'y') ) THEN
447 
448             DO j = 1, bdyzone
449             DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
450               dat(i,jds-j) = dat(i,jds+j-1) 
451             ENDDO
452             ENDDO
453 
454           ELSE
455 
456             IF (variable == 'v') THEN
457 
458               DO j = 1, bdyzone
459               DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
460                 dat(i,jds-j) = - dat(i,jds+j) 
461               ENDDO              
462               ENDDO
463 
464             ELSE
465 
466               DO j = 1, bdyzone
467               DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
468                 dat(i,jds-j) = dat(i,jds+j) 
469               ENDDO              
470               ENDDO
471 
472             END IF
473 
474           ENDIF
475 
476         ENDIF symmetry_ys
477 
478 !  now the symmetry boundary at ye
479 
480         symmetry_ye: IF( ( config_flags%symmetric_ye ) .and.  &
481                          ( jte == jde )                  )  THEN
482 
483           IF ( (variable /= 'v') .and. (variable /= 'y') ) THEN
484 
485             DO j = 1, bdyzone
486             DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
487               dat(i,jde+j-1) = dat(i,jde-j) 
488             ENDDO                               
489             ENDDO
490 
491           ELSE
492 
493             IF (variable == 'v' ) THEN
494 
495               DO j = 1, bdyzone
496               DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
497                 dat(i,jde+j) = - dat(i,jde-j)    ! bugfix: changed jds on rhs to jde , JM 20020410
498               ENDDO                               
499               ENDDO
500 
501             ELSE
502 
503               DO j = 1, bdyzone
504               DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
505                 dat(i,jde+j) = dat(i,jde-j)
506               ENDDO                               
507               ENDDO
508 
509             END IF
510 
511           ENDIF
512 
513         END IF symmetry_ye
514 
515 !  set open b.c in Y copy into boundary zone here.  WCS, 19 March 2000
516 
517         open_ys: IF( ( config_flags%open_ys   .or. &
518                        config_flags%specified .or. &
519                        config_flags%nested            ) .and.  &
520                          ( jts == jds) .and. open_bc_copy )  THEN
521 
522             DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
523               dat(i,jds-1) = dat(i,jds) 
524               dat(i,jds-2) = dat(i,jds) 
525               dat(i,jds-3) = dat(i,jds) 
526             ENDDO
527 
528         ENDIF open_ys
529 
530 !  now the open boundary copy at ye
531 
532         open_ye: IF( ( config_flags%open_ye   .or. &
533                        config_flags%specified .or. &
534                        config_flags%nested            ) .and.  &
535                          ( jte == jde ) .and. open_bc_copy )  THEN
536 
537           IF  (variable /= 'v' .and. variable /= 'y' ) THEN
538 
539             DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
540               dat(i,jde  ) = dat(i,jde-1) 
541               dat(i,jde+1) = dat(i,jde-1) 
542               dat(i,jde+2) = dat(i,jde-1) 
543             ENDDO                               
544 
545           ELSE
546 
547             DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
548               dat(i,jde+1) = dat(i,jde) 
549               dat(i,jde+2) = dat(i,jde) 
550               dat(i,jde+3) = dat(i,jde) 
551             ENDDO                               
552 
553           ENDIF
554 
555         END IF open_ye
556       
557 !  end open b.c in Y copy into boundary zone addition.  WCS, 19 March 2000
558 
559       END IF periodicity_y
560 
561 !  fix corners for doubly periodic domains
562 
563       IF ( config_flags%periodic_x .and. config_flags%periodic_y &
564            .and. (ids == ips) .and. (ide == ipe)                 &
565            .and. (jds == jps) .and. (jde == jpe)                   ) THEN
566 
567          IF ( (its == ids) .and. (jts == jds) ) THEN  ! lower left corner fill
568            DO j = 0, -(bdyzone-1), -1
569            DO i = 0, -(bdyzone-1), -1
570              dat(ids+i-1,jds+j-1) = dat(ide+i-1,jde+j-1)
571            ENDDO
572            ENDDO
573          END IF
574 
575          IF ( (ite == ide) .and. (jts == jds) ) THEN  ! lower right corner fill
576            DO j = 0, -(bdyzone-1), -1
577            DO i = 1, bdyzone
578              dat(ide+i+istag,jds+j-1) = dat(ids+i+istag,jde+j-1)
579            ENDDO
580            ENDDO
581          END IF
582 
583          IF ( (ite == ide) .and. (jte == jde) ) THEN  ! upper right corner fill
584            DO j = 1, bdyzone
585            DO i = 1, bdyzone
586              dat(ide+i+istag,jde+j+jstag) = dat(ids+i+istag,jds+j+jstag)
587            ENDDO
588            ENDDO
589          END IF
590 
591          IF ( (its == ids) .and. (jte == jde) ) THEN  ! upper left corner fill
592            DO j = 1, bdyzone
593            DO i = 0, -(bdyzone-1), -1
594              dat(ids+i-1,jde+j+jstag) = dat(ide+i-1,jds+j+jstag)
595            ENDDO
596            ENDDO
597          END IF
598 
599        END IF
600 
601    END SUBROUTINE set_physical_bc2d
602 
603 !-----------------------------------
604 
605    SUBROUTINE set_physical_bc3d( dat, variable_in,        &
606                                config_flags,                   & 
607                                ids,ide, jds,jde, kds,kde,  & ! domain dims
608                                ims,ime, jms,jme, kms,kme,  & ! memory dims
609                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
610                                its,ite, jts,jte, kts,kte )
611 
612 !  This subroutine sets the data in the boundary region, by direct
613 !  assignment if possible, for periodic and symmetric (wall)
614 !  boundary conditions.  Currently, we are only doing 1 variable
615 !  at a time - lots of overhead, so maybe this routine can be easily 
616 !  inlined later or we could pass multiple variables -
617 !  would probably want a largestep and smallstep version.
618 
619 !  15 Jan 99, Dave
620 !  Modified the incoming its,ite,jts,jte to truly be the tile size.
621 !  This required modifying the loop limits when the "istag" or "jstag"
622 !  is used, as this is only required at the end of the domain.
623 
624       IMPLICIT NONE
625 
626       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
627       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
628       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
629       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
630       CHARACTER,    INTENT(IN   )    :: variable_in
631 
632       CHARACTER                      :: variable
633 
634       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ) :: dat
635       TYPE( grid_config_rec_type ) config_flags
636 
637       INTEGER  :: i, j, k, istag, jstag, itime, k_end
638 
639       LOGICAL  :: debug, open_bc_copy
640 
641 !------------
642 
643       debug = .false.
644 
645       open_bc_copy = .false.
646 
647       variable = variable_in
648       IF ( variable_in .ge. 'A' .and. variable_in .le. 'Z' ) THEN
649         variable = CHAR( ICHAR(variable_in) - ICHAR('A') + ICHAR('a') )
650       ENDIF
651 
652       IF ((variable == 'u') .or. (variable == 'v') .or.     &
653           (variable == 'w') .or. (variable == 't') .or.     &
654           (variable == 'd') .or. (variable == 'e') .or. &
655           (variable == 'x') .or. (variable == 'y') .or. &
656           (variable == 'f') .or. (variable == 'r') .or. &
657           (variable == 'p')                        ) open_bc_copy = .true.
658 
659 !  begin, first set a staggering variable
660 
661       istag = -1
662       jstag = -1
663       k_end = max(1,min(kde-1,kte))
664 
665 
666       IF ((variable == 'u') .or. (variable == 'x')) istag = 0
667       IF ((variable == 'v') .or. (variable == 'y')) jstag = 0
668       IF ((variable == 'd') .or. (variable == 'xy')) then
669          istag = 0
670          jstag = 0
671       ENDIF
672       IF ((variable == 'e') ) then
673          istag = 0
674          k_end = min(kde,kte)
675       ENDIF
676 
677       IF ((variable == 'f') ) then
678          jstag = 0
679          k_end = min(kde,kte)
680       ENDIF
681 
682       IF ( variable == 'w')  k_end = min(kde,kte)
683 
684 !      k_end = kte
685 
686       if(debug) then
687         write(6,*) ' in bc, var is ',variable, istag, jstag, kte, k_end
688         write(6,*) ' b.cs are ',  &
689       config_flags%periodic_x,  &
690       config_flags%periodic_y
691       end if
692       
693 
694 
695 !  periodic conditions.
696 !  note, patch must cover full range in periodic dir, or else
697 !  its intra-patch communication that is handled elsewheres.
698 !  symmetry conditions can always be handled here, because no
699 !  outside patch communication is needed
700 
701       periodicity_x:  IF( ( config_flags%periodic_x ) ) THEN
702 
703         IF ( ( ids == ips ) .and. ( ide == ipe ) ) THEN  ! test if both east and west on-processor
704           IF ( its == ids ) THEN
705 
706             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
707             DO k = kts, k_end
708             DO i = 0,-(bdyzone-1),-1
709               dat(ids+i-1,k,j) = dat(ide+i-1,k,j)
710             ENDDO
711             ENDDO
712             ENDDO
713 
714           ENDIF
715 
716 
717           IF ( ite == ide ) THEN
718 
719             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
720             DO k = kts, k_end
721             DO i = -istag , bdyzone
722               dat(ide+i+istag,k,j) = dat(ids+i+istag,k,j)
723             ENDDO
724             ENDDO
725             ENDDO
726 
727           ENDIF
728 
729         ENDIF
730 
731       ELSE 
732 
733         symmetry_xs: IF( ( config_flags%symmetric_xs ) .and.  &
734                          ( its == ids )                  )  THEN
735 
736           IF ( (variable /= 'u') .and. (variable /= 'x') ) THEN
737 
738             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
739             DO k = kts, k_end
740             DO i = 1, bdyzone
741               dat(ids-i,k,j) = dat(ids+i-1,k,j) !  here, dat(0) = dat(1), etc
742             ENDDO                                 !  symmetry about dat(0.5) (u = 0 pt)
743             ENDDO
744             ENDDO
745 
746           ELSE
747 
748             IF ( variable == 'u' ) THEN
749 
750               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
751               DO k = kts, k_end
752               DO i = 1, bdyzone
753                 dat(ids-i,k,j) = - dat(ids+i,k,j) ! here, u(0) = - u(2), etc
754               ENDDO                                 !  normal b.c symmetry at u(1)
755               ENDDO
756               ENDDO
757 
758             ELSE
759 
760               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
761               DO k = kts, k_end
762               DO i = 1, bdyzone
763                 dat(ids-i,k,j) = dat(ids+i,k,j) ! here, phi(0) = phi(2), etc
764               ENDDO                               !  normal b.c symmetry at phi(1)
765               ENDDO
766               ENDDO
767 
768             END IF
769 
770           ENDIF
771 
772         ENDIF symmetry_xs
773 
774 
775 !  now the symmetry boundary at xe
776 
777         symmetry_xe: IF( ( config_flags%symmetric_xe ) .and.  &
778                          ( ite == ide )                  )  THEN
779 
780           IF ( (variable /= 'u') .and. (variable /= 'x') ) THEN
781 
782             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
783             DO k = kts, k_end
784             DO i = 1, bdyzone
785               dat(ide+i-1,k,j) = dat(ide-i,k,j)  !  sym. about dat(ide-0.5)
786             ENDDO
787             ENDDO
788             ENDDO
789 
790           ELSE
791 
792             IF (variable == 'u') THEN
793 
794               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
795               DO k = kts, k_end
796               DO i = 1, bdyzone
797                 dat(ide+i,k,j) = - dat(ide-i,k,j)  ! u(ide+1) = - u(ide-1), etc.
798               ENDDO
799               ENDDO
800               ENDDO
801 
802             ELSE
803 
804               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
805               DO k = kts, k_end
806               DO i = 1, bdyzone
807                 dat(ide+i,k,j) = dat(ide-i,k,j)  ! phi(ide+1) = - phi(ide-1), etc.
808               ENDDO
809               ENDDO
810               ENDDO
811 
812              END IF
813 
814           END IF 
815 
816         END IF symmetry_xe
817 
818 !  set open b.c in X copy into boundary zone here.  WCS, 19 March 2000
819 
820         open_xs: IF( ( config_flags%open_xs   .or. &
821                        config_flags%specified .or. &
822                        config_flags%nested            ) .and.  &
823                          ( its == ids ) .and. open_bc_copy  )  THEN
824 
825             DO j = jts-bdyzone, MIN(jte,jde+jstag)+bdyzone
826             DO k = kts, k_end
827               dat(ids-1,k,j) = dat(ids,k,j) !  here, dat(0) = dat(1), etc
828               dat(ids-2,k,j) = dat(ids,k,j)
829               dat(ids-3,k,j) = dat(ids,k,j)
830             ENDDO
831             ENDDO
832 
833         ENDIF open_xs
834 
835 
836 !  now the open_xe boundary copy 
837 
838         open_xe: IF( ( config_flags%open_xe   .or. &
839                        config_flags%specified .or. &
840                        config_flags%nested            ) .and.  &
841                          ( ite == ide ) .and. open_bc_copy )  THEN
842 
843           IF (variable /= 'u' .and. variable /= 'x' ) THEN
844 
845             DO j = jts-bdyzone, MIN(jte,jde+jstag)+bdyzone
846             DO k = kts, k_end
847               dat(ide  ,k,j) = dat(ide-1,k,j)
848               dat(ide+1,k,j) = dat(ide-1,k,j)
849               dat(ide+2,k,j) = dat(ide-1,k,j)
850             ENDDO
851             ENDDO
852 
853           ELSE
854 
855 !!!!!!! I am not sure about this one!  JM 20020402
856             DO j = MAX(jds,jts-1)-bdyzone, MIN(jte+1,jde+jstag)+bdyzone
857             DO k = kts, k_end
858               dat(ide+1,k,j) = dat(ide,k,j)
859               dat(ide+2,k,j) = dat(ide,k,j)
860               dat(ide+3,k,j) = dat(ide,k,j)
861             ENDDO
862             ENDDO
863 
864           END IF 
865 
866         END IF open_xe
867 
868 !  end open b.c in X copy into boundary zone addition.  WCS, 19 March 2000
869 
870       END IF periodicity_x
871 
872 !  same procedure in y
873 
874       periodicity_y:  IF( ( config_flags%periodic_y ) ) THEN
875         IF ( ( jds == jps ) .and. ( jde == jpe ) )  THEN      ! test if both north and south on processor
876           IF( jts == jds ) then
877 
878             DO j = 0, -(bdyzone-1), -1
879             DO k = kts, k_end
880             DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
881               dat(i,k,jds+j-1) = dat(i,k,jde+j-1)
882             ENDDO
883             ENDDO
884             ENDDO
885 
886           END IF
887 
888           IF( jte == jde ) then
889 
890             DO j = -jstag, bdyzone
891             DO k = kts, k_end
892             DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
893               dat(i,k,jde+j+jstag) = dat(i,k,jds+j+jstag)
894             ENDDO
895             ENDDO
896             ENDDO
897 
898           END IF
899 
900         END IF
901 
902       ELSE
903 
904         symmetry_ys: IF( ( config_flags%symmetric_ys ) .and.  &
905                          ( jts == jds)                   )  THEN
906 
907           IF ( (variable /= 'v') .and. (variable /= 'y') ) THEN
908 
909             DO j = 1, bdyzone
910             DO k = kts, k_end
911             DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
912               dat(i,k,jds-j) = dat(i,k,jds+j-1) 
913             ENDDO                               
914             ENDDO
915             ENDDO
916 
917           ELSE
918 
919             IF (variable == 'v') THEN
920 
921               DO j = 1, bdyzone
922               DO k = kts, k_end
923               DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
924                 dat(i,k,jds-j) = - dat(i,k,jds+j) 
925               ENDDO              
926               ENDDO
927               ENDDO
928 
929             ELSE
930 
931               DO j = 1, bdyzone
932               DO k = kts, k_end
933               DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
934                 dat(i,k,jds-j) = dat(i,k,jds+j) 
935               ENDDO              
936               ENDDO
937               ENDDO
938 
939             END IF
940 
941           ENDIF
942 
943         ENDIF symmetry_ys
944 
945 !  now the symmetry boundary at ye
946 
947         symmetry_ye: IF( ( config_flags%symmetric_ye ) .and.  &
948                          ( jte == jde )                  )  THEN
949 
950           IF ( (variable /= 'v') .and. (variable /= 'y') ) THEN
951 
952             DO j = 1, bdyzone
953             DO k = kts, k_end
954             DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
955               dat(i,k,jde+j-1) = dat(i,k,jde-j) 
956             ENDDO                               
957             ENDDO
958             ENDDO
959 
960           ELSE
961 
962             IF ( variable == 'v' ) THEN
963 
964               DO j = 1, bdyzone
965               DO k = kts, k_end
966               DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
967                 dat(i,k,jde+j) = - dat(i,k,jde-j) 
968               ENDDO                               
969               ENDDO
970               ENDDO
971 
972             ELSE
973 
974               DO j = 1, bdyzone
975               DO k = kts, k_end
976               DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
977                 dat(i,k,jde+j) = dat(i,k,jde-j) 
978               ENDDO                               
979               ENDDO
980               ENDDO
981 
982             END IF
983 
984           ENDIF
985 
986         END IF symmetry_ye
987       
988 !  set open b.c in Y copy into boundary zone here.  WCS, 19 March 2000
989 
990         open_ys: IF( ( config_flags%open_ys   .or. &
991                        config_flags%specified .or. &
992                        config_flags%nested            ) .and.  &
993                          ( jts == jds) .and. open_bc_copy )  THEN
994 
995             DO k = kts, k_end
996             DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
997               dat(i,k,jds-1) = dat(i,k,jds) 
998               dat(i,k,jds-2) = dat(i,k,jds) 
999               dat(i,k,jds-3) = dat(i,k,jds) 
1000             ENDDO
1001             ENDDO
1002 
1003         ENDIF open_ys
1004 
1005 !  now the open boundary copy at ye
1006 
1007         open_ye: IF( ( config_flags%open_ye   .or. &
1008                        config_flags%specified .or. &
1009                        config_flags%nested            ) .and.  &
1010                          ( jte == jde ) .and. open_bc_copy )  THEN
1011 
1012           IF (variable /= 'v' .and. variable /= 'y' ) THEN
1013 
1014             DO k = kts, k_end
1015             DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
1016               dat(i,k,jde  ) = dat(i,k,jde-1) 
1017               dat(i,k,jde+1) = dat(i,k,jde-1) 
1018               dat(i,k,jde+2) = dat(i,k,jde-1) 
1019             ENDDO                               
1020             ENDDO
1021 
1022           ELSE
1023 
1024             DO k = kts, k_end
1025             DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
1026               dat(i,k,jde+1) = dat(i,k,jde) 
1027               dat(i,k,jde+2) = dat(i,k,jde) 
1028               dat(i,k,jde+3) = dat(i,k,jde) 
1029             ENDDO                               
1030             ENDDO
1031 
1032           ENDIF
1033 
1034       END IF open_ye
1035 
1036 !  end open b.c in Y copy into boundary zone addition.  WCS, 19 March 2000
1037 
1038       END IF periodicity_y
1039 
1040 !  fix corners for doubly periodic domains
1041 
1042       IF ( config_flags%periodic_x .and. config_flags%periodic_y &
1043            .and. (ids == ips) .and. (ide == ipe)                 &
1044            .and. (jds == jps) .and. (jde == jpe)                   ) THEN
1045 
1046          IF ( (its == ids) .and. (jts == jds) ) THEN  ! lower left corner fill
1047            DO j = 0, -(bdyzone-1), -1
1048            DO k = kts, k_end
1049            DO i = 0, -(bdyzone-1), -1
1050              dat(ids+i-1,k,jds+j-1) = dat(ide+i-1,k,jde+j-1)
1051            ENDDO
1052            ENDDO
1053            ENDDO
1054          END IF
1055 
1056          IF ( (ite == ide) .and. (jts == jds) ) THEN  ! lower right corner fill
1057            DO j = 0, -(bdyzone-1), -1
1058            DO k = kts, k_end
1059            DO i = 1, bdyzone
1060              dat(ide+i+istag,k,jds+j-1) = dat(ids+i+istag,k,jde+j-1)
1061            ENDDO
1062            ENDDO
1063            ENDDO
1064          END IF
1065 
1066          IF ( (ite == ide) .and. (jte == jde) ) THEN  ! upper right corner fill
1067            DO j = 1, bdyzone
1068            DO k = kts, k_end
1069            DO i = 1, bdyzone
1070              dat(ide+i+istag,k,jde+j+jstag) = dat(ids+i+istag,k,jds+j+jstag)
1071            ENDDO
1072            ENDDO
1073            ENDDO
1074          END IF
1075 
1076          IF ( (its == ids) .and. (jte == jde) ) THEN  ! upper left corner fill
1077            DO j = 1, bdyzone
1078            DO k = kts, k_end
1079            DO i = 0, -(bdyzone-1), -1
1080              dat(ids+i-1,k,jde+j+jstag) = dat(ide+i-1,k,jds+j+jstag)
1081            ENDDO
1082            ENDDO
1083            ENDDO
1084          END IF
1085 
1086        END IF
1087 
1088    END SUBROUTINE set_physical_bc3d
1089 
1090    SUBROUTINE init_module_bc
1091    END SUBROUTINE init_module_bc
1092 
1093 !------------------------------------------------------------------------
1094    SUBROUTINE relax_bdytend   ( field, field_tend,                     &
1095                                 field_bdy_xs, field_bdy_xe,            &
1096                                 field_bdy_ys, field_bdy_ye,            &
1097                                 field_bdy_tend_xs, field_bdy_tend_xe,  &
1098                                 field_bdy_tend_ys, field_bdy_tend_ye,  &
1099                                 variable_in, config_flags,             &
1100                                 spec_bdy_width, spec_zone, relax_zone, &
1101                                 dtbc, fcx, gcx,             &
1102                                 ids,ide, jds,jde, kds,kde,  & ! domain dims
1103                                 ims,ime, jms,jme, kms,kme,  & ! memory dims
1104                                 ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1105                                 its,ite, jts,jte, kts,kte )
1106 
1107 !  This subroutine adds the tendencies in the boundary relaxation region, for specified
1108 !  boundary conditions.  
1109 !  spec_bdy_width is only used to dimension the boundary arrays.
1110 !  relax_zone is the inner edge of the boundary relaxation zone treated here.
1111 !  spec_zone is the width of the outer specified b.c.s that are not changed here.
1112 !  (JD July 2000)
1113 
1114       IMPLICIT NONE
1115 
1116       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1117       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1118       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
1119       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1120       INTEGER,      INTENT(IN   )    :: spec_bdy_width, spec_zone, relax_zone
1121       REAL,         INTENT(IN   )    :: dtbc
1122       CHARACTER,    INTENT(IN   )    :: variable_in
1123 
1124 
1125       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: field
1126       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field_tend
1127       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_xs, field_bdy_xe
1128       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_ys, field_bdy_ye
1129       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_tend_xs, field_bdy_tend_xe
1130       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_tend_ys, field_bdy_tend_ye
1131       REAL,  DIMENSION( spec_bdy_width ), INTENT(IN   ) :: fcx, gcx
1132       TYPE( grid_config_rec_type ) config_flags
1133 
1134       CHARACTER  :: variable
1135       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf, im1, ip1
1136       INTEGER    :: b_dist, b_limit
1137       REAL       :: fls0, fls1, fls2, fls3, fls4
1138       LOGICAL    :: periodic_x
1139 
1140       periodic_x = config_flags%periodic_x
1141       variable = variable_in
1142 
1143       IF (variable == 'U') variable = 'u'
1144       IF (variable == 'V') variable = 'v'
1145       IF (variable == 'M') variable = 'm'
1146       IF (variable == 'H') variable = 'h'
1147 
1148       ibs = ids
1149       ibe = ide-1
1150       itf = min(ite,ide-1)
1151       jbs = jds
1152       jbe = jde-1
1153       jtf = min(jte,jde-1)
1154       ktf = kde-1
1155       IF (variable == 'u') ibe = ide
1156       IF (variable == 'u') itf = min(ite,ide)
1157       IF (variable == 'v') jbe = jde
1158       IF (variable == 'v') jtf = min(jte,jde)
1159       IF (variable == 'm') ktf = kte
1160       IF (variable == 'h') ktf = kte
1161 
1162       IF (jts - jbs .lt. relax_zone) THEN
1163 ! Y-start boundary
1164         DO j = max(jts,jbs+spec_zone), min(jtf,jbs+relax_zone-1)
1165           b_dist = j - jbs
1166           b_limit = b_dist
1167           IF(periodic_x)b_limit = 0
1168           DO k = kts, ktf
1169             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1170               im1 = max(i-1,ibs)
1171               ip1 = min(i+1,ibe)
1172               fls0 = field_bdy_ys(i, k, b_dist+1) &
1173                    + dtbc * field_bdy_tend_ys(i, k, b_dist+1) &
1174                    - field(i,k,j)
1175               fls1 = field_bdy_ys(im1, k, b_dist+1) &
1176                    + dtbc * field_bdy_tend_ys(im1, k, b_dist+1) &
1177                    - field(im1,k,j)
1178               fls2 = field_bdy_ys(ip1, k, b_dist+1) &
1179                    + dtbc * field_bdy_tend_ys(ip1, k, b_dist+1) &
1180                    - field(ip1,k,j)
1181               fls3 = field_bdy_ys(i, k, b_dist) &
1182                    + dtbc * field_bdy_tend_ys(i, k, b_dist) &
1183                    - field(i,k,j-1)
1184               fls4 = field_bdy_ys(i, k, b_dist+2) &
1185                    + dtbc * field_bdy_tend_ys(i, k, b_dist+2) &
1186                    - field(i,k,j+1)
1187               field_tend(i,k,j) = field_tend(i,k,j) &
1188                                 + fcx(b_dist+1)*fls0 &
1189                                 - gcx(b_dist+1)*(fls1+fls2+fls3+fls4-4.*fls0)
1190             ENDDO
1191           ENDDO
1192         ENDDO
1193       ENDIF
1194 
1195       IF (jbe - jtf .lt. relax_zone) THEN
1196 ! Y-end boundary
1197         DO j = max(jts,jbe-relax_zone+1), min(jtf,jbe-spec_zone)
1198           b_dist = jbe - j
1199           b_limit = b_dist
1200           IF(periodic_x)b_limit = 0
1201           DO k = kts, ktf
1202             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1203               im1 = max(i-1,ibs)
1204               ip1 = min(i+1,ibe)
1205               fls0 = field_bdy_ye(i, k, b_dist+1) &
1206                    + dtbc * field_bdy_tend_ye(i, k, b_dist+1) &
1207                    - field(i,k,j)
1208               fls1 = field_bdy_ye(im1, k, b_dist+1) &
1209                    + dtbc * field_bdy_tend_ye(im1, k, b_dist+1) &
1210                    - field(im1,k,j)
1211               fls2 = field_bdy_ye(ip1, k, b_dist+1) &
1212                    + dtbc * field_bdy_tend_ye(ip1, k, b_dist+1) &
1213                    - field(ip1,k,j)
1214               fls3 = field_bdy_ye(i, k, b_dist) &
1215                    + dtbc * field_bdy_tend_ye(i, k, b_dist) &
1216                    - field(i,k,j+1)
1217               fls4 = field_bdy_ye(i, k, b_dist+2) &
1218                    + dtbc * field_bdy_tend_ye(i, k, b_dist+2) &
1219                    - field(i,k,j-1)
1220               field_tend(i,k,j) = field_tend(i,k,j) &
1221                                 + fcx(b_dist+1)*fls0 &
1222                                 - gcx(b_dist+1)*(fls1+fls2+fls3+fls4-4.*fls0)
1223 
1224             ENDDO
1225           ENDDO
1226         ENDDO
1227       ENDIF
1228 
1229     IF(.NOT.periodic_x)THEN
1230       IF (its - ibs .lt. relax_zone) THEN
1231 ! X-start boundary
1232         DO i = max(its,ibs+spec_zone), min(itf,ibs+relax_zone-1)
1233           b_dist = i - ibs
1234           DO k = kts, ktf
1235             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1236               fls0 = field_bdy_xs(j, k, b_dist+1) &
1237                    + dtbc * field_bdy_tend_xs(j, k, b_dist+1) &
1238                    - field(i,k,j)
1239               fls1 = field_bdy_xs(j-1, k, b_dist+1) &
1240                    + dtbc * field_bdy_tend_xs(j-1, k, b_dist+1) &
1241                    - field(i,k,j-1)
1242               fls2 = field_bdy_xs(j+1, k, b_dist+1) &
1243                    + dtbc * field_bdy_tend_xs(j+1, k, b_dist+1) &
1244                    - field(i,k,j+1)
1245               fls3 = field_bdy_xs(j, k, b_dist) &
1246                    + dtbc * field_bdy_tend_xs(j, k, b_dist) &
1247                    - field(i-1,k,j)
1248               fls4 = field_bdy_xs(j, k, b_dist+2) &
1249                    + dtbc * field_bdy_tend_xs(j, k, b_dist+2) &
1250                    - field(i+1,k,j)
1251               field_tend(i,k,j) = field_tend(i,k,j) &
1252                                 + fcx(b_dist+1)*fls0 &
1253                                 - gcx(b_dist+1)*(fls1+fls2+fls3+fls4-4.*fls0)
1254 
1255             ENDDO
1256           ENDDO
1257         ENDDO
1258       ENDIF
1259 
1260       IF (ibe - itf .lt. relax_zone) THEN
1261 ! X-end boundary
1262         DO i = max(its,ibe-relax_zone+1), min(itf,ibe-spec_zone)
1263           b_dist = ibe - i
1264           DO k = kts, ktf
1265             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1266               fls0 = field_bdy_xe(j, k, b_dist+1) &
1267                    + dtbc * field_bdy_tend_xe(j, k, b_dist+1) &
1268                    - field(i,k,j)
1269               fls1 = field_bdy_xe(j-1, k, b_dist+1) &
1270                    + dtbc * field_bdy_tend_xe(j-1, k, b_dist+1) &
1271                    - field(i,k,j-1)
1272               fls2 = field_bdy_xe(j+1, k, b_dist+1) &
1273                    + dtbc * field_bdy_tend_xe(j+1, k, b_dist+1) &
1274                    - field(i,k,j+1)
1275               fls3 = field_bdy_xe(j, k, b_dist) &
1276                    + dtbc * field_bdy_tend_xe(j, k, b_dist) &
1277                    - field(i+1,k,j)
1278               fls4 = field_bdy_xe(j, k, b_dist+2) &
1279                    + dtbc * field_bdy_tend_xe(j, k, b_dist+2) &
1280                    - field(i-1,k,j)
1281               field_tend(i,k,j) = field_tend(i,k,j) &
1282                                 + fcx(b_dist+1)*fls0 &
1283                                 - gcx(b_dist+1)*(fls1+fls2+fls3+fls4-4.*fls0)
1284             ENDDO
1285           ENDDO
1286         ENDDO
1287       ENDIF
1288     ENDIF
1289 
1290 
1291    END SUBROUTINE relax_bdytend
1292 !------------------------------------------------------------------------
1293 
1294    SUBROUTINE spec_bdytend   ( field_tend,                           &
1295                                field_bdy_xs, field_bdy_xe,           &
1296                                field_bdy_ys, field_bdy_ye,           &
1297                                field_bdy_tend_xs, field_bdy_tend_xe, &
1298                                field_bdy_tend_ys, field_bdy_tend_ye, &
1299                                variable_in, config_flags, & 
1300                                spec_bdy_width, spec_zone, &
1301                                ids,ide, jds,jde, kds,kde,  & ! domain dims
1302                                ims,ime, jms,jme, kms,kme,  & ! memory dims
1303                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1304                                its,ite, jts,jte, kts,kte )
1305 
1306 !  This subroutine sets the tendencies in the boundary specified region.
1307 !  spec_bdy_width is only used to dimension the boundary arrays.
1308 !  spec_zone is the width of the outer specified b.c.s that are set here.
1309 !  (JD July 2000)
1310 
1311       IMPLICIT NONE
1312 
1313       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1314       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1315       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
1316       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1317       INTEGER,      INTENT(IN   )    :: spec_bdy_width, spec_zone
1318       CHARACTER,    INTENT(IN   )    :: variable_in
1319 
1320 
1321       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(OUT  ) :: field_tend
1322       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_xs, field_bdy_xe
1323       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_ys, field_bdy_ye 
1324       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_tend_xs, field_bdy_tend_xe
1325       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_tend_ys, field_bdy_tend_ye 
1326       TYPE( grid_config_rec_type ) config_flags
1327 
1328       CHARACTER  :: variable
1329       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf
1330       INTEGER    :: b_dist, b_limit
1331       LOGICAL    :: periodic_x
1332 
1333       periodic_x = config_flags%periodic_x
1334 
1335       variable = variable_in
1336 
1337       IF (variable == 'U') variable = 'u'
1338       IF (variable == 'V') variable = 'v'
1339       IF (variable == 'M') variable = 'm'
1340       IF (variable == 'H') variable = 'h'
1341 
1342       ibs = ids
1343       ibe = ide-1
1344       itf = min(ite,ide-1)
1345       jbs = jds
1346       jbe = jde-1
1347       jtf = min(jte,jde-1)
1348       ktf = kde-1
1349       IF (variable == 'u') ibe = ide
1350       IF (variable == 'u') itf = min(ite,ide)
1351       IF (variable == 'v') jbe = jde
1352       IF (variable == 'v') jtf = min(jte,jde)
1353       IF (variable == 'm') ktf = kte
1354       IF (variable == 'h') ktf = kte
1355 
1356       IF (jts - jbs .lt. spec_zone) THEN
1357 ! Y-start boundary
1358         DO j = jts, min(jtf,jbs+spec_zone-1)
1359           b_dist = j - jbs
1360           b_limit = b_dist
1361           IF(periodic_x)b_limit = 0
1362           DO k = kts, ktf
1363             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1364               field_tend(i,k,j) = field_bdy_tend_ys(i, k, b_dist+1)
1365             ENDDO
1366           ENDDO
1367         ENDDO
1368       ENDIF 
1369       IF (jbe - jtf .lt. spec_zone) THEN 
1370 ! Y-end boundary 
1371         DO j = max(jts,jbe-spec_zone+1), jtf 
1372           b_dist = jbe - j 
1373           b_limit = b_dist
1374           IF(periodic_x)b_limit = 0
1375           DO k = kts, ktf 
1376             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1377               field_tend(i,k,j) = field_bdy_tend_ye(i, k, b_dist+1)
1378             ENDDO
1379           ENDDO
1380         ENDDO
1381       ENDIF 
1382 
1383     IF(.NOT.periodic_x)THEN
1384       IF (its - ibs .lt. spec_zone) THEN
1385 ! X-start boundary
1386         DO i = its, min(itf,ibs+spec_zone-1)
1387           b_dist = i - ibs
1388           DO k = kts, ktf
1389             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1390               field_tend(i,k,j) = field_bdy_tend_xs(j, k, b_dist+1)
1391             ENDDO
1392           ENDDO
1393         ENDDO
1394       ENDIF 
1395 
1396       IF (ibe - itf .lt. spec_zone) THEN
1397 ! X-end boundary
1398         DO i = max(its,ibe-spec_zone+1), itf
1399           b_dist = ibe - i
1400           DO k = kts, ktf
1401             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1402               field_tend(i,k,j) = field_bdy_tend_xe(j, k, b_dist+1)
1403             ENDDO
1404           ENDDO
1405         ENDDO
1406       ENDIF 
1407     ENDIF
1408 
1409    END SUBROUTINE spec_bdytend
1410 !------------------------------------------------------------------------
1411 
1412    SUBROUTINE spec_bdyupdate(  field,      &
1413                                field_tend, dt,            &
1414                                variable_in, config_flags, & 
1415                                spec_zone,                  &
1416                                ids,ide, jds,jde, kds,kde,  & ! domain dims
1417                                ims,ime, jms,jme, kms,kme,  & ! memory dims
1418                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1419                                its,ite, jts,jte, kts,kte )
1420 
1421 !  This subroutine adds the tendencies in the boundary specified region.
1422 !  spec_zone is the width of the outer specified b.c.s that are set here.
1423 !  (JD August 2000)
1424 
1425       IMPLICIT NONE
1426 
1427       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1428       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1429       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
1430       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1431       INTEGER,      INTENT(IN   )    :: spec_zone
1432       CHARACTER,    INTENT(IN   )    :: variable_in
1433       REAL,         INTENT(IN   )    :: dt
1434 
1435 
1436       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field
1437       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: field_tend
1438       TYPE( grid_config_rec_type ) config_flags
1439 
1440       CHARACTER  :: variable
1441       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf
1442       INTEGER    :: b_dist, b_limit
1443       LOGICAL    :: periodic_x
1444 
1445       periodic_x = config_flags%periodic_x
1446 
1447       variable = variable_in
1448 
1449       IF (variable == 'U') variable = 'u'
1450       IF (variable == 'V') variable = 'v'
1451       IF (variable == 'M') variable = 'm'
1452       IF (variable == 'H') variable = 'h'
1453 
1454       ibs = ids
1455       ibe = ide-1
1456       itf = min(ite,ide-1)
1457       jbs = jds
1458       jbe = jde-1
1459       jtf = min(jte,jde-1)
1460       ktf = kde-1
1461       IF (variable == 'u') ibe = ide
1462       IF (variable == 'u') itf = min(ite,ide)
1463       IF (variable == 'v') jbe = jde
1464       IF (variable == 'v') jtf = min(jte,jde)
1465       IF (variable == 'm') ktf = kte
1466       IF (variable == 'h') ktf = kte
1467 
1468       IF (jts - jbs .lt. spec_zone) THEN
1469 ! Y-start boundary
1470         DO j = jts, min(jtf,jbs+spec_zone-1)
1471           b_dist = j - jbs
1472           b_limit = b_dist
1473           IF(periodic_x)b_limit = 0
1474           DO k = kts, ktf
1475             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1476               field(i,k,j) = field(i,k,j) + dt*field_tend(i,k,j) 
1477             ENDDO
1478           ENDDO
1479         ENDDO
1480       ENDIF 
1481       IF (jbe - jtf .lt. spec_zone) THEN 
1482 ! Y-end boundary 
1483         DO j = max(jts,jbe-spec_zone+1), jtf 
1484           b_dist = jbe - j 
1485           b_limit = b_dist
1486           IF(periodic_x)b_limit = 0
1487           DO k = kts, ktf 
1488             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1489               field(i,k,j) = field(i,k,j) + dt*field_tend(i,k,j) 
1490             ENDDO
1491           ENDDO
1492         ENDDO
1493       ENDIF 
1494 
1495     IF(.NOT.periodic_x)THEN
1496       IF (its - ibs .lt. spec_zone) THEN
1497 ! X-start boundary
1498         DO i = its, min(itf,ibs+spec_zone-1)
1499           b_dist = i - ibs
1500           DO k = kts, ktf
1501             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1502               field(i,k,j) = field(i,k,j) + dt*field_tend(i,k,j) 
1503             ENDDO
1504           ENDDO
1505         ENDDO
1506       ENDIF 
1507 
1508       IF (ibe - itf .lt. spec_zone) THEN
1509 ! X-end boundary
1510         DO i = max(its,ibe-spec_zone+1), itf
1511           b_dist = ibe - i
1512           DO k = kts, ktf
1513             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1514               field(i,k,j) = field(i,k,j) + dt*field_tend(i,k,j) 
1515             ENDDO
1516           ENDDO
1517         ENDDO
1518       ENDIF 
1519     ENDIF
1520 
1521    END SUBROUTINE spec_bdyupdate
1522 !------------------------------------------------------------------------
1523 
1524    SUBROUTINE zero_grad_bdy (  field,                     &
1525                                variable_in, config_flags, & 
1526                                spec_zone,                  &
1527                                ids,ide, jds,jde, kds,kde,  & ! domain dims
1528                                ims,ime, jms,jme, kms,kme,  & ! memory dims
1529                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1530                                its,ite, jts,jte, kts,kte )
1531 
1532 !  This subroutine sets zero gradient conditions in the boundary specified region.
1533 !  spec_zone is the width of the outer specified b.c.s that are set here.
1534 !  (JD August 2000)
1535 
1536       IMPLICIT NONE
1537 
1538       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1539       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1540       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
1541       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1542       INTEGER,      INTENT(IN   )    :: spec_zone
1543       CHARACTER,    INTENT(IN   )    :: variable_in
1544 
1545 
1546       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field
1547       TYPE( grid_config_rec_type ) config_flags
1548 
1549       CHARACTER  :: variable
1550       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf, i_inner, j_inner
1551       INTEGER    :: b_dist, b_limit
1552       LOGICAL    :: periodic_x
1553 
1554       periodic_x = config_flags%periodic_x
1555 
1556       variable = variable_in
1557 
1558       IF (variable == 'U') variable = 'u'
1559       IF (variable == 'V') variable = 'v'
1560 
1561       ibs = ids
1562       ibe = ide-1
1563       itf = min(ite,ide-1)
1564       jbs = jds
1565       jbe = jde-1
1566       jtf = min(jte,jde-1)
1567       ktf = kde-1
1568       IF (variable == 'u') ibe = ide
1569       IF (variable == 'u') itf = min(ite,ide)
1570       IF (variable == 'v') jbe = jde
1571       IF (variable == 'v') jtf = min(jte,jde)
1572       IF (variable == 'w') ktf = kde
1573 
1574       IF (jts - jbs .lt. spec_zone) THEN
1575 ! Y-start boundary
1576         DO j = jts, min(jtf,jbs+spec_zone-1)
1577           b_dist = j - jbs
1578           b_limit = b_dist
1579           IF(periodic_x)b_limit = 0
1580           DO k = kts, ktf
1581             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1582               i_inner = max(i,ibs+spec_zone)
1583               i_inner = min(i_inner,ibe-spec_zone)
1584               IF(periodic_x)i_inner = i
1585               field(i,k,j) = field(i_inner,k,jbs+spec_zone)
1586             ENDDO
1587           ENDDO
1588         ENDDO
1589       ENDIF 
1590       IF (jbe - jtf .lt. spec_zone) THEN 
1591 ! Y-end boundary 
1592         DO j = max(jts,jbe-spec_zone+1), jtf 
1593           b_dist = jbe - j 
1594           b_limit = b_dist
1595           IF(periodic_x)b_limit = 0
1596           DO k = kts, ktf 
1597             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1598               i_inner = max(i,ibs+spec_zone)
1599               i_inner = min(i_inner,ibe-spec_zone)
1600               IF(periodic_x)i_inner = i
1601               field(i,k,j) = field(i_inner,k,jbe-spec_zone)
1602             ENDDO
1603           ENDDO
1604         ENDDO
1605       ENDIF 
1606 
1607     IF(.NOT.periodic_x)THEN
1608       IF (its - ibs .lt. spec_zone) THEN
1609 ! X-start boundary
1610         DO i = its, min(itf,ibs+spec_zone-1)
1611           b_dist = i - ibs
1612           DO k = kts, ktf
1613             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1614               j_inner = max(j,jbs+spec_zone)
1615               j_inner = min(j_inner,jbe-spec_zone)
1616               field(i,k,j) = field(ibs+spec_zone,k,j_inner)
1617             ENDDO
1618           ENDDO
1619         ENDDO
1620       ENDIF 
1621 
1622       IF (ibe - itf .lt. spec_zone) THEN
1623 ! X-end boundary
1624         DO i = max(its,ibe-spec_zone+1), itf
1625           b_dist = ibe - i
1626           DO k = kts, ktf
1627             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1628               j_inner = max(j,jbs+spec_zone)
1629               j_inner = min(j_inner,jbe-spec_zone)
1630               field(i,k,j) = field(ibe-spec_zone,k,j_inner)
1631             ENDDO
1632           ENDDO
1633         ENDDO
1634       ENDIF 
1635     ENDIF
1636 
1637    END SUBROUTINE zero_grad_bdy
1638 !------------------------------------------------------------------------
1639 
1640    SUBROUTINE flow_dep_bdy  (  field,                     &
1641                                u, v, config_flags, & 
1642                                spec_zone,                  &
1643                                ids,ide, jds,jde, kds,kde,  & ! domain dims
1644                                ims,ime, jms,jme, kms,kme,  & ! memory dims
1645                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1646                                its,ite, jts,jte, kts,kte )
1647 
1648 !  This subroutine sets zero gradient conditions for outflow and zero value
1649 !  for inflow in the boundary specified region. Note that field must be unstaggered.
1650 !  The velocities, u and v, will only be used to check their sign (coupled vels OK)
1651 !  spec_zone is the width of the outer specified b.c.s that are set here.
1652 !  (JD August 2000)
1653 
1654       IMPLICIT NONE
1655 
1656       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1657       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1658       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
1659       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1660       INTEGER,      INTENT(IN   )    :: spec_zone
1661 
1662 
1663       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field
1664       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: u
1665       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: v
1666       TYPE( grid_config_rec_type ) config_flags
1667 
1668       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf, i_inner, j_inner
1669       INTEGER    :: b_dist, b_limit
1670       LOGICAL    :: periodic_x
1671 
1672       periodic_x = config_flags%periodic_x
1673 
1674       ibs = ids
1675       ibe = ide-1
1676       itf = min(ite,ide-1)
1677       jbs = jds
1678       jbe = jde-1
1679       jtf = min(jte,jde-1)
1680       ktf = kde-1
1681 
1682       IF (jts - jbs .lt. spec_zone) THEN
1683 ! Y-start boundary
1684         DO j = jts, min(jtf,jbs+spec_zone-1)
1685           b_dist = j - jbs
1686           b_limit = b_dist
1687           IF(periodic_x)b_limit = 0
1688           DO k = kts, ktf
1689             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1690               i_inner = max(i,ibs+spec_zone)
1691               i_inner = min(i_inner,ibe-spec_zone)
1692               IF(periodic_x)i_inner = i
1693               IF(v(i,k,j) .lt. 0.)THEN
1694                 field(i,k,j) = field(i_inner,k,jbs+spec_zone)
1695               ELSE
1696                 field(i,k,j) = 0.
1697               ENDIF
1698             ENDDO
1699           ENDDO
1700         ENDDO
1701       ENDIF 
1702       IF (jbe - jtf .lt. spec_zone) THEN 
1703 ! Y-end boundary 
1704         DO j = max(jts,jbe-spec_zone+1), jtf 
1705           b_dist = jbe - j 
1706           b_limit = b_dist
1707           IF(periodic_x)b_limit = 0
1708           DO k = kts, ktf 
1709             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1710               i_inner = max(i,ibs+spec_zone)
1711               i_inner = min(i_inner,ibe-spec_zone)
1712               IF(periodic_x)i_inner = i
1713               IF(v(i,k,j+1) .gt. 0.)THEN
1714                 field(i,k,j) = field(i_inner,k,jbe-spec_zone)
1715               ELSE
1716                 field(i,k,j) = 0.
1717               ENDIF
1718             ENDDO
1719           ENDDO
1720         ENDDO
1721       ENDIF 
1722 
1723     IF(.NOT.periodic_x)THEN
1724       IF (its - ibs .lt. spec_zone) THEN
1725 ! X-start boundary
1726         DO i = its, min(itf,ibs+spec_zone-1)
1727           b_dist = i - ibs
1728           DO k = kts, ktf
1729             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1730               j_inner = max(j,jbs+spec_zone)
1731               j_inner = min(j_inner,jbe-spec_zone)
1732               IF(u(i,k,j) .lt. 0.)THEN
1733                 field(i,k,j) = field(ibs+spec_zone,k,j_inner)
1734               ELSE
1735                 field(i,k,j) = 0.
1736               ENDIF
1737             ENDDO
1738           ENDDO
1739         ENDDO
1740       ENDIF 
1741 
1742       IF (ibe - itf .lt. spec_zone) THEN
1743 ! X-end boundary
1744         DO i = max(its,ibe-spec_zone+1), itf
1745           b_dist = ibe - i
1746           DO k = kts, ktf
1747             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1748               j_inner = max(j,jbs+spec_zone)
1749               j_inner = min(j_inner,jbe-spec_zone)
1750               IF(u(i+1,k,j) .gt. 0.)THEN
1751                 field(i,k,j) = field(ibe-spec_zone,k,j_inner)
1752               ELSE
1753                 field(i,k,j) = 0.
1754               ENDIF
1755             ENDDO
1756           ENDDO
1757         ENDDO
1758       ENDIF 
1759     ENDIF
1760 
1761    END SUBROUTINE flow_dep_bdy
1762 
1763 !------------------------------------------------------------------------------
1764 
1765  SUBROUTINE stuff_bdy_new ( data3d , space_bdy_xs, space_bdy_xe, space_bdy_ys, space_bdy_ye, &
1766                              char_stagger , &
1767                              spec_bdy_width , &
1768                              ids, ide, jds, jde, kds, kde , &
1769                              ims, ime, jms, jme, kms, kme , & 
1770                              its, ite, jts, jte, kts, kte )
1771  
1772  !  This routine puts the data in the 3d arrays into the proper locations
1773  !  for the lateral boundary arrays.
1774  
1775     USE module_state_description
1776     
1777     IMPLICIT NONE
1778  
1779     INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde
1780     INTEGER , INTENT(IN) :: ims, ime, jms, jme, kms, kme
1781     INTEGER , INTENT(IN) :: its, ite, jts, jte, kts, kte
1782     INTEGER , INTENT(IN) :: spec_bdy_width
1783     REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: data3d
1784     REAL , DIMENSION(jms:jme,kds:kde,spec_bdy_width) , INTENT(OUT) :: space_bdy_xs, space_bdy_xe
1785     REAL , DIMENSION(ims:ime,kds:kde,spec_bdy_width) , INTENT(OUT) :: space_bdy_ys, space_bdy_ye
1786     CHARACTER (LEN=1) , INTENT(IN) :: char_stagger
1787  
1788     INTEGER :: i , ii , j , jj , k
1789  
1790     !  There are four lateral boundary locations that are stored.
1791  
1792     !  X start boundary
1793  
1794     IF ( char_stagger .EQ. 'W' ) THEN
1795        DO j = MAX(jds,jts) , MIN(jde-1,jte)
1796        DO k = kds , kde
1797        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
1798           space_bdy_xs(j,k,i) = data3d(i,k,j)
1799        END DO
1800        END DO
1801        END DO
1802     ELSE IF ( char_stagger .EQ. 'M' ) THEN
1803        DO j = MAX(jds,jts) , MIN(jde-1,jte)
1804        DO k = kds , kde
1805        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
1806           space_bdy_xs(j,k,i) = data3d(i,k,j)
1807        END DO
1808        END DO
1809        END DO
1810     ELSE IF ( char_stagger .EQ. 'V' ) THEN
1811        DO j = MAX(jds,jts) , MIN(jde,jte)
1812        DO k = kds , kde - 1
1813        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
1814           space_bdy_xs(j,k,i) = data3d(i,k,j)
1815        END DO
1816        END DO
1817        END DO
1818     ELSE
1819        DO j = MAX(jds,jts) , MIN(jde-1,jte)
1820        DO k = kds , kde - 1
1821        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
1822           space_bdy_xs(j,k,i) = data3d(i,k,j)
1823        END DO
1824        END DO
1825        END DO
1826     END IF
1827  
1828     !  X end boundary
1829  
1830     IF      ( char_stagger .EQ. 'U' ) THEN
1831        DO j = MAX(jds,jts) , MIN(jde-1,jte)
1832        DO k = kds , kde - 1
1833        DO i = MIN(ide,ite) , MAX(ide - spec_bdy_width + 1,its) , -1
1834           ii = ide - i + 1
1835           space_bdy_xe(j,k,ii) = data3d(i,k,j)
1836        END DO
1837        END DO
1838        END DO
1839     ELSE IF ( char_stagger .EQ. 'V' ) THEN
1840        DO j = MAX(jds,jts) , MIN(jde,jte)
1841        DO k = kds , kde - 1
1842        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
1843           ii = ide - i
1844           space_bdy_xe(j,k,ii) = data3d(i,k,j)
1845        END DO
1846        END DO
1847        END DO
1848     ELSE IF ( char_stagger .EQ. 'W' ) THEN
1849        DO j = MAX(jds,jts) , MIN(jde-1,jte)
1850        DO k = kds , kde
1851        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
1852           ii = ide - i
1853           space_bdy_xe(j,k,ii) = data3d(i,k,j)
1854        END DO
1855        END DO
1856        END DO
1857     ELSE IF ( char_stagger .EQ. 'M' ) THEN
1858        DO j = MAX(jds,jts) , MIN(jde-1,jte)
1859        DO k = kds , kde
1860        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
1861           ii = ide - i
1862           space_bdy_xe(j,k,ii) = data3d(i,k,j)
1863        END DO
1864        END DO
1865        END DO
1866     ELSE
1867        DO j = MAX(jds,jts) , MIN(jde-1,jte)
1868        DO k = kds , kde - 1
1869        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
1870           ii = ide - i
1871           space_bdy_xe(j,k,ii) = data3d(i,k,j)
1872        END DO
1873        END DO
1874        END DO
1875     END IF
1876  
1877     !  Y start boundary
1878  
1879     IF ( char_stagger .EQ. 'W' ) THEN
1880        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
1881        DO k = kds , kde
1882        DO i = MAX(ids,its) , MIN(ide-1,ite)
1883           space_bdy_ys(i,k,j) = data3d(i,k,j)
1884        END DO
1885        END DO
1886        END DO
1887     ELSE IF ( char_stagger .EQ. 'M' ) THEN
1888        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
1889        DO k = kds , kde
1890        DO i = MAX(ids,its) , MIN(ide-1,ite)
1891           space_bdy_ys(i,k,j) = data3d(i,k,j)
1892        END DO
1893        END DO
1894        END DO
1895     ELSE IF ( char_stagger .EQ. 'U' ) THEN
1896        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
1897        DO k = kds , kde - 1
1898        DO i = MAX(ids,its) , MIN(ide,ite)
1899           space_bdy_ys(i,k,j) = data3d(i,k,j)
1900        END DO
1901        END DO
1902        END DO
1903     ELSE
1904        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
1905        DO k = kds , kde - 1
1906        DO i = MAX(ids,its) , MIN(ide-1,ite)
1907           space_bdy_ys(i,k,j) = data3d(i,k,j)
1908        END DO
1909        END DO
1910        END DO
1911     END IF
1912  
1913     !  Y end boundary
1914  
1915     IF      ( char_stagger .EQ. 'V' ) THEN
1916        DO j = MIN(jde,jte) , MAX(jde - spec_bdy_width + 1,jts) , -1
1917        DO k = kds , kde - 1
1918        DO i = MAX(ids,its) , MIN(ide-1,ite)
1919           jj = jde - j + 1
1920           space_bdy_ye(i,k,jj) = data3d(i,k,j)
1921        END DO
1922        END DO
1923        END DO
1924     ELSE IF ( char_stagger .EQ. 'U' ) THEN
1925        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
1926        DO k = kds , kde - 1
1927        DO i = MAX(ids,its) , MIN(ide,ite)
1928           jj = jde - j
1929           space_bdy_ye(i,k,jj) = data3d(i,k,j)
1930        END DO
1931        END DO
1932        END DO
1933     ELSE IF ( char_stagger .EQ. 'W' ) THEN
1934        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
1935        DO k = kds , kde
1936        DO i = MAX(ids,its) , MIN(ide-1,ite)
1937           jj = jde - j
1938           space_bdy_ye(i,k,jj) = data3d(i,k,j)
1939        END DO
1940        END DO
1941        END DO
1942     ELSE IF ( char_stagger .EQ. 'M' ) THEN
1943        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
1944        DO k = kds , kde
1945        DO i = MAX(ids,its) , MIN(ide-1,ite)
1946           jj = jde - j
1947           space_bdy_ye(i,k,jj) = data3d(i,k,j)
1948        END DO
1949        END DO
1950        END DO
1951     ELSE
1952        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
1953        DO k = kds , kde - 1
1954        DO i = MAX(ids,its) , MIN(ide-1,ite)
1955           jj = jde - j
1956           space_bdy_ye(i,k,jj) = data3d(i,k,j)
1957        END DO
1958        END DO
1959        END DO
1960     END IF
1961     
1962  END SUBROUTINE stuff_bdy_new
1963  
1964  SUBROUTINE stuff_bdytend_new ( data3dnew , data3dold , time_diff , &
1965                              space_bdy_xs, space_bdy_xe, space_bdy_ys, space_bdy_ye, &
1966                              char_stagger , &
1967                              spec_bdy_width , &
1968                              ids, ide, jds, jde, kds, kde , &
1969                              ims, ime, jms, jme, kms, kme , & 
1970                              its, ite, jts, jte, kts, kte )
1971  
1972  !  This routine puts the tendency data into the proper locations
1973  !  for the lateral boundary arrays.
1974  
1975     USE module_state_description
1976     
1977     IMPLICIT NONE
1978  
1979     INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde
1980     INTEGER , INTENT(IN) :: ims, ime, jms, jme, kms, kme
1981     INTEGER , INTENT(IN) :: its, ite, jts, jte, kts, kte
1982     INTEGER , INTENT(IN) :: spec_bdy_width
1983     REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: data3dnew , data3dold
1984     REAL , DIMENSION(jms:jme,kds:kde,spec_bdy_width) , INTENT(OUT) :: space_bdy_xs, space_bdy_xe
1985     REAL , DIMENSION(ims:ime,kds:kde,spec_bdy_width) , INTENT(OUT) :: space_bdy_ys, space_bdy_ye
1986     CHARACTER (LEN=1) , INTENT(IN) :: char_stagger
1987     REAL , INTENT(IN) :: time_diff ! seconds
1988  
1989     INTEGER :: i , ii , j , jj , k
1990  
1991     !  There are four lateral boundary locations that are stored.
1992  
1993     !  X start boundary
1994  
1995     IF ( char_stagger .EQ. 'W' ) THEN
1996        DO j = MAX(jds,jts) , MIN(jde-1,jte)
1997        DO k = kds , kde
1998        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
1999           space_bdy_xs(j,k,i) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2000        END DO
2001        END DO
2002        END DO
2003     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2004        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2005        DO k = kds , kde
2006        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2007           space_bdy_xs(j,k,i) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2008        END DO
2009        END DO
2010        END DO
2011     ELSE IF ( char_stagger .EQ. 'V' ) THEN
2012        DO j = MAX(jds,jts) , MIN(jde,jte)
2013        DO k = kds , kde - 1
2014        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2015           space_bdy_xs(j,k,i) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2016        END DO
2017        END DO
2018        END DO
2019     ELSE
2020        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2021        DO k = kds , kde - 1
2022        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2023           space_bdy_xs(j,k,i) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2024        END DO
2025        END DO
2026        END DO
2027     END IF
2028  
2029     !  X end boundary
2030  
2031     IF      ( char_stagger .EQ. 'U' ) THEN
2032        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2033        DO k = kds , kde - 1
2034        DO i = MIN(ide,ite) , MAX(ide - spec_bdy_width + 1,its) , -1
2035           ii = ide - i + 1
2036           space_bdy_xe(j,k,ii) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2037        END DO
2038        END DO
2039        END DO
2040     ELSE IF ( char_stagger .EQ. 'V' ) THEN
2041        DO j = MAX(jds,jts) , MIN(jde,jte)
2042        DO k = kds , kde - 1
2043        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2044           ii = ide - i
2045           space_bdy_xe(j,k,ii) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2046        END DO
2047        END DO
2048        END DO
2049     ELSE IF ( char_stagger .EQ. 'W' ) THEN
2050        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2051        DO k = kds , kde
2052        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2053           ii = ide - i
2054           space_bdy_xe(j,k,ii) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2055        END DO
2056        END DO
2057        END DO
2058     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2059        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2060        DO k = kds , kde
2061        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2062           ii = ide - i
2063           space_bdy_xe(j,k,ii) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2064        END DO
2065        END DO
2066        END DO
2067     ELSE
2068        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2069        DO k = kds , kde - 1
2070        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2071           ii = ide - i
2072           space_bdy_xe(j,k,ii) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2073        END DO
2074        END DO
2075        END DO
2076     END IF
2077  
2078     !  Y start boundary
2079  
2080     IF ( char_stagger .EQ. 'W' ) THEN
2081        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2082        DO k = kds , kde
2083        DO i = MAX(ids,its) , MIN(ide-1,ite)
2084           space_bdy_ys(i,k,j) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2085        END DO
2086        END DO
2087        END DO
2088     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2089        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2090        DO k = kds , kde
2091        DO i = MAX(ids,its) , MIN(ide-1,ite)
2092           space_bdy_ys(i,k,j) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2093        END DO
2094        END DO
2095        END DO
2096     ELSE IF ( char_stagger .EQ. 'U' ) THEN
2097        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2098        DO k = kds , kde - 1
2099        DO i = MAX(ids,its) , MIN(ide,ite)
2100           space_bdy_ys(i,k,j) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2101        END DO
2102        END DO
2103        END DO
2104     ELSE
2105        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2106        DO k = kds , kde - 1
2107        DO i = MAX(ids,its) , MIN(ide-1,ite)
2108           space_bdy_ys(i,k,j) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2109        END DO
2110        END DO
2111        END DO
2112     END IF
2113  
2114     !  Y end boundary
2115  
2116     IF      ( char_stagger .EQ. 'V' ) THEN
2117        DO j = MIN(jde,jte) , MAX(jde - spec_bdy_width + 1,jts) , -1
2118        DO k = kds , kde - 1
2119        DO i = MAX(ids,its) , MIN(ide-1,ite)
2120           jj = jde - j + 1
2121           space_bdy_ye(i,k,jj) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2122        END DO
2123        END DO
2124        END DO
2125     ELSE IF ( char_stagger .EQ. 'U' ) THEN
2126        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2127        DO k = kds , kde - 1
2128        DO i = MAX(ids,its) , MIN(ide,ite)
2129           jj = jde - j
2130           space_bdy_ye(i,k,jj) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2131        END DO
2132        END DO
2133        END DO
2134     ELSE IF ( char_stagger .EQ. 'W' ) THEN
2135        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2136        DO k = kds , kde
2137        DO i = MAX(ids,its) , MIN(ide-1,ite)
2138           jj = jde - j
2139           space_bdy_ye(i,k,jj) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2140        END DO
2141        END DO
2142        END DO
2143     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2144        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2145        DO k = kds , kde
2146        DO i = MAX(ids,its) , MIN(ide-1,ite)
2147           jj = jde - j
2148           space_bdy_ye(i,k,jj) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2149        END DO
2150        END DO
2151        END DO
2152     ELSE
2153        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2154        DO k = kds , kde - 1
2155        DO i = MAX(ids,its) , MIN(ide-1,ite)
2156           jj = jde - j
2157           space_bdy_ye(i,k,jj) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2158        END DO
2159        END DO
2160        END DO
2161     END IF
2162     
2163  END SUBROUTINE stuff_bdytend_new
2164 
2165 !--- old versions for use with modules that use the old bdy data structures ---
2166 
2167  SUBROUTINE stuff_bdy_old ( data3d , space_bdy , char_stagger , &
2168                              ijds , ijde , spec_bdy_width , &
2169                              ids, ide, jds, jde, kds, kde , &
2170                              ims, ime, jms, jme, kms, kme , & 
2171                              its, ite, jts, jte, kts, kte )
2172  
2173  !  This routine puts the data in the 3d arrays into the proper locations
2174  !  for the lateral boundary arrays.
2175  
2176     USE module_state_description
2177     
2178     IMPLICIT NONE
2179  
2180     INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde
2181     INTEGER , INTENT(IN) :: ims, ime, jms, jme, kms, kme
2182     INTEGER , INTENT(IN) :: its, ite, jts, jte, kts, kte
2183     INTEGER , INTENT(IN) :: ijds , ijde , spec_bdy_width
2184     REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: data3d
2185     REAL , DIMENSION(ijds:ijde,kds:kde,spec_bdy_width,4) , INTENT(OUT) :: space_bdy
2186     CHARACTER (LEN=1) , INTENT(IN) :: char_stagger
2187  
2188     INTEGER :: i , ii , j , jj , k
2189  
2190     !  There are four lateral boundary locations that are stored.
2191  
2192     !  X start boundary
2193  
2194     IF ( char_stagger .EQ. 'W' ) THEN
2195        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2196        DO k = kds , kde
2197        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2198           space_bdy(j,k,i,P_XSB) = data3d(i,k,j)
2199        END DO
2200        END DO
2201        END DO
2202     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2203        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2204        DO k = kds , kde
2205        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2206           space_bdy(j,k,i,P_XSB) = data3d(i,k,j)
2207        END DO
2208        END DO
2209        END DO
2210     ELSE IF ( char_stagger .EQ. 'V' ) THEN
2211        DO j = MAX(jds,jts) , MIN(jde,jte)
2212        DO k = kds , kde - 1
2213        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2214           space_bdy(j,k,i,P_XSB) = data3d(i,k,j)
2215        END DO
2216        END DO
2217        END DO
2218     ELSE
2219        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2220        DO k = kds , kde - 1
2221        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2222           space_bdy(j,k,i,P_XSB) = data3d(i,k,j)
2223        END DO
2224        END DO
2225        END DO
2226     END IF
2227  
2228     !  X end boundary
2229  
2230     IF      ( char_stagger .EQ. 'U' ) THEN
2231        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2232        DO k = kds , kde - 1
2233        DO i = MIN(ide,ite) , MAX(ide - spec_bdy_width + 1,its) , -1
2234           ii = ide - i + 1
2235           space_bdy(j,k,ii,P_XEB) = data3d(i,k,j)
2236        END DO
2237        END DO
2238        END DO
2239     ELSE IF ( char_stagger .EQ. 'V' ) THEN
2240        DO j = MAX(jds,jts) , MIN(jde,jte)
2241        DO k = kds , kde - 1
2242        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2243           ii = ide - i
2244           space_bdy(j,k,ii,P_XEB) = data3d(i,k,j)
2245        END DO
2246        END DO
2247        END DO
2248     ELSE IF ( char_stagger .EQ. 'W' ) THEN
2249        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2250        DO k = kds , kde
2251        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2252           ii = ide - i
2253           space_bdy(j,k,ii,P_XEB) = data3d(i,k,j)
2254        END DO
2255        END DO
2256        END DO
2257     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2258        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2259        DO k = kds , kde
2260        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2261           ii = ide - i
2262           space_bdy(j,k,ii,P_XEB) = data3d(i,k,j)
2263        END DO
2264        END DO
2265        END DO
2266     ELSE
2267        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2268        DO k = kds , kde - 1
2269        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2270           ii = ide - i
2271           space_bdy(j,k,ii,P_XEB) = data3d(i,k,j)
2272        END DO
2273        END DO
2274        END DO
2275     END IF
2276  
2277     !  Y start boundary
2278  
2279     IF ( char_stagger .EQ. 'W' ) THEN
2280        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2281        DO k = kds , kde
2282        DO i = MAX(ids,its) , MIN(ide-1,ite)
2283           space_bdy(i,k,j,P_YSB) = data3d(i,k,j)
2284        END DO
2285        END DO
2286        END DO
2287     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2288        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2289        DO k = kds , kde
2290        DO i = MAX(ids,its) , MIN(ide-1,ite)
2291           space_bdy(i,k,j,P_YSB) = data3d(i,k,j)
2292        END DO
2293        END DO
2294        END DO
2295     ELSE IF ( char_stagger .EQ. 'U' ) THEN
2296        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2297        DO k = kds , kde - 1
2298        DO i = MAX(ids,its) , MIN(ide,ite)
2299           space_bdy(i,k,j,P_YSB) = data3d(i,k,j)
2300        END DO
2301        END DO
2302        END DO
2303     ELSE
2304        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2305        DO k = kds , kde - 1
2306        DO i = MAX(ids,its) , MIN(ide-1,ite)
2307           space_bdy(i,k,j,P_YSB) = data3d(i,k,j)
2308        END DO
2309        END DO
2310        END DO
2311     END IF
2312  
2313     !  Y end boundary
2314  
2315     IF      ( char_stagger .EQ. 'V' ) THEN
2316        DO j = MIN(jde,jte) , MAX(jde - spec_bdy_width + 1,jts) , -1
2317        DO k = kds , kde - 1
2318        DO i = MAX(ids,its) , MIN(ide-1,ite)
2319           jj = jde - j + 1
2320           space_bdy(i,k,jj,P_YEB) = data3d(i,k,j)
2321        END DO
2322        END DO
2323        END DO
2324     ELSE IF ( char_stagger .EQ. 'U' ) THEN
2325        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2326        DO k = kds , kde - 1
2327        DO i = MAX(ids,its) , MIN(ide,ite)
2328           jj = jde - j
2329           space_bdy(i,k,jj,P_YEB) = data3d(i,k,j)
2330        END DO
2331        END DO
2332        END DO
2333     ELSE IF ( char_stagger .EQ. 'W' ) THEN
2334        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2335        DO k = kds , kde
2336        DO i = MAX(ids,its) , MIN(ide-1,ite)
2337           jj = jde - j
2338           space_bdy(i,k,jj,P_YEB) = data3d(i,k,j)
2339        END DO
2340        END DO
2341        END DO
2342     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2343        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2344        DO k = kds , kde
2345        DO i = MAX(ids,its) , MIN(ide-1,ite)
2346           jj = jde - j
2347           space_bdy(i,k,jj,P_YEB) = data3d(i,k,j)
2348        END DO
2349        END DO
2350        END DO
2351     ELSE
2352        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2353        DO k = kds , kde - 1
2354        DO i = MAX(ids,its) , MIN(ide-1,ite)
2355           jj = jde - j
2356           space_bdy(i,k,jj,P_YEB) = data3d(i,k,j)
2357        END DO
2358        END DO
2359        END DO
2360     END IF
2361     
2362  END SUBROUTINE stuff_bdy_old
2363  
2364  SUBROUTINE stuff_bdytend_old ( data3dnew , data3dold , time_diff , space_bdy , char_stagger , &
2365                              ijds , ijde , spec_bdy_width , &
2366                              ids, ide, jds, jde, kds, kde , &
2367                              ims, ime, jms, jme, kms, kme , & 
2368                              its, ite, jts, jte, kts, kte )
2369  
2370  !  This routine puts the tendency data into the proper locations
2371  !  for the lateral boundary arrays.
2372  
2373     USE module_state_description
2374     
2375     IMPLICIT NONE
2376  
2377     INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde
2378     INTEGER , INTENT(IN) :: ims, ime, jms, jme, kms, kme
2379     INTEGER , INTENT(IN) :: its, ite, jts, jte, kts, kte
2380     INTEGER , INTENT(IN) :: ijds , ijde , spec_bdy_width
2381     REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: data3dnew , data3dold
2382 !    REAL , DIMENSION(:,:,:,:) , INTENT(OUT) :: space_bdy
2383     REAL , DIMENSION(ijds:ijde,kds:kde,spec_bdy_width,4) , INTENT(OUT) :: space_bdy
2384     CHARACTER (LEN=1) , INTENT(IN) :: char_stagger
2385     REAL , INTENT(IN) :: time_diff ! seconds
2386  
2387     INTEGER :: i , ii , j , jj , k
2388  
2389     !  There are four lateral boundary locations that are stored.
2390  
2391     !  X start boundary
2392  
2393     IF ( char_stagger .EQ. 'W' ) THEN
2394        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2395        DO k = kds , kde
2396        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2397           space_bdy(j,k,i,P_XSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2398 !         space_bdy(j,k,i,P_XSB) = 0. ! zeroout
2399        END DO
2400        END DO
2401        END DO
2402     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2403        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2404        DO k = kds , kde
2405        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2406           space_bdy(j,k,i,P_XSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2407 !         space_bdy(j,k,i,P_XSB) = 0. ! zeroout
2408        END DO
2409        END DO
2410        END DO
2411     ELSE IF ( char_stagger .EQ. 'V' ) THEN
2412        DO j = MAX(jds,jts) , MIN(jde,jte)
2413        DO k = kds , kde - 1
2414        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2415           space_bdy(j,k,i,P_XSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2416 !         space_bdy(j,k,i,P_XSB) = 0. ! zeroout
2417        END DO
2418        END DO
2419        END DO
2420     ELSE
2421        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2422        DO k = kds , kde - 1
2423        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2424           space_bdy(j,k,i,P_XSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2425 !         space_bdy(j,k,i,P_XSB) = 0. ! zeroout
2426        END DO
2427        END DO
2428        END DO
2429     END IF
2430  
2431     !  X end boundary
2432  
2433     IF      ( char_stagger .EQ. 'U' ) THEN
2434        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2435        DO k = kds , kde - 1
2436        DO i = MIN(ide,ite) , MAX(ide - spec_bdy_width + 1,its) , -1
2437           ii = ide - i + 1
2438           space_bdy(j,k,ii,P_XEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2439 !         space_bdy(j,k,ii,P_XEB) = 0. ! zeroout
2440        END DO
2441        END DO
2442        END DO
2443     ELSE IF ( char_stagger .EQ. 'V' ) THEN
2444        DO j = MAX(jds,jts) , MIN(jde,jte)
2445        DO k = kds , kde - 1
2446        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2447           ii = ide - i
2448           space_bdy(j,k,ii,P_XEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2449 !         space_bdy(j,k,ii,P_XEB) = 0. ! zeroout
2450        END DO
2451        END DO
2452        END DO
2453     ELSE IF ( char_stagger .EQ. 'W' ) THEN
2454        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2455        DO k = kds , kde
2456        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2457           ii = ide - i
2458           space_bdy(j,k,ii,P_XEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2459 !         space_bdy(j,k,ii,P_XEB) = 0. ! zeroout
2460        END DO
2461        END DO
2462        END DO
2463     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2464        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2465        DO k = kds , kde
2466        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2467           ii = ide - i
2468           space_bdy(j,k,ii,P_XEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2469 !         space_bdy(j,k,ii,P_XEB) = 0. ! zeroout
2470        END DO
2471        END DO
2472        END DO
2473     ELSE
2474        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2475        DO k = kds , kde - 1
2476        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2477           ii = ide - i
2478           space_bdy(j,k,ii,P_XEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2479 !         space_bdy(j,k,ii,P_XEB) = 0. ! zeroout
2480        END DO
2481        END DO
2482        END DO
2483     END IF
2484  
2485     !  Y start boundary
2486  
2487     IF ( char_stagger .EQ. 'W' ) THEN
2488        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2489        DO k = kds , kde
2490        DO i = MAX(ids,its) , MIN(ide-1,ite)
2491           space_bdy(i,k,j,P_YSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2492 !         space_bdy(i,k,j,P_YSB) = 0. ! zeroout
2493        END DO
2494        END DO
2495        END DO
2496     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2497        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2498        DO k = kds , kde
2499        DO i = MAX(ids,its) , MIN(ide-1,ite)
2500           space_bdy(i,k,j,P_YSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2501 !         space_bdy(i,k,j,P_YSB) = 0. ! zeroout
2502        END DO
2503        END DO
2504        END DO
2505     ELSE IF ( char_stagger .EQ. 'U' ) THEN
2506        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2507        DO k = kds , kde - 1
2508        DO i = MAX(ids,its) , MIN(ide,ite)
2509           space_bdy(i,k,j,P_YSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2510 !         space_bdy(i,k,j,P_YSB) = 0. ! zeroout
2511        END DO
2512        END DO
2513        END DO
2514     ELSE
2515        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2516        DO k = kds , kde - 1
2517        DO i = MAX(ids,its) , MIN(ide-1,ite)
2518           space_bdy(i,k,j,P_YSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2519 !         space_bdy(i,k,j,P_YSB) = 0. ! zeroout
2520        END DO
2521        END DO
2522        END DO
2523     END IF
2524  
2525     !  Y end boundary
2526  
2527     IF      ( char_stagger .EQ. 'V' ) THEN
2528        DO j = MIN(jde,jte) , MAX(jde - spec_bdy_width + 1,jts) , -1
2529        DO k = kds , kde - 1
2530        DO i = MAX(ids,its) , MIN(ide-1,ite)
2531           jj = jde - j + 1
2532           space_bdy(i,k,jj,P_YEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2533 !         space_bdy(i,k,jj,P_YEB) = 0. ! zeroout
2534        END DO
2535        END DO
2536        END DO
2537     ELSE IF ( char_stagger .EQ. 'U' ) THEN
2538        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2539        DO k = kds , kde - 1
2540        DO i = MAX(ids,its) , MIN(ide,ite)
2541           jj = jde - j
2542           space_bdy(i,k,jj,P_YEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2543 !         space_bdy(i,k,jj,P_YEB) = 0. ! zeroout
2544        END DO
2545        END DO
2546        END DO
2547     ELSE IF ( char_stagger .EQ. 'W' ) THEN
2548        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2549        DO k = kds , kde
2550        DO i = MAX(ids,its) , MIN(ide-1,ite)
2551           jj = jde - j
2552           space_bdy(i,k,jj,P_YEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2553 !         space_bdy(i,k,jj,P_YEB) = 0. ! zeroout
2554        END DO
2555        END DO
2556        END DO
2557     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2558        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2559        DO k = kds , kde
2560        DO i = MAX(ids,its) , MIN(ide-1,ite)
2561           jj = jde - j
2562           space_bdy(i,k,jj,P_YEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2563 !         space_bdy(i,k,jj,P_YEB) = 0. ! zeroout
2564        END DO
2565        END DO
2566        END DO
2567     ELSE
2568        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2569        DO k = kds , kde - 1
2570        DO i = MAX(ids,its) , MIN(ide-1,ite)
2571           jj = jde - j
2572           space_bdy(i,k,jj,P_YEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2573 !         space_bdy(i,k,jj,P_YEB) = 0. ! zeroout
2574        END DO
2575        END DO
2576        END DO
2577     END IF
2578     
2579  END SUBROUTINE stuff_bdytend_old
2580 
2581 END MODULE module_bc
2582 
2583 SUBROUTINE get_bdyzone_x ( bzx )
2584   USE module_bc
2585   IMPLICIT NONE
2586   INTEGER bzx
2587   bzx = bdyzone_x
2588 END SUBROUTINE get_bdyzone_x
2589 
2590 SUBROUTINE get_bdyzone_y ( bzy)
2591   USE module_bc
2592   IMPLICIT NONE
2593   INTEGER bzy
2594   bzy = bdyzone_y
2595 END SUBROUTINE get_bdyzone_y
2596 
2597 SUBROUTINE get_bdyzone ( bz)
2598   USE module_bc
2599   IMPLICIT NONE
2600   INTEGER bz
2601   bz = bdyzone
2602 END SUBROUTINE get_bdyzone
2603