#include "options.inc"

c
c     prepare forcing for BB here
c
#define prepare_for_BB

c#define prepare_for_BB_for_D3
c correct for a strange point in Bay of Fundi with this option

#define prepare_for_BB_for_E1





      subroutine prep_BB
      use prep_module
      implicit none
#ifdef prepare_for_BB
      integer nx,ny,nt
      real, allocatable :: tau(:,:,:,:),x(:),y(:)
      real, allocatable :: vec(:,:,:),nvec(:,:,:)
      real, allocatable :: ntr_flux(:,:,:)
      real, allocatable :: tr_flux(:,:,:)
      logical verbose
      integer :: n


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

       allocate(tau(nx,ny,nt,2),x(nx),y(ny))
       tau=0.;x=0.;y=0.
       call read_tauBB(tau,x,y,nx,ny,nt,spval)

       allocate( vec(nx,ny,2), nvec(imt,jmt,2) )
       vec=0.;nvec=0.
       print*,' interpolating data '
       verbose=.true.
       do n=1,nt
        print*,'n=',n
        vec(:,:,1)=tau(:,:,n,1)
        vec(:,:,2)=tau(:,:,n,2)
        call rot_intrp_vctr(vec, x, y, nx, ny, 
     &                      nvec, xu, yu ,
     &              imt,jmt, psir, thetar, phir,spval,verbose)
        verbose=.false.

#ifdef prepare_for_BB_for_D3
c
c       Bay of Fundi
c
        nvec(104,208,:)=0.
#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.)

        print*,' writing wind stress to output file for n=',n
        call write_tauBB_to_cdf(nvec,nt,n)
       enddo
       deallocate( tau,vec,nvec)
       print*,' done '

       print*,' reading heat flux data '
       allocate(  ntr_flux(imt,jmt,nt) )
       ntr_flux = 0.
       allocate(tr_flux(nx,ny,nt) )
       x=0;y=0;tr_flux=0.
       call read_trBB(tr_flux,x,y,nx,ny,nt,spval)

       print*,' interpolating data '
       verbose=.true.
       do n=1,nt
        print*,'nt=',n
         call rot_intrp_sclr (tr_flux(:,:,n), x, y, nx, ny, 
     &                        ntr_flux(:,:,n), xt, yt,
     &               imt, jmt, psir, thetar, phir,spval,verbose)
         verbose=.false.
         call fillgaps(imt, jmt, km, 1, kmt, ntr_flux(:,:,n), 
     &                     spval, 2,.true.)
       enddo
       deallocate(tr_flux,x,y)
       print*,' done '

       print*,' writing heat flux to output file '
       call write_trBB_to_cdf(ntr_flux,nt)
       deallocate( ntr_flux)

      print*,' reading sst feedback pattern data '
      call read_sstf_BB_dim(nx,ny,nt)
      print*,' nx=',nx,' ny=',ny,' nt=',nt

       allocate(  ntr_flux(imt,jmt,nt) )
       ntr_flux = 0.
       allocate(tr_flux(nx,ny,nt),x(nx),y(ny) )
       x=0;y=0;tr_flux=0.
       call read_sstf_BB(tr_flux,x,y,nx,ny,nt,spval)

       print*,' interpolating data '
       verbose=.true.
       do n=1,nt
        print*,'nt=',n
         call rot_intrp_sclr (tr_flux(:,:,n), x, y, nx, ny, 
     &                        ntr_flux(:,:,n), xt, yt,
     &               imt, jmt, psir, thetar, phir,spval,verbose)
         verbose=.false.
         call fillgaps(imt, jmt, km, 1, kmt, ntr_flux(:,:,n), 
     &                     spval, 2,.true.)
       enddo
       deallocate(tr_flux,x,y)
       print*,' done '

       print*,' writing sst feeback pattern to output file '
       call write_sstf_BB_to_cdf(ntr_flux,nt)
       deallocate( ntr_flux)



      print*,' reading mean sst data '
      call read_sstm_BB_dim(nx,ny,nt)
      print*,' nx=',nx,' ny=',ny,' nt=',nt

       allocate(  ntr_flux(imt,jmt,nt) )
       ntr_flux = 0.
       allocate(tr_flux(nx,ny,nt),x(nx),y(ny) )
       x=0;y=0;tr_flux=0.
       call read_sstm_BB(tr_flux,x,y,nx,ny,nt,spval)

       print*,' interpolating data '
       verbose=.true.
       do n=1,nt
        print*,'nt=',n
         call rot_intrp_sclr (tr_flux(:,:,n), x, y, nx, ny, 
     &                        ntr_flux(:,:,n), xt, yt,
     &               imt, jmt, psir, thetar, phir,spval,verbose)
         verbose=.false.
         call fillgaps(imt, jmt, km, 1, kmt, ntr_flux(:,:,n), 
     &                     spval, 2,.true.)
       enddo
       deallocate(tr_flux,x,y)
       print*,' done '

       print*,' writing mean sst pattern to output file '
       call write_sstm_BB_to_cdf(ntr_flux,nt)
       deallocate( ntr_flux)

