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

      subroutine grid_template()
c-----------------------------------------------------------------------
c     set up a grid here
c-----------------------------------------------------------------------
      use spflame_module
      implicit none
      real :: dx=0.,dy=0.,dz=0.,x0=0.,y0=0.,sumt,sumu
      integer i,k

#ifdef notdef
c
c     vertical grid same as in FLAME
c
      logical :: flame_grid = .true.
      integer,parameter :: km_in = 45
      real ddz(km_in)
      data ddz/ 
     &  1000.0,   1044.051, 1102.236, 1169.269, 1246.635, 1336.169,
     &  1440.147, 1561.417, 1703.569, 1871.169, 2070.068, 2307.849,
     &  2594.425, 2992.892, 3518.793, 4100.175, 4863.114, 5680.496,
     &  6513.772, 7038.668, 7397.827, 7758.212, 8130.866, 8726.499,
     &  9576.454,10908.060,12732.562,15187.976,18426.63,
     &  22000.,  24182.902,24817.098, 13*250.00e2 /
#endif

c#ifdef notdef
c
c     equidistant vertical grid
c
      logical :: flame_grid = .false.
      integer,parameter :: km_in = 50
      real ddz(km_in)
      data ddz/ 50*20.e2/
c#endif
c
c     horizontal grid
c
c      midlatidude basin
c      x0=280.; y0=30.; dx = 50./jmt; dy = dx*cos(y0/180.*pi)
c
c     equatorial basin
c      x0=280.; dx=2.;dy=2.; y0=-jmt/2.*dy + dy/2.
c
c     small eddy resolving region
c      x0=280.; y0=30.; dx = 1./3.; dy = dx*cos(y0/180.*pi)
      x0=280.; y0=-15.0; dx = 1./3.; dy = dx

      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 (my_pe==0 .and. flame_grid) then
        print*,''
        print*,' using vertical FLAME grid'
        print*,''
      endif
      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
      integer :: i,j,indp
      real :: hy,hcm

c     flat bottom
      kmt_big=km

#ifdef notdef
      do i=1,imt
       do j=1,jmt
          hy = (j-1.)/(jmt-1.)
          hcm=( zw(km)-zw(km)/2.* exp(- (xt(i)-xt(imt/2+1))**2/2.**2 ))
          kmt_big(i,j) = indp (hcm, zw, km)
       enddo
      enddo
#endif

#ifdef notdef
c     a ridge 
      do i=1,imt
       do j=1,jmt
          hy = (j-1.)/(jmt-1.)
          hcm=( 5500-3000.* exp(- (xt(i)-xt(imt/2+1))**2/15.**2 ) )*100.
          kmt_big(i,j) = indp (hcm, zw, km)
       enddo
      enddo
#endif

#ifdef notdef
c     a  north-south slope
      do i=1,imt
       do j=1,jmt
          hy = (j-1.)/(jmt-1.)
          hcm=(5500.0-5500.0/2*hy)*100.
          kmt_big(i,j) = indp (hcm, zw, km)
       enddo
      enddo
#endif

#ifdef notdef
c       a slope at the western boundary
        kmt_big(2,:)=indp(200*100.,zw,km)
        kmt_big(3,:)=indp(200*100.,zw,km)
        kmt_big(4,:)=indp(1000*100.,zw,km)
        kmt_big(5,:)=indp(2000*100.,zw,km)
        kmt_big(6,:)=indp(3000*100.,zw,km)
        kmt_big(7,:)=indp(4000*100.,zw,km)
#endif

      end subroutine kmt_template




      subroutine htp_template
c-----------------------------------------------------------------------
c     set up a simple topography for partial cells
c-----------------------------------------------------------------------
      use spflame_module
      implicit none
      integer :: i,j,indp
      real :: hy,hcm

       htp=0.
       do i=1,imt
        do j=1,jmt
           if (kmt_big(i,j)>0)  htp(i,j)=zw(kmt_big(i,j))
        enddo
       enddo

#ifdef notdef
c       a ridge 
        htp=0.; kmt_big=0
        do i=2,imt-1
         do j=2,jmt-1
          htp(i,j)=5500-3000.* exp(- (xt(i)-xt(imt/2+1))**2/15.**2 )
          htp(i,j)=htp(i,j)*100.
          kmt_big(i,j) = indp (htp(i,j), zw, km)
          if (zw( kmt_big(i,j) ) < htp(i,j) )
     &       kmt_big(i,j) = kmt_big(i,j)+1
         enddo
        enddo
#endif
c
#ifdef notdef
c       flat bottom
        kmt_big(2:imt-1,2:jmt-1)=km
        htp(2:imt-1,2:jmt-1) = zw(km)
