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

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

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

revision 1.10 by adcroft, Mon Jun 22 15:26:25 1998 UTC revision 1.20 by heimbach, Fri Jun 9 14:26:30 2000 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_GS(        SUBROUTINE CALC_GS(
7       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
8       I           xA,yA,uTrans,vTrans,wTrans,maskup,maskC,       I           xA,yA,uTrans,vTrans,rTrans,maskup,maskC,
9       I           K13,K23,KappaZS,KapGM,       I           K13,K23,KappaRS,KapGM,
10       U           af,df,fZon,fMer,fVerS,       U           af,df,fZon,fMer,fVerS,
11       I           myThid )       I           myCurrentTime, myThid )
12  C     /==========================================================\  C     /==========================================================\
13  C     | SUBROUTINE CALC_GS                                       |  C     | SUBROUTINE CALC_GS                                       |
14  C     | o Calculate the salt tendency terms.                     |  C     | o Calculate the salt 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  C     == Routine arguments ==  C     == Routine arguments ==
51  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 57  C     xA      - Tracer cell face area no Line 60  C     xA      - Tracer cell face area no
60  C     yA      - Tracer cell face area normal to X  C     yA      - Tracer cell face area normal to X
61  C     uTrans  - Zonal volume transport through cell face  C     uTrans  - Zonal volume transport through cell face
62  C     vTrans  - Meridional volume transport through cell face  C     vTrans  - Meridional volume transport through cell face
63  C     wTrans  - Vertical volume transport through cell face  C     rTrans  - Vertical volume transport through cell face
64  C     af      - Advective flux component work array  C     af      - Advective flux component work array
65  C     df      - Diffusive flux component work array  C     df      - Diffusive flux component work array
66  C     bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation  C     bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation
# Line 70  C     myThid - Instance number for this Line 73  C     myThid - Instance number for this
73        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76        _RL wTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78        _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79        _RL K13   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)        _RL K13   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
80        _RL K23   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)        _RL K23   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
81        _RL KappaZS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)        _RL KappaRS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
82        _RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85        INTEGER k,kUp,kDown,kM1        INTEGER k,kUp,kDown,kM1
86        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
87          _RL     myCurrentTime
88        INTEGER myThid        INTEGER myThid
89  CEndOfInterface  CEndOfInterface
90    
# Line 91  C     I, J, K - Loop counters Line 95  C     I, J, K - Loop counters
95        _RL afFacS, dfFacS        _RL afFacS, dfFacS
96        _RL dSdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dSdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _RL dSdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dSdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98          _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99    
100    #ifdef ALLOW_AUTODIFF_TAMC
101    C--   only the kUp part of fverS is set in this subroutine
102    C--   the kDown is still required
103    
104          fVerS(1,1,kDown) = fVerS(1,1,kDown)
105          DO j=1-OLy,sNy+OLy
106           DO i=1-OLx,sNx+OLx
107            fZon(i,j)      = 0.0
108            fMer(i,j)      = 0.0
109            fVerS(i,j,kUp) = 0.0
110           ENDDO
111          ENDDO
112    #endif
113    
114        afFacS = 1. _d 0        afFacS = 1. _d 0
115        dfFacS = 1. _d 0        dfFacS = 1. _d 0
# Line 98  C     I, J, K - Loop counters Line 117  C     I, J, K - Loop counters
117    
118  C---  Calculate advective and diffusive fluxes between cells.  C---  Calculate advective and diffusive fluxes between cells.
119    
120    #ifdef INCLUDE_T_DIFFUSION_CODE
121    C     o Zonal tracer gradient
122          DO j=1-Oly,sNy+Oly
123           DO i=1-Olx+1,sNx+Olx
124            dSdx(i,j) = _recip_dxC(i,j,bi,bj)*
125         &  (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))
126           ENDDO
127          ENDDO
128    C     o Meridional tracer gradient
129          DO j=1-Oly+1,sNy+Oly
130           DO i=1-Olx,sNx+Olx
131            dSdy(i,j) = _recip_dyC(i,j,bi,bj)*
132         &  (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
133           ENDDO
134          ENDDO
135    
136    C--   del^2 of S, needed for bi-harmonic (del^4) term
137          IF (diffK4S .NE. 0.) THEN
138           DO j=1-Oly+1,sNy+Oly-1
139            DO i=1-Olx+1,sNx+Olx-1
140             df4(i,j)= _recip_hFacC(i,j,k,bi,bj)
141         &             *recip_drF(k)/_rA(i,j,bi,bj)
142         &            *(
143         &             +( xA(i+1,j)*dSdx(i+1,j)-xA(i,j)*dSdx(i,j) )
144         &             +( yA(i,j+1)*dSdy(i,j+1)-yA(i,j)*dSdy(i,j) )
145         &             )
146            ENDDO
147           ENDDO
148          ENDIF
149    #endif
150    
151  C--   Zonal flux (fZon is at west face of "salt" cell)  C--   Zonal flux (fZon is at west face of "salt" cell)
152  C     Advective component of zonal flux  C     Advective component of zonal flux
153        DO j=jMin,jMax        DO j=jMin,jMax
# Line 106  C     Advective component of zonal flux Line 156  C     Advective component of zonal flux
156       &   uTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i-1,j,k,bi,bj))*0.5 _d 0       &   uTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i-1,j,k,bi,bj))*0.5 _d 0
157         ENDDO         ENDDO
158        ENDDO        ENDDO
159  C     Zonal tracer gradient  C     o Diffusive component of zonal flux
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         dSdx(i,j) = _rdxC(i,j,bi,bj)*  
      &  (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))  
        ENDDO  
       ENDDO  
 C     Diffusive component of zonal flux  