#endif
      end subroutine prep_BB

#ifdef prepare_for_BB

      subroutine read_dims_tauBB(nx,ny,nt)
      implicit none
      integer nx,ny,nt
#include "netcdf.inc"
      character (len=80) :: infile,name
      integer :: iret,ncid,txid,dims(4)
      infile = '/users/model/ceden/dov3/nao_reg/reg_tau_nao.cdf'

      iret=nf_open(infile,NF_nowrite,ncid)
      iret=nf_inq_varid(ncid,'taux_reg',txid)
      iret=NF_INQ_VARDIMID (NCID,txid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nt)
      call ncclos (ncid, iret)
      end subroutine


      subroutine  read_tauBB(tau,x,y,nx,ny,nt,spval)
      implicit none
      integer nx,ny,nt
      real :: tau(nx,ny,nt,2),x(nx),y(ny),spval
#include "netcdf.inc"
      character (len=80) :: infile,name
      integer :: iret,ncid,txid,tyid,latid,lonid,n
      integer :: corner(4),edges(4),dims(4)
#ifdef netcdf_real4
      real (kind=4) :: tau4(nx,ny,nt),x4(nx),y4(ny), spv
#else
      real          :: tau4(nx,ny,nt),x4(nx),y4(ny), spv
#endif

      infile = '/users/model/ceden/dov3/nao_reg/reg_tau_nao.cdf'
      iret=nf_open(infile,NF_nowrite,ncid)
      iret=nf_inq_varid(ncid,'taux_reg',txid)
      iret=nf_inq_varid(ncid,'tauy_reg',tyid)
      iret=NF_INQ_VARDIMID (NCID,txid,DIMS)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),name) 
      iret=NF_INQ_VARID (NCID, name, lonid)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),name) 
      iret=NF_INQ_VARID (NCID, name, latid)

      corner=1; edges(1)=nx; 
      iret= nf_get_vara_real (ncid,lonid ,corner,edges,x4);x=x4
      corner=1; edges(1)=ny; 
      iret= nf_get_vara_real (ncid,latid ,corner,edges,y4);y=y4

      iret = nf_get_att_real(ncid, txid, 'missing_value', spv)

      corner=(/1,1,1,1/); edges=(/nx,ny,nt,1/); 
      iret= nf_get_vara_real (ncid,txid ,corner,edges,tau4)
      tau(:,:,:,1)=tau4
      iret= nf_get_vara_real (ncid,tyid ,corner,edges,tau4)
      tau(:,:,:,2)=tau4
      where( tau == spv) tau=spval
!      convert N/m**2 to dyn/cm**2 
      where( tau /= spval) tau=tau*10.
      call ncclos (ncid, iret)
      end subroutine




      subroutine  read_trBB(hflx,x,y,nx,ny,nt,spval)
      implicit none
      integer nx,ny,nt
      real :: hflx(nx,ny,nt),x(nx),y(ny),spval
