! $Id: wake.F90 1999 2014-03-20 09:57:19Z fairhead $ SUBROUTINE wake(p, ph, pi, dtime, sigd_con, te0, qe0, omgb, dtdwn, dqdwn, & amdwn, amup, dta, dqa, wdtpbl, wdqpbl, udtpbl, udqpbl, deltatw, deltaqw, & dth, hw, sigmaw, wape, fip, gfl, dtls, dqls, ktopw, omgbdth, dp_omgb, & wdens, tu, qu, dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, spread, cstar, & d_deltat_gw, d_deltatw2, d_deltaqw2) ! ************************************************************** ! * ! WAKE * ! retour a un Pupper fixe * ! * ! written by : GRANDPEIX Jean-Yves 09/03/2000 * ! modified by : ROEHRIG Romain 01/29/2007 * ! ************************************************************** USE dimphy 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. ! Variables d'etat : deltatw : ecart de temperature wake-undisturbed ! area ! deltaqw : ecart d'humidite wake-undisturbed area ! sigmaw : fraction d'aire occupee par la poche. ! 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 : hauteur de la poche ! dp_omgb : vertical gradient of large scale omega ! 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) ! spread : spreading term in dt_wake and dq_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 ! te0 : 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) ! 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) ! 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 ! te : 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 : hw initial ! sigmaw0: sigmaw initial ! 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 période de onde de gravité ! Cgw : vitesse de propagation de onde de gravité ! LL : distance entre 2 poches ! ------------------------------------------------------------------------- ! Déclaration de variables ! ------------------------------------------------------------------------- include "dimensions.h" include "YOMCST.h" include "cvthermo.h" include "iniprint.h" ! Arguments en entree ! -------------------- REAL, DIMENSION (klon, klev) :: p, pi REAL, DIMENSION (klon, klev+1) :: ph, omgb REAL dtime REAL, DIMENSION (klon, klev) :: te0, qe0 REAL, DIMENSION (klon, klev) :: dtdwn, dqdwn REAL, DIMENSION (klon, klev) :: wdtpbl, wdqpbl REAL, DIMENSION (klon, klev) :: udtpbl, udqpbl REAL, DIMENSION (klon, klev) :: amdwn, amup REAL, DIMENSION (klon, klev) :: dta, dqa REAL, DIMENSION (klon) :: sigd_con ! Sorties ! -------- REAL, DIMENSION (klon, klev) :: deltatw, deltaqw, dth REAL, DIMENSION (klon, klev) :: tu, qu REAL, DIMENSION (klon, klev) :: dtls, dqls REAL, DIMENSION (klon, klev) :: dtke, dqke REAL, DIMENSION (klon, klev) :: dtpbl, dqpbl REAL, DIMENSION (klon, klev) :: spread REAL, DIMENSION (klon, klev) :: d_deltatgw REAL, DIMENSION (klon, klev) :: d_deltatw2, d_deltaqw2 REAL, DIMENSION (klon, klev+1) :: omgbdth, omg REAL, DIMENSION (klon, klev) :: dp_omgb, dp_deltomg REAL, DIMENSION (klon, klev) :: d_deltat_gw REAL, DIMENSION (klon) :: hw, sigmaw, wape, fip, gfl, cstar REAL, DIMENSION (klon) :: wdens INTEGER, DIMENSION (klon) :: ktopw ! Variables internes ! ------------------- ! Variables à fixer REAL alon REAL coefgw REAL :: wdens_ref REAL stark REAL alpk REAL delta_t_min INTEGER nsub REAL dtimesub REAL sigmad, hwmin, wapecut REAL :: sigmaw_max REAL :: dens_rate REAL wdens0 ! IM 080208 LOGICAL, DIMENSION (klon) :: gwake ! Variables de sauvegarde REAL, DIMENSION (klon, klev) :: deltatw0 REAL, DIMENSION (klon, klev) :: deltaqw0 REAL, DIMENSION (klon, klev) :: te, qe REAL, DIMENSION (klon) :: sigmaw0, sigmaw1 ! Variables pour les GW REAL, DIMENSION (klon) :: ll REAL, DIMENSION (klon, klev) :: n2 REAL, DIMENSION (klon, klev) :: cgw REAL, DIMENSION (klon, klev) :: tgw ! Variables liées au calcul de hw REAL, DIMENSION (klon) :: ptop_provis, ptop, ptop_new REAL, DIMENSION (klon) :: sum_dth REAL, DIMENSION (klon) :: dthmin REAL, DIMENSION (klon) :: z, dz, hw0 INTEGER, DIMENSION (klon) :: ktop, kupper ! Sub-timestep tendencies and related variables REAL d_deltatw(klon, klev), d_deltaqw(klon, klev) REAL d_te(klon, klev), d_qe(klon, klev) REAL d_sigmaw(klon), alpha(klon) REAL q0_min(klon), q1_min(klon) LOGICAL wk_adv(klon), ok_qx_qw(klon) REAL epsilon DATA epsilon/1.E-15/ ! Autres variables internes INTEGER isubstep, k, i REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu REAL, DIMENSION (klon) :: sum_dq, sum_rho REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn REAL, DIMENSION (klon) :: av_thu, av_tu, av_qu, av_thvu REAL, DIMENSION (klon) :: av_dth, av_dq, av_rho REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn REAL, DIMENSION (klon, klev) :: rho, rhow REAL, DIMENSION (klon, klev+1) :: rhoh REAL, DIMENSION (klon, klev) :: rhow_moyen REAL, DIMENSION (klon, klev) :: zh REAL, DIMENSION (klon, klev+1) :: zhh REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2 REAL, DIMENSION (klon, klev) :: the, thu ! REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw REAL, DIMENSION (klon, klev+1) :: omgbw REAL, DIMENSION (klon) :: pupper REAL, DIMENSION (klon) :: omgtop REAL, DIMENSION (klon, klev) :: dp_omgbw REAL, DIMENSION (klon) :: ztop, dztop REAL, DIMENSION (klon, klev) :: alpha_up REAL, DIMENSION (klon) :: rre1, rre2 REAL :: rrd1, rrd2 REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2 REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq REAL, DIMENSION (klon, klev) :: omgbdq REAL, DIMENSION (klon) :: ff, gg REAL, DIMENSION (klon) :: wape2, cstar2, heff REAL, DIMENSION (klon, klev) :: crep REAL crep_upper, crep_sol REAL, DIMENSION (klon, klev) :: ppi ! cc nrlmd REAL, DIMENSION (klon) :: death_rate, nat_rate REAL, DIMENSION (klon, klev) :: entr REAL, DIMENSION (klon, klev) :: detr ! ------------------------------------------------------------------------- ! Initialisations ! ------------------------------------------------------------------------- ! print*, 'wake initialisations' ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10. ! ------------------------------------------------------------------------- DATA wapecut, sigmad, hwmin/5., .02, 10./ ! cc nrlmd DATA sigmaw_max/0.4/ DATA dens_rate/0.1/ ! cc ! Longueur de maille (en m) ! ------------------------------------------------------------------------- ! ALON = 3.e5 alon = 1.E6 ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei) ! coefgw : Coefficient pour les ondes de gravité ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) ! wdens : Densité de poche froide par maille ! ------------------------------------------------------------------------- ! cc nrlmd coefgw=10 ! coefgw=1 ! wdens0 = 1.0/(alon**2) ! cc nrlmd wdens = 1.0/(alon**2) ! cc nrlmd stark = 0.50 ! CRtest ! cc nrlmd alpk=0.1 ! alpk = 1.0 ! alpk = 0.5 ! alpk = 0.05 stark = 0.33 alpk = 0.25 wdens_ref = 8.E-12 coefgw = 4. crep_upper = 0.9 crep_sol = 1.0 ! cc nrlmd Lecture du fichier wake_param.data OPEN (99, FILE='wake_param.data', STATUS='old', FORM='formatted', ERR=9999) READ (99, *, END=9998) stark READ (99, *, END=9998) alpk READ (99, *, END=9998) wdens_ref READ (99, *, END=9998) coefgw 9998 CONTINUE CLOSE (99) 9999 CONTINUE ! Initialisation de toutes des densites a wdens_ref. ! Les densites peuvent evoluer si les poches debordent ! (voir au tout debut de la boucle sur les substeps) wdens = wdens_ref ! print*,'stark',stark ! print*,'alpk',alpk ! print*,'wdens',wdens ! print*,'coefgw',coefgw ! cc ! Minimum value for |T_wake - T_undist|. Used for wake top definition ! ------------------------------------------------------------------------- delta_t_min = 0.2 ! 1. - Save initial values and initialize tendencies ! -------------------------------------------------- DO k = 1, klev DO i = 1, klon ppi(i, k) = pi(i, k) deltatw0(i, k) = deltatw(i, k) deltaqw0(i, k) = deltaqw(i, k) te(i, k) = te0(i, k) qe(i, k) = qe0(i, k) dtls(i, k) = 0. dqls(i, k) = 0. d_deltat_gw(i, k) = 0. d_te(i, k) = 0. d_qe(i, k) = 0. d_deltatw(i, k) = 0. d_deltaqw(i, k) = 0. ! IM 060508 beg d_deltatw2(i, k) = 0. d_deltaqw2(i, k) = 0. ! IM 060508 end END DO END DO ! sigmaw1=sigmaw ! IF (sigd_con.GT.sigmaw1) THEN ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con ! ENDIF DO i = 1, klon ! c sigmaw(i) = amax1(sigmaw(i),sigd_con(i)) sigmaw(i) = amax1(sigmaw(i), sigmad) sigmaw(i) = amin1(sigmaw(i), 0.99) sigmaw0(i) = sigmaw(i) wape(i) = 0. wape2(i) = 0. d_sigmaw(i) = 0. ktopw(i) = 0 END DO ! 2. - Prognostic part ! -------------------- ! 2.1 - Undisturbed area and Wake integrals ! --------------------------------------------------------- DO i = 1, klon z(i) = 0. ktop(i) = 0 kupper(i) = 0 sum_thu(i) = 0. sum_tu(i) = 0. sum_qu(i) = 0. sum_thvu(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_rho(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thu(i) = 0. av_tu(i) = 0. av_qu(i) = 0. av_thvu(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_rho(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END DO ! Distance between wakes DO i = 1, klon ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i)) END DO ! Potential temperatures and humidity ! ---------------------------------------------------------- DO k = 1, klev DO i = 1, klon ! write(*,*)'wake 1',i,k,rd,te(i,k) rho(i, k) = p(i, k)/(rd*te(i,k)) ! write(*,*)'wake 2',rho(i,k) IF (k==1) THEN ! write(*,*)'wake 3',i,k,rd,te(i,k) rhoh(i, k) = ph(i, k)/(rd*te(i,k)) ! write(*,*)'wake 4',i,k,rd,te(i,k) zhh(i, k) = 0 ELSE ! write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1)) rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1))) ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1) zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1) END IF ! write(*,*)'wake 7',ppi(i,k) the(i, k) = te(i, k)/ppi(i, k) thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i) qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k))) rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k))) dth(i, k) = deltatw(i, k)/ppi(i, k) END DO END DO DO k = 1, klev - 1 DO i = 1, klon IF (k==1) THEN n2(i, k) = 0 ELSE n2(i, k) = amax1(0., -rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i, & k-1))/(p(i,k+1)-p(i,k-1))) END IF zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2 cgw(i, k) = sqrt(n2(i,k))*zh(i, k) tgw(i, k) = coefgw*cgw(i, k)/ll(i) END DO END DO DO i = 1, klon n2(i, klev) = 0 zh(i, klev) = 0 cgw(i, klev) = 0 tgw(i, klev) = 0 END DO ! Calcul de la masse volumique moyenne de la colonne (bdlmd) ! ----------------------------------------------------------------- DO k = 1, klev DO i = 1, klon epaisseur1(i, k) = 0. epaisseur2(i, k) = 0. END DO END DO DO i = 1, klon epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1. epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1. rhow_moyen(i, 1) = rhow(i, 1) END DO DO k = 2, klev DO i = 1, klon epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1. epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k) rhow_moyen(i, k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+rhow(i,k)* & epaisseur1(i,k))/epaisseur2(i, k) END DO END DO ! Choose an integration bound well above wake top ! ----------------------------------------------------------------- ! Pupper = 50000. ! melting level ! Pupper = 60000. ! Pupper = 80000. ! essais pour case_e DO i = 1, klon pupper(i) = 0.6*ph(i, 1) pupper(i) = max(pupper(i), 45000.) ! cc Pupper(i) = 60000. END DO ! Determine Wake top pressure (Ptop) from buoyancy integral ! -------------------------------------------------------- ! -1/ Pressure of the level where dth becomes less than delta_t_min. DO i = 1, klon ptop_provis(i) = ph(i, 1) END DO DO k = 2, klev DO i = 1, klon ! IM v3JYG; ptop_provis(i).LT. ph(i,1) IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. & ptop_provis(i)==ph(i,1)) THEN ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO ! -2/ dth integral DO i = 1, klon sum_dth(i) = 0. dthmin(i) = -delta_t_min z(i) = 0. END DO DO k = 1, klev DO i = 1, klon dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) dthmin(i) = amin1(dthmin(i), dth(i,k)) END IF END DO END DO ! -3/ height of triangle with area= sum_dth and base = dthmin DO i = 1, klon hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) hw0(i) = amax1(hwmin, hw0(i)) END DO ! -4/ now, get Ptop DO i = 1, klon z(i) = 0. ptop(i) = ph(i, 1) END DO DO k = 1, klev DO i = 1, klon dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i)) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i) END IF END DO END DO ! -5/ Determination de ktop et kupper DO k = klev, 1, -1 DO i = 1, klon IF (ph(i,k+1)-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)=kupper(i)) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. END IF END DO END DO ! Vertical gradient of LS omega DO k = 1, klev DO i = 1, klon IF (k<=kupper(i)) THEN dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k)) END IF END DO END DO ! Integrals (and wake top level number) ! -------------------------------------- ! Initialize sum_thvu to 1st level virt. pot. temp. DO i = 1, klon z(i) = 1. dz(i) = 1. sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i) sum_dth(i) = 0. END DO DO k = 1, klev DO i = 1, klon dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) END IF END DO END DO DO i = 1, klon hw0(i) = z(i) END DO ! 2.1 - WAPE and mean forcing computation ! --------------------------------------- ! --------------------------------------- ! Means DO i = 1, klon av_thu(i) = sum_thu(i)/hw0(i) av_tu(i) = sum_tu(i)/hw0(i) av_qu(i) = sum_qu(i)/hw0(i) av_thvu(i) = sum_thvu(i)/hw0(i) ! av_thve = sum_thve/hw0 av_dth(i) = sum_dth(i)/hw0(i) av_dq(i) = sum_dq(i)/hw0(i) av_rho(i) = sum_rho(i)/hw0(i) av_dtdwn(i) = sum_dtdwn(i)/hw0(i) av_dqdwn(i) = sum_dqdwn(i)/hw0(i) wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i & )+av_dth(i)*av_dq(i)))/av_thvu(i) END DO ! 2.2 Prognostic variable update ! ------------------------------ ! Filter out bad wakes DO k = 1, klev DO i = 1, klon IF (wape(i)<0.) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. END IF END DO END DO DO i = 1, klon IF (wape(i)<0.) THEN wape(i) = 0. cstar(i) = 0. hw(i) = hwmin sigmaw(i) = amax1(sigmad, sigd_con(i)) fip(i) = 0. gwake(i) = .FALSE. ELSE cstar(i) = stark*sqrt(2.*wape(i)) gwake(i) = .TRUE. END IF END DO ! Check qx and qw positivity ! -------------------------- DO i = 1, klon q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), (qe(i, & 1)+(1.-sigmaw(i))*deltaqw(i,1))) END DO DO k = 2, klev DO i = 1, klon q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), (qe(i, & k)+(1.-sigmaw(i))*deltaqw(i,k))) IF (q1_min(i)<=q0_min(i)) THEN q0_min(i) = q1_min(i) END IF END DO END DO DO i = 1, klon ok_qx_qw(i) = q0_min(i) >= 0. alpha(i) = 1. END DO ! C ----------------------------------------------------------------- ! Sub-time-stepping ! ----------------- nsub = 10 dtimesub = dtime/nsub ! ------------------------------------------------------------ DO isubstep = 1, nsub ! ------------------------------------------------------------ ! wk_adv is the logical flag enabling wake evolution in the time advance ! loop DO i = 1, klon wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1. END DO ! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement ! négatif de ktop à kupper -------- ! cc On calcule pour cela une densité wdens0 pour laquelle on ! aurait un entrainement nul --- DO i = 1, klon ! c print *,' isubstep,wk_adv(i),cstar(i),wape(i) ', ! c $ isubstep,wk_adv(i),cstar(i),wape(i) IF (wk_adv(i) .AND. cstar(i)>0.01) THEN omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + & rg*amup(i, kupper(i)+1)/(1.-sigmaw(i)) wdens0 = (sigmaw(i)/(4.*3.14))*((1.-sigmaw(i))*omg(i,kupper(i)+1)/(( & ph(i,1)-pupper(i))*cstar(i)))**(2) IF (wdens(i)<=wdens0*1.1) THEN wdens(i) = wdens0 END IF ! c print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i) ! c $ ,ph(i,1)-pupper(i)', ! c $ omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i) ! c $ ,ph(i,1)-pupper(i) END IF END DO ! cc nrlmd DO i = 1, klon IF (wk_adv(i)) THEN gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i)) sigmaw(i) = amin1(sigmaw(i), sigmaw_max) END IF END DO DO i = 1, klon IF (wk_adv(i)) THEN ! cc nrlmd Introduction du taux de mortalité des poches et ! test sur sigmaw_max=0.4 ! cc d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub IF (sigmaw(i)>=sigmaw_max) THEN death_rate(i) = gfl(i)*cstar(i)/sigmaw(i) ELSE death_rate(i) = 0. END IF d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* & dtimesub ! $ - nat_rate(i)*sigmaw(i)*dtimesub ! c print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), ! c $ death_rate(i),ktop(i),kupper(i)', ! c $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), ! c $ death_rate(i),ktop(i),kupper(i) ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub ! sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!! ! wdens = wdens0/(10.*sigmaw) ! sigmaw =max(sigmaw,sigd_con) ! sigmaw =max(sigmaw,sigmad) END IF END DO ! calcul de la difference de vitesse verticale poche - zone non perturbee ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on ! definit ! IM 060208 au niveau k=1..? DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd dp_deltomg(i, k) = 0. END IF END DO END DO DO k = 1, klev + 1 DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd omg(i, k) = 0. END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN z(i) = 0. omg(i, 1) = 0. dp_deltomg(i, 1) = -(gfl(i)*cstar(i))/(sigmaw(i)*(1-sigmaw(i))) END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=ktop(i)) THEN dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg) z(i) = z(i) + dz(i) dp_deltomg(i, k) = dp_deltomg(i, 1) omg(i, k) = dp_deltomg(i, 1)*z(i) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg) ztop(i) = z(i) + dztop(i) omgtop(i) = dp_deltomg(i, 1)*ztop(i) END IF END DO ! ----------------- ! From m/s to Pa/s ! ----------------- DO i = 1, klon IF (wk_adv(i)) THEN omgtop(i) = -rho(i, ktop(i))*rg*omgtop(i) dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1)) END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=ktop(i)) THEN omg(i, k) = -rho(i, k)*rg*omg(i, k) dp_deltomg(i, k) = dp_deltomg(i, 1) END IF END DO END DO ! raccordement lineaire de omg de ptop a pupper DO i = 1, klon IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + & rg*amup(i, kupper(i)+1)/(1.-sigmaw(i)) dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ & (ptop(i)-pupper(i)) END IF END DO ! c DO i=1,klon ! c print*,'Pente entre 0 et kupper (référence)' ! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1)) ! c print*,'Pente entre ktop et kupper' ! c $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i)) ! c ENDDO ! c DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN dp_deltomg(i, k) = dp_deltomg(i, kupper(i)) omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i)) END IF END DO END DO ! cc nrlmd ! c DO i=1,klon ! c print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1) ! c END DO ! cc ! -- Compute wake average vertical velocity omgbw DO k = 1, klev + 1 DO i = 1, klon IF (wk_adv(i)) THEN omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k) END IF END DO END DO ! -- and its vertical gradient dp_omgbw DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k)) END IF END DO END DO ! -- Upstream coefficients for omgb velocity ! -- (alpha_up(k) is the coefficient of the value at level k) ! -- (1-alpha_up(k) is the coefficient of the value at level k-1) DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN alpha_up(i, k) = 0. IF (omgb(i,k)>0.) alpha_up(i, k) = 1. END IF END DO END DO ! Matrix expressing [The,deltatw] from [Th1,Th2] DO i = 1, klon IF (wk_adv(i)) THEN rre1(i) = 1. - sigmaw(i) rre2(i) = sigmaw(i) END IF END DO rrd1 = -1. rrd2 = 1. ! -- Get [Th1,Th2], dth and [q1,q2] DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN dth(i, k) = deltatw(i, k)/ppi(i, k) th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd d_th1(i, 1) = 0. d_th2(i, 1) = 0. d_dth(i, 1) = 0. d_q1(i, 1) = 0. d_q2(i, 1) = 0. d_dq(i, 1) = 0. END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN d_th1(i, k) = th1(i, k-1) - th1(i, k) d_th2(i, k) = th2(i, k-1) - th2(i, k) d_dth(i, k) = dth(i, k-1) - dth(i, k) d_q1(i, k) = q1(i, k-1) - q1(i, k) d_q2(i, k) = q2(i, k-1) - q2(i, k) d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN omgbdth(i, 1) = 0. omgbdq(i, 1) = 0. END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN ! loop on interfaces omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k)) omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k)) END IF END DO END DO ! ----------------------------------------------------------------- DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN ! ----------------------------------------------------------------- ! Compute redistribution (advective) term d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* & (rrd1*omg(i,k)*sigmaw(i)*d_th1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( & i))*d_th2(i,k+1)-(1.-alpha_up(i,k))*omgbdth(i,k)-alpha_up(i,k+1)* & omgbdth(i,k+1))*ppi(i, k) ! print*,'d_deltatw=',d_deltatw(i,k) d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* & (rrd1*omg(i,k)*sigmaw(i)*d_q1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( & i))*d_q2(i,k+1)-(1.-alpha_up(i,k))*omgbdq(i,k)-alpha_up(i,k+1)* & omgbdq(i,k+1)) ! print*,'d_deltaqw=',d_deltaqw(i,k) ! and increment large scale tendencies ! C ! ----------------------------------------------------------------- d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, & k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/(ph(i,k)-ph(i, & k+1)) & ! cc nrlmd $ ! -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k) -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/(ph(i, & k)-ph(i,k+1)) & ! cc )*ppi(i, k) d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, & k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/(ph(i,k)-ph(i, & k+1)) & ! cc nrlmd $ ! -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k) -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i, & k+1))/(ph(i,k)-ph(i,k+1)) & ! cc ) ! cc nrlmd ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, & k)/(ph(i,k)-ph(i,k+1))))*ppi(i, k) d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, & k)/(ph(i,k)-ph(i,k+1)))) END IF ! cc END DO END DO ! ------------------------------------------------------------------ ! Increment state variables DO k = 1, klev DO i = 1, klon ! cc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN IF (wk_adv(i) .AND. k<=kupper(i)) THEN ! cc ! Coefficient de répartition crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ & (ph(i,kupper(i))-ph(i,1)) crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-ph(i & ,kupper(i))) ! Reintroduce compensating subsidence term. ! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)) ! . /(1-sigmaw) ! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)) ! . /(1-sigmaw) ! dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k)) ! . /(1-sigmaw) ! dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k)) ! . /(1-sigmaw) dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i))) dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i))) ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k) dtpbl(i, k) = (wdtpbl(i,k)/sigmaw(i)-udtpbl(i,k)/(1.-sigmaw(i))) dqpbl(i, k) = (wdqpbl(i,k)/sigmaw(i)-udqpbl(i,k)/(1.-sigmaw(i))) ! print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k) ! cc nrlmd Prise en compte du taux de mortalité ! cc Définitions de entr, detr detr(i, k) = 0. entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + & sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k) spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i) ! cc spread(i,k) = ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/ ! cc $ sigmaw(i) ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU ! Jingmei ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k), ! & Tgw(i,k),deltatw(i,k) d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)*deltatw(i, k)* & dtimesub ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k) ff(i) = d_deltatw(i, k)/dtimesub ! Sans GW ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k)) ! GW formule 1 ! deltatw(k) = deltatw(k)+dtimesub* ! $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k)) ! GW formule 2 IF (dtimesub*tgw(i,k)<1.E-10) THEN d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc ! $ ! -spread(i,k)*deltatw(i,k) -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( & i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc -tgw(i,k)*deltatw(i,k)) ELSE d_deltatw(i, k) = 1/tgw(i, k)*(1-exp(-dtimesub*tgw(i, & k)))*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc $ ! -spread(i,k)*deltatw(i,k) -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( & i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc -tgw(i,k)*deltatw(i,k)) END IF dth(i, k) = deltatw(i, k)/ppi(i, k) gg(i) = d_deltaqw(i, k)/dtimesub d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k)+dqpbl(i,k) & ! cc $ ! -spread(i,k)*deltaqw(i,k)) -entr(i,k)*deltaqw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw(i)+detr( & i,k))*deltaqw(i,k)/(1.-sigmaw(i))) ! cc ! cc nrlmd ! cc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k) ! cc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k) ! cc END IF END DO END DO ! Scale tendencies so that water vapour remains positive in w and x. CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, & d_deltaqw, sigmaw, d_sigmaw, alpha) ! cc nrlmd ! c print*,'alpha' ! c do i=1,klon ! c print*,alpha(i) ! c end do ! cc DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN d_te(i, k) = alpha(i)*d_te(i, k) d_qe(i, k) = alpha(i)*d_qe(i, k) d_deltatw(i, k) = alpha(i)*d_deltatw(i, k) d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k) d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN d_sigmaw(i) = alpha(i)*d_sigmaw(i) END IF END DO ! Update large scale variables and wake variables ! IM 060208 manque DO i + remplace DO k=1,kupper(i) ! IM 060208 DO k = 1,kupper(i) DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN dtls(i, k) = dtls(i, k) + d_te(i, k) dqls(i, k) = dqls(i, k) + d_qe(i, k) ! cc nrlmd d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k) d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k) ! cc END IF END DO END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i) .AND. k<=kupper(i)) THEN te(i, k) = te0(i, k) + dtls(i, k) qe(i, k) = qe0(i, k) + dqls(i, k) the(i, k) = te(i, k)/ppi(i, k) deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k) deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k) dth(i, k) = deltatw(i, k)/ppi(i, k) ! c print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k) ! c $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN sigmaw(i) = sigmaw(i) + d_sigmaw(i) END IF END DO ! Determine Ptop from buoyancy integral ! --------------------------------------- ! - 1/ Pressure of the level where dth changes sign. DO i = 1, klon IF (wk_adv(i)) THEN ptop_provis(i) = ph(i, 1) END IF END DO DO k = 2, klev DO i = 1, klon IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. & dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO ! - 2/ dth integral DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd sum_dth(i) = 0. dthmin(i) = -delta_t_min z(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) dthmin(i) = amin1(dthmin(i), dth(i,k)) END IF END IF END DO END DO ! - 3/ height of triangle with area= sum_dth and base = dthmin DO i = 1, klon IF (wk_adv(i)) THEN hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) hw(i) = amax1(hwmin, hw(i)) END IF END DO ! - 4/ now, get Ptop DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd ktop(i) = 0 z(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i)) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i) ktop(i) = k END IF END IF END DO END DO ! 4.5/Correct ktop and ptop DO i = 1, klon IF (wk_adv(i)) THEN ptop_new(i) = ptop(i) END IF END DO DO k = klev, 2, -1 DO i = 1, klon ! IM v3JYG; IF (k .GE. ktop(i) IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN ptop(i) = ptop_new(i) END IF END DO DO k = klev, 1, -1 DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (ph(i,k+1)=kupper(i)) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. END IF END DO END DO ! -------------Cstar computation--------------------------------- DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd sum_thu(i) = 0. sum_tu(i) = 0. sum_qu(i) = 0. sum_thvu(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_rho(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thu(i) = 0. av_tu(i) = 0. av_qu(i) = 0. av_thvu(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_rho(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END IF END DO ! Integrals (and wake top level number) ! -------------------------------------- ! Initialize sum_thvu to 1st level virt. pot. temp. DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd z(i) = 1. dz(i) = 1. sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i) sum_dth(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) END IF END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd hw0(i) = z(i) END IF END DO ! - WAPE and mean forcing computation ! --------------------------------------- ! --------------------------------------- ! Means DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd av_thu(i) = sum_thu(i)/hw0(i) av_tu(i) = sum_tu(i)/hw0(i) av_qu(i) = sum_qu(i)/hw0(i) av_thvu(i) = sum_thvu(i)/hw0(i) av_dth(i) = sum_dth(i)/hw0(i) av_dq(i) = sum_dq(i)/hw0(i) av_rho(i) = sum_rho(i)/hw0(i) av_dtdwn(i) = sum_dtdwn(i)/hw0(i) av_dqdwn(i) = sum_dqdwn(i)/hw0(i) wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* & av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i) END IF END DO ! Filter out bad wakes DO k = 1, klev DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (wape(i)<0.) THEN deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. END IF END IF END DO END DO DO i = 1, klon IF (wk_adv(i)) THEN !!! nrlmd IF (wape(i)<0.) THEN wape(i) = 0. cstar(i) = 0. hw(i) = hwmin sigmaw(i) = max(sigmad, sigd_con(i)) fip(i) = 0. gwake(i) = .FALSE. ELSE cstar(i) = stark*sqrt(2.*wape(i)) gwake(i) = .TRUE. END IF END IF END DO END DO ! end sub-timestep loop ! ----------------------------------------------------------------- ! Get back to tendencies per second DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN ! cc dtls(i, k) = dtls(i, k)/dtime dqls(i, k) = dqls(i, k)/dtime d_deltatw2(i, k) = d_deltatw2(i, k)/dtime d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime ! c print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k) ! c $ ,death_rate(i)*sigmaw(i) END IF END DO END DO ! ---------------------------------------------------------- ! Determine wake final state; recompute wape, cstar, ktop; ! filter out bad wakes. ! ---------------------------------------------------------- ! 2.1 - Undisturbed area and Wake integrals ! --------------------------------------------------------- DO i = 1, klon ! cc nrlmd if (wk_adv(i)) then !!! nrlmd IF (ok_qx_qw(i)) THEN ! cc z(i) = 0. sum_thu(i) = 0. sum_tu(i) = 0. sum_qu(i) = 0. sum_thvu(i) = 0. sum_dth(i) = 0. sum_dq(i) = 0. sum_rho(i) = 0. sum_dtdwn(i) = 0. sum_dqdwn(i) = 0. av_thu(i) = 0. av_tu(i) = 0. av_qu(i) = 0. av_thvu(i) = 0. av_dth(i) = 0. av_dq(i) = 0. av_rho(i) = 0. av_dtdwn(i) = 0. av_dqdwn(i) = 0. END IF END DO ! Potential temperatures and humidity ! ---------------------------------------------------------- DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc rho(i, k) = p(i, k)/(rd*te(i,k)) IF (k==1) THEN rhoh(i, k) = ph(i, k)/(rd*te(i,k)) zhh(i, k) = 0 ELSE rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1))) zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1) END IF the(i, k) = te(i, k)/ppi(i, k) thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i) qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k))) dth(i, k) = deltatw(i, k)/ppi(i, k) END IF END DO END DO ! Integrals (and wake top level number) ! ----------------------------------------------------------- ! Initialize sum_thvu to 1st level virt. pot. temp. DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc z(i) = 1. dz(i) = 1. sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i) sum_dth(i) = 0. END IF END DO DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i)>0) THEN z(i) = z(i) + dz(i) sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i) sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i) sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i) sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i) sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) END IF END IF END DO END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc hw0(i) = z(i) END IF END DO ! - WAPE and mean forcing computation ! ------------------------------------------------------------- ! Means DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc av_thu(i) = sum_thu(i)/hw0(i) av_tu(i) = sum_tu(i)/hw0(i) av_qu(i) = sum_qu(i)/hw0(i) av_thvu(i) = sum_thvu(i)/hw0(i) av_dth(i) = sum_dth(i)/hw0(i) av_dq(i) = sum_dq(i)/hw0(i) av_rho(i) = sum_rho(i)/hw0(i) av_dtdwn(i) = sum_dtdwn(i)/hw0(i) av_dqdwn(i) = sum_dqdwn(i)/hw0(i) wape2(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* & av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i) END IF END DO ! Prognostic variable update ! ------------------------------------------------------------ ! Filter out bad wakes DO k = 1, klev DO i = 1, klon ! cc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN ! cc deltatw(i, k) = 0. deltaqw(i, k) = 0. dth(i, k) = 0. END IF END DO END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc IF (wape2(i)<0.) THEN wape2(i) = 0. cstar2(i) = 0. hw(i) = hwmin sigmaw(i) = amax1(sigmad, sigd_con(i)) fip(i) = 0. gwake(i) = .FALSE. ELSE IF (prt_level>=10) PRINT *, 'wape2>0' cstar2(i) = stark*sqrt(2.*wape2(i)) gwake(i) = .TRUE. END IF END IF END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc ktopw(i) = ktop(i) END IF END DO DO i = 1, klon ! cc nrlmd IF ( wk_adv(i)) THEN IF (ok_qx_qw(i)) THEN ! cc IF (ktopw(i)>0 .AND. gwake(i)) THEN ! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer) ! cc heff = 600. ! Utilisation de la hauteur hw ! c heff = 0.7*hw heff(i) = hw(i) fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* & sqrt(sigmaw(i)*wdens(i)*3.14) fip(i) = alpk*fip(i) ! jyg2 ELSE fip(i) = 0. END IF END IF END DO ! Limitation de sigmaw ! cc nrlmd ! DO i=1,klon ! IF (OK_qx_qw(i)) THEN ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max ! ENDIF ! ENDDO ! cc DO k = 1, klev DO i = 1, klon ! cc nrlmd On maintient désormais constant sigmaw en régime ! permanent ! cc IF ((sigmaw(i).GT.sigmaw_max).or. IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. & .NOT. ok_qx_qw(i)) THEN ! cc dtls(i, k) = 0. dqls(i, k) = 0. deltatw(i, k) = 0. deltaqw(i, k) = 0. END IF END DO END DO ! cc nrlmd On maintient désormais constant sigmaw en régime permanent DO i = 1, klon IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. & .NOT. ok_qx_qw(i)) THEN wape(i) = 0. cstar(i) = 0. hw(i) = hwmin sigmaw(i) = sigmad fip(i) = 0. ELSE wape(i) = wape2(i) cstar(i) = cstar2(i) END IF ! c print*,'wape wape2 ktopw OK_qx_qw =', ! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i) END DO RETURN END SUBROUTINE wake SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, & d_deltaqw, sigmaw, d_sigmaw, alpha) ! ------------------------------------------------------ ! Dtermination du coefficient alpha tel que les tendances ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent ! a une humidite positive dans la zone (x) et dans la zone (w). ! ------------------------------------------------------ ! Input REAL qe(nlon, nl), d_qe(nlon, nl) REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl) REAL sigmaw(nlon), d_sigmaw(nlon) LOGICAL wk_adv(nlon) INTEGER nl, nlon ! Output REAL alpha(nlon) ! Internal variables REAL zeta(nlon, nl) REAL alpha1(nlon) REAL x, a, b, c, discrim REAL epsilon ! DATA epsilon/1.e-15/ DO k = 1, nl DO i = 1, nlon IF (wk_adv(i)) THEN IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN zeta(i, k) = 0. ELSE zeta(i, k) = 1. END IF END IF END DO DO i = 1, nlon IF (wk_adv(i)) THEN x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + & (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i)*(deltaqw(i,k)+ & d_deltaqw(i,k)) a = -d_sigmaw(i)*d_deltaqw(i, k) b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - & deltaqw(i, k)*d_sigmaw(i) c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon discrim = b*b - 4.*a*c ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap alpha1(i) = 1. ELSE IF (x>=0.) THEN alpha1(i) = 1. ELSE IF (a>0.) THEN alpha1(i) = 0.9*min((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim & ))/(2.*a)) ELSE IF (a==0.) THEN alpha1(i) = 0.9*(-c/b) ELSE ! print*,'a,b,c discrim',a,b,c discrim alpha1(i) = 0.9*max((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim & ))/(2.*a)) END IF END IF END IF alpha(i) = min(alpha(i), alpha1(i)) END IF END DO END DO RETURN END SUBROUTINE wake_vec_modulation SUBROUTINE wake_scal(p, ph, ppi, dtime, sigd_con, te0, qe0, omgb, dtdwn, & dqdwn, amdwn, amup, dta, dqa, wdtpbl, wdqpbl, udtpbl, udqpbl, deltatw, & deltaqw, dth, hw, sigmaw, wape, fip, gfl, dtls, dqls, ktopw, omgbdth, & dp_omgb, wdens, tu, qu, dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, & spread, cstar, d_deltat_gw, d_deltatw2, d_deltaqw2) ! ************************************************************** ! * ! WAKE * ! retour a un Pupper fixe * ! * ! written by : GRANDPEIX Jean-Yves 09/03/2000 * ! modified by : ROEHRIG Romain 01/29/2007 * ! ************************************************************** USE dimphy 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. ! Variables d'etat : deltatw : ecart de temperature wake-undisturbed ! area ! deltaqw : ecart d'humidite wake-undisturbed area ! sigmaw : fraction d'aire occupee par la poche. ! 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 : hauteur de la poche ! dp_omgb : vertical gradient of large scale omega ! 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) ! spread : spreading term in dt_wake and dq_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 ! te0 : 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) ! 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) ! p : pressions aux milieux des couches (Pa) ! ph : pressions aux interfaces (Pa) ! ppi : (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 ! te : 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 : hw initial ! sigmaw0: sigmaw initial ! amflux : horizontal mass flux through wake boundary ! wdens : number of wakes per unit area (3D) or per ! unit length (2D) ! Tgw : 1 sur la période de onde de gravité ! Cgw : vitesse de propagation de onde de gravité ! LL : distance entre 2 poches ! ------------------------------------------------------------------------- ! Déclaration de variables ! ------------------------------------------------------------------------- include "dimensions.h" ! ccc include "dimphy.h" include "YOMCST.h" include "cvthermo.h" include "iniprint.h" ! Arguments en entree ! -------------------- REAL p(klev), ph(klev+1), ppi(klev) REAL dtime REAL te0(klev), qe0(klev) REAL omgb(klev+1) REAL dtdwn(klev), dqdwn(klev) REAL wdtpbl(klev), wdqpbl(klev) REAL udtpbl(klev), udqpbl(klev) REAL amdwn(klev), amup(klev) REAL dta(klev), dqa(klev) REAL sigd_con ! Sorties ! -------- REAL deltatw(klev), deltaqw(klev), dth(klev) REAL tu(klev), qu(klev) REAL dtls(klev), dqls(klev) REAL dtke(klev), dqke(klev) REAL dtpbl(klev), dqpbl(klev) REAL spread(klev) REAL d_deltatgw(klev) REAL d_deltatw2(klev), d_deltaqw2(klev) REAL omgbdth(klev+1), omg(klev+1) REAL dp_omgb(klev), dp_deltomg(klev) REAL d_deltat_gw(klev) REAL hw, sigmaw, wape, fip, gfl, cstar INTEGER ktopw ! Variables internes ! ------------------- ! Variables à fixer REAL alon REAL coefgw REAL wdens0, wdens REAL stark REAL alpk REAL delta_t_min REAL pupper INTEGER nsub REAL dtimesub REAL sigmad, hwmin ! Variables de sauvegarde REAL deltatw0(klev) REAL deltaqw0(klev) REAL te(klev), qe(klev) REAL sigmaw0, sigmaw1 ! Variables pour les GW REAL ll REAL n2(klev) REAL cgw(klev) REAL tgw(klev) ! Variables liées au calcul de hw REAL ptop_provis, ptop, ptop_new REAL sum_dth REAL dthmin REAL z, dz, hw0 INTEGER ktop, kupper ! Autres variables internes INTEGER isubstep, k REAL sum_thu, sum_tu, sum_qu, sum_thvu REAL sum_dq, sum_rho REAL sum_dtdwn, sum_dqdwn REAL av_thu, av_tu, av_qu, av_thvu REAL av_dth, av_dq, av_rho REAL av_dtdwn, av_dqdwn REAL rho(klev), rhoh(klev+1), rhow(klev) REAL rhow_moyen(klev) REAL zh(klev), zhh(klev+1) REAL epaisseur1(klev), epaisseur2(klev) REAL the(klev), thu(klev) REAL d_deltatw(klev), d_deltaqw(klev) REAL omgbw(klev+1), omgtop REAL dp_omgbw(klev) REAL ztop, dztop REAL alpha_up(klev) REAL rre1, rre2, rrd1, rrd2 REAL th1(klev), th2(klev), q1(klev), q2(klev) REAL d_th1(klev), d_th2(klev), d_dth(klev) REAL d_q1(klev), d_q2(klev), d_dq(klev) REAL omgbdq(klev) REAL ff, gg REAL wape2, cstar2, heff REAL crep(klev) REAL crep_upper, crep_sol ! ------------------------------------------------------------------------- ! Initialisations ! ------------------------------------------------------------------------- ! print*, 'wake initialisations' ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10. ! ------------------------------------------------------------------------- DATA sigmad, hwmin/.02, 10./ ! Longueur de maille (en m) ! ------------------------------------------------------------------------- ! ALON = 3.e5 alon = 1.E6 ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei) ! coefgw : Coefficient pour les ondes de gravité ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) ! wdens : Densité de poche froide par maille ! ------------------------------------------------------------------------- coefgw = 10 ! coefgw=1 ! wdens0 = 1.0/(alon**2) wdens = 1.0/(alon**2) stark = 0.50 ! CRtest alpk = 0.1 ! alpk = 1.0 ! alpk = 0.5 ! alpk = 0.05 crep_upper = 0.9 crep_sol = 1.0 ! Minimum value for |T_wake - T_undist|. Used for wake top definition ! ------------------------------------------------------------------------- delta_t_min = 0.2 ! 1. - Save initial values and initialize tendencies ! -------------------------------------------------- DO k = 1, klev deltatw0(k) = deltatw(k) deltaqw0(k) = deltaqw(k) te(k) = te0(k) qe(k) = qe0(k) dtls(k) = 0. dqls(k) = 0. d_deltat_gw(k) = 0. d_deltatw2(k) = 0. d_deltaqw2(k) = 0. END DO ! sigmaw1=sigmaw ! IF (sigd_con.GT.sigmaw1) THEN ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con ! ENDIF sigmaw = max(sigmaw, sigd_con) sigmaw = max(sigmaw, sigmad) sigmaw = min(sigmaw, 0.99) sigmaw0 = sigmaw ! wdens=wdens0/(10.*sigmaw) ! IF (sigd_con.GT.sigmaw1) THEN ! print*, 'sigmaw1,sigd1', sigmaw, sigd_con ! ENDIF ! 2. - Prognostic part ! ========================================================= ! print *, 'prognostic wake computation' ! 2.1 - Undisturbed area and Wake integrals ! --------------------------------------------------------- z = 0. ktop = 0 kupper = 0 sum_thu = 0. sum_tu = 0. sum_qu = 0. sum_thvu = 0. sum_dth = 0. sum_dq = 0. sum_rho = 0. sum_dtdwn = 0. sum_dqdwn = 0. av_thu = 0. av_tu = 0. av_qu = 0. av_thvu = 0. av_dth = 0. av_dq = 0. av_rho = 0. av_dtdwn = 0. av_dqdwn = 0. ! Potential temperatures and humidity ! ---------------------------------------------------------- DO k = 1, klev rho(k) = p(k)/(rd*te(k)) IF (k==1) THEN rhoh(k) = ph(k)/(rd*te(k)) zhh(k) = 0 ELSE rhoh(k) = ph(k)*2./(rd*(te(k)+te(k-1))) zhh(k) = (ph(k)-ph(k-1))/(-rhoh(k)*rg) + zhh(k-1) END IF the(k) = te(k)/ppi(k) thu(k) = (te(k)-deltatw(k)*sigmaw)/ppi(k) tu(k) = te(k) - deltatw(k)*sigmaw qu(k) = qe(k) - deltaqw(k)*sigmaw rhow(k) = p(k)/(rd*(te(k)+deltatw(k))) dth(k) = deltatw(k)/ppi(k) ll = (1-sqrt(sigmaw))/sqrt(wdens) END DO DO k = 1, klev - 1 IF (k==1) THEN n2(k) = 0 ELSE n2(k) = max(0., -rg**2/the(k)*rho(k)*(the(k+1)-the(k-1))/(p(k+ & 1)-p(k-1))) END IF zh(k) = (zhh(k)+zhh(k+1))/2 cgw(k) = sqrt(n2(k))*zh(k) tgw(k) = coefgw*cgw(k)/ll END DO n2(klev) = 0 zh(klev) = 0 cgw(klev) = 0 tgw(klev) = 0 ! Calcul de la masse volumique moyenne de la colonne ! ----------------------------------------------------------------- DO k = 1, klev epaisseur1(k) = 0. epaisseur2(k) = 0. END DO epaisseur1(1) = -(ph(2)-ph(1))/(rho(1)*rg) + 1. epaisseur2(1) = -(ph(2)-ph(1))/(rho(1)*rg) + 1. rhow_moyen(1) = rhow(1) DO k = 2, klev epaisseur1(k) = -(ph(k+1)-ph(k))/(rho(k)*rg) + 1. epaisseur2(k) = epaisseur2(k-1) + epaisseur1(k) rhow_moyen(k) = (rhow_moyen(k-1)*epaisseur2(k-1)+rhow(k)*epaisseur1(k))/ & epaisseur2(k) END DO ! Choose an integration bound well above wake top ! ----------------------------------------------------------------- ! Pupper = 50000. ! melting level pupper = 60000. ! Pupper = 70000. ! Determine Wake top pressure (Ptop) from buoyancy integral ! ----------------------------------------------------------------- ! -1/ Pressure of the level where dth becomes less than delta_t_min. ptop_provis = ph(1) DO k = 2, klev IF (dth(k)>-delta_t_min .AND. dth(k-1)<-delta_t_min) THEN ptop_provis = ((dth(k)+delta_t_min)*p(k-1)-(dth(k- & 1)+delta_t_min)*p(k))/(dth(k)-dth(k-1)) GO TO 25 END IF END DO 25 CONTINUE ! -2/ dth integral sum_dth = 0. dthmin = -delta_t_min z = 0. DO k = 1, klev dz = -(max(ph(k+1),ptop_provis)-ph(k))/(rho(k)*rg) IF (dz<=0) GO TO 40 z = z + dz sum_dth = sum_dth + dth(k)*dz dthmin = min(dthmin, dth(k)) END DO 40 CONTINUE ! -3/ height of triangle with area= sum_dth and base = dthmin hw0 = 2.*sum_dth/min(dthmin, -0.5) hw0 = max(hwmin, hw0) ! -4/ now, get Ptop z = 0. ptop = ph(1) DO k = 1, klev dz = min(-(ph(k+1)-ph(k))/(rho(k)*rg), hw0-z) IF (dz<=0) GO TO 45 z = z + dz ptop = ph(k) - rho(k)*rg*dz END DO 45 CONTINUE ! -5/ Determination de ktop et kupper DO k = klev, 1, -1 IF (ph(k+1)-delta_t_min .AND. dth(k-1)<-delta_t_min) THEN ptop_new = ((dth(k)+delta_t_min)*p(k-1)-(dth(k-1)+delta_t_min)*p(k))/ & (dth(k)-dth(k-1)) GO TO 225 END IF END DO 225 CONTINUE ptop = ptop_new DO k = klev, 1, -1 IF (ph(k+1)=10) PRINT *, 'wape<0' wape = 0. hw = hwmin sigmaw = max(sigmad, sigd_con) fip = 0. DO k = 1, klev deltatw(k) = 0. deltaqw(k) = 0. dth(k) = 0. END DO ELSE IF (prt_level>=10) PRINT *, 'wape>0' cstar = stark*sqrt(2.*wape) END IF ! ------------------------------------------------------------------ ! Sub-time-stepping ! ------------------------------------------------------------------ ! nsub=36 nsub = 10 dtimesub = dtime/nsub ! ------------------------------------------------------------ DO isubstep = 1, nsub ! ------------------------------------------------------------ ! print*,'---------------','substep=',isubstep,'-------------' ! Evolution of sigmaw gfl = 2.*sqrt(3.14*wdens*sigmaw) sigmaw = sigmaw + gfl*cstar*dtimesub sigmaw = min(sigmaw, 0.99) !!!!!!!! ! wdens = wdens0/(10.*sigmaw) ! sigmaw =max(sigmaw,sigd_con) ! sigmaw =max(sigmaw,sigmad) ! calcul de la difference de vitesse verticale poche - zone non perturbee z = 0. dp_deltomg(1:klev) = 0. omg(1:klev+1) = 0. omg(1) = 0. dp_deltomg(1) = -(gfl*cstar)/(sigmaw*(1-sigmaw)) DO k = 2, ktop dz = -(ph(k)-ph(k-1))/(rho(k-1)*rg) z = z + dz dp_deltomg(k) = dp_deltomg(1) omg(k) = dp_deltomg(1)*z END DO dztop = -(ptop-ph(ktop))/(rho(ktop)*rg) ztop = z + dztop omgtop = dp_deltomg(1)*ztop ! Conversion de la vitesse verticale de m/s a Pa/s omgtop = -rho(ktop)*rg*omgtop dp_deltomg(1) = omgtop/(ptop-ph(1)) DO k = 1, ktop omg(k) = -rho(k)*rg*omg(k) dp_deltomg(k) = dp_deltomg(1) END DO ! raccordement lineaire de omg de ptop a pupper IF (kupper>ktop) THEN omg(kupper+1) = -rg*amdwn(kupper+1)/sigmaw + rg*amup(kupper+1)/(1.- & sigmaw) dp_deltomg(kupper) = (omgtop-omg(kupper+1))/(ptop-pupper) DO k = ktop + 1, kupper dp_deltomg(k) = dp_deltomg(kupper) omg(k) = omgtop + (ph(k)-ptop)*dp_deltomg(kupper) END DO END IF ! Compute wake average vertical velocity omgbw DO k = 1, klev + 1 omgbw(k) = omgb(k) + (1.-sigmaw)*omg(k) END DO ! and its vertical gradient dp_omgbw DO k = 1, klev dp_omgbw(k) = (omgbw(k+1)-omgbw(k))/(ph(k+1)-ph(k)) 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 alpha_up(k) = 0. IF (omgb(k)>0.) alpha_up(k) = 1. END DO ! Matrix expressing [The,deltatw] from [Th1,Th2] rre1 = 1. - sigmaw rre2 = sigmaw rrd1 = -1. rrd2 = 1. ! Get [Th1,Th2], dth and [q1,q2] DO k = 1, kupper + 1 dth(k) = deltatw(k)/ppi(k) th1(k) = the(k) - sigmaw*dth(k) ! undisturbed area th2(k) = the(k) + (1.-sigmaw)*dth(k) ! wake q1(k) = qe(k) - sigmaw*deltaqw(k) ! undisturbed area q2(k) = qe(k) + (1.-sigmaw)*deltaqw(k) ! wake END DO d_th1(1) = 0. d_th2(1) = 0. d_dth(1) = 0. d_q1(1) = 0. d_q2(1) = 0. d_dq(1) = 0. DO k = 2, kupper + 1 ! loop on interfaces d_th1(k) = th1(k-1) - th1(k) d_th2(k) = th2(k-1) - th2(k) d_dth(k) = dth(k-1) - dth(k) d_q1(k) = q1(k-1) - q1(k) d_q2(k) = q2(k-1) - q2(k) d_dq(k) = deltaqw(k-1) - deltaqw(k) END DO omgbdth(1) = 0. omgbdq(1) = 0. DO k = 2, kupper + 1 ! loop on interfaces omgbdth(k) = omgb(k)*(dth(k-1)-dth(k)) omgbdq(k) = omgb(k)*(deltaqw(k-1)-deltaqw(k)) END DO ! ----------------------------------------------------------------- DO k = 1, kupper - 1 ! ----------------------------------------------------------------- ! Compute redistribution (advective) term d_deltatw(k) = dtimesub/(ph(k)-ph(k+1))*(rrd1*omg(k)*sigmaw*d_th1(k)- & rrd2*omg(k+1)*(1.-sigmaw)*d_th2(k+1)-(1.-alpha_up( & k))*omgbdth(k)-alpha_up(k+1)*omgbdth(k+1))*ppi(k) ! print*,'d_deltatw=',d_deltatw(k) d_deltaqw(k) = dtimesub/(ph(k)-ph(k+1))*(rrd1*omg(k)*sigmaw*d_q1(k)- & rrd2*omg(k+1)*(1.-sigmaw)*d_q2(k+1)-(1.-alpha_up( & k))*omgbdq(k)-alpha_up(k+1)*omgbdq(k+1)) ! print*,'d_deltaqw=',d_deltaqw(k) ! and increment large scale tendencies dtls(k) = dtls(k) + dtimesub*((rre1*omg(k)*sigmaw*d_th1(k)-rre2*omg(k+ & 1)*(1.-sigmaw)*d_th2(k+1))/(ph(k)-ph(k+1))-sigmaw*(1.-sigmaw)*dth(k)* & dp_deltomg(k))*ppi(k) ! print*,'dtls=',dtls(k) dqls(k) = dqls(k) + dtimesub*((rre1*omg(k)*sigmaw*d_q1(k)-rre2*omg(k+ & 1)*(1.-sigmaw)*d_q2(k+1))/(ph(k)-ph(k+1))-sigmaw*(1.-sigmaw)*deltaqw( & k)*dp_deltomg(k)) ! print*,'dqls=',dqls(k) ! ------------------------------------------------------------------- END DO ! ------------------------------------------------------------------ ! Increment state variables DO k = 1, kupper - 1 ! Coefficient de répartition crep(k) = crep_sol*(ph(kupper)-ph(k))/(ph(kupper)-ph(1)) crep(k) = crep(k) + crep_upper*(ph(1)-ph(k))/(p(1)-ph(kupper)) ! 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(k) = (dtdwn(k)/sigmaw-dta(k)/(1.-sigmaw)) dqke(k) = (dqdwn(k)/sigmaw-dqa(k)/(1.-sigmaw)) ! print*,'dtKE=',dtKE(k) ! print*,'dqKE=',dqKE(k) dtpbl(k) = (wdtpbl(k)/sigmaw-udtpbl(k)/(1.-sigmaw)) dqpbl(k) = (wdqpbl(k)/sigmaw-udqpbl(k)/(1.-sigmaw)) spread(k) = (1.-sigmaw)*dp_deltomg(k) + gfl*cstar/sigmaw ! print*,'spread=',spread(k) ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU ! Jingmei d_deltat_gw(k) = d_deltat_gw(k) - tgw(k)*deltatw(k)*dtimesub ! print*,'d_delta_gw=',d_deltat_gw(k) ff = d_deltatw(k)/dtimesub ! Sans GW ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k)) ! GW formule 1 ! deltatw(k) = deltatw(k)+dtimesub* ! $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k)) ! GW formule 2 IF (dtimesub*tgw(k)<1.E-10) THEN deltatw(k) = deltatw(k) + dtimesub*(ff+dtke(k)+dtpbl(k)-spread(k)* & deltatw(k)-tgw(k)*deltatw(k)) ELSE deltatw(k) = deltatw(k) + 1/tgw(k)*(1-exp(-dtimesub*tgw(k)))*(ff+dtke & (k)+dtpbl(k)-spread(k)*deltatw(k)-tgw(k)*deltatw(k)) END IF dth(k) = deltatw(k)/ppi(k) gg = d_deltaqw(k)/dtimesub deltaqw(k) = deltaqw(k) + dtimesub*(gg+dqke(k)+dqpbl(k)-spread(k)* & deltaqw(k)) d_deltatw2(k) = d_deltatw2(k) + d_deltatw(k) d_deltaqw2(k) = d_deltaqw2(k) + d_deltaqw(k) END DO ! And update large scale variables DO k = 1, kupper te(k) = te0(k) + dtls(k) qe(k) = qe0(k) + dqls(k) the(k) = te(k)/ppi(k) END DO ! Determine Ptop from buoyancy integral ! ---------------------------------------------------------------------- ! -1/ Pressure of the level where dth changes sign. ptop_provis = ph(1) DO k = 2, klev IF (dth(k)>-delta_t_min .AND. dth(k-1)<-delta_t_min) THEN ptop_provis = ((dth(k)+delta_t_min)*p(k-1)-(dth(k- & 1)+delta_t_min)*p(k))/(dth(k)-dth(k-1)) GO TO 65 END IF END DO 65 CONTINUE ! -2/ dth integral sum_dth = 0. dthmin = -delta_t_min z = 0. DO k = 1, klev dz = -(max(ph(k+1),ptop_provis)-ph(k))/(rho(k)*rg) IF (dz<=0) GO TO 70 z = z + dz sum_dth = sum_dth + dth(k)*dz dthmin = min(dthmin, dth(k)) END DO 70 CONTINUE ! -3/ height of triangle with area= sum_dth and base = dthmin hw = 2.*sum_dth/min(dthmin, -0.5) hw = max(hwmin, hw) ! -4/ now, get Ptop ktop = 0 z = 0. DO k = 1, klev dz = min(-(ph(k+1)-ph(k))/(rho(k)*rg), hw-z) IF (dz<=0) GO TO 75 z = z + dz ptop = ph(k) - rho(k)*rg*dz ktop = k END DO 75 CONTINUE ! -5/Correct ktop and ptop ptop_new = ptop DO k = ktop, 2, -1 IF (dth(k)>-delta_t_min .AND. dth(k-1)<-delta_t_min) THEN ptop_new = ((dth(k)+delta_t_min)*p(k-1)-(dth(k-1)+delta_t_min)*p(k))/ & (dth(k)-dth(k-1)) GO TO 275 END IF END DO 275 CONTINUE ptop = ptop_new DO k = klev, 1, -1 IF (ph(k+1)=10) PRINT *, 'wape2<0' wape2 = 0. hw = hwmin sigmaw = max(sigmad, sigd_con) fip = 0. DO k = 1, klev deltatw(k) = 0. deltaqw(k) = 0. dth(k) = 0. END DO ELSE IF (prt_level>=10) PRINT *, 'wape2>0' cstar2 = stark*sqrt(2.*wape2) END IF ktopw = ktop IF (ktopw>0) THEN ! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer) ! cc heff = 600. ! Utilisation de la hauteur hw ! c heff = 0.7*hw heff = hw fip = 0.5*rho(ktopw)*cstar2**3*heff*2*sqrt(sigmaw*wdens*3.14) fip = alpk*fip ! jyg2 ELSE fip = 0. END IF ! Limitation de sigmaw ! sécurité : si le wake occuppe plus de 90 % de la surface de la maille, ! alors il disparait en se mélangeant à la partie undisturbed ! correction NICOLAS . ((wape.ge.wape2).and.(wape2.le.1.0))) THEN IF ((sigmaw>0.9) .OR. ((wape>=wape2) .AND. (wape2<= & 1.0)) .OR. (ktopw<=2)) THEN ! IM cf NR/JYG 251108 . ((wape.ge.wape2).and.(wape2.le.1.0))) THEN ! IF (sigmaw.GT.0.9) THEN DO k = 1, klev dtls(k) = 0. dqls(k) = 0. deltatw(k) = 0. deltaqw(k) = 0. END DO wape = 0. hw = hwmin sigmaw = sigmad fip = 0. END IF RETURN END SUBROUTINE wake_scal