/[MITgcm]/MITgcm/pkg/generic_advdiff/gad_dst3_adv_x.F
ViewVC logotype

Annotation of /MITgcm/pkg/generic_advdiff/gad_dst3_adv_x.F

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


Revision 1.5 - (hide annotations) (download)
Sun Apr 3 16:05:34 2005 UTC (19 years, 2 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint57m_post, checkpoint57g_pre, checkpoint57s_post, checkpoint57g_post, checkpoint57r_post, checkpoint57i_post, checkpoint57n_post, checkpoint57l_post, checkpoint57t_post, checkpoint57v_post, checkpoint57f_post, checkpoint57h_pre, checkpoint57h_post, checkpoint57p_post, checkpint57u_post, checkpoint57q_post, checkpoint57h_done, checkpoint57j_post, checkpoint57f_pre, checkpoint57o_post, checkpoint57k_post
Changes since 1.4: +30 -14 lines
Modifications to make DST3(=30) adjointable
(seems to work at 1/4 deg. for at least 1200 timesteps).

1 heimbach 1.5 C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_dst3_adv_x.F,v 1.4 2004/09/24 16:53:45 jmc Exp $
2 adcroft 1.2 C $Name: $
3 adcroft 1.1
4     #include "GAD_OPTIONS.h"
5    
6     SUBROUTINE GAD_DST3_ADV_X(
7 heimbach 1.5 I bi,bj,k,deltaTloc,
8 adcroft 1.1 I uTrans, uVel,
9 jmc 1.4 I maskLocW, tracer,
10 adcroft 1.1 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 heimbach 1.5 #include "EEPARAMS.h"
23     #include "PARAMS.h"
24 adcroft 1.1 #include "GAD.h"
25    
26     C == Routine arguments ==
27     INTEGER bi,bj,k
28 heimbach 1.5 _RL deltaTloc
29 adcroft 1.1 _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 jmc 1.4 _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
32 adcroft 1.1 _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 jmc 1.3 C uFld :: velocity [m/s], zonal component
38 adcroft 1.1 INTEGER i,j
39     _RL Rjm,Rj,Rjp,cfl,d0,d1
40 adcroft 1.2 _RL psiP,psiM,thetaP,thetaM
41 jmc 1.3 _RL uFld
42 heimbach 1.5 _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 adcroft 1.1
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 jmc 1.4 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 adcroft 1.1
60 jmc 1.3 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 heimbach 1.5 cfl=abs(uFld*deltaTloc*recip_dxC(i,j,bi,bj))
64 adcroft 1.2 d0=(2.-cfl)*(1.-cfl)*oneSixth
65     d1=(1.-cfl*cfl)*oneSixth
66 heimbach 1.5 IF ( ABS(Rj).LT.smallNo .OR.
67     & ABS(Rjm).LT.smallNo ) THEN
68     thetaP=0.
69     psiP=0.
70     ELSE
71     thetaP=(Rjm+smallNo)/(smallNo+Rj)
72     psiP=d0+d1*thetaP
73     ENDIF
74     IF ( ABS(Rj).LT.smallNo .OR.
75     & ABS(Rjp).LT.smallNo ) THEN
76     thetaM=0.
77     psiM=0.
78     ELSE
79     thetaM=(Rjp+smallNo)/(smallNo+Rj)
80     psiM=d0+d1*thetaM
81     ENDIF
82 adcroft 1.1 uT(i,j)=
83     & 0.5*(uTrans(i,j)+abs(uTrans(i,j)))
84 adcroft 1.2 & *( Tracer(i-1,j) + psiP*Rj )
85 adcroft 1.1 & +0.5*(uTrans(i,j)-abs(uTrans(i,j)))
86 adcroft 1.2 & *( Tracer( i ,j) - psiM*Rj )
87 adcroft 1.1
88     ENDDO
89     ENDDO
90    
91     RETURN
92     END

  ViewVC Help
Powered by ViewVC 1.1.22