MODULE lmdz_wake_popdyn_3
PUBLIC wake_popdyn_3
CONTAINS

    SUBROUTINE wake_popdyn_3 ( klon, klev, phys_sub, wk_adv, dtimesub, wgen, &
                             wdensmin, &
                             sigmaw, asigmaw, wdens, awdens, &                       !! state variables
                             gfl, agfl, cstar, cin, wape, &
                             rad_wk, arad_wk, irad_wk, &
                             d_sigmaw, d_asigmaw, d_wdens, d_awdens, &               !! tendencies 
                             d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
                             d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_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 : CPPKEY_IOPHYS_WK
  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
  USE lmdz_wake_ini , ONLY : smallestreal
  
IMPLICIT NONE

  INTEGER, INTENT(IN)                                   :: klon,klev
  LOGICAL,                          INTENT(IN)          :: phys_sub 
  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)       :: asigmaw   !! sigma = fractional area of active 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 = mean wake radius
  REAL, DIMENSION (klon),           INTENT(OUT)         :: arad_wk    !! r_A = wake radius of active wakes
  REAL, DIMENSION (klon),           INTENT(OUT)         :: irad_wk    !! r_I = wake radius of inactive wakes
  REAL, DIMENSION (klon),           INTENT(OUT)         :: gfl       !! Lg = gust front length per unit area
  REAL, DIMENSION (klon),           INTENT(OUT)         :: agfl      !! LgA = gust front length of active wakes
                                                                     !!  per unit area
  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw, d_asigmaw, d_wdens, d_awdens
  ! 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_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_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_acol, d_adens_icol, d_adens_bnd


!! internal variables
  
  INTEGER                                               :: i, k
  REAL, DIMENSION (klon)                                :: iwdens, isigmaw !! inactive wake density and fractional area
!!  REAL, DIMENSION (klon)                                :: d_arad, d_irad
  REAL, DIMENSION (klon)                                :: igfl            !! LgI = gust front length of inactive wakes
                                                                           !!  per unit area
  REAL, DIMENSION (klon)                                :: s_wk            !! mean area of individual wakes
  REAL, DIMENSION (klon)                                :: as_wk           !! mean area of individual active wakes
  REAL, DIMENSION (klon)                                :: is_wk           !! mean area of individual inactive wakes
  REAL, DIMENSION (klon)                                :: tau_wk_inv      !! tau = life time of wakes
  REAL, DIMENSION (klon)                                :: tau_prime       !! tau_prime = life time of actives wakes
  REAL                                                  :: d_wdens_targ, d_sigmaw_targ


!! Equations
!! ---------
!! Gust fronts:
!! Lg_A = 2 pi r_A A
!! Lg_I = 2 pi r_I I
!! Lg   = 2 pi r   D
!!
!! Areas:
!! s = pi r^2
!! s_A = pi r_A^2
!! s_I = pi r_I^2
!!
!! Life expectancy:
!! tau_I = 3 C* ((C*/C*t)^3/2 - 1) / r_I
!!
!! Time deratives:
!! dD/dt = B - (D-A)/tau_I - 2 Lg C* D
!! dA/dt = B - A/tau_A     + 2 Lg_I C* (D-A) - 2 Lg_A C* A
!! dsigma/dt = B a0 - sigma_I/tau_I + Lg C* - 2 Lg_I C* (D-A) (2 s_I - a0)
!! dsigma_A/dt = B a0 - sigma_A/tau_A + Lg_A C* + (Lg_A I + Lg_I A) C* s_I + 2 Lg_I C* I a0
!!

! Initialization
 tau_wk_inv(:) = 0.
! Initialization of output variables
 rad_wk(:) = 0.
 arad_wk(:) = 0.
 irad_wk(:) = 0.
 gfl(:) = 0.
 agfl(:) = 0.
!
 d_wdens(:) = 0.
 d_awdens(:) = 0.
 d_sigmaw(:) = 0.
 d_asigmaw (:) = 0.
