! ! $Id: wake.F 1669 2012-10-16 12:41:50Z 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 * !*************************************************************** !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 wdens : densite de poches !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 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 pi : (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_ref: initial number of wakes per unit area (3D) or per !C unit length (2D), at the beginning of each time step !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" #include "YOMCST.h" #include "cvthermo.h" #include "iniprint.h" !c Arguments en entree !c-------------------- 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 !c Sorties !c-------- 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 !c Variables internes !c------------------- !c 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 !c Variables de sauvegarde REAL, dimension(klon,klev) :: deltatw0 REAL, dimension(klon,klev) :: deltaqw0 REAL, dimension(klon,klev) :: te, qe REAL, dimension(klon) :: sigmaw0, sigmaw1 !c Variables pour les GW REAL, DIMENSION(klon) :: LL REAL, dimension(klon,klev) :: N2 REAL, dimension(klon,klev) :: Cgw REAL, dimension(klon,klev) :: Tgw !c 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 !c 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/ !c 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 !ccc nrlmd real, dimension(klon) :: death_rate,nat_rate real, dimension(klon,klev) :: entr real, dimension(klon,klev) :: detr !C------------------------------------------------------------------------- !c Initialisations !c------------------------------------------------------------------------- !c print*, 'wake initialisations' !c Essais d'initialisation avec sigmaw = 0.02 et hw = 10. !c------------------------------------------------------------------------- DATA wapecut,sigmad, hwmin /5.,.02,10./ !ccc nrlmd DATA sigmaw_max /0.4/ DATA dens_rate /0.1/ !ccc !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------------------------------------------------------------------------- !ccc nrlmd coefgw=10 !c coefgw=1 !c wdens0 = 1.0/(alon**2) !ccc nrlmd wdens = 1.0/(alon**2) !ccc nrlmd stark = 0.50 !cCRtest !ccc nrlmd alpk=0.1 !c alpk = 1.0 !c alpk = 0.5 !c alpk = 0.05 !c stark = 0.33 Alpk = 0.25 wdens_ref = 8.e-12 coefgw = 4. Crep_upper=0.9 Crep_sol=1.0 !ccc 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 !c !c Initialisation de toutes des densites a wdens_ref. !c Les densites peuvent evoluer si les poches debordent !c (voir au tout debut de la boucle sur les substeps) wdens = wdens_ref !c !c print*,'stark',stark !c print*,'alpk',alpk !c print*,'wdens',wdens !c print*,'coefgw',coefgw !ccc !C Minimum value for |T_wake - T_undist|. Used for wake top definition !c------------------------------------------------------------------------- delta_t_min = 0.2 !C 1. - Save initial values and initialize tendencies !C -------------------------------------------------- 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 ENDDO ENDDO !c sigmaw1=sigmaw !c IF (sigd_con.GT.sigmaw1) THEN !c print*, 'sigmaw,sigd_con', sigmaw, sigd_con !c ENDIF DO i=1, klon !cc 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 ENDDO !C !C !C 2. - Prognostic part !C -------------------- !C !C !C 2.1 - Undisturbed area and Wake integrals !C --------------------------------------------------------- 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. ENDDO !c !c Distance between wakes DO i = 1,klon LL(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i)) ENDDO !C Potential temperatures and humidity !c---------------------------------------------------------- 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 .eq. 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) ENDIF ! 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) ENDDO ENDDO DO k = 1, klev-1 DO i=1, klon IF(k.eq.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))) ENDIF 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) ENDDO ENDDO DO i=1, klon N2(i,klev)=0 ZH(i,klev)=0 Cgw(i,klev)=0 Tgw(i,klev)=0 ENDDO !c Calcul de la masse volumique moyenne de la colonne (bdlmd) !c----------------------------------------------------------------- DO k=1,klev DO i=1, klon epaisseur1(i,k)=0. epaisseur2(i,k)=0. ENDDO ENDDO 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) ENDDO 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) ENDDO ENDDO !C !C Choose an integration bound well above wake top !c----------------------------------------------------------------- !c !C Pupper = 50000. ! melting level !c Pupper = 60000. !c Pupper = 80000. ! essais pour case_e DO i = 1,klon Pupper(i) = 0.6*ph(i,1) Pupper(i) = max(Pupper(i), 45000.) !ccc Pupper(i) = 60000. ENDDO !C !C Determine Wake top pressure (Ptop) from buoyancy integral !C -------------------------------------------------------- !c !c-1/ Pressure of the level where dth becomes less than delta_t_min. DO i=1,klon ptop_provis(i)=ph(i,1) ENDDO DO k= 2,klev DO i=1,klon !c !IM v3JYG; ptop_provis(i).LT. ph(i,1) !c IF (dth(i,k) .GT. -delta_t_min .and. & & dth(i,k-1).LT. -delta_t_min .and. & & ptop_provis(i).EQ. 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)) ENDIF ENDDO ENDDO !c-2/ dth integral DO i=1,klon sum_dth(i) = 0. dthmin(i) = -delta_t_min z(i) = 0. ENDDO 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) .gt. 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)) ENDIF ENDDO ENDDO !c-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)) ENDDO !c-4/ now, get Ptop DO i=1,klon z(i) = 0. ptop(i) = ph(i,1) ENDDO 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) .gt. 0) THEN z(i) = z(i)+dz(i) ptop(i) = ph(i,k)-rho(i,k)*rg*dz(i) ENDIF ENDDO ENDDO !C-5/ Determination de ktop et kupper DO k=klev,1,-1 DO i=1,klon IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k IF (ph(i,k+1) .lt. pupper(i)) kupper(i)=k ENDDO ENDDO !c On evite kupper = 1 et kupper = klev DO i=1,klon kupper(i) = max(kupper(i),2) kupper(i) = min(kupper(i),klev-1) ENDDO !c-6/ Correct ktop and ptop DO i = 1,klon ptop_new(i)=ptop(i) ENDDO DO k= klev,2,-1 DO i=1,klon IF (k .LE. ktop(i) .and. & & ptop_new(i) .EQ. ptop(i) .and. & & dth(i,k) .GT. -delta_t_min .and. & & dth(i,k-1).LT. -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)) ENDIF ENDDO ENDDO DO i=1,klon ptop(i) = ptop_new(i) ENDDO DO k=klev,1,-1 DO i=1,klon IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k ENDDO ENDDO !c !c-5/ Set deltatw & deltaqw to 0 above kupper !c DO k = 1,klev DO i=1,klon IF (k.GE. kupper(i)) THEN deltatw(i,k) = 0. deltaqw(i,k) = 0. ENDIF ENDDO ENDDO !c !C !C Vertical gradient of LS omega !C DO k = 1,klev DO i=1,klon IF (k.LE. kupper(i)) THEN dp_omgb(i,k) = (omgb(i,k+1) - omgb(i,k))/(ph(i,k+1)-ph(i,k)) ENDIF ENDDO ENDDO !C !C Integrals (and wake top level number) !C -------------------------------------- !C !C 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. ENDDO 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) .GT. 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) ENDIF ENDDO ENDDO !c DO i=1,klon hw0(i) = z(i) ENDDO !c !C !C 2.1 - WAPE and mean forcing computation !C --------------------------------------- !C !C --------------------------------------- !C !C 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) !c 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) ENDDO !C !C 2.2 Prognostic variable update !C ------------------------------ !C !C Filter out bad wakes DO k = 1,klev DO i=1,klon IF ( wape(i) .LT. 0.) THEN deltatw(i,k) = 0. deltaqw(i,k) = 0. dth(i,k) = 0. ENDIF ENDDO ENDDO !c DO i=1,klon IF ( wape(i) .LT. 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. ENDIF ENDDO !c !c Check qx and qw positivity !c -------------------------- 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)) ) ENDDO 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).le.q0_min(i)) THEN q0_min(i)=q1_min(i) ENDIF ENDDO ENDDO !c DO i = 1,klon OK_qx_qw(i) = q0_min(i) .GE. 0. alpha(i) = 1. ENDDO !c !CC ----------------------------------------------------------------- !C Sub-time-stepping !C ----------------- !C nsub=10 dtimesub=dtime/nsub !c !c------------------------------------------------------------ DO isubstep = 1,nsub !c------------------------------------------------------------ !c !c 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) .GE. 1. ENDDO !c !ccc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement négatif de ktop à kupper -------- !ccc On calcule pour cela une densité wdens0 pour laquelle on aurait un entrainement nul --- DO i=1,klon !cc print *,' isubstep,wk_adv(i),cstar(i),wape(i) ', !cc $ isubstep,wk_adv(i),cstar(i),wape(i) IF (wk_adv(i) .AND. cstar(i).GT.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) .LE. wdens0*1.1 ) THEN wdens(i) = wdens0 ENDIF !cc print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i) !cc $ ,ph(i,1)-pupper(i)', !cc $ omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i) !cc $ ,ph(i,1)-pupper(i) ENDIF ENDDO !ccc 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) ENDIF ENDDO DO i=1,klon IF (wk_adv(i)) THEN !ccc nrlmd Introduction du taux de mortalité des poches et test sur sigmaw_max=0.4 !ccc d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub IF (sigmaw(i).ge.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 !c $ - nat_rate(i)*sigmaw(i)*dtimesub !cc print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), !cc $ death_rate(i),ktop(i),kupper(i)', !cc $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), !cc $ death_rate(i),ktop(i),kupper(i) !c sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub !c sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!! !c wdens = wdens0/(10.*sigmaw) !c sigmaw =max(sigmaw,sigd_con) !c sigmaw =max(sigmaw,sigmad) ENDIF ENDDO !C !C !c 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 ENDDO ENDDO DO k= 1,klev+1 DO i = 1,klon if (wk_adv(i)) THEN !!! nrlmd omg(i,k)=0. end if ENDDO ENDDO !c 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))) ENDIF ENDDO !c DO k= 2,klev DO i = 1,klon IF( wk_adv(i) .AND. k .LE. 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) ENDIF ENDDO ENDDO !c 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) ENDIF ENDDO !c !c ----------------- !c From m/s to Pa/s !c ----------------- !c 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)) ENDIF ENDDO !c DO k= 1,klev DO i = 1,klon IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN omg(i,k) = - rho(i,k)*rg*omg(i,k) dp_deltomg(i,k) = dp_deltomg(i,1) ENDIF ENDDO ENDDO !c !c raccordement lineaire de omg de ptop a pupper DO i=1,klon IF ( wk_adv(i) .AND. kupper(i) .GT. 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)) ENDIF ENDDO !c !cc DO i=1,klon !cc print*,'Pente entre 0 et kupper (référence)' !cc $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1)) !cc print*,'Pente entre ktop et kupper' !cc $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i)) !cc ENDDO !cc DO k= 1,klev DO i = 1,klon IF( wk_adv(i) .AND. k .GT. ktop(i) .AND. k .LE. 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)) ENDIF ENDDO ENDDO !ccc nrlmd !cc DO i=1,klon !cc print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1) !cc END DO !ccc !c !c !c-- Compute wake average vertical velocity omgbw !c !c 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) ENDIF ENDDO ENDDO !c-- and its vertical gradient dp_omgbw !c 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)) ENDIF ENDDO ENDDO !C !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 DO i=1,klon IF ( wk_adv(i)) THEN alpha_up(i,k) = 0. IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1. ENDIF ENDDO ENDDO !c 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) ENDIF ENDDO RRd1 = -1. RRd2 = 1. !c !c-- Get [Th1,Th2], dth and [q1,q2] !c DO k= 1,klev DO i = 1,klon IF( wk_adv(i) .AND. k .LE. 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 ENDIF ENDDO ENDDO 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 ENDDO DO k= 2,klev DO i = 1,klon IF( wk_adv(i) .AND. k .LE. 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) ENDIF ENDDO ENDDO DO i=1,klon IF( wk_adv(i)) THEN omgbdth(i,1) = 0. omgbdq(i,1) = 0. ENDIF ENDDO DO k= 2,klev DO i = 1,klon IF( wk_adv(i) .AND. k .LE. 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)) ENDIF ENDDO ENDDO !c !c----------------------------------------------------------------- DO k= 1,klev DO i = 1,klon IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN !c----------------------------------------------------------------- !c !c Compute redistribution (advective) term !c 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) !c print*,'d_deltatw=',d_deltatw(i,k) !c 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)) !c print*,'d_deltaqw=',d_deltaqw(i,k) !c !c and increment large scale tendencies !c !c !C !CC ----------------------------------------------------------------- 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)) & !ccc 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)) & !ccc & & )*ppi(i,k) !c 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)) & !ccc 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)) & !ccc & & ) !ccc nrlmd ELSE IF(wk_adv(i) .AND. k .EQ. 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))) & & ) ENDIF !ccc ENDDO ENDDO !c------------------------------------------------------------------ !C !C Increment state variables DO k= 1,klev DO i = 1,klon !ccc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN !ccc !c !c 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))) !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(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))) !c print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k) !c 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))) !c print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k) !c !ccc nrlmd Prise en compte du taux de mortalité !ccc 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) !ccc spread(i,k) = (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/ !ccc $ sigmaw(i) !c 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 !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(i,k).lt.1.e-10) THEN d_deltatw(i,k) = dtimesub* & & (ff(i)+dtKE(i,k)+dtPBL(i,k) & !ccc $ -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)) & !ccc & & -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) & !ccc $ -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)) & !ccc & & -Tgw(i,k)*deltatw(i,k)) ENDIF 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) & !ccc $ -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))) !ccc !ccc nrlmd !ccc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k) !ccc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k) !ccc ENDIF ENDDO ENDDO !C !C Scale tendencies so that water vapour remains positive in w and x. !C call wake_vec_modulation(klon,klev,wk_adv,epsilon,qe,d_qe,deltaqw, & & d_deltaqw,sigmaw,d_sigmaw,alpha) !c !ccc nrlmd !cc print*,'alpha' !cc do i=1,klon !cc print*,alpha(i) !cc end do !ccc DO k = 1,klev DO i = 1,klon IF( wk_adv(i) .AND. k .LE. 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) ENDIF ENDDO ENDDO DO i = 1,klon IF( wk_adv(i)) THEN d_sigmaw(i)=alpha(i)*d_sigmaw(i) ENDIF ENDDO !C 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 .LE. kupper(i)) THEN dtls(i,k)=dtls(i,k)+d_te(i,k) dqls(i,k)=dqls(i,k)+d_qe(i,k) !ccc 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) !ccc ENDIF ENDDO ENDDO DO k= 1,klev DO i = 1,klon IF( wk_adv(i) .AND. k .LE. 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) !cc print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k) !cc $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k) ENDIF ENDDO ENDDO DO i = 1,klon IF( wk_adv(i)) THEN sigmaw(i) = sigmaw(i)+d_sigmaw(i) ENDIF ENDDO !c !C !c Determine Ptop from buoyancy integral !c --------------------------------------- !c !c- 1/ Pressure of the level where dth changes sign. !c DO i=1,klon IF ( wk_adv(i)) THEN Ptop_provis(i)=ph(i,1) ENDIF ENDDO !c DO k= 2,klev DO i=1,klon IF ( wk_adv(i) .AND. & & Ptop_provis(i) .EQ. ph(i,1) .AND. & & dth(i,k) .GT. -delta_t_min .and. & & dth(i,k-1).LT. -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)) ENDIF ENDDO ENDDO !c !c- 2/ dth integral !c DO i=1,klon if (wk_adv(i)) then !!! nrlmd sum_dth(i) = 0. dthmin(i) = -delta_t_min z(i) = 0. end if ENDDO 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) .gt. 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)) ENDIF ENDIF ENDDO ENDDO !c !c- 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)) ENDIF ENDDO !c !c- 4/ now, get Ptop !c DO i=1,klon if (wk_adv(i)) then !!! nrlmd ktop(i) = 0 z(i)=0. end if ENDDO !c 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) .gt. 0) THEN z(i) = z(i)+dz(i) Ptop(i) = Ph(i,k)-rho(i,k)*rg*dz(i) ktop(i) = k ENDIF ENDIF ENDDO ENDDO !c !c 4.5/Correct ktop and ptop !c DO i=1,klon IF ( wk_adv(i)) THEN Ptop_new(i)=ptop(i) ENDIF ENDDO !c DO k= klev,2,-1 DO i=1,klon !IM v3JYG; IF (k .GE. ktop(i) IF ( wk_adv(i) .AND. & & k .LE. ktop(i) .AND. & & ptop_new(i) .EQ. ptop(i) .AND. & & dth(i,k) .GT. -delta_t_min .and. & & dth(i,k-1).LT. -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)) ENDIF ENDDO ENDDO !c !c DO i=1,klon IF ( wk_adv(i)) THEN ptop(i) = ptop_new(i) ENDIF ENDDO DO k=klev,1,-1 DO i=1,klon if (wk_adv(i)) then !!! nrlmd IF (ph(i,k+1) .LT. ptop(i)) ktop(i)=k end if ENDDO ENDDO !c !c 5/ Set deltatw & deltaqw to 0 above kupper !c DO k = 1,klev DO i=1,klon IF ( wk_adv(i) .AND. k .GE. kupper(i)) THEN deltatw(i,k) = 0. deltaqw(i,k) = 0. ENDIF ENDDO ENDDO !c !C !c-------------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 ENDDO !C !C Integrals (and wake top level number) !C -------------------------------------- !C !C 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 ENDDO 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) .GT. 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) ENDIF end if ENDDO ENDDO !c DO i=1,klon if (wk_adv(i)) then !!! nrlmd hw0(i) = z(i) end if ENDDO !c !C !C - WAPE and mean forcing computation !C --------------------------------------- !C !C --------------------------------------- !C !C 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) !c 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 ENDDO !C !C Filter out bad wakes DO k = 1,klev DO i=1,klon if (wk_adv(i)) then !!! nrlmd IF ( wape(i) .LT. 0.) THEN deltatw(i,k) = 0. deltaqw(i,k) = 0. dth(i,k) = 0. ENDIF end if ENDDO ENDDO !c DO i=1,klon if (wk_adv(i)) then !!! nrlmd IF ( wape(i) .LT. 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. ENDIF end if ENDDO ENDDO ! end sub-timestep loop !C !C ----------------------------------------------------------------- !c Get back to tendencies per second !c DO k = 1,klev DO i=1,klon !ccc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN IF ( OK_qx_qw(i) .AND. k .LE. kupper(i)) THEN !ccc 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 !cc print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k) !cc $ ,death_rate(i)*sigmaw(i) ENDIF ENDDO ENDDO !c !c---------------------------------------------------------- !c Determine wake final state; recompute wape, cstar, ktop; !c filter out bad wakes. !c---------------------------------------------------------- !c !C 2.1 - Undisturbed area and Wake integrals !C --------------------------------------------------------- DO i=1,klon !ccc nrlmd if (wk_adv(i)) then !!! nrlmd if (OK_qx_qw(i)) then !ccc 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 ENDDO !C Potential temperatures and humidity !c---------------------------------------------------------- DO k =1,klev DO i=1,klon !ccc nrlmd IF ( wk_adv(i)) THEN if (OK_qx_qw(i)) then !ccc rho(i,k) = p(i,k)/(rd*te(i,k)) IF(k .eq. 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) ENDIF 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) ENDIF ENDDO ENDDO !C Integrals (and wake top level number) !C ----------------------------------------------------------- !C Initialize sum_thvu to 1st level virt. pot. temp. DO i=1,klon !ccc nrlmd IF ( wk_adv(i)) THEN if (OK_qx_qw(i)) then !ccc z(i) = 1. dz(i) = 1. sum_thvu(i) = thu(i,1)*(1.+eps*qu(i,1))*dz(i) sum_dth(i) = 0. ENDIF ENDDO DO k = 1,klev DO i=1,klon !ccc nrlmd IF ( wk_adv(i)) THEN if (OK_qx_qw(i)) then !ccc dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) IF (dz(i) .GT. 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) ENDIF ENDIF ENDDO ENDDO !c DO i=1,klon !ccc nrlmd IF ( wk_adv(i)) THEN if (OK_qx_qw(i)) then !ccc hw0(i) = z(i) ENDIF ENDDO !c !C - WAPE and mean forcing computation !C------------------------------------------------------------- !C Means DO i=1, klon !ccc nrlmd IF ( wk_adv(i)) THEN if (OK_qx_qw(i)) then !ccc 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) ENDIF ENDDO !C Prognostic variable update !C ------------------------------------------------------------ !C Filter out bad wakes !c DO k = 1,klev DO i=1,klon !ccc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN if (OK_qx_qw(i) .AND. wape2(i) .LT. 0.) then !ccc deltatw(i,k) = 0. deltaqw(i,k) = 0. dth(i,k) = 0. ENDIF ENDDO ENDDO !c DO i=1, klon !ccc nrlmd IF ( wk_adv(i)) THEN if (OK_qx_qw(i)) then !ccc IF ( wape2(i) .LT. 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.ge.10) print*,'wape2>0' Cstar2(i) = stark*sqrt(2.*wape2(i)) gwake(i) = .TRUE. ENDIF ENDIF ENDDO !c DO i=1, klon !ccc nrlmd IF ( wk_adv(i)) THEN if (OK_qx_qw(i)) then !ccc ktopw(i) = ktop(i) ENDIF ENDDO !c DO i=1, klon !ccc nrlmd IF ( wk_adv(i)) THEN if (OK_qx_qw(i)) then !ccc IF (ktopw(i) .gt. 0 .and. gwake(i)) 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(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) !Cjyg2 ELSE FIP(i) = 0. ENDIF ENDIF ENDDO !c !C Limitation de sigmaw !ccc nrlmd !c DO i=1,klon !c IF (OK_qx_qw(i)) THEN !c IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max !c ENDIF !c ENDDO !ccc DO k = 1,klev DO i=1, klon !ccc nrlmd On maintient désormais constant sigmaw en régime permanent !ccc IF ((sigmaw(i).GT.sigmaw_max).or. IF ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or. & & (ktopw(i).le.2) .OR. & & .not. OK_qx_qw(i) ) THEN !ccc dtls(i,k) = 0. dqls(i,k) = 0. deltatw(i,k) = 0. deltaqw(i,k) = 0. ENDIF ENDDO ENDDO !c !ccc nrlmd On maintient désormais constant sigmaw en régime permanent DO i=1, klon IF ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or. & & (ktopw(i).le.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) ENDIF !cc print*,'wape wape2 ktopw OK_qx_qw =', !cc $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i) ENDDO !c !c RETURN END SUBROUTINE WAKE SUBROUTINE wake_vec_modulation(nlon,nl,wk_adv,epsilon,qe,d_qe, & & deltaqw,d_deltaqw,sigmaw,d_sigmaw,alpha) !c------------------------------------------------------ !cDtermination du coefficient alpha tel que les tendances !c corriges alpha*d_G, pour toutes les grandeurs G, correspondent !c a une humidite positive dans la zone (x) et dans la zone (w). !c------------------------------------------------------ !c !c 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 !c Output REAL alpha(nlon) !c Internal variables REAL zeta(nlon,nl) REAL alpha1(nlon) REAL x,a,b,c,discrim REAL epsilon ! DATA epsilon/1.e-15/ !c DO k=1,nl DO i = 1,nlon IF (wk_adv(i)) THEN IF ((deltaqw(i,k)+d_deltaqw(i,k)).ge.0.) then zeta(i,k)=0. ELSE zeta(i,k)=1. END IF ENDIF ENDDO 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 !c print*, 'x, a, b, c, discrim', x, a, b, c, discrim IF (a+b .GE. 0.) THEN !! Condition suffisante pour la positivité de ovap alpha1(i)=1. ELSE IF (x .GE. 0.) THEN alpha1(i)=1. ELSE IF (a .GT. 0.) THEN alpha1(i)=0.9*min( (2.*c)/(-b+sqrt(discrim)), & & (-b+sqrt(discrim))/(2.*a) ) ELSE IF (a .eq. 0.) then alpha1(i)=0.9*(-c/b) ELSE !c print*,'a,b,c discrim',a,b,c discrim alpha1(i)=0.9*max( (2.*c)/(-b+sqrt(discrim)), & & (-b+sqrt(discrim))/(2.*a) ) ENDIF ENDIF ENDIF alpha(i) = min(alpha(i),alpha1(i)) ENDIF ENDDO ENDDO ! 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 * !*************************************************************** !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 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" #include "iniprint.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) 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 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 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 if(prt_level.ge.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. ENDDO ELSE if(prt_level.ge.10) 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 if(prt_level.ge.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. ENDDO ELSE if(prt_level.ge.10) 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 ! correction NICOLAS . ((wape.ge.wape2).and.(wape2.le.1.0))) THEN IF ((sigmaw.GT.0.9).or. & & ((wape.ge.wape2).and.(wape2.le.1.0)).or.(ktopw.le.2)) THEN !IM cf NR/JYG 251108 . ((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 SUBROUTINE WAKE_scal