module tendency_module

use fc_params, only : phys_kind
use domain_params, only : my_task

implicit none
private

   include "mpif.h"
   integer mpierr

integer n, np, npp, ntem   ! time level indices
integer mstep
integer::  flag_counter = 0

integer, public, parameter :: max_tend = 20
type tendency
   real(kind=phys_kind), allocatable, dimension(:,:,:,:) :: tt
end type tendency

type (tendency), dimension(max_tend), target :: tend3

public :: tendency_init, tendency_update

contains 

    subroutine tendency_init(t,tp,tpp,tflag, idim, jdim, kdim)
    integer, intent(in) :: idim, jdim, kdim
    integer, intent(out) :: tflag
    real(kind=phys_kind), intent(out), dimension(:,:,:), pointer :: &
        t, tp, tpp

       if(flag_counter == max_tend) then
         if(my_task==0) print *, 'too many tendencies, increase max_tend'
         call mpi_abort(MPI_COMM_WORLD, mpierr)
       endif
    
       flag_counter = flag_counter + 1
       tflag = flag_counter
       allocate(tend3(tflag)%tt(idim, jdim, kdim, 3))
       n = 1; np = 2; npp = 3
       tend3(tflag)%tt(idim, jdim, kdim, 3) = 0
       t => tend3(tflag)%tt(:, :, :, n)
       tp => tend3(tflag)%tt(:, :, :, np)
       tpp => tend3(tflag)%tt(:, :, :, npp)
       
    end subroutine tendency_init

    subroutine tendency_update(t,tp,tpp,tflag,istep)
    
    integer, intent(in) :: tflag, istep
    real(kind=phys_kind), intent(inout), dimension(:,:,:), pointer :: &
        t, tp, tpp

    if(.not.allocated(tend3(tflag)%tt)) then
         if(my_task==0) print *, ' tendency_update bad flag passed'
         call mpi_abort(MPI_COMM_WORLD, mpierr)
    endif
    
! nullify the pointers   
       nullify(t) 
       nullify(tp) 
       nullify(tpp) 

! update the timelevel indices
    if(istep .ne. mstep) then
       ntem = npp
       npp = np
       np = n
       n = ntem
       mstep = istep
    endif
      
! reset the pointers   
       t => tend3(tflag)%tt(:, :, :, n)
       tp => tend3(tflag)%tt(:, :, :, np)
       tpp => tend3(tflag)%tt(:, :, :, npp)
    
    end subroutine tendency_update

end module tendency_module


! below is a skeleton main program to illustrate the module above
program fc_model

  use fc_params, only : phys_kind
  use microphysics_module, only : micro_scheme_name
  use tendency_module, only : tendency_init, tendency_update

  implicit none

    integer itheta; real(kind=phys_kind), pointer, dimension(:,:,:) :: tendency_theta, tendency_theta_p, tendency_theta_pp
    integer iqvl  ; real(kind=phys_kind), pointer, dimension(:,:,:) :: tendency_qvl, tendency_qvl_p, tendency_qvl_pp
    integer iqcl  ; real(kind=phys_kind), pointer, dimension(:,:,:) :: tendency_qcl, tendency_qcl_p, tendency_qcl_pp
    integer inqcl  ; real(kind=phys_kind), pointer, dimension(:,:,:) :: tendency_nqcl, tendency_nqcl_p, tendency_nqcl_pp

!...
        call tendency_init(tendency_theta,tendency_theta_p,tendency_theta_pp,itheta,imng2,jmng2,km)
        call tendency_init(tendency_qvl,tendency_qvl_p,tendency_qvl_pp,iqvl,imng2,jmng2,km)
        call tendency_init(tendency_qcl,tendency_qcl_p,tendency_qcl_pp,iqcl,imng2,jmng2,km)
     if(trim(micro_scheme_name()) == 'm2005') then
        call tendency_init(tendency_nqcl,tendency_nqcl_p,tendency_nqcl_pp,inqcl,imng2,jmng2,km)
     endif
!...
! the following is in an integrated time loop
        call tendency_update(tendency_theta,tendency_theta_p,tendency_theta_pp,itheta,istep3)
        call tendency_update(tendency_qvl,tendency_qvl_p,tendency_qvl_pp,iqvl,istep3)
        call tendency_update(tendency_qcl,tendency_qcl_p,tendency_qcl_pp,iqcl,istep3)
     if(trim(micro_scheme_name()) == 'm2005') then
        call tendency_update(tendency_nqcl,tendency_nqcl_p,tendency_nqcl_pp,inqcl,istep3)
     endif

end program fc_model