integer i, j, k

     &                 j, ishfts, maxitr, mode, nconv

            do j=1, nconv

               print *, 'ph-ev ', j, s(j,1), s(j,2)

            do 20 j=1, nconv

               s(j,1) = sqrt(s(j,1))

                  tmpvec1(l) = v(l,j)

               call DCOPY( m, tmpvec2, 1, u(l,j), 1 )

               temp = one / dnrm2( m, u(l,j), 1 )

               print *, 'ph-print ', j, dnrm2( m, u(l,j), 1 )

               call DSCAL(m , temp, u(l,j), 1 )

               call DAXPY(m, -s(j,1), u(1,j), 1, tmpvec2, 1)

               s(j,2) = dnrm2(m, tmpvec2, 1)

      integer i, j, k

      integer i, j, k

        DO j=1,sNy

     &      + g_objf_state_final(i,j,bi,bj,k+0*Nr)**2/stdDevtheta**2

     &      + g_objf_state_final(i,j,bi,bj,k+1*Nr)**2/stdDevsalt**2

     &      + g_objf_state_final(i,j,bi,bj,k+2*Nr)**2/stdDevuvel**2

     &      + g_objf_state_final(i,j,bi,bj,k+3*Nr)**2/stdDevvvel**2

     &      + g_objf_state_final(i,j,bi,bj,1+4*Nr)**2/stdDevetan**2

        DO j=1,sNy

            do j = jmin, jmax

                tmpfld3d(i,j,k,bi,bj) = tmpfld3d(i,j,k,bi,bj)

     &                + advvel(i,j,k,bi,bj)

            do j = jmin, jmax

                tmpfld3d(i,j,k,bi,bj) = tmpfld3d(i,j,k,bi,bj)

     &                + aduvel(i,j,k,bi,bj)

            do j = jmin, jmax

                tmpfld3d(i,j,k,bi,bj) = tmpfld3d(i,j,k,bi,bj)

     &                + adsalt(i,j,k,bi,bj)*fac

            do j = jmin, jmax

                tmpfld3d(i,j,k,bi,bj) = tmpfld3d(i,j,k,bi,bj)

     &                + adtheta(i,j,k,bi,bj)*fac

      integer i, j, k

        DO j = jmin, jmax

            adtheta(i,j,k,bi,bj) = 

     &       g_objf_state_final(i,j,bi,bj,k)

            adsalt(i,j,k,bi,bj) = 

     &       g_objf_state_final(i,j,bi,bj,1*Nr+k)

            aduvel(i,j,k,bi,bj) = 

     &       g_objf_state_final(i,j,bi,bj,2*Nr+k)

            advvel(i,j,k,bi,bj) = 

     &       g_objf_state_final(i,j,bi,bj,3*Nr+k)

          adetan(i,j,bi,bj) = 

     &       g_objf_state_final(i,j,bi,bj,4*Nr+1)

          do j = jmin, jmax

              tmpfld2d(i,j,bi,bj) = tmpfld2d(i,j,bi,bj)

     &              + adetan(i,j,bi,bj)

      integer i, j, k

        DO j=1,sNy

            adobjf_state_final(i,j,bi,bj,k) = 

     &       g_objf_state_final(i,j,bi,bj,k)

           adobjf_state_final(i,j,bi,bj,4*Nr+1) = 

     &      g_objf_state_final(i,j,bi,bj,4*Nr+1)

        DO j = 1,sNy

         ReCountX(j,bj) = CountX(j,bj)

         IF (yC(1,j,bi,bj) .LE. -62.0) THEN

            ReCountX(j,bj) = 90.

         ELSE IF (yC(1,j,bi,bj) .EQ. 74.0) THEN

            ReCountX(j,bj) = CountX(j,bj) + 9.0

         ELSE IF (yC(1,j,bi,bj) .EQ. 70.0) THEN

            ReCountX(j,bj) = CountX(j,bj) + 8.0

         ELSE IF (yC(1,j,bi,bj) .EQ. 66.0) THEN

            ReCountX(j,bj) = CountX(j,bj) + 5.0

         ELSE IF (yC(1,j,bi,bj) .EQ. 62.0) THEN

            ReCountX(j,bj) = CountX(j,bj) + 1.0

        DO j=1,sNy

          S(i,j,bj) = sin(yC(i,j,bi,bj)*deg2rad)

          P2(i,j,bj) = 0.5*(3.*S(i,j,bj)**2 - 1.)

          P4(i,j,bj) = 0.12*(35.*S(i,j,bj)**4 - 30.*S(i,j,bj)**2 + 3.)

        DO j=1,sNy

          IF (yC(i,j,bi,bj) .LT. 0.) THEN

          SW(j,bj) = 0.25*Q0*(1 + Q2*P2(i,j,bj))*

     <         (1 - A0 - A2*P2(i,j,bj) - A4*P4(i,j,bj) )

          T(j,bj) = T0(no_so) + T2(no_so)*P2(i,j,bj)

          LW(j,bj) = LW0 + LW1*(T(j,bj)-t_mlt)

             LW(j,bj) = LW(j,bj) - 

         Qnet(i,j,bi,bj) = -1.0*( SW(j,bj) - LW(j,bj) - 

     <        0.000728e4      - 0.00678e4*S(i,j,bj) + 

     <        0.0955e4*S(i,j,bj)**2 + 0.0769e4*S(i,j,bj)**3 - 

     <        0.8508e4*S(i,j,bj)**4 - 0.3581e4*S(i,j,bj)**5 + 

     <        2.9240e4*S(i,j,bj)**6 + 0.8311e4*S(i,j,bj)**7 -   

     <        4.9548e4*S(i,j,bj)**8 - 0.8808e4*S(i,j,bj)**9 + 

     <        4.0644e4*S(i,j,bj)**10 +0.3409e4*S(i,j,bj)**11 - 

     <        1.2893e4*S(i,j,bj)**12 )

          IF (ReCountX(j,bj) .GT. 0.) THEN

             Qnet(i,j,bi,bj) = (-90./ReCountX(j,bj))*

     <            ( SW(j,bj) - LW(j,bj) - 

     <            Hd35(no_so)*(3.054e1 - 3.763e1*S(i,j,bj) + 

     <        1.892e2*S(i,j,bj)**2 + 3.041e2*S(i,j,bj)**3 - 

     <        1.540e3*S(i,j,bj)**4 - 9.586e2*S(i,j,bj)**5 + 

     <        2.939e3*S(i,j,bj)**6 + 1.219e3*S(i,j,bj)**7 -   

     <        2.550e3*S(i,j,bj)**8 - 5.396e2*S(i,j,bj)**9 + 

     <        8.119e2*S(i,j,bj)**10)

             Qnet(i,j,bi,bj) = 0.

          EmPmR(i,j,bi,bj) = -1.e-3*Fw35(no_so)

     <    *(-0.8454e5*S(i,j,bj)**14 + 0.5367e5*S(i,j,bj)**13 

     <    +3.3173e5*S(i,j,bj)**12 - 1.8965e5*S(i,j,bj)**11 

     <    -5.1701e5*S(i,j,bj)**10

     <    +2.6240e5*S(i,j,bj)**9 + 4.077e5*S(i,j,bj)**8 

     <    -1.791e5*S(i,j,bj)**7

     <    -1.7231e5*S(i,j,bj)**6 + 0.6229e5*S(i,j,bj)**5 

     <    +0.3824e5*S(i,j,bj)**4

     <    -0.1017e5*S(i,j,bj)**3 - 0.0387e5*S(i,j,bj)**2 

     <    +0.00562e5*S(i,j,bj)  + 0.0007743e5)

          IF (yC(i,j,bi,bj) .LT. -40.) THEN

           EmPmR(i,j,bi,bj) = -1.e-3*(Fw35(no_so)*

     <            (-6.5 + 35.3 + 71.7*S(i,j,bj)  

     <           - 1336.3*S(i,j,bj)**2 - 425.8*S(i,j,bj)**3 

     <           + 5434.8*S(i,j,bj)**4 + 707.9*S(i,j,bj)**5 

     <           - 6987.7*S(i,j,bj)**6 - 360.4*S(i,j,bj)**7 

     <           + 2855.0*S(i,j,bj)**8)

           IF (xC(i,j,bi,bj) .GT. 284. 

     <      .OR. xC(i,j,bi,bj) .LT. 28.) THEN

              EmPmR(i,j,bi,bj) = -1.e-3*(Fw35(no_so)*

     <             (-6.5 -2.878 + 3.157e2*S(i,j,bj) -

     <             2.388e3*S(i,j,bj)**2 - 4.101e3*S(i,j,bj)**3 + 

     <             1.963e4*S(i,j,bj)**4 + 1.534e4*S(i,j,bj)**5 - 

     <             6.556e4*S(i,j,bj)**6 - 2.478e4*S(i,j,bj)**7 + 

     <             1.083e5*S(i,j,bj)**8 + 1.85e4*S(i,j,bj)**9 - 

     <             8.703e4*S(i,j,bj)**10 - 5.276e3*S(i,j,bj)**11 + 

     <             2.703e4*S(i,j,bj)**12)

              EmPmR(i,j,bi,bj) = -1.e-3*(Fw35(no_so)

     <             *(-6.5 +51.89 + 4.916e2*S(i,j,bj) - 

     <             1.041e3*S(i,j,bj)**2 - 7.546e3*S(i,j,bj)**3 + 

     <             2.335e3*S(i,j,bj)**4 + 3.449e4*S(i,j,bj)**5 + 

     <             6.702e3*S(i,j,bj)**6 - 6.601e4*S(i,j,bj)**7 - 

     <             2.594e4*S(i,j,bj)**8 + 5.652e4*S(i,j,bj)**9 + 

     <             2.738e4*S(i,j,bj)**10 - 1.795e4*S(i,j,bj)**11 - 

     <             9.486e3*S(i,j,bj)**12)

      INTEGER i, j

        DO j=1-oLy,sNy+oLy

          S(i,j,bj) = 0.0

          P2(i,j,bj) = 0.0

          P4(i,j,bj) = 0.0

         SW(j,bj) = 0.0

         LW(j,bj) = 0.0

         Hd(j,bj) = 0.0

         Fw(j,bj) = 0.0

         T(j,bj) = 0.0

         ReCountX(j,bj) = 0.0

      INTEGER i,j

      DO j = jMin, jMax

          surfaceForcingU(i,j,bi,bj) = 

     &      fu(i,j,bi,bj)*horiVertRatio*recip_rhoConst

     &           + winPert(i,j,bi,bj)

     &            *drF(kSurface)*hFacW(i,j,kSurface,bi,bj)

          surfaceForcingV(i,j,bi,bj) = 

     &      fv(i,j,bi,bj)*horiVertRatio*recip_rhoConst

          surfaceForcingT(i,j,bi,bj) = 

     &      -Qnet(i,j,bi,bj)*recip_Cp*horiVertRatio*recip_rhoConst

     &      -lambdaThetaZonRelax*maskC(i,j,kSurface,bi,bj)*

     &         (theta(i,j,kSurface,bi,bj)-ZonalMeanSST(j,bj))

     &        *drF(kSurface)*hFacC(i,j,kSurface,bi,bj)

          surfaceForcingS(i,j,bi,bj) = 

     &      EmPmR(i,j,bi,bj)*convertFW2Salt*convertEmP2rUnit

     &      +Run(i,j,bi,bj)*scale_runoff

     &      -lambdaSaltClimRelax(i,j,bi,bj)*maskC(i,j,kSurface,bi,bj)

     &        *(salt(i,j,kSurface,bi,bj)-SSS(i,j,bi,bj))

     &        *drF(kSurface)*hFacC(i,j,kSurface,bi,bj)

      INTEGER i, j, k, bi, bj

       DO j=1-OLy,sNy+OLy

         ZonalMeanSST(j,bj) = 0.0

         CountX(j,bj) = 0.0

        DO j = 1-OLy, sNy+OLy

         CountX_tile(j,bi,bj) = 0.0

         ZonalMean_tile(j,bi,bj) = 0.0

          ZonalMean_tile(j,bi,bj) = ZonalMean_tile(j,bi,bj) +

     &           theta(i,j,k,bi,bj)

          CountX_tile(j,bi,bj) = CountX_tile(j,bi,bj) +  

     &         maskC(i,j,k,bi,bj)

         ZonalMeanSST(j,bj) = ZonalMeanSST(j,bj) + 

     &    ZonalMean_tile(j,bi,bj)

         CountX(j,bj) = CountX(j,bj) + CountX_tile(j,bi,bj)

       DO j=1-OLy,sNy+OLy

        _GLOBAL_SUM_R8( CountX(j,bj), myThid )

        _GLOBAL_SUM_R8( ZonalMeanSST(j,bj), myThid )

       DO j=1-OLy,sNy+OLy

        IF ( CountX(j,bj) .GT. 0.0) THEN

          ZonalMeanSST(j,bj) = ZonalMeanSST(j,bj)/CountX(j,bj)

       DO j = jMin, jMax

         y(j) = j

         ya(j) = y(j) - MOD(y(j),6.0)

         ya2(j) = ya(j)+6.0

         if ( ya2(j) .gt. sNy+Oly ) then

            ya2(j) = 0.0

         winPert(i,j,bi,bj) = maskW(i,j,k,bi,bj)*

     &        (1.0/1.66)*(0.75*winPert(i,j,bi,bj) + 

     &        stdev(j)*(1.0/6.0)*

     &        ((y(j)-ya(j))*temp2(i,INT(ya2(j)))+

     &        (6.0-y(j)+ya(j))*temp2(i,INT(ya(j)))))

      INTEGER i, j, bi, bj

      DO j = 1-OLy, sNy+OLy

         y(j) = 0.0

         ya(j) = 0.0

         ya2(j) =  0.0

         stdev(j) = 0.0

        DO j = 1-OLy, sNy+OLy

         temp(i,j) = 0.0

         temp2(i,j) = 0.0

      DO j = 1, sNy

         stdev(j) = std(j)

      DO j = jMin, jMax

         temp(i,j) = 1.73*(2.0*temp(i,j) - 1.0)

      DO j = jMin, jMax

         temp2(i,j) = 0.1*( (x(i)-xa(i))*temp(INT(xa2(i)),j)+

     &        (10.0-x(i)+xa(i))*temp(INT(xa(i)),j) )

      INTEGER bi,bj,i,j

       DO j=1,Ny-1

        lY = lY + delY(j)

         DO j=1,sNy

           distY = (yC(i,j,bi,bj)-(yC0))/lY

           fu(i,j,bi,bj) = tauX

         DO j=1,sNy

           fv(i,j,bi,bj) = 0.0

      INTEGER bi,bj,i,j,k

        DO j=1-OLy,sNy+OLy

          Run   (i,j,bi,bj) = 0.

          winPert(i,j,bi,bj) = 0.

# define _recip_hFacS(i,j,k,bi,bj) recip_hFacS(i,j,k,bi,bj)*maskS(i,j,k,bi,bj)

          ELSEIF ( yC(i,j,bi,bj) .GE. lat(2) .AND. 

     &       yC(i,j,bi,bj) .LT. lat(3)) THEN

     &           theta(i,j,k,bi,bj)

     &           maskC(i,j,k,bi,bj)

      INTEGER i, j, k, bi, bj

        DO j = 1,sNy

          IF ( yC(i,j,bi,bj) .GE. -lat(3) .AND. 

     &         yC(i,j,bi,bj) .LE. -lat(2)) THEN

     &           theta(i,j,k,bi,bj)

     &           maskC(i,j,k,bi,bj)

          ELSEIF ( yC(i,j,bi,bj) .GT. -lat(2) .AND. 

     &       yC(i,j,bi,bj) .LT. -lat(1)) THEN 

     &           theta(i,j,k,bi,bj)

     &           maskC(i,j,k,bi,bj)

          ELSEIF ( yC(i,j,bi,bj) .GT. lat(1) .AND. 

     &       yC(i,j,bi,bj) .LT. lat(2)) THEN

     &           theta(i,j,k,bi,bj)

     &           maskC(i,j,k,bi,bj)

          if ( ABS(uwind(i,j,bi,bj)) .GT. 100. .AND.

     &         maskW(i,j,1,bi,bj) .NE. 0. ) then

     &          bi, bj, i, j, myiter, uwind(i,j,bi,bj)