#include "netcdf.inc"
      character (len=80) :: infile,name
      integer :: iret,ncid,hflxid,latid,lonid,n
      integer :: corner(4),edges(4),dims(4)
#ifdef netcdf_real4
      real (kind=4) :: hflx4(nx,ny,nt),x4(nx),y4(ny),spv
#else
      real          :: hflx4(nx,ny,nt),x4(nx),y4(ny),spv
#endif
      real,parameter :: dimfac=0.2389e-4

      infile = '/users/model/ceden/dov3/nao_reg/reg_hflx_nao.cdf'
      iret=nf_open(infile,NF_nowrite,ncid)
      iret=nf_inq_varid(ncid,'HFLAN_reg',hflxid)
      iret=NF_INQ_VARDIMID (NCID,hflxid,DIMS)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),name) 
      iret=NF_INQ_VARID (NCID, name, lonid)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),name) 
      iret=NF_INQ_VARID (NCID, name, latid)

      corner=1; edges(1)=nx; 
      iret= nf_get_vara_real (ncid,lonid ,corner,edges,x4);x=x4
      corner=1; edges(1)=ny; 
      iret= nf_get_vara_real (ncid,latid ,corner,edges,y4);y=y4

      iret = nf_get_att_real(ncid,hflxid, 'missing_value', spv)

      corner=(/1,1,1,1/); edges=(/nx,ny,nt,1/); 
      iret= nf_get_vara_real (ncid,hflxid ,corner,edges,hflx4)
      hflx=hflx4
      where( hflx == spv) hflx=spval
      ! convert W/m^2 to cal/s/cm**2
      where( hflx /= spval) hflx=hflx*dimfac
      call ncclos (ncid, iret)
      end subroutine



      subroutine  read_sstf_BB_dim(nx,ny,nt)
      implicit none
      integer nx,ny,nt
#include "netcdf.inc"
      character (len=80) :: infile,name
      integer :: iret,ncid,sstfid,latid,lonid,n
      integer :: corner(4),edges(4),dims(4)
      infile = '/users/model/ceden/dov3/nao_reg/sst_pattern2_ufo.cdf'
      iret=nf_open(infile,NF_nowrite,ncid)
      iret=nf_inq_varid(ncid,'sst',sstfid)
      iret=NF_INQ_VARDIMID (NCID,sstfid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
      iret=NF_INQ_DIMLEN (NCID, DIMs(3), nt)
      call ncclos (ncid, iret)
      end subroutine  read_sstf_BB_dim





      subroutine  read_sstf_BB(sstf,x,y,nx,ny,nt,spval)
      implicit none
      integer nx,ny,nt
      real :: sstf(nx,ny,nt),x(nx),y(ny),spval
#include "netcdf.inc"
      character (len=80) :: infile,name
      integer :: iret,ncid,sstfid,latid,lonid,n
      integer :: corner(4),edges(4),dims(4)
#ifdef netcdf_real4
      real (kind=4) :: sstf4(nx,ny,nt),x4(nx),y4(ny),spv
#else
      real          :: sstf4(nx,ny,nt),x4(nx),y4(ny),spv
#endif

      infile = '/users/model/ceden/dov3/nao_reg/sst_pattern2_ufo.cdf'
      iret=nf_open(infile,NF_nowrite,ncid)
      iret=nf_inq_varid(ncid,'sst',sstfid)
      iret=NF_INQ_VARDIMID (NCID,sstfid,DIMS)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),name) 
      iret=NF_INQ_VARID (NCID, name, lonid)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),name) 
      iret=NF_INQ_VARID (NCID, name, latid)

      corner=1; edges(1)=nx; 
      iret= nf_get_vara_real (ncid,lonid ,corner,edges,x4);x=x4
      corner=1; edges(1)=ny; 
      iret= nf_get_vara_real (ncid,latid ,corner,edges,y4);y=y4

      iret = nf_get_att_real(ncid,sstfid, 'missing_value', spv)

      corner=(/1,1,1,1/); edges=(/nx,ny,1,1/); 
      iret= nf_get_vara_real (ncid,sstfid ,corner,edges,sstf4)
      sstf=sstf4
      do n=2,nt
       sstf(:,:,n)=sstf(:,:,1)
      enddo
      where( sstf == spv) sstf=spval
      call ncclos (ncid, iret)
      end subroutine




      subroutine  read_sstm_BB_dim(nx,ny,nt)
      implicit none
      integer nx,ny,nt
