#include "CPP_OPTIONS.h"

      subroutine smooth_diff3D (fld_in,nbt_in,mythid)

      IMPLICIT NONE
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "FFIELDS.h"
c#include "ctrl.h"
c#include "ctrl_dummy.h"
c#include "optim.h"
#ifdef ALLOW_AUTODIFF_TAMC
#include "tamc.h"
#include "tamc_keys.h"
#endif /* ALLOW_AUTODIFF_TAMC */

#include "GAD.h"
#include "GMREDI.h"
#include "smooth.h"

c declarations:
      integer i,j,k, bi,bj, iMin,iMax,jMin,jMax 
      integer itlo,ithi, jtlo,jthi
      integer myThid, myIter(nSx,nSy)
      integer mykkey, wc01key
      _RL fld_in(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL fld_tmp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL gT_in(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL gTm1_in(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      integer nbt_in
      integer wc01_di, wc01_dj,  wc01_dk, ii, jj, kk
      integer iloop, ilev_1, ilev_2, ilev_3
      integer max_lev2, max_lev3, key_in
      _RL dTtracerLev_bak(nr)

c for now: useless, because level 3 is recomputed anyway
c but : if level3 was computed during the fwd loop by callung
c 	mdsmooth_diff3D (assumes that it would be called
c	directly by the_main_loop) then I would need to pass key_in
c	as a parameter, with different values for T, S, ...
c	in order not to overwrite the same tape
      key_in=0

c computations:
      jtlo = mybylo(mythid)
      jthi = mybyhi(mythid)
      itlo = mybxlo(mythid)
      ithi = mybxhi(mythid)

      _BEGIN_MASTER( mythid )
       DO k=1,nr
       dTtracerLev_bak(k)=dTtracerLev(k)
       dTtracerLev(k)=1.
       ENDDO
      _END_MASTER( mythid )

c      call smooth_print (0,fld_in,mythid)

      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1,sNy
          DO i=1,sNx
c division by ~sqrt(volume):
c now done in wrapper routines:
c           fld_in(i,j,k,bi,bj)=fld_in(i,j,k,bi,bj)
c     & *sqrt(recip_rA(i,j,bi,bj)*recip_drF(k))
c initialization to 0:
           gT_in(i,j,k,bi,bj)   = 0. _d 0
           gTm1_in(i,j,k,bi,bj)   = 0. _d 0
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      _EXCH_XYZ_RL ( fld_in , myThid )
      _EXCH_XYZ_RL ( gt_in , myThid )
      _EXCH_XYZ_RL ( gtm1_in , myThid )

#ifdef ALLOW_TAMC_CHECKPOINTING

c checkpointing:
      max_lev3=nbt_in/(nchklev_1*nchklev_2)+1
      max_lev2=nbt_in/nchklev_1+1
#ifdef ALLOW_AUTODIFF_TAMC
CADJ INIT tape_wc01_lev3 = USER
#endif /* ALLOW_AUTODIFF_TAMC */
      do ilev_3 = 1,nchklev_3
        if(ilev_3.le.max_lev3) then
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE fld_in = tape_wc01_lev3 ,
CADJ & key = key_in*max_lev3 + ilev_3
CADJ STORE gTm1_in = tape_wc01_lev3 ,
CADJ & key = key_in*max_lev3 + ilev_3
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ INIT tape_wc01_lev2 = USER
#endif /* ALLOW_AUTODIFF_TAMC */
      do ilev_2 = 1,nchklev_2
         if(ilev_2.le.max_lev2) then
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE fld_in = tape_wc01_lev2 , 
CADJ & key = key_in*nchklev_2 + ilev_2
CADJ STORE gTm1_in = tape_wc01_lev2 , 
CADJ & key = key_in*nchklev_2 + ilev_2
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ INIT tape_wc01_lev1  = COMMON,
CADJ & nchklev_1*nsx*nsy*nthreads_chkpt
#endif /* ALLOW_AUTODIFF_TAMC */
      do ilev_1 = 1,nchklev_1
      iloop = (ilev_2 - 1)*nchklev_1 + ilev_1
     &            + (ilev_3 - 1)*nchklev_2*nchklev_1
      if ( iloop .le. nbt_in ) then

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE gTm1_in(:,:,:,bi,bj) = tape_wc01_lev1 ,  
CADJ & key = key_in*nchklev_1 + ilev_1
#endif /* ALLOW_AUTODIFF_TAMC */


#else /* ALLOW_TAMC_CHECKPOINTING */
      do iloop=1,nbt_in
#endif

      DO bj=jtlo,jthi
       DO bi=itlo,ithi
        DO k=1,Nr
         DO j=1,sNy
          DO i=1,sNx
           gT_in(i,j,k,bi,bj)   = 0. _d 0
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      _EXCH_XYZ_RL ( gt_in , myThid )

c compute gT_in:
cc     ikey_dynamics=iloop 
      CALL smooth_rhs( fld_in, gT_in, myThid )

      DO bj=jtlo,jthi
       DO bi=itlo,ithi
c adams bashfort on gT_in:
      myIter(bi,bj)=iloop-1
c      print*,"myIter",myIter(bi,bj)
        DO k=1,Nr
        CALL ADAMS_BASHFORTH2(
     &                        bi, bj, k,
     &                        gT_in, gTm1_in,
     &                        myIter(bi,bj) , myThid )
        ENDDO
c time stepping: 
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           fld_in(i,j,k,bi,bj)=fld_in(i,j,k,bi,bj)
     &            +wc01_dt*gT_in(i,j,k,bi,bj)
           gT_in(i,j,k,bi,bj)=0
          ENDDO
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE fld_in(:,:,:,bi,bj) = tape_wc01_lev1,
CADJ & key = key_in*nchklev_1 + ilev_1
#endif /* ALLOW_AUTODIFF_TAMC */

          CALL IMPLDIFF(
     I         bi, bj, 1, sNx, 1, sNy ,
     I         1, kappaRwc01(1-OLx,1-OLy,1,bi,bj),
     I         recip_hFacC,
     U         fld_in,
     I         myThid )

       ENDDO
      ENDDO

      _EXCH_XYZ_RL ( fld_in , myThid )
      _EXCH_XYZ_RL ( gt_in , myThid )
      _EXCH_XYZ_RL ( gtm1_in , myThid )

c      call smooth_print (iloop,fld_in,mythid)

#ifdef ALLOW_TAMC_CHECKPOINTING
      endif
      enddo
      endif
      enddo
      endif
      enddo
#else /* ALLOW_TAMC_CHECKPOINTING */
      enddo 
#endif 

      _BEGIN_MASTER( mythid )
       DO k=1,nr
       dTtracerLev(k)=dTtracerLev_bak(k)
       ENDDO
      _END_MASTER( mythid )

      END