/[MITgcm]/MITgcm/pkg/ggl90/ggl90_calc.F
ViewVC logotype

Diff of /MITgcm/pkg/ggl90/ggl90_calc.F

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

revision 1.20 by jmc, Thu Mar 15 15:23:22 2012 UTC revision 1.25 by mlosch, Mon Nov 5 13:11:11 2012 UTC
# Line 8  C !ROUTINE: GGL90_CALC Line 8  C !ROUTINE: GGL90_CALC
8    
9  C !INTERFACE: ======================================================  C !INTERFACE: ======================================================
10        SUBROUTINE GGL90_CALC(        SUBROUTINE GGL90_CALC(
11       I     bi, bj, myTime, myIter, myThid )       I                 bi, bj, sigmaR, myTime, myIter, myThid )
12    
13    
14  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
15  C     *==========================================================*  C     *==========================================================*
# Line 41  C !USES: =============================== Line 42  C !USES: ===============================
42    
43  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
44  C Routine arguments  C Routine arguments
45  C     bi, bj :: array indices on which to apply calculations  C     bi, bj :: Current tile indices
46    C     sigmaR :: Vertical gradient of iso-neutral density
47  C     myTime :: Current time in simulation  C     myTime :: Current time in simulation
48  C     myIter :: Current time-step number  C     myIter :: Current time-step number
49  C     myThid :: My Thread Id number  C     myThid :: My Thread Id number
50        INTEGER bi, bj        INTEGER bi, bj
51          _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
52        _RL     myTime        _RL     myTime
53        INTEGER myIter        INTEGER myIter
54        INTEGER myThid        INTEGER myThid
# Line 72  C     GGL90mixingLength:: mixing length Line 75  C     GGL90mixingLength:: mixing length
75  C         rMixingLength:: inverse of mixing length  C         rMixingLength:: inverse of mixing length
76  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)
77  C     TKEPrandtlNumber :: here, an empirical function of the Richardson number  C     TKEPrandtlNumber :: here, an empirical function of the Richardson number
 C     rhoK, rhoKm1     :: density at layer k and km1 (relative to k)  
78        INTEGER iMin ,iMax ,jMin ,jMax        INTEGER iMin ,iMax ,jMin ,jMax
79        INTEGER i, j, k, kp1, km1, kSurf, kBottom        INTEGER i, j, k, kp1, km1, kSurf, kBottom
80        _RL     explDissFac, implDissFac        _RL     explDissFac, implDissFac
# Line 93  c     _RL     SQRTTKE Line 95  c     _RL     SQRTTKE
95        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
96        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
98        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
100  C-    tri-diagonal matrix  C-    tri-diagonal matrix
# Line 103  C-    tri-diagonal matrix Line 103  C-    tri-diagonal matrix
103        _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
104        INTEGER errCode        INTEGER errCode
105  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
106    C     hFac     :: fractional thickness of W-cell
107  C     xA, yA   :: area of lateral faces  C     xA, yA   :: area of lateral faces
108  C     dfx, dfy :: diffusive flux across lateral faces  C     dfx, dfy :: diffusive flux across lateral faces
109  C     gTKE     :: right hand side of diffusion equation  C     gTKE     :: right hand side of diffusion equation
110          _RL     hFac
111        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 133  C     explicit/implicit timestepping wei Line 135  C     explicit/implicit timestepping wei
135    
136  C     Initialize local fields  C     Initialize local fields
137        DO k = 1, Nr        DO k = 1, Nr
138         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
139          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
140           KappaE(i,j,k)            = 0. _d 0           KappaE(i,j,k)            = 0. _d 0
141           TKEPrandtlNumber(i,j,k)  = 1. _d 0           TKEPrandtlNumber(i,j,k)  = 1. _d 0
142           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
# Line 147  C     Initialize local fields Line 149  C     Initialize local fields
149          ENDDO          ENDDO
150         ENDDO         ENDDO
151        ENDDO        ENDDO
152        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
153         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
         rhoK(i,j)          = 0. _d 0  
         rhoKm1(i,j)        = 0. _d 0  
154          totalDepth(i,j)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)          totalDepth(i,j)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
155          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
156          mxLength_Dn(i,j,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
# Line 160  C     Initialize local fields Line 160  C     Initialize local fields
160    
161  C     start k-loop  C     start k-loop
162        DO k = 2, Nr        DO k = 2, Nr
163         km1 = k-1  c      km1 = k-1
164  c      kp1 = MIN(Nr,k+1)  c      kp1 = MIN(Nr,k+1)
        CALL FIND_RHO_2D(  
      I      iMin, iMax, jMin, jMax, k,  
      I      theta(1-OLx,1-OLy,km1,bi,bj), salt(1-OLx,1-OLy,km1,bi,bj),  
      O      rhoKm1,  
      I      km1, bi, bj, myThid )  
   
        CALL FIND_RHO_2D(  
      I      iMin, iMax, jMin, jMax, k,  
      I      theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),  
      O      rhoK,  
      I      k, bi, bj, myThid )  
