#include "options.inc"

c
c ------------------------------------------------------
c   Templates to create restoring zones (sponge layers) 
c   for tracers in the model domain.
c
c   Notes:
c    uses  "write_spg_to_cdf" (in prep_out.F) to
c    write out the sponge layer to netcdf output file.
c    temperature should be potential temperature in
c    deg C, salinity in (psu-35)/1000.
c    A single subroutine is defined "prep_sponge"
c    Header of subroutines should not be touched.
c    CPP directives are set in prep_template.F
c
c                            c.eden  nov 2002
c ------------------------------------------------------
c

       subroutine prep_sponge
       use prep_module
       implicit none
       real, allocatable :: time_scale(:,:,:,:)
       real, allocatable :: tr(:,:,:,:)
       real, allocatable :: trbuf(:,:,:,:)
       integer is,ie,js,je,js2,j1,i,j,k
       real :: dstnce
       integer :: n,indp,m
       real, allocatable :: buf(:,:,:,:)

       call init_forcing_file('forcing_sponge.cdf')

       allocate(  buf(imt,jmt,km,number_tr) )


#ifdef prep_sponge_for_K1

c sponge layer number 1 in the Greenland Sea
c use the same restoring timescales as in DYNAMO I
c (from file interpol.f written by P.Herrmann) 
c and also the same data !!
      js=indp(60.5,yt,jmt); je=jmt-1 
      is=indp(320.0,xt,imt); ie=imt
      allocate(time_scale(is:ie,js:je,km,number_tr) )
      time_scale=0. ! set inverse time scale to zero for passive tracers
      call peters_rest_time_scale(time_scale,is,ie,js,je)
      time_scale(:,:,:,2)=time_scale(:,:,:,1)
      allocate(tr(is:ie,js:je,km,number_tr) ); tr=0.
      do n=1,12
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)
c sponge layer 2 for the med sea
c search for point (i,j) of med sea outflow
      ie=indp(354.0,xt,imt)
      j1=indp(36.0,yt,jmt)
c search for southernmost and northernmost jrows for rest. zone
      js=indp(yt(j1)-7.,yt,jmt) 
      je=indp(yt(j1)+7.,yt,jmt)
      is=ie-3
      allocate(time_scale(is:ie,js:je,km,number_tr) ); time_scale=0.
      do i=is,ie; do j=js,je
        if (dstnce(yt(j),xt(i),yt(j1),xt(ie))<=600.0) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,1:2)=1.0/(10.0*86400.0)
     &             *exp(- (zt(k)/100.0-1200.0)**2/400.0**2 )
         enddo
        endif
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr) )
      do n=1,12
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)
#endif




#if defined prep_sponge_for_E1 || defined prep_sponge_for_E0

c sponge layer number 1 in the Greenland Sea
c use the same restoring timescales as in DYNAMO I
c (from file interpol.f written by P.Herrmann) 
c and also the same data !!
      print*,' setting up sponge 1 for E1/E0 setup'
      js=indp(60.5,yt,jmt); je=jmt-1 
      is=indp(320.0,xt,imt); ie=imt
      allocate(time_scale(is:ie,js:je,km,number_tr) )
      time_scale=0. ! set inverse time scale to zero for passive tracers
      call peters_rest_time_scale(time_scale,is,ie,js,je)
      time_scale(:,:,:,2)=time_scale(:,:,:,1)
      allocate(tr(is:ie,js:je,km,number_tr) )
      allocate(trbuf(imt,jmt,km,number_tr) )
      tr=0.;trbuf=0.
      do n=1,12
       call read_ts_from_file(n,trbuf)
       call read_schiller_sponges(trbuf(:,:,:,1),trbuf(:,:,:,2),n)
       tr=trbuf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale,trbuf)
      print*,'done'
c sponge layer 2 for the med sea
c search for point (i,j) of med sea outflow
      print*,' setting up sponge 2 for E1/E0 setup'
      ie=indp(354.0,xt,imt)
      j1=indp(36.0,yt,jmt)
