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

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

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

revision 1.19 by jmc, Fri Dec 5 17:39:39 2014 UTC revision 1.20 by jmc, Tue Jan 20 20:52:12 2015 UTC
# Line 25  C     == Global variables == Line 25  C     == Global variables ==
25  #include "PARAMS.h"  #include "PARAMS.h"
26  #include "GRID.h"  #include "GRID.h"
27  #include "DYNVARS.h"  #include "DYNVARS.h"
28    #include "FFIELDS.h"
29  #include "GMREDI.h"  #include "GMREDI.h"
30    
31  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
# Line 208  C     psistar :: eddy induced streamfunc Line 209  C     psistar :: eddy induced streamfunc
209        _RL psistar(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)        _RL psistar(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)
210  #endif  #endif
211    
   
212  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
213    
214  C     ======================================  C     ======================================
# Line 257  C     Gradient of Coriolis Line 257  C     Gradient of Coriolis
257         ENDDO         ENDDO
258        ENDDO        ENDDO
259    
260  C     Coriolis at C points enforcing a minimum value so  C     Coriolis at C points enforcing a minimum value so
261  C     that it is defined at the equator  C     that it is defined at the equator
262        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
263         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
# Line 270  C     Coriolis at U and V points Line 270  C     Coriolis at U and V points
270         DO i=1-Olx+1,sNx+Olx         DO i=1-Olx+1,sNx+Olx
271  C       Limited so that the inverse is defined at the equator  C       Limited so that the inverse is defined at the equator
272          coriU(i,j) = op5*( cori(i,j)+cori(i-1,j) )          coriU(i,j) = op5*( cori(i,j)+cori(i-1,j) )
273          coriU(i,j) = SIGN( MAX( ABS(coriU(i,j)),GM_K3D_minCori ),          coriU(i,j) = SIGN( MAX( ABS(coriU(i,j)),GM_K3D_minCori ),
274       &       coriU(i,j) )       &       coriU(i,j) )
275    
276  C       Not limited  C       Not limited
# Line 281  C       Not limited Line 281  C       Not limited
281         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
282  C       Limited so that the inverse is defined at the equator  C       Limited so that the inverse is defined at the equator
283          coriV(i,j) = op5*( cori(i,j)+cori(i,j-1) )          coriV(i,j) = op5*( cori(i,j)+cori(i,j-1) )
284          coriV(i,j) = SIGN( MAX( ABS(coriV(i,j)),GM_K3D_minCori ),          coriV(i,j) = SIGN( MAX( ABS(coriV(i,j)),GM_K3D_minCori ),
285       &       coriV(i,j) )       &       coriV(i,j) )
286    
287  C       Not limited  C       Not limited
# Line 355  C     initialise remaining 3d variables Line 355  C     initialise remaining 3d variables
355  C     Find the zonal velocity at the cell centre  C     Find the zonal velocity at the cell centre
356  C     The logicals here are, in order: 1/ go from grid to north/east directions  C     The logicals here are, in order: 1/ go from grid to north/east directions
357  C     2/ go from C to A grid and 3/ apply the mask  C     2/ go from C to A grid and 3/ apply the mask
358  #ifdef ALLOW_EDDYPSI        #ifdef ALLOW_EDDYPSI
359        IF (GM_InMomAsStress) THEN        IF (GM_InMomAsStress) THEN
360          CALL rotate_uv2en_rl(uMean, vMean, ubar, vbar, .TRUE., .TRUE.,          CALL ROTATE_UV2EN_RL( uEulerMean, vEulerMean, ubar, vbar,
361       &       .TRUE.,Nr,mythid)       &                        .TRUE., .TRUE., .TRUE., Nr, mythid )
362        ELSE        ELSE
363  #endif  #endif
364          CALL rotate_uv2en_rl(uVel, vVel, ubar, vbar, .TRUE., .TRUE.,          CALL ROTATE_UV2EN_RL( uVel, vVel, ubar, vbar,
365       &       .TRUE.,Nr,mythid)       &                        .TRUE., .TRUE., .TRUE., Nr, mythid )
366  #ifdef ALLOW_EDDYPSI        #ifdef ALLOW_EDDYPSI
367        ENDIF        ENDIF
368  #endif  #endif
369    
# Line 383  C     N2(k=1) is always zero Line 383  C     N2(k=1) is always zero
383          N(i,j,k)  = zeroRL          N(i,j,k)  = zeroRL
384         ENDDO         ENDDO
385        ENDDO        ENDDO
386  C     Enforce a minimum N2  C     Enforce a minimum N2
387        DO k=2,Nr        DO k=2,Nr
388         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
389          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
# Line 395  C     Enforce a minimum N2 Line 395  C     Enforce a minimum N2
395  C     Calculate the minimum drho/dz  C     Calculate the minimum drho/dz
396        maxDRhoDz = -rhoConst*GM_K3D_minN2/gravity        maxDRhoDz = -rhoConst*GM_K3D_minN2/gravity
397    
398  C     Calculate the barotropic velocity by vertically integrating  C     Calculate the barotropic velocity by vertically integrating
399  C     and the dividing by the depth of the water column  C     and the dividing by the depth of the water column
400  C     Note that Ubaro is on the U grid.  C     Note that Ubaro is on the U grid.
401        DO k=1,Nr        DO k=1,Nr
402         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
403          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
404           Ubaro(i,j) = Ubaro(i,j) +           Ubaro(i,j) = Ubaro(i,j) +
405       &        maskW(i,j,k,bi,bj)*drF(k)*hfacC(i,j,k,bi,bj)       &        maskW(i,j,k,bi,bj)*drF(k)*hfacC(i,j,k,bi,bj)
406       &                *ubar(i,j,k,bi,bj)       &                *ubar(i,j,k,bi,bj)
407          ENDDO          ENDDO
# Line 455  C       Calculate the Rossby Radius Line 455  C       Calculate the Rossby Radius
455        ENDIF        ENDIF
456    
457  C     Average dsigma/dx and dsigma/dy onto the centre points  C     Average dsigma/dx and dsigma/dy onto the centre points
458          
459        DO k=1,Nr        DO k=1,Nr
460         DO j=1-Oly,sNy+Oly-1         DO j=1-Oly,sNy+Oly-1
461          DO i=1-Olx,sNx+Olx-1          DO i=1-Olx,sNx+Olx-1
# Line 472  C     =============================== Line 472  C     ===============================
472    
473         DO j=1-Oly,sNy+Oly-1         DO j=1-Oly,sNy+Oly-1
474          DO i=1-Olx,sNx+Olx-1          DO i=1-Olx,sNx+Olx-1
475           M4loc(i,j,k) = g_reciprho_sq*( dSigmaDx(i,j,k)**2           M4loc(i,j,k) = g_reciprho_sq*( dSigmaDx(i,j,k)**2
476       &                                 +dSigmaDy(i,j,k)**2 )       &                                 +dSigmaDy(i,j,k)**2 )
477           IF (k.NE.kLow_C(i,j)) THEN           IF (k.NE.kLow_C(i,j)) THEN
478             N2loc(i,j,k) = op5*(N2(i,j,k)+N2(i,j,k+1))             N2loc(i,j,k) = op5*(N2(i,j,k)+N2(i,j,k+1))
# Line 492  C      cell is deeper than the bottom in Line 492  C      cell is deeper than the bottom in
492  C      We are in the integration depth range  C      We are in the integration depth range
493         DO j=1-Oly,sNy+Oly-1         DO j=1-Oly,sNy+Oly-1
494          DO i=1-Olx,sNx+Olx-1          DO i=1-Olx,sNx+Olx-1
495           IF ( (kLow_C(i,j).GE.k) .AND.           IF ( (kLow_C(i,j).GE.k) .AND.
496       &        (-hMixLayer(i,j,bi,bj).LE.-rC(k)) ) THEN       &        (-hMixLayer(i,j,bi,bj).LE.-rC(k)) ) THEN
497    
498             slopeC(i,j,k) = SQRT(M4loc(i,j,k))/N2loc(i,j,k)             slopeC(i,j,k) = SQRT(M4loc(i,j,k))/N2loc(i,j,k)
# Line 503  C          Limit the slope.  Note, this Line 503  C          Limit the slope.  Note, this
503               slopeC(i,j,k) = GM_maxslope               slopeC(i,j,k) = GM_maxslope
504               M4onN2(i,j,k) = SQRT(M4loc(i,j,k))*GM_maxslope               M4onN2(i,j,k) = SQRT(M4loc(i,j,k))*GM_maxslope
505             ENDIF             ENDIF
506             eady(i,j)   = eady(i,j)               eady(i,j)   = eady(i,j)
507       &                   + hfacC(i,j,k,bi,bj)*drF(k)*M4onN2(i,j,k)       &                   + hfacC(i,j,k,bi,bj)*drF(k)*M4onN2(i,j,k)
508             deltaH(i,j) = deltaH(i,j) + drF(k)             deltaH(i,j) = deltaH(i,j) + drF(k)
509           ENDIF           ENDIF
# Line 514  C          Limit the slope.  Note, this Line 514  C          Limit the slope.  Note, this
514        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
515         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
516  C       If the minimum depth for the integration is deeper than the ocean  C       If the minimum depth for the integration is deeper than the ocean
517  C       bottom OR the mixed layer is deeper than the maximum depth of  C       bottom OR the mixed layer is deeper than the maximum depth of
518  C       integration, we set the Eady growth rate to something small  C       integration, we set the Eady growth rate to something small
519  C       to avoid floating point exceptions.  C       to avoid floating point exceptions.
520  C       Later, these areas will be given a small diffusivity.  C       Later, these areas will be given a small diffusivity.
# Line 525  C       Otherwise, divide over the integ Line 525  C       Otherwise, divide over the integ
525  C       to actually find the Eady growth rate.  C       to actually find the Eady growth rate.
526          ELSE          ELSE
527            eady(i,j) = SQRT(eady(i,j)/deltaH(i,j))            eady(i,j) = SQRT(eady(i,j)/deltaH(i,j))
528              
529          ENDIF          ENDIF
530    
531         ENDDO         ENDDO
# Line 600  C     ================================== Line 600  C     ==================================
600  C     Find the PV gradient  C     Find the PV gradient
601  C     ======================================  C     ======================================
602  C     Calculate the surface layer thickness.  C     Calculate the surface layer thickness.
603  C     Use hMixLayer (calculated in model/src/calc_oce_mxlayer)  C     Use hMixLayer (calculated in model/src/calc_oce_mxlayer)
604  C     for the mixed layer depth.  C     for the mixed layer depth.
605    
606  C     Enforce a minimum surface layer depth  C     Enforce a minimum surface layer depth
# Line 615  C     Enforce a minimum surface layer de Line 615  C     Enforce a minimum surface layer de
615        DO k=1,Nr        DO k=1,Nr
616         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
617          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
618           IF (rF(k).GT.surfkz(i,j) .AND. surfkz(i,j).GE.rF(k+1))           IF (rF(k).GT.surfkz(i,j) .AND. surfkz(i,j).GE.rF(k+1))
619       &        surfk(i,j) = k       &        surfk(i,j) = k
620          ENDDO          ENDDO
621         ENDDO         ENDDO
# Line 645  C          Calculate the zonal slope at Line 645  C          Calculate the zonal slope at
645             IF(ABS(slope).GT.GM_maxSlope)             IF(ABS(slope).GT.GM_maxSlope)
646       &          slope = SIGN(GM_maxSlope,slope)       &          slope = SIGN(GM_maxSlope,slope)
647             SlopeX(i,j,k)=-maskW(i,j,k-1,bi,bj)*maskW(i,j,k,bi,bj)*slope             SlopeX(i,j,k)=-maskW(i,j,k-1,bi,bj)*maskW(i,j,k,bi,bj)*slope
648              
649  C          Calculate the meridional slope at the southern cell face (V grid)  C          Calculate the meridional slope at the southern cell face (V grid)
650             sigz  = MIN( op5*(sigmaR(i,j,k)+sigmaR(i,j-1,k)), maxDRhoDz )             sigz  = MIN( op5*(sigmaR(i,j,k)+sigmaR(i,j-1,k)), maxDRhoDz )
651             sigy  = op5*( sigmaY(i,j,k) + sigmaY(i,j,k-1) )             sigy  = op5*( sigmaY(i,j,k) + sigmaY(i,j,k-1) )
# Line 658  C          Calculate the meridional slop Line 658  C          Calculate the meridional slop
658         ENDDO         ENDDO
659        ENDDO        ENDDO
660    
661  C     Calculate the thickness flux and diffusivity.  These may be altered later  C     Calculate the thickness flux and diffusivity.  These may be altered later
662  C     depending on namelist options.  C     depending on namelist options.
663  C     Enforce a zero slope bottom boundary condition for the bottom most cells (k=Nr)  C     Enforce a zero slope bottom boundary condition for the bottom most cells (k=Nr)
664        k=Nr        k=Nr
# Line 671  C          Meridional thickness flux at Line 671  C          Meridional thickness flux at
671             tfluxY(i,j,k) = -fCoriV(i,j)*SlopeY(i,j,k)             tfluxY(i,j,k) = -fCoriV(i,j)*SlopeY(i,j,k)
672       &                     *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)       &                     *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
673    
674  C          Use the interior diffusivity. Note: if we are using a  C          Use the interior diffusivity. Note: if we are using a
675  C          constant diffusivity KPV is overwritten later  C          constant diffusivity KPV is overwritten later
676             KPV(i,j,k) = K3D(i,j,k,bi,bj)             KPV(i,j,k) = K3D(i,j,k,bi,bj)
677    
# Line 686  C        Zonal thickness flux at the wes Line 686  C        Zonal thickness flux at the wes
686           tfluxX(i,j,k)=-fCoriU(i,j)*(SlopeX(i,j,k)-SlopeX(i,j,k+1))           tfluxX(i,j,k)=-fCoriU(i,j)*(SlopeX(i,j,k)-SlopeX(i,j,k+1))
687       &        *recip_drF(k)*recip_hFacW(i,j,k,bi,bj)       &        *recip_drF(k)*recip_hFacW(i,j,k,bi,bj)
688       &        *maskW(i,j,k,bi,bj)       &        *maskW(i,j,k,bi,bj)
689            
690  C        Meridional thickness flux at the southern cell face  C        Meridional thickness flux at the southern cell face
691           tfluxY(i,j,k)=-fCoriV(i,j)*(SlopeY(i,j,k)-SlopeY(i,j,k+1))           tfluxY(i,j,k)=-fCoriV(i,j)*(SlopeY(i,j,k)-SlopeY(i,j,k+1))
692       &        *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)       &        *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
693       &        *maskS(i,j,k,bi,bj)       &        *maskS(i,j,k,bi,bj)
694            
695  C        Use the interior diffusivity. Note: if we are using a  C        Use the interior diffusivity. Note: if we are using a
696  C        constant diffusivity KPV is overwritten later  C        constant diffusivity KPV is overwritten later
697           KPV(i,j,k) = K3D(i,j,k,bi,bj)           KPV(i,j,k) = K3D(i,j,k,bi,bj)
698          ENDDO          ENDDO
699         ENDDO         ENDDO
700        ENDDO        ENDDO
701    
702    C     Special treatment for the surface layer (if necessary), which overwrites
 C     Special treatment for the surface layer (if necessary), which overwrites  
