Changeset 11111 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90
- Timestamp:
- 2019-06-14T15:55:28+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_ecmwf.F90
r10069 r11111 1 1 MODULE sbcblk_algo_ecmwf 2 2 !!====================================================================== 3 !! *** MODULE sbcblk_algo_ecmwf *** 4 !! Computes turbulent components of surface fluxes 5 !! according to the method in IFS of the ECMWF model 6 !! 3 !! *** MODULE sbcblk_algo_ecmwf *** 4 !! Computes: 7 5 !! * bulk transfer coefficients C_D, C_E and C_H 8 6 !! * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed … … 10 8 !! => all these are used in bulk formulas in sbcblk.F90 11 9 !! 12 !! Using the bulk formulation/param. of IFS of ECMWF (cycle 31r2)10 !! Using the bulk formulation/param. of IFS of ECMWF (cycle 40r1) 13 11 !! based on IFS doc (avaible online on the ECMWF's website) 14 12 !! 13 !! Routine turb_ecmwf maintained and developed in AeroBulk 14 !! (https://github.com/brodeau/aerobulk) 15 15 !! 16 !! Routine turb_ecmwf maintained and developed in AeroBulk 17 !! (http://aerobulk.sourceforge.net/) 18 !! 19 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 16 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk) 20 17 !!---------------------------------------------------------------------- 21 18 !! History : 4.0 ! 2016-02 (L.Brodeau) Original code … … 41 38 42 39 USE sbc_oce ! Surface boundary condition: ocean fields 40 USE sbcblk_phymbl !LB: all thermodynamics functions, rho_air, q_sat, etc... 41 USE sbcblk_skin ! cool-skin/warm layer scheme (CSWL_ECMWF) 43 42 44 43 IMPLICIT NONE … … 52 51 REAL(wp), PARAMETER :: zi0 = 1000. ! scale height of the atmospheric boundary layer...1 53 52 REAL(wp), PARAMETER :: Beta0 = 1. ! gustiness parameter ( = 1.25 in COAREv3) 54 REAL(wp), PARAMETER :: rctv0 = 0.608 ! constant to obtain virtual temperature...55 REAL(wp), PARAMETER :: Cp_dry = 1005.0 ! Specic heat of dry air, constant pressure [J/K/kg]56 REAL(wp), PARAMETER :: Cp_vap = 1860.0 ! Specic heat of water vapor, constant pressure [J/K/kg]57 53 REAL(wp), PARAMETER :: alpha_M = 0.11 ! For roughness length (smooth surface term) 58 54 REAL(wp), PARAMETER :: alpha_H = 0.40 ! (Chapter 3, p.34, IFS doc Cy31r1) 59 55 REAL(wp), PARAMETER :: alpha_Q = 0.62 ! 56 57 INTEGER , PARAMETER :: nb_itt = 5 ! number of itterations 58 59 !! Cool-Skin / Warm-Layer related parameters: 60 REAL(wp), PARAMETER :: & 61 & rdt0 = 3600.*1.5, & !: time step 62 & rd0 = 3. , & !: Depth scale [m], "d" in Eq.11 (Zeng & Beljaars 2005) 63 & rNu0 = 0.5 !: Nu (exponent of temperature profile) Eq.11 64 ! !: (Zeng & Beljaars 2005) !: set to 0.5 instead of 65 ! !: 0.3 to respect a warming of +3 K in calm 66 ! !: condition for the insolation peak of +1000W/m^2 67 INTEGER, PARAMETER :: & 68 & nb_itt_wl = 10 !: number of sub-itterations for solving the differential equation in warm-layer part 69 ! !: => use "nb_itt_wl = 1" for No itteration! => way cheaper !!! 70 ! !: => assumes balance between the last 2 terms of Eq.11 (lhs of eq.11 = 0) 71 ! !: => in that case no need for sub-itterations ! 72 ! !: => ACCEPTABLE IN MOST CONDITIONS ! (UNLESS: sunny + very calm/low-wind conditions) 73 ! !: => Otherwize use "nb_itt_wl = 10" 74 75 76 77 60 78 !!---------------------------------------------------------------------- 61 79 CONTAINS 62 80 63 SUBROUTINE TURB_ECMWF( zt, zu, sst, t_zt, ssq , q_zt , U_zu, & 64 & Cd, Ch, Ce , t_zu, q_zu, U_blk, & 65 & Cdn, Chn, Cen ) 81 SUBROUTINE TURB_ECMWF( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, & 82 & Cd, Ch, Ce, t_zu, q_zu, U_blk, & 83 & Cdn, Chn, Cen, & 84 & Qsw, rad_lw, slp ) 66 85 !!---------------------------------------------------------------------------------- 67 86 !! *** ROUTINE turb_ecmwf *** 68 87 !! 69 !! 2015: L. Brodeau (brodeau@gmail.com)70 !!71 88 !! ** Purpose : Computes turbulent transfert coefficients of surface 72 !! fluxes according to IFS doc. (cycle 31)89 !! fluxes according to IFS doc. (cycle 40) 73 90 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 74 !! 75 !! ** Method : Monin Obukhov Similarity Theory 91 !! Returns the effective bulk wind speed at 10m to be used in the bulk formulas 92 !! 93 !! Applies the cool-skin warm-layer correction of the SST to T_s 94 !! if the net shortwave flux at the surface (Qsw), the downwelling longwave 95 !! radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp) 96 !! are provided as (optional) arguments! 76 97 !! 77 98 !! INPUT : … … 80 101 !! * zu : height for wind speed (generally 10m) [m] 81 102 !! * U_zu : scalar wind speed at 10m [m/s] 82 !! * sst : SST [K]83 103 !! * t_zt : potential air temperature at zt [K] 84 !! * ssq : specific humidity at saturation at SST [kg/kg]85 104 !! * q_zt : specific humidity of air at zt [kg/kg] 86 105 !! 106 !! INPUT/OUTPUT: 107 !! ------------- 108 !! * T_s : SST or skin temperature [K] 109 !! * q_s : SSQ aka saturation specific humidity at temp. T_s [kg/kg] 110 !! -> doesn't need to be given a value if skin temp computed (in case l_use_skin=True) 111 !! -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False) 112 !! 113 !! OPTIONAL INPUT (will trigger l_use_skin=TRUE if present!): 114 !! --------------- 115 !! * Qsw : net solar flux (after albedo) at the surface (>0) [W/m^2] 116 !! * rad_lw : downwelling longwave radiation at the surface (>0) [W/m^2] 117 !! * slp : sea-level pressure [Pa] 87 118 !! 88 119 !! OUTPUT : … … 93 124 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 94 125 !! * q_zu : specific humidity of air // [kg/kg] 95 !! * U_blk : bulk wind at 10m[m/s]96 !! 97 !! 98 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)126 !! * U_blk : bulk wind speed at 10m [m/s] 127 !! 128 !! 129 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 99 130 !!---------------------------------------------------------------------------------- 100 131 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] 101 132 REAL(wp), INTENT(in ) :: zu ! height for U_zu [m] 102 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: sst! sea surface temperature [Kelvin]133 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: T_s ! sea surface temperature [Kelvin] 103 134 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin] 104 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 135 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: q_s ! sea surface specific humidity [kg/kg] 136 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg] 106 137 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s] 107 138 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau) … … 113 144 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients 114 145 ! 146 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: Qsw ! [W/m^2] 147 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: rad_lw ! [W/m^2] 148 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: slp ! [Pa] 149 ! 115 150 INTEGER :: j_itt 116 LOGICAL :: 117 INTEGER , PARAMETER :: nb_itt = 4 ! number of itterations118 !119 REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star,&151 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 152 ! 153 REAL(wp), DIMENSION(jpi,jpj) :: & 154 & u_star, t_star, q_star, & 120 155 & dt_zu, dq_zu, & 121 156 & znu_a, & !: Nu_air, Viscosity of air 122 157 & Linv, & !: 1/L (inverse of Monin Obukhov length... 123 158 & z0, z0t, z0q 159 ! 160 ! Cool skin: 161 LOGICAL :: l_use_skin = .FALSE. 162 REAL(wp), DIMENSION(jpi,jpj) :: zsst ! to back up the initial bulk SST 163 164 ! 124 165 REAL(wp), DIMENSION(jpi,jpj) :: func_m, func_h 125 166 REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 126 167 !!---------------------------------------------------------------------------------- 127 168 ! 169 ! Cool skin ? 170 IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN 171 l_use_skin = .TRUE. 172 END IF 173 IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin 174 128 175 ! Identical first gess as in COARE, with IFS parameter values though 129 176 ! … … 131 178 IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 132 179 180 !! Initialization for cool skin: 181 IF( l_use_skin ) THEN 182 zsst = T_s ! save the bulk SST 183 T_s = T_s - 0.25 ! First guess of correction 184 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 185 END IF 133 186 134 187 !! First guess of temperature and humidity at height zu: 135 t_zu = MAX( t_zt , 0.0 ) ! who knows what's given on masked-continental regions...136 q_zu = MAX( q_zt , 1.e-6 ) ! "188 t_zu = MAX( t_zt , 0.0_wp ) ! who knows what's given on masked-continental regions... 189 q_zu = MAX( q_zt , 1.e-6_wp) ! " 137 190 138 191 !! Pot. temp. difference (and we don't want it to be 0!) 139 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.e-6), dt_zu )140 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.e-9), dq_zu )192 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 193 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 141 194 142 195 znu_a = visc_air(t_zt) ! Air viscosity (m^2/s) at zt given from temperature in (K) 143 196 144 ztmp2 = 0.5 * 0.5! initial guess for wind gustiness contribution197 ztmp2 = 0.5_wp*0.5_wp ! initial guess for wind gustiness contribution 145 198 U_blk = SQRT(U_zu*U_zu + ztmp2) 146 199 147 ! z0 = 0.0001 148 ztmp2 = 10000. ! optimization: ztmp2 == 1/z0 200 ztmp2 = 10000._wp ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001) 149 201 ztmp0 = LOG(zu*ztmp2) 150 202 ztmp1 = LOG(10.*ztmp2) 151 u_star = 0.035*U_blk*ztmp1/ztmp0 ! (u* = 0.035*Un10) 152 153 z0 = charn0*u_star*u_star/grav + 0.11*znu_a/u_star 154 z0t = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ! WARNING: 1/z0t ! 155 156 Cd = (vkarmn/ztmp0)**2 ! first guess of Cd 157 158 ztmp0 = vkarmn*vkarmn/LOG(zt*z0t)/Cd 203 u_star = 0.035_wp*U_blk*ztmp1/ztmp0 ! (u* = 0.035*Un10) 204 205 z0 = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star 206 z0 = MIN(ABS(z0), 0.001_wp) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 207 z0t = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 208 z0t = MIN(ABS(z0t), 0.001_wp) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 209 210 ztmp2 = vkarmn/ztmp0 211 Cd = ztmp2*ztmp2 ! first guess of Cd 212 213 ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 159 214 160 215 ztmp2 = Ri_bulk( zu, t_zu, dt_zu, q_zu, dq_zu, U_blk ) ! Ribu = Bulk Richardson number 161 216 162 217 !! First estimate of zeta_u, depending on the stability, ie sign of Ribu (ztmp2): 163 ztmp1 = 0.5 + SIGN( 0.5 , ztmp2 )218 ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 164 219 func_m = ztmp0*ztmp2 ! temporary array !! 165 !! Ribu < 0 Ribu > 0 Beta = 1.25166 func_h = (1.-ztmp1)*(func_m/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & ! temporary array !!! func_h == zeta_u167 & + ztmp1*(func_m*(1. + 27./9.*ztmp2/ztmp0))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 222 !#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" ! 168 223 169 224 !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 170 ztmp0 = vkarmn/(LOG(zu*z0t) - psi_h_ecmwf(func_h))225 ztmp0 = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h)) 171 226 172 227 u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_ecmwf(func_h)) … … 176 231 ! What's need to be done if zt /= zu: 177 232 IF( .NOT. l_zt_equal_zu ) THEN 178 !179 233 !! First update of values at zu (or zt for wind) 180 234 ztmp0 = psi_h_ecmwf(func_h) - psi_h_ecmwf(zt*func_h/zu) ! zt*func_h/zu == zeta_t 181 ztmp1 = log(zt/zu) + ztmp0235 ztmp1 = LOG(zt/zu) + ztmp0 182 236 t_zu = t_zt - t_star/vkarmn*ztmp1 183 237 q_zu = q_zt - q_star/vkarmn*ztmp1 184 q_zu = (0.5 + sign(0.5,q_zu))*q_zu !Makes it impossible to have negative humidity : 185 186 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 187 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 238 q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity : 188 239 ! 189 ENDIF 240 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 241 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 242 END IF 190 243 191 244 … … 195 248 !! First guess of inverse of Monin-Obukov length (1/L) : 196 249 ztmp0 = (1. + rctv0*q_zu) ! the factor to apply to temp. to get virt. temp... 197 Linv = grav*vkarmn*(t_star*ztmp0 + rctv0*t_zu*q_star) / ( u_star*u_star * t_zu*ztmp0 ) 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 198 252 199 253 !! Functions such as u* = U_blk*vkarmn/func_m 200 ztmp1 = zu + z0 201 ztmp0 = ztmp1*Linv 202 func_m = LOG(ztmp1) -LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0*Linv) 203 func_h = LOG(ztmp1*z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(1./z0t*Linv) 204 254 ztmp0 = zu*Linv 255 func_m = LOG(zu) - LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv) 256 func_h = LOG(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 205 257 206 258 !! ITERATION BLOCK 207 !! ***************208 209 259 DO j_itt = 1, nb_itt 210 260 … … 213 263 214 264 !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) : 215 Linv = ztmp0*func_m*func_m/func_h / zu ! From Eq. 3.23, Chap.3, p.33, IFS doc - Cy31r1 265 Linv = ztmp0*func_m*func_m/func_h / zu ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1 266 !! Note: it is slightly different that the L we would get with the usual 267 !! expression, as in coare algorithm or in 'mod_phymbl.f90' (One_on_L_MO()) 268 Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 216 269 217 270 !! Update func_m with new Linv: 218 ztmp1 = zu + z0 219 func_m = LOG(ztmp1) -LOG(z0) - psi_m_ecmwf(ztmp1*Linv) + psi_m_ecmwf(z0*Linv) 271 func_m = LOG(zu) -LOG(z0) - psi_m_ecmwf(zu*Linv) + psi_m_ecmwf(z0*Linv) ! LB: should be "zu+z0" rather than "zu" alone, but z0 is tiny wrt zu! 220 272 221 273 !! Need to update roughness lengthes: … … 223 275 ztmp2 = u_star*u_star 224 276 ztmp1 = znu_a/u_star 225 z0 = alpha_M*ztmp1 + charn0*ztmp2/grav226 z0t = alpha_H*ztmp1! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1227 z0q = alpha_Q*ztmp1228 277 z0 = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp) 278 z0t = MIN( ABS( alpha_H*ztmp1 ) , 0.001_wp) ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 279 z0q = MIN( ABS( alpha_Q*ztmp1 ) , 0.001_wp) 280 229 281 !! Update wind at 10m taking into acount convection-related wind gustiness: 282 !! => Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8 230 283 ! Only true when unstable (L<0) => when ztmp0 < 0 => - !!! 231 ztmp2 = ztmp2 * ( MAX(-zi0*Linv/vkarmn,0.))**(2./3.) ! => w*^2 (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1)284 ztmp2 = ztmp2 * ( MAX(-zi0*Linv/vkarmn , 0._wp))**(2._wp/3._wp) ! => w*^2 (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1) 232 285 !! => equivalent using Beta=1 (gustiness parameter, 1.25 for COARE, also zi0=600 in COARE..) 233 U_blk = MAX( sqrt(U_zu*U_zu + ztmp2), 0.2) ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1286 U_blk = MAX( SQRT(U_zu*U_zu + ztmp2) , 0.2_wp ) ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1 234 287 ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 235 288 … … 238 291 !! as well the air-sea differences: 239 292 IF( .NOT. l_zt_equal_zu ) THEN 240 241 293 !! Arrays func_m and func_h are free for a while so using them as temporary arrays... 242 func_h = psi_h_ecmwf( (zu+z0)*Linv) ! temporary array !!!243 func_m = psi_h_ecmwf( (zt+z0)*Linv) ! temporary array !!!294 func_h = psi_h_ecmwf(zu*Linv) ! temporary array !!! 295 func_m = psi_h_ecmwf(zt*Linv) ! temporary array !!! 244 296 245 297 ztmp2 = psi_h_ecmwf(z0t*Linv) 246 298 ztmp0 = func_h - ztmp2 247 ztmp1 = vkarmn/(LOG(zu +z0) - LOG(z0t) - ztmp0)299 ztmp1 = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) 248 300 t_star = dt_zu*ztmp1 249 301 ztmp2 = ztmp0 - func_m + ztmp2 … … 253 305 ztmp2 = psi_h_ecmwf(z0q*Linv) 254 306 ztmp0 = func_h - ztmp2 255 ztmp1 = vkarmn/(LOG(zu +z0) - LOG(z0q) - ztmp0)307 ztmp1 = vkarmn/(LOG(zu) - LOG(z0q) - ztmp0) 256 308 q_star = dq_zu*ztmp1 257 309 ztmp2 = ztmp0 - func_m + ztmp2 258 ztmp1 = log(zt/zu) + ztmp2310 ztmp1 = LOG(zt/zu) + ztmp2 259 311 q_zu = q_zt - q_star/vkarmn*ztmp1 260 261 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )262 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )263 264 312 END IF 265 313 266 314 !! Updating because of updated z0 and z0t and new Linv... 267 ztmp1 = zu + z0 268 ztmp0 = ztmp1*Linv 269 func_m = log(ztmp1) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv) 270 func_h = log(ztmp1) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 271 272 END DO 315 ztmp0 = zu*Linv 316 func_m = log(zu) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv) 317 func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 318 319 !! SKIN related part 320 !! ----------------- 321 IF( l_use_skin ) THEN 322 !! compute transfer coefficients at zu : lolo: verifier... 323 Cd = vkarmn*vkarmn/(func_m*func_m) 324 Ch = vkarmn*vkarmn/(func_m*func_h) 325 ztmp1 = LOG(zu) - LOG(z0q) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0q*Linv) ! func_q 326 Ce = vkarmn*vkarmn/(func_m*ztmp1) 327 ! Non-Solar heat flux to the ocean: 328 ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp) ! rho*U10 329 ztmp2 = T_s*T_s 330 ztmp1 = ztmp1 * ( Ce*rLevap*(q_zu - q_s) + Ch*rCp_dry*(t_zu - T_s) ) & ! Total turb. heat flux 331 & + (rad_lw - emiss_w*stefan*ztmp2*ztmp2) ! Net longwave flux 332 !! Updating the values of the skin temperature T_s and q_s : 333 CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) 334 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! 200 -> just to avoid numerics problem on masked regions if silly values are given 335 END IF 336 337 IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN 338 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 339 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 340 END IF 341 342 END DO !DO j_itt = 1, nb_itt 273 343 274 344 Cd = vkarmn*vkarmn/(func_m*func_m) 275 345 Ch = vkarmn*vkarmn/(func_m*func_h) 276 ztmp1 = log((zu + z0)/z0q) - psi_h_ecmwf((zu + z0)*Linv) + psi_h_ecmwf(z0q*Linv) ! func_q 277 Ce = vkarmn*vkarmn/(func_m*ztmp1) 278 279 ztmp1 = zu + z0 280 Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 )) 281 Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t)) 282 Cen = vkarmn*vkarmn / (log(ztmp1/z0q)*log(ztmp1/z0q)) 346 ztmp2 = log(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv) ! func_q 347 Ce = vkarmn*vkarmn/(func_m*ztmp2) 348 349 Cdn = vkarmn*vkarmn / (log(zu/z0 )*log(zu/z0 )) 350 Chn = vkarmn*vkarmn / (log(zu/z0t)*log(zu/z0t)) 351 Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q)) 283 352 284 353 END SUBROUTINE TURB_ECMWF … … 294 363 !! and L is M-O length 295 364 !! 296 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)365 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 297 366 !!---------------------------------------------------------------------------------- 298 367 REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf … … 306 375 DO ji = 1, jpi 307 376 ! 308 zzeta = MIN( pzeta(ji,jj) , 5. ) !! Very stable conditions (L positif and big!):377 zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 309 378 ! 310 379 ! Unstable (Paulson 1970): 311 380 ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 312 zx = SQRT(ABS(1. - 16.*zzeta))313 ztmp = 1. + SQRT(zx)381 zx = SQRT(ABS(1._wp - 16._wp*zzeta)) 382 ztmp = 1._wp + SQRT(zx) 314 383 ztmp = ztmp*ztmp 315 psi_unst = LOG( 0.125 *ztmp*(1.+ zx) ) &316 & -2. *ATAN( SQRT(zx) ) + 0.5*rpi384 psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) ) & 385 & -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi 317 386 ! 318 387 ! Unstable: 319 388 ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 320 psi_stab = -2. /3.*(zzeta - 5./0.35)*EXP(-0.35*zzeta) &321 & - zzeta - 2. /3.*5./0.35389 psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & 390 & - zzeta - 2._wp/3._wp*5._wp/0.35_wp 322 391 ! 323 392 ! Combining: 324 stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1325 ! 326 psi_m_ecmwf(ji,jj) = (1. - stab) * psi_unst & ! (zzeta < 0) Unstable327 & + stab * psi_stab ! (zzeta > 0) Stable393 stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 394 ! 395 psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 396 & + stab * psi_stab ! (zzeta > 0) Stable 328 397 ! 329 398 END DO … … 332 401 END FUNCTION psi_m_ecmwf 333 402 334 403 335 404 FUNCTION psi_h_ecmwf( pzeta ) 336 405 !!---------------------------------------------------------------------------------- … … 342 411 !! and L is M-O length 343 412 !! 344 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)413 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 345 414 !!---------------------------------------------------------------------------------- 346 415 REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf … … 354 423 DO ji = 1, jpi 355 424 ! 356 zzeta = MIN(pzeta(ji,jj) , 5. ) ! Very stable conditions (L positif and big!):357 ! 358 zx = ABS(1. - 16.*zzeta)**.25 ! this is actually (1/phi_m)**2 !!!425 zzeta = MIN(pzeta(ji,jj) , 5._wp) ! Very stable conditions (L positif and big!): 426 ! 427 zx = ABS(1._wp - 16._wp*zzeta)**.25 ! this is actually (1/phi_m)**2 !!! 359 428 ! ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 360 429 ! Unstable (Paulson 1970) : 361 psi_unst = 2. *LOG(0.5*(1.+ zx*zx)) ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1430 psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx)) ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 362 431 ! 363 432 ! Stable: 364 psi_stab = -2. /3.*(zzeta - 5./0.35)*EXP(-0.35*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1365 & - ABS(1. + 2./3.*zzeta)**1.5 - 2./3.*5./0.35 + 1.433 psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 434 & - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp 366 435 ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 367 436 ! 368 stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1369 ! 370 ! 371 psi_h_ecmwf(ji,jj) = (1. - stab) * psi_unst & ! (zzeta < 0) Unstable372 & + stab * psi_stab ! (zzeta > 0) Stable437 stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 438 ! 439 ! 440 psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 441 & + stab * psi_stab ! (zzeta > 0) Stable 373 442 ! 374 443 END DO … … 382 451 !! Bulk Richardson number (Eq. 3.25 IFS doc) 383 452 !! 384 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)453 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 385 454 !!---------------------------------------------------------------------------------- 386 455 REAL(wp), DIMENSION(jpi,jpj) :: Ri_bulk ! … … 394 463 !!---------------------------------------------------------------------------------- 395 464 ! 396 Ri_bulk = grav*pz/(pub*pub) 397 & * ( pdt/(ptz - 0.5_wp*(pdt + grav*pz/( Cp_dry+Cp_vap*pqz)))&465 Ri_bulk = grav*pz/(pub*pub) & 466 & * ( pdt/(ptz - 0.5_wp*(pdt + grav*pz/(rCp_dry + rCp_vap*pqz))) & 398 467 & + rctv0*pdq ) 399 468 ! … … 401 470 402 471 403 FUNCTION visc_air(ptak)404 !!----------------------------------------------------------------------------------405 !! Air kinetic viscosity (m^2/s) given from temperature in degrees...406 !!407 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)408 !!----------------------------------------------------------------------------------409 REAL(wp), DIMENSION(jpi,jpj) :: visc_air !410 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature in (K)411 !412 INTEGER :: ji, jj ! dummy loop indices413 REAL(wp) :: ztc, ztc2 ! local scalar414 !!----------------------------------------------------------------------------------415 !416 DO jj = 1, jpj417 DO ji = 1, jpi418 ztc = ptak(ji,jj) - rt0 ! air temp, in deg. C419 ztc2 = ztc*ztc420 visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc)421 END DO422 END DO423 !424 END FUNCTION visc_air425 472 426 473 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
