Changeset 4039 for palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
 Timestamp:
 Jun 18, 2019 10:32:41 AM (4 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
r3998 r4039 25 25 !  26 26 ! $Id$ 27 !  Add output of uu, vv, ww to enable variance calculation according temporal 28 ! EC method 29 !  Allocate arrays only when they are required 30 !  Formatting adjustment 31 !  Rename subroutines 32 !  Further modularization 33 ! 34 ! 3998 20190523 13:38:11Z suehring 27 35 ! Bugfix in gathering all output strings 28 36 ! … … 43 51 MODULE diagnostic_output_quantities_mod 44 52 45 46 ! USE arrays_3d, & 47 ! ONLY: dzw, e, heatflux_output_conversion, nc, nr, p, pt, & 48 ! precipitation_amount, prr, q, qc, ql, ql_c, ql_v, qr, s, tend, & 49 ! u, v, vpt, w, zu, zw, waterflux_output_conversion, hyrho, d_exner 53 USE arrays_3d, & 54 ONLY: ddzu, u, v, w 50 55 ! 51 56 ! USE averaging … … 62 67 ! USE cpulog, & 63 68 ! ONLY: cpu_log, log_point 64 ! 65 ! USE indices, & 66 ! ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, & 67 ! nzb, nzt, wall_flags_0 69 70 USE grid_variables, & 71 ONLY: ddx, ddy 72 ! 73 USE indices, & 74 ONLY: nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0 68 75 ! 69 76 USE kinds … … 103 110 104 111 INTEGER(iwp) :: timestep_number_at_prev_calc = 0 !< ...at previous diagnostic output calculation 105 106 LOGICAL :: initialized_diagnostic_output_quantities = .FALSE. 107 LOGICAL :: prepared_diagnostic_output_quantities = .FALSE. 112 113 LOGICAL :: initialized_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is initialized 114 LOGICAL :: prepared_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is p 108 115 109 116 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti !< rotation(u,v,w) aka turbulence intensity 110 117 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti_av !< avg. rotation(u,v,w) aka turbulence intensity 118 119 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu !< uu 120 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu_av !< mean of uu 121 122 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv !< vv 123 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv_av !< mean of vv 124 125 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww !< ww 126 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww_av !< mean of ww 111 127 112 128 … … 117 133 ! 118 134 ! Public variables 119 PUBLIC do_all, & 120 initialized_diagnostic_output_quantities, & 121 prepared_diagnostic_output_quantities, & 122 ti, ti_av, & 123 timestep_number_at_prev_calc 124 ! 125 ! Public routines 126 PUBLIC diagnostic_output_quantities_calculate, & 127 diagnostic_output_quantities_init, & 128 diagnostic_output_quantities_prepare 129 130 131 INTERFACE diagnostic_output_quantities_init 132 MODULE PROCEDURE diagnostic_output_quantities_init 133 END INTERFACE diagnostic_output_quantities_init 134 135 INTERFACE diagnostic_output_quantities_calculate 136 MODULE PROCEDURE diagnostic_output_quantities_calculate 137 END INTERFACE diagnostic_output_quantities_calculate 138 139 INTERFACE diagnostic_output_quantities_prepare 140 MODULE PROCEDURE diagnostic_output_quantities_prepare 141 END INTERFACE diagnostic_output_quantities_prepare 135 PUBLIC do_all, & 136 initialized_diagnostic_output_quantities, & 137 prepared_diagnostic_output_quantities, & 138 timestep_number_at_prev_calc, & 139 ti_av, & 140 uu_av, & 141 vv_av, & 142 ww_av 143 ! 144 ! Public routines 145 PUBLIC doq_3d_data_averaging, & 146 doq_calculate, & 147 doq_check_data_output, & 148 doq_define_netcdf_grid, & 149 doq_output_2d, & 150 doq_output_3d, & 151 doq_output_mask, & 152 doq_init, & 153 doq_prepare, & 154 doq_wrd_local 155 ! doq_rrd_local, & 156 157 158 INTERFACE doq_3d_data_averaging 159 MODULE PROCEDURE doq_3d_data_averaging 160 END INTERFACE doq_3d_data_averaging 161 162 INTERFACE doq_calculate 163 MODULE PROCEDURE doq_calculate 164 END INTERFACE doq_calculate 165 166 INTERFACE doq_check_data_output 167 MODULE PROCEDURE doq_check_data_output 168 END INTERFACE doq_check_data_output 169 170 INTERFACE doq_define_netcdf_grid 171 MODULE PROCEDURE doq_define_netcdf_grid 172 END INTERFACE doq_define_netcdf_grid 173 174 INTERFACE doq_output_2d 175 MODULE PROCEDURE doq_output_2d 176 END INTERFACE doq_output_2d 177 178 INTERFACE doq_output_3d 179 MODULE PROCEDURE doq_output_3d 180 END INTERFACE doq_output_3d 181 182 INTERFACE doq_output_mask 183 MODULE PROCEDURE doq_output_mask 184 END INTERFACE doq_output_mask 185 186 INTERFACE doq_init 187 MODULE PROCEDURE doq_init 188 END INTERFACE doq_init 189 190 INTERFACE doq_prepare 191 MODULE PROCEDURE doq_prepare 192 END INTERFACE doq_prepare 193 194 ! INTERFACE doq_rrd_local 195 ! MODULE PROCEDURE doq_rrd_local 196 ! END INTERFACE doq_rrd_local 197 198 INTERFACE doq_wrd_local 199 MODULE PROCEDURE doq_wrd_local 200 END INTERFACE doq_wrd_local 142 201 143 202 144 203 CONTAINS 145 204 146 205 !! 147 206 ! Description: 148 207 !  149 !> ... 150 !! 151 SUBROUTINE diagnostic_output_quantities_init 152 153 154 USE indices, & 155 ONLY: nxl, nxr, nyn, nys, nzb, nzt 208 !> Sum up and timeaverage diagnostic output quantities as well as allocate 209 !> the array necessary for storing the average. 210 !! 211 SUBROUTINE doq_3d_data_averaging( mode, variable ) 212 213 USE control_parameters, & 214 ONLY: average_count_3d 215 216 CHARACTER (LEN=*) :: mode !< 217 CHARACTER (LEN=*) :: variable !< 218 219 INTEGER(iwp) :: i !< 220 INTEGER(iwp) :: j !< 221 INTEGER(iwp) :: k !< 222 223 IF ( mode == 'allocate' ) THEN 224 225 SELECT CASE ( TRIM( variable ) ) 226 227 CASE ( 'ti' ) 228 IF ( .NOT. ALLOCATED( ti_av ) ) THEN 229 ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 230 ENDIF 231 ti_av = 0.0_wp 232 233 CASE ( 'uu' ) 234 IF ( .NOT. ALLOCATED( uu_av ) ) THEN 235 ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 236 ENDIF 237 uu_av = 0.0_wp 238 239 CASE ( 'vv' ) 240 IF ( .NOT. ALLOCATED( vv_av ) ) THEN 241 ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 242 ENDIF 243 vv_av = 0.0_wp 244 245 CASE ( 'ww' ) 246 IF ( .NOT. ALLOCATED( ww_av ) ) THEN 247 ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 248 ENDIF 249 ww_av = 0.0_wp 250 251 CASE DEFAULT 252 CONTINUE 253 254 END SELECT 255 256 ELSEIF ( mode == 'sum' ) THEN 257 258 SELECT CASE ( TRIM( variable ) ) 259 260 CASE ( 'ti' ) 261 IF ( ALLOCATED( ti_av ) ) THEN 262 DO i = nxl, nxr 263 DO j = nys, nyn 264 DO k = nzb, nzt+1 265 ti_av(k,j,i) = ti_av(k,j,i) + ti(k,j,i) 266 ENDDO 267 ENDDO 268 ENDDO 269 ENDIF 270 271 CASE ( 'uu' ) 272 IF ( ALLOCATED( uu_av ) ) THEN 273 DO i = nxl, nxr 274 DO j = nys, nyn 275 DO k = nzb, nzt+1 276 uu_av(k,j,i) = uu_av(k,j,i) + uu(k,j,i) 277 ENDDO 278 ENDDO 279 ENDDO 280 ENDIF 281 282 CASE ( 'vv' ) 283 IF ( ALLOCATED( vv_av ) ) THEN 284 DO i = nxl, nxr 285 DO j = nys, nyn 286 DO k = nzb, nzt+1 287 vv_av(k,j,i) = vv_av(k,j,i) + vv(k,j,i) 288 ENDDO 289 ENDDO 290 ENDDO 291 ENDIF 292 293 CASE ( 'ww' ) 294 IF ( ALLOCATED( ww_av ) ) THEN 295 DO i = nxl, nxr 296 DO j = nys, nyn 297 DO k = nzb, nzt+1 298 ww_av(k,j,i) = ww_av(k,j,i) + ww(k,j,i) 299 ENDDO 300 ENDDO 301 ENDDO 302 ENDIF 303 304 CASE DEFAULT 305 CONTINUE 306 307 END SELECT 308 309 ELSEIF ( mode == 'average' ) THEN 310 311 SELECT CASE ( TRIM( variable ) ) 312 313 CASE ( 'ti' ) 314 IF ( ALLOCATED( ti_av ) ) THEN 315 DO i = nxl, nxr 316 DO j = nys, nyn 317 DO k = nzb, nzt+1 318 ti_av(k,j,i) = ti_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 319 ENDDO 320 ENDDO 321 ENDDO 322 ENDIF 323 324 CASE ( 'uu' ) 325 IF ( ALLOCATED( uu_av ) ) THEN 326 DO i = nxl, nxr 327 DO j = nys, nyn 328 DO k = nzb, nzt+1 329 uu_av(k,j,i) = uu_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 330 ENDDO 331 ENDDO 332 ENDDO 333 ENDIF 334 335 CASE ( 'vv' ) 336 IF ( ALLOCATED( vv_av ) ) THEN 337 DO i = nxl, nxr 338 DO j = nys, nyn 339 DO k = nzb, nzt+1 340 vv_av(k,j,i) = vv_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 341 ENDDO 342 ENDDO 343 ENDDO 344 ENDIF 345 346 CASE ( 'ww' ) 347 IF ( ALLOCATED( ww_av ) ) THEN 348 DO i = nxl, nxr 349 DO j = nys, nyn 350 DO k = nzb, nzt+1 351 ww_av(k,j,i) = ww_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 352 ENDDO 353 ENDDO 354 ENDDO 355 ENDIF 356 357 END SELECT 358 359 ENDIF 360 361 362 END SUBROUTINE doq_3d_data_averaging 363 364 !! 365 ! Description: 366 !  367 !> Check data output for diagnostic output 368 !! 369 SUBROUTINE doq_check_data_output( var, unit ) 156 370 157 371 IMPLICIT NONE 372 373 CHARACTER (LEN=*) :: unit !< 374 CHARACTER (LEN=*) :: var !< 375 376 SELECT CASE ( TRIM( var ) ) 377 378 CASE ( 'ti' ) 379 unit = '1/s' 380 381 CASE ( 'uu' ) 382 unit = 'm2/s2' 383 384 CASE ( 'vv' ) 385 unit = 'm2/s2' 386 387 CASE ( 'ww' ) 388 unit = 'm2/s2' 389 390 CASE DEFAULT 391 unit = 'illegal' 392 393 END SELECT 394 395 396 END SUBROUTINE doq_check_data_output 397 398 !! 399 ! 400 ! Description: 401 !  402 !> Subroutine defining appropriate grid for netcdf variables. 403 !! 404 SUBROUTINE doq_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z ) 405 406 IMPLICIT NONE 407 408 CHARACTER (LEN=*), INTENT(IN) :: variable !< 409 LOGICAL, INTENT(OUT) :: found !< 410 CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< 411 CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< 412 CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< 413 414 found = .TRUE. 415 416 SELECT CASE ( TRIM( variable ) ) 417 ! 418 ! s grid 419 CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz' ) 420 421 grid_x = 'x' 422 grid_y = 'y' 423 grid_z = 'zu' 424 ! 425 ! u grid 426 CASE ( 'uu', 'uu_xy', 'uu_xz', 'uu_yz' ) 427 428 grid_x = 'xu' 429 grid_y = 'y' 430 grid_z = 'zu' 431 ! 432 ! v grid 433 CASE ( 'vv', 'vv_xy', 'vv_xz', 'vv_yz' ) 434 435 grid_x = 'x' 436 grid_y = 'yv' 437 grid_z = 'zu' 438 ! 439 ! w grid 440 CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz' ) 441 442 grid_x = 'x' 443 grid_y = 'y' 444 grid_z = 'zw' 445 446 CASE DEFAULT 447 found = .FALSE. 448 grid_x = 'none' 449 grid_y = 'none' 450 grid_z = 'none' 451 452 END SELECT 453 454 455 END SUBROUTINE doq_define_netcdf_grid 456 457 !! 458 ! 459 ! Description: 460 !  461 !> Subroutine defining 2D output variables 462 !! 463 SUBROUTINE doq_output_2d( av, variable, found, grid, & 464 mode, local_pf, two_d, nzb_do, nzt_do, fill_value ) 465 466 467 IMPLICIT NONE 468 469 CHARACTER (LEN=*) :: grid !< 470 CHARACTER (LEN=*) :: mode !< 471 CHARACTER (LEN=*) :: variable !< 472 473 INTEGER(iwp) :: av !< value indicating averaged or nonaveraged output 474 INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w) 475 INTEGER(iwp) :: i !< grid index xdirection 476 INTEGER(iwp) :: j !< grid index ydirection 477 INTEGER(iwp) :: k !< grid index zdirection 478 INTEGER(iwp) :: nzb_do !< 479 INTEGER(iwp) :: nzt_do !< 480 481 LOGICAL :: found !< true if variable is in list 482 LOGICAL :: resorted !< true if array is resorted 483 LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections) 484 485 REAL(wp) :: fill_value !< value for the _FillValue attribute 486 487 REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< 488 REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output 489 490 flag_nr = 0 491 found = .TRUE. 492 resorted = .FALSE. 493 two_d = .FALSE. 494 495 SELECT CASE ( TRIM( variable ) ) 496 497 CASE ( 'ti_xy', 'ti_xz', 'ti_yz' ) 498 IF ( av == 0 ) THEN 499 to_be_resorted => ti 500 ELSE 501 IF ( .NOT. ALLOCATED( ti_av ) ) THEN 502 ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 503 ti_av = REAL( fill_value, KIND = wp ) 504 ENDIF 505 to_be_resorted => ti_av 506 ENDIF 507 flag_nr = 0 508 509 IF ( mode == 'xy' ) grid = 'zu' 510 511 CASE ( 'uu_xy', 'uu_xz', 'uu_yz' ) 512 IF ( av == 0 ) THEN 513 to_be_resorted => uu 514 ELSE 515 IF ( .NOT. ALLOCATED( uu_av ) ) THEN 516 ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 517 uu_av = REAL( fill_value, KIND = wp ) 518 ENDIF 519 to_be_resorted => uu_av 520 ENDIF 521 flag_nr = 1 522 523 IF ( mode == 'xy' ) grid = 'zu' 524 525 CASE ( 'vv_xy', 'vv_xz', 'vv_yz' ) 526 IF ( av == 0 ) THEN 527 to_be_resorted => vv 528 ELSE 529 IF ( .NOT. ALLOCATED( vv_av ) ) THEN 530 ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 531 vv_av = REAL( fill_value, KIND = wp ) 532 ENDIF 533 to_be_resorted => vv_av 534 ENDIF 535 flag_nr = 2 536 537 IF ( mode == 'xy' ) grid = 'zu' 538 539 CASE ( 'ww_xy', 'ww_xz', 'ww_yz' ) 540 IF ( av == 0 ) THEN 541 to_be_resorted => ww 542 ELSE 543 IF ( .NOT. ALLOCATED( ww_av ) ) THEN 544 ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 545 ww_av = REAL( fill_value, KIND = wp ) 546 ENDIF 547 to_be_resorted => ww_av 548 ENDIF 549 flag_nr = 3 550 551 IF ( mode == 'xy' ) grid = 'zw' 552 553 CASE DEFAULT 554 found = .FALSE. 555 grid = 'none' 556 557 END SELECT 558 559 IF ( found .AND. .NOT. resorted ) THEN 560 DO i = nxl, nxr 561 DO j = nys, nyn 562 DO k = nzb_do, nzt_do 563 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), & 564 REAL( fill_value, KIND = wp ), & 565 BTEST( wall_flags_0(k,j,i), flag_nr ) ) 566 ENDDO 567 ENDDO 568 ENDDO 569 ENDIF 570 571 END SUBROUTINE doq_output_2d 572 573 574 !! 575 ! 576 ! Description: 577 !  578 !> Subroutine defining 3D output variables 579 !! 580 SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do, & 581 nzt_do ) 582 583 IMPLICIT NONE 584 585 CHARACTER (LEN=*) :: variable !< 586 587 INTEGER(iwp) :: av !< index indicating averaged or instantaneous output 588 INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w) 589 INTEGER(iwp) :: i !< index variable along xdirection 590 INTEGER(iwp) :: j !< index variable along ydirection 591 INTEGER(iwp) :: k !< index variable along zdirection 592 INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0) 593 INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d) 594 595 LOGICAL :: found !< true if variable is in list 596 LOGICAL :: resorted !< true if array is resorted 597 598 REAL(wp) :: fill_value !< value for the _FillValue attribute 599 600 REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< 601 REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output 602 603 flag_nr = 0 604 found = .TRUE. 605 resorted = .FALSE. 606 607 SELECT CASE ( TRIM( variable ) ) 608 609 CASE ( 'ti' ) 610 IF ( av == 0 ) THEN 611 to_be_resorted => ti 612 ELSE 613 IF ( .NOT. ALLOCATED( ti_av ) ) THEN 614 ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 615 ti_av = REAL( fill_value, KIND = wp ) 616 ENDIF 617 to_be_resorted => ti_av 618 ENDIF 619 flag_nr = 0 620 621 CASE ( 'uu' ) 622 IF ( av == 0 ) THEN 623 to_be_resorted => uu 624 ELSE 625 IF ( .NOT. ALLOCATED( uu_av ) ) THEN 626 ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 627 uu_av = REAL( fill_value, KIND = wp ) 628 ENDIF 629 to_be_resorted => uu_av 630 ENDIF 631 flag_nr = 1 632 633 CASE ( 'vv' ) 634 IF ( av == 0 ) THEN 635 to_be_resorted => vv 636 ELSE 637 IF ( .NOT. ALLOCATED( vv_av ) ) THEN 638 ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 639 vv_av = REAL( fill_value, KIND = wp ) 640 ENDIF 641 to_be_resorted => vv_av 642 ENDIF 643 flag_nr = 2 644 645 CASE ( 'ww' ) 646 IF ( av == 0 ) THEN 647 to_be_resorted => ww 648 ELSE 649 IF ( .NOT. ALLOCATED( ww_av ) ) THEN 650 ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 651 ww_av = REAL( fill_value, KIND = wp ) 652 ENDIF 653 to_be_resorted => ww_av 654 ENDIF 655 flag_nr = 3 656 657 CASE DEFAULT 658 found = .FALSE. 659 660 END SELECT 661 662 IF ( found .AND. .NOT. resorted ) THEN 663 DO i = nxl, nxr 664 DO j = nys, nyn 665 DO k = nzb_do, nzt_do 666 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), & 667 REAL( fill_value, KIND = wp ), & 668 BTEST( wall_flags_0(k,j,i), flag_nr ) ) 669 ENDDO 670 ENDDO 671 ENDDO 672 ENDIF 673 674 END SUBROUTINE doq_output_3d 675 676 ! Description: 677 !  678 !> Resorts the userdefined output quantity with indices (k,j,i) to a 679 !> temporary array with indices (i,j,k) for masked data output. 680 !! 681 SUBROUTINE doq_output_mask( av, variable, found, local_pf ) 682 683 USE control_parameters 684 685 USE indices 686 687 USE surface_mod, & 688 ONLY: get_topography_top_index_ji 689 690 IMPLICIT NONE 691 692 CHARACTER (LEN=*) :: variable !< 693 CHARACTER (LEN=5) :: grid !< flag to distinquish between staggered grids 694 695 INTEGER(iwp) :: av !< index indicating averaged or instantaneous output 696 INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w) 697 INTEGER(iwp) :: i !< index variable along xdirection 698 INTEGER(iwp) :: j !< index variable along ydirection 699 INTEGER(iwp) :: k !< index variable along zdirection 700 INTEGER(iwp) :: topo_top_ind !< k index of highest horizontal surface 701 702 LOGICAL :: found !< true if variable is in list 703 LOGICAL :: resorted !< true if array is resorted 704 705 REAL(wp), & 706 DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: & 707 local_pf !< 708 REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output 709 710 flag_nr = 0 711 found = .TRUE. 712 resorted = .FALSE. 713 grid = 's' 714 715 SELECT CASE ( TRIM( variable ) ) 716 717 CASE ( 'ti' ) 718 IF ( av == 0 ) THEN 719 to_be_resorted => ti 720 ELSE 721 to_be_resorted => ti_av 722 ENDIF 723 grid = 's' 724 flag_nr = 0 725 726 CASE ( 'uu' ) 727 IF ( av == 0 ) THEN 728 to_be_resorted => uu 729 ELSE 730 to_be_resorted => uu_av 731 ENDIF 732 grid = 'u' 733 flag_nr = 1 734 735 CASE ( 'vv' ) 736 IF ( av == 0 ) THEN 737 to_be_resorted => vv 738 ELSE 739 to_be_resorted => vv_av 740 ENDIF 741 grid = 'v' 742 flag_nr = 2 743 744 CASE ( 'ww' ) 745 IF ( av == 0 ) THEN 746 to_be_resorted => ww 747 ELSE 748 to_be_resorted => ww_av 749 ENDIF 750 grid = 'w' 751 flag_nr = 3 752 753 754 CASE DEFAULT 755 found = .FALSE. 756 757 END SELECT 758 759 IF ( found .AND. .NOT. resorted ) THEN 760 IF ( .NOT. mask_surface(mid) ) THEN 761 ! 762 ! Default masked output 763 DO i = 1, mask_size_l(mid,1) 764 DO j = 1, mask_size_l(mid,2) 765 DO k = 1, mask_size_l(mid,3) 766 local_pf(i,j,k) = to_be_resorted(mask_k(mid,k), & 767 mask_j(mid,j), & 768 mask_i(mid,i)) 769 ENDDO 770 ENDDO 771 ENDDO 772 773 ELSE 774 ! 775 ! Terrainfollowing masked output 776 DO i = 1, mask_size_l(mid,1) 777 DO j = 1, mask_size_l(mid,2) 778 ! 779 ! Get k index of highest horizontal surface 780 topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), & 781 mask_i(mid,i), & 782 grid ) 783 ! 784 ! Save output array 785 DO k = 1, mask_size_l(mid,3) 786 local_pf(i,j,k) = to_be_resorted( & 787 MIN( topo_top_ind+mask_k(mid,k), & 788 nzt+1 ), & 789 mask_j(mid,j), & 790 mask_i(mid,i) ) 791 ENDDO 792 ENDDO 793 ENDDO 794 795 ENDIF 796 ENDIF 797 798 END SUBROUTINE doq_output_mask 799 800 !! 801 ! Description: 802 !  803 !> Allocate required arrays 804 !! 805 SUBROUTINE doq_init 806 807 IMPLICIT NONE 808 809 INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities 158 810 ! 159 811 ! Next line is to avoid compiler warnings about unused variables … … 161 813 162 814 initialized_diagnostic_output_quantities = .FALSE. 163 164 IF ( .NOT. ALLOCATED( ti ) ) THEN 165 ALLOCATE( ti(nzb:nzt+1,nys:nyn,nxl:nxr) ) 166 ti = 0.0_wp 167 ENDIF 815 816 ivar = 1 817 818 DO WHILE ( ivar <= SIZE( do_all ) ) 819 820 SELECT CASE ( TRIM( do_all(ivar) ) ) 821 ! 822 ! Allocate array for 'turbulence intensity' 823 CASE ( 'ti' ) 824 IF ( .NOT. ALLOCATED( ti ) ) THEN 825 ALLOCATE( ti(nzb:nzt+1,nys:nyn,nxl:nxr) ) 826 ti = 0.0_wp 827 ENDIF 828 ! 829 ! Allocate array for uu 830 CASE ( 'uu' ) 831 IF ( .NOT. ALLOCATED( uu ) ) THEN 832 ALLOCATE( uu(nzb:nzt+1,nys:nyn,nxl:nxr) ) 833 uu = 0.0_wp 834 ENDIF 835 836 ! 837 ! Allocate array for vv 838 CASE ( 'vv' ) 839 IF ( .NOT. ALLOCATED( vv ) ) THEN 840 ALLOCATE( vv(nzb:nzt+1,nys:nyn,nxl:nxr) ) 841 vv = 0.0_wp 842 ENDIF 843 844 ! 845 ! Allocate array for ww 846 CASE ( 'ww' ) 847 IF ( .NOT. ALLOCATED( ww ) ) THEN 848 ALLOCATE( ww(nzb:nzt+1,nys:nyn,nxl:nxr) ) 849 ww = 0.0_wp 850 ENDIF 851 852 END SELECT 853 854 ivar = ivar + 1 855 ENDDO 168 856 169 857 initialized_diagnostic_output_quantities = .TRUE. 170 171 END SUBROUTINE d iagnostic_output_quantities_init858 859 END SUBROUTINE doq_init 172 860 173 861 … … 175 863 ! Description: 176 864 !  177 !> ...865 !> Calculate diagnostic quantities 178 866 !! 179 SUBROUTINE diagnostic_output_quantities_calculate 180 181 182 USE arrays_3d, & 183 ONLY: ddzu, u, v, w 184 185 USE grid_variables, & 186 ONLY: ddx, ddy 187 188 USE indices, & 189 ONLY: nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0 867 SUBROUTINE doq_calculate 190 868 191 869 IMPLICIT NONE 192 870 193 871 INTEGER(iwp) :: i, j, k !< grid loop index in x, y, zdirection 194 INTEGER(iwp) :: ivar _all!< loop index over all 2d/3d/mask output quantities872 INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities 195 873 196 874 … … 198 876 ! 199 877 ! Preparatory steps and initialization of output arrays 200 IF ( .NOT. prepared_diagnostic_output_quantities ) CALL d iagnostic_output_quantities_prepare201 IF ( .NOT. initialized_diagnostic_output_quantities ) CALL d iagnostic_output_quantities_init202 ! 203 ! Save timestep number to check in time_integration if d iagnostic_output_quantities_calculate878 IF ( .NOT. prepared_diagnostic_output_quantities ) CALL doq_prepare 879 IF ( .NOT. initialized_diagnostic_output_quantities ) CALL doq_init 880 ! 881 ! Save timestep number to check in time_integration if doq_calculate 204 882 ! has been called already, since the CALL occurs at two locations, but the calculations need to be 205 883 ! done only once per timestep. 206 884 timestep_number_at_prev_calc = current_timestep_number 207 885 208 209 ivar_all = 1 210 211 DO WHILE ( do_all(ivar_all)(1:1) /= ' ' ) 212 213 SELECT CASE ( TRIM( do_all(ivar_all) ) ) 886 ivar = 1 887 888 DO WHILE ( ivar <= SIZE( do_all ) ) 889 890 SELECT CASE ( TRIM( do_all(ivar) ) ) 214 891 ! 215 892 ! Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point … … 218 895 DO j = nys, nyn 219 896 DO k = nzb+1, nzt 220 221 ti(k,j,i) = 0.25_wp * SQRT( & 222 ( ( w(k,j+1,i) + w(k1,j+1,i) & 223  w(k,j1,i)  w(k1,j1,i) ) * ddy & 224  ( v(k+1,j,i) + v(k+1,j+1,i) & 225  v(k1,j,i)  v(k1,j+1,i) ) * ddzu(k) )**2 & 226 + ( ( u(k+1,j,i) + u(k+1,j,i+1) & 227  u(k1,j,i)  u(k1,j,i+1) ) * ddzu(k) & 228  ( w(k,j,i+1) + w(k1,j,i+1) & 229  w(k,j,i1)  w(k1,j,i1) ) * ddx )**2 & 230 + ( ( v(k,j,i+1) + v(k,j+1,i+1) & 231  v(k,j,i1)  v(k,j+1,i1) ) * ddx & 232  ( u(k,j+1,i) + u(k,j+1,i+1) & 233  u(k,j1,i)  u(k,j1,i+1) ) * ddy )**2 ) & 234 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0) ) 897 ti(k,j,i) = 0.25_wp * SQRT( & 898 ( ( w(k,j+1,i) + w(k1,j+1,i) & 899  w(k,j1,i)  w(k1,j1,i) ) * ddy & 900  ( v(k+1,j,i) + v(k+1,j+1,i) & 901  v(k1,j,i)  v(k1,j+1,i) ) * ddzu(k) )**2 & 902 + ( ( u(k+1,j,i) + u(k+1,j,i+1) & 903  u(k1,j,i)  u(k1,j,i+1) ) * ddzu(k) & 904  ( w(k,j,i+1) + w(k1,j,i+1) & 905  w(k,j,i1)  w(k1,j,i1) ) * ddx )**2 & 906 + ( ( v(k,j,i+1) + v(k,j+1,i+1) & 907  v(k,j,i1)  v(k,j+1,i1) ) * ddx & 908  ( u(k,j+1,i) + u(k,j+1,i+1) & 909  u(k,j1,i)  u(k,j1,i+1) ) * ddy )**2 ) & 910 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0) ) 235 911 ENDDO 236 912 ENDDO 913 ENDDO 914 ! 915 ! uu 916 CASE ( 'uu' ) 917 DO i = nxl, nxr 918 DO j = nys, nyn 919 DO k = nzb+1, nzt 920 uu(k,j,i) = u(k,j,i) * u(k,j,i) & 921 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 1) ) 922 ENDDO 923 ENDDO 237 924 ENDDO 238 239 END SELECT 240 241 ivar_all = ivar_all + 1 925 ! 926 ! vv 927 CASE ( 'vv' ) 928 DO i = nxl, nxr 929 DO j = nys, nyn 930 DO k = nzb+1, nzt 931 vv(k,j,i) = v(k,j,i) * v(k,j,i) & 932 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 2) ) 933 ENDDO 934 ENDDO 935 ENDDO 936 ! 937 ! ww 938 CASE ( 'ww' ) 939 DO i = nxl, nxr 940 DO j = nys, nyn 941 DO k = nzb+1, nzt1 942 ww(k,j,i) = w(k,j,i) * w(k,j,i) & 943 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 3) ) 944 ENDDO 945 ENDDO 946 ENDDO 947 948 END SELECT 949 950 ivar = ivar + 1 242 951 ENDDO 243 952 244 953 ! CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' ) 245 954 246 END SUBROUTINE d iagnostic_output_quantities_calculate955 END SUBROUTINE doq_calculate 247 956 248 957 … … 252 961 !> ... 253 962 !! 254 SUBROUTINE d iagnostic_output_quantities_prepare255 256 257 USE control_parameters, 963 SUBROUTINE doq_prepare 964 965 966 USE control_parameters, & 258 967 ONLY: do2d, do3d, domask, masks, mid 259 968 … … 268 977 INTEGER(iwp) :: l !< index for cutting string 269 978 270 271 979 prepared_diagnostic_output_quantities = .FALSE. 272 980 … … 280 988 do2d_var(av,ivar)(1:l3) = do2d(av,ivar)(1:l3) 281 989 ! 282 ! Gather 2d output quantity names, check for double occurrence of output quantity 990 ! Gather 2d output quantity names. 991 ! Check for double occurrence of output quantity, e.g. by _xy, 992 ! _yz, _xz. 283 993 DO WHILE ( do2d_var(av,ivar)(1:1) /= ' ' ) 284 994 IF ( .NOT. ANY( do_all == do2d_var(av,ivar) ) ) THEN … … 293 1003 ivar = 1 294 1004 ! 295 ! Gather 3d output quantity names , check for double occurrence of output quantity1005 ! Gather 3d output quantity names 296 1006 DO WHILE ( do3d(av,ivar)(1:1) /= ' ' ) 297 1007 do_all(ivar_all) = do3d(av,ivar) 298 299 1008 ivar = ivar + 1 300 1009 ivar_all = ivar_all + 1 … … 303 1012 ivar = 1 304 1013 ! 305 ! Gather masked output quantity names, check for double occurrence of output quantity 1014 ! Gather masked output quantity names. Also check for double output 1015 ! e.g. by different masks. 306 1016 DO mid = 1, masks 307 1017 DO WHILE ( domask(mid,av,ivar)(1:1) /= ' ' ) 308 do_all(ivar_all) = domask(mid,av,ivar) 1018 IF ( .NOT. ANY( do_all == domask(mid,av,ivar) ) ) THEN 1019 do_all(ivar_all) = do2d_var(av,ivar) 1020 ENDIF 309 1021 310 1022 ivar = ivar + 1 … … 318 1030 prepared_diagnostic_output_quantities = .TRUE. 319 1031 320 END SUBROUTINE diagnostic_output_quantities_prepare 1032 END SUBROUTINE doq_prepare 1033 1034 !! 1035 ! Description: 1036 !  1037 !> Subroutine reads local (subdomain) restart data 1038 !> Note: With the current structure reading of nonstandard array is not 1039 !> possible 1040 !! 1041 ! SUBROUTINE doq_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, & 1042 ! nxr_on_file, nynf, nync, nyn_on_file, nysf, & 1043 ! nysc, nys_on_file, tmp_3d_non_standard, found ) 1044 ! 1045 ! 1046 ! USE control_parameters 1047 ! 1048 ! USE indices 1049 ! 1050 ! USE kinds 1051 ! 1052 ! USE pegrid 1053 ! 1054 ! 1055 ! IMPLICIT NONE 1056 ! 1057 ! INTEGER(iwp) :: k !< 1058 ! INTEGER(iwp) :: nxlc !< 1059 ! INTEGER(iwp) :: nxlf !< 1060 ! INTEGER(iwp) :: nxl_on_file !< 1061 ! INTEGER(iwp) :: nxrc !< 1062 ! INTEGER(iwp) :: nxrf !< 1063 ! INTEGER(iwp) :: nxr_on_file !< 1064 ! INTEGER(iwp) :: nync !< 1065 ! INTEGER(iwp) :: nynf !< 1066 ! INTEGER(iwp) :: nyn_on_file !< 1067 ! INTEGER(iwp) :: nysc !< 1068 ! INTEGER(iwp) :: nysf !< 1069 ! INTEGER(iwp) :: nys_on_file !< 1070 ! 1071 ! LOGICAL, INTENT(OUT) :: found 1072 ! 1073 ! REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmp_3d_non_standard !< temporary array for storing 3D data with non standard dimensions 1074 ! ! 1075 ! ! If temporary nonstandard array for reading is already allocated, 1076 ! ! deallocate it. 1077 ! IF ( ALLOCATED( tmp_3d_non_standard ) ) DEALLOCATE( tmp_3d_non_standard ) 1078 ! 1079 ! found = .TRUE. 1080 ! 1081 ! SELECT CASE ( restart_string(1:length) ) 1082 ! 1083 ! CASE ( 'ti_av' ) 1084 ! IF ( .NOT. ALLOCATED( ti_av ) ) THEN 1085 ! ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 1086 ! ENDIF 1087 ! IF ( k == 1 ) THEN 1088 ! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, & 1089 ! nxl_on_file:nxr_on_file) ) 1090 ! READ ( 13 ) tmp_3d_non_standard 1091 ! ENDIF 1092 ! ti_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf) 1093 ! 1094 ! CASE ( 'uu_av' ) 1095 ! IF ( .NOT. ALLOCATED( uu_av ) ) THEN 1096 ! ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 1097 ! ENDIF 1098 ! IF ( k == 1 ) THEN 1099 ! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, & 1100 ! nxl_on_file:nxr_on_file) ) 1101 ! READ ( 13 ) tmp_3d_non_standard 1102 ! ENDIF 1103 ! uu_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf) 1104 ! 1105 ! CASE ( 'vv_av' ) 1106 ! IF ( .NOT. ALLOCATED( vv_av ) ) THEN 1107 ! ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 1108 ! ENDIF 1109 ! IF ( k == 1 ) THEN 1110 ! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, & 1111 ! nxl_on_file:nxr_on_file) ) 1112 ! READ ( 13 ) tmp_3d_non_standard 1113 ! ENDIF 1114 ! vv_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf) 1115 ! 1116 ! CASE ( 'ww_av' ) 1117 ! IF ( .NOT. ALLOCATED( ww_av ) ) THEN 1118 ! ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 1119 ! ENDIF 1120 ! IF ( k == 1 ) THEN 1121 ! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, & 1122 ! nxl_on_file:nxr_on_file) ) 1123 ! READ ( 13 ) tmp_3d_non_standard 1124 ! ENDIF 1125 ! ww_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf) 1126 ! 1127 ! 1128 ! CASE DEFAULT 1129 ! 1130 ! found = .FALSE. 1131 ! 1132 ! END SELECT 1133 ! 1134 ! END SUBROUTINE doq_rrd_local 1135 1136 !! 1137 ! Description: 1138 !  1139 !> Subroutine writes local (subdomain) restart data 1140 !! 1141 SUBROUTINE doq_wrd_local 1142 1143 1144 IMPLICIT NONE 1145 1146 IF ( ALLOCATED( ti_av ) ) THEN 1147 CALL wrd_write_string( 'ti_av' ) 1148 WRITE ( 14 ) ti_av 1149 ENDIF 1150 1151 IF ( ALLOCATED( uu_av ) ) THEN 1152 CALL wrd_write_string( 'uu_av' ) 1153 WRITE ( 14 ) uu_av 1154 ENDIF 1155 1156 IF ( ALLOCATED( vv_av ) ) THEN 1157 CALL wrd_write_string( 'vv_av' ) 1158 WRITE ( 14 ) vv_av 1159 ENDIF 1160 1161 IF ( ALLOCATED( ww_av ) ) THEN 1162 CALL wrd_write_string( 'ww_av' ) 1163 WRITE ( 14 ) ww_av 1164 ENDIF 1165 1166 1167 END SUBROUTINE doq_wrd_local 1168 1169 321 1170 322 1171 END MODULE diagnostic_output_quantities_mod
Note: See TracChangeset
for help on using the changeset viewer.