MODULE lmdz_wake_popdyn_2
PUBLIC wake_popdyn_2
CONTAINS

    SUBROUTINE wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, &
                             wdensmin, &
                             sigmaw, wdens, awdens, &   !! states variables
                             gfl, cstar, cin, wape, rad_wk, &
                             d_sigmaw, d_wdens, d_awdens, &  !! tendences 
                             cont_fact, &
                             d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
                             d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
                             d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
                             
                                             

  USE lmdz_wake_ini , ONLY : wake_ini
  USE lmdz_wake_ini , ONLY : prt_level,RG
  USE lmdz_wake_ini , ONLY : stark, wdens_ref
  USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0
!!  USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin
  USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn
  USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max
  
IMPLICIT NONE

  INTEGER, INTENT(IN)                                   :: klon,klev
  LOGICAL, DIMENSION (klon),        INTENT(IN)          :: wk_adv
  REAL,                             INTENT(IN)          :: dtimesub
  REAL,                             INTENT(IN)          :: wdensmin
  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen      !! B = birth rate of wakes
  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw    !! sigma = fractional area of wakes 
  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens     !! D = number of wakes per unit area 
  REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens    !! A = number of active wakes per unit area
  REAL, DIMENSION (klon),           INTENT(IN)          :: cstar     !! C* = spreading velocity of wakes
  REAL, DIMENSION (klon),           INTENT(IN)          :: cin, wape  ! RM : A Faire disparaitre

!  
  REAL, DIMENSION (klon),           INTENT(OUT)         :: rad_wk    !! r = wake radius
  REAL, DIMENSION (klon),           INTENT(OUT)         :: gfl       !! Lg = gust front lenght per unit area
  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw, d_wdens, d_awdens
  REAL, DIMENSION (klon),           INTENT(OUT)         :: cont_fact  !! RM facteur de contact = 2 pi * rad * C*
  ! Some components of the tendencies of state variables  
  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd


!! internal variables
  
  INTEGER                                               :: i, k
  REAL, DIMENSION (klon)                                :: tau_wk_inv      !! tau = life time of wakes
  REAL                                                  :: tau_wk_inv_min
  REAL, DIMENSION (klon)                                :: tau_prime       !! tau_prime = life time of actives wakes
  REAL                                                  :: d_wdens_targ, d_sigmaw_targ
  

!! Equations
!! dD/dt = B - (D-A)/tau - f D^2
!! dA/dt = B - A/tau_prime + f (D-A)^2 - f A^2
!! dsigma/dt = B a0 - sigma/D (D-A)/tau + Lg C* - f (D-A)^2 (sigma/D-a0)
!!
!! f = 2 (B (a0-sigma/D) + Lg C*) / (2 (D-A)^2 (2 sigma/D-a0) + D (1-2 sigma))


      DO i = 1, klon
        IF (wk_adv(i)) THEN
          rad_wk(i) = max( sqrt(sigmaw(i)/(3.14*wdens(i))) , rzero)
          gfl(i)  = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
        END IF
      END DO


      DO i = 1, klon
        IF (wk_adv(i)) THEN
!!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
          tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
          tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
          tau_prime(i) = tau_cv
!!          cont_fact(i) = 2.*(wgen(i)*(aa0-sigmaw(i)/wdens(i)) + gfl(i)*cstar(i)) / &
!!                             (2.*(wdens(i)-awdens(i))**2*(2.*sigmaw(i)/wdens(i) - aa0) + wdens(i)*(1.-2.*sigmaw(i)))
!!          cont_fact(i) = 2.*3.14*rad_wk(i)*cstar(i)     ! bug
!!          cont_fact(i) = 4.*3.14*rad_wk(i)*cstar(i)
          cont_fact(i) = 2.*gfl(i)*cstar(i)/wdens(i)

          d_sig_gen(i) = wgen(i)*aa0
          d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min
          d_sig_col(i) = - cont_fact(i)*(wdens(i)-awdens(i))**2*(2.*sigmaw(i)/wdens(i)-aa0)
          d_sig_spread(i) = gfl(i)*cstar(i)
!
          d_sig_gen(i) =  d_sig_gen(i)*dtimesub
          d_sig_death(i) = d_sig_death(i)*dtimesub
          d_sig_col(i) =  d_sig_col(i)*dtimesub 
          d_sig_spread(i) =  d_sig_spread(i)*dtimesub 
          d_sigmaw(i) =  d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)

         
          d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
!!          d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
!!          d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
          d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
          d_sigmaw(i) = d_sigmaw_targ
!!          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
         
         
          d_dens_gen(i) = wgen(i)
          d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min 
          d_dens_col(i) =  - cont_fact(i)*wdens(i)**2
! 
          d_dens_gen(i) =  d_dens_gen(i)*dtimesub
          d_dens_death(i) = d_dens_death(i)*dtimesub
          d_dens_col(i) =  d_dens_col(i)*dtimesub 
          d_wdens(i) = d_dens_gen(i) + d_dens_death(i) + d_dens_col(i)


          d_adens_death(i) = -awdens(i)/tau_prime(i)
          d_adens_icol(i) = cont_fact(i)*(wdens(i)-awdens(i))**2
          d_adens_acol(i) = - cont_fact(i)*awdens(i)**2
!
          d_adens_death(i) =  d_adens_death(i)*dtimesub
          d_adens_icol(i) =   d_adens_icol(i)*dtimesub
          d_adens_acol(i) =   d_adens_acol(i)*dtimesub
          d_awdens(i) =   d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i)     
            
!!
          d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
!!          d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
          d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
          d_wdens(i) = d_wdens_targ
          
          d_wdens_targ = min(max(d_awdens(i),-awdens(i)), wdens(i)-awdens(i))
!!          d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
          d_adens_bnd(i) = d_wdens_targ - d_awdens(i)
          d_awdens(i) = d_wdens_targ



        ENDIF
      ENDDO

      IF (prt_level >= 10) THEN
        print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1) ', &
                       cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1)
        print *,'wake, wdens(1), awdens(1), d_awdens(1) ', &
                       wdens(1), awdens(1), d_awdens(1)
        print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
                       d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
      ENDIF
sigmaw=sigmaw+d_sigmaw
wdens=wdens+d_wdens
awdens=awdens+d_awdens

    RETURN 
    END SUBROUTINE wake_popdyn_2  
END MODULE lmdz_wake_popdyn_2
