Subroutine WAKE (p,ph,ppi,dtime,sigd_con : ,te0,qe0,omgb,ibas : ,dtdwn,dqdwn,amdwn,amup,dta,dqa : ,wdtPBL,wdqPBL,udtPBL,udqPBL o ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl o ,dtls,dqls o ,ktopw,omgbdth,dp_omgb,wdens o ,tu,qu o ,dtKE,dqKE o ,dtPBL,dqPBL o ,omg,dp_deltomg,spread o ,Cstar,d_deltat_gw o ,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 * *************************************************************** c USE dimphy IMPLICIT none c============================================================================ C C C But : Decrire le comportement des poches froides apparaissant dans les C grands systemes convectifs, et fournir l'energie disponible pour C le declenchement de nouvelles colonnes convectives. C C Variables d'etat : deltatw : ecart de temperature wake-undisturbed area C deltaqw : ecart d'humidite wake-undisturbed area C sigmaw : fraction d'aire occupee par la poche. C C Variable de sortie : c c wape : WAke Potential Energy c fip : Front Incident Power (W/m2) - ALP c gfl : Gust Front Length per unit area (m-1) C dtls : large scale temperature tendency due to wake C dqls : large scale humidity tendency due to wake C hw : hauteur de la poche C dp_omgb : vertical gradient of large scale omega C omgbdth: flux of Delta_Theta transported by LS omega C dtKE : differential heating (wake - unpertubed) C dqKE : differential moistening (wake - unpertubed) C omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s) C dp_deltomg : vertical gradient of omg (s-1) C spread : spreading term in dt_wake and dq_wake C deltatw : updated temperature difference (T_w-T_u). C deltaqw : updated humidity difference (q_w-q_u). C sigmaw : updated wake fractional area. C d_deltat_gw : delta T tendency due to GW c C Variables d'entree : c c aire : aire de la maille c te0 : temperature dans l'environnement (K) C qe0 : humidite dans l'environnement (kg/kg) C omgb : vitesse verticale moyenne sur la maille (Pa/s) C ibas : cloud base level number C dtdwn: source de chaleur due aux descentes (K/s) C dqdwn: source d'humidite due aux descentes (kg/kg/s) C dta : source de chaleur due courants satures et detrain (K/s) C dqa : source d'humidite due aux courants satures et detra (kg/kg/s) C amdwn: flux de masse total des descentes, par unite de C surface de la maille (kg/m2/s) C amup : flux de masse total des ascendances, par unite de C surface de la maille (kg/m2/s) C p : pressions aux milieux des couches (Pa) C ph : pressions aux interfaces (Pa) C ppi : (p/p_0)**kapa (adim) C dtime: increment temporel (s) c C Variables internes : c c rhow : masse volumique de la poche froide C rho : environment density at P levels C rhoh : environment density at Ph levels C te : environment temperature | may change within C qe : environment humidity | sub-time-stepping C the : environment potential temperature C thu : potential temperature in undisturbed area C tu : temperature in undisturbed area C qu : humidity in undisturbed area C dp_omgb: vertical gradient og LS omega C omgbw : wake average vertical omega C dp_omgbw: vertical gradient of omgbw C omgbdq : flux of Delta_q transported by LS omega C dth : potential temperature diff. wake-undist. C th1 : first pot. temp. for vertical advection (=thu) C th2 : second pot. temp. for vertical advection (=thw) C q1 : first humidity for vertical advection C q2 : second humidity for vertical advection C d_deltatw : terme de redistribution pour deltatw C d_deltaqw : terme de redistribution pour deltaqw C deltatw0 : deltatw initial C deltaqw0 : deltaqw initial C hw0 : hw initial C sigmaw0: sigmaw initial C amflux : horizontal mass flux through wake boundary C wdens : number of wakes per unit area (3D) or per C unit length (2D) C Tgw : 1 sur la période de onde de gravité c Cgw : vitesse de propagation de onde de gravité c LL : distance entre 2 poches c------------------------------------------------------------------------- c Déclaration de variables c------------------------------------------------------------------------- #include "dimensions.h" cccc#include "dimphy.h" #include "YOMCST.h" #include "cvthermo.h" c Arguments en entree c-------------------- REAL p(klev),ph(klev+1),ppi(klev) REAL dtime REAL te0(klev),qe0(klev) REAL omgb(klev+1) INTEGER ibas 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 c Sorties c-------- 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 c Variables internes c------------------- c 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 c Variables de sauvegarde REAL deltatw0(klev) REAL deltaqw0(klev) REAL te(klev), qe(klev) REAL sigmaw0, sigmaw1 c Variables pour les GW REAL LL REAL N2(klev) REAL Cgw(klev) REAL Tgw(klev) c 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 c 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 pbase 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 C------------------------------------------------------------------------- c Initialisations c------------------------------------------------------------------------- c print*, 'wake initialisations' c Essais d'initialisation avec sigmaw = 0.02 et hw = 10. c------------------------------------------------------------------------- DATA sigmad, hwmin /.02,10./ C Longueur de maille (en m) c------------------------------------------------------------------------- c ALON = 3.e5 ALON = 1.e6 C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei) c c coefgw : Coefficient pour les ondes de gravité c stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) c wdens : Densité de poche froide par maille c------------------------------------------------------------------------- coefgw=10 c coefgw=1 c wdens0 = 1.0/(alon**2) wdens = 1.0/(alon**2) stark = 0.50 cCRtest alpk=0.1 c alpk = 1.0 c alpk = 0.5 c alpk = 0.05 Crep_upper=0.9 Crep_sol=1.0 C Minimum value for |T_wake - T_undist|. Used for wake top definition c------------------------------------------------------------------------- delta_t_min = 0.2 C Cloud base c------------------------------------------------------------------------- Pbase = P(ibas) C 1. - Save initial values and initialize tendencies C -------------------------------------------------- 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. ENDDO c sigmaw1=sigmaw c IF (sigd_con.GT.sigmaw1) THEN c print*, 'sigmaw,sigd_con', sigmaw, sigd_con c ENDIF sigmaw = max(sigmaw,sigd_con) sigmaw = max(sigmaw,sigmad) sigmaw = min(sigmaw,0.99) sigmaw0 = sigmaw c wdens=wdens0/(10.*sigmaw) c IF (sigd_con.GT.sigmaw1) THEN c print*, 'sigmaw1,sigd1', sigmaw, sigd_con c ENDIF C 2. - Prognostic part C ========================================================= c print *, 'prognostic wake computation' C 2.1 - Undisturbed area and Wake integrals C --------------------------------------------------------- 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. C Potential temperatures and humidity c---------------------------------------------------------- DO k =1,klev rho(k) = p(k)/(rd*te(k)) IF(k .eq. 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) ENDIF 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) ENDDO DO k = 1, klev-1 IF(k.eq.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))) ENDIF ZH(k)=(zhh(k)+zhh(k+1))/2 Cgw(k)=sqrt(N2(k))*ZH(k) Tgw(k)=coefgw*Cgw(k)/LL ENDDO N2(klev)=0 ZH(klev)=0 Cgw(klev)=0 Tgw(klev)=0 c Calcul de la masse volumique moyenne de la colonne c----------------------------------------------------------------- DO k=1,klev epaisseur1(k)=0. epaisseur2(k)=0. ENDDO 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) ENDDO C Choose an integration bound well above wake top c----------------------------------------------------------------- c Pupper = 50000. ! melting level Pupper = 60000. c Pupper = 70000. C Determine Wake top pressure (Ptop) from buoyancy integral C----------------------------------------------------------------- c-1/ Pressure of the level where dth becomes less than delta_t_min. Ptop_provis=ph(1) DO k= 2,klev IF (dth(k) .GT. -delta_t_min .and. $ dth(k-1).LT. -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 ENDIF ENDDO 25 CONTINUE c-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 .le. 0) GO TO 40 z = z+dz sum_dth = sum_dth + dth(k)*dz dthmin = min(dthmin,dth(k)) ENDDO 40 CONTINUE c-3/ height of triangle with area= sum_dth and base = dthmin hw0 = 2.*sum_dth/min(dthmin,-0.5) hw0 = max(hwmin,hw0) c-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 .le. 0) GO TO 45 z = z+dz Ptop = Ph(k)-rho(k)*rg*dz ENDDO 45 CONTINUE C-5/ Determination de ktop et kupper DO k=klev,1,-1 IF (ph(k+1) .lt. ptop) ktop=k IF (ph(k+1) .lt. pupper) kupper=k ENDDO c-6/ Correct ktop and ptop Ptop_new=ptop DO k= ktop,2,-1 IF (dth(k) .GT. -delta_t_min .and. $ dth(k-1).LT. -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 ENDIF ENDDO 225 CONTINUE ptop = ptop_new DO k=klev,1,-1 IF (ph(k+1) .lt. ptop) ktop=k ENDDO c Set deltatw & deltaqw to 0 above kupper c----------------------------------------------------------- DO k = kupper,klev deltatw(k) = 0. deltaqw(k) = 0. ENDDO C Vertical gradient of LS omega C------------------------------------------------------------ DO k = 1,kupper dp_omgb(k) = (omgb(k+1) - omgb(k))/(ph(k+1)-ph(k)) ENDDO C Integrals (and wake top level number) C ----------------------------------------------------------- C Initialize sum_thvu to 1st level virt. pot. temp. z = 1. dz = 1. sum_thvu = thu(1)*(1.+eps*qu(1))*dz sum_dth = 0. DO k = 1,klev dz = -(max(ph(k+1),Ptop)-Ph(k))/(rho(k)*rg) IF (dz .LE. 0) GO TO 50 z = z+dz sum_thu = sum_thu + thu(k)*dz sum_tu = sum_tu + tu(k)*dz sum_qu = sum_qu + qu(k)*dz sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz sum_dth = sum_dth + dth(k)*dz sum_dq = sum_dq + deltaqw(k)*dz sum_rho = sum_rho + rhow(k)*dz sum_dtdwn = sum_dtdwn + dtdwn(k)*dz sum_dqdwn = sum_dqdwn + dqdwn(k)*dz ENDDO 50 CONTINUE hw0 = z C 2.1 - WAPE and mean forcing computation C------------------------------------------------------------- C Means av_thu = sum_thu/hw0 av_tu = sum_tu/hw0 av_qu = sum_qu/hw0 av_thvu = sum_thvu/hw0 c av_thve = sum_thve/hw0 av_dth = sum_dth/hw0 av_dq = sum_dq/hw0 av_rho = sum_rho/hw0 av_dtdwn = sum_dtdwn/hw0 av_dqdwn = sum_dqdwn/hw0 wape = - rg*hw0*(av_dth $ + eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq ))/av_thvu C 2.2 Prognostic variable update C ------------------------------------------------------------ C Filter out bad wakes IF ( wape .LT. 0.) THEN 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. ENDDO ELSE print*,'wape>0' Cstar = stark*sqrt(2.*wape) ENDIF C------------------------------------------------------------------ C Sub-time-stepping C------------------------------------------------------------------ c nsub=36 nsub=10 dtimesub=dtime/nsub c------------------------------------------------------------ DO isubstep = 1,nsub c------------------------------------------------------------ c print*,'---------------','substep=',isubstep,'-------------' c Evolution of sigmaw gfl = 2.*sqrt(3.14*wdens*sigmaw) sigmaw =sigmaw + gfl*Cstar*dtimesub sigmaw =min(sigmaw,0.99) !!!!!!!! c wdens = wdens0/(10.*sigmaw) c sigmaw =max(sigmaw,sigd_con) c sigmaw =max(sigmaw,sigmad) c 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 ENDDO dztop=-(Ptop-Ph(ktop))/(rho(ktop)*rg) ztop = z+dztop omgtop=dp_deltomg(1)*ztop c 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) ENDDO c raccordement lineaire de omg de ptop a pupper IF (kupper .GT. 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) ENDDO ENDIF c Compute wake average vertical velocity omgbw DO k = 1,klev+1 omgbw(k) = omgb(k)+(1.-sigmaw)*omg(k) ENDDO c and its vertical gradient dp_omgbw DO k = 1,klev dp_omgbw(k) = (omgbw(k+1)-omgbw(k))/(ph(k+1)-ph(k)) ENDDO c Upstream coefficients for omgb velocity c-- (alpha_up(k) is the coefficient of the value at level k) c-- (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) .GT. 0.) alpha_up(k) = 1. ENDDO c Matrix expressing [The,deltatw] from [Th1,Th2] RRe1 = 1.-sigmaw RRe2 = sigmaw RRd1 = -1. RRd2 = 1. c 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 ENDDO 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) ENDDO 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)) ENDDO c----------------------------------------------------------------- DO k=1,kupper-1 c----------------------------------------------------------------- c c Compute redistribution (advective) term c 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) c print*,'d_deltatw=',d_deltatw(k) c 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) $ ) c print*,'d_deltaqw=',d_deltaqw(k) c c and increment large scale tendencies c 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) c print*,'dtls=',dtls(k) c 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) $ ) c print*,'dqls=',dqls(k) c------------------------------------------------------------------- ENDDO c------------------------------------------------------------------ C Increment state variables DO k = 1,kupper-1 c 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)) c Reintroduce compensating subsidence term. c dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)) c . /(1-sigmaw) c dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)) c . /(1-sigmaw) c c dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k)) c . /(1-sigmaw) c dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k)) c . /(1-sigmaw) dtKE(k)=(dtdwn(k)/sigmaw - dta(k)/(1.-sigmaw)) dqKE(k)=(dqdwn(k)/sigmaw - dqa(k)/(1.-sigmaw)) c print*,'dtKE=',dtKE(k) c print*,'dqKE=',dqKE(k) c dtPBL(k)=(wdtPBL(k)/sigmaw - udtPBL(k)/(1.-sigmaw)) dqPBL(k)=(wdqPBL(k)/sigmaw - udqPBL(k)/(1.-sigmaw)) c spread(k) = (1.-sigmaw)*dp_deltomg(k)+gfl*Cstar/sigmaw c print*,'spread=',spread(k) c 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 c print*,'d_delta_gw=',d_deltat_gw(k) ff=d_deltatw(k)/dtimesub c Sans GW c c deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k)) c c GW formule 1 c c deltatw(k) = deltatw(k)+dtimesub* c $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k)) c c GW formule 2 IF (dtimesub*Tgw(k).lt.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)) ENDIF 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) ENDDO C 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) ENDDO c Determine Ptop from buoyancy integral c---------------------------------------------------------------------- c-1/ Pressure of the level where dth changes sign. Ptop_provis=ph(1) DO k= 2,klev IF (dth(k) .GT. -delta_t_min .and. $ dth(k-1).LT. -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 ENDIF ENDDO 65 CONTINUE c-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 .le. 0) GO TO 70 z = z+dz sum_dth = sum_dth + dth(k)*dz dthmin = min(dthmin,dth(k)) ENDDO 70 CONTINUE c-3/ height of triangle with area= sum_dth and base = dthmin hw = 2.*sum_dth/min(dthmin,-0.5) hw = max(hwmin,hw) c-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 .le. 0) GO TO 75 z = z+dz Ptop = Ph(k)-rho(k)*rg*dz ktop = k ENDDO 75 CONTINUE c-5/Correct ktop and ptop Ptop_new=ptop DO k= ktop,2,-1 IF (dth(k) .GT. -delta_t_min .and. $ dth(k-1).LT. -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 ENDIF ENDDO 275 CONTINUE ptop = ptop_new DO k=klev,1,-1 IF (ph(k+1) .LT. ptop) ktop=k ENDDO c-6/ Set deltatw & deltaqw to 0 above kupper DO k = kupper,klev deltatw(k) = 0. deltaqw(k) = 0. ENDDO c------------------------------------------------------------------ ENDDO ! end sub-timestep loop C ----------------------------------------------------------------- c Get back to tendencies per second DO k = 1,kupper-1 dtls(k) = dtls(k)/dtime dqls(k) = dqls(k)/dtime d_deltatw2(k)=d_deltatw2(k)/dtime d_deltaqw2(k)=d_deltaqw2(k)/dtime d_deltat_gw(k) = d_deltat_gw(k)/dtime ENDDO C 2.1 - Undisturbed area and Wake integrals C --------------------------------------------------------- z = 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. C Potential temperatures and humidity c---------------------------------------------------------- DO k =1,klev rho(k) = p(k)/(rd*te(k)) IF(k .eq. 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) ENDIF 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) ENDDO C Integrals (and wake top level number) C ----------------------------------------------------------- C Initialize sum_thvu to 1st level virt. pot. temp. z = 1. dz = 1. sum_thvu = thu(1)*(1.+eps*qu(1))*dz sum_dth = 0. DO k = 1,klev dz = -(max(ph(k+1),Ptop)-Ph(k))/(rho(k)*rg) IF (dz .LE. 0) GO TO 51 z = z+dz sum_thu = sum_thu + thu(k)*dz sum_tu = sum_tu + tu(k)*dz sum_qu = sum_qu + qu(k)*dz sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz sum_dth = sum_dth + dth(k)*dz sum_dq = sum_dq + deltaqw(k)*dz sum_rho = sum_rho + rhow(k)*dz sum_dtdwn = sum_dtdwn + dtdwn(k)*dz sum_dqdwn = sum_dqdwn + dqdwn(k)*dz ENDDO 51 CONTINUE hw0 = z C 2.1 - WAPE and mean forcing computation C------------------------------------------------------------- C Means av_thu = sum_thu/hw0 av_tu = sum_tu/hw0 av_qu = sum_qu/hw0 av_thvu = sum_thvu/hw0 av_dth = sum_dth/hw0 av_dq = sum_dq/hw0 av_rho = sum_rho/hw0 av_dtdwn = sum_dtdwn/hw0 av_dqdwn = sum_dqdwn/hw0 wape2 = - rg*hw0*(av_dth $ + eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq ))/av_thvu C 2.2 Prognostic variable update C ------------------------------------------------------------ C Filter out bad wakes IF ( wape2 .LT. 0.) THEN 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. ENDDO ELSE print*,'wape2>0' Cstar2 = stark*sqrt(2.*wape2) ENDIF ktopw = ktop IF (ktopw .gt. 0) then Cjyg1 Utilisation d'un h_efficace constant ( ~ feeding layer) ccc heff = 600. C Utilisation de la hauteur hw cc heff = 0.7*hw heff = hw FIP = 0.5*rho(ktopw)*Cstar2**3*heff*2*sqrt(sigmaw*wdens*3.14) FIP = alpk * FIP Cjyg2 ELSE FIP = 0. ENDIF C Limitation de sigmaw c C sécurité : si le wake occuppe plus de 90 % de la surface de la maille, C alors il disparait en se mélangeant à la partie undisturbed IF ((sigmaw.GT.0.9).or. . ((wape.ge.wape2).and.(wape2.le.1.0))) THEN c IF (sigmaw.GT.0.9) THEN DO k = 1,klev dtls(k) = 0. dqls(k) = 0. deltatw(k) = 0. deltaqw(k) = 0. ENDDO wape = 0. hw = hwmin sigmaw = sigmad fip = 0. ENDIF RETURN END