Changeset 11215 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90
- Timestamp:
- 2019-07-04T12:34: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
r11209 r11215 14 14 !! (https://github.com/brodeau/aerobulk) 15 15 !! 16 !! ** Author: L. Brodeau, june 2016/ AeroBulk (https://github.com/brodeau/aerobulk)16 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 17 17 !!---------------------------------------------------------------------- 18 18 !! History : 4.0 ! 2016-02 (L.Brodeau) Original code … … 38 38 39 39 USE sbc_oce ! Surface boundary condition: ocean fields 40 USE sbcblk_phy ! LB: all thermodynamics functions, rho_air, q_sat, etc... #LB41 USE sbcblk_skin ! cool-skin/warm layer scheme (CSWL_ECMWF) 40 USE sbcblk_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB 41 USE sbcblk_skin ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB 42 42 43 43 IMPLICIT NONE … … 57 57 INTEGER , PARAMETER :: nb_itt = 5 ! number of itterations 58 58 59 !! Cool-Skin / Warm-Layer related parameters:60 REAL(wp), PARAMETER :: &61 & rdt0 = 3600.*1.5, & !: time step62 & rd0 = 3. , & !: Depth scale [m], "d" in Eq.11 (Zeng & Beljaars 2005)63 & rNu0 = 0.5 !: Nu (exponent of temperature profile) Eq.1164 ! !: (Zeng & Beljaars 2005) !: set to 0.5 instead of65 ! !: 0.3 to respect a warming of +3 K in calm66 ! !: condition for the insolation peak of +1000W/m^267 INTEGER, PARAMETER :: &68 & nb_itt_wl = 10 !: number of sub-itterations for solving the differential equation in warm-layer part69 ! !: => 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 59 !!---------------------------------------------------------------------- 75 60 CONTAINS … … 123 108 !! 124 109 !! 125 !! ** Author: L. Brodeau, june 2016/ AeroBulk (https://github.com/brodeau/aerobulk/)110 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 126 111 !!---------------------------------------------------------------------------------- 127 112 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] … … 157 142 LOGICAL :: l_use_skin = .FALSE. 158 143 REAL(wp), DIMENSION(jpi,jpj) :: zsst ! to back up the initial bulk SST 159 160 144 ! 161 145 REAL(wp), DIMENSION(jpi,jpj) :: func_m, func_h 162 146 REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 163 147 !!---------------------------------------------------------------------------------- 164 ! 148 165 149 ! Cool skin ? 166 150 IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN … … 169 153 IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin 170 154 171 ! Identical first gess as in COARE, with IFS parameter values though 155 ! Identical first gess as in COARE, with IFS parameter values though... 172 156 ! 173 157 l_zt_equal_zu = .FALSE. … … 182 166 183 167 !! First guess of temperature and humidity at height zu: 184 t_zu = MAX( t_zt , 0.0_wp) ! who knows what's given on masked-continental regions...185 q_zu = MAX( q_zt , 1.e-6_wp ) ! "168 t_zu = MAX( t_zt , 180._wp ) ! who knows what's given on masked-continental regions... 169 q_zu = MAX( q_zt , 1.e-6_wp ) ! " 186 170 187 171 !! Pot. temp. difference (and we don't want it to be 0!) … … 189 173 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 190 174 191 znu_a = visc_air(t_z t) ! Air viscosity (m^2/s) at zt given from temperature in (K)175 znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 192 176 193 177 ztmp2 = 0.5_wp*0.5_wp ! initial guess for wind gustiness contribution … … 214 198 ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 215 199 func_m = ztmp0*ztmp2 ! temporary array !! 216 func_h = (1. -ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & ! BRN < 0 ! temporary array !!! func_h == zeta_u217 & + ztmp1 * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m)) ! BRN > 0200 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 201 & + ztmp1 * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m)) ! BRN > 0 218 202 !#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" ! 219 203 … … 259 243 Linv = ztmp0*func_m*func_m/func_h / zu ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1 260 244 !! Note: it is slightly different that the L we would get with the usual 261 !! expression, as in coare algorithm or in 'mod_phy.f90' (One_on_L_MO())262 245 Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 263 246 … … 315 298 IF( l_use_skin ) THEN 316 299 !! compute transfer coefficients at zu : lolo: verifier... 317 Cd = vkarmn*vkarmn/(func_m*func_m)318 300 Ch = vkarmn*vkarmn/(func_m*func_h) 319 301 ztmp1 = LOG(zu) - LOG(z0q) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0q*Linv) ! func_q … … 357 339 !! and L is M-O length 358 340 !! 359 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)341 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 360 342 !!---------------------------------------------------------------------------------- 361 343 REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf … … 405 387 !! and L is M-O length 406 388 !! 407 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)389 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 408 390 !!---------------------------------------------------------------------------------- 409 391 REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf
Note: See TracChangeset
for help on using the changeset viewer.