wrftladj/module_small_step_em_ad.F

   1
   2 ! ======================================================================================
   3 ! This file was generated by the version 4.3.6 of ADG on 06/26/2010. The Adjoint Code
   4 ! Generator (ADG) was developed and sponsored by LASG of IAP (1999-2010)
   5 ! The Copyright of the ADG system was declared by Walls at LASG, 1999-2010
   6 ! ======================================================================================
   7 ! Ning Pan, 2010-07-10: Change Adj_ to a_
   8 ! Ning Pan, 2010-07-11: Change the order of arguments
   9
  10 MODULE a_module_small_step_em
  11
  12    USE module_model_constants
  13    USE module_configure
  14
  15 CONTAINS
  16
  17 SUBROUTINE a_advance_all (a_u, ru_tend, a_ru_tend, a_v, rv_tend, a_rv_tend,     &
  18                           a_w, rw_tend, a_rw_tend, a_t, t_tend, a_t_tend,       &
  19                           a_mu, a_mu_tend, a_ph, ph_tend, a_ph_tend,   &
  20                           muu, a_muu, muv, a_muv, mut, a_mut,                   &
  21                           msfuy, msfvx, msfty,           &
  22                           dts,                           &
  23                           config_flags, spec_zone,       &
  24                           ids, ide, jds, jde, kds, kde,  &
  25                           ims, ime, jms, jme, kms, kme,  &
  26                           ips, ipe, jps, jpe, kps, kpe,  &
  27                           its, ite, jts, jte, kts, kte  )
  28
  29       IMPLICIT NONE  ! religion first
  30
  31       TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
  32
  33       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
  34       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
  35       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
  36       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
  37       INTEGER,      INTENT(IN   )    :: spec_zone
  38
  39       REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),  &
  40             INTENT(IN   ) ::                          &
  41                                                   a_u,  &
  42                                                   a_v,  &
  43                                                   a_w,  &
  44                                                   a_t,  &
  45                                                   a_ph
  46
  47       REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
  48             INTENT(INOUT) ::                          &
  49                                              a_ru_tend, &
  50                                              a_rv_tend, &
  51                                              a_rw_tend, &
  52                                               a_t_tend, &
  53                                              a_ph_tend
  54       REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
  55             INTENT(IN   ) ::                          &
  56                                              ru_tend, &
  57                                              rv_tend, &
  58                                              rw_tend, &
  59                                               t_tend, &
  60                                              ph_tend
  61
  62       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: a_mu
  63       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(INOUT) :: a_mu_tend
  64
  65       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(INOUT) :: a_muu,  &
  66                                                                 a_muv,  &
  67                                                                 a_mut
  68       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: muu,  &
  69                                                                 muv,  &
  70                                                                 mut
  71
  72       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: msfuy,  &
  73                                                                 msfvx,  &
  74                                                                 msfty
  75
  76       REAL,                                    INTENT(IN   ) :: dts
  77
  78
  79       INTEGER :: i,j,k, i_start, i_end, j_start, j_end, k_start, k_end
  80       INTEGER :: i_endu, j_endv, k_endw
  81
  82       INTEGER :: i_start_u_tend, i_end_u_tend, j_start_v_tend, j_end_v_tend
  83
  84 !<DESCRIPTION>
  85 !
  86 !  advance_all advances the explicit perturbation horizontal momentum
  87 !  equations (u,v) by adding in the large-timestep tendency along with
  88 !  the small timestep pressure gradient tendency.
  89
  90 !</DESCRIPTION>
  91
  92 !  now, the real work.
  93 !  set the loop bounds taking into account boundary conditions.
  94
  95    i_start = its
  96    i_end   = min(ite,ide-1)
  97    j_start = jts
  98    j_end   = min(jte,jde-1)
  99    k_start = kts
 100    k_end   = kte-1
 101    IF ( .NOT. config_flags%periodic_x )THEN
 102       IF ( config_flags%specified .or. config_flags%nested ) then
 103         i_start = max(its,ids+1)
 104         i_end   = min(ite,ide-2)
 105       ENDIF
 106    ENDIF
 107    IF ( config_flags%specified .or. config_flags%nested ) then
 108       j_start = max(jts,jds+1)
 109       j_end   = min(jte,jde-2)
 110    ENDIF
 111
 112 IF ( config_flags%non_hydrostatic ) THEN
 113   j_loop_w:  DO j = j_start, j_end
 114     DO k=2,k_end+1
 115       DO i=i_start, i_end
 116         a_ph_tend(i,k,j) = a_ph_tend(i,k,j) &
 117                          + dts*msfty(i,j)/mut(i,j) * a_ph(i,k,j)
 118         a_mut(i,j) = a_mut(i,j) &
 119                    - dts*msfty(i,j)*ph_tend(i,k,j)/(mut(i,j)*mut(i,j)) * a_ph(i,k,j)
 120
 121         a_rw_tend(i,k,j) = a_rw_tend(i,k,j) &
 122                          + dts*msfty(i,j)/mut(i,j) * a_w(i,k,j)
 123         a_mut(i,j) = a_mut(i,j) &
 124                    - dts*msfty(i,j)*rw_tend(i,k,j)/(mut(i,j)*mut(i,j)) * a_w(i,k,j)
 125
 126       ENDDO
 127     ENDDO
 128   ENDDO j_loop_w
 129 ENDIF
 130
 131    i_start = its
 132    i_end   = min(ite,ide-1)
 133    j_start = jts
 134    j_end   = min(jte,jde-1)
 135    k_start = kts
 136    k_end   = kte-1
 137    IF ( .NOT. config_flags%periodic_x )THEN
 138      IF ( config_flags%specified .or. config_flags%nested ) then
 139        i_start = max(its,ids+1)
 140        i_end   = min(ite,ide-2)
 141      ENDIF
 142    ENDIF
 143    IF ( config_flags%specified .or. config_flags%nested ) then
 144      j_start = max(jts,jds+1)
 145      j_end   = min(jte,jde-2)
 146    ENDIF
 147
 148    DO j=j_start, j_end
 149      DO k=1,k_end
 150      DO i=i_start, i_end
 151        a_t_tend(i,k,j) = a_t_tend(i,k,j) &
 152                        + dts*msfty(i,j)/mut(i,j) * a_t(i,k,j)
 153        a_mut(i,j) = a_mut(i,j) &
 154                   - dts*msfty(i,j)*t_tend(i,k,j)/(mut(i,j)*mut(i,j)) * a_t(i,k,j)
 155      END DO
 156      END DO
 157    ENDDO
 158
 159    DO j = j_start, j_end
 160    DO i=i_start, i_end
 161      a_mu_tend(i,j) = a_mu_tend(i,j)+dts * a_mu(i,j)
 162    ENDDO
 163    ENDDO
 164
 165    IF( config_flags%nested .or. config_flags%specified ) THEN
 166       i_start = max( its,ids+spec_zone )
 167       i_end   = min( ite,ide-spec_zone-1 )
 168       j_start = max( jts,jds+spec_zone )
 169       j_end   = min( jte,jde-spec_zone-1 )
 170       k_start = kts
 171       k_end   = min( kte, kde-1 )
 172
 173       i_endu = min( ite,ide-spec_zone )
 174       j_endv = min( jte,jde-spec_zone )
 175       k_endw = k_end
 176
 177       IF( config_flags%periodic_x) THEN
 178         i_start = its
 179         i_end   = min(ite,ide-1)
 180         i_endu = ite
 181       ENDIF
 182    ELSE
 183       i_start = its
 184       i_end   = min(ite,ide-1)
 185       j_start = jts
 186       j_end   = min(jte,jde-1)
 187       k_start = kts
 188       k_end   = kte-1
 189
 190       i_endu = ite
 191       j_endv = jte
 192       k_endw = k_end
 193    ENDIF
 194
 195    i_start_u_tend = i_start
 196    i_end_u_tend   = i_endu
 197    j_start_v_tend = j_start
 198    j_end_v_tend   = j_endv
 199
 200    IF ( config_flags%symmetric_xs .and. (its == ids) ) &
 201                   i_start_u_tend = i_start_u_tend+1
 202    IF ( config_flags%symmetric_xe .and. (ite == ide) ) &
 203                   i_end_u_tend = i_end_u_tend-1
 204    IF ( config_flags%symmetric_ys .and. (jts == jds) ) &
 205                   j_start_v_tend = j_start_v_tend+1
 206    IF ( config_flags%symmetric_ye .and. (jte == jde) ) &
 207                   j_end_v_tend = j_end_v_tend-1
 208
 209    v_outer_j_loop: DO j = j_start_v_tend, j_end_v_tend
 210       DO k = k_start, k_end
 211       DO i = i_start, i_end
 212         a_rv_tend(i,k,j) = a_rv_tend(i,k,j) &
 213                          + dts*msfvx(i,j)/muv(i,j) * a_v(i,k,j)
 214         a_muv(i,j) = a_muv(i,j) &
 215                    - dts*msfvx(i,j)*rv_tend(i,k,j)/(muv(i,j)*muv(i,j)) * a_v(i,k,j)
 216       ENDDO
 217       ENDDO
 218    ENDDO v_outer_j_loop
 219
 220    u_outer_j_loop: DO j = j_start, j_end
 221       DO k = k_start, k_end
 222       DO i = i_start_u_tend, i_end_u_tend
 223         a_ru_tend(i,k,j) = a_ru_tend(i,k,j) &
 224                          + dts*msfuy(i,j)/muu(i,j) * a_u(i,k,j)
 225         a_muu(i,j) = a_muu(i,j) &
 226                    - dts*msfuy(i,j)*ru_tend(i,k,j)/(muu(i,j)*muu(i,j)) * a_u(i,k,j)
 227       ENDDO
 228       ENDDO
 229    ENDDO u_outer_j_loop
 230
 231 END SUBROUTINE a_advance_all
 232
 233 SUBROUTINE a_save_ph_mu ( a_ph_1, a_ph_2, a_ph_save, &
 234                           a_mu_1, a_mu_2, a_mu_save, &
 235                           rk_step,                      &
 236                           ids,ide, jds,jde, kds,kde,    &
 237                           ims,ime, jms,jme, kms,kme,    &
 238                           its,ite, jts,jte, kts,kte    )
 239
 240   IMPLICIT NONE  ! religion first
 241
 242 ! declarations for the stuff coming in
 243
 244   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
 245   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
 246   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
 247
 248   INTEGER,      INTENT(IN   )    :: rk_step
 249
 250   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: a_ph_1
 251   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: a_ph_save
 252   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: a_ph_2
 253   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: a_mu_1
 254   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: a_mu_save
 255   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: a_mu_2
 256
 257 ! local variables
 258
 259   INTEGER :: i, j, k
 260   INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
 261   INTEGER :: i_endu, j_endv
 262
 263
 264     i_start = its
 265     i_end   = min(ite,ide-1)
 266     j_start = jts
 267     j_end   = min(jte,jde-1)
 268     k_start = kts
 269     k_end = min(kte,kde-1)
 270
 271     DO j=j_start, j_end
 272     DO i=i_start, i_end
 273       a_mu_1(i,j)=a_mu_1(i,j)+a_mu_2(i,j)
 274       a_mu_2(i,j)=-a_mu_2(i,j)
 275       a_mu_2(i,j)=a_mu_2(i,j) + a_mu_save(i,j)
 276       a_mu_save(i,j)=0.
 277     ENDDO
 278     ENDDO
 279
 280     DO j=j_start, j_end
 281 !    DO k=k_start, min(kde,kte)
 282     DO k=k_start, kde
 283     DO i=i_start, i_end
 284       a_ph_1(i,k,j) = a_ph_1(i,k,j) + a_ph_2(i,k,j)
 285       a_ph_2(i,k,j) = -a_ph_2(i,k,j)
 286       a_ph_2(i,k,j) = a_ph_2(i,k,j) + a_ph_save(i,k,j)
 287       a_ph_save(i,k,j) = 0.
 288     ENDDO
 289     ENDDO
 290     ENDDO
 291
 292     IF ((rk_step == 1) ) THEN
 293
 294       DO j=j_start, j_end
 295       DO i=i_start, i_end
 296         a_mu_2(i,j) = a_mu_2(i,j) + a_mu_1(i,j)
 297         a_mu_1(i,j) = 0.
 298       ENDDO
 299       ENDDO
 300
 301       DO j=j_start, j_end
 302       DO k=k_start, min(kde,kte)
 303       DO i=i_start, i_end
 304         a_ph_2(i,k,j) = a_ph_2(i,k,j) + a_ph_1(i,k,j)
 305         a_ph_1(i,k,j) = 0.
 306       ENDDO
 307       ENDDO
 308       ENDDO
 309
 310     END IF
 311
 312 END SUBROUTINE a_save_ph_mu
 313
 314 !----------------------------------------------------------------------
 315
 316 SUBROUTINE a_restore_ph_mu ( a_ph_2, a_ph_save, &
 317                              a_mu_2, a_mu_save, &
 318                              ids,ide, jds,jde, kds,kde,    &
 319                              ims,ime, jms,jme, kms,kme,    &
 320                              its,ite, jts,jte, kts,kte )
 321
 322   IMPLICIT NONE  ! religion first
 323
 324 ! declarations for the stuff coming in
 325
 326   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
 327   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
 328   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
 329
 330   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: a_ph_save
 331   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: a_ph_2
 332   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: a_mu_2
 333   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: a_mu_save
 334
 335 ! local variables
 336
 337   INTEGER :: i, j, k
 338   INTEGER :: i_start, i_end, j_start, j_end
 339
 340
 341   i_start = its
 342   i_end   = min(ite,ide-1)
 343   j_start = jts
 344   j_end   = min(jte,jde-1)
 345
 346   DO j = j_start, j_end
 347   DO i = i_start, i_end
 348     a_mu_save(i,j) = a_mu_save(i,j) + a_mu_2(i,j)
 349   ENDDO
 350   ENDDO
 351
 352   DO j = j_start, j_end
 353   DO k = kds, kde
 354   DO i = i_start, i_end
 355     a_ph_save(i,k,j) = a_ph_save(i,k,j) + a_ph_2(i,k,j)
 356   ENDDO
 357   ENDDO
 358   ENDDO
 359
 360 END SUBROUTINE a_restore_ph_mu
 361
 362    SUBROUTINE a_small_step_prep(u_1,a_u_1,u_2,a_u_2,v_1,a_v_1,v_2,a_v_2, &
 363    w_1,a_w_1,w_2,a_w_2,t_1,a_t_1,t_2,a_t_2,ph_1,a_ph_1,ph_2,a_ph_2, &
 364    mub,mu_1,a_mu_1,mu_2,a_mu_2,muu,a_muu,muus,a_muus,muv,a_muv,muvs, &
 365    a_muvs,mut,a_mut,muts,a_muts,mudf,a_mudf,u_save,a_u_save,v_save, &
 366    a_v_save,w_save,a_w_save,t_save,a_t_save,ph_save,a_ph_save,mu_save, &
 367    a_mu_save,ww,a_ww,ww_save,a_ww_save,c2a,a_c2a,pb,p,a_p,alt,a_alt, &
 368    msfux,msfuy,msfvx,msfvx_inv,msfvy,msftx,msfty,rdx,rdy,rk_step,ids,ide,jds,jde,kds, &
 369    kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte)
 370
 371
 372 !PART I: DECLARATION OF VARIABLES
 373
 374    IMPLICIT NONE
 375
 376    INTEGER :: K0_ADJ,K1_ADJ,K2_ADJ,K3_ADJ
 377    INTEGER :: ids,ide,jds,jde,kds,kde
 378    INTEGER :: ims,ime,jms,jme,kms,kme
 379    INTEGER :: its,ite,jts,jte,kts,kte
 380    INTEGER :: rk_step
 381    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_u_1,u_1,a_v_1,v_1,a_w_1,w_1, &
 382    a_t_1,t_1,a_ph_1,ph_1
 383    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_u_save,u_save,a_v_save,v_save, &
 384    a_w_save,w_save,a_t_save,t_save,a_ph_save,ph_save
 385    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_u_2,u_2,a_v_2,v_2,a_w_2,w_2, &
 386    a_t_2,t_2,a_ph_2,ph_2
 387    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_c2a,c2a,a_ww_save,ww_save
 388    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: pb,a_p,p,a_alt,alt,a_ww,ww
 389    REAL,DIMENSION(ims:ime,jms:jme) :: a_mu_1,mu_1,a_mu_2,mu_2
 390    REAL,DIMENSION(ims:ime,jms:jme) :: mub,a_muu,muu,a_muv,muv,a_mut,mut,msfux, &
 391    msfuy,msfvx,msfvx_inv,msfvy,msftx,msfty
 392    REAL,DIMENSION(ims:ime,jms:jme) :: a_muus,muus,a_muvs,muvs,a_muts,muts,a_mudf,mudf
 393    REAL,DIMENSION(ims:ime,jms:jme) :: a_mu_save,mu_save
 394    REAL :: rdx,rdy
 395    INTEGER :: i,j,k
 396    INTEGER :: i_start,i_end,j_start,j_end,k_start,k_end
 397    INTEGER :: i_endu,j_endv
 398
 399    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Keep_Lpb4_u_2
 400    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Keep_Lpb5_v_2
 401    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Keep_Lpb6_t_2
 402 !  REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Keep_Lpb7_w_2
 403    REAL :: a_Tmpv1,Tmpv001,a_Tmpv2,Tmpv002,a_Tmpv3,Tmpv003,a_Tmpv4,Tmpv004
 404    REAL,DIMENSION(its:max0(ite,min(ite,ide-1)),kts:max0(min(kte,kde-1),kde)) :: Tmpv300
 405
 406
 407 !PART II: CALCULATIONS OF B. S. TRAJECTORY
 408
 409 !LPB[0]
 410        i_start = its
 411        i_end   = min(ite,ide-1)
 412        j_start = jts
 413        j_end   = min(jte,jde-1)
 414        k_start = kts
 415        k_end = min(kte,kde-1)
 416        i_endu = ite
 417        j_endv = jte
 418
 419 !LPB[1]
 420     IF ((rk_step == 1) ) THEN
 421
 422          DO j=j_start, j_end
 423          DO i=i_start, i_end
 424            mu_1(i,j)=mu_2(i,j)
 425            ww_save(i,kde,j) = 0.
 426            ww_save(i,1,j) = 0.
 427            mudf(i,j) = 0.
 428          ENDDO
 429          ENDDO
 430
 431          DO j=j_start, j_end
 432          DO k=k_start, k_end
 433          DO i=i_start, i_endu
 434            u_1(i,k,j) = u_2(i,k,j)
 435          ENDDO
 436          ENDDO
 437          ENDDO
 438
 439          DO j=j_start, j_endv
 440          DO k=k_start, k_end
 441          DO i=i_start, i_end
 442            v_1(i,k,j) = v_2(i,k,j)
 443          ENDDO
 444          ENDDO
 445          ENDDO
 446
 447          DO j=j_start, j_end
 448          DO k=k_start, k_end
 449          DO i=i_start, i_end
 450            t_1(i,k,j) = t_2(i,k,j)
 451          ENDDO
 452          ENDDO
 453          ENDDO
 454
 455          DO j=j_start, j_end
 456          DO k=k_start, min(kde,kte)
 457          DO i=i_start, i_end
 458            w_1(i,k,j)  = w_2(i,k,j)
 459            ph_1(i,k,j) = ph_2(i,k,j)
 460          ENDDO
 461          ENDDO
 462          ENDDO
 463
 464        DO j=j_start, j_end
 465          DO i=i_start, i_end
 466            muts(i,j)=mub(i,j)+mu_2(i,j)
 467          ENDDO
 468
 469          DO i=i_start, i_endu
 470            muus(i,j) = muu(i,j)
 471          ENDDO
 472        ENDDO
 473
 474        DO j=j_start, j_endv
 475        DO i=i_start, i_end
 476            muvs(i,j) = muv(i,j)
 477        ENDDO
 478        ENDDO
 479
 480        DO j=j_start, j_end
 481          DO i=i_start, i_end
 482            mu_save(i,j)=mu_2(i,j)
 483            mu_2(i,j)=mu_2(i,j)-mu_2(i,j)
 484          ENDDO
 485        ENDDO
 486        ELSE
 487
 488        DO j=j_start, j_end
 489          DO i=i_start, i_end
 490            muts(i,j)=mub(i,j)+mu_1(i,j)
 491          ENDDO
 492
 493          DO i=i_start, i_endu
 494            muus(i,j)=0.5*(mub(i,j)+mu_1(i,j)+mub(i-1,j)+mu_1(i-1,j))
 495          ENDDO
 496        ENDDO
 497
 498        DO j=j_start, j_endv
 499        DO i=i_start, i_end
 500          muvs(i,j)=0.5*(mub(i,j)+mu_1(i,j)+mub(i,j-1)+mu_1(i,j-1))
 501        ENDDO
 502        ENDDO
 503
 504        DO j=j_start, j_end
 505          DO i=i_start, i_end
 506            mu_save(i,j)=mu_2(i,j)
 507            mu_2(i,j)=mu_1(i,j)-mu_2(i,j)
 508          ENDDO
 509        ENDDO
 510
 511    END IF
 512
 513 !LPB[2]
 514          DO j=j_start, j_end
 515
 516          DO i=i_start, i_end
 517            ww_save(i,kde,j) = 0.
 518            ww_save(i,1,j) = 0.
 519          ENDDO
 520
 521          ENDDO
 522
 523 !LPB[3]
 524        DO j=j_start, j_end
 525
 526        DO k=k_start, k_end
 527        DO i=i_start, i_end
 528          c2a(i,k,j) = cpovcv*(pb(i,k,j)+p(i,k,j))/alt(i,k,j)
 529        ENDDO
 530        ENDDO
 531
 532        ENDDO
 533
 534 !LPB[4]
 535        DO j=j_start, j_end
 536
 537        DO k=k_start, k_end
 538        DO i=i_start, i_endu
 539        Keep_Lpb4_u_2(i,k,j) =u_2(i,k,j)
 540        END DO
 541        END DO
 542
 543        DO k=k_start, k_end
 544        DO i=i_start, i_endu
 545          u_save(i,k,j) = u_2(i,k,j)
 546          u_2(i,k,j) = (muus(i,j)*u_1(i,k,j)-muu(i,j)*u_2(i,k,j))/msfuy(i,j)
 547        ENDDO
 548        ENDDO
 549
 550        ENDDO
 551
 552 !LPB[5]
 553        DO j=j_start, j_endv
 554
 555        DO k=k_start, k_end
 556        DO i=i_start, i_end
 557        Keep_Lpb5_v_2(i,k,j) =v_2(i,k,j)
 558        END DO
 559        END DO
 560
 561        DO k=k_start, k_end
 562        DO i=i_start, i_end
 563          v_save(i,k,j) = v_2(i,k,j)
 564          v_2(i,k,j) = (muvs(i,j)*v_1(i,k,j)-muv(i,j)*v_2(i,k,j))*msfvx_inv(i,j)
 565        ENDDO
 566        ENDDO
 567
 568        ENDDO
 569
 570 !LPB[6]
 571        DO j=j_start, j_end
 572
 573        DO k=k_start, k_end
 574        DO i=i_start, i_end
 575        Keep_Lpb6_t_2(i,k,j) =t_2(i,k,j)
 576        END DO
 577        END DO
 578
 579        DO k=k_start, k_end
 580        DO i=i_start, i_end
 581          t_save(i,k,j) = t_2(i,k,j)
 582          t_2(i,k,j) = muts(i,j)*t_1(i,k,j)-mut(i,j)*t_2(i,k,j)
 583        ENDDO
 584        ENDDO
 585
 586        ENDDO
 587
 588 !!LPB[7]
 589 !       DO j=j_start, j_end
 590
 591 !!      DO k=k_start, kde
 592 !!      DO i=i_start, i_end
 593 !    !  Keep_Lpb7_w_2(i,k,j) =w_2(i,k,j)
 594 !!      END DO
 595 !!      END DO
 596
 597 !       DO k=k_start, kde
 598 !       DO i=i_start, i_end
 599 !         w_save(i,k,j) = w_2(i,k,j)
 600 !         w_2(i,k,j)  = (muts(i,j)* w_1(i,k,j)-mut(i,j)* w_2(i,k,j))/msfty(i,j)
 601 !         ph_save(i,k,j) = ph_2(i,k,j)
 602 !         ph_2(i,k,j) = ph_1(i,k,j)-ph_2(i,k,j)
 603 !       ENDDO
 604 !       ENDDO
 605
 606 !       ENDDO
 607
 608 !!LPB[8]
 609 !         DO j=j_start, j_end
 610
 611 !         DO k=k_start, kde
 612 !         DO i=i_start, i_end
 613 !           ww_save(i,k,j) = ww(i,k,j)
 614 !         ENDDO
 615 !         ENDDO
 616
 617 !         ENDDO
 618
 619 !PART IV: REVERSE/BACKWARD ACCUMULATIONS
 620
 621 !LPB[8]
 622    DO j =j_end, j_start, -1
 623
 624 !  DO k =k_start, kde
 625 !  DO i =i_start, i_end
 626 !  ww_save(i,k,j) =ww(i,k,j)
 627
 628 !  ENDDO
 629 !  ENDDO
 630
 631    DO k =kde, k_start, -1
 632    DO i =i_end, i_start, -1
 633    a_ww(i,k,j) =a_ww(i,k,j) +a_ww_save(i,k,j)
 634    a_ww_save(i,k,j) =0.0
 635    ENDDO
 636    ENDDO
 637
 638    ENDDO
 639
 640 !LPB[7]
 641    DO j =j_end, j_start, -1
 642
 643 !  DO k=k_start, kde
 644 !  DO i=i_start, i_end
 645 !  w_2(i,k,j) =Keep_Lpb7_w_2(i,k,j)
 646 !  END DO
 647 !  END DO
 648
 649    DO k =k_start, kde
 650    DO i =i_start, i_end
 651 !  w_save(i,k,j) =w_2(i,k,j)
 652
 653    Tmpv001 =muts(i,j)*w_1(i,k,j)
 654    Tmpv002 =mut(i,j)*w_2(i,k,j)
 655    Tmpv003 =Tmpv001 -Tmpv002
 656    Tmpv004 =Tmpv003/msfty(i,j)
 657    Tmpv300(i,k) =w_2(i,k,j)
 658    w_2(i,k,j) =Tmpv004
 659
 660 !  ph_save(i,k,j) =ph_2(i,k,j)
 661
 662    Tmpv001 =ph_1(i,k,j) -ph_2(i,k,j)
 663 !  ph_2(i,k,j) =Tmpv001
 664
 665    ENDDO
 666    ENDDO
 667
 668    DO k =kde, k_start, -1
 669    DO i =i_end, i_start, -1
 670    a_Tmpv1 =a_ph_2(i,k,j)
 671    a_ph_2(i,k,j) =0.0
 672    a_ph_1(i,k,j) =a_ph_1(i,k,j) +a_Tmpv1
 673    a_ph_2(i,k,j) =a_ph_2(i,k,j) -a_Tmpv1
 674    a_ph_2(i,k,j) =a_ph_2(i,k,j) +a_ph_save(i,k,j)
 675    a_ph_save(i,k,j) =0.0
 676
 677    w_2(i,k,j) =Tmpv300(i,k)
 678
 679    a_Tmpv4 =a_w_2(i,k,j)
 680    a_w_2(i,k,j) =0.0
 681    a_Tmpv3 =a_Tmpv4/msfty(i,j)
 682    a_Tmpv1 =a_Tmpv3
 683    a_Tmpv2 =-a_Tmpv3
 684    a_mut(i,j) =a_mut(i,j) +w_2(i,k,j)*a_Tmpv2
 685    a_w_2(i,k,j) =a_w_2(i,k,j) +mut(i,j)*a_Tmpv2
 686    a_muts(i,j) =a_muts(i,j) +w_1(i,k,j)*a_Tmpv1
 687    a_w_1(i,k,j) =a_w_1(i,k,j) +muts(i,j)*a_Tmpv1
 688    a_w_2(i,k,j) =a_w_2(i,k,j) +a_w_save(i,k,j)
 689    a_w_save(i,k,j) =0.0
 690    ENDDO
 691    ENDDO
 692
 693    ENDDO
 694
 695 !LPB[6]
 696    DO j =j_end, j_start, -1
 697
 698    DO k=k_start, k_end
 699    DO i=i_start, i_end
 700    t_2(i,k,j) =Keep_Lpb6_t_2(i,k,j)
 701    END DO
 702    END DO
 703
 704    DO k =k_start, k_end
 705    DO i =i_start, i_end
 706 !  t_save(i,k,j) =t_2(i,k,j)
 707
 708    Tmpv001 =muts(i,j)*t_1(i,k,j)
 709    Tmpv002 =mut(i,j)*t_2(i,k,j)
 710    Tmpv003 =Tmpv001 -Tmpv002
 711    Tmpv300(i,k) =t_2(i,k,j)
 712    t_2(i,k,j) =Tmpv003
 713
 714    ENDDO
 715    ENDDO
 716
 717    DO k =k_end, k_start, -1
 718    DO i =i_end, i_start, -1
 719
 720    t_2(i,k,j) =Tmpv300(i,k)
 721
 722    a_Tmpv3 =a_t_2(i,k,j)
 723    a_t_2(i,k,j) =0.0
 724    a_Tmpv1 =a_Tmpv3
 725    a_Tmpv2 =-a_Tmpv3
 726    a_mut(i,j) =a_mut(i,j) +t_2(i,k,j)*a_Tmpv2
 727    a_t_2(i,k,j) =a_t_2(i,k,j) +mut(i,j)*a_Tmpv2
 728    a_muts(i,j) =a_muts(i,j) +t_1(i,k,j)*a_Tmpv1
 729    a_t_1(i,k,j) =a_t_1(i,k,j) +muts(i,j)*a_Tmpv1
 730    a_t_2(i,k,j) =a_t_2(i,k,j) +a_t_save(i,k,j)
 731    a_t_save(i,k,j) =0.0
 732    ENDDO
 733    ENDDO
 734
 735    ENDDO
 736
 737 !LPB[5]
 738    DO j =j_endv, j_start, -1
 739
 740    DO k=k_start, k_end
 741    DO i=i_start, i_end
 742    v_2(i,k,j) =Keep_Lpb5_v_2(i,k,j)
 743 !  IF(a_v_2(i,k,j).NE.0.0) PRINT*, 'a_v_2(i,k,j)=', a_v_2(i,k,j)
 744    END DO
 745    END DO
 746
 747    DO k =k_start, k_end
 748    DO i =i_start, i_end
 749 !  v_save(i,k,j) =v_2(i,k,j)
 750
 751    Tmpv001 =muvs(i,j)*v_1(i,k,j)
 752    Tmpv002 =muv(i,j)*v_2(i,k,j)
 753    Tmpv003 =Tmpv001 -Tmpv002
 754    Tmpv004 =Tmpv003*msfvx_inv(i,j)
 755    Tmpv300(i,k) =v_2(i,k,j)
 756    v_2(i,k,j) =Tmpv004
 757
 758    ENDDO
 759    ENDDO
 760
 761    DO k =k_end, k_start, -1
 762    DO i =i_end, i_start, -1
 763
 764    v_2(i,k,j) =Tmpv300(i,k)
 765
 766    a_Tmpv4 =a_v_2(i,k,j)
 767    a_v_2(i,k,j) =0.0
 768    a_Tmpv3 =msfvx_inv(i,j)*a_Tmpv4
 769    a_Tmpv1 =a_Tmpv3
 770    a_Tmpv2 =-a_Tmpv3
 771    a_muv(i,j) =a_muv(i,j) +v_2(i,k,j)*a_Tmpv2
 772    a_v_2(i,k,j) =a_v_2(i,k,j) +muv(i,j)*a_Tmpv2
 773    a_muvs(i,j) =a_muvs(i,j) +v_1(i,k,j)*a_Tmpv1
 774    a_v_1(i,k,j) =a_v_1(i,k,j) +muvs(i,j)*a_Tmpv1
 775    a_v_2(i,k,j) =a_v_2(i,k,j) +a_v_save(i,k,j)
 776    a_v_save(i,k,j) =0.0
 777    ENDDO
 778    ENDDO
 779
 780    ENDDO
 781
 782 !LPB[4]
 783    DO j =j_end, j_start, -1
 784
 785    DO k=k_start, k_end
 786    DO i=i_start, i_endu
 787    u_2(i,k,j) =Keep_Lpb4_u_2(i,k,j)
 788    END DO
 789    END DO
 790
 791    DO k =k_start, k_end
 792    DO i =i_start, i_endu
 793 !  u_save(i,k,j) =u_2(i,k,j)
 794
 795    Tmpv001 =muus(i,j)*u_1(i,k,j)
 796    Tmpv002 =muu(i,j)*u_2(i,k,j)
 797    Tmpv003 =Tmpv001 -Tmpv002
 798    Tmpv004 =Tmpv003/msfuy(i,j)
 799    Tmpv300(i,k) =u_2(i,k,j)
 800    u_2(i,k,j) =Tmpv004
 801
 802    ENDDO
 803    ENDDO
 804
 805    DO k =k_end, k_start, -1
 806    DO i =i_endu, i_start, -1
 807
 808    u_2(i,k,j) =Tmpv300(i,k)
 809
 810    a_Tmpv4 =a_u_2(i,k,j)
 811    a_u_2(i,k,j) =0.0
 812    a_Tmpv3 =a_Tmpv4/msfuy(i,j)
 813    a_Tmpv1 =a_Tmpv3
 814    a_Tmpv2 =-a_Tmpv3
 815    a_muu(i,j) =a_muu(i,j) +u_2(i,k,j)*a_Tmpv2
 816    a_u_2(i,k,j) =a_u_2(i,k,j) +muu(i,j)*a_Tmpv2
 817    a_muus(i,j) =a_muus(i,j) +u_1(i,k,j)*a_Tmpv1
 818    a_u_1(i,k,j) =a_u_1(i,k,j) +muus(i,j)*a_Tmpv1
 819    a_u_2(i,k,j) =a_u_2(i,k,j) +a_u_save(i,k,j)
 820    a_u_save(i,k,j) =0.0
 821    ENDDO
 822    ENDDO
 823
 824    ENDDO
 825
 826 !LPB[3]
 827    DO j =j_end, j_start, -1
 828
 829 !  DO k =k_start, k_end
 830 !  DO i =i_start, i_end
 831 !  Tmpv001 =cpovcv*(pb(i,k,j) +p(i,k,j))/alt(i,k,j)
 832 !  c2a(i,k,j) =Tmpv001
 833
 834 !  ENDDO
 835 !  ENDDO
 836
 837    DO k =k_end, k_start, -1
 838    DO i =i_end, i_start, -1
 839    a_Tmpv1 =a_c2a(i,k,j)
 840    a_c2a(i,k,j) =0.0
 841    a_p(i,k,j) =a_p(i,k,j) +cpovcv/alt(i,k,j)*a_Tmpv1
 842    a_alt(i,k,j) =a_alt(i,k,j) -cpovcv*(pb(i,k,j) +p(i,k,j))/(alt(i,k,j)  &
 843    *alt(i,k,j))*a_Tmpv1
 844    ENDDO
 845    ENDDO
 846
 847    ENDDO
 848
 849 !LPB[2]
 850    DO j =j_end, j_start, -1
 851
 852 !  DO i =i_start, i_end
 853 !  ww_save(i,kde,j) =0.
 854
 855 !  ww_save(i,1,j) =0.
 856
 857 !  ENDDO
 858
 859    DO i =i_end, i_start, -1
 860    a_ww_save(i,1,j) =0.0
 861    a_ww_save(i,kde,j) =0.0
 862    ENDDO
 863
 864    ENDDO
 865
 866 !LPB[1]
 867
 868 !  IF((rk_step == 1) ) THEN
 869 !  DO j =j_start, j_end
 870 !  DO i =i_start, i_end
 871 !  mu_1(i,j) =mu_2(i,j)
 872
 873 !  ww_save(i,kde,j) =0.
 874
 875 !  ww_save(i,1,j) =0.
 876
 877 !  mudf(i,j) =0.
 878
 879 !  ENDDO
 880 !  ENDDO
 881 !  DO j =j_start, j_end
 882 !  DO k =k_start, k_end
 883 !  DO i =i_start, i_endu
 884 !  u_1(i,k,j) =u_2(i,k,j)
 885
 886 !  ENDDO
 887 !  ENDDO
 888 !  ENDDO
 889 !  DO j =j_start, j_endv
 890 !  DO k =k_start, k_end
 891 !  DO i =i_start, i_end
 892 !  v_1(i,k,j) =v_2(i,k,j)
 893
 894 !  ENDDO
 895 !  ENDDO
 896 !  ENDDO
 897 !  DO j =j_start, j_end
 898 !  DO k =k_start, k_end
 899 !  DO i =i_start, i_end
 900 !  t_1(i,k,j) =t_2(i,k,j)
 901
 902 !  ENDDO
 903 !  ENDDO
 904 !  ENDDO
 905 !  DO j =j_start, j_end
 906 !  DO k =k_start, min(kde, kte)
 907 !  DO i =i_start, i_end
 908 !  w_1(i,k,j) =w_2(i,k,j)
 909
 910 !  ph_1(i,k,j) =ph_2(i,k,j)
 911
 912 !  ENDDO
 913 !  ENDDO
 914 !  ENDDO
 915 !  DO j =j_start, j_end
 916 !  DO i =i_start, i_end
 917 !  muts(i,j) =mub(i,j) +mu_2(i,j)
 918
 919 !  ENDDO
 920
 921 !  DO i =i_start, i_endu
 922 !  muus(i,j) =muu(i,j)
 923
 924 !  ENDDO
 925 !  ENDDO
 926 !  DO j =j_start, j_endv
 927 !  DO i =i_start, i_end
 928 !  muvs(i,j) =muv(i,j)
 929
 930 !  ENDDO
 931 !  ENDDO
 932 !  DO j =j_start, j_end
 933 !  DO i =i_start, i_end
 934 !  mu_save(i,j) =mu_2(i,j)
 935
 936 !  mu_2(i,j) =mu_2(i,j) -mu_2(i,j)
 937
 938 !  ENDDO
 939 !  ENDDO
 940 !  ELSE
 941 !  DO j =j_start, j_end
 942 !  DO i =i_start, i_end
 943 !  muts(i,j) =mub(i,j) +mu_1(i,j)
 944
 945 !  ENDDO
 946
 947 !  DO i =i_start, i_endu
 948 !  Tmpv001 =mub(i,j) +mu_1(i,j) +mub(i-1,j) +mu_1(i-1,j)
 949 !  Tmpv002 =0.5*Tmpv001
 950 !  muus(i,j) =Tmpv002
 951
 952 !  ENDDO
 953 !  ENDDO
 954 !  DO j =j_start, j_endv
 955 !  DO i =i_start, i_end
 956 !  Tmpv001 =mub(i,j) +mu_1(i,j) +mub(i,j-1) +mu_1(i,j-1)
 957 !  Tmpv002 =0.5*Tmpv001
 958 !  muvs(i,j) =Tmpv002
 959
 960 !  ENDDO
 961 !  ENDDO
 962 !  DO j =j_start, j_end
 963 !  DO i =i_start, i_end
 964 !  mu_save(i,j) =mu_2(i,j)
 965
 966 !  Tmpv001 =mu_1(i,j) -mu_2(i,j)
 967 !  mu_2(i,j) =Tmpv001
 968
 969 !  ENDDO
 970 !  ENDDO
 971 !  END IF
 972
 973    IF((rk_step == 1) ) THEN
 974
 975    DO j =j_end, j_start, -1
 976    DO i =i_end, i_start, -1
 977 !BIG ERROR HERE, REVISED BY WALLS
 978 !  a_mu_2(i,j) =1.0 -1.0*a_mu_2(i,j)
 979    a_mu_2(i,j) =0.0  !REVISED BY WALLS
 980
 981    a_mu_2(i,j) =a_mu_2(i,j) +a_mu_save(i,j)
 982    a_mu_save(i,j) =0.0
 983    ENDDO
 984    ENDDO
 985    DO j =j_endv, j_start, -1
 986    DO i =i_end, i_start, -1
 987    a_muv(i,j) =a_muv(i,j) +a_muvs(i,j)
 988    a_muvs(i,j) =0.0
 989    ENDDO
 990    ENDDO
 991    DO j =j_end, j_start, -1
 992    DO i =i_endu, i_start, -1
 993    a_muu(i,j) =a_muu(i,j) +a_muus(i,j)
 994    a_muus(i,j) =0.0
 995    ENDDO
 996    DO i =i_end, i_start, -1
 997    a_mu_2(i,j) =a_mu_2(i,j) +a_muts(i,j)
 998    a_muts(i,j) =0.0
 999    ENDDO
