Changeset 5744 for LMDZ6/trunk/libf/phylmd

Timestamp:

Jul 1, 2025, 5:07:30 PM (4 days ago)

Author:

jyg

Message:

Correction of a bug found by Catherine Rio in pbl_surface_mod.F90:
modif = use of smallestreal in order to prevent division by 0.

File:

• LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 (modified) (5 diffs)

LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90

@@ -76 +76 @@
   !$OMP THREADPRIVATE(ok_bug_zg_wk_pbl)
 
+
+!JYG<
+  REAL, SAVE, PROTECTED     :: smallestreal
+  !$OMP THREADPRIVATE(smallestreal)
+
 !FC
 !  integer, save :: iflag_frein
@@ -116 +121 @@
     CHARACTER(len=80)             :: abort_message
     CHARACTER(len = 20)           :: modname = 'pbl_surface_init'
+
+!****************************************************************************************
+! Initialize some module variables
+!****************************************************************************************    
+    smallestreal = tiny(smallestreal)
     
 !****************************************************************************************
@@ -3852 +3862 @@
      ! 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 +3877 @@
      ! 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 +3891 @@
      ! 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))