c search for southernmost and northernmost jrows for rest. zone
      js=indp(yt(j1)-7.,yt,jmt) 
      je=indp(yt(j1)+7.,yt,jmt)
      is=ie-3
      allocate(time_scale(is:ie,js:je,km,number_tr) ); time_scale=0.
      do i=is,ie; do j=js,je
        if (dstnce(yt(j),xt(i),yt(j1),xt(ie))<=600.0) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,1:2)=1.0/(10.0*86400.0)
     &             *exp(- (zt(k)/100.0-1200.0)**2/400.0**2 )
         enddo
        endif
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr) )
      do n=1,12
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)
      print*,'done'
#endif

#if defined prep_sponge_for_E0
c
c sponge layer 3 for southern boundary
c
      print*,' setting up sponge 3 for E0 setup'
      is=2; ie=imt-1; js=2; je=5
      allocate(time_scale(is:ie,js:je,km,number_tr) ); 
      time_scale=0.
      do i=is,ie; do j=js,je
        do k=1,kmt(i,j)
          time_scale(i,j,k,1:2)=1.0/(3.0*86400.0)
        enddo
#ifdef  prep_for_NPZD
        if (number_tr >=3) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,3)=1.0/(3.0*86400.0)
         enddo
        endif
        if (number_tr >=7) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,7:numbert_tr)=1.0/(3.0*86400.0)
         enddo
        endif
c        if (number_tr >=8) then
c         do k=1,kmt(i,j)
c          time_scale(i,j,k,8)=1.0/(3.0*86400.0)
c         enddo
c        endif
#endif
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr)  )
      do n=1,12
       print*,' month=',n
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)
      print*,'done'
#endif

#ifdef prep_sponge_for_E7

c sponge layer 1 for the med sea
      print*,'generating sponge for E7: Med outflow'
c search for point (i,j) of med sea outflow
      ie=indp(354.0,xt,imt)
      j1=indp(36.0,yt,jmt)
c search for southernmost and northernmost jrows for rest. zone
      js=indp(yt(j1)-7.,yt,jmt) 
      je=indp(yt(j1)+7.,yt,jmt)
      is=ie-3
      allocate(time_scale(is:ie,js:je,km,number_tr) ); time_scale=0.
      do i=is,ie; do j=js,je
        if (dstnce(yt(j),xt(i),yt(j1),xt(ie))<=600.0) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,1:2)=1.0/(10.0*86400.0)
     &             *exp(- (zt(k)/100.0-1200.0)**2/400.0**2 )
         enddo
        endif
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr) )
      do n=1,12
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)
#endif



#ifdef prep_sponge_for_D3

c sponge layer number 1 in the Greenland Sea
c use the same restoring timescales as in DYNAMO I
c (from file interpol.f written by P.Herrmann) 
      js=indp(60.5,yt,jmt); je=jmt-1 
      is=indp(320.0,xt,imt); ie=imt
      allocate(tr(is:ie,js:je,km,1) ); tr=0.
      call peters_rest_time_scale(tr,is,ie,js,je)
      js2=jmt-10
      allocate(time_scale(is:ie,js2:je,km,number_tr) )
      time_scale=0.
      time_scale(is:ie,js2:je,1:km,1) = tr(is:ie,js2:je,1:km,1)
      time_scale(is:ie,js2:je,1:km,2) = tr(is:ie,js2:je,1:km,1)
      deallocate( tr) 
      js=js2
      allocate(tr(is:ie,js:je,km,number_tr) )
      allocate(trbuf(imt,jmt,km,number_tr)  )
      tr=0.;trbuf=0.
      do n=1,12
       call read_ts_from_file(n,trbuf)
       call read_schiller_sponges(trbuf(1,1,1,1),trbuf(1,1,1,2),n)
       tr=trbuf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale,trbuf)
c sponge layer 2 is for the med sea
c search for point (i,j) of med sea outflow
      ie=282
      j1=172
c search for southernmost and northernmost jrows for rest. zone
      js=152
      je=192
      is=263 ! = imtf
      allocate(time_scale(is:ie,js:je,km,number_tr) ); time_scale=0.
      do i=is,ie; do j=js,je
        if (dstnce(yt(j),xt(i),yt(j1),xt(ie))<=600.0) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,1:2)=1.0/(10.0*86400.0)
     &             *exp(- (zt(k)/100.0-1200.0)**2/400.0**2 )
         enddo
        endif
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr)  )
      do n=1,12
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)
#endif


