Index: trunk/LMDZ.MARS/README =================================================================== --- trunk/LMDZ.MARS/README (revision 266) +++ trunk/LMDZ.MARS/README (revision 267) @@ -847,2 +847,5 @@ - vdifc.F Now takes into account sub-grid gustiness, evaluated from thermals activity (it's proxy being the maximum vertical velocity) - Minor modification in meso_inc_les.F: u* is now taken from the new vdifc and not recomputed from a simple law +== 08/08/2011 == AC +- Formula correction for the Reynolds in vdif_cd : will affect the thermal roughness length. Results obtained with the precedent formula (previous commit) probably underestimate by about 1K the temperature in the mixing layer. +- Added a new subroutine to libf/phymars called surflayer_interpol.F. This routine can be called as a diagnostic for interpolation of horizontal velocity norm and temperature in the surface layer, in 1D or 3D. The output height is controlled by a parameter z_out, to be set in physiq.F. This might be later added to a .def file. z_out is set to 0. by default, resulting in no call to this subroutine. Index: trunk/LMDZ.MARS/libf/phymars/physiq.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/physiq.F (revision 266) +++ trunk/LMDZ.MARS/libf/phymars/physiq.F (revision 267) @@ -318,9 +318,9 @@ REAL time_phys -c Variables from thermal +c Variables for PBL REAL lmax_th_out(ngridmx),zmax_th(ngridmx) REAL, SAVE :: wmax_th(ngridmx) - REAL, SAVE :: hfmax_th(ngridmx) !! AS: pourquoi en save ? + REAL hfmax_th(ngridmx) REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx) REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx) @@ -328,4 +328,8 @@ REAL dtke_th(ngridmx,nlayermx+1) REAL zcdv(ngridmx) + REAL Teta_out(ngridmx),u_out(ngridmx) ! Interpolated teta and u at z_out + REAL z_out ! height of interpolation between z0 and z1 + REAL ustar(ngridmx),tstar(ngridmx) ! friction velocity and friction potential temp + REAL zu2(ngridmx) c======================================================================= @@ -1553,4 +1557,47 @@ c ---------------------------------------------------------- +c ---------------------------------------------------------- +c PBL OUTPUS +c ---------------------------------------------------------- +c ---------------------------------------------------------- + + +c ---------------------------------------------------------- +c Outputs of surface layer +c ---------------------------------------------------------- + + + z_out=0. + if (calltherm .and. (z_out .gt. 0.)) then + call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th + & ,tsurf,zt(:,:)*(zplay(:,:)/zplev(:,:))**rcp + & ,z_out,Teta_out,u_out,ustar,tstar) + + zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1)) + call WRITEDIAGFI(ngridmx,'sqrt(zu2)', + & 'horizontal velocity norm','m/s', + & 2,zu2) + + call WRITEDIAGFI(ngridmx,'Teta_out', + & 'potential temperature at z_out','K', + & 2,Teta_out) + call WRITEDIAGFI(ngridmx,'u_out', + & 'horizontal velocity norm at z_out','m/s', + & 2,u_out) + call WRITEDIAGFI(ngridmx,'u*', + & 'friction velocity','m/s', + & 2,ustar) + call WRITEDIAGFI(ngridmx,'teta*', + & 'friction potential temperature','K', + & 2,tstar) + else + if((.not. calltherm).and.(z_out .gt. 0.)) then + print*, 'WARNING : no interpolation in surface-layer :' + print*, 'Outputing surface-layer quantities without thermals + & does not make sense' + endif + endif + +c ---------------------------------------------------------- c Outputs of thermals c ---------------------------------------------------------- @@ -1589,4 +1636,11 @@ c ---------------------------------------------------------- +c ---------------------------------------------------------- +c END OF PBL OUTPUS +c ---------------------------------------------------------- +c ---------------------------------------------------------- + + +c ---------------------------------------------------------- c Output in netcdf file "diagsoil.nc" for subterranean c variables (output every "ecritphy", as for writediagfi) Index: trunk/LMDZ.MARS/libf/phymars/surflayer_interpol.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/surflayer_interpol.F (revision 267) +++ trunk/LMDZ.MARS/libf/phymars/surflayer_interpol.F (revision 267) @@ -0,0 +1,299 @@ + SUBROUTINE surflayer_interpol(ngrid,nlay,pz0, + & pg,pz,pu,pv,wmax,pts,ph,z_out,Teta_out,u_out,ustar,tstar) + IMPLICIT NONE +!======================================================================= +! +! Subject: interpolation of temperature and velocity norm in the surface layer +! by recomputation of surface layer quantities (Richardson, Prandtl, u*, teta*) +! ------- +! +! Author: Arnaud Colaitis 05/08/11 +! ------- +! +! Arguments: +! ---------- +! +! inputs: +! ------ +! ngrid size of the horizontal grid +! pg gravity (m s -2) +! pz(ngrid,nlay) height of layers +! pu(ngrid,nlay) u component of the wind +! pv(ngrid,nlay) v component of the wind +! pts(ngrid) surface temperature +! ph(ngrid) potential temperature T*(p/ps)^kappa +! +! +!======================================================================= +! +!----------------------------------------------------------------------- +! Declarations: +! ------------- + +#include "comcstfi.h" + +! Arguments: +! ---------- + + INTEGER, INTENT(IN) :: ngrid,nlay + REAL, INTENT(IN) :: pz0(ngrid) + REAL, INTENT(IN) :: pg,pz(ngrid,nlay) + REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay) + REAL, INTENT(IN) :: wmax(ngrid) + REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay) + REAL, INTENT(INOUT) :: z_out ! output height (in m above surface) + REAL, INTENT(OUT) :: Teta_out(ngrid),u_out(ngrid) ! interpolated fields at z_out : potential temperature and norm(uv) + REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid) ! u* and teta* + +! Local: +! ------ + + INTEGER ig + + REAL karman,nu + DATA karman,nu/.41,0.001/ + SAVE karman,nu + +! Local(2): +! --------- + + REAL rib(ngrid) ! Bulk Richardson number + REAL fm(ngrid) ! stability function for momentum + REAL fh(ngrid) ! stability function for heat + REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T + +! phim = 1+betam*zeta or (1-bm*zeta)**am +! phih = alphah + betah*zeta or alphah(1.-bh*zeta)**ah + REAL betam, betah, alphah, bm, bh, lambda +! ah and am are assumed to be -0.25 and -0.5 respectively + + REAL cdn(ngrid),chn(ngrid) ! neutral momentum and heat drag coefficient + REAL pz0t ! initial thermal roughness length. (local) + REAL ric ! critical richardson number + REAL reynolds(ngrid) ! reynolds number for UBL + REAL prandtl(ngrid) ! prandtl number for UBL + REAL pz0tcomp(ngrid) ! computed z0t + REAL ite + REAL residual + REAL pcdv(ngrid),pcdh(ngrid) +! For output : + + REAL zu2 ! Large-scale wind at first layer + REAL L_mo(ngrid) ! Monin-Obukhov length +!----------------------------------------------------------------------- +! couche de surface: +! ------------------ + + tstar(:)=0. + ustar(:)=0. + reynolds(:)=0. + +! Original formulation : + +! DO ig=1,ngrid +! z1=1.E+0 + pz(ig,1)/pz0(ig) +! zcd0=karman/log(z1) +! zcd0=zcd0*zcd0 +! pcdv(ig)=zcd0 +! pcdh(ig)=zcd0 +! ENDDO + +! print*,'old : cd,ch; ',pcdv,pcdh + +! New formulation (AC) : + +! phim = 1+betam*zeta or (1-bm*zeta)**am +! phih = alphah + betah*zeta or alphah(1.-bh*zeta)**ah +! am=-0.25, ah=-0.5 + + pz0t = 0. ! for the sake of simplicity + pz0tcomp(:) = 0.1*pz0(:) + rib(:)=0. + pcdv(:)=0. + pcdh(:)=0. + +! this formulation assumes alphah=1., implying betah=betam +! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers : + bm=16. !UBL + bh=16. !UBL + alphah=1. + betam=5. !SBL + betah=5. !SBL + lambda=(sqrt(bh/bm))/alphah + ric=betah/(betam**2) + + DO ig=1,ngrid + + ite=0. + residual=abs(pz0tcomp(ig)-pz0t) + + do while((residual .gt. 0.01*pz0(ig)) .and. (ite .lt. 10.)) + + pz0t=pz0tcomp(ig) + + if ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then + +! Classical Richardson number formulation + +! rib(ig) = (pg/ph(ig,1))*pz(ig,1)*(ph(ig,1)-pts(ig)) +! & /(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)) + +! Richardson number formulation proposed by D.E. England et al. (1995) + + rib(ig) = (pg/ph(ig,1)) + & *pz(ig,1)*pz0(ig)/sqrt(pz(ig,1)*pz0t) +! & *sqrt(pz(ig,1)*pz0(ig)) + & *(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t))) + & *(ph(ig,1)-pts(ig)) + & /(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)) + + else + print*,'warning, infinite Richardson at surface' + print*,pu(ig,1),pv(ig,1) + rib(ig) = ric ! traiter ce cas ! + endif + + z1z0=pz(ig,1)/pz0(ig) + z1z0t=pz(ig,1)/pz0t + + cdn(ig)=karman/log(z1z0) + cdn(ig)=cdn(ig)*cdn(ig) + chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t) + +! Stable case : + if (rib(ig) .gt. 0.) then +! From D.E. England et al. (95) + prandtl(ig)=1. + if(rib(ig) .lt. ric) then +! Assuming alphah=1. and bh=bm for stable conditions : + fm(ig)=((ric-rib(ig))/ric)**2 + fh(ig)=((ric-rib(ig))/ric)**2 + else +! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface + fm(ig)=0. + fh(ig)=0. + endif +! Unstable case : + else +! From D.E. England et al. (95) + fm(ig)=sqrt(1.-lambda*bm*rib(ig)) + fh(ig)=(1./alphah)*((1.-lambda*bh*rib(ig))**0.5)* + & (1.-lambda*bm*rib(ig))**0.25 + prandtl(ig)=alphah*((1.-lambda*bm*rib(ig))**0.25)/ + & ((1.-lambda*bh*rib(ig))**0.5) + endif + + reynolds(ig)=karman*sqrt(fm(ig))*sqrt(pu(ig,1)**2 + pv(ig,1)**2) + & *pz0(ig)/(log(z1z0)*nu) + pz0tcomp(ig)=pz0(ig)*exp(-karman*7.3* + & (reynolds(ig)**0.25)*(prandtl(ig)**0.5)) + + + residual = abs(pz0t-pz0tcomp(ig)) + ite = ite+1 +! if(ite .eq. 10) then +! print*, 'iteration max reached' +! endif +! print*, "iteration nnumber, residual",ite,residual + + enddo ! of while + + pz0t=0. + +! Drag computation : + + pcdv(ig)=cdn(ig)*fm(ig) + pcdh(ig)=chn(ig)*fh(ig) + + ENDDO + +! Large-scale wind at first layer (with gustiness) and +! u* theta* computation + DO ig=1,ngrid + + if (rib(ig) .gt. ric) then + ustar(ig)=0. + tstar(ig)=0. + else + zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)+(0.*wmax(ig))**2 + ustar(ig)=karman*sqrt(fm(ig)*zu2)/(log(pz(ig,1)/pz0(ig))) + tstar(ig)=karman*fh(ig)*(ph(ig,1)-pts(ig)) + & /(log(pz(ig,1)/pz0tcomp(ig))*sqrt(fm(ig))) + endif + ENDDO + +! Monin Obukhov length : + + DO ig=1,ngrid + if (rib(ig) .gt. ric) then + L_mo(ig)=0. + else + L_mo(ig)=pts(ig)*(ustar(ig)**2)/(pg*karman*tstar(ig)) !as defined here, L is positive for SBL, negative for UBL + endif + ENDDO + + IF(z_out .ge. pz(ig,1)) THEN + z_out=1. + print*, 'Warning, z_out should be less than the first + & vertical grid-point' + print*, 'z1 =',pz(ig,1) + print*, 'z_out=',z_out + print*, 'z_out has been changed to 1m + & and computation is resumed' + ENDIF + + IF(z_out .lt. pz0(ig)) THEN + z_out=1. + print*, 'Warning, z_out should be more than the roughness + & length' + print*, 'z0 =',pz0(ig) + print*, 'z_out=',z_out + print*, 'z_out has been changed to z0 + & and computation is resumed' + ENDIF + + + print*, 'interpolation of u and teta at z_out=',z_out + + DO ig=1,ngrid + IF (L_mo(ig) .gt. 0.) THEN + u_out(ig)=(ustar(ig)/karman)*log(z_out/pz0(ig)) + + & 5.*(ustar(ig)/(karman*L_mo(ig)))*(z_out-pz0(ig)) + Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman)) + & *log(z_out/pz0tcomp(ig)) + + & 5.*(tstar(ig)/(prandtl(ig)*karman*L_mo(ig))) + & *(z_out-pz0tcomp(ig)) + ELSEIF (L_mo(ig) .lt. 0.) THEN + u_out(ig)=(ustar(ig)/karman)*(( + & 2.*atan((1.-16.*z_out/L_mo(ig))**0.25) + & +log((-1.+(1.-16.*z_out/L_mo(ig))**0.25) + & /(-1.+(1.-16.*z_out/L_mo(ig))**0.25)))-( + & 2.*atan((1.-16.*pz0(ig)/L_mo(ig))**0.25) + & +log((-1.+(1.-16.*pz0(ig)/L_mo(ig))**0.25) + & /(-1.+(1.-16.*pz0(ig)/L_mo(ig))**0.25))) + & ) + Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*(( + & log((-1.+sqrt(1.-16.*z_out/L_mo(ig))) + & /(1.+sqrt(1.-16.*z_out/L_mo(ig)))))-( + & log((-1.+sqrt(1.-16.*pz0tcomp(ig)/L_mo(ig))) + & /(1.+sqrt(1.-16.*pz0tcomp(ig)/L_mo(ig))))) + & ) + ELSE + u_out(ig)=0. + Teta_out(ig)=0. + ENDIF + ENDDO + + call WRITEDIAGFI(ngrid,'L', + & 'Monin Obukhov length','m', + & 2,L_mo) + call WRITEDIAGFI(ngrid,'z0T', + & 'thermal roughness length','m', + & 2,pz0tcomp) + call WRITEDIAGFI(ngrid,'z0', + & 'roughness length','m', + & 2,pz0) + + + RETURN + END Index: trunk/LMDZ.MARS/libf/phymars/vdif_cd.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/vdif_cd.F (revision 266) +++ trunk/LMDZ.MARS/libf/phymars/vdif_cd.F (revision 267) @@ -52,9 +52,9 @@ INTEGER ig - REAL karman + REAL karman,nu ! Von Karman constant and fluid kinematic viscosity LOGICAL firstcal - DATA karman/.41/ + DATA karman,nu/.41,0.001/ DATA firstcal/.true./ - SAVE karman + SAVE karman,nu c Local(2): @@ -179,6 +179,6 @@ endif - reynolds(ig)=sqrt(fm(ig))*sqrt(pu(ig,1)**2 + pv(ig,1)**2)*pz0(ig) - & /(log(z1z0)) + reynolds(ig)=karman*sqrt(fm(ig))*sqrt(pu(ig,1)**2 + pv(ig,1)**2) + & *pz0(ig)/(log(z1z0)*nu) pz0tcomp(ig)=pz0(ig)*exp(-karman*7.3* & (reynolds(ig)**0.25)*(prandtl(ig)**0.5)) Index: trunk/LMDZ.MARS/libf/phymars/vdifc.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/vdifc.F (revision 266) +++ trunk/LMDZ.MARS/libf/phymars/vdifc.F (revision 267) @@ -238,6 +238,6 @@ zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) - & +(1.2*wmax(ig))**2 - + & +(1.2*wmax(ig))**2 !subgrid gustiness is used to enhance surface flux only, and not u*,t* computations + zcdv(ig)=zcdv_true(ig)*sqrt(zu2) zcdh(ig)=zcdh_true(ig)*sqrt(zu2)