#include "options.inc"

c
c-----------------------------------------------------------------------
c   Templates to read streamfunction for open boundaries
c   from external files.
c
c   Notes:
c      two subroutines are defined at the moment: 
c       "read_psi_south" and "read_psi_north"
c      Borth called by prep_obc (in prep_main.F)
c      psi should be in cm^3/s = 1e-12 Sv
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 read_psi_south(psi,kmt,imt)
       implicit none
       integer, intent(in) :: imt, kmt(imt)
       real                :: psi(12,imt)
       integer             :: n,i

#if defined read_psi_for_E1 ||defined read_psi_for_E7
       character (len=80) :: name
       character (len=32) :: stamp
       integer,parameter :: imtf=263,imtuf=348
       integer :: j,k
       real (kind=4) psiD3(imtuf),dpm
       name='data/PSI/D3_psi_south.mom'
       print*,' reading streamfct data from file ',name
       open(10,file=name,form='unformatted')
c      read data from FLAME D grid
       read(10) ! awi_id 
       do n=1,12
        read(10) ! iotext 
        read(10) stamp,dpm,i,j,k,psiD3(1:imtf)
c       unfold
        do i=imtf-1,imtuf
         psiD3(i)=psiD3(i-imtf+2)
        enddo
c       transfer to FLAME E grid
        do i=1,imt
         psi(n,i)=(psiD3(i*4)+psiD3(i*4-1)+psiD3(i*4-2)+psiD3(i*4-3))/4.
         if (kmt(i)==0) psi(n,i)=0.
        enddo
       enddo
       close(10)
#elif defined read_psi_for_D3 || defined read_psi_for_D4
       character (len=80) :: name
       character (len=32) :: stamp
       integer,parameter :: imtf=263,imtuf=348
       integer :: j,k
       real (kind=4) psiD3(imtuf),dpm
       name='data/PSI/D3_psi_south.mom'
       print*,' reading streamfct data from file ',name
       open(10,file=name,form='unformatted')
c      read data from FLAME D grid
       read(10) ! awi_id 
       do n=1,12
        read(10) ! iotext 
        read(10) stamp,dpm,i,j,k,psiD3(1:imtf)
c       unfold
        do i=imtf-1,imtuf
         psiD3(i)=psiD3(i-imtf+2)
        enddo
c       transfer
        do i=1,imt
         psi(n,i)=psiD3(i)
         if (kmt(i)==0) psi(n,i)=0.
        enddo
       enddo
       close(10)
#elif defined read_psi_for_HS2 || defined read_psi_for_HS3 || defined read_psi_for_HS4
       character (len=80) :: name
       character (len=32) :: stamp
       integer,parameter :: imtf=1021,imtuf=1394
       integer :: j,k
       real (kind=4) psiD3(imtuf),dpm
       name='data/PSI/H2_psi_south.mom'
       print*,' reading streamfct data from file ',name
       call sub_flush(6)
       open(10,file=name,form='unformatted')
c      read data from FLAME D grid
       read(10) ! awi_id 
       do n=1,12
        read(10) ! iotext 
        read(10) stamp,dpm,i,j,k,psiD3(1:imtf)
c       unfold
        do i=imtf-1,imtuf
         psiD3(i)=psiD3(i-imtf+2)
        enddo
c       transfer
        do i=1,imt
         psi(n,i)=psiD3(i)
         if (kmt(i)==0) psi(n,i)=0.
        enddo
       enddo
       close(10)
#elif defined read_psi_for_K1
       character (len=80) :: name
       name='data/PSI/K1_psi_south.dat'
       print*,' reading streamfct data from file ',name
       open(10,file=name,form='formatted')
       do n=1,12
        read(10,*) (psi(n,i),i=1,imt)
       enddo
       close(10)
c      convert to model units
       psi=psi*1e12
#elif defined read_psi_for_Q1
       character (len=80) :: name
       name='data/PSI/Q1_psi_south.dat'
       print*,' reading streamfct data from file ',name
       open(10,file=name,form='formatted')
       do n=1,12
        read(10,*) (psi(n,i),i=1,imt)
       enddo
       close(10)
c      convert to model units
       psi=psi*1e12
#endif
c
c      make sure that boundary cond. for psi  is met
c
       do n=1,12
        do i=2,imt-1
         if (kmt(i)>0 .and. kmt(i-1) ==0) psi(n,i)=0.
         if (kmt(i)>0 .and. kmt(i+1) ==0) psi(n,i)=0.
        enddo
       enddo
       end subroutine read_psi_south




       subroutine read_psi_north(psi,kmt,imt)
       implicit none
       integer, intent(in) :: imt, kmt(imt)
       real                :: psi(12,imt)
       integer             :: n,i

#ifdef read_psi_for_Q1
       character (len=80) :: name
       name='data/PSI/Q1_psi_north.dat'
       print*,' reading streamfct data from file ',name
       open(10,file=name,form='formatted')
       do n=1,12
        read(10,*) (psi(n,i),i=1,imt)
       enddo
       close(10)
