C-----------------------------------------------------------------------
Copyright (1999,2005) K. Bergeron, E.A. Coutsias, A.H. Nielsen, J.P. Lynov
C Research Software; use freely - we make no guarantees!
C-----------------------------------------------------------------------
      subroutine init_field(omega,psi,m_act,n_act,lu_log)

      implicit none
        
      include 'limits.h' 
      include 'indata.h'

        integer i, j, ier, m_act, n_act, lu_log, lu_oldbc, k

        parameter(lu_oldbc    = 32)

        double precision lambda, r0, psi0, theta0, a1, c0, epsilon
        double precision factor, costheta
        double precision omega(m_max+2,n_max+2)
        double precision psi(m_max+2,n_max+2)
        double precision disp, sqrtlam, bessel0, bes0, bes1, gamma11
        double precision r, ran2, seed, amp, b, eps, mode
        double precision u, lg, ro, d, po, gamma, co, omegao
        double precision bes01, rcos

	double precision x0, y0, old_x, old_y, new_x, new_y, delta
	double precision new_theta, new_r
   
        double precision aa, pin, ee, ddd, sum, phase(n_act_max)
        integer          per, msize, fcase
        character*60 s


         external disp, ran2

         theta0 = pi/dble(period)
 
c.. put a asymmetric dipole in a annulus 
         if (initcase .eq. 1) then
           u  = init_parval(1)
           lg = init_parval(2)
           ro = init_parval(3)
           d  = init_parval(4)
           po = init_parval(5)
           po = pi*po/180
        
           gamma11 = 3.83170597020751d0
           lambda  = gamma11/ro
           call besj0( gamma11 , bes0)

           co      = 2.0d0*u/lg
           gamma   = -pi*co*ro**2
           omegao  = -gamma/(4*pi*aspect)
           amp     = 2.0d0*u*lambda/bes0

           do i=1,m+1
             do j=1,n+1
               rcos = (lg+d+(x(i)+aspect)*cos(theta(j)))*cos(po)+
     #                (x(i)+aspect)*sin(theta(j))*sin(po)
               r = sqrt((x(i)+aspect)**2 + (lg+d)**2 + 
     #                  2*(x(i)+aspect)*(lg+d)*cos(theta(j)) 
     #                   + lg**2 - 2.0d0*lg*rcos )
               if (r .ge. ro) then
                  omega(i,j) = omegao
               else
                  call besj0( lambda*r , bes01)
                  omega(i,j) = omegao + co*(1.0d0 - bes01/bes0) 
                  if (r .ge. 0.0001) then
                    call besj1(lambda*r,bes1)
                    omega(i,j) = omega(i,j) - amp*bes1*(lg - rcos)/r
                  endif
               endif
            end do
          end do
         end if


c    put a double-vortex in the middle of the flow.
c...      init_parval(1)   :  amplituede of vortex (velocity)
c...      init_parval(2)   :  Not used
c...      init_parval(3)   :  diameter of vortex
c...      init_parval(4)   :  Asymetric dipole, cos(theta)->cos(theta)+epsilon
         if ( initcase .eq. 2) then 

          amp     = init_parval(1)
          eps     = init_parval(3)
          epsilon = init_parval(4) 

          gamma11 = 3.83170597020751d0
          call besj0(gamma11,bes0)
          amp = amp *2.0d0 * gamma11 /(eps*bes0)
          do i=1 , m+1
             do j=1 , n+1
               r = sqrt( (x(i)+aspect)**2 + (aspect)**2 - 
     #                  2.0d0*(aspect)*(x(i)+aspect)* cos(theta(j)-theta0) )
               if (r .gt. eps .or. r.eq.0.d0) then
                  omega(i,j) = 0.0d0
               else 
                  call besj1(gamma11*r/eps , bes1)
                  omega(i,j) = amp*bes1*(x(i)+aspect)*
     -            (cos(theta(j)-theta0-pi/2.0d0)-epsilon)/r 
               endif
             end do
          end do


        endif

c    put a double-vortex in the middle of the flow
c...      init_parval(1)   :  amplituede of vortex (velocity)
c...      init_parval(2)   :  impact on inner cylindre
c...      init_parval(3)   :  diameter of vortex
         if ( initcase .eq. 3) then 

          amp     = init_parval(1)
          delta   = asin(init_parval(2)/aspect)
	  eps     = init_parval(3)
          epsilon = init_parval(4) 

          gamma11 = 3.83170597020751d0
          call besj0(gamma11,bes0)
          amp = -amp *2.0d0 * gamma11 /(eps*bes0)
          do i=1 , m+1
             do j=1 , n+1

	       x0    = aspect*dcos(theta0)
	       y0    = aspect*dsin(theta0)

               old_x = (aspect + x(i))*dcos(theta(j))
               old_y = (aspect + x(i))*dsin(theta(j))

	       new_x = (old_x - x0)*dcos(theta0) + (old_y - y0)*dsin(theta0)
               new_y =-(old_x - x0)*dsin(theta0) + (old_y - y0)*dcos(theta0)

	       if (new_x .ge. 1d-20) then
                 new_theta = atan(new_y/new_x)
               else
                 if (new_x .le. -1d-20) then
                   new_theta = atan(new_y/new_x) + pi
                 else
                   new_theta = pi/2.0d0
                 endif
               endif

               new_r = dsqrt(new_x**2 + new_y**2)

               if (new_r .gt. eps .or. new_r.eq.0.d0) then
                  omega(i,j) = 0.0d0
               else 
                  call besj1(gamma11*new_r/eps , bes1)
                  omega(i,j) = amp*bes1*dsin(new_theta - theta0 - delta)
               endif



             end do
          end do


        endif


         if (initcase .eq. 4) then
          seed = 0.314152965d0
          amp = init_parval(1)
          b   = init_parval(2)
          eps = init_parval(3)
          mode= init_parval(4)
           do i=1,m+1
             do j=1,n+1
                omega(i,j) = -b*(eps/(x(i)+aspect) -
     -                        2.0d0*eps*eps*dtanh(eps*x(i)))/
     -                    (2.0d0*(dcosh(eps*x(i))**2)) *
     -                    (1 + amp*0.02d0*(1-ran2(seed)))*cos(mode*theta(j))
             end do
           end do
         end if

c***********************************************
         if (initcase .eq. 5) then
           seed = 0.314152965d0
           amp = init_parval(1)
           b   = init_parval(2)
           eps = init_parval(3)
           do i=1,m+1
             do j=1,n+1
                omega(i,j) = -b*(4.0d0*x(i)**2*eps**2 - 4.0d0*eps +
     -                    2.0d0*aspect*eps/(x(i)+aspect) - amp**2) *
     -                    dexp(-x(i)**2*eps)*dcos(amp*theta(j)) *
     -                    (1 + amp*0.02d0*(1-ran2(seed)))

             end do
          end do
         end if

         if (initcase .eq. 6) then
           amp = init_parval(1)
           b   = init_parval(2)
           seed = 0.314152965d0
           open(unit=33, file = 'data/jac.dat', form='formatted',
     _          status = 'old')
               read(33,*) aa,per,ee, pin,pin,ddd,fcase,msize
               ekman_coef = 8.d0*visc/ee/ee
               force_coef = ekman_coef
               if (aa .ne. aspect) stop 'a in in wrong'
               if (per .ne. period) stop 'period in in is wrong'
               if (msize .lt. m_act) stop 'too few coef. in data-file'
               do i = 1, m_act
                 read(33,*) pin, pin, omega(i,1), pin
                 do j = 2, n_act
                   omega(i,j) = 0.0d0
                 end do
               end do

           close(unit = 33, status = 'keep')
	   open(unit=lu_log,file=logfile,status='old',access='append',
     -          form='formatted')
           write(lu_log,*) 'initial condition: shear flow'
           write(lu_log,*) 'ekman_coef   =',ekman_coef
           write(lu_log,*) 'force_coef   =',force_coef
           close(unit=lu_log)
           call cosfou(omega,m,n,-1)

c... find max of the setup
           pin = 0.0d0
           do i=1,m+1
            do j=1,n
              if (pin .le. abs(omega(i,j))) then
                pin = abs(omega(i,j))
              end if
            end do
           end do

c... initialize a random phase for the dif. pertubated modes
           do i=1,b
              phase(i) = 2*pi*ran2(seed)
           end do

           do i=1,m+1
            do j=1,n
              sum = 0.0d0
              do k=1,b
                sum = sum + amp*dcos(dble(k)*theta(j)+phase(k))
              end do
              omega(i,j) = omega(i,j)*(1.0d0 + sum/pin)
            end do
           end do
           

         end if


c... use data from a old_run as a initia condition 
         if (initcase .eq. 7) then
            open(unit=lu_oldbc,file='data/save.999',status='old',
     -           form='unformatted')
             read(lu_oldbc)
             read(lu_oldbc)
             read(lu_oldbc)
             read(lu_oldbc)
             do j=1,n_act
               read(lu_oldbc) (psi(i,j) ,   i=1, m_act)
               read(lu_oldbc)
               read(lu_oldbc)
               read(lu_oldbc)
               read(lu_oldbc)
             end do
             close(unit=lu_oldbc)
             call rsqdiv(omega,psi,m_act,n_act)
             call remmod(psi,m_act,n_act)
             call cosfou(omega,m,n,-1)

         end if

c... Initialyze omega with init_parval(1)
         if (initcase .eq. 99) then

           do i=1,m+1
            do j=1,n
              omega(i,j) =init_parval(1)
            end do
           end do
           
         end if

 

        end


 

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

c**************************************************************************
c                     -----------------------
c          ---------<   p r o j e c t i o n   >---------
c                     -----------------------
c   a subroutine that projects initial vorticity to the hyperplane
c   of vorticities compatible with noslip constraints.
c**************************************************************************
        subroutine projection(omega,psi,boop,m_act,n_act,lu_log,logfile,bc_case)

        implicit none

        include 'limits.h'
        include 'helmholtz.h'

        double precision oldchange, newchange, newnorm
        double precision omega(m_max+2,n_max+2)
        double precision psi(m_max+2,n_max+2)
        double precision boop(m_max+2,n_max+2)
        double precision ipcp, ipcm,ipsp, ipsm, np, nm
        double precision pps, pms, ppc, pmc
        integer*4        i, j, k, n_act, m_act, l, lu_log, bc_case
        double precision ipb, ppb

  
        character*50 logfile

      do i=1,m_act
        do j=1,n_act/2
          psi(i,j) = c2r(i,j)
        end do
      end do

      do i=1,m_act
        do j=1,n_act
          boop(i,j) = omega(i,j)
        end do
      end do

c      if (bc_case .eq. 2 .or. bc_case .eq. 4) then
      do 7 j = 3, n_act, 2
         k = 2*(j/2)+1
	 l = (j/2) + 1
         np = 0.0d0
         ipcp = 0.0d0
         ipsp = 0.0d0
         ipb = 0.0d0
         do 12 i = 1, m_act
            ipcp = ipcp + omega(i,k  ) * c1r(i,l)
            ipsp = ipsp + omega(i,k+1) * c1r(i,l)
            np   = np   + c1r(i,l) * c1r(i,l)
            ipb  = ipb  + c2r(i,l) * c1r(i,l)
c           ipb  = ipb  + psi(i,l) * c1r(i,l)
12       continue
         ppc = ipcp / np
         pps = ipsp / np
         ppb = ipb  / np
         do 22 i = 1, m_act
            omega(i,k  ) = omega(i,k  ) - ppc*c1r(i,l)
            omega(i,k+1) = omega(i,k+1) - pps*c1r(i,l)
            c2r(  i,l  ) = c2r(  i,l  ) - ppb*c1r(i,l)
c           psi(  i,l  ) = psi(  i,l  ) - ppb*c1r(i,l)
22       continue
7     continue
c      end if

c      if (bc_case .eq. 2 .or. bc_case .eq. 5) then
      do 5 j = 3, n_act, 2
         k = 2*(j/2)+1
	 l = (j/2) + 1
         nm = 0.0d0
         ipcm = 0.0d0
         ipsm = 0.0d0
         do 10 i = 1, m_act
            ipcm = ipcm + omega(i,k  ) * c2r(i,l) 
c           ipcm = ipcm + omega(i,k  ) * psi(i,l) 
            ipsm = ipsm + omega(i,k+1) * c2r(i,l)
c           ipsm = ipsm + omega(i,k+1) * psi(i,l)
            nm   = nm   + c2r(  i,l)   * c2r(i,l)
c           nm   = nm   + psi(  i,l)   * psi(i,l)
10       continue
         pmc = ipcm / nm
         pms = ipsm / nm
         do 20 i = 1, m_act
            omega(i,k  ) = omega(i,k  ) - pmc*c2r(i,l)
            omega(i,k+1) = omega(i,k+1) - pms*c2r(i,l)
c           omega(i,k  ) = omega(i,k  ) - pmc*psi(i,l)
c           omega(i,k+1) = omega(i,k+1) - pms*psi(i,l)
20       continue
5     continue
c      end if


c... Culculate chage in omega after projection
      newnorm  = 0.0d0
      do i=1,m_act
        do j=1,n_act
          newnorm = newnorm + (boop(i,j) - omega(i,j))**2 
        end do
      end do

c... Write change in projection to log-file
      open(unit=lu_log,file=logfile,status='old',access='append')
       write(lu_log,*) 'change in omega after projection:',newnorm
      close (unit=lu_log,status='keep')
      end


c+++++++++++++++ Initialyzed forcing +++++++++++++++++++
      subroutine init_force(m_act,n_act,force_array,psi,dy2,lu_log)
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++

      implicit none

      include 'limits.h'
      include 'indata.h'
      include 'outdata.h'

      integer*4         m_act, n_act,lu_log
 
      double precision  force_array(m_act_max,n_act_max)
      double precision  psi(m_max+2,n_max+2)
      double precision  dy2(m_max+2,n_max+2)

        double precision aa, pin, ee, ddd
        integer          per, msize, fcase

      integer*4         i, j


      if (init_force_case .eq. 1) then
c... initialized forcing in real space in array psi
           open(unit=33, file = 'data/jac.dat', form='formatted',
     _          status = 'old')
               read(33,*) aa,per,ee, ee,ee,ddd,fcase,msize
               do i = 1, m_act
                 read(33,*) pin, pin, pin, force_array(i,1)
                 do j = 2, n_act
                   force_array(i,j) = 0.d0
                 end do
               end do
           close(unit = 33, status = 'keep')

        open(unit=lu_log,file=logfile,status='old',access='append',
     -       form='formatted')

 	write(lu_log,*) 'in init_force'
        call rsqdiv(dy2,force_array,m_act,n_act)
        ddd = sur_int(m_act,dy2)
        write(lu_log,*) 'circulation of forcing   =', ddd
        write(lu_log,*) 'p_bc_in    = ',p_bc_in
        write(lu_log,*) 'p_bc_out   = ',p_bc_out
	write(lu_log,*) 'force_coef = ',force_coef
        close(unit=lu_log)

        do 30 i=1,m_act
          do 40 j=1,n_act
            force_array(i,j) = force_coef*force_array(i,j)
 40       continue
 30     continue
       else
         stop 'Illegal init_force_case'
       endif

       end 


c++++++++++++++ Random generator +++++++++++++
        double precision function ran2(seed)
c+++++++++++++++++++++++++++++++++++++++++++++

        double precision seed


        ran2 = seed + dsqrt(3.0d0) - dint(seed + dsqrt(3.0d0))


        seed = ran2

        return
        end 

