/[MITgcm]/MITgcm/model/src/calc_gt.F
ViewVC logotype

Diff of /MITgcm/model/src/calc_gt.F

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

revision 1.31 by adcroft, Wed Feb 7 16:27:29 2001 UTC revision 1.32 by adcroft, Tue May 29 14:01:36 2001 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6    #define COSINEMETH_III
7    #undef  ISOTROPIC_COS_SCALING
8    
9  CStartOfInterFace  CStartOfInterFace
10        SUBROUTINE CALC_GT(        SUBROUTINE CALC_GT(
11       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
12       I           xA,yA,uTrans,vTrans,rTrans,maskup,maskC,       I           xA,yA,uTrans,vTrans,rTrans,maskUp,
13       I           KappaRT,       I           KappaRT,
14       U           fVerT,       U           fVerT,
15       I           myCurrentTime, myThid )       I           myCurrentTime, myThid )
# Line 51  C     == Routine arguments == Line 54  C     == Routine arguments ==
54  C     fVerT   - Flux of temperature (T) in the vertical  C     fVerT   - Flux of temperature (T) in the vertical
55  C               direction at the upper(U) and lower(D) faces of a cell.  C               direction at the upper(U) and lower(D) faces of a cell.
56  C     maskUp  - Land mask used to denote base of the domain.  C     maskUp  - Land mask used to denote base of the domain.
 C     maskC   - Land mask for theta cells (used in TOP_LAYER only)  
57  C     xA      - Tracer cell face area normal to X  C     xA      - Tracer cell face area normal to X
58  C     yA      - Tracer cell face area normal to X  C     yA      - Tracer cell face area normal to X
59  C     uTrans  - Zonal volume transport through cell face  C     uTrans  - Zonal volume transport through cell face
# Line 67  C     myThid - Instance number for this Line 69  C     myThid - Instance number for this
69        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
72        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
73        INTEGER k,kUp,kDown,kM1        INTEGER k,kUp,kDown,kM1
74        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
# Line 80  C     I, J, K - Loop counters Line 81  C     I, J, K - Loop counters
81        INTEGER i,j        INTEGER i,j
82        LOGICAL TOP_LAYER        LOGICAL TOP_LAYER
83        _RL afFacT, dfFacT        _RL afFacT, dfFacT
       _RL dTdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL dTdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
84        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 111  C---  Calculate advective and diffusive Line 110  C---  Calculate advective and diffusive
110  C     o Zonal tracer gradient  C     o Zonal tracer gradient
111        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
112         DO i=1-Olx+1,sNx+Olx         DO i=1-Olx+1,sNx+Olx
113          dTdx(i,j) = _recip_dxC(i,j,bi,bj)*          fZon(i,j) = _recip_dxC(i,j,bi,bj)*xA(i,j)
114       &  (theta(i,j,k,bi,bj)-theta(i-1,j,k,bi,bj))       &  *(theta(i,j,k,bi,bj)-theta(i-1,j,k,bi,bj))
115    #ifdef COSINEMETH_III
116         &   *sqCosFacU(j,bi,bj)
117    #endif
118         ENDDO         ENDDO
119        ENDDO        ENDDO
120  C     o Meridional tracer gradient  C     o Meridional tracer gradient
121        DO j=1-Oly+1,sNy+Oly        DO j=1-Oly+1,sNy+Oly
122         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
123          dTdy(i,j) = _recip_dyC(i,j,bi,bj)*          fMer(i,j) = _recip_dyC(i,j,bi,bj)*yA(i,j)
124       &  (theta(i,j,k,bi,bj)-theta(i,j-1,k,bi,bj))       &  *(theta(i,j,k,bi,bj)-theta(i,j-1,k,bi,bj))
125    #ifdef ISOTROPIC_COS_SCALING
126    #ifdef COSINEMETH_III
127         &   *sqCosFacV(j,bi,bj)
128    #endif
129    #endif
130         ENDDO         ENDDO
131        ENDDO        ENDDO
   