!
 d_sig_gen(:) = 0.
 d_sig_death(:) = 0.
 d_sig_col(:) = 0.
 d_sig_spread(:) = 0.
 d_sig_bnd(:) = 0.
 d_asig_death(:) = 0.
 d_asig_aicol(:) = 0.
 d_asig_iicol(:) = 0.
 d_asig_spread(:) = 0.
 d_asig_bnd(:) = 0.
 d_dens_gen(:) = 0.
 d_dens_death(:) = 0.
 d_dens_col(:) = 0.
 d_dens_bnd(:) = 0.
 d_adens_death(:) = 0.
 d_adens_icol(:) = 0.
 d_adens_acol(:) = 0.
 d_adens_bnd(:) = 0.


      DO i = 1, klon
        IF (wk_adv(i)) THEN
         iwdens(i) = wdens(i) - awdens(i)
         isigmaw(i) = sigmaw(i) - asigmaw(i)
!
         arad_wk(i) = max( sqrt(asigmaw(i)/(3.14*awdens(i))) , rzero)
         irad_wk(i) = max( sqrt((sigmaw(i)-asigmaw(i))/  &
                           (3.14*max(smallestreal,(wdens(i)-awdens(i))))), rzero)
         rad_wk(i) = (awdens(i)*arad_wk(i)+(wdens(i)-awdens(i))*irad_wk(i))/wdens(i)
!
         s_wk(i) = 3.14*rad_wk(i)**2
         as_wk(i) = 3.14*arad_wk(i)**2
         is_wk(i) = 3.14*irad_wk(i)**2
!
         gfl(i)  = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
         agfl(i) = 2.*sqrt(3.14*awdens(i)*asigmaw(i))
         igfl(i) = gfl(i) - agfl(i)
        ENDIF
      ENDDO


      DO i = 1, klon
        IF (wk_adv(i)) THEN