#include "netcdf.inc"
      character (len=80) :: infile,name
      integer :: iret,ncid,sstfid,latid,lonid,n
      integer :: corner(4),edges(4),dims(4)
      infile = 
#ifdef prepare_for_BB_for_D3
     &'/users/model/ceden/dov3/N_atl/DS1_1/sst_mmean_ufo.cdf'
#endif
#ifdef prepare_for_BB_for_E1
     &'/users/model/ceden/dov3/N_atl/E5_bbl/sst_mmean_ufo.cdf'
#endif
      iret=nf_open(infile,NF_nowrite,ncid)
      iret=nf_inq_varid(ncid,'temp',sstfid)
      iret=NF_INQ_VARDIMID (NCID,sstfid,DIMS)
      iret=NF_INQ_DIMLEN (NCID, DIMs(1), nx)
      iret=NF_INQ_DIMLEN (NCID, DIMs(2), ny)
      iret=NF_INQ_DIMLEN (NCID, DIMs(4), nt)
      call ncclos (ncid, iret)
      end subroutine  read_sstm_BB_dim




      subroutine  read_sstm_BB(sstf,x,y,nx,ny,nt,spval)
      implicit none
      integer nx,ny,nt
      real :: sstf(nx,ny,nt),x(nx),y(ny),spval
#include "netcdf.inc"
      character (len=80) :: infile,name
      integer :: iret,ncid,sstfid,latid,lonid,n
      integer :: corner(4),edges(4),dims(4)
#ifdef netcdf_real4
      real (kind=4) :: sstf4(nx,ny,nt),x4(nx),y4(ny),spv
#else
      real          :: sstf4(nx,ny,nt),x4(nx),y4(ny),spv
#endif

      infile = 
#ifdef prepare_for_BB_for_D3
     &'/users/model/ceden/dov3/N_atl/DS1_1/sst_mmean_ufo.cdf'
#endif
#ifdef prepare_for_BB_for_E1
     &'/users/model/ceden/dov3/N_atl/E5_bbl/sst_mmean_ufo.cdf'
#endif
      iret=nf_open(infile,NF_nowrite,ncid)
      iret=nf_inq_varid(ncid,'temp',sstfid)
      iret=NF_INQ_VARDIMID (NCID,sstfid,DIMS)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(1),name) 
      iret=NF_INQ_VARID (NCID, name, lonid)
      name='';iret=NF_INQ_DIMNAME (NCID, DIMs(2),name) 
      iret=NF_INQ_VARID (NCID, name, latid)

      corner=1; edges(1)=nx; 
      iret= nf_get_vara_real (ncid,lonid ,corner,edges,x4);x=x4
      corner=1; edges(1)=ny; 
      iret= nf_get_vara_real (ncid,latid ,corner,edges,y4);y=y4

      iret = nf_get_att_real(ncid,sstfid, 'missing_value', spv)

      corner=(/1,1,1,1/); edges=(/nx,ny,1,nt/); 
      iret= nf_get_vara_real (ncid,sstfid ,corner,edges,sstf4)
      sstf=sstf4
      where( sstf == spv) sstf=spval
      call ncclos (ncid, iret)
      end subroutine




      subroutine write_tauBB_to_cdf(tau,nt,k)
!     write the k.th time level
!     if k=1 then define variables
      use prep_module
      implicit none
