pkg/ecco/cost_theta.F

C $Header:  $
C $Name:  $

#include "COST_CPPOPTIONS.h"


      subroutine cost_theta( myiter, mytime, mythid )

c     ==================================================================
c     SUBROUTINE cost_theta
c     ==================================================================
c
c     o Evaluate cost function contribution of temperature.
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     changed: Patrick Heimbach heimbach@mit.edu 27-May-2000
c
c              - set ladinit to .true. to initialise adtbar file
c
c     ==================================================================
c     SUBROUTINE cost_theta
c     ==================================================================

      implicit none

c     == global variables ==

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

      integer bi,bj
      integer i,j,k
      integer itlo,ithi
      integer jtlo,jthi
      integer jmin,jmax
      integer imin,imax
      integer irec, irectmp
      integer levmon
      integer levoff
      integer iltheta

      _RL fctile
      _RL fcthread

      _RL cmask (1-olx:snx+olx,1-oly:sny+oly)
      _RL spval
      _RL spmax

      character*(80) fnametheta

      logical doglobalread
      logical ladinit

      character*(MAX_LEN_MBUF) msgbuf
#ifdef GENERIC_BAR_MONTH
      integer mrec, nyears, iyear
#endif

      _RL diagnosfld3d(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)

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

      spval = -1.8
      spmax = 40.

c--   Read tiled data.
      doglobalread = .false.
      ladinit      = .false.

#ifdef ALLOW_THETA_COST_CONTRIBUTION

      if (optimcycle .ge. 0) then
        iltheta = ilnblnk( tbarfile )
        write(fnametheta(1:80),'(2a,i10.10)')
     &    tbarfile(1:iltheta),'.',optimcycle
      endif

      fcthread = 0. _d 0

#ifdef GENERIC_BAR_MONTH
c--   Loop over month
      do irec = 1,min(nmonsrec,12)
         nyears=int((nmonsrec-irec)/12)+1
            do iyear=1,nyears
            mrec=irec+(iyear-1)*12
            irectmp=mrec
c--     Read time averages and the monthly mean data.
            call active_read_xyz( fnametheta, tbar, mrec,
     &           doglobalread, ladinit,
     &           optimcycle, mythid,
     &           xx_tbar_mean_dummy )
            do bj = jtlo,jthi
               do bi = itlo,ithi
                  do k = 1,nr
                     do j = jmin,jmax
                        do i = imin,imax
                           if(iyear.eq.1) then
                              tbar_gen(i,j,k,bi,bj) =tbar(i,j,k,bi,bj)
                           elseif(iyear.eq.nyears) then
                              tbar(i,j,k,bi,bj) =(tbar_gen(i,j,k,bi,bj)
     $                             +tbar(i,j,k,bi,bj))/float(nyears)
                           else
                           tbar_gen(i,j,k,bi,bj) =tbar_gen(i,j,k,bi,bj)
     $                             +tbar(i,j,k,bi,bj)
                           endif
                        enddo
                     enddo
                  enddo
               enddo
            enddo
            enddo
#else
c--   Loop over records.
      do irec = 1,nmonsrec

        irectmp = irec
c--     Read time averages and the monthly mean data.
        call active_read_xyz( fnametheta, tbar, irec,
     &                        doglobalread, ladinit,
     &                        optimcycle, mythid,
     &                        xx_tbar_mean_dummy )
#endif
c--     Determine the month to be read.
        levoff = mod(modelstartdate(1)/100,100)
        levmon = (irectmp-1) + levoff
        levmon = mod(levmon-1,12)+1

        call mdsreadfield( tdatfile, cost_iprec, cost_yftype,
     &                     nr, tdat, levmon, mythid)

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

c--         Loop over the model layers
            fctile = 0. _d 0
            do k = 1,nr

c--           Determine the mask on weights
              do j = jmin,jmax
                do i = imin,imax
                  cmask(i,j) = cosphi(i,j,bi,bj)
                  if (tdat(i,j,k,bi,bj) .eq. 0.) then
                    cmask(i,j) = 0. _d 0
                  else if (tdat(i,j,k,bi,bj) .lt. spval) then
                    cmask(i,j) = 0. _d 0
                  else if (tdat(i,j,k,bi,bj) .gt. spmax) then
                    cmask(i,j) = 0. _d 0
                  endif
                enddo
              enddo