#ifdef prep_sponge_for_D4

      print*,' setting up Med Sea sponge for D4 setup'
c sponge layer 1 is for the med sea
c search for point (i,j) of med sea outflow
      ie=282
      j1=172
c search for southernmost and northernmost jrows for rest. zone
      js=152
      je=192
      is=263 ! = imtf
      allocate(time_scale(is:ie,js:je,km,number_tr) ); time_scale=0.
      do i=is,ie; do j=js,je
        if (dstnce(yt(j),xt(i),yt(j1),xt(ie))<=600.0) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,1:2)=1.0/(10.0*86400.0)
     &             *exp(- (zt(k)/100.0-1200.0)**2/400.0**2 )
         enddo
        endif
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr)  )
      do n=1,12
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)

#ifdef  prep_for_NPZD
      print*,' setting up northern NO3 sponge for D4 setup'
      is=2; ie=imt-1; js=jmt-1-10; je=jmt-1
      allocate(time_scale(is:ie,js:je,km,number_tr) ); 
      time_scale=0.
      do i=is,ie; do j=js,je
        if (number_tr >=3) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,3)=1.0/(3.0*86400.0)
         enddo
        endif
        if (number_tr >=7) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,7)=1.0/(3.0*86400.0)
         enddo
        endif
        if (number_tr >=8) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,8)=1.0/(3.0*86400.0)
         enddo
        endif
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr)  )
      do n=1,12
       print*,' month=',n
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)


      print*,' setting up southern NO3 sponge for D4 setup'
      is=2; ie=imt-1; js=2; je=5*2
      allocate(time_scale(is:ie,js:je,km,number_tr) ); 
      time_scale=0.
      do i=is,ie; do j=js,je
        if (number_tr >=3) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,3)=1.0/(3.0*86400.0)
         enddo
        endif
        if (number_tr >=7) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,7)=1.0/(3.0*86400.0)
         enddo
        endif
        if (number_tr >=8) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,8)=1.0/(3.0*86400.0)
         enddo
        endif
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr)  )
      do n=1,12
       print*,' month=',n
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)
#endif
      print*,'done'
#endif



#ifdef prep_sponge_for_HS2

c sponge layer number 1 in the Greenland Sea
c use the same restoring timescales as in DYNAMO I
c (from file interpol.f written by P.Herrmann) 
      js=indp(60.5,yt,jmt); je=jmt-1 
      is=indp(320.0,xt,imt); ie=imt
      allocate(tr(is:ie,js:je,km,1) ); tr=0.
      call peters_rest_time_scale(tr,is,ie,js,je)
      js2=jmt-37
      allocate(time_scale(is:ie,js2:je,km,number_tr) )
      time_scale=0.
      time_scale(is:ie,js2:je,1:km,1) = tr(is:ie,js2:je,1:km,1)
      time_scale(is:ie,js2:je,1:km,2) = tr(is:ie,js2:je,1:km,1)
      deallocate( tr) 
      js=js2
      allocate(tr(is:ie,js:je,km,number_tr) )
      allocate(trbuf(imt,jmt,km,number_tr)  )
      tr=0.;trbuf=0.
      do n=1,12
       call read_ts_from_file(n,trbuf)
       call read_schiller_sponges(trbuf(1,1,1,1),trbuf(1,1,1,2),n)
       tr=trbuf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale,trbuf)

c sponge layer 2 is for the med sea
c search for point (i,j) of med sea outflow
      ie=1127
      j1=indp(36.0,yt,jmt)
c search for southernmost and northernmost jrows for rest. zone
      js=603
      je=763
      is=1051 ! = imtf
      allocate(time_scale(is:ie,js:je,km,number_tr) ); time_scale=0.
      do i=is,ie; do j=js,je
        if (dstnce(yt(j),xt(i),yt(j1),xt(ie))<=600.0) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,1:2)=1.0/(10.0*86400.0)
     &             *exp(- (zt(k)/100.0-1200.0)**2/400.0**2 )
         enddo
        endif
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr)  )
      do n=1,12
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)
#endif

