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

Annotation of /MITgcm/verification/matrix_example/code/gad_dst3_adv_y.F

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


Revision 1.1 - (hide annotations) (download)
Mon Apr 18 20:17:29 2005 UTC (19 years, 2 months ago) by spk
Branch: MAIN
CVS Tags: checkpint57u_post, checkpoint57t_post, checkpoint57o_post, checkpoint57h_post, checkpoint57r_post, checkpoint57s_post, checkpoint57l_post, checkpoint57v_post, checkpoint57q_post, checkpoint57j_post, checkpoint57k_post, checkpoint57h_done, checkpoint57i_post, checkpoint57m_post, checkpoint57p_post, checkpoint57h_pre
Changes to model code to implement pkg/matrix

1 spk 1.1 C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_dst3_adv_y.F,v 1.5 2005/04/03 16:05:34 heimbach Exp $
2     C $Name: $
3    
4     #include "GAD_OPTIONS.h"
5    
6     SUBROUTINE GAD_DST3_ADV_Y(
7     I bi,bj,k,deltaTloc,
8     I vTrans, vVel,
9     I maskLocS, tracer,
10     O vT,
11     I myThid )
12     C /==========================================================\
13     C | SUBROUTINE GAD_DST3_ADV_Y |
14     C | o Compute Meridional 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 vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
30     _RL vVel(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
31     _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
32     _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
33     _RL vT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
34     INTEGER myThid
35    
36     C == Local variables ==
37     C vFld :: velocity [m/s], meridional component
38     INTEGER i,j
39     _RL Rjm,Rj,Rjp,cfl,d0,d1
40     _RL psiP,psiM,thetaP,thetaM
41     _RL vFld
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 i=1-Olx,sNx+Olx
52     vT(i,1-Oly)=0.
53     vT(i,2-Oly)=0.
54     vT(i,sNy+Oly)=0.
55     ENDDO
56     DO j=1-Oly+2,sNy+Oly-1
57     DO i=1-Olx,sNx+Olx
58     Rjp=(tracer(i,j+1)-tracer(i, j ))*maskLocS(i,j+1)
59     Rj =(tracer(i, j )-tracer(i,j-1))*maskLocS(i, j )
60     Rjm=(tracer(i,j-1)-tracer(i,j-2))*maskLocS(i,j-1)
61    
62     c vFld = vVel(i,j,k,bi,bj)
63     vFld = vTrans(i,j)*recip_dxG(i,j,bi,bj)
64     & *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
65     cfl=abs(vFld*deltaTloc*recip_dyC(i,j,bi,bj))
66     d0=(2.-cfl)*(1.-cfl)*oneSixth
67     d1=(1.-cfl*cfl)*oneSixth
68     #ifdef ALLOW_MATRIX
69     IF (.NOT.useMATRIX) THEN
70     #endif /* ALLOW_MATRIX */
71     IF ( ABS(Rj).LT.smallNo .OR.
72     & ABS(Rjm).LT.smallNo ) THEN
73     thetaP=0.
74     psiP=0.
75     ELSE
76     thetaP=(Rjm+smallNo)/(smallNo+Rj)
77     psiP=d0+d1*thetaP
78     ENDIF
79     IF ( ABS(Rj).LT.smallNo .OR.
80     & ABS(Rjp).LT.smallNo ) THEN
81     thetaM=0.
82     psiM=0.
83     ELSE
84     thetaM=(Rjp+smallNo)/(smallNo+Rj)
85     psiM=d0+d1*thetaM
86     ENDIF
87     vT(i,j)=
88     & 0.5*(vTrans(i,j)+abs(vTrans(i,j)))
89     & *( Tracer(i,j-1) + psiP*Rj )
90     & +0.5*(vTrans(i,j)-abs(vTrans(i,j)))
91     & *( Tracer(i, j ) - psiM*Rj )
92     #ifdef ALLOW_MATRIX
93     ELSE
94     vT(i,j)=
95     & 0.5*(vTrans(i,j)+abs(vTrans(i,j)))
96     & *( Tracer(i,j-1) + (d0*Rj+d1*Rjm) )
97     & +0.5*(vTrans(i,j)-abs(vTrans(i,j)))
98     & *( Tracer(i, j ) - (d0*Rj+d1*Rjp) )
99     ENDIF
100     #endif /* ALLOW_MATRIX */
101    
102     ENDDO
103     ENDDO
104    
105     RETURN
106     END

  ViewVC Help
Powered by ViewVC 1.1.22