Ignore:
Timestamp:
Feb 16, 2018, 12:42:18 PM (6 years ago)
Author:
jyg
Message:

Implementation of a first crude model of the
dynamic of wake population.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/wake.F90

    r2922 r3208  
    44SUBROUTINE wake(znatsurf, p, ph, pi, dtime, &
    55                te0, qe0, omgb, &
    6                 dtdwn, dqdwn, amdwn, amup, dta, dqa, &
    7                 sigd_con, &
    8                 deltatw, deltaqw, sigmaw, wdens, &                          ! state variables
     6                dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
     7                sigd_con, Cin, &
     8                deltatw, deltaqw, sigmaw, awdens, wdens, &                  ! state variables
    99                dth, hw, wape, fip, gfl, &
    1010                dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, &
    1111                dtke, dqke, omg, dp_deltomg, spread, cstar, &
    1212                d_deltat_gw, &
    13                 d_deltatw2, d_deltaqw2, d_sigmaw2, d_wdens2)                ! tendencies
     13                d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2)     ! tendencies
    1414
    1515
    … …  
    4848  ! dtls : large scale temperature tendency due to wake
    4949  ! dqls : large scale humidity tendency due to wake
    50   ! hw   : hauteur de la poche
     50  ! hw   : wake top hight (given by hw*deltatw(1)/2=wape)
    5151  ! dp_omgb : vertical gradient of large scale omega
     52  ! awdens  : densite de poches actives
    5253  ! wdens   : densite de poches
    5354  ! omgbdth: flux of Delta_Theta transported by LS omega
    … …  
    7273  ! dta  : source de chaleur due courants satures et detrain  (K/s)
    7374  ! dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
     75  ! wgen : number of wakes generated per unit area and per sec (/m^2/s)
    7476  ! amdwn: flux de masse total des descentes, par unite de
    75   ! surface de la maille (kg/m2/s)
     77  !        surface de la maille (kg/m2/s)
    7678  ! amup : flux de masse total des ascendances, par unite de
    77   ! surface de la maille (kg/m2/s)
     79  !        surface de la maille (kg/m2/s)
     80  ! sigd_con:
     81  ! Cin  : convective inhibition
    7882  ! p    : pressions aux milieux des couches (Pa)
    7983  ! ph   : pressions aux interfaces (Pa)
    … …  
    105109  ! deltatw0   : deltatw initial
    106110  ! deltaqw0   : deltaqw initial
    107   ! hw0    : hw initial
    108   ! sigmaw0: sigmaw initial
     111  ! hw0    : wake top hight (defined as the altitude at which deltatw=0)
    109112  ! amflux : horizontal mass flux through wake boundary
    110113  ! wdens_ref: initial number of wakes per unit area (3D) or per
    … …  
    133136  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: amdwn, amup
    134137  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dta, dqa
     138  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen
    135139  REAL, DIMENSION (klon),           INTENT(IN)          :: sigd_con
     140  REAL, DIMENSION (klon),           INTENT(IN)          :: Cin
    136141
    137142  !
    … …  
    140145  REAL, DIMENSION (klon, klev),     INTENT(INOUT)       :: deltatw, deltaqw
    141146  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw
     147  REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens
    142148  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens
    143149
    … …  
    149155  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls
    150156  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
    151   REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: spread
     157  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: spread    !  unused (jyg)
    152158  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: omgbdth, omg
    153159  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dp_omgb, dp_deltomg
    … …  
    157163  ! Tendencies of state variables
    158164  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
    159   REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_wdens2
     165  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_awdens2, d_wdens2
    160166
    161167  ! Variables internes
    … …  
    165171  INTEGER, SAVE                                         :: igout
    166172  !$OMP THREADPRIVATE(igout)
    167   REAL                                                  :: alon
    168173  LOGICAL, SAVE                                         :: first = .TRUE.
    169174  !$OMP THREADPRIVATE(first)
    … …  
    176181  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
    177182
     183  REAL, SAVE                                            :: tau_cv
     184  !$OMP THREADPRIVATE(tau_cv)
     185  REAL, SAVE                                            :: rzero, aa0 ! minimal wake radius and area
     186  !$OMP THREADPRIVATE(rzero, aa0)
     187
    178188  LOGICAL, SAVE                                         :: flag_wk_check_trgl
    179189  !$OMP THREADPRIVATE(flag_wk_check_trgl)
    180190  INTEGER, SAVE                                         :: iflag_wk_check_trgl
    181191  !$OMP THREADPRIVATE(iflag_wk_check_trgl)
     192  INTEGER, SAVE                                         :: iflag_wk_pop_dyn
     193  !$OMP THREADPRIVATE(iflag_wk_pop_dyn)
    182194
    183195  REAL                                                  :: delta_t_min
    184196  INTEGER                                               :: nsub
    185197  REAL                                                  :: dtimesub
    186   REAL                                                  :: sigmad, hwmin, wapecut
     198  REAL                                                  :: wdensmin
     199  REAL, SAVE                                            :: sigmad, hwmin, wapecut, cstart
     200  !$OMP THREADPRIVATE(sigmad, hwmin, wapecut, cstart)
    187201  REAL                                                  :: sigmaw_max
    188202  REAL                                                  :: dens_rate
    … …  
    195209  REAL, DIMENSION (klon, klev)                          :: deltaqw0
    196210  REAL, DIMENSION (klon, klev)                          :: te, qe
    197   REAL, DIMENSION (klon)                                :: sigmaw0
    198211!!  REAL, DIMENSION (klon)                                :: sigmaw1
     212
     213  ! Variables liees a la dynamique de population
     214  REAL, DIMENSION(klon)                                 :: act
     215  REAL, DIMENSION(klon)                                 :: rad_wk, tau_wk_inv
     216  REAL, DIMENSION(klon)                                 :: f_shear
     217  REAL, DIMENSION(klon)                                 :: drdt
     218  REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col
     219  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
     220  LOGICAL, DIMENSION (klon)                             :: kill_wake
     221  INTEGER, SAVE                                         :: iflag_wk_act
     222  !$OMP THREADPRIVATE(iflag_wk_act)
     223  REAL                                                  :: drdt_pos
     224  REAL                                                  :: tau_wk_inv_min
    199225
    200226  ! Variables pour les GW
    … …  
    204230  REAL, DIMENSION (klon, klev)                          :: tgw
    205231
    206   ! Variables liées au calcul de hw
     232  ! Variables liees au calcul de hw
    207233  REAL, DIMENSION (klon)                                :: ptop_provis, ptop, ptop_new
    208234  REAL, DIMENSION (klon)                                :: sum_dth
    … …  
    211237  INTEGER, DIMENSION (klon)                             :: ktop, kupper
    212238
    213   ! Variables liées au test de la forme triangulaire du profil de Delta_theta
     239  ! Variables liees au test de la forme triangulaire du profil de Delta_theta
    214240  REAL, DIMENSION (klon)                                :: sum_half_dth
    215241  REAL, DIMENSION (klon)                                :: dz_half
    … …  
    218244  REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
    219245  REAL, DIMENSION (klon, klev)                          :: d_te, d_qe
     246  REAL, DIMENSION (klon)                                :: d_awdens, d_wdens
    220247  REAL, DIMENSION (klon)                                :: d_sigmaw, alpha
    221248  REAL, DIMENSION (klon)                                :: q0_min, q1_min
    … …  
    228255  INTEGER                                               ::isubstep, k, i
    229256
     257  REAL                                                  :: wdens_targ
    230258  REAL                                                  :: sigmaw_targ
    231259
    … …  
    273301  REAL, DIMENSION (klon, klev)                          :: detr
    274302
    275   REAL, DIMENSION(klon)                                 :: sigmaw_in   ! pour les prints
     303  REAL, DIMENSION(klon)                                 :: sigmaw_in             ! pour les prints
     304  REAL, DIMENSION(klon)                                 :: awdens_in, wdens_in   ! pour les prints
    276305
    277306  ! -------------------------------------------------------------------------
    … …  
    284313  ! -------------------------------------------------------------------------
    285314
    286   DATA wapecut, sigmad, hwmin/5., .02, 10./
     315!!  DATA wapecut, sigmad, hwmin/5., .02, 10./
     316  DATA wapecut, sigmad, hwmin/1., .02, 10./
     317  DATA wdensmin/1.e-12/
    287318  ! cc nrlmd
    288319  DATA sigmaw_max/0.4/
    289320  DATA dens_rate/0.1/
     321  DATA rzero /5000./
    290322  ! cc
    291323  ! Longueur de maille (en m)
    … …  
    293325
    294326  ! ALON = 3.e5
    295   alon = 1.E6
     327  ! alon = 1.E6
     328
     329  !  Provisionnal; to be suppressed when f_shear is parameterized
     330  f_shear(:) = 1.       ! 0. for strong shear, 1. for weak shear
    296331
    297332
    … …  
    300335  ! coefgw : Coefficient pour les ondes de gravité
    301336  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
    302   ! wdens : Densité de poche froide par maille
     337  ! wdens : Densité surfacique de poche froide
    303338  ! -------------------------------------------------------------------------
    304339
    … …  
    321356  crep_sol = 1.0
    322357
     358  aa0 = 3.14*rzero*rzero
     359
     360  tau_cv = 4000.
     361
    323362  ! cc nrlmd Lecture du fichier wake_param.data
    324363  stark=0.33
    325364  CALL getin_p('stark',stark)
     365  cstart = stark*sqrt(2.*wapecut)
     366
    326367  alpk=0.25
    327368  CALL getin_p('alpk',alpk)
    … …  
    334375  CALL getin_p('wdens_ref_l',wdens_ref(2))    !wake number per unit area ; land
    335376!>jyg
     377  iflag_wk_pop_dyn = 0
     378  CALL getin_p('iflag_wk_pop_dyn',iflag_wk_pop_dyn) ! switch between wdens prescribed
     379                                                    ! and wdens prognostic
     380  iflag_wk_act = 0
     381  CALL getin_p('iflag_wk_act',iflag_wk_act) ! 0: act(:)=0.
     382                                            ! 1: act(:)=1.
     383                                            ! 2: act(:)=f(Wape)
    336384  coefgw=4.
    337385  CALL getin_p('coefgw',coefgw)
    … …  
    344392  WRITE(*,*) 'wdens_ref_l=', wdens_ref(2)
    345393!>jyg
     394  WRITE(*,*) 'iflag_wk_pop_dyn=',iflag_wk_pop_dyn
     395  WRITE(*,*) 'iflag_wk_act',iflag_wk_act
    346396  WRITE(*,*) 'coefgw=', coefgw
    347397
    … …  
    357407 endif
    358408
     409 IF (iflag_wk_pop_dyn == 0) THEN
    359410  ! Initialisation de toutes des densites a wdens_ref.
    360411  ! Les densites peuvent evoluer si les poches debordent
    361412  ! (voir au tout debut de la boucle sur les substeps)
    362 !jyg<
    363 !!  wdens(:) = wdens_ref
    364   DO i = 1,klon
    365     wdens(i) = wdens_ref(znatsurf(i)+1)
    366   ENDDO
    367 !>jyg
     413  !jyg<
     414  !!  wdens(:) = wdens_ref
     415   DO i = 1,klon
     416     wdens(i) = wdens_ref(znatsurf(i)+1)
     417   ENDDO
     418  !>jyg
     419 ENDIF  ! (iflag_wk_pop_dyn == 0)
    368420
    369421  ! print*,'stark',stark
    … …  
    415467      d_deltatw2(:,:) = 0.
    416468      d_deltaqw2(:,:) = 0.
     469
     470      IF (iflag_wk_act == 0) THEN
     471        act(:) = 0.
     472      ELSEIF (iflag_wk_act == 1) THEN
     473        act(:) = 1.
     474      ENDIF
     475
    417476!!  DO i = 1, klon
    418477!!   sigmaw_in(i) = sigmaw(i)
    … …  
    425484  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
    426485  ! ENDIF
     486  IF (iflag_wk_pop_dyn >=1) THEN
     487    DO i = 1, klon
     488      wdens_targ = max(wdens(i),wdensmin)
     489      d_wdens2(i) = wdens_targ - wdens(i)
     490      wdens(i) = wdens_targ
     491    END DO
     492  ELSE
     493    DO i = 1, klon
     494      d_awdens2(i) = 0.
     495      d_wdens2(i) = 0.
     496    END DO
     497  ENDIF  ! (iflag_wk_pop_dyn >=1)
     498!
    427499  DO i = 1, klon
    428500    ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
    … …  
    434506    sigmaw(i) = sigmaw_targ
    435507!>jyg
    436     sigmaw0(i) = sigmaw(i)
    437     wape(i) = 0.
    438     wape2(i) = 0.
    439     d_sigmaw(i) = 0.
    440     d_wdens2(i) = 0.
    441     ktopw(i) = 0
    442   END DO
     508  END DO
     509
     510!
     511  IF (iflag_wk_pop_dyn >= 1) THEN
     512    awdens_in(:) = awdens(:)
     513    wdens_in(:) = wdens(:)
     514!!    wdens(:) = wdens(:) + wgen(:)*dtime
     515!!    d_wdens2(:) = wgen(:)*dtime
     516!!  ELSE
     517  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     518
     519  wape(:) = 0.
     520  wape2(:) = 0.
     521  d_sigmaw(:) = 0.
     522  ktopw(:) = 0
    443523!
    444524!<jyg
    … …  
    833913      gwake(i) = .FALSE.
    834914    ELSE
     915      hw(i) = hw0(i)
    835916      cstar(i) = stark*sqrt(2.*wape(i))
    836917      gwake(i) = .TRUE.
    … …  
    891972    ! cc           On calcule pour cela une densité wdens0 pour laquelle on
    892973    ! aurait un entrainement nul ---
     974    !jyg<
     975    ! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou
     976    ! les poches sont insuffisantes pour accueillir tout le flux de masse
     977    ! des descentes unsaturees. Nous faisons alors l'hypothese que la
     978    ! convection profonde cree directement de nouvelles poches, sans passer
     979    ! par les thermiques. La nouvelle valeur de wdens est alors imposée.
     980
    893981    DO i = 1, klon
    894982      ! c       print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
    … …  
    900988          ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2)
    901989        IF (wdens(i)<=wdens0*1.1) THEN
     990          IF (iflag_wk_pop_dyn >= 1) THEN
     991             d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
     992          ENDIF
    902993          wdens(i) = wdens0
    903994        END IF
    … …  
    9091000    END DO
    9101001
    911     ! cc nrlmd
    912 
    9131002    DO i = 1, klon
    9141003      IF (wk_adv(i)) THEN
    9151004        gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     1005        rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
    9161006!jyg<
    9171007!!        sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
    … …  
    9231013    END DO
    9241014
    925     DO i = 1, klon
    926       IF (wk_adv(i)) THEN
    927         ! cc nrlmd          Introduction du taux de mortalité des poches et
    928         ! test sur sigmaw_max=0.4
    929         ! cc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
    930         IF (sigmaw(i)>=sigmaw_max) THEN
    931           death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
    932         ELSE
    933           death_rate(i) = 0.
    934         END IF
    935 
    936         d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* &
    937           dtimesub
    938         ! $              - nat_rate(i)*sigmaw(i)*dtimesub
    939         ! c        print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
    940         ! c     $  death_rate(i),ktop(i),kupper(i)',
    941         ! c     $                d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
    942         ! c     $  death_rate(i),ktop(i),kupper(i)
    943 
    944         ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
    945         ! sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
    946         ! wdens = wdens0/(10.*sigmaw)
    947         ! sigmaw =max(sigmaw,sigd_con)
    948         ! sigmaw =max(sigmaw,sigmad)
    949       END IF
    950     END DO
     1015    IF (iflag_wk_pop_dyn >= 1) THEN
     1016
     1017      IF (iflag_wk_act ==2) THEN
     1018      DO i = 1, klon
     1019        IF (wk_adv(i)) THEN
     1020          wape1_act(i) = abs(cin(i))
     1021          wape2_act(i) = 2.*wape1_act(i) + 1.
     1022          act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
     1023        ENDIF  ! (wk_adv(i))
     1024      ENDDO
     1025      ENDIF  ! (iflag_wk_act ==2)
     1026
     1027
     1028      DO i = 1, klon
     1029        IF (wk_adv(i)) THEN
     1030!!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
     1031          tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
     1032          tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
     1033          drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
     1034                    (1 - 2*sigmaw(i)*(1.-f_shear(i)))
     1035          drdt_pos=max(drdt(i),0.)
     1036
     1037!!          d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) &
     1038!!                     - wdens(i)*tau_wk_inv_min &
     1039!!                     - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub
     1040          d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub
     1041          d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min -  &
     1042                         2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub
     1043          d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
     1044
     1045!!          d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) &
     1046!!                      + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) &
     1047!!                      - sigmaw(i)*tau_wk_inv_min )*dtimesub
     1048          d_sig_gen(i) = wgen(i)*aa0
     1049          d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min
     1050!!          d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos
     1051          d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos
     1052          d_sigmaw(i) = ( d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + gfl(i)*cstar(i) )*dtimesub
     1053          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
     1054        ENDIF
     1055      ENDDO
     1056
     1057      IF (prt_level >= 10) THEN
     1058        print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
     1059                       cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
     1060        print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
     1061                       wdens(1), awdens(1), act(1), d_awdens(1)
     1062        print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', &
     1063                       wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1)
     1064        print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
     1065                       d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
     1066      ENDIF
     1067   
     1068    ELSE  !  (iflag_wk_pop_dyn >= 1)
     1069
     1070    ! cc nrlmd
     1071
     1072      DO i = 1, klon
     1073        IF (wk_adv(i)) THEN
     1074          ! cc nrlmd          Introduction du taux de mortalité des poches et
     1075          ! test sur sigmaw_max=0.4
     1076          ! cc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
     1077          IF (sigmaw(i)>=sigmaw_max) THEN
     1078            death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
     1079          ELSE
     1080            death_rate(i) = 0.
     1081          END IF
     1082   
     1083          d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* &
     1084            dtimesub
     1085          ! $              - nat_rate(i)*sigmaw(i)*dtimesub
     1086          ! c        print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
     1087          ! c     $  death_rate(i),ktop(i),kupper(i)',
     1088          ! c     $              d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
     1089          ! c     $  death_rate(i),ktop(i),kupper(i)
     1090   
     1091          ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
     1092          ! sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
     1093          ! wdens = wdens0/(10.*sigmaw)
     1094          ! sigmaw =max(sigmaw,sigd_con)
     1095          ! sigmaw =max(sigmaw,sigmad)
     1096        END IF
     1097      END DO
     1098
     1099    ENDIF   !  (iflag_wk_pop_dyn >= 1)
     1100
    9511101
    9521102    ! calcul de la difference de vitesse verticale poche - zone non perturbee
    … …  
    12231373
    12241374    ! Increment state variables
     1375!jyg<
     1376    IF (iflag_wk_pop_dyn >= 1) THEN
     1377      DO k = 1, klev
     1378        DO i = 1, klon
     1379          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     1380            detr(i,k) = - d_sig_death(i) - d_sig_col(i)     
     1381            entr(i,k) = d_sig_gen(i)
     1382          ENDIF
     1383        ENDDO
     1384      ENDDO
     1385      ELSE  ! (iflag_wk_pop_dyn >= 1)
     1386      DO k = 1, klev
     1387        DO i = 1, klon
     1388          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     1389            detr(i, k) = 0.
     1390   
     1391            entr(i, k) = 0.
     1392          ENDIF
     1393        ENDDO
     1394      ENDDO
     1395    ENDIF  ! (iflag_wk_pop_dyn >= 1)
     1396
     1397   
    12251398
    12261399    DO k = 1, klev
    … …  
    12641437          ! cc nrlmd          Prise en compte du taux de mortalité
    12651438          ! cc               Définitions de entr, detr
    1266           detr(i, k) = 0.
    1267 
    1268           entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
    1269             sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)
    1270 
    1271           spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
     1439!jyg<
     1440!!            detr(i, k) = 0.
     1441!!   
     1442!!            entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
     1443!!              sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)
     1444!!
     1445            entr(i, k) = entr(i,k) + gfl(i)*cstar(i) + &
     1446                         sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)   
     1447!>jyg
     1448            spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
     1449
    12721450          ! cc        spread(i,k) =
    12731451          ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
    … …  
    13841562      END DO
    13851563    END DO
     1564!
    13861565    DO i = 1, klon
    13871566      IF (wk_adv(i)) THEN
    13881567        sigmaw(i) = sigmaw(i) + d_sigmaw(i)
     1568        d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
     1569      END IF
     1570    END DO
    13891571!jyg<
    1390         d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
     1572    IF (iflag_wk_pop_dyn >= 1) THEN
     1573      DO i = 1, klon
     1574        IF (wk_adv(i)) THEN
     1575          awdens(i) = awdens(i) + d_awdens(i)
     1576          wdens(i) = wdens(i) + d_wdens(i)
     1577          d_awdens2(i) = d_awdens2(i) + d_awdens(i)
     1578          d_wdens2(i) = d_wdens2(i) + d_wdens(i)
     1579        END IF
     1580      END DO
     1581      DO i = 1, klon
     1582        IF (wk_adv(i)) THEN
     1583          wdens_targ = max(wdens(i),wdensmin)
     1584          d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
     1585          wdens(i) = wdens_targ
     1586!
     1587          wdens_targ = min( max(awdens(i),0.), wdens(i) )
     1588          d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
     1589          awdens(i) = wdens_targ
     1590        END IF
     1591      END DO
     1592      DO i = 1, klon
     1593        IF (wk_adv(i)) THEN
     1594          sigmaw_targ = max(sigmaw(i),sigmad)
     1595          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     1596          sigmaw(i) = sigmaw_targ
     1597        END IF
     1598      END DO
     1599    ENDIF  ! (iflag_wk_pop_dyn >= 1)
    13911600!>jyg
    1392       END IF
    1393     END DO
    13941601
    13951602
    … …  
    19012108  ! ENDDO
    19022109  ! cc
     2110
     2111  !jyg<
     2112  IF (iflag_wk_pop_dyn >= 1) THEN
     2113    DO i = 1, klon
     2114      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     2115          .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
     2116    ENDDO
     2117  ELSE  ! (iflag_wk_pop_dyn >= 1)
     2118    DO i = 1, klon
     2119      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     2120          .NOT. ok_qx_qw(i)
     2121    ENDDO
     2122  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     2123  !>jyg
     2124
    19032125  DO k = 1, klev
    19042126    DO i = 1, klon
    1905 
    1906       ! cc nrlmd      On maintient désormais constant sigmaw en régime
    1907       ! permanent
    1908       ! cc      IF ((sigmaw(i).GT.sigmaw_max).or.
    1909       IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
    1910           .NOT. ok_qx_qw(i)) THEN
     2127!!jyg      IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     2128!!jyg          .NOT. ok_qx_qw(i)) THEN
     2129      IF (kill_wake(i)) THEN
    19112130        ! cc
    19122131        dtls(i, k) = 0.
    … …  
    19162135        d_deltatw2(i,k) = -deltatw0(i,k)
    19172136        d_deltaqw2(i,k) = -deltaqw0(i,k)
    1918       END IF
    1919     END DO
    1920   END DO
    1921 
    1922   ! cc nrlmd      On maintient désormais constant sigmaw en régime permanent
    1923   DO i = 1, klon
    1924     IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
    1925         .NOT. ok_qx_qw(i)) THEN
     2137      END IF  ! (kill_wake(i))
     2138    END DO
     2139  END DO
     2140
     2141  DO i = 1, klon
     2142!!jyg    IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     2143!!jyg        .NOT. ok_qx_qw(i)) THEN
     2144      IF (kill_wake(i)) THEN
    19262145      ktopw(i) = 0
    19272146      wape(i) = 0.
    19282147      cstar(i) = 0.
    1929 !!jyg   Outside subroutine "Wake" hw and sigmaw are zero when there are no wakes
     2148!!jyg   Outside subroutine "Wake" hw, wdens and sigmaw are zero when there are no wakes
    19302149!!      hw(i) = hwmin                       !jyg
    19312150!!      sigmaw(i) = sigmad                  !jyg
    19322151      hw(i) = 0.                            !jyg
    1933       sigmaw(i) = 0.                        !jyg
    19342152      fip(i) = 0.
    1935     ELSE
     2153!!      sigmaw(i) = 0.                        !jyg
     2154      sigmaw_targ = 0.
     2155      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     2156      sigmaw(i) = sigmaw_targ
     2157      IF (iflag_wk_pop_dyn >= 1) THEN
     2158!!        awdens(i) = 0.
     2159!!        wdens(i) = 0.
     2160        wdens_targ = 0.
     2161        d_wdens2(i) = wdens_targ - wdens(i)
     2162        wdens(i) = wdens_targ
     2163        wdens_targ = 0.
     2164        d_awdens2(i) = wdens_targ - awdens(i)
     2165        awdens(i) = wdens_targ
     2166      ENDIF  ! (iflag_wk_pop_dyn >= 1)
     2167    ELSE  ! (kill_wake(i))
    19362168      wape(i) = wape2(i)
    19372169      cstar(i) = cstar2(i)
    1938     END IF
     2170    END IF  ! (kill_wake(i))
    19392171    ! c        print*,'wape wape2 ktopw OK_qx_qw =',
    19402172    ! c     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
    … …  
    19722204  DO i = 1, klon
    19732205    d_sigmaw2(i) = d_sigmaw2(i)/dtime
     2206    d_awdens2(i) = d_awdens2(i)/dtime
    19742207    d_wdens2(i) = d_wdens2(i)/dtime
    19752208  ENDDO
Note: See TracChangeset for help on using the changeset viewer.