/[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.17 by cnh, Wed Oct 28 03:11:36 1998 UTC revision 1.20 by adcroft, Tue May 18 18:01:12 1999 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    
3  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
4    
5  CStartOfInterFace  CStartOfInterFace
6        SUBROUTINE CALC_GT(        SUBROUTINE CALC_GT(
# Line 8  CStartOfInterFace Line 8  CStartOfInterFace
8       I           xA,yA,uTrans,vTrans,rTrans,maskup,maskC,       I           xA,yA,uTrans,vTrans,rTrans,maskup,maskC,
9       I           K13,K23,KappaRT,KapGM,       I           K13,K23,KappaRT,KapGM,
10       U           af,df,fZon,fMer,fVerT,       U           af,df,fZon,fMer,fVerT,
11       I           myThid )       I           myCurrentTime, myThid )
12  C     /==========================================================\  C     /==========================================================\
13  C     | SUBROUTINE CALC_GT                                       |  C     | SUBROUTINE CALC_GT                                       |
14  C     | o Calculate the temperature tendency terms.              |  C     | o Calculate the temperature tendency terms.              |
# Line 43  C     == GLobal variables == Line 43  C     == GLobal variables ==
43  #include "PARAMS.h"  #include "PARAMS.h"
44  #include "GRID.h"  #include "GRID.h"
45  #include "FFIELDS.h"  #include "FFIELDS.h"
46    #ifdef ALLOW_KPP
47    #include "KPPMIX.h"
48    #endif
49    
50    
51  C     == Routine arguments ==  C     == Routine arguments ==
52  C     fZon    - Work array for flux of temperature in the east-west  C     fZon    - Work array for flux of temperature in the east-west
# Line 82  C     myThid - Instance number for this Line 86  C     myThid - Instance number for this
86        INTEGER k,kUp,kDown,kM1        INTEGER k,kUp,kDown,kM1
87        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
88        INTEGER myThid        INTEGER myThid
89          _RL     myCurrentTime
90  CEndOfInterface  CEndOfInterface
91    
92  C     == Local variables ==  C     == Local variables ==
# Line 91  C     I, J, K - Loop counters Line 96  C     I, J, K - Loop counters
96        _RL afFacT, dfFacT        _RL afFacT, dfFacT
97        _RL dTdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dTdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RL dTdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dTdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99    #ifdef ALLOW_KPP
100          _RL hbl  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)    ! used by KPP mixing scheme
101          _RL frac (1-OLx:sNx+OLx,1-OLy:sNy+OLy)    ! used by KPP mixing scheme
102          integer jwtype                            ! index for Jerlov water type
103    #endif
104    
105        afFacT = 1. _d 0        afFacT = 1. _d 0
106        dfFacT = 1. _d 0        dfFacT = 1. _d 0
# Line 99  C     I, J, K - Loop counters Line 109  C     I, J, K - Loop counters
109  C---  Calculate advective and diffusive fluxes between cells.  C---  Calculate advective and diffusive fluxes between cells.
110    
111  C--   Zonal flux (fZon is at west face of "theta" cell)  C--   Zonal flux (fZon is at west face of "theta" cell)
112  C     Advective component of zonal flux  #ifdef INCLUDE_T_ADVECTION_CODE
113    C     o Advective component of zonal flux
114        DO j=jMin,jMax        DO j=jMin,jMax
115         DO i=iMin,iMax         DO i=iMin,iMax
116          af(i,j) =          af(i,j) =
117       &   uTrans(i,j)*(theta(i,j,k,bi,bj)+theta(i-1,j,k,bi,bj))*0.5 _d 0       &   uTrans(i,j)*(theta(i,j,k,bi,bj)+theta(i-1,j,k,bi,bj))*0.5 _d 0
118         ENDDO         ENDDO
119        ENDDO        ENDDO
120  C     Zonal tracer gradient  #endif /* INCLUDE_T_ADVECTION_CODE */
121    #ifdef INCLUDE_T_DIFFUSION_CODE
122    C     o Zonal tracer gradient
123        DO j=jMin,jMax        DO j=jMin,jMax
124         DO i=iMin,iMax         DO i=iMin,iMax
125          dTdx(i,j) = _recip_dxC(i,j,bi,bj)*          dTdx(i,j) = _recip_dxC(i,j,bi,bj)*
126       &  (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))
127         ENDDO         ENDDO
128        ENDDO        ENDDO
129  C     Diffusive component of zonal flux  C     o Diffusive component of zonal flux
130        DO j=jMin,jMax        DO j=jMin,jMax
131         DO i=iMin,iMax         DO i=iMin,iMax
132          df(i,j) = -(diffKhT+0.5*(KapGM(i,j)+KapGM(i-1,j)))*          df(i,j) = -(diffKhT+0.5*(KapGM(i,j)+KapGM(i-1,j)))*
133       &            xA(i,j)*dTdx(i,j)       &            xA(i,j)*dTdx(i,j)
134         ENDDO         ENDDO
135        ENDDO        ENDDO
136  C     Net zonal flux  #endif /* INCLUDE_T_DIFFUSION_CODE */
137    C     o Net zonal flux
138        DO j=jMin,jMax        DO j=jMin,jMax
139         DO i=iMin,iMax         DO i=iMin,iMax
140          fZon(i,j) = afFacT*af(i,j) + dfFacT*df(i,j)          fZon(i,j) = 0.
141         _ADT(&            + afFacT*af(i,j) )
142         _LPT(&            + dfFacT*df(i,j) )
143         ENDDO         ENDDO
144        ENDDO        ENDDO
145    
146  C--   Meridional flux (fMer is at south face of "theta" cell)  C--   Meridional flux (fMer is at south face of "theta" cell)
147  C     Advective component of meridional flux  #ifdef INCLUDE_T_ADVECTION_CODE
148    C     o Advective component of meridional flux
149        DO j=jMin,jMax        DO j=jMin,jMax
150         DO i=iMin,iMax         DO i=iMin,iMax
 C       Advective component of meridional flux  