#include "netcdf.inc"
      integer nt,n,k
      real tau(imt,jmt,2)
      integer ncid,iret,i
      integer itauxid,itauyid,timedim,timeid,tstrtid,taveid
      integer lon_udim,lat_udim
      integer  corner(4), edges(4),dims(4)
#ifdef netcdf_real4
      real (kind=4) , allocatable :: var(:,:)
#else
      real , allocatable :: var(:,:)
#endif
      integer :: nn,dpm(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
      integer :: toffset = 0*365
      character(len=80) :: name,unit

      iret=nf_open('forcing.cdf',NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)

      if (k==1) then

      print*,' defining variable for tauxBB,tauyBB in NetCDF file'
      iret=nf_redef(ncid)

      i=nt; timedim  = ncddef(ncid, 'time_tauBB', i, iret)

      dims(1)=timedim;i=1
      timeid = ncvdef (ncid,'time_tauBB', NCFLOAT,i,dims,iret)
      tstrtid = ncvdef (ncid,'tstrt_tauBB', NCFLOAT,i,dims,iret)
      taveid = ncvdef (ncid,'tave_tauBB', NCFLOAT,i,dims,iret)

      name = 'Time'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'long_name', NCCHAR, i, name, iret) 

      name = 'Start of averaging period'; i=len_trim(name)
      call ncaptc(ncid, tstrtid, 'long_name', NCCHAR, i, name, iret) 

      name = 'averaging period'; i=len_trim(name)
      call ncaptc(ncid, taveid, 'long_name', NCCHAR, i, name, iret) 

      name = 'days'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid, 'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, taveid, 'units', NCCHAR, i, name, iret) 
      name = '01-JAN-0000 00:00:00'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'time_origin', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid, 'time_origin', NCCHAR, i, name, iret) 

      name = 'no_leap'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'calendar_type',NCCHAR,i,name,iret) 
      name = 'tstrt_tauBB'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'start_ave_period',NCCHAR,i,name,iret) 
      name = 'tave_tauBB'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'ave_period',NCCHAR,i,name,iret) 
c
      iret=nf_inq_dimid(ncid,'Longitude_u',lon_udim)
      iret=nf_inq_dimid(ncid,'Latitude_u',lat_udim)
      dims(1) = Lon_udim; dims(2) = Lat_udim; dims(3) = timedim; i=3
      itauxid   = ncvdef (ncid,'tauxBB', NCFLOAT,i,dims,iret)
      itauyid   = ncvdef (ncid,'tauyBB', NCFLOAT,i,dims,iret)
      name = 'Zonal Wind Stress     '
      unit = 'dyn/cm^2'
      call dvcdf(ncid,itauxid,name,24,unit,16,spval)
      name = 'Meridional Wind Stress     '
      unit = 'dyn/cm^2'
      call dvcdf(ncid,itauyid,name,24,unit,16,spval)
      iret=nf_enddef(ncid)

      allocate(var(nt,1)); var(1,1)=dpm(1)+toffset
      do n=2,nt
       nn=mod(n-1,12)+1; var(n,1)=var(n-1,1)+dpm(nn)
      enddo
      do n=1,nt
       nn=mod(n-1,12)+1; var(n,1)=var(n,1)-dpm(nn)/2.
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, timeid, corner, edges,var, iret)

      var(1,1)=toffset
      do n=2,nt
       nn=mod(n-2,12)+1; 
       var(n,1)=var(n-1,1)+dpm(nn)
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, tstrtid, corner, edges,var, iret)

      do n=1,nt
       nn=mod(n-1,12)+1; 
       var(n,1)=dpm(nn)
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, taveid, corner, edges,var, iret)
 
      deallocate(var)
      endif !k==1

      iret=nf_inq_varid(ncid,'tauxBB',itauxid)
      iret=nf_inq_varid(ncid,'tauyBB',itauyid)

      edges(1)=imt; edges(2)=jmt; corner=1
      allocate(var(imt,jmt) )
      corner(3)=k; edges(3)=1
      var=tau(:,:,1)
      where( var == spval ) var = spval+100.
      where( kmu==0 ) var=spval
      call ncvpt(ncid, itauxid, corner, edges,var, iret)
      var=tau(:,:,2)
      where( var == spval ) var = spval+100.
      where( kmu==0 ) var=spval
      call ncvpt(ncid, itauyid, corner, edges,var, iret)
      call ncclos (ncid, iret)
      deallocate(var)
      end subroutine



      subroutine write_trBB_to_cdf(hflx,nt)
      use prep_module
      implicit none