!  print *,'ZZZZpopdyn3 wgen(1) ',wgen(1)
!  print *,'ZZZZpopdyn3 cstar(1) ',cstar(1)
!  print *,'ZZZZpopdyn3 isigmaw(1) ',isigmaw(1)
!  print *,'ZZZZpopdyn3 gfl(1) ',gfl(1)
!!          tau_wk_inv(i) = max( (3.*cstar(i))/(irad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
          tau_wk_inv(i) = min(max( (3.*cstar(i))/(irad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.), 1./dtimesub)
          tau_prime(i) = tau_cv

          d_sig_gen(i) = wgen(i)*aa0
          d_sig_death(i) = - isigmaw(i)*tau_wk_inv(i)
          d_sig_col(i) = - 2.*igfl(i)*cstar(i)*iwdens(i)*(2.*is_wk(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_asig_death(i) = - asigmaw(i)/tau_prime(i)
!!  Bug : factor 2 omitted by mistake (bug found by Lamine Thiam)
!!          d_asig_aicol(i) = (agfl(i)*iwdens(i) + igfl(i)*awdens(i))*cstar(i)*is_wk(i)
          d_asig_aicol(i) = 2.*(agfl(i)*iwdens(i) + igfl(i)*awdens(i))*cstar(i)*is_wk(i)
          d_asig_iicol(i) = 2.*igfl(i)*cstar(i)*iwdens(i)*aa0
          d_asig_spread(i) = agfl(i)*cstar(i)
!
          d_asig_death(i) = d_asig_death(i)*dtimesub
          d_asig_aicol(i) =  d_asig_aicol(i)*dtimesub 
          d_asig_iicol(i) =  d_asig_iicol(i)*dtimesub 
          d_asig_spread(i) =  d_asig_spread(i)*dtimesub 
          d_asigmaw(i) =  d_sig_gen(i) + d_asig_death(i) + d_asig_aicol(i) + d_asig_iicol(i) + d_asig_spread(i)
!
          d_sigmaw_targ = min(max(d_asigmaw(i),-asigmaw(i)), sigmaw(i)-asigmaw(i))
!!          d_dens_bnd(i) = d_dens_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
          d_asig_bnd(i) = d_sigmaw_targ - d_asigmaw(i)
          d_asigmaw(i) = d_sigmaw_targ
          d_dens_gen(i) = wgen(i)
          d_dens_death(i) = - iwdens(i)*tau_wk_inv(i)
          d_dens_col(i) =  - 2.*gfl(i)*cstar(i)*wdens(i)
! 
          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_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_adens_death(i) = -awdens(i)/tau_prime(i)
          d_adens_icol(i) =   2.*igfl(i)*cstar(i)*iwdens(i)
          d_adens_acol(i)  = - 2.*agfl(i)*cstar(i)*awdens(i)
!
          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 = 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

!!          d_irad(i) = (d_sigmaw(i)-d_asigmaw(i)-isigmaw(i)*(d_wdens(i)-awdens(i))/iwdens(i)) / &
!!                      max(smallestreal,(2.*3.14*iwdens(i)*irad_wk(i)))
!!          d_arad(i) = (d_asigmaw(i)-asigmaw(i)*d_awdens(i)/awdens(i)) / &
!!                      max(smallestreal,(2.*3.14*awdens(i)*arad_wk(i)))
!!          d_irad(i) = d_irad(i)*dtimesub
!!          d_arad(i) = d_arad(i)*dtimesub
!!          call iophys_ecrit('d_irad',1,'d_irad','m',d_irad)
!!          call iophys_ecrit('d_airad',1,'d_arad','m',d_arad)
!!
        ENDIF
      ENDDO
IF (CPPKEY_IOPHYS_WK) THEN
    IF (phys_sub) THEN
       call iophys_ecrit('d_sigmaw0',1,'d_sigmaw0','',d_sigmaw)
!
       call iophys_ecrit('cstar',1,'cstar','',cstar)
       call iophys_ecrit('wgen_pd3',1,'wgen_popdyn3','',wgen)
       call iophys_ecrit('tauwk_inv',1,'tau_wk_inv','',tau_wk_inv)
       call iophys_ecrit('d_sigmaw',1,'d_sigmaw','',d_sigmaw)
       call iophys_ecrit('d_sig_gen',1,'d_sig_gen','',d_sig_gen)
       call iophys_ecrit('d_sig_death',1,'d_sig_death','',d_sig_death)
       call iophys_ecrit('d_sig_col',1,'d_sig_col','',d_sig_col)
       call iophys_ecrit('d_sig_spread',1,'d_sig_spread','',d_sig_spread)
       call iophys_ecrit('d_sig_bnd',1,'d_sig_bnd','',d_sig_bnd)
!
       call iophys_ecrit('d_asigmaw0',1,'d_asigmaw0','',d_asigmaw)
!
       call iophys_ecrit('d_asigmaw',1,'d_asigmaw','',d_asigmaw)
       call iophys_ecrit('d_asig_death',1,'d_asig_death','',d_asig_death)
       call iophys_ecrit('d_asig_aicol',1,'d_asig_aicol','',d_asig_aicol)
       call iophys_ecrit('d_asig_iicol',1,'d_asig_iicol','',d_asig_iicol)
       call iophys_ecrit('d_asig_spread',1,'d_asig_spread','',d_asig_spread)
       call iophys_ecrit('d_asig_bnd',1,'d_asig_bnd','',d_asig_bnd)
!
       call iophys_ecrit('d_wdens',1,'d_wdens','',d_wdens)
       call iophys_ecrit('d_dens_gen',1,'d_dens_gen','',d_dens_gen)
       call iophys_ecrit('d_dens_death',1,'d_dens_death','',d_dens_death)
       call iophys_ecrit('d_dens_col',1,'d_dens_col','',d_dens_col)
!  
       call iophys_ecrit('d_awdens',1,'d_awdens','',d_awdens)
       call iophys_ecrit('d_adens_death',1,'d_adens_death','',d_adens_death)
       call iophys_ecrit('d_adens_icol',1,'d_adens_icol','',d_adens_icol)
       call iophys_ecrit('d_adens_acol',1,'d_adens_acol','',d_adens_acol)
    ENDIF
END IF


      IF (prt_level >= 10) THEN
        print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), gfl(1) ', &
                       cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), gfl(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
asigmaw=asigmaw+d_asigmaw
wdens=wdens+d_wdens
awdens=awdens+d_awdens

    RETURN 
    END SUBROUTINE wake_popdyn_3  
END MODULE lmdz_wake_popdyn_3
