pkg/ecco/cost_generic.F

C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_generic.F,v 1.14 2012/02/28 00:51:07 gforget Exp $
C $Name:  $

#include "COST_CPPOPTIONS.h"


      subroutine cost_generic(
     &     nnzbar, localbarfile, localbar, xx_localbar_mean_dummy,
     &     nnzobs, localobsfile, mult_local,
     &     nrecloc, localstartdate, localperiod,
     &     ylocmask, localweight,
     &     spminloc, spmaxloc, spzeroloc,
     &     objf_local, num_local,
     &     myiter, mytime, mythid )

c     ==================================================================
c     SUBROUTINE cost_generic
c     ==================================================================
c
c     o Generic routine for evaluating time-dependent
c       cost function contribution
c
c     ==================================================================
c     SUBROUTINE cost_generic
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_COST
# include "ecco_cost.h"
# include "optim.h"
# ifdef ALLOW_SEAICE
#  include "SEAICE_COST.h"
# endif
#endif

c     == routine arguments ==

      integer nnzbar
      integer nnzobs
      integer nrecloc
      integer myiter
      integer mythid
      integer localstartdate(4)

      _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 xx_localbar_mean_dummy
      _RL mult_local
      _RL mytime
      _RL localperiod
      _RL spminloc
      _RL spmaxloc
      _RL spzeroloc
      _RL objf_local(nsx,nsy)
      _RL num_local(nsx,nsy)

      character*(1) ylocmask
      character*(MAX_LEN_FNAM) localbarfile
      character*(MAX_LEN_FNAM) localobsfile

#ifdef ALLOW_COST
c     == local variables ==

      integer bi,bj
      integer i,j,k
      integer itlo,ithi
      integer jtlo,jthi
      integer jmin,jmax
      integer imin,imax
      integer irec
      integer  il
      integer localrec
      integer obsrec

      logical doglobalread
      logical ladinit

      _RL spval
      parameter (spval = -9999. )
      _RL localwww
      _RL localcost
      _RL junk

      _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 cmask (1-olx:snx+olx,1-oly:sny+oly,nnzobs)

      character*(128) fname1, fname2
      character*(MAX_LEN_MBUF) msgbuf

cnew(
      _RL daytime
      _RL diffsecs
      integer dayiter
      integer daydate(4)
      integer difftime(4)
      integer tempDate_1
      integer middate(4)
      integer yday, ymod
      integer md, dd, sd, ld, wd
      integer mody, modm
      integer beginmodel, beginlocal
      logical exst
cnew)

c     == external functions ==

      integer  ilnblnk
      external ilnblnk

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.

      localwww = 0. _d 0

      do bj = jtlo,jthi
        do bi = itlo,ithi
          objf_local(bi,bj) = 0. _d 0
          num_local(bi,bj) = 0. _d 0
          do k = 1,nnzobs
            do j = jmin,jmax
              do i = imin,imax
                localobs(i,j,k,bi,bj) = 0. _d 0
              enddo
            enddo
          enddo
        enddo
      enddo

