module condensation_generic_mod
    implicit none
    
contains
    
    subroutine condensation_generic(ngrid,nlayer,nq,ptimestep, pplev, pplay,    &
                pt, pq, pdt, pdq, pdtlsc, pdqvaplsc, pdqliqlsc, rneb)
        use ioipsl_getin_p_mod, only: getin_p !-> to get the metallicity 
        use generic_cloud_common_h
        USE tracer_h
        USE mod_phys_lmdz_para, only: is_master
        use generic_tracer_index_mod, only: generic_tracer_index
        IMPLICIT none

!=======================================================================
!     
!     Purpose
!     -------
!     Calculates large-scale condensation of generic tracer "tname".
!     By convention, tname ends with the suffix "_gas", as it represents the 
!     gas phase of the generic tracer
!     
!     Authors
!     -------
!     Adapted from largescale.F90 by Lucas Teinturier & Noé Clément (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_gas     ! name of the tracer we consider. BY convention, it ends with _gas !!!
        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
        REAL, intent(out) :: rneb(ngrid,nlayer,nq)    ! fraction nuageuse

!       Options : 
        real, save :: metallicity !metallicity of planet
        REAL, SAVE :: qvap_deep   ! deep mixing ratio of vapor when simulating bottom less planets
        REAL, SAVE :: qvap_top   ! top mixing ratio of vapor when simulating bottom less planets
        logical, save :: perfect_gas_profile
!$OMP THREADPRIVATE(metallicity, qvap_deep, qvap_top, perfect_gas_profile)

!       Local variables

        ! to call only one time the ice/vap pair of a tracer
        logical call_ice_gas_generic

        INTEGER i, k , nn, iq
        INTEGER,PARAMETER :: nitermax=5000
        REAL tau ! tau is in seconds and must not be lower than the physical step time.
        integer k_cold_trap
        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
        real zqs_temp_1, zqs_temp_2, zqs_temp_3
        integer igcm_generic_gas, 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)
        DOUBLE PRECISION zy_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

                write(*,*) "Deep generic tracer vapor mixing ratio ? (no effect if negative) "
                qvap_deep=-1. ! default value
                call getin_p("qvap_deep",qvap_deep)
                write(*,*) " qvap_deep = ",qvap_deep   

                write(*,*) "top generic tracer vapor mixing ratio ? (no effect if negative) "
                qvap_top=-1. ! default value
                call getin_p("qvap_top",qvap_top)
                write(*,*) " qvap_top = ",qvap_top 
                
                write(*,*) " perfect_gas_profile ? "
                perfect_gas_profile=.false. ! default value
                call getin_p("perfect_gas_profile",perfect_gas_profile)
                write(*,*) " perfect_gas_profile = ",perfect_gas_profile

                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
        rneb(1:ngrid,1:nlayer,1:nq) = 0.0
        ! Let's loop on tracers 
        do iq=1,nq

                call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic)

                if(call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer
                        m=constants_mass(iq)
                        delta_gasH=constants_delta_gasH(iq)
                        Tref=constants_Tref(iq)
                        Pref=constants_Pref(iq)
                        epsi_generic=constants_epsi_generic(iq)
                        RLVTT_generic=constants_RLVTT_generic(iq)
                        metallicity_coeff=constants_metallicity_coeff(iq)

                        Lcp=RLVTT_generic/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_gas)+pdq(i,k,igcm_generic_gas)*ptimestep

                                        call Psat_generic(zt(i),local_p,metallicity,psat_tmp,zqs_temp)
                                        zy_q(i) = pq(i,k,igcm_generic_ice)+pdq(i,k,igcm_generic_ice)*ptimestep
                                        
                                        if ((zx_q(i) .le. zqs_temp) .and. (zy_q(i) .eq. 0.)) then
                                                ! if we are are not saturated and if there is no ice
                                                ! then no change

                                                zcond(i) = 0.0d0
                                        
                                        else    ! if we are saturated : we must evaporate
                                                ! if there is ice : we must check if we can condensate

                                                ! 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

                                                ! if zcond(i) > 0, zcond(i) is the amount of vapor that we can condensate
                                                !                       because we are saturated. zcond(i) will not change below
                                                ! if zcond(i) < 0, zcond(i) is the amount of ice that we can evaporate.
                                                !                       We can not evaporate more than the existing ice, 
                                                !                       so the line below is to check how much we can evaporate.
                                                !                       If there is no ice available, zcond(i) will be 0. (first condidition of "if")

                                                zcond(i)=MAX(zcond(i),-(pq(i,k,igcm_generic_ice)+pdq(i,k,igcm_generic_ice)*ptimestep))
                                        
                                        endif

                                        if (zcond(i) .gt. 0.) then
                                                rneb(i,k,iq)=1
                                        else
                                                rneb(i,k,iq)=0.
                                        endif

                                        zcond(i) = zcond(i)/ptimestep
                                ENDDO ! i=1,ngrid

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

                        Enddo ! k= nlayer, 1, -1

                        if ((perfect_gas_profile) .and. (ngrid.eq.1)) then
                                ! perfect_gas_profile is a mode that should a priori only be used in 1D:
                                ! as it is written below, it aims to force the vap profile:
                                ! - below condensation, it is forced to qvap_deep
                                ! - at condensation levels, it is forced to 99% of sat
                                ! - above the cold trap, it is forced to the value allowed by the cold trap

                                ! perfect_gas_profile can be customed as you want
                                
                                tau = 10. * ptimestep ! tau must not be lower than the physical step time. 

                                k_cold_trap = 2
                                do k=2,nlayer-1

                                        zt(1)=pt(1,k-1)+pdt(1,k-1)*ptimestep
                                        call Psat_generic(zt(1),pplay(1,k-1),metallicity,psat_tmp,zqs_temp_1)
                                        zt(1)=pt(1,k)+pdt(1,k)*ptimestep
                                        call Psat_generic(zt(1),pplay(1,k),metallicity,psat_tmp,zqs_temp_2)
                                        zt(1)=pt(1,k+1)+pdt(1,k+1)*ptimestep
                                        call Psat_generic(zt(1),pplay(1,k+1),metallicity,psat_tmp,zqs_temp_3)

                                        if ((zqs_temp_1 .gt. zqs_temp_2) .and. (zqs_temp_3 .gt. zqs_temp_2)) then
                                                k_cold_trap = k
                                                exit
                                        endif
                                enddo
                                if (k_cold_trap .lt. nlayer) then
                                        do k=k_cold_trap+1,nlayer
                                                pdqvaplsc(1,k,igcm_generic_gas) = (pq(1,k_cold_trap,igcm_generic_gas) - pq(1,k,igcm_generic_gas))/tau - pdq(1,k,igcm_generic_gas)
                                        enddo
                                endif

                                do k=1,k_cold_trap
                                        zt(1)=pt(1,k)+pdt(1,k)*ptimestep
                                        call Psat_generic(zt(1),pplay(1,k),metallicity,psat_tmp,zqs_temp)
                                        if (zqs_temp .gt. qvap_deep) then
                                                pdqvaplsc(1,k,igcm_generic_gas)  = (qvap_deep - pq(1,k,igcm_generic_gas))/tau - pdq(1,k,igcm_generic_gas)
                                        endif
                                        if (zqs_temp .lt. qvap_deep) then
                                                pdqvaplsc(1,k,igcm_generic_gas)  = (0.99*zqs_temp - pq(1,k,igcm_generic_gas))/tau - pdq(1,k,igcm_generic_gas)
                                        endif
                                enddo

                                pdqliqlsc(1:ngrid,:,igcm_generic_ice) = 0.
                                pdtlsc(1:ngrid,:)  = 0.
                        endif

                        if (qvap_deep >= 0.) then
                                tau = 10. * ptimestep ! tau must not be lower than the physical step time. 
                                ! brings lower generic vapor ratio to a fixed value.
                                ! tau is in seconds and must not be lower than the physical step time.
                                pdqvaplsc(1:ngrid,1,igcm_generic_gas) = (qvap_deep - pq(1:ngrid,1,igcm_generic_gas))/tau - pdq(1:ngrid,1,igcm_generic_gas)
                        endif
                        if (qvap_top >= 0.) then
                                tau = 10. * ptimestep ! tau must not be lower than the physical step time. 
                                ! brings lower generic vapor ratio to a fixed value.
                                ! tau is in seconds and must not be lower than the physical step time.
                                pdqvaplsc(1:ngrid,nlayer,igcm_generic_gas) = (qvap_top - pq(1:ngrid,nlayer,igcm_generic_gas))/tau - pdq(1:ngrid,nlayer,igcm_generic_gas)
                        endif
                endif !if(call_ice_gas_generic)
        enddo ! iq=1,nq

    end subroutine condensation_generic
end module condensation_generic_mod