forked from altMITgcm/MITgcm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
cost_vector.F
91 lines (73 loc) · 2.56 KB
/
cost_vector.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
C $Header: /u/gcmpack/MITgcm/pkg/cost/cost_vector.F,v 1.5 2012/08/10 19:36:02 jmc Exp $
C $Name: $
#include "COST_OPTIONS.h"
subroutine cost_vector( myThid )
C /==========================================================\
C | subroutine cost_vector |
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_VECTOR
C ========= Local variables =========================
integer isecbeg , isecend , jsec
integer kmaxdepth
integer i, j, k
integer ig, jg
integer bi, bj
_RL locfc
_RL vVel_bar(Nr), theta_bar(Nr), count(Nr)
_RL petawatt
_RL sum
parameter( petawatt = 1.e+15 )
C 80W - 0W at 24N
parameter( isecbeg = 70, isecend = 90, jsec = 27 )
parameter ( kmaxdepth = 15 )
C 80W - 0W at 48N
C parameter( isecbeg = 70, isecend = 90, jsec = 33 )
C parameter ( kmaxdepth = 14 )
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
do j=1,sNy
jg = myYGlobalLo-1+(bj-1)*sNy+j
if (jg .eq. jsec) then
do i=1,sNx
ig = myXGlobalLo-1+(bi-1)*sNx+i
if ((ig .ge. isecbeg) .and. (ig .le. isecend)) then
sum = 0. _d 0
do k = 1, kmaxdepth
sum = sum
& + vVel(i,j,k,bi,bj) * maskS(i,j,k,bi,bj)
& * drF(k)
C & * 0.5*(theta(i,j,k,bi,bj)+theta(i,j-1,k,bi,bj))
end do
objf_vector(i,bi,bj) = sum*dxG(i,j,bi,bj)
end if
end do
end if
end do
do i = 1,sNx
print*,' --> objf_vector(i,bi,bj) = ',
& objf_vector(i,bi,bj)
end do
END DO
END DO
#endif
end