#include "netcdf.inc"
      integer nt,n
      real hflx(imt,jmt,nt)
      integer ncid,iret,i
      integer hflxid,timedim,timeid,tstrtid,taveid
      integer lon_tdim,lat_tdim
      integer corner(4), edges(4),dims(4)
#ifdef netcdf_real4
      real (kind=4) , allocatable :: var(:,:)
#else
      real , allocatable :: var(:,:)
#endif
      integer :: nn,dpm(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
      integer :: toffset = 0*365
      character(len=80) :: name,unit

      iret=nf_open('forcing.cdf',NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      iret=nf_redef(ncid)

      i=nt; timedim  = ncddef(ncid, 'time_hflxBB', i, iret)

      dims(1)=timedim;i=1
      timeid = ncvdef (ncid,'time_hflxBB', NCFLOAT,i,dims,iret)
      tstrtid = ncvdef (ncid,'tstrt_hflxBB', NCFLOAT,i,dims,iret)
      taveid = ncvdef (ncid,'tave_hflxBB', NCFLOAT,i,dims,iret)

      name = 'Time'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'long_name', NCCHAR, i, name, iret) 

      name = 'Start of averaging period'; i=len_trim(name)
      call ncaptc(ncid, tstrtid, 'long_name', NCCHAR, i, name, iret) 

      name = 'averaging period'; i=len_trim(name)
      call ncaptc(ncid, taveid, 'long_name', NCCHAR, i, name, iret) 

      name = 'days'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid, 'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, taveid, 'units', NCCHAR, i, name, iret) 
      name = '01-JAN-0000 00:00:00'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'time_origin', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid, 'time_origin', NCCHAR, i, name, iret) 

      name = 'no_leap'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'calendar_type',NCCHAR,i,name,iret) 
      name = 'tstrt_hflxBB'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'start_ave_period',NCCHAR,i,name,iret) 
      name = 'tave_hflxBB'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'ave_period',NCCHAR,i,name,iret) 
c
      iret=nf_inq_dimid(ncid,'Longitude_t',lon_tdim)
      iret=nf_inq_dimid(ncid,'Latitude_t',lat_tdim)
      dims(1) = Lon_tdim; dims(2) = Lat_tdim; dims(3) = timedim; i=3
      hflxid   = ncvdef (ncid,'hflxBB', NCFLOAT,i,dims,iret)

      name = 'Heat flux anomaly'
      unit = 'cal/s/cm**2'
      call dvcdf(ncid,hflxid,name,24,unit,16,spval)

      iret=nf_enddef(ncid)

      allocate(var(nt,1)); var(1,1)=dpm(1)+toffset
      do n=2,nt
       nn=mod(n-1,12)+1; var(n,1)=var(n-1,1)+dpm(nn)
      enddo
      do n=1,nt
       nn=mod(n-1,12)+1; var(n,1)=var(n,1)-dpm(nn)/2.
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, timeid, corner, edges,var, iret)

      var(1,1)=toffset
      do n=2,nt
       nn=mod(n-2,12)+1; 
       var(n,1)=var(n-1,1)+dpm(nn)
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, tstrtid, corner, edges,var, iret)

      do n=1,nt
       nn=mod(n-1,12)+1; 
       var(n,1)=dpm(nn)
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, taveid, corner, edges,var, iret)

      edges(1)=imt; edges(2)=jmt; corner=1
      deallocate(var);allocate(var(imt,jmt) )
      do n=1,nt
       corner(3)=n; edges(3)=1
       var=hflx(:,:,n)
       call ncvpt(ncid,hflxid, corner, edges,var, iret)
      enddo
      call ncclos (ncid, iret)
      deallocate(var)
      end




      subroutine write_sstf_BB_to_cdf(hflx,nt)
      use prep_module
      implicit none
