! ======================================================================================
! 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
!
!
! 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.
!
!
! 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
!
!
! small_step_finish reconstructs the full uncoupled prognostic variables
! from the coupled perturbation variables used in the small timesteps.
!
!
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
!
!
! 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.
!
!
! 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
!
!
! advance_w advances the implicit w and geopotential equations.
!
!
! 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