#include "options.inc"
c
c-----------------------------------------------------------------------
c      Templates for hardwired grid, topography, surface boundary
c      and initial concitions.
c-----------------------------------------------------------------------
c

      module eddy1b_module
c      real :: Bv0 = 0.004, g=9.81 
c      real :: fac = -20.0
c      real :: tscl = 1./(86400.0*20.0)
      real :: Bv0 = 0.004, g=9.81 
      real :: fac = -10.0
      real :: tscl = 1./(86400.0*1.0)
      end module eddy1b_module

      subroutine grid_template()
c-----------------------------------------------------------------------
c     set up a grid here
c-----------------------------------------------------------------------
      use spflame_module
      use eddy1b_module
      implicit none
      real :: dx=0.,dy=0.,dz=0.,x0=0.,y0=0.,sumt,sumu
      integer i,k
c     equidistant vertical grid
      integer,parameter :: km_in = 20
      real ddz(km_in)
      data ddz/ 20*100.e2/
c     small eddy resolving region
      x0=280.; y0=30.; dx = 1./3.; dy = dx*cos(y0/180.*pi)

      if (my_pe==0) then
        print*,''
        print*,'setting up a simple grid '
        print*,' dx= ',dx,' dy= ',dy
        print*,''
      endif
      dxtdeg =  dx  ; dytdeg =  dy  ; dxudeg =  dx  ; dyudeg =  dy  
      xt(1)=x0; xu(1)=x0+dxtdeg(1)/2.
      do i=2,imt
        xt(i)=xt(i-1)+dxtdeg(i); xu(i)=xu(i-1)+dxudeg(i)
      enddo
      yt(1)=y0; yu(1)=y0+dytdeg(1)/2.
      do i=2,jmt
        yt(i)=yt(i-1)+dytdeg(i); yu(i)=yu(i-1)+dyudeg(i)
      enddo
      if (km /= km_in ) then
        print*,' km must be ',km_in
        stop
      endif
      dzt=ddz; dzw(1) = dzt(1); sumt = dzt(1)
      sumu = 0.5*dzt(1) + dzw(1)
      do i=2,km
       sumt = sumt + dzt(i)
       dzw(i) = 2*(sumt-sumu)
       sumu = sumu + dzw(i)
      enddo
      zt(1) = 0.5*dzt(1); zw(1) = dzt(1)
      do i=2,km
        zw(i) = zw(i-1) + dzt(i); zt(i) = zt(i-1) + dzw(i-1)
      enddo
      dzw(0)  = zt(1); dzw(km) = zw(km) - zt(km)
      end subroutine grid_template



      subroutine kmt_template
c-----------------------------------------------------------------------
c     set up a simple topography
c-----------------------------------------------------------------------
      use spflame_module
      implicit none
      kmt_big=km
      end subroutine kmt_template



      subroutine htp_template
c-----------------------------------------------------------------------
c     set up a simple topography for partial cells
c-----------------------------------------------------------------------
      use spflame_module
      implicit none
      htp=zw(km)
      end subroutine htp_template



      subroutine init_cond_template
c-----------------------------------------------------------------------
c     setup simple initial conditions
c     for temperature, salinity and velocity
c     Note that streamfunction is diagnosed from full velocities.
c-----------------------------------------------------------------------
      use spflame_module
      use freesurf_module
      use eddy1b_module
      implicit none
      integer :: i,j,k,n
      real :: x,y,z,h0,h_w,T0
c-----------------------------------------------------------------------
c        A mesoscale front
c-----------------------------------------------------------------------
      do j=js_pe,je_pe
       do k=1,km
        do i=is_pe,ie_pe
            t(i,k,j,2,:) = 0.   ! salinity equals 35 PSU
            x=float(i)/float(imt)
            y=float(j)/float(jmt)
            T0=20+y*fac
            t(i,k,j,1,:) = T0+1/(g*alpha_lin)*zt(k)/100.*Bv0**2
     &       + 0.01*sin( x*2.0*2*pi)
            t(i,k,j,1,:) = t(i,k,j,1,:)*tmask(i,k,j)
        enddo 
       enddo
      enddo

      do n=1,nt
       call set_cyclic_4th_order(t(:,:,:,n,0),km,2)
       call set_cyclic_4th_order(t(:,:,:,n,1),km,2)
      enddo
      end subroutine init_cond_template


      subroutine sbc_template
c-----------------------------------------------------------------------
c     setup simple surface boundary conditions
c-----------------------------------------------------------------------
      use spflame_module
      implicit none
c     no wind stress      
      smf=0.
c       no heat or salt fluxes
      stf_clim=0.
      stf_rest=0.
c       no wind stirring
      if (enable_ktmix.or.enable_tkemix) ustar=0.     ! no wind stirring
      end subroutine sbc_template



      subroutine blue_template
c-----------------------------------------------------------------------
c     setup simple temp / salinity arrays for BLUE
c-----------------------------------------------------------------------
      use spflame_module
      use blue_module
      implicit none
      end subroutine blue_template



      subroutine sponge_template(n)
c-----------------------------------------------------------------------
c     setup simple sponge layers. 
c-----------------------------------------------------------------------
      use spflame_module
      use eddy1b_module
      implicit none
      integer, intent(in) :: n ! number of the tracer
      integer :: i,j,k
      real :: x,y,z,h0,h_w,T0
      integer :: spg_width = 3

      if (n==1) then

       do j=js_pe,je_pe
         do k=1,km
          do i=is_pe,ie_pe
           x=float(i)/float(imt)
           y=float(j)/float(jmt)
           T0=20+y*fac+1/(g*alpha_lin)*zt(k)/100.*Bv0**2
           source(i,k,j)=source(i,k,j)+
     &        tscl*(T0-t(i,k,j,n,taum1))*tmask(i,k,j)
          enddo
         enddo
       enddo

      elseif (n==2) then

       if (my_blk_j==1) then
        do j=2,(spg_width+1)
         do k=1,km
          T0=1.0
          do i=is_pe,ie_pe
           source(i,k,j)=source(i,k,j)+
     &         tscl/(j-1.)*(T0-t(i,k,j,n,taum1))*tmask(i,k,j)
     &        *exp( - (zt(k)-500.e2)**2/ (250.e2)**2 )
          enddo
         enddo
        enddo
       endif
       if (my_blk_j==n_pes_j) then
        do j=jmt-1,jmt-spg_width,-1
         do k=1,km
          T0 = -1.0
          do i=is_pe,ie_pe
           source(i,k,j)=source(i,k,j)+
     &          tscl/(-1*(j-jmt))*
     &        (T0-t(i,k,j,n,taum1))*tmask(i,k,j)
     &        *exp( - (zt(k)-1500.e2)**2/ (250.e2)**2 )
          enddo
         enddo
        enddo
       endif


      endif
      end subroutine sponge_template
