subroutine latbdt_river(n,lll)
use mod_xc ! HYCOM communication interface
use mod_cb_arrays ! HYCOM saved arrays
use mod_nc ! Netcdf output (debug)
implicit none
c
integer n,lll,il,jl,imid,jmid
c
c --- Apply lateral boundary conditions to barotropic flow field
c
c --- River version:
c --- Uses the 'Browning and Kreiss' MICOM open boundary condition.
c
c --- Note that 'speed' is in fact thref*SQRT(gH) which represents
c --- c1*thref/g in the notation of (Bleck and Sun, Open boundary conditions
c --- for MICOM). The thref/g allows for the use of pressure fields.
c
c --- The first call is made during initialization.
c
c --- Alan J. Wallcraft, NRL, July, 2001.
c --- Alexandra Bozec, COAPS, based on Remy Baraille, SHOM, Apr, 2016
c
integer, parameter :: npas=18
integer, parameter :: nchar=120
c
character*10 cfall
logical lfatal,lfatalp,rarwave, boundary_mpi,debitf
integer i,j,k,isec,ifrst,ilast,l,np,npf,npi,npl
integer ite,ios,jname,nrec,margin,nts_day,urivr
real aline(nchar)
real crs,fin,fatal,xf,xfder,dtime
character*3 char3
character*100 cline
c
integer ll1, ll2
logical ifirst
c --- Ports properties and locations
integer, save :: nports
integer, save, allocatable :: kdport(:),
& ifport(:),ilport(:),
& jfport(:),jlport(:),lnport(:)
real , save, allocatable :: xpport(:)
character*10, save, allocatable :: nameriver(:)
c --- River properties
real :: xtemp,xsaln,xdeb
real, save, allocatable :: xtemp1(:),xsaln1(:),
& xtemp2(:),xsaln2(:),
& xtime1(:),xdeb1(:),
& xtime2(:),xdeb2(:)
real, save, allocatable :: pline(:),uline(:,:),xline(:)
c --- Boundary conditions properties
real, save, allocatable :: speedw(:,:),rspedw(:,:),speede(:,:),
& rspede(:,:),speedn(:,:),rspedn(:,:),
& speeds(:,:),rspeds(:,:)
c
c --- Loop indices of ports
integer, save, allocatable :: ibeg(:), iend(:),jbeg(:), jend(:)
integer, save, allocatable :: indport(:)
integer, save, allocatable :: ifport_loc(:), ilport_loc(:),
& lnport_loc(:),
& jfport_loc(:), jlport_loc(:)
c
character*13 fmt
save fmt
data fmt / '(i4,1x,120i1)' /
c
integer lcount
save lcount
data lcount / 0 /
if (n.eq.0.and.lcount.ne.0) return
c --- Conservation of transport
c rarwave = true for invariants, false for characteristics
c debitf = true for transport conservation, false for SSH conservation
rarwave = .true.
boundary_mpi = .true.
debitf = .true.
c
lcount = lcount + 1
urivr = 1000+19
margin = 5
c --- Get time
nts_day = nint(86400.0d0/baclin) !no. time steps per day
dtime=(nstep/nts_day)+mod(nstep,nts_day)*(baclin/86400.0d0)
c
c --- The first call just initializes data structures.
c
if (lcount.eq.1) then
c
open(unit=uoff+99,file=trim(flnminp)//'rivers.list')
c
c --- 'nports' = number of boundary port sections.
call blkini(nports,'nports')
c --- Array allocations
allocate(
& nameriver(nports),
& kdport(nports),
& ifport(nports),
& ilport(nports),
& jfport(nports),
& jlport(nports),
& xpport(nports),
& lnport(nports) )
c --- Read the rivers.list
do l=1,nports
call blkinc(nameriver(l),'namriv','(a6," = ",a10)')
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(xpport(l),'xpport','(a6," =",f10.4," %")')
enddo
c
c --- Sanity check.
c
do l= 1,nports
c
lnport(l) = ilport(l)-ifport(l)+jlport(l)-jfport(l)+1
if (kdport(l).gt.2) then
if (ifport(l).ne.ilport(l)) then
if (mnproc.eq.1) then
write(lp,*)
write(lp,*) ' port #', l
write(lp,*) 'error in latbdt - port direction',
& ' and orientation are not consistent'
write(lp,*)
call flush(lp)
endif
call xcstop('(latbdt)')
stop '(latbdt)'
endif
else
if (jfport(l).ne.jlport(l)) then
if (mnproc.eq.1) then
write(lp,*)
write(lp,*) ' port #', l
write(lp,*) 'error in latbdt - port direction',
& ' and orientation are not consistent'
write(lp,*)
call flush(lp)
endif
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,*) ' port #', l
write(lp,*) 'error in latbdt - port',
& ' location is not consistent'
write(lp,*)
call flush(lp)
endif
call xcstop('(latbdt)')
stop '(latbdt)'
endif
enddo
c
close(unit=99)
c
c --- Check ports against masks,
c --- Mark the port locations on masks and print them out.
c
lfatal = .false.
do l= 1,nports
lfatalp = .false.
c
if (kdport(l).eq.4) then
c
c --- Western port
c
i = ifport(l)
do j= jfport(l),jlport(l)
if (i.lt.1 .or. i.gt.itdm-2 .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
c
elseif (kdport(l).eq.3) then
c
c --- Eastern port
c
i = ifport(l)
do j= jfport(l),jlport(l)
if (i.lt.3 .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
c
elseif (kdport(l).eq.1) then
c
c --- Northern port
c
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
c
elseif (kdport(l).eq.2) then
c
c --- Southern port
c
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
c
endif
c
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
c
c --- Local lfatal to global lfatal
c
if (lfatal) then
fatal = 1.0
else
fatal = 0.0
endif
call xcmaxr(fatal)
lfatal = fatal.gt.0.5
c
c --- Write out -iu- and -iv- arrays, if they are not too big
c --- Data are written in strips nchar points wide
c
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)
c
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
c
if (lfatal) then
write(lp,*)
write(lp,*) 'error in latbdt - bad port(s)'
write(lp,*)
call flush(lp)
call xchalt('(latbdt)')
stop '(latbdt)'
endif
c
c --- Restore iu and iv, and zero iuopn and ivopn.
c
!$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
c
c --- Modification of the storage of loop indices in case of complex boundaries
c
c --- Array allocations
allocate(
& jbeg(nports),ibeg(nports),
& jend(nports),iend(nports),
& indport(nports))
c --- Western port
jbeg(:) = 0
jend(:) = 0
ll1 = 0
do l = 1, nports
if (kdport(l).eq.4) then
i = ifport(l)
j = jfport(l)
do ll2 = 1, ll1
if ((j.ge.jbeg(ll1).or.j+lnport(l).ge.jbeg(ll1)).and.
& (j.le.jend(ll1).or.j+lnport(l).le.jend(ll1))) then
ll1 = ll1 + 1
endif
enddo
ll1 = max(1,ll1)
jbeg(ll1) = j
jend(ll1) = j+lnport(l)
indport(l) = ll1
endif
enddo
c --- Eastern port
jbeg(:) = 0
jend(:) = 0
ll1 = 0
do l = 1, nports
if (kdport(l).eq.3) then
i = ifport(l)
j = jfport(l)
do ll2 = 1, ll1
if ((j.ge.jbeg(ll1).or.j+lnport(l).ge.jbeg(ll1)).and.
& (j.le.jend(ll1).or.j+lnport(l).le.jend(ll1))) then
ll1 = ll1 + 1
endif
enddo
ll1 = max(1,ll1)
jbeg(ll1) = j
jend(ll1) = j+lnport(l)
indport(l) = ll1
endif
enddo
c --- Southern port
ibeg(:) = 0
iend(:) = 0
ll1 = 0
do l = 1, nports
if (kdport(l).eq.2) then
i = ifport(l)
j = jfport(l)
do ll2 = 1, ll1
if ((i.ge.ibeg(ll1).or.i+lnport(l).ge.ibeg(ll1)).and.
& (i.le.iend(ll1).or.i+lnport(l).le.iend(ll1))) then
ll1 = ll1 + 1
endif
enddo
ll1 = max(1,ll1)
ibeg(ll1) = i
iend(ll1) = i+lnport(l)
indport(l) = ll1
endif
enddo
c --- Northern port
ibeg(:) = 0
iend(:) = 0
ll1 = 0
do l = 1, nports
if (kdport(l).eq.1) then
i = ifport(l)
j = jfport(l)
do ll2 = 1, ll1
if ((i.ge.ibeg(ll1).or.i+lnport(l).ge.ibeg(ll1)).and.
& (i.le.iend(ll1).or.i+lnport(l).le.iend(ll1))) then
ll1 = ll1 + 1
endif
enddo
ll1 = max(1,ll1)
ibeg(ll1) = i
iend(ll1) = i+lnport(l)
indport(l) = ll1
endif
enddo
c
i = -10
do l = 1, nports
i = max(i,indport(l))
enddo
c
c --- Initialize the ports
c
c --- Array allocations
allocate (speedw(jtdm,i) , rspedw(jtdm,i) , speede(jtdm,i) ,
& rspede(jtdm,i) , speedn(itdm,i) , rspedn(itdm,i) ,
& speeds(itdm,i) , rspeds(itdm,i))
allocate(
& pline(itdm+jtdm),
& uline(itdm+jtdm,2),
& xline(itdm+jtdm))
c
do l= 1,nports
ll1 = indport(l)
if (kdport(l).eq.4) then
c
c --- Western port
c
xf= 0.
i = ifport(l)
j = jfport(l)
if (boundary_mpi) then
call xclget(pline(j) ,lnport(l),depths,i,j,0,1,0)
call xclget(rspedw(j,ll1),lnport(l),onetai,i,j,0,1,0)
call xclget(xline(j) ,lnport(l), scuy,i,j,0,1,0) !!Alex
else
if (i.ge.i0+1-nbdy.and.i.le.i0+ii+nbdy) then
do k = max(1-nbdy,j-j0), min(jj+nbdy,j-j0+lnport(l)-1)
pline(k+j0) = depths(i-i0,k)
rspedw(k+j0,ll1) = onetai(i-i0,k)
xline(k+j0) = scuy(i-i0,k) !!Alex
enddo
endif
endif
do j= jfport(l),jlport(l)
xf = xf + xline(j)
enddo
do j= jfport(l),jlport(l)
c
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
speedw(j,ll1) = rspedw(j,ll1)*pline(j)*onem
rspedw(j,ll1) = sqrt(thref/(onem*pline(j)))
iuopn(i-i0,j-j0) = 2 ! river
else
speedw(j,ll1) = onem
rspedw(j,ll1) = sqrt(thref/onem)
endif
enddo
#if defined(MPI)
call gdsum(xf,xfder) !! send xf and xfder to all processors
#else
xfder=xf
#endif
xpport(l) = 0.01*onem*xpport(l)/max(1.,xf)
c
elseif (kdport(l).eq.3) then
c
c --- Eastern port
c
xf= 0.
i = ifport(l)-1
j = jfport(l)
if (boundary_mpi) then
call xclget(pline(j) ,lnport(l),depths,i ,j,0,1,0)
call xclget(rspede(j,ll1),lnport(l),onetai,i ,j,0,1,0)
call xclget(xline(j) ,lnport(l), scuy,i+1,j,0,1,0) !!Alex
!!Alex call xclget(xline(j) ,lnport(l), scuy,i,j,0,1,0)
else
if (i.ge.i0+1-nbdy.and.i.le.i0+ii+nbdy) then
do k = max(1-nbdy,j-j0), min(jj+nbdy,j-j0+lnport(l)-1)
pline(k+j0) = depths(i-i0,k)
rspede(k+j0,ll1) = onetai(i-i0,k)
xline(k+j0) = scuy(i+1-i0,k) !!Alex
!!Alex xline(k+j0) = scuy(i-i0,k)
enddo
endif
endif
do j= jfport(l),jlport(l)
xf = xf + xline(j)
enddo
do j= jfport(l),jlport(l)
c
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
speede(j,ll1) = rspede(j,ll1)*pline(j)*onem
rspede(j,ll1) = sqrt(thref/(onem*pline(j)))
iuopn(i-i0+1,j-j0) = 2 ! river
else
speede(j,ll1) = onem
rspede(j,ll1) = sqrt(thref/onem)
endif
enddo
#if defined(MPI)
call gdsum(xf,xfder) !! send xf and xfder to all processors
#else
xfder=xf
#endif
xpport(l) = -0.01*onem*xpport(l)/max(1.,xf)
c
elseif (kdport(l).eq.1) then
c
c --- Northern port
c
xf= 0.
j = jfport(l)-1
i = ifport(l)
if (boundary_mpi) then
call xclget(pline(i) ,lnport(l),depths,i,j ,1,0,0)
call xclget(rspedn(i,ll1),lnport(l),onetai,i,j ,1,0,0)
call xclget(xline(i) ,lnport(l), scvx,i,j+1,1,0,0) !!Alex
!!Alex call xclget(xline(i) ,lnport(l), scvx,i,j,1,0,0)
else
if (j.ge.j0+1-nbdy.and.j.le.j0+jj+nbdy) then
do k = max(1-nbdy,i-i0), min(ii+nbdy,i-i0+lnport(l)-1)
pline(k+i0) = depths(k,j-j0)
rspedn(k+i0,ll1) = onetai(k,j-j0)
xline(k+i0) = scvx(k,j+1-j0) !!Alex
!!Alex xline(k+i0) = scvx(k,j-j0)
enddo
endif
endif
do i= ifport(l),ilport(l)
xf = xf + xline(i)
enddo
do i= ifport(l),ilport(l)
c
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
speedn(i,ll1) = rspedn(i,ll1)*pline(i)*onem
rspedn(i,ll1) = sqrt(thref/(onem*pline(i)))
ivopn(i-i0,j-j0+1) = 2 ! river
else
speedn(i,ll1) = onem
rspedn(i,ll1) = sqrt(thref/onem)
endif
enddo
#if defined(MPI)
call gdsum(xf,xfder) !! send xf and xfder to all processors
#else
xfder=xf
#endif
xpport(l) = -0.01*onem*xpport(l)/max(1.,xf)
c
elseif (kdport(l).eq.2) then
c
c --- Southern port
c
xf= 0.
j = jfport(l)
i = ifport(l)
if (boundary_mpi) then
call xclget(pline(i) ,lnport(l),depths,i,j,1,0,0)
call xclget(rspeds(i,ll1),lnport(l),onetai,i,j,1,0,0)
call xclget(xline(i) ,lnport(l), scvx,i,j,1,0,0) !!Alex
else
if (j.ge.j0+1-nbdy.and.j.le.j0+jj+nbdy) then
do k = max(1-nbdy,i-i0), min(ii+nbdy,i-i0+lnport(l)-1)
pline(k+i0) = depths(k,j-j0)
rspeds(k+i0,ll1) = onetai(k,j-j0)
xline(k+i0) = scvx(k,j-j0) !!Alex
enddo
endif
endif
do i= ifport(l),ilport(l)
xf = xf + xline(i)
enddo
do i= ifport(l),ilport(l)
c
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
speeds(i,ll1) = rspeds(i,ll1)*pline(i)*onem
rspeds(i,ll1) = sqrt(thref/(onem*pline(i)))
ivopn(i-i0,j-j0) = 2 ! river
else
speeds(i,ll1) = onem
rspeds(i,ll1) = sqrt(thref/onem)
endif
enddo
c
#if defined(MPI)
call gdsum(xf,xfder) !! send xf and xfder to all processors
#else
xfder=xf
#endif
xpport(l) = 0.01*onem*xpport(l)/max(1.,xf)
endif
enddo !l=1,nports
c
c -- Update halo of iuopn and ivopn
util1(:,:) = iuopn(:,:)
util2(:,:) = ivopn(:,:)
call xctilr(util1,1,1, nbdy,nbdy, halo_us)
call xctilr(util2,1,1, nbdy,nbdy, halo_vs)
iuopn(:,:) = util1(:,:)
ivopn(:,:) = util2(:,:)
c
c --- Modification of loop indices for ports
c
c --- Arrays allocation
allocate(
& ifport_loc(nports),ilport_loc(nports),
& jfport_loc(nports),jlport_loc(nports),
& lnport_loc(nports))
allocate(
& xtemp1(nports),xtemp2(nports),
& xsaln1(nports),xsaln2(nports),
& xdeb1(nports), xdeb2(nports),
& xtime1(nports),xtime2(nports))
do l = 1,nports
c --- Western/Eastern port
if (kdport(l).eq.4.or.kdport(l).eq.3) then
i = ifport(l)
ifirst = .false.
if (i.ge.1+i0-nbdy.and.i.le.ii+i0+nbdy) then
do j = jfport(l),jlport(l)
if (j.ge.1+j0-nbdy.and.j.le.jj+j0+nbdy) then
if (.not.ifirst) then
ifrst = j
ifirst = .true.
endif
ilast = j
endif
enddo
endif
if (ifirst) then
jfport_loc(l) = ifrst
jlport_loc(l) = ilast
lnport_loc(l) = ilast-ifrst+1
else
lnport_loc(l) = 0
jfport_loc(l) = 10000
jlport_loc(l) = -10000
endif
endif
c --- Southern/Northern port
if (kdport(l).eq.1.or.kdport(l).eq.2) then
j = jfport(l)
ifirst = .false.
if (j.ge.1+j0-nbdy.and.j.le.jj+j0+nbdy) then
do i = ifport(l),ilport(l)
if (i.ge.1+i0-nbdy.and.i.le.ii+i0+nbdy) then
if (.not.ifirst) then
ifrst = i
ifirst = .true.
endif
ilast = i
endif
enddo
endif
if (ifirst) then
ifport_loc(l) = ifrst
ilport_loc(l) = ilast
lnport_loc(l) = ilast-ifrst+1
else
lnport_loc(l) = 0
ifport_loc(l) = 10000
ilport_loc(l) = -10000
endif
endif
enddo
c
c --- End of initialization
call xcsync(flush_lp)
xtime1(:) = 0.0
xtime2(:) = 0.0
return
endif ! lcount.eq.1
c
c --- 'wellposed' treatment of pressure and velocity fields.
c --- Alternate order of ports in case corners are open.
c
if (mod(lll,2).eq.1) then
npf = 1
npl = nports
npi = 1
else
npf = nports
npl = 1
npi = -1
endif
c
do l = npf, npl, npi
c
ll1 = indport(l)
c
if (lll.eq.1.and. (xtime1(l).gt.dtime .or.
& dtime .gt.xtime2(l)) ) then
open(unit=urivr,file=trim(nameriver(l))//'.r')
c
j = 0
do nrec = 1, 9999999
read(urivr,'(a)',iostat=ios) cline
if (j.eq.0) then
if (cline(1:3).eq.'---') j=1
else
if (ios.ne.0) then
if (mnproc.eq.1) write(lp,*)'missing data '//
& 'for river:',nameriver(l)
call xcstop('(latbdt)')
stop '(latbdt)'
endif
read (cline(1:),*) xtime2(l), xdeb2(l),
& xtemp2(l),xsaln2(l)
if (j.eq.1.and.
& dtime-0.5*baclin/86400..le.xtime2(l)) then
if (mnproc.eq.1) write(lp,*)'missing data '//
& 'for river:',nameriver(l)
call xcstop('(latbdt)')
stop '(latbdt)'
endif
if (dtime-0.5*baclin/86400..le.xtime2(l)) then
close(urivr)
goto 200
endif
xtime1(l) = xtime2(l)
xdeb1(l) = xdeb2(l)
xtemp1(l) = xtemp2(l)
xsaln1(l) = xsaln2(l)
j = j + 1
endif
enddo
200 continue
endif ! lll.eq.1 .and. xtime
c --- Calculate the transport of the river for each point of the port
xdeb = ( xdeb2(l)+( dtime-xtime2(l)) *(xdeb1(l)-xdeb2(l))/
& (xtime1(l)-xtime2(l)) )* xpport(l)
if (lll.eq.1) then
c --- Calculate the T and S of the river for each point of the port
xtemp = xtemp2(l) + (dtime-xtime2(l))*(xtemp1(l)-xtemp2(l))
& /(xtime1(l)-xtime2(l))
xsaln = xsaln2(l) + (dtime-xtime2(l))*(xsaln1(l)-xsaln2(l))
& /(xtime1(l)-xtime2(l))
c
c --- Western/Eastern port
c --- Fill triver and sriver arrays for tsadvc (advem_fct2c)
c
if (kdport(l).eq.4.or.kdport(l).eq.3) then
i = ifport(l)
do j = jfport(l),jlport(l)
! pline(j) = xtemp
if (xtemp.gt.-99.9) then
pline(j) = xtemp
else ! no temp data for river
if (kdport(l).eq.4) pline(j) = temp(i ,j,1,n)
if (kdport(l).eq.3) pline(j) = temp(i-1,j,1,n)
endif
enddo
j = jfport(l)
call xclput(pline(j),lnport(l),
& triver(1-nbdy,1-nbdy), i, j, 0, 1)
do j = jfport(l),jlport(l)
! pline(j) = xsaln
if (xsaln.gt.-99.9) then
pline(j) = xsaln
else ! no saln data for river
if (kdport(l).eq.4) pline(j) = saln(i ,j,1,n)
if (kdport(l).eq.3) pline(j) = saln(i-1,j,1,n)
endif
enddo
j = jfport(l)
call xclput(pline(j),lnport(l),
& sriver(1-nbdy,1-nbdy), i, j, 0, 1)
endif
c
c --- Southern/Northern port
c
if (kdport(l).eq.1.or.kdport(l).eq.2) then
j = jfport(l)
do i = ifport(l),ilport(l)
! pline(i) = xtemp
if (xtemp.gt.-99.9) then
pline(i) = xtemp
else ! no temp data for river
if (kdport(l).eq.1) pline(i) = temp(i,j-1,1,n)
if (kdport(l).eq.2) pline(i) = temp(i,j ,1,n)
endif
enddo
i = ifport(l)
call xclput(pline(i),lnport(l),
& triver(1-nbdy,1-nbdy), i, j, 1, 0)
do i = ifport(l),ilport(l)
! pline(i) = xsaln
if (xsaln.gt.-99.9) then
pline(i) = xsaln
else ! no saln data for river
if (kdport(l).eq.1) pline(i) = saln(i,j-1,1,n)
if (kdport(l).eq.2) pline(i) = saln(i,j ,1,n)
endif
enddo
i = ifport(l)
call xclput(pline(i),lnport(l),
& sriver(1-nbdy,1-nbdy), i, j, 1, 0)
endif
endif ! lll.eq.1
c
if (kdport(l).eq.4) then
c
c --- Western port
c
i = ifport(l)
j = jfport(l)
if (boundary_mpi) then
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)
call xclget(pline(j), lnport(l),
& pbavg(1-nbdy,1-nbdy,n), i, j, 0,1, 0)
else
if (i.ge.i0+1-margin.and.i.le.i0+ii+margin) then
do k=max(1 -margin,j-j0),
& min(jj+margin,j-j0+lnport(l)-1)
pline(k+j0) = pbavg(i-i0,k,n)
enddo
endif
if (i.ge.i0-margin.and.i+1.le.i0+ii+margin) then
do k=max(1 -margin,j-j0),
& min(jj+margin,j-j0+lnport(l)-1)
uline(k+j0,1) = ubavg(i+1-i0,k,n)
enddo
endif
if (i+1.ge.i0-margin.and.i+2.le.i0+ii+margin) then
do k=max(1 -margin,j-j0),
& min(jj+margin,j-j0+lnport(l)-1)
uline(k+j0,2) = ubavg(i+2-i0,k,n)
enddo
endif
endif ! boundary_mpi
if (rarwave) then
do j= max(j0+1 -margin,jfport_loc(l)),
& min(jj+j0+margin,jlport_loc(l))
fin = uline(j,1)
& - 2.*sqrt(max(1.e-6,(speedw(j,ll1)+pline(j))*thref))
do ite = 1, npas
xf = 0.25*uline(j,1)*(fin-uline(j,1))*
& (fin-uline(j,1))/thref - xdeb
xfder = 0.25*(fin- uline(j,1))*
& (fin-3.*uline(j,1))/thref
uline(j,1) = uline(j,1) - xf/xfder
enddo
pline(j) = 0.25*(fin-uline(j,1))
& *(fin-uline(j,1))/thref - speedw(j,ll1)
c
c --- Solver convergence problem
c
if (debitf.and.abs(xf/xfder).ge.1e-3) then
pline(j) = -speedw(j,ll1)
uline(j,1) = 0.
if (mnproc.eq.1) then
write(6,*) 'Solver convergence pb in ',i,',',j
write(6,*)'xf/xfder = ', xf/xfder, 'xf=',xf,
& 'xfder=',xfder
endif
endif
enddo
else
do j= max(j0+1 -margin,jfport_loc(l)),
& min(j0+jj+margin,jlport_loc(l))
fin = 1.5*uline(j,1)-0.5*uline(j,2)
& -rspedw(j,ll1)*pline( j)
do ite = 1, npas
xf = (speedw(j,ll1)-(fin-uline(j,1))/rspedw(j,ll1))
& *uline(j,1) - xdeb
xfder = speedw(j,ll1)-(fin-2.*uline(j,1))
& /rspedw(j,ll1)
uline(j,1) = uline(j,1) - xf/xfder
enddo
pline(j) = -(fin - uline(j,1))/rspedw(j,ll1)
c
c --- Solver convergence problem
c
if (debitf.and.abs(xf/xfder).ge.1e-3) then
pline(j) = -speedw(j,ll1)
uline(j,1) = 0.
if (mnproc.eq.1) then
write(6,*) 'Solver convergence pb in ',i,',',j
write(6,*)'xf/xfder = ', xf/xfder, 'xf=',xf,
& 'xfder=',xfder
endif
endif
enddo
endif ! rarwave
j = jfport(l)
if (boundary_mpi) then
call xclput(pline(j), lnport(l),
& pbsavu(1-nbdy,1-nbdy ), i, j, 0,1)
call xclput(uline(j,1),lnport(l),
& ubavg (1-nbdy,1-nbdy,n), i, j, 0,1) ! normal
else
if (i.ge.i0+1-margin.and.i.le.i0+ii+margin) then
do k=max(1-margin,j-j0),min(jj+margin,j-j0+lnport(l)-1)
pbsavu(i-i0,k ) = pline(k+j0)
ubavg (i-i0,k,n) = uline(k+j0,1)
enddo
endif
endif ! boundary_mpi
c
c
elseif (kdport(l).eq.3) then
c
c --- Eastern port
c
i = ifport(l)-1
j = jfport(l)
if (boundary_mpi) then
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)
call xclget(pline(j), lnport(l),
& pbavg(1-nbdy,1-nbdy,n), i, j, 0,1, 0)
else
if (i.ge.i0+1-margin.and.i.le.i0+ii+margin) then
do k = max(1 -margin,j-j0),
& min(jj+margin,j-j0+lnport(l)-1)
pline(k+j0) = pbavg(i-i0,k,n)
uline(k+j0,1) = ubavg(i-i0,k,n)
enddo
endif
if (i.ge.i0+2-margin.and.i-1.le.i0+ii+margin) then
do k = max(1 -margin,j-j0),
& min(jj+margin,j-j0+lnport(l)-1)
uline(k+j0,2) = ubavg(i-1-i0,k,n)
enddo
endif
endif ! boundary_mpi
if (rarwave) then
do j= max(j0+1 -margin,jfport_loc(l)),
& min(j0+jj+margin,jlport_loc(l))
fin = uline(j,1)
& + 2.*sqrt(max(1.e-6,(speede(j,ll1)+pline(j))*thref))
do ite = 1, npas
xf = 0.25*uline(j,1)*(fin-uline(j,1))*
& (fin-uline(j,1))/thref - xdeb
xfder = 0.25*(fin- uline(j,1))*
& (fin-3.*uline(j,1))/thref
uline(j,1) = uline(j,1) - xf/xfder
enddo
pline(j) = 0.25*(fin-uline(j,1))
& *(fin-uline(j,1))/thref - speede(j,ll1)
c
c --- Solver convergence problem
c
if (debitf.and.abs(xf/xfder).ge.1e-3) then
if (mnproc.eq.1) then
write(6,*) 'Solver convergence pb in ',i,',',j
write(6,*)'xf/xfder = ', xf/xfder, 'xf=',xf,
& 'xfder=',xfder
endif
pline(j) = -speede(j,ll1)
uline(j,1) = 0.
endif
enddo
else
do j= max(j0+1 -margin,jfport_loc(l)),
& min(j0+jj+margin,jlport_loc(l))
fin = 1.5*uline(j,1)-0.5*uline(j,2)
& +rspede(j,ll1)*pline( j)
do ite = 1, npas
xf = (speede(j,ll1)+(fin-uline(j,1))/rspede(j,ll1))
& *uline(j,1) - xdeb
xfder = speede(j,ll1)+(fin-2.*uline(j,1))
& /rspede(j,ll1)
uline(j,1) = uline(j,1) - xf/xfder
enddo
pline(j) = (fin - uline(j,1))/rspede(j,ll1)
c
c --- Solver convergence problem
c
if (debitf.and.abs(xf/xfder).ge.1e-3) then
pline(j) = -speede(j,ll1)
uline(j,1) = 0.
if (mnproc.eq.1) then
write(6,*) 'Solver convergence pb in ',i,',',j
write(6,*)'xf/xfder = ', xf/xfder,
& 'xf =',xf,
& 'xfder =',xfder
endif
endif
enddo
endif ! rarwave
j = jfport(l)
if (boundary_mpi) then
call xclput(pline(j), lnport(l),
& pbsavu(1-nbdy,1-nbdy ), i+1,j,0,1)
call xclput(uline(j,1),lnport(l),
& ubavg (1-nbdy,1-nbdy,n), i+1,j,0,1) ! normal
else
if (i.ge.i0-margin.and.i+1.le.i0+ii+margin) then
do k=max(1 -margin,j-j0),
& min(jj+margin,j-j0+lnport(l)-1)
pbsavu(i+1-i0,k ) = pline(k+j0)
ubavg (i+1-i0,k,n) = uline(k+j0,1)
enddo
endif
endif ! boundary_mpi
c
elseif (kdport(l).eq.1) then
c
c --- Northern port
c
j = jfport(l)-1
i = ifport(l)
if (boundary_mpi) then
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)
call xclget(pline(i), lnport(l),
& pbavg(1-nbdy,1-nbdy,n), i, j, 1,0, 0)
else
if (j.ge.j0+1-margin.and.j.le.j0+jj+margin) then
do k = max(1 -margin,i-i0),
& min(ii+margin,i-i0+lnport(l)-1)
pline(k+i0) = pbavg(k,j-j0,n)
uline(k+i0,1) = vbavg(k,j-j0,n)
enddo
endif
if (j.ge.j0+2-margin.and.j-1.le.j0+jj+margin) then
do k = max(1 -margin,i-i0),
& min(ii+margin,i-i0+lnport(l)-1)
uline(k+i0,2) = vbavg(k,j-1-j0,n)
enddo
endif
endif ! boundary_mpi
if (rarwave) then
do i= max(i0+1 -margin,ifport_loc(l)),
& min(i0+ii+margin,ilport_loc(l))
fin = uline(i,1)
& + 2.*sqrt(max(1.e-6,(speedn(i,ll1)+pline(i))*thref))
do ite = 1, npas
xf = 0.25*uline(i,1)*(fin-uline(i,1))*
& (fin-uline(i,1))/thref -xdeb
xfder = 0.25*(fin-uline(i,1))*
& (fin-3.*uline(i,1))/thref
uline(i,1) = uline(i,1) - xf/xfder
enddo
pline(i) = 0.25*(fin-uline(i,1))*
& (fin-uline(i,1))/thref - speedn(i,ll1)
c
c --- Solver convergence problem
c
if (debitf.and.abs(xf/xfder).ge.1e-3) then
pline(i) = -speedn(i,ll1)
uline(i,1) = 0.
if (mnproc.eq.1) then
write(6,*) 'Solver convergence pb in ',i,',',j
write(6,*) 'xf/xfder = ', xf/xfder, 'xf=',xf,
& 'xfder=',xfder
endif
endif
enddo
else
do i= max(i0+1 -margin,ifport_loc(l)),
& min(i0+ii+margin,ilport_loc(l))
fin = 1.5*uline(i,1)-0.5*uline(i,2)
& +rspedn(i,ll1)*pline( i)
do ite = 1, npas
xf = (speedn(i,ll1)+(fin-uline(i,1))/rspedn(i,ll1))
& *uline(i,1) - xdeb
xfder = speedn(i,ll1)+(fin-2.*uline(i,1))
& /rspedn(i,ll1)
uline(i,1) = uline(i,1) - xf/xfder
enddo
pline(i) = (fin - uline(i,1))/rspedn(i,ll1)
c
c --- Solver convergence problem
c
if (debitf.and.abs(xf/xfder).ge.1e-3) then
pline(i) = -speedn(i,ll1)
uline(i,1) = 0.
if (mnproc.eq.1) then
write(6,*) 'Solver convergence pb in ',i,',',j
write(6,*)'xf/xfder = ', xf/xfder, 'xf=',xf,
& 'xfder=',xfder
endif
endif
enddo
endif ! rarwave
i = ifport(l)
if (boundary_mpi) then
call xclput(pline(i), lnport(l),
& pbsavv(1-nbdy,1-nbdy ), i, j+1,1,0)
call xclput(uline(i,1),lnport(l),
& vbavg (1-nbdy,1-nbdy,n), i, j+1,1,0) ! normal
else
if (j.ge.j0-margin.and.j+1.le.j0+jj+margin) then
do k=max(1 -margin,i-i0),
& min(ii+margin,i-i0+lnport(l)-1)
pbsavv(k,j+1-j0 ) = pline(k+i0)
vbavg (k,j+1-j0,n) = uline(k+i0,1)
enddo
endif
endif ! boundary_mpi
elseif (kdport(l).eq.2) then
c
c --- Southern port
c
j = jfport(l)
i = ifport(l)
if (boundary_mpi) then
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)
call xclget(pline(i), lnport(l),
& pbavg(1-nbdy,1-nbdy,n), i, j, 1,0, 0)
else
if (j.ge.j0+1-margin.and.j.le.j0+jj+margin) then
do k = max(1 -margin,i-i0),
& min(ii+margin,i-i0+lnport(l)-1)
pline(k+i0) = pbavg(k,j-j0,n)
enddo
endif
if (j.ge.j0-margin.and.j+1.le.j0+jj+margin) then
do k = max(1 -margin,i-i0),
& min(ii+margin,i-i0+lnport(l)-1)
uline(k+i0,1) = vbavg(k,j+1-j0,n)
enddo
endif
if (j+1.ge.j0-margin.and.j+2.le.j0+jj+margin) then
do k = max(1 -margin,i-i0),
& min(ii+margin,i-i0+lnport(l)-1)
uline(k+i0,2) = vbavg(k,j+2-j0,n)
enddo
endif
endif ! boundary_mpi
if (rarwave) then
do i= max(i0+1 -margin,ifport_loc(l)),
& min(i0+ii+margin,ilport_loc(l))
fin = uline(i,1)
& - 2.*sqrt(max(1.e-6,(speeds(i,ll1)+pline(i))*thref))
do ite = 1, npas
xf = 0.25*uline(i,1)*(fin-uline(i,1))*
& (fin-uline(i,1))/thref - xdeb
xfder = 0.25*(fin-uline(i,1))*
& (fin-3.*uline(i,1))/thref
uline(i,1) = uline(i,1) - xf/xfder
enddo
pline(i) = 0.25*(fin-uline(i,1))*
& (fin-uline(i,1))/thref - speeds(i,ll1)
c
c --- Solver convergence problem
c
if (debitf.and.abs(xf/xfder).ge.1e-3) then
pline(i) = -speeds(i,ll1)
uline(i,1) = 0.
if (mnproc.eq.1) then
write(6,*) 'Solver convergence pb in ',i,',',j
write(6,*)'xf/xfder = ', xf/xfder, 'xf=',xf,
& 'xfder=',xfder
endif
endif
enddo
else
do i= max(i0+1 -margin,ifport_loc(l)),
& min(i0+ii+margin,ilport_loc(l))
fin = 1.5*uline(i,1)-0.5*uline(i,2)
& -rspeds(i,ll1)*pline( i)
do ite = 1, npas
xf = (speeds(i,ll1)-(fin-uline(i,1))/rspeds(i,ll1))
& *uline(i,1) - xdeb
xfder = speeds(i,ll1)-(fin-2.*uline(i,1))
& /rspeds(i,ll1)
uline(i,1) = uline(i,1) - xf/xfder
enddo
pline(i) = -(fin - uline(i,1))/rspeds(i,ll1)
c
c --- Solver convergence problem
c
if (debitf.and.abs(xf/xfder).ge.1e-3) then
pline(i) = -speeds(i,ll1)
uline(i,1) = 0.
if (mnproc.eq.1) then
write(6,*) 'Solver convergence pb in ',i,',',j
write(6,*)'xf/xfder = ', xf/xfder, 'xf=',xf,
& 'xfder=',xfder
endif
endif
enddo
endif ! rarwave
i = ifport(l)
if (boundary_mpi) then
call xclput(pline(i), lnport(l),
& pbsavv(1-nbdy,1-nbdy), i, j, 1, 0)
call xclput(uline(i,1),lnport(l),
& vbavg (1-nbdy,1-nbdy,n), i, j, 1, 0) ! normal
else
if (j.ge.j0+1-margin.and.j.le.j0+jj+margin) then
do k=max(1 -margin,i-i0),
& min(ii+margin,i-i0+lnport(l)-1)
pbsavv(k,j-j0 ) = pline(k+i0)
vbavg (k,j-j0,n) = uline(k+i0,1)
enddo
endif
endif ! boundary_mpi
c
endif ! kdports
c
enddo !l=1,nports
c
return
end