!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!MNH_LIC for details. version 1.
!     #######################
      MODULE MODI_PROGNOS
!     #######################
!
INTERFACE
!
SUBROUTINE PROGNOS(HLUOUT,PDT,PDZ,PLV,PCPH,PPRES,PRHOD,PRR,PTT,PRV,PRC,PS0,PCN,PCL)
!
CHARACTER(LEN=*),         INTENT(IN)    :: HLUOUT   ! Output-listing name for
                                                    ! model n
REAL,                     INTENT(IN)    :: PDT
REAL, DIMENSION(:),       INTENT(IN)    :: PPRES
REAL, DIMENSION(:),       INTENT(IN)    :: PDZ
REAL, DIMENSION(:),       INTENT(IN)    :: PLV
REAL, DIMENSION(:),       INTENT(IN)    :: PCPH
REAL, DIMENSION(:),       INTENT(IN)    :: PRHOD
REAL, DIMENSION(:),       INTENT(IN)    :: PRR
REAL, DIMENSION(:),       INTENT(INOUT) :: PTT
REAL, DIMENSION(:),       INTENT(INOUT) :: PRV
REAL, DIMENSION(:),       INTENT(INOUT) :: PRC
REAL, DIMENSION(:),       INTENT(INOUT) :: PS0
REAL, DIMENSION(:),       INTENT(INOUT) :: PCN
REAL, DIMENSION(:),       INTENT(INOUT) :: PCL
!
END SUBROUTINE PROGNOS
!
END INTERFACE
!
END MODULE MODI_PROGNOS
!
!     ###################################################################################
      SUBROUTINE PROGNOS(HLUOUT,PDT,PDZ,PLV,PCPH,PPRES,PRHOD,PRR,PTT,PRV,PRC,PS0,PCN,PCL)
!     ###################################################################################
!
!!****  * -  compute pseudo-prognostic of supersaturation according to Thouron
!                                                                     et al. 2012
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!
!!    REFERENCE
!!    ---------
!!
!!      Thouron, O., J.-L. Brenguier, and F. Burnet, Supersaturation calculation
!!      in large eddy simulation models for prediction of the droplet number
!!      concentration, Geosci. Model Dev., 5, 761-772, 2012.
!!
!!    AUTHOR
!!    ------
!!
!!      O. Thouron   * CNRM Meteo-France* : 
!!
!!    MODIFICATIONS
!!    -------------
!!     2014 G.Delautier : remplace MODD_RAIN_C2R2_PARAM par MODD_RAIN_C2R2_KHKO_PARAM
!!     2015 M.Mazoyer and O.Thouron : Physical tunings
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!
USE MODD_CST
USE MODD_RAIN_C2R2_KHKO_PARAM
USE MODD_PARAM_C2R2
USE MODI_GAMMA
USE MODE_IO_ll
!
IMPLICIT NONE
!
!*       0.1   Declarations of dummy arguments :
!
!
!
CHARACTER(LEN=*),         INTENT(IN)    :: HLUOUT   ! Output-listing name for
                                                    ! model n
REAL,                     INTENT(IN)    :: PDT
REAL, DIMENSION(:),       INTENT(IN)    :: PPRES
REAL, DIMENSION(:),       INTENT(IN)    :: PDZ
REAL, DIMENSION(:),       INTENT(IN)    :: PLV
REAL, DIMENSION(:),       INTENT(IN)    :: PCPH
REAL, DIMENSION(:),       INTENT(IN)    :: PRHOD
REAL, DIMENSION(:),       INTENT(IN)    :: PRR
REAL, DIMENSION(:),       INTENT(INOUT) :: PTT
REAL, DIMENSION(:),       INTENT(INOUT) :: PRV
REAL, DIMENSION(:),       INTENT(INOUT) :: PRC
REAL, DIMENSION(:),       INTENT(INOUT) :: PS0
REAL, DIMENSION(:),       INTENT(INOUT) :: PCN
REAL, DIMENSION(:),       INTENT(INOUT) :: PCL
!
!
!*       0.2   Declarations of local variables :
!
!
REAL, DIMENSION(SIZE(PRHOD,1))   ::  ZZW1,ZZW2,ZDZRC2,ZDZRC,ZCPH
REAL, DIMENSION(SIZE(PRHOD,1))   :: ZA1,ZA2,ZB,ZC,ZG
REAL, DIMENSION(SIZE(PRHOD,1))   :: ZLV,ZTT1,ZRT,ZTL,ZTT1_TEMP,ZTT2_TEMP
REAL, DIMENSION(SIZE(PRHOD,1))   :: ZRMOY,ZRVSAT1,ZRVSAT2
REAL, DIMENSION(SIZE(PRHOD,1))   :: ZVEC2  ! Work vectors forinterpolations
INTEGER, DIMENSION(SIZE(PRHOD,1)):: IVEC2   ! Vectors of indices for interpolations
INTEGER :: J1,J2
REAL,DIMENSION(SIZE(PS0,1))      ::MEM_PS0,ADJU2
REAL::AER_RAD
REAL, DIMENSION(SIZE(PRHOD,1))   :: ZFLAG_ACT   !Flag for activation
!
INTEGER                          :: IRESP      ! Return code of FM routines
INTEGER                          :: ILUOUT     ! Logical unit of output listing

