#include "options.inc"


c
c ------------------------------------------------------
c     fold the forcing for SPFLAME with this programm
c     fully automatic if you want. 
c     However, there is one command line parameter: imtf
c                              aug 2001  c eden
c ------------------------------------------------------
c

      program prep_fold
      implicit none
      character (len=80) :: infile,outfile
      integer, allocatable ::  kmt(:,:)
      integer :: imt,jmt,km,imtf, imtf_given,m
      integer   :: ncid,iret,n,ncout
      integer   :: xtid,ytid,ztid,kmtid
      integer   :: start(4),count(4),dims(4)
      character (len=80) :: name
      integer :: ret
#include "netcdf.inc"

      print*,' Folding forcing data'

      imtf_given=-1
      call  get_commandline_arg_int(1,imtf_given,iret)

      if (iret/=0) then
       print*,' you might want to preset imtf '
       imtf_given=-1
      else
       print*,' read imtf =',imtf_given,' from command line'
      endif

      call  get_commandline_arg(2,infile,iret)
      if (iret/=0) infile='forcing.cdf'
      call  get_commandline_arg(3,outfile,iret)
      if (iret/=0) outfile='forcing_folded.cdf'

      print*,' reading from file ',infile(1:len_trim(infile))
      print*,' writing to file ',outfile(1:len_trim(outfile))

      iret=nf_open(infile,NF_noWRITE,ncid)
      ncout = nccre (outfile, NCCLOB, iret)
      iret=nf_set_fill(ncout, NF_NOFILL, iret)

c     read grid 

      iret=nf_inq_varid(ncid,'xt',xtid)
      iret=nf_inq_varid(ncid,'yt',ytid)
      iret=nf_inq_varid(ncid,'zt',ztid)

      iret=NF_INQ_VARDIMID (NCID,xtid,dims)
      iret=NF_INQ_DIMLEN (NCID, dims(1), imt)
      iret=NF_INQ_VARDIMID (NCID,ytid,dims)
      iret=NF_INQ_DIMLEN (NCID, dims(1), jmt)
      iret=NF_INQ_VARDIMID (NCID,ztid,dims)
      iret=NF_INQ_DIMLEN (NCID, dims(1), km)

      print*,'imt=',imt,' jmt=',jmt,' km=',km

c     determine possible folding


      iret=nf_inq_varid(ncid,'kmt',kmtid)
      if (iret == 0) then
       allocate( kmt(imt,jmt) )
       start=(/1,1,1,1/); count=(/imt,jmt,1,1/)
       iret= nf_get_vara_int(ncid,kmtid ,start,count,kmt)
       call try_fold(kmt,imt,jmt,imtf)
      else
        print*,' Cannot find variable kmt in file'
        print*,' NetCDF error message:', nf_strerror(iret)
      endif

      print*,' determined minimal possible folding =',imtf
      if (imtf_given > 0) then
        print*,' however setting imtf to ',imtf_given
        imtf=imtf_given
      endif

c
c     define x-axis in new file
c
      iret = ncddef(ncout, 'Longitude_t',imtf,iret)
      iret = ncddef(ncout, 'Longitude_u',imtf,iret)
      call move_var_cdf(ncid,ncout,'Longitude_t','spv',ret)
      call move_var_cdf(ncid,ncout,'Longitude_u','spv',ret)
c
c     grid
c
      call move_var_cdf(ncid,ncout,'xt','spv',ret)
      call move_var_cdf(ncid,ncout,'xu','spv',ret)
      call move_var_cdf(ncid,ncout,'yt','spv',ret)
      call move_var_cdf(ncid,ncout,'yu','spv',ret)
      call move_var_cdf(ncid,ncout,'zt','spv',ret)
      call move_var_cdf(ncid,ncout,'zw','spv',ret)
      call move_var_cdf(ncid,ncout,'dxtdeg','spv',ret)
      call move_var_cdf(ncid,ncout,'dxudeg','spv',ret)
      call move_var_cdf(ncid,ncout,'dytdeg','spv',ret)
      call move_var_cdf(ncid,ncout,'dyudeg','spv',ret)
      call move_var_cdf(ncid,ncout,'dzt','spv',ret)
      call move_var_cdf(ncid,ncout,'dzw','spv',ret)
      call move_var_cdf(ncid,ncout,'glat','spv',ret)
