- Timestamp:
- Nov 21, 2014 5:14:03 PM (10 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/prandtl_fluxes.f90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/prandtl_fluxes.f90
r1362 r1494 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Bugfixes: qs is now calculated before calculation of Rif. Ccalculation of 23 ! buoyancy flux in Rif corrected (added missing humidity term), allow use of 24 ! topography for coupled runs (not tested) 23 25 ! 24 26 ! Former revisions: … … 129 131 130 132 IF ( large_scale_forcing .AND. lsf_surf ) THEN 131 pt( 0,:,:) = pt_surface133 pt(nzb_s_inner(j,i),:,:) = pt_surface 132 134 ENDIF 133 135 … … 153 155 b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p ) 154 156 155 ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( &156 LOG( z_p / z0h(j,i) ) - &157 ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( & 158 LOG( z_p / z0h(j,i) ) - & 157 159 2.0_wp * LOG( ( 1.0_wp + a ) / ( 1.0_wp + b ) ) ) 158 160 ENDIF … … 160 162 ENDDO 161 163 ENDDO 164 ENDIF 165 166 ! 167 !-- If required compute q* 168 IF ( humidity .OR. passive_scalar ) THEN 169 IF ( constant_waterflux ) THEN 170 ! 171 !-- For a given water flux in the Prandtl layer: 172 !$OMP PARALLEL DO 173 !$acc kernels loop 174 DO i = nxlg, nxrg 175 DO j = nysg, nyng 176 qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp ) 177 ENDDO 178 ENDDO 179 180 ELSE 181 coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. run_coupled ) 182 183 IF ( large_scale_forcing .AND. lsf_surf ) THEN 184 q(nzb_s_inner(j,i),:,:) = q_surface 185 ENDIF 186 187 !$OMP PARALLEL DO PRIVATE( a, b, k, z_p ) 188 !$acc kernels loop independent 189 DO i = nxlg, nxrg 190 !$acc loop independent 191 DO j = nysg, nyng 192 193 k = nzb_s_inner(j,i) 194 z_p = zu(k+1) - zw(k) 195 196 ! 197 !-- Assume saturation for atmosphere coupled to ocean (but not 198 !-- in case of precursor runs) 199 IF ( coupled_run ) THEN 200 e_q = 6.1_wp * & 201 EXP( 0.07_wp * ( MIN(pt(k,j,i),pt(k+1,j,i)) - 273.15_wp ) ) 202 q(k,j,i) = 0.622_wp * e_q / ( surface_pressure - e_q ) 203 ENDIF 204 IF ( rif(j,i) >= 0.0_wp ) THEN 205 ! 206 !-- Stable stratification 207 qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( & 208 LOG( z_p / z0h(j,i) ) + & 209 5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p & 210 ) 211 ELSE 212 ! 213 !-- Unstable stratification 214 a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) ) 215 b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p ) 216 217 qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( & 218 LOG( z_p / z0h(j,i) ) - & 219 2.0_wp * LOG( (1.0_wp + a ) / ( 1.0_wp + b ) ) ) 220 ENDIF 221 222 ENDDO 223 ENDDO 224 ENDIF 162 225 ENDIF 163 226 … … 190 253 k = nzb_s_inner(j,i) 191 254 z_p = zu(k+1) - zw(k) 192 rif(j,i) = z_p * kappa * g * & 193 ( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i) ) / & 255 rif(j,i) = z_p * kappa * g * & 256 ( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i) + 0.61_wp & 257 * q(k+1,j,i) * ts(j,i)) / & 194 258 ( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30_wp ) ) 195 259 ! … … 218 282 !-- Compute the absolute value of the horizontal velocity 219 283 !-- (relative to the surface) 220 uv_total = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) &221 - u(k,j,i) - u(k,j,i+1) ) )**2 + &222 ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) &284 uv_total = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) & 285 - u(k,j,i) - u(k,j,i+1) ) )**2 + & 286 ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) & 223 287 - v(k,j,i) - v(k,j+1,i) ) )**2 ) 224 288 … … 227 291 ! 228 292 !-- Stable stratification 229 us(j,i) = kappa * uv_total / ( &230 LOG( z_p / z0(j,i) ) + &231 5.0_wp * rif(j,i) * ( z_p - z0(j,i) ) / z_p &293 us(j,i) = kappa * uv_total / ( & 294 LOG( z_p / z0(j,i) ) + & 295 5.0_wp * rif(j,i) * ( z_p - z0(j,i) ) / z_p & 232 296 ) 233 297 ELSE … … 237 301 b = SQRT( SQRT( 1.0_wp - 16.0_wp * rif(j,i) / z_p * z0(j,i) ) ) 238 302 239 us(j,i) = kappa * uv_total / ( &240 LOG( z_p / z0(j,i) ) - &241 LOG( ( 1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / ( 242 ( 1.0_wp + b )**2 * ( 1.0_wp + b**2 ) ) ) + 243 2.0_wp * ( ATAN( a ) - ATAN( b ) ) 303 us(j,i) = kappa * uv_total / ( & 304 LOG( z_p / z0(j,i) ) - & 305 LOG( ( 1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / ( & 306 ( 1.0_wp + b )**2 * ( 1.0_wp + b**2 ) ) ) + & 307 2.0_wp * ( ATAN( a ) - ATAN( b ) ) & 244 308 ) 245 309 ENDIF … … 271 335 ! 272 336 !-- Stable stratification 273 usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) )/ ( &274 LOG( z_p / z0(j,i) ) + &275 5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p &337 usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) )/ ( & 338 LOG( z_p / z0(j,i) ) + & 339 5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p & 276 340 ) 277 341 ELSE … … 281 345 b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) ) 282 346 283 usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) ) / ( &284 LOG( z_p / z0(j,i) ) - &285 LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / ( 286 (1.0_wp + b )**2 * ( 1.0_wp + b**2 ) ) ) + 287 2.0_wp * ( ATAN( a ) - ATAN( b ) ) 347 usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) ) / ( & 348 LOG( z_p / z0(j,i) ) - & 349 LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / ( & 350 (1.0_wp + b )**2 * ( 1.0_wp + b**2 ) ) ) + & 351 2.0_wp * ( ATAN( a ) - ATAN( b ) ) & 288 352 ) 289 353 ENDIF … … 309 373 ! 310 374 !-- Stable stratification 311 vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) / ( &312 LOG( z_p / z0(j,i) ) + &313 5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p &375 vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) / ( & 376 LOG( z_p / z0(j,i) ) + & 377 5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p & 314 378 ) 315 379 ELSE … … 319 383 b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) ) 320 384 321 vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) / ( &322 LOG( z_p / z0(j,i) ) - &323 LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / ( 324 (1.0_wp + b )**2 * ( 1.0_wp + b**2 ) ) ) + 325 2.0_wp * ( ATAN( a ) - ATAN( b ) ) 385 vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) / ( & 386 LOG( z_p / z0(j,i) ) - & 387 LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / ( & 388 (1.0_wp + b )**2 * ( 1.0_wp + b**2 ) ) ) + & 389 2.0_wp * ( ATAN( a ) - ATAN( b ) ) & 326 390 ) 327 391 ENDIF … … 331 395 332 396 ! 333 !-- If required compute q* 334 IF ( humidity .OR. passive_scalar ) THEN 335 IF ( constant_waterflux ) THEN 336 ! 337 !-- For a given water flux in the Prandtl layer: 338 !$OMP PARALLEL DO 339 !$acc kernels loop 340 DO i = nxlg, nxrg 341 DO j = nysg, nyng 342 qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp ) 343 ENDDO 344 ENDDO 345 346 ELSE 347 coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. run_coupled ) 348 349 IF ( large_scale_forcing .AND. lsf_surf ) THEN 350 q(0,:,:) = q_surface 351 ENDIF 352 353 !$OMP PARALLEL DO PRIVATE( a, b, k, z_p ) 354 !$acc kernels loop independent 355 DO i = nxlg, nxrg 356 !$acc loop independent 357 DO j = nysg, nyng 358 359 k = nzb_s_inner(j,i) 360 z_p = zu(k+1) - zw(k) 361 362 ! 363 !-- Assume saturation for atmosphere coupled to ocean (but not 364 !-- in case of precursor runs) 365 IF ( coupled_run ) THEN 366 e_q = 6.1_wp * & 367 EXP( 0.07_wp * ( MIN(pt(0,j,i),pt(1,j,i)) - 273.15_wp ) ) 368 q(k,j,i) = 0.622_wp * e_q / ( surface_pressure - e_q ) 369 ENDIF 370 IF ( rif(j,i) >= 0.0_wp ) THEN 371 ! 372 !-- Stable stratification 373 qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( & 374 LOG( z_p / z0h(j,i) ) + & 375 5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p & 376 ) 377 ELSE 378 ! 379 !-- Unstable stratification 380 a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) ) 381 b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p ) 397 !-- If required compute qr* and nr* 398 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation ) THEN 399 400 !$OMP PARALLEL DO PRIVATE( a, b, k, z_p ) 401 !$acc kernels loop independent 402 DO i = nxlg, nxrg 403 !$acc loop independent 404 DO j = nysg, nyng 405 406 k = nzb_s_inner(j,i) 407 z_p = zu(k+1) - zw(k) 408 409 IF ( rif(j,i) >= 0.0 ) THEN 410 ! 411 !-- Stable stratification 412 qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) / ( & 413 LOG( z_p / z0h(j,i) ) + & 414 5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p ) 415 nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) / ( & 416 LOG( z_p / z0h(j,i) ) + & 417 5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p ) 418 419 ELSE 420 ! 421 !-- Unstable stratification 422 a = SQRT( 1.0 - 16.0 * rif(j,i) ) 423 b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p ) 382 424 383 qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( & 384 LOG( z_p / z0h(j,i) ) - & 385 2.0_wp * LOG( (1.0_wp + a ) / ( 1.0_wp + b ) ) ) 386 ENDIF 387 388 ENDDO 389 ENDDO 390 ENDIF 391 392 IF ( cloud_physics .AND. icloud_scheme == 0 & 393 .AND. precipitation ) THEN 394 395 !$OMP PARALLEL DO PRIVATE( a, b, k, z_p ) 396 !$acc kernels loop independent 397 DO i = nxlg, nxrg 398 !$acc loop independent 399 DO j = nysg, nyng 400 401 k = nzb_s_inner(j,i) 402 z_p = zu(k+1) - zw(k) 403 404 IF ( rif(j,i) >= 0.0 ) THEN 405 ! 406 !-- Stable stratification 407 qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) / ( & 408 LOG( z_p / z0h(j,i) ) + & 409 5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p ) 410 nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) / ( & 411 LOG( z_p / z0h(j,i) ) + & 412 5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p ) 413 414 ELSE 415 ! 416 !-- Unstable stratification 417 a = SQRT( 1.0 - 16.0 * rif(j,i) ) 418 b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p ) 419 420 qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) / ( & 421 LOG( z_p / z0h(j,i) ) - & 422 2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) ) 423 nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) / ( & 424 LOG( z_p / z0h(j,i) ) - & 425 2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) ) 426 427 ENDIF 428 429 ENDDO 430 ENDDO 431 432 ENDIF 425 qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) / ( & 426 LOG( z_p / z0h(j,i) ) - & 427 2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) ) 428 nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) / ( & 429 LOG( z_p / z0h(j,i) ) - & 430 2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) ) 431 432 ENDIF 433 434 ENDDO 435 ENDDO 433 436 434 437 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
