pkg/ecco/cost_theta0.F

C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_theta0.F,v 1.7 2007/05/14 22:06:55 heimbach Exp $

#include "COST_CPPOPTIONS.h"


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

c     ==================================================================
c     SUBROUTINE cost_zonstress
c     ==================================================================
c
c     o Calculate the zonal wind stress contribution to the cost function.
c
c     started: Christian Eckert eckert@mit.edu 30-Jun-1999
c
c     changed: Christian Eckert eckert@mit.edu 25-Feb-2000
c
c              - Restructured the code in order to create a package
c                for the MITgcmUV.
c
c     ==================================================================
c     SUBROUTINE cost_zonstress
c     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "GRID.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 ==

      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 tmpx
      _RL lengthscale

      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

      lengthscale = 1. _d 0

c--   Read state record from global file.
      doglobalread = .false.
      ladinit      = .false.
      
      irec = 1

#ifdef ALLOW_THETA0_COST_CONTRIBUTION

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

      call active_read_xyz( 
     &     fnamefld, tmpfld3d, irec, doglobalread,
     &     ladinit, optimcycle, mythid, xx_theta_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
              do i = imin,imax
                if (_hFacC(i,j,k,bi,bj) .ne. 0.) then
                  tmpx = tmpfld3d(i,j,k,bi,bj)
#ifndef ALLOW_SMOOTH_CORREL3D
                  if ( ABS(R_low(i,j,bi,bj)) .LT. 100. )
     &              tmpx = tmpx*ABS(R_low(i,j,bi,bj))/100.
                  fctile = fctile
     &                 + wthetaLev(i,j,k,bi,bj)*cosphi(i,j,bi,bj)
     &                 *tmpx*tmpx
#else
                  fctile = fctile + tmpx*tmpx
#endif
                  if ( wthetaLev(i,j,k,bi,bj)*cosphi(i,j,bi,bj).ne.0. )
     &                 num_temp0(bi,bj) = num_temp0(bi,bj) + 1. _d 0
                endif
              enddo
             enddo
            enddo

            objf_temp0(bi,bj) = objf_temp0(bi,bj) + fctile

          enddo
        enddo

#ifndef ALLOW_SMOOTH_CORREL3D
#ifdef ALLOW_SMOOTH_IC_COST_CONTRIBUTION

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

        _EXCH_XYZ_R8(tmpfld3d, mythid)

c--     Loop over this thread's tiles.
        do bj = jtlo,jthi
          do bi = itlo,ithi

            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
                  tmpx = 
     &              ( tmpfld3d(i+2,j,k,bi,bj)-tmpfld3d(i+1,j,k,bi,bj) )
     &                *maskW(i+1,j,k,bi,bj)*maskW(i+2,j,k,bi,bj)
     &            + ( tmpfld3d(i+1,j,k,bi,bj)-tmpfld3d(i,j,k,bi,bj) )
     &                *maskW(i+1,j,k,bi,bj)
     &            + ( tmpfld3d(i,j+2,k,bi,bj)-tmpfld3d(i,j+1,k,bi,bj) )
     &                *maskS(i,j+1,k,bi,bj)*maskS(i,j+2,k,bi,bj)
     &            + ( tmpfld3d(i,j+1,k,bi,bj)-tmpfld3d(i,j,k,bi,bj) )
     &                *maskS(i,j+1,k,bi,bj)
                  if ( ABS(R_low(i,j,bi,bj)) .LT. 100. )
     &              tmpx = tmpx*ABS(R_low(i,j,bi,bj))/100.
                  fctile = fctile
     &               + wthetaLev(i,j,k,bi,bj)*cosphi(i,j,bi,bj)
     *                 *0.0161*lengthscale/4.0
     &                 *tmpx*tmpx
                endif
              enddo
             enddo
            enddo

            objf_temp0smoo(bi,bj) = objf_temp0smoo(bi,bj) + fctile

          enddo
        enddo
#endif
#endif

#endif

      return
      end
      
      
1	gforget	1.8	C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_theta0.F,v 1.7 2007/05/14 22:06:55 heimbach Exp $
2	heimbach	1.1
3			#include "COST_CPPOPTIONS.h"
4
5
6			subroutine cost_theta0(
7			I myiter,
8			I mytime,
9			I mythid
10			& )
11
12			c ==================================================================
13			c SUBROUTINE cost_zonstress
14			c ==================================================================
15			c
16			c o Calculate the zonal wind stress contribution to the cost function.
17			c
18			c started: Christian Eckert eckert@mit.edu 30-Jun-1999
19			c
20			c changed: Christian Eckert eckert@mit.edu 25-Feb-2000
21			c
22			c - Restructured the code in order to create a package
23			c for the MITgcmUV.
24			c
25			c ==================================================================
26			c SUBROUTINE cost_zonstress
27			c ==================================================================
28
29			implicit none
30
31			c == global variables ==
32
33			#include "EEPARAMS.h"
34			#include "SIZE.h"
35			#include "GRID.h"
36
37			#include "ecco_cost.h"
38			#include "ctrl.h"
39			#include "ctrl_dummy.h"
40			#include "optim.h"
41
42			c == routine arguments ==
43
44			integer myiter
45			_RL mytime
46			integer mythid
47
48			c == local variables ==
49
50			integer bi,bj
51			integer i,j,k
52			integer itlo,ithi
53			integer jtlo,jthi
54			integer jmin,jmax
55			integer imin,imax
56			integer nrec
57			integer irec
58			integer ilfld
59
60			_RL fctile
61			_RL tmpx
62	heimbach	1.5	_RL lengthscale
63	heimbach	1.1
64			logical doglobalread
65			logical ladinit
66
67			character*(80) fnamefld
68
69			character*(MAX_LEN_MBUF) msgbuf
70
71			c == external functions ==
72
73			integer ilnblnk
74			external ilnblnk
75
76			c == end of interface ==
77
78			jtlo = mybylo(mythid)
79			jthi = mybyhi(mythid)
80			itlo = mybxlo(mythid)
81			ithi = mybxhi(mythid)
82			jmin = 1
83			jmax = sny
84			imin = 1
85			imax = snx
86
87	heimbach	1.5	lengthscale = 1. _d 0
88
89	heimbach	1.1	c-- Read state record from global file.
90			doglobalread = .false.
91			ladinit = .false.
92
93			irec = 1
94
95			#ifdef ALLOW_THETA0_COST_CONTRIBUTION
96
97			if (optimcycle .ge. 0) then
98			ilfld = ilnblnk( xx_theta_file )
99			write(fnamefld(1:80),'(2a,i10.10)')
100			& xx_theta_file(1:ilfld),'.',optimcycle
101			endif
102
103	heimbach	1.7	call active_read_xyz(
104	heimbach	1.5	& fnamefld, tmpfld3d, irec, doglobalread,
105			& ladinit, optimcycle, mythid, xx_theta_dummy )
106	heimbach	1.1
107			c-- Loop over this thread's tiles.
108			do bj = jtlo,jthi
109			do bi = itlo,ithi
110
111			c-- Determine the weights to be used.
112
113			fctile = 0. _d 0
114			do k = 1,nr
115	heimbach	1.6	do j = jmin,jmax
116	heimbach	1.1	do i = imin,imax
117			if (_hFacC(i,j,k,bi,bj) .ne. 0.) then
118			tmpx = tmpfld3d(i,j,k,bi,bj)
119	gforget	1.8	#ifndef ALLOW_SMOOTH_CORREL3D
120			if ( ABS(R_low(i,j,bi,bj)) .LT. 100. )
121	heimbach	1.4	& tmpx = tmpx*ABS(R_low(i,j,bi,bj))/100.
122	heimbach	1.1	fctile = fctile
123	heimbach	1.4	& + wthetaLev(i,j,k,bi,bj)*cosphi(i,j,bi,bj)
124	heimbach	1.1	& tmpxtmpx
125	gforget	1.8	#else
126			fctile = fctile + tmpx*tmpx
127			#endif
128	heimbach	1.4	if ( wthetaLev(i,j,k,bi,bj)*cosphi(i,j,bi,bj).ne.0. )
129	heimbach	1.3	& num_temp0(bi,bj) = num_temp0(bi,bj) + 1. _d 0
130	heimbach	1.1	endif
131			enddo
132	heimbach	1.6	enddo
133	heimbach	1.1	enddo
134
135			objf_temp0(bi,bj) = objf_temp0(bi,bj) + fctile
136
137			enddo
138			enddo
139
140	gforget	1.8	#ifndef ALLOW_SMOOTH_CORREL3D
141	heimbach	1.5	#ifdef ALLOW_SMOOTH_IC_COST_CONTRIBUTION
142
143	heimbach	1.7	call active_read_xyz(
144	heimbach	1.5	& fnamefld, tmpfld3d, irec, doglobalread,
145			& ladinit, optimcycle, mythid, xx_theta_dummy )
146
147			_EXCH_XYZ_R8(tmpfld3d, mythid)
148
149			c-- Loop over this thread's tiles.
150			do bj = jtlo,jthi
151			do bi = itlo,ithi
152
153			fctile = 0. _d 0
154			do k = 1,nr
155	heimbach	1.6	do j = jmin,jmax
156	heimbach	1.5	do i = imin,imax
157			if (_hFacC(i,j,k,bi,bj) .ne. 0.) then
158			tmpx =
159			& ( tmpfld3d(i+2,j,k,bi,bj)-tmpfld3d(i+1,j,k,bi,bj) )
160			& maskW(i+1,j,k,bi,bj)maskW(i+2,j,k,bi,bj)
161			& + ( tmpfld3d(i+1,j,k,bi,bj)-tmpfld3d(i,j,k,bi,bj) )
162			& *maskW(i+1,j,k,bi,bj)
163			& + ( tmpfld3d(i,j+2,k,bi,bj)-tmpfld3d(i,j+1,k,bi,bj) )
164			& maskS(i,j+1,k,bi,bj)maskS(i,j+2,k,bi,bj)
165			& + ( tmpfld3d(i,j+1,k,bi,bj)-tmpfld3d(i,j,k,bi,bj) )
166			& *maskS(i,j+1,k,bi,bj)
167	gforget	1.8	if ( ABS(R_low(i,j,bi,bj)) .LT. 100. )
168	heimbach	1.5	& tmpx = tmpx*ABS(R_low(i,j,bi,bj))/100.
169			fctile = fctile
170			& + wthetaLev(i,j,k,bi,bj)*cosphi(i,j,bi,bj)
171			* 0.0161lengthscale/4.0
172			& tmpxtmpx
173			endif
174			enddo
175	heimbach	1.6	enddo
176	heimbach	1.5	enddo
177
178	heimbach	1.6	objf_temp0smoo(bi,bj) = objf_temp0smoo(bi,bj) + fctile
179	heimbach	1.5
180			enddo
181			enddo
182	heimbach	1.1	#endif
183	gforget	1.8	#endif
184	heimbach	1.1
185			#endif
186
187			return
188			end
189
190