/[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.70 by mlosch, Mon Dec 12 15:33:53 2011 UTC revision 1.71 by jmc, Fri Mar 9 20:13:44 2012 UTC
# Line 10  C !ROUTINE: GAD_ADVECTION Line 10  C !ROUTINE: GAD_ADVECTION
10  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
11        SUBROUTINE GAD_ADVECTION(        SUBROUTINE GAD_ADVECTION(
12       I     implicitAdvection, advectionScheme, vertAdvecScheme,       I     implicitAdvection, advectionScheme, vertAdvecScheme,
13       I     tracerIdentity, deltaTLev,       I     trIdentity, deltaTLev,
14       I     uVel, vVel, wVel, tracer,       I     uVel, vVel, wVel, tracer,
15       O     gTracer,       O     gTracer,
16       I     bi,bj, myTime,myIter,myThid)       I     bi,bj, myTime,myIter,myThid)
# Line 57  C !INPUT PARAMETERS: =================== Line 57  C !INPUT PARAMETERS: ===================
57  C  implicitAdvection :: implicit vertical advection (later on)  C  implicitAdvection :: implicit vertical advection (later on)
58  C  advectionScheme   :: advection scheme to use (Horizontal plane)  C  advectionScheme   :: advection scheme to use (Horizontal plane)
59  C  vertAdvecScheme   :: advection scheme to use (vertical direction)  C  vertAdvecScheme   :: advection scheme to use (vertical direction)
60  C  tracerIdentity    :: tracer identifier (required only for OBCS)  C  trIdentity        :: tracer identifier
61  C  uVel              :: velocity, zonal component  C  uVel              :: velocity, zonal component
62  C  vVel              :: velocity, meridional component  C  vVel              :: velocity, meridional component
63  C  wVel              :: velocity, vertical component  C  wVel              :: velocity, vertical component
# Line 68  C  myIter            :: iteration number Line 68  C  myIter            :: iteration number
68  C  myThid            :: thread number  C  myThid            :: thread number
69        LOGICAL implicitAdvection        LOGICAL implicitAdvection
70        INTEGER advectionScheme, vertAdvecScheme        INTEGER advectionScheme, vertAdvecScheme
71        INTEGER tracerIdentity        INTEGER trIdentity
72        _RL deltaTLev(Nr)        _RL deltaTLev(Nr)
73        _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)
74        _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)
75        _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)
76        _RL tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
77        INTEGER bi,bj        INTEGER bi,bj
78        _RL myTime        _RL myTime
79        INTEGER myIter        INTEGER myIter
# Line 81  C  myThid            :: thread number Line 81  C  myThid            :: thread number
81    
82  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
83  C  gTracer           :: tendency 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: ====================================================
87  C  maskUp        :: 2-D array for mask at W points  C  maskUp        :: 2-D array for mask at W points
# Line 162  C-    Functions: Line 162  C-    Functions:
162  CEOP  CEOP
163    
164  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
165            act0 = tracerIdentity            act0 = trIdentity
166            max0 = maxpass            max0 = maxpass
167            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
168            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
# Line 176  CEOP Line 176  CEOP
176       &                      + act2*max0*max1       &                      + act2*max0*max1
177       &                      + act3*max0*max1*max2       &                      + act3*max0*max1*max2
178       &                      + act4*max0*max1*max2*max3       &                      + act4*max0*max1*max2*max3
179            IF (tracerIdentity.GT.maxpass) THEN            IF (trIdentity.GT.maxpass) THEN
180             WRITE(msgBuf,'(A,2I3)')             WRITE(msgBuf,'(A,2I3)')
181       &      'GAD_ADVECTION: maxpass < tracerIdentity ',       &      'GAD_ADVECTION: maxpass < trIdentity ',
182       &      maxpass, tracerIdentity       &      maxpass, trIdentity
183             CALL PRINT_ERROR( msgBuf, myThid )             CALL PRINT_ERROR( msgBuf, myThid )
184             STOP 'ABNORMAL END: S/R GAD_ADVECTION'             STOP 'ABNORMAL END: S/R GAD_ADVECTION'
185            ENDIF            ENDIF
# Line 191  C--   Set diagnostics flags and suffix f Line 191  C--   Set diagnostics flags and suffix f
191        doDiagAdvY = .FALSE.        doDiagAdvY = .FALSE.
192        doDiagAdvR = .FALSE.        doDiagAdvR = .FALSE.
193        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
194          diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )          diagSufx = GAD_DIAG_SUFX( trIdentity, myThid )
195          diagName = 'ADVx'//diagSufx          diagName = 'ADVx'//diagSufx
196          doDiagAdvX = DIAGNOSTICS_IS_ON( diagName, myThid )          doDiagAdvX = DIAGNOSTICS_IS_ON( diagName, myThid )
197          diagName = 'ADVy'//diagSufx          diagName = 'ADVy'//diagSufx
# Line 277  C--   Residual transp = Bolus transp + E Line 277  C--   Residual transp = Bolus transp + E
277  C--   Make local copy of tracer array and mask West & South  C--   Make local copy of tracer array and mask West & South
278        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
279         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
280           localTij(i,j)=tracer(i,j,k,bi,bj)           localTij(i,j) = tracer(i,j,k,bi,bj)
281  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
282           maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)           maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
283           maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)           maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
# Line 336  C-    not CubedSphere Line 336  C-    not CubedSphere
336  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
337  C--   X direction  C--   X direction
338  C-     Advective flux in X  C-     Advective flux in X
339          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
340           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
341            af(i,j) = 0.            af(i,j) = 0.
342           ENDDO           ENDDO
343          ENDDO          ENDDO
# Line 373  CADJ &     comlev1_bibj_k_gad_pass, key= Line 373  CADJ &     comlev1_bibj_k_gad_pass, key=
373    
374          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
375       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
376            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,            CALL GAD_DST2U1_ADV_X(    bi,bj,k, advectionScheme, .TRUE.,
377       I                           deltaTLev(k),uTrans,uFld,localTij,       I                              deltaTLev(k),uTrans,uFld,localTij,
378       O                           af, myThid )       O                              af, myThid )
379          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
380            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
381       I                              uTrans, uFld, maskLocW, localTij,       I                              uTrans, uFld, maskLocW, localTij,
# Line 398  CADJ &     comlev1_bibj_k_gad_pass, key= Line 398  CADJ &     comlev1_bibj_k_gad_pass, key=
398           STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'           STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
399          ENDIF          ENDIF
400    
401    #ifdef ALLOW_OBCS
402            IF ( useOBCS ) THEN
403    C-      replace advective flux with 1st order upwind scheme estimate
404              CALL OBCS_U1_ADV_TRACER( .TRUE., trIdentity, bi, bj, k,
405         I                             maskLocW, uTrans, localTij,
406         U                             af, myThid )
407            ENDIF
408    #endif /* ALLOW_OBCS */
409    
410  C-     Internal exchange for next calculations in Y  C-     Internal exchange for next calculations in Y
411          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
412           CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,           CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
# Line 412  C      use "maksInC" to prevent updating Line 421  C      use "maksInC" to prevent updating
421    
422  C      update in overlap-Only  C      update in overlap-Only
423         IF ( overlapOnly ) THEN         IF ( overlapOnly ) THEN
424          iMinUpd = 1-Olx+1          iMinUpd = 1-OLx+1
425          iMaxUpd = sNx+Olx-1          iMaxUpd = sNx+OLx-1
426  C- notes: these 2 lines below have no real effect (because recip_hFac=0  C- notes: these 2 lines below have no real effect (because recip_hFac=0
427  C         in corner region) but safer to keep them.  C         in corner region) but safer to keep them.
428          IF ( W_edge ) iMinUpd = 1          IF ( W_edge ) iMinUpd = 1
429          IF ( E_edge ) iMaxUpd = sNx          IF ( E_edge ) iMaxUpd = sNx
430    
431          IF ( S_edge ) THEN          IF ( S_edge ) THEN
432           DO j=1-Oly,0           DO j=1-OLy,0
433            DO i=iMinUpd,iMaxUpd            DO i=iMinUpd,iMaxUpd
434             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
435       &      -deltaTLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
# Line 433  C         in corner region) but safer to Line 442  C         in corner region) but safer to
442           ENDDO           ENDDO
443          ENDIF          ENDIF
444          IF ( N_edge ) THEN          IF ( N_edge ) THEN
445           DO j=sNy+1,sNy+Oly           DO j=sNy+1,sNy+OLy
446            DO i=iMinUpd,iMaxUpd            DO i=iMinUpd,iMaxUpd
447             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
448       &      -deltaTLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
# Line 448  C         in corner region) but safer to Line 457  C         in corner region) but safer to
457    
458         ELSE         ELSE
459  C      do not only update the overlap  C      do not only update the overlap
460          jMinUpd = 1-Oly          jMinUpd = 1-OLy
461          jMaxUpd = sNy+Oly          jMaxUpd = sNy+OLy
462          IF ( interiorOnly .AND. S_edge ) jMinUpd = 1          IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
463          IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy          IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
464          DO j=jMinUpd,jMaxUpd          DO j=jMinUpd,jMaxUpd
465           DO i=1-Olx+1,sNx+Olx-1           DO i=1-OLx+1,sNx+OLx-1
466             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
467       &      -deltaTLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
468       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
# Line 464  C      do not only update the overlap Line 473  C      do not only update the overlap
473           ENDDO           ENDDO
474          ENDDO          ENDDO
475  C-      keep advective flux (for diagnostics)  C-      keep advective flux (for diagnostics)
476          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
477           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
478            afx(i,j) = af(i,j)            afx(i,j) = af(i,j)
479           ENDDO           ENDDO
480          ENDDO          ENDDO
# Line 480  C---+----1----+----2----+----3----+----4 Line 489  C---+----1----+----2----+----3----+----4
489  C--   Y direction  C--   Y direction
490  cph-test  cph-test
491  C-     Advective flux in Y  C-     Advective flux in Y
492          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
493           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
494            af(i,j) = 0.            af(i,j) = 0.
495           ENDDO           ENDDO
496          ENDDO          ENDDO
# Line 509  C-     Internal exchange for calculation Line 518  C-     Internal exchange for calculation
518          ENDIF          ENDIF
519    
520  C-     Advective flux in Y  C-     Advective flux in Y
521          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
522           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
523            af(i,j) = 0.            af(i,j) = 0.
524           ENDDO           ENDDO
525          ENDDO          ENDDO
# Line 524  CADJ &     comlev1_bibj_k_gad_pass, key= Line 533  CADJ &     comlev1_bibj_k_gad_pass, key=
533    
534          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
535       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
536            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,            CALL GAD_DST2U1_ADV_Y(    bi,bj,k, advectionScheme, .TRUE.,
537       I                           deltaTLev(k),vTrans,vFld,localTij,       I                              deltaTLev(k),vTrans,vFld,localTij,
538       O                           af, myThid )       O                              af, myThid )
539          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
540            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
541       I                              vTrans, vFld, maskLocS, localTij,       I                              vTrans, vFld, maskLocS, localTij,
# Line 549  CADJ &     comlev1_bibj_k_gad_pass, key= Line 558  CADJ &     comlev1_bibj_k_gad_pass, key=
558           STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'           STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
559          ENDIF          ENDIF
560    
561    #ifdef ALLOW_OBCS
562            IF ( useOBCS ) THEN
563    C-      replace advective flux with 1st order upwind scheme estimate
564              CALL OBCS_U1_ADV_TRACER( .FALSE., trIdentity, bi, bj, k,
565         I                             maskLocS, vTrans, localTij,
566         U                             af, myThid )
567            ENDIF
568    #endif /* ALLOW_OBCS */
569    
570  C-     Internal exchange for next calculations in X  C-     Internal exchange for next calculations in X
571          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
572           CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,           CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
# Line 563  C      use "maksInC" to prevent updating Line 581  C      use "maksInC" to prevent updating
581    
582  C      update in overlap-Only  C      update in overlap-Only
583         IF ( overlapOnly ) THEN         IF ( overlapOnly ) THEN
584          jMinUpd = 1-Oly+1          jMinUpd = 1-OLy+1
585          jMaxUpd = sNy+Oly-1          jMaxUpd = sNy+OLy-1
586  C- notes: these 2 lines below have no real effect (because recip_hFac=0  C- notes: these 2 lines below have no real effect (because recip_hFac=0
587  C         in corner region) but safer to keep them.  C         in corner region) but safer to keep them.
588          IF ( S_edge ) jMinUpd = 1          IF ( S_edge ) jMinUpd = 1
# Line 572  C         in corner region) but safer to Line 590  C         in corner region) but safer to
590    
591          IF ( W_edge ) THEN          IF ( W_edge ) THEN
592           DO j=jMinUpd,jMaxUpd           DO j=jMinUpd,jMaxUpd
593            DO i=1-Olx,0            DO i=1-OLx,0
594             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
595       &      -deltaTLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
596       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
# Line 585  C         in corner region) but safer to Line 603  C         in corner region) but safer to
603          ENDIF          ENDIF
604          IF ( E_edge ) THEN          IF ( E_edge ) THEN
605           DO j=jMinUpd,jMaxUpd           DO j=jMinUpd,jMaxUpd
606            DO i=sNx+1,sNx+Olx            DO i=sNx+1,sNx+OLx
607             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
608       &      -deltaTLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
609       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
# Line 599  C         in corner region) but safer to Line 617  C         in corner region) but safer to
617    
618         ELSE         ELSE
619  C      do not only update the overlap  C      do not only update the overlap
620          iMinUpd = 1-Olx          iMinUpd = 1-OLx
621          iMaxUpd = sNx+Olx          iMaxUpd = sNx+OLx
622          IF ( interiorOnly .AND. W_edge ) iMinUpd = 1          IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
623          IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx          IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
624          DO j=1-Oly+1,sNy+Oly-1          DO j=1-OLy+1,sNy+OLy-1
625           DO i=iMinUpd,iMaxUpd           DO i=iMinUpd,iMaxUpd
626             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
627       &      -deltaTLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
# Line 615  C      do not only update the overlap Line 633  C      do not only update the overlap
633           ENDDO           ENDDO
634          ENDDO          ENDDO
635  C-      keep advective flux (for diagnostics)  C-      keep advective flux (for diagnostics)
636          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
637           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
638            afy(i,j) = af(i,j)            afy(i,j) = af(i,j)
639           ENDDO           ENDDO
640          ENDDO          ENDDO
# Line 632  C--   End of ipass loop Line 650  C--   End of ipass loop
650    
651        IF ( implicitAdvection ) THEN        IF ( implicitAdvection ) THEN
652  C-    explicit advection is done ; store tendency in gTracer:  C-    explicit advection is done ; store tendency in gTracer:
653          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
654           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
655            gTracer(i,j,k,bi,bj)=            gTracer(i,j,k,bi,bj)=
656       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)
657           ENDDO           ENDDO
658          ENDDO          ENDDO
659        ELSE        ELSE
660  C-    horizontal advection done; store intermediate result in 3D array:  C-    horizontal advection done; store intermediate result in 3D array:
661          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
662           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
663            localTijk(i,j,k)=localTij(i,j)            localTijk(i,j,k) = localTij(i,j)
664           ENDDO           ENDDO
665          ENDDO          ENDDO
666        ENDIF        ENDIF
# Line 650  C-    horizontal advection done; store i Line 668  C-    horizontal advection done; store i
668  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
669          IF ( doDiagAdvX ) THEN          IF ( doDiagAdvX ) THEN
670            diagName = 'ADVx'//diagSufx            diagName = 'ADVx'//diagSufx
671            CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( afx, diagName, k,1, 2,bi,bj, myThid )
672          ENDIF          ENDIF
673          IF ( doDiagAdvY ) THEN          IF ( doDiagAdvY ) THEN
674            diagName = 'ADVy'//diagSufx            diagName = 'ADVy'//diagSufx
675            CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( afy, diagName, k,1, 2,bi,bj, myThid )
676          ENDIF          ENDIF
677  #endif  #endif
678    
679  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
680        IF ( debugLevel .GE. debLevC        IF ( debugLevel .GE. debLevC
681       &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE       &   .AND. trIdentity.EQ.GAD_TEMPERATURE
682       &   .AND. k.LE.3 .AND. myIter.EQ.1+nIter0       &   .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
683       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
684       &   .AND. useCubedSphereExchange ) THEN       &   .AND. useCubedSphereExchange ) THEN
# Line 699  C-- Compute Vertical transport Line 717  C-- Compute Vertical transport
717  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
718  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
719          IF ( k.EQ.1 .OR.          IF ( k.EQ.1 .OR.
720       &     (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)       &     (useAIM .AND. trIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
721       &              ) THEN       &              ) THEN
722  #else  #else
723          IF ( k.EQ.1 ) THEN          IF ( k.EQ.1 ) THEN
# Line 711  cphmultiCADJ &     comlev1_bibj_k_gad, k Line 729  cphmultiCADJ &     comlev1_bibj_k_gad, k
729  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
730    
731  C- Surface interface :  C- Surface interface :
732           DO j=1-Oly,sNy+Oly           DO j=1-OLy,sNy+OLy
733            DO i=1-Olx,sNx+Olx            DO i=1-OLx,sNx+OLx
734             rTransKp1(i,j) = kp1Msk*rTrans(i,j)             rTransKp1(i,j) = kp1Msk*rTrans(i,j)
735             wFld(i,j)   = 0.             wFld(i,j)   = 0.
736             rTrans(i,j) = 0.             rTrans(i,j) = 0.
# Line 728  cphmultiCADJ &     comlev1_bibj_k_gad, k Line 746  cphmultiCADJ &     comlev1_bibj_k_gad, k
746  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
747    
748  C- Interior interface :  C- Interior interface :
749           DO j=1-Oly,sNy+Oly           DO j=1-OLy,sNy+OLy
750            DO i=1-Olx,sNx+Olx            DO i=1-OLx,sNx+OLx
751             rTransKp1(i,j) = kp1Msk*rTrans(i,j)             rTransKp1(i,j) = kp1Msk*rTrans(i,j)
752             wFld(i,j)   = wVel(i,j,k,bi,bj)             wFld(i,j)   = wVel(i,j,k,bi,bj)
753             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)
# Line 757  cphmultiCADJ &     = comlev1_bibj_k_gad, Line 775  cphmultiCADJ &     = comlev1_bibj_k_gad,
775  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
776           IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST           IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
777       &      .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN       &      .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
778             CALL GAD_DST2U1_ADV_R( bi,bj,k, advectionScheme,             CALL GAD_DST2U1_ADV_R(    bi,bj,k, advectionScheme,
779       I                            deltaTLev(k),rTrans,wFld,localTijk,       I                          deltaTLev(k),rTrans,wFld,localTijk,
780       O                            fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-OLx,1-OLy,kUp), myThid )
781           ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
782             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
783       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localTijk,
784       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-OLx,1-OLy,kUp), myThid )
785           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
786             CALL GAD_DST3_ADV_R(      bi,bj,k, deltaTLev(k),             CALL GAD_DST3_ADV_R(      bi,bj,k, deltaTLev(k),
787       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localTijk,
788       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-OLx,1-OLy,kUp), myThid )
789           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
790             CALL GAD_DST3FL_ADV_R(    bi,bj,k, deltaTLev(k),             CALL GAD_DST3FL_ADV_R(    bi,bj,k, deltaTLev(k),
791       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localTijk,
792       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-OLx,1-OLy,kUp), myThid )
793  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
794           ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN           ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
795             CALL GAD_OS7MP_ADV_R(     bi,bj,k, deltaTLev(k),             CALL GAD_OS7MP_ADV_R(     bi,bj,k, deltaTLev(k),
796       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localTijk,
797       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-OLx,1-OLy,kUp), myThid )
798  #endif  #endif
799           ELSE           ELSE
800            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
# Line 799  CADJ &     = comlev1_bibj_k_gad, key=kke Line 817  CADJ &     = comlev1_bibj_k_gad, key=kke
817  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
818    
819  C--   Divergence of vertical fluxes  C--   Divergence of vertical fluxes
820          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
821           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
822            localTij(i,j) = localTijk(i,j,k)            localTij(i,j) = localTijk(i,j,k)
823       &      -deltaTLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
824       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
# Line 816  C--   Divergence of vertical fluxes Line 834  C--   Divergence of vertical fluxes
834  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
835          IF ( doDiagAdvR ) THEN          IF ( doDiagAdvR ) THEN
836            diagName = 'ADVr'//diagSufx            diagName = 'ADVr'//diagSufx
837            CALL DIAGNOSTICS_FILL( fVerT(1-Olx,1-Oly,kUp),            CALL DIAGNOSTICS_FILL( fVerT(1-OLx,1-OLy,kUp),
838       &                           diagName, k,1, 2,bi,bj, myThid)       &                           diagName, k,1, 2,bi,bj, myThid )
839          ENDIF          ENDIF
840  #endif  #endif
841    
Legend:
Removed from v.1.70  
changed lines
  Added in v.1.71
  ViewVC Help
Powered by ViewVC 1.1.22