c23456789012345678901234567890123456789012345678901234567890123456789012


      SUBROUTINE x_prior(opt) 

c     ******************************************************************
c
c     This routine sets an a priori estimate of the state xo
c
c     ******************************************************************

c     -------------
c     Declarations
c     -------------

c     ''''''''''''''''''''''
      implicit none
      include 'commonvar.f'      
c     ''''''''''''''''''''''

c---- Argument variables

      character	opt*5


c---- Local variables

      integer	i,j,k,kcomp,ic,n,n_sel

      real*8	conse_m,
     >		therm_U,
     >		therm_V

      real*8	fac

      character*12 stringx

      logical	flag 




c     ---------------------------------
c     Initialization of state variables
c     ---------------------------------

      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1
         r(i,j,k) = 0.d0 
         U(i,j,k) = 0.d0 
         V(i,j,k) = 0.d0 
         W(i,j,k) = 0.d0 
         p(i,j,k) = 0.d0
      enddo
      enddo 
      enddo




c     --------------------------------------------
c     Notes
c     -----
c     1) r-field from WOCE climatology  
c     2) U-field from thermal wind or observations
c     3) V-field from thermal wind or observations 
c     4) W-field from mass conservation
c     5) p-field from WOCE climatology
c     ---------------------------------------------




c     -------------
c     Density field
c     -------------

c---- Read - anomalies of in situ density   r(,,)
c          - vertical eddy density fluxes F_r(,,) 

c     --------------------
      call wocexxx('rval') 
c     --------------------


c---- Include compression effects
c
c      do i=1,itot+1
c      do j=1,jtot+1
c      do k=1,ktot+1
c         if (flag_r(i,j,k)) r(i,j,k) =
c     >			    r(i,j,k) - compr*z_l(k) 
c      enddo
c      enddo
c      enddo


c---- Density anomaly

      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1
         if (flag_r(i,j,k)) r(i,j,k) =
     >                      r(i,j,k) - rho_0
      enddo
      enddo
      enddo




c     ---------------------
c     Transport component U 
c     ---------------------

c---- Level of no motion at 1000 m

c++
      if (lnm.eq.'1p00km') then
c++

c     In progress: Algo below requires dz(4) in therm_U
 

      do   i=1,itot+1
      do 1 j=1,jtot

         if (abs(lat(j  )*deg).le.5.d0) goto 1 
         if (abs(lat(j+1)*deg).le.5.d0) goto 1 

         fac = 1./z_l(ktot+1) + 1./dz(ktot)
c        -----------------------------------------------------
         U(i,j,3) =  -therm_U(i,j,4)                   / fac
         U(i,j,2) = (-therm_U(i,j,3) + U(i,j,3)/dz(3)) * dz(2) 
         U(i,j,1) = (-therm_U(i,j,2) + U(i,j,2)/dz(2)) * dz(1)
c        -----------------------------------------------------

         if (.not.flag_U(i,j,1)) then
           U(i,j,1) = 0.d0
         endif

         if (.not.flag_U(i,j,2)) then
           U(i,j,1) = 0.d0
           U(i,j,2) = 0.d0
         endif

         if (.not.flag_U(i,j,3)) then
           U(i,j,1) = 0.d0
           U(i,j,2) = 0.d0
           U(i,j,3) = 0.d0
         endif

1     continue
      enddo

c++
      endif
c++


c---- Level of no motion at 3750 m

c++
      if (lnm.eq.'3p75km') then
c++

      do   i=1,itot+1
      do 2 j=1,jtot

         if (abs(lat(j  )*deg).le.5.d0) goto 2 
         if (abs(lat(j+1)*deg).le.5.d0) goto 2 

         fac = 1./dz(1) + 1./dz(2)
c        ---------------------------------------------------
         U(i,j,1) = -therm_U(i,j,2) / fac
         U(i,j,2) =  therm_U(i,j,2) / fac
         U(i,j,3) = (therm_U(i,j,3) + U(i,j,2)/dz(2)) * dz(3)
c        ----------------------------------------------------

         if (.not.flag_U(i,j,1)) then
           U(i,j,1) = 0.d0
         endif

         if (.not.flag_U(i,j,2)) then
           U(i,j,1) = 0.d0
           U(i,j,2) = 0.d0
           fac      = 1./dz(2) + 1./dz(3)
c          -------------------------------
           U(i,j,3) = therm_U(i,j,3) / fac 
c          -------------------------------
         endif

         if (.not.flag_U(i,j,3)) then
           U(i,j,1) = 0.d0
           U(i,j,2) = 0.d0
           U(i,j,3) = 0.d0 
         endif

