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

Timestamp:

Aug 16, 2011, 11:06:45 AM (13 years ago)

Author:

acolaitis

Message:

--- AC 03/08/2011 ---
M 267 libf/phymars/physiq.F
<> Minor modification to pass Ch from vdifc to meso_inc_les

M 267 libf/phymars/surflayer_interpol.F
<> Major modification to the formulation of integrals

Now stable for most cases. Some cases with highly negative Monin Obukhov length
remain to be explored.

M 267 libf/phymars/vdif_cd.F
<> Added gustiness to the Richardson computation. Gustiness factor is for now of beta=1., after

several comparisons with LES aerodynamic conductances. May be subject to a minor change (+/- 0.1)
in the near future. (almost transparent for the user)

M 267 libf/phymars/vdifc.F
<> Minor modifications relative to variables.

M 267 libf/phymars/calltherm_mars.F90
<> Added a comment on a "sensitive" parameter that should not be changed without knowing the consequence !

M 267 libf/phymars/meso_inc/meso_inc_les.F
<> Changed the definition for HFX computation in the LES (to be discussed with Aymeric). New definition yields

very similar results to old one and follows a strict definition of what HFX should be.

---

File:

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

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

@@ -5 +5 @@
 !
 !   Subject: interpolation of temperature and velocity norm in the surface layer
-!   by recomputation of surface layer quantities (Richardson, Prandtl, u*, teta*)
+!   by recomputation of surface layer quantities (Richardson, Prandtl, Reynolds, u*, teta*)
 !   -------  
 !
@@ -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(OUT) :: Teta_out(ngrid),u_out(ngrid) ! interpolated fields at z_out : potential temperature and norm(uv)
+      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*
 
@@ -79 +79 @@
 ! For output :
 
-      REAL zu2                 ! Large-scale wind at first layer
+      REAL zu2(ngrid)                  ! Large-scale wind at first layer
       REAL L_mo(ngrid)                ! Monin-Obukhov length
 !-----------------------------------------------------------------------
@@ -88 +88 @@
       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) :
@@ -141 +129 @@
 ! 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) + wmax(ig)**2)
+!     &  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
+
           rib(ig) = (pg/ph(ig,1))
-     &      *pz(ig,1)*pz0(ig)/sqrt(pz(ig,1)*pz0t)
-!     &      *sqrt(pz(ig,1)*pz0(ig))
+!     &      *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))
+     &      *(ph(ig,1)-pts(ig))/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) + wmax(ig)**2)
 
          else
@@ -208 +208 @@
       ENDDO
 
-! Large-scale wind at first layer (with gustiness) and
+! Large-scale wind at first layer (without gustiness) and
 ! u* theta* computation
       DO ig=1,ngrid
@@ -216 +216 @@
            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)))
+
+!           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) then
+             print*, fh(ig),rib(ig),(ph(ig,1)-pts(ig))
+     &               ,log(pz(ig,1)/pz0tcomp(ig)),sqrt(fm(ig))
+           endif 
          endif
       ENDDO
@@ -233 +243 @@
       ENDDO
 
+      DO ig=1,ngrid
       IF(z_out .ge. pz(ig,1)) THEN
            z_out=1.
@@ -252 +263 @@
      &                  and computation is resumed'
       ENDIF
-
+      ENDDO
 
       print*, 'interpolation of u and teta at z_out=',z_out
@@ -265 +276 @@
      &                            *(z_out-pz0tcomp(ig))
         ELSEIF (L_mo(ig) .lt. 0.) THEN
-           u_out(ig)=(ustar(ig)/karman)*((
+
+          IF(L_mo(ig) .gt. -1000.) 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)))
+     & -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.*pz0(ig)/L_mo(ig))**0.25)
+     & -   log(1.+sqrt(1.-16.*z_out/L_mo(ig)))
+     & +   log(1.+sqrt(1.-16.*pz0(ig)/L_mo(ig)))
+     & +   log(z_out/pz0(ig))
      &                                  )
-           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)))))
+
+           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))
      &                                                        )
+
+          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))*(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))
+     &                                   )
+
+           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))
+     &                                                           )
+
+          ENDIF
         ELSE
            u_out(ig)=0.
-           Teta_out(ig)=0.
+           Teta_out(ig)=pts(ig)
         ENDIF
       ENDDO
+
+! Usefull diagnostics for the interpolation routine :
 
          call WRITEDIAGFI(ngrid,'L',