pkg/cost/cost_atlantic_heat.F

C $Header: /u/gcmpack/MITgcm/pkg/cost/cost_atlantic_heat.F,v 1.8 2008/01/31 16:55:29 heimbach Exp $
C $Name:  $

#include "COST_CPPOPTIONS.h"

      subroutine cost_atlantic_heat( myThid )
C     /==========================================================\
C     | subroutine cost_atlantic_heat                            |
C     | o This routine computes the meridional heat transport.   |
C     |   The current indices are for North Atlantic 29N         |
C     |   2x2 global setup.                                      |
C     \==========================================================/
       implicit none

C     == Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "cost.h"

C     ======== Routine arguments ======================
C     myThid - Thread number for this instance of the routine.
      integer myThid

#ifdef ALLOW_COST_ATLANTIC_HEAT
C     ========= Local variables =========================
      integer    isecbeg      , isecend      , jsec
      integer    jsecbeg      , jsecend      , isec
      integer    kmaxdepth
      integer i, j, k
      integer ig, jg
      integer bi, bj
      _RL     locfc
      _RL     uVel_bar(Nr), vVel_bar(Nr), theta_bar(Nr)
      _RL     thetaUvel_bar(Nr), thetaVvel_bar(Nr)
      _RL     countU(Nr), countV(Nr), countT(Nr)
      _RL     countTU(Nr), countTV(Nr)
      _RL     petawatt
      _RL     sum
      parameter( petawatt = 1. _d +15 )

C     80W - 0W at 24N
      parameter( isecbeg = 69, isecend = 87, jsec = 28 )
cph      parameter( isecbeg = 70, isecend = 90, jsec = 30 )
      parameter( jsecbeg = 10, jsecend = 27, isec = 59 )
#ifdef ALLOW_COST_ATLANTIC_HEAT_DOMASS
      parameter ( kmaxdepth = 8 )
#else
      parameter ( kmaxdepth = 14 )
#endif

C------------------------------------------------------
C     Accumulate meridionally integrated transports
C     Note bar(V)*bar(T) not bar(VT)
C     Attention pYFaceA [m^2*gravity*rhoConst]
C------------------------------------------------------

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)

        locfc = 0.0
        sum = 0.0

#define MERID_TRANSPORT

#ifdef MERID_TRANSPORT

#undef ENERGYNORM
#ifdef ENERGYNORM

        do i=1,sNx
         ig = myXGlobalLo-1+(bi-1)*sNx+i
         if ((ig .ge. isecbeg) .and. (ig .le. isecend)) then
          sum = 0.0
          do k = 1, kmaxdepth
           sum = sum
     &            + (vVel(i,j,k,bi,bj) * maskS(i,j,k,bi,bj)
     &            * drF(k))**2
          end do
          locfc = locfc + sum*dxG(i,j,bi,bj)
         end if
        end do

#else

        do j=1,sNy
         jg = myYGlobalLo-1+(bj-1)*sNy+j
         if (jg .eq. jsec) then

          do k = 1, Nr
           vVel_bar(k) = 0.0
           thetaVvel_bar(k) = 0.0
           countV(k) = 0.0
           countTV(k) = 0.0
           do i=1,sNx
            ig = myXGlobalLo-1+(bi-1)*sNx+i
c
            if ((ig .ge. isecbeg) .and. (ig .le. isecend)) then
               vVel_bar(k) = vVel_bar(k)
     &              + cMeanVVel(i,j,k,bi,bj)*dxG(i,j,bi,bj)
     &                 *maskS(i,j,k,bi,bj)

               thetaVvel_bar(k) = thetaVvel_bar(k)
     &            + cMeanThetaVVel(i,j,k,bi,bj)*dxG(i,j,bi,bj)
     &                 *maskS(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)

             countTV(k) = countTV(k) +
     &            maskS(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)
             countV(k) = countV(k) +
     &            maskS(i,j,k,bi,bj)
            end if

           enddo
          enddo
c
          do k = 1, Nr
