source: LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90 @ 4804

MODULE lmdz_lscp_ini

IMPLICIT NONE

  ! PARAMETERS for lscp:
  !--------------------
 
  REAL RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
  !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG)

  REAL, SAVE, PROTECTED :: seuil_neb=0.001      ! cloud fraction threshold: a cloud really exists when exceeded
  !$OMP THREADPRIVATE(seuil_neb)

  INTEGER, SAVE :: lunout, prt_level            ! Logical unit number and level for standard output
  !$OMP THREADPRIVATE(lunout,prt_level)

  INTEGER, SAVE, PROTECTED :: niter_lscp=5      ! number of iterations to calculate autoconversion to precipitation
  !$OMP THREADPRIVATE(niter_lscp)

  INTEGER, SAVE, PROTECTED :: iflag_evap_prec=1 ! precipitation evaporation flag. 0: nothing, 1: "old way", 
                                                ! 2: Max cloud fraction above to calculate the max of reevaporation
                                                ! >=4: LTP'method i.e. evaporation in the clear-sky fraction of the mesh only
                                                ! pay attention that iflag_evap_prec=4 may lead to unrealistic and overestimated
                                                ! evaporation. Use 5 instead
  !$OMP THREADPRIVATE(iflag_evap_prec)

  REAL, SAVE, PROTECTED :: t_coup=234.0         ! temperature threshold which determines the phase 
                                                ! for which the saturation vapor pressure is calculated
  !$OMP THREADPRIVATE(t_coup)
  REAL, SAVE, PROTECTED :: DDT0=0.01            ! iteration parameter 
  !$OMP THREADPRIVATE(DDT0)

  REAL, SAVE, PROTECTED :: ztfondue=278.15      ! parameter to calculate melting fraction of precipitation
  !$OMP THREADPRIVATE(ztfondue)

  REAL, SAVE, PROTECTED :: temp_nowater=233.15  ! temperature below which liquid water no longer exists
  !$OMP THREADPRIVATE(temp_nowater)

  REAL, SAVE, PROTECTED :: a_tr_sca(4)          ! Variables for tracers temporary: alpha parameter for scavenging, 4 possible scavenging processes
  !$OMP THREADPRIVATE(a_tr_sca)
  
  INTEGER, SAVE, PROTECTED ::  iflag_mpc_bl=0   ! flag to activate boundary layer mixed phase cloud param
  !$OMP THREADPRIVATE(iflag_mpc_bl)
  
  LOGICAL, SAVE, PROTECTED :: ok_radocond_snow=.false. ! take into account the mass of ice precip in the cloud ice content seen by radiation
  !$OMP THREADPRIVATE(ok_radocond_snow)

  REAL, SAVE, PROTECTED :: t_glace_min=258.0    ! lower-bound temperature parameter for cloud phase determination
  !$OMP THREADPRIVATE(t_glace_min)

  REAL, SAVE, PROTECTED :: t_glace_max=273.15   ! upper-bound temperature parameter for cloud phase determination
  !$OMP THREADPRIVATE(t_glace_max)

  REAL, SAVE, PROTECTED :: exposant_glace=1.0   ! parameter for cloud phase determination
  !$OMP THREADPRIVATE(exposant_glace)

  INTEGER, SAVE, PROTECTED :: iflag_vice=0      ! which expression for ice crystall fall velocity
  !$OMP THREADPRIVATE(iflag_vice)

  INTEGER, SAVE, PROTECTED :: iflag_t_glace=0   ! which expression for cloud phase partitioning
  !$OMP THREADPRIVATE(iflag_t_glace)

  INTEGER, SAVE, PROTECTED :: iflag_cloudth_vert=0         ! option for determining cloud fraction and content in convective boundary layers
  !$OMP THREADPRIVATE(iflag_cloudth_vert)

  INTEGER, SAVE, PROTECTED :: iflag_gammasat=0             ! which threshold for homogeneous nucleation below -40oC
  !$OMP THREADPRIVATE(iflag_gammasat)

  INTEGER, SAVE, PROTECTED :: iflag_rain_incloud_vol=0     ! use of volume cloud fraction for rain autoconversion 
  !$OMP THREADPRIVATE(iflag_rain_incloud_vol)

  INTEGER, SAVE, PROTECTED :: iflag_bergeron=0             ! bergeron effect for liquid precipitation treatment  
  !$OMP THREADPRIVATE(iflag_bergeron)

  INTEGER, SAVE, PROTECTED :: iflag_fisrtilp_qsat=0        ! qsat adjustment (iterative) during autoconversion
  !$OMP THREADPRIVATE(iflag_fisrtilp_qsat)

  INTEGER, SAVE, PROTECTED :: iflag_pdf=0                  ! type of subgrid scale qtot pdf
  !$OMP THREADPRIVATE(iflag_pdf)

  INTEGER, SAVE, PROTECTED :: iflag_autoconversion=0       ! autoconversion option
  !$OMP THREADPRIVATE(iflag_autoconversion)

  LOGICAL, SAVE, PROTECTED :: reevap_ice=.false.           ! no liquid precip for T< threshold
  !$OMP THREADPRIVATE(reevap_ice)

  REAL, SAVE, PROTECTED :: cld_lc_lsc=2.6e-4               ! liquid autoconversion coefficient, stratiform rain
  !$OMP THREADPRIVATE(cld_lc_lsc)

  REAL, SAVE, PROTECTED :: cld_lc_con=2.6e-4                ! liquid autoconversion coefficient, convective rain
  !$OMP THREADPRIVATE(cld_lc_con)

  REAL, SAVE, PROTECTED :: cld_tau_lsc=3600.                ! liquid autoconversion timescale, stratiform rain
  !$OMP THREADPRIVATE(cld_tau_lsc)

  REAL, SAVE, PROTECTED :: cld_tau_con=3600.                ! liquid autoconversion timescale, convective rain
  !$OMP THREADPRIVATE(cld_tau_con)

  REAL, SAVE, PROTECTED :: cld_expo_lsc=2.                  ! liquid autoconversion threshold exponent, stratiform rain
  !$OMP THREADPRIVATE(cld_expo_lsc)

  REAL, SAVE, PROTECTED :: cld_expo_con=2.                  ! liquid autoconversion threshold exponent, convective rain
  !$OMP THREADPRIVATE(cld_expo_con)

  REAL, SAVE, PROTECTED :: ffallv_lsc=1.                    ! tuning coefficient crystal fall velocity, stratiform 
  !$OMP THREADPRIVATE(ffallv_lsc)

  REAL, SAVE, PROTECTED :: ffallv_con=1.                    ! tuning coefficient crystal fall velocity, convective
  !$OMP THREADPRIVATE(ffallv_con)

  REAL, SAVE, PROTECTED :: coef_eva=2e-5                    ! tuning coefficient liquid precip evaporation
  !$OMP THREADPRIVATE(coef_eva)

  REAL, SAVE, PROTECTED :: coef_eva_i                       ! tuning coefficient ice precip sublimation
  !$OMP THREADPRIVATE(coef_eva_i)

  REAL, SAVE, PROTECTED :: expo_eva=0.5                     ! tuning coefficient liquid precip evaporation
  !$OMP THREADPRIVATE(expo_eva)

  REAL, SAVE, PROTECTED :: expo_eva_i                       ! tuning coefficient ice precip sublimation
  !$OMP THREADPRIVATE(expo_eva_i)

  REAL, SAVE, PROTECTED :: cice_velo=1.645                  ! factor in the ice fall velocity formulation
  !$OMP THREADPRIVATE(cice_velo)

  REAL, SAVE, PROTECTED ::  dice_velo=0.16                  ! exponent in the ice fall velocity formulation
  !$OMP THREADPRIVATE(dice_velo)

  REAL, SAVE, PROTECTED :: dist_liq=300.                    ! typical deph of cloud-top liquid layer in mpcs
  !$OMP THREADPRIVATE(dist_liq)

  REAL, SAVE, PROTECTED    :: tresh_cl=0.0                  ! cloud fraction threshold for cloud top search
  !$OMP THREADPRIVATE(tresh_cl)

  !--Parameters for poprecip
  LOGICAL, SAVE, PROTECTED :: ok_poprecip=.FALSE.           ! use the processes-oriented formulation of precipitations
  !$OMP THREADPRIVATE(ok_poprecip)

  REAL, SAVE, PROTECTED :: rain_int_min=0.001               ! Minimum local rain intensity [mm/s] before the decrease in associated precipitation fraction
  !$OMP THREADPRIVATE(rain_int_min)

  REAL, SAVE, PROTECTED :: thresh_precip_frac=1.E-6         ! precipitation fraction threshold TODO [-]
  !$OMP THREADPRIVATE(thresh_precip_frac)

  REAL, SAVE, PROTECTED :: gamma_col=1.                     ! A COMMENTER TODO [-]
  !$OMP THREADPRIVATE(gamma_col)

  REAL, SAVE, PROTECTED :: gamma_agg=1.                     ! A COMMENTER TODO [-]
  !$OMP THREADPRIVATE(gamma_agg)

  REAL, SAVE, PROTECTED :: gamma_rim=1.                     ! A COMMENTER TODO [-]
  !$OMP THREADPRIVATE(gamma_rim)

  REAL, SAVE, PROTECTED :: rho_rain=1000.                    ! A COMMENTER TODO [kg/m3]
  !$OMP THREADPRIVATE(rho_rain)

  REAL, SAVE, PROTECTED :: rho_snow=500.                   ! A COMMENTER TODO [kg/m3]
  !$OMP THREADPRIVATE(rho_snow)

  REAL, SAVE, PROTECTED :: r_rain=100.E-6                   ! A COMMENTER TODO [m]
  !$OMP THREADPRIVATE(r_rain)

  REAL, SAVE, PROTECTED :: r_snow=100.E-6                    ! A COMMENTER TODO [m]
  !$OMP THREADPRIVATE(r_snow)

  REAL, SAVE, PROTECTED :: Eff_rain_liq=1.0                  ! A COMMENTER TODO [-]
  !$OMP THREADPRIVATE(Eff_rain_liq)

  REAL, SAVE, PROTECTED :: Eff_snow_ice=0.5                ! A COMMENTER TODO [-]
  !$OMP THREADPRIVATE(Eff_snow_ice)

  REAL, SAVE, PROTECTED :: Eff_snow_liq=1.0              ! A COMMENTER TODO [-]
  !$OMP THREADPRIVATE(Eff_snow_liq)

  REAL, SAVE, PROTECTED :: tau_auto_snow_min=1800.          ! A COMMENTER TODO [s]
  !$OMP THREADPRIVATE(tau_auto_snow_min)

  REAL, SAVE, PROTECTED :: tau_auto_snow_max=7200.          ! A COMMENTER TODO [s]
  !$OMP THREADPRIVATE(tau_auto_snow_max)

  REAL, SAVE, PROTECTED :: eps=1.E-10                       ! A COMMENTER TODO [-]
  !$OMP THREADPRIVATE(eps)
  !--End of the parameters for poprecip

