pkg/ecco/cost_obcse.F

C $Header: $
C $Name:  $

#include "COST_CPPOPTIONS.h"
#ifdef ALLOW_OBCS
# include "OBCS_OPTIONS.h"
#endif

CBOP
C     !ROUTINE: COST_OBCSE
C     !INTERFACE:
      subroutine cost_obcse(
     I                       myiter,
     I                       mytime,
     I                       startrec,
     I                       endrec,
     I                       mythid
     &                     )

C     !DESCRIPTION: \bv
c     ==================================================================
c     SUBROUTINE cost_obcse
c     ==================================================================
c
c     o cost function contribution obc
c
c     ==================================================================
c     SUBROUTINE cost_obcse
c     ==================================================================
C     \ev

C     !USES:

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#ifdef ALLOW_OBCS
# include "OBCS.h"
#endif

#include "cal.h"
#include "ecco_cost.h"
#include "ctrl.h"
#include "ctrl_dummy.h"
#include "optim.h"

C     !INPUT/OUTPUT PARAMETERS:
c     == routine arguments ==

      integer myiter
      _RL     mytime
      integer mythid
      integer startrec
      integer endrec

C     !LOCAL VARIABLES:
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 iobcs
      integer ip1
      integer nrec
      integer ilfld
      integer igg

      _RL fctile
      _RL fcthread
      _RL dummy
      _RL gg
      _RL tmpx
cgg(
      _RL tmpfield (1-oly:sny+oly,nr,nsx,nsy)
      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)

      character*(80) fnamefld

      logical doglobalread
      logical ladinit

#ifdef ECCO_VERBOSE
      character*(MAX_LEN_MBUF) msgbuf
#endif

c     == external functions ==

      integer  ilnblnk
      external ilnblnk

c     == end of interface ==
CEOP

      jtlo = mybylo(mythid)
      jthi = mybyhi(mythid)
      itlo = mybxlo(mythid)
      ithi = mybxhi(mythid)
      jmin = 1
      jmax = sny
      imin = 1
      imax = snx

c--   Read tiled data.
      doglobalread = .false.
      ladinit      = .false.

c     Number of records to be used.
      nrec = endrec-startrec+1

#ifdef ALLOW_OBCSE_COST_CONTRIBUTION

      ip1 = 1
      fcthread = 0. _d 0

#ifdef ECCO_VERBOSE
      _BEGIN_MASTER( mythid )
      write(msgbuf,'(a)') ' '
      call print_message( msgbuf, standardmessageunit,
     &                    SQUEEZE_RIGHT , mythid)
      write(msgbuf,'(a)') ' '
      call print_message( msgbuf, standardmessageunit,
     &                    SQUEEZE_RIGHT , mythid)
      write(msgbuf,'(a,i9.8)')
     &  ' cost_obcse: number of records to process: ',nrec
      call print_message( msgbuf, standardmessageunit,
     &                    SQUEEZE_RIGHT , mythid)
      write(msgbuf,'(a)') ' '
      call print_message( msgbuf, standardmessageunit,
     &                    SQUEEZE_RIGHT , mythid)
      _END_MASTER( mythid )
#endif

      if (optimcycle .ge. 0) then
        ilfld=ilnblnk( xx_obcse_file )
        write(fnamefld(1:80),'(2a,i10.10)') 
     &       xx_obcse_file(1:ilfld), '.', optimcycle
      endif

c--   Loop over records.
      do irec = 1,nrec

        call active_read_yz( fnamefld, tmpfield, irec, doglobalread,
     &                       ladinit, optimcycle, mythid
     &        , xx_obcse_dummy )

cgg    Need to solve for iobcs would have been.
          gg    = (irec-1)/nobcs
          igg   = int(gg)
          iobcs = irec - igg*nobcs

          call active_read_yz( 'maskobcse', maskyz, 
     &                         iobcs,
     &                         doglobalread, ladinit, 0,
     &                         mythid, dummy )

c--     Loop over this thread's tiles.
        do bj = jtlo,jthi
          do bi = itlo,ithi

c--         Determine the weights to be used.
            fctile = 0. _d 0

            do k = 1, Nr
              do j = jmin,jmax
                i = OB_Iw(j,bi,bj)
cgg                if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
                  tmpx = tmpfield(j,k,bi,bj)
CMM                  fctile = fctile + wobcse2(j,k,bi,bj,iobcs)
                  fctile = fctile + wobcse(k,iobcs)
     &                        *tmpx*tmpx*maskyz(j,k,bi,bj)
