Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 1376)
+++ trunk/LMDZ.MARS/README	(revision 1377)
@@ -2150,2 +2150,5 @@
 == 20/10/2014 == JYC
 - Update and bug correction on the escape fluxes of H and H2.
+
+== 13/02/2015 == EM
+- Update pbl_parameter with corrected version by Anne-Flore Moreau.
Index: trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F	(revision 1376)
+++ trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F	(revision 1377)
@@ -145,6 +145,5 @@
        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
@@ -158,5 +157,5 @@
      &        *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)
+     &        *(ph(ig,1)-pts(ig))/(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1))
          ELSE
             print*,'warning, infinite Richardson at surface'
@@ -220,6 +219,5 @@
 
 
-! u* theta* computation
-! and Monin Obukhov length:
+! u* theta* computation:
 
       DO ig=1,ngrid
@@ -227,5 +225,5 @@
            ustar(ig)=0.
            tstar(ig)=0.
-           L_mo(ig)=0.
+
          ELSE
            ustar(ig)=sqrt(pcdv(ig))
@@ -233,19 +231,7 @@
            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:
+
+         ENDIF
+      ENDDO
 
       DO ig=1,ngrid
@@ -253,124 +239,17 @@
            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)
+         ELSE
+           u_out(ig,n)= ustar(ig)*log(zout/pz0(ig))/
+     &(karman*sqrt(fm(ig)))
+
+           Teta_out(ig,n)=pts(ig)+(tstar(ig)*sqrt(fm(ig))*log(zout/
+     & (pz0tcomp(ig)))/
+     &(karman*fh(ig)))
+
         ENDIF
-        IF(zout .lt. pz0(ig)) THEN
+
+        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
+        ENDIF 
 
       ENDDO
@@ -488,8 +367,8 @@
      &   'momentum stab function','m/s',2,fh)
             call WRITEDIAGFI(ngrid,'B_pbl',
-     &   'flottabilité','m/',2,(ph(:,1)-pts(:))/pts(:))
+     &   'buoyancy','m/',2,(ph(:,1)-pts(:))/pts(:))
             call WRITEDIAGFI(ngrid,'zot_pbl',
-     &   'flottabilité','ms',2,pz0tcomp)
-       call WRITEDIAGFI(ngrid,'zz1','flottabilité','m',2,zzlay(:,1))
+     &   'buoyancy','ms',2,pz0tcomp)
+       call WRITEDIAGFI(ngrid,'zz1','buoyancy','m',2,zzlay(:,1))
 #endif
 
Index: trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 1376)
+++ trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 1377)
@@ -518,5 +518,7 @@
       zdtsurf(:)=0
       dqsurf(:,:)=0
+#ifdef DUSTSTORM
       pq_tmp(:,:,:)=0
+#endif
       igout=ngrid/2+1 
 
@@ -1567,6 +1569,6 @@
             call WRITEDIAGFI(ngrid,'teta_star',
      &   'friction potential temperature','K',dimout,tstar)
-            call WRITEDIAGFI(ngrid,'L',
-     &   'Monin Obukhov length','m',dimout,L_mo)
+!            call WRITEDIAGFI(ngrid,'L',
+!     &   'Monin Obukhov length','m',dimout,L_mo)
             call WRITEDIAGFI(ngrid,'vvv',
      &   'Vertical velocity variance at zout','m',dimout,vvv)
