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

Annotation of /MITgcm/pkg/generic_advdiff/gad_dst3_adv_y.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, 1 month 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_y.F,v 1.4 2004/09/24 16:53:46 jmc Exp $
2 jmc 1.3 C $Name: $
3 adcroft 1.1
4     #include "GAD_OPTIONS.h"
5    
6     SUBROUTINE GAD_DST3_ADV_Y(
7 heimbach 1.5 I bi,bj,k,deltaTloc,
8 adcroft 1.1 I vTrans, vVel,
9 jmc 1.4 I maskLocS, tracer,
10 adcroft 1.1 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 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 vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
30     _RL vVel(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
31 jmc 1.4 _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
32 adcroft 1.1 _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 jmc 1.3 C vFld :: velocity [m/s], meridional 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 vFld
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 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 jmc 1.4 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 adcroft 1.1
62 jmc 1.3 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 heimbach 1.5 cfl=abs(vFld*deltaTloc*recip_dyC(i,j,bi,bj))
66 adcroft 1.2 d0=(2.-cfl)*(1.-cfl)*oneSixth
67     d1=(1.-cfl*cfl)*oneSixth
68 heimbach 1.5 IF ( ABS(Rj).LT.smallNo .OR.
69     & ABS(Rjm).LT.smallNo ) THEN
70     thetaP=0.
71     psiP=0.
72     ELSE
73     thetaP=(Rjm+smallNo)/(smallNo+Rj)
74     psiP=d0+d1*thetaP
75     ENDIF
76     IF ( ABS(Rj).LT.smallNo .OR.
77     & ABS(Rjp).LT.smallNo ) THEN
78     thetaM=0.
79     psiM=0.
80     ELSE
81     thetaM=(Rjp+smallNo)/(smallNo+Rj)
82     psiM=d0+d1*thetaM
83     ENDIF
84 adcroft 1.1 vT(i,j)=
85     & 0.5*(vTrans(i,j)+abs(vTrans(i,j)))
86 adcroft 1.2 & *( Tracer(i,j-1) + psiP*Rj )
87 adcroft 1.1 & +0.5*(vTrans(i,j)-abs(vTrans(i,j)))
88 adcroft 1.2 & *( Tracer(i, j ) - psiM*Rj )
89 adcroft 1.1
90     ENDDO
91     ENDDO
92    
93     RETURN
94     END

  ViewVC Help
Powered by ViewVC 1.1.22