/[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.27 by adcroft, Wed Feb 7 16:27:29 2001 UTC revision 1.28 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_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           KappaRS,       I           KappaRS,
14       U           fVerS,       U           fVerS,
15       I           myCurrentTime, myThid )       I           myCurrentTime, myThid )
# Line 49  C     == Routine arguments == Line 52  C     == Routine arguments ==
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
# Line 65  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)  
70        _RL KappaRS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaRS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
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
# Line 78  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
       _RL dSdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL dSdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
82        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83        _RL fZon  (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)        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 105  C--   the kDown is still required Line 104  C--   the kDown is still required
104    
105  C---  Calculate advective and diffusive fluxes between cells.  C---  Calculate advective and diffusive fluxes between cells.
106    
 #ifdef INCLUDE_T_DIFFUSION_CODE  
107  C     o Zonal tracer gradient  C     o Zonal tracer gradient
108        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
109         DO i=1-Olx+1,sNx+Olx         DO i=1-Olx+1,sNx+Olx
110          dSdx(i,j) = _recip_dxC(i,j,bi,bj)*          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))       &   *(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         ENDDO
116        ENDDO        ENDDO
117  C     o Meridional tracer gradient  C     o Meridional tracer gradient
118        DO j=1-Oly+1,sNy+Oly        DO j=1-Oly+1,sNy+Oly
119         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
120          dSdy(i,j) = _recip_dyC(i,j,bi,bj)*          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))       &   *(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         ENDDO
128        ENDDO        ENDDO
129    
# Line 128  C--   del^2 of S, needed for bi-harmonic Line 134  C--   del^2 of S, needed for bi-harmonic
134           df4(i,j)= _recip_hFacC(i,j,k,bi,bj)           df4(i,j)= _recip_hFacC(i,j,k,bi,bj)
135       &             *recip_drF(k)/_rA(i,j,bi,bj)       &             *recip_drF(k)/_rA(i,j,bi,bj)
136       &            *(       &            *(
137       &             +( xA(i+1,j)*dSdx(i+1,j)-xA(i,j)*dSdx(i,j) )       &             +( fZon(i+1,j)-fZon(i,j) )
138       &             +( yA(i,j+1)*dSdy(i,j+1)-yA(i,j)*dSdy(i,j) )       &             +( fMer(i,j+1)-fMer(i,j) )
139       &             )       &             )
140          ENDDO          ENDDO
141         ENDDO         ENDDO
142        ENDIF        ENDIF
 #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
# Line 147  C     Advective component of zonal flux Line 152  C     Advective component of zonal flux
152  C     o Diffusive component of zonal flux  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          df(i,j) = -diffKhS*xA(i,j)*dSdx(i,j)          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))
157         &   *CosFacU(j,bi,bj)
158         ENDDO         ENDDO
159        ENDDO        ENDDO
160  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
# Line 163  C     o Add the bi-harmonic contribution Line 170  C     o Add the bi-harmonic contribution
170          DO i=iMin,iMax          DO i=iMin,iMax
171           df(i,j) = df(i,j) + xA(i,j)*           df(i,j) = df(i,j) + xA(i,j)*
172       &    diffK4S*(df4(i,j)-df4(i-1,j))*_recip_dxC(i,j,bi,bj)       &    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          ENDDO
179         ENDDO         ENDDO
180        ENDIF        ENDIF
# Line 185  C       Advective component of meridiona Line 197  C       Advective component of meridiona
197  C     Diffusive component of meridional flux  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          df(i,j) = -diffKhS*yA(i,j)*dSdy(i,j)          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))
202         &   *CosFacV(j,bi,bj)
203         ENDDO         ENDDO
204        ENDDO        ENDDO
205  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
# Line 201  C     o Add the bi-harmonic contribution Line 215  C     o Add the bi-harmonic contribution
215          DO i=iMin,iMax          DO i=iMin,iMax
216           df(i,j) = df(i,j) + yA(i,j)*           df(i,j) = df(i,j) + yA(i,j)*
217       &    diffK4S*(df4(i,j)-df4(i,j-1))*_recip_dyC(i,j,bi,bj)       &    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          ENDDO
226         ENDDO         ENDDO
227        ENDIF        ENDIF
# Line 212  C     Net meridional flux Line 233  C     Net meridional flux
233         ENDDO         ENDDO
234        ENDDO        ENDDO
235    
236  C--   Vertical flux (fVerS) above  C--   Vertical flux ( fVerS(,,kUp) is at upper face of "Tracer" cell )
237  C     Advective component of vertical flux  C     o Advective component of vertical flux : assume W_bottom=0 (mask)
238  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(S) = S
239  C     (this plays the role of the free-surface correction)  C     (this plays the role of the free-surface correction)
240        DO j=jMin,jMax        IF ( rigidLid .AND. TOP_LAYER) THEN
241         DO i=iMin,iMax         DO j=jMin,jMax
242          af(i,j) =          DO i=iMin,iMax
243       &   rTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0           af(i,j) = 0.
244            ENDDO
245         ENDDO         ENDDO
246        ENDDO        ELSEIF ( rigidLid ) THEN
247           DO j=jMin,jMax
248            DO i=iMin,iMax
249             af(i,j) = rTrans(i,j)*
250         &       (salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0
251            ENDDO
252           ENDDO
253          ELSE
254           DO j=jMin,jMax
255            DO i=iMin,iMax
256             af(i,j) = rTrans(i,j)*(
257         &      maskC(i,j,kM1,bi,bj)*
258         &       (salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0
259         &    +(maskC(i,j,k,bi,bj)-maskC(i,j,kM1,bi,bj))*
260         &        salt(i,j,k,bi,bj) )
261            ENDDO
262           ENDDO
263          ENDIF
264  C     o Diffusive component of vertical flux  C     o Diffusive component of vertical flux
265  C     Note: For K=1 then KM1=1 and this gives a dS/dr = 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.
# Line 254  C           boundary condition. Line 293  C           boundary condition.
293  C--   Add non-local KPP transport term (ghat) to diffusive salt flux.  C--   Add non-local KPP transport term (ghat) to diffusive salt flux.
294        IF (useKPP) CALL KPP_TRANSPORT_S(        IF (useKPP) CALL KPP_TRANSPORT_S(
295       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
296       I     maskC,KappaRS,       I     KappaRS,
297       U     df )       U     df )
298  #endif  #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 278  C           i=iMin+1:iMax-1, j=jMin+1:jM Line 310  C           i=iMin+1:iMax-1, j=jMin+1:jM
310  C           will contain valid floating point numbers but  C           will contain valid floating point numbers but
311  C           they are not algorithmically correct. These points  C           they are not algorithmically correct. These points
312  C           are not used.  C           are not used.
313        DO j=jMin,jMax-1        DO j=jMin,jMax
314         DO i=iMin,iMax-1         DO i=iMin,iMax
 C     DO j=1-2,OLy+2  
 C      DO i=1-2,OLx+2  
 #define _recip_VolS1(i,j,k,bi,bj) _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)  
 #define _recip_VolS2(i,j,k,bi,bj) /_rA(i,j,bi,bj)  
315          gS(i,j,k,bi,bj)=          gS(i,j,k,bi,bj)=
316       &   -_recip_VolS1(i,j,k,bi,bj)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
317       &    _recip_VolS2(i,j,k,bi,bj)       &    *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 298  C      DO i=1-2,OLx+2 Line 326  C      DO i=1-2,OLx+2
326  C--   External forcing term(s)  C--   External forcing term(s)
327        CALL EXTERNAL_FORCING_S(        CALL EXTERNAL_FORCING_S(
328       I     iMin,iMax,jMin,jMax,bi,bj,k,       I     iMin,iMax,jMin,jMax,bi,bj,k,
      I     maskC,  
329       I     myCurrentTime,myThid)       I     myCurrentTime,myThid)
330    
331        RETURN        RETURN
Legend:
Removed from v.1.27  
changed lines
  Added in v.1.28
  ViewVC Help
Powered by ViewVC 1.1.22