/[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.21 by adcroft, Wed Jun 21 19:15:26 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           KappaRS,
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 57  C     xA      - Tracer cell face area no Line 57  C     xA      - Tracer cell face area no
57  C     yA      - Tracer cell face area normal to X  C     yA      - Tracer cell face area normal to X
58  C     uTrans  - Zonal volume transport through cell face  C     uTrans  - Zonal volume transport through cell face
59  C     vTrans  - Meridional volume transport through cell face  C     vTrans  - Meridional volume transport through cell face
60  C     wTrans  - Vertical volume transport through cell face  C     rTrans  - Vertical volume transport through cell face
61  C     af      - Advective flux component work array  C     af      - Advective flux component work array
62  C     df      - Diffusive flux component work array  C     df      - Diffusive flux component work array
63  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 70  C     myThid - Instance number for this
70        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73        _RL wTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75        _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76        _RL K13   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)        _RL KappaRS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL K23   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)  
       _RL KappaZS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)  
       _RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
77        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79        INTEGER k,kUp,kDown,kM1        INTEGER k,kUp,kDown,kM1
80        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
81          _RL     myCurrentTime
82        INTEGER myThid        INTEGER myThid
83  CEndOfInterface  CEndOfInterface
84    
# Line 91  C     I, J, K - Loop counters Line 89  C     I, J, K - Loop counters
89        _RL afFacS, dfFacS        _RL afFacS, dfFacS
90        _RL dSdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dSdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91        _RL dSdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dSdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92          _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93    
94    #ifdef ALLOW_AUTODIFF_TAMC
95    C--   only the kUp part of fverS is set in this subroutine
96    C--   the kDown is still required
97    
98          fVerS(1,1,kDown) = fVerS(1,1,kDown)
99          DO j=1-OLy,sNy+OLy
100           DO i=1-OLx,sNx+OLx
101            fZon(i,j)      = 0.0
102            fMer(i,j)      = 0.0
103            fVerS(i,j,kUp) = 0.0
104           ENDDO
105          ENDDO
106    #endif
107    
108        afFacS = 1. _d 0        afFacS = 1. _d 0
109        dfFacS = 1. _d 0        dfFacS = 1. _d 0
# Line 98  C     I, J, K - Loop counters Line 111  C     I, J, K - Loop counters
111    
112  C---  Calculate advective and diffusive fluxes between cells.  C---  Calculate advective and diffusive fluxes between cells.
113    
114    #ifdef INCLUDE_T_DIFFUSION_CODE
115    C     o Zonal tracer gradient
116          DO j=1-Oly,sNy+Oly
117           DO i=1-Olx+1,sNx+Olx
118            dSdx(i,j) = _recip_dxC(i,j,bi,bj)*
119         &  (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))
120           ENDDO
121          ENDDO
122    C     o Meridional tracer gradient
123          DO j=1-Oly+1,sNy+Oly
124           DO i=1-Olx,sNx+Olx
125            dSdy(i,j) = _recip_dyC(i,j,bi,bj)*
126         &  (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
127           ENDDO
128          ENDDO
129    
130    C--   del^2 of S, needed for bi-harmonic (del^4) term
131          IF (diffK4S .NE. 0.) THEN
132           DO j=1-Oly+1,sNy+Oly-1
133            DO i=1-Olx+1,sNx+Olx-1
134             df4(i,j)= _recip_hFacC(i,j,k,bi,bj)
135         &             *recip_drF(k)/_rA(i,j,bi,bj)
136         &            *(
137         &             +( xA(i+1,j)*dSdx(i+1,j)-xA(i,j)*dSdx(i,j) )
138         &             +( yA(i,j+1)*dSdy(i,j+1)-yA(i,j)*dSdy(i,j) )
139         &             )
140            ENDDO
141           ENDDO
142          ENDIF
143    #endif
144    
145  C--   Zonal flux (fZon is at west face of "salt" cell)  C--   Zonal flux (fZon is at west face of "salt" cell)
146  C     Advective component of zonal flux  C     Advective component of zonal flux
147        DO j=jMin,jMax        DO j=jMin,jMax
# Line 106  C     Advective component of zonal flux Line 150  C     Advective component of zonal flux
150       &   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
151         ENDDO         ENDDO
152        ENDDO        ENDDO
153  C     Zonal tracer gradient  C     o Diffusive component of zonal flux
154        DO j=jMin,jMax        DO j=jMin,jMax
155         DO i=iMin,iMax         DO i=iMin,iMax
156          dSdx(i,j) = _rdxC(i,j,bi,bj)*          df(i,j) = -diffKhS*xA(i,j)*dSdx(i,j)
      &  (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))  
