Changeset 12154 for NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/sbcblk_algo_ecmwf.F90
- Timestamp:
- 2019-12-10T15:44:23+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/sbcblk_algo_ecmwf.F90
r10069 r12154 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 2019 / 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_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB 41 USE sbcblk_skin_ecmwf ! cool-skin/warm layer scheme !LB 43 42 44 43 IMPLICIT NONE 45 44 PRIVATE 46 45 47 PUBLIC :: TURB_ECMWF ! called by sbcblk.F9048 49 ! !! ECMWF own values for given constants, taken form IFS documentation...46 PUBLIC :: SBCBLK_ALGO_ECMWF_INIT, TURB_ECMWF 47 48 !! ECMWF own values for given constants, taken form IFS documentation... 50 49 REAL(wp), PARAMETER :: charn0 = 0.018 ! Charnock constant (pretty high value here !!! 51 50 ! ! => Usually 0.011 for moderate winds) 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 = 10 ! number of itterations 58 60 59 !!---------------------------------------------------------------------- 61 60 CONTAINS 62 61 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 ) 66 !!---------------------------------------------------------------------------------- 67 !! *** ROUTINE turb_ecmwf *** 68 !! 69 !! 2015: L. Brodeau (brodeau@gmail.com) 70 !! 71 !! ** Purpose : Computes turbulent transfert coefficients of surface 72 !! fluxes according to IFS doc. (cycle 31) 73 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 74 !! 75 !! ** Method : Monin Obukhov Similarity Theory 62 63 SUBROUTINE sbcblk_algo_ecmwf_init(l_use_cs, l_use_wl) 64 !!--------------------------------------------------------------------- 65 !! *** FUNCTION sbcblk_algo_ecmwf_init *** 76 66 !! 77 67 !! INPUT : 78 68 !! ------- 69 !! * l_use_cs : use the cool-skin parameterization 70 !! * l_use_wl : use the warm-layer parameterization 71 !!--------------------------------------------------------------------- 72 LOGICAL , INTENT(in) :: l_use_cs ! use the cool-skin parameterization 73 LOGICAL , INTENT(in) :: l_use_wl ! use the warm-layer parameterization 74 INTEGER :: ierr 75 !!--------------------------------------------------------------------- 76 IF( l_use_wl ) THEN 77 ierr = 0 78 ALLOCATE ( dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr ) 79 IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_wl & Hz_wl failed!' ) 80 dT_wl(:,:) = 0._wp 81 Hz_wl(:,:) = rd0 ! (rd0, constant, = 3m is default for Zeng & Beljaars) 82 ENDIF 83 IF( l_use_cs ) THEN 84 ierr = 0 85 ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 86 IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_cs failed!' ) 87 dT_cs(:,:) = -0.25_wp ! First guess of skin correction 88 ENDIF 89 END SUBROUTINE sbcblk_algo_ecmwf_init 90 91 92 93 SUBROUTINE turb_ecmwf( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 94 & Cd, Ch, Ce, t_zu, q_zu, U_blk, & 95 & Cdn, Chn, Cen, & 96 & Qsw, rad_lw, slp, pdT_cs, & ! optionals for cool-skin (and warm-layer) 97 & pdT_wl, pHz_wl ) ! optionals for warm-layer only 98 !!---------------------------------------------------------------------- 99 !! *** ROUTINE turb_ecmwf *** 100 !! 101 !! ** Purpose : Computes turbulent transfert coefficients of surface 102 !! fluxes according to IFS doc. (cycle 45r1) 103 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 104 !! Returns the effective bulk wind speed at zu to be used in the bulk formulas 105 !! 106 !! Applies the cool-skin warm-layer correction of the SST to T_s 107 !! if the net shortwave flux at the surface (Qsw), the downwelling longwave 108 !! radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp) 109 !! are provided as (optional) arguments! 110 !! 111 !! INPUT : 112 !! ------- 113 !! * kt : current time step (starts at 1) 79 114 !! * zt : height for temperature and spec. hum. of air [m] 80 !! * zu : height for wind speed (generally 10m) [m] 81 !! * U_zu : scalar wind speed at 10m [m/s] 82 !! * sst : SST [K] 115 !! * zu : height for wind speed (usually 10m) [m] 83 116 !! * t_zt : potential air temperature at zt [K] 84 !! * ssq : specific humidity at saturation at SST [kg/kg]85 117 !! * q_zt : specific humidity of air at zt [kg/kg] 86 !! 118 !! * U_zu : scalar wind speed at zu [m/s] 119 !! * l_use_cs : use the cool-skin parameterization 120 !! * l_use_wl : use the warm-layer parameterization 121 !! 122 !! INPUT/OUTPUT: 123 !! ------------- 124 !! * T_s : always "bulk SST" as input [K] 125 !! -> unchanged "bulk SST" as output if CSWL not used [K] 126 !! -> skin temperature as output if CSWL used [K] 127 !! 128 !! * q_s : SSQ aka saturation specific humidity at temp. T_s [kg/kg] 129 !! -> doesn't need to be given a value if skin temp computed (in case l_use_cs=True or l_use_wl=True) 130 !! -> MUST be given the correct value if not computing skint temp. (in case l_use_cs=False or l_use_wl=False) 131 !! 132 !! OPTIONAL INPUT: 133 !! --------------- 134 !! * Qsw : net solar flux (after albedo) at the surface (>0) [W/m^2] 135 !! * rad_lw : downwelling longwave radiation at the surface (>0) [W/m^2] 136 !! * slp : sea-level pressure [Pa] 137 !! 138 !! OPTIONAL OUTPUT: 139 !! ---------------- 140 !! * pdT_cs : SST increment "dT" for cool-skin correction [K] 141 !! * pdT_wl : SST increment "dT" for warm-layer correction [K] 142 !! * pHz_wl : thickness of warm-layer [m] 87 143 !! 88 144 !! OUTPUT : … … 93 149 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 94 150 !! * 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) 99 !!---------------------------------------------------------------------------------- 151 !! * U_blk : bulk wind speed at zu [m/s] 152 !! 153 !! 154 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 155 !!---------------------------------------------------------------------------------- 156 INTEGER, INTENT(in ) :: kt ! current time step 100 157 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] 101 158 REAL(wp), INTENT(in ) :: zu ! height for U_zu [m] 102 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: sst! sea surface temperature [Kelvin]159 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: T_s ! sea surface temperature [Kelvin] 103 160 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 161 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: q_s ! sea surface specific humidity [kg/kg] 162 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg] 106 163 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s] 164 LOGICAL , INTENT(in ) :: l_use_cs ! use the cool-skin parameterization 165 LOGICAL , INTENT(in ) :: l_use_wl ! use the warm-layer parameterization 107 166 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau) 108 167 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ch ! transfer coefficient for sensible heat (Q_sens) … … 110 169 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K] 111 170 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]171 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind speed at zu [m/s] 113 172 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients 114 173 ! 174 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: Qsw ! [W/m^2] 175 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: rad_lw ! [W/m^2] 176 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: slp ! [Pa] 177 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_cs 178 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_wl ! [K] 179 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pHz_wl ! [m] 180 ! 115 181 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 118 ! 119 REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star, & 120 & dt_zu, dq_zu, & 121 & znu_a, & !: Nu_air, Viscosity of air 122 & Linv, & !: 1/L (inverse of Monin Obukhov length... 123 & z0, z0t, z0q 124 REAL(wp), DIMENSION(jpi,jpj) :: func_m, func_h 125 REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 126 !!---------------------------------------------------------------------------------- 127 ! 128 ! Identical first gess as in COARE, with IFS parameter values though 129 ! 182 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 183 ! 184 REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star 185 REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu 186 REAL(wp), DIMENSION(jpi,jpj) :: znu_a !: Nu_air, Viscosity of air 187 REAL(wp), DIMENSION(jpi,jpj) :: Linv !: 1/L (inverse of Monin Obukhov length... 188 REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t, z0q 189 ! 190 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst ! to back up the initial bulk SST 191 ! 192 REAL(wp), DIMENSION(jpi,jpj) :: func_m, func_h 193 REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 194 CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90' 195 !!---------------------------------------------------------------------------------- 196 197 IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl) 198 130 199 l_zt_equal_zu = .FALSE. 131 IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 132 133 200 IF( ABS(zu - zt) < 0.01_wp ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 201 202 !! Initializations for cool skin and warm layer: 203 IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 204 & CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' ) 205 206 IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 207 & CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' ) 208 209 IF( l_use_cs .OR. l_use_wl ) THEN 210 ALLOCATE ( zsst(jpi,jpj) ) 211 zsst = T_s ! backing up the bulk SST 212 IF( l_use_cs ) T_s = T_s - 0.25_wp ! First guess of correction 213 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 214 ENDIF 215 216 217 ! Identical first gess as in COARE, with IFS parameter values though... 218 ! 134 219 !! 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 ) ! "220 t_zu = MAX( t_zt , 180._wp ) ! who knows what's given on masked-continental regions... 221 q_zu = MAX( q_zt , 1.e-6_wp ) ! " 137 222 138 223 !! 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 )141 142 znu_a = visc_air(t_z t) ! Air viscosity (m^2/s) at zt given from temperature in (K)143 144 ztmp2 = 0.5 * 0.5! initial guess for wind gustiness contribution145 U_blk = SQRT(U_zu*U_zu + ztmp2) 146 147 ! z0 = 0.0001148 ztmp2 = 10000. ! optimization: ztmp2 == 1/z0149 ztmp0 = LOG(zu*ztmp2) 150 z tmp1 = 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_star154 z0t = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ! WARNING: 1/z0t !224 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 225 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 226 227 znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 228 229 U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution 230 231 ztmp0 = LOG( zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001) 232 ztmp1 = LOG(10._wp*10000._wp) ! " " " 233 u_star = 0.035_wp*U_blk*ztmp1/ztmp0 ! (u* = 0.035*Un10) 234 235 z0 = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star 236 z0 = MIN( MAX(ABS(z0), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 237 238 z0t = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 239 z0t = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 155 240 156 241 Cd = (vkarmn/ztmp0)**2 ! first guess of Cd 157 242 158 ztmp0 = vkarmn*vkarmn/LOG(zt *z0t)/Cd159 160 ztmp2 = Ri_bulk( zu, t_zu, dt_zu, q_zu, dq_zu, U_blk ) ! Ribu = Bulk Richardson number161 162 !! First estimate of zeta_u, depending on the stability, ie sign of Ribu(ztmp2):163 ztmp1 = 0.5 + SIGN( 0.5 , ztmp2 )243 ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 244 245 ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 246 247 !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 248 ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 164 249 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))250 func_h = (1._wp-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & ! BRN < 0 ! temporary array !!! func_h == zeta_u 251 & + ztmp1 * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m)) ! BRN > 0 252 !#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 253 169 254 !! 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))171 172 u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_ecmwf(func_h))255 ztmp0 = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h)) 256 257 u_star = MAX ( U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_ecmwf(func_h)) , 1.E-9 ) ! (MAX => prevents FPE from stupid values from masked region later on) 173 258 t_star = dt_zu*ztmp0 174 259 q_star = dq_zu*ztmp0 175 260 176 ! What 's needto be done if zt /= zu:261 ! What needs to be done if zt /= zu: 177 262 IF( .NOT. l_zt_equal_zu ) THEN 178 !179 263 !! First update of values at zu (or zt for wind) 180 264 ztmp0 = psi_h_ecmwf(func_h) - psi_h_ecmwf(zt*func_h/zu) ! zt*func_h/zu == zeta_t 181 ztmp1 = log(zt/zu) + ztmp0265 ztmp1 = LOG(zt/zu) + ztmp0 182 266 t_zu = t_zt - t_star/vkarmn*ztmp1 183 267 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 ) 268 q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity : 188 269 ! 270 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 271 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 189 272 ENDIF 190 273 … … 194 277 195 278 !! First guess of inverse of Monin-Obukov length (1/L) : 196 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 ) 279 Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star ) 198 280 199 281 !! 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 282 ztmp0 = zu*Linv 283 func_m = LOG(zu) - LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv) 284 func_h = LOG(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 205 285 206 286 !! ITERATION BLOCK 207 !! ***************208 209 287 DO j_itt = 1, nb_itt 210 288 211 289 !! Bulk Richardson Number at z=zu (Eq. 3.25) 212 ztmp0 = Ri_bulk( zu, t_zu, dt_zu, q_zu, dq_zu, U_blk)290 ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 213 291 214 292 !! 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 293 Linv = ztmp0*func_m*func_m/func_h / zu ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1 294 !! Note: it is slightly different that the L we would get with the usual 295 Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) 216 296 217 297 !! 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) 298 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 299 221 300 !! Need to update roughness lengthes: … … 223 302 ztmp2 = u_star*u_star 224 303 ztmp1 = znu_a/u_star 225 z0 = alpha_M*ztmp1 + charn0*ztmp2/grav 226 z0t = alpha_H*ztmp1 ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 227 z0q = alpha_Q*ztmp1 228 229 !! Update wind at 10m taking into acount convection-related wind gustiness: 230 ! 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) 232 !! => 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 - Cy31r1 304 z0 = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp) 305 z0t = MIN( ABS( alpha_H*ztmp1 ) , 0.001_wp) ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 306 z0q = MIN( ABS( alpha_Q*ztmp1 ) , 0.001_wp) 307 308 !! Update wind at zu with convection-related wind gustiness in unstable conditions (Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8) 309 ztmp2 = Beta0*Beta0*ztmp2*(MAX(-zi0*Linv/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1) 310 !! ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0 311 U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp) ! include gustiness in bulk wind speed 234 312 ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 235 313 … … 238 316 !! as well the air-sea differences: 239 317 IF( .NOT. l_zt_equal_zu ) THEN 240 241 318 !! 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 !!!319 func_h = psi_h_ecmwf(zu*Linv) ! temporary array !!! 320 func_m = psi_h_ecmwf(zt*Linv) ! temporary array !!! 244 321 245 322 ztmp2 = psi_h_ecmwf(z0t*Linv) 246 323 ztmp0 = func_h - ztmp2 247 ztmp1 = vkarmn/(LOG(zu +z0) - LOG(z0t) - ztmp0)324 ztmp1 = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) 248 325 t_star = dt_zu*ztmp1 249 326 ztmp2 = ztmp0 - func_m + ztmp2 … … 253 330 ztmp2 = psi_h_ecmwf(z0q*Linv) 254 331 ztmp0 = func_h - ztmp2 255 ztmp1 = vkarmn/(LOG(zu +z0) - LOG(z0q) - ztmp0)332 ztmp1 = vkarmn/(LOG(zu) - LOG(z0q) - ztmp0) 256 333 q_star = dq_zu*ztmp1 257 334 ztmp2 = ztmp0 - func_m + ztmp2 258 ztmp1 = log(zt/zu) + ztmp2335 ztmp1 = LOG(zt/zu) + ztmp2 259 336 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 END IF 337 ENDIF 265 338 266 339 !! 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 340 ztmp0 = zu*Linv 341 func_m = log(zu) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv) 342 func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 343 344 345 IF( l_use_cs ) THEN 346 !! Cool-skin contribution 347 348 CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 349 & ztmp1, ztmp0, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp0 350 351 CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst ) ! Qnsol -> ztmp1 352 353 T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1) 354 IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1) 355 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 356 357 ENDIF 358 359 IF( l_use_wl ) THEN 360 !! Warm-layer contribution 361 CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 362 & ztmp1, ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp2 363 CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst ) 364 !! Updating T_s and q_s !!! 365 T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1) ! 366 IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 367 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 368 ENDIF 369 370 IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 371 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 372 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 373 ENDIF 374 375 END DO !DO j_itt = 1, nb_itt 273 376 274 377 Cd = vkarmn*vkarmn/(func_m*func_m) 275 378 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)) 283 284 END SUBROUTINE TURB_ECMWF 379 ztmp2 = log(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv) ! func_q 380 Ce = vkarmn*vkarmn/(func_m*ztmp2) 381 382 Cdn = vkarmn*vkarmn / (log(zu/z0 )*log(zu/z0 )) 383 Chn = vkarmn*vkarmn / (log(zu/z0t)*log(zu/z0t)) 384 Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q)) 385 386 IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 387 IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 388 IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 389 390 IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 391 392 END SUBROUTINE turb_ecmwf 285 393 286 394 … … 294 402 !! and L is M-O length 295 403 !! 296 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)404 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 297 405 !!---------------------------------------------------------------------------------- 298 406 REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf … … 302 410 REAL(wp) :: zzeta, zx, ztmp, psi_unst, psi_stab, stab 303 411 !!---------------------------------------------------------------------------------- 304 !305 412 DO jj = 1, jpj 306 413 DO ji = 1, jpi 307 414 ! 308 zzeta = MIN( pzeta(ji,jj) , 5. ) !! Very stable conditions (L positif and big!):415 zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 309 416 ! 310 417 ! Unstable (Paulson 1970): 311 418 ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 312 zx = SQRT(ABS(1. - 16.*zzeta))313 ztmp = 1. + SQRT(zx)419 zx = SQRT(ABS(1._wp - 16._wp*zzeta)) 420 ztmp = 1._wp + SQRT(zx) 314 421 ztmp = ztmp*ztmp 315 psi_unst = LOG( 0.125 *ztmp*(1.+ zx) ) &316 & -2. *ATAN( SQRT(zx) ) + 0.5*rpi422 psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) ) & 423 & -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi 317 424 ! 318 425 ! Unstable: 319 426 ! 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.35427 psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & 428 & - zzeta - 2._wp/3._wp*5._wp/0.35_wp 322 429 ! 323 430 ! 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) Stable431 stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 432 ! 433 psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 434 & + stab * psi_stab ! (zzeta > 0) Stable 328 435 ! 329 436 END DO 330 437 END DO 331 !332 438 END FUNCTION psi_m_ecmwf 333 439 334 440 335 441 FUNCTION psi_h_ecmwf( pzeta ) 336 442 !!---------------------------------------------------------------------------------- … … 342 448 !! and L is M-O length 343 449 !! 344 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)450 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 345 451 !!---------------------------------------------------------------------------------- 346 452 REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf … … 354 460 DO ji = 1, jpi 355 461 ! 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 !!!462 zzeta = MIN(pzeta(ji,jj) , 5._wp) ! Very stable conditions (L positif and big!): 463 ! 464 zx = ABS(1._wp - 16._wp*zzeta)**.25 ! this is actually (1/phi_m)**2 !!! 359 465 ! ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 360 466 ! Unstable (Paulson 1970) : 361 psi_unst = 2. *LOG(0.5*(1.+ zx*zx)) ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1467 psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx)) ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 362 468 ! 363 469 ! 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.470 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 471 & - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp 366 472 ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 367 473 ! 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) Stable474 stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 475 ! 476 ! 477 psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 478 & + stab * psi_stab ! (zzeta > 0) Stable 373 479 ! 374 480 END DO 375 481 END DO 376 !377 482 END FUNCTION psi_h_ecmwf 378 483 379 380 FUNCTION Ri_bulk( pz, ptz, pdt, pqz, pdq, pub )381 !!----------------------------------------------------------------------------------382 !! Bulk Richardson number (Eq. 3.25 IFS doc)383 !!384 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)385 !!----------------------------------------------------------------------------------386 REAL(wp), DIMENSION(jpi,jpj) :: Ri_bulk !387 !388 REAL(wp) , INTENT(in) :: pz ! height above the sea [m]389 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptz ! air temperature at pz m [K]390 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pdt ! ptz - sst [K]391 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqz ! air temperature at pz m [kg/kg]392 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pdq ! pqz - ssq [kg/kg]393 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub ! bulk wind speed [m/s]394 !!----------------------------------------------------------------------------------395 !396 Ri_bulk = grav*pz/(pub*pub) &397 & * ( pdt/(ptz - 0.5_wp*(pdt + grav*pz/(Cp_dry+Cp_vap*pqz))) &398 & + rctv0*pdq )399 !400 END FUNCTION Ri_bulk401 402 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 484 426 485 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.