c      convert to model units
       psi=psi*1e12
#elif defined read_psi_for_D4
       character (len=80) :: name
       real :: x
       name='data/PSI/D4_psi_north.dat'
       print*,' reading streamfct data from file ',name
       open(10,file=name,form='formatted')
       do n=1,12
        read(10,*) (psi(n,i),i=1,imt)
       enddo
       close(10)
c      care about transports near coast
       do i=232,244
        x= (i-232.0)/12.0
c        psi(:,i) = psi(:,i)*exp( - (1.0-x)*10.0 )
        psi(:,i) = psi(:,i)*x
       enddo
c      enhance transports artifically
       psi=psi*2.0

       do i=335,346
        x= (i-335.0)/11.0
        psi(:,i) = psi(:,i)*(1.0-x)
       enddo
c      convert to model units
       psi=psi*1e12
#elif defined read_psi_for_E7
       character (len=80) :: name
       real :: x
       name='data/PSI/E7_psi_north.dat'
       print*,' reading streamfct data from file ',name
       open(10,file=name,form='formatted')
       do n=1,12
        read(10,*) (psi(n,i),i=1,imt)
       enddo
       close(10)
c      care about transports near coast
       psi(:,59)= 0.0
       psi(:,60)=-2.0
c      enhance transports artifically
       psi=psi*2.0
c      convert to model units
       psi=psi*1e12
#elif defined read_psi_for_HS3
       character (len=80) :: name
#include "netcdf.inc"
       integer :: iret,ncid,vid
#ifdef netcdf_real4
      real (kind=4) :: sp,v4(imt)
#else
      real          :: sp,v4(imt)
#endif
      real :: x

       name='data/PSI/HS3_psi_north.cdf'
       print*,' reading streamfct data from file ',name
       call sub_flush(6)
       iret=nf_open(name,NF_nowrite,ncid)
       if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       iret=nf_inq_varid(ncid,'psi',vid)
       iret= nf_get_att_real (ncid,vid ,'missing_value',sp)
       do n=1,12
        iret= nf_get_vara_real(ncid,vid,(/1,9,n/),(/imt,1,1/),v4)
        where ( v4 == sp ) v4=0.0
        psi(n,:)=v4
       enddo
       call ncclos (ncid, iret)
c      care about transports near coast
       psi=min(0.0,psi) 
       do i=1350,1390
        x= (i-1350.0)/40.0
        psi(:,i) = psi(:,i)*(1.0-x)
       enddo
       psi(:,1390:imt)=0.0
c      convert to model units
       psi=psi*1e12

#elif defined read_psi_for_LS1
       character (len=80) :: name
       real :: x
       name='data/PSI/LS1_psi_north.dat'
       print*,' reading streamfct data from file ',name
       open(10,file=name,form='formatted')
       do n=1,12
        read(10,*) (psi(n,i),i=1,imt)
       enddo
       close(10)
c      care about transports near coast
       do i=582,627
        x= (i-582.0)/(627.0-582.0)
        psi(:,i) = psi(:,i)*x
       enddo
c      enhance transports artifically
       psi=psi*2.0

       do i=995,1040
        x= (i-995.0)/(1040.-995)
        psi(:,i) = psi(:,i)*(1.0-x)
       enddo
c      convert to model units
       psi=psi*1e12

#elif defined read_psi_for_HS4
       character (len=80) :: name
       real :: x
       name='data/PSI/HS4_psinorth_j2.dat'
       print*,' reading streamfct data from file ',name
       open(10,file=name,form='formatted')
       do n=1,12
        read(10,*) (psi(n,i),i=1,imt)
       enddo
       close(10)
c      care about transports near coast
       do i=931,976
        x= (i-931.0)/(976.0-931.0)
        psi(:,i) = psi(:,i)*x
       enddo
c      enhance transports artifically
       psi=psi*2.0

       do i=1344,1389
        x= (i-1344.0)/(1389.0-1344.0)
        psi(:,i) = psi(:,i)*(1.0-x)
       enddo
c      convert to model units
       psi=psi*1e12
#elif defined read_psi_for_A1
c     ACC model
       character (len=80) :: name
       real :: x
       name='data/PSI/nonmae.dat'
       print*,' reading streamfct data from file ',name
c       open(10,file=name,form='formatted')
c       do n=1,12
c        read(10,*) (psi(n,i),i=1,imt)
c       enddo
c       close(10)
c      care about transports near coast

c      convert to model units
       psi=0*psi*1e12

#endif
c
c      make sure that boundary cond. for psi  is met
c
       do n=1,12
        do i=2,imt-1
         if (kmt(i)>0 .and. kmt(i-1) ==0) psi(n,i)=0.
         if (kmt(i)>0 .and. kmt(i+1) ==0) psi(n,i)=0.
        enddo
       enddo
       end subroutine read_psi_north