c
c     topography
c
      call move_var_cdf(ncid,ncout,'kmt','max',ret)
      call move_var_cdf(ncid,ncout,'topo','max',ret)
#ifdef partial_cell
      call move_var_cdf(ncid,ncout,'htp','max',ret)
#endif
      call move_var_cdf(ncid,ncout,'rose','max',ret)
      call move_var_cdf(ncid,ncout,'raw_rose','max',ret)
c
c     initial cond.
c
      call move_var_cdf(ncid,ncout,'temp_ic','spv',ret)
      call move_var_cdf(ncid,ncout,'salt_ic','spv',ret)
      call move_var_cdf(ncid,ncout,'temp_mean','spv',ret)
      call move_var_cdf(ncid,ncout,'salt_mean','spv',ret)
      m=3
 20    write(name,'("tracer_",i2,"_ic")') m
       call replace_space_zero(name)
       call move_var_cdf(ncid,ncout,name,'spv',ret)
      if (ret==0) then; m=m+1; goto 20; endif
      m=3
 25   write(name,'("tracer_",i2,"_mean")') m
       call replace_space_zero(name)
       call move_var_cdf(ncid,ncout,name,'spv',ret)
      if (ret==0) then; m=m+1; goto 25; endif
c
c     s b c
c
      call move_var_cdf(ncid,ncout,'taux','spv',ret)
      call move_var_cdf(ncid,ncout,'tauy','spv',ret)
      call move_var_cdf(ncid,ncout,'sst_clim','spv',ret)
      call move_var_cdf(ncid,ncout,'sst_rest','spv',ret)
      call move_var_cdf(ncid,ncout,'sst_flux','spv',ret)
      call move_var_cdf(ncid,ncout,'sss_clim','spv',ret)
      call move_var_cdf(ncid,ncout,'sss_rest','spv',ret)
      call move_var_cdf(ncid,ncout,'sss_flux','spv',ret)
      call move_var_cdf(ncid,ncout,'ustar','spv',ret)
      m=3
 30    write(name,'("str_",i2,"_clim")') m
       call replace_space_zero(name)
       call move_var_cdf(ncid,ncout,name,'spv',ret)
       write(name,'("str_",i2,"_rest")') m
       call replace_space_zero(name)
       call move_var_cdf(ncid,ncout,name,'spv',ret)
       write(name,'("str_",i2,"_flux")') m
       call replace_space_zero(name)
       call move_var_cdf(ncid,ncout,name,'spv',ret)
      if (ret==0) then; m=m+1; goto 30; endif
c
c     obc
c
      call move_var_cdf(ncid,ncout,'temp_obc_south','spv',ret)
      call move_var_cdf(ncid,ncout,'salt_obc_south','spv',ret)
      call move_var_cdf(ncid,ncout,'psi_obc_south','spv',ret)
      call move_var_cdf(ncid,ncout,'temp_obc_north','spv',ret)
      call move_var_cdf(ncid,ncout,'salt_obc_north','spv',ret)
      call move_var_cdf(ncid,ncout,'psi_obc_north','spv',ret)
      call move_var_cdf(ncid,ncout,'temp_obc_east','spv',ret)
      call move_var_cdf(ncid,ncout,'salt_obc_east','spv',ret)
      call move_var_cdf(ncid,ncout,'psi_obc_east','spv',ret)
      call move_var_cdf(ncid,ncout,'temp_obc_west','spv',ret)
      call move_var_cdf(ncid,ncout,'salt_obc_west','spv',ret)
      call move_var_cdf(ncid,ncout,'psi_obc_west','spv',ret)
      m=3
 40    write(name,'("tr",i2,"_obc_north")') m
       call replace_space_zero(name)
       call move_var_cdf(ncid,ncout,name,'spv',ret)
      if (ret==0) then; m=m+1; goto 40; endif
      m=3
 50    write(name,'("tr",i2,"_obc_south")') m
       call replace_space_zero(name)
       call move_var_cdf(ncid,ncout,name,'spv',ret)
      if (ret==0) then; m=m+1; goto 50; endif
      m=3
 60    write(name,'("tr",i2,"_obc_east")') m
       call replace_space_zero(name)
       call move_var_cdf(ncid,ncout,name,'spv',ret)
      if (ret==0) then; m=m+1; goto 60; endif
      m=3
 70    write(name,'("tr",i2,"_obc_west")') m
       call replace_space_zero(name)
       call move_var_cdf(ncid,ncout,name,'spv',ret)
      if (ret==0) then; m=m+1; goto 70; endif
