module mod_mercatorgrid integer :: dimlenx, dimleny, dimlenz real, allocatable, dimension(:,:) :: merclon , merclat ,mercbathy , mercmask ! T-cell centers real, allocatable, dimension(:,:) :: merclonu, merclatu, mercmasku ! U-cell real, allocatable, dimension(:,:) :: merclonv, merclatv, mercmaskv ! V-cell character(len=*), parameter :: clon ='gphit', clat ='glamt' character(len=*), parameter :: clonu ='gphiu', clatu='glamu' character(len=*), parameter :: clonv ='gphiv', clatv='glamv' character(len=*), parameter :: cbathy='hdept' character(len=*), parameter :: cmask ='tmask' character(len=*), parameter :: cmasku='umask' character(len=*), parameter :: cmaskv='vmask' character(len=*), parameter :: bathyfile='ext-mesh_mask.nc' contains subroutine mercgrid_init() use netcdf use m_spherdist implicit none integer :: ierr integer :: k integer :: ncid integer :: dimidx, dimidy, dimidz integer :: tempid,zvarid,lonid,latid, tmpid integer, allocatable :: ip(:,:),tmpint(:,:) integer, dimension(nf90_max_var_dims) :: dimIDs integer :: i2,j2,i,j integer :: ndims real :: zlev(1), amax, amin logical :: skipip ! Read Mercator netcdf file ierr=NF90_OPEN(bathyfile,NF90_NOWRITE,ncid) ! Get dims ierr=nf90_inq_dimid(ncid,'x',dimidx) ierr=nf90_inquire_dimension(ncid,dimidx,len=dimlenx) if (ierr/=NF90_NOERR) stop '(mercgrid:error getting x dim)' print *,dimidx,dimlenx ierr=nf90_inq_dimid(ncid,'y',dimidy) ierr=nf90_inquire_dimension(ncid,dimidy,len=dimleny) if (ierr/=NF90_NOERR) stop '(mercgrid:error getting y dim)' print *,dimidy,dimleny ierr=nf90_inq_dimid(ncid,'z',dimidz) ierr=nf90_inquire_dimension(ncid,dimidz,len=dimlenz) if (ierr/=NF90_NOERR) stop '(mercgrid:error getting z dim)' print *,dimidz,dimlenz ! Start allocating allocate(merclon (dimlenx,dimleny)) allocate(merclonu (dimlenx,dimleny)) allocate(merclonv (dimlenx,dimleny)) allocate(merclat (dimlenx,dimleny)) allocate(merclatu (dimlenx,dimleny)) allocate(merclatv (dimlenx,dimleny)) allocate(mercmask (dimlenx,dimleny)) allocate(mercmasku(dimlenx,dimleny)) ! NB - no z dim allocate(mercmaskv(dimlenx,dimleny)) ! NB - no z dim !!!!!!!!!! Longitudes ! Get Mercator lon (T) ierr=nf90_inq_varid(ncid,trim(clon),lonid) if (ierr==NF90_NOERR) then ierr=nf90_get_var(ncid,lonid,merclon) if (ierr/=NF90_NOERR) stop '(mercgrid:error getting t lon)' print *,'merc lon (T):',minval(merclon),maxval(merclon) else print *,'Error on reading mercator grid - t lon' call exit(1) end if ! Get Mercator lon (U) ierr=nf90_inq_varid(ncid,trim(clonu),lonid) if (ierr==NF90_NOERR) then ierr=nf90_get_var(ncid,lonid,merclonu) if (ierr/=NF90_NOERR) stop '(mercgrid:error getting u lon)' print *,'merc lon (U):',minval(merclonu),maxval(merclonu) else print *,'Error on reading mercator grid - u lon' call exit(1) end if ! Get Mercator lon (V) ierr=nf90_inq_varid(ncid,trim(clonv),lonid) if (ierr==NF90_NOERR) then ierr=nf90_get_var(ncid,lonid,merclonv) if (ierr/=NF90_NOERR) stop '(mercgrid:error getting V lon)' print *,'merc lon (V):',minval(merclonv),maxval(merclonv) else print *,'Error on reading mercator grid - V lon' call exit(1) end if !!!!!!!!!! Latitudes ! Get Mercator lat (T) ierr=nf90_inq_varid(ncid,trim(clat),latid) if (ierr==NF90_NOERR) then ierr=nf90_get_var(ncid,latid,merclat) if (ierr/=NF90_NOERR) stop '(mercgrid:error getting t lat)' print *,'merc lat (T):',minval(merclat),maxval(merclat) else print *,'Error on reading mercator grid - t lat' call exit(1) end if ! Get Mercator lat (U) ierr=nf90_inq_varid(ncid,trim(clatu),latid) if (ierr==NF90_NOERR) then ierr=nf90_get_var(ncid,latid,merclatu) if (ierr/=NF90_NOERR) stop '(mercgrid:error getting U lat)' print *,'merc lat (U):',minval(merclatu),maxval(merclatu) else print *,'Error on reading mercator grid - U lat' call exit(1) end if ! Get Mercator lat (V) ierr=nf90_inq_varid(ncid,trim(clatv),latid) if (ierr==NF90_NOERR) then ierr=nf90_get_var(ncid,latid,merclatv) if (ierr/=NF90_NOERR) stop '(mercgrid:error getting V lat)' print *,'merc lat (V):',minval(merclatv),maxval(merclatv) else print *,'Error on reading mercator grid - V lat' call exit(1) end if !!!!!!!!!! The rest ! Get Mercator bathymetry (T) ierr=nf90_inq_varid(ncid,trim(cbathy),tmpid) if (ierr==NF90_NOERR) then ierr=nf90_get_var(ncid,tmpid,mercbathy) if (ierr/=NF90_NOERR) stop '(mercgrid:error getting t bathy)' print *,'merc bathy (T):',minval(mercbathy),maxval(mercbathy) else print *,'Error on reading mercator grid - t bathy' call exit(1) end if ierr=nf90_close(ncid) end subroutine subroutine mercgrid_mask(k) use netcdf implicit none integer, intent(in) :: k ! vertical level integer :: ierr, tmpid, ncid ! Read Mercator netcdf file ierr=NF90_OPEN(bathyfile,NF90_NOWRITE,ncid) ! Get Mercator mask (T) ierr=nf90_inq_varid(ncid,trim(cmask),tmpid) if (ierr==NF90_NOERR) then ierr=nf90_get_var(ncid,tmpid,mercmask,start=(/1,1,k,1/)) if (ierr/=NF90_NOERR) stop '(mercgrid_mask:error getting t mask)' print *,'merc mask (T):',minval(mercmask),maxval(mercmask) else print *,'Error on reading mercator grid - T mask' call exit(1) end if ! Get Mercator mask (U) ierr=nf90_inq_varid(ncid,trim(cmasku),tmpid) if (ierr==NF90_NOERR) then ierr=nf90_get_var(ncid,tmpid,mercmasku,start=(/1,1,k,1/)) if (ierr/=NF90_NOERR) stop '(mercgrid_mask:error getting u mask)' print *,'merc mask (U):',minval(mercmasku),maxval(mercmasku) else print *,'Error on reading mercator grid - T mask' call exit(1) end if ! Get Mercator mask (V) ierr=nf90_inq_varid(ncid,trim(cmaskv),tmpid) if (ierr==NF90_NOERR) then ierr=nf90_get_var(ncid,tmpid,mercmaskv,start=(/1,1,k,1/)) if (ierr/=NF90_NOERR) stop '(mercgrid_mask:error getting v mask)' print *,'merc mask (V):',minval(mercmaskv),maxval(mercmaskv) else print *,'Error on reading mercator grid - V mask' call exit(1) end if end subroutine mercgrid_mask subroutine mercgrid_pivots(merclon,merclat,hyclon,hyclat,ipiv,jpiv, & i4,j4,dist4,idm,jdm,usefile) use m_spherdist use mod_za, only : zaiowr, zaiocl, zaiopf, zaiord implicit none integer, intent(in) :: idm,jdm real, dimension(dimlenx,dimleny) :: merclon, merclat real, dimension(idm,jdm) :: hyclon, hyclat integer, dimension(idm,jdm) :: ipiv,jpiv logical, intent(in) :: usefile real, dimension(idm,jdm,4) :: dist4 integer, dimension(idm,jdm,4) :: i4,j4 logical ::skipip, match integer :: i,j,i2,j2 real :: dist, mindist,tmpr(idm,jdm),amin,amax integer :: ipivtmp,jpivtmp,ip(idm,jdm) real :: mindist4(4) integer :: idist, idist2 integer :: mini(4), minj(4) integer :: k ! Go through hycom grid - find point closest point on mercator grid if (usefile) then ipiv=1 jpiv=1 print *,'NB!! Using stored pivots !!!' print *,'NB!! Using stored pivots !!!' print *,'NB!! Using stored pivots !!!' print *,'NB!! Using stored pivots !!!' print *,'NB!! Using stored pivots !!!' print *,'NB!! Using stored pivots !!!' print *,'NB!! Using stored pivots !!!' print *,'NB!! Using stored pivots !!!' print *,'NB!! Using stored pivots !!!' print *,'NB!! Using stored pivots !!!' print *,'NB!! Using stored pivots !!!' call zaiopf('pivots.a','old',66) call zaiord(tmpr,ip,.false.,amin,amax,66) !plon ipiv=int(tmpr) call zaiord(tmpr,ip,.false.,amin,amax,66) !plon jpiv=int(tmpr) do k=1,4 call zaiord(tmpr,ip,.false.,amin,amax,66) !plon i4(:,:,k)=nint(tmpr) call zaiord(tmpr,ip,.false.,amin,amax,66) !plon j4(:,:,k)=nint(tmpr) call zaiord(tmpr,ip,.false.,amin,amax,66) !plon dist4(:,:,k)=tmpr end do call zaiocl(66) else do iedge=1,3 ! 1=t,2=u,3=v (not exactly edge but... do j=1,jdm print *,j do i=1,idm mindist=1e14 mindist4=1e14 ipiv(i,j)=-1 jpiv(i,j)=-1 do j2=1,dimleny do i2=1,dimlenx if (iedge==1) then dist=spherdist(hyclon(i,j),hyclat(i,j),merclon(i2,j2),merclat(i2,j2)) if (dist < mindist) then ipivtmp=i2 jpivtmp=j2 mindist=dist !print *,i2,j2,mindist end if if (any(mindist4>dist)) then match=.false. idist=0 do while (.not.match) idist=idist+1 if (dist < mindist4(idist)) match=.true. end do !print *,'idist,dist is ',idist,dist !print *,'mindist4 is ',mindist4 do idist2=4,max(idist,2),-1 ! Move stuff down a step if (idist2<4) then mini (idist2)=mini (idist2-1) minj (idist2)=minj (idist2-1) mindist4(idist2)=mindist4(idist2-1) end if end do ! This one should be placed here mini (idist)=i2 minj (idist)=j2 mindist4(idist)=dist end if end do end do !print *,'mindist4 is ',mindist4 !print *,'sum mindist4 is ',sum(mindist4) !print *,'mini is ',mini !print *,'minj is ',minj ipiv(i,j)=ipivtmp jpiv(i,j)=jpivtmp !print *,i,j,ipivtmp,jpivtmp,hyclon(i,j),hyclat(i,j),merclon(ipivtmp,jpivtmp),merclat(ipivtmp,jpivtmp) dist4(i,j,:)=mindist4 i4 (i,j,:)=mini j4 (i,j,:)=minj end do end do call zaiopf('pivots.a','replace',66) tmpr=real(ipiv) call zaiowr(tmpr,ip,.false.,amin,amax,66,.false.) !plon tmpr=real(jpiv) call zaiowr(tmpr,ip,.false.,amin,amax,66,.false.) !plon do k=1,4 tmpr=real(i4(:,:,k)) call zaiowr(tmpr,ip,.false.,amin,amax,66,.false.) !plon tmpr=real(j4(:,:,k)) call zaiowr(tmpr,ip,.false.,amin,amax,66,.false.) !plon tmpr=dist4(:,:,k) call zaiowr(tmpr,ip,.false.,amin,amax,66,.false.) !plon end do call zaiocl(66) endif end subroutine mercgrid_pivots end module mod_mercatorgrid