
! ======================================================================================
! This file was generated by the version 5.3.5 of DFT on 06/19/2010. The differentiation
! transforming system(DFT) was jointly developed and sponsored by LASG of IAP(1998-2010)
! and LSEC of ICMSEC, AMSS(2001-2003)
! The copyright of the DFT system was declared by Walls at LASG, 1998-2010
! ======================================================================================
! Ning Pan, 2010-07-10: Change Diff_ to g_ 
! Ning Pan, 2010-07-11: Change the order of arguments

 MODULE g_module_small_step_em

 USE module_configure
 USE module_model_constants

 CONTAINS

SUBROUTINE g_advance_all (u, g_u, ru_tend, g_ru_tend, v, g_v, rv_tend, g_rv_tend,     &
                          w, g_w, rw_tend, g_rw_tend, t, g_t, t_tend, g_t_tend,       &
                          mu, g_mu, mu_tend, g_mu_tend, ph, g_ph, ph_tend, g_ph_tend, &
                          muu, g_muu, muv, g_muv, mut, g_mut,                         &
                          msfuy, msfvx, msfty,           &
                          dts,                           &
                          config_flags, spec_zone,       &
                          ids, ide, jds, jde, kds, kde,  &
                          ims, ime, jms, jme, kms, kme,  &
                          ips, ipe, jps, jpe, kps, kpe,  &
                          its, ite, jts, jte, kts, kte  )

      IMPLICIT NONE  ! religion first

      TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags

      INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
      INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
      INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
      INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
      INTEGER,      INTENT(IN   )    :: spec_zone

      REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),  &
            INTENT(INOUT) ::                          &
                                                  g_u,  &
                                                  g_v,  &
                                                  g_w,  &
                                                  g_t,  &
                                                  g_ph
      REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),  &
            INTENT(INOUT) ::                          &
                                                  u,  &
                                                  v,  &
                                                  w,  &
                                                  t,  &
                                                  ph

      REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
            INTENT(IN   ) ::                          &
                                             g_ru_tend, &
                                             g_rv_tend, &
                                             g_rw_tend, &
                                              g_t_tend, &
                                             g_ph_tend
      REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
            INTENT(IN   ) ::                          &
                                             ru_tend, &
                                             rv_tend, &
                                             rw_tend, &
                                              t_tend, &
                                             ph_tend

      REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(INOUT) :: g_mu
      REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(INOUT) :: mu
      REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: g_mu_tend
      REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: mu_tend

      REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: g_muu,  &
                                                                g_muv,  &
                                                                g_mut
      REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: muu,  &
                                                                muv,  &
                                                                mut

      REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: msfuy,  &
                                                                msfvx,  &
                                                                msfty

      REAL,                                    INTENT(IN   ) :: dts
    

      INTEGER :: i,j,k, i_start, i_end, j_start, j_end, k_start, k_end
      INTEGER :: i_endu, j_endv, k_endw

      INTEGER :: i_start_u_tend, i_end_u_tend, j_start_v_tend, j_end_v_tend

!<DESCRIPTION>
!
!  advance_all advances the explicit perturbation horizontal momentum 
!  equations (u,v) by adding in the large-timestep tendency along with 
!  the small timestep pressure gradient tendency.
!
!</DESCRIPTION>

!  now, the real work.
!  set the loop bounds taking into account boundary conditions.

   IF( config_flags%nested .or. config_flags%specified ) THEN
      i_start = max( its,ids+spec_zone )
      i_end   = min( ite,ide-spec_zone-1 )
      j_start = max( jts,jds+spec_zone )
      j_end   = min( jte,jde-spec_zone-1 )
      k_start = kts
      k_end   = min( kte, kde-1 )

      i_endu = min( ite,ide-spec_zone )
      j_endv = min( jte,jde-spec_zone )
      k_endw = k_end

      IF( config_flags%periodic_x) THEN
        i_start = its
        i_end   = min(ite,ide-1)
        i_endu = ite
      ENDIF
   ELSE
      i_start = its
      i_end   = min(ite,ide-1)
      j_start = jts
      j_end   = min(jte,jde-1)
      k_start = kts
      k_end   = kte-1

      i_endu = ite
      j_endv = jte
      k_endw = k_end
   ENDIF

   i_start_u_tend = i_start
   i_end_u_tend   = i_endu
   j_start_v_tend = j_start
   j_end_v_tend   = j_endv

   IF ( config_flags%symmetric_xs .and. (its == ids) ) &
                  i_start_u_tend = i_start_u_tend+1
   IF ( config_flags%symmetric_xe .and. (ite == ide) ) &
                  i_end_u_tend = i_end_u_tend-1
   IF ( config_flags%symmetric_ys .and. (jts == jds) ) &
                  j_start_v_tend = j_start_v_tend+1
   IF ( config_flags%symmetric_ye .and. (jte == jde) ) &
                  j_end_v_tend = j_end_v_tend-1

   u_outer_j_loop: DO j = j_start, j_end
      DO k = k_start, k_end
      DO i = i_start_u_tend, i_end_u_tend
        g_u(i,k,j) = g_u(i,k,j) &
                   + dts*msfuy(i,j)*( g_ru_tend(i,k,j)/muu(i,j) - &
                                      ru_tend(i,k,j)*g_muu(i,j)/(muu(i,j)*muu(i,j)) )
        u(i,k,j) = u(i,k,j) + dts*ru_tend(i,k,j)*msfuy(i,j)/muu(i,j)
      ENDDO
      ENDDO
   ENDDO u_outer_j_loop

   v_outer_j_loop: DO j = j_start_v_tend, j_end_v_tend
      DO k = k_start, k_end
      DO i = i_start, i_end
        g_v(i,k,j) = g_v(i,k,j) &
                   + dts*msfvx(i,j)*( g_rv_tend(i,k,j)/muv(i,j) - &
                                      rv_tend(i,k,j)*g_muv(i,j)/(muv(i,j)*muv(i,j)) )
        v(i,k,j) = v(i,k,j) + dts*rv_tend(i,k,j)*msfvx(i,j)/muv(i,j)
      ENDDO
      ENDDO
   ENDDO v_outer_j_loop

   i_start = its
   i_end   = min(ite,ide-1)
   j_start = jts
   j_end   = min(jte,jde-1)
   k_start = kts
   k_end   = kte-1
   IF ( .NOT. config_flags%periodic_x )THEN
     IF ( config_flags%specified .or. config_flags%nested ) then
       i_start = max(its,ids+1)
       i_end   = min(ite,ide-2)
     ENDIF
   ENDIF
   IF ( config_flags%specified .or. config_flags%nested ) then
     j_start = max(jts,jds+1)
     j_end   = min(jte,jde-2)
   ENDIF

   DO j = j_start, j_end
   DO i=i_start, i_end
     g_mu(i,j) = g_mu(i,j)+dts*g_mu_tend(i,j)
     mu(i,j) = mu(i,j)+dts*mu_tend(i,j)
   ENDDO
   ENDDO

   DO j=j_start, j_end
     DO k=1,k_end
     DO i=i_start, i_end
       g_t(i,k,j) = g_t(i,k,j) &
                  + dts*msfty(i,j)*( g_t_tend(i,k,j)/mut(i,j) - &
                                     t_tend(i,k,j)*g_mut(i,j)/(mut(i,j)*mut(i,j)) )
       t(i,k,j) = t(i,k,j) + msfty(i,j)*dts*t_tend(i,k,j)/mut(i,j)
     END DO
     END DO
   ENDDO   

   i_start = its
   i_end   = min(ite,ide-1)
   j_start = jts
   j_end   = min(jte,jde-1)
   k_start = kts
   k_end   = kte-1
   IF ( .NOT. config_flags%periodic_x )THEN
      IF ( config_flags%specified .or. config_flags%nested ) then
        i_start = max(its,ids+1)
        i_end   = min(ite,ide-2)
      ENDIF
   ENDIF
   IF ( config_flags%specified .or. config_flags%nested ) then
      j_start = max(jts,jds+1)
      j_end   = min(jte,jde-2)
   ENDIF

IF ( config_flags%non_hydrostatic ) THEN
  j_loop_w:  DO j = j_start, j_end
    DO k=2,k_end+1
      DO i=i_start, i_end
        g_w(i,k,j) = g_w(i,k,j) &
                   + dts*msfty(i,j)*( g_rw_tend(i,k,j)/mut(i,j) - &
                                      rw_tend(i,k,j)*g_mut(i,j)/(mut(i,j)*mut(i,j)) )
        w(i,k,j)=w(i,k,j)+dts*rw_tend(i,k,j)*msfty(i,j)/mut(i,j)

        g_ph(i,k,j) = g_ph(i,k,j) &
                    + dts*msfty(i,j)*( g_ph_tend(i,k,j)/mut(i,j) - &
                                       ph_tend(i,k,j)*g_mut(i,j)/(mut(i,j)*mut(i,j)) )
        ph(i,k,j) = ph(i,k,j)+dts*ph_tend(i,k,j)*msfty(i,j)/mut(i,j)
      ENDDO
    ENDDO
  ENDDO j_loop_w
ENDIF

END SUBROUTINE g_advance_all


SUBROUTINE g_save_ph_mu ( ph_1, g_ph_1, ph_2, g_ph_2, ph_save, g_ph_save, &
                          mu_1, g_mu_1, mu_2, g_mu_2, mu_save, g_mu_save, &
                          rk_step,                      &
                          ids,ide, jds,jde, kds,kde,    & 
                          ims,ime, jms,jme, kms,kme,    &
                          its,ite, jts,jte, kts,kte    )

  IMPLICIT NONE  ! religion first

! declarations for the stuff coming in

  INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
  INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
  INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte

  INTEGER,      INTENT(IN   )    :: rk_step

  REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: g_ph_1
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: ph_1

  REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(  OUT) :: g_ph_save
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(  OUT) :: ph_save

  REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: g_ph_2
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: ph_2

  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: g_mu_1,g_mu_2
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mu_1,mu_2
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(  OUT) :: g_mu_save
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(  OUT) :: mu_save

! local variables

  INTEGER :: i, j, k
  INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end


    i_start = its
    i_end   = min(ite,ide-1)
    j_start = jts
    j_end   = min(jte,jde-1)
    k_start = kts
    k_end = min(kte,kde-1)

    !  if this is the first RK step, reset *_1 to *_2
    !  (we are replacing the t-dt fields with the time t fields)

    IF ((rk_step == 1) ) THEN

      DO j=j_start, j_end
      DO k=k_start, min(kde,kte)
      DO i=i_start, i_end
        g_ph_1(i,k,j) = g_ph_2(i,k,j)
        ph_1(i,k,j) = ph_2(i,k,j)
      ENDDO
      ENDDO
      ENDDO

      DO j=j_start, j_end
      DO i=i_start, i_end
        g_mu_1(i,j) = g_mu_2(i,j)
        mu_1(i,j) = mu_2(i,j)
      ENDDO
      ENDDO

    END IF

    ! set up the small timestep variables

    DO j=j_start, j_end
!    DO k=k_start, min(kde,kte)
    DO k=k_start, kde
    DO i=i_start, i_end
      g_ph_save(i,k,j) = g_ph_2(i,k,j)
      ph_save(i,k,j) = ph_2(i,k,j)
      g_ph_2(i,k,j) = g_ph_1(i,k,j)-g_ph_2(i,k,j)
      ph_2(i,k,j) = ph_1(i,k,j)-ph_2(i,k,j)
    ENDDO
    ENDDO
    ENDDO

    DO j=j_start, j_end
    DO i=i_start, i_end
      g_mu_save(i,j)=g_mu_2(i,j)
      mu_save(i,j)=mu_2(i,j)
      g_mu_2(i,j)=g_mu_1(i,j)-g_mu_2(i,j)
      mu_2(i,j)=mu_1(i,j)-mu_2(i,j)
    ENDDO
    ENDDO

END SUBROUTINE g_save_ph_mu

!----------------------------------------------------------------------