#ifdef prep_sponge_for_HS3

      print*,' setting up sponge 1 for HS3 setup'
      call sub_flush(6)

c the only sponge layer is in the med sea

      is=1280;ie=1315; js=695; je=817 
      allocate(time_scale(is:ie,js:je,km,number_tr) ); 
      time_scale=0.
      do i=is,ie; do j=js,je
        do k=1,kmt(i,j)
          time_scale(i,j,k,1:2)=1.0/(3.0*86400.0)
        enddo
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr)  )
      do n=1,12
       print*,' month=',n
       call sub_flush(6)
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)

#ifdef  prep_for_NPZD
      print*,' setting up northern NO3 sponge for HS3 setup'
      is=2; ie=imt-1; js=jmt-1-10; je=jmt-1
      allocate(time_scale(is:ie,js:je,km,number_tr) ); 
      time_scale=0.
      do i=is,ie; do j=js,je
        if (number_tr >=3) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,3)=1.0/(3.0*86400.0)
         enddo
        endif
        if (number_tr >=7) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,7)=1.0/(3.0*86400.0)
         enddo
        endif
        if (number_tr >=8) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,8)=1.0/(3.0*86400.0)
         enddo
        endif
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr)  )
      do n=1,12
       print*,' month=',n
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)


      print*,' setting up southern NO3 sponge for HS3 setup'
      is=2; ie=imt-1; js=2; je=5*2
      allocate(time_scale(is:ie,js:je,km,number_tr) ); 
      time_scale=0.
      do i=is,ie; do j=js,je
        if (number_tr >=3) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,3)=1.0/(3.0*86400.0)
         enddo
        endif
        if (number_tr >=7) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,7)=1.0/(3.0*86400.0)
         enddo
        endif
        if (number_tr >=8) then
         do k=1,kmt(i,j)
          time_scale(i,j,k,8)=1.0/(3.0*86400.0)
         enddo
        endif
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr)  )
      do n=1,12
       print*,' month=',n
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)
#endif
      print*,'done'

#endif

#ifdef prep_sponge_for_PAC

c sponge layer 1 is for Indonesian throughflow

      print*,' setting up sponge 1 for PACific setup'
      is=2; ie=4; js=45; je=66
      allocate(time_scale(is:ie,js:je,km,number_tr) ); 
      time_scale=0.
      do i=is,ie; do j=js,je
        do k=1,kmt(i,j)
          time_scale(i,j,k,1:2)=1.0/(3.0*86400.0)
        enddo
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr)  )
      do n=1,12
       print*,' month=',n
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)
      print*,'done'

      print*,' setting up sponge 2 for PACific setup'
      is=5; ie=15; js=48; je=50 !/i=5:15/j=48:50
      allocate(time_scale(is:ie,js:je,km,number_tr) ); 
      time_scale=0.
      do i=is,ie; do j=js,je
        do k=1,kmt(i,j)
          time_scale(i,j,k,1:2)=1.0/(3.0*86400.0)
        enddo
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr)  )
      do n=1,12
       print*,' month=',n
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)
      print*,'done'
#endif

#ifdef prep_sponge_for_I1

c sponge layer 1 is for Red Sea

      print*,' setting up sponge 1 for Indic setup'
      is=22; ie=35; js=143; je=154
      allocate(time_scale(is:ie,js:je,km,number_tr) ); 
      time_scale=0.
      do i=is,ie; do j=js,je
        do k=1,kmt(i,j)
          time_scale(i,j,k,1:2)=1.0/(3.0*86400.0)
        enddo
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr)  )
      do n=1,12
       print*,' month=',n
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)
      print*,'done'

      print*,' setting up sponge 2 for Indic setup'
      is=75; ie=93; js=202; je=206 !/i=5:15/j=48:50
      allocate(time_scale(is:ie,js:je,km,number_tr) ); 
      time_scale=0.
      do i=is,ie; do j=js,je
        do k=1,kmt(i,j)
          time_scale(i,j,k,1:2)=1.0/(3.0*86400.0)
        enddo
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr)  )
      do n=1,12
       print*,' month=',n
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)
      print*,'done'



      print*,' setting up sponge 3 for Indic setup'
      is=306; ie=337; js=71; je=77 
      allocate(time_scale(is:ie,js:je,km,number_tr) ); 
      time_scale=0.
      do i=is,ie; do j=js,je
        do k=1,kmt(i,j)
          time_scale(i,j,k,1:2)=1.0/(3.0*86400.0)
        enddo
      enddo; enddo
      allocate(tr(is:ie,js:je,km,number_tr)  )
      do n=1,12
       print*,' month=',n
       call read_ts_from_file(n,buf)
       tr(:,:,:,:)=buf(is:ie,js:je,:,:)
       do m=1,number_tr
        call write_spg_to_cdf(time_scale(:,:,:,m),
     &         tr(:,:,:,m),is,ie,js,je,12,n,m)
       enddo
      enddo
      deallocate(tr,time_scale)
      print*,'done'