160        DO j=jMin,jMax        DO j=jMin,jMax
161         DO i=iMin,iMax         DO i=iMin,iMax
162          df(i,j) = -(diffKhS+0.5*(KapGM(i,j)+KapGM(i-1,j)))*          df(i,j) = -(diffKhS+0.5*(KapGM(i,j)+KapGM(i-1,j)))*
163       &            xA(i,j)*dSdx(i,j)       &            xA(i,j)*dSdx(i,j)
164         ENDDO         ENDDO
165        ENDDO        ENDDO
166    C     o Add the bi-harmonic contribution
167          IF (diffK4S .NE. 0.) THEN
168           DO j=jMin,jMax
169            DO i=iMin,iMax
170             df(i,j) = df(i,j) + xA(i,j)*
171         &    diffK4S*(df4(i,j)-df4(i-1,j))*_recip_dxC(i,j,bi,bj)
172            ENDDO
173           ENDDO
174          ENDIF
175  C     Net zonal flux  C     Net zonal flux
176        DO j=jMin,jMax        DO j=jMin,jMax
177         DO i=iMin,iMax         DO i=iMin,iMax
# Line 136  C       Advective component of meridiona Line 188  C       Advective component of meridiona
188       &   vTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i,j-1,k,bi,bj))*0.5 _d 0       &   vTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i,j-1,k,bi,bj))*0.5 _d 0
189         ENDDO         ENDDO
190        ENDDO        ENDDO
 C     Zonal tracer gradient  
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         dSdy(i,j) = _rdyC(i,j,bi,bj)*  
      &  (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))  
        ENDDO  
       ENDDO  
