- Timestamp:
- 2019-07-05T16:44:30+02:00 (5 years ago)
- 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 26 26 USE dom_oce ! ocean space and time domain 27 27 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 29 32 USE sbcwave, ONLY : cdn_wave ! wave module 30 33 #if defined key_si3 || defined key_cice 31 34 USE sbc_ice ! Surface boundary condition: ice fields 32 35 #endif 33 !34 USE iom ! I/O manager library35 USE lib_mpp ! distribued memory computing library36 USE in_out_manager ! I/O manager37 USE prtctl ! Print control38 36 USE lib_fortran ! to use key_nosignedzero 39 37 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 41 40 42 41 IMPLICIT NONE … … 45 44 PUBLIC :: TURB_NCAR ! called by sbcblk.F90 46 45 46 INTEGER , PARAMETER :: nb_itt = 4 ! number of itterations 47 47 48 !!---------------------------------------------------------------------- 48 49 CONTAINS … … 59 60 !! Returns the effective bulk wind speed at 10m to be used in the bulk formulas 60 61 !! 61 !! ** Method : Monin Obukhov Similarity Theory62 !! + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10)63 !!64 !! ** References : Large & Yeager, 2004 / Large & Yeager, 200865 !!66 !! ** Last update: Laurent Brodeau, June 2014:67 !! - handles both cases zt=zu and zt/=zu68 !! - optimized: less 2D arrays allocated and less operations69 !! - better first guess of stability by checking air-sea difference of virtual temperature70 !! rather than temperature difference only...71 !! - added function "cd_neutral_10m" that uses the improved parametrization of72 !! Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions!73 !! - using code-wide physical constants defined into "phycst.mod" rather than redifining them74 !! => 'vkarmn' and 'grav'75 !!76 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)77 62 !! 78 63 !! INPUT : … … 95 80 !! * q_zu : specific humidity of air // [kg/kg] 96 81 !! * U_blk : bulk wind speed at 10m [m/s] 82 !! 83 !! 84 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 97 85 !!---------------------------------------------------------------------------------- 98 86 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] … … 111 99 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients 112 100 ! 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 116 103 ! 117 104 REAL(wp), DIMENSION(jpi,jpj) :: Cx_n10 ! 10m neutral latent/sensible coefficient … … 124 111 ! 125 112 l_zt_equal_zu = .FALSE. 126 IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision113 IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 127 114 128 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 … … 137 124 ztmp0 (:,:) = cdn_wave(:,:) 138 125 ELSE 139 126 ztmp0 = cd_neutral_10m( U_blk ) 140 127 ENDIF 141 128 … … 144 131 !! Initializing transf. coeff. with their first guess neutral equivalents : 145 132 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)) 148 135 stab = sqrt_Cd_n10 ! Temporaty array !!! stab == SQRT(Cd) 149 136 150 IF( ln_cdgw ) Cen = Ce ; Chn = Ch 137 IF( ln_cdgw ) THEN 138 Cen = Ce 139 Chn = Ch 140 ENDIF 151 141 152 142 !! Initializing values at z_u with z_t values: … … 169 159 !! Stability parameters : 170 160 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 ) 172 162 zpsi_h_u = psi_h( zeta_u ) 173 163 … … 176 166 !! Array 'stab' is free for the moment so using it to store 'zeta_t' 177 167 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 !!! 179 169 stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab) ! stab just used as temp array again! 180 170 t_zu = t_zt - ztmp1/vkarmn*stab ! ztmp1 is still theta* L&Y 2004 eq.(9b) … … 187 177 stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 ) ! (stab == SQRT(Cd)) 188 178 Cd = stab * stab 189 ztmp0 = (LOG(zu/10. ) - zpsi_h_u) / vkarmn / sqrt_Cd_n10179 ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 190 180 ztmp2 = stab / sqrt_Cd_n10 ! (stab == SQRT(Cd)) 191 ztmp1 = 1. + Chn * ztmp0181 ztmp1 = 1._wp + Chn * ztmp0 192 182 Ch = Chn * ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 193 ztmp1 = 1. + Cen * ztmp0183 ztmp1 = 1._wp + Cen * ztmp0 194 184 Ce = Cen * ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 195 185 … … 205 195 206 196 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) 208 198 Chn(:,:) = Cx_n10 209 199 210 200 !! 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)) 212 202 Cd = ztmp0 / ( ztmp1*ztmp1 ) 213 203 stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 214 204 215 ztmp0 = (LOG(zu/10. ) - zpsi_h_u) / vkarmn / sqrt_Cd_n10205 ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 216 206 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) 218 208 Ch = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 219 209 220 Cx_n10 = 1.e-3 * (34.6* sqrt_Cd_n10) ! L&Y 2004 eq. (6b) ! Cx_n10 == Ce_n10210 Cx_n10 = 1.e-3_wp * (34.6_wp * sqrt_Cd_n10) ! L&Y 2004 eq. (6b) ! Cx_n10 == Ce_n10 221 211 Cen(:,:) = Cx_n10 222 ztmp1 = 1. + Cx_n10*ztmp0212 ztmp1 = 1._wp + Cx_n10*ztmp0 223 213 Ce = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 224 214 ENDIF … … 255 245 zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) ) ! If pw10 < 33. => 0, else => 1 256 246 ! 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/s259 & + zgt33 * 2.34 )! wind >= 33 m/s247 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 260 250 ! 261 251 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.