pkg/ecco/cost_ssh_mean.F

C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_ssh_mean.F,v 1.4 2005/05/05 23:47:57 heimbach Exp $

#include "COST_CPPOPTIONS.h"


      subroutine cost_ssh_mean(
     I                     psmean, offset
     O                   , costmean
     I                   , mythid
     &                   )

c     ==================================================================
c     SUBROUTINE cost_ssh_mean
c     ==================================================================
c
c     o Evaluate cost function contribution of sea surface height.
c       using of geoid error covariances requires regular model grid
c     o TODO: interpolate model grid to regular grid
c             check mask
c
c     started: Detlef Stammer, Ralf Giering Jul-1996
c     changed: Ralf Giering Ralf.Giering@FastOpt.de 12-Jun-2001
c              - totally rewrite for parallel processing
c              heimbach@mit.edu 13-Mar-2002
c              - several wrong i-loop boundaries (spotted by G. Gebbie)
c              - nprocs need to be replaced by nPx*nPy
c              - geoid error does not work as of now
c              heimbach@mit.edu 05-May-2005
c              - debugged and restructuted
c
c     ==================================================================
c     SUBROUTINE cost_ssh_mean
c     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"

#include "ecco_cost.h"
#include "optim.h"

#ifdef ALLOW_EGM96_ERROR_COV
# if (defined (ALLOW_USE_MPI) || defined (ALWAYS_USE_MPI))
#  include "EESUPPORT.h"
# else
      INTEGER nProcs
      PARAMETER (nProcs=1)
# endif
# include "sphere.h"
#endif

c     == routine arguments ==

      _RL     psmean ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL     offset
      _RL     costmean
      integer mythid

c     == local variables ==

      integer bi,bj
      integer i,j
      integer itlo,ithi
      integer jtlo,jthi
      integer jmin,jmax
      integer imin,imax