#endif
      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
      implicit none
      integer :: i,j,k,io
      real    :: baru(is_pe:ie_pe)
      real    :: Bv0,T0,g,stoe,x,y,z,h_W
      integer :: ii,jj,kk,indp
      real :: delta = 0.0e-13 ! 1/s /cm
      real :: u_y   = 5.0e-6  ! 1/s
      real :: u_0   = 0.0     ! cm/s
      real :: u_z   = 1.0e-3  ! 1/s
      real :: amp   = 0.1     ! cm/s
      real :: tscl   = 1.0/(3.0*86400.0)   ! 1/s
      real :: Bvsqr(km)
      namelist /equat/ delta,u_y,u_0,u_z,Bv0,amp,h_W,tscl

c#ifdef notdef
c-----------------------------------------------------------------------
c     inertial instability
c-----------------------------------------------------------------------
      g=9.81e2 ! cm/s^2
      Bv0 = 0.002 ! 1/s
      h_W = 300.0e2

      call getunit(io,'namelist.equat','fsr' )
      read(io, nml = equat )
      close(io)

      Bvsqr(1)= Bv0**2*dzw(0)
      do k=2,km
       Bvsqr(k)= Bvsqr(k-1)+Bv0**2*exp(-zw(k-1)/h_W) *dzw(k-1)
      enddo

      do j=js_pe,je_pe
        do k=1,km
         do i=is_pe,ie_pe
           y=yu(j)* radius/180.0*pi  ! in centimeters
           u(i,k,j,1,:) = u_0 + u_y*y+ delta*y*y/2.0 
c     &       + u_z*zt(k)
     &       + amp*(sin(  zt(k)/120.0e2  *2*pi)
     &       +      sin(  zt(k)/220.0e2  *2*pi)
     &       +      sin(  zt(k)/320.0e2  *2*pi)
     &       +      sin(  zt(k)/420.0e2  *2*pi)
     &       +      sin(  zt(k)/520.0e2  *2*pi)
     &       +      sin(  zt(k)/620.0e2  *2*pi)  )
c           u(i,k,j,1,:) = min(u(i,k,j,1,:),20.0)
c           u(i,k,j,1,:) = max(u(i,k,j,1,:),-20.0)

           u(i,k,j,2,:) = 0.   

           y     =  yt(j)* radius/180.0*pi  ! in centimeters
           t(i,k,j,1,:) = 20.0+rho0/(g*alpha_lin)*Bvsqr(k)
           t(i,k,j,2,:) = 0.   ! salinity equals 35 PSU

           t(i,k,j,1,:) = t(i,k,j,1,:)*tmask(i,k,j)
           t(i,k,j,2,:) = t(i,k,j,2,:)*tmask(i,k,j)
           u(i,k,j,1,:) = u(i,k,j,1,:)*umask(i,k,j)
           u(i,k,j,2,:) = u(i,k,j,2,:)*umask(i,k,j)

         enddo 
        enddo
      enddo

c#endif

#ifdef notdef
c
c-----------------------------------------------------------------------
c       a constant background stratification (const N^2)
c-----------------------------------------------------------------------
c
         Bv0 = 0.0026; T0=20.0; g=9.81; 
c         Bv0 = 0.0015; T0=20.0; g=9.81; 
c         Bv0=0. ! no stratification

         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
            h_W = 500.
            x=(float(i)/float(imt)-0.5)*2.
            y=(float(j)/float(jmt)-0.5)*2.

c        no pertubation
            stoe = 0. 
c        a small but large scale gaussian pertubation
c            stoe = 0.5*exp(-(x-0.5)**2*10.)*exp(-(y)**2*10.)
c        a greater large scale gaussian pertubation
c            stoe = 20.*exp(-(x-0.5)**2*10.)*exp(-(y)**2*10.)
c        a meridional gradient with a small pertubation
c            stoe = 4.*y+0.2*exp(-(x-0.5)**2/0.01)*exp(-(y)**2/0.01)
            T0=20.0
#ifdef partial_cell
            T0 = T0+stoe*exp(-ztp(i,k,j)/100. /h_W )
            t(i,k,j,1,:) = T0+1/(g*alpha_lin)*ztp(i,k,j)/100.*Bv0**2
#else
            T0 = T0+stoe*exp(-zt(k)/100. /h_W )
            t(i,k,j,1,:) = T0+1/(g*alpha_lin)*zt(k)/100.*Bv0**2
#endif
            t(i,k,j,1,:) = t(i,k,j,1,:)*tmask(i,k,j)

c            u(i,k,j,1,:) = 30. ! cm/s
c            u(i,k,j,2,:) = 0.  ! cm/s

           enddo 
          enddo
         enddo
