MODULE timeloop_sw_mod
  USE genmod
  USE transfert_mod
  USE etat0_mod
  
  INTEGER,PARAMETER :: euler=1, leapfrog=2, leapfrog_matsuno=3, adam_bashforth=4
  
CONTAINS
  
  SUBROUTINE timeloop
  USE field_mod
  USE domain_mod
!  USE wave_mod
  USE caldyn_sw_mod
  USE write_field
  USE dimensions
  USE geometry
  USE transfert_mod
  USE metric
  USE dissip_sw_mod
  USE ioipsl

  USE etat0_mod
  IMPLICIT NONE
  TYPE(t_field),POINTER :: f_h(:),f_hm1(:), f_hm2(:)
  TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:)
  TYPE(t_field),POINTER :: f_dh(:),f_dhm1(:), f_dhm2(:)
  TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:)

  REAL(rstd),POINTER :: h(:) ,hm1(:), hm2(:)
  REAL(rstd),POINTER :: u(:) , um1(:), um2(:)
  REAL(rstd),POINTER :: dh(:), dhm1(:), dhm2(:)
  REAL(rstd),POINTER :: du(:), dum1(:), dum2(:)
  REAL(rstd) :: dt
  INTEGER :: ind
  INTEGER :: it,i,j,n
  INTEGER :: scheme
  INTEGER :: matsuno_period
  INTEGER :: itaumax
  

  dt=90.
  CALL getin('dt',dt)

  itaumax=100
  CALL getin('itaumax',itaumax)
  
  scheme=leapfrog_matsuno
  CALL getin('scheme',scheme)
  
  matsuno_period=5
  CALL getin('matsuno_period',matsuno_period)
  IF (scheme==leapfrog) matsuno_period=itaumax+1
 
 
  CALL allocate_field(f_h,field_t,type_real)
  CALL allocate_field(f_hm1,field_t,type_real)
  CALL allocate_field(f_hm2,field_t,type_real)
  CALL allocate_field(f_dh,field_t,type_real)
  CALL allocate_field(f_dhm1,field_t,type_real)
  CALL allocate_field(f_dhm2,field_t,type_real)

  CALL allocate_field(f_u,field_u,type_real)
  CALL allocate_field(f_um1,field_u,type_real)
  CALL allocate_field(f_um2,field_u,type_real)
  CALL allocate_field(f_du,field_u,type_real)
  CALL allocate_field(f_dum1,field_u,type_real)
  CALL allocate_field(f_dum2,field_u,type_real)


  CALL allocate_caldyn
  CALL init_dissip(dt)

  CALL etat0_williamson(f_h,f_u)
  
  DO it=0,itaumax
    PRINT *,"SW: It No :",It

    CALL caldyn(f_h, f_u, f_dh, f_du)
    

    IF (scheme==Euler) THEN
      CALL  euler_scheme
    ELSE IF (scheme==leapfrog) THEN
      CALL leapfrog_scheme
    ELSE IF (scheme==leapfrog_matsuno) THEN
      CALL  leapfrog_matsuno_scheme
    ELSE IF (scheme==adam_bashforth) THEN
      CALL dissip(f_u,f_du)   
      CALL adam_bashforth_scheme
    ENDIF

  ENDDO
  
  CONTAINS

    SUBROUTINE Euler_scheme
    IMPLICIT NONE
      INTEGER :: ind
      
      DO ind=1,ndomain
        h=f_h(ind) ; u=f_u(ind) ; 
        dh=f_dh(ind) ; du=f_du(ind) ;
        h(:)=h(:)+dt*dh(:)
        u(:)=u(:)+dt*du(:)
      ENDDO
      
    END SUBROUTINE Euler_scheme
    
    SUBROUTINE leapfrog_scheme
    IMPLICIT NONE
      INTEGER :: ind

      DO ind=1,ndomain
        h=f_h(ind)   ; u=f_u(ind)   ; 
        hm1=f_hm1(ind) ; um1=f_um1(ind) ; 
        hm2=f_hm2(ind) ; um2=f_um2(ind) ; 
        dh=f_dh(ind) ; du=f_du(ind) ;

        IF (it==0) THEN
          hm2(:)=h(:) ; um2(:)=u(:) 

          h(:)=h(:)+dt*dh(:)
          u(:)=u(:)+dt*du(:)
          hm1(:)=h(:) ;  um1(:)=u(:) 

        ELSE
        
          h(:)=hm2(:)+2*dt*dh(:)
          u(:)=um2(:)+2*dt*du(:)

          hm2(:)=hm1(:) ; um2(:)=um1(:) 
          hm1(:)=h(:) ; um1(:)=u(:) 
        ENDIF
          
      ENDDO
    END SUBROUTINE leapfrog_scheme  
 
    SUBROUTINE leapfrog_matsuno_scheme
    IMPLICIT NONE
      INTEGER :: ind

      DO ind=1,ndomain
        h=f_h(ind)   ; u=f_u(ind)  
        hm1=f_hm1(ind) ; um1=f_um1(ind)  
        hm2=f_hm2(ind) ; um2=f_um2(ind) 
        dh=f_dh(ind) ; du=f_du(ind) ; 

        
        IF (MOD(it,matsuno_period+1)==0) THEN
          hm1(:)=h(:) ; um1(:)=u(:) 
          h(:)=hm1(:)+dt*dh(:)
          u(:)=um1(:)+dt*du(:)

        ELSE IF (MOD(it,matsuno_period+1)==1) THEN

          h(:)=hm1(:)+dt*dh(:)
          u(:)=um1(:)+dt*du(:)

          hm2(:)=hm1(:) ;  um2(:)=um1(:) 
          hm1(:)=h(:) ;  um1(:)=u(:) 

        ELSE 

          h(:)=hm2(:)+2*dt*dh(:)
          u(:)=um2(:)+2*dt*du(:)

          hm2(:)=hm1(:) ; um2(:)=um1(:) 
          hm1(:)=h(:) ;  um1(:)=u(:) 

        ENDIF
      
      ENDDO
      
    END SUBROUTINE leapfrog_matsuno_scheme  
         
    SUBROUTINE adam_bashforth_scheme
    IMPLICIT NONE
      INTEGER :: ind

      DO ind=1,ndomain
        h=f_h(ind)   ; u=f_u(ind)   
        dh=f_dh(ind) ; du=f_du(ind) 
        dhm1=f_dhm1(ind) ; dum1=f_dum1(ind) 
        dhm2=f_dhm2(ind) ; dum2=f_dum2(ind) 

        IF (it==0) THEN
          dhm1(:)=dh(:) ; dum1(:)=du(:) 
          dhm2(:)=dhm1(:) ; dum2(:)=dum1(:)
        ENDIF
              
        h(:)=h(:)+dt*(23*dh(:)-16*dhm1(:)+5*dhm2(:))/12
        u(:)=u(:)+dt*(23*du(:)-16*dum1(:)+5*dum2(:))/12

        dhm2(:)=dhm1(:) ; dum2(:)=dum1(:) 
        dhm1(:)=dh(:) ; dum1(:)=du(:) 

      ENDDO
     
      
    END SUBROUTINE adam_bashforth_scheme

  END SUBROUTINE timeloop
    
    
  
END MODULE timeloop_sw_mod