Changeset 4581 for palm/trunk/SOURCE/chemistry_model_mod.f90
- Timestamp:
- Jun 29, 2020 8:49:58 AM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/chemistry_model_mod.f90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/chemistry_model_mod.f90
r4577 r4581 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Enable output of resolved-scale vertical fluxes of chemical species. 29 ! 30 ! 4577 2020-06-25 09:53:58Z raasch 28 31 ! further re-formatting to follow the PALM coding standard 29 32 ! … … 436 439 ! 437 440 !-- Interface section 438 INTERFACE chem_3d_data_averaging439 MODULE PROCEDURE chem_3d_data_averaging440 END INTERFACE chem_3d_data_averaging441 442 INTERFACE chem_boundary_conditions443 MODULE PROCEDURE chem_boundary_conditions444 END INTERFACE chem_boundary_conditions445 446 INTERFACE chem_check_data_output447 MODULE PROCEDURE chem_check_data_output448 END INTERFACE chem_check_data_output449 450 INTERFACE chem_data_output_2d451 MODULE PROCEDURE chem_data_output_2d452 END INTERFACE chem_data_output_2d453 454 INTERFACE chem_data_output_3d455 MODULE PROCEDURE chem_data_output_3d456 END INTERFACE chem_data_output_3d457 458 INTERFACE chem_data_output_mask459 MODULE PROCEDURE chem_data_output_mask460 END INTERFACE chem_data_output_mask461 462 INTERFACE chem_check_data_output_pr463 MODULE PROCEDURE chem_check_data_output_pr464 END INTERFACE chem_check_data_output_pr465 466 INTERFACE chem_check_parameters467 MODULE PROCEDURE chem_check_parameters468 END INTERFACE chem_check_parameters469 470 INTERFACE chem_define_netcdf_grid471 MODULE PROCEDURE chem_define_netcdf_grid472 END INTERFACE chem_define_netcdf_grid473 474 INTERFACE chem_header475 MODULE PROCEDURE chem_header476 END INTERFACE chem_header477 478 INTERFACE chem_init_arrays479 MODULE PROCEDURE chem_init_arrays480 END INTERFACE chem_init_arrays481 482 INTERFACE chem_init483 MODULE PROCEDURE chem_init484 END INTERFACE chem_init485 486 INTERFACE chem_init_profiles487 MODULE PROCEDURE chem_init_profiles488 END INTERFACE chem_init_profiles489 490 INTERFACE chem_integrate491 MODULE PROCEDURE chem_integrate_ij492 END INTERFACE chem_integrate493 494 INTERFACE chem_parin495 MODULE PROCEDURE chem_parin496 END INTERFACE chem_parin497 498 441 INTERFACE chem_actions 499 442 MODULE PROCEDURE chem_actions 500 443 MODULE PROCEDURE chem_actions_ij 501 444 END INTERFACE chem_actions 445 446 INTERFACE chem_3d_data_averaging 447 MODULE PROCEDURE chem_3d_data_averaging 448 END INTERFACE chem_3d_data_averaging 449 450 INTERFACE chem_boundary_conditions 451 MODULE PROCEDURE chem_boundary_conditions 452 END INTERFACE chem_boundary_conditions 453 454 INTERFACE chem_check_data_output 455 MODULE PROCEDURE chem_check_data_output 456 END INTERFACE chem_check_data_output 457 458 INTERFACE chem_data_output_2d 459 MODULE PROCEDURE chem_data_output_2d 460 END INTERFACE chem_data_output_2d 461 462 INTERFACE chem_data_output_3d 463 MODULE PROCEDURE chem_data_output_3d 464 END INTERFACE chem_data_output_3d 465 466 INTERFACE chem_data_output_mask 467 MODULE PROCEDURE chem_data_output_mask 468 END INTERFACE chem_data_output_mask 469 470 INTERFACE chem_check_data_output_pr 471 MODULE PROCEDURE chem_check_data_output_pr 472 END INTERFACE chem_check_data_output_pr 473 474 INTERFACE chem_check_parameters 475 MODULE PROCEDURE chem_check_parameters 476 END INTERFACE chem_check_parameters 477 478 INTERFACE chem_define_netcdf_grid 479 MODULE PROCEDURE chem_define_netcdf_grid 480 END INTERFACE chem_define_netcdf_grid 481 482 INTERFACE chem_header 483 MODULE PROCEDURE chem_header 484 END INTERFACE chem_header 485 486 INTERFACE chem_init_arrays 487 MODULE PROCEDURE chem_init_arrays 488 END INTERFACE chem_init_arrays 489 490 INTERFACE chem_init 491 MODULE PROCEDURE chem_init 492 END INTERFACE chem_init 493 494 INTERFACE chem_init_profiles 495 MODULE PROCEDURE chem_init_profiles 496 END INTERFACE chem_init_profiles 497 498 INTERFACE chem_integrate 499 MODULE PROCEDURE chem_integrate_ij 500 END INTERFACE chem_integrate 501 502 INTERFACE chem_parin 503 MODULE PROCEDURE chem_parin 504 END INTERFACE chem_parin 502 505 503 506 INTERFACE chem_non_advective_processes … … 1069 1072 1070 1073 1071 CHARACTER (LEN=*) :: unit !< 1072 CHARACTER (LEN=*) :: variable !< 1073 CHARACTER (LEN=*) :: dopr_unit 1074 CHARACTER (LEN=16) :: spec_name 1075 1076 INTEGER(iwp) :: var_count, lsp !< 1077 1078 1079 spec_name = TRIM( variable(4:) ) 1080 1081 IF ( .NOT. air_chemistry ) THEN 1082 message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // & 1083 ' is not imp' // 'lemented for air_chemistry = .FALSE.' 1084 CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 ) 1085 1086 ELSE 1074 CHARACTER (LEN=*) :: unit !< unit 1075 CHARACTER (LEN=*) :: variable !< variable name 1076 CHARACTER (LEN=*) :: dopr_unit !< unit 1077 CHARACTER (LEN=16) :: spec_name !< species name extracted from output string 1078 1079 INTEGER(iwp) :: index_start !< start index of the species name in data-output string 1080 INTEGER(iwp) :: index_end !< end index of the species name in data-output string 1081 INTEGER(iwp) :: var_count !< number of data-output quantity 1082 INTEGER(iwp) :: lsp !< running index over species 1083 1084 SELECT CASE ( TRIM( variable(1:3 ) ) ) 1085 1086 CASE ( 'kc_' ) 1087 1088 IF ( .NOT. air_chemistry ) THEN 1089 message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // & 1090 ' is not implemented for air_chemistry = .FALSE.' 1091 CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 ) 1092 1093 ENDIF 1094 ! 1095 !-- Output of total fluxes is not allowed to date. 1096 IF ( TRIM( variable(1:4) ) == 'kc_w' ) THEN 1097 IF ( TRIM( variable(1:5) ) /= 'kc_w*' .AND. TRIM( variable(1:5) ) /= 'kc_w"' ) THEN 1098 message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // & 1099 ' is currently not implemented. Please output resolved- and '// & 1100 'subgrid-scale fluxes individually to obtain the total flux.' 1101 CALL message( 'chem_check_parameters', 'CM0498', 1, 2, 0, 6, 0 ) 1102 ENDIF 1103 ENDIF 1104 ! 1105 !-- Check for profile output of first-order moments, i.e. variable(4:) equals a species name. 1106 spec_name = TRIM( variable(4:) ) 1087 1107 DO lsp = 1, nspec 1088 1108 IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN 1089 cs_pr_count = cs_pr_count+11090 cs_pr_index (cs_pr_count) = lsp1091 dopr_index(var_count) = pr_palm+cs_pr_count1092 dopr_unit = 'ppm'1109 cs_pr_count_sp = cs_pr_count_sp + 1 1110 cs_pr_index_sp(cs_pr_count_sp) = lsp 1111 dopr_index(var_count) = pr_palm + cs_pr_count_sp + cs_pr_count_fl_sgs + cs_pr_count_fl_res 1112 dopr_unit = 'ppm' 1093 1113 IF ( spec_name(1:2) == 'PM') THEN 1094 1114 dopr_unit = 'kg m-3' 1095 1115 ENDIF 1096 hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 ) 1097 unit = dopr_unit 1116 hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 ) 1117 unit = dopr_unit 1118 1119 hom_index_spec(cs_pr_count_sp) = dopr_index(var_count) 1098 1120 ENDIF 1099 1121 ENDDO 1100 ENDIF 1122 ! 1123 !-- Check for profile output of fluxes. variable(index_start:index_end) equals a species name. 1124 !-- Start with SGS components. 1125 IF ( TRIM( variable(1:5) ) == 'kc_w"' ) THEN 1126 DO lsp = 1, nspec 1127 index_start = INDEX( TRIM( variable ), TRIM( chem_species(lsp)%name ) ) 1128 index_end = LEN( TRIM( variable ) ) - 1 1129 IF ( index_start /= 0 ) THEN 1130 spec_name = TRIM( variable(index_start:index_end) ) 1131 IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN 1132 cs_pr_count_fl_sgs = cs_pr_count_fl_sgs + 1 1133 cs_pr_index_fl_sgs(cs_pr_count_fl_sgs) = lsp 1134 dopr_index(var_count) = pr_palm + cs_pr_count_sp + & 1135 cs_pr_count_fl_sgs + & 1136 cs_pr_count_fl_res 1137 dopr_unit = 'm ppm s-1' 1138 IF ( spec_name(1:2) == 'PM') THEN 1139 dopr_unit = 'kg m-2 s-1' 1140 ENDIF 1141 hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 ) 1142 unit = dopr_unit 1143 1144 hom_index_fl_sgs(cs_pr_count_fl_sgs) = dopr_index(var_count) 1145 ENDIF 1146 ENDIF 1147 ENDDO 1148 ENDIF 1149 ! 1150 !-- Proceed with resolved-scale fluxes. 1151 IF ( TRIM( variable(1:5) ) == 'kc_w*' ) THEN 1152 spec_name = TRIM( variable(6:) ) 1153 DO lsp = 1, nspec 1154 index_start = INDEX( TRIM( variable ), TRIM( chem_species(lsp)%name ) ) 1155 index_end = LEN( TRIM( variable ) ) - 1 1156 IF ( index_start /= 0 ) THEN 1157 spec_name = TRIM( variable(index_start:index_end) ) 1158 IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN 1159 cs_pr_count_fl_res = cs_pr_count_fl_res + 1 1160 cs_pr_index_fl_res(cs_pr_count_fl_res) = lsp 1161 dopr_index(var_count) = pr_palm + cs_pr_count_sp + & 1162 cs_pr_count_fl_sgs + & 1163 cs_pr_count_fl_res 1164 dopr_unit = 'm ppm s-1' 1165 IF ( spec_name(1:2) == 'PM') THEN 1166 dopr_unit = 'kg m-2 s-1' 1167 ENDIF 1168 hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 ) 1169 unit = dopr_unit 1170 1171 hom_index_fl_res(cs_pr_count_fl_res) = dopr_index(var_count) 1172 ENDIF 1173 ENDIF 1174 ENDDO 1175 ENDIF 1176 CASE DEFAULT 1177 unit = 'illegal' 1178 1179 END SELECT 1101 1180 1102 1181 END SUBROUTINE chem_check_data_output_pr … … 1976 2055 spec_conc_3 (:,:,:,:) = 0.0_wp 1977 2056 2057 ! 2058 !-- Allocate array to store locally summed-up resolved-scale vertical fluxes. 2059 IF ( scalar_advec == 'ws-scheme' ) THEN 2060 ALLOCATE( sums_ws_l(nzb:nzt+1,0:threads_per_task-1,nspec) ) 2061 sums_ws_l = 0.0_wp 2062 ENDIF 1978 2063 1979 2064 DO lsp = 1, nspec … … 2070 2155 ENDDO 2071 2156 ENDIF 2157 2072 2158 ENDIF 2073 2159 ! … … 2620 2706 ! 2621 2707 !-- Determine the number of user-defined profiles and append them to the standard data output 2622 ! >>(data_output_pr)2708 !-- (data_output_pr) 2623 2709 max_pr_cs_tmp = 0 2624 2710 i = 1 2625 2711 2626 2712 DO WHILE ( data_output_pr(i) /= ' ' .AND. i <= 100 ) 2713 2627 2714 IF ( TRIM( data_output_pr(i)(1:3) ) == 'kc_' ) THEN 2628 2715 max_pr_cs_tmp = max_pr_cs_tmp+1 … … 2708 2795 ENDIF 2709 2796 2797 CASE ( 'before_timestep' ) 2798 ! 2799 !-- Set array used to sum-up resolved scale fluxes to zero. 2800 IF ( ws_scheme_sca ) THEN 2801 sums_ws_l = 0.0_wp 2802 ENDIF 2803 2710 2804 CASE DEFAULT 2711 2805 CONTINUE … … 2876 2970 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2877 2971 IF ( ws_scheme_sca ) THEN 2972 sums_wschs_ws_l(nzb:,0:) => sums_ws_l(:,:,ilsp) 2878 2973 CALL advec_s_ws( cs_advc_flags_s, chem_species(ilsp)%conc, 'kc', & 2879 2974 bc_dirichlet_cs_l .OR. bc_radiation_cs_l, & … … 2986 3081 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2987 3082 IF ( ws_scheme_sca ) THEN 3083 sums_wschs_ws_l(nzb:,0:) => sums_ws_l(nzb:nzt+1,0:threads_per_task-1,ilsp) 2988 3084 CALL advec_s_ws( cs_advc_flags_s, & 2989 3085 i, & … … 3029 3125 * ( chem_species(ilsp)%conc(k,j,i) - chem_species(ilsp)%conc_pr_init(k) ) & 3030 3126 ) & 3031 * MERGE( 1.0_wp, 0.0_wp, & 3032 BTEST( wall_flags_total_0(k,j,i), 0 ) & 3127 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) & 3033 3128 ) 3034 3129 … … 3170 3265 CHARACTER (LEN=*) :: mode !< 3171 3266 3172 INTEGER(iwp) :: i !< running index on x-axis 3173 INTEGER(iwp) :: j !< running index on y-axis 3174 INTEGER(iwp) :: k !< vertical index counter 3175 INTEGER(iwp) :: lpr !< running index chem spcs 3176 INTEGER(iwp) :: sr !< statistical region 3177 INTEGER(iwp) :: tn !< thread number 3267 INTEGER(iwp) :: i !< running index on x-axis 3268 INTEGER(iwp) :: j !< running index on y-axis 3269 INTEGER(iwp) :: k !< vertical index counter 3270 INTEGER(iwp) :: lpr !< running index chem spcs 3271 INTEGER(iwp) :: m !< running index for surface elements 3272 INTEGER(iwp) :: sr !< statistical region 3273 INTEGER(iwp) :: surf_e !< end surface index 3274 INTEGER(iwp) :: surf_s !< start surface index 3275 INTEGER(iwp) :: tn !< thread number 3276 3277 REAL(wp) :: flag !< topography masking flag 3278 REAL(wp), DIMENSION(nzb:nzt+1,0:threads_per_task-1) :: sums_tmp !< temporary array used to sum-up profiles 3178 3279 3179 3280 IF ( mode == 'profiles' ) THEN 3180 ! 3181 ! 3182 !-- Sample on how to calculate horizontally averaged profiles of user-defined quantities. Each 3183 !-- quantity is identified by the index "pr_palm+#" where "#" is an integer starting from 1. 3184 !-- These user-profile-numbers must also be assigned to the respective strings given by 3185 !-- data_output_pr_cs in routine user_check_data_output_pr. 3186 !-- hom(:,:,:,:) = dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g. 3187 !-- w*pt*), dim-4 = statistical region. 3188 3189 !$OMP DO 3190 DO i = nxl, nxr 3191 DO j = nys, nyn 3192 DO k = nzb, nzt+1 3193 DO lpr = 1, cs_pr_count 3194 3195 sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) + & 3196 chem_species(cs_pr_index(lpr))%conc(k,j,i) * & 3281 ! 3282 !-- Sum-up profiles for the species 3283 tn = 0 3284 !$OMP PARALLEL PRIVATE( i, j, k, tn, lpr, sums_tmp ) 3285 !$ tn = omp_get_thread_num() 3286 !$OMP DO 3287 DO lpr = 1, cs_pr_count_sp 3288 sums_tmp(:,tn) = 0.0_wp 3289 DO i = nxl, nxr 3290 DO j = nys, nyn 3291 DO k = nzb, nzt+1 3292 sums_tmp(k,tn) = sums_tmp(k,tn) + & 3293 chem_species(cs_pr_index_sp(lpr))%conc(k,j,i) * & 3197 3294 rmask(j,i,sr) * & 3198 MERGE( 1.0_wp, 0.0_wp, & 3199 BTEST( wall_flags_total_0(k,j,i), 22 ) ) 3295 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 22 ) ) 3200 3296 ENDDO 3201 3297 ENDDO 3202 3298 ENDDO 3299 sums_l(nzb:nzt+1,hom_index_spec(lpr),tn) = sums_tmp(nzb:nzt+1,tn) 3203 3300 ENDDO 3301 ! 3302 !-- Sum-up profiles for vertical fluxes of the the species. Note, in case of WS5 scheme the 3303 !-- profiles of resolved-scale fluxes have been already summed-up, while resolved-scale fluxes 3304 !-- need to be calculated in case of PW scheme. 3305 !-- For summation employ a temporary array. 3306 !$OMP PARALLEL PRIVATE( i, j, k, tn, lpr, sums_tmp, flag ) 3307 !$ tn = omp_get_thread_num() 3308 !$OMP DO 3309 DO lpr = 1, cs_pr_count_fl_sgs 3310 sums_tmp(:,tn) = 0.0_wp 3311 DO i = nxl, nxr 3312 DO j = nys, nyn 3313 DO k = nzb, nzt 3314 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 23 ) ) * & 3315 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 9 ) ) 3316 sums_tmp(k,tn) = sums_tmp(k,tn) - & 3317 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) & 3318 * ( chem_species(cs_pr_index_fl_sgs(lpr))%conc(k+1,j,i) - & 3319 chem_species(cs_pr_index_fl_sgs(lpr))%conc(k,j,i) ) * & 3320 ddzu(k+1) * rmask(j,i,sr) * flag 3321 ENDDO 3322 ! 3323 !-- Add surface fluxes 3324 surf_s = surf_def_h(0)%start_index(j,i) 3325 surf_e = surf_def_h(0)%end_index(j,i) 3326 DO m = surf_s, surf_e 3327 k = surf_def_h(0)%k(m) + surf_def_h(0)%koff 3328 sums_tmp(k,tn) = sums_tmp(k,tn) + surf_def_h(0)%cssws(cs_pr_index_fl_sgs(lpr),m) 3329 ENDDO 3330 surf_s = surf_def_h(1)%start_index(j,i) 3331 surf_e = surf_def_h(1)%end_index(j,i) 3332 DO m = surf_s, surf_e 3333 k = surf_def_h(1)%k(m) + surf_def_h(1)%koff 3334 sums_tmp(k,tn) = sums_tmp(k,tn) + surf_def_h(1)%cssws(cs_pr_index_fl_sgs(lpr),m) 3335 ENDDO 3336 surf_s = surf_lsm_h%start_index(j,i) 3337 surf_e = surf_lsm_h%end_index(j,i) 3338 DO m = surf_s, surf_e 3339 k = surf_lsm_h%k(m) + surf_lsm_h%koff 3340 sums_tmp(k,tn) = sums_tmp(k,tn) + surf_lsm_h%cssws(cs_pr_index_fl_sgs(lpr),m) 3341 ENDDO 3342 surf_s = surf_usm_h%start_index(j,i) 3343 surf_e = surf_usm_h%end_index(j,i) 3344 DO m = surf_s, surf_e 3345 k = surf_usm_h%k(m) + surf_usm_h%koff 3346 sums_tmp(k,tn) = sums_tmp(k,tn) + surf_usm_h%cssws(cs_pr_index_fl_sgs(lpr),m) 3347 ENDDO 3348 ENDDO 3349 ENDDO 3350 sums_l(nzb:nzt+1,hom_index_fl_sgs(lpr),tn) = sums_tmp(nzb:nzt+1,tn) 3351 ENDDO 3352 ! 3353 !-- Resolved-scale fluxes from the WS5 scheme 3354 IF ( ws_scheme_sca ) THEN 3355 !$OMP PARALLEL PRIVATE( tn, lpr ) 3356 !$ tn = omp_get_thread_num() 3357 !$OMP DO 3358 DO lpr = 1, cs_pr_count_fl_res 3359 sums_l(nzb:nzt+1,hom_index_fl_res(lpr),tn) = & 3360 sums_ws_l(nzb:nzt+1,tn,cs_pr_index_fl_res(lpr)) 3361 ENDDO 3362 ENDIF 3363 3204 3364 ELSEIF ( mode == 'time_series' ) THEN 3205 3365 ! @todo
Note: See TracChangeset
for help on using the changeset viewer.