703  C     values in the previous loops.  C     values in the previous loops.
704        IF (GM_K3D_ThickSheet .OR. GM_K3D_surfK) THEN        IF (GM_K3D_ThickSheet .OR. GM_K3D_surfK) THEN
705          DO k=Nr-1,1,-1          DO k=Nr-1,1,-1
706           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
707            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
708             IF(k.LE.surfk(i,j)) THEN             IF(k.LE.surfk(i,j)) THEN
709  C            We are in the surface layer.  Change the thickness flux  C            We are in the surface layer.  Change the thickness flux
710  C            and diffusivity as necessary.  C            and diffusivity as necessary.
711    
712               IF (GM_K3D_ThickSheet) THEN               IF (GM_K3D_ThickSheet) THEN
713  C              We are in the surface layer, so set the thickness flux  C              We are in the surface layer, so set the thickness flux
714  C              based on the average slope over the surface layer  C              based on the average slope over the surface layer
715  C              If we are on the edge of a "cliff" the surface layer at the  C              If we are on the edge of a "cliff" the surface layer at the
716  C              centre of the grid point could be deeper than the U or V point.    C              centre of the grid point could be deeper than the U or V point.
717  C              So, we ensure that we always take a sensible slope.  C              So, we ensure that we always take a sensible slope.
718                 IF(kLow_U(i,j).LT.surfk(i,j)) THEN                 IF(kLow_U(i,j).LT.surfk(i,j)) THEN
719                   kk=kLow_U(i,j)                   kk=kLow_U(i,j)
# Line 734  C              So, we ensure that we alw Line 733  C              So, we ensure that we alw
733                 ELSE                 ELSE
734                   tfluxX(i,j,k) = zeroRL                   tfluxX(i,j,k) = zeroRL
735                 ENDIF                 ENDIF
736                  
737                 IF(kLow_V(i,j).LT.surfk(i,j)) THEN                 IF(kLow_V(i,j).LT.surfk(i,j)) THEN
738                   kk=kLow_V(i,j)                   kk=kLow_V(i,j)
739                   hsurf = -rLowS(i,j,bi,bj)                   hsurf = -rLowS(i,j,bi,bj)
# Line 776  C       Ignore beta in the calculation o Line 775  C       Ignore beta in the calculation o
775            ENDDO            ENDDO
776           ENDDO           ENDDO
777          ENDDO          ENDDO
778            
779        ELSE        ELSE
780  C       Do not ignore beta  C       Do not ignore beta
781          DO k=1,Nr          DO k=1,Nr
# Line 821  C     If GM_K3D_constRedi=.TRUE. K3D wil Line 820  C     If GM_K3D_constRedi=.TRUE. K3D wil
820    
821        IF (.NOT. GM_K3D_smooth) THEN        IF (.NOT. GM_K3D_smooth) THEN
822  C       Do not expand K grad(q) => no smoothing  C       Do not expand K grad(q) => no smoothing
823  C       May only be done with a constant K, otherwise the  C       May only be done with a constant K, otherwise the
824  C       integral constraint is violated.  C       integral constraint is violated.
825          DO k=1,Nr          DO k=1,Nr
826           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
# Line 847  C       Calculate the eigenvectors at th Line 846  C       Calculate the eigenvectors at th
846       I         rLowW(:,:,bi,bj),GM_K3D_NModes,.FALSE.,       I         rLowW(:,:,bi,bj),GM_K3D_NModes,.FALSE.,
847       O         dummy,modesW(:,:,:,:,bi,bj))       O         dummy,modesW(:,:,:,:,bi,bj))
848          ENDIF          ENDIF
849            
850  C       Calculate Xi_m at the west face of a cell  C       Calculate Xi_m at the west face of a cell
851          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
852           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
# Line 861  C       Calculate Xi_m at the west face Line 860  C       Calculate Xi_m at the west face
860            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
861             DO m=1,GM_K3D_NModes             DO m=1,GM_K3D_NModes
862              Kdqdx(i,j,k) = KPV(i,j,k)*gradqx(i,j,k)              Kdqdx(i,j,k) = KPV(i,j,k)*gradqx(i,j,k)
863              XimX(m,i,j) = XimX(m,i,j)              XimX(m,i,j) = XimX(m,i,j)
864       &           - maskW(i,j,k,bi,bj)*drF(k)*hfacW(i,j,k,bi,bj)       &           - maskW(i,j,k,bi,bj)*drF(k)*hfacW(i,j,k,bi,bj)
865       &           *Kdqdx(i,j,k)*modesW(m,i,j,k,bi,bj)       &           *Kdqdx(i,j,k)*modesW(m,i,j,k,bi,bj)
866             ENDDO             ENDDO
867            ENDDO            ENDDO
868           ENDDO           ENDDO
869          ENDDO          ENDDO
870              
871  C     Calculate Xi in the X direction at the west face  C     Calculate Xi in the X direction at the west face
872          DO k=1,Nr          DO k=1,Nr
873           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
# Line 881  C     Calculate Xi in the X direction at Line 880  C     Calculate Xi in the X direction at
880           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
881            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
882             DO m=1,GM_K3D_NModes             DO m=1,GM_K3D_NModes
883              Xix(i,j,k) = Xix(i,j,k)              Xix(i,j,k) = Xix(i,j,k)
884       &           + maskW(i,j,k,bi,bj)*XimX(m,i,j)*modesW(m,i,j,k,bi,bj)       &           + maskW(i,j,k,bi,bj)*XimX(m,i,j)*modesW(m,i,j,k,bi,bj)
885             ENDDO             ENDDO
886            ENDDO            ENDDO
# Line 919  C     Calculate the eigenvectors at the Line 918  C     Calculate the eigenvectors at the
918            ENDDO            ENDDO
919           ENDDO           ENDDO
920          ENDDO          ENDDO
921            
922  C     Calculate Xi for Y direction at the south face  C     Calculate Xi for Y direction at the south face
923          DO k=1,Nr          DO k=1,Nr
924           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
# Line 932  C     Calculate Xi for Y direction at th Line 931  C     Calculate Xi for Y direction at th
931           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
932            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
933             DO m=1,GM_K3D_NModes             DO m=1,GM_K3D_NModes
934              Xiy(i,j,k) = Xiy(i,j,k)              Xiy(i,j,k) = Xiy(i,j,k)
935       &           + maskS(i,j,k,bi,bj)*XimY(m,i,j)*modesS(m,i,j,k,bi,bj)       &           + maskS(i,j,k,bi,bj)*XimY(m,i,j)*modesS(m,i,j,k,bi,bj)
936             ENDDO             ENDDO
937            ENDDO            ENDDO
# Line 942  C     Calculate Xi for Y direction at th Line 941  C     Calculate Xi for Y direction at th
941  C     ENDIF (.NOT. GM_K3D_smooth)  C     ENDIF (.NOT. GM_K3D_smooth)
942        ENDIF        ENDIF
943    
   