cgg                endif
CMM                  if (wobcsw2(j,k,bi,bj,iobcs)*maskyz(j,k,bi,bj).ne.0.)
                  if (wobcse(k,iobcs)*maskyz(j,k,bi,bj).ne.0.)
     &                    num_obcse(bi,bj) = num_obcse(bi,bj) + 1. _d 0
              enddo
            enddo

            objf_obcse(bi,bj) = objf_obcse(bi,bj) + fctile
            fcthread         = fcthread + fctile
          enddo
        enddo

#ifdef ECCO_VERBOSE
c--     Print cost function for all tiles.
        _GLOBAL_SUM_RL( fcthread , myThid )
        write(msgbuf,'(a)') ' '
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
        write(msgbuf,'(a,i8.8)')
     &    ' cost_obcse: irec = ',irec
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
        write(msgbuf,'(a,a,d22.15)')
     &    ' global cost function value',
     &    ' (obcse) = ',fcthread
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
        write(msgbuf,'(a)') ' '
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
#endif

      enddo
c--   End of loop over records.

#endif

      return
      end


1	mlosch	1.5	C $Header: $
2	jmc	1.3	C $Name: $
3	heimbach	1.1
4			#include "COST_CPPOPTIONS.h"
5			#ifdef ALLOW_OBCS
6			# include "OBCS_OPTIONS.h"
7			#endif
8
9	mlosch	1.5	CBOP
10			C !ROUTINE: COST_OBCSE
11			C !INTERFACE:
12	heimbach	1.1	subroutine cost_obcse(
13			I myiter,
14			I mytime,
15	mlosch	1.5	I startrec,
16			I endrec,
17	heimbach	1.1	I mythid
18			& )
19
20	mlosch	1.5	C !DESCRIPTION: \bv
21	heimbach	1.1	c ==================================================================
22			c SUBROUTINE cost_obcse
23			c ==================================================================
24			c
25			c o cost function contribution obc
26			c
27			c ==================================================================
28			c SUBROUTINE cost_obcse
29			c ==================================================================
30	mlosch	1.5	C \ev
31
32			C !USES:
33	heimbach	1.1
34			implicit none
35
36			c == global variables ==
37
38			#include "EEPARAMS.h"
39			#include "SIZE.h"
40			#include "PARAMS.h"
41			#include "GRID.h"
42			#include "DYNVARS.h"
43			#ifdef ALLOW_OBCS
44			# include "OBCS.h"
45			#endif
46
47			#include "cal.h"
48			#include "ecco_cost.h"
49			#include "ctrl.h"
50			#include "ctrl_dummy.h"
51			#include "optim.h"
52
53	mlosch	1.5	C !INPUT/OUTPUT PARAMETERS:
54	heimbach	1.1	c == routine arguments ==
55
56			integer myiter
57			_RL mytime
58			integer mythid
59	mlosch	1.5	integer startrec
60			integer endrec
61	heimbach	1.1
62	mlosch	1.5	C !LOCAL VARIABLES:
63	heimbach	1.1	c == local variables ==
64
65			integer bi,bj
66			integer i,j,k
67			integer itlo,ithi
68			integer jtlo,jthi
69			integer jmin,jmax
70			integer imin,imax
71			integer irec
72			integer il
73			integer iobcs
74			integer ip1
75	mlosch	1.5	integer nrec
76			integer ilfld
77			integer igg
78	heimbach	1.1
79			_RL fctile
80			_RL fcthread
81			_RL dummy
82	mlosch	1.5	_RL gg
83			_RL tmpx
84			cgg(
85			_RL tmpfield (1-oly:sny+oly,nr,nsx,nsy)
86			_RL maskyz (1-oly:sny+oly,nr,nsx,nsy)
87	heimbach	1.1
88	mlosch	1.5	character*(80) fnamefld
89	heimbach	1.1
90			logical doglobalread
91			logical ladinit
92
93			#ifdef ECCO_VERBOSE
94			character*(MAX_LEN_MBUF) msgbuf
95			#endif
96
97			c == external functions ==
98
99			integer ilnblnk
100			external ilnblnk
101
102			c == end of interface ==
103	mlosch	1.5	CEOP
104	heimbach	1.1
105			jtlo = mybylo(mythid)
106			jthi = mybyhi(mythid)
107			itlo = mybxlo(mythid)
108			ithi = mybxhi(mythid)
109			jmin = 1
110			jmax = sny
111			imin = 1
112			imax = snx
113
114			c-- Read tiled data.
115			doglobalread = .false.
116			ladinit = .false.
117
118	mlosch	1.5	c Number of records to be used.
119			nrec = endrec-startrec+1
120
121	heimbach	1.1	#ifdef ALLOW_OBCSE_COST_CONTRIBUTION
122
123	mlosch	1.5	ip1 = 1
124	heimbach	1.1	fcthread = 0. _d 0
125
126	mlosch	1.5	#ifdef ECCO_VERBOSE
127			_BEGIN_MASTER( mythid )
128			write(msgbuf,'(a)') ' '
129			call print_message( msgbuf, standardmessageunit,
130			& SQUEEZE_RIGHT , mythid)
131			write(msgbuf,'(a)') ' '
132			call print_message( msgbuf, standardmessageunit,
133			& SQUEEZE_RIGHT , mythid)
134			write(msgbuf,'(a,i9.8)')
135			& ' cost_obcse: number of records to process: ',nrec
136			call print_message( msgbuf, standardmessageunit,
137			& SQUEEZE_RIGHT , mythid)
138			write(msgbuf,'(a)') ' '
139			call print_message( msgbuf, standardmessageunit,
140			& SQUEEZE_RIGHT , mythid)
141			_END_MASTER( mythid )
142			#endif
143
144			if (optimcycle .ge. 0) then
145			ilfld=ilnblnk( xx_obcse_file )
146			write(fnamefld(1:80),'(2a,i10.10)')
147			& xx_obcse_file(1:ilfld), '.', optimcycle
148			endif
149
150	heimbach	1.2	c-- Loop over records.
151	mlosch	1.5	do irec = 1,nrec
152	heimbach	1.1
153	mlosch	1.5	call active_read_yz( fnamefld, tmpfield, irec, doglobalread,
154			& ladinit, optimcycle, mythid
155			& , xx_obcse_dummy )
156
157			cgg Need to solve for iobcs would have been.
158			gg = (irec-1)/nobcs
159			igg = int(gg)
160			iobcs = irec - igg*nobcs
161
162			call active_read_yz( 'maskobcse', maskyz,
163			& iobcs,
164			& doglobalread, ladinit, 0,
165			& mythid, dummy )
166	heimbach	1.1
167	mlosch	1.5	c-- Loop over this thread's tiles.
168	heimbach	1.1	do bj = jtlo,jthi
169			do bi = itlo,ithi
170	heimbach	1.2
171	mlosch	1.5	c-- Determine the weights to be used.
172	heimbach	1.2	fctile = 0. _d 0
173	heimbach	1.1
174	mlosch	1.5	do k = 1, Nr
175			do j = jmin,jmax
176			i = OB_Iw(j,bi,bj)
177			cgg if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
178			tmpx = tmpfield(j,k,bi,bj)
179			CMM fctile = fctile + wobcse2(j,k,bi,bj,iobcs)
180			fctile = fctile + wobcse(k,iobcs)
181			& tmpxtmpx*maskyz(j,k,bi,bj)
182			cgg endif
183			CMM if (wobcsw2(j,k,bi,bj,iobcs)*maskyz(j,k,bi,bj).ne.0.)
184			if (wobcse(k,iobcs)*maskyz(j,k,bi,bj).ne.0.)
185			& num_obcse(bi,bj) = num_obcse(bi,bj) + 1. _d 0
186			enddo
187	heimbach	1.1	enddo
188	heimbach	1.2
189	heimbach	1.1	objf_obcse(bi,bj) = objf_obcse(bi,bj) + fctile
190	mlosch	1.5	fcthread = fcthread + fctile
191	heimbach	1.1	enddo
192			enddo
193
194			#ifdef ECCO_VERBOSE
195			c-- Print cost function for all tiles.
196	jmc	1.4	_GLOBAL_SUM_RL( fcthread , myThid )
197	heimbach	1.1	write(msgbuf,'(a)') ' '
198			call print_message( msgbuf, standardmessageunit,
199			& SQUEEZE_RIGHT , mythid)
200			write(msgbuf,'(a,i8.8)')
201			& ' cost_obcse: irec = ',irec
202			call print_message( msgbuf, standardmessageunit,
203			& SQUEEZE_RIGHT , mythid)
204			write(msgbuf,'(a,a,d22.15)')
205			& ' global cost function value',
206			& ' (obcse) = ',fcthread
207			call print_message( msgbuf, standardmessageunit,
208			& SQUEEZE_RIGHT , mythid)
209			write(msgbuf,'(a)') ' '
210			call print_message( msgbuf, standardmessageunit,
211			& SQUEEZE_RIGHT , mythid)
212			#endif
213
214			enddo
215			c-- End of loop over records.
216
217			#endif
218
219			return
220			end
221
222
223
224
225
226
227