pkg/ecco/cost_generic.F

C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_generic.F,v 1.27 2014/10/22 13:28:27 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, localbar, dummy,
     &     nnzobs, localobsfile, mult_local,
     &     nrecloc, nrecobs,
     &     localstartdate, localperiod,
     &     ylocmask, localweight,
     &     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 localbar   (1-olx:snx+olx,1-oly:sny+oly,nnzbar,nsx,nsy)
      _RL localweight(1-olx:snx+olx,1-oly:sny+oly,nnzobs,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) 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
      integer itlo,ithi
      integer jtlo,jthi
      integer irec
      integer  il
      integer localrec, obsrec
      integer nrecloop

      _RL localmask  (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
      _RL localobs   (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)
      _RL localdif   (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)
      _RL difmask    (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)

      character*(128) fname1, fname2, fname3

      logical exst

c     == external functions ==

      integer  ilnblnk
      external ilnblnk

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(localobs,nnzobs,zeroRL,myThid)
      call ecco_zero(localdif,nnzobs,zeroRL,myThid)
      call ecco_zero(difmask,nnzobs,zeroRL,myThid)

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 ( preproc(1) .EQ. 'clim') nrecloop=MIN(nrecloop,preproc_i(1))

      if ( .NOT. ( localobsfile.EQ.' ' ) ) then

c--   loop over obsfile records
      do irec = 1, nrecloop

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, irec, nnzbar,
     &       nrecloc, preproc, preproc_c, preproc_i, preproc_r, 
     &       dummy, 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 )

#ifdef ALLOW_SMOOTH
      if ( (useSMOOTH).AND.(posproc(1).EQ.'smooth') )
     &  call smooth_hetero2d(localdif,maskc,
     &     posproc_c(1),posproc_i(1),mythid)
#endif

