MODULE lmdz_wake_ini
IMPLICIT NONE 

  ! ============================================================================


  ! But : Decrire le comportement des poches froides apparaissant dans les
  ! grands systemes convectifs, et fournir l'energie disponible pour
  ! le declenchement de nouvelles colonnes convectives.

  ! State variables : 
  ! deltatw    : temperature difference between wake and off-wake regions
  ! deltaqw    : specific humidity difference between wake and off-wake regions
  ! sigmaw     : fractional area covered by wakes.
  ! wdens      : number of wakes per unit area

  ! -------------------------------------------------------------------------
  ! Declaration de variables
  ! -------------------------------------------------------------------------

  ! Variables a fixer
!jyg<
!!  REAL, SAVE                                            :: stark, wdens_ref, coefgw, alpk
  INTEGER, SAVE, PROTECTED                                    :: prt_level
  REAL, SAVE, PROTECTED, DIMENSION(2)                         :: wdens_ref
  REAL, SAVE, PROTECTED                                       :: stark, coefgw, alpk, wk_pupper
!>jyg
  REAL, SAVE, PROTECTED                                       :: crep_upper, crep_sol  
  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, wk_pupper, crep_upper, crep_sol)

  REAL, SAVE, PROTECTED                                       :: tau_cv
  !$OMP THREADPRIVATE(tau_cv)

   REAL, SAVE, PROTECTED                                       :: wk_delta_t_min
  !$OMP THREADPRIVATE(wk_delta_t_min)

   REAL, SAVE, PROTECTED                                       :: wk_frac_int_delta_t
  !$OMP THREADPRIVATE(wk_frac_int_delta_t)

  REAL, SAVE, PROTECTED                                       :: rzero, aa0 ! minimal wake radius and area
  !$OMP THREADPRIVATE(rzero, aa0)

  LOGICAL, SAVE, PROTECTED                                    :: ok_bug_gfl
  !$OMP THREADPRIVATE(ok_bug_gfl)
  LOGICAL, SAVE, PROTECTED                                    :: flag_wk_check_trgl
  !$OMP THREADPRIVATE(flag_wk_check_trgl)
  INTEGER, SAVE, PROTECTED                                    :: iflag_wk_act
  !$OMP THREADPRIVATE(iflag_wk_act)

  INTEGER, SAVE, PROTECTED                                    :: iflag_wk_check_trgl
  !$OMP THREADPRIVATE(iflag_wk_check_trgl)
  INTEGER, SAVE, PROTECTED                                    :: iflag_wk_pop_dyn
  !$OMP THREADPRIVATE(iflag_wk_pop_dyn)

  INTEGER, SAVE, PROTECTED                                    :: iflag_wk_profile
  !$OMP THREADPRIVATE(iflag_wk_profile)

  INTEGER, SAVE, PROTECTED                                    :: wk_nsub
  !$OMP THREADPRIVATE(wk_nsub)

  INTEGER, SAVE, PROTECTED                                    :: iflag_wk_new_ptop
  !$OMP THREADPRIVATE(iflag_wk_new_ptop)

  REAL, SAVE, PROTECTED                                       :: wdensinit ! Minimum wake density used to restart wakes from a wake-free state
  !$OMP THREADPRIVATE(wdensinit)
  REAL, SAVE, PROTECTED                                       :: wdensthreshold ! Threshold wake density below which wakes are killed
  !$OMP THREADPRIVATE(wdensthreshold)
  REAL, SAVE, PROTECTED                                       :: sigmad, hwmin, wapecut, cstart
  !$OMP THREADPRIVATE(sigmad, hwmin, wapecut, cstart)
  REAL, SAVE, PROTECTED                                       :: sigmaw_max
  !$OMP THREADPRIVATE(sigmaw_max)  
  REAL, SAVE, PROTECTED                                       :: dens_rate
  !$OMP THREADPRIVATE(dens_rate)
  REAL, SAVE, PROTECTED                                       :: epsilon_loc
  !$OMP THREADPRIVATE(epsilon_loc)
  REAL, SAVE, PROTECTED                                       :: epsim1,RG,RD
  !$OMP THREADPRIVATE(epsim1,RG,RD)
  REAL, SAVE, PROTECTED                                        ::smallestreal
  !$OMP THREADPRIVATE(smallestreal)
  REAL, SAVE, PROTECTED                                        :: wk_int_delta_t_min
  !$OMP THREADPRIVATE(wk_int_delta_t_min)



