Changeset 12928 for NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/SBC/sbcblk_algo_ecmwf.F90
- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/SBC/sbcblk_algo_ecmwf.F90
r10069 r12928 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.F90 48 49 ! !! ECMWF own values for given constants, taken form IFS documentation... 46 PUBLIC :: SBCBLK_ALGO_ECMWF_INIT, TURB_ECMWF 47 !! * Substitutions 48 # include "do_loop_substitute.h90" 49 50 !! ECMWF own values for given constants, taken form IFS documentation... 50 51 REAL(wp), PARAMETER :: charn0 = 0.018 ! Charnock constant (pretty high value here !!! 51 52 ! ! => Usually 0.011 for moderate winds) 52 53 REAL(wp), PARAMETER :: zi0 = 1000. ! scale height of the atmospheric boundary layer...1 53 54 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 55 REAL(wp), PARAMETER :: alpha_M = 0.11 ! For roughness length (smooth surface term) 58 56 REAL(wp), PARAMETER :: alpha_H = 0.40 ! (Chapter 3, p.34, IFS doc Cy31r1) 59 57 REAL(wp), PARAMETER :: alpha_Q = 0.62 ! 58 59 INTEGER , PARAMETER :: nb_itt = 10 ! number of itterations 60 60 61 !!---------------------------------------------------------------------- 61 62 CONTAINS 62 63 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 64 65 SUBROUTINE sbcblk_algo_ecmwf_init(l_use_cs, l_use_wl) 66 !!--------------------------------------------------------------------- 67 !! *** FUNCTION sbcblk_algo_ecmwf_init *** 76 68 !! 77 69 !! INPUT : 78 70 !! ------- 71 !! * l_use_cs : use the cool-skin parameterization 72 !! * l_use_wl : use the warm-layer parameterization 73 !!--------------------------------------------------------------------- 74 LOGICAL , INTENT(in) :: l_use_cs ! use the cool-skin parameterization 75 LOGICAL , INTENT(in) :: l_use_wl ! use the warm-layer parameterization 76 INTEGER :: ierr 77 !!--------------------------------------------------------------------- 78 IF( l_use_wl ) THEN 79 ierr = 0 80 ALLOCATE ( dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr ) 81 IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_wl & Hz_wl failed!' ) 82 dT_wl(:,:) = 0._wp 83 Hz_wl(:,:) = rd0 ! (rd0, constant, = 3m is default for Zeng & Beljaars) 84 ENDIF 85 IF( l_use_cs ) THEN 86 ierr = 0 87 ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 88 IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_cs failed!' ) 89 dT_cs(:,:) = -0.25_wp ! First guess of skin correction 90 ENDIF 91 END SUBROUTINE sbcblk_algo_ecmwf_init 92 93 94 95 SUBROUTINE turb_ecmwf( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 96 & Cd, Ch, Ce, t_zu, q_zu, U_blk, & 97 & Cdn, Chn, Cen, & 98 & Qsw, rad_lw, slp, pdT_cs, & ! optionals for cool-skin (and warm-layer) 99 & pdT_wl, pHz_wl ) ! optionals for warm-layer only 100 !!---------------------------------------------------------------------------------- 101 !! *** ROUTINE turb_ecmwf *** 102 !! 103 !! ** Purpose : Computes turbulent transfert coefficients of surface 104 !! fluxes according to IFS doc. (cycle 45r1) 105 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 106 !! Returns the effective bulk wind speed at zu to be used in the bulk formulas 107 !! 108 !! Applies the cool-skin warm-layer correction of the SST to T_s 109 !! if the net shortwave flux at the surface (Qsw), the downwelling longwave 110 !! radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp) 111 !! are provided as (optional) arguments! 112 !! 113 !! INPUT : 114 !! ------- 115 !! * kt : current time step (starts at 1) 79 116 !! * 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] 117 !! * zu : height for wind speed (usually 10m) [m] 83 118 !! * t_zt : potential air temperature at zt [K] 84 !! * ssq : specific humidity at saturation at SST [kg/kg]85 119 !! * q_zt : specific humidity of air at zt [kg/kg] 86 !! 120 !! * U_zu : scalar wind speed at zu [m/s] 121 !! * l_use_cs : use the cool-skin parameterization 122 !! * l_use_wl : use the warm-layer parameterization 123 !! 124 !! INPUT/OUTPUT: 125 !! ------------- 126 !! * T_s : always "bulk SST" as input [K] 127 !! -> unchanged "bulk SST" as output if CSWL not used [K] 128 !! -> skin temperature as output if CSWL used [K] 129 !! 130 !! * q_s : SSQ aka saturation specific humidity at temp. T_s [kg/kg] 131 !! -> doesn't need to be given a value if skin temp computed (in case l_use_cs=True or l_use_wl=True) 132 !! -> MUST be given the correct value if not computing skint temp. (in case l_use_cs=False or l_use_wl=False) 133 !! 134 !! OPTIONAL INPUT: 135 !! --------------- 136 !! * Qsw : net solar flux (after albedo) at the surface (>0) [W/m^2] 137 !! * rad_lw : downwelling longwave radiation at the surface (>0) [W/m^2] 138 !! * slp : sea-level pressure [Pa] 139 !! 140 !! OPTIONAL OUTPUT: 141 !! ---------------- 142 !! * pdT_cs : SST increment "dT" for cool-skin correction [K] 143 !! * pdT_wl : SST increment "dT" for warm-layer correction [K] 144 !! * pHz_wl : thickness of warm-layer [m] 87 145 !! 88 146 !! OUTPUT : … … 93 151 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 94 152 !! * 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 !!---------------------------------------------------------------------------------- 153 !! * U_blk : bulk wind speed at zu [m/s] 154 !! 155 !! 156 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 157 !!---------------------------------------------------------------------------------- 158 INTEGER, INTENT(in ) :: kt ! current time step 100 159 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] 101 160 REAL(wp), INTENT(in ) :: zu ! height for U_zu [m] 102 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: sst! sea surface temperature [Kelvin]161 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: T_s ! sea surface temperature [Kelvin] 103 162 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 163 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: q_s ! sea surface specific humidity [kg/kg] 164 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg] 106 165 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s] 166 LOGICAL , INTENT(in ) :: l_use_cs ! use the cool-skin parameterization 167 LOGICAL , INTENT(in ) :: l_use_wl ! use the warm-layer parameterization 107 168 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau) 108 169 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ch ! transfer coefficient for sensible heat (Q_sens) … … 110 171 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K] 111 172 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]173 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind speed at zu [m/s] 113 174 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients 114 175 ! 176 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: Qsw ! [W/m^2] 177 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: rad_lw ! [W/m^2] 178 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: slp ! [Pa] 179 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_cs 180 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_wl ! [K] 181 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pHz_wl ! [m] 182 ! 115 183 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 ! 130 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 184 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 185 ! 186 REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star 187 REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu 188 REAL(wp), DIMENSION(jpi,jpj) :: znu_a !: Nu_air, Viscosity of air 189 REAL(wp), DIMENSION(jpi,jpj) :: Linv !: 1/L (inverse of Monin Obukhov length... 190 REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t, z0q 191 ! 192 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst ! to back up the initial bulk SST 193 ! 194 REAL(wp), DIMENSION(jpi,jpj) :: func_m, func_h 195 REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 196 CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90' 197 !!---------------------------------------------------------------------------------- 198 IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl) 199 200 l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! 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 DO jj = 1, jpj 306 DO ji = 1, jpi 307 ! 308 zzeta = MIN( pzeta(ji,jj) , 5. ) !! Very stable conditions (L positif and big!): 309 ! 310 ! Unstable (Paulson 1970): 311 ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 312 zx = SQRT(ABS(1. - 16.*zzeta)) 313 ztmp = 1. + SQRT(zx) 314 ztmp = ztmp*ztmp 315 psi_unst = LOG( 0.125*ztmp*(1. + zx) ) & 316 & -2.*ATAN( SQRT(zx) ) + 0.5*rpi 317 ! 318 ! Unstable: 319 ! 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.35 322 ! 323 ! Combining: 324 stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1 325 ! 326 psi_m_ecmwf(ji,jj) = (1. - stab) * psi_unst & ! (zzeta < 0) Unstable 327 & + stab * psi_stab ! (zzeta > 0) Stable 328 ! 329 END DO 330 END DO 331 ! 412 DO_2D_11_11 413 ! 414 zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 415 ! 416 ! Unstable (Paulson 1970): 417 ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 418 zx = SQRT(ABS(1._wp - 16._wp*zzeta)) 419 ztmp = 1._wp + SQRT(zx) 420 ztmp = ztmp*ztmp 421 psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) ) & 422 & -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi 423 ! 424 ! Unstable: 425 ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 426 psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & 427 & - zzeta - 2._wp/3._wp*5._wp/0.35_wp 428 ! 429 ! Combining: 430 stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 431 ! 432 psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 433 & + stab * psi_stab ! (zzeta > 0) Stable 434 ! 435 END_2D 332 436 END FUNCTION psi_m_ecmwf 333 437 334 438 335 439 FUNCTION psi_h_ecmwf( pzeta ) 336 440 !!---------------------------------------------------------------------------------- … … 342 446 !! and L is M-O length 343 447 !! 344 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)448 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 345 449 !!---------------------------------------------------------------------------------- 346 450 REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf … … 351 455 !!---------------------------------------------------------------------------------- 352 456 ! 353 DO jj = 1, jpj 354 DO ji = 1, jpi 355 ! 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 !!! 359 ! ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 360 ! Unstable (Paulson 1970) : 361 psi_unst = 2.*LOG(0.5*(1. + zx*zx)) ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 362 ! 363 ! Stable: 364 psi_stab = -2./3.*(zzeta - 5./0.35)*EXP(-0.35*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 365 & - ABS(1. + 2./3.*zzeta)**1.5 - 2./3.*5./0.35 + 1. 366 ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 367 ! 368 stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1 369 ! 370 ! 371 psi_h_ecmwf(ji,jj) = (1. - stab) * psi_unst & ! (zzeta < 0) Unstable 372 & + stab * psi_stab ! (zzeta > 0) Stable 373 ! 374 END DO 375 END DO 376 ! 457 DO_2D_11_11 458 ! 459 zzeta = MIN(pzeta(ji,jj) , 5._wp) ! Very stable conditions (L positif and big!): 460 ! 461 zx = ABS(1._wp - 16._wp*zzeta)**.25 ! this is actually (1/phi_m)**2 !!! 462 ! ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 463 ! Unstable (Paulson 1970) : 464 psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx)) ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 465 ! 466 ! Stable: 467 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 468 & - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp 469 ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 470 ! 471 stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 472 ! 473 ! 474 psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 475 & + stab * psi_stab ! (zzeta > 0) Stable 476 ! 477 END_2D 377 478 END FUNCTION psi_h_ecmwf 378 479 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 480 426 481 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.