pkg/smooth/smooth_impldiff.F

0c3d35c9cd Gael*0001 #include "SMOOTH_OPTIONS.h"
                0002
                0003       SUBROUTINE SMOOTH_IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax,
                0004      I                     deltaTX, KappaRX, recip_hFac,
                0005      U                     gXNm1,
                0006      I                     myThid )
                0007
                0008 C     *==========================================================*
                0009 C     | SUBROUTINE smooth_impldiff
                0010 C     | o Copy of model/src/impldiff
                0011 C     |    (simplified and with specified time step)
                0012 C     *==========================================================*
                0013
                0014 C     !USES:
                0015       IMPLICIT NONE
                0016 C     == Global data ==
                0017 #include "SIZE.h"
                0018 #include "EEPARAMS.h"
                0019 #include "PARAMS.h"
                0020 #include "GRID.h"
                0021
                0022 C     !INPUT/OUTPUT PARAMETERS:
9f5240b52a Jean*0023       INTEGER bi, bj, iMin, iMax, jMin, jMax
                0024       _RL deltaTX
                0025       _RL KappaRX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0026       _RS recip_hFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0027       _RL gXnm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
0c3d35c9cd Gael*0028       INTEGER myThid
                0029
                0030 C     !LOCAL VARIABLES:
                0031       INTEGER i,j,k
9f5240b52a Jean*0032       _RL gYnm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0033       _RL a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0034       _RL b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0035       _RL c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0036       _RL bet(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0037       _RL gam(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
0c3d35c9cd Gael*0038 CEOP
                0039
                0040 C--   Initialise
                0041       DO k=1,Nr
                0042        DO j=jMin,jMax
                0043         DO i=iMin,iMax
                0044          gYNm1(i,j,k,bi,bj) = 0. _d 0
                0045         ENDDO
                0046        ENDDO
                0047       ENDDO
                0048
                0049 C--   Old aLower
                0050       DO j=jMin,jMax
                0051        DO i=iMin,iMax
f3ec8d3513 Jean*0052          a(i,j,1) = 0. _d 0
0c3d35c9cd Gael*0053        ENDDO
                0054       ENDDO
                0055       DO k=2,Nr
                0056        DO j=jMin,jMax
                0057         DO i=iMin,iMax
                0058           a(i,j,k) = -deltaTX*recip_hFac(i,j,k,bi,bj)*recip_drF(k)
                0059      &               *recip_deepFac2C(k)*recip_rhoFacC(k)
                0060      &               *KappaRX(i,j, k )*recip_drC( k )
                0061      &               *deepFac2F(k)*rhoFacF(k)
                0062           IF (recip_hFac(i,j,k-1,bi,bj).EQ.0.) a(i,j,k)=0.
                0063         ENDDO
                0064        ENDDO
                0065       ENDDO
                0066
                0067 C--   Old aUpper
                0068       DO k=1,Nr-1
                0069        DO j=jMin,jMax
                0070         DO i=iMin,iMax
                0071           c(i,j,k) = -deltaTX*recip_hFac(i,j,k,bi,bj)*recip_drF(k)
                0072      &               *recip_deepFac2C(k)*recip_rhoFacC(k)
                0073      &               *KappaRX(i,j,k+1)*recip_drC(k+1)
                0074      &               *deepFac2F(k+1)*rhoFacF(k+1)
                0075           IF (recip_hFac(i,j,k+1,bi,bj).EQ.0.) c(i,j,k)=0.
                0076         ENDDO
                0077        ENDDO
                0078       ENDDO
                0079       DO j=jMin,jMax
                0080        DO i=iMin,iMax
                0081          c(i,j,Nr) = 0. _d 0
                0082        ENDDO
                0083       ENDDO
                0084
                0085 C--   Old aCenter
                0086       DO k=1,Nr
                0087        DO j=jMin,jMax
                0088         DO i=iMin,iMax
                0089           b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)
                0090         ENDDO
                0091        ENDDO
                0092       ENDDO
                0093
                0094 C--   Old and new gam, bet are the same
                0095       DO k=1,Nr
                0096        DO j=jMin,jMax
                0097         DO i=iMin,iMax
                0098           bet(i,j,k) = 1. _d 0
                0099           gam(i,j,k) = 0. _d 0
                0100         ENDDO
                0101        ENDDO
                0102       ENDDO
                0103
                0104 C--   Only need do anything if Nr>1
                0105       IF (Nr.GT.1) THEN
                0106
                0107        k = 1
                0108 C--    Beginning of forward sweep (top level)
                0109        DO j=jMin,jMax
                0110         DO i=iMin,iMax
                0111          IF (b(i,j,1).NE.0.) bet(i,j,1) = 1. _d 0 / b(i,j,1)
                0112         ENDDO
                0113        ENDDO
                0114
                0115       ENDIF
                0116
                0117 C--   Middle of forward sweep
                0118       IF (Nr.GE.2) THEN
                0119
                0120 CADJ loop = sequential
                0121        DO k=2,Nr
                0122
                0123         DO j=jMin,jMax
                0124          DO i=iMin,iMax
                0125           gam(i,j,k) = c(i,j,k-1)*bet(i,j,k-1)
f3ec8d3513 Jean*0126           IF ( ( b(i,j,k) - a(i,j,k)*gam(i,j,k) ) .NE. 0.)
0c3d35c9cd Gael*0127      &        bet(i,j,k) = 1. _d 0 / ( b(i,j,k) - a(i,j,k)*gam(i,j,k) )
                0128          ENDDO
                0129         ENDDO
                0130
                0131        ENDDO
                0132
                0133       ENDIF
                0134
                0135       DO j=jMin,jMax
                0136        DO i=iMin,iMax
                0137         gYNm1(i,j,1,bi,bj) = gXNm1(i,j,1,bi,bj)*bet(i,j,1)
                0138        ENDDO
                0139       ENDDO
                0140       DO k=2,Nr
                0141        DO j=jMin,jMax
                0142         DO i=iMin,iMax
                0143          gYnm1(i,j,k,bi,bj) = bet(i,j,k)*
                0144      &        (gXnm1(i,j,k,bi,bj) - a(i,j,k)*gYnm1(i,j,k-1,bi,bj))
                0145         ENDDO
                0146        ENDDO
                0147       ENDDO
                0148
                0149 C--    Backward sweep
                0150 CADJ loop = sequential
                0151        DO k=Nr-1,1,-1
                0152         DO j=jMin,jMax
                0153          DO i=iMin,iMax
                0154           gYnm1(i,j,k,bi,bj)=gYnm1(i,j,k,bi,bj)
                0155      &              -gam(i,j,k+1)*gYnm1(i,j,k+1,bi,bj)
                0156          ENDDO
                0157         ENDDO
                0158        ENDDO
                0159
                0160        DO k=1,Nr
                0161         DO j=jMin,jMax
                0162          DO i=iMin,iMax
                0163           gXnm1(i,j,k,bi,bj)=gYnm1(i,j,k,bi,bj)
                0164          ENDDO
                0165         ENDDO
                0166        ENDDO
                0167
                0168       RETURN
                0169       END