module_bc.F

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