pkg/ecco/cost_driftw.F

C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_driftw.F,v 1.5 2007/10/09 00:02:50 jmc Exp $
C $Name:  $

#include "COST_CPPOPTIONS.h"


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

c     ==================================================================
c     SUBROUTINE cost_Driftw
c     ==================================================================
c
c     o Evaluate cost function contribution of the w difference
c       between the first and the last year.
c
c     started: from cost_drift
c
c     Armin Koehl  akoehl@ucsd.edu  26-Feb-2002
c
c     ==================================================================
c     SUBROUTINE cost_Driftw
c     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "GRID.h"
#include "DYNVARS.h"

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

c     == routine arguments ==

      integer myiter
      _RL     mytime
      integer mythid

c     == local variables ==

      _RS        one_rs
      parameter( one_rs = 1. )

      integer bi,bj
      integer i,j,k
      integer itlo,ithi
      integer jtlo,jthi
      integer jmin,jmax
      integer imin,imax
      integer irec
      integer ilw
      integer nf, nl
      integer minrec

      _RL fctilew
      _RL fcthread_wdrift

      character*(80) fnamew

      logical doglobalread
      logical ladinit

      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 tiled data.
      doglobalread = .false.
      ladinit      = .false.

#ifdef ALLOW_DRIFTW_COST_CONTRIBUTION

      if (optimcycle .ge. 0) then
        ilw = ilnblnk( wbarfile )
        write(fnamew(1:80),'(2a,i10.10)')
     &    wbarfile(1:ilw),'.',optimcycle
      endif

      fcthread_wdrift = 0. _d 0

      do bj = jtlo,jthi
        do bi = itlo,ithi
          do k = 1,nr
            do j = jmin,jmax
              do i = imin,imax
                wfmean(i,j,k,bi,bj) = 0.0
                wlmean(i,j,k,bi,bj) = 0.0
              enddo
            enddo
          enddo
        enddo
      enddo

      nf = 0
      nl = 0
c--   Number of full years
      nfmin = MAX(INT(FLOAT(nmonsrec)/12.),1)
c--   Prevent code from crashing if integrated for less than a year
      minrec = MIN(nmonsrec,12)

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

c--     Read time averages and the monthly mean data.
        call active_read_xyz( fnamew, wbar, irec,
     &                        doglobalread, ladinit,
     &                        optimcycle, mythid,
     &                        xx_wbar_mean_dummy )

        nf = nf + 1
        do bj = jtlo,jthi
          do bi = itlo,ithi
            do k = 1,nr
              do j = jmin,jmax
                do i = imin,imax
                  wfmean(i,j,k,bi,bj) = wfmean(i,j,k,bi,bj) +
     &                  wbar(i,j,k,bi,bj)
                enddo
              enddo
            enddo
          enddo
        enddo

      enddo

      do irec = nmonsrec-minrec+1, nmonsrec
 
c--     Read time averages and the monthly mean data.
        call active_read_xyz( fnamew, wbar, irec,
     &                        doglobalread, ladinit,
     &                        optimcycle, mythid,
     &                        xx_wbar_mean_dummy )

        nl = nl + 1

        do bj = jtlo,jthi
          do bi = itlo,ithi
            do k = 1,nr
              do j = jmin,jmax
                do i = imin,imax
                  wlmean(i,j,k,bi,bj) = wlmean(i,j,k,bi,bj) +
     &                  wbar(i,j,k,bi,bj)
                enddo
              enddo
            enddo
          enddo
        enddo

      enddo


      do bj = jtlo,jthi
        do bi = itlo,ithi

c--       Loop over the model layers
          fctilew = 0. _d 0

          do k = 1,nr

c--         Compute model misfit and cost function term for
c           the vertical velovity field. The error is 1e-4 m/s.
            do j = jmin,jmax
             do i = imin,imax
              if (_hFacC(i,j,k,bi,bj) .ne. 0.) then
                     fctilew = fctilew +
     &                  (2.5e11*cosphi(i,j,bi,bj)*
     &          (wlmean(i,j,k,bi,bj)/nl - wfmean(i,j,k,bi,bj)/nf)*
     &          (wlmean(i,j,k,bi,bj)/nl - wfmean(i,j,k,bi,bj)/nf))
                     if ( cosphi(i,j,bi,bj) .ne. 0. )
     &                 num_wdrift(bi,bj) = num_wdrift(bi,bj) + 1. _d 0
              endif
             enddo
            enddo

          enddo
c--       End of loop over layers.

          fcthread_wdrift    = fcthread_wdrift    + fctilew
          objf_wdrift(bi,bj) = objf_wdrift(bi,bj) + fctilew

#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_Driftw: irec,bi,bj          =  ',irec,bi,bj
          call print_message( msgbuf, standardmessageunit,
     &                          SQUEEZE_RIGHT , mythid)
          write(msgbuf,'(a,d22.15)')
     &        '     cost function (wvel) = ',
     &        fcthread_wdrift
          call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
          write(msgbuf,'(a)') ' '
          call print_message( msgbuf, standardmessageunit,
     &                          SQUEEZE_RIGHT , mythid)