151          af(i,j) =          af(i,j) =
152       &   vTrans(i,j)*(theta(i,j,k,bi,bj)+theta(i,j-1,k,bi,bj))*0.5 _d 0       &   vTrans(i,j)*(theta(i,j,k,bi,bj)+theta(i,j-1,k,bi,bj))*0.5 _d 0
153         ENDDO         ENDDO
154        ENDDO        ENDDO
155  C     Zonal tracer gradient  #endif /* INCLUDE_T_ADVECTION_CODE */
156    #ifdef INCLUDE_T_DIFFUSION_CODE
157    C     o Meridional tracer gradient
158        DO j=jMin,jMax        DO j=jMin,jMax
159         DO i=iMin,iMax         DO i=iMin,iMax
160          dTdy(i,j) = _recip_dyC(i,j,bi,bj)*          dTdy(i,j) = _recip_dyC(i,j,bi,bj)*
161       &  (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))
162         ENDDO         ENDDO
163        ENDDO        ENDDO
164  C     Diffusive component of meridional flux  C     o Diffusive component of meridional flux
165        DO j=jMin,jMax        DO j=jMin,jMax
166         DO i=iMin,iMax         DO i=iMin,iMax
167          df(i,j) = -(diffKhT+0.5*(KapGM(i,j)+KapGM(i,j-1)))*          df(i,j) = -(diffKhT+0.5*(KapGM(i,j)+KapGM(i,j-1)))*
168       &            yA(i,j)*dTdy(i,j)       &            yA(i,j)*dTdy(i,j)
169         ENDDO         ENDDO
170        ENDDO        ENDDO
171  C     Net meridional flux  #endif /* INCLUDE_T_DIFFUSION_CODE */
172    C     o Net meridional flux
173        DO j=jMin,jMax        DO j=jMin,jMax
174         DO i=iMin,iMax         DO i=iMin,iMax
175          fMer(i,j) = afFacT*af(i,j) + dfFacT*df(i,j)          fMer(i,j) = 0.
176         _ADT(&  + afFacT*af(i,j) )
177         _LPT(&  + dfFacT*df(i,j) )
178         ENDDO         ENDDO
179        ENDDO        ENDDO
180    
181  C--   Interpolate terms for Redi/GM scheme  #ifdef INCLUDE_T_DIFFUSION_CODE
182    C--   Terms that diffusion tensor projects onto z
183        DO j=jMin,jMax        DO j=jMin,jMax
184         DO i=iMin,iMax         DO i=iMin,iMax
185          dTdx(i,j) = 0.5*(          dTdx(i,j) = 0.5*(
# Line 194  C--   Interpolate terms for Redi/GM sche Line 216  C--   Interpolate terms for Redi/GM sche
216       &       )       &       )
217         ENDDO         ENDDO
218        ENDDO        ENDDO
219    #endif /* INCLUDE_T_DIFFUSION_CODE */
220    
221  C--   Vertical flux (fVerT) above  C--   Vertical flux ( fVerT(,,kUp) is at upper face of "theta" cell )
222  C     Advective component of vertical flux  #ifdef INCLUDE_T_ADVECTION_CODE
223    C     o Advective component of vertical flux
224  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
225  C     (this plays the role of the free-surface correction)  C     (this plays the role of the free-surface correction)
226        DO j=jMin,jMax        DO j=jMin,jMax
# Line 205  C     (this plays the role of the free-s Line 229  C     (this plays the role of the free-s
229       &   rTrans(i,j)*(theta(i,j,k,bi,bj)+theta(i,j,kM1,bi,bj))*0.5 _d 0       &   rTrans(i,j)*(theta(i,j,k,bi,bj)+theta(i,j,kM1,bi,bj))*0.5 _d 0
230         ENDDO         ENDDO
231        ENDDO        ENDDO
232  C     Diffusive component of vertical flux  #endif /* INCLUDE_T_ADVECTION_CODE */
233  C     Note: For K=1 then KM1=1 this gives a dT/dr = 0 upper  #ifdef INCLUDE_T_DIFFUSION_CODE
234    C     o Diffusive component of vertical flux
235    C     Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
236  C           boundary condition.  C           boundary condition.
237        DO j=jMin,jMax        DO j=jMin,jMax
238         DO i=iMin,iMax         DO i=iMin,iMax
# Line 226  C           boundary condition. Line 252  C           boundary condition.
252          ENDDO          ENDDO
253         ENDDO         ENDDO
254        ENDIF        ENDIF
255  C     Net vertical flux  #endif /* INCLUDE_T_DIFFUSION_CODE */
256    
257    #ifdef ALLOW_KPP
258          IF (usingKPPmixing) THEN
259    C--   Compute fraction of solar short-wave flux penetrating to
260    C     the bottom of the mixing layer
261           DO j=jMin,jMax
262            DO i=iMin,iMax
263             hbl(i,j) = KPPhbl(i,j,bi,bj)
264            ENDDO
265           ENDDO
266           j=(sNx+2*OLx)*(sNy+2*OLy)
267           jwtype = 3
268           CALL SWFRAC(
269         I     j, -1., hbl, jwtype,
270         O     frac )
271    
272    C     Add non local transport coefficient (ghat term) to right-hand-side
273    C     The nonlocal transport term is noNrero only for scalars in unstable
274    C     (convective) forcing conditions.
275    C     Note: -[Qnet * delZ(1) + Qsw * (1-frac) / KPPhbl] * 4000 * rho
276    C     is the total heat flux
277    C     penetrating the mixed layer from the surface in (deg C / s)
278           IF ( TOP_LAYER ) THEN
279            DO j=jMin,jMax
280             DO i=iMin,iMax
281              df(i,j) = df(i,j) + _rA(i,j,bi,bj) *
282         &           ( Qnet(i,j,bi,bj) * delZ(1) +
283         &           Qsw(i,j,bi,bj) * (1.-frac(i,j))
284         &           / KPPhbl(i,j,bi,bj) ) *
285         &           ( KappaRT(i,j,k) * KPPghat(i,j,k,  bi,bj) )
286             ENDDO
287            ENDDO
288           ELSE
289            DO j=jMin,jMax
290             DO i=iMin,iMax
291              df(i,j) = df(i,j) + _rA(i,j,bi,bj) *
292         &           ( Qnet(i,j,bi,bj) * delZ(1) +
293         &           Qsw(i,j,bi,bj)  * (1.-frac(i,j))
294         &           / KPPhbl(i,j,bi,bj) ) *
295         &           ( KappaRT(i,j,k)   * KPPghat(i,j,k,  bi,bj)
296         &           - KappaRT(i,j,k-1) * KPPghat(i,j,k-1,bi,bj) )
297             ENDDO
298            ENDDO
299           ENDIF
300          ENDIF
301    #endif /* ALLOW_KPP */
302    
303    C     o Net vertical flux
304        DO j=jMin,jMax        DO j=jMin,jMax
305         DO i=iMin,iMax         DO i=iMin,iMax
306          fVerT(i,j,kUp) = ( afFacT*af(i,j)+  dfFacT*df(i,j) )*maskUp(i,j)          fVerT(i,j,kUp) = 0.
307         _ADT(&  +afFacT*af(i,j)*maskUp(i,j) )
308         _LPT(&  +dfFacT*df(i,j)*maskUp(i,j) )
309         ENDDO         ENDDO
310        ENDDO        ENDDO
311    #ifdef INCLUDE_T_ADVECTION_CODE
312        IF ( TOP_LAYER ) THEN        IF ( TOP_LAYER ) THEN
313         DO j=jMin,jMax         DO j=jMin,jMax
314          DO i=iMin,iMax          DO i=iMin,iMax
# Line 239  C     Net vertical flux Line 316  C     Net vertical flux
316          ENDDO          ENDDO
317         ENDDO         ENDDO
318        ENDIF        ENDIF
319    #endif /* INCLUDE_T_ADVECTION_CODE */
320    
321  C--   Tendency is minus divergence of the fluxes.  C--   Tendency is minus divergence of the fluxes.
322  C     Note. Tendency terms will only be correct for range  C     Note. Tendency terms will only be correct for range
# Line 261  C           are not used. Line 339  C           are not used.
339         ENDDO         ENDDO
340        ENDDO        ENDDO
341    
342    #ifdef INCLUDE_T_FORCING_CODE
343  C--   External thermal forcing term(s)  C--   External thermal forcing term(s)
344  C     o Surface relaxation term        CALL EXTERNAL_FORCING_T(
345        IF ( TOP_LAYER ) THEN       I     iMin,iMax,jMin,jMax,bi,bj,k,
346         DO j=jMin,jMax       I     maskC,
347          DO i=iMin,iMax       I     myCurrentTime,myThid)
348           gT(i,j,k,bi,bj)=gT(i,j,k,bi,bj)  #endif /*  INCLUDE_T_FORCING_CODE */
349       &  +maskC(i,j)*(  
350       &   -lambdaThetaClimRelax*(theta(i,j,k,bi,bj)-SST(i,j,bi,bj))  #ifdef INCLUDE_LAT_CIRC_FFT_FILTER_CODE
351       &   -Qnet(i,j,bi,bj) )  C--   Zonal FFT filter of tendency
352          ENDDO        CALL FILTER_LATCIRCS_FFT_APPLY(
353         ENDDO       U     gT,
354        ENDIF       I     1, sNy, k, k, bi, bj, 1, myThid)
355    #endif /* INCLUDE_LAT_CIRC_FFT_FILTER_CODE */
356    
357    
358        RETURN        RETURN
359        END        END
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22