SUBROUTINE g_restore_ph_mu ( ph_2, g_ph_2, ph_save, g_ph_save, &
                             mu_2, g_mu_2, mu_save, g_mu_save, &
                             ids,ide, jds,jde, kds,kde,    & 
                             ims,ime, jms,jme, kms,kme,    &
                             its,ite, jts,jte, kts,kte )

  IMPLICIT NONE  ! religion first

! declarations for the stuff coming in

  INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
  INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
  INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte

  REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(IN   ) :: g_ph_save
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(IN   ) :: ph_save
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: g_ph_2
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: ph_2
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: g_mu_2
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mu_2
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN   ) :: g_mu_save
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN   ) :: mu_save

! local variables

  INTEGER :: i, j, k
  INTEGER :: i_start, i_end, j_start, j_end


  i_start = its
  i_end   = min(ite,ide-1)
  j_start = jts
  j_end   = min(jte,jde-1)

  DO j = j_start, j_end
  DO k = kds, kde
  DO i = i_start, i_end
    g_ph_2(i,k,j) = g_ph_2(i,k,j) + g_ph_save(i,k,j)
    ph_2(i,k,j) = ph_2(i,k,j) + ph_save(i,k,j)
  ENDDO
  ENDDO
  ENDDO

  DO j = j_start, j_end
  DO i = i_start, i_end
    g_mu_2(i,j) = g_mu_2(i,j) + g_mu_save(i,j)
    mu_2(i,j) = mu_2(i,j) + mu_save(i,j)
  ENDDO
  ENDDO

END SUBROUTINE g_restore_ph_mu

 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, &
 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, &
 mub,mu_1,g_mu_1,mu_2,g_mu_2,muu,g_muu,muus,g_muus,muv,g_muv,muvs, &
 g_muvs,mut,g_mut,muts,g_muts,mudf,g_mudf,u_save,g_u_save,v_save, &
 g_v_save,w_save,g_w_save,t_save,g_t_save,ph_save,g_ph_save,mu_save, &
 g_mu_save,ww,g_ww,ww_save,g_ww_save,c2a,g_c2a,pb,p,g_p,alt,g_alt, &
 msfux,msfuy,msfvx,msfvx_inv,msfvy,msftx,msfty,rdx,rdy,rk_step,ids,ide,jds,jde,kds, &
 kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte)

 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ids,ide,jds,jde,kds,kde
 INTEGER,INTENT(IN) :: ims,ime,jms,jme,kms,kme
 INTEGER,INTENT(IN) :: its,ite,jts,jte,kts,kte
 INTEGER,INTENT(IN) :: rk_step
 REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(INOUT) :: g_u_1,u_1,g_v_1,v_1, &
 g_w_1,w_1,g_t_1,t_1,g_ph_1,ph_1
 REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(OUT) :: g_u_save,u_save, &
 g_v_save,v_save,g_w_save,w_save,g_t_save,t_save,g_ph_save,ph_save
 REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(INOUT) :: g_u_2,u_2,g_v_2,v_2, &
 g_w_2,w_2,g_t_2,t_2,g_ph_2,ph_2
 REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(OUT) :: g_c2a,c2a,g_ww_save,ww_save
 REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(IN) :: pb,g_p,p,g_alt,alt,g_ww,ww
 REAL,DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: g_mu_1,mu_1,g_mu_2,mu_2
 REAL,DIMENSION(ims:ime,jms:jme),INTENT(IN) :: mub,g_muu,muu,g_muv,muv,g_mut, &
 mut,msfux,msfuy,msfvx,msfvx_inv,msfvy,msftx,msfty
 REAL,DIMENSION(ims:ime,jms:jme),INTENT(OUT) :: g_muus,muus,g_muvs,muvs, &
 g_muts,muts,g_mudf,mudf
 REAL,DIMENSION(ims:ime,jms:jme),INTENT(OUT) :: g_mu_save,mu_save
 REAL,INTENT(IN) :: rdx,rdy
 INTEGER :: i,j,k
 INTEGER :: i_start,i_end,j_start,j_end,k_start,k_end
 INTEGER :: i_endu,j_endv

 i_start =its

 i_end =min(ite,ide-1)

 j_start =jts

 j_end =min(jte,jde-1)

 k_start =kts

 k_end =min(kte,kde-1)

 i_endu =ite

 j_endv =jte

 IF((rk_step == 1) ) THEN

 DO j =j_start,j_end
 DO i =i_start,i_end

 g_mu_1(i,j) =g_mu_2(i,j)
 mu_1(i,j) =mu_2(i,j)

 g_ww_save(i,kde,j) =0.0
 ww_save(i,kde,j) =0.

 g_ww_save(i,1,j) =0.0
 ww_save(i,1,j) =0.

 g_mudf(i,j) =0.0
 mudf(i,j) =0.

 ENDDO
 ENDDO

 DO j =j_start,j_end
 DO k =k_start,k_end
 DO i =i_start,i_endu

 g_u_1(i,k,j) =g_u_2(i,k,j)
 u_1(i,k,j) =u_2(i,k,j)

 ENDDO
 ENDDO
 ENDDO

 DO j =j_start,j_endv
 DO k =k_start,k_end
 DO i =i_start,i_end

 g_v_1(i,k,j) =g_v_2(i,k,j)
 v_1(i,k,j) =v_2(i,k,j)

 ENDDO
 ENDDO
 ENDDO

 DO j =j_start,j_end
 DO k =k_start,k_end
 DO i =i_start,i_end

 g_t_1(i,k,j) =g_t_2(i,k,j)
 t_1(i,k,j) =t_2(i,k,j)

 ENDDO
 ENDDO
 ENDDO

 DO j =j_start,j_end
 DO k =k_start,min(kde,kte)
 DO i =i_start,i_end

 g_w_1(i,k,j) =g_w_2(i,k,j)
 w_1(i,k,j) =w_2(i,k,j)

 g_ph_1(i,k,j) =g_ph_2(i,k,j)
 ph_1(i,k,j) =ph_2(i,k,j)

 ENDDO
 ENDDO
 ENDDO

 DO j =j_start,j_end
 DO i =i_start,i_end

 g_muts(i,j) =g_mu_2(i,j)
 muts(i,j) =mub(i,j) +mu_2(i,j)

 ENDDO

 DO i =i_start,i_endu

 g_muus(i,j) =g_muu(i,j)
 muus(i,j) =muu(i,j)

 ENDDO
 ENDDO

 DO j =j_start,j_endv
 DO i =i_start,i_end

 g_muvs(i,j) =g_muv(i,j)
 muvs(i,j) =muv(i,j)

 ENDDO
 ENDDO

 DO j =j_start,j_end
 DO i =i_start,i_end

 g_mu_save(i,j) =g_mu_2(i,j)
 mu_save(i,j) =mu_2(i,j)

!g_mu_2(i,j) =g_mu_2(i,j) -g_mu_2(i,j)
!mu_2(i,j) =mu_2(i,j) -mu_2(i,j)

 g_mu_2(i,j) =0.0  !REVISED BY WALLS
 mu_2(i,j) =0.0  !REVISED BY WALLS

 ENDDO
 ENDDO
 ELSE

 DO j =j_start,j_end
 DO i =i_start,i_end

 g_muts(i,j) =g_mu_1(i,j)
 muts(i,j) =mub(i,j) +mu_1(i,j)

 ENDDO

 DO i =i_start,i_endu

 g_muus(i,j) =0.5*(g_mu_1(i,j) +g_mu_1(i-1,j))
 muus(i,j) =0.5*(mub(i,j) +mu_1(i,j) +mub(i-1,j) +mu_1(i-1,j))

 ENDDO
 ENDDO

 DO j =j_start,j_endv
 DO i =i_start,i_end

 g_muvs(i,j) =0.5*(g_mu_1(i,j) +g_mu_1(i,j-1))
 muvs(i,j) =0.5*(mub(i,j) +mu_1(i,j) +mub(i,j-1) +mu_1(i,j-1))

 ENDDO
 ENDDO

 DO j =j_start,j_end
 DO i =i_start,i_end

 g_mu_save(i,j) =g_mu_2(i,j)
 mu_save(i,j) =mu_2(i,j)

 g_mu_2(i,j) =g_mu_1(i,j) -g_mu_2(i,j)
 mu_2(i,j) =mu_1(i,j) -mu_2(i,j)

 ENDDO
 ENDDO
 END IF

 DO j =j_start,j_end
 DO i =i_start,i_end

 g_ww_save(i,kde,j) =0.0
 ww_save(i,kde,j) =0.

 g_ww_save(i,1,j) =0.0
 ww_save(i,1,j) =0.

 ENDDO
 ENDDO

 DO j =j_start,j_end
 DO k =k_start,k_end
 DO i =i_start,i_end

 g_c2a(i,k,j) =(cpovcv*(g_p(i,k,j))*alt(i,k,j) -g_alt(i,k,j)*cpovcv*(pb(i,k, &
 j) +p(i,k,j)))/(alt(i,k,j)*alt(i,k,j))
 c2a(i,k,j) =cpovcv*(pb(i,k,j) +p(i,k,j))/alt(i,k,j)

 ENDDO
 ENDDO
 ENDDO

 DO j =j_start,j_end
 DO k =k_start,k_end
 DO i =i_start,i_endu

 g_u_save(i,k,j) =g_u_2(i,k,j)
 u_save(i,k,j) =u_2(i,k,j)

 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) &
*g_u_2(i,k,j) +g_muu(i,j)*u_2(i,k,j)))/msfuy(i,j)
 u_2(i,k,j) =(muus(i,j)*u_1(i,k,j) -muu(i,j)*u_2(i,k,j))/msfuy(i,j)

 ENDDO
 ENDDO
 ENDDO

 DO j =j_start,j_endv
 DO k =k_start,k_end
 DO i =i_start,i_end

 g_v_save(i,k,j) =g_v_2(i,k,j)
 v_save(i,k,j) =v_2(i,k,j)

 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) &
*g_v_2(i,k,j) +g_muv(i,j)*v_2(i,k,j)))*msfvx_inv(i,j)
 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)

 ENDDO
 ENDDO
 ENDDO

!LPB[6]
 DO j =j_start,j_end
 DO k =k_start,k_end
 DO i =i_start,i_end

 g_t_save(i,k,j) =g_t_2(i,k,j)
 t_save(i,k,j) =t_2(i,k,j)

 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) &
*g_t_2(i,k,j) +g_mut(i,j)*t_2(i,k,j))
 t_2(i,k,j) =muts(i,j)*t_1(i,k,j) -mut(i,j)*t_2(i,k,j)

 ENDDO
 ENDDO
 ENDDO

!LPB[7]
 DO j =j_start,j_end
 DO k =k_start,kde
 DO i =i_start,i_end

 g_w_save(i,k,j) =g_w_2(i,k,j)
 w_save(i,k,j) =w_2(i,k,j)

 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) &
*g_w_2(i,k,j) +g_mut(i,j)*w_2(i,k,j)))/msfty(i,j)
 w_2(i,k,j) =(muts(i,j)*w_1(i,k,j) -mut(i,j)*w_2(i,k,j))/msfty(i,j)

 g_ph_save(i,k,j) =g_ph_2(i,k,j)
 ph_save(i,k,j) =ph_2(i,k,j)

 g_ph_2(i,k,j) =g_ph_1(i,k,j) -g_ph_2(i,k,j)
 ph_2(i,k,j) =ph_1(i,k,j) -ph_2(i,k,j)

 ENDDO
 ENDDO
 ENDDO

 DO j =j_start,j_end
 DO k =k_start,kde
 DO i =i_start,i_end

 g_ww_save(i,k,j) =g_ww(i,k,j)
 ww_save(i,k,j) =ww(i,k,j)

 ENDDO
 ENDDO
 ENDDO

 END SUBROUTINE g_small_step_prep

