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 11221 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90 – NEMO

Ignore:
Timestamp:
2019-07-05T16:44:30+02:00 (5 years ago)
Author:
laurent
Message:

LB: "sbcblk_algo_ncar.F90" now uses functions "virt_temp()" and "One_on_L()" defined in "sbcblk_phy.F90".

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90

    r11209 r11221  
    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 
    40    USE sbcblk_phy      !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB 
     38   USE sbc_oce         ! Surface boundary condition: ocean fields 
     39   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB 
    4140 
    4241   IMPLICIT NONE 
    … …  
    4544   PUBLIC ::   TURB_NCAR   ! called by sbcblk.F90 
    4645    
     46   INTEGER , PARAMETER ::   nb_itt = 4        ! number of itterations 
     47 
    4748   !!---------------------------------------------------------------------- 
    4849CONTAINS 
    … …  
    5960      !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas 
    6061      !! 
    61       !! ** Method : Monin Obukhov Similarity Theory 
    62       !!             + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10) 
    63       !! 
    64       !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008 
    65       !! 
    66       !! ** Last update: Laurent Brodeau, June 2014: 
    67       !!    - handles both cases zt=zu and zt/=zu 
    68       !!    - optimized: less 2D arrays allocated and less operations 
    69       !!    - better first guess of stability by checking air-sea difference of virtual temperature 
    70       !!       rather than temperature difference only... 
    71       !!    - added function "cd_neutral_10m" that uses the improved parametrization of 
    72       !!      Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions! 
    73       !!    - using code-wide physical constants defined into "phycst.mod" rather than redifining them 
    74       !!      => 'vkarmn' and 'grav' 
    75       !! 
    76       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    7762      !! 
    7863      !! INPUT : 
    … …  
    9580      !!    *  q_zu   : specific humidity of air        //                    [kg/kg] 
    9681      !!    *  U_blk  : bulk wind speed at 10m                                [m/s] 
     82      !! 
     83      !! 
     84      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    9785      !!---------------------------------------------------------------------------------- 
    9886      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m] 
    … …  
    11199      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    112100      ! 
    113       INTEGER ::   j_itt 
    114       LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
    115       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 
    116103      ! 
    117104      REAL(wp), DIMENSION(jpi,jpj) ::   Cx_n10        ! 10m neutral latent/sensible coefficient 
    … …  
    124111      ! 
    125112      l_zt_equal_zu = .FALSE. 
    126       IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
     113      IF( ABS(zu - zt) < 0.01 )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
    127114 
    128115      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 
    … …  
    137124         ztmp0   (:,:) = cdn_wave(:,:) 
    138125      ELSE 
    139          ztmp0 = cd_neutral_10m( U_blk ) 
     126      ztmp0 = cd_neutral_10m( U_blk ) 
    140127      ENDIF 
    141128 
    … …  
    144131      !! Initializing transf. coeff. with their first guess neutral equivalents : 
    145132      Cd = ztmp0 
    146       Ce = 1.e-3*( 34.6 * sqrt_Cd_n10 ) 
    147       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)) 
    148135      stab = sqrt_Cd_n10   ! Temporaty array !!! stab == SQRT(Cd) 
    149136  
    150       IF( ln_cdgw )   Cen = Ce  ; Chn = Ch 
     137      IF( ln_cdgw ) THEN 
     138   Cen = Ce 
     139   Chn = Ch 
     140      ENDIF 
    151141 
    152142      !! Initializing values at z_u with z_t values: 
    … …  
    169159         !! Stability parameters : 
    170160         zeta_u   = zu*ztmp0 
    171          zeta_u = sign( min(abs(zeta_u),10.0_wp), zeta_u ) 
     161         zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u ) 
    172162         zpsi_h_u = psi_h( zeta_u ) 
    173163 
    … …  
    176166            !! Array 'stab' is free for the moment so using it to store 'zeta_t' 
    177167            stab = zt*ztmp0 
    178             stab = SIGN( MIN(ABS(stab),10.0_wp), stab )  ! Temporaty array stab == zeta_t !!! 
     168            stab = SIGN( MIN(ABS(stab),10._wp), stab )  ! Temporaty array stab == zeta_t !!! 
    179169            stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab)                   ! stab just used as temp array again! 
    180170            t_zu = t_zt - ztmp1/vkarmn*stab    ! ztmp1 is still theta*  L&Y 2004 eq.(9b) 
    … …  
    187177            stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 )  ! (stab == SQRT(Cd)) 
    188178            Cd   = stab * stab 
    189             ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
     179            ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
    190180            ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd)) 
    191             ztmp1 = 1. + Chn * ztmp0      
     181            ztmp1 = 1._wp + Chn * ztmp0      
    192182            Ch    = Chn * ztmp2 / ztmp1  ! L&Y 2004 eq. (10b) 
    193             ztmp1 = 1. + Cen * ztmp0 
     183            ztmp1 = 1._wp + Cen * ztmp0 
    194184            Ce    = Cen * ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
    195185 
    … …  
    205195 
    206196         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) 
     197         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) 
    208198         Chn(:,:) = Cx_n10 
    209199 
    210200         !! 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)) 
     201         ztmp1 = 1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u)) 
    212202         Cd      = ztmp0 / ( ztmp1*ztmp1 ) 
    213203         stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 
    214204 
    215          ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
     205         ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
    216206         ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd)) 
    217          ztmp1 = 1. + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10) 
     207         ztmp1 = 1._wp + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10) 
    218208         Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b) 
    219209 
    220          Cx_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10 
     210         Cx_n10  = 1.e-3_wp * (34.6_wp * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10 
    221211         Cen(:,:) = Cx_n10 
    222          ztmp1 = 1. + Cx_n10*ztmp0 
     212         ztmp1 = 1._wp + Cx_n10*ztmp0 
    223213         Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
    224214         ENDIF 
    … …  
    255245            zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) )   ! If pw10 < 33. => 0, else => 1 
    256246            ! 
    257             cd_neutral_10m(ji,jj) = 1.e-3 * ( & 
    258                &       (1. - zgt33)*( 2.7/zw + 0.142 + zw/13.09 - 3.14807E-10*zw6) & ! wind <  33 m/s 
    259                &      +    zgt33   *      2.34 )                                     ! wind >= 33 m/s 
     247            cd_neutral_10m(ji,jj) = 1.e-3_wp * ( & 
     248               &       (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind <  33 m/s 
     249               &      +    zgt33   *      2.34_wp )                                                 ! wind >= 33 m/s 
    260250            ! 
    261251            cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp) 
Note: See TracChangeset for help on using the changeset viewer.