Changeset 3848 for palm/trunk/SOURCE
- Timestamp:
- Apr 1, 2019 3:25:48 PM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/chemistry_model_mod.f90 (modified) (233 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/chemistry_model_mod.f90
r3833 r3848 27 27 ! ----------------- 28 28 ! $Id: chemistry_model_mod.f90 3784 2019-03-05 14:16:20Z banzhafs 29 ! some formatting 30 ! 31 ! 3784 2019-03-05 14:16:20Z banzhafs 29 32 ! added cs_mech to USE chem_gasphase_mod 30 33 ! … … 260 263 ! ------------ 261 264 !> Chemistry model for PALM-4U 262 !> @todo Extend chem_species type by nspec and nvar as addititional elements 265 !> @todo Extend chem_species type by nspec and nvar as addititional elements (RF) 266 !> @todo Check possibility to reduce dimension of chem_species%conc from nspec to nvar (RF) 263 267 !> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not 264 268 !> allowed to use the chemistry model in a precursor run and additionally 265 269 !> not using it in a main run 266 !> @todo Update/clean-up todo list! (FK)267 !> @todo Add routine chem_check_parameters, add checks for inconsistent or268 !> unallowed parameter settings.269 !> CALL of chem_check_parameters from check_parameters. (FK)270 270 !> @todo Implement turbulent inflow of chem spcs in inflow_turbulence. 271 271 !> @todo Separate boundary conditions for each chem spcs to be implemented … … 277 277 !> @todo test nesting for chem spcs, was implemented by suehring (kanani) 278 278 !> @todo chemistry error messages 279 !> @todo Format this module according to PALM coding standard (see e.g. module280 !> template under http://palm.muk.uni-hannover.de/mosaik/downloads/8 or281 !> D3_coding_standard.pdf under https://palm.muk.uni-hannover.de/trac/downloads/16)282 279 ! 283 280 !------------------------------------------------------------------------------! … … 318 315 319 316 LOGICAL :: nest_chemistry = .TRUE. !< flag for nesting mode of chemical species, independent on parent or not 320 321 !- Define chemical variables322 TYPE 323 CHARACTER(LEN=15) :: name 324 CHARACTER(LEN=15) :: unit 325 REAL(kind=wp), POINTER, DIMENSION(:,:,:) :: conc 326 REAL(kind=wp), POINTER, DIMENSION(:,:,:) :: conc_av 327 REAL(kind=wp), POINTER, DIMENSION(:,:,:) :: conc_p 328 REAL(kind=wp), POINTER, DIMENSION(:,:,:) :: tconc_m 329 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:) :: cssws_av 330 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:) :: flux_s_cs 331 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:) :: diss_s_cs 332 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) :: flux_l_cs 333 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) :: diss_l_cs 334 REAL(kind=wp), ALLOCATABLE, DIMENSION(:) :: conc_pr_init 317 ! 318 !- Define chemical variables within chem_species 319 TYPE species_def 320 CHARACTER(LEN=15) :: name !< name of chemical species 321 CHARACTER(LEN=15) :: unit !< unit (ppm for gases, kg m^-3 for aerosol tracers) 322 REAL(kind=wp), POINTER, DIMENSION(:,:,:) :: conc !< concentrations of trace gases 323 REAL(kind=wp), POINTER, DIMENSION(:,:,:) :: conc_av !< averaged concentrations 324 REAL(kind=wp), POINTER, DIMENSION(:,:,:) :: conc_p !< conc at prognostic time level 325 REAL(kind=wp), POINTER, DIMENSION(:,:,:) :: tconc_m !< weighted tendency of conc for previous sub-timestep (Runge-Kutta) 326 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:) :: cssws_av !< averaged fluxes of trace gases at surface 327 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:) :: flux_s_cs !< 6th-order advective flux at south face of grid box of chemical species (='cs') 328 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:) :: diss_s_cs !< artificial numerical dissipation flux at south face of grid box of chemical species 329 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) :: flux_l_cs !< 6th-order advective flux at left face of grid box of chemical species (='cs') 330 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) :: diss_l_cs !< artificial numerical dissipation flux at left face of grid box of chemical species 331 REAL(kind=wp), ALLOCATABLE, DIMENSION(:) :: conc_pr_init !< initial profile of chemical species 335 332 END TYPE species_def 336 337 338 TYPE 339 CHARACTER(LEN=15) :: name 340 CHARACTER(LEN=15) :: unit 341 REAL(kind=wp), POINTER, DIMENSION(:,:,:) :: freq 333 ! 334 !-- Define photolysis frequencies in phot_frequen 335 TYPE photols_def 336 CHARACTER(LEN=15) :: name !< name of pgotolysis frequency 337 CHARACTER(LEN=15) :: unit !< unit (1/s) 338 REAL(kind=wp), POINTER, DIMENSION(:,:,:) :: freq !< photolysis frequency 342 339 END TYPE photols_def 343 344 345 PUBLIC species_def346 PUBLIC photols_def347 340 348 341 … … 350 343 TYPE(photols_def), ALLOCATABLE, DIMENSION(:), TARGET, PUBLIC :: phot_frequen 351 344 352 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_1 353 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_2 354 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_3 355 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_av 356 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: freq_1 357 358 INTEGER, DIMENSION(nkppctrl) :: icntrl !< Fine tuning kpp 359 REAL(kind=wp), DIMENSION(nkppctrl) :: rcntrl !< Fine tuning kpp 360 LOGICAL :: decycle_chem_lr = .FALSE. 361 LOGICAL :: decycle_chem_ns = .FALSE. 345 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_1 !< pointer for swapping of timelevels for conc 346 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_2 !< pointer for swapping of timelevels for conc 347 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_3 !< pointer for swapping of timelevels for conc 348 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_av !< averaged concentrations of chemical species 349 REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: freq_1 !< pointer for phtolysis frequncies 350 !< (only 1 timelevel required) 351 INTEGER, DIMENSION(nkppctrl) :: icntrl !< 20 integer parameters for fine tuning KPP code 352 !< (e.g. solver type) 353 REAL(kind=wp), DIMENSION(nkppctrl) :: rcntrl !< 20 real parameters for fine tuning of KPP code 354 !< (e.g starting internal timestep of solver) 355 ! 356 !-- Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not 357 !-- approproate for chemicals compounds since they may accumulate too much. 358 !-- If no proper boundary conditions from a DYNAMIC input file are available, 359 !-- de-cycling applies the initial profiles at the inflow boundaries for 360 !-- Dirichlet boundary conditions 361 LOGICAL :: decycle_chem_lr = .FALSE. !< switch for setting decycling in left-right direction 362 LOGICAL :: decycle_chem_ns = .FALSE. !< switch for setting decycling in south-norht direction 362 363 CHARACTER (LEN=20), DIMENSION(4) :: decycle_method = & 363 364 (/'dirichlet','dirichlet','dirichlet','dirichlet'/) … … 376 377 377 378 378 379 379 ! 380 !-- Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010) 380 381 ! 381 382 INTEGER(iwp), PARAMETER :: nlu_dep = 15 !< Number of DEPAC landuse classes (lu's) 382 383 INTEGER(iwp), PARAMETER :: ncmp = 10 !< Number of DEPAC gas components 383 384 INTEGER(iwp), PARAMETER :: nposp = 69 !< Number of possible species for deposition 384 385 385 ! 386 !-- DEPAC landuse classes as defined in LOTOS-EUROS model v2.1 386 387 INTEGER(iwp) :: ilu_grass = 1 387 388 INTEGER(iwp) :: ilu_arable = 2 … … 400 401 INTEGER(iwp) :: ilu_semi_natural_veg = 15 401 402 402 403 403 ! 404 !-- NH3/SO2 ratio regimes: 404 405 INTEGER(iwp), PARAMETER :: iratns_low = 1 !< low ratio NH3/SO2 405 406 INTEGER(iwp), PARAMETER :: iratns_high = 2 !< high ratio NH3/SO2 406 407 INTEGER(iwp), PARAMETER :: iratns_very_low = 3 !< very low ratio NH3/SO2 407 408 408 ! 409 !-- Default: 409 410 INTEGER, PARAMETER :: iratns_default = iratns_low 410 411 411 ! 412 !-- Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) to W m-2 412 413 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: alpha =(/ 0.009, 0.009, 0.009, 0.006, 0.006, -999., -999., 0.009, -999., & 413 414 -999., 0.009, 0.006, -999., 0.009, 0.008/)*4.57 414 415 415 ! 416 !-- Set temperatures per land use for f_temp 416 417 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: tmin = (/ 12.0, 12.0, 12.0, 0.0, 0.0, -999., -999., 12.0, -999., -999., & 417 418 12.0, 0.0, -999., 12.0, 8.0/) … … 420 421 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: tmax = (/ 40.0, 40.0, 40.0, 36.0, 35.0, -999., -999., 40.0, -999., -999., & 421 422 40.0, 35.0, -999., 40.0, 39.0 /) 422 423 423 ! 424 !-- Set f_min: 424 425 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: f_min = (/ 0.01, 0.01, 0.01, 0.1, 0.1, -999., -999., 0.01, -999., -999., 0.01, & 425 426 0.1, -999., 0.01, 0.04/) 426 427 427 428 429 428 ! 429 !-- Set maximal conductance (m/s) 430 !-- (R T/P) = 1/41000 mmol/m3 is given for 20 deg C to go from mmol O3/m2/s to m/s 430 431 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: g_max = (/ 270., 300., 300., 140., 150., -999., -999., 270., -999., -999., & 431 270., 150., -999., 300., 422./)/41000 432 433 432 270., 150., -999., 300., 422./)/41000. 433 ! 434 !-- Set max, min for vapour pressure deficit vpd 434 435 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: vpd_max = (/1.3, 0.9, 0.9, 0.5, 1.0, -999., -999., 1.3, -999., -999., 1.3, & 435 436 1.0, -999., 0.9, 2.8/) … … 437 438 3.25, -999., 2.8, 4.5/) 438 439 439 440 440 PUBLIC nest_chemistry 441 PUBLIC nspec442 441 PUBLIC nreact 443 PUBLIC nvar 444 PUBLIC spc_names 442 PUBLIC nspec !< number of gas phase chemical species including constant compound (e.g. N2) 443 PUBLIC nvar !< number of variable gas phase chemical species (nvar <= nspec) 444 PUBLIC photols_def !< type defining phot_frequen 445 PUBLIC species_def !< type defining chem_species 446 PUBLIC spc_names !< names of gas phase chemical species (come from KPP) (come from KPP) 445 447 PUBLIC spec_conc_2 446 447 ! 448 !-- Interface section 448 ! 449 !-- Interface section 449 450 INTERFACE chem_3d_data_averaging 450 451 MODULE PROCEDURE chem_3d_data_averaging … … 607 608 chem_statistics, chem_swap_timelevel, chem_wrd_local, chem_depo 608 609 609 610 611 610 CONTAINS 612 611 613 ! 614 !------------------------------------------------------------------------------! 615 ! 616 ! Description: 617 ! ------------ 618 !> Subroutine for averaging 3D data of chemical species. Due to the fact that 619 !> the averaged chem arrays are allocated in chem_init, no if-query concerning 620 !> the allocation is required (in any mode). Attention: If you just specify an 621 !> averaged output quantity in the _p3dr file during restarts the first output 622 !> includes the time between the beginning of the restart run and the first 623 !> output time (not necessarily the whole averaging_interval you have 624 !> specified in your _p3d/_p3dr file ) 625 !------------------------------------------------------------------------------! 626 612 613 !------------------------------------------------------------------------------! 614 ! Description: 615 ! ------------ 616 !> Subroutine for averaging 3D data of chemical species. Due to the fact that 617 !> the averaged chem arrays are allocated in chem_init, no if-query concerning 618 !> the allocation is required (in any mode). Attention: If you just specify an 619 !> averaged output quantity in the _p3dr file during restarts the first output 620 !> includes the time between the beginning of the restart run and the first 621 !> output time (not necessarily the whole averaging_interval you have 622 !> specified in your _p3d/_p3dr file ) 623 !------------------------------------------------------------------------------! 627 624 SUBROUTINE chem_3d_data_averaging( mode, variable ) 628 625 … … 638 635 IMPLICIT NONE 639 636 640 CHARACTER (LEN=*) :: mode !<641 CHARACTER (LEN=*) :: variable !<642 643 LOGICAL :: match_def !< flag indicating natural-type surface644 LOGICAL :: match_lsm!< flag indicating natural-type surface645 LOGICAL :: match_usm!< flag indicating urban-type surface637 CHARACTER (LEN=*) :: mode !< 638 CHARACTER (LEN=*) :: variable !< 639 640 LOGICAL :: match_def !< flag indicating default-type surface 641 LOGICAL :: match_lsm !< flag indicating natural-type surface 642 LOGICAL :: match_usm !< flag indicating urban-type surface 646 643 647 644 INTEGER(iwp) :: i !< grid index x direction … … 649 646 INTEGER(iwp) :: k !< grid index z direction 650 647 INTEGER(iwp) :: m !< running index surface type 651 INTEGER(iwp) :: lsp!< running index for chem spcs648 INTEGER(iwp) :: lsp !< running index for chem spcs 652 649 653 650 IF ( (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_') ) THEN 654 651 655 IF ( mode == 'allocate' ) THEN 656 DO lsp = 1, nspec 657 IF ( TRIM(variable(1:3)) == 'kc_' .AND. & 658 TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN 659 chem_species(lsp)%conc_av = 0.0_wp 660 ENDIF 661 ENDDO 662 663 ELSEIF ( mode == 'sum' ) THEN 664 665 DO lsp = 1, nspec 666 IF ( TRIM(variable(1:3)) == 'kc_' .AND. & 667 TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN 668 DO i = nxlg, nxrg 669 DO j = nysg, nyng 670 DO k = nzb, nzt+1 671 chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) + & 672 chem_species(lsp)%conc(k,j,i) 652 IF ( mode == 'allocate' ) THEN 653 654 DO lsp = 1, nspec 655 IF ( TRIM( variable(1:3) ) == 'kc_' .AND. & 656 TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) ) THEN 657 chem_species(lsp)%conc_av = 0.0_wp 658 ENDIF 659 ENDDO 660 661 ELSEIF ( mode == 'sum' ) THEN 662 663 DO lsp = 1, nspec 664 IF ( TRIM( variable(1:3) ) == 'kc_' .AND. & 665 TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) ) THEN 666 DO i = nxlg, nxrg 667 DO j = nysg, nyng 668 DO k = nzb, nzt+1 669 chem_species(lsp)%conc_av(k,j,i) = & 670 chem_species(lsp)%conc_av(k,j,i) + & 671 chem_species(lsp)%conc(k,j,i) 672 ENDDO 673 673 ENDDO 674 674 ENDDO 675 ENDDO 676 ELSEIF ( TRIM(variable(4:)) == TRIM('cssws*') ) THEN 677 DO i = nxl, nxr 678 DO j = nys, nyn 679 match_def = surf_def_h(0)%start_index(j,i) <= & 680 surf_def_h(0)%end_index(j,i) 681 match_lsm = surf_lsm_h%start_index(j,i) <= & 682 surf_lsm_h%end_index(j,i) 683 match_usm = surf_usm_h%start_index(j,i) <= & 684 surf_usm_h%end_index(j,i) 685 686 IF ( match_def ) THEN 687 m = surf_def_h(0)%end_index(j,i) 688 chem_species(lsp)%cssws_av(j,i) = & 689 chem_species(lsp)%cssws_av(j,i) + & 690 surf_def_h(0)%cssws(lsp,m) 691 ELSEIF ( match_lsm .AND. .NOT. match_usm ) THEN 692 m = surf_lsm_h%end_index(j,i) 693 chem_species(lsp)%cssws_av(j,i) = & 694 chem_species(lsp)%cssws_av(j,i) + & 695 surf_lsm_h%cssws(lsp,m) 696 ELSEIF ( match_usm ) THEN 697 m = surf_usm_h%end_index(j,i) 698 chem_species(lsp)%cssws_av(j,i) = & 699 chem_species(lsp)%cssws_av(j,i) + & 700 surf_usm_h%cssws(lsp,m) 701 ENDIF 702 ENDDO 703 ENDDO 704 ENDIF 705 ENDDO 706 707 ELSEIF ( mode == 'average' ) THEN 708 DO lsp = 1, nspec 709 IF ( TRIM(variable(1:3)) == 'kc_' .AND. & 710 TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN 711 DO i = nxlg, nxrg 712 DO j = nysg, nyng 713 DO k = nzb, nzt+1 714 chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 675 ELSEIF ( TRIM( variable(4:) ) == TRIM( 'cssws*' ) ) THEN 676 DO i = nxl, nxr 677 DO j = nys, nyn 678 match_def = surf_def_h(0)%start_index(j,i) <= & 679 surf_def_h(0)%end_index(j,i) 680 match_lsm = surf_lsm_h%start_index(j,i) <= & 681 surf_lsm_h%end_index(j,i) 682 match_usm = surf_usm_h%start_index(j,i) <= & 683 surf_usm_h%end_index(j,i) 684 685 IF ( match_def ) THEN 686 m = surf_def_h(0)%end_index(j,i) 687 chem_species(lsp)%cssws_av(j,i) = & 688 chem_species(lsp)%cssws_av(j,i) + & 689 surf_def_h(0)%cssws(lsp,m) 690 ELSEIF ( match_lsm .AND. .NOT. match_usm ) THEN 691 m = surf_lsm_h%end_index(j,i) 692 chem_species(lsp)%cssws_av(j,i) = & 693 chem_species(lsp)%cssws_av(j,i) + & 694 surf_lsm_h%cssws(lsp,m) 695 ELSEIF ( match_usm ) THEN 696 m = surf_usm_h%end_index(j,i) 697 chem_species(lsp)%cssws_av(j,i) = & 698 chem_species(lsp)%cssws_av(j,i) + & 699 surf_usm_h%cssws(lsp,m) 700 ENDIF 715 701 ENDDO 716 702 ENDDO 717 ENDDO 718 719 ELSEIF (TRIM(variable(4:)) == TRIM('cssws*') ) THEN 720 DO i = nxlg, nxrg 721 DO j = nysg, nyng 722 chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp ) 703 ENDIF 704 ENDDO 705 706 ELSEIF ( mode == 'average' ) THEN 707 708 DO lsp = 1, nspec 709 IF ( TRIM( variable(1:3) ) == 'kc_' .AND. & 710 TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) ) THEN 711 DO i = nxlg, nxrg 712 DO j = nysg, nyng 713 DO k = nzb, nzt+1 714 chem_species(lsp)%conc_av(k,j,i) = & 715 chem_species(lsp)%conc_av(k,j,i) / & 716 REAL( average_count_3d, KIND=wp ) 717 ENDDO 718 ENDDO 723 719 ENDDO 724 ENDDO 725 CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp ) 726 ENDIF 727 ENDDO 720 721 ELSEIF ( TRIM( variable(4:) ) == TRIM( 'cssws*' ) ) THEN 722 DO i = nxlg, nxrg 723 DO j = nysg, nyng 724 chem_species(lsp)%cssws_av(j,i) = & 725 chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp ) 726 ENDDO 727 ENDDO 728 CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp ) 729 ENDIF 730 ENDDO 731 ENDIF 728 732 729 733 ENDIF 730 734 731 ENDIF732 733 735 END SUBROUTINE chem_3d_data_averaging 734 736 735 ! 736 !------------------------------------------------------------------------------! 737 ! 737 738 !------------------------------------------------------------------------------! 738 739 ! Description: 739 740 ! ------------ 740 741 !> Subroutine to initialize and set all boundary conditions for chemical species 741 742 !------------------------------------------------------------------------------! 742 743 743 SUBROUTINE chem_boundary_conds( mode ) 744 744 … … 755 755 ONLY: bc_h 756 756 757 CHARACTER ( len=*), INTENT(IN) :: mode757 CHARACTER (LEN=*), INTENT(IN) :: mode 758 758 INTEGER(iwp) :: i !< grid index x direction. 759 759 INTEGER(iwp) :: j !< grid index y direction. … … 768 768 CASE ( 'init' ) 769 769 770 IF ( bc_cs_b == 'dirichlet' ) 770 IF ( bc_cs_b == 'dirichlet' ) THEN 771 771 ibc_cs_b = 0 772 772 ELSEIF ( bc_cs_b == 'neumann' ) THEN … … 774 774 ELSE 775 775 message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"' 776 CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 ) !<776 CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 ) 777 777 ENDIF 778 778 ! 779 !-- 780 !-- 781 IF ( bc_cs_t == 'dirichlet' ) 779 !-- Set Integer flags and check for possible erroneous settings for top 780 !-- boundary condition. 781 IF ( bc_cs_t == 'dirichlet' ) THEN 782 782 ibc_cs_t = 0 783 ELSEIF ( bc_cs_t == 'neumann' ) 783 ELSEIF ( bc_cs_t == 'neumann' ) THEN 784 784 ibc_cs_t = 1 785 785 ELSEIF ( bc_cs_t == 'initial_gradient' ) THEN 786 786 ibc_cs_t = 2 787 ELSEIF ( bc_cs_t == 'nested' ) 787 ELSEIF ( bc_cs_t == 'nested' ) THEN 788 788 ibc_cs_t = 3 789 789 ELSE … … 792 792 ENDIF 793 793 794 795 794 CASE ( 'set_bc_bottomtop' ) 795 ! 796 796 !-- Bottom boundary condtions for chemical species 797 DO lsp = 1, nspec 798 IF ( ibc_cs_b == 0 ) THEN 799 DO l = 0, 1 797 DO lsp = 1, nspec 798 IF ( ibc_cs_b == 0 ) THEN 799 DO l = 0, 1 800 ! 800 801 !-- Set index kb: For upward-facing surfaces (l=0), kb=-1, i.e. 801 802 !-- the chem_species(nspec)%conc_p value at the topography top (k-1) … … 805 806 kb = MERGE( -1, 1, l == 0 ) 806 807 !$OMP PARALLEL DO PRIVATE( i, j, k ) 807 DO m = 1, bc_h(l)%ns808 i = bc_h(l)%i(m)809 j = bc_h(l)%j(m)810 k = bc_h(l)%k(m)808 DO m = 1, bc_h(l)%ns 809 i = bc_h(l)%i(m) 810 j = bc_h(l)%j(m) 811 k = bc_h(l)%k(m) 811 812 chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc(k+kb,j,i) 812 813 ENDDO 813 814 ENDDO 814 815 815 ELSEIF ( ibc_cs_b == 1 ) THEN 816 ELSEIF ( ibc_cs_b == 1 ) THEN 817 ! 816 818 !-- in boundary_conds there is som extra loop over m here for passive tracer 817 DO l = 0, 1819 DO l = 0, 1 818 820 kb = MERGE( -1, 1, l == 0 ) 819 821 !$OMP PARALLEL DO PRIVATE( i, j, k ) … … 828 830 ENDIF 829 831 ENDDO ! end lsp loop 830 832 ! 831 833 !-- Top boundary conditions for chemical species - Should this not be done for all species? 832 834 IF ( ibc_cs_t == 0 ) THEN 833 DO lsp = 1, nspec835 DO lsp = 1, nspec 834 836 chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:) 835 837 ENDDO 836 838 ELSEIF ( ibc_cs_t == 1 ) THEN 837 DO lsp = 1, nspec839 DO lsp = 1, nspec 838 840 chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) 839 841 ENDDO 840 842 ELSEIF ( ibc_cs_t == 2 ) THEN 841 DO lsp = 1, nspec843 DO lsp = 1, nspec 842 844 chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1) 843 845 ENDDO 844 846 ENDIF 845 ! 847 846 848 CASE ( 'set_bc_lateral' ) 849 ! 847 850 !-- Lateral boundary conditions for chem species at inflow boundary 848 851 !-- are automatically set when chem_species concentration is … … 854 857 855 858 IF ( bc_radiation_s ) THEN 856 DO lsp = 1, nspec859 DO lsp = 1, nspec 857 860 chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:) 858 861 ENDDO 859 862 ELSEIF ( bc_radiation_n ) THEN 860 DO lsp = 1, nspec863 DO lsp = 1, nspec 861 864 chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:) 862 865 ENDDO 863 866 ELSEIF ( bc_radiation_l ) THEN 864 DO lsp = 1, nspec867 DO lsp = 1, nspec 865 868 chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl) 866 869 ENDDO 867 870 ELSEIF ( bc_radiation_r ) THEN 868 DO lsp = 1, nspec871 DO lsp = 1, nspec 869 872 chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr) 870 873 ENDDO … … 875 878 END SUBROUTINE chem_boundary_conds 876 879 877 ! 880 878 881 !------------------------------------------------------------------------------! 879 882 ! Description: … … 881 884 !> Boundary conditions for prognostic variables in chemistry: decycling in the 882 885 !> x-direction 883 !----------------------------------------------------------------------------- 886 !------------------------------------------------------------------------------! 884 887 SUBROUTINE chem_boundary_conds_decycle( cs_3d, cs_pr_init ) 888 ! 889 !-- Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not 890 !-- approproate for chemicals compounds since they may accumulate too much. 891 !-- If no proper boundary conditions from a DYNAMIC input file are available, 892 !-- de-cycling applies the initial profiles at the inflow boundaries for 893 !-- Dirichlet boundary conditions 885 894 886 895 IMPLICIT NONE … … 897 906 898 907 flag = 0.0_wp 899 900 908 ! 901 909 !-- Left and right boundaries 902 910 IF ( decycle_chem_lr .AND. bc_lr_cyc ) THEN … … 962 970 ENDDO 963 971 ENDIF 964 965 972 ! 966 973 !-- South and north boundaries 967 974 IF ( decycle_chem_ns .AND. bc_ns_cyc ) THEN … … 1029 1036 ENDIF 1030 1037 END SUBROUTINE chem_boundary_conds_decycle 1031 ! 1032 !------------------------------------------------------------------------------! 1033 ! 1038 1039 1040 !------------------------------------------------------------------------------! 1034 1041 ! Description: 1035 1042 ! ------------ 1036 1043 !> Subroutine for checking data output for chemical species 1037 1044 !------------------------------------------------------------------------------! 1038 1039 1045 SUBROUTINE chem_check_data_output( var, unit, i, ilen, k ) 1040 1041 1046 1042 1047 IMPLICIT NONE … … 1050 1055 INTEGER(iwp) :: k 1051 1056 1052 CHARACTER( len=16) :: spec_name1057 CHARACTER(LEN=16) :: spec_name 1053 1058 1054 1059 ! … … 1058 1063 unit = 'illegal' 1059 1064 1060 spec_name = TRIM( var(4:)) !< var 1:3 is 'kc_' or 'em_'.1061 1062 IF ( TRIM( var(1:3)) == 'em_' ) THEN1063 DO lsp=1,nspec1064 IF (TRIM( spec_name) == TRIM(chem_species(lsp)%name))THEN1065 spec_name = TRIM( var(4:) ) !< var 1:3 is 'kc_' or 'em_'. 1066 1067 IF ( TRIM( var(1:3) ) == 'em_' ) THEN 1068 DO lsp=1,nspec 1069 IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN 1065 1070 unit = 'mol m-2 s-1' 1066 1071 ENDIF 1067 1068 !--It is possible to plant PM10 and PM25 into the gasphase chemistry code1069 !--as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):1070 !--set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)1071 IF (spec_name(1:2) == 'PM') 1072 ! 1073 !-- It is possible to plant PM10 and PM25 into the gasphase chemistry code 1074 !-- as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms): 1075 !-- set unit to micrograms per m**3 for PM10 and PM25 (PM2.5) 1076 IF (spec_name(1:2) == 'PM') THEN 1072 1077 unit = 'kg m-2 s-1' 1073 1078 ENDIF … … 1076 1081 ELSE 1077 1082 1078 DO lsp=1,nspec1079 IF (TRIM( spec_name) == TRIM(chem_species(lsp)%name))THEN1083 DO lsp=1,nspec 1084 IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN 1080 1085 unit = 'ppm' 1081 1086 ENDIF … … 1084 1089 !-- as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms): 1085 1090 !-- set unit to kilograms per m**3 for PM10 and PM25 (PM2.5) 1086 IF (spec_name(1:2) == 'PM') 1091 IF (spec_name(1:2) == 'PM') THEN 1087 1092 unit = 'kg m-3' 1088 1093 ENDIF 1089 1094 ENDDO 1090 1095 1091 DO lsp=1,nphot1092 IF (TRIM( spec_name) == TRIM(phot_frequen(lsp)%name))THEN1096 DO lsp=1,nphot 1097 IF (TRIM( spec_name ) == TRIM( phot_frequen(lsp)%name ) ) THEN 1093 1098 unit = 'sec-1' 1094 1099 ENDIF … … 1099 1104 RETURN 1100 1105 END SUBROUTINE chem_check_data_output 1101 ! 1102 !------------------------------------------------------------------------------! 1103 ! 1106 1107 1108 !------------------------------------------------------------------------------! 1104 1109 ! Description: 1105 1110 ! ------------ 1106 1111 !> Subroutine for checking data output of profiles for chemistry model 1107 1112 !------------------------------------------------------------------------------! 1108 1109 1113 SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit ) 1110 1114 … … 1122 1126 CHARACTER (LEN=*) :: variable !< 1123 1127 CHARACTER (LEN=*) :: dopr_unit 1124 CHARACTER (len=16) :: spec_name1128 CHARACTER (LEN=16) :: spec_name 1125 1129 1126 1130 INTEGER(iwp) :: var_count, lsp !< 1127 1131 1128 1132 1129 spec_name = TRIM( variable(4:))1133 spec_name = TRIM( variable(4:) ) 1130 1134 1131 1135 IF ( .NOT. air_chemistry ) THEN … … 1136 1140 1137 1141 ELSE 1138 DO lsp = 1, nspec1139 IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN1142 DO lsp = 1, nspec 1143 IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN 1140 1144 cs_pr_count = cs_pr_count+1 1141 1145 cs_pr_index(cs_pr_count) = lsp 1142 1146 dopr_index(var_count) = pr_palm+cs_pr_count 1143 1147 dopr_unit = 'ppm' 1144 IF (spec_name(1:2) == 'PM') 1148 IF (spec_name(1:2) == 'PM') THEN 1145 1149 dopr_unit = 'kg m-3' 1146 1150 ENDIF … … 1153 1157 END SUBROUTINE chem_check_data_output_pr 1154 1158 1155 ! 1159 1156 1160 !------------------------------------------------------------------------------! 1157 1161 ! Description: … … 1163 1167 IMPLICIT NONE 1164 1168 1165 LOGICAL 1169 LOGICAL :: found 1166 1170 INTEGER (iwp) :: lsp_usr !< running index for user defined chem spcs 1167 1171 INTEGER (iwp) :: lsp !< running index for chem spcs. … … 1171 1175 message_string = 'Chemical reactions: ON' 1172 1176 CALL message( 'chem_check_parameters', 'CM0421', 0, 0, 0, 6, 0 ) 1173 ELSEIF ( . not. (chem_gasphase_on) )THEN1177 ELSEIF ( .NOT. (chem_gasphase_on) ) THEN 1174 1178 message_string = 'Chemical reactions: OFF' 1175 1179 CALL message( 'chem_check_parameters', 'CM0422', 0, 0, 0, 6, 0 ) 1176 1180 ENDIF 1177 1181 ! 1178 1182 !-- check for chemistry time-step 1179 1183 IF ( call_chem_at_all_substeps ) THEN 1180 1184 message_string = 'Chemistry is calculated at all meteorology time-step' 1181 1185 CALL message( 'chem_check_parameters', 'CM0423', 0, 0, 0, 6, 0 ) 1182 ELSEIF ( .not. (call_chem_at_all_substeps) ) THEN1186 ELSEIF ( .not. (call_chem_at_all_substeps) ) THEN 1183 1187 message_string = 'Sub-time-steps are skipped for chemistry time-steps' 1184 1188 CALL message( 'chem_check_parameters', 'CM0424', 0, 0, 0, 6, 0 ) 1185 1189 ENDIF 1186 1190 ! 1187 1191 !-- check for photolysis scheme 1188 1192 IF ( (photolysis_scheme /= 'simple') .AND. (photolysis_scheme /= 'constant') ) THEN … … 1190 1194 CALL message( 'chem_check_parameters', 'CM0425', 1, 2, 0, 6, 0 ) 1191 1195 ENDIF 1192 1196 ! 1193 1197 !-- check for decycling of chem species 1194 IF ( (. not. any(decycle_method == 'neumann') ) .AND. (.not. any(decycle_method == 'dirichlet') ) )THEN1198 IF ( (.NOT. any(decycle_method == 'neumann') ) .AND. (.NOT. any(decycle_method == 'dirichlet') ) ) THEN 1195 1199 message_string = 'Incorrect boundary conditions. Only neumann or ' & 1196 1200 // 'dirichlet &available for decycling chemical species ' 1197 1201 CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 ) 1198 1202 ENDIF 1203 ! 1199 1204 !-- check for chemical mechanism used 1200 1205 CALL get_mechanism_name 1201 IF ( chem_mechanism /= trim(cs_mech) ) THEN1206 IF ( chem_mechanism /= TRIM( cs_mech ) ) THEN 1202 1207 message_string = 'Incorrect chemistry mechanism selected, check spelling in namelist and/or chem_gasphase_mod' 1203 1208 CALL message( 'chem_check_parameters', 'CM0462', 1, 2, 0, 6, 0 ) 1204 1209 ENDIF 1205 1206 1207 !--------------------- 1210 ! 1208 1211 !-- chem_check_parameters is called before the array chem_species is allocated! 1209 1212 !-- temporary switch of this part of the check 1210 1213 ! RETURN !bK commented 1214 1211 1215 CALL chem_init_internal 1212 !--------------------- 1213 1216 ! 1214 1217 !-- check for initial chem species input 1215 1218 lsp_usr = 1 … … 1217 1220 DO WHILE ( cs_name (lsp_usr) /= 'novalue') 1218 1221 found = .FALSE. 1219 DO lsp = 1, nvar1220 IF ( TRIM( cs_name (lsp_usr)) == TRIM(chem_species(lsp)%name) )THEN1222 DO lsp = 1, nvar 1223 IF ( TRIM( cs_name (lsp_usr) ) == TRIM( chem_species(lsp)%name) ) THEN 1221 1224 found = .TRUE. 1222 1225 EXIT 1223 1226 ENDIF 1224 1227 ENDDO 1225 IF ( .not. found ) THEN 1226 message_string = 'Unused/incorrect input for initial surface value: ' // TRIM(cs_name(lsp_usr)) 1228 IF ( .NOT. found ) THEN 1229 message_string = 'Unused/incorrect input for initial surface value: ' // & 1230 TRIM( cs_name(lsp_usr) ) 1227 1231 CALL message( 'chem_check_parameters', 'CM0427', 1, 2, 0, 6, 0 ) 1228 1232 ENDIF 1229 1233 lsp_usr = lsp_usr + 1 1230 1234 ENDDO 1231 1235 ! 1232 1236 !-- check for surface emission flux chem species 1233 1234 1237 lsp_usr = 1 1235 1238 lsp = 1 1236 1239 DO WHILE ( surface_csflux_name (lsp_usr) /= 'novalue') 1237 1240 found = .FALSE. 1238 DO lsp = 1, nvar1239 IF ( TRIM( surface_csflux_name (lsp_usr)) == TRIM(chem_species(lsp)%name) )THEN1241 DO lsp = 1, nvar 1242 IF ( TRIM( surface_csflux_name (lsp_usr) ) == TRIM( chem_species(lsp)%name ) ) THEN 1240 1243 found = .TRUE. 1241 1244 EXIT 1242 1245 ENDIF 1243 1246 ENDDO 1244 IF ( . not. found )THEN1247 IF ( .NOT. found ) THEN 1245 1248 message_string = 'Unused/incorrect input of chemical species for surface emission fluxes: ' & 1246 // TRIM( surface_csflux_name(lsp_usr))1249 // TRIM( surface_csflux_name(lsp_usr) ) 1247 1250 CALL message( 'chem_check_parameters', 'CM0428', 1, 2, 0, 6, 0 ) 1248 1251 ENDIF … … 1252 1255 END SUBROUTINE chem_check_parameters 1253 1256 1254 ! 1255 !------------------------------------------------------------------------------! 1256 ! 1257 1258 !------------------------------------------------------------------------------! 1257 1259 ! Description: 1258 1260 ! ------------ … … 1260 1262 !> @todo: Remove "mode" from argument list, not used. 1261 1263 !------------------------------------------------------------------------------! 1262 1263 1264 SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf, & 1264 1265 two_d, nzb_do, nzt_do, fill_value ) … … 1281 1282 REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) :: local_pf 1282 1283 1283 1284 1285 CHARACTER( len=16) :: spec_name1284 ! 1285 !-- local variables. 1286 CHARACTER(LEN=16) :: spec_name 1286 1287 INTEGER(iwp) :: lsp 1287 1288 INTEGER(iwp) :: i !< grid index along x-direction … … 1289 1290 INTEGER(iwp) :: k !< grid index along z-direction 1290 1291 INTEGER(iwp) :: char_len !< length of a character string 1291 1292 1292 ! 1293 1293 !-- Next statement is to avoid compiler warnings about unused variables … … 1295 1295 1296 1296 found = .FALSE. 1297 char_len = LEN_TRIM( variable)1297 char_len = LEN_TRIM( variable ) 1298 1298 1299 1299 spec_name = TRIM( variable(4:char_len-3) ) 1300 1300 1301 DO lsp=1,nspec1302 IF (TRIM( spec_name) == TRIM(chem_species(lsp)%name) .AND.&1303 ( (variable(char_len-2:) == '_xy') .OR.&1304 (variable(char_len-2:) == '_xz') .OR.&1305 (variable(char_len-2:) == '_yz') ) ) 1301 DO lsp=1,nspec 1302 IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) .AND. & 1303 ( (variable(char_len-2:) == '_xy') .OR. & 1304 (variable(char_len-2:) == '_xz') .OR. & 1305 (variable(char_len-2:) == '_yz') ) ) THEN 1306 1306 ! 1307 1307 !-- todo: remove or replace by "CALL message" mechanism (kanani) 1308 ! IF(myid == 0) WRITE(6,*) 'Output of species ' // TRIM( variable) //&1309 ! TRIM( chem_species(lsp)%name)1310 IF (av == 0) THEN1308 ! IF(myid == 0) WRITE(6,*) 'Output of species ' // TRIM( variable ) // & 1309 ! TRIM( chem_species(lsp)%name ) 1310 IF (av == 0) THEN 1311 1311 DO i = nxl, nxr 1312 1312 DO j = nys, nyn … … 1316 1316 REAL( fill_value, KIND = wp ), & 1317 1317 BTEST( wall_flags_0(k,j,i), 0 ) ) 1318 1319 1320 1318 ENDDO 1321 1319 ENDDO … … 1334 1332 ENDDO 1335 1333 ENDIF 1336 1337 1334 grid = 'zu' 1335 found = .TRUE. 1338 1336 ENDIF 1339 1337 ENDDO … … 1343 1341 END SUBROUTINE chem_data_output_2d 1344 1342 1345 ! 1346 !------------------------------------------------------------------------------! 1347 ! 1343 1344 !------------------------------------------------------------------------------! 1348 1345 ! Description: 1349 1346 ! ------------ 1350 1347 !> Subroutine defining 3D output variables for chemical species 1351 1348 !------------------------------------------------------------------------------! 1352 1353 1349 SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do ) 1354 1355 1350 1356 1351 USE indices … … 1372 1367 REAL(wp) :: fill_value !< 1373 1368 REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf 1374 1375 1376 !-- local variables 1377 CHARACTER(len=16) :: spec_name 1369 ! 1370 !-- local variables 1371 CHARACTER(LEN=16) :: spec_name 1378 1372 INTEGER(iwp) :: i 1379 1373 INTEGER(iwp) :: j … … 1385 1379 1386 1380 found = .FALSE. 1387 IF ( .NOT. (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_' ) ) THEN1381 IF ( .NOT. (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_' ) ) THEN 1388 1382 RETURN 1389 1383 ENDIF 1390 1384 1391 spec_name = TRIM( variable(4:))1392 1393 IF ( variable(1:3) == 'em_' ) THEN1394 1395 local_pf = 0.0_wp1396 1397 DOlsp = 1, nvar !!! cssws - nvar species, chem_species - nspec species !!!1398 IF ( TRIM(spec_name) == TRIM(chem_species(lsp)%name) )THEN1399 1400 !--no average for now1401 DOm = 1, surf_usm_h%ns1402 local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) = &1403 local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) + surf_usm_h%cssws(lsp,m)1404 ENDDO1405 DOm = 1, surf_lsm_h%ns1406 local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) = &1385 spec_name = TRIM( variable(4:) ) 1386 1387 IF ( variable(1:3) == 'em_' ) THEN 1388 1389 local_pf = 0.0_wp 1390 1391 DO lsp = 1, nvar !!! cssws - nvar species, chem_species - nspec species !!! 1392 IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) ) THEN 1393 ! 1394 !-- no average for now 1395 DO m = 1, surf_usm_h%ns 1396 local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) = & 1397 local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) + surf_usm_h%cssws(lsp,m) 1398 ENDDO 1399 DO m = 1, surf_lsm_h%ns 1400 local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) = & 1407 1401 local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) + surf_lsm_h%cssws(lsp,m) 1408 ENDDO1409 DOl = 0, 31410 DOm = 1, surf_usm_v(l)%ns1411 local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) = &1402 ENDDO 1403 DO l = 0, 3 1404 DO m = 1, surf_usm_v(l)%ns 1405 local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) = & 1412 1406 local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) + surf_usm_v(l)%cssws(lsp,m) 1413 ENDDO1414 DOm = 1, surf_lsm_v(l)%ns1415 local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) = &1416 local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) + surf_lsm_v(l)%cssws(lsp,m)1417 ENDDO1418 ENDDO1419 found = .TRUE.1420 ENDIF1421 ENDDO1407 ENDDO 1408 DO m = 1, surf_lsm_v(l)%ns 1409 local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) = & 1410 local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) + surf_lsm_v(l)%cssws(lsp,m) 1411 ENDDO 1412 ENDDO 1413 found = .TRUE. 1414 ENDIF 1415 ENDDO 1422 1416 ELSE 1423 DO lsp=1,nspec1424 IF (TRIM( spec_name) == TRIM(chem_species(lsp)%name))THEN1417 DO lsp = 1, nspec 1418 IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) ) THEN 1425 1419 ! 1426 1420 !-- todo: remove or replace by "CALL message" mechanism (kanani) 1427 ! IF(myid == 0 .AND. chem_debug0 ) WRITE(6,*) 'Output of species ' // TRIM( variable) //&1428 ! TRIM( chem_species(lsp)%name)1429 IF (av == 0) THEN1421 ! IF(myid == 0 .AND. chem_debug0 ) WRITE(6,*) 'Output of species ' // TRIM( variable ) // & 1422 ! TRIM( chem_species(lsp)%name ) 1423 IF (av == 0) THEN 1430 1424 DO i = nxl, nxr 1431 1425 DO j = nys, nyn … … 1440 1434 1441 1435 ELSE 1436 1442 1437 DO i = nxl, nxr 1443 1438 DO j = nys, nyn … … 1459 1454 1460 1455 END SUBROUTINE chem_data_output_3d 1461 ! 1462 !------------------------------------------------------------------------------! 1463 ! 1456 1457 1458 !------------------------------------------------------------------------------! 1464 1459 ! Description: 1465 1460 ! ------------ 1466 1461 !> Subroutine defining mask output variables for chemical species 1467 1462 !------------------------------------------------------------------------------! 1468 1469 1463 SUBROUTINE chem_data_output_mask( av, variable, found, local_pf ) 1470 1464 … … 1477 1471 IMPLICIT NONE 1478 1472 1479 CHARACTER(LEN=5) :: grid !< flag to distinquish between staggered grids 1480 1481 CHARACTER (LEN=*):: variable !< 1482 INTEGER(iwp) :: av !< flag to control data output of instantaneous or time-averaged data 1483 LOGICAL :: found !< 1473 CHARACTER(LEN=5) :: grid !< flag to distinquish between staggered grids 1474 CHARACTER(LEN=*) :: variable !< 1475 INTEGER(iwp) :: av !< flag to control data output of instantaneous or time-averaged data 1476 LOGICAL :: found 1484 1477 REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: & 1485 1478 local_pf !< 1486 1487 1488 !-- local variables. 1489 CHARACTER(len=16) :: spec_name 1479 ! 1480 !-- local variables. 1481 CHARACTER(LEN=16) :: spec_name 1490 1482 INTEGER(iwp) :: lsp 1491 1483 INTEGER(iwp) :: i !< grid index along x-direction … … 1495 1487 1496 1488 found = .TRUE. 1497 grid 1489 grid = 's' 1498 1490 1499 1491 spec_name = TRIM( variable(4:) ) 1500 1492 1501 DO lsp=1,nspec1502 IF (TRIM( spec_name) == TRIM(chem_species(lsp)%name) )THEN1493 DO lsp=1,nspec 1494 IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) ) THEN 1503 1495 ! 1504 1496 !-- todo: remove or replace by "CALL message" mechanism (kanani) 1505 ! IF(myid == 0 .AND. chem_debug0 ) WRITE(6,*) 'Output of species ' // TRIM( variable) // &1506 ! TRIM( chem_species(lsp)%name)1507 IF (av == 0) THEN1497 ! IF(myid == 0 .AND. chem_debug0 ) WRITE(6,*) 'Output of species ' // TRIM( variable ) // & 1498 ! TRIM( chem_species(lsp)%name ) 1499 IF (av == 0) THEN 1508 1500 IF ( .NOT. mask_surface(mid) ) THEN 1509 1501 … … 1582 1574 ENDIF 1583 1575 1584 1585 1576 ENDIF 1586 1577 found = .FALSE. … … 1592 1583 END SUBROUTINE chem_data_output_mask 1593 1584 1594 ! 1595 !------------------------------------------------------------------------------! 1596 ! 1585 1586 !------------------------------------------------------------------------------! 1597 1587 ! Description: 1598 1588 ! ------------ … … 1612 1602 found = .TRUE. 1613 1603 1614 IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' ) 1604 IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' ) THEN !< always the same grid for chemistry variables 1615 1605 grid_x = 'x' 1616 1606 grid_y = 'y' … … 1625 1615 1626 1616 END SUBROUTINE chem_define_netcdf_grid 1627 ! 1628 !------------------------------------------------------------------------------! 1629 ! 1617 1618 1619 !------------------------------------------------------------------------------! 1630 1620 ! Description: 1631 1621 ! ------------ … … 1637 1627 1638 1628 INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file 1639 INTEGER(iwp) :: lsp!< running index for chem spcs1640 INTEGER(iwp) 1641 CHARACTER (LEN=80) 1642 CHARACTER (LEN=80) 1629 INTEGER(iwp) :: lsp !< running index for chem spcs 1630 INTEGER(iwp) :: cs_fixed 1631 CHARACTER (LEN=80) :: docsflux_chr 1632 CHARACTER (LEN=80) :: docsinit_chr 1643 1633 ! 1644 1634 ! Get name of chemical mechanism from chem_gasphase_mod 1645 1635 CALL get_mechanism_name 1646 1636 ! 1647 1637 !-- Write chemistry model header 1648 1638 WRITE( io, 1 ) 1649 1639 ! 1650 1640 !-- Gasphase reaction status 1651 1641 IF ( chem_gasphase_on ) THEN … … 1654 1644 WRITE( io, 3 ) 1655 1645 ENDIF 1656 1657 1646 ! 1658 1647 !-- Chemistry time-step 1659 1648 WRITE ( io, 4 ) cs_time_step 1660 1649 ! 1661 1650 !-- Emission mode info 1662 1651 IF ( mode_emis == "DEFAULT" ) THEN … … 1667 1656 WRITE( io, 7 ) 1668 1657 ENDIF 1669 1658 ! 1670 1659 !-- Photolysis scheme info 1671 IF ( photolysis_scheme == "simple" ) 1660 IF ( photolysis_scheme == "simple" ) THEN 1672 1661 WRITE( io, 8 ) 1673 ELSEIF (photolysis_scheme == "constant" ) THEN1662 ELSEIF (photolysis_scheme == "constant" ) THEN 1674 1663 WRITE( io, 9 ) 1675 1664 ENDIF 1676 1665 ! 1677 1666 !-- Emission flux info 1678 1667 lsp = 1 … … 1681 1670 docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ',' 1682 1671 IF ( LEN_TRIM( docsflux_chr ) >= 75 ) THEN 1683 WRITE ( io, 10 ) docsflux_chr1684 docsflux_chr = ' '1672 WRITE ( io, 10 ) docsflux_chr 1673 docsflux_chr = ' ' 1685 1674 ENDIF 1686 1675 lsp = lsp + 1 … … 1690 1679 WRITE ( io, 10 ) docsflux_chr 1691 1680 ENDIF 1692 1693 1681 ! 1694 1682 !-- initializatoin of Surface and profile chemical species 1695 1683 … … 1708 1696 WRITE ( io, 11 ) docsinit_chr 1709 1697 ENDIF 1710 1698 ! 1711 1699 !-- number of variable and fix chemical species and number of reactions 1712 1700 cs_fixed = nspec - nvar … … 1734 1722 END SUBROUTINE chem_header 1735 1723 1736 ! 1737 !------------------------------------------------------------------------------! 1738 ! 1724 1725 !------------------------------------------------------------------------------! 1739 1726 ! Description: 1740 1727 ! ------------ … … 1742 1729 !------------------------------------------------------------------------------! 1743 1730 SUBROUTINE chem_init_arrays 1744 1731 ! 1745 1732 !-- Please use this place to allocate required arrays 1746 1733 1747 1734 END SUBROUTINE chem_init_arrays 1748 1735 1749 ! 1750 !------------------------------------------------------------------------------! 1751 ! 1736 1737 !------------------------------------------------------------------------------! 1752 1738 ! Description: 1753 1739 ! ------------ … … 1767 1753 INTEGER(iwp) :: j !< running index y dimension 1768 1754 INTEGER(iwp) :: n !< running index for chemical species 1769 1770 1755 ! 1771 1756 !-- Next statement is to avoid compiler warning about unused variables … … 1774 1759 ilu_urban ) == 0 ) CONTINUE 1775 1760 1776 1777 1778 1761 IF ( emissions_anthropogenic ) CALL chem_emissions_init 1779 1762 ! … … 1796 1779 END SUBROUTINE chem_init 1797 1780 1798 ! 1799 !------------------------------------------------------------------------------! 1800 ! 1781 1782 !------------------------------------------------------------------------------! 1801 1783 ! Description: 1802 1784 ! ------------ … … 1813 1795 1814 1796 IMPLICIT NONE 1815 1816 1797 ! 1817 1798 !-- Local variables … … 1821 1802 INTEGER(iwp) :: lpr_lev !< running index for chem spcs profile level 1822 1803 1823 IF ( emissions_anthropogenic ) THEN1804 IF ( emissions_anthropogenic ) THEN 1824 1805 CALL netcdf_data_input_chemistry_data( chem_emis_att, chem_emis ) 1825 1806 ENDIF … … 1874 1855 1875 1856 ENDDO 1876 1877 1857 ! 1878 1858 !-- Initial concentration of profiles is prescribed by parameters cs_profile 1879 1859 !-- and cs_heights in the namelist &chemistry_parameters 1860 1880 1861 CALL chem_init_profiles 1881 1862 ! … … 1892 1873 ENDDO 1893 1874 ENDIF 1894 1895 1896 1897 1875 ! 1898 1876 !-- Initialize model variables 1899 1877 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. & 1900 1878 TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN 1901 1902 1879 ! 1903 1880 !-- First model run of a possible job queue. 1904 1881 !-- Initial profiles of the variables must be computed. … … 1907 1884 ! 1908 1885 !-- Transfer initial profiles to the arrays of the 3D model 1909 DO lsp = 1, nspec1886 DO lsp = 1, nspec 1910 1887 DO i = nxlg, nxrg 1911 1888 DO j = nysg, nyng … … 1921 1898 CALL location_message( 'initializing with constant chemistry profiles', .FALSE. ) 1922 1899 1923 DO lsp = 1, nspec1900 DO lsp = 1, nspec 1924 1901 DO i = nxlg, nxrg 1925 1902 DO j = nysg, nyng … … 1930 1907 1931 1908 ENDIF 1932 1933 1909 ! 1934 1910 !-- If required, change the surface chem spcs at the start of the 3D run 1935 IF ( cs_surface_initial_change(1) /= 0.0_wp ) THEN1936 DO lsp = 1, nspec1911 IF ( cs_surface_initial_change(1) /= 0.0_wp ) THEN 1912 DO lsp = 1, nspec 1937 1913 chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) + & 1938 1914 cs_surface_initial_change(lsp) … … 1941 1917 ! 1942 1918 !-- Initiale old and new time levels. 1943 DO lsp = 1, nvar1919 DO lsp = 1, nvar 1944 1920 chem_species(lsp)%tconc_m = 0.0_wp 1945 1921 chem_species(lsp)%conc_p = chem_species(lsp)%conc … … 1948 1924 ENDIF 1949 1925 1950 1951 1952 !--- new code add above this line 1953 DO lsp = 1, nphot 1926 DO lsp = 1, nphot 1954 1927 phot_frequen(lsp)%name = phot_names(lsp) 1955 1928 ! 1956 1929 !-- todo: remove or replace by "CALL message" mechanism (kanani) 1957 ! 1958 ! WRITE(6,'(a,i4,3x,a)') 'Photolysis: ',lsp,TRIM(phot_names(lsp))1959 ! 1960 phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg) =>freq_1(:,:,:,lsp)1930 !-- IF( myid == 0 ) THEN 1931 !-- WRITE(6,'(a,i4,3x,a)') 'Photolysis: ',lsp,TRIM( phot_names(lsp) ) 1932 !-- ENDIF 1933 phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => freq_1(:,:,:,lsp) 1961 1934 ENDDO 1962 1935 … … 1967 1940 END SUBROUTINE chem_init_internal 1968 1941 1969 ! 1970 !------------------------------------------------------------------------------! 1971 ! 1942 1943 !------------------------------------------------------------------------------! 1972 1944 ! Description: 1973 1945 ! ------------ … … 1976 1948 !> analogue to parameters u_profile, v_profile and uv_heights) 1977 1949 !------------------------------------------------------------------------------! 1978 SUBROUTINE chem_init_profiles !< SUBROUTINE is called from chem_init in case of 1979 !< TRIM( initializing_actions ) /= 'read_restart_data' 1980 !< We still need to see what has to be done in case of restart run 1950 SUBROUTINE chem_init_profiles 1951 ! 1952 !-- SUBROUTINE is called from chem_init in case of TRIM( initializing_actions ) /= 'read_restart_data' 1953 !< We still need to see what has to be done in case of restart run 1954 1981 1955 USE chem_modules 1982 1956 1983 1957 IMPLICIT NONE 1984 1985 1958 ! 1959 !-- Local variables 1986 1960 INTEGER :: lsp !< running index for number of species in derived data type species_def 1987 1961 INTEGER :: lsp_usr !< running index for number of species (user defined) in cs_names, cs_profiles etc 1988 1962 INTEGER :: lpr_lev !< running index for profile level for each chem spcs. 1989 1963 INTEGER :: npr_lev !< the next available profile lev 1990 1991 !----------------- 1992 !-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles 1993 !-- and heights. If parameter "cs_profile" is not prescribed then initial surface values 1994 !-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and 1995 !-- "cs_heights" are prescribed, their values will!override the constant profile given by 1996 !-- "cs_surface". 1997 ! 1964 ! 1965 !-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles 1966 !-- and heights. If parameter "cs_profile" is not prescribed then initial surface values 1967 !-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and 1968 !-- "cs_heights" are prescribed, their values will!override the constant profile given by 1969 !-- "cs_surface". 1998 1970 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 1999 1971 lsp_usr = 1 2000 1972 DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' ) !'novalue' is the default 2001 1973 DO lsp = 1, nspec ! 2002 !-- create initial profile (conc_pr_init) for each chemical species 1974 ! 1975 !-- create initial profile (conc_pr_init) for each chemical species 2003 1976 IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) ) THEN ! 2004 IF ( cs_profile(lsp_usr,1) == 9999999.9_wp ) THEN 2005 !-- set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species 1977 IF ( cs_profile(lsp_usr,1) == 9999999.9_wp ) THEN 1978 ! 1979 !-- set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species 2006 1980 DO lpr_lev = 0, nzt+1 2007 1981 chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr) … … 2016 1990 2017 1991 npr_lev = 1 2018 !chem_species(lsp)%conc_pr_init(0) = 0.0_wp1992 ! chem_species(lsp)%conc_pr_init(0) = 0.0_wp 2019 1993 DO lpr_lev = 1, nz+1 2020 1994 IF ( npr_lev < 100 ) THEN … … 2028 2002 ENDDO 2029 2003 ENDIF 2030 IF ( npr_lev < 100 .AND. cs_heights(lsp_usr, npr_lev +1) /= 9999999.9_wp ) THEN2004 IF ( npr_lev < 100 .AND. cs_heights(lsp_usr,npr_lev+1) /= 9999999.9_wp ) THEN 2031 2005 chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) + & 2032 2006 ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) / & … … 2038 2012 ENDDO 2039 2013 ENDIF 2040 !-- If a profile is prescribed explicity using cs_profiles and cs_heights, then 2041 !-- chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based 2042 !-- on the cs_profiles(lsp_usr,:) and cs_heights(lsp_usr,:). 2014 ! 2015 !-- If a profile is prescribed explicity using cs_profiles and cs_heights, then 2016 !-- chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based 2017 !-- on the cs_profiles(lsp_usr,:) and cs_heights(lsp_usr,:). 2043 2018 ENDIF 2044 2019 ENDDO … … 2048 2023 2049 2024 END SUBROUTINE chem_init_profiles 2050 ! 2051 !------------------------------------------------------------------------------! 2052 ! 2053 ! Description: 2054 ! ------------ 2055 !> Subroutine to integrate chemical species in the given chemical mechanism 2056 !------------------------------------------------------------------------------! 2057 2025 2026 2027 !------------------------------------------------------------------------------! 2028 ! Description: 2029 ! ------------ 2030 !> Subroutine to integrate chemical species in the given chemical mechanism 2031 !------------------------------------------------------------------------------! 2058 2032 SUBROUTINE chem_integrate_ij( i, j ) 2059 2033 … … 2068 2042 INTEGER,INTENT(IN) :: i 2069 2043 INTEGER,INTENT(IN) :: j 2070 2071 2044 ! 2045 !-- local variables 2072 2046 INTEGER(iwp) :: lsp !< running index for chem spcs. 2073 2047 INTEGER(iwp) :: lph !< running index for photolysis frequencies … … 2096 2070 REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local 2097 2071 2098 2099 2072 REAL(kind=wp) :: dt_chem 2100 2101 2073 ! 2102 2074 !-- Set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry 2103 IF (chem_gasphase_on) THEN2075 IF (chem_gasphase_on) THEN 2104 2076 nacc = 0 2105 2077 nrej = 0 2106 2078 2107 2079 tmp_temp(:) = pt(nzb+1:nzt,j,i) * exner(nzb+1:nzt) 2108 ! 2109 !-- ppm to molecules/cm**3 2110 !-- tmp_fact = 1.e-6_wp*6.022e23_wp/(22.414_wp*1000._wp) * 273.15_wp * hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp ) 2080 ! 2081 !-- convert ppm to molecules/cm**3 2082 !-- tmp_fact = 1.e-6_wp*6.022e23_wp/(22.414_wp*1000._wp) * 273.15_wp * 2083 !-- hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp ) 2111 2084 conv = ppm2fr * xna / vmolcm 2112 2085 tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std) 2113 2086 tmp_fact_i = 1.0_wp/tmp_fact 2114 2087 2115 IF ( humidity ) THEN2088 IF ( humidity ) THEN 2116 2089 IF ( bulk_cloud_model ) THEN 2117 2090 tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) * & … … 2124 2097 ENDIF 2125 2098 2126 DO lsp = 1,nspec2099 DO lsp = 1,nspec 2127 2100 tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:) 2128 2101 ENDDO … … 2131 2104 tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i) 2132 2105 ENDDO 2133 ! 2134 !-- todo: remove (kanani) 2135 ! IF(myid == 0 .AND. chem_debug0 ) THEN 2136 ! IF (i == 10 .AND. j == 10) WRITE(0,*) 'begin chemics step ',dt_3d 2137 ! ENDIF 2138 2139 !-- Compute length of time step 2106 ! 2107 !-- Compute length of time step 2140 2108 IF ( call_chem_at_all_substeps ) THEN 2141 2109 dt_chem = dt_3d * weight_pres(intermediate_timestep_count) … … 2146 2114 cs_time_step = dt_chem 2147 2115 2148 IF(maxval(rcntrl) > 0.0) 2116 IF(maxval(rcntrl) > 0.0) THEN ! Only if rcntrl is set 2149 2117 IF( time_since_reference_point <= 2*dt_3d) THEN 2150 2118 rcntrl_local = 0 … … 2159 2127 icntrl_i = icntrl, rcntrl_i = rcntrl_local, xnacc = nacc, xnrej = nrej, istatus=istatus ) 2160 2128 2161 DO lsp = 1,nspec2129 DO lsp = 1,nspec 2162 2130 chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:) 2163 2131 ENDDO … … 2168 2136 RETURN 2169 2137 END SUBROUTINE chem_integrate_ij 2170 ! 2171 !------------------------------------------------------------------------------! 2172 2173 2174 2175 2176 2138 2139 2140 !------------------------------------------------------------------------------! 2141 ! Description: 2142 ! ------------ 2143 !> Subroutine defining parin for &chemistry_parameters for chemistry model 2144 !------------------------------------------------------------------------------! 2177 2145 SUBROUTINE chem_parin 2178 2146 … … 2230 2198 surface_csflux_name, & 2231 2199 time_fac_type 2232 2233 2234 2200 ! 2201 !-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj) 2202 !-- so this way we could prescribe a specific flux value for each species 2235 2203 !> chemistry_parameters for initial profiles 2236 2204 !> cs_names = 'O3', 'NO2', 'NO', ... to set initial profiles) … … 2240 2208 !> If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)" 2241 2209 !> then write these cs_surface values to chem_species(lsp)%conc_pr_init(:) 2242 2243 ! 2244 !-- Read chem namelist 2210 ! 2211 !-- Read chem namelist 2245 2212 2246 2213 CHARACTER(LEN=8) :: solver_type … … 2252 2219 atol = 1.0_wp 2253 2220 rtol = 0.01_wp 2254 2255 2221 ! 2222 !-- Try to find chemistry package 2256 2223 REWIND ( 11 ) 2257 2224 line = ' ' … … 2260 2227 ENDDO 2261 2228 BACKSPACE ( 11 ) 2262 2263 2229 ! 2230 !-- Read chemistry namelist 2264 2231 READ ( 11, chemistry_parameters, ERR = 10, END = 20 ) 2265 2266 2232 ! 2233 !-- Enable chemistry model 2267 2234 air_chemistry = .TRUE. 2268 2235 GOTO 20 … … 2273 2240 2274 2241 20 CONTINUE 2275 2276 ! 2277 !-- check for emission mode for chem species 2242 ! 2243 !-- check for emission mode for chem species 2278 2244 IF ( (mode_emis /= 'PARAMETERIZED') .AND. ( mode_emis /= 'DEFAULT' ) .AND. ( mode_emis /= 'PRE-PROCESSED' ) ) THEN 2279 2245 message_string = 'Incorrect mode_emiss option select. Please check spelling' … … 2282 2248 2283 2249 t_steps = my_steps 2284 2285 2286 2250 ! 2251 !-- Determine the number of user-defined profiles and append them to the 2252 !-- standard data output (data_output_pr) 2287 2253 max_pr_cs_tmp = 0 2288 2254 i = 1 2289 2255 2290 2256 DO WHILE ( data_output_pr(i) /= ' ' .AND. i <= 100 ) 2291 IF ( TRIM( data_output_pr(i)(1:3) ) == 'kc_' )THEN2257 IF ( TRIM( data_output_pr(i)(1:3) ) == 'kc_' ) THEN 2292 2258 max_pr_cs_tmp = max_pr_cs_tmp+1 2293 2259 ENDIF … … 2295 2261 ENDDO 2296 2262 2297 IF ( max_pr_cs_tmp > 0 ) THEN2263 IF ( max_pr_cs_tmp > 0 ) THEN 2298 2264 cs_pr_namelist_found = .TRUE. 2299 2265 max_pr_cs = max_pr_cs_tmp … … 2301 2267 2302 2268 ! Set Solver Type 2303 IF(icntrl(3) == 0) 2269 IF(icntrl(3) == 0) THEN 2304 2270 solver_type = 'rodas3' !Default 2305 ELSE IF(icntrl(3) == 1) 2271 ELSE IF(icntrl(3) == 1) THEN 2306 2272 solver_type = 'ros2' 2307 ELSE IF(icntrl(3) == 2) 2273 ELSE IF(icntrl(3) == 2) THEN 2308 2274 solver_type = 'ros3' 2309 ELSE IF(icntrl(3) == 3) 2275 ELSE IF(icntrl(3) == 3) THEN 2310 2276 solver_type = 'ro4' 2311 ELSE IF(icntrl(3) == 4) 2277 ELSE IF(icntrl(3) == 4) THEN 2312 2278 solver_type = 'rodas3' 2313 ELSE IF(icntrl(3) == 5) 2279 ELSE IF(icntrl(3) == 5) THEN 2314 2280 solver_type = 'rodas4' 2315 ELSE IF(icntrl(3) == 6) 2281 ELSE IF(icntrl(3) == 6) THEN 2316 2282 solver_type = 'Rang3' 2317 2283 ELSE … … 2320 2286 END IF 2321 2287 2322 2323 2324 ! write(text,*) 'gas_phase chemistry: solver_type = ',TRIM(solver_type)2325 2326 !kk CALL location_message( TRIM(text), .FALSE. )2327 ! IF(myid == 0)THEN2328 2329 2330 2331 ! write(9,*) ' gas_phase chemistry: solver_type = ',TRIM(solver_type)2332 2333 2334 2335 2336 2337 ! IF(vl_dim > 1)THEN2338 2339 2340 2341 2342 2343 2288 ! 2289 !-- todo: remove or replace by "CALL message" mechanism (kanani) 2290 ! write(text,*) 'gas_phase chemistry: solver_type = ',TRIM( solver_type ) 2291 !kk Has to be changed to right calling sequence 2292 !kk CALL location_message( TRIM( text ), .FALSE. ) 2293 ! IF(myid == 0) THEN 2294 ! write(9,*) ' ' 2295 ! write(9,*) 'kpp setup ' 2296 ! write(9,*) ' ' 2297 ! write(9,*) ' gas_phase chemistry: solver_type = ',TRIM( solver_type ) 2298 ! write(9,*) ' ' 2299 ! write(9,*) ' Hstart = ',rcntrl(3) 2300 ! write(9,*) ' FacMin = ',rcntrl(4) 2301 ! write(9,*) ' FacMax = ',rcntrl(5) 2302 ! write(9,*) ' ' 2303 ! IF(vl_dim > 1) THEN 2304 ! write(9,*) ' Vector mode vektor length = ',vl_dim 2305 ! ELSE 2306 ! write(9,*) ' Scalar mode' 2307 ! ENDIF 2308 ! write(9,*) ' ' 2309 ! END IF 2344 2310 2345 2311 RETURN … … 2347 2313 END SUBROUTINE chem_parin 2348 2314 2349 ! 2350 !------------------------------------------------------------------------------! 2351 ! 2352 ! Description: 2353 ! ------------ 2354 !> Subroutine calculating prognostic equations for chemical species 2355 !> (vector-optimized). 2356 !> Routine is called separately for each chemical species over a loop from 2357 !> prognostic_equations. 2358 !------------------------------------------------------------------------------! 2315 2316 !------------------------------------------------------------------------------! 2317 ! Description: 2318 ! ------------ 2319 !> Subroutine calculating prognostic equations for chemical species 2320 !> (vector-optimized). 2321 !> Routine is called separately for each chemical species over a loop from 2322 !> prognostic_equations. 2323 !------------------------------------------------------------------------------! 2359 2324 SUBROUTINE chem_prognostic_equations( cs_scalar_p, cs_scalar, tcs_scalar_m, & 2360 2325 pr_init_cs, ilsp ) … … 2390 2355 2391 2356 2392 2393 2357 ! 2358 !-- Tendency terms for chemical species 2394 2359 tend = 0.0_wp 2395 2396 2360 ! 2361 !-- Advection terms 2397 2362 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2398 2363 IF ( ws_scheme_sca ) THEN … … 2404 2369 CALL advec_s_up( cs_scalar ) 2405 2370 ENDIF 2406 2407 2371 ! 2372 !-- Diffusion terms (the last three arguments are zero) 2408 2373 CALL diffusion_s( cs_scalar, & 2409 2374 surf_def_h(0)%cssws(ilsp,:), & … … 2424 2389 surf_usm_v(2)%cssws(ilsp,:), & 2425 2390 surf_usm_v(3)%cssws(ilsp,:) ) 2426 2427 2391 ! 2392 !-- Prognostic equation for chemical species 2428 2393 DO i = nxl, nxr 2429 2394 DO j = nys, nyn … … 2443 2408 ENDDO 2444 2409 ENDDO 2445 2446 2410 ! 2411 !-- Calculate tendencies for the next Runge-Kutta step 2447 2412 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2448 2413 IF ( intermediate_timestep_count == 1 ) THEN … … 2469 2434 END SUBROUTINE chem_prognostic_equations 2470 2435 2471 !------------------------------------------------------------------------------! 2472 2473 2474 2475 2476 2477 2478 2479 2480 SUBROUTINE chem_prognostic_equations_ij( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs,&2481 i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs,&2436 2437 !------------------------------------------------------------------------------! 2438 ! Description: 2439 ! ------------ 2440 !> Subroutine calculating prognostic equations for chemical species 2441 !> (cache-optimized). 2442 !> Routine is called separately for each chemical species over a loop from 2443 !> prognostic_equations. 2444 !------------------------------------------------------------------------------! 2445 SUBROUTINE chem_prognostic_equations_ij( cs_scalar_p, cs_scalar, tcs_scalar_m, & 2446 pr_init_cs, i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs, & 2482 2447 flux_l_cs, diss_l_cs ) 2483 2448 USE pegrid … … 2507 2472 REAL(wp), DIMENSION(0:nz+1) :: pr_init_cs !< 2508 2473 2509 2510 2474 ! 2475 !-- local variables 2511 2476 2512 2477 INTEGER :: k 2513 2514 2478 ! 2479 !-- Tendency-terms for chem spcs. 2515 2480 tend(:,j,i) = 0.0_wp 2516 2517 2481 ! 2482 !-- Advection terms 2518 2483 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2519 2484 IF ( ws_scheme_sca ) THEN … … 2526 2491 CALL advec_s_up( i, j, cs_scalar ) 2527 2492 ENDIF 2528 2529 ! 2530 2531 !-- Diffusion terms (the last three arguments are zero) 2493 ! 2494 !-- Diffusion terms (the last three arguments are zero) 2532 2495 2533 2496 CALL diffusion_s( i, j, cs_scalar, & … … 2541 2504 surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:), & 2542 2505 surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) ) 2543 2544 ! 2545 !-- Prognostic equation for chem spcs 2546 DO k = nzb+1, nzt 2506 ! 2507 !-- Prognostic equation for chem spcs 2508 DO k = nzb+1, nzt 2547 2509 cs_scalar_p(k,j,i) = cs_scalar(k,j,i) + ( dt_3d * & 2548 2510 ( tsc(2) * tend(k,j,i) + & … … 2557 2519 IF ( cs_scalar_p(k,j,i) < 0.0_wp ) cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i) !FKS6 2558 2520 ENDDO 2559 2560 ! 2561 !-- Calculate tendencies for the next Runge-Kutta step 2521 ! 2522 !-- Calculate tendencies for the next Runge-Kutta step 2562 2523 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2563 2524 IF ( intermediate_timestep_count == 1 ) THEN … … 2576 2537 END SUBROUTINE chem_prognostic_equations_ij 2577 2538 2578 ! 2579 !------------------------------------------------------------------------------! 2580 ! 2581 ! Description: 2582 ! ------------ 2583 !> Subroutine to read restart data of chemical species 2584 !------------------------------------------------------------------------------! 2585 2539 2540 !------------------------------------------------------------------------------! 2541 ! Description: 2542 ! ------------ 2543 !> Subroutine to read restart data of chemical species 2544 !------------------------------------------------------------------------------! 2586 2545 SUBROUTINE chem_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, & 2587 2546 nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc, & … … 2623 2582 IF ( ALLOCATED(chem_species) ) THEN 2624 2583 2625 DO lsp = 1, nspec2584 DO lsp = 1, nspec 2626 2585 2627 2586 !< for time-averaged chemical conc. 2628 spc_name_av = TRIM( chem_species(lsp)%name)//'_av'2629 2630 IF ( restart_string(1:length) == TRIM( chem_species(lsp)%name) ) &2631 2587 spc_name_av = TRIM( chem_species(lsp)%name )//'_av' 2588 2589 IF ( restart_string(1:length) == TRIM( chem_species(lsp)%name) ) & 2590 THEN 2632 2591 !< read data into tmp_3d 2633 2592 IF ( k == 1 ) READ ( 13 ) tmp_3d 2634 2593 !< fill ..%conc in the restart run 2635 chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp, &2594 chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp, & 2636 2595 nxlc-nbgp:nxrc+nbgp) = & 2637 2596 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) … … 2639 2598 ELSEIF (restart_string(1:length) == spc_name_av ) THEN 2640 2599 IF ( k == 1 ) READ ( 13 ) tmp_3d 2641 chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp, &2600 chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp, & 2642 2601 nxlc-nbgp:nxrc+nbgp) = & 2643 2602 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) … … 2651 2610 2652 2611 END SUBROUTINE chem_rrd_local 2653 ! 2654 !-------------------------------------------------------------------------------! 2655 !> Description: 2656 !> Calculation of horizontally averaged profiles 2657 !> This routine is called for every statistic region (sr) defined by the user, 2658 !> but at least for the region "total domain" (sr=0). 2659 !> quantities. 2660 !------------------------------------------------------------------------------! 2612 2613 2614 !-------------------------------------------------------------------------------! 2615 !> Description: 2616 !> Calculation of horizontally averaged profiles 2617 !> This routine is called for every statistic region (sr) defined by the user, 2618 !> but at least for the region "total domain" (sr=0). 2619 !> quantities. 2620 !-------------------------------------------------------------------------------! 2661 2621 SUBROUTINE chem_statistics( mode, sr, tn ) 2662 2622 … … 2681 2641 IF ( mode == 'profiles' ) THEN 2682 2642 ! 2683 !-- Sample on how to calculate horizontally averaged profiles of user- 2684 !-- defined quantities. Each quantity is identified by the index 2685 !-- "pr_palm+#" where "#" is an integer starting from 1. These 2686 !-- user-profile-numbers must also be assigned to the respective strings 2687 !-- given by data_output_pr_cs in routine user_check_data_output_pr. 2688 !-- hom(:,:,:,:) = dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g. 2689 ! w*pt*), dim-4 = statistical region. 2690 2691 !$OMP DO 2643 ! 2644 !-- Sample on how to calculate horizontally averaged profiles of user- 2645 !-- defined quantities. Each quantity is identified by the index 2646 !-- "pr_palm+#" where "#" is an integer starting from 1. These 2647 !-- user-profile-numbers must also be assigned to the respective strings 2648 !-- given by data_output_pr_cs in routine user_check_data_output_pr. 2649 !-- hom(:,:,:,:) = dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g. 2650 !-- w*pt*), dim-4 = statistical region. 2651 2652 !$OMP DO 2692 2653 DO i = nxl, nxr 2693 2654 DO j = nys, nyn … … 2710 2671 END SUBROUTINE chem_statistics 2711 2672 2712 ! 2713 2714 ! 2715 ! Description: 2716 ! ------------ 2717 !> Subroutine for swapping of timelevels for chemical species 2718 !> called out from subroutine swap_timelevel 2719 !------------------------------------------------------------------------------! 2673 2674 !------------------------------------------------------------------------------! 2675 ! Description: 2676 ! ------------ 2677 !> Subroutine for swapping of timelevels for chemical species 2678 !> called out from subroutine swap_timelevel 2679 !------------------------------------------------------------------------------! 2680 2720 2681 2721 2682 SUBROUTINE chem_swap_timelevel( level ) … … 2724 2685 2725 2686 INTEGER(iwp), INTENT(IN) :: level 2726 2727 2687 ! 2688 !-- local variables 2728 2689 INTEGER(iwp) :: lsp 2729 2690 2730 2691 2731 IF ( level == 0 ) THEN2732 DO lsp=1, nvar2692 IF ( level == 0 ) THEN 2693 DO lsp=1, nvar 2733 2694 chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1(:,:,:,lsp) 2734 2695 chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2(:,:,:,lsp) 2735 2696 ENDDO 2736 2697 ELSE 2737 DO lsp=1, nvar2698 DO lsp=1, nvar 2738 2699 chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2(:,:,:,lsp) 2739 2700 chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1(:,:,:,lsp) … … 2743 2704 RETURN 2744 2705 END SUBROUTINE chem_swap_timelevel 2745 ! 2746 !------------------------------------------------------------------------------! 2747 2748 2749 2750 2751 2706 2707 2708 !------------------------------------------------------------------------------! 2709 ! Description: 2710 ! ------------ 2711 !> Subroutine to write restart data for chemistry model 2712 !------------------------------------------------------------------------------! 2752 2713 SUBROUTINE chem_wrd_local 2714 2753 2715 2754 2716 IMPLICIT NONE … … 2766 2728 2767 2729 2768 2769 !-------------------------------------------------------------------------------! 2770 ! 2771 ! Description: 2772 ! ------------ 2773 !> Subroutine to calculate the deposition of gases and PMs. For now deposition 2774 !> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT 2775 !> considered. The deposition of particles is derived following Zhang et al., 2776 !> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010). 2777 !> 2778 !> @TODO: Consider deposition on vertical surfaces 2779 !> @TODO: Consider overlaying horizontal surfaces 2780 !> @TODO: Consider resolved vegetation 2781 !> @TODO: Check error messages 2782 !-------------------------------------------------------------------------------! 2783 2730 !-------------------------------------------------------------------------------! 2731 ! Description: 2732 ! ------------ 2733 !> Subroutine to calculate the deposition of gases and PMs. For now deposition 2734 !> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT 2735 !> considered. The deposition of particles is derived following Zhang et al., 2736 !> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010). 2737 !> 2738 !> @TODO: Consider deposition on vertical surfaces 2739 !> @TODO: Consider overlaying horizontal surfaces 2740 !> @TODO: Consider resolved vegetation 2741 !> @TODO: Check error messages 2742 !-------------------------------------------------------------------------------! 2784 2743 SUBROUTINE chem_depo( i, j ) 2785 2744 … … 2801 2760 2802 2761 2803 2804 2762 IMPLICIT NONE 2805 2763 2806 2764 INTEGER(iwp), INTENT(IN) :: i 2807 2765 INTEGER(iwp), INTENT(IN) :: j 2808 2809 2766 INTEGER(iwp) :: k !< matching k to surface m at i,j 2810 2767 INTEGER(iwp) :: lsp !< running index for chem spcs. 2811 ! INTEGER(iwp) :: lu 2768 ! INTEGER(iwp) :: lu !< running index for landuse classes 2812 2769 INTEGER(iwp) :: lu_palm !< index of PALM LSM vegetation_type at current surface element 2813 2770 INTEGER(iwp) :: lup_palm !< index of PALM LSM pavement_type at current surface element … … 2826 2783 INTEGER(iwp) :: pspec !< running index 2827 2784 INTEGER(iwp) :: i_pspec !< index for matching depac gas component 2828 2829 ! 2830 !-- Vegetation !<Match to DEPAC 2785 ! 2786 !-- Vegetation !<Match to DEPAC 2831 2787 INTEGER(iwp) :: ind_lu_user = 0 !< --> ERROR 2832 2788 INTEGER(iwp) :: ind_lu_b_soil = 1 !< --> ilu_desert … … 2848 2804 INTEGER(iwp) :: ind_lu_mixed_forest = 17 !< --> ilu_coniferous_forest (ave(decid+conif)) 2849 2805 INTEGER(iwp) :: ind_lu_intrup_forest = 18 !< --> ilu_other (ave(other+decid)) 2850 2851 ! 2852 !-- Water 2806 ! 2807 !-- Water 2853 2808 INTEGER(iwp) :: ind_luw_user = 0 !< --> ERROR 2854 2809 INTEGER(iwp) :: ind_luw_lake = 1 !< --> ilu_water_inland … … 2857 2812 INTEGER(iwp) :: ind_luw_pond = 4 !< --> ilu_water_inland 2858 2813 INTEGER(iwp) :: ind_luw_fountain = 5 !< --> ilu_water_inland 2859 2860 ! 2861 !-- Pavement 2814 ! 2815 !-- Pavement 2862 2816 INTEGER(iwp) :: ind_lup_user = 0 !< --> ERROR 2863 2817 INTEGER(iwp) :: ind_lup_asph_conc = 1 !< --> ilu_desert … … 2876 2830 INTEGER(iwp) :: ind_lup_art_turf = 14 !< --> ilu_desert 2877 2831 INTEGER(iwp) :: ind_lup_clay = 15 !< --> ilu_desert 2878 2879 2880 ! 2881 !-- Particle parameters according to the respective aerosol classes (PM25, PM10) 2832 ! 2833 !-- Particle parameters according to the respective aerosol classes (PM25, PM10) 2882 2834 INTEGER(iwp) :: ind_p_size = 1 !< index for partsize in particle_pars 2883 2835 INTEGER(iwp) :: ind_p_dens = 2 !< index for rhopart in particle_pars … … 2929 2881 REAL(wp) :: ts !< surface temperatur in degrees celsius 2930 2882 REAL(wp) :: qv_tmp !< surface mixing ratio at current surface element 2931 2932 ! 2933 !-- Particle parameters (PM10 (1), PM25 (2)) 2934 !-- partsize (diameter in m), rhopart (density in kg/m3), slipcor 2935 !-- (slip correction factor dimensionless, Seinfeld and Pandis 2006, Table 9.3) 2883 ! 2884 !-- Particle parameters (PM10 (1), PM25 (2)) 2885 !-- partsize (diameter in m), rhopart (density in kg/m3), slipcor 2886 !-- (slip correction factor dimensionless, Seinfeld and Pandis 2006, Table 9.3) 2936 2887 REAL(wp), DIMENSION(1:3,1:2), PARAMETER :: particle_pars = RESHAPE( (/ & 2937 2888 8.0e-6_wp, 1.14e3_wp, 1.016_wp, & !< 1 … … 2939 2890 /), (/ 3, 2 /) ) 2940 2891 2941 2942 2892 LOGICAL :: match_lsm !< flag indicating natural-type surface 2943 2893 LOGICAL :: match_usm !< flag indicating urban-type surface 2944 2945 2946 ! 2947 !-- List of names of possible tracers 2948 CHARACTER(len=*), PARAMETER :: pspecnames(nposp) = (/ & 2894 ! 2895 !-- List of names of possible tracers 2896 CHARACTER(LEN=*), PARAMETER :: pspecnames(nposp) = (/ & 2949 2897 'NO2 ', & !< NO2 2950 2898 'NO ', & !< NO … … 3016 2964 'tpm25_biascorr', & !< tpm25_biascorr 3017 2965 'O3_biascorr ' /) !< o3_biascorr 3018 3019 3020 ! 3021 !-- tracer mole mass: 2966 ! 2967 !-- tracer mole mass: 3022 2968 REAL(wp), PARAMETER :: specmolm(nposp) = (/ & 3023 2969 xm_O * 2 + xm_N, & !< NO2 … … 3090 3036 xm_dummy, & !< tpm25_biascorr 3091 3037 xm_O * 3 /) !< o3_biascorr 3092 3093 3094 ! 3095 !-- Initialize surface element m 3038 ! 3039 !-- Initialize surface element m 3096 3040 m = 0 3097 3041 k = 0 3098 3099 3100 3042 ! 3043 !-- LSM or USM surface present at i,j: 3044 !-- Default surfaces are NOT considered for deposition 3101 3045 match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i) 3102 3046 match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i) 3103 3104 3105 ! 3106 !--For LSM surfaces 3047 ! 3048 !--For LSM surfaces 3107 3049 3108 3050 IF ( match_lsm ) THEN 3109 3110 ! 3111 !-- Get surface element information at i,j: 3051 ! 3052 !-- Get surface element information at i,j: 3112 3053 m = surf_lsm_h%start_index(j,i) 3113 3054 k = surf_lsm_h%k(m) 3114 3115 ! 3116 !-- Get needed variables for surface element m 3055 ! 3056 !-- Get needed variables for surface element m 3117 3057 ustar_lu = surf_lsm_h%us(m) 3118 3058 z0h_lu = surf_lsm_h%z0h(m) … … 3121 3061 lai = surf_lsm_h%lai(m) 3122 3062 sai = lai + 1 3123 3124 ! 3125 !-- For small grid spacing neglect R_a 3063 ! 3064 !-- For small grid spacing neglect R_a 3126 3065 IF ( dzw(k) <= 1.0 ) THEN 3127 3066 r_aero_lu = 0.0_wp 3128 3067 ENDIF 3129 3130 ! 3131 !--Initialize lu's 3068 ! 3069 !-- Initialize lu's 3132 3070 lu_palm = 0 3133 3071 lu_dep = 0 … … 3136 3074 luw_palm = 0 3137 3075 luw_dep = 0 3138 3139 ! 3140 !--Initialize budgets 3076 ! 3077 !-- Initialize budgets 3141 3078 bud_now_lu = 0.0_wp 3142 3079 bud_now_lup = 0.0_wp 3143 3080 bud_now_luw = 0.0_wp 3144 3145 3146 ! 3147 !-- Get land use for i,j and assign to DEPAC lu 3081 ! 3082 !-- Get land use for i,j and assign to DEPAC lu 3148 3083 IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 ) THEN 3149 3084 lu_palm = surf_lsm_h%vegetation_type(m) … … 3245 3180 ENDIF 3246 3181 ENDIF 3247 3248 3249 ! 3250 !-- Set wetness indicator to dry or wet for lsm vegetation or pavement 3182 ! 3183 !-- Set wetness indicator to dry or wet for lsm vegetation or pavement 3251 3184 IF ( surf_lsm_h%c_liq(m) > 0 ) THEN 3252 3185 nwet = 1 … … 3254 3187 nwet = 0 3255 3188 ENDIF 3256 3257 ! 3258 !--Compute length of time step 3189 ! 3190 !-- Compute length of time step 3259 3191 IF ( call_chem_at_all_substeps ) THEN 3260 3192 dt_chem = dt_3d * weight_pres(intermediate_timestep_count) … … 3263 3195 ENDIF 3264 3196 3265 3266 3197 dh = dzw(k) 3267 3198 inv_dh = 1.0_wp / dh 3268 3199 dt_dh = dt_chem/dh 3269 3270 ! 3271 !-- Concentration at i,j,k 3200 ! 3201 !-- Concentration at i,j,k 3272 3202 DO lsp = 1, nspec 3273 3203 cc(lsp) = chem_species(lsp)%conc(k,j,i) 3274 3204 ENDDO 3275 3205 3276 3277 ! 3278 !-- Temperature at i,j,k 3206 !-- Temperature at i,j,k 3279 3207 ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp 3280 3208 3281 3209 ts = ttemp - 273.15 !< in degrees celcius 3282 3283 ! 3284 !-- Viscosity of air 3210 ! 3211 !-- Viscosity of air 3285 3212 visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0) 3286 3287 ! 3288 !-- Air density at k 3213 ! 3214 !-- Air density at k 3289 3215 dens = rho_air_zw(k) 3290 3291 ! 3292 !-- Calculate relative humidity from specific humidity for DEPAC 3216 ! 3217 !-- Calculate relative humidity from specific humidity for DEPAC 3293 3218 qv_tmp = MAX(q(k,j,i),0.0_wp) 3294 3219 relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(k) ) 3295 3296 3297 3298 ! 3299 !-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget 3300 !-- for each surface fraction. Then derive overall budget taking into account the surface fractions. 3301 3302 ! 3303 !-- Vegetation 3220 ! 3221 !-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget 3222 !-- for each surface fraction. Then derive overall budget taking into account the surface fractions. 3223 ! 3224 !-- Vegetation 3304 3225 IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 ) THEN 3305 3226 3306 3307 3227 slinnfac = 1.0_wp 3308 3309 ! 3310 !-- Get deposition velocity vd 3228 ! 3229 !-- Get deposition velocity vd 3311 3230 DO lsp = 1, nvar 3312 3313 !--Initialize3231 ! 3232 !-- Initialize 3314 3233 vs = 0.0_wp 3315 3234 vd_lu = 0.0_wp … … 3319 3238 IF ( spc_names(lsp) == 'PM10' ) THEN 3320 3239 part_type = 1 3321 3322 !--Sedimentation velocity3240 ! 3241 !-- Sedimentation velocity 3323 3242 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3324 3243 particle_pars(ind_p_size, part_type), & … … 3340 3259 ELSEIF ( spc_names(lsp) == 'PM25' ) THEN 3341 3260 part_type = 2 3342 3343 !--Sedimentation velocity3261 ! 3262 !-- Sedimentation velocity 3344 3263 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3345 3264 particle_pars(ind_p_size, part_type), & 3346 3265 particle_pars(ind_p_slip, part_type), & 3347 3266 visc) 3348 3349 3267 3350 3268 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & … … 3359 3277 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3360 3278 3361 3362 3279 ELSE !< GASES 3363 3364 ! 3365 !-- Read spc_name of current species for gas parameter 3366 3280 ! 3281 !-- Read spc_name of current species for gas parameter 3367 3282 IF ( ANY( pspecnames(:) == spc_names(lsp) ) ) THEN 3368 3283 i_pspec = 0 … … 3374 3289 3375 3290 ELSE 3376 3377 !--For now species not deposited3291 ! 3292 !-- For now species not deposited 3378 3293 CYCLE 3379 3294 ENDIF 3380 3381 ! 3382 !-- Factor used for conversion from ppb to ug/m3 : 3383 !-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 3384 ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 3385 ! c 1e-9 xm_tracer 1e9 / xm_air dens 3386 !-- thus: 3387 ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 3388 !-- Use density at k: 3295 ! 3296 !-- Factor used for conversion from ppb to ug/m3 : 3297 !-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 3298 !-- (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 3299 !-- c 1e-9 xm_tracer 1e9 / xm_air dens 3300 !-- thus: 3301 !-- c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 3302 !-- Use density at k: 3389 3303 3390 3304 ppm_to_ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m3 3391 3392 ! 3393 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 3305 ! 3306 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 3394 3307 ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 3395 3308 catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 3396 3397 ! 3398 !-- Diffusivity for DEPAC relevant gases 3399 !-- Use default value 3309 ! 3310 !-- Diffusivity for DEPAC relevant gases 3311 !-- Use default value 3400 3312 diffc = 0.11e-4 3401 3402 !--overwrite with known coefficients of diffusivity from Massman (1998)3313 ! 3314 !-- overwrite with known coefficients of diffusivity from Massman (1998) 3403 3315 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 3404 3316 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 … … 3408 3320 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 3409 3321 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 3410 3411 3412 ! 3413 !-- Get quasi-laminar boundary layer resistance Rb: 3322 ! 3323 !-- Get quasi-laminar boundary layer resistance Rb: 3414 3324 CALL get_rb_cell( (lu_dep == ilu_water_sea) .OR. (lu_dep == ilu_water_inland), & 3415 3325 z0h_lu, ustar_lu, diffc, & 3416 3326 Rb ) 3417 3418 ! 3419 !-- Get Rc_tot 3327 ! 3328 !-- Get Rc_tot 3420 3329 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, & 3421 3330 relh, lai, sai, nwet, lu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 3422 3331 r_aero_lu , Rb ) 3423 3424 3425 ! 3426 !-- Calculate budget 3332 ! 3333 !-- Calculate budget 3427 3334 IF ( Rc_tot <= 0.0 ) THEN 3428 3335 … … 3440 3347 ENDDO 3441 3348 ENDIF 3442 3443 3444 ! 3445 !-- Pavement 3349 ! 3350 !-- Pavement 3446 3351 IF ( surf_lsm_h%frac(ind_pav_green,m) > 0 ) THEN 3447 3448 3449 ! 3450 !-- No vegetation on pavements: 3352 ! 3353 !-- No vegetation on pavements: 3451 3354 lai = 0.0_wp 3452 3355 sai = 0.0_wp 3453 3356 3454 3357 slinnfac = 1.0_wp 3455 3456 ! 3457 !-- Get vd 3358 ! 3359 !-- Get vd 3458 3360 DO lsp = 1, nvar 3459 3460 !--Initialize3361 ! 3362 !-- Initialize 3461 3363 vs = 0.0_wp 3462 3364 vd_lu = 0.0_wp … … 3466 3368 IF ( spc_names(lsp) == 'PM10' ) THEN 3467 3369 part_type = 1 3468 3469 !--Sedimentation velocity:3370 ! 3371 !-- Sedimentation velocity: 3470 3372 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3471 3373 particle_pars(ind_p_size, part_type), & … … 3487 3389 ELSEIF ( spc_names(lsp) == 'PM25' ) THEN 3488 3390 part_type = 2 3489 3490 !--Sedimentation velocity:3391 ! 3392 !-- Sedimentation velocity: 3491 3393 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3492 3394 particle_pars(ind_p_size, part_type), & 3493 3395 particle_pars(ind_p_slip, part_type), & 3494 3396 visc) 3495 3496 3397 3497 3398 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & … … 3506 3407 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3507 3408 3508 3509 3409 ELSE !<GASES 3510 3511 ! 3512 !-- Read spc_name of current species for gas parameter 3410 ! 3411 !-- Read spc_name of current species for gas parameter 3513 3412 3514 3413 IF ( ANY(pspecnames(:) == spc_names(lsp) ) ) THEN … … 3521 3420 3522 3421 ELSE 3523 3524 !--For now species not deposited3422 ! 3423 !-- For now species not deposited 3525 3424 CYCLE 3526 3425 ENDIF 3527 3528 !--Factor used for conversion from ppb to ug/m3 :3529 !ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &3530 !(mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)3531 !c 1e-9 xm_tracer 1e9 / xm_air dens3532 !--thus:3533 !c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm33534 !--Use density at lowest layer:3426 ! 3427 !-- Factor used for conversion from ppb to ug/m3 : 3428 !-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 3429 !-- (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 3430 !-- c 1e-9 xm_tracer 1e9 / xm_air dens 3431 !-- thus: 3432 !-- c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 3433 !-- Use density at lowest layer: 3535 3434 3536 3435 ppm_to_ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m3 3537 3538 ! 3539 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 3436 ! 3437 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 3540 3438 ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 3541 3439 catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 3542 3543 ! 3544 !-- Diffusivity for DEPAC relevant gases 3545 !-- Use default value 3440 ! 3441 !-- Diffusivity for DEPAC relevant gases 3442 !-- Use default value 3546 3443 diffc = 0.11e-4 3547 3548 !--overwrite with known coefficients of diffusivity from Massman (1998)3444 ! 3445 !-- overwrite with known coefficients of diffusivity from Massman (1998) 3549 3446 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 3550 3447 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 … … 3554 3451 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 3555 3452 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 3556 3557 3558 ! 3559 !-- Get quasi-laminar boundary layer resistance Rb: 3560 CALL get_rb_cell( (lup_dep == ilu_water_sea) .OR. (lup_dep == ilu_water_inland), & 3561 z0h_lu, ustar_lu, diffc, & 3562 Rb ) 3563 3564 3565 ! 3566 !-- Get Rc_tot 3567 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, & 3568 relh, lai, sai, nwet, lup_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 3569 r_aero_lu , Rb ) 3570 3571 3572 ! 3573 !-- Calculate budget 3453 ! 3454 !-- Get quasi-laminar boundary layer resistance Rb: 3455 CALL get_rb_cell( (lup_dep == ilu_water_sea) .OR. (lup_dep == ilu_water_inland), & 3456 z0h_lu, ustar_lu, diffc, Rb ) 3457 ! 3458 !-- Get Rc_tot 3459 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, & 3460 glrad, cos_zenith, relh, lai, sai, nwet, lup_dep, 2, & 3461 rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 3462 r_aero_lu , Rb ) 3463 ! 3464 !-- Calculate budget 3574 3465 IF ( Rc_tot <= 0.0 ) THEN 3575 3576 3466 bud_now_lup(lsp) = 0.0_wp 3577 3578 3467 ELSE 3579 3580 3468 vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot ) 3581 3582 3469 bud_now_lup(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * & 3583 3470 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3584 3471 ENDIF 3585 3472 3586 3587 3473 ENDIF 3588 3474 ENDDO 3589 3475 ENDIF 3590 3591 3592 ! 3593 !-- Water 3476 ! 3477 !-- Water 3594 3478 IF ( surf_lsm_h%frac(ind_wat_win,m) > 0 ) THEN 3595 3596 3597 ! 3598 !-- No vegetation on water: 3479 ! 3480 !-- No vegetation on water: 3599 3481 lai = 0.0_wp 3600 3482 sai = 0.0_wp 3601 3602 3483 slinnfac = 1.0_wp 3603 3604 ! 3605 !-- Get vd 3484 ! 3485 !-- Get vd 3606 3486 DO lsp = 1, nvar 3607 3608 !--Initialize3487 ! 3488 !-- Initialize 3609 3489 vs = 0.0_wp 3610 3490 vd_lu = 0.0_wp … … 3614 3494 IF ( spc_names(lsp) == 'PM10' ) THEN 3615 3495 part_type = 1 3616 3617 !--Sedimentation velocity:3496 ! 3497 !-- Sedimentation velocity: 3618 3498 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3619 3499 particle_pars(ind_p_size, part_type), & … … 3632 3512 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3633 3513 3634 3635 3514 ELSEIF ( spc_names(lsp) == 'PM25' ) THEN 3636 3515 part_type = 2 3637 3638 !--Sedimentation velocity:3516 ! 3517 !-- Sedimentation velocity: 3639 3518 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 3640 3519 particle_pars(ind_p_size, part_type), & 3641 3520 particle_pars(ind_p_slip, part_type), & 3642 3521 visc) 3643 3644 3522 3645 3523 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & … … 3654 3532 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 3655 3533 3656 3657 3534 ELSE !<GASES 3658 3659 ! 3660 !-- Read spc_name of current species for gas PARAMETER 3535 ! 3536 !-- Read spc_name of current species for gas PARAMETER 3661 3537 3662 3538 IF ( ANY(pspecnames(:) == spc_names(lsp) ) ) THEN … … 3669 3545 3670 3546 ELSE 3671 3672 !--For now species not deposited3547 ! 3548 !-- For now species not deposited 3673 3549 CYCLE 3674 3550 ENDIF 3675 3676 !-- Factor used for conversion from ppb to ug/m3 : 3677 ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 3678 ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 3679 ! c 1e-9 xm_tracer 1e9 / xm_air dens 3680 !-- thus: 3681 ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 3682 !-- Use density at lowest layer: 3683 3684 ppm_to_ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m3 3685 3686 ! 3687 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 3688 ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 3689 catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 3690 3691 ! 3692 !-- Diffusivity for DEPAC relevant gases 3693 !-- Use default value 3551 ! 3552 !-- Factor used for conversion from ppb to ug/m3 : 3553 !-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 3554 !-- (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 3555 !-- c 1e-9 xm_tracer 1e9 / xm_air dens 3556 !-- thus: 3557 !-- c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 3558 !-- Use density at lowest layer: 3559 3560 ppm_to_ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m3 3561 ! 3562 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 3563 !-- ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 3564 catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 3565 ! 3566 !-- Diffusivity for DEPAC relevant gases 3567 !-- Use default value 3694 3568 diffc = 0.11e-4 3695 3696 !--overwrite with known coefficients of diffusivity from Massman (1998)3569 ! 3570 !-- overwrite with known coefficients of diffusivity from Massman (1998) 3697 3571 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 3698 3572 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 … … 3702 3576 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 3703 3577 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 3704 3705 3706 ! 3707 !-- Get quasi-laminar boundary layer resistance Rb: 3708 CALL get_rb_cell( (luw_dep == ilu_water_sea) .OR. (luw_dep == ilu_water_inland), & 3709 z0h_lu, ustar_lu, diffc, & 3710 Rb ) 3711 3712 ! 3713 !-- Get Rc_tot 3714 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, & 3715 relh, lai, sai, nwet, luw_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 3716 r_aero_lu , Rb ) 3717 3718 3719 ! Calculate budget 3578 ! 3579 !-- Get quasi-laminar boundary layer resistance Rb: 3580 CALL get_rb_cell( (luw_dep == ilu_water_sea) .OR. (luw_dep == ilu_water_inland), & 3581 z0h_lu, ustar_lu, diffc, rb ) 3582 3583 !-- Get Rc_tot 3584 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, & 3585 glrad, cos_zenith, relh, lai, sai, nwet, luw_dep, 2, & 3586 rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 3587 r_aero_lu , Rb ) 3588 ! 3589 !-- Calculate budget 3720 3590 IF ( Rc_tot <= 0.0 ) THEN 3721 3591 … … 3736 3606 3737 3607 bud_now = 0.0_wp 3738 3739 !--Calculate overall budget for surface m and adapt concentration3608 ! 3609 !-- Calculate overall budget for surface m and adapt concentration 3740 3610 DO lsp = 1, nspec 3741 3742 3611 3743 3612 bud_now(lsp) = surf_lsm_h%frac(ind_veg_wall,m) * bud_now_lu(lsp) + & 3744 3613 surf_lsm_h%frac(ind_pav_green,m) * bud_now_lup(lsp) + & 3745 3614 surf_lsm_h%frac(ind_wat_win,m) * bud_now_luw(lsp) 3746 3747 ! 3748 !-- Compute new concentration: 3615 ! 3616 !-- Compute new concentration: 3749 3617 cc(lsp) = cc(lsp) + bud_now(lsp) * inv_dh 3750 3618 … … 3754 3622 3755 3623 ENDIF 3756 3757 3758 3759 3760 ! 3761 !-- For USM surfaces 3624 ! 3625 !-- For USM surfaces 3762 3626 3763 3627 IF ( match_usm ) THEN 3764 3765 ! 3766 !-- Get surface element information at i,j: 3628 ! 3629 !-- Get surface element information at i,j: 3767 3630 m = surf_usm_h%start_index(j,i) 3768 3631 k = surf_usm_h%k(m) 3769 3770 ! 3771 !-- Get needed variables for surface element m 3632 ! 3633 !-- Get needed variables for surface element m 3772 3634 ustar_lu = surf_usm_h%us(m) 3773 3635 z0h_lu = surf_usm_h%z0h(m) … … 3776 3638 lai = surf_usm_h%lai(m) 3777 3639 sai = lai + 1 3778 3779 ! 3780 !-- For small grid spacing neglect R_a 3640 ! 3641 !-- For small grid spacing neglect R_a 3781 3642 IF ( dzw(k) <= 1.0 ) THEN 3782 3643 r_aero_lu = 0.0_wp 3783 3644 ENDIF 3784 3785 ! 3786 !-- Initialize lu's 3645 ! 3646 !-- Initialize lu's 3787 3647 luu_palm = 0 3788 3648 luu_dep = 0 … … 3791 3651 lud_palm = 0 3792 3652 lud_dep = 0 3793 3794 ! 3795 !-- Initialize budgets 3653 ! 3654 !-- Initialize budgets 3796 3655 bud_now_luu = 0.0_wp 3797 3656 bud_now_lug = 0.0_wp 3798 3657 bud_now_lud = 0.0_wp 3799 3800 3801 ! 3802 !-- Get land use for i,j and assign to DEPAC lu 3658 ! 3659 !-- Get land use for i,j and assign to DEPAC lu 3803 3660 IF ( surf_usm_h%frac(ind_pav_green,m) > 0 ) THEN 3804 3805 !--For green urban surfaces (e.g. green roofs3806 !--assume LU short grass3661 ! 3662 !-- For green urban surfaces (e.g. green roofs 3663 !-- assume LU short grass 3807 3664 lug_palm = ind_lu_s_grass 3808 3665 IF ( lug_palm == ind_lu_user ) THEN … … 3849 3706 3850 3707 IF ( surf_usm_h%frac(ind_veg_wall,m) > 0 ) THEN 3851 3852 !--For walls in USM assume concrete walls/roofs,3853 !--assumed LU class desert as also assumed for3854 !--pavements in LSM3708 ! 3709 !-- For walls in USM assume concrete walls/roofs, 3710 !-- assumed LU class desert as also assumed for 3711 !-- pavements in LSM 3855 3712 luu_palm = ind_lup_conc 3856 3713 IF ( luu_palm == ind_lup_user ) THEN … … 3891 3748 3892 3749 IF ( surf_usm_h%frac(ind_wat_win,m) > 0 ) THEN 3893 3894 !--For windows in USM assume metal as this is3895 !--as close as we get, assumed LU class desert3896 !--as also assumed for pavements in LSM3750 ! 3751 !-- For windows in USM assume metal as this is 3752 !-- as close as we get, assumed LU class desert 3753 !-- as also assumed for pavements in LSM 3897 3754 lud_palm = ind_lup_metal 3898 3755 IF ( lud_palm == ind_lup_user ) THEN … … 3931 3788 ENDIF 3932 3789 ENDIF 3933 3934 3935 ! 3936 !-- @TODO: Activate these lines as soon as new ebsolver branch is merged: 3937 !-- Set wetness indicator to dry or wet for usm vegetation or pavement 3790 ! 3791 !-- @TODO: Activate these lines as soon as new ebsolver branch is merged: 3792 !-- Set wetness indicator to dry or wet for usm vegetation or pavement 3938 3793 !IF ( surf_usm_h%c_liq(m) > 0 ) THEN 3939 3794 ! nwet = 1 … … 3941 3796 nwet = 0 3942 3797 !ENDIF 3943 3944 ! 3945 !-- Compute length of time step 3798 ! 3799 !-- Compute length of time step 3946 3800 IF ( call_chem_at_all_substeps ) THEN 3947 3801 dt_chem = dt_3d * weight_pres(intermediate_timestep_count) … … 3950 3804 ENDIF 3951 3805 3952 3953 3806 dh = dzw(k) 3954 3807 inv_dh = 1.0_wp / dh 3955 3808 dt_dh = dt_chem/dh 3956 3957 ! 3958 !-- Concentration at i,j,k 3809 ! 3810 !-- Concentration at i,j,k 3959 3811 DO lsp = 1, nspec 3960 3812 cc(lsp) = chem_species(lsp)%conc(k,j,i) 3961 3813 ENDDO 3962 3963 ! 3964 !-- Temperature at i,j,k 3814 ! 3815 !-- Temperature at i,j,k 3965 3816 ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp 3966 3817 3967 3818 ts = ttemp - 273.15 !< in degrees celcius 3968 3969 ! 3970 !-- Viscosity of air 3819 ! 3820 !-- Viscosity of air 3971 3821 visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0) 3972 3973 ! 3974 !-- Air density at k 3822 ! 3823 !-- Air density at k 3975 3824 dens = rho_air_zw(k) 3976 3977 ! 3978 !-- Calculate relative humidity from specific humidity for DEPAC 3825 ! 3826 !-- Calculate relative humidity from specific humidity for DEPAC 3979 3827 qv_tmp = MAX(q(k,j,i),0.0_wp) 3980 3828 relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(nzb) ) 3981 3982 3983 3984 ! 3985 !-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget 3986 !-- for each surface fraction. Then derive overall budget taking into account the surface fractions. 3987 3988 ! 3989 !-- Walls/roofs 3829 ! 3830 !-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget 3831 !-- for each surface fraction. Then derive overall budget taking into account the surface fractions. 3832 ! 3833 !-- Walls/roofs 3990 3834 IF ( surf_usm_h%frac(ind_veg_wall,m) > 0 ) THEN 3991 3992 3993 ! 3994 !-- No vegetation on non-green walls: 3835 ! 3836 !-- No vegetation on non-green walls: 3995 3837 lai = 0.0_wp 3996 3838 sai = 0.0_wp 3997 3839 3998 3840 slinnfac = 1.0_wp 3999 4000 ! 4001 !-- Get vd 3841 ! 3842 !-- Get vd 4002 3843 DO lsp = 1, nvar 4003 4004 !--Initialize3844 ! 3845 !-- Initialize 4005 3846 vs = 0.0_wp 4006 3847 vd_lu = 0.0_wp … … 4010 3851 IF (spc_names(lsp) == 'PM10' ) THEN 4011 3852 part_type = 1 4012 4013 !--Sedimentation velocity3853 ! 3854 !-- Sedimentation velocity 4014 3855 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 4015 3856 particle_pars(ind_p_size, part_type), & … … 4028 3869 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 4029 3870 4030 4031 3871 ELSEIF ( spc_names(lsp) == 'PM25' ) THEN 4032 3872 part_type = 2 4033 4034 !--Sedimentation velocity3873 ! 3874 !-- Sedimentation velocity 4035 3875 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 4036 3876 particle_pars(ind_p_size, part_type), & 4037 3877 particle_pars(ind_p_slip, part_type), & 4038 3878 visc) 4039 4040 3879 4041 3880 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & … … 4050 3889 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 4051 3890 4052 4053 3891 ELSE !< GASES 4054 4055 ! 4056 !-- Read spc_name of current species for gas parameter 3892 ! 3893 !-- Read spc_name of current species for gas parameter 4057 3894 4058 3895 IF ( ANY( pspecnames(:) == spc_names(lsp) ) ) THEN … … 4064 3901 ENDDO 4065 3902 ELSE 4066 4067 !--For now species not deposited3903 ! 3904 !-- For now species not deposited 4068 3905 CYCLE 4069 3906 ENDIF 3907 ! 3908 !-- Factor used for conversion from ppb to ug/m3 : 3909 !-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 3910 !-- (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 3911 !-- c 1e-9 xm_tracer 1e9 / xm_air dens 3912 !-- thus: 3913 !-- c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 3914 !-- Use density at k: 3915 3916 ppm_to_ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m3 4070 3917 4071 3918 ! 4072 !-- Factor used for conversion from ppb to ug/m3 : 4073 ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 4074 ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 4075 ! c 1e-9 xm_tracer 1e9 / xm_air dens 4076 !-- thus: 4077 ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 4078 !-- Use density at k: 4079 4080 ppm_to_ugm3 = (dens/xm_air) * 0.001_wp !< (mole air)/m3 4081 4082 ! 4083 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 4084 ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 4085 catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 4086 4087 ! 4088 !-- Diffusivity for DEPAC relevant gases 4089 !-- Use default value 3919 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 3920 !-- ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 3921 !-- catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 3922 ! 3923 !-- Diffusivity for DEPAC relevant gases 3924 !-- Use default value 4090 3925 diffc = 0.11e-4 4091 4092 !--overwrite with known coefficients of diffusivity from Massman (1998)3926 ! 3927 !-- overwrite with known coefficients of diffusivity from Massman (1998) 4093 3928 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 4094 3929 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 … … 4098 3933 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 4099 3934 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 4100 4101 4102 ! 4103 !-- Get quasi-laminar boundary layer resistance Rb: 4104 CALL get_rb_cell( (luu_dep == ilu_water_sea) .OR. (luu_dep == ilu_water_inland), & 3935 ! 3936 !-- Get quasi-laminar boundary layer resistance Rb: 3937 CALL get_rb_cell( (luu_dep == ilu_water_sea) .OR. (luu_dep == ilu_water_inland), & 4105 3938 z0h_lu, ustar_lu, diffc, & 4106 3939 Rb ) 4107 4108 ! 4109 !-- Get Rc_tot 4110 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, & 4111 relh, lai, sai, nwet, luu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 4112 r_aero_lu , Rb ) 4113 4114 4115 ! 4116 !-- Calculate budget 3940 ! 3941 !-- Get Rc_tot 3942 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, & 3943 glrad, cos_zenith, relh, lai, sai, nwet, luu_dep, 2, & 3944 rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 3945 r_aero_lu, rb ) 3946 ! 3947 !-- Calculate budget 4117 3948 IF ( Rc_tot <= 0.0 ) THEN 4118 3949 … … 4130 3961 ENDDO 4131 3962 ENDIF 4132 4133 4134 ! 4135 !-- Green usm surfaces 3963 ! 3964 !-- Green usm surfaces 4136 3965 IF ( surf_usm_h%frac(ind_pav_green,m) > 0 ) THEN 4137 3966 4138 4139 3967 slinnfac = 1.0_wp 4140 4141 ! 4142 !-- Get vd 3968 ! 3969 !-- Get vd 4143 3970 DO lsp = 1, nvar 4144 4145 !--Initialize3971 ! 3972 !-- Initialize 4146 3973 vs = 0.0_wp 4147 3974 vd_lu = 0.0_wp … … 4151 3978 IF ( spc_names(lsp) == 'PM10' ) THEN 4152 3979 part_type = 1 4153 ! Sedimentation velocity 3980 ! 3981 !-- Sedimentation velocity 4154 3982 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 4155 3983 particle_pars(ind_p_size, part_type), & … … 4168 3996 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 4169 3997 4170 4171 3998 ELSEIF ( spc_names(lsp) == 'PM25' ) THEN 4172 3999 part_type = 2 4173 ! Sedimentation velocity 4000 ! 4001 !-- Sedimentation velocity 4174 4002 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 4175 4003 particle_pars(ind_p_size, part_type), & 4176 4004 particle_pars(ind_p_slip, part_type), & 4177 4005 visc) 4178 4179 4006 4180 4007 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & … … 4189 4016 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 4190 4017 4191 4192 4018 ELSE !< GASES 4193 4194 ! 4195 !-- Read spc_name of current species for gas parameter 4019 ! 4020 !-- Read spc_name of current species for gas parameter 4196 4021 4197 4022 IF ( ANY( pspecnames(:) == spc_names(lsp) ) ) THEN … … 4203 4028 ENDDO 4204 4029 ELSE 4205 4206 !--For now species not deposited4030 ! 4031 !-- For now species not deposited 4207 4032 CYCLE 4208 4033 ENDIF 4209 4210 ! 4211 !-- Factor used for conversion from ppb to ug/m3 : 4212 ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 4213 ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 4214 ! c 1e-9 xm_tracer 1e9 / xm_air dens 4215 !-- thus: 4216 ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 4217 !-- Use density at k: 4034 ! 4035 !-- Factor used for conversion from ppb to ug/m3 : 4036 !-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 4037 !-- (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 4038 !-- c 1e-9 xm_tracer 1e9 / xm_air dens 4039 !-- thus: 4040 !-- c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 4041 !-- Use density at k: 4218 4042 4219 4043 ppm_to_ugm3 = (dens/xm_air) * 0.001_wp ! (mole air)/m3 4220 4221 ! 4222 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 4044 ! 4045 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 4223 4046 ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 4224 4047 catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 4225 4226 ! 4227 !-- Diffusivity for DEPAC relevant gases 4228 !-- Use default value 4048 ! 4049 !-- Diffusivity for DEPAC relevant gases 4050 !-- Use default value 4229 4051 diffc = 0.11e-4 4230 4231 !--overwrite with known coefficients of diffusivity from Massman (1998)4052 ! 4053 !-- overwrite with known coefficients of diffusivity from Massman (1998) 4232 4054 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 4233 4055 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 … … 4237 4059 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 4238 4060 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 4239 4240 4241 ! 4242 !-- Get quasi-laminar boundary layer resistance Rb: 4243 CALL get_rb_cell( (lug_dep == ilu_water_sea) .OR. (lug_dep == ilu_water_inland), & 4061 ! 4062 !-- Get quasi-laminar boundary layer resistance Rb: 4063 CALL get_rb_cell( (lug_dep == ilu_water_sea) .OR. (lug_dep == ilu_water_inland), & 4244 4064 z0h_lu, ustar_lu, diffc, & 4245 4065 Rb ) 4246 4247 4248 ! 4249 !-- Get Rc_tot 4250 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, & 4251 relh, lai, sai, nwet, lug_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 4252 r_aero_lu , Rb ) 4253 4254 4255 ! 4256 !-- Calculate budget 4066 ! 4067 !-- Get Rc_tot 4068 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, & 4069 glrad, cos_zenith, relh, lai, sai, nwet, lug_dep, 2, & 4070 rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 4071 r_aero_lu , Rb ) 4072 ! 4073 !-- Calculate budget 4257 4074 IF ( Rc_tot <= 0.0 ) THEN 4258 4075 … … 4267 4084 ENDIF 4268 4085 4269 4270 4086 ENDIF 4271 4087 ENDDO 4272 4088 ENDIF 4273 4274 4275 ! 4276 !-- Windows 4089 ! 4090 !-- Windows 4277 4091 IF ( surf_usm_h%frac(ind_wat_win,m) > 0 ) THEN 4278 4279 4280 ! 4281 !-- No vegetation on windows: 4092 ! 4093 !-- No vegetation on windows: 4282 4094 lai = 0.0_wp 4283 4095 sai = 0.0_wp 4284 4096 4285 4097 slinnfac = 1.0_wp 4286 4287 ! 4288 !-- Get vd 4098 ! 4099 !-- Get vd 4289 4100 DO lsp = 1, nvar 4290 4291 !--Initialize4101 ! 4102 !-- Initialize 4292 4103 vs = 0.0_wp 4293 4104 vd_lu = 0.0_wp … … 4297 4108 IF ( spc_names(lsp) == 'PM10' ) THEN 4298 4109 part_type = 1 4299 4300 !--Sedimentation velocity4110 ! 4111 !-- Sedimentation velocity 4301 4112 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 4302 4113 particle_pars(ind_p_size, part_type), & … … 4304 4115 visc) 4305 4116 4306 CALL drydepo_aero_zhang_vd( vd_lu, Rs, & 4307 vs, & 4117 CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, & 4118 particle_pars(ind_p_size, part_type), & 4119 particle_pars(ind_p_slip, part_type), & 4120 nwet, ttemp, dens, visc, & 4121 lud_dep, r_aero_lu, ustar_lu ) 4122 4123 bud_now_lud(lsp) = - cc(lsp) * & 4124 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 4125 4126 ELSEIF ( spc_names(lsp) == 'PM25' ) THEN 4127 part_type = 2 4128 ! 4129 !-- Sedimentation velocity 4130 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), & 4131 particle_pars(ind_p_size, part_type), & 4132 particle_pars(ind_p_slip, part_type), & 4133 visc) 4134 4135 CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, & 4308 4136 particle_pars(ind_p_size, part_type), & 4309 4137 particle_pars(ind_p_slip, part_type), & … … 4315 4143 (1.0_wp - exp(-vd_lu * dt_dh )) * dh 4316 4144 4317 4318 ELSEIF ( spc_names(lsp) == 'PM25' ) THEN4319 part_type = 24320 !4321 !-- Sedimentation velocity4322 vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &4323 particle_pars(ind_p_size, part_type), &4324 particle_pars(ind_p_slip, part_type), &4325 visc)4326 4327 4328 CALL drydepo_aero_zhang_vd( vd_lu, Rs, &4329 vs, &4330 particle_pars(ind_p_size, part_type), &4331 particle_pars(ind_p_slip, part_type), &4332 nwet, ttemp, dens, visc, &4333 lud_dep, &4334 r_aero_lu, ustar_lu )4335 4336 bud_now_lud(lsp) = - cc(lsp) * &4337 (1.0_wp - exp(-vd_lu * dt_dh )) * dh4338 4339 4340 4145 ELSE !< GASES 4341 4342 ! 4343 !-- Read spc_name of current species for gas PARAMETER 4146 ! 4147 !-- Read spc_name of current species for gas PARAMETER 4344 4148 4345 4149 IF ( ANY( pspecnames(:) == spc_names(lsp) ) ) THEN … … 4351 4155 ENDDO 4352 4156 ELSE 4353 ! For now species not deposited 4157 ! 4158 !-- For now species not deposited 4354 4159 CYCLE 4355 4160 ENDIF 4356 4357 ! 4358 !-- Factor used for conversion from ppb to ug/m3 : 4359 ! ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 4360 ! (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 4361 ! c 1e-9 xm_tracer 1e9 / xm_air dens 4362 !-- thus: 4363 ! c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 4364 !-- Use density at k: 4161 ! 4162 !-- Factor used for conversion from ppb to ug/m3 : 4163 !-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) & 4164 !-- (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole) 4165 !-- c 1e-9 xm_tracer 1e9 / xm_air dens 4166 !-- thus: 4167 !-- c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3 4168 !-- Use density at k: 4365 4169 4366 4170 ppm_to_ugm3 = (dens/xm_air) * 0.001_wp ! (mole air)/m3 4367 4368 ! 4369 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 4370 ! ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 4171 ! 4172 !-- Atmospheric concentration in DEPAC is requested in ug/m3: 4173 !-- ug/m3 ppm (ug/m3)/ppm/(kg/mole) kg/mole 4371 4174 catm = cc(lsp) * ppm_to_ugm3 * specmolm(i_pspec) ! in ug/m3 4372 4373 ! 4374 !-- Diffusivity for DEPAC relevant gases 4375 !-- Use default value 4376 diffc = 0.11e-4 4377 ! 4378 !-- overwrite with known coefficients of diffusivity from Massman (1998) 4175 ! 4176 !-- Diffusivity for DEPAC relevant gases 4177 !-- Use default value 4178 diffc = 0.11e-4 4179 ! 4180 !-- overwrite with known coefficients of diffusivity from Massman (1998) 4379 4181 IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 4380 4182 IF ( spc_names(lsp) == 'NO' ) diffc = 0.199e-4 … … 4384 4186 IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4 4385 4187 IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4 4386 4387 4388 ! 4389 !-- Get quasi-laminar boundary layer resistance Rb: 4390 CALL get_rb_cell( (lud_dep == ilu_water_sea) .OR. (lud_dep == ilu_water_inland), & 4391 z0h_lu, ustar_lu, diffc, & 4392 Rb ) 4393 4394 ! 4395 !-- Get Rc_tot 4396 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, & 4397 relh, lai, sai, nwet, lud_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 4398 r_aero_lu , Rb ) 4399 4400 4401 ! 4402 !-- Calculate budget 4188 ! 4189 !-- Get quasi-laminar boundary layer resistance Rb: 4190 CALL get_rb_cell( (lud_dep == ilu_water_sea) .OR. (lud_dep == ilu_water_inland), & 4191 z0h_lu, ustar_lu, diffc, rb ) 4192 ! 4193 !-- Get Rc_tot 4194 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, & 4195 glrad, cos_zenith, relh, lai, sai, nwet, lud_dep, 2, & 4196 rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 4197 r_aero_lu , rb ) 4198 ! 4199 !-- Calculate budget 4403 4200 IF ( Rc_tot <= 0.0 ) THEN 4404 4201 … … 4419 4216 4420 4217 bud_now = 0.0_wp 4421 4422 !--Calculate overall budget for surface m and adapt concentration4218 ! 4219 !-- Calculate overall budget for surface m and adapt concentration 4423 4220 DO lsp = 1, nspec 4424 4221 … … 4427 4224 surf_usm_h%frac(ind_pav_green,m) * bud_now_lug(lsp) + & 4428 4225 surf_usm_h%frac(ind_wat_win,m) * bud_now_lud(lsp) 4429 4430 ! 4431 !-- Compute new concentration 4226 ! 4227 !-- Compute new concentration 4432 4228 cc(lsp) = cc(lsp) + bud_now(lsp) * inv_dh 4433 4229 … … 4439 4235 4440 4236 4441 4442 4237 END SUBROUTINE chem_depo 4443 4238 4444 4239 4445 4446 !---------------------------------------------------------------------------------- 4447 ! 4448 !> DEPAC: 4449 !> Code of the DEPAC routine and corresponding subroutines below from the DEPAC 4450 !> module of the LOTOS-EUROS model (Manders et al., 2017) 4451 ! 4452 !> Original DEPAC routines by RIVM and TNO (2015), for Documentation see 4453 !> van Zanten et al., 2010. 4454 !--------------------------------------------------------------------------------- 4455 ! 4456 !---------------------------------------------------------------------------------- 4457 ! 4458 !> depac: compute total canopy (or surface) resistance Rc for gases 4459 !---------------------------------------------------------------------------------- 4460 4240 !------------------------------------------------------------------------------! 4241 ! Description: 4242 ! ------------ 4243 !> Subroutine to compute total canopy (or surface) resistance Rc for gases 4244 !> 4245 !> DEPAC: 4246 !> Code of the DEPAC routine and corresponding subroutines below from the DEPAC 4247 !> module of the LOTOS-EUROS model (Manders et al., 2017) 4248 !> 4249 !> Original DEPAC routines by RIVM and TNO (2015), for Documentation see 4250 !> van Zanten et al., 2010. 4251 !------------------------------------------------------------------------------! 4461 4252 SUBROUTINE drydepos_gas_depac( compnam, day_of_year, lat, t, ust, glrad, sinphi, & 4462 4253 rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, catm, diffc, & 4463 4254 ra, rb ) 4464 4465 4466 ! 4467 !--Some of depac arguments are OPTIONAL: 4468 ! 4469 !-- A. compute Rc_tot without compensation points (ccomp_tot will be zero): 4470 !-- CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi]) 4471 ! 4472 !-- B. compute Rc_tot with compensation points (used for LOTOS-EUROS): 4473 !-- CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & 4474 ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water) 4475 ! 4476 !-- C. compute effective Rc based on compensation points (used for OPS): 4477 !-- CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & 4478 ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, & 4479 ! ra, rb, rc_eff) 4480 !-- X1. Extra (OPTIONAL) output variables: 4481 !-- CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & 4482 ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, & 4483 ! ra, rb, rc_eff, & 4484 ! gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out) 4485 !-- X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves): 4486 !-- CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & 4487 ! c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, & 4488 ! ra, rb, rc_eff, & 4489 ! gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out, & 4490 ! calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA) 4491 4255 ! 4256 !-- Some of depac arguments are OPTIONAL: 4257 !-- A. compute Rc_tot without compensation points (ccomp_tot will be zero): 4258 !-- CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi]) 4259 !-- B. compute Rc_tot with compensation points (used for LOTOS-EUROS): 4260 !-- CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & 4261 !-- c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water) 4262 !-- 4263 !-- C. compute effective Rc based on compensation points (used for OPS): 4264 !-- CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & 4265 !-- c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, & 4266 !-- ra, rb, rc_eff) 4267 !-- X1. Extra (OPTIONAL) output variables: 4268 !-- CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & 4269 !-- c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, & 4270 !-- ra, rb, rc_eff, & 4271 !-- gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out) 4272 !-- X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves): 4273 !-- CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], & 4274 !-- c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, & 4275 !-- ra, rb, rc_eff, & 4276 !-- gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out, & 4277 !-- calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA) 4492 4278 4493 4279 IMPLICIT NONE 4494 4280 4495 CHARACTER( len=*), INTENT(IN) :: compnam !< component name4281 CHARACTER(LEN=*), INTENT(IN) :: compnam !< component name 4496 4282 !< 'HNO3','NO','NO2','O3','SO2','NH3' 4497 4283 INTEGER(iwp), INTENT(IN) :: day_of_year !< day of year, 1 ... 365 (366) … … 4518 4304 REAL(wp), INTENT(OUT) :: rc_tot !< total canopy resistance Rc (s/m) 4519 4305 REAL(wp), INTENT(OUT) :: ccomp_tot !< total compensation point (ug/m3) 4520 !< [= 0 for species that don't have a compensation 4521 ! 4522 !-- Local variables: 4523 ! 4524 !-- Component number taken from component name, paramteres matched with include files 4306 ! !< [= 0 for species that don't have a compensation 4307 !-- Local variables: 4308 ! 4309 !-- Component number taken from component name, paramteres matched with include files 4525 4310 INTEGER(iwp) :: icmp 4526 4527 ! 4528 !-- Component numbers: 4311 ! 4312 !-- Component numbers: 4529 4313 INTEGER(iwp), PARAMETER :: icmp_o3 = 1 4530 4314 INTEGER(iwp), PARAMETER :: icmp_so2 = 2 … … 4541 4325 !< = 1 -> constant Rc 4542 4326 !< = 2 -> temperature dependent Rc 4543 4544 4327 ! 4328 !-- Vegetation indicators: 4545 4329 LOGICAL :: lai_present !< leaves are present for current land use type 4546 4330 LOGICAL :: sai_present !< vegetation is present for current land use type … … 4554 4338 REAL(wp) :: cstom !< stomatal compensation point (ug/m3) 4555 4339 REAL(wp) :: csoil !< soil compensation point (ug/m3) 4556 4557 4558 4340 ! 4559 4341 !-- Next statement is just to avoid compiler warning about unused variable 4560 4342 IF ( day_of_year == 0 .OR. ( catm + lat + ra + rb ) > 0.0_wp ) CONTINUE 4561 4562 4563 ! 4564 !-- Define component number 4343 ! 4344 !-- Define component number 4565 4345 SELECT CASE ( TRIM( compnam ) ) 4566 4346 … … 4596 4376 4597 4377 CASE default 4598 4599 !--Component not part of DEPAC --> not deposited4378 ! 4379 !-- Component not part of DEPAC --> not deposited 4600 4380 RETURN 4601 4381 4602 4382 END SELECT 4603 4383 4604 4605 4384 ! 4385 !-- Inititalize 4606 4386 gw = 0.0_wp 4607 4387 gstom = 0.0_wp … … 4611 4391 cstom = 0.0_wp 4612 4392 csoil = 0.0_wp 4613 4614 4615 ! 4616 !-- Check whether vegetation is present: 4393 ! 4394 !-- Check whether vegetation is present: 4617 4395 lai_present = ( lai > 0.0 ) 4618 4396 sai_present = ( sai > 0.0 ) 4619 4620 ! 4621 !-- Set Rc (i.e. rc_tot) in special cases: 4397 ! 4398 !-- Set Rc (i.e. rc_tot) in special cases: 4622 4399 CALL rc_special( icmp, compnam, lu, t, nwet, rc_tot, ready, ccomp_tot ) 4623 4624 ! 4625 !-- If Rc is not set: 4400 ! 4401 !-- If Rc is not set: 4626 4402 IF ( .NOT. ready ) then 4627 4628 ! 4629 !-- External conductance: 4403 ! 4404 !-- External conductance: 4630 4405 CALL rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai,gw ) 4631 4632 ! 4633 !-- Stomatal conductance: 4406 ! 4407 !-- Stomatal conductance: 4634 4408 CALL rc_gstom( icmp, compnam, lu, lai_present, lai, glrad, sinphi, t, rh, diffc, gstom, p ) 4635 4636 !--Effective soil conductance:4409 ! 4410 !-- Effective soil conductance: 4637 4411 CALL rc_gsoil_eff( icmp, lu, sai, ust, nwet, t, gsoil_eff ) 4638 4639 ! 4640 !-- Total canopy conductance (gc_tot) and resistance Rc (rc_tot): 4412 ! 4413 !-- Total canopy conductance (gc_tot) and resistance Rc (rc_tot): 4641 4414 CALL rc_rctot( gstom, gsoil_eff, gw, gc_tot, rc_tot ) 4642 4643 ! 4644 !-- Needed to include compensation point for NH3 4645 !-- Compensation points (always returns ccomp_tot; currently ccomp_tot != 0 only for NH3): 4646 !-- CALL rc_comp_point( compnam,lu,day_of_year,t,gw,gstom,gsoil_eff,gc_tot,&4647 ! lai_present, sai_present, &4648 ! ccomp_tot, &4649 ! catm=catm,c_ave_prev_nh3=c_ave_prev_nh3, &4650 ! c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, & 4651 ! tsea=tsea,cw=cw,cstom=cstom,csoil=csoil ) 4652 ! 4653 !-- Effective Rc based on compensation points: 4654 ! IF ( present(rc_eff) ) then 4655 ! check on required arguments: 4656 ! IF ( (.not. present(catm)) .OR. (.not. present(ra)) .OR. (.not. present(rb)) ) then 4657 ! stop 'output argument rc_eff requires input arguments catm, ra and rb' 4658 ! END IF 4659 !--Compute rc_eff :4415 ! 4416 !-- Needed to include compensation point for NH3 4417 !-- Compensation points (always returns ccomp_tot; currently ccomp_tot != 0 only for NH3): 4418 !-- CALL rc_comp_point( compnam,lu,day_of_year,t,gw,gstom,gsoil_eff,gc_tot,& 4419 !-- lai_present, sai_present, & 4420 !-- ccomp_tot, & 4421 !-- catm=catm,c_ave_prev_nh3=c_ave_prev_nh3, & 4422 !-- c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, & 4423 !-- tsea=tsea,cw=cw,cstom=cstom,csoil=csoil ) 4424 ! 4425 !-- Effective Rc based on compensation points: 4426 !-- IF ( present(rc_eff) ) then 4427 !-- check on required arguments: 4428 !-- IF ( (.not. present(catm)) .OR. (.not. present(ra)) .OR. (.not. present(rb)) ) then 4429 !-- stop 'output argument rc_eff requires input arguments catm, ra and rb' 4430 !-- END IF 4431 ! 4432 !-- Compute rc_eff : 4660 4433 ! CALL rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff) 4661 4434 ! ENDIF … … 4665 4438 4666 4439 4667 4668 !------------------------------------------------------------------- 4669 !> rc_special: compute total canopy resistance in special CASEs 4670 !------------------------------------------------------------------- 4440 !------------------------------------------------------------------------------! 4441 ! Description: 4442 ! ------------ 4443 !> Subroutine to compute total canopy resistance in special cases 4444 !------------------------------------------------------------------------------! 4671 4445 SUBROUTINE rc_special( icmp, compnam, lu, t, nwet, rc_tot, ready, ccomp_tot ) 4672 4446 4673 4447 IMPLICIT NONE 4674 4448 4675 CHARACTER( len=*), INTENT(IN):: compnam !< component name4676 4677 INTEGER(iwp), INTENT(IN) :: icmp !< component index4678 INTEGER(iwp), INTENT(IN) :: lu !< land use type, lu = 1,...,nlu4679 INTEGER(iwp), INTENT(IN) :: nwet !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow4680 4681 REAL(wp), INTENT(IN) :: t !< temperature (C)4449 CHARACTER(LEN=*), INTENT(IN) :: compnam !< component name 4450 4451 INTEGER(iwp), INTENT(IN) :: icmp !< component index 4452 INTEGER(iwp), INTENT(IN) :: lu !< land use type, lu = 1,...,nlu 4453 INTEGER(iwp), INTENT(IN) :: nwet !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow 4454 4455 REAL(wp), INTENT(IN) :: t !< temperature (C) 4682 4456 4683 4457 REAL(wp), INTENT(OUT) :: rc_tot !< total canopy resistance Rc (s/m) … … 4686 4460 LOGICAL, INTENT(OUT) :: ready !< Rc has been set 4687 4461 !< = 1 -> constant Rc 4688 !< = 2 -> temperature dependent Rc4689 4690 4462 ! 4691 4463 !-- Next line is to avoid compiler warning about unused variable 4692 4464 IF ( icmp == 0 ) CONTINUE 4693 4694 ! 4695 !-- rc_tot is not yet set: 4465 ! 4466 !-- rc_tot is not yet set: 4696 4467 ready = .FALSE. 4697 4698 ! 4699 !-- Default compensation point in special CASEs = 0: 4468 ! 4469 !-- Default compensation point in special CASEs = 0: 4700 4470 ccomp_tot = 0.0_wp 4701 4471 4702 4472 SELECT CASE( TRIM( compnam ) ) 4703 4473 CASE( 'HNO3', 'N2O5', 'NO3', 'H2O2' ) 4704 4705 !--No separate resistances for HNO3; just one total canopy resistance:4474 ! 4475 !-- No separate resistances for HNO3; just one total canopy resistance: 4706 4476 IF ( t < -5.0_wp .AND. nwet == 9 ) THEN 4707 4708 !--T < 5 C and snow:4477 ! 4478 !-- T < 5 C and snow: 4709 4479 rc_tot = 50.0_wp 4710 4480 ELSE 4711 4712 !--all other circumstances:4481 ! 4482 !-- all other circumstances: 4713 4483 rc_tot = 10.0_wp 4714 4484 ENDIF … … 4724 4494 ENDIF 4725 4495 CASE( 'NO2', 'O3', 'SO2', 'NH3' ) 4726 4727 !--snow surface:4496 ! 4497 !-- snow surface: 4728 4498 IF ( nwet == 9 ) THEN 4729 4730 !--To be activated when snow is implemented4499 ! 4500 !-- To be activated when snow is implemented 4731 4501 !CALL rc_snow(ipar_snow(icmp),t,rc_tot) 4732 4502 ready = .TRUE. … … 4740 4510 4741 4511 4742 4743 !------------------------------------------------------------------- 4744 !> rc_gw: compute external conductance 4745 !------------------------------------------------------------------- 4512 !------------------------------------------------------------------------------! 4513 ! Description: 4514 ! ------------ 4515 !> Subroutine to compute external conductance 4516 !------------------------------------------------------------------------------! 4746 4517 SUBROUTINE rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai, gw ) 4747 4518 4748 4519 IMPLICIT NONE 4749 4750 ! 4751 !-- Input/output variables: 4752 CHARACTER(len=*), INTENT(IN) :: compnam !< component name 4520 ! 4521 !-- Input/output variables: 4522 CHARACTER(LEN=*), INTENT(IN) :: compnam !< component name 4753 4523 4754 4524 INTEGER(iwp), INTENT(IN) :: nwet !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow … … 4757 4527 !< iratns = 2: high NH3/SO2 4758 4528 !< iratns = 3: very low NH3/SO2 4759 4760 4529 LOGICAL, INTENT(IN) :: sai_present 4761 4530 … … 4766 4535 REAL(wp), INTENT(OUT) :: gw !< external leaf conductance (m/s) 4767 4536 4768 4769 4537 SELECT CASE( TRIM( compnam ) ) 4770 4538 … … 4783 4551 CASE( 'NH3' ) 4784 4552 CALL rw_nh3_sutton( t, rh, sai_present, gw ) 4785 4786 ! 4787 !-- conversion from leaf resistance to canopy resistance by multiplying with sai: 4553 ! 4554 !-- conversion from leaf resistance to canopy resistance by multiplying with sai: 4788 4555 gw = sai * gw 4789 4556 4790 4557 CASE default 4791 message_string = 'Component '// TRIM( compnam) // ' not supported'4558 message_string = 'Component '// TRIM( compnam ) // ' not supported' 4792 4559 CALL message( 'rc_gw', 'CM0458', 1, 2, 0, 6, 0 ) 4793 4560 END SELECT … … 4796 4563 4797 4564 4798 4799 !------------------------------------------------------------------- 4800 !> rw_so2: compute external leaf conductance for SO2 4801 !------------------------------------------------------------------- 4565 !------------------------------------------------------------------------------! 4566 ! Description: 4567 ! ------------ 4568 !> Subroutine to compute external leaf conductance for SO2 4569 !------------------------------------------------------------------------------! 4802 4570 SUBROUTINE rw_so2( t, nwet, rh, iratns, sai_present, gw ) 4803 4571 4804 4572 IMPLICIT NONE 4805 4806 ! 4807 !-- Input/output variables: 4573 ! 4574 !-- Input/output variables: 4808 4575 INTEGER(iwp), INTENT(IN) :: nwet !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow 4809 4576 INTEGER(iwp), INTENT(IN) :: iratns !< index for NH3/SO2 ratio: … … 4817 4584 4818 4585 REAL(wp), INTENT(OUT) :: gw !< external leaf conductance (m/s) 4819 4820 ! 4821 !-- Local variables: 4586 ! 4587 !-- Local variables: 4822 4588 REAL(wp) :: rw !< external leaf resistance (s/m) 4823 4824 4825 ! 4826 !-- Check if vegetation present: 4589 ! 4590 !-- Check if vegetation present: 4827 4591 IF ( sai_present ) THEN 4828 4592 4829 4593 IF ( nwet == 0 ) THEN 4830 !-------------------------- 4831 ! dry surface 4832 !-------------------------- 4833 ! T > -1 C 4594 ! 4595 !-- ------------------------ 4596 !-- dry surface 4597 !-- ------------------------ 4598 !-- T > -1 C 4834 4599 IF ( t > -1.0_wp ) THEN 4835 4600 IF ( rh < 81.3_wp ) THEN … … 4848 4613 ENDIF 4849 4614 ELSE 4850 !-------------------------- 4851 ! wet surface 4852 !-------------------------- 4615 ! 4616 !-- ------------------------ 4617 !-- wet surface 4618 !-- ------------------------ 4853 4619 rw = 10.0_wp !see Table 5, Erisman et al, 1994 Atm. Environment, 0 is impl. as 10 4854 4620 ENDIF 4855 4856 !very low NH3/SO2 ratio:4621 ! 4622 !-- very low NH3/SO2 ratio: 4857 4623 IF ( iratns == iratns_very_low ) rw = rw + 50.0_wp 4858 4859 ! Conductance: 4624 ! 4625 !-- Conductance: 4626 gw = 1.0_wp / rw 4627 ELSE 4628 ! 4629 !-- no vegetation: 4630 gw = 0.0_wp 4631 ENDIF 4632 4633 END SUBROUTINE rw_so2 4634 4635 4636 !------------------------------------------------------------------------------! 4637 ! Description: 4638 ! ------------ 4639 !> Subroutine to compute external leaf conductance for NH3, 4640 !> following Sutton & Fowler, 1993 4641 !------------------------------------------------------------------------------! 4642 SUBROUTINE rw_nh3_sutton( tsurf, rh,sai_present, gw ) 4643 4644 IMPLICIT NONE 4645 ! 4646 !-- Input/output variables: 4647 LOGICAL, INTENT(IN) :: sai_present 4648 4649 REAL(wp), INTENT(IN) :: tsurf !< surface temperature (C) 4650 REAL(wp), INTENT(IN) :: rh !< relative humidity (%) 4651 4652 REAL(wp), INTENT(OUT) :: gw !< external leaf conductance (m/s) 4653 ! 4654 !-- Local variables: 4655 REAL(wp) :: rw !< external leaf resistance (s/m) 4656 REAL(wp) :: sai_grass_haarweg !< surface area index at experimental site Haarweg 4657 ! 4658 !-- Fix sai_grass at value valid for Haarweg data for which gamma_w parametrization is derived 4659 sai_grass_haarweg = 3.5_wp 4660 ! 4661 !-- Calculation rw: 4662 !-- 100 - rh 4663 !-- rw = 2.0 * exp(----------) 4664 !-- 12 4665 4666 IF ( sai_present ) THEN 4667 ! 4668 !-- External resistance according to Sutton & Fowler, 1993 4669 rw = 2.0_wp * exp( ( 100.0_wp - rh ) / 12.0_wp ) 4670 rw = sai_grass_haarweg * rw 4671 ! 4672 !-- Frozen soil (from Depac v1): 4673 IF ( tsurf < 0.0_wp ) rw = 200.0_wp 4674 ! 4675 !-- Conductance: 4860 4676 gw = 1.0_wp / rw 4861 4677 ELSE … … 4864 4680 ENDIF 4865 4681 4866 END SUBROUTINE rw_ so24867 4868 4869 4870 !------------------------------------------------------------------- 4871 !> rw_nh3_sutton: compute external leaf conductance for NH3, 4872 !> following Sutton & Fowler, 1993 4873 !------------------------------------------------------------------- 4874 SUBROUTINE rw_ nh3_sutton( tsurf, rh,sai_present, gw )4682 END SUBROUTINE rw_nh3_sutton 4683 4684 4685 !------------------------------------------------------------------------------! 4686 ! Description: 4687 ! ------------ 4688 !> Subroutine to compute external leaf conductance 4689 !------------------------------------------------------------------------------! 4690 SUBROUTINE rw_constant( rw_val, sai_present, gw ) 4875 4691 4876 4692 IMPLICIT NONE 4877 4878 ! 4879 !-- Input/output variables: 4693 ! 4694 !-- Input/output variables: 4880 4695 LOGICAL, INTENT(IN) :: sai_present 4881 4696 4882 REAL(wp), INTENT(IN) :: tsurf !< surface temperature (C)4883 REAL(wp), INTENT(IN) :: rh !< relative humidity (%)4884 4885 REAL(wp), INTENT(OUT) :: gw !< external leaf conductance (m/s)4886 4887 !4888 !-- Local variables:4889 REAL(wp) :: rw !< external leaf resistance (s/m)4890 REAL(wp) :: sai_grass_haarweg !< surface area index at experimental site Haarweg4891 4892 !4893 !-- Fix sai_grass at value valid for Haarweg data for which gamma_w parametrization is derived4894 sai_grass_haarweg = 3.5_wp4895 4896 !4897 !-- Calculation rw:4898 ! 100 - rh4899 ! rw = 2.0 * exp(----------)4900 ! 124901 4902 IF ( sai_present ) THEN4903 4904 ! External resistance according to Sutton & Fowler, 19934905 rw = 2.0_wp * exp( ( 100.0_wp - rh ) / 12.0_wp )4906 rw = sai_grass_haarweg * rw4907 4908 ! Frozen soil (from Depac v1):4909 IF ( tsurf < 0.0_wp ) rw = 200.0_wp4910 4911 ! Conductance:4912 gw = 1.0_wp / rw4913 ELSE4914 ! no vegetation:4915 gw = 0.0_wp4916 ENDIF4917 4918 END SUBROUTINE rw_nh3_sutton4919 4920 4921 4922 !-------------------------------------------------------------------4923 !> rw_constant: compute constant external leaf conductance4924 !-------------------------------------------------------------------4925 SUBROUTINE rw_constant( rw_val, sai_present, gw )4926 4927 IMPLICIT NONE4928 4929 !4930 !-- Input/output variables:4931 LOGICAL, INTENT(IN) :: sai_present4932 4933 4697 REAL(wp), INTENT(IN) :: rw_val !< constant value of Rw 4934 4698 4935 4699 REAL(wp), INTENT(OUT) :: gw !< wernal leaf conductance (m/s) 4936 4937 ! 4938 !-- Compute conductance: 4700 ! 4701 !-- Compute conductance: 4939 4702 IF ( sai_present .AND. .NOT.missing(rw_val) ) THEN 4940 4703 gw = 1.0_wp / rw_val … … 4946 4709 4947 4710 4948 4949 !------------------------------------------------------------------- 4950 !> rc_gstom: compute stomatal conductance 4951 !------------------------------------------------------------------- 4711 !------------------------------------------------------------------------------! 4712 ! Description: 4713 ! ------------ 4714 !> Subroutine to compute stomatal conductance 4715 !------------------------------------------------------------------------------! 4952 4716 SUBROUTINE rc_gstom( icmp, compnam, lu, lai_present, lai, glrad, sinphi, t, rh, diffc, gstom, p ) 4953 4717 4954 4718 IMPLICIT NONE 4955 4956 ! 4957 !-- input/output variables: 4958 CHARACTER(len=*), INTENT(IN) :: compnam !< component name 4719 ! 4720 !-- input/output variables: 4721 CHARACTER(LEN=*), INTENT(IN) :: compnam !< component name 4959 4722 4960 4723 INTEGER(iwp), INTENT(IN) :: icmp !< component index … … 4973 4736 4974 4737 REAL(wp), INTENT(OUT) :: gstom !< stomatal conductance (m/s) 4975 4976 ! 4977 !-- Local variables 4738 ! 4739 !-- Local variables 4978 4740 REAL(wp) :: vpd !< vapour pressure deficit (kPa) 4979 4741 4980 4742 REAL(wp), PARAMETER :: dO3 = 0.13e-4 !< diffusion coefficient of ozon (m2/s) 4981 4982 4743 ! 4983 4744 !-- Next line is to avoid compiler warning about unused variables 4984 4745 IF ( icmp == 0 ) CONTINUE 4985 4746 4986 SELECT CASE( TRIM( compnam) )4747 SELECT CASE( TRIM( compnam ) ) 4987 4748 4988 4749 CASE( 'NO', 'CO' ) 4989 4990 !--For no stomatal uptake is neglected:4750 ! 4751 !-- For no stomatal uptake is neglected: 4991 4752 gstom = 0.0_wp 4992 4753 4993 4754 CASE( 'NO2', 'O3', 'SO2', 'NH3' ) 4994 4995 ! 4996 !-- if vegetation present: 4755 ! 4756 !-- if vegetation present: 4997 4757 IF ( lai_present ) THEN 4998 4758 … … 5005 4765 ENDIF 5006 4766 ELSE 5007 5008 !--no vegetation; zero conductance (infinite resistance):4767 ! 4768 !-- no vegetation; zero conductance (infinite resistance): 5009 4769 gstom = 0.0_wp 5010 4770 ENDIF 5011 4771 5012 4772 CASE default 5013 message_string = 'Component '// TRIM( compnam) // ' not supported'4773 message_string = 'Component '// TRIM( compnam ) // ' not supported' 5014 4774 CALL message( 'rc_gstom', 'CM0459', 1, 2, 0, 6, 0 ) 5015 4775 END SELECT … … 5018 4778 5019 4779 5020 5021 !------------------------------------------------------------------- 5022 !> rc_gstom_emb: stomatal conductance according to Emberson 5023 !------------------------------------------------------------------- 4780 !------------------------------------------------------------------------------! 4781 ! Description: 4782 ! ------------ 4783 !> Subroutine to compute stomatal conductance according to Emberson 4784 !------------------------------------------------------------------------------! 5024 4785 SUBROUTINE rc_gstom_emb( lu, glrad, T, vpd, lai_present, lai, sinp, Gsto, p ) 5025 ! 5026 !-- History 5027 !> Original code from Lotos-Euros, TNO, M. Schaap 5028 !> 2009-08, M.C. van Zanten, Rivm 5029 !> Updated and extended. 5030 !> 2009-09, Arjo Segers, TNO 5031 !> Limitted temperature influence to range to avoid 5032 !> floating point exceptions. 5033 ! 5034 !> Method 5035 ! 5036 !> Code based on Emberson et al, 2000, Env. Poll., 403-413 5037 !> Notation conform Unified EMEP Model Description Part 1, ch 8 5038 ! 5039 !> In the calculation of f_light the modification of L. Zhang 2001, AE to the PARshade and PARsun 5040 !> parametrizations of Norman 1982 are applied 5041 ! 5042 !> f_phen and f_SWP are set to 1 5043 ! 5044 !> Land use types DEPAC versus Emberson (Table 5.1, EMEP model description) 5045 !> DEPAC Emberson 5046 !> 1 = grass GR = grassland 5047 !> 2 = arable land TC = temperate crops ( lai according to RC = rootcrops) 5048 !< 3 = permanent crops TC = temperate crops ( lai according to RC = rootcrops) 5049 !< 4 = coniferous forest CF = tempareate/boREAL(wp) coniferous forest 5050 !> 5 = deciduous forest DF = temperate/boREAL(wp) deciduous forest 5051 !> 6 = water W = water 5052 !> 7 = urban U = urban 5053 !> 8 = other GR = grassland 5054 !> 9 = desert DE = desert 4786 ! 4787 !> History 4788 !> Original code from Lotos-Euros, TNO, M. Schaap 4789 !> 2009-08, M.C. van Zanten, Rivm 4790 !> Updated and extended. 4791 !> 2009-09, Arjo Segers, TNO 4792 !> Limitted temperature influence to range to avoid 4793 !> floating point exceptions. 4794 4795 !> Method 4796 4797 !> Code based on Emberson et al, 2000, Env. Poll., 403-413 4798 !> Notation conform Unified EMEP Model Description Part 1, ch 8 4799 ! 4800 !> In the calculation of f_light the modification of L. Zhang 2001, AE to the PARshade and PARsun 4801 !> parametrizations of Norman 1982 are applied 4802 !> f_phen and f_SWP are set to 1 4803 ! 4804 !> Land use types DEPAC versus Emberson (Table 5.1, EMEP model description) 4805 !> DEPAC Emberson 4806 !> 1 = grass GR = grassland 4807 !> 2 = arable land TC = temperate crops ( lai according to RC = rootcrops) 4808 !> 3 = permanent crops TC = temperate crops ( lai according to RC = rootcrops) 4809 !> 4 = coniferous forest CF = tempareate/boREAL(wp) coniferous forest 4810 !> 5 = deciduous forest DF = temperate/boREAL(wp) deciduous forest 4811 !> 6 = water W = water 4812 !> 7 = urban U = urban 4813 !> 8 = other GR = grassland 4814 !> 9 = desert DE = desert 5055 4815 5056 4816 IMPLICIT NONE 5057 5058 ! 5059 !-- Emberson specific declarations 5060 5061 ! 5062 !-- Input/output variables: 4817 ! 4818 !-- Emberson specific declarations 4819 ! 4820 !-- Input/output variables: 5063 4821 INTEGER(iwp), INTENT(IN) :: lu !< land use type, lu = 1,...,nlu 5064 4822 … … 5075 4833 5076 4834 REAL(wp), INTENT(OUT) :: gsto !< stomatal conductance (m/s) 5077 5078 ! 5079 !-- Local variables: 4835 ! 4836 !-- Local variables: 5080 4837 REAL(wp) :: f_light 5081 4838 REAL(wp) :: f_phen … … 5094 4851 REAL(wp) :: pres 5095 4852 REAL(wp), PARAMETER :: p_sealevel = 1.01325e05 !< Pa 5096 5097 5098 ! 5099 !-- Check whether vegetation is present: 4853 ! 4854 !-- Check whether vegetation is present: 5100 4855 IF ( lai_present ) THEN 5101 4856 … … 5106 4861 sinphi = sinp 5107 4862 END IF 5108 5109 ! 5110 !-- ratio between actual and sea-level pressure is used 5111 !-- to correct for height in the computation of par; 5112 !-- should not exceed sea-level pressure therefore ... 4863 ! 4864 !-- ratio between actual and sea-level pressure is used 4865 !-- to correct for height in the computation of par; 4866 !-- should not exceed sea-level pressure therefore ... 5113 4867 IF ( present(p) ) THEN 5114 4868 pres = min( p, p_sealevel ) … … 5116 4870 pres = p_sealevel 5117 4871 ENDIF 5118 5119 ! 5120 !-- direct and diffuse par, Photoactive (=visible) radiation: 4872 ! 4873 !-- direct and diffuse par, Photoactive (=visible) radiation: 5121 4874 CALL par_dir_diff( glrad, sinphi, pres, p_sealevel, pardir, pardiff ) 5122 5123 ! 5124 !-- par for shaded leaves (canopy averaged): 4875 ! 4876 !-- par for shaded leaves (canopy averaged): 5125 4877 parshade = pardiff * exp( -0.5 * lai**0.7 ) + 0.07 * pardir * ( 1.1 - 0.1 * lai ) * exp( -sinphi ) !< Norman,1982 5126 4878 IF ( glrad > 200.0_wp .AND. lai > 2.5_wp ) THEN 5127 4879 parshade = pardiff * exp( -0.5 * lai**0.8 ) + 0.07 * pardir * ( 1.1 - 0.1 * lai ) * exp( -sinphi ) !< Zhang et al., 2001 5128 4880 END IF 5129 5130 ! 5131 !-- par for sunlit leaves (canopy averaged): 5132 !-- alpha -> mean angle between leaves and the sun is fixed at 60 deg -> i.e. cos alpha = 0.5 4881 ! 4882 !-- par for sunlit leaves (canopy averaged): 4883 !-- alpha -> mean angle between leaves and the sun is fixed at 60 deg -> i.e. cos alpha = 0.5 5133 4884 parsun = pardir * 0.5/sinphi + parshade !< Norman, 1982 5134 4885 IF ( glrad > 200.0_wp .AND. lai > 2.5_wp ) THEN 5135 4886 parsun = pardir**0.8 * 0.5 / sinphi + parshade !< Zhang et al., 2001 5136 4887 END IF 5137 5138 ! 5139 !-- leaf area index for sunlit and shaded leaves: 4888 ! 4889 !-- leaf area index for sunlit and shaded leaves: 5140 4890 IF ( sinphi > 0 ) THEN 5141 4891 laisun = 2 * sinphi * ( 1 - exp( -0.5 * lai / sinphi ) ) … … 5150 4900 5151 4901 f_light = MAX(f_light,f_min(lu)) 5152 5153 ! 5154 !-- temperature influence; only non-zero within range [tmin,tmax]: 4902 ! 4903 !-- temperature influence; only non-zero within range [tmin,tmax]: 5155 4904 IF ( ( tmin(lu) < t ) .AND. ( t < tmax(lu) ) ) THEN 5156 4905 bt = ( tmax(lu) - topt(lu) ) / ( topt(lu) - tmin(lu) ) … … 5160 4909 END IF 5161 4910 f_temp = MAX( f_temp, f_min(lu) ) 5162 5163 !vapour pressure deficit influence4911 ! 4912 !-- vapour pressure deficit influence 5164 4913 f_vpd = MIN( 1.0_wp, ( ( 1.0_wp - f_min(lu) ) * ( vpd_min(lu) - vpd ) / ( vpd_min(lu) - vpd_max(lu) ) + f_min(lu) ) ) 5165 4914 f_vpd = MAX( f_vpd, f_min(lu) ) 5166 4915 5167 4916 f_swp = 1.0_wp 5168 5169 !influence of phenology on stom. conductance5170 !ignored for now in DEPAC since influence of f_phen on lu classes in use is negligible.5171 !When other EMEP classes (e.g. med. broadleaf) are used f_phen might be too important to ignore4917 ! 4918 !-- influence of phenology on stom. conductance 4919 !-- ignored for now in DEPAC since influence of f_phen on lu classes in use is negligible. 4920 !-- When other EMEP classes (e.g. med. broadleaf) are used f_phen might be too important to ignore 5172 4921 f_phen = 1.0_wp 5173 5174 !evaluate total stomatal conductance4922 ! 4923 !-- evaluate total stomatal conductance 5175 4924 f_env = f_temp * f_vpd * f_swp 5176 4925 f_env = MAX( f_env,f_min(lu) ) 5177 4926 gsto = g_max(lu) * f_light * f_phen * f_env 5178 5179 !gstom expressed per m2 leafarea;5180 !this is converted with lai to m2 surface.4927 ! 4928 !-- gstom expressed per m2 leafarea; 4929 !-- this is converted with lai to m2 surface. 5181 4930 gsto = lai * gsto ! in m/s 5182 4931 … … 5188 4937 5189 4938 5190 5191 4939 !------------------------------------------------------------------- 5192 4940 !> par_dir_diff 5193 !5194 4941 !> Weiss, A., Norman, J.M. (1985) Partitioning solar radiation into direct and 5195 4942 !> diffuse, visible and near-infrared components. Agric. Forest Meteorol. 5196 4943 !> 34, 205-213. 5197 !5198 4944 !> From a SUBROUTINE obtained from Leiming Zhang, 5199 4945 !> Meteorological Service of Canada 5200 !5201 4946 !> Leiming uses solar irradiance. This should be equal to global radiation and 5202 4947 !> Willem Asman set it to global radiation … … 5228 4973 REAL(wp) :: ww !< water absorption in the near infrared for 10 mm of precipitable water 5229 4974 5230 5231 5232 ! 5233 !-- Calculate visible (PAR) direct beam radiation 5234 !-- 600 W m-2 represents average amount of par (400-700 nm wavelength) 5235 !-- at the top of the atmosphere; this is roughly 0.45*solar constant (solar constant=1320 Wm-2) 4975 ! 4976 !-- Calculate visible (PAR) direct beam radiation 4977 !-- 600 W m-2 represents average amount of par (400-700 nm wavelength) 4978 !-- at the top of the atmosphere; this is roughly 0.45*solar constant (solar constant=1320 Wm-2) 5236 4979 rdu = 600.0_wp* exp( -0.185_wp * ( pres / pres_0 ) / sinphi ) * sinphi 5237 5238 ! 5239 !-- Calculate potential visible diffuse radiation 4980 ! 4981 !-- Calculate potential visible diffuse radiation 5240 4982 rdv = 0.4_wp * ( 600.0_wp - rdu ) * sinphi 5241 5242 ! 5243 !-- Calculate the water absorption in the-near infrared 4983 ! 4984 !-- Calculate the water absorption in the-near infrared 5244 4985 ww = 1320 * 10**( -1.195_wp + 0.4459_wp * log10( 1.0_wp / sinphi ) - 0.0345_wp * ( log10( 1.0_wp / sinphi ) )**2 ) 5245 5246 ! 5247 !-- Calculate potential direct beam near-infrared radiation 4986 ! 4987 !-- Calculate potential direct beam near-infrared radiation 5248 4988 rdm = (720.0_wp * exp(-0.06_wp * (pres / pres_0) / sinphi ) - ww ) * sinphi !< 720 = solar constant - 600 5249 5250 ! 5251 !-- Calculate potential diffuse near-infrared radiation 4989 ! 4990 !-- Calculate potential diffuse near-infrared radiation 5252 4991 rdn = 0.6_wp * ( 720 - rdm - ww ) * sinphi 5253 5254 ! 5255 !-- Compute visible and near-infrared radiation 4992 ! 4993 !-- Compute visible and near-infrared radiation 5256 4994 rv = MAX( 0.1_wp, rdu + rdv ) 5257 4995 rn = MAX( 0.01_wp, rdm + rdn ) 5258 5259 ! 5260 !-- Compute ratio between input global radiation and total radiation computed here 4996 ! 4997 !-- Compute ratio between input global radiation and total radiation computed here 5261 4998 ratio = MIN( 0.89_wp, glrad / ( rv + rn ) ) 5262 5263 ! 5264 !-- Calculate total visible radiation 4999 ! 5000 !-- Calculate total visible radiation 5265 5001 sv = ratio * rv 5266 5267 ! 5268 !-- Calculate fraction of par in the direct beam 5002 ! 5003 !-- Calculate fraction of par in the direct beam 5269 5004 fv = MIN( 0.99_wp, ( 0.9_wp - ratio ) / 0.7_wp ) !< help variable 5270 5005 fv = MAX( 0.01_wp, rdu / rv * ( 1.0_wp - fv**0.6667_wp ) ) !< fraction of par in the direct beam 5271 5272 ! 5273 !-- Compute direct and diffuse parts of par 5006 ! 5007 !-- Compute direct and diffuse parts of par 5274 5008 par_dir = fv * sv 5275 5009 par_diff = sv - par_dir … … 5284 5018 5285 5019 IMPLICIT NONE 5286 5287 ! 5288 !-- Input/output variables: 5020 ! 5021 !-- Input/output variables: 5289 5022 REAL(wp), INTENT(IN) :: temp !< temperature (C) 5290 5023 REAL(wp), INTENT(IN) :: relh !< relative humidity (%) 5291 5024 5292 5025 REAL(wp), INTENT(OUT) :: vpd !< vapour pressure deficit (kPa) 5293 5294 ! 5295 !-- Local variables: 5026 ! 5027 !-- Local variables: 5296 5028 REAL(wp) :: esat 5297 5298 ! 5299 !-- fit parameters: 5029 ! 5030 !-- fit parameters: 5300 5031 REAL(wp), PARAMETER :: a1 = 6.113718e-01 5301 5032 REAL(wp), PARAMETER :: a2 = 4.43839e-02 … … 5304 5035 REAL(wp), PARAMETER :: a5 = 2.16e-07 5305 5036 REAL(wp), PARAMETER :: a6 = 3.0e-09 5306 5307 ! 5308 !-- esat is saturation vapour pressure (kPa) at temp(C) following Monteith(1973) 5037 ! 5038 !-- esat is saturation vapour pressure (kPa) at temp(C) following Monteith(1973) 5309 5039 esat = a1 + a2 * temp + a3 * temp**2 + a4 * temp**3 + a5 * temp**4 + a6 * temp**5 5310 5040 vpd = esat * ( 1 - relh / 100 ) 5311 5041 5312 5042 END SUBROUTINE rc_get_vpd 5313 5314 5043 5315 5044 … … 5320 5049 5321 5050 IMPLICIT NONE 5322 5323 ! 5324 !-- Input/output variables: 5051 ! 5052 !-- Input/output variables: 5325 5053 INTEGER(iwp), INTENT(IN) :: icmp !< component index 5326 5054 INTEGER(iwp), INTENT(IN) :: lu !< land use type, lu = 1,..., nlu … … 5329 5057 !< N.B. this routine cannot be called with nwet = 9, 5330 5058 !< nwet = 9 should be handled outside this routine. 5331 5332 5059 REAL(wp), INTENT(IN) :: sai !< surface area index 5333 5060 REAL(wp), INTENT(IN) :: ust !< friction velocity (m/s) 5334 5061 REAL(wp), INTENT(IN) :: t !< temperature (C) 5335 5336 5062 REAL(wp), INTENT(OUT) :: gsoil_eff !< effective soil conductance (m/s) 5337 5338 ! 5339 !-- local variables: 5063 ! 5064 !-- local variables: 5340 5065 REAL(wp) :: rinc !< in canopy resistance (s/m) 5341 5066 REAL(wp) :: rsoil_eff !< effective soil resistance (s/m) 5342 5343 5344 ! 5345 !-- Soil resistance (numbers matched with lu_classes and component numbers) 5067 ! 5068 !-- Soil resistance (numbers matched with lu_classes and component numbers) 5346 5069 ! grs ara crp cnf dec wat urb oth des ice sav trf wai med sem 5347 5070 REAL(wp), PARAMETER :: rsoil(nlu_dep,ncmp) = reshape( (/ & … … 5357 5080 -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999. /),& !< H2O2 5358 5081 (/nlu_dep,ncmp/) ) 5359 5360 ! 5361 !-- For o3 so2 no2 no nh3 co no3 hno3 n2o5 h2o2 5082 ! 5083 !-- For o3 so2 no2 no nh3 co no3 hno3 n2o5 h2o2 5362 5084 REAL(wp), PARAMETER :: rsoil_wet(ncmp) = (/2000., 10. , 2000., -999., 10. , -999., -999., -999., -999., -999./) 5363 5085 REAL(wp), PARAMETER :: rsoil_frozen(ncmp) = (/2000., 500., 2000., -999., 1000., -999., -999., -999., -999., -999./) 5364 5365 5366 5367 ! 5368 !-- Compute in canopy (in crop) resistance: 5086 ! 5087 !-- Compute in canopy (in crop) resistance: 5369 5088 CALL rc_rinc( lu, sai, ust, rinc ) 5370 5371 ! 5372 !-- Check for missing deposition path: 5089 ! 5090 !-- Check for missing deposition path: 5373 5091 IF ( missing(rinc) ) THEN 5374 5092 rsoil_eff = -9999.0_wp 5375 5093 ELSE 5376 5377 !--Frozen soil (temperature below 0):5094 ! 5095 !-- Frozen soil (temperature below 0): 5378 5096 IF ( t < 0.0_wp ) THEN 5379 5097 IF ( missing( rsoil_frozen( icmp ) ) ) THEN … … 5383 5101 ENDIF 5384 5102 ELSE 5385 5386 !--Non-frozen soil; dry:5103 ! 5104 !-- Non-frozen soil; dry: 5387 5105 IF ( nwet == 0 ) THEN 5388 5106 IF ( missing( rsoil( lu, icmp ) ) ) THEN … … 5391 5109 rsoil_eff = rsoil( lu, icmp ) + rinc 5392 5110 ENDIF 5393 5394 ! 5395 !-- Non-frozen soil; wet: 5111 ! 5112 !-- Non-frozen soil; wet: 5396 5113 ELSEIF ( nwet == 1 ) THEN 5397 5114 IF ( missing( rsoil_wet( icmp ) ) ) THEN … … 5406 5123 ENDIF 5407 5124 ENDIF 5408 5409 ! 5410 !-- Compute conductance: 5125 ! 5126 !-- Compute conductance: 5411 5127 IF ( rsoil_eff > 0.0_wp ) THEN 5412 5128 gsoil_eff = 1.0_wp / rsoil_eff … … 5426 5142 IMPLICIT NONE 5427 5143 5428 5429 5144 ! 5145 !-- Input/output variables: 5430 5146 INTEGER(iwp), INTENT(IN) :: lu !< land use class, lu = 1, ..., nlu 5431 5147 … … 5434 5150 5435 5151 REAL(wp), INTENT(OUT) :: rinc !< in canopy resistance (s/m) 5436 5437 ! 5438 !-- b = empirical constant for computation of rinc (in canopy resistance) (= 14 m-1 or -999 if not applicable) 5439 !-- h = vegetation height (m) gra ara crop con dec wat urb oth des ice sav trf wai med semi 5152 ! 5153 !-- b = empirical constant for computation of rinc (in canopy resistance) (= 14 m-1 or -999 if not applicable) 5154 !-- h = vegetation height (m) gra ara crop con dec wat urb oth des ice sav trf wai med semi 5440 5155 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: b = (/ -999, 14, 14, 14, 14, -999, -999, -999, -999, -999, -999, 14, -999, & 5441 5156 14, 14 /) 5442 5157 REAL(wp), DIMENSION(nlu_dep), PARAMETER :: h = (/ -999, 1, 1, 20, 20, -999, -999, -999, -999, -999, -999, 20, -999, & 5443 5158 1 , 1 /) 5444 5445 ! 5446 !-- Compute Rinc only for arable land, perm. crops, forest; otherwise Rinc = 0: 5159 ! 5160 !-- Compute Rinc only for arable land, perm. crops, forest; otherwise Rinc = 0: 5447 5161 IF ( b(lu) > 0.0_wp ) THEN 5448 !5449 !--Check for u* > 0 (otherwise denominator = 0):5162 ! ! 5163 !-- Check for u* > 0 (otherwise denominator = 0): 5450 5164 IF ( ust > 0.0_wp ) THEN 5451 5165 rinc = b(lu) * h(lu) * sai/ust … … 5464 5178 5465 5179 5466 5467 5180 !------------------------------------------------------------------- 5468 5181 !> rc_rctot: compute total canopy (or surface) resistance Rc … … 5472 5185 IMPLICIT NONE 5473 5186 5474 5475 5187 ! 5188 !-- Input/output variables: 5476 5189 REAL(wp), INTENT(IN) :: gstom !< stomatal conductance (s/m) 5477 5190 REAL(wp), INTENT(IN) :: gsoil_eff !< effective soil conductance (s/m) … … 5480 5193 REAL(wp), INTENT(OUT) :: gc_tot !< total canopy conductance (m/s) 5481 5194 REAL(wp), INTENT(OUT) :: rc_tot !< total canopy resistance Rc (s/m) 5482 5483 ! 5484 !-- Total conductance: 5195 ! 5196 !-- Total conductance: 5485 5197 gc_tot = gstom + gsoil_eff + gw 5486 5487 ! 5488 !-- Total resistance (note: gw can be negative, but no total emission allowed here): 5198 ! 5199 !-- Total resistance (note: gw can be negative, but no total emission allowed here): 5489 5200 IF ( gc_tot <= 0.0_wp .OR. gw < 0.0_wp ) THEN 5490 5201 rc_tot = -9999.0_wp … … 5496 5207 5497 5208 5498 5499 5209 !------------------------------------------------------------------- 5500 5210 !> rc_comp_point_rc_eff: calculate the effective resistance Rc 5501 5211 !> based on one or more compensation points 5502 5212 !------------------------------------------------------------------- 5503 !5504 5213 !> NH3rc (see depac v3.6 is based on Avero workshop Marc Sutton. p. 173. 5505 5214 !> Sutton 1998 AE 473-480) … … 5558 5267 !> 5559 5268 ! ------------------------------------------------------------------------------------------- 5560 5561 5269 ! SUBROUTINE rc_comp_point_rc_eff( ccomp_tot, catm, ra, rb, rc_tot, rc_eff ) 5562 5270 ! 5563 5271 ! IMPLICIT NONE 5564 5272 ! 5565 ! ! 5566 ! !-- Input/output variables: 5273 !!-- Input/output variables: 5567 5274 ! REAL(wp), INTENT(IN) :: ccomp_tot !< total compensation point (weighed average of separate compensation points) (ug/m3) 5568 5275 ! REAL(wp), INTENT(IN) :: catm !< atmospheric concentration (ug/m3) … … 5574 5281 ! 5575 5282 ! ! 5576 ! 5283 !!-- Compute effective resistance: 5577 5284 ! IF ( ccomp_tot == 0.0_wp ) THEN 5578 5285 ! ! 5579 ! !--trace with no compensiation point ( or compensation point equal to zero)5286 !!-- trace with no compensiation point ( or compensation point equal to zero) 5580 5287 ! rc_eff = rc_tot 5581 5288 ! 5582 5289 ! ELSE IF ( ccomp_tot > 0.0_wp .AND. ( abs( catm - ccomp_tot ) < 1.e-8 ) ) THEN 5583 5290 ! ! 5584 ! !--surface concentration (almost) equal to atmospheric concentration5585 ! !--no exchange between surface and atmosphere, infinite RC --> vd=05291 !!-- surface concentration (almost) equal to atmospheric concentration 5292 !!-- no exchange between surface and atmosphere, infinite RC --> vd=0 5586 5293 ! rc_eff = 9999999999.0_wp 5587 5294 ! 5588 5295 ! ELSE IF ( ccomp_tot > 0.0_wp ) THEN 5589 5296 ! ! 5590 ! !--compensation point available, calculate effective resistance5297 !!-- compensation point available, calculate effective resistance 5591 5298 ! rc_eff = ( ( ra + rb ) * ccomp_tot + rc_tot * catm ) / ( catm - ccomp_tot ) 5592 5299 ! … … 5602 5309 5603 5310 5604 5605 5311 !------------------------------------------------------------------- 5606 5312 !> missing: check for data that correspond with a missing deposition path 5607 5313 !> this data is represented by -999 5608 5314 !------------------------------------------------------------------- 5609 5610 5315 LOGICAL function missing( x ) 5611 5316 5612 5317 REAL(wp), INTENT(IN) :: x 5613 5318 5614 5615 5319 ! 5320 !-- bandwidth for checking (in)equalities of floats 5616 5321 REAL(wp), PARAMETER :: eps = 1.0e-5 5617 5322 … … 5621 5326 5622 5327 5623 5624 5328 ELEMENTAL FUNCTION sedimentation_velocity( rhopart, partsize, slipcor, visc ) RESULT( vs ) 5625 5329 5626 5627 5330 IMPLICIT NONE 5628 5331 5629 5630 5332 ! 5333 !-- in/out 5631 5334 5632 5335 REAL(wp), INTENT(IN) :: rhopart !< particle density (kg/m3) … … 5636 5339 5637 5340 REAL(wp) :: vs 5638 5639 ! 5640 !-- acceleration of gravity: 5341 ! 5342 !-- acceleration of gravity: 5641 5343 REAL(wp), PARAMETER :: grav = 9.80665 !< m/s2 5642 5344 5643 5644 5645 ! 5646 !-- sedimentation velocity 5345 !-- sedimentation velocity 5647 5346 vs = rhopart * ( partsize**2.0_wp ) * grav * slipcor / ( 18.0_wp * visc ) 5648 5347 … … 5656 5355 luc, ftop_lu, ustar_lu ) 5657 5356 5658 5659 5357 IMPLICIT NONE 5660 5358 5661 ! -- in/out 5359 ! 5360 !-- in/out 5662 5361 5663 5362 INTEGER(iwp), INTENT(IN) :: nwet !< 1=rain, 9=snowcover … … 5675 5374 REAL(wp), INTENT(OUT) :: vd !< deposition velocity (m/s) 5676 5375 REAL(wp), INTENT(OUT) :: Rs !< sedimentaion resistance (s/m) 5677 5678 5679 ! 5680 !-- constants 5376 ! 5377 !-- constants 5681 5378 5682 5379 REAL(wp), PARAMETER :: grav = 9.80665 !< acceleration of gravity (m/s2) … … 5693 5390 REAL(wp), PARAMETER ::A_lu(nlu_dep) = & 5694 5391 (/3.0, 3.0, 2.0, 2.0, 7.0, -99., 10.0, 3.0, -99., -99., 3.0, 7.0, -99., 2.0, -99./) 5695 ! 5696 !-- grass arabl crops conif decid water urba othr desr ice sav trf wai med sem 5697 5698 ! 5699 !-- local 5700 5392 ! 5393 !-- grass arabl crops conif decid water urba othr desr ice sav trf wai med sem 5394 ! 5395 !-- local 5701 5396 REAL(wp) :: kinvisc 5702 5397 REAL(wp) :: diff_part … … 5707 5402 REAL(wp) :: Einterc 5708 5403 REAL(wp) :: Reffic 5709 5710 5711 ! 5712 !-- kinetic viscosity & diffusivity 5404 ! 5405 !-- kinetic viscosity & diffusivity 5713 5406 kinvisc = viscos1 / dens1 !< only needed at surface 5714 5407 5715 5408 diff_part = kb * tsurf * slipcor / ( 3 * pi * viscos1 * partsize ) 5716 5717 ! 5718 !-- Schmidt number 5409 ! 5410 !-- Schmidt number 5719 5411 schmidt = kinvisc / diff_part 5720 5721 ! 5722 !-- calculate collection efficiencie E 5412 ! 5413 !-- calculate collection efficiencie E 5723 5414 Ebrown = Schmidt**( -gamma_lu(luc) ) !< Brownian diffusion 5724 5725 5726 5415 ! 5416 !-- determine Stokes number, interception efficiency 5417 !-- and sticking efficiency R (1 = no rebound) 5727 5418 IF ( luc == ilu_ice .OR. nwet==9 .OR. luc == ilu_water_sea .OR. luc == ilu_water_inland ) THEN 5728 5419 stokes = vs1 * ustar_lu**2 / ( grav * kinvisc ) … … 5738 5429 Reffic = exp( -Stokes**0.5_wp ) 5739 5430 END IF 5740 5741 5431 ! 5432 !-- when surface is wet all particles do not rebound: 5742 5433 IF ( nwet==1 ) Reffic = 1.0_wp 5743 5744 5434 ! 5435 !-- determine impaction efficiency: 5745 5436 Eimpac = ( stokes / ( alfa_lu(luc) + stokes ) )**beta 5746 5747 ! 5748 !-- sedimentation resistance: 5437 ! 5438 !-- sedimentation resistance: 5749 5439 Rs = 1.0_wp / ( epsilon0 * ustar_lu * ( Ebrown + Eimpac + Einterc ) * Reffic ) 5750 5440 5751 ! 5752 !-- deposition velocity according to Seinfeld and Pandis (2006; eq 19.7): 5753 ! 5754 ! 1 5755 ! vd = ------------------ + vs 5756 ! Ra + Rs + Ra*Rs*vs 5757 ! 5758 !-- where: Rs = Rb (in Seinfeld and Pandis, 2006) 5759 ! 5760 ! 5441 !-- deposition velocity according to Seinfeld and Pandis (2006; eq 19.7): 5442 !-- 5443 !-- 1 5444 !-- vd = ------------------ + vs 5445 !-- Ra + Rs + Ra*Rs*vs 5446 !-- 5447 !-- where: Rs = Rb (in Seinfeld and Pandis, 2006) 5761 5448 5762 5449 vd = 1.0_wp / ( ftop_lu + Rs + ftop_lu * Rs * vs1) + vs1 … … 5772 5459 SUBROUTINE get_rb_cell( is_water, z0h, ustar, diffc, Rb ) 5773 5460 5774 5461 5775 5462 IMPLICIT NONE 5776 5463 5777 !-- in/out 5464 ! 5465 !-- in/out 5778 5466 5779 5467 LOGICAL , INTENT(IN) :: is_water … … 5784 5472 5785 5473 REAL(wp), INTENT(OUT) :: Rb !< boundary layer resistance 5786 5787 ! 5788 !-- const 5474 ! 5475 !-- const 5789 5476 5790 5477 REAL(wp), PARAMETER :: thk = 0.19e-4 !< thermal diffusivity of dry air 20 C 5791 5478 REAL(wp), PARAMETER :: kappa_stab = 0.35 !< von Karman constant 5792 5793 5479 ! 5794 5480 !-- Next line is to avoid compiler warning about unused variable 5795 5481 IF ( is_water .OR. ( z0h + kappa_stab ) > 0.0_wp ) CONTINUE 5796 5797 ! 5798 !-- Use Simpson et al. (2003) 5799 !-- @TODO: Check Rb over water calculation, until then leave commented lines 5800 !IF ( is_water ) THEN 5801 ! org: Rb = 1.0_wp / (kappa_stab*MAX(0.01_wp,ustar)) * log(z0h/diffc*kappa_stab*MAX(0.01_wp,ustar)) 5802 ! Rb = 1.0_wp / (kappa_stab*MAX(0.1_wp,ustar)) * log(z0h/diffc*kappa_stab*MAX(0.1_wp,ustar)) 5803 !ELSE 5482 ! 5483 !-- Use Simpson et al. (2003) 5484 !-- @TODO: Check Rb over water calculation, until then leave commented lines 5485 !-- IF ( is_water ) THEN 5486 !-- org: Rb = 1.0_wp / (kappa_stab*MAX(0.01_wp,ustar)) * log(z0h/diffc*kappa_stab*MAX(0.01_wp,ustar)) 5487 !-- Rb = 1.0_wp / (kappa_stab*MAX(0.1_wp,ustar)) * log(z0h/diffc*kappa_stab*MAX(0.1_wp,ustar)) 5488 !-- ELSE 5804 5489 Rb = 5.0_wp / MAX( 0.01_wp, ustar ) * ( thk / diffc )**0.67_wp 5805 !END IF5490 !-- END IF 5806 5491 5807 5492 END SUBROUTINE get_rb_cell … … 5835 5520 IMPLICIT NONE 5836 5521 5837 5838 5522 ! 5523 !-- in/out 5839 5524 5840 5525 REAL(wp), INTENT(IN) :: q !< specific humidity [(kg water)/(kg air)] … … 5842 5527 5843 5528 REAL(wp) :: p_w !< water vapor partial pressure [Pa] 5844 5845 ! 5846 !-- const 5529 ! 5530 !-- const 5847 5531 5848 5532 REAL(wp), PARAMETER :: eps = xm_h2o / xm_air !< mole mass ratio ~ 0.622 5849 5850 ! 5851 !-- partial pressure of water vapor: 5533 ! 5534 !-- partial pressure of water vapor: 5852 5535 p_w = p * q / eps 5853 5536 5854 5537 END function watervaporpartialpressure 5855 5856 5538 5857 5539 … … 5880 5562 IMPLICIT NONE 5881 5563 5882 5883 5564 ! 5565 !-- in/out 5884 5566 5885 5567 REAL(wp), INTENT(IN) :: t !< temperature [K] 5886 5568 5887 5569 REAL(wp) :: e_sat_w !< saturation vapor pressure [Pa] 5888 5889 ! 5890 !-- const 5570 ! 5571 !-- const 5891 5572 REAL(wp), PARAMETER :: p0 = 611.2 !< base pressure [Pa] 5892 5893 ! 5894 !-- saturation vapor pressure: 5573 ! 5574 !-- saturation vapor pressure: 5895 5575 e_sat_w = p0 * exp( 17.67_wp * ( t - 273.16_wp ) / ( t - 29.66_wp ) ) !< [Pa] 5896 5576 5897 5577 END FUNCTION saturationvaporpressure 5898 5899 5578 5900 5579 … … 5912 5591 IMPLICIT NONE 5913 5592 5914 5915 5593 ! 5594 !-- in/out 5916 5595 5917 5596 REAL(wp), INTENT(IN) :: q !< specific humidity [(kg water)/(kg air)] … … 5920 5599 5921 5600 REAL(wp) :: rh !< relative humidity [%] 5922 5923 ! 5924 !-- relative humidity: 5601 ! 5602 !-- relative humidity: 5925 5603 rh = watervaporpartialpressure( q, p ) / saturationvaporpressure( t ) * 100.0_wp 5926 5604 5927 5605 END FUNCTION relativehumidity_from_specifichumidity 5928 5929 5606 5930 5607
Note: See TracChangeset
for help on using the changeset viewer.