!        Generated by TAPENADE     (INRIA, Tropics team)
!  Tapenade 3.5 (r3785) - 22 Mar 2011 18:35
!
!  Differentiation of small_step_finish in forward (tangent) mode:
!   variations   of useful results: ww ph_2 u_2 mu_2 w_2 t_2 v_2
!   with respect to varying inputs: u_save ph_save ww w_save mu_save
!                muvs ph_2 u_2 h_diabatic mu_2 w_2 t_save v_save
!                muts t_2 mut muu v_2 muv ww1 muus
!   RW status of diff variables: u_save:in ph_save:in ww:in-out
!                w_save:in mu_save:in muvs:in ph_2:in-out u_2:in-out
!                h_diabatic:in mu_2:in-out w_2:in-out t_save:in
!                v_save:in muts:in t_2:in-out mut:in muu:in v_2:in-out
!                muv:in ww1:in muus:in
SUBROUTINE G_SMALL_STEP_FINISH(u_2, u_2d, u_1, v_2, v_2d, v_1, w_2, w_2d&
&  , w_1, t_2, t_2d, t_1, ph_2, ph_2d, ph_1, ww, wwd, ww1, ww1d, mu_2, &
&  mu_2d, mu_1, mut, mutd, muts, mutsd, muu, muud, muus, muusd, muv, muvd&
&  , muvs, muvsd, u_save, u_saved, v_save, v_saved, w_save, w_saved, &
&  t_save, t_saved, ph_save, ph_saved, mu_save, mu_saved, msfux, msfuy, &
&  msfvx, msfvy, msftx, msfty, h_diabatic, h_diabaticd, &
&  number_of_small_timesteps, dts, rk_step, rk_order, ids, ide, jds, jde&
&  , kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte&
&)
  IMPLICIT NONE
! religion first
!  stuff passed in
  INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde
  INTEGER, INTENT(IN) :: ims, ime, jms, jme, kms, kme
  INTEGER, INTENT(IN) :: its, ite, jts, jte, kts, kte
  INTEGER, INTENT(IN) :: number_of_small_timesteps
  INTEGER, INTENT(IN) :: rk_step, rk_order
  REAL, INTENT(IN) :: dts
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u_1, v_1, &
&  w_1, t_1, ww1, ph_1
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: ww1d
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: u_2, v_2&
&  , w_2, t_2, ww, ph_2
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: u_2d, &
&  v_2d, w_2d, t_2d, wwd, ph_2d
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u_save, &
&  v_save, w_save, t_save, ph_save, h_diabatic
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u_saved, &
&  v_saved, w_saved, t_saved, ph_saved, h_diabaticd
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: muus, muvs
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: muusd, muvsd
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mu_2, mu_1
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mu_2d
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mut, muts, muu, &
&  muv, mu_save
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mutd, mutsd, muud&
&  , muvd, mu_saved
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: msfux, msfuy, msfvx, &
&  msfvy, msftx, msfty
! local stuff
  INTEGER :: i, j, k
  INTEGER :: i_start, i_end, j_start, j_end, i_endu, j_endv
  INTRINSIC MIN
!<DESCRIPTION>
!
!  small_step_finish reconstructs the full uncoupled prognostic variables
!  from the coupled perturbation variables used in the small timesteps.
!
!</DESCRIPTION>
  i_start = its
  IF (ite .GT. ide - 1) THEN
    i_end = ide - 1
  ELSE
    i_end = ite
  END IF
  j_start = jts
  IF (jte .GT. jde - 1) THEN
    j_end = jde - 1
  ELSE
    j_end = jte
  END IF
  i_endu = ite
  j_endv = jte
!    addition of time level t back into variables
  DO j=j_start,j_endv
    DO k=kds,kde-1
      DO i=i_start,i_end
! v coupled with mx
        v_2d(i, k, j) = ((msfvx(i, j)*v_2d(i, k, j)+v_saved(i, k, j)*muv&
&          (i, j)+v_save(i, k, j)*muvd(i, j))*muvs(i, j)-(msfvx(i, j)*v_2&
&          (i, k, j)+v_save(i, k, j)*muv(i, j))*muvsd(i, j))/muvs(i, j)**&
&          2
        v_2(i, k, j) = (msfvx(i, j)*v_2(i, k, j)+v_save(i, k, j)*muv(i, &
&          j))/muvs(i, j)
      END DO
    END DO
  END DO
  DO j=j_start,j_end
    DO k=kds,kde-1
      DO i=i_start,i_endu
! u coupled with my
        u_2d(i, k, j) = ((msfuy(i, j)*u_2d(i, k, j)+u_saved(i, k, j)*muu&
&          (i, j)+u_save(i, k, j)*muud(i, j))*muus(i, j)-(msfuy(i, j)*u_2&
&          (i, k, j)+u_save(i, k, j)*muu(i, j))*muusd(i, j))/muus(i, j)**&
&          2
        u_2(i, k, j) = (msfuy(i, j)*u_2(i, k, j)+u_save(i, k, j)*muu(i, &
&          j))/muus(i, j)
      END DO
    END DO
  END DO
  DO j=j_start,j_end
    DO k=kds,kde
      DO i=i_start,i_end
! w coupled with my
        w_2d(i, k, j) = ((msfty(i, j)*w_2d(i, k, j)+w_saved(i, k, j)*mut&
&          (i, j)+w_save(i, k, j)*mutd(i, j))*muts(i, j)-(msfty(i, j)*w_2&
&          (i, k, j)+w_save(i, k, j)*mut(i, j))*mutsd(i, j))/muts(i, j)**&
&          2
        w_2(i, k, j) = (msfty(i, j)*w_2(i, k, j)+w_save(i, k, j)*mut(i, &
&          j))/muts(i, j)
        ph_2d(i, k, j) = ph_2d(i, k, j) + ph_saved(i, k, j)
        ph_2(i, k, j) = ph_2(i, k, j) + ph_save(i, k, j)
        wwd(i, k, j) = wwd(i, k, j) + ww1d(i, k, j)
        ww(i, k, j) = ww(i, k, j) + ww1(i, k, j)
      END DO
    END DO
  END DO
  IF (rk_step .LT. rk_order) THEN
    DO j=j_start,j_end
      DO k=kds,kde-1
        DO i=i_start,i_end
          t_2d(i, k, j) = ((t_2d(i, k, j)+t_saved(i, k, j)*mut(i, j)+&
&            t_save(i, k, j)*mutd(i, j))*muts(i, j)-(t_2(i, k, j)+t_save(&
&            i, k, j)*mut(i, j))*mutsd(i, j))/muts(i, j)**2
          t_2(i, k, j) = (t_2(i, k, j)+t_save(i, k, j)*mut(i, j))/muts(i&
&            , j)
        END DO
      END DO
    END DO
  ELSE
    DO j=j_start,j_end
      DO k=kds,kde-1
        DO i=i_start,i_end
          t_2d(i, k, j) = ((t_2d(i, k, j)-dts*number_of_small_timesteps*&
&            (mutd(i, j)*h_diabatic(i, k, j)+mut(i, j)*h_diabaticd(i, k, &
&            j))+t_saved(i, k, j)*mut(i, j)+t_save(i, k, j)*mutd(i, j))*&
&            muts(i, j)-(t_2(i, k, j)-dts*number_of_small_timesteps*mut(i&
&            , j)*h_diabatic(i, k, j)+t_save(i, k, j)*mut(i, j))*mutsd(i&
&            , j))/muts(i, j)**2
          t_2(i, k, j) = (t_2(i, k, j)-dts*number_of_small_timesteps*mut&
&            (i, j)*h_diabatic(i, k, j)+t_save(i, k, j)*mut(i, j))/muts(i&
&            , j)
        END DO
      END DO
    END DO
  END IF
  DO j=j_start,j_end
    DO i=i_start,i_end
      mu_2d(i, j) = mu_2d(i, j) + mu_saved(i, j)
      mu_2(i, j) = mu_2(i, j) + mu_save(i, j)
    END DO
  END DO
END SUBROUTINE G_SMALL_STEP_FINISH

 SUBROUTINE g_calc_p_rho(al,g_al,p,g_p,ph,g_ph,alt,g_alt,t_2,g_t_2, &
 t_1,g_t_1,c2a,g_c2a,pm1,g_pm1,mu,g_mu,muts,g_muts,znu,t0,rdnw,dnw, &
 smdiv,non_hydrostatic,step,ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite, &
 jts,jte,kts,kte)

 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ids,ide,jds,jde,kds,kde
 INTEGER,INTENT(IN) :: ims,ime,jms,jme,kms,kme
 INTEGER,INTENT(IN) :: its,ite,jts,jte,kts,kte
 INTEGER,INTENT(IN) :: step
 REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(OUT) :: g_al,al,g_p,p
 REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(IN) :: g_alt,alt,g_t_2,t_2, &
 g_t_1,t_1,g_c2a,c2a
 REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(INOUT) :: g_ph,ph,g_pm1,pm1
 REAL,DIMENSION(ims:ime,jms:jme),INTENT(IN) :: g_mu,mu,g_muts,muts
 REAL,DIMENSION(kms:kme),INTENT(IN) :: dnw,rdnw,znu
 REAL,INTENT(IN) :: t0,smdiv
 LOGICAL,INTENT(IN) :: non_hydrostatic
 INTEGER :: i,j,k
 INTEGER :: i_start,i_end,j_start,j_end,k_start,k_end
 REAL :: g_ptmp,ptmp

 i_start =its

 i_end =min(ite,ide-1)

 j_start =jts

 j_end =min(jte,jde-1)

 k_start =kts

 k_end =min(kte,kde-1)

 IF(non_hydrostatic) THEN

 DO j =j_start,j_end
 DO k =k_start,k_end
 DO i =i_start,i_end

 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) &
 +rdnw(k)*(g_ph(i,k+1,j) -g_ph(i,k,j))) +(-(-1.)*g_muts(i,j)/(muts(i,j) &
*muts(i,j)))*(alt(i,k,j)*mu(i,j) +rdnw(k)*(ph(i,k+1,j) -ph(i,k,j)))
 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)))

 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) &
 +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))) &
*(muts(i,j)*(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)))*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)))*(muts(i,j)*(t0 +t_1(i,k,j)))) -g_al(i,k,j)) +g_c2a(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))) -al(i,k,j))
 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) &
*(t0 +t_1(i,k,j))) -al(i,k,j))

 ENDDO
 ENDDO
 ENDDO

 ELSE

 DO j =j_start,j_end
 DO k =k_start,k_end
 DO i =i_start,i_end

 g_p(i,k,j) =g_mu(i,j)*znu(k)
 p(i,k,j) =mu(i,j)*znu(k)

 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) &
*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) &
*(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))) &
*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))) &
*(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)) &
/(c2a(i,k,j)*c2a(i,k,j)))
 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))) &
 -p(i,k,j)/c2a(i,k,j)

 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) &
*al(i,k,j) +(mu(i,j)*g_alt(i,k,j) +g_mu(i,j)*alt(i,k,j)))
 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))

 ENDDO
 ENDDO
 ENDDO
 END IF

 IF(step == 0) THEN

 DO j =j_start,j_end
 DO k =k_start,k_end
 DO i =i_start,i_end

 g_pm1(i,k,j) =g_p(i,k,j)
 pm1(i,k,j) =p(i,k,j)

 ENDDO
 ENDDO
 ENDDO

 ELSE

 DO j =j_start,j_end
 DO k =k_start,k_end
 DO i =i_start,i_end

 g_ptmp =g_p(i,k,j)
 ptmp =p(i,k,j)

 g_p(i,k,j) =g_p(i,k,j) +smdiv*(g_p(i,k,j) -g_pm1(i,k,j))
 p(i,k,j) =p(i,k,j) +smdiv*(p(i,k,j) -pm1(i,k,j))

 g_pm1(i,k,j) =g_ptmp
 pm1(i,k,j) =ptmp

 ENDDO
 ENDDO
 ENDDO
 END IF

 END SUBROUTINE g_calc_p_rho

 SUBROUTINE g_calc_coef_w(a,g_a,alpha,g_alpha,gamma,g_gamma,mut,g_mut, &
 cqw,g_cqw,rdn,rdnw,c2a,g_c2a,dts,g,epssm,top_lid,ids,ide,jds,jde,kds,kde,ims, &
 ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte)

 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ids,ide,jds,jde,kds,kde
 INTEGER,INTENT(IN) :: ims,ime,jms,jme,kms,kme
 INTEGER,INTENT(IN) :: its,ite,jts,jte,kts,kte
 LOGICAL,INTENT(IN) :: top_lid
 REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(IN) :: g_c2a,c2a,g_cqw,cqw
 REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(INOUT) :: g_alpha,alpha,g_gamma, &
 gamma,g_a,a
 REAL,DIMENSION(ims:ime,jms:jme),INTENT(IN) :: g_mut,mut
 REAL,DIMENSION(kms:kme),INTENT(IN) :: rdn,rdnw
 REAL,INTENT(IN) :: epssm,dts,g
 REAL,DIMENSION(ims:ime) :: g_cof,cof
 REAL :: g_b,b,g_c,c
 INTEGER :: i,j,k,i_start,i_end,j_start,j_end,k_start,k_end
 INTEGER :: ij,ijp,ijm,lid_flag

 i_start =its

 i_end =min(ite,ide-1)

 j_start =jts

 j_end =min(jte,jde-1)

 k_start =kts

 k_end =kte-1

 lid_flag =1

 IF(top_lid) lid_flag =0

 DO j =j_start,j_end
 DO i =i_start,i_end

 g_cof(i) =2.0*(-.5 *dts *g *(1.+epssm)*g_mut(i,j)/(mut(i,j)*mut(i,j))) &
