C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_gencost_bpv4.F,v 1.9 2017/04/03 23:16:38 ou.wang Exp $
C $Name: $
#include "ECCO_OPTIONS.h"
subroutine COST_GENCOST_BPV4(
I mythid
& )
c ==================================================================
c SUBROUTINE cost_gencost_bpv4
c ==================================================================
c
c o Evaluate cost function contribution of bottom pressure anoamlies
c => GRACE data
c
c started: Gael Forget Oct-2009
c
c ==================================================================
c SUBROUTINE cost_bp
c ==================================================================
implicit none
c == global variables ==
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#ifdef ALLOW_ECCO
# include "ecco.h"
#endif
c == routine arguments ==
integer mythid
#ifdef ALLOW_ECCO
#ifdef ALLOW_GENCOST_CONTRIBUTION
c == local variables ==
integer bi,bj
integer i,j
integer itlo,ithi
integer jtlo,jthi
integer irec
integer il
logical doglobalread
logical ladinit
_RL locbpbar(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
_RL locbpdat(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
_RL locbpmask(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
_RL locwbp(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
_RL bpdifmean ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
_RL bpdifanom ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
_RL bpdatmean ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
_RL bpdatanom ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
_RL bpcount ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
_RL junk
character*(80) fname
character*(80) fname4test
_RL fac
_RL offset
_RL offset_sum
integer k, kgen
c == external functions ==
integer ilnblnk
external
c == end of interface ==
jtlo = mybylo(mythid)
jthi = mybyhi(mythid)
itlo = mybxlo(mythid)
ithi = mybxhi(mythid)
kgen=0
do k=1,NGENCOST
if ( (gencost_name(k).EQ.'bpv4-grace').AND.
& (.NOT.gencost_is1d(k)).AND.
& (using_gencost(k)) ) kgen=k
enddo
if (kgen.GT.0) then
call ECCO_ZERO(gencost_weight(:,:,1,1,kgen),1,zeroRL,myThid)
if ( gencost_errfile(kgen) .NE. ' ' )
& call ECCO_READWEI(gencost_errfile(kgen),
& gencost_weight(:,:,1,1,kgen),1,1,mythid)
c-- initialise local variables
cgf convert phibot from m2/s2 to cm
fac = 1. _d 2 / 9.81 _d 0
do bj = jtlo,jthi
do bi = itlo,ithi
do j = 1,sny
do i = 1,snx
bpdifmean(i,j,bi,bj) = 0. _d 0
bpdifanom(i,j,bi,bj) = 0. _d 0
bpdatmean(i,j,bi,bj) = 0. _d 0
bpdatanom(i,j,bi,bj) = 0. _d 0
bpcount(i,j,bi,bj) = 0. _d 0
locwbp(i,j,bi,bj) = 0. _d 0
locbpbar(i,j,bi,bj) = 0. _d 0
locbpdat(i,j,bi,bj) = 0. _d 0
locbpmask(i,j,bi,bj) = 0. _d 0
enddo
enddo
enddo
enddo
doglobalread = .false.
ladinit = .false.
c-- map global variable to local variables
do bj = jtlo,jthi
do bi = itlo,ithi
do j = 1,sny
do i = 1,snx
locwbp(i,j,bi,bj) = gencost_weight(i,j,bi,bj,kgen)
enddo
enddo
enddo
enddo
c--
#ifdef ALLOW_CTRL
write(fname(1:80),'(80a)') ' '
il=ilnblnk( gencost_barfile(kgen) )
write(fname(1:80),'(2a,i10.10)')
& gencost_barfile(kgen)(1:il),'.',eccoiter
#endif
c-- ============
c-- Mean values.
c-- ============
do irec = 1, nmonsrec
c-- Compute the mean over all bpdat records.
#ifdef ALLOW_AUTODIFF
call ACTIVE_READ_XY( fname, locbpbar, irec, doglobalread,
& ladinit, eccoiter, mythid,
& gencost_dummy(kgen) )
#else
CALL READ_REC_XY_RL( fname, locbpbar,
& iRec, 1, myThid )
#endif
call COST_BP_READ( gencost_datafile(kgen),
& gencost_startdate(1,kgen),
& locbpdat, locbpmask, irec, mythid )
do bj = jtlo,jthi
do bi = itlo,ithi
do j = 1,sny
do i = 1,snx
if ( (locbpmask(i,j,bi,bj).NE. 0. _d 0).AND.
& (maskc(i,j,1,bi,bj).NE. 0. _d 0) ) then
bpdifmean(i,j,bi,bj) = bpdifmean(i,j,bi,bj) +
& ( fac*locbpbar(i,j,bi,bj) - locbpdat(i,j,bi,bj) )
bpdatmean(i,j,bi,bj) = bpdatmean(i,j,bi,bj) +
& locbpdat(i,j,bi,bj)
bpcount(i,j,bi,bj) = bpcount(i,j,bi,bj) + 1. _d 0
endif
enddo
enddo
enddo
enddo
enddo
do bj = jtlo,jthi
do bi = itlo,ithi
do j = 1,sny
do i = 1,snx
if (bpcount(i,j,bi,bj).GT. 0. _d 0) then
bpdifmean(i,j,bi,bj) =
& bpdifmean(i,j,bi,bj)/bpcount(i,j,bi,bj)
bpdatmean(i,j,bi,bj) =
& bpdatmean(i,j,bi,bj)/bpcount(i,j,bi,bj)
endif
enddo
enddo
enddo
enddo
c-- ==========
c-- Anomalies.
c-- ==========
c-- Loop over records for the second time.
do irec = 1, nmonsrec
#ifdef ALLOW_AUTODIFF
call ACTIVE_READ_XY( fname, locbpbar, irec, doglobalread,
& ladinit, eccoiter, mythid,
& gencost_dummy(kgen) )
#else
CALL READ_REC_XY_RL( fname, locbpbar,
& iRec, 1, myThid )
#endif
call COST_BP_READ( gencost_datafile(kgen),
& gencost_startdate(1,kgen),
& locbpdat, locbpmask, irec, mythid )
c-- Compute field of anomalies
do bj = jtlo,jthi
do bi = itlo,ithi
do j = 1,sny
do i = 1,snx
if ( (locbpmask(i,j,bi,bj).NE. 0. _d 0).AND.
& (maskc(i,j,1,bi,bj).NE. 0. _d 0) ) then
bpdifanom(i,j,bi,bj) =
& ( fac*locbpbar(i,j,bi,bj) - locbpdat(i,j,bi,bj) )
& - bpdifmean(i,j,bi,bj)
bpdatanom(i,j,bi,bj) =
& locbpdat(i,j,bi,bj) - bpdatmean(i,j,bi,bj)
else
bpdifanom(i,j,bi,bj) = 0. _d 0
bpdatanom(i,j,bi,bj) = 0. _d 0
endif
enddo
enddo
enddo
enddo
c-- Remove global mean value
offset = 0. _d 0
offset_sum = 0. _d 0
do bj = jtlo,jthi
do bi = itlo,ithi
do j = 1,sny
do i = 1,snx
if ( (locbpmask(i,j,bi,bj).NE. 0. _d 0).AND.
& (maskc(i,j,1,bi,bj).NE. 0. _d 0) ) then
offset = offset + RA(i,j,bi,bj)*bpdifanom(i,j,bi,bj)
offset_sum = offset_sum + RA(i,j,bi,bj)
endif
enddo
enddo
enddo
enddo
_GLOBAL_SUM_RL( offset , mythid )
_GLOBAL_SUM_RL( offset_sum , mythid )
do bj = jtlo,jthi
do bi = itlo,ithi
do j = 1,sny
do i = 1,snx
if ((offset_sum.GT. 0. _d 0).AND.
& (locbpmask(i,j,bi,bj).NE. 0. _d 0).AND.
& (maskc(i,j,1,bi,bj).NE. 0. _d 0)) then
bpdifanom(i,j,bi,bj) = bpdifanom(i,j,bi,bj)
& - offset/offset_sum
endif
enddo
enddo
enddo
enddo
c-- Smooth field of anomalies
if (gencost_outputlevel(kgen).GT.0) then
write(fname4test(1:80),'(1a)') 'bpdifanom_raw'
call MDSWRITEFIELD(fname4test,32,.false.,'RL',
& 1,bpdifanom,irec,1,mythid)
write(fname4test(1:80),'(1a)') 'bpdatanom_raw'
call MDSWRITEFIELD(fname4test,32,.false.,'RL',
& 1,bpdatanom,irec,1,mythid)
endif
#ifdef ALLOW_SMOOTH
if ( useSMOOTH )
& call SMOOTH_BASIC2D(bpdifanom,maskc,300000. _d 0,3000,mythid)
#endif
if (gencost_outputlevel(kgen).GT.0) then
#ifdef ALLOW_SMOOTH
if ( useSMOOTH )
& call SMOOTH_BASIC2D(bpdatanom,maskc,300000. _d 0,3000,mythid)
#endif
write(fname4test(1:80),'(1a)') 'bpdifanom_smooth'
call MDSWRITEFIELD(fname4test,32,.false.,'RL',
& 1,bpdifanom,irec,1,mythid)
write(fname4test(1:80),'(1a)') 'bpdatanom_smooth'
call MDSWRITEFIELD(fname4test,32,.false.,'RL',
& 1,bpdatanom,irec,1,mythid)
endif
c-- Compute cost function
do bj = jtlo,jthi
do bi = itlo,ithi
do j = 1,sny
do i = 1,snx
c-- map to global cost variables
if ( (locwbp(i,j,bi,bj).NE. 0. _d 0).AND.
& (locbpmask(i,j,bi,bj).NE. 0. _d 0).AND.
& (maskc(i,j,1,bi,bj).NE. 0. _d 0) ) then
junk = bpdifanom(i,j,bi,bj)
objf_gencost(kgen,bi,bj) = objf_gencost(kgen,bi,bj)
& + junk*junk*locwbp(i,j,bi,bj)
num_gencost(kgen,bi,bj) = num_gencost(kgen,bi,bj)
& + 1. _d 0
endif
enddo
enddo
enddo
enddo
enddo
endif !if (kgen.GT.0) then
#endif /* ifdef ALLOW_GENCOST_CONTRIBUTION */
#endif /* ifdef ALLOW_ECCO */
end