157         ENDDO         ENDDO
158        ENDDO        ENDDO
159  C     Diffusive component of zonal flux  #ifdef ALLOW_GMREDI
160        DO j=jMin,jMax        IF (use_GMRedi) CALL GMREDI_XTRANSPORT(
161         DO i=iMin,iMax       I     iMin,iMax,jMin,jMax,bi,bj,K,
162          df(i,j) = -(diffKhS+0.5*(KapGM(i,j)+KapGM(i-1,j)))*       I     xA,salt,
163       &            xA(i,j)*dSdx(i,j)       U     df,
164         I     myThid)
165    #endif
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         ENDDO
174        ENDDO        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
194          df(i,j) = -(diffKhS+0.5*(KapGM(i,j)+KapGM(i,j-1)))*          df(i,j) = -diffKhS*yA(i,j)*dSdy(i,j)
      &            yA(i,j)*dSdy(i,j)  
195         ENDDO         ENDDO
196        ENDDO        ENDDO
197  C     Net meridional flux  #ifdef ALLOW_GMREDI
198        DO j=jMin,jMax        IF (use_GMRedi) CALL GMREDI_YTRANSPORT(
199         DO i=iMin,iMax       I     iMin,iMax,jMin,jMax,bi,bj,K,
200          fMer(i,j) = afFacS*af(i,j) + dfFacS*df(i,j)       I     yA,salt,
201         U     df,
202         I     myThid)
203    #endif
204    C     o Add the bi-harmonic contribution
205          IF (diffK4S .NE. 0.) THEN
206           DO j=jMin,jMax
207            DO i=iMin,iMax
208             df(i,j) = df(i,j) + yA(i,j)*
209         &    diffK4S*(df4(i,j)-df4(i,j-1))*_recip_dyC(i,j,bi,bj)
210            ENDDO
211         ENDDO         ENDDO
212        ENDDO        ENDIF
213    
214  C--   Interpolate terms for Redi/GM scheme  C     Net meridional flux
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         dSdx(i,j) = 0.5*(  
      &   +0.5*(_maskW(i+1,j,k,bi,bj)*_rdxC(i+1,j,bi,bj)*  
      &           (salt(i+1,j,k,bi,bj)-salt(i,j,k,bi,bj))  
      &        +_maskW(i,j,k,bi,bj)*_rdxC(i,j,bi,bj)*  
      &           (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj)))  
      &   +0.5*(_maskW(i+1,j,km1,bi,bj)*_rdxC(i+1,j,bi,bj)*  
      &           (salt(i+1,j,km1,bi,bj)-salt(i,j,km1,bi,bj))  
      &        +_maskW(i,j,km1,bi,bj)*_rdxC(i,j,bi,bj)*  
      &           (salt(i,j,km1,bi,bj)-salt(i-1,j,km1,bi,bj)))  
      &       )  
        ENDDO  
       ENDDO  
215        DO j=jMin,jMax        DO j=jMin,jMax
216         DO i=iMin,iMax         DO i=iMin,iMax
217          dSdy(i,j) = 0.5*(          fMer(i,j) = afFacS*af(i,j) + dfFacS*df(i,j)
      &   +0.5*(_maskS(i,j,k,bi,bj)*_rdyC(i,j,bi,bj)*  
      &           (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))  
      &        +_maskS(i,j+1,k,bi,bj)*_rdyC(i,j+1,bi,bj)*  
      &           (salt(i,j+1,k,bi,bj)-salt(i,j,k,bi,bj)))  
      &   +0.5*(_maskS(i,j,km1,bi,bj)*_rdyC(i,j,bi,bj)*  
      &           (salt(i,j,km1,bi,bj)-salt(i,j-1,km1,bi,bj))  
      &        +_maskS(i,j+1,km1,bi,bj)*_rdyC(i,j+1,bi,bj)*  
      &           (salt(i,j+1,km1,bi,bj)-salt(i,j,km1,bi,bj)))  
      &       )  
