New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 11785 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90 – NEMO

Ignore:
Timestamp:
2019-10-24T12:38:43+02:00 (4 years ago)
Author:
laurent
Message:

Fixes to prevent FPE-related errors at runtime (masked regions).

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90

    r11772 r11785  
    195195      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p6@sbcblk_algo_coare3p6' 
    196196      !!---------------------------------------------------------------------------------- 
    197  
    198197      IF ( kt == nit000 ) CALL COARE3P6_INIT(l_use_cs, l_use_wl) 
    199198 
     
    219218         zsst = T_s ! backing up the bulk SST 
    220219         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction 
    221          q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!! 
     220         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 
    222221      END IF 
    223222 
     
    240239 
    241240      z0     = alfa_charn_3p6(U_zu)*u_star*u_star/grav + 0.11_wp*znu_a/u_star 
    242       z0     = MIN(ABS(z0), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
     241      z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
     242 
    243243      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 
    244       z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
     244      z0t    = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
    245245 
    246246      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd 
     
    260260      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u)) 
    261261 
    262       u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u)) 
     262      u_star = MAX ( U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u)) , 1.E-9 )  !  (MAX => prevents FPE from stupid values from masked region later on) 
    263263      t_star = dt_zu*ztmp0 
    264264      q_star = dq_zu*ztmp0 
     
    283283         !!Inverse of Monin-Obukov length (1/L) : 
    284284         ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Monin-Obukhov length] 
    285          ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...) !#LOLO 
     285         ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...) 
    286286 
    287287         ztmp1 = u_star*u_star   ! u*^2 
     
    307307         ztmp2 = u_star/vkarmn*LOG(10./z0)                                 ! Neutral wind speed at 10m 
    308308         z0    = alfa_charn_3p6(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6) [ ztmp1==u*^2 ] 
     309         z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
     310 
    309311         ztmp1 = ( znu_a / (z0*u_star) )**0.72_wp     ! COARE3.6-specific! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific! 
    310312         z0t   = MIN( 1.6E-4_wp , 5.8E-5_wp*ztmp1 )   ! COARE3.6-specific! 
     313         z0t   = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
    311314 
    312315         !! Turbulent scales at zu : 
    313316         ztmp0   = psi_h_coare(zeta_u) 
    314          ztmp1   = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) ! #LOLO: in ztmp0, some use psi_h_coare(zeta_t) rather than psi_h_coare(zeta_t) ??? 
     317         ztmp1   = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) ! #LB: in ztmp0, some use psi_h_coare(zeta_t) rather than psi_h_coare(zeta_t) ??? 
    315318 
    316319         t_star = dt_zu*ztmp1 
    317320         q_star = dq_zu*ztmp1 
    318          u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u)) 
     321         u_star = MAX( U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u)) , 1.E-9 )  !  (MAX => prevents FPE from stupid values from masked region later on) 
    319322 
    320323         IF( .NOT. l_zt_equal_zu ) THEN 
Note: See TracChangeset for help on using the changeset viewer.