C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_generic.F,v 1.39 2016/03/22 22:29:25 gforget Exp $
C $Name: $
#include "ECCO_OPTIONS.h"
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: cost_generic
C !INTERFACE:
subroutine COST_GENERIC(
& nnzbar, localbarfile, dummy,
& nnzobs, localobsfile, localerrfile,
& mult_local, nrecloc, nrecobs,
& localstartdate, localperiod,
& ylocmask, spminloc, spmaxloc, spzeroloc,
& preproc, preproc_c, preproc_i, preproc_r,
& posproc, posproc_c, posproc_i, posproc_r,
& outlev, outname,
& objf_local, num_local,
& myiter, mytime, mythid )
C !DESCRIPTION: \bv
C Generic routine for evaluating time-dependent
c cost function contribution
C \ev
C !USES:
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
#ifdef ALLOW_SEAICE
# include "SEAICE_COST.h"
#endif
c == routine arguments ==
integer myiter
integer mythid
integer nnzbar, nnzobs
integer nrecloc, nrecobs
integer localstartdate(4)
integer outlev
integer preproc_i(NGENPPROC)
integer posproc_i(NGENPPROC)
_RL objf_local(nsx,nsy)
_RL num_local(nsx,nsy)
_RL dummy
_RL mult_local
_RL mytime
_RL localperiod
_RL spminloc
_RL spmaxloc
_RL spzeroloc
_RL preproc_r(NGENPPROC)
_RL posproc_r(NGENPPROC)
character*(1) ylocmask
character*(MAX_LEN_FNAM) localbarfile
character*(MAX_LEN_FNAM) localobsfile
character*(MAX_LEN_FNAM) localerrfile
character*(MAX_LEN_FNAM) preproc(NGENPPROC)
character*(MAX_LEN_FNAM) preproc_c(NGENPPROC)
character*(MAX_LEN_FNAM) posproc(NGENPPROC)
character*(MAX_LEN_FNAM) posproc_c(NGENPPROC)
character*(MAX_LEN_FNAM) outname
#ifdef ALLOW_ECCO
c == local variables ==
integer bi,bj,k2
integer itlo,ithi
integer jtlo,jthi
logical domean, doanom
_RL localdifmean1 (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
_RL localdifmean2 (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
CEOP
jtlo = mybylo(mythid)
jthi = mybyhi(mythid)
itlo = mybxlo(mythid)
ithi = mybxhi(mythid)
c-- Initialise local variables.
do bj = jtlo,jthi
do bi = itlo,ithi
objf_local(bi,bj) = 0. _d 0
num_local(bi,bj) = 0. _d 0
enddo
enddo
call ECCO_ZERO(localdifmean1,Nr,zeroRL,myThid)
call ECCO_ZERO(localdifmean2,Nr,zeroRL,myThid)
domean=.FALSE.
doanom=.FALSE.
do k2 = 1, NGENPPROC
if (preproc(k2).EQ.'mean') domean=.TRUE.
if (preproc(k2).EQ.'anom') doanom=.TRUE.
enddo
C Extra time loop to compute time-mean fields and costs
if ( (.NOT. ( localobsfile.EQ.' ' ) )
& .AND. ( domean .OR. doanom ) ) then
call COST_GENLOOP(
& localdifmean1,localdifmean2,.FALSE.,
& nnzbar, localbarfile, dummy,
& nnzobs, localobsfile, localerrfile,
& mult_local, nrecloc, nrecobs,
& localstartdate, localperiod,
& ylocmask, spminloc, spmaxloc, spzeroloc,
& preproc, preproc_c, preproc_i, preproc_r,
& posproc, posproc_c, posproc_i, posproc_r,
& outlev, outname,
& objf_local, num_local,
& myiter, mytime, mythid )
endif
call ECCO_ZERO(localdifmean1,Nr,zeroRL,myThid)
if ((.NOT.(localobsfile.EQ.' ')).AND.(.NOT.domean)) then
call COST_GENLOOP(
& localdifmean2,localdifmean1,.TRUE.,
& nnzbar, localbarfile, dummy,
& nnzobs, localobsfile, localerrfile,
& mult_local, nrecloc, nrecobs,
& localstartdate, localperiod,
& ylocmask, spminloc, spmaxloc, spzeroloc,
& preproc, preproc_c, preproc_i, preproc_r,
& posproc, posproc_c, posproc_i, posproc_r,
& outlev, outname,
& objf_local, num_local,
& myiter, mytime, mythid )
endif
#endif /* ALLOW_ECCO */
return
end
C--------------
subroutine COST_GENLOOP(
& localdifmeanIn,localdifmeanOut, addVariaCost,
& nnzbar, localbarfile, dummy,
& nnzobs, localobsfile, localerrfile,
& mult_local, nrecloc, nrecobs,
& localstartdate, localperiod,
& ylocmask, spminloc, spmaxloc, spzeroloc,
& preproc, preproc_c, preproc_i, preproc_r,
& posproc, posproc_c, posproc_i, posproc_r,
& outlev, outname,
& objf_local, num_local,
& myiter, mytime, mythid )
C !DESCRIPTION: \bv
C Generic routine for evaluating time-dependent
c cost function contribution
C \ev
C !USES:
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
#ifdef ALLOW_SEAICE
# include "SEAICE_COST.h"
#endif
c == routine arguments ==
integer myiter
integer mythid
integer nnzbar, nnzobs
integer nrecloc, nrecobs
integer localstartdate(4)
integer outlev
integer preproc_i(NGENPPROC)
integer posproc_i(NGENPPROC)
_RL objf_local(nsx,nsy)
_RL num_local(nsx,nsy)
_RL dummy
_RL mult_local
_RL mytime
_RL localperiod
_RL spminloc
_RL spmaxloc
_RL spzeroloc
_RL preproc_r(NGENPPROC)
_RL posproc_r(NGENPPROC)
character*(1) ylocmask
character*(MAX_LEN_FNAM) localbarfile
character*(MAX_LEN_FNAM) localobsfile
character*(MAX_LEN_FNAM) localerrfile
character*(MAX_LEN_FNAM) preproc(NGENPPROC)
character*(MAX_LEN_FNAM) preproc_c(NGENPPROC)
character*(MAX_LEN_FNAM) posproc(NGENPPROC)
character*(MAX_LEN_FNAM) posproc_c(NGENPPROC)
character*(MAX_LEN_FNAM) outname
logical addVariaCost
_RL localdifmeanIn (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
_RL localdifmeanOut (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
#ifdef ALLOW_ECCO
c == local variables ==
integer bi,bj
integer itlo,ithi
integer jtlo,jthi
integer irec, jrec
integer il, k2
integer localrec, obsrec
integer nrecloop, nrecclim, k2smooth
logical domean, doanom, dovarwei, doclim, dosmooth, dosumsq
_RL localmask (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
_RL localbar (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 localobs (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 localdifmsk (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
_RL localdifsum (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
_RL localdifnum (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
_RL fac
character*(128) fname1, fname2, fname3
character*200 msgbuf
logical exst
c == external functions ==
integer ilnblnk
external
CEOP
call ECCO_ZERO(localbar,Nr,zeroRL,myThid)
call ECCO_ZERO(localweight,Nr,zeroRL,myThid)
call ECCO_ZERO(localtmp,Nr,zeroRL,myThid)
call ECCO_ZERO(localmask,Nr,zeroRL,myThid)
call ECCO_ZERO(localobs,Nr,zeroRL,myThid)
call ECCO_ZERO(localdif,Nr,zeroRL,myThid)
call ECCO_ZERO(difmask,Nr,zeroRL,myThid)
call ECCO_ZERO(localdifmsk,Nr,zeroRL,myThid)
call ECCO_ZERO(localdifsum,Nr,zeroRL,myThid)
call ECCO_ZERO(localdifnum,Nr,zeroRL,myThid)
dosumsq=.TRUE.
domean=.FALSE.
doanom=.FALSE.
dovarwei=.FALSE.
dosmooth=.FALSE.
k2smooth=1
doclim=.FALSE.
nrecclim=nrecloc
fac=oneRL
do k2 = 1, NGENPPROC
if (preproc(k2).EQ.'mean') domean=.TRUE.
if (preproc(k2).EQ.'anom') doanom=.TRUE.
if (preproc(k2).EQ.'variaweight') dovarwei=.TRUE.
if (preproc(k2).EQ.'nosumsq') dosumsq=.FALSE.
if (posproc(k2).EQ.'smooth') then
dosmooth=.TRUE.
k2smooth=k2
endif
if (preproc(k2).EQ.'clim') then
doclim=.TRUE.
nrecclim=preproc_i(k2)
endif
if (preproc(k2).EQ.'factor') then
fac=preproc_r(k2)
endif
enddo
c-- Assign mask
if ( ylocmask .EQ. 'C' .OR. ylocmask .EQ. 'c' ) then
call ECCO_CPRSRL(maskC,nr,localmask,nr,myThid)
elseif ( ylocmask .EQ. 'S' .OR. ylocmask .EQ. 's' ) then
call ECCO_CPRSRL(maskS,nr,localmask,nr,myThid)
elseif ( ylocmask .EQ. 'W' .OR. ylocmask .EQ. 'w' ) then
call ECCO_CPRSRL(maskW,nr,localmask,nr,myThid)
else
STOP 'cost_generic: wrong ylocmask'
endif
c-- set nrecloop to nrecloc
nrecloop=nrecloc
c-- reset nrecloop, if needed, according to preproc
if ( doclim ) nrecloop=MIN(nrecloop,nrecclim)
c-- loop over obsfile records
do irec = 1, nrecloop
c-- load weights
exst=.FALSE.
jrec=1
if( dovarwei ) jrec = irec
call COST_GENCAL(localbarfile, localerrfile,
& jrec, localstartdate, localperiod, fname1,
& fname3, localrec, obsrec, exst, mythid )
call ECCO_ZERO(localweight,nnzobs,zeroRL,myThid)
if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) )
& call ECCO_READWEI(fname3,localweight,localrec,nnzobs,mythid)
c-- determine records and file names
exst=.FALSE.
call COST_GENCAL(localbarfile, localobsfile,
& irec, localstartdate, localperiod, fname1,
& fname2, localrec, obsrec, exst, mythid )
c-- load model average and observed average
call ECCO_ZERO(localbar,nnzbar,zeroRL,myThid)
call COST_GENREAD( fname1, localbar, localtmp, irec, nnzbar,
& nrecloc, preproc, preproc_c, preproc_i, preproc_r,
& dummy, mythid )
call ECCO_MULT(localbar,nnzbar,fac,myThid)
call ECCO_ZERO(localobs,nnzobs,spzeroloc,myThid)
if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) )
& call MDSREADFIELD( fname2, cost_iprec, cost_yftype, nnzobs,
& localobs, localrec, mythid )
c-- Compute masked model-data difference
call ECCO_DIFFMSK( localbar, nnzbar, localobs, nnzobs,
& localmask, spminloc, spmaxloc, spzeroloc,
& localdif, difmask, myThid )
if ( doanom ) call ECCO_SUBTRACT( localdifmeanIn,
& nnzobs, localdif, nnzobs, myThid )
if ( domean.OR.doanom )
& call ECCO_ADDMASK(localdif,difmask, nnzobs,localdifsum,
& localdifnum, nnzobs,myThid)
if (addVariaCost) then
#ifdef ALLOW_SMOOTH
if ( useSMOOTH.AND.dosmooth.AND.
& (nnzbar.EQ.1).AND.(nnzobs.EQ.1) )
& call SMOOTH_HETERO2D(localdif,maskc,
& posproc_c(k2smooth),posproc_i(k2smooth),mythid)
#endif
c-- Compute normalized model-obs cost function
call ECCO_ADDCOST(
I localdif, localweight, difmask, nnzobs, dosumsq,
U objf_local, num_local,
I myThid
& )
c-- output model-data difference to disk
if ( outlev.GT.0 ) then
il=ilnblnk(outname)
write(fname3(1:128),'(2a)') 'misfit_', outname(1:il)
if ( nnzobs.EQ.1 ) CALL
& WRITE_REC_XY_RL( fname3, localdif,irec, eccoiter, mythid )
if ( nnzobs.EQ.nr ) CALL
& WRITE_REC_XYZ_RL( fname3, localdif,irec, eccoiter, mythid )
endif
endif
enddo
c-- End of loop over obsfile records.
call ECCO_ZERO(localdifmeanOut,Nr,zeroRL,myThid)
call ECCO_CP(localdifsum,nnzobs,localdifmeanOut,nnzobs,myThid)
call ECCO_DIVFIELD(localdifmeanOut,nnzobs,localdifnum,myThid)
call ECCO_CP(localdifnum,nnzobs,localdifmsk,nnzobs,myThid)
call ECCO_DIVFIELD(localdifmsk,nnzobs,localdifnum,myThid)
if ( domean ) then
c-- Compute normalized model-obs cost function
call ECCO_ADDCOST(
I localdifmeanOut, localweight, localdifmsk, nnzobs, dosumsq,
U objf_local, num_local, myThid)
c-- output model-data difference to disk
if ( outlev.GT.0 ) then
il=ilnblnk(outname)
write(fname3(1:128),'(2a)') 'misfit_', outname(1:il)
if ( nnzobs.EQ.1 ) CALL
& WRITE_REC_XY_RL(fname3,localdifmeanOut,1,eccoiter,mythid)
if ( nnzobs.EQ.nr ) CALL
& WRITE_REC_XYZ_RL(fname3,localdifmeanOut,1,eccoiter,mythid)
endif
endif
if ( outlev.GT.1 ) then
il=ilnblnk(outname)
write(fname3(1:128),'(2a)') 'weight_', outname(1:il)
if ( nnzobs.EQ.1 ) CALL
& WRITE_REC_XY_RL( fname3, localweight,irec, eccoiter, mythid )
if ( nnzobs.EQ.nr ) CALL
& WRITE_REC_XYZ_RL( fname3, localweight,irec, eccoiter, mythid )
endif
#endif /* ALLOW_ECCO */
return
end