#endif
c
#ifdef notdef
c-----------------------------------------------------------------------
c        A mesoscale front
c-----------------------------------------------------------------------
c
         Bv0 = 0.0000; T0=20.0; g=9.81; 
         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=zt(km)/100.
            h_w = y/2.+tanh( (x-0.5)/0.3 )*y/5.
            T0=20 + 1.-tanh( (zt(k)/100.-h_w)/(y/5.)  )
            t(i,k,j,1,:) = T0+1/(g*alpha_lin)*zt(k)/100.*Bv0**2
            t(i,k,j,1,:) = t(i,k,j,1,:)*tmask(i,k,j)
           enddo 
          enddo
         enddo
#endif

#ifdef notdef
c
c-----------------------------------------------------------------------
c        Dammbruch exp.
c-----------------------------------------------------------------------
c
         ii = indp(330.5,xt,imt) 
         jj = indp(66.0,yt,jmt) 
         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
            if (j>=jj .and. i>=ii .and. k>=17) then
c            if (k>=17) then
              t(i,k,j,1,:) = 0.   
            else
              t(i,k,j,1,:) = 5.
            endif
           enddo
          enddo
         enddo
#endif


c-----------------------------------------------------------------------
c        barotropic mode
c-----------------------------------------------------------------------

         if (enable_freesurf) then
c
c         set the free surface
c
#ifdef notdef
          ubar=0.
          do j=js_pe,je_pe
           do i=is_pe,ie_pe
            y=(float(j)/float(jmt)-0.5)*2.
            etat(i,j,:)=0.
            if (j<jmt/2) etat(i,j,:)= 50.
            T0= model_dens_scalar(
     &           t(i,1,j,1,tau), t(i,1,j,2,tau),1
#ifdef partial_cell
     &                       ,ztp(i-1,1,j)
#endif
     &                       )
            ps(i,j,:) = etat(i,j,1)*grav*T0/rho0
           enddo
          enddo
#endif

         else ! enable_freesurf
c
c        diagnose streamfunction from full velocities
c        (works only for du/dx=0)
c
         psi = 0.0
         do j=js_pe,je_pe
c         depth integrated velocity
          baru=0.
          do k=1,km
           do i=is_pe,ie_pe
            baru(i) = baru(i) + u(i,k,j,1,0)*umask(i,k,j)
#ifdef partial_cell
     &                                    *dhu(i,k,j)
#else
     &                                    *dzt(k)
#endif
           enddo
          enddo
c         local integrant
          do i=is_pe,ie_pe
            psi(i,j+1,1) = -dyu(j)*baru(i)*umask(i,1,j)
          enddo
         enddo
         call border_exchg(psi(:,:,1),1,1)
         call set_cyclic(  psi(:,:,1),1,1)
c        Meridionally integrate psiu.
         do j=js_pe,je_pe
          do i=is_pe,ie_pe
            psi(i,j,1) = psi(i,j,1) + psi(i,j-1,1)
            psi(i,j,2) = psi(i,j,1) 
          enddo
         enddo
        endif ! enable_freesurf

      end subroutine init_cond_template



      subroutine sbc_template
c-----------------------------------------------------------------------
c     setup simple surface boundary conditions
c-----------------------------------------------------------------------
      use spflame_module
      implicit none
      integer :: i,j,k,is,ie,js,je
      real            :: y,x

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

c     no wind stress      
      smf=0.

#ifdef notdef
c     cosine shaped zonal wind stress
      do j=js,je
       do i=is,ie
          y = (yu(j)-yu(1))/ (yu(jmt)*1.0-yu(1)) ! cosine shaped wind stress
          smf(i,j,1) =-0.6*cos(2*pi*y)*umask(i,1,j) ! in taux only
          smf(i,j,2) =0.
       enddo
      enddo
#endif

#ifdef notdef
c       zonal wind stress in western part of domain, El Nino setup
	is=1;   is=max(is,is_pe)
	ie = nint(imt/6.) ; ie = min(ie,ie_pe)
 	if (is<=ie) smf(is:ie,js:je,1)=-0.5*umask(is:ie,1,js:je)
        is=max(is_pe,2); ie=min(ie_pe,imt-1)
#endif

c#ifdef notdef
c       no heat or salt fluxes
        stf_clim=0.
        stf_rest=0.
c#endif


#ifdef notdef
c       restoring to a profile
        do j=js,je
         do i=is,ie
c          y = (yt(j)-yt(1))/ (yt(jmt)*1.0-yt(1))
          y = ((j-1.)/(jmt-1.)-0.5)*2/1000.*2 ! this is a density -2:2, rho/alpha=  T
          stf_clim(i,j,1)=20.+y/alpha_lin
          stf_clim(i,j,2)=0.        ! 35 psu for flame density
          stf_rest(i,j,:)=1./(10.*86400.)*dzt(1)  ! 10 days
         enddo
        enddo