2     continue
      enddo

c++
      endif
c++


c---- Local volume transports

c     -----------------
      call observa(opt)
c     -----------------




c     ---------------------
c     Transport component V 
c     ---------------------

c---- Level of no motion at 1000 m

c++
      if (lnm.eq.'1p00km') then
c++

c     In progress: Algo below requires dz4) in therm_V

      do   i=1,itot
      do 3 j=1,jtot+1

         if (abs(lat(j)*deg).le.5.d0) goto 3 

         fac = 1./z_l(ktot+1) + 1./dz(ktot)
c        -----------------------------------------------------
         V(i,j,3) =  -therm_V(i,j,4)                   / fac
         V(i,j,2) = (-therm_V(i,j,3) + V(i,j,3)/dz(3)) * dz(2) 
         V(i,j,1) = (-therm_V(i,j,2) + V(i,j,2)/dz(2)) * dz(1)
c        -----------------------------------------------------

         if (.not.flag_V(i,j,1)) then
           V(i,j,1) = 0.d0
         endif

         if (.not.flag_V(i,j,2)) then
           V(i,j,1) = 0.d0
           V(i,j,2) = 0.d0
         endif

         if (.not.flag_V(i,j,3)) then
           V(i,j,1) = 0.d0
           V(i,j,2) = 0.d0
           V(i,j,3) = 0.d0
         endif

3     continue
      enddo

c++
      endif
c++


c---- Level of no motion at 3750 m 

c++
      if (lnm.eq.'3p75km') then
c++

      do   i=1,itot
      do 4 j=1,jtot+1

         if (abs(lat(j)*deg).le.5.d0) goto 4 

         fac = 1./dz(1) + 1./dz(2) 
c        ----------------------------------------------------
         V(i,j,1) = -therm_V(i,j,2) / fac
         V(i,j,2) =  therm_V(i,j,2) / fac
         V(i,j,3) = (therm_V(i,j,3) + V(i,j,2)/dz(2)) * dz(3)
c        ----------------------------------------------------

         if (.not.flag_V(i,j,1)) then
           V(i,j,1) = 0.d0
         endif

         if (.not.flag_V(i,j,2)) then
           V(i,j,1) = 0.d0
           V(i,j,2) = 0.d0
           fac      = 1./dz(2) + 1./dz(3)
c          -------------------------------
           V(i,j,3) = therm_V(i,j,3) / fac 
c          -------------------------------
         endif

         if (.not.flag_V(i,j,3)) then
           V(i,j,1) = 0.d0
           V(i,j,2) = 0.d0
           V(i,j,3) = 0.d0 
         endif

4     continue
      enddo

c++
      endif
c++


c---- Local volume transports

c     -----------------
      call observa(opt)
c     -----------------




c     ---------------------
c     Transport component W
c     ---------------------

      do i=1,itot
      do j=1,jtot

         W(i,j,1) = 0.d0
         do k=1+1,ktot+1
c           -------------------------------
            W(i,j,k) = conse_m('W',i,j,k-1)
c           -------------------------------
            if (.not.flag_W(i,j,k)) W(i,j,k) = 0.d0
         enddo

      enddo
      enddo          




c     ---------
c     PO4 field
c     ---------

c---- WOCE

c++
      if (iflag.eq.1) then
c++

c       --------------------
        call wocexxx('pval') 
c       --------------------

c++
      endif
c++




c     -------------------
c     End of the routine
c     -------------------

      return

      end

c     ******************************************************************


      REAL*8 FUNCTION rr(t,s,p)

c     ******************************************************************
c
c     This function calculates the density of seawater as a function of
c     its in situ temperature, practical scale salinity, and pressure
c     (Ref: UNESCO, 1981/1983) 
c
c     ******************************************************************

c     -------------
c     Declarations
c     -------------

c     ''''''''''''''
      implicit none
c     ''''''''''''''

c---- Argument variables

      real*8    t,s,p


c---- Internal variables

      real*8    a0,a1,a2,a3,a4,a5,
     >          b0,b1,b2,b3,b4,
     >          c0,c1,c2,
     >          d0

      real*8    f0,f1,f2,f3,
     >          g0,g1,g2,
     >          i0,i1,i2,
     >          j0,
     >          m0,m1,m2,
     >          e0,e1,e2,e3,e4,
     >          h0,h1,h2,h3,
     >          k0,k1,k2

      real*8    rW,AW,BW,KW,B,A,Kst0,Kstp

      real*8    T1,T2,T3




