#include "options.inc"
c
c-----------------------------------------------------------------------
c     preparation module
c-----------------------------------------------------------------------
c
      module prep_module
      implicit none
!     domain size
      integer imt,jmt,km
      integer :: number_tr = 2  ! numer of tracer

!     parameter which determines the grid, see MOM2 grids.F for explanation
      integer :: nxlon, nylat,nzdepth           
      integer,parameter :: max_regions = 100
      real :: x_lon(max_regions), dx_lon(max_regions)
      real :: y_lat(max_regions), dy_lat(max_regions)
      real :: z_depth(max_regions), dz_depth(max_regions)

!     the grid
      real, allocatable :: xt(:),xu(:)
      real, allocatable :: yt(:),yu(:)
      real, allocatable :: zt(:),zw(:),dzt(:),dzw(:)
      real, allocatable :: dxtdeg(:),dxudeg(:)
      real, allocatable :: dytdeg(:),dyudeg(:)
!     topography
      integer, allocatable :: kmt(:,:),kmu(:,:)
      real, allocatable :: htp(:,:),hup(:,:)  ! for partial cells

      logical :: cyclic= .false.   ! cyclic boundary condition
      logical :: enable_obc_north = .false.  ! open boundary at the north?
      logical :: enable_obc_south = .false.
      logical :: enable_obc_east = .false.
      logical :: enable_obc_west = .false.
      logical :: read_grid_topo = .false.    ! read in grid.dta and kmt.dta
      logical :: stop_after_topo = .false.   ! stop after generating 
                                             ! the grid and the
                                             ! topography
      logical :: isotropic = .false.         ! make an mercator projection
      logical :: flame_vert_grid = .false.   ! use the flame vertical grid
      logical :: pop_vert_grid   = .false.   ! use the POP 1/10 deg vertical grid


      logical :: rotated_grid = .true.       ! use other interpolation routines
      real (kind=8) phir,thetar,psir         ! rotation angles
      real (kind=8) gltnp,glnnp,gltpp,glnpp  ! latitude/long. of the 
                                             ! rotated north pole

      integer :: topo_smooth_ntimes=0, kmt_min=4 ! smoothing of topography,
                                                 ! minimal depth

      logical :: enable_blue      = .false.  ! prepare for blue
      logical :: enable_blue_mean = .false.  ! prepare for blue second approach

      real, parameter :: spval = -1.e20

      character (len=80) :: old_grid_file = 'old_grid.dta'
      character (len=80) :: new_grid_file = 'new_grid.dta'
      character (len=80) :: old_kmt_file  = 'old_kmt.dta'
      character (len=80) :: new_kmt_file  = 'new_kmt.dta'
      character (len=80) :: new_htp_file  = 'new_htp.dta'
      character (len=80) :: old_htp_file  = 'old_htp.dta'

      real :: tol_gcell = 1.0e-3  ! a factor used for the grid generation
      logical :: centered_t = .false. ! the grid can be t centered

      logical :: apply_landmask = .true.

      end module prep_module



      program prep_main
      use prep_module
      use time_type_module
      implicit none

      namelist /prep/ nxlon, nylat,nzdepth,
     &          x_lon,dx_lon,y_lat,dy_lat,z_depth,dz_depth,
     &          isotropic,flame_vert_grid,pop_vert_grid,cyclic,
     &          gltnp,glnnp,gltpp,glnpp,
     &          topo_smooth_ntimes,kmt_min,
     &          enable_obc_north, enable_obc_south,
     &          enable_obc_east,  enable_obc_west,
     &          read_grid_topo,stop_after_topo,
     &          rotated_grid, enable_blue, enable_blue_mean,
     &          number_tr,tol_gcell,centered_t

      call set_calendar_type(no_leap)

      open(10,file='namelist.prep',form='formatted')
      read(10, nml = prep )
      close(10)

      if (.not. rotated_grid) then
