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

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

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

revision 1.60 by jmc, Tue May 12 19:56:35 2009 UTC revision 1.61 by jahn, Fri Jun 26 23:10:09 2009 UTC
# Line 11  C !ROUTINE: GAD_ADVECTION Line 11  C !ROUTINE: GAD_ADVECTION
11  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
12        SUBROUTINE GAD_ADVECTION(        SUBROUTINE GAD_ADVECTION(
13       I     implicitAdvection, advectionScheme, vertAdvecScheme,       I     implicitAdvection, advectionScheme, vertAdvecScheme,
14       I     tracerIdentity,       I     tracerIdentity, deltaTLev,
15       I     uVel, vVel, wVel, tracer,       I     uVel, vVel, wVel, tracer,
16       O     gTracer,       O     gTracer,
17       I     bi,bj, myTime,myIter,myThid)       I     bi,bj, myTime,myIter,myThid)
# Line 70  C  myThid            :: thread number Line 70  C  myThid            :: thread number
70        LOGICAL implicitAdvection        LOGICAL implicitAdvection
71        INTEGER advectionScheme, vertAdvecScheme        INTEGER advectionScheme, vertAdvecScheme
72        INTEGER tracerIdentity        INTEGER tracerIdentity
73          _RL deltaTLev(Nr)
74        _RL uVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL uVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
75        _RL vVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL vVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
76        _RL wVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL wVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
# Line 376  CADJ &     comlev1_bibj_k_gad_pass, key= Line 377  CADJ &     comlev1_bibj_k_gad_pass, key=
377          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
378       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
379            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
380       I                           dTtracerLev(k),uTrans,uFld,localTij,       I                           deltaTLev(k),uTrans,uFld,localTij,
381       O                           af, myThid )       O                           af, myThid )
382          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
383            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
384       I                              uTrans, uFld, maskLocW, localTij,       I                              uTrans, uFld, maskLocW, localTij,
385       O                              af, myThid )       O                              af, myThid )
386          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
387            CALL GAD_DST3_ADV_X(      bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_DST3_ADV_X(      bi,bj,k, .TRUE., deltaTLev(k),
388       I                              uTrans, uFld, maskLocW, localTij,       I                              uTrans, uFld, maskLocW, localTij,
389       O                              af, myThid )       O                              af, myThid )
390          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
391            CALL GAD_DST3FL_ADV_X(    bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_DST3FL_ADV_X(    bi,bj,k, .TRUE., deltaTLev(k),
392       I                              uTrans, uFld, maskLocW, localTij,       I                              uTrans, uFld, maskLocW, localTij,
393       O                              af, myThid )       O                              af, myThid )
394  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
395          ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN          ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
396            CALL GAD_OS7MP_ADV_X(     bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_OS7MP_ADV_X(     bi,bj,k, .TRUE., deltaTLev(k),
397       I                              uTrans, uFld, maskLocW, localTij,       I                              uTrans, uFld, maskLocW, localTij,
398       O                              af, myThid )       O                              af, myThid )
399  #endif  #endif
# Line 424  C         in corner region) but safer to Line 425  C         in corner region) but safer to
425           DO j=1-Oly,0           DO j=1-Oly,0
426            DO i=iMinUpd,iMaxUpd            DO i=iMinUpd,iMaxUpd
427             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
428       &      -dTtracerLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
429       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
430       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
431       &       *( af(i+1,j)-af(i,j)       &       *( af(i+1,j)-af(i,j)
# Line 437  C         in corner region) but safer to Line 438  C         in corner region) but safer to
438           DO j=sNy+1,sNy+Oly           DO j=sNy+1,sNy+Oly
439            DO i=iMinUpd,iMaxUpd            DO i=iMinUpd,iMaxUpd
440             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
441       &      -dTtracerLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
442       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
443       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
444       &       *( af(i+1,j)-af(i,j)       &       *( af(i+1,j)-af(i,j)
# Line 456  C      do not only update the overlap Line 457  C      do not only update the overlap
457          DO j=jMinUpd,jMaxUpd          DO j=jMinUpd,jMaxUpd
458           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
459             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
460       &      -dTtracerLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
461       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
462       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
463       &       *( af(i+1,j)-af(i,j)       &       *( af(i+1,j)-af(i,j)
# Line 546  CADJ &     comlev1_bibj_k_gad_pass, key= Line 547  CADJ &     comlev1_bibj_k_gad_pass, key=
547          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
548       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
549            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
550       I                           dTtracerLev(k),vTrans,vFld,localTij,       I                           deltaTLev(k),vTrans,vFld,localTij,
551       O                           af, myThid )       O                           af, myThid )
552          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
553            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
554       I                              vTrans, vFld, maskLocS, localTij,       I                              vTrans, vFld, maskLocS, localTij,
555       O                              af, myThid )       O                              af, myThid )
556          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
557            CALL GAD_DST3_ADV_Y(      bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_DST3_ADV_Y(      bi,bj,k, .TRUE., deltaTLev(k),
558       I                              vTrans, vFld, maskLocS, localTij,       I                              vTrans, vFld, maskLocS, localTij,
559       O                              af, myThid )       O                              af, myThid )
560          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
561            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, .TRUE., deltaTLev(k),
562       I                              vTrans, vFld, maskLocS, localTij,       I                              vTrans, vFld, maskLocS, localTij,
563       O                              af, myThid )       O                              af, myThid )
564  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
565          ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN          ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
566            CALL GAD_OS7MP_ADV_Y(     bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_OS7MP_ADV_Y(     bi,bj,k, .TRUE., deltaTLev(k),
567       I                              vTrans, vFld, maskLocS, localTij,       I                              vTrans, vFld, maskLocS, localTij,
568       O                              af, myThid )       O                              af, myThid )
569  #endif  #endif
# Line 594  C         in corner region) but safer to Line 595  C         in corner region) but safer to
595           DO j=jMinUpd,jMaxUpd           DO j=jMinUpd,jMaxUpd
596            DO i=1-Olx,0            DO i=1-Olx,0
597             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
598       &      -dTtracerLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
599       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
600       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
601       &       *( af(i,j+1)-af(i,j)       &       *( af(i,j+1)-af(i,j)
# Line 607  C         in corner region) but safer to Line 608  C         in corner region) but safer to
608           DO j=jMinUpd,jMaxUpd           DO j=jMinUpd,jMaxUpd
609            DO i=sNx+1,sNx+Olx            DO i=sNx+1,sNx+Olx
610             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
611       &      -dTtracerLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
612       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
613       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
614       &       *( af(i,j+1)-af(i,j)       &       *( af(i,j+1)-af(i,j)
# Line 626  C      do not only update the overlap Line 627  C      do not only update the overlap
627          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
628           DO i=iMinUpd,iMaxUpd           DO i=iMinUpd,iMaxUpd
629             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
630       &      -dTtracerLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
631       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
632       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
633       &       *( af(i,j+1)-af(i,j)       &       *( af(i,j+1)-af(i,j)
# Line 671  C-    explicit advection is done ; store Line 672  C-    explicit advection is done ; store
672          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
673           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
674            gTracer(i,j,k,bi,bj)=            gTracer(i,j,k,bi,bj)=
675       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)
676           ENDDO           ENDDO
677          ENDDO          ENDDO
678        ELSE        ELSE
# Line 794  C-    Compute vertical advective flux in Line 795  C-    Compute vertical advective flux in
795           IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST           IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
796       &      .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN       &      .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
797             CALL GAD_DST2U1_ADV_R( bi,bj,k, advectionScheme,             CALL GAD_DST2U1_ADV_R( bi,bj,k, advectionScheme,
798       I                            dTtracerLev(k),rTrans,wFld,localTijk,       I                            deltaTLev(k),rTrans,wFld,localTijk,
799       O                            fVerT(1-Olx,1-Oly,kUp), myThid )       O                            fVerT(1-Olx,1-Oly,kUp), myThid )
800           ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
801             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, dTtracerLev(k),             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
802       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localTijk,
803       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-Olx,1-Oly,kUp), myThid )
804           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
805             CALL GAD_DST3_ADV_R(      bi,bj,k, dTtracerLev(k),             CALL GAD_DST3_ADV_R(      bi,bj,k, deltaTLev(k),
806       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localTijk,
807       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-Olx,1-Oly,kUp), myThid )
808           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
809             CALL GAD_DST3FL_ADV_R(    bi,bj,k, dTtracerLev(k),             CALL GAD_DST3FL_ADV_R(    bi,bj,k, deltaTLev(k),
810       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localTijk,
811       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-Olx,1-Oly,kUp), myThid )
812  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
813           ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN           ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
814             CALL GAD_OS7MP_ADV_R(     bi,bj,k, dTtracerLev(k),             CALL GAD_OS7MP_ADV_R(     bi,bj,k, deltaTLev(k),
815       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localTijk,
816       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-Olx,1-Oly,kUp), myThid )
817  #endif  #endif
# Line 838  C--   Divergence of vertical fluxes Line 839  C--   Divergence of vertical fluxes
839          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
840           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
841            localTij(i,j) = localTijk(i,j,k)            localTij(i,j) = localTijk(i,j,k)
842       &      -dTtracerLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
843       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
844       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
845       &       *( fVerT(i,j,kDown)-fVerT(i,j,kUp)       &       *( fVerT(i,j,kDown)-fVerT(i,j,kUp)
846       &         -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))
847       &        )*rkSign       &        )*rkSign
848            gTracer(i,j,k,bi,bj)=            gTracer(i,j,k,bi,bj)=
849       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)
850           ENDDO           ENDDO
851          ENDDO          ENDDO
852    
Legend:
Removed from v.1.60  
changed lines
  Added in v.1.61
  ViewVC Help
Powered by ViewVC 1.1.22