/[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.45 by jmc, Wed Jan 10 18:53:25 2007 UTC revision 1.54 by mlosch, Thu Feb 7 08:52:12 2008 UTC
# Line 156  c     _RS maskUp  (1-OLx:sNx+OLx,1-OLy:s Line 156  c     _RS maskUp  (1-OLx:sNx+OLx,1-OLy:s
156  CEOP  CEOP
157    
158  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
159            act0 = tracerIdentity - 1            act0 = tracerIdentity
160            max0 = maxpass            max0 = maxpass
161            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
162            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
# Line 165  CEOP Line 165  CEOP
165            act3 = myThid - 1            act3 = myThid - 1
166            max3 = nTx*nTy            max3 = nTx*nTy
167            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
168            igadkey = (act0 + 1)            igadkey = act0
169       &                      + act1*max0       &                      + act1*max0
170       &                      + act2*max0*max1       &                      + act2*max0*max1
171       &                      + act3*max0*max1*max2       &                      + act3*max0*max1*max2
# Line 265  C--   Make local copy of tracer array an Line 265  C--   Make local copy of tracer array an
265        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
266         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
267           localTij(i,j)=tracer(i,j,k,bi,bj)           localTij(i,j)=tracer(i,j,k,bi,bj)
268           maskLocW(i,j)=maskW(i,j,k,bi,bj)           maskLocW(i,j)=_maskW(i,j,k,bi,bj)
269           maskLocS(i,j)=maskS(i,j,k,bi,bj)           maskLocS(i,j)=_maskS(i,j,k,bi,bj)
270         ENDDO         ENDDO
271        ENDDO        ENDDO
272    
273  #ifndef ALLOW_AUTODIFF_TAMC  cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
274        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
275          withSigns = .FALSE.          withSigns = .FALSE.
276          CALL FILL_CS_CORNER_UV_RS(          CALL FILL_CS_CORNER_UV_RS(
277       &            withSigns, maskLocW,maskLocS, bi,bj, myThid )       &            withSigns, maskLocW,maskLocS, bi,bj, myThid )
278        ENDIF        ENDIF
279  #endif  cph-exch2#endif
280    
281  C--   Multiple passes for different directions on different tiles  C--   Multiple passes for different directions on different tiles
282  C--   For cube need one pass for each of red, green and blue axes.  C--   For cube need one pass for each of red, green and blue axes.
283        DO ipass=1,nipass        DO ipass=1,nipass
284  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
285           passkey = ipass + (k-1)      *maxcube           passkey = ipass
286       &                   + (igadkey-1)*maxcube*Nr       &                   + (k-1)      *maxpass
287         &                   + (igadkey-1)*maxpass*Nr
288           IF (nipass .GT. maxpass) THEN           IF (nipass .GT. maxpass) THEN
289            STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'            STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'
290           ENDIF           ENDIF
# Line 347  C       a) needed in overlap only Line 348  C       a) needed in overlap only
348  C   and b) the overlap of myTile are not cube-face Edges  C   and b) the overlap of myTile are not cube-face Edges
349         IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN         IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN
350    
351  #ifndef ALLOW_AUTODIFF_TAMC  cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
352  C-     Internal exchange for calculations in X  C-     Internal exchange for calculations in X
353  #ifdef MULTIDIM_OLD_VERSION  #ifdef MULTIDIM_OLD_VERSION
354          IF ( useCubedSphereExchange ) THEN          IF ( useCubedSphereExchange ) THEN
# Line 355  C-     Internal exchange for calculation Line 356  C-     Internal exchange for calculation
356          IF ( useCubedSphereExchange .AND.          IF ( useCubedSphereExchange .AND.
357       &       ( overlapOnly .OR. ipass.EQ.1 ) ) THEN       &       ( overlapOnly .OR. ipass.EQ.1 ) ) THEN
358  #endif  #endif
359           CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid )           CALL FILL_CS_CORNER_TR_RL( .TRUE., .FALSE.,
360         &                              localTij, bi,bj, myThid )
361          ENDIF          ENDIF
362  #endif  cph-exch2#endif
363    
364  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
365  # ifndef DISABLE_MULTIDIM_ADVECTION  # ifndef DISABLE_MULTIDIM_ADVECTION
# Line 368  CADJ &     comlev1_bibj_k_gad_pass, key= Line 370  CADJ &     comlev1_bibj_k_gad_pass, key=
370    
371          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
372       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
373            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
374       I                           dTtracerLev(k),uTrans,uFld,localTij,       I                           dTtracerLev(k),uTrans,uFld,localTij,
375       O                           af, myThid )       O                           af, myThid )
376          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
377            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
378       I                              uTrans, uFld, maskLocW, localTij,       I                              uTrans, uFld, maskLocW, localTij,
379       O                              af, myThid )       O                              af, myThid )
380          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
381            CALL GAD_DST3_ADV_X(      bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_X(      bi,bj,k, .TRUE., dTtracerLev(k),
382       I                              uTrans, uFld, maskLocW, localTij,       I                              uTrans, uFld, maskLocW, localTij,
383       O                              af, myThid )       O                              af, myThid )
384          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
385            CALL GAD_DST3FL_ADV_X(    bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_X(    bi,bj,k, .TRUE., dTtracerLev(k),
386         I                              uTrans, uFld, maskLocW, localTij,
387         O                              af, myThid )
388    #ifndef ALLOW_AUTODIFF_TAMC
389            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
390              CALL GAD_OS7MP_ADV_X(     bi,bj,k, .TRUE., dTtracerLev(k),
391       I                              uTrans, uFld, maskLocW, localTij,       I                              uTrans, uFld, maskLocW, localTij,
392       O                              af, myThid )       O                              af, myThid )
393    #endif
394          ELSE          ELSE
395           STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'           STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
396          ENDIF          ENDIF
# Line 390  CADJ &     comlev1_bibj_k_gad_pass, key= Line 398  CADJ &     comlev1_bibj_k_gad_pass, key=
398  C-     Advective flux in X : done  C-     Advective flux in X : done
399         ENDIF         ENDIF
400    
401  #ifndef ALLOW_AUTODIFF_TAMC  cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
402  C-     Internal exchange for next calculations in Y  C-     Internal exchange for next calculations in Y
403         IF ( overlapOnly .AND. ipass.EQ.1 ) THEN         IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
404           CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid )           CALL FILL_CS_CORNER_TR_RL(.FALSE., .FALSE.,
405         &                              localTij, bi,bj, myThid )
406         ENDIF         ENDIF
407  #endif  cph-exch2#endif
408    
409  C-     Update the local tracer field where needed:  C-     Update the local tracer field where needed:
410    
# Line 507  C       a) needed in overlap only Line 516  C       a) needed in overlap only
516  C   and b) the overlap of myTile are not cube-face edges  C   and b) the overlap of myTile are not cube-face edges
517         IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN         IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
518    
519  #ifndef ALLOW_AUTODIFF_TAMC  cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
520  C-     Internal exchange for calculations in Y  C-     Internal exchange for calculations in Y
521  #ifdef MULTIDIM_OLD_VERSION  #ifdef MULTIDIM_OLD_VERSION
522          IF ( useCubedSphereExchange ) THEN          IF ( useCubedSphereExchange ) THEN
# Line 515  C-     Internal exchange for calculation Line 524  C-     Internal exchange for calculation
524          IF ( useCubedSphereExchange .AND.          IF ( useCubedSphereExchange .AND.
525       &       ( overlapOnly .OR. ipass.EQ.1 ) ) THEN       &       ( overlapOnly .OR. ipass.EQ.1 ) ) THEN
526  #endif  #endif
527           CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid )           CALL FILL_CS_CORNER_TR_RL(.FALSE., .FALSE.,
528         &                              localTij, bi,bj, myThid )
529          ENDIF          ENDIF
530  #endif  cph-exch2#endif
531    
532  C-     Advective flux in Y  C-     Advective flux in Y
533          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
# Line 535  CADJ &     comlev1_bibj_k_gad_pass, key= Line 545  CADJ &     comlev1_bibj_k_gad_pass, key=
545    
546          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
547       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
548            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
549       I                           dTtracerLev(k),vTrans,vFld,localTij,       I                           dTtracerLev(k),vTrans,vFld,localTij,
550       O                           af, myThid )       O                           af, myThid )
551          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
552            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
553       I                              vTrans, vFld, maskLocS, localTij,       I                              vTrans, vFld, maskLocS, localTij,
554       O                              af, myThid )       O                              af, myThid )
555          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
556            CALL GAD_DST3_ADV_Y(      bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_Y(      bi,bj,k, .TRUE., dTtracerLev(k),
557       I                              vTrans, vFld, maskLocS, localTij,       I                              vTrans, vFld, maskLocS, localTij,
558       O                              af, myThid )       O                              af, myThid )
559          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
560            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, .TRUE., dTtracerLev(k),
561         I                              vTrans, vFld, maskLocS, localTij,
562         O                              af, myThid )
563    #ifndef ALLOW_AUTODIFF_TAMC
564            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
565              CALL GAD_OS7MP_ADV_Y(     bi,bj,k, .TRUE., dTtracerLev(k),
566       I                              vTrans, vFld, maskLocS, localTij,       I                              vTrans, vFld, maskLocS, localTij,
567       O                              af, myThid )       O                              af, myThid )
568    #endif
569          ELSE          ELSE
570           STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'           STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
571          ENDIF          ENDIF
# Line 557  CADJ &     comlev1_bibj_k_gad_pass, key= Line 573  CADJ &     comlev1_bibj_k_gad_pass, key=
573  C-     Advective flux in Y : done  C-     Advective flux in Y : done
574         ENDIF         ENDIF
575    
576  #ifndef ALLOW_AUTODIFF_TAMC  cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
577  C-     Internal exchange for next calculations in X  C-     Internal exchange for next calculations in X
578         IF ( overlapOnly .AND. ipass.EQ.1 ) THEN         IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
579           CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid )           CALL FILL_CS_CORNER_TR_RL( .TRUE., .FALSE.,
580         &                              localTij, bi,bj, myThid )
581         ENDIF         ENDIF
582  #endif  cph-exch2#endif
583    
584  C-     Update the local tracer field where needed:  C-     Update the local tracer field where needed:
585    
# Line 697  C---+----1----+----2----+----3----+----4 Line 714  C---+----1----+----2----+----3----+----4
714  C--   Start of k loop for vertical flux  C--   Start of k loop for vertical flux
715         DO k=Nr,1,-1         DO k=Nr,1,-1
716  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
717           kkey = (igadkey-1)*Nr + k           kkey = (igadkey-1)*Nr + (Nr-k+1)
718  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
719  C--   kUp    Cycles through 1,2 to point to w-layer above  C--   kUp    Cycles through 1,2 to point to w-layer above
720  C--   kDown  Cycles through 2,1 to point to w-layer below  C--   kDown  Cycles through 2,1 to point to w-layer below
# Line 707  c       kp1=min(Nr,k+1) Line 724  c       kp1=min(Nr,k+1)
724          kp1Msk=1.          kp1Msk=1.
725          if (k.EQ.Nr) kp1Msk=0.          if (k.EQ.Nr) kp1Msk=0.
726    
727    #ifdef ALLOW_AUTODIFF_TAMC
728    CADJ STORE rtrans(:,:)  =
729    CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte
730    cphCADJ STORE wfld(:,:)  =
731    cphCADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte
732    #endif
733    
734  C-- Compute Vertical transport  C-- Compute Vertical transport
735  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
736  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
# Line 718  C- a hack to prevent Water-Vapor vert.tr Line 742  C- a hack to prevent Water-Vapor vert.tr
742  #endif  #endif
743    
744  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
745  CADJ STORE rtrans(:,:)  =  cphmultiCADJ STORE wfld(:,:)  =
746  CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte  cphmultiCADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte
 CADJ STORE wfld(:,:)  =  
 CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte  