CONTAINS

  ! =========================================================================
  SUBROUTINE wake_ini(rg_in,rd_in,rv_in,prt_lev)
  ! =========================================================================

  ! **************************************************************
  ! *
  ! WAKE                                                        *
  ! retour a un Pupper fixe                                *
  ! *
  ! written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
  ! modified by :   ROEHRIG Romain        01/29/2007            *
  ! **************************************************************

  ! -------------------------------------------------------------------------
  ! Initialisations
  ! -------------------------------------------------------------------------

  USE ioipsl_getin_p_mod, ONLY: getin_p
  real eps
  integer, intent(in) :: prt_lev
  real, intent(in) :: rg_in,rd_in,rv_in

  smallestreal=tiny(smallestreal)

  prt_level=prt_lev
  epsilon_loc=1.E-15
  wapecut=1. ! previously 5.

  rzero = 5000.
  CALL getin_p('rzero_wk', rzero)
  aa0 = 3.14*rzero*rzero

  ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
!!  sigmad=0.005
  sigmad=0.02
  CALL getin_p('sigmad', sigmad)
  hwmin=10.

!!wdensthreshold=1.e-12
  wdensthreshold=1.e-14
  wdensthreshold=2.e-14
  CALL getin_p('wdensthreshold', wdensthreshold)

  IF (sigmad < 0.) THEN
    sigmad = abs(sigmad)
!!    wdensmin=sigmad/(3.14*rzero**2)
    wdensinit=sigmad/(3.14*rzero**2)
  ELSE
    wdensinit = wdensthreshold/2.
  ENDIF


  ! cc nrlmd
  sigmaw_max=0.4
  dens_rate=0.1

  eps = rd_in/rv_in
  epsim1 = 1.0/eps - 1.0
  RG=rg_in
  RD=rd_in


  ! cc
  ! Longueur de maille (en m)
  ! -------------------------------------------------------------------------

  ! ALON = 3.e5
  ! alon = 1.E6

  ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)

  ! coefgw : Coefficient pour les ondes de gravite
  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
  ! wdens : Densite surfacique de poche froide
  ! -------------------------------------------------------------------------

  ! cc nrlmd      coefgw=10
  ! coefgw=1
  ! wdens0 = 1.0/(alon**2)
  ! cc nrlmd      wdens = 1.0/(alon**2)
  ! cc nrlmd      stark = 0.50
  ! CRtest
  ! cc nrlmd      alpk=0.1
  ! alpk = 1.0
  ! alpk = 0.5
  ! alpk = 0.05



  crep_upper = 0.9
  crep_sol = 1.0

  ! Flag concerning the bug in gfl computation
  ok_bug_gfl = .True.
  call getin_p('ok_bug_gfl', ok_bug_gfl)

  ! Get wapecut from parameter file
  wapecut = 1.

print*,'wapecut',wapecut
  CALL getin_p('wapecut', wapecut)
print*,'wapecut',wapecut

  ! cc nrlmd Lecture du fichier wake_param.data


  ! cc nrlmd Lecture du fichier wake_param.data
  stark=0.33
  CALL getin_p('stark',stark)
  cstart = stark*sqrt(2.*wapecut)

  alpk=0.25
  CALL getin_p('alpk',alpk)
  
  wk_pupper=0.6
  CALL getin_p('wk_pupper',wk_pupper)


