pkg/ecco/cost_obcse.F

C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_obcse.F,v 1.7 2011/06/24 19:39:32 jmc Exp $
C $Name:  $

#include "COST_CPPOPTIONS.h"

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#include "DYNVARS.h"
#ifdef ALLOW_OBCS
c# include "OBCS_PARAMS.h"
# include "OBCS_GRID.h"
#endif

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