c
c     if blue density is in the file
c
      call move_var_cdf(ncid,ncout,'temp_blue','spv',ret)
      call move_var_cdf(ncid,ncout,'salt_blue','spv',ret)
c
c     if blue mean density is in the file
c
      call move_var_cdf(ncid,ncout,'temp_blue_mean','spv',ret)
      call move_var_cdf(ncid,ncout,'salt_blue_mean','spv',ret)
c
c     check if BB forcing files are present
c
      call move_var_cdf(ncid,ncout,'tauxBB','spv',ret)
      call move_var_cdf(ncid,ncout,'tauyBB','spv',ret)
      call move_var_cdf(ncid,ncout,'hflxBB','spv',ret)
      call move_var_cdf(ncid,ncout,'sstf_BB','spv',ret)
      call move_var_cdf(ncid,ncout,'sstm_BB','spv',ret)
c
c     sponge layers need special treament
c 
      call move_sponges(ncid,ncout,imt,imtf,jmt,km)

      call ncclos (ncid, iret)
      call ncclos (ncout, iret)

      end program prep_fold


      subroutine move_sponges(ncin,ncout,imt,imtf,jmt,km)
      implicit none
      integer :: ncin,ncout,imt,imtf,jmt,km,nt
      integer :: n,nr_spg,dims(4),nr_nspg,id,natts,k,m
      character (len=2) :: number,nr2
      character (len=80) :: name,unit,s
      integer :: is,ie,js,je,tid,sid,tsid,iret
      integer :: t1id,s1id,ts1id,t2id,s2id,ts2id
      integer :: start(4),count(4),x1dim,x2dim,ydim,zdim,tdim
      real (kind=4), allocatable  :: buf(:,:,:,:,:)
      integer :: ret ,i
      character (len=4) :: tr_name(99)
#include "netcdf.inc"

      iret= nf_get_att_int(ncin,nf_global,'number_of_sponge_layers',
     &                     nr_spg  )
      if (iret/=0) nr_spg = 0
c
c     set tracer names
c
      do m=1,99
        if (m==1) then
         tr_name(m) = 'temp'
        elseif (m==2) then
         tr_name(m) = 'salt'
        else
         write(tr_name(m),'("tr",i2)') m
         do i=1,len_trim(tr_name(m))
          if (tr_name(m)(i:i)==' ')tr_name(m)(i:i)='0'
         enddo
        endif
      enddo

      nr_nspg=nr_spg
      do n=1,nr_spg

       print*,' processing sponge layer #',n

       write(number,'(i2)') n
       if (number(1:1) == ' ') number(1:1)='0'