c     ------------------------
c     Polynomial coefficients
c     ------------------------

      a0 = 999.842594
      a1 =   6.793952d-2
      a2 =  -9.095290d-3
      a3 =   1.001685d-4
      a4 =  -1.120083d-6
      a5 =   6.536332d-9

      b0 = 8.24493d-1
      b1 =-4.0899d-3
      b2 = 7.6438d-5
      b3 =-8.2467d-7
      b4 = 5.3875d-9

      c0 =-5.72466d-3
      c1 = 1.0227d-4
      c2 =-1.6546d-6

      d0 = 4.8314d-4

      f0 = 54.6746d0
      f1 = -0.603459
      f2 =  1.09987d-2
      f3 = -6.1670d-5

      g0 = 7.944d-2
      g1 = 1.6483d-2
      g2 =-5.3009d-4

      i0 = 2.2838d-3
      i1 =-1.0981d-5
      i2 =-1.6078d-6

      j0 = 1.91075d-4

      m0 = -9.9348d-7
      m1 = 2.0816d-8
      m2 = 9.1697d-10

      e0 = 19652.21d0
      e1 =   148.4206d0
      e2 =    -2.327105
      e3 =     1.360477d-2
      e4 =    -5.155288d-5

      h0 = 3.239908d0
      h1 = 1.43713d-3
      h2 = 1.16092d-4
      h3 =-5.77905d-7

      k0 = 8.50935d-5
      k1 =-6.12293d-6
      k2 = 5.2787d-8




c     ----------------------------
c     Density of sea water sample 
c     ----------------------------

c---- Density of reference pure water at atmospheric pressure

      rW = a0      +
     >     a1*t    +
     >     a2*t**2 +
     >     a3*t**3 +
     >     a4*t**4 +
     >     a5*t**5


c---- Density of sea water at atmospheric pressure

      T1 = b0      +
     >     b1*t    +
     >     b2*t**2 +
     >     b3*t**3 +
     >     b4*t**4

      T2 = c0     +
     >     c1*t   +
     >     c2*t**2

      T3 = d0

      rr = rW + T1*s + T2*s**(3.d0/2.d0) + T3*s**2


c---- Secant bulk modulus of reference pure water at atmospheric pressure

      KW = e0      + 
     >     e1*t    + 
     >     e2*t**2 + 
     >     e3*t**3 + 
     >     e4*t**4


c---- Secant bulk modulus of sea water at atmospheric pressure

      Kst0 = 
     >KW   + (f0+f1*t+f2*t**2+f3*t**3)*s 
     >     + (g0+g1*t+g2*t**2        )*s**(3.d0/2.d0)


c---- Secant bulk modulus of sea water

      AW = h0 + h1*t + h2*t**2 + h3*t**3
      BW = k0 + k1*t + k2*t**2

      A  = AW + (i0+i1*t+i2*t**2)*s + j0*s**(3.d0/2.d0) 
      B  = BW + (m0+m1*t+m2*t**2)*s

      Kstp = Kst0 + A*p + B*p**2


c---- Density of sea water at in situ pressure

      rr = rr / (1.d0-p/Kstp)




c     --------------------
c     End of the function 
c     --------------------

      return

      end

c     ******************************************************************


      SUBROUTINE observa(opt) 

c     ******************************************************************
c
c     This routine assigns - the prior value of U,V according to 
c                            an observational estimate 
c                          - the error in this prior value
c
c     ******************************************************************

c     -------------
c     Declarations
c     -------------

c     ''''''''''''''''''''''
      implicit none
      include 'commonvar.f'      
c     ''''''''''''''''''''''

c---- Argument variable

      character	opt*5


c---- Local variables

      integer	n,i,j,k,kco2

      character*12 stringx

      integer	kcomp,j0(jtot+1)
 
      real*8	V_ave





c     --------------------------------------------
c     Local volume transports: 3-degree resolution
c     --------------------------------------------

c---- Test

c++
      if (opt.eq.'AT003') then
c++


c---- WBC: Transport component U

      do i=1,itot+1
      do j=1,jtot
      do k=1,ktot

         do n=1,n_wbc
            if (stringx('UU',i,j,k).eq.str_wbc(n)) U(i,j,k) = b_wbc(n)
         enddo

      enddo
      enddo
      enddo


c---- WBC: Transport component V 

      do i=1,itot
      do j=1,jtot+1
      do k=1,ktot

         do n=1,n_wbc
            if (stringx('VV',i,j,k).eq.str_wbc(n)) V(i,j,k) = b_wbc(n)
         enddo

      enddo
      enddo
      enddo


c---- MOW

      U(23,28,3) = b_mow0


c---- AABW through Vema Channel

      U(14,7,1) = b_vem0


