Changeset 12154 for NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/sbcblk_algo_ncar.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_ncar.F90
r10190 r12154 11 11 !! 12 12 !! Routine turb_ncar maintained and developed in AeroBulk 13 !! (http ://aerobulk.sourceforge.net/)13 !! (https://github.com/brodeau/aerobulk/) 14 14 !! 15 15 !! L. Brodeau, 2015 … … 38 38 USE lib_fortran ! to use key_nosignedzero 39 39 40 USE sbcblk_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB 40 41 41 42 IMPLICIT NONE 42 43 PRIVATE 43 44 44 PUBLIC :: TURB_NCAR ! called by sbcblk.F90 45 46 ! ! NCAR own values for given constants: 47 REAL(wp), PARAMETER :: rctv0 = 0.608 ! constant to obtain virtual temperature... 48 45 PUBLIC :: TURB_NCAR ! called by sbcblk.F90 46 47 INTEGER , PARAMETER :: nb_itt = 5 ! number of itterations 48 49 49 !!---------------------------------------------------------------------- 50 50 CONTAINS … … 61 61 !! Returns the effective bulk wind speed at 10m to be used in the bulk formulas 62 62 !! 63 !! ** Method : Monin Obukhov Similarity Theory64 !! + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10)65 !!66 !! ** References : Large & Yeager, 2004 / Large & Yeager, 200867 !!68 !! ** Last update: Laurent Brodeau, June 2014:69 !! - handles both cases zt=zu and zt/=zu70 !! - optimized: less 2D arrays allocated and less operations71 !! - better first guess of stability by checking air-sea difference of virtual temperature72 !! rather than temperature difference only...73 !! - added function "cd_neutral_10m" that uses the improved parametrization of74 !! Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions!75 !! - using code-wide physical constants defined into "phycst.mod" rather than redifining them76 !! => 'vkarmn' and 'grav'77 !!78 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)79 63 !! 80 64 !! INPUT : 81 65 !! ------- 82 66 !! * zt : height for temperature and spec. hum. of air [m] 83 !! * zu : height for wind speed (generally 10m) [m] 84 !! * U_zu : scalar wind speed at 10m [m/s] 85 !! * sst : SST [K] 67 !! * zu : height for wind speed (usually 10m) [m] 68 !! * sst : bulk SST [K] 86 69 !! * t_zt : potential air temperature at zt [K] 87 70 !! * ssq : specific humidity at saturation at SST [kg/kg] 88 71 !! * q_zt : specific humidity of air at zt [kg/kg] 72 !! * U_zu : scalar wind speed at zu [m/s] 89 73 !! 90 74 !! … … 96 80 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 97 81 !! * q_zu : specific humidity of air // [kg/kg] 98 !! * U_blk : bulk wind at 10m [m/s] 82 !! * U_blk : bulk wind speed at zu [m/s] 83 !! 84 !! 85 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 99 86 !!---------------------------------------------------------------------------------- 100 87 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] … … 103 90 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin] 104 91 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 92 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg] 106 93 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s] 107 94 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau) … … 110 97 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K] 111 98 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]99 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind speed at zu [m/s] 113 100 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients 114 101 ! 115 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 102 INTEGER :: j_itt 103 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 118 104 ! 119 105 REAL(wp), DIMENSION(jpi,jpj) :: Cx_n10 ! 10m neutral latent/sensible coefficient … … 126 112 ! 127 113 l_zt_equal_zu = .FALSE. 128 IF( ABS(zu - zt) < 0.01 )l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision129 130 U_blk = MAX( 0.5 , U_zu ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s114 IF( ABS(zu - zt) < 0.01_wp ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 115 116 U_blk = MAX( 0.5_wp , U_zu ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 131 117 132 118 !! First guess of stability: 133 ztmp0 = t_zt*(1. + rctv0*q_zt) - sst*(1. + rctv0*ssq) ! air-sea difference of virtual pot. temp. at zt134 stab = 0.5 + sign(0.5,ztmp0) ! stab = 1 if dTv > 0 => STABLE, 0 if unstable119 ztmp0 = virt_temp(t_zt, q_zt) - virt_temp(sst, ssq) ! air-sea difference of virtual pot. temp. at zt 120 stab = 0.5_wp + sign(0.5_wp,ztmp0) ! stab = 1 if dTv > 0 => STABLE, 0 if unstable 135 121 136 122 !! Neutral coefficients at 10m: … … 139 125 ztmp0 (:,:) = cdn_wave(:,:) 140 126 ELSE 141 127 ztmp0 = cd_neutral_10m( U_blk ) 142 128 ENDIF 143 129 … … 146 132 !! Initializing transf. coeff. with their first guess neutral equivalents : 147 133 Cd = ztmp0 148 Ce = 1.e-3 *( 34.6* sqrt_Cd_n10 )149 Ch = 1.e-3 *sqrt_Cd_n10*(18.*stab + 32.7*(1.- stab))134 Ce = 1.e-3_wp*( 34.6_wp * sqrt_Cd_n10 ) 135 Ch = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab)) 150 136 stab = sqrt_Cd_n10 ! Temporaty array !!! stab == SQRT(Cd) 151 137 152 IF( ln_cdgw ) Cen = Ce ; Chn = Ch 138 IF( ln_cdgw ) THEN 139 Cen = Ce 140 Chn = Ch 141 ENDIF 153 142 154 143 !! Initializing values at z_u with z_t values: 155 144 t_zu = t_zt ; q_zu = q_zt 156 145 157 !! * Now starting iteration loop158 DO j_itt =1, nb_itt146 !! ITERATION BLOCK 147 DO j_itt = 1, nb_itt 159 148 ! 160 149 ztmp1 = t_zu - sst ! Updating air/sea differences … … 162 151 163 152 ! Updating turbulent scales : (L&Y 2004 eq. (7)) 164 ztmp1 = Ch/stab*ztmp1 ! theta* (stab == SQRT(Cd)) 165 ztmp2 = Ce/stab*ztmp2 ! q* (stab == SQRT(Cd)) 166 167 ztmp0 = 1. + rctv0*q_zu ! multiply this with t and you have the virtual temperature 153 ztmp0 = stab*U_blk ! u* (stab == SQRT(Cd)) 154 ztmp1 = Ch/stab*ztmp1 ! theta* (stab == SQRT(Cd)) 155 ztmp2 = Ce/stab*ztmp2 ! q* (stab == SQRT(Cd)) 168 156 169 157 ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 170 ztmp0 = (grav*vkarmn/(t_zu*ztmp0)*(ztmp1*ztmp0 + rctv0*t_zu*ztmp2)) / (Cd*U_blk*U_blk) 171 ! ( Cd*U_blk*U_blk is U*^2 at zu ) 172 158 ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 ) 159 173 160 !! Stability parameters : 174 zeta_u = zu*ztmp0 ; zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 161 zeta_u = zu*ztmp0 162 zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u ) 175 163 zpsi_h_u = psi_h( zeta_u ) 176 164 … … 178 166 IF( .NOT. l_zt_equal_zu ) THEN 179 167 !! Array 'stab' is free for the moment so using it to store 'zeta_t' 180 stab = zt*ztmp0 ; stab = SIGN( MIN(ABS(stab),10.0), stab ) ! Temporaty array stab == zeta_t !!! 168 stab = zt*ztmp0 169 stab = SIGN( MIN(ABS(stab),10._wp), stab ) ! Temporaty array stab == zeta_t !!! 181 170 stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab) ! stab just used as temp array again! 182 171 t_zu = t_zt - ztmp1/vkarmn*stab ! ztmp1 is still theta* L&Y 2004 eq.(9b) 183 172 q_zu = q_zt - ztmp2/vkarmn*stab ! ztmp2 is still q* L&Y 2004 eq.(9c) 184 q_zu = max(0., q_zu) 185 END IF 186 173 q_zu = max(0._wp, q_zu) 174 ENDIF 175 176 ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 177 ! In very rare low-wind conditions, the old way of estimating the 178 ! neutral wind speed at 10m leads to a negative value that causes the code 179 ! to crash. To prevent this a threshold of 0.25m/s is imposed. 187 180 ztmp2 = psi_m(zeta_u) 188 181 IF( ln_cdgw ) THEN ! surface wave case 189 182 stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 ) ! (stab == SQRT(Cd)) 190 183 Cd = stab * stab 191 ztmp0 = (LOG(zu/10. ) - zpsi_h_u) / vkarmn / sqrt_Cd_n10184 ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 192 185 ztmp2 = stab / sqrt_Cd_n10 ! (stab == SQRT(Cd)) 193 ztmp1 = 1. + Chn * ztmp0186 ztmp1 = 1._wp + Chn * ztmp0 194 187 Ch = Chn * ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 195 ztmp1 = 1. + Cen * ztmp0188 ztmp1 = 1._wp + Cen * ztmp0 196 189 Ce = Cen * ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 197 190 198 191 ELSE 199 200 201 202 203 ztmp0 = MAX( 0.25 , U_blk/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u))204 205 206 207 208 stab = 0.5 + sign(0.5,zeta_u)! update stability209 Cx_n10 = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.- stab)) ! L&Y 2004 eq. (6c-6d) (Cx_n10 == Ch_n10)210 211 212 213 ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2) ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u))214 215 216 217 ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10218 219 ztmp1 = 1.+ Cx_n10*ztmp0 ! (Cx_n10 == Ch_n10)220 221 222 Cx_n10 = 1.e-3 * (34.6* sqrt_Cd_n10) ! L&Y 2004 eq. (6b) ! Cx_n10 == Ce_n10223 224 ztmp1 = 1.+ Cx_n10*ztmp0225 226 227 ! 228 END DO 229 ! 192 ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 193 ! In very rare low-wind conditions, the old way of estimating the 194 ! neutral wind speed at 10m leads to a negative value that causes the code 195 ! to crash. To prevent this a threshold of 0.25m/s is imposed. 196 ztmp0 = MAX( 0.25_wp , U_blk/(1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u)) 197 ztmp0 = cd_neutral_10m(ztmp0) ! Cd_n10 198 Cdn(:,:) = ztmp0 199 sqrt_Cd_n10 = sqrt(ztmp0) 200 201 stab = 0.5_wp + sign(0.5_wp,zeta_u) ! update stability 202 Cx_n10 = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab)) ! L&Y 2004 eq. (6c-6d) (Cx_n10 == Ch_n10) 203 Chn(:,:) = Cx_n10 204 205 !! Update of transfer coefficients: 206 ztmp1 = 1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2) ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u)) 207 Cd = ztmp0 / ( ztmp1*ztmp1 ) 208 stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 209 210 ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 211 ztmp2 = stab / sqrt_Cd_n10 ! (stab == SQRT(Cd)) 212 ztmp1 = 1._wp + Cx_n10*ztmp0 ! (Cx_n10 == Ch_n10) 213 Ch = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 214 215 Cx_n10 = 1.e-3_wp * (34.6_wp * sqrt_Cd_n10) ! L&Y 2004 eq. (6b) ! Cx_n10 == Ce_n10 216 Cen(:,:) = Cx_n10 217 ztmp1 = 1._wp + Cx_n10*ztmp0 218 Ce = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 219 ENDIF 220 221 END DO !DO j_itt = 1, nb_itt 222 230 223 END SUBROUTINE turb_ncar 231 224 … … 238 231 !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b) 239 232 !! 240 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)233 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 241 234 !!---------------------------------------------------------------------------------- 242 235 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10 ! scalar wind speed at 10m (m/s) … … 255 248 ! 256 249 ! When wind speed > 33 m/s => Cyclone conditions => special treatment 257 zgt33 = 0.5 + SIGN( 0.5, (zw - 33.) ) ! If pw10 < 33. => 0, else => 1258 ! 259 cd_neutral_10m(ji,jj) = 1.e-3 * ( &260 & (1. - zgt33)*( 2.7/zw + 0.142 + zw/13.09 - 3.14807E-10*zw6) & ! wind < 33 m/s261 & + zgt33 * 2.34 )! wind >= 33 m/s262 ! 263 cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6 )250 zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) ) ! If pw10 < 33. => 0, else => 1 251 ! 252 cd_neutral_10m(ji,jj) = 1.e-3_wp * ( & 253 & (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind < 33 m/s 254 & + zgt33 * 2.34_wp ) ! wind >= 33 m/s 255 ! 256 cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp) 264 257 ! 265 258 END DO … … 273 266 !! Universal profile stability function for momentum 274 267 !! !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 275 !! 276 !! pzet 0 : stability paramenter, z/L where z is altitude measurement268 !! 269 !! pzeta : stability paramenter, z/L where z is altitude measurement 277 270 !! and L is M-O length 278 271 !! 279 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)280 !!---------------------------------------------------------------------------------- 281 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: pzeta282 REAL(wp), DIMENSION(jpi,jpj) :: psi_m283 ! 284 INTEGER :: ji, jj 272 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 273 !!---------------------------------------------------------------------------------- 274 REAL(wp), DIMENSION(jpi,jpj) :: psi_m 275 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 276 ! 277 INTEGER :: ji, jj ! dummy loop indices 285 278 REAL(wp) :: zx2, zx, zstab ! local scalars 286 279 !!---------------------------------------------------------------------------------- 287 !288 280 DO jj = 1, jpj 289 281 DO ji = 1, jpi 290 zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) )291 zx2 = MAX ( zx2 , 1.)282 zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 283 zx2 = MAX( zx2 , 1._wp ) 292 284 zx = SQRT( zx2 ) 293 zstab = 0.5 + SIGN( 0.5, pzeta(ji,jj) )294 ! 295 psi_m(ji,jj) = zstab * (-5. *pzeta(ji,jj)) & ! Stable296 & + (1. - zstab) * (2.*LOG((1. + zx)*0.5) & ! Unstable297 & + LOG((1. + zx2)*0.5) - 2.*ATAN(zx) + rpi*0.5) ! "285 zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 286 ! 287 psi_m(ji,jj) = zstab * (-5._wp*pzeta(ji,jj)) & ! Stable 288 & + (1._wp - zstab) * (2._wp*LOG((1._wp + zx)*0.5_wp) & ! Unstable 289 & + LOG((1._wp + zx2)*0.5_wp) - 2._wp*ATAN(zx) + rpi*0.5_wp) ! " 298 290 ! 299 291 END DO 300 292 END DO 301 !302 293 END FUNCTION psi_m 303 294 … … 308 299 !! !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 309 300 !! 310 !! pzet 0 : stability paramenter, z/L where z is altitude measurement301 !! pzeta : stability paramenter, z/L where z is altitude measurement 311 302 !! and L is M-O length 312 303 !! 313 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 314 !!---------------------------------------------------------------------------------- 304 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 305 !!---------------------------------------------------------------------------------- 306 REAL(wp), DIMENSION(jpi,jpj) :: psi_h 315 307 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 316 REAL(wp), DIMENSION(jpi,jpj) :: psi_h 317 ! 318 INTEGER :: ji, jj ! dummy loop indices 308 ! 309 INTEGER :: ji, jj ! dummy loop indices 319 310 REAL(wp) :: zx2, zstab ! local scalars 320 311 !!---------------------------------------------------------------------------------- … … 322 313 DO jj = 1, jpj 323 314 DO ji = 1, jpi 324 zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) )325 zx2 = MAX ( zx2 , 1.)326 zstab = 0.5 + SIGN( 0.5, pzeta(ji,jj) )327 ! 328 psi_h(ji,jj) = zstab * (-5. *pzeta(ji,jj)) & ! Stable329 & + (1. - zstab) * (2.*LOG( (1. + zx2)*0.5)) ! Unstable315 zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 316 zx2 = MAX( zx2 , 1._wp ) 317 zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 318 ! 319 psi_h(ji,jj) = zstab * (-5._wp*pzeta(ji,jj)) & ! Stable 320 & + (1._wp - zstab) * (2._wp*LOG( (1._wp + zx2)*0.5_wp )) ! Unstable 330 321 ! 331 322 END DO 332 323 END DO 333 !334 324 END FUNCTION psi_h 335 325
Note: See TracChangeset
for help on using the changeset viewer.