pkg/cost/cost_ssh_mean.F

C $Header: /u/gcmpack/development/heimbach/ecco_env/pkg/cost/cost_ssh_mean.F,v 1.1 2001/06/12 19:13:37 ralf 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
c     changed: Ralf Giering Ralf.Giering@FastOpt.de 12-Jun-2001
c
c              - totally rewrite for parallel processing
c
c     ==================================================================
c     SUBROUTINE cost_ssh_mean
c     ==================================================================

      implicit none

c     == global variables ==

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

#include "EESUPPORT.h"

#include "cost.h"

#ifdef ALLOW_EGM96_ERROR_COV
# include "cost_sph.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

#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

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 = jmin,jmax
               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
      call exch_all_2d_rl( misfit, misfitgl, mythid )

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

      _BEGIN_MASTER( mythid )
      print *
      print *, ' cost_SSH: grid starts at point ', joffset
     &       , ' of full earth.'
      _END_MASTER( 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

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)

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

            do j = jmin,jmax
              do i = jmin,jmax
                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 = jmin,jmax
              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)
ce            --> check masks!
            enddo
          enddo
        enddo
      enddo

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

      end
1	C $Header: /u/gcmpack/development/heimbach/ecco_env/pkg/cost/cost_ssh_mean.F,v 1.1 2001/06/12 19:13:37 ralf Exp $
2
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
23	c changed: Ralf Giering Ralf.Giering@FastOpt.de 12-Jun-2001
24	c
25	c - totally rewrite for parallel processing
26	c
27	c ==================================================================
28	c SUBROUTINE cost_ssh_mean
29	c ==================================================================
30
31	implicit none
32
33	c == global variables ==
34
35	#include "EEPARAMS.h"
36	#include "SIZE.h"
37	#include "PARAMS.h"
38
39	#include "EESUPPORT.h"
40
41	#include "cost.h"
42
43	#ifdef ALLOW_EGM96_ERROR_COV
44	# include "cost_sph.h"
45	#endif
46
47	c == routine arguments ==
48
49	_RL psmean ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
50	_RL offset
51	_RL costmean
52	integer mythid
53
54	c == local variables ==
55
56	integer bi,bj
57	integer i,j
58	integer itlo,ithi
59	integer jtlo,jthi
60	integer jmin,jmax
61	integer imin,imax
62
63	#ifdef ALLOW_EGM96_ERROR_COV
64	_RL cphi
65	_RL xmean
66	_RL shc( ncshc )
67	_RL misfit ( 1-olx:snx+olx, 1-oly:sny+oly, nsx,nsy )
68	_RL misfitgl( 1-olx:snx+olx, 1-oly:sny+oly, nsx,nsy, npx,npy )
69
70	_RL pwork ( (lshcmax+1)*(lshcmax+2)/2 )
71	_RL diffearth ( lonshc, latshc )
72	integer joffset
73	integer ipx, ipy
74	integer iglobe, jglobe
75	integer iproc
76	integer mpirc
77	integer mpicrd(2)
78	#endif
79
80	_RL diff
81	Real*8 sumc
82
83	c == end of interface ==
84
85	c-- Initialise local variables.
86	jtlo = mybylo(mythid)
87	jthi = mybyhi(mythid)
88	itlo = mybxlo(mythid)
89	ithi = mybxhi(mythid)
90
91	jmin = 1
92	jmax = sny
93	imin = 1
94	imax = snx
95
96	costmean = 0.
97
98	#ifdef ALLOW_EGM96_ERROR_COV
99
100	c-- compute misfit on local domain
101	do bj = jtlo,jthi
102	do bi = itlo,ithi
103	do j = jmin,jmax
104	do i = jmin,jmax
105	if ( tpmeanmask(i,j,bi,bj) .gt. 0. ) then
106	misfit(i,j,bi,bj) =
107	& tpmean(i,j,bi,bj) - offset - psmean(i,j,bi,bj)
108	else
109	misfit(i,j,bi,bj) = 0.
110	endif
111	enddo
112	enddo
113	enddo
114	enddo
115
116	c-- communicate to get misfit on full 2d model surface grid
117	call exch_all_2d_rl( misfit, misfitgl, mythid )
118
119	c-- set meridional index offset,
120	c-- it is non zero if the grid does no reach the poles
121	joffset = (latshc - ny)/2
122
123	_BEGIN_MASTER( mythid )
124	print *
125	print *, ' cost_SSH: grid starts at point ', joffset
126	& , ' of full earth.'
127	_END_MASTER( mythid )
128
129	c-- preset field
130	c-- necessary if the grid does no reach the poles
131	do j = 1, latshc
132	do i = 1, lonshc
133	diffearth(i,j) = 0.
134	enddo
135	enddo
136
137	c-- -interpolate- from model grid to regular grid
138	c-- so far we assume that the model grid is already regular
139	_BEGIN_MASTER( mythid )
140
141	do iproc = 1, nprocs
142
143	c-- get coordinates of processor (iporc-1)
144	call MPI_Cart_coords(
145	I MPI_COMM_MODEL, iproc-1, 2, mpicrd
146	O , mpirc
147	& )
148
149	ipx = 1 + mpicrd(1)
150	ipy = 1 + mpicrd(2)
151
152	do bj = jtlo,jthi
153	do bi = itlo,ithi
154
155	do j = jmin,jmax
156	do i = jmin,jmax
157	jglobe = joffset+ j + sNy(bj-1) + sNynSy*(ipy-1)
158	iglobe = i + sNx(bi-1) + sNxnSx*(ipx-1)
159
160	diffearth(iglobe,jglobe) = misfitgl(i,j,bi,bj,ipx,ipy)
161	enddo
162	enddo
163
164	enddo
165	enddo
166
167	enddo
168
169	c-- Project regular grid misfit onto sperical harmonics
170	call shc4grid(
171	I lshcmax
172	O , shc
173	I , latshc, lonshc
174	I , diffearth
175	W , pwork
176	& )
177
178	c-- Remove the C(0,0) component, i.e. the global mean.
179	shc(1) = 0.
180
181	C-- Compute the cost function for the mean SSH.
182	call cost_geoid(
183	O costmean
184	I , shc
185	I , mythid
186	& )
187
188	_END_MASTER( mythid )
189
190	_BARRIER
191
192	#else
193
194	c-- Compute cost function for SSH by using the diagonal
195	c-- of the error covariance only.
196	c--
197	c-- Note: wp is assumed to include latitude dependence
198	c-- of error due to meridian convergence;
199	c-- --> no weighting by cosphi.
200
201	sumc = 0.
202	do bj = jtlo,jthi
203	do bi = itlo,ithi
204	do j = jmin,jmax
205	do i = jmin,jmax
206	diff = psmean(i,j,bi,bj) - tpmean(i,j,bi,bj) - offset
207	sumc = sumc + diff*diff
208	& * wp(i,j,bi,bj)*tpmeanmask(i,j,bi,bj)
209	ce --> check masks!
210	enddo
211	enddo
212	enddo
213	enddo
214
215	c-- Do the global summation.
216	_GLOBAL_SUM_R8( sumc, mythid )
217	costmean = sumc
218	#endif
219
220	end