LMDZ
wake.F90
Go to the documentation of this file.
1 
2 ! $Id: wake.F90 2346 2015-08-21 15:13:46Z emillour $
3 
4 SUBROUTINE wake(p, ph, pi, dtime, sigd_con, &
5  te0, qe0, omgb, &
6  dtdwn, dqdwn, amdwn, amup, dta, dqa, &
7  wdtpbl, wdqpbl, udtpbl, udqpbl, &
8  deltatw, deltaqw, dth, hw, sigmaw, wape, fip, gfl, &
9  dtls, dqls, ktopw, omgbdth, dp_omgb, wdens, tu, qu, &
10  dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, spread, cstar, &
11  d_deltat_gw, d_deltatw2, d_deltaqw2)
12 
13 
14  ! **************************************************************
15  ! *
16  ! WAKE *
17  ! retour a un Pupper fixe *
18  ! *
19  ! written by : GRANDPEIX Jean-Yves 09/03/2000 *
20  ! modified by : ROEHRIG Romain 01/29/2007 *
21  ! **************************************************************
22 
23  USE dimphy
24  use mod_phys_lmdz_para
25  USE print_control_mod, ONLY: prt_level
26  IMPLICIT NONE
27  ! ============================================================================
28 
29 
30  ! But : Decrire le comportement des poches froides apparaissant dans les
31  ! grands systemes convectifs, et fournir l'energie disponible pour
32  ! le declenchement de nouvelles colonnes convectives.
33 
34  ! Variables d'etat : deltatw : ecart de temperature wake-undisturbed
35  ! area
36  ! deltaqw : ecart d'humidite wake-undisturbed area
37  ! sigmaw : fraction d'aire occupee par la poche.
38 
39  ! Variable de sortie :
40 
41  ! wape : WAke Potential Energy
42  ! fip : Front Incident Power (W/m2) - ALP
43  ! gfl : Gust Front Length per unit area (m-1)
44  ! dtls : large scale temperature tendency due to wake
45  ! dqls : large scale humidity tendency due to wake
46  ! hw : hauteur de la poche
47  ! dp_omgb : vertical gradient of large scale omega
48  ! wdens : densite de poches
49  ! omgbdth: flux of Delta_Theta transported by LS omega
50  ! dtKE : differential heating (wake - unpertubed)
51  ! dqKE : differential moistening (wake - unpertubed)
52  ! omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
53  ! dp_deltomg : vertical gradient of omg (s-1)
54  ! spread : spreading term in dt_wake and dq_wake
55  ! deltatw : updated temperature difference (T_w-T_u).
56  ! deltaqw : updated humidity difference (q_w-q_u).
57  ! sigmaw : updated wake fractional area.
58  ! d_deltat_gw : delta T tendency due to GW
59 
60  ! Variables d'entree :
61 
62  ! aire : aire de la maille
63  ! te0 : temperature dans l'environnement (K)
64  ! qe0 : humidite dans l'environnement (kg/kg)
65  ! omgb : vitesse verticale moyenne sur la maille (Pa/s)
66  ! dtdwn: source de chaleur due aux descentes (K/s)
67  ! dqdwn: source d'humidite due aux descentes (kg/kg/s)
68  ! dta : source de chaleur due courants satures et detrain (K/s)
69  ! dqa : source d'humidite due aux courants satures et detra (kg/kg/s)
70  ! amdwn: flux de masse total des descentes, par unite de
71  ! surface de la maille (kg/m2/s)
72  ! amup : flux de masse total des ascendances, par unite de
73  ! surface de la maille (kg/m2/s)
74  ! p : pressions aux milieux des couches (Pa)
75  ! ph : pressions aux interfaces (Pa)
76  ! pi : (p/p_0)**kapa (adim)
77  ! dtime: increment temporel (s)
78 
79  ! Variables internes :
80 
81  ! rhow : masse volumique de la poche froide
82  ! rho : environment density at P levels
83  ! rhoh : environment density at Ph levels
84  ! te : environment temperature | may change within
85  ! qe : environment humidity | sub-time-stepping
86  ! the : environment potential temperature
87  ! thu : potential temperature in undisturbed area
88  ! tu : temperature in undisturbed area
89  ! qu : humidity in undisturbed area
90  ! dp_omgb: vertical gradient og LS omega
91  ! omgbw : wake average vertical omega
92  ! dp_omgbw: vertical gradient of omgbw
93  ! omgbdq : flux of Delta_q transported by LS omega
94  ! dth : potential temperature diff. wake-undist.
95  ! th1 : first pot. temp. for vertical advection (=thu)
96  ! th2 : second pot. temp. for vertical advection (=thw)
97  ! q1 : first humidity for vertical advection
98  ! q2 : second humidity for vertical advection
99  ! d_deltatw : terme de redistribution pour deltatw
100  ! d_deltaqw : terme de redistribution pour deltaqw
101  ! deltatw0 : deltatw initial
102  ! deltaqw0 : deltaqw initial
103  ! hw0 : hw initial
104  ! sigmaw0: sigmaw initial
105  ! amflux : horizontal mass flux through wake boundary
106  ! wdens_ref: initial number of wakes per unit area (3D) or per
107  ! unit length (2D), at the beginning of each time step
108  ! Tgw : 1 sur la période de onde de gravité ! Cgw : vitesse de propagation de onde de gravité ! LL : distance entre 2 poches ! ------------------------------------------------------------------------- ! Déclaration de variables ! ------------------------------------------------------------------------- include "YOMCST.h" include "cvthermo.h" ! Arguments en entree ! -------------------- REAL, DIMENSION (klon, klev), INTENT(IN) :: p, pi REAL, DIMENSION (klon, klev+1), INTENT(IN) :: ph, omgb REAL, INTENT(IN) :: dtime REAL, DIMENSION (klon, klev), INTENT(IN) :: te0, qe0 REAL, DIMENSION (klon, klev), INTENT(IN) :: dtdwn, dqdwn REAL, DIMENSION (klon, klev), INTENT(IN) :: wdtpbl, wdqpbl, udtpbl, udqpbl ! UNUSED REAL, DIMENSION (klon, klev), INTENT(IN) :: amdwn, amup REAL, DIMENSION (klon, klev), INTENT(IN) :: dta, dqa REAL, DIMENSION (klon), INTENT(IN) :: sigd_con ! ! Input/Output REAL, DIMENSION (klon, klev), INTENT(INOUT) :: deltatw, deltaqw REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw ! Sorties ! -------- REAL, DIMENSION (klon, klev), INTENT(OUT) :: dth REAL, DIMENSION (klon, klev), INTENT(OUT) :: tu, qu REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtls, dqls REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtke, dqke REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtpbl, dqpbl REAL, DIMENSION (klon, klev), INTENT(OUT) :: spread REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltatw2, d_deltaqw2 REAL, DIMENSION (klon, klev+1), INTENT(OUT) :: omgbdth, omg REAL, DIMENSION (klon, klev), INTENT(OUT) :: dp_omgb, dp_deltomg REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_gw REAL, DIMENSION (klon), INTENT(OUT) :: hw, wape, fip, gfl, cstar REAL, DIMENSION (klon), INTENT(OUT) :: wdens INTEGER, DIMENSION (klon), INTENT(OUT) :: ktopw ! Variables internes ! ------------------- ! Variables à fixer REAL alon LOGICAL, SAVE :: first = .TRUE. !$OMP THREADPRIVATE(first) REAL, SAVE :: stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol) REAL delta_t_min INTEGER nsub REAL dtimesub REAL sigmad, hwmin, wapecut REAL :: sigmaw_max REAL :: dens_rate REAL wdens0 ! IM 080208 LOGICAL, DIMENSION (klon) :: gwake ! Variables de sauvegarde REAL, DIMENSION (klon, klev) :: deltatw0 REAL, DIMENSION (klon, klev) :: deltaqw0 REAL, DIMENSION (klon, klev) :: te, qe REAL, DIMENSION (klon) :: sigmaw0, sigmaw1 ! Variables pour les GW REAL, DIMENSION (klon) :: ll REAL, DIMENSION (klon, klev) :: n2 REAL, DIMENSION (klon, klev) :: cgw REAL, DIMENSION (klon, klev) :: tgw ! Variables liées au calcul de hw REAL, DIMENSION (klon) :: ptop_provis, ptop, ptop_new REAL, DIMENSION (klon) :: sum_dth REAL, DIMENSION (klon) :: dthmin REAL, DIMENSION (klon) :: z, dz, hw0 INTEGER, DIMENSION (klon) :: ktop, kupper ! Sub-timestep tendencies and related variables REAL d_deltatw(klon, klev), d_deltaqw(klon, klev) REAL d_te(klon, klev), d_qe(klon, klev) REAL d_sigmaw(klon), alpha(klon) REAL q0_min(klon), q1_min(klon) LOGICAL wk_adv(klon), ok_qx_qw(klon) REAL epsilon DATA epsilon/1.E-15/ ! Autres variables internes INTEGER isubstep, k, i REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu REAL, DIMENSION (klon) :: sum_dq, sum_rho REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn REAL, DIMENSION (klon) :: av_thu, av_tu, av_qu, av_thvu REAL, DIMENSION (klon) :: av_dth, av_dq, av_rho REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn REAL, DIMENSION (klon, klev) :: rho, rhow REAL, DIMENSION (klon, klev+1) :: rhoh REAL, DIMENSION (klon, klev) :: rhow_moyen REAL, DIMENSION (klon, klev) :: zh REAL, DIMENSION (klon, klev+1) :: zhh REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2 REAL, DIMENSION (klon, klev) :: the, thu ! REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw REAL, DIMENSION (klon, klev+1) :: omgbw REAL, DIMENSION (klon) :: pupper REAL, DIMENSION (klon) :: omgtop REAL, DIMENSION (klon, klev) :: dp_omgbw REAL, DIMENSION (klon) :: ztop, dztop REAL, DIMENSION (klon, klev) :: alpha_up REAL, DIMENSION (klon) :: rre1, rre2 REAL :: rrd1, rrd2 REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2 REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq REAL, DIMENSION (klon, klev) :: omgbdq REAL, DIMENSION (klon) :: ff, gg REAL, DIMENSION (klon) :: wape2, cstar2, heff REAL, DIMENSION (klon, klev) :: crep REAL, DIMENSION (klon, klev) :: ppi ! cc nrlmd REAL, DIMENSION (klon) :: death_rate, nat_rate REAL, DIMENSION (klon, klev) :: entr REAL, DIMENSION (klon, klev) :: detr ! ------------------------------------------------------------------------- ! Initialisations ! ------------------------------------------------------------------------- ! print*, 'wake initialisations' ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10. ! ------------------------------------------------------------------------- DATA wapecut, sigmad, hwmin/5., .02, 10./ ! cc nrlmd DATA sigmaw_max/0.4/ DATA dens_rate/0.1/ ! cc ! Longueur de maille (en m) ! ------------------------------------------------------------------------- ! ALON = 3.e5 alon = 1.E6 ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei) ! coefgw : Coefficient pour les ondes de gravité ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) ! wdens : Densité de poche froide par maille ! ------------------------------------------------------------------------- ! cc nrlmd coefgw=10 ! coefgw=1 ! wdens0 = 1.0/(alon**2) ! cc nrlmd wdens = 1.0/(alon**2) ! cc nrlmd stark = 0.50 ! CRtest ! cc nrlmd alpk=0.1 ! alpk = 1.0 ! alpk = 0.5 ! alpk = 0.05 if (first) then stark = 0.33 alpk = 0.25 wdens_ref = 8.E-12 coefgw = 4. crep_upper = 0.9 crep_sol = 1.0 ! cc nrlmd Lecture du fichier wake_param.data !$OMP MASTER OPEN (99, FILE='wake_param.data', STATUS='old', FORM='formatted', ERR=9999) READ (99, *, END=9998) stark READ (99, *, END=9998) alpk READ (99, *, END=9998) wdens_ref READ (99, *, END=9998) coefgw 9998 CONTINUE CLOSE (99) 9999 CONTINUE !$OMP END MASTER CALL bcast(stark) CALL bcast(alpk) CALL bcast(wdens_ref) CALL bcast(coefgw) first=.false. endif ! Initialisation de toutes des densites a wdens_ref. ! Les densites peuvent evoluer si les poches debordent ! (voir au tout debut de la boucle sur les substeps) wdens = wdens_ref ! print*,'stark',stark ! print*,'alpk',alpk ! print*,'wdens',wdens ! print*,'coefgw',coefgw ! cc ! Minimum value for |T_wake - T_undist|. Used for wake top definition ! ------------------------------------------------------------------------- delta_t_min = 0.2 ! 1. - Save initial values and initialize tendencies ! -------------------------------------------------- DO k = 1, klev DO i = 1, klon ppi(i, k) = pi(i, k) deltatw0(i, k) = deltatw(i, k) deltaqw0(i, k) = deltaqw(i, k) te(i, k) = te0(i, k) qe(i, k) = qe0(i, k) dtls(i, k) = 0. dqls(i, k) = 0. d_deltat_gw(i, k) = 0. d_te(i, k) = 0. d_qe(i, k) = 0. d_deltatw(i, k) = 0. d_deltaqw(i, k) = 0. ! IM 060508 beg d_deltatw2(i, k) = 0. d_deltaqw2(i, k) = 0. ! IM 060508 end END DO END DO ! sigmaw1=sigmaw ! IF (sigd_con.GT.sigmaw1) THEN ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con ! ENDIF DO i = 1, klon ! c sigmaw(i) = amax1(sigmaw(i),sigd_con(i)) sigmaw(i) = amax1(sigmaw(i), sigmad) sigmaw(i) = amin1(sigmaw(i), 0.99) sigmaw0(i) = sigmaw(i) wape(i) = 0. wape2(i) = 0. d_sigmaw(i) = 0. ktopw(i) = 0 END DO ! 2. - Prognostic part ! -------------------- ! 2.1 - Undisturbed area and Wake integrals ! --------------------------------------------------------- DO i = 1, klon z(i) = 0. ktop(i) = 0 kupper(i) = 0 sum_thu(i) = 0. sum_tu(i) = 0. sum_qu(i) = 0. sum_thvu(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_rho(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thu(i) = 0. av_tu(i) = 0. av_qu(i) = 0. av_thvu(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_rho(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END DO ! Distance between wakes DO i = 1, klon ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i)) END DO ! Potential temperatures and humidity ! ---------------------------------------------------------- DO k = 1, klev DO i = 1, klon ! write(*,*)'wake 1',i,k,rd,te(i,k) rho(i, k) = p(i, k)/(rd*te(i,k)) ! write(*,*)'wake 2',rho(i,k) IF (k==1) THEN ! write(*,*)'wake 3',i,k,rd,te(i,k) rhoh(i, k) = ph(i, k)/(rd*te(i,k)) ! write(*,*)'wake 4',i,k,rd,te(i,k) zhh(i, k) = 0 ELSE ! write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1)) rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1))) ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1) zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1) END IF ! write(*,*)'wake 7',ppi(i,k) the(i, k) = te(i, k)/ppi(i, k) thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i) qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k))) rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k))) dth(i, k) = deltatw(i, k)/ppi(i, k) END DO END DO DO k = 1, klev - 1 DO i = 1, klon IF (k==1) THEN n2(i, k) = 0 ELSE n2(i, k) = amax1(0., -rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i, & k-1))/(p(i,k+1)-p(i,k-1))) END IF zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2 cgw(i, k) = sqrt(n2(i,k))*zh(i, k) tgw(i, k) = coefgw*cgw(i, k)/ll(i) END DO END DO DO i = 1, klon n2(i, klev) = 0 zh(i, klev) = 0 cgw(i, klev) = 0 tgw(i, klev) = 0 END DO ! Calcul de la masse volumique moyenne de la colonne (bdlmd) ! ----------------------------------------------------------------- DO k = 1, klev DO i = 1, klon epaisseur1(i, k) = 0. epaisseur2(i, k) = 0. END DO END DO DO i = 1, klon epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1. epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1. rhow_moyen(i, 1) = rhow(i, 1) END DO DO k = 2, klev DO i = 1, klon epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1. epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k) rhow_moyen(i, k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+rhow(i,k)* & epaisseur1(i,k))/epaisseur2(i, k) END DO END DO ! Choose an integration bound well above wake top ! ----------------------------------------------------------------- ! Pupper = 50000. ! melting level ! Pupper = 60000. ! Pupper = 80000. ! essais pour case_e DO i = 1, klon pupper(i) = 0.6*ph(i, 1) pupper(i) = max(pupper(i), 45000.) ! cc Pupper(i) = 60000. END DO ! Determine Wake top pressure (Ptop) from buoyancy integral ! -------------------------------------------------------- ! -1/ Pressure of the level where dth becomes less than delta_t_min. DO i = 1, klon ptop_provis(i) = ph(i, 1) END DO DO k = 2, klev DO i = 1, klon ! IM v3JYG; ptop_provis(i).LT. ph(i,1) IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. & ptop_provis(i)==ph(i,1)) THEN ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO ! -2/ dth integral DO i = 1, klon sum_dth(i) = 0. dthmin(i) = -delta_t_min z(i) = 0. END DO DO k = 1, klev DO i = 1, klon dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) dthmin(i) = amin1(dthmin(i), dth(i,k)) END IF END DO END DO ! -3/ height of triangle with area= sum_dth and base = dthmin DO i = 1, klon hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) hw0(i) = amax1(hwmin, hw0(i)) END DO ! -4/ now, get Ptop DO i = 1, klon z(i) = 0. ptop(i) = ph(i, 1) END DO DO k = 1, klev DO i = 1, klon dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i)) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i) END IF END DO END DO ! -5/ Determination de ktop et kupper DO k = klev, 1, -1 DO i = 1, klon IF (ph(i,k+1)<ptop(i)) ktop(i) = k IF (ph(i,k+1)<pupper(i)) kupper(i) = k END DO END DO ! On evite kupper = 1 et kupper = klev DO i = 1, klon kupper(i) = max(kupper(i), 2) kupper(i) = min(kupper(i), klev-1) END DO ! -6/ Correct ktop and ptop DO i = 1, klon ptop_new(i) = ptop(i) END DO DO k = klev, 2, -1 DO i = 1, klon IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO DO i = 1, klon ptop(i) = ptop_new(i) END DO DO k = klev, 1, -1 DO i = 1, klon IF (ph(i,k+1)<ptop(i)) ktop(i) = k END DO END DO ! -5/ Set deltatw & deltaqw to 0 above kupper DO k = 1, klev DO i = 1, klon IF (k>=kupper(i)) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. END IF END DO END DO ! Vertical gradient of LS omega DO k = 1, klev DO i = 1, klon IF (k<=kupper(i)) THEN dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k)) END IF END DO END DO ! Integrals (and wake top level number) ! -------------------------------------- ! Initialize sum_thvu to 1st level virt. pot. temp. DO i = 1, klon z(i) = 1. dz(i) = 1. sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i) sum_dth(i) = 0. END DO DO k = 1, klev DO i = 1, klon dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) END IF END DO END DO DO i = 1, klon hw0(i) = z(i) END DO ! 2.1 - WAPE and mean forcing computation ! --------------------------------------- ! --------------------------------------- ! Means DO i = 1, klon av_thu(i) = sum_thu(i)/hw0(i) av_tu(i) = sum_tu(i)/hw0(i) av_qu(i) = sum_qu(i)/hw0(i) av_thvu(i) = sum_thvu(i)/hw0(i) ! av_thve = sum_thve/hw0 av_dth(i) = sum_dth(i)/hw0(i) av_dq(i) = sum_dq(i)/hw0(i) av_rho(i) = sum_rho(i)/hw0(i) av_dtdwn(i) = sum_dtdwn(i)/hw0(i) av_dqdwn(i) = sum_dqdwn(i)/hw0(i) wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i & )+av_dth(i)*av_dq(i)))/av_thvu(i) END DO ! 2.2 Prognostic variable update ! ------------------------------ ! Filter out bad wakes DO k = 1, klev DO i = 1, klon IF (wape(i)<0.) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. END IF END DO END DO DO i = 1, klon IF (wape(i)<0.) THEN wape(i) = 0. cstar(i) = 0. hw(i) = hwmin sigmaw(i) = amax1(sigmad, sigd_con(i)) fip(i) = 0. gwake(i) = .FALSE. ELSE cstar(i) = stark*sqrt(2.*wape(i)) gwake(i) = .TRUE. END IF END DO ! Check qx and qw positivity ! -------------------------- DO i = 1, klon q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), (qe(i, & 1)+(1.-sigmaw(i))*deltaqw(i,1))) END DO DO k = 2, klev DO i = 1, klon q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), (qe(i, & k)+(1.-sigmaw(i))*deltaqw(i,k))) IF (q1_min(i)<=q0_min(i)) THEN q0_min(i) = q1_min(i) END IF END DO END DO DO i = 1, klon ok_qx_qw(i) = q0_min(i) >= 0. alpha(i) = 1. END DO ! C ----------------------------------------------------------------- ! Sub-time-stepping ! ----------------- nsub = 10 dtimesub = dtime/nsub ! ------------------------------------------------------------ DO isubstep = 1, nsub ! ------------------------------------------------------------ ! wk_adv is the logical flag enabling wake evolution in the time advance ! loop DO i = 1, klon wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1. END DO ! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement ! négatif de ktop à kupper -------- ! cc On calcule pour cela une densité wdens0 pour laquelle on ! aurait un entrainement nul --- DO i = 1, klon ! c print *,' isubstep,wk_adv(i),cstar(i),wape(i) ', ! c $ isubstep,wk_adv(i),cstar(i),wape(i) IF (wk_adv(i) .AND. cstar(i)>0.01) THEN omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + & rg*amup(i, kupper(i)+1)/(1.-sigmaw(i)) wdens0 = (sigmaw(i)/(4.*3.14))*((1.-sigmaw(i))*omg(i,kupper(i)+1)/(( & ph(i,1)-pupper(i))*cstar(i)))**(2) IF (wdens(i)<=wdens0*1.1) THEN wdens(i) = wdens0 END IF ! c print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i) ! c $ ,ph(i,1)-pupper(i)', ! c $ omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i) ! c $ ,ph(i,1)-pupper(i) END IF END DO ! cc nrlmd DO i = 1, klon IF (wk_adv(i)) THEN gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i)) sigmaw(i) = amin1(sigmaw(i), sigmaw_max) END IF END DO DO i = 1, klon IF (wk_adv(i)) THEN ! cc nrlmd Introduction du taux de mortalité des poches et ! test sur sigmaw_max=0.4 ! cc d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub IF (sigmaw(i)>=sigmaw_max) THEN death_rate(i) = gfl(i)*cstar(i)/sigmaw(i) ELSE death_rate(i) = 0. END IF d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* & dtimesub ! $ - nat_rate(i)*sigmaw(i)*dtimesub ! c print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), ! c $ death_rate(i),ktop(i),kupper(i)', ! c $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), ! c $ death_rate(i),ktop(i),kupper(i) ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub ! sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!! ! wdens = wdens0/(10.*sigmaw) ! sigmaw =max(sigmaw,sigd_con) ! sigmaw =max(sigmaw,sigmad) END IF END DO ! calcul de la difference de vitesse verticale poche - zone non perturbee ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on ! definit ! IM 060208 au niveau k=1..? DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd dp_deltomg(i, k) = 0. END IF END DO END DO DO k = 1, klev + 1 DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd omg(i, k) = 0. END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN z(i) = 0. omg(i, 1) = 0. dp_deltomg(i, 1) = -(gfl(i)*cstar(i))/(sigmaw(i)*(1-sigmaw(i))) END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=ktop(i)) THEN dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg) z(i) = z(i) + dz(i) dp_deltomg(i, k) = dp_deltomg(i, 1) omg(i, k) = dp_deltomg(i, 1)*z(i) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg) ztop(i) = z(i) + dztop(i) omgtop(i) = dp_deltomg(i, 1)*ztop(i) END IF END DO ! ----------------- ! From m/s to Pa/s ! ----------------- DO i = 1, klon IF (wk_adv(i)) THEN omgtop(i) = -rho(i, ktop(i))*rg*omgtop(i) dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1)) END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=ktop(i)) THEN omg(i, k) = -rho(i, k)*rg*omg(i, k) dp_deltomg(i, k) = dp_deltomg(i, 1) END IF END DO END DO ! raccordement lineaire de omg de ptop a pupper DO i = 1, klon IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + & rg*amup(i, kupper(i)+1)/(1.-sigmaw(i)) dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ & (ptop(i)-pupper(i)) END IF END DO ! c DO i=1,klon ! c print*,'Pente entre 0 et kupper (référence)' ! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1)) ! c print*,'Pente entre ktop et kupper' ! c $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i)) ! c ENDDO ! c DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN dp_deltomg(i, k) = dp_deltomg(i, kupper(i)) omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i)) END IF END DO END DO ! cc nrlmd ! c DO i=1,klon ! c print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1) ! c END DO ! cc ! -- Compute wake average vertical velocity omgbw DO k = 1, klev + 1 DO i = 1, klon IF (wk_adv(i)) THEN omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k) END IF END DO END DO ! -- and its vertical gradient dp_omgbw DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k)) END IF END DO END DO ! -- Upstream coefficients for omgb velocity ! -- (alpha_up(k) is the coefficient of the value at level k) ! -- (1-alpha_up(k) is the coefficient of the value at level k-1) DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN alpha_up(i, k) = 0. IF (omgb(i,k)>0.) alpha_up(i, k) = 1. END IF END DO END DO ! Matrix expressing [The,deltatw] from [Th1,Th2] DO i = 1, klon IF (wk_adv(i)) THEN rre1(i) = 1. - sigmaw(i) rre2(i) = sigmaw(i) END IF END DO rrd1 = -1. rrd2 = 1. ! -- Get [Th1,Th2], dth and [q1,q2] DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN dth(i, k) = deltatw(i, k)/ppi(i, k) th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd d_th1(i, 1) = 0. d_th2(i, 1) = 0. d_dth(i, 1) = 0. d_q1(i, 1) = 0. d_q2(i, 1) = 0. d_dq(i, 1) = 0. END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN d_th1(i, k) = th1(i, k-1) - th1(i, k) d_th2(i, k) = th2(i, k-1) - th2(i, k) d_dth(i, k) = dth(i, k-1) - dth(i, k) d_q1(i, k) = q1(i, k-1) - q1(i, k) d_q2(i, k) = q2(i, k-1) - q2(i, k) d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN omgbdth(i, 1) = 0. omgbdq(i, 1) = 0. END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN ! loop on interfaces omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k)) omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k)) END IF END DO END DO ! ----------------------------------------------------------------- DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN ! ----------------------------------------------------------------- ! Compute redistribution (advective) term d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* & (rrd1*omg(i,k)*sigmaw(i)*d_th1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( & i))*d_th2(i,k+1)-(1.-alpha_up(i,k))*omgbdth(i,k)-alpha_up(i,k+1)* & omgbdth(i,k+1))*ppi(i, k) ! print*,'d_deltatw=',d_deltatw(i,k) d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* & (rrd1*omg(i,k)*sigmaw(i)*d_q1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( & i))*d_q2(i,k+1)-(1.-alpha_up(i,k))*omgbdq(i,k)-alpha_up(i,k+1)* & omgbdq(i,k+1)) ! print*,'d_deltaqw=',d_deltaqw(i,k) ! and increment large scale tendencies ! C ! ----------------------------------------------------------------- d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, & k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/(ph(i,k)-ph(i, & k+1)) & ! cc nrlmd $ ! -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k) -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/(ph(i, & k)-ph(i,k+1)) & ! cc )*ppi(i, k) d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, & k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/(ph(i,k)-ph(i, & k+1)) & ! cc nrlmd $ ! -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k) -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i, & k+1))/(ph(i,k)-ph(i,k+1)) & ! cc ) ! cc nrlmd ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, & k)/(ph(i,k)-ph(i,k+1))))*ppi(i, k) d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, & k)/(ph(i,k)-ph(i,k+1)))) END IF ! cc END DO END DO ! ------------------------------------------------------------------ ! Increment state variables DO k = 1, klev DO i = 1, klon ! cc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN IF (wk_adv(i) .AND. k<=kupper(i)) THEN ! cc ! Coefficient de répartition crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ & (ph(i,kupper(i))-ph(i,1)) crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-ph(i & ,kupper(i))) ! Reintroduce compensating subsidence term. ! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)) ! . /(1-sigmaw) ! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)) ! . /(1-sigmaw) ! dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k)) ! . /(1-sigmaw) ! dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k)) ! . /(1-sigmaw) dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i))) dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i))) ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k) !jyg< !! !!--------------------------------------------------------------- !! The change of delta_T due to PBL (vertical diffusion plus thermal plumes) !! is accounted for by the PBL and the Thermals schemes. It is now set to zero !! within the Wake scheme. !!--------------------------------------------------------------- dtPBL(i,k) = 0. dqPBL(i,k) = 0. ! !! dtPBL(i,k)=wdtPBL(i,k) - udtPBL(i,k) !! dqPBL(i,k)=wdqPBL(i,k) - udqPBL(i,k) ! !! dtpbl(i, k) = (wdtpbl(i,k)/sigmaw(i)-udtpbl(i,k)/(1.-sigmaw(i))) !! dqpbl(i, k) = (wdqpbl(i,k)/sigmaw(i)-udqpbl(i,k)/(1.-sigmaw(i))) ! print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k) !>jyg ! ! cc nrlmd Prise en compte du taux de mortalité ! cc Définitions de entr, detr detr(i, k) = 0. entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + & sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k) spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i) ! cc spread(i,k) = ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/ ! cc $ sigmaw(i) ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU ! Jingmei ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k), ! & Tgw(i,k),deltatw(i,k) d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)*deltatw(i, k)* & dtimesub ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k) ff(i) = d_deltatw(i, k)/dtimesub ! Sans GW ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k)) ! GW formule 1 ! deltatw(k) = deltatw(k)+dtimesub* ! $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k)) ! GW formule 2 IF (dtimesub*tgw(i,k)<1.E-10) THEN d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc ! $ ! -spread(i,k)*deltatw(i,k) -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( & i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc -tgw(i,k)*deltatw(i,k)) ELSE d_deltatw(i, k) = 1/tgw(i, k)*(1-exp(-dtimesub*tgw(i, & k)))*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc $ ! -spread(i,k)*deltatw(i,k) -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( & i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc -tgw(i,k)*deltatw(i,k)) END IF dth(i, k) = deltatw(i, k)/ppi(i, k) gg(i) = d_deltaqw(i, k)/dtimesub d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k)+dqpbl(i,k) & ! cc $ ! -spread(i,k)*deltaqw(i,k)) -entr(i,k)*deltaqw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw(i)+detr( & i,k))*deltaqw(i,k)/(1.-sigmaw(i))) ! cc ! cc nrlmd ! cc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k) ! cc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k) ! cc END IF END DO END DO ! Scale tendencies so that water vapour remains positive in w and x. CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, & d_deltaqw, sigmaw, d_sigmaw, alpha) ! cc nrlmd ! c print*,'alpha' ! c do i=1,klon ! c print*,alpha(i) ! c end do ! cc DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN d_te(i, k) = alpha(i)*d_te(i, k) d_qe(i, k) = alpha(i)*d_qe(i, k) d_deltatw(i, k) = alpha(i)*d_deltatw(i, k) d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k) d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN d_sigmaw(i) = alpha(i)*d_sigmaw(i) END IF END DO ! Update large scale variables and wake variables ! IM 060208 manque DO i + remplace DO k=1,kupper(i) ! IM 060208 DO k = 1,kupper(i) DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN dtls(i, k) = dtls(i, k) + d_te(i, k) dqls(i, k) = dqls(i, k) + d_qe(i, k) ! cc nrlmd d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k) d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k) ! cc END IF END DO END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN te(i, k) = te0(i, k) + dtls(i, k) qe(i, k) = qe0(i, k) + dqls(i, k) the(i, k) = te(i, k)/ppi(i, k) deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k) deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k) dth(i, k) = deltatw(i, k)/ppi(i, k) ! c print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k) ! c $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN sigmaw(i) = sigmaw(i) + d_sigmaw(i) END IF END DO ! Determine Ptop from buoyancy integral ! --------------------------------------- ! - 1/ Pressure of the level where dth changes sign. DO i = 1, klon IF (wk_adv(i)) THEN ptop_provis(i) = ph(i, 1) END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. & dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO ! - 2/ dth integral DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd sum_dth(i) = 0. dthmin(i) = -delta_t_min z(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) dthmin(i) = amin1(dthmin(i), dth(i,k)) END IF END IF END DO END DO ! - 3/ height of triangle with area= sum_dth and base = dthmin DO i = 1, klon IF (wk_adv(i)) THEN hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) hw(i) = amax1(hwmin, hw(i)) END IF END DO ! - 4/ now, get Ptop DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd ktop(i) = 0 z(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i)) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i) ktop(i) = k END IF END IF END DO END DO ! 4.5/Correct ktop and ptop DO i = 1, klon IF (wk_adv(i)) THEN ptop_new(i) = ptop(i) END IF END DO DO k = klev, 2, -1 DO i = 1, klon ! IM v3JYG; IF (k .GE. ktop(i) IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN ptop(i) = ptop_new(i) END IF END DO DO k = klev, 1, -1 DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (ph(i,k+1)<ptop(i)) ktop(i) = k END IF END DO END DO ! 5/ Set deltatw & deltaqw to 0 above kupper DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k>=kupper(i)) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. END IF END DO END DO ! -------------Cstar computation--------------------------------- DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd sum_thu(i) = 0. sum_tu(i) = 0. sum_qu(i) = 0. sum_thvu(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_rho(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thu(i) = 0. av_tu(i) = 0. av_qu(i) = 0. av_thvu(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_rho(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END IF END DO ! Integrals (and wake top level number) ! -------------------------------------- ! Initialize sum_thvu to 1st level virt. pot. temp. DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd z(i) = 1. dz(i) = 1. sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i) sum_dth(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) END IF END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd hw0(i) = z(i) END IF END DO ! - WAPE and mean forcing computation ! --------------------------------------- ! --------------------------------------- ! Means DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd av_thu(i) = sum_thu(i)/hw0(i) av_tu(i) = sum_tu(i)/hw0(i) av_qu(i) = sum_qu(i)/hw0(i) av_thvu(i) = sum_thvu(i)/hw0(i) av_dth(i) = sum_dth(i)/hw0(i) av_dq(i) = sum_dq(i)/hw0(i) av_rho(i) = sum_rho(i)/hw0(i) av_dtdwn(i) = sum_dtdwn(i)/hw0(i) av_dqdwn(i) = sum_dqdwn(i)/hw0(i) wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* & av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i) END IF END DO ! Filter out bad wakes DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (wape(i)<0.) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. END IF END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (wape(i)<0.) THEN wape(i) = 0. cstar(i) = 0. hw(i) = hwmin sigmaw(i) = max(sigmad, sigd_con(i)) fip(i) = 0. gwake(i) = .FALSE. ELSE cstar(i) = stark*sqrt(2.*wape(i)) gwake(i) = .TRUE. END IF END IF END DO END DO ! end sub-timestep loop ! ----------------------------------------------------------------- ! Get back to tendencies per second DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN ! cc dtls(i, k) = dtls(i, k)/dtime dqls(i, k) = dqls(i, k)/dtime d_deltatw2(i, k) = d_deltatw2(i, k)/dtime d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime ! c print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k) ! c $ ,death_rate(i)*sigmaw(i) END IF END DO END DO ! ---------------------------------------------------------- ! Determine wake final state; recompute wape, cstar, ktop; ! filter out bad wakes. ! ---------------------------------------------------------- ! 2.1 - Undisturbed area and Wake integrals ! --------------------------------------------------------- DO i = 1, klon ! cc nrlmd if (wk_adv(i)) then !!! nrlmd IF (ok_qx_qw(i)) THEN ! cc z(i) = 0. sum_thu(i) = 0. sum_tu(i) = 0. sum_qu(i) = 0. sum_thvu(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_rho(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thu(i) = 0. av_tu(i) = 0. av_qu(i) = 0. av_thvu(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_rho(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END IF END DO ! Potential temperatures and humidity ! ---------------------------------------------------------- DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc rho(i, k) = p(i, k)/(rd*te(i,k)) IF (k==1) THEN rhoh(i, k) = ph(i, k)/(rd*te(i,k)) zhh(i, k) = 0 ELSE rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1))) zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1) END IF the(i, k) = te(i, k)/ppi(i, k) thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i) qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k))) dth(i, k) = deltatw(i, k)/ppi(i, k) END IF END DO END DO ! Integrals (and wake top level number) ! ----------------------------------------------------------- ! Initialize sum_thvu to 1st level virt. pot. temp. DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc z(i) = 1. dz(i) = 1. sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i) sum_dth(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) END IF END IF END DO END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc hw0(i) = z(i) END IF END DO ! - WAPE and mean forcing computation ! ------------------------------------------------------------- ! Means DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc av_thu(i) = sum_thu(i)/hw0(i) av_tu(i) = sum_tu(i)/hw0(i) av_qu(i) = sum_qu(i)/hw0(i) av_thvu(i) = sum_thvu(i)/hw0(i) av_dth(i) = sum_dth(i)/hw0(i) av_dq(i) = sum_dq(i)/hw0(i) av_rho(i) = sum_rho(i)/hw0(i) av_dtdwn(i) = sum_dtdwn(i)/hw0(i) av_dqdwn(i) = sum_dqdwn(i)/hw0(i) wape2(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* & av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i) END IF END DO ! Prognostic variable update ! ------------------------------------------------------------ ! Filter out bad wakes DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN ! cc deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. END IF END DO END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc IF (wape2(i)<0.) THEN wape2(i) = 0. cstar2(i) = 0. hw(i) = hwmin sigmaw(i) = amax1(sigmad, sigd_con(i)) fip(i) = 0. gwake(i) = .FALSE. ELSE IF (prt_level>=10) PRINT *, 'wape2>0' cstar2(i) = stark*sqrt(2.*wape2(i)) gwake(i) = .TRUE. END IF END IF END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc ktopw(i) = ktop(i) END IF END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc IF (ktopw(i)>0 .AND. gwake(i)) THEN ! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer) ! cc heff = 600. ! Utilisation de la hauteur hw ! c heff = 0.7*hw heff(i) = hw(i) fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* & sqrt(sigmaw(i)*wdens(i)*3.14) fip(i) = alpk*fip(i) ! jyg2 ELSE fip(i) = 0. END IF END IF END DO ! Limitation de sigmaw ! cc nrlmd ! DO i=1,klon ! IF (OK_qx_qw(i)) THEN ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max ! ENDIF ! ENDDO ! cc DO k = 1, klev DO i = 1, klon ! cc nrlmd On maintient désormais constant sigmaw en régime ! permanent ! cc IF ((sigmaw(i).GT.sigmaw_max).or. IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. & .NOT. ok_qx_qw(i)) THEN ! cc dtls(i, k) = 0. dqls(i, k) = 0. deltatw(i, k) = 0. deltaqw(i, k) = 0. END IF END DO END DO ! cc nrlmd On maintient désormais constant sigmaw en régime permanent DO i = 1, klon IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. & .NOT. ok_qx_qw(i)) THEN wape(i) = 0. cstar(i) = 0. !!jyg Outside subroutine "Wake" hw and sigmaw are zero when there are no wakes !! hw(i) = hwmin !jyg !! sigmaw(i) = sigmad !jyg hw(i) = 0. !jyg sigmaw(i) = 0. !jyg fip(i) = 0. ELSE wape(i) = wape2(i) cstar(i) = cstar2(i) END IF ! c print*,'wape wape2 ktopw OK_qx_qw =', ! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i) END DO RETURN END SUBROUTINE wake SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, & d_deltaqw, sigmaw, d_sigmaw, alpha) ! ------------------------------------------------------ ! Dtermination du coefficient alpha tel que les tendances ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent ! a une humidite positive dans la zone (x) et dans la zone (w). ! ------------------------------------------------------ IMPLICIT NONE ! Input REAL qe(nlon, nl), d_qe(nlon, nl) REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl) REAL sigmaw(nlon), d_sigmaw(nlon) LOGICAL wk_adv(nlon) INTEGER nl, nlon ! Output REAL alpha(nlon) ! Internal variables REAL zeta(nlon, nl) REAL alpha1(nlon) REAL x, a, b, c, discrim REAL epsilon ! DATA epsilon/1.e-15/ INTEGER i,k DO k = 1, nl DO i = 1, nlon IF (wk_adv(i)) THEN IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN zeta(i, k) = 0. ELSE zeta(i, k) = 1. END IF END IF END DO DO i = 1, nlon IF (wk_adv(i)) THEN x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + & (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i)*(deltaqw(i,k)+ & d_deltaqw(i,k)) a = -d_sigmaw(i)*d_deltaqw(i, k) b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - & deltaqw(i, k)*d_sigmaw(i) c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon discrim = b*b - 4.*a*c ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap alpha1(i) = 1. ELSE IF (x>=0.) THEN alpha1(i) = 1. ELSE IF (a>0.) THEN alpha1(i) = 0.9*min((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim & ))/(2.*a)) ELSE IF (a==0.) THEN alpha1(i) = 0.9*(-c/b) ELSE ! print*,'a,b,c discrim',a,b,c discrim alpha1(i) = 0.9*max((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim & ))/(2.*a)) END IF END IF END IF alpha(i) = min(alpha(i), alpha1(i)) END IF END DO END DO RETURN END SUBROUTINE wake_vec_modulation
109  ! Cgw : vitesse de propagation de onde de gravité ! LL : distance entre 2 poches ! ------------------------------------------------------------------------- ! Déclaration de variables ! ------------------------------------------------------------------------- include "YOMCST.h" include "cvthermo.h" ! Arguments en entree ! -------------------- REAL, DIMENSION (klon, klev), INTENT(IN) :: p, pi REAL, DIMENSION (klon, klev+1), INTENT(IN) :: ph, omgb REAL, INTENT(IN) :: dtime REAL, DIMENSION (klon, klev), INTENT(IN) :: te0, qe0 REAL, DIMENSION (klon, klev), INTENT(IN) :: dtdwn, dqdwn REAL, DIMENSION (klon, klev), INTENT(IN) :: wdtpbl, wdqpbl, udtpbl, udqpbl ! UNUSED REAL, DIMENSION (klon, klev), INTENT(IN) :: amdwn, amup REAL, DIMENSION (klon, klev), INTENT(IN) :: dta, dqa REAL, DIMENSION (klon), INTENT(IN) :: sigd_con ! ! Input/Output REAL, DIMENSION (klon, klev), INTENT(INOUT) :: deltatw, deltaqw REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw ! Sorties ! -------- REAL, DIMENSION (klon, klev), INTENT(OUT) :: dth REAL, DIMENSION (klon, klev), INTENT(OUT) :: tu, qu REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtls, dqls REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtke, dqke REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtpbl, dqpbl REAL, DIMENSION (klon, klev), INTENT(OUT) :: spread REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltatw2, d_deltaqw2 REAL, DIMENSION (klon, klev+1), INTENT(OUT) :: omgbdth, omg REAL, DIMENSION (klon, klev), INTENT(OUT) :: dp_omgb, dp_deltomg REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_gw REAL, DIMENSION (klon), INTENT(OUT) :: hw, wape, fip, gfl, cstar REAL, DIMENSION (klon), INTENT(OUT) :: wdens INTEGER, DIMENSION (klon), INTENT(OUT) :: ktopw ! Variables internes ! ------------------- ! Variables à fixer REAL alon LOGICAL, SAVE :: first = .TRUE. !$OMP THREADPRIVATE(first) REAL, SAVE :: stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol) REAL delta_t_min INTEGER nsub REAL dtimesub REAL sigmad, hwmin, wapecut REAL :: sigmaw_max REAL :: dens_rate REAL wdens0 ! IM 080208 LOGICAL, DIMENSION (klon) :: gwake ! Variables de sauvegarde REAL, DIMENSION (klon, klev) :: deltatw0 REAL, DIMENSION (klon, klev) :: deltaqw0 REAL, DIMENSION (klon, klev) :: te, qe REAL, DIMENSION (klon) :: sigmaw0, sigmaw1 ! Variables pour les GW REAL, DIMENSION (klon) :: ll REAL, DIMENSION (klon, klev) :: n2 REAL, DIMENSION (klon, klev) :: cgw REAL, DIMENSION (klon, klev) :: tgw ! Variables liées au calcul de hw REAL, DIMENSION (klon) :: ptop_provis, ptop, ptop_new REAL, DIMENSION (klon) :: sum_dth REAL, DIMENSION (klon) :: dthmin REAL, DIMENSION (klon) :: z, dz, hw0 INTEGER, DIMENSION (klon) :: ktop, kupper ! Sub-timestep tendencies and related variables REAL d_deltatw(klon, klev), d_deltaqw(klon, klev) REAL d_te(klon, klev), d_qe(klon, klev) REAL d_sigmaw(klon), alpha(klon) REAL q0_min(klon), q1_min(klon) LOGICAL wk_adv(klon), ok_qx_qw(klon) REAL epsilon DATA epsilon/1.E-15/ ! Autres variables internes INTEGER isubstep, k, i REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu REAL, DIMENSION (klon) :: sum_dq, sum_rho REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn REAL, DIMENSION (klon) :: av_thu, av_tu, av_qu, av_thvu REAL, DIMENSION (klon) :: av_dth, av_dq, av_rho REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn REAL, DIMENSION (klon, klev) :: rho, rhow REAL, DIMENSION (klon, klev+1) :: rhoh REAL, DIMENSION (klon, klev) :: rhow_moyen REAL, DIMENSION (klon, klev) :: zh REAL, DIMENSION (klon, klev+1) :: zhh REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2 REAL, DIMENSION (klon, klev) :: the, thu ! REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw REAL, DIMENSION (klon, klev+1) :: omgbw REAL, DIMENSION (klon) :: pupper REAL, DIMENSION (klon) :: omgtop REAL, DIMENSION (klon, klev) :: dp_omgbw REAL, DIMENSION (klon) :: ztop, dztop REAL, DIMENSION (klon, klev) :: alpha_up REAL, DIMENSION (klon) :: rre1, rre2 REAL :: rrd1, rrd2 REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2 REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq REAL, DIMENSION (klon, klev) :: omgbdq REAL, DIMENSION (klon) :: ff, gg REAL, DIMENSION (klon) :: wape2, cstar2, heff REAL, DIMENSION (klon, klev) :: crep REAL, DIMENSION (klon, klev) :: ppi ! cc nrlmd REAL, DIMENSION (klon) :: death_rate, nat_rate REAL, DIMENSION (klon, klev) :: entr REAL, DIMENSION (klon, klev) :: detr ! ------------------------------------------------------------------------- ! Initialisations ! ------------------------------------------------------------------------- ! print*, 'wake initialisations' ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10. ! ------------------------------------------------------------------------- DATA wapecut, sigmad, hwmin/5., .02, 10./ ! cc nrlmd DATA sigmaw_max/0.4/ DATA dens_rate/0.1/ ! cc ! Longueur de maille (en m) ! ------------------------------------------------------------------------- ! ALON = 3.e5 alon = 1.E6 ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei) ! coefgw : Coefficient pour les ondes de gravité ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) ! wdens : Densité de poche froide par maille ! ------------------------------------------------------------------------- ! cc nrlmd coefgw=10 ! coefgw=1 ! wdens0 = 1.0/(alon**2) ! cc nrlmd wdens = 1.0/(alon**2) ! cc nrlmd stark = 0.50 ! CRtest ! cc nrlmd alpk=0.1 ! alpk = 1.0 ! alpk = 0.5 ! alpk = 0.05 if (first) then stark = 0.33 alpk = 0.25 wdens_ref = 8.E-12 coefgw = 4. crep_upper = 0.9 crep_sol = 1.0 ! cc nrlmd Lecture du fichier wake_param.data !$OMP MASTER OPEN (99, FILE='wake_param.data', STATUS='old', FORM='formatted', ERR=9999) READ (99, *, END=9998) stark READ (99, *, END=9998) alpk READ (99, *, END=9998) wdens_ref READ (99, *, END=9998) coefgw 9998 CONTINUE CLOSE (99) 9999 CONTINUE !$OMP END MASTER CALL bcast(stark) CALL bcast(alpk) CALL bcast(wdens_ref) CALL bcast(coefgw) first=.false. endif ! Initialisation de toutes des densites a wdens_ref. ! Les densites peuvent evoluer si les poches debordent ! (voir au tout debut de la boucle sur les substeps) wdens = wdens_ref ! print*,'stark',stark ! print*,'alpk',alpk ! print*,'wdens',wdens ! print*,'coefgw',coefgw ! cc ! Minimum value for |T_wake - T_undist|. Used for wake top definition ! ------------------------------------------------------------------------- delta_t_min = 0.2 ! 1. - Save initial values and initialize tendencies ! -------------------------------------------------- DO k = 1, klev DO i = 1, klon ppi(i, k) = pi(i, k) deltatw0(i, k) = deltatw(i, k) deltaqw0(i, k) = deltaqw(i, k) te(i, k) = te0(i, k) qe(i, k) = qe0(i, k) dtls(i, k) = 0. dqls(i, k) = 0. d_deltat_gw(i, k) = 0. d_te(i, k) = 0. d_qe(i, k) = 0. d_deltatw(i, k) = 0. d_deltaqw(i, k) = 0. ! IM 060508 beg d_deltatw2(i, k) = 0. d_deltaqw2(i, k) = 0. ! IM 060508 end END DO END DO ! sigmaw1=sigmaw ! IF (sigd_con.GT.sigmaw1) THEN ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con ! ENDIF DO i = 1, klon ! c sigmaw(i) = amax1(sigmaw(i),sigd_con(i)) sigmaw(i) = amax1(sigmaw(i), sigmad) sigmaw(i) = amin1(sigmaw(i), 0.99) sigmaw0(i) = sigmaw(i) wape(i) = 0. wape2(i) = 0. d_sigmaw(i) = 0. ktopw(i) = 0 END DO ! 2. - Prognostic part ! -------------------- ! 2.1 - Undisturbed area and Wake integrals ! --------------------------------------------------------- DO i = 1, klon z(i) = 0. ktop(i) = 0 kupper(i) = 0 sum_thu(i) = 0. sum_tu(i) = 0. sum_qu(i) = 0. sum_thvu(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_rho(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thu(i) = 0. av_tu(i) = 0. av_qu(i) = 0. av_thvu(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_rho(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END DO ! Distance between wakes DO i = 1, klon ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i)) END DO ! Potential temperatures and humidity ! ---------------------------------------------------------- DO k = 1, klev DO i = 1, klon ! write(*,*)'wake 1',i,k,rd,te(i,k) rho(i, k) = p(i, k)/(rd*te(i,k)) ! write(*,*)'wake 2',rho(i,k) IF (k==1) THEN ! write(*,*)'wake 3',i,k,rd,te(i,k) rhoh(i, k) = ph(i, k)/(rd*te(i,k)) ! write(*,*)'wake 4',i,k,rd,te(i,k) zhh(i, k) = 0 ELSE ! write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1)) rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1))) ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1) zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1) END IF ! write(*,*)'wake 7',ppi(i,k) the(i, k) = te(i, k)/ppi(i, k) thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i) qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k))) rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k))) dth(i, k) = deltatw(i, k)/ppi(i, k) END DO END DO DO k = 1, klev - 1 DO i = 1, klon IF (k==1) THEN n2(i, k) = 0 ELSE n2(i, k) = amax1(0., -rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i, & k-1))/(p(i,k+1)-p(i,k-1))) END IF zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2 cgw(i, k) = sqrt(n2(i,k))*zh(i, k) tgw(i, k) = coefgw*cgw(i, k)/ll(i) END DO END DO DO i = 1, klon n2(i, klev) = 0 zh(i, klev) = 0 cgw(i, klev) = 0 tgw(i, klev) = 0 END DO ! Calcul de la masse volumique moyenne de la colonne (bdlmd) ! ----------------------------------------------------------------- DO k = 1, klev DO i = 1, klon epaisseur1(i, k) = 0. epaisseur2(i, k) = 0. END DO END DO DO i = 1, klon epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1. epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1. rhow_moyen(i, 1) = rhow(i, 1) END DO DO k = 2, klev DO i = 1, klon epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1. epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k) rhow_moyen(i, k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+rhow(i,k)* & epaisseur1(i,k))/epaisseur2(i, k) END DO END DO ! Choose an integration bound well above wake top ! ----------------------------------------------------------------- ! Pupper = 50000. ! melting level ! Pupper = 60000. ! Pupper = 80000. ! essais pour case_e DO i = 1, klon pupper(i) = 0.6*ph(i, 1) pupper(i) = max(pupper(i), 45000.) ! cc Pupper(i) = 60000. END DO ! Determine Wake top pressure (Ptop) from buoyancy integral ! -------------------------------------------------------- ! -1/ Pressure of the level where dth becomes less than delta_t_min. DO i = 1, klon ptop_provis(i) = ph(i, 1) END DO DO k = 2, klev DO i = 1, klon ! IM v3JYG; ptop_provis(i).LT. ph(i,1) IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. & ptop_provis(i)==ph(i,1)) THEN ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO ! -2/ dth integral DO i = 1, klon sum_dth(i) = 0. dthmin(i) = -delta_t_min z(i) = 0. END DO DO k = 1, klev DO i = 1, klon dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) dthmin(i) = amin1(dthmin(i), dth(i,k)) END IF END DO END DO ! -3/ height of triangle with area= sum_dth and base = dthmin DO i = 1, klon hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) hw0(i) = amax1(hwmin, hw0(i)) END DO ! -4/ now, get Ptop DO i = 1, klon z(i) = 0. ptop(i) = ph(i, 1) END DO DO k = 1, klev DO i = 1, klon dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i)) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i) END IF END DO END DO ! -5/ Determination de ktop et kupper DO k = klev, 1, -1 DO i = 1, klon IF (ph(i,k+1)<ptop(i)) ktop(i) = k IF (ph(i,k+1)<pupper(i)) kupper(i) = k END DO END DO ! On evite kupper = 1 et kupper = klev DO i = 1, klon kupper(i) = max(kupper(i), 2) kupper(i) = min(kupper(i), klev-1) END DO ! -6/ Correct ktop and ptop DO i = 1, klon ptop_new(i) = ptop(i) END DO DO k = klev, 2, -1 DO i = 1, klon IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO DO i = 1, klon ptop(i) = ptop_new(i) END DO DO k = klev, 1, -1 DO i = 1, klon IF (ph(i,k+1)<ptop(i)) ktop(i) = k END DO END DO ! -5/ Set deltatw & deltaqw to 0 above kupper DO k = 1, klev DO i = 1, klon IF (k>=kupper(i)) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. END IF END DO END DO ! Vertical gradient of LS omega DO k = 1, klev DO i = 1, klon IF (k<=kupper(i)) THEN dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k)) END IF END DO END DO ! Integrals (and wake top level number) ! -------------------------------------- ! Initialize sum_thvu to 1st level virt. pot. temp. DO i = 1, klon z(i) = 1. dz(i) = 1. sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i) sum_dth(i) = 0. END DO DO k = 1, klev DO i = 1, klon dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) END IF END DO END DO DO i = 1, klon hw0(i) = z(i) END DO ! 2.1 - WAPE and mean forcing computation ! --------------------------------------- ! --------------------------------------- ! Means DO i = 1, klon av_thu(i) = sum_thu(i)/hw0(i) av_tu(i) = sum_tu(i)/hw0(i) av_qu(i) = sum_qu(i)/hw0(i) av_thvu(i) = sum_thvu(i)/hw0(i) ! av_thve = sum_thve/hw0 av_dth(i) = sum_dth(i)/hw0(i) av_dq(i) = sum_dq(i)/hw0(i) av_rho(i) = sum_rho(i)/hw0(i) av_dtdwn(i) = sum_dtdwn(i)/hw0(i) av_dqdwn(i) = sum_dqdwn(i)/hw0(i) wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i & )+av_dth(i)*av_dq(i)))/av_thvu(i) END DO ! 2.2 Prognostic variable update ! ------------------------------ ! Filter out bad wakes DO k = 1, klev DO i = 1, klon IF (wape(i)<0.) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. END IF END DO END DO DO i = 1, klon IF (wape(i)<0.) THEN wape(i) = 0. cstar(i) = 0. hw(i) = hwmin sigmaw(i) = amax1(sigmad, sigd_con(i)) fip(i) = 0. gwake(i) = .FALSE. ELSE cstar(i) = stark*sqrt(2.*wape(i)) gwake(i) = .TRUE. END IF END DO ! Check qx and qw positivity ! -------------------------- DO i = 1, klon q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), (qe(i, & 1)+(1.-sigmaw(i))*deltaqw(i,1))) END DO DO k = 2, klev DO i = 1, klon q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), (qe(i, & k)+(1.-sigmaw(i))*deltaqw(i,k))) IF (q1_min(i)<=q0_min(i)) THEN q0_min(i) = q1_min(i) END IF END DO END DO DO i = 1, klon ok_qx_qw(i) = q0_min(i) >= 0. alpha(i) = 1. END DO ! C ----------------------------------------------------------------- ! Sub-time-stepping ! ----------------- nsub = 10 dtimesub = dtime/nsub ! ------------------------------------------------------------ DO isubstep = 1, nsub ! ------------------------------------------------------------ ! wk_adv is the logical flag enabling wake evolution in the time advance ! loop DO i = 1, klon wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1. END DO ! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement ! négatif de ktop à kupper -------- ! cc On calcule pour cela une densité wdens0 pour laquelle on ! aurait un entrainement nul --- DO i = 1, klon ! c print *,' isubstep,wk_adv(i),cstar(i),wape(i) ', ! c $ isubstep,wk_adv(i),cstar(i),wape(i) IF (wk_adv(i) .AND. cstar(i)>0.01) THEN omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + & rg*amup(i, kupper(i)+1)/(1.-sigmaw(i)) wdens0 = (sigmaw(i)/(4.*3.14))*((1.-sigmaw(i))*omg(i,kupper(i)+1)/(( & ph(i,1)-pupper(i))*cstar(i)))**(2) IF (wdens(i)<=wdens0*1.1) THEN wdens(i) = wdens0 END IF ! c print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i) ! c $ ,ph(i,1)-pupper(i)', ! c $ omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i) ! c $ ,ph(i,1)-pupper(i) END IF END DO ! cc nrlmd DO i = 1, klon IF (wk_adv(i)) THEN gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i)) sigmaw(i) = amin1(sigmaw(i), sigmaw_max) END IF END DO DO i = 1, klon IF (wk_adv(i)) THEN ! cc nrlmd Introduction du taux de mortalité des poches et ! test sur sigmaw_max=0.4 ! cc d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub IF (sigmaw(i)>=sigmaw_max) THEN death_rate(i) = gfl(i)*cstar(i)/sigmaw(i) ELSE death_rate(i) = 0. END IF d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* & dtimesub ! $ - nat_rate(i)*sigmaw(i)*dtimesub ! c print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), ! c $ death_rate(i),ktop(i),kupper(i)', ! c $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), ! c $ death_rate(i),ktop(i),kupper(i) ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub ! sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!! ! wdens = wdens0/(10.*sigmaw) ! sigmaw =max(sigmaw,sigd_con) ! sigmaw =max(sigmaw,sigmad) END IF END DO ! calcul de la difference de vitesse verticale poche - zone non perturbee ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on ! definit ! IM 060208 au niveau k=1..? DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd dp_deltomg(i, k) = 0. END IF END DO END DO DO k = 1, klev + 1 DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd omg(i, k) = 0. END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN z(i) = 0. omg(i, 1) = 0. dp_deltomg(i, 1) = -(gfl(i)*cstar(i))/(sigmaw(i)*(1-sigmaw(i))) END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=ktop(i)) THEN dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg) z(i) = z(i) + dz(i) dp_deltomg(i, k) = dp_deltomg(i, 1) omg(i, k) = dp_deltomg(i, 1)*z(i) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg) ztop(i) = z(i) + dztop(i) omgtop(i) = dp_deltomg(i, 1)*ztop(i) END IF END DO ! ----------------- ! From m/s to Pa/s ! ----------------- DO i = 1, klon IF (wk_adv(i)) THEN omgtop(i) = -rho(i, ktop(i))*rg*omgtop(i) dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1)) END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=ktop(i)) THEN omg(i, k) = -rho(i, k)*rg*omg(i, k) dp_deltomg(i, k) = dp_deltomg(i, 1) END IF END DO END DO ! raccordement lineaire de omg de ptop a pupper DO i = 1, klon IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + & rg*amup(i, kupper(i)+1)/(1.-sigmaw(i)) dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ & (ptop(i)-pupper(i)) END IF END DO ! c DO i=1,klon ! c print*,'Pente entre 0 et kupper (référence)' ! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1)) ! c print*,'Pente entre ktop et kupper' ! c $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i)) ! c ENDDO ! c DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN dp_deltomg(i, k) = dp_deltomg(i, kupper(i)) omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i)) END IF END DO END DO ! cc nrlmd ! c DO i=1,klon ! c print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1) ! c END DO ! cc ! -- Compute wake average vertical velocity omgbw DO k = 1, klev + 1 DO i = 1, klon IF (wk_adv(i)) THEN omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k) END IF END DO END DO ! -- and its vertical gradient dp_omgbw DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k)) END IF END DO END DO ! -- Upstream coefficients for omgb velocity ! -- (alpha_up(k) is the coefficient of the value at level k) ! -- (1-alpha_up(k) is the coefficient of the value at level k-1) DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN alpha_up(i, k) = 0. IF (omgb(i,k)>0.) alpha_up(i, k) = 1. END IF END DO END DO ! Matrix expressing [The,deltatw] from [Th1,Th2] DO i = 1, klon IF (wk_adv(i)) THEN rre1(i) = 1. - sigmaw(i) rre2(i) = sigmaw(i) END IF END DO rrd1 = -1. rrd2 = 1. ! -- Get [Th1,Th2], dth and [q1,q2] DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN dth(i, k) = deltatw(i, k)/ppi(i, k) th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd d_th1(i, 1) = 0. d_th2(i, 1) = 0. d_dth(i, 1) = 0. d_q1(i, 1) = 0. d_q2(i, 1) = 0. d_dq(i, 1) = 0. END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN d_th1(i, k) = th1(i, k-1) - th1(i, k) d_th2(i, k) = th2(i, k-1) - th2(i, k) d_dth(i, k) = dth(i, k-1) - dth(i, k) d_q1(i, k) = q1(i, k-1) - q1(i, k) d_q2(i, k) = q2(i, k-1) - q2(i, k) d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN omgbdth(i, 1) = 0. omgbdq(i, 1) = 0. END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN ! loop on interfaces omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k)) omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k)) END IF END DO END DO ! ----------------------------------------------------------------- DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN ! ----------------------------------------------------------------- ! Compute redistribution (advective) term d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* & (rrd1*omg(i,k)*sigmaw(i)*d_th1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( & i))*d_th2(i,k+1)-(1.-alpha_up(i,k))*omgbdth(i,k)-alpha_up(i,k+1)* & omgbdth(i,k+1))*ppi(i, k) ! print*,'d_deltatw=',d_deltatw(i,k) d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* & (rrd1*omg(i,k)*sigmaw(i)*d_q1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( & i))*d_q2(i,k+1)-(1.-alpha_up(i,k))*omgbdq(i,k)-alpha_up(i,k+1)* & omgbdq(i,k+1)) ! print*,'d_deltaqw=',d_deltaqw(i,k) ! and increment large scale tendencies ! C ! ----------------------------------------------------------------- d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, & k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/(ph(i,k)-ph(i, & k+1)) & ! cc nrlmd $ ! -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k) -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/(ph(i, & k)-ph(i,k+1)) & ! cc )*ppi(i, k) d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, & k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/(ph(i,k)-ph(i, & k+1)) & ! cc nrlmd $ ! -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k) -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i, & k+1))/(ph(i,k)-ph(i,k+1)) & ! cc ) ! cc nrlmd ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, & k)/(ph(i,k)-ph(i,k+1))))*ppi(i, k) d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, & k)/(ph(i,k)-ph(i,k+1)))) END IF ! cc END DO END DO ! ------------------------------------------------------------------ ! Increment state variables DO k = 1, klev DO i = 1, klon ! cc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN IF (wk_adv(i) .AND. k<=kupper(i)) THEN ! cc ! Coefficient de répartition crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ & (ph(i,kupper(i))-ph(i,1)) crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-ph(i & ,kupper(i))) ! Reintroduce compensating subsidence term. ! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)) ! . /(1-sigmaw) ! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)) ! . /(1-sigmaw) ! dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k)) ! . /(1-sigmaw) ! dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k)) ! . /(1-sigmaw) dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i))) dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i))) ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k) !jyg< !! !!--------------------------------------------------------------- !! The change of delta_T due to PBL (vertical diffusion plus thermal plumes) !! is accounted for by the PBL and the Thermals schemes. It is now set to zero !! within the Wake scheme. !!--------------------------------------------------------------- dtPBL(i,k) = 0. dqPBL(i,k) = 0. ! !! dtPBL(i,k)=wdtPBL(i,k) - udtPBL(i,k) !! dqPBL(i,k)=wdqPBL(i,k) - udqPBL(i,k) ! !! dtpbl(i, k) = (wdtpbl(i,k)/sigmaw(i)-udtpbl(i,k)/(1.-sigmaw(i))) !! dqpbl(i, k) = (wdqpbl(i,k)/sigmaw(i)-udqpbl(i,k)/(1.-sigmaw(i))) ! print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k) !>jyg ! ! cc nrlmd Prise en compte du taux de mortalité ! cc Définitions de entr, detr detr(i, k) = 0. entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + & sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k) spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i) ! cc spread(i,k) = ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/ ! cc $ sigmaw(i) ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU ! Jingmei ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k), ! & Tgw(i,k),deltatw(i,k) d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)*deltatw(i, k)* & dtimesub ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k) ff(i) = d_deltatw(i, k)/dtimesub ! Sans GW ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k)) ! GW formule 1 ! deltatw(k) = deltatw(k)+dtimesub* ! $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k)) ! GW formule 2 IF (dtimesub*tgw(i,k)<1.E-10) THEN d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc ! $ ! -spread(i,k)*deltatw(i,k) -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( & i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc -tgw(i,k)*deltatw(i,k)) ELSE d_deltatw(i, k) = 1/tgw(i, k)*(1-exp(-dtimesub*tgw(i, & k)))*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc $ ! -spread(i,k)*deltatw(i,k) -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( & i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc -tgw(i,k)*deltatw(i,k)) END IF dth(i, k) = deltatw(i, k)/ppi(i, k) gg(i) = d_deltaqw(i, k)/dtimesub d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k)+dqpbl(i,k) & ! cc $ ! -spread(i,k)*deltaqw(i,k)) -entr(i,k)*deltaqw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw(i)+detr( & i,k))*deltaqw(i,k)/(1.-sigmaw(i))) ! cc ! cc nrlmd ! cc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k) ! cc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k) ! cc END IF END DO END DO ! Scale tendencies so that water vapour remains positive in w and x. CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, & d_deltaqw, sigmaw, d_sigmaw, alpha) ! cc nrlmd ! c print*,'alpha' ! c do i=1,klon ! c print*,alpha(i) ! c end do ! cc DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN d_te(i, k) = alpha(i)*d_te(i, k) d_qe(i, k) = alpha(i)*d_qe(i, k) d_deltatw(i, k) = alpha(i)*d_deltatw(i, k) d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k) d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN d_sigmaw(i) = alpha(i)*d_sigmaw(i) END IF END DO ! Update large scale variables and wake variables ! IM 060208 manque DO i + remplace DO k=1,kupper(i) ! IM 060208 DO k = 1,kupper(i) DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN dtls(i, k) = dtls(i, k) + d_te(i, k) dqls(i, k) = dqls(i, k) + d_qe(i, k) ! cc nrlmd d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k) d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k) ! cc END IF END DO END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN te(i, k) = te0(i, k) + dtls(i, k) qe(i, k) = qe0(i, k) + dqls(i, k) the(i, k) = te(i, k)/ppi(i, k) deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k) deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k) dth(i, k) = deltatw(i, k)/ppi(i, k) ! c print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k) ! c $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN sigmaw(i) = sigmaw(i) + d_sigmaw(i) END IF END DO ! Determine Ptop from buoyancy integral ! --------------------------------------- ! - 1/ Pressure of the level where dth changes sign. DO i = 1, klon IF (wk_adv(i)) THEN ptop_provis(i) = ph(i, 1) END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. & dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO ! - 2/ dth integral DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd sum_dth(i) = 0. dthmin(i) = -delta_t_min z(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) dthmin(i) = amin1(dthmin(i), dth(i,k)) END IF END IF END DO END DO ! - 3/ height of triangle with area= sum_dth and base = dthmin DO i = 1, klon IF (wk_adv(i)) THEN hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) hw(i) = amax1(hwmin, hw(i)) END IF END DO ! - 4/ now, get Ptop DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd ktop(i) = 0 z(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i)) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i) ktop(i) = k END IF END IF END DO END DO ! 4.5/Correct ktop and ptop DO i = 1, klon IF (wk_adv(i)) THEN ptop_new(i) = ptop(i) END IF END DO DO k = klev, 2, -1 DO i = 1, klon ! IM v3JYG; IF (k .GE. ktop(i) IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN ptop(i) = ptop_new(i) END IF END DO DO k = klev, 1, -1 DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (ph(i,k+1)<ptop(i)) ktop(i) = k END IF END DO END DO ! 5/ Set deltatw & deltaqw to 0 above kupper DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k>=kupper(i)) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. END IF END DO END DO ! -------------Cstar computation--------------------------------- DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd sum_thu(i) = 0. sum_tu(i) = 0. sum_qu(i) = 0. sum_thvu(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_rho(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thu(i) = 0. av_tu(i) = 0. av_qu(i) = 0. av_thvu(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_rho(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END IF END DO ! Integrals (and wake top level number) ! -------------------------------------- ! Initialize sum_thvu to 1st level virt. pot. temp. DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd z(i) = 1. dz(i) = 1. sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i) sum_dth(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) END IF END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd hw0(i) = z(i) END IF END DO ! - WAPE and mean forcing computation ! --------------------------------------- ! --------------------------------------- ! Means DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd av_thu(i) = sum_thu(i)/hw0(i) av_tu(i) = sum_tu(i)/hw0(i) av_qu(i) = sum_qu(i)/hw0(i) av_thvu(i) = sum_thvu(i)/hw0(i) av_dth(i) = sum_dth(i)/hw0(i) av_dq(i) = sum_dq(i)/hw0(i) av_rho(i) = sum_rho(i)/hw0(i) av_dtdwn(i) = sum_dtdwn(i)/hw0(i) av_dqdwn(i) = sum_dqdwn(i)/hw0(i) wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* & av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i) END IF END DO ! Filter out bad wakes DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (wape(i)<0.) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. END IF END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (wape(i)<0.) THEN wape(i) = 0. cstar(i) = 0. hw(i) = hwmin sigmaw(i) = max(sigmad, sigd_con(i)) fip(i) = 0. gwake(i) = .FALSE. ELSE cstar(i) = stark*sqrt(2.*wape(i)) gwake(i) = .TRUE. END IF END IF END DO END DO ! end sub-timestep loop ! ----------------------------------------------------------------- ! Get back to tendencies per second DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN ! cc dtls(i, k) = dtls(i, k)/dtime dqls(i, k) = dqls(i, k)/dtime d_deltatw2(i, k) = d_deltatw2(i, k)/dtime d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime ! c print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k) ! c $ ,death_rate(i)*sigmaw(i) END IF END DO END DO ! ---------------------------------------------------------- ! Determine wake final state; recompute wape, cstar, ktop; ! filter out bad wakes. ! ---------------------------------------------------------- ! 2.1 - Undisturbed area and Wake integrals ! --------------------------------------------------------- DO i = 1, klon ! cc nrlmd if (wk_adv(i)) then !!! nrlmd IF (ok_qx_qw(i)) THEN ! cc z(i) = 0. sum_thu(i) = 0. sum_tu(i) = 0. sum_qu(i) = 0. sum_thvu(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_rho(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thu(i) = 0. av_tu(i) = 0. av_qu(i) = 0. av_thvu(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_rho(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END IF END DO ! Potential temperatures and humidity ! ---------------------------------------------------------- DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc rho(i, k) = p(i, k)/(rd*te(i,k)) IF (k==1) THEN rhoh(i, k) = ph(i, k)/(rd*te(i,k)) zhh(i, k) = 0 ELSE rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1))) zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1) END IF the(i, k) = te(i, k)/ppi(i, k) thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i) qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k))) dth(i, k) = deltatw(i, k)/ppi(i, k) END IF END DO END DO ! Integrals (and wake top level number) ! ----------------------------------------------------------- ! Initialize sum_thvu to 1st level virt. pot. temp. DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc z(i) = 1. dz(i) = 1. sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i) sum_dth(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) END IF END IF END DO END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc hw0(i) = z(i) END IF END DO ! - WAPE and mean forcing computation ! ------------------------------------------------------------- ! Means DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc av_thu(i) = sum_thu(i)/hw0(i) av_tu(i) = sum_tu(i)/hw0(i) av_qu(i) = sum_qu(i)/hw0(i) av_thvu(i) = sum_thvu(i)/hw0(i) av_dth(i) = sum_dth(i)/hw0(i) av_dq(i) = sum_dq(i)/hw0(i) av_rho(i) = sum_rho(i)/hw0(i) av_dtdwn(i) = sum_dtdwn(i)/hw0(i) av_dqdwn(i) = sum_dqdwn(i)/hw0(i) wape2(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* & av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i) END IF END DO ! Prognostic variable update ! ------------------------------------------------------------ ! Filter out bad wakes DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN ! cc deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. END IF END DO END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc IF (wape2(i)<0.) THEN wape2(i) = 0. cstar2(i) = 0. hw(i) = hwmin sigmaw(i) = amax1(sigmad, sigd_con(i)) fip(i) = 0. gwake(i) = .FALSE. ELSE IF (prt_level>=10) PRINT *, 'wape2>0' cstar2(i) = stark*sqrt(2.*wape2(i)) gwake(i) = .TRUE. END IF END IF END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc ktopw(i) = ktop(i) END IF END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc IF (ktopw(i)>0 .AND. gwake(i)) THEN ! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer) ! cc heff = 600. ! Utilisation de la hauteur hw ! c heff = 0.7*hw heff(i) = hw(i) fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* & sqrt(sigmaw(i)*wdens(i)*3.14) fip(i) = alpk*fip(i) ! jyg2 ELSE fip(i) = 0. END IF END IF END DO ! Limitation de sigmaw ! cc nrlmd ! DO i=1,klon ! IF (OK_qx_qw(i)) THEN ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max ! ENDIF ! ENDDO ! cc DO k = 1, klev DO i = 1, klon ! cc nrlmd On maintient désormais constant sigmaw en régime ! permanent ! cc IF ((sigmaw(i).GT.sigmaw_max).or. IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. & .NOT. ok_qx_qw(i)) THEN ! cc dtls(i, k) = 0. dqls(i, k) = 0. deltatw(i, k) = 0. deltaqw(i, k) = 0. END IF END DO END DO ! cc nrlmd On maintient désormais constant sigmaw en régime permanent DO i = 1, klon IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. & .NOT. ok_qx_qw(i)) THEN wape(i) = 0. cstar(i) = 0. !!jyg Outside subroutine "Wake" hw and sigmaw are zero when there are no wakes !! hw(i) = hwmin !jyg !! sigmaw(i) = sigmad !jyg hw(i) = 0. !jyg sigmaw(i) = 0. !jyg fip(i) = 0. ELSE wape(i) = wape2(i) cstar(i) = cstar2(i) END IF ! c print*,'wape wape2 ktopw OK_qx_qw =', ! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i) END DO RETURN END SUBROUTINE wake SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, & d_deltaqw, sigmaw, d_sigmaw, alpha) ! ------------------------------------------------------ ! Dtermination du coefficient alpha tel que les tendances ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent ! a une humidite positive dans la zone (x) et dans la zone (w). ! ------------------------------------------------------ IMPLICIT NONE ! Input REAL qe(nlon, nl), d_qe(nlon, nl) REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl) REAL sigmaw(nlon), d_sigmaw(nlon) LOGICAL wk_adv(nlon) INTEGER nl, nlon ! Output REAL alpha(nlon) ! Internal variables REAL zeta(nlon, nl) REAL alpha1(nlon) REAL x, a, b, c, discrim REAL epsilon ! DATA epsilon/1.e-15/ INTEGER i,k DO k = 1, nl DO i = 1, nlon IF (wk_adv(i)) THEN IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN zeta(i, k) = 0. ELSE zeta(i, k) = 1. END IF END IF END DO DO i = 1, nlon IF (wk_adv(i)) THEN x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + & (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i)*(deltaqw(i,k)+ & d_deltaqw(i,k)) a = -d_sigmaw(i)*d_deltaqw(i, k) b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - & deltaqw(i, k)*d_sigmaw(i) c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon discrim = b*b - 4.*a*c ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap alpha1(i) = 1. ELSE IF (x>=0.) THEN alpha1(i) = 1. ELSE IF (a>0.) THEN alpha1(i) = 0.9*min((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim & ))/(2.*a)) ELSE IF (a==0.) THEN alpha1(i) = 0.9*(-c/b) ELSE ! print*,'a,b,c discrim',a,b,c discrim alpha1(i) = 0.9*max((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim & ))/(2.*a)) END IF END IF END IF alpha(i) = min(alpha(i), alpha1(i)) END IF END DO END DO RETURN END SUBROUTINE wake_vec_modulation
110  ! LL : distance entre 2 poches
111 
112  ! -------------------------------------------------------------------------
113  ! Déclaration de variables
114  ! -------------------------------------------------------------------------
115 
116  include "YOMCST.h"
117  include "cvthermo.h"
118 
119  ! Arguments en entree
120  ! --------------------
121 
122  REAL, DIMENSION (klon, klev), INTENT(IN) :: p, pi
123  REAL, DIMENSION (klon, klev+1), INTENT(IN) :: ph, omgb
124  REAL, INTENT(IN) :: dtime
125  REAL, DIMENSION (klon, klev), INTENT(IN) :: te0, qe0
126  REAL, DIMENSION (klon, klev), INTENT(IN) :: dtdwn, dqdwn
127  REAL, DIMENSION (klon, klev), INTENT(IN) :: wdtpbl, wdqpbl, udtpbl, udqpbl ! UNUSED
128  REAL, DIMENSION (klon, klev), INTENT(IN) :: amdwn, amup
129  REAL, DIMENSION (klon, klev), INTENT(IN) :: dta, dqa
130  REAL, DIMENSION (klon), INTENT(IN) :: sigd_con
131 
132  !
133  ! Input/Output
134  REAL, DIMENSION (klon, klev), INTENT(INOUT) :: deltatw, deltaqw
135  REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw
136 
137  ! Sorties
138  ! --------
139 
140  REAL, DIMENSION (klon, klev), INTENT(OUT) :: dth
141  REAL, DIMENSION (klon, klev), INTENT(OUT) :: tu, qu
142  REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtls, dqls
143  REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtke, dqke
144  REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtpbl, dqpbl
145  REAL, DIMENSION (klon, klev), INTENT(OUT) :: spread
146  REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltatw2, d_deltaqw2
147  REAL, DIMENSION (klon, klev+1), INTENT(OUT) :: omgbdth, omg
148  REAL, DIMENSION (klon, klev), INTENT(OUT) :: dp_omgb, dp_deltomg
149  REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_gw
150  REAL, DIMENSION (klon), INTENT(OUT) :: hw, wape, fip, gfl, cstar
151  REAL, DIMENSION (klon), INTENT(OUT) :: wdens
152  INTEGER, DIMENSION (klon), INTENT(OUT) :: ktopw
153 
154  ! Variables internes
155  ! -------------------
156 
157  ! Variables à fixer
158  REAL alon
159  LOGICAL, SAVE :: first = .true.
160  !$OMP THREADPRIVATE(first)
161  REAL, SAVE :: stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol
162  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
163  REAL delta_t_min
164  INTEGER nsub
165  REAL dtimesub
166  REAL sigmad, hwmin, wapecut
167  REAL :: sigmaw_max
168  REAL :: dens_rate
169  REAL wdens0
170  ! IM 080208
171  LOGICAL, DIMENSION (klon) :: gwake
172 
173  ! Variables de sauvegarde
174  REAL, DIMENSION (klon, klev) :: deltatw0
175  REAL, DIMENSION (klon, klev) :: deltaqw0
176  REAL, DIMENSION (klon, klev) :: te, qe
177  REAL, DIMENSION (klon) :: sigmaw0, sigmaw1
178 
179  ! Variables pour les GW
180  REAL, DIMENSION (klon) :: ll
181  REAL, DIMENSION (klon, klev) :: n2
182  REAL, DIMENSION (klon, klev) :: cgw
183  REAL, DIMENSION (klon, klev) :: tgw
184 
185  ! Variables liées au calcul de hw
186  REAL, DIMENSION (klon) :: ptop_provis, ptop, ptop_new
187  REAL, DIMENSION (klon) :: sum_dth
188  REAL, DIMENSION (klon) :: dthmin
189  REAL, DIMENSION (klon) :: z, dz, hw0
190  INTEGER, DIMENSION (klon) :: ktop, kupper
191 
192  ! Sub-timestep tendencies and related variables
193  REAL d_deltatw(klon, klev), d_deltaqw(klon, klev)
194  REAL d_te(klon, klev), d_qe(klon, klev)
195  REAL d_sigmaw(klon), alpha(klon)
196  REAL q0_min(klon), q1_min(klon)
197  LOGICAL wk_adv(klon), ok_qx_qw(klon)
198  REAL epsilon
199  DATA epsilon/1.e-15/
200 
201  ! Autres variables internes
202  INTEGER isubstep, k, i
203 
204  REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu
205  REAL, DIMENSION (klon) :: sum_dq, sum_rho
206  REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn
207  REAL, DIMENSION (klon) :: av_thu, av_tu, av_qu, av_thvu
208  REAL, DIMENSION (klon) :: av_dth, av_dq, av_rho
209  REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn
210 
211  REAL, DIMENSION (klon, klev) :: rho, rhow
212  REAL, DIMENSION (klon, klev+1) :: rhoh
213  REAL, DIMENSION (klon, klev) :: rhow_moyen
214  REAL, DIMENSION (klon, klev) :: zh
215  REAL, DIMENSION (klon, klev+1) :: zhh
216  REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2
217 
218  REAL, DIMENSION (klon, klev) :: the, thu
219 
220  ! REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
221 
222  REAL, DIMENSION (klon, klev+1) :: omgbw
223  REAL, DIMENSION (klon) :: pupper
224  REAL, DIMENSION (klon) :: omgtop
225  REAL, DIMENSION (klon, klev) :: dp_omgbw
226  REAL, DIMENSION (klon) :: ztop, dztop
227  REAL, DIMENSION (klon, klev) :: alpha_up
228 
229  REAL, DIMENSION (klon) :: rre1, rre2
230  REAL :: rrd1, rrd2
231  REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2
232  REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth
233  REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq
234  REAL, DIMENSION (klon, klev) :: omgbdq
235 
236  REAL, DIMENSION (klon) :: ff, gg
237  REAL, DIMENSION (klon) :: wape2, cstar2, heff
238 
239  REAL, DIMENSION (klon, klev) :: crep
240 
241  REAL, DIMENSION (klon, klev) :: ppi
242 
243  ! cc nrlmd
244  REAL, DIMENSION (klon) :: death_rate, nat_rate
245  REAL, DIMENSION (klon, klev) :: entr
246  REAL, DIMENSION (klon, klev) :: detr
247 
248  ! -------------------------------------------------------------------------
249  ! Initialisations
250  ! -------------------------------------------------------------------------
251 
252  ! print*, 'wake initialisations'
253 
254  ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
255  ! -------------------------------------------------------------------------
256 
257  DATA wapecut, sigmad, hwmin/5., .02, 10./
258  ! cc nrlmd
259  DATA sigmaw_max/0.4/
260  DATA dens_rate/0.1/
261  ! cc
262  ! Longueur de maille (en m)
263  ! -------------------------------------------------------------------------
264 
265  ! ALON = 3.e5
266  alon = 1.e6
267 
268 
269  ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
270 
271  ! coefgw : Coefficient pour les ondes de gravité ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) ! wdens : Densité de poche froide par maille ! ------------------------------------------------------------------------- ! cc nrlmd coefgw=10 ! coefgw=1 ! wdens0 = 1.0/(alon**2) ! cc nrlmd wdens = 1.0/(alon**2) ! cc nrlmd stark = 0.50 ! CRtest ! cc nrlmd alpk=0.1 ! alpk = 1.0 ! alpk = 0.5 ! alpk = 0.05 if (first) then stark = 0.33 alpk = 0.25 wdens_ref = 8.E-12 coefgw = 4. crep_upper = 0.9 crep_sol = 1.0 ! cc nrlmd Lecture du fichier wake_param.data !$OMP MASTER OPEN (99, FILE='wake_param.data', STATUS='old', FORM='formatted', ERR=9999) READ (99, *, END=9998) stark READ (99, *, END=9998) alpk READ (99, *, END=9998) wdens_ref READ (99, *, END=9998) coefgw 9998 CONTINUE CLOSE (99) 9999 CONTINUE !$OMP END MASTER CALL bcast(stark) CALL bcast(alpk) CALL bcast(wdens_ref) CALL bcast(coefgw) first=.false. endif ! Initialisation de toutes des densites a wdens_ref. ! Les densites peuvent evoluer si les poches debordent ! (voir au tout debut de la boucle sur les substeps) wdens = wdens_ref ! print*,'stark',stark ! print*,'alpk',alpk ! print*,'wdens',wdens ! print*,'coefgw',coefgw ! cc ! Minimum value for |T_wake - T_undist|. Used for wake top definition ! ------------------------------------------------------------------------- delta_t_min = 0.2 ! 1. - Save initial values and initialize tendencies ! -------------------------------------------------- DO k = 1, klev DO i = 1, klon ppi(i, k) = pi(i, k) deltatw0(i, k) = deltatw(i, k) deltaqw0(i, k) = deltaqw(i, k) te(i, k) = te0(i, k) qe(i, k) = qe0(i, k) dtls(i, k) = 0. dqls(i, k) = 0. d_deltat_gw(i, k) = 0. d_te(i, k) = 0. d_qe(i, k) = 0. d_deltatw(i, k) = 0. d_deltaqw(i, k) = 0. ! IM 060508 beg d_deltatw2(i, k) = 0. d_deltaqw2(i, k) = 0. ! IM 060508 end END DO END DO ! sigmaw1=sigmaw ! IF (sigd_con.GT.sigmaw1) THEN ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con ! ENDIF DO i = 1, klon ! c sigmaw(i) = amax1(sigmaw(i),sigd_con(i)) sigmaw(i) = amax1(sigmaw(i), sigmad) sigmaw(i) = amin1(sigmaw(i), 0.99) sigmaw0(i) = sigmaw(i) wape(i) = 0. wape2(i) = 0. d_sigmaw(i) = 0. ktopw(i) = 0 END DO ! 2. - Prognostic part ! -------------------- ! 2.1 - Undisturbed area and Wake integrals ! --------------------------------------------------------- DO i = 1, klon z(i) = 0. ktop(i) = 0 kupper(i) = 0 sum_thu(i) = 0. sum_tu(i) = 0. sum_qu(i) = 0. sum_thvu(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_rho(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thu(i) = 0. av_tu(i) = 0. av_qu(i) = 0. av_thvu(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_rho(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END DO ! Distance between wakes DO i = 1, klon ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i)) END DO ! Potential temperatures and humidity ! ---------------------------------------------------------- DO k = 1, klev DO i = 1, klon ! write(*,*)'wake 1',i,k,rd,te(i,k) rho(i, k) = p(i, k)/(rd*te(i,k)) ! write(*,*)'wake 2',rho(i,k) IF (k==1) THEN ! write(*,*)'wake 3',i,k,rd,te(i,k) rhoh(i, k) = ph(i, k)/(rd*te(i,k)) ! write(*,*)'wake 4',i,k,rd,te(i,k) zhh(i, k) = 0 ELSE ! write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1)) rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1))) ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1) zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1) END IF ! write(*,*)'wake 7',ppi(i,k) the(i, k) = te(i, k)/ppi(i, k) thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i) qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k))) rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k))) dth(i, k) = deltatw(i, k)/ppi(i, k) END DO END DO DO k = 1, klev - 1 DO i = 1, klon IF (k==1) THEN n2(i, k) = 0 ELSE n2(i, k) = amax1(0., -rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i, & k-1))/(p(i,k+1)-p(i,k-1))) END IF zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2 cgw(i, k) = sqrt(n2(i,k))*zh(i, k) tgw(i, k) = coefgw*cgw(i, k)/ll(i) END DO END DO DO i = 1, klon n2(i, klev) = 0 zh(i, klev) = 0 cgw(i, klev) = 0 tgw(i, klev) = 0 END DO ! Calcul de la masse volumique moyenne de la colonne (bdlmd) ! ----------------------------------------------------------------- DO k = 1, klev DO i = 1, klon epaisseur1(i, k) = 0. epaisseur2(i, k) = 0. END DO END DO DO i = 1, klon epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1. epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1. rhow_moyen(i, 1) = rhow(i, 1) END DO DO k = 2, klev DO i = 1, klon epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1. epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k) rhow_moyen(i, k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+rhow(i,k)* & epaisseur1(i,k))/epaisseur2(i, k) END DO END DO ! Choose an integration bound well above wake top ! ----------------------------------------------------------------- ! Pupper = 50000. ! melting level ! Pupper = 60000. ! Pupper = 80000. ! essais pour case_e DO i = 1, klon pupper(i) = 0.6*ph(i, 1) pupper(i) = max(pupper(i), 45000.) ! cc Pupper(i) = 60000. END DO ! Determine Wake top pressure (Ptop) from buoyancy integral ! -------------------------------------------------------- ! -1/ Pressure of the level where dth becomes less than delta_t_min. DO i = 1, klon ptop_provis(i) = ph(i, 1) END DO DO k = 2, klev DO i = 1, klon ! IM v3JYG; ptop_provis(i).LT. ph(i,1) IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. & ptop_provis(i)==ph(i,1)) THEN ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO ! -2/ dth integral DO i = 1, klon sum_dth(i) = 0. dthmin(i) = -delta_t_min z(i) = 0. END DO DO k = 1, klev DO i = 1, klon dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) dthmin(i) = amin1(dthmin(i), dth(i,k)) END IF END DO END DO ! -3/ height of triangle with area= sum_dth and base = dthmin DO i = 1, klon hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) hw0(i) = amax1(hwmin, hw0(i)) END DO ! -4/ now, get Ptop DO i = 1, klon z(i) = 0. ptop(i) = ph(i, 1) END DO DO k = 1, klev DO i = 1, klon dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i)) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i) END IF END DO END DO ! -5/ Determination de ktop et kupper DO k = klev, 1, -1 DO i = 1, klon IF (ph(i,k+1)<ptop(i)) ktop(i) = k IF (ph(i,k+1)<pupper(i)) kupper(i) = k END DO END DO ! On evite kupper = 1 et kupper = klev DO i = 1, klon kupper(i) = max(kupper(i), 2) kupper(i) = min(kupper(i), klev-1) END DO ! -6/ Correct ktop and ptop DO i = 1, klon ptop_new(i) = ptop(i) END DO DO k = klev, 2, -1 DO i = 1, klon IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO DO i = 1, klon ptop(i) = ptop_new(i) END DO DO k = klev, 1, -1 DO i = 1, klon IF (ph(i,k+1)<ptop(i)) ktop(i) = k END DO END DO ! -5/ Set deltatw & deltaqw to 0 above kupper DO k = 1, klev DO i = 1, klon IF (k>=kupper(i)) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. END IF END DO END DO ! Vertical gradient of LS omega DO k = 1, klev DO i = 1, klon IF (k<=kupper(i)) THEN dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k)) END IF END DO END DO ! Integrals (and wake top level number) ! -------------------------------------- ! Initialize sum_thvu to 1st level virt. pot. temp. DO i = 1, klon z(i) = 1. dz(i) = 1. sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i) sum_dth(i) = 0. END DO DO k = 1, klev DO i = 1, klon dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) END IF END DO END DO DO i = 1, klon hw0(i) = z(i) END DO ! 2.1 - WAPE and mean forcing computation ! --------------------------------------- ! --------------------------------------- ! Means DO i = 1, klon av_thu(i) = sum_thu(i)/hw0(i) av_tu(i) = sum_tu(i)/hw0(i) av_qu(i) = sum_qu(i)/hw0(i) av_thvu(i) = sum_thvu(i)/hw0(i) ! av_thve = sum_thve/hw0 av_dth(i) = sum_dth(i)/hw0(i) av_dq(i) = sum_dq(i)/hw0(i) av_rho(i) = sum_rho(i)/hw0(i) av_dtdwn(i) = sum_dtdwn(i)/hw0(i) av_dqdwn(i) = sum_dqdwn(i)/hw0(i) wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i & )+av_dth(i)*av_dq(i)))/av_thvu(i) END DO ! 2.2 Prognostic variable update ! ------------------------------ ! Filter out bad wakes DO k = 1, klev DO i = 1, klon IF (wape(i)<0.) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. END IF END DO END DO DO i = 1, klon IF (wape(i)<0.) THEN wape(i) = 0. cstar(i) = 0. hw(i) = hwmin sigmaw(i) = amax1(sigmad, sigd_con(i)) fip(i) = 0. gwake(i) = .FALSE. ELSE cstar(i) = stark*sqrt(2.*wape(i)) gwake(i) = .TRUE. END IF END DO ! Check qx and qw positivity ! -------------------------- DO i = 1, klon q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), (qe(i, & 1)+(1.-sigmaw(i))*deltaqw(i,1))) END DO DO k = 2, klev DO i = 1, klon q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), (qe(i, & k)+(1.-sigmaw(i))*deltaqw(i,k))) IF (q1_min(i)<=q0_min(i)) THEN q0_min(i) = q1_min(i) END IF END DO END DO DO i = 1, klon ok_qx_qw(i) = q0_min(i) >= 0. alpha(i) = 1. END DO ! C ----------------------------------------------------------------- ! Sub-time-stepping ! ----------------- nsub = 10 dtimesub = dtime/nsub ! ------------------------------------------------------------ DO isubstep = 1, nsub ! ------------------------------------------------------------ ! wk_adv is the logical flag enabling wake evolution in the time advance ! loop DO i = 1, klon wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1. END DO ! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement ! négatif de ktop à kupper -------- ! cc On calcule pour cela une densité wdens0 pour laquelle on ! aurait un entrainement nul --- DO i = 1, klon ! c print *,' isubstep,wk_adv(i),cstar(i),wape(i) ', ! c $ isubstep,wk_adv(i),cstar(i),wape(i) IF (wk_adv(i) .AND. cstar(i)>0.01) THEN omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + & rg*amup(i, kupper(i)+1)/(1.-sigmaw(i)) wdens0 = (sigmaw(i)/(4.*3.14))*((1.-sigmaw(i))*omg(i,kupper(i)+1)/(( & ph(i,1)-pupper(i))*cstar(i)))**(2) IF (wdens(i)<=wdens0*1.1) THEN wdens(i) = wdens0 END IF ! c print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i) ! c $ ,ph(i,1)-pupper(i)', ! c $ omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i) ! c $ ,ph(i,1)-pupper(i) END IF END DO ! cc nrlmd DO i = 1, klon IF (wk_adv(i)) THEN gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i)) sigmaw(i) = amin1(sigmaw(i), sigmaw_max) END IF END DO DO i = 1, klon IF (wk_adv(i)) THEN ! cc nrlmd Introduction du taux de mortalité des poches et ! test sur sigmaw_max=0.4 ! cc d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub IF (sigmaw(i)>=sigmaw_max) THEN death_rate(i) = gfl(i)*cstar(i)/sigmaw(i) ELSE death_rate(i) = 0. END IF d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* & dtimesub ! $ - nat_rate(i)*sigmaw(i)*dtimesub ! c print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), ! c $ death_rate(i),ktop(i),kupper(i)', ! c $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), ! c $ death_rate(i),ktop(i),kupper(i) ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub ! sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!! ! wdens = wdens0/(10.*sigmaw) ! sigmaw =max(sigmaw,sigd_con) ! sigmaw =max(sigmaw,sigmad) END IF END DO ! calcul de la difference de vitesse verticale poche - zone non perturbee ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on ! definit ! IM 060208 au niveau k=1..? DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd dp_deltomg(i, k) = 0. END IF END DO END DO DO k = 1, klev + 1 DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd omg(i, k) = 0. END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN z(i) = 0. omg(i, 1) = 0. dp_deltomg(i, 1) = -(gfl(i)*cstar(i))/(sigmaw(i)*(1-sigmaw(i))) END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=ktop(i)) THEN dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg) z(i) = z(i) + dz(i) dp_deltomg(i, k) = dp_deltomg(i, 1) omg(i, k) = dp_deltomg(i, 1)*z(i) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg) ztop(i) = z(i) + dztop(i) omgtop(i) = dp_deltomg(i, 1)*ztop(i) END IF END DO ! ----------------- ! From m/s to Pa/s ! ----------------- DO i = 1, klon IF (wk_adv(i)) THEN omgtop(i) = -rho(i, ktop(i))*rg*omgtop(i) dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1)) END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=ktop(i)) THEN omg(i, k) = -rho(i, k)*rg*omg(i, k) dp_deltomg(i, k) = dp_deltomg(i, 1) END IF END DO END DO ! raccordement lineaire de omg de ptop a pupper DO i = 1, klon IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + & rg*amup(i, kupper(i)+1)/(1.-sigmaw(i)) dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ & (ptop(i)-pupper(i)) END IF END DO ! c DO i=1,klon ! c print*,'Pente entre 0 et kupper (référence)' ! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1)) ! c print*,'Pente entre ktop et kupper' ! c $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i)) ! c ENDDO ! c DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN dp_deltomg(i, k) = dp_deltomg(i, kupper(i)) omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i)) END IF END DO END DO ! cc nrlmd ! c DO i=1,klon ! c print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1) ! c END DO ! cc ! -- Compute wake average vertical velocity omgbw DO k = 1, klev + 1 DO i = 1, klon IF (wk_adv(i)) THEN omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k) END IF END DO END DO ! -- and its vertical gradient dp_omgbw DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k)) END IF END DO END DO ! -- Upstream coefficients for omgb velocity ! -- (alpha_up(k) is the coefficient of the value at level k) ! -- (1-alpha_up(k) is the coefficient of the value at level k-1) DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN alpha_up(i, k) = 0. IF (omgb(i,k)>0.) alpha_up(i, k) = 1. END IF END DO END DO ! Matrix expressing [The,deltatw] from [Th1,Th2] DO i = 1, klon IF (wk_adv(i)) THEN rre1(i) = 1. - sigmaw(i) rre2(i) = sigmaw(i) END IF END DO rrd1 = -1. rrd2 = 1. ! -- Get [Th1,Th2], dth and [q1,q2] DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN dth(i, k) = deltatw(i, k)/ppi(i, k) th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd d_th1(i, 1) = 0. d_th2(i, 1) = 0. d_dth(i, 1) = 0. d_q1(i, 1) = 0. d_q2(i, 1) = 0. d_dq(i, 1) = 0. END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN d_th1(i, k) = th1(i, k-1) - th1(i, k) d_th2(i, k) = th2(i, k-1) - th2(i, k) d_dth(i, k) = dth(i, k-1) - dth(i, k) d_q1(i, k) = q1(i, k-1) - q1(i, k) d_q2(i, k) = q2(i, k-1) - q2(i, k) d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN omgbdth(i, 1) = 0. omgbdq(i, 1) = 0. END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN ! loop on interfaces omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k)) omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k)) END IF END DO END DO ! ----------------------------------------------------------------- DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN ! ----------------------------------------------------------------- ! Compute redistribution (advective) term d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* & (rrd1*omg(i,k)*sigmaw(i)*d_th1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( & i))*d_th2(i,k+1)-(1.-alpha_up(i,k))*omgbdth(i,k)-alpha_up(i,k+1)* & omgbdth(i,k+1))*ppi(i, k) ! print*,'d_deltatw=',d_deltatw(i,k) d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* & (rrd1*omg(i,k)*sigmaw(i)*d_q1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( & i))*d_q2(i,k+1)-(1.-alpha_up(i,k))*omgbdq(i,k)-alpha_up(i,k+1)* & omgbdq(i,k+1)) ! print*,'d_deltaqw=',d_deltaqw(i,k) ! and increment large scale tendencies ! C ! ----------------------------------------------------------------- d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, & k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/(ph(i,k)-ph(i, & k+1)) & ! cc nrlmd $ ! -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k) -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/(ph(i, & k)-ph(i,k+1)) & ! cc )*ppi(i, k) d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, & k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/(ph(i,k)-ph(i, & k+1)) & ! cc nrlmd $ ! -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k) -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i, & k+1))/(ph(i,k)-ph(i,k+1)) & ! cc ) ! cc nrlmd ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, & k)/(ph(i,k)-ph(i,k+1))))*ppi(i, k) d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, & k)/(ph(i,k)-ph(i,k+1)))) END IF ! cc END DO END DO ! ------------------------------------------------------------------ ! Increment state variables DO k = 1, klev DO i = 1, klon ! cc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN IF (wk_adv(i) .AND. k<=kupper(i)) THEN ! cc ! Coefficient de répartition crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ & (ph(i,kupper(i))-ph(i,1)) crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-ph(i & ,kupper(i))) ! Reintroduce compensating subsidence term. ! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)) ! . /(1-sigmaw) ! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)) ! . /(1-sigmaw) ! dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k)) ! . /(1-sigmaw) ! dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k)) ! . /(1-sigmaw) dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i))) dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i))) ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k) !jyg< !! !!--------------------------------------------------------------- !! The change of delta_T due to PBL (vertical diffusion plus thermal plumes) !! is accounted for by the PBL and the Thermals schemes. It is now set to zero !! within the Wake scheme. !!--------------------------------------------------------------- dtPBL(i,k) = 0. dqPBL(i,k) = 0. ! !! dtPBL(i,k)=wdtPBL(i,k) - udtPBL(i,k) !! dqPBL(i,k)=wdqPBL(i,k) - udqPBL(i,k) ! !! dtpbl(i, k) = (wdtpbl(i,k)/sigmaw(i)-udtpbl(i,k)/(1.-sigmaw(i))) !! dqpbl(i, k) = (wdqpbl(i,k)/sigmaw(i)-udqpbl(i,k)/(1.-sigmaw(i))) ! print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k) !>jyg ! ! cc nrlmd Prise en compte du taux de mortalité ! cc Définitions de entr, detr detr(i, k) = 0. entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + & sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k) spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i) ! cc spread(i,k) = ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/ ! cc $ sigmaw(i) ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU ! Jingmei ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k), ! & Tgw(i,k),deltatw(i,k) d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)*deltatw(i, k)* & dtimesub ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k) ff(i) = d_deltatw(i, k)/dtimesub ! Sans GW ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k)) ! GW formule 1 ! deltatw(k) = deltatw(k)+dtimesub* ! $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k)) ! GW formule 2 IF (dtimesub*tgw(i,k)<1.E-10) THEN d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc ! $ ! -spread(i,k)*deltatw(i,k) -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( & i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc -tgw(i,k)*deltatw(i,k)) ELSE d_deltatw(i, k) = 1/tgw(i, k)*(1-exp(-dtimesub*tgw(i, & k)))*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc $ ! -spread(i,k)*deltatw(i,k) -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( & i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc -tgw(i,k)*deltatw(i,k)) END IF dth(i, k) = deltatw(i, k)/ppi(i, k) gg(i) = d_deltaqw(i, k)/dtimesub d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k)+dqpbl(i,k) & ! cc $ ! -spread(i,k)*deltaqw(i,k)) -entr(i,k)*deltaqw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw(i)+detr( & i,k))*deltaqw(i,k)/(1.-sigmaw(i))) ! cc ! cc nrlmd ! cc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k) ! cc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k) ! cc END IF END DO END DO ! Scale tendencies so that water vapour remains positive in w and x. CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, & d_deltaqw, sigmaw, d_sigmaw, alpha) ! cc nrlmd ! c print*,'alpha' ! c do i=1,klon ! c print*,alpha(i) ! c end do ! cc DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN d_te(i, k) = alpha(i)*d_te(i, k) d_qe(i, k) = alpha(i)*d_qe(i, k) d_deltatw(i, k) = alpha(i)*d_deltatw(i, k) d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k) d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN d_sigmaw(i) = alpha(i)*d_sigmaw(i) END IF END DO ! Update large scale variables and wake variables ! IM 060208 manque DO i + remplace DO k=1,kupper(i) ! IM 060208 DO k = 1,kupper(i) DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN dtls(i, k) = dtls(i, k) + d_te(i, k) dqls(i, k) = dqls(i, k) + d_qe(i, k) ! cc nrlmd d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k) d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k) ! cc END IF END DO END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN te(i, k) = te0(i, k) + dtls(i, k) qe(i, k) = qe0(i, k) + dqls(i, k) the(i, k) = te(i, k)/ppi(i, k) deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k) deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k) dth(i, k) = deltatw(i, k)/ppi(i, k) ! c print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k) ! c $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN sigmaw(i) = sigmaw(i) + d_sigmaw(i) END IF END DO ! Determine Ptop from buoyancy integral ! --------------------------------------- ! - 1/ Pressure of the level where dth changes sign. DO i = 1, klon IF (wk_adv(i)) THEN ptop_provis(i) = ph(i, 1) END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. & dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO ! - 2/ dth integral DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd sum_dth(i) = 0. dthmin(i) = -delta_t_min z(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) dthmin(i) = amin1(dthmin(i), dth(i,k)) END IF END IF END DO END DO ! - 3/ height of triangle with area= sum_dth and base = dthmin DO i = 1, klon IF (wk_adv(i)) THEN hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) hw(i) = amax1(hwmin, hw(i)) END IF END DO ! - 4/ now, get Ptop DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd ktop(i) = 0 z(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i)) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i) ktop(i) = k END IF END IF END DO END DO ! 4.5/Correct ktop and ptop DO i = 1, klon IF (wk_adv(i)) THEN ptop_new(i) = ptop(i) END IF END DO DO k = klev, 2, -1 DO i = 1, klon ! IM v3JYG; IF (k .GE. ktop(i) IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN ptop(i) = ptop_new(i) END IF END DO DO k = klev, 1, -1 DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (ph(i,k+1)<ptop(i)) ktop(i) = k END IF END DO END DO ! 5/ Set deltatw & deltaqw to 0 above kupper DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k>=kupper(i)) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. END IF END DO END DO ! -------------Cstar computation--------------------------------- DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd sum_thu(i) = 0. sum_tu(i) = 0. sum_qu(i) = 0. sum_thvu(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_rho(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thu(i) = 0. av_tu(i) = 0. av_qu(i) = 0. av_thvu(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_rho(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END IF END DO ! Integrals (and wake top level number) ! -------------------------------------- ! Initialize sum_thvu to 1st level virt. pot. temp. DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd z(i) = 1. dz(i) = 1. sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i) sum_dth(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) END IF END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd hw0(i) = z(i) END IF END DO ! - WAPE and mean forcing computation ! --------------------------------------- ! --------------------------------------- ! Means DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd av_thu(i) = sum_thu(i)/hw0(i) av_tu(i) = sum_tu(i)/hw0(i) av_qu(i) = sum_qu(i)/hw0(i) av_thvu(i) = sum_thvu(i)/hw0(i) av_dth(i) = sum_dth(i)/hw0(i) av_dq(i) = sum_dq(i)/hw0(i) av_rho(i) = sum_rho(i)/hw0(i) av_dtdwn(i) = sum_dtdwn(i)/hw0(i) av_dqdwn(i) = sum_dqdwn(i)/hw0(i) wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* & av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i) END IF END DO ! Filter out bad wakes DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (wape(i)<0.) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. END IF END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (wape(i)<0.) THEN wape(i) = 0. cstar(i) = 0. hw(i) = hwmin sigmaw(i) = max(sigmad, sigd_con(i)) fip(i) = 0. gwake(i) = .FALSE. ELSE cstar(i) = stark*sqrt(2.*wape(i)) gwake(i) = .TRUE. END IF END IF END DO END DO ! end sub-timestep loop ! ----------------------------------------------------------------- ! Get back to tendencies per second DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN ! cc dtls(i, k) = dtls(i, k)/dtime dqls(i, k) = dqls(i, k)/dtime d_deltatw2(i, k) = d_deltatw2(i, k)/dtime d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime ! c print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k) ! c $ ,death_rate(i)*sigmaw(i) END IF END DO END DO ! ---------------------------------------------------------- ! Determine wake final state; recompute wape, cstar, ktop; ! filter out bad wakes. ! ---------------------------------------------------------- ! 2.1 - Undisturbed area and Wake integrals ! --------------------------------------------------------- DO i = 1, klon ! cc nrlmd if (wk_adv(i)) then !!! nrlmd IF (ok_qx_qw(i)) THEN ! cc z(i) = 0. sum_thu(i) = 0. sum_tu(i) = 0. sum_qu(i) = 0. sum_thvu(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_rho(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thu(i) = 0. av_tu(i) = 0. av_qu(i) = 0. av_thvu(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_rho(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END IF END DO ! Potential temperatures and humidity ! ---------------------------------------------------------- DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc rho(i, k) = p(i, k)/(rd*te(i,k)) IF (k==1) THEN rhoh(i, k) = ph(i, k)/(rd*te(i,k)) zhh(i, k) = 0 ELSE rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1))) zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1) END IF the(i, k) = te(i, k)/ppi(i, k) thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i) qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k))) dth(i, k) = deltatw(i, k)/ppi(i, k) END IF END DO END DO ! Integrals (and wake top level number) ! ----------------------------------------------------------- ! Initialize sum_thvu to 1st level virt. pot. temp. DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc z(i) = 1. dz(i) = 1. sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i) sum_dth(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) END IF END IF END DO END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc hw0(i) = z(i) END IF END DO ! - WAPE and mean forcing computation ! ------------------------------------------------------------- ! Means DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc av_thu(i) = sum_thu(i)/hw0(i) av_tu(i) = sum_tu(i)/hw0(i) av_qu(i) = sum_qu(i)/hw0(i) av_thvu(i) = sum_thvu(i)/hw0(i) av_dth(i) = sum_dth(i)/hw0(i) av_dq(i) = sum_dq(i)/hw0(i) av_rho(i) = sum_rho(i)/hw0(i) av_dtdwn(i) = sum_dtdwn(i)/hw0(i) av_dqdwn(i) = sum_dqdwn(i)/hw0(i) wape2(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* & av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i) END IF END DO ! Prognostic variable update ! ------------------------------------------------------------ ! Filter out bad wakes DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN ! cc deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. END IF END DO END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc IF (wape2(i)<0.) THEN wape2(i) = 0. cstar2(i) = 0. hw(i) = hwmin sigmaw(i) = amax1(sigmad, sigd_con(i)) fip(i) = 0. gwake(i) = .FALSE. ELSE IF (prt_level>=10) PRINT *, 'wape2>0' cstar2(i) = stark*sqrt(2.*wape2(i)) gwake(i) = .TRUE. END IF END IF END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc ktopw(i) = ktop(i) END IF END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc IF (ktopw(i)>0 .AND. gwake(i)) THEN ! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer) ! cc heff = 600. ! Utilisation de la hauteur hw ! c heff = 0.7*hw heff(i) = hw(i) fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* & sqrt(sigmaw(i)*wdens(i)*3.14) fip(i) = alpk*fip(i) ! jyg2 ELSE fip(i) = 0. END IF END IF END DO ! Limitation de sigmaw ! cc nrlmd ! DO i=1,klon ! IF (OK_qx_qw(i)) THEN ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max ! ENDIF ! ENDDO ! cc DO k = 1, klev DO i = 1, klon ! cc nrlmd On maintient désormais constant sigmaw en régime ! permanent ! cc IF ((sigmaw(i).GT.sigmaw_max).or. IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. & .NOT. ok_qx_qw(i)) THEN ! cc dtls(i, k) = 0. dqls(i, k) = 0. deltatw(i, k) = 0. deltaqw(i, k) = 0. END IF END DO END DO ! cc nrlmd On maintient désormais constant sigmaw en régime permanent DO i = 1, klon IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. & .NOT. ok_qx_qw(i)) THEN wape(i) = 0. cstar(i) = 0. !!jyg Outside subroutine "Wake" hw and sigmaw are zero when there are no wakes !! hw(i) = hwmin !jyg !! sigmaw(i) = sigmad !jyg hw(i) = 0. !jyg sigmaw(i) = 0. !jyg fip(i) = 0. ELSE wape(i) = wape2(i) cstar(i) = cstar2(i) END IF ! c print*,'wape wape2 ktopw OK_qx_qw =', ! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i) END DO RETURN END SUBROUTINE wake SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, & d_deltaqw, sigmaw, d_sigmaw, alpha) ! ------------------------------------------------------ ! Dtermination du coefficient alpha tel que les tendances ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent ! a une humidite positive dans la zone (x) et dans la zone (w). ! ------------------------------------------------------ IMPLICIT NONE ! Input REAL qe(nlon, nl), d_qe(nlon, nl) REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl) REAL sigmaw(nlon), d_sigmaw(nlon) LOGICAL wk_adv(nlon) INTEGER nl, nlon ! Output REAL alpha(nlon) ! Internal variables REAL zeta(nlon, nl) REAL alpha1(nlon) REAL x, a, b, c, discrim REAL epsilon ! DATA epsilon/1.e-15/ INTEGER i,k DO k = 1, nl DO i = 1, nlon IF (wk_adv(i)) THEN IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN zeta(i, k) = 0. ELSE zeta(i, k) = 1. END IF END IF END DO DO i = 1, nlon IF (wk_adv(i)) THEN x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + & (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i)*(deltaqw(i,k)+ & d_deltaqw(i,k)) a = -d_sigmaw(i)*d_deltaqw(i, k) b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - & deltaqw(i, k)*d_sigmaw(i) c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon discrim = b*b - 4.*a*c ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap alpha1(i) = 1. ELSE IF (x>=0.) THEN alpha1(i) = 1. ELSE IF (a>0.) THEN alpha1(i) = 0.9*min((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim & ))/(2.*a)) ELSE IF (a==0.) THEN alpha1(i) = 0.9*(-c/b) ELSE ! print*,'a,b,c discrim',a,b,c discrim alpha1(i) = 0.9*max((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim & ))/(2.*a)) END IF END IF END IF alpha(i) = min(alpha(i), alpha1(i)) END IF END DO END DO RETURN END SUBROUTINE wake_vec_modulation
272  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
273  ! wdens : Densité de poche froide par maille
274  ! -------------------------------------------------------------------------
275 
276  ! cc nrlmd coefgw=10
277  ! coefgw=1
278  ! wdens0 = 1.0/(alon**2)
279  ! cc nrlmd wdens = 1.0/(alon**2)
280  ! cc nrlmd stark = 0.50
281  ! CRtest
282  ! cc nrlmd alpk=0.1
283  ! alpk = 1.0
284  ! alpk = 0.5
285  ! alpk = 0.05
286 
287  if (first) then
288  stark = 0.33
289  alpk = 0.25
290  wdens_ref = 8.e-12
291  coefgw = 4.
292  crep_upper = 0.9
293  crep_sol = 1.0
294 
295  ! cc nrlmd Lecture du fichier wake_param.data
296  !$OMP MASTER
297  OPEN (99, file='wake_param.data', status='old', form='formatted', err=9999)
298  READ (99, *, end=9998) stark
299  READ (99, *, end=9998) alpk
300  READ (99, *, end=9998) wdens_ref
301  READ (99, *, end=9998) coefgw
302 9998 CONTINUE
303  CLOSE (99)
304 9999 CONTINUE
305  !$OMP END MASTER
306  CALL bcast(stark)
307  CALL bcast(alpk)
308  CALL bcast(wdens_ref)
309  CALL bcast(coefgw)
310 
311  first=.false.
312  endif
313 
314  ! Initialisation de toutes des densites a wdens_ref.
315  ! Les densites peuvent evoluer si les poches debordent
316  ! (voir au tout debut de la boucle sur les substeps)
317  wdens = wdens_ref
318 
319  ! print*,'stark',stark
320  ! print*,'alpk',alpk
321  ! print*,'wdens',wdens
322  ! print*,'coefgw',coefgw
323  ! cc
324  ! Minimum value for |T_wake - T_undist|. Used for wake top definition
325  ! -------------------------------------------------------------------------
326 
327  delta_t_min = 0.2
328 
329  ! 1. - Save initial values and initialize tendencies
330  ! --------------------------------------------------
331 
332  DO k = 1, klev
333  DO i = 1, klon
334  ppi(i, k) = pi(i, k)
335  deltatw0(i, k) = deltatw(i, k)
336  deltaqw0(i, k) = deltaqw(i, k)
337  te(i, k) = te0(i, k)
338  qe(i, k) = qe0(i, k)
339  dtls(i, k) = 0.
340  dqls(i, k) = 0.
341  d_deltat_gw(i, k) = 0.
342  d_te(i, k) = 0.
343  d_qe(i, k) = 0.
344  d_deltatw(i, k) = 0.
345  d_deltaqw(i, k) = 0.
346  ! IM 060508 beg
347  d_deltatw2(i, k) = 0.
348  d_deltaqw2(i, k) = 0.
349  ! IM 060508 end
350  END DO
351  END DO
352  ! sigmaw1=sigmaw
353  ! IF (sigd_con.GT.sigmaw1) THEN
354  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
355  ! ENDIF
356  DO i = 1, klon
357  ! c sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
358  sigmaw(i) = amax1(sigmaw(i), sigmad)
359  sigmaw(i) = amin1(sigmaw(i), 0.99)
360  sigmaw0(i) = sigmaw(i)
361  wape(i) = 0.
362  wape2(i) = 0.
363  d_sigmaw(i) = 0.
364  ktopw(i) = 0
365  END DO
366 
367 
368  ! 2. - Prognostic part
369  ! --------------------
370 
371 
372  ! 2.1 - Undisturbed area and Wake integrals
373  ! ---------------------------------------------------------
374 
375  DO i = 1, klon
376  z(i) = 0.
377  ktop(i) = 0
378  kupper(i) = 0
379  sum_thu(i) = 0.
380  sum_tu(i) = 0.
381  sum_qu(i) = 0.
382  sum_thvu(i) = 0.
383  sum_dth(i) = 0.
384  sum_dq(i) = 0.
385  sum_rho(i) = 0.
386  sum_dtdwn(i) = 0.
387  sum_dqdwn(i) = 0.
388 
389  av_thu(i) = 0.
390  av_tu(i) = 0.
391  av_qu(i) = 0.
392  av_thvu(i) = 0.
393  av_dth(i) = 0.
394  av_dq(i) = 0.
395  av_rho(i) = 0.
396  av_dtdwn(i) = 0.
397  av_dqdwn(i) = 0.
398  END DO
399 
400  ! Distance between wakes
401  DO i = 1, klon
402  ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
403  END DO
404  ! Potential temperatures and humidity
405  ! ----------------------------------------------------------
406  DO k = 1, klev
407  DO i = 1, klon
408  ! write(*,*)'wake 1',i,k,rd,te(i,k)
409  rho(i, k) = p(i, k)/(rd*te(i,k))
410  ! write(*,*)'wake 2',rho(i,k)
411  IF (k==1) THEN
412  ! write(*,*)'wake 3',i,k,rd,te(i,k)
413  rhoh(i, k) = ph(i, k)/(rd*te(i,k))
414  ! write(*,*)'wake 4',i,k,rd,te(i,k)
415  zhh(i, k) = 0
416  ELSE
417  ! write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1))
418  rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
419  ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
420  zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
421  END IF
422  ! write(*,*)'wake 7',ppi(i,k)
423  the(i, k) = te(i, k)/ppi(i, k)
424  thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
425  tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
426  qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
427  ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
428  rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
429  dth(i, k) = deltatw(i, k)/ppi(i, k)
430  END DO
431  END DO
432 
433  DO k = 1, klev - 1
434  DO i = 1, klon
435  IF (k==1) THEN
436  n2(i, k) = 0
437  ELSE
438  n2(i, k) = amax1(0., -rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i, &
439  k-1))/(p(i,k+1)-p(i,k-1)))
440  END IF
441  zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
442 
443  cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
444  tgw(i, k) = coefgw*cgw(i, k)/ll(i)
445  END DO
446  END DO
447 
448  DO i = 1, klon
449  n2(i, klev) = 0
450  zh(i, klev) = 0
451  cgw(i, klev) = 0
452  tgw(i, klev) = 0
453  END DO
454 
455  ! Calcul de la masse volumique moyenne de la colonne (bdlmd)
456  ! -----------------------------------------------------------------
457 
458  DO k = 1, klev
459  DO i = 1, klon
460  epaisseur1(i, k) = 0.
461  epaisseur2(i, k) = 0.
462  END DO
463  END DO
464 
465  DO i = 1, klon
466  epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
467  epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
468  rhow_moyen(i, 1) = rhow(i, 1)
469  END DO
470 
471  DO k = 2, klev
472  DO i = 1, klon
473  epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1.
474  epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k)
475  rhow_moyen(i, k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+rhow(i,k)* &
476  epaisseur1(i,k))/epaisseur2(i, k)
477  END DO
478  END DO
479 
480 
481  ! Choose an integration bound well above wake top
482  ! -----------------------------------------------------------------
483 
484  ! Pupper = 50000. ! melting level
485  ! Pupper = 60000.
486  ! Pupper = 80000. ! essais pour case_e
487  DO i = 1, klon
488  pupper(i) = 0.6*ph(i, 1)
489  pupper(i) = max(pupper(i), 45000.)
490  ! cc Pupper(i) = 60000.
491  END DO
492 
493 
494  ! Determine Wake top pressure (Ptop) from buoyancy integral
495  ! --------------------------------------------------------
496 
497  ! -1/ Pressure of the level where dth becomes less than delta_t_min.
498 
499  DO i = 1, klon
500  ptop_provis(i) = ph(i, 1)
501  END DO
502  DO k = 2, klev
503  DO i = 1, klon
504 
505  ! IM v3JYG; ptop_provis(i).LT. ph(i,1)
506 
507  IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. &
508  ptop_provis(i)==ph(i,1)) THEN
509  ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
510  k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
511  END IF
512  END DO
513  END DO
514 
515  ! -2/ dth integral
516 
517  DO i = 1, klon
518  sum_dth(i) = 0.
519  dthmin(i) = -delta_t_min
520  z(i) = 0.
521  END DO
522 
523  DO k = 1, klev
524  DO i = 1, klon
525  dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
526  IF (dz(i)>0) THEN
527  z(i) = z(i) + dz(i)
528  sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
529  dthmin(i) = amin1(dthmin(i), dth(i,k))
530  END IF
531  END DO
532  END DO
533 
534  ! -3/ height of triangle with area= sum_dth and base = dthmin
535 
536  DO i = 1, klon
537  hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
538  hw0(i) = amax1(hwmin, hw0(i))
539  END DO
540 
541  ! -4/ now, get Ptop
542 
543  DO i = 1, klon
544  z(i) = 0.
545  ptop(i) = ph(i, 1)
546  END DO
547 
548  DO k = 1, klev
549  DO i = 1, klon
550  dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i))
551  IF (dz(i)>0) THEN
552  z(i) = z(i) + dz(i)
553  ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i)
554  END IF
555  END DO
556  END DO
557 
558 
559  ! -5/ Determination de ktop et kupper
560 
561  DO k = klev, 1, -1
562  DO i = 1, klon
563  IF (ph(i,k+1)<ptop(i)) ktop(i) = k
564  IF (ph(i,k+1)<pupper(i)) kupper(i) = k
565  END DO
566  END DO
567 
568  ! On evite kupper = 1 et kupper = klev
569  DO i = 1, klon
570  kupper(i) = max(kupper(i), 2)
571  kupper(i) = min(kupper(i), klev-1)
572  END DO
573 
574 
575  ! -6/ Correct ktop and ptop
576 
577  DO i = 1, klon
578  ptop_new(i) = ptop(i)
579  END DO
580  DO k = klev, 2, -1
581  DO i = 1, klon
582  IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
583  dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
584  ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
585  k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
586  END IF
587  END DO
588  END DO
589 
590  DO i = 1, klon
591  ptop(i) = ptop_new(i)
592  END DO
593 
594  DO k = klev, 1, -1
595  DO i = 1, klon
596  IF (ph(i,k+1)<ptop(i)) ktop(i) = k
597  END DO
598  END DO
599 
600  ! -5/ Set deltatw & deltaqw to 0 above kupper
601 
602  DO k = 1, klev
603  DO i = 1, klon
604  IF (k>=kupper(i)) THEN
605  deltatw(i, k) = 0.
606  deltaqw(i, k) = 0.
607  END IF
608  END DO
609  END DO
610 
611 
612  ! Vertical gradient of LS omega
613 
614  DO k = 1, klev
615  DO i = 1, klon
616  IF (k<=kupper(i)) THEN
617  dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
618  END IF
619  END DO
620  END DO
621 
622  ! Integrals (and wake top level number)
623  ! --------------------------------------
624 
625  ! Initialize sum_thvu to 1st level virt. pot. temp.
626 
627  DO i = 1, klon
628  z(i) = 1.
629  dz(i) = 1.
630  sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i)
631  sum_dth(i) = 0.
632  END DO
633 
634  DO k = 1, klev
635  DO i = 1, klon
636  dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
637  IF (dz(i)>0) THEN
638  z(i) = z(i) + dz(i)
639  sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
640  sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
641  sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
642  sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i)
643  sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
644  sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
645  sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
646  sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
647  sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
648  END IF
649  END DO
650  END DO
651 
652  DO i = 1, klon
653  hw0(i) = z(i)
654  END DO
655 
656 
657  ! 2.1 - WAPE and mean forcing computation
658  ! ---------------------------------------
659 
660  ! ---------------------------------------
661 
662  ! Means
663 
664  DO i = 1, klon
665  av_thu(i) = sum_thu(i)/hw0(i)
666  av_tu(i) = sum_tu(i)/hw0(i)
667  av_qu(i) = sum_qu(i)/hw0(i)
668  av_thvu(i) = sum_thvu(i)/hw0(i)
669  ! av_thve = sum_thve/hw0
670  av_dth(i) = sum_dth(i)/hw0(i)
671  av_dq(i) = sum_dq(i)/hw0(i)
672  av_rho(i) = sum_rho(i)/hw0(i)
673  av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
674  av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
675 
676  wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i &
677  )+av_dth(i)*av_dq(i)))/av_thvu(i)
678  END DO
679 
680  ! 2.2 Prognostic variable update
681  ! ------------------------------
682 
683  ! Filter out bad wakes
684 
685  DO k = 1, klev
686  DO i = 1, klon
687  IF (wape(i)<0.) THEN
688  deltatw(i, k) = 0.
689  deltaqw(i, k) = 0.
690  dth(i, k) = 0.
691  END IF
692  END DO
693  END DO
694 
695  DO i = 1, klon
696  IF (wape(i)<0.) THEN
697  wape(i) = 0.
698  cstar(i) = 0.
699  hw(i) = hwmin
700  sigmaw(i) = amax1(sigmad, sigd_con(i))
701  fip(i) = 0.
702  gwake(i) = .false.
703  ELSE
704  cstar(i) = stark*sqrt(2.*wape(i))
705  gwake(i) = .true.
706  END IF
707  END DO
708 
709 
710  ! Check qx and qw positivity
711  ! --------------------------
712  DO i = 1, klon
713  q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), (qe(i, &
714  1)+(1.-sigmaw(i))*deltaqw(i,1)))
715  END DO
716  DO k = 2, klev
717  DO i = 1, klon
718  q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), (qe(i, &
719  k)+(1.-sigmaw(i))*deltaqw(i,k)))
720  IF (q1_min(i)<=q0_min(i)) THEN
721  q0_min(i) = q1_min(i)
722  END IF
723  END DO
724  END DO
725 
726  DO i = 1, klon
727  ok_qx_qw(i) = q0_min(i) >= 0.
728  alpha(i) = 1.
729  END DO
730 
731  ! C -----------------------------------------------------------------
732  ! Sub-time-stepping
733  ! -----------------
734 
735  nsub = 10
736  dtimesub = dtime/nsub
737 
738  ! ------------------------------------------------------------
739  DO isubstep = 1, nsub
740  ! ------------------------------------------------------------
741 
742  ! wk_adv is the logical flag enabling wake evolution in the time advance
743  ! loop
744  DO i = 1, klon
745  wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
746  END DO
747 
748  ! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement
749  ! négatif de ktop à kupper --------
750  ! cc On calcule pour cela une densité wdens0 pour laquelle on
751  ! aurait un entrainement nul ---
752  DO i = 1, klon
753  ! c print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
754  ! c $ isubstep,wk_adv(i),cstar(i),wape(i)
755  IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
756  omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
757  rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
758  wdens0 = (sigmaw(i)/(4.*3.14))*((1.-sigmaw(i))*omg(i,kupper(i)+1)/(( &
759  ph(i,1)-pupper(i))*cstar(i)))**(2)
760  IF (wdens(i)<=wdens0*1.1) THEN
761  wdens(i) = wdens0
762  END IF
763  ! c print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
764  ! c $ ,ph(i,1)-pupper(i)',
765  ! c $ omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
766  ! c $ ,ph(i,1)-pupper(i)
767  END IF
768  END DO
769 
770  ! cc nrlmd
771 
772  DO i = 1, klon
773  IF (wk_adv(i)) THEN
774  gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
775  sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
776  END IF
777  END DO
778  DO i = 1, klon
779  IF (wk_adv(i)) THEN
780  ! cc nrlmd Introduction du taux de mortalité des poches et
781  ! test sur sigmaw_max=0.4
782  ! cc d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
783  IF (sigmaw(i)>=sigmaw_max) THEN
784  death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
785  ELSE
786  death_rate(i) = 0.
787  END IF
788  d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* &
789  dtimesub
790  ! $ - nat_rate(i)*sigmaw(i)*dtimesub
791  ! c print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
792  ! c $ death_rate(i),ktop(i),kupper(i)',
793  ! c $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
794  ! c $ death_rate(i),ktop(i),kupper(i)
795 
796  ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
797  ! sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!!
798  ! wdens = wdens0/(10.*sigmaw)
799  ! sigmaw =max(sigmaw,sigd_con)
800  ! sigmaw =max(sigmaw,sigmad)
801  END IF
802  END DO
803 
804 
805  ! calcul de la difference de vitesse verticale poche - zone non perturbee
806  ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
807  ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on
808  ! definit
809  ! IM 060208 au niveau k=1..?
810  DO k = 1, klev
811  DO i = 1, klon
812  IF (wk_adv(i)) THEN !!! nrlmd
813  dp_deltomg(i, k) = 0.
814  END IF
815  END DO
816  END DO
817  DO k = 1, klev + 1
818  DO i = 1, klon
819  IF (wk_adv(i)) THEN !!! nrlmd
820  omg(i, k) = 0.
821  END IF
822  END DO
823  END DO
824 
825  DO i = 1, klon
826  IF (wk_adv(i)) THEN
827  z(i) = 0.
828  omg(i, 1) = 0.
829  dp_deltomg(i, 1) = -(gfl(i)*cstar(i))/(sigmaw(i)*(1-sigmaw(i)))
830  END IF
831  END DO
832 
833  DO k = 2, klev
834  DO i = 1, klon
835  IF (wk_adv(i) .AND. k<=ktop(i)) THEN
836  dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
837  z(i) = z(i) + dz(i)
838  dp_deltomg(i, k) = dp_deltomg(i, 1)
839  omg(i, k) = dp_deltomg(i, 1)*z(i)
840  END IF
841  END DO
842  END DO
843 
844  DO i = 1, klon
845  IF (wk_adv(i)) THEN
846  dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
847  ztop(i) = z(i) + dztop(i)
848  omgtop(i) = dp_deltomg(i, 1)*ztop(i)
849  END IF
850  END DO
851 
852  ! -----------------
853  ! From m/s to Pa/s
854  ! -----------------
855 
856  DO i = 1, klon
857  IF (wk_adv(i)) THEN
858  omgtop(i) = -rho(i, ktop(i))*rg*omgtop(i)
859  dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
860  END IF
861  END DO
862 
863  DO k = 1, klev
864  DO i = 1, klon
865  IF (wk_adv(i) .AND. k<=ktop(i)) THEN
866  omg(i, k) = -rho(i, k)*rg*omg(i, k)
867  dp_deltomg(i, k) = dp_deltomg(i, 1)
868  END IF
869  END DO
870  END DO
871 
872  ! raccordement lineaire de omg de ptop a pupper
873 
874  DO i = 1, klon
875  IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
876  omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
877  rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
878  dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
879  (ptop(i)-pupper(i))
880  END IF
881  END DO
882 
883  ! c DO i=1,klon
884  ! c print*,'Pente entre 0 et kupper (référence)'
885  ! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
886  ! c print*,'Pente entre ktop et kupper'
887  ! c $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
888  ! c ENDDO
889  ! c
890  DO k = 1, klev
891  DO i = 1, klon
892  IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
893  dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
894  omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
895  END IF
896  END DO
897  END DO
898  ! cc nrlmd
899  ! c DO i=1,klon
900  ! c print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
901  ! c END DO
902  ! cc
903 
904 
905  ! -- Compute wake average vertical velocity omgbw
906 
907 
908  DO k = 1, klev + 1
909  DO i = 1, klon
910  IF (wk_adv(i)) THEN
911  omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
912  END IF
913  END DO
914  END DO
915  ! -- and its vertical gradient dp_omgbw
916 
917  DO k = 1, klev
918  DO i = 1, klon
919  IF (wk_adv(i)) THEN
920  dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
921  END IF
922  END DO
923  END DO
924 
925  ! -- Upstream coefficients for omgb velocity
926  ! -- (alpha_up(k) is the coefficient of the value at level k)
927  ! -- (1-alpha_up(k) is the coefficient of the value at level k-1)
928  DO k = 1, klev
929  DO i = 1, klon
930  IF (wk_adv(i)) THEN
931  alpha_up(i, k) = 0.
932  IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
933  END IF
934  END DO
935  END DO
936 
937  ! Matrix expressing [The,deltatw] from [Th1,Th2]
938 
939  DO i = 1, klon
940  IF (wk_adv(i)) THEN
941  rre1(i) = 1. - sigmaw(i)
942  rre2(i) = sigmaw(i)
943  END IF
944  END DO
945  rrd1 = -1.
946  rrd2 = 1.
947 
948  ! -- Get [Th1,Th2], dth and [q1,q2]
949 
950  DO k = 1, klev
951  DO i = 1, klon
952  IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
953  dth(i, k) = deltatw(i, k)/ppi(i, k)
954  th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
955  th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
956  q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
957  q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
958  END IF
959  END DO
960  END DO
961 
962  DO i = 1, klon
963  IF (wk_adv(i)) THEN !!! nrlmd
964  d_th1(i, 1) = 0.
965  d_th2(i, 1) = 0.
966  d_dth(i, 1) = 0.
967  d_q1(i, 1) = 0.
968  d_q2(i, 1) = 0.
969  d_dq(i, 1) = 0.
970  END IF
971  END DO
972 
973  DO k = 2, klev
974  DO i = 1, klon
975  IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
976  d_th1(i, k) = th1(i, k-1) - th1(i, k)
977  d_th2(i, k) = th2(i, k-1) - th2(i, k)
978  d_dth(i, k) = dth(i, k-1) - dth(i, k)
979  d_q1(i, k) = q1(i, k-1) - q1(i, k)
980  d_q2(i, k) = q2(i, k-1) - q2(i, k)
981  d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
982  END IF
983  END DO
984  END DO
985 
986  DO i = 1, klon
987  IF (wk_adv(i)) THEN
988  omgbdth(i, 1) = 0.
989  omgbdq(i, 1) = 0.
990  END IF
991  END DO
992 
993  DO k = 2, klev
994  DO i = 1, klon
995  IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN ! loop on interfaces
996  omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
997  omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
998  END IF
999  END DO
1000  END DO
1001 
1002  ! -----------------------------------------------------------------
1003  DO k = 1, klev
1004  DO i = 1, klon
1005  IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
1006  ! -----------------------------------------------------------------
1007 
1008  ! Compute redistribution (advective) term
1009 
1010  d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1011  (rrd1*omg(i,k)*sigmaw(i)*d_th1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( &
1012  i))*d_th2(i,k+1)-(1.-alpha_up(i,k))*omgbdth(i,k)-alpha_up(i,k+1)* &
1013  omgbdth(i,k+1))*ppi(i, k)
1014  ! print*,'d_deltatw=',d_deltatw(i,k)
1015 
1016  d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1017  (rrd1*omg(i,k)*sigmaw(i)*d_q1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( &
1018  i))*d_q2(i,k+1)-(1.-alpha_up(i,k))*omgbdq(i,k)-alpha_up(i,k+1)* &
1019  omgbdq(i,k+1))
1020  ! print*,'d_deltaqw=',d_deltaqw(i,k)
1021 
1022  ! and increment large scale tendencies
1023 
1024 
1025 
1026 
1027  ! C
1028  ! -----------------------------------------------------------------
1029  d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, &
1030  k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/(ph(i,k)-ph(i, &
1031  k+1)) & ! cc nrlmd $
1032  ! -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k)
1033  -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/(ph(i, &
1034  k)-ph(i,k+1)) & ! cc
1035  )*ppi(i, k)
1036 
1037  d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, &
1038  k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/(ph(i,k)-ph(i, &
1039  k+1)) & ! cc nrlmd $
1040  ! -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k)
1041  -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i, &
1042  k+1))/(ph(i,k)-ph(i,k+1)) & ! cc
1043  )
1044  ! cc nrlmd
1045  ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
1046  d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, &
1047  k)/(ph(i,k)-ph(i,k+1))))*ppi(i, k)
1048 
1049  d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, &
1050  k)/(ph(i,k)-ph(i,k+1))))
1051 
1052  END IF
1053  ! cc
1054  END DO
1055  END DO
1056  ! ------------------------------------------------------------------
1057 
1058  ! Increment state variables
1059 
1060  DO k = 1, klev
1061  DO i = 1, klon
1062  ! cc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
1063  IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1064  ! cc
1065 
1066 
1067 
1068  ! Coefficient de répartition
1069 
1070  crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
1071  (ph(i,kupper(i))-ph(i,1))
1072  crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-ph(i &
1073  ,kupper(i)))
1074 
1075 
1076  ! Reintroduce compensating subsidence term.
1077 
1078  ! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
1079  ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
1080  ! . /(1-sigmaw)
1081  ! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
1082  ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
1083  ! . /(1-sigmaw)
1084 
1085  ! dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
1086  ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
1087  ! . /(1-sigmaw)
1088  ! dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
1089  ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
1090  ! . /(1-sigmaw)
1091 
1092  dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
1093  dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
1094  ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
1095 
1096 !jyg<
1097 !!
1098 !!---------------------------------------------------------------
1099 !! The change of delta_T due to PBL (vertical diffusion plus thermal plumes)
1100 !! is accounted for by the PBL and the Thermals schemes. It is now set to zero
1101 !! within the Wake scheme.
1102 !!---------------------------------------------------------------
1103  dtpbl(i,k) = 0.
1104  dqpbl(i,k) = 0.
1105 !
1106 !! dtPBL(i,k)=wdtPBL(i,k) - udtPBL(i,k)
1107 !! dqPBL(i,k)=wdqPBL(i,k) - udqPBL(i,k)
1108 !
1109 !! dtpbl(i, k) = (wdtpbl(i,k)/sigmaw(i)-udtpbl(i,k)/(1.-sigmaw(i)))
1110 !! dqpbl(i, k) = (wdqpbl(i,k)/sigmaw(i)-udqpbl(i,k)/(1.-sigmaw(i)))
1111  ! print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k)
1113 !
1114 
1115  ! cc nrlmd Prise en compte du taux de mortalité ! cc Définitions de entr, detr detr(i, k) = 0. entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + & sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k) spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i) ! cc spread(i,k) = ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/ ! cc $ sigmaw(i) ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU ! Jingmei ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k), ! & Tgw(i,k),deltatw(i,k) d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)*deltatw(i, k)* & dtimesub ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k) ff(i) = d_deltatw(i, k)/dtimesub ! Sans GW ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k)) ! GW formule 1 ! deltatw(k) = deltatw(k)+dtimesub* ! $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k)) ! GW formule 2 IF (dtimesub*tgw(i,k)<1.E-10) THEN d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc ! $ ! -spread(i,k)*deltatw(i,k) -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( & i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc -tgw(i,k)*deltatw(i,k)) ELSE d_deltatw(i, k) = 1/tgw(i, k)*(1-exp(-dtimesub*tgw(i, & k)))*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc $ ! -spread(i,k)*deltatw(i,k) -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( & i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc -tgw(i,k)*deltatw(i,k)) END IF dth(i, k) = deltatw(i, k)/ppi(i, k) gg(i) = d_deltaqw(i, k)/dtimesub d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k)+dqpbl(i,k) & ! cc $ ! -spread(i,k)*deltaqw(i,k)) -entr(i,k)*deltaqw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw(i)+detr( & i,k))*deltaqw(i,k)/(1.-sigmaw(i))) ! cc ! cc nrlmd ! cc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k) ! cc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k) ! cc END IF END DO END DO ! Scale tendencies so that water vapour remains positive in w and x. CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, & d_deltaqw, sigmaw, d_sigmaw, alpha) ! cc nrlmd ! c print*,'alpha' ! c do i=1,klon ! c print*,alpha(i) ! c end do ! cc DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN d_te(i, k) = alpha(i)*d_te(i, k) d_qe(i, k) = alpha(i)*d_qe(i, k) d_deltatw(i, k) = alpha(i)*d_deltatw(i, k) d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k) d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN d_sigmaw(i) = alpha(i)*d_sigmaw(i) END IF END DO ! Update large scale variables and wake variables ! IM 060208 manque DO i + remplace DO k=1,kupper(i) ! IM 060208 DO k = 1,kupper(i) DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN dtls(i, k) = dtls(i, k) + d_te(i, k) dqls(i, k) = dqls(i, k) + d_qe(i, k) ! cc nrlmd d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k) d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k) ! cc END IF END DO END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN te(i, k) = te0(i, k) + dtls(i, k) qe(i, k) = qe0(i, k) + dqls(i, k) the(i, k) = te(i, k)/ppi(i, k) deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k) deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k) dth(i, k) = deltatw(i, k)/ppi(i, k) ! c print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k) ! c $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN sigmaw(i) = sigmaw(i) + d_sigmaw(i) END IF END DO ! Determine Ptop from buoyancy integral ! --------------------------------------- ! - 1/ Pressure of the level where dth changes sign. DO i = 1, klon IF (wk_adv(i)) THEN ptop_provis(i) = ph(i, 1) END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. & dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO ! - 2/ dth integral DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd sum_dth(i) = 0. dthmin(i) = -delta_t_min z(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) dthmin(i) = amin1(dthmin(i), dth(i,k)) END IF END IF END DO END DO ! - 3/ height of triangle with area= sum_dth and base = dthmin DO i = 1, klon IF (wk_adv(i)) THEN hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) hw(i) = amax1(hwmin, hw(i)) END IF END DO ! - 4/ now, get Ptop DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd ktop(i) = 0 z(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i)) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i) ktop(i) = k END IF END IF END DO END DO ! 4.5/Correct ktop and ptop DO i = 1, klon IF (wk_adv(i)) THEN ptop_new(i) = ptop(i) END IF END DO DO k = klev, 2, -1 DO i = 1, klon ! IM v3JYG; IF (k .GE. ktop(i) IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN ptop(i) = ptop_new(i) END IF END DO DO k = klev, 1, -1 DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (ph(i,k+1)<ptop(i)) ktop(i) = k END IF END DO END DO ! 5/ Set deltatw & deltaqw to 0 above kupper DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k>=kupper(i)) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. END IF END DO END DO ! -------------Cstar computation--------------------------------- DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd sum_thu(i) = 0. sum_tu(i) = 0. sum_qu(i) = 0. sum_thvu(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_rho(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thu(i) = 0. av_tu(i) = 0. av_qu(i) = 0. av_thvu(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_rho(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END IF END DO ! Integrals (and wake top level number) ! -------------------------------------- ! Initialize sum_thvu to 1st level virt. pot. temp. DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd z(i) = 1. dz(i) = 1. sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i) sum_dth(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) END IF END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd hw0(i) = z(i) END IF END DO ! - WAPE and mean forcing computation ! --------------------------------------- ! --------------------------------------- ! Means DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd av_thu(i) = sum_thu(i)/hw0(i) av_tu(i) = sum_tu(i)/hw0(i) av_qu(i) = sum_qu(i)/hw0(i) av_thvu(i) = sum_thvu(i)/hw0(i) av_dth(i) = sum_dth(i)/hw0(i) av_dq(i) = sum_dq(i)/hw0(i) av_rho(i) = sum_rho(i)/hw0(i) av_dtdwn(i) = sum_dtdwn(i)/hw0(i) av_dqdwn(i) = sum_dqdwn(i)/hw0(i) wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* & av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i) END IF END DO ! Filter out bad wakes DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (wape(i)<0.) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. END IF END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (wape(i)<0.) THEN wape(i) = 0. cstar(i) = 0. hw(i) = hwmin sigmaw(i) = max(sigmad, sigd_con(i)) fip(i) = 0. gwake(i) = .FALSE. ELSE cstar(i) = stark*sqrt(2.*wape(i)) gwake(i) = .TRUE. END IF END IF END DO END DO ! end sub-timestep loop ! ----------------------------------------------------------------- ! Get back to tendencies per second DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN ! cc dtls(i, k) = dtls(i, k)/dtime dqls(i, k) = dqls(i, k)/dtime d_deltatw2(i, k) = d_deltatw2(i, k)/dtime d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime ! c print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k) ! c $ ,death_rate(i)*sigmaw(i) END IF END DO END DO ! ---------------------------------------------------------- ! Determine wake final state; recompute wape, cstar, ktop; ! filter out bad wakes. ! ---------------------------------------------------------- ! 2.1 - Undisturbed area and Wake integrals ! --------------------------------------------------------- DO i = 1, klon ! cc nrlmd if (wk_adv(i)) then !!! nrlmd IF (ok_qx_qw(i)) THEN ! cc z(i) = 0. sum_thu(i) = 0. sum_tu(i) = 0. sum_qu(i) = 0. sum_thvu(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_rho(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thu(i) = 0. av_tu(i) = 0. av_qu(i) = 0. av_thvu(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_rho(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END IF END DO ! Potential temperatures and humidity ! ---------------------------------------------------------- DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc rho(i, k) = p(i, k)/(rd*te(i,k)) IF (k==1) THEN rhoh(i, k) = ph(i, k)/(rd*te(i,k)) zhh(i, k) = 0 ELSE rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1))) zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1) END IF the(i, k) = te(i, k)/ppi(i, k) thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i) qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k))) dth(i, k) = deltatw(i, k)/ppi(i, k) END IF END DO END DO ! Integrals (and wake top level number) ! ----------------------------------------------------------- ! Initialize sum_thvu to 1st level virt. pot. temp. DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc z(i) = 1. dz(i) = 1. sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i) sum_dth(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) END IF END IF END DO END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc hw0(i) = z(i) END IF END DO ! - WAPE and mean forcing computation ! ------------------------------------------------------------- ! Means DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc av_thu(i) = sum_thu(i)/hw0(i) av_tu(i) = sum_tu(i)/hw0(i) av_qu(i) = sum_qu(i)/hw0(i) av_thvu(i) = sum_thvu(i)/hw0(i) av_dth(i) = sum_dth(i)/hw0(i) av_dq(i) = sum_dq(i)/hw0(i) av_rho(i) = sum_rho(i)/hw0(i) av_dtdwn(i) = sum_dtdwn(i)/hw0(i) av_dqdwn(i) = sum_dqdwn(i)/hw0(i) wape2(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* & av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i) END IF END DO ! Prognostic variable update ! ------------------------------------------------------------ ! Filter out bad wakes DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN ! cc deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. END IF END DO END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc IF (wape2(i)<0.) THEN wape2(i) = 0. cstar2(i) = 0. hw(i) = hwmin sigmaw(i) = amax1(sigmad, sigd_con(i)) fip(i) = 0. gwake(i) = .FALSE. ELSE IF (prt_level>=10) PRINT *, 'wape2>0' cstar2(i) = stark*sqrt(2.*wape2(i)) gwake(i) = .TRUE. END IF END IF END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc ktopw(i) = ktop(i) END IF END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc IF (ktopw(i)>0 .AND. gwake(i)) THEN ! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer) ! cc heff = 600. ! Utilisation de la hauteur hw ! c heff = 0.7*hw heff(i) = hw(i) fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* & sqrt(sigmaw(i)*wdens(i)*3.14) fip(i) = alpk*fip(i) ! jyg2 ELSE fip(i) = 0. END IF END IF END DO ! Limitation de sigmaw ! cc nrlmd ! DO i=1,klon ! IF (OK_qx_qw(i)) THEN ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max ! ENDIF ! ENDDO ! cc DO k = 1, klev DO i = 1, klon ! cc nrlmd On maintient désormais constant sigmaw en régime ! permanent ! cc IF ((sigmaw(i).GT.sigmaw_max).or. IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. & .NOT. ok_qx_qw(i)) THEN ! cc dtls(i, k) = 0. dqls(i, k) = 0. deltatw(i, k) = 0. deltaqw(i, k) = 0. END IF END DO END DO ! cc nrlmd On maintient désormais constant sigmaw en régime permanent DO i = 1, klon IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. & .NOT. ok_qx_qw(i)) THEN wape(i) = 0. cstar(i) = 0. !!jyg Outside subroutine "Wake" hw and sigmaw are zero when there are no wakes !! hw(i) = hwmin !jyg !! sigmaw(i) = sigmad !jyg hw(i) = 0. !jyg sigmaw(i) = 0. !jyg fip(i) = 0. ELSE wape(i) = wape2(i) cstar(i) = cstar2(i) END IF ! c print*,'wape wape2 ktopw OK_qx_qw =', ! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i) END DO RETURN END SUBROUTINE wake SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, & d_deltaqw, sigmaw, d_sigmaw, alpha) ! ------------------------------------------------------ ! Dtermination du coefficient alpha tel que les tendances ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent ! a une humidite positive dans la zone (x) et dans la zone (w). ! ------------------------------------------------------ IMPLICIT NONE ! Input REAL qe(nlon, nl), d_qe(nlon, nl) REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl) REAL sigmaw(nlon), d_sigmaw(nlon) LOGICAL wk_adv(nlon) INTEGER nl, nlon ! Output REAL alpha(nlon) ! Internal variables REAL zeta(nlon, nl) REAL alpha1(nlon) REAL x, a, b, c, discrim REAL epsilon ! DATA epsilon/1.e-15/ INTEGER i,k DO k = 1, nl DO i = 1, nlon IF (wk_adv(i)) THEN IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN zeta(i, k) = 0. ELSE zeta(i, k) = 1. END IF END IF END DO DO i = 1, nlon IF (wk_adv(i)) THEN x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + & (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i)*(deltaqw(i,k)+ & d_deltaqw(i,k)) a = -d_sigmaw(i)*d_deltaqw(i, k) b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - & deltaqw(i, k)*d_sigmaw(i) c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon discrim = b*b - 4.*a*c ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap alpha1(i) = 1. ELSE IF (x>=0.) THEN alpha1(i) = 1. ELSE IF (a>0.) THEN alpha1(i) = 0.9*min((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim & ))/(2.*a)) ELSE IF (a==0.) THEN alpha1(i) = 0.9*(-c/b) ELSE ! print*,'a,b,c discrim',a,b,c discrim alpha1(i) = 0.9*max((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim & ))/(2.*a)) END IF END IF END IF alpha(i) = min(alpha(i), alpha1(i)) END IF END DO END DO RETURN END SUBROUTINE wake_vec_modulation
1116  ! cc Définitions de entr, detr
1117  detr(i, k) = 0.
1118 
1119  entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
1120  sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)
1121 
1122  spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
1123  ! cc spread(i,k) =
1124  ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
1125  ! cc $ sigmaw(i)
1126 
1127 
1128  ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
1129  ! Jingmei
1130 
1131  ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
1132  ! & Tgw(i,k),deltatw(i,k)
1133  d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)*deltatw(i, k)* &
1134  dtimesub
1135  ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
1136  ff(i) = d_deltatw(i, k)/dtimesub
1137 
1138  ! Sans GW
1139 
1140  ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
1141 
1142  ! GW formule 1
1143 
1144  ! deltatw(k) = deltatw(k)+dtimesub*
1145  ! $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
1146 
1147  ! GW formule 2
1148 
1149  IF (dtimesub*tgw(i,k)<1.e-10) THEN
1150  d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc
1151  ! $
1152  ! -spread(i,k)*deltatw(i,k)
1153  -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( &
1154  i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
1155  -tgw(i,k)*deltatw(i,k))
1156  ELSE
1157  d_deltatw(i, k) = 1/tgw(i, k)*(1-exp(-dtimesub*tgw(i, &
1158  k)))*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc $
1159  ! -spread(i,k)*deltatw(i,k)
1160  -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( &
1161  i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
1162  -tgw(i,k)*deltatw(i,k))
1163  END IF
1164 
1165  dth(i, k) = deltatw(i, k)/ppi(i, k)
1166 
1167  gg(i) = d_deltaqw(i, k)/dtimesub
1168 
1169  d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k)+dqpbl(i,k) & ! cc $
1170  ! -spread(i,k)*deltaqw(i,k))
1171  -entr(i,k)*deltaqw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw(i)+detr( &
1172  i,k))*deltaqw(i,k)/(1.-sigmaw(i)))
1173  ! cc
1174 
1175  ! cc nrlmd
1176  ! cc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
1177  ! cc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
1178  ! cc
1179  END IF
1180  END DO
1181  END DO
1182 
1183 
1184  ! Scale tendencies so that water vapour remains positive in w and x.
1185 
1186  CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, &
1187  d_deltaqw, sigmaw, d_sigmaw, alpha)
1188 
1189  ! cc nrlmd
1190  ! c print*,'alpha'
1191  ! c do i=1,klon
1192  ! c print*,alpha(i)
1193  ! c end do
1194  ! cc
1195  DO k = 1, klev
1196  DO i = 1, klon
1197  IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1198  d_te(i, k) = alpha(i)*d_te(i, k)
1199  d_qe(i, k) = alpha(i)*d_qe(i, k)
1200  d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
1201  d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
1202  d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
1203  END IF
1204  END DO
1205  END DO
1206  DO i = 1, klon
1207  IF (wk_adv(i)) THEN
1208  d_sigmaw(i) = alpha(i)*d_sigmaw(i)
1209  END IF
1210  END DO
1211 
1212  ! Update large scale variables and wake variables
1213  ! IM 060208 manque DO i + remplace DO k=1,kupper(i)
1214  ! IM 060208 DO k = 1,kupper(i)
1215  DO k = 1, klev
1216  DO i = 1, klon
1217  IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1218  dtls(i, k) = dtls(i, k) + d_te(i, k)
1219  dqls(i, k) = dqls(i, k) + d_qe(i, k)
1220  ! cc nrlmd
1221  d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
1222  d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
1223  ! cc
1224  END IF
1225  END DO
1226  END DO
1227  DO k = 1, klev
1228  DO i = 1, klon
1229  IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1230  te(i, k) = te0(i, k) + dtls(i, k)
1231  qe(i, k) = qe0(i, k) + dqls(i, k)
1232  the(i, k) = te(i, k)/ppi(i, k)
1233  deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
1234  deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
1235  dth(i, k) = deltatw(i, k)/ppi(i, k)
1236  ! c print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k)
1237  ! c $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
1238  END IF
1239  END DO
1240  END DO
1241  DO i = 1, klon
1242  IF (wk_adv(i)) THEN
1243  sigmaw(i) = sigmaw(i) + d_sigmaw(i)
1244  END IF
1245  END DO
1246 
1247 
1248  ! Determine Ptop from buoyancy integral
1249  ! ---------------------------------------
1250 
1251  ! - 1/ Pressure of the level where dth changes sign.
1252 
1253  DO i = 1, klon
1254  IF (wk_adv(i)) THEN
1255  ptop_provis(i) = ph(i, 1)
1256  END IF
1257  END DO
1258 
1259  DO k = 2, klev
1260  DO i = 1, klon
1261  IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
1262  dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1263  ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
1264  k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1265  END IF
1266  END DO
1267  END DO
1268 
1269  ! - 2/ dth integral
1270 
1271  DO i = 1, klon
1272  IF (wk_adv(i)) THEN !!! nrlmd
1273  sum_dth(i) = 0.
1274  dthmin(i) = -delta_t_min
1275  z(i) = 0.
1276  END IF
1277  END DO
1278 
1279  DO k = 1, klev
1280  DO i = 1, klon
1281  IF (wk_adv(i)) THEN
1282  dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
1283  IF (dz(i)>0) THEN
1284  z(i) = z(i) + dz(i)
1285  sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1286  dthmin(i) = amin1(dthmin(i), dth(i,k))
1287  END IF
1288  END IF
1289  END DO
1290  END DO
1291 
1292  ! - 3/ height of triangle with area= sum_dth and base = dthmin
1293 
1294  DO i = 1, klon
1295  IF (wk_adv(i)) THEN
1296  hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
1297  hw(i) = amax1(hwmin, hw(i))
1298  END IF
1299  END DO
1300 
1301  ! - 4/ now, get Ptop
1302 
1303  DO i = 1, klon
1304  IF (wk_adv(i)) THEN !!! nrlmd
1305  ktop(i) = 0
1306  z(i) = 0.
1307  END IF
1308  END DO
1309 
1310  DO k = 1, klev
1311  DO i = 1, klon
1312  IF (wk_adv(i)) THEN
1313  dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i))
1314  IF (dz(i)>0) THEN
1315  z(i) = z(i) + dz(i)
1316  ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i)
1317  ktop(i) = k
1318  END IF
1319  END IF
1320  END DO
1321  END DO
1322 
1323  ! 4.5/Correct ktop and ptop
1324 
1325  DO i = 1, klon
1326  IF (wk_adv(i)) THEN
1327  ptop_new(i) = ptop(i)
1328  END IF
1329  END DO
1330 
1331  DO k = klev, 2, -1
1332  DO i = 1, klon
1333  ! IM v3JYG; IF (k .GE. ktop(i)
1334  IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
1335  dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1336  ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
1337  k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1338  END IF
1339  END DO
1340  END DO
1341 
1342 
1343  DO i = 1, klon
1344  IF (wk_adv(i)) THEN
1345  ptop(i) = ptop_new(i)
1346  END IF
1347  END DO
1348 
1349  DO k = klev, 1, -1
1350  DO i = 1, klon
1351  IF (wk_adv(i)) THEN !!! nrlmd
1352  IF (ph(i,k+1)<ptop(i)) ktop(i) = k
1353  END IF
1354  END DO
1355  END DO
1356 
1357  ! 5/ Set deltatw & deltaqw to 0 above kupper
1358 
1359  DO k = 1, klev
1360  DO i = 1, klon
1361  IF (wk_adv(i) .AND. k>=kupper(i)) THEN
1362  deltatw(i, k) = 0.
1363  deltaqw(i, k) = 0.
1364  END IF
1365  END DO
1366  END DO
1367 
1368 
1369  ! -------------Cstar computation---------------------------------
1370  DO i = 1, klon
1371  IF (wk_adv(i)) THEN !!! nrlmd
1372  sum_thu(i) = 0.
1373  sum_tu(i) = 0.
1374  sum_qu(i) = 0.
1375  sum_thvu(i) = 0.
1376  sum_dth(i) = 0.
1377  sum_dq(i) = 0.
1378  sum_rho(i) = 0.
1379  sum_dtdwn(i) = 0.
1380  sum_dqdwn(i) = 0.
1381 
1382  av_thu(i) = 0.
1383  av_tu(i) = 0.
1384  av_qu(i) = 0.
1385  av_thvu(i) = 0.
1386  av_dth(i) = 0.
1387  av_dq(i) = 0.
1388  av_rho(i) = 0.
1389  av_dtdwn(i) = 0.
1390  av_dqdwn(i) = 0.
1391  END IF
1392  END DO
1393 
1394  ! Integrals (and wake top level number)
1395  ! --------------------------------------
1396 
1397  ! Initialize sum_thvu to 1st level virt. pot. temp.
1398 
1399  DO i = 1, klon
1400  IF (wk_adv(i)) THEN !!! nrlmd
1401  z(i) = 1.
1402  dz(i) = 1.
1403  sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i)
1404  sum_dth(i) = 0.
1405  END IF
1406  END DO
1407 
1408  DO k = 1, klev
1409  DO i = 1, klon
1410  IF (wk_adv(i)) THEN !!! nrlmd
1411  dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
1412  IF (dz(i)>0) THEN
1413  z(i) = z(i) + dz(i)
1414  sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
1415  sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
1416  sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
1417  sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i)
1418  sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1419  sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1420  sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
1421  sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
1422  sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
1423  END IF
1424  END IF
1425  END DO
1426  END DO
1427 
1428  DO i = 1, klon
1429  IF (wk_adv(i)) THEN !!! nrlmd
1430  hw0(i) = z(i)
1431  END IF
1432  END DO
1433 
1434 
1435  ! - WAPE and mean forcing computation
1436  ! ---------------------------------------
1437 
1438  ! ---------------------------------------
1439 
1440  ! Means
1441 
1442  DO i = 1, klon
1443  IF (wk_adv(i)) THEN !!! nrlmd
1444  av_thu(i) = sum_thu(i)/hw0(i)
1445  av_tu(i) = sum_tu(i)/hw0(i)
1446  av_qu(i) = sum_qu(i)/hw0(i)
1447  av_thvu(i) = sum_thvu(i)/hw0(i)
1448  av_dth(i) = sum_dth(i)/hw0(i)
1449  av_dq(i) = sum_dq(i)/hw0(i)
1450  av_rho(i) = sum_rho(i)/hw0(i)
1451  av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1452  av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1453 
1454  wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* &
1455  av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
1456  END IF
1457  END DO
1458 
1459  ! Filter out bad wakes
1460 
1461  DO k = 1, klev
1462  DO i = 1, klon
1463  IF (wk_adv(i)) THEN !!! nrlmd
1464  IF (wape(i)<0.) THEN
1465  deltatw(i, k) = 0.
1466  deltaqw(i, k) = 0.
1467  dth(i, k) = 0.
1468  END IF
1469  END IF
1470  END DO
1471  END DO
1472 
1473  DO i = 1, klon
1474  IF (wk_adv(i)) THEN !!! nrlmd
1475  IF (wape(i)<0.) THEN
1476  wape(i) = 0.
1477  cstar(i) = 0.
1478  hw(i) = hwmin
1479  sigmaw(i) = max(sigmad, sigd_con(i))
1480  fip(i) = 0.
1481  gwake(i) = .false.
1482  ELSE
1483  cstar(i) = stark*sqrt(2.*wape(i))
1484  gwake(i) = .true.
1485  END IF
1486  END IF
1487  END DO
1488 
1489  END DO ! end sub-timestep loop
1490 
1491  ! -----------------------------------------------------------------
1492  ! Get back to tendencies per second
1493 
1494  DO k = 1, klev
1495  DO i = 1, klon
1496 
1497  ! cc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
1498  IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
1499  ! cc
1500  dtls(i, k) = dtls(i, k)/dtime
1501  dqls(i, k) = dqls(i, k)/dtime
1502  d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
1503  d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
1504  d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
1505  ! c print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
1506  ! c $ ,death_rate(i)*sigmaw(i)
1507  END IF
1508  END DO
1509  END DO
1510 
1511 
1512  ! ----------------------------------------------------------
1513  ! Determine wake final state; recompute wape, cstar, ktop;
1514  ! filter out bad wakes.
1515  ! ----------------------------------------------------------
1516 
1517  ! 2.1 - Undisturbed area and Wake integrals
1518  ! ---------------------------------------------------------
1519 
1520  DO i = 1, klon
1521  ! cc nrlmd if (wk_adv(i)) then !!! nrlmd
1522  IF (ok_qx_qw(i)) THEN
1523  ! cc
1524  z(i) = 0.
1525  sum_thu(i) = 0.
1526  sum_tu(i) = 0.
1527  sum_qu(i) = 0.
1528  sum_thvu(i) = 0.
1529  sum_dth(i) = 0.
1530  sum_dq(i) = 0.
1531  sum_rho(i) = 0.
1532  sum_dtdwn(i) = 0.
1533  sum_dqdwn(i) = 0.
1534 
1535  av_thu(i) = 0.
1536  av_tu(i) = 0.
1537  av_qu(i) = 0.
1538  av_thvu(i) = 0.
1539  av_dth(i) = 0.
1540  av_dq(i) = 0.
1541  av_rho(i) = 0.
1542  av_dtdwn(i) = 0.
1543  av_dqdwn(i) = 0.
1544  END IF
1545  END DO
1546  ! Potential temperatures and humidity
1547  ! ----------------------------------------------------------
1548 
1549  DO k = 1, klev
1550  DO i = 1, klon
1551  ! cc nrlmd IF ( wk_adv(i)) THEN
1552  IF (ok_qx_qw(i)) THEN
1553  ! cc
1554  rho(i, k) = p(i, k)/(rd*te(i,k))
1555  IF (k==1) THEN
1556  rhoh(i, k) = ph(i, k)/(rd*te(i,k))
1557  zhh(i, k) = 0
1558  ELSE
1559  rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
1560  zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
1561  END IF
1562  the(i, k) = te(i, k)/ppi(i, k)
1563  thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
1564  tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
1565  qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
1566  rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
1567  dth(i, k) = deltatw(i, k)/ppi(i, k)
1568  END IF
1569  END DO
1570  END DO
1571 
1572  ! Integrals (and wake top level number)
1573  ! -----------------------------------------------------------
1574 
1575  ! Initialize sum_thvu to 1st level virt. pot. temp.
1576 
1577  DO i = 1, klon
1578  ! cc nrlmd IF ( wk_adv(i)) THEN
1579  IF (ok_qx_qw(i)) THEN
1580  ! cc
1581  z(i) = 1.
1582  dz(i) = 1.
1583  sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i)
1584  sum_dth(i) = 0.
1585  END IF
1586  END DO
1587 
1588  DO k = 1, klev
1589  DO i = 1, klon
1590  ! cc nrlmd IF ( wk_adv(i)) THEN
1591  IF (ok_qx_qw(i)) THEN
1592  ! cc
1593  dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
1594  IF (dz(i)>0) THEN
1595  z(i) = z(i) + dz(i)
1596  sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
1597  sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
1598  sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
1599  sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i)
1600  sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1601  sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1602  sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
1603  sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
1604  sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
1605  END IF
1606  END IF
1607  END DO
1608  END DO
1609 
1610  DO i = 1, klon
1611  ! cc nrlmd IF ( wk_adv(i)) THEN
1612  IF (ok_qx_qw(i)) THEN
1613  ! cc
1614  hw0(i) = z(i)
1615  END IF
1616  END DO
1617 
1618  ! - WAPE and mean forcing computation
1619  ! -------------------------------------------------------------
1620 
1621  ! Means
1622 
1623  DO i = 1, klon
1624  ! cc nrlmd IF ( wk_adv(i)) THEN
1625  IF (ok_qx_qw(i)) THEN
1626  ! cc
1627  av_thu(i) = sum_thu(i)/hw0(i)
1628  av_tu(i) = sum_tu(i)/hw0(i)
1629  av_qu(i) = sum_qu(i)/hw0(i)
1630  av_thvu(i) = sum_thvu(i)/hw0(i)
1631  av_dth(i) = sum_dth(i)/hw0(i)
1632  av_dq(i) = sum_dq(i)/hw0(i)
1633  av_rho(i) = sum_rho(i)/hw0(i)
1634  av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1635  av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1636 
1637  wape2(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* &
1638  av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
1639  END IF
1640  END DO
1641 
1642  ! Prognostic variable update
1643  ! ------------------------------------------------------------
1644 
1645  ! Filter out bad wakes
1646 
1647  DO k = 1, klev
1648  DO i = 1, klon
1649  ! cc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
1650  IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
1651  ! cc
1652  deltatw(i, k) = 0.
1653  deltaqw(i, k) = 0.
1654  dth(i, k) = 0.
1655  END IF
1656  END DO
1657  END DO
1658 
1659 
1660  DO i = 1, klon
1661  ! cc nrlmd IF ( wk_adv(i)) THEN
1662  IF (ok_qx_qw(i)) THEN
1663  ! cc
1664  IF (wape2(i)<0.) THEN
1665  wape2(i) = 0.
1666  cstar2(i) = 0.
1667  hw(i) = hwmin
1668  sigmaw(i) = amax1(sigmad, sigd_con(i))
1669  fip(i) = 0.
1670  gwake(i) = .false.
1671  ELSE
1672  IF (prt_level>=10) print *, 'wape2>0'
1673  cstar2(i) = stark*sqrt(2.*wape2(i))
1674  gwake(i) = .true.
1675  END IF
1676  END IF
1677  END DO
1678 
1679  DO i = 1, klon
1680  ! cc nrlmd IF ( wk_adv(i)) THEN
1681  IF (ok_qx_qw(i)) THEN
1682  ! cc
1683  ktopw(i) = ktop(i)
1684  END IF
1685  END DO
1686 
1687  DO i = 1, klon
1688  ! cc nrlmd IF ( wk_adv(i)) THEN
1689  IF (ok_qx_qw(i)) THEN
1690  ! cc
1691  IF (ktopw(i)>0 .AND. gwake(i)) THEN
1692 
1693  ! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer)
1694  ! cc heff = 600.
1695  ! Utilisation de la hauteur hw
1696  ! c heff = 0.7*hw
1697  heff(i) = hw(i)
1698 
1699  fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* &
1700  sqrt(sigmaw(i)*wdens(i)*3.14)
1701  fip(i) = alpk*fip(i)
1702  ! jyg2
1703  ELSE
1704  fip(i) = 0.
1705  END IF
1706  END IF
1707  END DO
1708 
1709  ! Limitation de sigmaw
1710 
1711  ! cc nrlmd
1712  ! DO i=1,klon
1713  ! IF (OK_qx_qw(i)) THEN
1714  ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
1715  ! ENDIF
1716  ! ENDDO
1717  ! cc
1718  DO k = 1, klev
1719  DO i = 1, klon
1720 
1721  ! cc nrlmd On maintient désormais constant sigmaw en régime
1722  ! permanent
1723  ! cc IF ((sigmaw(i).GT.sigmaw_max).or.
1724  IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
1725  .NOT. ok_qx_qw(i)) THEN
1726  ! cc
1727  dtls(i, k) = 0.
1728  dqls(i, k) = 0.
1729  deltatw(i, k) = 0.
1730  deltaqw(i, k) = 0.
1731  END IF
1732  END DO
1733  END DO
1734 
1735  ! cc nrlmd On maintient désormais constant sigmaw en régime permanent
1736  DO i = 1, klon
1737  IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
1738  .NOT. ok_qx_qw(i)) THEN
1739  wape(i) = 0.
1740  cstar(i) = 0.
1741 !!jyg Outside subroutine "Wake" hw and sigmaw are zero when there are no wakes
1742 !! hw(i) = hwmin !jyg
1743 !! sigmaw(i) = sigmad !jyg
1744  hw(i) = 0. !jyg
1745  sigmaw(i) = 0. !jyg
1746  fip(i) = 0.
1747  ELSE
1748  wape(i) = wape2(i)
1749  cstar(i) = cstar2(i)
1750  END IF
1751  ! c print*,'wape wape2 ktopw OK_qx_qw =',
1752  ! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
1753  END DO
1754 
1755 
1756  RETURN
1757 END SUBROUTINE wake
1758 
1759 SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, &
1760  d_deltaqw, sigmaw, d_sigmaw, alpha)
1761  ! ------------------------------------------------------
1762  ! Dtermination du coefficient alpha tel que les tendances
1763  ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
1764  ! a une humidite positive dans la zone (x) et dans la zone (w).
1765  ! ------------------------------------------------------
1766  IMPLICIT NONE
1767 
1768  ! Input
1769  REAL qe(nlon, nl), d_qe(nlon, nl)
1770  REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
1771  REAL sigmaw(nlon), d_sigmaw(nlon)
1772  LOGICAL wk_adv(nlon)
1773  INTEGER nl, nlon
1774  ! Output
1775  REAL alpha(nlon)
1776  ! Internal variables
1777  REAL zeta(nlon, nl)
1778  REAL alpha1(nlon)
1779  REAL x, a, b, c, discrim
1780  REAL epsilon
1781  ! DATA epsilon/1.e-15/
1782  INTEGER i,k
1783 
1784  DO k = 1, nl
1785  DO i = 1, nlon
1786  IF (wk_adv(i)) THEN
1787  IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
1788  zeta(i, k) = 0.
1789  ELSE
1790  zeta(i, k) = 1.
1791  END IF
1792  END IF
1793  END DO
1794  DO i = 1, nlon
1795  IF (wk_adv(i)) THEN
1796  x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
1797  (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i)*(deltaqw(i,k)+ &
1798  d_deltaqw(i,k))
1799  a = -d_sigmaw(i)*d_deltaqw(i, k)
1800  b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
1801  deltaqw(i, k)*d_sigmaw(i)
1802  c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon
1803  discrim = b*b - 4.*a*c
1804  ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
1805  IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap
1806  alpha1(i) = 1.
1807  ELSE
1808  IF (x>=0.) THEN
1809  alpha1(i) = 1.
1810  ELSE
1811  IF (a>0.) THEN
1812  alpha1(i) = 0.9*min((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
1813  ))/(2.*a))
1814  ELSE IF (a==0.) THEN
1815  alpha1(i) = 0.9*(-c/b)
1816  ELSE
1817  ! print*,'a,b,c discrim',a,b,c discrim
1818  alpha1(i) = 0.9*max((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
1819  ))/(2.*a))
1820  END IF
1821  END IF
1822  END IF
1823  alpha(i) = min(alpha(i), alpha1(i))
1824  END IF
1825  END DO
1826  END DO
1827 
1828  RETURN
1829 END SUBROUTINE wake_vec_modulation
1830 
1831 
bcast
!IM Implemente en modes sequentiel et parallele CALL rlon_glo CALL bcast(rlon_glo)!$OMP MASTER if(is_mpi_root) then!zstophy
dimphy::klon
integer, save klon
Definition: dimphy.F90:3
dimphy::klev
integer, save klev
Definition: dimphy.F90:7
eps
!$Id!Thermodynamical constants for t0 real clmci real eps
Definition: cvthermo.h:6
hw
!$Id hw
Definition: 1Dconv.h:4
err
subroutine err(ierr, typ, nam)
Definition: dynetat0.f90:189
false
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
Definition: calcul_STDlev.h:26
wake_vec_modulation
subroutine wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, d_deltaqw, sigmaw, d_sigmaw, alpha)
Definition: wake.F90:1761
nl
!$Id!Parameters for nl
Definition: cv30param.h:5
mod_phys_lmdz_para
Definition: mod_phys_lmdz_para.F90:4
prt_level
!FH On elimine toutes les clefs physiques dans la dynamique prt_level
Definition: write_paramLMDZ_dyn.h:9
true
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
Definition: calcul_STDlev.h:139
alpha
!$Id!Parameters for nlm real spfac!IM cf epmax real ptcrit real omtrain real dttrig real alpha real delta real betad COMMON cv30param nlm spfac &!IM cf ptcrit omtrain dttrig alpha
Definition: cv30param.h:5
print_control_mod
Definition: print_control_mod.F90:2
wake
subroutine wake(p, ph, pi, dtime, sigd_con, te0, qe0, omgb, dtdwn, dqdwn, amdwn, amup, dta, dqa, wdtpbl, wdqpbl, udtpbl, udqpbl, deltatw, deltaqw, dth, hw, sigmaw, wape, fip, gfl, dtls, dqls, ktopw, omgbdth, dp_omgb, wdens, tu, qu, dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, spread, cstar, d_deltat_gw, d_deltatw2, d_deltaqw2)
Definition: wake.F90:12
dimphy
Definition: dimphy.F90:1
rg
real rg
Definition: comcstphy.h:1