!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS 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 General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module fill_module !----------------------------------------------------------------------- !BOP ! implicit none ! !MODULE: fill_module --- utilities for filling in nagative values ! ! !PUBLIC MEMBER FUNCTIONS: public filew, fillxy, fillz, pfix ! ! !DESCRIPTION: ! ! !REVISION HISTORY: ! 99.03.01 Lin Creation ! 01.02.14 Lin Routines coalesced into this module ! 01.03.26 Sawyer Added ProTeX documentation ! !EOP !----------------------------------------------------------------------- contains !----------------------------------------------------------------------- !BOP ! !IROUTINE: fillxy --- Fill from east, west, north and south neighbors ! ! !INTERFACE: subroutine fillxy(q, im, jm, jfirst, jlast, acap, cosp, acosp) integer, intent(in):: im ! Longitudes integer, intent(in):: jm ! Total latitudes integer, intent(in):: jfirst ! Starting latitude integer, intent(in):: jlast ! Finishing latitude real , intent(in):: acap real , intent(in):: cosp(jm) real , intent(in):: acosp(jm) real , intent(inout):: q(im,jfirst:jlast) ! Field to adjust ! ! LOCAL VARIABLES: integer ipx, i,j real :: tiny = 1.e-40 call filew(q,im,jm,jfirst,jlast,acap,ipx,tiny, cosp(2), cosp(3)) ! if(ipx.ne.0) then ! Fill zonally ! call xfix(q,IM,JM,tiny,qt) ! endif !EOC end subroutine fillxy !----------------------------------------------------------------------- !BOP ! !IROUTINE: filew --- Fill from east and west neighbors; essentially ! performing local flux adjustment ! ! !INTERFACE: subroutine filew(q, im, jm, jfirst, jlast, acap, ipx, tiny, cosp2, cosp3) ! !INPUT PARAMETERS: integer im ! Longitudes integer jm ! Total latitudes integer jfirst ! Starting latitude integer jlast ! Finishing latitude real tiny ! A small number to pump up value real acap ! 1/(polar cap area) real cosp2 ! cosine(lat) at j=2 real cosp3 ! cosine(lat) at j=3 ! !INPUT/OUTPUT PARAMETERS: real q(im,jfirst:jlast) ! Field to adjust ! !OUTPUT PARAMETERS: integer ipx ! Flag: 0 if Q not change, 1 if changed ! !DESCRIPTION: ! Check for "bad" data and fill from east and west neighbors ! ! !REVISION HISTORY: ! 99.10.01 Lin Creation ! 07.30.01 Lin Improvement ! !EOP !----------------------------------------------------------------------- !BOC ! !LOCAL VARIABLES: real d0, d1, d2 real qtmp(jfirst:jlast,im) integer i, j, jm1, ip2 integer j1, j2 jm1 = jm-1 j1 = max( jfirst, 2 ) j2 = min( jlast, jm1 ) ipx = 0 do j=jfirst,jlast do i=1,im if( q(i,j) < 0. ) then ipx = 1 go to 444 endif enddo enddo if ( ipx == 0 ) return ! Copy & swap direction for vectorization. 444 do i=1,im do j=j1,j2 qtmp(j,i) = q(i,j) enddo enddo do i=2,im-1 do j=j1,j2 if(qtmp(j,i) < 0.) then ! west d0 = max(0.,qtmp(j,i-1)) d1 = min(-qtmp(j,i),d0) qtmp(j,i-1) = qtmp(j,i-1) - d1 qtmp(j,i) = qtmp(j,i) + d1 ! east d0 = max(0.,qtmp(j,i+1)) d2 = min(-qtmp(j,i),d0) qtmp(j,i+1) = qtmp(j,i+1) - d2 qtmp(j,i) = qtmp(j,i) + d2 + tiny endif enddo enddo i=1 do j=j1,j2 if(qtmp(j,i) < 0.) then ! west d0 = max(0.,qtmp(j,im)) d1 = min(-qtmp(j,i),d0) qtmp(j,im) = qtmp(j,im) - d1 qtmp(j,i) = qtmp(j,i) + d1 ! east d0 = max(0.,qtmp(j,i+1)) d2 = min(-qtmp(j,i),d0) qtmp(j,i+1) = qtmp(j,i+1) - d2 qtmp(j,i) = qtmp(j,i) + d2 + tiny endif enddo i=im do j=j1,j2 if(qtmp(j,i) < 0.) then ! west d0 = max(0.,qtmp(j,i-1)) d1 = min(-qtmp(j,i),d0) qtmp(j,i-1) = qtmp(j,i-1) - d1 qtmp(j,i) = qtmp(j,i) + d1 ! east d0 = max(0.,qtmp(j,1)) d2 = min(-qtmp(j,i),d0) qtmp(j,1) = qtmp(j,1) - d2 qtmp(j,i) = qtmp(j,i) + d2 + tiny endif enddo !----------- ! Final pass !----------- do i=1,im-1 do j=j1,j2 if (qtmp(j,i) < 0. ) then ! Take mass from east (essentially adjusting fx(i+1,j)) qtmp(j,i+1) = qtmp(j,i+1) + qtmp(j,i) qtmp(j,i) = 0. endif enddo enddo do i=im,2,-1 do j=j1,j2 if (qtmp(j,i) < 0. ) then ! Take mass from west (essentially adjusting fx(i,j)) qtmp(j,i-1) = qtmp(j,i-1) + qtmp(j,i) qtmp(j,i) = 0. endif enddo enddo do j=j1,j2 do i=1,im q(i,j) = qtmp(j,i) enddo enddo ! Check Poles. if ( jfirst == 1 ) then ip2 = 0 if(q(1,1) < 0.) then call pfix(q(1,2),q(1,1),im,ipx,acap,cosp2) else ! Check the latitude next to the SP do i=1,im if(q(i,2) < 0.) then ip2 = 1 go to 322 endif enddo 322 continue if(ip2.ne.0) call pfix(q(1,2),q(1,1),im,ipx,acap,cosp2) endif endif if ( jlast == jm ) then if(q(1,jm) < 0.) then call pfix(q(1,jm1),q(1,jm),im,ipx,acap,cosp2) else ! Check the latitude next to the NP ip2 = 0 do i=1,im if(q(i,jm1) < 0.) then ip2 = 1 go to 323 endif enddo 323 continue if(ip2.ne.0) call pfix(q(1,jm1),q(1,jm),im,ipx,acap,cosp2) endif endif !EOC end subroutine filew !----------------------------------------------------------------------- !----------------------------------------------------------------------- !BOP ! !IROUTINE: fillz --- Fill from neighbors below and above ! ! !INTERFACE: subroutine fillz(im, km, nq, q, dp) ! !INPUT PARAMETERS: integer, intent(in):: im ! No. of longitudes integer, intent(in):: km ! No. of levels integer, intent(in):: nq ! Total number of tracers real , intent(in):: dp(im,km) ! pressure thickness ! !INPUT/OUTPUT PARAMETERS: real , intent(inout) :: q(im,km,nq) ! tracer mixing ratio ! !DESCRIPTION: ! Check for "bad" data and fill from east and west neighbors ! ! !BUGS: ! Currently this routine only performs the east-west fill algorithm. ! This is because the N-S fill is very hard to do in a reproducible ! fashion when the problem is decomposed by latitudes. ! ! !REVISION HISTORY: ! 00.04.01 Lin Creation ! !EOP !----------------------------------------------------------------------- !BOC ! ! !LOCAL VARIABLES: integer i, k, ic real qup, qly, dup do ic=1,nq ! Top layer do i=1,im if( q(i,1,ic) < 0.) then q(i,2,ic) = q(i,2,ic) + q(i,1,ic)*dp(i,1)/dp(i,2) q(i,1,ic) = 0. endif enddo ! Interior do k=2,km-1 do i=1,im if( q(i,k,ic) < 0. ) then ! Borrow from above qup = q(i,k-1,ic)*dp(i,k-1) qly = -q(i,k ,ic)*dp(i,k ) dup = min( 0.75*qly, qup ) !borrow no more than 75% from top q(i,k-1,ic) = q(i,k-1,ic) - dup/dp(i,k-1) ! Borrow from below: q(i,k,ic) is still negative at this stage q(i,k+1,ic) = q(i,k+1,ic) + (dup-qly)/dp(i,k+1) q(i,k ,ic) = 0. endif enddo enddo ! Bottom layer k = km do i=1,im if( q(i,k,ic) < 0.) then ! Borrow from above qup = q(i,k-1,ic)*dp(i,k-1) qly = -q(i,k ,ic)*dp(i,k ) dup = min( qly, qup ) q(i,k-1,ic) = q(i,k-1,ic) - dup/dp(i,k-1) q(i,k,ic) = 0. endif enddo enddo end subroutine fillz !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! subroutine pfix(q, qp, im, ipx, acap, cosp2) integer im ! Longitudes real acap ! ??? real cosp2 ! ??? real q(im) ! Latitude-level field to adjust real qp(im) ! Second latitude-level field to adjust (usually pole) ! !OUTPUT PARAMETERS: integer ipx ! Flag: 0 if Q not change, 1 if changed ! !LOCAL VARIABLES: integer i real summ, sump, pmean summ = 0. sump = 0. do i=1,im summ = summ + q(i) sump = sump + qp(i) enddo sump = sump/im pmean = (sump*acap + summ*cosp2) / (acap + cosp2*im) do i=1,im q(i) = pmean qp(i) = pmean enddo if( qp(1) < 0. ) then ipx = 1 endif !EOC end subroutine pfix !----------------------------------------------------------------------- end module fill_module