#include "netcdf.inc"
      integer nt,n
      real hflx(imt,jmt,nt)
      integer ncid,iret,i
      integer hflxid,timedim,timeid,tstrtid,taveid
      integer lon_tdim,lat_tdim
      integer corner(4), edges(4),dims(4)
#ifdef netcdf_real4
      real (kind=4) , allocatable :: var(:,:)
#else
      real , allocatable :: var(:,:)
#endif
      integer :: nn,dpm(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
      integer :: toffset = 0*365
      character(len=80) :: name,unit

      iret=nf_open('forcing.cdf',NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      iret=nf_redef(ncid)

      i=nt; timedim  = ncddef(ncid, 'time_sstfBB', i, iret)

      dims(1)=timedim;i=1
      timeid = ncvdef (ncid,'time_sstfBB', NCFLOAT,i,dims,iret)
      tstrtid = ncvdef (ncid,'tstrt_sstfBB', NCFLOAT,i,dims,iret)
      taveid = ncvdef (ncid,'tave_sstfBB', NCFLOAT,i,dims,iret)

      name = 'Time'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'long_name', NCCHAR, i, name, iret) 

      name = 'Start of averaging period'; i=len_trim(name)
      call ncaptc(ncid, tstrtid, 'long_name', NCCHAR, i, name, iret) 

      name = 'averaging period'; i=len_trim(name)
      call ncaptc(ncid, taveid, 'long_name', NCCHAR, i, name, iret) 

      name = 'days'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid, 'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, taveid, 'units', NCCHAR, i, name, iret) 
      name = '01-JAN-0000 00:00:00'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'time_origin', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid, 'time_origin', NCCHAR, i, name, iret) 

      name = 'no_leap'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'calendar_type',NCCHAR,i,name,iret) 
      name = 'tstrt_sstfBB'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'start_ave_period',NCCHAR,i,name,iret) 
      name = 'tave_sstfBB'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'ave_period',NCCHAR,i,name,iret) 
c
      iret=nf_inq_dimid(ncid,'Longitude_t',lon_tdim)
      iret=nf_inq_dimid(ncid,'Latitude_t',lat_tdim)
      dims(1) = Lon_tdim; dims(2) = Lat_tdim; dims(3) = timedim; i=3
      hflxid   = ncvdef (ncid,'sstf_BB', NCFLOAT,i,dims,iret)

      name = 'SST feedback pattern'
      unit = 'degC'
      call dvcdf(ncid,hflxid,name,24,unit,16,spval)

      iret=nf_enddef(ncid)

      allocate(var(nt,1)); var(1,1)=dpm(1)+toffset
      do n=2,nt
       nn=mod(n-1,12)+1; var(n,1)=var(n-1,1)+dpm(nn)
      enddo
      do n=1,nt
       nn=mod(n-1,12)+1; var(n,1)=var(n,1)-dpm(nn)/2.
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, timeid, corner, edges,var, iret)

      var(1,1)=toffset
      do n=2,nt
       nn=mod(n-2,12)+1; 
       var(n,1)=var(n-1,1)+dpm(nn)
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, tstrtid, corner, edges,var, iret)

      do n=1,nt
       nn=mod(n-1,12)+1; 
       var(n,1)=dpm(nn)
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, taveid, corner, edges,var, iret)

      edges(1)=imt; edges(2)=jmt; corner=1
      deallocate(var);allocate(var(imt,jmt) )
      do n=1,nt
       corner(3)=n; edges(3)=1
       var=hflx(:,:,n)
       call ncvpt(ncid,hflxid, corner, edges,var, iret)
      enddo
      call ncclos (ncid, iret)
      deallocate(var)
      end




      subroutine write_sstm_BB_to_cdf(hflx,nt)
      use prep_module
      implicit none
