! source file: /Users/jfidler/work/UVic_ESCM/2.9/source/mom/poisson.F subroutine border (v, sym) !----------------------------------------------------------------------- ! adjust borders of an array for cyclic and symmetry settings ! symmetry conditions are: ! 't odd': "t" grid variable: asymmetric reflection at north ! 't even': "t" grid variable: symmetric reflection at north ! 'u odd': "u" grid variable: asymmetric reflection at north ! 'u even': "u" grid variable: symmetric reflection at north !----------------------------------------------------------------------- implicit none character(*) :: sym integer i, j include "size.h" include "stdunits.h" real v(imt,jmt) ! set southern border do i=2,imt-1 v(i,1) = 0.0 enddo do i=2,imt-1 v(i,jmt) = 0.0 enddo do j=1,jmt v(1,j) = v(imt-1,j) v(imt,j) = v(2,j) enddo return end subroutine iborder (iv, sym) !----------------------------------------------------------------------- ! adjust borders of an array for cyclic and symmetry settings ! symmetry conditions are: ! 't odd': "t" grid variable: asymmetric reflection at north ! 't even': "t" grid variable: symmetric reflection at north ! 'u odd': "u" grid variable: asymmetric reflection at north ! 'u even': "u" grid variable: symmetric reflection at north !----------------------------------------------------------------------- implicit none character(*) :: sym integer i, j include "size.h" include "stdunits.h" real iv(imt,jmt) ! set southern border do i=2,imt-1 iv(i,1) = 0 enddo do i=2,imt-1 iv(i,jmt) = 0 enddo do j=1,jmt iv(1,j) = iv(imt-1,j) iv(imt,j) = iv(2,j) enddo return end subroutine checkerboard (solution, map) !----------------------------------------------------------------------- ! removes "checkboard" null space from an array "solution" !----------------------------------------------------------------------- implicit none integer noceansum(0:1,0:1), i1, j1, j, i real redsum, blacksum, nred, nblack, diff, c real sum(0:1,0:1), correction(0:1,0:1) include "size.h" integer map(imt,jmt) real solution(imt,jmt) do i1=0,1 do j1=0,1 sum(i1,j1) = 0.0 noceansum(i1,j1) = 0 enddo enddo do i1=0,1 do j1=0,1 do j=2+j1,jmt-1,2 do i=2+i1,imt-1,2 sum(i1,j1) = sum(i1,j1) + solution(i,j) enddo enddo enddo enddo do i1=0,1 do j1=0,1 do j=2+j1,jmt-1,2 do i=2+i1,imt-1,2 if (map(i,j) .le. 0) then noceansum(i1,j1) = noceansum(i1,j1) + 1 endif enddo enddo enddo enddo redsum = sum(0,0) + sum(1,1) blacksum = sum(1,0) + sum(0,1) nred = noceansum(0,0) + noceansum(1,1) nblack = noceansum(1,0) + noceansum(0,1) diff = redsum/nred - blacksum/nblack c = diff / 2.0 print *, ' ' print '(a,i6,a,i6,a,e14.7)' &, '=> checkerboard: nred = ',nred, ', nblack = ',nblack &, ', removing a checkerboard correction of ', c correction (0,0) = -c correction (1,1) = -c correction (1,0) = c correction (0,1) = c do i1=0,1 do j1=0,1 do j=2+j1,jmt-1,2 do i=2+i1,imt-1,2 if (map(i,j) .le. 0) then solution(i,j) = solution(i,j) + correction(i1,j1) endif enddo enddo enddo enddo return end subroutine fill_land (solution, map, noslip, & nisle, iperm, jperm, iofs, nippts) implicit none integer nisle logical noslip include "size.h" integer map(imt,jmt), iperm(maxipp), jperm(maxipp) integer nippts(mnisle), iofs(mnisle) real solution(imt,jmt) if (noslip) then call fill_land1 (solution, map, & nisle, iperm, jperm, iofs, nippts) else call fill_land2 (solution, map, & nisle, iperm, jperm, iofs, nippts) endif return end subroutine fill_land1 (solution, map, & nisle, iperm, jperm, iofs, nippts) !======================================================================= ! fills each land area with the [presumed constant] value ! of solution along its ocean perimeter. !======================================================================= implicit none integer i, j, isle, nisle real fill include "size.h" integer map(imt,jmt), iperm(maxipp), jperm(maxipp) integer nippts(mnisle), iofs(mnisle) real solution(imt,jmt) do i=2,imt-1 do j=2,jmt-1 if (map(i,j) .gt. 0) then isle = map(i,j) fill = solution(iperm(iofs(isle)+1),jperm(iofs(isle)+1)) solution(i,j) = fill endif enddo enddo call mirror_adjust (solution) return end subroutine fill_land2 (solution, map, & nisle, iperm, jperm, iofs, nippts) !======================================================================= ! fills the boundary cells of each land area with the value ! of solution at the adjacent ocean perimeter point. ! only NSEW directions are searched [no diagonal directions] ! in case of multiple ocean perimeter points, their average ! is used. !======================================================================= implicit none integer i, j, isle, nbrs, nisle logical last_pass real sum include "size.h" integer map(imt,jmt), iperm(maxipp), jperm(maxipp), nippts(mnisle) integer iofs(mnisle) real solution(imt,jmt) 1000 continue last_pass = .true. do i=2,imt-1 do j=2,jmt-1 if (map(i,j) .gt. 0) then isle = map(i,j) nbrs = 0 sum = 0.0 if (map(i,j+1) .eq. -isle) then nbrs = nbrs + 1 sum = sum + solution(i,j+1) endif if (map(i+1,j) .eq. -isle) then nbrs = nbrs + 1 sum = sum + solution(i+1,j) endif if (map(i,j-1) .eq. -isle) then nbrs = nbrs + 1 sum = sum + solution(i,j-1) endif if (map(i-1,j) .eq. -isle) then nbrs = nbrs + 1 sum = sum + solution(i-1,j) endif if (nbrs .gt. 0) then solution(i,j) = sum / nbrs ! last_pass = .false. endif endif enddo enddo call mirror_adjust (solution) if (.not. last_pass) goto 1000 return end subroutine mirror_adjust (solution) !======================================================================= ! fills each border cell with the value adjacent to it !======================================================================= implicit none integer i include "size.h" real solution(imt,jmt) call border(solution, 't odd') do i=1,imt solution(i,1) = solution(i,2) solution(i,jmt) = solution(i,jmt-1) enddo return end subroutine zero_level (surfpres, variable, map, dxt, dyt, cst) implicit none character(*) :: variable integer i, j real sum, area_ocean, area, surfpres0 include "size.h" integer map(imt,jmt) real surfpres(imt,jmt), dxt(imt), dyt(jmt), cst(jmt) ! this does not correctly handle multiple basins sum = 0.0 area_ocean = 0.0 do i=2,imt-1 do j=2,jmt-1 if (map(i,j) .le. 0) then area = dxt(i)*cst(j)*dyt(j) sum = sum + surfpres(i,j)*area area_ocean = area_ocean + area endif enddo enddo surfpres0 = sum / area_ocean call con_adjust (surfpres, surfpres0, map) print '(a,e14.7,a,a/)' &, '=> zero_level: removing a mean of ', surfpres0, ' from ' &, variable return end subroutine ddxu (tquant, uquant, dxu, cosu) !======================================================================= ! Calculates x partial derivative of field tquant ! Answer is centered at u/v points !======================================================================= implicit none integer i, j include "size.h" real tquant(imt,jmt), uquant(imt,jmt) real dxu(imt), cosu(jmt) ! calculate partial derivative = ddx (tquant) call diffdxu (tquant, uquant) do i=1,imt-1 do j=1,jmt-1 uquant(i,j) = uquant(i,j) / (dxu(i)*cosu(j)) enddo enddo return end subroutine diffdxu (tquant, uquant) !======================================================================= ! Calculates x partial difference of field tquant ! Answer is centered at u/v points !======================================================================= implicit none integer i, j, i1, j1 real c0, p5, cddxu(0:1,0:1), cddyu(0:1,0:1) real cddxt(-1:0,-1:0), cddyt(-1:0,-1:0) include "size.h" real tquant(imt,jmt), uquant(imt,jmt) !----------------------------------------------------------------------- ! set locally needed constants !----------------------------------------------------------------------- c0 = 0.0 p5 = 0.5 !----------------------------------------------------------------------- ! construct coefficients for partial differences. a partial ! difference in "x" is defined as an "x" difference of a quantity ! which is averaged in "y". (and symmetrically for "y" differences). ! Note that this is an x difference and NOT an x derivitive. ! partial differences of quantities on the "t" grid are defined on ! the "u" grid and visa versa. ! therefore partial differences at: ! u/v points (i,j), involve nearby t/s points with subscripts: ! (i ,j+1) (i+1,j+1) ! (i ,j ) (i+1,j ) ! t/s points (i,j), involve nearby u/v points with subscripts: ! (i-1,j ) (i ,j ) ! (i-1,j-1) (i ,j-1) ! thus if qu(i,j) is defined on u/v points, its partial ! difference ddxqt = ddxt(qu) is defined on t/s points and has the ! value ! ddxqt(i,j) = cddxt(-1,-1)*qu(i-1,j-1) + cddxt(-1,0)*qu(i-1,j+0) ! + cddxt( 0,-1)*qu(i+0,j-1) + cddxt( 0,0)*qu(i+0,j+0) !----------------------------------------------------------------------- cddxu( 0, 0) = -p5 cddxu( 0, 1) = -p5 cddxu( 1, 0) = p5 cddxu( 1, 1) = p5 cddxt(-1,-1) = -p5 cddxt(-1, 0) = -p5 cddxt( 0,-1) = p5 cddxt( 0, 0) = p5 cddyu( 0, 0) = -p5 cddyu( 0, 1) = p5 cddyu( 1, 0) = -p5 cddyu( 1, 1) = p5 cddyt(-1,-1) = -p5 cddyt(-1, 0) = p5 cddyt( 0,-1) = -p5 cddyt( 0, 0) = p5 !----------------------------------------------------------------------- ! calculate partial difference ! diffdx (tquant) = deltax (bary (tquant)) !----------------------------------------------------------------------- do i=1,imt do j=1,jmt uquant(i,j) = 0.0 enddo enddo do i1=0,1 do j1=0,1 do i=1,imt-1 do j=1,jmt-1 uquant(i,j) = uquant(i,j) + cddxu(i1,j1)*tquant(i+i1,j+j1) enddo enddo enddo enddo return end subroutine ddyu (tquant, uquant, dyu) !======================================================================= ! Calculates y partial derivative of field tquant ! Answer is centered at u/v points !======================================================================= implicit none integer i, j include "size.h" real tquant(imt,jmt), uquant(imt,jmt), dyu(jmt) call diffdyu (tquant, uquant) do i=1,imt-1 do j=1,jmt-1 uquant(i,j) = uquant(i,j) / dyu(j) enddo enddo return end subroutine diffdyu (tquant, uquant) !======================================================================= ! Calculates y partial difference of field tquant ! Answer is centered at u/v points !======================================================================= implicit none integer i, j, i1, j1 real c0, p5, cddxu(0:1,0:1), cddyu(0:1,0:1) real cddxt(-1:0,-1:0), cddyt(-1:0,-1:0) include "size.h" real tquant(imt,jmt), uquant(imt,jmt) !----------------------------------------------------------------------- ! set locally needed constants !----------------------------------------------------------------------- c0 = 0.0 p5 = 0.5 !----------------------------------------------------------------------- ! construct coefficients for partial differences. a partial ! difference in "x" is defined as an "x" difference of a quantity ! which is averaged in "y". (and symmetrically for "y" differences). ! Note that this is an x difference and NOT an x derivitive. ! partial differences of quantities on the "t" grid are defined on ! the "u" grid and visa versa. ! therefore partial differences at: ! u/v points (i,j), involve nearby t/s points with subscripts: ! (i ,j+1) (i+1,j+1) ! (i ,j ) (i+1,j ) ! t/s points (i,j), involve nearby u/v points with subscripts: ! (i-1,j ) (i ,j ) ! (i-1,j-1) (i ,j-1) ! thus if qu(i,j) is defined on u/v points, its partial ! difference ddxqt = ddxt(qu) is defined on t/s points and has the ! value ! ddxqt(i,j) = cddxt(-1,-1)*qu(i-1,j-1) + cddxt(-1,0)*qu(i-1,j+0) ! + cddxt( 0,-1)*qu(i+0,j-1) + cddxt( 0,0)*qu(i+0,j+0) !----------------------------------------------------------------------- cddxu( 0, 0) = -p5 cddxu( 0, 1) = -p5 cddxu( 1, 0) = p5 cddxu( 1, 1) = p5 cddxt(-1,-1) = -p5 cddxt(-1, 0) = -p5 cddxt( 0,-1) = p5 cddxt( 0, 0) = p5 cddyu( 0, 0) = -p5 cddyu( 0, 1) = p5 cddyu( 1, 0) = -p5 cddyu( 1, 1) = p5 cddyt(-1,-1) = -p5 cddyt(-1, 0) = p5 cddyt( 0,-1) = -p5 cddyt( 0, 0) = p5 !----------------------------------------------------------------------- ! calculate partial difference ! diffdy (tquant) = deltay (barx (tquant)) !----------------------------------------------------------------------- do i=1,imt do j=1,jmt uquant(i,j) = 0.0 enddo enddo do i1=0,1 do j1=0,1 do i=1,imt-1 do j=1,jmt-1 uquant(i,j) = uquant(i,j) + cddyu(i1,j1)*tquant(i+i1,j+j1) enddo enddo enddo enddo return end subroutine con_adjust (dpsi, dpsi1, map) !----------------------------------------------------------------------- ! the constant dpsi1 is subtracted from dpsi(i,j) at all ! ocean points (i.e., where map(i,j) .le. 0) !----------------------------------------------------------------------- implicit none integer i, j real dpsi1 include "size.h" integer map(imt,jmt) real dpsi(imt,jmt) do i=1,imt do j=1,jmt if (map(i,j) .le. 0) then dpsi(i,j) = dpsi(i,j) - dpsi1 endif enddo enddo return end function conv (converged) !----------------------------------------------------------------------- ! converts logical to character form for printing !----------------------------------------------------------------------- implicit none character(11) :: conv logical converged if (converged) then conv = '[converged]' else conv = '[diverged] ' endif return end subroutine compare2(x, y, ax, ay, imax, jmax) !----------------------------------------------------------------------- ! compare arrays "x" and "y" ! ax = alphabetical identifier for "x" ! ay = alphabetical identifier for "y" ! prints count of relative differences greater than threshhold ! which is set in parameter statement below !----------------------------------------------------------------------- implicit none character(*) :: ax, ay integer imax, jmax, numneq, numbadrel, i, j real threshhold, relerrormax, badabserrormax, abserrormax real relerror, relerr, x(imax,jmax), y(imax,jmax) parameter (threshhold= 1.0e-6) write (6,*) 'comparing ',ax,' with ',ay numneq = 0 numbadrel = 0 relerrormax = 0.0 badabserrormax = 0.0 abserrormax = 0.0 do i=1,imax do j=1,jmax if (x(i,j) .ne. y(i,j) ) then numneq = numneq + 1 relerror = relerr (x(i,j), y(i,j)) relerrormax = max (relerror, relerrormax) if (numneq .le. 20 .and. & relerror .gt. threshhold) then write(6,'(3a10,2i4,tr1,2e25.18,a,e25.18)') & '-->',ax,ay,i,j,x(i,j),y(i,j), & ' rel error = ', relerror endif if (relerror .gt. threshhold) then numbadrel = numbadrel + 1 badabserrormax = max(badabserrormax, abs(x(i,j)-y(i,j))) endif abserrormax = max(abserrormax, abs(x(i,j)-y(i,j))) endif enddo enddo if (numneq .ne. 0) then write(6,*) numneq, ' entries differ' write(6,*) numbadrel, ' have relative error >', threshhold write(6,*) ' bad max absolute error = ',badabserrormax write(6,*) ' maximum relative error = ',relerrormax write(6,*) ' maximum absolute error = ',abserrormax else write(6,*) 'no differences detected' endif return end function relerr (x, y) implicit none real x, y, relerr, xymin, xymax if (x .eq. y) then relerr = 0.0 else xymin = min(abs(x), abs(y)) xymax = max(abs(x), abs(y)) if (xymin .gt. 1.0e-23 * xymax) then relerr = abs(x-y)/xymin else relerr = xymax endif endif return end