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

Changeset 13017


Ignore:
Timestamp:
2020-06-03T12:55:26+02:00 (5 years ago)
Author:
laurent
Message:

Bugfix to prevent "TURB_NCAR()" to return inconsistsent Cd and CdN when "ln_wave==T".

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbcblk_algo_ncar.F90

    r10190 r13017  
    1111   !! 
    1212   !!       Routine turb_ncar maintained and developed in AeroBulk 
    13    !!                     (http://aerobulk.sourceforge.net/) 
     13   !!                     (https://github.com/brodeau/aerobulk/) 
    1414   !! 
    15    !!                         L. Brodeau, 2015 
     15   !!                         L. Brodeau, 2020 
    1616   !!===================================================================== 
    17    !! History :  3.6  !  2016-02  (L.Brodeau) successor of old turb_ncar of former sbcblk_core.F90 
     17   !! History :  4.0  !  2020-06  (L.Brodeau) successor of old turb_ncar of former sbcblk_core.F90 
    1818   !!---------------------------------------------------------------------- 
    1919 
    … …  
    4242   PRIVATE 
    4343 
    44    PUBLIC ::   TURB_NCAR   ! called by sbcblk.F90 
     44   PUBLIC :: TURB_NCAR   ! called by sbcblk.F90 
    4545 
    4646   !                              ! NCAR own values for given constants: 
    4747   REAL(wp), PARAMETER ::   rctv0 = 0.608   ! constant to obtain virtual temperature... 
    48     
     48 
     49   INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations 
     50 
    4951   !!---------------------------------------------------------------------- 
    5052CONTAINS 
    5153 
    5254   SUBROUTINE turb_ncar( zt, zu, sst, t_zt, ssq, q_zt, U_zu, & 
    53       &                  Cd, Ch, Ce, t_zu, q_zu, U_blk,      & 
    54       &                  Cdn, Chn, Cen                       ) 
    55       !!---------------------------------------------------------------------------------- 
     55      &                  Cd, Ch, Ce, t_zu, q_zu, Ub,         & 
     56      &                  CdN, ChN, CeN                       ) 
     57      !!---------------------------------------------------------------------- 
    5658      !!                      ***  ROUTINE  turb_ncar  *** 
    5759      !! 
    … …  
    5961      !!                fluxes according to Large & Yeager (2004) and Large & Yeager (2008) 
    6062      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 
    61       !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas 
    62       !! 
    63       !! ** Method : Monin Obukhov Similarity Theory 
    64       !!             + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10) 
    65       !! 
    66       !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008 
    67       !! 
    68       !! ** Last update: Laurent Brodeau, June 2014: 
    69       !!    - handles both cases zt=zu and zt/=zu 
    70       !!    - optimized: less 2D arrays allocated and less operations 
    71       !!    - better first guess of stability by checking air-sea difference of virtual temperature 
    72       !!       rather than temperature difference only... 
    73       !!    - added function "cd_neutral_10m" that uses the improved parametrization of 
    74       !!      Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions! 
    75       !!    - using code-wide physical constants defined into "phycst.mod" rather than redifining them 
    76       !!      => 'vkarmn' and 'grav' 
    77       !! 
    78       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     63      !!                Returns the effective bulk wind speed at zu to be used in the bulk formulas 
    7964      !! 
    8065      !! INPUT : 
    8166      !! ------- 
    8267      !!    *  zt   : height for temperature and spec. hum. of air            [m] 
    83       !!    *  zu   : height for wind speed (generally 10m)                   [m] 
    84       !!    *  U_zu : scalar wind speed at 10m                                [m/s] 
    85       !!    *  sst  : SST                                                     [K] 
     68      !!    *  zu   : height for wind speed (usually 10m)                     [m] 
     69      !!    *  sst  : bulk SST                                                [K] 
    8670      !!    *  t_zt : potential air temperature at zt                         [K] 
    8771      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg] 
    8872      !!    *  q_zt : specific humidity of air at zt                          [kg/kg] 
    89       !! 
     73      !!    *  U_zu : scalar wind speed at zu                                 [m/s] 
    9074      !! 
    9175      !! OUTPUT : 
    … …  
    9680      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K] 
    9781      !!    *  q_zu   : specific humidity of air        //                    [kg/kg] 
    98       !!    *  U_blk  : bulk wind at 10m                                      [m/s] 
     82      !!    *  Ub  : bulk wind speed at zu                                 [m/s] 
     83      !! 
     84      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    9985      !!---------------------------------------------------------------------------------- 
    10086      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m] 
    … …  
    10389      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin] 
    10490      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity           [kg/kg] 
    105       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity                   [kg/kg] 
     91      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg] 
    10692      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s] 
    10793      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau) 
    … …  
    11096      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K] 
    11197      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg] 
    112       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s] 
    113       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    114       ! 
    115       INTEGER ::   j_itt 
    116       LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
    117       INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations 
     98      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ub    ! bulk wind speed at zu                     [m/s] 
     99      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   CdN, ChN, CeN ! neutral transfer coefficients 
     100      ! 
     101      INTEGER :: j_itt 
     102      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
    118103      ! 
    119104      REAL(wp), DIMENSION(jpi,jpj) ::   Cx_n10        ! 10m neutral latent/sensible coefficient 
    120       REAL(wp), DIMENSION(jpi,jpj) ::   sqrt_Cd_n10   ! root square of Cd_n10 
     105      REAL(wp), DIMENSION(jpi,jpj) ::   sqrtCdN10   ! root square of Cd_n10 
    121106      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu 
    122107      REAL(wp), DIMENSION(jpi,jpj) ::   zpsi_h_u 
    123108      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2 
    124       REAL(wp), DIMENSION(jpi,jpj) ::   stab          ! stability test integer 
    125       !!---------------------------------------------------------------------------------- 
    126       ! 
    127       l_zt_equal_zu = .FALSE. 
    128       IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
    129  
    130       U_blk = MAX( 0.5 , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 
    131  
    132       !! First guess of stability: 
    133       ztmp0 = t_zt*(1. + rctv0*q_zt) - sst*(1. + rctv0*ssq) ! air-sea difference of virtual pot. temp. at zt 
    134       stab  = 0.5 + sign(0.5,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable 
    135  
    136       !! Neutral coefficients at 10m: 
     109      REAL(wp), DIMENSION(jpi,jpj) ::   sqrtCd 
     110      !!---------------------------------------------------------------------------------- 
     111 
     112      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) 
     113 
     114      Ub = MAX( 0.5_wp , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 
     115 
     116      !! Neutral drag coefficient at zu: 
    137117      IF( ln_cdgw ) THEN      ! wave drag case 
    138          cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) 
    139          ztmp0   (:,:) = cdn_wave(:,:) 
     118         CdN = MAX( cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) , 0.1E-3_wp ) 
    140119      ELSE 
    141          ztmp0 = cd_neutral_10m( U_blk ) 
     120         CdN = CD_N10_NCAR( Ub ) 
    142121      ENDIF 
    143  
    144       sqrt_Cd_n10 = SQRT( ztmp0 ) 
     122      sqrtCdN10 = SQRT( CdN ) 
    145123 
    146124      !! Initializing transf. coeff. with their first guess neutral equivalents : 
    147       Cd = ztmp0 
    148       Ce = 1.e-3*( 34.6 * sqrt_Cd_n10 ) 
    149       Ch = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) 
    150       stab = sqrt_Cd_n10   ! Temporaty array !!! stab == SQRT(Cd) 
    151   
    152       IF( ln_cdgw )   Cen = Ce  ; Chn = Ch 
     125      Cd = CdN 
     126      Ce = CE_N10_NCAR( sqrtCdN10 ) 
     127      ztmp0 = 0.5_wp + SIGN(0.5_wp, virt_temp(t_zt, q_zt) - virt_temp(sst, ssq)) ! we guess stability based on delta of virt. pot. temp. 
     128      Ch = CH_N10_NCAR( sqrtCdN10 , ztmp0 ) 
     129      sqrtCd = sqrtCdN10 
    153130 
    154131      !! Initializing values at z_u with z_t values: 
    155       t_zu = t_zt   ;   q_zu = q_zt 
    156  
    157       !!  * Now starting iteration loop 
    158       DO j_itt=1, nb_itt 
     132      t_zu = t_zt 
     133      q_zu = q_zt 
     134 
     135      !! ITERATION BLOCK 
     136      DO j_itt = 1, nb_itt 
    159137         ! 
    160138         ztmp1 = t_zu - sst   ! Updating air/sea differences 
    161139         ztmp2 = q_zu - ssq 
    162140 
    163          ! Updating turbulent scales :   (L&Y 2004 eq. (7)) 
    164          ztmp1  = Ch/stab*ztmp1    ! theta*   (stab == SQRT(Cd)) 
    165          ztmp2  = Ce/stab*ztmp2    ! q*       (stab == SQRT(Cd)) 
    166  
    167          ztmp0 = 1. + rctv0*q_zu      ! multiply this with t and you have the virtual temperature 
    168  
    169          ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 
    170          ztmp0 =  (grav*vkarmn/(t_zu*ztmp0)*(ztmp1*ztmp0 + rctv0*t_zu*ztmp2)) / (Cd*U_blk*U_blk) 
    171          !                                                      ( Cd*U_blk*U_blk is U*^2 at zu ) 
     141         ! Updating turbulent scales :   (L&Y 2004 Eq. (7)) 
     142         ztmp0 = sqrtCd*Ub          ! u* 
     143         ztmp1 = Ch/sqrtCd*ztmp1    ! theta* 
     144         ztmp2 = Ce/sqrtCd*ztmp2    ! q* 
     145 
     146         ! Estimate the inverse of Obukov length (1/L) at height zu: 
     147         ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 ) 
    172148 
    173149         !! Stability parameters : 
    174          zeta_u   = zu*ztmp0   ;  zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 
    175          zpsi_h_u = psi_h( zeta_u ) 
    176  
    177          !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c)) 
     150         zeta_u   = zu*ztmp0 
     151         zeta_u   = sign( min(abs(zeta_u),10._wp), zeta_u ) 
     152 
     153         !! Shifting temperature and humidity at zu (L&Y 2004 Eq. (9b-9c)) 
    178154         IF( .NOT. l_zt_equal_zu ) THEN 
    179             !! Array 'stab' is free for the moment so using it to store 'zeta_t' 
    180             stab = zt*ztmp0 ;  stab = SIGN( MIN(ABS(stab),10.0), stab )  ! Temporaty array stab == zeta_t !!! 
    181             stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab)                   ! stab just used as temp array again! 
    182             t_zu = t_zt - ztmp1/vkarmn*stab    ! ztmp1 is still theta*  L&Y 2004 eq.(9b) 
    183             q_zu = q_zt - ztmp2/vkarmn*stab    ! ztmp2 is still q*      L&Y 2004 eq.(9c) 
    184             q_zu = max(0., q_zu) 
     155            ztmp0 = zt*ztmp0 ! zeta_t ! 
     156            ztmp0 = SIGN( MIN(ABS(ztmp0),10._wp), ztmp0 )  ! Temporaty array ztmp0 == zeta_t !!! 
     157            ztmp0 = LOG(zt/zu) + psi_h_ncar(zeta_u) - psi_h_ncar(ztmp0)                   ! ztmp0 just used as temp array again! 
     158            t_zu = t_zt - ztmp1/vkarmn*ztmp0    ! ztmp1 is still theta*  L&Y 2004 Eq. (9b) 
     159            !! 
     160            q_zu = q_zt - ztmp2/vkarmn*ztmp0    ! ztmp2 is still q*      L&Y 2004 Eq. (9c) 
     161            q_zu = MAX(0._wp, q_zu) 
    185162         END IF 
    186163 
    187          ztmp2 = psi_m(zeta_u) 
    188          IF( ln_cdgw ) THEN      ! surface wave case 
    189             stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 )  ! (stab == SQRT(Cd)) 
    190             Cd   = stab * stab 
    191             ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
    192             ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd)) 
    193             ztmp1 = 1. + Chn * ztmp0      
    194             Ch    = Chn * ztmp2 / ztmp1  ! L&Y 2004 eq. (10b) 
    195             ztmp1 = 1. + Cen * ztmp0 
    196             Ce    = Cen * ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
    197  
     164         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 Eq. 9a)... 
     165         !   In very rare low-wind conditions, the old way of estimating the 
     166         !   neutral wind speed at 10m leads to a negative value that causes the code 
     167         !   to crash. To prevent this a threshold of 0.25m/s is imposed. 
     168         ztmp2 = psi_m_ncar(zeta_u) 
     169         ztmp0 = MAX( 0.25_wp , UN10_from_CD(zu, Ub, Cd, ppsi=ztmp2) ) ! U_n10 (ztmp2 == psi_m_ncar(zeta_u)) 
     170 
     171         IF( ln_cdgw ) THEN      ! wave drag case 
     172            CdN = MAX( cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) , 0.1E-3_wp ) 
    198173         ELSE 
    199             ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 
    200             !   In very rare low-wind conditions, the old way of estimating the 
    201             !   neutral wind speed at 10m leads to a negative value that causes the code 
    202             !   to crash. To prevent this a threshold of 0.25m/s is imposed. 
    203             ztmp0 = MAX( 0.25 , U_blk/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u)) 
    204             ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10 
    205             Cdn(:,:) = ztmp0 
    206             sqrt_Cd_n10 = sqrt(ztmp0) 
    207  
    208             stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability 
    209             Cx_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10) 
    210             Chn(:,:) = Cx_n10 
    211  
    212             !! Update of transfer coefficients: 
    213             ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u)) 
    214             Cd      = ztmp0 / ( ztmp1*ztmp1 ) 
    215             stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 
    216  
    217             ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
    218             ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd)) 
    219             ztmp1 = 1. + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10) 
    220             Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b) 
    221  
    222             Cx_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10 
    223             Cen(:,:) = Cx_n10 
    224             ztmp1 = 1. + Cx_n10*ztmp0 
    225             Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
    226             ENDIF 
    227          ! 
    228       END DO 
    229       ! 
     174            CdN   = CD_N10_NCAR(ztmp0)                                       ! Cd_n10 
     175         END IF 
     176         sqrtCdN10 = SQRT(CdN) 
     177 
     178         !! Update of transfer coefficients: 
     179         ztmp1  = 1._wp + sqrtCdN10/vkarmn*(LOG(zu/10._wp) - ztmp2)   ! L&Y 2004 Eq. (10a) (ztmp2 == psi_m(zeta_u)) 
     180         Cd     = MAX( CdN / ( ztmp1*ztmp1 ) , 0.1E-3_wp ) 
     181         sqrtCd = SQRT( Cd ) 
     182 
     183         ztmp0  = ( LOG(zu/10._wp) - psi_h_ncar(zeta_u) ) / vkarmn / sqrtCdN10 
     184         ztmp2  = sqrtCd / sqrtCdN10 
     185 
     186         ztmp1  = 0.5_wp + sign(0.5_wp,zeta_u) ! stability flag 
     187         ChN    = CH_N10_NCAR( sqrtCdN10 , ztmp1 ) 
     188         ztmp1  = 1._wp + ChN*ztmp0 
     189         Ch     = MAX( ChN*ztmp2 / ztmp1 , 0.1E-3_wp )   ! L&Y 2004 Eq. (10b) 
     190 
     191         CeN = CE_N10_NCAR( sqrtCdN10 ) 
     192         ztmp1  = 1._wp + CeN*ztmp0 
     193         Ce     = MAX( CeN*ztmp2 / ztmp1 , 0.1E-3_wp )  ! L&Y 2004 Eq. (10c) 
     194 
     195      END DO !DO j_itt = 1, nb_itt 
     196 
    230197   END SUBROUTINE turb_ncar 
    231198 
    232199 
    233    FUNCTION cd_neutral_10m( pw10 ) 
    234       !!----------------------------------------------------------------------------------       
     200   FUNCTION CD_N10_NCAR( pw10 ) 
     201      !!---------------------------------------------------------------------------------- 
    235202      !! Estimate of the neutral drag coefficient at 10m as a function 
    236203      !! of neutral wind  speed at 10m 
    237204      !! 
    238       !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b) 
    239       !! 
    240       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     205      !! Origin: Large & Yeager 2008, Eq. (11) 
     206      !! 
     207      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    241208      !!---------------------------------------------------------------------------------- 
    242209      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10           ! scalar wind speed at 10m (m/s) 
    243       REAL(wp), DIMENSION(jpi,jpj)             :: cd_neutral_10m 
     210      REAL(wp), DIMENSION(jpi,jpj)             :: CD_N10_NCAR 
    244211      ! 
    245212      INTEGER  ::     ji, jj     ! dummy loop indices 
    … …  
    255222            ! 
    256223            ! When wind speed > 33 m/s => Cyclone conditions => special treatment 
    257             zgt33 = 0.5 + SIGN( 0.5, (zw - 33.) )   ! If pw10 < 33. => 0, else => 1 
    258             ! 
    259             cd_neutral_10m(ji,jj) = 1.e-3 * ( & 
    260                &       (1. - zgt33)*( 2.7/zw + 0.142 + zw/13.09 - 3.14807E-10*zw6) & ! wind <  33 m/s 
    261                &      +    zgt33   *      2.34 )                                     ! wind >= 33 m/s 
    262             ! 
    263             cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6) 
     224            zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) )   ! If pw10 < 33. => 0, else => 1 
     225            ! 
     226            CD_N10_NCAR(ji,jj) = 1.e-3_wp * ( & 
     227               &       (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind <  33 m/s 
     228               &      +    zgt33   *      2.34_wp )                                                 ! wind >= 33 m/s 
     229            ! 
     230            CD_N10_NCAR(ji,jj) = MAX( CD_N10_NCAR(ji,jj), 0.1E-3_wp ) 
    264231            ! 
    265232         END DO 
    266233      END DO 
    267234      ! 
    268    END FUNCTION cd_neutral_10m 
    269  
    270  
    271    FUNCTION psi_m( pzeta ) 
     235   END FUNCTION CD_N10_NCAR 
     236 
     237 
     238 
     239   FUNCTION CH_N10_NCAR( psqrtcdn10 , pstab ) 
     240      !!---------------------------------------------------------------------------------- 
     241      !! Estimate of the neutral heat transfer coefficient at 10m      !! 
     242      !! Origin: Large & Yeager 2008, Eq. (9) and (12) 
     243 
     244      !!---------------------------------------------------------------------------------- 
     245      REAL(wp), DIMENSION(jpi,jpj)             :: ch_n10_ncar 
     246      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 ) 
     247      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pstab      ! stable ABL => 1 / unstable ABL => 0 
     248      !!---------------------------------------------------------------------------------- 
     249      ! 
     250      ch_n10_ncar = MAX( 1.e-3_wp * psqrtcdn10*( 18._wp*pstab + 32.7_wp*(1._wp - pstab) )  , 0.1E-3_wp )   ! Eq. (9) & (12) Large & Yeager, 2008 
     251      ! 
     252   END FUNCTION CH_N10_NCAR 
     253 
     254   FUNCTION CE_N10_NCAR( psqrtcdn10 ) 
     255      !!---------------------------------------------------------------------------------- 
     256      !! Estimate of the neutral heat transfer coefficient at 10m      !! 
     257      !! Origin: Large & Yeager 2008, Eq. (9) and (13) 
     258      !!---------------------------------------------------------------------------------- 
     259      REAL(wp), DIMENSION(jpi,jpj)             :: ce_n10_ncar 
     260      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 ) 
     261      !!---------------------------------------------------------------------------------- 
     262      ce_n10_ncar = MAX( 1.e-3_wp * ( 34.6_wp * psqrtcdn10 ) , 0.1E-3_wp ) 
     263      ! 
     264   END FUNCTION CE_N10_NCAR 
     265 
     266 
     267   FUNCTION psi_m_ncar( pzeta ) 
    272268      !!---------------------------------------------------------------------------------- 
    273269      !! Universal profile stability function for momentum 
    274       !!    !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 
    275       !!      
    276       !! pzet0 : stability paramenter, z/L where z is altitude measurement                                           
     270      !!    !! Psis, L&Y 2004, Eq. (8c), (8d), (8e) 
     271      !! 
     272      !! pzeta : stability paramenter, z/L where z is altitude measurement 
    277273      !!         and L is M-O length 
    278274      !! 
    279       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    280       !!---------------------------------------------------------------------------------- 
    281       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pzeta 
    282       REAL(wp), DIMENSION(jpi,jpj)             ::   psi_m 
    283       ! 
    284       INTEGER  ::   ji, jj         ! dummy loop indices 
    285       REAL(wp) :: zx2, zx, zstab   ! local scalars 
    286       !!---------------------------------------------------------------------------------- 
    287       ! 
     275      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     276      !!---------------------------------------------------------------------------------- 
     277      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ncar 
     278      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 
     279      ! 
     280      INTEGER  ::   ji, jj    ! dummy loop indices 
     281      REAL(wp) :: zzeta, zx2, zx, zpsi_unst, zpsi_stab,  zstab   ! local scalars 
     282      !!---------------------------------------------------------------------------------- 
    288283      DO jj = 1, jpj 
    289284         DO ji = 1, jpi 
    290             zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) ) 
    291             zx2 = MAX ( zx2 , 1. ) 
    292             zx  = SQRT( zx2 ) 
    293             zstab = 0.5 + SIGN( 0.5 , pzeta(ji,jj) ) 
    294             ! 
    295             psi_m(ji,jj) =        zstab  * (-5.*pzeta(ji,jj))       &          ! Stable 
    296                &          + (1. - zstab) * (2.*LOG((1. + zx)*0.5)   &          ! Unstable 
    297                &               + LOG((1. + zx2)*0.5) - 2.*ATAN(zx) + rpi*0.5)  !    " 
     285 
     286            zzeta = pzeta(ji,jj) 
     287            ! 
     288            zx2 = SQRT( ABS(1._wp - 16._wp*zzeta) )  ! (1 - 16z)^0.5 
     289            zx2 = MAX( zx2 , 1._wp ) 
     290            zx  = SQRT(zx2)                          ! (1 - 16z)^0.25 
     291            zpsi_unst = 2._wp*LOG( (1._wp + zx )*0.5_wp )   & 
     292               &            + LOG( (1._wp + zx2)*0.5_wp )   & 
     293               &          - 2._wp*ATAN(zx) + rpi*0.5_wp 
     294            ! 
     295            zpsi_stab = -5._wp*zzeta 
     296            ! 
     297            zstab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => zstab = 1 
     298            ! 
     299            psi_m_ncar(ji,jj) =          zstab  * zpsi_stab &  ! (zzeta > 0) Stable 
     300               &              + (1._wp - zstab) * zpsi_unst    ! (zzeta < 0) Unstable 
    298301            ! 
    299302         END DO 
    300303      END DO 
    301       ! 
    302    END FUNCTION psi_m 
    303  
    304  
    305    FUNCTION psi_h( pzeta ) 
     304   END FUNCTION psi_m_ncar 
     305 
     306 
     307   FUNCTION psi_h_ncar( pzeta ) 
    306308      !!---------------------------------------------------------------------------------- 
    307309      !! Universal profile stability function for temperature and humidity 
    308       !!    !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 
    309       !! 
    310       !! pzet0 : stability paramenter, z/L where z is altitude measurement                                           
     310      !!    !! Psis, L&Y 2004, Eq. (8c), (8d), (8e) 
     311      !! 
     312      !! pzeta : stability paramenter, z/L where z is altitude measurement 
    311313      !!         and L is M-O length 
    312314      !! 
    313       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    314       !!---------------------------------------------------------------------------------- 
     315      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     316      !!---------------------------------------------------------------------------------- 
     317      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ncar 
    315318      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 
    316       REAL(wp), DIMENSION(jpi,jpj)             :: psi_h 
    317       ! 
    318       INTEGER  ::   ji, jj    ! dummy loop indices 
    319       REAL(wp) :: zx2, zstab  ! local scalars 
     319      ! 
     320      INTEGER  ::   ji, jj     ! dummy loop indices 
     321      REAL(wp) :: zzeta, zx2, zpsi_unst, zpsi_stab, zstab  ! local scalars 
    320322      !!---------------------------------------------------------------------------------- 
    321323      ! 
    322324      DO jj = 1, jpj 
    323325         DO ji = 1, jpi 
    324             zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) ) 
    325             zx2 = MAX ( zx2 , 1. ) 
    326             zstab = 0.5 + SIGN( 0.5 , pzeta(ji,jj) ) 
    327             ! 
    328             psi_h(ji,jj) =         zstab  * (-5.*pzeta(ji,jj))        &  ! Stable 
    329                &           + (1. - zstab) * (2.*LOG( (1. + zx2)*0.5 ))   ! Unstable 
     326            ! 
     327            zzeta = pzeta(ji,jj) 
     328            ! 
     329            zx2 = SQRT( ABS(1._wp - 16._wp*zzeta) )  ! (1 -16z)^0.5 
     330            zx2 = MAX( zx2 , 1._wp ) 
     331            zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) ) 
     332            ! 
     333            zpsi_stab = -5._wp*zzeta 
     334            ! 
     335            zstab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => zstab = 1 
     336            ! 
     337            psi_h_ncar(ji,jj) =          zstab  * zpsi_stab &  ! (zzeta > 0) Stable 
     338               &              + (1._wp - zstab) * zpsi_unst    ! (zzeta < 0) Unstable 
    330339            ! 
    331340         END DO 
    332341      END DO 
    333       ! 
    334    END FUNCTION psi_h 
     342   END FUNCTION psi_h_ncar 
     343 
     344 
     345 
     346 
     347   FUNCTION UN10_from_CD( pzu, pUb, pCd, ppsi ) 
     348      !!---------------------------------------------------------------------------------- 
     349      !!  Provides the neutral-stability wind speed at 10 m 
     350      !!---------------------------------------------------------------------------------- 
     351      REAL(wp), DIMENSION(jpi,jpj)             :: UN10_from_CD  !: [m/s] 
     352      REAL(wp),                     INTENT(in) :: pzu  !: measurement heigh of bulk wind speed 
     353      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb  !: bulk wind speed at height pzu m   [m/s] 
     354      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd  !: drag coefficient 
     355      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppsi !: "Psi_m(pzu/L)" stability correction profile for momentum [] 
     356      !!---------------------------------------------------------------------------------- 
     357      !! Reminder: UN10 = u*/vkarmn * log(10/z0) 
     358      !!     and: u* = sqrt(Cd) * Ub 
     359      !!                                  u*/vkarmn * log(   10   /       z0    ) 
     360      UN10_from_CD(:,:) = SQRT(pCd(:,:))*pUb/vkarmn * LOG( 10._wp / z0_from_Cd( pzu, pCd(:,:), ppsi=ppsi(:,:) ) ) 
     361      !! 
     362   END FUNCTION UN10_from_CD 
     363 
     364 
     365   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs ) 
     366      !!------------------------------------------------------------------------ 
     367      !! 
     368      !! Evaluates the 1./(Obukhov length) from air temperature, 
     369      !! air specific humidity, and frictional scales u*, t* and q* 
     370      !! 
     371      !! Author: L. Brodeau, June 2019 / AeroBulk 
     372      !!         (https://github.com/brodeau/aerobulk/) 
     373      !!------------------------------------------------------------------------ 
     374      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L     !: 1./(Obukhov length) [m^-1] 
     375      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha         !: reference potential temperature of air [K] 
     376      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa          !: reference specific humidity of air   [kg/kg] 
     377      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus          !: u*: friction velocity [m/s] 
     378      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pts, pqs     !: \theta* and q* friction aka turb. scales for temp. and spec. hum. 
     379      ! 
     380      INTEGER  ::   ji, jj         ! dummy loop indices 
     381      REAL(wp) ::     zqa          ! local scalar 
     382      !!------------------------------------------------------------------- 
     383      ! 
     384      DO jj = 1, jpj 
     385         DO ji = 1, jpi 
     386            ! 
     387            zqa = (1._wp + rctv0*pqa(ji,jj)) 
     388            ! 
     389            ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with: 
     390            !  a/  -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one! 
     391            !                      or 
     392            !  b/  -u* [ theta*              + 0.61 theta q* ] 
     393            ! 
     394            One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) & 
     395               &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp ) 
     396            ! 
     397         END DO 
     398      END DO 
     399      ! 
     400      One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...) 
     401      ! 
     402   END FUNCTION One_on_L 
     403 
     404 
     405   FUNCTION z0_from_Cd( pzu, pCd,  ppsi ) 
     406      REAL(wp), DIMENSION(jpi,jpj) :: z0_from_Cd        !: roughness length [m] 
     407      REAL(wp)                    , INTENT(in) :: pzu   !: reference height zu [m] 
     408      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd   !: (neutral or non-neutral) drag coefficient [] 
     409      REAL(wp), DIMENSION(jpi,jpj), INTENT(in), OPTIONAL :: ppsi !: "Psi_m(pzu/L)" stability correction profile for momentum [] 
     410      !! 
     411      !! If pCd is the NEUTRAL-STABILITY drag coefficient then ppsi must be 0 or not given 
     412      !! If pCd is the drag coefficient (in stable or unstable conditions) then pssi must be provided 
     413      !!---------------------------------------------------------------------------------- 
     414      IF ( PRESENT(ppsi) ) THEN 
     415         !! Cd provided is the actual Cd (not the neutral-stability CdN) : 
     416         z0_from_Cd = pzu * EXP( - ( vkarmn/SQRT(pCd(:,:)) + ppsi(:,:) ) ) !LB: ok, double-checked! 
     417      ELSE 
     418         !! Cd provided is the neutral-stability Cd, aka CdN : 
     419         z0_from_Cd = pzu * EXP( - vkarmn/SQRT(pCd(:,:)) )            !LB: ok, double-checked! 
     420      END IF 
     421   END FUNCTION z0_from_Cd 
     422 
     423   FUNCTION virt_temp( pta, pqa ) 
     424      REAL(wp), DIMENSION(jpi,jpj)             :: virt_temp !: virtual temperature [K] 
     425      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute or potential air temperature [K] 
     426      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa !: specific humidity of air   [kg/kg] 
     427      virt_temp(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:)) 
     428   END FUNCTION virt_temp 
    335429 
    336430   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.