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

Contents of /MITgcm/pkg/generic_advdiff/gad_fluxlimit_adv_r.F

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


Revision 1.6 - (show annotations) (download)
Thu Jan 17 16:48:59 2002 UTC (22 years, 4 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint44e_post, checkpoint44f_post, chkpt44d_post, checkpoint44e_pre, chkpt44a_post, chkpt44c_pre, checkpoint44g_post, release1_final_v1, checkpoint44b_post, chkpt44a_pre, checkpoint44b_pre, checkpoint44, chkpt44c_post, checkpoint44f_pre
Branch point for: release1_final
Changes since 1.5: +3 -3 lines
Bug fix:

The mask arrays had the wrong tile indices
 o only affected multi-tile per thread calculations

1 C $Header: /u/gcmpack/models/MITgcmUV/pkg/generic_advdiff/gad_fluxlimit_adv_r.F,v 1.5 2001/09/21 13:11:43 adcroft Exp $
2 C $Name: $
3
4 #include "GAD_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: GAD_FLUXLIMIT_ADV_R
8
9 C !INTERFACE: ==========================================================
10 SUBROUTINE GAD_FLUXLIMIT_ADV_R(
11 I bi_arg,bj_arg,k,dTarg,
12 I rTrans, wVel,
13 I tracer,
14 O wT,
15 I myThid )
16
17 C !DESCRIPTION:
18 C Calculates the area integrated vertical flux due to advection of a tracer
19 C using second-order interpolation with a flux limiter:
20 C \begin{equation*}
21 C F^x_{adv} = W \overline{ \theta }^k
22 C - \frac{1}{2} \left(
23 C [ 1 - \psi(C_r) ] |W|
24 C + W \frac{w \Delta t}{\Delta r_c} \psi(C_r)
25 C \right) \delta_k \theta
26 C \end{equation*}
27 C where the $\psi(C_r)$ is the limiter function and $C_r$ is
28 C the slope ratio.
29
30 C !USES: ===============================================================
31 IMPLICIT NONE
32 #include "SIZE.h"
33 #include "GRID.h"
34 #include "EEPARAMS.h"
35 #include "PARAMS.h"
36
37 C !INPUT PARAMETERS: ===================================================
38 C bi_arg,bj_arg :: tile indices
39 C k :: vertical level
40 C rTrans :: vertical volume transport
41 C wVel :: vertical flow
42 C tracer :: tracer field
43 C myThid :: thread number
44 INTEGER bi_arg,bj_arg,k
45 _RL dTarg
46 _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47 _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
48 _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
49 INTEGER myThid
50
51 C !OUTPUT PARAMETERS: ==================================================
52 C wT :: vertical advective flux
53 _RL wT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54
55 C !LOCAL VARIABLES: ====================================================
56 C i,j :: loop indices
57 C kp1 :: =min( k+1 , Nr )
58 C km1 :: =max( k-1 , 1 )
59 C km2 :: =max( k-2 , 1 )
60 C bi,bj :: tile indices or (1,1) depending on use
61 C Cr :: slope ratio
62 C Rjm,Rj,Rjp :: differences at i-1,i,i+1
63 INTEGER i,j,kp1,km1,km2,bi,bj
64 _RL Cr,Rjm,Rj,Rjp
65 C Statement function provides Limiter(Cr)
66 #include "GAD_FLUX_LIMITER.h"
67 CEOP
68
69 IF (.NOT. multiDimAdvection) THEN
70 C If using the standard time-stepping/advection schemes (ie. AB-II)
71 C then the data-structures are all global arrays
72 bi=bi_arg
73 bj=bj_arg
74 ELSE
75 C otherwise if using the multi-dimensional advection schemes
76 C then the data-structures are all local arrays except
77 C for maskC(...) and wVel(...)
78 bi=1
79 bj=1
80 ENDIF
81
82 km2=MAX(1,k-2)
83 km1=MAX(1,k-1)
84 kp1=MIN(Nr,k+1)
85
86 IF ( k.GT.Nr) THEN
87 DO j=1-Oly,sNy+Oly
88 DO i=1-Olx,sNx+Olx
89 wT(i,j) = 0.
90 ENDDO
91 ENDDO
92 ELSE
93 DO j=1-Oly,sNy+Oly
94 DO i=1-Olx,sNx+Olx
95 Rjp=(tracer(i,j,kp1,bi,bj)-tracer(i,j,k,bi,bj))
96 & *maskC(i,j,kp1,bi_arg,bj_arg)
97 Rj=(tracer(i,j,k,bi,bj)-tracer(i,j,kM1,bi,bj))
98 Rjm=(tracer(i,j,km1,bi,bj)-tracer(i,j,kM2,bi,bj))
99 & *maskC(i,j,km2,bi_arg,bj_arg)
100 IF (Rj.NE.0.) THEN
101 IF (rTrans(i,j).LT.0.) THEN
102 Cr=Rjm/Rj
103 ELSE
104 Cr=Rjp/Rj
105 ENDIF
106 ELSE
107 IF (rTrans(i,j).LT.0.) THEN
108 Cr=Rjm*1.E20
109 ELSE
110 Cr=Rjp*1.E20
111 ENDIF
112 ENDIF
113 Cr=Limiter(Cr)
114 wT(i,j) = maskC(i,j,kM1,bi_arg,bj_arg)*(
115 & rTrans(i,j)*
116 & (Tracer(i,j,k,bi,bj)+Tracer(i,j,kM1,bi,bj))*0.5 _d 0
117 & +(ABS(rTrans(i,j))*(1-Cr)
118 & +rTrans(i,j)*wVel(i,j,k,bi_arg,bj_arg)*dTarg*recip_drC(k)
119 & *Cr
120 & )*Rj*0.5 _d 0 )
121 ENDDO
122 ENDDO
123 ENDIF
124
125 RETURN
126 END

  ViewVC Help
Powered by ViewVC 1.1.22