/[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.23 by gforget, Wed Aug 8 22:22:42 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 133  C     explicit/implicit timestepping wei Line 133  C     explicit/implicit timestepping wei
133    
134  C     Initialize local fields  C     Initialize local fields
135        DO k = 1, Nr        DO k = 1, Nr
136         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
137          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
138           KappaE(i,j,k)            = 0. _d 0           KappaE(i,j,k)            = 0. _d 0
139           TKEPrandtlNumber(i,j,k)  = 1. _d 0           TKEPrandtlNumber(i,j,k)  = 1. _d 0
140           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
# Line 147  C     Initialize local fields Line 147  C     Initialize local fields
147          ENDDO          ENDDO
148         ENDDO         ENDDO
149        ENDDO        ENDDO
150        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
151         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
         rhoK(i,j)          = 0. _d 0  
         rhoKm1(i,j)        = 0. _d 0  
152          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)
153          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
154          mxLength_Dn(i,j,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
# Line 160  C     Initialize local fields Line 158  C     Initialize local fields
158    
159  C     start k-loop  C     start k-loop
160        DO k = 2, Nr        DO k = 2, Nr
161         km1 = k-1  c      km1 = k-1
162  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 )  
163         DO j=jMin,jMax         DO j=jMin,jMax
164          DO i=iMin,iMax          DO i=iMin,iMax
165           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
166    
167  C     buoyancy frequency  C     buoyancy frequency
168           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(k)           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
169       &        * ( rhoKm1(i,j) - rhoK(i,j) )*maskC(i,j,k,bi,bj)       &                  * sigmaR(i,j,k)
170  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
171  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)
172  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 187  C     mixing length
187         ENDDO         ENDDO
188        ENDDO        ENDDO
189    
190    C- ensure mixing between first and second level
191          IF (mxlSurfFlag) THEN
192           DO j=jMin,jMax
193            DO i=iMin,iMax
194             GGL90mixingLength(i,j,2)=drF(1)
195            ENDDO
196           ENDDO
197          ENDIF
198    
199  C- Impose upper and lower bound for mixing length  C- Impose upper and lower bound for mixing length
200        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
201  C-  
202         DO k=2,Nr         DO k=2,Nr
203          DO j=jMin,jMax          DO j=jMin,jMax
204           DO i=iMin,iMax           DO i=iMin,iMax
# Line 224  C- Line 220  C-
220         ENDDO         ENDDO
221    
222        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
223  C-  
224         DO k=2,Nr         DO k=2,Nr
225          DO j=jMin,jMax          DO j=jMin,jMax
226           DO i=iMin,iMax           DO i=iMin,iMax
# Line 247  c         MaxLength=MAX(MaxLength,20. _d Line 243  c         MaxLength=MAX(MaxLength,20. _d
243         ENDDO         ENDDO
244    
245        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
246  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  
247         DO k=2,Nr         DO k=2,Nr
248          DO j=jMin,jMax          DO j=jMin,jMax
249           DO i=iMin,iMax           DO i=iMin,iMax
# Line 289  cgf Line 278  cgf
278         ENDDO         ENDDO
279    
280        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
281  C-  
282         DO k=2,Nr         DO k=2,Nr
283          DO j=jMin,jMax          DO j=jMin,jMax
284           DO i=iMin,iMax           DO i=iMin,iMax
# Line 348  c     ENDDO Line 337  c     ENDDO
337  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
338  C      do_fields_blocking_exchanges)  C      do_fields_blocking_exchanges)
339  C     common factors  C     common factors
340          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
341           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
342            xA(i,j) = _dyG(i,j,bi,bj)            xA(i,j) = _dyG(i,j,bi,bj)
343       &         *drF(k)*_hFacW(i,j,k,bi,bj)       &         *drF(k)*_hFacW(i,j,k,bi,bj)
344            yA(i,j) = _dxG(i,j,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)
# Line 358  C     common factors Line 347  C     common factors
347          ENDDO          ENDDO
348  C     Compute diffusive fluxes  C     Compute diffusive fluxes
349  C     ... across x-faces  C     ... across x-faces
350          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
351           dfx(1-Olx,j)=0. _d 0           dfx(1-OLx,j)=0. _d 0
352           DO i=1-Olx+1,sNx+Olx           DO i=1-OLx+1,sNx+OLx
353            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
354       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
355       &      *(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))
# Line 368  C     ... across x-faces Line 357  C     ... across x-faces
357           ENDDO           ENDDO
358          ENDDO          ENDDO
359  C     ... across y-faces  C     ... across y-faces
360          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
361           dfy(i,1-Oly)=0. _d 0           dfy(i,1-OLy)=0. _d 0
362          ENDDO          ENDDO
363          DO j=1-Oly+1,sNy+Oly          DO j=1-OLy+1,sNy+OLy
364           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
365            dfy(i,j) = -GGL90diffTKEh*yA(i,j)            dfy(i,j) = -GGL90diffTKEh*yA(i,j)
366       &      *_recip_dyC(i,j,bi,bj)       &      *_recip_dyC(i,j,bi,bj)
367       &      *(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 371  C     ... across y-faces
371           ENDDO           ENDDO
372          ENDDO          ENDDO
373  C     Compute divergence of fluxes  C     Compute divergence of fluxes
374          DO j=1-Oly,sNy+Oly-1          DO j=1-OLy,sNy+OLy-1
375           DO i=1-Olx,sNx+Olx-1           DO i=1-OLx,sNx+OLx-1
376            gTKE(i,j) =            gTKE(i,j) =
377       &    -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &    -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
378       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))

Legend:
Removed from v.1.20  
changed lines
  Added in v.1.23

  ViewVC Help
Powered by ViewVC 1.1.22