pkg/generic_advdiff/gad_dst3fl_adv_x.F

6d046cd062 Alis*0001 #include "GAD_OPTIONS.h"
                0002
0af3073e4e Jean*0003       SUBROUTINE GAD_DST3FL_ADV_X(
692dd30681 Jean*0004      I           bi,bj,k, calcCFL, deltaTloc,
0af3073e4e Jean*0005      I           uTrans, uFld,
38a970fb5d Jean*0006      I           maskLocW, tracer,
6d046cd062 Alis*0007      O           uT,
                0008      I           myThid )
                0009 C     /==========================================================\
                0010 C     | SUBROUTINE GAD_DST3FL_ADV_X                              |
                0011 C     | o Compute Zonal 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 uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0af3073e4e Jean*0026       _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
38a970fb5d Jean*0027       _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6d046cd062 Alis*0028       _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0029       _RL uT    (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,uCFL,d0,d1,psiP,psiM,thetaP,thetaM
20eb3a9e6b Jean*0035       _RL thetaMax
                0036       PARAMETER( thetaMax = 1.D+20 )
                0037
                0038 C- jmc: an alternative would be to compute directly psiM*Rj & psiP*Rj
                0039 C       (if Rj*Rjm < 0 => psiP*Rj = 0 , elsef Rj > 0 ... , else  ... )
                0040 C       with no need to compute thetaM (might be easier to differentiate)
6d046cd062 Alis*0041
64442af1fe Jean*0042       DO j=1-OLy,sNy+OLy
                0043        uT(1-OLx,j)=0. _d 0
                0044        uT(2-OLx,j)=0. _d 0
                0045        uT(sNx+OLx,j)=0. _d 0
e1940c5eb1 Mart*0046       ENDDO
64442af1fe Jean*0047       DO j=1-OLy,sNy+OLy
                0048        DO i=1-OLx+2,sNx+OLx-1
                0049 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
bf0222e9e9 Mart*0050 C     These lines make TAF create vectorizable code
                0051         thetaP = 0. _d 0
                0052         thetaM = 0. _d 0
                0053 #endif
38a970fb5d Jean*0054         Rjp=(tracer(i+1,j)-tracer( i ,j))*maskLocW(i+1,j)
                0055         Rj =(tracer( i ,j)-tracer(i-1,j))*maskLocW( i ,j)
                0056         Rjm=(tracer(i-1,j)-tracer(i-2,j))*maskLocW(i-1,j)
6d046cd062 Alis*0057
692dd30681 Jean*0058         uCFL = uFld(i,j)
                0059         IF ( calcCFL ) uCFL = ABS( uFld(i,j)*deltaTloc
62743db5c5 Jean*0060      &                  *recip_dxC(i,j,bi,bj)*recip_deepFacC(k) )
                0061         d0=(2. _d 0 -uCFL)*(1. _d 0 -uCFL)*oneSixth
                0062         d1=(1. _d 0 -uCFL*uCFL)*oneSixth
20eb3a9e6b Jean*0063
                0064 C-      the old version: can produce overflow, division by zero,
                0065 c       and is wrong for tracer with low concentration:
                0066 c       thetaP=Rjm/(1.D-20+Rj)
                0067 c       thetaM=Rjp/(1.D-20+Rj)
                0068 C-      the right expression, but not bounded:
51862bbf1d Patr*0069 c       thetaP=0.D0
                0070 c       thetaM=0.D0
20eb3a9e6b Jean*0071 c       IF (Rj.NE.0.D0) thetaP=Rjm/Rj
51862bbf1d Patr*0072 c       IF (Rj.NE.0.D0) thetaM=Rjp/Rj
20eb3a9e6b Jean*0073 C-      prevent |thetaP,M| to reach too big value:
                0074         IF ( ABS(Rj)*thetaMax .LE. ABS(Rjm) ) THEN
                0075           thetaP=SIGN(thetaMax,Rjm*Rj)
                0076         ELSE
                0077           thetaP=Rjm/Rj
                0078         ENDIF
                0079         IF ( ABS(Rj)*thetaMax .LE. ABS(Rjp) ) THEN
                0080           thetaM=SIGN(thetaMax,Rjp*Rj)
                0081         ELSE
                0082           thetaM=Rjp/Rj
                0083         ENDIF
                0084
                0085         psiP=d0+d1*thetaP
62743db5c5 Jean*0086         psiP=MAX(0. _d 0,MIN(MIN(1. _d 0,psiP),
                0087      &                       thetaP*(1. _d 0 -uCFL)/(uCFL+1. _d -20) ))
6d046cd062 Alis*0088         psiM=d0+d1*thetaM
62743db5c5 Jean*0089         psiM=MAX(0. _d 0,MIN(MIN(1. _d 0,psiM),
                0090      &                       thetaM*(1. _d 0 -uCFL)/(uCFL+1. _d -20) ))
20eb3a9e6b Jean*0091
6d046cd062 Alis*0092         uT(i,j)=
62743db5c5 Jean*0093      &   0.5*(uTrans(i,j)+ABS(uTrans(i,j)))
6d046cd062 Alis*0094      &      *( Tracer(i-1,j) + psiP*Rj )
62743db5c5 Jean*0095      &  +0.5*(uTrans(i,j)-ABS(uTrans(i,j)))
6d046cd062 Alis*0096      &      *( Tracer( i ,j) - psiM*Rj )
                0097
                0098        ENDDO
                0099       ENDDO
                0100
                0101       RETURN
                0102       END