#endif

        enddo
      enddo

#ifdef ECCO_VERBOSE
c--     Print cost function for all tiles.
        _GLOBAL_SUM_R8( fcthread_wdrift , myThid )
        write(msgbuf,'(a)') ' '
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
        write(msgbuf,'(a,i8.8)')
     &    ' cost_Driftw: irec = ',irec
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
        write(msgbuf,'(a,a,d22.15)')
     &    '  cost function value',
     &    ' (wvel) = ',fcthread_wdrift
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
        write(msgbuf,'(a)') ' '
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
#endif

#endif

      return
      end

1	heimbach	1.6	C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_driftw.F,v 1.5 2007/10/09 00:02:50 jmc Exp $
2	jmc	1.5	C $Name: $
3	heimbach	1.1
4			#include "COST_CPPOPTIONS.h"
5
6
7			subroutine cost_driftw(
8			I myiter,
9			I mytime,
10			I mythid
11			& )
12
13			c ==================================================================
14			c SUBROUTINE cost_Driftw
15			c ==================================================================
16			c
17			c o Evaluate cost function contribution of the w difference
18			c between the first and the last year.
19			c
20			c started: from cost_drift
21			c
22			c Armin Koehl akoehl@ucsd.edu 26-Feb-2002
23			c
24			c ==================================================================
25			c SUBROUTINE cost_Driftw
26			c ==================================================================
27
28			implicit none
29
30			c == global variables ==
31
32			#include "EEPARAMS.h"
33			#include "SIZE.h"
34			#include "GRID.h"
35			#include "DYNVARS.h"
36
37			#include "cal.h"
38			#include "ecco_cost.h"
39			#include "ctrl.h"
40			#include "ctrl_dummy.h"
41			#include "optim.h"
42
43			c == routine arguments ==
44
45			integer myiter
46			_RL mytime
47			integer mythid
48
49			c == local variables ==
50
51			_RS one_rs
52			parameter( one_rs = 1. )
53
54			integer bi,bj
55			integer i,j,k
56			integer itlo,ithi
57			integer jtlo,jthi
58			integer jmin,jmax
59			integer imin,imax
60			integer irec
61			integer ilw
62			integer nf, nl
63	heimbach	1.6	integer minrec
64	heimbach	1.1
65			_RL fctilew
66			_RL fcthread_wdrift
67
68			character*(80) fnamew
69
70			logical doglobalread
71			logical ladinit
72
73			character*(MAX_LEN_MBUF) msgbuf
74
75			c == external functions ==
76
77			integer ilnblnk
78			external ilnblnk
79
80			c == end of interface ==
81
82			jtlo = mybylo(mythid)
83			jthi = mybyhi(mythid)
84			itlo = mybxlo(mythid)
85			ithi = mybxhi(mythid)
86			jmin = 1
87			jmax = sny
88			imin = 1
89			imax = snx
90
91			c-- Read tiled data.
92			doglobalread = .false.
93			ladinit = .false.
94
95			#ifdef ALLOW_DRIFTW_COST_CONTRIBUTION
96
97			if (optimcycle .ge. 0) then
98			ilw = ilnblnk( wbarfile )
99			write(fnamew(1:80),'(2a,i10.10)')
100			& wbarfile(1:ilw),'.',optimcycle
101			endif
102	jmc	1.5
103	heimbach	1.1	fcthread_wdrift = 0. _d 0
104	jmc	1.5
105	heimbach	1.1	do bj = jtlo,jthi
106			do bi = itlo,ithi
107	jmc	1.5	do k = 1,nr
108	heimbach	1.1	do j = jmin,jmax
109			do i = imin,imax
110			wfmean(i,j,k,bi,bj) = 0.0
111			wlmean(i,j,k,bi,bj) = 0.0
112			enddo
113	jmc	1.5	enddo
114	heimbach	1.1	enddo
115			enddo
116			enddo
117	jmc	1.5
118	heimbach	1.1	nf = 0
119			nl = 0
120	heimbach	1.6	c-- Number of full years
121			nfmin = MAX(INT(FLOAT(nmonsrec)/12.),1)
122			c-- Prevent code from crashing if integrated for less than a year
123			minrec = MIN(nmonsrec,12)
124	heimbach	1.1
125			c-- Loop over records.
126	heimbach	1.6	do irec = 1,minrec
127	heimbach	1.1
128			c-- Read time averages and the monthly mean data.
129			call active_read_xyz( fnamew, wbar, irec,
130			& doglobalread, ladinit,
131			& optimcycle, mythid,
132			& xx_wbar_mean_dummy )
133	jmc	1.5
134	heimbach	1.1	nf = nf + 1
135			do bj = jtlo,jthi
136			do bi = itlo,ithi
137	jmc	1.5	do k = 1,nr
138	heimbach	1.1	do j = jmin,jmax
139			do i = imin,imax
140	jmc	1.5	wfmean(i,j,k,bi,bj) = wfmean(i,j,k,bi,bj) +
141	heimbach	1.1	& wbar(i,j,k,bi,bj)
142			enddo
143	jmc	1.5	enddo
144	heimbach	1.1	enddo
145			enddo
146			enddo
147	jmc	1.5
148	heimbach	1.1	enddo
149	jmc	1.5
150	heimbach	1.6	do irec = nmonsrec-minrec+1, nmonsrec
151
152	heimbach	1.1	c-- Read time averages and the monthly mean data.
153			call active_read_xyz( fnamew, wbar, irec,
154			& doglobalread, ladinit,
155			& optimcycle, mythid,
156			& xx_wbar_mean_dummy )
157	jmc	1.5
158			nl = nl + 1
159
160	heimbach	1.1	do bj = jtlo,jthi
161			do bi = itlo,ithi
162	jmc	1.5	do k = 1,nr
163	heimbach	1.1	do j = jmin,jmax
164			do i = imin,imax
165	jmc	1.5	wlmean(i,j,k,bi,bj) = wlmean(i,j,k,bi,bj) +
166	heimbach	1.1	& wbar(i,j,k,bi,bj)
167			enddo
168	jmc	1.5	enddo
169	heimbach	1.1	enddo
170			enddo
171			enddo
172	jmc	1.5
173	heimbach	1.1	enddo
174	jmc	1.5
175
176	heimbach	1.1	do bj = jtlo,jthi
177			do bi = itlo,ithi
178
179			c-- Loop over the model layers
180			fctilew = 0. _d 0
181	jmc	1.5
182	heimbach	1.1	do k = 1,nr
183
184			c-- Compute model misfit and cost function term for
185			c the vertical velovity field. The error is 1e-4 m/s.
186			do j = jmin,jmax
187	heimbach	1.3	do i = imin,imax
188			if (_hFacC(i,j,k,bi,bj) .ne. 0.) then
189	heimbach	1.1	fctilew = fctilew +
190			& (2.5e11cosphi(i,j,bi,bj)
191			& (wlmean(i,j,k,bi,bj)/nl - wfmean(i,j,k,bi,bj)/nf)*
192			& (wlmean(i,j,k,bi,bj)/nl - wfmean(i,j,k,bi,bj)/nf))
193	heimbach	1.3	if ( cosphi(i,j,bi,bj) .ne. 0. )
194			& num_wdrift(bi,bj) = num_wdrift(bi,bj) + 1. _d 0
195			endif
196			enddo
197	heimbach	1.1	enddo
198
199			enddo
200			c-- End of loop over layers.
201
202			fcthread_wdrift = fcthread_wdrift + fctilew
203			objf_wdrift(bi,bj) = objf_wdrift(bi,bj) + fctilew
204
205			#ifdef ECCO_VERBOSE
206			c-- Print cost function for each tile in each thread.
207			write(msgbuf,'(a)') ' '
208			call print_message( msgbuf, standardmessageunit,
209			& SQUEEZE_RIGHT , mythid)
210			write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
211			& ' cost_Driftw: irec,bi,bj = ',irec,bi,bj
212			call print_message( msgbuf, standardmessageunit,
213			& SQUEEZE_RIGHT , mythid)
214			write(msgbuf,'(a,d22.15)')
215			& ' cost function (wvel) = ',
216	heimbach	1.4	& fcthread_wdrift
217	heimbach	1.1	call print_message( msgbuf, standardmessageunit,
218			& SQUEEZE_RIGHT , mythid)
219			write(msgbuf,'(a)') ' '
220			call print_message( msgbuf, standardmessageunit,
221			& SQUEEZE_RIGHT , mythid)
222			#endif
223
224			enddo
225			enddo
226
227			#ifdef ECCO_VERBOSE
228			c-- Print cost function for all tiles.
229	heimbach	1.4	_GLOBAL_SUM_R8( fcthread_wdrift , myThid )
230	heimbach	1.1	write(msgbuf,'(a)') ' '
231			call print_message( msgbuf, standardmessageunit,
232			& SQUEEZE_RIGHT , mythid)
233			write(msgbuf,'(a,i8.8)')
234	heimbach	1.4	& ' cost_Driftw: irec = ',irec
235	heimbach	1.1	call print_message( msgbuf, standardmessageunit,
236			& SQUEEZE_RIGHT , mythid)
237			write(msgbuf,'(a,a,d22.15)')
238			& ' cost function value',
239	heimbach	1.4	& ' (wvel) = ',fcthread_wdrift
240	heimbach	1.1	call print_message( msgbuf, standardmessageunit,
241			& SQUEEZE_RIGHT , mythid)
242			write(msgbuf,'(a)') ' '
243			call print_message( msgbuf, standardmessageunit,
244			& SQUEEZE_RIGHT , mythid)
245			#endif
246
247			#endif
248
249			return
250			end
251