Index: LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 5743)
+++ LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 5744)
@@ -76,4 +76,9 @@
   !$OMP THREADPRIVATE(ok_bug_zg_wk_pbl)
 
+
+!JYG<
+  REAL, SAVE, PROTECTED     :: smallestreal
+  !$OMP THREADPRIVATE(smallestreal)
+
 !FC
 !  integer, save :: iflag_frein
@@ -116,4 +121,9 @@
     CHARACTER(len=80)             :: abort_message
     CHARACTER(len = 20)           :: modname = 'pbl_surface_init'
+
+!****************************************************************************************
+! Initialize some module variables
+!****************************************************************************************    
+    smallestreal = tiny(smallestreal)
     
 !****************************************************************************************
@@ -3852,6 +3862,6 @@
      ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
           ustar(i,nsrf)=yustar(j)
-          u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)**2+vmer(j)**2)
-          v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)**2+vmer(j)**2)
+          u10m(i,nsrf)=(yu10m(j) * uzon(j))/max(SQRT(uzon(j)**2+vmer(j)**2), smallestreal)
+          v10m(i,nsrf)=(yu10m(j) * vmer(j))/max(SQRT(uzon(j)**2+vmer(j)**2), smallestreal)
 !
           DO k = 1, 6
@@ -3867,6 +3877,6 @@
      ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
           ustar_x(i,nsrf)=yustar_x(j)
-          u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/SQRT(uzon_x(j)**2+vmer_x(j)**2)
-          v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/SQRT(uzon_x(j)**2+vmer_x(j)**2)
+          u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/max(SQRT(uzon_x(j)**2+vmer_x(j)**2), smallestreal)
+          v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/max(SQRT(uzon_x(j)**2+vmer_x(j)**2), smallestreal)
 !
           DO k = 1, 6
@@ -3881,6 +3891,6 @@
      ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
           ustar_w(i,nsrf)=yustar_w(j)
-          u10m_w(i,nsrf)=(yu10m_w(j) * uzon_w(j))/SQRT(uzon_w(j)**2+vmer_w(j)**2)
-          v10m_w(i,nsrf)=(yu10m_w(j) * vmer_w(j))/SQRT(uzon_w(j)**2+vmer_w(j)**2)
+          u10m_w(i,nsrf)=(yu10m_w(j) * uzon_w(j))/max(SQRT(uzon_w(j)**2+vmer_w(j)**2), smallestreal)
+          v10m_w(i,nsrf)=(yu10m_w(j) * vmer_w(j))/max(SQRT(uzon_w(j)**2+vmer_w(j)**2), smallestreal)
 !
           ustar(i,nsrf) = ustar_x(i,nsrf) + wake_s(i)*(ustar_w(i,nsrf)-ustar_x(i,nsrf))