!minimum radius of cloud droplet
AER_RAD=1.0E-6
!
ZFLAG_ACT(:)=0.0
!ACTIVATION
ZVEC2(:) =0.0
IVEC2(:) =0.0
!
DO J1 = 1,4
  WHERE (PS0(:).GT.0.0)
   ZVEC2(:) = MAX( 1.00001, MIN( FLOAT(NHYP)-0.00001,      &
                   XHYPINTP1*LOG(PS0(:))+XHYPINTP2 ) )
   IVEC2(:) = INT( ZVEC2(:) )
   ZVEC2(:) = ZVEC2(:) - FLOAT( IVEC2(:) )
  END WHERE
END DO
ZZW1(:) =0.0
WHERE (PS0(:).GT.0.0)
ZZW1(:) =   XHYPF12( IVEC2(:)+1 )* ZVEC2(:)      &
             - XHYPF12( IVEC2(:)   )*(ZVEC2(:) - 1.0)
END WHERE
!
! the CCN spectra formula uses ZSMAX in percent
!
ZZW2(:)=0.0
IF (XCONC_CCN > 0) THEN
WHERE (PS0(:).GT.0.0)
   ZZW2(:) = MIN( XCONC_CCN,XCHEN * (100.0*PS0(:))**XKHEN * ZZW1(:) ) 
END WHERE
ELSE
WHERE (PS0(:).GT.0.0)
   ZZW2(:) = XCHEN * (100.0*PS0(:))**XKHEN * ZZW1(:) 
END WHERE
ENDIF
ZZW2(:)=MAX( (ZZW2(:)-PCN(:)),0.0 )
!
WHERE (ZZW2(:).LT.1.0)   !Non physique d'activer moins d'une particule
ZZW2(:)=0
END WHERE
!
!
!FLAG ACTIVE A TRUE (1.0) si on active pas
DO J2=1,SIZE(PRC,1)
 IF (ZZW2(J2).EQ.0.0) THEN
 ZFLAG_ACT(J2)=1.0
 ENDIF
ENDDO
!
! Mean radius
ZRMOY(:)=0.0
DO J2=1,SIZE(PRC,1)
 IF (PRC(J2).NE.0.0) THEN
 ZRMOY(J2)=(MOMG(XALPHAC,XNUC,3.0)*4.0*XPI*PCL(J2)*XRHOLW/&
        (3.0*PRC(J2)*PRHOD(J2)))**(1.0/3.0)
 ZRMOY(J2)=(PCL(J2)*MOMG(XALPHAC,XNUC,1.0)/ZRMOY(J2))
 ENDIF
 ZRMOY(J2)=ZRMOY(J2)+(ZZW2(J2)*AER_RAD)
