C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_ssh.F,v 1.24 2015/10/24 18:49:43 gforget Exp $
C $Name: $
#include "ECCO_OPTIONS.h"
subroutine COST_SSH(
I myiter,
I mytime,
I mythid
& )
c ==================================================================
c SUBROUTINE cost_ssh
c ==================================================================
c
c o Evaluate cost function contribution of sea surface height.
c using of geoid error covariances requires regular model grid
c
c started: Detlef Stammer, Ralf Giering Jul-1996
c
c changed: Christian Eckert eckert@mit.edu 25-Feb-2000
c
c - Restructured the code in order to create a package
c for the MITgcmUV.
c
c changed: Ralf Giering Ralf.Giering@FastOpt.de 12-Jun-2001
c
c - totally rewrite for parallel processing
c
c ==================================================================
c SUBROUTINE cost_ssh
c ==================================================================
implicit none
c == global variables ==
#ifdef ALLOW_SSH_COST_CONTRIBUTION
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"
#include "ecco_cost.h"
#include "CTRL_SIZE.h"
#include "ctrl.h"
#include "ctrl_dummy.h"
#include "optim.h"
#include "DYNVARS.h"
#ifdef ALLOW_PROFILES
# include "PROFILES_SIZE.h"
# include "profiles.h"
#endif
#endif
c == routine arguments ==
integer myiter
_RL mytime
integer mythid
#ifdef ALLOW_SSH_COST_CONTRIBUTION
c == local variables ==
integer bi,bj
integer i,j
integer itlo,ithi
integer jtlo,jthi
integer jmin,jmax
integer imin,imax
integer irec
integer ilps
logical doglobalread
logical ladinit
_RL offset
_RL costmean
_RL offset_sum
_RL psmean ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
_RL wwwtp ( 1-olx:snx+olx, 1-oly:sny+oly )
_RL wwwers ( 1-olx:snx+olx, 1-oly:sny+oly )
_RL wwwgfo ( 1-olx:snx+olx, 1-oly:sny+oly )
_RL junk
character*(80) fname
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-- Initialise local variables.
costmean = 0. _d 0
do j = jmin, jmax
do i = imin, imax
wwwtp(i,j) = 0. _d 0
wwwers(i,j) = 0. _d 0
wwwgfo(i,j) = 0. _d 0
enddo
enddo
c-- First, read tiled data.
doglobalread = .false.
ladinit = .false.
write(fname(1:80),'(80a)') ' '
ilps=ilnblnk( psbarfile )
write(fname(1:80),'(2a,i10.10)')
& psbarfile(1:ilps),'.',optimcycle
c-- ============
c-- Mean values.
c-- ============
do bj = jtlo,jthi
do bi = itlo,ithi
do j = jmin,jmax
do i = imin,imax
psmean(i,j,bi,bj) = 0. _d 0
enddo
enddo
enddo
enddo
c-- Read mean field and generate mask
c-- Convert mean ssh from cm to m
#ifdef ALLOW_SSH_MEAN_COST_CONTRIBUTION
call MDSREADFIELD( mdtdatfile, cost_iprec, cost_yftype, 1,
& mdt, 1, mythid )
do bj = jtlo,jthi
do bi = itlo,ithi
do j = jmin,jmax
do i = imin,imax
mdtmask(i,j,bi,bj)=maskC(i,j,1,bi,bj)
if (mdt(i,j,bi,bj) .lt. -9990. _d 0) then
mdtmask(i,j,bi,bj) = 0. _d 0
endif
if ( R_low(i,j,bi,bj) .GT. -200. ) then
mdtmask(i,j,bi,bj) = 0. _d 0
endif
mdt(i,j,bi,bj) = mdt(i,j,bi,bj)*
& mdtmask(i,j,bi,bj)*0.01 _d 0
enddo
enddo
enddo
enddo
#endif /* ALLOW_SSH_MEAN_COST_CONTRIBUTION */
c-- Loop over records for the first time.
do irec = 1, ndaysrec
c-- Compute the mean over all psbar records.
call ACTIVE_READ_XY( fname, psbar, irec, doglobalread,
& ladinit, optimcycle, mythid,
& xx_psbar_mean_dummy )
do bj = jtlo,jthi
do bi = itlo,ithi
do j = jmin,jmax
do i = imin,imax
psmean(i,j,bi,bj) = psmean(i,j,bi,bj) +
& psbar(i,j,bi,bj)/
& float(ndaysrec)
enddo
enddo
enddo
enddo
enddo
c-- Compute and remove offset for current tile and sum over all
c-- tiles of this instance.
offset = 0. _d 0
offset_sum = 0. _d 0
c-- Sum over this thread tiles.
do bj = jtlo,jthi
do bi = itlo,ithi
do j = 1,sny
do i = 1,snx
offset = offset +
& mdtmask(i,j,bi,bj)*cosphi(i,j,bi,bj)*
& (mdt(i,j,bi,bj) - psmean(i,j,bi,bj))
offset_sum = offset_sum +
& mdtmask(i,j,bi,bj)*cosphi(i,j,bi,bj)
enddo
enddo
enddo
enddo
c-- Do a global summation.
_GLOBAL_SUM_RL( offset , mythid )
_GLOBAL_SUM_RL( offset_sum , mythid )
if (offset_sum .eq. 0.0) then
_BEGIN_MASTER( mythid )
write(msgbuf,'(a)') ' cost_ssh: offset_sum = zero!'
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
_END_MASTER( mythid )
cph stop ' ... stopped in cost_ssh.'
else
_BEGIN_MASTER( mythid )
write(msgbuf,'(a,d22.15)')
& ' cost_ssh: offset_sum = ',offset_sum
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
_END_MASTER( mythid )
endif
offset = offset / offset_sum
#ifdef ALLOW_PROFILES
if ( usePROFILES ) then
do bj = jtlo,jthi
do bi = itlo,ithi
do j = 1,sny
do i = 1,snx
prof_etan_mean(i,j,bi,bj)=offset+psmean(i,j,bi,bj)
enddo
enddo
enddo
enddo
_EXCH_XY_RL( prof_etan_mean, mythid )
endif
#endif
#ifdef ALLOW_SSH_MEAN_COST_CONTRIBUTION
#ifndef ALLOW_SSH_TOT
c-- ==========
c-- Mean
c-- ==========
c-- compute mean ssh difference cost contribution
objf_hmean = 0. _d 0
do bj = jtlo,jthi
do bi = itlo,ithi
do j = jmin,jmax
do i = imin,imax
junk = psmean(i,j,bi,bj) - mdt(i,j,bi,bj) + offset
objf_hmean = objf_hmean + junk*junk
& * wp(i,j,bi,bj)*mdtmask(i,j,bi,bj)
if ( wp(i,j,bi,bj)*mdtmask(i,j,bi,bj) .ne. 0. )
& num_hmean = num_hmean + 1.D0
enddo
enddo
enddo
enddo
CALL GLOBAL_SUM_R8 ( objf_hmean, mythid )
CALL GLOBAL_SUM_R8 ( num_hmean, mythid )
#endif /* ALLOW_SSH_TOT */
#endif /* ALLOW_SSH_MEAN_COST_CONTRIBUTION */
c-- ==========
c-- Anomalies.
c-- ==========
c-- Loop over records for the second time.
do irec = 1, ndaysrec
call ACTIVE_READ_XY( fname, psbar, irec, doglobalread,
& ladinit, optimcycle, mythid,
& xx_psbar_mean_dummy )
#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
call COST_SLA_READ( topexfile, topexstartdate, topexperiod,
& topexintercept, topexslope,
& tpobs, tpmask,
& irec, mythid )
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
call COST_SLA_READ( ersfile, ersstartdate, ersperiod,
& ersintercept, ersslope,
& ersobs, ersmask,
& irec, mythid )
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
call COST_SLA_READ( gfofile, gfostartdate, gfoperiod,
& gfointercept, gfoslope,
& gfoobs, gfomask,
& irec, mythid )
#endif
do bj = jtlo,jthi
do bi = itlo,ithi
#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
do j = jmin,jmax
do i = imin,imax
c-- The array psobs contains SSH anomalies.
wwwtp(i,j) = wtp(i,j,bi,bj) *cosphi(i,j,bi,bj)
#ifndef ALLOW_SSH_TOT
junk = ((psbar(i,j,bi,bj) - psmean(i,j,bi,bj)) -
& tpobs(i,j,bi,bj))
& *tpmask(i,j,bi,bj)
objf_tp(bi,bj) = objf_tp(bi,bj)
& + junk*junk*wwwtp(i,j)
if ( wwwtp(i,j)*junk .ne. 0. )
& num_tp(bi,bj) = num_tp(bi,bj) + 1. _d 0
#else
if (mdtmask(i,j,bi,bj)*tpmask(i,j,bi,bj)
& *wp(i,j,bi,bj)*wwwtp(i,j) .ne.0.) then
junk = ( psbar(i,j,bi,bj) -
& (tpobs(i,j,bi,bj)+mdt(i,j,bi,bj)-offset) )
objf_tp(bi,bj) = objf_tp(bi,bj)
& +junk*junk/( 1/wp(i,j,bi,bj)+1/wwwtp(i,j) )
num_tp(bi,bj) = num_tp(bi,bj) + 1. _d 0
endif
#endif /* ALLOW_SSH_TOT */
enddo
enddo
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
do j = jmin,jmax
do i = imin,imax
c-- The array ersobs contains SSH anomalies.
wwwers(i,j) = wers(i,j,bi,bj)*cosphi(i,j,bi,bj)
#ifndef ALLOW_SSH_TOT
junk = ((psbar(i,j,bi,bj) - psmean(i,j,bi,bj)) -
& ersobs(i,j,bi,bj))
& *ersmask(i,j,bi,bj)
objf_ers(bi,bj) = objf_ers(bi,bj)
& + junk*junk*wwwers(i,j)
if ( wwwers(i,j)*junk .ne. 0. )
& num_ers(bi,bj) = num_ers(bi,bj) + 1. _d 0
#else
if (mdtmask(i,j,bi,bj)*ersmask(i,j,bi,bj)
& *wp(i,j,bi,bj)*wwwers(i,j) .ne.0.) then
junk = ( psbar(i,j,bi,bj) -
& (ersobs(i,j,bi,bj)+mdt(i,j,bi,bj)-offset) )
objf_ers(bi,bj) = objf_ers(bi,bj)
& +junk*junk/( 1/wp(i,j,bi,bj)+1/wwwers(i,j) )
num_ers(bi,bj) = num_ers(bi,bj) + 1. _d 0
endif
#endif /* ALLOW_SSH_TOT */
enddo
enddo
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
do j = jmin,jmax
do i = imin,imax
c-- The array gfoobs contains SSH anomalies.
wwwgfo(i,j) = wgfo(i,j,bi,bj)*cosphi(i,j,bi,bj)
#ifndef ALLOW_SSH_TOT
junk = ((psbar(i,j,bi,bj) - psmean(i,j,bi,bj)) -
& gfoobs(i,j,bi,bj))
& *gfomask(i,j,bi,bj)
objf_gfo(bi,bj) = objf_gfo(bi,bj)
& + junk*junk*wwwgfo(i,j)
if ( wwwgfo(i,j)*junk .ne. 0. )
& num_gfo(bi,bj) = num_gfo(bi,bj) + 1. _d 0
#else
if (mdtmask(i,j,bi,bj)*gfomask(i,j,bi,bj)
& *wp(i,j,bi,bj)*wwwgfo(i,j) .ne.0.) then
junk = ( psbar(i,j,bi,bj) -
& (gfoobs(i,j,bi,bj)+mdt(i,j,bi,bj)-offset) )
objf_gfo(bi,bj) = objf_gfo(bi,bj)
& +junk*junk/( 1/wp(i,j,bi,bj)+1/wwwgfo(i,j) )
num_gfo(bi,bj) = num_gfo(bi,bj) + 1. _d 0
endif
#endif /* ALLOW_SSH_TOT */
enddo
enddo
#endif
enddo
enddo
enddo
c-- End of second loop over records.
do bj = jtlo,jthi
do bi = itlo,ithi
objf_h(bi,bj) = objf_h(bi,bj) +
& objf_tp(bi,bj) + objf_ers(bi,bj) + objf_gfo(bi,bj)
num_h(bi,bj) = num_h(bi,bj) +
& num_tp(bi,bj) + num_ers(bi,bj) + num_gfo(bi,bj)
enddo
enddo
#endif /* ifdef ALLOW_SSH_COST_CONTRIBUTION */
end