/[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.10 by dfer, Tue Oct 7 19:35:10 2008 UTC revision 1.11 by dfer, Fri Jan 30 02:23:56 2009 UTC
# Line 78  C     gTKE             - right hand side Line 78  C     gTKE             - right hand side
78        _RL     uStarSquare        _RL     uStarSquare
79        _RL     verticalShear        _RL     verticalShear
80        _RL     KappaM, KappaH        _RL     KappaM, KappaH
81        _RL     Nsquare  c     _RL     Nsquare
82          _RL     Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
83        _RL     deltaTggl90        _RL     deltaTggl90
84        _RL     SQRTTKE  c     _RL     SQRTTKE
85          _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
86        _RL     RiNumber        _RL     RiNumber
87        _RL     TKEdissipation        _RL     TKEdissipation
88        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
89          _RL     MaxLength
90        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
91        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
92        _RL         rMixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL         rMixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 104  C     dfx, dfy - diffusive flux across l Line 107  C     dfx, dfy - diffusive flux across l
107        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
110    #ifdef ALLOW_GGL90_SMOOTH
111          _RL p4, p8, p16, tmpdiffKrS
112          p4=0.25 _d 0
113          p8=0.125 _d 0
114          p16=0.0625 _d 0
115    #endif
116  CEOP  CEOP
117        iMin = 2-OLx        iMin = 2-OLx
118        iMax = sNx+OLx-1        iMax = sNx+OLx-1
# Line 136  C     Initialize local fields Line 145  C     Initialize local fields
145         DO I=1-Olx,sNx+Olx         DO I=1-Olx,sNx+Olx
146          rhoK    (I,J) = 0. _d 0          rhoK    (I,J) = 0. _d 0
147          rhoKm1  (I,J) = 0. _d 0          rhoKm1  (I,J) = 0. _d 0
148          totalDepth(I,J)   = 0. _d 0          totalDepth(I,J) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
         IF ( recip_Rcol(I,J,bi,bj) .NE. 0. _d 0 )  
      &       totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)  
149         ENDDO         ENDDO
150        ENDDO        ENDDO
151    
152  C     start k-loop  C     start k-loop
153        DO K = 2, Nr        DO K = 2, Nr
154         Km1 = K-1         Km1 = K-1
155         Kp1 = MIN(Nr,K+1)  c      Kp1 = MIN(Nr,K+1)
156         CALL FIND_RHO_2D(         CALL FIND_RHO_2D(
157       I      iMin, iMax, jMin, jMax, K,       I      iMin, iMax, jMin, jMax, K,
158       I      theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),       I      theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),
# Line 159  C     start k-loop Line 166  C     start k-loop
166       I      K, bi, bj, myThid )       I      K, bi, bj, myThid )
167         DO J=jMin,jMax         DO J=jMin,jMax
168          DO I=iMin,iMax          DO I=iMin,iMax
169           SQRTTKE=SQRT( GGL90TKE(I,J,K,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) )
170  C  C
171  C     buoyancy frequency  C     buoyancy frequency
172  C  C
173           Nsquare = - gravity*recip_rhoConst*recip_drC(K)           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)
174       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)
175    cC     vertical shear term (dU/dz)^2+(dV/dz)^2
176    c         tempU= .5 _d 0*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)
177    c     &                 -( uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )
178    c     &        *recip_drC(K)
179    c         tempV= .5 _d 0*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)
180    c     &                 -( vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )
181    c     &        *recip_drC(K)
182    c         verticalShear = tempU*tempU + tempV*tempV
183    c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
184    cC     compute Prandtl number (always greater than 0)
185    c         prTemp = 1. _d 0
186    c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
187    c         TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)
188    C     mixing length
189             GGL90mixingLength(I,J,K) = SQRTTWO *
190         &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
191            ENDDO
192           ENDDO
193          ENDDO
194    
195    C- Impose upper bound for mixing length (total depth)
196          IF ( mxlMaxFlag .EQ. 0 ) THEN
197           DO k=2,Nr
198            DO J=jMin,jMax
199             DO I=iMin,iMax
200              MaxLength=totalDepth(I,J)
201              GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
202         &        MaxLength)
203             ENDDO
204            ENDDO
205           ENDDO
206          ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
207           DO k=2,Nr
208            DO J=jMin,jMax
209             DO I=iMin,iMax
210              MaxLength=MIN(Ro_surf(I,J,bi,bj)-rF(k),rF(k)-R_low(I,J,bi,bj))
211    c         MaxLength=MAX(MaxLength,20. _d 0)
212              GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
213         &        MaxLength)
214             ENDDO
215            ENDDO
216           ENDDO
217          ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
218           DO k=2,Nr
219            DO J=jMin,jMax
220             DO I=iMin,iMax
221              GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
222         &        GGL90mixingLength(I,J,K-1)+drF(k-1))
223             ENDDO
224            ENDDO
225           ENDDO
226           DO J=jMin,jMax
227            DO I=iMin,iMax
228              GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),
229         &       GGL90mixingLengthMin+drF(Nr))
230            ENDDO
231           ENDDO
232           DO k=Nr-1,2,-1
233            DO J=jMin,jMax
234             DO I=iMin,iMax
235              GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
236         &        GGL90mixingLength(I,J,K+1)+drF(k))
237             ENDDO
238            ENDDO  
239           ENDDO
240          ELSE
241           STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing lenght limit)'
242          ENDIF
243    
244    C- Impose minimum mixing length (to avoid division by zero)
245          DO k=2,Nr
246           DO J=jMin,jMax
247            DO I=iMin,iMax
248             GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
249         &        GGL90mixingLengthMin)
250             rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)
251            ENDDO
252           ENDDO
253          ENDDO
254    
255          DO k=2,Nr
256           Km1 = K-1
257           DO J=jMin,jMax
258            DO I=iMin,iMax
259  C     vertical shear term (dU/dz)^2+(dV/dz)^2  C     vertical shear term (dU/dz)^2+(dV/dz)^2
260           tempu= .5 _d 0*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)           tempU= .5 _d 0*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)
261       &                 -( uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )       &                 -( uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )
262       &        *recip_drC(K)       &        *recip_drC(K)
263           tempv= .5 _d 0*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)           tempV= .5 _d 0*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)
264       &                 -( vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )       &                 -( vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )
265       &        *recip_drC(K)       &        *recip_drC(K)
266           verticalShear = tempU*tempU + tempV*tempV           verticalShear = tempU*tempU + tempV*tempV
267           RiNumber   = MAX(Nsquare,0. _d 0)/(verticalShear+GGL90eps)           RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
268  C     compute Prandtl number (always greater than 0)  C     compute Prandtl number (always greater than 0)
269           prTemp = 1. _d 0           prTemp = 1. _d 0
270           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
271           TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)           TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)
272  C     mixing length  
273           GGL90mixingLength(I,J,K) = SQRTTWO *  C     viscosity and diffusivity
274       &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE(i,j,k)
 C     impose upper bound for mixing length (total depth)  
          GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),  
      &        totalDepth(I,J))  
 C     impose minimum mixing length (to avoid division by zero)  
          GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),  
      &        GGL90mixingLengthMin)  
          rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)  
 C     viscosity of last timestep  
          KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE  
