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 – NEMO

Changeset 11785 for NEMO


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

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

Location:
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90

    r11772 r11785  
    4545   PUBLIC :: COARE3P0_INIT, TURB_COARE3P0 
    4646 
    47    !                                              !! COARE own values for given constants: 
    48    REAL(wp), PARAMETER ::   zi0     = 600._wp     ! scale height of the atmospheric boundary layer... 
    49    REAL(wp), PARAMETER ::   Beta0   =   1.25_wp   ! gustiness parameter 
     47   !! COARE own values for given constants: 
     48   REAL(wp), PARAMETER :: zi0   = 600._wp     ! scale height of the atmospheric boundary layer... 
     49   REAL(wp), PARAMETER :: Beta0 =  1.25_wp    ! gustiness parameter 
    5050 
    5151   INTEGER , PARAMETER ::   nb_itt = 10             ! number of itterations 
    … …  
    7171         ierr = 0 
    7272         ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), Hz_wl(jpi,jpj), dT_wl(jpi,jpj), STAT=ierr ) 
    73          !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of Tau_ac and Qnt_ac failed!' 
     73         !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of Tau_ac, Qnt_ac, dT_wl & Hz_wl failed!' 
    7474         Tau_ac(:,:) = 0._wp 
    7575         Qnt_ac(:,:) = 0._wp 
     76         dT_wl(:,:)  = 0._wp 
    7677         Hz_wl(:,:)  = Hwl_max 
    77          dT_wl(:,:)  = 0._wp 
    7878      END IF 
    7979      !! 
    … …  
    8181         ierr = 0 
    8282         ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 
    83          !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of dT_cs and Qnt_ac failed!' 
     83         !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of dT_cs failed!' 
    8484         dT_cs(:,:) = -0.25_wp  ! First guess of skin correction 
    8585      END IF 
    … …  
    132132      !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2] 
    133133      !!    *  slp    : sea-level pressure                                    [Pa] 
     134      !! 
     135      !! OPTIONAL OUTPUT: 
     136      !! ---------------- 
    134137      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K] 
    135138      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K] 
    … …  
    170173      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa] 
    171174      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs 
    172       ! 
    173175      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K] 
    174176      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pHz_wl   !             [m] 
    … …  
    193195      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p0@sbcblk_algo_coare3p0' 
    194196      !!---------------------------------------------------------------------------------- 
    195  
    196197      IF ( kt == nit000 ) CALL COARE3P0_INIT(l_use_cs, l_use_wl) 
    197198 
    … …  
    208209 
    209210      IF ( l_use_wl ) THEN 
    210          IF(.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) )) THEN 
     211         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN 
    211212            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use warm-layer param!'; STOP 
    212213         END IF 
    … …  
    217218         zsst = T_s ! backing up the bulk SST 
    218219         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction 
    219          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 
    220221      END IF 
    221222 
    … …  
    238239 
    239240      z0     = alfa_charn_3p0(U_zu)*u_star*u_star/grav + 0.11_wp*znu_a/u_star 
    240       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 
    241243      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 
    242       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) 
    243245 
    244246      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd 
    … …  
    258260      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u)) 
    259261 
    260       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) 
    261263      t_star = dt_zu*ztmp0 
    262264      q_star = dq_zu*ztmp0 
    … …  
    281283         !!Inverse of Monin-Obukov length (1/L) : 
    282284         ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Monin-Obukhov length] 
    283          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...) 
    284286 
    285287         ztmp1 = u_star*u_star   ! u*^2 
    … …  
    305307         ztmp2 = u_star/vkarmn*LOG(10./z0)                                 ! Neutral wind speed at 10m 
    306308         z0    = alfa_charn_3p0(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6) 
     309         z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
     310 
    307311         ztmp1 = ( znu_a / (z0*u_star) )**0.6_wp    ! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific! 
    308312         z0t   = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1 ) ! Scalar roughness for both theta and q (eq.28) #LOLO: some use 1.15 not 1.1 !!! 
     313         z0t   = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
    309314 
    310315         !! Turbulent scales at zu : 
    311316         ztmp0   = psi_h_coare(zeta_u) 
    312          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) ??? 
    313318 
    314319         t_star = dt_zu*ztmp1 
    315320         q_star = dq_zu*ztmp1 
    316          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) 
    317322 
    318323         IF( .NOT. l_zt_equal_zu ) THEN 

  • 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 

  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90

    r11772 r11785  
    4646   PUBLIC :: ECMWF_INIT, TURB_ECMWF 
    4747 
    48    !                   !! ECMWF own values for given constants, taken form IFS documentation... 
     48   !! ECMWF own values for given constants, taken form IFS documentation... 
    4949   REAL(wp), PARAMETER ::   charn0 = 0.018    ! Charnock constant (pretty high value here !!! 
    5050   !                                          !    =>  Usually 0.011 for moderate winds) 
    … …  
    213213 
    214214      IF ( l_use_wl ) THEN 
    215          IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) THEN 
     215         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN 
    216216            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use warm-layer param!'; STOP 
    217217         END IF 
    … …  
    222222         zsst = T_s ! backing up the bulk SST 
    223223         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction 
    224          q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!! 
     224         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 
    225225      END IF 
    226226 
    … …  
    245245 
    246246      z0     = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star 
    247       z0     = MIN(ABS(z0), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
     247      z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
     248 
    248249      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 
    249       z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
     250      z0t    = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
    250251 
    251252      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd 
    … …  
    265266      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h)) 
    266267 
    267       u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_ecmwf(func_h)) 
     268      u_star = MAX ( U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_ecmwf(func_h)) , 1.E-9 )  !  (MAX => prevents FPE from stupid values from masked region later on) 
    268269      t_star = dt_zu*ztmp0 
    269270      q_star = dq_zu*ztmp0 
Note: See TracChangeset for help on using the changeset viewer.