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

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

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

Revision 1.3 - (hide annotations) (download)
Wed Jun 7 01:55:14 2006 UTC (17 years, 11 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint58l_post, mitgcm_mapl_00, checkpoint58u_post, checkpoint58w_post, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint60, checkpoint61, checkpoint62, checkpoint63, checkpoint58r_post, checkpoint58n_post, checkpoint58x_post, checkpoint58t_post, checkpoint58h_post, checkpoint58q_post, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint58j_post, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint59j, checkpoint59, checkpoint58i_post, checkpoint58o_post, checkpoint62c, checkpoint62b, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint58y_post, checkpoint58k_post, checkpoint58v_post, checkpoint58s_post, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint58p_post, checkpoint61a, checkpoint61n, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61q, checkpoint61z, checkpoint61x, checkpoint61y, checkpoint58m_post
Changes since 1.2: +9 -9 lines 
Modifications for bottom topography control
o replace hFacC by _hFacC at various places
o replace ALLOW_HFACC_CONTROL by ALLOW_DEPTH_CONTROL
o add non-self-adjoint cg2d_nsa
o update autodiff support routines
o re-initialise hfac after ctrl_depth_ini
o works for 5x5 box, doesnt work for global_ocean.90x40x15
1 heimbach 1.3 C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_dst3fl_impl_r.F,v 1.2 2005/11/04 01:31:36 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "GAD_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: GAD_DST3FL_IMPL_R
8     C !INTERFACE:
9     SUBROUTINE GAD_DST3FL_IMPL_R(
10     I bi,bj,k, iMin,iMax,jMin,jMax,
11     I deltaTarg, rTrans, tFld,
12     O a5d, b5d, c5d, d5d, e5d,
13     I myThid )
14    
15     C !DESCRIPTION:
16    
17     C Compute matrix element to solve vertical advection implicitly
18     C using 3rd order Direct Space and Time (DST) advection scheme
19     C with Flux-Limiter.
20     C Method:
21     C contribution of vertical transport at interface k is added
22     C to matrix lines k and k-1
23    
24     C !USES:
25     IMPLICIT NONE
26    
27     C == Global variables ===
28     #include "SIZE.h"
29     #include "GRID.h"
30     #include "EEPARAMS.h"
31     #include "PARAMS.h"
32     #include "GAD.h"
33    
34     C !INPUT/OUTPUT PARAMETERS:
35     C == Routine Arguments ==
36     C bi,bj :: tile indices
37     C k :: vertical level
38     C iMin,iMax :: computation domain
39     C jMin,jMax :: computation domain
40     C deltaTarg :: time step
41     C rTrans :: vertical volume transport
42     C tFld :: tracer field
43     C a5d :: 2nd lower diag of pentadiagonal matrix
44     C b5d :: 1rst lower diag of pentadiagonal matrix
45     C c5d :: main diag of pentadiagonal matrix
46     C d5d :: 1rst upper diag of pentadiagonal matrix
47     C e5d :: 2nd upper diag of pentadiagonal matrix
48     C myThid :: thread number
49     INTEGER bi,bj,k
50     INTEGER iMin,iMax,jMin,jMax
51     _RL deltaTarg(Nr)
52     _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
53     _RL tFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
54     _RL a5d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
55     _RL b5d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
56     _RL c5d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
57     _RL d5d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
58     _RL e5d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
59     INTEGER myThid
60    
61     C == Local Variables ==
62     C i,j :: loop indices
63     C kp1 :: =min( k+1 , Nr )
64     C km2 :: =max( k-2 , 1 )
65     C wCFL :: Courant-Friedrich-Levy number
66     C lowFac :: low order term factor
67     C highFac :: high order term factor
68     C rCenter :: centered contribution
69     C rUpwind :: upwind contribution
70     C rC4km, rC4kp :: high order contributions
71     INTEGER i,j,kp1,km2
72     _RL wCFL, rCenter, rUpwind
73     _RL lowFac (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74     _RL highFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75     _RL rC4km, rC4kp
76     _RL mskM, mskP, maskM2, maskP1
77     _RL Rj, Rjh, cL1, cH3, cM2, th1, th2
78     _RL deltaTcfl
79     CEOP
80    
81     C-- process interior interface only:
82     IF ( k.GT.1 .AND. k.LE.Nr ) THEN
83    
84     km2=MAX(1,k-2)
85     kp1=MIN(Nr,k+1)
86     maskP1 = 1. _d 0
87     maskM2 = 1. _d 0
88     IF ( k.LE.2 ) maskM2 = 0. _d 0
89     IF ( k.GE.Nr) maskP1 = 0. _d 0
90    
91     C-- Compute the low-order term & high-order term fractions :
92     deltaTcfl = deltaTarg(k)
93     C DST-3 Flux-Limiter Advection Scheme:
94     C- Limiter: Psi=max(0,min(1,cL1+theta*cH1,theta*(1-cfl)/cfl) )
95     C with theta=Rjh/Rj ;
96     C is linearize arround the current value of theta(tFld) & cfl:
97     C lowFac & highFac are set such as Psi*Rj = lowFac*Rj + highFac*Rjh
98     DO j=jMin,jMax
99     DO i=iMin,iMax
100     wCFL = deltaTcfl*ABS(rTrans(i,j))
101     & *recip_rA(i,j,bi,bj)*recip_drC(k)
102     cL1 = (2. _d 0 -wCFL)*(1. _d 0 -wCFL)*oneSixth
103     cH3 = (1. _d 0 -wCFL*wCFL)*oneSixth
104     c cM2 = (1. _d 0 - wCFL)/( wCFL +1. _d -20)
105     cM2 = (1. _d 0 + wCFL)/( wCFL +1. _d -20)
106    
107     Rj =(tFld(i,j,k) -tFld(i,j,k-1))
108     IF ( rTrans(i,j).GT.0. _d 0 ) THEN
109     Rjh = (tFld(i,j,k-1)-tFld(i,j,km2))*maskC(i,j,km2,bi,bj)
110     ELSE
111     Rjh = (tFld(i,j,kp1)-tFld(i,j,k) )*maskC(i,j,kp1,bi,bj)
112     ENDIF
113     IF ( Rj*Rjh.LE.0. _d 0 ) THEN
114     C- 1rst case: theta < 0 (Rj & Rjh opposite sign) => Psi = 0
115     lowFac(i,j) = 0. _d 0
116     highFac(i,j)= 0. _d 0
117     ELSE
118     Rj = ABS(Rj)
119     Rjh = ABS(Rjh)
120     th1 = cL1*Rj+cH3*Rjh
121     th2 = cM2*Rjh
122     IF ( th1.LE.th2 .AND. th1.LE.Rj ) THEN
123     C- 2nd case: cL1+theta*cH3 = min of the three = Psi
124     lowFac(i,j) = cL1
125     highFac(i,j)= cH3
126     ELSEIF ( th2.LT.th1 .AND. th2.LE.Rj ) THEN
127     C- 3rd case: theta*cM2 = min of the three = Psi
128     lowFac(i,j) = 0. _d 0
129     highFac(i,j)= cM2
130     ELSE
131     C- 4th case (Rj < th1 & Rj < th2) : 1 = min of the three = Psi
132     lowFac(i,j) = 1. _d 0
133     highFac(i,j)= 0. _d 0
134     ENDIF
135     ENDIF
136     ENDDO
137     ENDDO
138    
139     C-- Add centered & upwind contributions
140     DO j=jMin,jMax
141     DO i=iMin,iMax
142     rCenter= 0.5 _d 0 *rTrans(i,j)*recip_rA(i,j,bi,bj)*rkSign
143     mskM = maskC(i,j,km2,bi,bj)*maskM2
144     mskP = maskC(i,j,kp1,bi,bj)*maskP1
145     rUpwind= (0.5 _d 0 -lowFac(i,j))*ABS(rCenter)*2. _d 0
146     rC4km = highFac(i,j)*(rCenter+ABS(rCenter))*mskM
147     rC4kp = highFac(i,j)*(rCenter-ABS(rCenter))*mskP
148    
149     a5d(i,j,k) = a5d(i,j,k)
150     & + rC4km
151     & *deltaTarg(k)
152 heimbach 1.3 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
153 jmc 1.1 b5d(i,j,k) = b5d(i,j,k)
154     & - ( (rCenter+rUpwind) + rC4km )
155     & *deltaTarg(k)
156 heimbach 1.3 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
157 jmc 1.1 c5d(i,j,k) = c5d(i,j,k)
158     & - ( (rCenter-rUpwind) + rC4kp )
159     & *deltaTarg(k)
160 heimbach 1.3 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
161 jmc 1.1 d5d(i,j,k) = d5d(i,j,k)
162     & + rC4kp
163     & *deltaTarg(k)
164 heimbach 1.3 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
165 jmc 1.1 b5d(i,j,k-1) = b5d(i,j,k-1)
166     & - rC4km
167     & *deltaTarg(k-1)
168 heimbach 1.3 & *_recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
169 jmc 1.1 c5d(i,j,k-1) = c5d(i,j,k-1)
170     & + ( (rCenter+rUpwind) + rC4km )
171     & *deltaTarg(k-1)
172 heimbach 1.3 & *_recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
173 jmc 1.1 d5d(i,j,k-1) = d5d(i,j,k-1)
174     & + ( (rCenter-rUpwind) + rC4kp )
175     & *deltaTarg(k-1)
176 heimbach 1.3 & *_recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
177 jmc 1.1 e5d(i,j,k-1) = e5d(i,j,k-1)
178     & - rC4kp
179     & *deltaTarg(k-1)
180 heimbach 1.3 & *_recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
181 jmc 1.1 ENDDO
182     ENDDO
183    
184     C-- process interior interface only: end
185     ENDIF
186    
187     RETURN
188     END
  ViewVC Help
Powered by ViewVC 1.1.22