- Timestamp:
- Apr 30, 2020 12:20:40 PM (5 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 4 edited
-
check_parameters.f90 (modified) (3 diffs)
-
chem_modules.f90 (modified) (4 diffs)
-
chemistry_model_mod.f90 (modified) (33 diffs)
-
time_integration.f90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/check_parameters.f90
r4495 r4511 25 25 ! ----------------- 26 26 ! $Id$ 27 ! call of chem_boundary_conds removed (respective settings are now done in the chemistry module) 28 ! 29 ! 4495 2020-04-13 20:11:20Z raasch 27 30 ! check new restart_data_format parameters 28 31 ! … … 131 134 132 135 USE chem_modules 133 134 USE chemistry_model_mod, &135 ONLY: chem_boundary_conds136 136 137 137 USE control_parameters … … 1554 1554 1555 1555 ! 1556 !-- Boundary conditions for chemical species1557 IF ( air_chemistry ) CALL chem_boundary_conds( 'init' )1558 1559 !1560 1556 !-- Boundary conditions for horizontal components of wind speed 1561 1557 IF ( bc_uv_b == 'dirichlet' ) THEN -
palm/trunk/SOURCE/chem_modules.f90
r4481 r4511 22 22 ! Current revisions: 23 23 ! ----------------- 24 ! 24 ! 25 25 ! 26 26 ! Former revisions: 27 27 ! ----------------- 28 28 ! $Id$ 29 ! new variables for explicit settings of lateral boundary conditions introduced 30 ! 31 ! 4481 2020-03-31 18:55:54Z maronga 29 32 ! added namelist flag 'emiss_read_legacy_mode' to allow concurrent 30 33 ! functioning of new emission read mode under development (ECC) … … 92 95 93 96 CHARACTER (LEN=20) :: bc_cs_b = 'dirichlet' !< namelist parameter: surface boundary condition for concentration 97 CHARACTER (LEN=20) :: bc_cs_l = 'undefined' !< left boundary condition 98 CHARACTER (LEN=20) :: bc_cs_n = 'undefined' !< north boundary condition 99 CHARACTER (LEN=20) :: bc_cs_r = 'undefined' !< right boundary condition 100 CHARACTER (LEN=20) :: bc_cs_s = 'undefined' !< south boundary condition 94 101 CHARACTER (LEN=20) :: bc_cs_t = 'initial_gradient' !< namelist parameter: top boudary condition for concentration 95 102 CHARACTER (LEN=30) :: chem_mechanism = 'phstatp' !< namelist parameter: chemistry mechanism … … 105 112 CHARACTER (LEN=11), DIMENSION(99) :: data_output_pr_cs = 'novalue' !< namelist parameter: tbc...??? 106 113 CHARACTER (LEN=11), DIMENSION(99) :: surface_csflux_name = 'novalue' !< namelist parameter: tbc...??? 114 115 INTEGER(iwp) :: communicator_chem !< stores the number of the MPI communicator to be used 116 !< for ghost layer data exchange 117 !< 1: cyclic, 2: cyclic along x, 3: cyclic along y, 118 !< 4: non-cyclic 107 119 108 120 INTEGER(iwp) :: cs_pr_count = 0 !< counter for chemical species profiles … … 131 143 !< chemical species near walls and lateral boundaries 132 144 145 LOGICAL :: bc_dirichlet_cs_l = .FALSE. !< flag for indicating a dirichlet condition at 146 !< the left boundary 147 LOGICAL :: bc_dirichlet_cs_n = .FALSE. !< flag for indicating a dirichlet condition at 148 !< the north boundary 149 LOGICAL :: bc_dirichlet_cs_r = .FALSE. !< flag for indicating a dirichlet condition at 150 !< the right boundary 151 LOGICAL :: bc_dirichlet_cs_s = .FALSE. !< flag for indicating a dirichlet condition at 152 !< the south boundary 153 LOGICAL :: bc_radiation_cs_l = .FALSE. !< flag for indicating a radiation/neumann 154 !< condition at the left boundary 155 LOGICAL :: bc_radiation_cs_n = .FALSE. !< flag for indicating a radiation/neumann 156 !< condition at the north boundary 157 LOGICAL :: bc_radiation_cs_r = .FALSE. !< flag for indicating a radiation/neumann 158 !< condition at the right boundary 159 LOGICAL :: bc_radiation_cs_s = .FALSE. !< flag for indicating a radiation/neumann 160 !< condition at the south boundary 133 161 LOGICAL :: constant_top_csflux(99) = .TRUE. !< internal flag, set to .FALSE. if no top_csflux is prescribed 134 162 LOGICAL :: constant_csflux(99) = .TRUE. !< internal flag, set to .FALSE. if no surface_csflux is prescribed -
palm/trunk/SOURCE/chemistry_model_mod.f90
r4487 r4511 27 27 ! ----------------- 28 28 ! $Id$ 29 ! decycling replaced by explicit setting of lateral boundary conditions 30 ! 31 ! 4487 2020-04-03 09:38:20Z raasch 29 32 ! bugfix for subroutine calls that contain the decycle_chem switches as arguments 30 ! 33 ! 31 34 ! 4481 2020-03-31 18:55:54Z maronga 32 35 ! use statement for exchange horiz added, … … 347 350 REAL(kind=wp), DIMENSION(nkppctrl) :: rcntrl !< 20 real parameters for fine tuning of KPP code 348 351 !< (e.g starting internal timestep of solver) 349 !350 !-- Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not351 !-- approproate for chemicals compounds since they may accumulate too much.352 !-- If no proper boundary conditions from a DYNAMIC input file are available,353 !-- de-cycling applies the initial profiles at the inflow boundaries for354 !-- Dirichlet boundary conditions355 LOGICAL :: decycle_chem_lr = .FALSE. !< switch for setting decycling in left-right direction356 LOGICAL :: decycle_chem_ns = .FALSE. !< switch for setting decycling in south-norht direction357 CHARACTER (LEN=20), DIMENSION(4) :: decycle_method = &358 (/'dirichlet','dirichlet','dirichlet','dirichlet'/)359 !< Decycling method at horizontal boundaries,360 !< 1=left, 2=right, 3=south, 4=north361 !< dirichlet = initial size distribution and362 !< chemical composition set for the ghost and363 !< first three layers364 !< neumann = zero gradient365 352 366 353 REAL(kind=wp), PUBLIC :: cs_time_step = 0._wp … … 438 425 END INTERFACE chem_3d_data_averaging 439 426 440 INTERFACE chem_boundary_conds441 MODULE PROCEDURE chem_boundary_conds442 MODULE PROCEDURE chem_boundary_conds_decycle443 END INTERFACE chem_boundary_conds444 445 427 INTERFACE chem_boundary_conditions 446 428 MODULE PROCEDURE chem_boundary_conditions … … 590 572 END INTERFACE rc_rctot 591 573 592 ! INTERFACE rc_comp_point_rc_eff593 ! MODULE PROCEDURE rc_comp_point_rc_eff594 ! END INTERFACE rc_comp_point_rc_eff595 596 574 INTERFACE drydepo_aero_zhang_vd 597 575 MODULE PROCEDURE drydepo_aero_zhang_vd … … 604 582 605 583 606 PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_boundary_conditions, & 607 chem_boundary_conds_decycle, chem_check_data_output, & 608 chem_check_data_output_pr, chem_check_parameters, & 609 chem_data_output_2d, chem_data_output_3d, chem_data_output_mask, & 610 chem_define_netcdf_grid, chem_header, chem_init, chem_init_arrays, & 611 chem_init_profiles, chem_integrate, chem_parin, & 612 chem_actions, chem_prognostic_equations, chem_rrd_local, & 613 chem_statistics, chem_swap_timelevel, chem_wrd_local, chem_depo, & 614 chem_non_advective_processes, chem_exchange_horiz_bounds 584 PUBLIC chem_3d_data_averaging, chem_boundary_conditions, chem_check_data_output, & 585 chem_check_data_output_pr, chem_check_parameters, chem_data_output_2d, & 586 chem_data_output_3d, chem_data_output_mask, chem_define_netcdf_grid, chem_header, & 587 chem_init, chem_init_arrays, chem_init_profiles, chem_integrate, chem_parin, & 588 chem_actions, chem_prognostic_equations, chem_rrd_local, chem_statistics, & 589 chem_swap_timelevel, chem_wrd_local, chem_depo, chem_non_advective_processes, & 590 chem_exchange_horiz_bounds 615 591 616 592 CONTAINS … … 742 718 !> Subroutine to initialize and set all boundary conditions for chemical species 743 719 !------------------------------------------------------------------------------! 744 SUBROUTINE chem_boundary_conds( mode ) 745 746 USE control_parameters, & 747 ONLY: bc_radiation_l, bc_radiation_n, bc_radiation_r, bc_radiation_s 748 749 USE arrays_3d, & 750 ONLY: dzu 720 SUBROUTINE chem_boundary_conditions( horizontal_conditions_only ) 721 722 USE arrays_3d, & 723 ONLY: dzu 751 724 752 725 USE surface_mod, & 753 726 ONLY: bc_h 754 727 755 CHARACTER (LEN=*), INTENT(IN) :: mode756 728 INTEGER(iwp) :: i !< grid index x direction. 757 729 INTEGER(iwp) :: j !< grid index y direction. … … 761 733 INTEGER(iwp) :: lsp !< running index for chem spcs. 762 734 763 764 SELECT CASE ( TRIM( mode ) ) 765 CASE ( 'init' ) 766 767 IF ( bc_cs_b == 'dirichlet' ) THEN 768 ibc_cs_b = 0 769 ELSEIF ( bc_cs_b == 'neumann' ) THEN 770 ibc_cs_b = 1 771 ELSE 772 message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"' 773 CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 ) 774 ENDIF 775 ! 776 !-- Set Integer flags and check for possible erroneous settings for top 777 !-- boundary condition. 778 IF ( bc_cs_t == 'dirichlet' ) THEN 779 ibc_cs_t = 0 780 ELSEIF ( bc_cs_t == 'neumann' ) THEN 781 ibc_cs_t = 1 782 ELSEIF ( bc_cs_t == 'initial_gradient' ) THEN 783 ibc_cs_t = 2 784 ELSEIF ( bc_cs_t == 'nested' ) THEN 785 ibc_cs_t = 3 786 ELSE 787 message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"' 788 CALL message( 'check_parameters', 'CM0430', 1, 2, 0, 6, 0 ) 789 ENDIF 790 791 CASE ( 'set_bc_bottomtop' ) 792 ! 793 !-- Boundary condtions for chemical species at horizontal walls 794 DO lsp = 1, nspec 795 IF ( ibc_cs_b == 0 ) THEN 796 DO l = 0, 1 797 !$OMP PARALLEL DO PRIVATE( i, j, k ) 798 DO m = 1, bc_h(l)%ns 799 i = bc_h(l)%i(m) 800 j = bc_h(l)%j(m) 801 k = bc_h(l)%k(m) 802 chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) = & 735 LOGICAL, OPTIONAL :: horizontal_conditions_only !< switch to set horizontal bc only 736 737 738 IF ( .NOT. PRESENT( horizontal_conditions_only ) ) THEN 739 ! 740 !-- Boundary condtions for chemical species at horizontal walls 741 DO lsp = 1, nspec 742 743 IF ( ibc_cs_b == 0 ) THEN 744 DO l = 0, 1 745 !$OMP PARALLEL DO PRIVATE( i, j, k ) 746 DO m = 1, bc_h(l)%ns 747 i = bc_h(l)%i(m) 748 j = bc_h(l)%j(m) 749 k = bc_h(l)%k(m) 750 chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) = & 803 751 chem_species(lsp)%conc(k+bc_h(l)%koff,j,i) 804 ENDDO 805 ENDDO 806 807 ELSEIF ( ibc_cs_b == 1 ) THEN 808 ! 809 !-- in boundary_conds there is som extra loop over m here for passive tracer 810 DO l = 0, 1 811 !$OMP PARALLEL DO PRIVATE( i, j, k ) 812 DO m = 1, bc_h(l)%ns 813 i = bc_h(l)%i(m) 814 j = bc_h(l)%j(m) 815 k = bc_h(l)%k(m) 816 chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) = & 752 ENDDO 753 ENDDO 754 755 ELSEIF ( ibc_cs_b == 1 ) THEN 756 ! 757 !-- In boundary_conds there is som extra loop over m here for passive tracer 758 !> TODO: clarify the meaning of the above comment. Explain in more detail or remove it. (Siggi) 759 DO l = 0, 1 760 !$OMP PARALLEL DO PRIVATE( i, j, k ) 761 DO m = 1, bc_h(l)%ns 762 i = bc_h(l)%i(m) 763 j = bc_h(l)%j(m) 764 k = bc_h(l)%k(m) 765 chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) = & 817 766 chem_species(lsp)%conc_p(k,j,i) 818 819 ENDDO820 767 ENDDO 821 ENDIF822 ENDDO ! end lsp loop823 !824 !-- Top boundary conditions for chemical species - Should this not be done for all species?825 IF ( ibc_cs_t == 0 ) THEN826 DO lsp = 1, nspec827 chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)828 ENDDO829 ELSEIF ( ibc_cs_t == 1 ) THEN830 DO lsp = 1, nspec831 chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)832 ENDDO833 ELSEIF ( ibc_cs_t == 2 ) THEN834 DO lsp = 1, nspec835 chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1)836 768 ENDDO 837 769 ENDIF 838 770 839 CASE ( 'set_bc_lateral' ) 840 ! 841 !-- Lateral boundary conditions for chem species at inflow boundary 842 !-- are automatically set when chem_species concentration is 843 !-- initialized. The initially set value at the inflow boundary is not 844 !-- touched during time integration, hence, this boundary value remains 845 !-- at a constant value, which is the concentration that flows into the 846 !-- domain. 847 !-- Lateral boundary conditions for chem species at outflow boundary 848 849 IF ( bc_radiation_s ) THEN 850 DO lsp = 1, nspec 851 chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:) 852 ENDDO 853 ELSEIF ( bc_radiation_n ) THEN 854 DO lsp = 1, nspec 855 chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:) 856 ENDDO 857 ELSEIF ( bc_radiation_l ) THEN 858 DO lsp = 1, nspec 859 chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl) 860 ENDDO 861 ELSEIF ( bc_radiation_r ) THEN 862 DO lsp = 1, nspec 863 chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr) 864 ENDDO 865 ENDIF 866 867 END SELECT 868 869 END SUBROUTINE chem_boundary_conds 870 871 !------------------------------------------------------------------------------! 872 ! Description: 873 ! ------------ 874 !> Subroutine for boundary conditions 875 !------------------------------------------------------------------------------! 876 SUBROUTINE chem_boundary_conditions 877 878 IMPLICIT NONE 879 880 INTEGER(iwp) :: lsp !< 881 INTEGER(iwp) :: lsp_usr !< 882 883 ! 884 !-- Top/bottom boundary conditions for chemical species 885 CALL chem_boundary_conds( 'set_bc_bottomtop' ) 886 ! 887 !-- Lateral boundary conditions for chemical species 888 CALL chem_boundary_conds( 'set_bc_lateral' ) 889 890 ! 891 !-- Boundary conditions for prognostic quantitites of other modules: 892 !-- Here, only decycling is carried out 893 DO lsp = 1, nvar 894 lsp_usr = 1 895 DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' ) 896 IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) ) THEN 897 CALL chem_boundary_conds( chem_species(lsp)%conc_p, & 898 chem_species(lsp)%conc_pr_init ) 899 ENDIF 900 lsp_usr = lsp_usr + 1 901 ENDDO 902 ENDDO 903 904 905 END SUBROUTINE chem_boundary_conditions 906 907 !------------------------------------------------------------------------------! 908 ! Description: 909 ! ------------ 910 !> Boundary conditions for prognostic variables in chemistry: decycling in the 911 !> x-direction- 912 !> Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not 913 !> approproate for chemicals compounds since they may accumulate too much. 914 !> If no proper boundary conditions from a DYNAMIC input file are available, 915 !> de-cycling applies the initial profiles at the inflow boundaries for 916 !> Dirichlet boundary conditions 917 !------------------------------------------------------------------------------! 918 SUBROUTINE chem_boundary_conds_decycle( cs_3d, cs_pr_init ) 919 920 USE control_parameters, & 921 ONLY: nesting_offline 922 923 INTEGER(iwp) :: boundary !< 924 INTEGER(iwp) :: ee !< 925 INTEGER(iwp) :: copied !< 926 INTEGER(iwp) :: i !< 927 INTEGER(iwp) :: j !< 928 INTEGER(iwp) :: k !< 929 INTEGER(iwp) :: ss !< 930 931 REAL(wp), DIMENSION(nzb:nzt+1) :: cs_pr_init 932 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: cs_3d 933 REAL(wp) :: flag !< flag to mask topography grid points 934 935 936 flag = 0.0_wp 937 ! 938 !-- Skip input if forcing from a larger-scale model is applied 939 IF ( nesting_offline .AND. nesting_offline_chem ) RETURN 940 ! 941 !-- Left and right boundaries 942 IF ( decycle_chem_lr .AND. bc_lr_cyc ) THEN 943 944 DO boundary = 1, 2 945 946 IF ( decycle_method(boundary) == 'dirichlet' ) THEN 947 ! 948 !-- Initial profile is copied to ghost and first three layers 949 ss = 1 950 ee = 0 951 IF ( boundary == 1 .AND. nxl == 0 ) THEN 952 ss = nxlg 953 ee = nxl-1 954 ELSEIF ( boundary == 2 .AND. nxr == nx ) THEN 955 ss = nxr+1 956 ee = nxrg 957 ENDIF 958 959 DO i = ss, ee 960 DO j = nysg, nyng 961 DO k = nzb+1, nzt 962 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 963 cs_3d(k,j,i) = cs_pr_init(k) * flag 964 ENDDO 965 ENDDO 966 ENDDO 967 968 ELSEIF ( decycle_method(boundary) == 'neumann' ) THEN 969 ! 970 !-- The value at the boundary is copied to the ghost layers to simulate 971 !-- an outlet with zero gradient 972 ss = 1 973 ee = 0 974 IF ( boundary == 1 .AND. nxl == 0 ) THEN 975 ss = nxlg 976 ee = nxl-1 977 copied = nxl 978 ELSEIF ( boundary == 2 .AND. nxr == nx ) THEN 979 ss = nxr+1 980 ee = nxrg 981 copied = nxr 982 ENDIF 983 984 DO i = ss, ee 985 DO j = nysg, nyng 986 DO k = nzb+1, nzt 987 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 988 cs_3d(k,j,i) = cs_3d(k,j,copied) * flag 989 ENDDO 990 ENDDO 991 ENDDO 992 993 ELSE 994 WRITE(message_string,*) & 995 'unknown decycling method: decycle_method (', & 996 boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"' 997 CALL message( 'chem_boundary_conds_decycle', 'CM0431', & 998 1, 2, 0, 6, 0 ) 771 ENDDO ! end lsp loop 772 773 ! 774 !-- Top boundary conditions for chemical species - Should this not be done for all species? 775 !> TODO: This question also needs to be clarified. I guess it can be removed because the loop 776 !> already runs over all species? (Siggi) 777 DO lsp = 1, nspec 778 IF ( ibc_cs_t == 0 ) THEN 779 chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:) 780 ELSEIF ( ibc_cs_t == 1 ) THEN 781 chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) 782 ELSEIF ( ibc_cs_t == 2 ) THEN 783 chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + & 784 bc_cs_t_val(lsp) * dzu(nzt+1) 999 785 ENDIF 1000 786 ENDDO 787 788 ! 789 !-- Lateral boundary conditions. 790 !-- Dirichlet conditions have been already set when chem_species concentration is initialized. 791 !-- The initially set value is not touched during time integration, hence, this boundary value 792 !-- remains at a constant value. 793 !-- Neumann conditions: 794 IF ( bc_radiation_cs_s ) THEN 795 DO lsp = 1, nspec 796 chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:) 797 ENDDO 798 ENDIF 799 IF ( bc_radiation_cs_n ) THEN 800 DO lsp = 1, nspec 801 chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:) 802 ENDDO 803 ENDIF 804 IF ( bc_radiation_cs_l ) THEN 805 DO lsp = 1, nspec 806 chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl) 807 ENDDO 808 ENDIF 809 IF ( bc_radiation_cs_r ) THEN 810 DO lsp = 1, nspec 811 chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr) 812 ENDDO 813 ENDIF 814 815 !-- For testing: set Dirichlet conditions for all boundaries 816 ! IF ( bc_dirichlet_cs_s ) THEN 817 ! DO lsp = 1, nspec 818 ! DO i = nxlg, nxrg 819 ! DO j = nysg, nys-1 820 ! DO k = nzb, nzt 821 ! IF ( k /= nzb ) THEN 822 ! flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 823 ! ELSE 824 ! flag = 1.0 825 ! ENDIF 826 ! chem_species(lsp)%conc_p(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag 827 ! ENDDO 828 ! ENDDO 829 ! ENDDO 830 ! ENDDO 831 ! ENDIF 832 ! IF ( bc_dirichlet_cs_n ) THEN 833 ! DO lsp = 1, nspec 834 ! DO i = nxlg, nxrg 835 ! DO j = nyn+1, nyng 836 ! DO k = nzb, nzt 837 ! IF ( k /= nzb ) THEN 838 ! flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 839 ! ELSE 840 ! flag = 1.0 841 ! ENDIF 842 ! chem_species(lsp)%conc_p(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag 843 ! ENDDO 844 ! ENDDO 845 ! ENDDO 846 ! ENDDO 847 ! ENDIF 848 ! IF ( bc_dirichlet_cs_l ) THEN 849 ! DO lsp = 1, nspec 850 ! DO i = nxlg, nxl-1 851 ! DO j = nysg, nyng 852 ! DO k = nzb, nzt 853 ! IF ( k /= nzb ) THEN 854 ! flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 855 ! ELSE 856 ! flag = 1.0 857 ! ENDIF 858 ! chem_species(lsp)%conc_p(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag 859 ! ENDDO 860 ! ENDDO 861 ! ENDDO 862 ! ENDDO 863 ! ENDIF 864 ! IF ( bc_dirichlet_cs_r ) THEN 865 ! DO lsp = 1, nspec 866 ! DO i = nxr+1, nxrg 867 ! DO j = nysg, nyng 868 ! DO k = nzb, nzt 869 ! IF ( k /= nzb ) THEN 870 ! flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 871 ! ELSE 872 ! flag = 1.0 873 ! ENDIF 874 ! chem_species(lsp)%conc_p(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag 875 ! ENDDO 876 ! ENDDO 877 ! ENDDO 878 ! ENDDO 879 ! ENDIF 880 881 ELSE 882 ! 883 !-- Lateral Neumann booundary conditions for timelevel t. 884 !-- This branch is executed when routine is called after the non-advective processes / 885 !-- before the prognostic equations. 886 IF ( bc_radiation_cs_s ) THEN 887 DO lsp = 1, nspec 888 chem_species(lsp)%conc(:,nys-1,:) = chem_species(lsp)%conc(:,nys,:) 889 ENDDO 890 ENDIF 891 IF ( bc_radiation_cs_n ) THEN 892 DO lsp = 1, nspec 893 chem_species(lsp)%conc(:,nyn+1,:) = chem_species(lsp)%conc(:,nyn,:) 894 ENDDO 895 ENDIF 896 IF ( bc_radiation_cs_l ) THEN 897 DO lsp = 1, nspec 898 chem_species(lsp)%conc(:,:,nxl-1) = chem_species(lsp)%conc(:,:,nxl) 899 ENDDO 900 ENDIF 901 IF ( bc_radiation_cs_r ) THEN 902 DO lsp = 1, nspec 903 chem_species(lsp)%conc(:,:,nxr+1) = chem_species(lsp)%conc(:,:,nxr) 904 ENDDO 905 ENDIF 906 907 !-- For testing: set Dirichlet conditions for all boundaries 908 ! IF ( bc_dirichlet_cs_s ) THEN 909 ! DO lsp = 1, nspec 910 ! DO i = nxlg, nxrg 911 ! DO j = nysg, nys-1 912 ! DO k = nzb, nzt 913 ! IF ( k /= nzb ) THEN 914 ! flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 915 ! ELSE 916 ! flag = 1.0 917 ! ENDIF 918 ! chem_species(lsp)%conc(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag 919 ! ENDDO 920 ! ENDDO 921 ! ENDDO 922 ! ENDDO 923 ! ENDIF 924 ! IF ( bc_dirichlet_cs_n ) THEN 925 ! DO lsp = 1, nspec 926 ! DO i = nxlg, nxrg 927 ! DO j = nyn+1, nyng 928 ! DO k = nzb, nzt 929 ! IF ( k /= nzb ) THEN 930 ! flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 931 ! ELSE 932 ! flag = 1.0 933 ! ENDIF 934 ! chem_species(lsp)%conc(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag 935 ! ENDDO 936 ! ENDDO 937 ! ENDDO 938 ! ENDDO 939 ! ENDIF 940 ! IF ( bc_dirichlet_cs_l ) THEN 941 ! DO lsp = 1, nspec 942 ! DO i = nxlg, nxl-1 943 ! DO j = nysg, nyng 944 ! DO k = nzb, nzt 945 ! IF ( k /= nzb ) THEN 946 ! flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 947 ! ELSE 948 ! flag = 1.0 949 ! ENDIF 950 ! chem_species(lsp)%conc(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag 951 ! ENDDO 952 ! ENDDO 953 ! ENDDO 954 ! ENDDO 955 ! ENDIF 956 ! IF ( bc_dirichlet_cs_r ) THEN 957 ! DO lsp = 1, nspec 958 ! DO i = nxr+1, nxrg 959 ! DO j = nysg, nyng 960 ! DO k = nzb, nzt 961 ! IF ( k /= nzb ) THEN 962 ! flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 963 ! ELSE 964 ! flag = 1.0 965 ! ENDIF 966 ! chem_species(lsp)%conc(k,j,i) = chem_species(lsp)%conc_pr_init(k) * flag 967 ! ENDDO 968 ! ENDDO 969 ! ENDDO 970 ! ENDDO 971 ! ENDIF 972 1001 973 ENDIF 1002 ! 1003 !-- South and north boundaries 1004 IF ( decycle_chem_ns .AND. bc_ns_cyc ) THEN 1005 1006 DO boundary = 3, 4 1007 1008 IF ( decycle_method(boundary) == 'dirichlet' ) THEN 1009 ! 1010 !-- Initial profile is copied to ghost and first three layers 1011 ss = 1 1012 ee = 0 1013 IF ( boundary == 3 .AND. nys == 0 ) THEN 1014 ss = nysg 1015 ee = nys-1 1016 ELSEIF ( boundary == 4 .AND. nyn == ny ) THEN 1017 ss = nyn+1 1018 ee = nyng 1019 ENDIF 1020 1021 DO i = nxlg, nxrg 1022 DO j = ss, ee 1023 DO k = nzb+1, nzt 1024 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 1025 cs_3d(k,j,i) = cs_pr_init(k) * flag 1026 ENDDO 1027 ENDDO 1028 ENDDO 1029 1030 1031 ELSEIF ( decycle_method(boundary) == 'neumann' ) THEN 1032 ! 1033 !-- The value at the boundary is copied to the ghost layers to simulate 1034 !-- an outlet with zero gradient 1035 ss = 1 1036 ee = 0 1037 IF ( boundary == 3 .AND. nys == 0 ) THEN 1038 ss = nysg 1039 ee = nys-1 1040 copied = nys 1041 ELSEIF ( boundary == 4 .AND. nyn == ny ) THEN 1042 ss = nyn+1 1043 ee = nyng 1044 copied = nyn 1045 ENDIF 1046 1047 DO i = nxlg, nxrg 1048 DO j = ss, ee 1049 DO k = nzb+1, nzt 1050 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 1051 cs_3d(k,j,i) = cs_3d(k,copied,i) * flag 1052 ENDDO 1053 ENDDO 1054 ENDDO 1055 1056 ELSE 1057 WRITE(message_string,*) & 1058 'unknown decycling method: decycle_method (', & 1059 boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"' 1060 CALL message( 'chem_boundary_conds_decycle', 'CM0432', & 1061 1, 2, 0, 6, 0 ) 1062 ENDIF 1063 ENDDO 1064 ENDIF 1065 1066 1067 END SUBROUTINE chem_boundary_conds_decycle 974 975 END SUBROUTINE chem_boundary_conditions 1068 976 1069 977 … … 1194 1102 SUBROUTINE chem_check_parameters 1195 1103 1104 USE control_parameters, & 1105 ONLY: bc_lr, bc_ns 1196 1106 1197 1107 LOGICAL :: found … … 1223 1133 ENDIF 1224 1134 ! 1225 !-- check for decycling of chem species1226 IF ( (.NOT. any(decycle_method == 'neumann') ) .AND. (.NOT. any(decycle_method == 'dirichlet') ) ) THEN1227 message_string = 'Incorrect boundary conditions. Only neumann or ' &1228 // 'dirichlet &available for decycling chemical species '1229 CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 )1230 ENDIF1231 !1232 1135 !-- check for chemical mechanism used 1233 1136 CALL get_mechanism_name … … 1236 1139 CALL message( 'chem_check_parameters', 'CM0462', 1, 2, 0, 6, 0 ) 1237 1140 ENDIF 1141 1142 ! 1143 !-- Check bottom boundary condition and set internal steering parameter 1144 IF ( bc_cs_b == 'dirichlet' ) THEN 1145 ibc_cs_b = 0 1146 ELSEIF ( bc_cs_b == 'neumann' ) THEN 1147 ibc_cs_b = 1 1148 ELSE 1149 message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"' 1150 CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 ) 1151 ENDIF 1152 ! 1153 !-- Check top boundary condition and set internal steering parameter 1154 IF ( bc_cs_t == 'dirichlet' ) THEN 1155 ibc_cs_t = 0 1156 ELSEIF ( bc_cs_t == 'neumann' ) THEN 1157 ibc_cs_t = 1 1158 ELSEIF ( bc_cs_t == 'initial_gradient' ) THEN 1159 ibc_cs_t = 2 1160 ELSEIF ( bc_cs_t == 'nested' ) THEN 1161 ibc_cs_t = 3 1162 ELSE 1163 message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"' 1164 CALL message( 'check_parameters', 'CM0430', 1, 2, 0, 6, 0 ) 1165 ENDIF 1166 1238 1167 ! 1239 1168 !-- If nesting_chem = .F., set top boundary condition to its default value … … 1242 1171 bc_cs_t = 'initial_gradient' 1243 1172 ENDIF 1173 1174 ! 1175 !-- Check left and right boundary conditions. First set default value if not set by user. 1176 IF ( bc_cs_l == 'undefined' ) THEN 1177 IF ( bc_lr == 'cyclic' ) THEN 1178 bc_cs_l = 'cyclic' 1179 ELSEIF ( bc_lr == 'dirichlet/radiation' ) THEN 1180 bc_cs_l = 'dirichlet' 1181 ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN 1182 bc_cs_l = 'neumann' 1183 ENDIF 1184 ENDIF 1185 IF ( bc_cs_r == 'undefined' ) THEN 1186 IF ( bc_lr == 'cyclic' ) THEN 1187 bc_cs_r = 'cyclic' 1188 ELSEIF ( bc_lr == 'dirichlet/radiation' ) THEN 1189 bc_cs_r = 'neumann' 1190 ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN 1191 bc_cs_r = 'dirichlet' 1192 ENDIF 1193 ENDIF 1194 IF ( bc_cs_l /= 'dirichlet' .AND. bc_cs_l /= 'neumann' .AND. bc_cs_l /= 'cyclic' ) THEN 1195 message_string = 'unknown boundary condition: bc_cs_l = "' // TRIM( bc_cs_l ) // '"' 1196 CALL message( 'chem_check_parameters','PA0505', 1, 2, 0, 6, 0 ) 1197 ENDIF 1198 IF ( bc_cs_r /= 'dirichlet' .AND. bc_cs_r /= 'neumann' .AND. bc_cs_r /= 'cyclic' ) THEN 1199 message_string = 'unknown boundary condition: bc_cs_r = "' // TRIM( bc_cs_r ) // '"' 1200 CALL message( 'chem_check_parameters','PA0551', 1, 2, 0, 6, 0 ) 1201 ENDIF 1202 ! 1203 !-- Check north and south boundary conditions. First set default value if not set by user. 1204 IF ( bc_cs_n == 'undefined' ) THEN 1205 IF ( bc_ns == 'cyclic' ) THEN 1206 bc_cs_n = 'cyclic' 1207 ELSEIF ( bc_ns == 'dirichlet/radiation' ) THEN 1208 bc_cs_n = 'dirichlet' 1209 ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN 1210 bc_cs_n = 'neumann' 1211 ENDIF 1212 ENDIF 1213 IF ( bc_cs_s == 'undefined' ) THEN 1214 IF ( bc_ns == 'cyclic' ) THEN 1215 bc_cs_s = 'cyclic' 1216 ELSEIF ( bc_ns == 'dirichlet/radiation' ) THEN 1217 bc_cs_s = 'neumann' 1218 ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN 1219 bc_cs_s = 'dirichlet' 1220 ENDIF 1221 ENDIF 1222 IF ( bc_cs_n /= 'dirichlet' .AND. bc_cs_n /= 'neumann' .AND. bc_cs_n /= 'cyclic' ) THEN 1223 message_string = 'unknown boundary condition: bc_cs_n = "' // TRIM( bc_cs_n ) // '"' 1224 CALL message( 'chem_check_parameters','PA0714', 1, 2, 0, 6, 0 ) 1225 ENDIF 1226 IF ( bc_cs_s /= 'dirichlet' .AND. bc_cs_s /= 'neumann' .AND. bc_cs_s /= 'cyclic' ) THEN 1227 message_string = 'unknown boundary condition: bc_cs_s = "' // TRIM( bc_cs_s ) // '"' 1228 CALL message( 'chem_check_parameters','PA0715', 1, 2, 0, 6, 0 ) 1229 ENDIF 1230 ! 1231 !-- Cyclic conditions must be set identically at opposing boundaries 1232 IF ( ( bc_cs_l == 'cyclic' .AND. bc_cs_r /= 'cyclic' ) .OR. & 1233 ( bc_cs_r == 'cyclic' .AND. bc_cs_l /= 'cyclic' ) ) THEN 1234 message_string = 'boundary conditions bc_cs_l and bc_cs_r must both be cyclic or non-cyclic' 1235 CALL message( 'chem_check_parameters','PA0716', 1, 2, 0, 6, 0 ) 1236 ENDIF 1237 IF ( ( bc_cs_n == 'cyclic' .AND. bc_cs_s /= 'cyclic' ) .OR. & 1238 ( bc_cs_s == 'cyclic' .AND. bc_cs_n /= 'cyclic' ) ) THEN 1239 message_string = 'boundary conditions bc_cs_n and bc_cs_s must both be cyclic or non-cyclic' 1240 CALL message( 'chem_check_parameters','PA0717', 1, 2, 0, 6, 0 ) 1241 ENDIF 1242 ! 1243 !-- Set the switches that control application of horizontal boundary conditions at the boundaries 1244 !-- of the total domain 1245 IF ( bc_cs_n == 'dirichlet' .AND. nyn == ny ) bc_dirichlet_cs_n = .TRUE. 1246 IF ( bc_cs_n == 'neumann' .AND. nyn == ny ) bc_radiation_cs_n = .TRUE. 1247 IF ( bc_cs_s == 'dirichlet' .AND. nys == 0 ) bc_dirichlet_cs_s = .TRUE. 1248 IF ( bc_cs_s == 'neumann' .AND. nys == 0 ) bc_radiation_cs_s = .TRUE. 1249 IF ( bc_cs_l == 'dirichlet' .AND. nxl == 0 ) bc_dirichlet_cs_l = .TRUE. 1250 IF ( bc_cs_l == 'neumann' .AND. nxl == 0 ) bc_radiation_cs_l = .TRUE. 1251 IF ( bc_cs_r == 'dirichlet' .AND. nxr == nx ) bc_dirichlet_cs_r = .TRUE. 1252 IF ( bc_cs_r == 'neumann' .AND. nxr == nx ) bc_radiation_cs_r = .TRUE. 1253 1254 ! 1255 !-- Set the communicator to be used for ghost layer data exchange 1256 !-- 1: cyclic, 2: cyclic along x, 3: cyclic along y, 4: non-cyclic 1257 IF ( bc_cs_l == 'cyclic' ) THEN 1258 IF ( bc_cs_s == 'cyclic' ) THEN 1259 communicator_chem = 1 1260 ELSE 1261 communicator_chem = 2 1262 ENDIF 1263 ELSE 1264 IF ( bc_cs_s == 'cyclic' ) THEN 1265 communicator_chem = 3 1266 ELSE 1267 communicator_chem = 4 1268 ENDIF 1269 ENDIF 1270 1244 1271 ! 1245 1272 !-- chem_check_parameters is called before the array chem_species is allocated! 1246 1273 !-- temporary switch of this part of the check 1247 1274 ! RETURN !bK commented 1248 1275 !> TODO: this workaround definitely needs to be removed from here!!! 1249 1276 CALL chem_init_internal 1250 1277 ! … … 1732 1759 ! 1733 1760 !-- initializatoin of Surface and profile chemical species 1734 1735 1761 lsp = 1 1736 1762 docsinit_chr ='Chemical species for initial surface and profile emissions: ' … … 1750 1776 IF ( nesting_chem ) WRITE( io, 12 ) nesting_chem 1751 1777 IF ( nesting_offline_chem ) WRITE( io, 13 ) nesting_offline_chem 1778 1779 WRITE( io, 14 ) TRIM( bc_cs_b ), TRIM( bc_cs_t ), TRIM( bc_cs_s ), TRIM( bc_cs_n ), & 1780 TRIM( bc_cs_l ), TRIM( bc_cs_r ) 1781 1752 1782 ! 1753 1783 !-- number of variable and fix chemical species and number of reactions 1754 1784 cs_fixed = nspec - nvar 1755 1756 1785 WRITE ( io, * ) ' --> Chemical Mechanism : ', cs_mech 1757 1786 WRITE ( io, * ) ' --> Chemical species, variable: ', nvar … … 1774 1803 12 FORMAT (/' Nesting for chemistry variables: ', L1 ) 1775 1804 13 FORMAT (/' Offline nesting for chemistry variables: ', L1 ) 1776 ! 1777 ! 1805 14 FORMAT (/' Boundary conditions for chemical species:', / & 1806 ' bottom/top: ',A10,' / ',A10, / & 1807 ' north/south: ',A10,' / ',A10, / & 1808 ' left/right: ',A10,' / ',A10) 1809 1778 1810 END SUBROUTINE chem_header 1779 1811 … … 1884 1916 INTEGER(iwp) :: lpr_lev !< running index for chem spcs profile level 1885 1917 1918 REAL(wp) :: flag !< flag for masking topography/building grid points 1886 1919 ! 1887 1920 !-- 20200203 ECC … … 1948 1981 1949 1982 ENDDO 1950 ! 1951 !-- Set control flags for decycling only at lateral boundary cores, within the 1952 !-- inner cores the decycle flag is set to .False.. Even though it does not 1953 !-- affect the setting of chemistry boundary conditions, this flag is used to 1954 !-- set advection control flags appropriately. 1955 decycle_chem_lr = MERGE( decycle_chem_lr, .FALSE., & 1956 nxl == 0 .OR. nxr == nx ) 1957 decycle_chem_ns = MERGE( decycle_chem_ns, .FALSE., & 1958 nys == 0 .OR. nyn == ny ) 1983 1959 1984 ! 1960 1985 !-- For some passive scalars decycling may be enabled. This case, the lateral … … 1982 2007 !-- Initialize with flag 31 only. 1983 2008 cs_advc_flags_s = 0 1984 cs_advc_flags_s = MERGE( IBSET( cs_advc_flags_s, 31 ), 0, & 1985 BTEST( wall_flags_total_0, 31 ) ) 1986 1987 IF ( decycle_chem_ns ) THEN 2009 cs_advc_flags_s = MERGE( IBSET( cs_advc_flags_s, 31 ), 0, BTEST( wall_flags_total_0, 31 ) ) 2010 2011 IF ( bc_dirichlet_cs_n .OR. bc_dirichlet_cs_s ) THEN 1988 2012 IF ( nys == 0 ) THEN 1989 DO i = 1, nbgp 2013 DO i = 1, nbgp 1990 2014 cs_advc_flags_s(:,nys-i,:) = MERGE( & 1991 2015 IBSET( cs_advc_flags_s(:,nys,:), 31 ), & … … 1996 2020 ENDIF 1997 2021 IF ( nyn == ny ) THEN 1998 DO i = 1, nbgp 2022 DO i = 1, nbgp 1999 2023 cs_advc_flags_s(:,nyn+i,:) = MERGE( & 2000 2024 IBSET( cs_advc_flags_s(:,nyn,:), 31 ), & … … 2005 2029 ENDIF 2006 2030 ENDIF 2007 IF ( decycle_chem_lr ) THEN2031 IF ( bc_dirichlet_cs_l .OR. bc_dirichlet_cs_r ) THEN 2008 2032 IF ( nxl == 0 ) THEN 2009 DO i = 1, nbgp 2033 DO i = 1, nbgp 2010 2034 cs_advc_flags_s(:,:,nxl-i) = MERGE( & 2011 2035 IBSET( cs_advc_flags_s(:,:,nxl), 31 ), & … … 2015 2039 ENDDO 2016 2040 ENDIF 2017 IF ( nxr == nx ) THEN 2018 DO i = 1, nbgp 2041 IF ( nxr == nx ) THEN 2042 DO i = 1, nbgp 2019 2043 cs_advc_flags_s(:,:,nxr+i) = MERGE( & 2020 2044 IBSET( cs_advc_flags_s(:,:,nxr), 31 ), & … … 2023 2047 ) 2024 2048 ENDDO 2025 ENDIF 2049 ENDIF 2026 2050 ENDIF 2027 2051 ! … … 2039 2063 !-- oscillations, which are responsible for high concentration maxima that may 2040 2064 !-- appear under shear-free stable conditions. 2041 CALL ws_init_flags_scalar( & 2042 bc_dirichlet_l .OR. bc_radiation_l .OR. ( decycle_chem_lr .AND. nxl == 0 ), & 2043 bc_dirichlet_n .OR. bc_radiation_n .OR. ( decycle_chem_ns .AND. nyn == ny ), & 2044 bc_dirichlet_r .OR. bc_radiation_r .OR. ( decycle_chem_lr .AND. nxr == nx ), & 2045 bc_dirichlet_s .OR. bc_radiation_s .OR. ( decycle_chem_ns .AND. nys == 0 ), & 2046 cs_advc_flags_s, .TRUE. ) 2065 CALL ws_init_flags_scalar( bc_dirichlet_cs_l .OR. bc_radiation_cs_l, & 2066 bc_dirichlet_cs_n .OR. bc_radiation_cs_n, & 2067 bc_dirichlet_cs_r .OR. bc_radiation_cs_r, & 2068 bc_dirichlet_cs_s .OR. bc_radiation_cs_s, & 2069 cs_advc_flags_s, .TRUE. ) 2047 2070 ENDIF 2048 2071 ! … … 2074 2097 ! 2075 2098 !-- Transfer initial profiles to the arrays of the 3D model 2099 !-- Concentrations within buildings are set to zero. 2076 2100 DO lsp = 1, nspec 2077 2101 DO i = nxlg, nxrg 2078 2102 DO j = nysg, nyng 2079 2103 DO lpr_lev = 1, nz + 1 2080 chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev) 2104 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(lpr_lev,j,i), 0 ) ) 2105 chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)& 2106 * flag 2081 2107 ENDDO 2082 2108 ENDDO … … 2090 2116 DO i = nxlg, nxrg 2091 2117 DO j = nysg, nyng 2092 chem_species(lsp)%conc(:,j,i) = chem_species(lsp)%conc_pr_init 2118 DO lpr_lev = nzb, nz+1 2119 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(lpr_lev,j,i), 0 ) ) 2120 chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)& 2121 * flag 2122 ENDDO 2093 2123 ENDDO 2094 2124 ENDDO … … 2100 2130 IF ( cs_surface_initial_change(1) /= 0.0_wp ) THEN 2101 2131 DO lsp = 1, nspec 2102 chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) + & 2103 cs_surface_initial_change(lsp) 2132 DO i = nxlg, nxrg 2133 DO j = nysg, nyng 2134 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb,j,i), 0 ) ) 2135 chem_species(lsp)%conc(nzb,j,i) = chem_species(lsp)%conc(nzb,j,i) + & 2136 cs_surface_initial_change(lsp) * flag 2137 ENDDO 2138 ENDDO 2104 2139 ENDDO 2105 2140 ENDIF … … 2346 2381 2347 2382 NAMELIST /chemistry_parameters/ bc_cs_b, & 2383 bc_cs_l, & 2384 bc_cs_n, & 2385 bc_cs_r, & 2386 bc_cs_s, & 2348 2387 bc_cs_t, & 2349 2388 call_chem_at_all_substeps, & … … 2360 2399 cs_vertical_gradient_level, & 2361 2400 daytype_mdh, & 2362 decycle_chem_lr, &2363 decycle_chem_ns, &2364 decycle_method, &2365 2401 deposition_dry, & 2366 2402 emissions_anthropogenic, & … … 2590 2626 solver_type = 'Rang3' 2591 2627 ELSE 2592 message_string = 'illegal solver type'2628 message_string = 'illegal Rosenbrock-solver type' 2593 2629 CALL message( 'chem_parin', 'PA0506', 1, 2, 0, 6, 0 ) 2594 2630 END IF … … 2774 2810 2775 2811 INTEGER(iwp) :: lsp !< 2776 INTEGER(iwp) :: lsp_usr !<2777 2812 2778 2813 ! … … 2780 2815 CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'start' ) 2781 2816 DO lsp = 1, nvar 2782 CALL exchange_horiz( chem_species(lsp)%conc, nbgp ) 2783 lsp_usr = 1 2784 DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' ) 2785 IF ( TRIM(chem_species(lsp)%name) == TRIM(cs_name(lsp_usr)) ) THEN 2786 ! 2787 !-- As chem_exchange_horiz_bounds is called at the beginning 2788 !-- of prognostic_equations, boundary conditions are set on 2789 !-- %conc. 2790 CALL chem_boundary_conds( chem_species(lsp)%conc, & 2791 chem_species(lsp)%conc_pr_init ) 2792 2793 ENDIF 2794 lsp_usr = lsp_usr + 1 2795 ENDDO 2796 2797 2817 CALL exchange_horiz( chem_species(lsp)%conc, nbgp, & 2818 alternative_communicator = communicator_chem ) 2798 2819 ENDDO 2820 2821 CALL chem_boundary_conditions( horizontal_conditions_only = .TRUE. ) 2822 2799 2823 CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'stop' ) 2800 2824 … … 2832 2856 IF ( ws_scheme_sca ) THEN 2833 2857 CALL advec_s_ws( cs_advc_flags_s, chem_species(ilsp)%conc, 'kc', & 2834 bc_dirichlet_l .OR. bc_radiation_l .OR. ( decycle_chem_lr .AND. nxl == 0 ),&2835 bc_dirichlet_n .OR. bc_radiation_n .OR. ( decycle_chem_ns .AND. nyn == ny ),&2836 bc_dirichlet_r .OR. bc_radiation_r .OR. ( decycle_chem_lr .AND. nxr == nx ),&2837 bc_dirichlet_s .OR. bc_radiation_s .OR. ( decycle_chem_ns .AND. nys == 0 ))2858 bc_dirichlet_cs_l .OR. bc_radiation_cs_l, & 2859 bc_dirichlet_cs_n .OR. bc_radiation_cs_n, & 2860 bc_dirichlet_cs_r .OR. bc_radiation_cs_r, & 2861 bc_dirichlet_cs_s .OR. bc_radiation_cs_s ) 2838 2862 ELSE 2839 2863 CALL advec_s_pw( chem_species(ilsp)%conc ) … … 2942 2966 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2943 2967 IF ( ws_scheme_sca ) THEN 2944 CALL advec_s_ws( cs_advc_flags_s, i, j, chem_species(ilsp)%conc, 'kc', & 2968 CALL advec_s_ws( cs_advc_flags_s, & 2969 i, & 2970 j, & 2971 chem_species(ilsp)%conc, & 2972 'kc', & 2945 2973 chem_species(ilsp)%flux_s_cs, & 2946 2974 chem_species(ilsp)%diss_s_cs, & … … 2949 2977 i_omp_start, & 2950 2978 tn, & 2951 bc_dirichlet_l .OR. bc_radiation_l .OR. ( decycle_chem_lr .AND. nxl == 0 ),&2952 bc_dirichlet_n .OR. bc_radiation_n .OR. ( decycle_chem_ns .AND. nyn == ny ),&2953 bc_dirichlet_r .OR. bc_radiation_r .OR. ( decycle_chem_lr .AND. nxr == nx ),&2954 bc_dirichlet_s .OR. bc_radiation_s .OR. ( decycle_chem_ns .AND. nys == 0 ),&2955 monotonic_limiter_z )2979 bc_dirichlet_cs_l .OR. bc_radiation_cs_l, & 2980 bc_dirichlet_cs_n .OR. bc_radiation_cs_n, & 2981 bc_dirichlet_cs_r .OR. bc_radiation_cs_r, & 2982 bc_dirichlet_cs_s .OR. bc_radiation_cs_s, & 2983 monotonic_limiter_z ) 2956 2984 ELSE 2957 2985 CALL advec_s_pw( i, j, chem_species(ilsp)%conc ) -
palm/trunk/SOURCE/time_integration.f90
r4508 r4511 25 25 ! ----------------- 26 26 ! $Id$ 27 ! chemistry decycling replaced by explicit setting of lateral boundary conditions 28 ! 29 ! 4508 2020-04-24 13:32:20Z raasch 27 30 ! salsa decycling replaced by explicit setting of lateral boundary conditions 28 31 ! … … 242 245 243 246 USE chem_modules, & 244 ONLY: bc_cs_t_val, chem_species, emissions_anthropogenic, emiss_read_legacy_mode, & 245 n_matched_vars 246 247 #if defined( __parallel ) 248 USE chem_modules, & 249 ONLY: cs_name 250 #endif 251 252 USE chemistry_model_mod, & 253 ONLY: chem_boundary_conds 247 ONLY: bc_cs_t_val, chem_species, communicator_chem, emissions_anthropogenic, & 248 emiss_read_legacy_mode, n_matched_vars 254 249 255 250 USE control_parameters, & … … 469 464 INTEGER(iwp) :: icc !< additional index for aerosol mass bins 470 465 INTEGER(iwp) :: ig !< index for salsa gases 471 INTEGER(iwp) :: lsp !<472 466 INTEGER(iwp) :: mid !< masked output running index 473 #if defined( __parallel )474 INTEGER(iwp) :: lsp_usr !<475 467 INTEGER(iwp) :: n !< loop counter for chemistry species 476 #endif477 468 478 469 REAL(wp) :: dt_3d_old !< temporary storage of timestep to be used for … … 675 666 bc_q_t_val = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1) 676 667 IF ( air_chemistry ) THEN 677 DO lsp= 1, nvar678 bc_cs_t_val = ( chem_species( lsp)%conc_pr_init(nzt+1)&679 - chem_species( lsp)%conc_pr_init(nzt) )&668 DO n = 1, nvar 669 bc_cs_t_val = ( chem_species(n)%conc_pr_init(nzt+1) & 670 - chem_species(n)%conc_pr_init(nzt) ) & 680 671 / dzu(nzt+1) 681 672 ENDDO … … 824 815 IF ( passive_scalar ) CALL exchange_horiz( s_p, nbgp ) 825 816 IF ( air_chemistry ) THEN 826 DO lsp = 1, nvar 827 CALL exchange_horiz( chem_species(lsp)%conc_p, nbgp ) 817 DO n = 1, nvar 818 CALL exchange_horiz( chem_species(n)%conc_p, nbgp, & 819 alternative_communicator = communicator_chem ) 828 820 ENDDO 829 821 ENDIF … … 932 924 IF ( air_chemistry ) THEN 933 925 DO n = 1, nvar 934 CALL exchange_horiz( chem_species(n)%conc, nbgp ) 926 CALL exchange_horiz( chem_species(n)%conc, nbgp, & 927 alternative_communicator = communicator_chem ) 935 928 ENDDO 936 929 ENDIF … … 960 953 !-- Set boundary conditions again after interpolation and anterpolation. 961 954 CALL pmci_boundary_conds 962 963 !964 !-- Set chemistry boundary conditions (decycling)965 IF ( air_chemistry ) THEN966 DO lsp = 1, nvar967 lsp_usr = 1968 DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )969 IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) ) THEN970 CALL chem_boundary_conds( chem_species(lsp)%conc, &971 chem_species(lsp)%conc_pr_init )972 ENDIF973 lsp_usr = lsp_usr + 1974 ENDDO975 ENDDO976 ENDIF977 955 978 956 CALL cpu_log( log_point(60), 'nesting', 'stop' )
Note: See TracChangeset
for help on using the changeset viewer.
