updated top-level README and version_decl for V4.5 (#1847)
[WRF.git] / share / module_bc.F
blob473f9343a6c59a0aeecb841fa3608200766e93c4
2 !WRF:MODEL_LAYER:BOUNDARY
5 MODULE module_bc
7    USE module_configure
8    USE module_wrf_error
9    USE module_model_constants
10    IMPLICIT NONE
12 !   TYPE bcs
14 !     LOGICAL                     :: periodic_x
15 !     LOGICAL                     :: symmetric_xs
16 !     LOGICAL                     :: symmetric_xe
17 !     LOGICAL                     :: open_xs
18 !     LOGICAL                     :: open_xe
19 !     LOGICAL                     :: periodic_y
20 !     LOGICAL                     :: symmetric_ys
21 !     LOGICAL                     :: symmetric_ye
22 !     LOGICAL                     :: open_ys
23 !     LOGICAL                     :: open_ye
24 !     LOGICAL                     :: nested
25 !     LOGICAL                     :: specified
26 !     LOGICAL                     :: top_radiation
28 !   END TYPE bcs
30 !  set the bdyzone.  We are hardwiring this here and we'll
31 !  decide later where it should be set and stored
33    INTEGER, PARAMETER            :: bdyzone = 4
34    INTEGER, PARAMETER            :: bdyzone_x = bdyzone
35    INTEGER, PARAMETER            :: bdyzone_y = bdyzone
37    INTERFACE stuff_bdy
38      MODULE PROCEDURE stuff_bdy_new , stuff_bdy_old
39    END INTERFACE
41    INTERFACE stuff_bdytend
42      MODULE PROCEDURE stuff_bdytend_new , stuff_bdytend_old
43    END INTERFACE
45 CONTAINS
47   SUBROUTINE boundary_condition_check ( config_flags, bzone, error, gn )
49 !  this routine checks the boundary condition logicals
50 !  to make sure that the boundary conditions are not over
51 !  or under specified.  The routine also checks that the
52 !  boundary zone is sufficiently sized for the specified
53 !  boundary conditions
55   IMPLICIT NONE
57   TYPE( grid_config_rec_type ) config_flags
59   INTEGER, INTENT(IN   ) :: bzone, gn
60   INTEGER, INTENT(INOUT) :: error
62 ! local variables
64   INTEGER :: xs_bc, xe_bc, ys_bc, ye_bc, bzone_min
65   INTEGER :: nprocx, nprocy
67   CALL wrf_debug( 100 , ' checking boundary conditions for grid ' )
69   error = 0
70   xs_bc = 0
71   xe_bc = 0
72   ys_bc = 0
73   ye_bc = 0
75 !  sum the number of conditions specified for each lateral boundary.
76 !  obviously, this number should be 1
78   IF( config_flags%periodic_x ) THEN
79     xs_bc = xs_bc+1
80     xe_bc = xe_bc+1
81   ENDIF
83   IF( config_flags%periodic_y ) THEN
84     ys_bc = ys_bc+1
85     ye_bc = ye_bc+1
86   ENDIF
88   IF( config_flags%symmetric_xs ) xs_bc = xs_bc + 1
89   IF( config_flags%symmetric_xe ) xe_bc = xe_bc + 1
90   IF( config_flags%open_xs )      xs_bc = xs_bc + 1
91   IF( config_flags%open_xe )      xe_bc = xe_bc + 1
94   IF( config_flags%symmetric_ys ) ys_bc = ys_bc + 1
95   IF( config_flags%symmetric_ye ) ye_bc = ye_bc + 1
96   IF( config_flags%open_ys )      ys_bc = ys_bc + 1
97   IF( config_flags%open_ye )      ye_bc = ye_bc + 1
99   IF( config_flags%nested ) THEN
100      xs_bc = xs_bc + 1
101      xe_bc = xe_bc + 1
102      ys_bc = ys_bc + 1
103      ye_bc = ye_bc + 1
104    ENDIF
106   IF( config_flags%specified ) THEN
107      IF( .NOT. config_flags%periodic_x)xs_bc = xs_bc + 1
108      IF( .NOT. config_flags%periodic_x)xe_bc = xe_bc + 1
109      ys_bc = ys_bc + 1
110      ye_bc = ye_bc + 1
111    ENDIF
113   IF( config_flags%polar ) THEN
114      ys_bc = ys_bc + 1
115      ye_bc = ye_bc + 1
116    ENDIF
118 !  check the number of conditions for each boundary
120    IF( (xs_bc /= 1) .or. &
121        (xe_bc /= 1) .or. &
122        (ys_bc /= 1) .or. &
123        (ye_bc /= 1)         ) THEN
125      error = 1
127      write( wrf_err_message ,*) ' *** Error in boundary condition specification '
128      CALL wrf_message ( wrf_err_message )
129      write( wrf_err_message ,*) ' boundary conditions at xs ', xs_bc
130      CALL wrf_message ( wrf_err_message )
131      write( wrf_err_message ,*) ' boundary conditions at xe ', xe_bc
132      CALL wrf_message ( wrf_err_message )
133      write( wrf_err_message ,*) ' boundary conditions at ys ', ys_bc
134      CALL wrf_message ( wrf_err_message )
135      write( wrf_err_message ,*) ' boundary conditions at ye ', ye_bc
136      CALL wrf_message ( wrf_err_message )
137      write( wrf_err_message ,*) ' boundary conditions logicals are '
138      CALL wrf_message ( wrf_err_message )
139      write( wrf_err_message ,*) ' periodic_x   ',config_flags%periodic_x
140      CALL wrf_message ( wrf_err_message )
141      write( wrf_err_message ,*) ' periodic_y   ',config_flags%periodic_y
142      CALL wrf_message ( wrf_err_message )
143      write( wrf_err_message ,*) ' symmetric_xs ',config_flags%symmetric_xs
144      CALL wrf_message ( wrf_err_message )
145      write( wrf_err_message ,*) ' symmetric_xe ',config_flags%symmetric_xe
146      CALL wrf_message ( wrf_err_message )
147      write( wrf_err_message ,*) ' symmetric_ys ',config_flags%symmetric_ys
148      CALL wrf_message ( wrf_err_message )
149      write( wrf_err_message ,*) ' symmetric_ye ',config_flags%symmetric_ye
150      CALL wrf_message ( wrf_err_message )
151      write( wrf_err_message ,*) ' open_xs      ',config_flags%open_xs
152      CALL wrf_message ( wrf_err_message )
153      write( wrf_err_message ,*) ' open_xe      ',config_flags%open_xe
154      CALL wrf_message ( wrf_err_message )
155      write( wrf_err_message ,*) ' open_ys      ',config_flags%open_ys
156      CALL wrf_message ( wrf_err_message )
157      write( wrf_err_message ,*) ' open_ye      ',config_flags%open_ye
158      CALL wrf_message ( wrf_err_message )
159      write( wrf_err_message ,*) ' polar        ',config_flags%polar
160      CALL wrf_message ( wrf_err_message )
161      write( wrf_err_message ,*) ' nested       ',config_flags%nested
162      CALL wrf_message ( wrf_err_message )
163      write( wrf_err_message ,*) ' specified    ',config_flags%specified
164      CALL wrf_message ( wrf_err_message )
165      CALL wrf_error_fatal( ' *** Error in boundary condition specification ' )
167    ENDIF
169 !  now check to see if boundary zone size is sufficient.
170 !  we could have the necessary boundary zone size be returned
171 !  to the calling routine.
173    IF( config_flags%periodic_x   .or. &
174        config_flags%periodic_y   .or. &
175        config_flags%symmetric_xs .or. &
176        config_flags%symmetric_xe .or. &
177        config_flags%symmetric_ys .or. &
178        config_flags%symmetric_ye        )  THEN
180        bzone_min = MAX( 1,                                  &
181                         (config_flags%h_mom_adv_order+1)/2, &
182                         (config_flags%h_sca_adv_order+1)/2 )
184        IF( bzone < bzone_min) THEN  
186          error = 2
187          WRITE ( wrf_err_message , * ) ' boundary zone not large enough '
188          CALL wrf_message ( wrf_err_message )
189          WRITE ( wrf_err_message , * ) ' boundary zone specified      ',bzone
190          CALL wrf_message ( wrf_err_message )
191          WRITE ( wrf_err_message , * ) ' minimum boundary zone needed ',bzone_min
192          CALL wrf_error_fatal ( wrf_err_message )
194        ENDIF
195    ENDIF
197    CALL wrf_debug ( 100 , ' boundary conditions OK for grid ' )
199    END subroutine boundary_condition_check
201 !--------------------------------------------------------------------------
202    SUBROUTINE set_physical_bc2d( dat, variable_in,  &
203                                  config_flags,           &
204                                  ids,ide, jds,jde,   & ! domain dims
205                                  ims,ime, jms,jme,   & ! memory dims
206                                  ips,ipe, jps,jpe,   & ! patch  dims
207                                  its,ite, jts,jte   )      
209 !  This subroutine sets the data in the boundary region, by direct
210 !  assignment if possible, for periodic and symmetric (wall)
211 !  boundary conditions.  Currently, we are only doing 1 variable
212 !  at a time - lots of overhead, so maybe this routine can be easily
213 !  inlined later or we could pass multiple variables -
214 !  would probably want a largestep and smallstep version.
216 !  15 Jan 99, Dave
217 !  Modified the incoming its,ite,jts,jte to truly be the tile size.
218 !  This required modifying the loop limits when the "istag" or "jstag"
219 !  is used, as this is only required at the end of the domain.
221       IMPLICIT NONE
223       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde
224       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme
225       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe
226       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte
227       CHARACTER,    INTENT(IN   )    :: variable_in
229       CHARACTER                      :: variable
231       REAL,  DIMENSION( ims:ime , jms:jme ) :: dat
232       TYPE( grid_config_rec_type ) config_flags
234       INTEGER  :: i, j, istag, jstag, itime, istart, iend
236       LOGICAL  :: debug, open_bc_copy
238 !------------
240       debug = .false.
242       open_bc_copy = .false.
244       variable = variable_in
245       IF ( variable_in .ge. 'A' .and. variable_in .le. 'Z' ) THEN
246         variable = CHAR( ICHAR(variable_in) - ICHAR('A') + ICHAR('a') )
247       ENDIF
248       IF ((variable == 'u') .or. (variable == 'v') .or.  &
249           (variable == 'w') .or. (variable == 't') .or.  &
250           (variable == 'x') .or. (variable == 'y') .or.  &
251           (variable == 'r') .or. (variable == 'p') ) open_bc_copy = .true.
253 !  begin, first set a staggering variable
255       istag = -1
256       jstag = -1
258       IF ((variable == 'u') .or. (variable == 'x')) istag = 0
259       IF ((variable == 'v') .or. (variable == 'y')) jstag = 0
261       if(debug) then
262         write(6,*) ' in bc2d, var is ',variable, istag, jstag
263         write(6,*) ' b.cs are ',  &
264       config_flags%periodic_x,  &
265       config_flags%periodic_y
266       end if
267       
268       IF ( variable == 'd' ) then  !JDM
269          istag = 0
270          jstag = 0
271       ENDIF
272       IF ( variable == 'e' ) then  !JDM
273          istag = 0
274       ENDIF
275       IF ( variable == 'f' ) then  !JDM
276          jstag = 0
277       ENDIF
279 !  periodic conditions.
280 !  note, patch must cover full range in periodic dir, or else
281 !  its intra-patch communication that is handled elsewheres.
282 !  symmetry conditions can always be handled here, because no
283 !  outside patch communication is needed
285       periodicity_x:  IF( ( config_flags%periodic_x ) ) THEN
286         IF ( ( ids == ips ) .and.  ( ide == ipe ) ) THEN  ! test if east and west both on-processor
287           IF ( its == ids ) THEN
289             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
290             DO i = 0,-(bdyzone-1),-1
291               dat(ids+i-1,j) = dat(ide+i-1,j)
292             ENDDO
293             ENDDO
295           ENDIF
297           IF ( ite == ide ) THEN
299             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
300 !!          DO i = 1 , bdyzone
301             DO i = -istag , bdyzone
302               dat(ide+i+istag,j) = dat(ids+i+istag,j)
303             ENDDO
304             ENDDO
306           ENDIF
307         ENDIF
309       ELSE
311         symmetry_xs: IF( ( config_flags%symmetric_xs ) .and.  &
312                          ( its == ids )                  )  THEN
314           IF ( (variable /= 'u') .and. (variable /= 'x') ) THEN
316             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
317             DO i = 1, bdyzone
318               dat(ids-i,j) = dat(ids+i-1,j) !  here, dat(0) = dat(1), etc
319             ENDDO                             !  symmetry about dat(0.5) (u=0 pt)
320             ENDDO
322           ELSE
324             IF( variable == 'u' ) THEN
326               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
327               DO i = 0, bdyzone-1
328                 dat(ids-i,j) = - dat(ids+i,j) ! here, u(0) = - u(2), etc
329               ENDDO                             !  normal b.c symmetry at u(1)
330               ENDDO
332             ELSE
334               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
335               DO i = 0, bdyzone-1
336                 dat(ids-i,j) =   dat(ids+i,j) ! here, phi(0) = phi(2), etc
337               ENDDO                             !  normal b.c symmetry at phi(1)
338               ENDDO
340             END IF
342           ENDIF
344         ENDIF symmetry_xs
347 !  now the symmetry boundary at xe
349         symmetry_xe: IF( ( config_flags%symmetric_xe ) .and.  &
350                          ( ite == ide )                  )  THEN
352           IF ( (variable /= 'u') .and. (variable /= 'x') ) THEN
354             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
355             DO i = 1, bdyzone
356               dat(ide+i-1,j) = dat(ide-i,j)  !  sym. about dat(ide-0.5)
357             ENDDO
358             ENDDO
360           ELSE
362             IF (variable == 'u' ) THEN
364               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
365               DO i = 0, bdyzone-1
366                 dat(ide+i,j) = - dat(ide-i,j)  ! u(ide+1) = - u(ide-1), etc.
367               ENDDO
368               ENDDO
371             ELSE
373               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
374               DO i = 0, bdyzone-1
375                 dat(ide+i,j) = dat(ide-i,j)  !  phi(ide+1) = phi(ide-1), etc.
376               ENDDO
377               ENDDO
379             END IF
381           END IF
383         END IF symmetry_xe
386 !  set open b.c in X copy into boundary zone here.  WCS, 19 March 2000
388         open_xs: IF( ( config_flags%open_xs   .or. &
389                        config_flags%specified .or. &
390                        config_flags%nested            ) .and.  &
391                          ( its == ids ) .and. open_bc_copy  )  THEN
393             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
394               dat(ids-1,j) = dat(ids,j) !  here, dat(0) = dat(1)
395               dat(ids-2,j) = dat(ids,j)
396               dat(ids-3,j) = dat(ids,j)
397             ENDDO
399         ENDIF open_xs
402 !  now the open boundary copy at xe
404         open_xe: IF( ( config_flags%open_xe   .or. &
405                        config_flags%specified .or. &
406                        config_flags%nested            ) .and.  &
407                           ( ite == ide ) .and. open_bc_copy  )  THEN
409           IF ( variable /= 'u' .and. variable /= 'x') THEN
411             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
412               dat(ide  ,j) = dat(ide-1,j)
413               dat(ide+1,j) = dat(ide-1,j)
414               dat(ide+2,j) = dat(ide-1,j)
415             ENDDO
417           ELSE
419             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
420               dat(ide+1,j) = dat(ide,j)
421               dat(ide+2,j) = dat(ide,j)
422               dat(ide+3,j) = dat(ide,j)
423             ENDDO
425           END IF
427         END IF open_xe
429 !  end open b.c in X copy into boundary zone addition.  WCS, 19 March 2000
431       END IF periodicity_x
433 !  same procedure in y
435 !  Set the starting and ending loop indexes in the 'i' direction, so that
436 !  halo cells on the edge of the domain are also updated.  Begin with a default
437 !  start and end index for inner tiles, and then modify if the tile is on the
438 !  edge of the domain.  
440         istart = MAX(ids, its-1)
441         iend = MIN(ite+1, ide+istag)
442         IF ( its .eq. ids) THEN
443           istart = ims
444         END IF
445         IF ( ite .eq. ide) THEN
446           iend = ime
447         END IF
449       periodicity_y:  IF( ( config_flags%periodic_y ) ) THEN
450         IF ( ( jds == jps ) .and. ( jde == jpe ) )  THEN    ! test of both north and south on processor
452           IF( jts == jds ) then
454             DO j = 0, -(bdyzone-1), -1
455               !DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
456               DO i = istart, iend
457                 dat(i,jds+j-1) = dat(i,jde+j-1)
458               ENDDO
459             ENDDO
461           END IF
463           IF( jte == jde ) then
465             DO j = -jstag, bdyzone
466               !DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
467               DO i = istart, iend
468                 dat(i,jde+j+jstag) = dat(i,jds+j+jstag)
469               ENDDO
470             ENDDO
472           END IF
474         END IF
476       ELSE
478         symmetry_ys: IF( ( config_flags%symmetric_ys ) .and.  &
479                          ( jts == jds)                   )  THEN
481           IF ( (variable /= 'v') .and. (variable /= 'y') ) THEN
483             DO j = 1, bdyzone
484               !DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
485               DO i = istart, iend
486                 dat(i,jds-j) = dat(i,jds+j-1)
487               ENDDO
488             ENDDO
490           ELSE
492             IF (variable == 'v') THEN
494               DO j = 1, bdyzone
495                 !DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
496                 DO i = istart, iend
497                   dat(i,jds-j) = - dat(i,jds+j)
498                 ENDDO              
499               ENDDO
501             ELSE
503               DO j = 1, bdyzone
504                 !DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
505                 DO i = istart, iend
506                   dat(i,jds-j) = dat(i,jds+j)
507                 ENDDO              
508               ENDDO
510             END IF
512           ENDIF
514         ENDIF symmetry_ys
516 !  now the symmetry boundary at ye
518         symmetry_ye: IF( ( config_flags%symmetric_ye ) .and.  &
519                          ( jte == jde )                  )  THEN
521           IF ( (variable /= 'v') .and. (variable /= 'y') ) THEN
523             DO j = 1, bdyzone
524               !DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
525               DO i = istart, iend
526                 dat(i,jde+j-1) = dat(i,jde-j)
527               ENDDO                               
528             ENDDO
530           ELSE
532             IF (variable == 'v' ) THEN
534               DO j = 1, bdyzone
535                 !DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
536                 DO i = istart, iend
537                   dat(i,jde+j) = - dat(i,jde-j)    ! bugfix: changed jds on rhs to jde , JM 20020410
538                 ENDDO                               
539               ENDDO
541             ELSE
543               DO j = 1, bdyzone
544                 !DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
545                 DO i = istart, iend
546                   dat(i,jde+j) = dat(i,jde-j)
547                 ENDDO                               
548               ENDDO
550             END IF
552           ENDIF
554         END IF symmetry_ye
556 !  set open b.c in Y copy into boundary zone here.  WCS, 19 March 2000
558         open_ys: IF( ( config_flags%open_ys   .or. &
559                        config_flags%polar     .or. &
560                        config_flags%specified .or. &
561                        config_flags%nested            ) .and.  &
562                          ( jts == jds) .and. open_bc_copy )  THEN
564             !DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
565             DO i = istart, iend
566               dat(i,jds-1) = dat(i,jds)
567               dat(i,jds-2) = dat(i,jds)
568               dat(i,jds-3) = dat(i,jds)
569             ENDDO
571         ENDIF open_ys
573 !  now the open boundary copy at ye
575         open_ye: IF( ( config_flags%open_ye   .or. &
576                        config_flags%polar     .or. &
577                        config_flags%specified .or. &
578                        config_flags%nested            ) .and.  &
579                          ( jte == jde ) .and. open_bc_copy )  THEN
581           IF  (variable /= 'v' .and. variable /= 'y' ) THEN
583             !DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
584             DO i = istart, iend
585               dat(i,jde  ) = dat(i,jde-1)
586               dat(i,jde+1) = dat(i,jde-1)
587               dat(i,jde+2) = dat(i,jde-1)
588             ENDDO                               
590           ELSE
592             !DO i = MAX(ids,its-1), MIN(ite+1,ide+istag)
593             DO i = istart, iend
594               dat(i,jde+1) = dat(i,jde)
595               dat(i,jde+2) = dat(i,jde)
596               dat(i,jde+3) = dat(i,jde)
597             ENDDO                               
599           ENDIF
601         END IF open_ye
602       
603 !  end open b.c in Y copy into boundary zone addition.  WCS, 19 March 2000
605       END IF periodicity_y
607 !  fix corners for doubly periodic domains
609       IF ( config_flags%periodic_x .and. config_flags%periodic_y &
610            .and. (ids == ips) .and. (ide == ipe)                 &
611            .and. (jds == jps) .and. (jde == jpe)                   ) THEN
613          IF ( (its == ids) .and. (jts == jds) ) THEN  ! lower left corner fill
614            DO j = 0, -(bdyzone-1), -1
615            DO i = 0, -(bdyzone-1), -1
616              dat(ids+i-1,jds+j-1) = dat(ide+i-1,jde+j-1)
617            ENDDO
618            ENDDO
619          END IF
621          IF ( (ite == ide) .and. (jts == jds) ) THEN  ! lower right corner fill
622            DO j = 0, -(bdyzone-1), -1
623            DO i = 1, bdyzone
624              dat(ide+i+istag,jds+j-1) = dat(ids+i+istag,jde+j-1)
625            ENDDO
626            ENDDO
627          END IF
629          IF ( (ite == ide) .and. (jte == jde) ) THEN  ! upper right corner fill
630            DO j = 1, bdyzone
631            DO i = 1, bdyzone
632              dat(ide+i+istag,jde+j+jstag) = dat(ids+i+istag,jds+j+jstag)
633            ENDDO
634            ENDDO
635          END IF
637          IF ( (its == ids) .and. (jte == jde) ) THEN  ! upper left corner fill
638            DO j = 1, bdyzone
639            DO i = 0, -(bdyzone-1), -1
640              dat(ids+i-1,jde+j+jstag) = dat(ide+i-1,jds+j+jstag)
641            ENDDO
642            ENDDO
643          END IF
645        END IF
647    END SUBROUTINE set_physical_bc2d
649 !-----------------------------------
651    SUBROUTINE set_physical_bc3d( dat, variable_in,        &
652                                config_flags,                   &
653                                ids,ide, jds,jde, kds,kde,  & ! domain dims
654                                ims,ime, jms,jme, kms,kme,  & ! memory dims
655                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
656                                its,ite, jts,jte, kts,kte )
658 !  This subroutine sets the data in the boundary region, by direct
659 !  assignment if possible, for periodic and symmetric (wall)
660 !  boundary conditions.  Currently, we are only doing 1 variable
661 !  at a time - lots of overhead, so maybe this routine can be easily
662 !  inlined later or we could pass multiple variables -
663 !  would probably want a largestep and smallstep version.
665 !  15 Jan 99, Dave
666 !  Modified the incoming its,ite,jts,jte to truly be the tile size.
667 !  This required modifying the loop limits when the "istag" or "jstag"
668 !  is used, as this is only required at the end of the domain.
670       IMPLICIT NONE
672       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
673       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
674       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
675       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
676       CHARACTER,    INTENT(IN   )    :: variable_in
678       CHARACTER                      :: variable
680       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ) :: dat
681       TYPE( grid_config_rec_type ) config_flags
683       INTEGER  :: i, j, k, istag, jstag, itime, k_end, &
684                   i_start, i_end
686       LOGICAL  :: debug, open_bc_copy
688 !------------
690       debug = .false.
692       open_bc_copy = .false.
694       variable = variable_in
695       IF ( variable_in .ge. 'A' .and. variable_in .le. 'Z' ) THEN
696         variable = CHAR( ICHAR(variable_in) - ICHAR('A') + ICHAR('a') )
697       ENDIF
699       IF ((variable == 'u') .or. (variable == 'v') .or.     &
700           (variable == 'w') .or. (variable == 't') .or.     &
701           (variable == 'd') .or. (variable == 'e') .or. &
702           (variable == 'x') .or. (variable == 'y') .or. &
703           (variable == 'f') .or. (variable == 'r') .or. &
704           (variable == 'p')                        ) open_bc_copy = .true.
706 !  begin, first set a staggering variable
708       istag = -1
709       jstag = -1
710       k_end = max(1,min(kde-1,kte))
713       IF ((variable == 'u') .or. (variable == 'x')) istag = 0
714       IF ((variable == 'v') .or. (variable == 'y')) jstag = 0
715       IF ((variable == 'd') .or. (variable == 'xy')) then
716          istag = 0
717          jstag = 0
718       ENDIF
719       IF ((variable == 'e') ) then
720          istag = 0
721          k_end = min(kde,kte)
722       ENDIF
724       IF ((variable == 'f') ) then
725          jstag = 0
726          k_end = min(kde,kte)
727       ENDIF
729       IF ( variable == 'w')  k_end = min(kde,kte)
731 !      k_end = kte
733       if(debug) then
734         write(6,*) ' in bc, var is ',variable, istag, jstag, kte, k_end
735         write(6,*) ' b.cs are ',  &
736       config_flags%periodic_x,  &
737       config_flags%periodic_y
738       end if
739       
742 !  periodic conditions.
743 !  note, patch must cover full range in periodic dir, or else
744 !  its intra-patch communication that is handled elsewheres.
745 !  symmetry conditions can always be handled here, because no
746 !  outside patch communication is needed
748       periodicity_x:  IF( ( config_flags%periodic_x ) ) THEN
750         IF ( ( ids == ips ) .and. ( ide == ipe ) ) THEN  ! test if both east and west on-processor
751           IF ( its == ids ) THEN
753             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
754             DO k = kts, k_end
755             DO i = 0,-(bdyzone-1),-1
756               dat(ids+i-1,k,j) = dat(ide+i-1,k,j)
757             ENDDO
758             ENDDO
759             ENDDO
761           ENDIF
764           IF ( ite == ide ) THEN
766             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
767             DO k = kts, k_end
768             DO i = -istag , bdyzone
769               dat(ide+i+istag,k,j) = dat(ids+i+istag,k,j)
770             ENDDO
771             ENDDO
772             ENDDO
774           ENDIF
776         ENDIF
778       ELSE
780         symmetry_xs: IF( ( config_flags%symmetric_xs ) .and.  &
781                          ( its == ids )                  )  THEN
783           IF ( istag == -1 ) THEN
785             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
786             DO k = kts, k_end
787             DO i = 1, bdyzone
788               dat(ids-i,k,j) = dat(ids+i-1,k,j) !  here, dat(0) = dat(1), etc
789             ENDDO                                 !  symmetry about dat(0.5) (u = 0 pt)
790             ENDDO
791             ENDDO
793           ELSE
795             IF ( variable == 'u' ) THEN
797               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
798               DO k = kts, k_end
799               DO i = 1, bdyzone
800                 dat(ids-i,k,j) = - dat(ids+i,k,j) ! here, u(0) = - u(2), etc
801               ENDDO                                 !  normal b.c symmetry at u(1)
802               ENDDO
803               ENDDO
805             ELSE
807               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
808               DO k = kts, k_end
809               DO i = 1, bdyzone
810                 dat(ids-i,k,j) = dat(ids+i,k,j) ! here, phi(0) = phi(2), etc
811               ENDDO                               !  normal b.c symmetry at phi(1)
812               ENDDO
813               ENDDO
815             END IF
817           ENDIF
819         ENDIF symmetry_xs
822 !  now the symmetry boundary at xe
824         symmetry_xe: IF( ( config_flags%symmetric_xe ) .and.  &
825                          ( ite == ide )                  )  THEN
827           IF ( istag == -1 ) THEN
829             DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
830             DO k = kts, k_end
831             DO i = 1, bdyzone
832               dat(ide+i-1,k,j) = dat(ide-i,k,j)  !  sym. about dat(ide-0.5)
833             ENDDO
834             ENDDO
835             ENDDO
837           ELSE
839             IF (variable == 'u') THEN
841               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
842               DO k = kts, k_end
843               DO i = 1, bdyzone
844                 dat(ide+i,k,j) = - dat(ide-i,k,j)  ! u(ide+1) = - u(ide-1), etc.
845               ENDDO
846               ENDDO
847               ENDDO
849             ELSE
851               DO j = MAX(jds,jts-1), MIN(jte+1,jde+jstag)
852               DO k = kts, k_end
853               DO i = 1, bdyzone
854                 dat(ide+i,k,j) = dat(ide-i,k,j)  ! phi(ide+1) = - phi(ide-1), etc.
855               ENDDO
856               ENDDO
857               ENDDO
859              END IF
861           END IF
863         END IF symmetry_xe
865 !  set open b.c in X copy into boundary zone here.  WCS, 19 March 2000
867         open_xs: IF( ( config_flags%open_xs   .or. &
868                        config_flags%specified .or. &
869                        config_flags%nested            ) .and.  &
870                          ( its == ids ) .and. open_bc_copy  )  THEN
872             DO j = jts-bdyzone, MIN(jte,jde+jstag)+bdyzone
873 #if defined(INTEL_ALIGN64)
874 !DEC$ ASSUME_ALIGNED dat:64
875 !DEC$ UNROLL(4)
876 !DEC$ IVDEP
877 #endif
878             DO k = kts, k_end
879               dat(ids-1,k,j) = dat(ids,k,j) !  here, dat(0) = dat(1), etc
880               dat(ids-2,k,j) = dat(ids,k,j)
881               dat(ids-3,k,j) = dat(ids,k,j)
882             ENDDO
883             ENDDO
885         ENDIF open_xs
888 !  now the open_xe boundary copy
890         open_xe: IF( ( config_flags%open_xe   .or. &
891                        config_flags%specified .or. &
892                        config_flags%nested            ) .and.  &
893                          ( ite == ide ) .and. open_bc_copy )  THEN
895           IF (variable /= 'u' .and. variable /= 'x' ) THEN
897             DO j = jts-bdyzone, MIN(jte,jde+jstag)+bdyzone
898 #if defined(INTEL_ALIGN64)
899 !DEC$ ASSUME_ALIGNED dat:64
900 !DEC$ UNROLL(4)
901 !DEC$ IVDEP
902 #endif
903             DO k = kts, k_end
904               dat(ide  ,k,j) = dat(ide-1,k,j)
905               dat(ide+1,k,j) = dat(ide-1,k,j)
906               dat(ide+2,k,j) = dat(ide-1,k,j)
907             ENDDO
908             ENDDO
910           ELSE
912 !!!!!!! I am not sure about this one!  JM 20020402
913             DO j = MAX(jds,jts-1)-bdyzone, MIN(jte+1,jde+jstag)+bdyzone
914             DO k = kts, k_end
915               dat(ide+1,k,j) = dat(ide,k,j)
916               dat(ide+2,k,j) = dat(ide,k,j)
917               dat(ide+3,k,j) = dat(ide,k,j)
918             ENDDO
919             ENDDO
921           END IF
923         END IF open_xe
925 !  end open b.c in X copy into boundary zone addition.  WCS, 19 March 2000
927       END IF periodicity_x
929 !  same procedure in y
931 !  Set the starting and ending loop indexes in the 'i' direction, so that
932 !  halo cells on the edge of the domain are also updated.  Begin with a default
933 !  start and end index for inner tiles, and then modify if the tile is on the
934 !  edge of the domain.
936       i_start = MAX(ids, its-1)
937       i_end = MIN(ite+1, ide+istag)
938       IF ( its .eq. ids) THEN
939         i_start = ims
940       END IF
941       IF ( ite .eq. ide) THEN
942         i_end = ime
943       END IF
945       periodicity_y:  IF( ( config_flags%periodic_y ) ) THEN
946         IF ( ( jds == jps ) .and. ( jde == jpe ) )  THEN      ! test if both north and south on processor
947           IF( jts == jds ) then
949             DO j = 0, -(bdyzone-1), -1
950             DO k = kts, k_end
951             DO i = i_start, i_end
952               dat(i,k,jds+j-1) = dat(i,k,jde+j-1)
953             ENDDO
954             ENDDO
955             ENDDO
957           END IF
959           IF( jte == jde ) then
961             DO j = -jstag, bdyzone
962             DO k = kts, k_end
963             DO i = i_start, i_end
964               dat(i,k,jde+j+jstag) = dat(i,k,jds+j+jstag)
965             ENDDO
966             ENDDO
967             ENDDO
969           END IF
971         END IF
973       ELSE
975         symmetry_ys: IF( ( config_flags%symmetric_ys ) .and.  &
976                          ( jts == jds)                   )  THEN
978           IF ( jstag == -1 ) THEN
980             DO j = 1, bdyzone
981             DO k = kts, k_end
982             DO i = i_start, i_end
983               dat(i,k,jds-j) = dat(i,k,jds+j-1)
984             ENDDO                               
985             ENDDO
986             ENDDO
988           ELSE
990             IF (variable == 'v') THEN
992               DO j = 1, bdyzone
993               DO k = kts, k_end
994               DO i = i_start, i_end
995                 dat(i,k,jds-j) = - dat(i,k,jds+j)
996               ENDDO              
997               ENDDO
998               ENDDO
1000             ELSE
1002               DO j = 1, bdyzone
1003               DO k = kts, k_end
1004               DO i = i_start, i_end
1005                 dat(i,k,jds-j) = dat(i,k,jds+j)
1006               ENDDO              
1007               ENDDO
1008               ENDDO
1010             END IF
1012           ENDIF
1014         ENDIF symmetry_ys
1016 !  now the symmetry boundary at ye
1018         symmetry_ye: IF( ( config_flags%symmetric_ye ) .and.  &
1019                          ( jte == jde )                  )  THEN
1021           IF ( jstag == -1 ) THEN
1023             DO j = 1, bdyzone
1024             DO k = kts, k_end
1025             DO i = i_start, i_end
1026               dat(i,k,jde+j-1) = dat(i,k,jde-j)
1027             ENDDO                               
1028             ENDDO
1029             ENDDO
1031           ELSE
1033             IF ( variable == 'v' ) THEN
1035               DO j = 1, bdyzone
1036               DO k = kts, k_end
1037               DO i = i_start, i_end
1038                 dat(i,k,jde+j) = - dat(i,k,jde-j)
1039               ENDDO                               
1040               ENDDO
1041               ENDDO
1043             ELSE
1045               DO j = 1, bdyzone
1046               DO k = kts, k_end
1047               DO i = i_start, i_end
1048                 dat(i,k,jde+j) = dat(i,k,jde-j)
1049               ENDDO                               
1050               ENDDO
1051               ENDDO
1053             END IF
1055           ENDIF
1057         END IF symmetry_ye
1058       
1059 !  set open b.c in Y copy into boundary zone here.  WCS, 19 March 2000
1061         open_ys: IF( ( config_flags%open_ys   .or. &
1062                        config_flags%polar     .or. &
1063                        config_flags%specified .or. &
1064                        config_flags%nested            ) .and.  &
1065                          ( jts == jds) .and. open_bc_copy )  THEN
1067             DO k = kts, k_end
1068             DO i = i_start, i_end
1069               dat(i,k,jds-1) = dat(i,k,jds)
1070               dat(i,k,jds-2) = dat(i,k,jds)
1071               dat(i,k,jds-3) = dat(i,k,jds)
1072             ENDDO
1073             ENDDO
1075         ENDIF open_ys
1077 !  now the open boundary copy at ye
1079         open_ye: IF( ( config_flags%open_ye   .or. &
1080                        config_flags%polar     .or. &
1081                        config_flags%specified .or. &
1082                        config_flags%nested            ) .and.  &
1083                          ( jte == jde ) .and. open_bc_copy )  THEN
1085           IF (variable /= 'v' .and. variable /= 'y' ) THEN
1087             DO k = kts, k_end
1088             DO i = i_start, i_end
1089               dat(i,k,jde  ) = dat(i,k,jde-1)
1090               dat(i,k,jde+1) = dat(i,k,jde-1)
1091               dat(i,k,jde+2) = dat(i,k,jde-1)
1092             ENDDO                               
1093             ENDDO
1095           ELSE
1097             DO k = kts, k_end
1098             DO i = i_start, i_end
1099               dat(i,k,jde+1) = dat(i,k,jde)
1100               dat(i,k,jde+2) = dat(i,k,jde)
1101               dat(i,k,jde+3) = dat(i,k,jde)
1102             ENDDO                               
1103             ENDDO
1105           ENDIF
1107       END IF open_ye
1109 !  end open b.c in Y copy into boundary zone addition.  WCS, 19 March 2000
1111       END IF periodicity_y
1113    END SUBROUTINE set_physical_bc3d
1115    SUBROUTINE init_module_bc
1116    END SUBROUTINE init_module_bc
1118 !------------------------------------------------------------------------
1120 ! a couple versions of this call to allow a smaller-than-memory dimensioned field (e.g. tile sized)
1121 ! to be passed in as the first argument.  Both of these call the _core version defined below.
1122    SUBROUTINE relax_bdytend   ( field, field_tend,                     &
1123                                 field_bdy_xs, field_bdy_xe,            &
1124                                 field_bdy_ys, field_bdy_ye,            &
1125                                 field_bdy_tend_xs, field_bdy_tend_xe,  &
1126                                 field_bdy_tend_ys, field_bdy_tend_ye,  &
1127                                 variable_in, config_flags,             &
1128                                 spec_bdy_width, spec_zone, relax_zone, &
1129                                 dtbc, fcx, gcx,             &
1130                                 ids,ide, jds,jde, kds,kde,  & ! domain dims
1131                                 ims,ime, jms,jme, kms,kme,  & ! memory dims
1132                                 ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1133                                 its,ite, jts,jte, kts,kte   &
1134                                 )
1136       IMPLICIT NONE
1138       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1139       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1140       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
1141       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1142       INTEGER,      INTENT(IN   )    :: spec_bdy_width, spec_zone, relax_zone
1143       REAL,         INTENT(IN   )    :: dtbc
1144       CHARACTER,    INTENT(IN   )    :: variable_in
1146       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: field
1147       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field_tend
1148       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_xs, field_bdy_xe
1149       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_ys, field_bdy_ye
1150       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_tend_xs, field_bdy_tend_xe
1151       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_tend_ys, field_bdy_tend_ye
1152       REAL,  DIMENSION( spec_bdy_width ), INTENT(IN   ) :: fcx, gcx
1153       TYPE( grid_config_rec_type ) config_flags
1155       CALL relax_bdytend_core   ( field, field_tend,                     &
1156                                 field_bdy_xs, field_bdy_xe,            &
1157                                 field_bdy_ys, field_bdy_ye,            &
1158                                 field_bdy_tend_xs, field_bdy_tend_xe,  &
1159                                 field_bdy_tend_ys, field_bdy_tend_ye,  &
1160                                 variable_in, config_flags,             &
1161                                 spec_bdy_width, spec_zone, relax_zone, &
1162                                 dtbc, fcx, gcx,             &
1163                                 ids,ide, jds,jde, kds,kde,  & ! domain dims
1164                                 ims,ime, jms,jme, kms,kme,  & ! memory dims
1165                                 ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1166                                 its,ite, jts,jte, kts,kte,  & ! patch  dims
1167                                 ims,ime, jms,jme, kms,kme )  ! dimension of the field argument
1168    END SUBROUTINE relax_bdytend
1170 ! version that allows tile-sized version of field. Note, caller should define the
1171 ! field to be -+1 of tile size in each dimension because routine is going off onto halo
1172 ! for example, see relax_bdytend in dyn_em/module_bc_em.F
1173    SUBROUTINE relax_bdytend_tile   ( field, field_tend,                     &
1174                                 field_bdy_xs, field_bdy_xe,            &
1175                                 field_bdy_ys, field_bdy_ye,            &
1176                                 field_bdy_tend_xs, field_bdy_tend_xe,  &
1177                                 field_bdy_tend_ys, field_bdy_tend_ye,  &
1178                                 variable_in, config_flags,             &
1179                                 spec_bdy_width, spec_zone, relax_zone, &
1180                                 dtbc, fcx, gcx,             &
1181                                 ids,ide, jds,jde, kds,kde,  & ! domain dims
1182                                 ims,ime, jms,jme, kms,kme,  & ! memory dims
1183                                 ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1184                                 its,ite, jts,jte, kts,kte,  &
1185                                 iXs,iXe, jXs,jXe, kXs,kXe   &  ! dims of first argument
1186                                 )
1188       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1189       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1190       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
1191       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1192       INTEGER,      INTENT(IN   )    :: iXs,iXe, jXs,jXe, kXs,kXe
1193       INTEGER,      INTENT(IN   )    :: spec_bdy_width, spec_zone, relax_zone
1194       REAL,         INTENT(IN   )    :: dtbc
1195       CHARACTER,    INTENT(IN   )    :: variable_in
1197       REAL,  DIMENSION( iXs:iXe , kXs:kXe , jXs:jXe ), INTENT(IN   ) :: field
1198       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field_tend
1199       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_xs, field_bdy_xe
1200       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_ys, field_bdy_ye
1201       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_tend_xs, field_bdy_tend_xe
1202       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_tend_ys, field_bdy_tend_ye
1203       REAL,  DIMENSION( spec_bdy_width ), INTENT(IN   ) :: fcx, gcx
1204       TYPE( grid_config_rec_type ) config_flags
1206       CALL relax_bdytend_core   ( field, field_tend,                     &
1207                                 field_bdy_xs, field_bdy_xe,            &
1208                                 field_bdy_ys, field_bdy_ye,            &
1209                                 field_bdy_tend_xs, field_bdy_tend_xe,  &
1210                                 field_bdy_tend_ys, field_bdy_tend_ye,  &
1211                                 variable_in, config_flags,             &
1212                                 spec_bdy_width, spec_zone, relax_zone, &
1213                                 dtbc, fcx, gcx,             &
1214                                 ids,ide, jds,jde, kds,kde,  & ! domain dims
1215                                 ims,ime, jms,jme, kms,kme,  & ! memory dims
1216                                 ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1217                                 its,ite, jts,jte, kts,kte,  &
1218                                 iXs,iXe, jXs,jXe, kXs,kXe )  ! dimension of the field argument
1219    END SUBROUTINE relax_bdytend_tile
1221    SUBROUTINE relax_bdytend_core   ( field, field_tend,                     &
1222                                 field_bdy_xs, field_bdy_xe,            &
1223                                 field_bdy_ys, field_bdy_ye,            &
1224                                 field_bdy_tend_xs, field_bdy_tend_xe,  &
1225                                 field_bdy_tend_ys, field_bdy_tend_ye,  &
1226                                 variable_in, config_flags,             &
1227                                 spec_bdy_width, spec_zone, relax_zone, &
1228                                 dtbc, fcx, gcx,             &
1229                                 ids,ide, jds,jde, kds,kde,  & ! domain dims
1230                                 ims,ime, jms,jme, kms,kme,  & ! memory dims
1231                                 ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1232                                 its,ite, jts,jte, kts,kte,  & ! patch  dims
1233                                 iXs,iXe, jXs,jXe, kXs,kXe   & ! field (1st arg) dims; might be tile or patch
1234                                 )
1236 !  This subroutine adds the tendencies in the boundary relaxation region, for specified
1237 !  boundary conditions.  
1238 !  spec_bdy_width is only used to dimension the boundary arrays.
1239 !  relax_zone is the inner edge of the boundary relaxation zone treated here.
1240 !  spec_zone is the width of the outer specified b.c.s that are not changed here.
1241 !  (JD July 2000)
1243       IMPLICIT NONE
1245       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1246       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1247       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
1248       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1249       INTEGER,      INTENT(IN   )    :: iXs,iXe, jXs,jXe, kXs,kXe
1250       INTEGER,      INTENT(IN   )    :: spec_bdy_width, spec_zone, relax_zone
1251       REAL,         INTENT(IN   )    :: dtbc
1252       CHARACTER,    INTENT(IN   )    :: variable_in
1255       REAL,  DIMENSION( iXs:iXe , kXs:kXe , jXs:jXe ), INTENT(IN   ) :: field
1256       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field_tend
1257       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_xs, field_bdy_xe
1258       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_ys, field_bdy_ye
1259       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_tend_xs, field_bdy_tend_xe
1260       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_tend_ys, field_bdy_tend_ye
1261       REAL,  DIMENSION( spec_bdy_width ), INTENT(IN   ) :: fcx, gcx
1262       TYPE( grid_config_rec_type ) config_flags
1264       CHARACTER  :: variable
1265       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf, im1, ip1
1266       INTEGER    :: b_dist, b_limit
1267       REAL       :: fls0, fls1, fls2, fls3, fls4
1268       LOGICAL    :: periodic_x
1270       periodic_x = config_flags%periodic_x
1271       variable = variable_in
1273       IF (variable == 'U') variable = 'u'
1274       IF (variable == 'V') variable = 'v'
1275       IF (variable == 'M') variable = 'm'
1276       IF (variable == 'H') variable = 'h'
1278       ibs = ids
1279       ibe = ide-1
1280       itf = min(ite,ide-1)
1281       jbs = jds
1282       jbe = jde-1
1283       jtf = min(jte,jde-1)
1284       ktf = kde-1
1285       IF (variable == 'u') ibe = ide
1286       IF (variable == 'u') itf = min(ite,ide)
1287       IF (variable == 'v') jbe = jde
1288       IF (variable == 'v') jtf = min(jte,jde)
1289       IF (variable == 'm') ktf = kte
1290       IF (variable == 'h') ktf = kte
1293       IF (jts - jbs .lt. relax_zone) THEN
1294 ! Y-start boundary
1295         DO j = max(jts,jbs+spec_zone), min(jtf,jbs+relax_zone-1)
1296           b_dist = j - jbs
1297           b_limit = b_dist
1298           IF(periodic_x)b_limit = 0
1299           DO k = kts, ktf
1300             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1301               im1 = max(i-1,ibs)
1302               ip1 = min(i+1,ibe)
1303               fls0 = field_bdy_ys(i, k, b_dist+1) &
1304                    + dtbc * field_bdy_tend_ys(i, k, b_dist+1) &
1305                    - field(i,k,j)
1306               fls1 = field_bdy_ys(im1, k, b_dist+1) &
1307                    + dtbc * field_bdy_tend_ys(im1, k, b_dist+1) &
1308                    - field(im1,k,j)
1309               fls2 = field_bdy_ys(ip1, k, b_dist+1) &
1310                    + dtbc * field_bdy_tend_ys(ip1, k, b_dist+1) &
1311                    - field(ip1,k,j)
1312               fls3 = field_bdy_ys(i, k, b_dist) &
1313                    + dtbc * field_bdy_tend_ys(i, k, b_dist) &
1314                    - field(i,k,j-1)
1315               fls4 = field_bdy_ys(i, k, b_dist+2) &
1316                    + dtbc * field_bdy_tend_ys(i, k, b_dist+2) &
1317                    - field(i,k,j+1)
1318               field_tend(i,k,j) = field_tend(i,k,j) &
1319                                 + fcx(b_dist+1)*fls0 &
1320                                 - gcx(b_dist+1)*(fls1+fls2+fls3+fls4-4.*fls0)
1321             ENDDO
1322           ENDDO
1323         ENDDO
1324       ENDIF
1326       IF (jbe - jtf .lt. relax_zone) THEN
1329 ! Y-end boundary
1330         DO j = max(jts,jbe-relax_zone+1), min(jtf,jbe-spec_zone)
1331           b_dist = jbe - j
1332           b_limit = b_dist
1333           IF(periodic_x)b_limit = 0
1336           DO k = kts, ktf
1337             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1338               im1 = max(i-1,ibs)
1339               ip1 = min(i+1,ibe)
1340               fls0 = field_bdy_ye(i, k, b_dist+1) &
1341                    + dtbc * field_bdy_tend_ye(i, k, b_dist+1) &
1342                    - field(i,k,j)
1343               fls1 = field_bdy_ye(im1, k, b_dist+1) &
1344                    + dtbc * field_bdy_tend_ye(im1, k, b_dist+1) &
1345                    - field(im1,k,j)
1346               fls2 = field_bdy_ye(ip1, k, b_dist+1) &
1347                    + dtbc * field_bdy_tend_ye(ip1, k, b_dist+1) &
1348                    - field(ip1,k,j)
1349               fls3 = field_bdy_ye(i, k, b_dist) &
1350                    + dtbc * field_bdy_tend_ye(i, k, b_dist) &
1351                    - field(i,k,j+1)
1352               fls4 = field_bdy_ye(i, k, b_dist+2) &
1353                    + dtbc * field_bdy_tend_ye(i, k, b_dist+2) &
1354                    - field(i,k,j-1)
1355               field_tend(i,k,j) = field_tend(i,k,j) &
1356                                 + fcx(b_dist+1)*fls0 &
1357                                 - gcx(b_dist+1)*(fls1+fls2+fls3+fls4-4.*fls0)
1359             ENDDO
1360           ENDDO
1361         ENDDO
1362       ENDIF
1364     IF(.NOT.periodic_x)THEN
1365       IF (its - ibs .lt. relax_zone) THEN
1366 ! X-start boundary
1367         DO i = max(its,ibs+spec_zone), min(itf,ibs+relax_zone-1)
1368           b_dist = i - ibs
1369           DO k = kts, ktf
1370             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1371               fls0 = field_bdy_xs(j, k, b_dist+1) &
1372                    + dtbc * field_bdy_tend_xs(j, k, b_dist+1) &
1373                    - field(i,k,j)
1374               fls1 = field_bdy_xs(j-1, k, b_dist+1) &
1375                    + dtbc * field_bdy_tend_xs(j-1, k, b_dist+1) &
1376                    - field(i,k,j-1)
1377               fls2 = field_bdy_xs(j+1, k, b_dist+1) &
1378                    + dtbc * field_bdy_tend_xs(j+1, k, b_dist+1) &
1379                    - field(i,k,j+1)
1380               fls3 = field_bdy_xs(j, k, b_dist) &
1381                    + dtbc * field_bdy_tend_xs(j, k, b_dist) &
1382                    - field(i-1,k,j)
1383               fls4 = field_bdy_xs(j, k, b_dist+2) &
1384                    + dtbc * field_bdy_tend_xs(j, k, b_dist+2) &
1385                    - field(i+1,k,j)
1386               field_tend(i,k,j) = field_tend(i,k,j) &
1387                                 + fcx(b_dist+1)*fls0 &
1388                                 - gcx(b_dist+1)*(fls1+fls2+fls3+fls4-4.*fls0)
1390             ENDDO
1391           ENDDO
1392         ENDDO
1393       ENDIF
1395       IF (ibe - itf .lt. relax_zone) THEN
1396 ! X-end boundary
1397         DO i = max(its,ibe-relax_zone+1), min(itf,ibe-spec_zone)
1398           b_dist = ibe - i
1399           DO k = kts, ktf
1400             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1401               fls0 = field_bdy_xe(j, k, b_dist+1) &
1402                    + dtbc * field_bdy_tend_xe(j, k, b_dist+1) &
1403                    - field(i,k,j)
1404               fls1 = field_bdy_xe(j-1, k, b_dist+1) &
1405                    + dtbc * field_bdy_tend_xe(j-1, k, b_dist+1) &
1406                    - field(i,k,j-1)
1407               fls2 = field_bdy_xe(j+1, k, b_dist+1) &
1408                    + dtbc * field_bdy_tend_xe(j+1, k, b_dist+1) &
1409                    - field(i,k,j+1)
1410               fls3 = field_bdy_xe(j, k, b_dist) &
1411                    + dtbc * field_bdy_tend_xe(j, k, b_dist) &
1412                    - field(i+1,k,j)
1413               fls4 = field_bdy_xe(j, k, b_dist+2) &
1414                    + dtbc * field_bdy_tend_xe(j, k, b_dist+2) &
1415                    - field(i-1,k,j)
1416               field_tend(i,k,j) = field_tend(i,k,j) &
1417                                 + fcx(b_dist+1)*fls0 &
1418                                 - gcx(b_dist+1)*(fls1+fls2+fls3+fls4-4.*fls0)
1419             ENDDO
1420           ENDDO
1421         ENDDO
1422       ENDIF
1423     ENDIF
1427    END SUBROUTINE relax_bdytend_core
1428 !------------------------------------------------------------------------
1430    SUBROUTINE spec_bdytend   ( field_tend,                           &
1431                                field_bdy_xs, field_bdy_xe,           &
1432                                field_bdy_ys, field_bdy_ye,           &
1433                                field_bdy_tend_xs, field_bdy_tend_xe, &
1434                                field_bdy_tend_ys, field_bdy_tend_ye, &
1435                                variable_in, config_flags, &
1436                                spec_bdy_width, spec_zone, &
1437                                ids,ide, jds,jde, kds,kde,  & ! domain dims
1438                                ims,ime, jms,jme, kms,kme,  & ! memory dims
1439                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1440                                its,ite, jts,jte, kts,kte )
1442 !  This subroutine sets the tendencies in the boundary specified region.
1443 !  spec_bdy_width is only used to dimension the boundary arrays.
1444 !  spec_zone is the width of the outer specified b.c.s that are set here.
1445 !  (JD July 2000)
1447       IMPLICIT NONE
1449       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1450       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1451       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
1452       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1453       INTEGER,      INTENT(IN   )    :: spec_bdy_width, spec_zone
1454       CHARACTER,    INTENT(IN   )    :: variable_in
1457       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(OUT  ) :: field_tend
1458       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_xs, field_bdy_xe
1459       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_ys, field_bdy_ye
1460       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_tend_xs, field_bdy_tend_xe
1461       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_tend_ys, field_bdy_tend_ye
1462       TYPE( grid_config_rec_type ) config_flags
1464       CHARACTER  :: variable
1465       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf
1466       INTEGER    :: b_dist, b_limit
1467       LOGICAL    :: periodic_x
1469       periodic_x = config_flags%periodic_x
1471       variable = variable_in
1473       IF (variable == 'U') variable = 'u'
1474       IF (variable == 'V') variable = 'v'
1475       IF (variable == 'M') variable = 'm'
1476       IF (variable == 'H') variable = 'h'
1478       ibs = ids
1479       ibe = ide-1
1480       itf = min(ite,ide-1)
1481       jbs = jds
1482       jbe = jde-1
1483       jtf = min(jte,jde-1)
1484       ktf = kde-1
1485       IF (variable == 'u') ibe = ide
1486       IF (variable == 'u') itf = min(ite,ide)
1487       IF (variable == 'v') jbe = jde
1488       IF (variable == 'v') jtf = min(jte,jde)
1489       IF (variable == 'm') ktf = kte
1490       IF (variable == 'h') ktf = kte
1492       IF (jts - jbs .lt. spec_zone) THEN
1493 ! Y-start boundary
1494         DO j = jts, min(jtf,jbs+spec_zone-1)
1495           b_dist = j - jbs
1496           b_limit = b_dist
1497           IF(periodic_x)b_limit = 0
1498           DO k = kts, ktf
1499             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1500               field_tend(i,k,j) = field_bdy_tend_ys(i, k, b_dist+1)
1501             ENDDO
1502           ENDDO
1503         ENDDO
1504       ENDIF
1505       IF (jbe - jtf .lt. spec_zone) THEN
1508 ! Y-end boundary
1509         DO j = max(jts,jbe-spec_zone+1), jtf
1510           b_dist = jbe - j
1511           b_limit = b_dist
1512           IF(periodic_x)b_limit = 0
1515           DO k = kts, ktf
1516             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1517               field_tend(i,k,j) = field_bdy_tend_ye(i, k, b_dist+1)
1518             ENDDO
1519           ENDDO
1520         ENDDO
1521       ENDIF
1523     IF(.NOT.periodic_x)THEN
1524       IF (its - ibs .lt. spec_zone) THEN
1525 ! X-start boundary
1526         DO i = its, min(itf,ibs+spec_zone-1)
1527           b_dist = i - ibs
1528           DO k = kts, ktf
1529             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1530               field_tend(i,k,j) = field_bdy_tend_xs(j, k, b_dist+1)
1531             ENDDO
1532           ENDDO
1533         ENDDO
1534       ENDIF
1536       IF (ibe - itf .lt. spec_zone) THEN
1537 ! X-end boundary
1538         DO i = max(its,ibe-spec_zone+1), itf
1539           b_dist = ibe - i
1540           DO k = kts, ktf
1541             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1542               field_tend(i,k,j) = field_bdy_tend_xe(j, k, b_dist+1)
1543             ENDDO
1544           ENDDO
1545         ENDDO
1546       ENDIF
1547     ENDIF
1549    END SUBROUTINE spec_bdytend
1550    
1551 !------------------------------------------------------------------------
1552 ! KRS 9/13/2012: New subroutine spec_bdytend_perturb.
1553 ! Reads in field_tend_perturb from dyn_em/module_bc_em.F  
1554 ! field_tend_perturb=r*_tendf_stoch if perturb_bdy=1
1555 ! field_tend_perturb=User provided pattern if perturb_bdy=2
1556 ! User can fill array in this subroutine.   
1557 ! Adds perturbations to boundaries. Outputs field_tend (e.g. ru_tend).
1558 ! mu_2, mub, and msf* passed into subrountine couple variables.
1559 !------------------------------------------------------------------------
1560    SUBROUTINE spec_bdytend_perturb   ( field_tend,                   &
1561                                        field_tend_perturb,           &
1562                                        mu_2, mub, c1, c2,            &
1563                                        variable_in,                  &
1564                                        msf, config_flags,            &
1565                                        spec_bdy_width, spec_zone,    &
1566                                        kme_stoch,                    & ! stoch  dims
1567                                        ids,ide, jds,jde, kds,kde,    & ! domain dims
1568                                        ims,ime, jms,jme, kms,kme,    & ! memory dims
1569                                        ips,ipe, jps,jpe, kps,kpe,    & ! patch  dims
1570                                        its,ite, jts,jte, kts,kte )
1572 !  This subroutine sets the tendencies in the boundary specified region.
1573 !  spec_bdy_width is only used to dimension the boundary arrays.
1574 !  spec_zone is the width of the outer specified b.c.s that are set here.
1575 !  (JD July 2000)
1577       IMPLICIT NONE
1579       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1580       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme, kme_stoch
1581       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
1582       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1583       INTEGER,      INTENT(IN   )    :: spec_bdy_width, spec_zone
1584       CHARACTER,    INTENT(IN   )    :: variable_in
1587       REAL,  DIMENSION( ims:ime , kms:kme       , jms:jme ), INTENT(INOUT  ) :: field_tend
1588       REAL,  DIMENSION( ims:ime , kms:kme_stoch , jms:jme ), INTENT(IN   )   :: field_tend_perturb
1589       REAL,  DIMENSION( ims:ime ,                 jms:jme ), INTENT(IN   )   :: mu_2
1590       REAL,  DIMENSION( ims:ime ,                 jms:jme ), INTENT(IN   )   :: mub
1591       REAL,  DIMENSION( ims:ime ,                 jms:jme ), INTENT(IN   )   :: msf
1592       REAL,  DIMENSION(           kms:kme                 ), INTENT(IN   )   :: c1, c2
1593    
1594       TYPE( grid_config_rec_type ) config_flags
1596       CHARACTER  :: variable
1597       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf
1598       INTEGER    :: b_dist, b_limit
1599       LOGICAL    :: periodic_x
1601       periodic_x = config_flags%periodic_x
1603       variable = variable_in
1605       IF (variable == 'U') variable = 'u'
1606       IF (variable == 'V') variable = 'v'
1607       IF (variable == 'T') variable = 't'
1608       IF (variable == 'H') variable = 'h'
1610       ibs = ids
1611       ibe = ide-1
1612       itf = min(ite,ide-1)
1613       jbs = jds
1614       jbe = jde-1
1615       jtf = min(jte,jde-1)
1616       ktf = kde-1
1617       IF (variable == 'u') ibe = ide
1618       IF (variable == 'u') itf = min(ite,ide)
1619       IF (variable == 'v') jbe = jde
1620       IF (variable == 'v') jtf = min(jte,jde)
1621       IF (variable == 'h') ktf = kte
1623       IF (jts - jbs .lt. spec_zone) THEN
1625 ! Y-start boundary
1626         DO j = jts, min(jtf,jbs+spec_zone-1)
1627           b_dist = j - jbs
1628           b_limit = b_dist
1629           IF(periodic_x)b_limit = 0
1630           DO k = kts, ktf
1631             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1632             IF (variable == 't') THEN
1633               field_tend(i,k,j) = field_tend(i,k,j) + &
1634                                   field_tend_perturb(i,min(k,kme_stoch),j) * &
1635                                   ((c1(k)*mu_2(i,j))+(c1(k)*mub(i,j)+c2(k)))
1636             ENDIF
1637             IF (variable == 'u') THEN
1638               field_tend(i,k,j) = field_tend(i,k,j) + &
1639                                   field_tend_perturb(i,min(k,kme_stoch),j) * &
1640                                   0.5*((c1(k)*mu_2(i,j))+(c1(k)*mu_2(i-1,j)) + (c1(k)*mub(i,j)+c2(k))+(c1(k)*mub(i-1,j)+c2(k))) / msf(i,j)
1641             ENDIF
1642             IF (variable == 'v') THEN
1643               field_tend(i,k,j) = field_tend(i,k,j) + &
1644                                   field_tend_perturb(i,min(k,kme_stoch),j) * &
1645                                   0.5*((c1(k)*mu_2(i,j))+(c1(k)*mu_2(i,j-1)) + (c1(k)*mub(i,j)+c2(k))+(c1(k)*mub(i,j-1)+c2(k))) / msf(i,j)
1646             ENDIF
1647             ENDDO
1648           ENDDO
1649         ENDDO
1650       ENDIF
1651       IF (jbe - jtf .lt. spec_zone) THEN
1654 ! Y-end boundary
1655         DO j = max(jts,jbe-spec_zone+1), jtf
1656           b_dist = jbe - j
1657           b_limit = b_dist
1658           IF(periodic_x)b_limit = 0
1659           DO k = kts, ktf
1660             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1661             IF (variable == 't') THEN
1662               field_tend(i,k,j) = field_tend(i,k,j) + &
1663                                   field_tend_perturb(i,min(k,kme_stoch),j) * ((c1(k)*mu_2(i,j))+(c1(k)*mub(i,j)+c2(k)))
1664             ENDIF
1665             IF (variable == 'u') THEN
1666               field_tend(i,k,j) = field_tend(i,k,j) + &
1667                                   field_tend_perturb(i,min(k,kme_stoch),j) * 0.5*((c1(k)*mu_2(i,j))+(c1(k)*mu_2(i-1,j)) + (c1(k)*mub(i,j)+c2(k))+(c1(k)*mub(i-1,j)+c2(k))) / msf(i,j)
1668             ENDIF
1669             IF (variable == 'v') THEN
1670               field_tend(i,k,j) = field_tend(i,k,j) + &
1671                                   field_tend_perturb(i,min(k,kme_stoch),j) * 0.5*((c1(k)*mu_2(i,j))+(c1(k)*mu_2(i,j-1)) + (c1(k)*mub(i,j)+c2(k))+(c1(k)*mub(i,j-1)+c2(k))) / msf(i,j)
1672             ENDIF
1673      ENDDO
1674           ENDDO
1675         ENDDO
1676       ENDIF
1678     IF(.NOT.periodic_x)THEN
1679       IF (its - ibs .lt. spec_zone) THEN
1680 ! X-start boundary
1681         DO i = its, min(itf,ibs+spec_zone-1)
1682           b_dist = i - ibs
1683           DO k = kts, ktf
1684             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1685             IF (variable == 't') THEN
1686               field_tend(i,k,j) = field_tend(i,k,j) + &
1687                                   field_tend_perturb(i,min(k,kme_stoch),j) * ((c1(k)*mu_2(i,j))+(c1(k)*mub(i,j)+c2(k)))
1688             ENDIF
1689             IF (variable == 'u') THEN
1690               field_tend(i,k,j) = field_tend(i,k,j) + &
1691                                   field_tend_perturb(i,min(k,kme_stoch),j) * 0.5*((c1(k)*mu_2(i,j))+(c1(k)*mu_2(i-1,j)) + (c1(k)*mub(i,j)+c2(k))+(c1(k)*mub(i-1,j)+c2(k))) / msf(i,j)
1692             ENDIF
1693             IF (variable == 'v') THEN
1694               field_tend(i,k,j) = field_tend(i,k,j) + &
1695                                   field_tend_perturb(i,min(k,kme_stoch),j) * 0.5*((c1(k)*mu_2(i,j))+(c1(k)*mu_2(i,j-1)) + (c1(k)*mub(i,j)+c2(k))+(c1(k)*mub(i,j-1)+c2(k))) / msf(i,j)
1696             ENDIF
1697             ENDDO
1698           ENDDO
1699         ENDDO
1700       ENDIF
1702       IF (ibe - itf .lt. spec_zone) THEN
1703 ! X-end boundary
1704         DO i = max(its,ibe-spec_zone+1), itf
1705           b_dist = ibe - i
1706           DO k = kts, ktf
1707             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1708             IF (variable == 't') THEN
1709               field_tend(i,k,j) = field_tend(i,k,j) + &
1710                                   field_tend_perturb(i,min(k,kme_stoch),j) * ((c1(k)*mu_2(i,j))+(c1(k)*mub(i,j)+c2(k)))
1711             ENDIF
1712             IF (variable == 'u') THEN
1713               field_tend(i,k,j) = field_tend(i,k,j) + &
1714                                   field_tend_perturb(i,min(k,kme_stoch),j) * 0.5*((c1(k)*mu_2(i,j))+(c1(k)*mu_2(i-1,j)) + (c1(k)*mub(i,j)+c2(k))+(c1(k)*mub(i-1,j)+c2(k))) / msf(i,j)
1715             ENDIF
1716             IF (variable == 'v') THEN
1717               field_tend(i,k,j) = field_tend(i,k,j) + &
1718                                   field_tend_perturb(i,min(k,kme_stoch),j) * 0.5*((c1(k)*mu_2(i,j))+(c1(k)*mu_2(i,j-1)) + (c1(k)*mub(i,j)+c2(k))+(c1(k)*mub(i,j-1)+c2(k))) / msf(i,j)
1719             ENDIF
1720             ENDDO
1721           ENDDO
1722         ENDDO
1723       ENDIF
1724     ENDIF
1726    END SUBROUTINE spec_bdytend_perturb
1729 !------------------------------------------------------------------------
1731    SUBROUTINE spec_bdytend_perturb_chem   ( field_bdy_tend_xs, field_bdy_tend_xe, &
1732                                             field_bdy_tend_ys, field_bdy_tend_ye, &
1733                                             field_scalar_perturb,         &
1734                                             variable_in,                &
1735                                             periodic_x,                 &
1736                                             spec_bdy_width, spec_zone,  &
1737                                             kme_stoch,                  &
1738                                             ids,ide, jds,jde, kds,kde,  &
1739                                             ims,ime, jms,jme, kms,kme,  &
1740                                             ips,ipe, jps,jpe, kps,kpe,  &
1741                                             its,ite, jts,jte, kts,kte )
1743       IMPLICIT NONE
1745       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1746       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme, kme_stoch
1747       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
1748       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1749       INTEGER,      INTENT(IN   )    :: spec_bdy_width, spec_zone
1750       CHARACTER,    INTENT(IN   )    :: variable_in
1751       LOGICAL,      INTENT(IN   )    :: periodic_x
1753       REAL,  DIMENSION( ims:ime , kms:kme_stoch , jms:jme ) , INTENT(IN   ) :: field_scalar_perturb
1754       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(INOUT) :: field_bdy_tend_xs, field_bdy_tend_xe
1755       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(INOUT) :: field_bdy_tend_ys, field_bdy_tend_ye
1757       CHARACTER  :: variable
1758       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf
1759       INTEGER    :: b_dist, b_limit
1761       variable = variable_in
1762       IF (variable == 'C') variable = 'c'
1764       IF (variable /= 'c') THEN
1765           write( wrf_err_message ,*) ' *** Error in spec_bdytend_perturb_chem'
1766           CALL wrf_message ( wrf_err_message )
1767       ENDIF
1769       ibs = ids
1770       ibe = ide-1
1771       itf = min(ite,ide-1)
1772       jbs = jds
1773       jbe = jde-1
1774       jtf = min(jte,jde-1)
1775       ktf = kde-1
1776       !ktf = kte
1778       IF (jts - jbs .lt. spec_zone) THEN
1779 ! Y-start boundary
1780         DO j = jts, min(jtf,jbs+spec_zone-1)
1781            b_dist = j - jbs
1782            b_limit = b_dist
1783            IF(periodic_x)b_limit = 0
1784            DO k = kts, ktf
1785               DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1786                  field_bdy_tend_ys(i,k,b_dist+1) = field_bdy_tend_ys(i,k,b_dist+1) * &
1787                                                    (1.0 + field_scalar_perturb(i,min(k,kme_stoch),j))
1788               ENDDO
1789            ENDDO
1790         ENDDO
1791       ENDIF
1793       IF (jbe - jtf .lt. spec_zone) THEN
1794 ! Y-end boundary
1795         DO j = max(jts,jbe-spec_zone+1), jtf
1796            b_dist = jbe - j
1797            b_limit = b_dist
1798            IF(periodic_x)b_limit = 0
1799            DO k = kts, ktf
1800               DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1801                  field_bdy_tend_ye(i,k,b_dist+1) = field_bdy_tend_ye(i,k,b_dist+1) * &
1802                                                    (1.0 + field_scalar_perturb(i,min(k,kme_stoch),j))
1803               ENDDO
1804            ENDDO
1805         ENDDO
1806       ENDIF
1808     IF(.NOT.periodic_x)THEN
1810       IF (its - ibs .lt. spec_zone) THEN
1811 ! X-start boundary
1812         DO i = its, min(itf,ibs+spec_zone-1)
1813            b_dist = i - ibs
1814            DO k = kts, ktf
1815               DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1816                  field_bdy_tend_xs(j,k,b_dist+1) = field_bdy_tend_xs(j,k,b_dist+1) * &
1817                                                    (1.0 + field_scalar_perturb(i,min(k,kme_stoch),j))
1818               ENDDO
1819            ENDDO
1820         ENDDO
1821       ENDIF
1823       IF (ibe - itf .lt. spec_zone) THEN
1824 ! X-end boundary
1825         DO i = max(its,ibe-spec_zone+1), itf
1826            b_dist = ibe - i
1827            DO k = kts, ktf
1828              DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1829                  field_bdy_tend_xe(j,k,b_dist+1) = field_bdy_tend_xe(j,k,b_dist+1) * (1.0 + field_scalar_perturb(i,k,j))
1830               ENDDO
1831            ENDDO
1832         ENDDO
1833       ENDIF
1835     ENDIF  !  IF(.NOT.periodic_x)THEN
1837    END SUBROUTINE spec_bdytend_perturb_chem
1839 !------------------------------------------------------------------------
1841    SUBROUTINE spec_bdyfield   ( field,                     &
1842                                field_bdy_xs, field_bdy_xe,           &
1843                                field_bdy_ys, field_bdy_ye,           &
1844                                variable_in, config_flags,  &
1845                                spec_bdy_width, spec_zone, &
1846                                ids,ide, jds,jde, kds,kde,  & ! domain dims
1847                                ims,ime, jms,jme, kms,kme,  & ! memory dims
1848                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1849                                its,ite, jts,jte, kts,kte )
1851 !  This subroutine sets the tendencies in the boundary specified region.
1852 !  spec_bdy_width is only used to dimension the boundary arrays.
1853 !  spec_zone is the width of the outer specified b.c.s that are set here.
1854 !  (JD July 2000)
1856       IMPLICIT NONE
1858       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1859       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1860       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
1861       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1862       INTEGER,      INTENT(IN   )    :: spec_bdy_width, spec_zone
1863       CHARACTER,    INTENT(IN   )    :: variable_in
1866       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(OUT  ) :: field
1867       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_xs, field_bdy_xe
1868       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_ys, field_bdy_ye
1869       TYPE( grid_config_rec_type ) config_flags
1871       CHARACTER  :: variable
1872       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf
1873       INTEGER    :: b_dist, b_limit
1874       LOGICAL    :: periodic_x
1876       periodic_x = config_flags%periodic_x
1878       variable = variable_in
1880       IF (variable == 'U') variable = 'u'
1881       IF (variable == 'V') variable = 'v'
1882       IF (variable == 'M') variable = 'm'
1883       IF (variable == 'H') variable = 'h'
1885       ibs = ids
1886       ibe = ide-1
1887       itf = min(ite,ide-1)
1888       jbs = jds
1889       jbe = jde-1
1890       jtf = min(jte,jde-1)
1891       ktf = kde-1
1892       IF (variable == 'u') ibe = ide
1893       IF (variable == 'u') itf = min(ite,ide)
1894       IF (variable == 'v') jbe = jde
1895       IF (variable == 'v') jtf = min(jte,jde)
1896       IF (variable == 'm') ktf = kte
1897       IF (variable == 'h') ktf = kte
1899       IF (jts - jbs .lt. spec_zone) THEN
1900 ! Y-start boundary
1901         DO j = jts, min(jtf,jbs+spec_zone-1)
1902           b_dist = j - jbs
1903           b_limit = b_dist
1904           IF(periodic_x)b_limit = 0
1905           DO k = kts, ktf
1906             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1907               field(i,k,j) = field_bdy_ys(i, k, b_dist+1)
1908             ENDDO
1909           ENDDO
1910         ENDDO
1911       ENDIF
1912       IF (jbe - jtf .lt. spec_zone) THEN
1913 ! Y-end boundary
1914         DO j = max(jts,jbe-spec_zone+1), jtf
1915           b_dist = jbe - j
1916           b_limit = b_dist
1917           IF(periodic_x)b_limit = 0
1918           DO k = kts, ktf
1919             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1920               field(i,k,j) = field_bdy_ye(i, k, b_dist+1)
1921             ENDDO
1922           ENDDO
1923         ENDDO
1924       ENDIF
1926     IF(.NOT.periodic_x)THEN
1927       IF (its - ibs .lt. spec_zone) THEN
1928 ! X-start boundary
1929         DO i = its, min(itf,ibs+spec_zone-1)
1930           b_dist = i - ibs
1931           DO k = kts, ktf
1932             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1933               field(i,k,j) = field_bdy_xs(j, k, b_dist+1)
1934             ENDDO
1935           ENDDO
1936         ENDDO
1937       ENDIF
1939       IF (ibe - itf .lt. spec_zone) THEN
1940 ! X-end boundary
1941         DO i = max(its,ibe-spec_zone+1), itf
1942           b_dist = ibe - i
1943           DO k = kts, ktf
1944             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1945               field(i,k,j) = field_bdy_xe(j, k, b_dist+1)
1946             ENDDO
1947           ENDDO
1948         ENDDO
1949       ENDIF
1950     ENDIF
1952    END SUBROUTINE spec_bdyfield
1953 !------------------------------------------------------------------------
1955    SUBROUTINE spec_bdyupdate(  field,      &
1956                                field_tend, dt,            &
1957                                variable_in, config_flags, &
1958                                spec_zone,                  &
1959                                ids,ide, jds,jde, kds,kde,  & ! domain dims
1960                                ims,ime, jms,jme, kms,kme,  & ! memory dims
1961                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1962                                its,ite, jts,jte, kts,kte )
1964 !  This subroutine adds the tendencies in the boundary specified region.
1965 !  spec_zone is the width of the outer specified b.c.s that are set here.
1966 !  (JD August 2000)
1968       IMPLICIT NONE
1970       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1971       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1972       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
1973       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1974       INTEGER,      INTENT(IN   )    :: spec_zone
1975       CHARACTER,    INTENT(IN   )    :: variable_in
1976       REAL,         INTENT(IN   )    :: dt
1979       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field
1980       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: field_tend
1981       TYPE( grid_config_rec_type ) config_flags
1983       CHARACTER  :: variable
1984       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf
1985       INTEGER    :: b_dist, b_limit
1986       LOGICAL    :: periodic_x
1988       periodic_x = config_flags%periodic_x
1990       variable = variable_in
1992       IF (variable == 'U') variable = 'u'
1993       IF (variable == 'V') variable = 'v'
1994       IF (variable == 'M') variable = 'm'
1995       IF (variable == 'H') variable = 'h'
1997       ibs = ids
1998       ibe = ide-1
1999       itf = min(ite,ide-1)
2000       jbs = jds
2001       jbe = jde-1
2002       jtf = min(jte,jde-1)
2003       ktf = kde-1
2004       IF (variable == 'u') ibe = ide
2005       IF (variable == 'u') itf = min(ite,ide)
2006       IF (variable == 'v') jbe = jde
2007       IF (variable == 'v') jtf = min(jte,jde)
2008       IF (variable == 'm') ktf = kte
2009       IF (variable == 'h') ktf = kte
2011       IF (jts - jbs .lt. spec_zone) THEN
2012 ! Y-start boundary
2013         DO j = jts, min(jtf,jbs+spec_zone-1)
2014           b_dist = j - jbs
2015           b_limit = b_dist
2016           IF(periodic_x)b_limit = 0
2017           DO k = kts, ktf
2018             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2019               field(i,k,j) = field(i,k,j) + dt*field_tend(i,k,j)
2020             ENDDO
2021           ENDDO
2022         ENDDO
2023       ENDIF
2024       IF (jbe - jtf .lt. spec_zone) THEN
2025 ! Y-end boundary
2026         DO j = max(jts,jbe-spec_zone+1), jtf
2027           b_dist = jbe - j
2028           b_limit = b_dist
2029           IF(periodic_x)b_limit = 0
2030           DO k = kts, ktf
2031             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2032               field(i,k,j) = field(i,k,j) + dt*field_tend(i,k,j)
2033             ENDDO
2034           ENDDO
2035         ENDDO
2036       ENDIF
2038     IF(.NOT.periodic_x)THEN
2039       IF (its - ibs .lt. spec_zone) THEN
2040 ! X-start boundary
2041         DO i = its, min(itf,ibs+spec_zone-1)
2042           b_dist = i - ibs
2043           DO k = kts, ktf
2044             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2045               field(i,k,j) = field(i,k,j) + dt*field_tend(i,k,j)
2046             ENDDO
2047           ENDDO
2048         ENDDO
2049       ENDIF
2051       IF (ibe - itf .lt. spec_zone) THEN
2052 ! X-end boundary
2053         DO i = max(its,ibe-spec_zone+1), itf
2054           b_dist = ibe - i
2055           DO k = kts, ktf
2056             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2057               field(i,k,j) = field(i,k,j) + dt*field_tend(i,k,j)
2058             ENDDO
2059           ENDDO
2060         ENDDO
2061       ENDIF
2062     ENDIF
2064    END SUBROUTINE spec_bdyupdate
2065 !------------------------------------------------------------------------
2066    SUBROUTINE spec_bdy_final   ( field, mu, c1, c2, msf,               &
2067                                 field_bdy_xs, field_bdy_xe,            &
2068                                 field_bdy_ys, field_bdy_ye,            &
2069                                 field_bdy_tend_xs, field_bdy_tend_xe,  &
2070                                 field_bdy_tend_ys, field_bdy_tend_ye,  &
2071                                 variable_in, config_flags,             &
2072                                 spec_bdy_width, spec_zone,             &
2073                                 dtbc,                       &
2074                                 ids,ide, jds,jde, kds,kde,  & ! domain dims
2075                                 ims,ime, jms,jme, kms,kme,  & ! memory dims
2076                                 ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
2077                                 its,ite, jts,jte, kts,kte   )
2079 !  This subroutine forces the boundary to match the boundary file value for specified
2080 !  boundary conditions. Added to avoid drift due to round-off error using just tendencies.
2081 !  Boundary-file coupling is u,v,w:mu/msf  other fields:mu
2082 !  Correctly  staggered mu and msf are passed in (as seen in small_step_finish)
2083 !  spec_bdy_width is only used to dimension the boundary arrays.
2084 !  relax_zone is the inner edge of the boundary relaxation zone treated here.
2085 !  spec_zone is the width of the outer specified b.c.s that are not changed here.
2086 !  (JD Jan 2015)
2088       IMPLICIT NONE
2090       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
2091       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
2092       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
2093       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
2094       INTEGER,      INTENT(IN   )    :: spec_bdy_width, spec_zone
2095       REAL,         INTENT(IN   )    :: dtbc
2096       CHARACTER,    INTENT(IN   )    :: variable_in
2099       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field
2100       REAL,  DIMENSION( ims:ime , jms:jme), INTENT(IN   ) :: mu, msf
2101       REAL,  DIMENSION( kms:kme ), INTENT(IN   ) :: c1, c2
2102       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_xs, field_bdy_xe
2103       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_ys, field_bdy_ye
2104       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_tend_xs, field_bdy_tend_xe
2105       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width ), INTENT(IN   ) :: field_bdy_tend_ys, field_bdy_tend_ye
2106       TYPE( grid_config_rec_type ) config_flags
2108       CHARACTER  :: variable
2109       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf, im1, ip1
2110       INTEGER    :: b_dist, b_limit
2111       REAL       :: bfield, xmsf, xmu
2112       LOGICAL    :: periodic_x, msfcouple, mucouple
2114       periodic_x = config_flags%periodic_x
2115       variable = variable_in
2117       IF (variable == 'U') variable = 'u'
2118       IF (variable == 'V') variable = 'v'
2119       IF (variable == 'W') variable = 'w'
2120       IF (variable == 'M') variable = 'm'
2121       IF (variable == 'T') variable = 't'
2122       IF (variable == 'H') variable = 'h'
2124       ibs = ids
2125       ibe = ide-1
2126       itf = min(ite,ide-1)
2127       jbs = jds
2128       jbe = jde-1
2129       jtf = min(jte,jde-1)
2130       ktf = kde-1
2131       IF (variable == 'u') ibe = ide
2132       IF (variable == 'u') itf = min(ite,ide)
2133       IF (variable == 'v') jbe = jde
2134       IF (variable == 'v') jtf = min(jte,jde)
2135       IF (variable == 'm') ktf = kde
2136       IF (variable == 'h') ktf = kde
2137       IF (variable == 'w') ktf = kde
2139       msfcouple = .false.
2140       mucouple = .true.
2141       IF (variable == 'u' .OR. variable == 'v' .OR. variable == 'w')msfcouple = .true.
2142       IF (variable == 'm' )mucouple = .false.
2143       xmsf = 1.
2144       xmu = 1.
2146       IF (jts - jbs .lt. spec_zone) THEN
2147 ! Y-start boundary
2148         DO j = jts, min(jtf,jbs+spec_zone-1)
2149           b_dist = j - jbs
2150           b_limit = b_dist
2151           IF(periodic_x)b_limit = 0
2152           DO k = kts, ktf
2153             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2154               bfield = field_bdy_ys(i, k, b_dist+1) &
2155                    + dtbc * field_bdy_tend_ys(i, k, b_dist+1)
2156               if(msfcouple)xmsf = msf(i,j)
2157               if(mucouple)xmu = (c1(k)*mu(i,j)+c2(k))
2158               field(i,k,j) = xmsf*bfield/xmu
2159             ENDDO
2160           ENDDO
2161         ENDDO
2162       ENDIF
2164       IF (jbe - jtf .lt. spec_zone) THEN
2165 ! Y-end boundary
2166         DO j = max(jts,jbe-spec_zone+1), jtf
2167           b_dist = jbe - j
2168           b_limit = b_dist
2169           IF(periodic_x)b_limit = 0
2170           DO k = kts, ktf
2171             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2172               bfield = field_bdy_ye(i, k, b_dist+1) &
2173                    + dtbc * field_bdy_tend_ye(i, k, b_dist+1)
2174               if(msfcouple)xmsf = msf(i,j)
2175               if(mucouple)xmu = (c1(k)*mu(i,j)+c2(k))
2176               field(i,k,j) = xmsf*bfield/xmu
2177             ENDDO
2178           ENDDO
2179         ENDDO
2180       ENDIF
2182     IF(.NOT.periodic_x)THEN
2183       IF (its - ibs .lt. spec_zone) THEN
2184 ! X-start boundary
2185         DO i = its, min(itf,ibs+spec_zone-1)
2186           b_dist = i - ibs
2187           DO k = kts, ktf
2188             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2189               bfield = field_bdy_xs(j, k, b_dist+1) &
2190                    + dtbc * field_bdy_tend_xs(j, k, b_dist+1)
2191               if(msfcouple)xmsf = msf(i,j)
2192               if(mucouple)xmu = (c1(k)*mu(i,j)+c2(k))
2193               field(i,k,j) = xmsf*bfield/xmu
2194             ENDDO
2195           ENDDO
2196         ENDDO
2197       ENDIF
2199      IF (ibe - itf .lt. spec_zone) THEN
2200 ! X-end boundary
2201         DO i = max(its,ibe-spec_zone+1), itf
2202           b_dist = ibe - i
2203           DO k = kts, ktf
2204             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2205               bfield = field_bdy_xe(j, k, b_dist+1) &
2206                    + dtbc * field_bdy_tend_xe(j, k, b_dist+1)
2207               if(msfcouple)xmsf = msf(i,j)
2208               if(mucouple)xmu = (c1(k)*mu(i,j)+c2(k))
2209               field(i,k,j) = xmsf*bfield/xmu
2210             ENDDO
2211           ENDDO
2212         ENDDO
2213       ENDIF
2214     ENDIF
2216    END SUBROUTINE spec_bdy_final
2217 !------------------------------------------------------------------------
2219    SUBROUTINE zero_grad_bdy (  field,                     &
2220                                variable_in, config_flags, &
2221                                spec_zone,                  &
2222                                ids,ide, jds,jde, kds,kde,  & ! domain dims
2223                                ims,ime, jms,jme, kms,kme,  & ! memory dims
2224                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
2225                                its,ite, jts,jte, kts,kte )
2227 !  This subroutine sets zero gradient conditions in the boundary specified region.
2228 !  spec_zone is the width of the outer specified b.c.s that are set here.
2229 !  (JD August 2000)
2231       IMPLICIT NONE
2233       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
2234       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
2235       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
2236       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
2237       INTEGER,      INTENT(IN   )    :: spec_zone
2238       CHARACTER,    INTENT(IN   )    :: variable_in
2241       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field
2242       TYPE( grid_config_rec_type ) config_flags
2244       CHARACTER  :: variable
2245       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf, i_inner, j_inner
2246       INTEGER    :: b_dist, b_limit
2247       LOGICAL    :: periodic_x
2249       periodic_x = config_flags%periodic_x
2251       variable = variable_in
2253       IF (variable == 'U') variable = 'u'
2254       IF (variable == 'V') variable = 'v'
2256       ibs = ids
2257       ibe = ide-1
2258       itf = min(ite,ide-1)
2259       jbs = jds
2260       jbe = jde-1
2261       jtf = min(jte,jde-1)
2262       ktf = kde-1
2263       IF (variable == 'u') ibe = ide
2264       IF (variable == 'u') itf = min(ite,ide)
2265       IF (variable == 'v') jbe = jde
2266       IF (variable == 'v') jtf = min(jte,jde)
2267       IF (variable == 'w') ktf = kde
2269       IF (jts - jbs .lt. spec_zone) THEN
2270 ! Y-start boundary
2271         DO j = jts, min(jtf,jbs+spec_zone-1)
2272           b_dist = j - jbs
2273           b_limit = b_dist
2274           IF(periodic_x)b_limit = 0
2275           DO k = kts, ktf
2276             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2277               i_inner = max(i,ibs+spec_zone)
2278               i_inner = min(i_inner,ibe-spec_zone)
2279               IF(periodic_x)i_inner = i
2280               field(i,k,j) = field(i_inner,k,jbs+spec_zone)
2281             ENDDO
2282           ENDDO
2283         ENDDO
2284       ENDIF
2285       IF (jbe - jtf .lt. spec_zone) THEN
2286 ! Y-end boundary
2287         DO j = max(jts,jbe-spec_zone+1), jtf
2288           b_dist = jbe - j
2289           b_limit = b_dist
2290           IF(periodic_x)b_limit = 0
2291           DO k = kts, ktf
2292             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2293               i_inner = max(i,ibs+spec_zone)
2294               i_inner = min(i_inner,ibe-spec_zone)
2295               IF(periodic_x)i_inner = i
2296               field(i,k,j) = field(i_inner,k,jbe-spec_zone)
2297             ENDDO
2298           ENDDO
2299         ENDDO
2300       ENDIF
2302     IF(.NOT.periodic_x)THEN
2303       IF (its - ibs .lt. spec_zone) THEN
2304 ! X-start boundary
2305         DO i = its, min(itf,ibs+spec_zone-1)
2306           b_dist = i - ibs
2307           DO k = kts, ktf
2308             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2309               j_inner = max(j,jbs+spec_zone)
2310               j_inner = min(j_inner,jbe-spec_zone)
2311               field(i,k,j) = field(ibs+spec_zone,k,j_inner)
2312             ENDDO
2313           ENDDO
2314         ENDDO
2315       ENDIF
2317       IF (ibe - itf .lt. spec_zone) THEN
2318 ! X-end boundary
2319         DO i = max(its,ibe-spec_zone+1), itf
2320           b_dist = ibe - i
2321           DO k = kts, ktf
2322             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2323               j_inner = max(j,jbs+spec_zone)
2324               j_inner = min(j_inner,jbe-spec_zone)
2325               field(i,k,j) = field(ibe-spec_zone,k,j_inner)
2326             ENDDO
2327           ENDDO
2328         ENDDO
2329       ENDIF
2330     ENDIF
2332    END SUBROUTINE zero_grad_bdy
2333 !------------------------------------------------------------------------
2335    SUBROUTINE flow_dep_bdy  (  field,                     &
2336                                u, v, config_flags, &
2337                                spec_zone,                  &
2338                                ids,ide, jds,jde, kds,kde,  & ! domain dims
2339                                ims,ime, jms,jme, kms,kme,  & ! memory dims
2340                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
2341                                its,ite, jts,jte, kts,kte )
2343 !  This subroutine sets zero gradient conditions for outflow and zero value
2344 !  for inflow in the boundary specified region. Note that field must be unstaggered.
2345 !  The velocities, u and v, will only be used to check their sign (coupled vels OK)
2346 !  spec_zone is the width of the outer specified b.c.s that are set here.
2347 !  (JD August 2000)
2349       IMPLICIT NONE
2351       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
2352       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
2353       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
2354       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
2355       INTEGER,      INTENT(IN   )    :: spec_zone
2358       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field
2359       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: u
2360       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: v
2361       TYPE( grid_config_rec_type ) config_flags
2363       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf, i_inner, j_inner
2364       INTEGER    :: b_dist, b_limit
2365       LOGICAL    :: periodic_x
2367       periodic_x = config_flags%periodic_x
2369       ibs = ids
2370       ibe = ide-1
2371       itf = min(ite,ide-1)
2372       jbs = jds
2373       jbe = jde-1
2374       jtf = min(jte,jde-1)
2375       ktf = kde-1
2377       IF (jts - jbs .lt. spec_zone) THEN
2378 ! Y-start boundary
2379         DO j = jts, min(jtf,jbs+spec_zone-1)
2380           b_dist = j - jbs
2381           b_limit = b_dist
2382           IF(periodic_x)b_limit = 0
2383           DO k = kts, ktf
2384             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2385               i_inner = max(i,ibs+spec_zone)
2386               i_inner = min(i_inner,ibe-spec_zone)
2387               IF(periodic_x)i_inner = i
2388               IF(v(i,k,j) .lt. 0.)THEN
2389                 field(i,k,j) = field(i_inner,k,jbs+spec_zone)
2390               ELSE
2391                 field(i,k,j) = 0.
2392               ENDIF
2393             ENDDO
2394           ENDDO
2395         ENDDO
2396       ENDIF
2397       IF (jbe - jtf .lt. spec_zone) THEN
2398 ! Y-end boundary
2399         DO j = max(jts,jbe-spec_zone+1), jtf
2400           b_dist = jbe - j
2401           b_limit = b_dist
2402           IF(periodic_x)b_limit = 0
2403           DO k = kts, ktf
2404             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2405               i_inner = max(i,ibs+spec_zone)
2406               i_inner = min(i_inner,ibe-spec_zone)
2407               IF(periodic_x)i_inner = i
2408               IF(v(i,k,j+1) .gt. 0.)THEN
2409                 field(i,k,j) = field(i_inner,k,jbe-spec_zone)
2410               ELSE
2411                 field(i,k,j) = 0.
2412               ENDIF
2413             ENDDO
2414           ENDDO
2415         ENDDO
2416       ENDIF
2418     IF(.NOT.periodic_x)THEN
2419       IF (its - ibs .lt. spec_zone) THEN
2420 ! X-start boundary
2421         DO i = its, min(itf,ibs+spec_zone-1)
2422           b_dist = i - ibs
2423           DO k = kts, ktf
2424             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2425               j_inner = max(j,jbs+spec_zone)
2426               j_inner = min(j_inner,jbe-spec_zone)
2427               IF(u(i,k,j) .lt. 0.)THEN
2428                 field(i,k,j) = field(ibs+spec_zone,k,j_inner)
2429               ELSE
2430                 field(i,k,j) = 0.
2431               ENDIF
2432             ENDDO
2433           ENDDO
2434         ENDDO
2435       ENDIF
2437       IF (ibe - itf .lt. spec_zone) THEN
2438 ! X-end boundary
2439         DO i = max(its,ibe-spec_zone+1), itf
2440           b_dist = ibe - i
2441           DO k = kts, ktf
2442             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2443               j_inner = max(j,jbs+spec_zone)
2444               j_inner = min(j_inner,jbe-spec_zone)
2445               IF(u(i+1,k,j) .gt. 0.)THEN
2446                 field(i,k,j) = field(ibe-spec_zone,k,j_inner)
2447               ELSE
2448                 field(i,k,j) = 0.
2449               ENDIF
2450             ENDDO
2451           ENDDO
2452         ENDDO
2453       ENDIF
2454     ENDIF
2456    END SUBROUTINE flow_dep_bdy
2458 !------------------------------------------------------------------------------
2460    SUBROUTINE flow_dep_bdy_qnn  (  field,          &
2461                                u, v, config_flags, &
2462                                spec_zone,                  &
2463                                ccn_conc,                   & ! RAS
2464                                ids,ide, jds,jde, kds,kde,  & ! domain dims
2465                                ims,ime, jms,jme, kms,kme,  & ! memory dims
2466                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
2467                                its,ite, jts,jte, kts,kte )
2469 !  This subroutine sets zero gradient conditions for outflow and zero value
2470 !  for inflow in the boundary specified region. Note that field must be unstaggered.
2471 !  The velocities, u and v, will only be used to check their sign (coupled vels OK)
2472 !  spec_zone is the width of the outer specified b.c.s that are set here.
2473 !  (JD August 2000)
2475       IMPLICIT NONE
2477       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
2478       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
2479       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
2480       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
2481       INTEGER,      INTENT(IN   )    :: spec_zone
2482       REAL,         INTENT(IN   )    :: ccn_conc ! RAS
2485       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field
2486       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: u
2487       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: v
2488       TYPE( grid_config_rec_type ) config_flags
2490       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf, i_inner, j_inner
2491       INTEGER    :: b_dist, b_limit
2492       LOGICAL    :: periodic_x
2494       periodic_x = config_flags%periodic_x
2496       ibs = ids
2497       ibe = ide-1
2498       itf = min(ite,ide-1)
2499       jbs = jds
2500       jbe = jde-1
2501       jtf = min(jte,jde-1)
2502       ktf = kde-1
2504       IF (jts - jbs .lt. spec_zone) THEN
2505 ! Y-start boundary
2506         DO j = jts, min(jtf,jbs+spec_zone-1)
2507           b_dist = j - jbs
2508           b_limit = b_dist
2509           IF(periodic_x)b_limit = 0
2510           DO k = kts, ktf
2511             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2512               i_inner = max(i,ibs+spec_zone)
2513               i_inner = min(i_inner,ibe-spec_zone)
2514               IF(periodic_x)i_inner = i
2515               IF(v(i,k,j) .lt. 0.)THEN
2516                 field(i,k,j) = field(i_inner,k,jbs+spec_zone)
2517               ELSE
2518                 field(i,k,j) = ccn_conc ! RAS
2519               ENDIF
2520             ENDDO
2521           ENDDO
2522         ENDDO
2523       ENDIF   
2524       IF (jbe - jtf .lt. spec_zone) THEN
2525 ! Y-end boundary
2526         DO j = max(jts,jbe-spec_zone+1), jtf
2527           b_dist = jbe - j
2528           b_limit = b_dist
2529           IF(periodic_x)b_limit = 0
2530           DO k = kts, ktf
2531             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2532               i_inner = max(i,ibs+spec_zone)
2533               i_inner = min(i_inner,ibe-spec_zone)
2534               IF(periodic_x)i_inner = i
2535               IF(v(i,k,j+1) .gt. 0.)THEN
2536                 field(i,k,j) = field(i_inner,k,jbe-spec_zone)
2537               ELSE
2538                 field(i,k,j) = ccn_conc ! RAS
2539               ENDIF
2540             ENDDO
2541           ENDDO
2542         ENDDO
2543       ENDIF
2545     IF(.NOT.periodic_x)THEN
2546       IF (its - ibs .lt. spec_zone) THEN
2547 ! X-start boundary
2548         DO i = its, min(itf,ibs+spec_zone-1)
2549           b_dist = i - ibs
2550           DO k = kts, ktf
2551             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2552               j_inner = max(j,jbs+spec_zone)
2553               j_inner = min(j_inner,jbe-spec_zone)
2554               IF(u(i,k,j) .lt. 0.)THEN
2555                 field(i,k,j) = field(ibs+spec_zone,k,j_inner)
2556               ELSE
2557                 field(i,k,j) = ccn_conc ! RAS
2558               ENDIF
2559             ENDDO
2560           ENDDO
2561         ENDDO
2562       ENDIF
2564       IF (ibe - itf .lt. spec_zone) THEN
2565 ! X-end boundary
2566         DO i = max(its,ibe-spec_zone+1), itf
2567           b_dist = ibe - i
2568           DO k = kts, ktf
2569             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2570               j_inner = max(j,jbs+spec_zone)
2571               j_inner = min(j_inner,jbe-spec_zone)
2572               IF(u(i+1,k,j) .gt. 0.)THEN
2573                 field(i,k,j) = field(ibe-spec_zone,k,j_inner)
2574               ELSE
2575                 field(i,k,j) = ccn_conc ! RAS
2576               ENDIF
2577             ENDDO
2578           ENDDO
2579         ENDDO
2580       ENDIF
2581     ENDIF
2583    END SUBROUTINE flow_dep_bdy_qnn
2585 !-------------------------for ntu3m-------------------------------------------------
2586    SUBROUTINE flow_dep_bdy_fixed_inflow(field,u,v,config_flags,        &
2587                                spec_zone,ids,ide,jds,jde,kds,kde,ims,  &
2588                                ime,jms,jme,kms,kme,ips,ipe,jps,jpe,kps,&
2589                                kpe,its,ite,jts,jte,kts,kte)
2590 !-------------------------------------------------------------------------------
2591       IMPLICIT NONE
2592       INTEGER, INTENT(IN) :: ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,  &
2593                              kms,kme,ips,ipe,jps,jpe,kps,kpe,its,ite,  &
2594                              jts,jte,kts,kte,spec_zone
2595       REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: field
2596       REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN) :: u,v
2597       TYPE( grid_config_rec_type ) config_flags
2598       INTEGER :: i,j,k,ibs,ibe,jbs,jbe,itf,jtf,ktf,i_inner,j_inner,    &
2599                  b_dist,b_limit
2600       LOGICAL :: periodic_x
2602       periodic_x = config_flags%periodic_x
2603       ibs = ids
2604       ibe = ide-1
2605       itf = min(ite,ide-1)
2606       jbs = jds
2607       jbe = jde-1
2608       jtf = min(jte,jde-1)
2609       ktf = kde-1
2611       IF (jts - jbs .lt. spec_zone) THEN
2612 ! Y-start boundary
2613         DO j = jts, min(jtf,jbs+spec_zone-1)
2614           b_dist = j - jbs
2615           b_limit = b_dist
2616           IF (periodic_x) b_limit = 0
2617           DO k = kts, ktf
2618             DO i = max(its,b_limit+ibs),min(itf,ibe-b_limit)
2619               i_inner = max(i,ibs+spec_zone)
2620               i_inner = min(i_inner,ibe-spec_zone)
2621               IF (periodic_x) i_inner = i
2622               IF (v(i,k,j) .lt. 0.) THEN
2623                 field(i,k,j) = field(i_inner,k,jbs+spec_zone)
2624               ELSE
2625                 field(i,k,j) = field(i,k,j)
2626               ENDIF
2627             ENDDO
2628           ENDDO
2629         ENDDO
2630       ENDIF
2631       IF (jbe - jtf .lt. spec_zone) THEN
2632 ! Y-end boundary
2633         DO j = max(jts,jbe-spec_zone+1),jtf
2634           b_dist = jbe - j
2635           b_limit = b_dist
2636           IF (periodic_x) b_limit = 0
2637           DO k = kts,ktf
2638             DO i = max(its,b_limit+ibs),min(itf,ibe-b_limit)
2639               i_inner = max(i,ibs+spec_zone)
2640               i_inner = min(i_inner,ibe-spec_zone)
2641               IF (periodic_x) i_inner = i
2642               IF (v(i,k,j+1) .gt. 0.) THEN
2643                 field(i,k,j) = field(i_inner,k,jbe-spec_zone)
2644               ELSE
2645                 field(i,k,j) = field(i,k,j)
2646               ENDIF
2647             ENDDO
2648           ENDDO
2649         ENDDO
2650       ENDIF
2652       IF (.NOT.periodic_x) THEN
2653         IF (its - ibs .lt. spec_zone) THEN
2654 ! X-start boundary
2655           DO i = its, min(itf,ibs+spec_zone-1)
2656             b_dist = i - ibs
2657             DO k = kts, ktf
2658               DO j = max(jts,b_dist+jbs+1),min(jtf,jbe-b_dist-1)
2659                 j_inner = max(j,jbs+spec_zone)
2660                 j_inner = min(j_inner,jbe-spec_zone)
2661                 IF (u(i,k,j) .lt. 0.) THEN
2662                   field(i,k,j) = field(ibs+spec_zone,k,j_inner)
2663                 ELSE
2664                   field(i,k,j) = field(i,k,j)
2665                 ENDIF
2666               ENDDO
2667             ENDDO
2668           ENDDO
2669         ENDIF
2670         IF (ibe - itf .lt. spec_zone) THEN
2671 ! X-end boundary
2672           DO i = max(its,ibe-spec_zone+1),itf
2673             b_dist = ibe - i
2674             DO k = kts, ktf
2675               DO j = max(jts,b_dist+jbs+1),min(jtf,jbe-b_dist-1)
2676                 j_inner = max(j,jbs+spec_zone)
2677                 j_inner = min(j_inner,jbe-spec_zone)
2678                 IF (u(i+1,k,j) .gt. 0.) THEN
2679                   field(i,k,j) = field(ibe-spec_zone,k,j_inner)
2680                 ELSE
2681                   field(i,k,j) = field(i,k,j)
2682                 ENDIF
2683               ENDDO
2684             ENDDO
2685           ENDDO
2686         ENDIF
2687       ENDIF
2689    END SUBROUTINE flow_dep_bdy_fixed_inflow
2690 !----------------------------for ntu3m--------------------------------------------
2692 !------------------------------------------------------------------------------
2694  SUBROUTINE stuff_bdy_new ( data3d , space_bdy_xs, space_bdy_xe, space_bdy_ys, space_bdy_ye, &
2695                              char_stagger , &
2696                              spec_bdy_width , &
2697                              ids, ide, jds, jde, kds, kde , &
2698                              ims, ime, jms, jme, kms, kme , &
2699                              its, ite, jts, jte, kts, kte )
2701  !  This routine puts the data in the 3d arrays into the proper locations
2702  !  for the lateral boundary arrays.
2704     USE module_state_description
2705     
2706     IMPLICIT NONE
2708     INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde
2709     INTEGER , INTENT(IN) :: ims, ime, jms, jme, kms, kme
2710     INTEGER , INTENT(IN) :: its, ite, jts, jte, kts, kte
2711     INTEGER , INTENT(IN) :: spec_bdy_width
2712     REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: data3d
2713     REAL , DIMENSION(jms:jme,kds:kde,spec_bdy_width) , INTENT(OUT) :: space_bdy_xs, space_bdy_xe
2714     REAL , DIMENSION(ims:ime,kds:kde,spec_bdy_width) , INTENT(OUT) :: space_bdy_ys, space_bdy_ye
2715     CHARACTER (LEN=1) , INTENT(IN) :: char_stagger
2717     INTEGER :: i , ii , j , jj , k
2719     !  There are four lateral boundary locations that are stored.
2721     !  X start boundary
2723     IF ( char_stagger .EQ. 'W' ) THEN
2724        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2725        DO k = kds , kde
2726        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2727           space_bdy_xs(j,k,i) = data3d(i,k,j)
2728        END DO
2729        END DO
2730        END DO
2731     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2732        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2733        DO k = kds , kde
2734        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2735           space_bdy_xs(j,k,i) = data3d(i,k,j)
2736        END DO
2737        END DO
2738        END DO
2739     ELSE IF ( char_stagger .EQ. 'V' ) THEN
2740        DO j = MAX(jds,jts) , MIN(jde,jte)
2741        DO k = kds , kde - 1
2742        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2743           space_bdy_xs(j,k,i) = data3d(i,k,j)
2744        END DO
2745        END DO
2746        END DO
2747     ELSE
2748        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2749        DO k = kds , kde - 1
2750        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2751           space_bdy_xs(j,k,i) = data3d(i,k,j)
2752        END DO
2753        END DO
2754        END DO
2755     END IF
2757     !  X end boundary
2759     IF      ( char_stagger .EQ. 'U' ) THEN
2760        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2761        DO k = kds , kde - 1
2762        DO i = MIN(ide,ite) , MAX(ide - spec_bdy_width + 1,its) , -1
2763           ii = ide - i + 1
2764           space_bdy_xe(j,k,ii) = data3d(i,k,j)
2765        END DO
2766        END DO
2767        END DO
2768     ELSE IF ( char_stagger .EQ. 'V' ) THEN
2769        DO j = MAX(jds,jts) , MIN(jde,jte)
2770        DO k = kds , kde - 1
2771        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2772           ii = ide - i
2773           space_bdy_xe(j,k,ii) = data3d(i,k,j)
2774        END DO
2775        END DO
2776        END DO
2777     ELSE IF ( char_stagger .EQ. 'W' ) THEN
2778        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2779        DO k = kds , kde
2780        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2781           ii = ide - i
2782           space_bdy_xe(j,k,ii) = data3d(i,k,j)
2783        END DO
2784        END DO
2785        END DO
2786     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2787        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2788        DO k = kds , kde
2789        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2790           ii = ide - i
2791           space_bdy_xe(j,k,ii) = data3d(i,k,j)
2792        END DO
2793        END DO
2794        END DO
2795     ELSE
2796        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2797        DO k = kds , kde - 1
2798        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2799           ii = ide - i
2800           space_bdy_xe(j,k,ii) = data3d(i,k,j)
2801        END DO
2802        END DO
2803        END DO
2804     END IF
2806     !  Y start boundary
2808     IF ( char_stagger .EQ. 'W' ) THEN
2809        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2810        DO k = kds , kde
2811        DO i = MAX(ids,its) , MIN(ide-1,ite)
2812           space_bdy_ys(i,k,j) = data3d(i,k,j)
2813        END DO
2814        END DO
2815        END DO
2816     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2817        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2818        DO k = kds , kde
2819        DO i = MAX(ids,its) , MIN(ide-1,ite)
2820           space_bdy_ys(i,k,j) = data3d(i,k,j)
2821        END DO
2822        END DO
2823        END DO
2824     ELSE IF ( char_stagger .EQ. 'U' ) THEN
2825        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2826        DO k = kds , kde - 1
2827        DO i = MAX(ids,its) , MIN(ide,ite)
2828           space_bdy_ys(i,k,j) = data3d(i,k,j)
2829        END DO
2830        END DO
2831        END DO
2832     ELSE
2833        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
2834        DO k = kds , kde - 1
2835        DO i = MAX(ids,its) , MIN(ide-1,ite)
2836           space_bdy_ys(i,k,j) = data3d(i,k,j)
2837        END DO
2838        END DO
2839        END DO
2840     END IF
2842     !  Y end boundary
2844     IF      ( char_stagger .EQ. 'V' ) THEN
2845        DO j = MIN(jde,jte) , MAX(jde - spec_bdy_width + 1,jts) , -1
2846        DO k = kds , kde - 1
2847        DO i = MAX(ids,its) , MIN(ide-1,ite)
2848           jj = jde - j + 1
2849           space_bdy_ye(i,k,jj) = data3d(i,k,j)
2850        END DO
2851        END DO
2852        END DO
2853     ELSE IF ( char_stagger .EQ. 'U' ) THEN
2854        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2855        DO k = kds , kde - 1
2856        DO i = MAX(ids,its) , MIN(ide,ite)
2857           jj = jde - j
2858           space_bdy_ye(i,k,jj) = data3d(i,k,j)
2859        END DO
2860        END DO
2861        END DO
2862     ELSE IF ( char_stagger .EQ. 'W' ) THEN
2863        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2864        DO k = kds , kde
2865        DO i = MAX(ids,its) , MIN(ide-1,ite)
2866           jj = jde - j
2867           space_bdy_ye(i,k,jj) = data3d(i,k,j)
2868        END DO
2869        END DO
2870        END DO
2871     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2872        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2873        DO k = kds , kde
2874        DO i = MAX(ids,its) , MIN(ide-1,ite)
2875           jj = jde - j
2876           space_bdy_ye(i,k,jj) = data3d(i,k,j)
2877        END DO
2878        END DO
2879        END DO
2880     ELSE
2881        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
2882        DO k = kds , kde - 1
2883        DO i = MAX(ids,its) , MIN(ide-1,ite)
2884           jj = jde - j
2885           space_bdy_ye(i,k,jj) = data3d(i,k,j)
2886        END DO
2887        END DO
2888        END DO
2889     END IF
2890     
2891  END SUBROUTINE stuff_bdy_new
2893  SUBROUTINE stuff_bdytend_new ( data3dnew , data3dold , time_diff , &
2894                              space_bdy_xs, space_bdy_xe, space_bdy_ys, space_bdy_ye, &
2895                              char_stagger , &
2896                              spec_bdy_width , &
2897                              ids, ide, jds, jde, kds, kde , &
2898                              ims, ime, jms, jme, kms, kme , &
2899                              its, ite, jts, jte, kts, kte )
2901  !  This routine puts the tendency data into the proper locations
2902  !  for the lateral boundary arrays.
2904     USE module_state_description
2905     
2906     IMPLICIT NONE
2908     INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde
2909     INTEGER , INTENT(IN) :: ims, ime, jms, jme, kms, kme
2910     INTEGER , INTENT(IN) :: its, ite, jts, jte, kts, kte
2911     INTEGER , INTENT(IN) :: spec_bdy_width
2912     REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: data3dnew , data3dold
2913     REAL , DIMENSION(jms:jme,kds:kde,spec_bdy_width) , INTENT(OUT) :: space_bdy_xs, space_bdy_xe
2914     REAL , DIMENSION(ims:ime,kds:kde,spec_bdy_width) , INTENT(OUT) :: space_bdy_ys, space_bdy_ye
2915     CHARACTER (LEN=1) , INTENT(IN) :: char_stagger
2916     REAL , INTENT(IN) :: time_diff ! seconds
2918     INTEGER :: i , ii , j , jj , k
2920 #if 0
2921     !  Here is the easy way to zero out the boundary tendencies for
2922     !  time-dependent boundaries.
2924     space_bdy_xs = 0.
2925     space_bdy_xe = 0.
2926     space_bdy_ys = 0.
2927     space_bdy_ye = 0.
2928     return
2929 #else
2930     !  There are four lateral boundary locations that are stored.
2932     !  X start boundary
2934     IF ( char_stagger .EQ. 'W' ) THEN
2935        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2936        DO k = kds , kde
2937        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2938           space_bdy_xs(j,k,i) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2939        END DO
2940        END DO
2941        END DO
2942     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2943        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2944        DO k = kds , kde
2945        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2946           space_bdy_xs(j,k,i) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2947        END DO
2948        END DO
2949        END DO
2950     ELSE IF ( char_stagger .EQ. 'V' ) THEN
2951        DO j = MAX(jds,jts) , MIN(jde,jte)
2952        DO k = kds , kde - 1
2953        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2954           space_bdy_xs(j,k,i) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2955        END DO
2956        END DO
2957        END DO
2958     ELSE
2959        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2960        DO k = kds , kde - 1
2961        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
2962           space_bdy_xs(j,k,i) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2963        END DO
2964        END DO
2965        END DO
2966     END IF
2968     !  X end boundary
2970     IF      ( char_stagger .EQ. 'U' ) THEN
2971        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2972        DO k = kds , kde - 1
2973        DO i = MIN(ide,ite) , MAX(ide - spec_bdy_width + 1,its) , -1
2974           ii = ide - i + 1
2975           space_bdy_xe(j,k,ii) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2976        END DO
2977        END DO
2978        END DO
2979     ELSE IF ( char_stagger .EQ. 'V' ) THEN
2980        DO j = MAX(jds,jts) , MIN(jde,jte)
2981        DO k = kds , kde - 1
2982        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2983           ii = ide - i
2984           space_bdy_xe(j,k,ii) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2985        END DO
2986        END DO
2987        END DO
2988     ELSE IF ( char_stagger .EQ. 'W' ) THEN
2989        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2990        DO k = kds , kde
2991        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
2992           ii = ide - i
2993           space_bdy_xe(j,k,ii) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
2994        END DO
2995        END DO
2996        END DO
2997     ELSE IF ( char_stagger .EQ. 'M' ) THEN
2998        DO j = MAX(jds,jts) , MIN(jde-1,jte)
2999        DO k = kds , kde
3000        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3001           ii = ide - i
3002           space_bdy_xe(j,k,ii) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3003        END DO
3004        END DO
3005        END DO
3006     ELSE
3007        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3008        DO k = kds , kde - 1
3009        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3010           ii = ide - i
3011           space_bdy_xe(j,k,ii) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3012        END DO
3013        END DO
3014        END DO
3015     END IF
3017     !  Y start boundary
3019     IF ( char_stagger .EQ. 'W' ) THEN
3020        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3021        DO k = kds , kde
3022        DO i = MAX(ids,its) , MIN(ide-1,ite)
3023           space_bdy_ys(i,k,j) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3024        END DO
3025        END DO
3026        END DO
3027     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3028        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3029        DO k = kds , kde
3030        DO i = MAX(ids,its) , MIN(ide-1,ite)
3031           space_bdy_ys(i,k,j) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3032        END DO
3033        END DO
3034        END DO
3035     ELSE IF ( char_stagger .EQ. 'U' ) THEN
3036        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3037        DO k = kds , kde - 1
3038        DO i = MAX(ids,its) , MIN(ide,ite)
3039           space_bdy_ys(i,k,j) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3040        END DO
3041        END DO
3042        END DO
3043     ELSE
3044        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3045        DO k = kds , kde - 1
3046        DO i = MAX(ids,its) , MIN(ide-1,ite)
3047           space_bdy_ys(i,k,j) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3048        END DO
3049        END DO
3050        END DO
3051     END IF
3053     !  Y end boundary
3055     IF      ( char_stagger .EQ. 'V' ) THEN
3056        DO j = MIN(jde,jte) , MAX(jde - spec_bdy_width + 1,jts) , -1
3057        DO k = kds , kde - 1
3058        DO i = MAX(ids,its) , MIN(ide-1,ite)
3059           jj = jde - j + 1
3060           space_bdy_ye(i,k,jj) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3061        END DO
3062        END DO
3063        END DO
3064     ELSE IF ( char_stagger .EQ. 'U' ) THEN
3065        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3066        DO k = kds , kde - 1
3067        DO i = MAX(ids,its) , MIN(ide,ite)
3068           jj = jde - j
3069           space_bdy_ye(i,k,jj) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3070        END DO
3071        END DO
3072        END DO
3073     ELSE IF ( char_stagger .EQ. 'W' ) THEN
3074        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3075        DO k = kds , kde
3076        DO i = MAX(ids,its) , MIN(ide-1,ite)
3077           jj = jde - j
3078           space_bdy_ye(i,k,jj) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3079        END DO
3080        END DO
3081        END DO
3082     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3083        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3084        DO k = kds , kde
3085        DO i = MAX(ids,its) , MIN(ide-1,ite)
3086           jj = jde - j
3087           space_bdy_ye(i,k,jj) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3088        END DO
3089        END DO
3090        END DO
3091     ELSE
3092        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3093        DO k = kds , kde - 1
3094        DO i = MAX(ids,its) , MIN(ide-1,ite)
3095           jj = jde - j
3096           space_bdy_ye(i,k,jj) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3097        END DO
3098        END DO
3099        END DO
3100     END IF
3101 #endif
3102     
3103  END SUBROUTINE stuff_bdytend_new
3105 !--- old versions for use with modules that use the old bdy data structures ---
3107  SUBROUTINE stuff_bdy_old ( data3d , space_bdy , char_stagger , &
3108                              ijds , ijde , spec_bdy_width , &
3109                              ids, ide, jds, jde, kds, kde , &
3110                              ims, ime, jms, jme, kms, kme , &
3111                              its, ite, jts, jte, kts, kte )
3113  !  This routine puts the data in the 3d arrays into the proper locations
3114  !  for the lateral boundary arrays.
3116     USE module_state_description
3117     
3118     IMPLICIT NONE
3120     INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde
3121     INTEGER , INTENT(IN) :: ims, ime, jms, jme, kms, kme
3122     INTEGER , INTENT(IN) :: its, ite, jts, jte, kts, kte
3123     INTEGER , INTENT(IN) :: ijds , ijde , spec_bdy_width
3124     REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: data3d
3125     REAL , DIMENSION(ijds:ijde,kds:kde,spec_bdy_width,4) , INTENT(OUT) :: space_bdy
3126     CHARACTER (LEN=1) , INTENT(IN) :: char_stagger
3128     INTEGER :: i , ii , j , jj , k
3130     !  There are four lateral boundary locations that are stored.
3132     !  X start boundary
3134     IF ( char_stagger .EQ. 'W' ) THEN
3135        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3136        DO k = kds , kde
3137        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3138           space_bdy(j,k,i,P_XSB) = data3d(i,k,j)
3139        END DO
3140        END DO
3141        END DO
3142     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3143        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3144        DO k = kds , kde
3145        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3146           space_bdy(j,k,i,P_XSB) = data3d(i,k,j)
3147        END DO
3148        END DO
3149        END DO
3150     ELSE IF ( char_stagger .EQ. 'V' ) THEN
3151        DO j = MAX(jds,jts) , MIN(jde,jte)
3152        DO k = kds , kde - 1
3153        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3154           space_bdy(j,k,i,P_XSB) = data3d(i,k,j)
3155        END DO
3156        END DO
3157        END DO
3158     ELSE
3159        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3160        DO k = kds , kde - 1
3161        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3162           space_bdy(j,k,i,P_XSB) = data3d(i,k,j)
3163        END DO
3164        END DO
3165        END DO
3166     END IF
3168     !  X end boundary
3170     IF      ( char_stagger .EQ. 'U' ) THEN
3171        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3172        DO k = kds , kde - 1
3173        DO i = MIN(ide,ite) , MAX(ide - spec_bdy_width + 1,its) , -1
3174           ii = ide - i + 1
3175           space_bdy(j,k,ii,P_XEB) = data3d(i,k,j)
3176        END DO
3177        END DO
3178        END DO
3179     ELSE IF ( char_stagger .EQ. 'V' ) THEN
3180        DO j = MAX(jds,jts) , MIN(jde,jte)
3181        DO k = kds , kde - 1
3182        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3183           ii = ide - i
3184           space_bdy(j,k,ii,P_XEB) = data3d(i,k,j)
3185        END DO
3186        END DO
3187        END DO
3188     ELSE IF ( char_stagger .EQ. 'W' ) THEN
3189        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3190        DO k = kds , kde
3191        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3192           ii = ide - i
3193           space_bdy(j,k,ii,P_XEB) = data3d(i,k,j)
3194        END DO
3195        END DO
3196        END DO
3197     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3198        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3199        DO k = kds , kde
3200        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3201           ii = ide - i
3202           space_bdy(j,k,ii,P_XEB) = data3d(i,k,j)
3203        END DO
3204        END DO
3205        END DO
3206     ELSE
3207        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3208        DO k = kds , kde - 1
3209        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3210           ii = ide - i
3211           space_bdy(j,k,ii,P_XEB) = data3d(i,k,j)
3212        END DO
3213        END DO
3214        END DO
3215     END IF
3217     !  Y start boundary
3219     IF ( char_stagger .EQ. 'W' ) THEN
3220        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3221        DO k = kds , kde
3222        DO i = MAX(ids,its) , MIN(ide-1,ite)
3223           space_bdy(i,k,j,P_YSB) = data3d(i,k,j)
3224        END DO
3225        END DO
3226        END DO
3227     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3228        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3229        DO k = kds , kde
3230        DO i = MAX(ids,its) , MIN(ide-1,ite)
3231           space_bdy(i,k,j,P_YSB) = data3d(i,k,j)
3232        END DO
3233        END DO
3234        END DO
3235     ELSE IF ( char_stagger .EQ. 'U' ) THEN
3236        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3237        DO k = kds , kde - 1
3238        DO i = MAX(ids,its) , MIN(ide,ite)
3239           space_bdy(i,k,j,P_YSB) = data3d(i,k,j)
3240        END DO
3241        END DO
3242        END DO
3243     ELSE
3244        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3245        DO k = kds , kde - 1
3246        DO i = MAX(ids,its) , MIN(ide-1,ite)
3247           space_bdy(i,k,j,P_YSB) = data3d(i,k,j)
3248        END DO
3249        END DO
3250        END DO
3251     END IF
3253     !  Y end boundary
3255     IF      ( char_stagger .EQ. 'V' ) THEN
3256        DO j = MIN(jde,jte) , MAX(jde - spec_bdy_width + 1,jts) , -1
3257        DO k = kds , kde - 1
3258        DO i = MAX(ids,its) , MIN(ide-1,ite)
3259           jj = jde - j + 1
3260           space_bdy(i,k,jj,P_YEB) = data3d(i,k,j)
3261        END DO
3262        END DO
3263        END DO
3264     ELSE IF ( char_stagger .EQ. 'U' ) THEN
3265        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3266        DO k = kds , kde - 1
3267        DO i = MAX(ids,its) , MIN(ide,ite)
3268           jj = jde - j
3269           space_bdy(i,k,jj,P_YEB) = data3d(i,k,j)
3270        END DO
3271        END DO
3272        END DO
3273     ELSE IF ( char_stagger .EQ. 'W' ) THEN
3274        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3275        DO k = kds , kde
3276        DO i = MAX(ids,its) , MIN(ide-1,ite)
3277           jj = jde - j
3278           space_bdy(i,k,jj,P_YEB) = data3d(i,k,j)
3279        END DO
3280        END DO
3281        END DO
3282     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3283        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3284        DO k = kds , kde
3285        DO i = MAX(ids,its) , MIN(ide-1,ite)
3286           jj = jde - j
3287           space_bdy(i,k,jj,P_YEB) = data3d(i,k,j)
3288        END DO
3289        END DO
3290        END DO
3291     ELSE
3292        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3293        DO k = kds , kde - 1
3294        DO i = MAX(ids,its) , MIN(ide-1,ite)
3295           jj = jde - j
3296           space_bdy(i,k,jj,P_YEB) = data3d(i,k,j)
3297        END DO
3298        END DO
3299        END DO
3300     END IF
3301     
3302  END SUBROUTINE stuff_bdy_old
3304  SUBROUTINE stuff_bdytend_old ( data3dnew , data3dold , time_diff , space_bdy , char_stagger , &
3305                              ijds , ijde , spec_bdy_width , &
3306                              ids, ide, jds, jde, kds, kde , &
3307                              ims, ime, jms, jme, kms, kme , &
3308                              its, ite, jts, jte, kts, kte )
3310  !  This routine puts the tendency data into the proper locations
3311  !  for the lateral boundary arrays.
3313     USE module_state_description
3314     
3315     IMPLICIT NONE
3317     INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde
3318     INTEGER , INTENT(IN) :: ims, ime, jms, jme, kms, kme
3319     INTEGER , INTENT(IN) :: its, ite, jts, jte, kts, kte
3320     INTEGER , INTENT(IN) :: ijds , ijde , spec_bdy_width
3321     REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: data3dnew , data3dold
3322 !    REAL , DIMENSION(:,:,:,:) , INTENT(OUT) :: space_bdy
3323     REAL , DIMENSION(ijds:ijde,kds:kde,spec_bdy_width,4) , INTENT(OUT) :: space_bdy
3324     CHARACTER (LEN=1) , INTENT(IN) :: char_stagger
3325     REAL , INTENT(IN) :: time_diff ! seconds
3327     INTEGER :: i , ii , j , jj , k
3329     !  There are four lateral boundary locations that are stored.
3331     !  X start boundary
3333     IF ( char_stagger .EQ. 'W' ) THEN
3334        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3335        DO k = kds , kde
3336        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3337           space_bdy(j,k,i,P_XSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3338 !         space_bdy(j,k,i,P_XSB) = 0. ! zeroout
3339        END DO
3340        END DO
3341        END DO
3342     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3343        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3344        DO k = kds , kde
3345        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3346           space_bdy(j,k,i,P_XSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3347 !         space_bdy(j,k,i,P_XSB) = 0. ! zeroout
3348        END DO
3349        END DO
3350        END DO
3351     ELSE IF ( char_stagger .EQ. 'V' ) THEN
3352        DO j = MAX(jds,jts) , MIN(jde,jte)
3353        DO k = kds , kde - 1
3354        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3355           space_bdy(j,k,i,P_XSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3356 !         space_bdy(j,k,i,P_XSB) = 0. ! zeroout
3357        END DO
3358        END DO
3359        END DO
3360     ELSE
3361        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3362        DO k = kds , kde - 1
3363        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3364           space_bdy(j,k,i,P_XSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3365 !         space_bdy(j,k,i,P_XSB) = 0. ! zeroout
3366        END DO
3367        END DO
3368        END DO
3369     END IF
3371     !  X end boundary
3373     IF      ( char_stagger .EQ. 'U' ) THEN
3374        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3375        DO k = kds , kde - 1
3376        DO i = MIN(ide,ite) , MAX(ide - spec_bdy_width + 1,its) , -1
3377           ii = ide - i + 1
3378           space_bdy(j,k,ii,P_XEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3379 !         space_bdy(j,k,ii,P_XEB) = 0. ! zeroout
3380        END DO
3381        END DO
3382        END DO
3383     ELSE IF ( char_stagger .EQ. 'V' ) THEN
3384        DO j = MAX(jds,jts) , MIN(jde,jte)
3385        DO k = kds , kde - 1
3386        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3387           ii = ide - i
3388           space_bdy(j,k,ii,P_XEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3389 !         space_bdy(j,k,ii,P_XEB) = 0. ! zeroout
3390        END DO
3391        END DO
3392        END DO
3393     ELSE IF ( char_stagger .EQ. 'W' ) THEN
3394        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3395        DO k = kds , kde
3396        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3397           ii = ide - i
3398           space_bdy(j,k,ii,P_XEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3399 !         space_bdy(j,k,ii,P_XEB) = 0. ! zeroout
3400        END DO
3401        END DO
3402        END DO
3403     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3404        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3405        DO k = kds , kde
3406        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3407           ii = ide - i
3408           space_bdy(j,k,ii,P_XEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3409 !         space_bdy(j,k,ii,P_XEB) = 0. ! zeroout
3410        END DO
3411        END DO
3412        END DO
3413     ELSE
3414        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3415        DO k = kds , kde - 1
3416        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3417           ii = ide - i
3418           space_bdy(j,k,ii,P_XEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3419 !         space_bdy(j,k,ii,P_XEB) = 0. ! zeroout
3420        END DO
3421        END DO
3422        END DO
3423     END IF
3425     !  Y start boundary
3427     IF ( char_stagger .EQ. 'W' ) THEN
3428        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3429        DO k = kds , kde
3430        DO i = MAX(ids,its) , MIN(ide-1,ite)
3431           space_bdy(i,k,j,P_YSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3432 !         space_bdy(i,k,j,P_YSB) = 0. ! zeroout
3433        END DO
3434        END DO
3435        END DO
3436     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3437        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3438        DO k = kds , kde
3439        DO i = MAX(ids,its) , MIN(ide-1,ite)
3440           space_bdy(i,k,j,P_YSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3441 !         space_bdy(i,k,j,P_YSB) = 0. ! zeroout
3442        END DO
3443        END DO
3444        END DO
3445     ELSE IF ( char_stagger .EQ. 'U' ) THEN
3446        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3447        DO k = kds , kde - 1
3448        DO i = MAX(ids,its) , MIN(ide,ite)
3449           space_bdy(i,k,j,P_YSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3450 !         space_bdy(i,k,j,P_YSB) = 0. ! zeroout
3451        END DO
3452        END DO
3453        END DO
3454     ELSE
3455        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3456        DO k = kds , kde - 1
3457        DO i = MAX(ids,its) , MIN(ide-1,ite)
3458           space_bdy(i,k,j,P_YSB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3459 !         space_bdy(i,k,j,P_YSB) = 0. ! zeroout
3460        END DO
3461        END DO
3462        END DO
3463     END IF
3465     !  Y end boundary
3467     IF      ( char_stagger .EQ. 'V' ) THEN
3468        DO j = MIN(jde,jte) , MAX(jde - spec_bdy_width + 1,jts) , -1
3469        DO k = kds , kde - 1
3470        DO i = MAX(ids,its) , MIN(ide-1,ite)
3471           jj = jde - j + 1
3472           space_bdy(i,k,jj,P_YEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3473 !         space_bdy(i,k,jj,P_YEB) = 0. ! zeroout
3474        END DO
3475        END DO
3476        END DO
3477     ELSE IF ( char_stagger .EQ. 'U' ) THEN
3478        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3479        DO k = kds , kde - 1
3480        DO i = MAX(ids,its) , MIN(ide,ite)
3481           jj = jde - j
3482           space_bdy(i,k,jj,P_YEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3483 !         space_bdy(i,k,jj,P_YEB) = 0. ! zeroout
3484        END DO
3485        END DO
3486        END DO
3487     ELSE IF ( char_stagger .EQ. 'W' ) THEN
3488        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3489        DO k = kds , kde
3490        DO i = MAX(ids,its) , MIN(ide-1,ite)
3491           jj = jde - j
3492           space_bdy(i,k,jj,P_YEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3493 !         space_bdy(i,k,jj,P_YEB) = 0. ! zeroout
3494        END DO
3495        END DO
3496        END DO
3497     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3498        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3499        DO k = kds , kde
3500        DO i = MAX(ids,its) , MIN(ide-1,ite)
3501           jj = jde - j
3502           space_bdy(i,k,jj,P_YEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3503 !         space_bdy(i,k,jj,P_YEB) = 0. ! zeroout
3504        END DO
3505        END DO
3506        END DO
3507     ELSE
3508        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3509        DO k = kds , kde - 1
3510        DO i = MAX(ids,its) , MIN(ide-1,ite)
3511           jj = jde - j
3512           space_bdy(i,k,jj,P_YEB) = ( data3dnew(i,k,j) - data3dold(i,k,j) ) / time_diff
3513 !         space_bdy(i,k,jj,P_YEB) = 0. ! zeroout
3514        END DO
3515        END DO
3516        END DO
3517     END IF
3518     
3519  END SUBROUTINE stuff_bdytend_old
3521  SUBROUTINE stuff_bdy_ijk ( data3d , space_bdy_xs, space_bdy_xe, &
3522                              space_bdy_ys, space_bdy_ye, &
3523                              char_stagger , spec_bdy_width, &
3524                              ids, ide, jds, jde, kds, kde , &
3525                              ims, ime, jms, jme, kms, kme , &
3526                              its, ite, jts, jte, kts, kte )
3528  !  This routine puts the data in the 3d arrays into the proper locations
3529  !  for the lateral boundary arrays.
3531     USE module_state_description
3532     
3533     IMPLICIT NONE
3535     INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde
3536     INTEGER , INTENT(IN) :: ims, ime, jms, jme, kms, kme
3537     INTEGER , INTENT(IN) :: its, ite, jts, jte, kts, kte
3538     INTEGER , INTENT(IN) :: spec_bdy_width
3539     REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(IN) :: data3d
3540 !    REAL , DIMENSION(:,:,:,:) , INTENT(OUT) :: space_bdy
3541 !    REAL , DIMENSION(ijds:ijde,kds:kde,spec_bdy_width,4,1) , INTENT(OUT) :: space_bdy
3542     REAL , DIMENSION(jms:jme,kds:kde,spec_bdy_width) , INTENT(OUT) :: space_bdy_xs, space_bdy_xe
3543     REAL , DIMENSION(ims:ime,kds:kde,spec_bdy_width) , INTENT(OUT) :: space_bdy_ys, space_bdy_ye
3544     CHARACTER (LEN=1) , INTENT(IN) :: char_stagger
3546     INTEGER :: i , ii , j , jj , k
3548     !  There are four lateral boundary locations that are stored.
3550     !  X start boundary
3552     IF ( char_stagger .EQ. 'W' ) THEN
3553        DO k = kds , kde
3554        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3555        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3556           space_bdy_xs(j,k,i) = data3d(i,j,k)
3557        END DO
3558        END DO
3559        END DO
3560     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3561        DO k = kds , kde
3562        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3563        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3564           space_bdy_xs(j,k,i) = data3d(i,j,k)
3565        END DO
3566        END DO
3567        END DO
3568     ELSE IF ( char_stagger .EQ. 'V' ) THEN
3569        DO k = kds , kde - 1
3570        DO j = MAX(jds,jts) , MIN(jde,jte)
3571        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3572           space_bdy_xs(j,k,i) = data3d(i,j,k)
3573        END DO
3574        END DO
3575        END DO
3576     ELSE
3577        DO k = kds , kde - 1
3578        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3579        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3580           space_bdy_xs(j,k,i) = data3d(i,j,k)
3581        END DO
3582        END DO
3583        END DO
3584     END IF
3586     !  X end boundary
3588     IF      ( char_stagger .EQ. 'U' ) THEN
3589        DO k = kds , kde - 1
3590        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3591        DO i = MIN(ide,ite) , MAX(ide - spec_bdy_width + 1,its) , -1
3592           ii = ide - i + 1
3593           space_bdy_xe(j,k,ii) = data3d(i,j,k)
3594        END DO
3595        END DO
3596        END DO
3597     ELSE IF ( char_stagger .EQ. 'V' ) THEN
3598        DO k = kds , kde - 1
3599        DO j = MAX(jds,jts) , MIN(jde,jte)
3600        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3601           ii = ide - i
3602           space_bdy_xe(j,k,ii) = data3d(i,j,k)
3603        END DO
3604        END DO
3605        END DO
3606     ELSE IF ( char_stagger .EQ. 'W' ) THEN
3607        DO k = kds , kde
3608        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3609        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3610           ii = ide - i
3611           space_bdy_xe(j,k,ii) = data3d(i,j,k)
3612        END DO
3613        END DO
3614        END DO
3615     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3616        DO k = kds , kde
3617        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3618        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3619           ii = ide - i
3620           space_bdy_xe(j,k,ii) = data3d(i,j,k)
3621        END DO
3622        END DO
3623        END DO
3624     ELSE
3625        DO k = kds , kde - 1
3626        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3627        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3628           ii = ide - i
3629           space_bdy_xe(j,k,ii) = data3d(i,j,k)
3630        END DO
3631        END DO
3632        END DO
3633     END IF
3635     !  Y start boundary
3637     IF ( char_stagger .EQ. 'W' ) THEN
3638        DO k = kds , kde
3639        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3640        DO i = MAX(ids,its) , MIN(ide-1,ite)
3641           space_bdy_ys(i,k,j) = data3d(i,j,k)
3642        END DO
3643        END DO
3644        END DO
3645     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3646        DO k = kds , kde
3647        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3648        DO i = MAX(ids,its) , MIN(ide-1,ite)
3649           space_bdy_ys(i,k,j) = data3d(i,j,k)
3650        END DO
3651        END DO
3652        END DO
3653     ELSE IF ( char_stagger .EQ. 'U' ) THEN
3654        DO k = kds , kde - 1
3655        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3656        DO i = MAX(ids,its) , MIN(ide,ite)
3657           space_bdy_ys(i,k,j) = data3d(i,j,k)
3658        END DO
3659        END DO
3660        END DO
3661     ELSE
3662        DO k = kds , kde - 1
3663        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3664        DO i = MAX(ids,its) , MIN(ide-1,ite)
3665           space_bdy_ys(i,k,j) = data3d(i,j,k)
3666        END DO
3667        END DO
3668        END DO
3669     END IF
3671     !  Y end boundary
3673     IF      ( char_stagger .EQ. 'V' ) THEN
3674        DO k = kds , kde - 1
3675        DO j = MIN(jde,jte) , MAX(jde - spec_bdy_width + 1,jts) , -1
3676        DO i = MAX(ids,its) , MIN(ide-1,ite)
3677           jj = jde - j + 1
3678           space_bdy_ye(i,k,jj) = data3d(i,j,k)
3679        END DO
3680        END DO
3681        END DO
3682     ELSE IF ( char_stagger .EQ. 'U' ) THEN
3683        DO k = kds , kde - 1
3684        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3685        DO i = MAX(ids,its) , MIN(ide,ite)
3686           jj = jde - j
3687           space_bdy_ye(i,k,jj) = data3d(i,j,k)
3688        END DO
3689        END DO
3690        END DO
3691     ELSE IF ( char_stagger .EQ. 'W' ) THEN
3692        DO k = kds , kde
3693        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3694        DO i = MAX(ids,its) , MIN(ide-1,ite)
3695           jj = jde - j
3696           space_bdy_ye(i,k,jj) = data3d(i,j,k)
3697        END DO
3698        END DO
3699        END DO
3700     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3701        DO k = kds , kde
3702        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3703        DO i = MAX(ids,its) , MIN(ide-1,ite)
3704           jj = jde - j
3705           space_bdy_ye(i,k,jj) = data3d(i,j,k)
3706        END DO
3707        END DO
3708        END DO
3709     ELSE
3710        DO k = kds , kde - 1
3711        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3712        DO i = MAX(ids,its) , MIN(ide-1,ite)
3713           jj = jde - j
3714           space_bdy_ye(i,k,jj) = data3d(i,j,k)
3715 !        if (K .eq. 54 .and. I .eq. 369) then
3716 ! write(0,*) 'N bound i,k,jj,P_YEB,data3d,space_bdy: ', i,k,jj,P_YEB,data3d(I,j,k),space_bdy(i,k,jj,P_YEB,1)
3717 ! endif
3719        END DO
3720        END DO
3721        END DO
3722     END IF
3723     
3724  END SUBROUTINE stuff_bdy_ijk
3726  SUBROUTINE stuff_bdytend_ijk ( data3dnew , data3dold , time_diff , &
3727                              space_bdy_xs, space_bdy_xe, space_bdy_ys, space_bdy_ye, &
3728                              char_stagger , &
3729                              spec_bdy_width , &
3730                              ids, ide, jds, jde, kds, kde , &
3731                              ims, ime, jms, jme, kms, kme , &
3732                              its, ite, jts, jte, kts, kte )
3734  !  This routine puts the tendency data into the proper locations
3735  !  for the lateral boundary arrays.
3737     USE module_state_description
3738     
3739     IMPLICIT NONE
3741     INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde
3742     INTEGER , INTENT(IN) :: ims, ime, jms, jme, kms, kme
3743     INTEGER , INTENT(IN) :: its, ite, jts, jte, kts, kte
3744     INTEGER , INTENT(IN) :: spec_bdy_width
3745 !    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: data3dnew , data3dold
3746     REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(IN) :: data3dnew , data3dold
3747     REAL , DIMENSION(jms:jme,kds:kde,spec_bdy_width) , INTENT(OUT) :: space_bdy_xs, space_bdy_xe
3748     REAL , DIMENSION(ims:ime,kds:kde,spec_bdy_width) , INTENT(OUT) :: space_bdy_ys, space_bdy_ye
3750     CHARACTER (LEN=1) , INTENT(IN) :: char_stagger
3751     REAL , INTENT(IN) :: time_diff ! seconds
3753     INTEGER :: i , ii , j , jj , k
3755     !  There are four lateral boundary locations that are stored.
3757     !  X start boundary
3759     IF ( char_stagger .EQ. 'W' ) THEN
3760        DO k = kds , kde
3761        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3762        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3763           space_bdy_xs(j,k,i) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3764        END DO
3765        END DO
3766        END DO
3767     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3768        DO k = kds , kde
3769        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3770        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3771           space_bdy_xs(j,k,i) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3772        END DO
3773        END DO
3774        END DO
3775     ELSE IF ( char_stagger .EQ. 'V' ) THEN
3776        DO k = kds , kde - 1
3777        DO j = MAX(jds,jts) , MIN(jde,jte)
3778        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3779           space_bdy_xs(j,k,i) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3780        END DO
3781        END DO
3782        END DO
3783     ELSE
3784        DO k = kds , kde - 1
3785        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3786        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
3787           space_bdy_xs(j,k,i) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3788        END DO
3789        END DO
3790        END DO
3791     END IF
3793     !  X end boundary
3795     IF      ( char_stagger .EQ. 'U' ) THEN
3796        DO k = kds , kde - 1
3797        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3798        DO i = MIN(ide,ite) , MAX(ide - spec_bdy_width + 1,its) , -1
3799           ii = ide - i + 1
3800           space_bdy_xe(j,k,ii) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3801        END DO
3802        END DO
3803        END DO
3804     ELSE IF ( char_stagger .EQ. 'V' ) THEN
3805        DO k = kds , kde - 1
3806        DO j = MAX(jds,jts) , MIN(jde,jte)
3807        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3808           ii = ide - i
3809           space_bdy_xe(j,k,ii) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3810        END DO
3811        END DO
3812        END DO
3813     ELSE IF ( char_stagger .EQ. 'W' ) THEN
3814        DO k = kds , kde
3815        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3816        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3817           ii = ide - i
3818           space_bdy_xe(j,k,ii) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3819        END DO
3820        END DO
3821        END DO
3822     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3823        DO k = kds , kde
3824        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3825        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3826           ii = ide - i
3827           space_bdy_xe(j,k,ii) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3828        END DO
3829        END DO
3830        END DO
3831     ELSE
3832        DO k = kds , kde - 1
3833        DO j = MAX(jds,jts) , MIN(jde-1,jte)
3834        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
3835           ii = ide - i
3836           space_bdy_xe(j,k,ii) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3837        END DO
3838        END DO
3839        END DO
3840     END IF
3842     !  Y start boundary
3844     IF ( char_stagger .EQ. 'W' ) THEN
3845        DO k = kds , kde
3846        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3847        DO i = MAX(ids,its) , MIN(ide-1,ite)
3848           space_bdy_ys(i,k,j) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3849        END DO
3850        END DO
3851        END DO
3852     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3853        DO k = kds , kde
3854        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3855        DO i = MAX(ids,its) , MIN(ide-1,ite)
3856           space_bdy_ys(i,k,j) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3857        END DO
3858        END DO
3859        END DO
3860     ELSE IF ( char_stagger .EQ. 'U' ) THEN
3861        DO k = kds , kde - 1
3862        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3863        DO i = MAX(ids,its) , MIN(ide,ite)
3864           space_bdy_ys(i,k,j) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3865        END DO
3866        END DO
3867        END DO
3868     ELSE
3869        DO k = kds , kde - 1
3870        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
3871        DO i = MAX(ids,its) , MIN(ide-1,ite)
3872           space_bdy_ys(i,k,j) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3873        END DO
3874        END DO
3875        END DO
3876     END IF
3878     !  Y end boundary
3880     IF      ( char_stagger .EQ. 'V' ) THEN
3881        DO k = kds , kde - 1
3882        DO j = MIN(jde,jte) , MAX(jde - spec_bdy_width + 1,jts) , -1
3883        DO i = MAX(ids,its) , MIN(ide-1,ite)
3884           jj = jde - j + 1
3885           space_bdy_ye(i,k,jj) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3886        END DO
3887        END DO
3888        END DO
3889     ELSE IF ( char_stagger .EQ. 'U' ) THEN
3890        DO k = kds , kde - 1
3891        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3892        DO i = MAX(ids,its) , MIN(ide,ite)
3893           jj = jde - j
3894           space_bdy_ye(i,k,jj) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3895        END DO
3896        END DO
3897        END DO
3898     ELSE IF ( char_stagger .EQ. 'W' ) THEN
3899        DO k = kds , kde
3900        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3901        DO i = MAX(ids,its) , MIN(ide-1,ite)
3902           jj = jde - j
3903           space_bdy_ye(i,k,jj) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3904        END DO
3905        END DO
3906        END DO
3907     ELSE IF ( char_stagger .EQ. 'M' ) THEN
3908        DO k = kds , kde
3909        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3910        DO i = MAX(ids,its) , MIN(ide-1,ite)
3911           jj = jde - j
3912           space_bdy_ye(i,k,jj) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3913        END DO
3914        END DO
3915        END DO
3916     ELSE
3917        DO k = kds , kde - 1
3918        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
3919        DO i = MAX(ids,its) , MIN(ide-1,ite)
3920           jj = jde - j
3921           space_bdy_ye(i,k,jj) = ( data3dnew(i,j,k) - data3dold(i,j,k) ) / time_diff
3922 !        if (K .eq. 54 .and. I .eq. 369) then
3923 ! write(0,*) 'N bound i,k,jj,data3dnew,data3dold: ', i,k,jj,data3dnew(I,j,k),data3dold(i,j,k)
3924 ! endif
3925        END DO
3926        END DO
3927        END DO
3928     END IF
3929     
3930  END SUBROUTINE stuff_bdytend_ijk
3932 #if ( WRFPLUS == 1 )
3933 !------------------------------------------------------------------------
3935    SUBROUTINE couple_bdy ( field,      &
3936                            variable_in, config_flags, & 
3937                            spec_zone,  &
3938                            mu, msf,    &
3939                            ids,ide, jds,jde, kds,kde,  & ! domain dims
3940                            ims,ime, jms,jme, kms,kme,  & ! memory dims
3941                            ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
3942                            its,ite, jts,jte, kts,kte )
3944 !  This subroutine adds the tendencies in the boundary specified region.
3945 !  spec_zone is the width of the outer specified b.c.s that are set here.
3946 !  (JD August 2000)
3948       IMPLICIT NONE
3950       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
3951       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
3952       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
3953       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
3954       INTEGER,      INTENT(IN   )    :: spec_zone
3955       CHARACTER,    INTENT(IN   )    :: variable_in
3956       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: mu
3957       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: msf
3958       REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field
3959       TYPE( grid_config_rec_type ) config_flags
3961       CHARACTER  :: variable
3962       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf
3963       INTEGER    :: b_dist, b_limit
3964       LOGICAL    :: periodic_x
3966       periodic_x = config_flags%periodic_x
3968       variable = variable_in
3970       IF (variable == 'U') variable = 'u'
3971       IF (variable == 'V') variable = 'v'
3972       IF (variable == 'T') variable = 't'
3973       IF (variable == 'H') variable = 'h'
3974       IF (variable == 'W') variable = 'w'
3976       ibs = ids
3977       ibe = ide-1
3978       itf = min(ite,ide-1)
3979       jbs = jds
3980       jbe = jde-1
3981       jtf = min(jte,jde-1)
3982       ktf = kde-1
3983       IF (variable == 'u') ibe = ide
3984       IF (variable == 'u') itf = min(ite,ide)
3985       IF (variable == 'v') jbe = jde
3986       IF (variable == 'v') jtf = min(jte,jde)
3987       IF (variable == 'h') ktf = kte
3988       IF (variable == 'w') ktf = kte
3990       IF (jts - jbs .lt. spec_zone) THEN
3991 ! Y-start boundary
3992         DO j = jts, min(jtf,jbs+spec_zone-1)
3993           b_dist = j - jbs
3994           b_limit = b_dist
3995           IF(periodic_x)b_limit = 0
3996           DO k = kts, ktf
3997             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
3998               if (variable == 't' .or. variable == 'h') then 
3999                  field(i,k,j) = field(i,k,j)*mu(i,j)
4000               else 
4001                  field(i,k,j) = field(i,k,j)*mu(i,j)/msf(i,j)
4002               end if
4003             ENDDO
4004           ENDDO
4005         ENDDO
4006       ENDIF 
4007       IF (jbe - jtf .lt. spec_zone) THEN 
4008 ! Y-end boundary 
4009         DO j = max(jts,jbe-spec_zone+1), jtf 
4010           b_dist = jbe - j 
4011           b_limit = b_dist
4012           IF(periodic_x)b_limit = 0
4013           DO k = kts, ktf 
4014             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
4015               if (variable == 't' .or. variable == 'h') then 
4016                  field(i,k,j) = field(i,k,j)*mu(i,j)
4017               else 
4018                  field(i,k,j) = field(i,k,j)*mu(i,j)/msf(i,j)
4019               end if
4020             ENDDO
4021           ENDDO
4022         ENDDO
4023       ENDIF 
4025     IF(.NOT.periodic_x)THEN
4026       IF (its - ibs .lt. spec_zone) THEN
4027 ! X-start boundary
4028         DO i = its, min(itf,ibs+spec_zone-1)
4029           b_dist = i - ibs
4030           DO k = kts, ktf
4031             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
4032               if (variable == 't' .or. variable == 'h') then 
4033                  field(i,k,j) = field(i,k,j)*mu(i,j)
4034               else 
4035                  field(i,k,j) = field(i,k,j)*mu(i,j)/msf(i,j)
4036               end if
4037             ENDDO
4038           ENDDO
4039         ENDDO
4040       ENDIF 
4042       IF (ibe - itf .lt. spec_zone) THEN
4043 ! X-end boundary
4044         DO i = max(its,ibe-spec_zone+1), itf
4045           b_dist = ibe - i
4046           DO k = kts, ktf
4047             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
4048               if (variable == 't' .or. variable == 'h') then 
4049                  field(i,k,j) = field(i,k,j)*mu(i,j)
4050               else 
4051                  field(i,k,j) = field(i,k,j)*mu(i,j)/msf(i,j)
4052               end if
4053             ENDDO
4054           ENDDO
4055         ENDDO
4056       ENDIF 
4057     ENDIF
4059    END SUBROUTINE couple_bdy
4061 !------------------------------------------------------------------------
4063    SUBROUTINE uncouple_bdy ( field,      &
4064                              variable_in, config_flags, & 
4065                              spec_zone,  &
4066                              mu, msf,    &
4067                              ids,ide, jds,jde, kds,kde,  & ! domain dims
4068                              ims,ime, jms,jme, kms,kme,  & ! memory dims
4069                              ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
4070                              its,ite, jts,jte, kts,kte )
4072 !  This subroutine adds the tendencies in the boundary specified region.
4073 !  spec_zone is the width of the outer specified b.c.s that are set here.
4074 !  (JD August 2000)
4076       IMPLICIT NONE
4078       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
4079       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
4080       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
4081       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
4082       INTEGER,      INTENT(IN   )    :: spec_zone
4083       CHARACTER,    INTENT(IN   )    :: variable_in
4084       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: mu
4085       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: msf
4086       REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field
4087       TYPE( grid_config_rec_type ) config_flags
4089       CHARACTER  :: variable
4090       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf
4091       INTEGER    :: b_dist, b_limit
4092       LOGICAL    :: periodic_x
4094       periodic_x = config_flags%periodic_x
4096       variable = variable_in
4098       IF (variable == 'U') variable = 'u'
4099       IF (variable == 'V') variable = 'v'
4100       IF (variable == 'T') variable = 't'
4101       IF (variable == 'H') variable = 'h'
4102       IF (variable == 'W') variable = 'w'
4104       ibs = ids
4105       ibe = ide-1
4106       itf = min(ite,ide-1)
4107       jbs = jds
4108       jbe = jde-1
4109       jtf = min(jte,jde-1)
4110       ktf = kde-1
4111       IF (variable == 'u') ibe = ide
4112       IF (variable == 'u') itf = min(ite,ide)
4113       IF (variable == 'v') jbe = jde
4114       IF (variable == 'v') jtf = min(jte,jde)
4115       IF (variable == 'h') ktf = kte
4116       IF (variable == 'w') ktf = kte
4118       IF (jts - jbs .lt. spec_zone) THEN
4119 ! Y-start boundary
4120         DO j = jts, min(jtf,jbs+spec_zone-1)
4121           b_dist = j - jbs
4122           b_limit = b_dist
4123           IF(periodic_x)b_limit = 0
4124           DO k = kts, ktf
4125             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
4126               if (variable == 't' .or. variable == 'h') then 
4127                  field(i,k,j) = field(i,k,j)/mu(i,j)
4128               else 
4129                  field(i,k,j) = field(i,k,j)/mu(i,j)*msf(i,j)
4130               end if
4131             ENDDO
4132           ENDDO
4133         ENDDO
4134       ENDIF 
4135       IF (jbe - jtf .lt. spec_zone) THEN 
4136 ! Y-end boundary 
4137         DO j = max(jts,jbe-spec_zone+1), jtf 
4138           b_dist = jbe - j 
4139           b_limit = b_dist
4140           IF(periodic_x)b_limit = 0
4141           DO k = kts, ktf 
4142             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
4143               if (variable == 't' .or. variable == 'h') then 
4144                  field(i,k,j) = field(i,k,j)/mu(i,j)
4145               else 
4146                  field(i,k,j) = field(i,k,j)/mu(i,j)*msf(i,j)
4147               end if
4148             ENDDO
4149           ENDDO
4150         ENDDO
4151       ENDIF 
4153     IF(.NOT.periodic_x)THEN
4154       IF (its - ibs .lt. spec_zone) THEN
4155 ! X-start boundary
4156         DO i = its, min(itf,ibs+spec_zone-1)
4157           b_dist = i - ibs
4158           DO k = kts, ktf
4159             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
4160               if (variable == 't' .or. variable == 'h') then 
4161                  field(i,k,j) = field(i,k,j)/mu(i,j)
4162               else 
4163                  field(i,k,j) = field(i,k,j)/mu(i,j)*msf(i,j)
4164               end if
4165             ENDDO
4166           ENDDO
4167         ENDDO
4168       ENDIF 
4170       IF (ibe - itf .lt. spec_zone) THEN
4171 ! X-end boundary
4172         DO i = max(its,ibe-spec_zone+1), itf
4173           b_dist = ibe - i
4174           DO k = kts, ktf
4175             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
4176               if (variable == 't' .or. variable == 'h') then 
4177                  field(i,k,j) = field(i,k,j)/mu(i,j)
4178               else 
4179                  field(i,k,j) = field(i,k,j)/mu(i,j)*msf(i,j)
4180               end if
4181             ENDDO
4182           ENDDO
4183         ENDDO
4184       ENDIF 
4185     ENDIF
4187    END SUBROUTINE uncouple_bdy
4189 !---------------------------------------
4190  SUBROUTINE bdy_fields_halo (data3du, data3dv, data3dt, data3dph, data3dmu, &
4191                              data3dm, dir, xy, spec_bdy_width,              &
4192                              u_bxs,  u_bxe,  u_bys,  u_bye,                 &
4193                              v_bxs,  v_bxe,  v_bys,  v_bye,                 &
4194                              t_bxs,  t_bxe,  t_bys,  t_bye,                 &
4195                              ph_bxs, ph_bxe, ph_bys, ph_bye,                &
4196                              mu_bxs, mu_bxe, mu_bys, mu_bye,                &
4197                              moist_bxs, moist_bxe, moist_bys, moist_bye,    &
4198                              ids, ide, jds, jde, kds, kde ,                 &
4199                              ims, ime, jms, jme, kms, kme ,                 &
4200                              its, ite, jts, jte, kts, kte )
4202     USE module_state_description
4203        
4204     IMPLICIT NONE
4206     INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde
4207     INTEGER , INTENT(IN) :: ims, ime, jms, jme, kms, kme
4208     INTEGER , INTENT(IN) :: its, ite, jts, jte, kts, kte
4209     INTEGER , INTENT(IN) :: spec_bdy_width
4210     INTEGER , INTENT(IN) :: dir ! 0----pack ; 1----unpack
4211     INTEGER , INTENT(IN) :: xy  ! 0----X ; 1----Y
4212     REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: data3du , data3dv, data3dt, &
4213                                                                  data3dph, data3dm
4214     REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: data3dmu
4215     REAL , DIMENSION(jms:jme,kds:kde,spec_bdy_width) , INTENT(INOUT) :: u_bxs, u_bxe, v_bxs, v_bxe, &
4216                                                                         t_bxs, t_bxe, ph_bxs, ph_bxe, &
4217                                                                         moist_bxs, moist_bxe
4218     REAL , DIMENSION(jms:jme,1:1,spec_bdy_width) , INTENT(INOUT) :: mu_bxs, mu_bxe
4219     REAL , DIMENSION(ims:ime,kds:kde,spec_bdy_width) , INTENT(INOUT) :: u_bys, u_bye, v_bys, v_bye, &
4220                                                                         t_bys, t_bye, ph_bys, ph_bye, &
4221                                                                         moist_bys, moist_bye
4222     REAL , DIMENSION(ims:ime,1:1,spec_bdy_width) , INTENT(INOUT) :: mu_bys, mu_bye
4224    CALL bdy_fields_pack ( data3du, u_bxs, u_bxe, u_bys, u_bye,        &
4225                                     'U' , dir, xy, spec_bdy_width,    &
4226                                     ids, ide, jds, jde, kds, kde,     &
4227                                     ims, ime, jms, jme, kms, kme,     &
4228                                     its, ite, jts, jte, kts, kte      )
4230    CALL bdy_fields_pack ( data3dv, v_bxs, v_bxe, v_bys, v_bye,        &
4231                                     'V' , dir, xy, spec_bdy_width,    &
4232                                     ids, ide, jds, jde, kds, kde,     &
4233                                     ims, ime, jms, jme, kms, kme,     &
4234                                     its, ite, jts, jte, kts, kte      )
4236    CALL bdy_fields_pack ( data3dt , t_bxs, t_bxe, t_bys, t_bye,       &
4237                                     'T' , dir, xy, spec_bdy_width,    &
4238                                     ids, ide, jds, jde, kds, kde,     &
4239                                     ims, ime, jms, jme, kms, kme,     &
4240                                     its, ite, jts, jte, kts, kte      )
4241    
4242    CALL bdy_fields_pack ( data3dph , ph_bxs, ph_bxe, ph_bys, ph_bye,  &
4243                                     'W' , dir, xy, spec_bdy_width,    &
4244                                     ids, ide, jds, jde, kds, kde,     &
4245                                     ims, ime, jms, jme, kms, kme,     &
4246                                     its, ite, jts, jte, kts, kte      )
4248    CALL bdy_fields_pack ( data3dmu , mu_bxs, mu_bxe, mu_bys, mu_bye,  &
4249                                     'M' , dir, xy, spec_bdy_width  ,  &
4250                                     ids, ide, jds, jde, 1  , 1  ,     &
4251                                     ims, ime, jms, jme, 1  , 1  ,     &
4252                                     its, ite, jts, jte, 1  , 1        )
4254    CALL bdy_fields_pack ( data3dm , moist_bxs, moist_bxe, moist_bys, moist_bye,       &
4255                                     'T' , dir, xy, spec_bdy_width,    &
4256                                     ids, ide, jds, jde, kds, kde,     &
4257                                     ims, ime, jms, jme, kms, kme,     &
4258                                     its, ite, jts, jte, kts, kte      )
4261  END SUBROUTINE bdy_fields_halo
4263  SUBROUTINE bdy_fields_pack ( data3d , space_bdy_xs, space_bdy_xe, space_bdy_ys, space_bdy_ye, &
4264                              char_stagger , dir , xy ,&
4265                              spec_bdy_width , &         
4266                              ids, ide, jds, jde, kds, kde , &
4267                              ims, ime, jms, jme, kms, kme , &
4268                              its, ite, jts, jte, kts, kte )
4269                           
4270 !-------------------------------------------------------------------------
4271 !  Calculate the first guess at the end of the time window
4272 !  Author: Xin Zhang, 10/7/2010
4273 !-------------------------------------------------------------------------
4274       
4275     IMPLICIT NONE
4276       
4277     INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde
4278     INTEGER , INTENT(IN) :: ims, ime, jms, jme, kms, kme
4279     INTEGER , INTENT(IN) :: its, ite, jts, jte, kts, kte
4280     INTEGER , INTENT(IN) :: spec_bdy_width
4281     REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: data3d
4282     REAL , DIMENSION(jms:jme,kds:kde,spec_bdy_width) , INTENT(INOUT) :: space_bdy_xs, space_bdy_xe
4283     REAL , DIMENSION(ims:ime,kds:kde,spec_bdy_width) , INTENT(INOUT) :: space_bdy_ys, space_bdy_ye
4284     CHARACTER (LEN=1) , INTENT(IN) :: char_stagger
4285     INTEGER , INTENT(IN) :: dir ! 0----pack ; 1----unpack
4286     INTEGER , INTENT(IN) :: xy  ! 0----X ; 1----Y
4287       
4288     INTEGER :: i , ii , j , jj , k
4289          
4290     !  There are four lateral boundary locations that are stored.
4291          
4292 IF (dir == 0 ) THEN  ! ----Pack
4293    IF ( xy == 0 ) THEN  ! ----X 
4295     !  X start boundary
4296       
4297     IF ( char_stagger .EQ. 'W' ) THEN
4298        DO j = jms, jme
4299        DO k = kds , kde
4300        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
4301           data3d(i,k,j) = space_bdy_xs(j,k,i)
4302           space_bdy_xs(j,k,i) = 0.0
4303        END DO
4304        END DO
4305        END DO
4306     ELSE IF ( char_stagger .EQ. 'M' ) THEN
4307        DO j = jms, jme
4308        DO k = kds , kde
4309        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
4310           data3d(i,k,j) = space_bdy_xs(j,k,i)
4311           space_bdy_xs(j,k,i) = 0.0
4312        END DO
4313        END DO
4314        END DO
4315     ELSE IF ( char_stagger .EQ. 'V' ) THEN
4316        DO j = jms, jme
4317        DO k = kds , kde - 1
4318        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
4319           data3d(i,k,j) = space_bdy_xs(j,k,i)
4320           space_bdy_xs(j,k,i) = 0.0
4321        END DO
4322        END DO
4323        END DO
4324     ELSE
4325        DO j = jms, jme
4326        DO k = kds , kde - 1
4327        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
4328           data3d(i,k,j) = space_bdy_xs(j,k,i)
4329           space_bdy_xs(j,k,i) = 0.0
4330        END DO
4331        END DO
4332        END DO
4333     END IF
4335     !  X end boundary
4337     IF      ( char_stagger .EQ. 'U' ) THEN
4338        DO j = jms, jme
4339        DO k = kds , kde - 1
4340        DO i = MIN(ide,ite) , MAX(ide - spec_bdy_width + 1,its) , -1
4341           ii = ide - i + 1
4342           data3d(i,k,j) = space_bdy_xe(j,k,ii)
4343           space_bdy_xe(j,k,ii) = 0.0
4344        END DO
4345        END DO
4346        END DO
4347     ELSE IF ( char_stagger .EQ. 'V' ) THEN
4348        DO j = jms, jme
4349        DO k = kds , kde - 1
4350        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
4351           ii = ide - i
4352           data3d(i,k,j) = space_bdy_xe(j,k,ii)
4353           space_bdy_xe(j,k,ii) = 0.0
4354        END DO
4355        END DO
4356        END DO
4357     ELSE IF ( char_stagger .EQ. 'W' ) THEN
4358        DO j = jms, jme
4359        DO k = kds , kde
4360        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
4361           ii = ide - i
4362           data3d(i,k,j) = space_bdy_xe(j,k,ii)
4363           space_bdy_xe(j,k,ii) = 0.0
4364        END DO
4365        END DO
4366        END DO
4367     ELSE IF ( char_stagger .EQ. 'M' ) THEN
4368        DO j = jms, jme
4369        DO k = kds , kde
4370        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
4371           ii = ide - i
4372           data3d(i,k,j) = space_bdy_xe(j,k,ii)
4373           space_bdy_xe(j,k,ii) = 0.0
4374        END DO
4375        END DO
4376        END DO
4377     ELSE
4378        DO j = jms, jme
4379        DO k = kds , kde - 1
4380        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
4381           ii = ide - i
4382           data3d(i,k,j) = space_bdy_xe(j,k,ii)
4383           space_bdy_xe(j,k,ii) = 0.0
4384        END DO
4385        END DO
4386        END DO
4387     END IF
4388     
4389 ELSE    !    Y 
4390     !  Y start boundary
4391     
4392     IF ( char_stagger .EQ. 'W' ) THEN
4393        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
4394        DO k = kds , kde
4395        DO i = ims, ime
4396           data3d(i,k,j) = space_bdy_ys(i,k,j)
4397           space_bdy_ys(i,k,j) = 0.0
4398        END DO
4399        END DO
4400        END DO 
4401     ELSE IF ( char_stagger .EQ. 'M' ) THEN
4402        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
4403        DO k = kds , kde
4404        DO i = ims, ime
4405           data3d(i,k,j) = space_bdy_ys(i,k,j)
4406           space_bdy_ys(i,k,j) = 0.0
4407        END DO
4408        END DO
4409        END DO 
4410     ELSE IF ( char_stagger .EQ. 'U' ) THEN
4411        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
4412        DO k = kds , kde - 1
4413        DO i = ims, ime
4414           data3d(i,k,j) = space_bdy_ys(i,k,j)
4415           space_bdy_ys(i,k,j) = 0.0
4416        END DO
4417        END DO
4418        END DO
4419     ELSE
4420        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
4421        DO k = kds , kde - 1
4422        DO i = ims, ime
4423           data3d(i,k,j) = space_bdy_ys(i,k,j)
4424           space_bdy_ys(i,k,j) = 0.0
4425        END DO
4426        END DO
4427        END DO
4428     END IF
4430     !  Y end boundary
4432     IF      ( char_stagger .EQ. 'V' ) THEN
4433        DO j = MIN(jde,jte) , MAX(jde - spec_bdy_width + 1,jts) , -1
4434        DO k = kds , kde - 1
4435        DO i = ims, ime
4436           jj = jde - j + 1
4437           data3d(i,k,j) = space_bdy_ye(i,k,jj)
4438           space_bdy_ye(i,k,jj) = 0.0
4439        END DO
4440        END DO
4441        END DO
4442     ELSE IF ( char_stagger .EQ. 'U' ) THEN
4443        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
4444        DO k = kds , kde - 1
4445        DO i = ims, ime
4446           jj = jde - j
4447           data3d(i,k,j) = space_bdy_ye(i,k,jj)
4448           space_bdy_ye(i,k,jj) = 0.0
4449        END DO
4450        END DO
4451        END DO
4452     ELSE IF ( char_stagger .EQ. 'W' ) THEN
4453        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
4454        DO k = kds , kde
4455        DO i = ims, ime
4456           jj = jde - j
4457           data3d(i,k,j) = space_bdy_ye(i,k,jj)
4458           space_bdy_ye(i,k,jj) = 0.0
4459        END DO
4460        END DO
4461        END DO
4462     ELSE IF ( char_stagger .EQ. 'M' ) THEN
4463        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
4464        DO k = kds , kde
4465        DO i = ims, ime
4466           jj = jde - j
4467           data3d(i,k,j) = space_bdy_ye(i,k,jj)
4468           space_bdy_ye(i,k,jj) = 0.0
4469        END DO
4470        END DO
4471        END DO
4472     ELSE
4473        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
4474        DO k = kds , kde - 1
4475        DO i = ims, ime
4476           jj = jde - j
4477           data3d(i,k,j) = space_bdy_ye(i,k,jj)
4478           space_bdy_ye(i,k,jj) = 0.0
4479        END DO
4480        END DO
4481        END DO
4482     END IF
4484    END IF
4485 END IF
4487 IF ( dir == 1 ) THEN  ! ---- Unpack
4488    IF ( xy == 0 ) THEN !----- X
4490     !  X start boundary
4491       
4492     IF ( char_stagger .EQ. 'W' ) THEN
4493        DO j = MAX(jds,jts) , MIN(jde-1,jte)
4494        DO k = kds , kde
4495        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
4496           space_bdy_xs(j,k,i) = data3d(i,k,j)
4497        END DO
4498        END DO
4499        END DO
4500     ELSE IF ( char_stagger .EQ. 'M' ) THEN
4501        DO j = MAX(jds,jts) , MIN(jde-1,jte)
4502        DO k = kds , kde
4503        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
4504           space_bdy_xs(j,k,i) = data3d(i,k,j)
4505        END DO
4506        END DO
4507        END DO
4508     ELSE IF ( char_stagger .EQ. 'V' ) THEN
4509        DO j = MAX(jds,jts) , MIN(jde,jte)
4510        DO k = kds , kde - 1
4511        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
4512           space_bdy_xs(j,k,i) = data3d(i,k,j)
4513        END DO
4514        END DO
4515        END DO
4516     ELSE
4517        DO j = MAX(jds,jts) , MIN(jde-1,jte)
4518        DO k = kds , kde - 1
4519        DO i = MAX(ids,its) , MIN(ids + spec_bdy_width - 1,ite)
4520           space_bdy_xs(j,k,i) = data3d(i,k,j)
4521        END DO
4522        END DO
4523        END DO
4524     END IF
4526     !  X end boundary
4528     IF      ( char_stagger .EQ. 'U' ) THEN
4529        DO j = MAX(jds,jts) , MIN(jde-1,jte)
4530        DO k = kds , kde - 1
4531        DO i = MIN(ide,ite) , MAX(ide - spec_bdy_width + 1,its) , -1
4532           ii = ide - i + 1
4533           space_bdy_xe(j,k,ii) = data3d(i,k,j)
4534        END DO
4535        END DO
4536        END DO
4537     ELSE IF ( char_stagger .EQ. 'V' ) THEN
4538        DO j = MAX(jds,jts) , MIN(jde,jte)
4539        DO k = kds , kde - 1
4540        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
4541           ii = ide - i
4542           space_bdy_xe(j,k,ii) = data3d(i,k,j)
4543        END DO
4544        END DO
4545        END DO
4546     ELSE IF ( char_stagger .EQ. 'W' ) THEN
4547        DO j = MAX(jds,jts) , MIN(jde-1,jte)
4548        DO k = kds , kde
4549        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
4550           ii = ide - i
4551           space_bdy_xe(j,k,ii) = data3d(i,k,j)
4552        END DO
4553        END DO
4554        END DO
4555     ELSE IF ( char_stagger .EQ. 'M' ) THEN
4556        DO j = MAX(jds,jts) , MIN(jde-1,jte)
4557        DO k = kds , kde
4558        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
4559           ii = ide - i
4560           space_bdy_xe(j,k,ii) = data3d(i,k,j)
4561        END DO
4562        END DO
4563        END DO
4564     ELSE
4565        DO j = MAX(jds,jts) , MIN(jde-1,jte)
4566        DO k = kds , kde - 1
4567        DO i = MIN(ide - 1,ite) , MAX(ide - spec_bdy_width,its) , -1
4568           ii = ide - i
4569           space_bdy_xe(j,k,ii) = data3d(i,k,j)
4570        END DO
4571        END DO
4572        END DO
4573     END IF
4574     
4575 ELSE   !  Y  
4576     !  Y start boundary
4577     
4578     IF ( char_stagger .EQ. 'W' ) THEN
4579        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
4580        DO k = kds , kde
4581        DO i = MAX(ids,its) , MIN(ide-1,ite)
4582           space_bdy_ys(i,k,j) = data3d(i,k,j)
4583        END DO
4584        END DO
4585        END DO 
4586     ELSE IF ( char_stagger .EQ. 'M' ) THEN
4587        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
4588        DO k = kds , kde
4589        DO i = MAX(ids,its) , MIN(ide-1,ite)
4590           space_bdy_ys(i,k,j) = data3d(i,k,j)
4591        END DO
4592        END DO
4593        END DO 
4594     ELSE IF ( char_stagger .EQ. 'U' ) THEN
4595        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
4596        DO k = kds , kde - 1
4597        DO i = MAX(ids,its) , MIN(ide,ite)
4598           space_bdy_ys(i,k,j) = data3d(i,k,j)
4599        END DO
4600        END DO
4601        END DO
4602     ELSE
4603        DO j = MAX(jds,jts) , MIN(jds + spec_bdy_width - 1,jte)
4604        DO k = kds , kde - 1
4605        DO i = MAX(ids,its) , MIN(ide-1,ite)
4606           space_bdy_ys(i,k,j) = data3d(i,k,j)
4607        END DO
4608        END DO
4609        END DO
4610     END IF
4612     !  Y end boundary
4614     IF      ( char_stagger .EQ. 'V' ) THEN
4615        DO j = MIN(jde,jte) , MAX(jde - spec_bdy_width + 1,jts) , -1
4616        DO k = kds , kde - 1
4617        DO i = MAX(ids,its) , MIN(ide-1,ite)
4618           jj = jde - j + 1
4619           space_bdy_ye(i,k,jj) = data3d(i,k,j)
4620        END DO
4621        END DO
4622        END DO
4623     ELSE IF ( char_stagger .EQ. 'U' ) THEN
4624        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
4625        DO k = kds , kde - 1
4626        DO i = MAX(ids,its) , MIN(ide,ite)
4627           jj = jde - j
4628           space_bdy_ye(i,k,jj) = data3d(i,k,j)
4629        END DO
4630        END DO
4631        END DO
4632     ELSE IF ( char_stagger .EQ. 'W' ) THEN
4633        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
4634        DO k = kds , kde
4635        DO i = MAX(ids,its) , MIN(ide-1,ite)
4636           jj = jde - j
4637           space_bdy_ye(i,k,jj) = data3d(i,k,j)
4638        END DO
4639        END DO
4640        END DO
4641     ELSE IF ( char_stagger .EQ. 'M' ) THEN
4642        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
4643        DO k = kds , kde
4644        DO i = MAX(ids,its) , MIN(ide-1,ite)
4645           jj = jde - j
4646           space_bdy_ye(i,k,jj) = data3d(i,k,j)
4647        END DO
4648        END DO
4649        END DO
4650     ELSE
4651        DO j = MIN(jde-1,jte) , MAX(jde - spec_bdy_width,jts) , -1
4652        DO k = kds , kde - 1
4653        DO i = MAX(ids,its) , MIN(ide-1,ite)
4654           jj = jde - j
4655           space_bdy_ye(i,k,jj) = data3d(i,k,j)
4656        END DO
4657        END DO
4658        END DO
4659     END IF
4661    END IF
4662 END IF
4664  END SUBROUTINE bdy_fields_pack
4665 !----------------------------------------
4666 #endif
4668 END MODULE module_bc
4670 SUBROUTINE get_bdyzone_x ( bzx )
4671   USE module_bc
4672   IMPLICIT NONE
4673   INTEGER bzx
4674   bzx = bdyzone_x
4675 END SUBROUTINE get_bdyzone_x
4677 SUBROUTINE get_bdyzone_y ( bzy)
4678   USE module_bc
4679   IMPLICIT NONE
4680   INTEGER bzy
4681   bzy = bdyzone_y
4682 END SUBROUTINE get_bdyzone_y
4684 SUBROUTINE get_bdyzone ( bz)
4685   USE module_bc
4686   IMPLICIT NONE
4687   INTEGER bz
4688   bz = bdyzone
4689 END SUBROUTINE get_bdyzone