/[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.42 by heimbach, Tue Jul 18 03:23:50 2006 UTC revision 1.46 by adcroft, Sat Jan 20 21:20:11 2007 UTC
# Line 17  C !INTERFACE: ========================== Line 17  C !INTERFACE: ==========================
17       I     bi,bj, myTime,myIter,myThid)       I     bi,bj, myTime,myIter,myThid)
18    
19  C !DESCRIPTION:  C !DESCRIPTION:
20  C Calculates the tendancy of a tracer due to advection.  C Calculates the tendency of a tracer due to advection.
21  C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}  C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}
22  C and can only be used for the non-linear advection schemes such as the  C and can only be used for the non-linear advection schemes such as the
23  C direct-space-time method and flux-limiters.  C direct-space-time method and flux-limiters.
# Line 33  C      - \Delta t \partial_r (w\theta^{( Line 33  C      - \Delta t \partial_r (w\theta^{(
33  C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$}  C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$}
34  C \end{itemize}  C \end{itemize}
35  C  C
36  C The tendancy (output) is over-written by this routine.  C The tendency (output) is over-written by this routine.
37    
38  C !USES: ===============================================================  C !USES: ===============================================================
39        IMPLICIT NONE        IMPLICIT NONE
# Line 80  C  myThid            :: thread number Line 80  C  myThid            :: thread number
80        INTEGER myThid        INTEGER myThid
81    
82  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
83  C  gTracer           :: tendancy array  C  gTracer           :: tendency array
84        _RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
85    
86  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
# Line 383  CADJ &     comlev1_bibj_k_gad_pass, key= Line 383  CADJ &     comlev1_bibj_k_gad_pass, key=
383            CALL GAD_DST3FL_ADV_X(    bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_X(    bi,bj,k, dTtracerLev(k),
384       I                              uTrans, uFld, maskLocW, localTij,       I                              uTrans, uFld, maskLocW, localTij,
385       O                              af, myThid )       O                              af, myThid )
386            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
387              CALL GAD_OS7MP_ADV_X(     bi,bj,k, dTtracerLev(k),
388         I                              uTrans, uFld, maskLocW, localTij,
389         O                              af, myThid )
390          ELSE          ELSE
391           STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'           STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
392          ENDIF          ENDIF
# Line 411  C         in corner region) but safer to Line 415  C         in corner region) but safer to
415          IF ( S_edge ) THEN          IF ( S_edge ) THEN
416           DO j=1-Oly,0           DO j=1-Oly,0
417            DO i=iMinUpd,iMaxUpd            DO i=iMinUpd,iMaxUpd
418             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*             localTij(i,j) = localTij(i,j)
419       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &      -dTtracerLev(k)*recip_rhoFacC(k)
420       &       *recip_rA(i,j,bi,bj)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
421         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
422       &       *( af(i+1,j)-af(i,j)       &       *( af(i+1,j)-af(i,j)
423       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
424       &        )       &        )
# Line 423  C         in corner region) but safer to Line 428  C         in corner region) but safer to
428          IF ( N_edge ) THEN          IF ( N_edge ) THEN
429           DO j=sNy+1,sNy+Oly           DO j=sNy+1,sNy+Oly
430            DO i=iMinUpd,iMaxUpd            DO i=iMinUpd,iMaxUpd
431             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*             localTij(i,j) = localTij(i,j)
432       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &      -dTtracerLev(k)*recip_rhoFacC(k)
433       &       *recip_rA(i,j,bi,bj)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
434         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
435       &       *( af(i+1,j)-af(i,j)       &       *( af(i+1,j)-af(i,j)
436       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
437       &        )       &        )
# Line 441  C      do not only update the overlap Line 447  C      do not only update the overlap
447          IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy          IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
448          DO j=jMinUpd,jMaxUpd          DO j=jMinUpd,jMaxUpd
449           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
450             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*             localTij(i,j) = localTij(i,j)
451       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &      -dTtracerLev(k)*recip_rhoFacC(k)
452       &       *recip_rA(i,j,bi,bj)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
453         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
454       &       *( af(i+1,j)-af(i,j)       &       *( af(i+1,j)-af(i,j)
455       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
456       &        )       &        )
# Line 547  CADJ &     comlev1_bibj_k_gad_pass, key= Line 554  CADJ &     comlev1_bibj_k_gad_pass, key=
554            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, dTtracerLev(k),
555       I                              vTrans, vFld, maskLocS, localTij,       I                              vTrans, vFld, maskLocS, localTij,
556       O                              af, myThid )       O                              af, myThid )
557            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
558              CALL GAD_OS7MP_ADV_Y(     bi,bj,k, dTtracerLev(k),
559         I                              vTrans, vFld, maskLocS, localTij,
560         O                              af, myThid )
561          ELSE          ELSE
562           STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'           STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
563          ENDIF          ENDIF
# Line 575  C         in corner region) but safer to Line 586  C         in corner region) but safer to
586          IF ( W_edge ) THEN          IF ( W_edge ) THEN
587           DO j=jMinUpd,jMaxUpd           DO j=jMinUpd,jMaxUpd
588            DO i=1-Olx,0            DO i=1-Olx,0
589             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*             localTij(i,j) = localTij(i,j)
590       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &      -dTtracerLev(k)*recip_rhoFacC(k)
591       &       *recip_rA(i,j,bi,bj)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
592         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
593       &       *( af(i,j+1)-af(i,j)       &       *( af(i,j+1)-af(i,j)
594       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
595       &        )       &        )
# Line 587  C         in corner region) but safer to Line 599  C         in corner region) but safer to
599          IF ( E_edge ) THEN          IF ( E_edge ) THEN
600           DO j=jMinUpd,jMaxUpd           DO j=jMinUpd,jMaxUpd
601            DO i=sNx+1,sNx+Olx            DO i=sNx+1,sNx+Olx
602             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*             localTij(i,j) = localTij(i,j)
603       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &      -dTtracerLev(k)*recip_rhoFacC(k)
604       &       *recip_rA(i,j,bi,bj)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
605         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
606       &       *( af(i,j+1)-af(i,j)       &       *( af(i,j+1)-af(i,j)
607       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
608       &        )       &        )
# Line 605  C      do not only update the overlap Line 618  C      do not only update the overlap
618          IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx          IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
619          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
620           DO i=iMinUpd,iMaxUpd           DO i=iMinUpd,iMaxUpd
621             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*             localTij(i,j) = localTij(i,j)
622       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &      -dTtracerLev(k)*recip_rhoFacC(k)
623       &       *recip_rA(i,j,bi,bj)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
624         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
625       &       *( af(i,j+1)-af(i,j)       &       *( af(i,j+1)-af(i,j)
626       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
627       &        )       &        )
# Line 655  C-    explicit advection is done ; store Line 669  C-    explicit advection is done ; store
669          ENDDO          ENDDO
670        ELSE        ELSE
671  C-    horizontal advection done; store intermediate result in 3D array:  C-    horizontal advection done; store intermediate result in 3D array:
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           localTijk(i,j,k)=localTij(i,j)            localTijk(i,j,k)=localTij(i,j)
675             ENDDO
676          ENDDO          ENDDO
        ENDDO  
