pkg/generic_advdiff/gad_dst3fl_adv_y.F

6d046cd062 Alis*0001 #include "GAD_OPTIONS.h"
                0002
0af3073e4e Jean*0003       SUBROUTINE GAD_DST3FL_ADV_Y(
692dd30681 Jean*0004      I           bi,bj,k, calcCFL, deltaTloc,
0af3073e4e Jean*0005      I           vTrans, vFld,
38a970fb5d Jean*0006      I           maskLocS, tracer,
6d046cd062 Alis*0007      O           vT,
                0008      I           myThid )
                0009 C     /==========================================================\
                0010 C     | SUBROUTINE GAD_DST3FL_ADV_Y                              |
                0011 C     | o Compute Meridional advective Flux of Tracer using      |
                0012 C     |   3rd Order DST Sceheme with flux limiting               |
                0013 C     |==========================================================|
                0014       IMPLICIT NONE
                0015
                0016 C     == GLobal variables ==
                0017 #include "SIZE.h"
                0018 #include "GRID.h"
                0019 #include "GAD.h"
                0020
                0021 C     == Routine arguments ==
                0022       INTEGER bi,bj,k
692dd30681 Jean*0023       LOGICAL calcCFL
1816bb0a69 Patr*0024       _RL deltaTloc
6d046cd062 Alis*0025       _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0af3073e4e Jean*0026       _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
38a970fb5d Jean*0027       _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6d046cd062 Alis*0028       _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0029       _RL vT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0030       INTEGER myThid
                0031
                0032 C     == Local variables ==
                0033       INTEGER i,j
62743db5c5 Jean*0034       _RL Rjm,Rj,Rjp,vCFL,d0,d1,psiP,psiM,thetaP,thetaM
20eb3a9e6b Jean*0035       _RL thetaMax
                0036       PARAMETER( thetaMax = 1.D+20 )
6d046cd062 Alis*0037
64442af1fe Jean*0038       DO i=1-OLx,sNx+OLx
                0039        vT(i,1-OLy)=0. _d 0
                0040        vT(i,2-OLy)=0. _d 0
                0041        vT(i,sNy+OLy)=0. _d 0
6d046cd062 Alis*0042       ENDDO
64442af1fe Jean*0043       DO j=1-OLy+2,sNy+OLy-1
                0044        DO i=1-OLx,sNx+OLx
                0045 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
bf0222e9e9 Mart*0046 C     These lines make TAF create vectorizable code
                0047         thetaP = 0. _d 0
                0048         thetaM = 0. _d 0
                0049 #endif
38a970fb5d Jean*0050         Rjp=(tracer(i,j+1)-tracer(i, j ))*maskLocS(i,j+1)
                0051         Rj =(tracer(i, j )-tracer(i,j-1))*maskLocS(i, j )
                0052         Rjm=(tracer(i,j-1)-tracer(i,j-2))*maskLocS(i,j-1)
6d046cd062 Alis*0053
692dd30681 Jean*0054         vCFL = vFld(i,j)
                0055         IF ( calcCFL ) vCFL = ABS( vFld(i,j)*deltaTloc
62743db5c5 Jean*0056      &                  *recip_dyC(i,j,bi,bj)*recip_deepFacC(k) )
                0057         d0=(2. _d 0 -vCFL)*(1. _d 0 -vCFL)*oneSixth
                0058         d1=(1. _d 0 -vCFL*vCFL)*oneSixth
20eb3a9e6b Jean*0059
                0060 C-      the old version: can produce overflow, division by zero,
                0061 c       and is wrong for tracer with low concentration:
                0062 c       thetaP=Rjm/(1.D-20+Rj)
                0063 c       thetaM=Rjp/(1.D-20+Rj)
                0064 C-      the right expression, but not bounded:
51862bbf1d Patr*0065 c       thetaP=0.D0
                0066 c       thetaM=0.D0
20eb3a9e6b Jean*0067 c       IF (Rj.NE.0.D0) thetaP=Rjm/Rj
51862bbf1d Patr*0068 c       IF (Rj.NE.0.D0) thetaM=Rjp/Rj
20eb3a9e6b Jean*0069 C-      prevent |thetaP,M| to reach too big value:
                0070         IF ( ABS(Rj)*thetaMax .LE. ABS(Rjm) ) THEN
                0071           thetaP=SIGN(thetaMax,Rjm*Rj)
                0072         ELSE
                0073           thetaP=Rjm/Rj
                0074         ENDIF
                0075         IF ( ABS(Rj)*thetaMax .LE. ABS(Rjp) ) THEN
                0076           thetaM=SIGN(thetaMax,Rjp*Rj)
                0077         ELSE
                0078           thetaM=Rjp/Rj
                0079         ENDIF
                0080
                0081         psiP=d0+d1*thetaP
62743db5c5 Jean*0082         psiP=MAX(0. _d 0,MIN(MIN(1. _d 0,psiP),
                0083      &                       thetaP*(1. _d 0 -vCFL)/(vCFL+1. _d -20) ))
6d046cd062 Alis*0084         psiM=d0+d1*thetaM
62743db5c5 Jean*0085         psiM=MAX(0. _d 0,MIN(MIN(1. _d 0,psiM),
                0086      &                       thetaM*(1. _d 0 -vCFL)/(vCFL+1. _d -20) ))
20eb3a9e6b Jean*0087
6d046cd062 Alis*0088         vT(i,j)=
62743db5c5 Jean*0089      &   0.5*(vTrans(i,j)+ABS(vTrans(i,j)))
6d046cd062 Alis*0090      &      *( Tracer(i,j-1) + psiP*Rj )
62743db5c5 Jean*0091      &  +0.5*(vTrans(i,j)-ABS(vTrans(i,j)))
6d046cd062 Alis*0092      &      *( Tracer(i, j ) - psiM*Rj )
                0093
                0094        ENDDO
                0095       ENDDO
                0096
                0097       RETURN
                0098       END