#endif
#ifdef notdef
c       restoring to a profile
        do j=js,je
         do i=is,ie
c        stf=stf_rest*(stf_clim-t)
          y = (yt(j)-yt(1))/ (yt(jmt)*1.0-yt(1))
          stf_clim(i,j,1)=t(i,1,j,1,taum1)+1
          stf_rest(i,j,1)=-500.*(y-0.5)/41868.
c          stf_clim(i,j,2)=t(i,1,j,2,taum1)+1.
c          stf_rest(i,j,2)=-5e-6
         enddo
        enddo
#endif

c       no wind stirring
        do j=js,je
         do i=is,ie
          if (enable_ktmix.or.enable_tkemix) ustar(i,j)=0.     ! no wind stirring
         enddo
        enddo

      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
      integer :: i,j,k,is,ie,js,je
      real :: Bv0,T0,g
      real :: temp(is_pe-1:ie_pe+1), salt(is_pe-1:ie_pe+1)

#ifdef notdef
      is=max(1,is_pe-1); ie = min(imt,ie_pe+1)
      js=max(1,js_pe-1); je = min(jmt,je_pe+1)
c
c-----------------------------------------------------------------------
c     a constant background stratification (const N^2)
c-----------------------------------------------------------------------
c
      blue_p1 = 1
      blue_p2 = 2
      blue_f1 = 0.5
      blue_f2 = 0.5

      Bv0 = 0.0026; T0=20.0; g=9.81; 
      do j=js,je
       do k=1,km
         do i=is,ie
            salt(i) = 0.   ! salinity equals 35 PSU
#ifdef partial_cell
            temp(i) = T0+1/(g*alpha_lin)*ztp(i,k,j)/100.*Bv0**2
#else
            temp(i) = T0+1/(g*alpha_lin)*zt(k)/100.*Bv0**2
#endif
            temp(i) = temp(i)*tmask(i,k,j)
         enddo 
         call model_dens(temp(is), salt(is),
     &                  blue_rho(is,k,j,blue_p1),k,ie-is+1
#ifdef partial_cell
     &                       ,ztp(is,k,j)
#endif
     &                       )
       enddo
      enddo
      blue_rho(:,:,:,blue_p2)=blue_rho(:,:,:,blue_p1)
#endif

      end subroutine blue_template



      subroutine sponge_template(n)
c-----------------------------------------------------------------------
c     setup simple sponge layers. 
c-----------------------------------------------------------------------
      use spflame_module
      implicit none
      integer, intent(in) :: n ! number of the tracer
      logical, save :: init = .true.

      integer :: i,j,k,io,is,ie,js,je
      real :: y,z,fx,tstar
      real,save :: g   = 9.81e2 ! cm/s^2
      real,save :: Bv0 = 0.002 ! 1/s
      real,save :: h_W = 300.0e2
      real,save :: delta = 0.0e-13 ! 1/s /cm
      real,save :: u_y   = 5.0e-6  ! 1/s
      real,save :: u_0   = 0.0     ! cm/s
      real,save :: u_z   = 1.0e-3  ! 1/s
      real,save :: amp   = 0.1     ! cm/s
      real,save,allocatable :: Bvsqr(:)
      real,save :: tscl   = 1.0/(3.0*86400.0)   ! 1/s

      namelist /equat/ delta,u_y,u_0,u_z,Bv0,amp,h_W,tscl

c-----------------------------------------------------------------------
c     inertial instability
c-----------------------------------------------------------------------

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      if (init) then

       call getunit(io,'namelist.equat','fsr' )
       read(io, nml = equat )
       close(io)
       allocate(Bvsqr(km))

       Bvsqr(1)= Bv0**2*dzw(0)
       do k=2,km
        Bvsqr(k)= Bvsqr(k-1)+Bv0**2*exp(-zw(k-1)/h_W) *dzw(k-1)
       enddo
       init = .false.
      endif

      if (n==1) then
       do j=js,je
        do k=1,km
         do i=is,ie
           y     =  yt(j)* radius/180.0*pi  ! in centimeters
           fx    =  (cori(i,j-1,1)+cori(i,j,1))/2.0
           tstar = 20.0+rho0/(g*alpha_lin)*(Bvsqr(k)-fx*u_z*y)
           source(i,k,j)=source(i,k,j)+
     &        tscl*(tstar-t(i,k,j,n,taum1))*tmask(i,k,j)
         enddo 
        enddo
       enddo
      endif

      end subroutine sponge_template
