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 13655 for NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_ecmwf.F90 – NEMO

Ignore:
Timestamp:
2020-10-21T16:15:13+02:00 (3 years ago)
Author:
laurent
Message:

Commit all my dev of 2020!

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_ecmwf.F90

    r13460 r13655  
    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   !! 
     
    2424   !!                   returns the effective bulk wind speed at 10m 
    2525   !!---------------------------------------------------------------------- 
    26    USE oce             ! ocean dynamics and tracers 
    2726   USE dom_oce         ! ocean space and time domain 
    2827   USE phycst          ! physical constants 
    29    USE iom             ! I/O manager library 
    30    USE lib_mpp         ! distribued memory computing library 
    31    USE in_out_manager  ! I/O manager 
    32    USE prtctl          ! Print control 
    33    USE sbcwave, ONLY   :  cdn_wave ! wave module 
    34 #if defined key_si3 || defined key_cice 
    35    USE sbc_ice         ! Surface boundary condition: ice fields 
    36 #endif 
    37    USE lib_fortran     ! to use key_nosignedzero 
    38  
    39    USE sbc_oce         ! Surface boundary condition: ocean fields 
    40    USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB 
     28   USE lib_mpp,        ONLY: ctl_stop         ! distribued memory computing library 
     29   USE in_out_manager, ONLY: nit000  ! I/O manager 
     30   USE sbc_phy         ! all thermodynamics functions, rho_air, q_sat, etc... !LB 
    4131   USE sbcblk_skin_ecmwf ! cool-skin/warm layer scheme !LB 
    4232 
     
    4535 
    4636   PUBLIC :: SBCBLK_ALGO_ECMWF_INIT, TURB_ECMWF 
    47    !! * Substitutions 
    48 #  include "do_loop_substitute.h90" 
    4937 
    5038   !! ECMWF own values for given constants, taken form IFS documentation... 
    51    REAL(wp), PARAMETER ::   charn0 = 0.018    ! Charnock constant (pretty high value here !!! 
     39   REAL(wp), PARAMETER, PUBLIC :: charn0_ecmwf = 0.018_wp    ! Charnock constant (pretty high value here !!! 
    5240   !                                          !    =>  Usually 0.011 for moderate winds) 
    5341   REAL(wp), PARAMETER ::   zi0     = 1000.   ! scale height of the atmospheric boundary layer...1 
     
    5745   REAL(wp), PARAMETER ::   alpha_Q = 0.62    ! 
    5846 
    59    INTEGER , PARAMETER ::   nb_itt = 10             ! number of itterations 
     47   !! * Substitutions 
     48#  include "do_loop_substitute.h90" 
    6049 
    6150   !!---------------------------------------------------------------------- 
     
    9483 
    9584   SUBROUTINE turb_ecmwf( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 
    96       &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                           & 
    97       &                      Cdn, Chn, Cen,                                           & 
     85      &                      Cd, Ch, Ce, t_zu, q_zu, Ubzu,                            & 
     86      &                      nb_iter, Cdn, Chn, Cen,                                           & ! optional output 
    9887      &                      Qsw, rad_lw, slp, pdT_cs,                                & ! optionals for cool-skin (and warm-layer) 
    9988      &                      pdT_wl, pHz_wl )                                           ! optionals for warm-layer only 
     
    151140      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K] 
    152141      !!    *  q_zu   : specific humidity of air        //                    [kg/kg] 
    153       !!    *  U_blk  : bulk wind speed at zu                                 [m/s] 
     142      !!    *  Ubzu   : bulk wind speed at zu                                 [m/s] 
    154143      !! 
    155144      !! 
     
    171160      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K] 
    172161      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg] 
    173       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind speed at zu                     [m/s] 
    174       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    175       ! 
     162      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ubzu    ! bulk wind speed at zu                     [m/s] 
     163      ! 
     164      INTEGER , INTENT(in   ), OPTIONAL                     :: nb_iter  ! number of iterations 
     165      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CdN 
     166      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   ChN 
     167      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CeN 
    176168      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Qsw      !             [W/m^2] 
    177169      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2] 
     
    181173      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pHz_wl   !             [m] 
    182174      ! 
    183       INTEGER :: j_itt 
     175      INTEGER :: nbit, jit 
    184176      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
    185177      ! 
     
    197189      !!---------------------------------------------------------------------------------- 
    198190      IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl) 
     191 
     192      nbit = nb_iter0 
     193      IF( PRESENT(nb_iter) ) nbit = nb_iter 
    199194 
    200195      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision 
     
    227222      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 
    228223 
    229       U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution 
     224      Ubzu = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution 
    230225 
    231226      ztmp0   = LOG(    zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001) 
    232227      ztmp1   = LOG(10._wp*10000._wp) !       "                    "               " 
    233       u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10) 
    234  
    235       z0     = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star 
     228      u_star = 0.035_wp*Ubzu*ztmp1/ztmp0       ! (u* = 0.035*Un10) 
     229 
     230      z0     = charn0_ecmwf*u_star*u_star/grav + 0.11_wp*znu_a/u_star 
    236231      z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
    237232 
     
    239234      z0t    = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
    240235 
    241       Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd 
    242  
    243       ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 
    244  
    245       ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 
     236      Cd     = MAX( (vkarmn/ztmp0)**2 , Cx_min )   ! first guess of Cd 
     237 
     238      ztmp0 = vkarmn2/LOG(zt/z0t)/Cd 
     239 
     240      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, Ubzu ) ! Bulk Richardson Number (BRN) 
    246241 
    247242      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 
    248243      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 
    249       func_m = ztmp0*ztmp2 ! temporary array !! 
    250       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 
    251          &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))              !  BRN > 0 
    252       !#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" ! 
    253  
     244      func_h = (1._wp - ztmp1) *   ztmp0*ztmp2 / (1._wp - ztmp2*zi0*0.004_wp*Beta0**3/zu) & !  BRN < 0 
     245         &  +       ztmp1      * ( ztmp0*ztmp2 + 27._wp/9._wp*ztmp2*ztmp2 )                 !  BRN > 0 
     246       
    254247      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 
    255248      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h)) 
    256249 
    257       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) 
     250      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) 
    258251      t_star = dt_zu*ztmp0 
    259252      q_star = dq_zu*ztmp0 
     
    276269 
    277270 
    278       !! First guess of inverse of Monin-Obukov length (1/L) : 
     271      !! First guess of inverse of Obukov length (1/L) : 
    279272      Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star ) 
    280273 
    281       !! Functions such as  u* = U_blk*vkarmn/func_m 
     274      !! Functions such as  u* = Ubzu*vkarmn/func_m 
    282275      ztmp0 = zu*Linv 
    283276      func_m = LOG(zu) - LOG(z0)  - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv) 
     
    285278 
    286279      !! ITERATION BLOCK 
    287       DO j_itt = 1, nb_itt 
     280      DO jit = 1, nbit 
    288281 
    289282         !! Bulk Richardson Number at z=zu (Eq. 3.25) 
    290          ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 
    291  
    292          !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) : 
     283         ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, Ubzu ) ! Bulk Richardson Number (BRN) 
     284 
     285         !! New estimate of the inverse of the Obukhon length (Linv == zeta/zu) : 
    293286         Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1 
    294287         !! Note: it is slightly different that the L we would get with the usual 
     
    299292 
    300293         !! Need to update roughness lengthes: 
    301          u_star = U_blk*vkarmn/func_m 
     294         u_star = Ubzu*vkarmn/func_m 
    302295         ztmp2  = u_star*u_star 
    303296         ztmp1  = znu_a/u_star 
    304          z0     = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp) 
    305          z0t    = MIN( ABS( alpha_H*ztmp1                     ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 
    306          z0q    = MIN( ABS( alpha_Q*ztmp1                     ) , 0.001_wp) 
     297         z0     = MIN( ABS( alpha_M*ztmp1 + charn0_ecmwf*ztmp2/grav ) , 0.001_wp) 
     298         z0t    = MIN( ABS( alpha_H*ztmp1                           ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 
     299         z0q    = MIN( ABS( alpha_Q*ztmp1                           ) , 0.001_wp) 
    307300 
    308301         !! 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) 
    309302         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) 
    310303         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0 
    311          U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed 
    312          ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 
     304         Ubzu = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed 
     305         ! => 0.2 prevents Ubzu to be 0 in stable case when U_zu=0. 
    313306 
    314307 
     
    346339            !! Cool-skin contribution 
    347340 
    348             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, & 
     341            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, & 
    349342               &                   ztmp1, ztmp0,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp0 
    350343 
     
    359352         IF( l_use_wl ) THEN 
    360353            !! Warm-layer contribution 
    361             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, & 
     354            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, & 
    362355               &                   ztmp1, ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp2 
    363356            CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst ) 
     
    373366         ENDIF 
    374367 
    375       END DO !DO j_itt = 1, nb_itt 
    376  
    377       Cd = vkarmn*vkarmn/(func_m*func_m) 
    378       Ch = vkarmn*vkarmn/(func_m*func_h) 
    379       ztmp2 = log(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q 
    380       Ce = vkarmn*vkarmn/(func_m*ztmp2) 
    381  
    382       Cdn = vkarmn*vkarmn / (log(zu/z0 )*log(zu/z0 )) 
    383       Chn = vkarmn*vkarmn / (log(zu/z0t)*log(zu/z0t)) 
    384       Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q)) 
     368      END DO !DO jit = 1, nbit 
     369 
     370      Cd = MAX( vkarmn2/(func_m*func_m) , Cx_min ) 
     371      Ch = MAX( vkarmn2/(func_m*func_h) , Cx_min ) 
     372      ztmp2 = LOG(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q 
     373      Ce = MAX( vkarmn2/(func_m*ztmp2)  , Cx_min ) 
     374 
     375      IF(PRESENT(Cdn)) Cdn = MAX( vkarmn2 / (LOG(zu/z0 )*LOG(zu/z0 )) , Cx_min ) 
     376      IF(PRESENT(Chn)) Chn = MAX( vkarmn2 / (LOG(zu/z0t)*LOG(zu/z0t)) , Cx_min ) 
     377      IF(PRESENT(Cen)) Cen = MAX( vkarmn2 / (LOG(zu/z0q)*LOG(zu/z0q)) , Cx_min ) 
    385378 
    386379      IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 
     
    408401      ! 
    409402      INTEGER  ::   ji, jj    ! dummy loop indices 
    410       REAL(wp) :: zzeta, zx, ztmp, psi_unst, psi_stab, stab 
    411       !!---------------------------------------------------------------------------------- 
     403      REAL(wp) :: zta, zx2, zx, ztmp, zpsi_unst, zpsi_stab, zstab, zc 
     404      !!---------------------------------------------------------------------------------- 
     405      zc = 5._wp/0.35_wp 
    412406      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    413       ! 
    414       zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 
    415       ! 
    416       ! Unstable (Paulson 1970): 
    417       !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
    418       zx = SQRT(ABS(1._wp - 16._wp*zzeta)) 
    419       ztmp = 1._wp + SQRT(zx) 
    420       ztmp = ztmp*ztmp 
    421       psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   & 
    422          &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi 
    423       ! 
    424       ! Unstable: 
    425       ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
    426       psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & 
    427          &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp 
    428       ! 
    429       ! Combining: 
    430       stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
    431       ! 
    432       psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 
    433          &                +      stab  * psi_stab      ! (zzeta > 0) Stable 
    434       ! 
     407            ! 
     408            zta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 
     409 
     410            ! *** Unstable (Paulson 1970)    [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] : 
     411            zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 - 16z)^0.5 
     412            zx  = SQRT(zx2)                          ! (1 - 16z)^0.25 
     413            ztmp = 1._wp + zx 
     414            zpsi_unst = LOG( 0.125_wp*ztmp*ztmp*(1._wp + zx2) ) - 2._wp*ATAN( zx ) + 0.5_wp*rpi 
     415 
     416            ! *** Stable                   [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] : 
     417            zpsi_stab = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) & 
     418               &       - zta - 2._wp/3._wp*zc 
     419            ! 
     420            zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1 
     421            ! 
     422            psi_m_ecmwf(ji,jj) =         zstab  * zpsi_stab &  ! (zta > 0) Stable 
     423               &              + (1._wp - zstab) * zpsi_unst    ! (zta < 0) Unstable 
     424            ! 
    435425      END_2D 
    436426   END FUNCTION psi_m_ecmwf 
     
    452442      ! 
    453443      INTEGER  ::   ji, jj     ! dummy loop indices 
    454       REAL(wp) ::  zzeta, zx, psi_unst, psi_stab, stab 
    455       !!---------------------------------------------------------------------------------- 
     444      REAL(wp) ::  zta, zx2, zpsi_unst, zpsi_stab, zstab, zc 
     445      !!---------------------------------------------------------------------------------- 
     446      zc = 5._wp/0.35_wp 
    456447      ! 
    457448      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    458       ! 
    459       zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!): 
    460       ! 
    461       zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!! 
    462       !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 
    463       ! Unstable (Paulson 1970) : 
    464       psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
    465       ! 
    466       ! Stable: 
    467       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 
    468          &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp 
    469       ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 
    470       ! 
    471       stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
    472       ! 
    473       ! 
    474       psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable 
    475          &                +    stab    * psi_stab        ! (zzeta > 0) Stable 
    476       ! 
     449            ! 
     450            zta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!): 
     451            ! 
     452            ! *** Unstable (Paulson 1970)   [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] : 
     453            zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 -16z)^0.5 
     454            zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) ) 
     455            ! 
     456            ! *** Stable [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] : 
     457            zpsi_stab = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) & 
     458               &       - ABS(1._wp + 2._wp/3._wp*zta)**1.5_wp - 2._wp/3._wp*zc + 1._wp 
     459            ! 
     460            ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 
     461            ! 
     462            zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1 
     463            ! 
     464            psi_h_ecmwf(ji,jj) =         zstab  * zpsi_stab &  ! (zta > 0) Stable 
     465               &              + (1._wp - zstab) * zpsi_unst    ! (zta < 0) Unstable             
     466            ! 
    477467      END_2D 
    478468   END FUNCTION psi_h_ecmwf 
Note: See TracChangeset for help on using the changeset viewer.