#ifdef ALLOW_COST_ATLANTIC_HEAT_DOMASS
           if ( k .LE. kmaxdepth .AND. countV(k) .NE. 0) then
            sum = sum
     &            + vVel_bar(k)*drF(k)/countV(k)
           end if
#else
           if ( k .LE. kmaxdepth .AND. countTV(k) .NE. 0) then
            sum = sum
     &            + thetaVVel_bar(k)*drF(k)/countTV(k)
           end if
#endif
          end do

#endif /* ENERGYNORM */

#else

cph need to change this part to go from
cph \bar{u}*\bar{T} to \bar{u*T}
cph (required store dir. are now in place)

        do i=1,sNx
         ig = myXGlobalLo-1+(bi-1)*sNx+i
         if (ig .eq. isec) then

          do k = 1, Nr
           uVel_bar(k) = 0.0
           theta_bar(k) = 0.0
           countT(k) = 0.0
           countU(k) = 0.0
           do j=1,sNy
            jg = myYGlobalLo-1+(bj-1)*sNy+j
c
            if ((jg .ge. jsecbeg) .and. (jg .le. jsecend)) then
             uVel_bar(k) = uVel_bar(k)
     &                      + cMeanUVel(i,j,k,bi,bj)
     &                      *maskW(i,j,k,bi,bj)
     &                      *maskC(i,j,k,bi,bj)*maskC(i-1,j,k,bi,bj)
             theta_bar(k) = theta_bar(k) +
     &            0.5*( cMeanTheta(i,j,k,bi,bj)
     &                 +cMeanTheta(i-,j,k,bi,bj) )
     &                 *maskW(i,j,k,bi,bj)*dyG(i,j,bi,bj)
     &                 *maskC(i,j,k,bi,bj)*maskC(i-1,j,k,bi,bj)
             countT(k) = countT(k) + maskW(i,j,k,bi,bj)
     &                   *maskC(i,j,k,bi,bj)*maskC(i-1,j,k,bi,bj)
             countU(k) = countU(k) + maskW(i,j,k,bi,bj)
     &                   *maskC(i,j,k,bi,bj)*maskC(i-1,j,k,bi,bj)
            end if

           enddo
          enddo
c
          do k = 1, Nr
           if ( k .LE. kmaxdepth .AND.
     &          countT(k) .NE. 0 .AND. countU(k) .NE. 0) then
            sum = sum
     &            + uVel_bar(k) * theta_bar(k) * drF(k)
     &            / ( countT(k) * countU(k) )
           end if
          end do

#endif

         end if
        end do

#ifdef ALLOW_COST_ATLANTIC_HEAT_DOMASS
        objf_atl(bi,bj) =
     &     sum*1.E-6
#else
        objf_atl(bi,bj) =
     &     sum*HeatCapacity_Cp*rhoConst/petawatt
#endif

c--   end of bi,bj loop
       end do
      end do

#endif

      end