#include "netcdf.inc"
      integer nt,n
      real hflx(imt,jmt,nt)
      integer ncid,iret,i
      integer hflxid,timedim,timeid,tstrtid,taveid
      integer lon_tdim,lat_tdim
      integer corner(4), edges(4),dims(4)
#ifdef netcdf_real4
      real (kind=4) , allocatable :: var(:,:)
#else
      real , allocatable :: var(:,:)
#endif
c      integer, parameter :: lsstm = 4
c      integer :: nn,dpm(lsstm) = (/90,91,92,92/)
      integer, parameter :: lsstm = 12
      integer :: nn,dpm(lsstm) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
      integer :: toffset = 0*365
      character(len=80) :: name,unit

      iret=nf_open('forcing.cdf',NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      iret=nf_redef(ncid)

      i=nt; timedim  = ncddef(ncid, 'time_sstmBB', i, iret)

      dims(1)=timedim;i=1
      timeid = ncvdef (ncid,'time_sstmBB', NCFLOAT,i,dims,iret)
      tstrtid = ncvdef (ncid,'tstrt_sstmBB', NCFLOAT,i,dims,iret)
      taveid = ncvdef (ncid,'tave_sstmBB', NCFLOAT,i,dims,iret)

      name = 'Time'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'long_name', NCCHAR, i, name, iret) 

      name = 'Start of averaging period'; i=len_trim(name)
      call ncaptc(ncid, tstrtid, 'long_name', NCCHAR, i, name, iret) 

      name = 'averaging period'; i=len_trim(name)
      call ncaptc(ncid, taveid, 'long_name', NCCHAR, i, name, iret) 

      name = 'days'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid, 'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, taveid, 'units', NCCHAR, i, name, iret) 
      name = '01-JAN-0000 00:00:00'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'time_origin', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid, 'time_origin', NCCHAR, i, name, iret) 

      name = 'no_leap'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'calendar_type',NCCHAR,i,name,iret) 
      name = 'tstrt_sstmBB'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'start_ave_period',NCCHAR,i,name,iret) 
      name = 'tave_sstmBB'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'ave_period',NCCHAR,i,name,iret) 
c
      iret=nf_inq_dimid(ncid,'Longitude_t',lon_tdim)
      iret=nf_inq_dimid(ncid,'Latitude_t',lat_tdim)
      dims(1) = Lon_tdim; dims(2) = Lat_tdim; dims(3) = timedim; i=3
      hflxid   = ncvdef (ncid,'sstm_BB', NCFLOAT,i,dims,iret)

      name = 'mean SST'
      unit = 'degC'
      call dvcdf(ncid,hflxid,name,24,unit,16,spval)

      iret=nf_enddef(ncid)

      allocate(var(nt,1)); var(1,1)=dpm(1)+toffset
      do n=2,nt
       nn=mod(n-1,lsstm)+1; var(n,1)=var(n-1,1)+dpm(nn)
      enddo
      do n=1,nt
       nn=mod(n-1,lsstm)+1; var(n,1)=var(n,1)-dpm(nn)/2.
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, timeid, corner, edges,var, iret)

      var(1,1)=toffset
      do n=2,nt
       nn=mod(n-2,lsstm)+1; 
       var(n,1)=var(n-1,1)+dpm(nn)
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, tstrtid, corner, edges,var, iret)

      do n=1,nt
       nn=mod(n-1,lsstm)+1; 
       var(n,1)=dpm(nn)
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, taveid, corner, edges,var, iret)

      edges(1)=imt; edges(2)=jmt; corner=1
      deallocate(var);allocate(var(imt,jmt) )
      do n=1,nt
       corner(3)=n; edges(3)=1
       var=hflx(:,:,n)
       call ncvpt(ncid,hflxid, corner, edges,var, iret)
      enddo
      call ncclos (ncid, iret)
      deallocate(var)
      end
#endif
