Ignore:
Timestamp:
2019-11-29T16:59:07+01:00 (3 months ago)
Author:
gsamson
Message:

dev_ASINTER-01-05_merged: merge dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk branch @rev11988 with dev_r11265_ASINTER-01_Guillaume_ABL1D branch @rev11937 (tickets #2159 and #2131); ORCA2_ICE(_ABL) reproductibility OK

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk_algo_ncar.F90

    r10190 r12015  
    1111   !! 
    1212   !!       Routine turb_ncar maintained and developed in AeroBulk 
    13    !!                     (http://aerobulk.sourceforge.net/) 
     13   !!                     (https://github.com/brodeau/aerobulk/) 
    1414   !! 
    1515   !!                         L. Brodeau, 2015 
     
    2626   USE dom_oce         ! ocean space and time domain 
    2727   USE phycst          ! physical constants 
    28    USE sbc_oce         ! Surface boundary condition: ocean fields 
     28   USE iom             ! I/O manager library 
     29   USE lib_mpp         ! distribued memory computing library 
     30   USE in_out_manager  ! I/O manager 
     31   USE prtctl          ! Print control 
    2932   USE sbcwave, ONLY   :  cdn_wave ! wave module 
    3033#if defined key_si3 || defined key_cice 
    3134   USE sbc_ice         ! Surface boundary condition: ice fields 
    3235#endif 
    33    ! 
    34    USE iom             ! I/O manager library 
    35    USE lib_mpp         ! distribued memory computing library 
    36    USE in_out_manager  ! I/O manager 
    37    USE prtctl          ! Print control 
    3836   USE lib_fortran     ! to use key_nosignedzero 
    3937 
     38   USE sbc_oce         ! Surface boundary condition: ocean fields 
     39   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB 
    4040 
    4141   IMPLICIT NONE 
    4242   PRIVATE 
    4343 
    44    PUBLIC ::   TURB_NCAR   ! called by sbcblk.F90 
    45  
    46    !                              ! NCAR own values for given constants: 
    47    REAL(wp), PARAMETER ::   rctv0 = 0.608   ! constant to obtain virtual temperature... 
    48     
     44   PUBLIC :: TURB_NCAR   ! called by sbcblk.F90 
     45 
     46   INTEGER , PARAMETER ::   nb_itt = 5        ! number of itterations 
     47 
    4948   !!---------------------------------------------------------------------- 
    5049CONTAINS 
     
    5352      &                  Cd, Ch, Ce, t_zu, q_zu, U_blk,      & 
    5453      &                  Cdn, Chn, Cen                       ) 
    55       !!---------------------------------------------------------------------------------- 
     54      !!---------------------------------------------------------------------- 
    5655      !!                      ***  ROUTINE  turb_ncar  *** 
    5756      !! 
     
    6160      !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas 
    6261      !! 
    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) 
    7962      !! 
    8063      !! INPUT : 
    8164      !! ------- 
    8265      !!    *  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] 
     66      !!    *  zu   : height for wind speed (usually 10m)                     [m] 
     67      !!    *  sst  : bulk SST                                                [K] 
    8668      !!    *  t_zt : potential air temperature at zt                         [K] 
    8769      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg] 
    8870      !!    *  q_zt : specific humidity of air at zt                          [kg/kg] 
     71      !!    *  U_zu : scalar wind speed at zu                                 [m/s] 
    8972      !! 
    9073      !! 
     
    9679      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K] 
    9780      !!    *  q_zu   : specific humidity of air        //                    [kg/kg] 
    98       !!    *  U_blk  : bulk wind at 10m                                      [m/s] 
     81      !!    *  U_blk  : bulk wind speed at zu                                 [m/s] 
     82      !! 
     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] 
     98      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind speed at zu                     [m/s] 
    11399      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    114100      ! 
    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 
     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 
     
    126111      ! 
    127112      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 
     113      IF( ABS(zu - zt) < 0.01_wp )  l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
     114 
     115      U_blk = MAX( 0.5_wp , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 
    131116 
    132117      !! 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 
     118      ztmp0 = virt_temp(t_zt, q_zt) - virt_temp(sst, ssq) ! air-sea difference of virtual pot. temp. at zt 
     119      stab  = 0.5_wp + sign(0.5_wp,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable 
    135120 
    136121      !! Neutral coefficients at 10m: 
     
    139124         ztmp0   (:,:) = cdn_wave(:,:) 
    140125      ELSE 
    141          ztmp0 = cd_neutral_10m( U_blk ) 
     126      ztmp0 = cd_neutral_10m( U_blk ) 
    142127      ENDIF 
    143128 
     
    146131      !! Initializing transf. coeff. with their first guess neutral equivalents : 
    147132      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)) 
     133      Ce = 1.e-3_wp*( 34.6_wp * sqrt_Cd_n10 ) 
     134      Ch = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab)) 
    150135      stab = sqrt_Cd_n10   ! Temporaty array !!! stab == SQRT(Cd) 
    151136  
    152       IF( ln_cdgw )   Cen = Ce  ; Chn = Ch 
     137      IF( ln_cdgw ) THEN 
     138   Cen = Ce 
     139   Chn = Ch 
     140      ENDIF 
    153141 
    154142      !! Initializing values at z_u with z_t values: 
    155143      t_zu = t_zt   ;   q_zu = q_zt 
    156144 
    157       !!  * Now starting iteration loop 
    158       DO j_itt=1, nb_itt 
     145      !! ITERATION BLOCK 
     146      DO j_itt = 1, nb_itt 
    159147         ! 
    160148         ztmp1 = t_zu - sst   ! Updating air/sea differences 
     
    162150 
    163151         ! 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 
     152         ztmp0 = stab*U_blk       ! u*       (stab == SQRT(Cd)) 
     153         ztmp1 = Ch/stab*ztmp1    ! theta*   (stab == SQRT(Cd)) 
     154         ztmp2 = Ce/stab*ztmp2    ! q*       (stab == SQRT(Cd)) 
    168155 
    169156         ! 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 ) 
    172  
     157         ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 ) 
     158          
    173159         !! Stability parameters : 
    174          zeta_u   = zu*ztmp0   ;  zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 
     160         zeta_u   = zu*ztmp0 
     161         zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u ) 
    175162         zpsi_h_u = psi_h( zeta_u ) 
    176163 
     
    178165         IF( .NOT. l_zt_equal_zu ) THEN 
    179166            !! 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 !!! 
     167            stab = zt*ztmp0 
     168            stab = SIGN( MIN(ABS(stab),10._wp), stab )  ! Temporaty array stab == zeta_t !!! 
    181169            stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab)                   ! stab just used as temp array again! 
    182170            t_zu = t_zt - ztmp1/vkarmn*stab    ! ztmp1 is still theta*  L&Y 2004 eq.(9b) 
    183171            q_zu = q_zt - ztmp2/vkarmn*stab    ! ztmp2 is still q*      L&Y 2004 eq.(9c) 
    184             q_zu = max(0., q_zu) 
    185          END IF 
    186  
     172            q_zu = max(0._wp, q_zu) 
     173         ENDIF 
     174 
     175         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 
     176         !   In very rare low-wind conditions, the old way of estimating the 
     177         !   neutral wind speed at 10m leads to a negative value that causes the code 
     178         !   to crash. To prevent this a threshold of 0.25m/s is imposed. 
    187179         ztmp2 = psi_m(zeta_u) 
    188180         IF( ln_cdgw ) THEN      ! surface wave case 
    189181            stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 )  ! (stab == SQRT(Cd)) 
    190182            Cd   = stab * stab 
    191             ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
     183            ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
    192184            ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd)) 
    193             ztmp1 = 1. + Chn * ztmp0      
     185            ztmp1 = 1._wp + Chn * ztmp0      
    194186            Ch    = Chn * ztmp2 / ztmp1  ! L&Y 2004 eq. (10b) 
    195             ztmp1 = 1. + Cen * ztmp0 
     187            ztmp1 = 1._wp + Cen * ztmp0 
    196188            Ce    = Cen * ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
    197189 
    198190         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       ! 
     191         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 
     192         !   In very rare low-wind conditions, the old way of estimating the 
     193         !   neutral wind speed at 10m leads to a negative value that causes the code 
     194         !   to crash. To prevent this a threshold of 0.25m/s is imposed. 
     195         ztmp0 = MAX( 0.25_wp , U_blk/(1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u)) 
     196         ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10 
     197         Cdn(:,:) = ztmp0 
     198         sqrt_Cd_n10 = sqrt(ztmp0) 
     199 
     200         stab    = 0.5_wp + sign(0.5_wp,zeta_u)                        ! update stability 
     201         Cx_n10  = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10) 
     202         Chn(:,:) = Cx_n10 
     203 
     204         !! Update of transfer coefficients: 
     205         ztmp1 = 1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u)) 
     206         Cd      = ztmp0 / ( ztmp1*ztmp1 ) 
     207         stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 
     208 
     209         ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
     210         ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd)) 
     211         ztmp1 = 1._wp + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10) 
     212         Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b) 
     213 
     214         Cx_n10  = 1.e-3_wp * (34.6_wp * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10 
     215         Cen(:,:) = Cx_n10 
     216         ztmp1 = 1._wp + Cx_n10*ztmp0 
     217         Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
     218         ENDIF 
     219 
     220      END DO !DO j_itt = 1, nb_itt 
     221 
    230222   END SUBROUTINE turb_ncar 
    231223 
     
    238230      !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b) 
    239231      !! 
    240       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     232      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    241233      !!---------------------------------------------------------------------------------- 
    242234      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10           ! scalar wind speed at 10m (m/s) 
     
    255247            ! 
    256248            ! 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) 
     249            zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) )   ! If pw10 < 33. => 0, else => 1 
     250            ! 
     251            cd_neutral_10m(ji,jj) = 1.e-3_wp * ( & 
     252               &       (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind <  33 m/s 
     253               &      +    zgt33   *      2.34_wp )                                                 ! wind >= 33 m/s 
     254            ! 
     255            cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp) 
    264256            ! 
    265257         END DO 
     
    273265      !! Universal profile stability function for momentum 
    274266      !!    !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 
    275       !!      
    276       !! pzet0 : stability paramenter, z/L where z is altitude measurement                                           
     267      !! 
     268      !! pzeta : stability paramenter, z/L where z is altitude measurement 
    277269      !!         and L is M-O length 
    278270      !! 
    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 
     271      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     272      !!---------------------------------------------------------------------------------- 
     273      REAL(wp), DIMENSION(jpi,jpj) :: psi_m 
     274      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 
     275      ! 
     276      INTEGER  ::   ji, jj    ! dummy loop indices 
    285277      REAL(wp) :: zx2, zx, zstab   ! local scalars 
    286278      !!---------------------------------------------------------------------------------- 
    287       ! 
    288279      DO jj = 1, jpj 
    289280         DO ji = 1, jpi 
    290             zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) ) 
    291             zx2 = MAX ( zx2 , 1. ) 
     281            zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 
     282            zx2 = MAX( zx2 , 1._wp ) 
    292283            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)  !    " 
     284            zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 
     285            ! 
     286            psi_m(ji,jj) =        zstab  * (-5._wp*pzeta(ji,jj))       &          ! Stable 
     287               &          + (1._wp - zstab) * (2._wp*LOG((1._wp + zx)*0.5_wp)   &          ! Unstable 
     288               &               + LOG((1._wp + zx2)*0.5_wp) - 2._wp*ATAN(zx) + rpi*0.5_wp)  !    " 
    298289            ! 
    299290         END DO 
    300291      END DO 
    301       ! 
    302292   END FUNCTION psi_m 
    303293 
     
    308298      !!    !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 
    309299      !! 
    310       !! pzet0 : stability paramenter, z/L where z is altitude measurement                                           
     300      !! pzeta : stability paramenter, z/L where z is altitude measurement 
    311301      !!         and L is M-O length 
    312302      !! 
    313       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    314       !!---------------------------------------------------------------------------------- 
     303      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     304      !!---------------------------------------------------------------------------------- 
     305      REAL(wp), DIMENSION(jpi,jpj) :: psi_h 
    315306      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 
    316       REAL(wp), DIMENSION(jpi,jpj)             :: psi_h 
    317       ! 
    318       INTEGER  ::   ji, jj    ! dummy loop indices 
     307      ! 
     308      INTEGER  ::   ji, jj     ! dummy loop indices 
    319309      REAL(wp) :: zx2, zstab  ! local scalars 
    320310      !!---------------------------------------------------------------------------------- 
     
    322312      DO jj = 1, jpj 
    323313         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 
     314            zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 
     315            zx2 = MAX( zx2 , 1._wp ) 
     316            zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 
     317            ! 
     318            psi_h(ji,jj) =         zstab  * (-5._wp*pzeta(ji,jj))        &  ! Stable 
     319               &           + (1._wp - zstab) * (2._wp*LOG( (1._wp + zx2)*0.5_wp ))   ! Unstable 
    330320            ! 
    331321         END DO 
    332322      END DO 
    333       ! 
    334323   END FUNCTION psi_h 
    335324 
Note: See TracChangeset for help on using the changeset viewer.