191  C     Diffusive component of meridional flux  C     Diffusive component of meridional flux
192        DO j=jMin,jMax        DO j=jMin,jMax
193         DO i=iMin,iMax         DO i=iMin,iMax
# Line 150  C     Diffusive component of meridional Line 195  C     Diffusive component of meridional
195       &            yA(i,j)*dSdy(i,j)       &            yA(i,j)*dSdy(i,j)
196         ENDDO         ENDDO
197        ENDDO        ENDDO
198    C     o Add the bi-harmonic contribution
199          IF (diffK4S .NE. 0.) THEN
200           DO j=jMin,jMax
201            DO i=iMin,iMax
202             df(i,j) = df(i,j) + yA(i,j)*
203         &    diffK4S*(df4(i,j)-df4(i,j-1))*_recip_dyC(i,j,bi,bj)
204            ENDDO
205           ENDDO
206          ENDIF
207    
208  C     Net meridional flux  C     Net meridional flux
209        DO j=jMin,jMax        DO j=jMin,jMax
210         DO i=iMin,iMax         DO i=iMin,iMax
# Line 161  C--   Interpolate terms for Redi/GM sche Line 216  C--   Interpolate terms for Redi/GM sche
216        DO j=jMin,jMax        DO j=jMin,jMax
217         DO i=iMin,iMax         DO i=iMin,iMax
218          dSdx(i,j) = 0.5*(          dSdx(i,j) = 0.5*(
219       &   +0.5*(_maskW(i+1,j,k,bi,bj)*_rdxC(i+1,j,bi,bj)*       &   +0.5*(_maskW(i+1,j,k,bi,bj)
220         &         *_recip_dxC(i+1,j,bi,bj)*
221       &           (salt(i+1,j,k,bi,bj)-salt(i,j,k,bi,bj))       &           (salt(i+1,j,k,bi,bj)-salt(i,j,k,bi,bj))
222       &        +_maskW(i,j,k,bi,bj)*_rdxC(i,j,bi,bj)*       &        +_maskW(i,j,k,bi,bj)
223         &         *_recip_dxC(i,j,bi,bj)*
224       &           (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj)))       &           (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj)))
225       &   +0.5*(_maskW(i+1,j,km1,bi,bj)*_rdxC(i+1,j,bi,bj)*       &   +0.5*(_maskW(i+1,j,km1,bi,bj)
226         &         *_recip_dxC(i+1,j,bi,bj)*
227       &           (salt(i+1,j,km1,bi,bj)-salt(i,j,km1,bi,bj))       &           (salt(i+1,j,km1,bi,bj)-salt(i,j,km1,bi,bj))
228       &        +_maskW(i,j,km1,bi,bj)*_rdxC(i,j,bi,bj)*       &        +_maskW(i,j,km1,bi,bj)
229         &         *_recip_dxC(i,j,bi,bj)*
230       &           (salt(i,j,km1,bi,bj)-salt(i-1,j,km1,bi,bj)))       &           (salt(i,j,km1,bi,bj)-salt(i-1,j,km1,bi,bj)))
231       &       )       &       )
232         ENDDO         ENDDO
# Line 175  C--   Interpolate terms for Redi/GM sche Line 234  C--   Interpolate terms for Redi/GM sche
234        DO j=jMin,jMax        DO j=jMin,jMax
235         DO i=iMin,iMax         DO i=iMin,iMax
236          dSdy(i,j) = 0.5*(          dSdy(i,j) = 0.5*(
237       &   +0.5*(_maskS(i,j,k,bi,bj)*_rdyC(i,j,bi,bj)*       &   +0.5*(_maskS(i,j,k,bi,bj)
238         &         *_recip_dyC(i,j,bi,bj)*
239       &           (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))       &           (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
240       &        +_maskS(i,j+1,k,bi,bj)*_rdyC(i,j+1,bi,bj)*       &        +_maskS(i,j+1,k,bi,bj)
241         &         *_recip_dyC(i,j+1,bi,bj)*
242       &           (salt(i,j+1,k,bi,bj)-salt(i,j,k,bi,bj)))       &           (salt(i,j+1,k,bi,bj)-salt(i,j,k,bi,bj)))
243       &   +0.5*(_maskS(i,j,km1,bi,bj)*_rdyC(i,j,bi,bj)*       &   +0.5*(_maskS(i,j,km1,bi,bj)
244         &         *_recip_dyC(i,j,bi,bj)*
245       &           (salt(i,j,km1,bi,bj)-salt(i,j-1,km1,bi,bj))       &           (salt(i,j,km1,bi,bj)-salt(i,j-1,km1,bi,bj))
246       &        +_maskS(i,j+1,km1,bi,bj)*_rdyC(i,j+1,bi,bj)*       &        +_maskS(i,j+1,km1,bi,bj)
247         &         *_recip_dyC(i,j+1,bi,bj)*
248       &           (salt(i,j+1,km1,bi,bj)-salt(i,j,km1,bi,bj)))       &           (salt(i,j+1,km1,bi,bj)-salt(i,j,km1,bi,bj)))
249       &       )       &       )
250         ENDDO         ENDDO
# Line 194  C     (this plays the role of the free-s Line 257  C     (this plays the role of the free-s
257        DO j=jMin,jMax        DO j=jMin,jMax
258         DO i=iMin,iMax         DO i=iMin,iMax
259          af(i,j) =          af(i,j) =
260       &   wTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0       &   rTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0
261         ENDDO         ENDDO
262        ENDDO        ENDDO
263  C     Diffusive component of vertical flux  C     Diffusive component of vertical flux
# Line 202  C     Note: For K=1 then KM1=1 this give Line 265  C     Note: For K=1 then KM1=1 this give
265  C           boundary condition.  C           boundary condition.
266        DO j=jMin,jMax        DO j=jMin,jMax
267         DO i=iMin,iMax         DO i=iMin,iMax
268          df(i,j) = _zA(i,j,bi,bj)*(          df(i,j) = _rA(i,j,bi,bj)*(
269       &   -KapGM(i,j)*K13(i,j,k)*dSdx(i,j)       &   -KapGM(i,j)*K13(i,j,k)*dSdx(i,j)
270       &   -KapGM(i,j)*K23(i,j,k)*dSdy(i,j)       &   -KapGM(i,j)*K23(i,j,k)*dSdy(i,j)
271       &   )       &   )
# Line 211  C           boundary condition. Line 274  C           boundary condition.
274        IF (.NOT.implicitDiffusion) THEN        IF (.NOT.implicitDiffusion) THEN
275         DO j=jMin,jMax         DO j=jMin,jMax
276          DO i=iMin,iMax          DO i=iMin,iMax
277           df(i,j) = df(i,j) + _zA(i,j,bi,bj)*(           df(i,j) = df(i,j) + _rA(i,j,bi,bj)*(
278       &    -KappaZS(i,j,k)*rdzC(k)       &    -KappaRS(i,j,k)*recip_drC(k)
279       &    *(salt(i,j,kM1,bi,bj)-salt(i,j,k,bi,bj))       &    *(salt(i,j,kM1,bi,bj)-salt(i,j,k,bi,bj))*rkFac
280       &    )       &    )
281          ENDDO          ENDDO
282         ENDDO         ENDDO
283        ENDIF        ENDIF
284    #ifdef ALLOW_KPP
285          IF (usingKPPmixing) THEN
286    C--   Add non local transport coefficient (ghat term) to right-hand-side
287    C     The nonlocal transport term is noNrero only for scalars in unstable
288    C     (convective) forcing conditions.
289           IF ( TOP_LAYER ) THEN
290            DO j=jMin,jMax
291             DO i=iMin,iMax
292              df(i,j) = df(i,j) - _rA(i,j,bi,bj) *
293         &              EmPmR(i,j,bi,bj) * delZ(1) *
294         &              ( KappaRS(i,j,k)   * KPPghat(i,j,k,bi,bj)   )
295             ENDDO
296            ENDDO
297           ELSE
298            DO j=jMin,jMax
299             DO i=iMin,iMax
300              df(i,j) = df(i,j) - _rA(i,j,bi,bj) *
301         &              EmPmR(i,j,bi,bj) * delZ(1) *
302         &              ( KappaRS(i,j,k)   * KPPghat(i,j,k,bi,bj)
303         &              - KappaRS(i,j,k-1) * KPPghat(i,j,k-1,bi,bj) )
304             ENDDO
305            ENDDO
306           ENDIF
307          ENDIF
308    #endif /* ALLOW_KPP */
309    
310  C     Net vertical flux  C     Net vertical flux
311        DO j=jMin,jMax        DO j=jMin,jMax
312         DO i=iMin,iMax         DO i=iMin,iMax
# Line 240  C           they are not algorithmically Line 329  C           they are not algorithmically
329  C           are not used.  C           are not used.
330        DO j=jMin,jMax        DO j=jMin,jMax
331         DO i=iMin,iMax         DO i=iMin,iMax
332  C    &   -_rhFacC(i,j,k,bi,bj)*rdzF(k)*_rdxF(i,j,bi,bj)*_rdyF(i,j,bi,bj)  #define _recip_VolS1(i,j,k,bi,bj) _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
333  C    &   -_rhFacC(i,j,k,bi,bj)*rdzF(k)/_zA(i,j,bi,bj)  #define _recip_VolS2(i,j,k,bi,bj) /_rA(i,j,bi,bj)
 C #define _rVolS(i,j,k,bi,bj) _rhFacC(i,j,k,bi,bj)*rdzF(k)*_rdxF(i,j,bi,bj)*_rdyF(i,j,bi,bj)  
 #define _rVolS(i,j,k,bi,bj) _rhFacC(i,j,k,bi,bj)*rdzF(k)/_zA(i,j,bi,bj)  
