-
Notifications
You must be signed in to change notification settings - Fork 237
/
gad_biharm_r.F
103 lines (90 loc) · 3.09 KB
/
gad_biharm_r.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
92
93
94
95
96
97
98
99
100
101
102
103
#include "GAD_OPTIONS.h"
CBOP
C !ROUTINE: GAD_BIHARM_R
C !INTERFACE: ==========================================================
SUBROUTINE GAD_BIHARM_R(
I bi, bj, k,
I maskUp, diffKr4, tracer,
U d4f,
I myThid )
C !DESCRIPTION:
C Calculates the vertical flux due to bi-harmonic diffusion of a tracer.
C \begin{equation*}
C F^r_{diff} = \kappa_4 \partial_r ( \partial^2 / \partial_r^2 \theta )
C \end{equation*}
C !USES: ===============================================================
IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
C !INPUT PARAMETERS: ===================================================
C bi,bj :: tile indices
C k :: vertical level
C maskUp :: 2-D array for mask at W points
C diffKr4 :: vertical bi-harmonic diffusivity
C tracer :: tracer field
C myThid :: my thread Id number
INTEGER bi, bj, k
_RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL diffKr4(Nr)
_RL tracer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
INTEGER myThid
C !OUTPUT PARAMETERS: ==================================================
C dfy :: vertical diffusive flux
_RL d4f (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C !LOCAL VARIABLES: ====================================================
C i,j,n :: loop indices
C del2T :: vertical 2nd derivative of tracer
C gradR :: vertical 1rst derivative of tracer (*rkSign)
INTEGER i, j, n
INTEGER kl, km, kp
_RL del2T(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:2)
_RL gradR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:3)
_RL tmpFac
CEOP
IF ( k.GE.2 ) THEN
C Calculate vertical gradient @ interface k-1, k & k+1 (n=1:3)
DO n=1,3
km = k+n-3
kl = k+n-2
IF ( km.LT.1 .OR. kl.GT.Nr ) THEN
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
gradR(i,j,n) = 0.
ENDDO
ENDDO
ELSE
tmpFac = recip_drC(kl)*deepFac2F(kl)*rhoFacF(kl)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
gradR(i,j,n) = ( tracer(i,j,kl)-tracer(i,j,km) )
& *tmpFac*maskC(i,j,kl,bi,bj)*maskC(i,j,km,bi,bj)
ENDDO
ENDDO
ENDIF
ENDDO
C Calculate vertical 2nd derivative @ level k-1 & k (n=1:2)
DO n=1,2
kl = k+n-2
kp = k+n-1
tmpFac = recip_drF(kl)*recip_deepFac2C(kl)*recip_rhoFacC(kl)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
del2T(i,j,n) = ( gradR(i,j,n+1)-gradR(i,j,n) )
& *_recip_hFacC(i,j,kl,bi,bj)
ENDDO
ENDDO
ENDDO
C Add bi-harmonic flux (gradient of 2nd derivative)
tmpFac = rkSign*recip_drC(k)*deepFac2F(k)*rhoFacF(k)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
d4f(i,j) = d4f(i,j)
& + diffKr4(k)*( del2T(i,j,2)-del2T(i,j,1) )
& *tmpFac*_rA(i,j,bi,bj)*maskUp(i,j)
ENDDO
ENDDO
ENDIF
RETURN
END