cph(
      _RL diagnosfld(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
cph)

#ifdef ALLOW_EGM96_ERROR_COV
      _RL cphi
      _RL xmean
      _RL shc( ncshc )
      _RL misfit  ( 1-olx:snx+olx, 1-oly:sny+oly, nsx,nsy )
      _RL misfitgl( 1-olx:snx+olx, 1-oly:sny+oly, nsx,nsy, npx,npy )

      _RL pwork     ( (lshcmax+1)*(lshcmax+2)/2 )
      _RL diffearth ( lonshc, latshc )
      integer joffset
      integer ipx, ipy
      integer iglobe, jglobe
      integer iproc
      integer mpirc
      integer mpicrd(2)
#endif

      _RL diff
      Real*8 sumc
      character*(max_len_mbuf) msgbuf

c     == end of interface ==

c--   Initialise local variables.
      jtlo = mybylo(mythid)
      jthi = mybyhi(mythid)
      itlo = mybxlo(mythid)
      ithi = mybxhi(mythid)

      jmin = 1
      jmax = sny
      imin = 1
      imax = snx

      costmean   = 0.

#ifdef ALLOW_EGM96_ERROR_COV

c--   compute misfit on local domain
      do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = jmin,jmax
            do i = imin,imax
               if ( tpmeanmask(i,j,bi,bj) .gt. 0. ) then
                 misfit(i,j,bi,bj) =
     &           tpmean(i,j,bi,bj) - offset - psmean(i,j,bi,bj)
              else
                 misfit(i,j,bi,bj) = 0.
              endif
            enddo
          enddo
        enddo
      enddo

c--   communicate to get misfit on full 2d model surface grid
#if (defined (ALLOW_USE_MPI) || defined (ALWAYS_USE_MPI))
      call exch_all_2d_rl( misfit, misfitgl, mythid )
#else
      do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = jmin,jmax
            do i = imin,imax
               misfitgl(i,j,bi,bj,1,1) = misfit(i,j,bi,bj)
            enddo
          enddo
        enddo
      enddo
#endif

c--   set meridional index offset,
c--   it is non zero if the grid does no reach the poles
      joffset = (latshc - ny)/2

      write(msgbuf,'(a,I10)')
     &     'cost_ssh_mean: grid starts at point ', joffset
      call print_message( msgbuf, standardmessageunit,
     &                    SQUEEZE_RIGHT , mythid)

c--   preset field
c--   necessary if the grid does no reach the poles
      do j = 1, latshc
        do i = 1, lonshc
          diffearth(i,j) = 0.
        enddo
      enddo

c--   -interpolate- from model grid to regular grid
c--   so far we assume that the model grid is already regular
      _BEGIN_MASTER( mythid )

      do iproc = 1, nprocs

#if (defined (ALLOW_USE_MPI) || defined (ALWAYS_USE_MPI))
c--     get coordinates of processor (iporc-1)
        call MPI_Cart_coords(
     I          MPI_COMM_MODEL, iproc-1, 2, mpicrd
     O        , mpirc
     &        )
        ipx = 1 + mpicrd(1)
        ipy = 1 + mpicrd(2)
#else
        ipx = 1
        ipy = 1
#endif

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

            do j = jmin,jmax
              do i = imin,imax
                jglobe = joffset+ j + sNy*(bj-1) + sNy*nSy*(ipy-1)
                iglobe =          i + sNx*(bi-1) + sNx*nSx*(ipx-1)

                diffearth(iglobe,jglobe) = misfitgl(i,j,bi,bj,ipx,ipy)
              enddo
            enddo

          enddo
        enddo

      enddo

c--   Project regular grid misfit onto sperical harmonics
      call shc4grid(
     I               lshcmax
     O             , shc
     I             , latshc, lonshc
     I             , diffearth
     W             , pwork
     &             )

c--   Remove the C(0,0) component, i.e. the global mean.
      shc(1) = 0.

C--   Compute the cost function for the mean SSH.
      call cost_geoid(
     O                 costmean
     I               , shc
     I               , mythid
     &               )

      _END_MASTER( mythid )

      _BARRIER

#else

c--   Compute cost function for SSH by using the diagonal
c--   of the error covariance only.
c--
c--   Note: wp is assumed to include latitude dependence
c--         of error due to meridian convergence;
c--         --> no weighting by cosphi.

      sumc = 0.
      do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = jmin,jmax
            do i = imin,imax
              diff = psmean(i,j,bi,bj) - tpmean(i,j,bi,bj) + offset
              sumc = sumc + diff*diff
     &                    * wp(i,j,bi,bj)*tpmeanmask(i,j,bi,bj)
              if ( wp(i,j,bi,bj)*tpmeanmask(i,j,bi,bj) .ne. 0. )
     &             num_hmean = num_hmean + 1. _d 0
              diagnosfld(i,j,bi,bj) = diff*diff
     &                    * wp(i,j,bi,bj)*tpmeanmask(i,j,bi,bj)
            enddo
          enddo
        enddo
      enddo

cph(
      CALL WRITE_FLD_XY_RL( 'DiagnosSSHmean', ' ', diagnosfld,
     &                           optimcycle, mythid )
cph)

c--   Do the global summation.
      _GLOBAL_SUM_R8( sumc, mythid )
      _GLOBAL_SUM_R8( num_hmean, mythid )
      costmean = sumc

#endif

      end
1	heimbach	1.5	C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_ssh_mean.F,v 1.4 2005/05/05 23:47:57 heimbach Exp $
2	heimbach	1.1
3			#include "COST_CPPOPTIONS.h"
4
5
6			subroutine cost_ssh_mean(
7			I psmean, offset
8			O , costmean
9			I , mythid
10			& )
11
12			c ==================================================================
13			c SUBROUTINE cost_ssh_mean
14			c ==================================================================
15			c
16			c o Evaluate cost function contribution of sea surface height.
17			c using of geoid error covariances requires regular model grid
18			c o TODO: interpolate model grid to regular grid
19			c check mask
20			c
21			c started: Detlef Stammer, Ralf Giering Jul-1996
22			c changed: Ralf Giering Ralf.Giering@FastOpt.de 12-Jun-2001
23			c - totally rewrite for parallel processing
24			c heimbach@mit.edu 13-Mar-2002
25			c - several wrong i-loop boundaries (spotted by G. Gebbie)
26			c - nprocs need to be replaced by nPx*nPy
27			c - geoid error does not work as of now
28	heimbach	1.4	c heimbach@mit.edu 05-May-2005
29			c - debugged and restructuted
30	heimbach	1.1	c
31			c ==================================================================
32			c SUBROUTINE cost_ssh_mean
33			c ==================================================================
34
35			implicit none
36
37			c == global variables ==
38
39			#include "EEPARAMS.h"
40			#include "SIZE.h"
41			#include "PARAMS.h"
42
43			#include "ecco_cost.h"
44	heimbach	1.5	#include "optim.h"
45	heimbach	1.1
46			#ifdef ALLOW_EGM96_ERROR_COV
47	heimbach	1.4	# if (defined (ALLOW_USE_MPI) \|\| defined (ALWAYS_USE_MPI))
48			# include "EESUPPORT.h"
49			# else
50			INTEGER nProcs
51			PARAMETER (nProcs=1)
52			# endif
53			# include "sphere.h"
54	heimbach	1.1	#endif
55
56			c == routine arguments ==
57
58			_RL psmean ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
59			_RL offset
60			_RL costmean
61			integer mythid
62
63			c == local variables ==
64
65			integer bi,bj
66			integer i,j
67			integer itlo,ithi
68			integer jtlo,jthi
69			integer jmin,jmax
70			integer imin,imax
71
72	heimbach	1.5	cph(
73			_RL diagnosfld(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
74			cph)
75
76	heimbach	1.1	#ifdef ALLOW_EGM96_ERROR_COV
77			_RL cphi
78			_RL xmean
79			_RL shc( ncshc )
80			_RL misfit ( 1-olx:snx+olx, 1-oly:sny+oly, nsx,nsy )
81			_RL misfitgl( 1-olx:snx+olx, 1-oly:sny+oly, nsx,nsy, npx,npy )
82
83			_RL pwork ( (lshcmax+1)*(lshcmax+2)/2 )
84			_RL diffearth ( lonshc, latshc )
85			integer joffset
86			integer ipx, ipy
87			integer iglobe, jglobe
88			integer iproc
89			integer mpirc
90			integer mpicrd(2)
91			#endif
92
93			_RL diff
94			Real*8 sumc
95	heimbach	1.4	character*(max_len_mbuf) msgbuf
96	heimbach	1.1
97			c == end of interface ==
98
99			c-- Initialise local variables.
100			jtlo = mybylo(mythid)
101			jthi = mybyhi(mythid)
102			itlo = mybxlo(mythid)
103			ithi = mybxhi(mythid)
104
105			jmin = 1
106			jmax = sny
107			imin = 1
108			imax = snx
109
110			costmean = 0.
111
112			#ifdef ALLOW_EGM96_ERROR_COV
113
114			c-- compute misfit on local domain
115			do bj = jtlo,jthi
116			do bi = itlo,ithi
117			do j = jmin,jmax
118			do i = imin,imax
119			if ( tpmeanmask(i,j,bi,bj) .gt. 0. ) then
120			misfit(i,j,bi,bj) =
121			& tpmean(i,j,bi,bj) - offset - psmean(i,j,bi,bj)
122			else
123			misfit(i,j,bi,bj) = 0.
124			endif
125			enddo
126			enddo
127			enddo
128			enddo
129
130			c-- communicate to get misfit on full 2d model surface grid
131	heimbach	1.4	#if (defined (ALLOW_USE_MPI) \|\| defined (ALWAYS_USE_MPI))
132	heimbach	1.1	call exch_all_2d_rl( misfit, misfitgl, mythid )
133	heimbach	1.4	#else
134			do bj = jtlo,jthi
135			do bi = itlo,ithi
136			do j = jmin,jmax
137			do i = imin,imax
138			misfitgl(i,j,bi,bj,1,1) = misfit(i,j,bi,bj)
139			enddo
140			enddo
141			enddo
142			enddo
143			#endif
144	heimbach	1.1
145			c-- set meridional index offset,
146			c-- it is non zero if the grid does no reach the poles
147			joffset = (latshc - ny)/2
148
149	heimbach	1.4	write(msgbuf,'(a,I10)')
150			& 'cost_ssh_mean: grid starts at point ', joffset
151			call print_message( msgbuf, standardmessageunit,
152			& SQUEEZE_RIGHT , mythid)
153	heimbach	1.1
154			c-- preset field
155			c-- necessary if the grid does no reach the poles
156			do j = 1, latshc
157			do i = 1, lonshc
158			diffearth(i,j) = 0.
159			enddo
160			enddo
161
162			c-- -interpolate- from model grid to regular grid
163			c-- so far we assume that the model grid is already regular
164			_BEGIN_MASTER( mythid )
165
166			do iproc = 1, nprocs
167
168	heimbach	1.4	#if (defined (ALLOW_USE_MPI) \|\| defined (ALWAYS_USE_MPI))
169	heimbach	1.1	c-- get coordinates of processor (iporc-1)
170			call MPI_Cart_coords(
171			I MPI_COMM_MODEL, iproc-1, 2, mpicrd
172			O , mpirc
173			& )
174			ipx = 1 + mpicrd(1)
175			ipy = 1 + mpicrd(2)
176	heimbach	1.4	#else
177			ipx = 1
178			ipy = 1
179			#endif
180	heimbach	1.1
181			do bj = jtlo,jthi
182			do bi = itlo,ithi
183
184			do j = jmin,jmax
185			do i = imin,imax
186			jglobe = joffset+ j + sNy(bj-1) + sNynSy*(ipy-1)
187			iglobe = i + sNx(bi-1) + sNxnSx*(ipx-1)
188
189			diffearth(iglobe,jglobe) = misfitgl(i,j,bi,bj,ipx,ipy)
190			enddo
191			enddo
192
193			enddo
194			enddo
195
196			enddo
197
198			c-- Project regular grid misfit onto sperical harmonics
199			call shc4grid(
200			I lshcmax
201			O , shc
202			I , latshc, lonshc
203			I , diffearth
204			W , pwork
205			& )
206
207			c-- Remove the C(0,0) component, i.e. the global mean.
208			shc(1) = 0.
209
210			C-- Compute the cost function for the mean SSH.
211			call cost_geoid(
212			O costmean
213			I , shc
214			I , mythid
215			& )
216
217			_END_MASTER( mythid )
218
219			_BARRIER
220
221			#else
222
223			c-- Compute cost function for SSH by using the diagonal
224			c-- of the error covariance only.
225			c--
226			c-- Note: wp is assumed to include latitude dependence
227			c-- of error due to meridian convergence;
228			c-- --> no weighting by cosphi.
229
230			sumc = 0.
231			do bj = jtlo,jthi
232			do bi = itlo,ithi
233			do j = jmin,jmax
234			do i = imin,imax
235			diff = psmean(i,j,bi,bj) - tpmean(i,j,bi,bj) + offset
236			sumc = sumc + diff*diff
237			& * wp(i,j,bi,bj)*tpmeanmask(i,j,bi,bj)
238	heimbach	1.3	if ( wp(i,j,bi,bj)*tpmeanmask(i,j,bi,bj) .ne. 0. )
239			& num_hmean = num_hmean + 1. _d 0
240	heimbach	1.5	diagnosfld(i,j,bi,bj) = diff*diff
241			& * wp(i,j,bi,bj)*tpmeanmask(i,j,bi,bj)
242	heimbach	1.1	enddo
243			enddo
244			enddo
245			enddo
246
247	heimbach	1.5	cph(
248			CALL WRITE_FLD_XY_RL( 'DiagnosSSHmean', ' ', diagnosfld,
249			& optimcycle, mythid )
250			cph)
251
252	heimbach	1.1	c-- Do the global summation.
253			_GLOBAL_SUM_R8( sumc, mythid )
254	heimbach	1.3	_GLOBAL_SUM_R8( num_hmean, mythid )
255	heimbach	1.1	costmean = sumc
256	heimbach	1.5
257	heimbach	1.1	#endif
258
259			end