program setregionmask !------------------------------------------------------------------- ! PURPOSE: ! o Reads POP input binaries kmt.da & grid.pop.da for a given grid, ! and writes out a region.ieeei4 file. ! o For modern topography, generates full region set; for paleo ! topography, generates north/south mask. ! o USER must specify boundary segment start/end points which, ! together with land, will completely box in a region. A ! boundary set for one region also serves as a boundary ! for the region immediately adjacent. ! AND ! must provide a seed point (i_region,j_region) from which ! a given region will be filled. ! !------------------------------------------------------------------- ! ! Compile on tempest with f90 setregionmask.f90 ! ! 8/07/2006: settings appropriate for gx1v5 region masking ! !------------------------------------------------------------------- implicit none !--- local --- integer :: imt, jmt, i, j, k, recl_factor,nb,nbad integer :: imiss=99 logical :: paleo integer*4, allocatable :: kmt(:,:), mask(:,:), cnt(:,:) real*8, allocatable :: ulat(:,:),ulng(:,:),tlat(:,:),tlng(:,:) real :: atlantic(6,4),pacific(5,4),indian(2,4) real :: north(4),northsea(2,4),labrador(2,4) character(len=80) :: gridf, kmtf, maskf ! line segments which demarc region boundaries ! /x0, y0, x1, y1/ in degrees ! data north(1:4) /0., 0., 360., 0. / !N. Hemisphere data atlantic(1,1:4) /300.,-34., 22., -34. / !S. Atlantic data atlantic(2,1:4) /353., 32., 353., 39. / !Med. data atlantic(3,1:4) /8., 60., 341., 65. / !Finland-Iceland data atlantic(4,1:4) /341., 65., 325., 72. / !Iceland-Greenland data atlantic(5,1:4) /315., 62., 315., 50. / !Greenland-Atlantic data atlantic(6,1:4) /300., 51., 315., 50. / !Atlantic-NovaScotia data pacific(1,1:4) /143.,-34., 300., -34. / !S. Pacific data pacific(2,1:4) /143.,-20., 143., -6. / !Torres Strait (?) data pacific(3,1:4) /130.,-1.5, 116., -1.5 / !Borneo Strait (?) data pacific(4,1:4) /102.,-1.5, 112., -1.5 / !Java Strait (?) data pacific(5,1:4) /180.,66., 203., 66. / !Bering Strait data indian(1,1:4) /26.,-34., 117., -34. / !S. Pacific data indian(2,1:4) /56.,22., 56., 28. / !Persian Gulf data northsea(1,1:4) /26.,70., 16.0, 79. / !Norway-Svalbard data northsea(2,1:4) /16.6,79., 336., 79. / !Svalbard-Greenland data labrador(1,1:4) /295.,58.,295.,65./ data labrador(2,1:4) /279.,75.,279.,78./ print *, 'Enter grid dimensions : imt, jmt' read *, imt, jmt print *, 'Enter name of POP grid binary:' read *, gridf print *, 'Enter name of POP kmt binary:' read *, kmtf print *, 'Enter name of POP region mask binary (output):' read *, maskf print *, 'Is this a paleo grid (T/F)?' read *, paleo allocate(kmt(imt,jmt)) allocate(mask(imt,jmt)) allocate(cnt(imt,jmt)) allocate(ulat(imt,jmt)) allocate(ulng(imt,jmt)) allocate(tlat(imt,jmt)) allocate(tlng(imt,jmt)) open (1, file=kmtf, access='direct', & form='unformatted', recl=jmt*imt*4, status='old') open (2, file=gridf, access='direct', & form='unformatted', recl=jmt*imt*8, status='old') open (3, file=maskf, access='direct', & form='unformatted', recl=jmt*imt*4, status='unknown', & action='write') ! read in kmt and grid: read (1, rec=1) kmt read (2, rec=1) ulat read (2, rec=2) ulng close(1) close(2) ! get T grid locations call calc_tpoints(imt,jmt,ulng,ulat,tlng,tlat) ! mask 99 for ocean, 0 for land mask = imiss where(kmt == 0) mask=0 cnt = 0 ! set region boundaries (where necessary), then fill if (paleo) then ! For paleo grids, presumably do something simple like ! Northern Hemisphere & Southern Hemisphere nb = size(north,1) call set_bounds(imt,jmt,tlng,tlat,mask,nb,north,1,'north') call neighbor_fill(imt,jmt,mask,1,3,87,imiss) ! fill North call neighbor_fill(imt,jmt,mask,2,3,40,imiss) ! fill South else nb = size(atlantic,1) print *, 'Setting atlantic bounds...' call set_bounds(imt,jmt,tlng,tlat,mask,nb,atlantic,6,'atlantic') nb = size(pacific,1) print *, 'Setting pacific bounds...' call set_bounds(imt,jmt,tlng,tlat,mask,nb,pacific,2,'pacific') nb = size(indian,1) print *, 'Setting indian bounds...' call set_bounds(imt,jmt,tlng,tlat,mask,nb,indian,3,'indian') nb = size(northsea,1) print *, 'Setting northsea bounds...' call set_bounds(imt,jmt,tlng,tlat,mask,nb,northsea,9,'northsea') nb = size(labrador,1) print *, 'Setting labsea bounds...' call set_bounds(imt,jmt,tlng,tlat,mask,nb,labrador,8,'labrador') ! Filling regions from seed points. Note that some regions require ! more than one fill because they are not contiguous call neighbor_fill(imt,jmt,mask,6,309,272,imiss) ! fill Atlantic call neighbor_fill(imt,jmt,mask,2,199,212,imiss) ! fill Pacific call neighbor_fill(imt,jmt,mask,3,98,149,imiss) ! fill Indian call neighbor_fill(imt,jmt,mask,1,234,39,imiss) ! fill Southern Ocean call neighbor_fill(imt,jmt,mask,8,305,348,imiss) ! fill Labrador call neighbor_fill(imt,jmt,mask,9,57,370,imiss) ! fill North Sea call neighbor_fill(imt,jmt,mask,10,171,362,imiss) ! fill Arctic call neighbor_fill(imt,jmt,mask,7,51,288,imiss) ! fill Med call neighbor_fill(imt,jmt,mask,-13,65,302,imiss) ! fill Black call neighbor_fill(imt,jmt,mask,-14,81,306,imiss) ! fill Caspian call neighbor_fill(imt,jmt,mask,-5,71,247,imiss) ! fill Red call neighbor_fill(imt,jmt,mask,4,83,266,imiss) ! fill Persian call neighbor_fill(imt,jmt,mask,-12,57,334,imiss) ! fill Baltic call neighbor_fill(imt,jmt,mask,11,273,342,imiss) ! fill Hudson endif where(mask==imiss) cnt=1 nbad = sum(cnt) if (nbad /= 0) then write(6,*) 'Still have ',nbad,' points with unassigned value' endif write (unit=3, rec=1)mask write(6,*) ' ' write(6,*) 'DONE writing ',maskf end program setregionmask !================================================================= recursive SUBROUTINE neighbor_fill(imt,jmt,mask,value,i0,j0,miss) !----------------------------------------------------------------------- ! Recursively fills a region starting from seed point (i0,j0), and ! filling with "value" wherever north/south/east/west shifts indicate ! a missing value. !----------------------------------------------------------------------- integer :: i0, j0, imt, jmt, miss, value integer :: ie,iw,jn,js integer, dimension(imt,jmt) :: mask if (mask(i0,j0) == miss) then mask(i0,j0)=value iw = iadd(i0,-1,imt) call neighbor_fill(imt,jmt,mask,value,iw,j0,miss) ie = iadd(i0,1,imt) call neighbor_fill(imt,jmt,mask,value,ie,j0,miss) jn = iadd(j0,1,jmt) call neighbor_fill(imt,jmt,mask,value,i0,jn,miss) js = iadd(j0,-1,jmt) call neighbor_fill(imt,jmt,mask,value,i0,js,miss) endif END SUBROUTINE neighbor_fill !================================================================= SUBROUTINE set_bounds(imt,jmt,tlng,tlat,mask,nb,bounds,value,rstr) !----------------------------------------------------------------------- ! ! Sets MASK array to VALUE along linear segments specified in the ! BOUNDS array (there are NB segments to process for a given region). ! The line segments are meant to follow great circles which intersect ! the endpoints in BOUNDS, but in practice, a slope is computed in ! grid i,j space, and this slope is traced from a mesh starting cell ! to a mesh ending cell. ! !----------------------------------------------------------------------- integer :: i, j, k, imt, jmt, nb, value, kcnt, lcnt integer :: sti,stj,eni,enj,thisi,thisj integer, allocatable :: movei(:),movej(:) integer, dimension(imt,jmt) :: mask, iind, jind real, dimension(nb,4) :: bounds real :: x0,x1,y0,y1 real*8, dimension(imt,jmt) :: tlng,tlat,dist0,dist1 real*8 :: mindist0,mindist1 character(len=*) :: rstr rstr = trim(rstr) do i=1,imt iind(i,1:jmt)=i enddo do j=1,jmt jind(1:imt,j)=j enddo ! loop over number of boundary segments for this region do k=1,nb sti=0 stj=0 eni=0 enj=0 lcnt=0 kcnt=0 x0 = bounds(k,1) y0 = bounds(k,2) x1 = bounds(k,3) y1 = bounds(k,4) dist0 = sqrt((tlng-x0)**2+(tlat-y0)**2) dist1 = sqrt((tlng-x1)**2+(tlat-y1)**2) mindist0 = 1.e30 mindist1 = 1.e30 do j=1,jmt do i=1,imt if (dist0(i,j) imt/2) then if (idx > 0) idx = idx-imt if (idx < 0) idx = idx+imt endif nidx = abs(idx) nidy = abs(idy) ! factor out common denominators in slope idy/idx 10 continue do i=2,min(nidx,nidy) if (mod(idx,i)==0 .and. mod(idy,i)==0) then idx = idx/i idy = idy/i goto 10 endif enddo ! now that slope is reduced, determine best move sequence ! write(6,*) 'using slope DY = ',idy,' DX = ',idx nmove = nidx+nidy allocate(movei(nmove)) allocate(movej(nmove)) call get_moves(nmove,idx,idy,movei,movej) ! write(6,*) 'movei = ',(movei(i),i=1,nmove) ! write(6,*) 'movej = ',(movej(i),i=1,nmove) ! step from starting point on mesh (sti,stj) to adjacent cells until ! you get to the grid cell nearest the segment end (eni,enj) ! Set ocean point mask to "value", leave land point mask values 0 thisi = sti thisj = stj if (mask(thisi,thisj) /= 0) mask(thisi,thisj) = value do do i=1,nmove thisi=iadd(thisi,movei(i),imt) thisj=jadd(thisj,movej(i),jmt) kcnt = kcnt + 1 ! write(6,*) thisi,thisj,mask(thisi,thisj) if (mask(thisi,thisj) /= 0) then mask(thisi,thisj) = value lcnt = lcnt+1 endif if (thisi==eni .and. thisj==enj) exit enddo if (thisi==eni .and. thisj==enj) exit if (kcnt > 400) stop 'not getting anywhere' enddo ! write(6,*) 'assigned region value to ',lcnt,' cells' deallocate(movei) deallocate(movej) enddo ! k, loop over nb boundary segments END SUBROUTINE set_bounds !*********************************************************************** integer FUNCTION iadd(i,di,ni) implicit none !--- arguments --- integer :: i, di, ni if (di < 0) then iadd = mod(i+di-1+ni,ni)+1 else iadd = mod(i,ni)+di endif END FUNCTION iadd !=============================================================================== integer FUNCTION jadd(j,dj,nj) implicit none !--- arguments --- integer :: j, dj, nj if (dj < 0) then jadd = max(j+dj,1) else jadd = min(j+dj,nj) endif END FUNCTION jadd !=============================================================================== SUBROUTINE get_moves(nmove,idx,idy,movei,movej) !----------------------------------------------------------------------- ! ! For a slope in i,j grid space given by idy & idx, this ! routine determines a set of adjacent cell moves to ! accomplish this slope in the most linear fashion possible. ! The moves are stored as 1-D arrays movei, movej which ! hold values of -1, 0, or 1. ! !----------------------------------------------------------------------- integer :: i, j , nmove, idx, idy, nidx, nidy, cnt integer, dimension(nmove) :: tmp1, tmp2, movei, movej movei = 0 movej = 0 cnt = 0 nidx = abs(idx) nidy = abs(idy) ! cannot improve the linearity when either idx or idy <= 1 if (min(nidx,nidy) <= 1) then do j=1,nidy cnt = cnt + 1 movej(cnt)=idy/nidy end do do i=1,nidx cnt = cnt + 1 movei(cnt)=idx/nidx end do return end if ! if (nidy>nidx) then slp = nint(real(nidy)/real(nidx)) nfull = min(int(nidy/slp),nidx) nxleft = nidx-nfull nyleft = nidy-(slp*nfull) cnt = 0 do i=1,nfull cnt = cnt + 1 movei(cnt) = idx/nidx do j=1,slp cnt = cnt + 1 movej(cnt) = idy/nidy end do end do do i=1,nxleft cnt = cnt + 1 movei(cnt) = idx/nidx end do do j=1,nyleft cnt = cnt + 1 movej(cnt) = idy/nidy end do else ! nidx>nidy since they can't be equal slp = nint(real(nidx)/real(nidy)) nfull = min(int(nidx/slp),nidy) nyleft = nidy-nfull nxleft = nidx-(slp*nfull) cnt = 0 do i=1,nfull cnt = cnt + 1 movej(cnt) = idy/nidy do j=1,slp cnt = cnt + 1 movei(cnt) = idx/nidx end do end do do i=1,nyleft cnt = cnt + 1 movej(cnt) = idy/nidy end do do j=1,nxleft cnt = cnt + 1 movei(cnt) = idx/nidx end do end if END SUBROUTINE get_moves !*********************************************************************** SUBROUTINE calc_tpoints(imt,jmt,ux,uy,tlng,tlat) !----------------------------------------------------------------------- ! ! calculates lat/lon coordinates of T points from U points. ! this computation is identical to the method used in the ccsm ! ice model subroutine Tlatlon in module ice_gird.F as of ! 14 September 2001 ! !----------------------------------------------------------------------- integer :: i, j , im1, imt, jmt real*8, parameter :: pi=3.14159265358979323846 real*8 :: z1,x1,y1,z2,x2,y2,z3,x3,y3,z4,x4,y4,tx,ty,tz,da real*8 :: c4 = 4.0 real*8, intent(in), dimension(imt,jmt) :: ux, uy real*8, intent(out), dimension(imt,jmt) :: tlat, tlng write(6,*) 'imt, jmt = ',imt,jmt write(6,*) 'ux[50,50] = ',ux(50,50)*(180./pi) write(6,*) 'uy[50,50] = ',uy(50,50)*(180./pi) do j=2,jmt do i=1,imt if (i.eq.1) then im1=imt else im1=i-1 endif z1 = cos(uy(im1,j-1)) x1 = cos(ux(im1,j-1))*z1 y1 = sin(ux(im1,j-1))*z1 z1 = sin(uy(im1,j-1)) z2 = cos(uy(i,j-1)) x2 = cos(ux(i,j-1))*z2 y2 = sin(ux(i,j-1))*z2 z2 = sin(uy(i,j-1)) z3 = cos(uy(im1,j)) x3 = cos(ux(im1,j))*z3 y3 = sin(ux(im1,j))*z3 z3 = sin(uy(im1,j)) z4 = cos(uy(i,j)) x4 = cos(ux(i,j))*z4 y4 = sin(ux(i,j))*z4 z4 = sin(uy(i,j)) tx = (x1+x2+x3+x4)/c4 ty = (y1+y2+y3+y4)/c4 tz = (z1+z2+z3+z4)/c4 da = sqrt(tx**2+ty**2+tz**2) tz = tz/da ! tlng in radians tlng(i,j) = 0.0 if (tx.ne.0.0.or.ty.ne.0.0) & tlng(i,j) = atan2(ty,tx) ! tlat in radians tlat(i,j) = asin(tz) end do end do ! j=1: linear approximation do i=1,imt tlng(i,1) = tlng(i,2) tlat (i,1) = 2.0*tlat(i,2) - tlat(i,3) end do where (tlng > 2.*pi) tlng = tlng - 2.*pi where (tlng < 0.0 ) tlng = tlng + 2.*pi ! convert to degrees tlat = tlat*(180./pi) tlng = tlng*(180./pi) write(6,*) 'tx[50,50] = ',tlng(50,50) write(6,*) 'ty[50,50] = ',tlat(50,50) !----------------------------------------------------------------------- END SUBROUTINE calc_tpoints