!jyg<
!!  wdens_ref=8.E-12
!!  CALL getin_p('wdens_ref',wdens_ref)
  wdens_ref(1)=8.E-12
  wdens_ref(2)=8.E-12
  CALL getin_p('wdens_ref_o',wdens_ref(1))    !wake number per unit area ; ocean
  CALL getin_p('wdens_ref_l',wdens_ref(2))    !wake number per unit area ; land
!>jyg

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!  Population dynamics parameters    !!!!!!!!!!!!!!!!!!!!!!!!!!!!
!------------------------------------------------------------------------ 

  iflag_wk_pop_dyn = 0
  CALL getin_p('iflag_wk_pop_dyn',iflag_wk_pop_dyn) ! switch between wdens prescribed 
                                                    ! and wdens prognostic
  iflag_wk_act = 0
  CALL getin_p('iflag_wk_act',iflag_wk_act) ! 0: act(:)=0.
                                            ! 1: act(:)=1.
                                            ! 2: act(:)=f(Wape)

  iflag_wk_profile = 0
  CALL getin_p('iflag_wk_profile',iflag_wk_profile) ! switch between wdens prescribed 
                                                    ! and wdens prognostic
 ! iflag_wk_profile = 0
  iflag_wk_new_ptop = 0
  CALL getin_p('iflag_wk_new_ptop',iflag_wk_new_ptop)

  wk_nsub = 10
  CALL getin_p('wk_nsub',wk_nsub)

  tau_cv = 4000.
  CALL getin_p('tau_cv', tau_cv)
  
  wk_delta_t_min = 0.
  CALL getin_p('wk_delta_t_min', wk_delta_t_min)
  
  wk_int_delta_t_min = 10.
  CALL getin_p('wk_int_delta_t_min', wk_int_delta_t_min)
  
  wk_frac_int_delta_t = 0.9
  CALL getin_p('wk_frac_int_delta_t', wk_frac_int_delta_t)


!------------------------------------------------------------------------ 

  coefgw=4.
  CALL getin_p('coefgw',coefgw)

  WRITE(*,*) 'stark=', stark
  WRITE(*,*) 'alpk=', alpk
  WRITE(*,*) 'wk_pupper=', wk_pupper
!jyg<
!!  WRITE(*,*) 'wdens_ref=', wdens_ref
  WRITE(*,*) 'wdens_ref_o=', wdens_ref(1)
  WRITE(*,*) 'wdens_ref_l=', wdens_ref(2)
!>jyg
  WRITE(*,*) 'iflag_wk_pop_dyn=',iflag_wk_pop_dyn
  WRITE(*,*) 'iflag_wk_act',iflag_wk_act
  WRITE(*,*) 'coefgw=', coefgw 

  flag_wk_check_trgl=.false.
  CALL getin_p('flag_wk_check_trgl ', flag_wk_check_trgl)
  WRITE(*,*) 'flag_wk_check_trgl=', flag_wk_check_trgl
  WRITE(*,*) 'flag_wk_check_trgl OBSOLETE. Utilisr iflag_wk_check_trgl plutot'
  iflag_wk_check_trgl=0 ; IF (flag_wk_check_trgl) iflag_wk_check_trgl=1
  CALL getin_p('iflag_wk_check_trgl ', iflag_wk_check_trgl)
  WRITE(*,*) 'iflag_wk_check_trgl=', iflag_wk_check_trgl
  
  WRITE(*,*) 'wk_delta_t_min=', wk_delta_t_min
  WRITE(*,*) 'wk_int_delta_t_min=', wk_int_delta_t_min
  WRITE(*,*) 'wk_frac_int_delta_t=', wk_frac_int_delta_t
  WRITE(*,*) 'iflag_wk_new_ptop=', iflag_wk_new_ptop
  WRITE(*,*) 'wk_nsub=', wk_nsub

 RETURN

END SUBROUTINE wake_ini



END MODULE lmdz_wake_ini