*(.5 *dts *g *(1.+epssm)/mut(i,j))
 cof(i) =(.5 *dts *g *(1.+epssm)/mut(i,j))**2

 g_a(i,2,j) =0.0
 a(i,2,j) =0.

 g_a(i,kde,j) =(-2.*cof(i)*rdnw(kde-1)**2*g_c2a(i,kde-1,j) +(-2.)*g_cof(i) &
*rdnw(kde-1)**2*c2a(i,kde-1,j))*lid_flag
 a(i,kde,j) =-2.*cof(i)*rdnw(kde-1)**2*c2a(i,kde-1,j)*lid_flag

 g_gamma(i,1,j) =0.0
 gamma(i,1,j) =0.

 ENDDO

 DO k =3,kde-1
 DO i =i_start,i_end

 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) &
*g_cof(i) -g_cqw(i,k,j)*cof(i))*rdn(k)*rdnw(k-1)*c2a(i,k-1,j)
 a(i,k,j) =-cqw(i,k,j)*cof(i)*rdn(k)*rdnw(k-1)*c2a(i,k-1,j)

 ENDDO
 ENDDO

 DO k =2,kde-1
 DO i =i_start,i_end

 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, &
 j)) +(cqw(i,k,j)*g_cof(i) +g_cqw(i,k,j)*cof(i))*rdn(k)*(rdnw(k)*c2a(i,k,j) &
 +rdnw(k-1)*c2a(i,k-1,j))
 b =1. +cqw(i,k,j)*cof(i)*rdn(k)*(rdnw(k)*c2a(i,k,j) +rdnw(k-1)*c2a(i,k-1,j))

 g_c =-cqw(i,k,j)*cof(i)*rdn(k)*rdnw(k)*g_c2a(i,k,j) +(-cqw(i,k,j)*g_cof(i) &
 -g_cqw(i,k,j)*cof(i))*rdn(k)*rdnw(k)*c2a(i,k,j)
 c =-cqw(i,k,j)*cof(i)*rdn(k)*rdnw(k)*c2a(i,k,j)

 g_alpha(i,k,j) =-1.*(g_b -(a(i,k,j)*g_gamma(i,k-1,j) +g_a(i,k,j) &
*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)))
 alpha(i,k,j) =1./(b -a(i,k,j)*gamma(i,k-1,j))

 g_gamma(i,k,j) =c*g_alpha(i,k,j) +g_c*alpha(i,k,j)
 gamma(i,k,j) =c*alpha(i,k,j)

 ENDDO
 ENDDO

 DO i =i_start,i_end

 g_b =2.*cof(i)*rdnw(kde-1)**2*g_c2a(i,kde-1,j) +2.*g_cof(i)*rdnw(kde-1) &
**2*c2a(i,kde-1,j)
 b =1. +2.*cof(i)*rdnw(kde-1)**2*c2a(i,kde-1,j)

 g_c =0.0
 c =0.

 g_alpha(i,kde,j) =-1.*(g_b -(a(i,kde,j)*g_gamma(i,kde-1,j) +g_a(i,kde,j) &
*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)))
 alpha(i,kde,j) =1./(b -a(i,kde,j)*gamma(i,kde-1,j))

 g_gamma(i,kde,j) =c*g_alpha(i,kde,j) +g_c*alpha(i,kde,j)
 gamma(i,kde,j) =c*alpha(i,kde,j)

 ENDDO
 ENDDO

 END SUBROUTINE g_calc_coef_w

!        Generated by TAPENADE     (INRIA, Tropics team)
!  Tapenade 3.5 (r3785) - 22 Mar 2011 18:35
!
!  Differentiation of advance_uv in forward (tangent) mode:
!   variations   of useful results: u v
!   with respect to varying inputs: p al u v ru_tend mudf cqu cqv
!                php rv_tend ph alt muu muv mu
!   RW status of diff variables: p:in al:in u:in-out v:in-out ru_tend:in
!                mudf:in cqu:in cqv:in php:in rv_tend:in ph:in
!                alt:in muu:in muv:in mu:in
SUBROUTINE G_ADVANCE_UV(u, ud, ru_tend, ru_tendd, v, vd, rv_tend, &
&  rv_tendd, p, pd, pb, ph, phd, php, phpd, alt, altd, al, ald, mu, mud, &
&  muu, muud, cqu, cqud, muv, muvd, cqv, cqvd, mudf, mudfd, msfux, msfuy&
&  , msfvx, msfvx_inv, msfvy, rdx, rdy, dts, cf1, cf2, cf3, fnm, fnp, &
&  emdiv, rdnw, config_flags, spec_zone, non_hydrostatic, top_lid, ids, &
&  ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, &
&  jte, kts, kte)
  IMPLICIT NONE
! religion first
! stuff coming in
  TYPE(GRID_CONFIG_REC_TYPE), INTENT(IN) :: config_flags
  LOGICAL, INTENT(IN) :: non_hydrostatic, top_lid
  INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde
  INTEGER, INTENT(IN) :: ims, ime, jms, jme, kms, kme
  INTEGER, INTENT(IN) :: its, ite, jts, jte, kts, kte
  INTEGER, INTENT(IN) :: spec_zone
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: u, v
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: ud, vd
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: ru_tend, &
&  rv_tend, ph, php, p, pb, alt, al, cqu, cqv
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: ru_tendd, &
&  rv_tendd, phd, phpd, pd, altd, ald, cqud, cqvd
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: muu, muv, mu, mudf
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: muud, muvd, mud, &
&  mudfd
  REAL, DIMENSION(kms:kme), INTENT(IN) :: fnm, fnp, rdnw
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: msfux, msfuy, msfvx, &
&  msfvy, msfvx_inv
  REAL, INTENT(IN) :: rdx, rdy, dts, cf1, cf2, cf3, emdiv
!  Local 3d array from the stack (note tile size)
  REAL, DIMENSION(its:ite, kts:kte) :: dpn, dpxy
  REAL, DIMENSION(its:ite, kts:kte) :: dpnd, dpxyd
  REAL, DIMENSION(its:ite) :: mudf_xy
  REAL, DIMENSION(its:ite) :: mudf_xyd
  REAL :: dx, dy
  INTEGER :: i, j, k, i_start, i_end, j_start, j_end, k_start, k_end
  INTEGER :: i_endu, j_endv, k_endw
  INTEGER :: i_start_up, i_end_up, j_start_up, j_end_up
  INTEGER :: i_start_vp, i_end_vp, j_start_vp, j_end_vp
  INTEGER :: i_start_u_tend, i_end_u_tend, j_start_v_tend, j_end_v_tend
!<DESCRIPTION>
!
!  advance_uv advances the explicit perturbation horizontal momentum 
!  equations (u,v) by adding in the large-timestep tendency along with 
!  the small timestep pressure gradient tendency.
!
!</DESCRIPTION>
!  now, the real work.
!  set the loop bounds taking into account boundary conditions.
  IF (config_flags%nested .OR. config_flags%specified) THEN
    IF (its .LT. ids + spec_zone) THEN
      i_start = ids + spec_zone
    ELSE
      i_start = its
    END IF
    IF (ite .GT. ide - spec_zone - 1) THEN
      i_end = ide - spec_zone - 1
    ELSE
      i_end = ite
    END IF
    IF (jts .LT. jds + spec_zone) THEN
      j_start = jds + spec_zone
    ELSE
      j_start = jts
    END IF
    IF (jte .GT. jde - spec_zone - 1) THEN
      j_end = jde - spec_zone - 1
    ELSE
      j_end = jte
    END IF
    k_start = kts
    IF (kte .GT. kde - 1) THEN
      k_end = kde - 1
    ELSE
      k_end = kte
    END IF
    IF (ite .GT. ide - spec_zone) THEN
      i_endu = ide - spec_zone
    ELSE
      i_endu = ite
    END IF
    IF (jte .GT. jde - spec_zone) THEN
      j_endv = jde - spec_zone
    ELSE
      j_endv = jte
    END IF
    k_endw = k_end
    IF (config_flags%periodic_x) THEN
      i_start = its
      IF (ite .GT. ide - 1) THEN
        i_end = ide - 1
      ELSE
        i_end = ite
      END IF
      i_endu = ite
    END IF
  ELSE
    i_start = its
    IF (ite .GT. ide - 1) THEN
      i_end = ide - 1
    ELSE
      i_end = ite
    END IF
    j_start = jts
    IF (jte .GT. jde - 1) THEN
      j_end = jde - 1
    ELSE
      j_end = jte
    END IF
    k_start = kts
    k_end = kte - 1
    i_endu = ite
    j_endv = jte
    k_endw = k_end
  END IF
  i_start_up = i_start
  i_end_up = i_endu
  j_start_up = j_start
  j_end_up = j_end
  i_start_vp = i_start
  i_end_vp = i_end
  j_start_vp = j_start
  j_end_vp = j_endv
  IF ((config_flags%open_xs .OR. config_flags%symmetric_xs) .AND. its &
&      .EQ. ids) i_start_up = i_start_up + 1
  IF ((config_flags%open_xe .OR. config_flags%symmetric_xe) .AND. ite &
&      .EQ. ide) i_end_up = i_end_up - 1
  IF (((config_flags%open_ys .OR. config_flags%symmetric_ys) .OR. &
&      config_flags%polar) .AND. jts .EQ. jds) j_start_vp = j_start_vp + &
&      1
  IF (((config_flags%open_ye .OR. config_flags%symmetric_ye) .OR. &
&      config_flags%polar) .AND. jte .EQ. jde) j_end_vp = j_end_vp - 1
  i_start_u_tend = i_start
  i_end_u_tend = i_endu
  j_start_v_tend = j_start
  j_end_v_tend = j_endv
  IF (config_flags%symmetric_xs .AND. its .EQ. ids) i_start_u_tend = &
&      i_start_u_tend + 1
  IF (config_flags%symmetric_xe .AND. ite .EQ. ide) i_end_u_tend = &
&      i_end_u_tend - 1
  IF (config_flags%symmetric_ys .AND. jts .EQ. jds) j_start_v_tend = &
&      j_start_v_tend + 1
  IF (config_flags%symmetric_ye .AND. jte .EQ. jde) j_end_v_tend = &
&      j_end_v_tend - 1
  dx = 1./rdx
  dy = 1./rdy
  mudf_xyd = 0.0
  dpxyd = 0.0
  dpnd = 0.0
!  start real calculations.
!  first, u
u_outer_j_loop:DO j=j_start,j_end
    DO k=k_start,k_end
      DO i=i_start_u_tend,i_end_u_tend
        ud(i, k, j) = ud(i, k, j) + dts*ru_tendd(i, k, j)
        u(i, k, j) = u(i, k, j) + dts*ru_tend(i, k, j)
      END DO
    END DO
    DO i=i_start_up,i_end_up
      mudf_xyd(i) = -(emdiv*dx*(mudfd(i, j)-mudfd(i-1, j))/msfuy(i, j))
      mudf_xy(i) = -(emdiv*dx*(mudf(i, j)-mudf(i-1, j))/msfuy(i, j))
    END DO
    DO k=k_start,k_end
      DO i=i_start_up,i_end_up
