Changeset 5763


Ignore:
Timestamp:
Jul 8, 2025, 11:09:14 AM (10 hours ago)
Author:
jyg
Message:

Sequel of commit 5761: two versions of subroutine wake are implemented in the module
lmdz_wake.f90, one with the old computations for the vertical differential advection
(called "wake") and the other with the new computations introduced in commit 5761
(called wake2).

The switch between the two versions is controled by the flag iflag_wake.
iflag_wake is now made of two digits: the first one points to the version of wake

to be used (0 for wake or 2 for wake2); the second one defines the differential convective
heating.

+ iflag_wake= 01, 02, 03 ==> as before
+ iflag_wake= 21, 22, 23 ==> use wake2 instead of wake

Location:
LMDZ6/trunk/libf/phylmd
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/calwake.f90

    r5499 r5763  
    5757  USE indice_sol_mod, ONLY: is_oce
    5858  USE print_control_mod, ONLY: lunout, prt_level
    59   USE lmdz_wake, ONLY : wake
     59  USE lmdz_wake, ONLY : wake, wake2
     60  USE alpale_mod, ONLY: iflag_wake
    6061  USE yomcst_mod_h
    6162IMPLICIT NONE
    … …  
    240241
    241242
    242   CALL wake(klon,klev,znatsurf, p, ph, pi, dtime, &
    243     te, qe, omgbe, &
    244     dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
    245     sigd0, Cin, &
    246     dtw, dqw, sigmaw, asigmaw, wdens, awdens, &                      ! state variables
    247     dth, hw, wape, fip, gfl, &
    248     dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, &
    249     dtke, dqke, omg, dp_deltomg, spread, cstar, &
    250     d_deltat_gw, &
    251     d_deltatw, d_deltaqw, d_sigmaw, d_asigmaw, d_wdens, d_awdens)      ! tendencies
     243  IF (iflag_wake/10 == 0) THEN
     244    CALL wake(klon,klev,znatsurf, p, ph, pi, dtime, &
     245      te, qe, omgbe, &
     246      dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
     247      sigd0, Cin, &
     248      dtw, dqw, sigmaw, asigmaw, wdens, awdens, &                      ! state variables
     249      dth, hw, wape, fip, gfl, &
     250      dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, &
     251      dtke, dqke, omg, dp_deltomg, spread, cstar, &
     252      d_deltat_gw, &
     253      d_deltatw, d_deltaqw, d_sigmaw, d_asigmaw, d_wdens, d_awdens)      ! tendencies
     254
     255  ELSE IF (iflag_wake/10 == 2) THEN
     256    CALL wake2(klon,klev,znatsurf, p, ph, pi, dtime, &
     257      te, qe, omgbe, &
     258      dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
     259      sigd0, Cin, &
     260      dtw, dqw, sigmaw, asigmaw, wdens, awdens, &                      ! state variables
     261      dth, hw, wape, fip, gfl, &
     262      dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, &
     263      dtke, dqke, omg, dp_deltomg, spread, cstar, &
     264      d_deltat_gw, &
     265      d_deltatw, d_deltaqw, d_sigmaw, d_asigmaw, d_wdens, d_awdens)      ! tendencies
     266
     267  END IF  !(iflag_wake/10 == 0)
     268
    252269
    253270!

  • LMDZ6/trunk/libf/phylmd/lmdz_wake.f90

    r5761 r5763  
    55  IMPLICIT NONE; PRIVATE
    66 
    7 !!!  LOGICAL, PARAMETER :: phys_sub=.false.
    8   LOGICAL, PARAMETER :: phys_sub=.true.
     7  LOGICAL, PARAMETER :: phys_sub=.false.
    98  LOGICAL            :: first_call=.true.
    109  !$OMP THREADPRIVATE(first_call)
    1110
    12   PUBLIC wake, wake_first
     11  PUBLIC wake, wake2, wake_first
    1312
    1413CONTAINS
    … …  
    4342                dth, hw, wape, fip, gfl, &
    4443                dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, &
    45 !!!                dtke, dqke, omg, dp_deltomg, wkspread, cstar, &   ! changes in notation
    46                 d_deltat_dcv, d_deltaq_dcv, omg, dp_deltomg, wkspread, cstar, &
     44                dtke, dqke, omg, dp_deltomg, wkspread, cstar, &
    4745                d_deltat_gw, &                                                      ! tendencies
    4846                d_deltatw2, d_deltaqw2, d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2)             ! tendencies
    … …  
    9896  ! wdens   : densite de poches
    9997  ! omgbdth: flux of Delta_Theta transported by LS omega
    100   ! d_deltat_dcv   : differential heating (wake - unpertubed)
    101   ! d_deltat_dcv   : differential moistening (wake - unpertubed)
     98  ! dtKE   : differential heating (wake - unpertubed)
     99  ! dqKE   : differential moistening (wake - unpertubed)
    102100  ! omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
    103101  ! dp_deltomg  : vertical gradient of omg (s-1)
    … …  
    170168  ! --------------------
    171169
    172   INTEGER,                          INTENT(IN)          :: klon,klev
     170  INTEGER, INTENT(IN) :: klon,klev
    173171  INTEGER, DIMENSION (klon),        INTENT(IN)          :: znatsurf
    174172  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: p, pi
    … …  
    199197  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: tx, qx
    200198  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls
    201 !!  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
    202   REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_dcv, d_deltaq_dcv
     199  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
    203200  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: wkspread    !  unused (jyg)
    204201  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: omgbdth, omg
    … …  
    270267  ! Sub-timestep tendencies and related variables
    271268  REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
    272   REAL, DIMENSION (klon, klev)                          :: d_deltat_dadv, d_deltaq_dadv
    273   REAL, DIMENSION (klon, klev)                          :: d_deltat_lsadv, d_deltaq_lsadv
    274   REAL, DIMENSION (klon, klev)                          :: d_deltat_entrp
    275   REAL, DIMENSION (klon, klev)                          :: d_deltat_spread, d_deltaq_spread
    276 
    277269  REAL, DIMENSION (klon, klev)                          :: d_tb, d_qb
    278   REAL, DIMENSION (klon, klev)                          :: d_tb_dadv, d_qb_dadv
    279   REAL, DIMENSION (klon, klev)                          :: d_tb_spread, d_qb_spread
    280 
    281270  REAL, DIMENSION (klon)                                :: d_wdens, d_awdens, d_sigmaw, d_asigmaw
    282271  REAL, DIMENSION (klon)                                :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
    … …  
    290279  LOGICAL, DIMENSION (klon)                             :: wk_adv, ok_qx_qw
    291280
    292   ! Tests
    293   REAL, DIMENSION (klon, klev)                          :: c5p
    294 
    295281  ! Autres variables internes
    296282  INTEGER                                               ::isubstep, k, i, igout
    … …  
    302288  REAL                                                  :: d_sigmaw_targ
    303289  REAL                                                  :: d_wdens_targ
    304 
    305   REAL, DIMENSION (klon)                                :: dsigspread  !rate of change of sigmaw due to spreading
    306290
    307291  REAL, DIMENSION (klon)                                :: sum_thx, sum_tx, sum_qx, sum_thvx
    … …  
    333317  REAL, DIMENSION (klon, klev)                          :: omgbdq
    334318
     319  REAL, DIMENSION (klon)                                :: ff, gg
    335320  REAL, DIMENSION (klon)                                :: wape2, cstar2, heff
    336321                                                       
    … …  
    342327  REAL, DIMENSION (klon)                                :: death_rate
    343328!!  REAL, DIMENSION (klon)                                :: nat_rate
    344   REAL, DIMENSION (klon, klev)                          :: entr   ! total entrainment into wakes (spread + birth)
    345   REAL, DIMENSION (klon, klev)                          :: entr_p ! entrainment into wakes (due to births)
    346   REAL, DIMENSION (klon, klev)                          :: detr   ! detrainment from wakes (due to deaths)
     329  REAL, DIMENSION (klon, klev)                          :: entr
     330  REAL, DIMENSION (klon, klev)                          :: detr
    347331
    348332  REAL, DIMENSION(klon)                                 :: sigmaw_in, asigmaw_in ! pour les prints
    … …  
    357341  REAL, DIMENSION (klon)                                :: omega
    358342  REAL, DIMENSION (klon)                                :: h_zzz
    359 
    360   !! Bidouilles
    361   REAL                                                  :: iwkadv
    362   REAL                                                  :: iokqxqw
    363343
    364344!print*,'WAKE LJYFz'
    … …  
    458438      dqls(:,:) = 0.
    459439      d_deltat_gw(:,:) = 0.
     440      d_tb(:,:) = 0.
     441      d_qb(:,:) = 0.
     442      d_deltatw(:,:) = 0.
     443      d_deltaqw(:,:) = 0.
     444      d_deltatw2(:,:) = 0.
     445      d_deltaqw2(:,:) = 0.
     446
     447      d_sig_gen2(:)   = 0.
     448      d_sig_death2(:) = 0.
     449      d_sig_col2(:)   = 0.
     450      d_sig_spread2(:)= 0.
     451      d_asig_death2(:) = 0.
     452      d_asig_iicol2(:) = 0.
     453      d_asig_aicol2(:) = 0.
     454      d_asig_spread2(:)= 0.
     455      d_asig_bnd2(:) = 0.
     456      d_asigmaw2(:) = 0.
     457!
     458      d_dens_gen2(:)   = 0.
     459      d_dens_death2(:) = 0.
     460      d_dens_col2(:)   = 0.
     461      d_dens_bnd2(:) = 0.
     462      d_wdens2(:) = 0.
     463      d_adens_bnd2(:) = 0.
     464      d_awdens2(:) = 0.
     465      d_adens_death2(:) = 0.
     466      d_adens_icol2(:) = 0.
     467      d_adens_acol2(:) = 0.
     468
     469      IF (iflag_wk_act == 0) THEN
     470        act(:) = 0.
     471      ELSEIF (iflag_wk_act == 1) THEN
     472        act(:) = 1.
     473      ENDIF
     474
     475!!  DO i = 1, klon
     476!!   sigmaw_in(i) = sigmaw(i)
     477!!  END DO
     478   sigmaw_in(:)  = sigmaw(:)
     479   asigmaw_in(:) = asigmaw(:)
     480!>jyg
     481!
     482  IF (iflag_wk_pop_dyn >= 1) THEN
     483    awdens_in(:) = awdens(:)
     484    wdens_in(:) = wdens(:)
     485!!    wdens(:) = wdens(:) + wgen(:)*dtime
     486!!    d_wdens2(:) = wgen(:)*dtime
     487!!  ELSE
     488  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     489
     490
     491  ! sigmaw1=sigmaw
     492  ! IF (sigd_con.GT.sigmaw1) THEN
     493  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
     494  ! ENDIF
     495  IF (iflag_wk_pop_dyn >= 1) THEN
     496    DO i = 1, klon
     497      d_dens_gen2(i)   = 0.
     498      d_dens_death2(i) = 0.
     499      d_dens_col2(i)   = 0.
     500      d_awdens2(i) = 0.
     501      IF (wdens(i) < wdensthreshold) THEN
     502  !!      wdens_targ = max(wdens(i),wdensmin)
     503        wdens_targ = max(wdens(i),wdensinit)
     504        d_dens_bnd2(i) = wdens_targ - wdens(i)
     505        d_wdens2(i) = wdens_targ - wdens(i)
     506        wdens(i) = wdens_targ
     507      ELSE
     508        d_dens_bnd2(i) = 0.
     509        d_wdens2(i) = 0.
     510      ENDIF  !! (wdens(i) < wdensthreshold)
     511    END DO
     512    IF (iflag_wk_pop_dyn >= 2) THEN
     513      DO i = 1, klon 
     514        IF (awdens(i) < wdensthreshold) THEN
     515!!          wdens_targ = min(max(awdens(i),wdensmin),wdens(i))
     516            wdens_targ = min(max(awdens(i),wdensinit),wdens(i))
     517            d_adens_bnd2(i) = wdens_targ - awdens(i)
     518            d_awdens2(i) = wdens_targ - awdens(i)
     519            awdens(i) = wdens_targ
     520        ELSE
     521            wdens_targ = min(awdens(i), wdens(i))
     522            d_adens_bnd2(i) = wdens_targ - awdens(i)
     523            d_awdens2(i) = wdens_targ - awdens(i)
     524            awdens(i) = wdens_targ
     525        ENDIF
     526      END DO
     527    ENDIF ! (iflag_wk_pop_dyn >= 2)
     528  ELSE 
     529    DO i = 1, klon
     530      d_awdens2(i) = 0.
     531      d_wdens2(i) = 0.
     532    END DO
     533  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     534!
     535  DO i = 1, klon
     536    sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
     537    d_sig_bnd2(i) = sigmaw_targ - sigmaw(i)
     538    d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
     539    sigmaw(i) = sigmaw_targ
     540  END DO
     541!
     542  IF (iflag_wk_pop_dyn == 3) THEN
     543     DO i = 1, klon 
     544        IF ((wdens(i)-awdens(i)) <= smallestreal) THEN
     545          sigmaw_targ = sigmaw(i)
     546        ELSE
     547          sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i))
     548        ENDIF
     549        d_asig_bnd2(i) = sigmaw_targ - asigmaw(i)
     550        d_asigmaw2(i) = sigmaw_targ - asigmaw(i)
     551        asigmaw(i) = sigmaw_targ
     552     END DO
     553  ENDIF ! (iflag_wk_pop_dyn == 3)
     554
     555  wape(:) = 0.
     556  wape2(:) = 0.
     557  d_sigmaw(:) = 0.
     558  d_asigmaw(:) = 0.
     559  ktopw(:) = 0
     560!
     561!<jyg
     562dth(:,:) = 0.
     563tx(:,:) = 0.
     564qx(:,:) = 0.
     565dtke(:,:) = 0.
     566dqke(:,:) = 0.
     567wkspread(:,:) = 0.
     568omgbdth(:,:) = 0.
     569omg(:,:) = 0.
     570dp_omgb(:,:) = 0.
     571dp_deltomg(:,:) = 0.
     572hw(:) = 0.
     573wape(:) = 0.
     574fip(:) = 0.
     575gfl(:) = 0.
     576cstar(:) = 0.
     577ktopw(:) = 0
     578!
     579!  Vertical advection local variables
     580omgbw(:,:) = 0.
     581omgtop(:) = 0
     582dp_omgbw(:,:) = 0.
     583omgbdq(:,:) = 0.
     584
     585!>jyg
     586!
     587  IF (prt_level>=10) THEN
     588    PRINT *, 'wake-1, sigmaw(igout) ', sigmaw(igout)
     589    PRINT *, 'wake-1, deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev)
     590    PRINT *, 'wake-1, deltaqw(igout,k) ', (k,deltaqw(igout,k), k=1,klev)
     591    PRINT *, 'wake-1, dowwdraughts, amdwn(igout,k) ', (k,amdwn(igout,k), k=1,klev)
     592    PRINT *, 'wake-1, dowwdraughts, dtdwn(igout,k) ', (k,dtdwn(igout,k), k=1,klev)
     593    PRINT *, 'wake-1, dowwdraughts, dqdwn(igout,k) ', (k,dqdwn(igout,k), k=1,klev)
     594    PRINT *, 'wake-1, updraughts, amup(igout,k) ', (k,amup(igout,k), k=1,klev)
     595    PRINT *, 'wake-1, updraughts, dta(igout,k) ', (k,dta(igout,k), k=1,klev)
     596    PRINT *, 'wake-1, updraughts, dqa(igout,k) ', (k,dqa(igout,k), k=1,klev)
     597  ENDIF
     598
     599  ! 2. - Prognostic part
     600  ! --------------------
     601
     602
     603  ! 2.1 - Undisturbed area and Wake integrals
     604  ! ---------------------------------------------------------
     605
     606  DO i = 1, klon
     607    z(i) = 0.
     608    ktop(i) = 0
     609    kupper(i) = 0
     610    sum_thx(i) = 0.
     611    sum_tx(i) = 0.
     612    sum_qx(i) = 0.
     613    sum_thvx(i) = 0.
     614    sum_dth(i) = 0.
     615    sum_dq(i) = 0.
     616    sum_dtdwn(i) = 0.
     617    sum_dqdwn(i) = 0.
     618
     619    av_thx(i) = 0.
     620    av_tx(i) = 0.
     621    av_qx(i) = 0.
     622    av_thvx(i) = 0.
     623    av_dth(i) = 0.
     624    av_dq(i) = 0.
     625    av_dtdwn(i) = 0.
     626    av_dqdwn(i) = 0.
     627  END DO
     628
     629  ! Distance between wakes
     630  DO i = 1, klon
     631    ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
     632  END DO
     633  ! Potential temperatures and humidity
     634  ! ----------------------------------------------------------
     635  DO k = 1, klev
     636    DO i = 1, klon
     637      ! write(*,*)'wake 1',i,k,RD,tb(i,k)
     638      rho(i, k) = p(i, k)/(RD*tb(i,k))
     639      ! write(*,*)'wake 2',rho(i,k)
     640      IF (k==1) THEN
     641        ! write(*,*)'wake 3',i,k,rd,tb(i,k)
     642        rhoh(i, k) = ph(i, k)/(RD*tb(i,k))
     643        ! write(*,*)'wake 4',i,k,rd,tb(i,k)
     644        zhh(i, k) = 0
     645      ELSE
     646        ! write(*,*)'wake 5',rd,(tb(i,k)+tb(i,k-1))
     647        rhoh(i, k) = ph(i, k)*2./(RD*(tb(i,k)+tb(i,k-1)))
     648        ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
     649        zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
     650      END IF
     651      ! write(*,*)'wake 7',ppi(i,k)
     652      thb(i, k) = tb(i, k)/ppi(i, k)
     653      thx(i, k) = (tb(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
     654      tx(i, k) = tb(i, k) - deltatw(i, k)*sigmaw(i)
     655      qx(i, k) = qb(i, k) - deltaqw(i, k)*sigmaw(i)
     656      ! write(*,*)'wake 8',(RD*(tb(i,k)+deltatw(i,k)))
     657      dth(i, k) = deltatw(i, k)/ppi(i, k)
     658    END DO
     659  END DO
     660
     661  DO k = 1, klev - 1
     662    DO i = 1, klon
     663      IF (k==1) THEN
     664        n2(i, k) = 0
     665      ELSE
     666        n2(i, k) = amax1(0., -RG**2/thb(i,k)*rho(i,k)*(thb(i,k+1)-thb(i,k-1))/ &
     667                             (p(i,k+1)-p(i,k-1)))
     668      END IF
     669      zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
     670
     671      cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
     672      tgw(i, k) = coefgw*cgw(i, k)/ll(i)
     673    END DO
     674  END DO
     675
     676  DO i = 1, klon
     677    n2(i, klev) = 0
     678    zh(i, klev) = 0
     679    cgw(i, klev) = 0
     680    tgw(i, klev) = 0
     681  END DO
     682
     683 
     684  ! Choose an integration bound well above wake top
     685  ! -----------------------------------------------------------------
     686
     687  ! Determine Wake top pressure (Ptop) from buoyancy integral
     688  ! --------------------------------------------------------
     689
     690   Do i=1, klon
     691       wk_adv(i) = .True.
     692   Enddo
     693   Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
     694                    dth, hw0, rho, delta_t_min, &
     695                    ktop, wk_adv, h_zzz, ptop1, ktop1)
     696 
     697   !!print'("pkupper APPEL ",7i6)',0,int(ptop/100.),int(ptop1/100.),int(pupper/100.),ktop,ktop1,kupper
     698   
     699   IF (prt_level>=10) THEN
     700     PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
     701  ENDIF
     702
     703  ! -5/ Set deltatw & deltaqw to 0 above kupper
     704
     705  DO k = 1, klev
     706    DO i = 1, klon
     707      IF (k>=kupper(i)) THEN
     708        deltatw(i, k) = 0.
     709        deltaqw(i, k) = 0.
     710        d_deltatw2(i,k) = -deltatw0(i,k)
     711        d_deltaqw2(i,k) = -deltaqw0(i,k)
     712      END IF
     713    END DO
     714  END DO
     715
     716
     717  ! Vertical gradient of LS omega
     718
     719  DO k = 1, klev
     720    DO i = 1, klon
     721      IF (k<=kupper(i)) THEN
     722        dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
     723      END IF
     724    END DO
     725  END DO
     726
     727  ! Integrals (and wake top level number)
     728  ! --------------------------------------
     729
     730  ! Initialize sum_thvx to 1st level virt. pot. temp.
     731
     732  DO i = 1, klon
     733    z(i) = 1.
     734    dz(i) = 1.
     735    sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i)
     736    sum_dth(i) = 0.
     737  END DO
     738
     739  DO k = 1, klev
     740    DO i = 1, klon
     741      dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
     742      IF (dz(i)>0) THEN
     743              ! LJYF : ecriture pas sympa avec un tableau z(i) qui n'est pas utilise come tableau
     744        z(i) = z(i) + dz(i)
     745        sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i)
     746        sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i)
     747        sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i)
     748        sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i)
     749        sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
     750        sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
     751        sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
     752        sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
     753      END IF
     754    END DO
     755  END DO
     756
     757  DO i = 1, klon
     758    hw0(i) = z(i)
     759  END DO
     760
     761
     762  ! 2.1 - WAPE and mean forcing computation
     763  ! ---------------------------------------
     764
     765  ! ---------------------------------------
     766
     767  ! Means
     768
     769  DO i = 1, klon
     770    av_thx(i) = sum_thx(i)/hw0(i)
     771    av_tx(i) = sum_tx(i)/hw0(i)
     772    av_qx(i) = sum_qx(i)/hw0(i)
     773    av_thvx(i) = sum_thvx(i)/hw0(i)
     774    ! av_thve = sum_thve/hw0
     775    av_dth(i) = sum_dth(i)/hw0(i)
     776    av_dq(i) = sum_dq(i)/hw0(i)
     777    av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     778    av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     779
     780    wape(i) = -RG*hw0(i)*(av_dth(i)+ &
     781        epsim1*(av_thx(i)*av_dq(i)+av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i)
     782
     783  END DO
     784IF (CPPKEY_IOPHYS_WK) THEN
     785  IF (.not.phys_sub) CALL iophys_ecrit('wape_a',1,'wape_a','J/kg',wape)
     786END IF
     787
     788  ! 2.2 Prognostic variable update
     789  ! ------------------------------
     790
     791  ! Filter out bad wakes
     792
     793  DO k = 1, klev
     794    DO i = 1, klon
     795      IF (wape(i)<0.) THEN
     796        deltatw(i, k) = 0.
     797        deltaqw(i, k) = 0.
     798        dth(i, k) = 0.
     799        d_deltatw2(i,k) = -deltatw0(i,k)
     800        d_deltaqw2(i,k) = -deltaqw0(i,k)
     801      END IF
     802    END DO
     803  END DO
     804
     805  DO i = 1, klon
     806    IF (wape(i)<0.) THEN
     807!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
     808      sigmaw_targ = max(sigmad, sigd_con(i))
     809      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     810      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     811      sigmaw(i) = sigmaw_targ
     812    ENDIF  !!  (wape(i)<0.)
     813  ENDDO
     814  !
     815  IF (iflag_wk_pop_dyn == 3) THEN
     816    DO i = 1, klon
     817      IF (wape(i)<0.) THEN
     818        sigmaw_targ = max(sigmad, sigd_con(i))
     819        d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     820        d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     821        asigmaw(i) = sigmaw_targ
     822      ENDIF  !!  (wape(i)<0.)
     823    ENDDO
     824  ENDIF  !!  (iflag_wk_pop_dyn == 3)
     825
     826  DO i = 1, klon
     827    IF (wape(i)<0.) THEN
     828      wape(i) = 0.
     829      cstar(i) = 0.
     830      hw(i) = hwmin
     831      fip(i) = 0.
     832      gwake(i) = .FALSE.
     833    ELSE
     834      hw(i) = hw0(i)
     835      cstar(i) = stark*sqrt(2.*wape(i))
     836      gwake(i) = .TRUE.
     837    END IF
     838  END DO
     839  !
     840
     841  ! Check qx and qw positivity
     842  ! --------------------------
     843  DO i = 1, klon
     844    q0_min(i) = min((qb(i,1)-sigmaw(i)*deltaqw(i,1)),  &
     845                    (qb(i,1)+(1.-sigmaw(i))*deltaqw(i,1)))
     846  END DO
     847  DO k = 2, klev
     848    DO i = 1, klon
     849      q1_min(i) = min((qb(i,k)-sigmaw(i)*deltaqw(i,k)), &
     850                      (qb(i,k)+(1.-sigmaw(i))*deltaqw(i,k)))
     851      IF (q1_min(i)<=q0_min(i)) THEN
     852        q0_min(i) = q1_min(i)
     853      END IF
     854    END DO
     855  END DO
     856
     857  DO i = 1, klon
     858    ok_qx_qw(i) = q0_min(i) >= 0.
     859    alpha(i) = 1.
     860    alpha_tot(i) = 1.
     861  END DO
     862
     863  IF (prt_level>=10) THEN
     864    PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ', &
     865                      sigmaw(igout), cstar(igout), wape(igout), ktop(igout)
     866  ENDIF
     867
     868
     869  ! C -----------------------------------------------------------------
     870  ! Sub-time-stepping
     871  ! -----------------
     872
     873!    wk_nsub and dtimesub definitions moved to begining of routine.
     874!!  wk_nsub = 10
     875!!  dtimesub = dtime/wk_nsub
     876
     877 
     878  ! ------------------------------------------------------------------------
     879  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     880  ! ------------------------------------------------------------------------
     881  !
     882  DO isubstep = 1, wk_nsub
     883  !
     884  ! ------------------------------------------------------------------------
     885    ! wk_adv is the logical flag enabling wake evolution in the time advance
     886    ! loop
     887    DO i = 1, klon
     888      wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
     889    END DO
     890    IF (prt_level>=10) THEN
     891      PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', &
     892                          isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout)
     893     
     894    ENDIF
     895
     896    ! cc nrlmd   Ajout d'un recalcul de wdens dans le cas d'un entrainement
     897    ! negatif de ktop a kupper --------
     898    ! cc           On calcule pour cela une densite wdens0 pour laquelle on
     899    ! aurait un entrainement nul ---
     900    !jyg<
     901    ! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou
     902    ! les poches sont insuffisantes pour accueillir tout le flux de masse
     903    ! des descentes unsaturees. Nous faisons alors l'hypothese que la
     904    ! convection profonde cree directement de nouvelles poches, sans passer
     905    ! par les thermiques. La nouvelle valeur de wdens est alors imposee.
     906
     907    DO i = 1, klon
     908      ! c       print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
     909      ! c     $           isubstep,wk_adv(i),cstar(i),wape(i)
     910      IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
     911        IF ( iflag_wk_profile == 0 ) THEN
     912           omg(i, kupper(i)+1)=-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
     913                               RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
     914        ELSE
     915           omg(i, kupper(i)+1)=0.
     916        ENDIF
     917        wdens0 = (sigmaw(i)/(4.*3.14))* &
     918          ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2)
     919        IF (prt_level >= 10) THEN
     920             print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', &
     921                     omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)
     922        ENDIF
     923        IF (wdens(i)<=wdens0*1.1) THEN
     924          IF (iflag_wk_pop_dyn >= 1) THEN
     925             d_dens_bnd2(i) = d_dens_bnd2(i) + wdens0 - wdens(i)
     926             d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
     927          ENDIF
     928          wdens(i) = wdens0
     929        END IF
     930      END IF
     931    END DO
     932
     933    IF (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl) THEN
     934!!--------------------------------------------------------
     935!!Bug : computing gfl and rad_wk before changing sigmaw
     936!!      This bug exists only for iflag_wk_pop_dyn=0. Otherwise, gfl and rad_wk
     937!!      are computed within  wake_popdyn
     938!!--------------------------------------------------------
     939      DO i = 1, klon
     940        IF (wk_adv(i)) THEN
     941          gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     942          rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
     943        END IF
     944      END DO
     945    ENDIF   ! (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl)
     946!!--------------------------------------------------------
     947
     948    DO i = 1, klon
     949      IF (wk_adv(i)) THEN
     950        sigmaw_targ = min(sigmaw(i), sigmaw_max)
     951        d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     952        d_sigmaw2(i)  = d_sigmaw2(i)  + sigmaw_targ - sigmaw(i)
     953        sigmaw(i) = sigmaw_targ
     954      END IF
     955    END DO
     956
     957    IF (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl) THEN
     958!!--------------------------------------------------------
     959!!Fix : computing gfl and rad_wk after changing sigmaw
     960!!--------------------------------------------------------
     961      DO i = 1, klon
     962        IF (wk_adv(i)) THEN
     963          gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     964          rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
     965        END IF
     966      END DO
     967    ENDIF   ! (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl)
     968!!--------------------------------------------------------
     969
     970    IF (iflag_wk_pop_dyn >= 1) THEN
     971  !  The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
     972  !  Here, it has to be set to zero.
     973      death_rate(:) = 0.
     974    ENDIF
     975 
     976    IF (iflag_wk_pop_dyn >= 3) THEN
     977      DO i = 1, klon
     978        IF (wk_adv(i)) THEN
     979         sigmaw_targ = min(asigmaw(i), sigmaw_max)
     980         d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     981         d_asigmaw2(i)  = d_asigmaw2(i)  + sigmaw_targ - asigmaw(i)
     982         asigmaw(i) = sigmaw_targ
     983        ENDIF
     984      ENDDO
     985    ENDIF
     986 
     987!!--------------------------------------------------------
     988!!--------------------------------------------------------
     989    IF (iflag_wk_pop_dyn == 1) THEN
     990  !
     991     CALL wake_popdyn_1 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
     992                  wdensmin, &
     993                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
     994                  d_awdens, d_wdens, d_sigmaw, &
     995                  iflag_wk_act, wk_adv, cin, wape, &
     996                  drdt, &
     997                  d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     998                  d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     999                  d_wdens_targ, d_sigmaw_targ)
     1000                     
     1001   
     1002!!--------------------------------------------------------
     1003    ELSEIF (iflag_wk_pop_dyn == 2) THEN
     1004  !
     1005     CALL wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, &
     1006                             wdensmin, &
     1007                             sigmaw, wdens, awdens, &   !! state variables
     1008                             gfl, cstar, cin, wape, rad_wk, &
     1009                             d_sigmaw, d_wdens, d_awdens, &  !! tendencies                             
     1010                             cont_fact, &
     1011                             d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     1012                             d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     1013                             d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
     1014sigmaw=sigmaw-d_sigmaw
     1015wdens=wdens-d_wdens
     1016awdens=awdens-d_awdens
     1017
     1018!!--------------------------------------------------------
     1019    ELSEIF (iflag_wk_pop_dyn == 3) THEN
     1020IF (CPPKEY_IOPHYS_WK) THEN
     1021    IF (phys_sub) THEN
     1022     CALL iophys_ecrit('ptop',1,'ptop','Pa',ptop)
     1023     CALL iophys_ecrit('sigmaw',1,'sigmaw','',sigmaw)
     1024     CALL iophys_ecrit('asigmaw',1,'asigmaw','',asigmaw)
     1025     CALL iophys_ecrit('wdens',1,'wdens','1/m2',wdens)
     1026     CALL iophys_ecrit('awdens',1,'awdens','1/m2',awdens)
     1027     CALL iophys_ecrit('rad_wk',1,'rad_wk','m',rad_wk)
     1028     CALL iophys_ecrit('arad_wk',1,'arad_wk','m',arad_wk)
     1029     CALL iophys_ecrit('irad_wk',1,'irad_wk','m',irad_wk)
     1030    ENDIF
     1031END IF
     1032  !
     1033     CALL wake_popdyn_3 ( klon, klev, phys_sub, wk_adv, dtimesub, wgen, &
     1034                             wdensmin, &
     1035                             sigmaw, asigmaw, wdens, awdens, &   !! state variables
     1036                             gfl, agfl, cstar, cin, wape, &
     1037                             rad_wk, arad_wk, irad_wk, &
     1038                             d_sigmaw, d_asigmaw, d_wdens, d_awdens, &  !! tendencies                             
     1039                             d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     1040                             d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd, &
     1041                             d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     1042                             d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
     1043sigmaw=sigmaw-d_sigmaw
     1044asigmaw=asigmaw-d_asigmaw
     1045wdens=wdens-d_wdens
     1046awdens=awdens-d_awdens
     1047   
     1048!!--------------------------------------------------------
     1049    ELSEIF (iflag_wk_pop_dyn == 0) THEN
     1050   
     1051    ! cc nrlmd
     1052
     1053      DO i = 1, klon
     1054        IF (wk_adv(i)) THEN
     1055
     1056          ! cc nrlmd          Introduction du taux de mortalite des poches et
     1057          ! test sur sigmaw_max=0.4
     1058          ! cc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
     1059          IF (sigmaw(i)>=sigmaw_max) THEN
     1060            death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
     1061          ELSE
     1062            death_rate(i) = 0.
     1063          END IF
     1064   
     1065          d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* &
     1066            dtimesub
     1067          ! $              - nat_rate(i)*sigmaw(i)*dtimesub
     1068          ! c        print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
     1069          ! c     $  death_rate(i),ktop(i),kupper(i)',
     1070          ! c     $              d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
     1071          ! c     $  death_rate(i),ktop(i),kupper(i)
     1072   
     1073          ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
     1074          ! sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
     1075          ! wdens = wdens0/(10.*sigmaw)
     1076          ! sigmaw =max(sigmaw,sigd_con)
     1077          ! sigmaw =max(sigmaw,sigmad)
     1078        END IF
     1079      END DO
     1080
     1081    ENDIF   !  (iflag_wk_pop_dyn == 1) ... ELSEIF (iflag_wk_pop_dyn == 0)
     1082!!--------------------------------------------------------
     1083!!--------------------------------------------------------
     1084
     1085IF (CPPKEY_IOPHYS_WK) THEN
     1086    IF (phys_sub) THEN
     1087     CALL iophys_ecrit('wdensa',1,'wdensa','m',wdens)
     1088     CALL iophys_ecrit('awdensa',1,'awdensa','m',awdens)
     1089     CALL iophys_ecrit('sigmawa',1,'sigmawa','m',sigmaw)
     1090     CALL iophys_ecrit('asigmawa',1,'asigmawa','m',asigmaw)
     1091    ENDIF
     1092END IF
     1093    ! calcul de la difference de vitesse verticale poche - zone non perturbee
     1094    ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
     1095    ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
     1096    ! IM 060208 au niveau k=1...
     1097    !JYG 161013 Correction : maintenant omg est dimensionne a klev.
     1098    DO k = 1, klev
     1099      DO i = 1, klon
     1100        IF (wk_adv(i)) THEN !!! nrlmd
     1101          dp_deltomg(i, k) = 0.
     1102        END IF
     1103      END DO
     1104    END DO
     1105    DO k = 1, klev
     1106      DO i = 1, klon
     1107        IF (wk_adv(i)) THEN !!! nrlmd
     1108          omg(i, k) = 0.
     1109        END IF
     1110      END DO
     1111    END DO
     1112
     1113    DO i = 1, klon
     1114      IF (wk_adv(i)) THEN
     1115        z(i) = 0.
     1116        omg(i, 1) = 0.
     1117        dp_deltomg(i, 1) = -(gfl(i)*cstar(i))/(sigmaw(i)*(1-sigmaw(i)))
     1118      END IF
     1119    END DO
     1120
     1121    DO k = 2, klev
     1122      DO i = 1, klon
     1123        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
     1124          dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*RG)
     1125          z(i) = z(i) + dz(i)
     1126          dp_deltomg(i, k) = dp_deltomg(i, 1)
     1127          omg(i, k) = dp_deltomg(i, 1)*z(i)
     1128        END IF
     1129      END DO
     1130    END DO
     1131
     1132    DO i = 1, klon
     1133      IF (wk_adv(i)) THEN
     1134        dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*RG)
     1135        ztop(i) = z(i) + dztop(i)
     1136        omgtop(i) = dp_deltomg(i, 1)*ztop(i)
     1137      END IF
     1138    END DO
     1139
     1140    IF (prt_level>=10) THEN
     1141      PRINT *, 'wake-4.2, omg(igout,k) ', (k,omg(igout,k), k=1,klev)
     1142      PRINT *, 'wake-4.2, omgtop(igout), ptop(igout), ktop(igout) ', &
     1143                          omgtop(igout), ptop(igout), ktop(igout)
     1144    ENDIF
     1145
     1146    ! -----------------
     1147    ! From m/s to Pa/s
     1148    ! -----------------
     1149
     1150    DO i = 1, klon
     1151      IF (wk_adv(i)) THEN
     1152        omgtop(i) = -rho(i, ktop(i))*RG*omgtop(i)
     1153!! LJYF        dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
     1154        dp_deltomg(i, 1) = omgtop(i)/min(ptop(i)-ph(i,1),-smallestreal)
     1155      END IF
     1156    END DO
     1157
     1158    DO k = 1, klev
     1159      DO i = 1, klon
     1160        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
     1161          omg(i, k) = -rho(i, k)*RG*omg(i, k)
     1162          dp_deltomg(i, k) = dp_deltomg(i, 1)
     1163        END IF
     1164      END DO
     1165    END DO
     1166
     1167    ! raccordement lineaire de omg de ptop a pupper
     1168
     1169    DO i = 1, klon
     1170      IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
     1171        IF ( iflag_wk_profile == 0 ) THEN
     1172           omg(i, kupper(i)+1) =-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
     1173          RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
     1174        ELSE
     1175           omg(i, kupper(i)+1) = 0.
     1176        ENDIF
     1177        dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
     1178          (ptop(i)-pupper(i))
     1179      END IF
     1180    END DO
     1181
     1182    ! c      DO i=1,klon
     1183    ! c        print*,'Pente entre 0 et kupper (reference)'
     1184    ! c     $           ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
     1185    ! c        print*,'Pente entre ktop et kupper'
     1186    ! c     $   ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
     1187    ! c      ENDDO
     1188    ! c
     1189    DO k = 1, klev
     1190      DO i = 1, klon
     1191        IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
     1192          dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
     1193          omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
     1194        END IF
     1195      END DO
     1196    END DO
     1197!!    print *,'omg(igout,k) ', (k,omg(igout,k),k=1,klev)
     1198    ! cc nrlmd
     1199    ! c      DO i=1,klon
     1200    ! c      print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
     1201    ! c      END DO
     1202    ! cc
     1203
     1204
     1205    ! --    Compute wake average vertical velocity omgbw
     1206
     1207
     1208    DO k = 1, klev
     1209      DO i = 1, klon
     1210        IF (wk_adv(i)) THEN
     1211          omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
     1212        END IF
     1213      END DO
     1214    END DO
     1215    ! --    and its vertical gradient dp_omgbw
     1216
     1217    DO k = 1, klev-1
     1218      DO i = 1, klon
     1219        IF (wk_adv(i)) THEN
     1220          dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
     1221        END IF
     1222      END DO
     1223    END DO
     1224    DO i = 1, klon
     1225      IF (wk_adv(i)) THEN
     1226          dp_omgbw(i, klev) = 0.
     1227      END IF
     1228    END DO
     1229
     1230    ! --    Upstream coefficients for omgb velocity
     1231    ! --    (alpha_up(k) is the coefficient of the value at level k)
     1232    ! --    (1-alpha_up(k) is the coefficient of the value at level k-1)
     1233    DO k = 1, klev
     1234      DO i = 1, klon
     1235        IF (wk_adv(i)) THEN
     1236          alpha_up(i, k) = 0.
     1237          IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
     1238        END IF
     1239      END DO
     1240    END DO
     1241
     1242    ! Matrix expressing [The,deltatw] from  [Th1,Th2]
     1243
     1244    DO i = 1, klon
     1245      IF (wk_adv(i)) THEN
     1246        rre1(i) = 1. - sigmaw(i)
     1247        rre2(i) = sigmaw(i)
     1248      END IF
     1249    END DO
     1250    rrd1 = -1.
     1251    rrd2 = 1.
     1252
     1253    ! --    Get [Th1,Th2], dth and [q1,q2]
     1254
     1255    DO k = 1, klev
     1256      DO i = 1, klon
     1257        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
     1258          dth(i, k) = deltatw(i, k)/ppi(i, k)
     1259          th1(i, k) = thb(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
     1260          th2(i, k) = thb(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
     1261          q1(i, k) = qb(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
     1262          q2(i, k) = qb(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
     1263        END IF
     1264      END DO
     1265    END DO
     1266
     1267    DO i = 1, klon
     1268      IF (wk_adv(i)) THEN !!! nrlmd
     1269        d_th1(i, 1) = 0.
     1270        d_th2(i, 1) = 0.
     1271        d_dth(i, 1) = 0.
     1272        d_q1(i, 1) = 0.
     1273        d_q2(i, 1) = 0.
     1274        d_dq(i, 1) = 0.
     1275      END IF
     1276    END DO
     1277
     1278    DO k = 2, klev
     1279      DO i = 1, klon
     1280        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
     1281          d_th1(i, k) = th1(i, k-1) - th1(i, k)
     1282          d_th2(i, k) = th2(i, k-1) - th2(i, k)
     1283          d_dth(i, k) = dth(i, k-1) - dth(i, k)
     1284          d_q1(i, k) = q1(i, k-1) - q1(i, k)
     1285          d_q2(i, k) = q2(i, k-1) - q2(i, k)
     1286          d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
     1287        END IF
     1288      END DO
     1289    END DO
     1290
     1291    DO i = 1, klon
     1292      IF (wk_adv(i)) THEN
     1293        omgbdth(i, 1) = 0.
     1294        omgbdq(i, 1) = 0.
     1295      END IF
     1296    END DO
     1297
     1298    DO k = 2, klev
     1299      DO i = 1, klon
     1300        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN !   loop on interfaces
     1301          omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
     1302          omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
     1303        END IF
     1304      END DO
     1305    END DO
     1306
     1307!!    IF (prt_level>=10) THEN
     1308    IF (prt_level>=10 .and. wk_adv(igout)) THEN
     1309      PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,kupper(igout))
     1310      PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,kupper(igout))
     1311      PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,kupper(igout))
     1312      PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,kupper(igout))
     1313    ENDIF
     1314
     1315    ! -----------------------------------------------------------------
     1316    DO k = 1, klev-1
     1317      DO i = 1, klon
     1318        IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
     1319          ! -----------------------------------------------------------------
     1320
     1321          ! Compute redistribution (advective) term
     1322
     1323          d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
     1324            (rrd1*omg(i,k)*sigmaw(i)*d_th1(i,k) - &
     1325             rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1)- &
     1326             (1.-alpha_up(i,k))*omgbdth(i,k)- &
     1327             alpha_up(i,k+1)*omgbdth(i,k+1))*ppi(i, k)
     1328!           print*,'d_d,k_ptop_provis(i)eltatw=', k, d_deltatw(i,k)
     1329
     1330          d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
     1331            (rrd1*omg(i,k)*sigmaw(i)*d_q1(i,k)- &
     1332             rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1)- &
     1333             (1.-alpha_up(i,k))*omgbdq(i,k)- &
     1334             alpha_up(i,k+1)*omgbdq(i,k+1))
     1335!           print*,'d_deltaqw=', k, d_deltaqw(i,k)
     1336
     1337          ! and increment large scale tendencies
     1338
     1339
     1340
     1341
     1342          ! C
     1343          ! -----------------------------------------------------------------
     1344          d_tb(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)- &
     1345                                  rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/ &
     1346                                 (ph(i,k)-ph(i,k+1)) &
     1347                                 -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/ &
     1348                                 (ph(i,k)-ph(i,k+1)) )*ppi(i, k)
     1349
     1350          d_qb(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)- &
     1351                                  rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/ &
     1352                                 (ph(i,k)-ph(i,k+1)) &
     1353                                 -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i,k+1))/ &
     1354                                 (ph(i,k)-ph(i,k+1)) )
     1355        ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
     1356          d_tb(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k)
     1357
     1358          d_qb(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)/(ph(i,k)-ph(i,k+1)))
     1359
     1360        END IF
     1361        ! cc
     1362      END DO
     1363    END DO
     1364    ! ------------------------------------------------------------------
     1365
     1366    IF (prt_level>=10) THEN
     1367      PRINT *, 'wake-4.3, d_deltatw(igout,k) ', (k,d_deltatw(igout,k), k=1,klev)
     1368      PRINT *, 'wake-4.3, d_deltaqw(igout,k) ', (k,d_deltaqw(igout,k), k=1,klev)
     1369    ENDIF
     1370
     1371    ! Increment state variables
     1372!jyg<
     1373    IF (iflag_wk_pop_dyn >= 1) THEN
     1374      DO k = 1, klev
     1375        DO i = 1, klon
     1376          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     1377            detr(i,k) = - d_sig_death(i) - d_sig_col(i)     
     1378            entr(i,k) = d_sig_gen(i)
     1379          ENDIF
     1380        ENDDO
     1381      ENDDO
     1382      ELSE  ! (iflag_wk_pop_dyn >= 1)
     1383      DO k = 1, klev
     1384        DO i = 1, klon
     1385          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     1386            detr(i, k) = 0.
     1387   
     1388            entr(i, k) = 0.
     1389          ENDIF
     1390        ENDDO
     1391      ENDDO
     1392    ENDIF  ! (iflag_wk_pop_dyn >= 1)
     1393
     1394   
     1395
     1396    DO k = 1, klev
     1397      DO i = 1, klon
     1398        ! cc nrlmd       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
     1399        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     1400          ! cc
     1401
     1402
     1403
     1404          ! Coefficient de repartition
     1405
     1406          crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
     1407            (ph(i,kupper(i))-ph(i,1))
     1408          crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ &
     1409            (ph(i,1)-ph(i,kupper(i)))
     1410
     1411
     1412          ! Reintroduce compensating subsidence term.
     1413
     1414          ! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
     1415          ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
     1416          ! .                   /(1-sigmaw)
     1417          ! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
     1418          ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
     1419          ! .                   /(1-sigmaw)
     1420
     1421          ! dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
     1422          ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
     1423          ! .                   /(1-sigmaw)
     1424          ! dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
     1425          ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
     1426          ! .                   /(1-sigmaw)
     1427
     1428          dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
     1429          dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
     1430          ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
     1431
     1432!
     1433
     1434          ! cc nrlmd          Prise en compte du taux de mortalite
     1435          ! cc               Definitions de entr, detr
     1436!jyg<
     1437!!            detr(i, k) = 0.
     1438!!   
     1439!!            entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
     1440!!              sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)
     1441!!
     1442            entr(i, k) = entr(i,k) + gfl(i)*cstar(i) + &
     1443                         sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)   
     1444!>jyg
     1445            wkspread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
     1446
     1447          ! cc        wkspread(i,k) =
     1448          ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
     1449          ! cc     $  sigmaw(i)
     1450
     1451
     1452          ! ajout d'un effet onde de gravite -Tgw(k)*deltatw(k) 03/02/06 YU
     1453          ! Jingmei
     1454
     1455          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
     1456          ! &  Tgw(i,k),deltatw(i,k)
     1457          d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)*deltatw(i, k)* &
     1458            dtimesub
     1459          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
     1460          ff(i) = d_deltatw(i, k)/dtimesub
     1461
     1462          ! Sans GW
     1463
     1464          ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-wkspread(k)*deltatw(k))
     1465
     1466          ! GW formule 1
     1467
     1468          ! deltatw(k) = deltatw(k)+dtimesub*
     1469          ! $         (ff+dtKE(k) - wkspread(k)*deltatw(k)-Tgw(k)*deltatw(k))
     1470
     1471          ! GW formule 2
     1472
     1473          IF (dtimesub*tgw(i,k)<1.E-10) THEN
     1474            d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k) - &
     1475               entr(i,k)*deltatw(i,k)/sigmaw(i) - &
     1476               (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) - & ! cc
     1477               tgw(i,k)*deltatw(i,k) )
     1478          ELSE
     1479            d_deltatw(i, k) = 1/tgw(i, k)*(1-exp(-dtimesub*tgw(i,k)))* &
     1480               (ff(i)+dtke(i,k) - &
     1481                entr(i,k)*deltatw(i,k)/sigmaw(i) - &
     1482                (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) - &
     1483                tgw(i,k)*deltatw(i,k) )
     1484          END IF
     1485
     1486          dth(i, k) = deltatw(i, k)/ppi(i, k)
     1487
     1488          gg(i) = d_deltaqw(i, k)/dtimesub
     1489
     1490          d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k) - &
     1491            entr(i,k)*deltaqw(i,k)/sigmaw(i) - &
     1492            (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k)/(1.-sigmaw(i)))
     1493          ! cc
     1494
     1495          ! cc nrlmd
     1496          ! cc       d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
     1497          ! cc       d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
     1498          ! cc
     1499        END IF
     1500      END DO
     1501    END DO
     1502
     1503
     1504    ! Scale tendencies so that water vapour remains positive in w and x.
     1505
     1506    CALL wake_vec_modulation(klon, klev, wk_adv, epsilon_loc, qb, d_qb, deltaqw, &
     1507      d_deltaqw, sigmaw, d_sigmaw, alpha)
     1508    !
     1509    ! Alpha_tot = Product of all the alpha's
     1510    DO i = 1, klon
     1511      IF (wk_adv(i)) THEN
     1512        alpha_tot(i) = alpha_tot(i)*alpha(i)   
     1513      END IF
     1514    END DO
     1515
     1516    ! cc nrlmd
     1517    ! c      print*,'alpha'
     1518    ! c      do i=1,klon
     1519    ! c         print*,alpha(i)
     1520    ! c      end do
     1521    ! cc
     1522    DO k = 1, klev
     1523      DO i = 1, klon
     1524        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     1525          d_tb(i, k) = alpha(i)*d_tb(i, k)
     1526          d_qb(i, k) = alpha(i)*d_qb(i, k)
     1527          d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
     1528          d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
     1529          d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
     1530        END IF
     1531      END DO
     1532    END DO
     1533    DO i = 1, klon
     1534      IF (wk_adv(i)) THEN
     1535        d_sigmaw(i) = alpha(i)*d_sigmaw(i)
     1536      END IF
     1537    END DO
     1538
     1539    ! Update large scale variables and wake variables
     1540    ! IM 060208 manque DO i + remplace DO k=1,kupper(i)
     1541    ! IM 060208     DO k = 1,kupper(i)
     1542    DO k = 1, klev
     1543      DO i = 1, klon
     1544        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     1545          dtls(i, k) = dtls(i, k) + d_tb(i, k)
     1546          dqls(i, k) = dqls(i, k) + d_qb(i, k)
     1547          ! cc nrlmd
     1548          d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
     1549          d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
     1550          ! cc
     1551        END IF
     1552      END DO
     1553    END DO
     1554    DO k = 1, klev
     1555      DO i = 1, klon
     1556        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     1557          tb(i, k) = tb0(i, k) + dtls(i, k)
     1558          qb(i, k) = qb0(i, k) + dqls(i, k)
     1559          thb(i, k) = tb(i, k)/ppi(i, k)
     1560          deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
     1561          deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
     1562          dth(i, k) = deltatw(i, k)/ppi(i, k)
     1563          ! c      print*,'k,qx,qw',k,qb(i,k)-sigmaw(i)*deltaqw(i,k)
     1564          ! c     $        ,qb(i,k)+(1-sigmaw(i))*deltaqw(i,k)
     1565        END IF
     1566      END DO
     1567    END DO
     1568!
     1569    DO i = 1, klon
     1570      IF (wk_adv(i)) THEN
     1571        sigmaw(i) = sigmaw(i) + d_sigmaw(i)
     1572        d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
     1573      END IF
     1574    END DO
     1575!jyg<
     1576    IF (iflag_wk_pop_dyn >= 1) THEN
     1577!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! sigmaw !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1578!  Cumulatives
     1579      DO i = 1, klon
     1580        IF (wk_adv(i)) THEN
     1581          d_sig_gen2(i)   = d_sig_gen2(i)   + d_sig_gen(i)
     1582          d_sig_death2(i) = d_sig_death2(i) + d_sig_death(i)
     1583          d_sig_col2(i)   = d_sig_col2(i)   + d_sig_col(i)
     1584          d_sig_spread2(i)= d_sig_spread2(i)+ d_sig_spread(i)
     1585          d_sig_bnd2(i)   = d_sig_bnd2(i)   + d_sig_bnd(i)
     1586        END IF
     1587      END DO
     1588!  Bounds
     1589      DO i = 1, klon
     1590        IF (wk_adv(i)) THEN
     1591          sigmaw_targ = max(sigmaw(i),sigmad)
     1592          d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     1593          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     1594          sigmaw(i) = sigmaw_targ
     1595        END IF
     1596      END DO
     1597!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! wdens  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1598!  Cumulatives
     1599      DO i = 1, klon
     1600        IF (wk_adv(i)) THEN
     1601          wdens(i) = wdens(i) + d_wdens(i)
     1602          d_wdens2(i) = d_wdens2(i) + d_wdens(i)
     1603          d_dens_gen2(i)   = d_dens_gen2(i)   + d_dens_gen(i)
     1604          d_dens_death2(i) = d_dens_death2(i) + d_dens_death(i)
     1605          d_dens_col2(i)   = d_dens_col2(i)   + d_dens_col(i)
     1606          d_dens_bnd2(i)   = d_dens_bnd2(i)   + d_dens_bnd(i)
     1607        END IF
     1608      END DO
     1609!  Bounds
     1610      DO i = 1, klon
     1611        IF (wk_adv(i)) THEN
     1612          wdens_targ = max(wdens(i),wdensmin)
     1613          d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
     1614          d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
     1615          wdens(i) = wdens_targ
     1616        END IF
     1617      END DO
     1618!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1619!  Cumulatives
     1620      DO i = 1, klon
     1621        IF (wk_adv(i)) THEN
     1622          awdens(i) = awdens(i) + d_awdens(i)
     1623          d_awdens2(i) = d_awdens2(i) + d_awdens(i)
     1624        END IF
     1625      END DO
     1626!  Bounds
     1627      DO i = 1, klon
     1628        IF (wk_adv(i)) THEN
     1629          wdens_targ = min( max(awdens(i),0.), wdens(i) )
     1630          d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     1631          d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
     1632          awdens(i) = wdens_targ
     1633        END IF
     1634      END DO
     1635!
     1636      IF (iflag_wk_pop_dyn >= 2) THEN
     1637!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens again for iflag_wk_pop_dyn >= 2!!!!!!
     1638!  Cumulatives
     1639        DO i = 1, klon
     1640           IF (wk_adv(i)) THEN
     1641               d_adens_death2(i)   = d_adens_death2(i)   + d_adens_death(i)
     1642               d_adens_icol2(i)   = d_adens_icol2(i)   + d_adens_icol(i)
     1643               d_adens_acol2(i)   = d_adens_acol2(i)   + d_adens_acol(i)
     1644               d_adens_bnd2(i)   = d_adens_bnd2(i)   + d_adens_bnd(i)         
     1645           END IF
     1646        END DO
     1647!  Bounds
     1648        DO i = 1, klon
     1649           IF (wk_adv(i)) THEN
     1650               wdens_targ = min( max(awdens(i),0.), wdens(i) )
     1651               d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     1652               awdens(i) = wdens_targ
     1653           END IF
     1654        END DO
     1655!
     1656        IF (iflag_wk_pop_dyn == 3) THEN
     1657!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! asigmaw for iflag_wk_pop_dyn = 3!!!!!!
     1658!  Cumulatives
     1659          DO i = 1, klon
     1660             IF (wk_adv(i)) THEN
     1661                 asigmaw(i) = asigmaw(i) + d_asigmaw(i)
     1662                 d_asigmaw2(i) = d_asigmaw2(i) + d_asigmaw(i)
     1663                 d_asig_death2(i)   = d_asig_death2(i)   + d_asig_death(i)
     1664                 d_asig_spread2(i)  = d_asig_spread2(i)  + d_asig_spread(i)
     1665                 d_asig_iicol2(i)   = d_asig_iicol2(i)   + d_asig_iicol(i)
     1666                 d_asig_aicol2(i)   = d_asig_aicol2(i)   + d_asig_aicol(i)
     1667                 d_asig_bnd2(i)     = d_asig_bnd2(i)     + d_asig_bnd(i)         
     1668             END IF
     1669          END DO
     1670!  Bounds
     1671          DO i = 1, klon
     1672             IF (wk_adv(i)) THEN
     1673   !   asigmaw lower bound set to sigmad/2 in order to allow asigmaw values lower than sigmad.
     1674   !!             sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i))
     1675                sigmaw_targ = min(max(asigmaw(i),sigmad/2.),sigmaw(i))
     1676                d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     1677                d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     1678                asigmaw(i) = sigmaw_targ
     1679             END IF
     1680          END DO
     1681
     1682IF (CPPKEY_IOPHYS_WK) THEN
     1683    IF (phys_sub) THEN
     1684     CALL iophys_ecrit('wdensb',1,'wdensb','m',wdens)
     1685     CALL iophys_ecrit('awdensb',1,'awdensb','m',awdens)
     1686     CALL iophys_ecrit('sigmawb',1,'sigmawb','m',sigmaw)
     1687     CALL iophys_ecrit('asigmawb',1,'asigmawb','m',asigmaw)
     1688!
     1689     call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2)
     1690     call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2)
     1691     call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2)
     1692     call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2)
     1693     call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2)
     1694!
     1695     call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2)
     1696     call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2)
     1697     call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2)
     1698     call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2)
     1699     call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2)
     1700!
     1701     CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2)
     1702     CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2)
     1703     CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2)
     1704     CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2)
     1705     CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2)
     1706     CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2)
     1707!
     1708     CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2)
     1709     CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2)
     1710     CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2)
     1711     CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2)
     1712     CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2)
     1713     CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2)
     1714    ENDIF
     1715END IF
     1716        ENDIF ! (iflag_wk_pop_dyn == 3)
     1717      ENDIF ! (iflag_wk_pop_dyn >= 2)
     1718    ENDIF  ! (iflag_wk_pop_dyn >= 1)
     1719
     1720
     1721
     1722   Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
     1723                    dth, hw, rho, delta_t_min, &
     1724                    ktop, wk_adv, h_zzz, ptop1, ktop1)
     1725   !! print'("pkupper APPEL ",7i6)',isubstep,int(ptop/100.),int(ptop1/100.),int(pupper/100.),ktop,ktop1,kupper
     1726
     1727    ! 5/ Set deltatw & deltaqw to 0 above kupper
     1728
     1729    DO k = 1, klev
     1730      DO i = 1, klon
     1731        IF (wk_adv(i) .AND. k>=kupper(i)) THEN
     1732          deltatw(i, k) = 0.
     1733          deltaqw(i, k) = 0.
     1734          d_deltatw2(i,k) = -deltatw0(i,k)
     1735          d_deltaqw2(i,k) = -deltaqw0(i,k)
     1736        END IF
     1737      END DO
     1738    END DO
     1739
     1740
     1741    ! -------------Cstar computation---------------------------------
     1742    DO i = 1, klon
     1743      IF (wk_adv(i)) THEN !!! nrlmd
     1744        sum_thx(i) = 0.
     1745        sum_tx(i) = 0.
     1746        sum_qx(i) = 0.
     1747        sum_thvx(i) = 0.
     1748        sum_dth(i) = 0.
     1749        sum_dq(i) = 0.
     1750        sum_dtdwn(i) = 0.
     1751        sum_dqdwn(i) = 0.
     1752
     1753        av_thx(i) = 0.
     1754        av_tx(i) = 0.
     1755        av_qx(i) = 0.
     1756        av_thvx(i) = 0.
     1757        av_dth(i) = 0.
     1758        av_dq(i) = 0.
     1759        av_dtdwn(i) = 0.
     1760        av_dqdwn(i) = 0.
     1761      END IF
     1762    END DO
     1763
     1764    ! Integrals (and wake top level number)
     1765    ! --------------------------------------
     1766
     1767    ! Initialize sum_thvx to 1st level virt. pot. temp.
     1768
     1769    DO i = 1, klon
     1770      IF (wk_adv(i)) THEN !!! nrlmd
     1771        z(i) = 1.
     1772        dz(i) = 1.
     1773        sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i)
     1774        sum_dth(i) = 0.
     1775      END IF
     1776    END DO
     1777
     1778    DO k = 1, klev
     1779      DO i = 1, klon
     1780        IF (wk_adv(i)) THEN !!! nrlmd
     1781          dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
     1782          IF (dz(i)>0) THEN
     1783            z(i) = z(i) + dz(i)
     1784            sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i)
     1785            sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i)
     1786            sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i)
     1787            sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i)
     1788            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
     1789            sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
     1790            sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
     1791            sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
     1792          END IF
     1793        END IF
     1794      END DO
     1795    END DO
     1796
     1797    DO i = 1, klon
     1798      IF (wk_adv(i)) THEN !!! nrlmd
     1799        hw0(i) = z(i)
     1800      END IF
     1801    END DO
     1802
     1803
     1804    ! - WAPE and mean forcing computation
     1805    ! ---------------------------------------
     1806
     1807    ! ---------------------------------------
     1808
     1809    ! Means
     1810
     1811    DO i = 1, klon
     1812      IF (wk_adv(i)) THEN !!! nrlmd
     1813        av_thx(i) = sum_thx(i)/hw0(i)
     1814        av_tx(i) = sum_tx(i)/hw0(i)
     1815        av_qx(i) = sum_qx(i)/hw0(i)
     1816        av_thvx(i) = sum_thvx(i)/hw0(i)
     1817        av_dth(i) = sum_dth(i)/hw0(i)
     1818        av_dq(i) = sum_dq(i)/hw0(i)
     1819        av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     1820        av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     1821
     1822        wape(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_thx(i)*av_dq(i) + &
     1823                              av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i)
     1824      END IF
     1825    END DO
     1826
     1827
     1828    ! Filter out bad wakes
     1829
     1830    DO k = 1, klev
     1831      DO i = 1, klon
     1832        IF (wk_adv(i)) THEN !!! nrlmd
     1833          IF (wape(i)<0.) THEN
     1834            deltatw(i, k) = 0.
     1835            deltaqw(i, k) = 0.
     1836            dth(i, k) = 0.
     1837            d_deltatw2(i,k) = -deltatw0(i,k)
     1838            d_deltaqw2(i,k) = -deltaqw0(i,k)
     1839          END IF
     1840        END IF
     1841      END DO
     1842    END DO
     1843
     1844    DO i = 1, klon
     1845      IF (wk_adv(i)) THEN !!! nrlmd
     1846        IF (wape(i)<0.) THEN
     1847          wape(i) = 0.
     1848          cstar(i) = 0.
     1849          hw(i) = hwmin
     1850!jyg<
     1851!!          sigmaw(i) = max(sigmad, sigd_con(i))
     1852          sigmaw_targ = max(sigmad, sigd_con(i))
     1853          d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     1854          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     1855          sigmaw(i) = sigmaw_targ
     1856!
     1857          d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     1858          d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     1859          asigmaw(i) = sigmaw_targ
     1860!>jyg
     1861          fip(i) = 0.
     1862          gwake(i) = .FALSE.
     1863        ELSE
     1864          cstar(i) = stark*sqrt(2.*wape(i))
     1865          gwake(i) = .TRUE.
     1866        END IF
     1867      END IF
     1868    END DO
     1869  !
     1870  ! ------------------------------------------------------------------------
     1871  !
     1872  END DO   ! isubstep end sub-timestep loop
     1873  !
     1874  ! ------------------------------------------------------------------------
     1875  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     1876  ! ------------------------------------------------------------------------
     1877  !
     1878
     1879IF (CPPKEY_IOPHYS_WK) THEN
     1880    IF (.not.phys_sub) CALL iophys_ecrit('wape_b',1,'wape_b','J/kg',wape)
     1881END IF
     1882  IF (prt_level>=10) THEN
     1883    PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', &
     1884                      sigmaw(igout), cstar(igout), wape(igout), ptop(igout)
     1885  ENDIF
     1886
     1887
     1888  ! ----------------------------------------------------------
     1889  ! Determine wake final state; recompute wape, cstar, ktop;
     1890  ! filter out bad wakes.
     1891  ! ----------------------------------------------------------
     1892
     1893  ! 2.1 - Undisturbed area and Wake integrals
     1894  ! ---------------------------------------------------------
     1895
     1896  DO i = 1, klon
     1897    ! cc nrlmd       if (wk_adv(i)) then !!! nrlmd
     1898    IF (ok_qx_qw(i)) THEN
     1899      ! cc
     1900      z(i) = 0.
     1901      sum_thx(i) = 0.
     1902      sum_tx(i) = 0.
     1903      sum_qx(i) = 0.
     1904      sum_thvx(i) = 0.
     1905      sum_dth(i) = 0.
     1906      sum_half_dth(i) = 0.
     1907      sum_dq(i) = 0.
     1908      sum_dtdwn(i) = 0.
     1909      sum_dqdwn(i) = 0.
     1910
     1911      av_thx(i) = 0.
     1912      av_tx(i) = 0.
     1913      av_qx(i) = 0.
     1914      av_thvx(i) = 0.
     1915      av_dth(i) = 0.
     1916      av_dq(i) = 0.
     1917      av_dtdwn(i) = 0.
     1918      av_dqdwn(i) = 0.
     1919
     1920      dthmin(i) = -delta_t_min
     1921    END IF
     1922  END DO
     1923  ! Potential temperatures and humidity
     1924  ! ----------------------------------------------------------
     1925
     1926  DO k = 1, klev
     1927    DO i = 1, klon
     1928      ! cc nrlmd       IF ( wk_adv(i)) THEN
     1929      IF (ok_qx_qw(i)) THEN
     1930        ! cc
     1931        rho(i, k) = p(i, k)/(RD*tb(i,k))
     1932        IF (k==1) THEN
     1933          rhoh(i, k) = ph(i, k)/(RD*tb(i,k))
     1934          zhh(i, k) = 0
     1935        ELSE
     1936          rhoh(i, k) = ph(i, k)*2./(RD*(tb(i,k)+tb(i,k-1)))
     1937          zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
     1938        END IF
     1939        thb(i, k) = tb(i, k)/ppi(i, k)
     1940        thx(i, k) = (tb(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
     1941        tx(i, k) = tb(i, k) - deltatw(i, k)*sigmaw(i)
     1942        qx(i, k) = qb(i, k) - deltaqw(i, k)*sigmaw(i)
     1943        dth(i, k) = deltatw(i, k)/ppi(i, k)
     1944      END IF
     1945    END DO
     1946  END DO
     1947
     1948  ! Integrals (and wake top level number)
     1949  ! -----------------------------------------------------------
     1950
     1951  ! Initialize sum_thvx to 1st level virt. pot. temp.
     1952
     1953  DO i = 1, klon
     1954    ! cc nrlmd       IF ( wk_adv(i)) THEN
     1955    IF (ok_qx_qw(i)) THEN
     1956      ! cc
     1957      z(i) = 1.
     1958      dz(i) = 1.
     1959      dz_half(i) = 1.
     1960      sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i)
     1961      sum_dth(i) = 0.
     1962    END IF
     1963  END DO
     1964
     1965  DO k = 1, klev
     1966    DO i = 1, klon
     1967      ! cc nrlmd       IF ( wk_adv(i)) THEN
     1968      IF (ok_qx_qw(i)) THEN
     1969        ! cc
     1970        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
     1971        dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)*RG)
     1972        IF (dz(i)>0) THEN
     1973          z(i) = z(i) + dz(i)
     1974          sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i)
     1975          sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i)
     1976          sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i)
     1977          sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i)
     1978          sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
     1979          sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
     1980          sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
     1981          sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
     1982!
     1983          dthmin(i) = min(dthmin(i), dth(i,k))
     1984        END IF
     1985        IF (dz_half(i)>0) THEN
     1986          sum_half_dth(i) = sum_half_dth(i) + dth(i, k)*dz_half(i)
     1987        END IF
     1988      END IF
     1989    END DO
     1990  END DO
     1991
     1992  DO i = 1, klon
     1993    ! cc nrlmd       IF ( wk_adv(i)) THEN
     1994    IF (ok_qx_qw(i)) THEN
     1995      ! cc
     1996      hw0(i) = z(i)
     1997    END IF
     1998  END DO
     1999
     2000  ! - WAPE and mean forcing computation
     2001  ! -------------------------------------------------------------
     2002
     2003  ! Means
     2004
     2005  DO i = 1, klon
     2006    ! cc nrlmd       IF ( wk_adv(i)) THEN
     2007    IF (ok_qx_qw(i)) THEN
     2008      ! cc
     2009      av_thx(i) = sum_thx(i)/hw0(i)
     2010      av_tx(i) = sum_tx(i)/hw0(i)
     2011      av_qx(i) = sum_qx(i)/hw0(i)
     2012      av_thvx(i) = sum_thvx(i)/hw0(i)
     2013      av_dth(i) = sum_dth(i)/hw0(i)
     2014      av_dq(i) = sum_dq(i)/hw0(i)
     2015      av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     2016      av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     2017
     2018      wape2(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_thx(i)*av_dq(i) + &
     2019                             av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i)
     2020    END IF
     2021  END DO
     2022IF (CPPKEY_IOPHYS_WK) THEN
     2023  IF (.not.phys_sub) CALL iophys_ecrit('wape2_a',1,'wape2_a','J/kg',wape2)
     2024END IF
     2025
     2026
     2027  ! Prognostic variable update
     2028  ! ------------------------------------------------------------
     2029
     2030  ! Filter out bad wakes
     2031
     2032  IF (iflag_wk_check_trgl>=1) THEN
     2033    ! Check triangular shape of dth profile
     2034    DO i = 1, klon
     2035      IF (ok_qx_qw(i)) THEN
     2036        !! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
     2037        !! print *,'wake, 2.*sum_dth(i)/(hw0(i)*dthmin(i)) ', &
     2038        !!                2.*sum_dth(i)/(hw0(i)*dthmin(i))
     2039        !! print *,'wake, sum_half_dth(i), sum_dth(i) ', &
     2040        !!                sum_half_dth(i), sum_dth(i)
     2041        IF ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min) ) THEN
     2042          wape2(i) = -1.
     2043          !! print *,'wake, rej 1'
     2044        ELSE IF (iflag_wk_check_trgl==1.AND.abs(2.*sum_dth(i)/(hw0(i)*dthmin(i)) - 1.) > 0.5) THEN
     2045          wape2(i) = -1.
     2046          !! print *,'wake, rej 2'
     2047        ELSE IF (abs(sum_half_dth(i)) < 0.5*abs(sum_dth(i)) ) THEN
     2048          wape2(i) = -1.
     2049          !! print *,'wake, rej 3'
     2050        END IF
     2051      END IF
     2052    END DO
     2053  END IF
     2054IF (CPPKEY_IOPHYS_WK) THEN
     2055  IF (.not.phys_sub) CALL iophys_ecrit('wape2_b',1,'wape2_b','J/kg',wape2)
     2056END IF
     2057
     2058
     2059  DO k = 1, klev
     2060    DO i = 1, klon
     2061      ! cc nrlmd        IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
     2062      IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
     2063        ! cc
     2064        deltatw(i, k) = 0.
     2065        deltaqw(i, k) = 0.
     2066        dth(i, k) = 0.
     2067        d_deltatw2(i,k) = -deltatw0(i,k)
     2068        d_deltaqw2(i,k) = -deltaqw0(i,k)
     2069      END IF
     2070    END DO
     2071  END DO
     2072
     2073
     2074  DO i = 1, klon
     2075    ! cc nrlmd       IF ( wk_adv(i)) THEN
     2076    IF (ok_qx_qw(i)) THEN
     2077      ! cc
     2078      IF (wape2(i)<0.) THEN
     2079        wape2(i) = 0.
     2080        cstar2(i) = 0.
     2081        hw(i) = hwmin
     2082!jyg<
     2083!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
     2084      sigmaw_targ = max(sigmad, sigd_con(i))
     2085      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     2086      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     2087      sigmaw(i) = sigmaw_targ
     2088!
     2089      d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     2090      d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     2091      asigmaw(i) = sigmaw_targ
     2092!>jyg
     2093        fip(i) = 0.
     2094        gwake(i) = .FALSE.
     2095      ELSE
     2096        IF (prt_level>=10) PRINT *, 'wape2>0'
     2097        cstar2(i) = stark*sqrt(2.*wape2(i))
     2098        gwake(i) = .TRUE.
     2099      END IF
     2100IF (CPPKEY_IOPHYS_WK) THEN
     2101  IF (.not.phys_sub) CALL iophys_ecrit('cstar2',1,'cstar2','J/kg',cstar2)
     2102END IF
     2103    END IF  ! (ok_qx_qw(i))
     2104  END DO
     2105
     2106  DO i = 1, klon
     2107    ! cc nrlmd       IF ( wk_adv(i)) THEN
     2108    IF (ok_qx_qw(i)) THEN
     2109      ! cc
     2110      ktopw(i) = ktop(i)
     2111    END IF
     2112  END DO
     2113
     2114  DO i = 1, klon
     2115    ! cc nrlmd       IF ( wk_adv(i)) THEN
     2116    IF (ok_qx_qw(i)) THEN
     2117      ! cc
     2118      IF (ktopw(i)>0 .AND. gwake(i)) THEN
     2119
     2120        ! jyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
     2121        ! cc       heff = 600.
     2122        ! Utilisation de la hauteur hw
     2123        ! c       heff = 0.7*hw
     2124        heff(i) = hw(i)
     2125
     2126        fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* &
     2127          sqrt(sigmaw(i)*wdens(i)*3.14)
     2128        fip(i) = alpk*fip(i)
     2129        ! jyg2
     2130      ELSE
     2131        fip(i) = 0.
     2132      END IF
     2133    END IF
     2134  END DO
     2135    IF (iflag_wk_pop_dyn >= 3) THEN
     2136IF (CPPKEY_IOPHYS_WK) THEN
     2137      IF (.not.phys_sub) THEN
     2138       CALL iophys_ecrit('fip',1,'fip','J/kg',fip)
     2139       CALL iophys_ecrit('hw',1,'hw','J/kg',hw)
     2140       CALL iophys_ecrit('ptop',1,'ptop','J/kg',ptop)
     2141       CALL iophys_ecrit('wdens',1,'wdens','J/kg',wdens)
     2142       CALL iophys_ecrit('awdens',1,'awdens','m',awdens)
     2143       CALL iophys_ecrit('sigmaw',1,'sigmaw','m',sigmaw)
     2144       CALL iophys_ecrit('asigmaw',1,'asigmaw','m',asigmaw)
     2145!   
     2146       CALL iophys_ecrit('rad_wk',1,'rad_wk','J/kg',rad_wk)
     2147       CALL iophys_ecrit('arad_wk',1,'arad_wk','J/kg',arad_wk)
     2148       CALL iophys_ecrit('irad_wk',1,'irad_wk','J/kg',irad_wk)
     2149!   
     2150       call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2)
     2151       call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2)
     2152       call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2)
     2153       call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2)
     2154       call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2)
     2155!   
     2156       call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2)
     2157       call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2)
     2158       call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2)
     2159       call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2)
     2160       call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2)
     2161!   
     2162       CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2)
     2163       CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2)
     2164       CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2)
     2165       CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2)
     2166       CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2)
     2167       CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2)
     2168!   
     2169       CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2)
     2170       CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2)
     2171       CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2)
     2172       CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2)
     2173       CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2)
     2174       CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2)
     2175      ENDIF  ! (.not.phys_sub)
     2176END IF
     2177    ENDIF  ! (iflag_wk_pop_dyn >= 3)
     2178  ! Limitation de sigmaw
     2179
     2180  ! cc nrlmd
     2181  ! DO i=1,klon
     2182  ! IF (OK_qx_qw(i)) THEN
     2183  ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
     2184  ! ENDIF
     2185  ! ENDDO
     2186  ! cc
     2187
     2188  !jyg<
     2189  IF (iflag_wk_pop_dyn >= 1) THEN
     2190    DO i = 1, klon
     2191      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     2192          .NOT. ok_qx_qw(i) .OR. (wdens(i) < wdensthreshold)
     2193!!          .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
     2194    ENDDO
     2195  ELSE  ! (iflag_wk_pop_dyn >= 1)
     2196    DO i = 1, klon
     2197      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     2198          .NOT. ok_qx_qw(i)
     2199    ENDDO
     2200  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     2201  !>jyg
     2202
     2203  DO k = 1, klev
     2204    DO i = 1, klon
     2205!!jyg      IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     2206!!jyg          .NOT. ok_qx_qw(i)) THEN
     2207      IF (kill_wake(i)) THEN
     2208        ! cc
     2209        dtls(i, k) = 0.
     2210        dqls(i, k) = 0.
     2211        deltatw(i, k) = 0.
     2212        deltaqw(i, k) = 0.
     2213        d_deltatw2(i,k) = -deltatw0(i,k)
     2214        d_deltaqw2(i,k) = -deltaqw0(i,k)
     2215      END IF  ! (kill_wake(i))
     2216    END DO
     2217  END DO
     2218
     2219  DO i = 1, klon
     2220!!jyg    IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     2221!!jyg        .NOT. ok_qx_qw(i)) THEN
     2222      IF (kill_wake(i)) THEN
     2223      ktopw(i) = 0
     2224      wape(i) = 0.
     2225      cstar(i) = 0.
     2226!!jyg   Outside subroutine "Wake" hw, wdens sigmaw and asigmaw are zero when there are no wakes
     2227!!      hw(i) = hwmin                       !jyg
     2228!!      sigmaw(i) = sigmad                  !jyg
     2229      hw(i) = 0.                            !jyg
     2230      fip(i) = 0.
     2231!
     2232!!      sigmaw(i) = 0.                        !jyg
     2233      sigmaw_targ = 0.
     2234      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     2235!!      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     2236      d_sigmaw2(i) = sigmaw_targ - sigmaw_in(i)      ! _in = correction jyg 20220124
     2237      sigmaw(i) = sigmaw_targ
     2238!
     2239      IF (iflag_wk_pop_dyn >= 3) THEN
     2240        sigmaw_targ = 0.
     2241        d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     2242!!        d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     2243        d_asigmaw2(i) = sigmaw_targ - asigmaw_in(i)      ! _in = correction jyg 20220124
     2244        asigmaw(i) = sigmaw_targ
     2245      ELSE
     2246        asigmaw(i) = 0.
     2247      ENDIF ! (iflag_wk_pop_dyn >= 3)
     2248!
     2249      IF (iflag_wk_pop_dyn >= 1) THEN
     2250!!        awdens(i) = 0.
     2251!!        wdens(i) = 0.
     2252        wdens_targ = 0.
     2253        d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
     2254!!        d_wdens2(i) = wdens_targ - wdens(i)
     2255        d_wdens2(i) = wdens_targ - wdens_in(i)      ! jyg 20220916
     2256        wdens(i) = wdens_targ
     2257        wdens_targ = 0.
     2258!!jyg: bug fix : the d_adens_bnd2 computation must be before the update of awdens.
     2259        IF (iflag_wk_pop_dyn >= 2) THEN
     2260            d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     2261        ENDIF ! (iflag_wk_pop_dyn >= 2)
     2262!!        d_awdens2(i) = wdens_targ - awdens(i)
     2263        d_awdens2(i) = wdens_targ - awdens_in(i)    ! jyg 20220916
     2264        awdens(i) = wdens_targ
     2265!!        IF (iflag_wk_pop_dyn == 2) THEN
     2266!!            d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     2267!!        ENDIF ! (iflag_wk_pop_dyn == 2)
     2268      ENDIF  ! (iflag_wk_pop_dyn >= 1)
     2269    ELSE  ! (kill_wake(i))
     2270      wape(i) = wape2(i)
     2271      cstar(i) = cstar2(i)
     2272    END IF  ! (kill_wake(i))
     2273    ! c        print*,'wape wape2 ktopw OK_qx_qw =',
     2274    ! c     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
     2275  END DO
     2276
     2277  IF (prt_level>=10) THEN
     2278    PRINT *, 'wake-6, wape wape2 ktopw OK_qx_qw =', &
     2279                      wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout)
     2280  ENDIF
     2281IF (CPPKEY_IOPHYS_WK) THEN
     2282  IF (.not.phys_sub) CALL iophys_ecrit('wape_c',1,'wape_c','J/kg',wape)
     2283END IF
     2284
     2285
     2286  ! -----------------------------------------------------------------
     2287  ! Get back to tendencies per second
     2288
     2289  DO k = 1, klev
     2290    DO i = 1, klon
     2291
     2292      ! cc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
     2293!jyg<
     2294!!      IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
     2295      IF (ok_qx_qw(i)) THEN
     2296!>jyg
     2297        ! cc
     2298        dtls(i, k) = dtls(i, k)/dtime
     2299        dqls(i, k) = dqls(i, k)/dtime
     2300        d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
     2301        d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
     2302        d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
     2303        ! c      print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
     2304        ! c     $         ,death_rate(i)*sigmaw(i)
     2305      END IF
     2306    END DO
     2307  END DO
     2308!jyg<
     2309  IF (iflag_wk_pop_dyn >= 1) THEN
     2310    DO i = 1, klon
     2311        IF (ok_qx_qw(i)) THEN
     2312      d_sig_gen2(i) = d_sig_gen2(i)/dtime
     2313      d_sig_death2(i) = d_sig_death2(i)/dtime
     2314      d_sig_col2(i) = d_sig_col2(i)/dtime
     2315      d_sig_spread2(i) = d_sig_spread2(i)/dtime
     2316      d_sig_bnd2(i) = d_sig_bnd2(i)/dtime
     2317      d_sigmaw2(i) = d_sigmaw2(i)/dtime
     2318!
     2319      d_dens_gen2(i) = d_dens_gen2(i)/dtime
     2320      d_dens_death2(i) = d_dens_death2(i)/dtime
     2321      d_dens_col2(i) = d_dens_col2(i)/dtime
     2322      d_dens_bnd2(i) = d_dens_bnd2(i)/dtime
     2323      d_awdens2(i) = d_awdens2(i)/dtime
     2324      d_wdens2(i) = d_wdens2(i)/dtime
     2325        ENDIF
     2326    ENDDO
     2327    IF (iflag_wk_pop_dyn >= 2) THEN
     2328      DO i = 1, klon
     2329        IF (ok_qx_qw(i)) THEN
     2330        d_adens_death2(i) = d_adens_death2(i)/dtime
     2331        d_adens_icol2(i) = d_adens_icol2(i)/dtime
     2332        d_adens_acol2(i) = d_adens_acol2(i)/dtime
     2333        d_adens_bnd2(i) = d_adens_bnd2(i)/dtime
     2334        ENDIF
     2335      ENDDO
     2336      IF (iflag_wk_pop_dyn == 3) THEN
     2337       DO i = 1, klon
     2338          IF (ok_qx_qw(i)) THEN
     2339        d_asig_death2(i)  = d_asig_death2(i)/dtime
     2340        d_asig_iicol2(i)  = d_asig_iicol2(i)/dtime
     2341        d_asig_aicol2(i)  = d_asig_aicol2(i)/dtime
     2342        d_asig_spread2(i) = d_asig_spread2(i)/dtime
     2343        d_asig_bnd2(i) = d_asig_bnd2(i)/dtime
     2344          ENDIF
     2345       ENDDO
     2346      ENDIF ! (iflag_wk_pop_dyn == 3) 
     2347    ENDIF ! (iflag_wk_pop_dyn >= 2) 
     2348  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     2349 
     2350!>jyg
     2351
     2352 RETURN
     2353END SUBROUTINE wake
     2354
     2355SUBROUTINE wake2(klon,klev,znatsurf, p, ph, pi, dtime, &
     2356                tb0, qb0, omgb, &
     2357                dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
     2358                sigd_con, Cin, &
     2359                deltatw, deltaqw, sigmaw, asigmaw, wdens, awdens, &                  ! state variables           
     2360                dth, hw, wape, fip, gfl, &
     2361                dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, &
     2362!!!                dtke, dqke, omg, dp_deltomg, wkspread, cstar, &   ! changes in notation
     2363                d_deltat_dcv, d_deltaq_dcv, omg, dp_deltomg, wkspread, cstar, &
     2364                d_deltat_gw, &                                                      ! tendencies
     2365                d_deltatw2, d_deltaqw2, d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2)             ! tendencies
     2366
     2367
     2368  ! **************************************************************
     2369  ! *
     2370  ! WAKE                                                        *
     2371  ! retour a un Pupper fixe                                *
     2372  ! *
     2373  ! written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
     2374  ! modified by :   ROEHRIG Romain        01/29/2007            *
     2375  ! **************************************************************
     2376
     2377
     2378  USE lmdz_wake_ini , ONLY : wake_ini
     2379  USE lmdz_wake_ini , ONLY : prt_level,epsim1,RG,RD
     2380  USE lmdz_wake_ini , ONLY : stark, wdens_ref, coefgw, alpk, wk_pupper
     2381  USE lmdz_wake_ini , ONLY : crep_upper, crep_sol, tau_cv, rzero, aa0, flag_wk_check_trgl
     2382  USE lmdz_wake_ini , ONLY : ok_bug_gfl
     2383  USE lmdz_wake_ini , ONLY : iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensinit, wdensthreshold
     2384  USE lmdz_wake_ini , ONLY : sigmad, hwmin, wapecut, cstart, sigmaw_max, dens_rate, epsilon_loc
     2385  USE lmdz_wake_ini , ONLY : iflag_wk_profile
     2386  USE lmdz_wake_ini , ONLY : smallestreal,wk_nsub
     2387
     2388
     2389  IMPLICIT NONE
     2390  ! ============================================================================
     2391
     2392
     2393  ! But : Decrire le comportement des poches froides apparaissant dans les
     2394  ! grands systemes convectifs, et fournir l'energie disponible pour
     2395  ! le declenchement de nouvelles colonnes convectives.
     2396
     2397  ! State variables :
     2398  ! deltatw    : temperature difference between wake and off-wake regions
     2399  ! deltaqw    : specific humidity difference between wake and off-wake regions
     2400  ! sigmaw     : fractional area covered by wakes.
     2401  ! asigmaw    : fractional area covered by active wakes.
     2402  ! wdens      : number of wakes per unit area
     2403  ! awdens     : number of active wakes per unit area
     2404
     2405  ! Variable de sortie :
     2406
     2407  ! wape : WAke Potential Energy
     2408  ! fip  : Front Incident Power (W/m2) - ALP
     2409  ! gfl  : Gust Front Length per unit area (m-1)
     2410  ! dtls : large scale temperature tendency due to wake
     2411  ! dqls : large scale humidity tendency due to wake
     2412  ! hw   : wake top hight (given by hw*deltatw(1)/2=wape)
     2413  ! dp_omgb : vertical gradient of large scale omega
     2414  ! awdens  : densite de poches actives
     2415  ! wdens   : densite de poches
     2416  ! omgbdth: flux of Delta_Theta transported by LS omega
     2417  ! d_deltat_dcv   : differential heating (wake - unpertubed)
     2418  ! d_deltat_dcv   : differential moistening (wake - unpertubed)
     2419  ! omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
     2420  ! dp_deltomg  : vertical gradient of omg (s-1)
     2421  ! wkspread  : spreading term in d_t_wake and d_q_wake
     2422  ! deltatw     : updated temperature difference (T_w-T_u).
     2423  ! deltaqw     : updated humidity difference (q_w-q_u).
     2424  ! sigmaw      : updated wake fractional area.
     2425  ! asigmaw     : updated active wake fractional area.
     2426  ! d_deltat_gw : delta T tendency due to GW
     2427
     2428  ! Variables d'entree :
     2429
     2430  ! aire : aire de la maille
     2431  ! tb0  : horizontal average of temperature  (K)
     2432  ! qb0  : horizontal average of humidity   (kg/kg)
     2433  ! omgb : vitesse verticale moyenne sur la maille (Pa/s)
     2434  ! dtdwn: source de chaleur due aux descentes (K/s)
     2435  ! dqdwn: source d'humidite due aux descentes (kg/kg/s)
     2436  ! dta  : source de chaleur due courants satures et detrain  (K/s)
     2437  ! dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
     2438  ! wgen : number of wakes generated per unit area and per sec (/m^2/s)
     2439  ! amdwn: flux de masse total des descentes, par unite de
     2440  !        surface de la maille (kg/m2/s)
     2441  ! amup : flux de masse total des ascendances, par unite de
     2442  !        surface de la maille (kg/m2/s)
     2443  ! sigd_con:
     2444  ! Cin  : convective inhibition
     2445  ! p    : pressions aux milieux des couches (Pa)
     2446  ! ph   : pressions aux interfaces (Pa)
     2447  ! pi  : (p/p_0)**kapa (adim)
     2448  ! dtime: increment temporel (s)
     2449
     2450  ! Variables internes :
     2451
     2452  ! rho  : mean density at P levels
     2453  ! rhoh : mean density at Ph levels
     2454  ! tb   : mean temperature | may change within
     2455  ! qb   : mean humidity    | sub-time-stepping
     2456  ! thb  : mean potential temperature
     2457  ! thx  : potential temperature in (x) area
     2458  ! tx   : temperature  in (x) area
     2459  ! qx   : humidity in (x) area
     2460  ! dp_omgb: vertical gradient og LS omega
     2461  ! omgbw  : wake average vertical omega
     2462  ! dp_omgbw: vertical gradient of omgbw
     2463  ! omgbdq : flux of Delta_q transported by LS omega
     2464  ! dth  : potential temperature diff. wake-undist.
     2465  ! th1  : first pot. temp. for vertical advection (=thx)
     2466  ! th2  : second pot. temp. for vertical advection (=thw)
     2467  ! q1   : first humidity for vertical advection
     2468  ! q2   : second humidity for vertical advection
     2469  ! d_deltatw   : redistribution term for deltatw
     2470  ! d_deltaqw   : redistribution term for deltaqw
     2471  ! deltatw0   : initial deltatw
     2472  ! deltaqw0   : initial deltaqw
     2473  ! hw0    : wake top hight (defined as the altitude at which deltatw=0)
     2474  ! amflux : horizontal mass flux through wake boundary
     2475  ! wdens_ref: initial number of wakes per unit area (3D) or per
     2476  !            unit length (2D), at the beginning of each time step
     2477  ! Tgw    : 1 sur la periode de onde de gravite
     2478  ! Cgw    : vitesse de propagation de onde de gravite
     2479  ! LL     : distance between 2 wakes
     2480
     2481  ! -------------------------------------------------------------------------
     2482  ! Declaration de variables
     2483  ! -------------------------------------------------------------------------
     2484
     2485
     2486  ! Arguments en entree
     2487  ! --------------------
     2488
     2489  INTEGER,                          INTENT(IN)          :: klon,klev
     2490  INTEGER, DIMENSION (klon),        INTENT(IN)          :: znatsurf
     2491  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: p, pi
     2492  REAL, DIMENSION (klon, klev+1),   INTENT(IN)          :: ph
     2493  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: omgb
     2494  REAL,                             INTENT(IN)          :: dtime
     2495  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: tb0, qb0
     2496  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dtdwn, dqdwn
     2497  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: amdwn, amup
     2498  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dta, dqa
     2499  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen
     2500  REAL, DIMENSION (klon),           INTENT(IN)          :: sigd_con
     2501  REAL, DIMENSION (klon),           INTENT(IN)          :: Cin
     2502
     2503  !
     2504  ! Input/Output
     2505  ! State variables
     2506  REAL, DIMENSION (klon, klev),     INTENT(INOUT)       :: deltatw, deltaqw
     2507  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw
     2508  REAL, DIMENSION (klon),           INTENT(INOUT)       :: asigmaw
     2509  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens
     2510  REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens
     2511
     2512  ! Sorties
     2513  ! --------
     2514
     2515  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dth
     2516  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: tx, qx
     2517  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls
     2518!!  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
     2519  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_dcv, d_deltaq_dcv
     2520  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: wkspread    !  unused (jyg)
     2521  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: omgbdth, omg
     2522  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dp_omgb, dp_deltomg
     2523  REAL, DIMENSION (klon),           INTENT(OUT)         :: hw, wape, fip, gfl, cstar
     2524  INTEGER, DIMENSION (klon),        INTENT(OUT)         :: ktopw
     2525  ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields
     2526  !                                 computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
     2527  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw
     2528  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
     2529  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2
     2530
     2531  ! Variables internes
     2532  ! -------------------
     2533
     2534  ! Variables a fixer
     2535
     2536  REAL                                                  :: delta_t_min
     2537  REAL                                                  :: dtimesub
     2538  REAL                                                  :: wdens0
     2539  ! IM 080208
     2540  LOGICAL, DIMENSION (klon)                             :: gwake
     2541
     2542  ! Variables de sauvegarde
     2543  REAL, DIMENSION (klon, klev)                          :: deltatw0
     2544  REAL, DIMENSION (klon, klev)                          :: deltaqw0
     2545  REAL, DIMENSION (klon, klev)                          :: tb, qb
     2546
     2547  ! Variables liees a la dynamique de population 1
     2548  REAL, DIMENSION(klon)                                 :: act
     2549  REAL, DIMENSION(klon)                                 :: rad_wk, tau_wk_inv
     2550  REAL, DIMENSION(klon)                                 :: f_shear
     2551  REAL, DIMENSION(klon)                                 :: drdt
     2552 
     2553  ! Variables liees a la dynamique de population 2
     2554  REAL, DIMENSION(klon)                                 :: cont_fact 
     2555 
     2556  ! Variables liees a la dynamique de population 3
     2557  REAL, DIMENSION(klon)                                 :: arad_wk, irad_wk
     2558 
     2559!!  REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col
     2560  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
     2561  LOGICAL, DIMENSION (klon)                             :: kill_wake
     2562  REAL                                                  :: drdt_pos
     2563  REAL                                                  :: tau_wk_inv_min
     2564  ! Some components of the tendencies of state variables 
     2565  REAL, DIMENSION (klon)                                :: d_sig_gen2, d_sig_death2, d_sig_col2, d_sig_spread2, d_sig_bnd2
     2566  REAL, DIMENSION (klon)                                :: d_asig_death2, d_asig_aicol2, d_asig_iicol2, d_asig_spread2, d_asig_bnd2
     2567  REAL, DIMENSION (klon)                                :: d_dens_gen2, d_dens_death2, d_dens_col2, d_dens_bnd2
     2568  REAL, DIMENSION (klon)                                :: d_adens_death2, d_adens_icol2, d_adens_acol2, d_adens_bnd2
     2569
     2570  ! Variables pour les GW
     2571  REAL, DIMENSION (klon)                                :: ll
     2572  REAL, DIMENSION (klon, klev)                          :: n2
     2573  REAL, DIMENSION (klon, klev)                          :: cgw
     2574  REAL, DIMENSION (klon, klev)                          :: tgw
     2575
     2576  ! Variables liees au calcul de hw
     2577  REAL, DIMENSION (klon)                                :: ptop
     2578  REAL, DIMENSION (klon)                                :: sum_dth
     2579  REAL, DIMENSION (klon)                                :: dthmin
     2580  REAL, DIMENSION (klon)                                :: z, dz, hw0
     2581  INTEGER, DIMENSION (klon)                             :: ktop, kupper
     2582
     2583  ! Variables liees au test de la forme triangulaire du profil de Delta_theta
     2584  REAL, DIMENSION (klon)                                :: sum_half_dth
     2585  REAL, DIMENSION (klon)                                :: dz_half
     2586
     2587  ! Sub-timestep tendencies and related variables
     2588  REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
     2589  REAL, DIMENSION (klon, klev)                          :: d_deltat_dadv, d_deltaq_dadv
     2590  REAL, DIMENSION (klon, klev)                          :: d_deltat_lsadv, d_deltaq_lsadv
     2591  REAL, DIMENSION (klon, klev)                          :: d_deltat_entrp
     2592  REAL, DIMENSION (klon, klev)                          :: d_deltat_spread, d_deltaq_spread
     2593
     2594  REAL, DIMENSION (klon, klev)                          :: d_tb, d_qb
     2595  REAL, DIMENSION (klon, klev)                          :: d_tb_dadv, d_qb_dadv
     2596  REAL, DIMENSION (klon, klev)                          :: d_tb_spread, d_qb_spread
     2597
     2598  REAL, DIMENSION (klon)                                :: d_wdens, d_awdens, d_sigmaw, d_asigmaw
     2599  REAL, DIMENSION (klon)                                :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
     2600  REAL, DIMENSION (klon)                                :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd
     2601  REAL, DIMENSION (klon)                                :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
     2602  REAL, DIMENSION (klon)                                :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd
     2603  REAL, DIMENSION (klon)                                :: agfl              !! gust front length of active wakes
     2604                                                                             !!  per unit area
     2605  REAL, DIMENSION (klon)                                :: alpha, alpha_tot
     2606  REAL, DIMENSION (klon)                                :: q0_min, q1_min
     2607  LOGICAL, DIMENSION (klon)                             :: wk_adv, ok_qx_qw
     2608
     2609
     2610  ! Autres variables internes
     2611  INTEGER                                               ::isubstep, k, i, igout
     2612
     2613  REAL                                                  :: wdensmin
     2614
     2615  REAL                                                  :: sigmaw_targ
     2616  REAL                                                  :: wdens_targ
     2617  REAL                                                  :: d_sigmaw_targ
     2618  REAL                                                  :: d_wdens_targ
     2619
     2620  REAL, DIMENSION (klon)                                :: dsigspread  !rate of change of sigmaw due to spreading
     2621
     2622  REAL, DIMENSION (klon)                                :: sum_thx, sum_tx, sum_qx, sum_thvx
     2623  REAL, DIMENSION (klon)                                :: sum_dq
     2624  REAL, DIMENSION (klon)                                :: sum_dtdwn, sum_dqdwn
     2625  REAL, DIMENSION (klon)                                :: av_thx, av_tx, av_qx, av_thvx
     2626  REAL, DIMENSION (klon)                                :: av_dth, av_dq
     2627  REAL, DIMENSION (klon)                                :: av_dtdwn, av_dqdwn
     2628
     2629  REAL, DIMENSION (klon, klev)                          :: rho
     2630  REAL, DIMENSION (klon, klev+1)                        :: rhoh
     2631  REAL, DIMENSION (klon, klev)                          :: zh
     2632  REAL, DIMENSION (klon, klev+1)                        :: zhh
     2633
     2634  REAL, DIMENSION (klon, klev)                          :: thb, thx
     2635
     2636  REAL, DIMENSION (klon, klev)                          :: omgbw
     2637  REAL, DIMENSION (klon)                                :: pupper
     2638  REAL, DIMENSION (klon)                                :: omgtop
     2639  REAL, DIMENSION (klon, klev)                          :: dp_omgbw
     2640  REAL, DIMENSION (klon)                                :: ztop, dztop
     2641  REAL, DIMENSION (klon, klev)                          :: alpha_up
     2642
     2643  REAL, DIMENSION (klon)                                :: rre1, rre2
     2644  REAL                                                  :: rrd1, rrd2
     2645  REAL, DIMENSION (klon, klev)                          :: th1, th2, q1, q2
     2646  REAL, DIMENSION (klon, klev)                          :: d_th1, d_th2, d_dth
     2647  REAL, DIMENSION (klon, klev)                          :: d_q1, d_q2, d_dq
     2648  REAL, DIMENSION (klon, klev)                          :: omgbdq
     2649
     2650  REAL, DIMENSION (klon)                                :: wape2, cstar2, heff
     2651                                                       
     2652  REAL, DIMENSION (klon, klev)                          :: crep
     2653                                                       
     2654  REAL, DIMENSION (klon, klev)                          :: ppi
     2655
     2656  ! cc nrlmd
     2657  REAL, DIMENSION (klon)                                :: death_rate
     2658!!  REAL, DIMENSION (klon)                                :: nat_rate
     2659  REAL, DIMENSION (klon, klev)                          :: entr   ! total entrainment into wakes (spread + birth)
     2660  REAL, DIMENSION (klon, klev)                          :: entr_p ! entrainment into wakes (due to births)
     2661  REAL, DIMENSION (klon, klev)                          :: detr   ! detrainment from wakes (due to deaths)
     2662
     2663  REAL, DIMENSION(klon)                                 :: sigmaw_in, asigmaw_in ! pour les prints
     2664  REAL, DIMENSION(klon)                                 :: wdens_in, awdens_in   ! pour les prints
     2665
     2666  !!-- variables liees au nouveau calcul de ptop et hw
     2667  REAL, DIMENSION (klon, klev)                          :: int_dth
     2668  REAL, DIMENSION (klon, klev)                          :: zzz, dzzz
     2669  REAL                                                  :: epsil
     2670  REAL, DIMENSION (klon)                                :: ptop1
     2671  INTEGER, DIMENSION (klon)                             :: ktop1
     2672  REAL, DIMENSION (klon)                                :: omega
     2673  REAL, DIMENSION (klon)                                :: h_zzz
     2674
     2675  !! Bidouilles
     2676  REAL                                                  :: iwkadv
     2677  REAL                                                  :: iokqxqw
     2678
     2679!print*,'WAKE LJYFz'
     2680
     2681  ! -------------------------------------------------------------------------
     2682  ! Initialisations
     2683  ! -------------------------------------------------------------------------
     2684  ! ALON = 3.e5
     2685  ! alon = 1.E6
     2686
     2687  !  Provisionnal; to be suppressed when f_shear is parameterized
     2688  f_shear(:) = 1.       ! 0. for strong shear, 1. for weak shear
     2689
     2690  ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
     2691
     2692  ! coefgw : Coefficient pour les ondes de gravite
     2693  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
     2694  ! wdens : Densite surfacique de poche froide
     2695  ! -------------------------------------------------------------------------
     2696
     2697  ! cc nrlmd      coefgw=10
     2698  ! coefgw=1
     2699  ! wdens0 = 1.0/(alon**2)
     2700  ! cc nrlmd      wdens = 1.0/(alon**2)
     2701  ! cc nrlmd      stark = 0.50
     2702  ! CRtest
     2703  ! cc nrlmd      alpk=0.1
     2704  ! alpk = 1.0
     2705  ! alpk = 0.5
     2706  ! alpk = 0.05
     2707!
     2708 igout = klon/2+1/klon 
     2709!
     2710!   sub-time-stepping parameters
     2711  dtimesub = dtime/wk_nsub
     2712!
     2713 IF (iflag_wk_pop_dyn == 0) THEN
     2714  ! Initialisation de toutes des densites a wdens_ref.
     2715  ! Les densites peuvent evoluer si les poches debordent
     2716  ! (voir au tout debut de la boucle sur les substeps)
     2717  !jyg<
     2718  !!  wdens(:) = wdens_ref
     2719   DO i = 1,klon
     2720     wdens(i) = wdens_ref(znatsurf(i)+1)
     2721   ENDDO
     2722  !>jyg
     2723 ENDIF  ! (iflag_wk_pop_dyn == 0)
     2724!
     2725 IF (iflag_wk_pop_dyn >=1) THEN
     2726   IF (iflag_wk_pop_dyn == 3) THEN
     2727     wdensmin = wdensthreshold
     2728   ELSE
     2729     wdensmin = wdensinit
     2730   ENDIF
     2731 ENDIF
     2732
     2733  ! print*,'stark',stark
     2734  ! print*,'alpk',alpk
     2735  ! print*,'wdens',wdens
     2736  ! print*,'coefgw',coefgw
     2737  ! cc
     2738  ! Minimum value for |T_wake - T_undist|. Used for wake top definition
     2739  ! -------------------------------------------------------------------------
     2740
     2741   delta_t_min = 0.2
     2742
     2743  ! 1. - Save initial values, initialize tendencies, initialize output fields
     2744  ! ------------------------------------------------------------------------
     2745
     2746!jyg<
     2747!!  DO k = 1, klev
     2748!!    DO i = 1, klon
     2749!!      ppi(i, k) = pi(i, k)
     2750!!      deltatw0(i, k) = deltatw(i, k)
     2751!!      deltaqw0(i, k) = deltaqw(i, k)
     2752!!      tb(i, k) = tb0(i, k)
     2753!!      qb(i, k) = qb0(i, k)
     2754!!      dtls(i, k) = 0.
     2755!!      dqls(i, k) = 0.
     2756!!      d_deltat_gw(i, k) = 0.
     2757!!      d_tb(i, k) = 0.
     2758!!      d_qb(i, k) = 0.
     2759!!      d_deltatw(i, k) = 0.
     2760!!      d_deltaqw(i, k) = 0.
     2761!!      ! IM 060508 beg
     2762!!      d_deltatw2(i, k) = 0.
     2763!!      d_deltaqw2(i, k) = 0.
     2764!!      ! IM 060508 end
     2765!!    END DO
     2766!!  END DO
     2767      ppi(:,:) = pi(:,:)
     2768      deltatw0(:,:) = deltatw(:,:)
     2769      deltaqw0(:,:) = deltaqw(:,:)
     2770      tb(:,:) = tb0(:,:)
     2771      qb(:,:) = qb0(:,:)
     2772      dtls(:,:) = 0.
     2773      dqls(:,:) = 0.
     2774      d_deltat_gw(:,:) = 0.
    4602775
    4612776      detr(:,:) = 0.
    4622777      entr(:,:) = 0.
    4632778      entr_p(:,:) = 0.
    464       c5p(:,:) = 0.
    4652779
    4662780      th1(:,:) = 0.
    … …  
    15503864
    15513865          dth(i, k) = deltatw(i, k)/ppi(i, k)
    1552           c5p(i,k) = - entr_p(i,k)*deltatw(i,k)/sigmaw(i)
    15533866
    15543867          !! Entrainment due to spread is supposed to be included in the differential advection
    … …  
    15703883!!    IF (prt_level>=10) THEN
    15713884    IF (prt_level>=10 .and. wk_adv(igout)) THEN
    1572       PRINT *, 'wake-4.4, isubstep= ', isubstep,' c5p(igout,k) ', (k,c5p(igout,k), k=1,klev)
    15733885      PRINT *, 'wake-4.4, isubstep= ', isubstep,' deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev)
    15743886      PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_deltat_dcv(igout,k) ', (k,d_deltat_dcv(igout,k), k=1,klev)
    … …  
    24524764
    24534765 RETURN
    2454 END SUBROUTINE wake
     4766END SUBROUTINE wake2
    24554767
    24564768SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon_loc, qb, d_qb, deltaqw, &

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

    r5709 r5763  
    25272527          ENDIF
    25282528          CALL histwrite_phy(o_tntc, zx_tmp_fi3d)
    2529        ELSE IF(iflag_thermals.GE.1.AND.iflag_wake.EQ.1) THEN
     2529       ELSE IF(iflag_thermals.GE.1.AND.iflag_wake.GE.1) THEN
    25302530          IF (vars_defined) THEN
    25312531             zx_tmp_fi3d(1:klon,1:klev)=d_t_con(1:klon,1:klev)/pdtphys + &
    … …  
    25472547          IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_con(1:klon,1:klev)/pdtphys
    25482548          CALL histwrite_phy(o_tnhusc, zx_tmp_fi3d)
    2549        ELSE IF (iflag_thermals.GE.1.AND.iflag_wake.EQ.1) THEN
     2549       ELSE IF (iflag_thermals.GE.1.AND.iflag_wake.GE.1) THEN
    25502550          IF (vars_defined) THEN
    25512551             zx_tmp_fi3d(1:klon,1:klev)=d_q_con(1:klon,1:klev)/pdtphys + &

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

    r5715 r5763  
    35253525          ENDDO
    35263526
    3527           IF (iflag_wake==2) THEN
     3527          IF (mod(iflag_wake,10)==2) THEN
    35283528             ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
    35293529             DO k = 1,klev
    … …  
    35333533                     ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/phys_tstep
    35343534             ENDDO
    3535           ELSEIF (iflag_wake==3) THEN
     3535          ELSEIF (mod(iflag_wake,10)==3) THEN
    35363536             ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
    35373537             DO k = 1,klev
Note: See TracChangeset for help on using the changeset viewer.