pkg/ecco/cost_kapgm.F

C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_kapgm.F,v 1.11 2014/10/18 18:15:44 gforget Exp $
C $Name:  $

#include "ECCO_OPTIONS.h"


      subroutine cost_kapgm(
     I                           myiter,
     I                           mytime,
     I                           mythid
     &                         )

C     o==========================================================o
C     | subroutine cost_kapgm                                    |
C     | o GM coefficient adjustment penalization                 |
C     o==========================================================o

      implicit none

c     == global variables ==

#ifdef ALLOW_KAPGM_COST_CONTRIBUTION
#include "EEPARAMS.h"
#include "SIZE.h"
#include "GRID.h"
#include "DYNVARS.h"
#ifdef ALLOW_GMREDI
# include "GMREDI.h"
#endif

#include "ecco_cost.h"
#include "CTRL_SIZE.h"
#include "ctrl.h"
#include "ctrl_dummy.h"
#include "optim.h"
#endif

c     == routine arguments ==

      integer myiter
      _RL     mytime
      integer mythid

c     == local variables ==

#ifdef ALLOW_KAPGM_COST_CONTRIBUTION

      integer bi,bj
      integer i,j,k
      integer itlo,ithi
      integer jtlo,jthi
      integer jmin,jmax
      integer imin,imax
      integer nrec
      integer irec
      integer ilfld

      _RL fctile
      _RL fcthread
      _RL tmpx

      logical doglobalread
      logical ladinit

      character*(80) fnamefld

      character*(MAX_LEN_MBUF) msgbuf

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--   Read state record from global file.
      doglobalread = .false.
      ladinit      = .false.

      irec = 1

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

      fcthread = 0. _d 0

      call active_read_xyz( fnamefld, tmpfld3d, irec, doglobalread,
     &                       ladinit, optimcycle, mythid
     &        , xx_kapgm_dummy )

c--     Loop over this thread 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
              do i = imin,imax
                if (_hFacC(i,j,k,bi,bj) .ne. 0.) then
c                  tmpx = (tmpfld3d(i,j,k,bi,bj)-GM_background_K)
                  tmpx = tmpfld3d(i,j,k,bi,bj)
                IF ( .NOT.ctrlSmoothCorrel3D ) THEN
                  fctile = fctile
     &                 + wkapgmFld(i,j,k,bi,bj)*cosphi(i,j,bi,bj)
     &                 *tmpx*tmpx
                ELSE !IF ( .NOT.ctrlSmoothCorrel3D ) THEN
                  fctile = fctile + tmpx*tmpx
                ENDIF !IF ( .NOT.ctrlSmoothCorrel3D ) THEN
                endif
              enddo
            enddo
            enddo

            objf_kapgm(bi,bj) = objf_kapgm(bi,bj) + fctile
            fcthread          = fcthread + fctile

#ifdef ECCO_VERBOSE
c--         Print cost function for each tile in each thread.
            write(msgbuf,'(a)') ' '
            call print_message( msgbuf, standardmessageunit,
     &                          SQUEEZE_RIGHT , mythid)
            write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
     &        ' cost_kapgm: irec,bi,bj          =  ',irec,bi,bj
            call print_message( msgbuf, standardmessageunit,
     &                          SQUEEZE_RIGHT , mythid)
            write(msgbuf,'(a,d22.15)')
     &        '               cost function (dT(0)) = ',
     &        fctile
            call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
#endif
          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_kapgm:                       irec =  ',irec
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
        write(msgbuf,'(a,d22.15)')
     &    '                 global cost function value = ',
     &    fcthread
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
        write(msgbuf,'(a)') ' '
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
#endif

#endif

      return
      end


