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

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

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

revision 1.3 by adcroft, Mon Sep 10 01:22:48 2001 UTC revision 1.11 by jmc, Mon Jun 19 14:40:43 2006 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
5    
6        SUBROUTINE GAD_FLUXLIMIT_ADV_R(  CBOP
7       I           bi_arg,bj_arg,k,dTarg,  C !ROUTINE: GAD_FLUXLIMIT_ADV_R
8       I           rTrans, wVel,  
9    C !INTERFACE: ==========================================================
10          SUBROUTINE GAD_FLUXLIMIT_ADV_R(
11         I           bi,bj,k,dTarg,
12         I           rTrans, wFld,
13       I           tracer,       I           tracer,
14       O           wT,       O           wT,
15       I           myThid )       I           myThid )
 C     /==========================================================\  
 C     | SUBROUTINE GAD_FLUXLIMIT_ADV_R                           |  
 C     | o Compute vertical advective Flux of Tracer using        |  
 C     |   Flux Limiter Scheme                                    |  
 C     |==========================================================|  
       IMPLICIT NONE  
16    
17  C     == GLobal variables ==  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"  #include "SIZE.h"
33  #include "GRID.h"  #include "GRID.h"
34  #include "EEPARAMS.h"  #include "EEPARAMS.h"
35  #include "PARAMS.h"  #include "PARAMS.h"
36    
37  C     == Routine arguments ==  C !INPUT PARAMETERS: ===================================================
38        INTEGER bi_arg,bj_arg,k  C  bi,bj              :: tile indices
39    C  k                  :: vertical level
40    C  rTrans             :: vertical volume transport
41    C  wFld               :: vertical flow
42    C  tracer             :: tracer field
43    C  myThid             :: thread number
44          INTEGER bi,bj,k
45        _RL dTarg        _RL dTarg
46        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47        _RL wVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48        _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL wT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
49        INTEGER myThid        INTEGER myThid
50    
51  C     == Local variables ==  C !OUTPUT PARAMETERS: ==================================================
52        INTEGER i,j,kp1,km1,km2,bi,bj  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    C  wLoc               :: velocity, vertical component
64          INTEGER i,j,kp1,km1,km2
65        _RL Cr,Rjm,Rj,Rjp        _RL Cr,Rjm,Rj,Rjp
66          _RL wLoc
67    C Statement function provides Limiter(Cr)
68  #include "GAD_FLUX_LIMITER.h"  #include "GAD_FLUX_LIMITER.h"
69    CEOP
       IF (.NOT. multiDimAdvection) THEN  
 C      If using the standard time-stepping/advection schemes (ie. AB-II)  
 C      then the data-structures are all global arrays  
        bi=bi_arg  
        bj=bj_arg  
       ELSE  
 C      otherwise if using the multi-dimensional advection schemes  
 C      then the data-structures are all local arrays except  
 C      for maskC(...) and wVel(...)  
        bi=1  
        bj=1  
       ENDIF  
70    
71        km2=MAX(1,k-2)        km2=MAX(1,k-2)
72        km1=MAX(1,k-1)        km1=MAX(1,k-1)
# Line 62  C      for maskC(...) and wVel(...) Line 81  C      for maskC(...) and wVel(...)
81        ELSE        ELSE
82         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
83          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
84           Rjp=(tracer(i,j,kp1,bi,bj)-tracer(i,j,k,bi,bj))  
85             wLoc = wFld(i,j)
86    c        wLoc = rTrans(i,j)*recip_rA(i,j,bi,bj)
87             Rjp=(tracer(i,j,kp1)-tracer(i,j,k))
88       &        *maskC(i,j,kp1,bi,bj)       &        *maskC(i,j,kp1,bi,bj)
89           Rj=(tracer(i,j,k,bi,bj)-tracer(i,j,kM1,bi,bj))           Rj= (tracer(i,j,k)  -tracer(i,j,kM1))
90           Rjm=(tracer(i,j,km1,bi,bj)-tracer(i,j,kM2,bi,bj))           Rjm=(tracer(i,j,km1)-tracer(i,j,kM2))
91       &        *maskC(i,j,km2,bi,bj)       &        *maskC(i,j,km2,bi,bj)
92    
93           IF (Rj.NE.0.) THEN           IF (Rj.NE.0.) THEN
94            IF (rTrans(i,j).LT.0) THEN            IF (rTrans(i,j).LT.0.) THEN
95              Cr=Rjm/Rj              Cr=Rjm/Rj
96            ELSE            ELSE
97              Cr=Rjp/Rj              Cr=Rjp/Rj
98            ENDIF            ENDIF
99           ELSE           ELSE
100            IF (rTrans(i,j).LT.0) THEN            IF (rTrans(i,j).LT.0.) THEN
101              Cr=Rjm*1.E20              Cr=Rjm*1.E20
102            ELSE            ELSE
103              Cr=Rjp*1.E20              Cr=Rjp*1.E20
104            ENDIF            ENDIF
105           ENDIF           ENDIF
106           Cr=Limiter(Cr)           Cr=Limiter(Cr)
107           wT(i,j) = maskC(i,j,kM1,bi_arg,bj_arg)*(           wT(i,j) = maskC(i,j,kM1,bi,bj)*(
108       &     rTrans(i,j)*       &     rTrans(i,j)*
109       &        (Tracer(i,j,k,bi,bj)+Tracer(i,j,kM1,bi,bj))*0.5 _d 0       &        (tracer(i,j,k)+tracer(i,j,kM1))*0.5 _d 0
110       &    +(ABS(rTrans(i,j))*(1-Cr)       &    +(ABS(rTrans(i,j))*(1-Cr)
111       &      +rTrans(i,j)*wVel(i,j,k,bi_arg,bj_arg)*dTarg*recip_drC(k)       &      +rTrans(i,j)*wLoc*dTarg*recip_drC(k)
112       &                  *Cr       &                  *Cr
113       &     )*Rj*0.5 _d 0                )       &     )*Rj*0.5 _d 0                )
114          ENDDO          ENDDO
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22