Changeset 4757 for palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
- Timestamp:
- Oct 26, 2020 10:23:38 AM (4 years ago)
- File:
-
- 1 edited
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 )
Note: See TracChangeset
for help on using the changeset viewer.
