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

Contents 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 - (show annotations) (download)
Wed Jun 7 01:55:14 2006 UTC (17 years, 10 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 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 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 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
153 b5d(i,j,k) = b5d(i,j,k)
154 & - ( (rCenter+rUpwind) + rC4km )
155 & *deltaTarg(k)
156 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
157 c5d(i,j,k) = c5d(i,j,k)
158 & - ( (rCenter-rUpwind) + rC4kp )
159 & *deltaTarg(k)
160 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
161 d5d(i,j,k) = d5d(i,j,k)
162 & + rC4kp
163 & *deltaTarg(k)
164 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
165 b5d(i,j,k-1) = b5d(i,j,k-1)
166 & - rC4km
167 & *deltaTarg(k-1)
168 & *_recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
169 c5d(i,j,k-1) = c5d(i,j,k-1)
170 & + ( (rCenter+rUpwind) + rC4km )
171 & *deltaTarg(k-1)
172 & *_recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
173 d5d(i,j,k-1) = d5d(i,j,k-1)
174 & + ( (rCenter-rUpwind) + rC4kp )
175 & *deltaTarg(k-1)
176 & *_recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
177 e5d(i,j,k-1) = e5d(i,j,k-1)
178 & - rC4kp
179 & *deltaTarg(k-1)
180 & *_recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
181 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