Index: trunk/LMDZ.MARS/libf/phymars/surflayer_interpol.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/surflayer_interpol.F (revision 269) +++ trunk/LMDZ.MARS/libf/phymars/surflayer_interpol.F (revision 272) @@ -42,5 +42,5 @@ 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(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* @@ -57,4 +57,5 @@ ! Local(2): ! --------- + REAL zout REAL rib(ngrid) ! Bulk Richardson number @@ -84,5 +85,5 @@ ! couche de surface: ! ------------------ - + zout=z_out tstar(:)=0. ustar(:)=0. @@ -184,5 +185,7 @@ endif - reynolds(ig)=karman*sqrt(fm(ig))*sqrt(pu(ig,1)**2 + pv(ig,1)**2) + reynolds(ig)=karman*sqrt(fm(ig)) + & *sqrt(zu2(ig)) +! & *sqrt(pu(ig,1)**2 + pv(ig,1)**2) & *pz0(ig)/(log(z1z0)*nu) pz0tcomp(ig)=pz0(ig)*exp(-karman*7.3* @@ -244,35 +247,35 @@ DO ig=1,ngrid - IF(z_out .ge. pz(ig,1)) THEN - z_out=1. + IF(zout .ge. pz(ig,1)) THEN + zout=1. print*, 'Warning, z_out should be less than the first - & vertical grid-point' + & 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' + & 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) + 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 z0 - & and computation is resumed' + print*, 'z_out has been changed to z0t + & and computation is resumed' ENDIF ENDDO - print*, 'interpolation of u and teta at z_out=',z_out + 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(z_out/pz0(ig)) + - & 5.*(ustar(ig)/(karman*L_mo(ig)))*(z_out-pz0(ig)) + 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(z_out/pz0tcomp(ig)) + + & *log(zout/pz0tcomp(ig)) + & 5.*(tstar(ig)/(prandtl(ig)*karman*L_mo(ig))) - & *(z_out-pz0tcomp(ig)) + & *(zout-pz0tcomp(ig)) ELSEIF (L_mo(ig) .lt. 0.) THEN @@ -280,17 +283,17 @@ u_out(ig)=(ustar(ig)/karman)*( - & 2.*atan((1.-16.*z_out/L_mo(ig))**0.25) + & 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.*z_out/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.*z_out/L_mo(ig))) + & - log(1.+sqrt(1.-16.*zout/L_mo(ig))) & + log(1.+sqrt(1.-16.*pz0(ig)/L_mo(ig))) - & + log(z_out/pz0(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.*z_out/L_mo(ig))) - & + log(z_out/pz0tcomp(ig)) + & -2.*log(1.+sqrt(1.-16.*zout/L_mo(ig))) + & + log(zout/pz0tcomp(ig)) & ) @@ -305,15 +308,15 @@ u_out(ig)=(ustar(ig)/karman)*( - & (4./L_mo(ig))*(z_out-pz0(ig)) - & + (20./(L_mo(ig))**2)*(z_out**2-pz0(ig)**2) - & + (160./(L_mo(ig))**3)*(z_out**3-pz0(ig)**3) - & + log(z_out/pz0(ig)) + & (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))*(z_out-pz0tcomp(ig)) - & + (48./(L_mo(ig))**2)*(z_out**2-pz0tcomp(ig)**2) - & + (1280./(3.*(L_mo(ig))**3))*(z_out**3-pz0tcomp(ig)**3) - & + log(z_out/pz0tcomp(ig)) + & (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)) & ) @@ -327,13 +330,13 @@ ! Usefull diagnostics for the interpolation routine : - 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) +! 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)