pkg/ecco/cost_obcse.F

C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_obcse.F,v 1.11 2014/10/09 00:50:16 gforget Exp $
C $Name:  $

#include "ECCO_OPTIONS.h"
#ifdef ALLOW_CTRL
# include "CTRL_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"
c#ifdef ALLOW_OBCS
c# include "OBCS_GRID.h"
c#endif

#ifdef ALLOW_CAL
# include "cal.h"
#endif
#ifdef ALLOW_CTRL
# include "CTRL_SIZE.h"
# include "ctrl.h"
# include "ctrl_dummy.h"
# include "optim.h"
# include "CTRL_OBCS.h"
#endif

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

      integer myiter
      _RL     mytime
      integer mythid
      integer startrec
      integer endrec

#if (defined (ALLOW_CTRL) && defined (ALLOW_OBCS))

#ifdef ALLOW_OBCSE_COST_CONTRIBUTION

c     == external functions ==
      integer  ilnblnk
      external ilnblnk

C     !LOCAL VARIABLES:
c     == local variables ==

      integer bi,bj
      integer j,k
      integer itlo,ithi
      integer jtlo,jthi
      integer jmin,jmax
      integer imin,imax
      integer irec
      integer il
      integer iobcs
c     integer i, 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     == 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

c     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
c               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 /* ALLOW_OBCSE_COST_CONTRIBUTION */

#endif /* ALLOW_CTRL and ALLOW_OBCS */

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