677        ENDIF        ENDIF
678    
679  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
# Line 743  C- Interior interface : Line 757  C- Interior interface :
757             rTransKp1(i,j) = kp1Msk*rTrans(i,j)             rTransKp1(i,j) = kp1Msk*rTrans(i,j)
758             wFld(i,j)   = wVel(i,j,k,bi,bj)             wFld(i,j)   = wVel(i,j,k,bi,bj)
759             rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)             rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
760         &                 *deepFac2F(k)*rhoFacF(k)
761       &                 *maskC(i,j,k-1,bi,bj)       &                 *maskC(i,j,k-1,bi,bj)
762             fVerT(i,j,kUp) = 0.             fVerT(i,j,kUp) = 0.
763            ENDDO            ENDDO
# Line 764  CADJ &     = comlev1_bibj_k_gad, key=kke Line 779  CADJ &     = comlev1_bibj_k_gad, key=kke
779  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
780    
781  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
782           IF ( advectionScheme.EQ.ENUM_UPWIND_1RST           IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
783       &      .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &      .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
784             CALL GAD_DST2U1_ADV_R( bi,bj,k, advectionScheme,             CALL GAD_DST2U1_ADV_R( bi,bj,k, advectionScheme,
785       I                            dTtracerLev(k),rTrans,wFld,localTijk,       I                            dTtracerLev(k),rTrans,wFld,localTijk,
786       O                            fVerT(1-Olx,1-Oly,kUp), myThid )       O                            fVerT(1-Olx,1-Oly,kUp), myThid )
787           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
788             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, dTtracerLev(k),             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, dTtracerLev(k),
789       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localTijk,
790       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-Olx,1-Oly,kUp), myThid )
791           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
792             CALL GAD_DST3_ADV_R(      bi,bj,k, dTtracerLev(k),             CALL GAD_DST3_ADV_R(      bi,bj,k, dTtracerLev(k),
793       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localTijk,
794       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-Olx,1-Oly,kUp), myThid )
795           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
796             CALL GAD_DST3FL_ADV_R(    bi,bj,k, dTtracerLev(k),             CALL GAD_DST3FL_ADV_R(    bi,bj,k, dTtracerLev(k),
797       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localTijk,
798       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-Olx,1-Oly,kUp), myThid )
799             ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
800               CALL GAD_OS7MP_ADV_R(     bi,bj,k, dTtracerLev(k),
801         I                               rTrans, wFld, localTijk,
802         O                               fVerT(1-Olx,1-Oly,kUp), myThid )
803           ELSE           ELSE
804            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
805           ENDIF           ENDIF
# Line 798  CADJ &     = comlev1_bibj_k_gad, key=kke Line 817  CADJ &     = comlev1_bibj_k_gad, key=kke
817  C--   Divergence of vertical fluxes  C--   Divergence of vertical fluxes
818          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
819           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
820            localTij(i,j)=localTijk(i,j,k)-dTtracerLev(k)*            localTij(i,j) = localTijk(i,j,k)
821       &     _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &      -dTtracerLev(k)*recip_rhoFacC(k)
822       &     *recip_rA(i,j,bi,bj)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
823       &     *( fVerT(i,j,kDown)-fVerT(i,j,kUp)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
824       &       -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))       &       *( fVerT(i,j,kDown)-fVerT(i,j,kUp)
825       &      )*rkSign       &         -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))
826         &        )*rkSign
827            gTracer(i,j,k,bi,bj)=            gTracer(i,j,k,bi,bj)=
828       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)
829           ENDDO           ENDDO
Legend:
Removed from v.1.42  
changed lines
  Added in v.1.46
  ViewVC Help
Powered by ViewVC 1.1.22