!     Comments on map scale factors:   
!     x pressure gradient: ADT eqn 44, penultimate term on RHS
!        = -(mx/my)*(mu/rho)*partial dp/dx
!     [i.e., first rho->mu; 2nd still rho; alpha=1/rho]
!     Klemp et al. splits into 2 terms:
!        mu alpha partial dp/dx + partial dp/dnu * partial dphi/dx 
!     then into 4 terms:
!        mu alpha partial dp'/dx + nu mu alpha' partial dmubar/dx +
!      + mu partial dphi/dx + partial dphi'/dx * (partial dp'/dnu - mu')
!
!     first 3 terms:
!     ph, alt, p, al, pb not coupled
!     since we want tendency to fit ADT eqn 44 (coupled) we need to
!     multiply by (mx/my):
!
        dpxyd(i, k) = msfux(i, j)*.5*rdx*(muud(i, j)*(ph(i, k+1, j)-ph(i&
&          -1, k+1, j)+(ph(i, k, j)-ph(i-1, k, j))+(alt(i, k, j)+alt(i-1&
&          , k, j))*(p(i, k, j)-p(i-1, k, j))+(al(i, k, j)+al(i-1, k, j))&
&          *(pb(i, k, j)-pb(i-1, k, j)))+muu(i, j)*(phd(i, k+1, j)-phd(i-&
&          1, k+1, j)+phd(i, k, j)-phd(i-1, k, j)+(altd(i, k, j)+altd(i-1&
&          , k, j))*(p(i, k, j)-p(i-1, k, j))+(alt(i, k, j)+alt(i-1, k, j&
&          ))*(pd(i, k, j)-pd(i-1, k, j))+(pb(i, k, j)-pb(i-1, k, j))*(&
&          ald(i, k, j)+ald(i-1, k, j))))/msfuy(i, j)
        dpxy(i, k) = msfux(i, j)/msfuy(i, j)*.5*rdx*muu(i, j)*(ph(i, k+1&
&          , j)-ph(i-1, k+1, j)+(ph(i, k, j)-ph(i-1, k, j))+(alt(i, k, j)&
&          +alt(i-1, k, j))*(p(i, k, j)-p(i-1, k, j))+(al(i, k, j)+al(i-1&
&          , k, j))*(pb(i, k, j)-pb(i-1, k, j)))
      END DO
    END DO
    IF (non_hydrostatic) THEN
      DO i=i_start_up,i_end_up
        dpnd(i, 1) = .5*(cf1*(pd(i, 1, j)+pd(i-1, 1, j))+cf2*(pd(i, 2, j&
&          )+pd(i-1, 2, j))+cf3*(pd(i, 3, j)+pd(i-1, 3, j)))
        dpn(i, 1) = .5*(cf1*(p(i, 1, j)+p(i-1, 1, j))+cf2*(p(i, 2, j)+p(&
&          i-1, 2, j))+cf3*(p(i, 3, j)+p(i-1, 3, j)))
        dpnd(i, kde) = 0.0
        dpn(i, kde) = 0.
      END DO
      IF (top_lid) THEN
        DO i=i_start_up,i_end_up
          dpnd(i, kde) = .5*(cf1*(pd(i-1, kde-1, j)+pd(i, kde-1, j))+cf2&
&            *(pd(i-1, kde-2, j)+pd(i, kde-2, j))+cf3*(pd(i-1, kde-3, j)+&
&            pd(i, kde-3, j)))
          dpn(i, kde) = .5*(cf1*(p(i-1, kde-1, j)+p(i, kde-1, j))+cf2*(p&
&            (i-1, kde-2, j)+p(i, kde-2, j))+cf3*(p(i-1, kde-3, j)+p(i, &
&            kde-3, j)))
        END DO
      END IF
      DO k=k_start+1,k_end
        DO i=i_start_up,i_end_up
          dpnd(i, k) = .5*(fnm(k)*(pd(i, k, j)+pd(i-1, k, j))+fnp(k)*(pd&
&            (i, k-1, j)+pd(i-1, k-1, j)))
          dpn(i, k) = .5*(fnm(k)*(p(i, k, j)+p(i-1, k, j))+fnp(k)*(p(i, &
&            k-1, j)+p(i-1, k-1, j)))
        END DO
      END DO
!    Comments on map scale factors:
!    4th term:
!    php, dpn, mu not coupled
!    since we want tendency to fit ADT eqn 44 (coupled) we need to
!    multiply by (mx/my):
      DO k=k_start,k_end
        DO i=i_start_up,i_end_up
          dpxyd(i, k) = dpxyd(i, k) + msfux(i, j)*rdx*((phpd(i, k, j)-&
&            phpd(i-1, k, j))*(rdnw(k)*(dpn(i, k+1)-dpn(i, k))-.5*(mu(i-1&
&            , j)+mu(i, j)))+(php(i, k, j)-php(i-1, k, j))*(rdnw(k)*(dpnd&
&            (i, k+1)-dpnd(i, k))-.5*(mud(i-1, j)+mud(i, j))))/msfuy(i, j&
&            )
          dpxy(i, k) = dpxy(i, k) + msfux(i, j)/msfuy(i, j)*rdx*(php(i, &
&            k, j)-php(i-1, k, j))*(rdnw(k)*(dpn(i, k+1)-dpn(i, k))-.5*(&
&            mu(i-1, j)+mu(i, j)))
        END DO
      END DO
    END IF
    DO k=k_start,k_end
      DO i=i_start_up,i_end_up
        ud(i, k, j) = ud(i, k, j) - dts*(cqud(i, k, j)*dpxy(i, k)+cqu(i&
&          , k, j)*dpxyd(i, k)) + mudf_xyd(i)
        u(i, k, j) = u(i, k, j) - dts*cqu(i, k, j)*dpxy(i, k) + mudf_xy(&
&          i)
      END DO
    END DO
  END DO u_outer_j_loop
! now v
v_outer_j_loop:DO j=j_start_v_tend,j_end_v_tend
    DO k=k_start,k_end
      DO i=i_start,i_end
        vd(i, k, j) = vd(i, k, j) + dts*rv_tendd(i, k, j)
        v(i, k, j) = v(i, k, j) + dts*rv_tend(i, k, j)
      END DO
    END DO
    DO i=i_start,i_end
      mudf_xyd(i) = -(emdiv*dy*msfvx_inv(i, j)*(mudfd(i, j)-mudfd(i, j-1&
&        )))
      mudf_xy(i) = -(emdiv*dy*(mudf(i, j)-mudf(i, j-1))*msfvx_inv(i, j))
    END DO
    IF (j .GE. j_start_vp .AND. j .LE. j_end_vp) THEN
      DO k=k_start,k_end
        DO i=i_start,i_end
!      Comments on map scale factors:
!      y pressure gradient: ADT eqn 45, penultimate term on RHS
!         = -(my/mx)*(mu/rho)*partial dp/dy
!      [i.e., first rho->mu; 2nd still rho; alpha=1/rho]
!      Klemp et al. splits into 2 terms:
!         mu alpha partial dp/dy + partial dp/dnu * partial dphi/dy 
!      then into 4 terms:
!         mu alpha partial dp'/dy + nu mu alpha' partial dmubar/dy +
!       + mu partial dphi/dy + partial dphi'/dy * (partial dp'/dnu - mu')
!
!      first 3 terms:
!      ph, alt, p, al, pb not coupled
!      since we want tendency to fit ADT eqn 45 (coupled) we need to
!      multiply by (my/mx):
!      mudf_xy is NOT a map scale factor coupling
!      it is some sort of divergence damping
          dpxyd(i, k) = msfvy(i, j)*.5*rdy*(muvd(i, j)*(ph(i, k+1, j)-ph&
&            (i, k+1, j-1)+(ph(i, k, j)-ph(i, k, j-1))+(alt(i, k, j)+alt(&
&            i, k, j-1))*(p(i, k, j)-p(i, k, j-1))+(al(i, k, j)+al(i, k, &
&            j-1))*(pb(i, k, j)-pb(i, k, j-1)))+muv(i, j)*(phd(i, k+1, j)&
&            -phd(i, k+1, j-1)+phd(i, k, j)-phd(i, k, j-1)+(altd(i, k, j)&
&            +altd(i, k, j-1))*(p(i, k, j)-p(i, k, j-1))+(alt(i, k, j)+&
&            alt(i, k, j-1))*(pd(i, k, j)-pd(i, k, j-1))+(pb(i, k, j)-pb(&
&            i, k, j-1))*(ald(i, k, j)+ald(i, k, j-1))))/msfvx(i, j)
          dpxy(i, k) = msfvy(i, j)/msfvx(i, j)*.5*rdy*muv(i, j)*(ph(i, k&
&            +1, j)-ph(i, k+1, j-1)+(ph(i, k, j)-ph(i, k, j-1))+(alt(i, k&
&            , j)+alt(i, k, j-1))*(p(i, k, j)-p(i, k, j-1))+(al(i, k, j)+&
&            al(i, k, j-1))*(pb(i, k, j)-pb(i, k, j-1)))
        END DO
      END DO
      IF (non_hydrostatic) THEN
        DO i=i_start,i_end
          dpnd(i, 1) = .5*(cf1*(pd(i, 1, j)+pd(i, 1, j-1))+cf2*(pd(i, 2&
&            , j)+pd(i, 2, j-1))+cf3*(pd(i, 3, j)+pd(i, 3, j-1)))
          dpn(i, 1) = .5*(cf1*(p(i, 1, j)+p(i, 1, j-1))+cf2*(p(i, 2, j)+&
&            p(i, 2, j-1))+cf3*(p(i, 3, j)+p(i, 3, j-1)))
          dpnd(i, kde) = 0.0
          dpn(i, kde) = 0.
        END DO
        IF (top_lid) THEN
          DO i=i_start,i_end
            dpnd(i, kde) = .5*(cf1*(pd(i, kde-1, j-1)+pd(i, kde-1, j))+&
&              cf2*(pd(i, kde-2, j-1)+pd(i, kde-2, j))+cf3*(pd(i, kde-3, &
&              j-1)+pd(i, kde-3, j)))
            dpn(i, kde) = .5*(cf1*(p(i, kde-1, j-1)+p(i, kde-1, j))+cf2*&
&              (p(i, kde-2, j-1)+p(i, kde-2, j))+cf3*(p(i, kde-3, j-1)+p(&
&              i, kde-3, j)))
          END DO
        END IF
        DO k=k_start+1,k_end
          DO i=i_start,i_end
            dpnd(i, k) = .5*(fnm(k)*(pd(i, k, j)+pd(i, k, j-1))+fnp(k)*(&
&              pd(i, k-1, j)+pd(i, k-1, j-1)))
            dpn(i, k) = .5*(fnm(k)*(p(i, k, j)+p(i, k, j-1))+fnp(k)*(p(i&
&              , k-1, j)+p(i, k-1, j-1)))
          END DO
        END DO
!      Comments on map scale factors:
!      4th term:
!      php, dpn, mu not coupled
!      since we want tendency to fit ADT eqn 45 (coupled) we need to
!      multiply by (my/mx):
        DO k=k_start,k_end
          DO i=i_start,i_end
            dpxyd(i, k) = dpxyd(i, k) + msfvy(i, j)*rdy*((phpd(i, k, j)-&
&              phpd(i, k, j-1))*(rdnw(k)*(dpn(i, k+1)-dpn(i, k))-.5*(mu(i&
&              , j-1)+mu(i, j)))+(php(i, k, j)-php(i, k, j-1))*(rdnw(k)*(&
&              dpnd(i, k+1)-dpnd(i, k))-.5*(mud(i, j-1)+mud(i, j))))/&
&              msfvx(i, j)
            dpxy(i, k) = dpxy(i, k) + msfvy(i, j)/msfvx(i, j)*rdy*(php(i&
&              , k, j)-php(i, k, j-1))*(rdnw(k)*(dpn(i, k+1)-dpn(i, k))-&
&              .5*(mu(i, j-1)+mu(i, j)))
          END DO
        END DO
      END IF
      DO k=k_start,k_end
        DO i=i_start,i_end
          vd(i, k, j) = vd(i, k, j) - dts*(cqvd(i, k, j)*dpxy(i, k)+cqv(&
&            i, k, j)*dpxyd(i, k)) + mudf_xyd(i)
          v(i, k, j) = v(i, k, j) - dts*cqv(i, k, j)*dpxy(i, k) + &
