! bounds_forcing.f
! spcify variable boundary conditions, atmospheric forcing, restoring
!_______________________________________________________________________
subroutine bcond(idx)
! apply open boundary conditions
! closed boundary conditions are automatically enabled through
! specification of the masks, dum, dvm and fsm, in which case the open
! boundary conditions, included below, will be overwritten
use setvars
implicit none
integer idx
integer i,j,k
real ga,u1,wm !RMY: as_is(ecoast+seamt),+gae,utide,etide(stide)
if(idx.eq.1) then
! external (2-D) elevation boundary conditions
do j=1,jm
if(n_west.eq.-1) elf(1,j)=elf(2,j)
if(n_east.eq.-1) elf(im,j)=elf(imm1,j)
end do
do i=1,im
if(n_south.eq.-1) elf(i,1)=elf(i,2)
if(n_north.eq.-1) elf(i,jm)=elf(i,jmm1)
end do
do j=1,jm
do i=1,im
elf(i,j)=elf(i,j)*fsm(i,j)
end do
end do
return
else if(idx.eq.2) then
! external (2-D) velocity boundary conditions
do j=2,jmm1 !RMY: as_is(ecoast+seamt),j=1,jm(stide)
! west
if(n_west.eq.-1) then
uaf(2,j)=uabw(j)-rfw*sqrt(grav/h(2,j))*(el(2,j)-elw(j))
uaf(2,j)=ramp*uaf(2,j)
uaf(1,j)=uaf(2,j)
vaf(1,j)=0.e0
end if
! east
if(n_east.eq.-1) then
uaf(im,j)=uabe(j) !RMY: as_is(ecoast+seamt),u/etide(stide)
$ +rfe*sqrt(grav/h(imm1,j))*(el(imm1,j)-ele(j))
uaf(im,j)=ramp*uaf(im,j)
vaf(im,j)=0.e0 !RMY: as_is(ecoast+seamt),nonzero_gae(stide)
end if
end do
do i=2,imm1 !RMY: as_is(ecoast+seamt),i=1,im(stide)
! south
if(n_south.eq.-1) then
vaf(i,2)=vabs(i) !RMY: as_is(ecoast),uncomment(seamt+stide)
! $ -rfs*sqrt(grav/h(i,2))*(el(i,2)-els(i))
vaf(i,2)=ramp*vaf(i,2)
vaf(i,1)=vaf(i,2)
uaf(i,1)=0.e0
end if
! north
if(n_north.eq.-1) then
vaf(i,jm)=vabn(i)
$ +rfn*sqrt(grav/h(i,jmm1))*(el(i,jmm1)-eln(i))
vaf(i,jm)=ramp*vaf(i,jm)
uaf(i,jm)=0.e0
end if
end do
do j=1,jm
do i=1,im
uaf(i,j)=uaf(i,j)*dum(i,j)
vaf(i,j)=vaf(i,j)*dvm(i,j)
end do
end do
return
else if(idx.eq.3) then
! internal (3-D) velocity boundary conditions
! radiation conditions !RMY: as_is(ecoast),x_to_next_!RMY(seamt+stide)
! smoothing is used in the direction tangential to the boundaries
do k=1,kbm1
do j=2,jmm1
! east
if(n_east.eq.-1) then
ga=sqrt(h(im,j)/hmax)
uf(im,j,k)=ga*(.25e0*u(imm1,j-1,k)+.5e0*u(imm1,j,k)
$ +.25e0*u(imm1,j+1,k))
$ +(1.e0-ga)*(.25e0*u(im,j-1,k)+.5e0*u(im,j,k)
$ +.25e0*u(im,j+1,k))
vf(im,j,k)=0.e0
end if
! west
if(n_west.eq.-1) then
ga=sqrt(h(1,j)/hmax)
uf(2,j,k)=ga*(.25e0*u(3,j-1,k)+.5e0*u(3,j,k)
$ +.25e0*u(3,j+1,k))
$ +(1.e0-ga)*(.25e0*u(2,j-1,k)+.5e0*u(2,j,k)
$ +.25e0*u(2,j+1,k))
uf(1,j,k)=uf(2,j,k)
vf(1,j,k)=0.e0
end if
end do
end do
do k=1,kbm1
do i=2,imm1
! south
if(n_south.eq.-1) then
ga=sqrt(h(i,1)/hmax)
vf(i,2,k)=ga*(.25e0*v(i-1,3,k)+.5e0*v(i,3,k)
$ +.25e0*v(i+1,3,k))
$ +(1.e0-ga)*(.25e0*v(i-1,2,k)+.5e0*v(i,2,k)
$ +.25e0*v(i+1,2,k))
vf(i,1,k)=vf(i,2,k)
uf(i,1,k)=0.e0
end if
! north
if(n_north.eq.-1) then
ga=sqrt(h(i,jm)/hmax)
vf(i,jm,k)=ga*(.25e0*v(i-1,jmm1,k)+.5e0*v(i,jmm1,k)
$ +.25e0*v(i+1,jmm1,k))
$ +(1.e0-ga)*(.25e0*v(i-1,jm,k)+.5e0*v(i,jm,k)
$ +.25e0*v(i+1,jm,k))
uf(i,jm,k)=0.e0
end if
end do
end do
!RMY: as_is(ecoast),x_above_to_previous_!RMY(seamt+stide)
do k=1,kbm1
do j=1,jm
do i=1,im
uf(i,j,k)=uf(i,j,k)*dum(i,j)
vf(i,j,k)=vf(i,j,k)*dvm(i,j)
end do
end do
end do
return
else if(idx.eq.4) then
! temperature and salinity boundary conditions (using uf and vf,
! respectively)
do k=1,kbm1
do j=1,jm
! east
if(n_east.eq.-1) then
u1=2.e0*u(im,j,k)*dti/(dx(im,j)+dx(imm1,j))
if(u1.le.0.e0) then
uf(im,j,k)=t(im,j,k)-u1*(tbe(j,k)-t(im,j,k))
vf(im,j,k)=s(im,j,k)-u1*(sbe(j,k)-s(im,j,k))
else
uf(im,j,k)=t(im,j,k)-u1*(t(im,j,k)-t(imm1,j,k))
vf(im,j,k)=s(im,j,k)-u1*(s(im,j,k)-s(imm1,j,k))
if(k.ne.1.and.k.ne.kbm1) then
wm=.5e0*(w(imm1,j,k)+w(imm1,j,k+1))*dti
$ /((zz(k-1)-zz(k+1))*dt(imm1,j))
uf(im,j,k)=uf(im,j,k)-wm*(t(imm1,j,k-1)-t(imm1,j,k+1))
vf(im,j,k)=vf(im,j,k)-wm*(s(imm1,j,k-1)-s(imm1,j,k+1))
endif
end if
end if
! west
if(n_west.eq.-1) then
u1=2.e0*u(2,j,k)*dti/(dx(1,j)+dx(2,j))
if(u1.ge.0.e0) then
uf(1,j,k)=t(1,j,k)-u1*(t(1,j,k)-tbw(j,k))
vf(1,j,k)=s(1,j,k)-u1*(s(1,j,k)-sbw(j,k))
else
uf(1,j,k)=t(1,j,k)-u1*(t(2,j,k)-t(1,j,k))
vf(1,j,k)=s(1,j,k)-u1*(s(2,j,k)-s(1,j,k))
if(k.ne.1.and.k.ne.kbm1) then
wm=.5e0*(w(2,j,k)+w(2,j,k+1))*dti
$ /((zz(k-1)-zz(k+1))*dt(2,j))
uf(1,j,k)=uf(1,j,k)-wm*(t(2,j,k-1)-t(2,j,k+1))
vf(1,j,k)=vf(1,j,k)-wm*(s(2,j,k-1)-s(2,j,k+1))
end if
end if
end if
end do
do i=1,im
! south
if(n_south.eq.-1) then
u1=2.e0*v(i,2,k)*dti/(dy(i,1)+dy(i,2))
if(u1.ge.0.e0) then
uf(i,1,k)=t(i,1,k)-u1*(t(i,1,k)-tbs(i,k))
vf(i,1,k)=s(i,1,k)-u1*(s(i,1,k)-sbs(i,k))
else
uf(i,1,k)=t(i,1,k)-u1*(t(i,2,k)-t(i,1,k))
vf(i,1,k)=s(i,1,k)-u1*(s(i,2,k)-s(i,1,k))
if(k.ne.1.and.k.ne.kbm1) then
wm=.5e0*(w(i,2,k)+w(i,2,k+1))*dti
$ /((zz(k-1)-zz(k+1))*dt(i,2))
uf(i,1,k)=uf(i,1,k)-wm*(t(i,2,k-1)-t(i,2,k+1))
vf(i,1,k)=vf(i,1,k)-wm*(s(i,2,k-1)-s(i,2,k+1))
end if
end if
end if
! north
if(n_north.eq.-1) then
u1=2.e0*v(i,jm,k)*dti/(dy(i,jm)+dy(i,jmm1))
if(u1.le.0.e0) then
uf(i,jm,k)=t(i,jm,k)-u1*(tbn(i,k)-t(i,jm,k))
vf(i,jm,k)=s(i,jm,k)-u1*(sbn(i,k)-s(i,jm,k))
else
uf(i,jm,k)=t(i,jm,k)-u1*(t(i,jm,k)-t(i,jmm1,k))
vf(i,jm,k)=s(i,jm,k)-u1*(s(i,jm,k)-s(i,jmm1,k))
if(k.ne.1.and.k.ne.kbm1) then
wm=.5e0*(w(i,jmm1,k)+w(i,jmm1,k+1))*dti
$ /((zz(k-1)-zz(k+1))*dt(i,jmm1))
uf(i,jm,k)=uf(i,jm,k)-wm*(t(i,jmm1,k-1)-t(i,jmm1,k+1))
vf(i,jm,k)=vf(i,jm,k)-wm*(s(i,jmm1,k-1)-s(i,jmm1,k+1))
end if
end if
end if
end do
end do
do k=1,kbm1
do j=1,jm
do i=1,im
uf(i,j,k)=uf(i,j,k)*fsm(i,j)
vf(i,j,k)=vf(i,j,k)*fsm(i,j)
end do
end do
end do
return
else if(idx.eq.5) then
! vertical velocity boundary conditions
do k=1,kbm1
do j=1,jm
do i=1,im
w(i,j,k)=w(i,j,k)*fsm(i,j)
end do
end do
end do
return
else if(idx.eq.6) then
! q2 and q2l boundary conditions
do k=1,kb
do j=1,jm
! west
if(n_west.eq.-1) then
u1=2.e0*u(2,j,k)*dti/(dx(1,j)+dx(2,j))
if(u1.ge.0.e0) then
uf(1,j,k)=q2(1,j,k)-u1*(q2(1,j,k)-small)
vf(1,j,k)=q2l(1,j,k)-u1*(q2l(1,j,k)-small)
else
uf(1,j,k)=q2(1,j,k)-u1*(q2(2,j,k)-q2(1,j,k))
vf(1,j,k)=q2l(1,j,k)-u1*(q2l(2,j,k)-q2l(1,j,k))
end if
end if
! east
if(n_east.eq.-1) then
u1=2.e0*u(im,j,k)*dti/(dx(im,j)+dx(imm1,j))
if(u1.le.0.e0) then
uf(im,j,k)=q2(im,j,k)-u1*(small-q2(im,j,k))
vf(im,j,k)=q2l(im,j,k)-u1*(small-q2l(im,j,k))
else
uf(im,j,k)=q2(im,j,k)-u1*(q2(im,j,k)-q2(imm1,j,k))
vf(im,j,k)=q2l(im,j,k)-u1*(q2l(im,j,k)-q2l(imm1,j,k))
end if
end if
end do
do i=1,im
! south
if(n_south.eq.-1) then
u1=2.e0*v(i,2,k)*dti/(dy(i,1)+dy(i,2))
if(u1.ge.0.e0) then
uf(i,1,k)=q2(i,1,k)-u1*(q2(i,1,k)-small)
vf(i,1,k)=q2l(i,1,k)-u1*(q2l(i,1,k)-small)
else
uf(i,1,k)=q2(i,1,k)-u1*(q2(i,2,k)-q2(i,1,k))
vf(i,1,k)=q2l(i,1,k)-u1*(q2l(i,2,k)-q2l(i,1,k))
end if
end if
! north
if(n_north.eq.-1) then
u1=2.e0*v(i,jm,k)*dti/(dy(i,jm)+dy(i,jmm1))
if(u1.le.0.e0) then
uf(i,jm,k)=q2(i,jm,k)-u1*(small-q2(i,jm,k))
vf(i,jm,k)=q2l(i,jm,k)-u1*(small-q2l(i,jm,k))
else
uf(i,jm,k)=q2(i,jm,k)-u1*(q2(i,jm,k)-q2(i,jmm1,k))
vf(i,jm,k)=q2l(i,jm,k)-u1*(q2l(i,jm,k)-q2l(i,jmm1,k))
end if
end if
end do
end do
do k=1,kb
do j=1,jm
do i=1,im
uf(i,j,k)=uf(i,j,k)*fsm(i,j)+1.e-10
vf(i,j,k)=vf(i,j,k)*fsm(i,j)+1.e-10
end do
end do
end do
return
else if(idx.eq.7) then !RMY: new option for idx = 7
!RMY: t, s, u, and v boundary conditions for an a-grid
do k=1,kbm1
do j=1,jm
! west
if(n_west.eq.-1) then
uf(1,j,k)=uf(2,j,k)
vf(1,j,k)=vf(2,j,k)
end if
! east
if(n_east.eq.-1) then
uf(im,j,k)=uf(im-1,j,k)
vf(im,j,k)=vf(im-1,j,k)
end if
end do
do i=1,im
! south
if(n_south.eq.-1) then
uf(i,1,k)=uf(i,2,k)
vf(i,1,k)=vf(i,2,k)
end if
! north
if(n_north.eq.-1) then
uf(i,jm,k)=uf(i,jm-1,k)
vf(i,jm,k)=vf(i,jm-1,k)
end if
end do
end do
do k=1,kbm1
do j=1,jm
do i=1,im
uf(i,j,k)=uf(i,j,k)*fsm(i,j)
vf(i,j,k)=vf(i,j,k)*fsm(i,j)
end do
end do
end do
return
end if !RMY: end of idx = 7 option
end
!RMY: as_is(ecoast),end_of_file(seamt),as_is_to_next_!RMY(stide)
!_______________________________________________________________________
subroutine bcondorl(idx)
! this is an optional subroutine replacing bcond and using Orlanski's
! scheme (J. Comp. Phys. 21, 251-269, 1976), specialized for the
! seamount problem
use setvars
implicit none
integer idx
real cl,denom
real ar,eps
integer i,j,k
if(idx.eq.1) then
! external (2-D) elevation boundary conditions
do j=1,jm
if(n_west.eq.-1) elf(1,j)=elf(2,j)
if(n_east.eq.-1) elf(im,j)=elf(imm1,j)
end do
do j=1,jm
do i=1,im
elf(i,j)=elf(i,j)*fsm(i,j)
end do
end do
return
else if(idx.eq.2) then
! external (2-D) velocity boundary conditions
do j=2,jmm1
! east
if(n_east.eq.-1) then
denom=(uaf(im-1,j)+uab(im-1,j)-2.e0*ua(im-2,j))
if(denom.eq.0.0e0)denom=0.01e0
cl=(uab(im-1,j)-uaf(im-1,j))/denom
if(cl.gt.1.e0) cl=1.e0
if(cl.lt.0.e0) cl=0.e0
uaf(im,j)=(uab(im,j)*(1.e0-cl)+2.e0*cl*ua(im-1,j))
$ /(1.e0+cl)
vaf(im,j)=0.e0
end if
! west
if(n_west.eq.-1) then
denom=(uaf(3,j)+uab(3,j)-2.e0*ua(4,j))
if(denom.eq.0.0e0)denom=0.01e0
cl=(uab(3,j)-uaf(3,j))/denom
if(cl.gt.1.e0) cl=1.e0
if(cl.lt.0.e0) cl=0.e0
uaf(2,j)=(uab(2,j)*(1.e0-cl)+2.e0*cl*ua(3,j))
$ /(1.e0+cl)
uaf(1,j)=uaf(2,j)
vaf(1,j)=0.e0
end if
end do
do i=2,imm1 !RMY: as_is(ecoast),x_this_do_loop(stide)
! south
if(n_south.eq.-1) then
denom=(vaf(i,3)+vab(i,3)-2.e0*va(i,4))
if(denom.eq.0.0e0)denom=0.01e0
cl=(vab(i,3)-vaf(i,3))/denom
if(cl.gt.1.e0) cl=1.e0
if(cl.lt.0.e0) cl=0.e0
vaf(i,2)=(vab(i,2)*(1.e0-cl)+2.e0*cl*va(i,3))
$ /(1.e0+cl)
vaf(i,1)=vaf(i,2)
uaf(i,1)=0.e0
end if
! north
if(n_north.eq.-1) then
denom=(vaf(i,jm-1)+vab(i,jm-1)-2.e0*va(i,jm-2))
if(denom.eq.0.0e0)denom=0.01e0
cl=(vab(i,jm-1)-vaf(i,jm-1))/denom
if(cl.gt.1.e0) cl=1.e0
if(cl.lt.0.e0) cl=0.e0
vaf(i,jm)=(vab(i,jm)*(1.e0-cl)+2.e0*cl*va(i,jm-1))
$ /(1.e0+cl)
uaf(i,jm)=0.e0
end if
end do !RMY: as_is(ecoast),x_this_do_loop(stide)
do j=1,jm
do i=1,im
uaf(i,j)=uaf(i,j)*dum(i,j)
vaf(i,j)=vaf(i,j)*dvm(i,j)
end do
end do
return
else if(idx.eq.3) then
! internal (3-D) velocity boundary conditions
do k=1,kbm1
do j=2,jmm1
! east
if(n_east.eq.-1) then
denom=(uf(im-1,j,k)+ub(im-1,j,k)-2.e0*u(im-2,j,k))
if(denom.eq.0.e0)denom=0.01e0
cl=(ub(im-1,j,k)-uf(im-1,j,k))/denom
if(cl.gt.1.e0) cl=1.e0
if(cl.lt.0.e0) cl=0.e0
uf(im,j,k)=(ub(im,j,k)*(1.e0-cl)+2.e0*cl*u(im-1,j,k))
$ /(1.e0+cl)
vf(im,j,k)=0.e0
end if
! west
if(n_west.eq.-1) then
denom=(uf(3,j,k)+ub(3,j,k)-2.e0*u(4,j,k))
if(denom.eq.0.e0)denom=0.01e0
cl=(ub(3,j,k)-uf(3,j,k))/denom
if(cl.gt.1.e0) cl=1.e0
if(cl.lt.0.e0) cl=0.e0
uf(2,j,k)=(ub(2,j,k)*(1.e0-cl)+2.e0*cl*u(3,j,k))
$ /(1.e0+cl)
uf(1,j,k)=uf(2,j,k)
vf(1,j,k)=0.e0
end if
end do
do i=2,imm1 !RMY: as_is(ecoast),x_this_do_loop(stide)
! south
if(n_south.eq.-1) then
denom=(vf(i,3,k)+vb(i,3,k)-2.e0*v(i,4,k))
if(abs(denom).eq.0.0e0)denom=0.01e0
cl=(vb(i,3,k)-vf(i,3,k))/denom
if(cl.gt.1.e0) cl=1.e0
if(cl.lt.0.e0) cl=0.e0
vf(i,2,k)=(vb(i,2,k)*(1.e0-cl)+2.e0*cl*v(i,3,k))
$ /(1.e0+cl)
vf(i,1,k)=vf(i,2,k)
uf(i,1,k)=0.e0
end if
! north
if(n_north.eq.-1) then
denom=(vf(i,jm-1,k)+vb(i,jm-1,k)-2.e0*v(i,jm-2,k))
if(abs(denom).eq.0.0e0)denom=0.01e0
cl=(vb(i,jm-1,k)-vf(i,jm-1,k))/denom
if(cl.gt.1.e0) cl=1.e0
if(cl.lt.0.e0) cl=0.e0
vf(i,jm,k)=(vb(i,jm,k)*(1.e0-cl)+2.e0*cl*v(i,jm-1,k))
$ /(1.e0+cl)
uf(i,jm,k)=0.e0
end if
end do !RMY: as_is(ecoast),x_this_do_loop(stide)
end do
do k=1,kbm1
do j=1,jm
do i=1,im
uf(i,j,k)=uf(i,j,k)*dum(i,j)
vf(i,j,k)=vf(i,j,k)*dvm(i,j)
end do
end do
end do
return
else if(idx.eq.4) then
! temperature and salinity boundary conditions (using uf and vf,
! respectively)
do k=1,kbm1
do j=1,jm
! east
if(n_east.eq.-1) then
ube(j,k)=ub(im,j,k)
denom=(uf(im-1,j,k)+tb(im-1,j,k)-2.e0*t(im-2,j,k))
if(denom.eq.0.e0) denom=0.01e0
cl=(tb(im-1,j,k)-uf(im-1,j,k))/denom
if(cl.gt.1.e0) cl=1.e0
if(cl.lt.0.e0) cl=0.e0
uf(im,j,k)=(tb(im,j,k)*(1.e0-cl)+2.e0*cl*t(im-1,j,k))
$ /(1.e0+cl)
if(cl.eq.0.e0.and.ube(j,k).le.0.e0) uf(im,j,k)=tbe(j,k)
denom=(vf(im-1,j,k)+sb(im-1,j,k)-2.e0*s(im-2,j,k))
if(denom.eq.0.e0) denom=0.01e0
cl=(sb(im-1,j,k)-vf(im-1,j,k))/denom
if(cl.gt.1.e0) cl=1.e0
if(cl.lt.0.e0) cl=0.e0
vf(im,j,k)=(sb(im,j,k)*(1.e0-cl)+2.e0*cl*s(im-1,j,k))
$ /(1.e0+cl)
if(cl.eq.0.e0.and.ube(j,k).le.0.e0) vf(im,j,k)=sbe(j,k)
end if
! west
if(n_west.eq.-1) then
ubw(j,k)=ub(2,j,k)
denom=(uf(2,j,k)+tb(2,j,k)-2.e0*t(3,j,k))
if(denom.eq.0.e0) denom=0.01e0
cl=(tb(2,j,k)-uf(2,j,k))/denom
if(cl.gt.1.e0) cl=1.e0
if(cl.lt.0.e0) cl=0.e0
uf(1,j,k)=(tb(1,j,k)*(1.e0-cl)+2.e0*cl*t(2,j,k))/(1.e0+cl)
if(cl.eq.0.e0.and.ubw(j,k).ge.0.e0) uf(1,j,k)=tbw(j,k)
denom=(vf(2,j,k)+sb(2,j,k)-2.e0*s(3,j,k))
if(denom.eq.0.e0) denom=0.01e0
cl=(sb(2,j,k)-vf(2,j,k))/denom
if(cl.gt.1.e0) cl=1.e0
if(cl.lt.0.e0) cl=0.e0
vf(1,j,k)=(sb(1,j,k)*(1.e0-cl)+2.e0*cl*s(2,j,k))/(1.e0+cl)
if(cl.eq.0.e0.and.ubw(j,k).ge.0.e0) vf(1,j,k)=sbw(j,k)
end if
end do
end do
do k=1,kbm1
do j=1,jm
do i=1,im
uf(i,j,k)=uf(i,j,k)*fsm(i,j)
vf(i,j,k)=vf(i,j,k)*fsm(i,j)
end do
end do
end do
return
else if(idx.eq.5) then
! vertical velocity boundary conditions
do k=1,kbm1
do j=1,jm
do i=1,im
w(i,j,k)=w(i,j,k)*fsm(i,j)
end do
end do
end do
return
else if(idx.eq.6) then
! q2 and q2l boundary conditions
do k=1,kb
do j=1,jm
if(n_east.eq.-1) then
uf(im,j,k)=1.e-10
vf(im,j,k)=1.e-10
end if
if(n_west.eq.-1) then
uf(1,j,k)=1.e-10
vf(1,j,k)=1.e-10
end if
end do
do j=1,jm
do i=1,im
uf(i,j,k)=uf(i,j,k)*fsm(i,j)
vf(i,j,k)=vf(i,j,k)*fsm(i,j)
end do
end do
end do
return
endif
end
!RMY: as_is(ecoast),end_of_file(stide)
! _____________________________________________________________________
subroutine lateral_bc
! create variable lateral boundary conditions
use setvars
implicit none
integer nz
!RMY parameter(nz=33) !RMY: 40(ecoast),33(pomtc); use pom.nml's nl
integer i,j,k,ntime,ibc
real tbc,fold,fnew
real z0(nl),hs(im,jm) !RMY: use nl
real t_w(jm,nl),s_w(jm,nl),t_e(jm,nl),s_e(jm,nl),q, !RMY: use nl
$ t_n(im,nl),s_n(im,nl),t_s(im,nl),s_s(im,nl) !RMY: use nl
real ts1(im,jm,nl),ss1(im,jm,nl) !RMY: use nl
real ts2(im,jm,kb),ss2(im,jm,kb)
nz=nl !RMY: set nz=nl locally
tbc=30 ! time between bc files (days)
ibc=int(tbc*86400.e0/dti)
ntime=time/tbc
! read bc data
! read initial bc file
if (iint.eq.1) then
call read_boundary_conditions_pnetcdf(iint/ibc,nz,z0,
$ t_w,s_w,uabwf,t_e,s_e,uabef,t_n,s_n,vabnf,t_s,s_s,vabsf)
! south
do i=1,im
do j=1,jm
hs(i,j)=h(i,2)
do k=1,nz
ts1(i,j,k)=t_s(i,k)
ss1(i,j,k)=s_s(i,k)
end do
end do
end do
call ztosig(z0,ts1,zz,hs,ts2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
call ztosig(z0,ss1,zz,hs,ss2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
do i=1,im
do k=1,kb
tbsf(i,k)=ts2(i,jm/2,k)
sbsf(i,k)=ss2(i,jm/2,k)
end do
end do
! north
do i=1,im
do j=1,jm
hs(i,j)=h(i,jmm1)
do k=1,nz
ts1(i,j,k)=t_n(i,k)
ss1(i,j,k)=s_n(i,k)
end do
end do
end do
call ztosig(z0,ts1,zz,hs,ts2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
call ztosig(z0,ss1,zz,hs,ss2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
do i=1,im
do k=1,kb
tbnf(i,k)=ts2(i,jm/2,k)
sbnf(i,k)=ss2(i,jm/2,k)
end do
end do
! east
do i=1,im
do j=1,jm
hs(i,j)=h(imm1,j)
do k=1,nz
ts1(i,j,k)=t_e(j,k)
ss1(i,j,k)=s_e(j,k)
end do
end do
end do
call ztosig(z0,ts1,zz,hs,ts2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
call ztosig(z0,ss1,zz,hs,ss2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
do j=1,jm
do k=1,kb
tbef(j,k)=ts2(im/2,j,k)
sbef(j,k)=ss2(im/2,j,k)
end do
end do
! west
do i=1,im
do j=1,jm
hs(i,j)=h(2,j)
do k=1,nz
ts1(i,j,k)=t_w(j,k)
ss1(i,j,k)=s_w(j,k)
end do
end do
end do
call ztosig(z0,ts1,zz,hs,ts2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
call ztosig(z0,ss1,zz,hs,ss2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
do j=1,jm
do k=1,kb
tbwf(j,k)=ts2(im/2,j,k)
sbwf(j,k)=ss2(im/2,j,k)
end do
end do
end if
! read bc file corresponding to next time
if (iint.eq.1 .or. mod(iint,ibc).eq.0.) then
do j=1,jm
do k=1,kb
tbwb(j,k)=tbwf(j,k)
sbwb(j,k)=sbwf(j,k)
tbeb(j,k)=tbef(j,k)
sbeb(j,k)=sbef(j,k)
end do
uabwb(j)=uabwf(j)
uabeb(j)=uabef(j)
end do
do i=1,im
do k=1,kb
tbnb(i,k)=tbnf(i,k)
sbnb(i,k)=sbnf(i,k)
tbsb(i,k)=tbsf(i,k)
sbsb(i,k)=sbsf(i,k)
end do
vabnb(i)=vabnf(i)
vabsb(i)=vabsf(i)
end do
if (iint.ne.iend) then
call read_boundary_conditions_pnetcdf((iint+ibc)/ibc,nz,z0,
$ t_w,s_w,uabwf,t_e,s_e,uabef,t_n,s_n,vabnf,t_s,s_s,vabsf)
! south
do i=1,im
do j=1,jm
hs(i,j)=h(i,2)
do k=1,nz
ts1(i,j,k)=t_s(i,k)
ss1(i,j,k)=s_s(i,k)
end do
end do
end do
call ztosig(z0,ts1,zz,hs,ts2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
call ztosig(z0,ss1,zz,hs,ss2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
do i=1,im
do k=1,kb
tbsf(i,k)=ts2(i,jm/2,k)
sbsf(i,k)=ss2(i,jm/2,k)
end do
end do
! north
do i=1,im
do j=1,jm
hs(i,j)=h(i,jmm1)
do k=1,nz
ts1(i,j,k)=t_n(i,k)
ss1(i,j,k)=s_n(i,k)
end do
end do
end do
call ztosig(z0,ts1,zz,hs,ts2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
call ztosig(z0,ss1,zz,hs,ss2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
do i=1,im
do k=1,kb
tbnf(i,k)=ts2(i,jm/2,k)
sbnf(i,k)=ss2(i,jm/2,k)
end do
end do
! east
do i=1,im
do j=1,jm
hs(i,j)=h(imm1,j)
do k=1,nz
ts1(i,j,k)=t_e(j,k)
ss1(i,j,k)=s_e(j,k)
end do
end do
end do
call ztosig(z0,ts1,zz,hs,ts2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
call ztosig(z0,ss1,zz,hs,ss2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
do j=1,jm
do k=1,kb
tbef(j,k)=ts2(im/2,j,k)
sbef(j,k)=ss2(im/2,j,k)
end do
end do
! west
do i=1,im
do j=1,jm
hs(i,j)=h(2,j)
do k=1,nz
ts1(i,j,k)=t_w(j,k)
ss1(i,j,k)=s_w(j,k)
end do
end do
end do
call ztosig(z0,ts1,zz,hs,ts2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
call ztosig(z0,ss1,zz,hs,ss2,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
do j=1,jm
do k=1,kb
tbwf(j,k)=ts2(im/2,j,k)
sbwf(j,k)=ss2(im/2,j,k)
end do
end do
! end
end if
end if
! linear interpolation in time
fnew=time/tbc-ntime
fold=1.-fnew
do j=1,jm
do k=1,kb
tbw(j,k)=fold*tbwb(j,k)+fnew*tbwf(j,k)
sbw(j,k)=fold*sbwb(j,k)+fnew*sbwf(j,k)
tbe(j,k)=fold*tbeb(j,k)+fnew*tbef(j,k)
sbe(j,k)=fold*sbeb(j,k)+fnew*sbef(j,k)
end do
uabe(j)=fold*uabeb(j)+fnew*uabef(j)
end do
do i=1,im
do k=1,kb
tbn(i,k)=fold*tbnb(i,k)+fnew*tbnf(i,k)
sbn(i,k)=fold*sbnb(i,k)+fnew*sbnf(i,k)
tbs(i,k)=fold*tbsb(i,k)+fnew*tbsf(i,k)
sbs(i,k)=fold*sbsb(i,k)+fnew*sbsf(i,k)
end do
vabn(i)=fold*vabnb(i)+fnew*vabnf(i)
vabs(i)=fold*vabsb(i)+fnew*vabsf(i)
end do
return
end
!_______________________________________________________________________
subroutine wind
! read and interpolate (in time) wind stress
use setvars
implicit none
integer i,j,ntime,iwind
real twind,fold,fnew
twind=30 ! time between wind files (days) !RMY: 30(sbpom),?(pomtc)
iwind=int(twind*86400.e0/dti)
! read wind stress data
! read initial wind file
if (iint.eq.1) call read_wind_pnetcdf(iint/iwind,wusurff,wvsurff)
! read wind file corresponding to next time
if (iint.eq.1 .or. mod(iint,iwind).eq.0.) then
do i=1,im
do j=1,jm
wusurfb(i,j)=wusurff(i,j)
wvsurfb(i,j)=wvsurff(i,j)
end do
end do
if (iint.ne.iend) then
call read_wind_pnetcdf((iint+iwind)/iwind,wusurff,wvsurff)
end if
end if
! linear interpolation in time
ntime=time/twind
fnew=time/twind-ntime
fold=1.-fnew
do i=1,im
do j=1,jm
wusurf(i,j)=fold*wusurfb(i,j)+fnew*wusurff(i,j)
wvsurf(i,j)=fold*wvsurfb(i,j)+fnew*wvsurff(i,j)
end do
end do
return
end
!_______________________________________________________________________
subroutine heat
! read and interpolate heat flux in time
use setvars
implicit none
integer i,j,ntime,iheat
real theat,fold,fnew
theat=1. ! time between wind forcing (days)
iheat=int(theat*86400.e0/dti)
! read wind stress data
! read initial heat file
if (iint.eq.1) call read_heat_pnetcdf(iint/iheat,wtsurff,swradf)
! read heat forcing corresponding to next twind
if (iint.eq.1 .or. mod(iint,iheat).eq.0.) then
do i=1,im
do j=1,jm
wtsurfb(i,j)=wtsurff(i,j)
swradb(i,j)=swradf(i,j)
end do
end do
call read_heat_pnetcdf((iint+iheat)/iheat,wtsurff,swradf)
end if
! linear interpolation in time
ntime=time/theat
fnew=time/theat-ntime
fold=1.-fnew
do i=1,im
do j=1,jm
wtsurf(i,j)=fold*wtsurfb(i,j)+fnew*wtsurff(i,j)
swrad(i,j)=fold*swradb(i,j)+fnew*swradf(i,j)
end do
end do
return
end
!_______________________________________________________________________
subroutine water
! read and interpolate water flux in time
use setvars
implicit none
integer i,j,ntime,iwater
real twater,fold,fnew
twater=1. ! time between wind forcing (days)
iwater=int(twater*86400.e0/dti)
! read wind stress data
! read initial water file
if (iint.eq.1) call read_water_pnetcdf(iint/iwater,wssurff)
! read heat forcing corresponding to next twind
if (iint.eq.1 .or. mod(iint,iwater).eq.0.) then
do i=1,im
do j=1,jm
wssurfb(i,j)=wssurff(i,j)
end do
end do
call read_water_pnetcdf((iint+iwater)/iwater,wssurff)
end if
! linear interpolation in time
ntime=time/twater
fnew=time/twater-ntime
fold=1.-fnew
do i=1,im
do j=1,jm
wssurf(i,j)=fold*wssurfb(i,j)+fnew*wssurff(i,j)
end do
end do
return
end
! _____________________________________________________________________
subroutine restore_interior
! read, interpolate (in time) and apply restore interior data
use setvars
implicit none
integer nz
!RMY parameter(nz=33) !RMY: 40(ecoast),33(pomtc); use pom.nml's nl
real z0(nl),f0(im,jm,nl) !RMY: use nl
integer i,j,k,ntime,irst
real trst,fold,fnew
nz=nl !RMY: set nz=nl locally
trst=30 ! time between restore files (days)
irst=int(trst*86400.e0/dti)
ntime=time/trst
! read restore data
! read initial restore file
if (iint.eq.2) then
call read_restore_t_interior_pnetcdf(iint/irst,nz,z0,f0)
call ztosig(z0,f0,zz,h,trstrf,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
call read_restore_s_interior_pnetcdf(iint/irst,nz,z0,f0)
call ztosig(z0,f0,zz,h,srstrf,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
call read_restore_tau_interior_pnetcdf(iint/irst,nz,z0,f0)
call ztosig(z0,f0,zz,h,taurstrf,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
end if
! read restore file corresponding to next time
if (iint.eq.2 .or. mod(iint,irst).eq.0.) then
do k=1,kbm1
do i=1,im
do j=1,jm
trstrb(i,j,k)=trstrf(i,j,k)
srstrb(i,j,k)=srstrf(i,j,k)
taurstrb(i,j,k)=taurstrf(i,j,k)
end do
end do
end do
if (iint.ne.iend) then
call read_restore_t_interior_pnetcdf((iint+irst)/irst,nz,z0,
$ f0)
call ztosig(z0,f0,zz,h,trstrf,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
call read_restore_s_interior_pnetcdf((iint+irst)/irst,nz,z0,
$ f0)
call ztosig(z0,f0,zz,h,srstrf,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
call read_restore_tau_interior_pnetcdf((iint+irst)/irst,nz,
$ z0,f0)
call ztosig(z0,f0,zz,h,taurstrf,im,jm,nz,kb,
$ n_west,n_east,n_south,n_north)
end if
end if
! linear interpolation in time
fnew=time/trst-ntime
fold=1.-fnew
do k=1,kbm1
do i=1,im
do j=1,jm
trstr(i,j,k)=fold*trstrb(i,j,k)+fnew*trstrf(i,j,k)
srstr(i,j,k)=fold*srstrb(i,j,k)+fnew*srstrf(i,j,k)
taurstr(i,j,k)=fold*taurstrb(i,j,k)+fnew*taurstrf(i,j,k)
end do
end do
end do
! restore
do k=1,kbm1
do i=1,im
do j=1,jm
t(i,j,k)=t(i,j,k)+2.*dti/86400.*taurstr(i,j,k)*
$ (trstr(i,j,k)-t(i,j,k))
tb(i,j,k)=tb(i,j,k)+2.*dti/86400.*taurstr(i,j,k)*
$ (trstr(i,j,k)-tb(i,j,k))
s(i,j,k)=s(i,j,k)+2.*dti/86400.*taurstr(i,j,k)*
$ (srstr(i,j,k)-s(i,j,k))
sb(i,j,k)=sb(i,j,k)+2.*dti/86400.*taurstr(i,j,k)*
$ (srstr(i,j,k)-sb(i,j,k))
end do
end do
end do
! mask
do k=1,kbm1
do j=1,jm
do i=1,im
t(i,j,k)=t(i,j,k)*fsm(i,j)
tb(i,j,k)=tb(i,j,k)*fsm(i,j)
s(i,j,k)=s(i,j,k)*fsm(i,j)
sb(i,j,k)=sb(i,j,k)*fsm(i,j)
end do
end do
end do
return
end