module condensation_generic_mod
    implicit none
    
contains
    
    subroutine condensation_generic(ngrid,nlayer,nq,ptimestep, pplev, pplay,    &
                pt, pq, pdt, pdq, pdtlsc, pdqvaplsc, pdqliqlsc)
        use ioipsl_getin_p_mod, only: getin_p !-> to get the metallicity 
        use generic_cloud_common_h, only : RLVTT, cpp, &
        Psat_generic,Lcpdqsat_generic,specie_parameters,specie_parameters_table
        USE tracer_h
        IMPLICIT none

!=======================================================================
!     
!     Purpose
!     -------
!     Calculates large-scale condensation of generic tracer "tname".
!     By convention, tname ends with the suffix "_vap", as it represents the 
!     gas phase of the generic tracer
!     
!     Authors
!     -------
!     Adapted from largescale.F90 by Lucas Teinturier (2022)
!     largescale.F90 adapted from the LMDTERRE code by R. Wordsworth (2009)
!     Original author Z. X. Li (1993)
!     
!=========================================================================

        INTEGER, intent(in) :: ngrid,nlayer,nq

!       Arguments
        REAL, intent(in) :: ptimestep             ! intervalle du temps (s)
        REAL, intent(in) :: pplev(ngrid,nlayer+1) ! pression a inter-couche
        REAL, intent(in) :: pplay(ngrid,nlayer)   ! pression au milieu de couche
        REAL, intent(in) :: pt(ngrid,nlayer)      ! temperature (K)
        REAL, intent(in) :: pq(ngrid,nlayer,nq)   ! tracer mixing ratio (kg/kg)
        REAL, intent(in) :: pdt(ngrid,nlayer)     ! physical temperature tendency (K/s)
        REAL, intent(in) :: pdq(ngrid,nlayer,nq)  ! physical tracer tendency (K/s)
        ! CHARACTER(*), intent(in) :: tname_vap     ! name of the tracer we consider. BY convention, it ends with _vap !!!
        REAL, intent(out) :: pdtlsc(ngrid,nlayer)  ! incrementation de la temperature (K)
        REAL, intent(out) :: pdqvaplsc(ngrid,nlayer,nq) ! incrementation de la vapeur du traceur
        REAL, intent(out) :: pdqliqlsc(ngrid,nlayer,nq) ! incrementation du traceur liquide

!       Options : 
        real, save :: metallicity !metallicity of planet
!$OMP THREADPRIVATE(metallicity)

!       Local variables
        INTEGER i, k , nn, iq
        INTEGER,PARAMETER :: nitermax=5000
        DOUBLE PRECISION,PARAMETER :: alpha=.1,qthreshold=1.d-8
        ! JL13: if "careful, T<Tmin in psat water" appears often, you may want to stabilise the model by
        !                   decreasing alpha and increasing nitermax accordingly
        DOUBLE PRECISION zq(ngrid)
        DOUBLE PRECISION zcond(ngrid),zcond_iter
        DOUBLE PRECISION zqs(ngrid)
        real zt(ngrid),local_p,psat_tmp,dlnpsat_tmp,Lcp,zqs_temp,zdqs
        integer igcm_generic_vap, igcm_generic_ice! index of the vap and ice of generic_tracer
        ! CHARACTER(len=*) :: tname_ice
        ! evaporation calculations
        REAL dqevap(ngrid,nlayer),dtevap(ngrid,nlayer)     
        REAL qevap(ngrid,nlayer,nq)
        REAL tevap(ngrid,nlayer)

        DOUBLE PRECISION zx_q(ngrid)
        LOGICAL,SAVE :: firstcall=.true.
!$OMP THREADPRIVATE(firstcall)
        IF (firstcall) THEN
                write(*,*) "value for metallicity? "
                metallicity=0.0 ! default value
                call getin_p("metallicity",metallicity)
                write(*,*) " metallicity = ",metallicity
                firstcall = .false.
        ENDIF
!       Initialisation of outputs and local variables
        pdtlsc(1:ngrid,1:nlayer)  = 0.0
        pdqvaplsc(1:ngrid,1:nlayer,1:nq)  = 0.0
        pdqliqlsc(1:ngrid,1:nlayer,1:nq) = 0.0
        dqevap(1:ngrid,1:nlayer)=0.0
        dtevap(1:ngrid,1:nlayer)=0.0
        qevap(1:ngrid,1:nlayer,1:nq)=0.0
        tevap(1:ngrid,1:nlayer)=0.0
        igcm_generic_vap = -1
        igcm_generic_ice = -1
        ! Let's loop on tracers 
        do iq=1,nq
                if((is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0)) then
                        write(*,*) "There is a specie which is condensable, for generic condensation : ", noms(iq)

