module fill_module!-----------------------------------------------------------------------!BOP!! !MODULE: fill_module --- utilities for filling in "bad" data!! !PUBLIC MEMBER FUNCTIONS:      public filew, fillxy, fillz, filns, pfix!! !DESCRIPTION:!!    This module provides the basic utilities to fill in regions!    with bad "data", for example slightly negative values in fields!    which must be positive, like mixing ratios.  Generally this !    means borrowing positive values from neighboring cells.!! !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: filew --- Fill from east and west neighbors; essentially!                      performing local flux adjustment!! !INTERFACE:  subroutine filew(q, im, jm, jfirst, jlast, acap, ipx, tiny, cosp2)! !USES: use precision implicit none! !INPUT PARAMETERS: integer im                  ! Longitudes integer jm                  ! Total latitudes integer jfirst              ! Starting latitude integer jlast               ! Finishing latitude real(r8) tiny               ! A small number to pump up value real(r8) acap               ! 1/(polar cap area) real(r8) cosp2              ! cosine(lat) at j=2! !INPUT/OUTPUT PARAMETERS: real(r8) 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(r8) d0, d1, d2 real(r8) qtmp(jfirst:jlast,im) integer i, j, jm1, ip2 integer j1, j2 j1 = max( jfirst,   2 ) j2 = min( jlast, jm-1 ) jm1 = jm-1 ipx = 0! Copy & swap direction for vectorization.  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           ipx =  1! 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        ipx =  1! 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        ipx =  1! 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 if(ipx .ne. 0) then!-----------! 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 endif ! Check Poles. if ( jfirst == 1 ) then      if(q(1,1) < 0.) then         call pfix(q(1,2),q(1,1),im,ipx,acap,cosp2)      else!            Check j=2             ip2 = 0         do i=1,im            if(q(i,2).lt.0.) then               ip2 = 1               go to 322            endif         enddo322      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 j=jm1              ip2 = 0         do i=1,im            if(q(i,jm1) < 0.) then               ip2 = 1               go to 323            endif         enddo323      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: fillxy --- Fill from east, west, north and south neighbors!! !INTERFACE:  subroutine fillxy(q, im, jm, jfirst, jlast, acap, cosp, acosp)! !USES: use precision implicit none integer im                  ! Longitudes integer jm                  ! Total latitudes integer jfirst              ! Starting latitude integer jlast               ! Finishing latitude real(r8) acap               ! ??? real(r8) cosp(jm)           ! ??? real(r8) acosp(jm)          ! ???!! !INPUT/OUTPUT PARAMETERS: real(r8) q(im,jfirst:jlast) ! Field to adjust! !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:!   99.03.01   Lin        Creation!!EOP!-----------------------------------------------------------------------!BOC!! !LOCAL VARIABLES:  integer ipx, ipy, j1, j2  real(r8) tiny  parameter( tiny = 1.e-20 )    call filew(q,im,jm,jfirst,jlast,acap,ipx,tiny,cosp(2))! WS 99.08.03 : S.-J. can you clean up the j1, j2 stuff here?   if(ipx.ne.0) then      j1 = max( 2,    jfirst )      j2 = min( jm-1, jlast )!! WS 99.08.03 : see comments in "BUGS" above!!!      call filns(q,im,jm,j1,j2,cosp,acosp,ipy,tiny)!     if(ipy .ne. 0) then! do fill zonally! xfx is problematic!     call xfix(q,IM,JM,tiny,qt)!     endif   endif!EOC end subroutine fillxy!-----------------------------------------------------------------------!-----------------------------------------------------------------------!BOP! !IROUTINE: fillz --- Fill from neighbors below and above!! !INTERFACE:  subroutine fillz(im, i1, i2, km, nq, q, dp)! !USES: use precision implicit none! !INPUT PARAMETERS:   integer, intent(in) :: im                ! No. of longitudes   integer, intent(in) :: km                ! No. of levels   integer, intent(in) :: i1                ! Starting longitude   integer, intent(in) :: i2                ! Finishing longitude   integer, intent(in) :: nq                ! Total number of tracers   real(r8), intent(in) ::  dp(im,km)       ! pressure thickness! !INPUT/OUTPUT PARAMETERS:   real(r8), 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(r8) qup, qly, dup   do ic=1,nq! Top layer      do i=i1,i2         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=i1,i2         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.5*qly, qup )        !borrow no more than 50%             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=i1,i2         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!EOCend subroutine fillz!-----------------------------------------------------------------------!-----------------------------------------------------------------------!BOP! !IROUTINE: filns --- Fill from north and south neighbors!! !INTERFACE:  subroutine filns(q,im,jm,j1,j2,cosp,acosp,ipy,tiny)! !USES: use precision implicit none! !INPUT PARAMETERS: integer im                  ! Longitudes integer jm                  ! Total latitudes integer j1                  ! Starting latitude integer j2                  ! Finishing latitude real(r8) tiny               ! A small number to pump up value real(r8) cosp(*)            ! ??? real(r8) acosp(*)           ! ???! !INPUT/OUTPUT PARAMETERS: real(r8) q(im,*)            ! Field to adjust! !OUTPUT PARAMETERS: integer  ipy                ! Flag: 0 if no fill-in, 1 if fill-in! !DESCRIPTION:!   Check for "bad" data and fill from north and south neighbors!! !BUGS:!   Currently this routine can only be used performs when the!   problem is *not* distributed in latitude (i.e. j1=1, j2=jm).!   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:!   99.03.01   Lin        Creation!!EOP!-----------------------------------------------------------------------!BOC!! !LOCAL VARIABLES: integer  i, j real(r8) dp, cap1, dq, dn, ds, d0, d1, d2 logical first data first /.true./ save cap1    if(first) then      dp = 4.d0*datan(1.d0)/dble(jm-1)      cap1 = im*(1.-cos((j1-1.5)*dp))/dp      first = .false.    endif       ipy = 0    do j=j1+1,j2-1      do i=1,im      if(q(i,j).lt.0.) then         ipy =  1         dq  = - q(i,j)*cosp(j)! North         dn = q(i,j+1)*cosp(j+1)         d0 = max(0.,dn)         d1 = min(dq,d0)         q(i,j+1) = (dn - d1)*acosp(j+1)         dq = dq - d1! South         ds = q(i,j-1)*cosp(j-1)         d0 = max(0.,ds)         d2 = min(dq,d0)         q(i,j-1) = (ds - d2)*acosp(j-1)         q(i,j) = (d2 - dq)*acosp(j) + tiny      endif      enddo    enddo       do i=1,im      if(q(i,j1).lt.0.) then      ipy =  1      dq  = - q(i,j1)*cosp(j1)! North      dn = q(i,j1+1)*cosp(j1+1)      d0 = max(0.,dn)      d1 = min(dq,d0)      q(i,j1+1) = (dn - d1)*acosp(j1+1)      q(i,j1) = (d1 - dq)*acosp(j1) + tiny      endif      enddo       j = j2      do i=1,im      if(q(i,j).lt.0.) then      ipy =  1      dq  = - q(i,j)*cosp(j)! South      ds = q(i,j-1)*cosp(j-1)      d0 = max(0.,ds)      d2 = min(dq,d0)      q(i,j-1) = (ds - d2)*acosp(j-1)      q(i,j) = (d2 - dq)*acosp(j) + tiny      endif      enddo ! Check Poles.      if(q(1,1).lt.0.) then      dq = q(1,1)*cap1/float(im)*acosp(j1)      do i=1,im      q(i,1) = tiny      q(i,j1) = q(i,j1) + dq      q(i,j1) = max(tiny, q(i,j1) + dq )      enddo      endif       if(q(1,jm).lt.0.) then      dq = q(1,jm)*cap1/float(im)*acosp(j2)      do i=1,im      q(i,jm) = tiny      q(i,j2) = max(tiny,  q(i,j2) + dq )      enddo      endif!EOC  end subroutine filns!-----------------------------------------------------------------------!-----------------------------------------------------------------------!BOP! !IROUTINE: pfix --- fix an individual latitude-level!! !INTERFACE:  subroutine pfix(q, qp, im, ipx, acap, cosp2)! !USES: use precision implicit none! !INPUT PARAMETERS: integer im                  ! Longitudes real(r8) acap               ! ??? real(r8) cosp2              ! ???! !INPUT/OUTPUT PARAMETERS: real(r8) q(im)              ! Latitude-level field to adjust real(r8) qp(im)             ! Second latitude-level field to adjust (usually pole)! !OUTPUT PARAMETERS: integer ipx                 ! Flag:  0 if Q not change, 1 if changed! !DESCRIPTION:!   Fill one latitude-level from east and west neighbors!! !REVISION HISTORY:!   99.03.01   Lin        Creation!!EOP!-----------------------------------------------------------------------!BOC!! !LOCAL VARIABLES: integer i real(r8) 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