1	utke	1.9	C $Header: /u/gcmpack/MITgcm/pkg/cost/cost_atlantic_heat.F,v 1.8 2008/01/31 16:55:29 heimbach Exp $
2	jmc	1.7	C $Name: $
3	heimbach	1.1
4			#include "COST_CPPOPTIONS.h"
5
6			subroutine cost_atlantic_heat( myThid )
7			C /==========================================================\
8			C \| subroutine cost_atlantic_heat \|
9			C \| o This routine computes the meridional heat transport. \|
10			C \| The current indices are for North Atlantic 29N \|
11			C \| 2x2 global setup. \|
12			C \==========================================================/
13			implicit none
14
15			C == Global variables ===
16			#include "SIZE.h"
17			#include "EEPARAMS.h"
18			#include "PARAMS.h"
19			#include "GRID.h"
20			#include "DYNVARS.h"
21			#include "cost.h"
22
23			C ======== Routine arguments ======================
24			C myThid - Thread number for this instance of the routine.
25			integer myThid
26
27			#ifdef ALLOW_COST_ATLANTIC_HEAT
28			C ========= Local variables =========================
29			integer isecbeg , isecend , jsec
30	heimbach	1.4	integer jsecbeg , jsecend , isec
31	heimbach	1.1	integer kmaxdepth
32			integer i, j, k
33			integer ig, jg
34			integer bi, bj
35			_RL locfc
36	heimbach	1.6	_RL uVel_bar(Nr), vVel_bar(Nr), theta_bar(Nr)
37			_RL thetaUvel_bar(Nr), thetaVvel_bar(Nr)
38	heimbach	1.4	_RL countU(Nr), countV(Nr), countT(Nr)
39	heimbach	1.6	_RL countTU(Nr), countTV(Nr)
40	heimbach	1.1	_RL petawatt
41			_RL sum
42	utke	1.9	parameter( petawatt = 1. _d +15 )
43	heimbach	1.1
44			C 80W - 0W at 24N
45	heimbach	1.4	parameter( isecbeg = 69, isecend = 87, jsec = 28 )
46			cph parameter( isecbeg = 70, isecend = 90, jsec = 30 )
47			parameter( jsecbeg = 10, jsecend = 27, isec = 59 )
48	heimbach	1.8	#ifdef ALLOW_COST_ATLANTIC_HEAT_DOMASS
49			parameter ( kmaxdepth = 8 )
50			#else
51			parameter ( kmaxdepth = 14 )
52			#endif
53	heimbach	1.1
54			C------------------------------------------------------
55			C Accumulate meridionally integrated transports
56			C Note bar(V)*bar(T) not bar(VT)
57	mlosch	1.3	C Attention pYFaceA [m^2gravityrhoConst]
58	heimbach	1.1	C------------------------------------------------------
59
60			DO bj=myByLo(myThid),myByHi(myThid)
61			DO bi=myBxLo(myThid),myBxHi(myThid)
62
63			locfc = 0.0
64	heimbach	1.2	sum = 0.0
65
66	heimbach	1.4	#define MERID_TRANSPORT
67
68			#ifdef MERID_TRANSPORT
69	heimbach	1.1
70			#undef ENERGYNORM
71			#ifdef ENERGYNORM
72
73			do i=1,sNx
74			ig = myXGlobalLo-1+(bi-1)*sNx+i
75			if ((ig .ge. isecbeg) .and. (ig .le. isecend)) then
76			sum = 0.0
77			do k = 1, kmaxdepth
78			sum = sum
79			& + (vVel(i,j,k,bi,bj) * maskS(i,j,k,bi,bj)
80			& * drF(k))**2
81			end do
82			locfc = locfc + sum*dxG(i,j,bi,bj)
83			end if
84			end do
85
86			#else
87
88	heimbach	1.4	do j=1,sNy
89			jg = myYGlobalLo-1+(bj-1)*sNy+j
90			if (jg .eq. jsec) then
91
92	heimbach	1.1	do k = 1, Nr
93			vVel_bar(k) = 0.0
94	heimbach	1.6	thetaVvel_bar(k) = 0.0
95	heimbach	1.2	countV(k) = 0.0
96	heimbach	1.6	countTV(k) = 0.0
97	heimbach	1.1	do i=1,sNx
98			ig = myXGlobalLo-1+(bi-1)*sNx+i
99	heimbach	1.2	c
100	heimbach	1.1	if ((ig .ge. isecbeg) .and. (ig .le. isecend)) then
101	jmc	1.7	vVel_bar(k) = vVel_bar(k)
102	heimbach	1.8	& + cMeanVVel(i,j,k,bi,bj)*dxG(i,j,bi,bj)
103			& *maskS(i,j,k,bi,bj)
104	heimbach	1.6
105			thetaVvel_bar(k) = thetaVvel_bar(k)
106			& + cMeanThetaVVel(i,j,k,bi,bj)*dxG(i,j,bi,bj)
107			& maskS(i,j,k,bi,bj)maskC(i,j,k,bi,bj)
108
109	jmc	1.7	countTV(k) = countTV(k) +
110	heimbach	1.6	& maskS(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)
111	jmc	1.7	countV(k) = countV(k) +
112	heimbach	1.6	& maskS(i,j,k,bi,bj)
113	heimbach	1.1	end if
114	heimbach	1.2
115	heimbach	1.1	enddo
116			enddo
117	heimbach	1.2	c
118			do k = 1, Nr
119	heimbach	1.6	#ifdef ALLOW_COST_ATLANTIC_HEAT_DOMASS
120			if ( k .LE. kmaxdepth .AND. countV(k) .NE. 0) then
121			sum = sum
122			& + vVel_bar(k)*drF(k)/countV(k)
123			end if
124			#else
125			if ( k .LE. kmaxdepth .AND. countTV(k) .NE. 0) then
126	heimbach	1.2	sum = sum
127	heimbach	1.6	& + thetaVVel_bar(k)*drF(k)/countTV(k)
128	heimbach	1.2	end if
129	heimbach	1.6	#endif
130	heimbach	1.1	end do
131
132	heimbach	1.4	#endif /* ENERGYNORM */
133
134			#else
135
136	jmc	1.7	cph need to change this part to go from
137	heimbach	1.6	cph \bar{u}\bar{T} to \bar{uT}
138			cph (required store dir. are now in place)
139
140	heimbach	1.4	do i=1,sNx
141			ig = myXGlobalLo-1+(bi-1)*sNx+i
142			if (ig .eq. isec) then
143
144			do k = 1, Nr
145			uVel_bar(k) = 0.0
146			theta_bar(k) = 0.0
147			countT(k) = 0.0
148			countU(k) = 0.0
149			do j=1,sNy
150			jg = myYGlobalLo-1+(bj-1)*sNy+j
151			c
152			if ((jg .ge. jsecbeg) .and. (jg .le. jsecend)) then
153			uVel_bar(k) = uVel_bar(k)
154			& + cMeanUVel(i,j,k,bi,bj)
155			& *maskW(i,j,k,bi,bj)
156			& maskC(i,j,k,bi,bj)maskC(i-1,j,k,bi,bj)
157			theta_bar(k) = theta_bar(k) +
158			& 0.5*( cMeanTheta(i,j,k,bi,bj)
159			& +cMeanTheta(i-,j,k,bi,bj) )
160			& maskW(i,j,k,bi,bj)dyG(i,j,bi,bj)
161			& maskC(i,j,k,bi,bj)maskC(i-1,j,k,bi,bj)
162			countT(k) = countT(k) + maskW(i,j,k,bi,bj)
163			& maskC(i,j,k,bi,bj)maskC(i-1,j,k,bi,bj)
164			countU(k) = countU(k) + maskW(i,j,k,bi,bj)
165			& maskC(i,j,k,bi,bj)maskC(i-1,j,k,bi,bj)
166			end if
167
168			enddo
169			enddo
170			c
171			do k = 1, Nr
172			if ( k .LE. kmaxdepth .AND.
173			& countT(k) .NE. 0 .AND. countU(k) .NE. 0) then
174			sum = sum
175			& + uVel_bar(k) * theta_bar(k) * drF(k)
176			& / ( countT(k) * countU(k) )
177			end if
178			end do
179
180	heimbach	1.1	#endif
181
182			end if
183			end do
184
185	heimbach	1.6	#ifdef ALLOW_COST_ATLANTIC_HEAT_DOMASS
186	jmc	1.7	objf_atl(bi,bj) =
187	heimbach	1.6	& sum*1.E-6
188			#else
189	jmc	1.7	objf_atl(bi,bj) =
190	heimbach	1.4	& sumHeatCapacity_CprhoConst/petawatt
191	heimbach	1.6	#endif
192	heimbach	1.1
193			c-- end of bi,bj loop
194			end do
195			end do
196
197			#endif
198	jmc	1.7
199	heimbach	1.1	end