132  C--   del^2 of T, needed for bi-harmonic (del^4) term  C--   del^2 of T, needed for bi-harmonic (del^4) term
133        IF (diffK4T .NE. 0.) THEN        IF (diffK4T .NE. 0.) THEN
134         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
# Line 130  C--   del^2 of T, needed for bi-harmonic Line 136  C--   del^2 of T, needed for bi-harmonic
136           df4(i,j)= _recip_hFacC(i,j,k,bi,bj)           df4(i,j)= _recip_hFacC(i,j,k,bi,bj)
137       &             *recip_drF(k)/_rA(i,j,bi,bj)       &             *recip_drF(k)/_rA(i,j,bi,bj)
138       &            *(       &            *(
139       &             +( xA(i+1,j)*dTdx(i+1,j)-xA(i,j)*dTdx(i,j) )       &             +( fZon(i+1,j)-fZon(i,j) )
140       &             +( yA(i,j+1)*dTdy(i,j+1)-yA(i,j)*dTdy(i,j) )       &             +( fMer(i,j+1)-fMer(i,j) )
141       &             )       &             )
142          ENDDO          ENDDO
143         ENDDO         ENDDO
# Line 152  C     o Advective component of zonal flu Line 158  C     o Advective component of zonal flu
158  C     o Diffusive component of zonal flux  C     o Diffusive component of zonal flux
159        DO j=jMin,jMax        DO j=jMin,jMax
160         DO i=iMin,iMax         DO i=iMin,iMax
161          df(i,j) = -diffKhT*xA(i,j)*dTdx(i,j)          df(i,j) = -diffKhT*xA(i,j)*_recip_dxC(i,j,bi,bj)*
162         &  (theta(i,j,k,bi,bj)-theta(i-1,j,k,bi,bj))
163         &   *CosFacU(j,bi,bj)
164         ENDDO         ENDDO
165        ENDDO        ENDDO
166  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
# Line 168  C     o Add the bi-harmonic contribution Line 176  C     o Add the bi-harmonic contribution
176          DO i=iMin,iMax          DO i=iMin,iMax
177           df(i,j) = df(i,j) + xA(i,j)*           df(i,j) = df(i,j) + xA(i,j)*
178       &    diffK4T*(df4(i,j)-df4(i-1,j))*_recip_dxC(i,j,bi,bj)       &    diffK4T*(df4(i,j)-df4(i-1,j))*_recip_dxC(i,j,bi,bj)
179    #ifdef COSINEMETH_III
180         &   *sqCosFacU(j,bi,bj)
181    #else
182         &   *CosFacU(j,bi,bj)
183    #endif
184          ENDDO          ENDDO
185         ENDDO         ENDDO
186        ENDIF        ENDIF
# Line 195  C     o Advective component of meridiona Line 208  C     o Advective component of meridiona
208  C     o Diffusive component of meridional flux  C     o Diffusive component of meridional flux
209        DO j=jMin,jMax        DO j=jMin,jMax
210         DO i=iMin,iMax         DO i=iMin,iMax
211          df(i,j) = -diffKhT*yA(i,j)*dTdy(i,j)          df(i,j) = -diffKhT*yA(i,j)*_recip_dyC(i,j,bi,bj)*
212         &  (theta(i,j,k,bi,bj)-theta(i,j-1,k,bi,bj))
213    #ifdef ISOTROPIC_COS_SCALING
214         &   *CosFacV(j,bi,bj)
215    #endif
216         ENDDO         ENDDO
217        ENDDO        ENDDO
218  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
# Line 211  C     o Add the bi-harmonic contribution Line 228  C     o Add the bi-harmonic contribution
228          DO i=iMin,iMax          DO i=iMin,iMax
229           df(i,j) = df(i,j) + yA(i,j)*           df(i,j) = df(i,j) + yA(i,j)*
230       &    diffK4T*(df4(i,j)-df4(i,j-1))*_recip_dyC(i,j,bi,bj)       &    diffK4T*(df4(i,j)-df4(i,j-1))*_recip_dyC(i,j,bi,bj)
231    #ifdef ISOTROPIC_COS_SCALING
232    #ifdef COSINEMETH_III
233         &   *sqCosFacV(j,bi,bj)
234    #else
235         &   *CosFacV(j,bi,bj)
236    #endif
237    #endif
238          ENDDO          ENDDO
239         ENDDO         ENDDO
240        ENDIF        ENDIF
# Line 224  C     o Net meridional flux Line 248  C     o Net meridional flux
248         ENDDO         ENDDO
249        ENDDO        ENDDO
250    
251  #ifdef INCLUDE_T_DIFFUSION_CODE  C--   Vertical flux ( fVerT(,,kUp) is at upper face of "Tracer" cell )
 C--   Terms that diffusion tensor projects onto z  
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         dTdx(i,j) = 0.5*(  
      &   +0.5*(_maskW(i+1,j,k,bi,bj)  
      &         *_recip_dxC(i+1,j,bi,bj)*  
      &           (theta(i+1,j,k,bi,bj)-theta(i,j,k,bi,bj))  
      &        +_maskW(i,j,k,bi,bj)  
      &         *_recip_dxC(i,j,bi,bj)*  
      &           (theta(i,j,k,bi,bj)-theta(i-1,j,k,bi,bj)))  
      &   +0.5*(_maskW(i+1,j,km1,bi,bj)  
      &         *_recip_dxC(i+1,j,bi,bj)*  
      &           (theta(i+1,j,km1,bi,bj)-theta(i,j,km1,bi,bj))  
      &        +_maskW(i,j,km1,bi,bj)  
      &         *_recip_dxC(i,j,bi,bj)*  
      &           (theta(i,j,km1,bi,bj)-theta(i-1,j,km1,bi,bj)))  
      &       )  
        ENDDO  
       ENDDO  
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         dTdy(i,j) = 0.5*(  
      &   +0.5*(_maskS(i,j,k,bi,bj)  
      &         *_recip_dyC(i,j,bi,bj)*  
      &           (theta(i,j,k,bi,bj)-theta(i,j-1,k,bi,bj))  
      &        +_maskS(i,j+1,k,bi,bj)  
      &         *_recip_dyC(i,j+1,bi,bj)*  
      &           (theta(i,j+1,k,bi,bj)-theta(i,j,k,bi,bj)))  
      &   +0.5*(_maskS(i,j,km1,bi,bj)  
      &         *_recip_dyC(i,j,bi,bj)*  
      &           (theta(i,j,km1,bi,bj)-theta(i,j-1,km1,bi,bj))  
      &        +_maskS(i,j+1,km1,bi,bj)  
      &         *_recip_dyC(i,j+1,bi,bj)*  
      &           (theta(i,j+1,km1,bi,bj)-theta(i,j,km1,bi,bj)))  
      &       )  
        ENDDO  
       ENDDO  
 #endif /* INCLUDE_T_DIFFUSION_CODE */  
   
 C--   Vertical flux ( fVerT(,,kUp) is at upper face of "theta" cell )  