c--   Assign mask
      do bj = jtlo,jthi
        do bi = itlo,ithi
          do k = 1,Nr
            do j = 1-oly,sny+oly
              do i = 1-olx,snx+olx
      if ( ylocmask .EQ. 'C' .OR. ylocmask .EQ. 'c' ) then
         localmask(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
      elseif ( ylocmask .EQ. 'S' .OR. ylocmask .EQ. 's' ) then
         localmask(i,j,k,bi,bj) = maskS(i,j,k,bi,bj)
      elseif ( ylocmask .EQ. 'W' .OR. ylocmask .EQ. 'w' ) then
         localmask(i,j,k,bi,bj) = maskW(i,j,k,bi,bj)
      else
         STOP 'cost_generic: wrong ylocmask'
      endif
              enddo
            enddo
          enddo
        enddo
      enddo

c--   First, read tiled data.
      doglobalread = .false.
      ladinit      = .false.

      write(fname1(1:128),'(80a)') ' '
      il=ilnblnk( localbarfile )
      write(fname1(1:128),'(2a,i10.10)')
     &     localbarfile(1:il),'.',optimcycle

cph      if ( .NOT. ( mult_local.EQ.0. .OR. localobsfile.EQ.' ' ) ) then
      if ( .NOT. ( localobsfile.EQ.' ' ) ) then

c--   Loop over records for the second time.
      do irec = 1, nrecloc

        if ( nnzbar .EQ. 1 ) then
           call active_read_xy( fname1, localbar, irec, doglobalread,
     &                      ladinit, optimcycle, mythid,
     &                      xx_localbar_mean_dummy )
        else
           call active_read_xyz( fname1, localbar, irec, doglobalread,
     &                       ladinit, optimcycle, mythid,
     &                       xx_localbar_mean_dummy )
        endif

cnew(
        if ( localperiod .EQ. 86400. ) then
c-- assume daily fields
           obsrec = irec
           daytime = FLOAT(secondsperday*(irec-1)) + modelstart
           dayiter = hoursperday*(irec-1) + modeliter0
           call cal_getdate( dayiter, daytime, daydate, mythid )
           call cal_convdate( daydate,yday,md,dd,sd,ld,wd,mythid )
           ymod = localstartdate(1)/10000
           if ( ymod .GE. yday ) then
              call cal_FullDate( localstartdate(1), 0, middate, mythid)
           else
              tempDate_1 = yday*10000+100+1
              call cal_FullDate( tempDate_1, 0, middate, mythid)
           endif
           call cal_TimePassed( middate, daydate, difftime, mythid )
           call cal_ToSeconds( difftime, diffsecs, mythid )
c           localrec = floor(diffsecs/localperiod) + 1
           localrec = int(diffsecs/localperiod) + 1
        else
c-- assume monthly fields
           beginlocal = localstartdate(1)/10000
           beginmodel = modelstartdate(1)/10000
           obsrec =
     &           ( beginmodel - beginlocal )*nmonthyear
     &         + ( mod(modelstartdate(1)/100,100)
     &            -mod(localstartdate(1)/100,100) )
     &         + irec
           mody   = modelstartdate(1)/10000
           modm   = modelstartdate(1)/100 - mody*100
           yday   = mody + INT((modm-1+irec-1)/12)
           localrec = 1 + MOD(modm-1+irec-1,12)
        endif

        il=ilnblnk(localobsfile)
        write(fname2(1:128),'(2a,i4)')
     &       localobsfile(1:il), '_', yday
        inquire( file=fname2, exist=exst )
        if ( (.NOT. exst).AND.( localperiod .NE. 86400. ) ) then
           write(fname2(1:128),'(a)') localobsfile(1:il)
           inquire( file=fname2, exist=exst )
#ifndef COST_GENERIC_ASSUME_CYCLIC
c assume we have one big file, one year after the other
           localrec = obsrec
c otherwise assume climatology, used for each year
#endif
        endif

        if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) ) then
          call mdsreadfield( fname2, cost_iprec, cost_yftype, nnzobs,
     &         localobs, localrec, mythid )
        else
          do bj = jtlo,jthi
            do bi = itlo,ithi
              do k = 1,nnzobs
                do j = jmin,jmax
                  do i = imin,imax
                     localobs(i,j,k,bi,bj) = spval
                  enddo
                enddo
              enddo
            enddo
          enddo
        endif
cnew)

        do bj = jtlo,jthi
          do bi = itlo,ithi

            localcost    = 0. _d 0

c--           Determine the mask on weights
            do k = 1,nnzobs
             do j = jmin,jmax
              do i = imin,imax
               cmask(i,j,k) = cosphi(i,j,bi,bj)*localmask(i,j,k,bi,bj)
                if ( localobs(i,j,k,bi,bj) .lt. spminloc .or.
     &               localobs(i,j,k,bi,bj) .gt. spmaxloc .or.
     &               localobs(i,j,k,bi,bj) .eq. spzeroloc ) then
                   cmask(i,j,k) = 0. _d 0
                endif
              enddo
             enddo
            enddo
c--
            do k = 1,nnzobs
             do j = jmin,jmax
              do i = imin,imax
                localwww  = localweight(i,j,k,bi,bj)*cmask(i,j,k)
                junk      = ( localbar(i,j,k,bi,bj) -
     &                        localobs(i,j,k,bi,bj) )
                localcost = localcost + junk*junk*localwww
                if ( localwww .ne. 0. )
     &               num_local(bi,bj) = num_local(bi,bj) + 1. _d 0
              enddo
             enddo
            enddo

            objf_local(bi,bj) = objf_local(bi,bj) + localcost

          enddo
        enddo

      enddo
