Index: LMDZ6/trunk/libf/phylmd/calwake.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/calwake.F90 (revision 4586) +++ LMDZ6/trunk/libf/phylmd/calwake.F90 (revision 4588) @@ -29,4 +29,5 @@ USE indice_sol_mod, ONLY: is_oce USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level + USE lmdz_wake, ONLY : wake IMPLICIT NONE ! ====================================================================== Index: LMDZ6/trunk/libf/phylmd/lmdz_wake.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_wake.F90 (revision 4588) +++ LMDZ6/trunk/libf/phylmd/lmdz_wake.F90 (revision 4588) @@ -0,0 +1,2670 @@ +MODULE lmdz_wake + +! $Id$ + +CONTAINS + +SUBROUTINE wake(klon,klev,znatsurf, p, ph, pi, dtime, & + tenv0, qe0, omgb, & + dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, & + sigd_con, Cin, & + deltatw, deltaqw, sigmaw, awdens, wdens, & ! state variables + dth, hw, wape, fip, gfl, & + dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, & + dtke, dqke, omg, dp_deltomg, wkspread, cstar, & + d_deltat_gw, & ! tendencies + d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2) ! tendencies + + + ! ************************************************************** + ! * + ! WAKE * + ! retour a un Pupper fixe * + ! * + ! written by : GRANDPEIX Jean-Yves 09/03/2000 * + ! modified by : ROEHRIG Romain 01/29/2007 * + ! ************************************************************** + + + USE lmdz_wake_ini , ONLY : wake_ini + USE lmdz_wake_ini , ONLY : prt_level,epsim1,RG,RD + USE lmdz_wake_ini , ONLY : stark, wdens_ref, coefgw, alpk, wk_pupper + USE lmdz_wake_ini , ONLY : crep_upper, crep_sol, tau_cv, rzero, aa0, flag_wk_check_trgl + USE lmdz_wake_ini , ONLY : iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensmin + USE lmdz_wake_ini , ONLY : sigmad, hwmin, wapecut, cstart, sigmaw_max, dens_rate, epsilon_loc + USE lmdz_wake_ini , ONLY : iflag_wk_profile + + + IMPLICIT NONE + ! ============================================================================ + + + ! But : Decrire le comportement des poches froides apparaissant dans les + ! grands systemes convectifs, et fournir l'energie disponible pour + ! le declenchement de nouvelles colonnes convectives. + + ! State variables : + ! deltatw : temperature difference between wake and off-wake regions + ! deltaqw : specific humidity difference between wake and off-wake regions + ! sigmaw : fractional area covered by wakes. + ! wdens : number of wakes per unit area + + ! Variable de sortie : + + ! wape : WAke Potential Energy + ! fip : Front Incident Power (W/m2) - ALP + ! gfl : Gust Front Length per unit area (m-1) + ! dtls : large scale temperature tendency due to wake + ! dqls : large scale humidity tendency due to wake + ! hw : wake top hight (given by hw*deltatw(1)/2=wape) + ! dp_omgb : vertical gradient of large scale omega + ! awdens : densite de poches actives + ! wdens : densite de poches + ! omgbdth: flux of Delta_Theta transported by LS omega + ! dtKE : differential heating (wake - unpertubed) + ! dqKE : differential moistening (wake - unpertubed) + ! omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s) + ! dp_deltomg : vertical gradient of omg (s-1) + ! wkspread : spreading term in d_t_wake and d_q_wake + ! deltatw : updated temperature difference (T_w-T_u). + ! deltaqw : updated humidity difference (q_w-q_u). + ! sigmaw : updated wake fractional area. + ! d_deltat_gw : delta T tendency due to GW + + ! Variables d'entree : + + ! aire : aire de la maille + ! tenv0 : temperature dans l'environnement (K) + ! qe0 : humidite dans l'environnement (kg/kg) + ! omgb : vitesse verticale moyenne sur la maille (Pa/s) + ! dtdwn: source de chaleur due aux descentes (K/s) + ! dqdwn: source d'humidite due aux descentes (kg/kg/s) + ! dta : source de chaleur due courants satures et detrain (K/s) + ! dqa : source d'humidite due aux courants satures et detra (kg/kg/s) + ! wgen : number of wakes generated per unit area and per sec (/m^2/s) + ! amdwn: flux de masse total des descentes, par unite de + ! surface de la maille (kg/m2/s) + ! amup : flux de masse total des ascendances, par unite de + ! surface de la maille (kg/m2/s) + ! sigd_con: + ! Cin : convective inhibition + ! p : pressions aux milieux des couches (Pa) + ! ph : pressions aux interfaces (Pa) + ! pi : (p/p_0)**kapa (adim) + ! dtime: increment temporel (s) + + ! Variables internes : + + ! rhow : masse volumique de la poche froide + ! rho : environment density at P levels + ! rhoh : environment density at Ph levels + ! tenv : environment temperature | may change within + ! qe : environment humidity | sub-time-stepping + ! the : environment potential temperature + ! thu : potential temperature in undisturbed area + ! tu : temperature in undisturbed area + ! qu : humidity in undisturbed area + ! dp_omgb: vertical gradient og LS omega + ! omgbw : wake average vertical omega + ! dp_omgbw: vertical gradient of omgbw + ! omgbdq : flux of Delta_q transported by LS omega + ! dth : potential temperature diff. wake-undist. + ! th1 : first pot. temp. for vertical advection (=thu) + ! th2 : second pot. temp. for vertical advection (=thw) + ! q1 : first humidity for vertical advection + ! q2 : second humidity for vertical advection + ! d_deltatw : terme de redistribution pour deltatw + ! d_deltaqw : terme de redistribution pour deltaqw + ! deltatw0 : deltatw initial + ! deltaqw0 : deltaqw initial + ! hw0 : wake top hight (defined as the altitude at which deltatw=0) + ! amflux : horizontal mass flux through wake boundary + ! wdens_ref: initial number of wakes per unit area (3D) or per + ! unit length (2D), at the beginning of each time step + ! Tgw : 1 sur la periode de onde de gravite + ! Cgw : vitesse de propagation de onde de gravite + ! LL : distance entre 2 poches + + ! ------------------------------------------------------------------------- + ! Declaration de variables + ! ------------------------------------------------------------------------- + + + ! Arguments en entree + ! -------------------- + + INTEGER, INTENT(IN) :: klon,klev + INTEGER, DIMENSION (klon), INTENT(IN) :: znatsurf + REAL, DIMENSION (klon, klev), INTENT(IN) :: p, pi + REAL, DIMENSION (klon, klev+1), INTENT(IN) :: ph + REAL, DIMENSION (klon, klev), INTENT(IN) :: omgb + REAL, INTENT(IN) :: dtime + REAL, DIMENSION (klon, klev), INTENT(IN) :: tenv0, qe0 + REAL, DIMENSION (klon, klev), INTENT(IN) :: dtdwn, dqdwn + REAL, DIMENSION (klon, klev), INTENT(IN) :: amdwn, amup + REAL, DIMENSION (klon, klev), INTENT(IN) :: dta, dqa + REAL, DIMENSION (klon), INTENT(IN) :: wgen + REAL, DIMENSION (klon), INTENT(IN) :: sigd_con + REAL, DIMENSION (klon), INTENT(IN) :: Cin + + ! + ! Input/Output + ! State variables + REAL, DIMENSION (klon, klev), INTENT(INOUT) :: deltatw, deltaqw + REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw + REAL, DIMENSION (klon), INTENT(INOUT) :: awdens + REAL, DIMENSION (klon), INTENT(INOUT) :: wdens + + ! 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) :: wkspread ! unused (jyg) + REAL, DIMENSION (klon, klev), INTENT(OUT) :: omgbdth, omg + REAL, DIMENSION (klon, klev), INTENT(OUT) :: dp_omgb, dp_deltomg + REAL, DIMENSION (klon), INTENT(OUT) :: hw, wape, fip, gfl, cstar + INTEGER, DIMENSION (klon), INTENT(OUT) :: ktopw + ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields + ! computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens) + REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_gw + REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltatw2, d_deltaqw2 + REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw2, d_awdens2, d_wdens2 + + ! Variables internes + ! ------------------- + + ! Variables a fixer + + REAL :: delta_t_min + INTEGER :: nsub + REAL :: dtimesub + 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) :: tenv, qe +!! REAL, DIMENSION (klon) :: sigmaw1 + + ! Variables liees a la dynamique de population 1 + REAL, DIMENSION(klon) :: act + REAL, DIMENSION(klon) :: rad_wk, tau_wk_inv + REAL, DIMENSION(klon) :: f_shear + REAL, DIMENSION(klon) :: drdt + + ! Variables liees a la dynamique de population 2 + REAL, DIMENSION(klon) :: cont_fact + +!! REAL, DIMENSION(klon) :: d_sig_gen, d_sig_death, d_sig_col + REAL, DIMENSION(klon) :: wape1_act, wape2_act + LOGICAL, DIMENSION (klon) :: kill_wake + REAL :: drdt_pos + REAL :: tau_wk_inv_min + ! Some components of the tendencies of state variables + REAL, DIMENSION (klon) :: d_sig_gen2, d_sig_death2, d_sig_col2, d_sig_spread2, d_sig_bnd2 + REAL, DIMENSION (klon) :: d_dens_gen2, d_dens_death2, d_dens_col2, d_dens_bnd2 + REAL, DIMENSION (klon) :: d_adens_death2, d_adens_icol2, d_adens_acol2, d_adens_bnd2 + + ! Variables pour les GW + REAL, DIMENSION (klon) :: ll + REAL, DIMENSION (klon, klev) :: n2 + REAL, DIMENSION (klon, klev) :: cgw + REAL, DIMENSION (klon, klev) :: tgw + + ! Variables liees 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 + + ! Variables liees au test de la forme triangulaire du profil de Delta_theta + REAL, DIMENSION (klon) :: sum_half_dth + REAL, DIMENSION (klon) :: dz_half + + ! Sub-timestep tendencies and related variables + REAL, DIMENSION (klon, klev) :: d_deltatw, d_deltaqw + REAL, DIMENSION (klon, klev) :: d_tenv, d_qe + REAL, DIMENSION (klon) :: d_awdens, d_wdens, d_sigmaw + REAL, DIMENSION (klon) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd + REAL, DIMENSION (klon) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd + REAL, DIMENSION (klon) :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd + REAL, DIMENSION (klon) :: alpha, alpha_tot + REAL, DIMENSION (klon) :: q0_min, q1_min + LOGICAL, DIMENSION (klon) :: wk_adv, ok_qx_qw + + ! Autres variables internes + INTEGER ::isubstep, k, i, igout + + REAL :: sigmaw_targ + REAL :: wdens_targ + REAL :: d_sigmaw_targ + REAL :: d_wdens_targ + + 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) :: 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 +!! REAL, DIMENSION (klon) :: nat_rate + REAL, DIMENSION (klon, klev) :: entr + REAL, DIMENSION (klon, klev) :: detr + + REAL, DIMENSION(klon) :: sigmaw_in ! pour les prints + REAL, DIMENSION(klon) :: awdens_in, wdens_in ! pour les prints + + ! ------------------------------------------------------------------------- + ! Initialisations + ! ------------------------------------------------------------------------- + ! ALON = 3.e5 + ! alon = 1.E6 + + ! Provisionnal; to be suppressed when f_shear is parameterized + f_shear(:) = 1. ! 0. for strong shear, 1. for weak shear + + + ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei) + + ! coefgw : Coefficient pour les ondes de gravite + ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) + ! wdens : Densite surfacique de poche froide + ! ------------------------------------------------------------------------- + + ! 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 +!print *,'XXXX dtime input ', dtime + igout = klon/2+1/klon + + IF (iflag_wk_pop_dyn == 0) THEN + ! 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) + !jyg< + !! wdens(:) = wdens_ref + DO i = 1,klon + wdens(i) = wdens_ref(znatsurf(i)+1) + ENDDO + !>jyg + ENDIF ! (iflag_wk_pop_dyn == 0) + + ! 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, initialize tendencies, initialize output fields + ! ------------------------------------------------------------------------ + +!jyg< +!! 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) +!! tenv(i, k) = tenv0(i, k) +!! qe(i, k) = qe0(i, k) +!! dtls(i, k) = 0. +!! dqls(i, k) = 0. +!! d_deltat_gw(i, k) = 0. +!! d_tenv(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 + ppi(:,:) = pi(:,:) + deltatw0(:,:) = deltatw(:,:) + deltaqw0(:,:) = deltaqw(:,:) + tenv(:,:) = tenv0(:,:) + qe(:,:) = qe0(:,:) + dtls(:,:) = 0. + dqls(:,:) = 0. + d_deltat_gw(:,:) = 0. + d_tenv(:,:) = 0. + d_qe(:,:) = 0. + d_deltatw(:,:) = 0. + d_deltaqw(:,:) = 0. + d_deltatw2(:,:) = 0. + d_deltaqw2(:,:) = 0. + + IF (iflag_wk_act == 0) THEN + act(:) = 0. + ELSEIF (iflag_wk_act == 1) THEN + act(:) = 1. + ENDIF + +!! DO i = 1, klon +!! sigmaw_in(i) = sigmaw(i) +!! END DO + sigmaw_in(:) = sigmaw(:) +!>jyg +! + IF (iflag_wk_pop_dyn >= 1) THEN + awdens_in(:) = awdens(:) + wdens_in(:) = wdens(:) +!! wdens(:) = wdens(:) + wgen(:)*dtime +!! d_wdens2(:) = wgen(:)*dtime +!! ELSE + ENDIF ! (iflag_wk_pop_dyn >= 1) + + + ! sigmaw1=sigmaw + ! IF (sigd_con.GT.sigmaw1) THEN + ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con + ! ENDIF + IF (iflag_wk_pop_dyn >= 1) THEN + DO i = 1, klon + d_dens_gen2(i) = 0. + d_dens_death2(i) = 0. + d_dens_col2(i) = 0. + d_awdens2(i) = 0. +! + wdens_targ = max(wdens(i),wdensmin) + d_dens_bnd2(i) = wdens_targ - wdens(i) + d_wdens2(i) = wdens_targ - wdens(i) + wdens(i) = wdens_targ + END DO + IF (iflag_wk_pop_dyn == 2) THEN + DO i = 1, klon + d_adens_death2(i) = 0. + d_adens_icol2(i) = 0. + d_adens_acol2(i) = 0. +! + wdens_targ = min(max(awdens(i),0.),wdens(i)) + d_adens_bnd2(i) = wdens_targ - awdens(i) + d_awdens2(i) = wdens_targ - awdens(i) + awdens(i) = wdens_targ + END DO + ENDIF ! (iflag_wk_pop_dyn == 2) + ELSE + DO i = 1, klon + d_awdens2(i) = 0. + d_wdens2(i) = 0. + END DO + ENDIF ! (iflag_wk_pop_dyn >= 1) +! + DO i = 1, klon + ! c sigmaw(i) = amax1(sigmaw(i),sigd_con(i)) +!jyg< +!! sigmaw(i) = amax1(sigmaw(i), sigmad) +!! sigmaw(i) = amin1(sigmaw(i), 0.99) + d_sig_gen2(i) = 0. + d_sig_death2(i) = 0. + d_sig_col2(i) = 0. + d_sig_spread2(i)= 0. + sigmaw_targ = min(max(sigmaw(i), sigmad),0.99) + d_sig_bnd2(i) = sigmaw_targ - sigmaw(i) + d_sigmaw2(i) = sigmaw_targ - sigmaw(i) +! print *,'XXXX1 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) + sigmaw(i) = sigmaw_targ +!>jyg + END DO + + wape(:) = 0. + wape2(:) = 0. + d_sigmaw(:) = 0. + ktopw(:) = 0 +! +!jyg +! + IF (prt_level>=10) THEN + PRINT *, 'wake-1, sigmaw(igout) ', sigmaw(igout) + PRINT *, 'wake-1, deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev) + PRINT *, 'wake-1, deltaqw(igout,k) ', (k,deltaqw(igout,k), k=1,klev) + PRINT *, 'wake-1, dowwdraughts, amdwn(igout,k) ', (k,amdwn(igout,k), k=1,klev) + PRINT *, 'wake-1, dowwdraughts, dtdwn(igout,k) ', (k,dtdwn(igout,k), k=1,klev) + PRINT *, 'wake-1, dowwdraughts, dqdwn(igout,k) ', (k,dqdwn(igout,k), k=1,klev) + PRINT *, 'wake-1, updraughts, amup(igout,k) ', (k,amup(igout,k), k=1,klev) + PRINT *, 'wake-1, updraughts, dta(igout,k) ', (k,dta(igout,k), k=1,klev) + PRINT *, 'wake-1, updraughts, dqa(igout,k) ', (k,dqa(igout,k), k=1,klev) + ENDIF + + ! 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,tenv(i,k) + rho(i, k) = p(i, k)/(RD*tenv(i,k)) + ! write(*,*)'wake 2',rho(i,k) + IF (k==1) THEN + ! write(*,*)'wake 3',i,k,rd,tenv(i,k) + rhoh(i, k) = ph(i, k)/(RD*tenv(i,k)) + ! write(*,*)'wake 4',i,k,rd,tenv(i,k) + zhh(i, k) = 0 + ELSE + ! write(*,*)'wake 5',rd,(tenv(i,k)+tenv(i,k-1)) + rhoh(i, k) = ph(i, k)*2./(RD*(tenv(i,k)+tenv(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) = tenv(i, k)/ppi(i, k) + thu(i, k) = (tenv(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) + tu(i, k) = tenv(i, k) - deltatw(i, k)*sigmaw(i) + qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) + ! write(*,*)'wake 8',(RD*(tenv(i,k)+deltatw(i,k))) + rhow(i, k) = p(i, k)/(RD*(tenv(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 + ! ----------------------------------------------------------------- + + ! 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 + + IF (prt_level>=10) THEN + PRINT *, 'wake-2, ptop_provis(igout), ptop(igout) ', ptop_provis(igout), ptop(igout) + ENDIF + + + ! -5/ Determination de ktop et kupper + + CALL pkupper (klon, klev, ptop, ph, pupper, kupper) + + DO k = klev, 1, -1 + DO i = 1, klon + IF (ph(i,k+1)-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)=10) THEN + PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout) + ENDIF + + ! -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. + d_deltatw2(i,k) = -deltatw0(i,k) + d_deltaqw2(i,k) = -deltaqw0(i,k) + 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.+epsim1*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.+epsim1*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)+ & + epsim1*(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. + d_deltatw2(i,k) = -deltatw0(i,k) + d_deltaqw2(i,k) = -deltaqw0(i,k) + END IF + END DO + END DO + + DO i = 1, klon + IF (wape(i)<0.) THEN + wape(i) = 0. + cstar(i) = 0. + hw(i) = hwmin +!jyg< +!! sigmaw(i) = amax1(sigmad, sigd_con(i)) + sigmaw_targ = max(sigmad, sigd_con(i)) + d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) + d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) +! print *,'XXXX2 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) + sigmaw(i) = sigmaw_targ +!>jyg + fip(i) = 0. + gwake(i) = .FALSE. + ELSE + hw(i) = hw0(i) + 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. + alpha_tot(i) = 1. + END DO + + IF (prt_level>=10) THEN + PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ', & + sigmaw(igout), cstar(igout), wape(igout), ktop(igout) + ENDIF + + + ! C ----------------------------------------------------------------- + ! Sub-time-stepping + ! ----------------- + + nsub = 10 + dtimesub = dtime/nsub + + + + ! ------------------------------------------------------------ + DO isubstep = 1, nsub + ! ------------------------------------------------------------ + CALL pkupper (klon, klev, ptop, ph, pupper, kupper) + + !print*, 'ptop, pupper, ktop, kupper', ptop, pupper, ktop, kupper + + ! 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 + IF (prt_level>=10) THEN + PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', & + isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) + + ENDIF + + ! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement + ! negatif de ktop a kupper -------- + ! cc On calcule pour cela une densite wdens0 pour laquelle on + ! aurait un entrainement nul --- + !jyg< + ! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou + ! les poches sont insuffisantes pour accueillir tout le flux de masse + ! des descentes unsaturees. Nous faisons alors l'hypothese que la + ! convection profonde cree directement de nouvelles poches, sans passer + ! par les thermiques. La nouvelle valeur de wdens est alors imposee. + + 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 + IF ( iflag_wk_profile == 0 ) THEN + omg(i, kupper(i)+1)=-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + & + RG*amup(i, kupper(i)+1)/(1.-sigmaw(i)) + ELSE + omg(i, kupper(i)+1)=0. + ENDIF + wdens0 = (sigmaw(i)/(4.*3.14))* & + ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2) + IF (prt_level >= 10) THEN + print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', & + omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i) + ENDIF + IF (wdens(i)<=wdens0*1.1) THEN + IF (iflag_wk_pop_dyn >= 1) THEN + d_dens_bnd2(i) = d_dens_bnd2(i) + wdens0 - wdens(i) + d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i) + ENDIF + wdens(i) = wdens0 + END IF + END IF + END DO + + DO i = 1, klon + IF (wk_adv(i)) THEN + gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i)) + rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i))) +!jyg< +!! sigmaw(i) = amin1(sigmaw(i), sigmaw_max) + sigmaw_targ = min(sigmaw(i), sigmaw_max) + d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) + d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) +! print *,'XXXX3 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) + sigmaw(i) = sigmaw_targ +!>jyg + END IF + END DO + + IF (iflag_wk_pop_dyn == 1) THEN + + CALL wake_popdyn_1 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, & + dtimesub, gfl, rad_wk, f_shear, drdt_pos, & + d_awdens, d_wdens, d_sigmaw, & + iflag_wk_act, wk_adv, cin, wape, & + drdt, & + d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & + d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & + d_wdens_targ, d_sigmaw_targ) + +! The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0. +! Here, it has to be set to zero. + death_rate(:) = 0. + + ELSEIF (iflag_wk_pop_dyn == 2) THEN + CALL wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, & + sigmaw, wdens, awdens, & !! states variables + gfl, cstar, cin, wape, rad_wk, & + d_sigmaw, d_wdens, d_awdens, & !! tendences + cont_fact, & + d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & + d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & + d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd ) + death_rate(:) = 0. +sigmaw=sigmaw-d_sigmaw +wdens=wdens-d_wdens +awdens=awdens-d_awdens + + ELSEIF (iflag_wk_pop_dyn == 0) THEN + + ! cc nrlmd + + DO i = 1, klon + IF (wk_adv(i)) THEN + ! cc nrlmd Introduction du taux de mortalite 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 + + ENDIF ! (iflag_wk_pop_dyn >= 1) + + + ! 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... + !JYG 161013 Correction : maintenant omg est dimensionne a klev. + 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 + 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 + + IF (prt_level>=10) THEN + PRINT *, 'wake-4.2, omg(igout,k) ', (k,omg(igout,k), k=1,klev) + PRINT *, 'wake-4.2, omgtop(igout), ptop(igout), ktop(igout) ', & + omgtop(igout), ptop(igout), ktop(igout) + ENDIF + + ! ----------------- + ! 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 + IF ( iflag_wk_profile == 0 ) THEN + omg(i, kupper(i)+1) =-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + & + RG*amup(i, kupper(i)+1)/(1.-sigmaw(i)) + ELSE + omg(i, kupper(i)+1) = 0. + ENDIF + 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 (reference)' + ! 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 +!! print *,'omg(igout,k) ', (k,omg(igout,k),k=1,klev) + ! 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 + 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-1 + 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 + DO i = 1, klon + IF (wk_adv(i)) THEN + dp_omgbw(i, klev) = 0. + END IF + 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) +! print *, 'VVVVwake k, the(i,k), dth(i,k), sigmaw(i) ', k, the(i,k), dth(i,k), sigmaw(i) + 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 + +!! IF (prt_level>=10) THEN + IF (prt_level>=10 .and. wk_adv(igout)) THEN + PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,kupper(igout)) + PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,kupper(igout)) + PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,kupper(igout)) + PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,kupper(igout)) + ENDIF + + ! ----------------------------------------------------------------- + DO k = 1, klev-1 + 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=', k, 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=', k, d_deltaqw(i,k) + + ! and increment large scale tendencies + + + + + ! C + ! ----------------------------------------------------------------- + d_tenv(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)) & + -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/ & + (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)- & + rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/ & + (ph(i,k)-ph(i,k+1)) & + -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i,k+1))/ & + (ph(i,k)-ph(i,k+1)) ) + ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN + d_tenv(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 + ! ------------------------------------------------------------------ + + IF (prt_level>=10) THEN + PRINT *, 'wake-4.3, d_deltatw(igout,k) ', (k,d_deltatw(igout,k), k=1,klev) + PRINT *, 'wake-4.3, d_deltaqw(igout,k) ', (k,d_deltaqw(igout,k), k=1,klev) + ENDIF + + ! Increment state variables +!jyg< + IF (iflag_wk_pop_dyn >= 1) THEN + DO k = 1, klev + DO i = 1, klon + IF (wk_adv(i) .AND. k<=kupper(i)) THEN + detr(i,k) = - d_sig_death(i) - d_sig_col(i) + entr(i,k) = d_sig_gen(i) + ENDIF + ENDDO + ENDDO + ELSE ! (iflag_wk_pop_dyn >= 1) + DO k = 1, klev + DO i = 1, klon + IF (wk_adv(i) .AND. k<=kupper(i)) THEN + detr(i, k) = 0. + + entr(i, k) = 0. + ENDIF + ENDDO + ENDDO + ENDIF ! (iflag_wk_pop_dyn >= 1) + + + + 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 repartition + + 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) + +! + + ! cc nrlmd Prise en compte du taux de mortalite + ! cc Definitions de entr, detr +!jyg< +!! detr(i, k) = 0. +!! +!! entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + & +!! sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k) +!! + entr(i, k) = entr(i,k) + gfl(i)*cstar(i) + & + sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k) +!>jyg + wkspread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i) + + ! cc wkspread(i,k) = + ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/ + ! cc $ sigmaw(i) + + + ! ajout d'un effet onde de gravite -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)-wkspread(k)*deltatw(k)) + + ! GW formule 1 + + ! deltatw(k) = deltatw(k)+dtimesub* + ! $ (ff+dtKE(k) - wkspread(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) - & + 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) - & + entr(i,k)*deltatw(i,k)/sigmaw(i) - & + (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) - & + 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) - & + 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_loc, qe, d_qe, deltaqw, & + d_deltaqw, sigmaw, d_sigmaw, alpha) + ! + ! Alpha_tot = Product of all the alpha's + DO i = 1, klon + IF (wk_adv(i)) THEN + alpha_tot(i) = alpha_tot(i)*alpha(i) + END IF + END DO + + ! 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_tenv(i, k) = alpha(i)*d_tenv(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_tenv(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 + tenv(i, k) = tenv0(i, k) + dtls(i, k) + qe(i, k) = qe0(i, k) + dqls(i, k) + the(i, k) = tenv(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) + d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i) +! print *,'XXXX4 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) + END IF + END DO +!jyg< + IF (iflag_wk_pop_dyn >= 1) THEN +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! sigmaw !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Cumulatives + DO i = 1, klon + IF (wk_adv(i)) THEN + d_sig_gen2(i) = d_sig_gen2(i) + d_sig_gen(i) + d_sig_death2(i) = d_sig_death2(i) + d_sig_death(i) + d_sig_col2(i) = d_sig_col2(i) + d_sig_col(i) + d_sig_spread2(i)= d_sig_spread2(i)+ d_sig_spread(i) + d_sig_bnd2(i) = d_sig_bnd2(i) + d_sig_bnd(i) + END IF + END DO +! Bounds + DO i = 1, klon + IF (wk_adv(i)) THEN + sigmaw_targ = max(sigmaw(i),sigmad) + d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) + d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) +! print *,'XXXX5 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) + sigmaw(i) = sigmaw_targ + END IF + END DO +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! wdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Cumulatives + DO i = 1, klon + IF (wk_adv(i)) THEN + wdens(i) = wdens(i) + d_wdens(i) + d_wdens2(i) = d_wdens2(i) + d_wdens(i) + d_dens_gen2(i) = d_dens_gen2(i) + d_dens_gen(i) + d_dens_death2(i) = d_dens_death2(i) + d_dens_death(i) + d_dens_col2(i) = d_dens_col2(i) + d_dens_col(i) + d_dens_bnd2(i) = d_dens_bnd2(i) + d_dens_bnd(i) + END IF + END DO +! Bounds + DO i = 1, klon + IF (wk_adv(i)) THEN + wdens_targ = max(wdens(i),wdensmin) + d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i) + d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i) + wdens(i) = wdens_targ + END IF + END DO +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Cumulatives + DO i = 1, klon + IF (wk_adv(i)) THEN + awdens(i) = awdens(i) + d_awdens(i) + d_awdens2(i) = d_awdens2(i) + d_awdens(i) + END IF + END DO +! Bounds + DO i = 1, klon + IF (wk_adv(i)) THEN + wdens_targ = min( max(awdens(i),0.), wdens(i) ) + d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i) + awdens(i) = wdens_targ + END IF + END DO +! + IF (iflag_wk_pop_dyn == 2) THEN +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens again for iflag_wk_pop_dyn = 2!!!!!! +! Cumulatives + DO i = 1, klon + IF (wk_adv(i)) THEN + d_adens_death2(i) = d_adens_death2(i) + d_adens_death(i) + d_adens_icol2(i) = d_adens_icol2(i) + d_adens_icol(i) + d_adens_acol2(i) = d_adens_acol2(i) + d_adens_acol(i) + d_adens_bnd2(i) = d_adens_bnd2(i) + d_adens_bnd(i) + END IF + END DO +! Bounds + DO i = 1, klon + IF (wk_adv(i)) THEN + wdens_targ = min( max(awdens(i),0.), wdens(i) ) + d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i) + END IF + END DO + ENDIF ! (iflag_wk_pop_dyn == 2) + ENDIF ! (iflag_wk_pop_dyn >= 1) + + + ! 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)=kupper(i)) THEN + deltatw(i, k) = 0. + deltaqw(i, k) = 0. + d_deltatw2(i,k) = -deltatw0(i,k) + d_deltaqw2(i,k) = -deltaqw0(i,k) + 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.+epsim1*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.+epsim1*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)+epsim1*(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. + d_deltatw2(i,k) = -deltatw0(i,k) + d_deltaqw2(i,k) = -deltaqw0(i,k) + 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 +!jyg< +!! sigmaw(i) = max(sigmad, sigd_con(i)) + sigmaw_targ = max(sigmad, sigd_con(i)) + d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) + d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) +! print *,'XXXX6 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) + sigmaw(i) = sigmaw_targ +!>jyg + 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 + + IF (prt_level>=10) THEN + PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', & + sigmaw(igout), cstar(igout), wape(igout), ptop(igout) + ENDIF + + + ! ---------------------------------------------------------- + ! 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_half_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. + + dthmin(i) = -delta_t_min + 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*tenv(i,k)) + IF (k==1) THEN + rhoh(i, k) = ph(i, k)/(RD*tenv(i,k)) + zhh(i, k) = 0 + ELSE + rhoh(i, k) = ph(i, k)*2./(RD*(tenv(i,k)+tenv(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) = tenv(i, k)/ppi(i, k) + thu(i, k) = (tenv(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) + tu(i, k) = tenv(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*(tenv(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. + dz_half(i) = 1. + sum_thvu(i) = thu(i, 1)*(1.+epsim1*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) + dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-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.+epsim1*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) +! + dthmin(i) = min(dthmin(i), dth(i,k)) + END IF + IF (dz_half(i)>0) THEN + sum_half_dth(i) = sum_half_dth(i) + dth(i, k)*dz_half(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)+epsim1*(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 + + IF (iflag_wk_check_trgl>=1) THEN + ! Check triangular shape of dth profile + DO i = 1, klon + IF (ok_qx_qw(i)) THEN + !! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i) + !! print *,'wake, 2.*sum_dth(i)/(hw0(i)*dthmin(i)) ', & + !! 2.*sum_dth(i)/(hw0(i)*dthmin(i)) + !! print *,'wake, sum_half_dth(i), sum_dth(i) ', & + !! sum_half_dth(i), sum_dth(i) + IF ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min) ) THEN + wape2(i) = -1. + !! print *,'wake, rej 1' + ELSE IF (iflag_wk_check_trgl==1.AND.abs(2.*sum_dth(i)/(hw0(i)*dthmin(i)) - 1.) > 0.5) THEN + wape2(i) = -1. + !! print *,'wake, rej 2' + ELSE IF (abs(sum_half_dth(i)) < 0.5*abs(sum_dth(i)) ) THEN + wape2(i) = -1. + !! print *,'wake, rej 3' + END IF + END IF + END DO + END IF + + + 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. + d_deltatw2(i,k) = -deltatw0(i,k) + d_deltaqw2(i,k) = -deltaqw0(i,k) + 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 +!jyg< +!! sigmaw(i) = amax1(sigmad, sigd_con(i)) + sigmaw_targ = max(sigmad, sigd_con(i)) + d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) + d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) +! print *,'XXXX7 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) + sigmaw(i) = sigmaw_targ +!>jyg + 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 + + !jyg< + IF (iflag_wk_pop_dyn >= 1) THEN + DO i = 1, klon + kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. & + .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin) + ENDDO + ELSE ! (iflag_wk_pop_dyn >= 1) + DO i = 1, klon + kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. & + .NOT. ok_qx_qw(i) + ENDDO + ENDIF ! (iflag_wk_pop_dyn >= 1) + !>jyg + + DO k = 1, klev + DO i = 1, klon +!!jyg IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. & +!!jyg .NOT. ok_qx_qw(i)) THEN + IF (kill_wake(i)) THEN + ! cc + dtls(i, k) = 0. + dqls(i, k) = 0. + deltatw(i, k) = 0. + deltaqw(i, k) = 0. + d_deltatw2(i,k) = -deltatw0(i,k) + d_deltaqw2(i,k) = -deltaqw0(i,k) + END IF ! (kill_wake(i)) + END DO + END DO + + DO i = 1, klon +!!jyg IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. & +!!jyg .NOT. ok_qx_qw(i)) THEN + IF (kill_wake(i)) THEN + ktopw(i) = 0 + wape(i) = 0. + cstar(i) = 0. +!!jyg Outside subroutine "Wake" hw, wdens and sigmaw are zero when there are no wakes +!! hw(i) = hwmin !jyg +!! sigmaw(i) = sigmad !jyg + hw(i) = 0. !jyg + fip(i) = 0. +!! sigmaw(i) = 0. !jyg + sigmaw_targ = 0. + d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) +!! d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) + d_sigmaw2(i) = sigmaw_targ - sigmaw_in(i) ! _in = correction jyg 20220124 +! print *,'XXXX8 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) + sigmaw(i) = sigmaw_targ + IF (iflag_wk_pop_dyn >= 1) THEN +!! awdens(i) = 0. +!! wdens(i) = 0. + wdens_targ = 0. + d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i) +!! d_wdens2(i) = wdens_targ - wdens(i) + d_wdens2(i) = wdens_targ - wdens_in(i) ! jyg 20220916 + wdens(i) = wdens_targ + wdens_targ = 0. +!!jyg: bug fix : the d_adens_bnd2 computation must be before the update of awdens. + IF (iflag_wk_pop_dyn == 2) THEN + d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i) + ENDIF ! (iflag_wk_pop_dyn == 2) +!! d_awdens2(i) = wdens_targ - awdens(i) + d_awdens2(i) = wdens_targ - awdens_in(i) ! jyg 20220916 + awdens(i) = wdens_targ +!! IF (iflag_wk_pop_dyn == 2) THEN +!! d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i) +!! ENDIF ! (iflag_wk_pop_dyn == 2) + ENDIF ! (iflag_wk_pop_dyn >= 1) + ELSE ! (kill_wake(i)) + wape(i) = wape2(i) + cstar(i) = cstar2(i) + END IF ! (kill_wake(i)) + ! c print*,'wape wape2 ktopw OK_qx_qw =', + ! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i) + END DO + + IF (prt_level>=10) THEN + PRINT *, 'wake-6, wape wape2 ktopw OK_qx_qw =', & + wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout) + ENDIF + + + ! ----------------------------------------------------------------- + ! 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 +!jyg< +!! IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN + IF (ok_qx_qw(i)) THEN +!>jyg + ! 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 +!jyg< + IF (iflag_wk_pop_dyn >= 1) THEN + DO i = 1, klon + IF (ok_qx_qw(i)) THEN + d_sig_gen2(i) = d_sig_gen2(i)/dtime + d_sig_death2(i) = d_sig_death2(i)/dtime + d_sig_col2(i) = d_sig_col2(i)/dtime + d_sig_spread2(i) = d_sig_spread2(i)/dtime + d_sig_bnd2(i) = d_sig_bnd2(i)/dtime + d_sigmaw2(i) = d_sigmaw2(i)/dtime +! print *,'XXXX9 d_sigmaw2(i), sigmaw(i), dtime ', d_sigmaw2(i), sigmaw(i), dtime +! + d_dens_gen2(i) = d_dens_gen2(i)/dtime + d_dens_death2(i) = d_dens_death2(i)/dtime + d_dens_col2(i) = d_dens_col2(i)/dtime + d_dens_bnd2(i) = d_dens_bnd2(i)/dtime + d_awdens2(i) = d_awdens2(i)/dtime + d_wdens2(i) = d_wdens2(i)/dtime + ENDIF + ENDDO + IF (iflag_wk_pop_dyn == 2) THEN + DO i = 1, klon + IF (ok_qx_qw(i)) THEN + d_adens_death2(i) = d_adens_death2(i)/dtime + d_adens_icol2(i) = d_adens_icol2(i)/dtime + d_adens_acol2(i) = d_adens_acol2(i)/dtime + d_adens_bnd2(i) = d_adens_bnd2(i)/dtime + ENDIF + ENDDO + ENDIF ! (iflag_wk_pop_dyn == 2) + ENDIF ! (iflag_wk_pop_dyn >= 1) + +!>jyg + + RETURN +END SUBROUTINE wake + +SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon_loc, 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_loc + 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_loc + 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 positivite 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 + + + +SUBROUTINE pkupper (klon, klev, ptop, ph, pupper, kupper) + +USE lmdz_wake_ini , ONLY : wk_pupper +IMPLICIT NONE + +INTEGER, INTENT(IN) :: klon,klev +REAL, INTENT(IN), DIMENSION (klon,klev+1) :: ph +REAL, INTENT(IN), DIMENSION (klon) :: ptop +REAL, INTENT(OUT), DIMENSION (klon) :: pupper +INTEGER, INTENT(OUT), DIMENSION (klon) :: kupper +INTEGER :: i,k + + + kupper = 0 + +IF (wk_pupper<1.) THEN + ! 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) = wk_pupper*ph(i, 1) + pupper(i) = max(pupper(i), 45000.) + ! cc Pupper(i) = 60000. + END DO + +ELSE + + DO i=1, klon + ! pupper(i) = wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) + pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-5000.) + END DO +END IF + + ! -5/ Determination de kupper + + DO k = klev, 1, -1 + DO i = 1, klon + IF (ph(i,k+1)jyg+mlt +! +!jyg< + d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i)) +!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) + d_dens_bnd(i) = d_wdens_targ - d_wdens(i) + d_wdens(i) = d_wdens_targ +!! d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i)) +!>jyg + +!jyg+mlt< +!! d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) & +!! + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) & +!! - sigmaw(i)*tau_wk_inv_min )*dtimesub + d_sig_gen(i) = wgen(i)*aa0 +!! print*, 'XXX sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min', & +!! sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min + d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min +!! + + d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos + d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos + d_sig_spread(i) = gfl(i)*cstar(i) + d_sig_gen(i) = d_sig_gen(i)*dtimesub + d_sig_death(i) = d_sig_death(i)*dtimesub + d_sig_col(i) = d_sig_col(i)*dtimesub + d_sig_spread(i) = d_sig_spread(i)*dtimesub + d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i) +!>jyg+mlt +! +!jyg< + d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i)) +!! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i) +!! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i) + d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i) + d_sigmaw(i) = d_sigmaw_targ +!! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i)) +!>jyg + ENDIF + ENDDO + + IF (prt_level >= 10) THEN + print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', & + cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) + print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', & + wdens(1), awdens(1), act(1), d_awdens(1) + print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', & + wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1) + print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', & + d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) + ENDIF + + RETURN + END SUBROUTINE wake_popdyn_1 + + SUBROUTINE wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, & + sigmaw, wdens, awdens, & !! states variables + gfl, cstar, cin, wape, rad_wk, & + d_sigmaw, d_wdens, d_awdens, & !! tendences + cont_fact, & + d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & + d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & + d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd ) + + + + USE lmdz_wake_ini , ONLY : wake_ini + USE lmdz_wake_ini , ONLY : prt_level,RG + USE lmdz_wake_ini , ONLY : stark, wdens_ref + USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0 + USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin + USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max + +IMPLICIT NONE + + INTEGER, INTENT(IN) :: klon,klev + LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv + REAL, INTENT(IN) :: dtimesub + REAL, DIMENSION (klon), INTENT(IN) :: wgen !! B = birth rate of wakes + REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw !! sigma = fractional area of wakes + REAL, DIMENSION (klon), INTENT(INOUT) :: wdens !! D = number of wakes per unit area + REAL, DIMENSION (klon), INTENT(INOUT) :: awdens !! A = number of active wakes per unit area + REAL, DIMENSION (klon), INTENT(IN) :: gfl !! Lg = gust front lenght per unit area + REAL, DIMENSION (klon), INTENT(IN) :: cstar !! C* = spreading velocity of wakes + REAL, DIMENSION (klon), INTENT(IN) :: cin, wape ! RM : A Faire disparaitre + REAL, DIMENSION (klon), INTENT(IN) :: rad_wk !! r = wake radius + +! + REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_wdens, d_awdens + REAL, DIMENSION (klon), INTENT(OUT) :: cont_fact !! RM facteur de contact = 2 pi * rad * C* + ! Some components of the tendencies of state variables + REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd + REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd + REAL, DIMENSION (klon), INTENT(OUT) :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd + + +!! internal variables + + INTEGER :: i, k + REAL, DIMENSION (klon) :: tau_wk_inv !! tau = life time of wakes + REAL :: tau_wk_inv_min + REAL, DIMENSION (klon) :: tau_prime !! tau_prime = life time of actives wakes + REAL :: d_wdens_targ, d_sigmaw_targ + + +!! Equations +!! dD/dt = B - (D-A)/tau - f D^2 +!! dA/dt = B - A/tau_prime + f (D-A)^2 - f A^2 +!! dsigma/dt = B a0 - sigma/D (D-A)/tau + Lg C* - f (D-A)^2 (sigma/D-a0) +!! +!! f = 2 (B (a0-sigma/D) + Lg C*) / (2 (D-A)^2 (2 sigma/D-a0) + D (1-2 sigma)) + + + DO i = 1, klon + ! print*, 'XXX wk_adv(i)', wk_adv(i) + IF (wk_adv(i)) THEN +!! tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.) + tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.) + tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub) + tau_prime(i) = tau_cv +!! cont_fact(i) = 2.*(wgen(i)*(aa0-sigmaw(i)/wdens(i)) + gfl(i)*cstar(i)) / & +!! (2.*(wdens(i)-awdens(i))**2*(2.*sigmaw(i)/wdens(i) - aa0) + wdens(i)*(1.-2.*sigmaw(i))) +!! cont_fact(i) = 2.*3.14*rad_wk(i)*cstar(i) ! bug +!! cont_fact(i) = 4.*3.14*rad_wk(i)*cstar(i) + cont_fact(i) = 2.*gfl(i)*cstar(i)/wdens(i) + + d_sig_gen(i) = wgen(i)*aa0 + d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min + d_sig_col(i) = - cont_fact(i)*(wdens(i)-awdens(i))**2*(2.*sigmaw(i)/wdens(i)-aa0) + d_sig_spread(i) = gfl(i)*cstar(i) +! + d_sig_gen(i) = d_sig_gen(i)*dtimesub + d_sig_death(i) = d_sig_death(i)*dtimesub + d_sig_col(i) = d_sig_col(i)*dtimesub + d_sig_spread(i) = d_sig_spread(i)*dtimesub + d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i) + + + d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i)) +!! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i) +!! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i) + d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i) + d_sigmaw(i) = d_sigmaw_targ +!! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i)) + + + d_dens_gen(i) = wgen(i) + d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min + d_dens_col(i) = - cont_fact(i)*wdens(i)**2 +! + d_dens_gen(i) = d_dens_gen(i)*dtimesub + d_dens_death(i) = d_dens_death(i)*dtimesub + d_dens_col(i) = d_dens_col(i)*dtimesub + d_wdens(i) = d_dens_gen(i) + d_dens_death(i) + d_dens_col(i) + + + d_adens_death(i) = -awdens(i)/tau_prime(i) + d_adens_icol(i) = cont_fact(i)*(wdens(i)-awdens(i))**2 + d_adens_acol(i) = - cont_fact(i)*awdens(i)**2 +! + d_adens_death(i) = d_adens_death(i)*dtimesub + d_adens_icol(i) = d_adens_icol(i)*dtimesub + d_adens_acol(i) = d_adens_acol(i)*dtimesub + d_awdens(i) = d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i) + +!! + d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i)) +!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) + d_dens_bnd(i) = d_wdens_targ - d_wdens(i) + d_wdens(i) = d_wdens_targ + + d_wdens_targ = min(max(d_awdens(i),-awdens(i)), wdens(i)-awdens(i)) +!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) + d_adens_bnd(i) = d_wdens_targ - d_awdens(i) + d_awdens(i) = d_wdens_targ + + + + ENDIF + ENDDO + + IF (prt_level >= 10) THEN + print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1) ', & + cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1) + print *,'wake, wdens(1), awdens(1), d_awdens(1) ', & + wdens(1), awdens(1), d_awdens(1) + print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', & + d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) + ENDIF +sigmaw=sigmaw+d_sigmaw +wdens=wdens+d_wdens +awdens=awdens+d_awdens + + RETURN + END SUBROUTINE wake_popdyn_2 + +END MODULE lmdz_wake Index: LMDZ6/trunk/libf/phylmd/lmdz_wake_ini.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_wake_ini.F90 (revision 4588) +++ LMDZ6/trunk/libf/phylmd/lmdz_wake_ini.F90 (revision 4588) @@ -0,0 +1,219 @@ +MODULE lmdz_wake_ini +IMPLICIT NONE + + + ! ============================================================================ + + + ! But : Decrire le comportement des poches froides apparaissant dans les + ! grands systemes convectifs, et fournir l'energie disponible pour + ! le declenchement de nouvelles colonnes convectives. + + ! State variables : + ! deltatw : temperature difference between wake and off-wake regions + ! deltaqw : specific humidity difference between wake and off-wake regions + ! sigmaw : fractional area covered by wakes. + ! wdens : number of wakes per unit area + + ! ------------------------------------------------------------------------- + ! Declaration de variables + ! ------------------------------------------------------------------------- + + ! Variables a fixer +!jyg< +!! REAL, SAVE :: stark, wdens_ref, coefgw, alpk + INTEGER, SAVE, PROTECTED :: prt_level + REAL, SAVE, PROTECTED, DIMENSION(2) :: wdens_ref + REAL, SAVE, PROTECTED :: stark, coefgw, alpk, wk_pupper +!>jyg + REAL, SAVE, PROTECTED :: crep_upper, crep_sol + !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, wk_pupper, crep_upper, crep_sol) + + REAL, SAVE, PROTECTED :: tau_cv + !$OMP THREADPRIVATE(tau_cv) + REAL, SAVE, PROTECTED :: rzero, aa0 ! minimal wake radius and area + !$OMP THREADPRIVATE(rzero, aa0) + + LOGICAL, SAVE, PROTECTED :: flag_wk_check_trgl + !$OMP THREADPRIVATE(flag_wk_check_trgl) + INTEGER, SAVE, PROTECTED :: iflag_wk_act + !$OMP THREADPRIVATE(iflag_wk_act) + + INTEGER, SAVE, PROTECTED :: iflag_wk_check_trgl + !$OMP THREADPRIVATE(iflag_wk_check_trgl) + INTEGER, SAVE, PROTECTED :: iflag_wk_pop_dyn + !$OMP THREADPRIVATE(iflag_wk_pop_dyn) + + INTEGER, SAVE, PROTECTED :: iflag_wk_profile + !$OMP THREADPRIVATE(iflag_wk_profile) + + REAL, SAVE, PROTECTED :: wdensmin + !$OMP THREADPRIVATE(wdensmin) + REAL, SAVE, PROTECTED :: sigmad, hwmin, wapecut, cstart + !$OMP THREADPRIVATE(sigmad, hwmin, wapecut, cstart) + REAL, SAVE, PROTECTED :: sigmaw_max + !$OMP THREADPRIVATE(sigmaw_max) + REAL, SAVE, PROTECTED :: dens_rate + !$OMP THREADPRIVATE(dens_rate) + REAL, SAVE, PROTECTED :: epsilon_loc + !$OMP THREADPRIVATE(epsilon_loc) + REAL, SAVE, PROTECTED :: epsim1,RG,RD + !$OMP THREADPRIVATE(epsim1,RG,RD) + + + +CONTAINS + + ! ========================================================================= + SUBROUTINE wake_ini(rg_in,rd_in,rv_in,prt_lev) + ! ========================================================================= + + ! ************************************************************** + ! * + ! WAKE * + ! retour a un Pupper fixe * + ! * + ! written by : GRANDPEIX Jean-Yves 09/03/2000 * + ! modified by : ROEHRIG Romain 01/29/2007 * + ! ************************************************************** + + ! ------------------------------------------------------------------------- + ! Initialisations + ! ------------------------------------------------------------------------- + + USE ioipsl_getin_p_mod, ONLY : getin_p + real eps + integer, intent(in) :: prt_lev + real, intent(in) :: rg_in,rd_in,rv_in + + prt_level=prt_lev + epsilon_loc=1.E-15 + wapecut=1. ! previously 5. + ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10. + sigmad=0.02 + hwmin=10. +!! DATA wdensmin/1.e-12/ + wdensmin=1.e-14 + ! cc nrlmd + sigmaw_max=0.4 + dens_rate=0.1 + + eps = rd_in/rv_in + epsim1 = 1.0/eps - 1.0 + RG=rg_in + RD=rd_in + + + ! 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 gravite + ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) + ! wdens : Densite surfacique de poche froide + ! ------------------------------------------------------------------------- + + ! 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 + + + + crep_upper = 0.9 + crep_sol = 1.0 + + ! Get wapecut from parameter file + wapecut = 1. + +print*,'wapecut',wapecut + CALL getin_p('wapecut', wapecut) +print*,'wapecut',wapecut + + ! cc nrlmd Lecture du fichier wake_param.data + + + ! cc nrlmd Lecture du fichier wake_param.data + stark=0.33 + CALL getin_p('stark',stark) + cstart = stark*sqrt(2.*wapecut) + + alpk=0.25 + CALL getin_p('alpk',alpk) + + wk_pupper=0.6 + CALL getin_p('wk_pupper',wk_pupper) + + +!jyg< +!! wdens_ref=8.E-12 +!! CALL getin_p('wdens_ref',wdens_ref) + wdens_ref(1)=8.E-12 + wdens_ref(2)=8.E-12 + CALL getin_p('wdens_ref_o',wdens_ref(1)) !wake number per unit area ; ocean + CALL getin_p('wdens_ref_l',wdens_ref(2)) !wake number per unit area ; land +!>jyg +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!! Population dynamics parameters !!!!!!!!!!!!!!!!!!!!!!!!!!!! +!------------------------------------------------------------------------ + + iflag_wk_pop_dyn = 0 + CALL getin_p('iflag_wk_pop_dyn',iflag_wk_pop_dyn) ! switch between wdens prescribed + ! and wdens prognostic + iflag_wk_act = 0 + CALL getin_p('iflag_wk_act',iflag_wk_act) ! 0: act(:)=0. + ! 1: act(:)=1. + ! 2: act(:)=f(Wape) + + iflag_wk_profile = 0 + CALL getin_p('iflag_wk_profile',iflag_wk_profile) ! switch between wdens prescribed + ! and wdens prognostic + rzero = 5000. + CALL getin_p('rzero_wk', rzero) + aa0 = 3.14*rzero*rzero +! + tau_cv = 4000. + CALL getin_p('tau_cv', tau_cv) + +!------------------------------------------------------------------------ + + coefgw=4. + CALL getin_p('coefgw',coefgw) + + WRITE(*,*) 'stark=', stark + WRITE(*,*) 'alpk=', alpk + WRITE(*,*) 'wk_pupper=', wk_pupper +!jyg< +!! WRITE(*,*) 'wdens_ref=', wdens_ref + WRITE(*,*) 'wdens_ref_o=', wdens_ref(1) + WRITE(*,*) 'wdens_ref_l=', wdens_ref(2) +!>jyg + WRITE(*,*) 'iflag_wk_pop_dyn=',iflag_wk_pop_dyn + WRITE(*,*) 'iflag_wk_act',iflag_wk_act + WRITE(*,*) 'coefgw=', coefgw + + flag_wk_check_trgl=.false. + CALL getin_p('flag_wk_check_trgl ', flag_wk_check_trgl) + WRITE(*,*) 'flag_wk_check_trgl=', flag_wk_check_trgl + WRITE(*,*) 'flag_wk_check_trgl OBSOLETE. Utilisr iflag_wk_check_trgl plutot' + iflag_wk_check_trgl=0 ; IF (flag_wk_check_trgl) iflag_wk_check_trgl=1 + CALL getin_p('iflag_wk_check_trgl ', iflag_wk_check_trgl) + WRITE(*,*) 'iflag_wk_check_trgl=', iflag_wk_check_trgl + + RETURN + +END SUBROUTINE wake_ini + +END MODULE lmdz_wake_ini Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 4586) +++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 4588) @@ -80,5 +80,5 @@ #endif USE lscp_mod, ONLY : lscp - USE wake_ini_mod, ONLY : wake_ini + USE lmdz_wake_ini, ONLY : wake_ini USE yamada_ini_mod, ONLY : yamada_ini USE atke_turbulence_ini_mod, ONLY : atke_ini Index: LMDZ6/trunk/libf/phylmd/wake.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/wake.F90 (revision 4586) +++ (revision ) @@ -1,2667 +1,0 @@ - -! $Id$ - - -SUBROUTINE wake(klon,klev,znatsurf, p, ph, pi, dtime, & - tenv0, qe0, omgb, & - dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, & - sigd_con, Cin, & - deltatw, deltaqw, sigmaw, awdens, wdens, & ! state variables - dth, hw, wape, fip, gfl, & - dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, & - dtke, dqke, omg, dp_deltomg, wkspread, cstar, & - d_deltat_gw, & ! tendencies - d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2) ! tendencies - - - ! ************************************************************** - ! * - ! WAKE * - ! retour a un Pupper fixe * - ! * - ! written by : GRANDPEIX Jean-Yves 09/03/2000 * - ! modified by : ROEHRIG Romain 01/29/2007 * - ! ************************************************************** - - - USE wake_ini_mod , ONLY : wake_ini - USE wake_ini_mod , ONLY : prt_level,epsim1,RG,RD - USE wake_ini_mod , ONLY : stark, wdens_ref, coefgw, alpk, wk_pupper - USE wake_ini_mod , ONLY : crep_upper, crep_sol, tau_cv, rzero, aa0, flag_wk_check_trgl - USE wake_ini_mod , ONLY : iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensmin - USE wake_ini_mod , ONLY : sigmad, hwmin, wapecut, cstart, sigmaw_max, dens_rate, epsilon_loc - USE wake_ini_mod , ONLY : iflag_wk_profile - - - IMPLICIT NONE - ! ============================================================================ - - - ! But : Decrire le comportement des poches froides apparaissant dans les - ! grands systemes convectifs, et fournir l'energie disponible pour - ! le declenchement de nouvelles colonnes convectives. - - ! State variables : - ! deltatw : temperature difference between wake and off-wake regions - ! deltaqw : specific humidity difference between wake and off-wake regions - ! sigmaw : fractional area covered by wakes. - ! wdens : number of wakes per unit area - - ! Variable de sortie : - - ! wape : WAke Potential Energy - ! fip : Front Incident Power (W/m2) - ALP - ! gfl : Gust Front Length per unit area (m-1) - ! dtls : large scale temperature tendency due to wake - ! dqls : large scale humidity tendency due to wake - ! hw : wake top hight (given by hw*deltatw(1)/2=wape) - ! dp_omgb : vertical gradient of large scale omega - ! awdens : densite de poches actives - ! wdens : densite de poches - ! omgbdth: flux of Delta_Theta transported by LS omega - ! dtKE : differential heating (wake - unpertubed) - ! dqKE : differential moistening (wake - unpertubed) - ! omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s) - ! dp_deltomg : vertical gradient of omg (s-1) - ! wkspread : spreading term in d_t_wake and d_q_wake - ! deltatw : updated temperature difference (T_w-T_u). - ! deltaqw : updated humidity difference (q_w-q_u). - ! sigmaw : updated wake fractional area. - ! d_deltat_gw : delta T tendency due to GW - - ! Variables d'entree : - - ! aire : aire de la maille - ! tenv0 : temperature dans l'environnement (K) - ! qe0 : humidite dans l'environnement (kg/kg) - ! omgb : vitesse verticale moyenne sur la maille (Pa/s) - ! dtdwn: source de chaleur due aux descentes (K/s) - ! dqdwn: source d'humidite due aux descentes (kg/kg/s) - ! dta : source de chaleur due courants satures et detrain (K/s) - ! dqa : source d'humidite due aux courants satures et detra (kg/kg/s) - ! wgen : number of wakes generated per unit area and per sec (/m^2/s) - ! amdwn: flux de masse total des descentes, par unite de - ! surface de la maille (kg/m2/s) - ! amup : flux de masse total des ascendances, par unite de - ! surface de la maille (kg/m2/s) - ! sigd_con: - ! Cin : convective inhibition - ! p : pressions aux milieux des couches (Pa) - ! ph : pressions aux interfaces (Pa) - ! pi : (p/p_0)**kapa (adim) - ! dtime: increment temporel (s) - - ! Variables internes : - - ! rhow : masse volumique de la poche froide - ! rho : environment density at P levels - ! rhoh : environment density at Ph levels - ! tenv : environment temperature | may change within - ! qe : environment humidity | sub-time-stepping - ! the : environment potential temperature - ! thu : potential temperature in undisturbed area - ! tu : temperature in undisturbed area - ! qu : humidity in undisturbed area - ! dp_omgb: vertical gradient og LS omega - ! omgbw : wake average vertical omega - ! dp_omgbw: vertical gradient of omgbw - ! omgbdq : flux of Delta_q transported by LS omega - ! dth : potential temperature diff. wake-undist. - ! th1 : first pot. temp. for vertical advection (=thu) - ! th2 : second pot. temp. for vertical advection (=thw) - ! q1 : first humidity for vertical advection - ! q2 : second humidity for vertical advection - ! d_deltatw : terme de redistribution pour deltatw - ! d_deltaqw : terme de redistribution pour deltaqw - ! deltatw0 : deltatw initial - ! deltaqw0 : deltaqw initial - ! hw0 : wake top hight (defined as the altitude at which deltatw=0) - ! amflux : horizontal mass flux through wake boundary - ! wdens_ref: initial number of wakes per unit area (3D) or per - ! unit length (2D), at the beginning of each time step - ! Tgw : 1 sur la periode de onde de gravite - ! Cgw : vitesse de propagation de onde de gravite - ! LL : distance entre 2 poches - - ! ------------------------------------------------------------------------- - ! Declaration de variables - ! ------------------------------------------------------------------------- - - - ! Arguments en entree - ! -------------------- - - INTEGER, INTENT(IN) :: klon,klev - INTEGER, DIMENSION (klon), INTENT(IN) :: znatsurf - REAL, DIMENSION (klon, klev), INTENT(IN) :: p, pi - REAL, DIMENSION (klon, klev+1), INTENT(IN) :: ph - REAL, DIMENSION (klon, klev), INTENT(IN) :: omgb - REAL, INTENT(IN) :: dtime - REAL, DIMENSION (klon, klev), INTENT(IN) :: tenv0, qe0 - REAL, DIMENSION (klon, klev), INTENT(IN) :: dtdwn, dqdwn - REAL, DIMENSION (klon, klev), INTENT(IN) :: amdwn, amup - REAL, DIMENSION (klon, klev), INTENT(IN) :: dta, dqa - REAL, DIMENSION (klon), INTENT(IN) :: wgen - REAL, DIMENSION (klon), INTENT(IN) :: sigd_con - REAL, DIMENSION (klon), INTENT(IN) :: Cin - - ! - ! Input/Output - ! State variables - REAL, DIMENSION (klon, klev), INTENT(INOUT) :: deltatw, deltaqw - REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw - REAL, DIMENSION (klon), INTENT(INOUT) :: awdens - REAL, DIMENSION (klon), INTENT(INOUT) :: wdens - - ! 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) :: wkspread ! unused (jyg) - REAL, DIMENSION (klon, klev), INTENT(OUT) :: omgbdth, omg - REAL, DIMENSION (klon, klev), INTENT(OUT) :: dp_omgb, dp_deltomg - REAL, DIMENSION (klon), INTENT(OUT) :: hw, wape, fip, gfl, cstar - INTEGER, DIMENSION (klon), INTENT(OUT) :: ktopw - ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields - ! computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens) - REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_gw - REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltatw2, d_deltaqw2 - REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw2, d_awdens2, d_wdens2 - - ! Variables internes - ! ------------------- - - ! Variables a fixer - - REAL :: delta_t_min - INTEGER :: nsub - REAL :: dtimesub - 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) :: tenv, qe -!! REAL, DIMENSION (klon) :: sigmaw1 - - ! Variables liees a la dynamique de population 1 - REAL, DIMENSION(klon) :: act - REAL, DIMENSION(klon) :: rad_wk, tau_wk_inv - REAL, DIMENSION(klon) :: f_shear - REAL, DIMENSION(klon) :: drdt - - ! Variables liees a la dynamique de population 2 - REAL, DIMENSION(klon) :: cont_fact - -!! REAL, DIMENSION(klon) :: d_sig_gen, d_sig_death, d_sig_col - REAL, DIMENSION(klon) :: wape1_act, wape2_act - LOGICAL, DIMENSION (klon) :: kill_wake - REAL :: drdt_pos - REAL :: tau_wk_inv_min - ! Some components of the tendencies of state variables - REAL, DIMENSION (klon) :: d_sig_gen2, d_sig_death2, d_sig_col2, d_sig_spread2, d_sig_bnd2 - REAL, DIMENSION (klon) :: d_dens_gen2, d_dens_death2, d_dens_col2, d_dens_bnd2 - REAL, DIMENSION (klon) :: d_adens_death2, d_adens_icol2, d_adens_acol2, d_adens_bnd2 - - ! Variables pour les GW - REAL, DIMENSION (klon) :: ll - REAL, DIMENSION (klon, klev) :: n2 - REAL, DIMENSION (klon, klev) :: cgw - REAL, DIMENSION (klon, klev) :: tgw - - ! Variables liees 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 - - ! Variables liees au test de la forme triangulaire du profil de Delta_theta - REAL, DIMENSION (klon) :: sum_half_dth - REAL, DIMENSION (klon) :: dz_half - - ! Sub-timestep tendencies and related variables - REAL, DIMENSION (klon, klev) :: d_deltatw, d_deltaqw - REAL, DIMENSION (klon, klev) :: d_tenv, d_qe - REAL, DIMENSION (klon) :: d_awdens, d_wdens, d_sigmaw - REAL, DIMENSION (klon) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd - REAL, DIMENSION (klon) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd - REAL, DIMENSION (klon) :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd - REAL, DIMENSION (klon) :: alpha, alpha_tot - REAL, DIMENSION (klon) :: q0_min, q1_min - LOGICAL, DIMENSION (klon) :: wk_adv, ok_qx_qw - - ! Autres variables internes - INTEGER ::isubstep, k, i, igout - - REAL :: sigmaw_targ - REAL :: wdens_targ - REAL :: d_sigmaw_targ - REAL :: d_wdens_targ - - 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) :: 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 -!! REAL, DIMENSION (klon) :: nat_rate - REAL, DIMENSION (klon, klev) :: entr - REAL, DIMENSION (klon, klev) :: detr - - REAL, DIMENSION(klon) :: sigmaw_in ! pour les prints - REAL, DIMENSION(klon) :: awdens_in, wdens_in ! pour les prints - - ! ------------------------------------------------------------------------- - ! Initialisations - ! ------------------------------------------------------------------------- - ! ALON = 3.e5 - ! alon = 1.E6 - - ! Provisionnal; to be suppressed when f_shear is parameterized - f_shear(:) = 1. ! 0. for strong shear, 1. for weak shear - - - ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei) - - ! coefgw : Coefficient pour les ondes de gravite - ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) - ! wdens : Densite surfacique de poche froide - ! ------------------------------------------------------------------------- - - ! 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 -!print *,'XXXX dtime input ', dtime - igout = klon/2+1/klon - - IF (iflag_wk_pop_dyn == 0) THEN - ! 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) - !jyg< - !! wdens(:) = wdens_ref - DO i = 1,klon - wdens(i) = wdens_ref(znatsurf(i)+1) - ENDDO - !>jyg - ENDIF ! (iflag_wk_pop_dyn == 0) - - ! 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, initialize tendencies, initialize output fields - ! ------------------------------------------------------------------------ - -!jyg< -!! 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) -!! tenv(i, k) = tenv0(i, k) -!! qe(i, k) = qe0(i, k) -!! dtls(i, k) = 0. -!! dqls(i, k) = 0. -!! d_deltat_gw(i, k) = 0. -!! d_tenv(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 - ppi(:,:) = pi(:,:) - deltatw0(:,:) = deltatw(:,:) - deltaqw0(:,:) = deltaqw(:,:) - tenv(:,:) = tenv0(:,:) - qe(:,:) = qe0(:,:) - dtls(:,:) = 0. - dqls(:,:) = 0. - d_deltat_gw(:,:) = 0. - d_tenv(:,:) = 0. - d_qe(:,:) = 0. - d_deltatw(:,:) = 0. - d_deltaqw(:,:) = 0. - d_deltatw2(:,:) = 0. - d_deltaqw2(:,:) = 0. - - IF (iflag_wk_act == 0) THEN - act(:) = 0. - ELSEIF (iflag_wk_act == 1) THEN - act(:) = 1. - ENDIF - -!! DO i = 1, klon -!! sigmaw_in(i) = sigmaw(i) -!! END DO - sigmaw_in(:) = sigmaw(:) -!>jyg -! - IF (iflag_wk_pop_dyn >= 1) THEN - awdens_in(:) = awdens(:) - wdens_in(:) = wdens(:) -!! wdens(:) = wdens(:) + wgen(:)*dtime -!! d_wdens2(:) = wgen(:)*dtime -!! ELSE - ENDIF ! (iflag_wk_pop_dyn >= 1) - - - ! sigmaw1=sigmaw - ! IF (sigd_con.GT.sigmaw1) THEN - ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con - ! ENDIF - IF (iflag_wk_pop_dyn >= 1) THEN - DO i = 1, klon - d_dens_gen2(i) = 0. - d_dens_death2(i) = 0. - d_dens_col2(i) = 0. - d_awdens2(i) = 0. -! - wdens_targ = max(wdens(i),wdensmin) - d_dens_bnd2(i) = wdens_targ - wdens(i) - d_wdens2(i) = wdens_targ - wdens(i) - wdens(i) = wdens_targ - END DO - IF (iflag_wk_pop_dyn == 2) THEN - DO i = 1, klon - d_adens_death2(i) = 0. - d_adens_icol2(i) = 0. - d_adens_acol2(i) = 0. -! - wdens_targ = min(max(awdens(i),0.),wdens(i)) - d_adens_bnd2(i) = wdens_targ - awdens(i) - d_awdens2(i) = wdens_targ - awdens(i) - awdens(i) = wdens_targ - END DO - ENDIF ! (iflag_wk_pop_dyn == 2) - ELSE - DO i = 1, klon - d_awdens2(i) = 0. - d_wdens2(i) = 0. - END DO - ENDIF ! (iflag_wk_pop_dyn >= 1) -! - DO i = 1, klon - ! c sigmaw(i) = amax1(sigmaw(i),sigd_con(i)) -!jyg< -!! sigmaw(i) = amax1(sigmaw(i), sigmad) -!! sigmaw(i) = amin1(sigmaw(i), 0.99) - d_sig_gen2(i) = 0. - d_sig_death2(i) = 0. - d_sig_col2(i) = 0. - d_sig_spread2(i)= 0. - sigmaw_targ = min(max(sigmaw(i), sigmad),0.99) - d_sig_bnd2(i) = sigmaw_targ - sigmaw(i) - d_sigmaw2(i) = sigmaw_targ - sigmaw(i) -! print *,'XXXX1 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) - sigmaw(i) = sigmaw_targ -!>jyg - END DO - - wape(:) = 0. - wape2(:) = 0. - d_sigmaw(:) = 0. - ktopw(:) = 0 -! -!jyg -! - IF (prt_level>=10) THEN - PRINT *, 'wake-1, sigmaw(igout) ', sigmaw(igout) - PRINT *, 'wake-1, deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev) - PRINT *, 'wake-1, deltaqw(igout,k) ', (k,deltaqw(igout,k), k=1,klev) - PRINT *, 'wake-1, dowwdraughts, amdwn(igout,k) ', (k,amdwn(igout,k), k=1,klev) - PRINT *, 'wake-1, dowwdraughts, dtdwn(igout,k) ', (k,dtdwn(igout,k), k=1,klev) - PRINT *, 'wake-1, dowwdraughts, dqdwn(igout,k) ', (k,dqdwn(igout,k), k=1,klev) - PRINT *, 'wake-1, updraughts, amup(igout,k) ', (k,amup(igout,k), k=1,klev) - PRINT *, 'wake-1, updraughts, dta(igout,k) ', (k,dta(igout,k), k=1,klev) - PRINT *, 'wake-1, updraughts, dqa(igout,k) ', (k,dqa(igout,k), k=1,klev) - ENDIF - - ! 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,tenv(i,k) - rho(i, k) = p(i, k)/(RD*tenv(i,k)) - ! write(*,*)'wake 2',rho(i,k) - IF (k==1) THEN - ! write(*,*)'wake 3',i,k,rd,tenv(i,k) - rhoh(i, k) = ph(i, k)/(RD*tenv(i,k)) - ! write(*,*)'wake 4',i,k,rd,tenv(i,k) - zhh(i, k) = 0 - ELSE - ! write(*,*)'wake 5',rd,(tenv(i,k)+tenv(i,k-1)) - rhoh(i, k) = ph(i, k)*2./(RD*(tenv(i,k)+tenv(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) = tenv(i, k)/ppi(i, k) - thu(i, k) = (tenv(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) - tu(i, k) = tenv(i, k) - deltatw(i, k)*sigmaw(i) - qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) - ! write(*,*)'wake 8',(RD*(tenv(i,k)+deltatw(i,k))) - rhow(i, k) = p(i, k)/(RD*(tenv(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 - ! ----------------------------------------------------------------- - - ! 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 - - IF (prt_level>=10) THEN - PRINT *, 'wake-2, ptop_provis(igout), ptop(igout) ', ptop_provis(igout), ptop(igout) - ENDIF - - - ! -5/ Determination de ktop et kupper - - CALL pkupper (klon, klev, ptop, ph, pupper, kupper) - - DO k = klev, 1, -1 - DO i = 1, klon - IF (ph(i,k+1)-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)=10) THEN - PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout) - ENDIF - - ! -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. - d_deltatw2(i,k) = -deltatw0(i,k) - d_deltaqw2(i,k) = -deltaqw0(i,k) - 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.+epsim1*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.+epsim1*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)+ & - epsim1*(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. - d_deltatw2(i,k) = -deltatw0(i,k) - d_deltaqw2(i,k) = -deltaqw0(i,k) - END IF - END DO - END DO - - DO i = 1, klon - IF (wape(i)<0.) THEN - wape(i) = 0. - cstar(i) = 0. - hw(i) = hwmin -!jyg< -!! sigmaw(i) = amax1(sigmad, sigd_con(i)) - sigmaw_targ = max(sigmad, sigd_con(i)) - d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) - d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) -! print *,'XXXX2 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) - sigmaw(i) = sigmaw_targ -!>jyg - fip(i) = 0. - gwake(i) = .FALSE. - ELSE - hw(i) = hw0(i) - 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. - alpha_tot(i) = 1. - END DO - - IF (prt_level>=10) THEN - PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ', & - sigmaw(igout), cstar(igout), wape(igout), ktop(igout) - ENDIF - - - ! C ----------------------------------------------------------------- - ! Sub-time-stepping - ! ----------------- - - nsub = 10 - dtimesub = dtime/nsub - - - - ! ------------------------------------------------------------ - DO isubstep = 1, nsub - ! ------------------------------------------------------------ - CALL pkupper (klon, klev, ptop, ph, pupper, kupper) - - !print*, 'ptop, pupper, ktop, kupper', ptop, pupper, ktop, kupper - - ! 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 - IF (prt_level>=10) THEN - PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', & - isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) - - ENDIF - - ! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement - ! negatif de ktop a kupper -------- - ! cc On calcule pour cela une densite wdens0 pour laquelle on - ! aurait un entrainement nul --- - !jyg< - ! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou - ! les poches sont insuffisantes pour accueillir tout le flux de masse - ! des descentes unsaturees. Nous faisons alors l'hypothese que la - ! convection profonde cree directement de nouvelles poches, sans passer - ! par les thermiques. La nouvelle valeur de wdens est alors imposee. - - 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 - IF ( iflag_wk_profile == 0 ) THEN - omg(i, kupper(i)+1)=-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + & - RG*amup(i, kupper(i)+1)/(1.-sigmaw(i)) - ELSE - omg(i, kupper(i)+1)=0. - ENDIF - wdens0 = (sigmaw(i)/(4.*3.14))* & - ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2) - IF (prt_level >= 10) THEN - print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', & - omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i) - ENDIF - IF (wdens(i)<=wdens0*1.1) THEN - IF (iflag_wk_pop_dyn >= 1) THEN - d_dens_bnd2(i) = d_dens_bnd2(i) + wdens0 - wdens(i) - d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i) - ENDIF - wdens(i) = wdens0 - END IF - END IF - END DO - - DO i = 1, klon - IF (wk_adv(i)) THEN - gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i)) - rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i))) -!jyg< -!! sigmaw(i) = amin1(sigmaw(i), sigmaw_max) - sigmaw_targ = min(sigmaw(i), sigmaw_max) - d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) - d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) -! print *,'XXXX3 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) - sigmaw(i) = sigmaw_targ -!>jyg - END IF - END DO - - IF (iflag_wk_pop_dyn == 1) THEN - - CALL wake_popdyn_1 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, & - dtimesub, gfl, rad_wk, f_shear, drdt_pos, & - d_awdens, d_wdens, d_sigmaw, & - iflag_wk_act, wk_adv, cin, wape, & - drdt, & - d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & - d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & - d_wdens_targ, d_sigmaw_targ) - -! The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0. -! Here, it has to be set to zero. - death_rate(:) = 0. - - ELSEIF (iflag_wk_pop_dyn == 2) THEN - CALL wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, & - sigmaw, wdens, awdens, & !! states variables - gfl, cstar, cin, wape, rad_wk, & - d_sigmaw, d_wdens, d_awdens, & !! tendences - cont_fact, & - d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & - d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & - d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd ) - death_rate(:) = 0. -sigmaw=sigmaw-d_sigmaw -wdens=wdens-d_wdens -awdens=awdens-d_awdens - - ELSEIF (iflag_wk_pop_dyn == 0) THEN - - ! cc nrlmd - - DO i = 1, klon - IF (wk_adv(i)) THEN - ! cc nrlmd Introduction du taux de mortalite 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 - - ENDIF ! (iflag_wk_pop_dyn >= 1) - - - ! 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... - !JYG 161013 Correction : maintenant omg est dimensionne a klev. - 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 - 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 - - IF (prt_level>=10) THEN - PRINT *, 'wake-4.2, omg(igout,k) ', (k,omg(igout,k), k=1,klev) - PRINT *, 'wake-4.2, omgtop(igout), ptop(igout), ktop(igout) ', & - omgtop(igout), ptop(igout), ktop(igout) - ENDIF - - ! ----------------- - ! 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 - IF ( iflag_wk_profile == 0 ) THEN - omg(i, kupper(i)+1) =-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + & - RG*amup(i, kupper(i)+1)/(1.-sigmaw(i)) - ELSE - omg(i, kupper(i)+1) = 0. - ENDIF - 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 (reference)' - ! 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 -!! print *,'omg(igout,k) ', (k,omg(igout,k),k=1,klev) - ! 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 - 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-1 - 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 - DO i = 1, klon - IF (wk_adv(i)) THEN - dp_omgbw(i, klev) = 0. - END IF - 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) -! print *, 'VVVVwake k, the(i,k), dth(i,k), sigmaw(i) ', k, the(i,k), dth(i,k), sigmaw(i) - 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 - -!! IF (prt_level>=10) THEN - IF (prt_level>=10 .and. wk_adv(igout)) THEN - PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,kupper(igout)) - PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,kupper(igout)) - PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,kupper(igout)) - PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,kupper(igout)) - ENDIF - - ! ----------------------------------------------------------------- - DO k = 1, klev-1 - 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=', k, 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=', k, d_deltaqw(i,k) - - ! and increment large scale tendencies - - - - - ! C - ! ----------------------------------------------------------------- - d_tenv(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)) & - -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/ & - (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)- & - rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/ & - (ph(i,k)-ph(i,k+1)) & - -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i,k+1))/ & - (ph(i,k)-ph(i,k+1)) ) - ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN - d_tenv(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 - ! ------------------------------------------------------------------ - - IF (prt_level>=10) THEN - PRINT *, 'wake-4.3, d_deltatw(igout,k) ', (k,d_deltatw(igout,k), k=1,klev) - PRINT *, 'wake-4.3, d_deltaqw(igout,k) ', (k,d_deltaqw(igout,k), k=1,klev) - ENDIF - - ! Increment state variables -!jyg< - IF (iflag_wk_pop_dyn >= 1) THEN - DO k = 1, klev - DO i = 1, klon - IF (wk_adv(i) .AND. k<=kupper(i)) THEN - detr(i,k) = - d_sig_death(i) - d_sig_col(i) - entr(i,k) = d_sig_gen(i) - ENDIF - ENDDO - ENDDO - ELSE ! (iflag_wk_pop_dyn >= 1) - DO k = 1, klev - DO i = 1, klon - IF (wk_adv(i) .AND. k<=kupper(i)) THEN - detr(i, k) = 0. - - entr(i, k) = 0. - ENDIF - ENDDO - ENDDO - ENDIF ! (iflag_wk_pop_dyn >= 1) - - - - 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 repartition - - 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) - -! - - ! cc nrlmd Prise en compte du taux de mortalite - ! cc Definitions de entr, detr -!jyg< -!! detr(i, k) = 0. -!! -!! entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + & -!! sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k) -!! - entr(i, k) = entr(i,k) + gfl(i)*cstar(i) + & - sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k) -!>jyg - wkspread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i) - - ! cc wkspread(i,k) = - ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/ - ! cc $ sigmaw(i) - - - ! ajout d'un effet onde de gravite -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)-wkspread(k)*deltatw(k)) - - ! GW formule 1 - - ! deltatw(k) = deltatw(k)+dtimesub* - ! $ (ff+dtKE(k) - wkspread(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) - & - 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) - & - entr(i,k)*deltatw(i,k)/sigmaw(i) - & - (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) - & - 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) - & - 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_loc, qe, d_qe, deltaqw, & - d_deltaqw, sigmaw, d_sigmaw, alpha) - ! - ! Alpha_tot = Product of all the alpha's - DO i = 1, klon - IF (wk_adv(i)) THEN - alpha_tot(i) = alpha_tot(i)*alpha(i) - END IF - END DO - - ! 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_tenv(i, k) = alpha(i)*d_tenv(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_tenv(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 - tenv(i, k) = tenv0(i, k) + dtls(i, k) - qe(i, k) = qe0(i, k) + dqls(i, k) - the(i, k) = tenv(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) - d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i) -! print *,'XXXX4 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) - END IF - END DO -!jyg< - IF (iflag_wk_pop_dyn >= 1) THEN -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! sigmaw !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! Cumulatives - DO i = 1, klon - IF (wk_adv(i)) THEN - d_sig_gen2(i) = d_sig_gen2(i) + d_sig_gen(i) - d_sig_death2(i) = d_sig_death2(i) + d_sig_death(i) - d_sig_col2(i) = d_sig_col2(i) + d_sig_col(i) - d_sig_spread2(i)= d_sig_spread2(i)+ d_sig_spread(i) - d_sig_bnd2(i) = d_sig_bnd2(i) + d_sig_bnd(i) - END IF - END DO -! Bounds - DO i = 1, klon - IF (wk_adv(i)) THEN - sigmaw_targ = max(sigmaw(i),sigmad) - d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) - d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) -! print *,'XXXX5 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) - sigmaw(i) = sigmaw_targ - END IF - END DO -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! wdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! Cumulatives - DO i = 1, klon - IF (wk_adv(i)) THEN - wdens(i) = wdens(i) + d_wdens(i) - d_wdens2(i) = d_wdens2(i) + d_wdens(i) - d_dens_gen2(i) = d_dens_gen2(i) + d_dens_gen(i) - d_dens_death2(i) = d_dens_death2(i) + d_dens_death(i) - d_dens_col2(i) = d_dens_col2(i) + d_dens_col(i) - d_dens_bnd2(i) = d_dens_bnd2(i) + d_dens_bnd(i) - END IF - END DO -! Bounds - DO i = 1, klon - IF (wk_adv(i)) THEN - wdens_targ = max(wdens(i),wdensmin) - d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i) - d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i) - wdens(i) = wdens_targ - END IF - END DO -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! Cumulatives - DO i = 1, klon - IF (wk_adv(i)) THEN - awdens(i) = awdens(i) + d_awdens(i) - d_awdens2(i) = d_awdens2(i) + d_awdens(i) - END IF - END DO -! Bounds - DO i = 1, klon - IF (wk_adv(i)) THEN - wdens_targ = min( max(awdens(i),0.), wdens(i) ) - d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i) - awdens(i) = wdens_targ - END IF - END DO -! - IF (iflag_wk_pop_dyn == 2) THEN -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens again for iflag_wk_pop_dyn = 2!!!!!! -! Cumulatives - DO i = 1, klon - IF (wk_adv(i)) THEN - d_adens_death2(i) = d_adens_death2(i) + d_adens_death(i) - d_adens_icol2(i) = d_adens_icol2(i) + d_adens_icol(i) - d_adens_acol2(i) = d_adens_acol2(i) + d_adens_acol(i) - d_adens_bnd2(i) = d_adens_bnd2(i) + d_adens_bnd(i) - END IF - END DO -! Bounds - DO i = 1, klon - IF (wk_adv(i)) THEN - wdens_targ = min( max(awdens(i),0.), wdens(i) ) - d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i) - END IF - END DO - ENDIF ! (iflag_wk_pop_dyn == 2) - ENDIF ! (iflag_wk_pop_dyn >= 1) - - - ! 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)=kupper(i)) THEN - deltatw(i, k) = 0. - deltaqw(i, k) = 0. - d_deltatw2(i,k) = -deltatw0(i,k) - d_deltaqw2(i,k) = -deltaqw0(i,k) - 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.+epsim1*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.+epsim1*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)+epsim1*(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. - d_deltatw2(i,k) = -deltatw0(i,k) - d_deltaqw2(i,k) = -deltaqw0(i,k) - 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 -!jyg< -!! sigmaw(i) = max(sigmad, sigd_con(i)) - sigmaw_targ = max(sigmad, sigd_con(i)) - d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) - d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) -! print *,'XXXX6 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) - sigmaw(i) = sigmaw_targ -!>jyg - 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 - - IF (prt_level>=10) THEN - PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', & - sigmaw(igout), cstar(igout), wape(igout), ptop(igout) - ENDIF - - - ! ---------------------------------------------------------- - ! 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_half_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. - - dthmin(i) = -delta_t_min - 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*tenv(i,k)) - IF (k==1) THEN - rhoh(i, k) = ph(i, k)/(RD*tenv(i,k)) - zhh(i, k) = 0 - ELSE - rhoh(i, k) = ph(i, k)*2./(RD*(tenv(i,k)+tenv(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) = tenv(i, k)/ppi(i, k) - thu(i, k) = (tenv(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) - tu(i, k) = tenv(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*(tenv(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. - dz_half(i) = 1. - sum_thvu(i) = thu(i, 1)*(1.+epsim1*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) - dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-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.+epsim1*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) -! - dthmin(i) = min(dthmin(i), dth(i,k)) - END IF - IF (dz_half(i)>0) THEN - sum_half_dth(i) = sum_half_dth(i) + dth(i, k)*dz_half(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)+epsim1*(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 - - IF (iflag_wk_check_trgl>=1) THEN - ! Check triangular shape of dth profile - DO i = 1, klon - IF (ok_qx_qw(i)) THEN - !! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i) - !! print *,'wake, 2.*sum_dth(i)/(hw0(i)*dthmin(i)) ', & - !! 2.*sum_dth(i)/(hw0(i)*dthmin(i)) - !! print *,'wake, sum_half_dth(i), sum_dth(i) ', & - !! sum_half_dth(i), sum_dth(i) - IF ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min) ) THEN - wape2(i) = -1. - !! print *,'wake, rej 1' - ELSE IF (iflag_wk_check_trgl==1.AND.abs(2.*sum_dth(i)/(hw0(i)*dthmin(i)) - 1.) > 0.5) THEN - wape2(i) = -1. - !! print *,'wake, rej 2' - ELSE IF (abs(sum_half_dth(i)) < 0.5*abs(sum_dth(i)) ) THEN - wape2(i) = -1. - !! print *,'wake, rej 3' - END IF - END IF - END DO - END IF - - - 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. - d_deltatw2(i,k) = -deltatw0(i,k) - d_deltaqw2(i,k) = -deltaqw0(i,k) - 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 -!jyg< -!! sigmaw(i) = amax1(sigmad, sigd_con(i)) - sigmaw_targ = max(sigmad, sigd_con(i)) - d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) - d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) -! print *,'XXXX7 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) - sigmaw(i) = sigmaw_targ -!>jyg - 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 - - !jyg< - IF (iflag_wk_pop_dyn >= 1) THEN - DO i = 1, klon - kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. & - .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin) - ENDDO - ELSE ! (iflag_wk_pop_dyn >= 1) - DO i = 1, klon - kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. & - .NOT. ok_qx_qw(i) - ENDDO - ENDIF ! (iflag_wk_pop_dyn >= 1) - !>jyg - - DO k = 1, klev - DO i = 1, klon -!!jyg IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. & -!!jyg .NOT. ok_qx_qw(i)) THEN - IF (kill_wake(i)) THEN - ! cc - dtls(i, k) = 0. - dqls(i, k) = 0. - deltatw(i, k) = 0. - deltaqw(i, k) = 0. - d_deltatw2(i,k) = -deltatw0(i,k) - d_deltaqw2(i,k) = -deltaqw0(i,k) - END IF ! (kill_wake(i)) - END DO - END DO - - DO i = 1, klon -!!jyg IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. & -!!jyg .NOT. ok_qx_qw(i)) THEN - IF (kill_wake(i)) THEN - ktopw(i) = 0 - wape(i) = 0. - cstar(i) = 0. -!!jyg Outside subroutine "Wake" hw, wdens and sigmaw are zero when there are no wakes -!! hw(i) = hwmin !jyg -!! sigmaw(i) = sigmad !jyg - hw(i) = 0. !jyg - fip(i) = 0. -!! sigmaw(i) = 0. !jyg - sigmaw_targ = 0. - d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i) -!! d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i) - d_sigmaw2(i) = sigmaw_targ - sigmaw_in(i) ! _in = correction jyg 20220124 -! print *,'XXXX8 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i) - sigmaw(i) = sigmaw_targ - IF (iflag_wk_pop_dyn >= 1) THEN -!! awdens(i) = 0. -!! wdens(i) = 0. - wdens_targ = 0. - d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i) -!! d_wdens2(i) = wdens_targ - wdens(i) - d_wdens2(i) = wdens_targ - wdens_in(i) ! jyg 20220916 - wdens(i) = wdens_targ - wdens_targ = 0. -!!jyg: bug fix : the d_adens_bnd2 computation must be before the update of awdens. - IF (iflag_wk_pop_dyn == 2) THEN - d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i) - ENDIF ! (iflag_wk_pop_dyn == 2) -!! d_awdens2(i) = wdens_targ - awdens(i) - d_awdens2(i) = wdens_targ - awdens_in(i) ! jyg 20220916 - awdens(i) = wdens_targ -!! IF (iflag_wk_pop_dyn == 2) THEN -!! d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i) -!! ENDIF ! (iflag_wk_pop_dyn == 2) - ENDIF ! (iflag_wk_pop_dyn >= 1) - ELSE ! (kill_wake(i)) - wape(i) = wape2(i) - cstar(i) = cstar2(i) - END IF ! (kill_wake(i)) - ! c print*,'wape wape2 ktopw OK_qx_qw =', - ! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i) - END DO - - IF (prt_level>=10) THEN - PRINT *, 'wake-6, wape wape2 ktopw OK_qx_qw =', & - wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout) - ENDIF - - - ! ----------------------------------------------------------------- - ! 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 -!jyg< -!! IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN - IF (ok_qx_qw(i)) THEN -!>jyg - ! 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 -!jyg< - IF (iflag_wk_pop_dyn >= 1) THEN - DO i = 1, klon - IF (ok_qx_qw(i)) THEN - d_sig_gen2(i) = d_sig_gen2(i)/dtime - d_sig_death2(i) = d_sig_death2(i)/dtime - d_sig_col2(i) = d_sig_col2(i)/dtime - d_sig_spread2(i) = d_sig_spread2(i)/dtime - d_sig_bnd2(i) = d_sig_bnd2(i)/dtime - d_sigmaw2(i) = d_sigmaw2(i)/dtime -! print *,'XXXX9 d_sigmaw2(i), sigmaw(i), dtime ', d_sigmaw2(i), sigmaw(i), dtime -! - d_dens_gen2(i) = d_dens_gen2(i)/dtime - d_dens_death2(i) = d_dens_death2(i)/dtime - d_dens_col2(i) = d_dens_col2(i)/dtime - d_dens_bnd2(i) = d_dens_bnd2(i)/dtime - d_awdens2(i) = d_awdens2(i)/dtime - d_wdens2(i) = d_wdens2(i)/dtime - ENDIF - ENDDO - IF (iflag_wk_pop_dyn == 2) THEN - DO i = 1, klon - IF (ok_qx_qw(i)) THEN - d_adens_death2(i) = d_adens_death2(i)/dtime - d_adens_icol2(i) = d_adens_icol2(i)/dtime - d_adens_acol2(i) = d_adens_acol2(i)/dtime - d_adens_bnd2(i) = d_adens_bnd2(i)/dtime - ENDIF - ENDDO - ENDIF ! (iflag_wk_pop_dyn == 2) - ENDIF ! (iflag_wk_pop_dyn >= 1) - -!>jyg - - RETURN -END SUBROUTINE wake - -SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon_loc, 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_loc - 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_loc - 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 positivite 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 - - - -SUBROUTINE pkupper (klon, klev, ptop, ph, pupper, kupper) - -USE wake_ini_mod , ONLY : wk_pupper -IMPLICIT NONE - -INTEGER, INTENT(IN) :: klon,klev -REAL, INTENT(IN), DIMENSION (klon,klev+1) :: ph -REAL, INTENT(IN), DIMENSION (klon) :: ptop -REAL, INTENT(OUT), DIMENSION (klon) :: pupper -INTEGER, INTENT(OUT), DIMENSION (klon) :: kupper -INTEGER :: i,k - - - kupper = 0 - -IF (wk_pupper<1.) THEN - ! 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) = wk_pupper*ph(i, 1) - pupper(i) = max(pupper(i), 45000.) - ! cc Pupper(i) = 60000. - END DO - -ELSE - - DO i=1, klon - ! pupper(i) = wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) - pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-5000.) - END DO -END IF - - ! -5/ Determination de kupper - - DO k = klev, 1, -1 - DO i = 1, klon - IF (ph(i,k+1)jyg+mlt -! -!jyg< - d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i)) -!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) - d_dens_bnd(i) = d_wdens_targ - d_wdens(i) - d_wdens(i) = d_wdens_targ -!! d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i)) -!>jyg - -!jyg+mlt< -!! d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) & -!! + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) & -!! - sigmaw(i)*tau_wk_inv_min )*dtimesub - d_sig_gen(i) = wgen(i)*aa0 -!! print*, 'XXX sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min', & -!! sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min - d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min -!! - - d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos - d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos - d_sig_spread(i) = gfl(i)*cstar(i) - d_sig_gen(i) = d_sig_gen(i)*dtimesub - d_sig_death(i) = d_sig_death(i)*dtimesub - d_sig_col(i) = d_sig_col(i)*dtimesub - d_sig_spread(i) = d_sig_spread(i)*dtimesub - d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i) -!>jyg+mlt -! -!jyg< - d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i)) -!! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i) -!! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i) - d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i) - d_sigmaw(i) = d_sigmaw_targ -!! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i)) -!>jyg - ENDIF - ENDDO - - IF (prt_level >= 10) THEN - print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', & - cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) - print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', & - wdens(1), awdens(1), act(1), d_awdens(1) - print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', & - wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1) - print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', & - d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) - ENDIF - - RETURN - END SUBROUTINE wake_popdyn_1 - - SUBROUTINE wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, & - sigmaw, wdens, awdens, & !! states variables - gfl, cstar, cin, wape, rad_wk, & - d_sigmaw, d_wdens, d_awdens, & !! tendences - cont_fact, & - d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, & - d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, & - d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd ) - - - - USE wake_ini_mod , ONLY : wake_ini - USE wake_ini_mod , ONLY : prt_level,RG - USE wake_ini_mod , ONLY : stark, wdens_ref - USE wake_ini_mod , ONLY : tau_cv, rzero, aa0 - USE wake_ini_mod , ONLY : iflag_wk_pop_dyn, wdensmin - USE wake_ini_mod , ONLY : sigmad, cstart, sigmaw_max - -IMPLICIT NONE - - INTEGER, INTENT(IN) :: klon,klev - LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv - REAL, INTENT(IN) :: dtimesub - REAL, DIMENSION (klon), INTENT(IN) :: wgen !! B = birth rate of wakes - REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw !! sigma = fractional area of wakes - REAL, DIMENSION (klon), INTENT(INOUT) :: wdens !! D = number of wakes per unit area - REAL, DIMENSION (klon), INTENT(INOUT) :: awdens !! A = number of active wakes per unit area - REAL, DIMENSION (klon), INTENT(IN) :: gfl !! Lg = gust front lenght per unit area - REAL, DIMENSION (klon), INTENT(IN) :: cstar !! C* = spreading velocity of wakes - REAL, DIMENSION (klon), INTENT(IN) :: cin, wape ! RM : A Faire disparaitre - REAL, DIMENSION (klon), INTENT(IN) :: rad_wk !! r = wake radius - -! - REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_wdens, d_awdens - REAL, DIMENSION (klon), INTENT(OUT) :: cont_fact !! RM facteur de contact = 2 pi * rad * C* - ! Some components of the tendencies of state variables - REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd - REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd - REAL, DIMENSION (klon), INTENT(OUT) :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd - - -!! internal variables - - INTEGER :: i, k - REAL, DIMENSION (klon) :: tau_wk_inv !! tau = life time of wakes - REAL :: tau_wk_inv_min - REAL, DIMENSION (klon) :: tau_prime !! tau_prime = life time of actives wakes - REAL :: d_wdens_targ, d_sigmaw_targ - - -!! Equations -!! dD/dt = B - (D-A)/tau - f D^2 -!! dA/dt = B - A/tau_prime + f (D-A)^2 - f A^2 -!! dsigma/dt = B a0 - sigma/D (D-A)/tau + Lg C* - f (D-A)^2 (sigma/D-a0) -!! -!! f = 2 (B (a0-sigma/D) + Lg C*) / (2 (D-A)^2 (2 sigma/D-a0) + D (1-2 sigma)) - - - DO i = 1, klon - ! print*, 'XXX wk_adv(i)', wk_adv(i) - IF (wk_adv(i)) THEN -!! tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.) - tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.) - tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub) - tau_prime(i) = tau_cv -!! cont_fact(i) = 2.*(wgen(i)*(aa0-sigmaw(i)/wdens(i)) + gfl(i)*cstar(i)) / & -!! (2.*(wdens(i)-awdens(i))**2*(2.*sigmaw(i)/wdens(i) - aa0) + wdens(i)*(1.-2.*sigmaw(i))) -!! cont_fact(i) = 2.*3.14*rad_wk(i)*cstar(i) ! bug -!! cont_fact(i) = 4.*3.14*rad_wk(i)*cstar(i) - cont_fact(i) = 2.*gfl(i)*cstar(i)/wdens(i) - - d_sig_gen(i) = wgen(i)*aa0 - d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min - d_sig_col(i) = - cont_fact(i)*(wdens(i)-awdens(i))**2*(2.*sigmaw(i)/wdens(i)-aa0) - d_sig_spread(i) = gfl(i)*cstar(i) -! - d_sig_gen(i) = d_sig_gen(i)*dtimesub - d_sig_death(i) = d_sig_death(i)*dtimesub - d_sig_col(i) = d_sig_col(i)*dtimesub - d_sig_spread(i) = d_sig_spread(i)*dtimesub - d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i) - - - d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i)) -!! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i) -!! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i) - d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i) - d_sigmaw(i) = d_sigmaw_targ -!! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i)) - - - d_dens_gen(i) = wgen(i) - d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min - d_dens_col(i) = - cont_fact(i)*wdens(i)**2 -! - d_dens_gen(i) = d_dens_gen(i)*dtimesub - d_dens_death(i) = d_dens_death(i)*dtimesub - d_dens_col(i) = d_dens_col(i)*dtimesub - d_wdens(i) = d_dens_gen(i) + d_dens_death(i) + d_dens_col(i) - - - d_adens_death(i) = -awdens(i)/tau_prime(i) - d_adens_icol(i) = cont_fact(i)*(wdens(i)-awdens(i))**2 - d_adens_acol(i) = - cont_fact(i)*awdens(i)**2 -! - d_adens_death(i) = d_adens_death(i)*dtimesub - d_adens_icol(i) = d_adens_icol(i)*dtimesub - d_adens_acol(i) = d_adens_acol(i)*dtimesub - d_awdens(i) = d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i) - -!! - d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i)) -!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) - d_dens_bnd(i) = d_wdens_targ - d_wdens(i) - d_wdens(i) = d_wdens_targ - - d_wdens_targ = min(max(d_awdens(i),-awdens(i)), wdens(i)-awdens(i)) -!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i) - d_adens_bnd(i) = d_wdens_targ - d_awdens(i) - d_awdens(i) = d_wdens_targ - - - - ENDIF - ENDDO - - IF (prt_level >= 10) THEN - print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1) ', & - cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1) - print *,'wake, wdens(1), awdens(1), d_awdens(1) ', & - wdens(1), awdens(1), d_awdens(1) - print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', & - d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) - ENDIF -sigmaw=sigmaw+d_sigmaw -wdens=wdens+d_wdens -awdens=awdens+d_awdens - - RETURN - END SUBROUTINE wake_popdyn_2 - Index: LMDZ6/trunk/libf/phylmd/wake_ini_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/wake_ini_mod.F90 (revision 4586) +++ (revision ) @@ -1,219 +1,0 @@ -MODULE wake_ini_mod -IMPLICIT NONE - - - ! ============================================================================ - - - ! But : Decrire le comportement des poches froides apparaissant dans les - ! grands systemes convectifs, et fournir l'energie disponible pour - ! le declenchement de nouvelles colonnes convectives. - - ! State variables : - ! deltatw : temperature difference between wake and off-wake regions - ! deltaqw : specific humidity difference between wake and off-wake regions - ! sigmaw : fractional area covered by wakes. - ! wdens : number of wakes per unit area - - ! ------------------------------------------------------------------------- - ! Declaration de variables - ! ------------------------------------------------------------------------- - - ! Variables a fixer -!jyg< -!! REAL, SAVE :: stark, wdens_ref, coefgw, alpk - INTEGER, SAVE, PROTECTED :: prt_level - REAL, SAVE, PROTECTED, DIMENSION(2) :: wdens_ref - REAL, SAVE, PROTECTED :: stark, coefgw, alpk, wk_pupper -!>jyg - REAL, SAVE, PROTECTED :: crep_upper, crep_sol - !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, wk_pupper, crep_upper, crep_sol) - - REAL, SAVE, PROTECTED :: tau_cv - !$OMP THREADPRIVATE(tau_cv) - REAL, SAVE, PROTECTED :: rzero, aa0 ! minimal wake radius and area - !$OMP THREADPRIVATE(rzero, aa0) - - LOGICAL, SAVE, PROTECTED :: flag_wk_check_trgl - !$OMP THREADPRIVATE(flag_wk_check_trgl) - INTEGER, SAVE, PROTECTED :: iflag_wk_act - !$OMP THREADPRIVATE(iflag_wk_act) - - INTEGER, SAVE, PROTECTED :: iflag_wk_check_trgl - !$OMP THREADPRIVATE(iflag_wk_check_trgl) - INTEGER, SAVE, PROTECTED :: iflag_wk_pop_dyn - !$OMP THREADPRIVATE(iflag_wk_pop_dyn) - - INTEGER, SAVE, PROTECTED :: iflag_wk_profile - !$OMP THREADPRIVATE(iflag_wk_profile) - - REAL, SAVE, PROTECTED :: wdensmin - !$OMP THREADPRIVATE(wdensmin) - REAL, SAVE, PROTECTED :: sigmad, hwmin, wapecut, cstart - !$OMP THREADPRIVATE(sigmad, hwmin, wapecut, cstart) - REAL, SAVE, PROTECTED :: sigmaw_max - !$OMP THREADPRIVATE(sigmaw_max) - REAL, SAVE, PROTECTED :: dens_rate - !$OMP THREADPRIVATE(dens_rate) - REAL, SAVE, PROTECTED :: epsilon_loc - !$OMP THREADPRIVATE(epsilon_loc) - REAL, SAVE, PROTECTED :: epsim1,RG,RD - !$OMP THREADPRIVATE(epsim1,RG,RD) - - - -CONTAINS - - ! ========================================================================= - SUBROUTINE wake_ini(rg_in,rd_in,rv_in,prt_lev) - ! ========================================================================= - - ! ************************************************************** - ! * - ! WAKE * - ! retour a un Pupper fixe * - ! * - ! written by : GRANDPEIX Jean-Yves 09/03/2000 * - ! modified by : ROEHRIG Romain 01/29/2007 * - ! ************************************************************** - - ! ------------------------------------------------------------------------- - ! Initialisations - ! ------------------------------------------------------------------------- - - USE ioipsl_getin_p_mod, ONLY : getin_p - real eps - integer, intent(in) :: prt_lev - real, intent(in) :: rg_in,rd_in,rv_in - - prt_level=prt_lev - epsilon_loc=1.E-15 - wapecut=1. ! previously 5. - ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10. - sigmad=0.02 - hwmin=10. -!! DATA wdensmin/1.e-12/ - wdensmin=1.e-14 - ! cc nrlmd - sigmaw_max=0.4 - dens_rate=0.1 - - eps = rd_in/rv_in - epsim1 = 1.0/eps - 1.0 - RG=rg_in - RD=rd_in - - - ! 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 gravite - ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) - ! wdens : Densite surfacique de poche froide - ! ------------------------------------------------------------------------- - - ! 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 - - - - crep_upper = 0.9 - crep_sol = 1.0 - - ! Get wapecut from parameter file - wapecut = 1. - -print*,'wapecut',wapecut - CALL getin_p('wapecut', wapecut) -print*,'wapecut',wapecut - - ! cc nrlmd Lecture du fichier wake_param.data - - - ! cc nrlmd Lecture du fichier wake_param.data - stark=0.33 - CALL getin_p('stark',stark) - cstart = stark*sqrt(2.*wapecut) - - alpk=0.25 - CALL getin_p('alpk',alpk) - - wk_pupper=0.6 - CALL getin_p('wk_pupper',wk_pupper) - - -!jyg< -!! wdens_ref=8.E-12 -!! CALL getin_p('wdens_ref',wdens_ref) - wdens_ref(1)=8.E-12 - wdens_ref(2)=8.E-12 - CALL getin_p('wdens_ref_o',wdens_ref(1)) !wake number per unit area ; ocean - CALL getin_p('wdens_ref_l',wdens_ref(2)) !wake number per unit area ; land -!>jyg -! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!!!!!!!! Population dynamics parameters !!!!!!!!!!!!!!!!!!!!!!!!!!!! -!------------------------------------------------------------------------ - - iflag_wk_pop_dyn = 0 - CALL getin_p('iflag_wk_pop_dyn',iflag_wk_pop_dyn) ! switch between wdens prescribed - ! and wdens prognostic - iflag_wk_act = 0 - CALL getin_p('iflag_wk_act',iflag_wk_act) ! 0: act(:)=0. - ! 1: act(:)=1. - ! 2: act(:)=f(Wape) - - iflag_wk_profile = 0 - CALL getin_p('iflag_wk_profile',iflag_wk_profile) ! switch between wdens prescribed - ! and wdens prognostic - rzero = 5000. - CALL getin_p('rzero_wk', rzero) - aa0 = 3.14*rzero*rzero -! - tau_cv = 4000. - CALL getin_p('tau_cv', tau_cv) - -!------------------------------------------------------------------------ - - coefgw=4. - CALL getin_p('coefgw',coefgw) - - WRITE(*,*) 'stark=', stark - WRITE(*,*) 'alpk=', alpk - WRITE(*,*) 'wk_pupper=', wk_pupper -!jyg< -!! WRITE(*,*) 'wdens_ref=', wdens_ref - WRITE(*,*) 'wdens_ref_o=', wdens_ref(1) - WRITE(*,*) 'wdens_ref_l=', wdens_ref(2) -!>jyg - WRITE(*,*) 'iflag_wk_pop_dyn=',iflag_wk_pop_dyn - WRITE(*,*) 'iflag_wk_act',iflag_wk_act - WRITE(*,*) 'coefgw=', coefgw - - flag_wk_check_trgl=.false. - CALL getin_p('flag_wk_check_trgl ', flag_wk_check_trgl) - WRITE(*,*) 'flag_wk_check_trgl=', flag_wk_check_trgl - WRITE(*,*) 'flag_wk_check_trgl OBSOLETE. Utilisr iflag_wk_check_trgl plutot' - iflag_wk_check_trgl=0 ; IF (flag_wk_check_trgl) iflag_wk_check_trgl=1 - CALL getin_p('iflag_wk_check_trgl ', iflag_wk_check_trgl) - WRITE(*,*) 'iflag_wk_check_trgl=', iflag_wk_check_trgl - - RETURN - -END SUBROUTINE wake_ini - -END MODULE wake_ini_mod