From cef6c901d60b9ef99c796891cf9cfbfcaa45958d Mon Sep 17 00:00:00 2001 From: William Putman Date: Sat, 1 Apr 2023 13:25:42 -0400 Subject: [PATCH 1/5] updted moist_import to reduce rst size, alt FV3 sponge and new FV3 diagnostics --- model/dyn_core.F90 | 2 +- model/fv_dynamics.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index cd1c0de12..a91b2d5b6 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -697,7 +697,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, nord_v(k)=0; endif d_con_k = 0. - elseif ( (k==2) .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 diff --git a/model/fv_dynamics.F90 b/model/fv_dynamics.F90 index 15c77f853..8e0a2e258 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, & - 150., 335., bad_range) + 130., 335., bad_range) if ( .not. hydrostatic ) then call range_check('W_dyn', w, is, ie, js, je, ng, npz, gridstruct%agrid, & -100., 100., bad_range) From 07b0eef0132a62d718e669bc983507c045fba654 Mon Sep 17 00:00:00 2001 From: William Putman Date: Tue, 4 Apr 2023 10:32:25 -0400 Subject: [PATCH 2/5] stretched grid and helicity updates --- model/fv_control.F90 | 2 +- model/fv_dynamics.F90 | 2 +- tools/fv_diagnostics.F90 | 12 +++++++++--- 3 files changed, 11 insertions(+), 5 deletions(-) diff --git a/model/fv_control.F90 b/model/fv_control.F90 index 4c6dec079..311a54d16 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -773,7 +773,7 @@ subroutine run_setup(Atm, dt_atmos, grids_on_this_pe, p_split) if ( dimx >= 360 ) ns0 = 7 if ( dimx >= 1440 ) ns0 = 8 endif - if (stretch_fac > 1.0) ns0 = 7 + if (stretch_fac > 1.0) ns0 = 6 #else dimx = 4.0*(npx-1) if ( hydrostatic ) then 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/tools/fv_diagnostics.F90 b/tools/fv_diagnostics.F90 index a64040f91..3c07b6cf0 100644 --- a/tools/fv_diagnostics.F90 +++ b/tools/fv_diagnostics.F90 @@ -3180,7 +3180,7 @@ 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) @@ -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 From eade38ac94b871d5f5c430eeb77844f750357669 Mon Sep 17 00:00:00 2001 From: William Putman Date: Fri, 7 Apr 2023 10:40:24 -0400 Subject: [PATCH 3/5] updated FV# code for advection schemes and stretched grid defaults --- model/dyn_core.F90 | 93 +++++------ model/sw_core.F90 | 374 +++++++++++++++++++----------------------- model/tp_core.F90 | 393 +++++++++++++++++++++++++++++---------------- 3 files changed, 459 insertions(+), 401 deletions(-) diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index a91b2d5b6..3d0b09374 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -234,6 +234,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, real, allocatable, dimension(:,:,:):: pem, heat_source ! Auto 1D & 2D arrays: real, dimension(bd%isd:bd%ied,bd%jsd:bd%jed):: ws3, z_rat + real:: d_ext(npz) real:: dp_ref(npz) real:: zs(bd%isd:bd%ied,bd%jsd:bd%jed) !< surface height (m) real:: p1d(bd%is:bd%ie) @@ -650,7 +651,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, !$OMP is,ie,js,je,isd,ied,jsd,jed,omga,delp,gridstruct,npx,npy, & !$OMP ng,zh,vt,ptc,pt,u,v,w,uc,vc,ua,va,divgd,mfx,mfy,cx,cy, & !$OMP crx,cry,xfx,yfx,q_con,zvir,sphum,nq,q,dt,bd,rdt,iep1,jep1, & -!$OMP heat_source,diss_est,dpx) & +!$OMP heat_source,diss_est,dpx,d_ext) & !$OMP private(nord_k, nord_w, nord_t, damp_w, damp_t, d2_divg, kfac, & !$OMP d_con_k,kgb, hord_m, hord_v, hord_t, hord_p, wk, heat_s,diss_e, z_rat) do k=1,npz @@ -684,6 +685,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, if ( npz==1 .or. flagstruct%n_sponge<0 ) then d2_divg = flagstruct%d2_bg + d_ext(k) = flagstruct%d_ext else ! Sponge layers with del-2 damping on divergence, vorticity, w, z, and air mass (delp). ! no special damping of potential temperature in sponge layers @@ -697,6 +699,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, nord_v(k)=0; endif d_con_k = 0. + d_ext(k) = 0.02 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 @@ -704,10 +707,12 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, nord_v(k)=0; endif d_con_k = 0. + d_ext(k) = 0.01 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. + d_ext(k) = 0.01 endif endif @@ -721,9 +726,12 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, endif !--- external mode divergence damping --- - if ( flagstruct%d_ext > 0. ) & + if ( d_ext(k) > 0. ) then call a2b_ord2(delp(isd,jsd,k), wk, gridstruct, npx, npy, is, & ie, js, je, ng, .false.) + else + wk(:,:) = 0.0 + endif if ( .not.hydrostatic .and. flagstruct%do_f3d ) then ! Correction factor for 3D Coriolis force @@ -757,13 +765,12 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, enddo endif - if ( flagstruct%d_ext > 0. ) then - do j=js,jep1 - do i=is,iep1 - ptc(i,j,k) = wk(i,j) ! delp at cell corners - enddo - enddo - endif + do j=js,jep1 + do i=is,iep1 + ptc(i,j,k) = wk(i,j) ! delp at cell corners + enddo + enddo + if ( flagstruct%d_con > 1.0E-5 .OR. flagstruct%do_skeb ) then ! if ( flagstruct%d_con > 1.0E-5 ) then ! Average horizontal "convergence" to cell center @@ -790,27 +797,23 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, #endif call timing_off('COMM_TOTAL') - if ( flagstruct%d_ext > 0. ) then - d2_divg = flagstruct%d_ext * gridstruct%da_min_c -!$OMP parallel do default(none) shared(is,iep1,js,jep1,npz,wk,ptc,divg2,vt,d2_divg) - do j=js,jep1 - do i=is,iep1 - wk(i,j) = ptc(i,j,1) - divg2(i,j) = wk(i,j)*vt(i,j,1) - enddo - do k=2,npz - do i=is,iep1 - wk(i,j) = wk(i,j) + ptc(i,j,k) - divg2(i,j) = divg2(i,j) + ptc(i,j,k)*vt(i,j,k) - enddo - enddo - do i=is,iep1 - divg2(i,j) = d2_divg*divg2(i,j)/wk(i,j) - enddo - enddo - else - divg2(:,:) = 0. - endif +!$OMP parallel do default(none) shared(is,iep1,js,jep1,npz,wk,ptc,divg2,vt,d_ext) + do j=js,jep1 + do i=is,iep1 + wk(i,j) = ptc(i,j,1) + divg2(i,j) = wk(i,j)*vt(i,j,1)*d_ext(1) + enddo + do k=2,npz + do i=is,iep1 + wk(i,j) = wk(i,j) + ptc(i,j,k) + divg2(i,j) = divg2(i,j) + ptc(i,j,k)*vt(i,j,k)*d_ext(k) + enddo + enddo + do i=is,iep1 + divg2(i,j) = divg2(i,j)/wk(i,j) + enddo + enddo + divg2=divg2*gridstruct%da_min_c call timing_on('COMM_TOTAL') call complete_group_halo_update(i_pack(1), domain) @@ -938,13 +941,13 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, if ( beta > 0. ) then call grad1_p_update(divg2, u, v, pkc, gz, dt, ng, gridstruct, bd, npx, npy, npz, ptop, beta_d, flagstruct%a2b_ord) else - call one_grad_p(u, v, pkc, gz, divg2, delp, dt, ng, gridstruct, bd, npx, npy, npz, ptop, hydrostatic, flagstruct%a2b_ord, flagstruct%d_ext) + call one_grad_p(u, v, pkc, gz, divg2, delp, dt, ng, gridstruct, bd, npx, npy, npz, ptop, hydrostatic, flagstruct%a2b_ord, d_ext) endif else if ( beta > 0. ) then call split_p_grad( u, v, pkc, gz, delp, pk3, beta_d, dt, ng, gridstruct, bd, npx, npy, npz, flagstruct%use_logp) elseif ( beta < -0.1 ) then - call one_grad_p(u, v, pkc, gz, divg2, delp, dt, ng, gridstruct, bd, npx, npy, npz, ptop, hydrostatic, flagstruct%a2b_ord, flagstruct%d_ext) + call one_grad_p(u, v, pkc, gz, divg2, delp, dt, ng, gridstruct, bd, npx, npy, npz, ptop, hydrostatic, flagstruct%a2b_ord, d_ext) else call nh_p_grad(u, v, pkc, gz, delp, pk3, dt, ng, gridstruct, bd, npx, npy, npz, flagstruct%use_logp) endif @@ -1732,7 +1735,7 @@ subroutine one_grad_p(u, v, pk, gz, divg2, delp, dt, ng, gridstruct, bd, npx, np ptop, hydrostatic, a2b_ord, d_ext) integer, intent(IN) :: ng, npx, npy, npz, a2b_ord -real, intent(IN) :: dt, ptop, d_ext +real, intent(IN) :: dt, ptop, d_ext(npz) logical, intent(in) :: hydrostatic type(fv_grid_bounds_type), intent(IN) :: bd real, intent(in) :: divg2(bd%is:bd%ie+1,bd%js:bd%je+1) @@ -1796,12 +1799,10 @@ subroutine one_grad_p(u, v, pk, gz, divg2, delp, dt, ng, gridstruct, bd, npx, np endif enddo -if ( d_ext > 0. ) then - !$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 @@ -1812,22 +1813,8 @@ subroutine one_grad_p(u, v, pk, gz, divg2, delp, dt, ng, gridstruct, bd, npx, np enddo enddo -else - - !$OMP parallel do default(none) shared(is,ie,js,je,wk1,wk2) - do j=js,je+1 - do i=is,ie - wk2(i,j) = 0. - enddo - do i=is,ie+1 - wk1(i,j) = 0. - enddo - enddo - -endif - !$OMP parallel do default(none) shared(is,ie,js,je,npz,pk,delp,hydrostatic,a2b_ord,gridstruct, & -!$OMP npx,npy,isd,jsd,ng,u,v,wk2,dt,gz,wk1) & +!$OMP npx,npy,isd,jsd,ng,u,v,wk2,dt,gz,wk1,d_ext) & !$OMP private(wk) do k=1,npz @@ -1847,14 +1834,14 @@ subroutine one_grad_p(u, v, pk, gz, divg2, delp, dt, ng, gridstruct, bd, npx, np do j=js,je+1 do i=is,ie - u(i,j,k) = gridstruct%rdx(i,j)*(wk2(i,j)+u(i,j,k) + dt/(wk(i,j)+wk(i+1,j)) * & + u(i,j,k) = gridstruct%rdx(i,j)*(d_ext(k)*wk2(i,j)+u(i,j,k) + dt/(wk(i,j)+wk(i+1,j)) * & ((gz(i,j,k+1)-gz(i+1,j,k))*(pk(i+1,j,k+1)-pk(i,j,k)) & + (gz(i,j,k)-gz(i+1,j,k+1))*(pk(i,j,k+1)-pk(i+1,j,k)))) enddo enddo do j=js,je do i=is,ie+1 - v(i,j,k) = gridstruct%rdy(i,j)*(wk1(i,j)+v(i,j,k) + dt/(wk(i,j)+wk(i,j+1)) * & + v(i,j,k) = gridstruct%rdy(i,j)*(d_ext(k)*wk1(i,j)+v(i,j,k) + dt/(wk(i,j)+wk(i,j+1)) * & ((gz(i,j,k+1)-gz(i,j+1,k))*(pk(i,j+1,k+1)-pk(i,j,k)) & + (gz(i,j,k)-gz(i,j+1,k+1))*(pk(i,j,k+1)-pk(i,j+1,k)))) enddo 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 From 8742ac9e6b13b902becec3b61986145c7f950cc5 Mon Sep 17 00:00:00 2001 From: William Putman Date: Fri, 7 Apr 2023 13:11:32 -0400 Subject: [PATCH 4/5] more conservative splitting for stretched grids --- model/fv_control.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/model/fv_control.F90 b/model/fv_control.F90 index 311a54d16..4c6dec079 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -773,7 +773,7 @@ subroutine run_setup(Atm, dt_atmos, grids_on_this_pe, p_split) if ( dimx >= 360 ) ns0 = 7 if ( dimx >= 1440 ) ns0 = 8 endif - if (stretch_fac > 1.0) ns0 = 6 + if (stretch_fac > 1.0) ns0 = 7 #else dimx = 4.0*(npx-1) if ( hydrostatic ) then From 42c2c2634457dbd59d50e0946a4512a59fd06bf4 Mon Sep 17 00:00:00 2001 From: William Putman Date: Sun, 9 Apr 2023 17:45:14 -0400 Subject: [PATCH 5/5] dz_min patch based on Zhou & Juang 2023, sponge layer tweaks --- model/dyn_core.F90 | 109 +++++++++++++++++++++------------------ model/fv_arrays.F90 | 2 + model/fv_control.F90 | 5 +- model/nh_utils.F90 | 13 +++-- tools/fv_diagnostics.F90 | 2 +- 5 files changed, 73 insertions(+), 58 deletions(-) diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index 3d0b09374..2ec8a7469 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -234,7 +234,6 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, real, allocatable, dimension(:,:,:):: pem, heat_source ! Auto 1D & 2D arrays: real, dimension(bd%isd:bd%ied,bd%jsd:bd%jed):: ws3, z_rat - real:: d_ext(npz) real:: dp_ref(npz) real:: zs(bd%isd:bd%ied,bd%jsd:bd%jed) !< surface height (m) real:: p1d(bd%is:bd%ie) @@ -558,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') @@ -651,7 +650,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, !$OMP is,ie,js,je,isd,ied,jsd,jed,omga,delp,gridstruct,npx,npy, & !$OMP ng,zh,vt,ptc,pt,u,v,w,uc,vc,ua,va,divgd,mfx,mfy,cx,cy, & !$OMP crx,cry,xfx,yfx,q_con,zvir,sphum,nq,q,dt,bd,rdt,iep1,jep1, & -!$OMP heat_source,diss_est,dpx,d_ext) & +!$OMP heat_source,diss_est,dpx) & !$OMP private(nord_k, nord_w, nord_t, damp_w, damp_t, d2_divg, kfac, & !$OMP d_con_k,kgb, hord_m, hord_v, hord_t, hord_p, wk, heat_s,diss_e, z_rat) do k=1,npz @@ -685,7 +684,6 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, if ( npz==1 .or. flagstruct%n_sponge<0 ) then d2_divg = flagstruct%d2_bg - d_ext(k) = flagstruct%d_ext else ! Sponge layers with del-2 damping on divergence, vorticity, w, z, and air mass (delp). ! no special damping of potential temperature in sponge layers @@ -697,22 +695,21 @@ 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. - d_ext(k) = 0.02 - 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. - d_ext(k) = 0.01 - 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. - d_ext(k) = 0.01 endif endif @@ -726,12 +723,9 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, endif !--- external mode divergence damping --- - if ( d_ext(k) > 0. ) then + if ( flagstruct%d_ext > 0. ) & call a2b_ord2(delp(isd,jsd,k), wk, gridstruct, npx, npy, is, & ie, js, je, ng, .false.) - else - wk(:,:) = 0.0 - endif if ( .not.hydrostatic .and. flagstruct%do_f3d ) then ! Correction factor for 3D Coriolis force @@ -765,12 +759,13 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, enddo endif - do j=js,jep1 - do i=is,iep1 - ptc(i,j,k) = wk(i,j) ! delp at cell corners - enddo - enddo - + if ( flagstruct%d_ext > 0. ) then + do j=js,jep1 + do i=is,iep1 + ptc(i,j,k) = wk(i,j) ! delp at cell corners + enddo + enddo + endif if ( flagstruct%d_con > 1.0E-5 .OR. flagstruct%do_skeb ) then ! if ( flagstruct%d_con > 1.0E-5 ) then ! Average horizontal "convergence" to cell center @@ -797,23 +792,27 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, #endif call timing_off('COMM_TOTAL') -!$OMP parallel do default(none) shared(is,iep1,js,jep1,npz,wk,ptc,divg2,vt,d_ext) - do j=js,jep1 - do i=is,iep1 - wk(i,j) = ptc(i,j,1) - divg2(i,j) = wk(i,j)*vt(i,j,1)*d_ext(1) - enddo - do k=2,npz - do i=is,iep1 - wk(i,j) = wk(i,j) + ptc(i,j,k) - divg2(i,j) = divg2(i,j) + ptc(i,j,k)*vt(i,j,k)*d_ext(k) - enddo - enddo - do i=is,iep1 - divg2(i,j) = divg2(i,j)/wk(i,j) - enddo - enddo - divg2=divg2*gridstruct%da_min_c + if ( flagstruct%d_ext > 0. ) then + d2_divg = flagstruct%d_ext * gridstruct%da_min_c +!$OMP parallel do default(none) shared(is,iep1,js,jep1,npz,wk,ptc,divg2,vt,d2_divg) + do j=js,jep1 + do i=is,iep1 + wk(i,j) = ptc(i,j,1) + divg2(i,j) = wk(i,j)*vt(i,j,1) + enddo + do k=2,npz + do i=is,iep1 + wk(i,j) = wk(i,j) + ptc(i,j,k) + divg2(i,j) = divg2(i,j) + ptc(i,j,k)*vt(i,j,k) + enddo + enddo + do i=is,iep1 + divg2(i,j) = d2_divg*divg2(i,j)/wk(i,j) + enddo + enddo + else + divg2(:,:) = 0. + endif call timing_on('COMM_TOTAL') call complete_group_halo_update(i_pack(1), domain) @@ -846,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 @@ -941,13 +940,13 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap, if ( beta > 0. ) then call grad1_p_update(divg2, u, v, pkc, gz, dt, ng, gridstruct, bd, npx, npy, npz, ptop, beta_d, flagstruct%a2b_ord) else - call one_grad_p(u, v, pkc, gz, divg2, delp, dt, ng, gridstruct, bd, npx, npy, npz, ptop, hydrostatic, flagstruct%a2b_ord, d_ext) + call one_grad_p(u, v, pkc, gz, divg2, delp, dt, ng, gridstruct, bd, npx, npy, npz, ptop, hydrostatic, flagstruct%a2b_ord, flagstruct%d_ext) endif else if ( beta > 0. ) then call split_p_grad( u, v, pkc, gz, delp, pk3, beta_d, dt, ng, gridstruct, bd, npx, npy, npz, flagstruct%use_logp) elseif ( beta < -0.1 ) then - call one_grad_p(u, v, pkc, gz, divg2, delp, dt, ng, gridstruct, bd, npx, npy, npz, ptop, hydrostatic, flagstruct%a2b_ord, d_ext) + call one_grad_p(u, v, pkc, gz, divg2, delp, dt, ng, gridstruct, bd, npx, npy, npz, ptop, hydrostatic, flagstruct%a2b_ord, flagstruct%d_ext) else call nh_p_grad(u, v, pkc, gz, delp, pk3, dt, ng, gridstruct, bd, npx, npy, npz, flagstruct%use_logp) endif @@ -1168,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 @@ -1735,7 +1730,7 @@ subroutine one_grad_p(u, v, pk, gz, divg2, delp, dt, ng, gridstruct, bd, npx, np ptop, hydrostatic, a2b_ord, d_ext) integer, intent(IN) :: ng, npx, npy, npz, a2b_ord -real, intent(IN) :: dt, ptop, d_ext(npz) +real, intent(IN) :: dt, ptop, d_ext logical, intent(in) :: hydrostatic type(fv_grid_bounds_type), intent(IN) :: bd real, intent(in) :: divg2(bd%is:bd%ie+1,bd%js:bd%je+1) @@ -1799,6 +1794,8 @@ subroutine one_grad_p(u, v, pk, gz, divg2, delp, dt, ng, gridstruct, bd, npx, np endif enddo +if ( d_ext > 0. ) then + !$OMP parallel do default(none) shared(is,ie,js,je,wk2,divg2) do j=js,je+1 do i=is,ie @@ -1813,8 +1810,22 @@ subroutine one_grad_p(u, v, pk, gz, divg2, delp, dt, ng, gridstruct, bd, npx, np enddo enddo +else + + !$OMP parallel do default(none) shared(is,ie,js,je,wk1,wk2) + do j=js,je+1 + do i=is,ie + wk2(i,j) = 0. + enddo + do i=is,ie+1 + wk1(i,j) = 0. + enddo + enddo + +endif + !$OMP parallel do default(none) shared(is,ie,js,je,npz,pk,delp,hydrostatic,a2b_ord,gridstruct, & -!$OMP npx,npy,isd,jsd,ng,u,v,wk2,dt,gz,wk1,d_ext) & +!$OMP npx,npy,isd,jsd,ng,u,v,wk2,dt,gz,wk1) & !$OMP private(wk) do k=1,npz @@ -1834,14 +1845,14 @@ subroutine one_grad_p(u, v, pk, gz, divg2, delp, dt, ng, gridstruct, bd, npx, np do j=js,je+1 do i=is,ie - u(i,j,k) = gridstruct%rdx(i,j)*(d_ext(k)*wk2(i,j)+u(i,j,k) + dt/(wk(i,j)+wk(i+1,j)) * & + u(i,j,k) = gridstruct%rdx(i,j)*(wk2(i,j)+u(i,j,k) + dt/(wk(i,j)+wk(i+1,j)) * & ((gz(i,j,k+1)-gz(i+1,j,k))*(pk(i+1,j,k+1)-pk(i,j,k)) & + (gz(i,j,k)-gz(i+1,j,k+1))*(pk(i,j,k+1)-pk(i+1,j,k)))) enddo enddo do j=js,je do i=is,ie+1 - v(i,j,k) = gridstruct%rdy(i,j)*(d_ext(k)*wk1(i,j)+v(i,j,k) + dt/(wk(i,j)+wk(i,j+1)) * & + v(i,j,k) = gridstruct%rdy(i,j)*(wk1(i,j)+v(i,j,k) + dt/(wk(i,j)+wk(i,j+1)) * & ((gz(i,j,k+1)-gz(i,j+1,k))*(pk(i,j+1,k+1)-pk(i,j,k)) & + (gz(i,j,k)-gz(i,j+1,k+1))*(pk(i,j,k+1)-pk(i,j+1,k)))) 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/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/tools/fv_diagnostics.F90 b/tools/fv_diagnostics.F90 index 3c07b6cf0..b3e4f358b 100644 --- a/tools/fv_diagnostics.F90 +++ b/tools/fv_diagnostics.F90 @@ -3182,7 +3182,7 @@ subroutine range_check(qname, q, is, ie, js, je, n_g, km, pos, q_low, q_hi, bad_ enddo !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