/[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.14 by cnh, Sat Aug 22 17:51:07 1998 UTC revision 1.28 by adcroft, Tue May 29 14:01:36 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6    #define COSINEMETH_III
7    #undef  ISOTROPIC_COS_SCALING
8    
9  CStartOfInterFace  CStartOfInterFace
10        SUBROUTINE CALC_GS(        SUBROUTINE CALC_GS(
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           K13,K23,KappaRS,KapGM,       I           KappaRS,
14       U           af,df,fZon,fMer,fVerS,       U           fVerS,
15       I           myThid )       I           myCurrentTime, myThid )
16  C     /==========================================================\  C     /==========================================================\
17  C     | SUBROUTINE CALC_GS                                       |  C     | SUBROUTINE CALC_GS                                       |
18  C     | o Calculate the salt tendency terms.                     |  C     | o Calculate the salt tendency terms.                     |
# Line 45  C     == GLobal variables == Line 49  C     == GLobal variables ==
49  #include "FFIELDS.h"  #include "FFIELDS.h"
50    
51  C     == Routine arguments ==  C     == Routine arguments ==
 C     fZon    - Work array for flux of temperature in the east-west  
 C               direction at the west face of a cell.  
 C     fMer    - Work array for flux of temperature in the north-south  
 C               direction at the south face of a cell.  
52  C     fVerS   - Flux of salt (S) in the vertical  C     fVerS   - Flux of salt (S) in the vertical
53  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.
54  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 salt cells (used in TOP_LAYER only)  
55  C     xA      - Tracer cell face area normal to X  C     xA      - Tracer cell face area normal to X
56  C     yA      - Tracer cell face area normal to X  C     yA      - Tracer cell face area normal to X
57  C     uTrans  - Zonal volume transport through cell face  C     uTrans  - Zonal volume transport through cell face
58  C     vTrans  - Meridional volume transport through cell face  C     vTrans  - Meridional volume transport through cell face
59  C     wTrans  - Vertical volume transport through cell face  C     rTrans  - Vertical volume transport through cell face
 C     af      - Advective flux component work array  
 C     df      - Diffusive flux component work array  
60  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
61  C                                      results will be set.  C                                      results will be set.
62  C     myThid - Instance number for this innvocation of CALC_GT  C     myThid - Instance number for this innvocation of CALC_GT
       _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
63        _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
64        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 72  C     myThid - Instance number for this Line 67  C     myThid - Instance number for this
67        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69        _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)  
       _RL K13   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL K23   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
