Changeset 4331 for palm/trunk
- Timestamp:
- Dec 10, 2019 6:25:02 PM (5 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 15 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r4309 r4331 25 25 # ----------------- 26 26 # $Id$ 27 # Move diagnostic surface output to diagnostic_output_quantities 28 # 29 # 4309 2019-11-26 18:49:59Z suehring 27 30 # Add dependency to parallel random generator for synthetic turbulence generator 28 31 # … … 519 522 diagnostic_output_quantities_mod.o: \ 520 523 mod_kinds.o \ 521 modules.o 524 modules.o \ 525 surface_layer_fluxes_mod.o 522 526 diffusion_s.o: \ 523 527 mod_kinds.o \ … … 1128 1132 basic_constants_and_equations_mod.o \ 1129 1133 cpulog_mod.o \ 1134 diagnostic_output_quantities_mod.o \ 1130 1135 land_surface_model_mod.o \ 1131 1136 mod_kinds.o \ -
palm/trunk/SOURCE/average_3d_data.f90
r4182 r4331 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Move 2-m potential temperature output to diagnostic_output_quantities 28 ! 29 ! 4182 2019-08-22 15:20:23Z scharf 27 30 ! Corrected "Former revisions" section 28 31 ! … … 315 318 ENDIF 316 319 317 CASE ( 'theta_2m*' )318 IF ( ALLOCATED( pt_2m_av ) ) THEN319 DO i = nxlg, nxrg320 DO j = nysg, nyng321 pt_2m_av(j,i) = pt_2m_av(j,i) / REAL( average_count_3d, KIND=wp )322 ENDDO323 ENDDO324 CALL exchange_horiz_2d( pt_2m_av, nbgp )325 ENDIF326 327 320 CASE ( 't*' ) 328 321 IF ( ALLOCATED( ts_av ) ) THEN -
palm/trunk/SOURCE/check_parameters.f90
r4301 r4331 25 25 ! ----------------- 26 26 ! 4172 2019-08-20 11:55:33Z oliver.maas 27 ! Move 2-m potential temperature output to diagnostic_output_quantities 28 ! 29 ! 11:55:33Z oliver.maas 27 30 ! removed message PA0421, concerning old parameter recycling_yshift 28 31 ! … … 2385 2388 2386 2389 CASE ( 'ghf*', 'lwp*', 'ol*', 'qsws*', 'r_a*', & 2387 'shf*', 'ssws*', 't*', 't heta_2m*', 'tsurf*', 'us*',&2390 'shf*', 'ssws*', 't*', 'tsurf*', 'us*', & 2388 2391 'z0*', 'z0h*', 'z0q*' ) 2389 2392 IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN … … 2424 2427 'res passive_scalar = .TRUE.' 2425 2428 CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 ) 2426 ENDIF2427 !2428 !-- Activate calculation of 2m temperature if output is requested2429 IF ( TRIM( var ) == 'theta_2m*' ) THEN2430 do_output_at_2m = .TRUE.2431 unit = 'K'2432 2429 ENDIF 2433 2430 -
palm/trunk/SOURCE/data_output_2d.f90
r4329 r4331 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Move 2-m potential temperature output to diagnostic_output_quantities 28 ! 29 ! 4329 2019-12-10 15:46:36Z motisi 27 30 ! Renamed wall_flags_0 to wall_flags_static_0 28 31 ! … … 1001 1004 ENDIF 1002 1005 IF ( mode == 'xy' ) level_z = zu 1003 1004 CASE ( 'theta_2m*_xy' ) ! 2d-array1005 IF ( av == 0 ) THEN1006 DO m = 1, surf_def_h(0)%ns1007 i = surf_def_h(0)%i(m)1008 j = surf_def_h(0)%j(m)1009 local_pf(i,j,nzb+1) = surf_def_h(0)%pt_2m(m)1010 ENDDO1011 DO m = 1, surf_lsm_h%ns1012 i = surf_lsm_h%i(m)1013 j = surf_lsm_h%j(m)1014 local_pf(i,j,nzb+1) = surf_lsm_h%pt_2m(m)1015 ENDDO1016 DO m = 1, surf_usm_h%ns1017 i = surf_usm_h%i(m)1018 j = surf_usm_h%j(m)1019 local_pf(i,j,nzb+1) = surf_usm_h%pt_2m(m)1020 ENDDO1021 ELSE1022 IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN1023 ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )1024 pt_2m_av = REAL( fill_value, KIND = wp )1025 ENDIF1026 DO i = nxl, nxr1027 DO j = nys, nyn1028 local_pf(i,j,nzb+1) = pt_2m_av(j,i)1029 ENDDO1030 ENDDO1031 ENDIF1032 resorted = .TRUE.1033 two_d = .TRUE.1034 level_z(nzb+1) = zu(nzb+1)1035 1006 1036 1007 CASE ( 'w_xy', 'w_xz', 'w_yz' ) -
palm/trunk/SOURCE/data_output_mask.f90
r4329 r4331 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Formatting adjustment 28 ! 29 ! 4329 2019-12-10 15:46:36Z motisi 27 30 ! Renamed wall_flags_0 to wall_flags_static_0 28 31 ! … … 603 606 IF ( .NOT. found ) THEN 604 607 CALL doq_output_mask( av, domask(mid,av,ivar), found, local_pf, & 605 mid)608 mid) 606 609 ENDIF 607 610 ! -
palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
r4329 r4331 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Modularize 2-m potential temperature output 28 ! - New output for 10-m wind speed 29 ! 30 ! 4329 2019-12-10 15:46:36Z motisi 27 31 ! Renamed wall_flags_0 to wall_flags_static_0 28 32 ! … … 75 79 76 80 USE arrays_3d, & 77 ONLY: ddzu, u, v, w 78 ! 79 ! USE averaging 80 ! 81 ! USE basic_constants_and_equations_mod, & 82 ! ONLY: c_p, lv_d_cp, l_v 83 ! 84 ! USE bulk_cloud_model_mod, & 85 ! ONLY: bulk_cloud_model 86 ! 81 ONLY: ddzu, & 82 pt, & 83 u, & 84 v, & 85 w, & 86 zu, & 87 zw 88 89 USE basic_constants_and_equations_mod, & 90 ONLY: kappa 91 87 92 USE control_parameters, & 88 ONLY: current_timestep_number, varnamelength 93 ONLY: current_timestep_number, & 94 data_output, & 95 message_string, & 96 varnamelength 89 97 ! 90 98 ! USE cpulog, & … … 93 101 USE grid_variables, & 94 102 ONLY: ddx, ddy 95 ! 103 96 104 USE indices, & 97 ONLY: nxl, nxr, nyn, nys, nzb, nzt, wall_flags_static_0 98 ! 105 ONLY: nbgp, & 106 nxl, & 107 nxlg, & 108 nxr, & 109 nxrg, & 110 nyn, & 111 nyng, & 112 nys, & 113 nysg, & 114 nzb, & 115 nzt, & 116 wall_flags_static_0 117 99 118 USE kinds 100 ! 101 ! USE land_surface_model_mod, & 102 ! ONLY: zs 103 ! 104 ! USE module_interface, & 105 ! ONLY: module_interface_data_output_2d 106 ! 107 ! #if defined( __netcdf ) 108 ! USE NETCDF 109 ! #endif 110 ! 111 ! USE netcdf_interface, & 112 ! ONLY: fill_value, id_set_xy, id_set_xz, id_set_yz, id_var_do2d, & 113 ! id_var_time_xy, id_var_time_xz, id_var_time_yz, nc_stat, & 114 ! netcdf_data_format, netcdf_handle_error 115 ! 116 ! USE particle_attributes, & 117 ! ONLY: grid_particles, number_of_particles, particle_advection_start, & 118 ! particles, prt_count 119 ! 120 ! USE pegrid 121 ! 122 ! USE surface_mod, & 123 ! ONLY: ind_pav_green, ind_veg_wall, ind_wat_win, surf_def_h, & 124 ! surf_lsm_h, surf_usm_h 125 ! 126 ! USE turbulence_closure_mod, & 127 ! ONLY: tcm_data_output_2d 119 120 USE surface_mod, & 121 ONLY: surf_def_h, & 122 surf_lsm_h, & 123 surf_type, & 124 surf_usm_h 128 125 129 126 … … 137 134 LOGICAL :: prepared_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is p 138 135 136 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt_2m !< 2-m air potential temperature 137 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt_2m_av !< averaged 2-m air potential temperature 138 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uv_10m !< horizontal wind speed at 10m 139 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uv_10m_av !< averaged horizontal wind speed at 10m 140 139 141 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti !< rotation(u,v,w) aka turbulence intensity 140 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti_av !< avg. rotation(u,v,w) aka turbulence intensity 141 142 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti_av !< avg. rotation(u,v,w) aka turbulence intensity 142 143 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu !< uu 143 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu_av !< mean of uu 144 144 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu_av !< mean of uu 145 145 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv !< vv 146 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv_av !< mean of vv 147 146 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv_av !< mean of vv 148 147 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww !< ww 149 148 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww_av !< mean of ww … … 160 159 prepared_diagnostic_output_quantities, & 161 160 timestep_number_at_prev_calc, & 161 pt_2m_av, & 162 162 ti_av, & 163 163 uu_av, & 164 uv_10m_av, & 164 165 vv_av, & 165 166 ww_av … … 270 271 ENDIF 271 272 ww_av = 0.0_wp 273 274 CASE ( 'theta_2m*' ) 275 IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN 276 ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) ) 277 ENDIF 278 pt_2m_av = 0.0_wp 279 280 CASE ( 'wspeed_10m*' ) 281 IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN 282 ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) ) 283 ENDIF 284 uv_10m_av = 0.0_wp 272 285 273 286 CASE DEFAULT … … 323 336 ENDDO 324 337 ENDIF 338 339 CASE ( 'theta_2m*' ) 340 IF ( ALLOCATED( pt_2m_av ) ) THEN 341 DO i = nxl, nxr 342 DO j = nys, nyn 343 pt_2m_av(j,i) = pt_2m_av(j,i) + pt_2m(j,i) 344 ENDDO 345 ENDDO 346 ENDIF 347 348 CASE ( 'wspeed_10m*' ) 349 IF ( ALLOCATED( uv_10m_av ) ) THEN 350 DO i = nxl, nxr 351 DO j = nys, nyn 352 uv_10m_av(j,i) = uv_10m_av(j,i) + uv_10m(j,i) 353 ENDDO 354 ENDDO 355 ENDIF 325 356 326 357 CASE DEFAULT … … 377 408 ENDIF 378 409 410 CASE ( 'theta_2m*' ) 411 IF ( ALLOCATED( pt_2m_av ) ) THEN 412 DO i = nxlg, nxrg 413 DO j = nysg, nyng 414 pt_2m_av(j,i) = pt_2m_av(j,i) / REAL( average_count_3d, KIND=wp ) 415 ENDDO 416 ENDDO 417 CALL exchange_horiz_2d( pt_2m_av, nbgp ) 418 ENDIF 419 420 CASE ( 'wspeed_10m*' ) 421 IF ( ALLOCATED( uv_10m_av ) ) THEN 422 DO i = nxlg, nxrg 423 DO j = nysg, nyng 424 uv_10m_av(j,i) = uv_10m_av(j,i) / REAL( average_count_3d, KIND=wp ) 425 ENDDO 426 ENDDO 427 CALL exchange_horiz_2d( uv_10m_av, nbgp ) 428 ENDIF 429 379 430 END SELECT 380 431 … … 389 440 !> Check data output for diagnostic output 390 441 !------------------------------------------------------------------------------! 391 SUBROUTINE doq_check_data_output( var, unit )442 SUBROUTINE doq_check_data_output( var, unit, i, ilen, k ) 392 443 393 444 IMPLICIT NONE … … 395 446 CHARACTER (LEN=*) :: unit !< 396 447 CHARACTER (LEN=*) :: var !< 448 449 INTEGER(iwp), OPTIONAL, INTENT(IN) :: i !< Current element of data_output 450 INTEGER(iwp), OPTIONAL, INTENT(IN) :: ilen !< Length of current entry in data_output 451 INTEGER(iwp), OPTIONAL, INTENT(IN) :: k !< Output is xy mode? 0 = no, 1 = yes 397 452 398 453 SELECT CASE ( TRIM( var ) ) … … 409 464 CASE ( 'ww' ) 410 465 unit = 'm2/s2' 411 466 ! 467 !-- Treat horizotal cross-section output quanatities 468 CASE ( 'theta_2m*', 'wspeed_10m*' ) 469 ! 470 !-- Check if output quantity is _xy only. 471 IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN 472 message_string = 'illegal value for data_output: "' // & 473 TRIM( var ) // '" & only 2d-horizontal ' // & 474 'cross sections are allowed for this value' 475 CALL message( 'diagnostic_output', 'PA0111', 1, 2, 0, 6, 0 ) 476 ENDIF 477 478 IF ( TRIM( var ) == 'theta_2m*' ) unit = 'K' 479 IF ( TRIM( var ) == 'wspeed_10m*' ) unit = 'm/s' 480 412 481 CASE DEFAULT 413 482 unit = 'illegal' … … 443 512 grid_x = 'x' 444 513 grid_y = 'y' 445 grid_z = 'zu' 514 grid_z = 'zu' 515 ! 516 !-- s grid surface variables 517 CASE ( 'theta_2m*_xy', 'wspeed_10m*' ) 518 519 grid_x = 'x' 520 grid_y = 'y' 521 grid_z = 'zu' 446 522 ! 447 523 !-- u grid … … 572 648 573 649 IF ( mode == 'xy' ) grid = 'zw' 650 651 CASE ( 'theta_2m*_xy' ) ! 2d-array 652 IF ( av == 0 ) THEN 653 DO i = nxl, nxr 654 DO j = nys, nyn 655 local_pf(i,j,nzb+1) = pt_2m(j,i) 656 ENDDO 657 ENDDO 658 ELSE 659 IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN 660 ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) ) 661 pt_2m_av = REAL( fill_value, KIND = wp ) 662 ENDIF 663 DO i = nxl, nxr 664 DO j = nys, nyn 665 local_pf(i,j,nzb+1) = pt_2m_av(j,i) 666 ENDDO 667 ENDDO 668 ENDIF 669 resorted = .TRUE. 670 two_d = .TRUE. 671 grid = 'zu1' 672 673 CASE ( 'wspeed_10m*_xy' ) ! 2d-array 674 IF ( av == 0 ) THEN 675 DO i = nxl, nxr 676 DO j = nys, nyn 677 local_pf(i,j,nzb+1) = uv_10m(j,i) 678 ENDDO 679 ENDDO 680 ELSE 681 IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN 682 ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) ) 683 uv_10m_av = REAL( fill_value, KIND = wp ) 684 ENDIF 685 DO i = nxl, nxr 686 DO j = nys, nyn 687 local_pf(i,j,nzb+1) = uv_10m_av(j,i) 688 ENDDO 689 ENDDO 690 ENDIF 691 resorted = .TRUE. 692 two_d = .TRUE. 693 grid = 'zu1' 574 694 575 695 CASE DEFAULT … … 776 896 flag_nr = 3 777 897 778 779 898 CASE DEFAULT 780 899 found = .FALSE. … … 864 983 uu = 0.0_wp 865 984 ENDIF 866 867 985 ! 868 986 !-- Allocate array for vv … … 872 990 vv = 0.0_wp 873 991 ENDIF 874 875 992 ! 876 993 !-- Allocate array for ww … … 879 996 ALLOCATE( ww(nzb:nzt+1,nys:nyn,nxl:nxr) ) 880 997 ww = 0.0_wp 998 ENDIF 999 ! 1000 !-- Allocate array for 2-m potential temperature 1001 CASE ( 'theta_2m*' ) 1002 IF ( .NOT. ALLOCATED( pt_2m ) ) THEN 1003 ALLOCATE( pt_2m(nys:nyn,nxl:nxr) ) 1004 pt_2m = 0.0_wp 1005 ENDIF 1006 ! 1007 !-- Allocate array for 10-m wind speed 1008 CASE ( 'wspeed_10m*' ) 1009 IF ( .NOT. ALLOCATED( uv_10m ) ) THEN 1010 ALLOCATE( uv_10m(nys:nyn,nxl:nxr) ) 1011 uv_10m = 0.0_wp 881 1012 ENDIF 882 1013 … … 900 1031 IMPLICIT NONE 901 1032 902 INTEGER(iwp) :: i, j, k !< grid loop index in x-, y-, z-direction 1033 INTEGER(iwp) :: i !< grid index x-dimension 1034 INTEGER(iwp) :: j !< grid index y-dimension 1035 INTEGER(iwp) :: k !< grid index z-dimension 903 1036 INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities 1037 1038 TYPE(surf_type), POINTER :: surf !< surf-type array, used to generalize subroutines 904 1039 905 1040 … … 973 1108 ENDDO 974 1109 ENDDO 1110 ! 1111 !-- 2-m potential temperature 1112 CASE ( 'theta_2m*' ) 1113 ! 1114 !-- 2-m potential temperature is caluclated from surface arrays. In 1115 !-- case the 2m level is below the first grid point, MOST is applied, 1116 !-- else, linear interpolation between two vertical grid levels is 1117 !-- applied. To access all surfaces, iterate over all horizontally- 1118 !-- upward facing surface types. 1119 surf => surf_def_h(0) 1120 CALL calc_pt_2m 1121 surf => surf_lsm_h 1122 CALL calc_pt_2m 1123 surf => surf_usm_h 1124 CALL calc_pt_2m 1125 ! 1126 !-- 10-m wind speed 1127 CASE ( 'wspeed_10m*' ) 1128 ! 1129 !-- 10-m wind speed is caluclated from surface arrays. In 1130 !-- case the 10m level is below the first grid point, MOST is applied, 1131 !-- else, linear interpolation between two vertical grid levels is 1132 !-- applied. To access all surfaces, iterate over all horizontally- 1133 !-- upward facing surface types. 1134 surf => surf_def_h(0) 1135 CALL calc_wind_10m 1136 surf => surf_lsm_h 1137 CALL calc_wind_10m 1138 surf => surf_usm_h 1139 CALL calc_wind_10m 975 1140 976 1141 END SELECT … … 981 1146 ! CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' ) 982 1147 983 END SUBROUTINE doq_calculate 984 985 1148 ! 1149 !-- The following block contains subroutines to calculate diagnostic 1150 !-- surface quantities. 1151 CONTAINS 986 1152 !------------------------------------------------------------------------------! 987 1153 ! Description: 988 1154 ! ------------ 989 !> ... 1155 !> Calculation of 2-m potential temperature. 1156 !------------------------------------------------------------------------------! 1157 SUBROUTINE calc_pt_2m 1158 1159 USE surface_layer_fluxes_mod, & 1160 ONLY: psi_h 1161 1162 IMPLICIT NONE 1163 1164 INTEGER(iwp) :: kk !< running index along the z-dimension 1165 INTEGER(iwp) :: m !< running index for surface elements 1166 1167 DO m = 1, surf%ns 1168 1169 i = surf%i(m) 1170 j = surf%j(m) 1171 k = surf%k(m) 1172 ! 1173 !-- If 2-m level is below the first grid level, MOST is 1174 !-- used for calculation of 2-m temperature. 1175 IF ( surf%z_mo(m) > 2.0_wp ) THEN 1176 pt_2m(j,i) = surf%pt_surface(m) + surf%ts(m) / kappa & 1177 * ( LOG( 2.0_wp / surf%z0h(m) ) & 1178 - psi_h( 2.0_wp / surf%ol(m) ) & 1179 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1180 ! 1181 !-- If 2-m level is above the first grid level, 2-m temperature 1182 !-- is linearly interpolated between the two nearest vertical grid 1183 !-- levels. Note, since 2-m temperature is only computed for 1184 !-- horizontal upward-facing surfaces, only a vertical 1185 !-- interpolation is necessary. 1186 ELSE 1187 ! 1188 !-- zw(k-1) defines the height of the surface. 1189 kk = k 1190 DO WHILE ( zu(kk) - zw(k-1) < 2.0_wp .AND. kk <= nzt ) 1191 kk = kk + 1 1192 ENDDO 1193 ! 1194 !-- kk defines the index of the first grid level >= 2m. 1195 pt_2m(j,i) = pt(kk-1,j,i) + & 1196 ( zw(k-1) + 2.0_wp - zu(kk-1) ) * & 1197 ( pt(kk,j,i) - pt(kk-1,j,i) ) / & 1198 ( zu(kk) - zu(kk-1) ) 1199 ENDIF 1200 1201 ENDDO 1202 1203 END SUBROUTINE calc_pt_2m 1204 1205 !------------------------------------------------------------------------------! 1206 ! Description: 1207 ! ------------ 1208 !> Calculation of 10-m wind speed. 1209 !------------------------------------------------------------------------------! 1210 SUBROUTINE calc_wind_10m 1211 1212 USE surface_layer_fluxes_mod, & 1213 ONLY: psi_m 1214 1215 IMPLICIT NONE 1216 1217 INTEGER(iwp) :: kk !< running index along the z-dimension 1218 INTEGER(iwp) :: m !< running index for surface elements 1219 1220 REAL(wp) :: uv_l !< wind speed at lower grid point 1221 REAL(wp) :: uv_u !< wind speed at upper grid point 1222 1223 DO m = 1, surf%ns 1224 1225 i = surf%i(m) 1226 j = surf%j(m) 1227 k = surf%k(m) 1228 ! 1229 !-- If 10-m level is below the first grid level, MOST is 1230 !-- used for calculation of 10-m temperature. 1231 IF ( surf%z_mo(m) > 10.0_wp ) THEN 1232 uv_10m(j,i) = surf%us(m) / kappa & 1233 * ( LOG( 10.0_wp / surf%z0(m) ) & 1234 - psi_m( 10.0_wp / surf%ol(m) ) & 1235 + psi_m( surf%z0(m) / surf%ol(m) ) ) 1236 ! 1237 !-- If 10-m level is above the first grid level, 10-m wind speed 1238 !-- is linearly interpolated between the two nearest vertical grid 1239 !-- levels. Note, since 10-m temperature is only computed for 1240 !-- horizontal upward-facing surfaces, only a vertical 1241 !-- interpolation is necessary. 1242 ELSE 1243 ! 1244 !-- zw(k-1) defines the height of the surface. 1245 kk = k 1246 DO WHILE ( zu(kk) - zw(k-1) < 10.0_wp .AND. kk <= nzt ) 1247 kk = kk + 1 1248 ENDDO 1249 ! 1250 !-- kk defines the index of the first grid level >= 10m. 1251 uv_l = SQRT( ( 0.5_wp * ( u(kk-1,j,i) + u(kk-1,j,i+1) ) )**2 & 1252 + ( 0.5_wp * ( v(kk-1,j,i) + v(kk-1,j+1,i) ) )**2 ) 1253 1254 uv_u = SQRT( ( 0.5_wp * ( u(kk,j,i) + u(kk,j,i+1) ) )**2 & 1255 + ( 0.5_wp * ( v(kk,j,i) + v(kk,j+1,i) ) )**2 ) 1256 1257 uv_10m(j,i) = uv_l + ( zw(k-1) + 10.0_wp - zu(kk-1) ) * & 1258 ( uv_u - uv_l ) / & 1259 ( zu(kk) - zu(kk-1) ) 1260 1261 ENDIF 1262 1263 ENDDO 1264 1265 END SUBROUTINE calc_wind_10m 1266 1267 END SUBROUTINE doq_calculate 1268 1269 1270 !------------------------------------------------------------------------------! 1271 ! Description: 1272 ! ------------ 1273 !> Preparation of the diagnostic output, counting of the module-specific 1274 !> output quantities and gathering of the output names. 990 1275 !------------------------------------------------------------------------------! 991 1276 SUBROUTINE doq_prepare … … 1173 1458 IMPLICIT NONE 1174 1459 1460 IF ( ALLOCATED( pt_2m_av ) ) THEN 1461 CALL wrd_write_string( 'pt_2m_av' ) 1462 WRITE ( 14 ) pt_2m_av 1463 ENDIF 1464 1175 1465 IF ( ALLOCATED( ti_av ) ) THEN 1176 CALL wrd_write_string( 'ti_av' ) 1466 CALL wrd_write_string( 'ti_av' ) 1177 1467 WRITE ( 14 ) ti_av 1178 1468 ENDIF 1179 1469 1180 1470 IF ( ALLOCATED( uu_av ) ) THEN 1181 CALL wrd_write_string( 'uu_av' ) 1471 CALL wrd_write_string( 'uu_av' ) 1182 1472 WRITE ( 14 ) uu_av 1183 1473 ENDIF 1184 1474 1475 IF ( ALLOCATED( uv_10m_av ) ) THEN 1476 CALL wrd_write_string( 'uv_10m_av' ) 1477 WRITE ( 14 ) uv_10m_av 1478 ENDIF 1479 1185 1480 IF ( ALLOCATED( vv_av ) ) THEN 1186 CALL wrd_write_string( 'vv_av' ) 1481 CALL wrd_write_string( 'vv_av' ) 1187 1482 WRITE ( 14 ) vv_av 1188 1483 ENDIF 1189 1484 1190 1485 IF ( ALLOCATED( ww_av ) ) THEN 1191 CALL wrd_write_string( 'ww_av' ) 1486 CALL wrd_write_string( 'ww_av' ) 1192 1487 WRITE ( 14 ) ww_av 1193 1488 ENDIF -
palm/trunk/SOURCE/module_interface.f90
r4281 r4331 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Change interface for doq_check_data_output, in order to perform further 28 ! output checks. 29 ! 30 ! 4281 2019-10-29 15:15:39Z schwenkel 27 31 ! Added dynamics boundary conditions 28 32 ! … … 872 876 873 877 IF ( unit == 'illegal' ) THEN 874 CALL doq_check_data_output( variable, unit )878 CALL doq_check_data_output( variable, unit, i, ilen, k ) 875 879 ENDIF 876 880 -
palm/trunk/SOURCE/modules.f90
r4329 r4331 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - do_output_at_2m, pt_2m_av 28 ! 29 ! 4329 2019-12-10 15:46:36Z motisi 27 30 ! Renamed wall_flags_0 to wall_flags_static_0 28 31 ! … … 413 416 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: lwp_av !< avg. liquid water path 414 417 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ol_av !< avg. Obukhov length 415 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt_2m_av !< avg. 2m- air potential temperature416 418 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: qsws_av !< avg. surface moisture flux 417 419 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: r_a_av !< avg. resistance … … 707 709 LOGICAL :: do2d_at_begin = .FALSE. !< namelist parameter 708 710 LOGICAL :: do3d_at_begin = .FALSE. !< namelist parameter 709 LOGICAL :: do_output_at_2m = .FALSE. !< flag for activating calculation of potential temperature at z = 2 m710 711 LOGICAL :: do_sum = .FALSE. !< contribute to time average of profile data? 711 712 LOGICAL :: dp_external = .FALSE. !< namelist parameter -
palm/trunk/SOURCE/plant_canopy_model_mod.f90
r4329 r4331 27 27 ! ----------------- 28 28 ! $Id$ 29 ! Typo corrected 30 ! 31 ! 4329 2019-12-10 15:46:36Z motisi 29 32 ! Renamed wall_flags_0 to wall_flags_static_0 30 33 ! … … 1135 1138 'defined on top of an artificially '// & 1136 1139 'created building grid point ' // & 1137 '(emerged from the fil ering) - ' //&1140 '(emerged from the filtering) - ' // & 1138 1141 'LAD profile is omitted at this ' // & 1139 1142 'grid point: (i,j) = ', i, j -
palm/trunk/SOURCE/read_restart_data_mod.f90
r4301 r4331 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Enable restart data for 2-m potential temperature output 28 ! 29 ! 4301 2019-11-22 12:09:09Z oliver.maas 27 30 ! removed recycling_yshift 28 31 ! … … 101 104 102 105 USE diagnostic_output_quantities_mod, & 103 ONLY: ti_av, uu_av, vv_av, ww_av 106 ONLY: pt_2m_av, & 107 ti_av, & 108 uu_av, & 109 uv_10m_av, & 110 vv_av, & 111 ww_av 104 112 105 113 USE grid_variables, & … … 1525 1533 uu_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf) 1526 1534 1535 CASE ( 'uv_10m_av' ) 1536 IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN 1537 ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) ) 1538 ENDIF 1539 IF ( k == 1 ) READ ( 13 ) tmp_2d 1540 uv_10m_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 1541 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 1542 1527 1543 CASE ( 'u_m_l' ) 1528 1544 IF ( k == 1 ) THEN -
palm/trunk/SOURCE/sum_up_3d_data.f90
r4182 r4331 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Move 2-m potential temperature output to diagnostic_output_quantities 28 ! 29 ! 4182 2019-08-22 15:20:23Z scharf 27 30 ! Corrected "Former revisions" section 28 31 ! … … 64 67 65 68 USE arrays_3d, & 66 ONLY: dzw, d_exner, e, heatflux_output_conversion, p, &67 pt, q, ql, ql_c, ql_v, s, u, v, vpt, w, &69 ONLY: dzw, d_exner, e, heatflux_output_conversion, p, & 70 pt, q, ql, ql_c, ql_v, s, u, v, vpt, w, & 68 71 waterflux_output_conversion 69 72 70 73 USE averaging, & 71 74 ONLY: e_av, ghf_av, lpt_av, lwp_av, ol_av, p_av, pc_av, pr_av, pt_av, & 72 pt_2m_av, q_av, ql_av, ql_c_av, ql_v_av, ql_vp_av, qsws_av,&75 q_av, ql_av, ql_c_av, ql_v_av, ql_vp_av, qsws_av, & 73 76 qv_av, r_a_av, s_av, shf_av, ssws_av, ts_av, tsurf_av, u_av, & 74 77 us_av, v_av, vpt_av, w_av, z0_av, z0h_av, z0q_av … … 296 299 ENDIF 297 300 vpt_av = 0.0_wp 298 299 CASE ( 'theta_2m*' )300 IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN301 ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )302 ENDIF303 pt_2m_av = 0.0_wp304 301 305 302 CASE ( 'w' ) … … 730 727 ENDIF 731 728 732 CASE ( 'theta_2m*' )733 IF ( ALLOCATED( pt_2m_av ) ) THEN734 DO i = nxl, nxr735 DO j = nys, nyn736 match_def = surf_def_h(0)%start_index(j,i) <= &737 surf_def_h(0)%end_index(j,i)738 match_lsm = surf_lsm_h%start_index(j,i) <= &739 surf_lsm_h%end_index(j,i)740 match_usm = surf_usm_h%start_index(j,i) <= &741 surf_usm_h%end_index(j,i)742 743 IF ( match_def ) THEN744 m = surf_def_h(0)%end_index(j,i)745 pt_2m_av(j,i) = pt_2m_av(j,i) + &746 surf_def_h(0)%pt_2m(m)747 ELSEIF ( match_lsm .AND. .NOT. match_usm ) THEN748 m = surf_lsm_h%end_index(j,i)749 pt_2m_av(j,i) = pt_2m_av(j,i) + &750 surf_lsm_h%pt_2m(m)751 ELSEIF ( match_usm ) THEN752 m = surf_usm_h%end_index(j,i)753 pt_2m_av(j,i) = pt_2m_av(j,i) + &754 surf_usm_h%pt_2m(m)755 ENDIF756 ENDDO757 ENDDO758 ENDIF759 760 761 729 CASE ( 't*' ) 762 730 IF ( ALLOCATED( ts_av ) ) THEN -
palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
r4298 r4331 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Calculation of diagnostic-only 2-m potential temperature moved to 29 ! diagnostic_output_quantities. 30 ! 31 ! 4298 2019-11-21 15:59:16Z suehring 28 32 ! Calculation of 2-m temperature adjusted to the case the 2-m level is above 29 33 ! the first grid point. … … 108 112 constant_waterflux, coupling_mode, & 109 113 debug_output_timestep, & 110 do_output_at_2m, humidity,&114 humidity, & 111 115 ibc_e_b, ibc_pt_b, indoor_model, & 112 116 land_surface, large_scale_forcing, lsf_surf, message_string, & … … 162 166 PRIVATE 163 167 164 PUBLIC init_surface_layer_fluxes, phi_m, surface_layer_fluxes 168 PUBLIC init_surface_layer_fluxes, & 169 phi_m, & 170 psi_h, & 171 psi_m, & 172 surface_layer_fluxes 165 173 166 174 INTERFACE init_surface_layer_fluxes … … 171 179 MODULE PROCEDURE phi_m 172 180 END INTERFACE phi_m 181 182 INTERFACE psi_h 183 MODULE PROCEDURE psi_h 184 END INTERFACE psi_h 185 186 INTERFACE psi_m 187 MODULE PROCEDURE psi_m 188 END INTERFACE psi_m 173 189 174 190 INTERFACE surface_layer_fluxes … … 259 275 CALL calc_scaling_parameters 260 276 CALL calc_surface_fluxes 261 IF ( do_output_at_2m ) THEN262 CALL calc_pt_near_surface ( '2m' )263 ENDIF264 277 ENDIF 265 278 ! … … 272 285 CALL calc_scaling_parameters 273 286 CALL calc_surface_fluxes 274 IF ( do_output_at_2m ) THEN275 CALL calc_pt_near_surface ( '2m' )276 ENDIF277 287 ENDIF 278 288 ! … … 285 295 CALL calc_scaling_parameters 286 296 CALL calc_surface_fluxes 287 IF ( do_output_at_2m ) THEN288 CALL calc_pt_near_surface ( '2m' )289 ENDIF290 297 ! 291 298 !-- Calculate 10cm temperature, required in indoor model … … 1863 1870 INTEGER(iwp) :: j !< grid index y-dimension 1864 1871 INTEGER(iwp) :: k !< grid index z-dimension 1865 INTEGER(iwp) :: kk !< running index along the z-dimension1866 1872 INTEGER(iwp) :: m !< running index for surface elements 1867 1873 … … 1882 1888 - psi_h( 0.1_wp / surf%ol(m) ) & 1883 1889 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1884 1885 ENDDO 1886 1887 ! 1888 !-- 2-m temperature. Note, this is only calculated for output reasons at 1889 !-- horizontal upward-facing surfaces. 1890 CASE ( '2m' ) 1891 1892 DO m = 1, surf%ns 1893 1894 i = surf%i(m) 1895 j = surf%j(m) 1896 k = surf%k(m) 1897 ! 1898 !-- If 2-m level is below the first grid level, MOST is 1899 !-- used for calculation of 2-m temperature. 1900 IF ( surf%z_mo(m) > 2.0_wp ) THEN 1901 surf%pt_2m(m) = surf%pt_surface(m) + surf%ts(m) / kappa & 1902 * ( LOG( 2.0_wp / surf%z0h(m) ) & 1903 - psi_h( 2.0_wp / surf%ol(m) ) & 1904 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1905 ! 1906 !-- If 2-m level is above the first grid level, 2-m temperature 1907 !-- is linearly interpolated between the two nearest vertical grid 1908 !-- levels. Note, since 2-m temperature is only computed for 1909 !-- horizontal upward-facing surfaces, only a vertical 1910 !-- interpolation is necessary. 1911 ELSE 1912 ! 1913 !-- zw(k-1) defines the height of the surface. 1914 kk = k 1915 DO WHILE ( zu(kk) - zw(k-1) < 2.0_wp .AND. kk <= nzt ) 1916 kk = kk + 1 1917 ENDDO 1918 ! 1919 !-- kk defines the index of the first grid level >= 2m. 1920 surf%pt_2m(m) = pt(kk-1,j,i) + & 1921 ( zw(k-1) + 2.0_wp - zu(kk-1) ) * & 1922 ( pt(kk,j,i) - pt(kk-1,j,i) ) / & 1923 ( zu(kk) - zu(kk-1) ) 1924 ENDIF 1925 1926 ENDDO 1927 1928 1890 1891 ENDDO 1892 1929 1893 END SELECT 1930 1894 -
palm/trunk/SOURCE/surface_mod.f90
r4329 r4331 26 26 ! ----------------- 27 27 ! $Id$ 28 ! -pt_2m - array is moved to diagnostic_output_quantities 29 ! 30 ! 4329 2019-12-10 15:46:36Z motisi 28 31 ! Renamed wall_flags_0 to wall_flags_static_0 29 32 ! … … 289 292 290 293 REAL(wp), DIMENSION(:), ALLOCATABLE :: pt_10cm !< near surface air potential temperature at distance 10 cm from the surface (K) 291 REAL(wp), DIMENSION(:), ALLOCATABLE :: pt_2m !< near surface air potential temperature at distance 2 m from the surface (K)292 294 293 295 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: alpha_vg !< coef. of Van Genuchten … … 1281 1283 !-- Salinity surface flux 1282 1284 IF ( ocean_mode ) DEALLOCATE ( surfaces%sasws ) 1283 !1284 !-- 2-m potential temperature (for output quantity theta_2m*)1285 IF ( do_output_at_2m ) DEALLOCATE ( surfaces%pt_2m )1286 1285 1287 1286 END SUBROUTINE deallocate_surface_attributes_h … … 1412 1411 !-- Salinity surface flux 1413 1412 IF ( ocean_mode ) ALLOCATE ( surfaces%sasws(1:surfaces%ns) ) 1414 !1415 !-- 2-m potential temperature (for output quantity theta_2m*)1416 IF ( do_output_at_2m ) THEN1417 ALLOCATE ( surfaces%pt_2m(1:surfaces%ns) )1418 surfaces%pt_2m = -9999.0_wp !< output array (for theta_2m*) must be initialized here,1419 !< otherwise simulation crash at do2d_at_begin with spinup=.F.1420 ENDIF1421 1413 1422 1414 END SUBROUTINE allocate_surface_attributes_h -
palm/trunk/SOURCE/time_integration_spinup.f90
r4227 r4331 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Enable output of diagnostic quantities, e.g. 2-m temperature 28 ! 29 ! 4227 2019-09-10 18:04:34Z gronemeier 27 30 ! implement new palm_date_time_mod 28 31 ! … … 81 84 USE cpulog, & 82 85 ONLY: cpu_log, log_point_s 86 87 USE diagnostic_output_quantities_mod, & 88 ONLY: doq_calculate 83 89 84 90 USE indices, & … … 468 474 !-- 2d-data output (cross-sections) 469 475 IF ( time_do2d_xy >= dt_do2d_xy ) THEN 476 CALL doq_calculate 470 477 CALL data_output_2d( 'xy', 0 ) 471 478 time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) ) … … 475 482 !-- 3d-data output (volume data) 476 483 IF ( time_do3d >= dt_do3d ) THEN 484 CALL doq_calculate 477 485 CALL data_output_3d( 0 ) 478 486 time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) ) … … 560 568 561 569 USE basic_constants_and_equations_mod, & 562 ONLY: pi570 ONLY: pi 563 571 564 572 USE kinds -
palm/trunk/SOURCE/write_restart_data_mod.f90
r4301 r4331 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Move 2-m potential temperature to diagnostic_output_quantities. 28 ! 29 ! 4301 2019-11-22 12:09:09Z oliver.maas 27 30 ! removed recycling_yshift 28 31 ! … … 942 945 ENDIF 943 946 944 IF ( ALLOCATED( pt_2m_av ) ) THEN945 CALL wrd_write_string( 'pt_2m_av' )946 WRITE ( 14 ) pt_2m_av947 ENDIF948 949 947 IF ( humidity ) THEN 950 948
Note: See TracChangeset
for help on using the changeset viewer.