ENDDO
!
PCL(:)= PCL(:)+ZZW2(:)
PCN(:)= PCN(:)+ZZW2(:)
!
!CALCUL DE A1 => Estimation de (drs/dt)f
!T(=à determiner) avant forcage; T'(=PTT) apres forcage
!Calcul de ZTT1: calculé en inversant S0(T)jusqu'à T:
! l'erreur faite sur cette inversion est supérieur à la précision
! recherchée, on applique à rs(T') pour cxalculer le DT=T'-T qui
! correspond à la variation rs(T')-rs(T). Permet de recuperer une valeur
! correcte de DT et donc de determiner T comme T=T'-DT         
!ZRVSAT1=rs(T)
!
ZRVSAT1(:)=PRV(:)/(PS0(:)+1.0)
!ZTT1<--es(T) de rs(T)
ZTT1_TEMP(:)=PPRES(:)*((((XMV / XMD)/ZRVSAT1(:))+1.0)**(-1D0))
!ZTT1<--T de es(T)
ZTT1_TEMP(:)=LOG(ZTT1_TEMP(:)/610.8)
ZTT1_TEMP(:)=(31.25*ZTT1_TEMP(:) -17.5688*273.15)/(ZTT1_TEMP(:) - 17.5688)
!es(T')
ZZW1(:)=EXP(XALPW-XBETAW/PTT(:)-XGAMW*LOG(PTT(:)))
!ZRVSAT2=rs(T')
ZRVSAT2(:)=(XMV / XMD)*ZZW1(:)/(PPRES(:)-ZZW1(:))
!ZTT2<--es(T') de rs(T')
ZTT2_TEMP(:)=PPRES(:)*((((XMV / XMD)/ZRVSAT2(:))+1.0)**(-1D0))
!ZTT2<--T' de es(T')
IF (MINVAL(ZTT2_TEMP).LT.0.0) THEN
  PRINT*,'ZTT2_TEMP',MINVAL(ZTT2_TEMP),MINLOC(ZTT2_TEMP)
  CALL CLOSE_ll(HLUOUT,IOSTAT=IRESP)
  CALL ABORT
  STOP
ENDIF
!
ZTT2_TEMP(:)=LOG(ZZW1(:)/610.8)
ZTT2_TEMP(:)=(31.25*ZTT2_TEMP(:) -17.5688*273.15)/(ZTT2_TEMP(:) - 17.5688)
!ZTT1=T'-DT
ZTT1(:)=PTT(:)-(ZTT2_TEMP(:)-ZTT1_TEMP(:))
!Lv(T)
ZLV(:) = XLVTT+(XCPV-XCL)*(ZTT1(:)-XTT)                
!
ZA1(:)=-(((PS0(:)+1.0)**2.0)/PRV(:))*(ZRVSAT2(:)-(PRV(:)/(PS0(:)+1.0)))/PDT
!G
ZG(:)= 1.0/(XRHOLW*((XRV*ZTT1(:)/(XDIVA*EXP(XALPW-(XBETAW/ZTT1(:))-(XGAMW*LOG(ZTT1(:))))))     &
+((ZLV(:)/(XTHCO*ZTT1(:)))*((ZLV(:)/(ZTT1(:)*XRV))-1.0))))
!
ZC(:)=4.0*XPI*(XRHOLW/PRHOD(:))*ZG(:)
ZDZRC(:)=0.0
ZDZRC(:)=ZC(:)*PS0(:)*ZRMOY(:)
MEM_PS0(:)=PS0(:)
!CALCUL DE B => Estimation de (drs/dT)ce
!T(=PTT) avant condensation; T'(=à determiner) apres condensation
!Lv(T),Cph(T)
ZLV(:) = XLVTT+(XCPV-XCL)*(PTT(:)-XTT)                
ZCPH(:)= XCPD+XCPV*PRV(:)+XCL*(PRC(:)+PRR(:))
!T'=T+(DT)ce
ZTT1(:)=PTT(:)+(ZDZRC(:)*PDT*ZLV(:)/ZCPH(:))
!es(T')
ZZW1(:)=EXP(XALPW-XBETAW/PTT(:)-XGAMW*LOG(PTT(:)))
!rs(T')
ZZW1(:)=(XMV / XMD)*ZZW1(:)/(PPRES(:)-ZZW1(:))
!es(Tcond)
ZZW2(:)=EXP(XALPW-XBETAW/ZTT1(:)-XGAMW*LOG(ZTT1(:)))
!rs(Tcond)
ZZW2(:)=(XMV / XMD)*ZZW2(:)/(PPRES(:)-ZZW2(:))
!
WHERE (ZTT1(:).NE.PTT(:))
 ZB(:)=(ZLV(:)/ZCPH(:))*((ZZW2(:)-ZZW1(:))/(ZTT1(:)-PTT(:)))
ELSEWHERE
 ZB(:)=0.0
 ZDZRC(:)=0.0
ENDWHERE
!Calcul de S+dS
PS0(:)=PS0(:)+((ZA1(:)-(((ZB(:)*(PS0(:)+1.0)+1.0)*ZDZRC(:))/ZRVSAT1(:)))*PDT)
!
!Ajustement tel que rv=(s+1)*rvs
ZTL(:)=PTT(:)-(PLV(:)/PCPH(:))*PRC(:)
ZRT(:)=PRC(:)+PRV(:)
ZDZRC2(:)=PRC(:)
DO J2=1,SIZE(ZDZRC,1)
  IF ((ZDZRC(J2).NE.0.0).OR.(ZDZRC2(J2).NE.0.0)) THEN 
    DO J1=1,5
     ZLV(J2) = XLVTT+(XCPV-XCL)*(PTT(J2)-XTT)                
     ZCPH(J2)=XCPD+XCPV*PRV(J2)+XCL*(PRC(J2)+PRR(J2))
     ZZW1(J2)=EXP(XALPW-XBETAW/PTT(J2)-XGAMW*LOG(PTT(J2)))
     ZRVSAT1(J2)=(XMV / XMD)*ZZW1(J2)/(PPRES(J2)-ZZW1(J2))
     PRV(J2)=MIN(ZRT(J2),(PS0(J2)+1.0)*ZRVSAT1(J2))
     PRC(J2)=MAX(ZRT(J2)-PRV(J2),0.0)
     PTT(J2)=0.5*PTT(J2)+0.5*(ZTL(J2)+(ZLV(J2)*PRC(J2)/ZCPH(J2)))
    ENDDO
    ZLV(J2) = XLVTT+(XCPV-XCL)*(PTT(J2)-XTT)                
    ZCPH(J2)=XCPD+XCPV*PRV(J2)+XCL*(PRC(J2)+PRR(J2))
    PTT(J2)=ZTL(J2)+(ZLV(J2)*PRC(J2)/ZCPH(J2))
 ENDIF
ENDDO
ADJU2(:)=0.0
!
!Correction dans les mailles où ds a été surestimée
ZDZRC2(:)=PRC(:)-ZDZRC2(:)
WHERE ((MEM_PS0(:).LE.0.0).AND.(PS0(:).GT.0.0).AND.(ZDZRC2(:).LT.0.0))
  PS0(:)=0.0
  ADJU2(:)=1.0
ENDWHERE
!
WHERE ((MEM_PS0(:).GE.0.0).AND.(PS0(:).LT.0.0).AND.(ZDZRC2(:).GT.0.0))
  PS0(:)=0.0
  ADJU2(:)=1.0
ENDWHERE
!
DO J2=1,SIZE(ADJU2,1)
  IF (ADJU2(J2)==1) THEN 
   DO J1=1,5
    ZLV(J2) = XLVTT+(XCPV-XCL)*(PTT(J2)-XTT)                
    ZCPH(J2)=XCPD+XCPV*PRV(J2)+XCL*(PRC(J2)+PRR(J2))
    ZZW1(J2)=EXP(XALPW-XBETAW/PTT(J2)-XGAMW*LOG(PTT(J2)))
    ZRVSAT1(J2)=(XMV / XMD)*ZZW1(J2)/(PPRES(J2)-ZZW1(J2))
    PRV(J2)=MIN(ZRT(J2),(PS0(J2)+1.0)*ZRVSAT1(J2))
    PRC(J2)=MAX(ZRT(J2)-PRV(J2),0.0)
    PTT(J2)=0.5*PTT(J2)+0.5*(ZTL(J2)+(ZLV(J2)*PRC(J2)/ZCPH(J2)))
   ENDDO
   ZLV(J2) = XLVTT+(XCPV-XCL)*(PTT(J2)-XTT)                
   ZCPH(J2)=XCPD+XCPV*PRV(J2)+XCL*(PRC(J2)+PRR(J2))
   PTT(J2)=ZTL(J2)+(ZLV(J2)*PRC(J2)/ZCPH(J2))
  ENDIF
ENDDO
!
!Elimination de l'eau liquide dans les mailles où le rayon des gouttelettes est
!inférieur à AER_RAD
ZRMOY(:)=0.0
DO J2=1,SIZE(PRC,1)
 IF (PRC(J2).NE.0.0) THEN
  ZRMOY(J2)=(MOMG(XALPHAC,XNUC,3.0)*4.0*XPI*PCL(J2)*XRHOLW/&
        (3.0*PRC(J2)*PRHOD(J2)))**(1.0/3.0)
  ZRMOY(J2)=MOMG(XALPHAC,XNUC,1.0)/ZRMOY(J2)
IF ((ZFLAG_ACT(J2).EQ.1.0).AND.(MEM_PS0(J2).LT.0.0).AND.(ZRMOY(J2).LT.AER_RAD)) THEN
   PTT(J2)=ZTL(J2)
   PRV(J2)=ZRT(J2)
   PRC(J2)=0.0
  ENDIF
 ENDIF
ENDDO
!
!Calcul de S au regard de T et rv en fin de pas de temps
ZZW1=EXP(XALPW-XBETAW/PTT(:)-XGAMW*LOG(PTT(:)))
 !rvsat
ZRVSAT1(:)=(XMV / XMD)*ZZW1(:)/(PPRES-ZZW1(:))
!
WHERE (PRC(:)==0.0D0)
 PS0(:)=(PRV(:)/ZRVSAT1(:))-1D0
ENDWHERE
!
CONTAINS
!
FUNCTION MOMG (PALPHA,PNU,PP) RESULT (PMOMG)
USE MODI_GAMMA
IMPLICIT NONE
REAL     :: PALPHA ! first shape parameter of the DIMENSIONnal distribution
REAL     :: PNU    ! second shape parameter of the DIMENSIONnal distribution
REAL     :: PP     ! order of the moment
REAL     :: PMOMG  ! result: moment of order ZP
PMOMG = GAMMA(PNU+PP/PALPHA)/GAMMA(PNU)
!
END FUNCTION MOMG
!
END SUBROUTINE PROGNOS
