MODULE compute_dtau_mod

	IMPLICIT NONE

	REAL,PARAMETER :: ti_injection=10/24. ! time of beginning injection
	REAL,PARAMETER :: tf_injection=18/24. ! time of end injection
        REAL,PARAMETER :: t_scenario=14/24.   ! time of tauref_scenario

	CONTAINS

	SUBROUTINE compute_dtau(ngrid,nlayer,            		&
       & 			zday,pplev,tauref,			&
       & 			ptimestep,dustliftday,local_time)

	USE geometry_mod, only: longitude_deg
        USE time_phylmdz_mod, only: dtphys, daysec
        USE comcstfi_h, only: g
        USE tracer_mod, only: alpha_lift,igcm_dust_mass,igcm_dust_number
	
	IMPLICIT NONE
        
        include "callkeys.h"
	
        INTEGER, INTENT(in) :: ngrid
	INTEGER, INTENT(in) :: nlayer
        REAL, INTENT(in) :: zday ! date at lon=0, in fraction of sols
	REAL, INTENT(in) :: pplev(ngrid,nlayer+1) ! pressure (Pa)
	REAL, INTENT(in) :: tauref(ngrid) ! Computed dust opacity at 610Pa
	REAL, INTENT(in) :: ptimestep 
        REAL, INTENT(in) :: local_time(ngrid)
        REAL, INTENT(out) :: dustliftday(ngrid) ! Dust injection rate (s-1)
        
        INTEGER :: ig, l
	INTEGER, SAVE :: nb_daystep ! nomber of step a day
        REAL :: tauref_scenario(ngrid) ! from dust scenario
        REAL :: zday_scenario
        REAL,ALLOCATABLE,SAVE :: local_time_prev(:)
        REAL,ALLOCATABLE,SAVE :: dtau(:) ! Dust opacity increment (at 610Pa)
        REAL,PARAMETER :: odpref=610. !DOD reference pressure (Pa)
        
        LOGICAL, SAVE :: firstcall=.TRUE. ! signals first call to physics
        
        
        IF(firstcall)THEN
                ALLOCATE(local_time_prev(ngrid))
                ALLOCATE(dtau(ngrid))
                DO ig=1,ngrid
                   local_time_prev(ig)=modulo(1.+(zday-ptimestep/daysec)&
                                      -INT(zday-ptimestep/daysec) 	&
                                      +(longitude_deg(ig)/15)/24,1.)
                   dtau(ig)=0.
                ENDDO
                nb_daystep=(daysec/dtphys)
                firstcall=.FALSE.
        ENDIF
        
	! 1. Obtain tauref_scenario from dust scenario
	zday_scenario=zday-modulo(zday,1.)
        call read_dust_scenario(ngrid,nlayer,zday_scenario,pplev,	&
       & 			tauref_scenario)
       ! for diagnostics
        call WRITEDIAGFI(ngrid,"tauref_scenario","tauref_scenario",     &
       &                  "",2,tauref_scenario)

	! 2. Compute dtau
	DO ig=1,ngrid
         IF ((local_time(ig).ge.t_scenario).and. 			&
       &          (local_time_prev(ig).lt.(t_scenario)))THEN
	         dtau(ig)=tauref_scenario(ig)-tauref(ig)
         ENDIF
         IF (dtau(ig).LT.0) THEN
                  dtau(ig)=0.
         ENDIF
         
        ! 3. Use dtau to compute dustliftday
         dustliftday(ig)=coeff_injection*				&
       &		(dtau(ig)*pplev(ig,1)/odpref)			&
       &		/(daysec*(tf_injection-ti_injection))
        ENDDO 
         
        ! 4. Save local time 
        local_time_prev(1:ngrid)=local_time(1:ngrid)
        
        end subroutine compute_dtau

	end module compute_dtau_mod