wrftladj/module_small_step_em_tl.F

   1
   2 ! ======================================================================================
   3 ! This file was generated by the version 5.3.5 of DFT on 06/19/2010. The differentiation
   4 ! transforming system(DFT) was jointly developed and sponsored by LASG of IAP(1998-2010)
   5 ! and LSEC of ICMSEC, AMSS(2001-2003)
   6 ! The copyright of the DFT system was declared by Walls at LASG, 1998-2010
   7 ! ======================================================================================
   8 ! Ning Pan, 2010-07-10: Change Diff_ to g_
   9 ! Ning Pan, 2010-07-11: Change the order of arguments
  10
  11  MODULE g_module_small_step_em
  12
  13  USE module_configure
  14  USE module_model_constants
  15
  16  CONTAINS
  17
  18 SUBROUTINE g_advance_all (u, g_u, ru_tend, g_ru_tend, v, g_v, rv_tend, g_rv_tend,     &
  19                           w, g_w, rw_tend, g_rw_tend, t, g_t, t_tend, g_t_tend,       &
  20                           mu, g_mu, mu_tend, g_mu_tend, ph, g_ph, ph_tend, g_ph_tend, &
  21                           muu, g_muu, muv, g_muv, mut, g_mut,                         &
  22                           msfuy, msfvx, msfty,           &
  23                           dts,                           &
  24                           config_flags, spec_zone,       &
  25                           ids, ide, jds, jde, kds, kde,  &
  26                           ims, ime, jms, jme, kms, kme,  &
  27                           ips, ipe, jps, jpe, kps, kpe,  &
  28                           its, ite, jts, jte, kts, kte  )
  29
  30       IMPLICIT NONE  ! religion first
  31
  32       TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
  33
  34       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
  35       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
  36       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
  37       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
  38       INTEGER,      INTENT(IN   )    :: spec_zone
  39
  40       REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),  &
  41             INTENT(INOUT) ::                          &
  42                                                   g_u,  &
  43                                                   g_v,  &
  44                                                   g_w,  &
  45                                                   g_t,  &
  46                                                   g_ph
  47       REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),  &
  48             INTENT(INOUT) ::                          &
  49                                                   u,  &
  50                                                   v,  &
  51                                                   w,  &
  52                                                   t,  &
  53                                                   ph
  54
  55       REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
  56             INTENT(IN   ) ::                          &
  57                                              g_ru_tend, &
  58                                              g_rv_tend, &
  59                                              g_rw_tend, &
  60                                               g_t_tend, &
  61                                              g_ph_tend
  62       REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
  63             INTENT(IN   ) ::                          &
  64                                              ru_tend, &
  65                                              rv_tend, &
  66                                              rw_tend, &
  67                                               t_tend, &
  68                                              ph_tend
  69
  70       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(INOUT) :: g_mu
  71       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(INOUT) :: mu
  72       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: g_mu_tend
  73       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: mu_tend
  74
  75       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: g_muu,  &
  76                                                                 g_muv,  &
  77                                                                 g_mut
  78       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: muu,  &
  79                                                                 muv,  &
  80                                                                 mut
  81
  82       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: msfuy,  &
  83                                                                 msfvx,  &
  84                                                                 msfty
  85
  86       REAL,                                    INTENT(IN   ) :: dts
  87
  88
  89       INTEGER :: i,j,k, i_start, i_end, j_start, j_end, k_start, k_end
  90       INTEGER :: i_endu, j_endv, k_endw
  91
  92       INTEGER :: i_start_u_tend, i_end_u_tend, j_start_v_tend, j_end_v_tend
  93
  94 !<DESCRIPTION>
  95 !
  96 !  advance_all advances the explicit perturbation horizontal momentum
  97 !  equations (u,v) by adding in the large-timestep tendency along with
  98 !  the small timestep pressure gradient tendency.
  99 !
 100 !</DESCRIPTION>
 101
 102 !  now, the real work.
 103 !  set the loop bounds taking into account boundary conditions.
 104
 105    IF( config_flags%nested .or. config_flags%specified ) THEN
 106       i_start = max( its,ids+spec_zone )
 107       i_end   = min( ite,ide-spec_zone-1 )
 108       j_start = max( jts,jds+spec_zone )
 109       j_end   = min( jte,jde-spec_zone-1 )
 110       k_start = kts
 111       k_end   = min( kte, kde-1 )
 112
 113       i_endu = min( ite,ide-spec_zone )
 114       j_endv = min( jte,jde-spec_zone )
 115       k_endw = k_end
 116
 117       IF( config_flags%periodic_x) THEN
 118         i_start = its
 119         i_end   = min(ite,ide-1)
 120         i_endu = ite
 121       ENDIF
 122    ELSE
 123       i_start = its
 124       i_end   = min(ite,ide-1)
 125       j_start = jts
 126       j_end   = min(jte,jde-1)
 127       k_start = kts
 128       k_end   = kte-1
 129
 130       i_endu = ite
 131       j_endv = jte
 132       k_endw = k_end
 133    ENDIF
 134
 135    i_start_u_tend = i_start
 136    i_end_u_tend   = i_endu
 137    j_start_v_tend = j_start
 138    j_end_v_tend   = j_endv
 139
 140    IF ( config_flags%symmetric_xs .and. (its == ids) ) &
 141                   i_start_u_tend = i_start_u_tend+1
 142    IF ( config_flags%symmetric_xe .and. (ite == ide) ) &
 143                   i_end_u_tend = i_end_u_tend-1
 144    IF ( config_flags%symmetric_ys .and. (jts == jds) ) &
 145                   j_start_v_tend = j_start_v_tend+1
 146    IF ( config_flags%symmetric_ye .and. (jte == jde) ) &
 147                   j_end_v_tend = j_end_v_tend-1
 148
 149    u_outer_j_loop: DO j = j_start, j_end
 150       DO k = k_start, k_end
 151       DO i = i_start_u_tend, i_end_u_tend
 152         g_u(i,k,j) = g_u(i,k,j) &
 153                    + dts*msfuy(i,j)*( g_ru_tend(i,k,j)/muu(i,j) - &
 154                                       ru_tend(i,k,j)*g_muu(i,j)/(muu(i,j)*muu(i,j)) )
 155         u(i,k,j) = u(i,k,j) + dts*ru_tend(i,k,j)*msfuy(i,j)/muu(i,j)
 156       ENDDO
 157       ENDDO
 158    ENDDO u_outer_j_loop
 159
 160    v_outer_j_loop: DO j = j_start_v_tend, j_end_v_tend
 161       DO k = k_start, k_end
 162       DO i = i_start, i_end
 163         g_v(i,k,j) = g_v(i,k,j) &
 164                    + dts*msfvx(i,j)*( g_rv_tend(i,k,j)/muv(i,j) - &
 165                                       rv_tend(i,k,j)*g_muv(i,j)/(muv(i,j)*muv(i,j)) )
 166         v(i,k,j) = v(i,k,j) + dts*rv_tend(i,k,j)*msfvx(i,j)/muv(i,j)
 167       ENDDO
 168       ENDDO
 169    ENDDO v_outer_j_loop
 170
 171    i_start = its
 172    i_end   = min(ite,ide-1)
 173    j_start = jts
 174    j_end   = min(jte,jde-1)
 175    k_start = kts
 176    k_end   = kte-1
 177    IF ( .NOT. config_flags%periodic_x )THEN
 178      IF ( config_flags%specified .or. config_flags%nested ) then
 179        i_start = max(its,ids+1)
 180        i_end   = min(ite,ide-2)
 181      ENDIF
 182    ENDIF
 183    IF ( config_flags%specified .or. config_flags%nested ) then
 184      j_start = max(jts,jds+1)
 185      j_end   = min(jte,jde-2)
 186    ENDIF
 187
 188    DO j = j_start, j_end
 189    DO i=i_start, i_end
 190      g_mu(i,j) = g_mu(i,j)+dts*g_mu_tend(i,j)
 191      mu(i,j) = mu(i,j)+dts*mu_tend(i,j)
 192    ENDDO
 193    ENDDO
 194
 195    DO j=j_start, j_end
 196      DO k=1,k_end
 197      DO i=i_start, i_end
 198        g_t(i,k,j) = g_t(i,k,j) &
 199                   + dts*msfty(i,j)*( g_t_tend(i,k,j)/mut(i,j) - &
 200                                      t_tend(i,k,j)*g_mut(i,j)/(mut(i,j)*mut(i,j)) )
 201        t(i,k,j) = t(i,k,j) + msfty(i,j)*dts*t_tend(i,k,j)/mut(i,j)
 202      END DO
 203      END DO
 204    ENDDO
 205
 206    i_start = its
 207    i_end   = min(ite,ide-1)
 208    j_start = jts
 209    j_end   = min(jte,jde-1)
 210    k_start = kts
 211    k_end   = kte-1
 212    IF ( .NOT. config_flags%periodic_x )THEN
 213       IF ( config_flags%specified .or. config_flags%nested ) then
 214         i_start = max(its,ids+1)
 215         i_end   = min(ite,ide-2)
 216       ENDIF
 217    ENDIF
 218    IF ( config_flags%specified .or. config_flags%nested ) then
 219       j_start = max(jts,jds+1)
 220       j_end   = min(jte,jde-2)
 221    ENDIF
 222
 223 IF ( config_flags%non_hydrostatic ) THEN
 224   j_loop_w:  DO j = j_start, j_end
 225     DO k=2,k_end+1
 226       DO i=i_start, i_end
 227         g_w(i,k,j) = g_w(i,k,j) &
 228                    + dts*msfty(i,j)*( g_rw_tend(i,k,j)/mut(i,j) - &
 229                                       rw_tend(i,k,j)*g_mut(i,j)/(mut(i,j)*mut(i,j)) )
 230         w(i,k,j)=w(i,k,j)+dts*rw_tend(i,k,j)*msfty(i,j)/mut(i,j)
 231
 232         g_ph(i,k,j) = g_ph(i,k,j) &
 233                     + dts*msfty(i,j)*( g_ph_tend(i,k,j)/mut(i,j) - &
 234                                        ph_tend(i,k,j)*g_mut(i,j)/(mut(i,j)*mut(i,j)) )
 235         ph(i,k,j) = ph(i,k,j)+dts*ph_tend(i,k,j)*msfty(i,j)/mut(i,j)
 236       ENDDO
 237     ENDDO
 238   ENDDO j_loop_w
 239 ENDIF
 240
 241 END SUBROUTINE g_advance_all
 242
 243
 244 SUBROUTINE g_save_ph_mu ( ph_1, g_ph_1, ph_2, g_ph_2, ph_save, g_ph_save, &
 245                           mu_1, g_mu_1, mu_2, g_mu_2, mu_save, g_mu_save, &
 246                           rk_step,                      &
 247                           ids,ide, jds,jde, kds,kde,    &
 248                           ims,ime, jms,jme, kms,kme,    &
 249                           its,ite, jts,jte, kts,kte    )
 250
 251   IMPLICIT NONE  ! religion first
 252
 253 ! declarations for the stuff coming in
 254
 255   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
 256   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
 257   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
 258
 259   INTEGER,      INTENT(IN   )    :: rk_step
 260
 261   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: g_ph_1
 262   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: ph_1
 263
 264   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(  OUT) :: g_ph_save
 265   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(  OUT) :: ph_save
 266
 267   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: g_ph_2
 268   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: ph_2
 269
 270   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: g_mu_1,g_mu_2
 271   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mu_1,mu_2
 272   REAL, DIMENSION(ims:ime, jms:jme), INTENT(  OUT) :: g_mu_save
 273   REAL, DIMENSION(ims:ime, jms:jme), INTENT(  OUT) :: mu_save
 274
 275 ! local variables
 276
 277   INTEGER :: i, j, k
 278   INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
 279
 280
 281     i_start = its
 282     i_end   = min(ite,ide-1)
 283     j_start = jts
 284     j_end   = min(jte,jde-1)
 285     k_start = kts
 286     k_end = min(kte,kde-1)
 287
 288     !  if this is the first RK step, reset *_1 to *_2
 289     !  (we are replacing the t-dt fields with the time t fields)
 290
 291     IF ((rk_step == 1) ) THEN
 292
 293       DO j=j_start, j_end
 294       DO k=k_start, min(kde,kte)
 295       DO i=i_start, i_end
 296         g_ph_1(i,k,j) = g_ph_2(i,k,j)
 297         ph_1(i,k,j) = ph_2(i,k,j)
 298       ENDDO
 299       ENDDO
 300       ENDDO
 301
 302       DO j=j_start, j_end
 303       DO i=i_start, i_end
 304         g_mu_1(i,j) = g_mu_2(i,j)
 305         mu_1(i,j) = mu_2(i,j)
 306       ENDDO
 307       ENDDO
 308
 309     END IF
 310
 311     ! set up the small timestep variables
 312
 313     DO j=j_start, j_end
 314 !    DO k=k_start, min(kde,kte)
 315     DO k=k_start, kde
 316     DO i=i_start, i_end
 317       g_ph_save(i,k,j) = g_ph_2(i,k,j)
 318       ph_save(i,k,j) = ph_2(i,k,j)
 319       g_ph_2(i,k,j) = g_ph_1(i,k,j)-g_ph_2(i,k,j)
 320       ph_2(i,k,j) = ph_1(i,k,j)-ph_2(i,k,j)
 321     ENDDO
 322     ENDDO
 323     ENDDO
 324
 325     DO j=j_start, j_end
 326     DO i=i_start, i_end
 327       g_mu_save(i,j)=g_mu_2(i,j)
 328       mu_save(i,j)=mu_2(i,j)
 329       g_mu_2(i,j)=g_mu_1(i,j)-g_mu_2(i,j)
 330       mu_2(i,j)=mu_1(i,j)-mu_2(i,j)
 331     ENDDO
 332     ENDDO
 333
 334 END SUBROUTINE g_save_ph_mu
 335
 336 !----------------------------------------------------------------------
 337
 338 SUBROUTINE g_restore_ph_mu ( ph_2, g_ph_2, ph_save, g_ph_save, &
 339                              mu_2, g_mu_2, mu_save, g_mu_save, &
 340                              ids,ide, jds,jde, kds,kde,    &
 341                              ims,ime, jms,jme, kms,kme,    &
 342                              its,ite, jts,jte, kts,kte )
 343
 344   IMPLICIT NONE  ! religion first
 345
 346 ! declarations for the stuff coming in
 347
 348   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
 349   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
 350   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
 351
 352   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(IN   ) :: g_ph_save
 353   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(IN   ) :: ph_save
 354   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: g_ph_2
 355   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: ph_2
 356   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: g_mu_2
 357   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mu_2
 358   REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN   ) :: g_mu_save
 359   REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN   ) :: mu_save
 360
 361 ! local variables
 362
 363   INTEGER :: i, j, k
 364   INTEGER :: i_start, i_end, j_start, j_end
 365
 366
 367   i_start = its
 368   i_end   = min(ite,ide-1)
 369   j_start = jts
 370   j_end   = min(jte,jde-1)
 371
 372   DO j = j_start, j_end
 373   DO k = kds, kde
 374   DO i = i_start, i_end
 375     g_ph_2(i,k,j) = g_ph_2(i,k,j) + g_ph_save(i,k,j)
 376     ph_2(i,k,j) = ph_2(i,k,j) + ph_save(i,k,j)
 377   ENDDO
 378   ENDDO
 379   ENDDO
 380
 381   DO j = j_start, j_end
 382   DO i = i_start, i_end
 383     g_mu_2(i,j) = g_mu_2(i,j) + g_mu_save(i,j)
 384     mu_2(i,j) = mu_2(i,j) + mu_save(i,j)
 385   ENDDO
 386   ENDDO
 387
 388 END SUBROUTINE g_restore_ph_mu
 389
 390  SUBROUTINE g_small_step_prep(u_1,g_u_1,u_2,g_u_2,v_1,g_v_1,v_2,g_v_2, &
 391  w_1,g_w_1,w_2,g_w_2,t_1,g_t_1,t_2,g_t_2,ph_1,g_ph_1,ph_2,g_ph_2, &
 392  mub,mu_1,g_mu_1,mu_2,g_mu_2,muu,g_muu,muus,g_muus,muv,g_muv,muvs, &
 393  g_muvs,mut,g_mut,muts,g_muts,mudf,g_mudf,u_save,g_u_save,v_save, &
 394  g_v_save,w_save,g_w_save,t_save,g_t_save,ph_save,g_ph_save,mu_save, &
 395  g_mu_save,ww,g_ww,ww_save,g_ww_save,c2a,g_c2a,pb,p,g_p,alt,g_alt, &
 396  msfux,msfuy,msfvx,msfvx_inv,msfvy,msftx,msfty,rdx,rdy,rk_step,ids,ide,jds,jde,kds, &
 397  kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte)
 398
 399  IMPLICIT NONE
 400  INTEGER,INTENT(IN) :: ids,ide,jds,jde,kds,kde
 401  INTEGER,INTENT(IN) :: ims,ime,jms,jme,kms,kme
 402  INTEGER,INTENT(IN) :: its,ite,jts,jte,kts,kte
 403  INTEGER,INTENT(IN) :: rk_step
 404  REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(INOUT) :: g_u_1,u_1,g_v_1,v_1, &
 405  g_w_1,w_1,g_t_1,t_1,g_ph_1,ph_1
 406  REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(OUT) :: g_u_save,u_save, &
 407  g_v_save,v_save,g_w_save,w_save,g_t_save,t_save,g_ph_save,ph_save
 408  REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(INOUT) :: g_u_2,u_2,g_v_2,v_2, &
 409  g_w_2,w_2,g_t_2,t_2,g_ph_2,ph_2
 410  REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(OUT) :: g_c2a,c2a,g_ww_save,ww_save
 411  REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(IN) :: pb,g_p,p,g_alt,alt,g_ww,ww
 412  REAL,DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: g_mu_1,mu_1,g_mu_2,mu_2
 413  REAL,DIMENSION(ims:ime,jms:jme),INTENT(IN) :: mub,g_muu,muu,g_muv,muv,g_mut, &
 414  mut,msfux,msfuy,msfvx,msfvx_inv,msfvy,msftx,msfty
 415  REAL,DIMENSION(ims:ime,jms:jme),INTENT(OUT) :: g_muus,muus,g_muvs,muvs, &
 416  g_muts,muts,g_mudf,mudf
 417  REAL,DIMENSION(ims:ime,jms:jme),INTENT(OUT) :: g_mu_save,mu_save
 418  REAL,INTENT(IN) :: rdx,rdy
 419  INTEGER :: i,j,k
 420  INTEGER :: i_start,i_end,j_start,j_end,k_start,k_end
 421  INTEGER :: i_endu,j_endv
 422
 423  i_start =its
 424
 425  i_end =min(ite,ide-1)
 426
 427  j_start =jts
 428
 429  j_end =min(jte,jde-1)
 430
 431  k_start =kts
 432
 433  k_end =min(kte,kde-1)
 434
 435  i_endu =ite
 436
 437  j_endv =jte
 438
 439  IF((rk_step == 1) ) THEN
 440
 441  DO j =j_start,j_end
 442  DO i =i_start,i_end
 443
 444  g_mu_1(i,j) =g_mu_2(i,j)
 445  mu_1(i,j) =mu_2(i,j)
 446
 447  g_ww_save(i,kde,j) =0.0
 448  ww_save(i,kde,j) =0.
 449
 450  g_ww_save(i,1,j) =0.0
 451  ww_save(i,1,j) =0.
 452
 453  g_mudf(i,j) =0.0
 454  mudf(i,j) =0.
 455
 456  ENDDO
 457  ENDDO
 458
 459  DO j =j_start,j_end
 460  DO k =k_start,k_end
 461  DO i =i_start,i_endu
 462
 463  g_u_1(i,k,j) =g_u_2(i,k,j)
 464  u_1(i,k,j) =u_2(i,k,j)
 465
 466  ENDDO
 467  ENDDO
 468  ENDDO
 469
 470  DO j =j_start,j_endv
 471  DO k =k_start,k_end
 472  DO i =i_start,i_end
 473
 474  g_v_1(i,k,j) =g_v_2(i,k,j)
 475  v_1(i,k,j) =v_2(i,k,j)
 476
 477  ENDDO
 478  ENDDO
 479  ENDDO
 480
 481  DO j =j_start,j_end
 482  DO k =k_start,k_end
 483  DO i =i_start,i_end
 484
 485  g_t_1(i,k,j) =g_t_2(i,k,j)
 486  t_1(i,k,j) =t_2(i,k,j)
 487
 488  ENDDO
 489  ENDDO
 490  ENDDO
 491
 492  DO j =j_start,j_end
 493  DO k =k_start,min(kde,kte)
 494  DO i =i_start,i_end
 495
 496  g_w_1(i,k,j) =g_w_2(i,k,j)
 497  w_1(i,k,j) =w_2(i,k,j)
 498
 499  g_ph_1(i,k,j) =g_ph_2(i,k,j)
 500  ph_1(i,k,j) =ph_2(i,k,j)
 501
 502  ENDDO
 503  ENDDO
 504  ENDDO
 505
 506  DO j =j_start,j_end
 507  DO i =i_start,i_end
 508
 509  g_muts(i,j) =g_mu_2(i,j)
 510  muts(i,j) =mub(i,j) +mu_2(i,j)
 511
 512  ENDDO
 513
 514  DO i =i_start,i_endu
 515
 516  g_muus(i,j) =g_muu(i,j)
 517  muus(i,j) =muu(i,j)
 518
 519  ENDDO
 520  ENDDO
 521
 522  DO j =j_start,j_endv
 523  DO i =i_start,i_end
 524
 525  g_muvs(i,j) =g_muv(i,j)
 526  muvs(i,j) =muv(i,j)
 527
 528  ENDDO
 529  ENDDO
 530
 531  DO j =j_start,j_end
 532  DO i =i_start,i_end
 533
 534  g_mu_save(i,j) =g_mu_2(i,j)
 535  mu_save(i,j) =mu_2(i,j)
 536
 537 !g_mu_2(i,j) =g_mu_2(i,j) -g_mu_2(i,j)
 538 !mu_2(i,j) =mu_2(i,j) -mu_2(i,j)
 539
 540  g_mu_2(i,j) =0.0  !REVISED BY WALLS
 541  mu_2(i,j) =0.0  !REVISED BY WALLS
 542
 543  ENDDO
 544  ENDDO
 545  ELSE
 546
 547  DO j =j_start,j_end
 548  DO i =i_start,i_end
 549
 550  g_muts(i,j) =g_mu_1(i,j)
 551  muts(i,j) =mub(i,j) +mu_1(i,j)
 552
 553  ENDDO
 554
 555  DO i =i_start,i_endu
 556
 557  g_muus(i,j) =0.5*(g_mu_1(i,j) +g_mu_1(i-1,j))
 558  muus(i,j) =0.5*(mub(i,j) +mu_1(i,j) +mub(i-1,j) +mu_1(i-1,j))
 559
 560  ENDDO
 561  ENDDO
 562
 563  DO j =j_start,j_endv
 564  DO i =i_start,i_end
 565
 566  g_muvs(i,j) =0.5*(g_mu_1(i,j) +g_mu_1(i,j-1))
 567  muvs(i,j) =0.5*(mub(i,j) +mu_1(i,j) +mub(i,j-1) +mu_1(i,j-1))
 568
 569  ENDDO
 570  ENDDO
 571
 572  DO j =j_start,j_end
 573  DO i =i_start,i_end
 574
 575  g_mu_save(i,j) =g_mu_2(i,j)
 576  mu_save(i,j) =mu_2(i,j)
 577
 578  g_mu_2(i,j) =g_mu_1(i,j) -g_mu_2(i,j)
 579  mu_2(i,j) =mu_1(i,j) -mu_2(i,j)
 580
 581  ENDDO
 582  ENDDO
 583  END IF
 584
 585  DO j =j_start,j_end
 586  DO i =i_start,i_end
 587
 588  g_ww_save(i,kde,j) =0.0
 589  ww_save(i,kde,j) =0.
 590
 591  g_ww_save(i,1,j) =0.0
 592  ww_save(i,1,j) =0.
 593
 594  ENDDO
 595  ENDDO
 596
 597  DO j =j_start,j_end
 598  DO k =k_start,k_end
 599  DO i =i_start,i_end
 600
 601  g_c2a(i,k,j) =(cpovcv*(g_p(i,k,j))*alt(i,k,j) -g_alt(i,k,j)*cpovcv*(pb(i,k, &
 602  j) +p(i,k,j)))/(alt(i,k,j)*alt(i,k,j))
 603  c2a(i,k,j) =cpovcv*(pb(i,k,j) +p(i,k,j))/alt(i,k,j)
 604
 605  ENDDO
 606  ENDDO
 607  ENDDO
 608
 609  DO j =j_start,j_end
 610  DO k =k_start,k_end
 611  DO i =i_start,i_endu
 612
 613  g_u_save(i,k,j) =g_u_2(i,k,j)
 614  u_save(i,k,j) =u_2(i,k,j)
 615
 616  g_u_2(i,k,j) =(muus(i,j)*g_u_1(i,k,j) +g_muus(i,j)*u_1(i,k,j) -(muu(i,j) &
 617 *g_u_2(i,k,j) +g_muu(i,j)*u_2(i,k,j)))/msfuy(i,j)
 618  u_2(i,k,j) =(muus(i,j)*u_1(i,k,j) -muu(i,j)*u_2(i,k,j))/msfuy(i,j)
 619
 620  ENDDO
 621  ENDDO
 622  ENDDO
 623
 624  DO j =j_start,j_endv
 625  DO k =k_start,k_end
 626  DO i =i_start,i_end
 627
 628  g_v_save(i,k,j) =g_v_2(i,k,j)
 629  v_save(i,k,j) =v_2(i,k,j)
 630
 631  g_v_2(i,k,j) =(muvs(i,j)*g_v_1(i,k,j) +g_muvs(i,j)*v_1(i,k,j) -(muv(i,j) &
 632 *g_v_2(i,k,j) +g_muv(i,j)*v_2(i,k,j)))*msfvx_inv(i,j)
 633  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)
 634
 635  ENDDO
 636  ENDDO
 637  ENDDO
 638
 639 !LPB[6]
 640  DO j =j_start,j_end
 641  DO k =k_start,k_end
 642  DO i =i_start,i_end
 643
 644  g_t_save(i,k,j) =g_t_2(i,k,j)
 645  t_save(i,k,j) =t_2(i,k,j)
 646
 647  g_t_2(i,k,j) =muts(i,j)*g_t_1(i,k,j) +g_muts(i,j)*t_1(i,k,j) -(mut(i,j) &
 648 *g_t_2(i,k,j) +g_mut(i,j)*t_2(i,k,j))
 649  t_2(i,k,j) =muts(i,j)*t_1(i,k,j) -mut(i,j)*t_2(i,k,j)
 650
 651  ENDDO
 652  ENDDO
 653  ENDDO
 654
 655 !LPB[7]
 656  DO j =j_start,j_end
 657  DO k =k_start,kde
 658  DO i =i_start,i_end
 659
 660  g_w_save(i,k,j) =g_w_2(i,k,j)
 661  w_save(i,k,j) =w_2(i,k,j)
 662
 663  g_w_2(i,k,j) =(muts(i,j)*g_w_1(i,k,j) +g_muts(i,j)*w_1(i,k,j) -(mut(i,j) &
 664 *g_w_2(i,k,j) +g_mut(i,j)*w_2(i,k,j)))/msfty(i,j)
 665  w_2(i,k,j) =(muts(i,j)*w_1(i,k,j) -mut(i,j)*w_2(i,k,j))/msfty(i,j)
 666
 667  g_ph_save(i,k,j) =g_ph_2(i,k,j)
 668  ph_save(i,k,j) =ph_2(i,k,j)
 669
 670  g_ph_2(i,k,j) =g_ph_1(i,k,j) -g_ph_2(i,k,j)
 671  ph_2(i,k,j) =ph_1(i,k,j) -ph_2(i,k,j)
 672
 673  ENDDO
 674  ENDDO
 675  ENDDO
 676
 677  DO j =j_start,j_end
 678  DO k =k_start,kde
 679  DO i =i_start,i_end
 680
 681  g_ww_save(i,k,j) =g_ww(i,k,j)
 682  ww_save(i,k,j) =ww(i,k,j)
 683
 684  ENDDO
 685  ENDDO
 686  ENDDO
 687
 688  END SUBROUTINE g_small_step_prep
 689
 690 !        Generated by TAPENADE     (INRIA, Tropics team)
 691 !  Tapenade 3.5 (r3785) - 22 Mar 2011 18:35
 692 !
 693 !  Differentiation of small_step_finish in forward (tangent) mode:
 694 !   variations   of useful results: ww ph_2 u_2 mu_2 w_2 t_2 v_2
 695 !   with respect to varying inputs: u_save ph_save ww w_save mu_save
 696 !                muvs ph_2 u_2 h_diabatic mu_2 w_2 t_save v_save
 697 !                muts t_2 mut muu v_2 muv ww1 muus
 698 !   RW status of diff variables: u_save:in ph_save:in ww:in-out
 699 !                w_save:in mu_save:in muvs:in ph_2:in-out u_2:in-out
 700 !                h_diabatic:in mu_2:in-out w_2:in-out t_save:in
 701 !                v_save:in muts:in t_2:in-out mut:in muu:in v_2:in-out
 702 !                muv:in ww1:in muus:in
 703 SUBROUTINE G_SMALL_STEP_FINISH(u_2, u_2d, u_1, v_2, v_2d, v_1, w_2, w_2d&
 704 &  , w_1, t_2, t_2d, t_1, ph_2, ph_2d, ph_1, ww, wwd, ww1, ww1d, mu_2, &
 705 &  mu_2d, mu_1, mut, mutd, muts, mutsd, muu, muud, muus, muusd, muv, muvd&
 706 &  , muvs, muvsd, u_save, u_saved, v_save, v_saved, w_save, w_saved, &
 707 &  t_save, t_saved, ph_save, ph_saved, mu_save, mu_saved, msfux, msfuy, &
 708 &  msfvx, msfvy, msftx, msfty, h_diabatic, h_diabaticd, &
 709 &  number_of_small_timesteps, dts, rk_step, rk_order, ids, ide, jds, jde&
 710 &  , kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte&
 711 &)
 712   IMPLICIT NONE
 713 ! religion first
 714 !  stuff passed in
 715   INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde
 716   INTEGER, INTENT(IN) :: ims, ime, jms, jme, kms, kme
 717   INTEGER, INTENT(IN) :: its, ite, jts, jte, kts, kte
 718   INTEGER, INTENT(IN) :: number_of_small_timesteps
 719   INTEGER, INTENT(IN) :: rk_step, rk_order
 720   REAL, INTENT(IN) :: dts
 721   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u_1, v_1, &
 722 &  w_1, t_1, ww1, ph_1
 723   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: ww1d
 724   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: u_2, v_2&
 725 &  , w_2, t_2, ww, ph_2
 726   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: u_2d, &
 727 &  v_2d, w_2d, t_2d, wwd, ph_2d
 728   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u_save, &
 729 &  v_save, w_save, t_save, ph_save, h_diabatic
 730   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u_saved, &
 731 &  v_saved, w_saved, t_saved, ph_saved, h_diabaticd
 732   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: muus, muvs
 733   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: muusd, muvsd
 734   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mu_2, mu_1
 735   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mu_2d
 736   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mut, muts, muu, &
 737 &  muv, mu_save
 738   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mutd, mutsd, muud&
 739 &  , muvd, mu_saved
 740   REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: msfux, msfuy, msfvx, &
 741 &  msfvy, msftx, msfty
 742 ! local stuff
 743   INTEGER :: i, j, k
 744   INTEGER :: i_start, i_end, j_start, j_end, i_endu, j_endv
 745   INTRINSIC MIN
 746 !<DESCRIPTION>
 747 !
 748 !  small_step_finish reconstructs the full uncoupled prognostic variables
 749 !  from the coupled perturbation variables used in the small timesteps.
 750 !
 751 !</DESCRIPTION>
 752   i_start = its
 753   IF (ite .GT. ide - 1) THEN
 754     i_end = ide - 1
 755   ELSE
 756     i_end = ite
 757   END IF
 758   j_start = jts
 759   IF (jte .GT. jde - 1) THEN
 760     j_end = jde - 1
 761   ELSE
 762     j_end = jte
 763   END IF
 764   i_endu = ite
 765   j_endv = jte
 766 !    addition of time level t back into variables
 767   DO j=j_start,j_endv
 768     DO k=kds,kde-1
 769       DO i=i_start,i_end
 770 ! v coupled with mx
 771         v_2d(i, k, j) = ((msfvx(i, j)*v_2d(i, k, j)+v_saved(i, k, j)*muv&
 772 &          (i, j)+v_save(i, k, j)*muvd(i, j))*muvs(i, j)-(msfvx(i, j)*v_2&
 773 &          (i, k, j)+v_save(i, k, j)*muv(i, j))*muvsd(i, j))/muvs(i, j)**&
 774 &          2
 775         v_2(i, k, j) = (msfvx(i, j)*v_2(i, k, j)+v_save(i, k, j)*muv(i, &
 776 &          j))/muvs(i, j)
 777       END DO
 778     END DO
 779   END DO
 780   DO j=j_start,j_end
 781     DO k=kds,kde-1
 782       DO i=i_start,i_endu
 783 ! u coupled with my
 784         u_2d(i, k, j) = ((msfuy(i, j)*u_2d(i, k, j)+u_saved(i, k, j)*muu&
 785 &          (i, j)+u_save(i, k, j)*muud(i, j))*muus(i, j)-(msfuy(i, j)*u_2&
 786 &          (i, k, j)+u_save(i, k, j)*muu(i, j))*muusd(i, j))/muus(i, j)**&
 787 &          2
 788         u_2(i, k, j) = (msfuy(i, j)*u_2(i, k, j)+u_save(i, k, j)*muu(i, &
 789 &          j))/muus(i, j)
 790       END DO
 791     END DO
 792   END DO
 793   DO j=j_start,j_end
 794     DO k=kds,kde
 795       DO i=i_start,i_end
 796 ! w coupled with my
 797         w_2d(i, k, j) = ((msfty(i, j)*w_2d(i, k, j)+w_saved(i, k, j)*mut&
 798 &          (i, j)+w_save(i, k, j)*mutd(i, j))*muts(i, j)-(msfty(i, j)*w_2&
 799 &          (i, k, j)+w_save(i, k, j)*mut(i, j))*mutsd(i, j))/muts(i, j)**&
 800 &          2
 801         w_2(i, k, j) = (msfty(i, j)*w_2(i, k, j)+w_save(i, k, j)*mut(i, &
 802 &          j))/muts(i, j)
 803         ph_2d(i, k, j) = ph_2d(i, k, j) + ph_saved(i, k, j)
 804         ph_2(i, k, j) = ph_2(i, k, j) + ph_save(i, k, j)
 805         wwd(i, k, j) = wwd(i, k, j) + ww1d(i, k, j)
 806         ww(i, k, j) = ww(i, k, j) + ww1(i, k, j)
 807       END DO
 808     END DO
 809   END DO
 810   IF (rk_step .LT. rk_order) THEN
 811     DO j=j_start,j_end
 812       DO k=kds,kde-1
 813         DO i=i_start,i_end
 814           t_2d(i, k, j) = ((t_2d(i, k, j)+t_saved(i, k, j)*mut(i, j)+&
 815 &            t_save(i, k, j)*mutd(i, j))*muts(i, j)-(t_2(i, k, j)+t_save(&
 816 &            i, k, j)*mut(i, j))*mutsd(i, j))/muts(i, j)**2
 817           t_2(i, k, j) = (t_2(i, k, j)+t_save(i, k, j)*mut(i, j))/muts(i&
 818 &            , j)
 819         END DO
 820       END DO
 821     END DO
 822   ELSE
 823     DO j=j_start,j_end
 824       DO k=kds,kde-1
 825         DO i=i_start,i_end
 826           t_2d(i, k, j) = ((t_2d(i, k, j)-dts*number_of_small_timesteps*&
 827 &            (mutd(i, j)*h_diabatic(i, k, j)+mut(i, j)*h_diabaticd(i, k, &
 828 &            j))+t_saved(i, k, j)*mut(i, j)+t_save(i, k, j)*mutd(i, j))*&
 829 &            muts(i, j)-(t_2(i, k, j)-dts*number_of_small_timesteps*mut(i&
 830 &            , j)*h_diabatic(i, k, j)+t_save(i, k, j)*mut(i, j))*mutsd(i&
 831 &            , j))/muts(i, j)**2
 832           t_2(i, k, j) = (t_2(i, k, j)-dts*number_of_small_timesteps*mut&
 833 &            (i, j)*h_diabatic(i, k, j)+t_save(i, k, j)*mut(i, j))/muts(i&
 834 &            , j)
 835         END DO
 836       END DO
 837     END DO
 838   END IF
 839   DO j=j_start,j_end
 840     DO i=i_start,i_end
 841       mu_2d(i, j) = mu_2d(i, j) + mu_saved(i, j)
 842       mu_2(i, j) = mu_2(i, j) + mu_save(i, j)
 843     END DO
 844   END DO
 845 END SUBROUTINE G_SMALL_STEP_FINISH
 846
 847  SUBROUTINE g_calc_p_rho(al,g_al,p,g_p,ph,g_ph,alt,g_alt,t_2,g_t_2, &
 848  t_1,g_t_1,c2a,g_c2a,pm1,g_pm1,mu,g_mu,muts,g_muts,znu,t0,rdnw,dnw, &
 849  smdiv,non_hydrostatic,step,ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite, &
 850  jts,jte,kts,kte)
 851
 852  IMPLICIT NONE
 853  INTEGER,INTENT(IN) :: ids,ide,jds,jde,kds,kde
 854  INTEGER,INTENT(IN) :: ims,ime,jms,jme,kms,kme
 855  INTEGER,INTENT(IN) :: its,ite,jts,jte,kts,kte
 856  INTEGER,INTENT(IN) :: step
 857  REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(OUT) :: g_al,al,g_p,p
 858  REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(IN) :: g_alt,alt,g_t_2,t_2, &
 859  g_t_1,t_1,g_c2a,c2a
 860  REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(INOUT) :: g_ph,ph,g_pm1,pm1
 861  REAL,DIMENSION(ims:ime,jms:jme),INTENT(IN) :: g_mu,mu,g_muts,muts
 862  REAL,DIMENSION(kms:kme),INTENT(IN) :: dnw,rdnw,znu
 863  REAL,INTENT(IN) :: t0,smdiv
 864  LOGICAL,INTENT(IN) :: non_hydrostatic
 865  INTEGER :: i,j,k
 866  INTEGER :: i_start,i_end,j_start,j_end,k_start,k_end
 867  REAL :: g_ptmp,ptmp
 868
 869  i_start =its
 870
 871  i_end =min(ite,ide-1)
 872
 873  j_start =jts
 874
 875  j_end =min(jte,jde-1)
 876
 877  k_start =kts
 878
 879  k_end =min(kte,kde-1)
 880
 881  IF(non_hydrostatic) THEN
 882
 883  DO j =j_start,j_end
 884  DO k =k_start,k_end
 885  DO i =i_start,i_end
 886
 887  g_al(i,k,j) =-1./muts(i,j)*(alt(i,k,j)*g_mu(i,j) +g_alt(i,k,j)*mu(i,j) &
 888  +rdnw(k)*(g_ph(i,k+1,j) -g_ph(i,k,j))) +(-(-1.)*g_muts(i,j)/(muts(i,j) &
 889 *muts(i,j)))*(alt(i,k,j)*mu(i,j) +rdnw(k)*(ph(i,k+1,j) -ph(i,k,j)))
 890  al(i,k,j) =-1./muts(i,j)*(alt(i,k,j)*mu(i,j) +rdnw(k)*(ph(i,k+1,j) -ph(i,k,j)))
 891
 892  g_p(i,k,j) =c2a(i,k,j)*(((alt(i,k,j)*(g_t_2(i,k,j) -(mu(i,j)*g_t_1(i,k,j) &
 893  +g_mu(i,j)*t_1(i,k,j))) +g_alt(i,k,j)*(t_2(i,k,j) -mu(i,j)*t_1(i,k,j))) &
 894 *(muts(i,j)*(t0 +t_1(i,k,j))) -(muts(i,j)*(g_t_1(i,k,j)) +g_muts(i,j) &
 895 *(t0 +t_1(i,k,j)))*alt(i,k,j)*(t_2(i,k,j) -mu(i,j)*t_1(i,k,j)))/((muts(i,j) &
 896 *(t0 +t_1(i,k,j)))*(muts(i,j)*(t0 +t_1(i,k,j)))) -g_al(i,k,j)) +g_c2a(i,k,j) &
 897 *(alt(i,k,j)*(t_2(i,k,j) -mu(i,j)*t_1(i,k,j))/(muts(i,j)*(t0 +t_1(i,k,j))) -al(i,k,j))
 898  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))/(muts(i,j) &
 899 *(t0 +t_1(i,k,j))) -al(i,k,j))
 900
 901  ENDDO
 902  ENDDO
 903  ENDDO
 904
 905  ELSE
 906
 907  DO j =j_start,j_end
 908  DO k =k_start,k_end
 909  DO i =i_start,i_end
 910
 911  g_p(i,k,j) =g_mu(i,j)*znu(k)
 912  p(i,k,j) =mu(i,j)*znu(k)
 913
 914  g_al(i,k,j) =((alt(i,k,j)*(g_t_2(i,k,j) -(mu(i,j)*g_t_1(i,k,j) +g_mu(i,j) &
 915 *t_1(i,k,j))) +g_alt(i,k,j)*(t_2(i,k,j) -mu(i,j)*t_1(i,k,j)))*(muts(i,j) &
 916 *(t0 +t_1(i,k,j))) -(muts(i,j)*(g_t_1(i,k,j)) +g_muts(i,j)*(t0 +t_1(i,k,j))) &
 917 *alt(i,k,j)*(t_2(i,k,j) -mu(i,j)*t_1(i,k,j)))/((muts(i,j)*(t0 +t_1(i,k,j))) &
 918 *(muts(i,j)*(t0 +t_1(i,k,j)))) -((g_p(i,k,j)*c2a(i,k,j) -g_c2a(i,k,j)*p(i,k,j)) &
 919 /(c2a(i,k,j)*c2a(i,k,j)))
 920  al(i,k,j) =alt(i,k,j)*(t_2(i,k,j) -mu(i,j)*t_1(i,k,j))/(muts(i,j)*(t0 +t_1(i,k,j))) &
 921  -p(i,k,j)/c2a(i,k,j)
 922
 923  g_ph(i,k+1,j) =g_ph(i,k,j) -dnw(k)*(muts(i,j)*g_al(i,k,j) +g_muts(i,j) &
 924 *al(i,k,j) +(mu(i,j)*g_alt(i,k,j) +g_mu(i,j)*alt(i,k,j)))
 925  ph(i,k+1,j) =ph(i,k,j) -dnw(k)*(muts(i,j)*al(i,k,j) +mu(i,j)*alt(i,k,j))
 926
 927  ENDDO
 928  ENDDO
 929  ENDDO
 930  END IF
 931
 932  IF(step == 0) THEN
 933
 934  DO j =j_start,j_end
 935  DO k =k_start,k_end
 936  DO i =i_start,i_end
 937
 938  g_pm1(i,k,j) =g_p(i,k,j)
 939  pm1(i,k,j) =p(i,k,j)
 940
 941  ENDDO
 942  ENDDO
 943  ENDDO
 944
 945  ELSE
 946
 947  DO j =j_start,j_end
 948  DO k =k_start,k_end
 949  DO i =i_start,i_end
 950
 951  g_ptmp =g_p(i,k,j)
 952  ptmp =p(i,k,j)
 953
 954  g_p(i,k,j) =g_p(i,k,j) +smdiv*(g_p(i,k,j) -g_pm1(i,k,j))
 955  p(i,k,j) =p(i,k,j) +smdiv*(p(i,k,j) -pm1(i,k,j))
 956
 957  g_pm1(i,k,j) =g_ptmp
 958  pm1(i,k,j) =ptmp
 959
 960  ENDDO
 961  ENDDO
 962  ENDDO
 963  END IF
 964
 965  END SUBROUTINE g_calc_p_rho
 966
 967  SUBROUTINE g_calc_coef_w(a,g_a,alpha,g_alpha,gamma,g_gamma,mut,g_mut, &
 968  cqw,g_cqw,rdn,rdnw,c2a,g_c2a,dts,g,epssm,top_lid,ids,ide,jds,jde,kds,kde,ims, &
 969  ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte)
 970
 971  IMPLICIT NONE
 972  INTEGER,INTENT(IN) :: ids,ide,jds,jde,kds,kde
 973  INTEGER,INTENT(IN) :: ims,ime,jms,jme,kms,kme
 974  INTEGER,INTENT(IN) :: its,ite,jts,jte,kts,kte
 975  LOGICAL,INTENT(IN) :: top_lid
 976  REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(IN) :: g_c2a,c2a,g_cqw,cqw
 977  REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(INOUT) :: g_alpha,alpha,g_gamma, &
 978  gamma,g_a,a
 979  REAL,DIMENSION(ims:ime,jms:jme),INTENT(IN) :: g_mut,mut
 980  REAL,DIMENSION(kms:kme),INTENT(IN) :: rdn,rdnw
 981  REAL,INTENT(IN) :: epssm,dts,g
 982  REAL,DIMENSION(ims:ime) :: g_cof,cof
 983  REAL :: g_b,b,g_c,c
 984  INTEGER :: i,j,k,i_start,i_end,j_start,j_end,k_start,k_end
 985  INTEGER :: ij,ijp,ijm,lid_flag
 986
 987  i_start =its
 988
 989  i_end =min(ite,ide-1)
 990
 991  j_start =jts
 992
 993  j_end =min(jte,jde-1)
 994
 995  k_start =kts
 996
 997  k_end =kte-1
 998
 999  lid_flag =1
