pkg/ecco/cost_obcsw.F

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

#include "COST_CPPOPTIONS.h"

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

C     !DESCRIPTION: \bv
c     ==================================================================
c     SUBROUTINE cost_obcsw
c     ==================================================================
c
c     o cost function contribution obc
c
c     ==================================================================
c     SUBROUTINE cost_obcsw
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 iobcs
      integer ip1
      integer nrec
      integer ilfld
      integer igg

      _RL fctile
      _RL fcthread
      _RL dummy
      _RL gg
      _RL tmpx
      _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_OBCSW_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_obcsw: 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_obcsw_file )
        write(fnamefld(1:80),'(2a,i10.10)')
     &       xx_obcsw_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_obcsw_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( 'maskobcsw', 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 + wobcsw2(j,k,bi,bj,iobcs)
                  fctile = fctile + wobcsw(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 (wobcsw(k,iobcs)*maskyz(j,k,bi,bj).ne.0.)
     &                    num_obcsw(bi,bj) = num_obcsw(bi,bj) + 1. _d 0
              enddo
            enddo

            objf_obcsw(bi,bj) = objf_obcsw(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_obcsw: irec = ',irec
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
        write(msgbuf,'(a,a,d22.15)')
     &    ' global cost function value',
     &    ' (obcsw) = ',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.9	C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_obcsw.F,v 1.8 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_OBCSS
8			C !INTERFACE:
9	heimbach	1.1	subroutine cost_obcsw(
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_obcsw
20			c ==================================================================
21			c
22			c o cost function contribution obc
23			c
24			c ==================================================================
25			c SUBROUTINE cost_obcsw
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.7	c#include "DYNVARS.h"
40	heimbach	1.1	#ifdef ALLOW_OBCS
41	jmc	1.7	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.9	#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 iobcs
72			integer ip1
73	mlosch	1.5	integer nrec
74			integer ilfld
75			integer igg
76	heimbach	1.1
77			_RL fctile
78			_RL fcthread
79			_RL dummy
80	mlosch	1.5	_RL gg
81			_RL tmpx
82			_RL tmpfield (1-oly:sny+oly,nr,nsx,nsy)
83			_RL maskyz (1-oly:sny+oly,nr,nsx,nsy)
84	heimbach	1.1
85	mlosch	1.5	character*(80) fnamefld
86	heimbach	1.1
87			logical doglobalread
88			logical ladinit
89
90			#ifdef ECCO_VERBOSE
91			character*(MAX_LEN_MBUF) msgbuf
92			#endif
93
94			c == external functions ==
95
96			integer ilnblnk
97			external ilnblnk
98
99			c == end of interface ==
100	mlosch	1.5	CEOP
101	heimbach	1.1
102			jtlo = mybylo(mythid)
103			jthi = mybyhi(mythid)
104			itlo = mybxlo(mythid)
105			ithi = mybxhi(mythid)
106			jmin = 1
107			jmax = sny
108			imin = 1
109			imax = snx
110
111			c-- Read tiled data.
112			doglobalread = .false.
113			ladinit = .false.
114
115	mlosch	1.5	c Number of records to be used.
116			nrec = endrec-startrec+1
117
118	heimbach	1.1	#ifdef ALLOW_OBCSW_COST_CONTRIBUTION
119
120			ip1 = 1
121			fcthread = 0. _d 0
122
123	mlosch	1.5	#ifdef ECCO_VERBOSE
124			_BEGIN_MASTER( mythid )
125			write(msgbuf,'(a)') ' '
126			call print_message( msgbuf, standardmessageunit,
127			& SQUEEZE_RIGHT , mythid)
128			write(msgbuf,'(a)') ' '
129			call print_message( msgbuf, standardmessageunit,
130			& SQUEEZE_RIGHT , mythid)
131			write(msgbuf,'(a,i9.8)')
132			& ' cost_obcsw: number of records to process: ',nrec
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			_END_MASTER( mythid )
139			#endif
140
141			if (optimcycle .ge. 0) then
142			ilfld=ilnblnk( xx_obcsw_file )
143	jmc	1.7	write(fnamefld(1:80),'(2a,i10.10)')
144	mlosch	1.5	& xx_obcsw_file(1:ilfld), '.', optimcycle
145			endif
146
147	heimbach	1.2	c-- Loop over records.
148	mlosch	1.5	do irec = 1,nrec
149	heimbach	1.1
150	mlosch	1.5	call active_read_yz( fnamefld, tmpfield, irec, doglobalread,
151			& ladinit, optimcycle, mythid
152			& , xx_obcsw_dummy )
153
154			cgg Need to solve for iobcs would have been.
155			gg = (irec-1)/nobcs
156			igg = int(gg)
157			iobcs = irec - igg*nobcs
158
159	jmc	1.7	call active_read_yz( 'maskobcsw', maskyz,
160	mlosch	1.5	& iobcs,
161			& doglobalread, ladinit, 0,
162			& mythid, dummy )
163	heimbach	1.1
164	jmc	1.8	c-- Loop over this thread s tiles.
165	heimbach	1.1	do bj = jtlo,jthi
166			do bi = itlo,ithi
167	heimbach	1.2
168	mlosch	1.5	c-- Determine the weights to be used.
169	heimbach	1.2	fctile = 0. _d 0
170	heimbach	1.1
171	mlosch	1.5	do k = 1, Nr
172			do j = jmin,jmax
173			i = OB_Iw(j,bi,bj)
174			cgg if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
175			tmpx = tmpfield(j,k,bi,bj)
176			CMM fctile = fctile + wobcsw2(j,k,bi,bj,iobcs)
177			fctile = fctile + wobcsw(k,iobcs)
178			& tmpxtmpx*maskyz(j,k,bi,bj)
179			cgg endif
180			CMM if (wobcsw2(j,k,bi,bj,iobcs)*maskyz(j,k,bi,bj).ne.0.)
181			if (wobcsw(k,iobcs)*maskyz(j,k,bi,bj).ne.0.)
182			& num_obcsw(bi,bj) = num_obcsw(bi,bj) + 1. _d 0
183			enddo
184	heimbach	1.2	enddo
185
186	heimbach	1.1	objf_obcsw(bi,bj) = objf_obcsw(bi,bj) + fctile
187	mlosch	1.5	fcthread = fcthread + fctile
188	heimbach	1.1	enddo
189			enddo
190
191			#ifdef ECCO_VERBOSE
192			c-- Print cost function for all tiles.
193	jmc	1.4	_GLOBAL_SUM_RL( fcthread , myThid )
194	heimbach	1.1	write(msgbuf,'(a)') ' '
195			call print_message( msgbuf, standardmessageunit,
196			& SQUEEZE_RIGHT , mythid)
197			write(msgbuf,'(a,i8.8)')
198			& ' cost_obcsw: irec = ',irec
199			call print_message( msgbuf, standardmessageunit,
200			& SQUEEZE_RIGHT , mythid)
201			write(msgbuf,'(a,a,d22.15)')
202			& ' global cost function value',
203			& ' (obcsw) = ',fcthread
204			call print_message( msgbuf, standardmessageunit,
205			& SQUEEZE_RIGHT , mythid)
206			write(msgbuf,'(a)') ' '
207			call print_message( msgbuf, standardmessageunit,
208			& SQUEEZE_RIGHT , mythid)
209			#endif
210
211			enddo
212			c-- End of loop over records.
213
214			#endif
215
216			return
217			end