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

Annotation of /MITgcm/pkg/generic_advdiff/gad_fluxlimit_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, 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.7: +5 -5 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_fluxlimit_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     C !ROUTINE: GAD_FLUXLIMIT_IMPL_R
8     C !INTERFACE:
9     SUBROUTINE GAD_FLUXLIMIT_IMPL_R(
10     I bi,bj,k, iMin,iMax,jMin,jMax,
11     I deltaTarg, rTrans, tFld,
12     O a3d, b3d, c3d,
13     I myThid )
14    
15 edhill 1.3 C !DESCRIPTION:
16     C Compute matrix element to solve vertical advection implicitly
17 jmc 1.7 C using flux--limiter advection scheme.
18     C Method:
19     C contribution of vertical transport at interface k is added
20     C to matrix lines k and k-1.
21 jmc 1.1
22     C !USES:
23     IMPLICIT NONE
24    
25     C == Global variables ===
26     #include "SIZE.h"
27     #include "GRID.h"
28     #include "EEPARAMS.h"
29     #include "PARAMS.h"
30    
31     C !INPUT/OUTPUT PARAMETERS:
32     C == Routine Arguments ==
33 edhill 1.3 C bi,bj :: tile indices
34     C k :: vertical level
35     C iMin,iMax :: computation domain
36     C jMin,jMax :: computation domain
37     C deltaTarg :: time step
38     C rTrans :: vertical volume transport
39     C tFld :: tracer field
40     C a3d :: lower diagonal of the tridiagonal matrix
41     C b3d :: main diagonal of the tridiagonal matrix
42     C c3d :: upper diagonal of the tridiagonal matrix
43     C myThid :: thread number
44 jmc 1.1 INTEGER bi,bj,k
45     INTEGER iMin,iMax,jMin,jMax
46 jmc 1.5 _RL deltaTarg(Nr)
47 jmc 1.1 _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48     _RL tFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
49     _RL a3d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
50     _RL b3d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
51     _RL c3d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
52     INTEGER myThid
53    
54     C == Local Variables ==
55 edhill 1.3 C i,j :: loop indices
56     C kp1 :: =min( k+1 , Nr )
57     C km1 :: =max( k-1 , 1 )
58     C km2 :: =max( k-2 , 1 )
59     C Cr :: slope ratio
60     C Rjm,Rj,Rjp :: differences at i-1,i,i+1
61     C w_CFL :: Courant-Friedrich-Levy number
62     C upwindFac :: upwind factor
63     C rCenter :: centered contribution
64     C rUpwind :: upwind contribution
65 jmc 1.1 INTEGER i,j,kp1,km1,km2
66     _RL Cr,Rjm,Rj,Rjp, w_CFL
67     _RL upwindFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68     _RL rCenter, rUpwind
69 jmc 1.5 _RL deltaTcfl
70 jmc 1.1
71     C Statement function provides Limiter(Cr)
72     #include "GAD_FLUX_LIMITER.h"
73     CEOP
74    
75     km2=MAX(1,k-2)
76     km1=MAX(1,k-1)
77     kp1=MIN(Nr,k+1)
78    
79 jmc 1.5 C-- process interior interface only:
80     IF ( k.GT.1 .AND. k.LE.Nr ) THEN
81 jmc 1.1
82     C-- Compute the upwind fraction:
83 jmc 1.5 deltaTcfl = deltaTarg(k)
84 jmc 1.1 DO j=jMin,jMax
85     DO i=iMin,iMax
86 jmc 1.5 w_CFL = deltaTcfl*rTrans(i,j)*recip_rA(i,j,bi,bj)*recip_drC(k)
87 jmc 1.1 Rjp=(tFld(i,j,kp1)-tFld(i,j,k) )*maskC(i,j,kp1,bi,bj)
88     Rj =(tFld(i,j,k) -tFld(i,j,km1))
89     Rjm=(tFld(i,j,km1)-tFld(i,j,km2))*maskC(i,j,km2,bi,bj)
90    
91     IF ( Rj.NE.0. _d 0) THEN
92     IF (rTrans(i,j).LT.0. _d 0) THEN
93     Cr=Rjm/Rj
94     ELSE
95     Cr=Rjp/Rj
96     ENDIF
97     upwindFac(i,j) = 1. _d 0
98 jmc 1.6 & - Limiter(Cr) * ( 1. _d 0 + ABS(w_CFL) )
99     upwindFac(i,j) = MAX( -1. _d 0, upwindFac(i,j) )
100 jmc 1.1 ELSE
101     upwindFac(i,j) = 0. _d 0
102     ENDIF
103     ENDDO
104     ENDDO
105    
106     C-- Add centered & upwind contributions
107     DO j=jMin,jMax
108     DO i=iMin,iMax
109 jmc 1.6 rCenter = 0.5 _d 0 *rTrans(i,j)*recip_rA(i,j,bi,bj)*rkSign
110     rUpwind = ABS(rCenter)*upwindFac(i,j)
111 jmc 1.1 a3d(i,j,k) = a3d(i,j,k)
112 jmc 1.6 & - (rCenter+rUpwind)*deltaTarg(k)
113 heimbach 1.8 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
114 jmc 1.1 b3d(i,j,k) = b3d(i,j,k)
115 jmc 1.6 & - (rCenter-rUpwind)*deltaTarg(k)
116 heimbach 1.8 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
117 jmc 1.1 b3d(i,j,k-1) = b3d(i,j,k-1)
118 jmc 1.6 & + (rCenter+rUpwind)*deltaTarg(k-1)
119 heimbach 1.8 & *_recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
120 jmc 1.1 c3d(i,j,k-1) = c3d(i,j,k-1)
121 jmc 1.6 & + (rCenter-rUpwind)*deltaTarg(k-1)
122 heimbach 1.8 & *_recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
123 jmc 1.1 ENDDO
124     ENDDO
125    
126 jmc 1.5 C-- process interior interface only: end
127     ENDIF
128    
129 jmc 1.1 RETURN
130     END
  ViewVC Help
Powered by ViewVC 1.1.22