Changeset 4583 for palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
- Timestamp:
- Jun 29, 2020 12:36:47 PM (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
r4560 r4583 1 1 !> @file diagnostic_output_quantities_mod.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 4 ! 5 ! PALM is free software: you can redistribute it and/or modify it under the 6 ! terms of the GNU General Public License as published by the Free Software 7 ! Foundation, either version 3 of the License, or (at your option) any later 8 ! version. 9 ! 10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY 11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 12 ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. 13 ! 14 ! You should have received a copy of the GNU General Public License along with 15 ! PALM. If not, see <http://www.gnu.org/licenses/>. 5 ! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General 6 ! Public License as published by the Free Software Foundation, either version 3 of the License, or 7 ! (at your option) any later version. 8 ! 9 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the 10 ! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 11 ! Public License for more details. 12 ! 13 ! You should have received a copy of the GNU General Public License along with PALM. If not, see 14 ! <http://www.gnu.org/licenses/>. 16 15 ! 17 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 18 !------------------------------------------------------------------------------ !17 !--------------------------------------------------------------------------------------------------! 19 18 ! 20 19 ! Current revisions: … … 25 24 ! ----------------- 26 25 ! $Id$ 26 ! file re-formatted to follow the PALM coding standard 27 ! 28 ! 4560 2020-06-11 12:19:47Z suehring 27 29 ! - Bugfix in calculation of vertical momentum and scalar fluxes 28 30 ! - remove averaged output variables from PUBLIC list 29 ! 31 ! 30 32 ! 4535 2020-05-15 12:07:23Z raasch 31 33 ! bugfix for restart data format query 32 ! 34 ! 33 35 ! 4518 2020-05-04 15:44:28Z suehring 34 ! * Define arrays over ghost points in order to allow for standard mpi-io 35 ! treatment. By this modularization of restart-data input is possible with 36 ! the module interface. 36 ! * Define arrays over ghost points in order to allow for standard mpi-io treatment. By this 37 ! modularization of restart-data input is possible with the module interface. 37 38 ! * Bugfix: add missing restart input of wtheta_av, wq_av, wu_av, and wv_av. 38 ! 39 ! 39 40 ! 4517 2020-05-03 14:29:30Z raasch 40 41 ! use statement for exchange horiz added, 41 42 ! bugfix for call of exchange horiz 2d 42 ! 43 ! 43 44 ! 4431 2020-02-27 23:23:01Z gronemeier 44 45 ! added wspeed and wdir output; bugfix: set fill_value in case of masked output 45 46 ! 46 47 ! 4360 2020-01-07 11:25:50Z suehring 47 ! added output of wu, wv, wtheta and wq to enable covariance calculation 48 ! according to temporal ECmethod48 ! added output of wu, wv, wtheta and wq to enable covariance calculation according to temporal EC 49 ! method 49 50 ! 50 51 ! 4346 2019-12-18 11:55:56Z motisi 51 ! Introduction of wall_flags_total_0, which currently sets bits based on static 52 ! topographyinformation used in wall_flags_static_052 ! Introduction of wall_flags_total_0, which currently sets bits based on static topography 53 ! information used in wall_flags_static_0 53 54 ! 54 55 ! 4331 2019-12-10 18:25:02Z suehring … … 63 64 ! 64 65 ! 4167 2019-08-16 11:01:48Z suehring 65 ! Changed behaviour of masked output over surface to follow terrain and ignore 66 ! buildings(J.Resler, T.Gronemeier)66 ! Changed behaviour of masked output over surface to follow terrain and ignore buildings 67 ! (J.Resler, T.Gronemeier) 67 68 ! 68 69 ! 4157 2019-08-14 09:19:12Z suehring 69 ! Initialization restructured, in order to work also when data output during 70 ! spin-up is enabled. 70 ! Initialization restructured, in order to work also when data output during spin-up is enabled. 71 71 ! 72 72 ! 4132 2019-08-02 12:34:17Z suehring … … 74 74 ! 75 75 ! 4069 2019-07-01 14:05:51Z Giersch 76 ! Masked output running index mid has been introduced as a local variable to 77 ! avoid runtime error(Loop variable has been modified) in time_integration76 ! Masked output running index mid has been introduced as a local variable to avoid runtime error 77 ! (Loop variable has been modified) in time_integration 78 78 ! 79 79 ! 4039 2019-06-18 10:32:41Z suehring 80 ! - Add output of uu, vv, ww to enable variance calculation according temporal 81 ! EC method 80 ! - Add output of uu, vv, ww to enable variance calculation according temporal EC method 82 81 ! - Allocate arrays only when they are required 83 82 ! - Formatting adjustment … … 89 88 ! 90 89 ! 3995 2019-05-22 18:59:54Z suehring 91 ! Avoid compiler warnings about unused variable and fix string operation which 92 ! is not allowed withPGI compiler90 ! Avoid compiler warnings about unused variable and fix string operation which is not allowed with 91 ! PGI compiler 93 92 ! 94 93 ! 3994 2019-05-22 18:08:09Z suehring … … 103 102 ! ------------ 104 103 !> ... 105 !------------------------------------------------------------------------------ !104 !--------------------------------------------------------------------------------------------------! 106 105 MODULE diagnostic_output_quantities_mod 107 106 108 USE arrays_3d, &109 ONLY: ddzu, &110 pt, &111 q, &112 u, &113 v, &114 w, &115 zu, &107 USE arrays_3d, & 108 ONLY: ddzu, & 109 pt, & 110 q, & 111 u, & 112 v, & 113 w, & 114 zu, & 116 115 zw 117 116 118 USE basic_constants_and_equations_mod, &117 USE basic_constants_and_equations_mod, & 119 118 ONLY: kappa, pi 120 119 121 USE control_parameters, &122 ONLY: current_timestep_number, &123 data_output, &124 message_string, &125 restart_data_format_output, &120 USE control_parameters, & 121 ONLY: current_timestep_number, & 122 data_output, & 123 message_string, & 124 restart_data_format_output, & 126 125 varnamelength 127 126 ! 128 ! USE cpulog, &127 ! USE cpulog, & 129 128 ! ONLY: cpu_log, log_point 130 129 131 USE exchange_horiz_mod, &130 USE exchange_horiz_mod, & 132 131 ONLY: exchange_horiz_2d 133 132 134 USE grid_variables, &133 USE grid_variables, & 135 134 ONLY: ddx, ddy 136 135 137 USE indices, &138 ONLY: nbgp, &139 nxl, &140 nxlg, &141 nxr, &142 nxrg, &143 nyn, &144 nyng, &145 nys, &146 nysg, &147 nzb, &148 nzt, &136 USE indices, & 137 ONLY: nbgp, & 138 nxl, & 139 nxlg, & 140 nxr, & 141 nxrg, & 142 nyn, & 143 nyng, & 144 nys, & 145 nysg, & 146 nzb, & 147 nzt, & 149 148 wall_flags_total_0 150 149 … … 156 155 wrd_mpi_io 157 156 158 USE surface_mod, &159 ONLY: surf_def_h, &160 surf_lsm_h, &161 surf_type, &157 USE surface_mod, & 158 ONLY: surf_def_h, & 159 surf_lsm_h, & 160 surf_type, & 162 161 surf_usm_h 163 162 … … 209 208 ! 210 209 !-- Public variables 211 PUBLIC do_all, &212 initialized_diagnostic_output_quantities, &213 prepared_diagnostic_output_quantities, &210 PUBLIC do_all, & 211 initialized_diagnostic_output_quantities, & 212 prepared_diagnostic_output_quantities, & 214 213 timestep_number_at_prev_calc 215 214 ! 216 215 !-- Public routines 217 PUBLIC doq_3d_data_averaging, &218 doq_calculate, &219 doq_check_data_output, &220 doq_define_netcdf_grid, &221 doq_init, &222 doq_output_2d, &223 doq_output_3d, &224 doq_output_mask, &225 doq_rrd_local, &216 PUBLIC doq_3d_data_averaging, & 217 doq_calculate, & 218 doq_check_data_output, & 219 doq_define_netcdf_grid, & 220 doq_init, & 221 doq_output_2d, & 222 doq_output_3d, & 223 doq_output_mask, & 224 doq_rrd_local, & 226 225 doq_wrd_local 227 226 … … 275 274 CONTAINS 276 275 277 !------------------------------------------------------------------------------ !276 !--------------------------------------------------------------------------------------------------! 278 277 ! Description: 279 278 ! ------------ 280 !> Sum up and time-average diagnostic output quantities as well as allocate 281 !> the array necessary forstoring the average.282 !------------------------------------------------------------------------------ !279 !> Sum up and time-average diagnostic output quantities as well as allocate the array necessary for 280 !> storing the average. 281 !--------------------------------------------------------------------------------------------------! 283 282 SUBROUTINE doq_3d_data_averaging( mode, variable ) 284 283 285 USE control_parameters, &284 USE control_parameters, & 286 285 ONLY: average_count_3d 287 286 … … 293 292 INTEGER(iwp) :: k !< 294 293 294 295 295 IF ( mode == 'allocate' ) THEN 296 296 … … 383 383 384 384 CASE ( 'ti' ) 385 IF ( ALLOCATED( ti_av ) ) THEN385 IF ( ALLOCATED( ti_av ) ) THEN 386 386 DO i = nxl, nxr 387 387 DO j = nys, nyn … … 394 394 395 395 CASE ( 'uu' ) 396 IF ( ALLOCATED( uu_av ) ) THEN396 IF ( ALLOCATED( uu_av ) ) THEN 397 397 DO i = nxl, nxr 398 398 DO j = nys, nyn … … 405 405 406 406 CASE ( 'vv' ) 407 IF ( ALLOCATED( vv_av ) ) THEN407 IF ( ALLOCATED( vv_av ) ) THEN 408 408 DO i = nxl, nxr 409 409 DO j = nys, nyn … … 416 416 417 417 CASE ( 'ww' ) 418 IF ( ALLOCATED( ww_av ) ) THEN418 IF ( ALLOCATED( ww_av ) ) THEN 419 419 DO i = nxl, nxr 420 420 DO j = nys, nyn … … 427 427 428 428 CASE ( 'wu' ) 429 IF ( ALLOCATED( wu_av ) ) THEN429 IF ( ALLOCATED( wu_av ) ) THEN 430 430 DO i = nxl, nxr 431 431 DO j = nys, nyn … … 438 438 439 439 CASE ( 'wv' ) 440 IF ( ALLOCATED( wv_av ) ) THEN440 IF ( ALLOCATED( wv_av ) ) THEN 441 441 DO i = nxl, nxr 442 442 DO j = nys, nyn … … 449 449 450 450 CASE ( 'wtheta' ) 451 IF ( ALLOCATED( wtheta_av ) ) THEN451 IF ( ALLOCATED( wtheta_av ) ) THEN 452 452 DO i = nxl, nxr 453 453 DO j = nys, nyn … … 460 460 461 461 CASE ( 'wq' ) 462 IF ( ALLOCATED( wq_av ) ) THEN462 IF ( ALLOCATED( wq_av ) ) THEN 463 463 DO i = nxl, nxr 464 464 DO j = nys, nyn … … 471 471 472 472 CASE ( 'theta_2m*' ) 473 IF ( ALLOCATED( pt_2m_av ) ) THEN473 IF ( ALLOCATED( pt_2m_av ) ) THEN 474 474 DO i = nxl, nxr 475 475 DO j = nys, nyn … … 480 480 481 481 CASE ( 'wspeed_10m*' ) 482 IF ( ALLOCATED( uv_10m_av ) ) THEN482 IF ( ALLOCATED( uv_10m_av ) ) THEN 483 483 DO i = nxl, nxr 484 484 DO j = nys, nyn … … 489 489 490 490 CASE ( 'wspeed' ) 491 IF ( ALLOCATED( wspeed_av ) ) THEN491 IF ( ALLOCATED( wspeed_av ) ) THEN 492 492 DO i = nxl, nxr 493 493 DO j = nys, nyn … … 500 500 501 501 CASE ( 'wdir' ) 502 IF ( ALLOCATED( u_center_av ) .AND. ALLOCATED( v_center_av ) ) THEN502 IF ( ALLOCATED( u_center_av ) .AND. ALLOCATED( v_center_av ) ) THEN 503 503 DO i = nxl, nxr 504 504 DO j = nys, nyn … … 521 521 522 522 CASE ( 'ti' ) 523 IF ( ALLOCATED( ti_av ) ) THEN523 IF ( ALLOCATED( ti_av ) ) THEN 524 524 DO i = nxl, nxr 525 525 DO j = nys, nyn … … 532 532 533 533 CASE ( 'uu' ) 534 IF ( ALLOCATED( uu_av ) ) THEN534 IF ( ALLOCATED( uu_av ) ) THEN 535 535 DO i = nxl, nxr 536 536 DO j = nys, nyn … … 543 543 544 544 CASE ( 'vv' ) 545 IF ( ALLOCATED( vv_av ) ) THEN545 IF ( ALLOCATED( vv_av ) ) THEN 546 546 DO i = nxl, nxr 547 547 DO j = nys, nyn … … 554 554 555 555 CASE ( 'ww' ) 556 IF ( ALLOCATED( ww_av ) ) THEN556 IF ( ALLOCATED( ww_av ) ) THEN 557 557 DO i = nxl, nxr 558 558 DO j = nys, nyn … … 565 565 566 566 CASE ( 'wu' ) 567 IF ( ALLOCATED( wu_av ) ) THEN567 IF ( ALLOCATED( wu_av ) ) THEN 568 568 DO i = nxl, nxr 569 569 DO j = nys, nyn … … 576 576 577 577 CASE ( 'wv' ) 578 IF ( ALLOCATED( wv_av ) ) THEN578 IF ( ALLOCATED( wv_av ) ) THEN 579 579 DO i = nxl, nxr 580 580 DO j = nys, nyn … … 587 587 588 588 CASE ( 'wtheta' ) 589 IF ( ALLOCATED( wtheta_av ) ) THEN589 IF ( ALLOCATED( wtheta_av ) ) THEN 590 590 DO i = nxl, nxr 591 591 DO j = nys, nyn … … 598 598 599 599 CASE ( 'wq' ) 600 IF ( ALLOCATED( wq_av ) ) THEN600 IF ( ALLOCATED( wq_av ) ) THEN 601 601 DO i = nxl, nxr 602 602 DO j = nys, nyn … … 609 609 610 610 CASE ( 'theta_2m*' ) 611 IF ( ALLOCATED( pt_2m_av ) ) THEN611 IF ( ALLOCATED( pt_2m_av ) ) THEN 612 612 DO i = nxlg, nxrg 613 613 DO j = nysg, nyng … … 619 619 620 620 CASE ( 'wspeed_10m*' ) 621 IF ( ALLOCATED( uv_10m_av ) ) THEN621 IF ( ALLOCATED( uv_10m_av ) ) THEN 622 622 DO i = nxlg, nxrg 623 623 DO j = nysg, nyng … … 629 629 630 630 CASE ( 'wspeed' ) 631 IF ( ALLOCATED( wspeed_av ) ) THEN631 IF ( ALLOCATED( wspeed_av ) ) THEN 632 632 DO i = nxl, nxr 633 633 DO j = nys, nyn … … 640 640 641 641 CASE ( 'wdir' ) 642 IF ( ALLOCATED( u_center_av ) .AND. ALLOCATED( v_center_av ) ) THEN642 IF ( ALLOCATED( u_center_av ) .AND. ALLOCATED( v_center_av ) ) THEN 643 643 644 644 IF ( .NOT. ALLOCATED( wdir_av ) ) THEN … … 652 652 u_center_av(k,j,i) = u_center_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 653 653 v_center_av(k,j,i) = v_center_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 654 wdir_av(k,j,i) = ATAN2( u_center_av(k,j,i), v_center_av(k,j,i) ) &655 / pi * 180.0_wp + 180.0_wp654 wdir_av(k,j,i) = ATAN2( u_center_av(k,j,i), v_center_av(k,j,i) ) & 655 / pi * 180.0_wp + 180.0_wp 656 656 ENDDO 657 657 ENDDO … … 666 666 END SUBROUTINE doq_3d_data_averaging 667 667 668 !------------------------------------------------------------------------------ !668 !--------------------------------------------------------------------------------------------------! 669 669 ! Description: 670 670 ! ------------ 671 671 !> Check data output for diagnostic output 672 !------------------------------------------------------------------------------ !672 !--------------------------------------------------------------------------------------------------! 673 673 SUBROUTINE doq_check_data_output( var, unit, i, ilen, k ) 674 674 … … 682 682 INTEGER(iwp), OPTIONAL, INTENT(IN) :: k !< Output is xy mode? 0 = no, 1 = yes 683 683 684 684 685 SELECT CASE ( TRIM( var ) ) 685 686 … … 719 720 !-- Check if output quantity is _xy only. 720 721 IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN 721 message_string = 'illegal value for data_output: "' // &722 TRIM( var ) // '" & only 2d-horizontal ' // &722 message_string = 'illegal value for data_output: "' // & 723 TRIM( var ) // '" & only 2d-horizontal ' // & 723 724 'cross sections are allowed for this value' 724 725 CALL message( 'diagnostic_output', 'PA0111', 1, 2, 0, 6, 0 ) … … 736 737 END SUBROUTINE doq_check_data_output 737 738 738 !------------------------------------------------------------------------------ !739 !--------------------------------------------------------------------------------------------------! 739 740 ! 740 741 ! Description: 741 742 ! ------------ 742 743 !> Subroutine defining appropriate grid for netcdf variables. 743 !------------------------------------------------------------------------------ !744 !--------------------------------------------------------------------------------------------------! 744 745 SUBROUTINE doq_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z ) 745 746 … … 747 748 748 749 CHARACTER (LEN=*), INTENT(IN) :: variable !< 749 LOGICAL, INTENT(OUT) :: found !<750 750 CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< 751 751 CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< 752 752 CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< 753 753 754 LOGICAL, INTENT(OUT) :: found !< 755 756 754 757 found = .TRUE. 755 758 … … 757 760 ! 758 761 !-- s grid 759 CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz', &760 'wspeed', 'wspeed_xy', 'wspeed_xz', 'wspeed_yz', &762 CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz', & 763 'wspeed', 'wspeed_xy', 'wspeed_xz', 'wspeed_yz', & 761 764 'wdir', 'wdir_xy', 'wdir_xz', 'wdir_yz' ) 762 765 … … 787 790 ! 788 791 !-- w grid 789 CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz', &790 'wu', 'wu_xy', 'wu_xz', 'wu_yz', &791 'wv', 'wv_xy', 'wv_xz', 'wv_yz', &792 'wtheta', 'wtheta_xy', 'wtheta_xz', 'wtheta_yz', &793 'wq', 'wq_xy', 'wq_xz', 'wq_yz' 792 CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz', & 793 'wu', 'wu_xy', 'wu_xz', 'wu_yz', & 794 'wv', 'wv_xy', 'wv_xz', 'wv_yz', & 795 'wtheta', 'wtheta_xy', 'wtheta_xz', 'wtheta_yz', & 796 'wq', 'wq_xy', 'wq_xz', 'wq_yz' ) 794 797 795 798 grid_x = 'x' … … 808 811 END SUBROUTINE doq_define_netcdf_grid 809 812 810 !------------------------------------------------------------------------------ !813 !--------------------------------------------------------------------------------------------------! 811 814 ! 812 815 ! Description: 813 816 ! ------------ 814 817 !> Subroutine defining 2D output variables 815 !------------------------------------------------------------------------------ !816 SUBROUTINE doq_output_2d( av, variable, found, grid, 817 mode, local_pf, two_d, nzb_do, nzt_do,fill_value )818 !--------------------------------------------------------------------------------------------------! 819 SUBROUTINE doq_output_2d( av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do, & 820 fill_value ) 818 821 819 822 … … 841 844 REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output 842 845 846 843 847 flag_nr = 0 844 848 found = .TRUE. … … 852 856 to_be_resorted => ti 853 857 ELSE 854 IF ( .NOT. ALLOCATED( ti_av ) ) THEN858 IF ( .NOT. ALLOCATED( ti_av ) ) THEN 855 859 ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 856 860 ti_av = REAL( fill_value, KIND = wp ) … … 866 870 to_be_resorted => uu 867 871 ELSE 868 IF ( .NOT. ALLOCATED( uu_av ) ) THEN872 IF ( .NOT. ALLOCATED( uu_av ) ) THEN 869 873 ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 870 874 uu_av = REAL( fill_value, KIND = wp ) … … 880 884 to_be_resorted => vv 881 885 ELSE 882 IF ( .NOT. ALLOCATED( vv_av ) ) THEN886 IF ( .NOT. ALLOCATED( vv_av ) ) THEN 883 887 ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 884 888 vv_av = REAL( fill_value, KIND = wp ) … … 894 898 to_be_resorted => ww 895 899 ELSE 896 IF ( .NOT. ALLOCATED( ww_av ) ) THEN900 IF ( .NOT. ALLOCATED( ww_av ) ) THEN 897 901 ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 898 902 ww_av = REAL( fill_value, KIND = wp ) … … 908 912 to_be_resorted => wu 909 913 ELSE 910 IF ( .NOT. ALLOCATED( wu_av ) ) THEN914 IF ( .NOT. ALLOCATED( wu_av ) ) THEN 911 915 ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 912 916 wu_av = REAL( fill_value, KIND = wp ) … … 922 926 to_be_resorted => wv 923 927 ELSE 924 IF ( .NOT. ALLOCATED( wv_av ) ) THEN928 IF ( .NOT. ALLOCATED( wv_av ) ) THEN 925 929 ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 926 930 wv_av = REAL( fill_value, KIND = wp ) … … 936 940 to_be_resorted => wtheta 937 941 ELSE 938 IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN942 IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN 939 943 ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 940 944 wtheta_av = REAL( fill_value, KIND = wp ) … … 950 954 to_be_resorted => wq 951 955 ELSE 952 IF ( .NOT. ALLOCATED( wq_av ) ) THEN956 IF ( .NOT. ALLOCATED( wq_av ) ) THEN 953 957 ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 954 958 wq_av = REAL( fill_value, KIND = wp ) … … 968 972 ENDDO 969 973 ELSE 970 IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN974 IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN 971 975 ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) ) 972 976 pt_2m_av = REAL( fill_value, KIND = wp ) … … 990 994 ENDDO 991 995 ELSE 992 IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN996 IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN 993 997 ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) ) 994 998 uv_10m_av = REAL( fill_value, KIND = wp ) … … 1008 1012 to_be_resorted => wspeed 1009 1013 ELSE 1010 IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN1014 IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN 1011 1015 ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1012 1016 wspeed_av = REAL( fill_value, KIND = wp ) … … 1022 1026 to_be_resorted => wdir 1023 1027 ELSE 1024 IF ( .NOT. ALLOCATED( wdir_av ) ) THEN1028 IF ( .NOT. ALLOCATED( wdir_av ) ) THEN 1025 1029 ALLOCATE( wdir_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1026 1030 wdir_av = REAL( fill_value, KIND = wp ) … … 1042 1046 DO j = nys, nyn 1043 1047 DO k = nzb_do, nzt_do 1044 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &1045 REAL( fill_value, KIND = wp ),&1046 BTEST( wall_flags_total_0(k,j,i), flag_nr ) )1048 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), & 1049 REAL( fill_value, KIND = wp ), & 1050 BTEST( wall_flags_total_0(k,j,i), flag_nr ) ) 1047 1051 ENDDO 1048 1052 ENDDO … … 1053 1057 1054 1058 1055 !------------------------------------------------------------------------------ !1059 !--------------------------------------------------------------------------------------------------! 1056 1060 ! 1057 1061 ! Description: 1058 1062 ! ------------ 1059 1063 !> Subroutine defining 3D output variables 1060 !------------------------------------------------------------------------------! 1061 SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do, & 1062 nzt_do ) 1064 !--------------------------------------------------------------------------------------------------! 1065 SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do ) 1063 1066 1064 1067 IMPLICIT NONE … … 1082 1085 REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output 1083 1086 1087 1084 1088 flag_nr = 0 1085 1089 found = .TRUE. … … 1092 1096 to_be_resorted => ti 1093 1097 ELSE 1094 IF ( .NOT. ALLOCATED( ti_av ) ) THEN1098 IF ( .NOT. ALLOCATED( ti_av ) ) THEN 1095 1099 ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1096 1100 ti_av = REAL( fill_value, KIND = wp ) … … 1104 1108 to_be_resorted => uu 1105 1109 ELSE 1106 IF ( .NOT. ALLOCATED( uu_av ) ) THEN1110 IF ( .NOT. ALLOCATED( uu_av ) ) THEN 1107 1111 ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1108 1112 uu_av = REAL( fill_value, KIND = wp ) … … 1116 1120 to_be_resorted => vv 1117 1121 ELSE 1118 IF ( .NOT. ALLOCATED( vv_av ) ) THEN1122 IF ( .NOT. ALLOCATED( vv_av ) ) THEN 1119 1123 ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1120 1124 vv_av = REAL( fill_value, KIND = wp ) … … 1128 1132 to_be_resorted => ww 1129 1133 ELSE 1130 IF ( .NOT. ALLOCATED( ww_av ) ) THEN1134 IF ( .NOT. ALLOCATED( ww_av ) ) THEN 1131 1135 ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1132 1136 ww_av = REAL( fill_value, KIND = wp ) … … 1140 1144 to_be_resorted => wu 1141 1145 ELSE 1142 IF ( .NOT. ALLOCATED( wu_av ) ) THEN1146 IF ( .NOT. ALLOCATED( wu_av ) ) THEN 1143 1147 ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1144 1148 wu_av = REAL( fill_value, KIND = wp ) … … 1152 1156 to_be_resorted => wv 1153 1157 ELSE 1154 IF ( .NOT. ALLOCATED( wv_av ) ) THEN1158 IF ( .NOT. ALLOCATED( wv_av ) ) THEN 1155 1159 ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1156 1160 wv_av = REAL( fill_value, KIND = wp ) … … 1164 1168 to_be_resorted => wtheta 1165 1169 ELSE 1166 IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN1170 IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN 1167 1171 ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1168 1172 wtheta_av = REAL( fill_value, KIND = wp ) … … 1176 1180 to_be_resorted => wq 1177 1181 ELSE 1178 IF ( .NOT. ALLOCATED( wq_av ) ) THEN1182 IF ( .NOT. ALLOCATED( wq_av ) ) THEN 1179 1183 ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1180 1184 wq_av = REAL( fill_value, KIND = wp ) … … 1188 1192 to_be_resorted => wspeed 1189 1193 ELSE 1190 IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN1194 IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN 1191 1195 ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1192 1196 wspeed_av = REAL( fill_value, KIND = wp ) … … 1200 1204 to_be_resorted => wdir 1201 1205 ELSE 1202 IF ( .NOT. ALLOCATED( wdir_av ) ) THEN1206 IF ( .NOT. ALLOCATED( wdir_av ) ) THEN 1203 1207 ALLOCATE( wdir_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1204 1208 wdir_av = REAL( fill_value, KIND = wp ) … … 1217 1221 DO j = nys, nyn 1218 1222 DO k = nzb_do, nzt_do 1219 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &1220 REAL( fill_value, KIND = wp ),&1221 BTEST( wall_flags_total_0(k,j,i), flag_nr ) )1223 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), & 1224 REAL( fill_value, KIND = wp ), & 1225 BTEST( wall_flags_total_0(k,j,i), flag_nr ) ) 1222 1226 ENDDO 1223 1227 ENDDO … … 1227 1231 END SUBROUTINE doq_output_3d 1228 1232 1233 1234 !--------------------------------------------------------------------------------------------------! 1235 ! 1229 1236 ! Description: 1230 1237 ! ------------ 1231 !> Resorts the user-defined output quantity with indices (k,j,i) to a 1232 !> temporary array with indices(i,j,k) for masked data output.1233 !------------------------------------------------------------------------------ !1238 !> Resorts the user-defined output quantity with indices (k,j,i) to a temporary array with indices 1239 !> (i,j,k) for masked data output. 1240 !--------------------------------------------------------------------------------------------------! 1234 1241 SUBROUTINE doq_output_mask( av, variable, found, local_pf, mid ) 1235 1242 … … 1239 1246 1240 1247 IMPLICIT NONE 1248 1249 REAL(wp), PARAMETER :: fill_value = -9999.0_wp !< value for the _FillValue attribute 1241 1250 1242 1251 CHARACTER (LEN=*) :: variable !< … … 1257 1266 LOGICAL :: resorted !< true if array is resorted 1258 1267 1259 REAL(wp), & 1260 DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: & 1261 local_pf !< 1262 REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output 1263 1264 REAL(wp), PARAMETER :: fill_value = -9999.0_wp !< value for the _FillValue attribute 1268 REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: local_pf !< 1269 1270 REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output 1271 1265 1272 1266 1273 flag_nr = 0 … … 1373 1380 DO j = 1, mask_size_l(mid,2) 1374 1381 DO k = 1, mask_size_l(mid,3) 1375 local_pf(i,j,k) = MERGE( to_be_resorted(mask_k(mid,k), & 1376 mask_j(mid,j), & 1377 mask_i(mid,i)), & 1378 REAL( fill_value, KIND = wp ), & 1379 BTEST( wall_flags_total_0( & 1380 mask_k(mid,k), & 1381 mask_j(mid,j), & 1382 mask_i(mid,i)), & 1382 local_pf(i,j,k) = MERGE( to_be_resorted(mask_k(mid,k), & 1383 mask_j(mid,j), & 1384 mask_i(mid,i)), & 1385 REAL( fill_value, KIND = wp ), & 1386 BTEST( wall_flags_total_0(mask_k(mid,k), & 1387 mask_j(mid,j), & 1388 mask_i(mid,i)), & 1383 1389 flag_nr ) ) 1384 1390 ENDDO … … 1395 1401 im = mask_i(mid,i) 1396 1402 jm = mask_j(mid,j) 1397 ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), & 1398 DIM = 1 ) - 1 1403 ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 ) ), DIM=1 ) - 1 1399 1404 DO k = 1, mask_size_l(mid,3) 1400 1405 kk = MIN( ktt+mask_k(mid,k), nzt+1 ) … … 1415 1420 END SUBROUTINE doq_output_mask 1416 1421 1417 !------------------------------------------------------------------------------ !1422 !--------------------------------------------------------------------------------------------------! 1418 1423 ! Description: 1419 1424 ! ------------ 1420 1425 !> Allocate required arrays 1421 !------------------------------------------------------------------------------ !1426 !--------------------------------------------------------------------------------------------------! 1422 1427 SUBROUTINE doq_init 1423 1428 … … 1425 1430 1426 1431 INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities 1432 1427 1433 1428 1434 ! … … 1564 1570 1565 1571 ! 1566 !-- Save timestep number to check in time_integration if doq_calculate 1567 !-- has been called already, since the CALL occurs at two locations, but the calculations need to be1568 !-- done only once pertimestep.1572 !-- Save timestep number to check in time_integration if doq_calculate has been called already, 1573 !-- since the CALL occurs at two locations, but the calculations need to be done only once per 1574 !-- timestep. 1569 1575 timestep_number_at_prev_calc = current_timestep_number 1570 1576 … … 1580 1586 DO j = nys, nyn 1581 1587 DO k = nzb+1, nzt 1582 ti(k,j,i) = 0.25_wp * SQRT( &1583 ( ( w(k,j+1,i) + w(k-1,j+1,i)&1584 - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy&1585 - ( v(k+1,j,i) + v(k+1,j+1,i)&1586 - v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2&1587 + ( ( u(k+1,j,i) + u(k+1,j,i+1)&1588 - u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k)&1589 - ( w(k,j,i+1) + w(k-1,j,i+1)&1590 - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx )**2&1591 + ( ( v(k,j,i+1) + v(k,j+1,i+1)&1592 - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx&1593 - ( u(k,j+1,i) + u(k,j+1,i+1)&1594 - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy )**2 )&1595 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0) )1588 ti(k,j,i) = 0.25_wp * SQRT( & 1589 ( ( w(k,j+1,i) + w(k-1,j+1,i) & 1590 - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy & 1591 - ( v(k+1,j,i) + v(k+1,j+1,i) & 1592 - v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2 & 1593 + ( ( u(k+1,j,i) + u(k+1,j,i+1) & 1594 - u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k) & 1595 - ( w(k,j,i+1) + w(k-1,j,i+1) & 1596 - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx )**2 & 1597 + ( ( v(k,j,i+1) + v(k,j+1,i+1) & 1598 - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx & 1599 - ( u(k,j+1,i) + u(k,j+1,i+1) & 1600 - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy )**2 ) & 1601 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 1596 1602 ENDDO 1597 1603 ENDDO … … 1603 1609 DO j = nys, nyn 1604 1610 DO k = nzb+1, nzt 1605 uu(k,j,i) = u(k,j,i) * u(k,j,i) & 1606 * MERGE( 1.0_wp, 0.0_wp, & 1607 BTEST( wall_flags_total_0(k,j,i), 1) ) 1611 uu(k,j,i) = u(k,j,i) * u(k,j,i) & 1612 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 1608 1613 ENDDO 1609 1614 ENDDO … … 1615 1620 DO j = nys, nyn 1616 1621 DO k = nzb+1, nzt 1617 vv(k,j,i) = v(k,j,i) * v(k,j,i) & 1618 * MERGE( 1.0_wp, 0.0_wp, & 1619 BTEST( wall_flags_total_0(k,j,i), 2) ) 1622 vv(k,j,i) = v(k,j,i) * v(k,j,i) & 1623 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) ) 1620 1624 ENDDO 1621 1625 ENDDO … … 1627 1631 DO j = nys, nyn 1628 1632 DO k = nzb+1, nzt-1 1629 ww(k,j,i) = w(k,j,i) * w(k,j,i) & 1630 * MERGE( 1.0_wp, 0.0_wp, & 1631 BTEST( wall_flags_total_0(k,j,i), 3) ) 1633 ww(k,j,i) = w(k,j,i) * w(k,j,i) & 1634 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) ) 1632 1635 ENDDO 1633 1636 ENDDO … … 1639 1642 DO j = nys, nyn 1640 1643 DO k = nzb+1, nzt-1 1641 wu(k,j,i) = w(k,j,i) & 1642 * 0.25_wp * ( u(k,j,i) + u(k,j,i+1) & 1643 + u(k+1,j,i) + u(k+1,j,i+1) ) & 1644 * MERGE( 1.0_wp, 0.0_wp, & 1645 BTEST( wall_flags_total_0(k,j,i), 0) ) 1644 wu(k,j,i) = w(k,j,i) & 1645 * 0.25_wp * ( u(k,j,i) + u(k,j,i+1) + u(k+1,j,i) + u(k+1,j,i+1) )& 1646 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 1646 1647 ENDDO 1647 1648 ENDDO … … 1653 1654 DO j = nys, nyn 1654 1655 DO k = nzb+1, nzt-1 1655 wv(k,j,i) = w(k,j,i) & 1656 * 0.25_wp * ( v(k,j,i) + v(k,j+1,i) & 1657 + v(k+1,j,i) + v(k+1,j+1,i) ) & 1658 * MERGE( 1.0_wp, 0.0_wp, & 1659 BTEST( wall_flags_total_0(k,j,i), 0) ) 1656 wv(k,j,i) = w(k,j,i) & 1657 * 0.25_wp * ( v(k,j,i) + v(k,j+1,i) + v(k+1,j,i) + v(k+1,j+1,i) )& 1658 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 1660 1659 ENDDO 1661 1660 ENDDO … … 1667 1666 DO j = nys, nyn 1668 1667 DO k = nzb+1, nzt-1 1669 wtheta(k,j,i) = w(k,j,i) & 1670 * 0.5_wp * ( pt(k,j,i) + pt(k+1,j,i) ) & 1671 * MERGE( 1.0_wp, 0.0_wp, & 1672 BTEST( wall_flags_total_0(k,j,i), 3) ) 1668 wtheta(k,j,i) = w(k,j,i) & 1669 * 0.5_wp * ( pt(k,j,i) + pt(k+1,j,i) ) & 1670 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 )) 1673 1671 ENDDO 1674 1672 ENDDO … … 1680 1678 DO j = nys, nyn 1681 1679 DO k = nzb+1, nzt-1 1682 wq(k,j,i) = w(k,j,i) * 0.5_wp * ( q(k,j,i) + q(k+1,j,i) )& 1683 * MERGE( 1.0_wp, 0.0_wp, & 1684 BTEST( wall_flags_total_0(k,j,i), 3) ) 1680 wq(k,j,i) = w(k,j,i) * 0.5_wp * ( q(k,j,i) + q(k+1,j,i) ) & 1681 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) ) 1685 1682 ENDDO 1686 1683 ENDDO … … 1690 1687 CASE ( 'theta_2m*' ) 1691 1688 ! 1692 !-- 2-m potential temperature is caluclated from surface arrays. In 1693 !-- case the 2m level is below the first grid point, MOST is applied, 1694 !-- else, linear interpolation between two vertical grid levels is 1695 !-- applied. To access all surfaces, iterate over all horizontally- 1689 !-- 2-m potential temperature is caluclated from surface arrays. In case the 2m level is 1690 !-- below the first grid point, MOST is applied, else, linear interpolation between two 1691 !-- vertical grid levels is applied. To access all surfaces, iterate over all horizontally- 1696 1692 !-- upward facing surface types. 1697 1693 surf => surf_def_h(0) … … 1705 1701 CASE ( 'wspeed_10m*' ) 1706 1702 ! 1707 !-- 10-m wind speed is caluclated from surface arrays. In 1708 !-- case the 10m level is below the first grid point, MOST is applied, 1709 !-- else, linear interpolation between two vertical grid levels is 1710 !-- applied. To access all surfaces, iterate over all horizontally- 1711 !-- upward facing surface types. 1703 !-- 10-m wind speed is caluclated from surface arrays. In case the 10m level is below the 1704 !-- first grid point, MOST is applied, else, linear interpolation between two vertical grid 1705 !-- levels is applied. To access all surfaces, iterate over all horizontally-upward facing 1706 !-- surface types. 1712 1707 surf => surf_def_h(0) 1713 1708 CALL calc_wind_10m … … 1717 1712 CALL calc_wind_10m 1718 1713 ! 1719 !-- horizontal wind speed1714 !-- Horizontal wind speed 1720 1715 CASE ( 'wspeed' ) 1721 1716 DO i = nxl, nxr 1722 1717 DO j = nys, nyn 1723 1718 DO k = nzb, nzt+1 1724 wspeed(k,j,i) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 &1725 + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 ) &1726 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0))1719 wspeed(k,j,i) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 & 1720 + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 ) & 1721 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 )) 1727 1722 ENDDO 1728 1723 ENDDO … … 1730 1725 1731 1726 ! 1732 !-- horizontal wind direction1727 !-- Horizontal wind direction 1733 1728 CASE ( 'wdir' ) 1734 1729 DO i = nxl, nxr … … 1738 1733 v_center(k,j,i) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) 1739 1734 1740 wdir(k,j,i) = ATAN2( u_center(k,j,i), v_center(k,j,i) ) &1741 / pi * 180.0_wp + 180.0_wp1735 wdir(k,j,i) = ATAN2( u_center(k,j,i), v_center(k,j,i) ) & 1736 / pi * 180.0_wp + 180.0_wp 1742 1737 ENDDO 1743 1738 ENDDO … … 1752 1747 1753 1748 ! 1754 !-- The following block contains subroutines to calculate diagnostic 1755 !-- surface quantities. 1749 !-- The following block contains subroutines to calculate diagnostic surface quantities. 1756 1750 CONTAINS 1757 !------------------------------------------------------------------------------ !1751 !--------------------------------------------------------------------------------------------------! 1758 1752 ! Description: 1759 1753 ! ------------ 1760 1754 !> Calculation of 2-m potential temperature. 1761 !------------------------------------------------------------------------------ !1755 !--------------------------------------------------------------------------------------------------! 1762 1756 SUBROUTINE calc_pt_2m 1763 1757 1764 USE surface_layer_fluxes_mod, &1758 USE surface_layer_fluxes_mod, & 1765 1759 ONLY: psi_h 1766 1760 … … 1776 1770 k = surf%k(m) 1777 1771 ! 1778 !-- If 2-m level is below the first grid level, MOST is 1779 !-- used for calculation of2-m temperature.1772 !-- If 2-m level is below the first grid level, MOST is used for calculation of 1773 !-- 2-m temperature. 1780 1774 IF ( surf%z_mo(m) > 2.0_wp ) THEN 1781 pt_2m(j,i) = surf%pt_surface(m) + surf%ts(m) / kappa & 1782 * ( LOG( 2.0_wp / surf%z0h(m) ) & 1783 - psi_h( 2.0_wp / surf%ol(m) ) & 1784 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1785 ! 1786 !-- If 2-m level is above the first grid level, 2-m temperature 1787 !-- is linearly interpolated between the two nearest vertical grid 1788 !-- levels. Note, since 2-m temperature is only computed for 1789 !-- horizontal upward-facing surfaces, only a vertical 1790 !-- interpolation is necessary. 1775 pt_2m(j,i) = surf%pt_surface(m) + surf%ts(m) / kappa & 1776 * ( LOG( 2.0_wp / surf%z0h(m) ) & 1777 - psi_h( 2.0_wp / surf%ol(m) ) & 1778 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1779 ! 1780 !-- If 2-m level is above the first grid level, 2-m temperature is linearly interpolated 1781 !-- between the two nearest vertical grid levels. Note, since 2-m temperature is only 1782 !-- computed for horizontal upward-facing surfaces, only a vertical interpolation is 1783 !-- necessary. 1791 1784 ELSE 1792 1785 ! … … 1798 1791 ! 1799 1792 !-- kk defines the index of the first grid level >= 2m. 1800 pt_2m(j,i) = pt(kk-1,j,i) + &1801 ( zw(k-1) + 2.0_wp - zu(kk-1) ) * &1802 ( pt(kk,j,i) - pt(kk-1,j,i) ) / &1793 pt_2m(j,i) = pt(kk-1,j,i) + & 1794 ( zw(k-1) + 2.0_wp - zu(kk-1) ) * & 1795 ( pt(kk,j,i) - pt(kk-1,j,i) ) / & 1803 1796 ( zu(kk) - zu(kk-1) ) 1804 1797 ENDIF … … 1808 1801 END SUBROUTINE calc_pt_2m 1809 1802 1810 !------------------------------------------------------------------------------ !1803 !--------------------------------------------------------------------------------------------------! 1811 1804 ! Description: 1812 1805 ! ------------ 1813 1806 !> Calculation of 10-m wind speed. 1814 !------------------------------------------------------------------------------ !1807 !--------------------------------------------------------------------------------------------------! 1815 1808 SUBROUTINE calc_wind_10m 1816 1809 1817 USE surface_layer_fluxes_mod, &1810 USE surface_layer_fluxes_mod, & 1818 1811 ONLY: psi_m 1819 1812 … … 1825 1818 REAL(wp) :: uv_l !< wind speed at lower grid point 1826 1819 REAL(wp) :: uv_u !< wind speed at upper grid point 1820 1827 1821 1828 1822 DO m = 1, surf%ns … … 1832 1826 k = surf%k(m) 1833 1827 ! 1834 !-- If 10-m level is below the first grid level, MOST is 1835 !-- used for calculation of 10-mtemperature.1828 !-- If 10-m level is below the first grid level, MOST is used for calculation of 10-m 1829 !-- temperature. 1836 1830 IF ( surf%z_mo(m) > 10.0_wp ) THEN 1837 uv_10m(j,i) = surf%us(m) / kappa & 1838 * ( LOG( 10.0_wp / surf%z0(m) ) & 1839 - psi_m( 10.0_wp / surf%ol(m) ) & 1840 + psi_m( surf%z0(m) / surf%ol(m) ) ) 1841 ! 1842 !-- If 10-m level is above the first grid level, 10-m wind speed 1843 !-- is linearly interpolated between the two nearest vertical grid 1844 !-- levels. Note, since 10-m temperature is only computed for 1845 !-- horizontal upward-facing surfaces, only a vertical 1846 !-- interpolation is necessary. 1831 uv_10m(j,i) = surf%us(m) / kappa & 1832 * ( LOG( 10.0_wp / surf%z0(m) ) & 1833 - psi_m( 10.0_wp / surf%ol(m) ) & 1834 + psi_m( surf%z0(m) / surf%ol(m) ) ) 1835 ! 1836 !-- If 10-m level is above the first grid level, 10-m wind speed is linearly interpolated 1837 !-- between the two nearest vertical grid levels. Note, since 10-m temperature is only 1838 !-- computed for horizontal upward-facing surfaces, only a vertical interpolation is 1839 !-- necessary. 1847 1840 ELSE 1848 1841 ! … … 1854 1847 ! 1855 1848 !-- kk defines the index of the first grid level >= 10m. 1856 uv_l = SQRT( ( 0.5_wp * ( u(kk-1,j,i) + u(kk-1,j,i+1) ) )**2 &1849 uv_l = SQRT( ( 0.5_wp * ( u(kk-1,j,i) + u(kk-1,j,i+1) ) )**2 & 1857 1850 + ( 0.5_wp * ( v(kk-1,j,i) + v(kk-1,j+1,i) ) )**2 ) 1858 1851 1859 uv_u = SQRT( ( 0.5_wp * ( u(kk,j,i) + u(kk,j,i+1) ) )**2 &1852 uv_u = SQRT( ( 0.5_wp * ( u(kk,j,i) + u(kk,j,i+1) ) )**2 & 1860 1853 + ( 0.5_wp * ( v(kk,j,i) + v(kk,j+1,i) ) )**2 ) 1861 1854 1862 uv_10m(j,i) = uv_l + ( zw(k-1) + 10.0_wp - zu(kk-1) ) * &1863 ( uv_u - uv_l ) / &1855 uv_10m(j,i) = uv_l + ( zw(k-1) + 10.0_wp - zu(kk-1) ) * & 1856 ( uv_u - uv_l ) / & 1864 1857 ( zu(kk) - zu(kk-1) ) 1865 1858 … … 1873 1866 1874 1867 1875 !------------------------------------------------------------------------------ !1868 !--------------------------------------------------------------------------------------------------! 1876 1869 ! Description: 1877 1870 ! ------------ 1878 !> Preparation of the diagnostic output, counting of the module-specific 1879 !> output quantities andgathering of the output names.1880 !------------------------------------------------------------------------------ !1871 !> Preparation of the diagnostic output, counting of the module-specific output quantities and 1872 !> gathering of the output names. 1873 !--------------------------------------------------------------------------------------------------! 1881 1874 SUBROUTINE doq_prepare 1882 1875 1883 USE control_parameters, &1876 USE control_parameters, & 1884 1877 ONLY: do2d, do3d, domask, masks 1885 1878 1886 1879 IMPLICIT NONE 1887 1880 1888 CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) :: do2d_var = ' ' !< 1889 !< label array for 2d output quantities 1881 CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) :: do2d_var = ' ' !< label array for 2d output quantities 1890 1882 1891 1883 INTEGER(iwp) :: av !< index defining type of output, av=0 instantaneous, av=1 averaged … … 1895 1887 INTEGER(iwp) :: mid !< masked output running index 1896 1888 1889 1897 1890 prepared_diagnostic_output_quantities = .FALSE. 1898 1891 … … 1907 1900 ! 1908 1901 !-- Gather 2d output quantity names. 1909 !-- Check for double occurrence of output quantity, e.g. by _xy, 1910 !-- _yz, _xz. 1902 !-- Check for double occurrence of output quantity, e.g. by _xy, _yz, _xz. 1911 1903 DO WHILE ( do2d_var(av,ivar)(1:1) /= ' ' ) 1912 1904 IF ( .NOT. ANY( do_all == do2d_var(av,ivar) ) ) THEN … … 1930 1922 ivar = 1 1931 1923 ! 1932 !-- Gather masked output quantity names. Also check for double output 1933 !-- e.g. by different masks. 1924 !-- Gather masked output quantity names. Also check for double output e.g. by different masks. 1934 1925 DO mid = 1, masks 1935 1926 DO WHILE ( domask(mid,av,ivar)(1:1) /= ' ' ) … … 1950 1941 END SUBROUTINE doq_prepare 1951 1942 1952 !------------------------------------------------------------------------------ !1943 !--------------------------------------------------------------------------------------------------! 1953 1944 ! Description: 1954 1945 ! ------------ 1955 1946 !> Subroutine reads local (subdomain) restart data 1956 !------------------------------------------------------------------------------! 1957 SUBROUTINE doq_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, & 1958 nxr_on_file, nynf, nync, nyn_on_file, nysf, & 1959 nysc, nys_on_file, tmp_2d, tmp_3d, found ) 1947 !--------------------------------------------------------------------------------------------------! 1948 SUBROUTINE doq_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync, & 1949 nyn_on_file, nysf, nysc, nys_on_file, tmp_2d, tmp_3d, found ) 1960 1950 1961 1951 USE control_parameters … … 2002 1992 pt_2m_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 2003 1993 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 2004 1994 2005 1995 CASE ( 'ti_av' ) 2006 1996 IF ( .NOT. ALLOCATED( ti_av ) ) THEN … … 2108 2098 END SUBROUTINE doq_rrd_local_ftn 2109 2099 2110 !------------------------------------------------------------------------------ !2100 !--------------------------------------------------------------------------------------------------! 2111 2101 ! Description: 2112 2102 ! ------------ 2113 2103 !> Read module-specific local restart data arrays (MPI-IO). 2114 !------------------------------------------------------------------------------ !2104 !--------------------------------------------------------------------------------------------------! 2115 2105 SUBROUTINE doq_rrd_local_mpi 2116 2106 … … 2198 2188 END SUBROUTINE doq_rrd_local_mpi 2199 2189 2200 !------------------------------------------------------------------------------ !2190 !--------------------------------------------------------------------------------------------------! 2201 2191 ! Description: 2202 2192 ! ------------ 2203 2193 !> Subroutine writes local (subdomain) restart data 2204 !------------------------------------------------------------------------------ !2194 !--------------------------------------------------------------------------------------------------! 2205 2195 SUBROUTINE doq_wrd_local 2206 2196 2207 2208 2197 IMPLICIT NONE 2198 2209 2199 2210 2200 IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN
Note: See TracChangeset
for help on using the changeset viewer.