&            mudf_xy(i)
        END DO
      END DO
    END IF
  END DO v_outer_j_loop
! The check for j_start_vp and j_end_vp makes sure that the edges in v
! are not updated.  Let's add a polar boundary condition check here for
! safety to ensure that v at the poles never gets to be non-zero in the
! small time steps.
  IF (config_flags%polar) THEN
    IF (jts .EQ. jds) THEN
      DO k=k_start,k_end
        DO i=i_start,i_end
          vd(i, k, jds) = 0.0
          v(i, k, jds) = 0.
        END DO
      END DO
    END IF
    IF (jte .EQ. jde) THEN
      DO k=k_start,k_end
        DO i=i_start,i_end
          vd(i, k, jde) = 0.0
          v(i, k, jde) = 0.
        END DO
      END DO
    END IF
  END IF
END SUBROUTINE G_ADVANCE_UV

 SUBROUTINE g_advance_mu_t(ww,g_ww,ww_1,g_ww_1,u,g_u,u_1,g_u_1,v, &
 g_v,v_1,g_v_1,mu,g_mu,mut,g_mut,muave,g_muave,muts,g_muts,muu,g_muu, &
 muv,g_muv,mudf,g_mudf,uam,g_uam,vam,g_vam,wwam,g_wwam,t,g_t, &
 t_1,g_t_1,t_ave,g_t_ave,ft,g_ft,mu_tend,g_mu_tend,rdx,rdy,dts,epssm,dnw, &
 fnm,fnp,rdnw,msfux,msfuy,msfvx,msfvx_inv,msfvy,msftx,msfty,step,config_flags,ids,ide, &
 jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte)

 IMPLICIT NONE
 TYPE(grid_config_rec_type),INTENT(IN) :: config_flags
 INTEGER,INTENT(IN) :: ids,ide,jds,jde,kds,kde
 INTEGER,INTENT(IN) :: ims,ime,jms,jme,kms,kme
 INTEGER,INTENT(IN) :: its,ite,jts,jte,kts,kte
 INTEGER,INTENT(IN) :: step
 REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(IN) :: g_u,u,g_v,v,g_u_1,u_1, &
 g_v_1,v_1,g_t_1,t_1,g_ft,ft
 REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(INOUT) :: g_ww,ww,g_ww_1,ww_1, &
 g_t,t,g_t_ave,t_ave,g_uam,uam,g_vam,vam,g_wwam,wwam
 REAL,DIMENSION(ims:ime,jms:jme),INTENT(IN) :: g_muu,muu,g_muv,muv,g_mut,mut, &
 msfux,msfuy,msfvx,msfvx_inv,msfvy,msftx,msfty,g_mu_tend,mu_tend
 REAL,DIMENSION(ims:ime,jms:jme),INTENT(OUT) :: g_muave,muave,g_muts,muts,g_mudf,mudf
 REAL,DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: g_mu,mu
 REAL,DIMENSION(kms:kme),INTENT(IN) :: fnm,fnp,dnw,rdnw
 REAL,INTENT(IN) :: rdx,rdy,dts,epssm
 REAL,DIMENSION(its:ite,kts:kte) :: g_wdtn,wdtn,g_dvdxi,dvdxi
 REAL,DIMENSION(its:ite) :: g_dmdt,dmdt
 INTEGER :: i,j,k,i_start,i_end,j_start,j_end,k_start,k_end
 INTEGER :: i_endu,j_endv
 REAL :: g_acc,acc

 i_start =its

 i_end =min(ite,ide-1)

 j_start =jts

 j_end =min(jte,jde-1)

 k_start =kts

 k_end =kte-1

 IF( .NOT. config_flags%periodic_x ) THEN

 IF( config_flags%specified .or. config_flags%nested ) THEN

 i_start =max(its,ids+1)

 i_end =min(ite,ide-2)
 ENDIF
 ENDIF

 IF( config_flags%specified .or. config_flags%nested ) THEN

 j_start =max(jts,jds+1)

 j_end =min(jte,jde-2)
 ENDIF

 i_endu =ite

 j_endv =jte

 DO j =j_start,j_end
 DO i =i_start,i_end

 g_dmdt(i) =0.0
 dmdt(i) =0.

 ENDDO

 DO k =k_start,k_end
 DO i =i_start,i_end

 g_dvdxi(i,k) =msftx(i,j) *msfty(i,j)*(rdy*((g_v(i,k,j+1) +(muv(i,j+1) &
*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) &
 +(muv(i,j)*g_v_1(i,k,j) +g_muv(i,j)*v_1(i,k,j))*msfvx_inv(i,j))) &
 +rdx*((g_u(i+1,k,j) +((muu(i+1,j)*g_u_1(i+1,k,j) +g_muu(i+1,j) &
*u_1(i+1,k,j))/msfuy(i+1,j))) -(g_u(i,k,j) +((muu(i,j)*g_u_1(i,k,j) &
 +g_muu(i,j)*u_1(i,k,j))/msfuy(i,j)))))
 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) &
*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) &
 +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))))

 g_dmdt(i) =g_dmdt(i) +dnw(k)*g_dvdxi(i,k)
 dmdt(i) =dmdt(i) +dnw(k)*dvdxi(i,k)

 ENDDO
 ENDDO

 DO i =i_start,i_end

 g_muave(i,j) =g_mu(i,j)
 muave(i,j) =mu(i,j)

 g_mu(i,j) =g_mu(i,j) +dts*(g_dmdt(i) +g_mu_tend(i,j))
 mu(i,j) =mu(i,j) +dts*(dmdt(i) +mu_tend(i,j))

 g_mudf(i,j) =(g_dmdt(i) +g_mu_tend(i,j))
 mudf(i,j) =(dmdt(i) +mu_tend(i,j))

 g_muts(i,j) =g_mut(i,j) +g_mu(i,j)
 muts(i,j) =mut(i,j) +mu(i,j)

 g_muave(i,j) =.5*((1.+epssm)*g_mu(i,j) +(1.-epssm)*g_muave(i,j))
 muave(i,j) =.5*((1.+epssm)*mu(i,j) +(1.-epssm)*muave(i,j))

 ENDDO

 DO k =2,k_end
 DO i =i_start,i_end

 g_ww(i,k,j) =g_ww(i,k-1,j) -(dnw(k-1)*(g_dmdt(i) +g_dvdxi(i,k-1) &
 +g_mu_tend(i,j))/msfty(i,j))
 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)

 ENDDO
 ENDDO

 DO k =1,k_end
 DO i =i_start,i_end

 g_ww(i,k,j) =g_ww(i,k,j) -g_ww_1(i,k,j)
 ww(i,k,j) =ww(i,k,j) -ww_1(i,k,j)

 ENDDO
 ENDDO
 ENDDO

 DO j =j_start,j_end
 DO k =1,k_end
 DO i =i_start,i_end

 g_t_ave(i,k,j) =g_t(i,k,j)
 t_ave(i,k,j) =t(i,k,j)

 g_t(i,k,j) =g_t(i,k,j) +msfty(i,j) *dts*g_ft(i,k,j)
 t(i,k,j) =t(i,k,j) +msfty(i,j) *dts*ft(i,k,j)

 ENDDO
 ENDDO
 ENDDO

 DO j =j_start,j_end
 DO i =i_start,i_end

 g_wdtn(i,1) =0.0
 wdtn(i,1) =0.

 g_wdtn(i,kde) =0.0
 wdtn(i,kde) =0.

 ENDDO

 DO k =2,k_end
 DO i =i_start,i_end

 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)) &
 +g_ww(i,k,j)*(fnm(k)*t_1(i,k,j) +fnp(k)*t_1(i,k-1,j))
 wdtn(i,k) =ww(i,k,j)*(fnm(k)*t_1(i,k,j) +fnp(k)*t_1(i,k-1,j))

 ENDDO
 ENDDO

 DO k =1,k_end
 DO i =i_start,i_end

 g_t(i,k,j) =g_t(i,k,j) -dts *msfty(i,j)*(msftx(i,j)*(.5 *rdy*(v(i,k,j+1) &
*(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)) &
 -(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) &
 +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)) &
 +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) &
 +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) &
*(g_wdtn(i,k+1) -g_wdtn(i,k)))
 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) &
 +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) &
*(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) &
*(wdtn(i,k+1) -wdtn(i,k)))

 ENDDO
 ENDDO
 ENDDO

 END SUBROUTINE g_advance_mu_t

!        Generated by TAPENADE     (INRIA, Tropics team)
!  Tapenade 3.5 (r3785) - 22 Mar 2011 18:35
!
!  Differentiation of advance_w in forward (tangent) mode:
!   variations   of useful results: w t_2ave ph
!   with respect to varying inputs: ww w_save muave gamma u alpha
!                v w rw_tend ph_1 t_2ave cqw muts t_1 t_2 ph alt
!                c2a mut ph_tend a
!   RW status of diff variables: ww:in w_save:in muave:in gamma:in
!                u:in alpha:in v:in w:in-out rw_tend:in ph_1:in
!                t_2ave:in-out cqw:in muts:in t_1:in t_2:in ph:in-out
!                alt:in c2a:in mut:in ph_tend:in a:in
! domain dims
! memory dims
SUBROUTINE G_ADVANCE_W(w, wd, rw_tend, rw_tendd, ww, wwd, w_save, &
&  w_saved, u, ud, v, vd, mu1, mut, mutd, muave, muaved, muts, mutsd, &
&  t_2ave, t_2aved, t_2, t_2d, t_1, t_1d, ph, phd, ph_1, ph_1d, phb, &
&  ph_tend, ph_tendd, ht, c2a, c2ad, cqw, cqwd, alt, altd, alb, a, ad, &
&  alpha, alphad, gamma, gammad, rdx, rdy, dts, t0, epssm, dnw, fnm, fnp&
&  , rdnw, rdn, cf1, cf2, cf3, msftx, msfty, config_flags, top_lid, ids, &
&  ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, &
&  jte, kts, kte)
  IMPLICIT NONE
! religion first
! stuff coming in
  TYPE(GRID_CONFIG_REC_TYPE), INTENT(IN) :: config_flags
  INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde
  INTEGER, INTENT(IN) :: ims, ime, jms, jme, kms, kme
  INTEGER, INTENT(IN) :: its, ite, jts, jte, kts, kte
  LOGICAL, INTENT(IN) :: top_lid
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: t_2ave, w&
&  , ph
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: t_2aved, &
&  wd, phd
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: rw_tend, ww&
&  , w_save, u, v, t_2, t_1, ph_1, phb, ph_tend, alpha, gamma, a, c2a, &
&  cqw, alb, alt
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: rw_tendd, &
&  wwd, w_saved, ud, vd, t_2d, t_1d, ph_1d, ph_tendd, alphad, gammad, ad&
&  , c2ad, cqwd, altd
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: mu1, mut, muave, muts&
&  , ht, msftx, msfty
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: mutd, muaved, mutsd
  REAL, DIMENSION(kms:kme), INTENT(IN) :: fnp, fnm, rdnw, rdn, dnw
  REAL, INTENT(IN) :: rdx, rdy, dts, cf1, cf2, cf3, t0, epssm
!  Stack based 3d data, tile size.
  REAL, DIMENSION(its:ite) :: mut_inv, msft_inv
  REAL, DIMENSION(its:ite) :: mut_invd, msft_invd
  REAL, DIMENSION(its:ite, kts:kte) :: rhs, wdwn
  REAL, DIMENSION(its:ite, kts:kte) :: rhsd, wdwnd
  INTEGER :: i, j, k, i_start, i_end, j_start, j_end, k_start, k_end
  REAL, DIMENSION(kts:kte) :: dampwt
  REAL, DIMENSION(kts:kte) :: dampwtd
  REAL :: htop, hbot, hdepth, hk
  REAL :: htopd, hbotd, hkd
  REAL :: pi, dampmag
  REAL :: arg1
  REAL :: arg1d
  REAL :: arg2
  REAL :: arg2d
