/[MITgcm]/MITgcm/pkg/mom_common/mom_calc_ke.F
ViewVC logotype

Contents of /MITgcm/pkg/mom_common/mom_calc_ke.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.5 - (show annotations) (download)
Sun Mar 4 23:13:07 2012 UTC (12 years, 2 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63k, checkpoint64
Changes since 1.4: +22 -14 lines 
Add full array initialisation of KE (ALLOW_AUTODIFF_TAMC):
Normally, the S/R which calls MOM_CALC_KE initialises KE but in some
recomputation part this initialisation is gone, resulting in Floating
Point Exception (caught by open64 compiler with debug option).
VS: ----------------------------------------------------------------------
1 C $Header: /u/gcmpack/MITgcm/pkg/mom_common/mom_calc_ke.F,v 1.4 2006/06/07 01:55:14 heimbach Exp $
2 C $Name: $
3
4 #include "MOM_COMMON_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: MOM_CALC_KE
8
9 C !INTERFACE: ==========================================================
10 SUBROUTINE MOM_CALC_KE(
11 I bi,bj,k,KEscheme,
12 I uFld, vFld,
13 O KE,
14 I myThid)
15
16 C !DESCRIPTION:
17 C Calculates the Kinetic energy of horizontal flow
18 C \begin{equation*}
19 C KE = \frac{1}{2} \left( h_w \overline{u^2}^i + h_s \overline{v^2}^j \right)
20 C \end{equation*}
21
22 C !USES: ===============================================================
23 IMPLICIT NONE
24 #include "SIZE.h"
25 #include "GRID.h"
26
27 C !INPUT PARAMETERS: ===================================================
28 C bi,bj :: tile indices
29 C k :: vertical level
30 C KEscheme :: spacial discretisation scheme for KE
31 C uFld :: zonal flow
32 C vFld :: meridional flow
33 C KE :: Kinetic Energy
34 C myThid :: thread number
35 INTEGER bi,bj,k
36 _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
37 _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
38 INTEGER KEscheme
39 INTEGER myThid
40
41 C !OUTPUT PARAMETERS: ==================================================
42 C KE :: Kinetic energy
43 _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
44
45 C !LOCAL VARIABLES: ====================================================
46 C i,j :: loop indices
47 INTEGER i,j
48 CEOP
49
50 C This defn of KE should not ever be used. Just to let you know.
51 C 2 2
52 C 1 / ___I ___J \
53 C KE = --- | U + V |
54 C 2 \ /
55 C
56 #ifdef ALLOW_AUTODIFF_TAMC
57 DO j=1-OLy,sNy+OLy
58 DO i=1-OLx,sNx+OLx
59 KE(i,j) = 0.
60 ENDDO
61 ENDDO
62 #endif
63
64 IF (KEscheme.EQ.-1) THEN
65 DO j=1-OLy,sNy+OLy-1
66 DO i=1-OLx,sNx+OLx-1
67 KE(i,j) = 0.125*(
68 & ( uFld(i,j)+uFld(i+1, j ) )**2
69 & +( vFld(i,j)+vFld( i ,j+1) )**2 )
70 ENDDO
71 ENDDO
72
73 ELSEIF (KEscheme.EQ.0) THEN
74 C This defn of KE should be used for the vector invariant equations.
75 C _____I _____J
76 C 1 / 2 2 \
77 C KE = --- | U + V |
78 C 2 \ /
79 C
80 DO j=1-OLy,sNy+OLy-1
81 DO i=1-OLx,sNx+OLx-1
82 KE(i,j) = 0.25*(
83 & ( uFld( i , j )*uFld( i , j )
84 & +uFld(i+1, j )*uFld(i+1, j ) )
85 & + ( vFld( i , j )*vFld( i , j )
86 & +vFld( i ,j+1)*vFld( i ,j+1) )
87 & )
88 ENDDO
89 ENDDO
90
91 ELSEIF (KEscheme.EQ.1) THEN
92 C As above but including the area
93 DO j=1-OLy,sNy+OLy-1
94 DO i=1-OLx,sNx+OLx-1
95 KE(i,j) = 0.25*(
96 & ( uFld(i, j )*uFld(i, j )*rAw(i ,j, bi,bj)
97 & +uFld(i+1,j)*uFld(i+1,j)*rAw(i+1,j,bi,bj) )
98 & + ( vFld(i, j )*vFld(i, j )*rAs(i ,j, bi,bj)
99 & +vFld(i,j+1)*vFld(i,j+1)*rAs(i,j+1,bi,bj) )
100 & )*recip_rA(i,j,bi,bj)
101 ENDDO
102 ENDDO
103
104 ELSEIF (KEscheme.EQ.2) THEN
105 C As KEscheme=0 but including the lopping factors and should be used
106 C for the conservative form of the momentum equations.
107 DO j=1-OLy,sNy+OLy-1
108 DO i=1-OLx,sNx+OLx-1
109 KE(i,j) = 0.25*(
110 & ( uFld( i , j )*uFld( i , j )*_hFacW(i,j,k,bi,bj)
111 & +uFld(i+1, j )*uFld(i+1, j )*_hFacW(i+1,j,k,bi,bj) )
112 & + ( vFld( i , j )*vFld( i , j )*_hFacS(i,j,k,bi,bj)
113 & +vFld( i ,j+1)*vFld( i ,j+1)*_hFacS(i,j+1,k,bi,bj) )
114 & )*_recip_hFacC(i,j,k,bi,bj)
115 ENDDO
116 ENDDO
117
118 ELSEIF (KEscheme.EQ.3) THEN
119 C As above but including the area
120 DO j=1-OLy,sNy+OLy-1
121 DO i=1-OLx,sNx+OLx-1
122 KE(i,j) = 0.25*(
123 & (
124 & uFld(i, j )*uFld(i, j )
125 & *_hFacW(i ,j, k,bi,bj)*rAw(i ,j, bi,bj)
126 & +uFld(i+1,j)*uFld(i+1,j)
127 & *_hFacW(i+1,j,k,bi,bj)*rAw(i+1,j,bi,bj)
128 & )
129 & + (
130 & vFld(i, j )*vFld(i, j )
131 & *_hFacS(i, j, k,bi,bj)*rAs(i ,j, bi,bj)
132 & +vFld(i,j+1)*vFld(i,j+1)
133 & *_hFacS(i,j+1,k,bi,bj)*rAs(i,j+1,bi,bj)
134 & ) )*_recip_hFacC(i,j,k,bi,bj)
135 & * recip_rA(i,j,bi,bj)
136 ENDDO
137 ENDDO
138
139 ELSE
140 STOP 'S/R MOM_CALC_KE: We should never reach this point!'
141 ENDIF
142
143 RETURN
144 END
  ViewVC Help
Powered by ViewVC 1.1.22