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 14072 for NEMO/trunk/src/OCE/SBC/sbcblk_algo_ecmwf.F90 – NEMO

Ignore:
Timestamp:
2020-12-04T08:48:38+01:00 (4 years ago)
Author:
laurent
Message:

Merging branch "2020/dev_r13648_ASINTER-04_laurent_bulk_ice", ticket #2369

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/src/OCE/SBC/sbcblk_algo_ecmwf.F90

    r14007 r14072  
    55   !!   * bulk transfer coefficients C_D, C_E and C_H 
    66   !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed 
    7    !!   * the effective bulk wind speed at 10m U_blk 
     7   !!   * the effective bulk wind speed at 10m Ubzu 
    88   !!   => all these are used in bulk formulas in sbcblk.F90 
    99   !! 
    … …  
    1717   !!---------------------------------------------------------------------- 
    1818   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code 
    19    !!            4.2  !  2020-12  (G. Madec, E. Clementi) Charnock coeff from wave model 
     19   !!            4.2  !  2020-12  (L. Brodeau) Introduction of various air-ice bulk parameterizations + improvements 
    2020   !!---------------------------------------------------------------------- 
    2121 
    … …  
    2525   !!                   returns the effective bulk wind speed at 10m 
    2626   !!---------------------------------------------------------------------- 
    27    USE oce             ! ocean dynamics and tracers 
    2827   USE dom_oce         ! ocean space and time domain 
    2928   USE phycst          ! physical constants 
    30    USE iom             ! I/O manager library 
    31    USE lib_mpp         ! distribued memory computing library 
    32    USE in_out_manager  ! I/O manager 
    33    USE prtctl          ! Print control 
    34    USE sbcwave, ONLY   : charn ! wave module 
    35 #if defined key_si3 || defined key_cice 
    36    USE sbc_ice         ! Surface boundary condition: ice fields 
    37 #endif 
    38    USE lib_fortran     ! to use key_nosignedzero 
    39  
    40    USE sbc_oce         ! Surface boundary condition: ocean fields 
    41    USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB 
     29   USE lib_mpp,        ONLY: ctl_stop         ! distribued memory computing library 
     30   USE in_out_manager, ONLY: nit000  ! I/O manager 
     31   USE sbc_phy         ! Catalog of functions for physical/meteorological parameters in the marine boundary layer 
    4232   USE sbcblk_skin_ecmwf ! cool-skin/warm layer scheme !LB 
    4333 
    … …  
    4636 
    4737   PUBLIC :: SBCBLK_ALGO_ECMWF_INIT, TURB_ECMWF 
    48    !! * Substitutions 
    49 #  include "do_loop_substitute.h90" 
    5038 
    5139   !! ECMWF own values for given constants, taken form IFS documentation... 
    52    REAL(wp), PARAMETER ::   charn0 = 0.018    ! Charnock constant (pretty high value here !!! 
     40   REAL(wp), PARAMETER, PUBLIC :: charn0_ecmwf = 0.018_wp    ! Charnock constant (pretty high value here !!! 
    5341   !                                          !    =>  Usually 0.011 for moderate winds) 
    5442   REAL(wp), PARAMETER ::   zi0     = 1000.   ! scale height of the atmospheric boundary layer...1 
    … …  
    5846   REAL(wp), PARAMETER ::   alpha_Q = 0.62    ! 
    5947 
    60    INTEGER , PARAMETER ::   nb_itt = 10             ! number of itterations 
     48   !! * Substitutions 
     49#  include "do_loop_substitute.h90" 
    6150 
    6251   !!---------------------------------------------------------------------- 
    … …  
    9584 
    9685   SUBROUTINE turb_ecmwf( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 
    97       &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                           & 
    98       &                      Cdn, Chn, Cen,                                           & 
     86      &                      Cd, Ch, Ce, t_zu, q_zu, Ubzu,                            & 
     87      &                      nb_iter, Cdn, Chn, Cen,                                           & ! optional output 
    9988      &                      Qsw, rad_lw, slp, pdT_cs,                                & ! optionals for cool-skin (and warm-layer) 
    10089      &                      pdT_wl, pHz_wl )                                           ! optionals for warm-layer only 
    … …  
    152141      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K] 
    153142      !!    *  q_zu   : specific humidity of air        //                    [kg/kg] 
    154       !!    *  U_blk  : bulk wind speed at zu                                 [m/s] 
     143      !!    *  Ubzu   : bulk wind speed at zu                                 [m/s] 
    155144      !! 
    156145      !! 
    … …  
    172161      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K] 
    173162      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg] 
    174       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind speed at zu                     [m/s] 
    175       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    176       ! 
     163      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ubzu    ! bulk wind speed at zu                     [m/s] 
     164      ! 
     165      INTEGER , INTENT(in   ), OPTIONAL                     :: nb_iter  ! number of iterations 
     166      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CdN 
     167      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   ChN 
     168      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CeN 
    177169      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Qsw      !             [W/m^2] 
    178170      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2] 
    … …  
    182174      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pHz_wl   !             [m] 
    183175      ! 
    184       INTEGER :: j_itt 
     176      INTEGER :: nbit, jit 
    185177      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
    186178      ! 
    … …  
    198190      !!---------------------------------------------------------------------------------- 
    199191      IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl) 
     192 
     193      nbit = nb_iter0 
     194      IF( PRESENT(nb_iter) ) nbit = nb_iter 
    200195 
    201196      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision 
    … …  
    228223      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 
    229224 
    230       U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution 
     225      Ubzu = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution 
    231226 
    232227      ztmp0   = LOG(    zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001) 
    233228      ztmp1   = LOG(10._wp*10000._wp) !       "                    "               " 
    234       u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10) 
    235  
    236       IF (ln_charn)  THEN          !  Charnock value if wave coupling 
    237          z0     = charn*u_star*u_star/grav + 0.11_wp*znu_a/u_star 
    238       ELSE 
    239          z0     = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star 
    240       ENDIF 
    241  
     229      u_star = 0.035_wp*Ubzu*ztmp1/ztmp0       ! (u* = 0.035*Un10) 
     230 
     231      z0     = charn0_ecmwf*u_star*u_star/grav + 0.11_wp*znu_a/u_star 
    242232      z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
    243233 
    … …  
    245235      z0t    = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
    246236 
    247       Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd 
    248  
    249       ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 
    250  
    251       ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 
     237      Cd     = MAX( (vkarmn/ztmp0)**2 , Cx_min )   ! first guess of Cd 
     238 
     239      ztmp0 = vkarmn2/LOG(zt/z0t)/Cd 
     240 
     241      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, Ubzu ) ! Bulk Richardson Number (BRN) 
    252242 
    253243      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 
    254244      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 
    255       func_m = ztmp0*ztmp2 ! temporary array !! 
    256       func_h = (1._wp-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0 ! temporary array !!! func_h == zeta_u 
    257          &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))              !  BRN > 0 
    258       !#LB: should make sure that the "func_m" of "27./9.*ztmp2/func_m" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 
     245      func_h = (1._wp - ztmp1) *   ztmp0*ztmp2 / (1._wp - ztmp2*zi0*0.004_wp*Beta0**3/zu) & !  BRN < 0 
     246         &  +       ztmp1      * ( ztmp0*ztmp2 + 27._wp/9._wp*ztmp2*ztmp2 )                 !  BRN > 0 
    259247 
    260248      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 
    261249      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h)) 
    262250 
    263       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) 
     251      u_star = MAX ( Ubzu*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_ecmwf(func_h)) , 1.E-9 )  !  (MAX => prevents FPE from stupid values from masked region later on) 
    264252      t_star = dt_zu*ztmp0 
    265253      q_star = dq_zu*ztmp0 
    … …  
    282270 
    283271 
    284       !! First guess of inverse of Monin-Obukov length (1/L) : 
     272      !! First guess of inverse of Obukov length (1/L) : 
    285273      Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star ) 
    286274 
    287       !! Functions such as  u* = U_blk*vkarmn/func_m 
     275      !! Functions such as  u* = Ubzu*vkarmn/func_m 
    288276      ztmp0 = zu*Linv 
    289277      func_m = LOG(zu) - LOG(z0)  - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv) 
    … …  
    291279 
    292280      !! ITERATION BLOCK 
    293       DO j_itt = 1, nb_itt 
     281      DO jit = 1, nbit 
    294282 
    295283         !! Bulk Richardson Number at z=zu (Eq. 3.25) 
    296          ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 
    297  
    298          !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) : 
     284         ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, Ubzu ) ! Bulk Richardson Number (BRN) 
     285 
     286         !! New estimate of the inverse of the Obukhon length (Linv == zeta/zu) : 
    299287         Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1 
    300288         !! Note: it is slightly different that the L we would get with the usual 
    … …  
    305293 
    306294         !! Need to update roughness lengthes: 
    307          u_star = U_blk*vkarmn/func_m 
     295         u_star = Ubzu*vkarmn/func_m 
    308296         ztmp2  = u_star*u_star 
    309297         ztmp1  = znu_a/u_star 
    310          IF (ln_charn) THEN     ! Charnock value if wave coupling 
    311             z0  = MIN( ABS( alpha_M*ztmp1 + charn*ztmp2/grav ) , 0.001_wp)          
    312          ELSE 
    313             z0     = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp) 
    314          ENDIF 
    315          z0t    = MIN( ABS( alpha_H*ztmp1                     ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 
    316          z0q    = MIN( ABS( alpha_Q*ztmp1                     ) , 0.001_wp) 
     298         z0     = MIN( ABS( alpha_M*ztmp1 + charn0_ecmwf*ztmp2/grav ) , 0.001_wp) 
     299         z0t    = MIN( ABS( alpha_H*ztmp1                           ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 
     300         z0q    = MIN( ABS( alpha_Q*ztmp1                           ) , 0.001_wp) 
    317301 
    318302         !! Update wind at zu with convection-related wind gustiness in unstable conditions (Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8) 
    319303         ztmp2 = Beta0*Beta0*ztmp2*(MAX(-zi0*Linv/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution  (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1) 
    320304         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0 
    321          U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed 
    322          ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 
     305         Ubzu = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed 
     306         ! => 0.2 prevents Ubzu to be 0 in stable case when U_zu=0. 
    323307 
    324308 
    … …  
    356340            !! Cool-skin contribution 
    357341 
    358             CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 
     342            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, & 
    359343               &                   ztmp1, ztmp0,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp0 
    360344 
    … …  
    369353         IF( l_use_wl ) THEN 
    370354            !! Warm-layer contribution 
    371             CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 
     355            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, & 
    372356               &                   ztmp1, ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp2 
    373357            CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst ) 
    … …  
    383367         ENDIF 
    384368 
    385       END DO !DO j_itt = 1, nb_itt 
    386  
    387       Cd = vkarmn*vkarmn/(func_m*func_m) 
    388       Ch = vkarmn*vkarmn/(func_m*func_h) 
    389       ztmp2 = log(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q 
    390       Ce = vkarmn*vkarmn/(func_m*ztmp2) 
    391  
    392       Cdn = vkarmn*vkarmn / (log(zu/z0 )*log(zu/z0 )) 
    393       Chn = vkarmn*vkarmn / (log(zu/z0t)*log(zu/z0t)) 
    394       Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q)) 
     369      END DO !DO jit = 1, nbit 
     370 
     371      Cd = MAX( vkarmn2/(func_m*func_m) , Cx_min ) 
     372      Ch = MAX( vkarmn2/(func_m*func_h) , Cx_min ) 
     373      ztmp2 = LOG(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q 
     374      Ce = MAX( vkarmn2/(func_m*ztmp2)  , Cx_min ) 
     375 
     376      IF(PRESENT(Cdn)) Cdn = MAX( vkarmn2 / (LOG(zu/z0 )*LOG(zu/z0 )) , Cx_min ) 
     377      IF(PRESENT(Chn)) Chn = MAX( vkarmn2 / (LOG(zu/z0t)*LOG(zu/z0t)) , Cx_min ) 
     378      IF(PRESENT(Cen)) Cen = MAX( vkarmn2 / (LOG(zu/z0q)*LOG(zu/z0q)) , Cx_min ) 
    395379 
    396380      IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 
    … …  
    418402      ! 
    419403      INTEGER  ::   ji, jj    ! dummy loop indices 
    420       REAL(wp) :: zzeta, zx, ztmp, psi_unst, psi_stab, stab 
    421       !!---------------------------------------------------------------------------------- 
     404      REAL(wp) :: zta, zx2, zx, ztmp, zpsi_unst, zpsi_stab, zstab, zc 
     405      !!---------------------------------------------------------------------------------- 
     406      zc = 5._wp/0.35_wp 
    422407      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    423       ! 
    424       zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 
    425       ! 
    426       ! Unstable (Paulson 1970): 
    427       !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
    428       zx = SQRT(ABS(1._wp - 16._wp*zzeta)) 
    429       ztmp = 1._wp + SQRT(zx) 
    430       ztmp = ztmp*ztmp 
    431       psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   & 
    432          &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi 
    433       ! 
    434       ! Unstable: 
    435       ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
    436       psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & 
    437          &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp 
    438       ! 
    439       ! Combining: 
    440       stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
    441       ! 
    442       psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 
    443          &                +      stab  * psi_stab      ! (zzeta > 0) Stable 
    444       ! 
     408            ! 
     409            zta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 
     410 
     411            ! *** Unstable (Paulson 1970)    [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] : 
     412            zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 - 16z)^0.5 
     413            zx  = SQRT(zx2)                          ! (1 - 16z)^0.25 
     414            ztmp = 1._wp + zx 
     415            zpsi_unst = LOG( 0.125_wp*ztmp*ztmp*(1._wp + zx2) ) - 2._wp*ATAN( zx ) + 0.5_wp*rpi 
     416 
     417            ! *** Stable                   [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] : 
     418            zpsi_stab = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) & 
     419               &       - zta - 2._wp/3._wp*zc 
     420            ! 
     421            zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1 
     422            ! 
     423            psi_m_ecmwf(ji,jj) =         zstab  * zpsi_stab &  ! (zta > 0) Stable 
     424               &              + (1._wp - zstab) * zpsi_unst    ! (zta < 0) Unstable 
     425            ! 
    445426      END_2D 
    446427   END FUNCTION psi_m_ecmwf 
    … …  
    462443      ! 
    463444      INTEGER  ::   ji, jj     ! dummy loop indices 
    464       REAL(wp) ::  zzeta, zx, psi_unst, psi_stab, stab 
    465       !!---------------------------------------------------------------------------------- 
     445      REAL(wp) ::  zta, zx2, zpsi_unst, zpsi_stab, zstab, zc 
     446      !!---------------------------------------------------------------------------------- 
     447      zc = 5._wp/0.35_wp 
    466448      ! 
    467449      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    468       ! 
    469       zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!): 
    470       ! 
    471       zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!! 
    472       !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 
    473       ! Unstable (Paulson 1970) : 
    474       psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
    475       ! 
    476       ! Stable: 
    477       psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
    478          &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp 
    479       ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 
    480       ! 
    481       stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
    482       ! 
    483       ! 
    484       psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable 
    485          &                +    stab    * psi_stab        ! (zzeta > 0) Stable 
    486       ! 
     450            ! 
     451            zta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!): 
     452            ! 
     453            ! *** Unstable (Paulson 1970)   [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] : 
     454            zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 -16z)^0.5 
     455            zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) ) 
     456            ! 
     457            ! *** Stable [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] : 
     458            zpsi_stab = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) & 
     459               &       - ABS(1._wp + 2._wp/3._wp*zta)**1.5_wp - 2._wp/3._wp*zc + 1._wp 
     460            ! 
     461            ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 
     462            ! 
     463            zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1 
     464            ! 
     465            psi_h_ecmwf(ji,jj) =         zstab  * zpsi_stab &  ! (zta > 0) Stable 
     466               &              + (1._wp - zstab) * zpsi_unst    ! (zta < 0) Unstable 
     467            ! 
    487468      END_2D 
    488469   END FUNCTION psi_h_ecmwf 
Note: See TracChangeset for help on using the changeset viewer.