wrftladj/module_bc_ad.F

   1 !WRF+/AD:MODEL_LAYER:BOUNDARY
   2 !Created by Ning Pan, 2010-08
   3 !
   4
   5 MODULE a_module_bc
   6
   7    USE module_configure
   8    USE module_wrf_error
   9    USE module_model_constants
  10
  11    IMPLICIT NONE
  12
  13 !  set the bdyzone.  We are hardwiring this here and we'll
  14 !  decide later where it should be set and stored
  15
  16    INTEGER, PARAMETER            :: bdyzone = 4
  17    INTEGER, PARAMETER            :: bdyzone_x = bdyzone
  18    INTEGER, PARAMETER            :: bdyzone_y = bdyzone
  19
  20 CONTAINS
  21
  22 !------------------------------------------------------------------------
  23
  24    SUBROUTINE a_set_physical_bc2d( a_dat, variable_in,  &
  25                                  config_flags,           &
  26                                  ids,ide, jds,jde,   & ! domain dims
  27                                  ims,ime, jms,jme,   & ! memory dims
  28                                  ips,ipe, jps,jpe,   & ! patch  dims
  29                                  its,ite, jts,jte   )
  30
  31       IMPLICIT NONE
  32
  33       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde
  34       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme
  35       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe
  36       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte
  37       CHARACTER,    INTENT(IN   )    :: variable_in
  38
  39       CHARACTER                      :: variable
  40
  41       REAL,  DIMENSION( ims:ime , jms:jme ) :: a_dat
  42       TYPE( grid_config_rec_type ) config_flags
  43
  44       INTEGER  :: i, j, istag, jstag, itime
  45
  46       LOGICAL  :: debug, open_bc_copy
  47
  48       real :: a_aux
  49
  50 !------------
  51
  52       a_aux = 0.0
  53
  54       debug = .false.
  55
  56       open_bc_copy = .false.
  57
  58       variable = variable_in
  59       IF ( variable_in .ge. 'A' .and. variable_in .le. 'Z' ) THEN
  60         variable = CHAR( ICHAR(variable_in) - ICHAR('A') + ICHAR('a') )
  61       ENDIF
  62       IF ((variable == 'u') .or. (variable == 'v') .or.  &
  63           (variable == 'w') .or. (variable == 't') .or.  &
  64           (variable == 'x') .or. (variable == 'y') .or.  &
  65           (variable == 'r') .or. (variable == 'p') ) open_bc_copy = .true.
  66
  67 !  begin, first set a staggering variable
  68
  69       istag = -1
  70       jstag = -1
  71
  72       IF ((variable == 'u') .or. (variable == 'x')) istag = 0
  73       IF ((variable == 'v') .or. (variable == 'y')) jstag = 0
  74
  75       if(debug) then
  76         write(6,*) ' in bc2d, var is ',variable, istag, jstag
  77         write(6,*) ' b.cs are ',  &
  78       config_flags%periodic_x,  &
  79       config_flags%periodic_y
  80       end if
  81
  82       IF ( variable == 'd' ) then  !JDM
  83          istag = 0
  84          jstag = 0
  85       ENDIF
  86       IF ( variable == 'e' ) then  !JDM
  87          istag = 0
  88       ENDIF
  89       IF ( variable == 'f' ) then  !JDM
  90          jstag = 0
  91       ENDIF
  92
  93 !  fix corners for doubly periodic domains
  94
  95       IF ( config_flags%periodic_x .and. config_flags%periodic_y &
  96            .and. (ids == ips) .and. (ide == ipe)                 &
  97            .and. (jds == jps) .and. (jde == jpe)                   ) THEN
  98
  99          IF ( (its == ids) .and. (jte == jde) ) THEN  ! upper left corner fill
 100            DO j = bdyzone,1,-1
 101            DO i = -(bdyzone-1),0,1
 102              a_aux = a_dat(ids+i-1,jde+j+jstag)
 103              a_dat(ids+i-1,jde+j+jstag) = 0.0
 104              a_dat(ide+i-1,jds+j+jstag) = a_dat(ide+i-1,jds+j+jstag) + a_aux
 105              a_aux = 0.0
 106            ENDDO
 107            ENDDO
 108          END IF
 109
 110          IF ( (ite == ide) .and. (jte == jde) ) THEN  ! upper right corner fill
 111            DO j = bdyzone,1,-1
 112            DO i = bdyzone,1,-1
 113              a_aux = a_dat(ide+i+istag,jde+j+jstag)
 114              a_dat(ide+i+istag,jde+j+jstag) = 0.0
 115              a_dat(ids+i+istag,jds+j+jstag) = a_dat(ids+i+istag,jds+j+jstag) + a_aux
 116              a_aux = 0.0
 117            ENDDO
 118            ENDDO
 119          END IF
 120
 121          IF ( (ite == ide) .and. (jts == jds) ) THEN  ! lower right corner fill
 122            DO j = -(bdyzone-1),0,1
 123            DO i = bdyzone,1,-1
 124              a_aux = a_dat(ide+i+istag,jds+j-1)
 125              a_dat(ide+i+istag,jds+j-1) = 0.0
 126              a_dat(ids+i+istag,jde+j-1) = a_dat(ids+i+istag,jde+j-1) + a_aux
 127              a_aux = 0.0
 128            ENDDO
 129            ENDDO
 130          END IF
 131
 132          IF ( (its == ids) .and. (jts == jds) ) THEN  ! lower left corner fill
 133            DO j = -(bdyzone-1),0,1
 134            DO i = -(bdyzone-1),0,1
 135              a_aux = a_dat(ids+i-1,jds+j-1)
 136              a_dat(ids+i-1,jds+j-1) = 0.0
 137              a_dat(ide+i-1,jde+j-1) = a_dat(ide+i-1,jde+j-1) + a_aux
 138              a_aux = 0.0
 139            ENDDO
 140            ENDDO
 141          END IF
 142
 143        END IF
 144
 145 !  same procedure in y
 146
 147       periodicity_y:  IF( ( config_flags%periodic_y ) ) THEN
 148
 149         IF ( ( jds == jps ) .and. ( jde == jpe ) )  THEN    ! test of both north and south on processor
 150
 151           IF( jte == jde ) then
 152
 153             DO j = bdyzone,-jstag,-1
 154             DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 155               a_aux = a_dat(i,jde+j+jstag)
 156               a_dat(i,jde+j+jstag) = 0.0
 157               a_dat(i,jds+j+jstag) = a_dat(i,jds+j+jstag) + a_aux
 158               a_aux = 0.0
 159             ENDDO
 160             ENDDO
 161
 162           END IF
 163
 164           IF( jts == jds ) then
 165
 166             DO j = -(bdyzone-1),0,1
 167             DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 168               a_aux = a_dat(i,jds+j-1)
 169               a_dat(i,jds+j-1) = 0.0
 170               a_dat(i,jde+j-1) = a_dat(i,jde+j-1) + a_aux
 171               a_aux = 0.0
 172             ENDDO
 173             ENDDO
 174
 175           END IF
 176
 177         END IF
 178
 179       ELSE
 180
 181 !  set open b.c in Y copy into boundary zone here.  WCS, 19 March 2000
 182
 183 !  now the open boundary copy at ye
 184
 185         open_ye: IF( ( config_flags%open_ye   .or. &
 186                        config_flags%polar     .or. &
 187                        config_flags%specified .or. &
 188                        config_flags%nested            ) .and.  &
 189                          ( jte == jde ) .and. open_bc_copy )  THEN
 190
 191           IF  (variable /= 'v' .and. variable /= 'y' ) THEN
 192
 193             DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 194               a_aux = a_dat(i,jde  )
 195               a_dat(i,jde  ) = 0.0
 196               a_dat(i,jde-1) = a_dat(i,jde-1) + a_aux
 197               a_aux = a_dat(i,jde+1)
 198               a_dat(i,jde+1) = 0.0
 199               a_dat(i,jde-1) = a_dat(i,jde-1) + a_aux
 200               a_aux = a_dat(i,jde+2)
 201               a_dat(i,jde+2) = 0.0
 202               a_dat(i,jde-1) = a_dat(i,jde-1) + a_aux
 203               a_aux = 0.0
 204             ENDDO
 205
 206           ELSE
 207
 208             DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 209               a_aux = a_dat(i,jde+1)
 210               a_dat(i,jde+1) = 0.0
 211               a_dat(i,jde  ) = a_dat(i,jde  ) + a_aux
 212               a_aux = a_dat(i,jde+2)
 213               a_dat(i,jde+2) = 0.0
 214               a_dat(i,jde  ) = a_dat(i,jde  ) + a_aux
 215               a_aux = a_dat(i,jde+3)
 216               a_dat(i,jde+3) = 0.0
 217               a_dat(i,jde  ) = a_dat(i,jde  ) + a_aux
 218               a_aux = 0.0
 219             ENDDO
 220
 221           ENDIF
 222
 223         END IF open_ye
 224
 225         open_ys: IF( ( config_flags%open_ys   .or. &
 226                        config_flags%polar     .or. &
 227                        config_flags%specified .or. &
 228                        config_flags%nested            ) .and.  &
 229                          ( jts == jds) .and. open_bc_copy )  THEN
 230
 231             DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 232               a_aux = a_dat(i,jds-1)
 233               a_dat(i,jds-1) = 0.0
 234               a_dat(i,jds) = a_dat(i,jds) + a_aux
 235               a_aux = a_dat(i,jds-2)
 236               a_dat(i,jds-2) = 0.0
 237               a_dat(i,jds) = a_dat(i,jds) + a_aux
 238               a_aux = a_dat(i,jds-3)
 239               a_dat(i,jds-3) = 0.0
 240               a_dat(i,jds) = a_dat(i,jds) + a_aux
 241               a_aux = 0.0
 242             ENDDO
 243
 244         ENDIF open_ys
 245
 246 !  now the symmetry boundary at ye
 247
 248         symmetry_ye: IF( ( config_flags%symmetric_ye ) .and.  &
 249                          ( jte == jde )                  )  THEN
 250
 251           IF ( (variable /= 'v') .and. (variable /= 'y') ) THEN
 252
 253             DO j = bdyzone,1,-1
 254             DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 255               a_aux = a_dat(i,jde+j-1)
 256               a_dat(i,jde+j-1) = 0.0
 257               a_dat(i,jde-j) = a_dat(i,jde-j) + a_aux
 258               a_aux = 0.0
 259             ENDDO
 260             ENDDO
 261
 262           ELSE
 263
 264             IF (variable == 'v' ) THEN
 265
 266               DO j = bdyzone,1,-1
 267               DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 268                 a_aux = - a_dat(i,jde+j)
 269                 a_dat(i,jde+j) = 0.0
 270                 a_dat(i,jde-j) = a_dat(i,jde-j) + a_aux
 271                 a_aux = 0.0
 272               ENDDO
 273               ENDDO
 274
 275             ELSE
 276
 277               DO j = bdyzone,1,-1
 278               DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 279                 a_aux = a_dat(i,jde+j)
 280                 a_dat(i,jde+j) = 0.0
 281                 a_dat(i,jde-j) = a_dat(i,jde-j) + a_aux
 282                 a_aux = 0.0
 283               ENDDO
 284               ENDDO
 285
 286             END IF
 287
 288           ENDIF
 289
 290         END IF symmetry_ye
 291
 292         symmetry_ys: IF( ( config_flags%symmetric_ys ) .and.  &
 293                          ( jts == jds)                   )  THEN
 294
 295           IF ( (variable /= 'v') .and. (variable /= 'y') ) THEN
 296
 297             DO j = bdyzone,1,-1
 298             DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 299               a_aux = a_dat(i,jds-j)
 300               a_dat(i,jds-j) = 0.0
 301               a_dat(i,jds+j-1) = a_dat(i,jds+j-1) + a_aux
 302               a_aux = 0.0
 303             ENDDO
 304             ENDDO
 305
 306           ELSE
 307
 308             IF (variable == 'v') THEN
 309
 310               DO j = bdyzone,1,-1
 311               DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 312                 a_aux = - a_dat(i,jds-j)
 313                 a_dat(i,jds-j) = 0.0
 314                 a_dat(i,jds+j) = a_dat(i,jds+j) + a_aux
 315                 a_aux = 0.0
 316               ENDDO
 317               ENDDO
 318
 319             ELSE
 320
 321               DO j = bdyzone,1,-1
 322               DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 323                 a_aux = a_dat(i,jds-j)
 324                 a_dat(i,jds-j) = 0.0
 325                 a_dat(i,jds+j) = a_dat(i,jds+j) + a_aux
 326                 a_aux = 0.0
 327               ENDDO
 328               ENDDO
 329
 330             END IF
 331
 332           ENDIF
 333
 334         ENDIF symmetry_ys
 335
 336       END IF periodicity_y
 337
 338       periodicity_x:  IF( ( config_flags%periodic_x ) ) THEN
 339
 340         IF ( ( ids == ips ) .and.  ( ide == ipe ) ) THEN  ! test if east and west both on-processor
 341           IF ( ite == ide ) THEN
 342
 343             DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
 344             DO i = bdyzone,-istag,-1
 345               a_aux = a_dat(ide+i+istag,j)
 346               a_dat(ide+i+istag,j) = 0.0
 347               a_dat(ids+i+istag,j) = a_dat(ids+i+istag,j) + a_aux
 348               a_aux = 0.0
 349             ENDDO
 350             ENDDO
 351
 352           ENDIF
 353
 354           IF ( its == ids ) THEN
 355
 356             DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
 357             DO i = -(bdyzone-1),0,1
 358               a_aux = a_dat(ids+i-1,j)
 359               a_dat(ids+i-1,j) = 0.0
 360               a_dat(ide+i-1,j) = a_dat(ide+i-1,j) + a_aux
 361               a_aux = 0.0
 362             ENDDO
 363             ENDDO
 364
 365           ENDIF
 366         ENDIF
 367
 368       ELSE
 369
 370 !  set open b.c in X copy into boundary zone here.  WCS, 19 March 2000
 371
 372 !  now the open boundary copy at xe
 373
 374         open_xe: IF( ( config_flags%open_xe   .or. &
 375                        config_flags%specified .or. &
 376                        config_flags%nested            ) .and.  &
 377                           ( ite == ide ) .and. open_bc_copy  )  THEN
 378
 379           IF ( variable /= 'u' .and. variable /= 'x') THEN
 380
 381             DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
 382               a_aux = a_dat(ide  ,j)
 383               a_dat(ide  ,j) = 0.0
 384               a_dat(ide-1,j) = a_dat(ide-1,j) + a_aux
 385               a_aux = a_dat(ide+1,j)
 386               a_dat(ide+1,j) = 0.0
 387               a_dat(ide-1,j) = a_dat(ide-1,j) + a_aux
 388               a_aux = a_dat(ide+2,j)
 389               a_dat(ide+2,j) = 0.0
 390               a_dat(ide-1,j) = a_dat(ide-1,j) + a_aux
 391               a_aux = 0.0
 392             ENDDO
 393
 394           ELSE
 395
 396             DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
 397               a_aux = a_dat(ide+1,j)
 398               a_dat(ide+1,j) = 0.0
 399               a_dat(ide,j) = a_dat(ide,j) + a_aux
 400               a_aux = a_dat(ide+2,j)
 401               a_dat(ide+2,j) = 0.0
 402               a_dat(ide,j) = a_dat(ide,j) + a_aux
 403               a_aux = a_dat(ide+3,j)
 404               a_dat(ide+3,j) = 0.0
 405               a_dat(ide,j) = a_dat(ide,j) + a_aux
 406               a_aux = 0.0
 407             ENDDO
 408
 409           END IF
 410
 411         END IF open_xe
 412
 413         open_xs: IF( ( config_flags%open_xs   .or. &
 414                        config_flags%specified .or. &
 415                        config_flags%nested            ) .and.  &
 416                          ( its == ids ) .and. open_bc_copy  )  THEN
 417
 418             DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
 419               a_aux = a_dat(ids-1,j)
 420               a_dat(ids-1,j) = 0.0
 421               a_dat(ids,j) = a_dat(ids,j) + a_aux
 422               a_aux = a_dat(ids-2,j)
 423               a_dat(ids-2,j) = 0.0
 424               a_dat(ids,j) = a_dat(ids,j) + a_aux
 425               a_aux = a_dat(ids-3,j)
 426               a_dat(ids-3,j) = 0.0
 427               a_dat(ids,j) = a_dat(ids,j) + a_aux
 428               a_aux = 0.0
 429             ENDDO
 430
 431         ENDIF open_xs
 432
 433 !  end open b.c in X copy into boundary zone addition.  WCS, 19 March 2000
 434
 435 !  now the symmetry boundary at xe
 436
 437         symmetry_xe: IF( ( config_flags%symmetric_xe ) .and.  &
 438                          ( ite == ide )                  )  THEN
 439
 440           IF ( (variable /= 'u') .and. (variable /= 'x') ) THEN
 441
 442             DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
 443             DO i = bdyzone,1,-1
 444               a_aux = a_dat(ide+i-1,j)
 445               a_dat(ide+i-1,j) = 0.0
 446               a_dat(ide-i,j) = a_dat(ide-i,j) + a_aux
 447               a_aux = 0.0
 448             ENDDO
 449             ENDDO
 450
 451           ELSE
 452
 453             IF (variable == 'u' ) THEN
 454
 455               DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
 456               DO i = bdyzone-1,0,-1
 457                 a_aux = - a_dat(ide+i,j)
 458                 a_dat(ide+i,j) = 0.0
 459                 a_dat(ide-i,j) = a_dat(ide-i,j) + a_aux
 460                 a_aux = 0.0
 461               ENDDO
 462               ENDDO
 463
 464
 465             ELSE
 466
 467               DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
 468               DO i = bdyzone-1,0,-1
 469                 a_aux = a_dat(ide+i,j)
 470                 a_dat(ide+i,j) = 0.0
 471                 a_dat(ide-i,j) = a_dat(ide-i,j) + a_aux
 472                 a_aux = 0.0
 473               ENDDO
 474               ENDDO
 475
 476             END IF
 477
 478           END IF
 479
 480         END IF symmetry_xe
 481
 482         symmetry_xs: IF( ( config_flags%symmetric_xs ) .and.  &
 483                          ( its == ids )                  )  THEN
 484
 485           IF ( (variable /= 'u') .and. (variable /= 'x') ) THEN
 486
 487             DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
 488             DO i = bdyzone,1,-1
 489               a_aux = a_dat(ids-i,j)
 490               a_dat(ids-i,j) = 0.0
 491               a_dat(ide+i-1,j) = a_dat(ide+i-1,j) + a_aux
 492               a_aux = 0.0
 493             ENDDO
 494             ENDDO
 495
 496           ELSE
 497
 498             IF( variable == 'u' ) THEN
 499
 500               DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
 501               DO i = bdyzone-1,0,-1
 502                 a_aux = - a_dat(ids-i,j)
 503                 a_dat(ids-i,j) = 0.0
 504                 a_dat(ids+i,j) = a_dat(ids+i,j) + a_aux
 505                 a_aux = 0.0
 506               ENDDO
 507               ENDDO
 508
 509             ELSE
 510
 511               DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
 512               DO i = bdyzone-1,0,-1
 513                 a_aux = a_dat(ids-i,j)
 514                 a_dat(ids-i,j) = 0.0
 515                 a_dat(ids+i,j) = a_dat(ids+i,j) + a_aux
 516                 a_aux = 0.0
 517               ENDDO
 518               ENDDO
 519
 520             END IF
 521
 522           ENDIF
 523
 524         ENDIF symmetry_xs
 525
 526       END IF periodicity_x
 527
 528    END SUBROUTINE a_set_physical_bc2d
 529
 530 !-----------------------------------
 531
 532    SUBROUTINE a_set_physical_bc3d( a_dat, variable_in,        &
 533                                config_flags,                   &
 534                                ids,ide, jds,jde, kds,kde,  & ! domain dims
 535                                ims,ime, jms,jme, kms,kme,  & ! memory dims
 536                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
 537                                its,ite, jts,jte, kts,kte )
 538
 539       IMPLICIT NONE
 540
 541       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
 542       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
 543       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
 544       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
 545       CHARACTER,    INTENT(IN   )    :: variable_in
 546
 547       CHARACTER                      :: variable
 548
 549       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ) :: a_dat
 550       TYPE( grid_config_rec_type ) config_flags
 551
 552       INTEGER  :: i, j, k, istag, jstag, itime, k_end
 553
 554       LOGICAL  :: debug, open_bc_copy
 555
 556       real :: a_aux
 557
 558 !------------
 559
 560       a_aux = 0.0
 561
 562       debug = .false.
 563
 564       open_bc_copy = .false.
 565
 566       variable = variable_in
 567       IF ( variable_in .ge. 'A' .and. variable_in .le. 'Z' ) THEN
 568         variable = CHAR( ICHAR(variable_in) - ICHAR('A') + ICHAR('a') )
 569       ENDIF
 570
 571       IF ((variable == 'u') .or. (variable == 'v') .or.     &
 572           (variable == 'w') .or. (variable == 't') .or.     &
 573           (variable == 'd') .or. (variable == 'e') .or. &
 574           (variable == 'x') .or. (variable == 'y') .or. &
 575           (variable == 'f') .or. (variable == 'r') .or. &
 576           (variable == 'p')                        ) open_bc_copy = .true.
 577
 578 !  begin, first set a staggering variable
 579
 580       istag = -1
 581       jstag = -1
 582       k_end = max(1,min(kde-1,kte))
 583
 584
 585       IF ((variable == 'u') .or. (variable == 'x')) istag = 0
 586       IF ((variable == 'v') .or. (variable == 'y')) jstag = 0
 587       IF ((variable == 'd') .or. (variable == 'xy')) then
 588          istag = 0
 589          jstag = 0
 590       ENDIF
 591       IF ((variable == 'e') ) then
 592          istag = 0
 593          k_end = min(kde,kte)
 594       ENDIF
 595
 596       IF ((variable == 'f') ) then
 597          jstag = 0
 598          k_end = min(kde,kte)
 599       ENDIF
 600
 601       IF ( variable == 'w')  k_end = min(kde,kte)
 602
 603 !      k_end = kte
 604
 605       if(debug) then
 606         write(6,*) ' in bc, var is ',variable, istag, jstag, kte, k_end
 607         write(6,*) ' b.cs are ',  &
 608       config_flags%periodic_x,  &
 609       config_flags%periodic_y
 610       end if
 611
 612 !  fix corners for doubly periodic domains
 613
 614       IF ( config_flags%periodic_x .and. config_flags%periodic_y &
 615            .and. (ids == ips) .and. (ide == ipe)                 &
 616            .and. (jds == jps) .and. (jde == jpe)                   ) THEN
 617
 618          IF ( (its == ids) .and. (jte == jde) ) THEN  ! upper left corner fill
 619            DO j = bdyzone,1,-1
 620            DO k = kts, k_end
 621            DO i = -(bdyzone-1),0,1
 622              a_aux = a_dat(ids+i-1,k,jde+j+jstag)
 623              a_dat(ids+i-1,k,jde+j+jstag) = 0.0
 624              a_dat(ide+i-1,k,jds+j+jstag) = a_dat(ide+i-1,k,jds+j+jstag) + a_aux
 625              a_aux = 0.0
 626            ENDDO
 627            ENDDO
 628            ENDDO
 629          END IF
 630
 631          IF ( (ite == ide) .and. (jte == jde) ) THEN  ! upper right corner fill
 632            DO j = bdyzone,1,-1
 633            DO k = kts, k_end
 634            DO i = bdyzone,1,-1
 635              a_aux = a_dat(ide+i+istag,k,jde+j+jstag)
 636              a_dat(ide+i+istag,k,jde+j+jstag) = 0.0
 637              a_dat(ids+i+istag,k,jds+j+jstag) = a_dat(ids+i+istag,k,jds+j+jstag) + a_aux
 638              a_aux = 0.0
 639            ENDDO
 640            ENDDO
 641            ENDDO
 642          END IF
 643
 644          IF ( (ite == ide) .and. (jts == jds) ) THEN  ! lower right corner fill
 645            DO j = -(bdyzone-1),0,1
 646            DO k = kts, k_end
 647            DO i = bdyzone,1,-1
 648              a_aux = a_dat(ide+i+istag,k,jds+j-1)
 649              a_dat(ide+i+istag,k,jds+j-1) = 0.0
 650              a_dat(ids+i+istag,k,jde+j-1) = a_dat(ids+i+istag,k,jde+j-1) + a_aux
 651              a_aux = 0.0
 652            ENDDO
 653            ENDDO
 654            ENDDO
 655          END IF
 656
 657          IF ( (its == ids) .and. (jts == jds) ) THEN  ! lower left corner fill
 658            DO j = -(bdyzone-1),0,1
 659            DO k = kts, k_end
 660            DO i = -(bdyzone-1),0,1
 661              a_aux = a_dat(ids+i-1,k,jds+j-1)
 662              a_dat(ids+i-1,k,jds+j-1) = 0.0
 663              a_dat(ide+i-1,k,jde+j-1) = a_dat(ide+i-1,k,jde+j-1) + a_aux
 664              a_aux = 0.0
 665            ENDDO
 666            ENDDO
 667            ENDDO
 668          END IF
 669
 670        END IF
 671
 672 !  same procedure in y
 673
 674       periodicity_y:  IF( ( config_flags%periodic_y ) ) THEN
 675         IF ( ( jds == jps ) .and. ( jde == jpe ) )  THEN      ! test if both north and south on processor
 676           IF( jte == jde ) then
 677
 678             DO j = bdyzone,-jstag,-1
 679             DO k = kts, k_end
 680             DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 681               a_aux = a_dat(i,k,jde+j+jstag)
 682               a_dat(i,k,jde+j+jstag) = 0.0
 683               a_dat(i,k,jds+j+jstag) = a_dat(i,k,jds+j+jstag) + a_aux
 684               a_aux = 0.0
 685             ENDDO
 686             ENDDO
 687             ENDDO
 688
 689           END IF
 690
 691           IF( jts == jds ) then
 692
 693             DO j = -(bdyzone-1),0,1
 694             DO k = kts, k_end
 695             DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 696               a_aux = a_dat(i,k,jds+j-1)
 697               a_dat(i,k,jds+j-1) = 0.0
 698               a_dat(i,k,jde+j-1) = a_dat(i,k,jde+j-1) + a_aux
 699               a_aux = 0.0
 700             ENDDO
 701             ENDDO
 702             ENDDO
 703
 704           END IF
 705
 706         END IF
 707
 708       ELSE
 709
 710 !  set open b.c in Y copy into boundary zone here.  WCS, 19 March 2000
 711
 712 !  now the open boundary copy at ye
 713
 714         open_ye: IF( ( config_flags%open_ye   .or. &
 715                        config_flags%polar     .or. &
 716                        config_flags%specified .or. &
 717                        config_flags%nested            ) .and.  &
 718                          ( jte == jde ) .and. open_bc_copy )  THEN
 719
 720           IF (variable /= 'v' .and. variable /= 'y' ) THEN
 721
 722             DO k = kts, k_end
 723             DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 724               a_aux = a_dat(i,k,jde  )
 725               a_dat(i,k,jde  ) = 0.0
 726               a_dat(i,k,jde-1) = a_dat(i,k,jde-1) + a_aux
 727               a_aux = a_dat(i,k,jde+1)
 728               a_dat(i,k,jde+1) = 0.0
 729               a_dat(i,k,jde-1) = a_dat(i,k,jde-1) + a_aux
 730               a_aux = a_dat(i,k,jde+2)
 731               a_dat(i,k,jde+2) = 0.0
 732               a_dat(i,k,jde-1) = a_dat(i,k,jde-1) + a_aux
 733               a_aux = 0.0
 734             ENDDO
 735             ENDDO
 736
 737           ELSE
 738
 739             DO k = kts, k_end
 740             DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 741               a_aux = a_dat(i,k,jde+1)
 742               a_dat(i,k,jde+1) = 0.0
 743               a_dat(i,k,jde) = a_dat(i,k,jde) + a_aux
 744               a_aux = a_dat(i,k,jde+2)
 745               a_dat(i,k,jde+2) = 0.0
 746               a_dat(i,k,jde) = a_dat(i,k,jde) + a_aux
 747               a_aux = a_dat(i,k,jde+3)
 748               a_dat(i,k,jde+3) = 0.0
 749               a_dat(i,k,jde) = a_dat(i,k,jde) + a_aux
 750               a_aux = 0.0
 751             ENDDO
 752             ENDDO
 753
 754           ENDIF
 755
 756         END IF open_ye
 757
 758         open_ys: IF( ( config_flags%open_ys   .or. &
 759                        config_flags%polar     .or. &
 760                        config_flags%specified .or. &
 761                        config_flags%nested            ) .and.  &
 762                          ( jts == jds) .and. open_bc_copy )  THEN
 763
 764             DO k = kts, k_end
 765             DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 766               a_aux = a_dat(i,k,jds-1)
 767               a_dat(i,k,jds-1) = 0.0
 768               a_dat(i,k,jds) = a_dat(i,k,jds) + a_aux
 769               a_aux = a_dat(i,k,jds-2)
 770               a_dat(i,k,jds-2) = 0.0
 771               a_dat(i,k,jds) = a_dat(i,k,jds) + a_aux
 772               a_aux = a_dat(i,k,jds-3)
 773               a_dat(i,k,jds-3) = 0.0
 774               a_dat(i,k,jds) = a_dat(i,k,jds) + a_aux
 775               a_aux = 0.0
 776             ENDDO
 777             ENDDO
 778
 779         ENDIF open_ys
 780
 781 !  end open b.c in Y copy into boundary zone addition.  WCS, 19 March 2000
 782
 783 !  now the symmetry boundary at ye
 784
 785         symmetry_ye: IF( ( config_flags%symmetric_ye ) .and.  &
 786                          ( jte == jde )                  )  THEN
 787
 788           IF ( (variable /= 'v') .and. (variable /= 'y') ) THEN
 789
 790             DO j = bdyzone,1,-1
 791             DO k = kts, k_end
 792             DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 793               a_aux = a_dat(i,k,jde+j-1)
 794               a_dat(i,k,jde+j-1) = 0.0
 795               a_dat(i,k,jde-j) = a_dat(i,k,jde-j) + a_aux
 796               a_aux = 0.0
 797             ENDDO
 798             ENDDO
 799             ENDDO
 800
 801           ELSE
 802
 803             IF ( variable == 'v' ) THEN
 804
 805               DO j = bdyzone,1,-1
 806               DO k = kts, k_end
 807               DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 808                 a_aux = - a_dat(i,k,jde+j)
 809                 a_dat(i,k,jde+j) = 0.0
 810                 a_dat(i,k,jde-j) = a_dat(i,k,jde-j) + a_aux
 811                 a_aux = 0.0
 812               ENDDO
 813               ENDDO
 814               ENDDO
 815
 816             ELSE
 817
 818               DO j = bdyzone,1,-1
 819               DO k = kts, k_end
 820               DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 821                 a_aux = a_dat(i,k,jde+j)
 822                 a_dat(i,k,jde+j) = 0.0
 823                 a_dat(i,k,jde-j) = a_dat(i,k,jde-j) + a_aux
 824                 a_aux = 0.0
 825               ENDDO
 826               ENDDO
 827               ENDDO
 828
 829             END IF
 830
 831           ENDIF
 832
 833         END IF symmetry_ye
 834
 835         symmetry_ys: IF( ( config_flags%symmetric_ys ) .and.  &
 836                          ( jts == jds)                   )  THEN
 837
 838           IF ( (variable /= 'v') .and. (variable /= 'y') ) THEN
 839
 840             DO j = bdyzone,1,-1
 841             DO k = kts, k_end
 842             DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 843               a_aux = a_dat(i,k,jds-j)
 844               a_dat(i,k,jds-j) = 0.0
 845               a_dat(i,k,jds+j-1) = a_dat(i,k,jds+j-1) + a_aux
 846               a_aux = 0.0
 847             ENDDO
 848             ENDDO
 849             ENDDO
 850
 851           ELSE
 852
 853             IF (variable == 'v') THEN
 854
 855               DO j = bdyzone,1,-1
 856               DO k = kts, k_end
 857               DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 858                 a_aux = - a_dat(i,k,jds-j)
 859                 a_dat(i,k,jds-j) = 0.0
 860                 a_dat(i,k,jds+j) = a_dat(i,k,jds+j) + a_aux
 861                 a_aux = 0.0
 862               ENDDO
 863               ENDDO
 864               ENDDO
 865
 866             ELSE
 867
 868               DO j = bdyzone,1,-1
 869               DO k = kts, k_end
 870               DO i = MIN(ite+1,ide+istag),MAX(ids,its-1),-1
 871                 a_aux = a_dat(i,k,jds-j)
 872                 a_dat(i,k,jds-j) = 0.0
 873                 a_dat(i,k,jds+j) = a_dat(i,k,jds+j) + a_aux
 874                 a_aux = 0.0
 875               ENDDO
 876               ENDDO
 877               ENDDO
 878
 879             END IF
 880
 881           ENDIF
 882
 883         ENDIF symmetry_ys
 884
 885       END IF periodicity_y
 886
 887       periodicity_x:  IF( ( config_flags%periodic_x ) ) THEN
 888
 889         IF ( ( ids == ips ) .and. ( ide == ipe ) ) THEN  ! test if both east and west on-processor
 890
 891           IF ( ite == ide ) THEN
 892
 893             DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
 894             DO k = kts, k_end
 895             DO i = bdyzone,-istag,-1
 896               a_aux = a_dat(ide+i+istag,k,j)
 897               a_dat(ide+i+istag,k,j) = 0.0
 898               a_dat(ids+i+istag,k,j) = a_dat(ids+i+istag,k,j) + a_aux
 899               a_aux = 0.0
 900             ENDDO
 901             ENDDO
 902             ENDDO
 903
 904           ENDIF
 905
 906           IF ( its == ids ) THEN
 907
 908             DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
 909             DO k = kts, k_end
 910             DO i = -(bdyzone-1),0,1
 911               a_aux = a_dat(ids+i-1,k,j)
 912               a_dat(ids+i-1,k,j) = 0.0
 913               a_dat(ide+i-1,k,j) = a_dat(ide+i-1,k,j) + a_aux
 914               a_aux = 0.0
 915             ENDDO
 916             ENDDO
 917             ENDDO
 918
 919           ENDIF
 920
 921         ENDIF
 922
 923       ELSE
 924
 925 !  set open b.c in X copy into boundary zone here.  WCS, 19 March 2000
 926
 927 !  now the open_xe boundary copy
 928
 929         open_xe: IF( ( config_flags%open_xe   .or. &
 930                        config_flags%specified .or. &
 931                        config_flags%nested            ) .and.  &
 932                          ( ite == ide ) .and. open_bc_copy )  THEN
 933
 934           IF (variable /= 'u' .and. variable /= 'x' ) THEN
 935
 936             DO j = MIN(jte,jde+jstag)+bdyzone,jts-bdyzone,-1
 937             DO k = kts, k_end
 938               a_aux = a_dat(ide  ,k,j)
 939               a_dat(ide  ,k,j) = 0.0
 940               a_dat(ide-1,k,j) = a_dat(ide-1,k,j) + a_aux
 941               a_aux = a_dat(ide+1,k,j)
 942               a_dat(ide+1,k,j) = 0.0
 943               a_dat(ide-1,k,j) = a_dat(ide-1,k,j) + a_aux
 944               a_aux = a_dat(ide+2,k,j)
 945               a_dat(ide+2,k,j) = 0.0
 946               a_dat(ide-1,k,j) = a_dat(ide-1,k,j) + a_aux
 947               a_aux = 0.0
 948             ENDDO
 949             ENDDO
 950
 951           ELSE
 952
 953 !!!!!!! I am not sure about this one!  JM 20020402
 954             DO j = MIN(jte+1,jde+jstag)+bdyzone,MAX(jds,jts-1)-bdyzone,-1
 955             DO k = kts, k_end
 956               a_aux = a_dat(ide+1,k,j)
 957               a_dat(ide+1,k,j) = 0.0
 958               a_dat(ide,k,j) = a_dat(ide,k,j) + a_aux
 959               a_aux = a_dat(ide+2,k,j)
 960               a_dat(ide+2,k,j) = 0.0
 961               a_dat(ide,k,j) = a_dat(ide,k,j) + a_aux
 962               a_aux = a_dat(ide+3,k,j)
 963               a_dat(ide+3,k,j) = 0.0
 964               a_dat(ide,k,j) = a_dat(ide,k,j) + a_aux
 965               a_aux = 0.0
 966             ENDDO
 967             ENDDO
 968
 969           END IF
 970
 971         END IF open_xe
 972
 973         open_xs: IF( ( config_flags%open_xs   .or. &
 974                        config_flags%specified .or. &
 975                        config_flags%nested            ) .and.  &
 976                          ( its == ids ) .and. open_bc_copy  )  THEN
 977
 978             DO j = MIN(jte,jde+jstag)+bdyzone,jts-bdyzone,-1
 979             DO k = kts, k_end
 980               a_aux = a_dat(ids-1,k,j)
 981               a_dat(ids-1,k,j) = 0.0
 982               a_dat(ids,k,j) = a_dat(ids,k,j) + a_aux
 983               a_aux = a_dat(ids-2,k,j)
 984               a_dat(ids-2,k,j) = 0.0
 985               a_dat(ids,k,j) = a_dat(ids,k,j) + a_aux
 986               a_aux = a_dat(ids-3,k,j)
 987               a_dat(ids-3,k,j) = 0.0
 988               a_dat(ids,k,j) = a_dat(ids,k,j) + a_aux
 989               a_aux = 0.0
 990             ENDDO
 991             ENDDO
 992
 993         ENDIF open_xs
 994
 995 !  end open b.c in X copy into boundary zone addition.  WCS, 19 March 2000
 996
 997 !  now the symmetry boundary at xe
 998
 999         symmetry_xe: IF( ( config_flags%symmetric_xe ) .and.  &
1000                          ( ite == ide )                  )  THEN
1001
1002           IF ( (variable /= 'u') .and. (variable /= 'x') ) THEN
1003
1004             DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
1005             DO k = kts, k_end
1006             DO i = bdyzone,1,-1
1007               a_aux = a_dat(ide+i-1,k,j)
1008               a_dat(ide+i-1,k,j) = 0.0
1009               a_dat(ide-i,k,j) = a_dat(ide-i,k,j) + a_aux
1010               a_aux = 0.0
1011             ENDDO
1012             ENDDO
1013             ENDDO
1014
1015           ELSE
1016
1017             IF (variable == 'u') THEN
1018
1019               DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
1020               DO k = kts, k_end
1021               DO i = bdyzone,1,-1
1022                 a_aux = - a_dat(ide+i,k,j)
1023                 a_dat(ide+i,k,j) = 0.0
1024                 a_dat(ide-i,k,j) = a_dat(ide-i,k,j) + a_aux
1025                 a_aux = 0.0
1026               ENDDO
1027               ENDDO
1028               ENDDO
1029
1030             ELSE
1031
1032               DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
1033               DO k = kts, k_end
1034               DO i = bdyzone,1,-1
1035                 a_aux = a_dat(ide+i,k,j)
1036                 a_dat(ide+i,k,j) = 0.0
1037                 a_dat(ide-i,k,j) = a_dat(ide-i,k,j) + a_aux
1038                 a_aux = 0.0
1039               ENDDO
1040               ENDDO
1041               ENDDO
1042
1043              END IF
1044
1045           END IF
1046
1047         END IF symmetry_xe
1048
1049         symmetry_xs: IF( ( config_flags%symmetric_xs ) .and.  &
1050                          ( its == ids )                  )  THEN
1051
1052           IF ( (variable /= 'u') .and. (variable /= 'x') ) THEN
1053
1054             DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
1055             DO k = kts, k_end
1056             DO i = bdyzone,1,-1
1057               a_aux = a_dat(ids-i,k,j)
1058               a_dat(ids-i,k,j) = 0.0
1059               a_dat(ids+i-1,k,j) = a_dat(ids+i-1,k,j) + a_aux
1060               a_aux = 0.0
1061             ENDDO
1062             ENDDO
1063             ENDDO
1064
1065           ELSE
1066
1067             IF ( variable == 'u' ) THEN
1068
1069               DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
1070               DO k = kts, k_end
1071               DO i = bdyzone,1,-1
1072                 a_aux = - a_dat(ids-i,k,j)
1073                 a_dat(ids-i,k,j) = 0.0
1074                 a_dat(ids+i,k,j) = a_dat(ids+i,k,j) + a_aux
1075                 a_aux = 0.0
1076               ENDDO
1077               ENDDO
1078               ENDDO
1079
1080             ELSE
1081
1082               DO j = MIN(jte+1,jde+jstag),MAX(jds,jts-1),-1
1083               DO k = kts, k_end
1084               DO i = bdyzone,1,-1
1085                 a_aux = a_dat(ids-i,k,j)
1086                 a_dat(ids-i,k,j) = 0.0
1087                 a_dat(ids+i,k,j) = a_dat(ids+i,k,j) + a_aux
1088                 a_aux = 0.0
1089               ENDDO
1090               ENDDO
1091               ENDDO
1092
1093             END IF
1094
1095           ENDIF
1096
1097         ENDIF symmetry_xs
1098
1099       END IF periodicity_x
1100
1101    END SUBROUTINE a_set_physical_bc3d
1102
1103 !------------------------------------------------------------------------
1104
1105    SUBROUTINE a_init_module_bc
1106    END SUBROUTINE a_init_module_bc
1107
1108 !------------------------------------------------------------------------
1109
1110 ! a couple versions of this call to allow a smaller-than-memory dimensioned field (e.g. tile sized) ! to be passed in as the first argument.  Both of these call the _core version defined below.
1111    SUBROUTINE a_relax_bdytend ( a_field, a_field_tend, &
1112                                 a_field_bdy_xs, a_field_bdy_xe,  &
1113                                 a_field_bdy_ys, a_field_bdy_ye,  &
1114                                 a_field_bdy_tend_xs, a_field_bdy_tend_xe,  &
1115                                 a_field_bdy_tend_ys, a_field_bdy_tend_ye,  &
1116                                 variable_in, config_flags,             &
1117                                 spec_bdy_width, spec_zone, relax_zone, &
1118                                 dtbc, fcx, gcx,             &
1119                                 ids,ide, jds,jde, kds,kde,  & ! domain dims
1120                                 ims,ime, jms,jme, kms,kme,  & ! memory dims
1121                                 ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1122                                 its,ite, jts,jte, kts,kte )
1123
1124    IMPLICIT NONE
1125
1126    INTEGER,   INTENT(IN) :: ids,ide, jds,jde, kds,kde
1127    INTEGER,   INTENT(IN) :: ims,ime, jms,jme, kms,kme
1128    INTEGER,   INTENT(IN) :: ips,ipe, jps,jpe, kps,kpe
1129    INTEGER,   INTENT(IN) :: its,ite, jts,jte, kts,kte
1130    INTEGER,   INTENT(IN) :: spec_bdy_width, spec_zone, relax_zone
1131    REAL,      INTENT(IN) :: dtbc
1132    CHARACTER, INTENT(IN) :: variable_in
1133
1134    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: a_field
1135    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN   ) :: a_field_tend
1136    REAL, DIMENSION(jms:jme, kds:kde, spec_bdy_width), INTENT(INOUT) :: a_field_bdy_xs, a_field_bdy_xe
1137    REAL, DIMENSION(ims:ime, kds:kde, spec_bdy_width), INTENT(INOUT) :: a_field_bdy_ys, a_field_bdy_ye
1138    REAL, DIMENSION(jms:jme, kds:kde, spec_bdy_width), INTENT(INOUT) :: a_field_bdy_tend_xs, a_field_bdy_tend_xe
1139    REAL, DIMENSION(ims:ime, kds:kde, spec_bdy_width), INTENT(INOUT) :: a_field_bdy_tend_ys, a_field_bdy_tend_ye
1140    REAL, DIMENSION(spec_bdy_width), INTENT(IN) :: fcx, gcx
1141    TYPE(grid_config_rec_type),      INTENT(IN) :: config_flags
1142
1143
1144    CALL a_relax_bdytend_core ( a_field, a_field_tend, &
1145                     a_field_bdy_xs, a_field_bdy_xe, &
1146                     a_field_bdy_ys, a_field_bdy_ye, &
1147                     a_field_bdy_tend_xs, a_field_bdy_tend_xe,  &
1148                     a_field_bdy_tend_ys, a_field_bdy_tend_ye,  &
1149                     variable_in, config_flags,             &
1150                     spec_bdy_width, spec_zone, relax_zone, &
1151                     dtbc, fcx, gcx,             &
1152                     ids,ide, jds,jde, kds,kde,  & ! domain dims
1153                     ims,ime, jms,jme, kms,kme,  & ! memory dims
1154                     ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1155                     its,ite, jts,jte, kts,kte,  & ! patch  dims
1156                     ims,ime, jms,jme, kms,kme )  ! dimension of the field argument
1157
1158    END SUBROUTINE a_relax_bdytend
1159
1160 ! version that allows tile-sized version of field. Note, caller should define the
1161 ! field to be -+1 of tile size in each dimension because routine is going off onto halo
1162 ! for example, see relax_bdytend in dyn_em/module_bc_em.F
1163    SUBROUTINE a_relax_bdytend_tile ( a_field, a_field_tend, &
1164                        a_field_bdy_xs, a_field_bdy_xe,  &
1165                        a_field_bdy_ys, a_field_bdy_ye,  &
1166                        a_field_bdy_tend_xs, a_field_bdy_tend_xe,  &
1167                        a_field_bdy_tend_ys, a_field_bdy_tend_ye,  &
1168                        variable_in, config_flags,             &
1169                        spec_bdy_width, spec_zone, relax_zone, &
1170                        dtbc, fcx, gcx,             &
1171                        ids,ide, jds,jde, kds,kde,  & ! domain dims
1172                        ims,ime, jms,jme, kms,kme,  & ! memory dims
1173                        ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1174                        its,ite, jts,jte, kts,kte,  &
1175                        iXs,iXe, jXs,jXe, kXs,kXe   &  ! dims of first argument
1176                        )
1177
1178    IMPLICIT NONE
1179
1180    INTEGER,   INTENT(IN) :: ids,ide, jds,jde, kds,kde
1181    INTEGER,   INTENT(IN) :: ims,ime, jms,jme, kms,kme
1182    INTEGER,   INTENT(IN) :: ips,ipe, jps,jpe, kps,kpe
1183    INTEGER,   INTENT(IN) :: its,ite, jts,jte, kts,kte
1184    INTEGER,   INTENT(IN) :: iXs,iXe, jXs,jXe, kXs,kXe
1185    INTEGER,   INTENT(IN) :: spec_bdy_width, spec_zone, relax_zone
1186    REAL,      INTENT(IN) :: dtbc
1187    CHARACTER, INTENT(IN) :: variable_in
1188
1189    REAL, DIMENSION(iXs:iXe, kXs:kXe, jXs:jXe), INTENT(INOUT) :: a_field
1190    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN   ) :: a_field_tend
1191    REAL, DIMENSION(jms:jme, kds:kde, spec_bdy_width), INTENT(INOUT) :: a_field_bdy_xs, a_field_bdy_xe
1192    REAL, DIMENSION(ims:ime, kds:kde, spec_bdy_width), INTENT(INOUT) :: a_field_bdy_ys, a_field_bdy_ye
1193    REAL, DIMENSION(jms:jme, kds:kde, spec_bdy_width), INTENT(INOUT) :: a_field_bdy_tend_xs, a_field_bdy_tend_xe
1194    REAL, DIMENSION(ims:ime, kds:kde, spec_bdy_width), INTENT(INOUT) :: a_field_bdy_tend_ys, a_field_bdy_tend_ye
1195    REAL, DIMENSION(spec_bdy_width), INTENT(IN) :: fcx, gcx
1196    TYPE(grid_config_rec_type),      INTENT(IN) :: config_flags
1197
1198
1199    CALL a_relax_bdytend_core ( a_field, a_field_tend, &
1200                     a_field_bdy_xs, a_field_bdy_xe,  &
1201                     a_field_bdy_ys, a_field_bdy_ye,  &
1202                     a_field_bdy_tend_xs, a_field_bdy_tend_xe,  &
1203                     a_field_bdy_tend_ys, a_field_bdy_tend_ye,  &
1204                     variable_in, config_flags,             &
1205                     spec_bdy_width, spec_zone, relax_zone, &
1206                     dtbc, fcx, gcx,             &
1207                     ids,ide, jds,jde, kds,kde,  & ! domain dims
1208                     ims,ime, jms,jme, kms,kme,  & ! memory dims
1209                     ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1210                     its,ite, jts,jte, kts,kte,  &
1211                     iXs,iXe, jXs,jXe, kXs,kXe )  ! dimension of the field argument
1212
1213    END SUBROUTINE a_relax_bdytend_tile
1214
1215 !        Generated by TAPENADE     (INRIA, Tropics team)
1216 !  Tapenade 3.5 (r3785) - 22 Mar 2011 18:35
1217 !
1218 !  Differentiation of relax_bdytend_core in reverse (adjoint) mode:
1219 !   gradient     of useful results: field field_bdy_xe field_bdy_tend_xe
1220 !                field_tend field_bdy_xs field_bdy_tend_xs field_bdy_ye
1221 !                field_bdy_tend_ye field_bdy_ys field_bdy_tend_ys
1222 !   with respect to varying inputs: field field_bdy_xe field_bdy_tend_xe
1223 !                field_tend field_bdy_xs field_bdy_tend_xs field_bdy_ye
1224 !                field_bdy_tend_ye field_bdy_ys field_bdy_tend_ys
1225 !   RW status of diff variables: field:incr field_bdy_xe:incr field_bdy_tend_xe:incr
1226 !                field_tend:in-out field_bdy_xs:incr field_bdy_tend_xs:incr
1227 !                field_bdy_ye:incr field_bdy_tend_ye:incr field_bdy_ys:incr
1228 !                field_bdy_tend_ys:incr
1229 ! domain dims
1230 ! memory dims
1231 ! patch  dims
1232 ! patch  dims
1233 ! field (1st arg) dims; might be tile or patch
1234 SUBROUTINE A_RELAX_BDYTEND_CORE(fieldb, field_tendb, &
1235 &  field_bdy_xsb, field_bdy_xeb&
1236 &  , field_bdy_ysb, field_bdy_yeb, &
1237 &  field_bdy_tend_xsb, field_bdy_tend_xeb, &
1238 &  field_bdy_tend_ysb, &
1239 &  field_bdy_tend_yeb, variable_in, config_flags, spec_bdy_width, &
1240 &  spec_zone, relax_zone, dtbc, fcx, gcx, ids, ide, jds, jde, kds, kde, &
1241 &  ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe, its, ite, &
1242 &  jts, jte, kts, kte, ixs, ixe, jxs, jxe, kxs, kxe)
1243   IMPLICIT NONE
1244   INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde
1245   INTEGER, INTENT(IN) :: ims, ime, jms, jme, kms, kme
1246   INTEGER, INTENT(IN) :: ips, ipe, jps, jpe, kps, kpe
1247   INTEGER, INTENT(IN) :: its, ite, jts, jte, kts, kte
1248   INTEGER, INTENT(IN) :: ixs, ixe, jxs, jxe, kxs, kxe
1249   INTEGER, INTENT(IN) :: spec_bdy_width, spec_zone, relax_zone
1250   REAL, INTENT(IN) :: dtbc
1251   CHARACTER, INTENT(IN) :: variable_in
1252   REAL, DIMENSION(ixs:ixe, kxs:kxe, jxs:jxe) :: fieldb
1253   REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: field_tendb
1254   REAL, DIMENSION(jms:jme, kds:kde, spec_bdy_width) :: field_bdy_xsb, &
1255 &  field_bdy_xeb
1256   REAL, DIMENSION(ims:ime, kds:kde, spec_bdy_width) :: field_bdy_ysb, &
1257 &  field_bdy_yeb
1258   REAL, DIMENSION(jms:jme, kds:kde, spec_bdy_width) :: &
1259 &  field_bdy_tend_xsb, field_bdy_tend_xeb
1260   REAL, DIMENSION(ims:ime, kds:kde, spec_bdy_width) :: &
1261 &  field_bdy_tend_ysb, field_bdy_tend_yeb
1262   REAL, DIMENSION(spec_bdy_width), INTENT(IN) :: fcx, gcx
1263   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   REAL :: fls0b, fls1b, fls2b, fls3b, fls4b
1269   LOGICAL :: periodic_x
1270   INTEGER :: branch
1271   INTEGER :: ad_from
1272   INTEGER :: ad_to
1273   INTEGER :: ad_from0
1274   INTEGER :: ad_to0
1275   INTEGER :: ad_from1
1276   INTEGER :: ad_to1
1277   INTEGER :: ad_from2
1278   INTEGER :: ad_to2
1279   INTEGER :: min8
1280   INTEGER :: min7
1281   INTEGER :: min6
1282   INTEGER :: min5
1283   INTEGER :: min4
1284   INTEGER :: min3
1285   INTEGER :: min2
1286   INTEGER :: min1
1287   REAL :: tempb2
1288   REAL :: tempb1
1289   REAL :: tempb0
1290   REAL :: tempb
1291   INTEGER :: max8
1292   INTEGER :: max7
1293   INTEGER :: max6
1294   INTEGER :: max5
1295   INTEGER :: max4
1296   INTEGER :: max3
1297   INTEGER :: max2
1298   INTEGER :: max1
1299   periodic_x = config_flags%periodic_x
1300   variable = variable_in
1301   IF (variable .EQ. 'U') variable = 'u'
1302   IF (variable .EQ. 'V') variable = 'v'
1303   IF (variable .EQ. 'M') variable = 'm'
1304   IF (variable .EQ. 'H') variable = 'h'
1305   ibs = ids
1306   ibe = ide - 1
1307   IF (ite .GT. ide - 1) THEN
1308     itf = ide - 1
1309   ELSE
1310     itf = ite
1311   END IF
1312   jbs = jds
1313   jbe = jde - 1
1314   IF (jte .GT. jde - 1) THEN
1315     jtf = jde - 1
1316   ELSE
1317     jtf = jte
1318   END IF
1319   ktf = kde - 1
1320   IF (variable .EQ. 'u') ibe = ide
1321   IF (variable .EQ. 'u') THEN
1322     IF (ite .GT. ide) THEN
1323       itf = ide
1324     ELSE
1325       itf = ite
1326     END IF
1327   END IF
1328   IF (variable .EQ. 'v') jbe = jde
1329   IF (variable .EQ. 'v') THEN
1330     IF (jte .GT. jde) THEN
1331       jtf = jde
1332     ELSE
1333       jtf = jte
1334     END IF
1335   END IF
1336   IF (variable .EQ. 'm') ktf = kte
1337   IF (variable .EQ. 'h') ktf = kte
1338   IF (jts - jbs .LT. relax_zone) THEN
1339     IF (jts .LT. jbs + spec_zone) THEN
1340       max1 = jbs + spec_zone
1341     ELSE
1342       max1 = jts
1343     END IF
1344     IF (jtf .GT. jbs + relax_zone - 1) THEN
1345       min1 = jbs + relax_zone - 1
1346     ELSE
1347       min1 = jtf
1348     END IF
1349 ! Y-start boundary
1350     DO j=max1,min1
1351       CALL PUSHINTEGER4(b_dist)
1352       b_dist = j - jbs
1353       b_limit = b_dist
1354       IF (periodic_x) b_limit = 0
1355       DO k=kts,ktf
1356         IF (its .LT. b_limit + ibs) THEN
1357           max2 = b_limit + ibs
1358         ELSE
1359           max2 = its
1360         END IF
1361         IF (itf .GT. ibe - b_limit) THEN
1362           min2 = ibe - b_limit
1363         ELSE
1364           min2 = itf
1365         END IF
1366         ad_from = max2
1367         DO i=ad_from,min2
1368           IF (i - 1 .LT. ibs) THEN
1369             CALL PUSHINTEGER4(im1)
1370             im1 = ibs
1371             CALL PUSHCONTROL1B(0)
1372           ELSE
1373             CALL PUSHINTEGER4(im1)
1374             im1 = i - 1
1375             CALL PUSHCONTROL1B(1)
1376           END IF
1377           IF (i + 1 .GT. ibe) THEN
1378             CALL PUSHINTEGER4(ip1)
1379             ip1 = ibe
1380             CALL PUSHCONTROL1B(0)
1381           ELSE
1382             CALL PUSHINTEGER4(ip1)
1383             ip1 = i + 1
1384             CALL PUSHCONTROL1B(1)
1385           END IF
1386         END DO
1387         CALL PUSHINTEGER4(i - 1)
1388         CALL PUSHINTEGER4(ad_from)
1389       END DO
1390     END DO
1391     CALL PUSHCONTROL1B(0)
1392   ELSE
1393     CALL PUSHCONTROL1B(1)
1394   END IF
1395   IF (jbe - jtf .LT. relax_zone) THEN
1396     IF (jts .LT. jbe - relax_zone + 1) THEN
1397       max3 = jbe - relax_zone + 1
1398     ELSE
1399       max3 = jts
1400     END IF
1401     IF (jtf .GT. jbe - spec_zone) THEN
1402       min3 = jbe - spec_zone
1403     ELSE
1404       min3 = jtf
1405     END IF
1406 ! Y-end boundary
1407     DO j=max3,min3
1408       CALL PUSHINTEGER4(b_dist)
1409       b_dist = jbe - j
1410       b_limit = b_dist
1411       IF (periodic_x) b_limit = 0
1412       DO k=kts,ktf
1413         IF (its .LT. b_limit + ibs) THEN
1414           max4 = b_limit + ibs
1415         ELSE
1416           max4 = its
1417         END IF
1418         IF (itf .GT. ibe - b_limit) THEN
1419           min4 = ibe - b_limit
1420         ELSE
1421           min4 = itf
1422         END IF
1423         ad_from0 = max4
1424         DO i=ad_from0,min4
1425           IF (i - 1 .LT. ibs) THEN
1426             CALL PUSHINTEGER4(im1)
1427             im1 = ibs
1428             CALL PUSHCONTROL1B(0)
1429           ELSE
1430             CALL PUSHINTEGER4(im1)
1431             im1 = i - 1
1432             CALL PUSHCONTROL1B(1)
1433           END IF
1434           IF (i + 1 .GT. ibe) THEN
1435             CALL PUSHINTEGER4(ip1)
1436             ip1 = ibe
1437             CALL PUSHCONTROL1B(0)
1438           ELSE
1439             CALL PUSHINTEGER4(ip1)
1440             ip1 = i + 1
1441             CALL PUSHCONTROL1B(1)
1442           END IF
1443         END DO
1444         CALL PUSHINTEGER4(i - 1)
1445         CALL PUSHINTEGER4(ad_from0)
1446       END DO
1447     END DO
1448     CALL PUSHCONTROL1B(0)
1449   ELSE
1450     CALL PUSHCONTROL1B(1)
1451   END IF
1452   IF (.NOT.periodic_x) THEN
1453     IF (its - ibs .LT. relax_zone) THEN
1454       IF (its .LT. ibs + spec_zone) THEN
1455         max5 = ibs + spec_zone
1456       ELSE
1457         max5 = its
1458       END IF
1459       IF (itf .GT. ibs + relax_zone - 1) THEN
1460         min5 = ibs + relax_zone - 1
1461       ELSE
1462         min5 = itf
1463       END IF
1464 ! X-start boundary
1465       DO i=max5,min5
1466         CALL PUSHINTEGER4(b_dist)
1467         b_dist = i - ibs
1468         DO k=kts,ktf
1469           IF (jts .LT. b_dist + jbs + 1) THEN
1470             max6 = b_dist + jbs + 1
1471           ELSE
1472             max6 = jts
1473           END IF
1474           IF (jtf .GT. jbe - b_dist - 1) THEN
1475             min6 = jbe - b_dist - 1
1476           ELSE
1477             min6 = jtf
1478           END IF
1479           ad_from1 = max6
1480           j = min6 + 1
1481           CALL PUSHINTEGER4(j - 1)
1482           CALL PUSHINTEGER4(ad_from1)
1483         END DO
1484       END DO
1485       CALL PUSHCONTROL1B(0)
1486     ELSE
1487       CALL PUSHCONTROL1B(1)
1488     END IF
1489     IF (ibe - itf .LT. relax_zone) THEN
1490       IF (its .LT. ibe - relax_zone + 1) THEN
1491         max7 = ibe - relax_zone + 1
1492       ELSE
1493         max7 = its
1494       END IF
1495       IF (itf .GT. ibe - spec_zone) THEN
1496         min7 = ibe - spec_zone
1497       ELSE
1498         min7 = itf
1499       END IF
1500 ! X-end boundary
1501       DO i=max7,min7
1502         CALL PUSHINTEGER4(b_dist)
1503         b_dist = ibe - i
1504         DO k=kts,ktf
1505           IF (jts .LT. b_dist + jbs + 1) THEN
1506             max8 = b_dist + jbs + 1
1507           ELSE
1508             max8 = jts
1509           END IF
1510           IF (jtf .GT. jbe - b_dist - 1) THEN
1511             min8 = jbe - b_dist - 1
1512           ELSE
1513             min8 = jtf
1514           END IF
1515           ad_from2 = max8
1516           j = min8 + 1
1517           CALL PUSHINTEGER4(j - 1)
1518           CALL PUSHINTEGER4(ad_from2)
1519         END DO
1520       END DO
1521       DO i=min7,max7,-1
1522         DO k=ktf,kts,-1
1523           CALL POPINTEGER4(ad_from2)
1524           CALL POPINTEGER4(ad_to2)
1525           DO j=ad_to2,ad_from2,-1
1526             tempb2 = -(gcx(b_dist+1)*field_tendb(i, k, j))
1527             fls0b = fcx(b_dist+1)*field_tendb(i, k, j) - 4.*tempb2
1528             fls1b = tempb2
1529             fls2b = tempb2
1530             fls3b = tempb2
1531             fls4b = tempb2
1532             field_bdy_xeb(j, k, b_dist+2) = field_bdy_xeb(j, k, b_dist+2&
1533 &              ) + fls4b
1534             field_bdy_tend_xeb(j, k, b_dist+2) = field_bdy_tend_xeb(j, k&
1535 &              , b_dist+2) + dtbc*fls4b
1536             fieldb(i-1, k, j) = fieldb(i-1, k, j) - fls4b
1537             field_bdy_xeb(j, k, b_dist) = field_bdy_xeb(j, k, b_dist) + &
1538 &              fls3b
1539             field_bdy_tend_xeb(j, k, b_dist) = field_bdy_tend_xeb(j, k, &
1540 &              b_dist) + dtbc*fls3b
1541             fieldb(i+1, k, j) = fieldb(i+1, k, j) - fls3b
1542             field_bdy_xeb(j+1, k, b_dist+1) = field_bdy_xeb(j+1, k, &
1543 &              b_dist+1) + fls2b
1544             field_bdy_tend_xeb(j+1, k, b_dist+1) = field_bdy_tend_xeb(j+&
1545 &              1, k, b_dist+1) + dtbc*fls2b
1546             fieldb(i, k, j+1) = fieldb(i, k, j+1) - fls2b
1547             field_bdy_xeb(j-1, k, b_dist+1) = field_bdy_xeb(j-1, k, &
1548 &              b_dist+1) + fls1b
1549             field_bdy_tend_xeb(j-1, k, b_dist+1) = field_bdy_tend_xeb(j-&
1550 &              1, k, b_dist+1) + dtbc*fls1b
1551             fieldb(i, k, j-1) = fieldb(i, k, j-1) - fls1b
1552             field_bdy_xeb(j, k, b_dist+1) = field_bdy_xeb(j, k, b_dist+1&
1553 &              ) + fls0b
1554             field_bdy_tend_xeb(j, k, b_dist+1) = field_bdy_tend_xeb(j, k&
1555 &              , b_dist+1) + dtbc*fls0b
1556             fieldb(i, k, j) = fieldb(i, k, j) - fls0b
1557           END DO
1558         END DO
1559         CALL POPINTEGER4(b_dist)
1560       END DO
1561     END IF
1562     CALL POPCONTROL1B(branch)
1563     IF (branch .EQ. 0) THEN
1564       DO i=min5,max5,-1
1565         DO k=ktf,kts,-1
1566           CALL POPINTEGER4(ad_from1)
1567           CALL POPINTEGER4(ad_to1)
1568           DO j=ad_to1,ad_from1,-1
1569             tempb1 = -(gcx(b_dist+1)*field_tendb(i, k, j))
1570             fls0b = fcx(b_dist+1)*field_tendb(i, k, j) - 4.*tempb1
1571             fls1b = tempb1
1572             fls2b = tempb1
1573             fls3b = tempb1
1574             fls4b = tempb1
1575             field_bdy_xsb(j, k, b_dist+2) = field_bdy_xsb(j, k, b_dist+2&
1576 &              ) + fls4b
1577             field_bdy_tend_xsb(j, k, b_dist+2) = field_bdy_tend_xsb(j, k&
1578 &              , b_dist+2) + dtbc*fls4b
1579             fieldb(i+1, k, j) = fieldb(i+1, k, j) - fls4b
1580             field_bdy_xsb(j, k, b_dist) = field_bdy_xsb(j, k, b_dist) + &
1581 &              fls3b
1582             field_bdy_tend_xsb(j, k, b_dist) = field_bdy_tend_xsb(j, k, &
1583 &              b_dist) + dtbc*fls3b
1584             fieldb(i-1, k, j) = fieldb(i-1, k, j) - fls3b
1585             field_bdy_xsb(j+1, k, b_dist+1) = field_bdy_xsb(j+1, k, &
1586 &              b_dist+1) + fls2b
1587             field_bdy_tend_xsb(j+1, k, b_dist+1) = field_bdy_tend_xsb(j+&
1588 &              1, k, b_dist+1) + dtbc*fls2b
1589             fieldb(i, k, j+1) = fieldb(i, k, j+1) - fls2b
1590             field_bdy_xsb(j-1, k, b_dist+1) = field_bdy_xsb(j-1, k, &
1591 &              b_dist+1) + fls1b
1592             field_bdy_tend_xsb(j-1, k, b_dist+1) = field_bdy_tend_xsb(j-&
1593 &              1, k, b_dist+1) + dtbc*fls1b
1594             fieldb(i, k, j-1) = fieldb(i, k, j-1) - fls1b
1595             field_bdy_xsb(j, k, b_dist+1) = field_bdy_xsb(j, k, b_dist+1&
1596 &              ) + fls0b
1597             field_bdy_tend_xsb(j, k, b_dist+1) = field_bdy_tend_xsb(j, k&
1598 &              , b_dist+1) + dtbc*fls0b
1599             fieldb(i, k, j) = fieldb(i, k, j) - fls0b
1600           END DO
1601         END DO
1602         CALL POPINTEGER4(b_dist)
1603       END DO
1604     END IF
1605   END IF
1606   CALL POPCONTROL1B(branch)
1607   IF (branch .EQ. 0) THEN
1608     DO j=min3,max3,-1
1609       DO k=ktf,kts,-1
1610         CALL POPINTEGER4(ad_from0)
1611         CALL POPINTEGER4(ad_to0)
1612         DO i=ad_to0,ad_from0,-1
1613           tempb0 = -(gcx(b_dist+1)*field_tendb(i, k, j))
1614           fls0b = fcx(b_dist+1)*field_tendb(i, k, j) - 4.*tempb0
1615           fls1b = tempb0
1616           fls2b = tempb0
1617           fls3b = tempb0
1618           fls4b = tempb0
1619           field_bdy_yeb(i, k, b_dist+2) = field_bdy_yeb(i, k, b_dist+2) &
1620 &            + fls4b
1621           field_bdy_tend_yeb(i, k, b_dist+2) = field_bdy_tend_yeb(i, k, &
1622 &            b_dist+2) + dtbc*fls4b
1623           fieldb(i, k, j-1) = fieldb(i, k, j-1) - fls4b
1624           field_bdy_yeb(i, k, b_dist) = field_bdy_yeb(i, k, b_dist) + &
1625 &            fls3b
1626           field_bdy_tend_yeb(i, k, b_dist) = field_bdy_tend_yeb(i, k, &
1627 &            b_dist) + dtbc*fls3b
1628           fieldb(i, k, j+1) = fieldb(i, k, j+1) - fls3b
1629           field_bdy_yeb(ip1, k, b_dist+1) = field_bdy_yeb(ip1, k, b_dist&
1630 &            +1) + fls2b
1631           field_bdy_tend_yeb(ip1, k, b_dist+1) = field_bdy_tend_yeb(ip1&
1632 &            , k, b_dist+1) + dtbc*fls2b
1633           fieldb(ip1, k, j) = fieldb(ip1, k, j) - fls2b
1634           field_bdy_yeb(im1, k, b_dist+1) = field_bdy_yeb(im1, k, b_dist&
1635 &            +1) + fls1b
1636           field_bdy_tend_yeb(im1, k, b_dist+1) = field_bdy_tend_yeb(im1&
1637 &            , k, b_dist+1) + dtbc*fls1b
1638           fieldb(im1, k, j) = fieldb(im1, k, j) - fls1b
1639           field_bdy_yeb(i, k, b_dist+1) = field_bdy_yeb(i, k, b_dist+1) &
1640 &            + fls0b
1641           field_bdy_tend_yeb(i, k, b_dist+1) = field_bdy_tend_yeb(i, k, &
1642 &            b_dist+1) + dtbc*fls0b
1643           fieldb(i, k, j) = fieldb(i, k, j) - fls0b
1644           CALL POPCONTROL1B(branch)
1645           IF (branch .EQ. 0) THEN
1646             CALL POPINTEGER4(ip1)
1647           ELSE
1648             CALL POPINTEGER4(ip1)
1649           END IF
1650           CALL POPCONTROL1B(branch)
1651           IF (branch .EQ. 0) THEN
1652             CALL POPINTEGER4(im1)
1653           ELSE
1654             CALL POPINTEGER4(im1)
1655           END IF
1656         END DO
1657       END DO
1658       CALL POPINTEGER4(b_dist)
1659     END DO
1660   END IF
1661   CALL POPCONTROL1B(branch)
1662   IF (branch .EQ. 0) THEN
1663     DO j=min1,max1,-1
1664       DO k=ktf,kts,-1
1665         CALL POPINTEGER4(ad_from)
1666         CALL POPINTEGER4(ad_to)
1667         DO i=ad_to,ad_from,-1
1668           tempb = -(gcx(b_dist+1)*field_tendb(i, k, j))
1669           fls0b = fcx(b_dist+1)*field_tendb(i, k, j) - 4.*tempb
1670           fls1b = tempb
1671           fls2b = tempb
1672           fls3b = tempb
1673           fls4b = tempb
1674           field_bdy_ysb(i, k, b_dist+2) = field_bdy_ysb(i, k, b_dist+2) &
1675 &            + fls4b
1676           field_bdy_tend_ysb(i, k, b_dist+2) = field_bdy_tend_ysb(i, k, &
1677 &            b_dist+2) + dtbc*fls4b
1678           fieldb(i, k, j+1) = fieldb(i, k, j+1) - fls4b
1679           field_bdy_ysb(i, k, b_dist) = field_bdy_ysb(i, k, b_dist) + &
1680 &            fls3b
1681           field_bdy_tend_ysb(i, k, b_dist) = field_bdy_tend_ysb(i, k, &
1682 &            b_dist) + dtbc*fls3b
1683           fieldb(i, k, j-1) = fieldb(i, k, j-1) - fls3b
1684           field_bdy_ysb(ip1, k, b_dist+1) = field_bdy_ysb(ip1, k, b_dist&
1685 &            +1) + fls2b
1686           field_bdy_tend_ysb(ip1, k, b_dist+1) = field_bdy_tend_ysb(ip1&
1687 &            , k, b_dist+1) + dtbc*fls2b
1688           fieldb(ip1, k, j) = fieldb(ip1, k, j) - fls2b
1689           field_bdy_ysb(im1, k, b_dist+1) = field_bdy_ysb(im1, k, b_dist&
1690 &            +1) + fls1b
1691           field_bdy_tend_ysb(im1, k, b_dist+1) = field_bdy_tend_ysb(im1&
1692 &            , k, b_dist+1) + dtbc*fls1b
1693           fieldb(im1, k, j) = fieldb(im1, k, j) - fls1b
1694           field_bdy_ysb(i, k, b_dist+1) = field_bdy_ysb(i, k, b_dist+1) &
1695 &            + fls0b
1696           field_bdy_tend_ysb(i, k, b_dist+1) = field_bdy_tend_ysb(i, k, &
1697 &            b_dist+1) + dtbc*fls0b
1698           fieldb(i, k, j) = fieldb(i, k, j) - fls0b
1699           CALL POPCONTROL1B(branch)
1700           IF (branch .EQ. 0) THEN
1701             CALL POPINTEGER4(ip1)
1702           ELSE
1703             CALL POPINTEGER4(ip1)
1704           END IF
1705           CALL POPCONTROL1B(branch)
1706           IF (branch .EQ. 0) THEN
1707             CALL POPINTEGER4(im1)
1708           ELSE
1709             CALL POPINTEGER4(im1)
1710           END IF
1711         END DO
1712       END DO
1713       CALL POPINTEGER4(b_dist)
1714     END DO
1715   END IF
1716 END SUBROUTINE A_RELAX_BDYTEND_CORE
1717 !------------------------------------------------------------------------
1718
1719    SUBROUTINE a_spec_bdytend ( a_field_tend,           &
1720                                a_field_bdy_tend_xs, a_field_bdy_tend_xe, &
1721                                a_field_bdy_tend_ys, a_field_bdy_tend_ye, &
1722                                variable_in, config_flags, &
1723                                spec_bdy_width, spec_zone, &
1724                                ids,ide, jds,jde, kds,kde,  & ! domain dims
1725                                ims,ime, jms,jme, kms,kme,  & ! memory dims
1726                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1727                                its,ite, jts,jte, kts,kte )
1728
1729 !  spec_bdy_width is only used to dimension the boundary arrays.
1730 !  spec_zone is the width of the outer specified b.c.s that are set here.
1731
1732       IMPLICIT NONE
1733
1734       INTEGER,   INTENT(IN) :: ids,ide, jds,jde, kds,kde
1735       INTEGER,   INTENT(IN) :: ims,ime, jms,jme, kms,kme
1736       INTEGER,   INTENT(IN) :: ips,ipe, jps,jpe, kps,kpe
1737       INTEGER,   INTENT(IN) :: its,ite, jts,jte, kts,kte
1738       INTEGER,   INTENT(IN) :: spec_bdy_width, spec_zone
1739       CHARACTER, INTENT(IN) :: variable_in
1740
1741       REAL, DIMENSION(ims:ime, kms:kme, jms:jme),        INTENT(INOUT) :: a_field_tend
1742       REAL, DIMENSION(jms:jme, kds:kde, spec_bdy_width), INTENT(INOUT) :: a_field_bdy_tend_xs, a_field_bdy_tend_xe
1743       REAL, DIMENSION(ims:ime, kds:kde, spec_bdy_width), INTENT(INOUT) :: a_field_bdy_tend_ys, a_field_bdy_tend_ye
1744       TYPE( grid_config_rec_type ) config_flags
1745
1746       CHARACTER  :: variable
1747       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf
1748       INTEGER    :: b_dist, b_limit
1749       LOGICAL    :: periodic_x
1750
1751
1752       periodic_x = config_flags%periodic_x
1753
1754       variable = variable_in
1755
1756       IF (variable == 'U') variable = 'u'
1757       IF (variable == 'V') variable = 'v'
1758       IF (variable == 'M') variable = 'm'
1759       IF (variable == 'H') variable = 'h'
1760
1761       ibs = ids
1762       ibe = ide-1
1763       itf = min(ite,ide-1)
1764       jbs = jds
1765       jbe = jde-1
1766       jtf = min(jte,jde-1)
1767       ktf = kde-1
1768       IF (variable == 'u') ibe = ide
1769       IF (variable == 'u') itf = min(ite,ide)
1770       IF (variable == 'v') jbe = jde
1771       IF (variable == 'v') jtf = min(jte,jde)
1772       IF (variable == 'm') ktf = kte
1773       IF (variable == 'h') ktf = kte
1774
1775
1776     IF(.NOT.periodic_x)THEN
1777       IF (ibe - itf .lt. spec_zone) THEN
1778 ! X-end boundary
1779         DO i = itf, max(its,ibe-spec_zone+1), -1
1780           b_dist = ibe - i
1781           DO k = kts, ktf
1782             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1783               a_field_bdy_tend_xe(j, k, b_dist+1) = a_field_bdy_tend_xe(j, k, b_dist+1) &
1784                                                   + a_field_tend(i,k,j)
1785               a_field_tend(i,k,j) = 0.
1786             ENDDO
1787           ENDDO
1788         ENDDO
1789       ENDIF
1790
1791       IF (its - ibs .lt. spec_zone) THEN
1792 ! X-start boundary
1793         DO i = min(itf,ibs+spec_zone-1), its, -1
1794           b_dist = i - ibs
1795           DO k = kts, ktf
1796             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1797               a_field_bdy_tend_xs(j, k, b_dist+1) = a_field_bdy_tend_xs(j, k, b_dist+1) &
1798                                                   + a_field_tend(i,k,j)
1799               a_field_tend(i,k,j) = 0.
1800             ENDDO
1801           ENDDO
1802         ENDDO
1803       ENDIF
1804
1805     ENDIF
1806
1807       IF (jbe - jtf .lt. spec_zone) THEN
1808 ! Y-end boundary
1809         DO j = jtf, max(jts,jbe-spec_zone+1), -1
1810           b_dist = jbe - j
1811           b_limit = b_dist
1812           IF(periodic_x)b_limit = 0
1813           DO k = kts, ktf
1814             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1815               a_field_bdy_tend_ye(i, k, b_dist+1) = a_field_bdy_tend_ye(i, k, b_dist+1) &
1816                                                   + a_field_tend(i,k,j)
1817               a_field_tend(i,k,j) = 0.
1818             ENDDO
1819           ENDDO
1820         ENDDO
1821       ENDIF
1822
1823       IF (jts - jbs .lt. spec_zone) THEN
1824 ! Y-start boundary
1825         DO j = min(jtf,jbs+spec_zone-1), jts, -1
1826           b_dist = j - jbs
1827           b_limit = b_dist
1828           IF(periodic_x)b_limit = 0
1829           DO k = kts, ktf
1830             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1831               a_field_bdy_tend_ys(i, k, b_dist+1) = a_field_bdy_tend_ys(i, k, b_dist+1) &
1832                                                   + a_field_tend(i,k,j)
1833               a_field_tend(i,k,j) = 0.
1834             ENDDO
1835           ENDDO
1836         ENDDO
1837       ENDIF
1838
1839    END SUBROUTINE a_spec_bdytend
1840
1841 !------------------------------------------------------------------------
1842
1843    SUBROUTINE a_spec_bdyupdate(  a_field,  &
1844                                  a_field_tend, dt,            &
1845                                  variable_in, config_flags, &
1846                                  spec_zone,                  &
1847                                  ids,ide, jds,jde, kds,kde,  & ! domain dims
1848                                  ims,ime, jms,jme, kms,kme,  & ! memory dims
1849                                  ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1850                                  its,ite, jts,jte, kts,kte )
1851
1852 !  This subroutine adds the tendencies in the boundary specified region.
1853 !  spec_zone is the width of the outer specified b.c.s that are set here.
1854 !  (JD August 2000)
1855
1856       IMPLICIT NONE
1857
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_zone
1863       CHARACTER,    INTENT(IN   )    :: variable_in
1864       REAL,         INTENT(IN   )    :: dt
1865
1866
1867       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: a_field
1868       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: a_field_tend
1869       TYPE( grid_config_rec_type ) config_flags
1870
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
1875
1876       periodic_x = config_flags%periodic_x
1877
1878       variable = variable_in
1879
1880       IF (variable == 'U') variable = 'u'
1881       IF (variable == 'V') variable = 'v'
1882       IF (variable == 'M') variable = 'm'
1883       IF (variable == 'H') variable = 'h'
1884
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
1898
1899     IF(.NOT.periodic_x)THEN
1900       IF (ibe - itf .lt. spec_zone) THEN
1901 ! X-end boundary
1902         DO i = max(its,ibe-spec_zone+1), itf
1903           b_dist = ibe - i
1904           DO k = kts, ktf
1905             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1906               a_field_tend(i,k,j) = a_field_tend(i,k,j) + dt * a_field(i,k,j)
1907             ENDDO
1908           ENDDO
1909         ENDDO
1910       ENDIF
1911
1912       IF (its - ibs .lt. spec_zone) THEN
1913 ! X-start boundary
1914         DO i = its, min(itf,ibs+spec_zone-1)
1915           b_dist = i - ibs
1916           DO k = kts, ktf
1917             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1918               a_field_tend(i,k,j) = a_field_tend(i,k,j) + dt * a_field(i,k,j)
1919             ENDDO
1920           ENDDO
1921         ENDDO
1922       ENDIF
1923
1924     ENDIF
1925
1926       IF (jbe - jtf .lt. spec_zone) THEN
1927 ! Y-end boundary
1928         DO j = max(jts,jbe-spec_zone+1), jtf
1929           b_dist = jbe - j
1930           b_limit = b_dist
1931           IF(periodic_x)b_limit = 0
1932           DO k = kts, ktf
1933             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1934               a_field_tend(i,k,j) = a_field_tend(i,k,j) + dt * a_field(i,k,j)
1935             ENDDO
1936           ENDDO
1937         ENDDO
1938       ENDIF
1939
1940       IF (jts - jbs .lt. spec_zone) THEN
1941 ! Y-start boundary
1942         DO j = jts, min(jtf,jbs+spec_zone-1)
1943           b_dist = j - jbs
1944           b_limit = b_dist
1945           IF(periodic_x)b_limit = 0
1946           DO k = kts, ktf
1947             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
1948               a_field_tend(i,k,j) = a_field_tend(i,k,j) + dt * a_field(i,k,j)
1949             ENDDO
1950           ENDDO
1951         ENDDO
1952       ENDIF
1953
1954    END SUBROUTINE a_spec_bdyupdate
1955 !------------------------------------------------------------------------
1956 !        Generated by TAPENADE     (INRIA, Tropics team)
1957 !  Tapenade 3.10 (r5498) - 20 Jan 2015 09:48
1958 !
1959 !  Differentiation of spec_bdy_final in reverse (adjoint) mode:
1960 !   gradient     of useful results: field field_bdy_xe field_bdy_tend_xe
1961 !                field_bdy_xs field_bdy_tend_xs field_bdy_ye field_bdy_tend_ye
1962 !                field_bdy_ys field_bdy_tend_ys mu
1963 !   with respect to varying inputs: field field_bdy_xe field_bdy_tend_xe
1964 !                field_bdy_xs field_bdy_tend_xs field_bdy_ye field_bdy_tend_ye
1965 !                field_bdy_ys field_bdy_tend_ys mu
1966 !   RW status of diff variables: field:in-out field_bdy_xe:incr
1967 !                field_bdy_tend_xe:incr field_bdy_xs:incr field_bdy_tend_xs:incr
1968 !                field_bdy_ye:incr field_bdy_tend_ye:incr field_bdy_ys:incr
1969 !                field_bdy_tend_ys:incr mu:incr
1970 ! domain dims
1971 ! memory dims
1972 ! patch  dims
1973 SUBROUTINE a_SPEC_BDY_FINAL(field, fieldb, mu, mub, msf, field_bdy_xs, &
1974 & field_bdy_xsb, field_bdy_xe, field_bdy_xeb, field_bdy_ys, &
1975 & field_bdy_ysb, field_bdy_ye, field_bdy_yeb, field_bdy_tend_xs, &
1976 & field_bdy_tend_xsb, field_bdy_tend_xe, field_bdy_tend_xeb, &
1977 & field_bdy_tend_ys, field_bdy_tend_ysb, field_bdy_tend_ye, &
1978 & field_bdy_tend_yeb, variable_in, config_flags, spec_bdy_width, &
1979 & spec_zone, dtbc, ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms&
1980 & , kme, ips, ipe, jps, jpe, kps, kpe, its, ite, jts, jte, kts, kte)
1981   IMPLICIT NONE
1982   INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde
1983   INTEGER, INTENT(IN) :: ims, ime, jms, jme, kms, kme
1984   INTEGER, INTENT(IN) :: ips, ipe, jps, jpe, kps, kpe
1985   INTEGER, INTENT(IN) :: its, ite, jts, jte, kts, kte
1986   INTEGER, INTENT(IN) :: spec_bdy_width, spec_zone
1987   REAL, INTENT(IN) :: dtbc
1988   CHARACTER, INTENT(IN) :: variable_in
1989   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: field
1990   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: fieldb
1991   REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: mu, msf
1992   REAL, DIMENSION(ims:ime, jms:jme) :: mub
1993   REAL, DIMENSION(jms:jme, kds:kde, spec_bdy_width), INTENT(IN) :: &
1994 & field_bdy_xs, field_bdy_xe
1995   REAL, DIMENSION(jms:jme, kds:kde, spec_bdy_width) :: field_bdy_xsb, &
1996 & field_bdy_xeb
1997   REAL, DIMENSION(ims:ime, kds:kde, spec_bdy_width), INTENT(IN) :: &
1998 & field_bdy_ys, field_bdy_ye
1999   REAL, DIMENSION(ims:ime, kds:kde, spec_bdy_width) :: field_bdy_ysb, &
2000 & field_bdy_yeb
2001   REAL, DIMENSION(jms:jme, kds:kde, spec_bdy_width), INTENT(IN) :: &
2002 & field_bdy_tend_xs, field_bdy_tend_xe
2003   REAL, DIMENSION(jms:jme, kds:kde, spec_bdy_width) :: &
2004 & field_bdy_tend_xsb, field_bdy_tend_xeb
2005   REAL, DIMENSION(ims:ime, kds:kde, spec_bdy_width), INTENT(IN) :: &
2006 & field_bdy_tend_ys, field_bdy_tend_ye
2007   REAL, DIMENSION(ims:ime, kds:kde, spec_bdy_width) :: &
2008 & field_bdy_tend_ysb, field_bdy_tend_yeb
2009   TYPE(GRID_CONFIG_REC_TYPE) :: config_flags
2010   CHARACTER :: variable
2011   INTEGER :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf, im1, ip1
2012   INTEGER :: b_dist, b_limit
2013   REAL :: bfield, xmsf, xmu
2014   REAL :: bfieldb, xmub
2015   LOGICAL :: periodic_x, msfcouple, mucouple
2016   INTEGER :: branch
2017   INTEGER :: ad_from
2018   INTEGER :: ad_to
2019   INTEGER :: ad_from0
2020   INTEGER :: ad_to0
2021   INTEGER :: ad_from1
2022   INTEGER :: ad_to1
2023   INTEGER :: ad_from2
2024   INTEGER :: ad_to2
2025   INTEGER :: min6
2026   INTEGER :: min5
2027   INTEGER :: min4
2028   INTEGER :: min3
2029   INTEGER :: min2
2030   INTEGER :: min1
2031   REAL :: tempb2
2032   REAL :: tempb1
2033   REAL :: tempb0
2034   REAL :: tempb
2035   INTEGER :: max6
2036   INTEGER :: max5
2037   INTEGER :: max4
2038   INTEGER :: max3
2039   INTEGER :: max2
2040   INTEGER :: max1
2041   periodic_x = config_flags%periodic_x
2042   variable = variable_in
2043   IF (variable .EQ. 'U') variable = 'u'
2044   IF (variable .EQ. 'V') variable = 'v'
2045   IF (variable .EQ. 'W') variable = 'w'
2046   IF (variable .EQ. 'M') variable = 'm'
2047   IF (variable .EQ. 'T') variable = 't'
2048   IF (variable .EQ. 'H') variable = 'h'
2049   ibs = ids
2050   ibe = ide - 1
2051   IF (ite .GT. ide - 1) THEN
2052     itf = ide - 1
2053   ELSE
2054     itf = ite
2055   END IF
2056   jbs = jds
2057   jbe = jde - 1
2058   IF (jte .GT. jde - 1) THEN
2059     jtf = jde - 1
2060   ELSE
2061     jtf = jte
2062   END IF
2063   ktf = kde - 1
2064   IF (variable .EQ. 'u') ibe = ide
2065   IF (variable .EQ. 'u') THEN
2066     IF (ite .GT. ide) THEN
2067       itf = ide
2068     ELSE
2069       itf = ite
2070     END IF
2071   END IF
2072   IF (variable .EQ. 'v') jbe = jde
2073   IF (variable .EQ. 'v') THEN
2074     IF (jte .GT. jde) THEN
2075       jtf = jde
2076     ELSE
2077       jtf = jte
2078     END IF
2079   END IF
2080   IF (variable .EQ. 'm') ktf = kde
2081   IF (variable .EQ. 'h') ktf = kde
2082   IF (variable .EQ. 'w') ktf = kde
2083   msfcouple = .false.
2084   mucouple = .true.
2085   IF ((variable .EQ. 'u' .OR. variable .EQ. 'v') .OR. variable .EQ. 'w'&
2086 & ) msfcouple = .true.
2087   IF (variable .EQ. 'm') mucouple = .false.
2088   xmsf = 1.
2089   xmu = 1.
2090   IF (jts - jbs .LT. spec_zone) THEN
2091     IF (jtf .GT. jbs + spec_zone - 1) THEN
2092       min1 = jbs + spec_zone - 1
2093     ELSE
2094       min1 = jtf
2095     END IF
2096 ! Y-start boundary
2097     DO j=jts,min1
2098       CALL PUSHINTEGER4(b_dist)
2099       b_dist = j - jbs
2100       b_limit = b_dist
2101       IF (periodic_x) b_limit = 0
2102       DO k=kts,ktf
2103         IF (its .LT. b_limit + ibs) THEN
2104           max1 = b_limit + ibs
2105         ELSE
2106           max1 = its
2107         END IF
2108         IF (itf .GT. ibe - b_limit) THEN
2109           min3 = ibe - b_limit
2110         ELSE
2111           min3 = itf
2112         END IF
2113         ad_from = max1
2114         DO i=ad_from,min3
2115           IF (msfcouple) THEN
2116             CALL PUSHREAL8(xmsf)
2117             xmsf = msf(i, j)
2118             CALL PUSHCONTROL1B(0)
2119           ELSE
2120             CALL PUSHCONTROL1B(1)
2121           END IF
2122           IF (mucouple) THEN
2123             CALL PUSHREAL8(xmu)
2124             xmu = mu(i, j)
2125             CALL PUSHCONTROL1B(0)
2126           ELSE
2127             CALL PUSHCONTROL1B(1)
2128           END IF
2129         END DO
2130         CALL PUSHINTEGER4(i - 1)
2131         CALL PUSHINTEGER4(ad_from)
2132       END DO
2133     END DO
2134     CALL PUSHCONTROL1B(0)
2135   ELSE
2136     CALL PUSHCONTROL1B(1)
2137   END IF
2138   IF (jbe - jtf .LT. spec_zone) THEN
2139     IF (jts .LT. jbe - spec_zone + 1) THEN
2140       max2 = jbe - spec_zone + 1
2141     ELSE
2142       max2 = jts
2143     END IF
2144 ! Y-end boundary
2145     DO j=max2,jtf
2146       CALL PUSHINTEGER4(b_dist)
2147       b_dist = jbe - j
2148       b_limit = b_dist
2149       IF (periodic_x) b_limit = 0
2150       DO k=kts,ktf
2151         IF (its .LT. b_limit + ibs) THEN
2152           max3 = b_limit + ibs
2153         ELSE
2154           max3 = its
2155         END IF
2156         IF (itf .GT. ibe - b_limit) THEN
2157           min4 = ibe - b_limit
2158         ELSE
2159           min4 = itf
2160         END IF
2161         ad_from0 = max3
2162         DO i=ad_from0,min4
2163           IF (msfcouple) THEN
2164             CALL PUSHREAL8(xmsf)
2165             xmsf = msf(i, j)
2166             CALL PUSHCONTROL1B(0)
2167           ELSE
2168             CALL PUSHCONTROL1B(1)
2169           END IF
2170           IF (mucouple) THEN
2171             CALL PUSHREAL8(xmu)
2172             xmu = mu(i, j)
2173             CALL PUSHCONTROL1B(0)
2174           ELSE
2175             CALL PUSHCONTROL1B(1)
2176           END IF
2177         END DO
2178         CALL PUSHINTEGER4(i - 1)
2179         CALL PUSHINTEGER4(ad_from0)
2180       END DO
2181     END DO
2182     CALL PUSHCONTROL1B(0)
2183   ELSE
2184     CALL PUSHCONTROL1B(1)
2185   END IF
2186   IF (.NOT.periodic_x) THEN
2187     IF (its - ibs .LT. spec_zone) THEN
2188       IF (itf .GT. ibs + spec_zone - 1) THEN
2189         min2 = ibs + spec_zone - 1
2190       ELSE
2191         min2 = itf
2192       END IF
2193 ! X-start boundary
2194       DO i=its,min2
2195         CALL PUSHINTEGER4(b_dist)
2196         b_dist = i - ibs
2197         DO k=kts,ktf
2198           IF (jts .LT. b_dist + jbs + 1) THEN
2199             max4 = b_dist + jbs + 1
2200           ELSE
2201             max4 = jts
2202           END IF
2203           IF (jtf .GT. jbe - b_dist - 1) THEN
2204             min5 = jbe - b_dist - 1
2205           ELSE
2206             min5 = jtf
2207           END IF
2208           ad_from1 = max4
2209           DO j=ad_from1,min5
2210             IF (msfcouple) THEN
2211               CALL PUSHREAL8(xmsf)
2212               xmsf = msf(i, j)
2213               CALL PUSHCONTROL1B(0)
2214             ELSE
2215               CALL PUSHCONTROL1B(1)
2216             END IF
2217             IF (mucouple) THEN
2218               CALL PUSHREAL8(xmu)
2219               xmu = mu(i, j)
2220               CALL PUSHCONTROL1B(0)
2221             ELSE
2222               CALL PUSHCONTROL1B(1)
2223             END IF
2224           END DO
2225           CALL PUSHINTEGER4(j - 1)
2226           CALL PUSHINTEGER4(ad_from1)
2227         END DO
2228       END DO
2229       CALL PUSHCONTROL1B(0)
2230     ELSE
2231       CALL PUSHCONTROL1B(1)
2232     END IF
2233     IF (ibe - itf .LT. spec_zone) THEN
2234       IF (its .LT. ibe - spec_zone + 1) THEN
2235         max5 = ibe - spec_zone + 1
2236       ELSE
2237         max5 = its
2238       END IF
2239 ! X-end boundary
2240       DO i=max5,itf
2241         CALL PUSHINTEGER4(b_dist)
2242         b_dist = ibe - i
2243         DO k=kts,ktf
2244           IF (jts .LT. b_dist + jbs + 1) THEN
2245             max6 = b_dist + jbs + 1
2246           ELSE
2247             max6 = jts
2248           END IF
2249           IF (jtf .GT. jbe - b_dist - 1) THEN
2250             min6 = jbe - b_dist - 1
2251           ELSE
2252             min6 = jtf
2253           END IF
2254           ad_from2 = max6
2255           DO j=ad_from2,min6
2256             IF (msfcouple) THEN
2257               CALL PUSHREAL8(xmsf)
2258               xmsf = msf(i, j)
2259               CALL PUSHCONTROL1B(0)
2260             ELSE
2261               CALL PUSHCONTROL1B(1)
2262             END IF
2263             IF (mucouple) THEN
2264               CALL PUSHREAL8(xmu)
2265               xmu = mu(i, j)
2266               CALL PUSHCONTROL1B(0)
2267             ELSE
2268               CALL PUSHCONTROL1B(1)
2269             END IF
2270           END DO
2271           CALL PUSHINTEGER4(j - 1)
2272           CALL PUSHINTEGER4(ad_from2)
2273         END DO
2274       END DO
2275       xmub = 0.0
2276       DO i=itf,max5,-1
2277         DO k=ktf,kts,-1
2278           CALL POPINTEGER4(ad_from2)
2279           CALL POPINTEGER4(ad_to2)
2280           DO j=ad_to2,ad_from2,-1
2281             bfield = field_bdy_xe(j, k, b_dist+1) + dtbc*&
2282 &             field_bdy_tend_xe(j, k, b_dist+1)
2283             tempb2 = xmsf*fieldb(i, k, j)/xmu
2284             bfieldb = tempb2
2285             xmub = xmub - bfield*tempb2/xmu
2286             fieldb(i, k, j) = 0.0
2287             CALL POPCONTROL1B(branch)
2288             IF (branch .EQ. 0) THEN
2289               CALL POPREAL8(xmu)
2290               mub(i, j) = mub(i, j) + xmub
2291               xmub = 0.0
2292             END IF
2293             CALL POPCONTROL1B(branch)
2294             IF (branch .EQ. 0) CALL POPREAL8(xmsf)
2295             field_bdy_xeb(j, k, b_dist+1) = field_bdy_xeb(j, k, b_dist+1&
2296 &             ) + bfieldb
2297             field_bdy_tend_xeb(j, k, b_dist+1) = field_bdy_tend_xeb(j, k&
2298 &             , b_dist+1) + dtbc*bfieldb
2299           END DO
2300         END DO
2301         CALL POPINTEGER4(b_dist)
2302       END DO
2303     ELSE
2304       xmub = 0.0
2305     END IF
2306     CALL POPCONTROL1B(branch)
2307     IF (branch .EQ. 0) THEN
2308       DO i=min2,its,-1
2309         DO k=ktf,kts,-1
2310           CALL POPINTEGER4(ad_from1)
2311           CALL POPINTEGER4(ad_to1)
2312           DO j=ad_to1,ad_from1,-1
2313             bfield = field_bdy_xs(j, k, b_dist+1) + dtbc*&
2314 &             field_bdy_tend_xs(j, k, b_dist+1)
2315             tempb1 = xmsf*fieldb(i, k, j)/xmu
2316             bfieldb = tempb1
2317             xmub = xmub - bfield*tempb1/xmu
2318             fieldb(i, k, j) = 0.0
2319             CALL POPCONTROL1B(branch)
2320             IF (branch .EQ. 0) THEN
2321               CALL POPREAL8(xmu)
2322               mub(i, j) = mub(i, j) + xmub
2323               xmub = 0.0
2324             END IF
2325             CALL POPCONTROL1B(branch)
2326             IF (branch .EQ. 0) CALL POPREAL8(xmsf)
2327             field_bdy_xsb(j, k, b_dist+1) = field_bdy_xsb(j, k, b_dist+1&
2328 &             ) + bfieldb
2329             field_bdy_tend_xsb(j, k, b_dist+1) = field_bdy_tend_xsb(j, k&
2330 &             , b_dist+1) + dtbc*bfieldb
2331           END DO
2332         END DO
2333         CALL POPINTEGER4(b_dist)
2334       END DO
2335     END IF
2336   ELSE
2337     xmub = 0.0
2338   END IF
2339   CALL POPCONTROL1B(branch)
2340   IF (branch .EQ. 0) THEN
2341     DO j=jtf,max2,-1
2342       DO k=ktf,kts,-1
2343         CALL POPINTEGER4(ad_from0)
2344         CALL POPINTEGER4(ad_to0)
2345         DO i=ad_to0,ad_from0,-1
2346           bfield = field_bdy_ye(i, k, b_dist+1) + dtbc*field_bdy_tend_ye&
2347 &           (i, k, b_dist+1)
2348           tempb0 = xmsf*fieldb(i, k, j)/xmu
2349           bfieldb = tempb0
2350           xmub = xmub - bfield*tempb0/xmu
2351           fieldb(i, k, j) = 0.0
2352           CALL POPCONTROL1B(branch)
2353           IF (branch .EQ. 0) THEN
2354             CALL POPREAL8(xmu)
2355             mub(i, j) = mub(i, j) + xmub
2356             xmub = 0.0
2357           END IF
2358           CALL POPCONTROL1B(branch)
2359           IF (branch .EQ. 0) CALL POPREAL8(xmsf)
2360           field_bdy_yeb(i, k, b_dist+1) = field_bdy_yeb(i, k, b_dist+1) &
2361 &           + bfieldb
2362           field_bdy_tend_yeb(i, k, b_dist+1) = field_bdy_tend_yeb(i, k, &
2363 &           b_dist+1) + dtbc*bfieldb
2364         END DO
2365       END DO
2366       CALL POPINTEGER4(b_dist)
2367     END DO
2368   END IF
2369   CALL POPCONTROL1B(branch)
2370   IF (branch .EQ. 0) THEN
2371     DO j=min1,jts,-1
2372       DO k=ktf,kts,-1
2373         CALL POPINTEGER4(ad_from)
2374         CALL POPINTEGER4(ad_to)
2375         DO i=ad_to,ad_from,-1
2376           bfield = field_bdy_ys(i, k, b_dist+1) + dtbc*field_bdy_tend_ys&
2377 &           (i, k, b_dist+1)
2378           tempb = xmsf*fieldb(i, k, j)/xmu
2379           bfieldb = tempb
2380           xmub = xmub - bfield*tempb/xmu
2381           fieldb(i, k, j) = 0.0
2382           CALL POPCONTROL1B(branch)
2383           IF (branch .EQ. 0) THEN
2384             CALL POPREAL8(xmu)
2385             mub(i, j) = mub(i, j) + xmub
2386             xmub = 0.0
2387           END IF
2388           CALL POPCONTROL1B(branch)
2389           IF (branch .EQ. 0) CALL POPREAL8(xmsf)
2390           field_bdy_ysb(i, k, b_dist+1) = field_bdy_ysb(i, k, b_dist+1) &
2391 &           + bfieldb
2392           field_bdy_tend_ysb(i, k, b_dist+1) = field_bdy_tend_ysb(i, k, &
2393 &           b_dist+1) + dtbc*bfieldb
2394         END DO
2395       END DO
2396       CALL POPINTEGER4(b_dist)
2397     END DO
2398   END IF
2399 END SUBROUTINE a_SPEC_BDY_FINAL
2400 !------------------------------------------------------------------------
2401
2402    SUBROUTINE a_zero_grad_bdy (  a_field,                     &
2403                                variable_in, config_flags, &
2404                                spec_zone,                  &
2405                                ids,ide, jds,jde, kds,kde,  & ! domain dims
2406                                ims,ime, jms,jme, kms,kme,  & ! memory dims
2407                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
2408                                its,ite, jts,jte, kts,kte )
2409
2410 !  This subroutine sets zero gradient conditions in the boundary specified region.
2411 !  spec_zone is the width of the outer specified b.c.s that are set here.
2412 !  (JD August 2000)
2413
2414       IMPLICIT NONE
2415
2416       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
2417       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
2418       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
2419       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
2420       INTEGER,      INTENT(IN   )    :: spec_zone
2421       CHARACTER,    INTENT(IN   )    :: variable_in
2422
2423
2424       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: a_field
2425       TYPE( grid_config_rec_type ) config_flags
2426
2427       CHARACTER  :: variable
2428       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf, i_inner, j_inner
2429       INTEGER    :: b_dist, b_limit
2430       LOGICAL    :: periodic_x
2431       REAL       :: a_aux
2432
2433
2434       a_aux = 0.
2435
2436       periodic_x = config_flags%periodic_x
2437
2438       variable = variable_in
2439
2440       IF (variable == 'U') variable = 'u'
2441       IF (variable == 'V') variable = 'v'
2442
2443       ibs = ids
2444       ibe = ide-1
2445       itf = min(ite,ide-1)
2446       jbs = jds
2447       jbe = jde-1
2448       jtf = min(jte,jde-1)
2449       ktf = kde-1
2450       IF (variable == 'u') ibe = ide
2451       IF (variable == 'u') itf = min(ite,ide)
2452       IF (variable == 'v') jbe = jde
2453       IF (variable == 'v') jtf = min(jte,jde)
2454       IF (variable == 'w') ktf = kde
2455
2456     IF(.NOT.periodic_x)THEN
2457
2458       IF (ibe - itf .lt. spec_zone) THEN
2459 ! X-end boundary
2460         DO i = itf, max(its,ibe-spec_zone+1), -1
2461           b_dist = ibe - i
2462           DO k = kts, ktf
2463             DO j = min(jtf,jbe-b_dist-1), max(jts,b_dist+jbs+1), -1
2464               j_inner = max(j,jbs+spec_zone)
2465               j_inner = min(j_inner,jbe-spec_zone)
2466               a_aux = a_aux + a_field(i,k,j)
2467               a_field(i,k,j) = 0.
2468               a_field(ibe-spec_zone,k,j_inner) = a_field(ibe-spec_zone,k,j_inner) + a_aux
2469               a_aux = 0.
2470             ENDDO
2471           ENDDO
2472         ENDDO
2473       ENDIF
2474
2475       IF (its - ibs .lt. spec_zone) THEN
2476 ! X-start boundary
2477         DO i = min(itf,ibs+spec_zone-1), its, -1
2478           b_dist = i - ibs
2479           DO k = kts, ktf
2480             DO j = min(jtf,jbe-b_dist-1), max(jts,b_dist+jbs+1), -1
2481               j_inner = max(j,jbs+spec_zone)
2482               j_inner = min(j_inner,jbe-spec_zone)
2483               a_aux = a_aux + a_field(i,k,j)
2484               a_field(i,k,j) = 0.
2485               a_field(ibs+spec_zone,k,j_inner) = a_field(ibs+spec_zone,k,j_inner) + a_aux
2486               a_aux = 0.
2487             ENDDO
2488           ENDDO
2489         ENDDO
2490       ENDIF
2491
2492     ENDIF
2493
2494       IF (jbe - jtf .lt. spec_zone) THEN
2495 ! Y-end boundary
2496         DO j = jtf, max(jts,jbe-spec_zone+1), -1
2497           b_dist = jbe - j
2498           b_limit = b_dist
2499           IF(periodic_x)b_limit = 0
2500           DO k = kts, ktf
2501             DO i = min(itf,ibe-b_limit), max(its,b_limit+ibs), -1
2502               i_inner = max(i,ibs+spec_zone)
2503               i_inner = min(i_inner,ibe-spec_zone)
2504               IF(periodic_x)i_inner = i
2505               a_aux = a_aux + a_field(i,k,j)
2506               a_field(i,k,j) = 0.
2507               a_field(i_inner,k,jbe-spec_zone) = a_field(i_inner,k,jbe-spec_zone) + a_aux
2508               a_aux = 0.
2509             ENDDO
2510           ENDDO
2511         ENDDO
2512       ENDIF
2513
2514       IF (jts - jbs .lt. spec_zone) THEN
2515 ! Y-start boundary
2516         DO j = min(jtf,jbs+spec_zone-1), jts, -1
2517           b_dist = j - jbs
2518           b_limit = b_dist
2519           IF(periodic_x)b_limit = 0
2520           DO k = kts, ktf
2521             DO i = min(itf,ibe-b_limit), max(its,b_limit+ibs), -1
2522               i_inner = max(i,ibs+spec_zone)
2523               i_inner = min(i_inner,ibe-spec_zone)
2524               IF(periodic_x)i_inner = i
2525               a_aux = a_aux + a_field(i,k,j)
2526               a_field(i,k,j) = 0.
2527               a_field(i_inner,k,jbs+spec_zone) = a_field(i_inner,k,jbs+spec_zone) + a_aux
2528               a_aux = 0.
2529             ENDDO
2530           ENDDO
2531         ENDDO
2532       ENDIF
2533
2534    END SUBROUTINE a_zero_grad_bdy
2535
2536 !------------------------------------------------------------------------
2537
2538    SUBROUTINE a_couple_bdy ( field, a_field,  &
2539                              variable_in, config_flags, &
2540                              spec_zone,       &
2541                              mu, a_mu, msf,   &
2542                              ids,ide, jds,jde, kds,kde,  & ! domain dims
2543                              ims,ime, jms,jme, kms,kme,  & ! memory dims
2544                              ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
2545                              its,ite, jts,jte, kts,kte )
2546
2547 !  This subroutine adds the tendencies in the boundary specified region.
2548 !  spec_zone is the width of the outer specified b.c.s that are set here.
2549 !  (JD August 2000)
2550
2551       IMPLICIT NONE
2552
2553       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
2554       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
2555       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
2556       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
2557       INTEGER,      INTENT(IN   )    :: spec_zone
2558       CHARACTER,    INTENT(IN   )    :: variable_in
2559       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(INOUT) :: a_mu
2560       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: mu
2561       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: msf
2562       REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: a_field
2563       REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: field
2564       TYPE( grid_config_rec_type ) config_flags
2565
2566       CHARACTER  :: variable
2567       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf
2568       INTEGER    :: b_dist, b_limit
2569       LOGICAL    :: periodic_x
2570
2571       periodic_x = config_flags%periodic_x
2572
2573       variable = variable_in
2574
2575       IF (variable == 'U') variable = 'u'
2576       IF (variable == 'V') variable = 'v'
2577       IF (variable == 'T') variable = 't'
2578       IF (variable == 'H') variable = 'h'
2579       IF (variable == 'W') variable = 'w'
2580
2581       ibs = ids
2582       ibe = ide-1
2583       itf = min(ite,ide-1)
2584       jbs = jds
2585       jbe = jde-1
2586       jtf = min(jte,jde-1)
2587       ktf = kde-1
2588       IF (variable == 'u') ibe = ide
2589       IF (variable == 'u') itf = min(ite,ide)
2590       IF (variable == 'v') jbe = jde
2591       IF (variable == 'v') jtf = min(jte,jde)
2592       IF (variable == 'h') ktf = kte
2593       IF (variable == 'w') ktf = kte
2594
2595     IF(.NOT.periodic_x)THEN
2596       IF (ibe - itf .lt. spec_zone) THEN
2597 ! X-end boundary
2598         DO i = max(its,ibe-spec_zone+1), itf
2599           b_dist = ibe - i
2600           DO k = kts, ktf
2601             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2602               if (variable == 't' .or. variable == 'h') then
2603                  a_mu(i,j) = a_mu(i,j) + field(i,k,j)*a_field(i,k,j)
2604                  a_field(i,k,j) = a_field(i,k,j)*mu(i,j)
2605               else
2606                  a_mu(i,j) = a_mu(i,j) + field(i,k,j)*a_field(i,k,j)/msf(i,j)
2607                  a_field(i,k,j) = a_field(i,k,j)*mu(i,j)/msf(i,j)
2608               end if
2609             ENDDO
2610           ENDDO
2611         ENDDO
2612       ENDIF
2613
2614        IF (its - ibs .lt. spec_zone) THEN
2615 ! X-start boundary
2616         DO i = its, min(itf,ibs+spec_zone-1)
2617           b_dist = i - ibs
2618           DO k = kts, ktf
2619             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2620               if (variable == 't' .or. variable == 'h') then
2621                  a_mu(i,j) = a_mu(i,j) + field(i,k,j)*a_field(i,k,j)
2622                  a_field(i,k,j) = a_field(i,k,j)*mu(i,j)
2623               else
2624                  a_mu(i,j) = a_mu(i,j) + field(i,k,j)*a_field(i,k,j)/msf(i,j)
2625                  a_field(i,k,j) = a_field(i,k,j)*mu(i,j)/msf(i,j)
2626               end if
2627             ENDDO
2628           ENDDO
2629         ENDDO
2630       ENDIF
2631     ENDIF
2632
2633       IF (jbe - jtf .lt. spec_zone) THEN
2634 ! Y-end boundary
2635         DO j = max(jts,jbe-spec_zone+1), jtf
2636           b_dist = jbe - j
2637           b_limit = b_dist
2638           IF(periodic_x)b_limit = 0
2639           DO k = kts, ktf
2640             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2641               if (variable == 't' .or. variable == 'h') then
2642                  a_mu(i,j) = a_mu(i,j) + field(i,k,j)*a_field(i,k,j)
2643                  a_field(i,k,j) = a_field(i,k,j)*mu(i,j)
2644               else
2645                  a_mu(i,j) = a_mu(i,j) + field(i,k,j)*a_field(i,k,j)/msf(i,j)
2646                  a_field(i,k,j) = a_field(i,k,j)*mu(i,j)/msf(i,j)
2647               end if
2648             ENDDO
2649           ENDDO
2650         ENDDO
2651       ENDIF
2652
2653       IF (jts - jbs .lt. spec_zone) THEN
2654 ! Y-start boundary
2655         DO j = jts, min(jtf,jbs+spec_zone-1)
2656           b_dist = j - jbs
2657           b_limit = b_dist
2658           IF(periodic_x)b_limit = 0
2659           DO k = kts, ktf
2660             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2661               if (variable == 't' .or. variable == 'h') then
2662                  a_mu(i,j) = a_mu(i,j) + field(i,k,j)*a_field(i,k,j)
2663                  a_field(i,k,j) = a_field(i,k,j)*mu(i,j)
2664               else
2665                  a_mu(i,j) = a_mu(i,j) + field(i,k,j)*a_field(i,k,j)/msf(i,j)
2666                  a_field(i,k,j) = a_field(i,k,j)*mu(i,j)/msf(i,j)
2667               end if
2668             ENDDO
2669           ENDDO
2670         ENDDO
2671       ENDIF
2672
2673    END SUBROUTINE a_couple_bdy
2674 !------------------------------------------------------------------------
2675
2676    SUBROUTINE a_uncouple_bdy(  field, a_field,  &
2677                                variable_in, config_flags, &
2678                                spec_zone,       &
2679                                mu, a_mu, msf,   &
2680                                ids,ide, jds,jde, kds,kde,  & ! domain dims
2681                                ims,ime, jms,jme, kms,kme,  & ! memory dims
2682                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
2683                                its,ite, jts,jte, kts,kte )
2684
2685 !  This subroutine adds the tendencies in the boundary specified region.
2686 !  spec_zone is the width of the outer specified b.c.s that are set here.
2687 !  (JD August 2000)
2688
2689       IMPLICIT NONE
2690
2691       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
2692       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
2693       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
2694       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
2695       INTEGER,      INTENT(IN   )    :: spec_zone
2696       CHARACTER,    INTENT(IN   )    :: variable_in
2697       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(INOUT) :: a_mu
2698       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: mu
2699       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: msf
2700       REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: a_field
2701       REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: field
2702       TYPE( grid_config_rec_type ) config_flags
2703
2704       CHARACTER  :: variable
2705       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf
2706       INTEGER    :: b_dist, b_limit
2707       LOGICAL    :: periodic_x
2708
2709       periodic_x = config_flags%periodic_x
2710
2711       variable = variable_in
2712
2713       IF (variable == 'U') variable = 'u'
2714       IF (variable == 'V') variable = 'v'
2715       IF (variable == 'T') variable = 't'
2716       IF (variable == 'H') variable = 'h'
2717       IF (variable == 'W') variable = 'w'
2718
2719       ibs = ids
2720       ibe = ide-1
2721       itf = min(ite,ide-1)
2722       jbs = jds
2723       jbe = jde-1
2724       jtf = min(jte,jde-1)
2725       ktf = kde-1
2726       IF (variable == 'u') ibe = ide
2727       IF (variable == 'u') itf = min(ite,ide)
2728       IF (variable == 'v') jbe = jde
2729       IF (variable == 'v') jtf = min(jte,jde)
2730       IF (variable == 'h') ktf = kte
2731       IF (variable == 'w') ktf = kte
2732
2733     IF(.NOT.periodic_x)THEN
2734       IF (ibe - itf .lt. spec_zone) THEN
2735 ! X-end boundary
2736         DO i = max(its,ibe-spec_zone+1), itf
2737           b_dist = ibe - i
2738           DO k = kts, ktf
2739             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2740               if (variable == 't' .or. variable == 'h') then
2741                  a_mu(i,j) = a_mu(i,j) &
2742                            - field(i,k,j)/(mu(i,j)*mu(i,j)) * a_field(i,k,j)
2743                  a_field(i,k,j) = a_field(i,k,j)/mu(i,j)
2744               else
2745                  a_mu(i,j) = a_mu(i,j) &
2746                            - field(i,k,j)/(mu(i,j)*mu(i,j))*msf(i,j) * a_field(i,k,j)
2747                  a_field(i,k,j) = a_field(i,k,j)/mu(i,j)*msf(i,j)
2748               end if
2749             ENDDO
2750           ENDDO
2751         ENDDO
2752       ENDIF
2753
2754       IF (its - ibs .lt. spec_zone) THEN
2755 ! X-start boundary
2756         DO i = its, min(itf,ibs+spec_zone-1)
2757           b_dist = i - ibs
2758           DO k = kts, ktf
2759             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2760               if (variable == 't' .or. variable == 'h') then
2761                  a_mu(i,j) = a_mu(i,j) &
2762                            - field(i,k,j)/(mu(i,j)*mu(i,j)) * a_field(i,k,j)
2763                  a_field(i,k,j) = a_field(i,k,j)/mu(i,j)
2764               else
2765                  a_mu(i,j) = a_mu(i,j) &
2766                            - field(i,k,j)/(mu(i,j)*mu(i,j))*msf(i,j) * a_field(i,k,j)
2767                  a_field(i,k,j) = a_field(i,k,j)/mu(i,j)*msf(i,j)
2768               end if
2769             ENDDO
2770           ENDDO
2771         ENDDO
2772       ENDIF
2773
2774     ENDIF
2775
2776       IF (jbe - jtf .lt. spec_zone) THEN
2777 ! Y-end boundary
2778         DO j = max(jts,jbe-spec_zone+1), jtf
2779           b_dist = jbe - j
2780           b_limit = b_dist
2781           IF(periodic_x)b_limit = 0
2782           DO k = kts, ktf
2783             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2784               if (variable == 't' .or. variable == 'h') then
2785                  a_mu(i,j) = a_mu(i,j) &
2786                            - field(i,k,j)/(mu(i,j)*mu(i,j)) * a_field(i,k,j)
2787                  a_field(i,k,j) = a_field(i,k,j)/mu(i,j)
2788               else
2789                  a_mu(i,j) = a_mu(i,j) &
2790                            - field(i,k,j)/(mu(i,j)*mu(i,j))*msf(i,j) * a_field(i,k,j)
2791                  a_field(i,k,j) = a_field(i,k,j)/mu(i,j)*msf(i,j)
2792               end if
2793             ENDDO
2794           ENDDO
2795         ENDDO
2796       ENDIF
2797
2798       IF (jts - jbs .lt. spec_zone) THEN
2799 ! Y-start boundary
2800         DO j = jts, min(jtf,jbs+spec_zone-1)
2801           b_dist = j - jbs
2802           b_limit = b_dist
2803           IF(periodic_x)b_limit = 0
2804           DO k = kts, ktf
2805             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2806               if (variable == 't' .or. variable == 'h') then
2807                  a_mu(i,j) = a_mu(i,j) &
2808                            - field(i,k,j)/(mu(i,j)*mu(i,j)) * a_field(i,k,j)
2809                  a_field(i,k,j) = a_field(i,k,j)/mu(i,j)
2810               else
2811                  a_mu(i,j) = a_mu(i,j) &
2812                            - field(i,k,j)/(mu(i,j)*mu(i,j))*msf(i,j) * a_field(i,k,j)
2813                  a_field(i,k,j) = a_field(i,k,j)/mu(i,j)*msf(i,j)
2814               end if
2815             ENDDO
2816           ENDDO
2817         ENDDO
2818       ENDIF
2819
2820    END SUBROUTINE a_uncouple_bdy
2821 !------------------------------------------------------------------------
2822
2823    SUBROUTINE a_flow_dep_bdy ( a_field,            &
2824                                u, v, config_flags, &
2825                                spec_zone,                  &
2826                                ids,ide, jds,jde, kds,kde,  & ! domain dims
2827                                ims,ime, jms,jme, kms,kme,  & ! memory dims
2828                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
2829                                its,ite, jts,jte, kts,kte )
2830
2831       IMPLICIT NONE
2832
2833       INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde
2834       INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme
2835       INTEGER, INTENT(IN) :: ips,ipe, jps,jpe, kps,kpe
2836       INTEGER, INTENT(IN) :: its,ite, jts,jte, kts,kte
2837       INTEGER, INTENT(IN) :: spec_zone
2838
2839
2840       REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: a_field
2841       REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN   ) :: u
2842       REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN   ) :: v
2843       TYPE(grid_config_rec_type),INTENT(IN) :: config_flags
2844
2845       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf, i_inner, j_inner
2846       INTEGER    :: b_dist, b_limit
2847       LOGICAL    :: periodic_x
2848       REAL       :: a_aux
2849
2850
2851       a_aux = 0.0
2852
2853       periodic_x = config_flags%periodic_x
2854
2855       ibs = ids
2856       ibe = ide-1
2857       itf = min(ite,ide-1)
2858       jbs = jds
2859       jbe = jde-1
2860       jtf = min(jte,jde-1)
2861       ktf = kde-1
2862
2863     IF(.NOT.periodic_x)THEN
2864       IF (ibe - itf .lt. spec_zone) THEN
2865 ! X-end boundary
2866         DO i = max(its,ibe-spec_zone+1), itf
2867           b_dist = ibe - i
2868           DO k = kts, ktf
2869             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2870               j_inner = max(j,jbs+spec_zone)
2871               j_inner = min(j_inner,jbe-spec_zone)
2872               IF(u(i+1,k,j) .gt. 0.)THEN
2873                 a_aux = a_aux + a_field(i,k,j)
2874                 a_field(i,k,j) = 0.
2875                 a_field(ibe-spec_zone,k,j_inner) = a_field(ibe-spec_zone,k,j_inner) + a_aux
2876                 a_aux = 0.
2877               ELSE
2878                 a_field(i,k,j) = 0.
2879               ENDIF
2880             ENDDO
2881           ENDDO
2882         ENDDO
2883       ENDIF
2884
2885       IF (its - ibs .lt. spec_zone) THEN
2886 ! X-start boundary
2887         DO i = its, min(itf,ibs+spec_zone-1)
2888           b_dist = i - ibs
2889           DO k = kts, ktf
2890             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2891               j_inner = max(j,jbs+spec_zone)
2892               j_inner = min(j_inner,jbe-spec_zone)
2893               IF(u(i,k,j) .lt. 0.)THEN
2894                 a_aux = a_aux + a_field(i,k,j)
2895                 a_field(i,k,j) = 0.
2896                 a_field(ibs+spec_zone,k,j_inner) = a_field(ibs+spec_zone,k,j_inner) + a_aux
2897                 a_aux = 0.
2898               ELSE
2899                 a_field(i,k,j) = 0.
2900               ENDIF
2901             ENDDO
2902           ENDDO
2903         ENDDO
2904       ENDIF
2905
2906     ENDIF
2907
2908       IF (jbe - jtf .lt. spec_zone) THEN
2909 ! Y-end boundary
2910         DO j = max(jts,jbe-spec_zone+1), jtf
2911           b_dist = jbe - j
2912           b_limit = b_dist
2913           IF(periodic_x)b_limit = 0
2914           DO k = kts, ktf
2915             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2916               i_inner = max(i,ibs+spec_zone)
2917               i_inner = min(i_inner,ibe-spec_zone)
2918               IF(periodic_x)i_inner = i
2919               IF(v(i,k,j+1) .gt. 0.)THEN
2920                 a_aux = a_aux + a_field(i,k,j)
2921                 a_field(i,k,j) = 0.
2922                 a_field(i_inner,k,jbe-spec_zone) = a_field(i_inner,k,jbe-spec_zone) + a_aux
2923                 a_aux = 0.
2924               ELSE
2925                 a_field(i,k,j) = 0.
2926               ENDIF
2927             ENDDO
2928           ENDDO
2929         ENDDO
2930       ENDIF
2931
2932       IF (jts - jbs .lt. spec_zone) THEN
2933 ! Y-start boundary
2934         DO j = jts, min(jtf,jbs+spec_zone-1)
2935           b_dist = j - jbs
2936           b_limit = b_dist
2937           IF(periodic_x)b_limit = 0
2938           DO k = kts, ktf
2939             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2940               i_inner = max(i,ibs+spec_zone)
2941               i_inner = min(i_inner,ibe-spec_zone)
2942               IF(periodic_x)i_inner = i
2943               IF(v(i,k,j) .lt. 0.)THEN
2944                 a_aux = a_aux + a_field(i,k,j)
2945                 a_field(i,k,j) = 0.
2946                 a_field(i_inner,k,jbs+spec_zone) = a_field(i_inner,k,jbs+spec_zone) + a_aux
2947                 a_aux = 0.
2948               ELSE
2949                 a_field(i,k,j) = 0.
2950               ENDIF
2951             ENDDO
2952           ENDDO
2953         ENDDO
2954       ENDIF
2955
2956    END SUBROUTINE a_flow_dep_bdy
2957
2958 !        Generated by TAPENADE     (INRIA, Tropics team)
2959 !  Tapenade 3.6 (r4165) - 21 sep 2011 20:54
2960 !
2961 !  Differentiation of flow_dep_bdy_qnn in reverse (adjoint) mode:
2962 !   gradient     of useful results: field
2963 !   with respect to varying inputs: field
2964 !   RW status of diff variables: field:in-out
2965 ! domain dims
2966 ! memory dims
2967 ! patch  dims
2968 SUBROUTINE A_FLOW_DEP_BDY_QNN(field, fieldb, u, v, config_flags, &
2969 &  spec_zone, ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme&
2970 &  , ips, ipe, jps, jpe, kps, kpe, its, ite, jts, jte, kts, kte)
2971   IMPLICIT NONE
2972   INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde
2973   INTEGER, INTENT(IN) :: ims, ime, jms, jme, kms, kme
2974   INTEGER, INTENT(IN) :: ips, ipe, jps, jpe, kps, kpe
2975   INTEGER, INTENT(IN) :: its, ite, jts, jte, kts, kte
2976   INTEGER, INTENT(IN) :: spec_zone
2977   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: field
2978   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: fieldb
2979   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u
2980   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: v
2981   TYPE(GRID_CONFIG_REC_TYPE) :: config_flags
2982   INTEGER :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf, i_inner, &
2983 &  j_inner
2984   INTEGER :: b_dist, b_limit
2985   LOGICAL :: periodic_x
2986   REAL :: tmp
2987   REAL :: tmp0
2988   REAL :: tmp1
2989   REAL :: tmp2
2990   INTEGER :: branch
2991   INTEGER :: ad_from
2992   INTEGER :: ad_to
2993   INTEGER :: ad_from0
2994   INTEGER :: ad_to0
2995   INTEGER :: ad_from1
2996   INTEGER :: ad_to1
2997   INTEGER :: ad_from2
2998   INTEGER :: ad_to2
2999   INTEGER :: min6
3000   INTEGER :: min5
3001   INTEGER :: min4
3002   INTEGER :: min3
3003   INTEGER :: min2
3004   INTEGER :: min1
3005   INTRINSIC MAX
3006   REAL :: tmpb
3007   REAL :: tmp0b
3008   REAL :: tmp2b
3009   INTRINSIC MIN
3010   INTEGER :: max6
3011   INTEGER :: max5
3012   INTEGER :: max4
3013   INTEGER :: max3
3014   INTEGER :: max2
3015   INTEGER :: max1
3016   REAL :: tmp1b
3017   periodic_x = config_flags%periodic_x
3018   ibs = ids
3019   ibe = ide - 1
3020   IF (ite .GT. ide - 1) THEN
3021     itf = ide - 1
3022   ELSE
3023     itf = ite
3024   END IF
3025   jbs = jds
3026   jbe = jde - 1
3027   IF (jte .GT. jde - 1) THEN
3028     jtf = jde - 1
3029   ELSE
3030     jtf = jte
3031   END IF
3032   ktf = kde - 1
3033   IF (jts - jbs .LT. spec_zone) THEN
3034     IF (jtf .GT. jbs + spec_zone - 1) THEN
3035       min1 = jbs + spec_zone - 1
3036     ELSE
3037       min1 = jtf
3038     END IF
3039 ! Y-start boundary
3040     DO j=jts,min1
3041       b_dist = j - jbs
3042       b_limit = b_dist
3043       IF (periodic_x) b_limit = 0
3044       DO k=kts,ktf
3045         IF (its .LT. b_limit + ibs) THEN
3046           max1 = b_limit + ibs
3047         ELSE
3048           max1 = its
3049         END IF
3050         IF (itf .GT. ibe - b_limit) THEN
3051           min3 = ibe - b_limit
3052         ELSE
3053           min3 = itf
3054         END IF
3055         ad_from = max1
3056         DO i=ad_from,min3
3057           IF (i .LT. ibs + spec_zone) THEN
3058             CALL PUSHINTEGER4(i_inner)
3059             i_inner = ibs + spec_zone
3060             CALL PUSHCONTROL1B(0)
3061           ELSE
3062             CALL PUSHINTEGER4(i_inner)
3063             i_inner = i
3064             CALL PUSHCONTROL1B(1)
3065           END IF
3066           IF (i_inner .GT. ibe - spec_zone) THEN
3067             i_inner = ibe - spec_zone
3068           ELSE
3069             i_inner = i_inner
3070           END IF
3071           IF (periodic_x) i_inner = i
3072           IF (v(i, k, j) .LT. 0.) THEN
3073             CALL PUSHCONTROL1B(1)
3074           ELSE
3075             CALL PUSHCONTROL1B(0)
3076           END IF
3077         END DO
3078         CALL PUSHINTEGER4(i - 1)
3079         CALL PUSHINTEGER4(ad_from)
3080       END DO
3081     END DO
3082     CALL PUSHCONTROL1B(0)
3083   ELSE
3084     CALL PUSHCONTROL1B(1)
3085   END IF
3086   IF (jbe - jtf .LT. spec_zone) THEN
3087     IF (jts .LT. jbe - spec_zone + 1) THEN
3088       max2 = jbe - spec_zone + 1
3089     ELSE
3090       max2 = jts
3091     END IF
3092 ! Y-end boundary
3093     DO j=max2,jtf
3094       b_dist = jbe - j
3095       b_limit = b_dist
3096       IF (periodic_x) b_limit = 0
3097       DO k=kts,ktf
3098         IF (its .LT. b_limit + ibs) THEN
3099           max3 = b_limit + ibs
3100         ELSE
3101           max3 = its
3102         END IF
3103         IF (itf .GT. ibe - b_limit) THEN
3104           min4 = ibe - b_limit
3105         ELSE
3106           min4 = itf
3107         END IF
3108         ad_from0 = max3
3109         DO i=ad_from0,min4
3110           IF (i .LT. ibs + spec_zone) THEN
3111             CALL PUSHINTEGER4(i_inner)
3112             i_inner = ibs + spec_zone
3113             CALL PUSHCONTROL1B(0)
3114           ELSE
3115             CALL PUSHINTEGER4(i_inner)
3116             i_inner = i
3117             CALL PUSHCONTROL1B(1)
3118           END IF
3119           IF (i_inner .GT. ibe - spec_zone) THEN
3120             i_inner = ibe - spec_zone
3121           ELSE
3122             i_inner = i_inner
3123           END IF
3124           IF (periodic_x) i_inner = i
3125           IF (v(i, k, j+1) .GT. 0.) THEN
3126             CALL PUSHCONTROL1B(1)
3127           ELSE
3128             CALL PUSHCONTROL1B(0)
3129           END IF
3130         END DO
3131         CALL PUSHINTEGER4(i - 1)
3132         CALL PUSHINTEGER4(ad_from0)
3133       END DO
3134     END DO
3135     CALL PUSHCONTROL1B(0)
3136   ELSE
3137     CALL PUSHCONTROL1B(1)
3138   END IF
3139   IF (.NOT.periodic_x) THEN
3140     IF (its - ibs .LT. spec_zone) THEN
3141       IF (itf .GT. ibs + spec_zone - 1) THEN
3142         min2 = ibs + spec_zone - 1
3143       ELSE
3144         min2 = itf
3145       END IF
3146 ! X-start boundary
3147       DO i=its,min2
3148         b_dist = i - ibs
3149         DO k=kts,ktf
3150           IF (jts .LT. b_dist + jbs + 1) THEN
3151             max4 = b_dist + jbs + 1
3152           ELSE
3153             max4 = jts
3154           END IF
3155           IF (jtf .GT. jbe - b_dist - 1) THEN
3156             min5 = jbe - b_dist - 1
3157           ELSE
3158             min5 = jtf
3159           END IF
3160           ad_from1 = max4
3161           DO j=ad_from1,min5
3162             IF (j .LT. jbs + spec_zone) THEN
3163               CALL PUSHINTEGER4(j_inner)
3164               j_inner = jbs + spec_zone
3165               CALL PUSHCONTROL1B(0)
3166             ELSE
3167               CALL PUSHINTEGER4(j_inner)
3168               j_inner = j
3169               CALL PUSHCONTROL1B(1)
3170             END IF
3171             IF (j_inner .GT. jbe - spec_zone) THEN
3172               j_inner = jbe - spec_zone
3173             ELSE
3174               j_inner = j_inner
3175             END IF
3176             IF (u(i, k, j) .LT. 0.) THEN
3177               CALL PUSHCONTROL1B(1)
3178             ELSE
3179               CALL PUSHCONTROL1B(0)
3180             END IF
3181           END DO
3182           CALL PUSHINTEGER4(j - 1)
3183           CALL PUSHINTEGER4(ad_from1)
3184         END DO
3185       END DO
3186       CALL PUSHCONTROL1B(0)
3187     ELSE
3188       CALL PUSHCONTROL1B(1)
3189     END IF
3190     IF (ibe - itf .LT. spec_zone) THEN
3191       IF (its .LT. ibe - spec_zone + 1) THEN
3192         max5 = ibe - spec_zone + 1
3193       ELSE
3194         max5 = its
3195       END IF
3196 ! X-end boundary
3197       DO i=max5,itf
3198         b_dist = ibe - i
3199         DO k=kts,ktf
3200           IF (jts .LT. b_dist + jbs + 1) THEN
3201             max6 = b_dist + jbs + 1
3202           ELSE
3203             max6 = jts
3204           END IF
3205           IF (jtf .GT. jbe - b_dist - 1) THEN
3206             min6 = jbe - b_dist - 1
3207           ELSE
3208             min6 = jtf
3209           END IF
3210           ad_from2 = max6
3211           DO j=ad_from2,min6
3212             IF (j .LT. jbs + spec_zone) THEN
3213               CALL PUSHINTEGER4(j_inner)
3214               j_inner = jbs + spec_zone
3215               CALL PUSHCONTROL1B(0)
3216             ELSE
3217               CALL PUSHINTEGER4(j_inner)
3218               j_inner = j
3219               CALL PUSHCONTROL1B(1)
3220             END IF
3221             IF (j_inner .GT. jbe - spec_zone) THEN
3222               j_inner = jbe - spec_zone
3223             ELSE
3224               j_inner = j_inner
3225             END IF
3226             IF (u(i+1, k, j) .GT. 0.) THEN
3227               CALL PUSHCONTROL1B(1)
3228             ELSE
3229               CALL PUSHCONTROL1B(0)
3230             END IF
3231           END DO
3232           CALL PUSHINTEGER4(j - 1)
3233           CALL PUSHINTEGER4(ad_from2)
3234         END DO
3235       END DO
3236       DO i=itf,max5,-1
3237         DO k=ktf,kts,-1
3238           CALL POPINTEGER4(ad_from2)
3239           CALL POPINTEGER4(ad_to2)
3240           DO j=ad_to2,ad_from2,-1
3241             CALL POPCONTROL1B(branch)
3242             IF (branch .EQ. 0) THEN
3243               fieldb(i, k, j) = 0.0
3244             ELSE
3245               tmp2b = fieldb(i, k, j)
3246               fieldb(i, k, j) = 0.0
3247               fieldb(ibe-spec_zone, k, j_inner) = fieldb(ibe-spec_zone, &
3248 &                k, j_inner) + tmp2b
3249             END IF
3250             CALL POPCONTROL1B(branch)
3251             IF (branch .EQ. 0) THEN
3252               CALL POPINTEGER4(j_inner)
3253             ELSE
3254               CALL POPINTEGER4(j_inner)
3255             END IF
3256           END DO
3257         END DO
3258       END DO
3259     END IF
3260     CALL POPCONTROL1B(branch)
3261     IF (branch .EQ. 0) THEN
3262       DO i=min2,its,-1
3263         DO k=ktf,kts,-1
3264           CALL POPINTEGER4(ad_from1)
3265           CALL POPINTEGER4(ad_to1)
3266           DO j=ad_to1,ad_from1,-1
3267             CALL POPCONTROL1B(branch)
3268             IF (branch .EQ. 0) THEN
3269               fieldb(i, k, j) = 0.0
3270             ELSE
3271               tmp1b = fieldb(i, k, j)
3272               fieldb(i, k, j) = 0.0
3273               fieldb(ibs+spec_zone, k, j_inner) = fieldb(ibs+spec_zone, &
3274 &                k, j_inner) + tmp1b
3275             END IF
3276             CALL POPCONTROL1B(branch)
3277             IF (branch .EQ. 0) THEN
3278               CALL POPINTEGER4(j_inner)
3279             ELSE
3280               CALL POPINTEGER4(j_inner)
3281             END IF
3282           END DO
3283         END DO
3284       END DO
3285     END IF
3286   END IF
3287   CALL POPCONTROL1B(branch)
3288   IF (branch .EQ. 0) THEN
3289     DO j=jtf,max2,-1
3290       DO k=ktf,kts,-1
3291         CALL POPINTEGER4(ad_from0)
3292         CALL POPINTEGER4(ad_to0)
3293         DO i=ad_to0,ad_from0,-1
3294           CALL POPCONTROL1B(branch)
3295           IF (branch .EQ. 0) THEN
3296             fieldb(i, k, j) = 0.0
3297           ELSE
3298             tmp0b = fieldb(i, k, j)
3299             fieldb(i, k, j) = 0.0
3300             fieldb(i_inner, k, jbe-spec_zone) = fieldb(i_inner, k, jbe-&
3301 &              spec_zone) + tmp0b
3302           END IF
3303           CALL POPCONTROL1B(branch)
3304           IF (branch .EQ. 0) THEN
3305             CALL POPINTEGER4(i_inner)
3306           ELSE
3307             CALL POPINTEGER4(i_inner)
3308           END IF
3309         END DO
3310       END DO
3311     END DO
3312   END IF
3313   CALL POPCONTROL1B(branch)
3314   IF (branch .EQ. 0) THEN
3315     DO j=min1,jts,-1
3316       DO k=ktf,kts,-1
3317         CALL POPINTEGER4(ad_from)
3318         CALL POPINTEGER4(ad_to)
3319         DO i=ad_to,ad_from,-1
3320           CALL POPCONTROL1B(branch)
3321           IF (branch .EQ. 0) THEN
3322             fieldb(i, k, j) = 0.0
3323           ELSE
3324             tmpb = fieldb(i, k, j)
3325             fieldb(i, k, j) = 0.0
3326             fieldb(i_inner, k, jbs+spec_zone) = fieldb(i_inner, k, jbs+&
3327 &              spec_zone) + tmpb
3328           END IF
3329           CALL POPCONTROL1B(branch)
3330           IF (branch .EQ. 0) THEN
3331             CALL POPINTEGER4(i_inner)
3332           ELSE
3333             CALL POPINTEGER4(i_inner)
3334           END IF
3335         END DO
3336       END DO
3337     END DO
3338   END IF
3339 END SUBROUTINE A_FLOW_DEP_BDY_QNN
3340
3341 !---------------------------------------ntu3m-------------------------
3342    SUBROUTINE a_flow_dep_bdy_fixed_inflow(a_field,u,v,config_flags,  &
3343                                  spec_zone,ids,ide,jds,jde,kds,kde,  &
3344                                  ims,ime,jms,jme,kms,kme,ips,ipe,jps,&
3345                                  jpe,kps,kpe,its,ite,jts,jte,kts,kte)
3346 !---------------------------------------------------------------------
3347       IMPLICIT NONE
3348       INTEGER, INTENT(IN) :: ids,ide,jds,jde,kds,kde,ims,ime,jms,jme, &
3349                              kms,kme,ips,ipe,jps,jpe,kps,kpe,its,ite, &
3350                              jts,jte,kts,kte,spec_zone
3351       REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) ::      &
3352                                                 a_field
3353       REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN) :: u,v
3354       TYPE(grid_config_rec_type),INTENT(IN) :: config_flags
3355       INTEGER :: i,j,k,ibs,ibe,jbs,jbe,itf,jtf,ktf,i_inner,j_inner,   &
3356                  b_dist,b_limit
3357       LOGICAL :: periodic_x
3358       REAL :: a_aux
3359
3360       a_aux = 0.0
3361       periodic_x = config_flags%periodic_x
3362       ibs = ids
3363       ibe = ide-1
3364       itf = min(ite,ide-1)
3365       jbs = jds
3366       jbe = jde-1
3367       jtf = min(jte,jde-1)
3368       ktf = kde-1
3369
3370       IF (.NOT.periodic_x) THEN
3371          IF (ibe - itf .lt. spec_zone) THEN
3372 ! X-end boundary
3373             DO i = max(its,ibe-spec_zone+1), itf
3374                b_dist = ibe - i
3375                DO k = kts, ktf
3376                   DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
3377                      j_inner = max(j,jbs+spec_zone)
3378                      j_inner = min(j_inner,jbe-spec_zone)
3379                      IF (u(i+1,k,j) .gt. 0.) THEN
3380                         a_aux = a_aux+a_field(i,k,j)
3381                         a_field(i,k,j) = 0.
3382                         a_field(ibe-spec_zone,k,j_inner) =             &
3383                         a_field(ibe-spec_zone,k,j_inner)+a_aux
3384                         a_aux = 0.
3385                      ELSE
3386                         a_field(i,k,j) = a_field(i,k,j)
3387                      ENDIF
3388                   ENDDO
3389                ENDDO
3390             ENDDO
3391          ENDIF
3392
3393          IF (its - ibs .lt. spec_zone) THEN
3394 ! X-start boundary
3395             DO i = its,min(itf,ibs+spec_zone-1)
3396                b_dist = i - ibs
3397                DO k = kts, ktf
3398                   DO j = max(jts,b_dist+jbs+1),min(jtf,jbe-b_dist-1)
3399                      j_inner = max(j,jbs+spec_zone)
3400                      j_inner = min(j_inner,jbe-spec_zone)
3401                      IF (u(i,k,j) .lt. 0.) THEN
3402                         a_aux = a_aux+a_field(i,k,j)
3403                         a_field(i,k,j) = 0.
3404                         a_field(ibs+spec_zone,k,j_inner) =             &
3405                         a_field(ibs+spec_zone,k,j_inner)+a_aux
3406                         a_aux = 0.
3407                      ELSE
3408                         a_field(i,k,j) = a_field(i,k,j)
3409                      ENDIF
3410                   ENDDO
3411                ENDDO
3412             ENDDO
3413          ENDIF
3414       ENDIF
3415
3416       IF (jbe - jtf .lt. spec_zone) THEN
3417 ! Y-end boundary
3418          DO j = max(jts,jbe-spec_zone+1), jtf
3419             b_dist = jbe - j
3420             b_limit = b_dist
3421             IF (periodic_x) b_limit = 0
3422             DO k = kts, ktf
3423                DO i = max(its,b_limit+ibs),min(itf,ibe-b_limit)
3424                   i_inner = max(i,ibs+spec_zone)
3425                   i_inner = min(i_inner,ibe-spec_zone)
3426                   IF (periodic_x) i_inner = i
3427                   IF (v(i,k,j+1).gt.0.) THEN
3428                      a_aux = a_aux+a_field(i,k,j)
3429                      a_field(i,k,j) = 0.
3430                      a_field(i_inner,k,jbe-spec_zone) =                &
3431                      a_field(i_inner,k,jbe-spec_zone)+a_aux
3432                      a_aux = 0.
3433                   ELSE
3434                      a_field(i,k,j) = a_field(i,k,j)
3435                   ENDIF
3436                ENDDO
3437             ENDDO
3438          ENDDO
3439       ENDIF
3440
3441       IF (jts - jbs .lt. spec_zone) THEN
3442 ! Y-start boundary
3443          DO j = jts, min(jtf,jbs+spec_zone-1)
3444             b_dist = j - jbs
3445             b_limit = b_dist
3446             IF (periodic_x) b_limit = 0
3447             DO k = kts, ktf
3448                DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
3449                   i_inner = max(i,ibs+spec_zone)
3450                   i_inner = min(i_inner,ibe-spec_zone)
3451                   IF (periodic_x) i_inner = i
3452                   IF (v(i,k,j) .lt. 0.) THEN
3453                      a_aux = a_aux+a_field(i,k,j)
3454                      a_field(i,k,j) = 0.
3455                      a_field(i_inner,k,jbs+spec_zone) =                &
3456                      a_field(i_inner,k,jbs+spec_zone)+a_aux
3457                      a_aux = 0.
3458                   ELSE
3459                      a_field(i,k,j) = a_field(i,k,j)
3460                   ENDIF
3461                ENDDO
3462             ENDDO
3463          ENDDO
3464       ENDIF
3465
3466    END SUBROUTINE a_flow_dep_bdy_fixed_inflow
3467 !--------------------------------------------------ntu3m-------------------------
3468
3469 END MODULE a_module_bc