Changeset 2232 for palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
 Timestamp:
 May 30, 2017 5:47:52 PM (4 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
r2119 r2232 20 20 ! Current revisions: 21 21 !  22 ! 22 ! Adjustments to new surface concept 23 ! OpenMP bugfix 23 24 ! 24 25 ! Former revisions: … … 171 172 172 173 USE arrays_3d, & 173 ONLY: e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, & 174 s, shf, ss, ssws, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, & 175 z0h, z0q, drho_air_zw, rho_air_zw 174 ONLY: e, kh, nr, pt, q, ql, qr, s, u, v, vpt, w, zu, zw, drho_air_zw, & 175 rho_air_zw 176 176 177 177 USE cloud_parameters, & … … 187 187 constant_waterflux, coupling_mode, g, humidity, ibc_e_b, & 188 188 ibc_pt_b, initializing_actions, kappa, & 189 intermediate_timestep_count, 190 intermediate_timestep_count_max, large_scale_forcing, lsf_surf,&189 intermediate_timestep_count, intermediate_timestep_count_max, & 190 land_surface, large_scale_forcing, lsf_surf, & 191 191 message_string, microphysics_seifert, most_method, neutral, & 192 192 passive_scalar, pt_surface, q_surface, run_coupled, & 193 193 surface_pressure, simulated_time, terminate_run, & 194 urban_surface, & 195 zeta_max, zeta_min 194 urban_surface, zeta_max, zeta_min 195 196 USE grid_variables, & 197 ONLY: dx, dy 196 198 197 199 USE indices, & 198 ONLY: nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_s_inner, & 199 nzb_u_inner, nzb_v_inner 200 ONLY: nxl, nxr, nys, nyn, nzb 200 201 201 202 USE kinds … … 204 205 205 206 USE land_surface_model_mod, & 206 ONLY: land_surface, skip_time_do_lsm 207 207 ONLY: aero_resist_kray, skip_time_do_lsm 208 209 USE surface_mod, & 210 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_type, & 211 surf_usm_h, surf_usm_v 208 212 209 213 … … 213 217 INTEGER(iwp) :: j !< loop index y direction 214 218 INTEGER(iwp) :: k !< loop index z direction 215 INTEGER(iwp) :: l_bnd = 7500 !< Lookup table index of the last time step 216 217 INTEGER(iwp), PARAMETER :: num_steps = 15000 !< number of steps in the lookup table 218 219 LOGICAL :: coupled_run !< Flag for coupled atmosphereocean runs 220 221 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt1, & !< Potential temperature at first grid level (required for cloud_physics = .T.) 222 qv1, & !< Specific humidity at first grid level (required for cloud_physics = .T.) 223 uv_total !< Total velocity at first grid level 219 INTEGER(iwp) :: l !< loop index for surf type 220 INTEGER(iwp) :: li_bnd = 7500 !< Lookup table index of the last time step 221 222 INTEGER(iwp), PARAMETER :: num_steps = 15000 !< number of steps in the lookup table 223 224 LOGICAL :: coupled_run !< Flag for coupled atmosphereocean runs 225 LOGICAL :: downward = .FALSE.!< Flag indicating downwardfacing horizontal surface 226 LOGICAL :: mom_uv = .FALSE. !< Flag indicating calculation of usvs and vsus at vertical surfaces 227 LOGICAL :: mom_w = .FALSE. !< Flag indicating calculation of wsus and wsvs at vertical surfaces 228 LOGICAL :: mom_tke = .FALSE. !< Flag indicating calculation of momentum fluxes at vertical surfaces used for TKE production 229 LOGICAL :: surf_vertical !< Flag indicating vertical surfaces 224 230 225 231 REAL(wp), DIMENSION(0:num_steps1) :: rib_tab, & !< Lookup table bulk Richardson number … … 232 238 z_mo !< Height of the constant flux layer where MOST is assumed 233 239 240 TYPE(surf_type), POINTER :: surf !< surftype array, used to generalize subroutines 241 234 242 235 243 SAVE … … 237 245 PRIVATE 238 246 239 PUBLIC init_surface_layer_fluxes, pt1, qv1, surface_layer_fluxes, uv_total247 PUBLIC init_surface_layer_fluxes, surface_layer_fluxes 240 248 241 249 INTERFACE init_surface_layer_fluxes … … 260 268 IMPLICIT NONE 261 269 270 surf_vertical = .FALSE. 271 downward = .FALSE. 262 272 ! 263 273 ! In case cloud physics is used, it is required to derive potential … … 265 275 ! and q 266 276 IF ( cloud_physics ) THEN 267 CALL calc_pt_q 277 ! 278 ! First call for horizontal defaulttype surfaces (l=0  upward facing, 279 ! l=1  downward facing) 280 DO l = 0, 1 281 IF ( surf_def_h(l)%ns >= 1 ) THEN 282 surf => surf_def_h(l) 283 CALL calc_pt_q 284 ENDIF 285 ENDDO 286 ! 287 ! Call for naturaltype surfaces 288 IF ( surf_lsm_h%ns >= 1 ) THEN 289 surf => surf_lsm_h 290 CALL calc_pt_q 291 ENDIF 292 ! 293 ! Call for urbantype surfaces 294 IF ( surf_usm_h%ns >= 1 ) THEN 295 surf => surf_usm_h 296 CALL calc_pt_q 297 ENDIF 268 298 ENDIF 269 299 … … 277 307 ! 278 308 ! Depending on setting of most_method use the "old" routine 309 ! Note, each routine is called for different surface types. 310 ! First call for defaulttype horizontal surfaces, for natural and 311 ! urbantype surfaces. Note, at this place only upwardfacing horizontal 312 ! surfaces are treted. 279 313 IF ( most_method == 'circular' ) THEN 280 281 CALL calc_scaling_parameters 282 283 CALL calc_uv_total 284 285 IF ( .NOT. neutral ) THEN 286 CALL calc_ol 287 ENDIF 288 289 CALL calc_us 290 291 CALL calc_surface_fluxes 292 314 ! 315 ! Defaulttype upwardfacing horizontal surfaces 316 IF ( surf_def_h(0)%ns >= 1 ) THEN 317 surf => surf_def_h(0) 318 CALL calc_scaling_parameters 319 CALL calc_uvw_abs 320 IF ( .NOT. neutral ) CALL calc_ol 321 CALL calc_us 322 CALL calc_surface_fluxes 323 ENDIF 324 ! 325 ! Naturaltype horizontal surfaces 326 IF ( surf_lsm_h%ns >= 1 ) THEN 327 surf => surf_lsm_h 328 CALL calc_scaling_parameters 329 CALL calc_uvw_abs 330 IF ( .NOT. neutral ) CALL calc_ol 331 CALL calc_us 332 CALL calc_surface_fluxes 333 ENDIF 334 ! 335 ! Urbantype horizontal surfaces 336 IF ( surf_usm_h%ns >= 1 ) THEN 337 surf => surf_usm_h 338 CALL calc_scaling_parameters 339 CALL calc_uvw_abs 340 IF ( .NOT. neutral ) CALL calc_ol 341 CALL calc_us 342 CALL calc_surface_fluxes 343 ENDIF 293 344 ! 294 345 ! Use either Newton iteration or a lookup table for the bulk Richardson 295 346 ! number to calculate the Obukhov length 296 347 ELSEIF ( most_method == 'newton' .OR. most_method == 'lookup' ) THEN 297 298 CALL calc_uv_total 299 300 IF ( .NOT. neutral ) THEN 301 CALL calc_ol 302 ENDIF 303 348 ! 349 ! Defaulttype upwardfacing horizontal surfaces 350 IF ( surf_def_h(0)%ns >= 1 ) THEN 351 surf => surf_def_h(0) 352 CALL calc_uvw_abs 353 IF ( .NOT. neutral ) CALL calc_ol 354 CALL calc_us 355 CALL calc_scaling_parameters 356 CALL calc_surface_fluxes 357 ENDIF 358 ! 359 ! Naturaltype horizontal surfaces 360 IF ( surf_lsm_h%ns >= 1 ) THEN 361 surf => surf_lsm_h 362 CALL calc_uvw_abs 363 IF ( .NOT. neutral ) CALL calc_ol 364 CALL calc_us 365 CALL calc_scaling_parameters 366 CALL calc_surface_fluxes 367 ENDIF 368 ! 369 ! Urbantype horizontal surfaces 370 IF ( surf_usm_h%ns >= 1 ) THEN 371 surf => surf_usm_h 372 CALL calc_uvw_abs 373 IF ( .NOT. neutral ) CALL calc_ol 374 CALL calc_us 375 CALL calc_scaling_parameters 376 CALL calc_surface_fluxes 377 ENDIF 378 379 ENDIF 380 ! 381 ! Treat downwardfacing horizontal surfaces. Note, so far, these are 382 ! always default type. Stratification is not considered 383 ! in this case, hence, no further distinction between different 384 ! most_method is required. 385 IF ( surf_def_h(1)%ns >= 1 ) THEN 386 downward = .TRUE. 387 surf => surf_def_h(1) 388 CALL calc_uvw_abs 304 389 CALL calc_us 305 306 CALL calc_scaling_parameters307 308 390 CALL calc_surface_fluxes 309 391 downward = .FALSE. 310 392 ENDIF 393 ! 394 ! Calculate surfaces fluxes at vertical surfaces for momentum 395 ! and subgridscale TKE. 396 ! No stability is considered. Therefore, scaling parameters and Obukhov 397 ! length do not need to be calculated and no distinction in 'circular', 398 ! 'Newton' or 'lookup' is necessary so far. 399 ! Note, this will change if stability is once considered. 400 surf_vertical = .TRUE. 401 ! 402 ! Calculate horizontal momentum fluxes at north and southfacing 403 ! surfaces(usvs). 404 ! For defaulttype surfaces 405 mom_uv = .TRUE. 406 DO l = 0, 1 407 IF ( surf_def_v(l)%ns >= 1 ) THEN 408 surf => surf_def_v(l) 409 ! 410 ! Compute surfaceparallel velocity 411 CALL calc_uvw_abs_v_ugrid 412 ! 413 ! Compute respective friction velocity on staggered grid 414 CALL calc_us 415 ! 416 ! Compute respective surface fluxes for momentum and TKE 417 CALL calc_surface_fluxes 418 ENDIF 419 ENDDO 420 ! 421 ! For naturaltype surfaces. Please note, even though stability is not 422 ! considered for the calculation of momentum fluxes at vertical surfaces, 423 ! scaling parameters and Obukov length are calculated nevertheless in this 424 ! case. This is due to the requirement of ts in parameterization of heat 425 ! flux in landsurface model in case of aero_resist_kray is not true. 426 IF ( .NOT. aero_resist_kray ) THEN 427 IF ( most_method == 'circular' ) THEN 428 DO l = 0, 1 429 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 430 surf => surf_lsm_v(l) 431 ! 432 ! Compute scaling parameters 433 CALL calc_scaling_parameters 434 ! 435 ! Compute surfaceparallel velocity 436 CALL calc_uvw_abs_v_ugrid 437 ! 438 ! Compute Obukhov length 439 IF ( .NOT. neutral ) CALL calc_ol 440 ! 441 ! Compute respective friction velocity on staggered grid 442 CALL calc_us 443 ! 444 ! Compute respective surface fluxes for momentum and TKE 445 CALL calc_surface_fluxes 446 ENDIF 447 ENDDO 448 ELSE 449 DO l = 0, 1 450 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 451 surf => surf_lsm_v(l) 452 ! 453 ! Compute surfaceparallel velocity 454 CALL calc_uvw_abs_v_ugrid 455 ! 456 ! Compute Obukhov length 457 IF ( .NOT. neutral ) CALL calc_ol 458 ! 459 ! Compute respective friction velocity on staggered grid 460 CALL calc_us 461 ! 462 ! Compute scaling parameters 463 CALL calc_scaling_parameters 464 ! 465 ! Compute respective surface fluxes for momentum and TKE 466 CALL calc_surface_fluxes 467 ENDIF 468 ENDDO 469 ENDIF 470 ! 471 ! No ts is required, so scaling parameters and Obukhov length do not need 472 ! to be computed. 473 ELSE 474 DO l = 0, 1 475 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 476 surf => surf_lsm_v(l) 477 ! 478 ! Compute surfaceparallel velocity 479 CALL calc_uvw_abs_v_ugrid 480 ! 481 ! Compute respective friction velocity on staggered grid 482 CALL calc_us 483 ! 484 ! Compute respective surface fluxes for momentum and TKE 485 CALL calc_surface_fluxes 486 ENDIF 487 ENDDO 488 ENDIF 489 ! 490 ! For urbantype surfaces 491 DO l = 0, 1 492 IF ( surf_usm_v(l)%ns >= 1 ) THEN 493 surf => surf_usm_v(l) 494 ! 495 ! Compute surfaceparallel velocity 496 CALL calc_uvw_abs_v_ugrid 497 ! 498 ! Compute respective friction velocity on staggered grid 499 CALL calc_us 500 ! 501 ! Compute respective surface fluxes for momentum and TKE 502 CALL calc_surface_fluxes 503 ENDIF 504 ENDDO 505 ! 506 ! Calculate horizontal momentum fluxes at east and westfacing 507 ! surfaces (vsus). 508 ! For defaulttype surfaces 509 DO l = 2, 3 510 IF ( surf_def_v(l)%ns >= 1 ) THEN 511 surf => surf_def_v(l) 512 ! 513 ! Compute surfaceparallel velocity 514 CALL calc_uvw_abs_v_vgrid 515 ! 516 ! Compute respective friction velocity on staggered grid 517 CALL calc_us 518 ! 519 ! Compute respective surface fluxes for momentum and TKE 520 CALL calc_surface_fluxes 521 ENDIF 522 ENDDO 523 ! 524 ! For naturaltype surfaces. Please note, even though stability is not 525 ! considered for the calculation of momentum fluxes at vertical surfaces, 526 ! scaling parameters and Obukov length are calculated nevertheless in this 527 ! case. This is due to the requirement of ts in parameterization of heat 528 ! flux in landsurface model in case of aero_resist_kray is not true. 529 IF ( .NOT. aero_resist_kray ) THEN 530 IF ( most_method == 'circular' ) THEN 531 DO l = 2, 3 532 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 533 surf => surf_lsm_v(l) 534 ! 535 ! Compute scaling parameters 536 CALL calc_scaling_parameters 537 ! 538 ! Compute surfaceparallel velocity 539 CALL calc_uvw_abs_v_vgrid 540 ! 541 ! Compute Obukhov length 542 IF ( .NOT. neutral ) CALL calc_ol 543 ! 544 ! Compute respective friction velocity on staggered grid 545 CALL calc_us 546 ! 547 ! Compute respective surface fluxes for momentum and TKE 548 CALL calc_surface_fluxes 549 ENDIF 550 ENDDO 551 ELSE 552 DO l = 2, 3 553 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 554 surf => surf_lsm_v(l) 555 ! 556 ! Compute surfaceparallel velocity 557 CALL calc_uvw_abs_v_vgrid 558 ! 559 ! Compute Obukhov length 560 IF ( .NOT. neutral ) CALL calc_ol 561 ! 562 ! Compute respective friction velocity on staggered grid 563 CALL calc_us 564 ! 565 ! Compute scaling parameters 566 CALL calc_scaling_parameters 567 ! 568 ! Compute respective surface fluxes for momentum and TKE 569 CALL calc_surface_fluxes 570 ENDIF 571 ENDDO 572 ENDIF 573 ELSE 574 DO l = 2, 3 575 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 576 surf => surf_lsm_v(l) 577 ! 578 ! Compute surfaceparallel velocity 579 CALL calc_uvw_abs_v_vgrid 580 ! 581 ! Compute respective friction velocity on staggered grid 582 CALL calc_us 583 ! 584 ! Compute respective surface fluxes for momentum and TKE 585 CALL calc_surface_fluxes 586 ENDIF 587 ENDDO 588 ENDIF 589 ! 590 ! For urbantype surfaces 591 DO l = 2, 3 592 IF ( surf_usm_v(l)%ns >= 1 ) THEN 593 surf => surf_usm_v(l) 594 ! 595 ! Compute surfaceparallel velocity 596 CALL calc_uvw_abs_v_vgrid 597 ! 598 ! Compute respective friction velocity on staggered grid 599 CALL calc_us 600 ! 601 ! Compute respective surface fluxes for momentum and TKE 602 CALL calc_surface_fluxes 603 ENDIF 604 ENDDO 605 mom_uv = .FALSE. 606 ! 607 ! Calculate horizontal momentum fluxes of w (wsus and wsvs) at vertial 608 ! surfaces. 609 mom_w = .TRUE. 610 ! 611 ! Defaulttype surfaces 612 DO l = 0, 3 613 IF ( surf_def_v(l)%ns >= 1 ) THEN 614 surf => surf_def_v(l) 615 CALL calc_uvw_abs_v_wgrid 616 CALL calc_us 617 CALL calc_surface_fluxes 618 ENDIF 619 ENDDO 620 ! 621 ! Naturaltype surfaces 622 DO l = 0, 3 623 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 624 surf => surf_lsm_v(l) 625 CALL calc_uvw_abs_v_wgrid 626 CALL calc_us 627 CALL calc_surface_fluxes 628 ENDIF 629 ENDDO 630 ! 631 ! Urbantype surfaces 632 DO l = 0, 3 633 IF ( surf_usm_v(l)%ns >= 1 ) THEN 634 surf => surf_usm_v(l) 635 CALL calc_uvw_abs_v_wgrid 636 CALL calc_us 637 CALL calc_surface_fluxes 638 ENDIF 639 ENDDO 640 mom_w = .FALSE. 641 ! 642 ! Calculate momentum fluxes usvs, vsus, wsus and wsvs at vertical 643 ! surfaces for TKE production. Note, here, momentum fluxes are defined 644 ! at grid center and are not staggered as before. 645 mom_tke = .TRUE. 646 ! 647 ! Defaulttype surfaces 648 DO l = 0, 3 649 IF ( surf_def_v(l)%ns >= 1 ) THEN 650 surf => surf_def_v(l) 651 CALL calc_uvw_abs_v_sgrid 652 CALL calc_us 653 CALL calc_surface_fluxes 654 ENDIF 655 ENDDO 656 ! 657 ! Naturaltype surfaces 658 DO l = 0, 3 659 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 660 surf => surf_lsm_v(l) 661 CALL calc_uvw_abs_v_sgrid 662 CALL calc_us 663 CALL calc_surface_fluxes 664 ENDIF 665 ENDDO 666 ! 667 ! Urbantype surfaces 668 DO l = 0, 3 669 IF ( surf_usm_v(l)%ns >= 1 ) THEN 670 surf => surf_usm_v(l) 671 CALL calc_uvw_abs_v_sgrid 672 CALL calc_us 673 CALL calc_surface_fluxes 674 ENDIF 675 ENDDO 676 mom_tke = .FALSE. 677 311 678 312 679 END SUBROUTINE surface_layer_fluxes … … 324 691 IMPLICIT NONE 325 692 326 INTEGER(iwp) :: l ,& !< Index for loop to create lookup table693 INTEGER(iwp) :: li, & !< Index for loop to create lookup table 327 694 num_steps_n !< Number of nonstretched zeta steps 328 695 … … 339 706 z0h_min = 0.0_wp, & !< Minimum value of z0h to create table 340 707 z0_min = 0.0_wp !< Minimum value of z0 to create table 341 ! 342 ! When cloud physics is used, arrays for storing potential temperature and 343 ! specific humidity at first grid level are required 344 IF ( cloud_physics ) THEN 345 ALLOCATE ( pt1(nysg:nyng,nxlg:nxrg) ) 346 ALLOCATE ( qv1(nysg:nyng,nxlg:nxrg) ) 347 ENDIF 348 349 ! 350 ! Allocate field for storing the horizontal velocity 351 ALLOCATE ( uv_total(nysg:nyng,nxlg:nxrg) ) 708 709 352 710 353 711 … … 355 713 ! In case of runs with neutral statification, set Obukhov length to a 356 714 ! large value 357 IF ( neutral ) ol = 1.0E10_wp 715 IF ( neutral ) THEN 716 IF ( surf_def_h(0)%ns >= 1 ) surf_def_h(0)%ol = 1.0E10_wp 717 IF ( surf_lsm_h%ns >= 1 ) surf_lsm_h%ol = 1.0E10_wp 718 IF ( surf_usm_h%ns >= 1 ) surf_usm_h%ol = 1.0E10_wp 719 ENDIF 358 720 359 721 IF ( most_method == 'lookup' ) THEN … … 361 723 ! 362 724 ! Check for roughness heterogeneity. In that case terminate run and 363 ! inform user 364 IF ( MINVAL( z0h ) /= MAXVAL( z0h ) .OR. & 365 MINVAL( z0 ) /= MAXVAL( z0 ) ) THEN 366 terminate_run_l = .TRUE. 367 ENDIF 725 ! inform user. Check for both, natural and nonnatural walls. 726 IF ( surf_def_h(0)%ns >= 1 ) THEN 727 IF ( MINVAL( surf_def_h(0)%z0h ) /= MAXVAL( surf_def_h(0)%z0h ) .OR. & 728 MINVAL( surf_def_h(0)%z0 ) /= MAXVAL( surf_def_h(0)%z0 ) ) THEN 729 terminate_run_l = .TRUE. 730 ENDIF 731 ENDIF 732 IF ( surf_lsm_h%ns >= 1 ) THEN 733 IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h ) .OR. & 734 MINVAL( surf_lsm_h%z0 ) /= MAXVAL( surf_lsm_h%z0 ) ) THEN 735 terminate_run_l = .TRUE. 736 ENDIF 737 ENDIF 738 IF ( surf_usm_h%ns >= 1 ) THEN 739 IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_usm_h%z0h ) .OR. & 740 MINVAL( surf_usm_h%z0 ) /= MAXVAL( surf_usm_h%z0 ) ) THEN 741 terminate_run_l = .TRUE. 742 ENDIF 743 ENDIF 744 ! 745 ! Check roughness homogeneity between differt surface types. 746 IF ( surf_lsm_h%ns >= 1 .AND. surf_def_h(0)%ns >= 1 ) THEN 747 IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h ) .OR. & 748 MINVAL( surf_lsm_h%z0 ) /= MAXVAL( surf_def_h(0)%z0 ) ) THEN 749 terminate_run_l = .TRUE. 750 ENDIF 751 ENDIF 752 IF ( surf_usm_h%ns >= 1 .AND. surf_def_h(0)%ns >= 1 ) THEN 753 IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h ) .OR. & 754 MINVAL( surf_usm_h%z0 ) /= MAXVAL( surf_def_h(0)%z0 ) ) THEN 755 terminate_run_l = .TRUE. 756 ENDIF 757 ENDIF 758 IF ( surf_usm_h%ns >= 1 .AND. surf_lsm_h%ns >= 1 ) THEN 759 IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h ) .OR. & 760 MINVAL( surf_usm_h%z0 ) /= MAXVAL( surf_lsm_h%z0 ) ) THEN 761 terminate_run_l = .TRUE. 762 ENDIF 763 ENDIF 764 368 765 369 766 #if defined( __parallel ) … … 389 786 ! 390 787 ! Use the lowest possible value for z_mo 391 k = MINVAL(nzb_s_inner)788 k = nzb 392 789 z_mo = zu(k+1)  zw(k) 393 790 … … 401 798 zeta_step = (zeta_max  zeta_stretch) / REAL(num_steps_n) 402 799 403 DO l = 1, num_steps_n404 zeta_tmp(l ) = zeta_tmp(l1) + zeta_step800 DO li = 1, num_steps_n 801 zeta_tmp(li) = zeta_tmp(li1) + zeta_step 405 802 ENDDO 406 803 … … 415 812 ! 416 813 ! Calculate z/L range from zeta_min to zeta_stretch 417 DO l = num_steps_n+1, num_steps1418 zeta_tmp(l ) = zeta_tmp(l1) + zeta_step814 DO li = num_steps_n+1, num_steps1 815 zeta_tmp(li) = zeta_tmp(li1) + zeta_step 419 816 zeta_step = zeta_step * regr 420 817 ENDDO … … 422 819 ! 423 820 ! Save roughness lengths to temporary variables 424 z0h_min = z0h(nys,nxl) 425 z0_min = z0(nys,nxl) 426 821 IF ( surf_def_h(0)%ns >= 1 ) THEN 822 z0h_min = surf_def_h(0)%z0h(1) 823 z0_min = surf_def_h(0)%z0(1) 824 ELSEIF ( surf_lsm_h%ns >= 1 ) THEN 825 z0h_min = surf_lsm_h%z0h(1) 826 z0_min = surf_lsm_h%z0(1) 827 ELSEIF ( surf_usm_h%ns >= 1 ) THEN 828 z0h_min = surf_usm_h%z0h(1) 829 z0_min = surf_usm_h%z0(1) 830 ENDIF 427 831 ! 428 832 ! Calculate lookup table for the Richardson number versus Obukhov length … … 430 834 ! boundary conditions for temperature 431 835 IF ( ibc_pt_b == 1 ) THEN 432 DO l = 0, num_steps1433 ol_tab(l ) =  z_mo / zeta_tmp(num_steps1l)434 rib_tab(l ) = z_mo / ol_tab(l) / ( LOG( z_mo / z0_min )&435  psi_m( z_mo / ol_tab(l ) ) &436 + psi_m( z0_min / ol_tab(l ) ) &836 DO li = 0, num_steps1 837 ol_tab(li) =  z_mo / zeta_tmp(num_steps1li) 838 rib_tab(li) = z_mo / ol_tab(li) / ( LOG( z_mo / z0_min ) & 839  psi_m( z_mo / ol_tab(li) ) & 840 + psi_m( z0_min / ol_tab(li) ) & 437 841 )**3 438 842 ENDDO 439 843 ELSE 440 DO l = 0, num_steps1441 ol_tab(l ) =  z_mo / zeta_tmp(num_steps1l)442 rib_tab(l ) = z_mo / ol_tab(l) * ( LOG( z_mo / z0h_min )&443  psi_h( z_mo / ol_tab(l ) )&444 + psi_h( z0h_min / ol_tab(l ) )&844 DO li = 0, num_steps1 845 ol_tab(li) =  z_mo / zeta_tmp(num_steps1li) 846 rib_tab(li) = z_mo / ol_tab(li) * ( LOG( z_mo / z0h_min ) & 847  psi_h( z_mo / ol_tab(li) ) & 848 + psi_h( z0h_min / ol_tab(li) ) & 445 849 ) & 446 850 / ( LOG( z_mo / z0_min ) & 447  psi_m( z_mo / ol_tab(l ) )&448 + psi_m( z0_min / ol_tab(l ) )&851  psi_m( z_mo / ol_tab(li) ) & 852 + psi_m( z0_min / ol_tab(li) ) & 449 853 )**2 450 854 ENDDO … … 467 871 !  468 872 !> Compute the absolute value of the horizontal velocity (relative to the 469 !> surface) . This is required by all methods873 !> surface) for horizontal surface elements. This is required by all methods. 470 874 !! 471 SUBROUTINE calc_uv _total875 SUBROUTINE calc_uvw_abs 472 876 473 877 IMPLICIT NONE 474 878 475 476 !$OMP PARALLEL DO PRIVATE( k ) 477 DO i = nxl, nxr 478 DO j = nys, nyn 479 480 k = nzb_s_inner(j,i) 481 uv_total(j,i) = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) & 482  u(k,j,i)  u(k,j,i+1) ) )**2 + & 483 ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) & 484  v(k,j,i)  v(k,j+1,i) ) )**2 ) 485 486 ! 487 ! For too small values of the local wind, MOST does not work. A 488 ! threshold value is thus set if required 489 ! uv_total(j,i) = MAX(0.01_wp,uv_total(j,i)) 490 879 INTEGER(iwp) :: i !< running index x direction 880 INTEGER(iwp) :: ibit !< flag to mask computation of relative velocity in case of downwardfacing surfaces 881 INTEGER(iwp) :: j !< running index y direction 882 INTEGER(iwp) :: k !< running index z direction 883 INTEGER(iwp) :: m !< running index surface elements 884 885 ! 886 ! ibit is 1 for upwardfacing surfaces, zero for downwardfacing surfaces. 887 ibit = MERGE( 1, 0, .NOT. downward ) 888 889 DO m = 1, surf%ns 890 891 i = surf%i(m) 892 j = surf%j(m) 893 k = surf%k(m) 894 ! 895 ! Compute the absolute value of the horizontal velocity. 896 ! (relative to the surface in case the lower surface is the ocean). 897 ! Please note, in new surface modelling concept the index values changed, 898 ! i.e. the reference grid point is not the surfacegrid point itself but 899 ! the first grid point outside of the topography. 900 ! Note, in case of coupled oceanatmosphere simulations relative velocity 901 ! with respect to the ocean surface is used, hence, (k1,j,i) values 902 ! are used to calculate the absolute velocity. 903 ! However, this do not apply for downwardfacing walls. To mask this, 904 ! use ibit, which checks for upward/downwardfacing surfaces. 905 906 surf%uvw_abs(m) = SQRT( & 907 ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) & 908  ( u(k1,j,i) + u(k1,j,i+1) & 909 ) * ibit & 910 ) & 911 )**2 + & 912 ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) & 913  ( v(k1,j,i) + v(k1,j+1,i) & 914 ) * ibit & 915 ) & 916 )**2 & 917 ) 918 919 ENDDO 920 921 END SUBROUTINE calc_uvw_abs 922 923 924 !! 925 ! Description: 926 !  927 !> Compute the absolute value of the horizontal velocity (relative to the 928 !> surface) for horizontal surface elements. This is required by all methods. 929 !! 930 SUBROUTINE calc_uvw_abs_v_ugrid 931 932 IMPLICIT NONE 933 934 INTEGER(iwp) :: i !< running index x direction 935 INTEGER(iwp) :: j !< running index y direction 936 INTEGER(iwp) :: k !< running index z direction 937 INTEGER(iwp) :: m !< running index surface elements 938 939 REAL(wp) :: u_i 940 REAL(wp) :: w_i 941 942 943 DO m = 1, surf%ns 944 i = surf%i(m) 945 j = surf%j(m) 946 k = surf%k(m) 947 ! 948 ! Compute the absolute value of the surface parallel velocity on ugrid. 949 u_i = u(k,j,i) 950 w_i = 0.25_wp * ( w(k1,j,i1) + w(k1,j,i) + & 951 w(k,j,i1) + w(k,j,i) ) 952 953 surf%uvw_abs(m) = SQRT( u_i**2 + w_i**2 ) 954 955 ENDDO 956 957 END SUBROUTINE calc_uvw_abs_v_ugrid 958 959 !! 960 ! Description: 961 !  962 !> Compute the absolute value of the horizontal velocity (relative to the 963 !> surface) for horizontal surface elements. This is required by all methods. 964 !! 965 SUBROUTINE calc_uvw_abs_v_vgrid 966 967 IMPLICIT NONE 968 969 INTEGER(iwp) :: i !< running index x direction 970 INTEGER(iwp) :: j !< running index y direction 971 INTEGER(iwp) :: k !< running index z direction 972 INTEGER(iwp) :: m !< running index surface elements 973 974 REAL(wp) :: v_i 975 REAL(wp) :: w_i 976 977 978 DO m = 1, surf%ns 979 i = surf%i(m) 980 j = surf%j(m) 981 k = surf%k(m) 982 983 v_i = u(k,j,i) 984 w_i = 0.25_wp * ( w(k1,j1,i) + w(k1,j,i) + & 985 w(k,j1,i) + w(k,j,i) ) 986 987 surf%uvw_abs(m) = SQRT( v_i**2 + w_i**2 ) 988 989 ENDDO 990 991 END SUBROUTINE calc_uvw_abs_v_vgrid 992 993 !! 994 ! Description: 995 !  996 !> Compute the absolute value of the horizontal velocity (relative to the 997 !> surface) for horizontal surface elements. This is required by all methods. 998 !! 999 SUBROUTINE calc_uvw_abs_v_wgrid 1000 1001 IMPLICIT NONE 1002 1003 INTEGER(iwp) :: i !< running index x direction 1004 INTEGER(iwp) :: j !< running index y direction 1005 INTEGER(iwp) :: k !< running index z direction 1006 INTEGER(iwp) :: m !< running index surface elements 1007 1008 REAL(wp) :: u_i 1009 REAL(wp) :: v_i 1010 REAL(wp) :: w_i 1011 ! 1012 ! North (l=0) and southfacing (l=1) surfaces 1013 IF ( l == 0 .OR. l == 1 ) THEN 1014 DO m = 1, surf%ns 1015 i = surf%i(m) 1016 j = surf%j(m) 1017 k = surf%k(m) 1018 1019 u_i = 0.25_wp * ( u(k+1,j,i+1) + u(k+1,j,i) + & 1020 u(k,j,i+1) + u(k,j,i) ) 1021 v_i = 0.0_wp 1022 w_i = w(k,j,i) 1023 1024 surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 491 1025 ENDDO 492 ENDDO 493 494 ! 495 ! Values of uv_total need to be exchanged at the ghost boundaries 496 CALL exchange_horiz_2d( uv_total ) 497 498 END SUBROUTINE calc_uv_total 1026 ! 1027 ! East (l=2) and westfacing (l=3) surfaces 1028 ELSE 1029 DO m = 1, surf%ns 1030 i = surf%i(m) 1031 j = surf%j(m) 1032 k = surf%k(m) 1033 1034 u_i = 0.0_wp 1035 v_i = 0.25_wp * ( v(k+1,j+1,i) + v(k+1,j,i) + & 1036 v(k,j+1,i) + v(k,j,i) ) 1037 w_i = w(k,j,i) 1038 1039 surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 1040 ENDDO 1041 ENDIF 1042 1043 END SUBROUTINE calc_uvw_abs_v_wgrid 1044 1045 !! 1046 ! Description: 1047 !  1048 !> Compute the absolute value of the horizontal velocity (relative to the 1049 !> surface) for horizontal surface elements. This is required by all methods. 1050 !! 1051 SUBROUTINE calc_uvw_abs_v_sgrid 1052 1053 IMPLICIT NONE 1054 1055 INTEGER(iwp) :: i !< running index x direction 1056 INTEGER(iwp) :: j !< running index y direction 1057 INTEGER(iwp) :: k !< running index z direction 1058 INTEGER(iwp) :: m !< running index surface elements 1059 1060 REAL(wp) :: u_i 1061 REAL(wp) :: v_i 1062 REAL(wp) :: w_i 1063 1064 ! 1065 ! North (l=0) and southfacing (l=1) walls 1066 IF ( l == 0 .OR. l == 1 ) THEN 1067 DO m = 1, surf%ns 1068 i = surf%i(m) 1069 j = surf%j(m) 1070 k = surf%k(m) 1071 1072 u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) 1073 v_i = 0.0_wp 1074 w_i = 0.5_wp * ( w(k,j,i) + w(k1,j,i) ) 1075 1076 surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 1077 ENDDO 1078 ! 1079 ! East (l=2) and westfacing (l=3) walls 1080 ELSE 1081 DO m = 1, surf%ns 1082 i = surf%i(m) 1083 j = surf%j(m) 1084 k = surf%k(m) 1085 1086 u_i = 0.0_wp 1087 v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) 1088 w_i = 0.5_wp * ( w(k,j,i) + w(k1,j,i) ) 1089 1090 surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 1091 ENDDO 1092 ENDIF 1093 1094 END SUBROUTINE calc_uvw_abs_v_sgrid 499 1095 500 1096 … … 508 1104 IMPLICIT NONE 509 1105 510 INTEGER(iwp) :: iter, & !< Newton iteration step 511 l !< look index 512 513 REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) :: rib !< Bulk Richardson number 1106 INTEGER(iwp) :: iter !< Newton iteration step 1107 INTEGER(iwp) :: li !< look index 1108 INTEGER(iwp) :: m !< loop variable over all horizontal wall elements 514 1109 515 1110 REAL(wp) :: f, & !< Function for Newton iteration: f = Ri  [...]/[...]^2 = 0 … … 521 1116 522 1117 IF ( TRIM( most_method ) /= 'circular' ) THEN 523 524 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 525 DO i = nxl, nxr 526 DO j = nys, nyn 527 528 k = nzb_s_inner(j,i) 529 z_mo = zu(k+1)  zw(k) 530 531 ! 532 ! Evaluate bulk Richardson number (calculation depends on 533 ! definition based on setting of boundary conditions 534 IF ( ibc_pt_b /= 1 ) THEN 535 IF ( humidity ) THEN 536 rib(j,i) = g * z_mo * ( vpt(k+1,j,i)  vpt(k,j,i) ) & 537 / ( uv_total(j,i)**2 * vpt(k+1,j,i) + 1.0E20_wp ) 538 ELSE 539 rib(j,i) = g * z_mo * ( pt(k+1,j,i)  pt(k,j,i) ) & 540 / ( uv_total(j,i)**2 * pt(k+1,j,i) + 1.0E20_wp ) 541 ENDIF 542 ELSE 543 ! 544 ! When using Neumann boundary conditions, the buoyancy flux 545 ! is required but cannot be calculated at the surface, as pt 546 ! and q are not known at the surface. Hence the values at 547 ! first grid level are used to estimate the buoyancy flux 548 IF ( humidity ) THEN 549 rib(j,i) =  g * z_mo * ( ( 1.0_wp + 0.61_wp & 550 * q(k+1,j,i) ) * shf(j,i) + 0.61_wp & 551 * pt(k+1,j,i) * qsws(j,i) ) * drho_air_zw(k) & 552 / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2& 553 + 1.0E20_wp) 554 ELSE 555 rib(j,i) =  g * z_mo * shf(j,i) * drho_air_zw(k) & 556 / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2 & 1118 ! 1119 ! Evaluate bulk Richardson number (calculation depends on 1120 ! definition based on setting of boundary conditions 1121 IF ( ibc_pt_b /= 1 ) THEN 1122 IF ( humidity ) THEN 1123 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1124 DO m = 1, surf%ns 1125 1126 i = surf%i(m) 1127 j = surf%j(m) 1128 k = surf%k(m) 1129 1130 z_mo = surf%z_mo(m) 1131 1132 surf%rib(m) = g * z_mo * & 1133 ( vpt(k,j,i)  vpt(k1,j,i) ) / & 1134 ( surf%uvw_abs(m)**2 * vpt(k,j,i) + 1.0E20_wp ) 1135 ENDDO 1136 ELSE 1137 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1138 DO m = 1, surf%ns 1139 1140 i = surf%i(m) 1141 j = surf%j(m) 1142 k = surf%k(m) 1143 1144 z_mo = surf%z_mo(m) 1145 1146 surf%rib(m) = g * z_mo * & 1147 ( pt(k,j,i)  pt(k1,j,i) ) / & 1148 ( surf%uvw_abs(m)**2 * pt(k,j,i) + 1.0E20_wp ) 1149 ENDDO 1150 ENDIF 1151 ELSE 1152 IF ( humidity ) THEN 1153 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1154 DO m = 1, surf%ns 1155 1156 i = surf%i(m) 1157 j = surf%j(m) 1158 k = surf%k(m) 1159 1160 z_mo = surf%z_mo(m) 1161 1162 surf%rib(m) =  g * z_mo * ( ( 1.0_wp + 0.61_wp & 1163 * q(k,j,i) ) * surf%shf(m) + 0.61_wp & 1164 * pt(k,j,i) * surf%qsws(m) ) * & 1165 drho_air_zw(k1) / & 1166 ( surf%uvw_abs(m)**3 * vpt(k,j,i) * kappa**2 & 557 1167 + 1.0E20_wp ) 558 ENDIF 559 ENDIF 560 561 ENDDO 562 ENDDO 1168 ENDDO 1169 ELSE 1170 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1171 DO m = 1, surf%ns 1172 1173 i = surf%i(m) 1174 j = surf%j(m) 1175 k = surf%k(m) 1176 1177 z_mo = surf%z_mo(m) 1178 1179 surf%rib(m) =  g * z_mo * surf%shf(m) * & 1180 drho_air_zw(k1) / & 1181 ( surf%uvw_abs(m)**3 * pt(k,j,i) * kappa**2 & 1182 + 1.0E20_wp ) 1183 ENDDO 1184 ENDIF 1185 ENDIF 563 1186 564 1187 ENDIF 1188 565 1189 566 1190 ! … … 569 1193 IF ( TRIM( most_method ) == 'newton' ) THEN 570 1194 571 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 572 DO i = nxl, nxr 573 DO j = nys, nyn 574 575 k = nzb_s_inner(j,i) 576 z_mo = zu(k+1)  zw(k) 577 578 ! 579 ! Store current value in case the Newton iteration fails 580 ol_old = ol(j,i) 1195 DO m = 1, surf%ns 1196 1197 i = surf%i(m) 1198 j = surf%j(m) 1199 k = surf%k(m) 1200 1201 z_mo = surf%z_mo(m) 1202 1203 ! 1204 ! Store current value in case the Newton iteration fails 1205 ol_old = surf%ol(m) 1206 1207 ! 1208 ! Ensure that the bulk Richardson number and the Obukhov 1209 ! length have the same sign 1210 IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. & 1211 ABS( surf%ol(m) ) == ol_max ) THEN 1212 IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp 1213 IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = 0.01_wp 1214 ENDIF 1215 ! 1216 ! Iteration to find Obukhov length 1217 iter = 0 1218 DO 1219 iter = iter + 1 1220 ! 1221 ! In case of divergence, use the value of the previous time step 1222 IF ( iter > 1000 ) THEN 1223 surf%ol(m) = ol_old 1224 EXIT 1225 ENDIF 1226 1227 ol_m = surf%ol(m) 1228 ol_l = ol_m  0.001_wp * ol_m 1229 ol_u = ol_m + 0.001_wp * ol_m 1230 1231 1232 IF ( ibc_pt_b /= 1 ) THEN 1233 ! 1234 ! Calculate f = Ri  [...]/[...]^2 = 0 1235 f = surf%rib(m)  ( z_mo / ol_m ) * ( & 1236 LOG( z_mo / surf%z0h(m) ) & 1237  psi_h( z_mo / ol_m ) & 1238 + psi_h( surf%z0h(m) / & 1239 ol_m ) & 1240 ) & 1241 / ( LOG( z_mo / surf%z0(m) ) & 1242  psi_m( z_mo / ol_m ) & 1243 + psi_m( surf%z0(m) / & 1244 ol_m ) & 1245 )**2 1246 1247 ! 1248 ! Calculate df/dL 1249 f_d_ol = (  ( z_mo / ol_u ) * ( LOG( z_mo / & 1250 surf%z0h(m) ) & 1251  psi_h( z_mo / ol_u ) & 1252 + psi_h( surf%z0h(m) / ol_u ) & 1253 ) & 1254 / ( LOG( z_mo / surf%z0(m) ) & 1255  psi_m( z_mo / ol_u ) & 1256 + psi_m( surf%z0(m) / ol_u ) & 1257 )**2 & 1258 + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) ) & 1259  psi_h( z_mo / ol_l ) & 1260 + psi_h( surf%z0h(m) / ol_l ) & 1261 ) & 1262 / ( LOG( z_mo / surf%z0(m) ) & 1263  psi_m( z_mo / ol_l ) & 1264 + psi_m( surf%z0(m) / ol_l ) & 1265 )**2 & 1266 ) / ( ol_u  ol_l ) 1267 ELSE 1268 ! 1269 ! Calculate f = Ri  1 /[...]^3 = 0 1270 f = surf%rib(m)  ( z_mo / ol_m ) / & 1271 ( LOG( z_mo / surf%z0(m) ) & 1272  psi_m( z_mo / ol_m ) & 1273 + psi_m( surf%z0(m) / ol_m ) & 1274 )**3 1275 1276 ! 1277 ! Calculate df/dL 1278 f_d_ol = (  ( z_mo / ol_u ) / ( LOG( z_mo / & 1279 surf%z0(m) ) & 1280  psi_m( z_mo / ol_u ) & 1281 + psi_m( surf%z0(m) / ol_u ) & 1282 )**3 & 1283 + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) ) & 1284  psi_m( z_mo / ol_l ) & 1285 + psi_m( surf%z0(m) / ol_l ) & 1286 )**3 & 1287 ) / ( ol_u  ol_l ) 1288 ENDIF 1289 ! 1290 ! Calculate new L 1291 surf%ol(m) = ol_m  f / f_d_ol 581 1292 582 1293 ! 583 1294 ! Ensure that the bulk Richardson number and the Obukhov 584 ! lengtH have the same sign 585 IF ( rib(j,i) * ol(j,i) < 0.0_wp .OR. & 586 ABS( ol(j,i) ) == ol_max ) THEN 587 IF ( rib(j,i) > 0.0_wp ) ol(j,i) = 0.01_wp 588 IF ( rib(j,i) < 0.0_wp ) ol(j,i) = 0.01_wp 1295 ! length have the same sign and ensure convergence. 1296 IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp 1297 ! 1298 ! If unrealistic value occurs, set L to the maximum 1299 ! value that is allowed 1300 IF ( ABS( surf%ol(m) ) > ol_max ) THEN 1301 surf%ol(m) = ol_max 1302 EXIT 589 1303 ENDIF 590 1304 ! 591 ! Iteration to find Obukhov length 592 iter = 0 593 DO 594 iter = iter + 1 595 ! 596 ! In case of divergence, use the value of the previous time step 597 IF ( iter > 1000 ) THEN 598 ol(j,i) = ol_old 599 EXIT 600 ENDIF 601 602 ol_m = ol(j,i) 603 ol_l = ol_m  0.001_wp * ol_m 604 ol_u = ol_m + 0.001_wp * ol_m 605 606 607 IF ( ibc_pt_b /= 1 ) THEN 608 ! 609 ! Calculate f = Ri  [...]/[...]^2 = 0 610 f = rib(j,i)  ( z_mo / ol_m ) * ( LOG( z_mo / z0h(j,i) )& 611  psi_h( z_mo / ol_m ) & 612 + psi_h( z0h(j,i) / ol_m ) & 613 ) & 614 / ( LOG( z_mo / z0(j,i) ) & 615  psi_m( z_mo / ol_m ) & 616 + psi_m( z0(j,i) / ol_m ) & 617 )**2 618 619 ! 620 ! Calculate df/dL 621 f_d_ol = (  ( z_mo / ol_u ) * ( LOG( z_mo / z0h(j,i) ) & 622  psi_h( z_mo / ol_u ) & 623 + psi_h( z0h(j,i) / ol_u ) & 624 ) & 625 / ( LOG( z_mo / z0(j,i) ) & 626  psi_m( z_mo / ol_u ) & 627 + psi_m( z0(j,i) / ol_u ) & 628 )**2 & 629 + ( z_mo / ol_l ) * ( LOG( z_mo / z0h(j,i) ) & 630  psi_h( z_mo / ol_l ) & 631 + psi_h( z0h(j,i) / ol_l ) & 632 ) & 633 / ( LOG( z_mo / z0(j,i) ) & 634  psi_m( z_mo / ol_l ) & 635 + psi_m( z0(j,i) / ol_l ) & 636 )**2 & 637 ) / ( ol_u  ol_l ) 638 ELSE 639 ! 640 ! Calculate f = Ri  1 /[...]^3 = 0 641 f = rib(j,i)  ( z_mo / ol_m ) / ( LOG( z_mo / z0(j,i) )& 642  psi_m( z_mo / ol_m ) & 643 + psi_m( z0(j,i) / ol_m ) & 644 )**3 645 646 ! 647 ! Calculate df/dL 648 f_d_ol = (  ( z_mo / ol_u ) / ( LOG( z_mo / z0(j,i) ) & 649  psi_m( z_mo / ol_u ) & 650 + psi_m( z0(j,i) / ol_u ) & 651 )**3 & 652 + ( z_mo / ol_l ) / ( LOG( z_mo / z0(j,i) ) & 653  psi_m( z_mo / ol_l ) & 654 + psi_m( z0(j,i) / ol_l ) & 655 )**3 & 656 ) / ( ol_u  ol_l ) 657 ENDIF 658 ! 659 ! Calculate new L 660 ol(j,i) = ol_m  f / f_d_ol 661 662 ! 663 ! Ensure that the bulk Richardson number and the Obukhov 664 ! length have the same sign and ensure convergence. 665 IF ( ol(j,i) * ol_m < 0.0_wp ) ol(j,i) = ol_m * 0.5_wp 666 667 ! 668 ! If unrealistic value occurs, set L to the maximum 669 ! value that is allowed 670 IF ( ABS( ol(j,i) ) > ol_max ) THEN 671 ol(j,i) = ol_max 672 EXIT 673 ENDIF 674 ! 675 ! Check for convergence 676 IF ( ABS( ( ol(j,i)  ol_m ) / ol(j,i) ) < 1.0E4_wp ) THEN 677 EXIT 678 ELSE 679 CYCLE 680 ENDIF 681 1305 ! Check for convergence 1306 IF ( ABS( ( surf%ol(m)  ol_m ) / & 1307 surf%ol(m) ) < 1.0E4_wp ) THEN 1308 EXIT 1309 ELSE 1310 CYCLE 1311 ENDIF 1312 1313 ENDDO 1314 ENDDO 1315 1316 ELSEIF ( TRIM( most_method ) == 'lookup' ) THEN 1317 1318 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo, li ) FIRSTPRIVATE( li_bnd ) LASTPRIVATE( li_bnd ) 1319 DO m = 1, surf%ns 1320 1321 i = surf%i(m) 1322 j = surf%j(m) 1323 k = surf%k(m) 1324 ! 1325 ! If the bulk Richardson number is outside the range of the lookup 1326 ! table, set it to the exceeding threshold value 1327 IF ( surf%rib(m) < rib_min ) surf%rib(m) = rib_min 1328 IF ( surf%rib(m) > rib_max ) surf%rib(m) = rib_max 1329 1330 ! 1331 ! Find the correct index bounds for linear interpolation. As the 1332 ! Richardson number will not differ very much from time step to 1333 ! time step , use the index from the last step and search in the 1334 ! correct direction 1335 li = li_bnd 1336 IF ( rib_tab(li)  surf%rib(m) > 1.0_wp ) THEN 1337 DO WHILE ( rib_tab(li1)  surf%rib(m) > 1.0_wp .AND. li > 1 ) 1338 li = li1 682 1339 ENDDO 683 684 ENDDO 1340 ELSE 1341 DO WHILE ( rib_tab(li)  surf%rib(m) < 0.0_wp & 1342 .AND. li < num_steps1 ) 1343 li = li+1 1344 ENDDO 1345 ENDIF 1346 li_bnd = li 1347 1348 ! 1349 ! Linear interpolation to find the correct value of z/L 1350 surf%ol(m) = ( ol_tab(li1) + ( ol_tab(li)  ol_tab(li1) ) & 1351 / ( rib_tab(li)  rib_tab(li1) ) & 1352 * ( surf%rib(m)  rib_tab(li1) ) ) 1353 685 1354 ENDDO 686 1355 687 ELSEIF ( TRIM( most_method ) == 'lookup' ) THEN688 689 !$OMP PARALLEL DO PRIVATE( k, l, z_mo ) FIRSTPRIVATE( l_bnd ) LASTPRIVATE( l_bnd )690 DO i = nxl, nxr691 DO j = nys, nyn692 693 !694 ! If the bulk Richardson number is outside the range of the lookup695 ! table, set it to the exceeding threshold value696 IF ( rib(j,i) < rib_min ) rib(j,i) = rib_min697 IF ( rib(j,i) > rib_max ) rib(j,i) = rib_max698 699 !700 ! Find the correct index bounds for linear interpolation. As the701 ! Richardson number will not differ very much from time step to702 ! time step , use the index from the last step and search in the703 ! correct direction704 l = l_bnd705 IF ( rib_tab(l)  rib(j,i) > 0.0_wp ) THEN706 DO WHILE ( rib_tab(l1)  rib(j,i) > 0.0_wp .AND. l > 0 )707 l = l1708 ENDDO709 ELSE710 DO WHILE ( rib_tab(l)  rib(j,i) < 0.0_wp &711 .AND. l < num_steps1 )712 l = l+1713 ENDDO714 ENDIF715 l_bnd = l716 717 !718 ! Linear interpolation to find the correct value of z/L719 ol(j,i) = ( ol_tab(l1) + ( ol_tab(l)  ol_tab(l1) ) &720 / ( rib_tab(l)  rib_tab(l1) ) &721 * ( rib(j,i)  rib_tab(l1) ) )722 723 ENDDO724 ENDDO725 726 1356 ELSEIF ( TRIM( most_method ) == 'circular' ) THEN 727 1357 728 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 729 DO i = nxl, nxr 730 DO j = nys, nyn 731 732 k = nzb_s_inner(j,i) 733 z_mo = zu(k+1)  zw(k) 734 735 IF ( .NOT. humidity ) THEN 736 ol(j,i) = ( pt(k+1,j,i) * us(j,i)**2 ) / ( kappa * g & 737 * ts(j,i) + 1E30_wp ) 738 ELSEIF ( cloud_physics ) THEN 739 740 ol(j,i) = ( vpt(k+1,j,i) * us(j,i)**2 ) / ( kappa * g & 741 * ( ts(j,i) + 0.61_wp * pt1(j,i) * qs(j,i) & 742 + 0.61_wp * qv1(j,i) * ts(j,i)  ts(j,i) & 743 * ql(k+1,j,i) ) + 1E30_wp ) 744 ELSE 745 ol(j,i) = ( vpt(k+1,j,i) * us(j,i)**2 ) / ( kappa * g & 746 * ( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i) & 747 + 0.61_wp * q(k+1,j,i) * ts(j,i) ) + 1E30_wp ) 748 ENDIF 1358 IF ( .NOT. humidity ) THEN 1359 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1360 DO m = 1, surf%ns 1361 1362 i = surf%i(m) 1363 j = surf%j(m) 1364 k = surf%k(m) 1365 1366 z_mo = surf%z_mo(m) 1367 1368 surf%ol(m) = ( pt(k,j,i) * surf%us(m)**2 ) / & 1369 ( kappa * g * & 1370 surf%ts(m) + 1E30_wp ) 749 1371 ! 750 1372 ! Limit the value range of the Obukhov length. 751 ! This is necessary for very small velocities (u,v > 0), because1373 ! This is necessary for very small velocities (u,v > 1), because 752 1374 ! the absolute value of ol can then become very small, which in 753 1375 ! consequence would result in very large shear stresses and very 754 1376 ! small momentum fluxes (both are generally unrealistic). 755 IF ( ( z_mo / ( ol(j,i) + 1E30_wp ) ) < zeta_min ) & 756 ol(j,i) = z_mo / zeta_min 757 IF ( ( z_mo / ( ol(j,i) + 1E30_wp ) ) > zeta_max ) & 758 ol(j,i) = z_mo / zeta_max 759 ENDDO 760 ENDDO 1377 IF ( ( z_mo / ( surf%ol(m) + 1E30_wp ) ) < zeta_min ) & 1378 surf%ol(m) = z_mo / zeta_min 1379 IF ( ( z_mo / ( surf%ol(m) + 1E30_wp ) ) > zeta_max ) & 1380 surf%ol(m) = z_mo / zeta_max 1381 1382 ENDDO 1383 ELSEIF ( cloud_physics ) THEN 1384 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1385 DO m = 1, surf%ns 1386 1387 i = surf%i(m) 1388 j = surf%j(m) 1389 k = surf%k(m) 1390 1391 z_mo = surf%z_mo(m) 1392 1393 1394 surf%ol(m) = ( vpt(k,j,i) * surf%us(m)**2 ) / & 1395 ( kappa * g * ( surf%ts(m) + & 1396 0.61_wp * surf%pt1(m) * surf%us(m) & 1397 + 0.61_wp * surf%qv1(m) * surf%ts(m)  & 1398 surf%ts(m) * ql(k,j,i) ) + 1E30_wp ) 1399 ! 1400 ! Limit the value range of the Obukhov length. 1401 ! This is necessary for very small velocities (u,v > 1), because 1402 ! the absolute value of ol can then become very small, which in 1403 ! consequence would result in very large shear stresses and very 1404 ! small momentum fluxes (both are generally unrealistic). 1405 IF ( ( z_mo / ( surf%ol(m) + 1E30_wp ) ) < zeta_min ) & 1406 surf%ol(m) = z_mo / zeta_min 1407 IF ( ( z_mo / ( surf%ol(m) + 1E30_wp ) ) > zeta_max ) & 1408 surf%ol(m) = z_mo / zeta_max 1409 1410 ENDDO 1411 ELSE 1412 1413 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1414 DO m = 1, surf%ns 1415 1416 i = surf%i(m) 1417 j = surf%j(m) 1418 k = surf%k(m) 1419 1420 z_mo = surf%z_mo(m) 1421 1422 surf%ol(m) = ( vpt(k,j,i) * surf%us(m)**2 ) / & 1423 ( kappa * g * ( surf%ts(m) + 0.61_wp * pt(k,j,i) * & 1424 surf%qs(m) + 0.61_wp * q(k,j,i) * & 1425 surf%ts(m) ) + 1E30_wp ) 1426 1427 ! 1428 ! Limit the value range of the Obukhov length. 1429 ! This is necessary for very small velocities (u,v > 1), because 1430 ! the absolute value of ol can then become very small, which in 1431 ! consequence would result in very large shear stresses and very 1432 ! small momentum fluxes (both are generally unrealistic). 1433 IF ( ( z_mo / ( surf%ol(m) + 1E30_wp ) ) < zeta_min ) & 1434 surf%ol(m) = z_mo / zeta_min 1435 IF ( ( z_mo / ( surf%ol(m) + 1E30_wp ) ) > zeta_max ) & 1436 surf%ol(m) = z_mo / zeta_max 1437 1438 ENDDO 1439 1440 ENDIF 761 1441 762 1442 ENDIF 763 764 !765 ! Values of ol at ghost point locations are needed for the evaluation766 ! of usws and vsws.767 CALL exchange_horiz_2d( ol )768 1443 769 1444 END SUBROUTINE calc_ol … … 775 1450 IMPLICIT NONE 776 1451 777 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 778 DO i = nxlg, nxrg 779 DO j = nysg, nyng 780 781 k = nzb_s_inner(j,i)+1 782 z_mo = zu(k+1)  zw(k) 783 784 ! 785 ! Compute u* at the scalars' grid points 786 us(j,i) = kappa * uv_total(j,i) / ( LOG( z_mo / z0(j,i) ) & 787  psi_m( z_mo / ol(j,i) ) & 788 + psi_m( z0(j,i) / ol(j,i) ) ) 1452 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1453 1454 ! 1455 ! Compute u* at horizontal surfaces at the scalars' grid points 1456 IF ( .NOT. surf_vertical ) THEN 1457 ! 1458 ! Compute u* at upwardfacing surfaces 1459 IF ( .NOT. downward ) THEN 1460 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 1461 DO m = 1, surf%ns 1462 1463 k = surf%k(m) 1464 z_mo = surf%z_mo(m) 1465 ! 1466 ! Compute u* at the scalars' grid points 1467 surf%us(m) = kappa * surf%uvw_abs(m) / & 1468 ( LOG( z_mo / surf%z0(m) ) & 1469  psi_m( z_mo / surf%ol(m) ) & 1470 + psi_m( surf%z0(m) / surf%ol(m) ) ) 1471 1472 ENDDO 1473 ! 1474 ! Compute u* at downwardfacing surfaces. This case, do not consider 1475 ! any stability. 1476 ELSE 1477 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 1478 DO m = 1, surf%ns 1479 1480 k = surf%k(m) 1481 z_mo = surf%z_mo(m) 1482 ! 1483 ! Compute u* at the scalars' grid points 1484 surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) ) 1485 1486 ENDDO 1487 ENDIF 1488 ! 1489 ! Compute u* at vertical surfaces at the u/v/v grid, respectively. 1490 ! No stability is considered in this case. 1491 ELSE 1492 !$OMP PARALLEL DO PRIVATE( z_mo ) 1493 DO m = 1, surf%ns 1494 z_mo = surf%z_mo(m) 1495 1496 surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) ) 789 1497 ENDDO 790 END DO1498 ENDIF 791 1499 792 1500 END SUBROUTINE calc_us … … 794 1502 ! 795 1503 ! Calculate potential temperature and specific humidity at first grid level 1504 ! ( only for upwardfacing surfs ) 796 1505 SUBROUTINE calc_pt_q 797 1506 798 1507 IMPLICIT NONE 799 1508 800 DO i = nxlg, nxrg 801 DO j = nysg, nyng 802 k = nzb_s_inner(j,i)+1 803 pt1(j,i) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i) 804 qv1(j,i) = q(k,j,i)  ql(k,j,i) 805 ENDDO 1509 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1510 1511 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1512 DO m = 1, surf%ns 1513 1514 i = surf%i(m) 1515 j = surf%j(m) 1516 k = surf%k(m) 1517 1518 surf%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i) 1519 surf%qv1(m) = q(k,j,i)  ql(k,j,i) 1520 806 1521 ENDDO 807 1522 … … 814 1529 IMPLICIT NONE 815 1530 1531 1532 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1533 816 1534 ! 817 ! Compute theta* 818 IF ( constant_heatflux ) THEN 819 1535 ! Compute theta* at horizontal surfaces 1536 IF ( constant_heatflux .AND. .NOT. surf_vertical ) THEN 820 1537 ! 821 1538 ! For a given heat flux in the surface layer: 822 !$OMP PARALLEL DO 823 DO i = nxlg, nxrg 824 DO j = nysg, nyng 825 k = nzb_s_inner(j,i) 826 ts(j,i) = shf(j,i) * drho_air_zw(k) / ( us(j,i) + 1E30_wp ) 827 ! 828 ! ts must be limited, because otherwise overflow may occur in case 829 ! of us=0 when computing ol further below 830 IF ( ts(j,i) < 1.05E5_wp ) ts(j,i) = 1.0E5_wp 831 IF ( ts(j,i) > 1.0E5_wp ) ts(j,i) = 1.0E5_wp 832 ENDDO 1539 1540 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1541 DO m = 1, surf%ns 1542 1543 i = surf%i(m) 1544 j = surf%j(m) 1545 k = surf%k(m) 1546 1547 surf%ts(m) = surf%shf(m) * drho_air_zw(k1) / & 1548 ( surf%us(m) + 1E30_wp ) 1549 1550 ! 1551 ! ts must be limited, because otherwise overflow may occur in case 1552 ! of us=0 when computing ol further below 1553 IF ( surf%ts(m) < 1.05E5_wp ) surf%ts(m) = 1.0E5_wp 1554 IF ( surf%ts(m) > 1.0E5_wp ) surf%ts(m) = 1.0E5_wp 1555 833 1556 ENDDO 834 1557 835 ELSE 1558 ELSEIF ( .NOT. surf_vertical ) THEN 836 1559 ! 837 1560 ! For a given surface temperature: 838 1561 IF ( large_scale_forcing .AND. lsf_surf ) THEN 839 !$OMP PARALLEL DO 840 DO i = nxlg, nxrg 841 DO j = nysg, nyng 842 k = nzb_s_inner(j,i) 843 pt(k,j,i) = pt_surface 844 ENDDO 845 ENDDO 846 ENDIF 847 848 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 849 DO i = nxlg, nxrg 850 DO j = nysg, nyng 851 852 k = nzb_s_inner(j,i) 853 z_mo = zu(k+1)  zw(k) 854 855 IF ( cloud_physics ) THEN 856 ts(j,i) = kappa * ( pt1(j,i)  pt(k,j,i) ) & 857 / ( LOG( z_mo / z0h(j,i) ) & 858  psi_h( z_mo / ol(j,i) ) & 859 + psi_h( z0h(j,i) / ol(j,i) ) ) 860 ELSE 861 ts(j,i) = kappa * ( pt(k+1,j,i)  pt(k,j,i) ) & 862 / ( LOG( z_mo / z0h(j,i) ) & 863  psi_h( z_mo / ol(j,i) ) & 864 + psi_h( z0h(j,i) / ol(j,i) ) ) 865 ENDIF 866 867 ENDDO 1562 1563 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1564 DO m = 1, surf%ns 1565 i = surf%i(m) 1566 j = surf%j(m) 1567 k = surf%k(m) 1568 1569 pt(k1,j,i) = pt_surface 1570 ENDDO 1571 ENDIF 1572 1573 IF ( cloud_physics ) THEN 1574 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1575 DO m = 1, surf%ns 1576 1577 i = surf%i(m) 1578 j = surf%j(m) 1579 k = surf%k(m) 1580 1581 z_mo = surf%z_mo(m) 1582 1583 surf%ts(m) = kappa * ( surf%pt1(m)  pt(k1,j,i) ) & 1584 / ( LOG( z_mo / surf%z0h(m) ) & 1585  psi_h( z_mo / surf%ol(m) ) & 1586 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1587 1588 ENDDO 1589 ELSE 1590 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1591 DO m = 1, surf%ns 1592 1593 i = surf%i(m) 1594 j = surf%j(m) 1595 k = surf%k(m) 1596 1597 z_mo = surf%z_mo(m) 1598 1599 surf%ts(m) = kappa * ( pt(k,j,i)  pt(k1,j,i) ) & 1600 / ( LOG( z_mo / surf%z0h(m) ) & 1601  psi_h( z_mo / surf%ol(m) ) & 1602 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1603 ENDDO 1604 ENDIF 1605 ENDIF 1606 ! 1607 ! Compute theta* at vertical surfaces. This is only required in case of 1608 ! landsurface model, in order to compute aerodynamical resistance. 1609 IF ( surf_vertical ) THEN 1610 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1611 DO m = 1, surf%ns 1612 1613 i = surf%i(m) 1614 j = surf%j(m) 1615 k = surf%k(m) 1616 surf%ts(m) = surf%shf(m) / ( surf%us(m) + 1E30_wp ) 1617 ! 1618 ! ts must be limited, because otherwise overflow may occur in case 1619 ! of us=0 when computing ol further below 1620 IF ( surf%ts(m) < 1.05E5_wp ) surf%ts(m) = 1.0E5_wp 1621 IF ( surf%ts(m) > 1.0E5_wp ) surf%ts(m) = 1.0E5_wp 1622 868 1623 ENDDO 869 1624 ENDIF 870 1625 871 1626 ! 872 ! If required compute q* 1627 ! If required compute q* at horizontal surfaces 873 1628 IF ( humidity ) THEN 874 IF ( constant_waterflux ) THEN1629 IF ( constant_waterflux .AND. .NOT. surf_vertical ) THEN 875 1630 ! 876 1631 ! For a given water flux in the surface layer 877 !$OMP PARALLEL DO 878 DO i = nxlg, nxrg 879 DO j = nysg, nyng 880 k = nzb_s_inner(j,i) 881 qs(j,i) = qsws(j,i) * drho_air_zw(k) / ( us(j,i) + 1E30_wp ) 882 ENDDO 883 ENDDO 884 885 ELSE 1632 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1633 DO m = 1, surf%ns 1634 1635 i = surf%i(m) 1636 j = surf%j(m) 1637 k = surf%k(m) 1638 surf%qs(m) = surf%qsws(m) * drho_air_zw(k1) / & 1639 ( surf%us(m) + 1E30_wp ) 1640 1641 ENDDO 1642 1643 ELSEIF ( .NOT. surf_vertical ) THEN 886 1644 coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. & 887 1645 run_coupled ) 888 1646 889 1647 IF ( large_scale_forcing .AND. lsf_surf ) THEN 890 !$OMP PARALLEL DO 891 DO i = nxlg, nxrg 892 DO j = nysg, nyng 893 k = nzb_s_inner(j,i) 894 q(k,j,i) = q_surface 895 ENDDO 1648 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1649 DO m = 1, surf%ns 1650 1651 i = surf%i(m) 1652 j = surf%j(m) 1653 k = surf%k(m) 1654 q(k1,j,i) = q_surface 1655 896 1656 ENDDO 897 1657 ENDIF 898 1658 899 !$OMP PARALLEL DO PRIVATE( e_s, k, z_mo ) 900 DO i = nxlg, nxrg 901 DO j = nysg, nyng 902 903 k = nzb_s_inner(j,i) 904 z_mo = zu(k+1)  zw(k) 905 906 ! 907 ! Assume saturation for atmosphere coupled to ocean (but not 908 ! in case of precursor runs) 909 IF ( coupled_run ) THEN 910 e_s = 6.1_wp * & 911 EXP( 0.07_wp * ( MIN(pt(k,j,i),pt(k+1,j,i)) & 1659 ! 1660 ! Assume saturation for atmosphere coupled to ocean (but not 1661 ! in case of precursor runs) 1662 IF ( coupled_run ) THEN 1663 !$OMP PARALLEL DO PRIVATE( i, j, k, e_s ) 1664 DO m = 1, surf%ns 1665 i = surf%i(m) 1666 j = surf%j(m) 1667 k = surf%k(m) 1668 e_s = 6.1_wp * & 1669 EXP( 0.07_wp * ( MIN(pt(k1,j,i),pt(k,j,i)) & 912 1670  273.15_wp ) ) 913 q(k,j,i) = 0.622_wp * e_s / ( surface_pressure  e_s ) 914 ENDIF 915 916 IF ( cloud_physics ) THEN 917 qs(j,i) = kappa * ( qv1(j,i)  q(k,j,i) ) & 918 / ( LOG( z_mo / z0q(j,i) ) & 919  psi_h( z_mo / ol(j,i) ) & 920 + psi_h( z0q(j,i) / ol(j,i) ) ) 921 922 ELSE 923 qs(j,i) = kappa * ( q(k+1,j,i)  q(k,j,i) ) & 924 / ( LOG( z_mo / z0q(j,i) ) & 925  psi_h( z_mo / ol(j,i) ) & 926 + psi_h( z0q(j,i) / ol(j,i) ) ) 927 ENDIF 928 1671 q(k1,j,i) = 0.622_wp * e_s / ( surface_pressure  e_s ) 929 1672 ENDDO 1673 ENDIF 1674 1675 IF ( cloud_physics ) THEN 1676 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1677 DO m = 1, surf%ns 1678 1679 i = surf%i(m) 1680 j = surf%j(m) 1681 k = surf%k(m) 1682 1683 z_mo = surf%z_mo(m) 1684 1685 surf%qs(m) = kappa * ( surf%qv1(m)  q(k1,j,i) ) & 1686 / ( LOG( z_mo / surf%z0q(m) ) & 1687  psi_h( z_mo / surf%ol(m) ) & 1688 + psi_h( surf%z0q(m) / & 1689 surf%ol(m) ) ) 1690 ENDDO 1691 ELSE 1692 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1693 DO m = 1, surf%ns 1694 1695 i = surf%i(m) 1696 j = surf%j(m) 1697 k = surf%k(m) 1698 1699 z_mo = surf%z_mo(m) 1700 1701 surf%qs(m) = kappa * ( q(k,j,i)  q(k1,j,i) ) & 1702 / ( LOG( z_mo / surf%z0q(m) ) & 1703  psi_h( z_mo / surf%ol(m) ) & 1704 + psi_h( surf%z0q(m) / & 1705 surf%ol(m) ) ) 1706 ENDDO 1707 ENDIF 1708 ENDIF 1709 ! 1710 ! Compute q* at vertical surfaces 1711 IF ( surf_vertical ) THEN 1712 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1713 DO m = 1, surf%ns 1714 1715 i = surf%i(m) 1716 j = surf%j(m) 1717 k = surf%k(m) 1718 surf%qs(m) = surf%qsws(m) / ( surf%us(m) + 1E30_wp ) 1719 930 1720 ENDDO 931 1721 ENDIF … … 935 1725 ! If required compute s* 936 1726 IF ( passive_scalar ) THEN 937 IF ( constant_scalarflux ) THEN 938 ! 939 ! For a given water flux in the surface layer 940 !$OMP PARALLEL DO 941 DO i = nxlg, nxrg 942 DO j = nysg, nyng 943 ss(j,i) = ssws(j,i) / ( us(j,i) + 1E30_wp ) 944 ENDDO 1727 ! 1728 ! At horizontal surfaces 1729 IF ( constant_scalarflux .AND. .NOT. surf_vertical ) THEN 1730 ! 1731 ! For a given scalar flux in the surface layer 1732 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1733 DO m = 1, surf%ns 1734 i = surf%i(m) 1735 j = surf%j(m) 1736 k = surf%k(m) 1737 surf%ss(m) = surf%ssws(m) / ( surf%us(m) + 1E30_wp ) 1738 ENDDO 1739 ENDIF 1740 ! 1741 ! At vertical surfaces 1742 IF ( surf_vertical ) THEN 1743 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1744 DO m = 1, surf%ns 1745 i = surf%i(m) 1746 j = surf%j(m) 1747 k = surf%k(m) 1748 surf%ss(m) = surf%ssws(m) / ( surf%us(m) + 1E30_wp ) 945 1749 ENDDO 946 1750 ENDIF … … 950 1754 ! 951 1755 ! If required compute qr* and nr* 952 IF ( cloud_physics .AND. microphysics_seifert )&953 THEN954 955 !$OMP PARALLEL DO PRIVATE( k, z_mo )956 DO i = nxlg, nxrg 957 DO j = nysg, nyng958 959 k = nzb_s_inner(j,i)960 z_mo = zu(k+1)  zw(k) 961 962 qrs(j,i) = kappa * ( qr(k+1,j,i)  qr(k,j,i) ) & 963 / ( LOG( z_mo / z0q(j,i) ) &964  psi_h( z_mo / ol(j,i) )&965 + psi_h( z0q(j,i) / ol(j,i) ) )966 967 nrs(j,i) = kappa * ( nr(k+1,j,i)  nr(k,j,i) ) & 968 / ( LOG( z_mo / z0q(j,i) ) &969  psi_h( z_mo / ol(j,i) )&970 + psi_h( z0q(j,i) / ol(j,i) ) )971 ENDDO1756 IF ( cloud_physics .AND. microphysics_seifert .AND. & 1757 .NOT. surf_vertical ) THEN 1758 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1759 DO m = 1, surf%ns 1760 1761 i = surf%i(m) 1762 j = surf%j(m) 1763 k = surf%k(m) 1764 1765 z_mo = surf%z_mo(m) 1766 1767 surf%qrs(m) = kappa * ( qr(k,j,i)  qr(k1,j,i) ) & 1768 / ( LOG( z_mo / surf%z0q(m) ) & 1769  psi_h( z_mo / surf%ol(m) ) & 1770 + psi_h( surf%z0q(m) / surf%ol(m) ) ) 1771 1772 surf%nrs(m) = kappa * ( nr(k,j,i)  nr(k1,j,i) ) & 1773 / ( LOG( z_mo / surf%z0q(m) ) & 1774  psi_h( z_mo / surf%ol(m) ) & 1775 + psi_h( surf%z0q(m) / surf%ol(m) ) ) 972 1776 ENDDO 973 1777 … … 984 1788 IMPLICIT NONE 985 1789 986 REAL(wp) :: ol_mid !< Gridinterpolated L 987 988 ! 989 ! Compute u'w' for the total model domain. 990 ! First compute the corresponding component of u* and square it. 991 !$OMP PARALLEL DO PRIVATE( k, ol_mid, z_mo ) 992 DO i = nxl, nxr 993 DO j = nys, nyn 994 995 k = nzb_u_inner(j,i) 996 z_mo = zu(k+1)  zw(k) 997 ! 998 ! Compute bulk Obukhov length for this point 999 ol_mid = 0.5_wp * ( ol(j,i1) + ol(j,i) ) 1000 1001 IF ( ol_mid == 0.0_wp ) THEN 1002 ol_mid = MIN(ol(j,i1), ol(j,i)) 1790 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1791 1792 REAL(wp) :: dum !< dummy to precalculate logarithm 1793 REAL(wp) :: flag_u !< flag indicating ugrid, used for calculation of horizontal momentum fluxes at vertical surfaces 1794 REAL(wp) :: flag_v !< flag indicating vgrid, used for calculation of horizontal momentum fluxes at vertical surfaces 1795 REAL(wp), DIMENSION(:), ALLOCATABLE :: u_i !< ucomponent interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces 1796 REAL(wp), DIMENSION(:), ALLOCATABLE :: v_i !< vcomponent interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces 1797 REAL(wp), DIMENSION(:), ALLOCATABLE :: w_i !< wcomponent interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces 1798 1799 ! 1800 ! Calcuate surface fluxes at horizontal walls 1801 IF ( .NOT. surf_vertical ) THEN 1802 ! 1803 ! Compute u'w' for the total model domain at upwardfacing surfaces. 1804 ! First compute the corresponding component of u* and square it. 1805 IF ( .NOT. downward ) THEN 1806 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1807 DO m = 1, surf%ns 1808 1809 i = surf%i(m) 1810 j = surf%j(m) 1811 k = surf%k(m) 1812 1813 z_mo = surf%z_mo(m) 1814 1815 surf%usws(m) = kappa * ( u(k,j,i)  u(k1,j,i) ) & 1816 / ( LOG( z_mo / surf%z0(m) ) & 1817  psi_m( z_mo / surf%ol(m) ) & 1818 + psi_m( surf%z0(m) / surf%ol(m) ) ) 1819 ! 1820 ! Please note, the computation of usws is not fully accurate. Actually 1821 ! a further interpolation of us onto the ugrid, where usws is defined, 1822 ! is required. However, this is not done as this would require several 1823 ! data transfers between 2Dgrid and the surftype. 1824 ! The impact of the missing interpolation is negligible as several 1825 ! tests had shown. 1826 ! Same also for ol. 1827 surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k1) 1828 1829 ENDDO 1830 ! 1831 ! At downwardfacing surfaces 1832 ELSE 1833 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1834 DO m = 1, surf%ns 1835 1836 i = surf%i(m) 1837 j = surf%j(m) 1838 k = surf%k(m) 1839 1840 z_mo = surf%z_mo(m) 1841 1842 surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) ) 1843 ! 1844 ! To Do: Is the sign correct??? 1845 surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k) 1846 1847 ENDDO 1848 ENDIF 1849 1850 ! 1851 ! Compute v'w' for the total model domain. 1852 ! First compute the corresponding component of u* and square it. 1853 ! Upwardfacing surfaces 1854 IF ( .NOT. downward ) THEN 1855 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1856 DO m = 1, surf%ns 1857 i = surf%i(m) 1858 j = surf%j(m) 1859 k = surf%k(m) 1860 1861 z_mo = surf%z_mo(m) 1862 1863 surf%vsws(m) = kappa * ( v(k,j,i)  v(k1,j,i) ) & 1864 / ( LOG( z_mo / surf%z0(m) ) & 1865  psi_m( z_mo / surf%ol(m) ) & 1866 + psi_m( surf%z0(m) / surf%ol(m) ) ) 1867 ! 1868 ! Please note, the computation of vsws is not fully accurate. Actually 1869 ! a further interpolation of us onto the vgrid, where vsws is defined, 1870 ! is required. However, this is not done as this would require several 1871 ! data transfers between 2Dgrid and the surftype. 1872 ! The impact of the missing interpolation is negligible as several 1873 ! tests had shown. 1874 ! Same also for ol. 1875 surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k1) 1876 ENDDO 1877 ! 1878 ! Downwardfacing surfaces 1879 ELSE 1880 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1881 DO m = 1, surf%ns 1882 i = surf%i(m) 1883 j = surf%j(m) 1884 k = surf%k(m) 1885 1886 z_mo = surf%z_mo(m) 1887 1888 surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) ) 1889 1890 surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k) 1891 ENDDO 1892 ENDIF 1893 ! 1894 ! Compute the vertical kinematic heat flux 1895 IF ( .NOT. constant_heatflux .AND. ( simulated_time <= & 1896 skip_time_do_lsm .OR. .NOT. land_surface ) .AND. & 1897 .NOT. urban_surface .AND. .NOT. downward ) THEN 1898 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1899 DO m = 1, surf%ns 1900 i = surf%i(m) 1901 j = surf%j(m) 1902 k = surf%k(m) 1903 surf%shf(m) = surf%ts(m) * surf%us(m) * rho_air_zw(k1) 1904 ENDDO 1905 ENDIF 1906 ! 1907 ! Compute the vertical water flux 1908 IF ( .NOT. constant_waterflux .AND. humidity .AND. & 1909 ( simulated_time <= skip_time_do_lsm & 1910 .OR. .NOT. land_surface ) .AND. .NOT. downward ) THEN 1911 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1912 DO m = 1, surf%ns 1913 i = surf%i(m) 1914 j = surf%j(m) 1915 k = surf%k(m) 1916 surf%qsws(m) = surf%qs(m) * surf%us(m) * rho_air_zw(k1) 1917 ENDDO 1918 ENDIF 1919 ! 1920 ! Compute the vertical scalar flux 1921 IF ( .NOT. constant_scalarflux .AND. passive_scalar .AND. & 1922 .NOT. downward ) THEN 1923 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1924 DO m = 1, surf%ns 1925 1926 i = surf%i(m) 1927 j = surf%j(m) 1928 k = surf%k(m) 1929 surf%ssws(m) = surf%ss(m) * surf%us(m) 1930 1931 ENDDO 1932 ENDIF 1933 ! 1934 ! Compute (turbulent) fluxes of rain water content and rain drop conc. 1935 IF ( cloud_physics .AND. microphysics_seifert .AND. & 1936 .NOT. downward) THEN 1937 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1938 DO m = 1, surf%ns 1939 1940 i = surf%i(m) 1941 j = surf%j(m) 1942 k = surf%k(m) 1943 1944 surf%qrsws(m) = surf%qrs(m) * surf%us(m) 1945 surf%nrsws(m) = surf%nrs(m) * surf%us(m) 1946 ENDDO 1947 ENDIF 1948 1949 ! 1950 ! Bottom boundary condition for the TKE. 1951 IF ( ibc_e_b == 2 ) THEN 1952 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1953 DO m = 1, surf%ns 1954 1955 i = surf%i(m) 1956 j = surf%j(m) 1957 k = surf%k(m) 1958 1959 e(k,j,i) = ( surf%us(m) / 0.1_wp )**2 1960 ! 1961 ! As a test: cm = 0.4 1962 ! e(k,j,i) = ( us(j,i) / 0.4_wp )**2 1963 e(k1,j,i) = e(k,j,i) 1964 1965 ENDDO 1966 ENDIF 1967 ! 1968 ! Calcuate surface fluxes at vertical surfaces. No stability is considered. 1969 ELSE 1970 ! 1971 ! Compute usvs l={0,1} and vsus l={2,3} 1972 IF ( mom_uv ) THEN 1973 ! 1974 ! Generalize computation by introducing flags. At north and south 1975 ! facing surfaces ucomponent is used, at east and westfacing 1976 ! surfaces vcomponent is used. 1977 flag_u = MERGE( 1.0_wp, 0.0_wp, l == 0 .OR. l == 1 ) 1978 flag_v = MERGE( 1.0_wp, 0.0_wp, l == 2 .OR. l == 3 ) 1979 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1980 DO m = 1, surf%ns 1981 i = surf%i(m) 1982 j = surf%j(m) 1983 k = surf%k(m) 1984 1985 z_mo = surf%z_mo(m) 1986 1987 surf%mom_flux_uv(m) = kappa * & 1988 ( flag_u * u(k,j,i) + flag_v * v(k,j,i) ) / & 1989 LOG( z_mo / surf%z0(m) ) 1990 1991 surf%mom_flux_uv(m) = & 1992  surf%mom_flux_uv(m) * surf%us(m) 1993 ENDDO 1994 ENDIF 1995 ! 1996 ! Compute wsus l={0,1} and wsvs l={2,3} 1997 IF ( mom_w ) THEN 1998 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1999 DO m = 1, surf%ns 2000 i = surf%i(m) 2001 j = surf%j(m) 2002 k = surf%k(m) 2003 2004 z_mo = surf%z_mo(m) 2005 2006 surf%mom_flux_w(m) = kappa * w(k,j,i) / LOG( z_mo / surf%z0(m) ) 2007 2008 surf%mom_flux_w(m) = & 2009  surf%mom_flux_w(m) * surf%us(m) 2010 ENDDO 2011 ENDIF 2012 ! 2013 ! Compute momentum fluxes used for subgridscale TKE production at 2014 ! vertical surfaces. In constrast to the calculated momentum fluxes at 2015 ! vertical surfaces before, which are defined on the u/v/wgrid, 2016 ! respectively), the TKE fluxes are defined at the scalar grid. 2017 ! 2018 IF ( mom_tke ) THEN 2019 ! 2020 ! Precalculate velocity components at scalar grid point. 2021 ALLOCATE( u_i(1:surf%ns) ) 2022 ALLOCATE( v_i(1:surf%ns) ) 2023 ALLOCATE( w_i(1:surf%ns) ) 2024 2025 IF ( l == 0 .OR. l == 1 ) THEN 2026 !$OMP PARALLEL DO PRIVATE( i, j, k ) 2027 DO m = 1, surf%ns 2028 i = surf%i(m) 2029 j = surf%j(m) 2030 k = surf%k(m) 2031 2032 u_i(m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) 2033 v_i(m) = 0.0_wp 2034 w_i(m) = 0.5_wp * ( w(k,j,i) + w(k1,j,i) ) 2035 ENDDO 2036 ELSE 2037 !$OMP PARALLEL DO PRIVATE( i, j, k ) 2038 DO m = 1, surf%ns 2039 i = surf%i(m) 2040 j = surf%j(m) 2041 k = surf%k(m) 2042 2043 u_i(m) = 0.0_wp 2044 v_i(m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) 2045 w_i(m) = 0.5_wp * ( w(k,j,i) + w(k1,j,i) ) 2046 ENDDO 1003 2047 ENDIF 1004 2048 1005 usws(j,i) = kappa * ( u(k+1,j,i)  u(k,j,i) ) & 1006 / ( LOG( z_mo / z0(j,i) ) & 1007  psi_m( z_mo / ol_mid ) & 1008 + psi_m( z0(j,i) / ol_mid ) ) 1009 1010 usws(j,i) = usws(j,i) * 0.5_wp * ( us(j,i1) + us(j,i) ) & 1011 * rho_air_zw(k) 1012 ENDDO 1013 ENDDO 1014 1015 ! 1016 ! Compute v'w' for the total model domain. 1017 ! First compute the corresponding component of u* and square it. 1018 !$OMP PARALLEL DO PRIVATE( k, ol_mid, z_mo ) 1019 DO i = nxl, nxr 1020 DO j = nys, nyn 1021 1022 k = nzb_v_inner(j,i) 1023 z_mo = zu(k+1)  zw(k) 1024 ! 1025 ! Compute bulk Obukhov length for this point 1026 ol_mid = 0.5_wp * ( ol(j1,i) + ol(j,i) ) 1027 1028 IF ( ol_mid == 0.0_wp ) THEN 1029 ol_mid = MIN(ol(j1,i), ol(j1,i)) 1030 ENDIF 1031 1032 vsws(j,i) = kappa * ( v(k+1,j,i)  v(k,j,i) ) & 1033 / ( LOG( z_mo / z0(j,i) ) & 1034  psi_m( z_mo / ol_mid ) & 1035 + psi_m( z0(j,i) / ol_mid ) ) 1036 1037 vsws(j,i) = vsws(j,i) * 0.5_wp * ( us(j1,i) + us(j,i) ) & 1038 * rho_air_zw(k) 1039 1040 ENDDO 1041 ENDDO 1042 1043 ! 1044 ! Exchange the boundaries for the momentum fluxes (is this still required?) 1045 CALL exchange_horiz_2d( usws ) 1046 CALL exchange_horiz_2d( vsws ) 1047 1048 ! 1049 ! Compute the vertical kinematic heat flux 1050 IF ( .NOT. constant_heatflux .AND. ( simulated_time <= & 1051 skip_time_do_lsm .OR. .NOT. land_surface ) .AND. & 1052 .NOT. urban_surface ) THEN 1053 !$OMP PARALLEL DO 1054 DO i = nxlg, nxrg 1055 DO j = nysg, nyng 1056 k = nzb_s_inner(j,i) 1057 shf(j,i) = ts(j,i) * us(j,i) * rho_air_zw(k) 1058 ENDDO 1059 ENDDO 1060 1061 ENDIF 1062 1063 ! 1064 ! Compute the vertical water flux 1065 IF ( .NOT. constant_waterflux .AND. humidity .AND. & 1066 ( simulated_time <= skip_time_do_lsm & 1067 .OR. .NOT. land_surface ) ) THEN 1068 !$OMP PARALLEL DO 1069 DO i = nxlg, nxrg 1070 DO j = nysg, nyng 1071 k = nzb_s_inner(j,i) 1072 qsws(j,i) = qs(j,i) * us(j,i) * rho_air_zw(k) 1073 ENDDO 1074 ENDDO 1075 1076 ENDIF 1077 1078 ! 1079 ! Compute the vertical scalar flux 1080 IF ( .NOT. constant_scalarflux .AND. passive_scalar .AND. & 1081 ( simulated_time <= skip_time_do_lsm & 1082 .OR. .NOT. land_surface ) ) THEN 1083 !$OMP PARALLEL DO 1084 DO i = nxlg, nxrg 1085 DO j = nysg, nyng 1086 ssws(j,i) = ss(j,i) * us(j,i) 1087 ENDDO 1088 ENDDO 1089 1090 ENDIF 1091 1092 ! 1093 ! Compute (turbulent) fluxes of rain water content and rain drop conc. 1094 IF ( cloud_physics .AND. microphysics_seifert ) THEN 1095 !$OMP PARALLEL DO 1096 DO i = nxlg, nxrg 1097 DO j = nysg, nyng 1098 qrsws(j,i) = qrs(j,i) * us(j,i) 1099 nrsws(j,i) = nrs(j,i) * us(j,i) 1100 ENDDO 1101 ENDDO 1102 ENDIF 1103 1104 ! 1105 ! Bottom boundary condition for the TKE 1106 IF ( ibc_e_b == 2 ) THEN 1107 !$OMP PARALLEL DO 1108 DO i = nxlg, nxrg 1109 DO j = nysg, nyng 1110 k = nzb_s_inner(j,i) 1111 e(k+1,j,i) = ( us(j,i) / 0.1_wp )**2 1112 ! 1113 ! As a test: cm = 0.4 1114 ! e(k+1,j,i) = ( us(j,i) / 0.4_wp )**2 1115 e(k,j,i) = e(k+1,j,i) 1116 ENDDO 1117 ENDDO 2049 !$OMP PARALLEL DO PRIVATE( i, j, k, dum, z_mo ) 2050 DO m = 1, surf%ns 2051 i = surf%i(m) 2052 j = surf%j(m) 2053 k = surf%k(m) 2054 2055 z_mo = surf%z_mo(m) 2056 2057 dum = kappa / LOG( z_mo / surf%z0(m) ) 2058 ! 2059 ! usvs (l=0,1) and vsus (l=2,3) 2060 surf%mom_flux_tke(0,m) = dum * ( u_i(m) + v_i(m) ) 2061 ! 2062 ! wsvs (l=0,1) and wsus (l=2,3) 2063 surf%mom_flux_tke(1,m) = dum * w_i(m) 2064 2065 surf%mom_flux_tke(0:1,m) = & 2066  surf%mom_flux_tke(0:1,m) * surf%us(m) 2067 ENDDO 2068 ! 2069 ! Deallocate temporary arrays 2070 DEALLOCATE( u_i ) 2071 DEALLOCATE( v_i ) 2072 DEALLOCATE( w_i ) 2073 ENDIF 1118 2074 ENDIF 1119 2075
Note: See TracChangeset
for help on using the changeset viewer.