1	gforget	1.12	C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_kapgm.F,v 1.11 2014/10/18 18:15:44 gforget Exp $
2	jmc	1.5	C $Name: $
3	heimbach	1.1
4	jmc	1.9	#include "ECCO_OPTIONS.h"
5	heimbach	1.1
6
7			subroutine cost_kapgm(
8			I myiter,
9			I mytime,
10			I mythid
11			& )
12
13	gforget	1.4	C o==========================================================o
14			C \| subroutine cost_kapgm \|
15			C \| o GM coefficient adjustment penalization \|
16			C o==========================================================o
17	heimbach	1.1
18			implicit none
19
20			c == global variables ==
21
22	gforget	1.11	#ifdef ALLOW_KAPGM_COST_CONTRIBUTION
23	heimbach	1.1	#include "EEPARAMS.h"
24			#include "SIZE.h"
25			#include "GRID.h"
26			#include "DYNVARS.h"
27	heimbach	1.2	#ifdef ALLOW_GMREDI
28			# include "GMREDI.h"
29			#endif
30	heimbach	1.1
31	gforget	1.4	#include "ecco_cost.h"
32	heimbach	1.8	#include "CTRL_SIZE.h"
33	heimbach	1.1	#include "ctrl.h"
34			#include "ctrl_dummy.h"
35			#include "optim.h"
36	gforget	1.11	#endif
37
38	heimbach	1.1	c == routine arguments ==
39
40			integer myiter
41			_RL mytime
42			integer mythid
43
44			c == local variables ==
45
46	gforget	1.11	#ifdef ALLOW_KAPGM_COST_CONTRIBUTION
47
48	heimbach	1.1	integer bi,bj
49			integer i,j,k
50			integer itlo,ithi
51			integer jtlo,jthi
52			integer jmin,jmax
53			integer imin,imax
54			integer nrec
55			integer irec
56			integer ilfld
57
58			_RL fctile
59			_RL fcthread
60			_RL tmpx
61
62			logical doglobalread
63			logical ladinit
64
65			character*(80) fnamefld
66
67			character*(MAX_LEN_MBUF) msgbuf
68
69			c == external functions ==
70
71			integer ilnblnk
72			external ilnblnk
73
74			c == end of interface ==
75
76			jtlo = mybylo(mythid)
77			jthi = mybyhi(mythid)
78			itlo = mybxlo(mythid)
79			ithi = mybxhi(mythid)
80			jmin = 1
81			jmax = sny
82			imin = 1
83			imax = snx
84
85			c-- Read state record from global file.
86			doglobalread = .false.
87			ladinit = .false.
88	jmc	1.5
89	heimbach	1.1	irec = 1
90
91			if (optimcycle .ge. 0) then
92			ilfld = ilnblnk( xx_kapgm_file )
93	jmc	1.5	write(fnamefld(1:80),'(2a,i10.10)')
94	heimbach	1.1	& xx_kapgm_file(1:ilfld),'.',optimcycle
95			endif
96
97			fcthread = 0. _d 0
98
99	heimbach	1.3	call active_read_xyz( fnamefld, tmpfld3d, irec, doglobalread,
100	heimbach	1.1	& ladinit, optimcycle, mythid
101			& , xx_kapgm_dummy )
102
103	jmc	1.7	c-- Loop over this thread tiles.
104	heimbach	1.1	do bj = jtlo,jthi
105			do bi = itlo,ithi
106
107			c-- Determine the weights to be used.
108	jmc	1.5
109	heimbach	1.1	fctile = 0. _d 0
110			do k = 1,nr
111			do j = jmin,jmax
112			do i = imin,imax
113			if (_hFacC(i,j,k,bi,bj) .ne. 0.) then
114	jmc	1.5	c tmpx = (tmpfld3d(i,j,k,bi,bj)-GM_background_K)
115	gforget	1.4	tmpx = tmpfld3d(i,j,k,bi,bj)
116	gforget	1.10	IF ( .NOT.ctrlSmoothCorrel3D ) THEN
117	heimbach	1.1	fctile = fctile
118	gforget	1.4	& + wkapgmFld(i,j,k,bi,bj)*cosphi(i,j,bi,bj)
119	heimbach	1.1	& tmpxtmpx
120	gforget	1.10	ELSE !IF ( .NOT.ctrlSmoothCorrel3D ) THEN
121	gforget	1.4	fctile = fctile + tmpx*tmpx
122	gforget	1.10	ENDIF !IF ( .NOT.ctrlSmoothCorrel3D ) THEN
123	heimbach	1.1	endif
124			enddo
125			enddo
126			enddo
127
128			objf_kapgm(bi,bj) = objf_kapgm(bi,bj) + fctile
129			fcthread = fcthread + fctile
130
131			#ifdef ECCO_VERBOSE
132			c-- Print cost function for each tile in each thread.
133			write(msgbuf,'(a)') ' '
134			call print_message( msgbuf, standardmessageunit,
135			& SQUEEZE_RIGHT , mythid)
136			write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
137			& ' cost_kapgm: irec,bi,bj = ',irec,bi,bj
138			call print_message( msgbuf, standardmessageunit,
139			& SQUEEZE_RIGHT , mythid)
140			write(msgbuf,'(a,d22.15)')
141			& ' cost function (dT(0)) = ',
142			& fctile
143			call print_message( msgbuf, standardmessageunit,
144			& SQUEEZE_RIGHT , mythid)
145			#endif
146			enddo
147			enddo
148
149			#ifdef ECCO_VERBOSE
150			c-- Print cost function for all tiles.
151	jmc	1.6	_GLOBAL_SUM_RL( fcthread , myThid )
152	heimbach	1.1	write(msgbuf,'(a)') ' '
153			call print_message( msgbuf, standardmessageunit,
154			& SQUEEZE_RIGHT , mythid)
155			write(msgbuf,'(a,i8.8)')
156	gforget	1.4	& ' cost_kapgm: irec = ',irec
157	heimbach	1.1	call print_message( msgbuf, standardmessageunit,
158			& SQUEEZE_RIGHT , mythid)
159			write(msgbuf,'(a,d22.15)')
160			& ' global cost function value = ',
161			& fcthread
162			call print_message( msgbuf, standardmessageunit,
163			& SQUEEZE_RIGHT , mythid)
164			write(msgbuf,'(a)') ' '
165			call print_message( msgbuf, standardmessageunit,
166			& SQUEEZE_RIGHT , mythid)
167			#endif
168
169			#endif
170
171			return
172			end
173	jmc	1.5
174