#endif

      deallocate(buf)
      end subroutine prep_sponge

c
c    local functions and subroutines used by prep_sponge
c

      function dstnce(phi1,rla1,phi2,rla2)
c returns distance (in km) between two geographical
c points
c        written by   p.herrmann 
c      phi1 : first latitude 
c      rla1 : first longitude
c      phi2 : second latitude
c      rla2 : second longitude

      real p1,p2,rl1,rl2,x,xx
      data s/0.0174533/
      p1=phi1*s; p2=phi2*s; rl1=rla1*s; rl2=rla2*s
      x=sin(p1)*sin(p2)+cos(p1)*cos(p2)*cos(rl2-rl1)
      xx=abs(x)
      if(xx.gt.1.0) x=1.0
      d=atan(sqrt((1-x)/(1+x)))*222.24/s; dstnce=d
      end function dstnce


#if defined prep_sponge_for_E1 || defined prep_sponge_for_E0
      subroutine read_schiller_sponges(temp,salt,nl)
c      read the old schiller dataset from binary file
c      it is assumed that data are already interpolated
c      and folded. read data into temp and salt from prep_module
      use prep_module
      implicit none
      integer :: nl
      real :: temp(imt,jmt,km), salt(imt,jmt,km)
      integer (kind=4) :: io=10,jsx,jnx,m,i,k
c      integer :: imtf =67 ! imtf of FLAME E
c      integer :: imtuf=67 ! imtuf of FLAME E
c      integer :: jmtflame=89 ! imtuf of FLAME E
      integer :: imtf =79 ! imtf of FLAME B
      integer :: imtuf=98 ! imtuf of FLAME B
      integer :: jmtflame=150 ! imtuf of FLAME B
      character (len=32) s
      real (kind=4) dpm
      real (kind=4), allocatable :: bufn(:,:,:,:)
      integer j,jj
      character (len=80) :: infile = 'data/schiller/grid_B_folded.mom'
      print*,' reading TS data from file ',
     &       infile(1:len_trim(infile)),' for month=',nl
      open (io,file=infile,form='unformatted',status='old')
      read (io); read (io) jsx,jnx 
      allocate( bufn(imtuf,km,jnx,2) )
      do m=1,nl-1
       read (io) ; read (io)   ! data
      enddo
      read (io)   ! iotext
c      does not work for pflame option sponge1_converted
      read (io) s,dpm,i,i,i,i,i,i,
     &             (dpm,i=1,jmtflame),(dpm,i=1,jmtflame),
     &             (dpm,i=1,imtf*km*jsx*2),bufn(1:imtf,:,:,:)
c     unfold and transfer to temp/salt fields
      do i=imtf-1,imtuf
         bufn(i,:,:,:)= bufn(i-imtf+2,:,:,:)
      enddo
      do j=1,jnx
       jj=j+jmt-jnx-1
       do i=1,min(imt,imtuf)
        temp(i,jj,:)= bufn(i,:,j,1)
        salt(i,jj,:)= bufn(i,:,j,2)
       enddo
      enddo
      deallocate( bufn)
      close(io)
      do k=1,km
       where( k>kmt) temp(:,:,k)=spval
       where( k>kmt) salt(:,:,k)=spval
      enddo
      end subroutine  read_schiller_sponges
#endif

#if defined prep_sponge_for_D3 || defined prep_sponge_for_HS2

      subroutine read_schiller_sponges(temp,salt,nl)