c--     Compute normalized model-obs cost function
        call ecco_addcost(
     I                   localdif, localweight, difmask, nnzobs,
     I                   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

      enddo
c--   End of loop over obsfile records.

      endif

#endif /* ALLOW_ECCO */

      end
1	gforget	1.28	C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_generic.F,v 1.27 2014/10/22 13:28:27 gforget Exp $
2	jmc	1.8	C $Name: $
3	heimbach	1.1
4	jmc	1.18	#include "ECCO_OPTIONS.h"
5	heimbach	1.1
6	gforget	1.25	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
7			CBOP
8			C !ROUTINE: cost_generic
9			C !INTERFACE:
10	jmc	1.8	subroutine cost_generic(
11	gforget	1.26	& nnzbar, localbarfile, localbar, dummy,
12	heimbach	1.9	& nnzobs, localobsfile, mult_local,
13	gforget	1.21	& nrecloc, nrecobs,
14			& localstartdate, localperiod,
15	heimbach	1.10	& ylocmask, localweight,
16	gforget	1.21	& spminloc, spmaxloc, spzeroloc,
17	gforget	1.24	& preproc, preproc_c, preproc_i, preproc_r,
18			& posproc, posproc_c, posproc_i, posproc_r,
19	gforget	1.22	& outlev, outname,
20	heimbach	1.1	& objf_local, num_local,
21			& myiter, mytime, mythid )
22
23	gforget	1.25	C !DESCRIPTION: \bv
24			C Generic routine for evaluating time-dependent
25	heimbach	1.1	c cost function contribution
26	gforget	1.25	C \ev
27	heimbach	1.1
28	gforget	1.25	C !USES:
29	heimbach	1.1	implicit none
30
31			c == global variables ==
32
33			#include "EEPARAMS.h"
34			#include "SIZE.h"
35			#include "PARAMS.h"
36	heimbach	1.10	#include "GRID.h"
37	heimbach	1.1	#ifdef ALLOW_CAL
38			# include "cal.h"
39			#endif
40	gforget	1.19	#ifdef ALLOW_ECCO
41			# include "ecco.h"
42	gforget	1.17	#endif
43			#ifdef ALLOW_SEAICE
44			# include "SEAICE_COST.h"
45	heimbach	1.1	#endif
46
47			c == routine arguments ==
48
49	gforget	1.24	integer myiter
50			integer mythid
51	gforget	1.21	integer nnzbar, nnzobs
52			integer nrecloc, nrecobs
53	heimbach	1.1	integer localstartdate(4)
54	gforget	1.21	integer outlev
55	gforget	1.24	integer preproc_i(NGENPPROC)
56			integer posproc_i(NGENPPROC)
57	heimbach	1.1
58	gforget	1.24	_RL objf_local(nsx,nsy)
59			_RL num_local(nsx,nsy)
60	heimbach	1.1	_RL localbar (1-olx:snx+olx,1-oly:sny+oly,nnzbar,nsx,nsy)
61			_RL localweight(1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)
62	gforget	1.26	_RL dummy
63	heimbach	1.2	_RL mult_local
64	heimbach	1.1	_RL mytime
65			_RL localperiod
66			_RL spminloc
67			_RL spmaxloc
68			_RL spzeroloc
69	gforget	1.24	_RL preproc_r(NGENPPROC)
70			_RL posproc_r(NGENPPROC)
71	heimbach	1.1
72	heimbach	1.10	character*(1) ylocmask
73	heimbach	1.1	character*(MAX_LEN_FNAM) localbarfile
74			character*(MAX_LEN_FNAM) localobsfile
75	gforget	1.24	character*(MAX_LEN_FNAM) preproc(NGENPPROC)
76			character*(MAX_LEN_FNAM) preproc_c(NGENPPROC)
77			character*(MAX_LEN_FNAM) posproc(NGENPPROC)
78			character*(MAX_LEN_FNAM) posproc_c(NGENPPROC)
79	gforget	1.22	character*(MAX_LEN_FNAM) outname
80	gforget	1.21
81			#ifdef ALLOW_ECCO
82	heimbach	1.1
83			c == local variables ==
84
85			integer bi,bj
86			integer itlo,ithi
87			integer jtlo,jthi
88			integer irec
89			integer il
90	gforget	1.23	integer localrec, obsrec
91	gforget	1.28	integer nrecloop
92	heimbach	1.1
93	heimbach	1.10	_RL localmask (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
94	heimbach	1.9	_RL localobs (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)
95	gforget	1.22	_RL localdif (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)
96	gforget	1.25	_RL difmask (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)
97	heimbach	1.1
98	gforget	1.22	character*(128) fname1, fname2, fname3
99	heimbach	1.1
100			logical exst
101
102			c == external functions ==
103
104			integer ilnblnk
105			external ilnblnk
106
107	gforget	1.25	CEOP
108	heimbach	1.1
109			jtlo = mybylo(mythid)
110			jthi = mybyhi(mythid)
111			itlo = mybxlo(mythid)
112			ithi = mybxhi(mythid)
113
114			c-- Initialise local variables.
115
116	heimbach	1.2	do bj = jtlo,jthi
117			do bi = itlo,ithi
118			objf_local(bi,bj) = 0. _d 0
119			num_local(bi,bj) = 0. _d 0
120			enddo
121			enddo
122	gforget	1.23	call ecco_zero(localobs,nnzobs,zeroRL,myThid)
123			call ecco_zero(localdif,nnzobs,zeroRL,myThid)
124			call ecco_zero(difmask,nnzobs,zeroRL,myThid)
125	heimbach	1.2
126	heimbach	1.10	c-- Assign mask
127			if ( ylocmask .EQ. 'C' .OR. ylocmask .EQ. 'c' ) then
128	gforget	1.26	call ecco_cprsrl(maskC,nr,localmask,nr,myThid)
129	heimbach	1.10	elseif ( ylocmask .EQ. 'S' .OR. ylocmask .EQ. 's' ) then
130	gforget	1.26	call ecco_cprsrl(maskS,nr,localmask,nr,myThid)
131	heimbach	1.10	elseif ( ylocmask .EQ. 'W' .OR. ylocmask .EQ. 'w' ) then
132	gforget	1.26	call ecco_cprsrl(maskW,nr,localmask,nr,myThid)
133	heimbach	1.10	else
134			STOP 'cost_generic: wrong ylocmask'
135			endif
136
137	gforget	1.28	c-- set nrecloop to nrecloc
138			nrecloop=nrecloc
139
140			c-- reset nrecloop, if needed, according to preproc
141			if ( preproc(1) .EQ. 'clim') nrecloop=MIN(nrecloop,preproc_i(1))
142	gforget	1.20
143	heimbach	1.7	if ( .NOT. ( localobsfile.EQ.' ' ) ) then
144	heimbach	1.2
145	gforget	1.22	c-- loop over obsfile records
146	gforget	1.28	do irec = 1, nrecloop
147	heimbach	1.1
148	gforget	1.23	c-- determine records and file names
149			exst=.FALSE.
150			call cost_gencal(localbarfile, localobsfile,
151			& irec, localstartdate, localperiod, fname1,
152			& fname2, localrec, obsrec, exst, mythid )
153
154			c-- load model average and observed average
155	gforget	1.28	call ecco_zero(localbar,nnzbar,zeroRL,myThid)
156	gforget	1.26	call cost_genread( fname1, localbar, irec, nnzbar,
157	gforget	1.28	& nrecloc, preproc, preproc_c, preproc_i, preproc_r,
158			& dummy, mythid )
159	heimbach	1.1
160	gforget	1.23	call ecco_zero(localobs,nnzobs,spzeroloc,myThid)
161			if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) )
162			& call mdsreadfield( fname2, cost_iprec, cost_yftype, nnzobs,
163	heimbach	1.3	& localobs, localrec, mythid )
164	heimbach	1.1
165	gforget	1.23	c-- Compute masked model-data difference
166	gforget	1.26	call ecco_diffmsk( localbar, nnzbar, localobs, nnzobs,
167	gforget	1.23	& localmask, spminloc, spmaxloc, spzeroloc,
168	gforget	1.25	& localdif, difmask, myThid )
169	gforget	1.22
170	gforget	1.27	#ifdef ALLOW_SMOOTH
171			if ( (useSMOOTH).AND.(posproc(1).EQ.'smooth') )
172			& call smooth_hetero2d(localdif,maskc,
173			& posproc_c(1),posproc_i(1),mythid)
174			#endif
175
176	gforget	1.22	c-- Compute normalized model-obs cost function
177	gforget	1.25	call ecco_addcost(
178			I localdif, localweight, difmask, nnzobs,
179			I objf_local, num_local,
180			I myThid
181			& )
182	gforget	1.22
183	gforget	1.23	c-- output model-data difference to disk
184			if ( outlev.GT.0 ) then
185			il=ilnblnk(outname)
186			write(fname3(1:128),'(2a)') 'misfit_', outname(1:il)
187			if ( nnzobs.EQ.1 ) CALL
188			& WRITE_REC_XY_RL( fname3, localdif,irec, eccoiter, mythid )
189			if ( nnzobs.EQ.nr ) CALL
190			& WRITE_REC_XYZ_RL( fname3, localdif,irec, eccoiter, mythid )
191			endif
192	gforget	1.22
193	heimbach	1.1	enddo
194	gforget	1.22	c-- End of loop over obsfile records.
195	heimbach	1.1
196	heimbach	1.2	endif
197
198	gforget	1.21	#endif /* ALLOW_ECCO */
199
200	heimbach	1.1	end