Changeset 272 for trunk/LMDZ.MARS/libf/phymars/surflayer_interpol.F

Timestamp:

Aug 26, 2011, 10:40:39 AM (13 years ago)

Author:

acolaitis

Message:

Modifications du surflayer_interpol.F : needed local duplicate for z_out, and included subgrid gustiness in friction quantities estimation.

File:

• trunk/LMDZ.MARS/libf/phymars/surflayer_interpol.F (modified) (8 diffs)

trunk/LMDZ.MARS/libf/phymars/surflayer_interpol.F

@@ -42 +42 @@
       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 +57 @@
 !    Local(2):
 !    ---------
+      REAL zout
 
       REAL rib(ngrid)  ! Bulk Richardson number
@@ -84 +85 @@
 !   couche de surface:
 !   ------------------
-      
+      zout=z_out 
       tstar(:)=0.
       ustar(:)=0.
@@ -184 +185 @@
       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 +247 @@
 
       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 +283 @@
            
            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 +308 @@
 
            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 +330 @@
 ! 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)