! Two parameters used for lmdz_lscp_old only
  INTEGER, SAVE, PROTECTED :: iflag_oldbug_fisrtilp=0, fl_cor_ebil 
  !$OMP THREADPRIVATE(iflag_oldbug_fisrtilp,fl_cor_ebil)

CONTAINS

SUBROUTINE lscp_ini(dtime,lunout_in,prt_level_in,ok_ice_sursat, iflag_ratqs, fl_cor_ebil_in, RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, &
                    RVTMP2_in, RTT_in,RD_in,RG_in)


   USE ioipsl_getin_p_mod, ONLY : getin_p
   USE ice_sursat_mod, ONLY: ice_sursat_init
   USE lmdz_cloudth_ini, ONLY : cloudth_ini

   REAL, INTENT(IN)      :: dtime
   INTEGER, INTENT(IN)   :: lunout_in,prt_level_in,iflag_ratqs,fl_cor_ebil_in
   LOGICAL, INTENT(IN)   :: ok_ice_sursat

   REAL, INTENT(IN)      :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in
   REAL, INTENT(IN)      ::  RVTMP2_in, RTT_in, RD_in, RG_in
   character (len=20) :: modname='lscp_ini_mod'
   character (len=80) :: abort_message


    lunout=lunout_in
    prt_level=prt_level_in
    fl_cor_ebil=fl_cor_ebil_in

    RG=RG_in
    RD=RD_in
    RCPD=RCPD_in
    RLVTT=RLVTT_in
    RLSTT=RLSTT_in
    RLMLT=RLMLT_in
    RTT=RTT_in
    RG=RG_in
    RVTMP2=RVTMP2_in



    CALL getin_p('niter_lscp',niter_lscp)
    CALL getin_p('iflag_evap_prec',iflag_evap_prec)
    CALL getin_p('seuil_neb',seuil_neb)
    CALL getin_p('iflag_mpc_bl',iflag_mpc_bl)
    CALL getin_p('ok_radocond_snow',ok_radocond_snow) 
    CALL getin_p('t_glace_max',t_glace_max)
    CALL getin_p('t_glace_min',t_glace_min)
    CALL getin_p('exposant_glace',exposant_glace) 
    CALL getin_p('iflag_vice',iflag_vice)
    CALL getin_p('iflag_t_glace',iflag_t_glace)
    CALL getin_p('iflag_cloudth_vert',iflag_cloudth_vert)
    CALL getin_p('iflag_gammasat',iflag_gammasat)
    CALL getin_p('iflag_rain_incloud_vol',iflag_rain_incloud_vol)
    CALL getin_p('iflag_bergeron',iflag_bergeron)
    CALL getin_p('iflag_fisrtilp_qsat',iflag_fisrtilp_qsat)
    CALL getin_p('iflag_pdf',iflag_pdf)
    CALL getin_p('reevap_ice',reevap_ice)
    CALL getin_p('cld_lc_lsc',cld_lc_lsc)
    CALL getin_p('cld_lc_con',cld_lc_con)
    CALL getin_p('cld_tau_lsc',cld_tau_lsc)
    CALL getin_p('cld_tau_con',cld_tau_con)
    CALL getin_p('cld_expo_lsc',cld_expo_lsc)
    CALL getin_p('cld_expo_con',cld_expo_con)
    CALL getin_p('ffallv_lsc',ffallv_lsc)
    CALL getin_p('ffallv_lsc',ffallv_con)
    CALL getin_p('coef_eva',coef_eva)
    coef_eva_i=coef_eva
    CALL getin_p('coef_eva_i',coef_eva_i)
    CALL getin_p('expo_eva',expo_eva)
    expo_eva_i=expo_eva
    CALL getin_p('expo_eva_i',expo_eva_i)
    CALL getin_p('iflag_autoconversion',iflag_autoconversion)
    CALL getin_p('dist_liq',dist_liq)
    CALL getin_p('tresh_cl',tresh_cl)
    CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
    CALL getin_p('ok_poprecip',ok_poprecip)
    CALL getin_p('rain_int_min',rain_int_min) 
    CALL getin_p('gamma_agg',gamma_agg)
    CALL getin_p('gamma_col',gamma_col)



    WRITE(lunout,*) 'lscp_ini, niter_lscp:', niter_lscp
    WRITE(lunout,*) 'lscp_ini, iflag_evap_prec:', iflag_evap_prec
    WRITE(lunout,*) 'lscp_ini, seuil_neb:', seuil_neb
    WRITE(lunout,*) 'lscp_ini, iflag_mpc_bl:', iflag_mpc_bl
    WRITE(lunout,*) 'lscp_ini, ok_radocond_snow:', ok_radocond_snow
    WRITE(lunout,*) 'lscp_ini, t_glace_max:', t_glace_max
    WRITE(lunout,*) 'lscp_ini, t_glace_min:', t_glace_min
    WRITE(lunout,*) 'lscp_ini, exposant_glace:', exposant_glace
    WRITE(lunout,*) 'lscp_ini, iflag_vice:', iflag_vice
    WRITE(lunout,*) 'lscp_ini, iflag_t_glace:', iflag_t_glace
    WRITE(lunout,*) 'lscp_ini, iflag_cloudth_vert:', iflag_cloudth_vert
    WRITE(lunout,*) 'lscp_ini, iflag_gammasat:', iflag_gammasat
    WRITE(lunout,*) 'lscp_ini, iflag_rain_incloud_vol:', iflag_rain_incloud_vol
    WRITE(lunout,*) 'lscp_ini, iflag_bergeron:', iflag_bergeron
    WRITE(lunout,*) 'lscp_ini, iflag_fisrtilp_qsat:', iflag_fisrtilp_qsat
    WRITE(lunout,*) 'lscp_ini, iflag_pdf', iflag_pdf
    WRITE(lunout,*) 'lscp_ini, reevap_ice', reevap_ice
    WRITE(lunout,*) 'lscp_ini, cld_lc_lsc', cld_lc_lsc
    WRITE(lunout,*) 'lscp_ini, cld_lc_con', cld_lc_con
    WRITE(lunout,*) 'lscp_ini, cld_tau_lsc', cld_tau_lsc
    WRITE(lunout,*) 'lscp_ini, cld_tau_con', cld_tau_con
    WRITE(lunout,*) 'lscp_ini, cld_expo_lsc', cld_expo_lsc
    WRITE(lunout,*) 'lscp_ini, cld_expo_con', cld_expo_con
    WRITE(lunout,*) 'lscp_ini, ffallv_lsc', ffallv_lsc
    WRITE(lunout,*) 'lscp_ini, ffallv_con', ffallv_con
    WRITE(lunout,*) 'lscp_ini, coef_eva', coef_eva
    WRITE(lunout,*) 'lscp_ini, coef_eva_i', coef_eva_i
    WRITE(lunout,*) 'lscp_ini, expo_eva', expo_eva
    WRITE(lunout,*) 'lscp_ini, expo_eva_i', expo_eva_i
    WRITE(lunout,*) 'lscp_ini, iflag_autoconversion', iflag_autoconversion
    WRITE(lunout,*) 'lscp_ini, dist_liq', dist_liq
    WRITE(lunout,*) 'lscp_ini, tresh_cl', tresh_cl
    WRITE(lunout,*) 'lscp_ini, iflag_oldbug_fisrtilp', iflag_oldbug_fisrtilp
    WRITE(lunout,*) 'lscp_ini, fl_cor_ebil', fl_cor_ebil
    WRITE(lunout,*) 'lscp_ini, ok_poprecip', ok_poprecip
    WRITE(lunout,*) 'lscp_ini, rain_int_min:', rain_int_min
    WRITE(lunout,*) 'lscp_ini, gamma_agg:', gamma_agg
    WRITE(lunout,*) 'lscp_ini, gamma_col:', gamma_col





    ! check for precipitation sub-time steps
    IF (ABS(dtime/REAL(niter_lscp)-360.0).GT.0.001) THEN
        WRITE(lunout,*) 'lscp: it is not expected, see Z.X.Li', dtime
        WRITE(lunout,*) 'I would prefer a 6 min sub-timestep'
    ENDIF

    ! check consistency between numerical resolution of autoconversion 
    ! and other options
   
    IF (iflag_autoconversion .EQ. 2) THEN
        IF ((iflag_vice .NE. 0) .OR. (niter_lscp .GT. 1)) THEN
           abort_message = 'in lscp, iflag_autoconversion=2 requires iflag_vice=0 and niter_lscp=1'
           CALL abort_physic (modname,abort_message,1) 
        ENDIF 
    ENDIF


    ! check consistency between the use of the processes-oriented precipitation formulations
    ! and other options
   
    IF (ok_poprecip) THEN
        IF ((iflag_evap_prec .LT. 4) .OR. (niter_lscp .GT. 1)) THEN
           abort_message = 'in lscp, ok_poprecip=y requires iflag_evap_prec >= 4 and niter_lscp=1'
           CALL abort_physic (modname,abort_message,1) 
        ENDIF 
    ENDIF


    !AA Temporary initialisation
    a_tr_sca(1) = -0.5
    a_tr_sca(2) = -0.5
    a_tr_sca(3) = -0.5
    a_tr_sca(4) = -0.5
    
    IF (ok_ice_sursat) CALL ice_sursat_init()

    CALL cloudth_ini(iflag_cloudth_vert,iflag_ratqs)

RETURN

END SUBROUTINE lscp_ini

END MODULE lmdz_lscp_ini