218         ENDDO         ENDDO
219        ENDDO        ENDDO
220    
# Line 194  C     (this plays the role of the free-s Line 225  C     (this plays the role of the free-s
225        DO j=jMin,jMax        DO j=jMin,jMax
226         DO i=iMin,iMax         DO i=iMin,iMax
227          af(i,j) =          af(i,j) =
228       &   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
229         ENDDO         ENDDO
230        ENDDO        ENDDO
231  C     Diffusive component of vertical flux  C     o Diffusive component of vertical flux
232  C     Note: For K=1 then KM1=1 this gives a dS/dz = 0 upper  C     Note: For K=1 then KM1=1 and this gives a dS/dr = 0 upper
233  C           boundary condition.  C           boundary condition.
234        DO j=jMin,jMax        IF (implicitDiffusion) THEN
235         DO i=iMin,iMax         DO j=jMin,jMax
236          df(i,j) = _zA(i,j,bi,bj)*(          DO i=iMin,iMax
237       &   -KapGM(i,j)*K13(i,j,k)*dSdx(i,j)           df(i,j) = 0.
238       &   -KapGM(i,j)*K23(i,j,k)*dSdy(i,j)          ENDDO
      &   )  
239         ENDDO         ENDDO
240        ENDDO        ELSE
       IF (.NOT.implicitDiffusion) THEN  
241         DO j=jMin,jMax         DO j=jMin,jMax
242          DO i=iMin,iMax          DO i=iMin,iMax
243           df(i,j) = df(i,j) + _zA(i,j,bi,bj)*(           df(i,j) = - _rA(i,j,bi,bj)*(
244       &    -KappaZS(i,j,k)*rdzC(k)       &    KappaRS(i,j,k)*recip_drC(k)
245       &    *(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
246       &    )       &    )
247          ENDDO          ENDDO
248         ENDDO         ENDDO
249        ENDIF        ENDIF
250    
251    #ifdef ALLOW_GMREDI
252          IF (use_GMRedi) CALL GMREDI_RTRANSPORT(
253         I     iMin,iMax,jMin,jMax,bi,bj,K,
254         I     maskUp,salt,
255         U     df,
256         I     myThid)
257    #endif
258    
259    #ifdef ALLOW_KPP
260    C--   Add non-local KPP transport term (ghat) to diffusive salt flux.
261          IF (use_KPPmixing) CALL KPP_TRANSPORT_S(
262         I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
263         I     maskC,KappaRS,
264         U     df )
265    #endif
266    
267  C     Net vertical flux  C     Net vertical flux
268        DO j=jMin,jMax        DO j=jMin,jMax
269         DO i=iMin,iMax         DO i=iMin,iMax
# Line 240  C           they are not algorithmically Line 286  C           they are not algorithmically
286  C           are not used.  C           are not used.
287        DO j=jMin,jMax        DO j=jMin,jMax
288         DO i=iMin,iMax         DO i=iMin,iMax
289  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)
290  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)  
291          gS(i,j,k,bi,bj)=          gS(i,j,k,bi,bj)=
292       &   -_rVolS(i,j,k,bi,bj)       &   -_recip_VolS1(i,j,k,bi,bj)
293         &    _recip_VolS2(i,j,k,bi,bj)
294       &   *(       &   *(
295       &    +( fZon(i+1,j)-fZon(i,j) )       &    +( fZon(i+1,j)-fZon(i,j) )
296       &    +( fMer(i,j+1)-fMer(i,j) )       &    +( fMer(i,j+1)-fMer(i,j) )
297       &    +( fVerS(i,j,kUp)-fVerS(i,j,kDown) )       &    +( fVerS(i,j,kUp)-fVerS(i,j,kDown) )*rkFac
298       &    )       &    )
299         ENDDO         ENDDO
300        ENDDO        ENDDO
301    
302  C--   External P-E forcing term(s)  C--   External forcing term(s)
303  C     o Surface relaxation term        CALL EXTERNAL_FORCING_S(
304        IF ( TOP_LAYER ) THEN       I     iMin,iMax,jMin,jMax,bi,bj,k,
305         DO j=jMin,jMax       I     maskC,
306          DO i=iMin,iMax       I     myCurrentTime,myThid)
307           gS(i,j,k,bi,bj)=gS(i,j,k,bi,bj)  
308       &   +maskC(i,j)*(  #ifdef INCLUDE_LAT_CIRC_FFT_FILTER_CODE
309       &   -lambdaSaltClimRelax*(salt(i,j,k,bi,bj)-SSS(i,j,bi,bj))  C--
310       &   -EmPpR(i,j,bi,bj) )        CALL FILTER_LATCIRCS_FFT_APPLY( gS, 1, sNy, k, k, bi, bj, 1, myThid)
311          ENDDO  #endif
        ENDDO  
       ENDIF  
   
312    
313        RETURN        RETURN
314        END        END
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.21
  ViewVC Help
Powered by ViewVC 1.1.22