Changeset 11209


Ignore:
Timestamp:
2019-07-03T09:09:17+02:00 (13 months ago)
Author:
laurent
Message:

LB: making more efficient and more centralized use of physical/thermodynamics functions defined into "OCE/SBC/sbcblk_phy.F90".

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

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90

    r11182 r11209  
    3131   USE sbc_ice         ! Surface boundary condition: ice fields 
    3232#endif 
    33    USE sbcblk_phy     !LB: all thermodynamics functions, rho_air, q_sat, etc... 
     33   USE sbcblk_phy     !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB 
    3434   USE in_out_manager ! I/O manager 
    3535   USE iom            ! I/O manager library 
     
    149149      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 
    150150 
    151       ztmp2  = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk)  !! Ribu Bulk Richardson number ;       !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth 
    152  
    153       !! First estimate of zeta_u, depending on the stability, ie sign of Ribu (ztmp2): 
     151      ztmp2 = Ri_bulk( zu, sst, t_zu, ssq, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 
     152      !ztmp2  = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk)  !! Bulk Richardson number ;       !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth 
     153 
     154      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 
    154155      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 
    155156      ztmp0 = ztmp0*ztmp2 
    156       zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & !  Ribu < 0 
    157          &  +     ztmp1   * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))               !  Ribu > 0 
     157      zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & !  BRN < 0 
     158         &  +     ztmp1   * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))               !  BRN > 0 
    158159      !#LOLO: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 
    159160 
     
    285286   END FUNCTION alfa_charn 
    286287 
    287  
    288    FUNCTION One_on_L( ptha, pqa, pus, pts, pqs ) 
    289       !!------------------------------------------------------------------------ 
    290       !! 
    291       !! Evaluates the 1./(Monin Obukhov length) from air temperature and 
    292       !!  specific humidity, and frictional scales u*, t* and q* 
    293       !! 
    294       !! Author: L. Brodeau, june 2016 / AeroBulk 
    295       !!         (https://github.com/brodeau/aerobulk/) 
    296       !!------------------------------------------------------------------------ 
    297       REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1] 
    298       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha,  &  !: average potetntial air temperature [K] 
    299          &                                        pqa,   &  !: average specific humidity of air   [kg/kg] 
    300          &                                      pus, pts, pqs   !: frictional velocity, temperature and humidity 
    301       ! 
    302       INTEGER  ::   ji, jj         ! dummy loop indices 
    303       REAL(wp) ::     zqa          ! local scalar 
    304       !!------------------------------------------------------------------- 
    305       ! 
    306       DO jj = 1, jpj 
    307          DO ji = 1, jpi 
    308             ! 
    309             zqa = (1. + rctv0*pqa(ji,jj)) 
    310             ! 
    311             One_on_L(ji,jj) =  grav*vkarmn*(pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj)) & 
    312                &                      / ( pus(ji,jj)*pus(ji,jj) * ptha(ji,jj)*zqa ) 
    313             ! 
    314          END DO 
    315       END DO 
    316       ! 
    317    END FUNCTION One_on_L 
    318  
    319  
    320288   FUNCTION psi_m_coare( pzeta ) 
    321289      !!---------------------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p5.F90

    r11111 r11209  
    3535   USE sbc_ice         ! Surface boundary condition: ice fields 
    3636#endif 
    37    ! 
     37   USE sbcblk_phy     !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB 
    3838   USE iom             ! I/O manager library 
    3939   USE lib_mpp         ! distribued memory computing library 
     
    255255   END SUBROUTINE turb_coare3p5 
    256256 
    257  
    258  
    259    FUNCTION One_on_L( ptha, pqa, pus, pts, pqs ) 
    260       !!------------------------------------------------------------------------ 
    261       !! 
    262       !! Evaluates the 1./(Monin Obukhov length) from air temperature and 
    263       !!  specific humidity, and frictional scales u*, t* and q* 
    264       !! 
    265       !! Author: L. Brodeau, june 2016 / AeroBulk 
    266       !!         (https://github.com/brodeau/aerobulk/) 
    267       !!------------------------------------------------------------------------ 
    268       REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1] 
    269       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha,  &  !: average potetntial air temperature [K] 
    270          &                                        pqa,   &  !: average specific humidity of air   [kg/kg] 
    271          &                                      pus, pts, pqs   !: frictional velocity, temperature and humidity 
    272       ! 
    273       INTEGER  ::   ji, jj         ! dummy loop indices 
    274       REAL(wp) ::     zqa          ! local scalar 
    275       ! 
    276       DO jj = 1, jpj 
    277          DO ji = 1, jpi 
    278             ! 
    279             zqa = (1. + rctv0*pqa(ji,jj)) 
    280             ! 
    281             One_on_L(ji,jj) =  grav*vkarmn*(pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj)) & 
    282                &                      / ( pus(ji,jj)*pus(ji,jj) * ptha(ji,jj)*zqa ) 
    283             ! 
    284          END DO 
    285       END DO 
    286       ! 
    287    END FUNCTION One_on_L 
    288  
    289  
    290257   FUNCTION psi_m_coare( pzeta ) 
    291258      !!---------------------------------------------------------------------------------- 
     
    389356   END FUNCTION psi_h_coare 
    390357 
    391  
    392    FUNCTION visc_air( ptak ) 
    393       !!--------------------------------------------------------------------- 
    394       !! Air kinetic viscosity (m^2/s) given from temperature in degrees... 
    395       !! 
    396       !! Author: L. Brodeau, june 2016 / AeroBulk 
    397       !!         (https://github.com/brodeau/aerobulk/) 
    398       !!--------------------------------------------------------------------- 
    399       REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air 
    400       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature   [K] 
    401       ! 
    402       INTEGER  ::   ji, jj         ! dummy loop indices 
    403       REAL(wp) ::   ztc, ztc2      ! local scalar 
    404       ! 
    405       DO jj = 1, jpj 
    406          DO ji = 1, jpi 
    407             ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C 
    408             ztc2 = ztc*ztc 
    409             visc_air(ji,jj) = 1.326E-5*(1. + 6.542E-3*ztc + 8.301E-6*ztc2 - 4.84E-9*ztc2*ztc) 
    410          END DO 
    411       END DO 
    412       ! 
    413    END FUNCTION visc_air 
    414  
    415358   !!====================================================================== 
    416359END MODULE sbcblk_algo_coare3p5 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90

    r11182 r11209  
    3838 
    3939   USE sbc_oce         ! Surface boundary condition: ocean fields 
    40    USE sbcblk_phy   !LB: all thermodynamics functions, rho_air, q_sat, etc... 
     40   USE sbcblk_phy      !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB 
    4141   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) 
    4242 
     
    7272   !                            !:    => ACCEPTABLE IN MOST CONDITIONS ! (UNLESS: sunny + very calm/low-wind conditions) 
    7373   !                            !:  => Otherwize use "nb_itt_wl = 10" 
    74  
    75  
    76  
    77  
    7874   !!---------------------------------------------------------------------- 
    7975CONTAINS 
     
    213209      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 
    214210 
    215       ztmp2 = Ri_bulk( zu, t_zu, dt_zu, q_zu, dq_zu, U_blk )   ! Ribu = Bulk Richardson number 
    216  
    217       !! First estimate of zeta_u, depending on the stability, ie sign of Ribu (ztmp2): 
     211      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 
     212 
     213      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 
    218214      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 
    219215      func_m = ztmp0*ztmp2 ! temporary array !! 
    220       func_h = (1.-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  Ribu < 0 ! temporary array !!! func_h == zeta_u 
    221          &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))              !  Ribu > 0 
     216      func_h = (1.-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0 ! temporary array !!! func_h == zeta_u 
     217         &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))           !  BRN > 0 
    222218      !#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" ! 
    223219 
     
    247243 
    248244      !! First guess of inverse of Monin-Obukov length (1/L) : 
    249       ztmp0 = (1. + rctv0*q_zu)  ! the factor to apply to temp. to get virt. temp... 
    250       Linv  =  grav*vkarmn*(t_star*ztmp0 + rctv0*t_zu*q_star) / MAX( u_star*u_star * t_zu*ztmp0 , 1.E-9 ) ! #LOLO 
    251       Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
     245      Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star ) 
    252246 
    253247      !! Functions such as  u* = U_blk*vkarmn/func_m 
     
    260254 
    261255         !! Bulk Richardson Number at z=zu (Eq. 3.25) 
    262          ztmp0 = Ri_bulk(zu, t_zu, dt_zu, q_zu, dq_zu, U_blk) 
     256         ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 
    263257 
    264258         !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) : 
     
    278272         z0t    = MIN( ABS( alpha_H*ztmp1                     ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 
    279273         z0q    = MIN( ABS( alpha_Q*ztmp1                     ) , 0.001_wp) 
    280           
     274 
    281275         !! Update wind at 10m taking into acount convection-related wind gustiness: 
    282276         !! => Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8 
     
    447441 
    448442 
    449    FUNCTION Ri_bulk( pz, ptz, pdt, pqz, pdq, pub ) 
    450       !!---------------------------------------------------------------------------------- 
    451       !! Bulk Richardson number (Eq. 3.25 IFS doc) 
    452       !! 
    453       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    454       !!---------------------------------------------------------------------------------- 
    455       REAL(wp), DIMENSION(jpi,jpj) ::   Ri_bulk   ! 
    456       ! 
    457       REAL(wp)                    , INTENT(in) ::   pz    ! height above the sea        [m] 
    458       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptz   ! air temperature at pz m     [K] 
    459       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pdt   ! ptz - sst                   [K] 
    460       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqz   ! air temperature at pz m [kg/kg] 
    461       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pdq   ! pqz - ssq               [kg/kg] 
    462       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pub   ! bulk wind speed           [m/s] 
    463       !!---------------------------------------------------------------------------------- 
    464       ! 
    465       Ri_bulk =   grav*pz/(pub*pub)   & 
    466          &      * ( pdt/(ptz - 0.5_wp*(pdt + grav*pz/(rCp_dry + rCp_vap*pqz))) & 
    467          &          + rctv0*pdq ) 
    468       ! 
    469    END FUNCTION Ri_bulk 
    470443 
    471444 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90

    r11111 r11209  
    3838   USE lib_fortran     ! to use key_nosignedzero 
    3939 
     40   USE sbcblk_phy      !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB 
    4041 
    4142   IMPLICIT NONE 
     
    9394      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K] 
    9495      !!    *  q_zu   : specific humidity of air        //                    [kg/kg] 
    95       !!    *  U_blk  : bulk wind at 10m                                      [m/s] 
     96      !!    *  U_blk  : bulk wind speed at 10m                                [m/s] 
    9697      !!---------------------------------------------------------------------------------- 
    9798      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m] 
     
    125126      IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
    126127 
    127       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 
     128      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 
    128129 
    129130      !! First guess of stability: 
    130       ztmp0 = t_zt*(1. + rctv0*q_zt) - sst*(1. + rctv0*ssq) ! air-sea difference of virtual pot. temp. at zt 
    131       stab  = 0.5 + sign(0.5,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable 
     131      ztmp0 = virt_temp(t_zt, q_zt) - virt_temp(sst, ssq) ! air-sea difference of virtual pot. temp. at zt 
     132      stab  = 0.5_wp + sign(0.5_wp,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable 
    132133 
    133134      !! Neutral coefficients at 10m: 
     
    159160 
    160161         ! Updating turbulent scales :   (L&Y 2004 eq. (7)) 
    161          ztmp1  = Ch/stab*ztmp1    ! theta*   (stab == SQRT(Cd)) 
    162          ztmp2  = Ce/stab*ztmp2    ! q*       (stab == SQRT(Cd)) 
    163  
    164          ztmp0 = 1. + rctv0*q_zu      ! multiply this with t and you have the virtual temperature 
     162         ztmp0 = stab*U_blk       ! u*       (stab == SQRT(Cd)) 
     163         ztmp1 = Ch/stab*ztmp1    ! theta*   (stab == SQRT(Cd)) 
     164         ztmp2 = Ce/stab*ztmp2    ! q*       (stab == SQRT(Cd)) 
    165165 
    166166         ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 
    167          ztmp0 =  (grav*vkarmn/(t_zu*ztmp0)*(ztmp1*ztmp0 + rctv0*t_zu*ztmp2)) / (Cd*U_blk*U_blk) 
    168          !                                                      ( Cd*U_blk*U_blk is U*^2 at zu ) 
    169  
     167         ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 ) 
     168          
    170169         !! Stability parameters : 
    171          zeta_u   = zu*ztmp0   ;  zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 
     170         zeta_u   = zu*ztmp0 
     171         zeta_u = sign( min(abs(zeta_u),10.0_wp), zeta_u ) 
    172172         zpsi_h_u = psi_h( zeta_u ) 
    173173 
     
    175175         IF( .NOT. l_zt_equal_zu ) THEN 
    176176            !! Array 'stab' is free for the moment so using it to store 'zeta_t' 
    177             stab = zt*ztmp0 ;  stab = SIGN( MIN(ABS(stab),10.0), stab )  ! Temporaty array stab == zeta_t !!! 
     177            stab = zt*ztmp0 
     178            stab = SIGN( MIN(ABS(stab),10.0_wp), stab )  ! Temporaty array stab == zeta_t !!! 
    178179            stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab)                   ! stab just used as temp array again! 
    179180            t_zu = t_zt - ztmp1/vkarmn*stab    ! ztmp1 is still theta*  L&Y 2004 eq.(9b) 
    180181            q_zu = q_zt - ztmp2/vkarmn*stab    ! ztmp2 is still q*      L&Y 2004 eq.(9c) 
    181             q_zu = max(0., q_zu) 
     182            q_zu = max(0._wp, q_zu) 
    182183         END IF 
    183184 
     
    194195 
    195196         ELSE 
    196             ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 
    197             !   In very rare low-wind conditions, the old way of estimating the 
    198             !   neutral wind speed at 10m leads to a negative value that causes the code 
    199             !   to crash. To prevent this a threshold of 0.25m/s is imposed. 
    200             ztmp0 = MAX( 0.25 , U_blk/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u)) 
    201             ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10 
    202             Cdn(:,:) = ztmp0 
    203             sqrt_Cd_n10 = sqrt(ztmp0) 
    204  
    205             stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability 
    206             Cx_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10) 
    207             Chn(:,:) = Cx_n10 
    208  
    209             !! Update of transfer coefficients: 
    210             ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u)) 
    211             Cd      = ztmp0 / ( ztmp1*ztmp1 ) 
    212             stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 
    213  
    214             ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
    215             ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd)) 
    216             ztmp1 = 1. + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10) 
    217             Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b) 
    218  
    219             Cx_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10 
    220             Cen(:,:) = Cx_n10 
    221             ztmp1 = 1. + Cx_n10*ztmp0 
    222             Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
    223             ENDIF 
     197         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 
     198         !   In very rare low-wind conditions, the old way of estimating the 
     199         !   neutral wind speed at 10m leads to a negative value that causes the code 
     200         !   to crash. To prevent this a threshold of 0.25m/s is imposed. 
     201         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)) 
     202         ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10 
     203         Cdn(:,:) = ztmp0 
     204         sqrt_Cd_n10 = sqrt(ztmp0) 
     205 
     206         stab    = 0.5_wp + sign(0.5_wp,zeta_u)                        ! update stability 
     207         Cx_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10) 
     208         Chn(:,:) = Cx_n10 
     209 
     210         !! Update of transfer coefficients: 
     211         ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u)) 
     212         Cd      = ztmp0 / ( ztmp1*ztmp1 ) 
     213         stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 
     214 
     215         ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
     216         ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd)) 
     217         ztmp1 = 1. + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10) 
     218         Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b) 
     219 
     220         Cx_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10 
     221         Cen(:,:) = Cx_n10 
     222         ztmp1 = 1. + Cx_n10*ztmp0 
     223         Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
     224         ENDIF 
    224225         ! 
    225226      END DO 
     
    252253            ! 
    253254            ! When wind speed > 33 m/s => Cyclone conditions => special treatment 
    254             zgt33 = 0.5 + SIGN( 0.5, (zw - 33.) )   ! If pw10 < 33. => 0, else => 1 
     255            zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) )   ! If pw10 < 33. => 0, else => 1 
    255256            ! 
    256257            cd_neutral_10m(ji,jj) = 1.e-3 * ( & 
     
    258259               &      +    zgt33   *      2.34 )                                     ! wind >= 33 m/s 
    259260            ! 
    260             cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6) 
     261            cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp) 
    261262            ! 
    262263         END DO 
     
    285286      DO jj = 1, jpj 
    286287         DO ji = 1, jpi 
    287             zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) ) 
    288             zx2 = MAX ( zx2 , 1. ) 
     288            zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 
     289            zx2 = MAX( zx2 , 1._wp ) 
    289290            zx  = SQRT( zx2 ) 
    290             zstab = 0.5 + SIGN( 0.5 , pzeta(ji,jj) ) 
    291             ! 
    292             psi_m(ji,jj) =        zstab  * (-5.*pzeta(ji,jj))       &          ! Stable 
    293                &          + (1. - zstab) * (2.*LOG((1. + zx)*0.5)   &          ! Unstable 
    294                &               + LOG((1. + zx2)*0.5) - 2.*ATAN(zx) + rpi*0.5)  !    " 
     291            zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 
     292            ! 
     293            psi_m(ji,jj) =        zstab  * (-5._wp*pzeta(ji,jj))       &          ! Stable 
     294               &          + (1._wp - zstab) * (2._wp*LOG((1._wp + zx)*0.5_wp)   &          ! Unstable 
     295               &               + LOG((1._wp + zx2)*0.5_wp) - 2._wp*ATAN(zx) + rpi*0.5_wp)  !    " 
    295296            ! 
    296297         END DO 
     
    319320      DO jj = 1, jpj 
    320321         DO ji = 1, jpi 
    321             zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) ) 
    322             zx2 = MAX ( zx2 , 1. ) 
    323             zstab = 0.5 + SIGN( 0.5 , pzeta(ji,jj) ) 
    324             ! 
    325             psi_h(ji,jj) =         zstab  * (-5.*pzeta(ji,jj))        &  ! Stable 
    326                &           + (1. - zstab) * (2.*LOG( (1. + zx2)*0.5 ))   ! Unstable 
     322            zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 
     323            zx2 = MAX( zx2 , 1._wp ) 
     324            zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 
     325            ! 
     326            psi_h(ji,jj) =         zstab  * (-5._wp*pzeta(ji,jj))        &  ! Stable 
     327               &           + (1._wp - zstab) * (2._wp*LOG( (1._wp + zx2)*0.5_wp ))   ! Unstable 
    327328            ! 
    328329         END DO 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90

    r11182 r11209  
    88   !!---------------------------------------------------------------------- 
    99 
     10   !!   virt_temp     : virtual (aka sensible) temperature (potential or absolute) 
    1011   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP 
     12   !!   visc_air      : kinematic viscosity (aka Nu_air) of air from temperature 
     13   !!   L_vap         : latent heat of vaporization of water as a function of temperature 
    1114   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air) 
     15   !!   gamma_moist   : adiabatic lapse-rate of moist air 
     16   !!   One_on_L      : 1. / ( Monin-Obukhov length ) 
     17   !!   Ri_bulk       : bulk Richardson number aka BRN 
    1218   !!   q_sat         : saturation humidity as a function of SLP and temperature 
    13    !!   gamma_moist   :  
    14    !!   L_vap         : latent heat of vaporization of water as a function of temperature 
    15    !!   visc_air      : kinematic viscosity (aka Nu_air) of air from temperature 
    16     
     19 
    1720   USE dom_oce        ! ocean space and time domain 
    1821   USE phycst         ! physical constants 
     
    2124   PRIVATE 
    2225 
     26   INTERFACE gamma_moist 
     27      MODULE PROCEDURE gamma_moist_vctr, gamma_moist_sclr 
     28   END INTERFACE gamma_moist 
     29 
     30   PUBLIC virt_temp 
    2331   PUBLIC rho_air 
     32   PUBLIC visc_air 
     33   PUBLIC L_vap 
    2434   PUBLIC cp_air 
     35   PUBLIC gamma_moist 
     36   PUBLIC One_on_L 
     37   PUBLIC Ri_bulk 
    2538   PUBLIC q_sat 
    26    PUBLIC gamma_moist 
    27    PUBLIC L_vap 
    28    PUBLIC visc_air 
    2939 
    3040   !!---------------------------------------------------------------------- 
     
    3545CONTAINS 
    3646 
     47   FUNCTION virt_temp( pta, pqa ) 
     48      !!------------------------------------------------------------------------ 
     49      !! 
     50      !! Compute the (absolute/potential) virtual temperature, knowing the 
     51      !! (absolute/potential) temperature and specific humidity 
     52      !! 
     53      !! If input temperature is absolute then output vitual temperature is absolute 
     54      !! If input temperature is potential then output vitual temperature is potential 
     55      !! 
     56      !! Author: L. Brodeau, June 2019 / AeroBulk 
     57      !!         (https://github.com/brodeau/aerobulk/) 
     58      !!------------------------------------------------------------------------ 
     59      REAL(wp), DIMENSION(jpi,jpj)             :: virt_temp         !: 1./(Monin Obukhov length) [m^-1] 
     60      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta,  &  !: absolute or potetntial air temperature [K] 
     61         &                                        pqa      !: specific humidity of air   [kg/kg] 
     62      !!------------------------------------------------------------------- 
     63      ! 
     64      virt_temp(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:)) 
     65      !! 
     66      !! This is exactly the same sing that: 
     67      !! virt_temp = pta * ( pwa + reps0) / (reps0*(1.+pwa)) 
     68      !! with wpa (mixing ration) defined as : pwa = pqa/(1.-pqa) 
     69      ! 
     70   END FUNCTION virt_temp 
     71 
    3772   FUNCTION rho_air( ptak, pqa, pslp ) 
    3873      !!------------------------------------------------------------------------------- 
     
    4176      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
    4277      !! 
    43       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     78      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    4479      !!------------------------------------------------------------------------------- 
    4580      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K] 
     
    5388   END FUNCTION rho_air 
    5489 
     90   FUNCTION visc_air(ptak) 
     91      !!---------------------------------------------------------------------------------- 
     92      !! Air kinetic viscosity (m^2/s) given from temperature in degrees... 
     93      !! 
     94      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     95      !!---------------------------------------------------------------------------------- 
     96      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air   ! 
     97      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature in (K) 
     98      ! 
     99      INTEGER  ::   ji, jj      ! dummy loop indices 
     100      REAL(wp) ::   ztc, ztc2   ! local scalar 
     101      !!---------------------------------------------------------------------------------- 
     102      ! 
     103      DO jj = 1, jpj 
     104         DO ji = 1, jpi 
     105            ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C 
     106            ztc2 = ztc*ztc 
     107            visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc) 
     108         END DO 
     109      END DO 
     110      ! 
     111   END FUNCTION visc_air 
     112 
     113   FUNCTION L_vap( psst ) 
     114      !!--------------------------------------------------------------------------------- 
     115      !!                           ***  FUNCTION L_vap  *** 
     116      !! 
     117      !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
     118      !! 
     119      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     120      !!---------------------------------------------------------------------------------- 
     121      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg] 
     122      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
     123      !!---------------------------------------------------------------------------------- 
     124      ! 
     125      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6 
     126      ! 
     127   END FUNCTION L_vap 
    55128 
    56129   FUNCTION cp_air( pqa ) 
     
    69142      ! 
    70143   END FUNCTION cp_air 
     144 
     145   FUNCTION gamma_moist_vctr( ptak, pqa ) 
     146      !!---------------------------------------------------------------------------------- 
     147      !!                           ***  FUNCTION gamma_moist_vctr  *** 
     148      !! 
     149      !! ** Purpose : Compute the moist adiabatic lapse-rate. 
     150      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 
     151      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 
     152      !! 
     153      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     154      !!---------------------------------------------------------------------------------- 
     155      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K] 
     156      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg] 
     157      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist_vctr   ! moist adiabatic lapse-rate 
     158      ! 
     159      INTEGER  ::   ji, jj         ! dummy loop indices 
     160      REAL(wp) :: zta, zqa, zwa, ziRT        ! local scalar 
     161      !!---------------------------------------------------------------------------------- 
     162      ! 
     163      DO jj = 1, jpj 
     164         DO ji = 1, jpi 
     165            zta = MAX( ptak(ji,jj),  180._wp) ! prevents screw-up over masked regions where field == 0. 
     166            zqa = MAX( pqa(ji,jj),  1.E-6_wp) !    "                   "                     " 
     167            ! 
     168            zwa = zqa / (1. - zqa)   ! w is mixing ratio w = q/(1-q) | q = w/(1+w) 
     169            ziRT = 1._wp/(R_dry*zta)    ! 1/RT 
     170            gamma_moist_vctr(ji,jj) = grav * ( 1._wp + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/zta ) 
     171         END DO 
     172      END DO 
     173      ! 
     174   END FUNCTION gamma_moist_vctr 
     175 
     176   FUNCTION gamma_moist_sclr( ptak, pqa ) 
     177      !!---------------------------------------------------------------------------------- 
     178      !! ** Purpose : Compute the moist adiabatic lapse-rate. 
     179      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 
     180      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 
     181      !! 
     182      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     183      !!---------------------------------------------------------------------------------- 
     184      REAL(wp)             :: gamma_moist_sclr 
     185      REAL(wp), INTENT(in) :: ptak, pqa ! air temperature (K) and specific humidity (kg/kg) 
     186      ! 
     187      REAL(wp) :: zta, zqa, zwa, ziRT        ! local scalar 
     188      !!---------------------------------------------------------------------------------- 
     189      zta = MAX( ptak,  180._wp) ! prevents screw-up over masked regions where field == 0. 
     190      zqa = MAX( pqa,  1.E-6_wp) !    "                   "                     " 
     191      !! 
     192      zwa = zqa / (1._wp - zqa)   ! w is mixing ratio w = q/(1-q) | q = w/(1+w) 
     193      ziRT = 1._wp / (R_dry*zta)    ! 1/RT 
     194      gamma_moist_sclr = grav * ( 1._wp + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/zta ) 
     195      !! 
     196   END FUNCTION gamma_moist_sclr 
     197 
     198   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs ) 
     199      !!------------------------------------------------------------------------ 
     200      !! 
     201      !! Evaluates the 1./(Monin Obukhov length) from air temperature and 
     202      !!  specific humidity, and frictional scales u*, t* and q* 
     203      !! 
     204      !! Author: L. Brodeau, June 2016 / AeroBulk 
     205      !!         (https://github.com/brodeau/aerobulk/) 
     206      !!------------------------------------------------------------------------ 
     207      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1] 
     208      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha,  &  !: average potetntial air temperature [K] 
     209         &                                        pqa,   &  !: average specific humidity of air   [kg/kg] 
     210         &                                      pus, pts, pqs   !: frictional velocity, temperature and humidity 
     211      ! 
     212      INTEGER  ::   ji, jj         ! dummy loop indices 
     213      REAL(wp) ::     zqa          ! local scalar 
     214      !!------------------------------------------------------------------- 
     215      ! 
     216      DO jj = 1, jpj 
     217         DO ji = 1, jpi 
     218            ! 
     219            zqa = (1._wp + rctv0*pqa(ji,jj)) 
     220            ! 
     221            One_on_L(ji,jj) = grav*vkarmn*(pts(ji,jj) + rctv0*ptha(ji,jj)*pqs(ji,jj)) & 
     222               &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp ) 
     223            ! 
     224         END DO 
     225      END DO 
     226      ! 
     227      One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...) 
     228      ! 
     229   END FUNCTION One_on_L 
     230 
     231   FUNCTION Ri_bulk( pz, psst, ptha, pssq, pqa, pub ) 
     232      !!---------------------------------------------------------------------------------- 
     233      !! Bulk Richardson number according to "wide-spread equation"... 
     234      !! 
     235      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     236      !!---------------------------------------------------------------------------------- 
     237      REAL(wp), DIMENSION(jpi,jpj)             :: Ri_bulk 
     238      REAL(wp)                    , INTENT(in) :: pz    ! height above the sea (aka "delta z")  [m] 
     239      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst  ! SST                                   [K] 
     240      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha  ! pot. air temp. at height "pz"         [K] 
     241      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssq  ! 0.98*q_sat(SST)                   [kg/kg] 
     242      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa   ! air spec. hum. at height "pz"     [kg/kg] 
     243      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub   ! bulk wind speed                     [m/s] 
     244      ! 
     245      INTEGER  ::   ji, jj                                ! dummy loop indices 
     246      REAL(wp) ::   zqa, zta, zgamma, zdth_v, ztv, zsstv  ! local scalars 
     247      !!------------------------------------------------------------------- 
     248      ! 
     249      DO jj = 1, jpj 
     250         DO ji = 1, jpi 
     251            ! 
     252            zqa = 0.5_wp*(pqa(ji,jj)+pssq(ji,jj))                                        ! ~ mean q within the layer... 
     253            zta = 0.5_wp*( psst(ji,jj) + ptha(ji,jj) - gamma_moist(ptha(ji,jj),zqa)*pz ) ! ~ mean absolute temperature of air within the layer 
     254            zta = 0.5_wp*( psst(ji,jj) + ptha(ji,jj) - gamma_moist(zta,        zqa)*pz ) ! ~ mean absolute temperature of air within the layer 
     255            zgamma =  gamma_moist(zta, zqa)                                              ! Adiabatic lapse-rate for moist air within the layer 
     256            ! 
     257            zsstv = psst(ji,jj)*(1._wp + rctv0*pssq(ji,jj)) ! absolute==potential virtual SST (absolute==potential because z=0!) 
     258            ! 
     259            zdth_v = ptha(ji,jj)*(1._wp + rctv0*pqa(ji,jj)) - zsstv ! air-sea delta of "virtual potential temperature" 
     260            ! 
     261            ztv = 0.5_wp*( zsstv + (ptha(ji,jj) - zgamma*pz)*(1._wp + rctv0*pqa(ji,jj)) )  ! ~ mean absolute virtual temp. within the layer 
     262            ! 
     263            Ri_bulk(ji,jj) = grav*zdth_v*pz / ( ztv*pub(ji,jj)*pub(ji,jj) )                            ! the usual definition of Ri_bulk 
     264            ! 
     265         END DO 
     266      END DO 
     267   END FUNCTION Ri_bulk 
    71268 
    72269 
     
    108305 
    109306 
    110    FUNCTION gamma_moist( ptak, pqa ) 
    111       !!---------------------------------------------------------------------------------- 
    112       !!                           ***  FUNCTION gamma_moist  *** 
    113       !! 
    114       !! ** Purpose : Compute the moist adiabatic lapse-rate. 
    115       !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 
    116       !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 
    117       !! 
    118       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    119       !!---------------------------------------------------------------------------------- 
    120       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K] 
    121       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg] 
    122       REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate 
    123       ! 
    124       INTEGER  ::   ji, jj         ! dummy loop indices 
    125       REAL(wp) :: zrv, ziRT        ! local scalar 
    126       !!---------------------------------------------------------------------------------- 
    127       ! 
    128       DO jj = 1, jpj 
    129          DO ji = 1, jpi 
    130             zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 
    131             ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT 
    132             gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( rCp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) ) 
    133          END DO 
    134       END DO 
    135       ! 
    136    END FUNCTION gamma_moist 
    137  
    138  
    139    FUNCTION L_vap( psst ) 
    140       !!--------------------------------------------------------------------------------- 
    141       !!                           ***  FUNCTION L_vap  *** 
    142       !! 
    143       !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
    144       !! 
    145       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    146       !!---------------------------------------------------------------------------------- 
    147       REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg] 
    148       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
    149       !!---------------------------------------------------------------------------------- 
    150       ! 
    151       L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6 
    152       ! 
    153    END FUNCTION L_vap 
    154  
    155  
    156  
    157    FUNCTION visc_air(ptak) 
    158       !!---------------------------------------------------------------------------------- 
    159       !! Air kinetic viscosity (m^2/s) given from temperature in degrees... 
    160       !! 
    161       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    162       !!---------------------------------------------------------------------------------- 
    163       REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air   ! 
    164       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature in (K) 
    165       ! 
    166       INTEGER  ::   ji, jj      ! dummy loop indices 
    167       REAL(wp) ::   ztc, ztc2   ! local scalar 
    168       !!---------------------------------------------------------------------------------- 
    169       ! 
    170       DO jj = 1, jpj 
    171          DO ji = 1, jpi 
    172             ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C 
    173             ztc2 = ztc*ztc 
    174             visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc) 
    175          END DO 
    176       END DO 
    177       ! 
    178    END FUNCTION visc_air 
    179     
    180  
    181307   !!====================================================================== 
    182308END MODULE sbcblk_phy 
Note: See TracChangeset for help on using the changeset viewer.