747  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
748    
749  C- Surface interface :  C- Surface interface :
# Line 737  C- Surface interface : Line 759  C- Surface interface :
759          ELSE          ELSE
760    
761  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
762  CADJ STORE rtrans(:,:)  =  cphmultiCADJ STORE wfld(:,:)  =
763  CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte  cphmultiCADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte
 CADJ STORE wfld(:,:)  =  
 CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte  
764  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
765    
766  C- Interior interface :  C- Interior interface :
# Line 764  C--   Residual transp = Bolus transp + E Line 784  C--   Residual transp = Bolus transp + E
784  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
785    
786  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
787  CADJ STORE localTijk(:,:,k)    cphmultiCADJ STORE localTijk(:,:,k)  
788  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
789  CADJ STORE rTrans(:,:)    cphmultiCADJ STORE rTrans(:,:)  
790  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
791  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
792    
793  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
# Line 788  C-    Compute vertical advective flux in Line 808  C-    Compute vertical advective flux in
808             CALL GAD_DST3FL_ADV_R(    bi,bj,k, dTtracerLev(k),             CALL GAD_DST3FL_ADV_R(    bi,bj,k, dTtracerLev(k),
809       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localTijk,
810       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-Olx,1-Oly,kUp), myThid )
811    #ifndef ALLOW_AUTODIFF_TAMC
812             ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
813               CALL GAD_OS7MP_ADV_R(     bi,bj,k, dTtracerLev(k),
814         I                               rTrans, wFld, localTijk,
815         O                               fVerT(1-Olx,1-Oly,kUp), myThid )
816    #endif
817           ELSE           ELSE
818            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
819           ENDIF           ENDIF
# Line 796  C- end Surface/Interior if bloc Line 822  C- end Surface/Interior if bloc
822          ENDIF          ENDIF
823    
824  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
825  CADJ STORE rTrans(:,:)    cphmultiCADJ STORE rTrans(:,:)  
826  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
827  CADJ STORE rTranskp1(:,:)    cphmultiCADJ STORE rTranskp1(:,:)  
828    cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
829    cph --- following storing of fVerT is critical for correct
830    cph --- gradient with multiDimAdvection
831    cph --- Without it, kDown component is not properly recomputed
832    cph --- This is a TAF bug (and no warning available)
833    CADJ STORE fVerT(:,:,:)  
834  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
835  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
836    
Legend:
Removed from v.1.45  
changed lines
  Added in v.1.54
  ViewVC Help
Powered by ViewVC 1.1.22