module lmdz_blosno_ini

   implicit none

   real, save, protected :: RCPD, RV, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI
   real, save, protected :: coef_sub_bs, fallv_bs, zeta_bs, c_esalt_bs
   real, save, protected :: prt_bs, pbst_bs, qbst_bs, r_bs
   integer, save, protected :: iflag_saltation_bs, iflag_sedim_bs, iflag_sublim_bs

   !$OMP THREADPRIVATE(RCPD, RV, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG, RPI)
   !$OMP THREADPRIVATE(coef_sub_bs, fallv_bs, r_bs, zeta_bs, c_esalt_bs)
   !$OMP THREADPRIVATE(pbst_bs, prt_bs, qbst_bs)
   !$OMP THREADPRIVATE(iflag_saltation_bs,iflag_sedim_bs, iflag_sublim_bs)

   real, save, protected :: tbsmelt = 278.15    ! parameter to calculate melting fraction of BS sedimentation
   real, save, protected :: taumeltbs0 = 600.0  ! Melting time scale of blowing snow at 273.15K
   real, save, protected :: qbmin = 1.E-10      ! Minimum blowing snow specific content
   !$OMP THREADPRIVATE(tbsmelt, taumeltbs0, qbmin)

   real, save, protected :: tau_dens0_bs = 864000.      ! 10 days by default, in s
   real, save, protected :: tau_densmin_bs = 21600.    ! 1/4 days according to in situ obs by C. Amory during blowing snow +
   ! Marshall et al. 1999 (snow densification during rain)
   real, save, protected :: tau_eqsalt_bs = 10.        ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt of about 10s
   real, save, protected :: rhofresh_bs = 300.0       ! fresh snow density kg/m3
   real, save, protected :: rhohard_bs = 450.0       ! hard snow density kg/m3
   real, save, protected :: rhoice_bs = 920.0         ! ice density kg/m3
   real, save, protected :: rhobs = 900.0               ! blowing snow density (kg/m3) following Bintanja et al. 2001 part I
   !$OMP THREADPRIVATE(rhoice_bs, rhofresh_bs, rhohard_bs, tau_dens0_bs, tau_densmin_bs, tau_eqsalt_bs, rhobs)

contains

   subroutine blosno_ini(RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, &
                               RVTMP2_in, RTT_in, RD_in, RG_in, RV_in, RPI_in)

      USE ioipsl_getin_p_mod, ONLY: getin_p

      real, intent(in) :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, RPI_in
      real, intent(in) ::  RVTMP2_in, RTT_in, RD_in, RG_in, RV_in

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

      c_esalt_bs = 3.25
      CALL getin_p('c_esalt_bs', c_esalt_bs)

      qbst_bs = 0.001
      CALL getin_p('qbst_bs', qbst_bs)

      pbst_bs = 0.00003
      CALL getin_p('pbst_bs', pbst_bs)

      prt_bs = 0.00003
      CALL getin_p('prt_bs', prt_bs)

      zeta_bs = 1.
      CALL getin_p('zeta_bs', zeta_bs)

      fallv_bs = 0.5
      CALL getin_p('fallv_bs', fallv_bs)

      coef_sub_bs = 0.01
      CALL getin_p('coef_sub_bs', coef_sub_bs)

      iflag_sublim_bs = 1
      CALL getin_p('iflag_sublim_bs', iflag_sublim_bs)

      iflag_sedim_bs = 1
      CALL getin_p('iflag_sedim_bs', iflag_sedim_bs)

      r_bs = 50.0e-6
      CALL getin_p('r_bs', r_bs)

   end subroutine blosno_ini

end module lmdz_blosno_ini