275           KappaH = KappaM/TKEPrandtlNumber(I,J,K)           KappaH = KappaM/TKEPrandtlNumber(I,J,K)
276    
277  C     Set a minium (= background) and maximum value  C     Set a minium (= background) and maximum value
# Line 198  C     Set a minium (= background) and ma Line 280  C     Set a minium (= background) and ma
280           KappaM = MIN(KappaM,GGL90viscMax)           KappaM = MIN(KappaM,GGL90viscMax)
281           KappaH = MIN(KappaH,GGL90diffMax)           KappaH = MIN(KappaH,GGL90diffMax)
282    
283  C     Mask land points and save  C     Mask land points and storage
          KappaE(I,J,K) = KappaM*GGL90alpha  
284           GGL90viscAr(I,J,K,bi,bj) = KappaM * maskC(I,J,K,bi,bj)           GGL90viscAr(I,J,K,bi,bj) = KappaM * maskC(I,J,K,bi,bj)
285           GGL90diffKr(I,J,K,bi,bj) = KappaH * maskC(I,J,K,bi,bj)           GGL90diffKr(I,J,K,bi,bj) = KappaH * maskC(I,J,K,bi,bj)
286             KappaE(I,J,K) = GGL90alpha * GGL90viscAr(I,J,K,bi,bj)
287    
288  C     dissipation term  C     dissipation term
289           TKEdissipation = ab05*GGL90ceps           TKEdissipation = ab05*GGL90ceps
290       &        *SQRTTKE*rMixingLength(I,J,K)       &        *SQRTTKE(i,j,k)*rMixingLength(I,J,K)
291       &        *GGL90TKE(I,J,K,bi,bj)       &        *GGL90TKE(I,J,K,bi,bj)
292  C     sum up contributions to form the right hand side  C     sum up contributions to form the right hand side
293           gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)           gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)
294       &        + deltaTggl90*(       &        + deltaTggl90*(
295       &        + KappaM*verticalShear       &        + KappaM*verticalShear
296       &        - KappaH*Nsquare       &        - KappaH*Nsquare(i,j,k)
 c    &        - KappaM*Nsquare/TKEPrandtlNumber(I,J,K)  
297       &        - TKEdissipation       &        - TKEdissipation
298       &        )       &        )
299          ENDDO          ENDDO
300         ENDDO         ENDDO
301        ENDDO        ENDDO
302    
303  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
304  C      do_fields_blocking_exchanges)  C      do_fields_blocking_exchanges)
305  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
# Line 285  C--   Lower diagonal Line 367  C--   Lower diagonal
367         ENDDO         ENDDO
368        ENDDO        ENDDO
369        DO k=2,Nr        DO k=2,Nr
370         km1=MAX(1,k-1)         km1=max(2,k-1)
371         DO j=jMin,jMax         DO j=jMin,jMax
372          DO i=iMin,iMax          DO i=iMin,iMax
373           a(i,j,k) = -deltaTggl90           a(i,j,k) = -deltaTggl90
374       &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)  c     &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)
375         &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
376       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
377       &        *recip_drC(k)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)
           IF (recip_hFacI(i,j,km1,bi,bj).EQ.0. _d 0) a(i,j,k)=0. _d 0  
