/[MITgcm]/MITgcm/pkg/gmredi/gmredi_calc_tensor.F
ViewVC logotype

Diff of /MITgcm/pkg/gmredi/gmredi_calc_tensor.F

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

revision 1.33 by jmc, Tue Oct 21 22:10:55 2008 UTC revision 1.34 by jmc, Mon Oct 27 22:13:12 2008 UTC
# Line 88  C     == Local variables == Line 88  C     == Local variables ==
88  #ifdef OLD_VISBECK_CALC  #ifdef OLD_VISBECK_CALC
89        _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
90  #else  #else
91        _RL dSigmaH        _RL dSigmaH, dSigmaR
92        _RL Sloc, M2loc        _RL Sloc, M2loc
93  #endif  #endif
94          _RL recipMaxSlope
95        _RL deltaH, integrDepth        _RL deltaH, integrDepth
96        _RL N2loc, SNloc        _RL N2loc, SNloc
97  #endif  #endif /* GM_VISBECK_VARIABLE_K */
98    
99  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
100        LOGICAL  doDiagRediFlx        LOGICAL  doDiagRediFlx
# Line 129  C---+----1----+----2----+----3----+----4 Line 130  C---+----1----+----2----+----3----+----4
130  #endif  #endif
131    
132  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
133          recipMaxSlope = 0. _d 0
134          IF ( GM_Visbeck_maxSlope.GT.0. _d 0 ) THEN
135            recipMaxSlope = 1. _d 0 / GM_Visbeck_maxSlope
136          ENDIF
137        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
138         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
139          VisbeckK(i,j,bi,bj) = 0. _d 0          VisbeckK(i,j,bi,bj) = 0. _d 0
# Line 255  C-- 1rst loop on k : compute Tensor Coef Line 260  C-- 1rst loop on k : compute Tensor Coef
260  # endif  # endif
261          ENDDO          ENDDO
262         ENDDO         ENDDO
263  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
264    
265         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
266          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
# Line 266  C      Gradient of Sigma at rVel points Line 271  C      Gradient of Sigma at rVel points
271           dSigmaDy(i,j)=op25*( sigmaY(i,j+1,k-1)+sigmaY(i,j,k-1)           dSigmaDy(i,j)=op25*( sigmaY(i,j+1,k-1)+sigmaY(i,j,k-1)
272       &                       +sigmaY(i,j+1, k )+sigmaY(i,j, k )       &                       +sigmaY(i,j+1, k )+sigmaY(i,j, k )
273       &                      )*maskC(i,j,k,bi,bj)       &                      )*maskC(i,j,k,bi,bj)
274           dSigmaDr(i,j)=sigmaR(i,j,k)  c        dSigmaDr(i,j)=sigmaR(i,j,k)
275          ENDDO          ENDDO
276         ENDDO         ENDDO
277    
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
278  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
279  #ifndef OLD_VISBECK_CALC  #ifndef OLD_VISBECK_CALC
280         IF ( GM_Visbeck_alpha.GT.0. .AND.         IF ( GM_Visbeck_alpha.GT.0. .AND.
281       &      -rC(k-1).LT.GM_Visbeck_depth ) THEN       &      -rC(k-1).LT.GM_Visbeck_depth ) THEN
282    
283            DO j=1-Oly,sNy+Oly
284             DO i=1-Olx,sNx+Olx
285              dSigmaDr(i,j) = MIN( sigmaR(i,j,k), 0. _d 0 )
286             ENDDO
287            ENDDO
288    
289  C--     Depth average of f/sqrt(Ri) = M^2/N^2 * N  C--     Depth average of f/sqrt(Ri) = M^2/N^2 * N
290  C       M^2 and N^2 are horizontal & vertical gradient of buoyancy.  C       M^2 and N^2 are horizontal & vertical gradient of buoyancy.
291    
# Line 298  C       GM_Visbeck_depth, whatever is th Line 300  C       GM_Visbeck_depth, whatever is th
300             integrDepth = -rC( kLowC(i,j,bi,bj) )             integrDepth = -rC( kLowC(i,j,bi,bj) )
301  C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments  C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
302             integrDepth = MIN( integrDepth, GM_Visbeck_depth )             integrDepth = MIN( integrDepth, GM_Visbeck_depth )
303    C-      to recover "old-visbeck" form with Visbeck_minDepth = Visbeck_depth
304               integrDepth = MAX( integrDepth, GM_Visbeck_minDepth )
305  C       Distance between level center above and the integration depth  C       Distance between level center above and the integration depth
306             deltaH = integrDepth + rC(k-1)             deltaH = integrDepth + rC(k-1)
307  C       If negative then we are below the integration level  C       If negative then we are below the integration level
# Line 308  C       Now we convert deltaH to a non-d Line 312  C       Now we convert deltaH to a non-d
312             deltaH = deltaH/( integrDepth+rC(1) )             deltaH = deltaH/( integrDepth+rC(1) )
313    
314  C--      compute: ( M^2 * S )^1/2   (= S*N since S=M^2/N^2 )  C--      compute: ( M^2 * S )^1/2   (= S*N since S=M^2/N^2 )
315    C        a 5 points average gives a more "homogeneous" formulation
316    C        (same stencil and same weights as for dSigmaH calculation)
317               dSigmaR = ( dSigmaDr(i,j)*4. _d 0
318         &               + dSigmaDr(i-1,j)
319         &               + dSigmaDr(i+1,j)
320         &               + dSigmaDr(i,j-1)
321         &               + dSigmaDr(i,j+1)
322         &               )/( 4. _d 0
323         &                 + maskC(i-1,j,k,bi,bj)
324         &                 + maskC(i+1,j,k,bi,bj)
325         &                 + maskC(i,j-1,k,bi,bj)
326         &                 + maskC(i,j+1,k,bi,bj)
327         &                 )
328             dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)             dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)
329       &             + dSigmaDy(i,j)*dSigmaDy(i,j)       &             + dSigmaDy(i,j)*dSigmaDy(i,j)
330             IF ( dSigmaH .GT. 0. _d 0 ) THEN             IF ( dSigmaH .GT. 0. _d 0 ) THEN
331               dSigmaH = SQRT( dSigmaH )               dSigmaH = SQRT( dSigmaH )
332  C-       compute slope, limited by GM_maxSlope:  C-       compute slope, limited by GM_Visbeck_maxSlope:
333               IF ( -dSigmaDr(i,j).GT.dSigmaH*GM_rMaxSlope ) THEN               IF ( -dSigmaR.GT.dSigmaH*recipMaxSlope ) THEN
334                Sloc = dSigmaH / ( -dSigmaDr(i,j) )                Sloc = dSigmaH / ( -dSigmaR )
335               ELSE               ELSE
336                Sloc = GM_maxSlope                Sloc = GM_Visbeck_maxSlope
337               ENDIF               ENDIF
338               M2loc = gravity*recip_rhoConst*dSigmaH               M2loc = gravity*recip_rhoConst*dSigmaH
339  c            SNloc = SQRT( Sloc*M2loc )  c            SNloc = SQRT( Sloc*M2loc )
340               N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)               N2loc = -gravity*recip_rhoConst*dSigmaR
341    c            N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)
342               IF ( N2loc.GT.0. _d 0 ) THEN               IF ( N2loc.GT.0. _d 0 ) THEN
343                 SNloc = Sloc*SQRT(N2loc)                 SNloc = Sloc*SQRT(N2loc)
344               ELSE               ELSE
# Line 338  c            SNloc = SQRT( Sloc*M2loc ) Line 356  c            SNloc = SQRT( Sloc*M2loc )
356         ENDIF         ENDIF
357  #endif /* ndef OLD_VISBECK_CALC */  #endif /* ndef OLD_VISBECK_CALC */
358  #endif /* GM_VISBECK_VARIABLE_K */  #endif /* GM_VISBECK_VARIABLE_K */
359           DO j=1-Oly,sNy+Oly
360            DO i=1-Olx,sNx+Olx
361             dSigmaDr(i,j)=sigmaR(i,j,k)
362            ENDDO
363           ENDDO
364    
365    #ifdef ALLOW_AUTODIFF_TAMC
366    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
367    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
368    CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
369    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
370    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
371    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
372    #endif /* ALLOW_AUTODIFF_TAMC */
373    
374  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
375         CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
# Line 425  C-     Limit range that KapGM can take Line 457  C-     Limit range that KapGM can take
457         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
458          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
459           VisbeckK(i,j,bi,bj)=           VisbeckK(i,j,bi,bj)=
460       &       MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)       &       MIN( MAX( VisbeckK(i,j,bi,bj), GM_Visbeck_minVal_K ),
461         &            GM_Visbeck_maxVal_K )
462          ENDDO          ENDDO
463         ENDDO         ENDDO
464        ENDIF        ENDIF
Legend:
Removed from v.1.33  
changed lines
  Added in v.1.34
  ViewVC Help
Powered by ViewVC 1.1.22