c       insure that no rotation at all is made
c       so that we can use simpler interpolation routines
        gltnp=90.;glnnp=90.
        gltpp=0.; glnpp=0.
      endif

      call prep_grid
      call sub_flush(6)
      call init_out_cdf
      call sub_flush(6)
      call prep_rot
      call sub_flush(6)
      call prep_topo
      call sub_flush(6)

      if (.not. stop_after_topo) then
       call prep_ic
       call sub_flush(6)
       call prep_sbc_tau
       call sub_flush(6)
       call prep_sbc_sflx
       call sub_flush(6)
       call prep_sponge
       call sub_flush(6)
       call prep_obc
       call sub_flush(6)
#ifdef BB_model
       call prep_BB()
       call sub_flush(6)
#endif
      endif

      end program prep_main



      subroutine prep_rot
      use prep_module
      implicit none
      integer i,j
      real (kind=8) rad
      real  (kind=8) rlnpp,rltpp,rltnp,rlnnp
      real  (kind=8) glt,gln,rlt,rln
      real,  allocatable :: glat(:,:)
      real :: p

c-----------------------------------------------------------------------
c     any spherical grid rotation can be achieved by specifying three, 
c     Euler angle, solid body rotations. the first rotation angle phir
c     is defined as a rotation about the original z axis. the second 
c     rotation angle thetar is defined as a rotation about the new x 
c     axis. the final rotation angle psir is defined as a rotation about
c     the new z axis. one method of deciding on the values of these 
c     angles is to define the geographic (unrotated) latitude and 
c     longitude of the north pole of the rotated grid. this defines 
c     the angles phir and thetar. The latitude of the new grid is now
c     specified but the longitude is still arbitrary. to define the new 
c     grid longitude, specify the geographic location of any point on 
c     the rotated grid`s prime meridian. the angle psir is then 
c     calculated as the rotation angle necessary to place this point on 
c     the prime meridian of the rotated grid.
c
c     for example, to set up a model with an equatorial grid over the 
c     Arctic and North Atlantic, one could specify the rotated grid pole
c     to be at 0 N, 110 W and a prime meridian point to be at 0 N O E. 
c     This specifies a grid rotation in which the new grid equator is 
c     defined along the 20 W and 160 E meridians. The rotated grid 
c     longitude would be east, north of the geographic equator and west 
c     to the south. The geographic pole would be at 0 N and 90 E.
c
c     latitudes and longitudes should be specified such that N and E are
c     positive and S and W are negative.
c 
c     gltnp = geographic latitude of the rotated grid north pole 
c     glnnp = geographic longitude of the rotated grid north pole
c     gltpp = geographic latitude of a prime point on the rotated grid
c     glnpp = geographic longitude of a prime point on the rotated grid
c-----------------------------------------------------------------------
c
c      if (precision(p) < 10) then
c       write(6,'(/a/)') '=>Warning: 64bit precision is recommended'
c       call halt_stop('')
c      endif
c
      rad = acos(-1.)/180.
c
c-----------------------------------------------------------------------
c     calculate the Euler angles required for the rotation
c-----------------------------------------------------------------------
c
      phir   = (glnnp - 90.)*rad
      thetar = (gltnp - 90.)*rad
c
c     set rotation angle psir to zero and rotate prime meridian point to
c     find the required angle psir to give the point a longitude of zero
c
      psir = 0.
      call rotate (gltpp, glnpp, phir, thetar, psir, rltpp, rlnpp)
      psir = rlnpp*rad

      print*,'psir=',psir/rad,' thetar=',thetar/rad,' phir=',phir/rad
c
c-----------------------------------------------------------------------
c     write out the grid definitions and rotation angles
c-----------------------------------------------------------------------
c
      call rotate (gltnp, glnnp, phir, thetar, psir, rltnp, rlnnp)
      call rotate (gltpp, glnpp, phir, thetar, psir, rltpp, rlnpp)
c
      write(6     ,'(/,a,/)') 'Rotated Model Grid Definitions:'
      write(6     ,'(16x,a)') '    Geographic              Rotated     '
      write(6     ,'(16x,a)')'Latitude  Longitude   Latitude  Longitude'
      write(6     ,'(a,4f11.3)') 'North Pole: ', gltnp, glnnp
     &,                           rltnp, rlnnp
      write(6     ,'(a,4f11.3)') 'Prime point:', gltpp, glnpp
     &,                           rltpp, rlnpp
c
      write(6     ,'(//,a,/)') 'Rotation Angles:'
      write(6     ,'(a,f11.3,a)') 'phir:       ', phir/rad, ' deg'
      write(6     ,'(a,f11.3,a)') 'thetar:     ', thetar/rad, ' deg'
      write(6     ,'(a,f11.3,a)') 'psir:       ', psir/rad, ' deg'
c
c-----------------------------------------------------------------------
c     calculate and write the geographic latitude of the model grid 
c-----------------------------------------------------------------------
c
c     latitude of velocity grid is needed to calculate Coriolis
c
      allocate( glat(imt,jmt) )
      do j=1,jmt
        do i=1,imt
          rlt=yu(j);rln=xu(i)
          call rotate (rlt, rln, -psir, -thetar, -phir, glt, gln)
          glat(i,j) = glt
        enddo
      enddo
      call write_glat_to_cdf(glat)
      deallocate(glat)
c
c-----------------------------------------------------------------------
c     write the locations of the corners of the rotated model grid
c-----------------------------------------------------------------------
c
      write(6     ,'(/,/,a,/)') 'Rotated Model Grid Corners:'
      write(6     ,'(16x,a)') '    Geographic              Rotated  '
      write(6     ,'(16x,a)') 
     &                'Latitude  Longitude   Latitude  Longitude'
      rlt=yt(jmt);rln=xt(1)
      call rotate (rlt, rln, -psir, -thetar, -phir, glt, gln)
      write(6     ,'(a,4f11.3)') 'NW corner:  ', glt, gln, yt(jmt)
     &,                             xt(1)
      rlt=yt(jmt);rln=xt(imt)
      call rotate (rlt, rln, -psir, -thetar, -phir, glt, gln)
      write(6     ,'(a,4f11.3)') 'NE corner:  ', glt, gln, yt(jmt)
     &,                             xt(imt)
      rlt=yt(1);rln=xt(1)
      call rotate (rlt, rln, -psir, -thetar, -phir, glt, gln)
      write(6     ,'(a,4f11.3)') 'SW corner:  ', glt, gln, yt(1), xt(1)
      rlt=yt(1);rln=xt(imt)
      call rotate (rlt, rln, -psir, -thetar, -phir, glt, gln)
      write(6     ,'(a,4f11.3,//)') 'SE corner:  ', glt, gln, yt(1)
     &,                             xt(imt)
      end subroutine prep_rot
      




       subroutine prep_ic
c
c--------------------------------------------------------
c      prepare the initial conditions
c      the data are also needed for sponge layers,
c      open boundaries and maybe for surface boundary conditions
c--------------------------------------------------------
c
       use prep_module
       implicit none
       integer nx,ny,nz,nt,n,k,j,m,ierr,i
       real, allocatable :: dtr(:,:,:,:),x(:),y(:),z(:),tx(:)
       real, allocatable :: ntr(:,:,:)
       real, allocatable :: tr(:,:,:,:)
       integer :: init_v,lenw_v
       real, allocatable :: work_v(:)
       integer :: init_h,lenw_h
       real, allocatable :: work_h(:)
       logical :: verbose=.true.

       call init_forcing_file('forcing_ic.cdf')
       if (number_tr>2) then
        call init_forcing_file('forcing_ic_tracer.cdf')
       endif

       print*,' reading global T/S fields from file'
       call read_dims_ic(nx,ny,nz,nt)
       print*,' nx=',nx,' ny=',ny,' nz=',nz,' nt=',nt
c      allocate working arrays
       allocate(dtr(nx,ny,nz,number_tr),x(nx),y(ny),z(nz+1),tx(nx)
     &       ,stat=ierr)
       if (ierr/=0) call halt_stop('allocate failed ')
       allocate(ntr(imt,jmt,nz+1))
       dtr=0.; x=0.; y=0.; z=0.; ntr=0.
       allocate(tr(imt,jmt,km,number_tr)); tr=0.

       print*,' interpolating data'
       do n=1,nt
c       read monthly means
        print*,' reading raw data for month =',n
        call sub_flush(6)
        call read_ic(nx,ny,nz,number_tr,n,x,y,z,dtr,spval)
        z(nz+1)=10000.
        if (.not.rotated_grid) then
         lenw_h=10*max(imt,jmt); allocate(work_h(lenw_h))
         init_h=1
        endif
        lenw_v=10*max(imt,km); allocate(work_v(lenw_v))
        init_v=1

        do m=1,number_tr
         print*,' working on tracer #',m
         print*,' interpolating horizontally '
         call sub_flush(6)
         do k=1,nz
          if (rotated_grid) then
           call rot_intrp_sclr(dtr(:,:,k,m), x, y, nx, ny, 
     &                        ntr(:,:,k), xt, yt,
     &               imt, jmt, psir, thetar, phir,spval,verbose)
           verbose=.false.
          else
c         use cheaper subroutine here
           call tranlon (dtr(:,:,k,m), nx, nx-2, ny, x, xt,imt, tx)
           call ctf (dtr(:,:,k,m), nx, ny, tx, y, 
     &                  ntr(:,:,k),imt, jmt, 1, imt, 1, jmt,
     &                  xt, yt, init_h, work_H, lenw_h,spval)
           init_h=0 ! we can use the same interpolation weights
          endif
         enddo
c
c       extrapolate to abyssal ocean
c
         print*,' interpolating vertically'
         call sub_flush(6)
         ntr(:,:,nz+1)=ntr(:,:,nz)
         do j=1,jmt
          call ctf (ntr(:,j,:), imt, nz+1, xt, z, 
     &              tr(:,j,:,m),imt, km, 1, imt, 1, km,
     &              xt, zt/100., init_v, work_v, lenw_v,spval)
          init_v=0
         enddo
        enddo ! m

        deallocate(work_v)
        if (.not.rotated_grid) deallocate(work_h)

        print*,' try to fill gaps horizontally  '
        call sub_flush(6)
        do k=1,km
         do m=1,number_tr
          call fillgaps(imt, jmt, km, k, kmt, tr(:,:,k,m), 
     &                     spval, 0,.false.)
         enddo
        enddo
        print*,' try to fill gaps vertically  '
        call sub_flush(6)
        do j=1,jmt
         do m=1,number_tr
          call fillholes(imt, jmt, km, j, kmt, tr(:,j,:,m), 
     &                      spval, 1,.true.)
         enddo
        enddo
        print*,' successfull interpolation for month ',n
        call sub_flush(6)

        ! store january values as initial conditions
        if (n==1) call write_ic_to_cdf(tr,'initial')

        ! write monthly values to a temporary file
        call write_ts_to_file(n,tr)

        ! write monthly mean values to netcdf file for BLUE
        if (enable_blue) call write_blue_to_cdf(
     &                  tr(1,1,1,1),tr(1,1,1,2),nt,n)

       enddo ! n

       ! store also mean values in netcdf file
       ! assuming equal months here, which is a good approximation
       deallocate(ntr)
       allocate(ntr(imt,jmt,km)); ntr=0.
       tr=0
       if (enable_blue_mean)  then
        do n=1,nt
         do m=1,number_tr
          call read_atr_from_file(n,m,ntr)
          do k=1,km; do j=1,jmt; do i=1,imt 
           if (ntr(i,j,k) /=spval ) then
            tr(i,j,k,m)=tr(i,j,k,m)+ntr(i,j,k)
           else
            tr(i,j,k,m)=spval
           endif
          enddo; enddo; enddo
         enddo
        enddo
        do k=1,km; do j=1,jmt; do i=1,imt 
         do m=1,number_tr
          if (tr(i,j,k,m) /=spval ) tr(i,j,k,m)=tr(i,j,k,m)/nt
         enddo
        enddo; enddo; enddo
c        call write_ic_to_cdf(tr,'mean')
        call sub_flush(6)
        ! write annual mean values to netcdf file for BLUE mean
        call write_blue_mean_to_cdf(tr(1,1,1,1),tr(1,1,1,2))
       endif

       deallocate(dtr,x,tx,y,z,ntr,tr )
       end subroutine prep_ic



       subroutine prep_sbc_tau
c
c--------------------------------------------------------
c      prepare surface boundary conditions for momentum
c--------------------------------------------------------
c
       use prep_module
       implicit none
       integer nx,ny,nt,ntr
       real, allocatable :: tau(:,:,:,:),x(:),y(:),mask(:,:)
       real, allocatable :: vec(:,:,:),nvec(:,:,:)
       real, allocatable :: time(:),tx(:)
       real :: toffset
       integer :: n,k,m
       real, allocatable :: ustar(:,:,:),nustar(:,:)
       logical verbose
       integer :: init_t = 1,lenw_t
       real , allocatable :: work_t(:)

       call init_forcing_file('forcing_wind.cdf')

       print*,' reading wind stress data '
       call read_dims_tau(nx,ny,nt)
       print*,' nx=',nx,' ny=',ny,' nt=',nt

       allocate(tau(nx,ny,nt,2),x(nx),y(ny),mask(nx,ny),time(nt) )
       tau=0.;x=0.;y=0.
       call read_tau(tau,x,y,mask,nx,ny,nt,spval,time,toffset)
       call sub_flush(6)
       allocate( vec(nx,ny,2), nvec(imt,jmt,2) ); vec=0.;nvec=0.

       print*,' interpolating data '
       call sub_flush(6)
       verbose=.true.
       if (.not.rotated_grid) then
          init_t=1
          allocate(tx(nx))
          lenw_t = 10*max(imt,jmt)
          allocate( work_t(lenw_t) )
       endif
       do n=1,nt
        call sub_flush(6)
        vec(:,:,1)=tau(:,:,n,1)
        vec(:,:,2)=tau(:,:,n,2)
        if (rotated_grid) then
          call rot_intrp_vctr(vec, x, y, nx, ny, 
     &                      nvec, xu, yu ,
     &              imt,jmt, psir, thetar, phir,spval,verbose)
           verbose=.false.
        else
c         use cheaper subroutine here
           call tranlon (vec(:,:,1), nx, nx-2, ny, x, xu,imt, tx)
           call tranlon (vec(:,:,2), nx, nx-2, ny, x, xu,imt, tx)
           call ctf (vec(:,:,1), nx, ny, tx, y, 
     &                  nvec(:,:,1),imt, jmt, 1, imt, 1, jmt,
     &                  xu, yu, init_t, work_t, lenw_t,spval)
           init_t=0 ! we can use the same interpolation weights
           call ctf (vec(:,:,2), nx, ny, tx, y, 
     &                  nvec(:,:,2),imt, jmt, 1, imt, 1, jmt,
     &                  xu, yu, init_t, work_t, lenw_t,spval)
        endif
        call fillgaps(imt, jmt, km, 1, kmu, nvec(1,1,1), 
     &                 spval, 2,.true.)
        call fillgaps(imt, jmt, km, 1, kmu, nvec(1,1,2), 
     &                 spval, 2,.true.)

c        print*,' writing wind stress to output file for n=',n
        call sub_flush(6)
        call write_tau_to_cdf(nvec,nt,n,time,toffset)
       enddo
       deallocate( tau,x,y,mask,vec,nvec, time)
       if (.not.rotated_grid) deallocate(tx,work_t)
       print*,' done '

c#ifdef notdef
       print*,' reading ustar data '
       call sub_flush(6)
       call read_dims_ustar(nx,ny,nt)
       print*,' nx=',nx,' ny=',ny,' nt=',nt
       allocate(ustar(nx,ny,nt),x(nx),y(ny),time(nt) )
       ustar=0.;x=0.;y=0.;time=0.0
       call read_ustar(ustar,x,y,nx,ny,nt,spval,time,toffset)
       allocate(  nustar(imt,jmt)); nustar=0.
       print*,' interpolating data '
       verbose=.true.
       if (.not.rotated_grid) then
          init_t=1
          allocate(tx(nx))
          lenw_t = 10*max(imt,jmt)
          allocate( work_t(lenw_t) )
       endif
       do n=1,nt
        if (rotated_grid) then
          call rot_intrp_sclr (ustar(1,1,n), x, y, nx, ny, 
     &                       nustar, xt, yt,
     &               imt, jmt, psir, thetar, phir,spval,verbose)
          verbose=.false.
        else
c         use cheaper subroutine here
           call tranlon (ustar(1,1,n), nx, nx-2, ny, x, xt,imt, tx)
           call ctf (ustar(1,1,n), nx, ny, tx, y, 
     &                  nustar,imt, jmt, 1, imt, 1, jmt,
     &                  xt, yt, init_t, work_t, lenw_t,spval)
           init_t=0 ! we can use the same interpolation weights
        endif
        call fillgaps(imt, jmt, km, 1, kmt, nustar, 
     &                     spval, 2,.true.)
        print*,' writing ustar to output file for n=',n
        call write_ustar_to_cdf(nustar,nt,n,time,toffset)
       enddo
       print*,' done '
       if (.not.rotated_grid) deallocate(tx,work_t)
       deallocate( ustar,x,y,nustar,time)
       call sub_flush(6)
c#endif
       end subroutine prep_sbc_tau

      
       subroutine prep_sbc_sflx
c
c--------------------------------------------------------
c      prepare surface boundary conditions for tracer
c--------------------------------------------------------
c
       use prep_module
       implicit none
       integer nx,ny,nt,ntr
       real, allocatable :: x(:),y(:)
       real, allocatable :: time(:),tx(:)
       real :: toffset
       integer :: n,k,m
       real, allocatable :: ntr_clim(:,:),ntr_rest(:,:)
       real, allocatable :: tr_clim(:,:,:),tr_rest(:,:,:)
       real, allocatable :: ntr_flux(:,:), tr_flux(:,:,:)
       logical verbose
       integer :: init_t = 1,lenw_t
       real , allocatable :: work_t(:)

       call init_forcing_file('forcing_sflx.cdf')
       if (number_tr>2) then
         call init_forcing_file('forcing_sflx_tracer.cdf')
       endif

       print*,' reading tracer flux data '
       call sub_flush(6)

       do m=1,number_tr

        print*,' Surface boundary conditions for tracer ',m
        call sub_flush(6)
        call read_dims_tr(nx,ny,nt,m)
        print*,' nx=',nx,' ny=',ny,' nt=',nt
        call sub_flush(6)

        allocate(  ntr_clim(imt,jmt) )
        allocate(  ntr_rest(imt,jmt) )
        allocate(  ntr_flux(imt,jmt) )
        ntr_clim=0.; ntr_rest=0.; ntr_flux = 0.
        allocate(tr_clim(nx,ny,nt),tr_rest(nx,ny,nt),
     &           x(nx),y(ny),tr_flux(nx,ny,nt),time(nt) )
        x=0;y=0;tr_clim=0;tr_rest=0; tr_flux=0.; time=0.0

        apply_landmask = .true.
        call read_tr(tr_clim,tr_rest,tr_flux,x,y,
     &                    time,toffset,nx,ny,nt,m)

        print*,' interpolating data '
        call sub_flush(6)
        verbose=.true.
        if (.not.rotated_grid) then
          init_t=1
          allocate(tx(nx))
          lenw_t = 10*max(imt,jmt)
          allocate( work_t(lenw_t) )
        endif
        do n=1,nt
         print*,'nt=',n
         call sub_flush(6)

         if (rotated_grid) then
          call rot_intrp_sclr (tr_clim(:,:,n), x, y, nx, ny, 
     &                        ntr_clim(:,:), xt, yt,
     &               imt, jmt, psir, thetar, phir,spval,verbose)
          verbose=.false.
          call rot_intrp_sclr (tr_rest(:,:,n), x, y, nx, ny, 
     &                        ntr_rest(:,:), xt, yt,
     &                imt, jmt, psir, thetar, phir,spval,verbose)
          call rot_intrp_sclr (tr_flux(:,:,n), x, y, nx, ny, 
     &                        ntr_flux(:,:), xt, yt,
     &                imt, jmt, psir, thetar, phir,spval,verbose)
         else
c         use cheaper subroutine here
           call tranlon (tr_clim(1,1,n), nx, nx-2, ny, x, xt,imt, tx)
           call ctf (tr_clim(1,1,n), nx, ny, tx, y, 
     &                  ntr_clim,imt, jmt, 1, imt, 1, jmt,
     &                  xt, yt, init_t, work_t, lenw_t,spval)
           init_t=0 ! we can use the same interpolation weights
           call tranlon (tr_rest(1,1,n), nx, nx-2, ny, x, xt,imt, tx)
           call ctf (tr_rest(1,1,n), nx, ny, tx, y, 
     &                  ntr_rest,imt, jmt, 1, imt, 1, jmt,
     &                  xt, yt, init_t, work_t, lenw_t,spval)
           call tranlon (tr_flux(1,1,n), nx, nx-2, ny, x, xt,imt, tx)
           call ctf (tr_flux(1,1,n), nx, ny, tx, y, 
     &                  ntr_flux,imt, jmt, 1, imt, 1, jmt,
     &                  xt, yt, init_t, work_t, lenw_t,spval)
         endif
         if (apply_landmask) then
         call fillgaps(imt, jmt, km, 1, kmt, ntr_clim(:,:), 
     &                     spval, 2,.true.)
         call fillgaps(imt, jmt, km, 1, kmt, ntr_rest(:,:), 
     &                     spval, 2,.true.)
         call fillgaps(imt, jmt, km, 1, kmt, ntr_flux(:,:), 
     &                     spval, 2,.true.)
         endif
         call write_tr_to_cdf(ntr_clim,ntr_rest,ntr_flux,
     &                                time,toffset,nt,m,n)
        enddo ! n: number of time step
        if (.not.rotated_grid) deallocate(tx,work_t)
        deallocate(tr_clim,tr_rest,tr_flux,x,y)
        deallocate( ntr_clim,ntr_rest,ntr_flux,time)
       enddo ! m=1,number_tr
       print*,' done '
       call sub_flush(6)
       print*,' done '
       call sub_flush(6)
       end subroutine prep_sbc_sflx

      
 
       subroutine prep_obc
c--------------------------------------------------------
c      prepare open boundary conditions
c--------------------------------------------------------
       use prep_module
       implicit none
       integer n

       real psi_south(12,imt)
       real psi_west(12,jmt)

       call init_forcing_file('forcing_obc.cdf')

       print*,' preparing open boundaries'
       call sub_flush(6)


       if (enable_obc_north) then

        print*,' preparing northern open boundaries'
        call sub_flush(6)
c
c      write northern obc,  no data for psi
c
        psi_south=0.
        call read_psi_north(psi_south,kmt(:,jmt-1),imt)
        do n=1,12
         where( kmt(:,jmt-1) ==0 ) psi_south(n,:)=spval
        enddo
        call write_obc_to_cdf(psi_south,'n',12)
       endif

       if (enable_obc_south) then

        print*,' preparing southern open boundaries'
        call sub_flush(6)
c
c      read streamfunction for southern boundary
c
        psi_south=0.
        call read_psi_south(psi_south,kmt(:,2),imt)
        do n=1,12
        where( kmt(:,2) ==0 ) psi_south(n,:)=spval
        enddo
c
c      write southern obc 
c
        call write_obc_to_cdf(psi_south,'s',12)
       endif

       if (enable_obc_east) then

        print*,' preparing eastern open boundaries'
        call sub_flush(6)
        psi_west=0.
c        print*,' not yeyt readdy ';stop
c        call read_psi_east(psi_west,kmt(imt-1,:),jmt)
        do n=1,12
         where( kmt(imt-1,:) ==0 ) psi_west(n,:)=spval
        enddo
        call write_obc_to_cdf(psi_west,'e',12)
       endif

       if (enable_obc_west) then

        print*,' preparing western open boundaries'
        call sub_flush(6)
        psi_west=0.
c        print*,' not yeyt readdy ';stop
c        call read_psi_west(psi_west,kmt(2,:),jmt)
        do n=1,12
         where( kmt(2,:) ==0 ) psi_west(n,:)=spval
        enddo
        call write_obc_to_cdf(psi_west,'w',12)
       endif
       print*,' done preparing open boundaries'
      
       end subroutine prep_obc

c
c     island_module need this mpp stuff
c     which is disabled here anyway
c
      subroutine barrier
      end

      subroutine halt_stop(s)
      character*(*) s
      print*,s
      stop
      end

