Ignore:
Timestamp:
Jun 29, 2020 8:49:58 AM (12 months ago)
Author:
suehring
Message:

mesoscale nesting: omit explicit pressure forcing via geostrophic wind components; chemistry: enable profile output of vertical fluxes; urban-surface: bugfix in initialization in case of cyclic_fill

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/chemistry_model_mod.f90

    r4577 r4581  
    2626! -----------------
    2727! $Id$
     28! Enable output of resolved-scale vertical fluxes of chemical species.
     29!
     30! 4577 2020-06-25 09:53:58Z raasch
    2831! further re-formatting to follow the PALM coding standard
    2932!
     
    436439!
    437440!-- Interface section
    438     INTERFACE chem_3d_data_averaging
    439        MODULE PROCEDURE chem_3d_data_averaging
    440     END INTERFACE chem_3d_data_averaging
    441 
    442     INTERFACE chem_boundary_conditions
    443        MODULE PROCEDURE chem_boundary_conditions
    444     END INTERFACE chem_boundary_conditions
    445 
    446     INTERFACE chem_check_data_output
    447        MODULE PROCEDURE chem_check_data_output
    448     END INTERFACE chem_check_data_output
    449 
    450     INTERFACE chem_data_output_2d
    451        MODULE PROCEDURE chem_data_output_2d
    452     END INTERFACE chem_data_output_2d
    453 
    454     INTERFACE chem_data_output_3d
    455        MODULE PROCEDURE chem_data_output_3d
    456     END INTERFACE chem_data_output_3d
    457 
    458     INTERFACE chem_data_output_mask
    459        MODULE PROCEDURE chem_data_output_mask
    460     END INTERFACE chem_data_output_mask
    461 
    462     INTERFACE chem_check_data_output_pr
    463        MODULE PROCEDURE chem_check_data_output_pr
    464     END INTERFACE chem_check_data_output_pr
    465 
    466     INTERFACE chem_check_parameters
    467        MODULE PROCEDURE chem_check_parameters
    468     END INTERFACE chem_check_parameters
    469 
    470     INTERFACE chem_define_netcdf_grid
    471        MODULE PROCEDURE chem_define_netcdf_grid
    472     END INTERFACE chem_define_netcdf_grid
    473 
    474     INTERFACE chem_header
    475        MODULE PROCEDURE chem_header
    476     END INTERFACE chem_header
    477 
    478     INTERFACE chem_init_arrays
    479        MODULE PROCEDURE chem_init_arrays
    480     END INTERFACE chem_init_arrays
    481 
    482     INTERFACE chem_init
    483        MODULE PROCEDURE chem_init
    484     END INTERFACE chem_init
    485 
    486     INTERFACE chem_init_profiles
    487        MODULE PROCEDURE chem_init_profiles
    488     END INTERFACE chem_init_profiles
    489 
    490     INTERFACE chem_integrate
    491        MODULE PROCEDURE chem_integrate_ij
    492     END INTERFACE chem_integrate
    493 
    494     INTERFACE chem_parin
    495        MODULE PROCEDURE chem_parin
    496     END INTERFACE chem_parin
    497 
    498441    INTERFACE chem_actions
    499442       MODULE PROCEDURE chem_actions
    500443       MODULE PROCEDURE chem_actions_ij
    501444    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
    502505
    503506    INTERFACE chem_non_advective_processes
     
    10691072
    10701073
    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:) )
    10871107          DO  lsp = 1, nspec
    10881108             IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
    1089                 cs_pr_count = cs_pr_count+1
    1090                 cs_pr_index(cs_pr_count) = lsp
    1091                 dopr_index(var_count) = pr_palm+cs_pr_count
    1092                 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'
    10931113                IF ( spec_name(1:2) == 'PM')  THEN
    1094                      dopr_unit = 'kg m-3'
     1114                   dopr_unit = 'kg m-3'
    10951115                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)
    10981120             ENDIF
    10991121          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
    11011180
    11021181 END SUBROUTINE chem_check_data_output_pr
     
    19762055    spec_conc_3 (:,:,:,:) = 0.0_wp
    19772056
     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
    19782063
    19792064    DO  lsp = 1, nspec
     
    20702155             ENDDO
    20712156          ENDIF
     2157
    20722158       ENDIF
    20732159!
     
    26202706!
    26212707!-- Determine the number of user-defined profiles and append them to the standard data output
    2622 !>> (data_output_pr)
     2708!-- (data_output_pr)
    26232709    max_pr_cs_tmp = 0
    26242710    i = 1
    26252711
    26262712    DO  WHILE ( data_output_pr(i)  /= ' '  .AND.  i <= 100 )
     2713
    26272714       IF ( TRIM( data_output_pr(i)(1:3) ) == 'kc_' )  THEN
    26282715          max_pr_cs_tmp = max_pr_cs_tmp+1
     
    27082795          ENDIF
    27092796
     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
    27102804       CASE DEFAULT
    27112805          CONTINUE
     
    28762970       IF ( timestep_scheme(1:5) == 'runge' )  THEN
    28772971          IF ( ws_scheme_sca )  THEN
     2972             sums_wschs_ws_l(nzb:,0:) => sums_ws_l(:,:,ilsp)
    28782973             CALL advec_s_ws( cs_advc_flags_s, chem_species(ilsp)%conc, 'kc',                      &
    28792974                              bc_dirichlet_cs_l  .OR.  bc_radiation_cs_l,                          &
     
    29863081       IF ( timestep_scheme(1:5) == 'runge' )  THEN
    29873082          IF ( ws_scheme_sca )  THEN
     3083             sums_wschs_ws_l(nzb:,0:) => sums_ws_l(nzb:nzt+1,0:threads_per_task-1,ilsp)
    29883084             CALL advec_s_ws( cs_advc_flags_s,                                                     &
    29893085                              i,                                                                   &
     
    30293125               * ( chem_species(ilsp)%conc(k,j,i) - chem_species(ilsp)%conc_pr_init(k) )           &
    30303126               )                                                                                   &
    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 )                      &
    30333128               )
    30343129
     
    31703265    CHARACTER (LEN=*) ::  mode   !<
    31713266
    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
    31783279
    31793280    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) *                            &
    31973294                        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 ) )
    32003296                ENDDO
    32013297             ENDDO
    32023298          ENDDO
     3299          sums_l(nzb:nzt+1,hom_index_spec(lpr),tn) = sums_tmp(nzb:nzt+1,tn)
    32033300       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
    32043364    ELSEIF ( mode == 'time_series' )  THEN
    32053365!      @todo
Note: See TracChangeset for help on using the changeset viewer.