Index: trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F (revision 528) +++ trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F (revision 529) @@ -1,5 +1,5 @@ - SUBROUTINE pbl_parameters(ngrid,nlay,pz0, - & pg,pz,pu,pv,wmax,hfmax,zmax,pts,ph,z_out, - & Teta_out,u_out,ustar,tstar,wstar,L_mo,vhf,vvv) + SUBROUTINE pbl_parameters(ngrid,nlay,ps,pplay,pz0, + & pg,zzlay,zzlev,pu,pv,wstar_in,hfmax,zmax,pts,ph,z_out, + & Teta_out,u_out,ustar,tstar,L_mo,vhf,vvv) IMPLICIT NONE !======================================================================= @@ -21,9 +21,9 @@ ! pz0(ngrid) surface roughness length ! pg gravity (m s -2) -! pz(ngrid,nlay) height of layers +! zzlay(ngrid,nlay) height of mid-layers +! zzlev(ngrid,nlay+1) height of layers interfaces ! pu(ngrid,nlay) u component of the wind ! pv(ngrid,nlay) v component of the wind -! wmax(ngrid) maximum vertical velocity in thermals (might not be needed -! if the computation of w* works) +! wstar_in(ngrid) free convection velocity in thermals ! hfmax(ngrid) maximum vertical turbulent heat flux in thermals ! zmax(ngrid) height reached by the thermals (pbl height) @@ -56,12 +56,18 @@ INTEGER, INTENT(IN) :: ngrid,nlay - REAL, INTENT(IN) :: pz0(ngrid) - REAL, INTENT(IN) :: pg,pz(ngrid,nlay) + REAL, INTENT(IN) :: pz0(ngrid),ps(ngrid),pplay(ngrid,nlay) + REAL, INTENT(IN) :: pg,zzlay(ngrid,nlay),zzlev(ngrid,nlay) REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay) - REAL, INTENT(IN) :: wmax(ngrid),hfmax(ngrid),zmax(ngrid) + REAL, INTENT(IN) :: wstar_in(ngrid),hfmax(ngrid),zmax(ngrid) REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay) REAL, INTENT(IN) :: z_out + +! Outputs: +! -------- + REAL, INTENT(OUT) :: Teta_out(ngrid),u_out(ngrid) - REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid),wstar(ngrid) + REAL T_out(ngrid) + REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid) + REAL wstar(ngrid) REAL, INTENT(OUT) :: L_mo(ngrid) @@ -139,5 +145,6 @@ ite=0. residual=abs(pz0tcomp(ig)-pz0t) - zu2(ig)=MAX(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1),wmax(ig)**2) + zu2(ig)=MAX(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) + & ,(0.3*wstar_in(ig))**2) DO WHILE((residual .gt. 0.01*pz0(ig)) .and. (ite .lt. 10.)) @@ -147,6 +154,6 @@ ! Richardson number formulation proposed by D.E. England et al. (1995) rib(ig) = (pg/ph(ig,1)) - & *sqrt(pz(ig,1)*pz0(ig)) - & *(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t))) + & *sqrt(zzlev(ig,2)*pz0(ig)) + & *(((log(zzlev(ig,2)/pz0(ig)))**2)/(log(zzlev(ig,2)/pz0t))) & *(ph(ig,1)-pts(ig))/zu2(ig) ELSE @@ -156,6 +163,6 @@ ENDIF - z1z0=pz(ig,1)/pz0(ig) - z1z0t=pz(ig,1)/pz0t + z1z0=zzlev(ig,2)/pz0(ig) + z1z0t=zzlev(ig,2)/pz0t cdn(ig)=karman/log(z1z0) @@ -310,4 +317,10 @@ ENDIF + T_out(:) = Teta_out(:)*(exp( + & (zout/zzlay(:,1))*(log(pplay(:,1)/ps)) + & ) + & )**rcp + + !------------------------------------------------------------------------ !------------------------------------------------------------------------ @@ -325,6 +338,6 @@ DO k=1,nlay-1 DO ig=1, ngrid - IF (abs(zmax(ig)-pz(ig,k)) .lt. - & abs(zmax(ig)-pz(ig,pbl_height_index(ig)))) THEN + IF (abs(zmax(ig)-zzlay(ig,k)) .lt. + & abs(zmax(ig)-zzlay(ig,pbl_height_index(ig)))) THEN pbl_height_index(ig)=k ENDIF @@ -337,5 +350,5 @@ DO k=1,nlay-1 DO ig=1, ngrid - dteta(ig,k) = (ph(ig,k+1)-ph(ig,k))/(pz(ig,k+1)-pz(ig,k)) + dteta(ig,k) = (ph(ig,k+1)-ph(ig,k))/(zzlay(ig,k+1)-zzlay(ig,k)) ENDDO ENDDO @@ -424,13 +437,13 @@ & 'Monin Obukhov length','m', & dimout,L_mo) - call WRITEDIAGFI(ngrid,'w_star', - & 'Free convection velocity','m', - & dimout,wstar) - call WRITEDIAGFI(ngrid,'z_i', - & 'PBL height','m', - & dimout,zmax) - call WRITEDIAGFI(ngrid,'hf_max', - & 'Maximum vertical heat flux','m', - & dimout,hfmax) +! call WRITEDIAGFI(ngrid,'w_star', +! & 'Free convection velocity','m', +! & dimout,wstar) +! call WRITEDIAGFI(ngrid,'z_i', +! & 'PBL height','m', +! & dimout,zmax) +! call WRITEDIAGFI(ngrid,'hf_max', +! & 'Maximum vertical heat flux','m', +! & dimout,hfmax) call WRITEDIAGFI(ngrid,'vvv', & 'Vertical velocity variance at zout','m', Index: trunk/LMDZ.MARS/libf/phymars/physiq.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/physiq.F (revision 528) +++ trunk/LMDZ.MARS/libf/phymars/physiq.F (revision 529) @@ -777,5 +777,10 @@ $ zdum1,zdum2,zdh,pdq,zflubid, $ zdudif,zdvdif,zdhdif,zdtsdif,q2, - & zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th) + & zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th +#ifdef MESOSCALE + & ,flag_LES +#endif + & ) + #ifdef MESOSCALE @@ -1848,6 +1853,6 @@ if (calltherm .and. (z_out .gt. 0.)) then - call pbl_parameters(ngrid,nlayer,z0, - & g,zzlay,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out, + call pbl_parameters(ngrid,nlayer,ps,zplay,z0, + & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out, & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv) @@ -1938,9 +1943,9 @@ ! THERMALS STUFF 1D - z_out=0. + z_out=1. if (calltherm .and. (z_out .gt. 0.)) then - call pbl_parameters(ngrid,nlayer,z0, - & g,zzlay,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out, + call pbl_parameters(ngrid,nlayer,ps,zplay,z0, + & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out, & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv) Index: trunk/LMDZ.MARS/libf/phymars/vdifc.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/vdifc.F (revision 528) +++ trunk/LMDZ.MARS/libf/phymars/vdifc.F (revision 529) @@ -5,5 +5,10 @@ $ pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf, $ pdudif,pdvdif,pdhdif,pdtsrf,pq2, - $ pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,hfmax) + $ pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,hfmax +#ifdef MESOSCALE + & ,flag_LES +#endif + & ) + USE ioipsl_getincom IMPLICIT NONE @@ -137,4 +142,7 @@ REAL qco2,mmean +#ifdef MESOSCALE + LOGICAL flag_LES !! pour LES avec isfflx!=0 +#endif c ** un petit test de coherence c -------------------------- @@ -242,5 +250,6 @@ DO ilev=1,nlay DO ig=1,ngrid - qco2=pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep + qco2=MAX(1.E-30 + & ,pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep) c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2) mmean=1/(A*qco2 +B) @@ -339,8 +348,7 @@ ENDIF - ! Some usefull diagnostics for the new surface layer parametrization : - -! call WRITEDIAGFI(ngridmx,'pcdv', + +! call WRITEDIAGFI(ngridmx,'pcdv', ! & 'momentum drag','no units', ! & 2,zcdv_true) @@ -361,9 +369,9 @@ c ** schema de diffusion turbulente dans la couche limite c ---------------------------------------------------- -! + CALL vdif_kc(ptimestep,g,pzlev,pzlay & ,pu,pv,ph,zcdv_true & ,pq2,zkv,zkh,zq) -! + ! pt(:,:)=ph(:,:)*ppopsk(:,:) ! CALL yamada4(ngrid,nlay,ptimestep,g,r,pplev,pt @@ -519,5 +527,4 @@ ENDDO - ! calcul de zc et zd pour la couche top en prenant en compte le terme ! de variation de masse (on fait une boucle pour que ça converge) @@ -526,5 +533,7 @@ llnt(:)=1 - +#ifdef MESOSCALE + IF (.not.flag_LES) THEN +#endif DO ig=1,ngrid DO l=1,nlay @@ -537,4 +546,8 @@ 5 continue ENDDO + +#ifdef MESOSCALE + ENDIF +#endif DO ig=1,ngrid @@ -589,5 +602,4 @@ ztsrf2(ig)=z1(ig)/z2(ig) ! pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep !incremented outside loop - zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig) @@ -633,9 +645,14 @@ c Using the wind modified by friction for lifting and sublimation c ---------------------------------------------------------------- -c DO ig=1,ngrid -c zu2(ig)=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1) -c zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig)) -c zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig)) -c ENDDO + +! This is computed above and takes into account surface-atmosphere flux +! enhancement by subgrid gustiness and atmospheric-stability related +! variations of transfer coefficients. + +! DO ig=1,ngrid +! zu2(ig)=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1) +! zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig)) +! zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig)) +! ENDDO c Calcul du flux vertical au bas de la premiere couche (dust) :