!<DESCRIPTION>
!
!  advance_w advances the implicit w and geopotential equations.
!
!</DESCRIPTION>
!  set loop limits.
!  Currently set for periodic boundary conditions
  i_start = its
  IF (ite .GT. ide - 1) THEN
    i_end = ide - 1
  ELSE
    i_end = ite
  END IF
  j_start = jts
  IF (jte .GT. jde - 1) THEN
    j_end = jde - 1
  ELSE
    j_end = jte
  END IF
  k_start = kts
  k_end = kte - 1
  IF (.NOT.config_flags%periodic_x) THEN
    IF (config_flags%specified .OR. config_flags%nested) THEN
      IF (its .LT. ids + 1) THEN
        i_start = ids + 1
      ELSE
        i_start = its
      END IF
      IF (ite .GT. ide - 2) THEN
        i_end = ide - 2
      ELSE
        i_end = ite
      END IF
    END IF
  END IF
  IF (config_flags%specified .OR. config_flags%nested) THEN
    IF (jts .LT. jds + 1) THEN
      j_start = jds + 1
    ELSE
      j_start = jts
    END IF
    IF (jte .GT. jde - 2) THEN
      j_end = jde - 2
    ELSE
      j_end = jte
    END IF
  END IF
  pi = 4.*ATAN(1.)
  dampmag = dts*config_flags%dampcoef
  hdepth = config_flags%zdamp
! calculation of phi and w equations
!   Comments on map scale factors:
!   phi equation is:
!    partial d/dt(rho phi/my) = -mx partial d/dx(phi rho u/my) 
!                               -mx partial d/dy(phi rho v/mx)
!                               - partial d/dz(phi rho w/my) + rho g w/my
!   as with scalar equation, use uncoupled value (here phi) to find the
!   coupled tendency (rho phi/my)
!   here as usual rho -> ~'mu'
!
!   w eqn [divided by my] is:
!     partial d/dt(rho w/my) = -mx partial d/dx(w rho u/my)
!                              -mx partial d/dy(v rho v/mx)
!                              - partial d/dz(w rho w/my)
!                              +rho[(u*u+v*v)/a + 2 u omega cos(lat) -
!                                   (1/rho) partial dp/dz - g + Fz]/my
!   here as usual rho -> ~'mu'
!
!  'u,v,w' sent here must be coupled variables (= rho w/my etc.) by their usage
  DO i=i_start,i_end
    rhsd(i, 1) = 0.0
    rhs(i, 1) = 0.
  END DO
  mut_invd = 0.0
  wdwnd = 0.0
  rhsd = 0.0
  dampwtd = 0.0
j_loop_w:DO j=j_start,j_end
    DO i=i_start,i_end
      mut_invd(i) = -(mutd(i, j)/mut(i, j)**2)
      mut_inv(i) = 1./mut(i, j)
      msft_invd(i) = 0.0
      msft_inv(i) = 1./msfty(i, j)
    END DO
    DO k=1,k_end
      DO i=i_start,i_end
        t_2aved(i, k, j) = .5*((1.+epssm)*t_2d(i, k, j)+(1.-epssm)*&
&          t_2aved(i, k, j))
        t_2ave(i, k, j) = .5*((1.+epssm)*t_2(i, k, j)+(1.-epssm)*t_2ave(&
&          i, k, j))
        t_2aved(i, k, j) = ((t_2aved(i, k, j)+t0*muaved(i, j))*muts(i, j&
&          )*(t0+t_1(i, k, j))-(t_2ave(i, k, j)+muave(i, j)*t0)*(mutsd(i&
&          , j)*(t0+t_1(i, k, j))+muts(i, j)*t_1d(i, k, j)))/(muts(i, j)*&
&          (t0+t_1(i, k, j)))**2
        t_2ave(i, k, j) = (t_2ave(i, k, j)+muave(i, j)*t0)/(muts(i, j)*(&
&          t0+t_1(i, k, j)))
        wdwnd(i, k+1) = .5*rdnw(k)*((wwd(i, k+1, j)+wwd(i, k, j))*(ph_1(&
&          i, k+1, j)-ph_1(i, k, j)+phb(i, k+1, j)-phb(i, k, j))+(ww(i, k&
&          +1, j)+ww(i, k, j))*(ph_1d(i, k+1, j)-ph_1d(i, k, j)))
        wdwn(i, k+1) = .5*(ww(i, k+1, j)+ww(i, k, j))*rdnw(k)*(ph_1(i, k&
&          +1, j)-ph_1(i, k, j)+phb(i, k+1, j)-phb(i, k, j))
        rhsd(i, k+1) = dts*(ph_tendd(i, k+1, j)+.5*g*(1.-epssm)*wd(i, k+&
&          1, j))
        rhs(i, k+1) = dts*(ph_tend(i, k+1, j)+.5*g*(1.-epssm)*w(i, k+1, &
&          j))
      END DO
    END DO
!   building up RHS of phi equation
!   ph_tend contains terms 1 and 2; now adding 3 and 4 in stages:
!   here rhs = delta t [ph_tend + ~g*w/2 - ~ww * partial d phi/dz]
    DO k=2,k_end
      DO i=i_start,i_end
        rhsd(i, k) = rhsd(i, k) - dts*(fnm(k)*wdwnd(i, k+1)+fnp(k)*wdwnd&
&          (i, k))
        rhs(i, k) = rhs(i, k) - dts*(fnm(k)*wdwn(i, k+1)+fnp(k)*wdwn(i, &
&          k))
      END DO
    END DO
!  NOTE:  phi'' is not coupled with the map-scale factor  (1/m), 
!         but it's tendency is, so must multiply by msft here
!   Comments on map scale factors:
!   building up RHS of phi equation
!   ph_tend contains terms 1 and 2; now adding 3 and 4 in stages:
!   here rhs = ph_previous + (msft/mu)*[rhs_previous - ~ww * delta t * 
!                                                      partial d phi/dz]
!            = ph_previous + (msft*delta t/mu)*[ph_tend + ~g*w/2 - ~ww *
!                                                         partial d phi/dz]
    DO k=2,k_end+1
      DO i=i_start,i_end
        rhsd(i, k) = phd(i, k, j) + msfty(i, j)*(rhsd(i, k)*mut_inv(i)+&
&          rhs(i, k)*mut_invd(i))
        rhs(i, k) = ph(i, k, j) + msfty(i, j)*rhs(i, k)*mut_inv(i)
        IF (top_lid .AND. k .EQ. k_end + 1) THEN
          rhsd(i, k) = 0.0
          rhs(i, k) = 0.
        END IF
      END DO
    END DO
!  lower boundary condition on w
!   Comments on map scale factors:
!   Chain rule: if Z=Z(X,Y) [true at the surface] then
!      dZ/dt = dZ/dX * dX/dt + dZ/dY * dY/dt, U=dX/dt, V=dY/dt
!   using capitals to denote actual values
!   In mapped values, u=U, v=V, z=Z, 1/dX=mx/dx, 1/dY=my/dy
!      w = dz/dt = mx u dz/dx + my v dz/dy
!   [where dz/dx is just the surface height change between x
!    gridpoints, and dz/dy is the change between y gridpoints]
!   [cf1, cf2 and cf3 do vertical weighting of u or v values nearest
!    the surface]
    DO i=i_start,i_end
      wd(i, 1, j) = msfty(i, j)*.5*rdy*((ht(i, j+1)-ht(i, j))*(cf1*vd(i&
&        , 1, j+1)+cf2*vd(i, 2, j+1)+cf3*vd(i, 3, j+1))+(ht(i, j)-ht(i, j&
&        -1))*(cf1*vd(i, 1, j)+cf2*vd(i, 2, j)+cf3*vd(i, 3, j))) + msftx(&
&        i, j)*.5*rdx*((ht(i+1, j)-ht(i, j))*(cf1*ud(i+1, 1, j)+cf2*ud(i+&
&        1, 2, j)+cf3*ud(i+1, 3, j))+(ht(i, j)-ht(i-1, j))*(cf1*ud(i, 1, &
&        j)+cf2*ud(i, 2, j)+cf3*ud(i, 3, j)))
      w(i, 1, j) = msfty(i, j)*.5*rdy*((ht(i, j+1)-ht(i, j))*(cf1*v(i, 1&
&        , j+1)+cf2*v(i, 2, j+1)+cf3*v(i, 3, j+1))+(ht(i, j)-ht(i, j-1))*&
&        (cf1*v(i, 1, j)+cf2*v(i, 2, j)+cf3*v(i, 3, j))) + msftx(i, j)*.5&
&        *rdx*((ht(i+1, j)-ht(i, j))*(cf1*u(i+1, 1, j)+cf2*u(i+1, 2, j)+&
&        cf3*u(i+1, 3, j))+(ht(i, j)-ht(i-1, j))*(cf1*u(i, 1, j)+cf2*u(i&
&        , 2, j)+cf3*u(i, 3, j)))
    END DO
!
! Jammed 3 doubly nested loops over k/i into 1 for slight improvement
! in efficiency.  No change in results (bit-for-bit). JM 20040514
! (left a blank line where the other two k/i-loops were)
!
!   above surface, begin by adding delta t * previous (coupled) w tendency
    DO k=2,k_end
      DO i=i_start,i_end
        wd(i, k, j) = wd(i, k, j) + dts*rw_tendd(i, k, j) + msft_inv(i)*&
&          .5*dts*g*rdn(k)*((cqwd(i, k, j)*mut_inv(i)+cqw(i, k, j)*&
&          mut_invd(i))*(c2a(i, k, j)*rdnw(k)*((1.+epssm)*(rhs(i, k+1)-&
&          rhs(i, k))+(1.-epssm)*(ph(i, k+1, j)-ph(i, k, j)))-c2a(i, k-1&
&          , j)*rdnw(k-1)*((1.+epssm)*(rhs(i, k)-rhs(i, k-1))+(1.-epssm)*&
&          (ph(i, k, j)-ph(i, k-1, j))))+cqw(i, k, j)*mut_inv(i)*(rdnw(k)&
&          *(c2ad(i, k, j)*((1.+epssm)*(rhs(i, k+1)-rhs(i, k))+(1.-epssm)&
&          *(ph(i, k+1, j)-ph(i, k, j)))+c2a(i, k, j)*((1.+epssm)*(rhsd(i&
&          , k+1)-rhsd(i, k))+(1.-epssm)*(phd(i, k+1, j)-phd(i, k, j))))-&
&          rdnw(k-1)*(c2ad(i, k-1, j)*((1.+epssm)*(rhs(i, k)-rhs(i, k-1))&
&          +(1.-epssm)*(ph(i, k, j)-ph(i, k-1, j)))+c2a(i, k-1, j)*((1.+&
&          epssm)*(rhsd(i, k)-rhsd(i, k-1))+(1.-epssm)*(phd(i, k, j)-phd(&
&          i, k-1, j)))))) + dts*g*msft_inv(i)*(rdn(k)*((c2ad(i, k, j)*&
&          alt(i, k, j)+c2a(i, k, j)*altd(i, k, j))*t_2ave(i, k, j)+c2a(i&
&          , k, j)*alt(i, k, j)*t_2aved(i, k, j)-(c2ad(i, k-1, j)*alt(i, &
&          k-1, j)+c2a(i, k-1, j)*altd(i, k-1, j))*t_2ave(i, k-1, j)-c2a(&
&          i, k-1, j)*alt(i, k-1, j)*t_2aved(i, k-1, j))-muaved(i, j))
        w(i, k, j) = w(i, k, j) + dts*rw_tend(i, k, j) + msft_inv(i)*cqw&
