/[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.14 by m_bates, Wed Jan 1 23:20:48 2014 UTC revision 1.15 by m_bates, Thu Mar 27 05:47:49 2014 UTC
# Line 150  C     surfkz   :: Depth of surface layer Line 150  C     surfkz   :: Depth of surface layer
150        _RL fCoriV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL fCoriV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
151        _RL surfkz(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL surfkz(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
152    
153          _RL centreX, centreY
154          _RL numerator, denominator
155          _RL uInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
156          _RL vInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
157          _RL KdqdxInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
158          _RL KdqdyInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
159          _RL uKdqdyInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
160          _RL vKdqdxInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
161          _RL uXiyInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
162          _RL vXixInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
163          _RL Renorm(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
164          _RL RenormU(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
165          _RL RenormV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
166    
167  C     gradqx   :: Potential vorticity gradient in x direction  C     gradqx   :: Potential vorticity gradient in x direction
168  C     gradqy   :: Potential vorticity gradient in y direction  C     gradqy   :: Potential vorticity gradient in y direction
169  C     XimX     :: Vertical integral of phi_m*K*gradqx  C     XimX     :: Vertical integral of phi_m*K*gradqx
# Line 881  C     ENDIF (.NOT. GM_K3D_smooth) Line 895  C     ENDIF (.NOT. GM_K3D_smooth)
895        ENDIF        ENDIF
896    
897    
898    C     Calculate the renormalisation factor
899          DO j=1-Oly,sNy+Oly
900           DO i=1-Olx,sNx+Olx
901            uInt(i,j)=zeroRL
902            vInt(i,j)=zeroRL
903            KdqdyInt(i,j)=zeroRL
904            KdqdxInt(i,j)=zeroRL
905            uKdqdyInt(i,j)=zeroRL
906            vKdqdxInt(i,j)=zeroRL
907            uXiyInt(i,j)=zeroRL
908            vXixInt(i,j)=zeroRL
909            Renorm(i,j)=zeroRL
910           ENDDO
911          ENDDO
912          DO k=1,Nr
913           DO j=1-Oly,sNy+Oly-1
914            DO i=1-Olx,sNx+Olx-1
915             centreX = op5*(uVel(i,j,k,bi,bj)+uVel(i+1,j,k,bi,bj))
916             centreY = op5*(Kdqdy(i,j,k)     +Kdqdy(i,j+1,k)     )
917    C        For the numerator
918             uInt(i,j) = uInt(i,j)
919         &        + centreX*hfacC(i,j,k,bi,bj)*drF(k)
920             KdqdyInt(i,j) = KdqdyInt(i,j)
921         &        + centreY*hfacC(i,j,k,bi,bj)*drF(k)
922             uKdqdyInt(i,j) = uKdqdyInt(i,j)
923         &        + centreX*centreY*hfacC(i,j,k,bi,bj)*drF(k)
924    C        For the denominator
925             centreY = op5*(Xiy(i,j,k) + Xiy(i,j+1,k))
926             uXiyInt(i,j) = uXiyInt(i,j)
927         &        + centreX*centreY*hfacC(i,j,k,bi,bj)*drF(k)
928    
929             centreX = op5*(Kdqdx(i,j,k)     +Kdqdx(i+1,j,k))
930             centreY = op5*(vVel(i,j,k,bi,bj)+vVel(i,j+1,k,bi,bj) )
931    C        For the numerator
932             vInt(i,j) = vInt(i,j)
933         &        + centreY*hfacC(i,j,k,bi,bj)*drF(k)
934             KdqdxInt(i,j) = KdqdxInt(i,j)
935         &        + CentreX*hfacC(i,j,k,bi,bj)*drF(k)
936             vKdqdxInt(i,j) = vKdqdxInt(i,j)
937         &        + centreY*centreX*hfacC(i,j,k,bi,bj)*drF(k)
938    C        For the denominator
939             centreX = op5*(Xix(i,j,k) + Xix(i+1,j,k))
940             vXixInt(i,j) = vXixInt(i,j)
941         &        + centreY*centreX*hfacC(i,j,k,bi,bj)*drF(k)
942    
943            ENDDO
944           ENDDO
945          ENDDO
946    
947          DO j=1-Oly,sNy+Oly-1
948           DO i=1-Olx,sNx+Olx-1
949            IF (kLowC(i,j,bi,bj).GT.0) THEN
950              numerator =
951         &         (uKdqdyInt(i,j)-uInt(i,j)*KdqdyInt(i,j)/R_low(i,j,bi,bj))
952         &        -(vKdqdxInt(i,j)-vInt(i,j)*KdqdxInt(i,j)/R_low(i,j,bi,bj))
953              denominator = uXiyInt(i,j) - vXixInt(i,j)
954    C         We can have troubles with floating point exceptions if the denominator
955    C         of the renormalisation if the ocean is resting (e.g. intial conditions).
956    C         So we make the renormalisation factor zero if the denominator is very small
957              IF (denominator.LE.small) THEN
958                Renorm(i,j)=zeroRL
959              ELSE
960                Renorm(i,j) = ABS(numerator/denominator)
961              ENDIF
962            ENDIF
963           ENDDO
964          ENDDO
965    C     Now put it back on to the velocity grids
966          DO j=1-Oly+1,sNy+Oly-1
967           DO i=1-Olx+1,sNx+Olx-1
968            RenormU(i,j) = op5*(Renorm(i-1,j)+Renorm(i,j))
969            RenormV(i,j) = op5*(Renorm(i,j-1)+Renorm(i,j))
970           ENDDO
971          ENDDO
972    
973  C     Calculate the eddy induced velocity in the X direction at the west face  C     Calculate the eddy induced velocity in the X direction at the west face
974        DO k=1,Nr        DO k=1,Nr
975         DO j=1-Oly+1,sNy+Oly         DO j=1-Oly+1,sNy+Oly
976          DO i=1-Olx+1,sNx+Olx          DO i=1-Olx+1,sNx+Olx
977           ustar(i,j,k) = -Xix(i,j,k)/coriU(i,j)           ustar(i,j,k) = -Renorm(i,j)*Xix(i,j,k)/coriU(i,j)
978          ENDDO          ENDDO
979         ENDDO         ENDDO
980        ENDDO              ENDDO      
# Line 894  C     Calculate the eddy induced velocit Line 983  C     Calculate the eddy induced velocit
983        DO k=1,Nr        DO k=1,Nr
984         DO j=1-Oly+1,sNy+Oly         DO j=1-Oly+1,sNy+Oly
985          DO i=1-Olx+1,sNx+Olx          DO i=1-Olx+1,sNx+Olx
986           vstar(i,j,k) = -Xiy(i,j,k)/coriV(i,j)           vstar(i,j,k) = -Renorm(i,j)*Xiy(i,j,k)/coriV(i,j)
987          ENDDO          ENDDO
988         ENDDO         ENDDO
989        ENDDO              ENDDO      
# Line 977  C     Diagnostics Line 1066  C     Diagnostics
1066          CALL DIAGNOSTICS_FILL(N2loc,  'GM_N2   ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(N2loc,  'GM_N2   ',0,Nr,0,1,1,myThid)
1067          CALL DIAGNOSTICS_FILL(M4onN2, 'GM_M4_N2',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(M4onN2, 'GM_M4_N2',0,Nr,0,1,1,myThid)
1068          CALL DIAGNOSTICS_FILL(slopeC, 'GM_SLOPE',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(slopeC, 'GM_SLOPE',0,Nr,0,1,1,myThid)
1069            CALL DIAGNOSTICS_FILL(Renorm, 'GM_RENRM',0, 1,0,1,1,myThid)
1070    
1071        ENDIF        ENDIF
1072  #endif  #endif
Legend:
Removed from v.1.14  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.22