252  #ifdef INCLUDE_T_ADVECTION_CODE  #ifdef INCLUDE_T_ADVECTION_CODE
253  C     o Advective component of vertical flux  C     o Advective component of vertical flux : assume W_bottom=0 (mask)
254  C     Note: For K=1 then KM1=1 this gives a barZ(T) = T  C     Note: For K=1 then KM1=1 this gives a barZ(T) = T
255  C     (this plays the role of the free-surface correction)  C     (this plays the role of the free-surface correction)
256        DO j=jMin,jMax        IF ( rigidLid .AND. TOP_LAYER) THEN
257         DO i=iMin,iMax         DO j=jMin,jMax
258          af(i,j) =          DO i=iMin,iMax
259       &   rTrans(i,j)*(theta(i,j,k,bi,bj)+theta(i,j,kM1,bi,bj))*0.5 _d 0           af(i,j) = 0.
260            ENDDO
261         ENDDO         ENDDO
262        ENDDO        ELSEIF ( rigidLid ) THEN
263           DO j=jMin,jMax
264            DO i=iMin,iMax
265             af(i,j) = rTrans(i,j)*
266         &       (theta(i,j,k,bi,bj)+theta(i,j,kM1,bi,bj))*0.5 _d 0
267            ENDDO
268           ENDDO
269          ELSE
270    C-  include "free-surface correction" :
271           DO j=jMin,jMax
272            DO i=iMin,iMax
273             af(i,j) = rTrans(i,j)*(
274         &      maskC(i,j,kM1,bi,bj)*
275         &       (theta(i,j,k,bi,bj)+theta(i,j,kM1,bi,bj))*0.5 _d 0
276         &    +(maskC(i,j,k,bi,bj)-maskC(i,j,kM1,bi,bj))*
277         &        theta(i,j,k,bi,bj) )
278            ENDDO
279           ENDDO
280          ENDIF
281  #endif /* INCLUDE_T_ADVECTION_CODE */  #endif /* INCLUDE_T_ADVECTION_CODE */
282  #ifdef INCLUDE_T_DIFFUSION_CODE  #ifdef INCLUDE_T_DIFFUSION_CODE
283  C     o Diffusive component of vertical flux  C     o Diffusive component of vertical flux
# Line 310  C           boundary condition. Line 313  C           boundary condition.
313  C--   Add non local KPP transport term (ghat) to diffusive T flux.  C--   Add non local KPP transport term (ghat) to diffusive T flux.
314        IF (useKPP) CALL KPP_TRANSPORT_T(        IF (useKPP) CALL KPP_TRANSPORT_T(
315       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
316       I     maskC,KappaRT,       I     KappaRT,
317       U     df )       U     df )
318  #endif  #endif
319    
320  C     o Net vertical flux  C     o Net vertical flux
321        DO j=jMin,jMax        DO j=jMin,jMax
322         DO i=iMin,iMax         DO i=iMin,iMax
323    c       fVerT(i,j,kUp) = afFacT*af(i,j) + dfFacT*df(i,j)*maskUp(i,j)
324          fVerT(i,j,kUp) = 0.          fVerT(i,j,kUp) = 0.
325       & _ADT( +afFacT*af(i,j)*maskUp(i,j) )       & _ADT( +afFacT*af(i,j) )
326       & _LPT( +dfFacT*df(i,j)*maskUp(i,j) )       & _LPT( +dfFacT*df(i,j)*maskUp(i,j) )
327         ENDDO         ENDDO
328        ENDDO        ENDDO
 #ifdef INCLUDE_T_ADVECTION_CODE  
       IF ( TOP_LAYER ) THEN  
        DO j=jMin,jMax  
         DO i=iMin,iMax  
          fVerT(i,j,kUp) = afFacT*af(i,j)*freeSurfFac  
         ENDDO  
        ENDDO  
       ENDIF  
 #endif /* INCLUDE_T_ADVECTION_CODE */  
