diff --git a/model/dyn_core.F90 b/model/dyn_core.F90
index a91b2d5b6..2ec8a7469 100644
--- a/model/dyn_core.F90
+++ b/model/dyn_core.F90
@@ -557,7 +557,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap,
enddo
endif
call timing_on('UPDATE_DZ_C')
- call update_dz_c(is, ie, js, je, npz, ng, dt2, dp_ref, zs, gridstruct%area, ut, vt, gz, ws3, &
+ call update_dz_c(is, ie, js, je, npz, ng, dt2, flagstruct%dz_min, dp_ref, zs, gridstruct%area, ut, vt, gz, ws3, &
npx, npy, gridstruct%sw_corner, gridstruct%se_corner, &
gridstruct%ne_corner, gridstruct%nw_corner, bd, gridstruct%grid_type)
call timing_off('UPDATE_DZ_C')
@@ -695,16 +695,18 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap,
if ( flagstruct%do_vort_damp ) then
! damping on delp and vorticity:
nord_v(k)=0;
+ damp_vt(k) = 0.5*d2_divg
endif
d_con_k = 0.
- elseif ( (k<=MAX(2,flagstruct%n_sponge-1)).and. flagstruct%d2_bg_k2>0.01 ) then
+ elseif ( k<=MAX(2,flagstruct%n_sponge-1) .and. flagstruct%d2_bg_k2>0.01 ) then
nord_k=0; d2_divg = max(flagstruct%d2_bg, flagstruct%d2_bg_k2)
nord_w=0; damp_w = d2_divg
if ( flagstruct%do_vort_damp ) then
nord_v(k)=0;
+ damp_vt(k) = 0.5*d2_divg
endif
d_con_k = 0.
- elseif ( (k<=MAX(3,flagstruct%n_sponge)) .and. flagstruct%d2_bg_k2>0.05 ) then
+ elseif ( k<=MAX(3,flagstruct%n_sponge) .and. flagstruct%d2_bg_k2>0.05 ) then
nord_k=0; d2_divg = max(flagstruct%d2_bg, 0.2*flagstruct%d2_bg_k2)
nord_w=0; damp_w = d2_divg
d_con_k = 0.
@@ -843,7 +845,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap,
#ifndef SW_DYNAMICS
call timing_on('UPDATE_DZ')
call update_dz_d(nord_v, damp_vt, flagstruct%hord_tm, is, ie, js, je, npz, ng, npx, npy, gridstruct%area, &
- gridstruct%rarea, dp_ref, zs, zh, crx, cry, xfx, yfx, delz, ws, rdt, gridstruct, bd, flagstruct%lim_fac)
+ gridstruct%rarea, dp_ref, zs, zh, crx, cry, xfx, yfx, delz, ws, rdt, flagstruct%dz_min, gridstruct, bd, flagstruct%lim_fac)
call timing_off('UPDATE_DZ')
if (idiag%id_ws>0 .and. last_step) then
@@ -1165,12 +1167,8 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap,
do k=1,n_con
delt = abs(bdt*flagstruct%delt_max)
! Sponge layers:
- if ( flagstruct%n_sponge == 0) then
- if ( k == 1 ) delt = 0.1*delt
- if ( k == 2 ) delt = 0.5*delt
- else
- delt = delt*MIN(1.0,FLOAT(k)/FLOAT(flagstruct%n_sponge))
- endif
+ if ( k == 1 ) delt = 0.1*delt
+ if ( k == 2 ) delt = 0.5*delt
do j=js,je
do i=is,ie
#ifdef MOIST_CAPPA
@@ -1801,7 +1799,7 @@ subroutine one_grad_p(u, v, pk, gz, divg2, delp, dt, ng, gridstruct, bd, npx, np
!$OMP parallel do default(none) shared(is,ie,js,je,wk2,divg2)
do j=js,je+1
do i=is,ie
- wk2(i,j) = divg2(i,j)-divg2(i+1,j)
+ wk2(i,j) = (divg2(i,j)-divg2(i+1,j))
enddo
enddo
diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90
index 663eaa174..bc3062084 100644
--- a/model/fv_arrays.F90
+++ b/model/fv_arrays.F90
@@ -614,6 +614,8 @@ module fv_arrays_mod
!< and 1. The default value is 0.75. Only used if
!< 'hydrostatic' = .false.
+ real :: dz_min = 2.0 !< Controls minimum thickness in NH solver
+
integer :: n_split = 0 !< The number of small dynamics (acoustic) time steps between
!< vertical remapping. 0 by default, in which case the model
!< produces a good first guess by examining the resolution,
diff --git a/model/fv_control.F90 b/model/fv_control.F90
index 4c6dec079..45499d58f 100644
--- a/model/fv_control.F90
+++ b/model/fv_control.F90
@@ -222,6 +222,7 @@ module fv_control_mod
logical , pointer :: reset_eta
real , pointer :: p_fac
real , pointer :: a_imp
+ real , pointer :: dz_min
integer , pointer :: n_split
! Default
integer , pointer :: m_split
@@ -656,7 +657,7 @@ subroutine run_setup(Atm, dt_atmos, grids_on_this_pe, p_split)
namelist /fv_grid_nml/ grid_name, grid_file
namelist /fv_core_nml/npx, npy, ntiles, npz, npz_rst, layout, io_layout, ncnst, nwat, &
- use_logp, p_fac, a_imp, k_split, n_split, m_split, q_split, print_freq, write_3d_diags, do_schmidt, &
+ use_logp, p_fac, a_imp, dz_min, k_split, n_split, m_split, q_split, print_freq, write_3d_diags, do_schmidt, &
hord_mt, hord_vt, hord_tm, hord_dp, hord_tr, shift_fac, stretch_fac, target_lat, target_lon, &
kord_mt, kord_wz, kord_tm, kord_tr, fv_debug, fv_land, nudge, do_sat_adj, do_f3d, &
external_ic, read_increment, ncep_ic, nggps_ic, ecmwf_ic, use_new_ncep, use_ncep_phy, fv_diag_ic, &
@@ -916,6 +917,7 @@ subroutine run_setup(Atm, dt_atmos, grids_on_this_pe, p_split)
if(is_master()) then
write(*,*) 'Off center implicit scheme param=', a_imp
write(*,*) ' p_fac=', p_fac
+ write(*,*) ' dz_min=', dz_min
endif
endif
@@ -1220,6 +1222,7 @@ subroutine setup_pointers(Atm)
reset_eta => Atm%flagstruct%reset_eta
p_fac => Atm%flagstruct%p_fac
a_imp => Atm%flagstruct%a_imp
+ dz_min => Atm%flagstruct%dz_min
n_split => Atm%flagstruct%n_split
m_split => Atm%flagstruct%m_split
k_split => Atm%flagstruct%k_split
diff --git a/model/fv_dynamics.F90 b/model/fv_dynamics.F90
index 8e0a2e258..23b0d1dc1 100644
--- a/model/fv_dynamics.F90
+++ b/model/fv_dynamics.F90
@@ -922,7 +922,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
call range_check('VA_dyn', ua, is, ie, js, je, ng, npz, gridstruct%agrid, &
-280., 280., bad_range)
call range_check('TA_dyn', pt, is, ie, js, je, ng, npz, gridstruct%agrid, &
- 130., 335., bad_range)
+ 150., 333., bad_range)
if ( .not. hydrostatic ) then
call range_check('W_dyn', w, is, ie, js, je, ng, npz, gridstruct%agrid, &
-100., 100., bad_range)
diff --git a/model/nh_utils.F90 b/model/nh_utils.F90
index 13d7c6963..5e17c5566 100644
--- a/model/nh_utils.F90
+++ b/model/nh_utils.F90
@@ -63,18 +63,17 @@ module nh_utils_mod
public sim3p0_solver, rim_2d
public Riem_Solver_c
- real, parameter:: dz_min = 2.
real, parameter:: r3 = 1./3.
CONTAINS
- subroutine update_dz_c(is, ie, js, je, km, ng, dt, dp0, zs, area, ut, vt, gz, ws, &
+ subroutine update_dz_c(is, ie, js, je, km, ng, dt, dz_min, dp0, zs, area, ut, vt, gz, ws, &
npx, npy, sw_corner, se_corner, ne_corner, nw_corner, bd, grid_type)
! !INPUT PARAMETERS:
type(fv_grid_bounds_type), intent(IN) :: bd
integer, intent(in):: is, ie, js, je, ng, km, npx, npy, grid_type
logical, intent(IN):: sw_corner, se_corner, ne_corner, nw_corner
- real, intent(in):: dt
+ real, intent(in):: dt, dz_min
real, intent(in):: dp0(km)
real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: ut, vt
real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng):: area
@@ -193,7 +192,7 @@ subroutine update_dz_c(is, ie, js, je, km, ng, dt, dp0, zs, area, ut, vt, gz, ws
6000 continue
! Enforce monotonicity of height to prevent blowup
-!$OMP parallel do default(none) shared(is1,ie1,js1,je1,ws,zs,gz,rdt,km)
+!$OMP parallel do default(none) shared(is1,ie1,js1,je1,ws,zs,gz,rdt,dz_min,km)
do j=js1, je1
do i=is1, ie1
ws(i,j) = ( zs(i,j) - gz(i,j,km+1) ) * rdt
@@ -209,12 +208,12 @@ end subroutine update_dz_c
subroutine update_dz_d(ndif, damp, hord, is, ie, js, je, km, ng, npx, npy, area, rarea, &
- dp0, zs, zh, crx, cry, xfx, yfx, delz, ws, rdt, gridstruct, bd, lim_fac)
+ dp0, zs, zh, crx, cry, xfx, yfx, delz, ws, rdt, dz_min, gridstruct, bd, lim_fac)
type(fv_grid_bounds_type), intent(IN) :: bd
integer, intent(in):: is, ie, js, je, ng, km, npx, npy
integer, intent(in):: hord
- real, intent(in) :: rdt
+ real, intent(in) :: rdt, dz_min
real, intent(in) :: dp0(km)
real, intent(in) :: area(is-ng:ie+ng,js-ng:je+ng)
real, intent(in) :: rarea(is-ng:ie+ng,js-ng:je+ng)
@@ -308,7 +307,7 @@ subroutine update_dz_d(ndif, damp, hord, is, ie, js, je, km, ng, npx, npy, area,
enddo
-!$OMP parallel do default(none) shared(is,ie,js,je,km,ws,zs,zh,rdt)
+!$OMP parallel do default(none) shared(is,ie,js,je,km,ws,zs,zh,rdt,dz_min)
do j=js, je
do i=is,ie
ws(i,j) = ( zs(i,j) - zh(i,j,km+1) ) * rdt
diff --git a/model/sw_core.F90 b/model/sw_core.F90
index 52ef71523..9a13170a7 100644
--- a/model/sw_core.F90
+++ b/model/sw_core.F90
@@ -1,21 +1,21 @@
!***********************************************************************
-!* GNU Lesser General Public License
+!* GNU Lesser General Public License
!*
!* This file is part of the FV3 dynamical core.
!*
-!* The FV3 dynamical core is free software: you can redistribute it
+!* The FV3 dynamical core is free software: you can redistribute it
!* and/or modify it under the terms of the
!* GNU Lesser General Public License as published by the
-!* Free Software Foundation, either version 3 of the License, or
+!* Free Software Foundation, either version 3 of the License, or
!* (at your option) any later version.
!*
-!* The FV3 dynamical core is distributed in the hope that it will be
-!* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
-!* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+!* The FV3 dynamical core is distributed in the hope that it will be
+!* useful, but WITHOUT ANY WARRANTY; without even the implied warranty
+!* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
!* See the GNU General Public License for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
-!* License along with the FV3 dynamical core.
+!* License along with the FV3 dynamical core.
!* If not, see .
!***********************************************************************
@@ -23,7 +23,7 @@
!! as described by \cite lin1997explicit, \cite lin2004vertically, and \cite harris2013two.
!>@details The step is applied to the cubed sphere.
- module sw_core_mod
+module sw_core_mod
! Modules Included:
!
@@ -119,7 +119,7 @@ subroutine c_sw(delpc, delp, ptc, pt, u,v, w, uc,vc, ua,va, wc, &
type(fv_flags_type), intent(IN), target :: flagstruct
! Local:
- logical:: sw_corner, se_corner, ne_corner, nw_corner
+ logical:: sw_corner, se_corner, ne_corner, nw_corner
real, dimension(bd%is-1:bd%ie+1,bd%js-1:bd%je+1):: vort, ke
real, dimension(bd%is-1:bd%ie+2,bd%js-1:bd%je+1):: fx, fx1, fx2
real, dimension(bd%is-1:bd%ie+1,bd%js-1:bd%je+2):: fy, fy1, fy2
@@ -183,7 +183,7 @@ subroutine c_sw(delpc, delp, ptc, pt, u,v, w, uc,vc, ua,va, wc, &
do j=js-1,jep1
do i=is-1,iep1+1
if (ut(i,j) > 0.) then
- ut(i,j) = dt2*ut(i,j)*dy(i,j)*sin_sg(i-1,j,3)
+ ut(i,j) = dt2*ut(i,j)*dy(i,j)*sin_sg(i-1,j,3)
else
ut(i,j) = dt2*ut(i,j)*dy(i,j)*sin_sg(i,j,1)
end if
@@ -192,7 +192,7 @@ subroutine c_sw(delpc, delp, ptc, pt, u,v, w, uc,vc, ua,va, wc, &
do j=js-1,je+2
do i=is-1,iep1
if (vt(i,j) > 0.) then
- vt(i,j) = dt2*vt(i,j)*dx(i,j)*sin_sg(i,j-1,4)
+ vt(i,j) = dt2*vt(i,j)*dx(i,j)*sin_sg(i,j-1,4)
else
vt(i,j) = dt2*vt(i,j)*dx(i,j)*sin_sg(i,j, 2)
end if
@@ -236,7 +236,7 @@ subroutine c_sw(delpc, delp, ptc, pt, u,v, w, uc,vc, ua,va, wc, &
if (flagstruct%grid_type < 3) &
call fill_4corners(w, 1, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
do j=js-1,je+1
- do i=is-1,ie+2
+ do i=is-1,ie+2
if ( ut(i,j) > 0. ) then
fx1(i,j) = delp(i-1,j)
fx(i,j) = pt(i-1,j)
@@ -257,7 +257,7 @@ subroutine c_sw(delpc, delp, ptc, pt, u,v, w, uc,vc, ua,va, wc, &
if (flagstruct%grid_type < 3 .and. .not. nested) call fill2_4corners(delp, pt, 2, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
if ( hydrostatic ) then
do j=js-1,jep1+1
- do i=is-1,iep1
+ do i=is-1,iep1
if ( vt(i,j) > 0. ) then
fy1(i,j) = delp(i,j-1)
fy(i,j) = pt(i,j-1)
@@ -270,20 +270,20 @@ subroutine c_sw(delpc, delp, ptc, pt, u,v, w, uc,vc, ua,va, wc, &
enddo
enddo
do j=js-1,jep1
- do i=is-1,iep1
- delpc(i,j) = delp(i,j) + ((fx1(i,j)-fx1(i+1,j))+(fy1(i,j)-fy1(i,j+1)))*gridstruct%rarea(i,j)
+ do i=is-1,iep1
+ delpc(i,j) = delp(i,j) + (fx1(i,j)-fx1(i+1,j)+fy1(i,j)-fy1(i,j+1))*gridstruct%rarea(i,j)
#ifdef SW_DYNAMICS
ptc(i,j) = pt(i,j)
#else
ptc(i,j) = (pt(i,j)*delp(i,j) + &
- ((fx(i,j)-fx(i+1,j))+(fy(i,j)-fy(i,j+1)))*gridstruct%rarea(i,j))/delpc(i,j)
+ (fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1))*gridstruct%rarea(i,j))/delpc(i,j)
#endif
enddo
enddo
else
if (flagstruct%grid_type < 3) call fill_4corners(w, 2, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
do j=js-1,je+2
- do i=is-1,ie+1
+ do i=is-1,ie+1
if ( vt(i,j) > 0. ) then
fy1(i,j) = delp(i,j-1)
fy(i,j) = pt(i,j-1)
@@ -299,12 +299,12 @@ subroutine c_sw(delpc, delp, ptc, pt, u,v, w, uc,vc, ua,va, wc, &
enddo
enddo
do j=js-1,je+1
- do i=is-1,ie+1
- delpc(i,j) = delp(i,j) + ((fx1(i,j)-fx1(i+1,j))+(fy1(i,j)-fy1(i,j+1)))*gridstruct%rarea(i,j)
+ do i=is-1,ie+1
+ delpc(i,j) = delp(i,j) + (fx1(i,j)-fx1(i+1,j)+fy1(i,j)-fy1(i,j+1))*gridstruct%rarea(i,j)
ptc(i,j) = (pt(i,j)*delp(i,j) + &
- ((fx(i,j)-fx(i+1,j))+(fy(i,j)-fy(i,j+1)))*gridstruct%rarea(i,j))/delpc(i,j)
- wc(i,j) = (w(i,j)*delp(i,j) + ((fx2(i,j)-fx2(i+1,j)) + &
- (fy2(i,j)-fy2(i,j+1)))*gridstruct%rarea(i,j))/delpc(i,j)
+ (fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1))*gridstruct%rarea(i,j))/delpc(i,j)
+ wc(i,j) = (w(i,j)*delp(i,j) + (fx2(i,j)-fx2(i+1,j) + &
+ fy2(i,j)-fy2(i,j+1))*gridstruct%rarea(i,j))/delpc(i,j)
enddo
enddo
endif
@@ -313,7 +313,7 @@ subroutine c_sw(delpc, delp, ptc, pt, u,v, w, uc,vc, ua,va, wc, &
! Compute KE:
!------------
-!Since uc = u*, i.e. the covariant wind perpendicular to the face edge, if we want to compute kinetic energy we will need the true coordinate-parallel covariant wind, computed through u = uc*sina + v*cosa.
+!Since uc = u*, i.e. the covariant wind perpendicular to the face edge, if we want to compute kinetic energy we will need the true coordinate-parallel covariant wind, computed through u = uc*sina + v*cosa.
!Use the alpha for the cell KE is being computed in.
!!! TO DO:
!!! Need separate versions for nesting/single-tile
@@ -385,7 +385,7 @@ subroutine c_sw(delpc, delp, ptc, pt, u,v, w, uc,vc, ua,va, wc, &
dt4 = 0.5*dt2
do j=js-1,jep1
do i=is-1,iep1
- ke(i,j) = dt4*(ua(i,j)*ke(i,j) + va(i,j)*vort(i,j))
+ ke(i,j) = dt4*(ua(i,j)*ke(i,j) + va(i,j)*vort(i,j))
enddo
enddo
@@ -407,7 +407,7 @@ subroutine c_sw(delpc, delp, ptc, pt, u,v, w, uc,vc, ua,va, wc, &
do j=js,je+1
do i=is,ie+1
- vort(i,j) = (fx(i,j-1) - fx(i,j)) + (fy(i,j) - fy(i-1,j))
+ vort(i,j) = fx(i,j-1) - fx(i,j) - fy(i-1,j) + fy(i,j)
enddo
enddo
@@ -513,7 +513,7 @@ end subroutine c_sw
-! d_sw :: D-Grid Shallow Water Routine
+! d_sw :: D-Grid Shallow Water Routine
!>@brief The subroutine 'd_sw' peforms a full-timestep advance of the D-grid winds
!! and other prognostic varaiables.
@@ -559,7 +559,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
type(fv_grid_type), intent(IN), target :: gridstruct
type(fv_flags_type), intent(IN), target :: flagstruct
! Local:
- logical:: sw_corner, se_corner, ne_corner, nw_corner
+ logical:: sw_corner, se_corner, ne_corner, nw_corner
real :: ut(bd%isd:bd%ied+1,bd%jsd:bd%jed)
real :: vt(bd%isd:bd%ied, bd%jsd:bd%jed+1)
!---
@@ -575,7 +575,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
real :: fy(bd%is:bd%ie ,bd%js:bd%je+1) !< 1-D Y-direction Fluxes
real :: ra_x(bd%is:bd%ie,bd%jsd:bd%jed)
real :: ra_y(bd%isd:bd%ied,bd%js:bd%je)
- real :: gx(bd%is:bd%ie+1,bd%js:bd%je )
+ real :: gx(bd%is:bd%ie+1,bd%js:bd%je )
real :: gy(bd%is:bd%ie ,bd%js:bd%je+1) !< work Y-dir flux array
logical :: fill_c
@@ -612,37 +612,37 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
npy = flagstruct%npy
nested = gridstruct%nested
- area => gridstruct%area
- rarea => gridstruct%rarea
- sin_sg => gridstruct%sin_sg
- cosa_u => gridstruct%cosa_u
- cosa_v => gridstruct%cosa_v
- cosa_s => gridstruct%cosa_s
- sina_u => gridstruct%sina_u
- sina_v => gridstruct%sina_v
- rsin_u => gridstruct%rsin_u
- rsin_v => gridstruct%rsin_v
- rsina => gridstruct%rsina
- f0 => gridstruct%f0
- rsin2 => gridstruct%rsin2
- divg_u => gridstruct%divg_u
- divg_v => gridstruct%divg_v
- cosa => gridstruct%cosa
- dx => gridstruct%dx
- dy => gridstruct%dy
- dxc => gridstruct%dxc
- dyc => gridstruct%dyc
- rdxa => gridstruct%rdxa
- rdya => gridstruct%rdya
- rdx => gridstruct%rdx
- rdy => gridstruct%rdy
+ area => gridstruct%area
+ rarea => gridstruct%rarea
+ sin_sg => gridstruct%sin_sg
+ cosa_u => gridstruct%cosa_u
+ cosa_v => gridstruct%cosa_v
+ cosa_s => gridstruct%cosa_s
+ sina_u => gridstruct%sina_u
+ sina_v => gridstruct%sina_v
+ rsin_u => gridstruct%rsin_u
+ rsin_v => gridstruct%rsin_v
+ rsina => gridstruct%rsina
+ f0 => gridstruct%f0
+ rsin2 => gridstruct%rsin2
+ divg_u => gridstruct%divg_u
+ divg_v => gridstruct%divg_v
+ cosa => gridstruct%cosa
+ dx => gridstruct%dx
+ dy => gridstruct%dy
+ dxc => gridstruct%dxc
+ dyc => gridstruct%dyc
+ rdxa => gridstruct%rdxa
+ rdya => gridstruct%rdya
+ rdx => gridstruct%rdx
+ rdy => gridstruct%rdy
sw_corner = gridstruct%sw_corner
se_corner = gridstruct%se_corner
nw_corner = gridstruct%nw_corner
ne_corner = gridstruct%ne_corner
-#ifdef SW_DYNAMICS
+#ifdef SW_DYNAMICS
if ( test_case == 1 ) then
do j=jsd,jed
do i=is,ie+1
@@ -674,12 +674,12 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
!!! TO DO: separate versions for nesting and for cubed-sphere
if (nested) then
do j=jsd,jed
- do i=is-1,ie+2
+ do i=is,ie+1
ut(i,j) = ( uc(i,j) - 0.25 * cosa_u(i,j) * &
(vc(i-1,j)+vc(i,j)+vc(i-1,j+1)+vc(i,j+1)))*rsin_u(i,j)
enddo
enddo
- do j=js-1,je+2
+ do j=js,je+1
do i=isd,ied
vt(i,j) = ( vc(i,j) - 0.25 * cosa_v(i,j) * &
(uc(i,j-1)+uc(i+1,j-1)+uc(i,j)+uc(i+1,j)))*rsin_v(i,j)
@@ -778,10 +778,10 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
enddo
endif
-! The following code solves a 2x2 system to get the interior parallel-to-edge uc,vc values
-! near the corners (ex: for the sw corner ut(2,1) and vt(1,2) are solved for simultaneously).
-! It then computes the halo uc, vc values so as to be consistent with the computations on
-! the facing panel.
+! The following code solves a 2x2 system to get the interior parallel-to-edge uc,vc values
+! near the corners (ex: for the sw corner ut(2,1) and vt(1,2) are solved for simultaneously).
+! It then computes the halo uc, vc values so as to be consistent with the computations on
+! the facing panel.
!The system solved is:
! ut(2,1) = uc(2,1) - avg(vt)*cosa_u(2,1)
@@ -871,10 +871,10 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
ut(i,j) = uc(i,j)
enddo
enddo
-
+
do j=js,je+1
do i=isd,ied
- vt(i,j) = vc(i,j)
+ vt(i,j) = vc(i,j)
enddo
enddo
endif ! end grid_type choices
@@ -892,7 +892,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
enddo
! Explanation of the following code:
-! xfx_adv = dt*ut*dy
+! xfx_adv = dt*ut*dy
! crx_adv = dt*ut/dx
do j=jsd,jed
@@ -915,7 +915,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
yfx_adv(i,j) = dx(i,j)*yfx_adv(i,j)*sin_sg(i,j-1,4)
else
cry_adv(i,j) = yfx_adv(i,j) * rdya(i,j)
- yfx_adv(i,j) = dx(i,j)*yfx_adv(i,j)*sin_sg(i,j,2)
+ yfx_adv(i,j) = dx(i,j)*yfx_adv(i,j)*sin_sg(i,j,2)
endif
enddo
enddo
@@ -926,12 +926,12 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
do j=jsd,jed
do i=is,ie
- ra_x(i,j) = area(i,j) + (xfx_adv(i,j) - xfx_adv(i+1,j))
+ ra_x(i,j) = area(i,j) + xfx_adv(i,j) - xfx_adv(i+1,j)
enddo
enddo
do j=js,je
do i=isd,ied
- ra_y(i,j) = area(i,j) + (yfx_adv(i,j) - yfx_adv(i,j+1))
+ ra_y(i,j) = area(i,j) + yfx_adv(i,j) - yfx_adv(i,j+1)
enddo
enddo
@@ -944,12 +944,12 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
do i=is,ie+1
cx(i,j) = cx(i,j) + crx_adv(i,j)
enddo
- enddo
+ enddo
do j=js,je
do i=is,ie+1
xflux(i,j) = xflux(i,j) + fx(i,j)
enddo
- enddo
+ enddo
do j=js,je+1
do i=isd,ied
cy(i,j) = cy(i,j) + cry_adv(i,j)
@@ -957,7 +957,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
do i=is,ie
yflux(i,j) = yflux(i,j) + fy(i,j)
enddo
- enddo
+ enddo
#ifndef SW_DYNAMICS
do j=js,je
@@ -969,16 +969,18 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
if ( .not. hydrostatic ) then
if ( damp_w>1.E-5 ) then
- dd8 = kgb*abs(dt)
- damp4 = (damp_w*gridstruct%da_min_c)**(nord_w+1)
- call del6_vt_flux(nord_w, npx, npy, damp4, w, wk, fx2, fy2, gridstruct, bd)
+ dd8 = kgb*abs(dt)
+ damp4 = (damp_w*gridstruct%da_min_c)**(nord_w+1)
+ call del6_vt_flux(nord_w, npx, npy, damp4, w, wk, fx2, fy2, gridstruct, bd)
do j=js,je
do i=is,ie
- dw(i,j) = ((fx2(i,j)-fx2(i+1,j))+(fy2(i,j)-fy2(i,j+1)))*rarea(i,j)
+ dw(i,j) = (fx2(i,j)-fx2(i+1,j)+fy2(i,j)-fy2(i,j+1))*rarea(i,j)
! 0.5 * [ (w+dw)**2 - w**2 ] = w*dw + 0.5*dw*dw
! heat_source(i,j) = -d_con*dw(i,j)*(w(i,j)+0.5*dw(i,j))
heat_source(i,j) = dd8 - dw(i,j)*(w(i,j)+0.5*dw(i,j))
- diss_est(i,j) = heat_source(i,j)
+ if ( flagstruct%do_skeb ) then
+ diss_est(i,j) = heat_source(i,j)
+ endif
enddo
enddo
endif
@@ -986,19 +988,19 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
gridstruct, bd, ra_x, ra_y, flagstruct%lim_fac, mfx=fx, mfy=fy)
do j=js,je
do i=is,ie
- w(i,j) = delp(i,j)*w(i,j) + ((gx(i,j)-gx(i+1,j))+(gy(i,j)-gy(i,j+1)))*rarea(i,j)
+ w(i,j) = delp(i,j)*w(i,j) + (gx(i,j)-gx(i+1,j)+gy(i,j)-gy(i,j+1))*rarea(i,j)
enddo
enddo
endif
#ifdef USE_COND
- call fv_tp_2d(q_con, crx_adv,cry_adv, npx, npy, hord_dp, gx, gy, &
+ call fv_tp_2d(q_con, crx_adv,cry_adv, npx, npy, hord_dp, gx, gy, &
xfx_adv,yfx_adv, gridstruct, bd, ra_x, ra_y, flagstruct%lim_fac, mfx=fx, mfy=fy, mass=delp, nord=nord_t, damp_c=damp_t)
- do j=js,je
- do i=is,ie
- q_con(i,j) = delp(i,j)*q_con(i,j) + ((gx(i,j)-gx(i+1,j))+(gy(i,j)-gy(i,j+1)))*rarea(i,j)
- enddo
- enddo
+ do j=js,je
+ do i=is,ie
+ q_con(i,j) = delp(i,j)*q_con(i,j) + (gx(i,j)-gx(i+1,j)+gy(i,j)-gy(i,j+1))*rarea(i,j)
+ enddo
+ enddo
#endif
call fv_tp_2d(pt, crx_adv,cry_adv, npx, npy, hord_tm, gx, gy, &
@@ -1010,12 +1012,12 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
do j=js,je
do i=is,ie
wk(i,j) = delp(i,j)
- delp(i,j) = wk(i,j) + ((fx(i,j)-fx(i+1,j))+(fy(i,j)-fy(i,j+1)))*rarea(i,j)
+ delp(i,j) = wk(i,j) + (fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1))*rarea(i,j)
#ifdef SW_DYNAMICS
ptc(i,j) = pt(i,j)
#else
pt(i,j) = (pt(i,j)*wk(i,j) + &
- ((gx(i,j)-gx(i+1,j))+(gy(i,j)-gy(i,j+1)))*rarea(i,j))/delp(i,j)
+ (gx(i,j)-gx(i+1,j)+gy(i,j)-gy(i,j+1))*rarea(i,j))/delp(i,j)
#endif
enddo
enddo
@@ -1026,19 +1028,27 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
do j=js,je
do i=is,ie
q(i,j,k,iq) = (q(i,j,k,iq)*wk(i,j) + &
- ((gx(i,j)-gx(i+1,j))+(gy(i,j)-gy(i,j+1)))*rarea(i,j))/delp(i,j)
+ (gx(i,j)-gx(i+1,j)+gy(i,j)-gy(i,j+1))*rarea(i,j))/delp(i,j)
enddo
enddo
enddo
+! if ( zvir>0.01 ) then
+! do j=js,je
+! do i=is,ie
+! pt(i,j) = pt(i,j)*(1.+zvir*q(i,j,k,sphum))
+! enddo
+! enddo
+! endif
+
else
do j=js,je
do i=is,ie
#ifndef SW_DYNAMICS
pt(i,j) = pt(i,j)*delp(i,j) + &
- ((gx(i,j)-gx(i+1,j))+(gy(i,j)-gy(i,j+1)))*rarea(i,j)
+ (gx(i,j)-gx(i+1,j)+gy(i,j)-gy(i,j+1))*rarea(i,j)
#endif
delp(i,j) = delp(i,j) + &
- ((fx(i,j)-fx(i+1,j))+(fy(i,j)-fy(i,j+1)))*rarea(i,j)
+ (fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1))*rarea(i,j)
#ifndef SW_DYNAMICS
pt(i,j) = pt(i,j) / delp(i,j)
@@ -1079,7 +1089,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
if (nested) then
do j=js2,je1
do i=is2,ie1
- vb(i,j) = dt5*((vc(i-1,j)+vc(i,j))-(uc(i,j-1)+uc(i,j))*cosa(i,j))*rsina(i,j)
+ vb(i,j) = dt5*(vc(i-1,j)+vc(i,j)-(uc(i,j-1)+uc(i,j))*cosa(i,j))*rsina(i,j)
enddo
enddo
else
@@ -1091,7 +1101,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
do j=js2,je1
do i=is2,ie1
- vb(i,j) = dt5*((vc(i-1,j)+vc(i,j))-(uc(i,j-1)+uc(i,j))*cosa(i,j))*rsina(i,j)
+ vb(i,j) = dt5*(vc(i-1,j)+vc(i,j)-(uc(i,j-1)+uc(i,j))*cosa(i,j))*rsina(i,j)
enddo
if ( is==1 ) then
@@ -1110,7 +1120,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
enddo
endif
endif
-
+
else
do j=js,je+1
do i=is,ie+1
@@ -1133,9 +1143,9 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
if (nested) then
do j=js,je+1
-
+
do i=is2,ie1
- ub(i,j) = dt5*((uc(i,j-1)+uc(i,j))-(vc(i-1,j)+vc(i,j))*cosa(i,j))*rsina(i,j)
+ ub(i,j) = dt5*(uc(i,j-1)+uc(i,j)-(vc(i-1,j)+vc(i,j))*cosa(i,j))*rsina(i,j)
enddo
enddo
@@ -1156,7 +1166,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
enddo
else
do i=is2,ie1
- ub(i,j) = dt5*((uc(i,j-1)+uc(i,j))-(vc(i-1,j)+vc(i,j))*cosa(i,j))*rsina(i,j)
+ ub(i,j) = dt5*(uc(i,j-1)+uc(i,j)-(vc(i-1,j)+vc(i,j))*cosa(i,j))*rsina(i,j)
enddo
endif
enddo
@@ -1167,7 +1177,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
enddo
endif
endif
-
+
else
do j=js,je+1
do i=is,ie+1
@@ -1230,7 +1240,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
! wk is "volume-mean" relative vorticity
do j=jsd,jed
do i=isd,ied
- wk(i,j) = rarea(i,j)*( (vt(i,j)-vt(i,j+1)) + (ut(i+1,j)-ut(i,j)) )
+ wk(i,j) = rarea(i,j)*(vt(i,j)-vt(i,j+1)-ut(i,j)+ut(i+1,j))
enddo
enddo
@@ -1325,9 +1335,9 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
vort(1, j) = v(1, j)*dxc(1, j)*sin_sg(1,j,1)
end if
end if
- if ( (ie+1)==npx ) then
+ if ( (ie+1)==npx ) then
if (uc(npx,j) > 0) then
- vort(npx,j) = v(npx,j)*dxc(npx,j)* &
+ vort(npx,j) = v(npx,j)*dxc(npx,j)* &
sin_sg(npx-1,j,3)
else
vort(npx,j) = v(npx,j)*dxc(npx,j)* &
@@ -1339,7 +1349,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
do j=js,je+1
do i=is,ie+1
- delpc(i,j) = (vort(i,j-1) - vort(i,j)) + (ptc(i-1,j) - ptc(i,j))
+ delpc(i,j) = vort(i,j-1) - vort(i,j) + ptc(i-1,j) - ptc(i,j)
enddo
enddo
@@ -1393,7 +1403,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
if ( fill_c ) call fill_corners(vc, uc, npx, npy, VECTOR=.true., DGRID=.true.)
do j=js-nt,je+1+nt
do i=is-nt,ie+1+nt
- divg_d(i,j) = (uc(i,j-1) - uc(i,j)) + (vc(i-1,j) - vc(i,j))
+ divg_d(i,j) = uc(i,j-1) - uc(i,j) + vc(i-1,j) - vc(i,j)
enddo
enddo
@@ -1487,12 +1497,12 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
xfx_adv,yfx_adv, gridstruct, bd, ra_x, ra_y, flagstruct%lim_fac)
do j=js,je+1
do i=is,ie
- u(i,j) = vt(i,j) + (ke(i,j) - ke(i+1,j)) + fy(i,j)
+ u(i,j) = vt(i,j) + ke(i,j) - ke(i+1,j) + fy(i,j)
enddo
enddo
do j=js,je
do i=is,ie+1
- v(i,j) = ut(i,j) + (ke(i,j) - ke(i,j+1)) - fx(i,j)
+ v(i,j) = ut(i,j) + ke(i,j) - ke(i,j+1) - fx(i,j)
enddo
enddo
@@ -1538,7 +1548,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, &
diss_est(i,j) = diss_est(i,j)-rsin2(i,j)*( &
(ub(i,j)**2 + ub(i,j+1)**2 + vb(i,j)**2 + vb(i+1,j)**2) &
+ 2.*(gy(i,j)+gy(i,j+1)+gx(i,j)+gx(i+1,j)) &
- - cosa_s(i,j)*(u2*dv2 + v2*du2 + du2*dv2))
+ - cosa_s(i,j)*(u2*dv2 + v2*du2 + du2*dv2))
endif
enddo
enddo
@@ -1595,10 +1605,10 @@ subroutine del6_vt_flux(nord, npx, npy, damp, q, d2, fx2, fy2, gridstruct, bd)
#ifdef USE_SG
sin_sg => gridstruct%sin_sg
- rdxc => gridstruct%rdxc
- rdyc => gridstruct%rdyc
- dx => gridstruct%dx
- dy => gridstruct%dy
+ rdxc => gridstruct%rdxc
+ rdyc => gridstruct%rdyc
+ dx => gridstruct%dx
+ dy => gridstruct%dy
#endif
nested = gridstruct%nested
@@ -1606,7 +1616,7 @@ subroutine del6_vt_flux(nord, npx, npy, damp, q, d2, fx2, fy2, gridstruct, bd)
ie = bd%ie
js = bd%js
je = bd%je
-
+
i1 = is-1-nord; i2 = ie+1+nord
j1 = js-1-nord; j2 = je+1+nord
@@ -1645,7 +1655,7 @@ subroutine del6_vt_flux(nord, npx, npy, damp, q, d2, fx2, fy2, gridstruct, bd)
nt = nord-n
do j=js-nt-1,je+nt+1
do i=is-nt-1,ie+nt+1
- d2(i,j) = ((fx2(i,j)-fx2(i+1,j))+(fy2(i,j)-fy2(i,j+1)))*gridstruct%rarea(i,j)
+ d2(i,j) = (fx2(i,j)-fx2(i+1,j)+fy2(i,j)-fy2(i,j+1))*gridstruct%rarea(i,j)
enddo
enddo
@@ -1711,10 +1721,10 @@ subroutine divergence_corner(u, v, ua, va, divg_d, gridstruct, flagstruct, bd)
npy = flagstruct%npy
nested = gridstruct%nested
- sin_sg => gridstruct%sin_sg
- cos_sg => gridstruct%cos_sg
- dxc => gridstruct%dxc
- dyc => gridstruct%dyc
+ sin_sg => gridstruct%sin_sg
+ cos_sg => gridstruct%cos_sg
+ dxc => gridstruct%dxc
+ dyc => gridstruct%dyc
if (nested) then
is2 = is; ie1 = ie+1
@@ -1722,7 +1732,7 @@ subroutine divergence_corner(u, v, ua, va, divg_d, gridstruct, flagstruct, bd)
is2 = max(2,is); ie1 = min(npx-1,ie+1)
end if
- if (flagstruct%grid_type==4) then
+ if (flagstruct%grid_type > 3) then
do j=js-1,je+2
do i=is-2,ie+2
uf(i,j) = u(i,j)*dyc(i,j)
@@ -1735,7 +1745,7 @@ subroutine divergence_corner(u, v, ua, va, divg_d, gridstruct, flagstruct, bd)
enddo
do j=js-1,je+2
do i=is-1,ie+2
- divg_d(i,j) = gridstruct%rarea_c(i,j)*((vf(i,j-1)-vf(i,j))+(uf(i-1,j)-uf(i,j)))
+ divg_d(i,j) = gridstruct%rarea_c(i,j)*(vf(i,j-1)-vf(i,j)+uf(i-1,j)-uf(i,j))
enddo
enddo
else
@@ -1768,7 +1778,7 @@ subroutine divergence_corner(u, v, ua, va, divg_d, gridstruct, flagstruct, bd)
do j=js,je+1
do i=is,ie+1
- divg_d(i,j) = (vf(i,j-1) - vf(i,j)) + (uf(i-1,j) - uf(i,j))
+ divg_d(i,j) = vf(i,j-1) - vf(i,j) + uf(i-1,j) - uf(i,j)
enddo
enddo
@@ -1824,18 +1834,18 @@ subroutine divergence_corner_nest(u, v, ua, va, divg_d, gridstruct, flagstruct,
nested = gridstruct%nested
rarea_c => gridstruct%rarea_c
- sin_sg => gridstruct%sin_sg
- cos_sg => gridstruct%cos_sg
- cosa_u => gridstruct%cosa_u
- cosa_v => gridstruct%cosa_v
- sina_u => gridstruct%sina_u
- sina_v => gridstruct%sina_v
- dxc => gridstruct%dxc
- dyc => gridstruct%dyc
+ sin_sg => gridstruct%sin_sg
+ cos_sg => gridstruct%cos_sg
+ cosa_u => gridstruct%cosa_u
+ cosa_v => gridstruct%cosa_v
+ sina_u => gridstruct%sina_u
+ sina_v => gridstruct%sina_v
+ dxc => gridstruct%dxc
+ dyc => gridstruct%dyc
divg_d = 1.e25
- if (flagstruct%grid_type==4) then
+ if (flagstruct%grid_type > 3) then
do j=jsd,jed
do i=isd,ied
uf(i,j) = u(i,j)*dyc(i,j)
@@ -1848,7 +1858,7 @@ subroutine divergence_corner_nest(u, v, ua, va, divg_d, gridstruct, flagstruct,
enddo
do j=jsd+1,jed
do i=isd+1,ied
- divg_d(i,j) = rarea_c(i,j)*((vf(i,j-1)-vf(i,j))+(uf(i-1,j)-uf(i,j)))
+ divg_d(i,j) = rarea_c(i,j)*(vf(i,j-1)-vf(i,j)+uf(i-1,j)-uf(i,j))
enddo
enddo
else
@@ -1869,30 +1879,10 @@ subroutine divergence_corner_nest(u, v, ua, va, divg_d, gridstruct, flagstruct,
do j=jsd+1,jed
do i=isd+1,ied
- divg_d(i,j) = ((vf(i,j-1) - vf(i,j)) + (uf(i-1,j) - uf(i,j)))*rarea_c(i,j)
+ divg_d(i,j) = (vf(i,j-1) - vf(i,j) + uf(i-1,j) - uf(i,j))*rarea_c(i,j)
enddo
enddo
-!!$ !Edges
-!!$
-!!$ !West, East
-!!$ do j=jsd+1,jed
-!!$ divg_d(isd ,j) = (vf(isd,j-1) - vf(isd,j) + uf(isd,j) - uf(isd+1,j))*rarea_c(isd,j)
-!!$ divg_d(ied+1,j) = (vf(ied+1,j-1) - vf(ied+1,j) + uf(ied-1,j) - uf(ied,j))*rarea_c(ied,j)
-!!$ end do
-!!$
-!!$ !North, South
-!!$ do i=isd+1,ied
-!!$ divg_d(i,jsd ) = (vf(i,jsd) - vf(i,jsd+1) + uf(i-1,jsd) - uf(i,jsd))*rarea_c(i,jsd)
-!!$ divg_d(i,jed+1) = (vf(i,jed-1) - vf(i,jed) + uf(i-1,jed+1) - uf(i,jed+1))*rarea_c(i,jed)
-!!$ end do
-!!$
-!!$ !Corners (just use next corner value)
-!!$ divg_d(isd,jsd) = divg_d(isd+1,jsd+1)
-!!$ divg_d(isd,jed+1) = divg_d(isd+1,jed)
-!!$ divg_d(ied+1,jsd) = divg_d(ied,jsd+1)
-!!$ divg_d(ied+1,jed+1) = divg_d(ied,jed)
-
endif
@@ -1922,7 +1912,7 @@ subroutine smag_corner(dt, u, v, ua, va, smag_c, bd, npx, npy, gridstruct, ng)
integer :: is, ie, js, je
integer :: isd, ied, jsd, jed
-
+
is = bd%is
ie = bd%ie
js = bd%js
@@ -1957,7 +1947,7 @@ subroutine smag_corner(dt, u, v, ua, va, smag_c, bd, npx, npy, gridstruct, ng)
enddo
do j=js,je+1
do i=is,ie+1
- smag_c(i,j) = rarea_c(i,j)*((vt(i,j-1)-vt(i,j)) + (ut(i,j)-ut(i-1,j)))
+ smag_c(i,j) = rarea_c(i,j)*(vt(i,j-1)-vt(i,j)-ut(i-1,j)+ut(i,j))
enddo
enddo
! Fix the corners?? if grid_type /= 4
@@ -1976,7 +1966,7 @@ subroutine smag_corner(dt, u, v, ua, va, smag_c, bd, npx, npy, gridstruct, ng)
do j=jsd,jed
do i=isd,ied
- wk(i,j) = rarea(i,j)*((vt(i,j)-vt(i,j+1))+(ut(i,j)-ut(i+1,j)))
+ wk(i,j) = rarea(i,j)*(vt(i,j)-vt(i,j+1)+ut(i,j)-ut(i+1,j))
enddo
enddo
call a2b_ord4(wk, sh, gridstruct, npx, npy, is, ie, js, je, ng, .false.)
@@ -2023,7 +2013,7 @@ subroutine xtp_u(is,ie,js,je,isd,ied,jsd,jed,c, u, v, flux, iord, dx, rdx, npx,
if ( iord < 8 ) then
-! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6
+! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6
do j=js,je+1
@@ -2182,19 +2172,6 @@ subroutine xtp_u(is,ie,js,je,isd,ied,jsd,jed,c, u, v, flux, iord, dx, rdx, npx,
do i=is-1, ie+1
smt5(i) = 3.*abs(b0(i)) < abs(bl(i)-br(i))
enddo
-!WMP
-! fix edge issues
- if ( (.not.nested) .and. grid_type < 3) then
- if( is==1 ) then
- smt5(0) = bl(0)*br(0) < 0.
- smt5(1) = bl(1)*br(1) < 0.
- endif
- if( (ie+1)==npx ) then
- smt5(npx-1) = bl(npx-1)*br(npx-1) < 0.
- smt5(npx ) = bl(npx )*br(npx ) < 0.
- endif
- endif
-!WMP
endif
!DEC$ VECTOR ALWAYS
@@ -2342,7 +2319,7 @@ subroutine xtp_u(is,ie,js,je,isd,ied,jsd,jed,c, u, v, flux, iord, dx, rdx, npx,
br(i) = min(max(0., pmp, lac), max(al(i+1)-u(i,j), min(0.,pmp, lac)))
enddo
endif
-
+
do i=is,ie+1
if( c(i,j)>0. ) then
cfl = c(i,j)*rdx(i-1,j)
@@ -2391,7 +2368,7 @@ subroutine ytp_v(is,ie,js,je,isd,ied,jsd,jed, c, u, v, flux, jord, dy, rdy, npx,
end if
if ( jord<8 ) then
-! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6
+! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6
do j=js3,je3+1
do i=is,ie+1
@@ -2604,23 +2581,6 @@ subroutine ytp_v(is,ie,js,je,isd,ied,jsd,jed, c, u, v, flux, jord, dy, rdy, npx,
smt6(i,j) = 3.*abs(b0(i,j)) < abs(bl(i,j)-br(i,j))
enddo
enddo
-!WMP
-! fix edge issues
- if ( (.not.nested) .and. grid_type < 3) then
- if( js==1 ) then
- do i=is,ie+1
- smt6(i,0) = bl(i,0)*br(i,0) < 0.
- smt6(i,1) = bl(i,1)*br(i,1) < 0.
- enddo
- endif
- if( (je+1)==npy ) then
- do i=is,ie+1
- smt6(i,npy-1) = bl(i,npy-1)*br(i,npy-1) < 0.
- smt6(i,npy ) = bl(i,npy )*br(i,npy ) < 0.
- enddo
- endif
- endif
-!WMP
do j=js,je+1
!DEC$ VECTOR ALWAYS
do i=is,ie+1
@@ -2663,7 +2623,7 @@ subroutine ytp_v(is,ie,js,je,isd,ied,jsd,jed, c, u, v, flux, jord, dy, rdy, npx,
al(i,j) = 0.5*(v(i,j-1)+v(i,j)) + r3*(dm(i,j-1)-dm(i,j))
enddo
enddo
-
+
if ( jord==8 ) then
do j=js3,je3
do i=is,ie+1
@@ -2675,7 +2635,7 @@ subroutine ytp_v(is,ie,js,je,isd,ied,jsd,jed, c, u, v, flux, jord, dy, rdy, npx,
elseif ( jord==9 ) then
do j=js3,je3
do i=is,ie+1
- pmp_1 = -2.*dq(i,j)
+ pmp_1 = -2.*dq(i,j)
lac_1 = pmp_1 + 1.5*dq(i,j+1)
bl(i,j) = min(max(0., pmp_1, lac_1), max(al(i,j)-v(i,j), min(0., pmp_1, lac_1)))
pmp_2 = 2.*dq(i,j-1)
@@ -2695,7 +2655,7 @@ subroutine ytp_v(is,ie,js,je,isd,ied,jsd,jed, c, u, v, flux, jord, dy, rdy, npx,
br(i,j) = 0.
endif
elseif( abs(3.*(bl(i,j)+br(i,j))) > abs(bl(i,j)-br(i,j)) ) then
- pmp_1 = -2.*dq(i,j)
+ pmp_1 = -2.*dq(i,j)
lac_1 = pmp_1 + 1.5*dq(i,j+1)
bl(i,j) = min(max(0., pmp_1, lac_1), max(bl(i,j), min(0., pmp_1, lac_1)))
pmp_2 = 2.*dq(i,j-1)
@@ -2713,7 +2673,7 @@ subroutine ytp_v(is,ie,js,je,isd,ied,jsd,jed, c, u, v, flux, jord, dy, rdy, npx,
enddo
enddo
endif
-
+
!--------------
! fix the edges
!--------------
@@ -2803,18 +2763,18 @@ subroutine ytp_v(is,ie,js,je,isd,ied,jsd,jed, c, u, v, flux, jord, dy, rdy, npx,
al(i,j) = 0.5*(v(i,j-1)+v(i,j)) + r3*(dm(i,j-1)-dm(i,j))
enddo
enddo
-
+
do j=js-1,je+1
do i=is,ie+1
pmp = 2.*dq(i,j-1)
lac = pmp - 1.5*dq(i,j-2)
br(i,j) = min(max(0.,pmp,lac), max(al(i,j+1)-v(i,j), min(0.,pmp,lac)))
- pmp = -2.*dq(i,j)
+ pmp = -2.*dq(i,j)
lac = pmp + 1.5*dq(i,j+1)
bl(i,j) = min(max(0.,pmp,lac), max(al(i,j)-v(i,j), min(0.,pmp,lac)))
enddo
enddo
-
+
endif
do j=js,je+1
@@ -2851,7 +2811,7 @@ subroutine d2a2c_vect(u, v, ua, va, uc, vc, ut, vt, dord4, gridstruct, &
integer, intent(IN) :: npx, npy, grid_type
logical, intent(IN) :: nested
type(fv_grid_type), intent(IN), target :: gridstruct
-! Local
+! Local
real, dimension(bd%isd:bd%ied,bd%jsd:bd%jed):: utmp, vtmp
integer npt, i, j, ifirst, ilast, id
integer :: is, ie, js, je
@@ -2871,15 +2831,15 @@ subroutine d2a2c_vect(u, v, ua, va, uc, vc, ut, vt, dord4, gridstruct, &
jsd = bd%jsd
jed = bd%jed
- sin_sg => gridstruct%sin_sg
- cosa_u => gridstruct%cosa_u
- cosa_v => gridstruct%cosa_v
- cosa_s => gridstruct%cosa_s
- rsin_u => gridstruct%rsin_u
- rsin_v => gridstruct%rsin_v
- rsin2 => gridstruct%rsin2
- dxa => gridstruct%dxa
- dya => gridstruct%dya
+ sin_sg => gridstruct%sin_sg
+ cosa_u => gridstruct%cosa_u
+ cosa_v => gridstruct%cosa_v
+ cosa_s => gridstruct%cosa_s
+ rsin_u => gridstruct%rsin_u
+ rsin_v => gridstruct%rsin_v
+ rsin2 => gridstruct%rsin2
+ dxa => gridstruct%dxa
+ dya => gridstruct%dya
if ( dord4 ) then
id = 1
@@ -2895,7 +2855,7 @@ subroutine d2a2c_vect(u, v, ua, va, uc, vc, ut, vt, dord4, gridstruct, &
! Initialize the non-existing corner regions
utmp(:,:) = big_number
- vtmp(:,:) = big_number
+ vtmp(:,:) = big_number
if ( nested) then
@@ -2916,7 +2876,7 @@ subroutine d2a2c_vect(u, v, ua, va, uc, vc, ut, vt, dord4, gridstruct, &
vtmp(i,j) = a2*(v(i-1,j)+v(i+2,j)) + a1*(v(i,j)+v(i+1,j))
enddo
i = isd
- vtmp(i,j) = 0.5*(v(i,j)+v(i+1,j))
+ vtmp(i,j) = 0.5*(v(i,j)+v(i+1,j))
i = ied
vtmp(i,j) = 0.5*(v(i,j)+v(i+1,j))
enddo
@@ -3041,24 +3001,24 @@ subroutine d2a2c_vect(u, v, ua, va, uc, vc, ut, vt, dord4, gridstruct, &
! Xdir:
if( gridstruct%sw_corner ) then
ua(-1,0) = -va(0,2)
- ua( 0,0) = -va(0,1)
+ ua( 0,0) = -va(0,1)
endif
if( gridstruct%se_corner ) then
ua(npx, 0) = va(npx,1)
- ua(npx+1,0) = va(npx,2)
+ ua(npx+1,0) = va(npx,2)
endif
if( gridstruct%ne_corner ) then
ua(npx, npy) = -va(npx,npy-1)
- ua(npx+1,npy) = -va(npx,npy-2)
+ ua(npx+1,npy) = -va(npx,npy-2)
endif
if( gridstruct%nw_corner ) then
ua(-1,npy) = va(0,npy-2)
- ua( 0,npy) = va(0,npy-1)
+ ua( 0,npy) = va(0,npy-1)
endif
if( is==1 .and. .not. nested ) then
do j=js-1,je+1
- uc(0,j) = c1*utmp(-2,j) + c2*utmp(-1,j) + c3*utmp(0,j)
+ uc(0,j) = c1*utmp(-2,j) + c2*utmp(-1,j) + c3*utmp(0,j)
ut(1,j) = edge_interpolate4(ua(-1:2,j), dxa(-1:2,j))
!Want to use the UPSTREAM value
if (ut(1,j) > 0.) then
@@ -3074,14 +3034,14 @@ subroutine d2a2c_vect(u, v, ua, va, uc, vc, ut, vt, dord4, gridstruct, &
if( (ie+1)==npx .and. .not. nested ) then
do j=js-1,je+1
- uc(npx-1,j) = c1*utmp(npx-3,j)+c2*utmp(npx-2,j)+c3*utmp(npx-1,j)
+ uc(npx-1,j) = c1*utmp(npx-3,j)+c2*utmp(npx-2,j)+c3*utmp(npx-1,j)
ut(npx, j) = edge_interpolate4(ua(npx-2:npx+1,j), dxa(npx-2:npx+1,j))
if (ut(npx,j) > 0.) then
uc(npx,j) = ut(npx,j)*sin_sg(npx-1,j,3)
else
uc(npx,j) = ut(npx,j)*sin_sg(npx,j,1)
end if
- uc(npx+1,j) = c3*utmp(npx,j) + c2*utmp(npx+1,j) + c1*utmp(npx+2,j)
+ uc(npx+1,j) = c3*utmp(npx,j) + c2*utmp(npx+1,j) + c1*utmp(npx+2,j)
ut(npx-1,j) = (uc(npx-1,j)-v(npx-1,j)*cosa_u(npx-1,j))*rsin_u(npx-1,j)
ut(npx+1,j) = (uc(npx+1,j)-v(npx+1,j)*cosa_u(npx+1,j))*rsin_u(npx+1,j)
enddo
@@ -3180,7 +3140,7 @@ subroutine d2a2c_vect(u, v, ua, va, uc, vc, ut, vt, dord4, gridstruct, &
end subroutine d2a2c_vect
-
+
real function edge_interpolate4(ua, dxa)
real, intent(in) :: ua(4)
diff --git a/model/tp_core.F90 b/model/tp_core.F90
index 479692dc6..6614eee62 100644
--- a/model/tp_core.F90
+++ b/model/tp_core.F90
@@ -1,21 +1,21 @@
!***********************************************************************
-!* GNU Lesser General Public License
+!* GNU Lesser General Public License
!*
!* This file is part of the FV3 dynamical core.
!*
-!* The FV3 dynamical core is free software: you can redistribute it
+!* The FV3 dynamical core is free software: you can redistribute it
!* and/or modify it under the terms of the
!* GNU Lesser General Public License as published by the
-!* Free Software Foundation, either version 3 of the License, or
+!* Free Software Foundation, either version 3 of the License, or
!* (at your option) any later version.
!*
-!* The FV3 dynamical core is distributed in the hope that it will be
-!* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
-!* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+!* The FV3 dynamical core is distributed in the hope that it will be
+!* useful, but WITHOUT ANY WARRANTY; without even the implied warranty
+!* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
!* See the GNU General Public License for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
-!* License along with the FV3 dynamical core.
+!* License along with the FV3 dynamical core.
!* If not, see .
!***********************************************************************
@@ -62,6 +62,7 @@ module tp_core_mod
real, parameter:: r3 = 1./3.
real, parameter:: near_zero = 1.E-25
real, parameter:: ppm_limiter = 2.0
+ real, parameter:: r12 = 1./12.
#ifdef WAVE_FORM
! Suresh & Huynh scheme 2.2 (purtabation form)
@@ -160,13 +161,13 @@ subroutine fv_tp_2d(q, crx, cry, npx, npy, hord, fx, fy, xfx, yfx, &
ord_ou = hord
if (.not. gridstruct%nested) call copy_corners(q, npx, npy, 2, gridstruct%nested, bd, &
- gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner)
+ gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner)
call yppm(fy2, q, cry, ord_in, isd,ied,isd,ied, js,je,jsd,jed, npx,npy, gridstruct%dya, gridstruct%nested, gridstruct%grid_type, lim_fac)
do j=js,je+1
do i=isd,ied
- fyy(i,j) = yfx(i,j) * fy2(i,j)
+ fyy(i,j) = yfx(i,j) * fy2(i,j)
enddo
enddo
do j=js,je
@@ -178,18 +179,18 @@ subroutine fv_tp_2d(q, crx, cry, npx, npy, hord, fx, fy, xfx, yfx, &
call xppm(fx, q_i, crx(is,js), ord_ou, is,ie,isd,ied, js,je,jsd,jed, npx,npy, gridstruct%dxa, gridstruct%nested, gridstruct%grid_type, lim_fac)
if (.not. gridstruct%nested) call copy_corners(q, npx, npy, 1, gridstruct%nested, bd, &
- gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner)
+ gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner)
call xppm(fx2, q, crx, ord_in, is,ie,isd,ied, jsd,jed,jsd,jed, npx,npy, gridstruct%dxa, gridstruct%nested, gridstruct%grid_type, lim_fac)
- do j=jsd,jed
- do i=is,ie+1
- fx1(i) = xfx(i,j) * fx2(i,j)
- enddo
- do i=is,ie
- q_j(i,j) = (q(i,j)*gridstruct%area(i,j) + fx1(i)-fx1(i+1))/ra_x(i,j)
- enddo
- enddo
+ do j=jsd,jed
+ do i=is,ie+1
+ fx1(i) = xfx(i,j) * fx2(i,j)
+ enddo
+ do i=is,ie
+ q_j(i,j) = (q(i,j)*gridstruct%area(i,j) + fx1(i)-fx1(i+1))/ra_x(i,j)
+ enddo
+ enddo
call yppm(fy, q_j, cry, ord_ou, is,ie,isd,ied, js,je,jsd,jed, npx, npy, gridstruct%dya, gridstruct%nested, gridstruct%grid_type, lim_fac)
@@ -317,7 +318,7 @@ subroutine copy_corners(q, npx, npy, dir, nested, bd, &
endif
endif
-
+
end subroutine copy_corners
subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy, dxa, nested, grid_type, lim_fac)
@@ -334,10 +335,10 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy,
!OUTPUT PARAMETERS:
real , INTENT(OUT) :: flux(is:ie+1,jfirst:jlast) !< Flux
! Local
- real, dimension(is-1:ie+1):: bl, br, b0
+ real, dimension(is-1:ie+1):: bl, br, b0, a4, da1
real:: q1(isd:ied)
real, dimension(is:ie+1):: fx0, fx1, xt1
- logical, dimension(is-1:ie+1):: smt5, smt6
+ logical, dimension(is-1:ie+1):: ext5, ext6, smt5, smt6
logical, dimension(is:ie+1):: hi5, hi6
real al(is-1:ie+2)
real dm(is-2:ie+2)
@@ -361,9 +362,9 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy,
q1(i) = q(i,j)
enddo
- if ( iord < 8 ) then
+ if ( iord < 7 ) then
! ord = 2: perfectly linear ppm scheme
-! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6
+! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6
do i=is1, ie3
al(i) = p1*(q1(i-1)+q1(i)) + p2*(q1(i-2)+q1(i+1))
@@ -406,7 +407,7 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy,
fx1(i) = (1.+c(i,j))*(bl(i) + c(i,j)*b0(i))
flux(i,j) = q1(i)
endif
- if (smt5(i-1).or.smt5(i)) flux(i,j) = flux(i,j) + fx1(i)
+ if (smt5(i-1).or.smt5(i)) flux(i,j) = flux(i,j) + fx1(i)
enddo
elseif ( mord==2 ) then ! perfectly linear scheme
@@ -416,10 +417,10 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy,
xt = c(i,j)
if ( xt > 0. ) then
qtmp = q1(i-1)
- flux(i,j) = qtmp + (1.-xt)*(al(i)-qtmp-xt*(al(i-1)+al(i)-(qtmp+qtmp)))
+ flux(i,j) = qtmp + (1.-xt)*(al(i)-qtmp-xt*(al(i-1)+al(i)-(qtmp+qtmp)))
else
qtmp = q1(i)
- flux(i,j) = qtmp + (1.+xt)*(al(i)-qtmp+xt*(al(i)+al(i+1)-(qtmp+qtmp)))
+ flux(i,j) = qtmp + (1.+xt)*(al(i)-qtmp+xt*(al(i)+al(i+1)-(qtmp+qtmp)))
endif
! x0 = sign(dim(xt, 0.), 1.)
! x1 = sign(dim(0., xt), 1.)
@@ -439,26 +440,19 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy,
smt6(i) = 3.*x0 < xt
enddo
do i=is,ie+1
- fx1(i) = 0.
xt1(i) = c(i,j)
- hi5(i) = smt5(i-1) .and. smt5(i) ! more diffusive
- hi6(i) = smt6(i-1) .or. smt6(i)
- enddo
- do i=is,ie+1
if ( xt1(i) > 0. ) then
- if ( hi6(i) ) then
- fx1(i) = br(i-1) - xt1(i)*b0(i-1)
- elseif ( hi5(i) ) then ! 2nd order, piece-wise linear
- fx1(i) = sign(min(abs(bl(i-1)),abs(br(i-1))), br(i-1))
- endif
- flux(i,j) = q1(i-1) + (1.-xt1(i))*fx1(i)
+ if ( smt5(i-1) .or. smt6(i) ) then
+ flux(i,j) = q1(i-1) + (1.-xt1(i))*(br(i-1) - xt1(i)*b0(i-1))
+ else
+ flux(i,j) = q1(i-1)
+ endif
else
- if ( hi6(i) ) then
- fx1(i) = bl(i) + xt1(i)*b0(i)
- elseif ( hi5(i) ) then ! 2nd order, piece-wise linear
- fx1(i) = sign(min(abs(bl(i)), abs(br(i))), bl(i))
- endif
- flux(i,j) = q1(i) + (1.+xt1(i))*fx1(i)
+ if ( smt6(i-1) .or. smt5(i) ) then
+ flux(i,j) = q1(i) + (1.+xt1(i))*(bl(i) + xt1(i)*b0(i))
+ else
+ flux(i,j) = q1(i)
+ endif
endif
enddo
@@ -481,25 +475,51 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy,
enddo
!DEC$ VECTOR ALWAYS
do i=is,ie+1
- if ( xt1(i) > 0. ) then
+ if ( xt1(i) > 0. ) then
fx1(i) = (1.-xt1(i))*(br(i-1) - xt1(i)*b0(i-1))
flux(i,j) = q1(i-1)
else
fx1(i) = (1.+xt1(i))*(bl(i) + xt1(i)*b0(i))
flux(i,j) = q1(i)
endif
- if ( hi5(i) ) flux(i,j) = flux(i,j) + fx1(i)
+ if ( hi5(i) ) flux(i,j) = flux(i,j) + fx1(i)
enddo
else
- if ( mord==5 ) then
+ if ( iord==5 ) then
do i=is-1,ie+1
bl(i) = al(i) - q1(i)
br(i) = al(i+1) - q1(i)
b0(i) = bl(i) + br(i)
smt5(i) = bl(i)*br(i) < 0.
enddo
+ elseif ( iord==-5 ) then
+ do i=is-1,ie+1
+ bl(i) = al(i) - q1(i)
+ br(i) = al(i+1) - q1(i)
+ b0(i) = bl(i) + br(i)
+ smt5(i) = bl(i)*br(i) < 0.
+ da1(i) = br(i) - bl(i)
+ a4(i) = -3.*b0(i)
+ enddo
+ do i=is-1,ie+1
+ if( abs(da1(i)) < -a4(i) ) then
+ if( q1(i)+0.25/a4(i)*da1(i)**2+a4(i)*r12 < 0. ) then
+ if( .not. smt5(i) ) then
+ br(i) = 0.
+ bl(i) = 0.
+ b0(i) = 0.
+ elseif( da1(i) > 0. ) then
+ br(i) = -2.*bl(i)
+ b0(i) = -bl(i)
+ else
+ bl(i) = -2.*br(i)
+ b0(i) = -br(i)
+ endif
+ endif
+ endif
+ enddo
else
do i=is-1,ie+1
bl(i) = al(i) - q1(i)
@@ -507,19 +527,6 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy,
b0(i) = bl(i) + br(i)
smt5(i) = 3.*abs(b0(i)) < abs(bl(i)-br(i))
enddo
-!WMP
-! fix edge issues
- if ( (.not.nested) .and. grid_type < 3) then
- if( is==1 ) then
- smt5(0) = bl(0)*br(0) < 0.
- smt5(1) = bl(1)*br(1) < 0.
- endif
- if( (ie+1)==npx ) then
- smt5(npx-1) = bl(npx-1)*br(npx-1) < 0.
- smt5(npx ) = bl(npx )*br(npx ) < 0.
- endif
- endif
-!WMP
endif
!DEC$ VECTOR ALWAYS
@@ -531,7 +538,7 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy,
fx1(i) = (1.+c(i,j))*(bl(i) + c(i,j)*b0(i))
flux(i,j) = q1(i)
endif
- if (smt5(i-1).or.smt5(i)) flux(i,j) = flux(i,j) + fx1(i)
+ if (smt5(i-1).or.smt5(i)) flux(i,j) = flux(i,j) + fx1(i)
enddo
endif
@@ -542,7 +549,7 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy,
! Monotonic constraints:
! ord = 8: PPM with Lin's PPM fast monotone constraint
! ord = 10: PPM with Lin's modification of Huynh 2nd constraint
-! ord = 13: 10 plus positive definite constraint
+! ord = 13: positive definite constraint
do i=is-2,ie+2
xt = 0.25*(q1(i+1) - q1(i-1))
@@ -559,14 +566,7 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy,
bl(i) = -sign(min(abs(xt), abs(al(i )-q1(i))), xt)
br(i) = sign(min(abs(xt), abs(al(i+1)-q1(i))), xt)
enddo
- elseif ( iord==11 ) then
-! This is emulation of 2nd van Leer scheme using PPM codes
- do i=is1, ie1
- xt = ppm_fac*dm(i)
- bl(i) = -sign(min(abs(xt), abs(al(i )-q1(i))), xt)
- br(i) = sign(min(abs(xt), abs(al(i+1)-q1(i))), xt)
- enddo
- else
+ elseif ( iord==10 ) then
do i=is1-2, ie1+1
dq(i) = 2.*(q1(i+1) - q1(i))
enddo
@@ -585,6 +585,41 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy,
bl(i) = min( max(0., pmp_1, lac_1), max(bl(i), min(0., pmp_1, lac_1)) )
endif
enddo
+ elseif ( iord==11 ) then
+! This is emulation of 2nd van Leer scheme using PPM codes
+ do i=is1, ie1
+ xt = ppm_fac*dm(i)
+ bl(i) = -sign(min(abs(xt), abs(al(i )-q1(i))), xt)
+ br(i) = sign(min(abs(xt), abs(al(i+1)-q1(i))), xt)
+ enddo
+ elseif ( iord==7 .or. iord==12 ) then ! positive definite (Lin & Rood 1996)
+ do i=is1, ie1
+ bl(i) = al(i) - q1(i)
+ br(i) = al(i+1) - q1(i)
+ a4(i) = -3.*(bl(i) + br(i))
+ da1(i) = br(i) - bl(i)
+ ext5(i) = br(i)*bl(i) > 0.
+ ext6(i) = abs(da1(i)) < -a4(i)
+ enddo
+ do i=is1, ie1
+ if( ext6(i) ) then
+ if( q1(i)+0.25/a4(i)*da1(i)**2+a4(i)*r12 < 0. ) then
+ if( ext5(i) ) then
+ br(i) = 0.
+ bl(i) = 0.
+ elseif( da1(i) > 0. ) then
+ br(i) = -2.*bl(i)
+ else
+ bl(i) = -2.*br(i)
+ endif
+ endif
+ endif
+ enddo
+ else
+ do i=is1, ie1
+ bl(i) = al(i ) - q1(i)
+ br(i) = al(i+1) - q1(i)
+ enddo
endif
! Positive definite constraint:
if(iord==9 .or. iord==13) call pert_ppm(ie1-is1+1, q1(is1), bl(is1), br(is1), 0)
@@ -631,13 +666,30 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy,
endif
- do i=is,ie+1
- if( c(i,j)>0. ) then
- flux(i,j) = q1(i-1) + (1.-c(i,j))*(br(i-1)-c(i,j)*(bl(i-1)+br(i-1)))
- else
- flux(i,j) = q1(i ) + (1.+c(i,j))*(bl(i )+c(i,j)*(bl(i)+br(i)))
- endif
- enddo
+ if ( iord==7 ) then
+ do i=is-1,ie+1
+ b0(i) = bl(i) + br(i)
+ smt5(i) = bl(i) * br(i) < 0.
+ enddo
+ do i=is,ie+1
+ if ( c(i,j) > 0. ) then
+ fx1(i) = (1.-c(i,j))*(br(i-1) - c(i,j)*b0(i-1))
+ flux(i,j) = q1(i-1)
+ else
+ fx1(i) = (1.+c(i,j))*(bl(i) + c(i,j)*b0(i))
+ flux(i,j) = q1(i)
+ endif
+ if ( smt5(i-1).or.smt5(i) ) flux(i,j) = flux(i,j) + fx1(i)
+ enddo
+ else
+ do i=is,ie+1
+ if( c(i,j)>0. ) then
+ flux(i,j) = q1(i-1) + (1.-c(i,j))*(br(i-1)-c(i,j)*(bl(i-1)+br(i-1)))
+ else
+ flux(i,j) = q1(i ) + (1.+c(i,j))*(bl(i )+c(i,j)*(bl(i)+br(i)))
+ endif
+ enddo
+ endif
666 continue
@@ -661,7 +713,7 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy
real:: al(ifirst:ilast,js-1:je+2)
real, dimension(ifirst:ilast,js-1:je+1):: bl, br, b0
real:: dq(ifirst:ilast,js-3:je+2)
- real, dimension(ifirst:ilast):: fx0, fx1, xt1
+ real, dimension(ifirst:ilast):: fx0, fx1, xt1, a4
logical, dimension(ifirst:ilast,js-1:je+1):: smt5, smt6
logical, dimension(ifirst:ilast):: hi5, hi6
real:: x0, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2, r1
@@ -679,7 +731,7 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy
mord = abs(jord)
-if ( jord < 8 ) then
+if ( jord < 7 ) then
do j=js1, je3
do i=ifirst,ilast
@@ -733,7 +785,7 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy
fx1(i) = (1.+c(i,j))*(bl(i,j) + c(i,j)*b0(i,j))
flux(i,j) = q(i,j)
endif
- if (smt5(i,j-1).or.smt5(i,j)) flux(i,j) = flux(i,j) + fx1(i)
+ if (smt5(i,j-1).or.smt5(i,j)) flux(i,j) = flux(i,j) + fx1(i)
enddo
enddo
@@ -761,7 +813,7 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy
bl(i,j) = al(i,j ) - q(i,j)
br(i,j) = al(i,j+1) - q(i,j)
b0(i,j) = bl(i,j) + br(i,j)
- x0 = abs(b0(i,j))
+ x0 = abs(b0(i,j))
xt = abs(bl(i,j)-br(i,j))
smt5(i,j) = x0 < xt
smt6(i,j) = 3.*x0 < xt
@@ -769,26 +821,21 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy
enddo
do j=js,je+1
do i=ifirst,ilast
- fx1(i) = 0.
xt1(i) = c(i,j)
- hi5(i) = smt5(i,j-1) .and. smt5(i,j)
- hi6(i) = smt6(i,j-1) .or. smt6(i,j)
enddo
do i=ifirst,ilast
if ( xt1(i) > 0. ) then
- if( hi6(i) ) then
- fx1(i) = br(i,j-1) - xt1(i)*b0(i,j-1)
- elseif ( hi5(i) ) then ! both up-downwind sides are noisy; 2nd order, piece-wise linear
- fx1(i) = sign(min(abs(bl(i,j-1)),abs(br(i,j-1))),br(i,j-1))
+ if( smt5(i,j-1) .or. smt6(i,j) ) then
+ flux(i,j) = q(i,j-1) + (1.-xt1(i))*(br(i,j-1) - xt1(i)*b0(i,j-1))
+ else
+ flux(i,j) = q(i,j-1)
endif
- flux(i,j) = q(i,j-1) + (1.-xt1(i))*fx1(i)
else
- if( hi6(i) ) then
- fx1(i) = bl(i,j) + xt1(i)*b0(i,j)
- elseif ( hi5(i) ) then ! both up-downwind sides are noisy; 2nd order, piece-wise linear
- fx1(i) = sign(min(abs(bl(i,j)),abs(br(i,j))), bl(i,j))
+ if( smt6(i,j-1) .or. smt5(i,j) ) then
+ flux(i,j) = q(i,j) + (1.+xt1(i))*(bl(i,j) + xt1(i)*b0(i,j))
+ else
+ flux(i,j) = q(i,j)
endif
- flux(i,j) = q(i,j) + (1.+xt1(i))*fx1(i)
endif
enddo
enddo
@@ -800,7 +847,7 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy
bl(i,j) = al(i,j ) - q(i,j)
br(i,j) = al(i,j+1) - q(i,j)
b0(i,j) = bl(i,j) + br(i,j)
- x0 = abs(b0(i,j))
+ x0 = abs(b0(i,j))
xt = abs(bl(i,j)-br(i,j))
smt5(i,j) = x0 < xt
smt6(i,j) = 3.*x0 < xt
@@ -822,19 +869,47 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy
fx1(i) = (1.+xt1(i))*(bl(i,j) + xt1(i)*b0(i,j))
flux(i,j) = q(i,j)
endif
- if ( hi5(i) ) flux(i,j) = flux(i,j) + fx1(i)
+ if ( hi5(i) ) flux(i,j) = flux(i,j) + fx1(i)
enddo
enddo
- else ! mord=5,6,7
- if ( mord==5 ) then
+ else ! mord=5,6
+ if ( jord==5 ) then
+ do j=js-1,je+1
+ do i=ifirst,ilast
+ bl(i,j) = al(i,j ) - q(i,j)
+ br(i,j) = al(i,j+1) - q(i,j)
+ b0(i,j) = bl(i,j) + br(i,j)
+ smt5(i,j) = bl(i,j)*br(i,j) < 0.
+ enddo
+ enddo
+ elseif ( jord==-5 ) then
do j=js-1,je+1
do i=ifirst,ilast
bl(i,j) = al(i,j ) - q(i,j)
br(i,j) = al(i,j+1) - q(i,j)
b0(i,j) = bl(i,j) + br(i,j)
+ xt1(i) = br(i,j) - bl(i,j)
+ a4(i) = -3.*b0(i,j)
smt5(i,j) = bl(i,j)*br(i,j) < 0.
enddo
+ do i=ifirst,ilast
+ if( abs(xt1(i)) < -a4(i) ) then
+ if( q(i,j)+0.25/a4(i)*xt1(i)**2+a4(i)*r12 < 0. ) then
+ if( .not. smt5(i,j) ) then
+ br(i,j) = 0.
+ bl(i,j) = 0.
+ b0(i,j) = 0.
+ elseif( xt1(i) > 0. ) then
+ br(i,j) = -2.*bl(i,j)
+ b0(i,j) = -bl(i,j)
+ else
+ bl(i,j) = -2.*br(i,j)
+ b0(i,j) = -br(i,j)
+ endif
+ endif
+ endif
+ enddo
enddo
else
do j=js-1,je+1
@@ -845,23 +920,6 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy
smt5(i,j) = 3.*abs(b0(i,j)) < abs(bl(i,j)-br(i,j))
enddo
enddo
-!WMP
-! fix edge issues
- if ( (.not.nested) .and. grid_type < 3) then
- if( js==1 ) then
- do i=ifirst,ilast
- smt5(i,0) = bl(i,0)*br(i,0) < 0.
- smt5(i,1) = bl(i,1)*br(i,1) < 0.
- enddo
- endif
- if( (je+1)==npy ) then
- do i=ifirst,ilast
- smt5(i,npy-1) = bl(i,npy-1)*br(i,npy-1) < 0.
- smt5(i,npy ) = bl(i,npy )*br(i,npy ) < 0.
- enddo
- endif
- endif
-!WMP
endif
do j=js,je+1
@@ -874,7 +932,7 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy
fx1(i) = (1.+c(i,j))*(bl(i,j) + c(i,j)*b0(i,j))
flux(i,j) = q(i,j)
endif
- if (smt5(i,j-1).or.smt5(i,j)) flux(i,j) = flux(i,j) + fx1(i)
+ if (smt5(i,j-1).or.smt5(i,j)) flux(i,j) = flux(i,j) + fx1(i)
enddo
enddo
@@ -885,7 +943,7 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy
! Monotonic constraints:
! ord = 8: PPM with Lin's PPM fast monotone constraint
! ord > 8: PPM with Lin's modification of Huynh 2nd constraint
-
+
do j=js-2,je+2
do i=ifirst,ilast
xt = 0.25*(q(i,j+1) - q(i,j-1))
@@ -907,15 +965,7 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy
br(i,j) = sign(min(abs(xt), abs(al(i,j+1)-q(i,j))), xt)
enddo
enddo
- elseif ( jord==11 ) then
- do j=js1,je1
- do i=ifirst,ilast
- xt = ppm_fac*dm(i,j)
- bl(i,j) = -sign(min(abs(xt), abs(al(i,j)-q(i,j))), xt)
- br(i,j) = sign(min(abs(xt), abs(al(i,j+1)-q(i,j))), xt)
- enddo
- enddo
- else
+ elseif ( jord==10 ) then
do j=js1-2,je1+1
do i=ifirst,ilast
dq(i,j) = 2.*(q(i,j+1) - q(i,j))
@@ -932,12 +982,52 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy
pmp_2 = dq(i,j-1)
lac_2 = pmp_2 - 0.75*dq(i,j-2)
br(i,j) = min(max(0.,pmp_2,lac_2), max(br(i,j), min(0.,pmp_2,lac_2)))
- pmp_1 = -dq(i,j)
+ pmp_1 = -dq(i,j)
lac_1 = pmp_1 + 0.75*dq(i,j+1)
bl(i,j) = min(max(0.,pmp_1,lac_1), max(bl(i,j), min(0.,pmp_1,lac_1)))
endif
enddo
enddo
+ elseif ( jord==11 ) then
+ do j=js1,je1
+ do i=ifirst,ilast
+ xt = ppm_fac*dm(i,j)
+ bl(i,j) = -sign(min(abs(xt), abs(al(i,j)-q(i,j))), xt)
+ br(i,j) = sign(min(abs(xt), abs(al(i,j+1)-q(i,j))), xt)
+ enddo
+ enddo
+ elseif ( jord==7 .or. jord==12 ) then
+ do j=js1,je1
+ do i=ifirst,ilast
+ bl(i,j) = al(i,j ) - q(i,j)
+ br(i,j) = al(i,j+1) - q(i,j)
+ xt1(i) = br(i,j) - bl(i,j)
+ a4(i) = -3.*(br(i,j) + bl(i,j))
+ hi5(i) = bl(i,j)*br(i,j) > 0.
+ hi6(i) = abs(xt1(i)) < -a4(i)
+ enddo
+ do i=ifirst,ilast
+ if( hi6(i) ) then
+ if( q(i,j)+0.25/a4(i)*xt1(i)**2+a4(i)*r12 < 0. ) then
+ if( hi5(i) ) then
+ br(i,j) = 0.
+ bl(i,j) = 0.
+ elseif( xt1(i) > 0. ) then
+ br(i,j) = -2.*bl(i,j)
+ else
+ bl(i,j) = -2.*br(i,j)
+ endif
+ endif
+ endif
+ enddo
+ enddo
+ else
+ do j=js1,je1
+ do i=ifirst,ilast
+ bl(i,j) = al(i,j ) - q(i,j)
+ br(i,j) = al(i,j+1) - q(i,j)
+ enddo
+ enddo
endif
if ( jord==9 .or. jord==13 ) then
! Positive definite constraint:
@@ -993,15 +1083,36 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy
endif
- do j=js,je+1
- do i=ifirst,ilast
- if( c(i,j)>0. ) then
- flux(i,j) = q(i,j-1) + (1.-c(i,j))*(br(i,j-1)-c(i,j)*(bl(i,j-1)+br(i,j-1)))
- else
- flux(i,j) = q(i,j ) + (1.+c(i,j))*(bl(i,j )+c(i,j)*(bl(i,j)+br(i,j)))
- endif
- enddo
- enddo
+ if ( jord==7 ) then
+ do j=js-1,je+1
+ do i=ifirst,ilast
+ b0(i,j) = bl(i,j) + br(i,j)
+ smt5(i,j) = bl(i,j) * br(i,j) < 0.
+ enddo
+ enddo
+ do j=js,je+1
+ do i=ifirst,ilast
+ if ( c(i,j) > 0. ) then
+ fx1(i) = (1.-c(i,j))*(br(i,j-1) - c(i,j)*b0(i,j-1))
+ flux(i,j) = q(i,j-1)
+ else
+ fx1(i) = (1.+c(i,j))*(bl(i,j) + c(i,j)*b0(i,j))
+ flux(i,j) = q(i,j)
+ endif
+ if ( smt5(i,j-1).or.smt5(i,j) ) flux(i,j) = flux(i,j) + fx1(i)
+ enddo
+ enddo
+ else
+ do j=js,je+1
+ do i=ifirst,ilast
+ if( c(i,j)>0. ) then
+ flux(i,j) = q(i,j-1) + (1.-c(i,j))*(br(i,j-1)-c(i,j)*(bl(i,j-1)+br(i,j-1)))
+ else
+ flux(i,j) = q(i,j ) + (1.+c(i,j))*(bl(i,j )+c(i,j)*(bl(i,j)+br(i,j)))
+ endif
+ enddo
+ enddo
+ endif
end subroutine yppm
@@ -1024,7 +1135,7 @@ subroutine mp_ghost_ew(im, jm, km, nq, ifirst, ilast, jfirst, jlast, &
!
! !DESCRIPTION:
!
-! Ghost 4d east/west
+! Ghost 4d east/west
!
! !REVISION HISTORY:
! 2005.08.22 Putman
@@ -1145,11 +1256,11 @@ subroutine deln_flux(nord,is,ie,js,je, npx, npy, damp, q, fx, fy, gridstruct, bd
#ifdef USE_SG
real, pointer, dimension(:,:) :: dx, dy, rdxc, rdyc
real, pointer, dimension(:,:,:) :: sin_sg
- dx => gridstruct%dx
- dy => gridstruct%dy
- rdxc => gridstruct%rdxc
- rdyc => gridstruct%rdyc
- sin_sg => gridstruct%sin_sg
+ dx => gridstruct%dx
+ dy => gridstruct%dy
+ rdxc => gridstruct%rdxc
+ rdyc => gridstruct%rdyc
+ sin_sg => gridstruct%sin_sg
#endif
i1 = is-1-nord; i2 = ie+1+nord
diff --git a/tools/fv_diagnostics.F90 b/tools/fv_diagnostics.F90
index a64040f91..b3e4f358b 100644
--- a/tools/fv_diagnostics.F90
+++ b/tools/fv_diagnostics.F90
@@ -3180,9 +3180,9 @@ subroutine range_check(qname, q, is, ie, js, je, n_g, km, pos, q_low, q_hi, bad_
enddo
enddo
enddo
- call mpp_error(NOTE,'==> Error from range_check: data out of bound')
+ !call mpp_error(NOTE,'==> Error from range_check: data out of bound')
endif
-106 format('Range_Warn: ',A,'(',i4,',',i4,',',i3,')=',f8.4,' at LON:',f8.4,' LAT:',f8.4)
+106 format('Range_Warn: ',A,'(',i4,',',i4,',',i3,')=',f10.4,' at LON:',f8.4,' LAT:',f8.4)
end subroutine range_check
@@ -3936,8 +3936,9 @@ subroutine helicity_relative(is, ie, js, je, ng, km, zvir, sphum, srh, &
srh(i,j) = 0.
below(i) = .true.
zh0(i) = 0.
+ k0 = -1
+ k1 = -1
-! if ( phis(i,j)/grav < 1.E3 ) then
do k=km,1,-1
if ( hydrostatic ) then
dz(i) = rdg*pt(i,j,k)*(1.+zvir*q(i,j,k,sphum))*(peln(i,k+1,j)-peln(i,k,j))
@@ -3966,6 +3967,7 @@ subroutine helicity_relative(is, ie, js, je, ng, km, zvir, sphum, srh, &
123 continue
! Lowest layer wind shear computed betw top edge and mid-layer
+ if ( (k0 .ne. -1) .and. (k1 .ne. -1) ) then
k = k1
srh(i,j) = 0.5*(va(i,j,k1)-vc(i))*(ua(i,j,k1-1)-ua(i,j,k1)) - &
0.5*(ua(i,j,k1)-uc(i))*(va(i,j,k1-1)-va(i,j,k1))
@@ -3973,7 +3975,7 @@ subroutine helicity_relative(is, ie, js, je, ng, km, zvir, sphum, srh, &
srh(i,j) = srh(i,j) + 0.5*(va(i,j,k)-vc(i))*(ua(i,j,k-1)-ua(i,j,k+1)) - &
0.5*(ua(i,j,k)-uc(i))*(va(i,j,k-1)-va(i,j,k+1))
enddo
-! endif
+ endif
enddo ! i-loop
enddo ! j-loop
@@ -4017,6 +4019,8 @@ subroutine helicity_relative_CAPS(is, ie, js, je, ng, km, zvir, sphum, srh, uc,
zh(i) = 0.
zh0 = 0.
below = .true.
+ k0 = -1
+ k1 = -1
do k=km,1,-1
if ( hydrostatic ) then
@@ -4042,6 +4046,7 @@ subroutine helicity_relative_CAPS(is, ie, js, je, ng, km, zvir, sphum, srh, uc,
123 continue
! Lowest layer wind shear computed betw top edge and mid-layer
+ if ( (k0 .ne. -1) .and. (k1 .ne. -1) ) then
k = k1
srh(i,j) = 0.5*(va(i,j,k1)-vc(i,j))*(ua(i,j,k1-1)-ua(i,j,k1)) - &
0.5*(ua(i,j,k1)-uc(i,j))*(va(i,j,k1-1)-va(i,j,k1))
@@ -4049,6 +4054,7 @@ subroutine helicity_relative_CAPS(is, ie, js, je, ng, km, zvir, sphum, srh, uc,
srh(i,j) = srh(i,j) + 0.5*(va(i,j,k)-vc(i,j))*(ua(i,j,k-1)-ua(i,j,k+1)) - &
0.5*(ua(i,j,k)-uc(i,j))*(va(i,j,k-1)-va(i,j,k+1))
enddo
+ endif
enddo ! i-loop
enddo ! j-loop