MODULE lmdz_wake_popdyn_1
PUBLIC wake_popdyn_1
CONTAINS

SUBROUTINE wake_popdyn_1(klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
                  wdensmin, &
                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
                  d_awdens, d_wdens, d_sigmaw, &
                  iflag_wk_act, wk_adv, cin, wape, &
                  drdt, &
                  d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
                  d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
                  d_wdens_targ, d_sigmaw_targ)
                

  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)          :: dtime
  REAL,                             INTENT(IN)          :: dtimesub
  REAL,                             INTENT(IN)          :: wdensmin
  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen
  REAL, DIMENSION (klon),           INTENT(IN)          :: wdens
  REAL, DIMENSION (klon),           INTENT(IN)          :: awdens
  REAL, DIMENSION (klon),           INTENT(IN)          :: sigmaw
  REAL, DIMENSION (klon),           INTENT(IN)          :: cstar
  REAL, DIMENSION (klon),           INTENT(IN)          :: cin, wape
  REAL, DIMENSION (klon),           INTENT(IN)          :: f_shear
  INTEGER,                          INTENT(IN)          :: iflag_wk_act

  
  !
 
  ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields 
  !                                 computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
  REAL, DIMENSION (klon),           INTENT(OUT)         :: rad_wk
  REAL, DIMENSION (klon),           INTENT(OUT)         :: gfl
  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw, d_awdens, d_wdens
  REAL, DIMENSION (klon),           INTENT(OUT)         :: drdt
  ! Some components of the tendencies of state variables  
  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_gen, d_sig_death, d_sig_col, d_sig_bnd
  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_spread
  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
  REAL,                             INTENT(OUT)         :: d_wdens_targ, d_sigmaw_targ
 
  
  REAL                                                  :: delta_t_min
  INTEGER                                               :: i, k
  REAL                                                  :: wdens0
  ! IM 080208
  LOGICAL, DIMENSION (klon)                             :: gwake
  
   ! Variables liees a la dynamique de population
  REAL, DIMENSION(klon)                                 :: act
  REAL, DIMENSION(klon)                                 :: tau_wk_inv
  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
  LOGICAL, DIMENSION (klon)                             :: kill_wake
  REAL                                                  :: drdt_pos
  REAL                                                  :: tau_wk_inv_min
  
     

      IF (iflag_wk_act == 0) THEN
        act(:) = 0.
      ELSEIF (iflag_wk_act == 1) THEN
        act(:) = 1.
      ELSEIF (iflag_wk_act ==2) THEN
      DO i = 1, klon
        IF (wk_adv(i)) THEN
          wape1_act(i) = abs(cin(i))
          wape2_act(i) = 2.*wape1_act(i) + 1.
          act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
        ENDIF  ! (wk_adv(i))
      ENDDO
      ENDIF  ! (iflag_wk_act ==2)

      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)
          drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
                    (1 + 2*f_shear(i)*(2.*sigmaw(i)-aa0*wdens(i)) - 2.*sigmaw(i))
!!                    (1 - 2*sigmaw(i)*(1.-f_shear(i)))
          drdt_pos=max(drdt(i),0.)

!!          d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) &
!!                     - wdens(i)*tau_wk_inv_min &
!!                     - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub
!jyg+mlt<
          d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub
          d_dens_gen(i) = wgen(i)
          d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min 
          d_dens_col(i) =  -2.*wdens(i)*gfl(i)*drdt_pos 
          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(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min -  &
!!                         2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub
!>jyg+mlt
!
!jyg<
          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(i) = max(d_wdens(i), wdensmin-wdens(i))
!>jyg

!jyg+mlt<
!!          d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) &
!!                      + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) &
!!                      - sigmaw(i)*tau_wk_inv_min )*dtimesub
          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) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos
          d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos
          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)
!>jyg+mlt
!
!jyg<
          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))
!>jyg
        ENDIF
      ENDDO

      IF (prt_level >= 10) THEN
        print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
                       cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
        print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
                       wdens(1), awdens(1), act(1), d_awdens(1)
        print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', &
                       wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(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
    
    RETURN 
    END SUBROUTINE wake_popdyn_1
END MODULE lmdz_wake_popdyn_1