c---- AABW through Hunter Channel

      V(16,5,1) = b_hun0


c---- AABW through Romanche & Chain FZ

      V(19,16,1) = b_rac0


c---- AABW through Equatorial Channel at 36W

      U(15,15,1) = b_eqc0


c---- AABW through Vema FZ

      V(12,19,1) = b_vfz0


c---- End of test

c++
      endif
c++




c     ------------------------------------------------------
c     Errors in local volume transports: 3-degree resolution
c     ------------------------------------------------------

c---- Test

c++
      if (opt.eq.'AT003') then
c++


c---- WBC: Transport component U

      do i=1,itot+1
      do j=1,jtot
      do k=1,ktot

         do n=1,n_wbc
            if (stringx('UU',i,j,k).eq.str_wbc(n)) e_U(i,j,k) = e_wbc0
         enddo

      enddo
      enddo
      enddo


c---- WBC: Transport component V 

      do i=1,itot
      do j=1,jtot+1
      do k=1,ktot

         do n=1,n_wbc
            if (stringx('VV',i,j,k).eq.str_wbc(n)) e_V(i,j,k) = e_wbc0
         enddo

      enddo
      enddo
      enddo


c---- MOW

      e_U(23,28,3) = e_mow0


c---- AABW through Vema Channel

      e_U(14,7,1) = e_vem0


c---- AABW through Hunter Channel

      e_V(16,5,1) = e_hun0


c---- AABW through Romanche & Chain FZ

      e_V(19,16,1) = e_rac0


c---- AABW through Equatorial Channel at 36W

      e_U(15,15,1) = e_eqc0


c---- AABW through Vema FZ

      e_V(12,19,1) = e_vfz0


c---- End of test

c++
      endif
c++




c     --------------------------------------------
c     Local volume transports: 2-degree resolution
c     --------------------------------------------

c---- Test

c++
      if (opt.eq.'AT002') then
c++


c---- WBC: Transport component U

      do i=1,itot+1
      do j=1,jtot
      do k=1,ktot

         do n=1,n_wbc
            if (stringx('UU',i,j,k).eq.str_wbc(n)) U(i,j,k) = b_wbc(n)
         enddo

      enddo
      enddo
      enddo


c---- WBC: Transport component V 

      do i=1,itot
      do j=1,jtot+1
      do k=1,ktot

         do n=1,n_wbc
            if (stringx('VV',i,j,k).eq.str_wbc(n)) V(i,j,k) = b_wbc(n)
         enddo

      enddo
      enddo
      enddo


c---- MOW

      U(35,42,3) = b_mow0

c      U(35,jtot-14,ktot) = b_mow0


c---- AABW through Vema Channel

      V(20,10,1) = b_vem0


c---- AABW through Hunter Channel

      V(25,7,1) = b_hun0


c---- AABW through Romanche & Chain FZ

      V(29,24,1) = b_rac0


c---- AABW through Equatorial Channel at 36W

      U(21,24,1) = b_eqc0


c---- AABW through Vema FZ

      U(18,29,1) = b_vfz0


c---- End of test

c++
      endif
c++




c     ------------------------------------------------------
c     Errors in local volume transports: 2-degree resolution
c     ------------------------------------------------------

c---- Test

c++
      if (opt.eq.'AT002') then
c++


c---- WBC: Transport component U

      do i=1,itot+1
      do j=1,jtot
      do k=1,ktot

         do n=1,n_wbc
            if (stringx('UU',i,j,k).eq.str_wbc(n)) e_U(i,j,k) = e_wbc0
         enddo

      enddo
      enddo
      enddo


c---- WBC: Transport component V 

      do i=1,itot
      do j=1,jtot+1
      do k=1,ktot

         do n=1,n_wbc
            if (stringx('VV',i,j,k).eq.str_wbc(n)) e_V(i,j,k) = e_wbc0
         enddo

      enddo
      enddo
      enddo


c---- MOW

      e_U(35,42,3) = e_mow0

c      e_U(35,jtot-14,ktot) = e_mow0


c---- AABW through Vema Channel

      e_V(20,10,1) = e_vem0


c---- AABW through Hunter Channel

      e_V(25,7,1) = e_hun0


c---- AABW through Romanche & Chain FZ

      e_V(29,24,1) = e_rac0


c---- AABW through Equatorial Channel at 36W

      e_U(21,24,1) = e_eqc0


c---- AABW through Vema FZ

      e_U(18,29,1) = e_vfz0


c---- End of test

c++
      endif
c++




c     ---------------------
c     End of the subroutine
c     ---------------------

      return

      end

c     ******************************************************************