!                       Let's get the index of our tracers (we look for igcm _generic_vap and igcm_generic_ice)
                        ! tname_ice = trim(noms(iq)(1:len(tname_ice)-3))//"ice" 
                        ! print*,trim(adjustl(trim(noms(iq))(9:len(trim(noms(iq)))-4))) !testing here, should go away
                        ! stop

                        igcm_generic_vap=iq

                        ! We look for the corresponding ice traceur (before or after in the list of traceurs, maybe could be generalised to the whole list)
                        if ((noms(iq)(9:len(trim(noms(iq)))-4) .eq. noms(iq+1)(9:len(trim(noms(iq+1)))-4)) .and. (index(noms(iq+1),"ice") .ne. 0)) then
                                igcm_generic_ice = iq+1
                        else if ((iq .gt. 1) .and. (noms(iq)(9:len(trim(noms(iq)))-4) .eq. noms(iq-1)(9:len(trim(noms(iq-1)))-4)) .and. (index(noms(iq-1),"ice") .ne. 0)) then
                                igcm_generic_ice = iq-1
                        else
                                write(*,*) "ERROR : You set a vap traceur but you forgot to set the corresponding ice traceur"
                        endif

                        ! Need to call specie_parameters of the considered specie 
                        call specie_parameters(noms(iq)(9:len(trim(noms(iq)))-4))
                        !call specie_parameters_table(noms(iq)(9:len(trim(noms(iq)))-4))
                        
                        Lcp=RLVTT/cpp ! need to be init here

                        !  Vertical loop (from top to bottom)
                        DO k = nlayer, 1, -1
                                zt(1:ngrid)=pt(1:ngrid,k)+pdt(1:ngrid,k)*ptimestep

                                ! Computes Psat and the partial condensation
                                DO i = 1, ngrid
                                        local_p=pplay(i,k)
                                        if(zt(i).le.15.) then
                                                print*,'in lsc',i,k,zt(i)
                                        !	    zt(i)=15.   ! check too low temperatures
                                        endif
                                        zx_q(i) = pq(i,k,igcm_generic_vap)+pdq(i,k,igcm_generic_vap)*ptimestep
                                        ! iterative process to stabilize the scheme when large water amounts JL12
                                        zcond(i) = 0.0d0
                                        Do nn=1,nitermax  
                                                call Psat_generic(zt(i),local_p,metallicity,psat_tmp,zqs_temp)
                                                zqs(i)=zqs_temp
                                                call Lcpdqsat_generic(zt(i),local_p,psat_tmp,zqs_temp,zdqs,dlnpsat_tmp)
                                                zcond_iter = alpha*(zx_q(i)-zqs(i))/(1.d0+zdqs)
                                                !zcond can be negative here
                                                zx_q(i) = zx_q(i) - zcond_iter
                                                zcond(i) = zcond(i) + zcond_iter
                                                zt(i) = zt(i) + zcond_iter*Lcp
                                                if (ABS(zcond_iter/alpha/zqs(i)).lt.qthreshold) exit
                                                if (nn.eq.nitermax) print*,'itermax in largescale'
                                        End do ! niter
                                        zcond(i)=MAX(zcond(i),-(pq(i,k,igcm_generic_ice)+pdq(i,k,igcm_generic_ice)*ptimestep))
                                        zcond(i) = zcond(i)/ptimestep
                                ENDDO ! i=1,ngrid

                                !Tendances de t et q
                                pdqvaplsc(1:ngrid,k,igcm_generic_vap)  = - zcond(1:ngrid) 
                                pdqliqlsc(1:ngrid,k,igcm_generic_ice) = - pdqvaplsc(1:ngrid,k,igcm_generic_vap)
                                pdtlsc(1:ngrid,k)  = pdtlsc(1:ngrid,k) + pdqliqlsc(1:ngrid,k,igcm_generic_ice)*Lcp

                        Enddo ! k= nlayer, 1, -1
                endif !(is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0))
        enddo ! iq=1,nq

    end subroutine condensation_generic
end module condensation_generic_mod