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

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

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

revision 1.42 by jmc, Sun Nov 19 22:04:09 2006 UTC revision 1.47 by jmc, Thu May 3 21:34:39 2007 UTC
# Line 210  C-    Advective flux in X Line 210  C-    Advective flux in X
210            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
211          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
212       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
213            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
214       I            dTtracerLev(k), uTrans, uFld, locABT,       I            dTtracerLev(k), uTrans, uFld, locABT,
215       O            af, myThid )       O            af, myThid )
216          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
217            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
218       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
219       O            af, myThid )       O            af, myThid )
220          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
# Line 222  C-    Advective flux in X Line 222  C-    Advective flux in X
222          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
223            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
224          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
225            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
226       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
227       O            af, myThid )       O            af, myThid )
228          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
# Line 230  C-    Advective flux in X Line 230  C-    Advective flux in X
230  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
231  cph IF inAdExact=.FALSE., we want to use DST3  cph IF inAdExact=.FALSE., we want to use DST3
232  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
233            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
234       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
235       O           af, myThid )       O           af, myThid )
236           ELSE           ELSE
237            CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
238       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
239       O           af, myThid )       O           af, myThid )
240           ENDIF           ENDIF
241    #ifndef ALLOW_AUTODIFF_TAMC
242            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
243              CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
244         I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
245         O            af, myThid )
246    #endif
247          ELSE          ELSE
248           STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'           STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
249          ENDIF          ENDIF
# Line 289  C *note* should update GMREDI_XTRANSPORT Line 295  C *note* should update GMREDI_XTRANSPORT
295          ENDIF          ENDIF
296        ENDIF        ENDIF
297  #endif  #endif
298    C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
299        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
300         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
301          fZon(i,j) = fZon(i,j) + df(i,j)          fZon(i,j) = fZon(i,j) + df(i,j)*rhoFacC(k)
302         ENDDO         ENDDO
303        ENDDO        ENDDO
304    
# Line 318  C-    Advective flux in Y Line 325  C-    Advective flux in Y
325            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
326          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
327       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
328            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
329       I            dTtracerLev(k), vTrans, vFld, locABT,       I            dTtracerLev(k), vTrans, vFld, locABT,
330       O            af, myThid )       O            af, myThid )
331          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
332            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
333       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
334       O            af, myThid )       O            af, myThid )
335          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
# Line 330  C-    Advective flux in Y Line 337  C-    Advective flux in Y
337          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
338            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
339          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
340            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
341       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
342       O            af, myThid )       O            af, myThid )
343          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
# Line 338  C-    Advective flux in Y Line 345  C-    Advective flux in Y
345  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
346  cph IF inAdExact=.FALSE., we want to use DST3  cph IF inAdExact=.FALSE., we want to use DST3
347  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
348            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
349       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
350       O           af, myThid )       O           af, myThid )
351           ELSE           ELSE
352            CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
353       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
354       O           af, myThid )       O           af, myThid )
355           ENDIF           ENDIF
356    #ifndef ALLOW_AUTODIFF_TAMC
357            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
358              CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
359         I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
360         O            af, myThid )
361    #endif
362          ELSE          ELSE
363            STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'            STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
364          ENDIF          ENDIF
# Line 397  C *note* should update GMREDI_YTRANSPORT Line 410  C *note* should update GMREDI_YTRANSPORT
410          ENDIF          ENDIF
411        ENDIF        ENDIF
412  #endif  #endif
413    C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
414        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
415         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
416          fMer(i,j) = fMer(i,j) + df(i,j)          fMer(i,j) = fMer(i,j) + df(i,j)*rhoFacC(k)
417         ENDDO         ENDDO
418        ENDDO        ENDDO
419    
# Line 456  cph with limiters in forward, but withou Line 470  cph with limiters in forward, but withou
470       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
471       O         af, myThid )       O         af, myThid )
472            ENDIF            ENDIF
473    #ifndef ALLOW_AUTODIFF_TAMC
474            ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
475               CALL GAD_OS7MP_ADV_R( bi,bj,k,
476         I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
477         O         af, myThid )
478    #endif
479          ELSE          ELSE
480            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
481          ENDIF          ENDIF
# Line 539  C-    Set non local KPP transport term ( Line 559  C-    Set non local KPP transport term (
559         ENDDO         ENDDO
560         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
561          CALL KPP_TRANSPORT_T(          CALL KPP_TRANSPORT_T(
562       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
563       O     df )       O           df,
564         I           myTime, myIter, myThid )
565         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
566          CALL KPP_TRANSPORT_S(          CALL KPP_TRANSPORT_S(
567       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
568       O     df )       O           df,
569         I           myTime, myIter, myThid )
570  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
571         ELSEIF (tracerIdentity .GE. GAD_TR1) THEN         ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
572          CALL KPP_TRANSPORT_PTR(          CALL KPP_TRANSPORT_PTR(
573       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
574       I     tracerIdentity-GAD_TR1+1,       I           tracerIdentity-GAD_TR1+1,
575       O     df )       O           df,
576         I           myTime, myIter, myThid )
577  #endif  #endif
578         ELSE         ELSE
579          PRINT*,'invalid tracer indentity: ', tracerIdentity          PRINT*,'invalid tracer indentity: ', tracerIdentity
# Line 558  C-    Set non local KPP transport term ( Line 581  C-    Set non local KPP transport term (
581         ENDIF         ENDIF
582         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
583          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
584           fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)           fVerT(i,j,kUp) = fVerT(i,j,kUp)
585         &                  + df(i,j)*maskUp(i,j)*rhoFacF(k)
586          ENDDO          ENDDO
587         ENDDO         ENDDO
588        ENDIF        ENDIF
589  #endif  #endif
590    
591  C--   Divergence of fluxes  C--   Divergence of fluxes
592    C     Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
593        DO j=1-Oly,sNy+Oly-1        DO j=1-Oly,sNy+Oly-1
594         DO i=1-Olx,sNx+Olx-1         DO i=1-Olx,sNx+Olx-1
595          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
596       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
597         &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
598       &   *( (fZon(i+1,j)-fZon(i,j))       &   *( (fZon(i+1,j)-fZon(i,j))
599       &     +(fMer(i,j+1)-fMer(i,j))       &     +(fMer(i,j+1)-fMer(i,j))
600       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
Legend:
Removed from v.1.42  
changed lines
  Added in v.1.47
  ViewVC Help
Powered by ViewVC 1.1.22