c      read the old schiller dataset from binary file
c      it is assumed that data are already interpolated
c      and folded. read data into temp and salt from prep_module
      use prep_module
      implicit none
      integer :: nl
      real :: temp(imt,jmt,km), salt(imt,jmt,km)
      integer (kind=4) :: io=10,jsx,jnx,m,i,k
#ifdef prep_sponge_for_HS2
      integer :: imtf =1021 ! imtf of FLAME H
      integer :: imtuf=1394 ! imtuf of FLAME H
      integer :: jmtflame=1416 ! imtuf of FLAME H
      character (len=80) :: infile = 'data/schiller/grid_H_folded.mom'
#else
      integer :: imtf =263 ! imtf of FLAME D
      integer :: imtuf=348 ! imtuf of FLAME D
      integer :: jmtflame=356 ! imtuf of FLAME D
      character (len=80) :: infile = 'data/schiller/grid_D_folded.mom'
#endif
      character (len=32) s
      character(len=80) iotext
      real (kind=4) dpm
      real (kind=4), allocatable :: bufn(:,:,:,:)
      integer j,jj
      print*,' reading TS data from file ',
     &       infile(1:len_trim(infile)),' for month=',nl
      open (io,file=infile,form='unformatted',status='old')
      read (io) ! FLAME_ID
      read (io) jsx,jnx 
      allocate( bufn(imtuf,km,jnx,2) )
      do m=1,nl-1
       read (io) iotext
       read (io)   ! stamp et al
       do j=1,jsx
        read (io)   ! souther sponge  temp
        read (io)   ! souther sponge  salt
       enddo
       do j=1,jnx
        read (io)   ! northern sponge  temp
        read (io)   ! northern sponge  salt
       enddo
      enddo
      read (io)   iotext
c      print*,'iot-',iotext
      read (io)   ! stamp et al
      do j=1,jsx
        read (io)   ! souther sponge  temp
        read (io)   ! souther sponge  salt
      enddo
      do j=1,jnx
        read (io) ((bufn(i,k,j,1),i=1,imtf),k=1,km)
        read (io) ((bufn(i,k,j,2),i=1,imtf),k=1,km)
      enddo
c     unfold and transfer to temp/salt fields
      do i=imtf-1,imtuf
         bufn(i,:,:,:)= bufn(i-imtf+2,:,:,:)
      enddo
      do j=1,jnx
       jj=j+jmt-jnx-1
       do i=1,min(imt,imtuf)
        temp(i,jj,:)= bufn(i,:,j,1)
        salt(i,jj,:)= bufn(i,:,j,2)
       enddo
      enddo
      deallocate( bufn)
      do k=1,km
       where( k>kmt) temp(:,:,k)=spval
       where( k>kmt) salt(:,:,k)=spval
      enddo
      close(io)
      end subroutine  read_schiller_sponges
#endif




      subroutine peters_rest_time_scale(rest_tn,i1,ie,jsn,jen)
c    creates the time scale for the northern restoring
c    zone in the DYNAMO models.
      use prep_module
      implicit none
      integer i1,ie,jsn,jen
      real rest_tn(i1:ie,jsn:jen,km)
      integer i,j,k,j1,indp,i2
      real dx1,phi67,factor,south,dstnce
      i2=indp(350.0,xt,imt)
      j1=indp(67.0,yt,jmt)
      dx1=yt(jen)-yt(j1)
      do i=i1,i2
       do j=j1,jen
        do k=1,kmt(i,j)
         rest_tn(i,j,k)=
     &        1.0/((97.0/dx1*(yt(jen)-yt(j))+3.0)*86400.0)
        enddo
       enddo
      enddo
      phi67=yt(j1)
      factor=(yt(jsn)-yt(j1))/(xt(imt)-xt(i2))
      do i=i2+1,imt
       south=phi67+factor*(xt(i)-xt(i2))
       j1=indp(south,yt,jmt)
       dx1=yt(jen)-yt(j1)
       do j=j1,jen
        do k=1,kmt(i,j)
         rest_tn(i,j,k)=
     &    1.0/((97.0/dx1*(yt(jen)-yt(j))+3.0)*86400.0)
        enddo
       enddo
      enddo
      end subroutine peters_rest_time_scale