165         DO j=jMin,jMax         DO j=jMin,jMax
166          DO i=iMin,iMax          DO i=iMin,iMax
167           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
168    
169  C     buoyancy frequency  C     buoyancy frequency
170           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(k)           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
171       &        * ( rhoKm1(i,j) - rhoK(i,j) )*maskC(i,j,k,bi,bj)       &                  * sigmaR(i,j,k)
172  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
173  c         tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)  c         tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
174  c     &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )  c     &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
# Line 200  C     mixing length Line 189  C     mixing length
189         ENDDO         ENDDO
190        ENDDO        ENDDO
191    
192    C- ensure mixing between first and second level
193          IF (mxlSurfFlag) THEN
194           DO j=jMin,jMax
195            DO i=iMin,iMax
196             GGL90mixingLength(i,j,2)=drF(1)
197            ENDDO
198           ENDDO
199          ENDIF
200    
201  C- Impose upper and lower bound for mixing length  C- Impose upper and lower bound for mixing length
202        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
203  C-  
204         DO k=2,Nr         DO k=2,Nr
205          DO j=jMin,jMax          DO j=jMin,jMax
206           DO i=iMin,iMax           DO i=iMin,iMax
# Line 224  C- Line 222  C-
222         ENDDO         ENDDO
223    
224        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
225  C-  
226         DO k=2,Nr         DO k=2,Nr
227          DO j=jMin,jMax          DO j=jMin,jMax
228           DO i=iMin,iMax           DO i=iMin,iMax
# Line 247  c         MaxLength=MAX(MaxLength,20. _d Line 245  c         MaxLength=MAX(MaxLength,20. _d
245         ENDDO         ENDDO
246    
247        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
248  C-  
 cgf ensure mixing between first and second level  
 c      DO j=jMin,jMax  
 c        DO i=iMin,iMax  
 c         GGL90mixingLength(i,j,2)=drF(1)  
 c        ENDDO  
 c      ENDDO  
 cgf  
249         DO k=2,Nr         DO k=2,Nr
250          DO j=jMin,jMax          DO j=jMin,jMax
251           DO i=iMin,iMax           DO i=iMin,iMax
# Line 289  cgf Line 280  cgf
280         ENDDO         ENDDO
281    
282        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
283  C-  
284         DO k=2,Nr         DO k=2,Nr
285          DO j=jMin,jMax          DO j=jMin,jMax
286           DO i=iMin,iMax           DO i=iMin,iMax
# Line 344  c     ENDDO Line 335  c     ENDDO
335         km1 = k-1         km1 = k-1
336    
337  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
338         IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
339  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
340  C      do_fields_blocking_exchanges)  C      do_fields_blocking_exchanges)
341  C     common factors  C     common factors
342          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
343           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
344            xA(i,j) = _dyG(i,j,bi,bj)            xA(i,j) = _dyG(i,j,bi,bj)*drC(k)*
345       &         *drF(k)*_hFacW(i,j,k,bi,bj)       &                 (min(.5 _d 0,_hFacW(i,j,k-1,bi,bj) ) +
346            yA(i,j) = _dxG(i,j,bi,bj)       &                  min(.5 _d 0,_hFacW(i,j,k  ,bi,bj) ) )
347       &         *drF(k)*_hFacS(i,j,k,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)*drC(k)*
348         &                 (min(.5 _d 0,_hFacS(i,j,k-1,bi,bj) ) +
349         &                  min(.5 _d 0,_hFacS(i,j,k  ,bi,bj) ) )
350           ENDDO           ENDDO
351          ENDDO          ENDDO
352  C     Compute diffusive fluxes  C     Compute diffusive fluxes
353  C     ... across x-faces  C     ... across x-faces
354          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
355           dfx(1-Olx,j)=0. _d 0           dfx(1-OLx,j)=0. _d 0
356           DO i=1-Olx+1,sNx+Olx           DO i=1-OLx+1,sNx+OLx
357            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
358       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
359       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
360    #ifdef ISOTROPIC_COS_SCALING
361       &      *CosFacU(j,bi,bj)       &      *CosFacU(j,bi,bj)
362    #endif /* ISOTROPIC_COS_SCALING */
363           ENDDO           ENDDO
364          ENDDO          ENDDO
365  C     ... across y-faces  C     ... across y-faces
366          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
367           dfy(i,1-Oly)=0. _d 0           dfy(i,1-OLy)=0. _d 0
368          ENDDO          ENDDO
369          DO j=1-Oly+1,sNy+Oly          DO j=1-OLy+1,sNy+OLy
370           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
371            dfy(i,j) = -GGL90diffTKEh*yA(i,j)            dfy(i,j) = -GGL90diffTKEh*yA(i,j)
372       &      *_recip_dyC(i,j,bi,bj)       &      *_recip_dyC(i,j,bi,bj)
373       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
# Line 382  C     ... across y-faces Line 377  C     ... across y-faces
377           ENDDO           ENDDO
378          ENDDO          ENDDO
379  C     Compute divergence of fluxes  C     Compute divergence of fluxes
380          DO j=1-Oly,sNy+Oly-1          DO j=1-OLy,sNy+OLy-1
381           DO i=1-Olx,sNx+Olx-1           DO i=1-OLx,sNx+OLx-1
382            gTKE(i,j) =            hFac = min(.5 _d 0,_hFacC(i,j,k-1,bi,bj) ) +
383       &    -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &          min(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )
384       &         *( (dfx(i+1,j)-dfx(i,j))            gTKE(i,j) = 0.0
385       &           +(dfy(i,j+1)-dfy(i,j))            if ( hFac .ne. 0.0 )
386       &          )       &      gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)/hFac
387         &         *((dfx(i+1,j)-dfx(i,j))
388         &          +(dfy(i,j+1)-dfy(i,j)) )
389           ENDDO           ENDDO
390          ENDDO          ENDDO
391  C      end if GGL90diffTKEh .eq. 0.  C      end if GGL90diffTKEh .eq. 0.
Legend:
Removed from v.1.20  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.22