- Timestamp:
- Feb 19, 2018 11:29:57 AM (7 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 2 edited
-
chemistry_model_mod.f90 (modified) (17 diffs)
-
prognostic_equations.f90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/chemistry_model_mod.f90
r2773 r2815 27 27 ! ----------------- 28 28 ! $Id$ 29 ! Bugfix in restart mechanism, 30 ! rename chem_tendency to chem_prognostic_equations, 31 ! implement vector-optimized version of chem_prognostic_equations, 32 ! some clean up (incl. todo list) 33 ! 34 ! 2773 2018-01-30 14:12:54Z suehring 29 35 ! Declare variables required for nesting as public 30 36 ! … … 60 66 !> Chemistry model for PALM-4U 61 67 !> @todo Update/clean-up todo list! (FK) 62 !> @todo Chemistry for prognostic_equations_vector! (FK)63 68 !> @todo Set proper fill values (/= 0) for chem output arrays! (FK) 64 69 !> @todo Add routine chem_check_parameters, add checks for inconsistent or 65 !> unallowed parameter settings (e.g. "only cache optimized version66 !> available").CALL of chem_check_parameters from check_parameters. (FK)70 !> unallowed parameter settings. 71 !> CALL of chem_check_parameters from check_parameters. (FK) 67 72 !> @todo Make routine chem_header available, CALL from header.f90 68 73 !> (see e.g. how it is done in routine lsm_header in … … 76 81 !> @todo slight differences in passive scalar and chem spcs when chem reactions 77 82 !> turned off. Need to be fixed. bK 78 !> @todo introduce nesting for chem spcs. bK83 !> @todo test nesting for chem spcs, was implemented by suehring (kanani) 79 84 !> @todo subroutine set_const_initial_values to be taken out from chemistry_model_mod !bK. 80 85 !> @todo chemistry error messages … … 211 216 END INTERFACE chem_read_restart_data 212 217 213 INTERFACE chem_tendency 214 MODULE PROCEDURE chem_tendency 215 END INTERFACE chem_tendency 218 INTERFACE chem_prognostic_equations 219 MODULE PROCEDURE chem_prognostic_equations 220 MODULE PROCEDURE chem_prognostic_equations_ij 221 END INTERFACE chem_prognostic_equations 216 222 217 223 INTERFACE chem_header … … 236 242 237 243 238 PUBLIC chem_ boundary_conds, chem_check_data_output,&244 PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_check_data_output, & 239 245 chem_check_data_output_pr, chem_data_output_3d, & 240 chem_define_netcdf_grid, chem_ header, chem_init,&241 chem_init_profiles, chem_integrate, chem_ parin,&242 chem_ swap_timelevel, chem_last_actions, chem_read_restart_data,&243 chem_ tendency, chem_3d_data_averaging, chem_emissions!, chem_write_var_list, &244 !chem_read_var_list, chem_skip_var_list246 chem_define_netcdf_grid, chem_emissions, chem_header, chem_init, & 247 chem_init_profiles, chem_integrate, chem_last_actions, & 248 chem_parin, chem_prognostic_equations, & 249 chem_read_restart_data, chem_swap_timelevel 250 245 251 246 252 … … 293 299 ELSE 294 300 ! message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"' ! bK commented 295 CALL message( 'chem_boundary_conds', 'C HEM001', 1, 2, 0, 6, 0 ) !< chemistry_model_mod should have special error numbers --> "CHEM###",301 CALL message( 'chem_boundary_conds', 'CM0010', 1, 2, 0, 6, 0 ) !< chemistry_model_mod should have special error numbers --> "CHEM###", 296 302 ENDIF 297 303 ! … … 308 314 ELSE 309 315 ! message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"' 310 CALL message( 'check_parameters', 'C HEM002', 1, 2, 0, 6, 0 )316 CALL message( 'check_parameters', 'CM0011', 1, 2, 0, 6, 0 ) 311 317 ENDIF 312 318 … … 449 455 ! IF ( cs_heights(lsp,1) /= 0.0_wp ) THEN 450 456 ! message_string = 'cs_heights(1,1) must be 0.0' 451 ! CALL message( 'chem_check_parameters', 'C HEM0004', 1, 2, 0, 6, 0 )457 ! CALL message( 'chem_check_parameters', 'CM0012', 1, 2, 0, 6, 0 ) 452 458 ! ENDIF 453 459 ! … … 658 664 chem_species(lsp)%conc_p = chem_species(lsp)%conc 659 665 ENDDO 660 ! CALL location_message( 'finished', .TRUE. ) 661 662 ! 663 !-- (todo (FK): Restarts not available yet. It is not correct to call 664 !-- read_3d_binary here) 665 ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' .OR. & 666 TRIM( initializing_actions ) == 'cyclic_fill' ) & 667 THEN 668 669 message_string = 'so far, restarts are not possible with chemistry' 670 CALL message( 'chem_init', 'CM0002', 1, 2, 0, 6, 0 ) 671 ! CALL location_message( 'initializing in case of restart / cyclic_fill', & 672 ! .FALSE. ) 673 ! ! 674 ! !-- When reading data for cyclic fill of 3D prerun data files, read 675 ! !-- some of the global variables from the restart file which are required 676 ! !-- for initializing the inflow 677 ! IF ( TRIM( initializing_actions ) == 'cyclic_fill' ) THEN 678 ! DO i = 0, io_blocks-1 679 ! IF ( i == io_group ) THEN 680 ! CALL read_parts_of_var_list 681 ! CALL close_file( 13 ) 682 ! ENDIF 683 ! #if defined( __parallel ) 684 ! CALL MPI_BARRIER( comm2d, ierr ) 685 ! #endif 686 ! ENDDO 687 ! ENDIF 688 ! ! 689 ! !-- Read binary data from restart file 690 ! DO i = 0, io_blocks-1 691 ! IF ( i == io_group ) THEN 692 ! CALL read_3d_binary 693 ! ENDIF 694 ! #if defined( __parallel ) 695 ! CALL MPI_BARRIER( comm2d, ierr ) 696 ! #endif 697 ! ENDDO 698 ! ! 699 ! !-- Initialization of the turbulence recycling method 700 ! IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. & 701 ! turbulent_inflow ) THEN 702 ! ! 703 ! !-- First store the profiles to be used at the inflow. 704 ! !-- These profiles are the (temporally) and horizontally averaged vertical 705 ! !-- profiles from the prerun. Alternatively, prescribed profiles 706 ! !-- for u,v-components can be used. 707 ! ALLOCATE( mean_inflow_profiles(nzb:nzt+1,8) ) 708 ! ! mean_inflow_profiles(:,8) = hom_sum(:,115,0) !cs 709 ! 710 ! ! IF ( inflow_l) THEN 711 ! ! DO j = nysg, nyng 712 ! ! DO k = nzb, nzt+1 713 ! ! DO lsp = 1, nvar 714 ! ! chem_species(lsp)%conc(k,j,nxlg:-1) = mean_inflow_profiles(k,8) ??? 715 ! ! ENDDO 716 ! ! ENDDO 717 ! ! ENDDO 718 ! ! ENDIF 719 ! ENDIF 720 ! 721 ! !-- Inside buildings set velocities and TKE back to zero 722 ! ! IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. & 723 ! ! topography /= 'flat' ) THEN 724 ! ! ...... 725 ! ! ENDIF 726 ! ! 727 ! !-- Initialize new time levels (only done in order to set boundary values 728 ! !-- including ghost points) 729 ! DO lsp = 1, nvar 730 ! chem_species(lsp)%conc_p = chem_species(lsp)%conc 731 ! ENDDO 732 ! ! 733 ! !-- Allthough tendency arrays are set in prognostic_equations, they have 734 ! !-- have to be predefined here because they are used (but multiplied with 0) 735 ! !-- there before they are set. 736 ! 737 ! DO lsp = 1, nvar 738 ! chem_species(lsp)%tconc_m = 0.0_wp 739 ! ENDDO 740 741 ELSE 742 ! 743 !-- Actually this part of the programm should not be reached 744 message_string = 'unknown initializing problem' 745 CALL message( 'chem_init', 'CM0003', 1, 2, 0, 6, 0 ) 666 746 667 ENDIF 747 668 … … 1397 1318 IMPLICIT NONE 1398 1319 1399 INTEGER(iwp) :: lsp !<1320 INTEGER(iwp) :: lsp !< 1400 1321 CHARACTER(LEN=20) :: cspcs_name 1401 1322 CHARACTER(LEN=20) :: cspcs_name_av 1402 1323 ! REAL(kind=wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: chems_conc 1403 1324 1325 1404 1326 IF ( write_binary ) THEN 1405 1327 DO lsp = 1, nspec … … 1428 1350 SUBROUTINE chem_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file, & 1429 1351 nynfa, nyn_on_file, nysfa, nys_on_file, & 1430 offset_xa, offset_ya, overlap_count, tmp_2d, tmp_3d) 1352 offset_xa, offset_ya, overlap_count, & 1353 tmp_2d, tmp_3d ) 1431 1354 1432 1355 USE control_parameters … … 1464 1387 INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nysfa !< 1465 1388 INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_xa !< 1466 INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_ya !< 1389 INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_ya !< 1390 1391 LOGICAL :: chem_found 1467 1392 !! 1468 REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp, & !< 2D array to temp store data 1469 nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d !< 1470 REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp, & !< 1471 nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !< 3D array to temp store data 1393 REAL(wp), & 1394 DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::& 1395 tmp_2d !< 2D array to temp store data 1396 1397 REAL(wp), & 1398 DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::& 1399 tmp_3d !< 3D array to temp store data 1472 1400 1473 1401 1474 1402 IF ( initializing_actions == 'read_restart_data' ) THEN 1475 1403 READ ( 13 ) field_char 1476 1477 DO WHILE ( TRIM( field_char ) /= '*** end chem ***') 1478 DO k = 1, overlap_count 1404 DO WHILE ( TRIM( field_char ) /= '*** end chem *** ' ) 1405 1406 DO k = 1, overlap_count 1407 1479 1408 nxlf = nxlfa(i,k) 1480 1409 nxlc = nxlfa(i,k) + offset_xa(i,k) … … 1485 1414 nynf = nynfa(i,k) 1486 1415 nync = nynfa(i,k) + offset_ya(i,k) 1487 1416 1417 1418 chem_found = .FALSE. 1419 1488 1420 DO lsp = 1, nspec 1489 spc_name_av = TRIM(chem_species(lsp)%name)//'_av' 1421 1422 spc_name_av = TRIM(chem_species(lsp)%name)//'_av' !< for time-averaged chemical conc. 1490 1423 IF (TRIM( field_char ) == TRIM(chem_species(lsp)%name) ) THEN 1491 IF ( k == 1 ) READ ( 13 ) tmp_3d !< read data into tmp_3d 1492 chem_species(lsp)%conc = tmp_3d !< fill ..%conc in the restart run 1493 ELSEIF (TRIM( field_char ) == spc_name_av ) THEN 1494 IF( k == 1 ) READ ( 13 ) tmp_3d 1495 chem_species(lsp)%conc_av = tmp_3d 1424 IF ( k == 1 ) READ ( 13 ) tmp_3d 1425 chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 1426 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 1427 chem_found = .TRUE. 1428 ELSEIF (TRIM( field_char ) == spc_name_av ) THEN 1429 IF ( k == 1 ) READ ( 13 ) tmp_3d 1430 chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 1431 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 1432 chem_found = .TRUE. 1496 1433 ENDIF 1497 1434 1498 1435 ENDDO 1499 1436 IF ( .NOT. chem_found ) THEN 1437 WRITE( message_string, * ) 'unknown variable named "', & 1438 TRIM( field_char ), '" found in', & 1439 '&data from prior run on PE ', myid 1440 CALL message( 'chem_read_restart_data', 'CM0008', 1, 2, 0, 6, 0 ) 1441 ENDIF 1500 1442 ENDDO 1501 1443 … … 1511 1453 1512 1454 END SUBROUTINE chem_read_restart_data 1513 ! 1514 !------------------------------------------------------------------------------! 1515 ! 1516 ! Description: 1517 ! ------------ 1518 !> Subroutine calculating tendencies for chemical species 1519 !------------------------------------------------------------------------------! 1520 1521 SUBROUTINE chem_tendency ( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs, & 1455 1456 1457 !------------------------------------------------------------------------------! 1458 ! 1459 ! Description: 1460 ! ------------ 1461 !> Subroutine calculating prognostic equations for chemical species 1462 !> (cache-optimized). 1463 !> Routine is called separately for each chemical species over a loop from 1464 !> prognostic_equations. 1465 !------------------------------------------------------------------------------! 1466 SUBROUTINE chem_prognostic_equations_ij ( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs, & 1522 1467 i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs, & 1523 1468 flux_l_cs, diss_l_cs ) … … 1530 1475 USE surface_mod, ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & 1531 1476 surf_usm_v 1532 !> Only one chem_spcs recieved prog_eqns at time. 1533 !> cem_tendency is called in prog_eqns over loop => nvar 1477 1534 1478 1535 1479 IMPLICIT NONE … … 1565 1509 ! 1566 1510 1567 !-- Diffusion terms ( ie hinteren 3 sind 0)1511 !-- Diffusion terms (the last three arguments are zero) 1568 1512 1569 1513 CALL diffusion_s( i, j, cs_scalar, & … … 1609 1553 ENDIF 1610 1554 ENDIF 1611 ! 1612 END SUBROUTINE chem_tendency 1555 1556 END SUBROUTINE chem_prognostic_equations_ij 1557 1558 1559 !------------------------------------------------------------------------------! 1560 ! 1561 ! Description: 1562 ! ------------ 1563 !> Subroutine calculating prognostic equations for chemical species 1564 !> (vector-optimized). 1565 !> Routine is called separately for each chemical species over a loop from 1566 !> prognostic_equations. 1567 !------------------------------------------------------------------------------! 1568 SUBROUTINE chem_prognostic_equations ( cs_scalar_p, cs_scalar, tcs_scalar_m, & 1569 pr_init_cs, ilsp ) 1570 1571 USE advec_s_pw_mod, & 1572 ONLY: advec_s_pw 1573 1574 USE advec_s_up_mod, & 1575 ONLY: advec_s_up 1576 1577 USE advec_ws, & 1578 ONLY: advec_s_ws 1579 1580 USE diffusion_s_mod, & 1581 ONLY: diffusion_s 1582 1583 USE indices, & 1584 ONLY: nxl, nxr, nyn, nys, wall_flags_0 1585 1586 USE pegrid 1587 1588 USE surface_mod, & 1589 ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & 1590 surf_usm_v 1591 1592 IMPLICIT NONE 1593 1594 INTEGER :: i !< running index 1595 INTEGER :: j !< running index 1596 INTEGER :: k !< running index 1597 1598 INTEGER(iwp),INTENT(IN) :: ilsp !< 1599 1600 REAL(wp), DIMENSION(0:nz+1) :: pr_init_cs !< 1601 1602 REAL(wp), DIMENSION(:,:,:), POINTER :: cs_scalar !< 1603 REAL(wp), DIMENSION(:,:,:), POINTER :: cs_scalar_p !< 1604 REAL(wp), DIMENSION(:,:,:), POINTER :: tcs_scalar_m !< 1605 1606 1607 ! 1608 !-- Tendency terms for chemical species 1609 tend = 0.0_wp 1610 ! 1611 !-- Advection terms 1612 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1613 IF ( ws_scheme_sca ) THEN 1614 CALL advec_s_ws( cs_scalar, 'kc' ) 1615 ELSE 1616 CALL advec_s_pw( cs_scalar ) 1617 ENDIF 1618 ELSE 1619 CALL advec_s_up( cs_scalar ) 1620 ENDIF 1621 ! 1622 !-- Diffusion terms (the last three arguments are zero) 1623 CALL diffusion_s( cs_scalar, & 1624 surf_def_h(0)%cssws(ilsp,:), & 1625 surf_def_h(1)%cssws(ilsp,:), & 1626 surf_def_h(2)%cssws(ilsp,:), & 1627 surf_lsm_h%cssws(ilsp,:), & 1628 surf_usm_h%cssws(ilsp,:), & 1629 surf_def_v(0)%cssws(ilsp,:), & 1630 surf_def_v(1)%cssws(ilsp,:), & 1631 surf_def_v(2)%cssws(ilsp,:), & 1632 surf_def_v(3)%cssws(ilsp,:), & 1633 surf_lsm_v(0)%cssws(ilsp,:), & 1634 surf_lsm_v(1)%cssws(ilsp,:), & 1635 surf_lsm_v(2)%cssws(ilsp,:), & 1636 surf_lsm_v(3)%cssws(ilsp,:), & 1637 surf_usm_v(0)%cssws(ilsp,:), & 1638 surf_usm_v(1)%cssws(ilsp,:), & 1639 surf_usm_v(2)%cssws(ilsp,:), & 1640 surf_usm_v(3)%cssws(ilsp,:) ) 1641 ! 1642 !-- Prognostic equation for chemical species 1643 DO i = nxl, nxr 1644 DO j = nys, nyn 1645 DO k = nzb+1, nzt 1646 cs_scalar_p(k,j,i) = cs_scalar(k,j,i) & 1647 + ( dt_3d * & 1648 ( tsc(2) * tend(k,j,i) & 1649 + tsc(3) * tcs_scalar_m(k,j,i) & 1650 ) & 1651 - tsc(5) * rdf_sc(k) & 1652 * ( cs_scalar(k,j,i) - pr_init_cs(k) ) & 1653 ) & 1654 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 1655 1656 IF ( cs_scalar_p(k,j,i) < 0.0_wp ) cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i) 1657 ENDDO 1658 ENDDO 1659 ENDDO 1660 ! 1661 !-- Calculate tendencies for the next Runge-Kutta step 1662 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1663 IF ( intermediate_timestep_count == 1 ) THEN 1664 DO i = nxl, nxr 1665 DO j = nys, nyn 1666 DO k = nzb+1, nzt 1667 tcs_scalar_m(k,j,i) = tend(k,j,i) 1668 ENDDO 1669 ENDDO 1670 ENDDO 1671 ELSEIF ( intermediate_timestep_count < & 1672 intermediate_timestep_count_max ) THEN 1673 DO i = nxl, nxr 1674 DO j = nys, nyn 1675 DO k = nzb+1, nzt 1676 tcs_scalar_m(k,j,i) = - 9.5625_wp * tend(k,j,i) & 1677 + 5.3125_wp * tcs_scalar_m(k,j,i) 1678 ENDDO 1679 ENDDO 1680 ENDDO 1681 ENDIF 1682 ENDIF 1683 1684 END SUBROUTINE chem_prognostic_equations 1685 1613 1686 1614 1687 !------------------------------------------------------------------------------! -
palm/trunk/SOURCE/prognostic_equations.f90
r2766 r2815 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Rename chem_tendency to chem_prognostic_equations, 28 ! implement vector version for air chemistry 29 ! 30 ! 2766 2018-01-22 17:17:47Z kanani 27 31 ! Removed preprocessor directive __chem 28 32 ! … … 270 274 271 275 USE chemistry_model_mod, & 272 ONLY: chem_integrate, chem_ species,&273 chem_ tendency, nspec, nvar, spc_names276 ONLY: chem_integrate, chem_prognostic_equations, & 277 chem_species, nspec, nvar, spc_names 274 278 275 279 USE chem_photolysis_mod, & … … 1301 1305 ! 1302 1306 !-- Loop over chemical species 1303 DO lsp = 1, nvar1304 CALL chem_ tendency ( chem_species(lsp)%conc_p,&1307 DO lsp = 1, nvar 1308 CALL chem_prognostic_equations ( chem_species(lsp)%conc_p, & 1305 1309 chem_species(lsp)%conc, & 1306 1310 chem_species(lsp)%tconc_m, & … … 1340 1344 IMPLICIT NONE 1341 1345 1342 INTEGER(iwp) :: i !< 1343 INTEGER(iwp) :: j !< 1344 INTEGER(iwp) :: k !< 1346 INTEGER(iwp) :: i !< 1347 INTEGER(iwp) :: j !< 1348 INTEGER(iwp) :: k !< 1349 INTEGER(iwp) :: lsp !< running index for chemical species 1345 1350 1346 1351 REAL(wp) :: sbt !< … … 2436 2441 CALL tcm_prognostic() 2437 2442 2443 ! 2444 !-- If required, compute prognostic equation for chemical quantites 2445 IF ( air_chemistry ) THEN 2446 CALL cpu_log( log_point(83), '(chem advec+diff+prog)', 'start' ) 2447 ! 2448 !-- Loop over chemical species 2449 DO lsp = 1, nvar 2450 CALL chem_prognostic_equations ( chem_species(lsp)%conc_p, & 2451 chem_species(lsp)%conc, & 2452 chem_species(lsp)%tconc_m, & 2453 chem_species(lsp)%conc_pr_init, & 2454 lsp ) 2455 ENDDO 2456 2457 CALL cpu_log( log_point(83), '(chem advec+diff+prog)', 'stop' ) 2458 ENDIF ! Chemicals equations 2459 2438 2460 2439 2461 END SUBROUTINE prognostic_equations_vector
Note: See TracChangeset
for help on using the changeset viewer.
