- Timestamp:
- Oct 26, 2020 10:23:38 AM (4 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 4 edited
-
diagnostic_output_quantities_mod.f90 (modified) (28 diffs)
-
dynamics_mod.f90 (modified) (5 diffs)
-
flow_statistics.f90 (modified) (8 diffs)
-
module_interface.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
r4671 r4757 24 24 ! ----------------- 25 25 ! $Id$ 26 ! Implement relative humidity as diagnostic output quantity 27 ! 28 ! 4671 2020-09-09 20:27:58Z pavelkrc 26 29 ! Implementation of downward facing USM and LSM surfaces 27 30 ! … … 110 113 USE arrays_3d, & 111 114 ONLY: ddzu, & 115 exner, & 116 hyp, & 112 117 pt, & 113 118 q, & … … 119 124 120 125 USE basic_constants_and_equations_mod, & 121 ONLY: kappa, pi 126 ONLY: kappa, & 127 pi, & 128 magnus, & 129 rd_d_rv 122 130 123 131 USE control_parameters, & 124 132 ONLY: current_timestep_number, & 125 133 data_output, & 134 data_output_pr, & 135 humidity, & 126 136 message_string, & 127 137 restart_data_format_output, & … … 153 163 USE kinds 154 164 165 USE profil_parameter, & 166 ONLY: dopr_index 167 155 168 USE restart_data_mpi_io_mod, & 156 169 ONLY: rd_mpi_io_check_array, & 157 170 rrd_mpi_io, & 158 171 wrd_mpi_io 172 173 USE statistics, & 174 ONLY: hom, & 175 pr_palm, & 176 rmask, & 177 statistic_regions, & 178 sums_l 159 179 160 180 USE surface_mod, & … … 164 184 surf_usm_h 165 185 166 167 186 IMPLICIT NONE 168 187 … … 179 198 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uv_10m_av !< averaged horizontal wind speed at 10m 180 199 200 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: rh !< relative humidity 201 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: rh_av !< avg. relative humidity 181 202 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti !< rotation(u,v,w) aka turbulence intensity 182 203 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti_av !< avg. rotation(u,v,w) aka turbulence intensity … … 220 241 doq_calculate, & 221 242 doq_check_data_output, & 243 doq_check_data_output_pr, & 222 244 doq_define_netcdf_grid, & 223 245 doq_init, & … … 225 247 doq_output_3d, & 226 248 doq_output_mask, & 249 doq_statistics, & 227 250 doq_rrd_local, & 228 251 doq_wrd_local … … 241 264 END INTERFACE doq_check_data_output 242 265 266 INTERFACE doq_check_data_output_pr 267 MODULE PROCEDURE doq_check_data_output_pr 268 END INTERFACE doq_check_data_output_pr 269 243 270 INTERFACE doq_define_netcdf_grid 244 271 MODULE PROCEDURE doq_define_netcdf_grid … … 260 287 MODULE PROCEDURE doq_init 261 288 END INTERFACE doq_init 289 290 INTERFACE doq_statistics 291 MODULE PROCEDURE doq_statistics 292 END INTERFACE doq_statistics 262 293 263 294 INTERFACE doq_prepare … … 274 305 END INTERFACE doq_wrd_local 275 306 276 277 307 CONTAINS 278 308 … … 300 330 SELECT CASE ( TRIM( variable ) ) 301 331 332 CASE ( 'rh' ) 333 IF ( .NOT. ALLOCATED( rh_av ) ) THEN 334 ALLOCATE( rh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 335 ENDIF 336 rh_av = 0.0_wp 337 302 338 CASE ( 'ti' ) 303 339 IF ( .NOT. ALLOCATED( ti_av ) ) THEN … … 384 420 385 421 SELECT CASE ( TRIM( variable ) ) 422 423 CASE ( 'rh' ) 424 IF ( ALLOCATED( rh_av ) ) THEN 425 DO i = nxl, nxr 426 DO j = nys, nyn 427 DO k = nzb, nzt+1 428 rh_av(k,j,i) = rh_av(k,j,i) + rh(k,j,i) 429 ENDDO 430 ENDDO 431 ENDDO 432 ENDIF 386 433 387 434 CASE ( 'ti' ) … … 523 570 SELECT CASE ( TRIM( variable ) ) 524 571 572 CASE ( 'rh' ) 573 IF ( ALLOCATED( rh_av ) ) THEN 574 DO i = nxl, nxr 575 DO j = nys, nyn 576 DO k = nzb, nzt+1 577 rh_av(k,j,i) = rh_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 578 ENDDO 579 ENDDO 580 ENDDO 581 ENDIF 582 525 583 CASE ( 'ti' ) 526 584 IF ( ALLOCATED( ti_av ) ) THEN … … 688 746 SELECT CASE ( TRIM( var ) ) 689 747 748 CASE ( 'rh' ) 749 unit = '%' 750 690 751 CASE ( 'ti' ) 691 752 unit = '1/s' … … 740 801 END SUBROUTINE doq_check_data_output 741 802 803 804 !--------------------------------------------------------------------------------------------------! 805 ! Description: 806 ! ------------ 807 !> Set the unit of user defined profile output quantities. For those variables not recognized by the 808 !> user, the parameter unit is set to "illegal", which tells the calling routine that the 809 !> output variable is not defined and leads to a program abort. 810 !--------------------------------------------------------------------------------------------------! 811 SUBROUTINE doq_check_data_output_pr( variable, var_count, unit, dopr_unit ) 812 813 CHARACTER (LEN=*) :: unit !< 814 CHARACTER (LEN=*) :: variable !< 815 CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit 816 817 INTEGER(iwp) :: pr_index !< 818 INTEGER(iwp) :: var_count !< 819 820 SELECT CASE ( TRIM( variable ) ) 821 822 CASE ( 'rh' ) 823 IF ( .NOT. humidity ) THEN 824 message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // & 825 ' requires humidity' 826 CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 ) 827 ENDIF 828 pr_index = 130 829 dopr_index(var_count) = pr_index 830 dopr_unit = '%' 831 unit = dopr_unit 832 hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 ) 833 834 CASE DEFAULT 835 unit = 'illegal' 836 837 END SELECT 838 839 840 END SUBROUTINE doq_check_data_output_pr 841 842 742 843 !--------------------------------------------------------------------------------------------------! 743 844 ! … … 763 864 ! 764 865 !-- s grid 765 CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz', & 866 CASE ( 'rh', 'rh_xy', 'rh_xz', 'rh_yz', & 867 'ti', 'ti_xy', 'ti_xz', 'ti_yz', & 766 868 'wspeed', 'wspeed_xy', 'wspeed_xz', 'wspeed_yz', & 767 869 'wdir', 'wdir_xy', 'wdir_xz', 'wdir_yz' ) … … 854 956 855 957 SELECT CASE ( TRIM( variable ) ) 958 959 CASE ( 'rh_xy', 'rh_xz', 'rh_yz' ) 960 IF ( av == 0 ) THEN 961 to_be_resorted => rh 962 ELSE 963 IF ( .NOT. ALLOCATED( rh_av ) ) THEN 964 ALLOCATE( rh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 965 rh_av = REAL( fill_value, KIND = wp ) 966 ENDIF 967 to_be_resorted => rh_av 968 ENDIF 969 flag_nr = 0 970 971 IF ( mode == 'xy' ) grid = 'zu' 856 972 857 973 CASE ( 'ti_xy', 'ti_xz', 'ti_yz' ) … … 1095 1211 SELECT CASE ( TRIM( variable ) ) 1096 1212 1213 CASE ( 'rh' ) 1214 IF ( av == 0 ) THEN 1215 to_be_resorted => rh 1216 ELSE 1217 IF ( .NOT. ALLOCATED( rh ) ) THEN 1218 ALLOCATE( rh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1219 rh = REAL( fill_value, KIND = wp ) 1220 ENDIF 1221 to_be_resorted => rh 1222 ENDIF 1223 flag_nr = 0 1224 1097 1225 CASE ( 'ti' ) 1098 1226 IF ( av == 0 ) THEN … … 1280 1408 1281 1409 SELECT CASE ( TRIM( variable ) ) 1410 1411 CASE ( 'rh' ) 1412 IF ( av == 0 ) THEN 1413 to_be_resorted => rh 1414 ELSE 1415 to_be_resorted => rh_av 1416 ENDIF 1417 grid = 's' 1418 flag_nr = 0 1282 1419 1283 1420 CASE ( 'ti' ) … … 1450 1587 SELECT CASE ( TRIM( do_all(ivar) ) ) 1451 1588 ! 1589 !-- Allocate array for 'relative humidity' 1590 CASE ( 'rh' ) 1591 IF ( .NOT. ALLOCATED( rh ) ) THEN 1592 ALLOCATE( rh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1593 rh = 0.0_wp 1594 ENDIF 1595 ! 1452 1596 !-- Allocate array for 'turbulence intensity' 1453 1597 CASE ( 'ti' ) … … 1556 1700 ! Description: 1557 1701 ! ------------ 1702 !> Calculation of module-specific statistics, i.e. horizontally averaged profiles and time series. 1703 !> This is called for every statistic region sr, but at least for the region "total domain" (sr=0). 1704 !--------------------------------------------------------------------------------------------------! 1705 SUBROUTINE doq_statistics( mode, sr, tn ) 1706 1707 1708 CHARACTER (LEN=*) :: mode !< 1709 1710 ! INTEGER(iwp) :: i !< 1711 ! INTEGER(iwp) :: j !< 1712 ! INTEGER(iwp) :: k !< 1713 INTEGER(iwp) :: sr !< 1714 INTEGER(iwp) :: tn !< 1715 ! 1716 !-- Next line is to avoid compiler warning about unused variables. Please remove. 1717 IF ( sr == 0 .OR. tn == 0 ) CONTINUE 1718 1719 IF ( mode == 'profiles' ) THEN 1720 1721 ELSEIF ( mode == 'time_series' ) THEN 1722 1723 ENDIF 1724 1725 END SUBROUTINE doq_statistics 1726 1727 1728 !--------------------------------------------------------------------------------------------------! 1729 ! Description: 1730 ! ------------ 1558 1731 !> Calculate diagnostic quantities 1559 1732 !--------------------------------------------------------------------------------------------------! … … 1567 1740 INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities 1568 1741 1742 REAL(wp) :: temp !< temperature 1743 REAL(wp) :: e_s !< saturation vapor pressure 1744 REAL(wp) :: q_s !< saturation mixing ratio 1745 1569 1746 TYPE(surf_type), POINTER :: surf !< surf-type array, used to generalize subroutines 1570 1747 … … 1583 1760 1584 1761 SELECT CASE ( TRIM( do_all(ivar) ) ) 1762 ! 1763 !-- rh (relative humidity) 1764 CASE ( 'rh' ) 1765 DO i = nxl, nxr 1766 DO j = nys, nyn 1767 DO k = nzb+1, nzt 1768 IF ( humidity ) THEN 1769 temp = exner(k) * pt(k,j,i) 1770 e_s = magnus( temp ) 1771 q_s = rd_d_rv * e_s / ( hyp(k) - e_s ) 1772 rh(k,j,i) = q(k,j,i) / q_s * 100.0_wp & 1773 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 1774 ENDIF 1775 ENDDO 1776 ENDDO 1777 ENDDO 1585 1778 ! 1586 1779 !-- Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point … … 1996 2189 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 1997 2190 2191 CASE ( 'rh_av' ) 2192 IF ( .NOT. ALLOCATED( rh_av ) ) THEN 2193 ALLOCATE( rh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 2194 ENDIF 2195 IF ( k == 1 ) READ ( 13 ) tmp_3d 2196 rh_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 2197 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 2198 1998 2199 CASE ( 'ti_av' ) 1999 2200 IF ( .NOT. ALLOCATED( ti_av ) ) THEN … … 2123 2324 ENDIF 2124 2325 2326 CALL rd_mpi_io_check_array( 'rh_av' , found = array_found ) 2327 IF ( array_found ) THEN 2328 IF ( .NOT. ALLOCATED( rh_av ) ) ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 2329 CALL rrd_mpi_io( 'rh_av', rh_av ) 2330 ENDIF 2331 2125 2332 CALL rd_mpi_io_check_array( 'u_center_av' , found = array_found ) 2126 2333 IF ( array_found ) THEN … … 2208 2415 ENDIF 2209 2416 2417 IF ( ALLOCATED( rh_av ) ) THEN 2418 CALL wrd_write_string( 'rh_av' ) 2419 WRITE ( 14 ) rh_av 2420 ENDIF 2421 2210 2422 IF ( ALLOCATED( ti_av ) ) THEN 2211 2423 CALL wrd_write_string( 'ti_av' ) … … 2271 2483 2272 2484 IF ( ALLOCATED( pt_2m_av ) ) CALL wrd_mpi_io( 'pt_2m_av', pt_2m_av ) 2485 IF ( ALLOCATED( rh_av ) ) CALL wrd_mpi_io( 'rh_av', rh_av ) 2273 2486 IF ( ALLOCATED( ti_av ) ) CALL wrd_mpi_io( 'ti_av', ti_av ) 2274 2487 IF ( ALLOCATED( u_center_av ) ) CALL wrd_mpi_io( 'u_center_av', u_center_av ) -
palm/trunk/SOURCE/dynamics_mod.f90
r4731 r4757 24 24 ! ----------------- 25 25 ! $Id$ 26 ! Implement relative humidity as diagnostic output quantity 27 ! 28 ! 4731 2020-10-07 13:25:11Z schwenkel 26 29 ! Move exchange_horiz from time_integration to modules 27 30 ! … … 82 85 u, u_1, u_2, u_init, u_p, u_m_l, u_m_n, u_m_r, u_m_s, & 83 86 v, v_1, v_2, v_p, v_init, v_m_l, v_m_n, v_m_r, v_m_s, & 84 w, w_1, w_2, w_p, w_m_l, w_m_n, w_m_r, w_m_s 87 w, w_1, w_2, w_p, w_m_l, w_m_n, w_m_r, w_m_s, & 88 zu 85 89 86 90 USE basic_constants_and_equations_mod, & … … 518 522 REAL(wp) :: q_s !< saturation mixing ratio 519 523 REAL(wp) :: t_l !< actual temperature 524 REAL(wp) :: rh = 9999999.9_wp !< relative humidity 525 REAL(wp) :: rh_min = 9999999.9_wp !< max relative humidity 526 REAL(wp) :: height = 9999999.9_wp !< height of supersaturated regions 527 REAL(wp) :: min_height = 9999999.9_wp !< height of supersaturated regions 520 528 521 529 ! … … 537 545 q_s = rd_d_rv * e_s / ( hyp(k) - e_s ) 538 546 539 IF ( q(k,j,i) > 1.02_wp * q_s ) realistic_q = .FALSE. 547 IF ( q(k,j,i) > 1.02_wp * q_s ) THEN 548 realistic_q = .FALSE. 549 rh = q(k,j,i) / q_s * 100.0_wp 550 height = zu(k) 551 ENDIF 540 552 ENDDO 541 553 ENDDO … … 546 558 #if defined( __parallel ) 547 559 CALL MPI_ALLREDUCE( MPI_IN_PLACE, realistic_q, 1, MPI_LOGICAL, MPI_LAND, comm2d, ierr) 560 CALL MPI_ALLREDUCE( rh, rh_min, 1, MPI_REAL, MPI_MIN, comm2d, ierr ) 561 CALL MPI_ALLREDUCE( height, min_height, 1, MPI_REAL, MPI_MIN, comm2d, ierr ) 548 562 #endif 549 563 550 564 IF ( .NOT. realistic_q .AND. & 551 565 TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 552 message_string = 'The initial mixing ratio exceeds the saturation mixing ratio.' 566 WRITE( message_string, * ) 'The initial mixing ratio exceeds the saturation mixing' // & 567 'ratio, with rh =', rh_min, '% in a height of', min_height, 'm for the first time' 553 568 CALL message( 'dynamic_init_checks', 'PA0697', 2, 2, 0, 6, 0 ) 554 569 ELSEIF ( .NOT. realistic_q .AND. & 555 570 TRIM( initializing_actions ) == 'read_restart_data' ) THEN 556 message_string = 'The mixing ratio exceeds the saturation mixing ratio.' 571 WRITE( message_string, * ) 'The initial mixing ratio exceeds the saturation mixing' // & 572 'ratio, with rh =', rh_min, '% in a height of', min_height, 'm for the first time' 557 573 CALL message( 'dynamic_init_checks', 'PA0697', 0, 1, 0, 6, 0 ) 558 574 ENDIF -
palm/trunk/SOURCE/flow_statistics.f90
r4742 r4757 24 24 ! ----------------- 25 25 ! $Id$ 26 ! Implement relative humidity as diagnostic output quantity 27 ! 28 ! 4742 2020-10-14 15:11:02Z schwenkel 26 29 ! Implement snow and graupel (bulk microphysics) 27 30 ! … … 105 108 106 109 USE arrays_3d, & 107 ONLY: ddzu, ddzw, d_exner, e, heatflux_output_conversion, hyp, km, kh,&110 ONLY: ddzu, ddzw, d_exner, e, exner, heatflux_output_conversion, hyp, km, kh, & 108 111 momentumflux_output_conversion, & 109 112 nc, ni, ng, nr, ns, p, prho, prr, pt, q, qc, qi, qg, ql, qr, qs, & … … 112 115 113 116 USE basic_constants_and_equations_mod, & 114 ONLY: g, lv_d_cp 117 ONLY: g, lv_d_cp, magnus, rd_d_rv 115 118 116 119 USE bulk_cloud_model_mod, & … … 220 223 REAL(wp) :: v2 !< 221 224 REAL(wp) :: w2 !< 225 226 REAL(wp) :: temp !< temperature 227 REAL(wp) :: e_s !< saturation vapor pressure 228 REAL(wp) :: q_s !< saturation mixing ratio 229 REAL(wp) :: rh_tmp !< local value of relative humidity 222 230 223 231 REAL(wp) :: dptdz(nzb+1:nzt+1) !< … … 1390 1398 END IF 1391 1399 END IF 1400 temp = exner(k) * pt(k,j,i) 1401 e_s = magnus( temp ) 1402 q_s = rd_d_rv * e_s / ( hyp(k) - e_s ) 1403 rh_tmp = q(k,j,i) / q_s * 100.0_wp 1404 sums_l(k,130,tn) = sums_l(k,130,tn) + rh_tmp * rmask(j,i,sr) * flag 1392 1405 ENDIF 1393 1406 ! … … 1747 1760 IF ( max_pr_cs > 0 ) THEN 1748 1761 sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+ max_pr_cs,0) = & 1749 sums_l(:,pr_palm+max_pr_user+1:pr_palm +max_pr_user+max_pr_cs,0) + &1750 sums_l(:,pr_palm+max_pr_user+1:pr_palm +max_pr_user+max_pr_cs,i)1762 sums_l(:,pr_palm+max_pr_user+1:pr_palm+max_pr_user+max_pr_cs,0) + & 1763 sums_l(:,pr_palm+max_pr_user+1:pr_palm+max_pr_user+max_pr_cs,i) 1751 1764 1752 1765 ENDIF … … 1832 1845 sums(k,116) = sums(k,116) / ngp_2dh_s_inner(k,sr) 1833 1846 sums(k,118:pr_palm-2) = sums(k,118:pr_palm-2) / ngp_2dh_s_inner(k,sr) 1834 sums(k,123:1 29) = sums(k,123:129) * ngp_2dh_s_inner(k,sr) / ngp_2dh(sr)1847 sums(k,123:130) = sums(k,123:130) * ngp_2dh_s_inner(k,sr) / ngp_2dh(sr) 1835 1848 ENDIF 1836 1849 ENDDO … … 1946 1959 hom(:,1,128,sr) = sums(:,128) ! ns 1947 1960 hom(:,1,129,sr) = sums(:,129) ! qs 1961 hom(:,1,130,sr) = sums(:,130) ! rh 1948 1962 hom(:,1,73,sr) = sums(:,73) ! nr 1949 1963 hom(:,1,74,sr) = sums(:,74) ! qr -
palm/trunk/SOURCE/module_interface.f90
r4753 r4757 24 24 ! ----------------- 25 25 ! $Id$ 26 ! Implement doq calls 27 ! 28 ! 4753 2020-10-21 14:55:41Z raasch 26 29 ! file re-formatted to follow the PALM coding standard 27 30 ! … … 342 345 ONLY: doq_3d_data_averaging, & 343 346 doq_check_data_output, & 347 doq_check_data_output_pr, & 344 348 doq_define_netcdf_grid, & 345 349 doq_init, & 346 350 doq_output_2d, & 347 351 doq_output_3d, & 352 doq_statistics, & 348 353 doq_rrd_local, & 349 354 doq_wrd_local … … 870 875 871 876 CALL dynamics_check_data_output_pr( variable, var_count, unit, dopr_unit ) 877 CALL doq_check_data_output_pr( variable, var_count, unit, dopr_unit ) 872 878 873 879 IF ( unit == 'illegal' .AND. bulk_cloud_model ) THEN … … 1701 1707 1702 1708 CALL dynamics_statistics( mode, sr, tn ) 1709 CALL doq_statistics( mode, sr, tn ) 1703 1710 1704 1711 IF ( gust_module_enabled ) CALL gust_statistics( mode, sr, tn, dots_max )
Note: See TracChangeset
for help on using the changeset viewer.