378          ENDDO          ENDDO
379         ENDDO         ENDDO
380        ENDDO        ENDDO
# Line 300  C--   Upper diagonal Line 382  C--   Upper diagonal
382        DO j=jMin,jMax        DO j=jMin,jMax
383         DO i=iMin,iMax         DO i=iMin,iMax
384           c(i,j,1)  = 0. _d 0           c(i,j,1)  = 0. _d 0
          c(i,j,Nr) = 0. _d 0  
385         ENDDO         ENDDO
386        ENDDO        ENDDO
387  CML      DO k=1,Nr-1        DO k=2,Nr
       DO k=2,Nr-1  
        kp1=min(Nr,k+1)  
388         DO j=jMin,jMax         DO j=jMin,jMax
389          DO i=iMin,iMax          DO i=iMin,iMax
390              kp1=min(klowC(i,j,bi,bj),k+1)
391            c(i,j,k) = -deltaTggl90            c(i,j,k) = -deltaTggl90
392       &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)  c     &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)
393       &               *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
394       &        *recip_drC(k)       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
395            IF (recip_hFacI(i,j,kp1,bi,bj).EQ.0. _d 0) c(i,j,k)=0. _d 0       &        *recip_drC(k)*maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)
396          ENDDO          ENDDO
397         ENDDO         ENDDO
398        ENDDO        ENDDO
# Line 322  C--   Center diagonal Line 402  C--   Center diagonal
402          DO i=iMin,iMax          DO i=iMin,iMax
403            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)
404       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))
405       &        *rMixingLength(I,J,K)       &        *rMixingLength(I,J,K)*maskC(i,j,k,bi,bj)
406           ENDDO           ENDDO
407         ENDDO         ENDDO
408        ENDDO        ENDDO
# Line 344  C     Dirichlet surface boundary conditi Line 424  C     Dirichlet surface boundary conditi
424          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
425       &                     *maskC(I,J,kSurf,bi,bj)       &                     *maskC(I,J,kSurf,bi,bj)
426  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
427          kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)          kBottom   = MAX(kLowC(I,J,bi,bj),1)
428          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)
429       &       - GGL90TKEbottom*c(I,J,kBottom)       &                    - GGL90TKEbottom*c(I,J,kBottom)
430            c(I,J,kBottom) = 0. _d 0
431         ENDDO         ENDDO
432        ENDDO        ENDDO
433  C  C
# Line 368  C     impose minimum TKE to avoid numeri Line 449  C     impose minimum TKE to avoid numeri
449          ENDDO          ENDDO
450         ENDDO         ENDDO
451        ENDDO        ENDDO
452  C  
453  C     end of time step  C     end of time step
454  C     ===============================  C     ===============================
455  C     compute viscosity coefficients  
456  C  #ifdef ALLOW_GGL90_SMOOTH
457  c      DO K=2,Nr        DO K=1,Nr
458  c       DO J=jMin,jMax         DO J=jMin,jMax
459  c        DO I=iMin,iMax          DO I=iMin,iMax
460  C     Eq. (11), (18) and (21)           tmpdiffKrS=
461  c         KappaM = GGL90ck*GGL90mixingLength(I,J,K)*       &  (
462  c     &                  SQRT( GGL90TKE(I,J,K,bi,bj) )       &   p4 *  GGL90viscAr(i  ,j  ,k,bi,bj) * mskCor(i  ,j  ,bi,bj)
463  c         KappaH = KappaM/TKEPrandtlNumber(I,J,K)       &  +p8 *( GGL90viscAr(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
464  C     Set a minium (= background) value       &       + GGL90viscAr(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
465  c         KappaM = MAX(KappaM,viscAr)       &       + GGL90viscAr(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
466  c         KappaH = MAX(KappaH,diffKrNrT(k))       &       + GGL90viscAr(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
467  C     Set a maximum and mask land point       &  +p16*( GGL90viscAr(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)
468  c        GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)       &       + GGL90viscAr(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)
469  c    &        * maskC(I,J,K,bi,bj)       &       + GGL90viscAr(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
470  c        GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)       &       + GGL90viscAr(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))
471  c    &        * maskC(I,J,K,bi,bj)       &  )
472  c        ENDDO       & /(p4
473  c       ENDDO       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
474  C     end third k-loop       &       +       maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
475  c      ENDDO       &       +       maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
476         &       +       maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
477         &  +p16*(       maskC(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)
478         &       +       maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)
479         &       +       maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
480         &       +       maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))
481         &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
482         &   /TKEPrandtlNumber(i,j,k)
483             GGL90diffKrS(I,J,K,bi,bj)= MAX( tmpdiffKrS , diffKrNrT(k) )
484            ENDDO
485           ENDDO
486          ENDDO
487    #endif
488    
489  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
490        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22