From 13c653fef3b9bbb75e6b8ed5824c8f5bc5f69e78 Mon Sep 17 00:00:00 2001 From: Xiaqiong Zhou Date: Thu, 10 Sep 2020 16:00:30 +0000 Subject: [PATCH 1/6] Add a positive-definite advection option; the reconstructed cell-interface value is not forced to be positive --- model/tp_core.F90 | 311 ++++++++++++++++++++++++---------------------- 1 file changed, 162 insertions(+), 149 deletions(-) diff --git a/model/tp_core.F90 b/model/tp_core.F90 index c5cc3b29e..6f8e810de 100644 --- a/model/tp_core.F90 +++ b/model/tp_core.F90 @@ -61,6 +61,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) @@ -336,15 +337,14 @@ 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 + real, dimension(is:ie+1):: fx0, fx1 logical, dimension(is-1:ie+1):: smt5, smt6 - logical, dimension(is:ie+1):: hi5, hi6 real al(is-1:ie+2) real dm(is-2:ie+2) real dq(is-3:ie+2) - integer:: i, j, ie3, is1, ie1, mord + integer:: i, j, ie3, is1, ie1 real:: x0, x1, xt, qtmp, pmp_1, lac_1, pmp_2, lac_2 if ( .not. bounded_domain .and. grid_type<3 ) then @@ -354,7 +354,6 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy, is1 = is-1; ie3 = ie+2 ie1 = ie+1 end if - mord = abs(iord) do 666 j=jfirst,jlast @@ -362,13 +361,18 @@ 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 < 8 ) 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)) enddo + if ( iord==7 ) then + do i=is1, ie3 + if ( al(i)<0. ) al(i) = 0.5*(q1(i-1)+q1(i)) + enddo + endif if ( .not. bounded_domain .and. grid_type<3 ) then if ( is==1 ) then @@ -376,51 +380,37 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy, al(1) = 0.5*(((2.*dxa(0,j)+dxa(-1,j))*q1(0)-dxa(0,j)*q1(-1))/(dxa(-1,j)+dxa(0,j)) & + ((2.*dxa(1,j)+dxa( 2,j))*q1(1)-dxa(1,j)*q1( 2))/(dxa(1, j)+dxa(2,j))) al(2) = c3*q1(1) + c2*q1(2) +c1*q1(3) + if(iord==7) then + al(0) = max(0., al(0)) + al(1) = max(0., al(1)) + al(2) = max(0., al(2)) + endif endif if ( (ie+1)==npx ) then al(npx-1) = c1*q1(npx-3) + c2*q1(npx-2) + c3*q1(npx-1) al(npx) = 0.5*(((2.*dxa(npx-1,j)+dxa(npx-2,j))*q1(npx-1)-dxa(npx-1,j)*q1(npx-2))/(dxa(npx-2,j)+dxa(npx-1,j)) & + ((2.*dxa(npx, j)+dxa(npx+1,j))*q1(npx )-dxa(npx, j)*q1(npx+1))/(dxa(npx, j)+dxa(npx+1,j))) al(npx+1) = c3*q1(npx) + c2*q1(npx+1) + c1*q1(npx+2) + if(iord==7) then + al(npx-1) = max(0., al(npx-1)) + al(npx ) = max(0., al(npx )) + al(npx+1) = max(0., al(npx+1)) + endif endif endif - if ( iord<0 ) then - do i=is-1, ie+2 - al(i) = max(0., al(i)) - enddo - endif - - if ( mord==1 ) then ! perfectly linear scheme - 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) = abs(lim_fac*b0(i)) < abs(bl(i)-br(i)) - enddo -!DEC$ VECTOR ALWAYS - 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 - - elseif ( mord==2 ) then ! perfectly linear scheme + if ( iord==2 ) then ! perfectly linear scheme +! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6 < ord7 !DEC$ VECTOR ALWAYS do i=is,ie+1 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.) @@ -428,8 +418,7 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy, ! + x1*(q1(i) +(1.+xt)*(al(i)-qtmp+xt*(al(i)+al(i+1)-(qtmp+qtmp)))) enddo - elseif ( mord==3 ) then - + elseif ( iord==3 ) then do i=is-1,ie+1 bl(i) = al(i) - q1(i) br(i) = al(i+1) - q1(i) @@ -441,30 +430,27 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy, 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 + xt = c(i,j) + if ( xt > 0. ) then + fx0(i) = q1(i-1) + if ( smt6(i-1).or.smt5(i) ) then + fx1(i) = br(i-1) - xt*b0(i-1) + elseif ( smt5(i-1) ) 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) else - if ( hi6(i) ) then - fx1(i) = bl(i) + xt1(i)*b0(i) - elseif ( hi5(i) ) then ! 2nd order, piece-wise linear + fx0(i) = q1(i) + if ( smt6(i).or.smt5(i-1) ) then + fx1(i) = bl(i) + xt*b0(i) + elseif ( smt5(i) ) then fx1(i) = sign(min(abs(bl(i)), abs(br(i))), bl(i)) endif - flux(i,j) = q1(i) + (1.+xt1(i))*fx1(i) endif + flux(i,j) = fx0(i) + (1.-abs(xt))*fx1(i) enddo - - elseif ( mord==4 ) then - + elseif ( iord==4 ) then do i=is-1,ie+1 bl(i) = al(i) - q1(i) br(i) = al(i+1) - q1(i) @@ -475,53 +461,74 @@ 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 - xt1(i) = c(i,j) - hi5(i) = smt5(i-1) .and. smt5(i) ! more diffusive - hi6(i) = smt6(i-1) .or. smt6(i) - hi5(i) = hi5(i) .or. hi6(i) + fx1(i) = 0. enddo !DEC$ VECTOR ALWAYS do i=is,ie+1 - if ( xt1(i) > 0. ) then - fx1(i) = (1.-xt1(i))*(br(i-1) - xt1(i)*b0(i-1)) - flux(i,j) = q1(i-1) + if ( c(i,j) > 0. ) then + fx0(i) = q1(i-1) + if ( smt6(i-1).or.smt5(i) ) fx1(i) = (1.-c(i,j))*(br(i-1) - c(i,j)*b0(i-1)) else - fx1(i) = (1.+xt1(i))*(bl(i) + xt1(i)*b0(i)) - flux(i,j) = q1(i) + fx0(i) = q1(i) + if ( smt6(i).or.smt5(i-1) ) fx1(i) = (1.+c(i,j))*(bl(i) + c(i,j)*b0(i)) endif - if ( hi5(i) ) flux(i,j) = flux(i,j) + fx1(i) + flux(i,j) = fx0(i) + fx1(i) enddo - else - - if ( mord==5 ) then +! iord = 5 & 6 + 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 +! A positive-definite scheme by S-J Lin; Harris et al (2020) JAMES + 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) br(i) = al(i+1) - q1(i) b0(i) = bl(i) + br(i) - smt5(i) = 3.*abs(b0(i)) < abs(bl(i)-br(i)) + smt5(i) = abs(3.*b0(i)) < abs(bl(i)-br(i)) enddo endif - !DEC$ VECTOR ALWAYS 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) + 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) + 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 goto 666 @@ -648,11 +655,10 @@ 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 - integer:: i, j, js1, je3, je1, mord + integer:: i, j, js1, je3, je1 if ( .not.bounded_domain .and. grid_type < 3 ) then ! Cubed-sphere: @@ -664,8 +670,6 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy je1 = je+1 endif - mord = abs(jord) - if ( jord < 8 ) then do j=js1, je3 @@ -673,6 +677,13 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy al(i,j) = p1*(q(i,j-1)+q(i,j)) + p2*(q(i,j-2)+q(i,j+1)) enddo enddo + if ( jord==7 ) then + do j=js1, je3 + do i=ifirst,ilast + if ( al(i,j)<0. ) al(i,j) = 0.5*(q(i,j)+q(i,j+1)) + enddo + enddo + endif if ( .not. bounded_domain .and. grid_type<3 ) then if( js==1 ) then @@ -682,6 +693,13 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy + ((2.*dya(i,1)+dya(i,2))*q(i,1)-dya(i,1)*q(i,2))/(dya(i,1)+dya(i,2))) al(i,2) = c3*q(i,1) + c2*q(i,2) + c1*q(i,3) enddo + if ( jord==7 ) then + do i=ifirst,ilast + al(i,0) = max(0., al(i,0)) + al(i,1) = max(0., al(i,1)) + al(i,2) = max(0., al(i,2)) + enddo + endif endif if( (je+1)==npy ) then do i=ifirst,ilast @@ -690,41 +708,17 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy + ((2.*dya(i,npy)+dya(i,npy+1))*q(i,npy)-dya(i,npy)*q(i,npy+1))/(dya(i,npy)+dya(i,npy+1))) al(i,npy+1) = c3*q(i,npy) + c2*q(i,npy+1) + c1*q(i,npy+2) enddo + if (jord==7 ) then + do i=ifirst,ilast + al(i,npy-1) = max(0., al(i,npy-1)) + al(i,npy ) = max(0., al(i,npy )) + al(i,npy+1) = max(0., al(i,npy+1)) + enddo + endif endif endif - if ( jord<0 ) then - do j=js-1, je+2 - do i=ifirst,ilast - al(i,j) = max(0., al(i,j)) - enddo - enddo - endif - - if ( mord==1 ) 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) = abs(lim_fac*b0(i,j)) < abs(bl(i,j)-br(i,j)) - enddo - enddo - do j=js,je+1 -!DEC$ VECTOR ALWAYS - 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 - - elseif ( mord==2 ) then ! Perfectly linear scheme + if ( jord==2 ) then ! Perfectly linear scheme ! Diffusivity: ord2 < ord5 < ord3 < ord4 < ord6 < ord7 do j=js,je+1 @@ -741,14 +735,13 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy enddo enddo - elseif ( mord==3 ) then - + elseif ( jord==3 ) 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) - 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 @@ -757,37 +750,35 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy 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 + xt = c(i,j) + if ( xt > 0. ) then + fx0(i) = q(i,j-1) + if( smt6(i,j-1).or.smt5(i,j) ) then + fx1(i) = br(i,j-1) - xt*b0(i,j-1) + elseif ( smt5(i,j-1) ) 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)) 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 + fx0(i) = q(i,j) + if( smt6(i,j).or.smt5(i,j-1) ) then + fx1(i) = bl(i,j) + xt*b0(i,j) + elseif ( smt5(i,j) ) then fx1(i) = sign(min(abs(bl(i,j)),abs(br(i,j))), bl(i,j)) endif - flux(i,j) = q(i,j) + (1.+xt1(i))*fx1(i) endif + flux(i,j) = fx0(i) + (1.-abs(xt))*fx1(i) enddo enddo - elseif ( mord==4 ) then - + elseif ( jord==4 ) 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) - 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 @@ -795,33 +786,58 @@ 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 - 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) - hi5(i) = hi5(i) .or. hi6(i) + fx1(i) = 0. enddo !DEC$ VECTOR ALWAYS do i=ifirst,ilast - if ( xt1(i) > 0. ) then - fx1(i) = (1.-xt1(i))*(br(i,j-1) - xt1(i)*b0(i,j-1)) - flux(i,j) = q(i,j-1) - else - 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 ( c(i,j) > 0. ) then + fx0(i) = q(i,j-1) + if( smt6(i,j-1).or.smt5(i,j) ) fx1(i) = (1.-c(i,j))*(br(i,j-1) - c(i,j)*b0(i,j-1)) + else + fx0(i) = q(i,j) + if( smt6(i,j).or.smt5(i,j-1) ) fx1(i) = (1.+c(i,j))*(bl(i,j) + c(i,j)*b0(i,j)) + endif + flux(i,j) = fx0(i) + fx1(i) enddo enddo - else ! mord=5,6,7 - if ( mord==5 ) then + else ! jord=5,6,7 + 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 @@ -829,11 +845,10 @@ 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) - smt5(i,j) = 3.*abs(b0(i,j)) < abs(bl(i,j)-br(i,j)) + smt5(i,j) = abs(3.*b0(i,j)) < abs(bl(i,j)-br(i,j)) enddo enddo endif - do j=js,je+1 !DEC$ VECTOR ALWAYS do i=ifirst,ilast @@ -844,10 +859,9 @@ 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 - endif return @@ -855,7 +869,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)) @@ -902,7 +916,7 @@ 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 @@ -1034,7 +1048,6 @@ subroutine pert_ppm(im, a0, al, ar, iv) ! Local: real a4, da1, da2, a6da, fmin integer i - real, parameter:: r12 = 1./12. !----------------------------------- ! Optimized PPM in perturbation form: From 272233afa871cd760d6ed5ca1f51b301b204374b Mon Sep 17 00:00:00 2001 From: "Chan-hoo.Jeon" Date: Wed, 23 Sep 2020 11:29:07 +0000 Subject: [PATCH 2/6] Modify the halo extents of u and vt in the regional_boundary_update call that caused a reproducibility issue. --- model/fv_regional_bc.F90 | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index 98285b7a3..e3b633256 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -4409,6 +4409,12 @@ subroutine regional_boundary_update(array & j1=jsd j2=jed ! +! CHJ --- s --- + if(trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='u')then + j2=jed+1 + endif +! CHJ --- e --- + i1=isd i2=is-1 ! @@ -4453,7 +4459,13 @@ subroutine regional_boundary_update(array & ! j1=jsd j2=jed -! + +! CHJ --- s --- + if(trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='u')then + j2=jed+1 + endif +! CHJ --- e --- + i1=ie+1 i2=ied if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v')then From 6f3d2034bacaadd4967feaf527b2884354290ce6 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 7 Oct 2020 08:51:57 -0600 Subject: [PATCH 3/6] model/fv_dynamics.F90: bugfix when debugging output is enabled --- model/fv_dynamics.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/model/fv_dynamics.F90 b/model/fv_dynamics.F90 index 16c0e5753..1b706681b 100644 --- a/model/fv_dynamics.F90 +++ b/model/fv_dynamics.F90 @@ -774,12 +774,12 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, if ( flagstruct%fv_debug ) then if (is_master()) write(*,'(A, I3, A1, I3)') 'finished k_split ', n_map, '/', k_split call prt_mxm('T_dyn_a3', pt, is, ie, js, je, ng, npz, 1., gridstruct%area_64, domain) - call prt_mxm('SPHUM_dyn', q(isd,jsd,1,sphum ), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain) - call prt_mxm('liq_wat_dyn', q(isd,jsd,1,liq_wat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain) - call prt_mxm('rainwat_dyn', q(isd,jsd,1,rainwat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain) - call prt_mxm('ice_wat_dyn', q(isd,jsd,1,ice_wat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain) - call prt_mxm('snowwat_dyn', q(isd,jsd,1,snowwat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain) - call prt_mxm('graupel_dyn', q(isd,jsd,1,graupel), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain) + if (sphum > 0) call prt_mxm('SPHUM_dyn', q(isd,jsd,1,sphum ), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain) + if (liq_wat > 0) call prt_mxm('liq_wat_dyn', q(isd,jsd,1,liq_wat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain) + if (rainwat > 0) call prt_mxm('rainwat_dyn', q(isd,jsd,1,rainwat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain) + if (ice_wat > 0) call prt_mxm('ice_wat_dyn', q(isd,jsd,1,ice_wat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain) + if (snowwat > 0) call prt_mxm('snowwat_dyn', q(isd,jsd,1,snowwat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain) + if (graupel > 0) call prt_mxm('graupel_dyn', q(isd,jsd,1,graupel), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain) endif #ifdef AVEC_TIMERS call avec_timer_stop(6) From c3b564f163ff97b2911735c8686dba7f1606f086 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Sat, 10 Oct 2020 15:05:46 -0600 Subject: [PATCH 4/6] Comment out test for ntrac > ntracers in tools/external_ic.F90 --- tools/external_ic.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tools/external_ic.F90 b/tools/external_ic.F90 index 432703501..44f76c22e 100644 --- a/tools/external_ic.F90 +++ b/tools/external_ic.F90 @@ -570,8 +570,10 @@ subroutine get_nggps_ic (Atm, fv_domain, dt_atmos ) !--- read in the number of tracers in the NCEP NGGPS ICs call read_data ('INPUT/'//trim(fn_gfs_ctl), 'ntrac', ntrac, no_domain=.TRUE.) - if (ntrac > ntracers) call mpp_error(FATAL,'==> External_ic::get_nggps_ic: more NGGPS tracers & - &than defined in field_table '//trim(fn_gfs_ctl)//' for NGGPS IC') + ! DH* 20200922 - this breaks Ferrier-Aligo MP runs + !if (ntrac > ntracers) call mpp_error(FATAL,'==> External_ic::get_nggps_ic: more NGGPS tracers & + ! &than defined in field_table '//trim(fn_gfs_ctl)//' for NGGPS IC') + ! *DH 20200922 ! call get_data_source(source,Atm%flagstruct%regional) From 1a24fa4254124c27a89a4257af2bbbd9ee3734de Mon Sep 17 00:00:00 2001 From: Xiaqiong Zhou Date: Fri, 16 Oct 2020 16:44:28 +0000 Subject: [PATCH 5/6] To avoid the error ( erroneous arithmetic operation) with the GNU compiler when DEBUG=Y --- model/tp_core.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/model/tp_core.F90 b/model/tp_core.F90 index 6f8e810de..0a19d524f 100644 --- a/model/tp_core.F90 +++ b/model/tp_core.F90 @@ -495,7 +495,8 @@ subroutine xppm(flux, q, c, iord, is,ie,isd,ied, jfirst,jlast,jsd,jed, npx, npy, 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( q1(i)+0.25/a4(i)*da1(i)**2+a4(i)*r12 < 0. ) then + if( q1(i)+0.25*da1(i)**2/a4(i)+a4(i)*r12 < 0. ) then if( .not. smt5(i) ) then br(i) = 0. bl(i) = 0. @@ -823,7 +824,8 @@ subroutine yppm(flux, q, c, jord, ifirst,ilast, isd,ied, js,je,jsd,jed, npx, npy 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( q(i,j)+0.25/a4(i)*xt1(i)**2+a4(i)*r12 < 0. ) then + if( q(i,j)+0.25*xt1(i)**2/a4(i)+a4(i)*r12 < 0. ) then if( .not. smt5(i,j) ) then br(i,j) = 0. bl(i,j) = 0. From 61875852b52951f6c6215603a19c826b952fc534 Mon Sep 17 00:00:00 2001 From: Jun Wang <37633869+junwang-noaa@users.noreply.github.com> Date: Tue, 17 Nov 2020 15:16:55 -0500 Subject: [PATCH 6/6] add GFSv16 dzmin change and changes for Dev/jcsda (#44) * add GFSv16 dzmin change * Add code changes in external_ic.F90 and fv_grid_tools.F90 for dev/jcsda, dycore PR #35 Co-authored-by: Jun Wang Co-authored-by: Dan Holdaway --- model/nh_utils.F90 | 22 +++++++++++----------- tools/external_ic.F90 | 13 +++++-------- tools/fv_grid_tools.F90 | 8 ++++---- 3 files changed, 20 insertions(+), 23 deletions(-) diff --git a/model/nh_utils.F90 b/model/nh_utils.F90 index c6ada26c0..893a8725c 100644 --- a/model/nh_utils.F90 +++ b/model/nh_utils.F90 @@ -66,7 +66,7 @@ module nh_utils_mod public sim3p0_solver, rim_2d public Riem_Solver_c - real, parameter:: dz_min = 2. + real, parameter:: dz_min = 6. real, parameter:: r3 = 1./3. CONTAINS @@ -198,14 +198,14 @@ subroutine update_dz_c(is, ie, js, je, km, ng, dt, dp0, zs, area, ut, vt, gz, ws ! Enforce monotonicity of height to prevent blowup !$OMP parallel do default(none) shared(is1,ie1,js1,je1,ws,zs,gz,rdt,km) do j=js1, je1 - do i=is1, ie1 - ws(i,j) = ( zs(i,j) - gz(i,j,km+1) ) * rdt - enddo - do k=km, 1, -1 + do k=2, km+1 do i=is1, ie1 - gz(i,j,k) = max( gz(i,j,k), gz(i,j,k+1) + dz_min ) + gz(i,j,k) = min( gz(i,j,k), gz(i,j,k-1) - dz_min ) enddo enddo + do i=is1, ie1 + ws(i,j) = ( zs(i,j) - gz(i,j,km+1) ) * rdt + enddo enddo end subroutine update_dz_c @@ -312,15 +312,15 @@ subroutine update_dz_d(ndif, damp, hord, is, ie, js, je, km, ng, npx, npy, area, !$OMP parallel do default(none) shared(is,ie,js,je,km,ws,zs,zh,rdt) do j=js, je - do i=is,ie - ws(i,j) = ( zs(i,j) - zh(i,j,km+1) ) * rdt - enddo - do k=km, 1, -1 + do k=2, km+1 do i=is, ie ! Enforce monotonicity of height to prevent blowup - zh(i,j,k) = max( zh(i,j,k), zh(i,j,k+1) + dz_min ) + zh(i,j,k) = min( zh(i,j,k), zh(i,j,k-1) - dz_min ) enddo enddo + do i=is,ie + ws(i,j) = ( zs(i,j) - zh(i,j,km+1) ) * rdt + enddo enddo end subroutine update_dz_d diff --git a/tools/external_ic.F90 b/tools/external_ic.F90 index 44f76c22e..69ef847bc 100644 --- a/tools/external_ic.F90 +++ b/tools/external_ic.F90 @@ -148,7 +148,7 @@ module external_ic_mod use external_sst_mod, only: i_sst, j_sst, sst_ncep use fms_mod, only: file_exist, read_data, field_exist, write_version_number use fms_mod, only: open_namelist_file, check_nml_error, close_file - use fms_mod, only: get_mosaic_tile_file, read_data, error_mesg + use fms_mod, only: get_mosaic_tile_file, error_mesg use fms_io_mod, only: get_tile_string, field_size, free_restart_type use fms_io_mod, only: restart_file_type, register_restart_field use fms_io_mod, only: save_restart, restore_state, set_filename_appendix, get_global_att_value @@ -197,10 +197,11 @@ module external_ic_mod real, parameter:: zvir = rvgas/rdgas - 1. real(kind=R_GRID), parameter :: cnst_0p20=0.20d0 - real :: deg2rad - character (len = 80) :: source ! This tells what the input source was for the data + real, parameter :: deg2rad = pi/180. + character (len = 80),public :: source ! This tells what the input source was for the data character(len=27), parameter :: source_fv3gfs = 'FV3GFS GAUSSIAN NEMSIO FILE' - public get_external_ic, get_cubed_sphere_terrain + public get_external_ic, get_cubed_sphere_terrain + public remap_scalar, remap_dwinds ! version number of this module ! Include variable "version" to be written to log file. @@ -999,8 +1000,6 @@ subroutine get_ncep_ic( Atm, fv_domain, nq ) jsd = Atm%bd%jsd jed = Atm%bd%jed - deg2rad = pi/180. - npz = Atm%npz call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers, num_prog=ntprog) if(is_master()) write(*,*) 'ntracers = ', ntracers, 'ntprog = ',ntprog @@ -1564,8 +1563,6 @@ subroutine get_ecmwf_ic( Atm, fv_domain ) jsd = Atm%bd%jsd jed = Atm%bd%jed - deg2rad = pi/180. - npz = Atm%npz call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers, num_prog=ntprog) if(is_master()) write(*,*) 'ntracers = ', ntracers, 'ntprog = ',ntprog diff --git a/tools/fv_grid_tools.F90 b/tools/fv_grid_tools.F90 index 5a3a964a4..45914aa17 100644 --- a/tools/fv_grid_tools.F90 +++ b/tools/fv_grid_tools.F90 @@ -629,7 +629,7 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, if (Atm%neststruct%nested .or. ANY(Atm%neststruct%child_grids)) then grid_global => Atm%grid_global - else if( trim(grid_file) .NE. 'INPUT/grid_spec.nc') then + else if( trim(grid_file) .EQ. 'Inline') then allocate(grid_global(1-ng:npx +ng,1-ng:npy +ng,ndims,1:nregions)) endif @@ -683,7 +683,7 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, ! still need to set up call setup_aligned_nest(Atm) else - if(trim(grid_file) == 'INPUT/grid_spec.nc') then + if(trim(grid_file) .NE. 'Inline') then call read_grid(Atm, grid_file, ndims, nregions, ng) else if (Atm%flagstruct%grid_type>=0) call gnomonic_grids(Atm%flagstruct%grid_type, npx-1, xs, ys) @@ -1180,8 +1180,8 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, enddo if (Atm%neststruct%nested .or. ANY(Atm%neststruct%child_grids)) then - nullify(grid_global) - else if( trim(grid_file) .NE. 'INPUT/grid_spec.nc') then + nullify(grid_global) + else if( trim(grid_file) .EQ. 'Inline') then deallocate(grid_global) endif