Ignore:
Timestamp:
Jul 28, 2025, 7:23:15 PM (8 days ago)
Author:
aborella
Message:

Merge with trunk r5789

Location:
LMDZ6/branches/contrails
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/contrails

  • LMDZ6/branches/contrails/libf/phylmd/lmdz_wake.f90

    r5618 r5791  
    99  !$OMP THREADPRIVATE(first_call)
    1010
    11   PUBLIC wake, wake_first
     11  PUBLIC wake, wake2, wake_first
    1212
    1313CONTAINS
    … …  
    23532353END SUBROUTINE wake
    23542354
     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  ! Tgen   : 1 sur le temps caracteristique d'amortissement par les naissances de poches
     2481
     2482  ! -------------------------------------------------------------------------
     2483  ! Declaration de variables
     2484  ! -------------------------------------------------------------------------
     2485
     2486
     2487  ! Arguments en entree
     2488  ! --------------------
     2489
     2490  INTEGER,                          INTENT(IN)          :: klon,klev
     2491  INTEGER, DIMENSION (klon),        INTENT(IN)          :: znatsurf
     2492  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: p, pi
     2493  REAL, DIMENSION (klon, klev+1),   INTENT(IN)          :: ph
     2494  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: omgb
     2495  REAL,                             INTENT(IN)          :: dtime
     2496  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: tb0, qb0
     2497  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dtdwn, dqdwn
     2498  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: amdwn, amup
     2499  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dta, dqa
     2500  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen
     2501  REAL, DIMENSION (klon),           INTENT(IN)          :: sigd_con
     2502  REAL, DIMENSION (klon),           INTENT(IN)          :: Cin
     2503
     2504  !
     2505  ! Input/Output
     2506  ! State variables
     2507  REAL, DIMENSION (klon, klev),     INTENT(INOUT)       :: deltatw, deltaqw
     2508  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw
     2509  REAL, DIMENSION (klon),           INTENT(INOUT)       :: asigmaw
     2510  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens
     2511  REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens
     2512
     2513  ! Sorties
     2514  ! --------
     2515
     2516  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dth
     2517  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: tx, qx
     2518  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls
     2519!!  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
     2520  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_dcv, d_deltaq_dcv
     2521  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: wkspread    !  unused (jyg)
     2522  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: omgbdth, omg
     2523  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dp_omgb, dp_deltomg
     2524  REAL, DIMENSION (klon),           INTENT(OUT)         :: hw, wape, fip, gfl, cstar
     2525  INTEGER, DIMENSION (klon),        INTENT(OUT)         :: ktopw
     2526  ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields
     2527  !                                 computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
     2528  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw
     2529  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
     2530  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2
     2531
     2532  ! Variables internes
     2533  ! -------------------
     2534
     2535  ! Variables a fixer
     2536
     2537  REAL                                                  :: delta_t_min
     2538  REAL                                                  :: dtimesub
     2539  REAL                                                  :: wdens0
     2540  ! IM 080208
     2541  LOGICAL, DIMENSION (klon)                             :: gwake
     2542
     2543  ! Variables de sauvegarde
     2544  REAL, DIMENSION (klon, klev)                          :: deltatw0
     2545  REAL, DIMENSION (klon, klev)                          :: deltaqw0
     2546  REAL, DIMENSION (klon, klev)                          :: tb, qb
     2547
     2548  ! Variables liees a la dynamique de population 1
     2549  REAL, DIMENSION(klon)                                 :: act
     2550  REAL, DIMENSION(klon)                                 :: rad_wk, tau_wk_inv
     2551  REAL, DIMENSION(klon)                                 :: f_shear
     2552  REAL, DIMENSION(klon)                                 :: drdt
     2553 
     2554  ! Variables liees a la dynamique de population 2
     2555  REAL, DIMENSION(klon)                                 :: cont_fact 
     2556 
     2557  ! Variables liees a la dynamique de population 3
     2558  REAL, DIMENSION(klon)                                 :: arad_wk, irad_wk
     2559 
     2560!!  REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col
     2561  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
     2562  LOGICAL, DIMENSION (klon)                             :: kill_wake
     2563  REAL                                                  :: drdt_pos
     2564  REAL                                                  :: tau_wk_inv_min
     2565  ! Some components of the tendencies of state variables 
     2566  REAL, DIMENSION (klon)                                :: d_sig_gen2, d_sig_death2, d_sig_col2, d_sig_spread2, d_sig_bnd2
     2567  REAL, DIMENSION (klon)                                :: d_asig_death2, d_asig_aicol2, d_asig_iicol2, d_asig_spread2, d_asig_bnd2
     2568  REAL, DIMENSION (klon)                                :: d_dens_gen2, d_dens_death2, d_dens_col2, d_dens_bnd2
     2569  REAL, DIMENSION (klon)                                :: d_adens_death2, d_adens_icol2, d_adens_acol2, d_adens_bnd2
     2570
     2571  ! Variables pour les GW
     2572  REAL, DIMENSION (klon)                                :: ll
     2573  REAL, DIMENSION (klon, klev)                          :: n2
     2574  REAL, DIMENSION (klon, klev)                          :: cgw
     2575  REAL, DIMENSION (klon, klev)                          :: tgw
     2576  REAL, DIMENSION (klon, klev)                          :: tgen
     2577
     2578  ! Variables liees au calcul de hw
     2579  REAL, DIMENSION (klon)                                :: ptop
     2580  REAL, DIMENSION (klon)                                :: sum_dth
     2581  REAL, DIMENSION (klon)                                :: dthmin
     2582  REAL, DIMENSION (klon)                                :: z, dz, hw0
     2583  INTEGER, DIMENSION (klon)                             :: ktop, kupper
     2584
     2585  ! Variables liees au test de la forme triangulaire du profil de Delta_theta
     2586  REAL, DIMENSION (klon)                                :: sum_half_dth
     2587  REAL, DIMENSION (klon)                                :: dz_half
     2588
     2589  ! Sub-timestep tendencies and related variables
     2590  REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
     2591  REAL, DIMENSION (klon, klev)                          :: d_deltat_dadv, d_deltaq_dadv
     2592  REAL, DIMENSION (klon, klev)                          :: d_deltat_lsadv, d_deltaq_lsadv
     2593  REAL, DIMENSION (klon, klev)                          :: d_deltat_entrp
     2594  REAL, DIMENSION (klon, klev)                          :: d_deltat_spread, d_deltaq_spread
     2595
     2596  REAL, DIMENSION (klon, klev)                          :: d_tb, d_qb
     2597  REAL, DIMENSION (klon, klev)                          :: d_tb_dadv, d_qb_dadv
     2598  REAL, DIMENSION (klon, klev)                          :: d_tb_spread, d_qb_spread
     2599
     2600  REAL, DIMENSION (klon)                                :: d_wdens, d_awdens, d_sigmaw, d_asigmaw
     2601  REAL, DIMENSION (klon)                                :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
     2602  REAL, DIMENSION (klon)                                :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd
     2603  REAL, DIMENSION (klon)                                :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
     2604  REAL, DIMENSION (klon)                                :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd
     2605  REAL, DIMENSION (klon)                                :: agfl              !! gust front length of active wakes
     2606                                                                             !!  per unit area
     2607  REAL, DIMENSION (klon)                                :: alpha, alpha_tot
     2608  REAL, DIMENSION (klon)                                :: q0_min, q1_min
     2609  LOGICAL, DIMENSION (klon)                             :: wk_adv, ok_qx_qw
     2610
     2611
     2612  ! Autres variables internes
     2613  INTEGER                                               ::isubstep, k, i, igout
     2614
     2615  REAL                                                  :: wdensmin
     2616
     2617  REAL                                                  :: sigmaw_targ
     2618  REAL                                                  :: wdens_targ
     2619  REAL                                                  :: d_sigmaw_targ
     2620  REAL                                                  :: d_wdens_targ
     2621
     2622  REAL, DIMENSION (klon)                                :: dsigspread  !rate of change of sigmaw due to spreading
     2623
     2624  REAL, DIMENSION (klon)                                :: sum_thx, sum_tx, sum_qx, sum_thvx
     2625  REAL, DIMENSION (klon)                                :: sum_dq
     2626  REAL, DIMENSION (klon)                                :: sum_dtdwn, sum_dqdwn
     2627  REAL, DIMENSION (klon)                                :: av_thx, av_tx, av_qx, av_thvx
     2628  REAL, DIMENSION (klon)                                :: av_dth, av_dq
     2629  REAL, DIMENSION (klon)                                :: av_dtdwn, av_dqdwn
     2630
     2631  REAL, DIMENSION (klon, klev)                          :: rho
     2632  REAL, DIMENSION (klon, klev+1)                        :: rhoh
     2633  REAL, DIMENSION (klon, klev)                          :: zh
     2634  REAL, DIMENSION (klon, klev+1)                        :: zhh
     2635
     2636  REAL, DIMENSION (klon, klev)                          :: thb, thx
     2637
     2638  REAL, DIMENSION (klon, klev)                          :: omgbw
     2639  REAL, DIMENSION (klon)                                :: pupper
     2640  REAL, DIMENSION (klon)                                :: omgtop
     2641  REAL, DIMENSION (klon, klev)                          :: dp_omgbw
     2642  REAL, DIMENSION (klon)                                :: ztop, dztop
     2643  REAL, DIMENSION (klon, klev)                          :: alpha_up
     2644
     2645  REAL, DIMENSION (klon)                                :: rre1, rre2
     2646  REAL                                                  :: rrd1, rrd2
     2647  REAL, DIMENSION (klon, klev)                          :: th1, th2, q1, q2
     2648  REAL, DIMENSION (klon, klev)                          :: d_th1, d_th2, d_dth
     2649  REAL, DIMENSION (klon, klev)                          :: d_q1, d_q2, d_dq
     2650  REAL, DIMENSION (klon, klev)                          :: omgbdq
     2651
     2652  REAL, DIMENSION (klon)                                :: wape2, cstar2, heff
     2653                                                       
     2654  REAL, DIMENSION (klon, klev)                          :: crep
     2655                                                       
     2656  REAL, DIMENSION (klon, klev)                          :: ppi
     2657
     2658  ! cc nrlmd
     2659  REAL, DIMENSION (klon)                                :: death_rate
     2660!!  REAL, DIMENSION (klon)                                :: nat_rate
     2661  REAL, DIMENSION (klon, klev)                          :: entr   ! total entrainment into wakes (spread + birth)
     2662  REAL, DIMENSION (klon, klev)                          :: entr_p ! entrainment into wakes (due to births)
     2663  REAL, DIMENSION (klon, klev)                          :: detr   ! detrainment from wakes (due to deaths)
     2664
     2665  REAL, DIMENSION(klon)                                 :: sigmaw_in, asigmaw_in ! pour les prints
     2666  REAL, DIMENSION(klon)                                 :: wdens_in, awdens_in   ! pour les prints
     2667
     2668  !!-- variables liees au nouveau calcul de ptop et hw
     2669  REAL, DIMENSION (klon, klev)                          :: int_dth
     2670  REAL, DIMENSION (klon, klev)                          :: zzz, dzzz
     2671  REAL                                                  :: epsil
     2672  REAL, DIMENSION (klon)                                :: ptop1
     2673  INTEGER, DIMENSION (klon)                             :: ktop1
     2674  REAL, DIMENSION (klon)                                :: omega
     2675  REAL, DIMENSION (klon)                                :: h_zzz
     2676
     2677  !! Bidouilles
     2678  REAL                                                  :: iwkadv
     2679  REAL                                                  :: iokqxqw
     2680
     2681!print*,'WAKE LJYFz'
     2682
     2683  ! -------------------------------------------------------------------------
     2684  ! Initialisations
     2685  ! -------------------------------------------------------------------------
     2686  ! ALON = 3.e5
     2687  ! alon = 1.E6
     2688
     2689  !  Provisionnal; to be suppressed when f_shear is parameterized
     2690  f_shear(:) = 1.       ! 0. for strong shear, 1. for weak shear
     2691
     2692  ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
     2693
     2694  ! coefgw : Coefficient pour les ondes de gravite
     2695  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
     2696  ! wdens : Densite surfacique de poche froide
     2697  ! -------------------------------------------------------------------------
     2698
     2699  ! cc nrlmd      coefgw=10
     2700  ! coefgw=1
     2701  ! wdens0 = 1.0/(alon**2)
     2702  ! cc nrlmd      wdens = 1.0/(alon**2)
     2703  ! cc nrlmd      stark = 0.50
     2704  ! CRtest
     2705  ! cc nrlmd      alpk=0.1
     2706  ! alpk = 1.0
     2707  ! alpk = 0.5
     2708  ! alpk = 0.05
     2709!
     2710 igout = klon/2+1/klon 
     2711!
     2712!   sub-time-stepping parameters
     2713  dtimesub = dtime/wk_nsub
     2714!
     2715 IF (iflag_wk_pop_dyn == 0) THEN
     2716  ! Initialisation de toutes des densites a wdens_ref.
     2717  ! Les densites peuvent evoluer si les poches debordent
     2718  ! (voir au tout debut de la boucle sur les substeps)
     2719  !jyg<
     2720  !!  wdens(:) = wdens_ref
     2721   DO i = 1,klon
     2722     wdens(i) = wdens_ref(znatsurf(i)+1)
     2723   ENDDO
     2724  !>jyg
     2725 ENDIF  ! (iflag_wk_pop_dyn == 0)
     2726!
     2727 IF (iflag_wk_pop_dyn >=1) THEN
     2728   IF (iflag_wk_pop_dyn == 3) THEN
     2729     wdensmin = wdensthreshold
     2730   ELSE
     2731     wdensmin = wdensinit
     2732   ENDIF
     2733 ENDIF
     2734
     2735  ! print*,'stark',stark
     2736  ! print*,'alpk',alpk
     2737  ! print*,'wdens',wdens
     2738  ! print*,'coefgw',coefgw
     2739  ! cc
     2740  ! Minimum value for |T_wake - T_undist|. Used for wake top definition
     2741  ! -------------------------------------------------------------------------
     2742
     2743   delta_t_min = 0.2
     2744
     2745  ! 1. - Save initial values, initialize tendencies, initialize output fields
     2746  ! ------------------------------------------------------------------------
     2747
     2748!jyg<
     2749!!  DO k = 1, klev
     2750!!    DO i = 1, klon
     2751!!      ppi(i, k) = pi(i, k)
     2752!!      deltatw0(i, k) = deltatw(i, k)
     2753!!      deltaqw0(i, k) = deltaqw(i, k)
     2754!!      tb(i, k) = tb0(i, k)
     2755!!      qb(i, k) = qb0(i, k)
     2756!!      dtls(i, k) = 0.
     2757!!      dqls(i, k) = 0.
     2758!!      d_deltat_gw(i, k) = 0.
     2759!!      d_tb(i, k) = 0.
     2760!!      d_qb(i, k) = 0.
     2761!!      d_deltatw(i, k) = 0.
     2762!!      d_deltaqw(i, k) = 0.
     2763!!      ! IM 060508 beg
     2764!!      d_deltatw2(i, k) = 0.
     2765!!      d_deltaqw2(i, k) = 0.
     2766!!      ! IM 060508 end
     2767!!    END DO
     2768!!  END DO
     2769      ppi(:,:) = pi(:,:)
     2770      deltatw0(:,:) = deltatw(:,:)
     2771      deltaqw0(:,:) = deltaqw(:,:)
     2772      tb(:,:) = tb0(:,:)
     2773      qb(:,:) = qb0(:,:)
     2774      dtls(:,:) = 0.
     2775      dqls(:,:) = 0.
     2776      d_deltat_gw(:,:) = 0.
     2777
     2778      detr(:,:) = 0.
     2779      entr(:,:) = 0.
     2780      entr_p(:,:) = 0.
     2781
     2782      th1(:,:) = 0.
     2783      th2(:,:) = 0.
     2784      q1(:,:) = 0.
     2785      q2(:,:) = 0.
     2786
     2787      d_tb(:,:) = 0.
     2788      d_tb_dadv(:,:) = 0.
     2789      d_tb_spread(:,:) = 0.
     2790
     2791      d_qb(:,:) = 0.
     2792      d_qb_dadv(:,:) = 0.
     2793      d_qb_spread(:,:) = 0.
     2794
     2795      d_deltatw(:,:) = 0.
     2796      d_deltat_dadv  (:,:) = 0.
     2797      d_deltat_lsadv (:,:) = 0.
     2798      d_deltat_dcv   (:,:) = 0.
     2799      d_deltat_spread(:,:) = 0.
     2800
     2801      d_deltaqw(:,:) = 0.
     2802      d_deltaq_dadv(:,:) = 0.
     2803      d_deltaq_lsadv(:,:) = 0.
     2804      d_deltaq_dcv(:,:) = 0.
     2805      d_deltaq_spread(:,:) = 0.
     2806
     2807      d_deltatw2(:,:) = 0.
     2808      d_deltaqw2(:,:) = 0.
     2809
     2810      d_sig_gen2(:)   = 0.
     2811      d_sig_death2(:) = 0.
     2812      d_sig_col2(:)   = 0.
     2813      d_sig_spread2(:)= 0.
     2814      d_asig_death2(:) = 0.
     2815      d_asig_iicol2(:) = 0.
     2816      d_asig_aicol2(:) = 0.
     2817      d_asig_spread2(:)= 0.
     2818      d_asig_bnd2(:) = 0.
     2819      d_asigmaw2(:) = 0.
     2820!
     2821      d_dens_gen2(:)   = 0.
     2822      d_dens_death2(:) = 0.
     2823      d_dens_col2(:)   = 0.
     2824      d_dens_bnd2(:) = 0.
     2825      d_wdens2(:) = 0.
     2826      d_adens_bnd2(:) = 0.
     2827      d_awdens2(:) = 0.
     2828      d_adens_death2(:) = 0.
     2829      d_adens_icol2(:) = 0.
     2830      d_adens_acol2(:) = 0.
     2831
     2832      IF (iflag_wk_act == 0) THEN
     2833        act(:) = 0.
     2834      ELSEIF (iflag_wk_act == 1) THEN
     2835        act(:) = 1.
     2836      ENDIF
     2837
     2838!!  DO i = 1, klon
     2839!!   sigmaw_in(i) = sigmaw(i)
     2840!!  END DO
     2841   sigmaw_in(:)  = sigmaw(:)
     2842   asigmaw_in(:) = asigmaw(:)
     2843!>jyg
     2844!
     2845  IF (iflag_wk_pop_dyn >= 1) THEN
     2846    awdens_in(:) = awdens(:)
     2847    wdens_in(:) = wdens(:)
     2848!!    wdens(:) = wdens(:) + wgen(:)*dtime
     2849!!    d_wdens2(:) = wgen(:)*dtime
     2850!!  ELSE
     2851  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     2852
     2853
     2854  ! sigmaw1=sigmaw
     2855  ! IF (sigd_con.GT.sigmaw1) THEN
     2856  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
     2857  ! ENDIF
     2858  IF (iflag_wk_pop_dyn >= 1) THEN
     2859    DO i = 1, klon
     2860      d_dens_gen2(i)   = 0.
     2861      d_dens_death2(i) = 0.
     2862      d_dens_col2(i)   = 0.
     2863      d_awdens2(i) = 0.
     2864      IF (wdens(i) < wdensthreshold) THEN
     2865  !!      wdens_targ = max(wdens(i),wdensmin)
     2866        wdens_targ = max(wdens(i),wdensinit)
     2867        d_dens_bnd2(i) = wdens_targ - wdens(i)
     2868        d_wdens2(i) = wdens_targ - wdens(i)
     2869        wdens(i) = wdens_targ
     2870      ELSE
     2871        d_dens_bnd2(i) = 0.
     2872        d_wdens2(i) = 0.
     2873      ENDIF  !! (wdens(i) < wdensthreshold)
     2874    END DO
     2875    IF (iflag_wk_pop_dyn >= 2) THEN
     2876      DO i = 1, klon 
     2877        IF (awdens(i) < wdensthreshold) THEN
     2878!!          wdens_targ = min(max(awdens(i),wdensmin),wdens(i))
     2879            wdens_targ = min(max(awdens(i),wdensinit),wdens(i))
     2880            d_adens_bnd2(i) = wdens_targ - awdens(i)
     2881            d_awdens2(i) = wdens_targ - awdens(i)
     2882            awdens(i) = wdens_targ
     2883        ELSE
     2884            wdens_targ = min(awdens(i), wdens(i))
     2885            d_adens_bnd2(i) = wdens_targ - awdens(i)
     2886            d_awdens2(i) = wdens_targ - awdens(i)
     2887            awdens(i) = wdens_targ
     2888        ENDIF
     2889      END DO
     2890    ENDIF ! (iflag_wk_pop_dyn >= 2)
     2891  ELSE 
     2892    DO i = 1, klon
     2893      d_awdens2(i) = 0.
     2894      d_wdens2(i) = 0.
     2895    END DO
     2896  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     2897!
     2898  DO i = 1, klon
     2899    sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
     2900    d_sig_bnd2(i) = sigmaw_targ - sigmaw(i)
     2901    d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
     2902    sigmaw(i) = sigmaw_targ
     2903  END DO
     2904!
     2905  IF (iflag_wk_pop_dyn == 3) THEN
     2906     DO i = 1, klon 
     2907        IF ((wdens(i)-awdens(i)) <= smallestreal) THEN
     2908          sigmaw_targ = sigmaw(i)
     2909        ELSE
     2910          sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i))
     2911        ENDIF
     2912        d_asig_bnd2(i) = sigmaw_targ - asigmaw(i)
     2913        d_asigmaw2(i) = sigmaw_targ - asigmaw(i)
     2914        asigmaw(i) = sigmaw_targ
     2915     END DO
     2916  ENDIF ! (iflag_wk_pop_dyn == 3)
     2917
     2918  wape(:) = 0.
     2919  wape2(:) = 0.
     2920  d_sigmaw(:) = 0.
     2921  d_asigmaw(:) = 0.
     2922  ktopw(:) = 0
     2923!
     2924!<jyg
     2925dth(:,:) = 0.
     2926tx(:,:) = 0.
     2927qx(:,:) = 0.
     2928d_deltat_dcv(:,:) = 0.
     2929d_deltaq_dcv(:,:) = 0.
     2930wkspread(:,:) = 0.
     2931omgbdth(:,:) = 0.
     2932omg(:,:) = 0.
     2933dp_omgb(:,:) = 0.
     2934dp_deltomg(:,:) = 0.
     2935tgen(:,:) = 0.
     2936hw(:) = 0.
     2937wape(:) = 0.
     2938fip(:) = 0.
     2939gfl(:) = 0.
     2940cstar(:) = 0.
     2941ktopw(:) = 0
     2942!
     2943!  Vertical advection local variables
     2944omgbw(:,:) = 0.
     2945omgtop(:) = 0
     2946dp_omgbw(:,:) = 0.
     2947omgbdq(:,:) = 0.
     2948
     2949!>jyg
     2950!
     2951  IF (prt_level>=10) THEN
     2952    PRINT *, 'wake-1, sigmaw(igout) ', sigmaw(igout)
     2953    PRINT *, 'wake-1, deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev)
     2954    PRINT *, 'wake-1, deltaqw(igout,k) ', (k,deltaqw(igout,k), k=1,klev)
     2955    PRINT *, 'wake-1, dowwdraughts, amdwn(igout,k) ', (k,amdwn(igout,k), k=1,klev)
     2956    PRINT *, 'wake-1, dowwdraughts, dtdwn(igout,k) ', (k,dtdwn(igout,k), k=1,klev)
     2957    PRINT *, 'wake-1, dowwdraughts, dqdwn(igout,k) ', (k,dqdwn(igout,k), k=1,klev)
     2958    PRINT *, 'wake-1, updraughts, amup(igout,k) ', (k,amup(igout,k), k=1,klev)
     2959    PRINT *, 'wake-1, updraughts, dta(igout,k) ', (k,dta(igout,k), k=1,klev)
     2960    PRINT *, 'wake-1, updraughts, dqa(igout,k) ', (k,dqa(igout,k), k=1,klev)
     2961  ENDIF
     2962
     2963  ! 2. - Prognostic part
     2964  ! --------------------
     2965
     2966
     2967  ! 2.1 - Undisturbed area and Wake integrals
     2968  ! ---------------------------------------------------------
     2969
     2970  DO i = 1, klon
     2971    z(i) = 0.
     2972    ktop(i) = 0
     2973    kupper(i) = 0
     2974    sum_thx(i) = 0.
     2975    sum_tx(i) = 0.
     2976    sum_qx(i) = 0.
     2977    sum_thvx(i) = 0.
     2978    sum_dth(i) = 0.
     2979    sum_dq(i) = 0.
     2980    sum_dtdwn(i) = 0.
     2981    sum_dqdwn(i) = 0.
     2982
     2983    av_thx(i) = 0.
     2984    av_tx(i) = 0.
     2985    av_qx(i) = 0.
     2986    av_thvx(i) = 0.
     2987    av_dth(i) = 0.
     2988    av_dq(i) = 0.
     2989    av_dtdwn(i) = 0.
     2990    av_dqdwn(i) = 0.
     2991  END DO
     2992
     2993  ! Distance between wakes
     2994  DO i = 1, klon
     2995    ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
     2996  END DO
     2997  ! Potential temperatures and humidity
     2998  ! ----------------------------------------------------------
     2999  DO k = 1, klev
     3000    DO i = 1, klon
     3001      ! write(*,*)'wake 1',i,k,RD,tb(i,k)
     3002      rho(i, k) = p(i, k)/(RD*tb(i,k))
     3003      ! write(*,*)'wake 2',rho(i,k)
     3004      IF (k==1) THEN
     3005        ! write(*,*)'wake 3',i,k,rd,tb(i,k)
     3006        rhoh(i, k) = ph(i, k)/(RD*tb(i,k))
     3007        ! write(*,*)'wake 4',i,k,rd,tb(i,k)
     3008        zhh(i, k) = 0
     3009      ELSE
     3010        ! write(*,*)'wake 5',rd,(tb(i,k)+tb(i,k-1))
     3011        rhoh(i, k) = ph(i, k)*2./(RD*(tb(i,k)+tb(i,k-1)))
     3012        ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
     3013        zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
     3014      END IF
     3015      ! write(*,*)'wake 7',ppi(i,k)
     3016      thb(i, k) = tb(i, k)/ppi(i, k)
     3017      thx(i, k) = (tb(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
     3018      tx(i, k) = tb(i, k) - deltatw(i, k)*sigmaw(i)
     3019      qx(i, k) = qb(i, k) - deltaqw(i, k)*sigmaw(i)
     3020      ! write(*,*)'wake 8',(RD*(tb(i,k)+deltatw(i,k)))
     3021      dth(i, k) = deltatw(i, k)/ppi(i, k)
     3022    END DO
     3023  END DO
     3024
     3025  DO k = 1, klev - 1
     3026    DO i = 1, klon
     3027      IF (k==1) THEN
     3028        n2(i, k) = 0
     3029      ELSE
     3030        n2(i, k) = amax1(0., -RG**2/thb(i,k)*rho(i,k)*(thb(i,k+1)-thb(i,k-1))/ &
     3031                             (p(i,k+1)-p(i,k-1)))
     3032      END IF
     3033      zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
     3034
     3035      cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
     3036      tgw(i, k) = coefgw*cgw(i, k)/ll(i)
     3037    END DO
     3038  END DO
     3039
     3040  DO i = 1, klon
     3041    n2(i, klev) = 0
     3042    zh(i, klev) = 0
     3043    cgw(i, klev) = 0
     3044    tgw(i, klev) = 0
     3045  END DO
     3046
     3047 
     3048  ! Choose an integration bound well above wake top
     3049  ! -----------------------------------------------------------------
     3050
     3051  ! Determine Wake top pressure (Ptop) from buoyancy integral
     3052  ! --------------------------------------------------------
     3053
     3054   Do i=1, klon
     3055       wk_adv(i) = .True.
     3056   Enddo
     3057   Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
     3058                    dth, hw0, rho, delta_t_min, &
     3059                    ktop, wk_adv, h_zzz, ptop1, ktop1)
     3060 
     3061   !!print'("pkupper APPEL ",7i6)',0,int(ptop/100.),int(ptop1/100.),int(pupper/100.),ktop,ktop1,kupper
     3062   
     3063   IF (prt_level>=10) THEN
     3064     PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
     3065  ENDIF
     3066
     3067  ! -5/ Set deltatw & deltaqw to 0 above kupper
     3068
     3069  DO k = 1, klev
     3070    DO i = 1, klon
     3071      IF (k>=kupper(i)) THEN
     3072        deltatw(i, k) = 0.
     3073        deltaqw(i, k) = 0.
     3074        d_deltatw2(i,k) = -deltatw0(i,k)
     3075        d_deltaqw2(i,k) = -deltaqw0(i,k)
     3076      END IF
     3077    END DO
     3078  END DO
     3079
     3080
     3081  ! Vertical gradient of LS omega
     3082
     3083  DO k = 1, klev
     3084    DO i = 1, klon
     3085      IF (k<=kupper(i)) THEN
     3086        dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
     3087      END IF
     3088    END DO
     3089  END DO
     3090
     3091  ! Integrals (and wake top level number)
     3092  ! --------------------------------------
     3093
     3094  ! Initialize sum_thvx to 1st level virt. pot. temp.
     3095
     3096  DO i = 1, klon
     3097    z(i) = 1.
     3098    dz(i) = 1.
     3099    sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i)
     3100    sum_dth(i) = 0.
     3101  END DO
     3102
     3103  DO k = 1, klev
     3104    DO i = 1, klon
     3105      dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
     3106      IF (dz(i)>0) THEN
     3107              ! LJYF : ecriture pas sympa avec un tableau z(i) qui n'est pas utilise come tableau
     3108        z(i) = z(i) + dz(i)
     3109        sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i)
     3110        sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i)
     3111        sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i)
     3112        sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i)
     3113        sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
     3114        sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
     3115        sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
     3116        sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
     3117      END IF
     3118    END DO
     3119  END DO
     3120
     3121  DO i = 1, klon
     3122    hw0(i) = z(i)
     3123  END DO
     3124
     3125
     3126  ! 2.1 - WAPE and mean forcing computation
     3127  ! ---------------------------------------
     3128
     3129  ! ---------------------------------------
     3130
     3131  ! Means
     3132
     3133  DO i = 1, klon
     3134    av_thx(i) = sum_thx(i)/hw0(i)
     3135    av_tx(i) = sum_tx(i)/hw0(i)
     3136    av_qx(i) = sum_qx(i)/hw0(i)
     3137    av_thvx(i) = sum_thvx(i)/hw0(i)
     3138    ! av_thve = sum_thve/hw0
     3139    av_dth(i) = sum_dth(i)/hw0(i)
     3140    av_dq(i) = sum_dq(i)/hw0(i)
     3141    av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     3142    av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     3143
     3144    wape(i) = -RG*hw0(i)*(av_dth(i)+ &
     3145        epsim1*(av_thx(i)*av_dq(i)+av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i)
     3146
     3147  END DO
     3148IF (CPPKEY_IOPHYS_WK) THEN
     3149  IF (.not.phys_sub) CALL iophys_ecrit('wape_a',1,'wape_a','J/kg',wape)
     3150END IF
     3151
     3152  ! 2.2 Prognostic variable update
     3153  ! ------------------------------
     3154
     3155  ! Filter out bad wakes
     3156
     3157  DO k = 1, klev
     3158    DO i = 1, klon
     3159      IF (wape(i)<0.) THEN
     3160        deltatw(i, k) = 0.
     3161        deltaqw(i, k) = 0.
     3162        dth(i, k) = 0.
     3163        d_deltatw2(i,k) = -deltatw0(i,k)
     3164        d_deltaqw2(i,k) = -deltaqw0(i,k)
     3165      END IF
     3166    END DO
     3167  END DO
     3168
     3169  DO i = 1, klon
     3170    IF (wape(i)<0.) THEN
     3171!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
     3172      sigmaw_targ = max(sigmad, sigd_con(i))
     3173      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     3174      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     3175      sigmaw(i) = sigmaw_targ
     3176    ENDIF  !!  (wape(i)<0.)
     3177  ENDDO
     3178  !
     3179  IF (iflag_wk_pop_dyn == 3) THEN
     3180    DO i = 1, klon
     3181      IF (wape(i)<0.) THEN
     3182        sigmaw_targ = max(sigmad, sigd_con(i))
     3183        d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     3184        d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     3185        asigmaw(i) = sigmaw_targ
     3186      ENDIF  !!  (wape(i)<0.)
     3187    ENDDO
     3188  ENDIF  !!  (iflag_wk_pop_dyn == 3)
     3189
     3190  DO i = 1, klon
     3191    IF (wape(i)<0.) THEN
     3192      wape(i) = 0.
     3193      cstar(i) = 0.
     3194      hw(i) = hwmin
     3195      fip(i) = 0.
     3196      gwake(i) = .FALSE.
     3197    ELSE
     3198      hw(i) = hw0(i)
     3199      cstar(i) = stark*sqrt(2.*wape(i))
     3200      gwake(i) = .TRUE.
     3201    END IF
     3202  END DO
     3203  !
     3204
     3205  ! Check qx and qw positivity
     3206  ! --------------------------
     3207  DO i = 1, klon
     3208    q0_min(i) = min((qb(i,1)-sigmaw(i)*deltaqw(i,1)),  &
     3209                    (qb(i,1)+(1.-sigmaw(i))*deltaqw(i,1)))
     3210  END DO
     3211  DO k = 2, klev
     3212    DO i = 1, klon
     3213      q1_min(i) = min((qb(i,k)-sigmaw(i)*deltaqw(i,k)), &
     3214                      (qb(i,k)+(1.-sigmaw(i))*deltaqw(i,k)))
     3215      IF (q1_min(i)<=q0_min(i)) THEN
     3216        q0_min(i) = q1_min(i)
     3217      END IF
     3218    END DO
     3219  END DO
     3220
     3221  DO i = 1, klon
     3222    ok_qx_qw(i) = q0_min(i) >= 0.
     3223    alpha(i) = 1.
     3224    alpha_tot(i) = 1.
     3225  END DO
     3226
     3227  IF (prt_level>=10) THEN
     3228    PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ', &
     3229                      sigmaw(igout), cstar(igout), wape(igout), ktop(igout)
     3230  ENDIF
     3231
     3232
     3233  ! C -----------------------------------------------------------------
     3234  ! Sub-time-stepping
     3235  ! -----------------
     3236
     3237!    wk_nsub and dtimesub definitions moved to begining of routine.
     3238!!  wk_nsub = 10
     3239!!  dtimesub = dtime/wk_nsub
     3240
     3241 
     3242  ! ------------------------------------------------------------------------
     3243  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     3244  ! ------------------------------------------------------------------------
     3245  !
     3246  DO isubstep = 1, wk_nsub
     3247  !
     3248  ! ------------------------------------------------------------------------
     3249  !
     3250    ! wk_adv is the logical flag enabling wake evolution in the time advance
     3251    ! loop
     3252    DO i = 1, klon
     3253      wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
     3254    END DO
     3255IF (CPPKEY_IOPHYS_WK) THEN
     3256    IF (phys_sub) THEN
     3257     iwkadv=0.
     3258     IF (wk_adv(1)) iwkadv=1.
     3259     iokqxqw=0.
     3260     IF (ok_qx_qw(1)) iokqxqw=1.
     3261     CALL iophys_ecrit('iwkadv',1,'iwkadv','',iwkadv)
     3262     CALL iophys_ecrit('iokqxqw',1,'iokqxqw','',iokqxqw)
     3263     CALL iophys_ecrit('alpha',1,'alpha','',alpha(1))
     3264    ENDIF
     3265END IF
     3266    IF (prt_level>=10) THEN
     3267      PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', &
     3268                          isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout)
     3269     
     3270    ENDIF
     3271
     3272    ! cc nrlmd   Ajout d'un recalcul de wdens dans le cas d'un entrainement
     3273    ! negatif de ktop a kupper --------
     3274    ! cc           On calcule pour cela une densite wdens0 pour laquelle on
     3275    ! aurait un entrainement nul ---
     3276    !jyg<
     3277    ! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou
     3278    ! les poches sont insuffisantes pour accueillir tout le flux de masse
     3279    ! des descentes unsaturees. Nous faisons alors l'hypothese que la
     3280    ! convection profonde cree directement de nouvelles poches, sans passer
     3281    ! par les thermiques. La nouvelle valeur de wdens est alors imposee.
     3282
     3283    DO i = 1, klon
     3284      ! c       print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
     3285      ! c     $           isubstep,wk_adv(i),cstar(i),wape(i)
     3286      IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
     3287        IF ( iflag_wk_profile == 0 ) THEN
     3288           omg(i, kupper(i)+1)=-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
     3289                               RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
     3290        ELSE
     3291           omg(i, kupper(i)+1)=0.
     3292        ENDIF
     3293        wdens0 = (sigmaw(i)/(4.*3.14))* &
     3294          ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2)
     3295        IF (prt_level >= 10) THEN
     3296             print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', &
     3297                     omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)
     3298        ENDIF
     3299        IF (wdens(i)<=wdens0*1.1) THEN
     3300          IF (iflag_wk_pop_dyn >= 1) THEN
     3301             d_dens_bnd2(i) = d_dens_bnd2(i) + wdens0 - wdens(i)
     3302             d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
     3303          ENDIF
     3304          wdens(i) = wdens0
     3305        END IF
     3306      END IF
     3307    END DO
     3308
     3309    IF (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl) THEN
     3310!!--------------------------------------------------------
     3311!!Bug : computing gfl and rad_wk before changing sigmaw
     3312!!      This bug exists only for iflag_wk_pop_dyn=0. Otherwise, gfl and rad_wk
     3313!!      are computed within  wake_popdyn
     3314!!--------------------------------------------------------
     3315      DO i = 1, klon
     3316        IF (wk_adv(i)) THEN
     3317          gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     3318          rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
     3319        END IF
     3320      END DO
     3321    ENDIF   ! (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl)
     3322!!--------------------------------------------------------
     3323
     3324    DO i = 1, klon
     3325      IF (wk_adv(i)) THEN
     3326        sigmaw_targ = min(sigmaw(i), sigmaw_max)
     3327        d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     3328        d_sigmaw2(i)  = d_sigmaw2(i)  + sigmaw_targ - sigmaw(i)
     3329        sigmaw(i) = sigmaw_targ
     3330      END IF
     3331    END DO
     3332
     3333    IF (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl) THEN
     3334!!--------------------------------------------------------
     3335!!Fix : computing gfl and rad_wk after changing sigmaw
     3336!!--------------------------------------------------------
     3337      DO i = 1, klon
     3338        IF (wk_adv(i)) THEN
     3339          gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     3340          rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
     3341        END IF
     3342      END DO
     3343    ENDIF   ! (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl)
     3344!!--------------------------------------------------------
     3345
     3346    IF (iflag_wk_pop_dyn >= 1) THEN
     3347  !  The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
     3348  !  Here, it has to be set to zero.
     3349      death_rate(:) = 0.
     3350    ENDIF
     3351 
     3352    IF (iflag_wk_pop_dyn >= 3) THEN
     3353      DO i = 1, klon
     3354        IF (wk_adv(i)) THEN
     3355         sigmaw_targ = min(asigmaw(i), sigmaw_max)
     3356         d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     3357         d_asigmaw2(i)  = d_asigmaw2(i)  + sigmaw_targ - asigmaw(i)
     3358         asigmaw(i) = sigmaw_targ
     3359        ENDIF
     3360      ENDDO
     3361    ENDIF
     3362 
     3363!!--------------------------------------------------------
     3364!!--------------------------------------------------------
     3365    IF (iflag_wk_pop_dyn == 1) THEN
     3366  !
     3367     CALL wake_popdyn_1 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
     3368                  wdensmin, &
     3369                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
     3370                  d_awdens, d_wdens, d_sigmaw, &
     3371                  iflag_wk_act, wk_adv, cin, wape, &
     3372                  drdt, &
     3373                  d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     3374                  d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     3375                  d_wdens_targ, d_sigmaw_targ)
     3376                     
     3377   
     3378!!--------------------------------------------------------
     3379    ELSEIF (iflag_wk_pop_dyn == 2) THEN
     3380  !
     3381     CALL wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, &
     3382                             wdensmin, &
     3383                             sigmaw, wdens, awdens, &   !! state variables
     3384                             gfl, cstar, cin, wape, rad_wk, &
     3385                             d_sigmaw, d_wdens, d_awdens, &  !! tendencies                             
     3386                             cont_fact, &
     3387                             d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     3388                             d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     3389                             d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
     3390sigmaw=sigmaw-d_sigmaw
     3391wdens=wdens-d_wdens
     3392awdens=awdens-d_awdens
     3393
     3394!!--------------------------------------------------------
     3395    ELSEIF (iflag_wk_pop_dyn == 3) THEN
     3396IF (CPPKEY_IOPHYS_WK) THEN
     3397    IF (phys_sub) THEN
     3398     CALL iophys_ecrit('ptop',1,'ptop','Pa',ptop)
     3399     CALL iophys_ecrit('wape',1,'wape','J/kg',wape)
     3400     CALL iophys_ecrit('wgen',1,'wgen','1/(s.m2)',wgen)
     3401     CALL iophys_ecrit('sigmaw',1,'sigmaw','',sigmaw)
     3402     CALL iophys_ecrit('asigmaw',1,'asigmaw','',asigmaw)
     3403     CALL iophys_ecrit('wdens',1,'wdens','1/m2',wdens)
     3404     CALL iophys_ecrit('awdens',1,'awdens','1/m2',awdens)
     3405    ENDIF
     3406END IF
     3407  !
     3408     CALL wake_popdyn_3 ( klon, klev, phys_sub, wk_adv, dtimesub, wgen, &
     3409                             wdensmin, &
     3410                             sigmaw, asigmaw, wdens, awdens, &   !! state variables
     3411                             gfl, agfl, cstar, cin, wape, &
     3412                             rad_wk, arad_wk, irad_wk, &
     3413                             d_sigmaw, d_asigmaw, d_wdens, d_awdens, &  !! tendencies                             
     3414                             d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
     3415                             d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd, &
     3416                             d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
     3417                             d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
     3418IF (CPPKEY_IOPHYS_WK) THEN
     3419    IF (phys_sub) THEN
     3420     CALL iophys_ecrit('rad_wk',1,'rad_wk','m',rad_wk)
     3421     CALL iophys_ecrit('arad_wk',1,'arad_wk','m',arad_wk)
     3422     CALL iophys_ecrit('irad_wk',1,'irad_wk','m',irad_wk)
     3423    ENDIF
     3424END IF
     3425sigmaw=sigmaw-d_sigmaw
     3426asigmaw=asigmaw-d_asigmaw
     3427wdens=wdens-d_wdens
     3428awdens=awdens-d_awdens
     3429   
     3430!!--------------------------------------------------------
     3431    ELSEIF (iflag_wk_pop_dyn == 0) THEN
     3432   
     3433    ! cc nrlmd
     3434
     3435      DO i = 1, klon
     3436        IF (wk_adv(i)) THEN
     3437
     3438          ! cc nrlmd          Introduction du taux de mortalite des poches et
     3439          ! test sur sigmaw_max=0.4
     3440          ! cc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
     3441          IF (sigmaw(i)>=sigmaw_max) THEN
     3442            death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
     3443          ELSE
     3444            death_rate(i) = 0.
     3445          END IF
     3446   
     3447          d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* &
     3448            dtimesub
     3449          ! $              - nat_rate(i)*sigmaw(i)*dtimesub
     3450          ! c        print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
     3451          ! c     $  death_rate(i),ktop(i),kupper(i)',
     3452          ! c     $              d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
     3453          ! c     $  death_rate(i),ktop(i),kupper(i)
     3454   
     3455          ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
     3456          ! sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
     3457          ! wdens = wdens0/(10.*sigmaw)
     3458          ! sigmaw =max(sigmaw,sigd_con)
     3459          ! sigmaw =max(sigmaw,sigmad)
     3460        END IF
     3461      END DO
     3462
     3463    ENDIF   !  (iflag_wk_pop_dyn == 1) ... ELSEIF (iflag_wk_pop_dyn == 0)
     3464!!--------------------------------------------------------
     3465!!--------------------------------------------------------
     3466
     3467IF (CPPKEY_IOPHYS_WK) THEN
     3468    IF (phys_sub) THEN
     3469     CALL iophys_ecrit('wdensa',1,'wdensa','m',wdens)
     3470     CALL iophys_ecrit('awdensa',1,'awdensa','m',awdens)
     3471     CALL iophys_ecrit('sigmawa',1,'sigmawa','m',sigmaw)
     3472     CALL iophys_ecrit('asigmawa',1,'asigmawa','m',asigmaw)
     3473    ENDIF
     3474END IF
     3475    ! calcul de la difference de vitesse verticale poche - zone non perturbee
     3476    ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
     3477    ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
     3478    ! IM 060208 au niveau k=1...
     3479    !JYG 161013 Correction : maintenant omg est dimensionne a klev.
     3480    DO k = 1, klev
     3481      DO i = 1, klon
     3482        IF (wk_adv(i)) THEN !!! nrlmd
     3483          dp_deltomg(i, k) = 0.
     3484        END IF
     3485      END DO
     3486    END DO
     3487    DO k = 1, klev
     3488      DO i = 1, klon
     3489        IF (wk_adv(i)) THEN !!! nrlmd
     3490          omg(i, k) = 0.
     3491        END IF
     3492      END DO
     3493    END DO
     3494
     3495    DO i = 1, klon
     3496      IF (wk_adv(i)) THEN
     3497        z(i) = 0.
     3498        omg(i, 1) = 0.
     3499        dp_deltomg(i, 1) = -(gfl(i)*cstar(i))/(sigmaw(i)*(1-sigmaw(i)))
     3500      END IF
     3501    END DO
     3502
     3503    DO k = 2, klev
     3504      DO i = 1, klon
     3505        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
     3506          dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*RG)
     3507          z(i) = z(i) + dz(i)
     3508          dp_deltomg(i, k) = dp_deltomg(i, 1)
     3509          omg(i, k) = dp_deltomg(i, 1)*z(i)
     3510        END IF
     3511      END DO
     3512    END DO
     3513
     3514    DO i = 1, klon
     3515      IF (wk_adv(i)) THEN
     3516        dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*RG)
     3517        ztop(i) = z(i) + dztop(i)
     3518        omgtop(i) = dp_deltomg(i, 1)*ztop(i)
     3519      END IF
     3520    END DO
     3521
     3522    IF (prt_level>=10) THEN
     3523      PRINT *, 'wake-4.2, omg(igout,k) ', (k,omg(igout,k), k=1,klev)
     3524      PRINT *, 'wake-4.2, omgtop(igout), ptop(igout), ktop(igout) ', &
     3525                          omgtop(igout), ptop(igout), ktop(igout)
     3526    ENDIF
     3527
     3528    ! -----------------
     3529    ! From m/s to Pa/s
     3530    ! -----------------
     3531
     3532    DO i = 1, klon
     3533      IF (wk_adv(i)) THEN
     3534        omgtop(i) = -rho(i, ktop(i))*RG*omgtop(i)
     3535!! LJYF        dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
     3536        dp_deltomg(i, 1) = omgtop(i)/min(ptop(i)-ph(i,1),-smallestreal)
     3537      END IF
     3538    END DO
     3539
     3540    DO k = 1, klev
     3541      DO i = 1, klon
     3542        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
     3543          omg(i, k) = -rho(i, k)*RG*omg(i, k)
     3544          dp_deltomg(i, k) = dp_deltomg(i, 1)
     3545        END IF
     3546      END DO
     3547    END DO
     3548
     3549    ! raccordement lineaire de omg de ptop a pupper
     3550
     3551    DO i = 1, klon
     3552      IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
     3553        IF ( iflag_wk_profile == 0 ) THEN
     3554           omg(i, kupper(i)+1) =-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
     3555          RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
     3556        ELSE
     3557           omg(i, kupper(i)+1) = 0.
     3558        ENDIF
     3559        dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
     3560          (ptop(i)-pupper(i))
     3561      END IF
     3562    END DO
     3563
     3564    ! c      DO i=1,klon
     3565    ! c        print*,'Pente entre 0 et kupper (reference)'
     3566    ! c     $           ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
     3567    ! c        print*,'Pente entre ktop et kupper'
     3568    ! c     $   ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
     3569    ! c      ENDDO
     3570    ! c
     3571    DO k = 1, klev
     3572      DO i = 1, klon
     3573        IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
     3574          dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
     3575          omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
     3576        END IF
     3577      END DO
     3578    END DO
     3579!!    print *,'omg(igout,k) ', (k,omg(igout,k),k=1,klev)
     3580    ! cc nrlmd
     3581    ! c      DO i=1,klon
     3582    ! c      print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
     3583    ! c      END DO
     3584    ! cc
     3585
     3586
     3587    ! --    Compute wake average vertical velocity omgbw
     3588
     3589
     3590    DO k = 1, klev
     3591      DO i = 1, klon
     3592        IF (wk_adv(i)) THEN
     3593          omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
     3594        END IF
     3595      END DO
     3596    END DO
     3597    ! --    and its vertical gradient dp_omgbw
     3598
     3599    DO k = 1, klev-1
     3600      DO i = 1, klon
     3601        IF (wk_adv(i)) THEN
     3602          dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
     3603        END IF
     3604      END DO
     3605    END DO
     3606    DO i = 1, klon
     3607      IF (wk_adv(i)) THEN
     3608          dp_omgbw(i, klev) = 0.
     3609      END IF
     3610    END DO
     3611
     3612    ! --    Upstream coefficients for omgb velocity
     3613    ! --    (alpha_up(k) is the coefficient of the value at level k)
     3614    ! --    (1-alpha_up(k) is the coefficient of the value at level k-1)
     3615    DO k = 1, klev
     3616      DO i = 1, klon
     3617        IF (wk_adv(i)) THEN
     3618          alpha_up(i, k) = 0.
     3619          IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
     3620        END IF
     3621      END DO
     3622    END DO
     3623
     3624    ! Matrix expressing [The,deltatw] from  [Th1,Th2]
     3625
     3626    DO i = 1, klon
     3627      IF (wk_adv(i)) THEN
     3628        rre1(i) = 1. - sigmaw(i)
     3629        rre2(i) = sigmaw(i)
     3630      END IF
     3631    END DO
     3632    rrd1 = -1.
     3633    rrd2 = 1.
     3634
     3635    ! --    Get [Th1,Th2], dth and [q1,q2]
     3636
     3637    DO k = 1, klev
     3638      DO i = 1, klon
     3639        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
     3640          dth(i, k) = deltatw(i, k)/ppi(i, k)
     3641          th1(i, k) = thb(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
     3642          th2(i, k) = thb(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
     3643          q1(i, k) = qb(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
     3644          q2(i, k) = qb(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
     3645        END IF
     3646      END DO
     3647    END DO
     3648
     3649    DO i = 1, klon
     3650      IF (wk_adv(i)) THEN !!! nrlmd
     3651        d_th1(i, 1) = 0.
     3652        d_th2(i, 1) = 0.
     3653        d_dth(i, 1) = 0.
     3654        d_q1(i, 1) = 0.
     3655        d_q2(i, 1) = 0.
     3656        d_dq(i, 1) = 0.
     3657      END IF
     3658    END DO
     3659
     3660    DO k = 2, klev
     3661      DO i = 1, klon
     3662        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
     3663          d_th1(i, k) = th1(i, k-1) - th1(i, k)
     3664          d_th2(i, k) = th2(i, k-1) - th2(i, k)
     3665          d_dth(i, k) = dth(i, k-1) - dth(i, k)
     3666          d_q1(i, k) = q1(i, k-1) - q1(i, k)
     3667          d_q2(i, k) = q2(i, k-1) - q2(i, k)
     3668          d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
     3669        END IF
     3670      END DO
     3671    END DO
     3672
     3673    DO i = 1, klon
     3674      IF (wk_adv(i)) THEN
     3675        omgbdth(i, 1) = 0.
     3676        omgbdq(i, 1) = 0.
     3677      END IF
     3678    END DO
     3679
     3680    DO k = 2, klev
     3681      DO i = 1, klon
     3682        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN !   loop on interfaces
     3683          omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
     3684          omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
     3685        END IF
     3686      END DO
     3687    END DO
     3688
     3689!!    IF (prt_level>=10) THEN
     3690    IF (prt_level>=10 .and. wk_adv(igout)) THEN
     3691      PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,kupper(igout))
     3692      PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,kupper(igout))
     3693      PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,kupper(igout))
     3694      PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,kupper(igout))
     3695    ENDIF
     3696
     3697
     3698    ! -----------------------------------------------------------------
     3699          ! Compute redistribution (advective) term
     3700
     3701!     rate of change of sigmaw due to spreading
     3702          dsigspread(:) = gfl(:)*cstar(:)
     3703
     3704          CALL wake_dadv(klon, klev, dtimesub, ph, ppi, wk_adv, kupper, &
     3705                    omg, dp_deltomg, sigmaw, dsigspread, &
     3706                    th2, th1, q2, q1, &
     3707                    d_deltat_dadv, d_deltaq_dadv, d_tb_dadv, d_qb_dadv)
     3708
     3709    ! For the difference fields: convert to change per second in order to combine with the
     3710    ! other terms (d_deltat_ls, d_deltat_cv, d_deltat_gw)
     3711    d_deltat_dadv(:,:) = d_deltat_dadv(:,:)/dtimesub
     3712    d_deltaq_dadv(:,:) = d_deltaq_dadv(:,:)/dtimesub
     3713!
     3714    !   For the mean fields: tb and qb the computation of the tendencies due to wakes is
     3715    !   already complete.
     3716    d_tb(:,:) = d_tb_dadv(:,:)
     3717    d_qb(:,:) = d_qb_dadv(:,:)
     3718
     3719    ! -----------------------------------------------------------------
     3720    DO k = 1, klev-1
     3721      DO i = 1, klon
     3722        IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
     3723          ! -----------------------------------------------------------------
     3724          d_deltat_lsadv(i, k) = 1./(ph(i,k)-ph(i,k+1))* &
     3725            (-(1.-alpha_up(i,k))*omgbdth(i,k)- &
     3726             alpha_up(i,k+1)*omgbdth(i,k+1))*ppi(i, k)
     3727
     3728          d_deltaq_lsadv(i, k) = 1./(ph(i,k)-ph(i,k+1))* &
     3729            (-(1.-alpha_up(i,k))*omgbdq(i,k)- &
     3730             alpha_up(i,k+1)*omgbdq(i,k+1))
     3731
     3732        END IF
     3733        ! cc
     3734      END DO
     3735    END DO
     3736    ! ------------------------------------------------------------------
     3737
     3738!!    IF (prt_level>=10) THEN
     3739    IF (prt_level>=10 .and. wk_adv(igout)) THEN
     3740      PRINT *, 'wake-4.3, d_deltat_dadv(igout,k) ', (k,d_deltat_dadv(igout,k), k=1,klev)
     3741      PRINT *, 'wake-4.3, d_deltat_lsadv(igout,k) ', (k,d_deltat_lsadv(igout,k), k=1,klev)
     3742      PRINT *, 'wake-4.3, d_deltaq_dadv(igout,k) ', (k,d_deltaq_dadv(igout,k), k=1,klev)
     3743      PRINT *, 'wake-4.3, d_deltaq_lsadv(igout,k) ', (k,d_deltaq_lsadv(igout,k), k=1,klev)
     3744    ENDIF
     3745
     3746    ! Increment state variables
     3747!jyg<
     3748    IF (iflag_wk_pop_dyn >= 1) THEN
     3749      DO k = 1, klev
     3750        DO i = 1, klon
     3751          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     3752            detr(i,k) = - d_sig_death(i) - d_sig_col(i)     
     3753            entr_p(i,k) = d_sig_gen(i)
     3754          ENDIF
     3755        ENDDO
     3756      ENDDO
     3757      ELSE  ! (iflag_wk_pop_dyn >= 1)
     3758      DO k = 1, klev
     3759        DO i = 1, klon
     3760          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     3761            detr(i, k) = 0.
     3762   
     3763            entr_p(i, k) = 0.
     3764          ENDIF
     3765        ENDDO
     3766      ENDDO
     3767    ENDIF  ! (iflag_wk_pop_dyn >= 1)
     3768
     3769   
     3770
     3771    DO k = 1, klev
     3772      DO i = 1, klon
     3773        ! cc nrlmd       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
     3774        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     3775          ! cc
     3776
     3777
     3778
     3779          ! Coefficient de repartition
     3780
     3781          crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
     3782            (ph(i,kupper(i))-ph(i,1))
     3783          crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ &
     3784            (ph(i,1)-ph(i,kupper(i)))
     3785
     3786
     3787          ! Reintroduce compensating subsidence term.
     3788
     3789          ! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
     3790          ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
     3791          ! .                   /(1-sigmaw)
     3792          ! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
     3793          ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
     3794          ! .                   /(1-sigmaw)
     3795
     3796          ! dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
     3797          ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
     3798          ! .                   /(1-sigmaw)
     3799          ! dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
     3800          ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
     3801          ! .                   /(1-sigmaw)
     3802
     3803          !! Differential heating (d_deltat_dcv) and moistening (d_deltaq_dcv) by deep convection
     3804          d_deltat_dcv(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
     3805!!          dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))        ! supprime
     3806          d_deltaq_dcv(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
     3807!!          dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))        ! supprime
     3808          ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
     3809
     3810!
     3811
     3812          ! cc nrlmd          Prise en compte du taux de mortalite
     3813          ! cc               Definitions de entr, detr
     3814!jyg<
     3815!!            detr(i, k) = 0.
     3816!!   
     3817!!            entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
     3818!!              sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)
     3819!!
     3820            entr(i, k) = entr_p(i,k) + gfl(i)*cstar(i) + &
     3821                         sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)   
     3822            tgen(i,k) = entr_p(i,k)/sigmaw(i)
     3823!>jyg
     3824            wkspread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
     3825
     3826          ! cc        wkspread(i,k) =
     3827          ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
     3828          ! cc     $  sigmaw(i)
     3829
     3830
     3831          ! ajout d'un effet onde de gravite -Tgw(k)*deltatw(k) 03/02/06 YU
     3832          ! Jingmei
     3833
     3834          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
     3835          ! &  Tgw(i,k),deltatw(i,k)
     3836          d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)*deltatw(i, k)* &
     3837            dtimesub
     3838          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
     3839
     3840          ! Sans GW
     3841
     3842          ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-wkspread(k)*deltatw(k))
     3843
     3844          ! GW formule 1
     3845
     3846          ! deltatw(k) = deltatw(k)+dtimesub*
     3847          ! $         (ff+dtKE(k) - wkspread(k)*deltatw(k)-Tgw(k)*deltatw(k))
     3848
     3849          ! GW formule 2
     3850
     3851          !! Entrainment due to spread is supposed to be included in the differential advection
     3852          !! term (d_deltat_dadv); hence only the entrainment due to population dynamics (entr_p)
     3853          !! appears in the expression of d_deltatw.
     3854          IF (dtimesub*(tgw(i,k)+tgen(i,k))<1.E-10) THEN
     3855!!!            d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k) - &        ! nouvelle notation
     3856            d_deltatw(i, k) = dtimesub*(d_deltat_dadv(i,k)+d_deltat_lsadv(i,k)+d_deltat_dcv(i,k) - &
     3857               (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) - & ! cc
     3858               (tgw(i,k)+tgen(i,k))*deltatw(i,k) )
     3859          ELSE
     3860            d_deltatw(i, k) = 1/(tgw(i,k)+tgen(i,k))*(1-exp(-dtimesub*(tgw(i,k)+tgen(i,k))))* &
     3861!!!               (ff(i)+dtke(i,k) - &                                ! nouvelle notation
     3862               (d_deltat_dadv(i,k)+d_deltat_lsadv(i,k)+d_deltat_dcv(i,k) - &
     3863                (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) - &
     3864                (tgw(i,k)+tgen(i,k))*deltatw(i,k) )
     3865          END IF
     3866
     3867          dth(i, k) = deltatw(i, k)/ppi(i, k)
     3868
     3869          !! Entrainment due to spread is supposed to be included in the differential advection
     3870          !! term (d_deltaq_dadv); hence only the entrainment due to population dynamics (entr_p)
     3871          !! appears in the expression of d_deltaqw.
     3872          IF (dtimesub*tgen(i,k)<1.E-10) THEN
     3873            d_deltaqw(i, k) = dtimesub*(d_deltaq_dadv(i,k)+d_deltaq_lsadv(i,k)+d_deltaq_dcv(i,k) - &
     3874               (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k)/(1.-sigmaw(i)) - &
     3875               tgen(i,k)*deltaqw(i,k))
     3876          ELSE
     3877            d_deltaqw(i, k) = 1/tgen(i,k)*(1-exp(-dtimesub*tgen(i,k))) * &
     3878               (d_deltaq_dadv(i,k)+d_deltaq_lsadv(i,k)+d_deltaq_dcv(i,k) - &
     3879               (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k)/(1.-sigmaw(i)) - &
     3880               tgen(i,k)*deltaqw(i,k))
     3881          END IF
     3882          ! cc
     3883
     3884          ! cc nrlmd
     3885          ! cc       d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
     3886          ! cc       d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
     3887          ! cc
     3888        END IF
     3889      END DO
     3890    END DO
     3891
     3892!!    IF (prt_level>=10) THEN
     3893    IF (prt_level>=10 .and. wk_adv(igout)) THEN
     3894      PRINT *, 'wake-4.4, isubstep= ', isubstep,' deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev)
     3895      PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_deltat_dcv(igout,k) ', (k,d_deltat_dcv(igout,k), k=1,klev)
     3896      PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_deltat_dadv(igout,k) ', (k,d_deltat_dadv(igout,k), k=1,klev)
     3897      PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_deltat_lsadv(igout,k) ', (k,d_deltat_lsadv(igout,k), k=1,klev)
     3898      PRINT *, 'wake-4.4, isubstep= ', isubstep,' tgen(igout,k)*deltatw(igout,k) ', (k,tgen(igout,k)*deltatw(igout,k), k=1,klev)
     3899      PRINT *, 'wake-4.4, isubstep= ', isubstep,' tgw(igout,k)*deltatw(igout,k) ', (k,tgw(igout,k)*deltatw(igout,k), k=1,klev)
     3900      PRINT *, 'wake-4.4, isubstep= ', isubstep,' death_rate(igout) ', death_rate(igout)
     3901      PRINT *, 'wake-4.4, isubstep= ', isubstep,' detr(igout,k) ',  (k,detr(igout,k), k=1,klev)
     3902      PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_deltatw(igout,k) ', (k,d_deltatw(igout,k), k=1,klev)
     3903      PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_deltaqw(igout,k) ', (k,d_deltaqw(igout,k), k=1,klev)
     3904      PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_tb(igout,k) ', (k,d_tb(igout,k), k=1,klev)
     3905      PRINT *, 'wake-4.4, isubstep= ', isubstep,' d_qb(igout,k) ', (k,d_qb(igout,k), k=1,klev)
     3906    ENDIF
     3907
     3908!
     3909IF (CPPKEY_IOPHYS_WK) THEN
     3910    IF (phys_sub) THEN
     3911     DO k = 1,klev
     3912      d_deltat_entrp(:,k) = - entr_p(:,k)*deltatw(:,k)/sigmaw(:)
     3913     ENDDO
     3914     CALL iophys_ecrit('d_deltatw',klev,'d_deltatw','K/s',d_deltatw(:,1:klev))
     3915     CALL iophys_ecrit('d_deltat_dadv',klev,'d_deltat_dadv','K/s',d_deltat_dadv(:,1:klev))
     3916     CALL iophys_ecrit('d_deltat_lsadv',klev,'d_deltat_lsadv','K/s',d_deltat_lsadv(:,1:klev))
     3917     CALL iophys_ecrit('d_deltat_dcv',klev,'d_deltat_dcv','K/s',d_deltat_dcv(:,1:klev))
     3918     CALL iophys_ecrit('d_deltat_entrp',klev,'d_deltat_entrp','K/s',d_deltat_entrp(:,1:klev))
     3919
     3920    ENDIF
     3921END IF
     3922
     3923    ! Scale tendencies so that water vapour remains positive in w and x.
     3924
     3925    CALL wake_vec_modulation(klon, klev, wk_adv, epsilon_loc, qb, d_qb, deltaqw, &
     3926      d_deltaqw, sigmaw, d_sigmaw, alpha)
     3927    !
     3928    ! Alpha_tot = Product of all the alpha's
     3929    DO i = 1, klon
     3930      IF (wk_adv(i)) THEN
     3931        alpha_tot(i) = alpha_tot(i)*alpha(i)   
     3932      END IF
     3933    END DO
     3934
     3935    ! cc nrlmd
     3936    ! c      print*,'alpha'
     3937    ! c      do i=1,klon
     3938    ! c         print*,alpha(i)
     3939    ! c      end do
     3940    ! cc
     3941    DO k = 1, klev
     3942      DO i = 1, klon
     3943        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     3944          d_tb(i, k) = alpha(i)*d_tb(i, k)
     3945          d_qb(i, k) = alpha(i)*d_qb(i, k)
     3946          d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
     3947          d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
     3948          d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
     3949        END IF
     3950      END DO
     3951    END DO
     3952    DO i = 1, klon
     3953      IF (wk_adv(i)) THEN
     3954        d_sigmaw(i) = alpha(i)*d_sigmaw(i)
     3955      END IF
     3956    END DO
     3957
     3958    ! Update large scale variables and wake variables
     3959    ! IM 060208 manque DO i + remplace DO k=1,kupper(i)
     3960    ! IM 060208     DO k = 1,kupper(i)
     3961    DO k = 1, klev
     3962      DO i = 1, klon
     3963        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     3964          dtls(i, k) = dtls(i, k) + d_tb(i, k)
     3965          dqls(i, k) = dqls(i, k) + d_qb(i, k)
     3966          ! cc nrlmd
     3967          d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
     3968          d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
     3969          ! cc
     3970        END IF
     3971      END DO
     3972    END DO
     3973    DO k = 1, klev
     3974      DO i = 1, klon
     3975        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     3976          tb(i, k) = tb0(i, k) + dtls(i, k)
     3977          qb(i, k) = qb0(i, k) + dqls(i, k)
     3978          thb(i, k) = tb(i, k)/ppi(i, k)
     3979          deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
     3980          deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
     3981          dth(i, k) = deltatw(i, k)/ppi(i, k)
     3982          ! c      print*,'k,qx,qw',k,qb(i,k)-sigmaw(i)*deltaqw(i,k)
     3983          ! c     $        ,qb(i,k)+(1-sigmaw(i))*deltaqw(i,k)
     3984        END IF
     3985      END DO
     3986    END DO
     3987!
     3988    DO i = 1, klon
     3989      IF (wk_adv(i)) THEN
     3990        sigmaw(i) = sigmaw(i) + d_sigmaw(i)
     3991        d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
     3992      END IF
     3993    END DO
     3994!jyg<
     3995    IF (iflag_wk_pop_dyn >= 1) THEN
     3996!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! sigmaw !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     3997!  Cumulatives
     3998      DO i = 1, klon
     3999        IF (wk_adv(i)) THEN
     4000          d_sig_gen2(i)   = d_sig_gen2(i)   + d_sig_gen(i)
     4001          d_sig_death2(i) = d_sig_death2(i) + d_sig_death(i)
     4002          d_sig_col2(i)   = d_sig_col2(i)   + d_sig_col(i)
     4003          d_sig_spread2(i)= d_sig_spread2(i)+ d_sig_spread(i)
     4004          d_sig_bnd2(i)   = d_sig_bnd2(i)   + d_sig_bnd(i)
     4005        END IF
     4006      END DO
     4007!  Bounds
     4008      DO i = 1, klon
     4009        IF (wk_adv(i)) THEN
     4010          sigmaw_targ = max(sigmaw(i),sigmad)
     4011          d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     4012          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     4013          sigmaw(i) = sigmaw_targ
     4014        END IF
     4015      END DO
     4016!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! wdens  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     4017!  Cumulatives
     4018      DO i = 1, klon
     4019        IF (wk_adv(i)) THEN
     4020          wdens(i) = wdens(i) + d_wdens(i)
     4021          d_wdens2(i) = d_wdens2(i) + d_wdens(i)
     4022          d_dens_gen2(i)   = d_dens_gen2(i)   + d_dens_gen(i)
     4023          d_dens_death2(i) = d_dens_death2(i) + d_dens_death(i)
     4024          d_dens_col2(i)   = d_dens_col2(i)   + d_dens_col(i)
     4025          d_dens_bnd2(i)   = d_dens_bnd2(i)   + d_dens_bnd(i)
     4026        END IF
     4027      END DO
     4028!  Bounds
     4029      DO i = 1, klon
     4030        IF (wk_adv(i)) THEN
     4031          wdens_targ = max(wdens(i),wdensmin)
     4032          d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
     4033          d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
     4034          wdens(i) = wdens_targ
     4035        END IF
     4036      END DO
     4037!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     4038!  Cumulatives
     4039      DO i = 1, klon
     4040        IF (wk_adv(i)) THEN
     4041          awdens(i) = awdens(i) + d_awdens(i)
     4042          d_awdens2(i) = d_awdens2(i) + d_awdens(i)
     4043        END IF
     4044      END DO
     4045!  Bounds
     4046      DO i = 1, klon
     4047        IF (wk_adv(i)) THEN
     4048          wdens_targ = min( max(awdens(i),0.), wdens(i) )
     4049          d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     4050          d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
     4051          awdens(i) = wdens_targ
     4052        END IF
     4053      END DO
     4054!
     4055      IF (iflag_wk_pop_dyn >= 2) THEN
     4056!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens again for iflag_wk_pop_dyn >= 2!!!!!!
     4057!  Cumulatives
     4058        DO i = 1, klon
     4059           IF (wk_adv(i)) THEN
     4060               d_adens_death2(i)   = d_adens_death2(i)   + d_adens_death(i)
     4061               d_adens_icol2(i)   = d_adens_icol2(i)   + d_adens_icol(i)
     4062               d_adens_acol2(i)   = d_adens_acol2(i)   + d_adens_acol(i)
     4063               d_adens_bnd2(i)   = d_adens_bnd2(i)   + d_adens_bnd(i)         
     4064           END IF
     4065        END DO
     4066!  Bounds
     4067        DO i = 1, klon
     4068           IF (wk_adv(i)) THEN
     4069               wdens_targ = min( max(awdens(i),0.), wdens(i) )
     4070               d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     4071               awdens(i) = wdens_targ
     4072           END IF
     4073        END DO
     4074!
     4075        IF (iflag_wk_pop_dyn == 3) THEN
     4076!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! asigmaw for iflag_wk_pop_dyn = 3!!!!!!
     4077!  Cumulatives
     4078          DO i = 1, klon
     4079             IF (wk_adv(i)) THEN
     4080                 asigmaw(i) = asigmaw(i) + d_asigmaw(i)
     4081                 d_asigmaw2(i) = d_asigmaw2(i) + d_asigmaw(i)
     4082                 d_asig_death2(i)   = d_asig_death2(i)   + d_asig_death(i)
     4083                 d_asig_spread2(i)  = d_asig_spread2(i)  + d_asig_spread(i)
     4084                 d_asig_iicol2(i)   = d_asig_iicol2(i)   + d_asig_iicol(i)
     4085                 d_asig_aicol2(i)   = d_asig_aicol2(i)   + d_asig_aicol(i)
     4086                 d_asig_bnd2(i)     = d_asig_bnd2(i)     + d_asig_bnd(i)         
     4087             END IF
     4088          END DO
     4089!  Bounds
     4090          DO i = 1, klon
     4091             IF (wk_adv(i)) THEN
     4092   !   asigmaw lower bound set to sigmad/2 in order to allow asigmaw values lower than sigmad.
     4093   !!             sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i))
     4094                sigmaw_targ = min(max(asigmaw(i),sigmad/2.),sigmaw(i))
     4095                d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     4096                d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     4097                asigmaw(i) = sigmaw_targ
     4098             END IF
     4099          END DO
     4100
     4101IF (CPPKEY_IOPHYS_WK) THEN
     4102    IF (phys_sub) THEN
     4103     CALL iophys_ecrit('wdensb',1,'wdensb','m',wdens)
     4104     CALL iophys_ecrit('awdensb',1,'awdensb','m',awdens)
     4105     CALL iophys_ecrit('sigmawb',1,'sigmawb','m',sigmaw)
     4106     CALL iophys_ecrit('asigmawb',1,'asigmawb','m',asigmaw)
     4107!
     4108     call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2)
     4109     call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2)
     4110     call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2)
     4111     call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2)
     4112     call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2)
     4113!
     4114     call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2)
     4115     call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2)
     4116     call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2)
     4117     call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2)
     4118     call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2)
     4119!
     4120     CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2)
     4121     CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2)
     4122     CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2)
     4123     CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2)
     4124     CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2)
     4125     CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2)
     4126!
     4127     CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2)
     4128     CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2)
     4129     CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2)
     4130     CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2)
     4131     CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2)
     4132     CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2)
     4133    ENDIF
     4134END IF
     4135        ENDIF ! (iflag_wk_pop_dyn == 3)
     4136      ENDIF ! (iflag_wk_pop_dyn >= 2)
     4137    ENDIF  ! (iflag_wk_pop_dyn >= 1)
     4138
     4139
     4140
     4141   Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
     4142                    dth, hw, rho, delta_t_min, &
     4143                    ktop, wk_adv, h_zzz, ptop1, ktop1)
     4144   !! print'("pkupper APPEL ",7i6)',isubstep,int(ptop/100.),int(ptop1/100.),int(pupper/100.),ktop,ktop1,kupper
     4145
     4146    ! 5/ Set deltatw & deltaqw to 0 above kupper
     4147
     4148    DO k = 1, klev
     4149      DO i = 1, klon
     4150        IF (wk_adv(i) .AND. k>=kupper(i)) THEN
     4151          deltatw(i, k) = 0.
     4152          deltaqw(i, k) = 0.
     4153          d_deltatw2(i,k) = -deltatw0(i,k)
     4154          d_deltaqw2(i,k) = -deltaqw0(i,k)
     4155        END IF
     4156      END DO
     4157    END DO
     4158
     4159
     4160    ! -------------Cstar computation---------------------------------
     4161    DO i = 1, klon
     4162      IF (wk_adv(i)) THEN !!! nrlmd
     4163        sum_thx(i) = 0.
     4164        sum_tx(i) = 0.
     4165        sum_qx(i) = 0.
     4166        sum_thvx(i) = 0.
     4167        sum_dth(i) = 0.
     4168        sum_dq(i) = 0.
     4169        sum_dtdwn(i) = 0.
     4170        sum_dqdwn(i) = 0.
     4171
     4172        av_thx(i) = 0.
     4173        av_tx(i) = 0.
     4174        av_qx(i) = 0.
     4175        av_thvx(i) = 0.
     4176        av_dth(i) = 0.
     4177        av_dq(i) = 0.
     4178        av_dtdwn(i) = 0.
     4179        av_dqdwn(i) = 0.
     4180      END IF
     4181    END DO
     4182
     4183    ! Integrals (and wake top level number)
     4184    ! --------------------------------------
     4185
     4186    ! Initialize sum_thvx to 1st level virt. pot. temp.
     4187
     4188    DO i = 1, klon
     4189      IF (wk_adv(i)) THEN !!! nrlmd
     4190        z(i) = 1.
     4191        dz(i) = 1.
     4192        sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i)
     4193        sum_dth(i) = 0.
     4194      END IF
     4195    END DO
     4196
     4197    DO k = 1, klev
     4198      DO i = 1, klon
     4199        IF (wk_adv(i)) THEN !!! nrlmd
     4200          dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
     4201          IF (dz(i)>0) THEN
     4202            z(i) = z(i) + dz(i)
     4203            sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i)
     4204            sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i)
     4205            sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i)
     4206            sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i)
     4207            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
     4208            sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
     4209            sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
     4210            sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
     4211          END IF
     4212        END IF
     4213      END DO
     4214    END DO
     4215
     4216    DO i = 1, klon
     4217      IF (wk_adv(i)) THEN !!! nrlmd
     4218        hw0(i) = z(i)
     4219      END IF
     4220    END DO
     4221
     4222
     4223    ! - WAPE and mean forcing computation
     4224    ! ---------------------------------------
     4225
     4226    ! ---------------------------------------
     4227
     4228    ! Means
     4229
     4230    DO i = 1, klon
     4231      IF (wk_adv(i)) THEN !!! nrlmd
     4232        av_thx(i) = sum_thx(i)/hw0(i)
     4233        av_tx(i) = sum_tx(i)/hw0(i)
     4234        av_qx(i) = sum_qx(i)/hw0(i)
     4235        av_thvx(i) = sum_thvx(i)/hw0(i)
     4236        av_dth(i) = sum_dth(i)/hw0(i)
     4237        av_dq(i) = sum_dq(i)/hw0(i)
     4238        av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     4239        av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     4240
     4241        wape(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_thx(i)*av_dq(i) + &
     4242                              av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i)
     4243!!print *,'XXXXwake wape(i), hw0(i), av_dth(i), av_thx(i), av_dq(i), av_qx(i), av_thvx(i) ', &
     4244!!                  wape(i), hw0(i), av_dth(i), av_thx(i), av_dq(i), av_qx(i), av_thvx(i)
     4245      END IF
     4246    END DO
     4247
     4248
     4249    ! Filter out bad wakes
     4250
     4251    DO k = 1, klev
     4252      DO i = 1, klon
     4253        IF (wk_adv(i)) THEN !!! nrlmd
     4254          IF (wape(i)<0.) THEN
     4255            deltatw(i, k) = 0.
     4256            deltaqw(i, k) = 0.
     4257            dth(i, k) = 0.
     4258            d_deltatw2(i,k) = -deltatw0(i,k)
     4259            d_deltaqw2(i,k) = -deltaqw0(i,k)
     4260          END IF
     4261        END IF
     4262      END DO
     4263    END DO
     4264
     4265    DO i = 1, klon
     4266      IF (wk_adv(i)) THEN !!! nrlmd
     4267        IF (wape(i)<0.) THEN
     4268          wape(i) = 0.
     4269          cstar(i) = 0.
     4270          hw(i) = hwmin
     4271!jyg<
     4272!!          sigmaw(i) = max(sigmad, sigd_con(i))
     4273          sigmaw_targ = max(sigmad, sigd_con(i))
     4274          d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     4275          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     4276          sigmaw(i) = sigmaw_targ
     4277!
     4278          d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     4279          d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     4280          asigmaw(i) = sigmaw_targ
     4281!>jyg
     4282          fip(i) = 0.
     4283          gwake(i) = .FALSE.
     4284        ELSE
     4285          cstar(i) = stark*sqrt(2.*wape(i))
     4286          gwake(i) = .TRUE.
     4287        END IF
     4288      END IF
     4289    END DO
     4290  !
     4291  ! ------------------------------------------------------------------------
     4292  !
     4293  END DO   ! isubstep end sub-timestep loop
     4294  !
     4295  ! ------------------------------------------------------------------------
     4296  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     4297  ! ------------------------------------------------------------------------
     4298  !
     4299
     4300IF (CPPKEY_IOPHYS_WK) THEN
     4301    IF (.not.phys_sub) CALL iophys_ecrit('wape_b',1,'wape_b','J/kg',wape)
     4302    IF (.not.phys_sub) CALL iophys_ecrit('alpha_mod',1,'alpha_modulation','-',alpha_tot)
     4303END IF
     4304  IF (prt_level>=10) THEN
     4305    PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', &
     4306                      sigmaw(igout), cstar(igout), wape(igout), ptop(igout)
     4307  ENDIF
     4308
     4309
     4310  ! ----------------------------------------------------------
     4311  ! Determine wake final state; recompute wape, cstar, ktop;
     4312  ! filter out bad wakes.
     4313  ! ----------------------------------------------------------
     4314
     4315  ! 2.1 - Undisturbed area and Wake integrals
     4316  ! ---------------------------------------------------------
     4317
     4318  DO i = 1, klon
     4319    ! cc nrlmd       if (wk_adv(i)) then !!! nrlmd
     4320    IF (ok_qx_qw(i)) THEN
     4321      ! cc
     4322      z(i) = 0.
     4323      sum_thx(i) = 0.
     4324      sum_tx(i) = 0.
     4325      sum_qx(i) = 0.
     4326      sum_thvx(i) = 0.
     4327      sum_dth(i) = 0.
     4328      sum_half_dth(i) = 0.
     4329      sum_dq(i) = 0.
     4330      sum_dtdwn(i) = 0.
     4331      sum_dqdwn(i) = 0.
     4332
     4333      av_thx(i) = 0.
     4334      av_tx(i) = 0.
     4335      av_qx(i) = 0.
     4336      av_thvx(i) = 0.
     4337      av_dth(i) = 0.
     4338      av_dq(i) = 0.
     4339      av_dtdwn(i) = 0.
     4340      av_dqdwn(i) = 0.
     4341
     4342      dthmin(i) = -delta_t_min
     4343    END IF
     4344  END DO
     4345  ! Potential temperatures and humidity
     4346  ! ----------------------------------------------------------
     4347
     4348  DO k = 1, klev
     4349    DO i = 1, klon
     4350      ! cc nrlmd       IF ( wk_adv(i)) THEN
     4351      IF (ok_qx_qw(i)) THEN
     4352        ! cc
     4353        rho(i, k) = p(i, k)/(RD*tb(i,k))
     4354        IF (k==1) THEN
     4355          rhoh(i, k) = ph(i, k)/(RD*tb(i,k))
     4356          zhh(i, k) = 0
     4357        ELSE
     4358          rhoh(i, k) = ph(i, k)*2./(RD*(tb(i,k)+tb(i,k-1)))
     4359          zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
     4360        END IF
     4361        thb(i, k) = tb(i, k)/ppi(i, k)
     4362        thx(i, k) = (tb(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
     4363        tx(i, k) = tb(i, k) - deltatw(i, k)*sigmaw(i)
     4364        qx(i, k) = qb(i, k) - deltaqw(i, k)*sigmaw(i)
     4365        dth(i, k) = deltatw(i, k)/ppi(i, k)
     4366      END IF
     4367    END DO
     4368  END DO
     4369
     4370  ! Integrals (and wake top level number)
     4371  ! -----------------------------------------------------------
     4372
     4373  ! Initialize sum_thvx to 1st level virt. pot. temp.
     4374
     4375  DO i = 1, klon
     4376    ! cc nrlmd       IF ( wk_adv(i)) THEN
     4377    IF (ok_qx_qw(i)) THEN
     4378      ! cc
     4379      z(i) = 1.
     4380      dz(i) = 1.
     4381      dz_half(i) = 1.
     4382      sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i)
     4383      sum_dth(i) = 0.
     4384    END IF
     4385  END DO
     4386
     4387  DO k = 1, klev
     4388    DO i = 1, klon
     4389      ! cc nrlmd       IF ( wk_adv(i)) THEN
     4390      IF (ok_qx_qw(i)) THEN
     4391        ! cc
     4392        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
     4393        dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)*RG)
     4394        IF (dz(i)>0) THEN
     4395          z(i) = z(i) + dz(i)
     4396          sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i)
     4397          sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i)
     4398          sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i)
     4399          sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i)
     4400          sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
     4401          sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
     4402          sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
     4403          sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
     4404!
     4405          dthmin(i) = min(dthmin(i), dth(i,k))
     4406        END IF
     4407        IF (dz_half(i)>0) THEN
     4408          sum_half_dth(i) = sum_half_dth(i) + dth(i, k)*dz_half(i)
     4409        END IF
     4410      END IF
     4411    END DO
     4412  END DO
     4413
     4414  DO i = 1, klon
     4415    ! cc nrlmd       IF ( wk_adv(i)) THEN
     4416    IF (ok_qx_qw(i)) THEN
     4417      ! cc
     4418      hw0(i) = z(i)
     4419    END IF
     4420  END DO
     4421
     4422  ! - WAPE and mean forcing computation
     4423  ! -------------------------------------------------------------
     4424
     4425  ! Means
     4426
     4427  DO i = 1, klon
     4428    ! cc nrlmd       IF ( wk_adv(i)) THEN
     4429    IF (ok_qx_qw(i)) THEN
     4430      ! cc
     4431      av_thx(i) = sum_thx(i)/hw0(i)
     4432      av_tx(i) = sum_tx(i)/hw0(i)
     4433      av_qx(i) = sum_qx(i)/hw0(i)
     4434      av_thvx(i) = sum_thvx(i)/hw0(i)
     4435      av_dth(i) = sum_dth(i)/hw0(i)
     4436      av_dq(i) = sum_dq(i)/hw0(i)
     4437      av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     4438      av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     4439
     4440      wape2(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_thx(i)*av_dq(i) + &
     4441                             av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i)
     4442    END IF
     4443  END DO
     4444IF (CPPKEY_IOPHYS_WK) THEN
     4445  IF (.not.phys_sub) CALL iophys_ecrit('wape2_a',1,'wape2_a','J/kg',wape2)
     4446END IF
     4447
     4448
     4449  ! Prognostic variable update
     4450  ! ------------------------------------------------------------
     4451
     4452  ! Filter out bad wakes
     4453
     4454  IF (iflag_wk_check_trgl>=1) THEN
     4455    ! Check triangular shape of dth profile
     4456    DO i = 1, klon
     4457      IF (ok_qx_qw(i)) THEN
     4458        !! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
     4459        !! print *,'wake, 2.*sum_dth(i)/(hw0(i)*dthmin(i)) ', &
     4460        !!                2.*sum_dth(i)/(hw0(i)*dthmin(i))
     4461        !! print *,'wake, sum_half_dth(i), sum_dth(i) ', &
     4462        !!                sum_half_dth(i), sum_dth(i)
     4463        IF ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min) ) THEN
     4464          wape2(i) = -1.
     4465          !! print *,'wake, rej 1'
     4466        ELSE IF (iflag_wk_check_trgl==1.AND.abs(2.*sum_dth(i)/(hw0(i)*dthmin(i)) - 1.) > 0.5) THEN
     4467          wape2(i) = -1.
     4468          !! print *,'wake, rej 2'
     4469        ELSE IF (abs(sum_half_dth(i)) < 0.5*abs(sum_dth(i)) ) THEN
     4470          wape2(i) = -1.
     4471          !! print *,'wake, rej 3'
     4472        END IF
     4473      END IF
     4474    END DO
     4475  END IF
     4476IF (CPPKEY_IOPHYS_WK) THEN
     4477  IF (.not.phys_sub) CALL iophys_ecrit('wape2_b',1,'wape2_b','J/kg',wape2)
     4478END IF
     4479
     4480
     4481  DO k = 1, klev
     4482    DO i = 1, klon
     4483      ! cc nrlmd        IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
     4484      IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
     4485        ! cc
     4486        deltatw(i, k) = 0.
     4487        deltaqw(i, k) = 0.
     4488        dth(i, k) = 0.
     4489        d_deltatw2(i,k) = -deltatw0(i,k)
     4490        d_deltaqw2(i,k) = -deltaqw0(i,k)
     4491      END IF
     4492    END DO
     4493  END DO
     4494
     4495
     4496  DO i = 1, klon
     4497    ! cc nrlmd       IF ( wk_adv(i)) THEN
     4498    IF (ok_qx_qw(i)) THEN
     4499      ! cc
     4500      IF (wape2(i)<0.) THEN
     4501        wape2(i) = 0.
     4502        cstar2(i) = 0.
     4503        hw(i) = hwmin
     4504!jyg<
     4505!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
     4506      sigmaw_targ = max(sigmad, sigd_con(i))
     4507      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     4508      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     4509      sigmaw(i) = sigmaw_targ
     4510!
     4511      d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     4512      d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
     4513      asigmaw(i) = sigmaw_targ
     4514!>jyg
     4515        fip(i) = 0.
     4516        gwake(i) = .FALSE.
     4517      ELSE
     4518        IF (prt_level>=10) PRINT *, 'wape2>0'
     4519        cstar2(i) = stark*sqrt(2.*wape2(i))
     4520        gwake(i) = .TRUE.
     4521      END IF
     4522IF (CPPKEY_IOPHYS_WK) THEN
     4523  IF (.not.phys_sub) CALL iophys_ecrit('cstar2',1,'cstar2','J/kg',cstar2)
     4524END IF
     4525    END IF  ! (ok_qx_qw(i))
     4526  END DO
     4527
     4528  DO i = 1, klon
     4529    ! cc nrlmd       IF ( wk_adv(i)) THEN
     4530    IF (ok_qx_qw(i)) THEN
     4531      ! cc
     4532      ktopw(i) = ktop(i)
     4533    END IF
     4534  END DO
     4535
     4536  DO i = 1, klon
     4537    ! cc nrlmd       IF ( wk_adv(i)) THEN
     4538    IF (ok_qx_qw(i)) THEN
     4539      ! cc
     4540      IF (ktopw(i)>0 .AND. gwake(i)) THEN
     4541
     4542        ! jyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
     4543        ! cc       heff = 600.
     4544        ! Utilisation de la hauteur hw
     4545        ! c       heff = 0.7*hw
     4546        heff(i) = hw(i)
     4547
     4548        fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* &
     4549          sqrt(sigmaw(i)*wdens(i)*3.14)
     4550        fip(i) = alpk*fip(i)
     4551        ! jyg2
     4552      ELSE
     4553        fip(i) = 0.
     4554      END IF
     4555    END IF
     4556  END DO
     4557    IF (iflag_wk_pop_dyn >= 3) THEN
     4558IF (CPPKEY_IOPHYS_WK) THEN
     4559      IF (.not.phys_sub) THEN
     4560       call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2)
     4561       call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2)
     4562       call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2)
     4563       call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2)
     4564       call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2)
     4565!   
     4566       call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2)
     4567       call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2)
     4568       call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2)
     4569       call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2)
     4570       call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2)
     4571!   
     4572       CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2)
     4573       CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2)
     4574       CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2)
     4575       CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2)
     4576       CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2)
     4577       CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2)
     4578!   
     4579       CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2)
     4580       CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2)
     4581       CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2)
     4582       CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2)
     4583       CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2)
     4584       CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2)
     4585      ENDIF  ! (.not.phys_sub)
     4586END IF
     4587    ENDIF  ! (iflag_wk_pop_dyn >= 3)
     4588  ! Limitation de sigmaw
     4589
     4590  ! cc nrlmd
     4591  ! DO i=1,klon
     4592  ! IF (OK_qx_qw(i)) THEN
     4593  ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
     4594  ! ENDIF
     4595  ! ENDDO
     4596  ! cc
     4597
     4598  !jyg<
     4599  IF (iflag_wk_pop_dyn >= 1) THEN
     4600    DO i = 1, klon
     4601      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     4602          .NOT. ok_qx_qw(i) .OR. (wdens(i) < wdensthreshold)
     4603!!          .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
     4604    ENDDO
     4605  ELSE  ! (iflag_wk_pop_dyn >= 1)
     4606    DO i = 1, klon
     4607      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     4608          .NOT. ok_qx_qw(i)
     4609    ENDDO
     4610  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     4611  !>jyg
     4612
     4613  DO k = 1, klev
     4614    DO i = 1, klon
     4615!!jyg      IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     4616!!jyg          .NOT. ok_qx_qw(i)) THEN
     4617      IF (kill_wake(i)) THEN
     4618        ! cc
     4619        dtls(i, k) = 0.
     4620        dqls(i, k) = 0.
     4621        deltatw(i, k) = 0.
     4622        deltaqw(i, k) = 0.
     4623        d_deltatw2(i,k) = -deltatw0(i,k)
     4624        d_deltaqw2(i,k) = -deltaqw0(i,k)
     4625      END IF  ! (kill_wake(i))
     4626    END DO
     4627  END DO
     4628
     4629  DO i = 1, klon
     4630!!jyg    IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
     4631!!jyg        .NOT. ok_qx_qw(i)) THEN
     4632      IF (kill_wake(i)) THEN
     4633      ktopw(i) = 0
     4634      wape(i) = 0.
     4635      cstar(i) = 0.
     4636!!jyg   Outside subroutine "Wake" hw, wdens sigmaw and asigmaw are zero when there are no wakes
     4637!!      hw(i) = hwmin                       !jyg
     4638!!      sigmaw(i) = sigmad                  !jyg
     4639      hw(i) = 0.                            !jyg
     4640      fip(i) = 0.
     4641!
     4642!!      sigmaw(i) = 0.                        !jyg
     4643      sigmaw_targ = 0.
     4644      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
     4645!!      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     4646      d_sigmaw2(i) = sigmaw_targ - sigmaw_in(i)      ! _in = correction jyg 20220124
     4647      sigmaw(i) = sigmaw_targ
     4648!
     4649      IF (iflag_wk_pop_dyn >= 3) THEN
     4650        sigmaw_targ = 0.
     4651        d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
     4652!!        d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
     4653        d_asigmaw2(i) = sigmaw_targ - asigmaw_in(i)      ! _in = correction jyg 20220124
     4654        asigmaw(i) = sigmaw_targ
     4655      ELSE
     4656        asigmaw(i) = 0.
     4657      ENDIF ! (iflag_wk_pop_dyn >= 3)
     4658!
     4659      IF (iflag_wk_pop_dyn >= 1) THEN
     4660!!        awdens(i) = 0.
     4661!!        wdens(i) = 0.
     4662        wdens_targ = 0.
     4663        d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
     4664!!        d_wdens2(i) = wdens_targ - wdens(i)
     4665        d_wdens2(i) = wdens_targ - wdens_in(i)      ! jyg 20220916
     4666        wdens(i) = wdens_targ
     4667        wdens_targ = 0.
     4668!!jyg: bug fix : the d_adens_bnd2 computation must be before the update of awdens.
     4669        IF (iflag_wk_pop_dyn >= 2) THEN
     4670            d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     4671        ENDIF ! (iflag_wk_pop_dyn >= 2)
     4672!!        d_awdens2(i) = wdens_targ - awdens(i)
     4673        d_awdens2(i) = wdens_targ - awdens_in(i)    ! jyg 20220916
     4674        awdens(i) = wdens_targ
     4675!!        IF (iflag_wk_pop_dyn == 2) THEN
     4676!!            d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
     4677!!        ENDIF ! (iflag_wk_pop_dyn == 2)
     4678      ENDIF  ! (iflag_wk_pop_dyn >= 1)
     4679    ELSE  ! (kill_wake(i))
     4680      wape(i) = wape2(i)
     4681      cstar(i) = cstar2(i)
     4682    END IF  ! (kill_wake(i))
     4683    ! c        print*,'wape wape2 ktopw OK_qx_qw =',
     4684    ! c     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
     4685  END DO
     4686
     4687  IF (prt_level>=10) THEN
     4688    PRINT *, 'wake-6, wape wape2 ktopw OK_qx_qw =', &
     4689                      wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout)
     4690  ENDIF
     4691IF (CPPKEY_IOPHYS_WK) THEN
     4692  IF (.not.phys_sub) CALL iophys_ecrit('wape_c',1,'wape_c','J/kg',wape)
     4693  IF (iflag_wk_pop_dyn >= 3) THEN
     4694    IF (.not.phys_sub) THEN
     4695     CALL iophys_ecrit('fip',1,'fip','J/kg',fip)
     4696     CALL iophys_ecrit('hw',1,'hw','J/kg',hw)
     4697     CALL iophys_ecrit('ptop',1,'ptop','J/kg',ptop)
     4698     CALL iophys_ecrit('wdens',1,'wdens','J/kg',wdens)
     4699     CALL iophys_ecrit('awdens',1,'awdens','m',awdens)
     4700     CALL iophys_ecrit('sigmaw',1,'sigmaw','m',sigmaw)
     4701     CALL iophys_ecrit('asigmaw',1,'asigmaw','m',asigmaw)
     4702!
     4703     CALL iophys_ecrit('rad_wk',1,'rad_wk','J/kg',rad_wk)
     4704     CALL iophys_ecrit('arad_wk',1,'arad_wk','J/kg',arad_wk)
     4705     CALL iophys_ecrit('irad_wk',1,'irad_wk','J/kg',irad_wk)
     4706    ENDIF  ! (.not.phys_sub)
     4707  ENDIF  ! (iflag_wk_pop_dyn >= 3)
     4708END IF  !(CPPKEY_IOPHYS_WK)
     4709
     4710
     4711  ! -----------------------------------------------------------------
     4712  ! Get back to tendencies per second
     4713
     4714  DO k = 1, klev
     4715    DO i = 1, klon
     4716
     4717      ! cc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
     4718!jyg<
     4719!!      IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
     4720      IF (ok_qx_qw(i)) THEN
     4721!>jyg
     4722        ! cc
     4723        dtls(i, k) = dtls(i, k)/dtime
     4724        dqls(i, k) = dqls(i, k)/dtime
     4725        d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
     4726        d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
     4727        d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
     4728        ! c      print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
     4729        ! c     $         ,death_rate(i)*sigmaw(i)
     4730      END IF
     4731    END DO
     4732  END DO
     4733!jyg<
     4734  IF (iflag_wk_pop_dyn >= 1) THEN
     4735    DO i = 1, klon
     4736        IF (ok_qx_qw(i)) THEN
     4737      d_sig_gen2(i) = d_sig_gen2(i)/dtime
     4738      d_sig_death2(i) = d_sig_death2(i)/dtime
     4739      d_sig_col2(i) = d_sig_col2(i)/dtime
     4740      d_sig_spread2(i) = d_sig_spread2(i)/dtime
     4741      d_sig_bnd2(i) = d_sig_bnd2(i)/dtime
     4742      d_sigmaw2(i) = d_sigmaw2(i)/dtime
     4743!
     4744      d_dens_gen2(i) = d_dens_gen2(i)/dtime
     4745      d_dens_death2(i) = d_dens_death2(i)/dtime
     4746      d_dens_col2(i) = d_dens_col2(i)/dtime
     4747      d_dens_bnd2(i) = d_dens_bnd2(i)/dtime
     4748      d_awdens2(i) = d_awdens2(i)/dtime
     4749      d_wdens2(i) = d_wdens2(i)/dtime
     4750        ENDIF
     4751    ENDDO
     4752    IF (iflag_wk_pop_dyn >= 2) THEN
     4753      DO i = 1, klon
     4754        IF (ok_qx_qw(i)) THEN
     4755        d_adens_death2(i) = d_adens_death2(i)/dtime
     4756        d_adens_icol2(i) = d_adens_icol2(i)/dtime
     4757        d_adens_acol2(i) = d_adens_acol2(i)/dtime
     4758        d_adens_bnd2(i) = d_adens_bnd2(i)/dtime
     4759        ENDIF
     4760      ENDDO
     4761      IF (iflag_wk_pop_dyn == 3) THEN
     4762       DO i = 1, klon
     4763          IF (ok_qx_qw(i)) THEN
     4764        d_asig_death2(i)  = d_asig_death2(i)/dtime
     4765        d_asig_iicol2(i)  = d_asig_iicol2(i)/dtime
     4766        d_asig_aicol2(i)  = d_asig_aicol2(i)/dtime
     4767        d_asig_spread2(i) = d_asig_spread2(i)/dtime
     4768        d_asig_bnd2(i) = d_asig_bnd2(i)/dtime
     4769          ENDIF
     4770       ENDDO
     4771      ENDIF ! (iflag_wk_pop_dyn == 3) 
     4772    ENDIF ! (iflag_wk_pop_dyn >= 2) 
     4773  ENDIF  ! (iflag_wk_pop_dyn >= 1)
     4774 
     4775!>jyg
     4776
     4777 RETURN
     4778END SUBROUTINE wake2
     4779
    23554780SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon_loc, qb, d_qb, deltaqw, &
    23564781    d_deltaqw, sigmaw, d_sigmaw, alpha)
    … …  
    23634788
    23644789  ! Input
    2365   REAL qb(nlon, nl), d_qb(nlon, nl)
    2366   REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
    2367   REAL sigmaw(nlon), d_sigmaw(nlon)
    2368   LOGICAL wk_adv(nlon)
    2369   INTEGER nl, nlon
     4790  INTEGER,                      INTENT(IN)               :: nl, nlon
     4791  REAL, DIMENSION(nlon, nl),    INTENT(IN)               :: qb, d_qb
     4792  REAL, DIMENSION(nlon, nl),    INTENT(IN)               :: deltaqw, d_deltaqw
     4793  REAL, DIMENSION(nlon),        INTENT(IN)               :: sigmaw, d_sigmaw
     4794  LOGICAL, DIMENSION(nlon),     INTENT(IN)               :: wk_adv
    23704795  ! Output
    2371   REAL alpha(nlon)
     4796  REAL, DIMENSION(nlon),        INTENT(INOUT)            :: alpha
    23724797  ! Internal variables
    23734798  REAL zeta(nlon, nl)
    … …  
    32025627  REAL, DIMENSION (klon)                                :: is_wk           !! mean area of individual inactive wakes
    32035628  REAL, DIMENSION (klon)                                :: tau_wk_inv      !! tau = life time of wakes
    3204   REAL                                                  :: tau_wk_inv_min
    32055629  REAL, DIMENSION (klon)                                :: tau_prime       !! tau_prime = life time of actives wakes
    32065630  REAL                                                  :: d_wdens_targ, d_sigmaw_targ
    … …  
    32295653!!
    32305654
     5655! Initialization
     5656 tau_wk_inv(:) = 0.
     5657! Initialization of output variables
     5658 rad_wk(:) = 0.
     5659 arad_wk(:) = 0.
     5660 irad_wk(:) = 0.
     5661 gfl(:) = 0.
     5662 agfl(:) = 0.
     5663!
     5664 d_wdens(:) = 0.
     5665 d_awdens(:) = 0.
     5666 d_sigmaw(:) = 0.
     5667 d_asigmaw (:) = 0.
     5668!
     5669 d_sig_gen(:) = 0.
     5670 d_sig_death(:) = 0.
     5671 d_sig_col(:) = 0.
     5672 d_sig_spread(:) = 0.
     5673 d_sig_bnd(:) = 0.
     5674 d_asig_death(:) = 0.
     5675 d_asig_aicol(:) = 0.
     5676 d_asig_iicol(:) = 0.
     5677 d_asig_spread(:) = 0.
     5678 d_asig_bnd(:) = 0.
     5679 d_dens_gen(:) = 0.
     5680 d_dens_death(:) = 0.
     5681 d_dens_col(:) = 0.
     5682 d_dens_bnd(:) = 0.
     5683 d_adens_death(:) = 0.
     5684 d_adens_icol(:) = 0.
     5685 d_adens_acol(:) = 0.
     5686 d_adens_bnd(:) = 0.
     5687
     5688
    32315689      DO i = 1, klon
    32325690        IF (wk_adv(i)) THEN
    … …  
    32525710      DO i = 1, klon
    32535711        IF (wk_adv(i)) THEN
    3254           tau_wk_inv(i) = max( (3.*cstar(i))/(irad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
    3255           tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
     5712!  print *,'ZZZZpopdyn3 wgen(1) ',wgen(1)
     5713!  print *,'ZZZZpopdyn3 cstar(1) ',cstar(1)
     5714!  print *,'ZZZZpopdyn3 isigmaw(1) ',isigmaw(1)
     5715!  print *,'ZZZZpopdyn3 gfl(1) ',gfl(1)
     5716!!          tau_wk_inv(i) = max( (3.*cstar(i))/(irad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
     5717          tau_wk_inv(i) = min(max( (3.*cstar(i))/(irad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.), 1./dtimesub)
    32565718          tau_prime(i) = tau_cv
    32575719
    32585720          d_sig_gen(i) = wgen(i)*aa0
    3259           d_sig_death(i) = - isigmaw(i)*tau_wk_inv_min
     5721          d_sig_death(i) = - isigmaw(i)*tau_wk_inv(i)
    32605722          d_sig_col(i) = - 2.*igfl(i)*cstar(i)*iwdens(i)*(2.*is_wk(i)-aa0)
    32615723          d_sig_spread(i) = gfl(i)*cstar(i)
    … …  
    32665728          d_sig_spread(i) =  d_sig_spread(i)*dtimesub
    32675729          d_sigmaw(i) =  d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
    3268 IF (CPPKEY_IOPHYS_WK) THEN
    3269           IF (phys_sub) call iophys_ecrit('d_sigmaw0',1,'d_sigmaw0','',d_sigmaw)
    3270 END IF
    3271 
    32725730         
    32735731          d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
    … …  
    32775735          d_sigmaw(i) = d_sigmaw_targ
    32785736!!          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
    3279 IF (CPPKEY_IOPHYS_WK) THEN
    3280           IF (phys_sub) THEN
    3281              call iophys_ecrit('tauwk_inv',1,'tau_wk_inv_min','',tau_wk_inv_min)
    3282              call iophys_ecrit('d_sigmaw',1,'d_sigmaw','',d_sigmaw)
    3283              call iophys_ecrit('d_sig_gen',1,'d_sig_gen','',d_sig_gen)
    3284              call iophys_ecrit('d_sig_death',1,'d_sig_death','',d_sig_death)
    3285              call iophys_ecrit('d_sig_col',1,'d_sig_col','',d_sig_col)
    3286              call iophys_ecrit('d_sig_spread',1,'d_sig_spread','',d_sig_spread)
    3287              call iophys_ecrit('d_sig_bnd',1,'d_sig_bnd','',d_sig_bnd)
    3288           ENDIF
    3289 END IF
    32905737          d_asig_death(i) = - asigmaw(i)/tau_prime(i)
    32915738          d_asig_aicol(i) = (agfl(i)*iwdens(i) + igfl(i)*awdens(i))*cstar(i)*is_wk(i)
    … …  
    32985745          d_asig_spread(i) =  d_asig_spread(i)*dtimesub
    32995746          d_asigmaw(i) =  d_sig_gen(i) + d_asig_death(i) + d_asig_aicol(i) + d_asig_iicol(i) + d_asig_spread(i)
    3300 IF (CPPKEY_IOPHYS_WK) THEN
    3301           IF (phys_sub) call iophys_ecrit('d_asigmaw0',1,'d_asigmaw0','',d_asigmaw)
    3302 END IF
    3303 
     5747!
    33045748          d_sigmaw_targ = min(max(d_asigmaw(i),-asigmaw(i)), sigmaw(i)-asigmaw(i))
    33055749!!          d_dens_bnd(i) = d_dens_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
    33065750          d_asig_bnd(i) = d_sigmaw_targ - d_asigmaw(i)
    33075751          d_asigmaw(i) = d_sigmaw_targ
    3308 IF (CPPKEY_IOPHYS_WK) THEN
    3309           IF (phys_sub) THEN
    3310              call iophys_ecrit('d_asigmaw',1,'d_asigmaw','',d_asigmaw)
    3311              call iophys_ecrit('d_asig_death',1,'d_asig_death','',d_asig_death)
    3312              call iophys_ecrit('d_asig_aicol',1,'d_asig_aicol','',d_asig_aicol)
    3313              call iophys_ecrit('d_asig_iicol',1,'d_asig_iicol','',d_asig_iicol)
    3314              call iophys_ecrit('d_asig_spread',1,'d_asig_spread','',d_asig_spread)
    3315              call iophys_ecrit('d_asig_bnd',1,'d_asig_bnd','',d_asig_bnd)
    3316           ENDIF
    3317 END IF
    33185752          d_dens_gen(i) = wgen(i)
    3319           d_dens_death(i) = - iwdens(i)*tau_wk_inv_min
     5753          d_dens_death(i) = - iwdens(i)*tau_wk_inv(i)
    33205754          d_dens_col(i) =  - 2.*gfl(i)*cstar(i)*wdens(i)
    33215755!
    … …  
    33295763          d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
    33305764          d_wdens(i) = d_wdens_targ
    3331 IF (CPPKEY_IOPHYS_WK) THEN
    3332     IF (phys_sub) THEN
    3333         call iophys_ecrit('d_wdens',1,'d_wdens','',d_wdens)
    3334         call iophys_ecrit('d_dens_gen',1,'d_dens_gen','',d_dens_gen)
    3335         call iophys_ecrit('d_dens_death',1,'d_dens_death','',d_dens_death)
    3336         call iophys_ecrit('d_dens_col',1,'d_dens_col','',d_dens_col)
    3337     ENDIF
    3338 END IF
    3339 
    33405765          d_adens_death(i) = -awdens(i)/tau_prime(i)
    33415766          d_adens_icol(i) =   2.*igfl(i)*cstar(i)*iwdens(i)
    … …  
    33465771          d_adens_acol(i)  =   d_adens_acol(i)*dtimesub
    33475772          d_awdens(i) =   d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i)     
    3348 IF (CPPKEY_IOPHYS_WK) THEN
    3349     IF (phys_sub) THEN
    3350         call iophys_ecrit('d_awdens',1,'d_awdens','',d_awdens)
    3351         call iophys_ecrit('d_adens_death',1,'d_adens_death','',d_adens_death)
    3352         call iophys_ecrit('d_adens_icol',1,'d_adens_icol','',d_adens_icol)
    3353         call iophys_ecrit('d_adens_acol',1,'d_adens_acol','',d_adens_acol)
    3354     ENDIF
    3355 END IF
    33565773          d_wdens_targ = min(max(d_awdens(i),-awdens(i)), wdens(i)-awdens(i))
    33575774!!          d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
    … …  
    33705787        ENDIF
    33715788      ENDDO
     5789IF (CPPKEY_IOPHYS_WK) THEN
     5790    IF (phys_sub) THEN
     5791       call iophys_ecrit('d_sigmaw0',1,'d_sigmaw0','',d_sigmaw)
     5792!
     5793       call iophys_ecrit('cstar',1,'cstar','',cstar)
     5794       call iophys_ecrit('wgen_pd3',1,'wgen_popdyn3','',wgen)
     5795       call iophys_ecrit('tauwk_inv',1,'tau_wk_inv','',tau_wk_inv)
     5796       call iophys_ecrit('d_sigmaw',1,'d_sigmaw','',d_sigmaw)
     5797       call iophys_ecrit('d_sig_gen',1,'d_sig_gen','',d_sig_gen)
     5798       call iophys_ecrit('d_sig_death',1,'d_sig_death','',d_sig_death)
     5799       call iophys_ecrit('d_sig_col',1,'d_sig_col','',d_sig_col)
     5800       call iophys_ecrit('d_sig_spread',1,'d_sig_spread','',d_sig_spread)
     5801       call iophys_ecrit('d_sig_bnd',1,'d_sig_bnd','',d_sig_bnd)
     5802!
     5803       call iophys_ecrit('d_asigmaw0',1,'d_asigmaw0','',d_asigmaw)
     5804!
     5805       call iophys_ecrit('d_asigmaw',1,'d_asigmaw','',d_asigmaw)
     5806       call iophys_ecrit('d_asig_death',1,'d_asig_death','',d_asig_death)
     5807       call iophys_ecrit('d_asig_aicol',1,'d_asig_aicol','',d_asig_aicol)
     5808       call iophys_ecrit('d_asig_iicol',1,'d_asig_iicol','',d_asig_iicol)
     5809       call iophys_ecrit('d_asig_spread',1,'d_asig_spread','',d_asig_spread)
     5810       call iophys_ecrit('d_asig_bnd',1,'d_asig_bnd','',d_asig_bnd)
     5811!
     5812       call iophys_ecrit('d_wdens',1,'d_wdens','',d_wdens)
     5813       call iophys_ecrit('d_dens_gen',1,'d_dens_gen','',d_dens_gen)
     5814       call iophys_ecrit('d_dens_death',1,'d_dens_death','',d_dens_death)
     5815       call iophys_ecrit('d_dens_col',1,'d_dens_col','',d_dens_col)
     5816! 
     5817       call iophys_ecrit('d_awdens',1,'d_awdens','',d_awdens)
     5818       call iophys_ecrit('d_adens_death',1,'d_adens_death','',d_adens_death)
     5819       call iophys_ecrit('d_adens_icol',1,'d_adens_icol','',d_adens_icol)
     5820       call iophys_ecrit('d_adens_acol',1,'d_adens_acol','',d_adens_acol)
     5821    ENDIF
     5822END IF
     5823
    33725824
    33735825      IF (prt_level >= 10) THEN
    … …  
    33865838    RETURN
    33875839    END SUBROUTINE wake_popdyn_3 
     5840   
     5841    SUBROUTINE wake_dadv(klon, klev, dtime, ph, ppi, wk_adv, kupper,  &
     5842                         deltomg, dp_deltomg, sigmaw, dsigspread,  &
     5843                         thw, thx, qw, qx, &
     5844                         d_deltat_dadv, d_deltaq_dadv, d_tb_dadv, d_qb_dadv)
     5845
     5846  USE lmdz_wake_ini , ONLY : flag_dadv_implicit
     5847
     5848IMPLICIT NONE
     5849
     5850  INTEGER, INTENT(IN)                                     :: klon, klev
     5851  REAL,                               INTENT(IN)          :: dtime
     5852  REAL, DIMENSION (klon, klev+1),     INTENT(IN)          :: ph
     5853  REAL, DIMENSION (klon, klev),       INTENT(IN)          :: ppi
     5854  LOGICAL, DIMENSION (klon),          INTENT(IN)          :: wk_adv
     5855  INTEGER, DIMENSION (klon),          INTENT(IN)          :: kupper
     5856  REAL, DIMENSION (klon, klev),       INTENT(IN)          :: deltomg
     5857  REAL, DIMENSION (klon, klev),       INTENT(IN)          :: dp_deltomg
     5858  REAL, DIMENSION (klon),             INTENT(IN)          :: sigmaw
     5859  REAL, DIMENSION (klon),             INTENT(IN)          :: dsigspread
     5860  REAL, DIMENSION (klon, klev),       INTENT(IN)          :: thw           ! component # 1
     5861  REAL, DIMENSION (klon, klev),       INTENT(IN)          :: thx           ! component # 2
     5862  REAL, DIMENSION (klon, klev),       INTENT(IN)          :: qw            ! component # 1
     5863  REAL, DIMENSION (klon, klev),       INTENT(IN)          :: qx            ! component # 2
     5864
     5865  REAL, DIMENSION (klon, klev),       INTENT(OUT)         :: d_deltat_dadv
     5866  REAL, DIMENSION (klon, klev),       INTENT(OUT)         :: d_deltaq_dadv
     5867  REAL, DIMENSION (klon, klev),       INTENT(OUT)         :: d_tb_dadv
     5868  REAL, DIMENSION (klon, klev),       INTENT(OUT)         :: d_qb_dadv
     5869
     5870! Internal variables
     5871  INTEGER                               :: i, k
     5872  REAL, DIMENSION (klon, klev)          :: entr_s    ! entrainment into wakes due to spread
     5873  REAL, DIMENSION (klon, klev)          :: thb, qb
     5874  REAL, DIMENSION (klon, klev)          :: delta_th, delta_q
     5875
     5876! Tests
     5877
     5878! Arrays used in the implicit scheme
     5879  REAL, DIMENSION (klon)                :: rr11, rr12, rr21, rr22
     5880
     5881  REAL, DIMENSION (klon, klev)          :: aa11, aa12, aa21, aa22
     5882  REAL, DIMENSION (klon, klev)          :: bb11, bb12, bb21, bb22
     5883  REAL, DIMENSION (klon, klev)          :: cc11, cc12, cc21, cc22
     5884
     5885  REAL, DIMENSION (klon, klev)          :: alpha11, alpha12, alpha21, alpha22
     5886  REAL, DIMENSION (klon, klev)          :: beta11, beta12, beta21, beta22
     5887  REAL, DIMENSION (klon, klev)          :: gamma11, gamma12, gamma21, gamma22
     5888  REAL, DIMENSION (klon, klev)          :: ai11, ai12, ai21, ai22             ! inverse of alpha
     5889
     5890  REAL, DIMENSION (klon, klev)          :: xt1, xt2, xq1, xq2
     5891  REAL, DIMENSION (klon, klev)          :: yt1, yt2, yq1, yq2
     5892  REAL, DIMENSION (klon, klev)          :: zt1, zt2, zq1, zq2
     5893  REAL, DIMENSION (klon, klev)          :: th1, th2, q1, q2
     5894
     5895  REAL                                  :: coef, det
     5896
     5897  REAL, DIMENSION (klon,klev)           :: xt1inv, xt2inv, xq1inv, xq2inv
     5898
     5899! Arrays used in the explicit scheme (vertical gradients)
     5900  REAL, DIMENSION (klon, klev)          :: d_thx, d_qx
     5901  REAL, DIMENSION (klon, klev)          :: d_thw, d_qw
     5902  REAL, DIMENSION (klon, klev)          :: d_dth, d_dq
     5903
     5904!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     5905
     5906    entr_s(:,:) = 0.
     5907    delta_th(:,:) = 0.
     5908   
     5909    d_deltat_dadv(:,:) = 0.
     5910    d_deltaq_dadv(:,:) = 0.
     5911    d_tb_dadv(:,:) = 0.
     5912    d_qb_dadv(:,:) = 0.
     5913
     5914
     5915    rr11(:) = sigmaw(:)
     5916    rr12(:) = 1.-sigmaw(:)
     5917    rr21(:) = 1.
     5918    rr22(:) = -1.
     5919
     5920    DO k = 1, klev
     5921      DO i = 1,klon
     5922        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
     5923         thb(i,k)      = rr11(i)*thw(i,k)+rr12(i)*thx(i,k)
     5924         delta_th(i,k) = rr21(i)*thw(i,k)+rr22(i)*thx(i,k)
     5925       
     5926         qb(i,k)      = rr11(i)*qw(i,k)+rr12(i)*qx(i,k)
     5927         delta_q(i,k) = rr21(i)*qw(i,k)+rr22(i)*qx(i,k)
     5928        ENDIF
     5929      ENDDO
     5930    ENDDO
     5931
     5932    DO i = 1, klon
     5933        entr_s(i,klev) = 0.
     5934    ENDDO
     5935
     5936    DO k = 1, klev-1
     5937      DO i = 1,klon
     5938        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
     5939!!        entr_s(i,k) = dsigspread(i) - sigmaw(i)*(1.-sigmaw(i))*(deltomg(i,k+1)-deltomg(i,k)) / &
     5940!!                     (ph(i,k)-ph(i,k+1))   
     5941
     5942          entr_s(i,k) = dsigspread(i) + sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i,k)
     5943!!  print *,'dadv, k, dp_deltomg(i,k), (deltomg(i,k)-deltomg(i,k+1))/(ph(i,k)-ph(i,k+1)) ', &
     5944!!                 k, dp_deltomg(i,k), (deltomg(i,k)-deltomg(i,k+1))/(ph(i,k)-ph(i,k+1))
     5945
     5946        ENDIF
     5947      ENDDO
     5948    ENDDO
     5949
     5950! -------------------------------------------------------------------------------------------
     5951!   Depending on flag_dadv_implicit, use implicit upstream scheme or explicit upstream scheme
     5952! -------------------------------------------------------------------------------------------
     5953
     5954  IF (flag_dadv_implicit) THEN
     5955
     5956!   Implicit scheme : solve for d_deltat_dadv and d_tb_dadv
     5957!                     (and similarly for d_deltaq_dadv and d_qb_dadv).
     5958!                     The system to be inverted is block-tridiagonal with 2x2 blocks.
     5959! -----------------------------------------------------------------------------------------
     5960
     5961!        Matrix indexing:                   Theta_w     Theta_x
     5962!
     5963!                                           /               
     5964!                           Theta_b        |  A11        A12 |
     5965!                                          |                 |
     5966!                           delta_theta    |  A21        A22 |
     5967!                                                           /
     5968!       Tridiagonal matrix
     5969!         /                                                         
     5970!         |   aa(1)   bb(1)  0                                       |
     5971!         |   cc(2)   aa(2)  bb(2)  0                                |
     5972!         |   0       cc(3)  aa(3)  bb(3)                            |
     5973!         |                                                          |
     5974!                     .      .      .       .                         
     5975!                                                                     
     5976!                            .      .       .       .                 
     5977!         |                                                          |
     5978!         |                         cc(n-2) aa(n-2) bb(n-2) 0        |   
     5979!         |                         0       cc(n-1) aa(n-1) bb(n-1)  |             
     5980!                                           0       cc(n)   aa(n)    /             
     5981! -----------------------------------------------------------------------------------------
     5982
     5983!! Building the tridiagonal matrix
     5984    DO i = 1,klon
     5985      IF (wk_adv(i)) THEN
     5986        k = kupper(i)
     5987        coef = dtime/(ph(i,k)-ph(i,k+1))
     5988        aa11(i,k) = rr11(i)+coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k)
     5989        aa12(i,k) = rr12(i)
     5990        aa21(i,k) = rr21(i)+coef*( dsigspread(i)/sigmaw(i)*(ph(i,k)-ph(i,k+1)) + &
     5991                                     (1.-sigmaw(i))*deltomg(i,k) )
     5992        aa22(i,k) = rr22(i)+coef*(-dsigspread(i)/sigmaw(i)*(ph(i,k)-ph(i,k+1)) - &
     5993                                     deltomg(i,k) )
     5994   
     5995        cc11(i,k) = 0.
     5996        cc12(i,k) = -coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k)
     5997        cc21(i,k) = 0.
     5998        cc22(i,k) = coef*sigmaw(i)*deltomg(i,k)
     5999      ENDIF  ! (wk_adv(i))
     6000    ENDDO
     6001    DO k = 2, klev-1
     6002      DO i = 1,klon
     6003        IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
     6004          coef = dtime/(ph(i,k)-ph(i,k+1))
     6005          aa11(i,k) = rr11(i)+coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k)
     6006          aa12(i,k) = rr12(i)+coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k+1)
     6007          aa21(i,k) = rr21(i)+coef*( dsigspread(i)/sigmaw(i)*(ph(i,k)-ph(i,k+1)) + &
     6008                                     (1.-sigmaw(i))*deltomg(i,k) )
     6009          aa22(i,k) = rr22(i)+coef*(-dsigspread(i)/sigmaw(i)*(ph(i,k)-ph(i,k+1)) + &
     6010                                     (1.-sigmaw(i))*(deltomg(i,k+1)-deltomg(i,k)) - &
     6011                                     sigmaw(i)*deltomg(i,k) )
     6012   
     6013          bb11(i,k) =  -coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k+1)
     6014          bb12(i,k) =  0.
     6015          bb21(i,k) =  -coef*(1.-sigmaw(i))*deltomg(i,k+1)
     6016          bb22(i,k) =  0.
     6017   
     6018          cc11(i,k) = 0.
     6019          cc12(i,k) = -coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,k)
     6020          cc21(i,k) = 0.
     6021          cc22(i,k) = coef*sigmaw(i)*deltomg(i,k)
     6022        ENDIF  ! (wk_adv(i) .AND. k<=kupper(i))
     6023      ENDDO
     6024    ENDDO
     6025    DO i = 1,klon
     6026      IF (wk_adv(i)) THEN
     6027        coef = dtime/(ph(i,1)-ph(i,2))
     6028        aa11(i,1) = rr11(i)+coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,1)
     6029        aa12(i,1) = rr12(i)+coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,2)
     6030        aa21(i,1) = rr21(i)+coef*( dsigspread(i)/sigmaw(i)*(ph(i,1)-ph(i,2)) + &
     6031                                   (1.-sigmaw(i))*deltomg(i,1) )
     6032        aa22(i,1) = rr22(i)+coef*(-dsigspread(i)/sigmaw(i)*(ph(i,1)-ph(i,2)) + &
     6033                                   (1.-sigmaw(i))*(deltomg(i,2)-deltomg(i,1)) - &
     6034                                   sigmaw(i)*deltomg(i,1) )
     6035   
     6036        bb11(i,1) =  -coef*sigmaw(i)*(1.-sigmaw(i))*deltomg(i,2)
     6037        bb12(i,1) =  0.
     6038        bb21(i,1) =  -coef*(1.-sigmaw(i))*deltomg(i,2)
     6039        bb22(i,1) =  0.
     6040      ENDIF  ! (wk_adv(i))
     6041    ENDDO
     6042
     6043!!  printing the tridiagonal matrix
     6044!!!  First row
     6045!!   k = 1
     6046!!   print 1789, k, aa11(1,1), aa12(1,1), bb11(1,1), bb12(1,1)
     6047!!   print 1789, k, aa21(1,1), aa22(1,1), bb21(1,1), bb22(1,1)
     6048!!1789 FORMAT(1X, I3, 3(4X, 2E13.5))
     6049!!        coef = dtime/(ph(1,k)-ph(1,k+1))
     6050!!   print *,'rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,1), deltomg(1,2) ', &
     6051!!            rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,1), deltomg(1,2)
     6052!!
     6053!!!  Rows 2 to klev-1
     6054!!   DO k = 2, klev-1
     6055!!     print 1789, k, cc11(1,k), cc12(1,k), aa11(1,k), aa12(1,k), bb11(1,k), bb12(1,k)
     6056!!     print 1789, k, cc21(1,k), cc22(1,k), aa21(1,k), aa22(1,k), bb21(1,k), bb22(1,k)
     6057!!        coef = dtime/(ph(1,k)-ph(1,k+1))
     6058!!   print *,'rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,k), deltomg(1,k+1) ', &
     6059!!            rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,k), deltomg(1,k+1)
     6060!!   ENDDO
     6061!!
     6062!!!  Row klev
     6063!!     print 1789, klev, cc11(1,klev), cc12(1,klev), aa11(1,klev), aa12(1,klev)
     6064!!     print 1789, klev, cc21(1,klev), cc22(1,klev), aa21(1,klev), aa22(1,klev)
     6065!!        coef = dtime/(ph(1,klev)-ph(1,klev+1))
     6066!!   print *,'rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,klev) ', &
     6067!!            rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,klev)
     6068
     6069
     6070!! Downward loop
     6071
     6072   xt1(:,:) = thb(:,:)     
     6073   xt2(:,:) = delta_th(:,:)
     6074   xq1(:,:) = qb(:,:)     
     6075   xq2(:,:) = delta_q(:,:)
     6076
     6077    DO i = 1,klon
     6078      IF (wk_adv(i)) THEN
     6079        k = kupper(i)
     6080        alpha11(:,k)=aa11(:,k)
     6081        alpha12(:,k)=aa12(:,k)
     6082        alpha21(:,k)=aa21(:,k)
     6083        alpha22(:,k)=aa22(:,k)
     6084        beta11(:,k)=0.
     6085        beta12(:,k)=0.
     6086        beta21(:,k)=0.
     6087        beta22(:,k)=0.
     6088        yt1(i,k) = xt1(i,k)
     6089        yt2(i,k) = xt2(i,k)
     6090        yq1(i,k) = xq1(i,k)
     6091        yq2(i,k) = xq2(i,k)
     6092      ENDIF  ! (wk_adv(i))
     6093    ENDDO
     6094    DO i = 1,klon
     6095      IF (wk_adv(i)) THEN
     6096        k = kupper(i)
     6097        det=alpha11(i,k)*alpha22(i,k) - alpha12(i,k)*alpha21(i,k)
     6098        ai11(i,k)= alpha22(i,k)/det
     6099        ai12(i,k)=-alpha12(i,k)/det
     6100        ai21(i,k)=-alpha21(i,k)/det
     6101        ai22(i,k)= alpha11(i,k)/det
     6102        zt1(i,k) = ai11(i,k)*yt1(i,k) + ai12(i,k)*yt2(i,k)
     6103        zt2(i,k) = ai21(i,k)*yt1(i,k) + ai22(i,k)*yt2(i,k)
     6104        zq1(i,k) = ai11(i,k)*yq1(i,k) + ai12(i,k)*yq2(i,k)
     6105        zq2(i,k) = ai21(i,k)*yq1(i,k) + ai22(i,k)*yq2(i,k)
     6106      ENDIF  ! (wk_adv(i))
     6107    ENDDO
     6108
     6109    DO k = klev, 2, -1
     6110      DO i = 1,klon
     6111        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     6112          gamma11(i,k) = ai11(i,k)*cc11(i,k) + ai12(i,k)*cc21(i,k)
     6113          gamma12(i,k) = ai11(i,k)*cc12(i,k) + ai12(i,k)*cc22(i,k)
     6114          gamma21(i,k) = ai21(i,k)*cc11(i,k) + ai22(i,k)*cc21(i,k)
     6115          gamma22(i,k) = ai21(i,k)*cc12(i,k) + ai22(i,k)*cc22(i,k)
     6116   
     6117          alpha11(i,k-1) = aa11(i,k-1) - ( bb11(i,k-1)*gamma11(i,k)+bb12(i,k-1)*gamma21(i,k) )
     6118          alpha12(i,k-1) = aa12(i,k-1) - ( bb11(i,k-1)*gamma12(i,k)+bb12(i,k-1)*gamma22(i,k) )
     6119          alpha21(i,k-1) = aa21(i,k-1) - ( bb21(i,k-1)*gamma11(i,k)+bb22(i,k-1)*gamma21(i,k) )
     6120          alpha22(i,k-1) = aa22(i,k-1) - ( bb21(i,k-1)*gamma12(i,k)+bb22(i,k-1)*gamma22(i,k) )
     6121   
     6122          beta11(i,k-1) = bb11(i,k-1)*ai11(i,k)+bb12(i,k-1)*ai21(i,k)
     6123          beta12(i,k-1) = bb11(i,k-1)*ai12(i,k)+bb12(i,k-1)*ai22(i,k)
     6124          beta21(i,k-1) = bb21(i,k-1)*ai11(i,k)+bb22(i,k-1)*ai21(i,k)
     6125          beta22(i,k-1) = bb21(i,k-1)*ai12(i,k)+bb22(i,k-1)*ai22(i,k)
     6126   
     6127          yt1(i,k-1) = xt1(i,k-1) - ( beta11(i,k-1)*yt1(i,k) +beta12(i,k-1)*yt2(i,k) )
     6128          yt2(i,k-1) = xt2(i,k-1) - ( beta21(i,k-1)*yt1(i,k) +beta22(i,k-1)*yt2(i,k) )
     6129          yq1(i,k-1) = xq1(i,k-1) - ( beta11(i,k-1)*yq1(i,k) +beta12(i,k-1)*yq2(i,k) )
     6130          yq2(i,k-1) = xq2(i,k-1) - ( beta21(i,k-1)*yq1(i,k) +beta22(i,k-1)*yq2(i,k) )
     6131   
     6132          det=alpha11(i,k-1)*alpha22(i,k-1) - alpha12(i,k-1)*alpha21(i,k-1)
     6133          ai11(i,k-1)= alpha22(i,k-1)/det
     6134          ai12(i,k-1)=-alpha12(i,k-1)/det
     6135          ai21(i,k-1)=-alpha21(i,k-1)/det
     6136          ai22(i,k-1)= alpha11(i,k-1)/det
     6137   
     6138          zt1(i,k-1) = ai11(i,k-1)*yt1(i,k-1)+ai12(i,k-1)*yt2(i,k-1)
     6139          zt2(i,k-1) = ai21(i,k-1)*yt1(i,k-1)+ai22(i,k-1)*yt2(i,k-1)
     6140          zq1(i,k-1) = ai11(i,k-1)*yq1(i,k-1)+ai12(i,k-1)*yq2(i,k-1)
     6141          zq2(i,k-1) = ai21(i,k-1)*yq1(i,k-1)+ai22(i,k-1)*yq2(i,k-1)
     6142        ENDIF  ! (wk_adv(i) .AND. k<=kupper(i))
     6143      ENDDO
     6144    ENDDO
     6145       
     6146!! Upward loop
     6147
     6148    DO i = 1,klon
     6149      IF (wk_adv(i)) THEN
     6150       th1(i,1) = zt1(i,1)
     6151       th2(i,1) = zt2(i,1)
     6152       q1(i,1)  = zq1(i,1)
     6153       q2(i,1)  = zq2(i,1)
     6154     
     6155       d_tb_dadv(i,1) =     ( rr11(i)*(th1(i,1)-thw(i,1))+rr12(i)*(th2(i,1)-thx(i,1)) )*ppi(i,1)
     6156       d_deltat_dadv(i,1) = ( rr21(i)*(th1(i,1)-thw(i,1))+rr22(i)*(th2(i,1)-thx(i,1)) )*ppi(i,1)
     6157       d_qb_dadv(i,1) =       rr11(i)*(q1(i,1) -qw(i,1)) +rr12(i)*(q2(i,1)-qx(i,1))
     6158       d_deltaq_dadv(i,1) =   rr21(i)*(q1(i,1) -qw(i,1)) +rr22(i)*(q2(i,1)-qx(i,1))
     6159      ENDIF  ! (wk_adv(i))
     6160    ENDDO
     6161
     6162    DO k = 2, klev
     6163      DO i = 1,klon
     6164        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
     6165          th1(i,k) = zt1(i,k) - ( gamma11(i,k)*th1(i,k-1)+gamma12(i,k)*th2(i,k-1) )
     6166          th2(i,k) = zt2(i,k) - ( gamma21(i,k)*th1(i,k-1)+gamma22(i,k)*th2(i,k-1) )
     6167          q1(i,k)  = zq1(i,k) - ( gamma11(i,k)*q1(i,k-1) +gamma12(i,k)*q2(i,k-1) )
     6168          q2(i,k)  = zq2(i,k) - ( gamma21(i,k)*q1(i,k-1) +gamma22(i,k)*q2(i,k-1) )
     6169   
     6170          d_tb_dadv(i,k) =     ( rr11(i)*(th1(i,k)-thw(i,k))+rr12(i)*(th2(i,k)-thx(i,k)) )*ppi(i,k)
     6171          d_deltat_dadv(i,k) = ( rr21(i)*(th1(i,k)-thw(i,k))+rr22(i)*(th2(i,k)-thx(i,k)) )*ppi(i,k)
     6172          d_qb_dadv(i,k) =       rr11(i)*(q1(i,k)-qw(i,k))  +rr12(i)*(q2(i,k)-qx(i,k))
     6173          d_deltaq_dadv(i,k) =   rr21(i)*(q1(i,k)-qw(i,k))  +rr22(i)*(q2(i,k)-qx(i,k))
     6174        ENDIF  ! (wk_adv(i) .AND. k<=kupper(i))
     6175      ENDDO
     6176    ENDDO
     6177
     6178!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     6179!!!!!!!!!!!!!!!!!!       Verification de l'inversion                !!!!!!!!!!!!!!!!!!!!!!!
     6180!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     6181!!
     6182!!    DO i = 1,klon
     6183!!        xt1inv(i,1) = aa11(i,1)*th1(i,1) + aa12(i,1)*th2(i,1) + bb11(i,1)*th1(i,2) + bb12(i,1)*th2(i,2)
     6184!!        xt2inv(i,1) = aa21(i,1)*th1(i,1) + aa22(i,1)*th2(i,1) + bb21(i,1)*th1(i,2) + bb22(i,1)*th2(i,2)
     6185!!        xq1inv(i,1)  = aa11(i,1)*q1(i,1)  + aa12(i,1)*q2(i,1)  + bb11(i,1)*q1(i,2)  + bb12(i,1)*q2(i,2)
     6186!!        xq2inv(i,1)  = aa21(i,1)*q1(i,1)  + aa22(i,1)*q2(i,1)  + bb21(i,1)*q1(i,2)  + bb22(i,1)*q2(i,2)
     6187!!    ENDDO
     6188!!   
     6189!!      DO k = 2, klev-1
     6190!!        DO i = 1,klon
     6191!!        xt1inv(i,k) = aa11(i,k)*th1(i,k) + aa12(i,k)*th2(i,k) + bb11(i,k)*th1(i,k+1) + bb12(i,k)*th2(i,k+1) &
     6192!!                                                              + cc11(i,k)*th1(i,k-1) + cc12(i,k)*th2(i,k-1)
     6193!!        xt2inv(i,k) = aa21(i,k)*th1(i,k) + aa22(i,k)*th2(i,k) + bb21(i,k)*th1(i,k+1) + bb22(i,k)*th2(i,k+1) &
     6194!!                                                              + cc21(i,k)*th1(i,k-1) + cc22(i,k)*th2(i,k-1)
     6195!!        xq1inv(i,k)  = aa11(i,k)*q1(i,k)  + aa12(i,k)*q2(i,k)  + bb11(i,k)*q1(i,k+1)  + bb12(i,k)*q2(i,k+1)  &
     6196!!                                                              + cc11(i,k)*q1(i,k-1)  + cc12(i,k)*q2(i,k-1)
     6197!!        xq2inv(i,k)  = aa21(i,k)*q1(i,k)  + aa22(i,k)*q2(i,k)  + bb21(i,k)*q1(i,k+1)  + bb22(i,k)*q2(i,k+1)  &
     6198!!                                                              + cc21(i,k)*q1(i,k-1)  + cc22(i,k)*q2(i,k-1)
     6199!!        ENDDO
     6200!!      ENDDO
     6201!!   
     6202!!    DO i = 1,klon
     6203!!        xt1inv(i,klev) = aa11(i,klev)*th1(i,klev) + aa12(i,klev)*th2(i,klev) + cc11(i,klev)*th1(i,klev-1) + cc12(i,klev)*th2(i,klev-1)
     6204!!        xt2inv(i,klev) = aa21(i,klev)*th1(i,klev) + aa22(i,klev)*th2(i,klev) + cc21(i,klev)*th1(i,klev-1) + cc22(i,klev)*th2(i,klev-1)
     6205!!        xq1inv(i,klev)  = aa11(i,klev)*q1(i,klev)  + aa12(i,klev)*q2(i,klev)  + cc11(i,klev)*q1(i,klev-1)  + cc12(i,klev)*q2(i,klev-1)
     6206!!        xq2inv(i,klev)  = aa21(i,klev)*q1(i,klev)  + aa22(i,klev)*q2(i,klev)  + cc21(i,klev)*q1(i,klev-1)  + cc22(i,klev)*q2(i,klev-1)
     6207!!    ENDDO
     6208!!   
     6209!!    DO k = 1, 20
     6210!!      IF (abs(xt1inv(1,k)-xt1(1,k)) .GT. 1.e-15*xt1(1,k) ) THEN
     6211!!        print *,'wake_dadv, k, xt1inv(1,k), xt1(1,k), xt1inv(1,k)-xt1(1,k) ', &
     6212!!                            k, xt1inv(1,k), xt1(1,k), xt1inv(1,k)-xt1(1,k)
     6213!!      ENDIF
     6214!!      IF (abs(xt2inv(1,k)-xt2(1,k)) .GT. 1.e-15*xt2(1,k) ) THEN
     6215!!        print *,'wake_dadv, k, xt2inv(1,k), xt2(1,k), xt2inv(1,k)-xt2(1,k) ', &
     6216!!                            k, xt2inv(1,k), xt2(1,k), xt2inv(1,k)-xt2(1,k)
     6217!!      ENDIF
     6218!!      IF (abs(xq1inv(1,k)-xq1(1,k)) .GT. 1.e-15*xq1(1,k) ) THEN
     6219!!        print *,'wake_dadv, k, xq1inv(1,k), xq1(1,k), xq1inv(1,k)-xq1(1,k) ', &
     6220!!                            k, xq1inv(1,k), xq1(1,k), xq1inv(1,k)-xq1(1,k)
     6221!!      ENDIF
     6222!!      IF (abs(xq2inv(1,k)-xq2(1,k)) .GT. 1.e-15*xq2(1,k) ) THEN
     6223!!        print *,'wake_dadv, k, xq2inv(1,k), xq2(1,k), xq2inv(1,k)-xq2(1,k) ', &
     6224!!                          k, xq2inv(1,k), xq2(1,k), xq2inv(1,k)-xq2(1,k)
     6225!!      ENDIF
     6226!!    ENDDO
     6227
     6228  ELSE  ! (flag_dadv_implicit)
     6229
     6230!   Explicit scheme : compute directly d_deltat_dadv and d_tb_dadv
     6231!                     (and similarly for d_deltaq_dadv and d_qb_dadv).
     6232! -----------------------------------------------------------------------------------------
     6233
     6234    DO i = 1, klon
     6235      IF (wk_adv(i)) THEN !!! nrlmd
     6236        d_thx(i, 1) = 0.
     6237        d_thw(i, 1) = 0.
     6238        d_dth(i, 1) = 0.
     6239        d_qx(i, 1) = 0.
     6240        d_qw(i, 1) = 0.
     6241        d_dq(i, 1) = 0.
     6242      END IF
     6243    END DO
     6244
     6245    DO k = 2, klev
     6246      DO i = 1, klon
     6247        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
     6248          d_thx(i, k) = thx(i, k-1) - thx(i, k)
     6249          d_thw(i, k) = thw(i, k-1) - thw(i, k)
     6250          d_dth(i, k) = delta_th(i, k-1) - delta_th(i, k)
     6251          d_qx(i, k) = qx(i, k-1) - qx(i, k)
     6252          d_qw(i, k) = qw(i, k-1) - qw(i, k)
     6253          d_dq(i, k) = delta_q(i, k-1) - delta_q(i, k)
     6254        END IF ! (wk_adv(i) .AND. k<=kupper(i)+1)
     6255      END DO
     6256    END DO
     6257
     6258    DO k = 1, klev-1
     6259      DO i = 1, klon
     6260        IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
     6261          d_deltat_dadv(i, k) = dtime/(ph(i,k)-ph(i,k+1))* &
     6262            (rr22(i)*deltomg(i,k)*sigmaw(i)*d_thx(i,k) - &
     6263             rr21(i)*deltomg(i,k+1)*(1.-sigmaw(i))*d_thw(i,k+1) )*ppi(i, k) - &
     6264             dtime*entr_s(i,k)*delta_th(i,k)/sigmaw(i)*ppi(i, k)
     6265!
     6266          d_deltaq_dadv(i, k) = dtime/(ph(i,k)-ph(i,k+1))* &
     6267            (rr22(i)*deltomg(i,k)*sigmaw(i)*d_qx(i,k)- &
     6268             rr21(i)*deltomg(i,k+1)*(1.-sigmaw(i))*d_qw(i,k+1) ) - &
     6269             dtime*entr_s(i,k)*delta_q(i,k)/sigmaw(i)
     6270
     6271          ! and increment large scale tendencies
     6272          d_tb_dadv(i, k) = dtime*((rr12(i)*deltomg(i,k)*sigmaw(i)*d_thx(i,k)- &
     6273                                  rr11(i)*deltomg(i,k+1)*(1.-sigmaw(i))*d_thw(i,k+1))/ &
     6274                                 (ph(i,k)-ph(i,k+1)) &
     6275                                 -sigmaw(i)*(1.-sigmaw(i))*delta_th(i,k)*(deltomg(i,k)-deltomg(i,k+1))/ &
     6276                                 (ph(i,k)-ph(i,k+1)) )*ppi(i, k)
     6277
     6278          d_qb_dadv(i, k) = dtime*((rr12(i)*deltomg(i,k)*sigmaw(i)*d_qx(i,k)- &
     6279                                  rr11(i)*deltomg(i,k+1)*(1.-sigmaw(i))*d_qw(i,k+1))/ &
     6280                                 (ph(i,k)-ph(i,k+1)) &
     6281                                 -sigmaw(i)*(1.-sigmaw(i))*delta_q(i,k)*(deltomg(i,k)-deltomg(i,k+1))/ &
     6282                                 (ph(i,k)-ph(i,k+1)) )
     6283        ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
     6284          d_tb_dadv(i, k) = dtime*(rr12(i)*deltomg(i, k)*sigmaw(i)*d_thx(i, k)/(ph(i, k)-ph(i, k+1)))*ppi(i, k)
     6285
     6286          d_qb_dadv(i, k) = dtime*(rr12(i)*deltomg(i, k)*sigmaw(i)*d_qx(i, k)/(ph(i, k)-ph(i, k+1)))
     6287        END IF ! (wk_adv(i) .AND. k<=kupper(i)-1)
     6288      END DO
     6289    END DO
     6290
     6291  ENDIF! (flag_dadv_implicit)
     6292
     6293    END SUBROUTINE wake_dadv
    33886294
    33896295END MODULE lmdz_wake
Note: See TracChangeset for help on using the changeset viewer.