C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_obcs_ageos.F,v 1.10 2014/10/20 03:16:12 gforget Exp $
C $Name: $
#include "ECCO_OPTIONS.h"
#ifdef ALLOW_CTRL
# include "CTRL_OPTIONS.h"
#endif
subroutine COST_OBCS_AGEOS( myiter, mytime, mythid )
c ==================================================================
c SUBROUTINE cost_obcs_ageos
c ==================================================================
c
c o cost function contribution obc -- Ageostrophic boundary flow.
c
c started: G. Gebbie gebbie@mit.edu 4-Feb-2003
c
c warning: masks redundantly assume no gradient of topography at
c boundary.
c
c ==================================================================
c SUBROUTINE cost_obcs_ageos
c ==================================================================
implicit none
c == global variables ==
#include "EEPARAMS.h"
#include "SIZE.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "PARAMS.h"
#ifdef ALLOW_OBCS
# include "OBCS_GRID.h"
#endif
#ifdef ALLOW_CAL
# include "cal.h"
#endif
#ifdef ALLOW_ECCO
# include "ecco_cost.h"
#endif
#ifdef ALLOW_CTRL
# include "CTRL_SIZE.h"
# include "ctrl.h"
# include "ctrl_dummy.h"
# include "optim.h"
# include "CTRL_OBCS.h"
#endif
c == routine arguments ==
integer myiter
_RL mytime
integer mythid
#if (defined (ALLOW_CTRL)
defined (ALLOW_OBCS)
defined (ALLOW_ECCO))
c == local variables ==
integer bi,bj
integer i,j,k
integer itlo,ithi
integer jtlo,jthi
integer jmin,jmax
integer imin,imax
integer irec
integer levmon
integer levoff
integer iltheta
integer ilsalt
integer iluvel
integer ilvvel
integer ip1, jp1
_RL fctile
_RL fcthread
_RL rholoc (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
_RL xzgrdrho(1-olx:snx+olx,Nr,nsx,nsy)
_RL yzgrdrho(1-oly:sny+oly,Nr,nsx,nsy)
_RL xzdvel1 (1-olx:snx+olx,nr,nsx,nsy)
_RL xzdvel2 (1-olx:snx+olx,nr,nsx,nsy)
_RL yzdvel1 (1-oly:sny+oly,nr,nsx,nsy)
_RL yzdvel2 (1-oly:sny+oly,nr,nsx,nsy)
_RL maskxzageos (1-olx:snx+olx,nr,nsx,nsy)
_RL maskyzageos (1-oly:sny+oly,nr,nsx,nsy)
_RL dummy
character*(80) fnametheta
character*(80) fnamesalt
character*(80) fnameuvel
character*(80) fnamevvel
logical doglobalread
logical ladinit
character*(MAX_LEN_MBUF) msgbuf
c == external functions ==
integer ilnblnk
external
c == end of interface ==
jtlo = mybylo(mythid)
jthi = mybyhi(mythid)
itlo = mybxlo(mythid)
ithi = mybxhi(mythid)
jmin = 1
jmax = sny
imin = 1
imax = snx
c-- Read tiled data.
doglobalread = .false.
ladinit = .false.
#ifdef OBCS_AGEOS_COST_CONTRIBUTION
#ifdef ECCO_VERBOSE
_BEGIN_MASTER( mythid )
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,i8.8)')
& ' cost_obcs_ageos: number of records to process = ',nmonsrec
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
_END_MASTER( mythid )
#endif
if (optimcycle .ge. 0) then
iltheta = ilnblnk( tbarfile )
write(fnametheta(1:80),'(2a,i10.10)')
& tbarfile(1:iltheta),'.',optimcycle
ilsalt = ilnblnk( sbarfile )
write(fnamesalt(1:80),'(2a,i10.10)')
& sbarfile(1:ilsalt),'.',optimcycle
iluvel = ilnblnk( ubarfile )
write(fnameuvel(1:80),'(2a,i10.10)')
& ubarfile(1:iluvel),'.',optimcycle
ilvvel = ilnblnk( vbarfile )
write(fnamevvel(1:80),'(2a,i10.10)')
& vbarfile(1:ilvvel),'.',optimcycle
endif
fcthread = 0. _d 0
fctile = 0. _d 0
cgg Code safe: always initialize to zero.
do bj = jtlo,jthi
do bi = itlo,ithi
do k = 1,nr
do i = 1-olx,snx+olx
maskxzageos(i,k,bi,bj) = 0. _d 0
xzdvel1(i,k,bi,bj) = 0. _d 0
xzdvel2(i,k,bi,bj) = 0. _d 0
xzgrdrho(i,k,bi,bj) = 0. _d 0
enddo
do j = 1-oly,sny+oly
maskyzageos(j,k,bi,bj) = 0. _d 0
yzdvel1(j,k,bi,bj) = 0. _d 0
yzdvel2(j,k,bi,bj) = 0. _d 0
yzgrdrho(j,k,bi,bj) = 0. _d 0
enddo
enddo
enddo
enddo
do bj = jtlo,jthi
do bi = itlo,ithi
do j = 1-oly,sny+oly
do i = 1-olx,snx+olx
rholoc(i,j,bi,bj) = 0. _d 0
enddo
enddo
enddo
enddo
c-- Loop over records.
do irec = 1,nmonsrec
c-- Read time averages and the monthly mean data.
call ACTIVE_READ_XYZ( fnametheta, tbar, irec,
& doglobalread, ladinit,
& optimcycle, mythid,
& xx_tbar_mean_dummy )
c-- Read time averages and the monthly mean data.
call ACTIVE_READ_XYZ( fnamesalt, sbar, irec,
& doglobalread, ladinit,
& optimcycle, mythid,
& xx_sbar_mean_dummy )
c-- Read time averages and the monthly mean data.
call ACTIVE_READ_XYZ( fnameuvel, ubar, irec,
& doglobalread, ladinit,
& optimcycle, mythid,
& xx_ubar_mean_dummy )
c-- Read time averages and the monthly mean data.
call ACTIVE_READ_XYZ( fnamevvel, vbar, irec,
& doglobalread, ladinit,
& optimcycle, mythid,
& xx_vbar_mean_dummy )
cgg Minor problem : grad T,S needs overlap values.
_BARRIER
_EXCH_XYZ_RL(tbar,myThid)
_EXCH_XYZ_RL(sbar,myThid)
do bj = jtlo,jthi
do bi = itlo,ithi
#ifdef ALLOW_OBCSN_CONTROL
jp1 = 0
cgg Make a mask for the velocity shear comparison.
do k = 1,nr-1
do i = imin, imax
j = OB_Jn(i,bi,bj)
cgg All these points need to be wet.
if ( j.eq.OB_indexNone ) then
maskxzageos(i,k,bi,bj) = 0.
else
maskxzageos(i,k,bi,bj) =
& hfacC(i,j+jp1,k,bi,bj)*hfacC(i+1,j+jp1,k,bi,bj) *
& hfacC(i-1,j+jp1,k,bi,bj)*hfacC(i,j+jp1,k+1,bi,bj)*
& hfacC(i-1,j+jp1,k+1,bi,bj)*hfacC(i+1,j+jp1,k+1,bi,bj)*
& hfacS(i,j+jp1,k,bi,bj)*hfacS(i,j+jp1,k+1,bi,bj)
endif
enddo
enddo
do k = 1,nr
C-- jmc: both calls below are wrong if more than 1 tile => stop here
IF ( bi.NE.1 .OR. bj.NE.1 )
& STOP 'COST_OBCS_AGEOS wrong with more than 1 tile/proc'
call FIND_RHO_2D(
I iMin, iMax, jMin, jMax, k,
I tbar(1-OLx,1-OLy,k,bi,bj),
I sbar(1-OLx,1-OLy,k,bi,bj),
O rholoc,
I k, bi, bj, myThid )
_EXCH_XY_RL(rholoc , myThid)
cgg Compute centered difference horizontal gradient on bdy.
do i = imin, imax
j = OB_Jn(i,bi,bj)
if ( j.eq.OB_indexNone ) j = 1
xzgrdrho(i,k,bi,bj) =
& (rholoc(i-1,j+jp1,bi,bj)-rholoc(i+1,j+jp1,bi,bj))
& /(2.*dxc(i,j+jp1,bi,bj))
enddo
enddo
cgg Compute vertical shear from geostrophy/thermal wind.
cgg Above level 4 needs not be geostrophic.
cgg Please get rid of the "4" ASAP. Ridiculous.
do k = 4,Nr-1
do i = imin,imax
j = OB_Jn(i,bi,bj)
if ( j.eq.OB_indexNone ) j = 1
xzdvel1(i,k,bi,bj) = vbar(i,j+jp1,k ,bi,bj)
& - vbar(i,j+jp1,k+1,bi,bj)
xzdvel2(i,k,bi,bj)=((xzgrdrho(i,k,bi,bj)*delz(k)/2.)+
& (xzgrdrho(i,k+1,bi,bj)*delz(k+1)/2.))
& * gravity / f0 / rhonil
fctile = fctile + 100*wcurrent(k,bi,bj) *
& maskxzageos(i,k,bi,bj)*
& (xzdvel2(i,k,bi,bj) - xzdvel1(i,k,bi,bj))*
& (xzdvel2(i,k,bi,bj) - xzdvel1(i,k,bi,bj))
if (maskxzageos(i,k,bi,bj) .ne. 0) then
cgg print*,'XX shear,geos shear',xzdvel1(i,k,bi,bj),xzdvel2(i,k,bi,bj)
cgg print*,'i,j,k,fctile N',i,j,k,fctile
endif
enddo
enddo
c-- End of loop over layers.
#endif /* ALLOW_OBCSN_CONTROL */
#ifdef ALLOW_OBCSS_CONTROL
jp1 = 1
cgg Make a mask for the velocity shear comparison.
do k = 1,nr-1
do i = imin, imax
j = OB_Js(i,bi,bj)
if ( j.eq.OB_indexNone ) then
maskxzageos(i,k,bi,bj) = 0.
else
cgg All these points need to be wet.
maskxzageos(i,k,bi,bj) =
& hfacC(i,j+jp1,k,bi,bj)*hfacC(i+1,j+jp1,k,bi,bj) *
& hfacC(i-1,j+jp1,k,bi,bj)*hfacC(i,j+jp1,k+1,bi,bj)*
& hfacC(i-1,j+jp1,k+1,bi,bj)*hfacC(i+1,j+jp1,k+1,bi,bj)*
& hfacS(i,j+jp1,k,bi,bj)*hfacS(i,j+jp1,k+1,bi,bj)
endif
enddo
enddo
do k = 1,nr
C-- jmc: both calls below are wrong if more than 1 tile => stop here
IF ( bi.NE.1 .OR. bj.NE.1 )
& STOP 'COST_OBCS_AGEOS wrong with more than 1 tile/proc'
call FIND_RHO_2D(
I iMin, iMax, jMin, jMax, k,
I tbar(1-OLx,1-OLy,k,bi,bj),
I sbar(1-OLx,1-OLy,k,bi,bj),
O rholoc,
I k, bi, bj, myThid )
_EXCH_XY_RL(rholoc , myThid)
cgg Compute centered difference horizontal gradient on bdy.
do i = imin, imax
j = OB_Js(i,bi,bj)
if ( j.eq.OB_indexNone ) j = 1
xzgrdrho(i,k,bi,bj) =
& (rholoc(i-1,j+jp1,bi,bj)-rholoc(i+1,j+jp1,bi,bj))
& /(2.*dxc(i,j+jp1,bi,bj))
enddo
enddo
cgg Compute vertical shear from geostrophy/thermal wind.
do k = 4,Nr-1
do i = imin,imax
j = OB_Js(i,bi,bj)
if ( j.eq.OB_indexNone ) j = 1
cgg Retrieve the model vertical shear.
xzdvel1(i,k,bi,bj) = vbar(i,j+jp1,k ,bi,bj)
& - vbar(i,j+jp1,k+1,bi,bj)
cgg Compute vertical shear from geostrophy/thermal wind.
xzdvel2(i,k,bi,bj) =((xzgrdrho(i,k,bi,bj)*delz(k)/2.)+
& (xzgrdrho(i,k+1,bi,bj)*delz(k+1)/2.))
& * gravity / f0 /rhonil
cgg Make a comparison.
fctile = fctile + 100*wcurrent(k,bi,bj) *
& maskxzageos(i,k,bi,bj)*
& (xzdvel2(i,k,bi,bj) - xzdvel1(i,k,bi,bj))*
& (xzdvel2(i,k,bi,bj) - xzdvel1(i,k,bi,bj))
cgg print*,'fctile S',fctile
enddo
enddo
c-- End of loop over layers.
#endif
#ifdef ALLOW_OBCSW_CONTROL
ip1 = 1
cgg Make a mask for the velocity shear comparison.
do k = 1,nr-1
do j = jmin, jmax
i = OB_Iw(j,bi,bj)
cgg All these points need to be wet.
if ( i.eq.OB_indexNone ) then
maskyzageos(j,k,bi,bj) = 0.
else
maskyzageos(j,k,bi,bj) =
& hfacC(i+ip1,j,k,bi,bj)*hfacC(i+ip1,j+1,k,bi,bj) *
& hfacC(i+ip1,j-1,k,bi,bj)*hfacC(i+ip1,j,k+1,bi,bj)*
& hfacC(i+ip1,j-1,k+1,bi,bj)*hfacC(i+ip1,j+1,k+1,bi,bj)*
& hfacW(i+ip1,j,k,bi,bj)*hfacW(i+ip1,j,k+1,bi,bj)
endif
enddo
enddo
do k = 1,nr
IF ( bi.NE.1 .OR. bj.NE.1 )
& STOP 'COST_OBCS_AGEOS wrong with more than 1 tile/proc'
call FIND_RHO_2D(
I iMin, iMax, jMin, jMax, k,
I tbar(1-OLx,1-OLy,k,bi,bj),
I sbar(1-OLx,1-OLy,k,bi,bj),
O rholoc,
I k, bi, bj, myThid )
_EXCH_XY_RL(rholoc , myThid)
cgg Compute centered difference horizontal gradient on bdy.
do j = jmin, jmax
i = OB_Iw(j,bi,bj)
if ( i.eq.OB_indexNone ) i = 1
cgg Negative sign due to geostrophy.
yzgrdrho(j,k,bi,bj) =
& (rholoc(i+ip1,j+1,bi,bj)-rholoc(i+ip1,j-1,bi,bj))
& /(2.*dyc(i+ip1,j,bi,bj))
enddo
enddo
cgg Compute vertical shear from geostrophy/thermal wind.
do k = 4,Nr-1
do j = jmin,jmax
i = OB_Iw(j,bi,bj)
if ( i.eq.OB_indexNone ) i = 1
cgg Retrieve the model vertical shear.
yzdvel1(j,k,bi,bj) = ubar(i+ip1,j,k ,bi,bj)
& - ubar(i+ip1,j,k+1,bi,bj)
cgg Compute vertical shear from geostrophy/thermal wind.
yzdvel2(j,k,bi,bj) =((yzgrdrho(j,k ,bi,bj)*delz(k)/2.)+
& (yzgrdrho(j,k+1,bi,bj)*delz(k+1)/2.))
& * gravity / f0 / rhonil
cgg Make a comparison.
fctile = fctile + 100*wcurrent(k,bi,bj) *
& maskyzageos(j,k,bi,bj) *
& (yzdvel2(j,k,bi,bj) - yzdvel1(j,k,bi,bj))*
& (yzdvel2(j,k,bi,bj) - yzdvel1(j,k,bi,bj))
enddo
enddo
c-- End of loop over layers.
#endif
#ifdef ALLOW_OBCSE_CONTROL
ip1 = 0
cgg Make a mask for the velocity shear comparison.
do k = 1,nr-1
do j = jmin, jmax
i = OB_Ie(j,bi,bj)
if ( i.eq.OB_indexNone ) then
maskyzageos(j,k,bi,bj) =0.
else
cgg All these points need to be wet.
maskyzageos(j,k,bi,bj) =
& hfacC(i+ip1,j,k,bi,bj)*hfacC(i+ip1,j+1,k,bi,bj) *
& hfacC(i+ip1,j-1,k,bi,bj)*hfacC(i+ip1,j,k+1,bi,bj)*
& hfacC(i+ip1,j-1,k+1,bi,bj)*hfacC(i+ip1,j+1,k+1,bi,bj)*
& hfacW(i+ip1,j,k,bi,bj)*hfacW(i+ip1,j,k+1,bi,bj)
endif
enddo
enddo
do k = 1,nr
IF ( bi.NE.1 .OR. bj.NE.1 )
& STOP 'COST_OBCS_AGEOS wrong with more than 1 tile/proc'
call FIND_RHO_2D(
I iMin, iMax, jMin, jMax, k,
I tbar(1-OLx,1-OLy,k,bi,bj),
I sbar(1-OLx,1-OLy,k,bi,bj),
O rholoc,
I k, bi, bj, myThid )
_EXCH_XY_RL(rholoc , myThid)
cgg Compute centered difference horizontal gradient on bdy.
do j = jmin, jmax
i = OB_Ie(j,bi,bj)
if ( i.eq.OB_indexNone ) i = 1
cgg Negative sign due to geostrophy.
yzgrdrho(j,k,bi,bj) =
& (rholoc(i+ip1,,j+1,bi,bj)-rholoc(i+ip1,j-1,bi,bj))
& /(2.*dyc(i+ip1,j,bi,bj))
enddo
enddo
cgg Compute vertical shear from geostrophy/thermal wind.
do k = 4,Nr-1
do j = jmin,jmax
i = OB_Ie(j,bi,bj)
if ( i.eq.OB_indexNone ) i = 1
cgg Retrieve the model vertical shear.
yzdvel1(j,k,bi,bj) = ubar(i+ip1,j,k ,bi,bj)
& - ubar(i+ip1,j,k+1,bi,bj)
cgg Compute vertical shear from geostrophy/thermal wind.
yzdvel2(j,k,bi,bj) =((yzgrdrho(j,k ,bi,bj)*delz(k)/2.)+
& (yzgrdrho(j,k+1,bi,bj)*delz(k+1)/2.))
& * gravity / f0 /rhonil
cgg Make a comparison.
fctile = fctile + 100*wcurrent(k,bi,bj) *
& maskyzageos(j,k,bi,bj) *
& (yzdvel2(j,k,bi,bj) - yzdvel1(j,k,bi,bj))*
& (yzdvel2(j,k,bi,bj) - yzdvel1(j,k,bi,bj))
enddo
enddo
c-- End of loop over layers.
#endif
fcthread = fcthread + fctile
objf_ageos(bi,bj) = objf_ageos(bi,bj) + fctile
#ifdef ECCO_VERBOSE
c-- Print cost function for each tile in each thread.
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
& ' cost_Theta: irec,bi,bj = ',irec,bi,bj
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,d22.15)')
& ' cost function (temperature) = ',
& fctile
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
#endif
enddo
enddo
#ifdef ECCO_VERBOSE
c-- Print cost function for all tiles.
_GLOBAL_SUM_RL( fcthread , myThid )
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,i8.8)')
& ' cost_Theta: irec = ',irec
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,a,d22.15)')
& ' global cost function value',
& ' (temperature) = ',fcthread
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
#endif
enddo
c-- End of loop over records.
#endif
#endif /* ALLOW_CTRL, ALLOW_OBCS and ALLOW_ECCO */
return
end