C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_gencost_transp.F,v 1.4 2015/11/12 13:22:14 atn Exp $
C $Name: $
#include "ECCO_OPTIONS.h"
subroutine COST_GENCOST_TRANSP(mythid)
c ==================================================================
c SUBROUTINE cost_gencost_transp
c ==================================================================
c
c o Evaluate cost function contributions of section transport.
c
c ==================================================================
c SUBROUTINE cost_gencost_boxmean
c ==================================================================
implicit none
c == global variables ==
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"
#ifdef ALLOW_CAL
# include "cal.h"
#endif
#ifdef ALLOW_ECCO
# include "ecco.h"
#endif
c == routine arguments ==
integer mythid
#ifdef ALLOW_GENCOST_CONTRIBUTION
c == local variables ==
integer nnzobs, nnzbar
parameter (nnzbar = Nr, nnzobs = Nr)
integer nrecloc, localrec
integer localstartdate(4)
_RL myobs (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
_RL mybar (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
_RL localdif (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
_RL difmask (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
_RL localweight(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
_RL localtmp (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
_RL dummy,facCost,facNum
c- facNum is 1 divided by the number of tiles in SIZE dot h
parameter ( facNum = 1./float(nPx) )
_RL localperiod
_RL spminloc, spmaxloc, spzeroloc
_RL tmpMeanTile(nSx,nSy),tmpNumTile(nSx,nSy)
_RL tmpMeanGlo,tmpNumGlo
character*(MAX_LEN_FNAM) mybarfile
character*(MAX_LEN_FNAM) myobsfile
integer kgen(NGENCOST3D)
integer bi,bj
integer i,j,k
integer itlo,ithi
integer jtlo,jthi
integer obsrec,irec,jrec
integer il,k2
integer icount,icount_transp
logical dosumsq, dovarwei, doreadobs
integer preproc_i(NGENPPROC)
_RL preproc_r(NGENPPROC)
character*(MAX_LEN_FNAM) preproc(NGENPPROC)
character*(MAX_LEN_FNAM) preproc_c(NGENPPROC)
logical doglobalread
logical ladinit
character*(MAX_LEN_MBUF) msgbuf
character*(128) fname1, fname0
logical exst
c == external functions ==
integer ilnblnk
external
c == end of interface ==
jtlo = mybylo(mythid)
jthi = mybyhi(mythid)
itlo = mybxlo(mythid)
ithi = mybxhi(mythid)
c=============== PART 0: initilization ===================
c-- detect the relevant gencost indices
do k=1,NGENCOST3D
kgen(k)=0
enddo
c-- write a report of how many transport costs
write(msgbuf,'(A)') 'Inside cost_gencost_transp:'
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT, myThid )
icount_transp=0
do k=1,NGENCOST
if ( (gencost_name(k)(1:6).EQ.'transp').AND.
& (using_gencost(k)) ) then
icount_transp=icount_transp+1
kgen(icount_transp)=k
il=ilnblnk(gencost_barfile(kgen(icount_transp)))
write(msgbuf,'(A,i4,A,A)') 'Cost ',kgen(icount_transp),
& ': ',gencost_barfile(kgen(icount_transp))(1:il)
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT, myThid )
endif
enddo
c-- write a report of how many transport costs
do icount=1,icount_transp
if (kgen(icount).NE.0) then
c ========
c-- initialize objf and num:
DO bj=jtlo,jthi
DO bi=itlo,ithi
objf_gencost(bi,bj,kgen(icount))=0. _d 0
num_gencost(bi,bj,kgen(icount))=0. _d 0
ENDDO
ENDDO
c-- Initialise local variables.
nrecloc=0
nrecloc=gencost_nrec(kgen(icount))
c-- only enters if there is at least 1 record
if(nrecloc.gt.0) then
facCost=1. _d 0 / float(nrecloc)
localperiod=0.
localperiod=gencost_period(kgen(icount))
dummy=gencost_dummy(kgen(icount))
spminloc=gencost_spmin(kgen(icount))
spmaxloc=gencost_spmax(kgen(icount))
spzeroloc=gencost_spzero(kgen(icount))
c prefer to have preproc match nosumsq but can not seem to get syntax
c to work for comparison of characters so match dosumsq for now.
dosumsq=.FALSE.
dovarwei=.FALSE.
do k2 = 1, NGENPPROC
preproc(k2)=gencost_preproc(k2,kgen(icount))
preproc_i(k2)=gencost_preproc_i(k2,kgen(icount))
preproc_c(k2)=gencost_preproc_c(k2,kgen(icount))
preproc_r(k2)=gencost_preproc_r(k2,kgen(icount))
if (preproc(k2).EQ.'variaweight') dovarwei=.TRUE.
if (preproc(k2)(1:7).EQ.'dosumsq') dosumsq=.TRUE.
enddo
c-- report of dosumsq flag to make sure it is false
il=ilnblnk(gencost_name(kgen(icount)))
write(msgbuf,'(3A,L5,2A)')
& 'Cost ',gencost_name(kgen(icount))(1:il),
& ' dosumsq: ',dosumsq,' preproc(1): ',preproc(1)(1:7)
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT, myThid )
c-- set mybarfile
mybarfile=gencost_barfile(kgen(icount))
c-- set obsfile if defined. Not use for now. send warning
doreadobs=.FALSE.
if( .not. gencost_datafile(kgen(icount)).eq.' ') then
c doreadobs=.TRUE.
write(msgbuf,'(A)')
& '**WARNING** S/R COST_GENCOST_TRANSP: gencost_datafile '
CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
& SQUEEZE_RIGHT, myThid )
write(msgbuf,'(A)')
& 'are currently ignored. Adjust the S/R to add code.'
CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
& SQUEEZE_RIGHT, myThid )
endif
c model mask[W,S]: already included in transp calc in ecco_phys
c--------- PART 0.1 read weights --------------------
c-- read mask in errfile, already stored under msktrVol[W,S]
c-- for now, assume non-time-variable mask
c=============== PART 1: main loop ===================
do irec = 1,nrecloc
c--------- PART 1.1 read barfiles ------------------
c-- set all bars to zeros:
call ECCO_ZERO(mybar,Nr,zeroRL,myThid)
c gencost_errfile and fname0 are dummy, get fname1 from mybarfile
exst=.FALSE.
call ECCO_ZERO(localtmp,Nr,zeroRL,myThid)
call COST_GENCAL(mybarfile,gencost_errfile(kgen(icount)),
& irec,localstartdate,localperiod,fname1,
& fname0,localrec,obsrec,exst,mythid)
call COST_GENREAD(fname1,mybar,localtmp,irec,nnzbar,
& nrecloc,preproc,preproc_c,preproc_i,preproc_r,
& dummy,mythid)
c--------- PART 1.2 read data --------------------
c-- ignore for now, but use doreadobs flag if needed
c-- be careful of recomputation when put inside if-end block
c if(doreadobs) then
call ECCO_ZERO(myobs,Nr,zeroRL,myThid)
c endif
c--------- PART 1.3 Cost calculation -------------
c-- keep total at each irec to print out for time-series
DO bj = jtlo,jthi
DO bi = itlo,ithi
tmpMeanTile(bi,bj) = 0. _d 0
tmpNumTile(bi,bj) = 0. _d 0
ENDDO
ENDDO
c compute obs minus bar (localdif) and mask (difmask)
c note localtmp is set to 1.
call ECCO_ZERO(localtmp,nnzobs,oneRL,mythid)
call ECCO_ZERO(localdif,nnzobs,zeroRL,mythid)
call ECCO_ZERO(difmask,nnzobs,zeroRL,mythid)
c take care to set sp[min,max,zero]loc carefully to not
c filter out signal. Can consider skip diffmsk step
call ECCO_DIFFMSK(
I mybar, nnzbar, myobs, nnzobs, localtmp,
I spminloc, spmaxloc, spzeroloc,
O localdif, difmask,
I myThid )
call ECCO_ZERO(localtmp,nnzobs,oneRL,mythid)
call ECCO_ADDCOST(
I localdif,localtmp,difmask,nnzobs,dosumsq,
O tmpMeanTile,tmpNumTile,
I mythid)
cc either use ecco_diffmsk and ecco_addcost from above or simplify
cc to call below. For now keep syntax consistent with gencost
c call ecco_zero(localtmp,nnzobs,oneRL,mythid)
c call ecco_addcost(
c I mybar,localtmp,localtmp,nnzobs,dosumsq,
c O tmpMeanTile,tmpNumTile,
c I mythid)
c global sums for display of time series
c note tmpNumGlo is the constant total wet points in the msktrVol[W,S]
tmpMeanGlo = 0. _d 0
tmpNumGlo = 0. _d 0
il=ilnblnk(gencost_barfile(kgen(icount)))
CALL GLOBAL_SUM_TILE_RL( tmpMeanTile, tmpMeanGlo, myThid )
CALL GLOBAL_SUM_TILE_RL( tmpNumTile, tmpNumGlo, myThid )
WRITE(msgBuf,'(2A,I3,A,1PE21.14,1PE21.14)')
& 'globalsum transp ',gencost_barfile(kgen(icount))(1:il),
& irec,' ',tmpMeanGlo,tmpNumGlo
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT, myThid )
c sum that is actually be used in cost function
DO bj = jtlo,jthi
DO bi = itlo,ithi
objf_gencost(bi,bj,kgen(icount))=
& objf_gencost(bi,bj,kgen(icount))+tmpMeanTile(bi,bj)
ENDDO
ENDDO
enddo !irec
c-- last step:
c-- divide by number of record to get mean transport:
c-- make num_gencost equals number of months/days used
do bj = jtlo,jthi
do bi = itlo,ithi
objf_gencost(bi,bj,kgen(icount))=
& objf_gencost(bi,bj,kgen(icount))*facCost
num_gencost(bi,bj,kgen(icount))=nrecloc*facNum
enddo
enddo
endif !if (nrecloc.gt.0)
endif !if (kgen.NE.0)
enddo !icount_transp
#endif /* ALLOW_GENCOST_CONTRIBUTION */
RETURN
end