c--           Compute model data misfit and cost function term for
c             the temperature field.
              do j = jmin,jmax
                do i = imin,imax
                  if ( _hFacC(i,j,k,bi,bj) .ne. 0. ) then
                     fctile = fctile +
     &                    (wthetaLev(i,j,k,bi,bj)*cmask(i,j)*
     &                    (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj))*
     &                    (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj)) )
                     if ( wthetaLev(i,j,k,bi,bj)*cmask(i,j) .ne. 0. )
     &                    num_temp(bi,bj) = num_temp(bi,bj) + 1. _d 0
                     diagnosfld3d(i,j,k,bi,bj) =
     &                    (wthetaLev(i,j,k,bi,bj)*cmask(i,j)*
     &                    (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj))*
     &                    (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj)) )
                  else
                     diagnosfld3d(i,j,k,bi,bj) = 0.
                  endif
                enddo
              enddo

            enddo
c--         End of loop over layers.

            call mdswritefield( 'DiagnosCost_ClimTheta',
     &           writeBinaryPrec, globalfiles, 'RL', Nr,
     &           diagnosfld3d, irec, optimcycle, mythid )

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

          enddo
        enddo

      enddo
c--   End of loop over records.

#endif

      return
      end

1	jmc	1.9	C $Header: $
2			C $Name: $
3	heimbach	1.1
4			#include "COST_CPPOPTIONS.h"
5
6
7	heimbach	1.7	subroutine cost_theta( myiter, mytime, mythid )
8	heimbach	1.1
9			c ==================================================================
10	heimbach	1.7	c SUBROUTINE cost_theta
11	heimbach	1.1	c ==================================================================
12			c
13			c o Evaluate cost function contribution of temperature.
14			c
15			c started: Christian Eckert eckert@mit.edu 30-Jun-1999
16			c
17			c changed: Christian Eckert eckert@mit.edu 25-Feb-2000
18			c
19			c - Restructured the code in order to create a package
20			c for the MITgcmUV.
21			c
22			c changed: Patrick Heimbach heimbach@mit.edu 27-May-2000
23			c
24			c - set ladinit to .true. to initialise adtbar file
25			c
26			c ==================================================================
27	heimbach	1.7	c SUBROUTINE cost_theta
28	heimbach	1.1	c ==================================================================
29
30			implicit none
31
32			c == global variables ==
33
34			#include "EEPARAMS.h"
35			#include "SIZE.h"
36	heimbach	1.8	#include "PARAMS.h"
37	heimbach	1.1	#include "GRID.h"
38			#include "DYNVARS.h"
39
40			#include "cal.h"
41			#include "ecco_cost.h"
42			#include "ctrl.h"
43			#include "ctrl_dummy.h"
44			#include "optim.h"
45
46			c == routine arguments ==
47
48			integer myiter
49			_RL mytime
50			integer mythid
51
52			c == local variables ==
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	heimbach	1.7	integer irec, irectmp
61	heimbach	1.1	integer levmon
62			integer levoff
63			integer iltheta
64
65			_RL fctile
66			_RL fcthread
67
68	jmc	1.9	_RL cmask (1-olx:snx+olx,1-oly:sny+oly)
69	heimbach	1.1	_RL spval
70	heimbach	1.5	_RL spmax
71	heimbach	1.1
72			character*(80) fnametheta
73
74			logical doglobalread
75			logical ladinit
76
77			character*(MAX_LEN_MBUF) msgbuf
78			#ifdef GENERIC_BAR_MONTH
79	heimbach	1.4	integer mrec, nyears, iyear
80	jmc	1.9	#endif
81	heimbach	1.8
82			_RL diagnosfld3d(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
83	jmc	1.9
84	heimbach	1.1	c == external functions ==
85
86			integer ilnblnk
87			external ilnblnk
88
89			c == end of interface ==
90
91			jtlo = mybylo(mythid)
92			jthi = mybyhi(mythid)
93			itlo = mybxlo(mythid)
94			ithi = mybxhi(mythid)
95			jmin = 1
96			jmax = sny
97			imin = 1
98			imax = snx
99
100	heimbach	1.4	spval = -1.8
101	heimbach	1.5	spmax = 40.
102	heimbach	1.1
103			c-- Read tiled data.
104			doglobalread = .false.
105			ladinit = .false.
106
107			#ifdef ALLOW_THETA_COST_CONTRIBUTION
108
109			if (optimcycle .ge. 0) then
110			iltheta = ilnblnk( tbarfile )
111			write(fnametheta(1:80),'(2a,i10.10)')
112			& tbarfile(1:iltheta),'.',optimcycle
113			endif
114
115			fcthread = 0. _d 0
116
117			#ifdef GENERIC_BAR_MONTH
118			c-- Loop over month
119	heimbach	1.7	do irec = 1,min(nmonsrec,12)
120	heimbach	1.1	nyears=int((nmonsrec-irec)/12)+1
121			do iyear=1,nyears
122			mrec=irec+(iyear-1)*12
123	heimbach	1.7	irectmp=mrec
124	heimbach	1.1	c-- Read time averages and the monthly mean data.
125			call active_read_xyz( fnametheta, tbar, mrec,
126			& doglobalread, ladinit,
127			& optimcycle, mythid,
128			& xx_tbar_mean_dummy )
129			do bj = jtlo,jthi
130			do bi = itlo,ithi
131			do k = 1,nr
132			do j = jmin,jmax
133			do i = imin,imax
134			if(iyear.eq.1) then
135			tbar_gen(i,j,k,bi,bj) =tbar(i,j,k,bi,bj)
136			elseif(iyear.eq.nyears) then
137			tbar(i,j,k,bi,bj) =(tbar_gen(i,j,k,bi,bj)
138			$ +tbar(i,j,k,bi,bj))/float(nyears)
139			else
140	heimbach	1.7	tbar_gen(i,j,k,bi,bj) =tbar_gen(i,j,k,bi,bj)
141	heimbach	1.1	$ +tbar(i,j,k,bi,bj)
142			endif
143			enddo
144			enddo
145			enddo
146			enddo
147			enddo
148	heimbach	1.7	enddo
149	heimbach	1.1	#else
150			c-- Loop over records.
151			do irec = 1,nmonsrec
152
153	heimbach	1.7	irectmp = irec
154	heimbach	1.1	c-- Read time averages and the monthly mean data.
155			call active_read_xyz( fnametheta, tbar, irec,
156			& doglobalread, ladinit,
157			& optimcycle, mythid,
158			& xx_tbar_mean_dummy )
159			#endif
160			c-- Determine the month to be read.
161			levoff = mod(modelstartdate(1)/100,100)
162	heimbach	1.7	levmon = (irectmp-1) + levoff
163	heimbach	1.1	levmon = mod(levmon-1,12)+1
164
165	jmc	1.9	call mdsreadfield( tdatfile, cost_iprec, cost_yftype,
166	heimbach	1.1	& nr, tdat, levmon, mythid)
167
168			do bj = jtlo,jthi
169			do bi = itlo,ithi
170
171			c-- Loop over the model layers
172			fctile = 0. _d 0
173			do k = 1,nr
174
175			c-- Determine the mask on weights
176			do j = jmin,jmax
177			do i = imin,imax
178	heimbach	1.5	cmask(i,j) = cosphi(i,j,bi,bj)
179	heimbach	1.1	if (tdat(i,j,k,bi,bj) .eq. 0.) then
180			cmask(i,j) = 0. _d 0
181	heimbach	1.5	else if (tdat(i,j,k,bi,bj) .lt. spval) then
182			cmask(i,j) = 0. _d 0
183			else if (tdat(i,j,k,bi,bj) .gt. spmax) then
184	heimbach	1.1	cmask(i,j) = 0. _d 0
185			endif
186			enddo
187			enddo
188
189			c-- Compute model data misfit and cost function term for
190			c the temperature field.
191			do j = jmin,jmax
192			do i = imin,imax
193	heimbach	1.6	if ( _hFacC(i,j,k,bi,bj) .ne. 0. ) then
194	heimbach	1.1	fctile = fctile +
195	heimbach	1.8	& (wthetaLev(i,j,k,bi,bj)cmask(i,j)
196	heimbach	1.1	& (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj))*
197	heimbach	1.5	& (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj)) )
198	heimbach	1.8	if ( wthetaLev(i,j,k,bi,bj)*cmask(i,j) .ne. 0. )
199	heimbach	1.5	& num_temp(bi,bj) = num_temp(bi,bj) + 1. _d 0
200	jmc	1.9	diagnosfld3d(i,j,k,bi,bj) =
201	heimbach	1.8	& (wthetaLev(i,j,k,bi,bj)cmask(i,j)
202			& (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj))*
203			& (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj)) )
204			else
205			diagnosfld3d(i,j,k,bi,bj) = 0.
206	heimbach	1.1	endif
207			enddo
208			enddo
209
210			enddo
211			c-- End of loop over layers.
212
213	heimbach	1.8	call mdswritefield( 'DiagnosCost_ClimTheta',
214			& writeBinaryPrec, globalfiles, 'RL', Nr,
215			& diagnosfld3d, irec, optimcycle, mythid )
216
217	heimbach	1.1	fcthread = fcthread + fctile
218			objf_temp(bi,bj) = objf_temp(bi,bj) + fctile
219
220			enddo
221			enddo
222
223			enddo
224			c-- End of loop over records.
225
226			#endif
227
228			return
229			end
230