c--   End of second loop over records.

c--   End of mult_local or localobsfile
      endif

#endif /* ifdef ALLOW_COST */

      end
1	jmc	1.15	C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_generic.F,v 1.14 2012/02/28 00:51:07 gforget Exp $
2	jmc	1.8	C $Name: $
3	heimbach	1.1
4			#include "COST_CPPOPTIONS.h"
5
6
7	jmc	1.8	subroutine cost_generic(
8	heimbach	1.1	& nnzbar, localbarfile, localbar, xx_localbar_mean_dummy,
9	heimbach	1.9	& nnzobs, localobsfile, mult_local,
10	heimbach	1.1	& nrecloc, localstartdate, localperiod,
11	heimbach	1.10	& ylocmask, localweight,
12	heimbach	1.1	& spminloc, spmaxloc, spzeroloc,
13			& objf_local, num_local,
14			& myiter, mytime, mythid )
15
16			c ==================================================================
17			c SUBROUTINE cost_generic
18			c ==================================================================
19			c
20	jmc	1.8	c o Generic routine for evaluating time-dependent
21	heimbach	1.1	c cost function contribution
22			c
23			c ==================================================================
24			c SUBROUTINE cost_generic
25			c ==================================================================
26
27			implicit none
28
29			c == global variables ==
30
31			#include "EEPARAMS.h"
32			#include "SIZE.h"
33			#include "PARAMS.h"
34	heimbach	1.10	#include "GRID.h"
35	heimbach	1.1	#ifdef ALLOW_CAL
36			# include "cal.h"
37			#endif
38			#ifdef ALLOW_COST
39			# include "ecco_cost.h"
40			# include "optim.h"
41	heimbach	1.4	# ifdef ALLOW_SEAICE
42			# include "SEAICE_COST.h"
43			# endif
44	heimbach	1.1	#endif
45
46			c == routine arguments ==
47
48			integer nnzbar
49			integer nnzobs
50			integer nrecloc
51			integer myiter
52			integer mythid
53			integer localstartdate(4)
54
55			_RL localbar (1-olx:snx+olx,1-oly:sny+oly,nnzbar,nsx,nsy)
56			_RL localweight(1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)
57			_RL xx_localbar_mean_dummy
58	heimbach	1.2	_RL mult_local
59	heimbach	1.1	_RL mytime
60			_RL localperiod
61			_RL spminloc
62			_RL spmaxloc
63			_RL spzeroloc
64			_RL objf_local(nsx,nsy)
65			_RL num_local(nsx,nsy)
66
67	heimbach	1.10	character*(1) ylocmask
68	heimbach	1.1	character*(MAX_LEN_FNAM) localbarfile
69			character*(MAX_LEN_FNAM) localobsfile
70
71			#ifdef ALLOW_COST
72			c == local variables ==
73
74			integer bi,bj
75			integer i,j,k
76			integer itlo,ithi
77			integer jtlo,jthi
78			integer jmin,jmax
79			integer imin,imax
80			integer irec
81			integer il
82			integer localrec
83	heimbach	1.3	integer obsrec
84	heimbach	1.1
85			logical doglobalread
86			logical ladinit
87
88	heimbach	1.3	_RL spval
89			parameter (spval = -9999. )
90	heimbach	1.1	_RL localwww
91			_RL localcost
92			_RL junk
93
94	heimbach	1.10	_RL localmask (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
95	heimbach	1.9	_RL localobs (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)
96	jmc	1.8	_RL cmask (1-olx:snx+olx,1-oly:sny+oly,nnzobs)
97	heimbach	1.1
98			character*(128) fname1, fname2
99			character*(MAX_LEN_MBUF) msgbuf
100
101			cnew(
102			_RL daytime
103			_RL diffsecs
104			integer dayiter
105			integer daydate(4)
106			integer difftime(4)
107	jmc	1.15	integer tempDate_1
108	heimbach	1.1	integer middate(4)
109			integer yday, ymod
110			integer md, dd, sd, ld, wd
111			integer mody, modm
112	heimbach	1.3	integer beginmodel, beginlocal
113	heimbach	1.1	logical exst
114			cnew)
115
116			c == external functions ==
117
118			integer ilnblnk
119			external ilnblnk
120
121			c == end of interface ==
122
123			jtlo = mybylo(mythid)
124			jthi = mybyhi(mythid)
125			itlo = mybxlo(mythid)
126			ithi = mybxhi(mythid)
127			jmin = 1
128			jmax = sny
129			imin = 1
130			imax = snx
131
132			c-- Initialise local variables.
133
134			localwww = 0. _d 0
135
136	heimbach	1.2	do bj = jtlo,jthi
137			do bi = itlo,ithi
138			objf_local(bi,bj) = 0. _d 0
139			num_local(bi,bj) = 0. _d 0
140	heimbach	1.9	do k = 1,nnzobs
141			do j = jmin,jmax
142			do i = imin,imax
143			localobs(i,j,k,bi,bj) = 0. _d 0
144			enddo
145			enddo
146			enddo
147	heimbach	1.2	enddo
148			enddo
149
150	heimbach	1.10	c-- Assign mask
151	gforget	1.11	do bj = jtlo,jthi
152			do bi = itlo,ithi
153			do k = 1,Nr
154			do j = 1-oly,sny+oly
155			do i = 1-olx,snx+olx
156	heimbach	1.10	if ( ylocmask .EQ. 'C' .OR. ylocmask .EQ. 'c' ) then
157	gforget	1.11	localmask(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
158	heimbach	1.10	elseif ( ylocmask .EQ. 'S' .OR. ylocmask .EQ. 's' ) then
159	gforget	1.11	localmask(i,j,k,bi,bj) = maskS(i,j,k,bi,bj)
160	heimbach	1.10	elseif ( ylocmask .EQ. 'W' .OR. ylocmask .EQ. 'w' ) then
161	gforget	1.11	localmask(i,j,k,bi,bj) = maskW(i,j,k,bi,bj)
162	heimbach	1.10	else
163			STOP 'cost_generic: wrong ylocmask'
164			endif
165	gforget	1.11	enddo
166			enddo
167			enddo
168			enddo
169			enddo
170	heimbach	1.10
171	heimbach	1.1	c-- First, read tiled data.
172			doglobalread = .false.
173			ladinit = .false.
174
175			write(fname1(1:128),'(80a)') ' '
176			il=ilnblnk( localbarfile )
177	jmc	1.8	write(fname1(1:128),'(2a,i10.10)')
178	heimbach	1.1	& localbarfile(1:il),'.',optimcycle
179
180	heimbach	1.7	cph if ( .NOT. ( mult_local.EQ.0. .OR. localobsfile.EQ.' ' ) ) then
181			if ( .NOT. ( localobsfile.EQ.' ' ) ) then
182	heimbach	1.2
183	heimbach	1.1	c-- Loop over records for the second time.
184			do irec = 1, nrecloc
185
186			if ( nnzbar .EQ. 1 ) then
187			call active_read_xy( fname1, localbar, irec, doglobalread,
188	jmc	1.8	& ladinit, optimcycle, mythid,
189	heimbach	1.1	& xx_localbar_mean_dummy )
190			else
191			call active_read_xyz( fname1, localbar, irec, doglobalread,
192	jmc	1.8	& ladinit, optimcycle, mythid,
193	heimbach	1.1	& xx_localbar_mean_dummy )
194			endif
195
196			cnew(
197			if ( localperiod .EQ. 86400. ) then
198			c-- assume daily fields
199	heimbach	1.4	obsrec = irec
200	gforget	1.13	daytime = FLOAT(secondsperday*(irec-1)) + modelstart
201			dayiter = hoursperday*(irec-1) + modeliter0
202	heimbach	1.1	call cal_getdate( dayiter, daytime, daydate, mythid )
203			call cal_convdate( daydate,yday,md,dd,sd,ld,wd,mythid )
204			ymod = localstartdate(1)/10000
205	gforget	1.14	if ( ymod .GE. yday ) then
206			call cal_FullDate( localstartdate(1), 0, middate, mythid)
207	heimbach	1.1	else
208	jmc	1.15	tempDate_1 = yday*10000+100+1
209			call cal_FullDate( tempDate_1, 0, middate, mythid)
210	heimbach	1.1	endif
211			call cal_TimePassed( middate, daydate, difftime, mythid )
212			call cal_ToSeconds( difftime, diffsecs, mythid )
213	gforget	1.14	c localrec = floor(diffsecs/localperiod) + 1
214	heimbach	1.1	localrec = int(diffsecs/localperiod) + 1
215			else
216			c-- assume monthly fields
217	heimbach	1.3	beginlocal = localstartdate(1)/10000
218			beginmodel = modelstartdate(1)/10000
219	jmc	1.8	obsrec =
220	heimbach	1.3	& ( beginmodel - beginlocal )*nmonthyear
221			& + ( mod(modelstartdate(1)/100,100)
222			& -mod(localstartdate(1)/100,100) )
223			& + irec
224			mody = modelstartdate(1)/10000
225			modm = modelstartdate(1)/100 - mody*100
226			yday = mody + INT((modm-1+irec-1)/12)
227	heimbach	1.1	localrec = 1 + MOD(modm-1+irec-1,12)
228			endif
229
230			il=ilnblnk(localobsfile)
231			write(fname2(1:128),'(2a,i4)')
232			& localobsfile(1:il), '_', yday
233			inquire( file=fname2, exist=exst )
234	gforget	1.14	if ( (.NOT. exst).AND.( localperiod .NE. 86400. ) ) then
235	heimbach	1.1	write(fname2(1:128),'(a)') localobsfile(1:il)
236	gforget	1.14	inquire( file=fname2, exist=exst )
237	jmc	1.15	#ifndef COST_GENERIC_ASSUME_CYCLIC
238	gforget	1.12	c assume we have one big file, one year after the other
239	heimbach	1.3	localrec = obsrec
240	jmc	1.15	c otherwise assume climatology, used for each year
241	gforget	1.12	#endif
242	heimbach	1.3	endif
243
244	gforget	1.14	if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) ) then
245	jmc	1.8	call mdsreadfield( fname2, cost_iprec, cost_yftype, nnzobs,
246	heimbach	1.3	& localobs, localrec, mythid )
247			else
248			do bj = jtlo,jthi
249			do bi = itlo,ithi
250			do k = 1,nnzobs
251			do j = jmin,jmax
252			do i = imin,imax
253			localobs(i,j,k,bi,bj) = spval
254			enddo
255			enddo
256			enddo
257			enddo
258			enddo
259	heimbach	1.1	endif
260			cnew)
261
262			do bj = jtlo,jthi
263			do bi = itlo,ithi
264
265			localcost = 0. _d 0
266
267			c-- Determine the mask on weights
268			do k = 1,nnzobs
269			do j = jmin,jmax
270			do i = imin,imax
271			cmask(i,j,k) = cosphi(i,j,bi,bj)*localmask(i,j,k,bi,bj)
272			if ( localobs(i,j,k,bi,bj) .lt. spminloc .or.
273			& localobs(i,j,k,bi,bj) .gt. spmaxloc .or.
274			& localobs(i,j,k,bi,bj) .eq. spzeroloc ) then
275			cmask(i,j,k) = 0. _d 0
276			endif
277			enddo
278			enddo
279			enddo
280			c--
281			do k = 1,nnzobs
282			do j = jmin,jmax
283			do i = imin,imax
284			localwww = localweight(i,j,k,bi,bj)*cmask(i,j,k)
285			junk = ( localbar(i,j,k,bi,bj) -
286			& localobs(i,j,k,bi,bj) )
287			localcost = localcost + junkjunklocalwww
288			if ( localwww .ne. 0. )
289			& num_local(bi,bj) = num_local(bi,bj) + 1. _d 0
290			enddo
291			enddo
292			enddo
293
294			objf_local(bi,bj) = objf_local(bi,bj) + localcost
295
296			enddo
297			enddo
298
299			enddo
300			c-- End of second loop over records.
301
302	heimbach	1.2	c-- End of mult_local or localobsfile
303			endif
304
305	heimbach	1.1	#endif /* ifdef ALLOW_COST */
306
307			end