70        _RL KappaRS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaRS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
71        INTEGER k,kUp,kDown,kM1        INTEGER k,kUp,kDown,kM1
72        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
73          _RL     myCurrentTime
74        INTEGER myThid        INTEGER myThid
75  CEndOfInterface  CEndOfInterface
76    
# Line 89  C     I, J, K - Loop counters Line 79  C     I, J, K - Loop counters
79        INTEGER i,j        INTEGER i,j
80        LOGICAL TOP_LAYER        LOGICAL TOP_LAYER
81        _RL afFacS, dfFacS        _RL afFacS, dfFacS
82        _RL dSdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83        _RL dSdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84          _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85          _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86          _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87    
88    #ifdef ALLOW_AUTODIFF_TAMC
89    C--   only the kUp part of fverS is set in this subroutine
90    C--   the kDown is still required
91          fVerS(1,1,kDown) = fVerS(1,1,kDown)
92    #endif
93          DO j=1-OLy,sNy+OLy
94           DO i=1-OLx,sNx+OLx
95            fZon(i,j)      = 0.0
96            fMer(i,j)      = 0.0
97            fVerS(i,j,kUp) = 0.0
98           ENDDO
99          ENDDO
100    
101        afFacS = 1. _d 0        afFacS = 1. _d 0
102        dfFacS = 1. _d 0        dfFacS = 1. _d 0
# Line 98  C     I, J, K - Loop counters Line 104  C     I, J, K - Loop counters
104    
105  C---  Calculate advective and diffusive fluxes between cells.  C---  Calculate advective and diffusive fluxes between cells.
106    
107    C     o Zonal tracer gradient
108          DO j=1-Oly,sNy+Oly
109           DO i=1-Olx+1,sNx+Olx
110            fZon(i,j) = _recip_dxC(i,j,bi,bj)*xA(i,j)
111         &   *(salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))
112    #ifdef COSINEMETH_III
113         &   *sqCosFacU(j,bi,bj)
114    #endif
115           ENDDO
116          ENDDO
117    C     o Meridional tracer gradient
118          DO j=1-Oly+1,sNy+Oly
119           DO i=1-Olx,sNx+Olx
120            fMer(i,j) = _recip_dyC(i,j,bi,bj)*yA(i,j)
121         &   *(salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
122    #ifdef ISOTROPIC_COS_SCALING
123    #ifdef COSINEMETH_III
124         &   *sqCosFacV(j,bi,bj)
125    #endif
126    #endif
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         &             +( fZon(i+1,j)-fZon(i,j) )
138         &             +( fMer(i,j+1)-fMer(i,j) )
139         &             )
140            ENDDO
141           ENDDO
142          ENDIF
143    
144  C--   Zonal flux (fZon is at west face of "salt" cell)  C--   Zonal flux (fZon is at west face of "salt" cell)
145  C     Advective component of zonal flux  C     Advective component of zonal flux
146        DO j=jMin,jMax        DO j=jMin,jMax
# Line 106  C     Advective component of zonal flux Line 149  C     Advective component of zonal flux
149       &   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
150         ENDDO         ENDDO
151        ENDDO        ENDDO
152  C     Zonal tracer gradient  C     o Diffusive component of zonal flux
153        DO j=jMin,jMax        DO j=jMin,jMax
154         DO i=iMin,iMax         DO i=iMin,iMax
155          dSdx(i,j) = _recip_dxC(i,j,bi,bj)*          df(i,j) = -diffKhS*xA(i,j)*_recip_dxC(i,j,bi,bj)*
156       &  (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))
157         &   *CosFacU(j,bi,bj)
158         ENDDO         ENDDO
159        ENDDO        ENDDO
160  C     Diffusive component of zonal flux  #ifdef ALLOW_GMREDI
161        DO j=jMin,jMax        IF (useGMRedi) CALL GMREDI_XTRANSPORT(
162         DO i=iMin,iMax       I     iMin,iMax,jMin,jMax,bi,bj,K,
163          df(i,j) = -(diffKhS+0.5*(KapGM(i,j)+KapGM(i-1,j)))*       I     xA,salt,
164       &            xA(i,j)*dSdx(i,j)       U     df,
165         I     myThid)
166    #endif
167    C     o Add the bi-harmonic contribution
168          IF (diffK4S .NE. 0.) THEN
169           DO j=jMin,jMax
170            DO i=iMin,iMax
171             df(i,j) = df(i,j) + xA(i,j)*
172         &    diffK4S*(df4(i,j)-df4(i-1,j))*_recip_dxC(i,j,bi,bj)
173    #ifdef COSINEMETH_III
174         &   *sqCosFacU(j,bi,bj)
175    #else
176         &   *CosFacU(j,bi,bj)
177    #endif
178            ENDDO
179         ENDDO         ENDDO
180        ENDDO        ENDIF
181  C     Net zonal flux  C     Net zonal flux
182        DO j=jMin,jMax        DO j=jMin,jMax
183         DO i=iMin,iMax         DO i=iMin,iMax
# Line 136  C       Advective component of meridiona Line 194  C       Advective component of meridiona
194       &   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
195         ENDDO         ENDDO
196        ENDDO        ENDDO
197  C     Zonal tracer gradient  C     Diffusive component of meridional flux
198        DO j=jMin,jMax        DO j=jMin,jMax
199         DO i=iMin,iMax         DO i=iMin,iMax
200          dSdy(i,j) = _recip_dyC(i,j,bi,bj)*          df(i,j) = -diffKhS*yA(i,j)*_recip_dyC(i,j,bi,bj)*
201       &  (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))
202         &   *CosFacV(j,bi,bj)
203         ENDDO         ENDDO
204        ENDDO        ENDDO
205  C     Diffusive component of meridional flux  #ifdef ALLOW_GMREDI
206        DO j=jMin,jMax        IF (useGMRedi) CALL GMREDI_YTRANSPORT(
207         DO i=iMin,iMax       I     iMin,iMax,jMin,jMax,bi,bj,K,
208          df(i,j) = -(diffKhS+0.5*(KapGM(i,j)+KapGM(i,j-1)))*       I     yA,salt,
209       &            yA(i,j)*dSdy(i,j)       U     df,
210         I     myThid)
211    #endif
212    C     o Add the bi-harmonic contribution
213          IF (diffK4S .NE. 0.) THEN
214           DO j=jMin,jMax
215            DO i=iMin,iMax
216             df(i,j) = df(i,j) + yA(i,j)*
217         &    diffK4S*(df4(i,j)-df4(i,j-1))*_recip_dyC(i,j,bi,bj)
218    #ifdef ISOTROPIC_COS_SCALING
219    #ifdef COSINEMETH_III
220         &   *sqCosFacV(j,bi,bj)
221    #else
222         &   *CosFacV(j,bi,bj)
223    #endif
224    #endif
225            ENDDO
226         ENDDO         ENDDO
227        ENDDO        ENDIF
228    
229  C     Net meridional flux  C     Net meridional flux
230        DO j=jMin,jMax        DO j=jMin,jMax
231         DO i=iMin,iMax         DO i=iMin,iMax
# Line 157  C     Net meridional flux Line 233  C     Net meridional flux
233         ENDDO         ENDDO
234        ENDDO        ENDDO
235    
236  C--   Interpolate terms for Redi/GM scheme  C--   Vertical flux ( fVerS(,,kUp) is at upper face of "Tracer" cell )
237        DO j=jMin,jMax  C     o Advective component of vertical flux : assume W_bottom=0 (mask)
238         DO i=iMin,iMax  C     Note: For K=1 then KM1=1 this gives a barZ(S) = S
239          dSdx(i,j) = 0.5*(  C     (this plays the role of the free-surface correction)
240       &   +0.5*(_maskW(i+1,j,k,bi,bj)*_recip_dxC(i+1,j,bi,bj)*        IF ( rigidLid .AND. TOP_LAYER) THEN
241       &           (salt(i+1,j,k,bi,bj)-salt(i,j,k,bi,bj))         DO j=jMin,jMax
242       &        +_maskW(i,j,k,bi,bj)*_recip_dxC(i,j,bi,bj)*          DO i=iMin,iMax
243       &           (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj)))           af(i,j) = 0.
244       &   +0.5*(_maskW(i+1,j,km1,bi,bj)*_recip_dxC(i+1,j,bi,bj)*          ENDDO
      &           (salt(i+1,j,km1,bi,bj)-salt(i,j,km1,bi,bj))  
      &        +_maskW(i,j,km1,bi,bj)*_recip_dxC(i,j,bi,bj)*  
      &           (salt(i,j,km1,bi,bj)-salt(i-1,j,km1,bi,bj)))  
      &       )  