329    
330  C--   Tendency is minus divergence of the fluxes.  C--   Tendency is minus divergence of the fluxes.
331  C     Note. Tendency terms will only be correct for range  C     Note. Tendency terms will only be correct for range
# Line 340  C           they are not algorithmically Line 335  C           they are not algorithmically
335  C           are not used.  C           are not used.
336        DO j=jMin,jMax        DO j=jMin,jMax
337         DO i=iMin,iMax         DO i=iMin,iMax
 #define _recip_VolT1(i,j,k,bi,bj) _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)  
 #define _recip_VolT2(i,j,k,bi,bj) /_rA(i,j,bi,bj)  
338          gT(i,j,k,bi,bj)=          gT(i,j,k,bi,bj)=
339       &   -_recip_VolT1(i,j,k,bi,bj)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
340       &    _recip_VolT2(i,j,k,bi,bj)       &    *recip_rA(i,j,bi,bj)
341       &   *(       &   *(
342       &    +( fZon(i+1,j)-fZon(i,j) )       &    +( fZon(i+1,j)-fZon(i,j) )
343       &    +( fMer(i,j+1)-fMer(i,j) )       &    +( fMer(i,j+1)-fMer(i,j) )
# Line 357  C           are not used. Line 350  C           are not used.
350  C--   External thermal forcing term(s)  C--   External thermal forcing term(s)
351        CALL EXTERNAL_FORCING_T(        CALL EXTERNAL_FORCING_T(
352       I     iMin,iMax,jMin,jMax,bi,bj,k,       I     iMin,iMax,jMin,jMax,bi,bj,k,
      I     maskC,  
353       I     myCurrentTime,myThid)       I     myCurrentTime,myThid)
354  #endif /*  INCLUDE_T_FORCING_CODE */  #endif /*  INCLUDE_T_FORCING_CODE */
355    
Legend:
Removed from v.1.31  
changed lines
  Added in v.1.32
  ViewVC Help
Powered by ViewVC 1.1.22