334          gS(i,j,k,bi,bj)=          gS(i,j,k,bi,bj)=
335       &   -_rVolS(i,j,k,bi,bj)       &   -_recip_VolS1(i,j,k,bi,bj)
336         &    _recip_VolS2(i,j,k,bi,bj)
337       &   *(       &   *(
338       &    +( fZon(i+1,j)-fZon(i,j) )       &    +( fZon(i+1,j)-fZon(i,j) )
339       &    +( fMer(i,j+1)-fMer(i,j) )       &    +( fMer(i,j+1)-fMer(i,j) )
340       &    +( fVerS(i,j,kUp)-fVerS(i,j,kDown) )       &    +( fVerS(i,j,kUp)-fVerS(i,j,kDown) )*rkFac
341       &    )       &    )
342         ENDDO         ENDDO
343        ENDDO        ENDDO
344    
345  C--   External P-E forcing term(s)  C--   External forcing term(s)
346  C     o Surface relaxation term        CALL EXTERNAL_FORCING_S(
347        IF ( TOP_LAYER ) THEN       I     iMin,iMax,jMin,jMax,bi,bj,k,
348         DO j=jMin,jMax       I     maskC,
349          DO i=iMin,iMax       I     myCurrentTime,myThid)
350           gS(i,j,k,bi,bj)=gS(i,j,k,bi,bj)  
351       &   +maskC(i,j)*(  #ifdef INCLUDE_LAT_CIRC_FFT_FILTER_CODE
352       &   -lambdaSaltClimRelax*(salt(i,j,k,bi,bj)-SSS(i,j,bi,bj))  C--
353       &   -EmPpR(i,j,bi,bj) )        CALL FILTER_LATCIRCS_FFT_APPLY( gS, 1, sNy, k, k, bi, bj, 1, myThid)
354          ENDDO  #endif
        ENDDO  
       ENDIF  
   
355    
356        RETURN        RETURN
357        END        END
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22