245         ENDDO         ENDDO
246        ENDDO        ELSEIF ( rigidLid ) THEN
247        DO j=jMin,jMax         DO j=jMin,jMax
248         DO i=iMin,iMax          DO i=iMin,iMax
249          dSdy(i,j) = 0.5*(           af(i,j) = rTrans(i,j)*
250       &   +0.5*(_maskS(i,j,k,bi,bj)*_recip_dyC(i,j,bi,bj)*       &       (salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0
251       &           (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))          ENDDO
      &        +_maskS(i,j+1,k,bi,bj)*_recip_dyC(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)*_recip_dyC(i,j,bi,bj)*  
      &           (salt(i,j,km1,bi,bj)-salt(i,j-1,km1,bi,bj))  
      &        +_maskS(i,j+1,km1,bi,bj)*_recip_dyC(i,j+1,bi,bj)*  
      &           (salt(i,j+1,km1,bi,bj)-salt(i,j,km1,bi,bj)))  
      &       )  
252         ENDDO         ENDDO
253        ENDDO        ELSE
254           DO j=jMin,jMax
255  C--   Vertical flux (fVerS) above          DO i=iMin,iMax
256  C     Advective component of vertical flux           af(i,j) = rTrans(i,j)*(
257  C     Note: For K=1 then KM1=1 this gives a barZ(T) = T       &      maskC(i,j,kM1,bi,bj)*
258  C     (this plays the role of the free-surface correction)       &       (salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0
259        DO j=jMin,jMax       &    +(maskC(i,j,k,bi,bj)-maskC(i,j,kM1,bi,bj))*
260         DO i=iMin,iMax       &        salt(i,j,k,bi,bj) )
261          af(i,j) =          ENDDO
      &   rTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0  
262         ENDDO         ENDDO
263        ENDDO        ENDIF
264  C     Diffusive component of vertical flux  C     o Diffusive component of vertical flux
265  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
266  C           boundary condition.  C           boundary condition.
267        DO j=jMin,jMax        IF (implicitDiffusion) THEN
268         DO i=iMin,iMax         DO j=jMin,jMax
269          df(i,j) = _rA(i,j,bi,bj)*(          DO i=iMin,iMax
270       &   -KapGM(i,j)*K13(i,j,k)*dSdx(i,j)           df(i,j) = 0.
271       &   -KapGM(i,j)*K23(i,j,k)*dSdy(i,j)          ENDDO
      &   )  
272         ENDDO         ENDDO
273        ENDDO        ELSE
       IF (.NOT.implicitDiffusion) THEN  
274         DO j=jMin,jMax         DO j=jMin,jMax
275          DO i=iMin,iMax          DO i=iMin,iMax
276           df(i,j) = df(i,j) + _rA(i,j,bi,bj)*(           df(i,j) = - _rA(i,j,bi,bj)*(
277       &    -KappaRS(i,j,k)*recip_drC(k)       &    KappaRS(i,j,k)*recip_drC(k)
278       &    *(salt(i,j,kM1,bi,bj)-salt(i,j,k,bi,bj))*rkFac       &    *(salt(i,j,kM1,bi,bj)-salt(i,j,k,bi,bj))*rkFac
279       &    )       &    )
280          ENDDO          ENDDO
281         ENDDO         ENDDO
282        ENDIF        ENDIF
283    
284    #ifdef ALLOW_GMREDI
285          IF (useGMRedi) CALL GMREDI_RTRANSPORT(
286         I     iMin,iMax,jMin,jMax,bi,bj,K,
287         I     maskUp,salt,
288         U     df,
289         I     myThid)
290    #endif
291    
292    #ifdef ALLOW_KPP
293    C--   Add non-local KPP transport term (ghat) to diffusive salt flux.
294          IF (useKPP) CALL KPP_TRANSPORT_S(
295         I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
296         I     KappaRS,
297         U     df )
298    #endif
299    
300  C     Net vertical flux  C     Net vertical flux
301        DO j=jMin,jMax        DO j=jMin,jMax
302         DO i=iMin,iMax         DO i=iMin,iMax
303          fVerS(i,j,kUp) = ( afFacS*af(i,j)+  dfFacS*df(i,j) )*maskUp(i,j)          fVerS(i,j,kUp) = afFacS*af(i,j) + dfFacS*df(i,j)*maskUp(i,j)
304         ENDDO         ENDDO
305        ENDDO        ENDDO
       IF ( TOP_LAYER ) THEN  
        DO j=jMin,jMax  
         DO i=iMin,iMax  
          fVerS(i,j,kUp) = afFacS*af(i,j)*freeSurfFac  
         ENDDO  
        ENDDO  
       ENDIF  
306    
307  C--   Tendency is minus divergence of the fluxes.  C--   Tendency is minus divergence of the fluxes.
308  C     Note. Tendency terms will only be correct for range  C     Note. Tendency terms will only be correct for range
# Line 240  C           they are not algorithmically Line 312  C           they are not algorithmically
312  C           are not used.  C           are not used.
313        DO j=jMin,jMax        DO j=jMin,jMax
314         DO i=iMin,iMax         DO i=iMin,iMax
 #define _recip_VolS(i,j,k,bi,bj) _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)/_rA(i,j,bi,bj)  
315          gS(i,j,k,bi,bj)=          gS(i,j,k,bi,bj)=
316       &   -_recip_VolS(i,j,k,bi,bj)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
317         &    *recip_rA(i,j,bi,bj)
318       &   *(       &   *(
319       &    +( fZon(i+1,j)-fZon(i,j) )       &    +( fZon(i+1,j)-fZon(i,j) )
320       &    +( fMer(i,j+1)-fMer(i,j) )       &    +( fMer(i,j+1)-fMer(i,j) )
# Line 251  C           are not used. Line 323  C           are not used.
323         ENDDO         ENDDO
324        ENDDO        ENDDO
325    
326  C--   External P-E forcing term(s)  C--   External forcing term(s)
327  C     o Surface relaxation term        CALL EXTERNAL_FORCING_S(
328        IF ( TOP_LAYER ) THEN       I     iMin,iMax,jMin,jMax,bi,bj,k,
329         DO j=jMin,jMax       I     myCurrentTime,myThid)
         DO i=iMin,iMax  
          gS(i,j,k,bi,bj)=gS(i,j,k,bi,bj)  
      &   +maskC(i,j)*(  
      &   -lambdaSaltClimRelax*(salt(i,j,k,bi,bj)-SSS(i,j,bi,bj))  
      &   +EmPmR(i,j,bi,bj) )  
         ENDDO  
        ENDDO  
       ENDIF  
   
330    
331        RETURN        RETURN
332        END        END
Legend:
Removed from v.1.14  
changed lines
  Added in v.1.28
  ViewVC Help
Powered by ViewVC 1.1.22