c
c      get range of sponge layer from temperature
c
       iret=nf_inq_varid(ncin,'temp_spg_'//number,tid)
       if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       iret=NF_INQ_VARDIMID (ncin,tid,dims)
       if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       iret=NF_INQ_DIMLEN (ncin, dims(3), km)
       if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       iret= nf_get_att_int(ncin,tid,'range',dims )
       if (iret/=0) then; print*,nf_strerror(iret);stop;endif
       is=dims(1);ie=dims(2);js=dims(3);je=dims(4)

       print*,' is=',is,' ie=',ie,' js=',js,' je=',je,'km=',km



       if (ie <= imtf ) then
c
c       just transfer
c
        print*,' we can just transfer this sponge layer'

        m=1
 20     continue
         call move_var_cdf(ncin,ncout,tr_name(m)//'_spg_'//number,
     &                    'spv',ret)
         call move_var_cdf(ncin,ncout,tr_name(m)//'_tscl_spg_'//number,
     &                    'spv',ret)
        if (ret==0) then; m=m+1; goto 20; endif

       elseif (is>=imtf) then
c
c       just transfer and adjust is and ie
c
        print*,' we can transfer this sponge layer'
        print*,' and adjust is and ie '

        m=1
 30     continue

         call move_var_cdf(ncin,ncout,tr_name(m)//'_spg_'//number,
     &                    'spv',ret)
         call move_var_cdf(ncin,ncout,tr_name(m)//'_tscl_spg_'//number,
     &                    'spv',ret)

         if (ret==0) then
          dims(1)=is-imtf+2
          dims(2)=ie-imtf+2
c
          iret=nf_inq_varid(ncout,tr_name(m)//'_spg_'//number,tid)
          if (iret/=0) then; print*,' cannot find var.';stop;endif
          iret=nf_inq_varid(ncout,tr_name(m)//'_tscl_spg_'//number,tsid)
          if (iret/=0) then; print*,' cannot find var.';stop;endif

          iret = nf_redef(ncout)
          iret=nf_set_fill(ncout, NF_NOFILL, iret)

          iret= nf_put_att_int(ncout,tid,'range',nf_int,4,dims)
          iret= nf_put_att_int(ncout,tsid,'range',nf_int,4,dims)

          iret = nf_enddef(ncout)
c
c         adjust maybe also the x-axis of sponge layer
c
         endif

        if (ret==0) then; m=m+1; goto 30; endif

       else
c
c       we have to split the sponge layer into two parts
c
        print*,' this sponge layer has to be to splitted '
        nr_nspg=nr_nspg+1
        write(nr2,'(i2)') nr_nspg
        if (nr2(1:1) == ' ') nr2(1:1)='0'

        m=1
 40     continue

c        time dimension and y-axis can be shared by both
         call move_var_cdf(ncin,ncout,
     &       'time_spg_'//number//'_'//tr_name(m),'spv',ret)
         if (ret/=0) goto 50  ! no more tracers

         if (m==1) then
          call move_var_cdf(ncin,ncout,'Lat_spg_'//number,'spv',ret)
         endif
c
         iret=nf_inq_varid(ncin,tr_name(m)//'_spg_'//number,tid)
         print*,' CHECK for depth axis'
c  check for depth axis
         iret=NF_INQ_VARDIMID(ncin,tid,dims)
         s='';iret=NF_INQ_DIMNAME(NCin,dims(3), s) ! depth axis
         print*,'name of depth axis',s
         iret=nf_inq_varid(ncout,s,tsid)
         if (iret/=0) then
           call move_var_cdf(ncin,ncout,s,'spv',ret)
         else
          print*,' depth axis already defined'
         endif

         iret=nf_inq_varid(ncin,tr_name(m)//'_tscl_spg_'//number,tsid)

         iret= NF_INQ_DIMLEN (NCin, DIMs(4), nt)

         allocate( buf(is:ie,js:je,km,nt,2) )

         start=1; count=(/ie-is+1,je-js+1,km,nt/)
         iret= nf_get_vara_real(ncin,tid,start,count,buf(:,:,:,:,1))
         iret= nf_get_vara_real(ncin,tsid,start,count,buf(:,:,:,:,2))

         iret = nf_redef(ncout)
         iret=nf_set_fill(ncout, NF_NOFILL, iret)
        
         if (m==1) then
          x1dim = ncddef(ncout, 'Lon_spg_'//number,imtf-1-is+1 , iret)
          x2dim = ncddef(ncout, 'Lon_spg_'//nr2   ,ie-imtf+1 , iret)
         else
          iret=NF_INQ_DIMID(ncout,'Lon_spg_'//number,x1dim)
          iret=NF_INQ_DIMID(ncout,'Lon_spg_'//nr2,   x2dim)
         endif
         iret=NF_INQ_DIMID(ncout,'Lat_spg_'//number,ydim)
         iret=NF_INQ_DIMID(ncout,
     &                'time_spg_'//number//'_'//tr_name(m),tdim)
         iret=NF_INQ_DIMID(ncout,'depth_t',zdim)

         dims=(/x1dim,ydim,zdim,tdim/)
         t1id = ncvdef (ncout,tr_name(m)//'_spg_'//number,
     &                  ncfloat,4,dims,iret)
         ts1id = ncvdef (ncout,tr_name(m)//'_tscl_spg_'//number,
     &                   ncfloat,4,dims,iret)

         dims=(/x2dim,ydim,zdim,tdim/)
         t2id = ncvdef (ncout,tr_name(m)//'_spg_'//nr2,
     &                  ncfloat,4,dims,iret)
         ts2id = ncvdef (ncout,tr_name(m)//'_tscl_spg_'//nr2,
     &                   ncfloat,4,dims,iret)

         iret =  NF_INQ_VARNATTS (NCin, tid, natts)
         do k=1,natts
          name='';iret= NF_INQ_ATTNAME (NCin, tid, k, name)
          iret= NF_COPY_ATT(NCIN, tid, name, NCOUT, t1id)
          iret= NF_COPY_ATT(NCIN, tid, name, NCOUT, t2id)
         enddo

         iret =  NF_INQ_VARNATTS (NCin, tsid, natts)
         do k=1,natts
          name='';iret= NF_INQ_ATTNAME (NCin, tsid, k, name)
          iret= NF_COPY_ATT(NCIN, tsid, name, NCOUT, ts1id)
          iret= NF_COPY_ATT(NCIN, tsid, name, NCOUT, ts2id)
         enddo

         dims=(/is,imtf-1,js,je/)
         iret= nf_put_att_int(ncout,t1id,'range',nf_int,4,dims)
         iret= nf_put_att_int(ncout,ts1id,'range',nf_int,4,dims)

         dims=(/imtf-imtf+2,ie-imtf+2,js,je/)
         iret= nf_put_att_int(ncout,t2id,'range',nf_int,4,dims)
         iret= nf_put_att_int(ncout,ts2id,'range',nf_int,4,dims)
c
c        maybe care also about x-axis but not really needed
c
         iret = nf_enddef(ncout)

         start=1; count=(/imtf-1-is+1,je-js+1,km,nt/)
         iret= nf_put_vara_real(ncout,t1id ,start,count,
     &                          buf(is:imtf-1,:,:,:,1)) 
         iret= nf_put_vara_real(ncout,ts1id ,start,count,
     &                          buf(is:imtf-1,:,:,:,2)) 

         start=1; count=(/ie-imtf+1,je-js+1,km,nt/)
         iret= nf_put_vara_real(ncout,t2id ,start,count,
     &                          buf(imtf:ie,:,:,:,1)) 
         iret= nf_put_vara_real(ncout,ts2id ,start,count,
     &                          buf(imtf:ie,:,:,:,2)) 

         deallocate(buf)

 50      continue

        if (ret==0) then; m=m+1; goto 40; endif

       endif ! transfer decision

       print*,' '
       print*,' Done processing sponge layer #',number
       print*,' '

      enddo ! loop over sponge layers
c
c     the new nr of sponges
c
      iret = nf_redef(ncout)
      iret= nf_put_att_int(ncout,nf_global,'number_of_sponge_layers',
     &                     nf_int,1,nr_nspg  )
      iret = nf_enddef(ncout)

      print*,' new number of sponge layer:',nr_nspg

      end subroutine move_sponges



      recursive subroutine  move_var_cdf(ncin,ncout,name,maxi,ret)
      implicit none
      integer :: ncin,ncout
      character*(*) name
      character*(*) maxi
      integer :: ret 
#include "netcdf.inc"
      integer :: idi,ido,xtype,ndims,dimin(4),dimout(4)
      integer :: natts,n,iret
      character (len=80) :: s
      integer :: start(4),count(4),count2(4),len,aid
      real (kind=4), allocatable  :: buff(:,:,:,:)
      integer, allocatable        :: bufi(:,:,:,:)
      real (kind=8), allocatable  :: bufd(:,:,:,:)
      logical :: axis = .false.
      real (kind=4) :: spval
      logical :: use_max = .false.  ! use maximum criterium to fold

      if (maxi(1:3) == 'max' ) use_max=.true.
      if (maxi(1:3) == 'spv' ) use_max=.false.

      iret=nf_inq_varid(ncin,name,idi)
      if (iret/=0) then
        print*,' WARING: cannot find variable ',
     &          name(1:len_trim(name)),' in input file'
        ret=-1
        return
      else
        print*,' moving variable ',
     &          name(1:len_trim(name)),' to output file'
        if (use_max) print*,' using maximum criterium to fold '
      endif

      iret=NF_INQ_VARTYPE(ncin,idi, xtype)
      iret=NF_INQ_VARNDIMS(ncin,idi, ndims)
      if (ndims>4) then
        print*,' cannot handle more than 4 dimensions'
        print*,' but found ',ndims
        stop
      endif
      iret=NF_INQ_VARDIMID(ncin,idi,dimin)

      do n=1,ndims
       s='';iret=NF_INQ_DIMNAME(NCin,dimin(n), s) 
       iret=NF_INQ_DIMID(NCout,s,dimout(n)) 
       if (iret/=0) then
        print*,' WARNING: cannot find dimension ',s(1:len_trim(s)),
     &        ' in output file'
        print*,' try to define it ...'
        iret= NF_INQ_DIMLEN (NCin, DIMin(n), len)
        print*,' dimension length = ',len
        iret = nf_redef(ncout)
        iret=nf_set_fill(ncout, NF_NOFILL, iret)
        dimout(n) = ncddef(ncout, s, len, iret)
        print*,' try to find axis variables for that dimension'
        iret=nf_inq_varid(ncin,s,aid)
        if (iret==0) then
         print*,' okay found it and transfering it'
         if (s(1:len_trim(s))==name(1:len_trim(name))) then
           print*,' that was already the variable which we want to '
           axis = .true.
         else
           call move_var_cdf(ncin,ncout,s,maxi,ret)
         endif

         s='';iret= nf_get_att_text(ncin,aid ,'ave_period',s)
         if (iret==0) then
           print*,' found also attribute ave_period, transfering also'
           call move_var_cdf(ncin,ncout,s,maxi,ret)
         endif
         s='';iret= nf_get_att_text(ncin,aid ,'start_ave_period',s)
         if (iret==0) then
           print*,' found also attribute start_ave_period,'
     &            //' transfering also'
           call move_var_cdf(ncin,ncout,s,maxi,ret)
         endif
         print*,'okay'         
         print*,''         
        else
         print*,' negative '
        endif
        iret = nf_enddef(ncout)
       endif
      enddo


      iret = nf_redef(ncout)
      iret=nf_set_fill(ncout, NF_NOFILL, iret)

      ido = ncvdef (ncout,name, xtype,ndims,dimout,iret)

      iret =  NF_INQ_VARNATTS (NCin, idi, natts)
      do n=1,natts
       s='';iret= NF_INQ_ATTNAME (NCin, idi, n, s)
       iret= NF_COPY_ATT(NCIN, idi, s, NCOUT, ido)
      enddo

      iret = nf_enddef(ncout)

      start=1; count=1; count2=1;
      do n=1,ndims
       iret= NF_INQ_DIMLEN (NCin, DIMin(n), count(n))
       iret= NF_INQ_DIMLEN (NCout, DIMout(n), count2(n))
      enddo

      if (count(1) > count2(1) .and. ndims>1 ) then
c      missing value must be of type ncfloat !!!
       iret = nf_get_att_real(ncout, ido, 'missing_value', spval) 
       if (iret/=0) then
        print*,' ERROR: cannot find attribute missing value'
        print*,' for variable'
        stop
       endif
      endif

      print*,'    (# of dimensions=',ndims,
     &       ' new dimensions = ',count2,')'

      if (xtype == nclong) then

       allocate( bufi(count(1),count(2),count(3),count(4)) )
       iret= nf_get_vara_int(ncin,idi ,start,count,bufi)

       if (count(1) > count2(1) .and. ndims>1 ) then
        call foldi(bufi,count2(1),count(1),count(2),
     &                  count(3),count(4),spval,use_max)
        print*,'    folding integer variable ',name(1:len_trim(name))
       endif
       iret= nf_put_vara_int(ncout,ido ,start,count2,
     &   bufi(1:count2(1),1:count2(2),1:count2(3),1:count2(4)) )

       deallocate(bufi)

      elseif (xtype == ncfloat) then


       allocate( buff(count(1),count(2),count(3),count(4)) )
       iret= nf_get_vara_real(ncin,idi ,start,count,buff)

       if (count(1) > count2(1) .and. ndims>1 ) then
        call foldf(buff,count2(1),count(1),count(2),
     &                  count(3),count(4),spval,use_max)
        print*,'    folding real variable ',name(1:len_trim(name))
       endif
       iret= nf_put_vara_real(ncout,ido ,start,count2,
     &   buff(1:count2(1),1:count2(2),1:count2(3),1:count2(4)) )

       deallocate(buff)


      elseif (xtype == ncdouble) then

       allocate( bufd(count(1),count(2),count(3),count(4)) )
       iret= nf_get_vara_double(ncin,idi ,start,count,bufd)

       if (count(1) > count2(1) .and. ndims>1 ) then
        call foldd(bufd,count2(1),count(1),count(2),
     &                  count(3),count(4),spval,use_max)
        print*,'    folding double variable ',name(1:len_trim(name))
       endif
       iret= nf_put_vara_double(ncout,ido ,start,count2,
     &   bufd(1:count2(1),1:count2(2),1:count2(3),1:count2(4)) )

       deallocate(bufd)

      else
       print*,' do not know how to handle NCtype ',xtype
       stop
      endif

      ret=0

      end



      subroutine try_fold(kmt,imtuf,jmt,maxfold)
      implicit none
      integer imtuf,jmt,kmt(imtuf,jmt)

      integer maxfold,i,j,l1(imtuf),l2(imtuf),l3(imtuf),mask(imtuf)
      integer m1,m2,imtm1,imt,ii
      logical twin
      character (len=1) :: mark
c
      maxfold=0
      do j=1,jmt-1
c
c transfer j-row and set mask
c
       do i=1,imtuf
        l1(i)=kmt(i,j)
        mask(i)=0
        if(kmt(i,j).gt.0) mask(i)=1
       end do
c
c look for eastern and western boundaries
c
       do i=1,imtuf
        if(mask(i).ne.0)goto 234
       end do
  234  m1=i-1
       do i=imtuf,1,-1
        if(mask(i).ne.0)goto 345
       end do
  345  m2=i+1
c
c determine imt
c
c       imt=imtuf-max(m1,imtuf-m2)+2   ! too soft
       imt=imtuf-(m1+imtuf-m2)+3  ! see Iceland for explanation
c       imt=imtuf-(m1+imtuf-m2)+2  ! without is(ce)land constraint
       imt=max(imt,imtuf/2+2)
       imtm1=imt-1
c
c clear buffers
c
       do i=1,imtuf
        l2(i)=0
        l3(i)=0
       end do
c
c fold l1 to l2
c
       do i=1,imtuf
        l2(i)=l1(i)
       end do
       do i=imtm1,imtuf
        ii=i-imt+2
        l2(ii)=max(l1(i),l1(ii))
       end do
       do i=imt+1,imtuf
        l2(i)=0
       end do
c
c unfold l2 to l3
c
       do i=1,imtuf
        l3(i)=l2(i)
       end do
       l3(1)=l2(imtm1)
       do i=imt,imtuf
        l3(i)=l2(i-imt+2)
       end do
c
c mask l3
c
       do i=1,imtuf
        l3(i)=l3(i)*mask(i)
       end do
c
c compare l1 and l3
c
       twin=.true.
       do i=1,imtuf
        if(l1(i).ne.l3(i))twin=.false.
       end do
       if (.not.twin) then
        write(6,'('' ...trying imt='',i5)')imt
        write(6,'(''     i  mask    l1    l2    l3'')')
        do i=1,imtuf
         mark=' '
         if(l1(i).ne.l3(i))mark='!'
         write(6,'(5i6,1x,a)')i,mask(i),l1(i),l2(i),l3(i),mark
        end do
        stop
       end if
c       write(6,'('' folding at j='',i5,'' possible with imt='',i5)')
c     &        j,imt
       maxfold=max(maxfold,imt)
      end do
c
      write(6,'(/'' overall folding possible with imt='',i5/)')maxfold
c

      end subroutine try_fold



      subroutine foldf(a,imtf,imtuf,jmt,km,nt,spval,use_max)
      implicit none
      integer imtf,imtuf,jmt,km,nt
      real (kind=4) :: a(imtuf,jmt,km,nt)
      integer i,j,k,ii,n
      real (kind=4) :: spval
      logical :: use_max
      do n=1,nt
      do k=1,km
       do j=1,jmt
        do i=imtf-1,imtuf
        ii=i-imtf+2
        if (use_max) then
         a(ii,j,k,n)=max(a(i,j,k,n),a(ii,j,k,n))
        else        
         if (a(ii,j,k,n) == spval) a(ii,j,k,n) = a(i,j,k,n)
        endif
        enddo
        do i=imtf+1,imtuf
         a(i,j,k,n)=0
        end do
        a(imtf,j,k,n)=a(2,j,k,n) ! apply cyclic boundary condition
       end do
      enddo
      enddo

      end



      subroutine foldd(a,imtf,imtuf,jmt,km,nt,spval,use_max)
      implicit none
      integer imtf,imtuf,jmt,km,nt
      real (kind=8) a(imtuf,jmt,km,nt)
      integer i,j,k,ii,n
      real (kind=4) :: spval
      logical :: use_max
      do n=1,nt
      do k=1,km
       do j=1,jmt
        do i=imtf-1,imtuf
        ii=i-imtf+2
        if (use_max) then
         a(ii,j,k,n)=max(a(i,j,k,n),a(ii,j,k,n))
        else
         if (a(ii,j,k,n) == spval) a(ii,j,k,n) = a(i,j,k,n)
        endif
        enddo
        do i=imtf+1,imtuf
         a(i,j,k,n)=0
        end do
        a(imtf,j,k,n)=a(2,j,k,n) ! apply cyclic boundary condition
       end do
      enddo
      enddo
      end



      subroutine foldi(a,imtf,imtuf,jmt,km,nt,spval,use_max)
      implicit none
      integer imtf,imtuf,jmt,km,nt
      integer a(imtuf,jmt,km,nt)
      integer i,j,k,ii,n
      real (kind=4) :: spval
      logical :: use_max
      do n=1,nt
      do k=1,km
       do j=1,jmt
        do i=imtf-1,imtuf
        ii=i-imtf+2
        if (use_max) then
         a(ii,j,k,n)=max(a(i,j,k,n),a(ii,j,k,n))
        else
         if (a(ii,j,k,n) == spval) a(ii,j,k,n) = a(i,j,k,n)
        endif 
        enddo
        do i=imtf+1,imtuf
         a(i,j,k,n)=0
        end do
        a(imtf,j,k,n)=a(2,j,k,n) ! apply cyclic boundary condition
       end do
      enddo
      enddo
      end

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

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

