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.)

     <        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 )

     <            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)

     <    *(-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)

     <            (-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)

     <             (-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)

     <             *(-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)

          S(i,j,bj) = 0.0

      _RL S(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSy)

      REAL    SINA, COSA, C, S

               S = SHC(K-M)

               SHC(K+M) = COSA*C + SINA*S

               SHC(K-M) =-SINA*C + COSA*S

      SUBROUTINE HELMHOLTZ(LMAX,S,P)

      REAL    S           ! INPUT  sin(latitude).

      C  = SNGL(DSQRT(1.0D0-DBLE(S)*DBLE(S)))

            P(K) = S*P(K-1)/AK

               P(K) =(S*P(K-1)-P(K-2)*A(M))/AK

      SR = ((((9.1697E-10*T+2.0816E-8)*T-9.9348E-7) * S + 

     &     S + ((-5.77905E-7*T+1.16092E-4)*T+1.43713E-3)*T-0.1194975) 

     &     ((-6.1670E-5*T+1.09987E-2)*T-0.603459)*T+54.6746) * S + 

      _RL DPT,LAT,S,T

     &     ((2.7759E-12*T-1.1351E-10)*(S-35.0)+

     &     1.8741E-8))*P+(-4.2393E-8*T+1.8932E-6)*(S-35.0)+

     &     ((2.7759E-12*T-1.1351E-10)*(S-35.0)+

     &     1.8741E-8))*P+(-4.2393E-8*T+1.8932E-6)*(S-35.0)+

     &     ((2.7759E-12*T-1.1351E-10)*(S-35.0)+

     &     1.8741E-8))*P+(-4.2393E-8*T+1.8932E-6)*(S-35.0)+

      _RL FUNCTION SBO_RHO( DPT, LAT, S, T )

     &     ((2.7759E-12*T-1.1351E-10)*(S-35.0)+

     &     1.8741E-8))*P+(-4.2393E-8*T+1.8932E-6)*(S-35.0)+

      SR = sqrt(abs(S))

      x = (4.8314E-4 * S +

     &     +8.24493E-1)*S + ((((6.536332E-9*T-1.120083E-6)

      _RL S,T,P,PR

      del_th = del_P*sw_adtg(S,T,P)

      del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)

      del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)

      del_th = del_P*sw_adtg(S,th,P+del_P)

      _RL function SW_PTMP  (S,T,P,PR)

      _RL S,T,P

     &     + (b0 + b1*T)*(S-sref)  

     &     + ( (c0 + (c1 + (c2 + c3*T)*T)*T) + (d0 + d1*T)*(S-sref) )*P 

      _RL function SW_ADTG  (S,T,P)