1000
1001  IF(top_lid) lid_flag =0
1002
1003  DO j =j_start,j_end
1004  DO i =i_start,i_end
1005
1006  g_cof(i) =2.0*(-.5 *dts *g *(1.+epssm)*g_mut(i,j)/(mut(i,j)*mut(i,j))) &
1007 *(.5 *dts *g *(1.+epssm)/mut(i,j))
1008  cof(i) =(.5 *dts *g *(1.+epssm)/mut(i,j))**2
1009
1010  g_a(i,2,j) =0.0
1011  a(i,2,j) =0.
1012
1013  g_a(i,kde,j) =(-2.*cof(i)*rdnw(kde-1)**2*g_c2a(i,kde-1,j) +(-2.)*g_cof(i) &
1014 *rdnw(kde-1)**2*c2a(i,kde-1,j))*lid_flag
1015  a(i,kde,j) =-2.*cof(i)*rdnw(kde-1)**2*c2a(i,kde-1,j)*lid_flag
1016
1017  g_gamma(i,1,j) =0.0
1018  gamma(i,1,j) =0.
1019
1020  ENDDO
1021
1022  DO k =3,kde-1
1023  DO i =i_start,i_end
1024
1025  g_a(i,k,j) =-cqw(i,k,j)*cof(i)*rdn(k)*rdnw(k-1)*g_c2a(i,k-1,j) +(-cqw(i,k,j) &
1026 *g_cof(i) -g_cqw(i,k,j)*cof(i))*rdn(k)*rdnw(k-1)*c2a(i,k-1,j)
1027  a(i,k,j) =-cqw(i,k,j)*cof(i)*rdn(k)*rdnw(k-1)*c2a(i,k-1,j)
1028
1029  ENDDO
1030  ENDDO
1031
1032  DO k =2,kde-1
1033  DO i =i_start,i_end
1034
1035  g_b =cqw(i,k,j)*cof(i)*rdn(k)*(rdnw(k)*g_c2a(i,k,j) +rdnw(k-1)*g_c2a(i,k-1, &
1036  j)) +(cqw(i,k,j)*g_cof(i) +g_cqw(i,k,j)*cof(i))*rdn(k)*(rdnw(k)*c2a(i,k,j) &
1037  +rdnw(k-1)*c2a(i,k-1,j))
1038  b =1. +cqw(i,k,j)*cof(i)*rdn(k)*(rdnw(k)*c2a(i,k,j) +rdnw(k-1)*c2a(i,k-1,j))
1039
1040  g_c =-cqw(i,k,j)*cof(i)*rdn(k)*rdnw(k)*g_c2a(i,k,j) +(-cqw(i,k,j)*g_cof(i) &
1041  -g_cqw(i,k,j)*cof(i))*rdn(k)*rdnw(k)*c2a(i,k,j)
1042  c =-cqw(i,k,j)*cof(i)*rdn(k)*rdnw(k)*c2a(i,k,j)
1043
1044  g_alpha(i,k,j) =-1.*(g_b -(a(i,k,j)*g_gamma(i,k-1,j) +g_a(i,k,j) &
1045 *gamma(i,k-1,j)))/((b -a(i,k,j)*gamma(i,k-1,j))*(b -a(i,k,j)*gamma(i,k-1,j)))
1046  alpha(i,k,j) =1./(b -a(i,k,j)*gamma(i,k-1,j))
1047
1048  g_gamma(i,k,j) =c*g_alpha(i,k,j) +g_c*alpha(i,k,j)
1049  gamma(i,k,j) =c*alpha(i,k,j)
1050
1051  ENDDO
1052  ENDDO
1053
1054  DO i =i_start,i_end
1055
1056  g_b =2.*cof(i)*rdnw(kde-1)**2*g_c2a(i,kde-1,j) +2.*g_cof(i)*rdnw(kde-1) &
1057 **2*c2a(i,kde-1,j)
1058  b =1. +2.*cof(i)*rdnw(kde-1)**2*c2a(i,kde-1,j)
1059
1060  g_c =0.0
1061  c =0.
1062
1063  g_alpha(i,kde,j) =-1.*(g_b -(a(i,kde,j)*g_gamma(i,kde-1,j) +g_a(i,kde,j) &
1064 *gamma(i,kde-1,j)))/((b -a(i,kde,j)*gamma(i,kde-1,j))*(b -a(i,kde,j)*gamma(i,kde-1,j)))
1065  alpha(i,kde,j) =1./(b -a(i,kde,j)*gamma(i,kde-1,j))
1066
1067  g_gamma(i,kde,j) =c*g_alpha(i,kde,j) +g_c*alpha(i,kde,j)
1068  gamma(i,kde,j) =c*alpha(i,kde,j)
1069
1070  ENDDO
1071  ENDDO
1072
1073  END SUBROUTINE g_calc_coef_w
1074
1075 !        Generated by TAPENADE     (INRIA, Tropics team)
1076 !  Tapenade 3.5 (r3785) - 22 Mar 2011 18:35
1077 !
1078 !  Differentiation of advance_uv in forward (tangent) mode:
1079 !   variations   of useful results: u v
1080 !   with respect to varying inputs: p al u v ru_tend mudf cqu cqv
1081 !                php rv_tend ph alt muu muv mu
1082 !   RW status of diff variables: p:in al:in u:in-out v:in-out ru_tend:in
1083 !                mudf:in cqu:in cqv:in php:in rv_tend:in ph:in
1084 !                alt:in muu:in muv:in mu:in
1085 SUBROUTINE G_ADVANCE_UV(u, ud, ru_tend, ru_tendd, v, vd, rv_tend, &
1086 &  rv_tendd, p, pd, pb, ph, phd, php, phpd, alt, altd, al, ald, mu, mud, &
1087 &  muu, muud, cqu, cqud, muv, muvd, cqv, cqvd, mudf, mudfd, msfux, msfuy&
1088 &  , msfvx, msfvx_inv, msfvy, rdx, rdy, dts, cf1, cf2, cf3, fnm, fnp, &
1089 &  emdiv, rdnw, config_flags, spec_zone, non_hydrostatic, top_lid, ids, &
1090 &  ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, &
1091 &  jte, kts, kte)
1092   IMPLICIT NONE
1093 ! religion first
1094 ! stuff coming in
1095   TYPE(GRID_CONFIG_REC_TYPE), INTENT(IN) :: config_flags
1096   LOGICAL, INTENT(IN) :: non_hydrostatic, top_lid
1097   INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde
1098   INTEGER, INTENT(IN) :: ims, ime, jms, jme, kms, kme
1099   INTEGER, INTENT(IN) :: its, ite, jts, jte, kts, kte
1100   INTEGER, INTENT(IN) :: spec_zone
1101   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: u, v
1102   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: ud, vd
1103   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: ru_tend, &
1104 &  rv_tend, ph, php, p, pb, alt, al, cqu, cqv
1105   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: ru_tendd, &
1106 &  rv_tendd, phd, phpd, pd, altd, ald, cqud, cqvd
1107   REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: muu, muv, mu, mudf
1108   REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: muud, muvd, mud, &
1109 &  mudfd
1110   REAL, DIMENSION(kms:kme), INTENT(IN) :: fnm, fnp, rdnw
1111   REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: msfux, msfuy, msfvx, &
1112 &  msfvy, msfvx_inv
1113   REAL, INTENT(IN) :: rdx, rdy, dts, cf1, cf2, cf3, emdiv
1114 !  Local 3d array from the stack (note tile size)
1115   REAL, DIMENSION(its:ite, kts:kte) :: dpn, dpxy
1116   REAL, DIMENSION(its:ite, kts:kte) :: dpnd, dpxyd
1117   REAL, DIMENSION(its:ite) :: mudf_xy
1118   REAL, DIMENSION(its:ite) :: mudf_xyd
1119   REAL :: dx, dy
1120   INTEGER :: i, j, k, i_start, i_end, j_start, j_end, k_start, k_end
1121   INTEGER :: i_endu, j_endv, k_endw
1122   INTEGER :: i_start_up, i_end_up, j_start_up, j_end_up
1123   INTEGER :: i_start_vp, i_end_vp, j_start_vp, j_end_vp
1124   INTEGER :: i_start_u_tend, i_end_u_tend, j_start_v_tend, j_end_v_tend
1125 !<DESCRIPTION>
1126 !
1127 !  advance_uv advances the explicit perturbation horizontal momentum
1128 !  equations (u,v) by adding in the large-timestep tendency along with
1129 !  the small timestep pressure gradient tendency.
1130 !
1131 !</DESCRIPTION>
1132 !  now, the real work.
1133 !  set the loop bounds taking into account boundary conditions.
1134   IF (config_flags%nested .OR. config_flags%specified) THEN
1135     IF (its .LT. ids + spec_zone) THEN
1136       i_start = ids + spec_zone
1137     ELSE
1138       i_start = its
1139     END IF
1140     IF (ite .GT. ide - spec_zone - 1) THEN
1141       i_end = ide - spec_zone - 1
1142     ELSE
1143       i_end = ite
1144     END IF
1145     IF (jts .LT. jds + spec_zone) THEN
1146       j_start = jds + spec_zone
1147     ELSE
1148       j_start = jts
1149     END IF
1150     IF (jte .GT. jde - spec_zone - 1) THEN
1151       j_end = jde - spec_zone - 1
1152     ELSE
1153       j_end = jte
1154     END IF
1155     k_start = kts
1156     IF (kte .GT. kde - 1) THEN
1157       k_end = kde - 1
1158     ELSE
1159       k_end = kte
1160     END IF
1161     IF (ite .GT. ide - spec_zone) THEN
1162       i_endu = ide - spec_zone
1163     ELSE
1164       i_endu = ite
1165     END IF
1166     IF (jte .GT. jde - spec_zone) THEN
1167       j_endv = jde - spec_zone
1168     ELSE
1169       j_endv = jte
1170     END IF
1171     k_endw = k_end
1172     IF (config_flags%periodic_x) THEN
1173       i_start = its
1174       IF (ite .GT. ide - 1) THEN
1175         i_end = ide - 1
1176       ELSE
1177         i_end = ite
1178       END IF
1179       i_endu = ite
1180     END IF
1181   ELSE
1182     i_start = its
1183     IF (ite .GT. ide - 1) THEN
1184       i_end = ide - 1
1185     ELSE
1186       i_end = ite
1187     END IF
1188     j_start = jts
1189     IF (jte .GT. jde - 1) THEN
1190       j_end = jde - 1
1191     ELSE
1192       j_end = jte
1193     END IF
1194     k_start = kts
1195     k_end = kte - 1
1196     i_endu = ite
1197     j_endv = jte
1198     k_endw = k_end
1199   END IF
1200   i_start_up = i_start
1201   i_end_up = i_endu
1202   j_start_up = j_start
1203   j_end_up = j_end
1204   i_start_vp = i_start
1205   i_end_vp = i_end
1206   j_start_vp = j_start
1207   j_end_vp = j_endv
1208   IF ((config_flags%open_xs .OR. config_flags%symmetric_xs) .AND. its &
1209 &      .EQ. ids) i_start_up = i_start_up + 1
1210   IF ((config_flags%open_xe .OR. config_flags%symmetric_xe) .AND. ite &
1211 &      .EQ. ide) i_end_up = i_end_up - 1
1212   IF (((config_flags%open_ys .OR. config_flags%symmetric_ys) .OR. &
1213 &      config_flags%polar) .AND. jts .EQ. jds) j_start_vp = j_start_vp + &
1214 &      1
1215   IF (((config_flags%open_ye .OR. config_flags%symmetric_ye) .OR. &
1216 &      config_flags%polar) .AND. jte .EQ. jde) j_end_vp = j_end_vp - 1
1217   i_start_u_tend = i_start
1218   i_end_u_tend = i_endu
1219   j_start_v_tend = j_start
1220   j_end_v_tend = j_endv
1221   IF (config_flags%symmetric_xs .AND. its .EQ. ids) i_start_u_tend = &
1222 &      i_start_u_tend + 1
1223   IF (config_flags%symmetric_xe .AND. ite .EQ. ide) i_end_u_tend = &
1224 &      i_end_u_tend - 1
1225   IF (config_flags%symmetric_ys .AND. jts .EQ. jds) j_start_v_tend = &
1226 &      j_start_v_tend + 1
1227   IF (config_flags%symmetric_ye .AND. jte .EQ. jde) j_end_v_tend = &
1228 &      j_end_v_tend - 1
1229   dx = 1./rdx
1230   dy = 1./rdy
1231   mudf_xyd = 0.0
1232   dpxyd = 0.0
1233   dpnd = 0.0
1234 !  start real calculations.
1235 !  first, u
1236 u_outer_j_loop:DO j=j_start,j_end
1237     DO k=k_start,k_end
1238       DO i=i_start_u_tend,i_end_u_tend
1239         ud(i, k, j) = ud(i, k, j) + dts*ru_tendd(i, k, j)
1240         u(i, k, j) = u(i, k, j) + dts*ru_tend(i, k, j)
1241       END DO
1242     END DO
1243     DO i=i_start_up,i_end_up
1244       mudf_xyd(i) = -(emdiv*dx*(mudfd(i, j)-mudfd(i-1, j))/msfuy(i, j))
1245       mudf_xy(i) = -(emdiv*dx*(mudf(i, j)-mudf(i-1, j))/msfuy(i, j))
1246     END DO
1247     DO k=k_start,k_end
1248       DO i=i_start_up,i_end_up
1249 !     Comments on map scale factors:
1250 !     x pressure gradient: ADT eqn 44, penultimate term on RHS
1251 !        = -(mx/my)*(mu/rho)*partial dp/dx
1252 !     [i.e., first rho->mu; 2nd still rho; alpha=1/rho]
1253 !     Klemp et al. splits into 2 terms:
1254 !        mu alpha partial dp/dx + partial dp/dnu * partial dphi/dx
1255 !     then into 4 terms:
1256 !        mu alpha partial dp'/dx + nu mu alpha' partial dmubar/dx +
1257 !      + mu partial dphi/dx + partial dphi'/dx * (partial dp'/dnu - mu')
1258 !
1259 !     first 3 terms:
1260 !     ph, alt, p, al, pb not coupled
1261 !     since we want tendency to fit ADT eqn 44 (coupled) we need to
1262 !     multiply by (mx/my):
1263 !
1264         dpxyd(i, k) = msfux(i, j)*.5*rdx*(muud(i, j)*(ph(i, k+1, j)-ph(i&
1265 &          -1, k+1, j)+(ph(i, k, j)-ph(i-1, k, j))+(alt(i, k, j)+alt(i-1&
1266 &          , k, j))*(p(i, k, j)-p(i-1, k, j))+(al(i, k, j)+al(i-1, k, j))&
1267 &          *(pb(i, k, j)-pb(i-1, k, j)))+muu(i, j)*(phd(i, k+1, j)-phd(i-&
1268 &          1, k+1, j)+phd(i, k, j)-phd(i-1, k, j)+(altd(i, k, j)+altd(i-1&
1269 &          , k, j))*(p(i, k, j)-p(i-1, k, j))+(alt(i, k, j)+alt(i-1, k, j&
1270 &          ))*(pd(i, k, j)-pd(i-1, k, j))+(pb(i, k, j)-pb(i-1, k, j))*(&
1271 &          ald(i, k, j)+ald(i-1, k, j))))/msfuy(i, j)
1272         dpxy(i, k) = msfux(i, j)/msfuy(i, j)*.5*rdx*muu(i, j)*(ph(i, k+1&
1273 &          , j)-ph(i-1, k+1, j)+(ph(i, k, j)-ph(i-1, k, j))+(alt(i, k, j)&
1274 &          +alt(i-1, k, j))*(p(i, k, j)-p(i-1, k, j))+(al(i, k, j)+al(i-1&
1275 &          , k, j))*(pb(i, k, j)-pb(i-1, k, j)))
1276       END DO
1277     END DO
1278     IF (non_hydrostatic) THEN
1279       DO i=i_start_up,i_end_up
1280         dpnd(i, 1) = .5*(cf1*(pd(i, 1, j)+pd(i-1, 1, j))+cf2*(pd(i, 2, j&
1281 &          )+pd(i-1, 2, j))+cf3*(pd(i, 3, j)+pd(i-1, 3, j)))
1282         dpn(i, 1) = .5*(cf1*(p(i, 1, j)+p(i-1, 1, j))+cf2*(p(i, 2, j)+p(&
1283 &          i-1, 2, j))+cf3*(p(i, 3, j)+p(i-1, 3, j)))
1284         dpnd(i, kde) = 0.0
1285         dpn(i, kde) = 0.
1286       END DO
1287       IF (top_lid) THEN
1288         DO i=i_start_up,i_end_up
1289           dpnd(i, kde) = .5*(cf1*(pd(i-1, kde-1, j)+pd(i, kde-1, j))+cf2&
1290 &            *(pd(i-1, kde-2, j)+pd(i, kde-2, j))+cf3*(pd(i-1, kde-3, j)+&
1291 &            pd(i, kde-3, j)))
1292           dpn(i, kde) = .5*(cf1*(p(i-1, kde-1, j)+p(i, kde-1, j))+cf2*(p&
1293 &            (i-1, kde-2, j)+p(i, kde-2, j))+cf3*(p(i-1, kde-3, j)+p(i, &
1294 &            kde-3, j)))
1295         END DO
1296       END IF
1297       DO k=k_start+1,k_end
1298         DO i=i_start_up,i_end_up
1299           dpnd(i, k) = .5*(fnm(k)*(pd(i, k, j)+pd(i-1, k, j))+fnp(k)*(pd&
1300 &            (i, k-1, j)+pd(i-1, k-1, j)))
1301           dpn(i, k) = .5*(fnm(k)*(p(i, k, j)+p(i-1, k, j))+fnp(k)*(p(i, &
1302 &            k-1, j)+p(i-1, k-1, j)))
1303         END DO
1304       END DO
1305 !    Comments on map scale factors:
1306 !    4th term:
1307 !    php, dpn, mu not coupled
1308 !    since we want tendency to fit ADT eqn 44 (coupled) we need to
1309 !    multiply by (mx/my):
1310       DO k=k_start,k_end
1311         DO i=i_start_up,i_end_up
1312           dpxyd(i, k) = dpxyd(i, k) + msfux(i, j)*rdx*((phpd(i, k, j)-&
1313 &            phpd(i-1, k, j))*(rdnw(k)*(dpn(i, k+1)-dpn(i, k))-.5*(mu(i-1&
1314 &            , j)+mu(i, j)))+(php(i, k, j)-php(i-1, k, j))*(rdnw(k)*(dpnd&
1315 &            (i, k+1)-dpnd(i, k))-.5*(mud(i-1, j)+mud(i, j))))/msfuy(i, j&
1316 &            )
1317           dpxy(i, k) = dpxy(i, k) + msfux(i, j)/msfuy(i, j)*rdx*(php(i, &
1318 &            k, j)-php(i-1, k, j))*(rdnw(k)*(dpn(i, k+1)-dpn(i, k))-.5*(&
1319 &            mu(i-1, j)+mu(i, j)))
1320         END DO
1321       END DO
1322     END IF
1323     DO k=k_start,k_end
1324       DO i=i_start_up,i_end_up
1325         ud(i, k, j) = ud(i, k, j) - dts*(cqud(i, k, j)*dpxy(i, k)+cqu(i&
1326 &          , k, j)*dpxyd(i, k)) + mudf_xyd(i)
1327         u(i, k, j) = u(i, k, j) - dts*cqu(i, k, j)*dpxy(i, k) + mudf_xy(&
1328 &          i)
1329       END DO
1330     END DO
1331   END DO u_outer_j_loop
1332 ! now v
1333 v_outer_j_loop:DO j=j_start_v_tend,j_end_v_tend
1334     DO k=k_start,k_end
1335       DO i=i_start,i_end
1336         vd(i, k, j) = vd(i, k, j) + dts*rv_tendd(i, k, j)
1337         v(i, k, j) = v(i, k, j) + dts*rv_tend(i, k, j)
1338       END DO
1339     END DO
1340     DO i=i_start,i_end
1341       mudf_xyd(i) = -(emdiv*dy*msfvx_inv(i, j)*(mudfd(i, j)-mudfd(i, j-1&
1342 &        )))
1343       mudf_xy(i) = -(emdiv*dy*(mudf(i, j)-mudf(i, j-1))*msfvx_inv(i, j))
1344     END DO
1345     IF (j .GE. j_start_vp .AND. j .LE. j_end_vp) THEN
1346       DO k=k_start,k_end
1347         DO i=i_start,i_end
1348 !      Comments on map scale factors:
1349 !      y pressure gradient: ADT eqn 45, penultimate term on RHS
1350 !         = -(my/mx)*(mu/rho)*partial dp/dy
1351 !      [i.e., first rho->mu; 2nd still rho; alpha=1/rho]
1352 !      Klemp et al. splits into 2 terms:
1353 !         mu alpha partial dp/dy + partial dp/dnu * partial dphi/dy
1354 !      then into 4 terms:
1355 !         mu alpha partial dp'/dy + nu mu alpha' partial dmubar/dy +
1356 !       + mu partial dphi/dy + partial dphi'/dy * (partial dp'/dnu - mu')
1357 !
1358 !      first 3 terms:
1359 !      ph, alt, p, al, pb not coupled
1360 !      since we want tendency to fit ADT eqn 45 (coupled) we need to
1361 !      multiply by (my/mx):
1362 !      mudf_xy is NOT a map scale factor coupling
1363 !      it is some sort of divergence damping
1364           dpxyd(i, k) = msfvy(i, j)*.5*rdy*(muvd(i, j)*(ph(i, k+1, j)-ph&
1365 &            (i, k+1, j-1)+(ph(i, k, j)-ph(i, k, j-1))+(alt(i, k, j)+alt(&
1366 &            i, k, j-1))*(p(i, k, j)-p(i, k, j-1))+(al(i, k, j)+al(i, k, &
1367 &            j-1))*(pb(i, k, j)-pb(i, k, j-1)))+muv(i, j)*(phd(i, k+1, j)&
1368 &            -phd(i, k+1, j-1)+phd(i, k, j)-phd(i, k, j-1)+(altd(i, k, j)&
1369 &            +altd(i, k, j-1))*(p(i, k, j)-p(i, k, j-1))+(alt(i, k, j)+&
1370 &            alt(i, k, j-1))*(pd(i, k, j)-pd(i, k, j-1))+(pb(i, k, j)-pb(&
1371 &            i, k, j-1))*(ald(i, k, j)+ald(i, k, j-1))))/msfvx(i, j)
1372           dpxy(i, k) = msfvy(i, j)/msfvx(i, j)*.5*rdy*muv(i, j)*(ph(i, k&
1373 &            +1, j)-ph(i, k+1, j-1)+(ph(i, k, j)-ph(i, k, j-1))+(alt(i, k&
1374 &            , j)+alt(i, k, j-1))*(p(i, k, j)-p(i, k, j-1))+(al(i, k, j)+&
1375 &            al(i, k, j-1))*(pb(i, k, j)-pb(i, k, j-1)))
1376         END DO
1377       END DO
1378       IF (non_hydrostatic) THEN
1379         DO i=i_start,i_end
1380           dpnd(i, 1) = .5*(cf1*(pd(i, 1, j)+pd(i, 1, j-1))+cf2*(pd(i, 2&
1381 &            , j)+pd(i, 2, j-1))+cf3*(pd(i, 3, j)+pd(i, 3, j-1)))
1382           dpn(i, 1) = .5*(cf1*(p(i, 1, j)+p(i, 1, j-1))+cf2*(p(i, 2, j)+&
1383 &            p(i, 2, j-1))+cf3*(p(i, 3, j)+p(i, 3, j-1)))
1384           dpnd(i, kde) = 0.0
1385           dpn(i, kde) = 0.
1386         END DO
1387         IF (top_lid) THEN
1388           DO i=i_start,i_end
1389             dpnd(i, kde) = .5*(cf1*(pd(i, kde-1, j-1)+pd(i, kde-1, j))+&
1390 &              cf2*(pd(i, kde-2, j-1)+pd(i, kde-2, j))+cf3*(pd(i, kde-3, &
1391 &              j-1)+pd(i, kde-3, j)))
1392             dpn(i, kde) = .5*(cf1*(p(i, kde-1, j-1)+p(i, kde-1, j))+cf2*&
1393 &              (p(i, kde-2, j-1)+p(i, kde-2, j))+cf3*(p(i, kde-3, j-1)+p(&
1394 &              i, kde-3, j)))
1395           END DO
1396         END IF
1397         DO k=k_start+1,k_end
1398           DO i=i_start,i_end
1399             dpnd(i, k) = .5*(fnm(k)*(pd(i, k, j)+pd(i, k, j-1))+fnp(k)*(&
1400 &              pd(i, k-1, j)+pd(i, k-1, j-1)))
1401             dpn(i, k) = .5*(fnm(k)*(p(i, k, j)+p(i, k, j-1))+fnp(k)*(p(i&
1402 &              , k-1, j)+p(i, k-1, j-1)))
1403           END DO
1404         END DO
1405 !      Comments on map scale factors:
1406 !      4th term:
1407 !      php, dpn, mu not coupled
1408 !      since we want tendency to fit ADT eqn 45 (coupled) we need to
1409 !      multiply by (my/mx):
1410         DO k=k_start,k_end
1411           DO i=i_start,i_end
1412             dpxyd(i, k) = dpxyd(i, k) + msfvy(i, j)*rdy*((phpd(i, k, j)-&
1413 &              phpd(i, k, j-1))*(rdnw(k)*(dpn(i, k+1)-dpn(i, k))-.5*(mu(i&
1414 &              , j-1)+mu(i, j)))+(php(i, k, j)-php(i, k, j-1))*(rdnw(k)*(&
1415 &              dpnd(i, k+1)-dpnd(i, k))-.5*(mud(i, j-1)+mud(i, j))))/&
1416 &              msfvx(i, j)
1417             dpxy(i, k) = dpxy(i, k) + msfvy(i, j)/msfvx(i, j)*rdy*(php(i&
1418 &              , k, j)-php(i, k, j-1))*(rdnw(k)*(dpn(i, k+1)-dpn(i, k))-&
1419 &              .5*(mu(i, j-1)+mu(i, j)))
1420           END DO
1421         END DO
1422       END IF
1423       DO k=k_start,k_end
1424         DO i=i_start,i_end
1425           vd(i, k, j) = vd(i, k, j) - dts*(cqvd(i, k, j)*dpxy(i, k)+cqv(&
1426 &            i, k, j)*dpxyd(i, k)) + mudf_xyd(i)
1427           v(i, k, j) = v(i, k, j) - dts*cqv(i, k, j)*dpxy(i, k) + &
1428 &            mudf_xy(i)
1429         END DO
1430       END DO
1431     END IF
1432   END DO v_outer_j_loop
1433 ! The check for j_start_vp and j_end_vp makes sure that the edges in v
1434 ! are not updated.  Let's add a polar boundary condition check here for
1435 ! safety to ensure that v at the poles never gets to be non-zero in the
1436 ! small time steps.
1437   IF (config_flags%polar) THEN
1438     IF (jts .EQ. jds) THEN
1439       DO k=k_start,k_end
1440         DO i=i_start,i_end
1441           vd(i, k, jds) = 0.0
1442           v(i, k, jds) = 0.
1443         END DO
1444       END DO
1445     END IF
1446     IF (jte .EQ. jde) THEN
1447       DO k=k_start,k_end
1448         DO i=i_start,i_end
1449           vd(i, k, jde) = 0.0
1450           v(i, k, jde) = 0.
1451         END DO
1452       END DO
1453     END IF
1454   END IF
1455 END SUBROUTINE G_ADVANCE_UV
1456
1457  SUBROUTINE g_advance_mu_t(ww,g_ww,ww_1,g_ww_1,u,g_u,u_1,g_u_1,v, &
1458  g_v,v_1,g_v_1,mu,g_mu,mut,g_mut,muave,g_muave,muts,g_muts,muu,g_muu, &
1459  muv,g_muv,mudf,g_mudf,uam,g_uam,vam,g_vam,wwam,g_wwam,t,g_t, &
1460  t_1,g_t_1,t_ave,g_t_ave,ft,g_ft,mu_tend,g_mu_tend,rdx,rdy,dts,epssm,dnw, &
1461  fnm,fnp,rdnw,msfux,msfuy,msfvx,msfvx_inv,msfvy,msftx,msfty,step,config_flags,ids,ide, &
1462  jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte)
1463
1464  IMPLICIT NONE
1465  TYPE(grid_config_rec_type),INTENT(IN) :: config_flags
1466  INTEGER,INTENT(IN) :: ids,ide,jds,jde,kds,kde
1467  INTEGER,INTENT(IN) :: ims,ime,jms,jme,kms,kme
1468  INTEGER,INTENT(IN) :: its,ite,jts,jte,kts,kte
1469  INTEGER,INTENT(IN) :: step
1470  REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(IN) :: g_u,u,g_v,v,g_u_1,u_1, &
1471  g_v_1,v_1,g_t_1,t_1,g_ft,ft
1472  REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(INOUT) :: g_ww,ww,g_ww_1,ww_1, &
1473  g_t,t,g_t_ave,t_ave,g_uam,uam,g_vam,vam,g_wwam,wwam
1474  REAL,DIMENSION(ims:ime,jms:jme),INTENT(IN) :: g_muu,muu,g_muv,muv,g_mut,mut, &
1475  msfux,msfuy,msfvx,msfvx_inv,msfvy,msftx,msfty,g_mu_tend,mu_tend
1476  REAL,DIMENSION(ims:ime,jms:jme),INTENT(OUT) :: g_muave,muave,g_muts,muts,g_mudf,mudf
1477  REAL,DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: g_mu,mu
1478  REAL,DIMENSION(kms:kme),INTENT(IN) :: fnm,fnp,dnw,rdnw
1479  REAL,INTENT(IN) :: rdx,rdy,dts,epssm
1480  REAL,DIMENSION(its:ite,kts:kte) :: g_wdtn,wdtn,g_dvdxi,dvdxi
1481  REAL,DIMENSION(its:ite) :: g_dmdt,dmdt
1482  INTEGER :: i,j,k,i_start,i_end,j_start,j_end,k_start,k_end
1483  INTEGER :: i_endu,j_endv
1484  REAL :: g_acc,acc
1485
1486  i_start =its
1487
1488  i_end =min(ite,ide-1)
1489
1490  j_start =jts
1491
1492  j_end =min(jte,jde-1)
1493
1494  k_start =kts
1495
1496  k_end =kte-1
1497
1498  IF( .NOT. config_flags%periodic_x ) THEN
1499
1500  IF( config_flags%specified .or. config_flags%nested ) THEN
1501
1502  i_start =max(its,ids+1)
1503
1504  i_end =min(ite,ide-2)
1505  ENDIF
1506  ENDIF
1507
1508  IF( config_flags%specified .or. config_flags%nested ) THEN
1509
1510  j_start =max(jts,jds+1)
1511
1512  j_end =min(jte,jde-2)
1513  ENDIF
1514
1515  i_endu =ite
1516
1517  j_endv =jte
1518
1519  DO j =j_start,j_end
1520  DO i =i_start,i_end
1521
1522  g_dmdt(i) =0.0
1523  dmdt(i) =0.
1524
1525  ENDDO
1526
1527  DO k =k_start,k_end
1528  DO i =i_start,i_end
1529
1530  g_dvdxi(i,k) =msftx(i,j) *msfty(i,j)*(rdy*((g_v(i,k,j+1) +(muv(i,j+1) &
1531 *g_v_1(i,k,j+1) +g_muv(i,j+1)*v_1(i,k,j+1))*msfvx_inv(i,j+1)) -(g_v(i,k,j) &
1532  +(muv(i,j)*g_v_1(i,k,j) +g_muv(i,j)*v_1(i,k,j))*msfvx_inv(i,j))) &
1533  +rdx*((g_u(i+1,k,j) +((muu(i+1,j)*g_u_1(i+1,k,j) +g_muu(i+1,j) &
1534 *u_1(i+1,k,j))/msfuy(i+1,j))) -(g_u(i,k,j) +((muu(i,j)*g_u_1(i,k,j) &
1535  +g_muu(i,j)*u_1(i,k,j))/msfuy(i,j)))))
1536  dvdxi(i,k) =msftx(i,j) *msfty(i,j)*(rdy*((v(i,k,j+1) +muv(i,j+1)*v_1(i,k,j+1) &
1537 *msfvx_inv(i,j+1)) -(v(i,k,j) +muv(i,j)*v_1(i,k,j)*msfvx_inv(i,j))) +rdx*((u(i+1,k,j) &
1538  +muu(i+1,j)*u_1(i+1,k,j)/msfuy(i+1,j)) -(u(i,k,j) +muu(i,j)*u_1(i,k,j)/msfuy(i,j))))
1539
1540  g_dmdt(i) =g_dmdt(i) +dnw(k)*g_dvdxi(i,k)
1541  dmdt(i) =dmdt(i) +dnw(k)*dvdxi(i,k)
1542
1543  ENDDO
1544  ENDDO
1545
1546  DO i =i_start,i_end
1547
1548  g_muave(i,j) =g_mu(i,j)
1549  muave(i,j) =mu(i,j)
1550
1551  g_mu(i,j) =g_mu(i,j) +dts*(g_dmdt(i) +g_mu_tend(i,j))
1552  mu(i,j) =mu(i,j) +dts*(dmdt(i) +mu_tend(i,j))
1553
1554  g_mudf(i,j) =(g_dmdt(i) +g_mu_tend(i,j))
1555  mudf(i,j) =(dmdt(i) +mu_tend(i,j))
1556
1557  g_muts(i,j) =g_mut(i,j) +g_mu(i,j)
1558  muts(i,j) =mut(i,j) +mu(i,j)
1559
1560  g_muave(i,j) =.5*((1.+epssm)*g_mu(i,j) +(1.-epssm)*g_muave(i,j))
1561  muave(i,j) =.5*((1.+epssm)*mu(i,j) +(1.-epssm)*muave(i,j))
1562
1563  ENDDO
1564
1565  DO k =2,k_end
1566  DO i =i_start,i_end
1567
1568  g_ww(i,k,j) =g_ww(i,k-1,j) -(dnw(k-1)*(g_dmdt(i) +g_dvdxi(i,k-1) &
1569  +g_mu_tend(i,j))/msfty(i,j))
1570  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)
1571
1572  ENDDO
1573  ENDDO
1574
1575  DO k =1,k_end
1576  DO i =i_start,i_end
1577
1578  g_ww(i,k,j) =g_ww(i,k,j) -g_ww_1(i,k,j)
1579  ww(i,k,j) =ww(i,k,j) -ww_1(i,k,j)
1580
1581  ENDDO
1582  ENDDO
1583  ENDDO
1584
1585  DO j =j_start,j_end
1586  DO k =1,k_end
1587  DO i =i_start,i_end
1588
1589  g_t_ave(i,k,j) =g_t(i,k,j)
1590  t_ave(i,k,j) =t(i,k,j)
1591
1592  g_t(i,k,j) =g_t(i,k,j) +msfty(i,j) *dts*g_ft(i,k,j)
1593  t(i,k,j) =t(i,k,j) +msfty(i,j) *dts*ft(i,k,j)
1594
1595  ENDDO
1596  ENDDO
1597  ENDDO
1598
1599  DO j =j_start,j_end
1600  DO i =i_start,i_end
1601
1602  g_wdtn(i,1) =0.0
1603  wdtn(i,1) =0.
1604
1605  g_wdtn(i,kde) =0.0
1606  wdtn(i,kde) =0.
1607
1608  ENDDO
1609
1610  DO k =2,k_end
1611  DO i =i_start,i_end
1612
1613  g_wdtn(i,k) =ww(i,k,j)*(fnm(k)*g_t_1(i,k,j) +fnp(k)*g_t_1(i,k-1,j)) &
1614  +g_ww(i,k,j)*(fnm(k)*t_1(i,k,j) +fnp(k)*t_1(i,k-1,j))
1615  wdtn(i,k) =ww(i,k,j)*(fnm(k)*t_1(i,k,j) +fnp(k)*t_1(i,k-1,j))
1616
1617  ENDDO
1618  ENDDO
1619
1620  DO k =1,k_end
1621  DO i =i_start,i_end
1622
1623  g_t(i,k,j) =g_t(i,k,j) -dts *msfty(i,j)*(msftx(i,j)*(.5 *rdy*(v(i,k,j+1) &
1624 *(g_t_1(i,k,j+1) +g_t_1(i,k,j)) +g_v(i,k,j+1)*(t_1(i,k,j+1) +t_1(i,k,j)) &
1625  -(v(i,k,j)*(g_t_1(i,k,j) +g_t_1(i,k,j-1)) +g_v(i,k,j)*(t_1(i,k,j) &
1626  +t_1(i,k,j-1)))) +.5 *rdx*(u(i+1,k,j)*(g_t_1(i+1,k,j) +g_t_1(i,k,j)) &
1627  +g_u(i+1,k,j)*(t_1(i+1,k,j) +t_1(i,k,j)) -(u(i,k,j)*(g_t_1(i,k,j) &
1628  +g_t_1(i-1,k,j)) +g_u(i,k,j)*(t_1(i,k,j) +t_1(i-1,k,j))))) +rdnw(k) &
1629 *(g_wdtn(i,k+1) -g_wdtn(i,k)))
1630  t(i,k,j) =t(i,k,j) -dts *msfty(i,j)*(msftx(i,j)*(.5 *rdy*(v(i,k,j+1)*(t_1(i,k,j+1) &
1631  +t_1(i,k,j)) -v(i,k,j)*(t_1(i,k,j) +t_1(i,k,j-1))) +.5 *rdx*(u(i+1,k,j) &
1632 *(t_1(i+1,k,j) +t_1(i,k,j)) -u(i,k,j)*(t_1(i,k,j) +t_1(i-1,k,j)))) +rdnw(k) &
1633 *(wdtn(i,k+1) -wdtn(i,k)))
1634
1635  ENDDO
1636  ENDDO
1637  ENDDO
1638
1639  END SUBROUTINE g_advance_mu_t
1640
1641 !        Generated by TAPENADE     (INRIA, Tropics team)
1642 !  Tapenade 3.5 (r3785) - 22 Mar 2011 18:35
1643 !
1644 !  Differentiation of advance_w in forward (tangent) mode:
1645 !   variations   of useful results: w t_2ave ph
1646 !   with respect to varying inputs: ww w_save muave gamma u alpha
1647 !                v w rw_tend ph_1 t_2ave cqw muts t_1 t_2 ph alt
1648 !                c2a mut ph_tend a
1649 !   RW status of diff variables: ww:in w_save:in muave:in gamma:in
1650 !                u:in alpha:in v:in w:in-out rw_tend:in ph_1:in
1651 !                t_2ave:in-out cqw:in muts:in t_1:in t_2:in ph:in-out
1652 !                alt:in c2a:in mut:in ph_tend:in a:in
1653 ! domain dims
1654 ! memory dims
1655 SUBROUTINE G_ADVANCE_W(w, wd, rw_tend, rw_tendd, ww, wwd, w_save, &
1656 &  w_saved, u, ud, v, vd, mu1, mut, mutd, muave, muaved, muts, mutsd, &
1657 &  t_2ave, t_2aved, t_2, t_2d, t_1, t_1d, ph, phd, ph_1, ph_1d, phb, &
1658 &  ph_tend, ph_tendd, ht, c2a, c2ad, cqw, cqwd, alt, altd, alb, a, ad, &
1659 &  alpha, alphad, gamma, gammad, rdx, rdy, dts, t0, epssm, dnw, fnm, fnp&
1660 &  , rdnw, rdn, cf1, cf2, cf3, msftx, msfty, config_flags, top_lid, ids, &
1661 &  ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, &
1662 &  jte, kts, kte)
1663   IMPLICIT NONE
1664 ! religion first
1665 ! stuff coming in
1666   TYPE(GRID_CONFIG_REC_TYPE), INTENT(IN) :: config_flags
1667   INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde
1668   INTEGER, INTENT(IN) :: ims, ime, jms, jme, kms, kme
1669   INTEGER, INTENT(IN) :: its, ite, jts, jte, kts, kte
1670   LOGICAL, INTENT(IN) :: top_lid
1671   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: t_2ave, w&
1672 &  , ph
1673   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: t_2aved, &
1674 &  wd, phd
1675   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: rw_tend, ww&
1676 &  , w_save, u, v, t_2, t_1, ph_1, phb, ph_tend, alpha, gamma, a, c2a, &
1677 &  cqw, alb, alt
1678   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: rw_tendd, &
1679 &  wwd, w_saved, ud, vd, t_2d, t_1d, ph_1d, ph_tendd, alphad, gammad, ad&
1680 &  , c2ad, cqwd, altd
1681   REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: mu1, mut, muave, muts&
1682 &  , ht, msftx, msfty
1683   REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: mutd, muaved, mutsd
1684   REAL, DIMENSION(kms:kme), INTENT(IN) :: fnp, fnm, rdnw, rdn, dnw
1685   REAL, INTENT(IN) :: rdx, rdy, dts, cf1, cf2, cf3, t0, epssm
1686 !  Stack based 3d data, tile size.
1687   REAL, DIMENSION(its:ite) :: mut_inv, msft_inv
1688   REAL, DIMENSION(its:ite) :: mut_invd, msft_invd
1689   REAL, DIMENSION(its:ite, kts:kte) :: rhs, wdwn
1690   REAL, DIMENSION(its:ite, kts:kte) :: rhsd, wdwnd
1691   INTEGER :: i, j, k, i_start, i_end, j_start, j_end, k_start, k_end
1692   REAL, DIMENSION(kts:kte) :: dampwt
1693   REAL, DIMENSION(kts:kte) :: dampwtd
1694   REAL :: htop, hbot, hdepth, hk
1695   REAL :: htopd, hbotd, hkd
1696   REAL :: pi, dampmag
1697   REAL :: arg1
1698   REAL :: arg1d
1699   REAL :: arg2
1700   REAL :: arg2d
1701 !<DESCRIPTION>
1702 !
1703 !  advance_w advances the implicit w and geopotential equations.
1704 !
1705 !</DESCRIPTION>
1706 !  set loop limits.
1707 !  Currently set for periodic boundary conditions
1708   i_start = its
1709   IF (ite .GT. ide - 1) THEN
1710     i_end = ide - 1
1711   ELSE
1712     i_end = ite
1713   END IF
1714   j_start = jts
1715   IF (jte .GT. jde - 1) THEN
1716     j_end = jde - 1
1717   ELSE
1718     j_end = jte
1719   END IF
1720   k_start = kts
1721   k_end = kte - 1
1722   IF (.NOT.config_flags%periodic_x) THEN
1723     IF (config_flags%specified .OR. config_flags%nested) THEN
1724       IF (its .LT. ids + 1) THEN
1725         i_start = ids + 1
1726       ELSE
1727         i_start = its
1728       END IF
1729       IF (ite .GT. ide - 2) THEN
1730         i_end = ide - 2
1731       ELSE
1732         i_end = ite
1733       END IF
1734     END IF
1735   END IF
1736   IF (config_flags%specified .OR. config_flags%nested) THEN
1737     IF (jts .LT. jds + 1) THEN
1738       j_start = jds + 1
1739     ELSE
1740       j_start = jts
1741     END IF
1742     IF (jte .GT. jde - 2) THEN
1743       j_end = jde - 2
1744     ELSE
1745       j_end = jte
1746     END IF
1747   END IF
1748   pi = 4.*ATAN(1.)
1749   dampmag = dts*config_flags%dampcoef
1750   hdepth = config_flags%zdamp
1751 ! calculation of phi and w equations
1752 !   Comments on map scale factors:
1753 !   phi equation is:
1754 !    partial d/dt(rho phi/my) = -mx partial d/dx(phi rho u/my)
1755 !                               -mx partial d/dy(phi rho v/mx)
1756 !                               - partial d/dz(phi rho w/my) + rho g w/my
1757 !   as with scalar equation, use uncoupled value (here phi) to find the
1758 !   coupled tendency (rho phi/my)
1759 !   here as usual rho -> ~'mu'
1760 !
1761 !   w eqn [divided by my] is:
1762 !     partial d/dt(rho w/my) = -mx partial d/dx(w rho u/my)
1763 !                              -mx partial d/dy(v rho v/mx)
1764 !                              - partial d/dz(w rho w/my)
1765 !                              +rho[(u*u+v*v)/a + 2 u omega cos(lat) -
1766 !                                   (1/rho) partial dp/dz - g + Fz]/my
1767 !   here as usual rho -> ~'mu'
1768 !
1769 !  'u,v,w' sent here must be coupled variables (= rho w/my etc.) by their usage
1770   DO i=i_start,i_end
1771     rhsd(i, 1) = 0.0
1772     rhs(i, 1) = 0.
1773   END DO
1774   mut_invd = 0.0
1775   wdwnd = 0.0
1776   rhsd = 0.0
1777   dampwtd = 0.0
1778 j_loop_w:DO j=j_start,j_end
1779     DO i=i_start,i_end
1780       mut_invd(i) = -(mutd(i, j)/mut(i, j)**2)
1781       mut_inv(i) = 1./mut(i, j)
1782       msft_invd(i) = 0.0
1783       msft_inv(i) = 1./msfty(i, j)
1784     END DO
1785     DO k=1,k_end
1786       DO i=i_start,i_end
1787         t_2aved(i, k, j) = .5*((1.+epssm)*t_2d(i, k, j)+(1.-epssm)*&
1788 &          t_2aved(i, k, j))
1789         t_2ave(i, k, j) = .5*((1.+epssm)*t_2(i, k, j)+(1.-epssm)*t_2ave(&
1790 &          i, k, j))
1791         t_2aved(i, k, j) = ((t_2aved(i, k, j)+t0*muaved(i, j))*muts(i, j&
1792 &          )*(t0+t_1(i, k, j))-(t_2ave(i, k, j)+muave(i, j)*t0)*(mutsd(i&
1793 &          , j)*(t0+t_1(i, k, j))+muts(i, j)*t_1d(i, k, j)))/(muts(i, j)*&
1794 &          (t0+t_1(i, k, j)))**2
1795         t_2ave(i, k, j) = (t_2ave(i, k, j)+muave(i, j)*t0)/(muts(i, j)*(&
1796 &          t0+t_1(i, k, j)))
1797         wdwnd(i, k+1) = .5*rdnw(k)*((wwd(i, k+1, j)+wwd(i, k, j))*(ph_1(&
1798 &          i, k+1, j)-ph_1(i, k, j)+phb(i, k+1, j)-phb(i, k, j))+(ww(i, k&
1799 &          +1, j)+ww(i, k, j))*(ph_1d(i, k+1, j)-ph_1d(i, k, j)))
1800         wdwn(i, k+1) = .5*(ww(i, k+1, j)+ww(i, k, j))*rdnw(k)*(ph_1(i, k&
1801 &          +1, j)-ph_1(i, k, j)+phb(i, k+1, j)-phb(i, k, j))
1802         rhsd(i, k+1) = dts*(ph_tendd(i, k+1, j)+.5*g*(1.-epssm)*wd(i, k+&
1803 &          1, j))
1804         rhs(i, k+1) = dts*(ph_tend(i, k+1, j)+.5*g*(1.-epssm)*w(i, k+1, &
1805 &          j))
1806       END DO
1807     END DO
1808 !   building up RHS of phi equation
1809 !   ph_tend contains terms 1 and 2; now adding 3 and 4 in stages:
1810 !   here rhs = delta t [ph_tend + ~g*w/2 - ~ww * partial d phi/dz]
1811     DO k=2,k_end
1812       DO i=i_start,i_end
1813         rhsd(i, k) = rhsd(i, k) - dts*(fnm(k)*wdwnd(i, k+1)+fnp(k)*wdwnd&
1814 &          (i, k))
1815         rhs(i, k) = rhs(i, k) - dts*(fnm(k)*wdwn(i, k+1)+fnp(k)*wdwn(i, &
1816 &          k))
1817       END DO
1818     END DO
1819 !  NOTE:  phi'' is not coupled with the map-scale factor  (1/m),
1820 !         but it's tendency is, so must multiply by msft here
1821 !   Comments on map scale factors:
1822 !   building up RHS of phi equation
1823 !   ph_tend contains terms 1 and 2; now adding 3 and 4 in stages:
1824 !   here rhs = ph_previous + (msft/mu)*[rhs_previous - ~ww * delta t *
1825 !                                                      partial d phi/dz]
1826 !            = ph_previous + (msft*delta t/mu)*[ph_tend + ~g*w/2 - ~ww *
1827 !                                                         partial d phi/dz]
1828     DO k=2,k_end+1
1829       DO i=i_start,i_end
1830         rhsd(i, k) = phd(i, k, j) + msfty(i, j)*(rhsd(i, k)*mut_inv(i)+&
1831 &          rhs(i, k)*mut_invd(i))
1832         rhs(i, k) = ph(i, k, j) + msfty(i, j)*rhs(i, k)*mut_inv(i)
1833         IF (top_lid .AND. k .EQ. k_end + 1) THEN
1834           rhsd(i, k) = 0.0
1835           rhs(i, k) = 0.
1836         END IF
1837       END DO
1838     END DO
1839 !  lower boundary condition on w
1840 !   Comments on map scale factors:
1841 !   Chain rule: if Z=Z(X,Y) [true at the surface] then
1842 !      dZ/dt = dZ/dX * dX/dt + dZ/dY * dY/dt, U=dX/dt, V=dY/dt
1843 !   using capitals to denote actual values
1844 !   In mapped values, u=U, v=V, z=Z, 1/dX=mx/dx, 1/dY=my/dy
1845 !      w = dz/dt = mx u dz/dx + my v dz/dy
1846 !   [where dz/dx is just the surface height change between x
1847 !    gridpoints, and dz/dy is the change between y gridpoints]
1848 !   [cf1, cf2 and cf3 do vertical weighting of u or v values nearest
1849 !    the surface]
1850     DO i=i_start,i_end
1851       wd(i, 1, j) = msfty(i, j)*.5*rdy*((ht(i, j+1)-ht(i, j))*(cf1*vd(i&
1852 &        , 1, j+1)+cf2*vd(i, 2, j+1)+cf3*vd(i, 3, j+1))+(ht(i, j)-ht(i, j&
1853 &        -1))*(cf1*vd(i, 1, j)+cf2*vd(i, 2, j)+cf3*vd(i, 3, j))) + msftx(&
1854 &        i, j)*.5*rdx*((ht(i+1, j)-ht(i, j))*(cf1*ud(i+1, 1, j)+cf2*ud(i+&
1855 &        1, 2, j)+cf3*ud(i+1, 3, j))+(ht(i, j)-ht(i-1, j))*(cf1*ud(i, 1, &
1856 &        j)+cf2*ud(i, 2, j)+cf3*ud(i, 3, j)))
1857       w(i, 1, j) = msfty(i, j)*.5*rdy*((ht(i, j+1)-ht(i, j))*(cf1*v(i, 1&
1858 &        , j+1)+cf2*v(i, 2, j+1)+cf3*v(i, 3, j+1))+(ht(i, j)-ht(i, j-1))*&
1859 &        (cf1*v(i, 1, j)+cf2*v(i, 2, j)+cf3*v(i, 3, j))) + msftx(i, j)*.5&
1860 &        *rdx*((ht(i+1, j)-ht(i, j))*(cf1*u(i+1, 1, j)+cf2*u(i+1, 2, j)+&
1861 &        cf3*u(i+1, 3, j))+(ht(i, j)-ht(i-1, j))*(cf1*u(i, 1, j)+cf2*u(i&
1862 &        , 2, j)+cf3*u(i, 3, j)))
1863     END DO
1864 !
1865 ! Jammed 3 doubly nested loops over k/i into 1 for slight improvement
1866 ! in efficiency.  No change in results (bit-for-bit). JM 20040514
1867 ! (left a blank line where the other two k/i-loops were)
1868 !
1869 !   above surface, begin by adding delta t * previous (coupled) w tendency
1870     DO k=2,k_end
1871       DO i=i_start,i_end
1872         wd(i, k, j) = wd(i, k, j) + dts*rw_tendd(i, k, j) + msft_inv(i)*&
1873 &          .5*dts*g*rdn(k)*((cqwd(i, k, j)*mut_inv(i)+cqw(i, k, j)*&
1874 &          mut_invd(i))*(c2a(i, k, j)*rdnw(k)*((1.+epssm)*(rhs(i, k+1)-&
1875 &          rhs(i, k))+(1.-epssm)*(ph(i, k+1, j)-ph(i, k, j)))-c2a(i, k-1&
1876 &          , j)*rdnw(k-1)*((1.+epssm)*(rhs(i, k)-rhs(i, k-1))+(1.-epssm)*&
1877 &          (ph(i, k, j)-ph(i, k-1, j))))+cqw(i, k, j)*mut_inv(i)*(rdnw(k)&
1878 &          *(c2ad(i, k, j)*((1.+epssm)*(rhs(i, k+1)-rhs(i, k))+(1.-epssm)&
1879 &          *(ph(i, k+1, j)-ph(i, k, j)))+c2a(i, k, j)*((1.+epssm)*(rhsd(i&
1880 &          , k+1)-rhsd(i, k))+(1.-epssm)*(phd(i, k+1, j)-phd(i, k, j))))-&
1881 &          rdnw(k-1)*(c2ad(i, k-1, j)*((1.+epssm)*(rhs(i, k)-rhs(i, k-1))&
1882 &          +(1.-epssm)*(ph(i, k, j)-ph(i, k-1, j)))+c2a(i, k-1, j)*((1.+&
1883 &          epssm)*(rhsd(i, k)-rhsd(i, k-1))+(1.-epssm)*(phd(i, k, j)-phd(&
1884 &          i, k-1, j)))))) + dts*g*msft_inv(i)*(rdn(k)*((c2ad(i, k, j)*&
1885 &          alt(i, k, j)+c2a(i, k, j)*altd(i, k, j))*t_2ave(i, k, j)+c2a(i&
1886 &          , k, j)*alt(i, k, j)*t_2aved(i, k, j)-(c2ad(i, k-1, j)*alt(i, &
1887 &          k-1, j)+c2a(i, k-1, j)*altd(i, k-1, j))*t_2ave(i, k-1, j)-c2a(&
1888 &          i, k-1, j)*alt(i, k-1, j)*t_2aved(i, k-1, j))-muaved(i, j))
1889         w(i, k, j) = w(i, k, j) + dts*rw_tend(i, k, j) + msft_inv(i)*cqw&
1890 &          (i, k, j)*(.5*dts*g*mut_inv(i)*rdn(k)*(c2a(i, k, j)*rdnw(k)*((&
1891 &          1.+epssm)*(rhs(i, k+1)-rhs(i, k))+(1.-epssm)*(ph(i, k+1, j)-ph&
1892 &          (i, k, j)))-c2a(i, k-1, j)*rdnw(k-1)*((1.+epssm)*(rhs(i, k)-&
1893 &          rhs(i, k-1))+(1.-epssm)*(ph(i, k, j)-ph(i, k-1, j))))) + dts*g&
1894 &          *msft_inv(i)*(rdn(k)*(c2a(i, k, j)*alt(i, k, j)*t_2ave(i, k, j&
1895 &          )-c2a(i, k-1, j)*alt(i, k-1, j)*t_2ave(i, k-1, j))-muave(i, j)&
1896 &          )
1897       END DO
1898     END DO
1899     k = k_end + 1
1900     DO i=i_start,i_end
1901       wd(i, k, j) = wd(i, k, j) + dts*rw_tendd(i, k, j) + msft_inv(i)*(-&
1902 &        (.5*dts*g*rdnw(k-1)**2*2.*((mut_invd(i)*c2a(i, k-1, j)+mut_inv(i&
1903 &        )*c2ad(i, k-1, j))*((1.+epssm)*(rhs(i, k)-rhs(i, k-1))+(1.-epssm&
1904 &        )*(ph(i, k, j)-ph(i, k-1, j)))+mut_inv(i)*c2a(i, k-1, j)*((1.+&
1905 &        epssm)*(rhsd(i, k)-rhsd(i, k-1))+(1.-epssm)*(phd(i, k, j)-phd(i&
1906 &        , k-1, j)))))-dts*g*(2.*rdnw(k-1)*((c2ad(i, k-1, j)*alt(i, k-1, &
1907 &        j)+c2a(i, k-1, j)*altd(i, k-1, j))*t_2ave(i, k-1, j)+c2a(i, k-1&
1908 &        , j)*alt(i, k-1, j)*t_2aved(i, k-1, j))+muaved(i, j)))
1909       w(i, k, j) = w(i, k, j) + dts*rw_tend(i, k, j) + msft_inv(i)*(-(.5&
1910 &        *dts*g*mut_inv(i)*rdnw(k-1)**2*2.*c2a(i, k-1, j)*((1.+epssm)*(&
1911 &        rhs(i, k)-rhs(i, k-1))+(1.-epssm)*(ph(i, k, j)-ph(i, k-1, j))))-&
1912 &        dts*g*(2.*rdnw(k-1)*c2a(i, k-1, j)*alt(i, k-1, j)*t_2ave(i, k-1&
1913 &        , j)+muave(i, j)))
1914       IF (top_lid) THEN
1915         wd(i, k, j) = 0.0
1916         w(i, k, j) = 0.
1917       END IF
1918     END DO
1919     DO k=2,k_end+1
1920       DO i=i_start,i_end
1921         wd(i, k, j) = (wd(i, k, j)-ad(i, k, j)*w(i, k-1, j)-a(i, k, j)*&
1922 &          wd(i, k-1, j))*alpha(i, k, j) + (w(i, k, j)-a(i, k, j)*w(i, k-&
1923 &          1, j))*alphad(i, k, j)
1924         w(i, k, j) = (w(i, k, j)-a(i, k, j)*w(i, k-1, j))*alpha(i, k, j)
1925       END DO
1926     END DO
1927     DO k=k_end,2,-1
1928       DO i=i_start,i_end
1929         wd(i, k, j) = wd(i, k, j) - gammad(i, k, j)*w(i, k+1, j) - gamma&
1930 &          (i, k, j)*wd(i, k+1, j)
1931         w(i, k, j) = w(i, k, j) - gamma(i, k, j)*w(i, k+1, j)
1932       END DO
1933     END DO
1934     IF (config_flags%damp_opt .EQ. 3) THEN
1935       DO k=k_end+1,2,-1
1936         DO i=i_start,i_end
1937           htopd = ph_1d(i, k_end+1, j)/g
1938           htop = (ph_1(i, k_end+1, j)+phb(i, k_end+1, j))/g
1939           hkd = ph_1d(i, k, j)/g
1940           hk = (ph_1(i, k, j)+phb(i, k, j))/g
1941           hbotd = htopd
1942           hbot = htop - hdepth
1943           dampwtd(k) = 0.0
1944           dampwt(k) = 0.
1945           IF (hk .GE. hbot) THEN
1946             arg1d = 0.5*pi*(hkd-hbotd)/hdepth
1947             arg1 = 0.5*pi*(hk-hbot)/hdepth
1948             arg2d = 0.5*pi*(hkd-hbotd)/hdepth
1949             arg2 = 0.5*pi*(hk-hbot)/hdepth
1950             dampwtd(k) = dampmag*(arg1d*COS(arg1)*SIN(arg2)+SIN(arg1)*&
1951 &              arg2d*COS(arg2))
1952             dampwt(k) = dampmag*SIN(arg1)*SIN(arg2)
1953           END IF
1954           wd(i, k, j) = ((wd(i, k, j)-(dampwtd(k)*mut(i, j)+dampwt(k)*&
1955 &            mutd(i, j))*w_save(i, k, j)-dampwt(k)*mut(i, j)*w_saved(i, k&
1956 &            , j))*(1.+dampwt(k))-(w(i, k, j)-dampwt(k)*mut(i, j)*w_save(&
1957 &            i, k, j))*dampwtd(k))/(1.+dampwt(k))**2
1958           w(i, k, j) = (w(i, k, j)-dampwt(k)*mut(i, j)*w_save(i, k, j))/&
1959 &            (1.+dampwt(k))
1960         END DO
1961       END DO
1962     END IF
1963     DO k=k_end+1,2,-1
1964       DO i=i_start,i_end
1965         phd(i, k, j) = rhsd(i, k) + (msfty(i, j)*.5*dts*g*(1.+epssm)*wd(&
1966 &          i, k, j)*muts(i, j)-msfty(i, j)*.5*dts*g*(1.+epssm)*w(i, k, j)&
1967 &          *mutsd(i, j))/muts(i, j)**2
1968         ph(i, k, j) = rhs(i, k) + msfty(i, j)*.5*dts*g*(1.+epssm)*w(i, k&
1969 &          , j)/muts(i, j)
1970       END DO
1971     END DO
1972   END DO j_loop_w
1973 END SUBROUTINE G_ADVANCE_W
1974
1975  SUBROUTINE g_sumflux(ru,g_ru,rv,g_rv,ww,g_ww,u_lin,g_u_lin,v_lin, &
1976  g_v_lin,ww_lin,g_ww_lin,muu,g_muu,muv,g_muv,ru_m,g_ru_m,rv_m,g_rv_m, &
1977  ww_m,g_ww_m,epssm,msfux,msfuy,msfvx,msfvx_inv,msfvy,iteration,number_of_small_timesteps, &
1978  ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte)
1979
1980  IMPLICIT NONE
1981  INTEGER,INTENT(IN) :: number_of_small_timesteps
1982  INTEGER,INTENT(IN) :: iteration
1983  INTEGER,INTENT(IN) :: ids,ide,jds,jde,kds,kde
1984  INTEGER,INTENT(IN) :: ims,ime,jms,jme,kms,kme
1985  INTEGER,INTENT(IN) :: its,ite,jts,jte,kts,kte
1986  REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(IN) :: g_ru,ru,g_rv,rv,g_ww, &
1987  ww,g_u_lin,u_lin,g_v_lin,v_lin,g_ww_lin,ww_lin
1988  REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(INOUT) :: g_ru_m,ru_m,g_rv_m, &
1989  rv_m,g_ww_m,ww_m
1990  REAL,DIMENSION(ims:ime,jms:jme),INTENT(IN) :: g_muu,muu,g_muv,muv,msfux,msfuy, &
1991  msfvx,msfvy,msfvx_inv
1992  INTEGER :: mini,minj,mink
1993  REAL,INTENT(IN) :: epssm
1994  INTEGER :: i,j,k
1995
1996  IF(iteration == 1 ) THEN
1997
1998  DO j =jts,jte
1999  DO k =kts,kte
2000  DO i =its,ite
2001
2002  g_ru_m(i,k,j) =0.0
2003  ru_m(i,k,j) =0.
2004
2005  g_rv_m(i,k,j) =0.0
2006  rv_m(i,k,j) =0.
2007
2008  g_ww_m(i,k,j) =0.0
2009  ww_m(i,k,j) =0.
2010
2011  ENDDO
2012  ENDDO
2013  ENDDO
2014  ENDIF
2015
2016  mini =min(ide-1,ite)
2017
2018  minj =min(jde-1,jte)
2019
2020  mink =min(kde-1,kte)
2021
2022  DO j =jts,minj
2023  DO k =kts,mink
2024  DO i =its,mini
2025
2026  g_ru_m(i,k,j) =g_ru_m(i,k,j) +g_ru(i,k,j)
2027  ru_m(i,k,j) =ru_m(i,k,j) +ru(i,k,j)
2028
2029  g_rv_m(i,k,j) =g_rv_m(i,k,j) +g_rv(i,k,j)
2030  rv_m(i,k,j) =rv_m(i,k,j) +rv(i,k,j)
2031
2032  g_ww_m(i,k,j) =g_ww_m(i,k,j) +g_ww(i,k,j)
2033  ww_m(i,k,j) =ww_m(i,k,j) +ww(i,k,j)
2034
2035  ENDDO
2036  ENDDO
2037  ENDDO
2038
2039  IF(ite .GT. mini) THEN
2040
2041  DO j =jts,minj
2042  DO k =kts,mink
2043  DO i =mini+1,ite
2044
2045  g_ru_m(i,k,j) =g_ru_m(i,k,j) +g_ru(i,k,j)
2046  ru_m(i,k,j) =ru_m(i,k,j) +ru(i,k,j)
2047
2048  ENDDO
2049  ENDDO
2050  ENDDO
2051  END IF
2052
2053  IF(jte .GT. minj) THEN
2054
2055  DO j =minj+1,jte
2056  DO k =kts,mink
2057  DO i =its,mini
2058
2059  g_rv_m(i,k,j) =g_rv_m(i,k,j) +g_rv(i,k,j)
2060  rv_m(i,k,j) =rv_m(i,k,j) +rv(i,k,j)
2061
2062  ENDDO
2063  ENDDO
2064  ENDDO
2065  END IF
2066
2067  IF( kte .GT. mink) THEN
2068
2069  DO j =jts,minj
2070  DO k =mink+1,kte
2071  DO i =its,mini
2072
2073  g_ww_m(i,k,j) =g_ww_m(i,k,j) +g_ww(i,k,j)
2074  ww_m(i,k,j) =ww_m(i,k,j) +ww(i,k,j)
2075
2076  ENDDO
2077  ENDDO
2078  ENDDO
2079  END IF
2080
2081  IF(iteration == number_of_small_timesteps) THEN
2082
2083  DO j =jts,minj
2084  DO k =kts,mink
2085  DO i =its,mini
2086
2087  g_ru_m(i,k,j) =g_ru_m(i,k,j)/number_of_small_timesteps +((muu(i,j) &
2088 *g_u_lin(i,k,j) +g_muu(i,j)*u_lin(i,k,j))/msfuy(i,j))
2089  ru_m(i,k,j) =ru_m(i,k,j)/number_of_small_timesteps +muu(i,j)*u_lin(i,k,j)/msfuy(i,j)
2090
2091  g_rv_m(i,k,j) =g_rv_m(i,k,j)/number_of_small_timesteps +(muv(i,j) &
2092 *g_v_lin(i,k,j) +g_muv(i,j)*v_lin(i,k,j))*msfvx_inv(i,j)
2093  rv_m(i,k,j) =rv_m(i,k,j)/number_of_small_timesteps +muv(i,j)*v_lin(i,k,j)*msfvx_inv(i,j)
2094
2095  g_ww_m(i,k,j) =g_ww_m(i,k,j)/number_of_small_timesteps +g_ww_lin(i,k,j)
2096  ww_m(i,k,j) =ww_m(i,k,j)/number_of_small_timesteps +ww_lin(i,k,j)
2097
2098  ENDDO
2099  ENDDO
2100  ENDDO
2101
2102  IF(ite .GT. mini) THEN
2103
2104  DO j =jts,minj
2105  DO k =kts,mink
2106  DO i =mini+1,ite
2107
2108  g_ru_m(i,k,j) =g_ru_m(i,k,j)/number_of_small_timesteps +((muu(i,j) &
2109 *g_u_lin(i,k,j) +g_muu(i,j)*u_lin(i,k,j))/msfuy(i,j))
2110  ru_m(i,k,j) =ru_m(i,k,j)/number_of_small_timesteps +muu(i,j)*u_lin(i,k,j)/msfuy(i,j)
2111
2112  ENDDO
2113  ENDDO
2114  ENDDO
2115  END IF
2116
2117  IF(jte .GT. minj) THEN
2118
2119  DO j =minj+1,jte
2120  DO k =kts,mink
2121  DO i =its,mini
2122
2123  g_rv_m(i,k,j) =g_rv_m(i,k,j)/number_of_small_timesteps +(muv(i,j) &
2124 *g_v_lin(i,k,j) +g_muv(i,j)*v_lin(i,k,j))*msfvx_inv(i,j)
2125  rv_m(i,k,j) =rv_m(i,k,j)/number_of_small_timesteps +muv(i,j)*v_lin(i,k,j)*msfvx_inv(i,j)
2126
2127  ENDDO
2128  ENDDO
2129  ENDDO
2130  END IF
2131
2132  IF( kte .GT. mink) THEN
2133
2134  DO j =jts,minj
2135  DO k =mink+1,kte
2136  DO i =its,mini
2137
2138  g_ww_m(i,k,j) =g_ww_m(i,k,j)/number_of_small_timesteps +g_ww_lin(i,k,j)
2139  ww_m(i,k,j) =ww_m(i,k,j)/number_of_small_timesteps +ww_lin(i,k,j)
2140
2141  ENDDO
2142  ENDDO
2143  ENDDO
2144  END IF
2145  ENDIF
2146
2147  END SUBROUTINE g_sumflux
2148
2149  SUBROUTINE g_init_module_small_step
2150
2151  END SUBROUTINE g_init_module_small_step
2152
2153  END MODULE g_module_small_step_em
2154