944  C     Calculate the renormalisation factor  C     Calculate the renormalisation factor
945        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
946         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
# Line 965  C     Calculate the renormalisation fact Line 963  C     Calculate the renormalisation fact
963           centreX = op5*(uVel(i,j,k,bi,bj)+uVel(i+1,j,k,bi,bj))           centreX = op5*(uVel(i,j,k,bi,bj)+uVel(i+1,j,k,bi,bj))
964           centreY = op5*(Kdqdy(i,j,k)     +Kdqdy(i,j+1,k)     )           centreY = op5*(Kdqdy(i,j,k)     +Kdqdy(i,j+1,k)     )
965  C        For the numerator  C        For the numerator
966           uInt(i,j) = uInt(i,j)           uInt(i,j) = uInt(i,j)
967       &        + centreX*hfacC(i,j,k,bi,bj)*drF(k)       &        + centreX*hfacC(i,j,k,bi,bj)*drF(k)
968           KdqdyInt(i,j) = KdqdyInt(i,j)           KdqdyInt(i,j) = KdqdyInt(i,j)
969       &        + centreY*hfacC(i,j,k,bi,bj)*drF(k)       &        + centreY*hfacC(i,j,k,bi,bj)*drF(k)
# Line 997  C        For the denominator Line 995  C        For the denominator
995        DO j=1-Oly,sNy+Oly-1        DO j=1-Oly,sNy+Oly-1
996         DO i=1-Olx,sNx+Olx-1         DO i=1-Olx,sNx+Olx-1
997          IF (kLowC(i,j,bi,bj).GT.0) THEN          IF (kLowC(i,j,bi,bj).GT.0) THEN
998            numerator =            numerator =
999       &         (uKdqdyInt(i,j)-uInt(i,j)*KdqdyInt(i,j)/R_low(i,j,bi,bj))       &         (uKdqdyInt(i,j)-uInt(i,j)*KdqdyInt(i,j)/R_low(i,j,bi,bj))
1000       &        -(vKdqdxInt(i,j)-vInt(i,j)*KdqdxInt(i,j)/R_low(i,j,bi,bj))       &        -(vKdqdxInt(i,j)-vInt(i,j)*KdqdxInt(i,j)/R_low(i,j,bi,bj))
1001            denominator = uXiyInt(i,j) - vXixInt(i,j)            denominator = uXiyInt(i,j) - vXixInt(i,j)
1002  C         We can have troubles with floating point exceptions if the denominator  C         We can have troubles with floating point exceptions if the denominator
1003  C         of the renormalisation if the ocean is resting (e.g. intial conditions).  C         of the renormalisation if the ocean is resting (e.g. intial conditions).
1004  C         So we make the renormalisation factor one if the denominator is very small  C         So we make the renormalisation factor one if the denominator is very small
1005  C         The renormalisation factor is supposed to correct the error in the extraction of  C         The renormalisation factor is supposed to correct the error in the extraction of
1006  C         potential energy associated with the truncation of the expansion. Thus, we  C         potential energy associated with the truncation of the expansion. Thus, we
# Line 1031  C     Calculate the eddy induced velocit Line 1029  C     Calculate the eddy induced velocit
1029           ustar(i,j,k) = -RenormU(i,j)*Xix(i,j,k)/coriU(i,j)           ustar(i,j,k) = -RenormU(i,j)*Xix(i,j,k)/coriU(i,j)
1030          ENDDO          ENDDO
1031         ENDDO         ENDDO
1032        ENDDO              ENDDO
1033    
1034  C     Calculate the eddy induced velocity in the Y direction at the south face  C     Calculate the eddy induced velocity in the Y direction at the south face
1035        DO k=1,Nr        DO k=1,Nr
# Line 1040  C     Calculate the eddy induced velocit Line 1038  C     Calculate the eddy induced velocit
1038           vstar(i,j,k) = -RenormV(i,j)*Xiy(i,j,k)/coriV(i,j)           vstar(i,j,k) = -RenormV(i,j)*Xiy(i,j,k)/coriV(i,j)
1039          ENDDO          ENDDO
1040         ENDDO         ENDDO
1041        ENDDO              ENDDO
1042    
1043  C     ======================================  C     ======================================
1044  C     Calculate the eddy induced overturning streamfunction  C     Calculate the eddy induced overturning streamfunction
# Line 1060  C     ================================== Line 1058  C     ==================================
1058          ENDDO          ENDDO
1059         ENDDO         ENDDO
1060        ENDDO        ENDDO
1061        
1062  #else  #else
1063    
1064        IF (GM_AdvForm) THEN        IF (GM_AdvForm) THEN
# Line 1088  C     ================================== Line 1086  C     ==================================
1086  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
1087  C     Diagnostics  C     Diagnostics
1088        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
1089          CALL DIAGNOSTICS_FILL(K3D,    'GM_K3D  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(K3D,    'GM_K3D  ',0,Nr,1,bi,bj,myThid)
1090          CALL DIAGNOSTICS_FILL(KPV,    'GM_KPV  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(KPV,    'GM_KPV  ',0,Nr,2,bi,bj,myThid)
1091          CALL DIAGNOSTICS_FILL(urms,   'GM_URMS ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(urms,   'GM_URMS ',0,Nr,2,bi,bj,myThid)
1092          CALL DIAGNOSTICS_FILL(Rdef,   'GM_RDEF ',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Rdef,   'GM_RDEF ',0, 1,1,bi,bj,myThid)
1093          CALL DIAGNOSTICS_FILL(Rurms,  'GM_RURMS',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Rurms,  'GM_RURMS',0, 1,2,bi,bj,myThid)
1094          CALL DIAGNOSTICS_FILL(RRhines,'GM_RRHNS',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(RRhines,'GM_RRHNS',0, 1,2,bi,bj,myThid)
1095          CALL DIAGNOSTICS_FILL(Rmix,   'GM_RMIX ',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Rmix,   'GM_RMIX ',0, 1,2,bi,bj,myThid)
1096          CALL DIAGNOSTICS_FILL(supp,   'GM_SUPP ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(supp,   'GM_SUPP ',0,Nr,2,bi,bj,myThid)
1097          CALL DIAGNOSTICS_FILL(Xix,    'GM_Xix  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Xix,    'GM_Xix  ',0,Nr,2,bi,bj,myThid)
1098          CALL DIAGNOSTICS_FILL(Xiy,    'GM_Xiy  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Xiy,    'GM_Xiy  ',0,Nr,2,bi,bj,myThid)
1099          CALL DIAGNOSTICS_FILL(cDopp,  'GM_C    ',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(cDopp,  'GM_C    ',0, 1,2,bi,bj,myThid)
1100          CALL DIAGNOSTICS_FILL(Ubaro,  'GM_UBARO',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Ubaro,  'GM_UBARO',0, 1,2,bi,bj,myThid)
1101          CALL DIAGNOSTICS_FILL(eady,   'GM_EADY ',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(eady,   'GM_EADY ',0, 1,2,bi,bj,myThid)
1102          CALL DIAGNOSTICS_FILL(SlopeX, 'GM_Sx   ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(SlopeX, 'GM_Sx   ',0,Nr,2,bi,bj,myThid)
1103          CALL DIAGNOSTICS_FILL(SlopeY, 'GM_Sy   ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(SlopeY, 'GM_Sy   ',0,Nr,2,bi,bj,myThid)
1104          CALL DIAGNOSTICS_FILL(tfluxX, 'GM_TFLXX',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(tfluxX, 'GM_TFLXX',0,Nr,2,bi,bj,myThid)
1105          CALL DIAGNOSTICS_FILL(tfluxY, 'GM_TFLXY',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(tfluxY, 'GM_TFLXY',0,Nr,2,bi,bj,myThid)
1106          CALL DIAGNOSTICS_FILL(gradqx, 'GM_dqdx ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(gradqx, 'GM_dqdx ',0,Nr,2,bi,bj,myThid)
1107          CALL DIAGNOSTICS_FILL(gradqy, 'GM_dqdy ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(gradqy, 'GM_dqdy ',0,Nr,2,bi,bj,myThid)
1108          CALL DIAGNOSTICS_FILL(Kdqdy,  'GM_Kdqdy',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Kdqdy,  'GM_Kdqdy',0,Nr,2,bi,bj,myThid)
1109          CALL DIAGNOSTICS_FILL(Kdqdx,  'GM_Kdqdx',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Kdqdx,  'GM_Kdqdx',0,Nr,2,bi,bj,myThid)
1110          CALL DIAGNOSTICS_FILL(surfkz, 'GM_SFLYR',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(surfkz, 'GM_SFLYR',0, 1,2,bi,bj,myThid)
1111          CALL DIAGNOSTICS_FILL(ustar,  'GM_USTAR',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(ustar,  'GM_USTAR',0,Nr,2,bi,bj,myThid)
1112          CALL DIAGNOSTICS_FILL(vstar,  'GM_VSTAR',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(vstar,  'GM_VSTAR',0,Nr,2,bi,bj,myThid)
1113          CALL DIAGNOSTICS_FILL(umc,    'GM_UMC  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(umc,    'GM_UMC  ',0,Nr,2,bi,bj,myThid)
1114          CALL DIAGNOSTICS_FILL(ubar,   'GM_UBAR ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(ubar,   'GM_UBAR ',0,Nr,2,bi,bj,myThid)
1115          CALL DIAGNOSTICS_FILL(modesC(1,:,:,:,bi,bj),          CALL DIAGNOSTICS_FILL(modesC, 'GM_MODEC',0,Nr,1,bi,bj,myThid)
1116       &                                'GM_MODEC',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(M4loc,  'GM_M4   ',0,Nr,2,bi,bj,myThid)
1117          CALL DIAGNOSTICS_FILL(M4loc,  'GM_M4   ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(N2loc,  'GM_N2   ',0,Nr,2,bi,bj,myThid)
1118          CALL DIAGNOSTICS_FILL(N2loc,  'GM_N2   ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(M4onN2, 'GM_M4_N2',0,Nr,2,bi,bj,myThid)
1119          CALL DIAGNOSTICS_FILL(M4onN2, 'GM_M4_N2',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(slopeC, 'GM_SLOPE',0,Nr,2,bi,bj,myThid)
1120          CALL DIAGNOSTICS_FILL(slopeC, 'GM_SLOPE',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Renorm, 'GM_RENRM',0, 1,2,bi,bj,myThid)
         CALL DIAGNOSTICS_FILL(Renorm, 'GM_RENRM',0, 1,0,1,1,myThid)  
   
1121        ENDIF        ENDIF
1122  #endif  #endif
1123    
1124  C     For the Redi diffusivity, we set K3D to a constant if  C     For the Redi diffusivity, we set K3D to a constant if
1125  C     GM_K3D_constRedi=.TRUE.  C     GM_K3D_constRedi=.TRUE.
1126        IF (GM_K3D_constRedi) THEN        IF (GM_K3D_constRedi) THEN
1127          DO k=1,Nr          DO k=1,Nr
# Line 1138  C     GM_K3D_constRedi=.TRUE. Line 1134  C     GM_K3D_constRedi=.TRUE.
1134        ENDIF        ENDIF
1135    
1136  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
1137        IF ( useDiagnostics )        IF ( useDiagnostics )
1138       &     CALL DIAGNOSTICS_FILL(K3D, 'GM_K3D_T',0,Nr,0,1,1,myThid)       &     CALL DIAGNOSTICS_FILL(K3D, 'GM_K3D_T',0,Nr,1,bi,bj,myThid)
1139  #endif  #endif
1140    
1141  #endif /* GM_K3D */  #endif /* GM_K3D */
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22