#if ! defined(RELO) subroutine latbdf(n,lll) use mod_xc ! HYCOM communication interface use mod_cb_arrays ! HYCOM saved arrays use mod_tides ! HYCOM tides implicit none ! integer n,lll ! ! --- apply lateral boundary conditions to barotropic flow field ! ! --- port flow version: ! --- NOT similar to the standard 'Browning and Kreiss' MICOM/HYCOM open ! --- boundary condition. This version uses algorithms based on a ! --- 1 invariant Flather boundary condition (setting the gradient ! --- of the incoming characteristic to zero). ! ! --- The tangential velocity is not constrained. ! ! --- see also: latbdp ! ! --- the code is as similar as possible to that for the standard case. ! --- so for example, 'speed' is in fact 1/SQRT(gH) which represents ! --- c1/g in the notation of (Bleck and Sun, Open boundary conditions ! --- for MICOM). The 1/g allows for the use of pressure fields. ! ! --- Note that East, West, North and South refers to the grid ! --- (i.e i,j points) and NOT geographic East, West, North and South ! ! --- the first call is made during initialization. ! ! --- Iris Lohmann, Carlos Lozano, NCEP, April 2006 ! logical, parameter :: ldebug_latbdf=.false. ! integer, parameter :: mports=1 !maximum number of ports ! integer, parameter :: nchar =120 ! logical lfatal,lfatalp integer i,j,isec,ifrst,ilast,l real aline(nchar), & dline(itdm+jtdm),xline(itdm+jtdm), & pline(itdm+jtdm),uline(itdm+jtdm) real sum,svspin,fatal character*3 char3 ! integer nports integer lnport(mports),kdport(mports) integer jfport(mports),jlport(mports), & ifport(mports),ilport(mports) real svpnow(mports),svport(mports) save lnport save svpnow,svport save nports,kdport,ifport,ilport,jfport,jlport ! real uportw(jtdm,mports),speedw(jtdm,mports),rspedw(jtdm,mports), & uporte(jtdm,mports),speede(jtdm,mports),rspede(jtdm,mports), & vportn(itdm,mports),speedn(itdm,mports),rspedn(itdm,mports), & vports(itdm,mports),speeds(itdm,mports),rspeds(itdm,mports) save uportw,speedw,rspedw,uporte,speede,rspede, & vportn,speedn,rspedn,vports,speeds,rspeds ! tides stuff integer npts_p,kdpt_p(mports), & ifpt_p(mports),ilpt_p(mports), & jfpt_p(mports),jlpt_p(mports),lnpt_p(mports) integer npts_v,kdpt_v(mports), & ifpt_v(mports),ilpt_v(mports), & jfpt_v(mports),jlpt_v(mports),lnpt_v(mports) integer npts_u,kdpt_u(mports), & ifpt_u(mports),ilpt_u(mports), & jfpt_u(mports),jlpt_u(mports),lnpt_u(mports) ! the max of itdm and jtdm is ok for any port ! number of tidal consituents (ncon) from mod_tides real z1r_p(mports,max(jtdm,itdm),ncon) real z1i_p(mports,max(jtdm,itdm),ncon) real z1r_u(mports,max(jtdm,itdm),ncon) real z1i_u(mports,max(jtdm,itdm),ncon) real z1r_v(mports,max(jtdm,itdm),ncon) real z1i_v(mports,max(jtdm,itdm),ncon) real tmp1, tmp2 real tmpr(max(itdm,jtdm),ncon), tmpi(max(itdm,jtdm),ncon) real upred(max(itdm,jtdm)),zpred(max(itdm,jtdm)) real udpred(max(itdm,jtdm)) real vpred(max(itdm,jtdm)) real ulow(max(itdm,jtdm),mports),plow(max(itdm,jtdm),mports) real uu(max(itdm,jtdm)),vv(max(itdm,jtdm)) integer bnd_init(mports) real*8 d_time real*8 timermp,frmp logical astroflag integer jn,in,ic save z1r_p, z1i_p save z1r_u, z1i_u save bnd_init save ulow,plow ! character*13 fmt save fmt data fmt / '(i4,1x,120i1)' / ! ! USER-INPUT: optimization coefficients for the 1 inv algorithm. real w_1,w_1c save w_1 data w_1 / 0.1 / ! integer lcount save lcount data lcount / 0 / ! ! Comment out the next line to run without ! boundary tides d_time = time_8 + lll*dlt/86400.d0 ! add 15384 for obtaining mjd lcount = lcount + 1 ! ! set 1-invariant coefficient w_1c=1.0-w_1 ! ! ! --- the first call just initializes data structures. ! if (lcount.eq.1) then ! open(unit=uoff+99,file=trim(flnminp)//'ports.input') ! ! --- 'nports' = number of boundary port sections. call blkini(nports,'nports') if (mnproc.eq.1) then write(lp,*) endif if (nports.eq.0) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdf - must have lbflag=0 for nports=0' write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbdf)') stop '(latbdf)' elseif (nports.lt.0 .or. nports.gt.mports) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdf - illegal nports value' if (nports.gt.mports) then write(lp,*) 'increase parameter mports to',nports endif write(lp,*) call flush(lp) endif call xcstop('(latbdf)') stop '(latbdf)' endif ! ! --- read in the ports one at a time ! do l= 1,nports ! ! --- port location is w.r.t. u (EW) or v (NS) grid ! --- and identifies the sea at the port ! --- the minimum index is 0 ! ! --- 'kdport' = port orientation (1=N, 2=S, 3=E, 4=W, 5=S(Fundy)) ! --- 'ifport' = first i-index ! --- 'ilport' = last i-index (=ifport for N or S orientation) ! --- 'jfport' = first j-index ! --- 'jlport' = last j-index (=jfport for E or W orientation) ! --- 'svpnow' = existing port transport in Sv (+ve towards E or S) ! --- 'svport' = target port transport in Sv (+ve towards E or S) ! --- 'lnport' = port length (calculated, not input) call blkini(kdport(l),'kdport') call blkini(ifport(l),'ifport') call blkini(ilport(l),'ilport') call blkini(jfport(l),'jfport') call blkini(jlport(l),'jlport') call blkinr(svpnow(l),'svpnow','(a6," =",f10.4," Sv")') call blkinr(svport(l),'svport','(a6," =",f10.4," Sv")') if (mnproc.eq.1) then write(lp,*) endif ! lnport(l) = ilport(l)-ifport(l)+jlport(l)-jfport(l)+1 ! ! --- sanity check. ! if (kdport(l).eq.3.or.kdport(l).eq.4) then if (ifport(l).ne.ilport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdp - port direction', & ' and orientation are not consistent' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif else if (jfport(l).ne.jlport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdp - port direction', & ' and orientation are not consistent' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif endif if (ifport(l).gt.ilport(l) .or. & jfport(l).gt.jlport(l) ) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdp - port', & ' location is not consistent' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif enddo ! close(unit=uoff+99) ! ********************************************** !------ read low frequency input-file ! initialize the low-frequency u and eta to zero do i = 1,max(itdm,jtdm) do j = 1,mports ulow(i,j) =0.0 plow(i,j) =0.0 enddo enddo open(unit=uoff+99,file=trim(flnminp)//'lowfreq.input') ! the file is read in order west, east, south, north do l= 1,nports if (kdport(l).eq.4) then ! western port do j= jfport(l),jlport(l) read(uoff+99,*) ulow(j,l),plow(j,l) enddo elseif (kdport(l).eq.3) then ! eastern port do j= jfport(l),jlport(l) read(uoff+99,*) ulow(j,l),plow(j,l) enddo elseif (kdport(l).eq.2) then ! southern port do i= ifport(l),ilport(l) read(uoff+99,*) ulow(i,l),plow(i,l) enddo elseif (kdport(l).eq.1) then ! northern port do i= ifport(l),ilport(l) read(uoff+99,*) ulow(i,l),plow(i,l) enddo endif !kdport enddo !l=1,nports close(uoff+99) !************************************************ ! ****** READ TIDAL CONSTITUENTS **************** if (tidflg.ge.1) then ! initialize the tidal constituents to zero do i = 1,mports do j = 1,max(itdm,jtdm) do ic = 1,ncon z1r_p(i,j,ic) = 0. z1i_p(i,j,ic) = 0. z1r_u(i,j,ic) = 0. z1i_u(i,j,ic) = 0. z1r_v(i,j,ic) = 0. z1i_v(i,j,ic) = 0. enddo enddo enddo do j = 1,max(itdm,jtdm) uu(j) = 0. vv(j) = 0. upred(j) = 0. vpred(j) = 0. zpred(j) = 0. enddo ! -------------------------------------------------------- ! Now the P points open(unit=uoff+99,file=trim(flnminp)//'tidalports_p.input') ! --- 'nports' = number of boundary port sections. call blkini(npts_p,'npts_p') if (mnproc.eq.1) then write(lp,*) endif if (npts_p.ne.nports) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error number of ports needs to be same' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif do l= 1,nports ! ! --- port location is w.r.t. u (EW) or v (NS) grid ! --- and identifies the sea at the port ! --- the minimum index is 0 ! ! --- 'kdport' = port orientation (1=N, 2=S, 3=E, 4=W, 5=Fundy(S)) ! --- 'ifport' = first i-index ! --- 'ilport' = last i-index (=ifport for N or S orientation) ! --- 'jfport' = first j-index ! --- 'jlport' = last j-index (=jfport for E or W orientation) ! --- 'lnport' = port length (calculated, not input) call blkini(kdpt_p(l),'kdpt_p') if (kdpt_p(l).ne.kdport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'Mismatch of the kdport' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif call blkini(ifpt_p(l),'ifpt_p') if (ifpt_p(l).ne.ifport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'Mismatch of the ifport' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif call blkini(ilpt_p(l),'ilpt_p') if (ilpt_p(l).ne.ilport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'Mismatch of the ilport' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif call blkini(jfpt_p(l),'jfpt_p') if (jfpt_p(l).ne.jfport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'Mismatch of the jfport' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif call blkini(jlpt_p(l),'jlpt_p') if (jlpt_p(l).ne.jlport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'Mismatch of the jlport' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif if (kdport(l).eq.4) then ! ! western port ! do j= jfport(l),jlport(l) do ic = 1, ncon read(uoff+99,*) tmp1, tmp2 z1r_p(l,j,ic) = tmp1 z1i_p(l,j,ic) = tmp2 enddo enddo ! elseif (kdport(l).eq.3) then ! ! eastern port do j= jfport(l),jlport(l) do ic = 1, ncon read(uoff+99,*) tmp1, tmp2 z1r_p(l,j,ic) = tmp1 z1i_p(l,j,ic) = tmp2 enddo enddo ! elseif (kdport(l).eq.1) then ! ! northern port ! do i= ifport(l),ilport(l) do ic = 1, ncon read(uoff+99,*) tmp1, tmp2 z1r_p(l,i,ic) = tmp1 z1i_p(l,i,ic) = tmp2 enddo enddo ! elseif (kdport(l).eq.2) then ! ! southern port ! do i= ifport(l),ilport(l) do ic = 1, ncon read(uoff+99,*) tmp1, tmp2 z1r_p(l,i,ic) = tmp1 z1i_p(l,i,ic) = tmp2 enddo enddo endif !kdport= ! Close the l = 1, nports loop enddo !nports close(uoff+99) ! ------------------------------------------------------------- ! Now the normal-velocity points open(unit=uoff+99,file=trim(flnminp)//'tidalports_v.input') ! --- 'nports' = number of boundary port sections. call blkini(npts_u,'npts_u') if (mnproc.eq.1) then write(lp,*) endif if (npts_u.ne.nports) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error number of ports needs to be same' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif do l= 1,nports ! ! --- port location is w.r.t. u (EW) or v (NS) grid ! --- and identifies the sea at the port ! --- the minimum index is 0 ! ! --- 'kdport' = port orientation (1=N, 2=S, 3=E, 4=W, 5=Fundy(S)) ! --- 'ifport' = first i-index ! --- 'ilport' = last i-index (=ifport for N or S orientation) ! --- 'jfport' = first j-index ! --- 'jlport' = last j-index (=jfport for E or W orientation) ! --- 'svpnow' = existing port transport in Sv (+ve towards E or S) ! --- 'svport' = target port transport in Sv (+ve towards E or S) ! --- 'lnport' = port length (calculated, not input) call blkini(kdpt_u(l),'kdpt_u') if (kdpt_u(l).ne.kdport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'Mismatch of the kdport' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif call blkini(ifpt_u(l),'ifpt_u') if (ifpt_u(l).ne.ifport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'Mismatch of the ifport' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif call blkini(ilpt_u(l),'ilpt_u') if (ilpt_u(l).ne.ilport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'Mismatch of the ilport' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif call blkini(jfpt_u(l),'jfpt_u') if (jfpt_u(l).ne.jfport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'Mismatch of the jfport' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif call blkini(jlpt_u(l),'jlpt_u') if (jlpt_u(l).ne.jlport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'Mismatch of the jlport' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif if (kdport(l).eq.4) then ! ! western port ! do j= jfport(l),jlport(l) do ic = 1, ncon read(uoff+99,*) tmp1, tmp2 z1r_u(l,j,ic) = tmp1 z1i_u(l,j,ic) = tmp2 enddo enddo ! elseif (kdport(l).eq.3) then ! ! eastern port do j= jfport(l),jlport(l) do ic = 1, ncon read(uoff+99,*) tmp1, tmp2 z1r_u(l,j,ic) = tmp1 z1i_u(l,j,ic) = tmp2 enddo enddo ! elseif (kdport(l).eq.1) then ! ! northern port ! do i= ifport(l),ilport(l) do ic = 1, ncon read(uoff+99,*) tmp1, tmp2 z1r_u(l,i,ic) = tmp1 z1i_u(l,i,ic) = tmp2 enddo enddo ! elseif (kdport(l).eq.2) then ! ! southern port ! do i= ifport(l),ilport(l) do ic = 1, ncon read(uoff+99,*) tmp1, tmp2 z1r_u(l,i,ic) = tmp1 z1i_u(l,i,ic) = tmp2 enddo enddo endif !kdport= ! Close the l = 1, nports loop enddo !nports close(uoff+99) !*****************END OF READING THE TIDAL CONSTITUENTS endif !tidflg.ge.1 ! ! --- check ports against masks, ! --- mark the port locations on masks and print them out. ! lfatal = .false. do l= 1,nports lfatalp = .false. ! if (kdport(l).eq.4) then ! ! western port ! i = ifport(l) do j= jfport(l),jlport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.1 .or. j.gt.jtdm ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iu(i-i0,j-j0).ne.0) then lfatalp = .true. iu(i-i0,j-j0) = 9 !indicate an error else iu(i-i0,j-j0) = -1 endif if (iu(i-i0+1,j-j0).ne.1 .or. & iu(i-i0+2,j-j0).ne.1 ) then lfatalp = .true. iu(i-i0,j-j0) = 7 !indicate an error endif enddo ! elseif (kdport(l).eq.3) then ! ! eastern port ! i = ifport(l) do j= jfport(l),jlport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.1 .or. j.gt.jtdm ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iu(i-i0,j-j0).ne.0) then lfatalp = .true. iu(i-i0,j-j0) = 9 !indicate an error else iu(i-i0,j-j0) = -1 endif if (iu(i-i0-1,j-j0).ne.1 .or. & iu(i-i0-2,j-j0).ne.1 ) then lfatalp = .true. iu(i-i0,j-j0) = 7 !indicate an error endif enddo ! elseif (kdport(l).eq.1) then ! ! northern port ! j = jfport(l) do i= ifport(l),ilport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.3 .or. j.gt.jtdm ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iv(i-i0,j-j0).ne.0) then lfatalp = .true. iv(i-i0,j-j0) = 9 !indicate an error else iv(i-i0,j-j0) = -1 endif if (iv(i-i0,j-j0-1).ne.1 .or. & iv(i-i0,j-j0-2).ne.1 ) then lfatalp = .true. iv(i-i0,j-j0) = 7 !indicate an error endif enddo ! elseif (kdport(l).eq.2.or.kdport(l).eq.5) then ! ! southern port ! j = jfport(l) do i= ifport(l),ilport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.1 .or. j.gt.jtdm-2 ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iv(i-i0,j-j0).ne.0) then lfatalp = .true. iv(i-i0,j-j0) = 9 !indicate an error else iv(i-i0,j-j0) = -1 endif if (iv(i-i0,j-j0+1).ne.1 .or. & iv(i-i0,j-j0+2).ne.1 ) then lfatalp = .true. iv(i-i0,j-j0) = 7 !indicate an error endif enddo ! endif !kdport= ! if (lfatalp) then write(lp,*) write(lp,*) 'error in latbdp - port ',l,' mislocated', & ' (mnproc = ',mnproc,')' write(lp,*) call flush(lp) endif lfatal = lfatal .or. lfatalp enddo !nports ! ! local lfatal to global lfatal ! if (lfatal) then fatal = 1.0 else fatal = 0.0 endif call xcmaxr(fatal) lfatal = fatal.gt.0.5 ! ! --- write out -iu- and -iv- arrays, if they are not too big ! --- data are written in strips nchar points wide if (lfatal .or. max(itdm,jtdm).le.2*nchar) then util1(1:ii,1:jj) = iu(1:ii,1:jj) ! xclget is for real arrays isec=(itdm-1)/nchar do ifrst=0,nchar*isec,nchar ilast=min(itdm,ifrst+nchar) write (char3,'(i3)') ilast-ifrst fmt(8:10)=char3 if (mnproc.eq.1) then write (lp,'(a,i5,a,i5)') & 'iu array, cols',ifrst+1,' --',ilast endif do j= jtdm,1,-1 call xclget(aline,ilast-ifrst, util1,ifrst+1,j,1,0, 1) if (mnproc.eq.1) then write (lp,fmt) j,(nint(aline(i)),i=1,ilast-ifrst) endif enddo enddo if (mnproc.eq.1) then write (lp,*) endif call xcsync(flush_lp) ! util1(1:ii,1:jj) = iv(1:ii,1:jj) ! xclget is for real arrays isec=(itdm-1)/nchar do ifrst=0,nchar*isec,nchar ilast=min(itdm,ifrst+nchar) write (char3,'(i3)') ilast-ifrst fmt(8:10)=char3 if (mnproc.eq.1) then write (lp,'(a,i5,a,i5)') & 'iv array, cols',ifrst+1,' --',ilast endif do j= jtdm,1,-1 call xclget(aline,ilast-ifrst, util1,ifrst+1,j,1,0, 1) if (mnproc.eq.1) then write (lp,fmt) j,(nint(aline(i)),i=1,ilast-ifrst) endif enddo enddo if (mnproc.eq.1) then write (lp,*) endif call xcsync(flush_lp) endif ! small region ! if (lfatal) then write(lp,*) write(lp,*) 'error in latbdp - bad port(s)' write(lp,*) call flush(lp) call xchalt('(latbdp)') stop '(latbdp)' endif ! ! --- restore iu and iv, and zero iuopn and ivopn. ! ! $OMP PARALLEL DO PRIVATE(j,i) ! $OMP& SCHEDULE(STATIC,jblk) do j= 1,jj do i= 1,ii iu(i,j) = max( iu(i,j), 0 ) iv(i,j) = max( iv(i,j), 0 ) enddo enddo ! $OMP PARALLEL DO PRIVATE(j,i) ! $OMP& SCHEDULE(STATIC,jblk) do j= 1-nbdy,jj+nbdy do i= 1-nbdy,ii+nbdy iuopn(i,j) = 0 ivopn(i,j) = 0 enddo enddo ! ! --- initialize the ports ! do l= 1,nports if (kdport(l).eq.4) then ! ! western port ! sum = 0.0 i = ifport(l) j = jfport(l) call xclget(dline(j),lnport(l), depths,i,j,0,1, 0) call xclget(xline(j),lnport(l), scuy, i, j,0,1, 0) do j= jfport(l),jlport(l) sum = sum + dline(j)*xline(j) enddo sum = 1.e6/sum do j= jfport(l),jlport(l) uportw(j,l) = sum speedw(j,l) = sqrt(1.0*svref/(onem*dline(j))) rspedw(j,l) = 1.0/speedw(j,l) if (i.ge.i0+ 1-nbdy .and. & i.le.i0+ii+nbdy .and. & j.ge.j0+ 1-nbdy .and. & j.le.j0+jj+nbdy ) then iuopn(i-i0,j-j0) = 1 endif enddo ! elseif (kdport(l).eq.3) then ! ! eastern port ! sum = 0.0 i = ifport(l)-1 j = jfport(l) call xclget(dline(j),lnport(l), depths,i, j,0,1, 0) call xclget(xline(j),lnport(l), scuy, i+1,j,0,1, 0) do j= jfport(l),jlport(l) sum = sum + dline(j)*xline(j) enddo sum = 1.e6/sum do j= jfport(l),jlport(l) uporte(j,l) = sum speede(j,l) = sqrt(1.0*svref/(onem*dline(j))) rspede(j,l) = 1.0/speede(j,l) if (i+1.ge.i0+ 1-nbdy .and. & i+1.le.i0+ii+nbdy .and. & j .ge.j0+ 1-nbdy .and. & j .le.j0+jj+nbdy ) then iuopn(i-i0+1,j-j0) = 1 endif enddo ! elseif (kdport(l).eq.1) then ! ! northern port ! sum = 0.0 j = jfport(l)-1 i = ifport(l) call xclget(dline(i),lnport(l), depths,i,j, 1,0, 0) call xclget(xline(i),lnport(l), scuy, i,j+1,1,0, 0) do i= ifport(l),ilport(l) sum = sum + dline(i)*xline(i) enddo sum = 1.e6/sum do i= ifport(l),ilport(l) vportn(i,l) = sum speedn(i,l) = sqrt(1.0*svref/(onem*dline(i))) rspedn(i,l) = 1.0/speedn(i,l) if (i .ge.i0+ 1-nbdy .and. & i .le.i0+ii+nbdy .and. & j+1.ge.j0+ 1-nbdy .and. & j+1.le.j0+jj+nbdy ) then ivopn(i-i0,j-j0+1) = 1 endif enddo ! elseif (kdport(l).eq.2.or.kdport(l).eq.5) then ! ! southern port ! sum = 0.0 j = jfport(l) i = ifport(l) call xclget(dline(i),lnport(l), depths,i,j,1,0, 0) call xclget(xline(i),lnport(l), scuy, i,j, 1,0, 0) do i= ifport(l),ilport(l) sum = sum + dline(i)*xline(i) enddo sum = 1.e6/sum do i= ifport(l),ilport(l) vports(i,l) = sum speeds(i,l) = sqrt(1.0*svref/(onem*dline(i))) rspeds(i,l) = 1.0/speeds(i,l) if (i.ge.i0+ 1-nbdy .and. & i.le.i0+ii+nbdy .and. & j.ge.j0+ 1-nbdy .and. & j.le.j0+jj+nbdy ) then ivopn(i-i0,j-j0) = 1 endif enddo ! endif ! enddo !nports if (mnproc.eq.1) then write(lp,*) call flush(lp) endif ! ! end of initialization ! call xcsync(flush_lp) return endif !lcount=1 if ((lll.eq.0) .or. (n.eq.0)) return !called from dpthuv ! ! --- 'wellposed' treatment of pressure and normal velocity fields ! --- not in fact wellposed with this exterior data ! ! --- set ramping factor ! when ramp_time=0:no ramping of tide (=full tide) ! when ramp_time>0: ramping of tide, if... ! ...ramp_orig<=timermp<=(ramp_orig+ramp_time) if (ramp_time.gt.0.0) then timermp=d_time if (timermp.ge.ramp_orig) then timermp=(timermp-ramp_orig)/ramp_time ! frmp=(1-exp(-10*timermp)) frmp=(1-exp(-5*timermp)) else frmp=0.0 endif else frmp=1.0 endif ! do l=1,nports if (kdport(l).eq.4) then ! ! western port ! i = ifport(l) j = jfport(l) call xclget(dline(j), lnport(l), & depthu, i,j,0,1, 0) call xclget(pline(j), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i, j,0,1, 0) call xclget(uline(j),lnport(l), & ubavg(1-nbdy,1-nbdy,n), i+1,j,0,1, 0) if (tidflg.ge.1) then !tide do jn= jfport(l),jlport(l) do ic = 1, ncon tmpr(jn,ic) = z1r_u(l,jn,ic) tmpi(jn,ic) = z1i_u(l,jn,ic) enddo enddo call tides_driver(tmpr,tmpi,d_time, & astroflag,uu,j,max(itdm,jtdm),lnport(l)) !normal ! Note!! uu and vv from tpx are transports do jn= jfport(l),jlport(l) upred(jn)=uu(jn)*frmp*onem/dline(jn) enddo do jn= jfport(l),jlport(l) do ic = 1, ncon tmpr(jn,ic) = z1r_p(l,jn,ic) tmpi(jn,ic) = z1i_p(l,jn,ic) enddo enddo call tides_driver(tmpr,tmpi,d_time, & astroflag,zpred,j,max(itdm,jtdm),lnport(l)) do j= jfport(l),jlport(l) zpred(j)=zpred(j)*onem*frmp enddo do j= jfport(l),jlport(l) ! ----set both u and eta at boundary; 1 invariant weighted: uline(j)=ulow(j,l)+upred(j) & +w_1*speedw(j,l)*((plow(j,l)+zpred(j))-pline(j)) pline(j)=w_1c*(plow(j,l)+zpred(j))+w_1*pline(j) enddo else !no bnd-tide do j= jfport(l),jlport(l) ! ---- set both u and eta; 1 invariant weighted: uline(j)=ulow(j,l)+ & w_1*speedw(j,l)*(plow(j,l)-pline(j)) pline(j)=w_1c*plow(j,l)+w_1*pline(j) enddo endif !tide/no tide j = jfport(l) call xclput(pline(j), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i, j,0,1) call xclput(uline(j),lnport(l), & ubavg(1-nbdy,1-nbdy,n), i, j,0,1) ! elseif (kdport(l).eq.3) then ! eastern port ! i = ifport(l)-1 j = jfport(l) call xclget(dline(j), lnport(l), & depthu, i+1,j,0,1, 0) call xclget(pline(j), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i, j,0,1, 0) call xclget(uline(j),lnport(l), & ubavg(1-nbdy,1-nbdy,n), i, j,0,1, 0) if (tidflg.ge.1) then !tide do jn= jfport(l),jlport(l) do ic = 1, ncon tmpr(jn,ic) = z1r_u(l,jn,ic) tmpi(jn,ic) = z1i_u(l,jn,ic) enddo enddo call tides_driver(tmpr,tmpi,d_time, & astroflag,uu,j,max(itdm,jtdm),lnport(l)) !normal ! Note!! uu and vv from tpx are transports do jn= jfport(l),jlport(l) upred(jn)=uu(jn)*frmp*onem/dline(jn) enddo do jn= jfport(l),jlport(l) do ic = 1, ncon tmpr(jn,ic) = z1r_p(l,jn,ic) tmpi(jn,ic) = z1i_p(l,jn,ic) enddo enddo call tides_driver(tmpr,tmpi,d_time, & astroflag,zpred,j,max(itdm,jtdm),lnport(l)) do j= jfport(l),jlport(l) zpred(j) = zpred(j)*onem*frmp enddo do j= jfport(l),jlport(l) ! ----set u and eta on boundary; 1 invariant weighted: uline(j)=ulow(j,l)+upred(j) & +w_1*speede(j,l)*(pline(j)-(plow(j,l)+zpred(j))) pline(j)=w_1c*(plow(j,l)+zpred(j))+w_1*pline(j) enddo else !no bnd-tide do j= jfport(l),jlport(l) ! ----set u and eta on boundary; 1 invariant weighted: uline(j)=ulow(j,l) & -w_1*speede(j,l)*(plow(j,l)-pline(j)) pline(j)=w_1c*plow(j,l)+w_1*pline(j) enddo endif !tide/no tide j = jfport(l) call xclput(pline(j), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i, j,0,1) call xclput(uline(j),lnport(l), & ubavg(1-nbdy,1-nbdy,n), i+1,j,0,1) elseif (kdport(l).eq.1) then ! ! northern port ! j = jfport(l)-1 i = ifport(l) call xclget(dline(i), lnport(l), & depthv, i,j+1,1,0, 0) call xclget(pline(i), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i,j, 1,0, 0) call xclget(uline(i),lnport(l), & vbavg(1-nbdy,1-nbdy,n), i,j, 1,0, 0) if (tidflg.ge.1) then !tide do in= ifport(l),ilport(l) do ic = 1, ncon tmpr(in,ic) = z1r_u(l,in,ic) tmpi(in,ic) = z1i_u(l,in,ic) enddo enddo call tides_driver(tmpr,tmpi,d_time, & astroflag,uu,i,max(itdm,jtdm),lnport(l)) !normal ! Note!! uu and vv from tpx are transports do in= ifport(l),ilport(l) vpred(in)=uu(in)*frmp*onem/dline(in) enddo do in= ifport(l),ilport(l) do ic = 1, ncon tmpr(in,ic) = z1r_p(l,in,ic) tmpi(in,ic) = z1i_p(l,in,ic) enddo enddo call tides_driver(tmpr,tmpi,d_time, & astroflag,zpred,i,max(itdm,jtdm),lnport(l)) do i= ifport(l),ilport(l) zpred(i) =zpred(i)*onem*frmp enddo do i= ifport(l),ilport(l) ! ----set u and eta at boundary; 1 invariant weighted: uline(i)=ulow(i,l)+vpred(i) & +w_1*speedn(i,l)*(pline(i)-(plow(i,l)+zpred(i))) pline(i)=w_1c*(plow(i,l)+zpred(i))+w_1*pline(i) enddo else !no bnd-tide do i= ifport(l),ilport(l) ! ----set u and eta at boundary; 1 invariant weighted: uline(i)=ulow(i,l)- & w_1*speedn(i,l)*(plow(i,l)-pline(i)) pline(i)=w_1c*plow(i,l)+w_1*pline(i) enddo endif !tide/no-tide i = ifport(l) call xclput(pline(i), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i,j, 1,0) call xclput(uline(i),lnport(l), & vbavg(1-nbdy,1-nbdy,n), i,j+1,1,0) ! elseif (kdport(l).eq.2.or.kdport(l).eq.5) then ! ! southern port ! j = jfport(l) i = ifport(l) call xclget(dline(i), lnport(l), & depthv, i,j, 1,0, 0) call xclget(pline(i), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i,j, 1,0, 0) call xclget(uline(i),lnport(l), & vbavg(1-nbdy,1-nbdy,n), i,j+1,1,0, 0) if (tidflg.ge.1) then !tide do in= ifport(l),ilport(l) do ic = 1, ncon tmpr(in,ic) = z1r_u(l,in,ic) tmpi(in,ic) = z1i_u(l,in,ic) enddo enddo call tides_driver(tmpr,tmpi,d_time, & astroflag,uu,i,max(itdm,jtdm),lnport(l)) !normal ! Note!! uu and vv from tpx are transports do in= ifport(l),ilport(l) vpred(in)=uu(in)*frmp*onem/dline(in) enddo do in= ifport(l),ilport(l) do ic = 1, ncon tmpr(in,ic) = z1r_p(l,in,ic) tmpi(in,ic) = z1i_p(l,in,ic) enddo enddo call tides_driver(tmpr,tmpi,d_time, & astroflag,zpred,i,max(itdm,jtdm),lnport(l)) do i= ifport(l),ilport(l) zpred(i) = zpred(i)*onem*frmp enddo do i= ifport(l),ilport(l) ! ----set u and eta at boundary; 1 invariant weighted: uline(i)=ulow(i,l)+vpred(i) & +w_1*speeds(i,l)*(pline(i)-(plow(i,l)+zpred(i))) pline(i)=w_1c*(plow(i,l)+zpred(i))+w_1*pline(i) enddo else !no bnd-tide do i= ifport(l),ilport(l) ! ----set u and eta at boundary; 1 invariant weighted: uline(i)=ulow(i,l)+ & w_1*speeds(i,l)*(plow(i,l)-pline(i)) pline(i)=w_1c*plow(i,l)+w_1*pline(i) enddo endif !tide/no-tide i = ifport(l) call xclput(pline(i), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i,j, 1,0) call xclput(uline(i),lnport(l), & vbavg(1-nbdy,1-nbdy,n), i,j, 1,0) ! endif !kdport enddo !nports ! return end subroutine latbdf subroutine latbdp(n) use mod_xc ! HYCOM communication interface use mod_cb_arrays ! HYCOM saved arrays implicit none ! integer n ! ! --- apply lateral boundary conditions to barotropic flow field ! ! --- port flow version: ! --- similar to the standard 'Browning and Kreiss' MICOM/HYCOM open ! --- boundary condition, except that the exterior normal velocity ! --- is constant in time and exterior pressure = interior pressure. ! --- tangential velocity is not constrained. ! ! --- see also: latbdp ! ! --- the code is as similar as possible to that for the standard case. ! --- so for example, 'speed' is in fact 1/SQRT(gH) which represents ! --- c1/g in the notation of (Bleck and Sun, Open boundary conditions ! --- for MICOM). The 1/g allows for the use of pressure fields. ! ! --- Note that East, West, North and South refers to the grid ! --- (i.e i,j points) and NOT geographic East, West, North and South ! ! --- the first call is made during initialization. ! ! --- Alan J. Wallcraft, NRL, November 1999. ! logical, parameter :: ldebug_latbdp=.false. ! integer, parameter :: mports=1 !maximum number of ports ! integer, parameter :: nchar=120 ! logical lfatal,lfatalp integer i,j,isec,ifrst,ilast,l real aline(nchar), & dline(itdm+jtdm),xline(itdm+jtdm), & pline(itdm+jtdm),uline(itdm+jtdm,2) real crs,fin,sum,svspin,uvscl,uvscl2,fatal real*8 tstep character*3 char3 ! integer nports,kdport(mports), & ifport(mports),ilport(mports), & jfport(mports),jlport(mports),lnport(mports) real pefold,svpnow(mports),svport(mports) real*8 refold save nports,kdport,ifport,ilport,jfport,jlport,lnport save pefold,svpnow,svport,refold ! real uportw(jtdm),speedw(jtdm),rspedw(jtdm), & uporte(jtdm),speede(jtdm),rspede(jtdm), & vportn(itdm),speedn(itdm),rspedn(itdm), & vports(itdm),speeds(itdm),rspeds(itdm) save uportw,speedw,rspedw,uporte,speede,rspede, & vportn,speedn,rspedn,vports,speeds,rspeds ! character*13 fmt save fmt data fmt / '(i4,1x,120i1)' / ! integer lcount save lcount data lcount / 0 / ! lcount = lcount + 1 ! ! --- the first call just initializes data structures. ! if (lcount.eq.1) then ! open(unit=uoff+99,file=trim(flnminp)//'ports.input') ! ! --- 'nports' = number of boundary port sections. call blkini(nports,'nports') if (mnproc.eq.1) then write(lp,*) endif if (nports.eq.0) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdp - must have lbflag=0 for nports=0' write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbdp)') stop '(latbdp)' elseif (nports.lt.0 .or. nports.gt.mports) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdp - illegal nports value' if (nports.gt.mports) then write(lp,*) 'increase parameter mports to',nports endif write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif ! ! --- 'pefold' = port transport e-folding time in days call blkinr(pefold,'pefold','(a6," =",f10.4," days")') if (mnproc.eq.1) then write(lp,*) endif ! ! --- switch units from days to baroclinic time steps ! --- shift lcount to prevent underflow (lcount*refold.ge.0.001) ! tstep = pefold*(86400.d0/batrop) refold = 1.d0/tstep lcount = lcount + int(tstep)/1000 ! ! --- read in the ports one at a time ! do l= 1,nports ! ! --- port location is w.r.t. u (EW) or v (NS) grid ! --- and identifies the sea at the port ! --- the minimum index is 0 ! ! --- 'kdport' = port orientation (1=N, 2=S, 3=E, 4=W) ! --- 'ifport' = first i-index ! --- 'ilport' = last i-index (=ifport for N or S orientation) ! --- 'jfport' = first j-index ! --- 'jlport' = last j-index (=jfport for E or W orientation) ! --- 'svpnow' = existing port transport in Sv (+ve towards E or S) ! --- 'svport' = target port transport in Sv (+ve towards E or S) ! --- 'lnport' = port length (calculated, not input) call blkini(kdport(l),'kdport') call blkini(ifport(l),'ifport') call blkini(ilport(l),'ilport') call blkini(jfport(l),'jfport') call blkini(jlport(l),'jlport') call blkinr(svpnow(l),'svpnow','(a6," =",f10.4," Sv")') call blkinr(svport(l),'svport','(a6," =",f10.4," Sv")') if (mnproc.eq.1) then write(lp,*) endif ! lnport(l) = ilport(l)-ifport(l)+jlport(l)-jfport(l)+1 ! ! --- sanity check. ! if (kdport(l).gt.2) then if (ifport(l).ne.ilport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdp - port direction', & ' and orientation are not consistent' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif else if (jfport(l).ne.jlport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdp - port direction', & ' and orientation are not consistent' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif endif if (ifport(l).gt.ilport(l) .or. & jfport(l).gt.jlport(l) ) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdp - port', & ' location is not consistent' write(lp,*) call flush(lp) endif call xcstop('(latbdp)') stop '(latbdp)' endif enddo ! close(unit=uoff+99) ! ! --- check ports against masks, ! --- mark the port locations on masks and print them out. ! lfatal = .false. do l= 1,nports lfatalp = .false. ! if (kdport(l).eq.4) then ! ! western port ! i = ifport(l) do j= jfport(l),jlport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.1 .or. j.gt.jtdm ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iu(i-i0,j-j0).ne.0) then lfatalp = .true. iu(i-i0,j-j0) = 9 !indicate an error else iu(i-i0,j-j0) = -1 endif if (iu(i-i0+1,j-j0).ne.1 .or. & iu(i-i0+2,j-j0).ne.1 ) then lfatalp = .true. iu(i-i0,j-j0) = 7 !indicate an error endif enddo ! elseif (kdport(l).eq.3) then ! ! eastern port ! i = ifport(l) do j= jfport(l),jlport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.1 .or. j.gt.jtdm ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iu(i-i0,j-j0).ne.0) then lfatalp = .true. iu(i-i0,j-j0) = 9 !indicate an error else iu(i-i0,j-j0) = -1 endif if (iu(i-i0-1,j-j0).ne.1 .or. & iu(i-i0-2,j-j0).ne.1 ) then lfatalp = .true. iu(i-i0,j-j0) = 7 !indicate an error endif enddo ! elseif (kdport(l).eq.1) then ! ! northern port ! j = jfport(l) do i= ifport(l),ilport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.3 .or. j.gt.jtdm ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iv(i-i0,j-j0).ne.0) then lfatalp = .true. iv(i-i0,j-j0) = 9 !indicate an error else iv(i-i0,j-j0) = -1 endif if (iv(i-i0,j-j0-1).ne.1 .or. & iv(i-i0,j-j0-2).ne.1 ) then lfatalp = .true. iv(i-i0,j-j0) = 7 !indicate an error endif enddo ! elseif (kdport(l).eq.2) then ! ! southern port ! j = jfport(l) do i= ifport(l),ilport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.1 .or. j.gt.jtdm-2 ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iv(i-i0,j-j0).ne.0) then lfatalp = .true. iv(i-i0,j-j0) = 9 !indicate an error else iv(i-i0,j-j0) = -1 endif if (iv(i-i0,j-j0+1).ne.1 .or. & iv(i-i0,j-j0+2).ne.1 ) then lfatalp = .true. iv(i-i0,j-j0) = 7 !indicate an error endif enddo ! endif ! if (lfatalp) then write(lp,*) write(lp,*) 'error in latbdp - port ',l,' mislocated', & ' (mnproc = ',mnproc,')' write(lp,*) call flush(lp) endif lfatal = lfatal .or. lfatalp enddo ! ! local lfatal to global lfatal ! if (lfatal) then fatal = 1.0 else fatal = 0.0 endif call xcmaxr(fatal) lfatal = fatal.gt.0.5 ! ! --- write out -iu- and -iv- arrays, if they are not too big ! --- data are written in strips nchar points wide if (lfatal .or. max(itdm,jtdm).le.2*nchar) then util1(1:ii,1:jj) = iu(1:ii,1:jj) ! xclget is for real arrays isec=(itdm-1)/nchar do ifrst=0,nchar*isec,nchar ilast=min(itdm,ifrst+nchar) write (char3,'(i3)') ilast-ifrst fmt(8:10)=char3 if (mnproc.eq.1) then write (lp,'(a,i5,a,i5)') & 'iu array, cols',ifrst+1,' --',ilast endif do j= jtdm,1,-1 call xclget(aline,ilast-ifrst, util1,ifrst+1,j,1,0, 1) if (mnproc.eq.1) then write (lp,fmt) j,(nint(aline(i)),i=1,ilast-ifrst) endif enddo enddo if (mnproc.eq.1) then write (lp,*) endif call xcsync(flush_lp) ! util1(1:ii,1:jj) = iv(1:ii,1:jj) ! xclget is for real arrays isec=(itdm-1)/nchar do ifrst=0,nchar*isec,nchar ilast=min(itdm,ifrst+nchar) write (char3,'(i3)') ilast-ifrst fmt(8:10)=char3 if (mnproc.eq.1) then write (lp,'(a,i5,a,i5)') & 'iv array, cols',ifrst+1,' --',ilast endif do j= jtdm,1,-1 call xclget(aline,ilast-ifrst, util1,ifrst+1,j,1,0, 1) if (mnproc.eq.1) then write (lp,fmt) j,(nint(aline(i)),i=1,ilast-ifrst) endif enddo enddo if (mnproc.eq.1) then write (lp,*) endif call xcsync(flush_lp) endif ! small region ! if (lfatal) then write(lp,*) write(lp,*) 'error in latbdp - bad port(s)' write(lp,*) call flush(lp) call xchalt('(latbdp)') stop '(latbdp)' endif ! ! --- restore iu and iv, and zero iuopn and ivopn. ! !$OMP PARALLEL DO PRIVATE(j,i) & !$OMP SCHEDULE(STATIC,jblk) do j= 1,jj do i= 1,ii iu(i,j) = max( iu(i,j), 0 ) iv(i,j) = max( iv(i,j), 0 ) enddo enddo !$OMP PARALLEL DO PRIVATE(j,i) & !$OMP SCHEDULE(STATIC,jblk) do j= 1-nbdy,jj+nbdy do i= 1-nbdy,ii+nbdy iuopn(i,j) = 0 ivopn(i,j) = 0 enddo enddo ! ! --- initialize the ports ! do l= 1,nports if (kdport(l).eq.4) then ! ! western port ! sum = 0.0 i = ifport(l) j = jfport(l) call xclget(dline(j),lnport(l), depths,i+1,j,0,1, 0) call xclget(xline(j),lnport(l), scuy, i, j,0,1, 0) do j= jfport(l),jlport(l) sum = sum + dline(j)*xline(j) enddo sum = 1.e6/sum do j= jfport(l),jlport(l) uportw(j) = sum speedw(j) = sqrt(svref/(onem*dline(j))) rspedw(j) = 1.0/speedw(j) if (mnproc.eq.1) then write(lp,'(a,i2,2i5,1p2e13.5)') & 'w port: ',l,i,j,uportw(j),speedw(j) endif ! if (i.ge.i0+ 1-nbdy .and. & i.le.i0+ii+nbdy .and. & j.ge.j0+ 1-nbdy .and. & j.le.j0+jj+nbdy ) then iuopn(i-i0,j-j0) = 1 endif enddo ! elseif (kdport(l).eq.3) then ! ! eastern port ! sum = 0.0 i = ifport(l)-1 j = jfport(l) call xclget(dline(j),lnport(l), depths,i, j,0,1, 0) call xclget(xline(j),lnport(l), scuy, i+1,j,0,1, 0) do j= jfport(l),jlport(l) sum = sum + dline(j)*xline(j) enddo sum = 1.e6/sum do j= jfport(l),jlport(l) uporte(j) = sum speede(j) = sqrt(svref/(onem*dline(j))) rspede(j) = 1.0/speede(j) if (mnproc.eq.1) then write(lp,'(a,i2,2i5,1p2e13.5)') & 'e port: ',l,i,j,uporte(j),speede(j) endif ! if (i+1.ge.i0+ 1-nbdy .and. & i+1.le.i0+ii+nbdy .and. & j .ge.j0+ 1-nbdy .and. & j .le.j0+jj+nbdy ) then iuopn(i-i0+1,j-j0) = 1 endif enddo ! elseif (kdport(l).eq.1) then ! ! northern port ! sum = 0.0 j = jfport(l)-1 i = ifport(l) call xclget(dline(i),lnport(l), depths,i,j, 1,0, 0) call xclget(xline(i),lnport(l), scuy, i,j+1,1,0, 0) do i= ifport(l),ilport(l) sum = sum + dline(i)*xline(i) enddo sum = 1.e6/sum do i= ifport(l),ilport(l) vportn(i) = sum speedn(i) = sqrt(svref/(onem*dline(i))) rspedn(i) = 1.0/speedn(i) if (mnproc.eq.1) then write(lp,'(a,i2,2i5,1p2e13.5)') & 'n port: ',l,i,j,vportn(i),speedn(i) endif ! if (i .ge.i0+ 1-nbdy .and. & i .le.i0+ii+nbdy .and. & j+1.ge.j0+ 1-nbdy .and. & j+1.le.j0+jj+nbdy ) then ivopn(i-i0,j-j0+1) = 1 endif enddo ! elseif (kdport(l).eq.2) then ! ! southern port ! sum = 0.0 j = jfport(l) i = ifport(l) call xclget(dline(i),lnport(l), depths,i,j+1,1,0, 0) call xclget(xline(i),lnport(l), scuy, i,j, 1,0, 0) do i= ifport(l),ilport(l) sum = sum + dline(i)*xline(i) enddo sum = 1.e6/sum do i= ifport(l),ilport(l) vports(i) = sum speeds(i) = sqrt(svref/(onem*dline(i))) rspeds(i) = 1.0/speeds(i) if (mnproc.eq.1) then write(lp,'(a,i2,2i5,1p2e13.5)') & 's port: ',l,i,j,vports(i),speeds(i) endif ! if (i.ge.i0+ 1-nbdy .and. & i.le.i0+ii+nbdy .and. & j.ge.j0+ 1-nbdy .and. & j.le.j0+jj+nbdy ) then ivopn(i-i0,j-j0) = 1 endif enddo ! endif ! if (mnproc.eq.1) then write(lp,*) 'port, now/target velocity = ', & l,svpnow(l)*sum,svport(l)*sum call flush(lp) endif enddo if (mnproc.eq.1) then write(lp,*) call flush(lp) endif ! ! end of initialization ! call xcsync(flush_lp) return endif if (n.eq.0) return !called from dpthuv ! ! --- 'wellposed' treatment of pressure and normal velocity fields ! --- not in fact wellposed with this exterior data ! tstep = lcount svspin = exp( -tstep*refold ) do l= 1,nports uvscl = svport(l) + svspin*(svpnow(l)-svport(l)) ! if (kdport(l).eq.4) then ! ! western port ! i = ifport(l) j = jfport(l) call xclget(dline(j), lnport(l), & depthu, i+1,j,0,1, 0) call xclget(xline(j), lnport(l), & scuy, i, j,0,1, 0) call xclget(pline(j), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i, j,0,1, 0) call xclget(uline(j,1),lnport(l), & ubavg(1-nbdy,1-nbdy,n), i+1,j,0,1, 0) call xclget(uline(j,2),lnport(l), & ubavg(1-nbdy,1-nbdy,n), i+2,j,0,1, 0) sum = 0.0 do j= jfport(l),jlport(l) crs=uvscl*uportw(j)+speedw(j)*pline(j) fin=1.5*uline(j,1)-.5*uline(j,2)-speedw(j)*pline(j) sum=sum+((crs+fin)-uline(j,1))*dline(j)*xline(j) enddo uvscl2 = uvscl + (uvscl - sum/(onem*1.e6)) sum = 0.0 do j= jfport(l),jlport(l) crs=uvscl2*uportw(j)+speedw(j)*pline(j) fin=1.5*uline(j,1)-.5*uline(j,2)-speedw(j)*pline(j) pline(j) =.5*(crs-fin)*rspedw(j) uline(j,1)=(crs+fin)-uline(j,1) sum=sum+uline(j,1)*dline(j)*xline(j) enddo j = jfport(l) call xclput(pline(j), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i, j,0,1) call xclput(uline(j,1),lnport(l), & ubavg(1-nbdy,1-nbdy,n), i, j,0,1) ! if (ldebug_latbdp .and. mnproc.eq.1) then write(lp,'(a,i2,3i5,2i2)') 'l,xclput(pb - ', & l,lnport(l),i, j,0,1 write(lp,'(a,i2,3i5,2i2)') 'l,xclput(ub - ', & l,lnport(l),i, j,0,1 call flush(lp) endif ! elseif (kdport(l).eq.3) then ! ! eastern port ! i = ifport(l)-1 j = jfport(l) call xclget(dline(j), lnport(l), & depthu, i+1,j,0,1, 0) call xclget(xline(j), lnport(l), & scuy, i+1,j,0,1, 0) call xclget(pline(j), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i, j,0,1, 0) call xclget(uline(j,1),lnport(l), & ubavg(1-nbdy,1-nbdy,n), i, j,0,1, 0) call xclget(uline(j,2),lnport(l), & ubavg(1-nbdy,1-nbdy,n), i-1,j,0,1, 0) sum = 0.0 do j= jfport(l),jlport(l) crs=uvscl*uporte(j)-speede(j)*pline(j) fin=1.5*uline(j,1)-.5*uline(j,2)+speede(j)*pline(j) sum=sum+((crs+fin)-uline(j,1))*dline(j)*xline(j) enddo uvscl2 = uvscl + (uvscl - sum/(onem*1.e6)) sum = 0.0 do j= jfport(l),jlport(l) crs=uvscl2*uporte(j)-speede(j)*pline(j) fin=1.5*uline(j,1)-.5*uline(j,2)+speede(j)*pline(j) pline(j) =.5*(fin-crs)*rspede(j) uline(j,1)=(fin+crs)-uline(j,1) sum=sum+uline(j,1)*dline(j)*xline(j) ! if (mnproc.eq.1) then ! write(lp,'(a,i2,2i5,1p2e13.5)') ! & 'e port: ',l,i,j,pline(j),uline(j,1) ! endif enddo j = jfport(l) call xclput(pline(j), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i, j,0,1) call xclput(uline(j,1),lnport(l), & ubavg(1-nbdy,1-nbdy,n), i+1,j,0,1) ! if (ldebug_latbdp .and. mnproc.eq.1) then write(lp,'(a,i2,3i5,2i2)') 'l,xclput(pb - ', & l,lnport(l),i, j,0,1 write(lp,'(a,i2,3i5,2i2)') 'l,xclput(ub - ', & l,lnport(l),i+1,j,0,1 call flush(lp) endif ! elseif (kdport(l).eq.1) then ! ! northern port ! j = jfport(l)-1 i = ifport(l) call xclget(dline(i), lnport(l), & depthv, i,j+1,1,0, 0) call xclget(xline(i), lnport(l), & scux, i,j+1,1,0, 0) call xclget(pline(i), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i,j, 1,0, 0) call xclget(uline(i,1),lnport(l), & vbavg(1-nbdy,1-nbdy,n), i,j, 1,0, 0) call xclget(uline(i,2),lnport(l), & vbavg(1-nbdy,1-nbdy,n), i,j-1,1,0, 0) sum = 0.0 do i= ifport(l),ilport(l) crs=uvscl*vportn(i)-speedn(i)*pline(i) fin=1.5*uline(i,1)-.5*uline(i,2)+speedn(i)*pline(i) sum=sum+((fin+crs)-uline(i,1))*dline(i)*xline(i) enddo uvscl2 = uvscl + (uvscl - sum/(onem*1.e6)) sum = 0.0 do i= ifport(l),ilport(l) crs=uvscl2*vportn(i)-speedn(i)*pline(i) fin=1.5*uline(i,1)-.5*uline(i,2)+speedn(i)*pline(i) pline(i) =.5*(fin-crs)*rspedn(i) uline(i,1)=(fin+crs)-uline(i,1) sum=sum+uline(i,1)*dline(i)*xline(i) enddo i = ifport(l) call xclput(pline(i), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i,j, 1,0) call xclput(uline(i,1),lnport(l), & vbavg(1-nbdy,1-nbdy,n), i,j+1,1,0) ! if (ldebug_latbdp .and. mnproc.eq.1) then write(lp,'(a,i2,3i5,2i2)') 'l,xclput(pb - ', & l,lnport(l),i,j, 1,0 write(lp,'(a,i2,3i5,2i2)') 'l,xclput(vb - ', & l,lnport(l),i,j+1,1,0 call flush(lp) endif ! elseif (kdport(l).eq.2) then ! ! southern port ! j = jfport(l) i = ifport(l) call xclget(dline(i), lnport(l), & depthv, i,j, 1,0, 0) call xclget(xline(i), lnport(l), & scux, i,j, 1,0, 0) call xclget(pline(i), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i,j, 1,0, 0) call xclget(uline(i,1),lnport(l), & vbavg(1-nbdy,1-nbdy,n), i,j+1,1,0, 0) call xclget(uline(i,2),lnport(l), & vbavg(1-nbdy,1-nbdy,n), i,j+2,1,0, 0) sum = 0.0 do i= ifport(l),ilport(l) crs=uvscl*vports(i)+speeds(i)*pline(i) fin=1.5*uline(i,1)-.5*uline(i,2)-speeds(i)*pline(i) sum=sum+((crs+fin)-uline(i,1))*dline(i)*xline(i) enddo uvscl2 = uvscl + (uvscl - sum/(onem*1.e6)) sum = 0.0 do i= ifport(l),ilport(l) crs=uvscl2*vports(i)+speeds(i)*pline(i) fin=1.5*uline(i,1)-.5*uline(i,2)-speeds(i)*pline(i) pline(i) =.5*(crs-fin)*rspeds(i) uline(i,1)=(crs+fin)-uline(i,1) sum=sum+uline(i,1)*dline(i)*xline(i) enddo i = ifport(l) call xclput(pline(i), lnport(l), & pbavg(1-nbdy,1-nbdy,n), i,j, 1,0) call xclput(uline(i,1),lnport(l), & vbavg(1-nbdy,1-nbdy,n), i,j, 1,0) ! if (ldebug_latbdp .and. mnproc.eq.1) then write(lp,'(a,i2,3i5,2i2)') 'l,xclput(pb - ', & l,lnport(l),i,j, 1,0 write(lp,'(a,i2,3i5,2i2)') 'l,xclput(vb - ', & l,lnport(l),i,j, 1,0 call flush(lp) endif ! endif ! ! if (mod(lcount,512).eq.0) then ! if (mnproc.eq.1) then ! write(lp,*) 'latbdp - l,sv,sum = ',l,uvscl,sum/(onem*1.e6) ! call flush(lp) ! endif ! endif enddo ! return end subroutine latbdp #endif /* ! RELO */ subroutine latbdt(n,lll) use mod_xc ! HYCOM communication interface use mod_cb_arrays ! HYCOM saved arrays use mod_tides ! HYCOM Tidal Data implicit none ! integer n,lll ! ! --- apply lateral boundary conditions to barotropic flow field ! ! --- nested sub-region version: ! --- Uses the 'Browning and Kreiss' MICOM open boundary condition. ! ! --- Tidal z,u & v components added to height and velocity fields ! --- at open boundary points when tidflg = 1 or 3. ! ! --- Note that 'speed' is in fact qonem*SQRT(g/H). ! --- The qonem allows for the use of pressure fields. ! ! --- Note that East, West, North and South refers to the grid ! --- (i.e i,j points) and NOT geographic East, West, North and South ! ! --- the first call is made during initialization. ! ! --- Alan J. Wallcraft, NRL, July, 2001. ! --- Dan Moore, QNA, August 2011 (tidal components). ! logical, parameter :: ldebug_latbdt=.false. !usually .false. logical, parameter :: ltrans_latbdt=.false. !usually .false. ! integer, parameter :: mnp=1 !mnflg for xciegt and xciput !0 == all tasks; 1 == task 1; #if defined(LATBDT_NPLINE3) integer, parameter :: npline=3 !set by a CPP macro #else integer, parameter :: npline=1 !usually 1 #endif !update pline every npline time steps !npline is odd, can be 1 integer, parameter :: nchar=120 ! logical lfatal,lfatalp integer i,im1,ip1,j,jm1,jp1,k,isec,ifrst,ilast,l,np real aline(nchar) real fatal,sb0,sb1 character*3 char3 ! integer, save :: & nportpts, & nptest, & nports !number of ports real, save :: & sports !scale factor for nested velocity integer, save, allocatable :: & nxport(:), & kdport(:), & ifport(:),ilport(:), & jfport(:),jlport(:),lnport(:) real*8, save, allocatable :: & trport(:), & trporp(:), & trnest(:) ! integer jline_start,jline,jtide_start,jtide,tidcon1 ! integer, save, allocatable :: & iaub(:),iavb(:),jaub(:),javb(:), & !boundary vel-point iaui(:),iavi(:),jaui(:),javi(:), & !1st sea vel-point inside boundary iau2(:),iav2(:),jau2(:),jav2(:), & !2nd sea vel-point inside boundary iapi(:), japi(:), & !1st sea p-point ndup(:) !1st sea p-point, duplicates real, save, allocatable :: pspeed(:),rspeed(:), & utrans(:),vtrans(:), & !velocity to transport pline(:),uline(:),vline(:), & ulin2(:),vlin2(:), & plnst(:),ulnst(:),vlnst(:), & crs(:), fin(:) ! real, save, allocatable :: z_R(:,:,:),z_I(:,:,:),z_A(:) real, save, allocatable :: port_tide(:,:) real*8, save :: d_time ! character*13 fmt save fmt data fmt / '(i4,1x,120i1)' / ! integer lcount save lcount data lcount / 0 / ! lcount = lcount + 1 ! ! --- the first call just initializes data structures. ! if (lcount.eq.1) then !call at start of dpthuv if (mnproc.eq.1) then write(lp,'(/a,i3/a)') & 'Browning and Kreiss barotropic nesting, npline =',npline, & ' now processing ports.input ...' endif !1st tile call xcsync(flush_lp) ! open(unit=uoff+99,file=trim(flnminp)//'ports.input') ! ! --- 'sports' = scale factor for nested velocity (optional, default 1.0) ! --- 'nports' = number of boundary port sections. call blkinr2(sports,i, & 'sports','(a6," =",f10.4," ")', & 'nports','(a6," =",f11.0)') !treat nports as real if (i.eq.1) then !sports call blkini(nports,'nports') else !nports nports = sports sports = 1.0 !default endif if (mnproc.eq.1) then write(lp,*) endif if (nports.eq.0) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdt - must have lbflag=0 for nports=0' write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbdt)') stop '(latbdt)' elseif (nports.lt.0) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdt - illegal nports value' write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbdt)') stop '(latbdt)' endif ! allocate( & nxport(nports), & kdport(nports), & ifport(nports), & ilport(nports), & jfport(nports), & jlport(nports), & lnport(nports) ) allocate( & trport(nports), & trporp(nports), & trnest(nports) ) ! ! --- read in the ports one at a time ! nportpts = 0 do l= 1,nports ! ! --- port location is w.r.t. u (EW) or v (NS) grid ! --- and identifies the sea at the port ! --- the minimum index is 0 ! ! --- 'kdport' = port orientation (1=N, 2=S, 3=E, 4=W) ! --- 'ifport' = first i-index ! --- 'ilport' = last i-index (=ifport for N or S orientation) ! --- 'jfport' = first j-index ! --- 'jlport' = last j-index (=jfport for E or W orientation) ! --- 'lnport' = port length (calculated, not input) call blkini(kdport(l),'kdport') call blkini(ifport(l),'ifport') call blkini(ilport(l),'ilport') call blkini(jfport(l),'jfport') call blkini(jlport(l),'jlport') if (mnproc.eq.1) then write(lp,*) endif ! nxport(l) = nportpts+1 lnport(l) = ilport(l)-ifport(l)+ & jlport(l)-jfport(l)+1 nportpts = nportpts+lnport(l) ! ! --- sanity check. ! if (kdport(l).gt.2) then if (ifport(l).ne.ilport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdt - port direction', & ' and orientation are not consistent' write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbdt)') stop '(latbdt)' endif else if (jfport(l).ne.jlport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdt - port direction', & ' and orientation are not consistent' write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbdt)') stop '(latbdt)' endif endif if (ifport(l).gt.ilport(l) .or. & jfport(l).gt.jlport(l) ) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdt - port', & ' location is not consistent' write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbdt)') stop '(latbdt)' endif enddo ! close(unit=uoff+99) ! ! --- check ports against masks, ! --- mark the port locations on masks and print them out. ! lfatal = .false. do l= 1,nports lfatalp = .false. ! if (kdport(l).eq.4) then ! ! western port ! i = ifport(l) do j= jfport(l),jlport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.1 .or. j.gt.jtdm ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iu(i-i0,j-j0).ne.0) then lfatalp = .true. iu(i-i0,j-j0) = 9 !indicate an error else iu(i-i0,j-j0) = -1 endif if (iu(i-i0+1,j-j0).ne.1 .or. & iu(i-i0+2,j-j0).ne.1 ) then lfatalp = .true. iu(i-i0,j-j0) = 7 !indicate an error endif enddo ! elseif (kdport(l).eq.3) then ! ! eastern port ! i = ifport(l) do j= jfport(l),jlport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.1 .or. j.gt.jtdm ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iu(i-i0,j-j0).ne.0) then lfatalp = .true. iu(i-i0,j-j0) = 9 !indicate an error else iu(i-i0,j-j0) = -1 endif if (iu(i-i0-1,j-j0).ne.1 .or. & iu(i-i0-2,j-j0).ne.1 ) then lfatalp = .true. iu(i-i0,j-j0) = 7 !indicate an error endif enddo ! elseif (kdport(l).eq.1) then ! ! northern port ! j = jfport(l) do i= ifport(l),ilport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.3 .or. j.gt.jtdm ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iv(i-i0,j-j0).ne.0) then lfatalp = .true. iv(i-i0,j-j0) = 9 !indicate an error else iv(i-i0,j-j0) = -1 endif if (iv(i-i0,j-j0-1).ne.1 .or. & iv(i-i0,j-j0-2).ne.1 ) then lfatalp = .true. iv(i-i0,j-j0) = 7 !indicate an error endif enddo ! elseif (kdport(l).eq.2) then ! ! southern port ! j = jfport(l) do i= ifport(l),ilport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.1 .or. j.gt.jtdm-2 ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iv(i-i0,j-j0).ne.0) then lfatalp = .true. iv(i-i0,j-j0) = 9 !indicate an error else iv(i-i0,j-j0) = -1 endif if (iv(i-i0,j-j0+1).ne.1 .or. & iv(i-i0,j-j0+2).ne.1 ) then lfatalp = .true. iv(i-i0,j-j0) = 7 !indicate an error endif enddo ! endif ! if (lfatalp) then write(lp,*) write(lp,*) 'error in latbdt - port ',l,' mislocated', & ' (mnproc = ',mnproc,')' write(lp,*) call flush(lp) endif lfatal = lfatal .or. lfatalp enddo !l=1,nports ! ! local lfatal to global lfatal ! if (lfatal) then fatal = 1.0 else fatal = 0.0 endif call xcmaxr(fatal) lfatal = fatal.gt.0.5 ! ! --- write out -iu- and -iv- arrays, if they are not too big ! --- data are written in strips nchar points wide ! if (lfatal .or. max(itdm,jtdm).le.2*nchar) then util1(1:ii,1:jj) = iu(1:ii,1:jj) ! xclget is for real arrays isec=(itdm-1)/nchar do ifrst=0,nchar*isec,nchar ilast=min(itdm,ifrst+nchar) write (char3,'(i3)') ilast-ifrst fmt(8:10)=char3 if (mnproc.eq.1) then write (lp,'(a,i5,a,i5)') & 'iu array, cols',ifrst+1,' --',ilast endif !1st tile do j= jtdm,1,-1 call xclget(aline,ilast-ifrst, util1,ifrst+1,j,1,0, 1) if (mnproc.eq.1) then write (lp,fmt) j,(nint(aline(i)),i=1,ilast-ifrst) endif !1st tile enddo enddo if (mnproc.eq.1) then write (lp,*) endif call xcsync(flush_lp) ! util1(1:ii,1:jj) = iv(1:ii,1:jj) ! xclget is for real arrays isec=(itdm-1)/nchar do ifrst=0,nchar*isec,nchar ilast=min(itdm,ifrst+nchar) write (char3,'(i3)') ilast-ifrst fmt(8:10)=char3 if (mnproc.eq.1) then write (lp,'(a,i5,a,i5)') & 'iv array, cols',ifrst+1,' --',ilast endif !1st tile do j= jtdm,1,-1 call xclget(aline,ilast-ifrst, util1,ifrst+1,j,1,0, 1) if (mnproc.eq.1) then write (lp,fmt) j,(nint(aline(i)),i=1,ilast-ifrst) endif !1st tile enddo enddo if (mnproc.eq.1) then write (lp,*) endif call xcsync(flush_lp) endif ! small region ! if (lfatal) then write(lp,*) write(lp,*) 'error in latbdt - bad port(s)' write(lp,*) call flush(lp) call xchalt('(latbdt)') stop '(latbdt)' endif ! ! --- restore iu and iv, and zero iuopn and ivopn. ! !$OMP PARALLEL DO PRIVATE(j,i) & !$OMP SCHEDULE(STATIC,jblk) do j= 1,jj do i= 1,ii iu(i,j) = max( iu(i,j), 0 ) iv(i,j) = max( iv(i,j), 0 ) enddo enddo !$OMP PARALLEL DO PRIVATE(j,i) & !$OMP SCHEDULE(STATIC,jblk) do j= 1-nbdy,jj+nbdy do i= 1-nbdy,ii+nbdy iuopn(i,j) = 0 ivopn(i,j) = 0 enddo enddo ! ! --- initialize the ports ! allocate(iaub(nportpts),jaub(nportpts), & iavb(nportpts),javb(nportpts), & iaui(nportpts),jaui(nportpts), & iavi(nportpts),javi(nportpts), & iau2(nportpts),jau2(nportpts), & iav2(nportpts),jav2(nportpts), & iapi(nportpts),japi(nportpts), & ndup(nportpts) ) ! do l= 1,nports if (kdport(l).eq.4) then !western port np = nxport(l)-jfport(l) i = ifport(l) if (l.eq.1) then nptest = np+jlport(l) endif do j= jfport(l),jlport(l) iapi(np+j) = i japi(np+j) = j ! --- u is normal to the port iaub(np+j) = i jaub(np+j) = j iaui(np+j) = i+1 jaui(np+j) = j iau2(np+j) = i+2 jau2(np+j) = j ! --- v is tangential to the port iavb(np+j) = i !j=jfport not active javb(np+j) = j !j=jfport not active iavi(np+j) = i javi(np+j) = j iav2(np+j) = i !not active jav2(np+j) = j !not active if (l.eq.1 .and. j.eq.jttest) then nptest = np+j endif if (ldebug_latbdt .and. np+j.eq.nptest & .and. mnproc.eq.1) then write(lp,'(a,5i5)') 'n,ia,ja w:',np+j, & iaub(np+j),jaub(np+j),iaui(np+j),jaui(np+j) endif !ldebug_latbdt enddo !j elseif (kdport(l).eq.3) then !eastern port np = nxport(l)-jfport(l) i = ifport(l)-1 if (l.eq.1) then nptest = np+jlport(l) endif do j= jfport(l),jlport(l) iapi(np+j) = i japi(np+j) = j ! --- u is normal to the port iaub(np+j) = i+1 jaub(np+j) = j iaui(np+j) = i jaui(np+j) = j iau2(np+j) = i-1 jau2(np+j) = j ! --- v is tangential to the port iavb(np+j) = i !j=jfport not active javb(np+j) = j !j=jfport not active iavi(np+j) = i javi(np+j) = j iav2(np+j) = i !not active jav2(np+j) = j !not active if (l.eq.1 .and. j.eq.jttest) then nptest = np+j endif if (ldebug_latbdt .and. np+j.eq.nptest & .and. mnproc.eq.1) then write(lp,'(a,5i5)') 'n,ia,ja e:',np+j, & iaub(np+j),jaub(np+j),iaui(np+j),jaui(np+j) endif !ldebug_latbdt enddo !j elseif (kdport(l).eq.1) then !northern port np = nxport(l)-ifport(l) j = jfport(l)-1 if (l.eq.1) then nptest = np+ilport(l) endif do i= ifport(l),ilport(l) iapi(np+i) = i japi(np+i) = j ! --- v is normal to the port iavb(np+i) = i javb(np+i) = j+1 iavi(np+i) = i javi(np+i) = j iav2(np+i) = i jav2(np+i) = j-1 ! --- u is tangential to the port iaub(np+i) = i !i=ifport not active jaub(np+i) = j !i=ifport not active iaui(np+i) = i jaui(np+i) = j iau2(np+i) = i !not active jau2(np+i) = j !not active if (l.eq.1 .and. i.eq.ittest) then nptest = np+i endif if (ldebug_latbdt .and. np+i.eq.nptest & .and. mnproc.eq.1) then write(lp,'(a,5i5)') 'n,ia,ja n:',np+i, & iavb(np+i),javb(np+i),iavi(np+i),javi(np+i) endif !ldebug_latbdt enddo !i elseif (kdport(l).eq.2) then !southern port np = nxport(l)-ifport(l) j = jfport(l) if (l.eq.1) then nptest = np+ilport(l) endif do i= ifport(l),ilport(l) iapi(np+i) = i japi(np+i) = j ! --- v is normal to the port iavb(np+i) = i javb(np+i) = j iavi(np+i) = i javi(np+i) = j+1 iav2(np+i) = i jav2(np+i) = j+2 ! --- u is tangential to the port iaub(np+i) = i !i=ifport not active jaub(np+i) = j !i=ifport not active iaui(np+i) = i jaui(np+i) = j iau2(np+i) = i !not active jau2(np+i) = j !not active if (l.eq.1 .and. i.eq.ittest) then nptest = np+i endif if (ldebug_latbdt .and. np+i.eq.nptest & .and. mnproc.eq.1) then write(lp,'(a,5i5)') 'n,ia,ja n:',np+i, & iavb(np+i),javb(np+i),iavi(np+i),javi(np+i) endif !ldebug_latbdt enddo !i endif !kdport enddo !l=1,nports do i= 1,nportpts ! --- find p-point duplicates, if any ndup(i) = 0 do j= 1,nportpts if (i.ne.j .and. iapi(i).eq.iapi(j) & .and. japi(i).eq.japi(j) ) then ndup(i) = j if (ldebug_latbdt .and. mnproc.eq.1) then write(lp,'(a,4i5)') 'n,ndup,ia,ja', & i,j,iapi(i),japi(i) endif !ldebug_latbdt endif enddo !j enddo !i ! allocate(pspeed(nportpts), & rspeed(nportpts), & pline(nportpts), & uline(nportpts), & ulin2(nportpts), & vline(nportpts), & vlin2(nportpts), & plnst(nportpts), & ulnst(nportpts), & vlnst(nportpts), & crs(nportpts), & fin(nportpts) ) ! ! --- Note that 'speed' is in fact qonem*SQRT(g/H). ! --- The qonem allows for the use of pressure fields. ! call xciget(pline,nportpts,depths,iapi,japi,0) ! do l= 1,nports if (kdport(l).eq.4) then ! ! western port ! np = nxport(l)-jfport(l) i = ifport(l) do j= jfport(l),jlport(l) pspeed(np+j) = qonem*sqrt(g/pline(np+j)) rspeed(np+j) = 1.0/pspeed(np+j) if (ldebug_latbdt .and. np+j.eq.nptest & .and. mnproc.eq.1) then write(lp,'(a,i2,4i5,1pe13.5)') & 'w port: ',l,i,j,iapi(np+j),japi(np+j),pspeed(np+j) endif !ldebug_latbdt ! if (i.ge.i0+ 1-nbdy .and. & i.le.i0+ii+nbdy .and. & j.ge.j0+ 1-nbdy .and. & j.le.j0+jj+nbdy ) then iuopn(i-i0,j-j0) = 1 endif enddo ! elseif (kdport(l).eq.3) then ! ! eastern port ! np = nxport(l)-jfport(l) i = ifport(l)-1 do j= jfport(l),jlport(l) pspeed(np+j) = qonem*sqrt(g/pline(np+j)) rspeed(np+j) = 1.0/pspeed(np+j) if (ldebug_latbdt .and. np+j.eq.nptest & .and. mnproc.eq.1) then write(lp,'(a,i2,4i5,1pe13.5)') & 'e port: ',l,i,j,iapi(np+j),japi(np+j),pspeed(np+j) endif !ldebug_latbdt ! if (i+1.ge.i0+ 1-nbdy .and. & i+1.le.i0+ii+nbdy .and. & j .ge.j0+ 1-nbdy .and. & j .le.j0+jj+nbdy ) then iuopn(i-i0+1,j-j0) = 1 endif enddo ! elseif (kdport(l).eq.1) then ! ! northern port ! np = nxport(l)-ifport(l) j = jfport(l)-1 do i= ifport(l),ilport(l) pspeed(np+i) = qonem*sqrt(g/pline(np+i)) rspeed(np+i) = 1.0/pspeed(np+i) if (ldebug_latbdt .and. np+i.eq.nptest & .and. mnproc.eq.1) then write(lp,'(a,i2,4i5,1pe13.5)') & 'n port: ',l,i,j,iapi(np+i),japi(np+i),pspeed(np+i) endif !ldebug_latbdt ! if (i .ge.i0+ 1-nbdy .and. & i .le.i0+ii+nbdy .and. & j+1.ge.j0+ 1-nbdy .and. & j+1.le.j0+jj+nbdy ) then ivopn(i-i0,j-j0+1) = 1 endif enddo ! elseif (kdport(l).eq.2) then ! ! southern port ! np = nxport(l)-ifport(l) j = jfport(l) do i= ifport(l),ilport(l) pspeed(np+i) = qonem*sqrt(g/pline(np+i)) rspeed(np+i) = 1.0/pspeed(np+i) if (ldebug_latbdt .and. np+i.eq.nptest & .and. mnproc.eq.1) then write(lp,'(a,i2,4i5,1pe13.5)') & 's port: ',l,i,j,iapi(np+i),japi(np+i),pspeed(np+i) endif !ldebug_latbdt ! if (i.ge.i0+ 1-nbdy .and. & i.le.i0+ii+nbdy .and. & j.ge.j0+ 1-nbdy .and. & j.le.j0+jj+nbdy ) then ivopn(i-i0,j-j0) = 1 endif enddo ! endif !kdport enddo !l=1,nports if (ldebug_latbdt .and. mnproc.eq.1) then write(lp,*) call flush(lp) endif !ldebug_latbdt if (tidflg.eq.1 .or. tidflg.eq.3) then ! --- last index is 1:3 for z,u,v allocate(port_tide(nportpts,3)) allocate( z_A( nportpts )) allocate( z_R(ncon,nportpts,3),z_I(ncon,nportpts,3)) ! ! --- read in Z, u & v tidal component at open boundary ports ! call latbd_tide(z_A,z_R,z_I,nportpts) endif !tidflg ! ! end of initialization ! call xcsync(flush_lp) return endif !lcount.eq.1 if (.not.allocated(utrans)) then allocate(utrans(nportpts), & vtrans(nportpts) ) utrans(:) = 0.0 vtrans(:) = 0.0 ! call xciget(uline,nportpts,depthu,iaub,jaub,0) call xciget(vline,nportpts,depthv,iavb,javb,0) call xciget(ulin2,nportpts,scuy, iaub,jaub,0) call xciget(vlin2,nportpts,scvx, iavb,javb,0) ! do l= 1,nports if (kdport(l).eq.4) then ! ! western port ! np = nxport(l)-jfport(l) i = ifport(l) do j= jfport(l),jlport(l) utrans(np+j) = uline(np+j)*ulin2(np+j)*qonem if (ltrans_latbdt .and. mnproc.eq.1) then write(lp,'(a,i2,4i5,1p3e13.5)') & 'w tran: ',l,i,j,iaub(np+j),jaub(np+j), & utrans(np+j),uline(np+j)*qonem,ulin2(np+j) endif !ltrans_latbdt enddo !j ! elseif (kdport(l).eq.3) then ! ! eastern port ! np = nxport(l)-jfport(l) i = ifport(l)-1 do j= jfport(l),jlport(l) utrans(np+j) = uline(np+j)*ulin2(np+j)*qonem if (ltrans_latbdt .and. mnproc.eq.1) then write(lp,'(a,i2,4i5,1p3e13.5)') & 'e tran: ',l,i,j,iaub(np+j),jaub(np+j), & utrans(np+j),uline(np+j)*qonem,ulin2(np+j) endif !ltrans_latbdt enddo !j ! elseif (kdport(l).eq.1) then ! ! northern port ! np = nxport(l)-ifport(l) j = jfport(l)-1 do i= ifport(l),ilport(l) vtrans(np+i) = vline(np+i)*vlin2(np+i)*qonem if (ltrans_latbdt .and. mnproc.eq.1) then write(lp,'(a,i2,4i5,1p3e13.5)') & 'n tran: ',l,i,j,iavb(np+i),javb(np+i), & vtrans(np+i),vline(np+i)*qonem,vlin2(np+i) endif !ltrans_latbdt enddo !i ! elseif (kdport(l).eq.2) then ! ! southern port ! np = nxport(l)-ifport(l) j = jfport(l) do i= ifport(l),ilport(l) vtrans(np+i) = vline(np+i)*vlin2(np+i)*qonem if (ltrans_latbdt .and. mnproc.eq.1) then write(lp,'(a,i2,4i5,1p3e13.5)') & 's tran: ',l,i,j,iavb(np+i),javb(np+i), & vtrans(np+i),vline(np+i)*qonem,vlin2(np+i) endif !ltrans_latbdt enddo ! endif !kdport enddo !l=1,nports if (ltrans_latbdt .and. mnproc.eq.1) then write(lp,*) call flush(lp) endif !ldebug_latbdt endif ![uv]trans ! if ((lll.eq.0) .or. (n.eq.0)) return !called from dpthuv ! ! --- nested input only required on first barotropic time step. ! if (lll.eq.1) then sb0 = sports*wb0 sb1 = sports*wb1 do j= 1,jj do i= 1,ii util3(i,j) = pbnest(i,j,lb0)*wb0+pbnest(i,j,lb1)*wb1 util4(i,j) = ubpnst(i,j,lb0)*sb0+ubpnst(i,j,lb1)*sb1 !scaled util5(i,j) = vbpnst(i,j,lb0)*sb0+vbpnst(i,j,lb1)*sb1 !scaled enddo enddo ! call xciget(plnst,nportpts,util3,iapi,japi,mnp) call xciget(ulnst,nportpts,util4,iapi,japi,mnp) call xciget(vlnst,nportpts,util5,iapi,japi,mnp) ! if (tidflg.eq.1 .or. tidflg.eq.3) then d_time = time_8 + lll*dlt/86400.d0 call tides_ports(d_time,nportpts,z_R,z_I,z_A, port_tide) if (mnp.eq.0 .or. mnp.eq.mnproc) then do i= 1,nportpts plnst(i)=plnst(i)+ onem*port_tide(i,1) ulnst(i)=ulnst(i)+0.01d0*port_tide(i,2) vlnst(i)=vlnst(i)+0.01d0*port_tide(i,3) enddo !i endif !mnp if (ldebug_latbdt .and. mnproc.eq.max(mnp,1)) then ! do l= 1,nports do l= 1,1 i=nptest if (kdport(l).eq.4) then write(lp,'(a,i4,a,i4,a,f12.4,3f14.5)') & 'W z,u,v at (',iapi(i),',',japi(i),')=', & d_time,(port_tide(i,k),k=1,3) elseif (kdport(l).eq.3) then write(lp,'(a,i4,a,i4,a,f12.4,3f14.5)') & 'E z,u,v at (',iapi(i),',',japi(i),')=', & d_time,(port_tide(i,k),k=1,3) elseif (kdport(l).eq.1) then write(lp,'(a,i4,a,i4,a,f12.4,3f14.5)') & 'N z,u,v at (',iapi(i),',',japi(i),')=', & d_time,(port_tide(i,k),k=1,3) elseif (kdport(l).eq.2) then write(lp,'(a,i4,a,i4,a,f12.4,3f14.5)') & 'S z,u,v at (',iapi(i),',',japi(i),')=', & d_time,(port_tide(i,k),k=1,3) endif !kdport enddo !l endif !ldebug_latbdt endif !tidflg endif !lll.eq.1 ! call xciget(uline,nportpts,ubavg(1-nbdy,1-nbdy,n),iaui,jaui,mnp) call xciget(ulin2,nportpts,ubavg(1-nbdy,1-nbdy,n),iau2,jau2,mnp) call xciget(vline,nportpts,vbavg(1-nbdy,1-nbdy,n),iavi,javi,mnp) call xciget(vlin2,nportpts,vbavg(1-nbdy,1-nbdy,n),iav2,jav2,mnp) call xciget(pline,nportpts,pbavg(1-nbdy,1-nbdy,n),iapi,japi,mnp) ! ! --- 'wellposed' treatment of pressure and velocity fields. ! --- https://www.hycom.org/attachments/067_boundary.pdf ! --- u* and p* (eqn. 5) are at 1st sea p-point ! --- nested velocity at p-point from 0.5*ui + 0.5*ub ([uv]pnst) ! --- inner velocity at p-point from 1.5*ui - 0.5*ui2 ! --- boundary velocity at vel-point from 2.0*u* - 1.0*ui ! do l= 1,nports ! if (kdport(l).eq.4) then ! ! western port ! if (mnp.eq.0 .or. mnp.eq.mnproc) then trporp(l) = 0.d0 trnest(l) = 0.d0 do j= nxport(l),nxport(l)-1+lnport(l) crs(j)= ulnst(j)+pspeed(j)*plnst(j) fin(j)=1.5*uline(j)-0.5*ulin2(j)-pspeed(j)*pline(j) pline(j)=0.5*(crs(j)-fin(j))*rspeed(j) ulin2(j)= (crs(j)+fin(j)) -uline(j) !temporary array trport(l) = trport(l) + uline(j)*utrans(j) trporp(l) = trporp(l) + 0.5*(crs(j)+fin(j))*utrans(j) trnest(l) = trnest(l) + ulnst(j)*utrans(j) enddo trport(l) = 0.d0 do j= nxport(l),nxport(l)-1+lnport(l) ! smooth the transport jm1 = max(j-1,nxport(l)) jp1 = min(j+1,nxport(l)-1+lnport(l)) uline(j) = (0.25*ulin2(jm1)*utrans(jm1) + & 0.50*ulin2(j )*utrans(j ) + & 0.25*ulin2(jp1)*utrans(jp1) ) / & (0.25* utrans(jm1) + & 0.50* utrans(j ) + & 0.25* utrans(jp1) ) trport(l) = trport(l) + uline(j)*utrans(j) enddo do j= nxport(l)+1,nxport(l)-1+lnport(l) vline(j)=0.5*(vlnst(j-1)+vlnst(j)) enddo endif !mnp ! elseif (kdport(l).eq.3) then ! ! eastern port ! if (ldebug_latbdt .and. & l.eq.1 .and. mnproc.eq.max(mnp,1)) then i=nptest write(lp,'(a,2i5,1pe13.5)') 'e port, uline:', & iaui(i),jaui(i),uline(i) write(lp,'(a,2i5,1pe13.5)') 'e port, ulin2:', & iau2(i),jau2(i),ulin2(i) write(lp,'(a,2i5,1pe13.5)') 'e port, pline:', & iapi(i),japi(i),pline(i) write(lp,'(a,2i5,1pe13.5)') 'e port, plnst:', & iapi(i),japi(i),plnst(i) write(lp,'(a,2i5,1pe13.5)') 'e port, ulnst:', & iapi(i),japi(i),ulnst(i) endif !ldebug_latbdt if (mnp.eq.0 .or. mnp.eq.mnproc) then trporp(l) = 0.d0 trnest(l) = 0.d0 do j= nxport(l),nxport(l)-1+lnport(l) crs(j)= ulnst(j)-pspeed(j)*plnst(j) fin(j)=1.5*uline(j)-0.5*ulin2(j)+pspeed(j)*pline(j) pline(j)=0.5*(fin(j)-crs(j))*rspeed(j) ulin2(j)= (fin(j)+crs(j)) -uline(j) !temporary array trporp(l) = trporp(l) + 0.5*(crs(j)+fin(j))*utrans(j) trnest(l) = trnest(l) + ulnst(j)*utrans(j) enddo trport(l) = 0.d0 do j= nxport(l),nxport(l)-1+lnport(l) ! smooth the transport jm1 = max(j-1,nxport(l)) jp1 = min(j+1,nxport(l)-1+lnport(l)) uline(j) = (0.25*ulin2(jm1)*utrans(jm1) + & 0.50*ulin2(j )*utrans(j ) + & 0.25*ulin2(jp1)*utrans(jp1) ) / & (0.25* utrans(jm1) + & 0.50* utrans(j ) + & 0.25* utrans(jp1) ) trport(l) = trport(l) + uline(j)*utrans(j) enddo do j= nxport(l)+1,nxport(l)-1+lnport(l) vline(j)=0.5*(vlnst(j-1)+vlnst(j)) enddo endif !mnp if (ldebug_latbdt .and. & l.eq.1 .and. mnproc.eq.max(mnp,1)) then i=nptest write(lp,'(a,2i5,1p2e13.5)') 'e port, crs:', & iapi(i),japi(i), crs(i),fin(i) write(lp,'(a,2i5,1p1e13.5)') 'e port, pbavg:', & iapi(i),japi(i),pline(i) write(lp,'(a,2i5,1p1e13.5)') 'e port, ubavg:', & iaub(i),jaub(i),uline(i) write(lp,'(a,2i5,1p1e13.5)') 'e port, vbavg:', & iavb(i),javb(i),vline(i) write(lp,*) call flush(lp) endif !ldebug_latbdt ! elseif (kdport(l).eq.1) then ! ! northern port ! if (ldebug_latbdt .and. & l.eq.1 .and. mnproc.eq.max(mnp,1)) then j=nptest write(lp,'(a,2i5,1pe13.5)') 'n port, vline:', & iavi(j),javi(j),vline(j) write(lp,'(a,2i5,1pe13.5)') 'n port, vlin2:', & iav2(j),jav2(j),vlin2(j) write(lp,'(a,2i5,1p1e13.5)') 'n port, pline:', & iapi(j),japi(j),pline(j) write(lp,'(a,2i5,1p1e13.5)') 'n port, plnst:', & iapi(j),japi(j),vlnst(j) write(lp,'(a,2i5,1p1e13.5)') 'n port, plnst:', & iapi(j),japi(j),vlnst(j) endif !ldebug_latbdt if (mnp.eq.0 .or. mnp.eq.mnproc) then trporp(l) = 0.d0 trnest(l) = 0.d0 do i= nxport(l),nxport(l)-1+lnport(l) crs(i)= vlnst(i)-pspeed(i)*plnst(i) fin(i)=1.5*vline(i)-0.5*vlin2(i)+pspeed(i)*pline(i) pline(i)=0.5*(fin(i)-crs(i))*rspeed(i) vlin2(i)= (fin(i)+crs(i)) -vline(i) !temporary array trporp(l) = trporp(l) + 0.5*(fin(i)+crs(i))*vtrans(i) trnest(l) = trnest(l) + vlnst(i)*vtrans(i) enddo trport(l) = 0.d0 do i= nxport(l),nxport(l)-1+lnport(l) ! smooth the transport im1 = max(i-1,nxport(l)) ip1 = min(i+1,nxport(l)-1+lnport(l)) vline(i) = (0.25*vlin2(im1)*vtrans(im1) + & 0.50*vlin2(i )*vtrans(i ) + & 0.25*vlin2(ip1)*vtrans(ip1) ) / & (0.25* vtrans(im1) + & 0.50* vtrans(i ) + & 0.25* vtrans(ip1) ) trport(l) = trport(l) + vline(i)*vtrans(i) enddo do i= nxport(l)+1,nxport(l)-1+lnport(l) uline(i)=0.5*(ulnst(i-1)+ulnst(i)) enddo endif !mnp if (ldebug_latbdt .and. & l.eq.1 .and. mnproc.eq.max(mnp,1)) then j=nptest write(lp,'(a,2i5,1p2e13.5)') 'n port, crs:', & iapi(j),japi(j), crs(j),fin(j) write(lp,'(a,2i5,1p1e13.5)') 'n port, pbavg:', & iapi(j),japi(j),pline(j) write(lp,'(a,2i5,1p1e13.5)') 'n port, vbavg:', & iavb(j),javb(j),vline(j) write(lp,'(a,2i5,1p1e13.5)') 'n port, ubavg:', & iaub(j),jaub(j),uline(j) write(lp,*) call flush(lp) endif !ldebug_latbdt ! elseif (kdport(l).eq.2) then ! ! southern port ! if (mnp.eq.0 .or. mnp.eq.mnproc) then trporp(l) = 0.d0 trnest(l) = 0.d0 do i= nxport(l),nxport(l)-1+lnport(l) crs(i)= vlnst(i)+pspeed(i)*plnst(i) fin(i)=1.5*vline(i)-0.5*vlin2(i)-pspeed(i)*pline(i) pline(i)=0.5*(crs(i)-fin(i))*rspeed(i) vlin2(i)= (crs(i)+fin(i)) -vline(i) !temporary array trporp(l) = trporp(l) + 0.5*(fin(i)+crs(i))*vtrans(i) trnest(l) = trnest(l) + vlnst(i)*vtrans(i) enddo trport(l) = 0.d0 do i= nxport(l),nxport(l)-1+lnport(l) ! smooth the transport im1 = max(i-1,nxport(l)) ip1 = min(i+1,nxport(l)-1+lnport(l)) vline(i) = (0.25*vlin2(im1)*vtrans(im1) + & 0.50*vlin2(i )*vtrans(i ) + & 0.25*vlin2(ip1)*vtrans(ip1) ) / & (0.25* vtrans(im1) + & 0.50* vtrans(i ) + & 0.25* vtrans(ip1) ) trport(l) = trport(l) + vline(i)*vtrans(i) enddo do i= nxport(l)+1,nxport(l)-1+lnport(l) uline(i)=0.5*(ulnst(i-1)+ulnst(i)) enddo endif !mnp ! endif !kdport ! enddo !l=1,nports ! if (ltrans_latbdt .and. mnproc.eq.1) then if (lll.eq.1) then write(lp,'(i9,i3,a,99f10.5)') & nstep,lll,' trnest =',trnest(1:nports)*1.0d-6 !Sv endif write(lp,'(i9,i3,a,99f10.5)') & nstep,lll,' trport =',trport(1:nports)*1.0d-6 !Sv write(lp,'(i9,i3,a,99f10.5)') & nstep,lll,' trporp =',trporp(1:nports)*1.0d-6 !Sv call flush(lp) endif !transport diagnostic ! ! --- update the boundary points. ! --- the same p-point can appear twice at a boundary corner, ! --- so average the two corner values. ! call xciput(uline,nportpts,ubavg(1-nbdy,1-nbdy,n),iaub,jaub,mnp) call xciput(vline,nportpts,vbavg(1-nbdy,1-nbdy,n),iavb,javb,mnp) if (mod(lll-1,npline).eq.0) then !every npline time steps ulin2(:) = pline(:) !temporary copy of pline do i= 1,nportpts j = ndup(i) if (j.ne.0) then pline(i) = 0.5*(ulin2(i)+ulin2(j)) endif enddo !i call xciput(pline,nportpts,pbavg(1-nbdy,1-nbdy,n),iapi,japi,mnp) endif !npline ! return end subroutine latbdt subroutine latbdtf(n,lll) use mod_xc ! HYCOM communication interface use mod_cb_arrays ! HYCOM saved arrays use mod_tides ! HYCOM Tidal Data implicit none ! integer n,lll ! ! --- apply lateral boundary conditions to barotropic flow field ! ! --- nested sub-region version: ! --- Uses a combination of Flather and clamped open boundary condition. ! ! --- The tangential velocity is not constrained. ! ! --- Tidal z,u & v components added to height and velocity fields ! --- at open boundary points when tidflg = 1 or 3. ! ! --- Note that 'speed' is in fact qonem*SQRT(g/H). ! --- The qonem allows for the use of pressure fields. ! ! --- Note that East, West, North and South refers to the grid ! --- (i.e i,j points) and NOT geographic East, West, North and South ! ! --- the first call is made during initialization. ! ! --- Based on latbdt, with Flather from latbdf. ! --- Alan J. Wallcraft, NRL, March 2012. ! logical, parameter :: ldebug_latbdtf=.false. !usually .false. ! integer, parameter :: mnp=1 !mnflg for xciget and xciput !0 == all tasks; 1 == task 1; integer, parameter :: nchar=120 ! real, parameter :: w_f=1.0 !fraction Flather real, parameter :: w_c=1.0-w_f !fraction clamped SSH and velocity ! logical lfatal,lfatalp integer i,j,k,isec,ifrst,ilast,l,np real aline(nchar) real fatal,sb0,sb1 character*3 char3 ! integer, save :: & nportpts, & nports !number of ports real, save :: & sports !scale factor for nested velocity integer, save, allocatable :: & nxport(:), & kdport(:), & ifport(:),ilport(:), & jfport(:),jlport(:),lnport(:) ! integer jline_start,jline,jtide_start,jtide,tidcon1 ! integer, save, allocatable :: & iaub(:),iavb(:),jaub(:),javb(:), & !boundary vel-point iaui(:),iavi(:),jaui(:),javi(:), & !1st sea vel-point inside boundary iapi(:), japi(:) !1st sea p-point real, save, allocatable :: pspeed(:), & pline(:),uline(:),vline(:), & plnst(:),ulnst(:),vlnst(:) real, save, allocatable :: z_R(:,:,:),z_I(:,:,:),z_A(:) real, save, allocatable :: port_tide(:,:) real*8, save :: d_time ! character*13 fmt save fmt data fmt / '(i4,1x,120i1)' / ! integer lcount save lcount data lcount / 0 / ! lcount = lcount + 1 ! ! --- the first call just initializes data structures. ! if (lcount.eq.1) then if (mnproc.eq.1) then write(lp,'(/a/a)') & 'Flather barotropic nesting', & ' now processing ports.input ...' endif !1st tile call xcsync(flush_lp) ! open(unit=uoff+99,file=trim(flnminp)//'ports.input') ! ! --- 'sports' = scale factor for nested velocity (optional, default 1.0) ! --- 'nports' = number of boundary port sections. call blkinr2(sports,i, & 'sports','(a6," =",f10.4," ")', & 'nports','(a6," =",f11.0)') !treat nports as real if (i.eq.1) then !sports call blkini(nports,'nports') else !nports nports = sports sports = 1.0 !default endif if (mnproc.eq.1) then write(lp,*) endif if (nports.eq.0) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdtf - must have lbflag=0 for nports=0' write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbdtf)') stop '(latbdtf)' elseif (nports.lt.0) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdtf - illegal nports value' write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbdtf)') stop '(latbdtf)' endif ! allocate( & nxport(nports), & kdport(nports), & ifport(nports), & ilport(nports), & jfport(nports), & jlport(nports), & lnport(nports) ) ! ! --- read in the ports one at a time ! nportpts = 0 do l= 1,nports ! ! --- port location is w.r.t. u (EW) or v (NS) grid ! --- and identifies the sea at the port ! --- the minimum index is 0 ! ! --- 'kdport' = port orientation (1=N, 2=S, 3=E, 4=W) ! --- 'ifport' = first i-index ! --- 'ilport' = last i-index (=ifport for N or S orientation) ! --- 'jfport' = first j-index ! --- 'jlport' = last j-index (=jfport for E or W orientation) ! --- 'lnport' = port length (calculated, not input) call blkini(kdport(l),'kdport') call blkini(ifport(l),'ifport') call blkini(ilport(l),'ilport') call blkini(jfport(l),'jfport') call blkini(jlport(l),'jlport') if (mnproc.eq.1) then write(lp,*) endif ! nxport(l) = nportpts+1 !start of this port lnport(l) = ilport(l)-ifport(l)+ & jlport(l)-jfport(l)+1 !length of this port nportpts = nportpts+lnport(l) !length of all ports ! ! --- sanity check. ! if (kdport(l).gt.2) then if (ifport(l).ne.ilport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdtf - port direction', & ' and orientation are not consistent' write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbdtf)') stop '(latbdtf)' endif else if (jfport(l).ne.jlport(l)) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdtf - port direction', & ' and orientation are not consistent' write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbdtf)') stop '(latbdtf)' endif endif if (ifport(l).gt.ilport(l) .or. & jfport(l).gt.jlport(l) ) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbdtf - port', & ' location is not consistent' write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbdtf)') stop '(latbdtf)' endif enddo !l ! close(unit=uoff+99) ! ! --- check ports against masks, ! --- mark the port locations on masks and print them out. ! lfatal = .false. do l= 1,nports lfatalp = .false. ! if (kdport(l).eq.4) then ! ! western port ! i = ifport(l) do j= jfport(l),jlport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.1 .or. j.gt.jtdm ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iu(i-i0,j-j0).ne.0) then lfatalp = .true. iu(i-i0,j-j0) = 9 !indicate an error else iu(i-i0,j-j0) = -1 endif if (iu(i-i0+1,j-j0).ne.1 .or. & iu(i-i0+2,j-j0).ne.1 ) then lfatalp = .true. iu(i-i0,j-j0) = 7 !indicate an error endif enddo ! elseif (kdport(l).eq.3) then ! ! eastern port ! i = ifport(l) do j= jfport(l),jlport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.1 .or. j.gt.jtdm ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iu(i-i0,j-j0).ne.0) then lfatalp = .true. iu(i-i0,j-j0) = 9 !indicate an error else iu(i-i0,j-j0) = -1 endif if (iu(i-i0-1,j-j0).ne.1 .or. & iu(i-i0-2,j-j0).ne.1 ) then lfatalp = .true. iu(i-i0,j-j0) = 7 !indicate an error endif enddo ! elseif (kdport(l).eq.1) then ! ! northern port ! j = jfport(l) do i= ifport(l),ilport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.3 .or. j.gt.jtdm ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iv(i-i0,j-j0).ne.0) then lfatalp = .true. iv(i-i0,j-j0) = 9 !indicate an error else iv(i-i0,j-j0) = -1 endif if (iv(i-i0,j-j0-1).ne.1 .or. & iv(i-i0,j-j0-2).ne.1 ) then lfatalp = .true. iv(i-i0,j-j0) = 7 !indicate an error endif enddo ! elseif (kdport(l).eq.2) then ! ! southern port ! j = jfport(l) do i= ifport(l),ilport(l) if (i.lt.1 .or. i.gt.itdm .or. & j.lt.1 .or. j.gt.jtdm-2 ) then lfatalp = .true. elseif (i.le.i0 .or. i.gt.i0+ii .or. & j.le.j0 .or. j.gt.j0+jj ) then cycle ! not on this tile. elseif (iv(i-i0,j-j0).ne.0) then lfatalp = .true. iv(i-i0,j-j0) = 9 !indicate an error else iv(i-i0,j-j0) = -1 endif if (iv(i-i0,j-j0+1).ne.1 .or. & iv(i-i0,j-j0+2).ne.1 ) then lfatalp = .true. iv(i-i0,j-j0) = 7 !indicate an error endif enddo ! endif ! if (lfatalp) then write(lp,*) write(lp,*) 'error in latbdtf - port ',l,' mislocated', & ' (mnproc = ',mnproc,')' write(lp,*) call flush(lp) endif lfatal = lfatal .or. lfatalp enddo !l=1,nports ! ! local lfatal to global lfatal ! if (lfatal) then fatal = 1.0 else fatal = 0.0 endif call xcmaxr(fatal) lfatal = fatal.gt.0.5 ! ! --- write out -iu- and -iv- arrays, if they are not too big ! --- data are written in strips nchar points wide ! if (lfatal .or. max(itdm,jtdm).le.2*nchar) then util1(1:ii,1:jj) = iu(1:ii,1:jj) ! xclget is for real arrays isec=(itdm-1)/nchar do ifrst=0,nchar*isec,nchar ilast=min(itdm,ifrst+nchar) write (char3,'(i3)') ilast-ifrst fmt(8:10)=char3 if (mnproc.eq.1) then write (lp,'(a,i5,a,i5)') & 'iu array, cols',ifrst+1,' --',ilast endif !1st tile do j= jtdm,1,-1 call xclget(aline,ilast-ifrst, util1,ifrst+1,j,1,0, 1) if (mnproc.eq.1) then write (lp,fmt) j,(nint(aline(i)),i=1,ilast-ifrst) endif !1st tile enddo enddo if (mnproc.eq.1) then write (lp,*) endif call xcsync(flush_lp) ! util1(1:ii,1:jj) = iv(1:ii,1:jj) ! xclget is for real arrays isec=(itdm-1)/nchar do ifrst=0,nchar*isec,nchar ilast=min(itdm,ifrst+nchar) write (char3,'(i3)') ilast-ifrst fmt(8:10)=char3 if (mnproc.eq.1) then write (lp,'(a,i5,a,i5)') & 'iv array, cols',ifrst+1,' --',ilast endif !1st tile do j= jtdm,1,-1 call xclget(aline,ilast-ifrst, util1,ifrst+1,j,1,0, 1) if (mnproc.eq.1) then write (lp,fmt) j,(nint(aline(i)),i=1,ilast-ifrst) endif !1st tile enddo enddo if (mnproc.eq.1) then write (lp,*) endif call xcsync(flush_lp) endif ! small region ! if (lfatal) then write(lp,*) write(lp,*) 'error in latbdtf - bad port(s)' write(lp,*) call flush(lp) call xchalt('(latbdtf)') stop '(latbdtf)' endif ! ! --- restore iu and iv, and zero iuopn and ivopn. ! !$OMP PARALLEL DO PRIVATE(j,i) & !$OMP SCHEDULE(STATIC,jblk) do j= 1,jj do i= 1,ii iu(i,j) = max( iu(i,j), 0 ) iv(i,j) = max( iv(i,j), 0 ) enddo enddo !$OMP PARALLEL DO PRIVATE(j,i) & !$OMP SCHEDULE(STATIC,jblk) do j= 1-nbdy,jj+nbdy do i= 1-nbdy,ii+nbdy iuopn(i,j) = 0 ivopn(i,j) = 0 enddo enddo ! ! --- initialize the ports ! allocate(iaub(nportpts),jaub(nportpts), & iavb(nportpts),javb(nportpts), & iaui(nportpts),jaui(nportpts), & iavi(nportpts),javi(nportpts), & iapi(nportpts),japi(nportpts) ) ! do l= 1,nports if (kdport(l).eq.4) then !western port np = nxport(l)-jfport(l) i = ifport(l) do j= jfport(l),jlport(l) iaub(np+j) = i jaub(np+j) = j iaui(np+j) = i+1 jaui(np+j) = j iapi(np+j) = i japi(np+j) = j ! --- v not active iavb(np+j) = i javb(np+j) = j iavi(np+j) = i javi(np+j) = j if (ldebug_latbdtf .and. mnproc.eq.1) then write(lp,'(a,5i5)') 'n,ia,ja w:',np+j, & iaub(np+j),jaub(np+j),iaui(np+j),jaui(np+j) endif !ldebug_latbdtf enddo !j elseif (kdport(l).eq.3) then !eastern port np = nxport(l)-jfport(l) i = ifport(l)-1 do j= jfport(l),jlport(l) iaub(np+j) = i+1 jaub(np+j) = j iaui(np+j) = i jaui(np+j) = j iapi(np+j) = i japi(np+j) = j ! --- v not active iavb(np+j) = i javb(np+j) = j iavi(np+j) = i javi(np+j) = j if (ldebug_latbdtf .and. mnproc.eq.1) then write(lp,'(a,5i5)') 'n,ia,ja e:',np+j, & iaub(np+j),jaub(np+j),iaui(np+j),jaui(np+j) endif !ldebug_latbdtf enddo !j elseif (kdport(l).eq.1) then !northern port np = nxport(l)-ifport(l) j = jfport(l)-1 do i= ifport(l),ilport(l) iavb(np+i) = i javb(np+i) = j+1 iavi(np+i) = i javi(np+i) = j iapi(np+i) = i japi(np+i) = j ! --- u not active iaub(np+i) = i jaub(np+i) = j iaui(np+i) = i jaui(np+i) = j if (ldebug_latbdtf .and. mnproc.eq.1) then write(lp,'(a,5i5)') 'n,ia,ja n:',np+i, & iavb(np+i),javb(np+i),iavi(np+i),javi(np+i) endif !ldebug_latbdtf enddo !i elseif (kdport(l).eq.2) then !southern port np = nxport(l)-ifport(l) j = jfport(l) do i= ifport(l),ilport(l) iavb(np+i) = i javb(np+i) = j iavi(np+i) = i javi(np+i) = j+1 iapi(np+i) = i japi(np+i) = j ! --- u not active iaub(np+i) = i jaub(np+i) = j iaui(np+i) = i jaui(np+i) = j if (ldebug_latbdtf .and. mnproc.eq.1) then write(lp,'(a,5i5)') 'n,ia,ja n:',np+i, & iavb(np+i),javb(np+i),iavi(np+i),javi(np+i) endif !ldebug_latbdtf enddo !i endif !kdport enddo !l=1,nports ! allocate(pspeed(nportpts), & pline(nportpts), & uline(nportpts), & vline(nportpts), & plnst(nportpts), & ulnst(nportpts), & vlnst(nportpts) ) ! ! --- Note that 'speed' is in fact qonem*SQRT(g/H). ! --- The qonem allows for the use of pressure fields. ! call xciget(pline,nportpts,depths,iapi,japi,0) ! do l= 1,nports if (kdport(l).eq.4) then ! ! western port ! np = nxport(l)-jfport(l) i = ifport(l) do j= jfport(l),jlport(l) pspeed(np+j) = qonem*sqrt(g/pline(np+j)) if (ldebug_latbdtf .and. mnproc.eq.1) then write(lp,'(a,i2,2i5,1pe13.5)') & 'w port: ',l,i,j,pspeed(np+j) endif !ldebug_latbdtf ! if (i.ge.i0+ 1-nbdy .and. & i.le.i0+ii+nbdy .and. & j.ge.j0+ 1-nbdy .and. & j.le.j0+jj+nbdy ) then iuopn(i-i0,j-j0) = 1 endif enddo ! elseif (kdport(l).eq.3) then ! ! eastern port ! np = nxport(l)-jfport(l) i = ifport(l)-1 do j= jfport(l),jlport(l) pspeed(np+j) = qonem*sqrt(g/pline(np+j)) if (ldebug_latbdtf .and. mnproc.eq.1) then write(lp,'(a,i2,2i5,1pe13.5)') & 'e port: ',l,i,j,pspeed(np+j) endif !ldebug_latbdtf ! if (i+1.ge.i0+ 1-nbdy .and. & i+1.le.i0+ii+nbdy .and. & j .ge.j0+ 1-nbdy .and. & j .le.j0+jj+nbdy ) then iuopn(i-i0+1,j-j0) = 1 endif enddo ! elseif (kdport(l).eq.1) then ! ! northern port ! np = nxport(l)-ifport(l) j = jfport(l)-1 do i= ifport(l),ilport(l) pspeed(np+i) = qonem*sqrt(g/pline(np+i)) if (ldebug_latbdtf .and. mnproc.eq.1) then write(lp,'(a,i2,2i5,1pe13.5)') & 'n port: ',l,i,j,pspeed(np+i) endif !ldebug_latbdtf ! if (i .ge.i0+ 1-nbdy .and. & i .le.i0+ii+nbdy .and. & j+1.ge.j0+ 1-nbdy .and. & j+1.le.j0+jj+nbdy ) then ivopn(i-i0,j-j0+1) = 1 endif enddo ! elseif (kdport(l).eq.2) then ! ! southern port ! np = nxport(l)-ifport(l) j = jfport(l) do i= ifport(l),ilport(l) pspeed(np+i) = qonem*sqrt(g/pline(np+i)) if (ldebug_latbdtf .and. mnproc.eq.1) then write(lp,'(a,i2,2i5,1pe13.5)') & 's port: ',l,i,j,pspeed(np+i) endif !ldebug_latbdtf ! if (i.ge.i0+ 1-nbdy .and. & i.le.i0+ii+nbdy .and. & j.ge.j0+ 1-nbdy .and. & j.le.j0+jj+nbdy ) then ivopn(i-i0,j-j0) = 1 endif enddo ! endif !kdport enddo !l=1,nports if (ldebug_latbdtf .and. mnproc.eq.1) then write(lp,*) call flush(lp) endif !ldebug_latbdtf if (tidflg.eq.1 .or. tidflg.eq.3) then ! --- last index is 1:3 for z,u,v allocate(port_tide(nportpts,3)) allocate( z_A( nportpts )) allocate( z_R(ncon,nportpts,3),z_I(ncon,nportpts,3)) ! ! --- read in Z, u & v tidal component at open boundary ports ! call latbd_tide(z_A,z_R,z_I,nportpts) endif !tidflg ! ! end of initialization ! call xcsync(flush_lp) return endif if ((lll.eq.0) .or. (n.eq.0)) return !called from dpthuv ! ! --- nested input only required on first barotropic time step. ! --- [uv]lnst on 1st sea p-point for compatibility with latbdt. ! if (lll.eq.1) then sb0 = sports*wb0 sb1 = sports*wb1 do j= 1,jj do i= 1,ii util1(i,j) = ubpnst(i,j,lb0)*sb0+ubpnst(i,j,lb1)*sb1 !scaled util2(i,j) = vbpnst(i,j,lb0)*sb0+vbpnst(i,j,lb1)*sb1 !scaled util3(i,j) = pbnest(i,j,lb0)*wb0+pbnest(i,j,lb1)*wb1 enddo enddo ! call xciget(ulnst,nportpts,util1,iapi,japi,mnp) call xciget(vlnst,nportpts,util2,iapi,japi,mnp) call xciget(plnst,nportpts,util3,iapi,japi,mnp) ! if (tidflg.eq.1 .or. tidflg.eq.3) then d_time = time_8 + lll*dlt/86400.d0 call tides_ports(d_time,nportpts,z_R,z_I,z_A, port_tide) if (mnp.eq.0 .or. mnp.eq.mnproc) then do i= 1,nportpts plnst(i)=plnst(i)+ onem*port_tide(i,1) ulnst(i)=ulnst(i)+0.01d0*port_tide(i,2) vlnst(i)=vlnst(i)+0.01d0*port_tide(i,3) enddo !i endif !mnp if (ldebug_latbdtf .and. mnproc.eq.max(mnp,1)) then do l= 1,nports i=nxport(l)-1 + lnport(l)/2 if (kdport(l).eq.4) then write(lp,'(a,i4,a,i4,a,f12.4,3f14.5)') & 'W z,u,v at (',iapi(i),',',japi(i),')=', & d_time,(port_tide(i,k),k=1,3) elseif (kdport(l).eq.3) then write(lp,'(a,i4,a,i4,a,f12.4,3f14.5)') & 'E z,u,v at (',iapi(i),',',japi(i),')=', & d_time,(port_tide(i,k),k=1,3) elseif (kdport(l).eq.1) then write(lp,'(a,i4,a,i4,a,f12.4,3f14.5)') & 'N z,u,v at (',iapi(i),',',japi(i),')=', & d_time,(port_tide(i,k),k=1,3) elseif (kdport(l).eq.2) then write(lp,'(a,i4,a,i4,a,f12.4,3f14.5)') & 'S z,u,v at (',iapi(i),',',japi(i),')=', & d_time,(port_tide(i,k),k=1,3) endif !kdport enddo !l=1,nports endif !ldebug_latbdtf endif !tidflg endif !lll.eq.1 ! call xciget(uline,nportpts,ubavg(1-nbdy,1-nbdy,n),iaui,jaui,mnp) call xciget(vline,nportpts,vbavg(1-nbdy,1-nbdy,n),iavi,javi,mnp) call xciget(pline,nportpts,pbavg(1-nbdy,1-nbdy,n),iapi,japi,mnp) ! ! --- 'wellposed' treatment of pressure and velocity fields. ! --- u* and p* are at 1st sea p-point ! --- boundary velocity at vel-point from 2.0*u* - 1.0*ui ! do l= 1,nports ! if (kdport(l).eq.4) then ! ! western port ! if (ldebug_latbdtf .and. mnproc.eq.max(mnp,1)) then j=jfport(l)-1 + lnport(l)/2 i=nxport(l)-1 + lnport(l)/2 write(lp,'(a,i4,e14.5)') 'w port, uline:',j,uline(i) write(lp,'(a,i4,e14.5)') 'w port, pline:',j,pline(i)*qonem write(lp,'(a,i4,e14.5)') 'w port, plnst:',j,plnst(i)*qonem write(lp,'(a,i4,e14.5)') 'w port, ulnst:',j,ulnst(i) endif !1st tile if (mnp.eq.0 .or. mnp.eq.mnproc) then do j= nxport(l),nxport(l)-1+lnport(l) pline(j)=w_f*pline(j)+w_c*plnst(j) uline(j)=2.0*(ulnst(j)+pspeed(j)*(plnst(j)-pline(j))) & !new pline -uline(j) enddo !j endif !mnp if (ldebug_latbdtf .and. mnproc.eq.max(mnp,1)) then j=jfport(l)-1 + lnport(l)/2 i=nxport(l)-1 + lnport(l)/2 write(lp,'(a,i4,e14.5)') 'e port, pbavg:',j,pline(i)*qonem write(lp,'(a,i4,e14.5)') 'e port, ubavg:',j,uline(i) write(lp,*) call flush(lp) endif !1st tile ! elseif (kdport(l).eq.3) then ! ! eastern port ! if (ldebug_latbdtf .and. mnproc.eq.max(mnp,1)) then j=jfport(l)-1 + lnport(l)/2 i=nxport(l)-1 + lnport(l)/2 write(lp,'(a,i4,e14.5)') 'e port, uline:',j,uline(i) write(lp,'(a,i4,e14.5)') 'e port, pline:',j,pline(i)*qonem write(lp,'(a,i4,e14.5)') 'e port, plnst:',j,plnst(i)*qonem write(lp,'(a,i4,e14.5)') 'e port, ulnst:',j,ulnst(i) endif !1st tile if (mnp.eq.0 .or. mnp.eq.mnproc) then do j= nxport(l),nxport(l)-1+lnport(l) pline(j)=w_f*pline(j)+w_c*plnst(j) uline(j)=2.0*(ulnst(j)+pspeed(j)*(pline(j)-plnst(j))) & !new pline -uline(j) enddo endif !mnp if (ldebug_latbdtf .and. mnproc.eq.max(mnp,1)) then j=jfport(l)-1 + lnport(l)/2 i=nxport(l)-1 + lnport(l)/2 write(lp,'(a,i4,e14.5)') 'e port, pbavg:',j,pline(i)*qonem write(lp,'(a,i4,e14.5)') 'e port, ubavg:',j,uline(i) write(lp,*) call flush(lp) endif !1st tile ! elseif (kdport(l).eq.1) then ! ! northern port ! if (ldebug_latbdtf .and. mnproc.eq.max(mnp,1)) then i=ifport(l)-1 + lnport(l)/2 j=nxport(l)-1 + lnport(l)/2 write(lp,'(a,i4,e14.5)') 'n port, vline:',i,vline(j) write(lp,'(a,i4,e14.5)') 'n port, pline:',i,pline(j)*qonem write(lp,'(a,i4,e14.5)') 'n port, plnst:',i,plnst(j)*qonem write(lp,'(a,i4,e14.5)') 'n port, vlnst:',i,vlnst(j) endif !ldebug_latbdtf if (mnp.eq.0 .or. mnp.eq.mnproc) then do i= nxport(l),nxport(l)-1+lnport(l) pline(i)=w_f*pline(i)+w_c*plnst(i) vline(i)=2.0*(vlnst(i)+pspeed(i)*(pline(i)-plnst(i))) & !new pline -vline(i) enddo endif !mnp if (ldebug_latbdtf .and. mnproc.eq.max(mnp,1)) then i=ifport(l)-1 + lnport(l)/2 j=nxport(l)-1 + lnport(l)/2 write(lp,'(a,i4,e14.5)') 'n port, pbavg:',i,pline(j)*qonem write(lp,'(a,i4,e14.5)') 'n port, vbavg:',i,vline(j) write(lp,*) call flush(lp) endif !ldebug_latbdtf ! elseif (kdport(l).eq.2) then ! ! southern port ! if (ldebug_latbdtf .and. mnproc.eq.max(mnp,1)) then i=ifport(l)-1 + lnport(l)/2 j=nxport(l)-1 + lnport(l)/2 write(lp,'(a,i4,e14.5)') 's port, vline:',i,vline(j) write(lp,'(a,i4,e14.5)') 's port, pline:',i,pline(j)*qonem write(lp,'(a,i4,e14.5)') 's port, plnst:',i,plnst(j)*qonem write(lp,'(a,i4,e14.5)') 's port, vlnst:',i,vlnst(j) endif !ldebug_latbdtf if (mnp.eq.0 .or. mnp.eq.mnproc) then do i= nxport(l),nxport(l)-1+lnport(l) pline(i)=w_f*pline(i)+w_c*plnst(i) vline(i)=2.0*(vlnst(i)+pspeed(i)*(plnst(i)-pline(i))) & !new pline -vline(i) enddo endif !mnp if (ldebug_latbdtf .and. mnproc.eq.max(mnp,1)) then i=ifport(l)-1 + lnport(l)/2 j=nxport(l)-1 + lnport(l)/2 write(lp,'(a,i4,e14.5)') 's port, pbavg:',i,pline(j)*qonem write(lp,'(a,i4,e14.5)') 's port, vbavg:',i,vline(j) write(lp,*) call flush(lp) endif !ldebug_latbdtf ! endif !kdport ! enddo !l=1,nports ! ! --- update the boundary points. ! if (w_f.ne.1.0) then call xciput(pline,nportpts,pbavg(1-nbdy,1-nbdy,n),iapi,japi,mnp) endif !clamped call xciput(uline,nportpts,ubavg(1-nbdy,1-nbdy,n),iaub,jaub,mnp) call xciput(vline,nportpts,vbavg(1-nbdy,1-nbdy,n),iavb,javb,mnp) ! return end subroutine latbdtf subroutine latbd_tide(z_A,z_R,z_I,nportpts) use mod_xc ! HYCOM communication interface use mod_cb_arrays ! HYCOM saved arrays use mod_tides ! HYCOM Tidal Data implicit none ! integer nportpts real z_A(nportpts),z_R(ncon,nportpts,3),z_I(ncon,nportpts,3) ! ! read in Z, u & v tidal component at open boundary ports ! ! These are generated using OTPSnc (or OTPS2) extract_HC ! http://volkov.oce.orst.edu/tides/otps.html ! http://volkov.oce.orst.edu/tides/tpxo8_atlas.html ! logical, parameter :: ldebug_latbd_t=.false. !usually .false. ! integer i,j real dum1,dum2 logical Tide_Modes_Correct character*156 Tide_Line ! Character*2 Tide_Names(8) save Tide_Names data Tide_Names/'m2','s2','k1','o1','n2','p1','k2','q1'/ ! ! Default is zero ! z_R(:,:,:) = 0.0 z_I(:,:,:) = 0.0 ! ! Read Pang (angle of xward wrt eward) for each port point ! open(unit=uoff+99,file=trim(flnminp)//'ports_a.input') read(uoff+99,'(a)')Tide_Line if (ldebug_latbd_t .and. mnproc.eq.1) then write(lp,'(a)') trim(Tide_Line) endif !ldebug_latbd_t if (Tide_Line(1:25).ne. & ' Lat Lon | Pang') then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbd_tide - ports_a.input' write(lp,*) 'Initial Line not correct!' write(lp,*) 'Expecting:'// & ' Lat Lon | Pang' write(lp,*) ' Got:'//Tide_Line(1:25) call flush(lp) endif !1st tile call xcstop('(latbd_tide)') stop '(latbd_tide)' endif !1st line do i= 1, nportpts read(uoff+99,'(2F9.3,F11.6)')dum1,dum2,z_A(i) if (ldebug_latbd_t .and. mnproc.eq.1) then write(6,'(a,i4,a,F11.6)') & 'Port Point ',i,' pang= ',z_A(i) endif !ldebug_latbd_t enddo !i close(unit=uoff+99) if (ldebug_latbd_t .and. mnproc.eq.1) then write(lp,*)'Pang read for',nportpts,' open boundary points' endif !ldebug_latbd_t ! ! Read z tidal components (real, imag) for each port point ! open(unit=uoff+99,file=trim(flnminp)//'ports_z.input') read(uoff+99,'(a)')Tide_Line if (mnproc.eq.1) then write(lp,'(a)') trim(Tide_Line) endif !1st tile read(uoff+99,'(a)')Tide_Line if (mnproc.eq.1) then write(lp,'(a)') trim(Tide_Line) endif !1st tile if (Tide_Line(1:15).ne.' Elevations (m)') then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbd_tide - ports_z.input' write(lp,*) 'Second Line not correct!' write(lp,*) 'Expecting: Elevations (m)' write(lp,*) ' Got:'//Tide_Line(1:15) write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbd_tide)') stop '(latbd_tide)' endif !2nd line read(uoff+99,'(a)')Tide_Line if (ldebug_latbd_t .and. mnproc.eq.1) then write(lp,'(a)') trim(Tide_Line) endif !ldebug_latbd_t ! --- always have an OTPSnc header, data can be from OTPS2 if (Tide_Line(1:50).ne. & ' Lat Lon | m2_Re m2_Im s2_Re s2_Im') then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbd_tide - ports_z.input' write(lp,*) 'Third Line not correct!' write(lp,*) 'Expecting:'// & ' Lat Lon | m2_Re m2_Im s2_Re s2_Im' write(lp,*) ' Got:'//Tide_Line(1:50) if (Tide_Line(22:27).eq.'m2_amp') then write(lp,*)'ports_z.input file is in Amplitude, Phase Mode!' endif write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbd_tide)') stop '(latbd_tide)' endif !3rd line Tide_Modes_Correct=.true. do i=1,8 Tide_Modes_Correct=Tide_Modes_Correct .and. & Tide_Line( 6+16*i:10+16*i).eq.Tide_Names(i)//'_Re' .and. & Tide_Line(14+16*i:18+16*i).eq.Tide_Names(i)//'_Im' enddo if (.not.Tide_Modes_Correct) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbd_tide - ports_z.input' write(lp,*) 'Tidal Modes may be in wrong order!' write(lp,*) 'Expecting: m2,s2,k1,o1,n2,p1,k2,q1' write(lp,*) ' Got: ',(Tide_Line(6+16*i:8+16*i),i=1,8) if (Tide_Line(22:27).eq.'m2_amp') then write(lp,*) & 'ports_z.input file is in Amplitude, Phase Mode!' endif write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbd_tide)') stop '(latbd_tide)' endif !.not.Tide_Modes_Correct do i= 1, nportpts read(uoff+99,'(a)')Tide_Line if (index(Tide_Line,'Site').eq.0) then if (len_trim(Tide_Line).eq.147) then ! --- from OTPS2 read(Tide_Line,'(F9.4,F10.4,16F8.3)')dum1,dum2, & (z_R(j,i,1),z_I(j,i,1),j=1,ncon) write(lp,*) Tide_Line write(lp,*) 'luana' else ! --- from OTPSnc read(Tide_Line,'(2F9.3,16F8.3)')dum1,dum2, & (z_R(j,i,1),z_I(j,i,1),j=1,ncon) endif elseif (mnproc.eq.1) then write(lp,'(a,i5,a / a)') & 'WARNING: port location',i,' treated as all zeros', & trim(Tide_Line) endif if (ldebug_latbd_t .and. mnproc.eq.1) then write(6,'(a,i4,a,8F8.3/23x,8F8.3)') & 'Port Point ',i,' z(j)= ', & (z_R(j,i,1),j=1,ncon), & (z_I(j,i,1),j=1,ncon) endif !ldebug_latbd_t enddo !i close(unit=uoff+99) if (ldebug_latbd_t .and. mnproc.eq.1) then write(lp,'(18a)')'z tidal Modes:', & (' '//Tide_Line(14+8*i:17+8*i),i=1,16) write(lp,*)'read for',nportpts,' open boundary points' endif !ldebug_latbd_t ! ! Read u tidal components (real, imag) for each port point open(unit=uoff+99,file=trim(flnminp)//'ports_u.input') read(uoff+99,'(a)')Tide_Line if (mnproc.eq.1) then write(lp,'(a)') trim(Tide_Line) endif !1st tile read(uoff+99,'(a)')Tide_Line if (mnproc.eq.1) then write(lp,'(a)') trim(Tide_Line) endif !1st tile if (Tide_Line(1:20).ne.' WE velocity (cm/s)') then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbd_tide - ports_u.input' write(lp,*) 'Second Line not correct!' write(lp,*) 'Expecting: WE velocity (cm/s)' write(lp,*) ' Got:'//Tide_Line(1:20) write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbd_tide)') stop '(latbd_tide)' endif !2nd line read(uoff+99,'(a)')Tide_Line if (ldebug_latbd_t .and. mnproc.eq.1) then write(lp,'(a)') trim(Tide_Line) endif !ldebug_latbd_t ! --- always have an OTPSnc header, data can be from OTPS2 if (Tide_Line(1:50).ne. & ' Lat Lon | m2_Re m2_Im s2_Re s2_Im') then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbd_tide - ports_u.input' write(lp,*) 'Third Line not correct!' write(lp,*) 'Expecting:'// & ' Lat Lon | m2_Re m2_Im s2_Re s2_Im' write(lp,*) ' Got:'//Tide_Line(1:50) if (Tide_Line(22:27).eq.'m2_amp') then write(lp,*)'ports_u.input file is in Amplitude, Phase Mode?' write(lp,*)'Expecting mode_Re mode_Im data' endif write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbd_tide)') stop '(latbd_tide)' endif !3rd line Tide_Modes_Correct=.true. do i=1,8 Tide_Modes_Correct=Tide_Modes_Correct .and. & Tide_Line( 6+16*i:10+16*i).eq.Tide_Names(i)//'_Re' .and. & Tide_Line(14+16*i:18+16*i).eq.Tide_Names(i)//'_Im' enddo !i if (.not.Tide_Modes_Correct) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbd_tide - ports_u.input' write(lp,*) 'Tidal Modes may be in wrong order!' write(lp,*) 'Expecting: m2,s2,k1,o1,n2,p1,k2,q1' write(lp,*) ' Got: ',(Tide_Line(6+16*i:8+16*i),i=1,8) if (Tide_Line(22:27).eq.'m2_amp') then write(lp,*) & 'ports_u.input file is in Amplitude, Phase Mode!' write(lp,*) & 'Expecting mode_Re mode_Im data' endif write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbd_tide)') stop '(latbd_tide)' endif !.not.Tide_Modes_Correct do i= 1,nportpts read(uoff+99,'(a)')Tide_Line if (index(Tide_Line,'Site').eq.0) then if (len_trim(Tide_Line).eq.147) then ! --- from OTPS2 read(Tide_Line,'(F9.4,F10.4,16F8.3)')dum1,dum2, & (z_R(j,i,2),z_I(j,i,2),j=1,ncon) else ! --- from OTPSnc read(Tide_Line,'(2F9.3,16F8.3)')dum1,dum2, & (z_R(j,i,2),z_I(j,i,2),j=1,ncon) endif elseif (mnproc.eq.1) then write(lp,'(a,i5,a / a)') & 'WARNING: port location',i,' treated as all zeros', & trim(Tide_Line) endif if (ldebug_latbd_t .and. mnproc.eq.1) then write(6,'(a,i4,a,8F8.3/23x,8F8.3)') & 'Port Point ',i,' z(j)= ', & (z_R(j,i,2),j=1,ncon), & (z_I(j,i,2),j=1,ncon) endif !ldebug_latbd_t enddo !i close(unit=uoff+99) if (ldebug_latbd_t .and. mnproc.eq.1) then write(lp,'(18a)')'u tidal Modes:', & (' '//Tide_Line(14+8*i:17+8*i),i=1,16) write(lp,*)'read for',nportpts,' open boundary points' endif !ldebug_latbd_t ! ! Read v tidal components (real, imag) for each port point ! open(unit=uoff+99,file=trim(flnminp)//'ports_v.input') read(uoff+99,'(a)')Tide_Line if (mnproc.eq.1) then write(lp,'(a /)') trim(Tide_Line) endif !1st tile read(uoff+99,'(a)')Tide_Line if (mnproc.eq.1) then write(lp,'(a)') trim(Tide_Line) endif !1st tile if (Tide_Line(1:20).ne.' SN velocity (cm/s)') then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbd_tide - ports_v.input' write(lp,*) 'Second Line not correct!' write(lp,*) 'Expecting: SN velocity (cm/s)' write(lp,*) ' Got:'//Tide_Line(1:20) write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbd_tide)') stop '(latbd_tide)' endif !2nd line read(uoff+99,'(a)')Tide_Line if (ldebug_latbd_t .and. mnproc.eq.1) then write(lp,'(a)') trim(Tide_Line) endif !ldebug_latbd_t ! --- always have an OTPSnc header, data can be from OTPS2 if (Tide_Line(1:50).ne. & ' Lat Lon | m2_Re m2_Im s2_Re s2_Im') then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbd_tide - ports_v.input' write(lp,*) 'Third Line not correct!' write(lp,*) 'Expecting:'// & ' Lat Lon | m2_Re m2_Im s2_Re s2_Im' write(lp,*) ' Got:'//Tide_Line(1:50) if (Tide_Line(22:27).eq.'m2_amp') then write(lp,*)'ports_v.input file is in Amplitude, Phase Mode?' write(lp,*)'Expecting mode_Re mode_Im data' endif write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbd_tide)') stop '(latbd_tide)' endif !3rd line ! write(6,*)'Modes = ',('/'//Tide_Line(14+8*i:18+8*i),i=1,16) Tide_Modes_Correct=.true. do i= 1,8 Tide_Modes_Correct=Tide_Modes_Correct .and. & Tide_Line( 6+16*i:10+16*i).eq.Tide_Names(i)//'_Re' .and. & Tide_Line(14+16*i:18+16*i).eq.Tide_Names(i)//'_Im' enddo !i if (.not.Tide_Modes_Correct) then if (mnproc.eq.1) then write(lp,*) write(lp,*) 'error in latbd_tide - ports_v.input' write(lp,*) 'Tidal Modes may be in wrong order!' write(lp,*) 'Expecting: m2,s2,k1,o1,n2,p1,k2,q1' write(lp,*) ' Got: ',(Tide_Line(6+16*i:8+16*i),i=1,8) if (Tide_Line(22:27).eq.'m2_amp') then write(lp,*) & 'ports_v.input file is in Amplitude, Phase Mode!' write(lp,*) & 'Expecting mode_Re mode_Im data' endif write(lp,*) call flush(lp) endif !1st tile call xcstop('(latbd_tide)') stop '(latbd_tide)' endif !.not.Tide_Modes_Correct do i= 1,nportpts read(uoff+99,'(a)')Tide_Line if (index(Tide_Line,'Site').eq.0) then if (len_trim(Tide_Line).eq.147) then ! --- from OTPS2 read(Tide_Line,'(F9.4,F10.4,16F8.3)')dum1,dum2, & (z_R(j,i,3),z_I(j,i,3),j=1,ncon) else ! --- from OTPSnc read(Tide_Line,'(2F9.3,16F8.3)')dum1,dum2, & (z_R(j,i,3),z_I(j,i,3),j=1,ncon) endif elseif (mnproc.eq.1) then write(lp,'(a,i5,a / a)') & 'WARNING: port location',i,' treated as all zeros', & trim(Tide_Line) endif if (ldebug_latbd_t .and. mnproc.eq.1) then write(6,'(a,i4,a,8F8.3/23x,8F8.3)') & 'Port Point ',i,' z(j)= ', & (z_R(j,i,3),j=1,ncon), & (z_I(j,i,3),j=1,ncon) endif !ldebug_latbd_t enddo !i close(unit=uoff+99) if (mnproc.eq.1 .and. ldebug_latbd_t) then write(lp,'(18a)')'v tidal Modes:', & (' '//Tide_Line(14+8*i:17+8*i),i=1,16) write(lp,*)'read for',nportpts,' open boundary points' endif !ldebug_latbd_t ! return end subroutine latbd_tide ! ! !> Revision history: !> !> Mar. 2004 -- fixed bug in latbdp's speed calculation !> Nov. 2006 -- added latbdf !> Aug. 2011 -- tidal response at ports added !> Mar. 2012 -- added latbdtf !> Mar. 2012 -- replaced speed[nsew] with pspeed, so that ports can "overlap" ! Mar. 2012 -- added mnp to select mnflg for xclegt and xclput !> Apr. 2012 -- added latbd_tide !> Apr. 2012 -- replaced xclget and xclput with xciget and xciput in latbdtf !> June 2012 -- replaced xclget and xclput with xciget and xciput in latbdt !> June 2012 -- removed parameter mports, arrays are allocated at run time !> June 2012 -- fixed bnest interpolation bug, use lb[01] not ln[01] !> Mar. 2013 -- latbd_tide: added support for OTPS2 data with OTPSnc headers !> June 2013 -- added latbdtc for clammped velocity obc. !> May 2014 -- removed latbdtc !> Feb 2019 -- error stop on nports=0, because lbflag>0 !> Oct 2019 -- optionally added sports: nested velocity scale factor !> Oct 2019 -- use [uv]pnst in Flather, as in Browning&Kreiss !> Oct 2019 -- added ltrans_latbdt for port transport diagnostic !> Oct 2019 -- update pline in latbdt every npline time steps !> Oct 2019 -- smooth the Browning&Kreiss normal transport !> Oct 2019 -- npline=3 via a CPP macro