&          (i, k, j)*(.5*dts*g*mut_inv(i)*rdn(k)*(c2a(i, k, j)*rdnw(k)*((&
&          1.+epssm)*(rhs(i, k+1)-rhs(i, k))+(1.-epssm)*(ph(i, k+1, j)-ph&
&          (i, k, j)))-c2a(i, k-1, j)*rdnw(k-1)*((1.+epssm)*(rhs(i, k)-&
&          rhs(i, k-1))+(1.-epssm)*(ph(i, k, j)-ph(i, k-1, j))))) + dts*g&
&          *msft_inv(i)*(rdn(k)*(c2a(i, k, j)*alt(i, k, j)*t_2ave(i, k, j&
&          )-c2a(i, k-1, j)*alt(i, k-1, j)*t_2ave(i, k-1, j))-muave(i, j)&
&          )
      END DO
    END DO
    k = k_end + 1
    DO i=i_start,i_end
      wd(i, k, j) = wd(i, k, j) + dts*rw_tendd(i, k, j) + msft_inv(i)*(-&
&        (.5*dts*g*rdnw(k-1)**2*2.*((mut_invd(i)*c2a(i, k-1, j)+mut_inv(i&
&        )*c2ad(i, k-1, j))*((1.+epssm)*(rhs(i, k)-rhs(i, k-1))+(1.-epssm&
&        )*(ph(i, k, j)-ph(i, k-1, j)))+mut_inv(i)*c2a(i, k-1, j)*((1.+&
&        epssm)*(rhsd(i, k)-rhsd(i, k-1))+(1.-epssm)*(phd(i, k, j)-phd(i&
&        , k-1, j)))))-dts*g*(2.*rdnw(k-1)*((c2ad(i, k-1, j)*alt(i, k-1, &
&        j)+c2a(i, k-1, j)*altd(i, k-1, j))*t_2ave(i, k-1, j)+c2a(i, k-1&
&        , j)*alt(i, k-1, j)*t_2aved(i, k-1, j))+muaved(i, j)))
      w(i, k, j) = w(i, k, j) + dts*rw_tend(i, k, j) + msft_inv(i)*(-(.5&
&        *dts*g*mut_inv(i)*rdnw(k-1)**2*2.*c2a(i, k-1, j)*((1.+epssm)*(&
&        rhs(i, k)-rhs(i, k-1))+(1.-epssm)*(ph(i, k, j)-ph(i, k-1, j))))-&
&        dts*g*(2.*rdnw(k-1)*c2a(i, k-1, j)*alt(i, k-1, j)*t_2ave(i, k-1&
&        , j)+muave(i, j)))
      IF (top_lid) THEN
        wd(i, k, j) = 0.0
        w(i, k, j) = 0.
      END IF
    END DO
    DO k=2,k_end+1
      DO i=i_start,i_end
        wd(i, k, j) = (wd(i, k, j)-ad(i, k, j)*w(i, k-1, j)-a(i, k, j)*&
&          wd(i, k-1, j))*alpha(i, k, j) + (w(i, k, j)-a(i, k, j)*w(i, k-&
&          1, j))*alphad(i, k, j)
        w(i, k, j) = (w(i, k, j)-a(i, k, j)*w(i, k-1, j))*alpha(i, k, j)
      END DO
    END DO
    DO k=k_end,2,-1
      DO i=i_start,i_end
        wd(i, k, j) = wd(i, k, j) - gammad(i, k, j)*w(i, k+1, j) - gamma&
&          (i, k, j)*wd(i, k+1, j)
        w(i, k, j) = w(i, k, j) - gamma(i, k, j)*w(i, k+1, j)
      END DO
    END DO
    IF (config_flags%damp_opt .EQ. 3) THEN
      DO k=k_end+1,2,-1
        DO i=i_start,i_end
          htopd = ph_1d(i, k_end+1, j)/g
          htop = (ph_1(i, k_end+1, j)+phb(i, k_end+1, j))/g
          hkd = ph_1d(i, k, j)/g
          hk = (ph_1(i, k, j)+phb(i, k, j))/g
          hbotd = htopd
          hbot = htop - hdepth
          dampwtd(k) = 0.0
          dampwt(k) = 0.
          IF (hk .GE. hbot) THEN
            arg1d = 0.5*pi*(hkd-hbotd)/hdepth
            arg1 = 0.5*pi*(hk-hbot)/hdepth
            arg2d = 0.5*pi*(hkd-hbotd)/hdepth
            arg2 = 0.5*pi*(hk-hbot)/hdepth
            dampwtd(k) = dampmag*(arg1d*COS(arg1)*SIN(arg2)+SIN(arg1)*&
&              arg2d*COS(arg2))
            dampwt(k) = dampmag*SIN(arg1)*SIN(arg2)
          END IF
          wd(i, k, j) = ((wd(i, k, j)-(dampwtd(k)*mut(i, j)+dampwt(k)*&
&            mutd(i, j))*w_save(i, k, j)-dampwt(k)*mut(i, j)*w_saved(i, k&
&            , j))*(1.+dampwt(k))-(w(i, k, j)-dampwt(k)*mut(i, j)*w_save(&
&            i, k, j))*dampwtd(k))/(1.+dampwt(k))**2
          w(i, k, j) = (w(i, k, j)-dampwt(k)*mut(i, j)*w_save(i, k, j))/&
&            (1.+dampwt(k))
        END DO
      END DO
    END IF
    DO k=k_end+1,2,-1
      DO i=i_start,i_end
        phd(i, k, j) = rhsd(i, k) + (msfty(i, j)*.5*dts*g*(1.+epssm)*wd(&
&          i, k, j)*muts(i, j)-msfty(i, j)*.5*dts*g*(1.+epssm)*w(i, k, j)&
&          *mutsd(i, j))/muts(i, j)**2
        ph(i, k, j) = rhs(i, k) + msfty(i, j)*.5*dts*g*(1.+epssm)*w(i, k&
&          , j)/muts(i, j)
      END DO
    END DO
  END DO j_loop_w
END SUBROUTINE G_ADVANCE_W

 SUBROUTINE g_sumflux(ru,g_ru,rv,g_rv,ww,g_ww,u_lin,g_u_lin,v_lin, &
 g_v_lin,ww_lin,g_ww_lin,muu,g_muu,muv,g_muv,ru_m,g_ru_m,rv_m,g_rv_m, &
 ww_m,g_ww_m,epssm,msfux,msfuy,msfvx,msfvx_inv,msfvy,iteration,number_of_small_timesteps, &
 ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte)

 IMPLICIT NONE
 INTEGER,INTENT(IN) :: number_of_small_timesteps
 INTEGER,INTENT(IN) :: iteration
 INTEGER,INTENT(IN) :: ids,ide,jds,jde,kds,kde
 INTEGER,INTENT(IN) :: ims,ime,jms,jme,kms,kme
 INTEGER,INTENT(IN) :: its,ite,jts,jte,kts,kte
 REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(IN) :: g_ru,ru,g_rv,rv,g_ww, &
 ww,g_u_lin,u_lin,g_v_lin,v_lin,g_ww_lin,ww_lin
 REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(INOUT) :: g_ru_m,ru_m,g_rv_m, &
 rv_m,g_ww_m,ww_m
 REAL,DIMENSION(ims:ime,jms:jme),INTENT(IN) :: g_muu,muu,g_muv,muv,msfux,msfuy, &
 msfvx,msfvy,msfvx_inv
 INTEGER :: mini,minj,mink
 REAL,INTENT(IN) :: epssm
 INTEGER :: i,j,k

 IF(iteration == 1 ) THEN

 DO j =jts,jte
 DO k =kts,kte
 DO i =its,ite

 g_ru_m(i,k,j) =0.0
 ru_m(i,k,j) =0.

 g_rv_m(i,k,j) =0.0
 rv_m(i,k,j) =0.

 g_ww_m(i,k,j) =0.0
 ww_m(i,k,j) =0.

 ENDDO
 ENDDO
 ENDDO
 ENDIF

 mini =min(ide-1,ite)

 minj =min(jde-1,jte)

 mink =min(kde-1,kte)

 DO j =jts,minj
 DO k =kts,mink
 DO i =its,mini

 g_ru_m(i,k,j) =g_ru_m(i,k,j) +g_ru(i,k,j)
 ru_m(i,k,j) =ru_m(i,k,j) +ru(i,k,j)

 g_rv_m(i,k,j) =g_rv_m(i,k,j) +g_rv(i,k,j)
 rv_m(i,k,j) =rv_m(i,k,j) +rv(i,k,j)

 g_ww_m(i,k,j) =g_ww_m(i,k,j) +g_ww(i,k,j)
 ww_m(i,k,j) =ww_m(i,k,j) +ww(i,k,j)

 ENDDO
 ENDDO
 ENDDO

 IF(ite .GT. mini) THEN

 DO j =jts,minj
 DO k =kts,mink
 DO i =mini+1,ite

 g_ru_m(i,k,j) =g_ru_m(i,k,j) +g_ru(i,k,j)
 ru_m(i,k,j) =ru_m(i,k,j) +ru(i,k,j)

 ENDDO
 ENDDO
 ENDDO
 END IF

 IF(jte .GT. minj) THEN

 DO j =minj+1,jte
 DO k =kts,mink
 DO i =its,mini

 g_rv_m(i,k,j) =g_rv_m(i,k,j) +g_rv(i,k,j)
 rv_m(i,k,j) =rv_m(i,k,j) +rv(i,k,j)

 ENDDO
 ENDDO
 ENDDO
 END IF

 IF( kte .GT. mink) THEN

 DO j =jts,minj
 DO k =mink+1,kte
 DO i =its,mini

 g_ww_m(i,k,j) =g_ww_m(i,k,j) +g_ww(i,k,j)
 ww_m(i,k,j) =ww_m(i,k,j) +ww(i,k,j)

 ENDDO
 ENDDO
 ENDDO
 END IF

 IF(iteration == number_of_small_timesteps) THEN

 DO j =jts,minj
 DO k =kts,mink
 DO i =its,mini

 g_ru_m(i,k,j) =g_ru_m(i,k,j)/number_of_small_timesteps +((muu(i,j) &
*g_u_lin(i,k,j) +g_muu(i,j)*u_lin(i,k,j))/msfuy(i,j))
 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)

 g_rv_m(i,k,j) =g_rv_m(i,k,j)/number_of_small_timesteps +(muv(i,j) &
*g_v_lin(i,k,j) +g_muv(i,j)*v_lin(i,k,j))*msfvx_inv(i,j)
 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)

 g_ww_m(i,k,j) =g_ww_m(i,k,j)/number_of_small_timesteps +g_ww_lin(i,k,j)
 ww_m(i,k,j) =ww_m(i,k,j)/number_of_small_timesteps +ww_lin(i,k,j)

 ENDDO
 ENDDO
 ENDDO

 IF(ite .GT. mini) THEN

 DO j =jts,minj
 DO k =kts,mink
 DO i =mini+1,ite

 g_ru_m(i,k,j) =g_ru_m(i,k,j)/number_of_small_timesteps +((muu(i,j) &
*g_u_lin(i,k,j) +g_muu(i,j)*u_lin(i,k,j))/msfuy(i,j))
 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)

 ENDDO
 ENDDO
 ENDDO
 END IF

 IF(jte .GT. minj) THEN

 DO j =minj+1,jte
 DO k =kts,mink
 DO i =its,mini

 g_rv_m(i,k,j) =g_rv_m(i,k,j)/number_of_small_timesteps +(muv(i,j) &
*g_v_lin(i,k,j) +g_muv(i,j)*v_lin(i,k,j))*msfvx_inv(i,j)
 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)

 ENDDO
 ENDDO
 ENDDO
 END IF

 IF( kte .GT. mink) THEN

 DO j =jts,minj
 DO k =mink+1,kte
 DO i =its,mini

 g_ww_m(i,k,j) =g_ww_m(i,k,j)/number_of_small_timesteps +g_ww_lin(i,k,j)
 ww_m(i,k,j) =ww_m(i,k,j)/number_of_small_timesteps +ww_lin(i,k,j)

 ENDDO
 ENDDO
 ENDDO
 END IF
 ENDIF

 END SUBROUTINE g_sumflux

 SUBROUTINE g_init_module_small_step

 END SUBROUTINE g_init_module_small_step

 END MODULE g_module_small_step_em

