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

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

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

Revision 1.8 - (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, 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.7: +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.8 C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_u3c4_impl_r.F,v 1.7 2005/10/22 20:17:44 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "GAD_OPTIONS.h"
5    
6     CBOP
7 jmc 1.6 C !ROUTINE: GAD_U3C4_IMPL_R
8 jmc 1.1 C !INTERFACE:
9 jmc 1.7 SUBROUTINE GAD_U3C4_IMPL_R(
10     I bi,bj,k, iMin,iMax,jMin,jMax,
11 jmc 1.1 I advectionScheme, deltaTarg, rTrans,
12     O a5d, b5d, c5d, d5d, e5d,
13     I myThid )
14    
15 edhill 1.2 C !DESCRIPTION:
16 jmc 1.6 C Compute matrix element to solve vertical advection implicitly
17     C using 3rd order upwind advection scheme,
18 jmc 1.7 C or 3rd order Direct Space and Time advection scheme,
19 jmc 1.6 C or 4th order Centered advection scheme.
20     C Method:
21     C contribution of vertical transport at interface k is added
22     C to matrix lines k and k-1
23 jmc 1.1
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 jmc 1.6 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 advectionScheme :: advection scheme to use
41     C deltaTarg :: time step
42     C rTrans :: vertical volume transport
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 jmc 1.1 INTEGER bi,bj,k
50     INTEGER iMin,iMax,jMin,jMax
51     INTEGER advectionScheme
52 jmc 1.4 _RL deltaTarg(Nr)
53 jmc 1.1 _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
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 jmc 1.6 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 rCenter :: centered contribution
66     C rUpwind :: upwind contribution
67 jmc 1.7 C rC4km, rC4kp :: high order contribution
68     C rHigh :: high order term factor
69 jmc 1.1 LOGICAL flagC4
70     INTEGER i,j,kp1,km2
71 jmc 1.7 _RL wCFL, rCenter, rUpwind
72     _RL rC4km, rC4kp, rHigh
73 jmc 1.1 _RL mskM, mskP, maskM2, maskP1
74 jmc 1.7 _RL deltaTcfl
75 jmc 1.1 CEOP
76    
77 jmc 1.4 C-- process interior interface only:
78     IF ( k.GT.1 .AND. k.LE.Nr ) THEN
79 jmc 1.1
80     km2=MAX(1,k-2)
81     kp1=MIN(Nr,k+1)
82 jmc 1.7 maskP1 = 1. _d 0
83     maskM2 = 1. _d 0
84 jmc 1.1 IF ( k.LE.2 ) maskM2 = 0. _d 0
85     IF ( k.GE.Nr) maskP1 = 0. _d 0
86 jmc 1.7 flagC4 = advectionScheme.EQ.ENUM_CENTERED_4TH
87 jmc 1.1 & .AND. k.GT.2 .AND. k.LT.Nr
88    
89 jmc 1.7 C-- Add centered, upwind and high-order contributions
90     deltaTcfl = deltaTarg(k)
91 jmc 1.1 DO j=jMin,jMax
92     DO i=iMin,iMax
93 jmc 1.5 rCenter= 0.5 _d 0 *rTrans(i,j)*recip_rA(i,j,bi,bj)*rkSign
94 jmc 1.1 mskM = maskC(i,j,km2,bi,bj)*maskM2
95     mskP = maskC(i,j,kp1,bi,bj)*maskP1
96     IF ( flagC4 .AND. mskM*mskP.GT.0. _d 0 ) THEN
97     rUpwind= 0. _d 0
98 jmc 1.7 rC4km = oneSixth*rCenter*mskM
99     rC4kp = oneSixth*rCenter*mskP
100     ELSEIF ( advectionScheme.EQ.ENUM_DST3 ) THEN
101     wCFL = deltaTcfl*ABS(rTrans(i,j))
102     & *recip_rA(i,j,bi,bj)*recip_drC(k)
103     rHigh = (1. _d 0 -wCFL*wCFL)*oneSixth
104     c rUpwind= (2. _d 0*rHigh - wCFL)*ABS(rCenter)
105     rUpwind= (2. _d 0*rHigh )*ABS(rCenter)
106     rC4km = rHigh * (rCenter+ABS(rCenter))*mskM
107     rC4kp = rHigh * (rCenter-ABS(rCenter))*mskP
108 jmc 1.1 ELSE
109 jmc 1.7 rUpwind= 2. _d 0*oneSixth*ABS(rCenter)
110     rC4km = oneSixth*(rCenter+ABS(rCenter))*mskM
111     rC4kp = oneSixth*(rCenter-ABS(rCenter))*mskP
112 jmc 1.1 ENDIF
113     a5d(i,j,k) = a5d(i,j,k)
114 jmc 1.7 & + rC4km
115 jmc 1.4 & *deltaTarg(k)
116 heimbach 1.8 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
117 jmc 1.1 b5d(i,j,k) = b5d(i,j,k)
118 jmc 1.7 & - ( (rCenter+rUpwind) + rC4km )
119 jmc 1.4 & *deltaTarg(k)
120 heimbach 1.8 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
121 jmc 1.1 c5d(i,j,k) = c5d(i,j,k)
122 jmc 1.7 & - ( (rCenter-rUpwind) + rC4kp )
123 jmc 1.4 & *deltaTarg(k)
124 heimbach 1.8 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
125 jmc 1.1 d5d(i,j,k) = d5d(i,j,k)
126 jmc 1.7 & + rC4kp
127 jmc 1.4 & *deltaTarg(k)
128 heimbach 1.8 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
129 jmc 1.1 b5d(i,j,k-1) = b5d(i,j,k-1)
130 jmc 1.7 & - rC4km
131 jmc 1.4 & *deltaTarg(k-1)
132 heimbach 1.8 & *_recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
133 jmc 1.1 c5d(i,j,k-1) = c5d(i,j,k-1)
134 jmc 1.7 & + ( (rCenter+rUpwind) + rC4km )
135 jmc 1.4 & *deltaTarg(k-1)
136 heimbach 1.8 & *_recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
137 jmc 1.1 d5d(i,j,k-1) = d5d(i,j,k-1)
138 jmc 1.7 & + ( (rCenter-rUpwind) + rC4kp )
139 jmc 1.4 & *deltaTarg(k-1)
140 heimbach 1.8 & *_recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
141 jmc 1.1 e5d(i,j,k-1) = e5d(i,j,k-1)
142 jmc 1.7 & - rC4kp
143 jmc 1.4 & *deltaTarg(k-1)
144 heimbach 1.8 & *_recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
145 jmc 1.1 ENDDO
146     ENDDO
147    
148 jmc 1.4 C-- process interior interface only: end
149     ENDIF
150    
151 jmc 1.1 RETURN
152     END
  ViewVC Help
Powered by ViewVC 1.1.22