SUBROUTINE pbl_parameters(ngrid,nlay,ps,pplay,pz0, & pg,zzlay,zzlev,pu,pv,wstar_in,hfmax,zmax,pts,ph,z_out,n_out, & T_out,u_out,ustar,tstar,L_mo,vhf,vvv) IMPLICIT NONE !======================================================================= ! ! Anlysis of the PBL from input temperature, wind field and thermals outputs. ! ! ------- ! ! Author: Arnaud Colaitis 09/01/12 ! ------- ! ! Arguments: ! ---------- ! ! inputs: ! ------ ! ngrid size of the horizontal grid ! nlay size of the vertical grid ! pz0(ngrid) surface roughness length ! pg gravity (m s -2) ! 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 ! wstar_in(ngrid) free convection velocity in PBL ! hfmax(ngrid) maximum vertical turbulent heat flux in thermals ! zmax(ngrid) height reached by the thermals (pbl height) ! pts(ngrid) surface temperature ! ph(ngrid,nlay) potential temperature T*(p/ps)^kappa ! z_out(n_out) heights of interpolation ! n_out number of points for interpolation ! ! outputs: ! ------ ! ! Teta_out(ngrid,n_out) interpolated teta ! u_out(ngrid,n_out) interpolated u ! ustar(ngrid) friction velocity ! tstar(ngrid) friction temperature ! L_mo(ngrid) monin_obukhov length ! ! !======================================================================= ! !----------------------------------------------------------------------- ! Declarations: ! ------------- #include "comcstfi.h" #include "callkeys.h" ! Arguments: ! ---------- INTEGER, INTENT(IN) :: ngrid,nlay,n_out 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) :: wstar_in(ngrid),hfmax(ngrid),zmax(ngrid) REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay) REAL, INTENT(IN) :: z_out(n_out) ! Outputs: ! -------- REAL, INTENT(OUT) :: T_out(ngrid,n_out),u_out(ngrid,n_out) REAL Teta_out(ngrid,n_out) REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid) REAL, INTENT(OUT) :: L_mo(ngrid) ! Local: ! ------ INTEGER ig,k,n REAL karman,nu DATA karman,nu/.41,0.001/ SAVE karman,nu ! Local(2): ! --------- REAL zout 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,zcd0,z1 REAL pcdv(ngrid),pcdh(ngrid) REAL zu2(ngrid) ! Large-scale wind at first layer REAL pbl_teta(ngrid) ! mixed-layer potential temperature INTEGER pbl_height_index(ngrid) ! index of nearest vertical grid point for zmax REAL dteta(ngrid,nlay),x(ngrid) ! potential temperature gradient and z/zi REAL dvhf(ngrid),dvvv(ngrid) ! dimensionless vertical heat flux and ! dimensionless vertical velocity variance REAL vhf(ngrid),vvv(ngrid) ! vertical heat flux and vertical velocity variance INTEGER ii(1) ! temporary INTEGER dimout !------------------------------------------------------------------------ !------------------------------------------------------------------------ ! PART I : RICHARDSON/REYNOLDS/THERMAL_ROUGHNESS/STABILITY_COEFFICIENTS !------------------------------------------------------------------------ !------------------------------------------------------------------------ DO n=1,n_out c Initialisation : L_mo(:)=0. ustar(:)=0. tstar(:)=0. zout=z_out(n) reynolds(:)=0. pz0t = 0. 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) ! zu2(ig)=MAX(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) ! & ,(0.3*wstar_in(ig))**2) zu2(ig)=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) & + (log(1.+0.7*wstar_in(ig) + 2.3*wstar_in(ig)**2))**2 DO WHILE((residual .gt. 0.01*pz0(ig)) .and. (ite .lt. 10.)) pz0t=pz0tcomp(ig) IF (zu2(ig) .ne. 0.) THEN ! Richardson number formulation proposed by D.E. England et al. (1995) rib(ig) = (pg/pts(ig)) & *sqrt(zzlay(ig,1)*pz0(ig)) & *(((log(zzlay(ig,1)/pz0(ig)))**2)/(log(zzlay(ig,1)/pz0t))) & *(ph(ig,1)-pts(ig))/zu2(ig) ELSE print*,'warning, infinite Richardson at surface' print*,pu(ig,1),pv(ig,1) rib(ig) = ric ENDIF z1z0=zzlay(ig,1)/pz0(ig) z1z0t=zzlay(ig,1)/pz0t cdn(ig)=karman/log(z1z0) cdn(ig)=cdn(ig)*cdn(ig) chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t) ! STABLE BOUNDARY LAYER : 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 BOUNDARY LAYER : 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(zu2(ig)) & *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 ENDDO ! of while pz0t=0. ! Drag computation: pcdv(ig)=cdn(ig)*fm(ig) pcdh(ig)=chn(ig)*fh(ig) ENDDO ! of ngrid !------------------------------------------------------------------------ !------------------------------------------------------------------------ ! PART II : USTAR/TSTAR/U_OUT/TETA_OUT COMPUTATIONS !------------------------------------------------------------------------ !------------------------------------------------------------------------ ! u* theta* computation ! and Monin Obukhov length: DO ig=1,ngrid IF (rib(ig) .ge. ric) THEN ustar(ig)=0. tstar(ig)=0. L_mo(ig)=0. ELSE ustar(ig)=sqrt(pcdv(ig)) & *sqrt(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)) tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1)) & /sqrt(pcdv(ig)) 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 ! Monin Obukhov length: ! DO ig=1,ngrid ! IF (rib(ig) .ge. 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 ! Interpolation: DO ig=1,ngrid IF(zout .lt. pz0tcomp(ig)) THEN u_out(ig,n)=0. Teta_out(ig,n)=pts(ig) ELSEIF (L_mo(ig) .gt. 0.) THEN ! u_out(ig,n)=(ustar(ig)/karman)*log(zout/pz0(ig)) + ! & 5.*(ustar(ig)/(karman*L_mo(ig)))*(zout-pz0(ig)) ! Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman)) ! & *log(zout/pz0tcomp(ig)) + ! & 5.*(tstar(ig)/(prandtl(ig)*karman*L_mo(ig))) ! & *(zout-pz0tcomp(ig)) IF ((zout/L_mo(ig).lt.ric).and.(pz0(ig).lt.ric)) THEN ! & ((zout/L_mo(ig).gt.ric).and.(pz0(ig).gt.ric)) ) THEN u_out(ig,n)=(ustar(ig)/karman)*( & log((ric-pz0(ig)/L_mo(ig))/(ric-zout/L_mo(ig))) & +log(zout/pz0(ig)) & ) ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0(ig).gt.ric)) THEN u_out(ig,n)=(ustar(ig)/karman)*( & log((zout-ric*L_mo(ig))/(pz0(ig)-ric*L_mo(ig))) & +log(pz0(ig)/zout) & ) ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0(ig).lt.ric)) THEN !integral of the stability function does not converge u_out(ig,n)=sqrt(pu(ig,1)**2 + pv(ig,1)**2) ENDIF IF((zout/L_mo(ig).lt.ric).and.(pz0tcomp(ig).lt.ric)) THEN ! & ((zout/L_mo(ig).gt.ric).and.(pz0tcomp(ig).gt.ric)) ) THEN Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*( & log((ric-pz0tcomp(ig)/L_mo(ig))/(ric-zout/L_mo(ig))) & +log(zout/pz0tcomp(ig)) & ) ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0tcomp(ig).gt.ric)) THEN Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*( & log((zout-ric*L_mo(ig))/(pz0tcomp(ig)-ric*L_mo(ig))) & +log(pz0tcomp(ig)/zout) & ) ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0tcomp(ig).lt.ric)) THEN !integral of the stability function does not converge Teta_out(ig,n)=ph(ig,1) ENDIF ELSEIF (L_mo(ig) .lt. 0.) THEN IF(L_mo(ig) .gt. -1000.) THEN u_out(ig,n)=(ustar(ig)/karman)*( & 2.*atan((1.-16.*zout/L_mo(ig))**0.25) & -2.*atan((1.-16.*pz0(ig)/L_mo(ig))**0.25) & -2.*log(1.+(1.-16.*zout/L_mo(ig))**0.25) & +2.*log(1.+(1.-16.*pz0(ig)/L_mo(ig))**0.25) & - log(1.+sqrt(1.-16.*zout/L_mo(ig))) & + log(1.+sqrt(1.-16.*pz0(ig)/L_mo(ig))) & + log(zout/pz0(ig)) & ) Teta_out(ig,n)=MAX(pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*( & 2.*log(1.+sqrt(1.-16.*pz0tcomp(ig)/L_mo(ig))) & -2.*log(1.+sqrt(1.-16.*zout/L_mo(ig))) & + log(zout/pz0tcomp(ig)) & ),ph(ig,1)) ELSE ! We have to treat the case where L is very large and negative, ! corresponding to a very nearly stable atmosphere (but not quite !) ! i.e. teta* <0 and almost zero. We choose the cutoff ! corresponding to L_mo < -1000 and use a 3rd order taylor expansion. The difference ! between the two expression for z/L = -1/1000 is -1.54324*10^-9 ! (we do that to avoid using r*4 precision, otherwise, we get -inf values) u_out(ig,n)=(ustar(ig)/karman)*( & (4./L_mo(ig))*(zout-pz0(ig)) & + (20./(L_mo(ig))**2)*(zout**2-pz0(ig)**2) & + (160./(L_mo(ig))**3)*(zout**3-pz0(ig)**3) & + log(zout/pz0(ig)) & ) Teta_out(ig,n)=MAX(pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*( & (8./L_mo(ig))*(zout-pz0tcomp(ig)) & + (48./(L_mo(ig))**2)*(zout**2-pz0tcomp(ig)**2) & + (1280./(3.*(L_mo(ig))**3))*(zout**3-pz0tcomp(ig)**3) & + log(zout/pz0tcomp(ig)) & ),ph(ig,1)) ENDIF ELSE u_out(ig,n)=sqrt(pu(ig,1)**2 + pv(ig,1)**2) Teta_out(ig,n)=ph(ig,1) ENDIF IF(zout .lt. pz0(ig)) THEN u_out(ig,n)=0. ENDIF ! Final check IF (L_mo(ig) .gt. 0) THEN IF (Teta_out(ig,n) .gt. ph(ig,1)) THEN Teta_out(ig,n)=ph(ig,1) ELSEIF (Teta_out(ig,n) .lt. pts(ig)) THEN Teta_out(ig,n)=pts(ig) ENDIF ELSEIF (L_mo(ig) .lt. 0) THEN IF (Teta_out(ig,n) .lt. ph(ig,1)) THEN Teta_out(ig,n)=ph(ig,1) ELSEIF (Teta_out(ig,n) .gt. pts(ig)) THEN Teta_out(ig,n)=pts(ig) ENDIF ENDIF IF (u_out(ig,n) .gt. sqrt(pu(ig,1)**2 + pv(ig,1)**2)) THEN u_out(ig,n)=sqrt(pu(ig,1)**2 + pv(ig,1)**2) ENDIF ENDDO ! when using convective adjustment without thermals, a vertical potential temperature ! profile is assumed up to the thermal roughness length. Hence, theoretically, theta ! interpolated at any height in the surface layer is theta at the first level. IF ((.not.calltherm).and.(calladj)) THEN Teta_out(:,n)=ph(:,1) u_out(:,n)=(sqrt(cdn(:))*sqrt(pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1)) & /karman)*log(zout/pz0(:)) ENDIF T_out(:,n) = Teta_out(:,n)*(exp( & (zout/zzlay(:,1))*(log(pplay(:,1)/ps)) & ) & )**rcp ENDDO !of n=1,n_out !------------------------------------------------------------------------ !------------------------------------------------------------------------ ! PART III : WSTAR COMPUTATION !------------------------------------------------------------------------ !------------------------------------------------------------------------ ! Detection of the mixed-layer potential temperature ! ------------ ! Nearest index for the pbl height IF (calltherm) THEN pbl_height_index(:)=1 DO k=1,nlay-1 DO ig=1, ngrid IF (abs(zmax(ig)-zzlay(ig,k)) .lt. & abs(zmax(ig)-zzlay(ig,pbl_height_index(ig)))) THEN pbl_height_index(ig)=k ENDIF ENDDO ENDDO ! Potential temperature gradient dteta(:,nlay)=0. DO k=1,nlay-1 DO ig=1, ngrid dteta(ig,k) = (ph(ig,k+1)-ph(ig,k))/(zzlay(ig,k+1)-zzlay(ig,k)) ENDDO ENDDO ! Computation of the pbl mixed layer temperature DO ig=1, ngrid ii=MINLOC(abs(dteta(ig,1:pbl_height_index(ig)))) pbl_teta(ig) = ph(ig,ii(1)) ENDDO !------------------------------------------------------------------------ !------------------------------------------------------------------------ ! PART IV : VERTICAL_VELOCITY_VARIANCE/VERTICAL_TURBULENT_FLUX PROFILES !------------------------------------------------------------------------ !------------------------------------------------------------------------ ! We follow Spiga et. al 2010 (QJRMS) ! ------------ DO ig=1, ngrid IF (zmax(ig) .gt. 0.) THEN x(ig) = zout/zmax(ig) ELSE x(ig) = 999. ENDIF ENDDO DO ig=1, ngrid ! dimensionless vertical heat flux IF (x(ig) .le. 0.3) THEN dvhf(ig) = ((-3.85/log(x(ig)))+0.07*log(x(ig))) & *exp(-4.61*x(ig)) ELSEIF (x(ig) .le. 1.) THEN dvhf(ig) = -1.52*x(ig) + 1.24 ELSE dvhf(ig) = 0. ENDIF ! dimensionless vertical velocity variance IF (x(ig) .le. 1.) THEN dvvv(ig) = 2.05*(x(ig)**(2./3.))*(1.-0.64*x(ig))**2 ELSE dvvv(ig) = 0. ENDIF ENDDO vhf(:) = dvhf(:)*hfmax(:) vvv(:) = dvvv(:)*(wstar_in(:))**2 ENDIF ! of if calltherm call WRITEDIAGFI(ngrid,'rib_pbl', & 'Richardson in pbl parameter','m/s',2,rib) call WRITEDIAGFI(ngrid,'cdn_pbl', & 'neutral momentum coef','m/s',2,cdn) call WRITEDIAGFI(ngrid,'fm_pbl', & 'momentum stab function','m/s',2,fm) call WRITEDIAGFI(ngrid,'uv', & 'wind norm first layer','m/s',2,sqrt(zu2)) call WRITEDIAGFI(ngrid,'uvtrue', & 'wind norm first layer','m/s',2,sqrt(pu(:,1)**2+pv(:,1)**2)) call WRITEDIAGFI(ngrid,'chn_pbl', & 'neutral momentum coef','m/s',2,chn) call WRITEDIAGFI(ngrid,'fh_pbl', & 'momentum stab function','m/s',2,fh) call WRITEDIAGFI(ngrid,'B_pbl', & 'flottabilité','m/',2,(ph(:,1)-pts(:))/pts(:)) call WRITEDIAGFI(ngrid,'zot_pbl', & 'flottabilité','ms',2,pz0tcomp) call WRITEDIAGFI(ngrid,'zz1','flottabilité','m',2,zzlay(:,1)) RETURN END