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

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

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


Revision 1.3 - (hide annotations) (download)
Sun Jun 18 23:31:35 2006 UTC (17 years, 11 months ago) by jmc
Branch: MAIN
Changes since 1.2: +9 -9 lines
change velocity argument from a 5-indices global array to a local 2-D array
 (but still not used)

1 jmc 1.3 C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_dst2u1_adv_r.F,v 1.2 2006/06/12 16:06:28 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "GAD_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: GAD_DST2U1_ADV_R
8    
9     C !INTERFACE: ==========================================================
10 jmc 1.2 SUBROUTINE GAD_DST2U1_ADV_R(
11 jmc 1.1 I bi,bj,k, advectionScheme,
12 jmc 1.3 I deltaTloc, rTrans, wFld,
13 jmc 1.1 I tracer,
14     O wT,
15     I myThid )
16    
17     C !DESCRIPTION:
18     C Calculates the area integrated vertical flux due to advection
19 jmc 1.2 C of a tracer using second-order Direct Space and Time (DST-2)
20 jmc 1.1 C interpolation (=Lax-Wendroff) or simple 1rst order upwind scheme.
21    
22     C !USES: ===============================================================
23     IMPLICIT NONE
24     #include "SIZE.h"
25     #include "GRID.h"
26     #include "GAD.h"
27    
28     C !INPUT PARAMETERS: ===================================================
29     C bi,bj :: tile indices
30     C k :: vertical level
31     C advectionScheme :: advection scheme to use: either 2nd Order DST
32     C or 1rst Order Upwind
33     C deltaTloc :: local time-step (s)
34     C rTrans :: vertical volume transport
35 jmc 1.3 C wFld :: vertical flow
36 jmc 1.1 C tracer :: tracer field
37     C myThid :: thread number
38     INTEGER bi,bj,k
39     INTEGER advectionScheme
40     _RL deltaTloc
41     _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42 jmc 1.3 _RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43 jmc 1.1 _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
44     INTEGER myThid
45    
46     C !OUTPUT PARAMETERS: ==================================================
47     C wT :: vertical advective flux
48     _RL wT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
49    
50     C !LOCAL VARIABLES: ====================================================
51     C i,j :: loop indices
52     C km1 :: =max( k-1 , 1 )
53     C rLimit :: centered (vs upwind) fraction
54 jmc 1.3 C wLoc :: velocity, vertical component
55 jmc 1.1 C wCFL :: Courant-Friedrich-Levy number
56     INTEGER i,j,km1
57 jmc 1.3 _RL wLoc, wCFL, rLimit, wAbs
58 jmc 1.1 CEOP
59    
60     rLimit = 0. _d 0
61     IF ( advectionScheme.EQ.ENUM_DST2 ) rLimit = 1. _d 0
62    
63     km1=MAX(1,k-1)
64    
65     IF ( k.LE.1 .OR. k.GT.Nr) THEN
66     DO j=1-Oly,sNy+Oly
67     DO i=1-Olx,sNx+Olx
68     wT(i,j) = 0.
69     ENDDO
70     ENDDO
71     ELSE
72     DO j=1-Oly,sNy+Oly
73     DO i=1-Olx,sNx+Olx
74    
75 jmc 1.3 c wLoc = wFld(i,j)
76     wLoc = rTrans(i,j)*recip_rA(i,j,bi,bj)
77     wCFL = ABS(wLoc*deltaTloc*recip_drC(k))
78 jmc 1.1
79 jmc 1.2 c wT(i,j) = maskC(i,j,km1,bi,bj)*(
80     c & rTrans(i,j)*(tracer(i,j,km1)+tracer(i,j,k))*0.5 _d 0
81     c & + ( 1. _d 0 - rLimit*(1. _d 0 - wCFL) )*ABS(rTrans(i,j))
82     c & *(tracer(i,j,km1)-tracer(i,j,k))*0.5 _d 0*rkSign
83     c & )
84     C-- above formulation produces large truncation error when:
85     C 1rst.O upWind and w*rkSign > 0 & |tracer(k-1)| << |tracer(k)|
86     C or w*rkSign < 0 & |tracer(k-1)| >> |tracer(k)|
87     C-- change to a more robust expression:
88     wAbs = ABS(rTrans(i,j))*rkSign
89     & *( 1. _d 0 - rLimit*(1. _d 0 - wCFL) )
90     wT(i,j) = maskC(i,j,km1,bi,bj)*(
91     & ( rTrans(i,j)+wAbs )* 0.5 _d 0 * tracer(i,j,km1)
92     & + ( rTrans(i,j)-wAbs )* 0.5 _d 0 * tracer(i,j,k)
93 jmc 1.1 & )
94     ENDDO
95     ENDDO
96     ENDIF
97    
98     RETURN
99     END

  ViewVC Help
Powered by ViewVC 1.1.22