1000    ENDDO
1001    DO j =j_end, j_start, -1
1002    DO k =min(kde, kte), k_start, -1
1003    DO i =i_end, i_start, -1
1004    a_ph_2(i,k,j) =a_ph_2(i,k,j) +a_ph_1(i,k,j)
1005    a_ph_1(i,k,j) =0.0
1006    a_w_2(i,k,j) =a_w_2(i,k,j) +a_w_1(i,k,j)
1007    a_w_1(i,k,j) =0.0
1008    ENDDO
1009    ENDDO
1010    ENDDO
1011    DO j =j_end, j_start, -1
1012    DO k =k_end, k_start, -1
1013    DO i =i_end, i_start, -1
1014    a_t_2(i,k,j) =a_t_2(i,k,j) +a_t_1(i,k,j)
1015    a_t_1(i,k,j) =0.0
1016    ENDDO
1017    ENDDO
1018    ENDDO
1019    DO j =j_endv, j_start, -1
1020    DO k =k_end, k_start, -1
1021    DO i =i_end, i_start, -1
1022    a_v_2(i,k,j) =a_v_2(i,k,j) +a_v_1(i,k,j)
1023    a_v_1(i,k,j) =0.0
1024    ENDDO
1025    ENDDO
1026    ENDDO
1027    DO j =j_end, j_start, -1
1028    DO k =k_end, k_start, -1
1029    DO i =i_endu, i_start, -1
1030    a_u_2(i,k,j) =a_u_2(i,k,j) +a_u_1(i,k,j)
1031    a_u_1(i,k,j) =0.0
1032    ENDDO
1033    ENDDO
1034    ENDDO
1035    DO j =j_end, j_start, -1
1036    DO i =i_end, i_start, -1
1037    a_mudf(i,j) =0.0
1038    a_ww_save(i,1,j) =0.0
1039    a_ww_save(i,kde,j) =0.0
1040    a_mu_2(i,j) =a_mu_2(i,j) +a_mu_1(i,j)
1041    a_mu_1(i,j) =0.0
1042    ENDDO
1043    ENDDO
1044
1045    ELSE
1046
1047    DO j =j_end, j_start, -1
1048    DO i =i_end, i_start, -1
1049    a_Tmpv1 =a_mu_2(i,j)
1050    a_mu_2(i,j) =0.0
1051    a_mu_1(i,j) =a_mu_1(i,j) +a_Tmpv1
1052    a_mu_2(i,j) =a_mu_2(i,j) -a_Tmpv1
1053    a_mu_2(i,j) =a_mu_2(i,j) +a_mu_save(i,j)
1054    a_mu_save(i,j) =0.0
1055    ENDDO
1056    ENDDO
1057    DO j =j_endv, j_start, -1
1058    DO i =i_end, i_start, -1
1059    a_Tmpv2 =a_muvs(i,j)
1060    a_muvs(i,j) =0.0
1061    a_Tmpv1 =0.5*a_Tmpv2
1062    a_mu_1(i,j) =a_mu_1(i,j) +a_Tmpv1
1063    a_mu_1(i,j-1) =a_mu_1(i,j-1) +a_Tmpv1
1064    ENDDO
1065    ENDDO
1066    DO j =j_end, j_start, -1
1067    DO i =i_endu, i_start, -1
1068    a_Tmpv2 =a_muus(i,j)
1069    a_muus(i,j) =0.0
1070    a_Tmpv1 =0.5*a_Tmpv2
1071    a_mu_1(i,j) =a_mu_1(i,j) +a_Tmpv1
1072    a_mu_1(i-1,j) =a_mu_1(i-1,j) +a_Tmpv1
1073    ENDDO
1074    DO i =i_end, i_start, -1
1075    a_mu_1(i,j) =a_mu_1(i,j) +a_muts(i,j)
1076    a_muts(i,j) =0.0
1077    ENDDO
1078    ENDDO
1079
1080    END IF
1081
1082 !LPB[0]
1083 !  i_start =its
1084 !  i_end =min(ite, ide-1)
1085 !  j_start =jts
1086 !  j_end =min(jte, jde-1)
1087 !  k_start =kts
1088 !  k_end =min(kte, kde-1)
1089 !  i_endu =ite
1090 !  j_endv =jte
1091
1092    END SUBROUTINE a_small_step_prep
1093
1094 !        Generated by TAPENADE     (INRIA, Tropics team)
1095 !  Tapenade 3.5 (r3785) - 22 Mar 2011 18:35
1096 !
1097 !  Differentiation of small_step_finish in reverse (adjoint) mode:
1098 !   gradient     of useful results: u_save ph_save ww w_save mu_save
1099 !                muvs ph_2 u_2 h_diabatic mu_2 w_2 t_save v_save
1100 !                muts t_2 mut muu v_2 muv ww1 muus
1101 !   with respect to varying inputs: u_save ph_save ww w_save mu_save
1102 !                muvs ph_2 u_2 h_diabatic mu_2 w_2 t_save v_save
1103 !                muts t_2 mut muu v_2 muv ww1 muus
1104 !   RW status of diff variables: u_save:incr ph_save:incr ww:in-out
1105 !                w_save:incr mu_save:incr muvs:incr ph_2:in-out
1106 !                u_2:in-out h_diabatic:incr mu_2:in-out w_2:in-out
1107 !                t_save:incr v_save:incr muts:incr t_2:in-out mut:incr
1108 !                muu:incr v_2:in-out muv:incr ww1:incr muus:incr
1109 SUBROUTINE A_SMALL_STEP_FINISH(u_2, u_2b, u_1, v_2, v_2b, v_1, w_2, w_2b&
1110 &  , w_1, t_2, t_2b, t_1, ph_2, ph_2b, ph_1, ww, wwb, ww1, ww1b, mu_2, &
1111 &  mu_2b, mu_1, mut, mutb, muts, mutsb, muu, muub, muus, muusb, muv, muvb&
1112 &  , muvs, muvsb, u_save, u_saveb, v_save, v_saveb, w_save, w_saveb, &
1113 &  t_save, t_saveb, ph_save, ph_saveb, mu_save, mu_saveb, msfux, msfuy, &
1114 &  msfvx, msfvy, msftx, msfty, h_diabatic, h_diabaticb, &
1115 &  number_of_small_timesteps, dts, rk_step, rk_order, ids, ide, jds, jde&
1116 &  , kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte&
1117 &)
1118   IMPLICIT NONE
1119 ! religion first
1120 !  stuff passed in
1121   INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde
1122   INTEGER, INTENT(IN) :: ims, ime, jms, jme, kms, kme
1123   INTEGER, INTENT(IN) :: its, ite, jts, jte, kts, kte
1124   INTEGER, INTENT(IN) :: number_of_small_timesteps
1125   INTEGER, INTENT(IN) :: rk_step, rk_order
1126   REAL, INTENT(IN) :: dts
1127   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u_1, v_1, &
1128 &  w_1, t_1, ww1, ph_1
1129   REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: ww1b
1130   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: u_2, v_2&
1131 &  , w_2, t_2, ww, ph_2
1132   REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: u_2b, v_2b, w_2b, t_2b, &
1133 &  wwb, ph_2b
1134   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u_save, &
1135 &  v_save, w_save, t_save, ph_save, h_diabatic
1136   REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: u_saveb, v_saveb, &
1137 &  w_saveb, t_saveb, ph_saveb, h_diabaticb
1138   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: muus, muvs
1139   REAL, DIMENSION(ims:ime, jms:jme) :: muusb, muvsb
1140   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mu_2, mu_1
1141   REAL, DIMENSION(ims:ime, jms:jme) :: mu_2b
1142   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mut, muts, muu, &
1143 &  muv, mu_save
1144   REAL, DIMENSION(ims:ime, jms:jme) :: mutb, mutsb, muub, muvb, mu_saveb
1145   REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: msfux, msfuy, msfvx, &
1146 &  msfvy, msftx, msfty
1147 ! local stuff
1148   INTEGER :: i, j, k
1149   INTEGER :: i_start, i_end, j_start, j_end, i_endu, j_endv
1150   INTEGER :: branch
1151   REAL :: tempb4
1152   REAL :: tempb3
1153   REAL :: tempb2
1154   REAL :: tempb1
1155   REAL :: tempb0
1156   REAL :: tempb
1157   INTRINSIC MIN
1158 !<DESCRIPTION>
1159 !
1160 !  small_step_finish reconstructs the full uncoupled prognostic variables
1161 !  from the coupled perturbation variables used in the small timesteps.
1162 !
1163 !</DESCRIPTION>
1164   i_start = its
1165   IF (ite .GT. ide - 1) THEN
1166     i_end = ide - 1
1167   ELSE
1168     i_end = ite
1169   END IF
1170   j_start = jts
1171   IF (jte .GT. jde - 1) THEN
1172     j_end = jde - 1
1173   ELSE
1174     j_end = jte
1175   END IF
1176   i_endu = ite
1177   j_endv = jte
1178   IF (rk_step .LT. rk_order) THEN
1179     CALL PUSHCONTROL1B(1)
1180   ELSE
1181     CALL PUSHCONTROL1B(0)
1182   END IF
1183   DO j=j_end,j_start,-1
1184     DO i=i_end,i_start,-1
1185       mu_saveb(i, j) = mu_saveb(i, j) + mu_2b(i, j)
1186     END DO
1187   END DO
1188   CALL POPCONTROL1B(branch)
1189   IF (branch .EQ. 0) THEN
1190     DO j=j_end,j_start,-1
1191       DO k=kde-1,kds,-1
1192         DO i=i_end,i_start,-1
1193           tempb3 = t_2b(i, k, j)/muts(i, j)
1194           tempb4 = -(dts*number_of_small_timesteps*tempb3)
1195           mutb(i, j) = mutb(i, j) + t_save(i, k, j)*tempb3 + h_diabatic(&
1196 &            i, k, j)*tempb4
1197           h_diabaticb(i, k, j) = h_diabaticb(i, k, j) + mut(i, j)*tempb4
1198           t_saveb(i, k, j) = t_saveb(i, k, j) + mut(i, j)*tempb3
1199           mutsb(i, j) = mutsb(i, j) - (t_2(i, k, j)-dts*&
1200 &            number_of_small_timesteps*(mut(i, j)*h_diabatic(i, k, j))+&
1201 &            t_save(i, k, j)*mut(i, j))*tempb3/muts(i, j)
1202           t_2b(i, k, j) = tempb3
1203         END DO
1204       END DO
1205     END DO
1206   ELSE
1207     DO j=j_end,j_start,-1
1208       DO k=kde-1,kds,-1
1209         DO i=i_end,i_start,-1
1210           tempb2 = t_2b(i, k, j)/muts(i, j)
1211           t_saveb(i, k, j) = t_saveb(i, k, j) + mut(i, j)*tempb2
1212           mutb(i, j) = mutb(i, j) + t_save(i, k, j)*tempb2
1213           mutsb(i, j) = mutsb(i, j) - (t_2(i, k, j)+t_save(i, k, j)*mut(&
1214 &            i, j))*tempb2/muts(i, j)
1215           t_2b(i, k, j) = tempb2
1216         END DO
1217       END DO
1218     END DO
1219   END IF
1220   DO j=j_end,j_start,-1
1221     DO k=kde,kds,-1
1222       DO i=i_end,i_start,-1
1223         ww1b(i, k, j) = ww1b(i, k, j) + wwb(i, k, j)
1224         ph_saveb(i, k, j) = ph_saveb(i, k, j) + ph_2b(i, k, j)
1225         tempb1 = w_2b(i, k, j)/muts(i, j)
1226         w_saveb(i, k, j) = w_saveb(i, k, j) + mut(i, j)*tempb1
1227         mutb(i, j) = mutb(i, j) + w_save(i, k, j)*tempb1
1228         mutsb(i, j) = mutsb(i, j) - (msfty(i, j)*w_2(i, k, j)+w_save(i, &
1229 &          k, j)*mut(i, j))*tempb1/muts(i, j)
1230         w_2b(i, k, j) = msfty(i, j)*tempb1
1231       END DO
1232     END DO
1233   END DO
1234   DO j=j_end,j_start,-1
1235     DO k=kde-1,kds,-1
1236       DO i=i_endu,i_start,-1
1237         tempb0 = u_2b(i, k, j)/muus(i, j)
1238         u_saveb(i, k, j) = u_saveb(i, k, j) + muu(i, j)*tempb0
1239         muub(i, j) = muub(i, j) + u_save(i, k, j)*tempb0
1240         muusb(i, j) = muusb(i, j) - (msfuy(i, j)*u_2(i, k, j)+u_save(i, &
1241 &          k, j)*muu(i, j))*tempb0/muus(i, j)
1242         u_2b(i, k, j) = msfuy(i, j)*tempb0
1243       END DO
1244     END DO
1245   END DO
1246   DO j=j_endv,j_start,-1
1247     DO k=kde-1,kds,-1
1248       DO i=i_end,i_start,-1
1249         tempb = v_2b(i, k, j)/muvs(i, j)
1250         v_saveb(i, k, j) = v_saveb(i, k, j) + muv(i, j)*tempb
1251         muvb(i, j) = muvb(i, j) + v_save(i, k, j)*tempb
1252         muvsb(i, j) = muvsb(i, j) - (msfvx(i, j)*v_2(i, k, j)+v_save(i, &
1253 &          k, j)*muv(i, j))*tempb/muvs(i, j)
1254         v_2b(i, k, j) = msfvx(i, j)*tempb
1255       END DO
1256     END DO
1257   END DO
1258 END SUBROUTINE A_SMALL_STEP_FINISH
1259
1260    SUBROUTINE a_calc_p_rho(al,a_al,p,a_p,ph,a_ph,alt,a_alt,t_2,a_t_2, &
1261    t_1,a_t_1,c2a,a_c2a,pm1,a_pm1,mu,a_mu,muts,a_muts,znu,t0,rdnw,dnw, &
1262    smdiv,non_hydrostatic,step,ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite, &
1263    jts,jte,kts,kte)
1264
1265 !PART I: DECLARATION OF VARIABLES
1266
1267    IMPLICIT NONE
1268
1269    INTEGER :: K0_ADJ,K1_ADJ,K2_ADJ,K3_ADJ
1270    INTEGER :: ids,ide,jds,jde,kds,kde
1271    INTEGER :: ims,ime,jms,jme,kms,kme
1272    INTEGER :: its,ite,jts,jte,kts,kte
1273    INTEGER :: step
1274    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_al,al,a_p,p
1275    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_alt,alt,a_t_2,t_2,a_t_1,t_1,a_c2a,c2a
1276    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_ph,ph,a_pm1,pm1
1277    REAL,DIMENSION(ims:ime,jms:jme) :: a_mu,mu,a_muts,muts
1278    REAL,DIMENSION(kms:kme) :: dnw,rdnw,znu
1279    REAL :: t0,smdiv
1280    LOGICAL :: non_hydrostatic
1281    INTEGER :: i,j,k
1282    INTEGER :: i_start,i_end,j_start,j_end,k_start,k_end
1283    REAL :: a_ptmp,ptmp
1284
1285 !  REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Keep_Lpb1_ph
1286    INTEGER :: IX1,IX2,IX3
1287
1288    REAL :: a_Tmpv1,Tmpv001,a_Tmpv2,Tmpv002,a_Tmpv3,Tmpv003,a_Tmpv4,Tmpv004, &
1289    a_Tmpv5,Tmpv005,a_Tmpv6,Tmpv006,a_Tmpv7,Tmpv007
1290    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Tmpv400
1291    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Tmpv401
1292    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Tmpv402
1293    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Tmpv403
1294    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Tmpv404
1295    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Tmpv405
1296    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Tmpv406
1297    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Tmpv407
1298    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Tmpv408
1299    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Tmpv409
1300    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Tmpv4010
1301    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Tmpv4011
1302
1303 !PART II: CALCULATIONS OF B. S. TRAJECTORY
1304
1305 !LPB[0]
1306       i_start = its
1307       i_end   = min(ite,ide-1)
1308       j_start = jts
1309       j_end   = min(jte,jde-1)
1310       k_start = kts
1311       k_end = min(kte,kde-1)
1312
1313 !!LPB[1]
1314 !!  DO IX3=jms,jme
1315 !!  DO IX2=kms,kme
1316 !!  DO IX1=ims,ime
1317 !    !  Keep_Lpb1_ph(IX1,IX2,IX3) =ph(IX1,IX2,IX3)
1318 !!  END DO
1319 !!  END DO
1320 !!  END DO
1321
1322 !   IF (non_hydrostatic) THEN
1323
1324 !        DO j=j_start, j_end
1325 !        DO k=k_start, k_end
1326 !        DO i=i_start, i_end
1327 !         al(i,k,j)=-1./muts(i,j)*(alt(i,k,j)*mu(i,j)    &
1328 !                +rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
1329 !         p(i,k,j)=c2a(i,k,j)*(alt(i,k,j)*(t_2(i,k,j)-mu(i,j)*t_1(i,k,j))    &
1330 !                          /(muts(i,j)*(t0+t_1(i,k,j)))-al (i,k,j))
1331 !        ENDDO
1332 !        ENDDO
1333 !        ENDDO
1334 !      ELSE
1335
1336 !          DO j=j_start, j_end
1337 !          DO k=k_start, k_end
1338 !          DO i=i_start, i_end
1339 !            p(i,k,j)=mu(i,j)*znu(k)
1340 !            al(i,k,j)=alt(i,k,j)*(t_2(i,k,j)-mu(i,j)*t_1(i,k,j))              &
1341 !                         /(muts(i,j)*(t0+t_1(i,k,j)))-p(i,k,j)/c2a(i,k,j)
1342 !            ph(i,k+1,j)=ph(i,k,j)-dnw(k)*(muts(i,j)*al (i,k,j)                &
1343 !                             +mu(i,j)*alt(i,k,j))
1344 !          ENDDO
1345 !          ENDDO
1346 !          ENDDO
1347
1348 !   END IF
1349
1350 !!LPB[2]
1351
1352 !!LPB[3]
1353 !
1354 !     IF (step == 0) then
1355
1356 !          DO j=j_start, j_end
1357 !          DO k=k_start, k_end
1358 !          DO i=i_start, i_end
1359 !            pm1(i,k,j)=p(i,k,j)
1360 !          ENDDO
1361 !          ENDDO
1362 !          ENDDO
1363 !        ELSE
1364
1365 !          DO j=j_start, j_end
1366 !          DO k=k_start, k_end
1367 !          DO i=i_start, i_end
1368 !            ptmp = p(i,k,j)
1369 !            p(i,k,j) = p(i,k,j) + smdiv*(p(i,k,j)-pm1(i,k,j))
1370 !            pm1(i,k,j) = ptmp
1371 !          ENDDO
1372 !          ENDDO
1373 !          ENDDO
1374
1375 !   END IF
1376
1377 !PART III: INITIALIZATION OF LOCAL ADJOINT PERTURBATIONS
1378
1379    a_ptmp =0.0
1380
1381 !PART IV: REVERSE/BACKWARD ACCUMULATIONS
1382
1383 !LPB[3]
1384
1385 !  IF(step == 0) THEN
1386 !  DO j =j_start, j_end
1387 !  DO k =k_start, k_end
1388 !  DO i =i_start, i_end
1389 !  pm1(i,k,j) =p(i,k,j)
1390
1391 !  ENDDO
1392 !  ENDDO
1393 !  ENDDO
1394 !  ELSE
1395 !  DO j =j_start, j_end
1396 !  DO k =k_start, k_end
1397 !  DO i =i_start, i_end
1398 !  ptmp =p(i,k,j)
1399
1400 !  Tmpv001 =p(i,k,j) -pm1(i,k,j)
1401 !  Tmpv002 =smdiv*Tmpv001
1402 !  Tmpv003 =p(i,k,j) +Tmpv002
1403 !  p(i,k,j) =Tmpv003
1404
1405 !  pm1(i,k,j) =ptmp
1406
1407 !  ENDDO
1408 !  ENDDO
1409 !  ENDDO
1410 !  END IF
1411
1412    IF(step == 0) THEN
1413
1414    DO j =j_end, j_start, -1
1415    DO k =k_end, k_start, -1
1416    DO i =i_end, i_start, -1
1417    a_p(i,k,j) =a_p(i,k,j) +a_pm1(i,k,j)
1418    a_pm1(i,k,j) =0.0
1419    ENDDO
1420    ENDDO
1421    ENDDO
1422
1423    ELSE
1424
1425    DO j =j_end, j_start, -1
1426    DO k =k_end, k_start, -1
1427    DO i =i_end, i_start, -1
1428    a_ptmp =a_ptmp +a_pm1(i,k,j)
1429    a_pm1(i,k,j) =0.0
1430    a_Tmpv3 =a_p(i,k,j)
1431    a_p(i,k,j) =0.0
1432    a_p(i,k,j) =a_p(i,k,j) +a_Tmpv3
1433    a_Tmpv2 =a_Tmpv3
1434    a_Tmpv1 =smdiv*a_Tmpv2
1435    a_p(i,k,j) =a_p(i,k,j) +a_Tmpv1
1436    a_pm1(i,k,j) =a_pm1(i,k,j) -a_Tmpv1
1437    a_p(i,k,j) =a_p(i,k,j) +a_ptmp
1438    a_ptmp =0.0
1439    ENDDO
1440    ENDDO
1441    ENDDO
1442
1443    END IF
1444
1445 !LPB[2]
1446
1447 !LPB[1]
1448 !  DO IX3=jms,jme
1449 !  DO IX2=kms,kme
1450 !  DO IX1=ims,ime
1451 !  ph(IX1,IX2,IX3) =Keep_Lpb1_ph(IX1,IX2,IX3)
1452 !  END DO
1453 !  END DO
1454 !  END DO
1455
1456    IF(non_hydrostatic) THEN
1457    DO j =j_start, j_end
1458    DO k =k_start, k_end
1459    DO i =i_start, i_end
1460    Tmpv001 =alt(i,k,j)*mu(i,j)
1461    Tmpv002 =ph(i,k+1,j) -ph(i,k,j)
1462    Tmpv003 =rdnw(k)*Tmpv002
1463    Tmpv004 =Tmpv001 +Tmpv003
1464    Tmpv400(i,k,j) =Tmpv004
1465    Tmpv005 =-1./muts(i,j)*Tmpv400(i,k,j)
1466    Tmpv401(i,k,j) =al(i,k,j)
1467    al(i,k,j) =Tmpv005
1468
1469    Tmpv001 =mu(i,j)*t_1(i,k,j)
1470    Tmpv002 =t_2(i,k,j) -Tmpv001
1471    Tmpv402(i,k,j) =Tmpv002
1472    Tmpv003 =alt(i,k,j)*Tmpv402(i,k,j)
1473    Tmpv004 =muts(i,j)*(t0 +t_1(i,k,j))
1474    Tmpv403(i,k,j) =Tmpv003
1475    Tmpv404(i,k,j) =Tmpv004
1476    Tmpv005 =Tmpv403(i,k,j)/Tmpv404(i,k,j)
1477    Tmpv006 =Tmpv005 -al(i,k,j)
1478    Tmpv405(i,k,j) =Tmpv006
1479    Tmpv007 =c2a(i,k,j)*Tmpv405(i,k,j)
1480    Tmpv406(i,k,j) =p(i,k,j)
1481    p(i,k,j) =Tmpv007
1482
1483    ENDDO
1484    ENDDO
1485    ENDDO
1486    ELSE
1487    DO j =j_start, j_end
1488    DO k =k_start, k_end
1489    DO i =i_start, i_end
1490    Tmpv407(i,k,j) =p(i,k,j)
1491    p(i,k,j) =mu(i,j)*znu(k)
1492
1493    Tmpv001 =mu(i,j)*t_1(i,k,j)
1494    Tmpv002 =t_2(i,k,j) -Tmpv001
1495    Tmpv408(i,k,j) =Tmpv002
1496    Tmpv003 =alt(i,k,j)*Tmpv408(i,k,j)
1497    Tmpv004 =muts(i,j)*(t0 +t_1(i,k,j))
1498    Tmpv409(i,k,j) =Tmpv003
1499    Tmpv4010(i,k,j) =Tmpv004
1500    Tmpv005 =Tmpv409(i,k,j)/Tmpv4010(i,k,j)
1501    Tmpv006 =p(i,k,j)/c2a(i,k,j)
1502    Tmpv007 =Tmpv005 -Tmpv006
1503    Tmpv4011(i,k,j) =al(i,k,j)
1504    al(i,k,j) =Tmpv007
1505
1506    Tmpv001 =muts(i,j)*al(i,k,j)
1507    Tmpv002 =mu(i,j)*alt(i,k,j)
1508    Tmpv003 =Tmpv001 +Tmpv002
1509    Tmpv004 =dnw(k)*Tmpv003
1510    Tmpv005 =ph(i,k,j) -Tmpv004
1511    ph(i,k+1,j) =Tmpv005
1512
1513    ENDDO
1514    ENDDO
1515    ENDDO
1516    END IF
1517
1518    IF(non_hydrostatic) THEN
1519
1520    DO j =j_end, j_start, -1
1521    DO k =k_end, k_start, -1
1522    DO i =i_end, i_start, -1
1523
1524    p(i,k,j) =Tmpv406(i,k,j)
1525
1526    a_Tmpv7 =a_p(i,k,j)
1527    a_p(i,k,j) =0.0
1528    a_c2a(i,k,j) =a_c2a(i,k,j) +Tmpv405(i,k,j)*a_Tmpv7
1529    a_Tmpv6 =c2a(i,k,j)*a_Tmpv7
1530    a_Tmpv5 =a_Tmpv6
1531    a_al(i,k,j) =a_al(i,k,j) -a_Tmpv6
1532    a_Tmpv3 =a_Tmpv5/Tmpv404(i,k,j)
1533    a_Tmpv4 =-Tmpv403(i,k,j)/(Tmpv404(i,k,j)*Tmpv404(i,k,j))*a_Tmpv5
1534    a_muts(i,j) =a_muts(i,j) +(t0 +t_1(i,k,j))*a_Tmpv4
1535    a_t_1(i,k,j) =a_t_1(i,k,j) +muts(i,j)*a_Tmpv4
1536    a_alt(i,k,j) =a_alt(i,k,j) +Tmpv402(i,k,j)*a_Tmpv3
1537    a_Tmpv2 =alt(i,k,j)*a_Tmpv3
1538    a_t_2(i,k,j) =a_t_2(i,k,j) +a_Tmpv2
1539    a_Tmpv1 =-a_Tmpv2
1540    a_mu(i,j) =a_mu(i,j) +t_1(i,k,j)*a_Tmpv1
1541    a_t_1(i,k,j) =a_t_1(i,k,j) +mu(i,j)*a_Tmpv1
1542
1543    al(i,k,j) =Tmpv401(i,k,j)
1544
1545    a_Tmpv5 =a_al(i,k,j)
1546    a_al(i,k,j) =0.0
1547    a_muts(i,j) =a_muts(i,j) +1./(muts(i,j)*muts(i,j))*Tmpv400(i,k,j)*a_Tmpv5
1548    a_Tmpv4 =-1./muts(i,j)*a_Tmpv5
1549    a_Tmpv1 =a_Tmpv4
1550    a_Tmpv3 =a_Tmpv4
1551    a_Tmpv2 =rdnw(k)*a_Tmpv3
1552    a_ph(i,k+1,j) =a_ph(i,k+1,j) +a_Tmpv2
1553    a_ph(i,k,j) =a_ph(i,k,j) -a_Tmpv2
1554    a_alt(i,k,j) =a_alt(i,k,j) +mu(i,j)*a_Tmpv1
1555    a_mu(i,j) =a_mu(i,j) +alt(i,k,j)*a_Tmpv1
1556    ENDDO
1557    ENDDO
1558    ENDDO
1559
1560    ELSE
1561
1562    DO j =j_end, j_start, -1
1563    DO k =k_end, k_start, -1
1564    DO i =i_end, i_start, -1
1565    a_Tmpv5 =a_ph(i,k+1,j)
1566    a_ph(i,k+1,j) =0.0
1567    a_ph(i,k,j) =a_ph(i,k,j) +a_Tmpv5
1568    a_Tmpv4 =-a_Tmpv5
1569    a_Tmpv3 =dnw(k)*a_Tmpv4
1570    a_Tmpv1 =a_Tmpv3
1571    a_Tmpv2 =a_Tmpv3
1572    a_mu(i,j) =a_mu(i,j) +alt(i,k,j)*a_Tmpv2
1573    a_alt(i,k,j) =a_alt(i,k,j) +mu(i,j)*a_Tmpv2
1574    a_muts(i,j) =a_muts(i,j) +al(i,k,j)*a_Tmpv1
1575    a_al(i,k,j) =a_al(i,k,j) +muts(i,j)*a_Tmpv1
1576
1577    al(i,k,j) =Tmpv4011(i,k,j)
1578
1579    a_Tmpv7 =a_al(i,k,j)
1580    a_al(i,k,j) =0.0
1581    a_Tmpv5 =a_Tmpv7
1582    a_Tmpv6 =-a_Tmpv7
1583    a_p(i,k,j) =a_p(i,k,j) +a_Tmpv6/c2a(i,k,j)
1584    a_c2a(i,k,j) =a_c2a(i,k,j) -p(i,k,j)/(c2a(i,k,j)*c2a(i,k,j))*a_Tmpv6
1585    a_Tmpv3 =a_Tmpv5/Tmpv4010(i,k,j)
1586    a_Tmpv4 =-Tmpv409(i,k,j)/(Tmpv4010(i,k,j)*Tmpv4010(i,k,j))*a_Tmpv5
1587    a_muts(i,j) =a_muts(i,j) +(t0 +t_1(i,k,j))*a_Tmpv4
1588    a_t_1(i,k,j) =a_t_1(i,k,j) +muts(i,j)*a_Tmpv4
1589    a_alt(i,k,j) =a_alt(i,k,j) +Tmpv408(i,k,j)*a_Tmpv3
1590    a_Tmpv2 =alt(i,k,j)*a_Tmpv3
1591    a_t_2(i,k,j) =a_t_2(i,k,j) +a_Tmpv2
1592    a_Tmpv1 =-a_Tmpv2
1593    a_mu(i,j) =a_mu(i,j) +t_1(i,k,j)*a_Tmpv1
1594    a_t_1(i,k,j) =a_t_1(i,k,j) +mu(i,j)*a_Tmpv1
1595
1596    p(i,k,j) =Tmpv407(i,k,j)
1597
1598    a_mu(i,j) =a_mu(i,j) +znu(k)*a_p(i,k,j)
1599    a_p(i,k,j) =0.0
1600    ENDDO
1601    ENDDO
1602    ENDDO
1603
1604    END IF
1605
1606 !LPB[0]
1607 !  i_start =its
1608 !  i_end =min(ite, ide-1)
1609 !  j_start =jts
1610 !  j_end =min(jte, jde-1)
1611 !  k_start =kts
1612 !  k_end =min(kte, kde-1)
1613
1614    END SUBROUTINE a_calc_p_rho
1615
1616    SUBROUTINE a_calc_coef_w(a,a_a,alpha,a_alpha,gamma,a_gamma,mut,a_mut, &
1617    cqw,a_cqw,rdn,rdnw,c2a,a_c2a,dts,g,epssm,top_lid,ids,ide,jds,jde,kds,kde,ims, &
1618    ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte)
1619
1620 !PART I: DECLARATION OF VARIABLES
1621
1622    IMPLICIT NONE
1623
1624    INTEGER :: K0_ADJ,K1_ADJ,K2_ADJ,K3_ADJ
1625    INTEGER :: ids,ide,jds,jde,kds,kde
1626    INTEGER :: ims,ime,jms,jme,kms,kme
1627    INTEGER :: its,ite,jts,jte,kts,kte
1628    LOGICAL :: top_lid
1629    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_c2a,c2a,a_cqw,cqw
1630    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_alpha,alpha,a_gamma,gamma,a_a,a
1631    REAL,DIMENSION(ims:ime,jms:jme) :: a_mut,mut
1632    REAL,DIMENSION(kms:kme) :: rdn,rdnw
1633    REAL :: epssm,dts,g
1634    REAL,DIMENSION(ims:ime) :: a_cof,cof
1635    REAL :: a_b,b,a_c,c
1636    INTEGER :: i,j,k,i_start,i_end,j_start,j_end,k_start,k_end
1637    INTEGER :: ij,ijp,ijm,lid_flag
1638
1639    REAL :: a_Tmpv1,Tmpv001,a_Tmpv2,Tmpv002,a_Tmpv3,Tmpv003,a_Tmpv4,Tmpv004, &
1640    a_Tmpv5,Tmpv005
1641    REAL,DIMENSION(ims:ime) :: Tmpv200
1642    REAL,DIMENSION(ims:ime) :: Tmpv201
1643    REAL,DIMENSION(ims:ime) :: Tmpv202
1644    REAL,DIMENSION(ims:ime) :: Tmpv203
1645    REAL,DIMENSION(ims:ime) :: Tmpv204
1646    REAL,DIMENSION(ims:ime) :: Tmpv205
1647    REAL,DIMENSION(ims:ime) :: Tmpv206
1648    REAL,DIMENSION(ims:ime) :: Tmpv207
1649    REAL,DIMENSION(its:min(ite,ide-1),3:kde-1) :: Tmpv300
1650    REAL,DIMENSION(its:min(ite,ide-1),3:kde-1) :: Tmpv301
1651    REAL,DIMENSION(its:min(ite,ide-1),2:kde-1) :: Tmpv302
1652    REAL,DIMENSION(its:min(ite,ide-1),2:kde-1) :: Tmpv303
1653    REAL,DIMENSION(its:min(ite,ide-1),2:kde-1) :: Tmpv304
1654    REAL,DIMENSION(its:min(ite,ide-1),2:kde-1) :: Tmpv305
1655    REAL,DIMENSION(its:min(ite,ide-1),2:kde-1) :: Tmpv306
1656    REAL,DIMENSION(its:min(ite,ide-1),2:kde-1) :: Tmpv307
1657    REAL,DIMENSION(its:min(ite,ide-1),2:kde-1) :: Tmpv308  !REVISED BY WALLS
1658
1659 !PART II: CALCULATIONS OF B. S. TRAJECTORY
1660
1661 !LPB[0]
1662
1663          i_start = its
1664          i_end   = min(ite,ide-1)
1665          j_start = jts
1666          j_end   = min(jte,jde-1)
1667          k_start = kts
1668          k_end   = kte-1
1669          lid_flag=1
1670
1671 !LPB[1]
1672       IF(top_lid)lid_flag=0
1673
1674 !!LPB[2]
1675 !        outer_j_loop:  DO j = j_start, j_end
1676
1677 !           DO i = i_start, i_end
1678 !             cof(i)  = (.5*dts*g*(1.+epssm)/mut(i,j))**2
1679 !             a(i, 2 ,j) = 0.
1680 !             a(i,kde,j) =-2.*cof(i)*rdnw(kde-1)**2*c2a(i,kde-1,j)*lid_flag
1681 !             gamma(i,1 ,j) = 0.
1682 !           ENDDO
1683
1684 !           DO k=3,kde-1
1685 !           DO i=i_start, i_end
1686 !             a(i,k,j) =   -cqw(i,k,j)*cof(i)*rdn(k)* rdnw(k-1)*c2a(i,k-1,j)
1687 !           ENDDO
1688 !           ENDDO
1689
1690 !           DO k=2,kde-1
1691 !           DO i=i_start, i_end
1692 !             b = 1.+cqw(i,k,j)*cof(i)*rdn(k)*(rdnw(k  )*c2a(i,k,j  )    &
1693 !                                          +rdnw(k-1)*c2a(i,k-1,j))
1694 !             c =   -cqw(i,k,j)*cof(i)*rdn(k)*rdnw(k  )*c2a(i,k,j  )
1695 !             alpha(i,k,j) = 1./(b-a(i,k,j)*gamma(i,k-1,j))
1696 !             gamma(i,k,j) = c*alpha(i,k,j)
1697 !           ENDDO
1698 !           ENDDO
1699
1700 !           DO i=i_start, i_end
1701 !             b = 1.+2.*cof(i)*rdnw(kde-1)**2*c2a(i,kde-1,j)
1702 !             c = 0.
1703 !             alpha(i,kde,j) = 1./(b-a(i,kde,j)*gamma(i,kde-1,j))
1704 !             gamma(i,kde,j) = c*alpha(i,kde,j)
1705 !           ENDDO
1706
1707 !         ENDDO outer_j_loop
1708
1709 !PART III: INITIALIZATION OF LOCAL ADJOINT PERTURBATIONS
1710
1711    Do K0_ADJ =ims, ime
1712    a_cof(K0_ADJ) =0.0
1713    End Do
1714
1715    a_b =0.0
1716    a_c =0.0
1717
1718 !PART IV: REVERSE/BACKWARD ACCUMULATIONS
1719
1720 !LPB[2]
1721    DO j =j_end, j_start, -1
1722
1723    DO i =i_start, i_end
1724    Tmpv200(i) =cof(i)
1725    cof(i) =(.5*dts*g*(1.+epssm)/mut(i,j))**2
1726
1727    Tmpv201(i) =a(i,2,j)
1728    a(i,2,j) =0.
1729
1730    Tmpv001 =-2.*cof(i)*rdnw(kde-1)**2*c2a(i,kde-1,j)
1731    Tmpv002 =Tmpv001*lid_flag
1732    Tmpv202(i) =a(i,kde,j)
1733    a(i,kde,j) =Tmpv002
1734
1735    Tmpv203(i) =gamma(i,1,j)
1736    gamma(i,1,j) =0.
1737
1738    ENDDO
1739
1740    DO k =3, kde-1
1741    DO i =i_start, i_end
1742    Tmpv001 =-cqw(i,k,j)*cof(i)
1743    Tmpv002 =Tmpv001*rdn(k)
1744    Tmpv003 =Tmpv002*rdnw(k-1)
1745    Tmpv300(i,k) =Tmpv003
1746    Tmpv004 =Tmpv300(i,k)*c2a(i,k-1,j)
1747    Tmpv301(i,k) =a(i,k,j)
1748    a(i,k,j) =Tmpv004
1749
1750    ENDDO
1751    ENDDO
1752    DO k =2, kde-1
1753    DO i =i_start, i_end
1754    Tmpv001 =cqw(i,k,j)*cof(i)
1755    Tmpv002 =Tmpv001*rdn(k)
1756    Tmpv003 =rdnw(k)*c2a(i,k,j) +rdnw(k-1)*c2a(i,k-1,j)
1757    Tmpv302(i,k) =Tmpv002
1758    Tmpv303(i,k) =Tmpv003
1759    Tmpv004 =Tmpv302(i,k)*Tmpv303(i,k)
1760    Tmpv005 =1. +Tmpv004
1761    b =Tmpv005
1762
1763    Tmpv001 =-cqw(i,k,j)*cof(i)
1764    Tmpv002 =Tmpv001*rdn(k)
1765    Tmpv003 =Tmpv002*rdnw(k)
1766    Tmpv304(i,k) =Tmpv003
1767    Tmpv004 =Tmpv304(i,k)*c2a(i,k,j)
1768 !  Tmpv305(i,k) =c
1769    c =Tmpv004
1770
1771    Tmpv305(i,k) =c  !REVISED BY WALLS
1772
1773    Tmpv001 =a(i,k,j)*gamma(i,k-1,j)
1774    Tmpv002 =b -Tmpv001
1775    Tmpv308(i,k) =Tmpv002  !REVISED BY WALLS
1776    Tmpv003 =1./Tmpv002
1777    Tmpv306(i,k) =alpha(i,k,j)
1778    alpha(i,k,j) =Tmpv003
1779
1780    Tmpv001 =c*alpha(i,k,j)
1781 !  Tmpv307(i,k) =gamma(i,k,j)
1782    Tmpv307(i,k) =gamma(i,k-1,j)
1783    gamma(i,k,j) =Tmpv001
1784
1785    ENDDO
1786    ENDDO
1787    DO i =i_start, i_end
1788    Tmpv001 =2.*cof(i)*rdnw(kde-1)**2*c2a(i,kde-1,j)
1789    Tmpv002 =1. +Tmpv001
1790    b =Tmpv002
1791
1792    Tmpv204(i) =c
1793    c =0.
1794
1795    Tmpv001 =a(i,kde,j)*gamma(i,kde-1,j)
1796    Tmpv002 =b -Tmpv001
1797    Tmpv207(i) = Tmpv002  ! Added by Ning Pan, 2010-08-02
1798    Tmpv003 =1./Tmpv002
1799    Tmpv205(i) =alpha(i,kde,j)
1800    alpha(i,kde,j) =Tmpv003
1801
1802    Tmpv001 =c*alpha(i,kde,j)
1803    Tmpv206(i) =gamma(i,kde,j)
1804    gamma(i,kde,j) =Tmpv001
1805
1806    ENDDO
1807
1808    DO i =i_end, i_start, -1
1809
1810    gamma(i,kde,j) =Tmpv206(i)
1811
1812    a_Tmpv1 =a_gamma(i,kde,j)
1813    a_gamma(i,kde,j) =0.0
1814    a_c =a_c +alpha(i,kde,j)*a_Tmpv1
1815    a_alpha(i,kde,j) =a_alpha(i,kde,j) +c*a_Tmpv1
1816
1817    alpha(i,kde,j) =Tmpv205(i)
1818
1819    a_Tmpv3 =a_alpha(i,kde,j)
1820    a_alpha(i,kde,j) =0.0
1821    Tmpv002 = Tmpv207(i)  ! Added by Ning Pan, 2010-08-02
1822    a_Tmpv2 =-(1.)*a_Tmpv3/(Tmpv002*Tmpv002)
1823    a_b =a_b +a_Tmpv2
1824    a_Tmpv1 =-a_Tmpv2
1825    a_a(i,kde,j) =a_a(i,kde,j) +gamma(i,kde-1,j)*a_Tmpv1
1826    a_gamma(i,kde-1,j) =a_gamma(i,kde-1,j) +a(i,kde,j)*a_Tmpv1
1827
1828    c =Tmpv204(i)
1829
1830    a_c =0.0
1831    a_Tmpv2 =a_b
1832    a_b =0.0
1833    a_Tmpv1 =a_Tmpv2
1834    a_cof(i) =a_cof(i) +2.*rdnw(kde-1)**2*c2a(i,kde-1,j)*a_Tmpv1
1835    a_c2a(i,kde-1,j) =a_c2a(i,kde-1,j) +2.*cof(i)*rdnw(kde-1)**2*a_Tmpv1
1836    ENDDO
1837
1838    DO k =kde-1, 2, -1
1839    DO i =i_end, i_start, -1
1840
1841    c =Tmpv305(i,k)  !REVISED BY WALLS
1842
1843 !  gamma(i,k,j) =Tmpv307(i,k)
1844    gamma(i,k-1,j) =Tmpv307(i,k)
1845
1846    a_Tmpv1 =a_gamma(i,k,j)
1847    a_gamma(i,k,j) =0.0
1848    a_c =a_c +alpha(i,k,j)*a_Tmpv1
1849    a_alpha(i,k,j) =a_alpha(i,k,j) +c*a_Tmpv1
1850
1851    alpha(i,k,j) =Tmpv306(i,k)
1852
1853    a_Tmpv3 =a_alpha(i,k,j)
1854    a_alpha(i,k,j) =0.0
1855    Tmpv002 =Tmpv308(i,k) !REVISED BY WALLS
1856    a_Tmpv2 =-(1.)*a_Tmpv3/(Tmpv002*Tmpv002)
1857    a_b =a_b +a_Tmpv2
1858    a_Tmpv1 =-a_Tmpv2
1859    a_a(i,k,j) =a_a(i,k,j) +gamma(i,k-1,j)*a_Tmpv1
1860    a_gamma(i,k-1,j) =a_gamma(i,k-1,j) +a(i,k,j)*a_Tmpv1
1861
1862 !  c =Tmpv305(i,k)
1863
1864    a_Tmpv4 =a_c
1865    a_c =0.0
1866    a_Tmpv3 =c2a(i,k,j)*a_Tmpv4
1867    a_c2a(i,k,j) =a_c2a(i,k,j) +Tmpv304(i,k)*a_Tmpv4
1868    a_Tmpv2 =rdnw(k)*a_Tmpv3
1869    a_Tmpv1 =rdn(k)*a_Tmpv2
1870    a_cqw(i,k,j) =a_cqw(i,k,j) -cof(i)*a_Tmpv1
1871    a_cof(i) =a_cof(i) -cqw(i,k,j)*a_Tmpv1
1872    a_Tmpv5 =a_b
1873    a_b =0.0
1874    a_Tmpv4 =a_Tmpv5
1875    a_Tmpv2 =Tmpv303(i,k)*a_Tmpv4
1876    a_Tmpv3 =Tmpv302(i,k)*a_Tmpv4
1877    a_c2a(i,k,j) =a_c2a(i,k,j) +rdnw(k)*a_Tmpv3
1878    a_c2a(i,k-1,j) =a_c2a(i,k-1,j) +rdnw(k-1)*a_Tmpv3
1879    a_Tmpv1 =rdn(k)*a_Tmpv2
1880    a_cqw(i,k,j) =a_cqw(i,k,j) +cof(i)*a_Tmpv1
1881    a_cof(i) =a_cof(i) +cqw(i,k,j)*a_Tmpv1
1882    ENDDO
1883    ENDDO
1884
1885    DO k =kde-1, 3, -1
1886    DO i =i_end, i_start, -1
1887
1888    a(i,k,j) =Tmpv301(i,k)
1889
1890    a_Tmpv4 =a_a(i,k,j)
1891    a_a(i,k,j) =0.0
1892    a_Tmpv3 =c2a(i,k-1,j)*a_Tmpv4
1893    a_c2a(i,k-1,j) =a_c2a(i,k-1,j) +Tmpv300(i,k)*a_Tmpv4
1894    a_Tmpv2 =rdnw(k-1)*a_Tmpv3
1895    a_Tmpv1 =rdn(k)*a_Tmpv2
1896    a_cqw(i,k,j) =a_cqw(i,k,j) -cof(i)*a_Tmpv1
1897    a_cof(i) =a_cof(i) -cqw(i,k,j)*a_Tmpv1
1898    ENDDO
1899    ENDDO
1900
1901    DO i =i_end, i_start, -1
1902
1903    gamma(i,1,j) =Tmpv203(i)
1904
1905    a_gamma(i,1,j) =0.0
1906
1907    a(i,kde,j) =Tmpv202(i)
1908
1909    a_Tmpv2 =a_a(i,kde,j)
1910    a_a(i,kde,j) =0.0
1911    a_Tmpv1 =lid_flag*a_Tmpv2
1912    a_cof(i) =a_cof(i) -2.*rdnw(kde-1)**2*c2a(i,kde-1,j)*a_Tmpv1
1913    a_c2a(i,kde-1,j) =a_c2a(i,kde-1,j) -2.*cof(i)*rdnw(kde-1)**2*a_Tmpv1
1914
1915    a(i,2,j) =Tmpv201(i)
1916
1917    a_a(i,2,j) =0.0
1918
1919    cof(i) =Tmpv200(i)
1920
1921    a_mut(i,j) =a_mut(i,j) -.5*dts*g*(1.+epssm)/(mut(i,j)*mut(i,j))*2.0*(.5*dts*g*  &
1922    (1.+epssm)/mut(i,j))*a_cof(i)
1923    a_cof(i) =0.0
1924    ENDDO
1925
1926    ENDDO
1927
1928 !LPB[1]
1929
1930 !  IF(top_lid) THEN
1931 !  lid_flag =0
1932 !  END IF
1933
1934 !  IF(top_lid) THEN
1935
1936 !  END IF
1937
1938 !LPB[0]
1939 !  i_start =its
1940 !  i_end =min(ite, ide-1)
1941 !  j_start =jts
1942 !  j_end =min(jte, jde-1)
1943 !  k_start =kts
1944 !  k_end =kte-1
1945 !  lid_flag =1
1946
1947    END SUBROUTINE a_calc_coef_w
1948
1949    SUBROUTINE a_advance_uv(u,a_u,ru_tend,a_ru_tend,v,a_v,rv_tend, &
1950    a_rv_tend,p,a_p,pb,ph,a_ph,php,a_php,alt,a_alt,al,a_al,mu,a_mu, &
1951    muu,a_muu,cqu,a_cqu,muv,a_muv,cqv,a_cqv,mudf,a_mudf,msfux,msfuy,msfvx, &
1952    msfvx_inv,msfvy,rdx,rdy,dts,cf1,cf2,cf3,fnm,fnp,emdiv,rdnw,config_flags,spec_zone, &
1953    non_hydrostatic,top_lid,ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts, &
1954    jte,kts,kte)
1955
1956 !PART I: DECLARATION OF VARIABLES
1957
1958    IMPLICIT NONE
1959
1960    INTEGER :: K0_ADJ,K1_ADJ,K2_ADJ,K3_ADJ
1961    TYPE(grid_config_rec_type) :: config_flags
1962    LOGICAL :: non_hydrostatic,top_lid
1963    INTEGER :: ids,ide,jds,jde,kds,kde
1964    INTEGER :: ims,ime,jms,jme,kms,kme
1965    INTEGER :: its,ite,jts,jte,kts,kte
1966    INTEGER :: spec_zone
1967    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_u,u,a_v,v
1968    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_ru_tend,ru_tend,a_rv_tend,rv_tend, &
1969    a_ph,ph,a_php,php,a_p,p,pb,a_alt,alt,a_al,al,a_cqu,cqu,a_cqv,cqv
1970    REAL,DIMENSION(ims:ime,jms:jme) :: a_muu,muu,a_muv,muv,a_mu,mu,a_mudf,mudf
1971    REAL,DIMENSION(kms:kme) :: fnm,fnp,rdnw
1972    REAL,DIMENSION(ims:ime,jms:jme) :: msfux,msfuy,msfvx,msfvy,msfvx_inv
1973    REAL :: rdx,rdy,dts,cf1,cf2,cf3,emdiv
1974    REAL,DIMENSION(its:ite,kts:kte) :: a_dpn,dpn,a_dpxy,dpxy
1975    REAL,DIMENSION(its:ite) :: a_mudf_xy,mudf_xy
1976    REAL :: dx,dy
1977    INTEGER :: i,j,k,i_start,i_end,j_start,j_end,k_start,k_end
1978    INTEGER :: i_endu,j_endv,k_endw
1979    INTEGER :: i_start_up,i_end_up,j_start_up,j_end_up
1980    INTEGER :: i_start_vp,i_end_vp,j_start_vp,j_end_vp
1981    INTEGER :: i_start_u_tend,i_end_u_tend,j_start_v_tend,j_end_v_tend
1982
1983 !  REAL,DIMENSION(its:ite,kts:kte) :: Keep_Lpb20_dpxy
1984 !  REAL,DIMENSION(its:ite,kts:kte) :: Keep_Lpb20_dpn
1985    REAL :: a_Tmpv1,Tmpv001,a_Tmpv2,Tmpv002,a_Tmpv3,Tmpv003,a_Tmpv4,Tmpv004, &
1986    a_Tmpv5,Tmpv005,a_Tmpv6,Tmpv006,a_Tmpv7,Tmpv007,a_Tmpv8,Tmpv008,a_Tmpv9, &
1987    Tmpv009,a_Tmpv10,Tmpv010,a_Tmpv11,Tmpv011
1988    REAL,DIMENSION(its:ite,kts:kte) :: Tmpv300
1989    REAL,DIMENSION(its:ite,kts:kte) :: Tmpv301
1990    REAL,DIMENSION(its:ite,kts:kte) :: Tmpv302
1991    REAL,DIMENSION(its:ite,kts:kte) :: Tmpv303
1992    REAL,DIMENSION(its:ite,kts:kte) :: Tmpv304
1993    REAL,DIMENSION(its:ite,kts:kte) :: Tmpv305
1994    REAL,DIMENSION(its:ite,kts:kte) :: Tmpv306
1995
1996 !PART II: CALCULATIONS OF B. S. TRAJECTORY
1997
1998 !LPB[0]
1999
2000 !LPB[1]
2001     IF( config_flags%nested .or. config_flags%specified ) THEN
2002
2003          i_start = max( its,ids+spec_zone )
2004          i_end   = min( ite,ide-spec_zone-1 )
2005          j_start = max( jts,jds+spec_zone )
2006          j_end   = min( jte,jde-spec_zone-1 )
2007          k_start = kts
2008          k_end   = min( kte, kde-1 )
2009          i_endu = min( ite,ide-spec_zone )
2010          j_endv = min( jte,jde-spec_zone )
2011          k_endw = k_end
2012       IF( config_flags%periodic_x) THEN
2013
2014            i_start = its
2015            i_end   = min(ite,ide-1)
2016            i_endu = ite
2017          ENDIF
2018        ELSE
2019          i_start = its
2020          i_end   = min(ite,ide-1)
2021          j_start = jts
2022          j_end   = min(jte,jde-1)
2023          k_start = kts
2024          k_end   = kte-1
2025          i_endu = ite
2026          j_endv = jte
2027          k_endw = k_end
2028
2029    ENDIF
2030
2031 !LPB[2]
2032
2033          i_start_up = i_start
2034          i_end_up   = i_endu
2035          j_start_up = j_start
2036          j_end_up   = j_end
2037          i_start_vp = i_start
2038          i_end_vp   = i_end
2039          j_start_vp = j_start
2040          j_end_vp   = j_endv
2041
2042 !LPB[3]
2043       IF ( (config_flags%open_xs   .or.       &
2044             config_flags%symmetric_xs   )     &
2045             .and. (its == ids) )              &
2046                  i_start_up = i_start_up + 1
2047
2048 !LPB[4]
2049
2050 !LPB[5]
2051       IF ( (config_flags%open_xe    .or.    &
2052             config_flags%symmetric_xe   )   &
2053              .and. (ite == ide) )           &
2054                  i_end_up   = i_end_up - 1
2055
2056 !LPB[6]
2057
2058 !LPB[7]
2059       IF ( (config_flags%open_ys    .or.     &
2060             config_flags%symmetric_ys  .or.     &
2061             config_flags%polar   )    &
2062                      .and. (jts == jds) )    &
2063                  j_start_vp = j_start_vp + 1
2064
2065 !LPB[8]
2066
2067 !LPB[9]
2068       IF ( (config_flags%open_ye     .or.   &
2069             config_flags%symmetric_ye  .or.     &
2070             config_flags%polar   )    &
2071             .and. (jte == jde) )            &
2072                  j_end_vp   = j_end_vp - 1
2073
2074 !LPB[10]
2075          i_start_u_tend = i_start
2076          i_end_u_tend   = i_endu
2077          j_start_v_tend = j_start
2078          j_end_v_tend   = j_endv
2079
2080 !LPB[11]
2081       IF ( config_flags%symmetric_xs .and. (its == ids) )   &
2082                      i_start_u_tend = i_start_u_tend+1
2083
2084 !LPB[12]
2085
2086 !LPB[13]
2087       IF ( config_flags%symmetric_xe .and. (ite == ide) )   &
2088                      i_end_u_tend = i_end_u_tend-1
2089
2090 !LPB[14]
2091
2092 !LPB[15]
2093       IF ( config_flags%symmetric_ys .and. (jts == jds) )   &
2094                      j_start_v_tend = j_start_v_tend+1
2095
2096 !LPB[16]
2097
2098 !LPB[17]
2099       IF ( config_flags%symmetric_ye .and. (jte == jde) )   &
2100                      j_end_v_tend = j_end_v_tend-1
2101
2102 !LPB[18]
2103       dx = 1./rdx
2104       dy = 1./rdy
2105
2106 !LPB[19]
2107      u_outer_j_loop: DO j = j_start, j_end
2108
2109       DO k = k_start, k_end
2110       DO i = i_start_u_tend, i_end_u_tend
2111         u(i,k,j) = u(i,k,j) + dts*ru_tend(i,k,j)
2112       ENDDO
2113       ENDDO
2114
2115       DO i = i_start_up, i_end_up
2116         mudf_xy(i)= -emdiv*dx*(mudf(i,j)-mudf(i-1,j))/msfuy(i,j)
2117       ENDDO
2118
2119       DO k = k_start, k_end
2120       DO i = i_start_up, i_end_up
2121         dpxy(i,k)= (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*(           &
2122           ((ph (i,k+1,j)-ph (i-1,k+1,j))+(ph (i,k,j)-ph (i-1,k,j)))    &
2123          +(alt(i,k  ,j)+alt(i-1,k  ,j))*(p  (i,k,j)-p  (i-1,k,j))    &
2124          +(al (i,k  ,j)+al (i-1,k  ,j))*(pb (i,k,j)-pb (i-1,k,j)) )
2125       ENDDO
2126       ENDDO
2127    IF (non_hydrostatic) THEN
2128
2129         DO i = i_start_up, i_end_up
2130           dpn(i,1) = .5*( cf1*(p(i,1,j)+p(i-1,1,j))    &
2131                          +cf2*(p(i,2,j)+p(i-1,2,j))    &
2132                          +cf3*(p(i,3,j)+p(i-1,3,j)) )
2133           dpn(i,kde) = 0.
2134         ENDDO
2135      IF (top_lid) THEN
2136
2137           DO i = i_start_up, i_end_up
2138             dpn(i,kde) =.5*( cf1*(p(i-1,kde-1,j)+p(i,kde-1,j))     &
2139                             +cf2*(p(i-1,kde-2,j)+p(i,kde-2,j))     &
2140                             +cf3*(p(i-1,kde-3,j)+p(i,kde-3,j))  )
2141           ENDDO
2142         ENDIF
2143
2144         DO k = k_start+1, k_end
2145           DO i = i_start_up, i_end_up
2146             dpn(i,k) = .5*( fnm(k)*(p(i,k  ,j)+p(i-1,k  ,j))     &
2147                            +fnp(k)*(p(i,k-1,j)+p(i-1,k-1,j)) )
2148           ENDDO
2149         ENDDO
2150
2151         DO k = k_start, k_end
2152           DO i = i_start_up, i_end_up
2153             dpxy(i,k)=dpxy(i,k) + (msfux(i,j)/msfuy(i,j))*rdx*(php(i,k,j)-php(i-1,k,j)) &
2154    *          &
2155               (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i-1,j)+mu(i,j)))
2156           ENDDO
2157         ENDDO
2158       END IF
2159
2160       DO k = k_start, k_end
2161         DO i = i_start_up, i_end_up
2162           u(i,k,j)=u(i,k,j)-dts*cqu(i,k,j)*dpxy(i,k)+mudf_xy(i)
2163         ENDDO
2164       ENDDO
2165
2166       ENDDO u_outer_j_loop
2167
2168 !!LPB[20]
2169 !     v_outer_j_loop: DO j = j_start_v_tend, j_end_v_tend
2170
2171 !!      DO k=k_start, k_end
2172 !    !  Keep_Lpb20_dpxy(i,k) =dpxy(i,k)
2173 !!      END DO
2174 !!      DO k=k_start+1, k_end
2175 !    !  Keep_Lpb20_dpn(i,k) =dpn(i,k)
2176 !!      END DO
2177
2178 !      DO k = k_start, k_end
2179 !      DO i = i_start, i_end
2180 !        v(i,k,j) = v(i,k,j) + dts*rv_tend(i,k,j)
2181 !      ENDDO
2182 !      ENDDO
2183
2184 !      DO i = i_start, i_end
2185 !        mudf_xy(i)= -emdiv*dy*(mudf(i,j)-mudf(i,j-1))*msfvx_inv(i,j)
2186 !      ENDDO
2187 !     IF (     ( j >= j_start_vp)    &
2188 !       .and.( j <= j_end_vp  ) )  THEN
2189
2190 !            DO k = k_start, k_end
2191 !            DO i = i_start, i_end
2192 !          dpxy(i,k)= (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*(         &
2193 !           ((ph(i,k+1,j)-ph(i,k+1,j-1))+(ph (i,k,j)-ph (i,k,j-1)))     &
2194 !           +(alt(i,k  ,j)+alt(i,k  ,j-1))*(p  (i,k,j)-p  (i,k,j-1))    &
2195 !           +(al (i,k  ,j)+al (i,k  ,j-1))*(pb (i,k,j)-pb (i,k,j-1)) )
2196 !               ENDDO
2197 !               ENDDO
2198 !     IF (non_hydrostatic) THEN
2199
2200 !          DO i = i_start, i_end
2201 !            dpn(i,1) = .5*( cf1*(p(i,1,j)+p(i,1,j-1))    &
2202 !                           +cf2*(p(i,2,j)+p(i,2,j-1))    &
2203 !                           +cf3*(p(i,3,j)+p(i,3,j-1)) )
2204 !            dpn(i,kde) = 0.
2205 !          ENDDO
2206 !     IF (top_lid) THEN
2207
2208 !          DO i = i_start, i_end
2209 !            dpn(i,kde) =.5*( cf1*(p(i,kde-1,j-1)+p(i,kde-1,j))     &
2210 !                            +cf2*(p(i,kde-2,j-1)+p(i,kde-2,j))     &
2211 !                            +cf3*(p(i,kde-3,j-1)+p(i,kde-3,j))  )
2212 !          ENDDO
2213 !        ENDIF
2214
2215 !          DO k = k_start+1, k_end
2216 !            DO i = i_start, i_end
2217 !              dpn(i,k) = .5*( fnm(k)*(p(i,k  ,j)+p(i,k  ,j-1))    &
2218 !                             +fnp(k)*(p(i,k-1,j)+p(i,k-1,j-1)) )
2219 !            ENDDO
2220 !          ENDDO
2221
2222 !          DO k = k_start, k_end
2223 !            DO i = i_start, i_end
2224 !              dpxy(i,k)=dpxy(i,k) + (msfvy(i,j)/msfvx(i,j))*rdy*(php(i,k,j) &
2225 !   -php(i,k,j-1))*      &
2226 !                (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i,j-1)+mu(i,j)))
2227 !            ENDDO
2228 !          ENDDO
2229 !        END IF
2230
2231 !        DO k = k_start, k_end
2232 !          DO i = i_start, i_end
2233 !            v(i,k,j)=v(i,k,j)-dts*cqv(i,k,j)*dpxy(i,k)+mudf_xy(i)
2234 !          ENDDO
2235 !        ENDDO
2236 !      END IF
2237
2238 !     ENDDO  v_outer_j_loop
2239
2240 !!LPB[21]
2241
2242 !!LPB[22]
2243 !  IF (config_flags%polar) THEN
2244
2245 !     IF (jts == jds) THEN
2246
2247 !           DO k = k_start, k_end
2248 !           DO i = i_start, i_end
2249 !              v(i,k,jds) = 0.
2250 !           ENDDO
2251 !           ENDDO
2252 !        END IF
2253 !     IF (jte == jde) THEN
2254
2255 !           DO k = k_start, k_end
2256 !           DO i = i_start, i_end
2257 !              v(i,k,jde) = 0.
2258 !           ENDDO
2259 !           ENDDO
2260 !        END IF
2261
2262 !   END IF
2263
2264 !PART III: INITIALIZATION OF LOCAL ADJOINT PERTURBATIONS
2265
2266    Do K1_ADJ =kts, kte
2267    Do K0_ADJ =its, ite
2268    a_dpn(K0_ADJ,K1_ADJ) =0.0
2269    End Do
2270    End Do
2271
2272    Do K1_ADJ =kts, kte
2273    Do K0_ADJ =its, ite
2274    a_dpxy(K0_ADJ,K1_ADJ) =0.0
2275    End Do
2276    End Do
2277
2278    Do K0_ADJ =its, ite
2279    a_mudf_xy(K0_ADJ) =0.0
2280    End Do
2281
2282 !PART IV: REVERSE/BACKWARD ACCUMULATIONS
2283
2284 !LPB[22]
2285
2286 !  IF(config_flags%polar) THEN
2287 !  IF(jts == jds) THEN
2288 !  DO k =k_start, k_end
2289 !  DO i =i_start, i_end
2290 !  v(i,k,jds) =0.
2291
2292 !  ENDDO
2293 !  ENDDO
2294 !  END IF
2295 !  IF(jte == jde) THEN
2296 !  DO k =k_start, k_end
2297 !  DO i =i_start, i_end
2298 !  v(i,k,jde) =0.
2299
2300 !  ENDDO
2301 !  ENDDO
2302 !  END IF
2303 !  END IF
2304
2305    IF(config_flags%polar) THEN
2306
2307    IF(jte == jde) THEN
2308
2309    DO k =k_end, k_start, -1
2310    DO i =i_end, i_start, -1
2311    a_v(i,k,jde) =0.0
2312    ENDDO
2313    ENDDO
2314
2315    END IF
2316
2317    IF(jts == jds) THEN
2318
2319    DO k =k_end, k_start, -1
2320    DO i =i_end, i_start, -1
2321    a_v(i,k,jds) =0.0
2322    ENDDO
2323    ENDDO
2324
2325    END IF
2326
2327    END IF
2328
2329 !LPB[21]
2330
2331 !LPB[20]
2332    DO j =j_end_v_tend, j_start_v_tend, -1
2333
2334 !  DO k=k_start, k_end
2335 !  dpxy(i,k) =Keep_Lpb20_dpxy(i,k)
2336 !  END DO
2337 !  DO k=k_start+1, k_end
2338 !  dpn(i,k) =Keep_Lpb20_dpn(i,k)
2339 !  END DO
2340
2341    DO k =k_start, k_end
2342    DO i =i_start, i_end
2343    Tmpv001 =v(i,k,j) +dts*rv_tend(i,k,j)
2344 !  v(i,k,j) =Tmpv001
2345
2346    ENDDO
2347    ENDDO
2348    DO i =i_start, i_end
2349    Tmpv001 =mudf(i,j) -mudf(i,j-1)
2350    Tmpv002 =-emdiv*dy*Tmpv001
2351    Tmpv003 =Tmpv002*msfvx_inv(i,j)
2352 !  mudf_xy(i) =Tmpv003
2353
2354    ENDDO
2355
2356    IF(     ( j >= j_start_vp)                  .and.( j <= j_end_vp  ) ) THEN
2357    DO k =k_start, k_end
2358    DO i =i_start, i_end
2359    Tmpv001 =ph(i,k+1,j) -ph(i,k+1,j-1)
2360    Tmpv002 =ph(i,k,j) -ph(i,k,j-1)
2361    Tmpv003 =Tmpv001 +Tmpv002
2362    Tmpv004 =alt(i,k,j) +alt(i,k,j-1)
2363    Tmpv005 =p(i,k,j) -p(i,k,j-1)
2364    Tmpv300(i,k) =Tmpv004
2365    Tmpv301(i,k) =Tmpv005
2366    Tmpv006 =Tmpv300(i,k)*Tmpv301(i,k)
2367    Tmpv007 =Tmpv003 +Tmpv006
2368    Tmpv008 =al(i,k,j) +al(i,k,j-1)
2369    Tmpv009 =Tmpv008*(pb(i,k,j)-pb(i,k,j-1))
2370    Tmpv010 =Tmpv007 +Tmpv009
2371    Tmpv302(i,k) =Tmpv010
2372    Tmpv011 =(msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*Tmpv302(i,k)
2373    Tmpv303(i,k) =dpxy(i,k)
2374    dpxy(i,k) =Tmpv011
2375
2376    ENDDO
2377    ENDDO
2378    IF(non_hydrostatic) THEN
2379    DO i =i_start, i_end
2380    Tmpv001 =p(i,1,j) +p(i,1,j-1)
2381    Tmpv002 =cf1*Tmpv001
2382    Tmpv003 =p(i,2,j) +p(i,2,j-1)
2383    Tmpv004 =cf2*Tmpv003
2384    Tmpv005 =Tmpv002 +Tmpv004
2385    Tmpv006 =p(i,3,j) +p(i,3,j-1)
2386    Tmpv007 =cf3*Tmpv006
2387    Tmpv008 =Tmpv005 +Tmpv007
2388    Tmpv009 =.5*Tmpv008
2389    dpn(i,1) =Tmpv009
2390
2391    dpn(i,kde) =0.
2392
2393    ENDDO
2394
2395    IF(top_lid) THEN
2396    DO i =i_start, i_end
2397    Tmpv001 =p(i,kde-1,j-1) +p(i,kde-1,j)
2398    Tmpv002 =cf1*Tmpv001
2399    Tmpv003 =p(i,kde-2,j-1) +p(i,kde-2,j)
2400    Tmpv004 =cf2*Tmpv003
2401    Tmpv005 =Tmpv002 +Tmpv004
2402    Tmpv006 =p(i,kde-3,j-1) +p(i,kde-3,j)
2403    Tmpv007 =cf3*Tmpv006
2404    Tmpv008 =Tmpv005 +Tmpv007
2405    Tmpv009 =.5*Tmpv008
2406    dpn(i,kde) =Tmpv009
2407
2408    ENDDO
2409
2410    ENDIF
2411    DO k =k_start+1, k_end
2412    DO i =i_start, i_end
2413    Tmpv001 =p(i,k,j) +p(i,k,j-1)
2414    Tmpv002 =fnm(k)*Tmpv001
2415    Tmpv003 =p(i,k-1,j) +p(i,k-1,j-1)
2416    Tmpv004 =fnp(k)*Tmpv003
2417    Tmpv005 =Tmpv002 +Tmpv004
2418    Tmpv006 =.5*Tmpv005
2419    dpn(i,k) =Tmpv006
2420
2421    ENDDO
2422    ENDDO
2423    DO k =k_start, k_end
2424    DO i =i_start, i_end
2425    Tmpv001 =php(i,k,j) -php(i,k,j-1)
2426    Tmpv002 =(msfvy(i,j)/msfvx(i,j))*rdy*Tmpv001
2427    Tmpv003 =dpn(i,k+1) -dpn(i,k)
2428    Tmpv004 =rdnw(k)*Tmpv003
2429    Tmpv005 =mu(i,j-1) +mu(i,j)
2430    Tmpv006 =.5*Tmpv005
2431    Tmpv007 =Tmpv004 -Tmpv006
2432    Tmpv304(i,k) =Tmpv002
2433    Tmpv305(i,k) =Tmpv007
2434    Tmpv008 =Tmpv304(i,k)*Tmpv305(i,k)
2435    Tmpv009 =dpxy(i,k) +Tmpv008
2436    Tmpv306(i,k) =dpxy(i,k)
2437    dpxy(i,k) =Tmpv009
2438
2439    ENDDO
2440    ENDDO
2441    END IF
2442    DO k =k_start, k_end
2443    DO i =i_start, i_end
2444    Tmpv001 =dts*cqv(i,k,j)*dpxy(i,k)
2445    Tmpv002 =v(i,k,j) -Tmpv001
2446    Tmpv003 =Tmpv002 +mudf_xy(i)
2447 !  v(i,k,j) =Tmpv003
2448
2449    ENDDO
2450    ENDDO
2451    END IF
2452
2453    IF(     ( j >= j_start_vp)    &
2454           .and.( j <= j_end_vp  ) ) THEN
2455
2456    DO k =k_end, k_start, -1
2457    DO i =i_end, i_start, -1
2458    a_Tmpv3 =a_v(i,k,j)
2459    a_v(i,k,j) =0.0
2460    a_Tmpv2 =a_Tmpv3
2461    a_mudf_xy(i) =a_mudf_xy(i) +a_Tmpv3
2462    a_v(i,k,j) =a_v(i,k,j) +a_Tmpv2
2463    a_Tmpv1 =-a_Tmpv2
2464    a_cqv(i,k,j) =a_cqv(i,k,j) +dts*dpxy(i,k)*a_Tmpv1
2465    a_dpxy(i,k) =a_dpxy(i,k) +dts*cqv(i,k,j)*a_Tmpv1
2466    ENDDO
2467    ENDDO
2468
2469    IF(non_hydrostatic) THEN
2470
2471    DO k =k_end, k_start, -1
2472    DO i =i_end, i_start, -1
2473
2474    dpxy(i,k) =Tmpv306(i,k)
2475
2476    a_Tmpv9 =a_dpxy(i,k)
2477    a_dpxy(i,k) =0.0
2478    a_dpxy(i,k) =a_dpxy(i,k) +a_Tmpv9
2479    a_Tmpv8 =a_Tmpv9
2480    a_Tmpv2 =Tmpv305(i,k)*a_Tmpv8
2481    a_Tmpv7 =Tmpv304(i,k)*a_Tmpv8
2482    a_Tmpv4 =a_Tmpv7
2483    a_Tmpv6 =-a_Tmpv7
2484    a_Tmpv5 =.5*a_Tmpv6
2485    a_mu(i,j-1) =a_mu(i,j-1) +a_Tmpv5
2486    a_mu(i,j) =a_mu(i,j) +a_Tmpv5
2487    a_Tmpv3 =rdnw(k)*a_Tmpv4
2488    a_dpn(i,k+1) =a_dpn(i,k+1) +a_Tmpv3
2489    a_dpn(i,k) =a_dpn(i,k) -a_Tmpv3
2490    a_Tmpv1 =(msfvy(i,j)/msfvx(i,j))*rdy*a_Tmpv2
2491    a_php(i,k,j) =a_php(i,k,j) +a_Tmpv1
2492    a_php(i,k,j-1) =a_php(i,k,j-1) -a_Tmpv1
2493    ENDDO
2494    ENDDO
2495    DO k =k_end, k_start+1, -1
2496    DO i =i_end, i_start, -1
2497    a_Tmpv6 =a_dpn(i,k)
2498    a_dpn(i,k) =0.0
2499    a_Tmpv5 =.5*a_Tmpv6
2500    a_Tmpv2 =a_Tmpv5
2501    a_Tmpv4 =a_Tmpv5
2502    a_Tmpv3 =fnp(k)*a_Tmpv4
2503    a_p(i,k-1,j) =a_p(i,k-1,j) +a_Tmpv3
2504    a_p(i,k-1,j-1) =a_p(i,k-1,j-1) +a_Tmpv3
2505    a_Tmpv1 =fnm(k)*a_Tmpv2
2506    a_p(i,k,j) =a_p(i,k,j) +a_Tmpv1
2507    a_p(i,k,j-1) =a_p(i,k,j-1) +a_Tmpv1
2508    ENDDO
2509    ENDDO
2510
2511    IF(top_lid) THEN
2512
2513    DO i =i_end, i_start, -1
2514    a_Tmpv9 =a_dpn(i,kde)
2515    a_dpn(i,kde) =0.0
2516    a_Tmpv8 =.5*a_Tmpv9
2517    a_Tmpv5 =a_Tmpv8
2518    a_Tmpv7 =a_Tmpv8
2519    a_Tmpv6 =cf3*a_Tmpv7
2520    a_p(i,kde-3,j-1) =a_p(i,kde-3,j-1) +a_Tmpv6
2521    a_p(i,kde-3,j) =a_p(i,kde-3,j) +a_Tmpv6
2522    a_Tmpv2 =a_Tmpv5
2523    a_Tmpv4 =a_Tmpv5
2524    a_Tmpv3 =cf2*a_Tmpv4
2525    a_p(i,kde-2,j-1) =a_p(i,kde-2,j-1) +a_Tmpv3
2526    a_p(i,kde-2,j) =a_p(i,kde-2,j) +a_Tmpv3
2527    a_Tmpv1 =cf1*a_Tmpv2
2528    a_p(i,kde-1,j-1) =a_p(i,kde-1,j-1) +a_Tmpv1
2529    a_p(i,kde-1,j) =a_p(i,kde-1,j) +a_Tmpv1
2530    ENDDO
2531
2532    ENDIF
2533    DO i =i_end, i_start, -1
2534    a_dpn(i,kde) =0.0
2535    a_Tmpv9 =a_dpn(i,1)
2536    a_dpn(i,1) =0.0
2537    a_Tmpv8 =.5*a_Tmpv9
2538    a_Tmpv5 =a_Tmpv8
2539    a_Tmpv7 =a_Tmpv8
2540    a_Tmpv6 =cf3*a_Tmpv7
2541    a_p(i,3,j) =a_p(i,3,j) +a_Tmpv6
2542    a_p(i,3,j-1) =a_p(i,3,j-1) +a_Tmpv6
2543    a_Tmpv2 =a_Tmpv5
2544    a_Tmpv4 =a_Tmpv5
2545    a_Tmpv3 =cf2*a_Tmpv4
2546    a_p(i,2,j) =a_p(i,2,j) +a_Tmpv3
2547    a_p(i,2,j-1) =a_p(i,2,j-1) +a_Tmpv3
2548    a_Tmpv1 =cf1*a_Tmpv2
2549    a_p(i,1,j) =a_p(i,1,j) +a_Tmpv1
2550    a_p(i,1,j-1) =a_p(i,1,j-1) +a_Tmpv1
2551    ENDDO
2552
2553    END IF
2554    DO k =k_end, k_start, -1
2555    DO i =i_end, i_start, -1
2556
2557    dpxy(i,k) =Tmpv303(i,k)
2558
2559    a_Tmpv11 =a_dpxy(i,k)
2560    a_dpxy(i,k) =0.0
2561    a_muv(i,j) =a_muv(i,j) +(msfvy(i,j)/msfvx(i,j))*.5*rdy*Tmpv302(i,k)*a_Tmpv11
2562    a_Tmpv10 =(msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*a_Tmpv11
2563    a_Tmpv7 =a_Tmpv10
2564    a_Tmpv9 =a_Tmpv10
2565    a_Tmpv8 =(pb(i,k,j)-pb(i,k,j-1))*a_Tmpv9
2566    a_al(i,k,j) =a_al(i,k,j) +a_Tmpv8
2567    a_al(i,k,j-1) =a_al(i,k,j-1) +a_Tmpv8
2568    a_Tmpv3 =a_Tmpv7
2569    a_Tmpv6 =a_Tmpv7
2570    a_Tmpv4 =Tmpv301(i,k)*a_Tmpv6
2571    a_Tmpv5 =Tmpv300(i,k)*a_Tmpv6
2572    a_p(i,k,j) =a_p(i,k,j) +a_Tmpv5
2573    a_p(i,k,j-1) =a_p(i,k,j-1) -a_Tmpv5
2574    a_alt(i,k,j) =a_alt(i,k,j) +a_Tmpv4
2575    a_alt(i,k,j-1) =a_alt(i,k,j-1) +a_Tmpv4
2576    a_Tmpv1 =a_Tmpv3
2577    a_Tmpv2 =a_Tmpv3
2578    a_ph(i,k,j) =a_ph(i,k,j) +a_Tmpv2
2579    a_ph(i,k,j-1) =a_ph(i,k,j-1) -a_Tmpv2
2580    a_ph(i,k+1,j) =a_ph(i,k+1,j) +a_Tmpv1
2581    a_ph(i,k+1,j-1) =a_ph(i,k+1,j-1) -a_Tmpv1
2582    ENDDO
2583    ENDDO
2584
2585    END IF
2586
2587    DO i =i_end, i_start, -1
2588    a_Tmpv3 =a_mudf_xy(i)
2589    a_mudf_xy(i) =0.0
2590    a_Tmpv2 =msfvx_inv(i,j)*a_Tmpv3
2591    a_Tmpv1 =-emdiv*dy*a_Tmpv2
2592    a_mudf(i,j) =a_mudf(i,j) +a_Tmpv1
2593    a_mudf(i,j-1) =a_mudf(i,j-1) -a_Tmpv1
2594    ENDDO
2595
2596    DO k =k_end, k_start, -1
2597    DO i =i_end, i_start, -1
2598    a_Tmpv1 =a_v(i,k,j)
2599    a_v(i,k,j) =0.0
2600    a_v(i,k,j) =a_v(i,k,j) +a_Tmpv1
2601    a_rv_tend(i,k,j) =a_rv_tend(i,k,j) +dts*a_Tmpv1
2602    ENDDO
2603    ENDDO
2604
2605    ENDDO
2606
2607 !LPB[19]
2608    DO j =j_end, j_start, -1
2609
2610    DO k =k_start, k_end
2611    DO i =i_start_u_tend, i_end_u_tend
2612    Tmpv001 =u(i,k,j) +dts*ru_tend(i,k,j)
2613 !  u(i,k,j) =Tmpv001
2614
2615    ENDDO
2616    ENDDO
2617    DO i =i_start_up, i_end_up
2618    Tmpv001 =mudf(i,j) -mudf(i-1,j)
2619    Tmpv002 =-emdiv*dx*Tmpv001
2620    Tmpv003 =Tmpv002/msfuy(i,j)
2621 !  mudf_xy(i) =Tmpv003
2622
2623    ENDDO
2624
2625    DO k =k_start, k_end
2626    DO i =i_start_up, i_end_up
2627    Tmpv001 =ph(i,k+1,j) -ph(i-1,k+1,j)
2628    Tmpv002 =ph(i,k,j) -ph(i-1,k,j)
2629    Tmpv003 =Tmpv001 +Tmpv002
2630    Tmpv004 =alt(i,k,j) +alt(i-1,k,j)
2631    Tmpv005 =p(i,k,j) -p(i-1,k,j)
2632    Tmpv300(i,k) =Tmpv004
2633    Tmpv301(i,k) =Tmpv005
2634    Tmpv006 =Tmpv300(i,k)*Tmpv301(i,k)
2635    Tmpv007 =Tmpv003 +Tmpv006
2636    Tmpv008 =al(i,k,j) +al(i-1,k,j)
2637    Tmpv009 =Tmpv008*(pb(i,k,j)-pb(i-1,k,j))
2638    Tmpv010 =Tmpv007 +Tmpv009
2639    Tmpv302(i,k) =Tmpv010
2640    Tmpv011 =(msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*Tmpv302(i,k)
2641    Tmpv303(i,k) =dpxy(i,k)
2642    dpxy(i,k) =Tmpv011
2643
2644    ENDDO
2645    ENDDO
2646    IF(non_hydrostatic) THEN
2647    DO i =i_start_up, i_end_up
2648    Tmpv001 =p(i,1,j) +p(i-1,1,j)
2649    Tmpv002 =cf1*Tmpv001
2650    Tmpv003 =p(i,2,j) +p(i-1,2,j)
2651    Tmpv004 =cf2*Tmpv003
2652    Tmpv005 =Tmpv002 +Tmpv004
2653    Tmpv006 =p(i,3,j) +p(i-1,3,j)
2654    Tmpv007 =cf3*Tmpv006
2655    Tmpv008 =Tmpv005 +Tmpv007
2656    Tmpv009 =.5*Tmpv008
2657    dpn(i,1) =Tmpv009
2658
2659    dpn(i,kde) =0.
2660
2661    ENDDO
2662
2663    IF(top_lid) THEN
2664    DO i =i_start_up, i_end_up
2665    Tmpv001 =p(i-1,kde-1,j) +p(i,kde-1,j)
2666    Tmpv002 =cf1*Tmpv001
2667    Tmpv003 =p(i-1,kde-2,j) +p(i,kde-2,j)
2668    Tmpv004 =cf2*Tmpv003
2669    Tmpv005 =Tmpv002 +Tmpv004
2670    Tmpv006 =p(i-1,kde-3,j) +p(i,kde-3,j)
2671    Tmpv007 =cf3*Tmpv006
2672    Tmpv008 =Tmpv005 +Tmpv007
2673    Tmpv009 =.5*Tmpv008
2674    dpn(i,kde) =Tmpv009
2675
2676    ENDDO
2677
2678    ENDIF
2679    DO k =k_start+1, k_end
2680    DO i =i_start_up, i_end_up
2681    Tmpv001 =p(i,k,j) +p(i-1,k,j)
2682    Tmpv002 =fnm(k)*Tmpv001
2683    Tmpv003 =p(i,k-1,j) +p(i-1,k-1,j)
2684    Tmpv004 =fnp(k)*Tmpv003
2685    Tmpv005 =Tmpv002 +Tmpv004
2686    Tmpv006 =.5*Tmpv005
2687    dpn(i,k) =Tmpv006
2688
2689    ENDDO
2690    ENDDO
2691    DO k =k_start, k_end
2692    DO i =i_start_up, i_end_up
2693    Tmpv001 =php(i,k,j) -php(i-1,k,j)
2694    Tmpv002 =(msfux(i,j)/msfuy(i,j))*rdx*Tmpv001
2695    Tmpv003 =dpn(i,k+1) -dpn(i,k)
2696    Tmpv004 =rdnw(k)*Tmpv003
2697    Tmpv005 =mu(i-1,j) +mu(i,j)
2698    Tmpv006 =.5*Tmpv005
2699    Tmpv007 =Tmpv004 -Tmpv006
2700    Tmpv304(i,k) =Tmpv002
2701    Tmpv305(i,k) =Tmpv007
2702    Tmpv008 =Tmpv304(i,k)*Tmpv305(i,k)
2703    Tmpv009 =dpxy(i,k) +Tmpv008
2704    Tmpv306(i,k) =dpxy(i,k)
2705    dpxy(i,k) =Tmpv009
2706
2707    ENDDO
2708    ENDDO
2709    END IF
2710    DO k =k_start, k_end
2711    DO i =i_start_up, i_end_up
2712    Tmpv001 =dts*cqu(i,k,j)*dpxy(i,k)
2713    Tmpv002 =u(i,k,j) -Tmpv001
2714    Tmpv003 =Tmpv002 +mudf_xy(i)
2715 !  u(i,k,j) =Tmpv003
2716
2717    ENDDO
2718    ENDDO
2719
2720    DO k =k_end, k_start, -1
2721    DO i =i_end_up, i_start_up, -1
2722    a_Tmpv3 =a_u(i,k,j)
2723    a_u(i,k,j) =0.0
2724    a_Tmpv2 =a_Tmpv3
2725    a_mudf_xy(i) =a_mudf_xy(i) +a_Tmpv3
2726    a_u(i,k,j) =a_u(i,k,j) +a_Tmpv2
2727    a_Tmpv1 =-a_Tmpv2
2728    a_cqu(i,k,j) =a_cqu(i,k,j) +dts*dpxy(i,k)*a_Tmpv1
2729    a_dpxy(i,k) =a_dpxy(i,k) +dts*cqu(i,k,j)*a_Tmpv1
2730    ENDDO
2731    ENDDO
2732
2733    IF(non_hydrostatic) THEN
2734
2735    DO k =k_end, k_start, -1
2736    DO i =i_end_up, i_start_up, -1
2737
2738    dpxy(i,k) =Tmpv306(i,k)
2739
2740    a_Tmpv9 =a_dpxy(i,k)
2741    a_dpxy(i,k) =0.0
2742    a_dpxy(i,k) =a_dpxy(i,k) +a_Tmpv9
2743    a_Tmpv8 =a_Tmpv9
2744    a_Tmpv2 =Tmpv305(i,k)*a_Tmpv8
2745    a_Tmpv7 =Tmpv304(i,k)*a_Tmpv8
2746    a_Tmpv4 =a_Tmpv7
2747    a_Tmpv6 =-a_Tmpv7
2748    a_Tmpv5 =.5*a_Tmpv6
2749    a_mu(i-1,j) =a_mu(i-1,j) +a_Tmpv5
2750    a_mu(i,j) =a_mu(i,j) +a_Tmpv5
2751    a_Tmpv3 =rdnw(k)*a_Tmpv4
2752    a_dpn(i,k+1) =a_dpn(i,k+1) +a_Tmpv3
2753    a_dpn(i,k) =a_dpn(i,k) -a_Tmpv3
2754    a_Tmpv1 =(msfux(i,j)/msfuy(i,j))*rdx*a_Tmpv2
2755    a_php(i,k,j) =a_php(i,k,j) +a_Tmpv1
2756    a_php(i-1,k,j) =a_php(i-1,k,j) -a_Tmpv1
2757    ENDDO
2758    ENDDO
2759    DO k =k_end, k_start+1, -1
2760    DO i =i_end_up, i_start_up, -1
2761    a_Tmpv6 =a_dpn(i,k)
2762    a_dpn(i,k) =0.0
2763    a_Tmpv5 =.5*a_Tmpv6
2764    a_Tmpv2 =a_Tmpv5
2765    a_Tmpv4 =a_Tmpv5
2766    a_Tmpv3 =fnp(k)*a_Tmpv4
2767    a_p(i,k-1,j) =a_p(i,k-1,j) +a_Tmpv3
2768    a_p(i-1,k-1,j) =a_p(i-1,k-1,j) +a_Tmpv3
2769    a_Tmpv1 =fnm(k)*a_Tmpv2
2770    a_p(i,k,j) =a_p(i,k,j) +a_Tmpv1
2771    a_p(i-1,k,j) =a_p(i-1,k,j) +a_Tmpv1
2772    ENDDO
2773    ENDDO
2774
2775    IF(top_lid) THEN
2776
2777    DO i =i_end_up, i_start_up, -1
2778    a_Tmpv9 =a_dpn(i,kde)
2779    a_dpn(i,kde) =0.0
2780    a_Tmpv8 =.5*a_Tmpv9
2781    a_Tmpv5 =a_Tmpv8
2782    a_Tmpv7 =a_Tmpv8
2783    a_Tmpv6 =cf3*a_Tmpv7
2784    a_p(i-1,kde-3,j) =a_p(i-1,kde-3,j) +a_Tmpv6
2785    a_p(i,kde-3,j) =a_p(i,kde-3,j) +a_Tmpv6
2786    a_Tmpv2 =a_Tmpv5
2787    a_Tmpv4 =a_Tmpv5
2788    a_Tmpv3 =cf2*a_Tmpv4
2789    a_p(i-1,kde-2,j) =a_p(i-1,kde-2,j) +a_Tmpv3
2790    a_p(i,kde-2,j) =a_p(i,kde-2,j) +a_Tmpv3
2791    a_Tmpv1 =cf1*a_Tmpv2
2792    a_p(i-1,kde-1,j) =a_p(i-1,kde-1,j) +a_Tmpv1
2793    a_p(i,kde-1,j) =a_p(i,kde-1,j) +a_Tmpv1
2794    ENDDO
2795
2796    ENDIF
2797    DO i =i_end_up, i_start_up, -1
2798    a_dpn(i,kde) =0.0
2799    a_Tmpv9 =a_dpn(i,1)
2800    a_dpn(i,1) =0.0
2801    a_Tmpv8 =.5*a_Tmpv9
2802    a_Tmpv5 =a_Tmpv8
2803    a_Tmpv7 =a_Tmpv8
2804    a_Tmpv6 =cf3*a_Tmpv7
2805    a_p(i,3,j) =a_p(i,3,j) +a_Tmpv6
2806    a_p(i-1,3,j) =a_p(i-1,3,j) +a_Tmpv6
2807    a_Tmpv2 =a_Tmpv5
2808    a_Tmpv4 =a_Tmpv5
2809    a_Tmpv3 =cf2*a_Tmpv4
2810    a_p(i,2,j) =a_p(i,2,j) +a_Tmpv3
2811    a_p(i-1,2,j) =a_p(i-1,2,j) +a_Tmpv3
2812    a_Tmpv1 =cf1*a_Tmpv2
2813    a_p(i,1,j) =a_p(i,1,j) +a_Tmpv1
2814    a_p(i-1,1,j) =a_p(i-1,1,j) +a_Tmpv1
2815    ENDDO
2816
2817    END IF
2818
2819    DO k =k_end, k_start, -1
2820    DO i =i_end_up, i_start_up, -1
2821
2822    dpxy(i,k) =Tmpv303(i,k)
2823
2824    a_Tmpv11 =a_dpxy(i,k)
2825    a_dpxy(i,k) =0.0
2826    a_muu(i,j) =a_muu(i,j) +(msfux(i,j)/msfuy(i,j))*.5*rdx*Tmpv302(i,k)*a_Tmpv11
2827    a_Tmpv10 =(msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*a_Tmpv11
2828    a_Tmpv7 =a_Tmpv10
2829    a_Tmpv9 =a_Tmpv10
2830    a_Tmpv8 =(pb(i,k,j)-pb(i-1,k,j))*a_Tmpv9
2831    a_al(i,k,j) =a_al(i,k,j) +a_Tmpv8
2832    a_al(i-1,k,j) =a_al(i-1,k,j) +a_Tmpv8
2833    a_Tmpv3 =a_Tmpv7
2834    a_Tmpv6 =a_Tmpv7
2835    a_Tmpv4 =Tmpv301(i,k)*a_Tmpv6
2836    a_Tmpv5 =Tmpv300(i,k)*a_Tmpv6
2837    a_p(i,k,j) =a_p(i,k,j) +a_Tmpv5
2838    a_p(i-1,k,j) =a_p(i-1,k,j) -a_Tmpv5
2839    a_alt(i,k,j) =a_alt(i,k,j) +a_Tmpv4
2840    a_alt(i-1,k,j) =a_alt(i-1,k,j) +a_Tmpv4
2841    a_Tmpv1 =a_Tmpv3
2842    a_Tmpv2 =a_Tmpv3
2843    a_ph(i,k,j) =a_ph(i,k,j) +a_Tmpv2
2844    a_ph(i-1,k,j) =a_ph(i-1,k,j) -a_Tmpv2
2845    a_ph(i,k+1,j) =a_ph(i,k+1,j) +a_Tmpv1
2846    a_ph(i-1,k+1,j) =a_ph(i-1,k+1,j) -a_Tmpv1
2847    ENDDO
2848    ENDDO
2849
2850    DO i =i_end_up, i_start_up, -1
2851    a_Tmpv3 =a_mudf_xy(i)
2852    a_mudf_xy(i) =0.0
2853    a_Tmpv2 =a_Tmpv3/msfuy(i,j)
2854    a_Tmpv1 =-emdiv*dx*a_Tmpv2
2855    a_mudf(i,j) =a_mudf(i,j) +a_Tmpv1
2856    a_mudf(i-1,j) =a_mudf(i-1,j) -a_Tmpv1
2857    ENDDO
2858
2859    DO k =k_end, k_start, -1
2860    DO i =i_end_u_tend, i_start_u_tend, -1
2861    a_Tmpv1 =a_u(i,k,j)
2862    a_u(i,k,j) =0.0
2863    a_u(i,k,j) =a_u(i,k,j) +a_Tmpv1
2864    a_ru_tend(i,k,j) =a_ru_tend(i,k,j) +dts*a_Tmpv1
2865    ENDDO
2866    ENDDO
2867
2868    ENDDO
2869
2870 !LPB[18]
2871 !  dx =1./rdx
2872 !  dy =1./rdy
2873
2874 !LPB[17]
2875
2876 !  IF( config_flags%symmetric_ye .and. (jte == jde) ) THEN
2877 !  j_end_v_tend =j_end_v_tend-1
2878 !  END IF
2879
2880 !  IF( config_flags%symmetric_ye .and. (jte == jde) ) THEN
2881
2882 !  END IF
2883
2884 !LPB[16]
2885
2886 !LPB[15]
2887
2888 !  IF( config_flags%symmetric_ys .and. (jts == jds) ) THEN
2889 !  j_start_v_tend =j_start_v_tend+1
2890 !  END IF
2891
2892 !  IF( config_flags%symmetric_ys .and. (jts == jds) ) THEN
2893
2894 !  END IF
2895
2896 !LPB[14]
2897
2898 !LPB[13]
2899
2900 !  IF( config_flags%symmetric_xe .and. (ite == ide) ) THEN
2901 !  i_end_u_tend =i_end_u_tend-1
2902 !  END IF
2903
2904 !  IF( config_flags%symmetric_xe .and. (ite == ide) ) THEN
2905
2906 !  END IF
2907
2908 !LPB[12]
2909
2910 !LPB[11]
2911
2912 !  IF( config_flags%symmetric_xs .and. (its == ids) ) THEN
2913 !  i_start_u_tend =i_start_u_tend+1
2914 !  END IF
2915
2916 !  IF( config_flags%symmetric_xs .and. (its == ids) ) THEN
2917
2918 !  END IF
2919
2920 !LPB[10]
2921 !  i_start_u_tend =i_start
2922 !  i_end_u_tend =i_endu
2923 !  j_start_v_tend =j_start
2924 !  j_end_v_tend =j_endv
2925
2926 !LPB[9]
2927
2928 !  IF( (config_flags%open_ye     .or.               config_flags%symmetric_ye  .or.                 config_flags%polar   )                  .and. (jte == jde) ) THEN
2929 !  j_end_vp =j_end_vp-1
2930 !  END IF
2931
2932 !  IF( (config_flags%open_ye     .or.   &
2933 !               config_flags%symmetric_ye  .or.     &
2934 !               config_flags%polar   )    &
2935 !               .and. (jte == jde) ) THEN
2936
2937 !  END IF
2938
2939 !LPB[8]
2940
2941 !LPB[7]
2942
2943 !  IF( (config_flags%open_ys    .or.                        config_flags%symmetric_ys  .or.                 config_flags%polar   )                           .and. (jts == jds) ) THEN
2944 !  j_start_vp =j_start_vp+1
2945 !  END IF
2946
2947 !  IF( (config_flags%open_ys    .or.     &
2948 !               config_flags%symmetric_ys  .or.     &
2949 !               config_flags%polar   )    &
2950 !                        .and. (jts == jds) ) THEN
2951
2952 !  END IF
2953
2954 !LPB[6]
2955
2956 !LPB[5]
2957
2958 !  IF( (config_flags%open_xe    .or.                config_flags%symmetric_xe   )                    .and. (ite == ide) ) THEN
2959 !  i_end_up =i_end_up-1
2960 !  END IF
2961
2962 !  IF( (config_flags%open_xe    .or.    &
2963 !               config_flags%symmetric_xe   )   &
2964 !                .and. (ite == ide) ) THEN
2965
2966 !  END IF
2967
2968 !LPB[4]
2969
2970 !LPB[3]
2971
2972 !  IF( (config_flags%open_xs   .or.                         config_flags%symmetric_xs   )                   .and. (its == ids) ) THEN
2973 !  i_start_up =i_start_up+1
2974 !  END IF
2975
2976 !  IF( (config_flags%open_xs   .or.       &
2977 !               config_flags%symmetric_xs   )     &
2978 !               .and. (its == ids) ) THEN
2979
2980 !  END IF
2981
2982 !LPB[2]
2983 !  i_start_up =i_start
2984 !  i_end_up =i_endu
2985 !  j_start_up =j_start
2986 !  j_end_up =j_end
2987 !  i_start_vp =i_start
2988 !  i_end_vp =i_end
2989 !  j_start_vp =j_start
2990 !  j_end_vp =j_endv
2991
2992 !LPB[1]
2993
2994 !  IF( config_flags%nested .or. config_flags%specified ) THEN
2995 !  i_start =max(its, ids+spec_zone)
2996 !  i_end =min(ite, ide-spec_zone-1)
2997 !  j_start =max(jts, jds+spec_zone)
2998 !  j_end =min(jte, jde-spec_zone-1)
2999 !  k_start =kts
3000 !  k_end =min(kte, kde-1)
3001 !  i_endu =min(ite, ide-spec_zone)
3002 !  j_endv =min(jte, jde-spec_zone)
3003 !  k_endw =k_end
3004 !  IF( config_flags%periodic_x) THEN
3005 !  i_start =its
3006 !  i_end =min(ite, ide-1)
3007 !  i_endu =ite
3008 !  ENDIF
3009 !  ELSE
3010 !  i_start =its
3011 !  i_end =min(ite, ide-1)
3012 !  j_start =jts
3013 !  j_end =min(jte, jde-1)
3014 !  k_start =kts
3015 !  k_end =kte-1
3016 !  i_endu =ite
3017 !  j_endv =jte
3018 !  k_endw =k_end
3019 !  ENDIF
3020
3021    IF( config_flags%nested .or. config_flags%specified ) THEN
3022
3023    IF( config_flags%periodic_x) THEN
3024
3025    ENDIF
3026
3027    ELSE
3028
3029    ENDIF
3030
3031 !LPB[0]
3032
3033    END SUBROUTINE a_advance_uv
3034
3035    SUBROUTINE a_advance_mu_t(ww,a_ww,ww_1,a_ww_1,u,a_u,u_1,a_u_1,v, &
3036    a_v,v_1,a_v_1,mu,a_mu,mut,a_mut,muave,a_muave,muts,a_muts,muu,a_muu, &
3037    muv,a_muv,mudf,a_mudf,uam,a_uam,vam,a_vam,wwam,a_wwam,t,a_t, &
3038    t_1,a_t_1,t_ave,a_t_ave,ft,a_ft,mu_tend,a_mu_tend,rdx,rdy,dts,epssm,dnw, &
3039    fnm,fnp,rdnw,msfux,msfuy,msfvx,msfvx_inv,msfvy,msftx,msfty,step,config_flags,ids,ide, &
3040    jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte)
3041
3042 !PART I: DECLARATION OF VARIABLES
3043
3044    IMPLICIT NONE
3045
3046    INTEGER :: K0_ADJ,K1_ADJ,K2_ADJ,K3_ADJ
3047    TYPE(grid_config_rec_type) :: config_flags
3048    INTEGER :: ids,ide,jds,jde,kds,kde
3049    INTEGER :: ims,ime,jms,jme,kms,kme
3050    INTEGER :: its,ite,jts,jte,kts,kte
3051    INTEGER :: step
3052    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_u,u,a_v,v,a_u_1,u_1,a_v_1,v_1, &
3053    a_t_1,t_1,a_ft,ft
3054    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_ww,ww,a_ww_1,ww_1,a_t,t,a_t_ave, &
3055    t_ave,a_uam,uam,a_vam,vam,a_wwam,wwam
3056    REAL,DIMENSION(ims:ime,jms:jme) :: a_muu,muu,a_muv,muv,a_mut,mut,msfux,msfuy, &
3057    msfvx,msfvx_inv,msfvy,msftx,msfty,a_mu_tend,mu_tend
3058    REAL,DIMENSION(ims:ime,jms:jme) :: a_muave,muave,a_muts,muts,a_mudf,mudf
3059    REAL,DIMENSION(ims:ime,jms:jme) :: a_mu,mu
3060    REAL,DIMENSION(kms:kme) :: fnm,fnp,dnw,rdnw
3061    REAL :: rdx,rdy,dts,epssm
3062    REAL,DIMENSION(its:ite,kts:kte) :: a_wdtn,wdtn,a_dvdxi,dvdxi
3063    REAL,DIMENSION(its:ite) :: a_dmdt,dmdt
3064    INTEGER :: i,j,k,i_start,i_end,j_start,j_end,k_start,k_end
3065    INTEGER :: i_endu,j_endv
3066    REAL :: a_acc,acc
3067
3068    REAL :: a_Tmpv1,Tmpv001,a_Tmpv2,Tmpv002,a_Tmpv3,Tmpv003,a_Tmpv4,Tmpv004, &
3069    a_Tmpv5,Tmpv005,a_Tmpv6,Tmpv006,a_Tmpv7,Tmpv007,a_Tmpv8,Tmpv008,a_Tmpv9, &
3070    Tmpv009,a_Tmpv10,Tmpv010,a_Tmpv11,Tmpv011,a_Tmpv12,Tmpv012,a_Tmpv13,Tmpv013, &
3071    a_Tmpv14,Tmpv014,a_Tmpv15,Tmpv015,a_Tmpv16,Tmpv016,a_Tmpv17,Tmpv017, &
3072    a_Tmpv18,Tmpv018,a_Tmpv19,Tmpv019
3073    REAL,DIMENSION(its:min(ite,ide-1),2:kte-1) :: Tmpv300
3074    REAL,DIMENSION(its:min(ite,ide-1),kte-1) :: Tmpv301
3075    REAL,DIMENSION(its:min(ite,ide-1),kte-1) :: Tmpv302
3076    REAL,DIMENSION(its:min(ite,ide-1),kte-1) :: Tmpv303
3077    REAL,DIMENSION(its:min(ite,ide-1),kte-1) :: Tmpv304
3078
3079 !PART II: CALCULATIONS OF B. S. TRAJECTORY
3080
3081 !LPB[0]
3082
3083      i_start = its
3084      i_end   = min(ite,ide-1)
3085      j_start = jts
3086      j_end   = min(jte,jde-1)
3087      k_start = kts
3088      k_end   = kte-1
3089
3090 !LPB[1]
3091   IF ( .NOT. config_flags%periodic_x )THEN
3092
3093      IF ( config_flags%specified .or. config_flags%nested ) then
3094
3095           i_start = max(its,ids+1)
3096           i_end   = min(ite,ide-2)
3097         ENDIF
3098
3099    ENDIF
3100
3101 !LPB[2]
3102
3103 !LPB[3]
3104
3105   IF ( config_flags%specified .or. config_flags%nested ) then
3106
3107         j_start = max(jts,jds+1)
3108         j_end   = min(jte,jde-2)
3109
3110    ENDIF
3111
3112 !LPB[4]
3113
3114      i_endu = ite
3115      j_endv = jte
3116
3117 !LPB[5]
3118       DO j = j_start, j_end
3119
3120         DO i=i_start, i_end
3121                dmdt(i) = 0.
3122         ENDDO
3123
3124         DO k=k_start, k_end
3125         DO i=i_start, i_end
3126             dvdxi(i,k) = msftx(i,j)*msfty(i,j)*(                                        &
3127                         rdy*( (v(i,k,j+1)+muv(i,j+1)*v_1(i,k,j+1)*msfvx_inv(i,j+1))     &
3128                              -(v(i,k,j  )+muv(i,j  )*v_1(i,k,j  )*msfvx_inv(i,j  )) )   &
3129                        +rdx*( (u(i+1,k,j)+muu(i+1,j)*u_1(i+1,k,j)/msfuy(i+1,j))         &
3130                              -(u(i,k,j  )+muu(i  ,j)*u_1(i,k,j  )/msfuy(i  ,j)) ))
3131            dmdt(i)    = dmdt(i) + dnw(k)*dvdxi(i,k)
3132         ENDDO
3133         ENDDO
3134
3135         DO i=i_start, i_end
3136           muave(i,j) = mu(i,j)
3137           mu(i,j) = mu(i,j)+dts*(dmdt(i)+mu_tend(i,j))
3138           mudf(i,j) = (dmdt(i)+mu_tend(i,j))
3139           muts(i,j) = mut(i,j)+mu(i,j)
3140           muave(i,j) =.5*((1.+epssm)*mu(i,j)+(1.-epssm)*muave(i,j))
3141         ENDDO
3142
3143         DO k=2,k_end
3144         DO i=i_start, i_end
3145           ww(i,k,j)=ww(i,k-1,j)-dnw(k-1)*(dmdt(i)+dvdxi(i,k-1)+mu_tend(i,j))/msfty(i,j)
3146         ENDDO
3147         END DO
3148
3149         DO k=1,k_end
3150         DO i=i_start, i_end
3151           ww(i,k,j)=ww(i,k,j)-ww_1(i,k,j)
3152         END DO
3153         END DO
3154
3155       ENDDO
3156
3157 ! Remarked by Ning Pan, 2010-08-31 : LPB[6]
3158 !LPB[6]
3159 !      DO j=j_start, j_end
3160
3161 !        DO k=1,k_end
3162 !        DO i=i_start, i_end
3163 !          t_ave(i,k,j) = t(i,k,j)
3164 !          t   (i,k,j) = t(i,k,j) + msfty(i,j)*dts*ft(i,k,j)
3165 !        END DO
3166 !        END DO
3167
3168 !      ENDDO
3169
3170 !!LPB[7]
3171 !      DO j=j_start, j_end
3172
3173 !        DO i=i_start, i_end
3174 !          wdtn(i,1  )=0.
3175 !          wdtn(i,kde)=0.
3176 !        ENDDO
3177
3178 !        DO k=2,k_end
3179 !        DO i=i_start, i_end
3180 !           wdtn(i,k)= ww(i,k,j)*(fnm(k)*t_1(i,k  ,j)+fnp(k)*t_1(i,k-1,j))
3181 !        ENDDO
3182 !        ENDDO
3183
3184 !        DO k=1,k_end
3185 !        DO i=i_start, i_end
3186 !          t(i,k,j) = t(i,k,j) - dts*msfty(i,j)*(                &
3187 !                              &
3188 !                             msftx(i,j)*(                       &
3189 !                  .5*rdy*                                       &
3190 !                 ( v(i,k,j+1)*(t_1(i,k,j+1)+t_1(i,k, j ))       &
3191 !                  -v(i,k,j  )*(t_1(i,k, j )+t_1(i,k,j-1)) )     &
3192 !                + .5*rdx*                                       &
3193 !                 ( u(i+1,k,j)*(t_1(i+1,k,j)+t_1(i  ,k,j))       &
3194 !                  -u(i  ,k,j)*(t_1(i  ,k,j)+t_1(i-1,k,j)) ) )   &
3195 !                + rdnw(k)*( wdtn(i,k+1)-wdtn(i,k) ) )
3196 !        ENDDO
3197 !        ENDDO
3198
3199 !      ENDDO
3200
3201 !PART III: INITIALIZATION OF LOCAL ADJOINT PERTURBATIONS
3202
3203    Do K1_ADJ =kts, kte
3204    Do K0_ADJ =its, ite
3205    a_wdtn(K0_ADJ,K1_ADJ) =0.0
3206    End Do
3207    End Do
3208
3209    Do K1_ADJ =kts, kte
3210    Do K0_ADJ =its, ite
3211    a_dvdxi(K0_ADJ,K1_ADJ) =0.0
3212    End Do
3213    End Do
3214
3215    Do K0_ADJ =its, ite
3216    a_dmdt(K0_ADJ) =0.0
3217    End Do
3218
3219    a_acc =0.0
3220
3221 !PART IV: REVERSE/BACKWARD ACCUMULATIONS
3222
3223 !LPB[7]
3224    DO j =j_end, j_start, -1
3225
3226 ! Remarked by Ning Pan, 2010-08-31
3227 !   DO i =i_start, i_end
3228 !   wdtn(i,1) =0.
3229
3230 !   wdtn(i,kde) =0.
3231
3232 !   ENDDO
3233
3234    DO k =2, k_end
3235    DO i =i_start, i_end
3236    Tmpv001 =fnm(k)*t_1(i,k,j) +fnp(k)*t_1(i,k-1,j)
3237    Tmpv300(i,k) =Tmpv001
3238 ! Remarked by Ning Pan, 2010-08-31
3239 !   Tmpv002 =ww(i,k,j)*Tmpv300(i,k)
3240 !!  wdtn(i,k) =Tmpv002
3241
3242    ENDDO
3243    ENDDO
3244    DO k =1, k_end
3245    DO i =i_start, i_end
3246    Tmpv001 =t_1(i,k,j+1) +t_1(i,k,j)
3247    Tmpv301(i,k) =Tmpv001
3248    Tmpv002 =v(i,k,j+1)*Tmpv301(i,k)
3249    Tmpv003 =t_1(i,k,j) +t_1(i,k,j-1)
3250    Tmpv302(i,k) =Tmpv003
3251    Tmpv004 =v(i,k,j)*Tmpv302(i,k)
3252    Tmpv005 =Tmpv002 -Tmpv004
3253    Tmpv006 =.5*rdy*Tmpv005
3254    Tmpv007 =t_1(i+1,k,j) +t_1(i,k,j)
3255    Tmpv303(i,k) =Tmpv007
3256    Tmpv008 =u(i+1,k,j)*Tmpv303(i,k)
3257    Tmpv009 =t_1(i,k,j) +t_1(i-1,k,j)
3258    Tmpv304(i,k) =Tmpv009
3259 ! Remarked by Ning Pan, 2010-08-31
3260 !   Tmpv010 =u(i,k,j)*Tmpv304(i,k)
3261 !   Tmpv011 =Tmpv008 -Tmpv010
3262 !   Tmpv012 =.5*rdx*Tmpv011
3263 !   Tmpv013 =Tmpv006 +Tmpv012
3264 !   Tmpv014 =msftx(i,j)*Tmpv013
3265 !   Tmpv015 =wdtn(i,k+1) -wdtn(i,k)
3266 !   Tmpv016 =rdnw(k)*Tmpv015
3267 !   Tmpv017 =Tmpv014 +Tmpv016
3268 !   Tmpv018 =dts*msfty(i,j)*Tmpv017
3269 !   Tmpv019 =t(i,k,j) -Tmpv018
3270 !!  t(i,k,j) =Tmpv019
3271
3272    ENDDO
3273    ENDDO
3274
3275    DO k =k_end, 1, -1
3276    DO i =i_end, i_start, -1
3277    a_Tmpv19 =a_t(i,k,j)
3278    a_t(i,k,j) =0.0
3279    a_t(i,k,j) =a_t(i,k,j) +a_Tmpv19
3280    a_Tmpv18 =-a_Tmpv19
3281    a_Tmpv17 =dts*msfty(i,j)*a_Tmpv18
3282    a_Tmpv14 =a_Tmpv17
3283    a_Tmpv16 =a_Tmpv17
3284    a_Tmpv15 =rdnw(k)*a_Tmpv16
3285    a_wdtn(i,k+1) =a_wdtn(i,k+1) +a_Tmpv15
3286    a_wdtn(i,k) =a_wdtn(i,k) -a_Tmpv15
3287    a_Tmpv13 =msftx(i,j)*a_Tmpv14
3288    a_Tmpv6 =a_Tmpv13
3289    a_Tmpv12 =a_Tmpv13
3290    a_Tmpv11 =.5*rdx*a_Tmpv12
3291    a_Tmpv8 =a_Tmpv11
3292    a_Tmpv10 =-a_Tmpv11
3293    a_u(i,k,j) =a_u(i,k,j) +Tmpv304(i,k)*a_Tmpv10
3294    a_Tmpv9 =u(i,k,j)*a_Tmpv10
3295    a_t_1(i,k,j) =a_t_1(i,k,j) +a_Tmpv9
3296    a_t_1(i-1,k,j) =a_t_1(i-1,k,j) +a_Tmpv9
3297    a_u(i+1,k,j) =a_u(i+1,k,j) +Tmpv303(i,k)*a_Tmpv8
3298    a_Tmpv7 =u(i+1,k,j)*a_Tmpv8
3299    a_t_1(i+1,k,j) =a_t_1(i+1,k,j) +a_Tmpv7
3300    a_t_1(i,k,j) =a_t_1(i,k,j) +a_Tmpv7
3301    a_Tmpv5 =.5*rdy*a_Tmpv6
3302    a_Tmpv2 =a_Tmpv5
3303    a_Tmpv4 =-a_Tmpv5
3304    a_v(i,k,j) =a_v(i,k,j) +Tmpv302(i,k)*a_Tmpv4
3305    a_Tmpv3 =v(i,k,j)*a_Tmpv4
3306    a_t_1(i,k,j) =a_t_1(i,k,j) +a_Tmpv3
3307    a_t_1(i,k,j-1) =a_t_1(i,k,j-1) +a_Tmpv3
3308    a_v(i,k,j+1) =a_v(i,k,j+1) +Tmpv301(i,k)*a_Tmpv2
3309    a_Tmpv1 =v(i,k,j+1)*a_Tmpv2
3310    a_t_1(i,k,j+1) =a_t_1(i,k,j+1) +a_Tmpv1
3311    a_t_1(i,k,j) =a_t_1(i,k,j) +a_Tmpv1
3312    ENDDO
3313    ENDDO
3314
3315    DO k =k_end, 2, -1
3316    DO i =i_end, i_start, -1
3317    a_Tmpv2 =a_wdtn(i,k)
3318    a_wdtn(i,k) =0.0
3319    a_ww(i,k,j) =a_ww(i,k,j) +Tmpv300(i,k)*a_Tmpv2
3320    a_Tmpv1 =ww(i,k,j)*a_Tmpv2
3321    a_t_1(i,k,j) =a_t_1(i,k,j) +fnm(k)*a_Tmpv1
3322    a_t_1(i,k-1,j) =a_t_1(i,k-1,j) +fnp(k)*a_Tmpv1
3323    ENDDO
3324    ENDDO
3325
3326    DO i =i_end, i_start, -1
3327    a_wdtn(i,kde) =0.0
3328    a_wdtn(i,1) =0.0
3329    ENDDO
3330
3331    ENDDO
3332
3333 !LPB[6]
3334    DO j =j_end, j_start, -1
3335
3336 !  DO k =1, k_end
3337 !  DO i =i_start, i_end
3338 !  t_ave(i,k,j) =t(i,k,j)
3339
3340 !  Tmpv001 =t(i,k,j) +msfty(i,j)*dts*ft(i,k,j)
3341 !  t(i,k,j) =Tmpv001
3342
3343 !  ENDDO
3344 !  ENDDO
3345
3346    DO k =k_end, 1, -1
3347    DO i =i_end, i_start, -1
3348    a_Tmpv1 =a_t(i,k,j)
3349    a_t(i,k,j) =0.0
3350    a_t(i,k,j) =a_t(i,k,j) +a_Tmpv1
3351    a_ft(i,k,j) =a_ft(i,k,j) +msfty(i,j)*dts*a_Tmpv1
3352    a_t(i,k,j) =a_t(i,k,j) +a_t_ave(i,k,j)
3353    a_t_ave(i,k,j) =0.0
3354    ENDDO
3355    ENDDO
3356
3357    ENDDO
3358
3359 !LPB[5]
3360    DO j =j_end, j_start, -1
3361
3362 ! Remarked by Ning Pan, 2010-08-31 : not need to recalculate
3363 !   DO i =i_start, i_end
3364 !   dmdt(i) =0.
3365
3366 !   ENDDO
3367
3368 !   DO k =k_start, k_end
3369 !   DO i =i_start, i_end
3370 !   Tmpv001 =muv(i,j+1)*v_1(i,k,j+1)
3371 !   Tmpv002 =Tmpv001*msfvx_inv(i,j+1)
3372 !   Tmpv003 =v(i,k,j+1) +Tmpv002
3373 !   Tmpv004 =muv(i,j)*v_1(i,k,j)
3374 !   Tmpv005 =Tmpv004*msfvx_inv(i,j)
3375 !   Tmpv006 =v(i,k,j) +Tmpv005
3376 !   Tmpv007 =Tmpv003 -Tmpv006
3377 !   Tmpv008 =rdy*Tmpv007
3378 !   Tmpv009 =muu(i+1,j)*u_1(i+1,k,j)
3379 !   Tmpv010 =Tmpv009/msfuy(i+1,j)
3380 !   Tmpv011 =u(i+1,k,j) +Tmpv010
3381 !   Tmpv012 =muu(i,j)*u_1(i,k,j)
3382 !   Tmpv013 =Tmpv012/msfuy(i,j)
3383 !   Tmpv014 =u(i,k,j) +Tmpv013
3384 !   Tmpv015 =Tmpv011 -Tmpv014
3385 !   Tmpv016 =rdx*Tmpv015
3386 !   Tmpv017 =Tmpv008 +Tmpv016
3387 !   Tmpv018 =msftx(i,j)*msfty(i,j)*Tmpv017
3388 !!  dvdxi(i,k) =Tmpv018
3389
3390 !   Tmpv001 =dmdt(i) +dnw(k)*dvdxi(i,k)
3391 !!  dmdt(i) =Tmpv001
3392
3393 !   ENDDO
3394 !   ENDDO
3395 !   DO i =i_start, i_end
3396 !!  muave(i,j) =mu(i,j)
3397
3398 !   Tmpv001 =dmdt(i) +mu_tend(i,j)
3399 !   Tmpv002 =dts*Tmpv001
3400 !   Tmpv003 =mu(i,j) +Tmpv002
3401 !!  mu(i,j) =Tmpv003
3402
3403 !   Tmpv001 =dmdt(i) +mu_tend(i,j)
3404 !!  mudf(i,j) =Tmpv001
3405
3406 !   Tmpv001 =mut(i,j) +mu(i,j)
3407 !!  muts(i,j) =Tmpv001
3408
3409 !   Tmpv001 =(1.+epssm)*mu(i,j) +(1.-epssm)*muave(i,j)
3410 !   Tmpv002 =.5*Tmpv001
3411 !!  muave(i,j) =Tmpv002
3412
3413 !   ENDDO
3414
3415 !   DO k =2, k_end
3416 !   DO i =i_start, i_end
3417 !   Tmpv001 =dmdt(i) +dvdxi(i,k-1)
3418 !   Tmpv002 =Tmpv001 +mu_tend(i,j)
3419 !   Tmpv003 =dnw(k-1)*Tmpv002
3420 !   Tmpv004 =Tmpv003/msfty(i,j)
3421 !   Tmpv005 =ww(i,k-1,j) -Tmpv004
3422 !!  ww(i,k,j) =Tmpv005
3423
3424 !   ENDDO
3425 !   ENDDO
3426 !   DO k =1, k_end
3427 !   DO i =i_start, i_end
3428 !   Tmpv001 =ww(i,k,j) -ww_1(i,k,j)
3429 !!  ww(i,k,j) =Tmpv001
3430
3431 !   ENDDO
3432 !   ENDDO
3433
3434    DO k =k_end, 1, -1
3435    DO i =i_end, i_start, -1
3436    a_Tmpv1 =a_ww(i,k,j)
3437    a_ww(i,k,j) =0.0
3438    a_ww(i,k,j) =a_ww(i,k,j) +a_Tmpv1
3439    a_ww_1(i,k,j) =a_ww_1(i,k,j) -a_Tmpv1
3440    ENDDO
3441    ENDDO
3442
3443    DO k =k_end, 2, -1
3444    DO i =i_end, i_start, -1
3445    a_Tmpv5 =a_ww(i,k,j)
3446    a_ww(i,k,j) =0.0
3447    a_ww(i,k-1,j) =a_ww(i,k-1,j) +a_Tmpv5
3448    a_Tmpv4 =-a_Tmpv5
3449    a_Tmpv3 =a_Tmpv4/msfty(i,j)
3450    a_Tmpv2 =dnw(k-1)*a_Tmpv3
3451    a_Tmpv1 =a_Tmpv2
3452    a_mu_tend(i,j) =a_mu_tend(i,j) +a_Tmpv2
3453    a_dmdt(i) =a_dmdt(i) +a_Tmpv1
3454    a_dvdxi(i,k-1) =a_dvdxi(i,k-1) +a_Tmpv1
3455    ENDDO
3456    ENDDO
3457
3458    DO i =i_end, i_start, -1
3459    a_Tmpv2 =a_muave(i,j)
3460    a_muave(i,j) =0.0
3461    a_Tmpv1 =.5*a_Tmpv2
3462    a_mu(i,j) =a_mu(i,j) +(1.+epssm)*a_Tmpv1
3463    a_muave(i,j) =a_muave(i,j) +(1.-epssm)*a_Tmpv1
3464    a_Tmpv1 =a_muts(i,j)
3465    a_muts(i,j) =0.0
3466    a_mut(i,j) =a_mut(i,j) +a_Tmpv1
3467    a_mu(i,j) =a_mu(i,j) +a_Tmpv1
3468    a_Tmpv1 =a_mudf(i,j)
3469    a_mudf(i,j) =0.0
3470    a_dmdt(i) =a_dmdt(i) +a_Tmpv1
3471    a_mu_tend(i,j) =a_mu_tend(i,j) +a_Tmpv1
3472    a_Tmpv3 =a_mu(i,j)
3473    a_mu(i,j) =0.0
3474    a_mu(i,j) =a_mu(i,j) +a_Tmpv3
3475    a_Tmpv2 =a_Tmpv3
3476    a_Tmpv1 =dts*a_Tmpv2
3477    a_dmdt(i) =a_dmdt(i) +a_Tmpv1
3478    a_mu_tend(i,j) =a_mu_tend(i,j) +a_Tmpv1
3479    a_mu(i,j) =a_mu(i,j) +a_muave(i,j)
3480    a_muave(i,j) =0.0
3481    ENDDO
3482
3483    DO k =k_end, k_start, -1
3484    DO i =i_end, i_start, -1
3485    a_Tmpv1 =a_dmdt(i)
3486    a_dmdt(i) =0.0
3487    a_dmdt(i) =a_dmdt(i) +a_Tmpv1
3488    a_dvdxi(i,k) =a_dvdxi(i,k) +dnw(k)*a_Tmpv1
3489    a_Tmpv18 =a_dvdxi(i,k)
3490    a_dvdxi(i,k) =0.0
3491    a_Tmpv17 =msftx(i,j)*msfty(i,j)*a_Tmpv18
3492    a_Tmpv8 =a_Tmpv17
3493    a_Tmpv16 =a_Tmpv17
3494    a_Tmpv15 =rdx*a_Tmpv16
3495    a_Tmpv11 =a_Tmpv15
3496    a_Tmpv14 =-a_Tmpv15
3497    a_u(i,k,j) =a_u(i,k,j) +a_Tmpv14
3498    a_Tmpv13 =a_Tmpv14
3499    a_Tmpv12 =a_Tmpv13/msfuy(i,j)
3500    a_muu(i,j) =a_muu(i,j) +u_1(i,k,j)*a_Tmpv12
3501    a_u_1(i,k,j) =a_u_1(i,k,j) +muu(i,j)*a_Tmpv12
3502    a_u(i+1,k,j) =a_u(i+1,k,j) +a_Tmpv11
3503    a_Tmpv10 =a_Tmpv11
3504    a_Tmpv9 =a_Tmpv10/msfuy(i+1,j)
3505    a_muu(i+1,j) =a_muu(i+1,j) +u_1(i+1,k,j)*a_Tmpv9
3506    a_u_1(i+1,k,j) =a_u_1(i+1,k,j) +muu(i+1,j)*a_Tmpv9
3507    a_Tmpv7 =rdy*a_Tmpv8
3508    a_Tmpv3 =a_Tmpv7
3509    a_Tmpv6 =-a_Tmpv7
3510    a_v(i,k,j) =a_v(i,k,j) +a_Tmpv6
3511    a_Tmpv5 =a_Tmpv6
3512    a_Tmpv4 =msfvx_inv(i,j)*a_Tmpv5
3513    a_muv(i,j) =a_muv(i,j) +v_1(i,k,j)*a_Tmpv4
3514    a_v_1(i,k,j) =a_v_1(i,k,j) +muv(i,j)*a_Tmpv4
3515    a_v(i,k,j+1) =a_v(i,k,j+1) +a_Tmpv3
3516    a_Tmpv2 =a_Tmpv3
3517    a_Tmpv1 =msfvx_inv(i,j+1)*a_Tmpv2
3518    a_muv(i,j+1) =a_muv(i,j+1) +v_1(i,k,j+1)*a_Tmpv1
3519    a_v_1(i,k,j+1) =a_v_1(i,k,j+1) +muv(i,j+1)*a_Tmpv1
3520    ENDDO
3521    ENDDO
3522
3523    DO i =i_end, i_start, -1
3524    a_dmdt(i) =0.0
3525    ENDDO
3526
3527    ENDDO
3528
3529 !LPB[4]
3530 !  i_endu =ite
3531 !  j_endv =jte
3532
3533 !LPB[3]
3534
3535 !  IF( config_flags%specified .or. config_flags%nested ) THEN
3536 !  j_start =max(jts, jds+1)
3537 !  j_end =min(jte, jde-2)
3538 !  ENDIF
3539
3540 ! Remarked by Ning Pan, 2010-08-31
3541 !   IF( config_flags%specified .or. config_flags%nested ) THEN
3542
3543 !   ENDIF
3544
3545 !LPB[2]
3546
3547 !LPB[1]
3548
3549 !  IF( .NOT. config_flags%periodic_x ) THEN
3550 !  IF( config_flags%specified .or. config_flags%nested ) THEN
3551 !  i_start =max(its, ids+1)
3552 !  i_end =min(ite, ide-2)
3553 !  ENDIF
3554 !  ENDIF
3555
3556 ! Remarked by Ning Pan, 2010-08-31
3557 !   IF( .NOT. config_flags%periodic_x ) THEN
3558
3559 !   IF( config_flags%specified .or. config_flags%nested ) THEN
3560
3561
3562 !   ENDIF
3563
3564 !   ENDIF
3565
3566 !LPB[0]
3567 !  i_start =its
3568 !  i_end =min(ite, ide-1)
3569 !  j_start =jts
3570 !  j_end =min(jte, jde-1)
3571 !  k_start =kts
3572 !  k_end =kte-1
3573
3574    END SUBROUTINE a_advance_mu_t
3575
3576    SUBROUTINE a_advance_w(w,a_w,rw_tend,a_rw_tend,ww,a_ww,w_save, &
3577    a_w_save,u,a_u,v,a_v,mu1,a_mu1,mut,a_mut,muave,a_muave,muts,a_muts, &
3578    t_2ave,a_t_2ave,t_2,a_t_2,t_1,a_t_1,ph,a_ph,ph_1,a_ph_1,phb, &
3579    !a_ph_tend,ph_tend,a_ht,ht,a_c2a,c2a,a_cqw,cqw,a_alt,alt,alb,a_a,a, & ! Remarked by Ning Pan, 2010-07-08
3580    ph_tend,a_ph_tend,ht,c2a,a_c2a,cqw,a_cqw,alt,a_alt,alb,a,a_a, &  ! Ning Pan, 2010-07-08
3581    alpha,a_alpha,gamma,a_gamma,rdx,rdy,dts,t0,epssm,dnw,fnm,fnp,rdnw,rdn,cf1,cf2, &
3582    cf3,msftx,msfty,config_flags,top_lid,ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme, &
3583    its,ite,jts,jte,kts,kte)
3584
3585 !PART I: DECLARATION OF VARIABLES
3586
3587    IMPLICIT NONE
3588
3589    INTEGER :: K0_ADJ,K1_ADJ,K2_ADJ,K3_ADJ
3590    TYPE(grid_config_rec_type) :: config_flags
3591    INTEGER :: ids,ide,jds,jde,kds,kde
3592    INTEGER :: ims,ime,jms,jme,kms,kme
3593    INTEGER :: its,ite,jts,jte,kts,kte
3594    LOGICAL :: top_lid
3595    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_t_2ave,t_2ave,a_w,w,a_ph,ph
3596    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_rw_tend,rw_tend,a_ww,ww,a_w_save, &
3597    w_save,a_u,u,a_v,v,a_t_2,t_2,a_t_1,t_1,a_ph_1,ph_1,phb,a_ph_tend,ph_tend, &
3598    a_alpha,alpha,a_gamma,gamma,a_a,a,a_c2a,c2a,a_cqw,cqw,alb,a_alt,alt
3599    REAL,DIMENSION(ims:ime,jms:jme) :: a_mu1,mu1,a_mut,mut,a_muave,muave,a_muts, &
3600 !   muts,a_ht,ht,msftx,msfty  ! Remarked by Ning Pan, 2010-07-09
3601    muts,a_ht,ht,msftx,msfty   ! Ning Pan, 2010-07-09
3602    REAL,DIMENSION(kms:kme) :: fnp,fnm,rdnw,rdn,dnw
3603    REAL :: rdx,rdy,dts,cf1,cf2,cf3,t0,epssm
3604    REAL,DIMENSION(its:ite) :: a_mut_inv,mut_inv,msft_inv
3605    REAL,DIMENSION(its:ite,kts:kte) :: a_rhs,rhs,a_wdwn,wdwn
3606    INTEGER :: i,j,k,i_start,i_end,j_start,j_end,k_start,k_end
3607    REAL,DIMENSION(kts:kte) :: a_dampwt,dampwt
3608    REAL :: a_htop,htop,a_hbot,hbot,hdepth,a_hk,hk
3609    REAL :: pi,dampmag
3610
3611 !  REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Keep_Lpb6_t_2ave
3612 !  REAL,DIMENSION(its:ite,kts:kte) :: Keep_Lpb6_rhs, REVISED BY WALLS
3613 !  REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Keep_Lpb6_w
3614 !  REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: Keep_Lpb6_ph
3615    REAL :: a_Tmpv1,Tmpv001,a_Tmpv2,Tmpv002,a_Tmpv3,Tmpv003,a_Tmpv4,Tmpv004, &
3616    a_Tmpv5,Tmpv005,a_Tmpv6,Tmpv006,a_Tmpv7,Tmpv007,a_Tmpv8,Tmpv008,a_Tmpv9, &
3617    Tmpv009,a_Tmpv10,Tmpv010,a_Tmpv11,Tmpv011,a_Tmpv12,Tmpv012,a_Tmpv13,Tmpv013, &
3618    a_Tmpv14,Tmpv014,a_Tmpv15,Tmpv015,a_Tmpv16,Tmpv016,a_Tmpv17,Tmpv017, &
3619    a_Tmpv18,Tmpv018,a_Tmpv19,Tmpv019,a_Tmpv20,Tmpv020,a_Tmpv21,Tmpv021, &
3620    a_Tmpv22,Tmpv022,a_Tmpv23,Tmpv023,a_Tmpv24,Tmpv024,a_Tmpv25,Tmpv025,a_Tmpv26,Tmpv026
3621 !  REAL,DIMENSION(ims:ime) :: Tmpv200, REVISED BY WALLS
3622    REAL,DIMENSION(ims:ime) :: Tmpv201
3623    REAL,DIMENSION(ims:ime) :: Tmpv202
3624    REAL,DIMENSION(ims:ime) :: Tmpv203
3625    REAL,DIMENSION(ims:ime) :: Tmpv204
3626    REAL,DIMENSION(ims:ime) :: Tmpv205
3627    REAL,DIMENSION(ims:ime) :: Tmpv206
3628    REAL,DIMENSION(ims:ime) :: Tmpv207
3629    REAL,DIMENSION(ims:ime) :: Tmpv208
3630    REAL,DIMENSION(ims:ime) :: Tmpv209
3631    REAL,DIMENSION(ims:ime) :: Tmpv2010
3632    REAL,DIMENSION(ims:ime) :: Tmpv2011
3633    REAL,DIMENSION(ims:ime) :: Tmpv2012
3634    REAL,DIMENSION(ims:ime) :: Tmpv2013
3635    REAL,DIMENSION(ims:ime) :: Tmpv2014
3636    REAL,DIMENSION(its:min(ite,ide-1),kte-1) :: Tmpv300
3637    REAL,DIMENSION(its:min(ite,ide-1),kte-1) :: Tmpv301
3638    REAL,DIMENSION(its:min(ite,ide-1),kte-1) :: Tmpv302
3639    REAL,DIMENSION(its:min(ite,ide-1),kte-1) :: Tmpv303
3640    REAL,DIMENSION(its:min(ite,ide-1),kte-1) :: Tmpv304
3641    REAL,DIMENSION(its:min(ite,ide-1),kte-1) :: Tmpv305
3642    REAL,DIMENSION(its:min(ite,ide-1),kte-1) :: Tmpv306
3643    REAL,DIMENSION(its:min(ite,ide-1),2:kte-1) :: Tmpv307
3644    REAL,DIMENSION(its:min(ite,ide-1),2:kte) :: Tmpv308
3645    REAL,DIMENSION(its:min(ite,ide-1),2:kte) :: Tmpv309
3646    REAL,DIMENSION(its:min(ite,ide-1),2:kte-1) :: Tmpv3010
3647    REAL,DIMENSION(its:min(ite,ide-1),2:kte-1) :: Tmpv3011
3648    REAL,DIMENSION(its:min(ite,ide-1),2:kte-1) :: Tmpv3012
3649    REAL,DIMENSION(its:min(ite,ide-1),2:kte-1) :: Tmpv3013
3650    REAL,DIMENSION(its:min(ite,ide-1),2:kte-1) :: Tmpv3014
3651    REAL,DIMENSION(its:min(ite,ide-1),2:kte-1) :: Tmpv3015
3652    REAL,DIMENSION(its:min(ite,ide-1),2:kte-1) :: Tmpv3016
3653    REAL,DIMENSION(its:min(ite,ide-1),2:kte) :: Tmpv3017
3654    REAL,DIMENSION(its:min(ite,ide-1),1:kte) :: Tmpv3018 !REVISED BY WALLS
3655    REAL,DIMENSION(its:min(ite,ide-1),2:kte) :: Tmpv3019 !REVISED BY WALLS
3656    REAL,DIMENSION(its:min(ite,ide-1),2:kte) :: Tmpv3020
3657    REAL,DIMENSION(its:min(ite,ide-1),2:kte) :: Tmpv3021
3658    REAL,DIMENSION(its:min(ite,ide-1),2:kte) :: Tmpv3022
3659    REAL,DIMENSION(its:min(ite,ide-1),2:kte) :: Tmpv3023
3660    REAL,DIMENSION(its:min(ite,ide-1),2:kte) :: Tmpv3024
3661    REAL,DIMENSION(its:min(ite,ide-1),2:kte) :: Tmpv3025
3662    REAL,DIMENSION(its:min(ite,ide-1),2:kte) :: Tmpv3026
3663    REAL,DIMENSION(its:min(ite,ide-1),2:kte) :: Tmpv3027
3664    REAL,DIMENSION(its:min(ite,ide-1),2:kte) :: Tmpv3028
3665    REAL,DIMENSION(its:min(ite,ide-1),2:kte) :: Tmpv3029  !REVISED BY WALLS
3666    REAL,DIMENSION(its:min(ite,ide-1),2:kte) :: Tmpv3030  !REVISED BY WALLS
3667
3668 !PART II: CALCULATIONS OF B. S. TRAJECTORY
3669
3670 !LPB[0]
3671          i_start = its
3672          i_end   = min(ite,ide-1)
3673          j_start = jts
3674          j_end   = min(jte,jde-1)
3675          k_start = kts
3676          k_end   = kte-1
3677
3678 !LPB[1]
3679       IF ( .NOT. config_flags%periodic_x )THEN
3680
3681          IF ( config_flags%specified .or. config_flags%nested ) then
3682
3683               i_start = max(its,ids+1)
3684               i_end   = min(ite,ide-2)
3685             ENDIF
3686
3687    ENDIF
3688
3689 !LPB[2]
3690
3691 !LPB[3]
3692
3693       IF ( config_flags%specified .or. config_flags%nested ) then
3694
3695             j_start = max(jts,jds+1)
3696             j_end   = min(jte,jde-2)
3697
3698    ENDIF
3699
3700 !LPB[4]
3701
3702       pi = 4.*atan(1.)
3703       dampmag = dts*config_flags%dampcoef
3704       hdepth=config_flags%zdamp
3705
3706 !LPB[5]
3707        DO i=i_start, i_end
3708
3709          rhs(i,1) = 0.
3710
3711        ENDDO
3712
3713 !LPB[6]
3714 !     j_loop_w:  DO j = j_start, j_end
3715 !
3716 !      DO k=1, k_end
3717 !      DO i=i_start, i_end
3718 !      Keep_Lpb6_t_2ave(i,k,j) =t_2ave(i,k,j)
3719 !      END DO
3720 !      END DO
3721 !      DO k=2, k_end
3722 !      Keep_Lpb6_rhs(i,k) =rhs(i,k)
3723 !      END DO
3724 !      DO k=2, k_end
3725 !      DO i=i_start, i_end
3726 !      Keep_Lpb6_w(i,k,j) =w(i,k,j)
3727 !      END DO
3728 !      END DO
3729 !      DO k=2, k_end+1
3730 !      DO i=i_start, i_end
3731 !      Keep_Lpb6_ph(i,k,j) =ph(i,k,j)
3732 !      END DO
3733 !      END DO
3734 !
3735 !       DO i=i_start, i_end
3736 !          mut_inv(i) = 1./mut(i,j)
3737 !          msft_inv(i) = 1./msfty(i,j)
3738 !       ENDDO
3739 !
3740 !       DO k=1, k_end
3741 !       DO i=i_start, i_end
3742 !         t_2ave(i,k,j)=.5*((1.+epssm)*t_2(i,k,j)         &
3743 !                       +(1.-epssm)*t_2ave(i,k,j))
3744 !         t_2ave(i,k,j)=(t_2ave(i,k,j) + muave(i,j)*t0)   &
3745 !                       /(muts(i,j)*(t0+t_1(i,k,j)))
3746 !         wdwn(i,k+1)=.5*(ww(i,k+1,j)+ww(i,k,j))*rdnw(k)      &
3747 !              *(ph_1(i,k+1,j)-ph_1(i,k,j)+phb(i,k+1,j)-phb(i,k,j))
3748 !         rhs(i,k+1) = dts*(ph_tend(i,k+1,j) + .5*g*(1.-epssm)*w(i,k+1,j))
3749 !       ENDDO
3750 !       ENDDO
3751 !
3752 !       DO k=2,k_end
3753 !       DO i=i_start, i_end
3754 !          rhs(i,k) = rhs(i,k)-dts*( fnm(k)*wdwn(i,k+1)    &
3755 !                                   +fnp(k)*wdwn(i,k  ) )
3756 !       ENDDO
3757 !       ENDDO
3758 !
3759 !       DO k=2,k_end+1
3760 !       DO i=i_start, i_end
3761 !         rhs(i,k) = ph(i,k,j) + msfty(i,j)*rhs(i,k)*mut_inv(i)
3762 !         if(top_lid .and. k.eq.k_end+1) rhs(i,k)=0.
3763 !
3764 !       ENDDO
3765 !       ENDDO
3766 !
3767 !       DO i=i_start, i_end
3768 !          w(i,1,j)=                                             &
3769 !                   msfty(i,j)*.5*rdy*(                          &
3770 !                            (ht(i,j+1)-ht(i,j  ))               &
3771 !           *(cf1*v(i,1,j+1)+cf2*v(i,2,j+1)+cf3*v(i,3,j+1))      &
3772 !                           +(ht(i,j  )-ht(i,j-1))               &
3773 !           *(cf1*v(i,1,j  )+cf2*v(i,2,j  )+cf3*v(i,3,j  ))  )   &
3774 !                  +msftx(i,j)*.5*rdx*(                          &
3775 !                            (ht(i+1,j)-ht(i,j  ))               &
3776 !           *(cf1*u(i+1,1,j)+cf2*u(i+1,2,j)+cf3*u(i+1,3,j))      &
3777 !                           +(ht(i,j  )-ht(i-1,j))               &
3778 !           *(cf1*u(i  ,1,j)+cf2*u(i  ,2,j)+cf3*u(i  ,3,j))  )
3779 !        ENDDO
3780 !
3781 !       DO k=2,k_end
3782 !         DO i=i_start, i_end
3783 !           w(i,k,j)=w(i,k,j)+dts*rw_tend(i,k,j)                         &
3784 !                    + msft_inv(i)*cqw(i,k,j)*(                          &
3785 !               +.5*dts*g*mut_inv(i)*rdn(k)*                             &
3786 !                (c2a(i,k  ,j)*rdnw(k  )                                 &
3787 !           *((1.+epssm)*(rhs(i,k+1  )-rhs(i,k    ))                     &
3788 !            +(1.-epssm)*(ph(i,k+1,j)-ph(i,k  ,j)))                      &
3789 !                -c2a(i,k-1,j)*rdnw(k-1)                                 &
3790 !           *((1.+epssm)*(rhs(i,k    )-rhs(i,k-1  ))                     &
3791 !            +(1.-epssm)*(ph(i,k  ,j)-ph(i,k-1,j)))))                    &
3792 !                   +dts*g*msft_inv(i)*(rdn(k)*                          &
3793 !                (c2a(i,k  ,j)*alt(i,k  ,j)*t_2ave(i,k  ,j)              &
3794 !                -c2a(i,k-1,j)*alt(i,k-1,j)*t_2ave(i,k-1,j))             &
3795 !                 - muave(i,j))
3796 !         ENDDO
3797 !       ENDDO
3798 !       K=k_end+1
3799 !
3800 !       DO i=i_start, i_end
3801 !          w(i,k,j)=w(i,k,j)+dts*rw_tend(i,k,j)                           &
3802 !              +msft_inv(i)*(                                          &
3803 !            -.5*dts*g*mut_inv(i)*rdnw(k-1)**2*2.*c2a(i,k-1,j)              &
3804 !                *((1.+epssm)*(rhs(i,k  )-rhs(i,k-1  ))                   &
3805 !                 +(1.-epssm)*(ph(i,k,j)-ph(i,k-1,j)))                    &
3806 !            -dts*g*(2.*rdnw(k-1)*                                        &
3807 !                      c2a(i,k-1,j)*alt(i,k-1,j)*t_2ave(i,k-1,j)          &
3808 !                 + muave(i,j)) )
3809 !       if(top_lid)w(i,k,j) = 0.
3810 !
3811 !       ENDDO
3812 !
3813 !       DO k=2,k_end+1
3814 !       DO i=i_start, i_end
3815 !          w(i,k,j)=(w(i,k,j)-a(i,k,j)*w(i,k-1,j))*alpha(i,k,j)
3816 !       ENDDO
3817 !      ENDDO
3818 !
3819 !      DO k=k_end,2,-1
3820 !       DO i=i_start, i_end
3821 !          w (i,k,j)=w (i,k,j)-gamma(i,k,j)*w(i,k+1,j)
3822 !       ENDDO
3823 !       ENDDO
3824 !      IF (config_flags%damp_opt .eq. 3) THEN
3825 !
3826 !        DO k=k_end+1,2,-1
3827 !         DO i=i_start, i_end
3828 !             htop=(ph_1(i,k_end+1,j)+phb(i,k_end+1,j))/g
3829 !             hk=(ph_1(i,k,j)+phb(i,k,j))/g
3830 !             hbot=htop-hdepth
3831 !             dampwt(k) = 0.
3832 !          if(hk .ge. hbot)then
3833 !
3834 !               dampwt(k) = dampmag*sin(0.5*pi*(hk-hbot)/hdepth)*sin(0.5*pi*(hk-hbot)/hdepth)
3835 !             endif
3836 !             w(i,k,j) = (w(i,k,j) - dampwt(k)*mut(i,j)*w_save(i,k,j))/(1.+dampwt(k))
3837 !         ENDDO
3838 !         ENDDO
3839 !       ENDIF
3840 !
3841 !       DO k=k_end+1,2,-1
3842 !       DO i=i_start, i_end
3843 !          ph(i,k,j) = rhs(i,k)+msfty(i,j)*.5*dts*g*(1.+epssm)   &
3844 !                         *w(i,k,j)/muts(i,j)
3845 !       ENDDO
3846 !       ENDDO
3847 !
3848 !     ENDDO j_loop_w
3849
3850 !PART III: INITIALIZATION OF LOCAL ADJOINT PERTURBATIONS
3851
3852    Do K0_ADJ =its, ite
3853    a_mut_inv(K0_ADJ) =0.0
3854    End Do
3855
3856    Do K1_ADJ =kts, kte
3857    Do K0_ADJ =its, ite
3858    a_rhs(K0_ADJ,K1_ADJ) =0.0
3859    End Do
3860    End Do
3861
3862    Do K1_ADJ =kts, kte
3863    Do K0_ADJ =its, ite
3864    a_wdwn(K0_ADJ,K1_ADJ) =0.0
3865    End Do
3866    End Do
3867
3868    Do K0_ADJ =kts, kte
3869    a_dampwt(K0_ADJ) =0.0
3870    End Do
3871
3872    a_htop =0.0
3873    a_hbot =0.0
3874    a_hk =0.0
3875
3876 !PART IV: REVERSE/BACKWARD ACCUMULATIONS
3877
3878 !LPB[6]
3879    DO j =j_end, j_start, -1
3880
3881 !  DO k=1, k_end
3882 !  DO i=i_start, i_end
3883 !  t_2ave(i,k,j) =Keep_Lpb6_t_2ave(i,k,j)
3884 !  END DO
3885 !  END DO
3886 !  DO k=2, k_end
3887 !  rhs(i,k) =Keep_Lpb6_rhs(i,k)
3888 !  END DO
3889 !  DO k=2, k_end
3890 !  DO i=i_start, i_end
3891 !  w(i,k,j) =Keep_Lpb6_w(i,k,j)
3892 !  END DO
3893 !  END DO
3894 !  DO k=2, k_end+1
3895 !  DO i=i_start, i_end
3896 !  ph(i,k,j) =Keep_Lpb6_ph(i,k,j)
3897 !  END DO
3898 !  END DO
3899
3900    DO i =i_start, i_end
3901 !THE KEEPING OPERATION IS NOT NECESSARY IF THE VARIABLE IS LOCALLY AN OUTPUT VARIABLE
3902 !REVISED BY WALLS
3903 !  Tmpv200(i) =mut_inv(i)
3904    mut_inv(i) =1./mut(i,j)
3905
3906    msft_inv(i) =1./msfty(i,j)
3907    ENDDO
3908
3909    DO k =1, k_end
3910    DO i =i_start, i_end
3911    Tmpv001 =(1.+epssm)*t_2(i,k,j) +(1.-epssm)*t_2ave(i,k,j)
3912    Tmpv002 =.5*Tmpv001
3913    Tmpv300(i,k) =t_2ave(i,k,j)
3914    t_2ave(i,k,j) =Tmpv002
3915
3916    Tmpv001 =t_2ave(i,k,j) +muave(i,j)*t0
3917    Tmpv002 =muts(i,j)*(t0 +t_1(i,k,j))
3918    Tmpv301(i,k) =Tmpv001
3919    Tmpv302(i,k) =Tmpv002
3920    Tmpv003 =Tmpv301(i,k)/Tmpv302(i,k)
3921    Tmpv303(i,k) =t_2ave(i,k,j)
3922    t_2ave(i,k,j) =Tmpv003
3923
3924    Tmpv001 =ww(i,k+1,j) +ww(i,k,j)
3925    Tmpv002 =.5*Tmpv001
3926    Tmpv003 =Tmpv002*rdnw(k)
3927    Tmpv004 =ph_1(i,k+1,j) -ph_1(i,k,j)
3928    Tmpv005 =Tmpv004 +phb(i,k+1,j)
3929    Tmpv006 =Tmpv005 -phb(i,k,j)
3930    Tmpv304(i,k) =Tmpv003
3931    Tmpv305(i,k) =Tmpv006
3932    Tmpv007 =Tmpv304(i,k)*Tmpv305(i,k)
3933   wdwn(i,k+1) =Tmpv007  ! Removed comment by Ning Pan, 2010-07-09
3934
3935    Tmpv001 =ph_tend(i,k+1,j) +.5*g*(1.-epssm)*w(i,k+1,j)
3936    Tmpv002 =dts*Tmpv001
3937    Tmpv306(i,k) =rhs(i,k+1)
3938    rhs(i,k+1) =Tmpv002
3939
3940    ENDDO
3941    ENDDO
3942
3943    DO k =2, k_end
3944    DO i =i_start, i_end
3945    Tmpv001 =fnm(k)*wdwn(i,k+1) +fnp(k)*wdwn(i,k)
3946    Tmpv002 =dts*Tmpv001
3947    Tmpv003 =rhs(i,k) -Tmpv002
3948    Tmpv307(i,k) =rhs(i,k)
3949    rhs(i,k) =Tmpv003
3950
3951    ENDDO
3952    ENDDO
3953
3954    DO k =2, k_end+1
3955    DO i =i_start, i_end
3956    Tmpv001 =msfty(i,j)*rhs(i,k)*mut_inv(i)
3957    Tmpv002 =ph(i,k,j) +Tmpv001
3958    Tmpv308(i,k) =rhs(i,k)
3959    rhs(i,k) =Tmpv002
3960
3961    IF(top_lid .and. k.eq.k_end+1) THEN
3962    Tmpv309(i,k) =rhs(i,k)
3963    rhs(i,k) =0.
3964    END IF
3965    ENDDO
3966    ENDDO
3967
3968    DO i =i_start, i_end
3969    Tmpv001 =ht(i,j+1) -ht(i,j)
3970    Tmpv002 =cf1*v(i,1,j+1) +cf2*v(i,2,j+1)
3971    Tmpv003 =Tmpv002 +cf3*v(i,3,j+1)
3972    Tmpv201(i) =Tmpv001
3973    Tmpv202(i) =Tmpv003
3974    Tmpv004 =Tmpv201(i)*Tmpv202(i)
3975    Tmpv005 =ht(i,j) -ht(i,j-1)
3976    Tmpv006 =cf1*v(i,1,j) +cf2*v(i,2,j)
3977    Tmpv007 =Tmpv006 +cf3*v(i,3,j)
3978    Tmpv203(i) =Tmpv005
3979    Tmpv204(i) =Tmpv007
3980    Tmpv008 =Tmpv203(i)*Tmpv204(i)
3981    Tmpv009 =Tmpv004 +Tmpv008
3982    Tmpv010 =msfty(i,j)*.5*rdy*Tmpv009
3983    Tmpv011 =ht(i+1,j) -ht(i,j)
3984    Tmpv012 =cf1*u(i+1,1,j) +cf2*u(i+1,2,j)
3985    Tmpv013 =Tmpv012 +cf3*u(i+1,3,j)
3986    Tmpv205(i) =Tmpv011
3987    Tmpv206(i) =Tmpv013
3988    Tmpv014 =Tmpv205(i)*Tmpv206(i)
3989    Tmpv015 =ht(i,j) -ht(i-1,j)
3990    Tmpv016 =cf1*u(i,1,j) +cf2*u(i,2,j)
3991    Tmpv017 =Tmpv016 +cf3*u(i,3,j)
3992    Tmpv207(i) =Tmpv015
3993    Tmpv208(i) =Tmpv017
3994    Tmpv018 =Tmpv207(i)*Tmpv208(i)
3995    Tmpv019 =Tmpv014 +Tmpv018
3996    Tmpv020 =msftx(i,j)*.5*rdx*Tmpv019
3997    Tmpv021 =Tmpv010 +Tmpv020
3998    Tmpv209(i) =w(i,1,j)
3999    w(i,1,j) =Tmpv021
4000
4001    ENDDO
4002
4003    DO k =2, k_end
4004    DO i =i_start, i_end
4005    Tmpv001 =w(i,k,j) +dts*rw_tend(i,k,j)
4006    Tmpv002 =rhs(i,k+1) -rhs(i,k)
4007    Tmpv003 =(1.+epssm)*Tmpv002
4008    Tmpv004 =ph(i,k+1,j) -ph(i,k,j)
4009    Tmpv005 =(1.-epssm)*Tmpv004
4010    Tmpv006 =Tmpv003 +Tmpv005
4011    Tmpv3010(i,k) =Tmpv006
4012    Tmpv007 =c2a(i,k,j)*rdnw(k)*Tmpv3010(i,k)
4013    Tmpv008 =rhs(i,k) -rhs(i,k-1)
4014    Tmpv009 =(1.+epssm)*Tmpv008
4015    Tmpv010 =ph(i,k,j) -ph(i,k-1,j)
4016    Tmpv011 =(1.-epssm)*Tmpv010
4017    Tmpv012 =Tmpv009 +Tmpv011
4018    Tmpv3011(i,k) =Tmpv012
4019    Tmpv013 =c2a(i,k-1,j)*rdnw(k-1)*Tmpv3011(i,k)
4020    Tmpv014 =Tmpv007 -Tmpv013
4021    Tmpv3012(i,k) =Tmpv014
4022    Tmpv015 =.5*dts*g*mut_inv(i)*rdn(k)*Tmpv3012(i,k)
4023    Tmpv3013(i,k) =+Tmpv015
4024    Tmpv016 =msft_inv(i)*cqw(i,k,j)*Tmpv3013(i,k)
4025    Tmpv017 =Tmpv001 +Tmpv016
4026    Tmpv018 =c2a(i,k,j)*alt(i,k,j)
4027    Tmpv3014(i,k) =Tmpv018
4028    Tmpv019 =Tmpv3014(i,k)*t_2ave(i,k,j)
4029    Tmpv020 =c2a(i,k-1,j)*alt(i,k-1,j)
4030    Tmpv3015(i,k) =Tmpv020
4031    Tmpv021 =Tmpv3015(i,k)*t_2ave(i,k-1,j)
4032    Tmpv022 =Tmpv019 -Tmpv021
4033    Tmpv023 =rdn(k)*Tmpv022
4034    Tmpv024 =Tmpv023 -muave(i,j)
4035    Tmpv025 =dts*g*msft_inv(i)*Tmpv024
4036    Tmpv026 =Tmpv017 +Tmpv025
4037    Tmpv3016(i,k) =w(i,k,j)
4038    w(i,k,j) =Tmpv026
4039
4040    ENDDO
4041    ENDDO
4042    K =k_end+1
4043    DO i =i_start, i_end
4044 !DELETE JUST FOR TUNNING
4045    Tmpv001 =w(i,k,j) +dts*rw_tend(i,k,j)
4046    Tmpv002 =-.5*dts*g*mut_inv(i)*rdnw(k-1)**2*2.*c2a(i,k-1,j)
4047    Tmpv003 =rhs(i,k) -rhs(i,k-1)
4048    Tmpv004 =(1.+epssm)*Tmpv003
4049    Tmpv005 =ph(i,k,j) -ph(i,k-1,j)
4050    Tmpv006 =(1.-epssm)*Tmpv005
4051    Tmpv007 =Tmpv004 +Tmpv006
4052    Tmpv2010(i) =Tmpv002
4053    Tmpv2011(i) =Tmpv007
4054    Tmpv008 =Tmpv2010(i)*Tmpv2011(i)
4055    Tmpv009 =2.*rdnw(k-1)*c2a(i,k-1,j)*alt(i,k-1,j)
4056    Tmpv2012(i) =Tmpv009
4057    Tmpv010 =Tmpv2012(i)*t_2ave(i,k-1,j)
4058    Tmpv011 =Tmpv010 +muave(i,j)
4059    Tmpv012 =dts*g*Tmpv011
4060    Tmpv013 =Tmpv008 -Tmpv012
4061    Tmpv014 =msft_inv(i)*Tmpv013
4062    Tmpv015 =Tmpv001 +Tmpv014
4063    Tmpv2013(i) =w(i,k,j)
4064    w(i,k,j) =Tmpv015
4065
4066    IF(top_lid) THEN
4067    Tmpv2014(i) =w(i,k,j)
4068    w(i,k,j) =0.
4069
4070    END IF
4071    ENDDO
4072
4073    DO k =2, k_end+1
4074    DO i =i_start, i_end
4075    Tmpv001 =a(i,k,j)*w(i,k-1,j)
4076    Tmpv002 =w(i,k,j) -Tmpv001
4077    Tmpv3017(i,k) =Tmpv002
4078    Tmpv003 =Tmpv3017(i,k)*alpha(i,k,j)
4079 !   Tmpv3018(i,k) =w(i,k,j)
4080    Tmpv3018(i,k) =w(i,k-1,j)  !REVISED BY WALLS
4081    w(i,k,j) =Tmpv003
4082
4083    ENDDO
4084    ENDDO
4085    DO k =k_end, 2, -1
4086    DO i =i_start, i_end
4087    Tmpv001 =gamma(i,k,j)*w(i,k+1,j)
4088    Tmpv002 =w(i,k,j) -Tmpv001
4089 !   Tmpv3019(i,k) =w(i,k,j)
4090    Tmpv3019(i,k) =w(i,k+1,j)  !REVISED BY WALLS
4091    w(i,k,j) =Tmpv002
4092
4093    ENDDO
4094    ENDDO
4095
4096    IF(config_flags%damp_opt .eq. 3) THEN
4097    DO k =k_end+1, 2, -1
4098    DO i =i_start, i_end
4099    htop =(ph_1(i,k_end+1,j) +phb(i,k_end+1,j))/g  !REVISED BY WALLS
4100
4101    hk =(ph_1(i,k,j) +phb(i,k,j))/g
4102
4103    hbot =htop -hdepth
4104
4105    Tmpv3020(i,k) =dampwt(k)
4106    dampwt(k) =0.
4107
4108    Tmpv3029(i,k) =hk  !REVISED BY WALLS
4109    Tmpv3030(i,k) =hbot  !REVISED BY WALLS
4110
4111    IF(hk .ge. hbot) THEN
4112    Tmpv001 =hk -hbot
4113    Tmpv002 =0.5*pi*Tmpv001
4114    Tmpv003 =Tmpv002/hdepth
4115    Tmpv3021(i,k) =Tmpv003
4116    Tmpv004 =sin(Tmpv3021(i,k))
4117    Tmpv005 =dampmag*Tmpv004
4118    Tmpv006 =hk -hbot
4119    Tmpv007 =0.5*pi*Tmpv006
4120    Tmpv008 =Tmpv007/hdepth
4121    Tmpv3022(i,k) =Tmpv008
4122    Tmpv009 =sin(Tmpv3022(i,k))
4123    Tmpv3023(i,k) =Tmpv005
4124    Tmpv3024(i,k) =Tmpv009
4125    Tmpv010 =Tmpv3023(i,k)*Tmpv3024(i,k)
4126    Tmpv3025(i,k) =dampwt(k)
4127    dampwt(k) =Tmpv010
4128
4129    endif
4130    Tmpv001 =dampwt(k)*mut(i,j)
4131    Tmpv3026(i,k) =Tmpv001
4132    Tmpv002 =Tmpv3026(i,k)*w_save(i,k,j)
4133    Tmpv003 =w(i,k,j) -Tmpv002
4134    Tmpv3027(i,k) =Tmpv003
4135    Tmpv004 =Tmpv3027(i,k)/(1. +dampwt(k))
4136 !  Tmpv3028(i,k) =w(i,k,j), REVISED BY WALLS
4137    Tmpv3028(i,k) =dampwt(k) !REVISED BY WALLS
4138    w(i,k,j) =Tmpv004
4139
4140    ENDDO
4141    ENDDO
4142    ENDIF
4143
4144 ! Ning Pan, 2010-07-09 - No need to recalculate ph
4145 !   DO k =k_end+1, 2, -1
4146 !   DO i =i_start, i_end
4147 !DELETE! JUST FOR TUNNING
4148 !!   Tmpv001 =msfty(i,j)*.5*dts*g*(1.+epssm)*w(i,k,j)/muts(i,j)
4149 !!   Tmpv002 =rhs(i,k) +Tmpv001
4150 !!   ph(i,k,j) =Tmpv002
4151
4152 !   ENDDO
4153 !   ENDDO
4154
4155    DO k =2, k_end+1, 1
4156    DO i =i_end, i_start, -1
4157 !DELETE JUST FOR TUNNING
4158    a_Tmpv2 =a_ph(i,k,j)
4159    a_ph(i,k,j) =0.0
4160    a_rhs(i,k) =a_rhs(i,k) +a_Tmpv2
4161    a_Tmpv1 =a_Tmpv2
4162    a_w(i,k,j) =a_w(i,k,j) +msfty(i,j)*.5*dts*g*(1.+epssm)/muts(i,j)*a_Tmpv1
4163    a_muts(i,j) =a_muts(i,j) -msfty(i,j)*.5*dts*g*(1.+epssm)*w(i,k,j)/(muts(i,j)  &
4164    *muts(i,j))*a_Tmpv1
4165    ENDDO
4166    ENDDO
4167
4168    IF(config_flags%damp_opt .eq. 3) THEN
4169
4170    DO k =2, k_end+1, 1
4171    DO i =i_end, i_start, -1
4172
4173 !  w(i,k,j) =Tmpv3028(i,k), REVISED BY WALLS
4174
4175    a_Tmpv4 =a_w(i,k,j)
4176    a_w(i,k,j) =0.0
4177    dampwt(k) =Tmpv3028(i,k) !REVISED BY WALLS
4178    a_Tmpv3 =a_Tmpv4/(1. +dampwt(k))
4179    a_dampwt(k) =a_dampwt(k) -Tmpv3027(i,k)/((1. +dampwt(k))*(1. +dampwt(k)))*a_Tmpv4
4180    a_w(i,k,j) =a_w(i,k,j) +a_Tmpv3
4181    a_Tmpv2 =-a_Tmpv3
4182    a_Tmpv1 =w_save(i,k,j)*a_Tmpv2
4183    a_w_save(i,k,j) =a_w_save(i,k,j) +Tmpv3026(i,k)*a_Tmpv2
4184    a_dampwt(k) =a_dampwt(k) +mut(i,j)*a_Tmpv1
4185    a_mut(i,j) =a_mut(i,j) +dampwt(k)*a_Tmpv1
4186
4187    hk =Tmpv3029(i,k)  !REVISED BY WALLS
4188    hbot =Tmpv3030(i,k)  !REVISED BY WALLS
4189    IF(hk .ge. hbot) THEN
4190
4191    dampwt(k) =Tmpv3025(i,k)
4192
4193    a_Tmpv10 =a_dampwt(k)
4194    a_dampwt(k) =0.0
4195    a_Tmpv5 =Tmpv3024(i,k)*a_Tmpv10
4196    a_Tmpv9 =Tmpv3023(i,k)*a_Tmpv10
4197    a_Tmpv8 =cos(Tmpv3022(i,k))*a_Tmpv9
4198    a_Tmpv7 =a_Tmpv8/hdepth
4199    a_Tmpv6 =0.5*pi*a_Tmpv7
4200    a_hk =a_hk +a_Tmpv6
4201    a_hbot =a_hbot -a_Tmpv6
4202    a_Tmpv4 =dampmag*a_Tmpv5
4203    a_Tmpv3 =cos(Tmpv3021(i,k))*a_Tmpv4
4204    a_Tmpv2 =a_Tmpv3/hdepth
4205    a_Tmpv1 =0.5*pi*a_Tmpv2
4206    a_hk =a_hk +a_Tmpv1
4207    a_hbot =a_hbot -a_Tmpv1
4208
4209    ENDIF !REVISED BY WALLS
4210 !  endif
4211
4212    dampwt(k) =Tmpv3020(i,k)
4213
4214    a_dampwt(k) =0.0
4215    a_htop =a_htop +a_hbot
4216    a_hbot =0.0
4217    a_ph_1(i,k,j) =a_ph_1(i,k,j) +1.0/g*a_hk
4218    a_hk =0.0
4219    a_ph_1(i,k_end+1,j) =a_ph_1(i,k_end+1,j) +1.0/g*a_htop
4220    a_htop =0.0
4221    ENDDO
4222    ENDDO
4223
4224    ENDIF
4225
4226    DO k =2, k_end, 1
4227    DO i =i_end, i_start, -1
4228
4229 !  w(i,k,j) =Tmpv3019(i,k)
4230    w(i,k+1,j) =Tmpv3019(i,k)  !REVISED BY WALLS
4231    a_Tmpv2 =a_w(i,k,j)
4232    a_w(i,k,j) =0.0
4233    a_w(i,k,j) =a_w(i,k,j) +a_Tmpv2
4234    a_Tmpv1 =-a_Tmpv2
4235    a_gamma(i,k,j) =a_gamma(i,k,j) +w(i,k+1,j)*a_Tmpv1
4236    a_w(i,k+1,j) =a_w(i,k+1,j) +gamma(i,k,j)*a_Tmpv1
4237    ENDDO
4238    ENDDO
4239
4240    DO k =k_end+1, 2, -1
4241    DO i =i_end, i_start, -1
4242
4243 !   w(i,k,j) =Tmpv3018(i,k)
4244    w(i,k-1,j) =Tmpv3018(i,k)
4245
4246    a_Tmpv3 =a_w(i,k,j)
4247    a_w(i,k,j) =0.0
4248    a_Tmpv2 =alpha(i,k,j)*a_Tmpv3
4249    a_alpha(i,k,j) =a_alpha(i,k,j) +Tmpv3017(i,k)*a_Tmpv3
4250    a_w(i,k,j) =a_w(i,k,j) +a_Tmpv2
4251    a_Tmpv1 =-a_Tmpv2
4252    a_a(i,k,j) =a_a(i,k,j) +w(i,k-1,j)*a_Tmpv1
4253    a_w(i,k-1,j) =a_w(i,k-1,j) +a(i,k,j)*a_Tmpv1
4254    ENDDO
4255    ENDDO
4256
4257    K=k_end+1      ! Added by Ning Pan, 2010-07-09
4258    DO i =i_end, i_start, -1
4259
4260 !DELETE JUST FOR TUNNING
4261    IF(top_lid) THEN
4262
4263    w(i,k,j) =Tmpv2014(i)
4264    a_w(i,k,j) =0.0
4265
4266    END IF
4267
4268    w(i,k,j) =Tmpv2013(i)
4269
4270    a_Tmpv15 =a_w(i,k,j)
4271    a_w(i,k,j) =0.0
4272    a_Tmpv1 =a_Tmpv15
4273    a_Tmpv14 =a_Tmpv15
4274    a_Tmpv13 =msft_inv(i)*a_Tmpv14
4275    a_Tmpv8 =a_Tmpv13
4276    a_Tmpv12 =-a_Tmpv13
4277    a_Tmpv11 =dts*g*a_Tmpv12
4278    a_Tmpv10 =a_Tmpv11
4279    a_muave(i,j) =a_muave(i,j) +a_Tmpv11
4280    a_Tmpv9 =t_2ave(i,k-1,j)*a_Tmpv10
4281    a_t_2ave(i,k-1,j) =a_t_2ave(i,k-1,j) +Tmpv2012(i)*a_Tmpv10
4282    a_c2a(i,k-1,j) =a_c2a(i,k-1,j) +2.*rdnw(k-1)*alt(i,k-1,j)*a_Tmpv9
4283    a_alt(i,k-1,j) =a_alt(i,k-1,j) +2.*rdnw(k-1)*c2a(i,k-1,j)*a_Tmpv9
4284    a_Tmpv2 =Tmpv2011(i)*a_Tmpv8
4285    a_Tmpv7 =Tmpv2010(i)*a_Tmpv8
4286    a_Tmpv4 =a_Tmpv7
4287    a_Tmpv6 =a_Tmpv7
4288    a_Tmpv5 =(1.-epssm)*a_Tmpv6
4289    a_ph(i,k,j) =a_ph(i,k,j) +a_Tmpv5
4290    a_ph(i,k-1,j) =a_ph(i,k-1,j) -a_Tmpv5
4291    a_Tmpv3 =(1.+epssm)*a_Tmpv4
4292    a_rhs(i,k) =a_rhs(i,k) +a_Tmpv3
4293    a_rhs(i,k-1) =a_rhs(i,k-1) -a_Tmpv3
4294    a_mut_inv(i) =a_mut_inv(i) -.5*dts*g*rdnw(k-1)**2*2.*c2a(i,k-1,j)*a_Tmpv2
4295    a_c2a(i,k-1,j) =a_c2a(i,k-1,j) -.5*dts*g*mut_inv(i)*rdnw(k-1)**2*2.*a_Tmpv2
4296    a_w(i,k,j) =a_w(i,k,j) +a_Tmpv1
4297    a_rw_tend(i,k,j) =a_rw_tend(i,k,j) +dts*a_Tmpv1
4298    ENDDO
4299
4300 !  K=k_end+1, REVISED BY WALLS
4301
4302    DO k =k_end, 2, -1
4303    DO i =i_end, i_start, -1
4304
4305    w(i,k,j) =Tmpv3016(i,k)
4306
4307    a_Tmpv26 =a_w(i,k,j)
4308    a_w(i,k,j) =0.0
4309    a_Tmpv17 =a_Tmpv26
4310    a_Tmpv25 =a_Tmpv26
4311    a_Tmpv24 =dts*g*msft_inv(i)*a_Tmpv25
4312    a_Tmpv23 =a_Tmpv24
4313    a_muave(i,j) =a_muave(i,j) -a_Tmpv24
4314    a_Tmpv22 =rdn(k)*a_Tmpv23
4315    a_Tmpv19 =a_Tmpv22
4316    a_Tmpv21 =-a_Tmpv22
4317    a_Tmpv20 =t_2ave(i,k-1,j)*a_Tmpv21
4318    a_t_2ave(i,k-1,j) =a_t_2ave(i,k-1,j) +Tmpv3015(i,k)*a_Tmpv21
4319    a_c2a(i,k-1,j) =a_c2a(i,k-1,j) +alt(i,k-1,j)*a_Tmpv20
4320    a_alt(i,k-1,j) =a_alt(i,k-1,j) +c2a(i,k-1,j)*a_Tmpv20
4321    a_Tmpv18 =t_2ave(i,k,j)*a_Tmpv19
4322    a_t_2ave(i,k,j) =a_t_2ave(i,k,j) +Tmpv3014(i,k)*a_Tmpv19
4323    a_c2a(i,k,j) =a_c2a(i,k,j) +alt(i,k,j)*a_Tmpv18
4324    a_alt(i,k,j) =a_alt(i,k,j) +c2a(i,k,j)*a_Tmpv18
4325    a_Tmpv1 =a_Tmpv17
4326    a_Tmpv16 =a_Tmpv17
4327    a_cqw(i,k,j) =a_cqw(i,k,j) +msft_inv(i)*Tmpv3013(i,k)*a_Tmpv16
4328    a_Tmpv15 =msft_inv(i)*cqw(i,k,j)*a_Tmpv16
4329    a_mut_inv(i) =a_mut_inv(i) +.5*dts*g*rdn(k)*Tmpv3012(i,k)*a_Tmpv15
4330    a_Tmpv14 =.5*dts*g*mut_inv(i)*rdn(k)*a_Tmpv15
4331    a_Tmpv7 =a_Tmpv14
4332    a_Tmpv13 =-a_Tmpv14
4333    a_c2a(i,k-1,j) =a_c2a(i,k-1,j) +rdnw(k-1)*Tmpv3011(i,k)*a_Tmpv13
4334    a_Tmpv12 =c2a(i,k-1,j)*rdnw(k-1)*a_Tmpv13
4335    a_Tmpv9 =a_Tmpv12
4336    a_Tmpv11 =a_Tmpv12
4337    a_Tmpv10 =(1.-epssm)*a_Tmpv11
4338    a_ph(i,k,j) =a_ph(i,k,j) +a_Tmpv10
4339    a_ph(i,k-1,j) =a_ph(i,k-1,j) -a_Tmpv10
4340    a_Tmpv8 =(1.+epssm)*a_Tmpv9
4341    a_rhs(i,k) =a_rhs(i,k) +a_Tmpv8
4342    a_rhs(i,k-1) =a_rhs(i,k-1) -a_Tmpv8
4343    a_c2a(i,k,j) =a_c2a(i,k,j) +rdnw(k)*Tmpv3010(i,k)*a_Tmpv7
4344    a_Tmpv6 =c2a(i,k,j)*rdnw(k)*a_Tmpv7
4345    a_Tmpv3 =a_Tmpv6
4346    a_Tmpv5 =a_Tmpv6
4347    a_Tmpv4 =(1.-epssm)*a_Tmpv5
4348    a_ph(i,k+1,j) =a_ph(i,k+1,j) +a_Tmpv4
4349    a_ph(i,k,j) =a_ph(i,k,j) -a_Tmpv4
4350    a_Tmpv2 =(1.+epssm)*a_Tmpv3
4351    a_rhs(i,k+1) =a_rhs(i,k+1) +a_Tmpv2
4352    a_rhs(i,k) =a_rhs(i,k) -a_Tmpv2
4353    a_w(i,k,j) =a_w(i,k,j) +a_Tmpv1
4354    a_rw_tend(i,k,j) =a_rw_tend(i,k,j) +dts*a_Tmpv1
4355    ENDDO
4356    ENDDO
4357
4358    DO i =i_end, i_start, -1
4359
4360    w(i,1,j) =Tmpv209(i)
4361
4362    a_Tmpv21 =a_w(i,1,j)
4363    a_w(i,1,j) =0.0
4364    a_Tmpv10 =a_Tmpv21
4365    a_Tmpv20 =a_Tmpv21
4366    a_Tmpv19 =msftx(i,j)*.5*rdx*a_Tmpv20
4367    a_Tmpv14 =a_Tmpv19
4368    a_Tmpv18 =a_Tmpv19
4369    a_Tmpv15 =Tmpv208(i)*a_Tmpv18
4370    a_Tmpv17 =Tmpv207(i)*a_Tmpv18
4371    a_Tmpv16 =a_Tmpv17
4372    a_u(i,3,j) =a_u(i,3,j) +cf3*a_Tmpv17
4373    a_u(i,1,j) =a_u(i,1,j) +cf1*a_Tmpv16
4374    a_u(i,2,j) =a_u(i,2,j) +cf2*a_Tmpv16
4375 !   a_ht(i,j) =a_ht(i,j) +a_Tmpv15     ! Remarked by Ning Pan, 2010-07-09
4376 !   a_ht(i-1,j) =a_ht(i-1,j) -a_Tmpv15 ! Remarked by Ning Pan, 2010-07-09
4377    a_Tmpv11 =Tmpv206(i)*a_Tmpv14
4378    a_Tmpv13 =Tmpv205(i)*a_Tmpv14
4379    a_Tmpv12 =a_Tmpv13
4380    a_u(i+1,3,j) =a_u(i+1,3,j) +cf3*a_Tmpv13
4381    a_u(i+1,1,j) =a_u(i+1,1,j) +cf1*a_Tmpv12
4382    a_u(i+1,2,j) =a_u(i+1,2,j) +cf2*a_Tmpv12
4383 !   a_ht(i+1,j) =a_ht(i+1,j) +a_Tmpv11 ! Remarked by Ning Pan, 2010-07-09
4384 !   a_ht(i,j) =a_ht(i,j) -a_Tmpv11     ! Remarked by Ning Pan, 2010-07-09
4385    a_Tmpv9 =msfty(i,j)*.5*rdy*a_Tmpv10
4386    a_Tmpv4 =a_Tmpv9
4387    a_Tmpv8 =a_Tmpv9
4388    a_Tmpv5 =Tmpv204(i)*a_Tmpv8
4389    a_Tmpv7 =Tmpv203(i)*a_Tmpv8
4390    a_Tmpv6 =a_Tmpv7
4391    a_v(i,3,j) =a_v(i,3,j) +cf3*a_Tmpv7
4392    a_v(i,1,j) =a_v(i,1,j) +cf1*a_Tmpv6
4393    a_v(i,2,j) =a_v(i,2,j) +cf2*a_Tmpv6
4394 !   a_ht(i,j) =a_ht(i,j) +a_Tmpv5     ! Remarked by Ning Pan, 2010-07-09
4395 !   a_ht(i,j-1) =a_ht(i,j-1) -a_Tmpv5 ! Remarked by Ning Pan, 2010-07-09
4396    a_Tmpv1 =Tmpv202(i)*a_Tmpv4
4397    a_Tmpv3 =Tmpv201(i)*a_Tmpv4
4398    a_Tmpv2 =a_Tmpv3
4399    a_v(i,3,j+1) =a_v(i,3,j+1) +cf3*a_Tmpv3
4400    a_v(i,1,j+1) =a_v(i,1,j+1) +cf1*a_Tmpv2
4401    a_v(i,2,j+1) =a_v(i,2,j+1) +cf2*a_Tmpv2
4402 !   a_ht(i,j+1) =a_ht(i,j+1) +a_Tmpv1 ! Remarked by Ning Pan, 2010-07-09
4403 !   a_ht(i,j) =a_ht(i,j) -a_Tmpv1     ! Remarked by Ning Pan, 2010-07-09
4404    ENDDO
4405
4406    DO k =k_end+1, 2, -1
4407    DO i =i_end, i_start, -1
4408
4409    IF(top_lid .and. k.eq.k_end+1) THEN
4410
4411    rhs(i,k) =Tmpv309(i,k)
4412
4413    a_rhs(i,k) =0.0
4414
4415    END IF
4416
4417    rhs(i,k) =Tmpv308(i,k)
4418
4419    a_Tmpv2 =a_rhs(i,k)
4420    a_rhs(i,k) =0.0
4421    a_ph(i,k,j) =a_ph(i,k,j) +a_Tmpv2
4422    a_Tmpv1 =a_Tmpv2
4423    a_rhs(i,k) =a_rhs(i,k) +msfty(i,j)*mut_inv(i)*a_Tmpv1
4424    a_mut_inv(i) =a_mut_inv(i) +msfty(i,j)*rhs(i,k)*a_Tmpv1
4425    ENDDO
4426    ENDDO
4427
4428    DO k =k_end, 2, -1
4429    DO i =i_end, i_start, -1
4430
4431    rhs(i,k) =Tmpv307(i,k)
4432
4433    a_Tmpv3 =a_rhs(i,k)
4434    a_rhs(i,k) =0.0
4435    a_rhs(i,k) =a_rhs(i,k) +a_Tmpv3
4436    a_Tmpv2 =-a_Tmpv3
4437    a_Tmpv1 =dts*a_Tmpv2
4438    a_wdwn(i,k+1) =a_wdwn(i,k+1) +fnm(k)*a_Tmpv1
4439    a_wdwn(i,k) =a_wdwn(i,k) +fnp(k)*a_Tmpv1
4440    ENDDO
4441    ENDDO
4442
4443    DO k =k_end, 1, -1
4444    DO i =i_end, i_start, -1
4445
4446    rhs(i,k+1) =Tmpv306(i,k)
4447
4448    a_Tmpv2 =a_rhs(i,k+1)
4449    a_rhs(i,k+1) =0.0
4450    a_Tmpv1 =dts*a_Tmpv2
4451    a_ph_tend(i,k+1,j) =a_ph_tend(i,k+1,j) +a_Tmpv1
4452    a_w(i,k+1,j) =a_w(i,k+1,j) +.5*g*(1.-epssm)*a_Tmpv1
4453    a_Tmpv7 =a_wdwn(i,k+1)
4454    a_wdwn(i,k+1) =0.0
4455    a_Tmpv3 =Tmpv305(i,k)*a_Tmpv7
4456    a_Tmpv6 =Tmpv304(i,k)*a_Tmpv7
4457    a_Tmpv5 =a_Tmpv6
4458    a_Tmpv4 =a_Tmpv5
4459    a_ph_1(i,k+1,j) =a_ph_1(i,k+1,j) +a_Tmpv4
4460    a_ph_1(i,k,j) =a_ph_1(i,k,j) -a_Tmpv4
4461    a_Tmpv2 =rdnw(k)*a_Tmpv3
4462    a_Tmpv1 =.5*a_Tmpv2
4463    a_ww(i,k+1,j) =a_ww(i,k+1,j) +a_Tmpv1
4464    a_ww(i,k,j) =a_ww(i,k,j) +a_Tmpv1
4465
4466    t_2ave(i,k,j) =Tmpv303(i,k)
4467
4468    a_Tmpv3 =a_t_2ave(i,k,j)
4469    a_t_2ave(i,k,j) =0.0
4470    a_Tmpv1 =a_Tmpv3/Tmpv302(i,k)
4471    a_Tmpv2 =-Tmpv301(i,k)/(Tmpv302(i,k)*Tmpv302(i,k))*a_Tmpv3
4472    a_muts(i,j) =a_muts(i,j) +(t0 +t_1(i,k,j))*a_Tmpv2
4473    a_t_1(i,k,j) =a_t_1(i,k,j) +muts(i,j)*a_Tmpv2
4474    a_t_2ave(i,k,j) =a_t_2ave(i,k,j) +a_Tmpv1
4475    a_muave(i,j) =a_muave(i,j) +t0*a_Tmpv1
4476
4477    t_2ave(i,k,j) =Tmpv300(i,k)
4478
4479    a_Tmpv2 =a_t_2ave(i,k,j)
4480    a_t_2ave(i,k,j) =0.0
4481    a_Tmpv1 =.5*a_Tmpv2
4482    a_t_2(i,k,j) =a_t_2(i,k,j) +(1.+epssm)*a_Tmpv1
4483    a_t_2ave(i,k,j) =a_t_2ave(i,k,j) +(1.-epssm)*a_Tmpv1
4484    ENDDO
4485    ENDDO
4486
4487    DO i =i_end, i_start, -1
4488
4489 !  mut_inv(i) =Tmpv200(i), REVISED BY WALLS
4490
4491    a_mut(i,j) =a_mut(i,j) -1./(mut(i,j)*mut(i,j))*a_mut_inv(i)
4492    a_mut_inv(i) =0.0
4493    ENDDO
4494
4495    ENDDO
4496
4497 !LPB[5]
4498    DO i =i_end, i_start, -1
4499
4500 !  rhs(i,1) =0.
4501
4502    a_rhs(i,1) =0.0
4503
4504    ENDDO
4505
4506 !LPB[4]
4507 !  pi =4.*Atan(1.)
4508 !  dampmag =dts*config_flags%dampcoef
4509 !  hdepth =config_flags%zdamp
4510
4511 !LPB[3]
4512
4513 !  IF( config_flags%specified .or. config_flags%nested ) THEN
4514 !  j_start =max(jts, jds+1)
4515 !  j_end =min(jte, jde-2)
4516 !  ENDIF
4517
4518    IF( config_flags%specified .or. config_flags%nested ) THEN
4519
4520    ENDIF
4521
4522 !LPB[2]
4523
4524 !LPB[1]
4525
4526 !  IF( .NOT. config_flags%periodic_x ) THEN
4527 !  IF( config_flags%specified .or. config_flags%nested ) THEN
4528 !  i_start =max(its, ids+1)
4529 !  i_end =min(ite, ide-2)
4530 !  ENDIF
4531 !  ENDIF
4532
4533    IF( .NOT. config_flags%periodic_x ) THEN
4534
4535    IF( config_flags%specified .or. config_flags%nested ) THEN
4536
4537    ENDIF
4538
4539    ENDIF
4540
4541 !LPB[0]
4542 !  i_start =its
4543 !  i_end =min(ite, ide-1)
4544 !  j_start =jts
4545 !  j_end =min(jte, jde-1)
4546 !  k_start =kts
4547 !  k_end =kte-1
4548
4549    END SUBROUTINE a_advance_w
4550
4551    SUBROUTINE a_sumflux(ru,a_ru,rv,a_rv,ww,a_ww,u_lin,a_u_lin,v_lin, &
4552    a_v_lin,ww_lin,a_ww_lin,muu,a_muu,muv,a_muv,ru_m,a_ru_m,rv_m,a_rv_m, &
4553    ww_m,a_ww_m,epssm,msfux,msfuy,msfvx,msfvx_inv,msfvy,iteration,number_of_small_timesteps, &
4554    ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte)
4555
4556 !PART I: DECLARATION OF VARIABLES
4557
4558    IMPLICIT NONE
4559
4560    INTEGER :: K0_ADJ,K1_ADJ,K2_ADJ,K3_ADJ
4561    INTEGER :: number_of_small_timesteps
4562    INTEGER :: iteration
4563    INTEGER :: ids,ide,jds,jde,kds,kde
4564    INTEGER :: ims,ime,jms,jme,kms,kme
4565    INTEGER :: its,ite,jts,jte,kts,kte
4566    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_ru,ru,a_rv,rv,a_ww,ww,a_u_lin, &
4567    u_lin,a_v_lin,v_lin,a_ww_lin,ww_lin
4568    REAL,DIMENSION(ims:ime,kms:kme,jms:jme) :: a_ru_m,ru_m,a_rv_m,rv_m,a_ww_m,ww_m
4569    REAL,DIMENSION(ims:ime,jms:jme) :: a_muu,muu,a_muv,muv,msfux,msfuy,msfvx,msfvy,msfvx_inv
4570    INTEGER :: mini,minj,mink
4571    REAL :: epssm
4572    INTEGER :: i,j,k
4573
4574    REAL :: a_Tmpv1,Tmpv001,a_Tmpv2,Tmpv002,a_Tmpv3,Tmpv003
4575
4576 !PART II: CALCULATIONS OF B. S. TRAJECTORY
4577
4578 !LPB[0]
4579
4580 !LPB[1]
4581     IF (iteration == 1 )THEN
4582
4583          DO  j = jts, jte
4584          DO  k = kts, kte
4585          DO  i = its, ite
4586            ru_m(i,k,j)  = 0.
4587            rv_m(i,k,j)  = 0.
4588            ww_m(i,k,j)  = 0.
4589          ENDDO
4590          ENDDO
4591          ENDDO
4592
4593    ENDIF
4594
4595 !LPB[2]
4596
4597      mini = min(ide-1,ite)
4598      minj = min(jde-1,jte)
4599      mink = min(kde-1,kte)
4600
4601 !LPB[3]
4602        DO  j = jts, minj
4603
4604        DO  k = kts, mink
4605        DO  i = its, mini
4606          ru_m(i,k,j)  = ru_m(i,k,j) + ru(i,k,j)
4607          rv_m(i,k,j)  = rv_m(i,k,j) + rv(i,k,j)
4608          ww_m(i,k,j)  = ww_m(i,k,j) + ww(i,k,j)
4609        ENDDO
4610        ENDDO
4611
4612        ENDDO
4613
4614 !LPB[4]
4615
4616 !LPB[5]
4617     IF (ite .GT. mini) THEN
4618
4619          DO  j = jts, minj
4620          DO  k = kts, mink
4621          DO  i = mini+1, ite
4622            ru_m(i,k,j)  = ru_m(i,k,j) + ru(i,k,j)
4623          ENDDO
4624          ENDDO
4625          ENDDO
4626
4627    END IF
4628
4629 !LPB[6]
4630
4631 !LPB[7]
4632
4633     IF (jte .GT. minj) THEN
4634
4635          DO  j = minj+1, jte
4636          DO  k = kts, mink
4637          DO  i = its, mini
4638            rv_m(i,k,j)  = rv_m(i,k,j) + rv(i,k,j)
4639          ENDDO
4640          ENDDO
4641          ENDDO
4642
4643    END IF
4644
4645 !LPB[8]
4646
4647 !LPB[9]
4648
4649     IF ( kte .GT. mink) THEN
4650
4651          DO  j = jts, minj
4652          DO  k = mink+1, kte
4653          DO  i = its, mini
4654            ww_m(i,k,j)  = ww_m(i,k,j) + ww(i,k,j)
4655          ENDDO
4656          ENDDO
4657          ENDDO
4658
4659    END IF
4660
4661 !LPB[10]
4662
4663 !!LPB[11]
4664 !
4665 !  IF (iteration == number_of_small_timesteps) THEN
4666
4667 !       DO  j = jts, minj
4668 !       DO  k = kts, mink
4669 !       DO  i = its, mini
4670 !         ru_m(i,k,j)  = ru_m(i,k,j) / number_of_small_timesteps     &
4671 !                        + muu(i,j)*u_lin(i,k,j)/msfuy(i,j)
4672 !         rv_m(i,k,j)  = rv_m(i,k,j) / number_of_small_timesteps     &
4673 !                        + muv(i,j)*v_lin(i,k,j)*msfvx_inv(i,j)
4674 !         ww_m(i,k,j)  = ww_m(i,k,j) / number_of_small_timesteps     &
4675 !                        + ww_lin(i,k,j)
4676 !       ENDDO
4677 !       ENDDO
4678 !       ENDDO
4679 !    IF (ite .GT. mini) THEN
4680
4681 !         DO  j = jts, minj
4682 !         DO  k = kts, mink
4683 !         DO  i = mini+1, ite
4684 !           ru_m(i,k,j)  = ru_m(i,k,j) / number_of_small_timesteps     &
4685 !                        + muu(i,j)*u_lin(i,k,j)/msfuy(i,j)
4686 !         ENDDO
4687 !         ENDDO
4688 !         ENDDO
4689 !       END IF
4690 !    IF (jte .GT. minj) THEN
4691
4692 !         DO  j = minj+1, jte
4693 !         DO  k = kts, mink
4694 !         DO  i = its, mini
4695 !           rv_m(i,k,j)  = rv_m(i,k,j) / number_of_small_timesteps     &
4696 !                        + muv(i,j)*v_lin(i,k,j)*msfvx_inv(i,j)
4697 !         ENDDO
4698 !         ENDDO
4699 !         ENDDO
4700 !       END IF
4701 !    IF ( kte .GT. mink) THEN
4702
4703 !         DO  j = jts, minj
4704 !         DO  k = mink+1, kte
4705 !         DO  i = its, mini
4706 !           ww_m(i,k,j)  = ww_m(i,k,j) / number_of_small_timesteps     &
4707 !                        + ww_lin(i,k,j)
4708 !         ENDDO
4709 !         ENDDO
4710 !         ENDDO
4711 !       END IF
4712
4713 !   ENDIF
4714
4715 !PART IV: REVERSE/BACKWARD ACCUMULATIONS
4716
4717 !LPB[11]
4718
4719    IF(iteration == number_of_small_timesteps) THEN
4720    DO j =jts, minj
4721    DO k =kts, mink
4722    DO i =its, mini
4723    Tmpv001 =muu(i,j)*u_lin(i,k,j)
4724    Tmpv002 =Tmpv001/msfuy(i,j)
4725    Tmpv003 =ru_m(i,k,j)/number_of_small_timesteps +Tmpv002
4726 !  ru_m(i,k,j) =Tmpv003
4727
4728    Tmpv001 =muv(i,j)*v_lin(i,k,j)
4729    Tmpv002 =Tmpv001*msfvx_inv(i,j)
4730    Tmpv003 =rv_m(i,k,j)/number_of_small_timesteps +Tmpv002
4731 !  rv_m(i,k,j) =Tmpv003
4732
4733    Tmpv001 =ww_m(i,k,j)/number_of_small_timesteps +ww_lin(i,k,j)
4734 !  ww_m(i,k,j) =Tmpv001
4735
4736    ENDDO
4737    ENDDO
4738    ENDDO
4739    IF(ite .GT. mini) THEN
4740    DO j =jts, minj
4741    DO k =kts, mink
4742    DO i =mini+1, ite
4743    Tmpv001 =muu(i,j)*u_lin(i,k,j)
4744    Tmpv002 =Tmpv001/msfuy(i,j)
4745    Tmpv003 =ru_m(i,k,j)/number_of_small_timesteps +Tmpv002
4746 !  ru_m(i,k,j) =Tmpv003
4747
4748    ENDDO
4749    ENDDO
4750    ENDDO
4751    END IF
4752    IF(jte .GT. minj) THEN
4753    DO j =minj+1, jte
4754    DO k =kts, mink
4755    DO i =its, mini
4756    Tmpv001 =muv(i,j)*v_lin(i,k,j)
4757    Tmpv002 =Tmpv001*msfvx_inv(i,j)
4758    Tmpv003 =rv_m(i,k,j)/number_of_small_timesteps +Tmpv002
4759 !  rv_m(i,k,j) =Tmpv003
4760
4761    ENDDO
4762    ENDDO
4763    ENDDO
4764    END IF
4765    IF( kte .GT. mink) THEN
4766    DO j =jts, minj
4767    DO k =mink+1, kte
4768    DO i =its, mini
4769    Tmpv001 =ww_m(i,k,j)/number_of_small_timesteps +ww_lin(i,k,j)
4770 !  ww_m(i,k,j) =Tmpv001
4771
4772    ENDDO
4773    ENDDO
4774    ENDDO
4775    END IF
4776    ENDIF
4777
4778    IF(iteration == number_of_small_timesteps) THEN
4779
4780    IF( kte .GT. mink) THEN
4781
4782    DO j =minj, jts, -1
4783    DO k =kte, mink+1, -1
4784    DO i =mini, its, -1
4785    a_Tmpv1 =a_ww_m(i,k,j)
4786    a_ww_m(i,k,j) =0.0
4787    a_ww_m(i,k,j) =a_ww_m(i,k,j) +1.0/number_of_small_timesteps*a_Tmpv1
4788    a_ww_lin(i,k,j) =a_ww_lin(i,k,j) +a_Tmpv1
4789    ENDDO
4790    ENDDO
4791    ENDDO
4792
4793    END IF
4794
4795    IF(jte .GT. minj) THEN
4796
4797    DO j =jte, minj+1, -1
4798    DO k =mink, kts, -1
4799    DO i =mini, its, -1
4800    a_Tmpv3 =a_rv_m(i,k,j)
4801    a_rv_m(i,k,j) =0.0
4802    a_rv_m(i,k,j) =a_rv_m(i,k,j) +1.0/number_of_small_timesteps*a_Tmpv3
4803    a_Tmpv2 =a_Tmpv3
4804    a_Tmpv1 =msfvx_inv(i,j)*a_Tmpv2
4805    a_muv(i,j) =a_muv(i,j) +v_lin(i,k,j)*a_Tmpv1
4806    a_v_lin(i,k,j) =a_v_lin(i,k,j) +muv(i,j)*a_Tmpv1
4807    ENDDO
4808    ENDDO
4809    ENDDO
4810
4811    END IF
4812
4813    IF(ite .GT. mini) THEN
4814
4815    DO j =minj, jts, -1
4816    DO k =mink, kts, -1
4817    DO i =ite, mini+1, -1
4818    a_Tmpv3 =a_ru_m(i,k,j)
4819    a_ru_m(i,k,j) =0.0
4820    a_ru_m(i,k,j) =a_ru_m(i,k,j) +1.0/number_of_small_timesteps*a_Tmpv3
4821    a_Tmpv2 =a_Tmpv3
4822    a_Tmpv1 =a_Tmpv2/msfuy(i,j)
4823    a_muu(i,j) =a_muu(i,j) +u_lin(i,k,j)*a_Tmpv1
4824    a_u_lin(i,k,j) =a_u_lin(i,k,j) +muu(i,j)*a_Tmpv1
4825    ENDDO
4826    ENDDO
4827    ENDDO
4828
4829    END IF
4830    DO j =minj, jts, -1
4831    DO k =mink, kts, -1
4832    DO i =mini, its, -1
4833    a_Tmpv1 =a_ww_m(i,k,j)
4834    a_ww_m(i,k,j) =0.0
4835    a_ww_m(i,k,j) =a_ww_m(i,k,j) +1.0/number_of_small_timesteps*a_Tmpv1
4836    a_ww_lin(i,k,j) =a_ww_lin(i,k,j) +a_Tmpv1
4837    a_Tmpv3 =a_rv_m(i,k,j)
4838    a_rv_m(i,k,j) =0.0
4839    a_rv_m(i,k,j) =a_rv_m(i,k,j) +1.0/number_of_small_timesteps*a_Tmpv3
4840    a_Tmpv2 =a_Tmpv3
4841    a_Tmpv1 =msfvx_inv(i,j)*a_Tmpv2
4842    a_muv(i,j) =a_muv(i,j) +v_lin(i,k,j)*a_Tmpv1
4843    a_v_lin(i,k,j) =a_v_lin(i,k,j) +muv(i,j)*a_Tmpv1
4844    a_Tmpv3 =a_ru_m(i,k,j)
4845    a_ru_m(i,k,j) =0.0
4846    a_ru_m(i,k,j) =a_ru_m(i,k,j) +1.0/number_of_small_timesteps*a_Tmpv3
4847    a_Tmpv2 =a_Tmpv3
4848    a_Tmpv1 =a_Tmpv2/msfuy(i,j)
4849    a_muu(i,j) =a_muu(i,j) +u_lin(i,k,j)*a_Tmpv1
4850    a_u_lin(i,k,j) =a_u_lin(i,k,j) +muu(i,j)*a_Tmpv1
4851    ENDDO
4852    ENDDO
4853    ENDDO
4854
4855    ENDIF
4856
4857 !LPB[10]
4858
4859 !LPB[9]
4860
4861 !  IF( kte .GT. mink) THEN
4862 !  DO j =jts, minj
4863 !  DO k =mink+1, kte
4864 !  DO i =its, mini
4865 !  Tmpv001 =ww_m(i,k,j) +ww(i,k,j)
4866 !  ww_m(i,k,j) =Tmpv001
4867
4868 !  ENDDO
4869 !  ENDDO
4870 !  ENDDO
4871 !  END IF
4872
4873    IF( kte .GT. mink) THEN
4874
4875    DO j =minj, jts, -1
4876    DO k =kte, mink+1, -1
4877    DO i =mini, its, -1
4878    a_Tmpv1 =a_ww_m(i,k,j)
4879    a_ww_m(i,k,j) =0.0
4880    a_ww_m(i,k,j) =a_ww_m(i,k,j) +a_Tmpv1
4881    a_ww(i,k,j) =a_ww(i,k,j) +a_Tmpv1
4882    ENDDO
4883    ENDDO
4884    ENDDO
4885
4886    END IF
4887
4888 !LPB[8]
4889
4890 !LPB[7]
4891
4892 !  IF(jte .GT. minj) THEN
4893 !  DO j =minj+1, jte
4894 !  DO k =kts, mink
4895 !  DO i =its, mini
4896 !  Tmpv001 =rv_m(i,k,j) +rv(i,k,j)
4897 !  rv_m(i,k,j) =Tmpv001
4898
4899 !  ENDDO
4900 !  ENDDO
4901 !  ENDDO
4902 !  END IF
4903
4904    IF(jte .GT. minj) THEN
4905
4906    DO j =jte, minj+1, -1
4907    DO k =mink, kts, -1
4908    DO i =mini, its, -1
4909    a_Tmpv1 =a_rv_m(i,k,j)
4910    a_rv_m(i,k,j) =0.0
4911    a_rv_m(i,k,j) =a_rv_m(i,k,j) +a_Tmpv1
4912    a_rv(i,k,j) =a_rv(i,k,j) +a_Tmpv1
4913    ENDDO
4914    ENDDO
4915    ENDDO
4916
4917    END IF
4918
4919 !LPB[6]
4920
4921 !LPB[5]
4922
4923 !  IF(ite .GT. mini) THEN
4924 !  DO j =jts, minj
4925 !  DO k =kts, mink
4926 !  DO i =mini+1, ite
4927 !  Tmpv001 =ru_m(i,k,j) +ru(i,k,j)
4928 !  ru_m(i,k,j) =Tmpv001
4929
4930 !  ENDDO
4931 !  ENDDO
4932 !  ENDDO
4933 !  END IF
4934
4935    IF(ite .GT. mini) THEN
4936
4937    DO j =minj, jts, -1
4938    DO k =mink, kts, -1
4939    DO i =ite, mini+1, -1
4940    a_Tmpv1 =a_ru_m(i,k,j)
4941    a_ru_m(i,k,j) =0.0
4942    a_ru_m(i,k,j) =a_ru_m(i,k,j) +a_Tmpv1
4943    a_ru(i,k,j) =a_ru(i,k,j) +a_Tmpv1
4944    ENDDO
4945    ENDDO
4946    ENDDO
4947
4948    END IF
4949
4950 !LPB[4]
4951
4952 !LPB[3]
4953    DO j =minj, jts, -1
4954
4955 !  DO k =kts, mink
4956 !  DO i =its, mini
4957 !  Tmpv001 =ru_m(i,k,j) +ru(i,k,j)
4958 !  ru_m(i,k,j) =Tmpv001
4959
4960 !  Tmpv001 =rv_m(i,k,j) +rv(i,k,j)
4961 !  rv_m(i,k,j) =Tmpv001
4962
4963 !  Tmpv001 =ww_m(i,k,j) +ww(i,k,j)
4964 !  ww_m(i,k,j) =Tmpv001
4965
4966 !  ENDDO
4967 !  ENDDO
4968
4969    DO k =mink, kts, -1
4970    DO i =mini, its, -1
4971    a_Tmpv1 =a_ww_m(i,k,j)
4972    a_ww_m(i,k,j) =0.0
4973    a_ww_m(i,k,j) =a_ww_m(i,k,j) +a_Tmpv1
4974    a_ww(i,k,j) =a_ww(i,k,j) +a_Tmpv1
4975    a_Tmpv1 =a_rv_m(i,k,j)
4976    a_rv_m(i,k,j) =0.0
4977    a_rv_m(i,k,j) =a_rv_m(i,k,j) +a_Tmpv1
4978    a_rv(i,k,j) =a_rv(i,k,j) +a_Tmpv1
4979    a_Tmpv1 =a_ru_m(i,k,j)
4980    a_ru_m(i,k,j) =0.0
4981    a_ru_m(i,k,j) =a_ru_m(i,k,j) +a_Tmpv1
4982    a_ru(i,k,j) =a_ru(i,k,j) +a_Tmpv1
4983    ENDDO
4984    ENDDO
4985
4986    ENDDO
4987
4988 !LPB[2]
4989 !  mini =min(ide-1, ite)
4990 !  minj =min(jde-1, jte)
4991 !  mink =min(kde-1, kte)
4992
4993 !LPB[1]
4994
4995 !  IF(iteration == 1 ) THEN
4996 !  DO j =jts, jte
4997 !  DO k =kts, kte
4998 !  DO i =its, ite
4999 !  ru_m(i,k,j) =0.
5000
5001 !  rv_m(i,k,j) =0.
5002
5003 !  ww_m(i,k,j) =0.
5004
5005 !  ENDDO
5006 !  ENDDO
5007 !  ENDDO
5008 !  ENDIF
5009
5010    IF(iteration == 1 ) THEN
5011
5012    DO j =jte, jts, -1
5013    DO k =kte, kts, -1
5014    DO i =ite, its, -1
5015    a_ww_m(i,k,j) =0.0
5016    a_rv_m(i,k,j) =0.0
5017    a_ru_m(i,k,j) =0.0
5018    ENDDO
5019    ENDDO
5020    ENDDO
5021
5022    ENDIF
5023
5024 !LPB[0]
5025
5026    END SUBROUTINE a_sumflux
5027
5028    SUBROUTINE a_init_module_small_step
5029
5030    END SUBROUTINE a_init_module_small_step
5031
5032 END MODULE a_module_small_step_em
5033