/[MITgcm]/MITgcm/verification/matrix_example/code/gad_dst3_adv_x.F
ViewVC logotype

Contents of /MITgcm/verification/matrix_example/code/gad_dst3_adv_x.F

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


Revision 1.2 - (show annotations) (download)
Fri Oct 14 22:14:56 2005 UTC (18 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +1 -1 lines
FILE REMOVED
now checked in pkg/generic_advdiff.

1 C $Header: /u/gcmpack/MITgcm/verification/matrix_example/code/gad_dst3_adv_x.F,v 1.1 2005/04/18 20:17:14 spk Exp $
2 C $Name: $
3
4 #include "GAD_OPTIONS.h"
5
6 SUBROUTINE GAD_DST3_ADV_X(
7 I bi,bj,k,deltaTloc,
8 I uTrans, uVel,
9 I maskLocW, tracer,
10 O uT,
11 I myThid )
12 C /==========================================================\
13 C | SUBROUTINE GAD_DST3_ADV_X |
14 C | o Compute Zonal advective Flux of Tracer using |
15 C | 3rd Order DST Sceheme |
16 C |==========================================================|
17 IMPLICIT NONE
18
19 C == GLobal variables ==
20 #include "SIZE.h"
21 #include "GRID.h"
22 #include "EEPARAMS.h"
23 #include "PARAMS.h"
24 #include "GAD.h"
25
26 C == Routine arguments ==
27 INTEGER bi,bj,k
28 _RL deltaTloc
29 _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
30 _RL uVel(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
31 _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
32 _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
33 _RL uT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
34 INTEGER myThid
35
36 C == Local variables ==
37 C uFld :: velocity [m/s], zonal component
38 INTEGER i,j
39 _RL Rjm,Rj,Rjp,cfl,d0,d1
40 _RL psiP,psiM,thetaP,thetaM
41 _RL uFld
42 _RL smallNo
43 _RL Rjjm,Rjjp
44
45 IF (inAdMode) THEN
46 smallNo = 1.0D-20
47 ELSE
48 smallNo = 1.0D-20
49 ENDIF
50
51 DO j=1-Oly,sNy+Oly
52 uT(1-Olx,j)=0.
53 uT(2-Olx,j)=0.
54 uT(sNx+Olx,j)=0.
55 DO i=1-Olx+2,sNx+Olx-1
56 Rjp=(tracer(i+1,j)-tracer( i ,j))*maskLocW(i+1,j)
57 Rj =(tracer( i ,j)-tracer(i-1,j))*maskLocW( i ,j)
58 Rjm=(tracer(i-1,j)-tracer(i-2,j))*maskLocW(i-1,j)
59
60 c uFld = uVel(i,j,k,bi,bj)
61 uFld = uTrans(i,j)*recip_dyG(i,j,bi,bj)
62 & *recip_drF(k)*recip_hFacW(i,j,k,bi,bj)
63 cfl=abs(uFld*deltaTloc*recip_dxC(i,j,bi,bj))
64 d0=(2.-cfl)*(1.-cfl)*oneSixth
65 d1=(1.-cfl*cfl)*oneSixth
66 #ifdef ALLOW_MATRIX
67 IF (.NOT.useMATRIX) THEN
68 #endif /* ALLOW_MATRIX */
69 IF ( ABS(Rj).LT.smallNo .OR.
70 & ABS(Rjm).LT.smallNo ) THEN
71 thetaP=0.
72 psiP=0.
73 ELSE
74 thetaP=(Rjm+smallNo)/(smallNo+Rj)
75 psiP=d0+d1*thetaP
76 ENDIF
77 IF ( ABS(Rj).LT.smallNo .OR.
78 & ABS(Rjp).LT.smallNo ) THEN
79 thetaM=0.
80 psiM=0.
81 ELSE
82 thetaM=(Rjp+smallNo)/(smallNo+Rj)
83 psiM=d0+d1*thetaM
84 ENDIF
85 uT(i,j)=
86 & 0.5*(uTrans(i,j)+abs(uTrans(i,j)))
87 & *( Tracer(i-1,j) + psiP*Rj )
88 & +0.5*(uTrans(i,j)-abs(uTrans(i,j)))
89 & *( Tracer( i ,j) - psiM*Rj )
90 #ifdef ALLOW_MATRIX
91 ELSE
92 uT(i,j)=
93 & 0.5*(uTrans(i,j)+abs(uTrans(i,j)))
94 & *( Tracer(i-1,j) + (d0*Rj+d1*Rjm) )
95 & +0.5*(uTrans(i,j)-abs(uTrans(i,j)))
96 & *( Tracer( i ,j) - (d0*Rj+d1*Rjp) )
97 ENDIF
98 #endif /* ALLOW_MATRIX */
99
100 ENDDO
101 ENDDO
102
103 RETURN
104 END

  ViewVC Help
Powered by ViewVC 1.1.22