diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index 0ffee531f..cd1c0de12 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -297,16 +297,16 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, iep1 = ie + 1 jep1 = je + 1 +!$OMP parallel do default(none) shared(npz,dp_ref,ak,bk) + do k=1,npz + dp_ref(k) = (ak(k+1)-ak(k)) + (bk(k+1)-bk(k))*1.E5 + enddo + if ( .not.hydrostatic ) then rgrav = 1.0/grav k1k = akap / (1.-akap) ! rg/Cv=0.4 -!$OMP parallel do default(none) shared(npz,dp_ref,ak,bk) - do k=1,npz - dp_ref(k) = (ak(k+1)-ak(k)) + (bk(k+1)-bk(k))*1.E5 - enddo - !$OMP parallel do default(none) shared(isd,ied,jsd,jed,zs,phis,rgrav) do j=jsd,jed do i=isd,ied diff --git a/model/fv_control.F90 b/model/fv_control.F90 index c4e687ca3..4c6dec079 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -767,13 +767,13 @@ subroutine run_setup(Atm, dt_atmos, grids_on_this_pe, p_split) ! Define n_split if not in namelist if (ntiles==6) then #ifdef MAPL_MODE - dimx = ceiling(stretch_fac)*4.0*(npx-1) - if (.not. hydrostatic ) then + dimx = stretch_fac*4.0*(npx-1) + if (.not. hydrostatic) then ns0 = 6 if ( dimx >= 360 ) ns0 = 7 if ( dimx >= 1440 ) ns0 = 8 - if ( stretch_fac > 1.0 ) ns0=ns0+1 - endif + endif + if (stretch_fac > 1.0) ns0 = 7 #else dimx = 4.0*(npx-1) if ( hydrostatic ) then diff --git a/model/fv_mapz.F90 b/model/fv_mapz.F90 index 13becf3eb..9b075c3f3 100644 --- a/model/fv_mapz.F90 +++ b/model/fv_mapz.F90 @@ -108,6 +108,10 @@ module fv_mapz_mod real, parameter:: cp_vap = cp_vapor !< 1846. real, parameter:: tice = 273.16 + logical, parameter :: w_limiter = .true. + real, parameter :: w_max = 90. + real, parameter :: w_min = -60. + real(kind=8) :: E_Flux = 0. private @@ -207,7 +211,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & !----------------------------------------------------------------------- real(kind=8), dimension(is:ie,js:je):: te_2d, zsum0, zsum1 real, dimension(is:ie,js:je):: dpln - real, dimension(is:ie,km) :: q2, dp2 + real, dimension(is:ie,km) :: q2, dp2, w2 real, dimension(is:ie,km+1):: pe1, pe2, pk1, pk2, pn2, phis real, dimension(is:ie+1,km+1):: pe0, pe3 real, dimension(is:ie):: gz, cvm @@ -332,7 +336,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & !$OMP ak,bk,nq,isd,ied,jsd,jed,kord_tr,fill, adiabatic, & !$OMP hs,w,ws,kord_wz,rrg,kord_mt,consv,remap_option,gmao_remap) & !$OMP private(gz,cvm,bkh,dp2, & -!$OMP pe0,pe1,pe2,pe3,pk1,pk2,pn2,phis,q2,dpln,dlnp) +!$OMP pe0,pe1,pe2,pe3,pk1,pk2,pn2,phis,q2,w2,dpln,dlnp) do 1000 j=js,je+1 do k=1,km+1 @@ -572,7 +576,62 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, & delz(i,j,k) = -delz(i,j,k)*dp2(i,k) enddo enddo - endif + + !Fix excessive w - momentum conserving --- sjl + ! gz(:) used here as a temporary array + if ( w_limiter ) then + do k=1,km + do i=is,ie + w2(i,k) = w(i,j,k) + enddo + enddo + do k=1, km-1 + do i=is,ie + if ( w2(i,k) > w_max ) then + gz(i) = (w2(i,k)-w_max) * dp2(i,k) + w2(i,k ) = w_max + w2(i,k+1) = w2(i,k+1) + gz(i)/dp2(i,k+1) + !print*, ' W_LIMITER down: ', i,j,k, w2(i,k:k+1), w(i,j,k:k+1) + elseif ( w2(i,k) < w_min ) then + gz(i) = (w2(i,k)-w_min) * dp2(i,k) + w2(i,k ) = w_min + w2(i,k+1) = w2(i,k+1) + gz(i)/dp2(i,k+1) + !print*, ' W_LIMITER down: ', i,j,k, w2(i,k:k+1), w(i,j,k:k+1) + endif + enddo + enddo + do k=km, 2, -1 + do i=is,ie + if ( w2(i,k) > w_max ) then + gz(i) = (w2(i,k)-w_max) * dp2(i,k) + w2(i,k ) = w_max + w2(i,k-1) = w2(i,k-1) + gz(i)/dp2(i,k-1) + !print*, ' W_LIMITER up: ', i,j,k, w2(i,k-1:k), w(i,j,k-1:k) + elseif ( w2(i,k) < w_min ) then + gz(i) = (w2(i,k)-w_min) * dp2(i,k) + w2(i,k ) = w_min + w2(i,k-1) = w2(i,k-1) + gz(i)/dp2(i,k-1) + !print*, ' W_LIMITER up: ', i,j,k, w2(i,k-1:k), w(i,j,k-1:k) + endif + enddo + enddo + do i=is,ie + if (w2(i,1) > w_max*2. ) then + w2(i,1) = w_max*2 ! sink out of the top of the domain + !print*, ' W_LIMITER top limited: ', i,j,1, w2(i,1), w(i,j,1) + elseif (w2(i,1) < w_min*2. ) then + w2(i,1) = w_min*2. + !print*, ' W_LIMITER top limited: ', i,j,1, w2(i,1), w(i,j,1) + endif + enddo + do k=1,km + do i=is,ie + w(i,j,k) = w2(i,k) + enddo + enddo + endif + + endif endif !(j < je+1)