Ignore:
Timestamp:
Jun 29, 2018, 12:31:11 PM (7 years ago)
Author:
Laurent Fairhead
Message:

First attempt at merging with trunk

Location:
LMDZ6/branches/DYNAMICO-conv
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/DYNAMICO-conv

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/wake.F90

    r2922 r3356  
    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/
     318  DATA wdensmin/1.e-14/
    287319  ! cc nrlmd
    288320  DATA sigmaw_max/0.4/
    289321  DATA dens_rate/0.1/
     322  DATA rzero /5000./
    290323  ! cc
    291324  ! Longueur de maille (en m)
    … …  
    293326
    294327  ! ALON = 3.e5
    295   alon = 1.E6
     328  ! alon = 1.E6
     329
     330  !  Provisionnal; to be suppressed when f_shear is parameterized
     331  f_shear(:) = 1.       ! 0. for strong shear, 1. for weak shear
    296332
    297333
    … …  
    300336  ! coefgw : Coefficient pour les ondes de gravité
    301337  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
    302   ! wdens : Densité de poche froide par maille
     338  ! wdens : Densité surfacique de poche froide
    303339  ! -------------------------------------------------------------------------
    304340
    … …  
    321357  crep_sol = 1.0
    322358
     359  aa0 = 3.14*rzero*rzero
     360
     361  tau_cv = 4000.
     362
    323363  ! cc nrlmd Lecture du fichier wake_param.data
    324364  stark=0.33
    325365  CALL getin_p('stark',stark)
     366  cstart = stark*sqrt(2.*wapecut)
     367
    326368  alpk=0.25
    327369  CALL getin_p('alpk',alpk)
    … …  
    334376  CALL getin_p('wdens_ref_l',wdens_ref(2))    !wake number per unit area ; land
    335377!>jyg
     378  iflag_wk_pop_dyn = 0
     379  CALL getin_p('iflag_wk_pop_dyn',iflag_wk_pop_dyn) ! switch between wdens prescribed
     380                                                    ! and wdens prognostic
     381  iflag_wk_act = 0
     382  CALL getin_p('iflag_wk_act',iflag_wk_act) ! 0: act(:)=0.
     383                                            ! 1: act(:)=1.
     384                                            ! 2: act(:)=f(Wape)
    336385  coefgw=4.
    337386  CALL getin_p('coefgw',coefgw)
    … …  
    344393  WRITE(*,*) 'wdens_ref_l=', wdens_ref(2)
    345394!>jyg
     395  WRITE(*,*) 'iflag_wk_pop_dyn=',iflag_wk_pop_dyn
     396  WRITE(*,*) 'iflag_wk_act',iflag_wk_act
    346397  WRITE(*,*) 'coefgw=', coefgw
    347398
    … …  
    357408 endif
    358409
     410 IF (iflag_wk_pop_dyn == 0) THEN
    359411  ! Initialisation de toutes des densites a wdens_ref.
    360412  ! Les densites peuvent evoluer si les poches debordent
    361413  ! (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
     414  !jyg<
     415  !!  wdens(:) = wdens_ref
     416   DO i = 1,klon
     417     wdens(i) = wdens_ref(znatsurf(i)+1)
     418   ENDDO
     419  !>jyg
     420 ENDIF  ! (iflag_wk_pop_dyn == 0)
    368421
    369422  ! print*,'stark',stark
    … …  
    415468      d_deltatw2(:,:) = 0.
    416469      d_deltaqw2(:,:) = 0.
     470
     471      IF (iflag_wk_act == 0) THEN
     472        act(:) = 0.
     473      ELSEIF (iflag_wk_act == 1) THEN
     474        act(:) = 1.
     475      ENDIF
     476
    417477!!  DO i = 1, klon
    418478!!   sigmaw_in(i) = sigmaw(i)
    … …  
    425485  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
    426486  ! ENDIF
     487  IF (iflag_wk_pop_dyn >=1) THEN
     488    DO i = 1, klon
     489      wdens_targ = max(wdens(i),wdensmin)
     490      d_wdens2(i) = wdens_targ - wdens(i)
     491      wdens(i) = wdens_targ
     492    END DO
     493  ELSE
     494    DO i = 1, klon
     495      d_awdens2(i) = 0.
     496      d_wdens2(i) = 0.
     497    END DO
     498  ENDIF  ! (iflag_wk_pop_dyn >=1)
     499!
    427500  DO i = 1, klon
    428501    ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
    … …  
    434507    sigmaw(i) = sigmaw_targ
    435508!>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
     509  END DO
     510
     511!
     512  IF (iflag_wk_pop_dyn >= 1) THEN
     513    awdens_in(:) = awdens(:)
     514    wdens_in(:) = wdens(:)
     515!!    wdens(:) = wdens(:) + wgen(:)*dtime
     516!!    d_wdens2(:) = wgen(:)*dtime
     517!!  ELSE
     518  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     519
     520  wape(:) = 0.
     521  wape2(:) = 0.
     522  d_sigmaw(:) = 0.
     523  ktopw(:) = 0
    443524!
    444525!<jyg
    … …  
    833914      gwake(i) = .FALSE.
    834915    ELSE
     916      hw(i) = hw0(i)
    835917      cstar(i) = stark*sqrt(2.*wape(i))
    836918      gwake(i) = .TRUE.
    … …  
    891973    ! cc           On calcule pour cela une densité wdens0 pour laquelle on
    892974    ! aurait un entrainement nul ---
     975    !jyg<
     976    ! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou
     977    ! les poches sont insuffisantes pour accueillir tout le flux de masse
     978    ! des descentes unsaturees. Nous faisons alors l'hypothese que la
     979    ! convection profonde cree directement de nouvelles poches, sans passer
     980    ! par les thermiques. La nouvelle valeur de wdens est alors imposée.
     981
    893982    DO i = 1, klon
    894983      ! c       print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
    … …  
    899988        wdens0 = (sigmaw(i)/(4.*3.14))* &
    900989          ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2)
     990        IF (prt_level >= 10) THEN
     991             print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', &
     992                     omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)
     993        ENDIF
    901994        IF (wdens(i)<=wdens0*1.1) THEN
     995          IF (iflag_wk_pop_dyn >= 1) THEN
     996             d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
     997          ENDIF
    902998          wdens(i) = wdens0
    903999        END IF
    904         ! c        print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
    905         ! c     $     ,ph(i,1)-pupper(i)',
    906         ! c     $             omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
    907         ! c     $     ,ph(i,1)-pupper(i)
    908       END IF
    909     END DO
    910 
    911     ! cc nrlmd
     1000      END IF
     1001    END DO
    9121002
    9131003    DO i = 1, klon
    9141004      IF (wk_adv(i)) THEN
    9151005        gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     1006        rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
    9161007!jyg<
    9171008!!        sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
    … …  
    9231014    END DO
    9241015
    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
     1016    IF (iflag_wk_pop_dyn >= 1) THEN
     1017
     1018      IF (iflag_wk_act ==2) THEN
     1019      DO i = 1, klon
     1020        IF (wk_adv(i)) THEN
     1021          wape1_act(i) = abs(cin(i))
     1022          wape2_act(i) = 2.*wape1_act(i) + 1.
     1023          act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
     1024        ENDIF  ! (wk_adv(i))
     1025      ENDDO
     1026      ENDIF  ! (iflag_wk_act ==2)
     1027
     1028
     1029      DO i = 1, klon
     1030        IF (wk_adv(i)) THEN
     1031!!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
     1032          tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
     1033          tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
     1034          drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
     1035                    (1 + 2*f_shear(i)*(2.*sigmaw(i)-aa0*wdens(i)) - 2.*sigmaw(i))
     1036!!                    (1 - 2*sigmaw(i)*(1.-f_shear(i)))
     1037          drdt_pos=max(drdt(i),0.)
     1038
     1039!!          d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) &
     1040!!                     - wdens(i)*tau_wk_inv_min &
     1041!!                     - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub
     1042          d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub
     1043          d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min -  &
     1044                         2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub
     1045          d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
     1046
     1047!!          d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) &
     1048!!                      + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) &
     1049!!                      - sigmaw(i)*tau_wk_inv_min )*dtimesub
     1050          d_sig_gen(i) = wgen(i)*aa0
     1051          d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min
     1052!!          d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos
     1053          d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos
     1054          d_sigmaw(i) = ( d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + gfl(i)*cstar(i) )*dtimesub
     1055          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
     1056        ENDIF
     1057      ENDDO
     1058
     1059      IF (prt_level >= 10) THEN
     1060        print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
     1061                       cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
     1062        print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
     1063                       wdens(1), awdens(1), act(1), d_awdens(1)
     1064        print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', &
     1065                       wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1)
     1066        print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
     1067                       d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
     1068      ENDIF
     1069   
     1070    ELSE  !  (iflag_wk_pop_dyn >= 1)
     1071
     1072    ! cc nrlmd
     1073
     1074      DO i = 1, klon
     1075        IF (wk_adv(i)) THEN
     1076          ! cc nrlmd          Introduction du taux de mortalité des poches et
     1077          ! test sur sigmaw_max=0.4
     1078          ! cc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
     1079          IF (sigmaw(i)>=sigmaw_max) THEN
     1080            death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
     1081          ELSE
     1082            death_rate(i) = 0.
     1083          END IF
     1084   
     1085          d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* &
     1086            dtimesub
     1087          ! $              - nat_rate(i)*sigmaw(i)*dtimesub
     1088          ! c        print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
     1089          ! c     $  death_rate(i),ktop(i),kupper(i)',
     1090          ! c     $              d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
     1091          ! c     $  death_rate(i),ktop(i),kupper(i)
     1092   
     1093          ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
     1094          ! sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
     1095          ! wdens = wdens0/(10.*sigmaw)
     1096          ! sigmaw =max(sigmaw,sigd_con)
     1097          ! sigmaw =max(sigmaw,sigmad)
     1098        END IF
     1099      END DO
     1100
     1101    ENDIF   !  (iflag_wk_pop_dyn >= 1)
     1102
    9511103
    9521104    ! calcul de la difference de vitesse verticale poche - zone non perturbee
    … …  
    12231375
    12241376    ! Increment state variables
     1377!jyg<
     1378    IF (iflag_wk_pop_dyn >= 1) THEN
     1379      DO k = 1, klev
     1380        DO i = 1, klon
     1381          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     1382            detr(i,k) = - d_sig_death(i) - d_sig_col(i)     
     1383            entr(i,k) = d_sig_gen(i)
     1384          ENDIF
     1385        ENDDO
     1386      ENDDO
     1387      ELSE  ! (iflag_wk_pop_dyn >= 1)
     1388      DO k = 1, klev
     1389        DO i = 1, klon
     1390          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     1391            detr(i, k) = 0.
     1392   
     1393            entr(i, k) = 0.
     1394          ENDIF
     1395        ENDDO
     1396      ENDDO
     1397    ENDIF  ! (iflag_wk_pop_dyn >= 1)
     1398
     1399   
    12251400
    12261401    DO k = 1, klev
    … …  
    12641439          ! cc nrlmd          Prise en compte du taux de mortalité
    12651440          ! 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)
     1441!jyg<
     1442!!            detr(i, k) = 0.
     1443!!   
     1444!!            entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
     1445!!              sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)
     1446!!
     1447            entr(i, k) = entr(i,k) + gfl(i)*cstar(i) + &
     1448                         sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)   
     1449!>jyg
     1450            spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
     1451
    12721452          ! cc        spread(i,k) =
    12731453          ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
    … …  
    13841564      END DO
    13851565    END DO
     1566!
    13861567    DO i = 1, klon
    13871568      IF (wk_adv(i)) THEN
    13881569        sigmaw(i) = sigmaw(i) + d_sigmaw(i)
     1570        d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
     1571      END IF
     1572    END DO
    13891573!jyg<
    1390         d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
     1574    IF (iflag_wk_pop_dyn >= 1) THEN
     1575      DO i = 1, klon
     1576        IF (wk_adv(i)) THEN
     1577          awdens(i) = awdens(i) + d_awdens(i)
     1578          wdens(i) = wdens(i) + d_wdens(i)
     1579          d_awdens2(i) = d_awdens2(i) + d_awdens(i)
     1580          d_wdens2(i) = d_wdens2(i) + d_wdens(i)
     1581        END IF
     1582      END DO
     1583      DO i = 1, klon
     1584        IF (wk_adv(i)) THEN
     1585          wdens_targ = max(wdens(i),wdensmin)
     1586          d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
     1587          wdens(i) = wdens_targ
     1588!
     1589          wdens_targ = min( max(awdens(i),0.), wdens(i) )
     1590          d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
     1591          awdens(i) = wdens_targ
     1592        END IF
     1593      END DO
     1594      DO i = 1, klon
     1595        IF (wk_adv(i)) THEN
     1596          sigmaw_targ = max(sigmaw(i),sigmad)
     1597          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     1598          sigmaw(i) = sigmaw_targ
     1599        END IF
     1600      END DO
     1601    ENDIF  ! (iflag_wk_pop_dyn >= 1)
    13911602!>jyg
    1392       END IF
    1393     END DO
    13941603
    13951604
    … …  
    19012110  ! ENDDO
    19022111  ! cc
     2112
     2113  !jyg<
     2114  IF (iflag_wk_pop_dyn >= 1) THEN
     2115    DO i = 1, klon
     2116      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     2117          .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
     2118    ENDDO
     2119  ELSE  ! (iflag_wk_pop_dyn >= 1)
     2120    DO i = 1, klon
     2121      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     2122          .NOT. ok_qx_qw(i)
     2123    ENDDO
     2124  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     2125  !>jyg
     2126
    19032127  DO k = 1, klev
    19042128    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
     2129!!jyg      IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     2130!!jyg          .NOT. ok_qx_qw(i)) THEN
     2131      IF (kill_wake(i)) THEN
    19112132        ! cc
    19122133        dtls(i, k) = 0.
    … …  
    19162137        d_deltatw2(i,k) = -deltatw0(i,k)
    19172138        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
     2139      END IF  ! (kill_wake(i))
     2140    END DO
     2141  END DO
     2142
     2143  DO i = 1, klon
     2144!!jyg    IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     2145!!jyg        .NOT. ok_qx_qw(i)) THEN
     2146      IF (kill_wake(i)) THEN
    19262147      ktopw(i) = 0
    19272148      wape(i) = 0.
    19282149      cstar(i) = 0.
    1929 !!jyg   Outside subroutine "Wake" hw and sigmaw are zero when there are no wakes
     2150!!jyg   Outside subroutine "Wake" hw, wdens and sigmaw are zero when there are no wakes
    19302151!!      hw(i) = hwmin                       !jyg
    19312152!!      sigmaw(i) = sigmad                  !jyg
    19322153      hw(i) = 0.                            !jyg
    1933       sigmaw(i) = 0.                        !jyg
    19342154      fip(i) = 0.
    1935     ELSE
     2155!!      sigmaw(i) = 0.                        !jyg
     2156      sigmaw_targ = 0.
     2157      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     2158      sigmaw(i) = sigmaw_targ
     2159      IF (iflag_wk_pop_dyn >= 1) THEN
     2160!!        awdens(i) = 0.
     2161!!        wdens(i) = 0.
     2162        wdens_targ = 0.
     2163        d_wdens2(i) = wdens_targ - wdens(i)
     2164        wdens(i) = wdens_targ
     2165        wdens_targ = 0.
     2166        d_awdens2(i) = wdens_targ - awdens(i)
     2167        awdens(i) = wdens_targ
     2168      ENDIF  ! (iflag_wk_pop_dyn >= 1)
     2169    ELSE  ! (kill_wake(i))
    19362170      wape(i) = wape2(i)
    19372171      cstar(i) = cstar2(i)
    1938     END IF
     2172    END IF  ! (kill_wake(i))
    19392173    ! c        print*,'wape wape2 ktopw OK_qx_qw =',
    19402174    ! c     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
    … …  
    19722206  DO i = 1, klon
    19732207    d_sigmaw2(i) = d_sigmaw2(i)/dtime
     2208    d_awdens2(i) = d_awdens2(i)/dtime
    19742209    d_wdens2(i) = d_wdens2(i)/dtime
    19752210  ENDDO
Note: See TracChangeset for help on using the changeset viewer.