SUBROUTINE surflayer_interpol(ngrid,nlay,pz0, & pg,pz,pu,pv,wmax,pts,ph,z_out,Teta_out,u_out,ustar,tstar,L_mo) IMPLICIT NONE !======================================================================= ! ! Subject: interpolation of temperature and velocity norm in the surface layer ! by recomputation of surface layer quantities (Richardson, Prandtl, Reynolds, 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" #include "callkeys.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(IN) :: 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* REAL, INTENT(OUT) :: L_mo(ngrid) ! Monin-Obukhov length ! Local: ! ------ INTEGER ig 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) ! For output : REAL zu2(ngrid) ! Large-scale wind at first layer !----------------------------------------------------------------------- ! couche de surface: ! ------------------ c Init : L_mo(:)=0. ustar(:)=0. tstar(:)=0. zout=z_out reynolds(:)=0. ! 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) ! IF((pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) .lt. 1.) ! & .and. (wmax(ig) .gt. 0.)) THEN zu2(ig)= ! & (MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),wmax(ig)**2)) & ( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)) ! & pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) ! ELSE ! zu2(ig)=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) ! ENDIF if(.not.callrichsl) then rib(ig) = (pg/ph(ig,1)) & *sqrt(pz(ig,1)*pz0(ig)) & *(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t))) & *(ph(ig,1)-pts(ig))/(zu2(ig)+6.) else 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))/(zu2(ig) + (0.5*wmax(ig))**2) ! & /(MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),wmax(ig)**2)) ! & /( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + wmax(ig)**2) endif 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(zu2(ig)+(0.5*wmax(ig))**2) ! & *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 ! endif !of if callrichsl ! Large-scale wind at first layer (without gustiness) and ! u* theta* computation DO ig=1,ngrid if (rib(ig) .ge. ric) then ustar(ig)=0. tstar(ig)=0. else ! ustar(ig)=karman*sqrt(fm(ig)*zu2(ig))/(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))) !simpler definition of u* and teta*, equivalent to the formula above : ustar(ig)=sqrt(pcdv(ig))*sqrt(zu2(ig)) tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1)) & /sqrt(pcdv(ig)) if ((tstar(ig) .lt. -50) .and. callrichsl) then print*, fh(ig),rib(ig),(ph(ig,1)-pts(ig)) & ,log(pz(ig,1)/pz0tcomp(ig)),sqrt(fm(ig)) endif 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 DO ig=1,ngrid IF(zout .ge. pz(ig,1)) THEN zout=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(zout .lt. pz0tcomp(ig)) THEN zout=pz0tcomp(ig) print*, 'Warning, z_out should be more than the thermal & roughness length' print*, 'z0 =',pz0tcomp(ig) print*, 'z_out=',z_out print*, 'z_out has been changed to z0t & and computation is resumed' ENDIF ENDDO print*, 'interpolation of u and teta at z_out=',zout DO ig=1,ngrid IF (L_mo(ig) .gt. 0.) THEN u_out(ig)=(ustar(ig)/karman)*log(zout/pz0(ig)) + & 5.*(ustar(ig)/(karman*L_mo(ig)))*(zout-pz0(ig)) Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman)) & *log(zout/pz0tcomp(ig)) + & 5.*(tstar(ig)/(prandtl(ig)*karman*L_mo(ig))) & *(zout-pz0tcomp(ig)) ELSEIF (L_mo(ig) .lt. 0.) THEN IF(L_mo(ig) .gt. -1000.) THEN u_out(ig)=(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)=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)) & ) 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)=(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)=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)) & ) ENDIF ELSE u_out(ig)=0. Teta_out(ig)=pts(ig) 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(:)=ph(:,1) ENDIF ! Usefull diagnostics for the interpolation routine : ! ! call WRITEDIAGFI(ngrid,'Ri', ! & 'Richardson','m', ! & 2,rib) ! ! call WRITEDIAGFI(ngrid,'z0T', ! & 'thermal roughness length','m', ! & 2,pz0tcomp) ! call WRITEDIAGFI(ngrid,'z0', ! & 'roughness length','m', ! & 2,pz0) RETURN END