MODULE timeloop_gcm_mod
  USE genmod
  USE transfert_mod
  USE etat0_mod
  
CONTAINS
  
  SUBROUTINE timeloop
  USE field_mod
  USE domain_mod
  USE write_field
  USE dimensions
  USE geometry
  USE transfert_mod
  USE metric
  USE dissip_gcm_mod
  USE ioipsl
  USE caldyn_mod
  USE theta2theta_rhodz_mod
  USE etat0_mod
  USE guided_mod
  USE advect_tracer_mod
  
  IMPLICIT NONE
  TYPE(t_field),POINTER :: f_phis(:)
  TYPE(t_field),POINTER :: f_theta(:)
  TYPE(t_field),POINTER :: f_q(:)
  TYPE(t_field),POINTER :: f_dtheta(:)
  TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:)
  TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:)
  TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:)
  TYPE(t_field),POINTER :: f_dps(:),f_dpsm1(:), f_dpsm2(:)
  TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:)
  TYPE(t_field),POINTER :: f_dtheta_rhodz(:),f_dtheta_rhodzm1(:),f_dtheta_rhodzm2(:)

  REAL(rstd),POINTER :: phis(:)
  REAL(rstd),POINTER :: q(:,:,:)
  REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:)
  REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:)
  REAL(rstd),POINTER :: theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:)
  REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:)
  REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:)
  REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
  REAL(rstd) :: dt
  INTEGER :: ind
  INTEGER :: it,i,j,n
  CHARACTER(len=255) :: scheme
  INTEGER :: matsuno_period
  INTEGER :: itaumax
  INTEGER :: write_period  
  INTEGER :: itau_out

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

  itaumax=100
  CALL getin('itaumax',itaumax)

  write_period=0
  CALL getin('write_period',write_period)
  itau_out=INT(write_period/dt)
  
  scheme='adam_bashforth'
  CALL getin('scheme',scheme)
  
  matsuno_period=5
  CALL getin('matsuno_period',matsuno_period)
  IF (TRIM(scheme)=='leapfrog') matsuno_period=itaumax+1

  CALL allocate_field(f_phis,field_t,type_real)
  
  CALL allocate_field(f_ps,field_t,type_real)
  CALL allocate_field(f_psm1,field_t,type_real)
  CALL allocate_field(f_psm2,field_t,type_real)
  CALL allocate_field(f_dps,field_t,type_real)
  CALL allocate_field(f_dpsm1,field_t,type_real)
  CALL allocate_field(f_dpsm2,field_t,type_real)

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

  CALL allocate_field(f_theta,field_t,type_real,llm)
  CALL allocate_field(f_dtheta,field_t,type_real,llm)

  CALL allocate_field(f_q,field_t,type_real,llm,nqtot)

  CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
  CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm)
  CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
  CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm)
  CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm)
  CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm)

  CALL init_dissip(dt)
  CALL init_caldyn(dt)
  CALL init_guided(dt)
  CALL init_advect_tracer(dt)
  
  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
  
  DO it=0,itaumax
    PRINT *,"It No :",It

    CALL guided(it,f_ps,f_theta_rhodz,f_u,f_q)
    CALL caldyn(it,f_phis,f_ps,f_theta_rhodz,f_u, f_dps, f_dtheta_rhodz, f_du)
    CALL advect_tracer(f_ps,f_u,f_q)
    
    SELECT CASE (TRIM(scheme))
      CASE('euler')
        CALL  euler_scheme

      CASE ('leapfrog')
        CALL leapfrog_scheme

      CASE ('leapfrog_matsuno')
        CALL  leapfrog_matsuno_scheme

      CASE ('adam_bashforth')
        CALL dissip(f_u,f_du,f_ps,f_theta_rhodz,f_dtheta_rhodz)
        CALL adam_bashforth_scheme

      CASE default
        PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),"> options are <euler>, <leapfrog>, <leapfrog_matsuno>, <adam_bashforth>" 
        STOP
        
    END SELECT

    
    IF ( itau_out>0 .AND. MOD(it,itau_out)==0) THEN
      CALL writefield("q",f_q)
      CALL writefield("ps",f_ps)
    ENDIF

  ENDDO
  
  CONTAINS

    SUBROUTINE Euler_scheme
    IMPLICIT NONE
      INTEGER :: ind
      
      DO ind=1,ndomain
        ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
        ps(:)=ps(:)+dt*dps(:)
        u(:,:)=u(:,:)+dt*du(:,:)
        theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)
      ENDDO
      
    END SUBROUTINE Euler_scheme
    
    SUBROUTINE leapfrog_scheme
    IMPLICIT NONE
      INTEGER :: ind

      DO ind=1,ndomain
        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)

        IF (it==0) THEN
          psm2(:)=ps(:) ; theta_rhodzm2(:,:)=theta_rhodz(:,:) ; um2(:,:)=u(:,:) 

          ps(:)=ps(:)+dt*dps(:)
          u(:,:)=u(:,:)+dt*du(:,:)
          theta_rhodz(:,:)=theta_rhodz(:,:)+dt*dtheta_rhodz(:,:)

          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
        ELSE
        
          ps(:)=psm2(:)+2*dt*dps(:)
          u(:,:)=um2(:,:)+2*dt*du(:,:)
          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*dt*dtheta_rhodz(:,:)

          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 
        ENDIF
          
      ENDDO
    END SUBROUTINE leapfrog_scheme  
 
    SUBROUTINE leapfrog_matsuno_scheme
    IMPLICIT NONE
      INTEGER :: ind

      DO ind=1,ndomain
        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
        psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)

        
        IF (MOD(it,matsuno_period+1)==0) THEN
          psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
          ps(:)=psm1(:)+dt*dps(:)
          u(:,:)=um1(:,:)+dt*du(:,:)
          theta_rhodz(:,:)=theta_rhodzm1(:,:)+dt*dtheta_rhodz(:,:)

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

          ps(:)=psm1(:)+dt*dps(:)
          u(:,:)=um1(:,:)+dt*du(:,:)
          theta_rhodz(:,:)=theta_rhodzm1(:,:)+dt*dtheta_rhodz(:,:)

          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 

        ELSE 

          ps(:)=psm2(:)+2*dt*dps(:)
          u(:,:)=um2(:,:)+2*dt*du(:,:)
          theta_rhodz(:,:)=theta_rhodzm2(:,:)+2*dt*dtheta_rhodz(:,:)

          psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:) 
          psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:) 

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

      DO ind=1,ndomain
        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
        dpsm1=f_dpsm1(ind) ; dum1=f_dum1(ind) ; dtheta_rhodzm1=f_dtheta_rhodzm1(ind)
        dpsm2=f_dpsm2(ind) ; dum2=f_dum2(ind) ; dtheta_rhodzm2=f_dtheta_rhodzm2(ind)

        IF (it==0) THEN
          dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)
          dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
        ENDIF
              
        ps(:)=ps(:)+dt*(23*dps(:)-16*dpsm1(:)+5*dpsm2(:))/12
        u(:,:)=u(:,:)+dt*(23*du(:,:)-16*dum1(:,:)+5*dum2(:,:))/12
        theta_rhodz(:,:)=theta_rhodz(:,:)+dt*(23*dtheta_rhodz(:,:)-16*dtheta_rhodzm1(:,:)+5*dtheta_rhodzm2(:,:))/12

        dpsm2(:)=dpsm1(:) ; dum2(:,:)=dum1(:,:) ; dtheta_rhodzm2(:,:)=dtheta_rhodzm1(:,:)
        dpsm1(:)=dps(:) ; dum1(:,:)=du(:,:) ; dtheta_rhodzm1(:,:)=dtheta_rhodz(:,:)

      ENDDO
     
      
    END SUBROUTINE adam_bashforth_scheme

  END SUBROUTINE timeloop
    
    
  
END MODULE timeloop_gcm_mod