Changeset 3298
- Timestamp:
- Oct 2, 2018 12:21:11 PM (6 years ago)
- Location:
- palm/trunk
- Files:
-
- 2 deleted
- 40 edited
- 13 copied
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE
- Property svn:mergeinfo changed
-
palm/trunk/SOURCE/Makefile
- Property svn:mergeinfo changed
r3294 r3298 25 25 # ----------------- 26 26 # $Id$ 27 # Added missing dependencies and replaced blanks with tabs (forkel) 28 # Added chemistry emission module (Russo) 29 # Added chemistry_model_mod to flow_statistics (basit) 30 # 31 # 3294 2018-10-01 02:37:10Z raasch 27 32 # changes related to modularization of the ocean mode, 28 33 # bugfix: dependency to advec_ws was missed in chemistry_model_mod … … 473 478 check_open.f90 \ 474 479 check_parameters.f90 \ 480 chem_emissions_mod.f90 \ 475 481 chem_gasphase_mod.f90 \ 476 482 chemistry_model_mod.f90 \ … … 576 582 random_gauss.f90 \ 577 583 random_generator_parallel_mod.f90 \ 578 584 read_restart_data_mod.f90 \ 579 585 run_control.f90 \ 580 586 set_slicer_attributes_dvrp.f90 \ … … 636 642 wind_turbine_model_mod.f90 \ 637 643 wrd_write_string.f90 \ 638 644 write_restart_data_mod.f90 639 645 640 646 … … 757 763 bulk_cloud_model_mod.o \ 758 764 chemistry_model_mod.o \ 765 chem_emissions_mod.o \ 759 766 gust_mod.o \ 760 767 init_vertical_profiles.o \ … … 778 785 vertical_nesting_mod.o \ 779 786 wind_turbine_model_mod.o 787 chem_emissions_mod.o: \ 788 chem_gasphase_mod.o \ 789 chem_modules.o \ 790 date_and_time_mod.o \ 791 mod_kinds.o \ 792 modules.o \ 793 netcdf_data_input_mod.o \ 794 surface_mod.o 780 795 chemistry_model_mod.o: \ 781 796 advec_ws.o \ 782 797 chem_gasphase_mod.o \ 783 798 chem_modules.o \ 799 date_and_time_mod.o \ 784 800 diffusion_s.o \ 785 801 mod_kinds.o \ 786 802 modules.o \ 787 803 netcdf_data_input_mod.o \ 788 surface_mod.o 804 surface_mod.o \ 805 user_module.o 789 806 chem_gasphase_mod.o: \ 790 807 mod_kinds.o \ … … 831 848 basic_constants_and_equations_mod.o \ 832 849 bulk_cloud_model_mod.o \ 850 chemistry_model_mod.o \ 833 851 cpulog_mod.o \ 834 852 mod_kinds.o \ … … 868 886 basic_constants_and_equations_mod.o \ 869 887 bulk_cloud_model_mod.o \ 888 chemistry_model_mod.o\ 870 889 cpulog_mod.o \ 871 890 gust_mod.o \ … … 942 961 basic_constants_and_equations_mod.o \ 943 962 bulk_cloud_model_mod.o \ 963 chemistry_model_mod.o \ 944 964 cpulog_mod.o \ 945 965 gust_mod.o \ … … 960 980 basic_constants_and_equations_mod.o \ 961 981 bulk_cloud_model_mod.o \ 982 chemistry_model_mod.o \ 962 983 cpulog_mod.o \ 963 984 date_and_time_mod.o \ … … 986 1007 basic_constants_and_equations_mod.o \ 987 1008 bulk_cloud_model_mod.o \ 988 chem istry_model_mod.o \1009 chem_emissions_mod.o \ 989 1010 cpulog_mod.o \ 990 1011 disturb_heatflux.o \ … … 1217 1238 basic_constants_and_equations_mod.o \ 1218 1239 mod_kinds.o \ 1219 modules.o 1240 modules.o \ 1241 time_to_string.o 1220 1242 modules.o: \ 1221 1243 mod_kinds.o … … 1233 1255 mod_kinds.o 1234 1256 netcdf_data_input_mod.o: \ 1257 chem_modules.o \ 1235 1258 cpulog_mod.o \ 1236 1259 mod_kinds.o \ … … 1274 1297 chem_photolysis_mod.o \ 1275 1298 cpulog_mod.o \ 1299 date_and_time_mod.o \ 1276 1300 land_surface_model_mod.o \ 1277 1301 mod_kinds.o \ … … 1282 1306 pmc_particle_interface.o \ 1283 1307 surface_layer_fluxes_mod.o \ 1308 time_to_string.o \ 1284 1309 write_restart_data_mod.o 1285 1310 parin.o: \ … … 1545 1570 buoyancy.o \ 1546 1571 calc_mean_profile.o \ 1572 chem_emissions_mod.o \ 1547 1573 chemistry_model_mod.o \ 1548 1574 chem_modules.o \ … … 1558 1584 modules.o \ 1559 1585 multi_agent_system_mod.o \ 1560 1586 ocean_mod.o \ 1561 1587 pmc_interface_mod.o \ 1562 1588 prognostic_equations.o \ … … 1567 1593 surface_mod.o \ 1568 1594 synthetic_turbulence_generator_mod.o \ 1595 time_to_string.o \ 1569 1596 turbulence_closure_mod.o\ 1570 1597 urban_surface_mod.o \ … … 1584 1611 radiation_model_mod.o \ 1585 1612 surface_layer_fluxes_mod.o \ 1613 time_to_string.o \ 1586 1614 urban_surface_mod.o 1587 1615 time_to_string.o: \ -
palm/trunk/SOURCE/advec_ws.f90
r3294 r3298 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Add formatted note from ketelsen about statistics for chemistry 28 ! 29 ! 3294 2018-10-01 02:37:10Z raasch 27 30 ! ocean renamed ocean_mode 28 31 ! … … 1697 1700 ) * weight_substep(intermediate_timestep_count) 1698 1701 ENDDO 1702 1703 ! CASE ( 'kc' ) 1704 !kk Has to be implemented for kpp chemistry 1705 1699 1706 1700 1707 END SELECT -
palm/trunk/SOURCE/check_parameters.f90
- Property svn:mergeinfo changed
r3294 r3298 25 25 ! ----------------- 26 26 ! $Id: check_parameters.f90 2520 2017-10-05 13:50:26Z gronemeier & 27 ! - Minor formatting and remarks 28 ! - Call of chem_emissions_check_parameters commented 29 ! (preliminary, as there is still a bug in this subroutine) (forkel) 30 ! - Corrected call for chemistry profile output added (basit) 31 ! - Call for init_vertical_profiles for chemistry removed (basit) 32 ! - Call for chem_init_profiles removed (basit) 33 ! - Setting of initial profiles for chemstry (constant and call of 34 ! init_vertical_profiles) commented (forkel) 35 ! 36 ! 2520 2017-10-05 13:50:26Z gronemeier 27 37 ! changes concerning modularization of ocean option, 28 38 ! init_vertical_profiles moved to separate file to avoid circular dependency … … 679 689 bcm_check_data_output_pr 680 690 691 USE chem_emissions_mod, & 692 ONLY: chem_emissions_check_parameters 693 694 USE chem_modules 695 681 696 USE chemistry_model_mod, & 682 697 ONLY: chem_boundary_conds, chem_check_data_output, & 683 chem_check_data_output_pr, chem_species 684 685 USE chem_modules 698 chem_check_data_output_pr, chem_check_parameters, chem_species 686 699 687 700 USE control_parameters … … 1443 1456 ENDIF 1444 1457 1445 ! 1458 ! Check for chem_emission_mod parameters setting 1459 ! IF ( air_chemistry ) CALL chem_emissions_check_parameters ! forkel preliminary 1460 1461 1462 ! Check for chemitry_model_mod parameters setting 1463 IF ( air_chemistry ) CALL chem_check_parameters 1464 1465 1446 1466 !-- Check the module settings 1447 1467 IF ( ocean_mode ) CALL ocean_check_parameters … … 1468 1488 IF ( humidity ) q_init = q_surface 1469 1489 IF ( passive_scalar ) s_init = s_surface 1470 IF ( air_chemistry ) THEN 1471 DO lsp = 1, nvar 1472 chem_species(lsp)%conc_pr_init = cs_surface(lsp) 1473 ENDDO 1474 ENDIF 1490 ! 1491 !-- TODO 1492 !-- Russo: Is done in chem_init and would overwrite what is done there 1493 !-- --> kanani: Revise 1494 ! IF ( air_chemistry ) THEN 1495 ! DO lsp = 1, nvar 1496 ! chem_species(lsp)%conc_pr_init = cs_surface(lsp) 1497 ! ENDDO 1498 ! ENDIF 1475 1499 ! 1476 1500 !-- … … 1683 1707 ENDIF 1684 1708 ! 1709 !-- TODO 1685 1710 !-- Compute initial chemistry profile using the given chemical species gradients 1686 IF ( air_chemistry ) THEN 1687 DO lsp = 1, nvar 1688 CALL init_vertical_profiles( cs_vertical_gradient_level_ind(lsp,:),& 1689 cs_vertical_gradient_level(lsp,:), & 1690 cs_vertical_gradient(lsp,:), & 1691 chem_species(lsp)%conc_pr_init, & 1692 cs_surface(lsp), bc_cs_t_val(lsp) ) 1693 ENDDO 1694 ENDIF 1711 !-- Russo: Is done in chem_init --> kanani: Revise 1695 1712 1696 1713 ENDIF -
palm/trunk/SOURCE/chem_gasphase_mod.f90
r2766 r3298 1 1 MODULE chem_gasphase_mod 2 3 ! Mechanism: passive 4 ! 2 5 !------------------------------------------------------------------------------! 3 6 ! 4 7 ! ******Module chem_gasphase_mod is automatically generated by kpp4palm ****** 5 8 ! 6 ! *********Please do NOT change this Code ********* 9 ! *********Please do NOT change this Code,it will be ovewritten ********* 10 ! 11 !------------------------------------------------------------------------------! 12 ! This file was created by KPP (http://people.cs.vt.edu/asandu/Software/Kpp/) 13 ! and kpp4palm (created by Klaus Ketelsen). kpp4palm is an adapted version 14 ! of KP4 (Jöckel,P.,Kerkweg,A.,Pozzer,A.,Sander,R.,Tost,H.,Riede, 15 ! H.,Baumgaertner,A.,Gromov,S.,and Kern,B.,2010: Development cycle 2 of 16 ! the Modular Earth Submodel System (MESSy2),Geosci. Model Dev.,3,717-752, 17 ! https://doi.org/10.5194/gmd-3-717-2010). KP4 is part of the Modular Earth 18 ! Submodel System (MESSy),which is is available under the GNU General Public 19 ! License (GPL). 20 ! 21 ! KPP is free software; you can redistribute it and/or modify it under the terms 22 ! of the General Public Licence as published by the Free Software Foundation; 23 ! either version 2 of the License,or (at your option) any later version. 24 ! KPP is distributed in the hope that it will be useful,but WITHOUT ANY WARRANTY; 25 ! without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 26 ! PURPOSE. See the GNU General Public Licence for more details. 7 27 ! 8 28 !------------------------------------------------------------------------------! … … 11 31 ! PALM is free software: you can redistribute it and/or modify it under the 12 32 ! terms of the GNU General Public License as published by the Free Software 13 ! Foundation, either version 3 of the License,or (at your option) any later33 ! Foundation,either version 3 of the License,or (at your option) any later 14 34 ! version. 15 35 ! 16 ! PALM is distributed in the hope that it will be useful, 36 ! PALM is distributed in the hope that it will be useful,but WITHOUT ANY 17 37 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 18 38 ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. 19 39 ! 20 40 ! You should have received a copy of the GNU General Public License along with 21 ! PALM. If not, 41 ! PALM. If not,see <http://www.gnu.org/licenses/>. 22 42 ! 23 43 ! Copyright 1997-2018 Leibniz Universitaet Hannover … … 30 50 ! Former revisions: 31 51 ! ----------------- 32 ! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $ 33 ! Removed preprocessor directive __chem 34 ! 35 ! 2460 2017-09-13 14:47:48Z forkel 36 ! 37 ! 38 ! Variables for photolyis added 39 ! 40 ! 41 ! 52 ! $Id: chem_gasphase_mod.f90 3287 2018-09-28 10:19:58Z forkel $ 53 ! forkel June 2018: qvap,fakt added 54 ! forkel June 2018: reset case in Initialize,Integrate,Update_rconst 55 ! 56 ! 57 ! 3287 2018-09-28 10:19:58Z forkel 58 ! 59 ! forkel Sept. 2017: Variables for photolyis added 42 60 ! 43 61 ! … … 52 70 USE kinds, ONLY: dp=>wp 53 71 54 USE pegrid, ONLY: myid, threads_per_task72 USE pegrid, ONLY: myid, threads_per_task 55 73 56 74 IMPLICIT NONE 57 75 PRIVATE 58 !SAVE ! NOTE: OCCURS AGAIN IN AUTOMATICALLY GENERATED CODE...76 !SAVE ! note: occurs again in automatically generated code ... 59 77 60 78 ! PUBLIC :: IERR_NAMES … … 63 81 ! ,REQ_MCFCT,IP_MAX,jname 64 82 65 PUBLIC :: eqn_names, phot_names, spc_names83 PUBLIC :: eqn_names, phot_names, spc_names 66 84 PUBLIC :: nmaxfixsteps 67 PUBLIC :: atol, rtol68 PUBLIC :: nspec, nreact85 PUBLIC :: atol, rtol 86 PUBLIC :: nspec, nreact 69 87 PUBLIC :: temp 88 PUBLIC :: qvap 89 PUBLIC :: fakt 70 90 PUBLIC :: phot 71 91 PUBLIC :: rconst 72 92 PUBLIC :: nvar 73 93 PUBLIC :: nphot 74 75 PUBLIC :: initialize,integrate,update_rconst 94 PUBLIC :: vl_dim ! PUBLIC to ebable other MODULEs to distiguish between scalar and vec 95 96 PUBLIC :: initialize, integrate, update_rconst 76 97 PUBLIC :: chem_gasphase_integrate 77 98 PUBLIC :: initialize_kpp_ctrl … … 82 103 83 104 LOGICAL, PARAMETER :: l_vector = .FALSE. 84 INTEGER, PARAMETER :: i_lu_di = 0105 INTEGER, PARAMETER :: i_lu_di = 2 85 106 INTEGER, PARAMETER :: vl_dim = 1 86 107 INTEGER :: vl 87 108 88 109 INTEGER :: vl_glo 89 INTEGER :: is, ie110 INTEGER :: is, ie 90 111 91 INTEGER, DIMENSION(vl_dim) :: kacc,krej 112 113 INTEGER, DIMENSION(vl_dim) :: kacc, krej 92 114 INTEGER, DIMENSION(vl_dim) :: ierrv 93 LOGICAL :: data_loaded = . false.115 LOGICAL :: data_loaded = .FALSE. 94 116 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 95 117 ! … … 107 129 ! 108 130 ! File : chem_gasphase_mod_Parameters.f90 109 ! Time : Fri Dec 8 11:54:15 2017110 ! Working directory : /home/forkel-r/palmstuff/work/chemistry201 71117/GASPHASE_PREPROC/tmp_kpp4palm131 ! Time : Tue Sep 25 18:34:57 2018 132 ! Working directory : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm 111 133 ! Equation file : chem_gasphase_mod.kpp 112 134 ! Output root filename : chem_gasphase_mod … … 120 142 121 143 ! NSPEC - Number of chemical species 122 INTEGER, PARAMETER :: nspec = 2144 INTEGER, PARAMETER :: nspec = 2 123 145 ! NVAR - Number of Variable species 124 INTEGER, PARAMETER :: nvar = 2146 INTEGER, PARAMETER :: nvar = 2 125 147 ! NVARACT - Number of Active species 126 INTEGER, PARAMETER :: nvaract = 2148 INTEGER, PARAMETER :: nvaract = 2 127 149 ! NFIX - Number of Fixed species 128 INTEGER, PARAMETER :: nfix = 1150 INTEGER, PARAMETER :: nfix = 1 129 151 ! NREACT - Number of reactions 130 INTEGER, PARAMETER :: nreact = 2152 INTEGER, PARAMETER :: nreact = 2 131 153 ! NVARST - Starting of variables in conc. vect. 132 INTEGER, PARAMETER :: nvarst = 1154 INTEGER, PARAMETER :: nvarst = 1 133 155 ! NFIXST - Starting of fixed in conc. vect. 134 INTEGER, PARAMETER :: nfixst = 3156 INTEGER, PARAMETER :: nfixst = 3 135 157 ! NONZERO - Number of nonzero entries in Jacobian 136 INTEGER, PARAMETER :: nonzero = 2158 INTEGER, PARAMETER :: nonzero = 2 137 159 ! LU_NONZERO - Number of nonzero entries in LU factoriz. of Jacobian 138 INTEGER, PARAMETER :: lu_nonzero = 2160 INTEGER, PARAMETER :: lu_nonzero = 2 139 161 ! CNVAR - (NVAR+1) Number of elements in compressed row format 140 INTEGER, PARAMETER :: cnvar = 3162 INTEGER, PARAMETER :: cnvar = 3 141 163 ! CNEQN - (NREACT+1) Number stoicm elements in compressed col format 142 INTEGER, PARAMETER :: cneqn = 3164 INTEGER, PARAMETER :: cneqn = 3 143 165 ! NHESS - Length of Sparse Hessian 144 INTEGER,PARAMETER :: nhess = 1 145 ! NLOOKAT - Number of species to look at 146 INTEGER,PARAMETER :: nlookat = 0 147 ! NMONITOR - Number of species to monitor 148 INTEGER,PARAMETER :: nmonitor = 0 166 INTEGER, PARAMETER :: nhess = 1 149 167 ! NMASS - Number of atoms to check mass balance 150 INTEGER, PARAMETER :: nmass = 1168 INTEGER, PARAMETER :: nmass = 1 151 169 152 170 ! Index declaration for variable species in C and VAR 153 171 ! VAR(ind_spc) = C(ind_spc) 154 172 155 INTEGER, PARAMETER,PUBLIC :: ind_pm10 = 1156 INTEGER, PARAMETER,PUBLIC :: ind_pm25 = 2173 INTEGER, PARAMETER, PUBLIC :: ind_pm10 = 1 174 INTEGER, PARAMETER, PUBLIC :: ind_pm25 = 2 157 175 158 176 ! Index declaration for fixed species in C … … 165 183 166 184 ! NJVRP - Length of sparse Jacobian JVRP 167 INTEGER, PARAMETER :: njvrp = 2185 INTEGER, PARAMETER :: njvrp = 2 168 186 169 187 ! NSTOICM - Length of Sparse Stoichiometric Matrix 170 INTEGER, PARAMETER :: nstoicm = 1188 INTEGER, PARAMETER :: nstoicm = 1 171 189 172 190 … … 186 204 ! 187 205 ! File : chem_gasphase_mod_Global.f90 188 ! Time : Fri Dec 8 11:54:15 2017189 ! Working directory : /home/forkel-r/palmstuff/work/chemistry201 71117/GASPHASE_PREPROC/tmp_kpp4palm206 ! Time : Tue Sep 25 18:34:57 2018 207 ! Working directory : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm 190 208 ! Equation file : chem_gasphase_mod.kpp 191 209 ! Output root filename : chem_gasphase_mod … … 207 225 REAL(kind=dp):: fix(nfix) 208 226 ! VAR,FIX are chunks of array C 209 equivalence( c(1),var(1))227 EQUIVALENCE( c(1), var(1)) 210 228 ! RCONST - Rate constants (global) 211 229 REAL(kind=dp):: rconst(nreact) 212 230 ! TIME - Current integration time 213 231 REAL(kind=dp):: time 214 ! SUN - Sunlight intensity between [0,1]215 REAL(kind=dp):: sun216 232 ! TEMP - Temperature 217 REAL(dp),dimension(:),allocatable :: temp 218 ! RTOLS - (scalar) Relative tolerance 219 REAL(kind=dp):: rtols 233 REAL(kind=dp):: temp 220 234 ! TSTART - Integration start time 221 235 REAL(kind=dp):: tstart 222 ! TEND - Integration end time223 REAL(kind=dp):: tend224 ! DT - Integration step225 REAL(kind=dp):: dt226 236 ! ATOL - Absolute tolerance 227 237 REAL(kind=dp):: atol(nvar) … … 230 240 ! STEPMIN - Lower bound for integration step 231 241 REAL(kind=dp):: stepmin 232 ! STEPMAX - Upper bound for integration step233 REAL(kind=dp):: stepmax234 242 ! CFACTOR - Conversion factor for concentration units 235 243 REAL(kind=dp):: cfactor 236 ! DDMTYPE - DDM sensitivity w.r.t.: 0=init.val.,1=params237 INTEGER :: ddmtype238 244 239 245 ! INLINED global variable declarations 240 246 241 ! declaration of global variable declarations for photolysis will come from 247 ! QVAP - Water vapor 248 REAL(kind=dp):: qvap 249 ! FAKT - Conversion factor 250 REAL(kind=dp):: fakt 251 242 252 243 253 ! INLINED global variable declarations … … 260 270 ! 261 271 ! File : chem_gasphase_mod_JacobianSP.f90 262 ! Time : Fri Dec 8 11:54:15 2017263 ! Working directory : /home/forkel-r/palmstuff/work/chemistry201 71117/GASPHASE_PREPROC/tmp_kpp4palm272 ! Time : Tue Sep 25 18:34:57 2018 273 ! Working directory : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm 264 274 ! Equation file : chem_gasphase_mod.kpp 265 275 ! Output root filename : chem_gasphase_mod … … 275 285 276 286 277 INTEGER, PARAMETER,DIMENSION(2):: lu_irow = (/ &287 INTEGER, PARAMETER, DIMENSION(2):: lu_irow = (/ & 278 288 1, 2 /) 279 289 280 INTEGER, PARAMETER,DIMENSION(2):: lu_icol = (/ &290 INTEGER, PARAMETER, DIMENSION(2):: lu_icol = (/ & 281 291 1, 2 /) 282 292 283 INTEGER, PARAMETER,DIMENSION(3):: lu_crow = (/ &293 INTEGER, PARAMETER, DIMENSION(3):: lu_crow = (/ & 284 294 1, 2, 3 /) 285 295 286 INTEGER, PARAMETER,DIMENSION(3):: lu_diag = (/ &296 INTEGER, PARAMETER, DIMENSION(3):: lu_diag = (/ & 287 297 1, 2, 3 /) 288 298 … … 304 314 ! 305 315 ! File : chem_gasphase_mod_Monitor.f90 306 ! Time : Fri Dec 8 11:54:15 2017307 ! Working directory : /home/forkel-r/palmstuff/work/chemistry201 71117/GASPHASE_PREPROC/tmp_kpp4palm316 ! Time : Tue Sep 25 18:34:57 2018 317 ! Working directory : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm 308 318 ! Equation file : chem_gasphase_mod.kpp 309 319 ! Output root filename : chem_gasphase_mod … … 315 325 316 326 317 CHARACTER(len=15), PARAMETER,DIMENSION(2):: spc_names = (/ &327 CHARACTER(len=15), PARAMETER, DIMENSION(2):: spc_names = (/ & 318 328 'PM10 ','PM25 ' /) 319 329 320 INTEGER,DIMENSION(1):: lookat 321 INTEGER,DIMENSION(1):: monitor 322 CHARACTER(len=15),DIMENSION(1):: smass 323 CHARACTER(len=100),PARAMETER,DIMENSION(2):: eqn_names = (/ & 330 CHARACTER(len=100), PARAMETER, DIMENSION(2):: eqn_names = (/ & 324 331 'PM10 --> PM10 ',& 325 332 'PM25 --> PM25 ' /) … … 329 336 ! inline f90_data: declaration of global variables for photolysis 330 337 ! REAL(kind=dp):: phot(nphot)must eventually be moved to global later for 331 INTEGER, PARAMETER :: nphot = 1338 INTEGER, PARAMETER :: nphot = 1 332 339 ! phot photolysis frequencies 333 340 REAL(kind=dp):: phot(nphot) 334 341 335 INTEGER, PARAMETER,PUBLIC :: j_no2 = 1336 337 CHARACTER(len=15), PARAMETER,DIMENSION(nphot):: phot_names = (/ &342 INTEGER, PARAMETER, PUBLIC :: j_no2 = 1 343 344 CHARACTER(len=15), PARAMETER, DIMENSION(nphot):: phot_names = (/ & 338 345 'J_NO2 '/) 339 346 … … 367 374 ! 368 375 ! File : chem_gasphase_mod_Initialize.f90 369 ! Time : Fri Dec 8 11:54:15 2017370 ! Working directory : /home/forkel-r/palmstuff/work/chemistry201 71117/GASPHASE_PREPROC/tmp_kpp4palm376 ! Time : Tue Sep 25 18:34:57 2018 377 ! Working directory : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm 371 378 ! Equation file : chem_gasphase_mod.kpp 372 379 ! Output root filename : chem_gasphase_mod … … 393 400 ! 394 401 ! File : chem_gasphase_mod_Integrator.f90 395 ! Time : Fri Dec 8 11:54:15 2017396 ! Working directory : /home/forkel-r/palmstuff/work/chemistry201 71117/GASPHASE_PREPROC/tmp_kpp4palm402 ! Time : Tue Sep 25 18:34:57 2018 403 ! Working directory : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm 397 404 ! Equation file : chem_gasphase_mod.kpp 398 405 ! Output root filename : chem_gasphase_mod … … 432 439 433 440 !~~~> statistics on the work performed by the rosenbrock method 434 INTEGER, PARAMETER :: nfun=1,njac=2,nstp=3,nacc=4,&435 nrej=5, ndec=6,nsol=7,nsng=8,&436 ntexit=1, nhexit=2,nhnew = 3441 INTEGER, PARAMETER :: nfun=1, njac=2, nstp=3, nacc=4, & 442 nrej=5, ndec=6, nsol=7, nsng=8, & 443 ntexit=1, nhexit=2, nhnew = 3 437 444 438 445 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ … … 451 458 ! 452 459 ! File : chem_gasphase_mod_LinearAlgebra.f90 453 ! Time : Fri Dec 8 11:54:15 2017454 ! Working directory : /home/forkel-r/palmstuff/work/chemistry201 71117/GASPHASE_PREPROC/tmp_kpp4palm460 ! Time : Tue Sep 25 18:34:57 2018 461 ! Working directory : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm 455 462 ! Equation file : chem_gasphase_mod.kpp 456 463 ! Output root filename : chem_gasphase_mod … … 478 485 ! 479 486 ! File : chem_gasphase_mod_Jacobian.f90 480 ! Time : Fri Dec 8 11:54:15 2017481 ! Working directory : /home/forkel-r/palmstuff/work/chemistry201 71117/GASPHASE_PREPROC/tmp_kpp4palm487 ! Time : Tue Sep 25 18:34:57 2018 488 ! Working directory : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm 482 489 ! Equation file : chem_gasphase_mod.kpp 483 490 ! Output root filename : chem_gasphase_mod … … 505 512 ! 506 513 ! File : chem_gasphase_mod_Function.f90 507 ! Time : Fri Dec 8 11:54:15 2017508 ! Working directory : /home/forkel-r/palmstuff/work/chemistry201 71117/GASPHASE_PREPROC/tmp_kpp4palm514 ! Time : Tue Sep 25 18:34:57 2018 515 ! Working directory : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm 509 516 ! Equation file : chem_gasphase_mod.kpp 510 517 ! Output root filename : chem_gasphase_mod … … 534 541 ! 535 542 ! File : chem_gasphase_mod_Rates.f90 536 ! Time : Fri Dec 8 11:54:15 2017537 ! Working directory : /home/forkel-r/palmstuff/work/chemistry201 71117/GASPHASE_PREPROC/tmp_kpp4palm543 ! Time : Tue Sep 25 18:34:57 2018 544 ! Working directory : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm 538 545 ! Equation file : chem_gasphase_mod.kpp 539 546 ! Output root filename : chem_gasphase_mod … … 560 567 ! 561 568 ! File : chem_gasphase_mod_Util.f90 562 ! Time : Fri Dec 8 11:54:15 2017563 ! Working directory : /home/forkel-r/palmstuff/work/chemistry201 71117/GASPHASE_PREPROC/tmp_kpp4palm569 ! Time : Tue Sep 25 18:34:57 2018 570 ! Working directory : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm 564 571 ! Equation file : chem_gasphase_mod.kpp 565 572 ! Output root filename : chem_gasphase_mod … … 575 582 576 583 ! notes: 577 ! - 578 ! - 579 ! - 580 ! - 581 ! - 582 ! - 583 ! - 584 ! - l_vector is automatically defined by kp4 585 ! - vl_dim is automatically defined by kp4 586 ! - i_lu_di is automatically defined by kp4 587 ! - wanted is automatically defined by xmecca 588 ! - icntrl rcntrl are automatically defined by kpp 589 ! - "USE messy_main_tools" is in MODULE_header of messy_mecca_kpp.f90 590 ! - SAVE will be automatically added by kp4 584 591 585 592 !SAVE … … 587 594 ! for fixed time step control 588 595 ! ... max. number of fixed time steps (sum must be 1) 589 INTEGER, PARAMETER :: nmaxfixsteps = 50596 INTEGER, PARAMETER :: nmaxfixsteps = 50 590 597 ! ... switch for fixed time stepping 591 LOGICAL, PUBLIC :: l_fixed_step = .false.592 INTEGER, PUBLIC :: nfsteps = 1598 LOGICAL, PUBLIC :: l_fixed_step = .FALSE. 599 INTEGER, PUBLIC :: nfsteps = 1 593 600 ! ... number of kpp control PARAMETERs 594 INTEGER, PARAMETER,PUBLIC :: nkppctrl = 20601 INTEGER, PARAMETER, PUBLIC :: nkppctrl = 20 595 602 ! 596 INTEGER, DIMENSION(nkppctrl), PUBLIC :: icntrl = 0597 REAL(dp), DIMENSION(nkppctrl),PUBLIC :: rcntrl = 0.0_dp598 REAL(dp), DIMENSION(nmaxfixsteps),PUBLIC :: t_steps = 0.0_dp603 INTEGER, DIMENSION(nkppctrl), PUBLIC :: icntrl = 0 604 REAL(dp), DIMENSION(nkppctrl), PUBLIC :: rcntrl = 0.0_dp 605 REAL(dp), DIMENSION(nmaxfixsteps), PUBLIC :: t_steps = 0.0_dp 599 606 600 607 ! END header MODULE initialize_kpp_ctrl_template … … 619 626 END INTERFACE kppsolve 620 627 628 INTERFACE jac_sp 629 MODULE PROCEDURE jac_sp 630 END INTERFACE jac_sp 631 632 INTERFACE k_arr 633 MODULE PROCEDURE k_arr 634 END INTERFACE k_arr 635 636 INTERFACE update_rconst 637 MODULE PROCEDURE update_rconst 638 END INTERFACE update_rconst 639 640 INTERFACE arr2 641 MODULE PROCEDURE arr2 642 END INTERFACE arr2 643 644 INTERFACE initialize_kpp_ctrl 645 MODULE PROCEDURE initialize_kpp_ctrl 646 END INTERFACE initialize_kpp_ctrl 647 648 INTERFACE error_output 649 MODULE PROCEDURE error_output 650 END INTERFACE error_output 651 652 INTERFACE wscal 653 MODULE PROCEDURE wscal 654 END INTERFACE wscal 655 656 !INTERFACE not working INTERFACE waxpy 657 !INTERFACE not working MODULE PROCEDURE waxpy 658 !INTERFACE not working END INTERFACE waxpy 659 660 INTERFACE rosenbrock 661 MODULE PROCEDURE rosenbrock 662 END INTERFACE rosenbrock 663 664 INTERFACE funtemplate 665 MODULE PROCEDURE funtemplate 666 END INTERFACE funtemplate 667 668 INTERFACE jactemplate 669 MODULE PROCEDURE jactemplate 670 END INTERFACE jactemplate 671 621 672 INTERFACE kppdecomp 622 673 MODULE PROCEDURE kppdecomp 623 674 END INTERFACE kppdecomp 624 675 625 INTERFACE wlamch626 MODULE PROCEDURE wlamch627 END INTERFACE wlamch628 629 INTERFACE wlamch_add630 MODULE PROCEDURE wlamch_add631 END INTERFACE wlamch_add632 633 INTERFACE jac_sp634 MODULE PROCEDURE jac_sp635 END INTERFACE jac_sp636 637 INTERFACE k_arr638 MODULE PROCEDURE k_arr639 END INTERFACE k_arr640 641 INTERFACE update_rconst642 MODULE PROCEDURE update_rconst643 END INTERFACE update_rconst644 645 INTERFACE arr2646 MODULE PROCEDURE arr2647 END INTERFACE arr2648 649 INTERFACE initialize_kpp_ctrl650 MODULE PROCEDURE initialize_kpp_ctrl651 END INTERFACE initialize_kpp_ctrl652 653 INTERFACE error_output654 MODULE PROCEDURE error_output655 END INTERFACE error_output656 657 !interface not working INTERFACE wcopy658 !interface not working MODULE PROCEDURE wcopy659 !interface not working END INTERFACE wcopy660 661 INTERFACE wscal662 MODULE PROCEDURE wscal663 END INTERFACE wscal664 665 !interface not working INTERFACE waxpy666 !interface not working MODULE PROCEDURE waxpy667 !interface not working END INTERFACE waxpy668 669 INTERFACE rosenbrock670 MODULE PROCEDURE rosenbrock671 END INTERFACE rosenbrock672 673 INTERFACE funtemplate674 MODULE PROCEDURE funtemplate675 END INTERFACE funtemplate676 677 INTERFACE jactemplate678 MODULE PROCEDURE jactemplate679 END INTERFACE jactemplate680 681 676 INTERFACE chem_gasphase_integrate 682 677 MODULE PROCEDURE chem_gasphase_integrate 683 678 END INTERFACE chem_gasphase_integrate 684 679 685 INTERFACE fill_temp686 MODULE PROCEDURE fill_temp687 END INTERFACE fill_temp688 PUBLIC fill_temp689 690 680 691 681 CONTAINS … … 694 684 695 685 686 INTEGER :: j, k 696 687 697 688 INTEGER :: i 698 689 REAL(kind=dp):: x 699 690 k = is 700 691 cfactor = 1.000000e+00_dp 701 692 702 x = (0.)*cfactor 703 DO i = 1,nvar 704 var(i) = x 693 x = (0.) * cfactor 694 DO i = 1 , nvar 705 695 ENDDO 706 696 707 x = (0.) *cfactor708 DO i = 1 ,nfix697 x = (0.) * cfactor 698 DO i = 1 , nfix 709 699 fix(i) = x 710 700 ENDDO … … 720 710 END SUBROUTINE initialize 721 711 722 SUBROUTINE integrate( tin, tout,&723 icntrl_u, rcntrl_u,istatus_u,rstatus_u,ierr_u)724 725 726 REAL(kind=dp), INTENT(in):: tin ! start time727 REAL(kind=dp), INTENT(in):: tout ! END time712 SUBROUTINE integrate( tin, tout, & 713 icntrl_u, rcntrl_u, istatus_u, rstatus_u, ierr_u) 714 715 716 REAL(kind=dp), INTENT(IN):: tin ! start time 717 REAL(kind=dp), INTENT(IN):: tout ! END time 728 718 ! OPTIONAL input PARAMETERs and statistics 729 INTEGER, INTENT( in), OPTIONAL :: icntrl_u(20)730 REAL(kind=dp), INTENT(in), OPTIONAL :: rcntrl_u(20)731 INTEGER, INTENT( out),OPTIONAL :: istatus_u(20)732 REAL(kind=dp), INTENT(out),OPTIONAL :: rstatus_u(20)733 INTEGER, INTENT( out),OPTIONAL :: ierr_u734 735 REAL(kind=dp):: rcntrl(20), rstatus(20)736 INTEGER :: icntrl(20), istatus(20),ierr737 738 INTEGER, SAVE :: ntotal = 0719 INTEGER, INTENT(IN), OPTIONAL :: icntrl_u(20) 720 REAL(kind=dp), INTENT(IN), OPTIONAL :: rcntrl_u(20) 721 INTEGER, INTENT(OUT), OPTIONAL :: istatus_u(20) 722 REAL(kind=dp), INTENT(OUT), OPTIONAL :: rstatus_u(20) 723 INTEGER, INTENT(OUT), OPTIONAL :: ierr_u 724 725 REAL(kind=dp):: rcntrl(20), rstatus(20) 726 INTEGER :: icntrl(20), istatus(20), ierr 727 728 INTEGER, SAVE :: ntotal = 0 739 729 740 730 icntrl(:) = 0 … … 744 734 745 735 !~~~> fine-tune the integrator: 746 icntrl(1) = 0 ! 0 - non- autonomous,1 -autonomous747 icntrl(2) = 0 ! 0 - vector tolerances,1 -scalars748 749 ! IF OPTIONAL PARAMETERs are given, and IF they are >0,736 icntrl(1) = 0 ! 0 - non- autonomous, 1 - autonomous 737 icntrl(2) = 0 ! 0 - vector tolerances, 1 - scalars 738 739 ! IF OPTIONAL PARAMETERs are given, and IF they are >0, 750 740 ! THEN they overwrite default settings. 751 IF ( present(icntrl_u))THEN752 where(icntrl_u(:)> 0)icntrl(:) = icntrl_u(:)741 IF (PRESENT(icntrl_u))THEN 742 WHERE(icntrl_u(:)> 0)icntrl(:) = icntrl_u(:) 753 743 ENDIF 754 IF ( present(rcntrl_u))THEN755 where(rcntrl_u(:)> 0)rcntrl(:) = rcntrl_u(:)744 IF (PRESENT(rcntrl_u))THEN 745 WHERE(rcntrl_u(:)> 0)rcntrl(:) = rcntrl_u(:) 756 746 ENDIF 757 747 758 748 759 CALL rosenbrock(nvar, var,tin,tout, &760 atol, rtol, &761 rcntrl, icntrl,rstatus,istatus,ierr)749 CALL rosenbrock(nvar, var, tin, tout, & 750 atol, rtol, & 751 rcntrl, icntrl, rstatus, istatus, ierr) 762 752 763 753 !~~~> debug option: show no of steps … … 768 758 ! IF OPTIONAL PARAMETERs are given for output they 769 759 ! are updated with the RETURN information 770 IF ( present(istatus_u))istatus_u(:) = istatus(:)771 IF ( present(rstatus_u))rstatus_u(:) = rstatus(:)772 IF ( present(ierr_u))ierr_u = ierr760 IF (PRESENT(istatus_u))istatus_u(:) = istatus(:) 761 IF (PRESENT(rstatus_u))rstatus_u(:) = rstatus(:) 762 IF (PRESENT(ierr_u)) ierr_u = ierr 773 763 774 764 END SUBROUTINE integrate 775 765 776 SUBROUTINE fun(v, f,rct,vdot)766 SUBROUTINE fun(v, f, rct, vdot) 777 767 778 768 ! V - Concentrations of variable species (local) … … 794 784 END SUBROUTINE fun 795 785 796 SUBROUTINE kppsolve(jvs, x)786 SUBROUTINE kppsolve(jvs, x) 797 787 798 788 ! JVS - sparse Jacobian of variables … … 801 791 REAL(kind=dp):: x(nvar) 802 792 803 x(2) = x(2) / jvs(2)804 x(1) = x(1) / jvs(1)793 x(2) = x(2) / jvs(2) 794 x(1) = x(1) / jvs(1) 805 795 806 796 END SUBROUTINE kppsolve 807 797 808 SUBROUTINE kppdecomp( jvs,ier) 809 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 810 ! Sparse LU factorization 811 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 812 813 814 INTEGER :: ier 815 REAL(kind=dp):: jvs(lu_nonzero),w(nvar),a 816 INTEGER :: k,kk,j,jj 817 818 a = 0. ! mz_rs_20050606 819 ier = 0 820 DO k=1,nvar 821 ! mz_rs_20050606: don't check if real value == 0 822 ! IF(jvs( lu_diag(k)).eq. 0.)THEN 823 IF(abs(jvs(lu_diag(k)))< tiny(a))THEN 824 ier = k 825 RETURN 826 ENDIF 827 DO kk = lu_crow(k),lu_crow(k+ 1)- 1 828 w( lu_icol(kk)) = jvs(kk) 829 ENDDO 830 DO kk = lu_crow(k),lu_diag(k)- 1 831 j = lu_icol(kk) 832 a = - w(j)/ jvs( lu_diag(j)) 833 w(j) = - a 834 DO jj = lu_diag(j)+ 1,lu_crow(j+ 1)- 1 835 w( lu_icol(jj)) = w( lu_icol(jj))+ a*jvs(jj) 836 ENDDO 837 ENDDO 838 DO kk = lu_crow(k),lu_crow(k+ 1)- 1 839 jvs(kk) = w( lu_icol(kk)) 840 ENDDO 841 ENDDO 842 843 END SUBROUTINE kppdecomp 844 845 REAL(kind=dp)FUNCTION wlamch( c) 846 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 847 ! returns epsilon machine 848 ! after LAPACK 849 ! replace this by the function from the optimized LAPACK implementation: 850 ! CALL SLAMCH('E') or CALL DLAMCH('E') 851 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 852 ! USE chem_gasphase_mod_Precision 853 854 CHARACTER :: c 855 INTEGER :: i 856 REAL(kind=dp),SAVE :: eps 857 REAL(kind=dp) :: suma 858 REAL(kind=dp),PARAMETER :: one=1.0_dp, half=0.5_dp 859 LOGICAL,SAVE :: first=.true. 860 861 IF (first)THEN 862 first = .false. 863 eps = half**(16) 864 DO i = 17,80 865 eps = eps*half 866 CALL wlamch_add(one,eps,suma) 867 IF (suma.le.one)goto 10 868 ENDDO 869 PRINT*,'ERROR IN WLAMCH. EPS < ',Eps 870 RETURN 871 10 eps = eps*2 872 i = i- 1 873 ENDIF 874 875 wlamch = eps 876 877 END FUNCTION wlamch 878 879 SUBROUTINE wlamch_add( a,b,suma) 880 ! USE chem_gasphase_mod_Precision 881 882 REAL(kind=dp)a,b,suma 883 suma = a + b 884 885 END SUBROUTINE wlamch_add 886 887 SUBROUTINE jac_sp(v,f,rct,jvs) 798 SUBROUTINE jac_sp(v, f, rct, jvs) 888 799 889 800 ! V - Concentrations of variable species (local) … … 914 825 END SUBROUTINE jac_sp 915 826 916 elemental REAL(kind=dp)FUNCTION k_arr (k_298, tdep,temp)827 elemental REAL(kind=dp)FUNCTION k_arr (k_298, tdep, temp) 917 828 ! arrhenius FUNCTION 918 829 919 REAL, INTENT( in):: k_298 ! k at t = 298.15k920 REAL, INTENT( in):: tdep ! temperature dependence921 REAL(kind=dp), INTENT(in):: temp ! temperature830 REAL, INTENT(IN):: k_298 ! k at t = 298.15k 831 REAL, INTENT(IN):: tdep ! temperature dependence 832 REAL(kind=dp), INTENT(IN):: temp ! temperature 922 833 923 834 intrinsic exp 924 835 925 k_arr = k_298 * exp(tdep*(1._dp/temp- 3.3540e-3_dp))! 1/298.15=3.3540e-3836 k_arr = k_298 * exp(tdep* (1._dp/temp- 3.3540e-3_dp))! 1/298.15=3.3540e-3 926 837 927 838 END FUNCTION k_arr 928 839 929 840 SUBROUTINE update_rconst() 930 INTEGER :: j,k841 INTEGER :: k 931 842 932 843 k = is … … 942 853 END SUBROUTINE update_rconst 943 854 944 REAL(kind=dp)FUNCTION arr2( a0,b0,temp) 855 ! END FUNCTION ARR2 856 REAL(kind=dp)FUNCTION arr2( a0, b0, temp) 945 857 REAL(kind=dp):: temp 946 REAL(kind=dp):: a0, b0947 arr2 = a0 * exp( - b0 / temp)858 REAL(kind=dp):: a0, b0 859 arr2 = a0 * exp( - b0 / temp) 948 860 END FUNCTION arr2 949 861 950 SUBROUTINE initialize_kpp_ctrl(status ,iou,modstr)862 SUBROUTINE initialize_kpp_ctrl(status) 951 863 952 864 953 865 ! i/o 954 INTEGER, INTENT(out):: status 955 INTEGER, INTENT(in) :: iou ! LOGICAL i/o unit 956 CHARACTER(len=*),INTENT(in) :: modstr ! read <modstr>.nml 866 INTEGER, INTENT(OUT):: status 957 867 958 868 ! local … … 962 872 ! check fixed time steps 963 873 tsum = 0.0_dp 964 DO i=1, nmaxfixsteps874 DO i=1, nmaxfixsteps 965 875 IF (t_steps(i)< tiny(0.0_dp))exit 966 876 tsum = tsum + t_steps(i) … … 982 892 WRITE(*,*) ' RCNTRL : ',rcntrl 983 893 ! 984 ! note: this is onlymeaningful for vectorized (kp4)rosenbrock- methods894 ! note: this is ONLY meaningful for vectorized (kp4)rosenbrock- methods 985 895 IF (l_vector)THEN 986 896 IF (l_fixed_step)THEN … … 1000 910 END SUBROUTINE initialize_kpp_ctrl 1001 911 1002 SUBROUTINE error_output(c, ierr,pe)1003 1004 1005 INTEGER, INTENT(in):: ierr1006 INTEGER, INTENT(in):: pe1007 REAL(dp), DIMENSION(:),INTENT(in):: c912 SUBROUTINE error_output(c, ierr, pe) 913 914 915 INTEGER, INTENT(IN):: ierr 916 INTEGER, INTENT(IN):: pe 917 REAL(dp), DIMENSION(:), INTENT(IN):: c 1008 918 1009 919 write(6,*) 'ERROR in chem_gasphase_mod ',ierr,C(1) … … 1012 922 END SUBROUTINE error_output 1013 923 1014 SUBROUTINE wcopy(n,x,incx,y,incy) 1015 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1016 ! copies a vector,x,to a vector,y: y <- x 1017 ! only for incX=incY=1 1018 ! after BLAS 1019 ! replace this by the function from the optimized BLAS implementation: 1020 ! CALL SCOPY(N,X,1,Y,1) or CALL DCOPY(N,X,1,Y,1) 1021 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1022 ! USE chem_gasphase_mod_Precision 1023 1024 INTEGER :: i,incx,incy,m,mp1,n 1025 REAL(kind=dp):: x(n),y(n) 1026 1027 IF (n.le.0)RETURN 1028 1029 m = mod(n,8) 1030 IF( m .ne. 0)THEN 1031 DO i = 1,m 1032 y(i) = x(i) 1033 ENDDO 1034 IF( n .lt. 8)RETURN 1035 ENDIF 1036 mp1 = m+ 1 1037 DO i = mp1,n,8 1038 y(i) = x(i) 1039 y(i + 1) = x(i + 1) 1040 y(i + 2) = x(i + 2) 1041 y(i + 3) = x(i + 3) 1042 y(i + 4) = x(i + 4) 1043 y(i + 5) = x(i + 5) 1044 y(i + 6) = x(i + 6) 1045 y(i + 7) = x(i + 7) 1046 ENDDO 1047 1048 END SUBROUTINE wcopy 1049 1050 SUBROUTINE wscal(n,alpha,x,incx) 924 SUBROUTINE wscal(n, alpha, x, incx) 1051 925 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1052 926 ! constant times a vector: x(1:N) <- Alpha*x(1:N) … … 1057 931 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1058 932 1059 INTEGER :: i, incx,m,mp1,n1060 REAL(kind=dp) :: x(n), alpha1061 REAL(kind=dp), PARAMETER :: zero=0.0_dp, one=1.0_dp933 INTEGER :: i, incx, m, mp1, n 934 REAL(kind=dp) :: x(n), alpha 935 REAL(kind=dp), PARAMETER :: zero=0.0_dp, one=1.0_dp 1062 936 1063 937 IF (alpha .eq. one)RETURN 1064 938 IF (n .le. 0)RETURN 1065 939 1066 m = mod(n, 5)1067 IF ( m .ne. 0)THEN940 m = mod(n, 5) 941 IF ( m .ne. 0)THEN 1068 942 IF (alpha .eq. (- one))THEN 1069 DO i = 1, m943 DO i = 1, m 1070 944 x(i) = - x(i) 1071 945 ENDDO 1072 946 ELSEIF (alpha .eq. zero)THEN 1073 DO i = 1, m947 DO i = 1, m 1074 948 x(i) = zero 1075 949 ENDDO 1076 950 ELSE 1077 DO i = 1, m1078 x(i) = alpha* x(i)951 DO i = 1, m 952 x(i) = alpha* x(i) 1079 953 ENDDO 1080 954 ENDIF 1081 IF ( n .lt. 5)RETURN955 IF ( n .lt. 5)RETURN 1082 956 ENDIF 1083 957 mp1 = m + 1 1084 958 IF (alpha .eq. (- one))THEN 1085 DO i = mp1, n,51086 x(i) 959 DO i = mp1, n, 5 960 x(i) = - x(i) 1087 961 x(i + 1) = - x(i + 1) 1088 962 x(i + 2) = - x(i + 2) … … 1091 965 ENDDO 1092 966 ELSEIF (alpha .eq. zero)THEN 1093 DO i = mp1, n,51094 x(i) 967 DO i = mp1, n, 5 968 x(i) = zero 1095 969 x(i + 1) = zero 1096 970 x(i + 2) = zero … … 1099 973 ENDDO 1100 974 ELSE 1101 DO i = mp1, n,51102 x(i) = alpha*x(i)1103 x(i + 1) = alpha* x(i + 1)1104 x(i + 2) = alpha* x(i + 2)1105 x(i + 3) = alpha* x(i + 3)1106 x(i + 4) = alpha* x(i + 4)975 DO i = mp1, n, 5 976 x(i) = alpha* x(i) 977 x(i + 1) = alpha* x(i + 1) 978 x(i + 2) = alpha* x(i + 2) 979 x(i + 3) = alpha* x(i + 3) 980 x(i + 4) = alpha* x(i + 4) 1107 981 ENDDO 1108 982 ENDIF … … 1110 984 END SUBROUTINE wscal 1111 985 1112 SUBROUTINE waxpy(n, alpha,x,incx,y,incy)986 SUBROUTINE waxpy(n, alpha, x, incx, y, incy) 1113 987 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1114 988 ! constant times a vector plus a vector: y <- y + Alpha*x … … 1119 993 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1120 994 1121 INTEGER :: i, incx,incy,m,mp1,n1122 REAL(kind=dp):: x(n), y(n),alpha1123 REAL(kind=dp), PARAMETER :: zero = 0.0_dp995 INTEGER :: i, incx, incy, m, mp1, n 996 REAL(kind=dp):: x(n), y(n), alpha 997 REAL(kind=dp), PARAMETER :: zero = 0.0_dp 1124 998 1125 999 IF (alpha .eq. zero)RETURN 1126 1000 IF (n .le. 0)RETURN 1127 1001 1128 m = mod(n, 4)1129 IF ( m .ne. 0)THEN1130 DO i = 1, m1131 y(i) = y(i) + alpha*x(i)1002 m = mod(n, 4) 1003 IF ( m .ne. 0)THEN 1004 DO i = 1, m 1005 y(i) = y(i) + alpha* x(i) 1132 1006 ENDDO 1133 IF ( n .lt. 4)RETURN1007 IF ( n .lt. 4)RETURN 1134 1008 ENDIF 1135 1009 mp1 = m + 1 1136 DO i = mp1, n,41137 y(i) = y(i) + alpha*x(i)1138 y(i + 1) = y(i + 1) + alpha*x(i + 1)1139 y(i + 2) = y(i + 2) + alpha*x(i + 2)1140 y(i + 3) = y(i + 3) + alpha*x(i + 3)1010 DO i = mp1, n, 4 1011 y(i) = y(i) + alpha* x(i) 1012 y(i + 1) = y(i + 1) + alpha* x(i + 1) 1013 y(i + 2) = y(i + 2) + alpha* x(i + 2) 1014 y(i + 3) = y(i + 3) + alpha* x(i + 3) 1141 1015 ENDDO 1142 1016 1143 1017 END SUBROUTINE waxpy 1144 1018 1145 SUBROUTINE rosenbrock(n, y,tstart,tend,&1146 abstol, reltol, &1147 rcntrl, icntrl,rstatus,istatus,ierr)1019 SUBROUTINE rosenbrock(n, y, tstart, tend, & 1020 abstol, reltol, & 1021 rcntrl, icntrl, rstatus, istatus, ierr) 1148 1022 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1149 1023 ! … … 1172 1046 !~~~> input arguments: 1173 1047 ! 1174 !- y(n)= vector of initial conditions (at t=tstart)1175 !- [tstart,tend] = time range of integration1048 !- y(n) = vector of initial conditions (at t=tstart) 1049 !- [tstart, tend] = time range of integration 1176 1050 ! (if Tstart>Tend the integration is performed backwards in time) 1177 !- reltol,abstol = user precribed accuracy1178 !- SUBROUTINE fun( t,y,ydot) = ode FUNCTION,1051 !- reltol, abstol = user precribed accuracy 1052 !- SUBROUTINE fun( t, y, ydot) = ode FUNCTION, 1179 1053 ! returns Ydot = Y' = F(T,Y) 1180 !- SUBROUTINE jac( t,y,jcb) = jacobian of the ode FUNCTION,1054 !- SUBROUTINE jac( t, y, jcb) = jacobian of the ode FUNCTION, 1181 1055 ! returns Jcb = dFun/dY 1182 !- icntrl(1:20)= INTEGER inputs PARAMETERs1183 !- rcntrl(1:20)= REAL inputs PARAMETERs1056 !- icntrl(1:20) = INTEGER inputs PARAMETERs 1057 !- rcntrl(1:20) = REAL inputs PARAMETERs 1184 1058 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1185 1059 ! 1186 1060 !~~~> output arguments: 1187 1061 ! 1188 !- y(n) - > vector of final states (at t- >tEND)1189 !- istatus(1:20)- > INTEGER output PARAMETERs1190 !- rstatus(1:20)- > REAL output PARAMETERs1191 !- 1062 !- y(n) - > vector of final states (at t- >tend) 1063 !- istatus(1:20) - > INTEGER output PARAMETERs 1064 !- rstatus(1:20) - > REAL output PARAMETERs 1065 !- ierr - > job status upon RETURN 1192 1066 ! success (positive value) or 1193 1067 ! failure (negative value) … … 1262 1136 1263 1137 !~~~> arguments 1264 INTEGER, INTENT( in):: n1265 REAL(kind=dp), INTENT(inout):: y(n)1266 REAL(kind=dp), INTENT(in) :: tstart,tend1267 REAL(kind=dp), INTENT(in) :: abstol(n),reltol(n)1268 INTEGER, INTENT( in):: icntrl(20)1269 REAL(kind=dp), INTENT(in):: rcntrl(20)1270 INTEGER, INTENT( inout):: istatus(20)1271 REAL(kind=dp), INTENT(inout):: rstatus(20)1272 INTEGER, INTENT(out):: ierr1273 !~~~> PARAMETERs of the rosenbrock method, up to 6 stages1274 INTEGER :: ros_s, rosmethod1275 INTEGER, PARAMETER :: rs2=1,rs3=2,rs4=3,rd3=4,rd4=5,rg3=61276 REAL(kind=dp):: ros_a(15), ros_c(15),ros_m(6),ros_e(6),&1277 ros_alpha(6), ros_gamma(6),ros_elo1138 INTEGER, INTENT(IN) :: n 1139 REAL(kind=dp), INTENT(INOUT):: y(n) 1140 REAL(kind=dp), INTENT(IN) :: tstart, tend 1141 REAL(kind=dp), INTENT(IN) :: abstol(n), reltol(n) 1142 INTEGER, INTENT(IN) :: icntrl(20) 1143 REAL(kind=dp), INTENT(IN) :: rcntrl(20) 1144 INTEGER, INTENT(INOUT):: istatus(20) 1145 REAL(kind=dp), INTENT(INOUT):: rstatus(20) 1146 INTEGER, INTENT(OUT) :: ierr 1147 !~~~> PARAMETERs of the rosenbrock method, up to 6 stages 1148 INTEGER :: ros_s, rosmethod 1149 INTEGER, PARAMETER :: rs2=1, rs3=2, rs4=3, rd3=4, rd4=5, rg3=6 1150 REAL(kind=dp):: ros_a(15), ros_c(15), ros_m(6), ros_e(6), & 1151 ros_alpha(6), ros_gamma(6), ros_elo 1278 1152 LOGICAL :: ros_newf(6) 1279 1153 CHARACTER(len=12):: ros_name 1280 1154 !~~~> local variables 1281 REAL(kind=dp):: roundoff, facmin,facmax,facrej,facsafe1282 REAL(kind=dp):: hmin, hmax,hstart1155 REAL(kind=dp):: roundoff, facmin, facmax, facrej, facsafe 1156 REAL(kind=dp):: hmin, hmax, hstart 1283 1157 REAL(kind=dp):: texit 1284 INTEGER :: i, uplimtol,max_no_steps1285 LOGICAL :: autonomous, vectortol1158 INTEGER :: i, uplimtol, max_no_steps 1159 LOGICAL :: autonomous, vectortol 1286 1160 !~~~> PARAMETERs 1287 REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp1288 REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp1161 REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp 1162 REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp 1289 1163 1290 1164 !~~~> initialize statistics … … 1298 1172 ! For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N) 1299 1173 IF (icntrl(2) == 0)THEN 1300 vectortol = . true.1174 vectortol = .TRUE. 1301 1175 uplimtol = n 1302 1176 ELSE 1303 vectortol = . false.1177 vectortol = .FALSE. 1304 1178 uplimtol = 1 1305 1179 ENDIF … … 1313 1187 CASE (3) 1314 1188 CALL ros4 1315 CASE (0, 4)1189 CASE (0, 4) 1316 1190 CALL rodas3 1317 1191 CASE (5) … … 1321 1195 CASE default 1322 1196 PRINT *,'Unknown Rosenbrock method: ICNTRL(3) =',ICNTRL(3) 1323 CALL ros_errormsg(- 2, tstart,zero,ierr)1197 CALL ros_errormsg(- 2, tstart, zero, ierr) 1324 1198 RETURN 1325 1199 END select … … 1332 1206 ELSE 1333 1207 PRINT *,'User-selected max no. of steps: ICNTRL(4) =',ICNTRL(4) 1334 CALL ros_errormsg(- 1, tstart,zero,ierr)1208 CALL ros_errormsg(- 1, tstart, zero, ierr) 1335 1209 RETURN 1336 1210 ENDIF 1337 1211 1338 1212 !~~~> unit roundoff (1+ roundoff>1) 1339 Roundoff = WLAMCH('E')1213 roundoff = epsilon(one) 1340 1214 1341 1215 !~~~> lower bound on the step size: (positive value) … … 1346 1220 ELSE 1347 1221 PRINT *,'User-selected Hmin: RCNTRL(1) =',RCNTRL(1) 1348 CALL ros_errormsg(- 3, tstart,zero,ierr)1222 CALL ros_errormsg(- 3, tstart, zero, ierr) 1349 1223 RETURN 1350 1224 ENDIF … … 1353 1227 hmax = abs(tend-tstart) 1354 1228 ELSEIF (rcntrl(2)> zero)THEN 1355 hmax = min(abs(rcntrl(2)), abs(tend-tstart))1229 hmax = min(abs(rcntrl(2)), abs(tend-tstart)) 1356 1230 ELSE 1357 1231 PRINT *,'User-selected Hmax: RCNTRL(2) =',RCNTRL(2) 1358 CALL ros_errormsg(- 3, tstart,zero,ierr)1232 CALL ros_errormsg(- 3, tstart, zero, ierr) 1359 1233 RETURN 1360 1234 ENDIF 1361 1235 !~~~> starting step size: (positive value) 1362 1236 IF (rcntrl(3) == zero)THEN 1363 hstart = max(hmin, deltamin)1237 hstart = max(hmin, deltamin) 1364 1238 ELSEIF (rcntrl(3)> zero)THEN 1365 hstart = min(abs(rcntrl(3)), abs(tend-tstart))1239 hstart = min(abs(rcntrl(3)), abs(tend-tstart)) 1366 1240 ELSE 1367 1241 PRINT *,'User-selected Hstart: RCNTRL(3) =',RCNTRL(3) 1368 CALL ros_errormsg(- 3, tstart,zero,ierr)1242 CALL ros_errormsg(- 3, tstart, zero, ierr) 1369 1243 RETURN 1370 1244 ENDIF … … 1376 1250 ELSE 1377 1251 PRINT *,'User-selected FacMin: RCNTRL(4) =',RCNTRL(4) 1378 CALL ros_errormsg(- 4, tstart,zero,ierr)1252 CALL ros_errormsg(- 4, tstart, zero, ierr) 1379 1253 RETURN 1380 1254 ENDIF … … 1385 1259 ELSE 1386 1260 PRINT *,'User-selected FacMax: RCNTRL(5) =',RCNTRL(5) 1387 CALL ros_errormsg(- 4, tstart,zero,ierr)1261 CALL ros_errormsg(- 4, tstart, zero, ierr) 1388 1262 RETURN 1389 1263 ENDIF … … 1395 1269 ELSE 1396 1270 PRINT *,'User-selected FacRej: RCNTRL(6) =',RCNTRL(6) 1397 CALL ros_errormsg(- 4, tstart,zero,ierr)1271 CALL ros_errormsg(- 4, tstart, zero, ierr) 1398 1272 RETURN 1399 1273 ENDIF … … 1405 1279 ELSE 1406 1280 PRINT *,'User-selected FacSafe: RCNTRL(7) =',RCNTRL(7) 1407 CALL ros_errormsg(- 4, tstart,zero,ierr)1281 CALL ros_errormsg(- 4, tstart, zero, ierr) 1408 1282 RETURN 1409 1283 ENDIF 1410 1284 !~~~> check IF tolerances are reasonable 1411 DO i=1, uplimtol1412 IF ((abstol(i)<= zero).or. (reltol(i)<= 10.0_dp*roundoff)&1285 DO i=1, uplimtol 1286 IF ((abstol(i)<= zero).or. (reltol(i)<= 10.0_dp* roundoff)& 1413 1287 .or. (reltol(i)>= 1.0_dp))THEN 1414 1288 PRINT *,' AbsTol(',i,') = ',AbsTol(i) 1415 1289 PRINT *,' RelTol(',i,') = ',RelTol(i) 1416 CALL ros_errormsg(- 5, tstart,zero,ierr)1290 CALL ros_errormsg(- 5, tstart, zero, ierr) 1417 1291 RETURN 1418 1292 ENDIF … … 1421 1295 1422 1296 !~~~> CALL rosenbrock method 1423 CALL ros_integrator(y, tstart,tend,texit, &1424 abstol, reltol, &1297 CALL ros_integrator(y, tstart, tend, texit, & 1298 abstol, reltol, & 1425 1299 ! Integration parameters 1426 autonomous, vectortol,max_no_steps, &1427 roundoff, hmin,hmax,hstart, &1428 facmin, facmax,facrej,facsafe, &1300 autonomous, vectortol, max_no_steps, & 1301 roundoff, hmin, hmax, hstart, & 1302 facmin, facmax, facrej, facsafe, & 1429 1303 ! Error indicator 1430 1304 ierr) … … 1435 1309 1436 1310 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1437 SUBROUTINE ros_errormsg(code, t,h,ierr)1311 SUBROUTINE ros_errormsg(code, t, h, ierr) 1438 1312 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1439 1313 ! Handles all error messages 1440 1314 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1441 1315 1442 REAL(kind=dp), INTENT(in):: t,h1443 INTEGER, INTENT(in) :: code1444 INTEGER, INTENT(out):: ierr1316 REAL(kind=dp), INTENT(IN):: t, h 1317 INTEGER, INTENT(IN) :: code 1318 INTEGER, INTENT(OUT):: ierr 1445 1319 1446 1320 ierr = code 1447 print * ,&1321 print * , & 1448 1322 'Forced exit from Rosenbrock due to the following error:' 1449 1323 1450 1324 select CASE (code) 1451 CASE (- 1) 1325 CASE (- 1) 1452 1326 PRINT *,'--> Improper value for maximal no of steps' 1453 CASE (- 2) 1327 CASE (- 2) 1454 1328 PRINT *,'--> Selected Rosenbrock method not implemented' 1455 CASE (- 3) 1329 CASE (- 3) 1456 1330 PRINT *,'--> Hmin/Hmax/Hstart must be positive' 1457 CASE (- 4) 1331 CASE (- 4) 1458 1332 PRINT *,'--> FacMin/FacMax/FacRej must be positive' 1459 1333 CASE (- 5) … … 1464 1338 PRINT *,'--> Step size too small: T + 10*H = T',& 1465 1339 ' or H < Roundoff' 1466 CASE (- 8) 1340 CASE (- 8) 1467 1341 PRINT *,'--> Matrix is repeatedly singular' 1468 1342 CASE default … … 1470 1344 END select 1471 1345 1472 print * ,"t=",t,"and h=",h1346 print * , "t=", t, "and h=", h 1473 1347 1474 1348 END SUBROUTINE ros_errormsg 1475 1349 1476 1350 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1477 SUBROUTINE ros_integrator (y, tstart,tend,t, &1478 abstol, reltol, &1351 SUBROUTINE ros_integrator (y, tstart, tend, t, & 1352 abstol, reltol, & 1479 1353 !~~~> integration PARAMETERs 1480 autonomous, vectortol,max_no_steps, &1481 roundoff, hmin,hmax,hstart, &1482 facmin, facmax,facrej,facsafe, &1354 autonomous, vectortol, max_no_steps, & 1355 roundoff, hmin, hmax, hstart, & 1356 facmin, facmax, facrej, facsafe, & 1483 1357 !~~~> error indicator 1484 1358 ierr) … … 1491 1365 1492 1366 !~~~> input: the initial condition at tstart; output: the solution at t 1493 REAL(kind=dp), INTENT(inout):: y(n)1367 REAL(kind=dp), INTENT(INOUT):: y(n) 1494 1368 !~~~> input: integration interval 1495 REAL(kind=dp), INTENT(in):: tstart,tend1496 !~~~> output: time at which the solution is RETURNed (t=t ENDIF success)1497 REAL(kind=dp), INTENT(out):: t1369 REAL(kind=dp), INTENT(IN):: tstart, tend 1370 !~~~> output: time at which the solution is RETURNed (t=tendIF success) 1371 REAL(kind=dp), INTENT(OUT):: t 1498 1372 !~~~> input: tolerances 1499 REAL(kind=dp), INTENT(in):: abstol(n),reltol(n)1373 REAL(kind=dp), INTENT(IN):: abstol(n), reltol(n) 1500 1374 !~~~> input: integration PARAMETERs 1501 LOGICAL, INTENT(in):: autonomous,vectortol1502 REAL(kind=dp), INTENT(in):: hstart,hmin,hmax1503 INTEGER, INTENT(in):: max_no_steps1504 REAL(kind=dp), INTENT(in):: roundoff,facmin,facmax,facrej,facsafe1375 LOGICAL, INTENT(IN):: autonomous, vectortol 1376 REAL(kind=dp), INTENT(IN):: hstart, hmin, hmax 1377 INTEGER, INTENT(IN):: max_no_steps 1378 REAL(kind=dp), INTENT(IN):: roundoff, facmin, facmax, facrej, facsafe 1505 1379 !~~~> output: error indicator 1506 INTEGER, INTENT(out):: ierr1380 INTEGER, INTENT(OUT):: ierr 1507 1381 ! ~~~~ Local variables 1508 REAL(kind=dp):: ynew(n), fcn0(n),fcn(n)1509 REAL(kind=dp):: k(n* ros_s),dfdt(n)1382 REAL(kind=dp):: ynew(n), fcn0(n), fcn(n) 1383 REAL(kind=dp):: k(n* ros_s), dfdt(n) 1510 1384 #ifdef full_algebra 1511 REAL(kind=dp):: jac0(n, n),ghimj(n,n)1385 REAL(kind=dp):: jac0(n, n), ghimj(n, n) 1512 1386 #else 1513 REAL(kind=dp):: jac0(lu_nonzero), ghimj(lu_nonzero)1387 REAL(kind=dp):: jac0(lu_nonzero), ghimj(lu_nonzero) 1514 1388 #endif 1515 REAL(kind=dp):: h, hnew,hc,hg,fac,tau1516 REAL(kind=dp):: err, yerr(n)1517 INTEGER :: pivot(n), direction,ioffset,j,istage1518 LOGICAL :: rejectlasth, rejectmoreh,singular1389 REAL(kind=dp):: h, hnew, hc, hg, fac, tau 1390 REAL(kind=dp):: err, yerr(n) 1391 INTEGER :: pivot(n), direction, ioffset, j, istage 1392 LOGICAL :: rejectlasth, rejectmoreh, singular 1519 1393 !~~~> local PARAMETERs 1520 REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp1521 REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp1394 REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp 1395 REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp 1522 1396 !~~~> locally called FUNCTIONs 1523 1397 ! REAL(kind=dp) WLAMCH … … 1529 1403 t = tstart 1530 1404 rstatus(nhexit) = zero 1531 h = min( max(abs(hmin), abs(hstart)),abs(hmax))1532 IF (abs(h)<= 10.0_dp* roundoff)h = deltamin1533 1534 IF (t END>= tstart)THEN1405 h = min( max(abs(hmin), abs(hstart)), abs(hmax)) 1406 IF (abs(h)<= 10.0_dp* roundoff)h = deltamin 1407 1408 IF (tend >= tstart)THEN 1535 1409 direction = + 1 1536 1410 ELSE 1537 1411 direction = - 1 1538 1412 ENDIF 1539 h = direction* h1540 1541 rejectlasth=. false.1542 rejectmoreh=. false.1413 h = direction* h 1414 1415 rejectlasth=.FALSE. 1416 rejectmoreh=.FALSE. 1543 1417 1544 1418 !~~~> time loop begins below 1545 1419 1546 timeloop: DO WHILE((direction > 0).and.((t- t END)+ roundoff <= zero)&1547 .or. (direction < 0).and.((tend-t) + roundoff <= zero))1548 1549 IF (istatus(nstp)> max_no_steps)THEN ! too many steps1550 CALL ros_errormsg(- 6, t,h,ierr)1420 timeloop: DO WHILE((direction > 0).and.((t- tend) + roundoff <= zero)& 1421 .or. (direction < 0).and.((tend-t) + roundoff <= zero)) 1422 1423 IF (istatus(nstp)> max_no_steps)THEN ! too many steps 1424 CALL ros_errormsg(- 6, t, h, ierr) 1551 1425 RETURN 1552 1426 ENDIF 1553 IF (((t+ 0.1_dp*h) == t).or.(h <= roundoff))THEN ! step size too small1554 CALL ros_errormsg(- 7, t,h,ierr)1427 IF (((t+ 0.1_dp* h) == t).or.(h <= roundoff))THEN ! step size too small 1428 CALL ros_errormsg(- 7, t, h, ierr) 1555 1429 RETURN 1556 1430 ENDIF 1557 1431 1558 1432 !~~~> limit h IF necessary to avoid going beyond tend 1559 h = min(h, abs(tend-t))1433 h = min(h, abs(tend-t)) 1560 1434 1561 1435 !~~~> compute the FUNCTION at current time 1562 CALL funtemplate(t, y,fcn0)1563 istatus(nfun) = istatus(nfun) + 11436 CALL funtemplate(t, y, fcn0) 1437 istatus(nfun) = istatus(nfun) + 1 1564 1438 1565 1439 !~~~> compute the FUNCTION derivative with respect to t 1566 1440 IF (.not.autonomous)THEN 1567 CALL ros_funtimederivative(t, roundoff,y,&1568 fcn0, dfdt)1441 CALL ros_funtimederivative(t, roundoff, y, & 1442 fcn0, dfdt) 1569 1443 ENDIF 1570 1444 1571 1445 !~~~> compute the jacobian at current time 1572 CALL jactemplate(t, y,jac0)1573 istatus(njac) = istatus(njac) + 11446 CALL jactemplate(t, y, jac0) 1447 istatus(njac) = istatus(njac) + 1 1574 1448 1575 1449 !~~~> repeat step calculation until current step accepted 1576 1450 untilaccepted: do 1577 1451 1578 CALL ros_preparematrix(h, direction,ros_gamma(1),&1579 jac0, ghimj,pivot,singular)1452 CALL ros_preparematrix(h, direction, ros_gamma(1), & 1453 jac0, ghimj, pivot, singular) 1580 1454 IF (singular)THEN ! more than 5 consecutive failed decompositions 1581 CALL ros_errormsg(- 8, t,h,ierr)1455 CALL ros_errormsg(- 8, t, h, ierr) 1582 1456 RETURN 1583 1457 ENDIF 1584 1458 1585 1459 !~~~> compute the stages 1586 stage: DO istage = 1, ros_s1460 stage: DO istage = 1, ros_s 1587 1461 1588 1462 ! current istage offset. current istage vector is k(ioffset+ 1:ioffset+ n) 1589 ioffset = n* (istage-1)1463 ioffset = n* (istage-1) 1590 1464 1591 1465 ! for the 1st istage the FUNCTION has been computed previously 1592 IF (istage == 1)THEN1593 !slim: CALL wcopy(n, fcn0,1,fcn,1)1594 1466 IF (istage == 1)THEN 1467 !slim: CALL wcopy(n, fcn0, 1, fcn, 1) 1468 fcn(1:n) = fcn0(1:n) 1595 1469 ! istage>1 and a new FUNCTION evaluation is needed at the current istage 1596 1470 ELSEIF(ros_newf(istage))THEN 1597 !slim: CALL wcopy(n, y,1,ynew,1)1598 1599 DO j = 1, istage-11600 CALL waxpy(n, ros_a((istage-1)*(istage-2)/2+ j),&1601 k(n* (j- 1)+ 1),1,ynew,1)1471 !slim: CALL wcopy(n, y, 1, ynew, 1) 1472 ynew(1:n) = y(1:n) 1473 DO j = 1, istage-1 1474 CALL waxpy(n, ros_a((istage-1) * (istage-2) /2+ j), & 1475 k(n* (j- 1) + 1), 1, ynew, 1) 1602 1476 ENDDO 1603 tau = t + ros_alpha(istage) *direction*h1604 CALL funtemplate(tau, ynew,fcn)1605 istatus(nfun) = istatus(nfun) + 11477 tau = t + ros_alpha(istage) * direction* h 1478 CALL funtemplate(tau, ynew, fcn) 1479 istatus(nfun) = istatus(nfun) + 1 1606 1480 ENDIF ! IF istage == 1 ELSEIF ros_newf(istage) 1607 !slim: CALL wcopy(n, fcn,1,k(ioffset+ 1),1)1481 !slim: CALL wcopy(n, fcn, 1, k(ioffset+ 1), 1) 1608 1482 k(ioffset+ 1:ioffset+ n) = fcn(1:n) 1609 DO j = 1, istage-11610 hc = ros_c((istage-1) *(istage-2)/2+ j)/(direction*h)1611 CALL waxpy(n, hc,k(n*(j- 1)+ 1),1,k(ioffset+ 1),1)1483 DO j = 1, istage-1 1484 hc = ros_c((istage-1) * (istage-2) /2+ j) /(direction* h) 1485 CALL waxpy(n, hc, k(n* (j- 1) + 1), 1, k(ioffset+ 1), 1) 1612 1486 ENDDO 1613 1487 IF ((.not. autonomous).and.(ros_gamma(istage).ne.zero))THEN 1614 hg = direction* h*ros_gamma(istage)1615 CALL waxpy(n, hg,dfdt,1,k(ioffset+ 1),1)1488 hg = direction* h* ros_gamma(istage) 1489 CALL waxpy(n, hg, dfdt, 1, k(ioffset+ 1), 1) 1616 1490 ENDIF 1617 CALL ros_solve(ghimj, pivot,k(ioffset+ 1))1491 CALL ros_solve(ghimj, pivot, k(ioffset+ 1)) 1618 1492 1619 1493 END DO stage … … 1621 1495 1622 1496 !~~~> compute the new solution 1623 !slim: CALL wcopy(n, y,1,ynew,1)1497 !slim: CALL wcopy(n, y, 1, ynew, 1) 1624 1498 ynew(1:n) = y(1:n) 1625 DO j=1, ros_s1626 CALL waxpy(n, ros_m(j),k(n*(j- 1)+ 1),1,ynew,1)1499 DO j=1, ros_s 1500 CALL waxpy(n, ros_m(j), k(n* (j- 1) + 1), 1, ynew, 1) 1627 1501 ENDDO 1628 1502 1629 1503 !~~~> compute the error estimation 1630 !slim: CALL wscal(n, zero,yerr,1)1504 !slim: CALL wscal(n, zero, yerr, 1) 1631 1505 yerr(1:n) = zero 1632 DO j=1, ros_s1633 CALL waxpy(n, ros_e(j),k(n*(j- 1)+ 1),1,yerr,1)1506 DO j=1, ros_s 1507 CALL waxpy(n, ros_e(j), k(n* (j- 1) + 1), 1, yerr, 1) 1634 1508 ENDDO 1635 err = ros_errornorm(y, ynew,yerr,abstol,reltol,vectortol)1509 err = ros_errornorm(y, ynew, yerr, abstol, reltol, vectortol) 1636 1510 1637 1511 !~~~> new step size is bounded by facmin <= hnew/h <= facmax 1638 fac = min(facmax, max(facmin,facsafe/err**(one/ros_elo)))1639 hnew = h* fac1512 fac = min(facmax, max(facmin, facsafe/err** (one/ros_elo))) 1513 hnew = h* fac 1640 1514 1641 1515 !~~~> check the error magnitude and adjust step size 1642 istatus(nstp) = istatus(nstp) + 11643 IF ((err <= one).or.(h <= hmin))THEN !~~~> accept step1644 istatus(nacc) = istatus(nacc) + 11645 !slim: CALL wcopy(n, ynew,1,y,1)1516 istatus(nstp) = istatus(nstp) + 1 1517 IF ((err <= one).or.(h <= hmin))THEN !~~~> accept step 1518 istatus(nacc) = istatus(nacc) + 1 1519 !slim: CALL wcopy(n, ynew, 1, y, 1) 1646 1520 y(1:n) = ynew(1:n) 1647 t = t + direction* h1648 hnew = max(hmin, min(hnew,hmax))1521 t = t + direction* h 1522 hnew = max(hmin, min(hnew, hmax)) 1649 1523 IF (rejectlasth)THEN ! no step size increase after a rejected step 1650 hnew = min(hnew, h)1524 hnew = min(hnew, h) 1651 1525 ENDIF 1652 1526 rstatus(nhexit) = h 1653 1527 rstatus(nhnew) = hnew 1654 1528 rstatus(ntexit) = t 1655 rejectlasth = . false.1656 rejectmoreh = . false.1529 rejectlasth = .FALSE. 1530 rejectmoreh = .FALSE. 1657 1531 h = hnew 1658 1532 exit untilaccepted ! exit the loop: WHILE step not accepted 1659 1533 ELSE !~~~> reject step 1660 1534 IF (rejectmoreh)THEN 1661 hnew = h* facrej1535 hnew = h* facrej 1662 1536 ENDIF 1663 1537 rejectmoreh = rejectlasth 1664 rejectlasth = . true.1538 rejectlasth = .TRUE. 1665 1539 h = hnew 1666 IF (istatus(nacc)>= 1) istatus(nrej) = istatus(nrej) + 11540 IF (istatus(nacc)>= 1) istatus(nrej) = istatus(nrej) + 1 1667 1541 ENDIF ! err <= 1 1668 1542 … … 1678 1552 1679 1553 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1680 REAL(kind=dp)FUNCTION ros_errornorm(y, ynew,yerr,&1681 abstol, reltol,vectortol)1554 REAL(kind=dp)FUNCTION ros_errornorm(y, ynew, yerr, & 1555 abstol, reltol, vectortol) 1682 1556 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1683 1557 !~~~> computes the "scaled norm" of the error vector yerr … … 1685 1559 1686 1560 ! Input arguments 1687 REAL(kind=dp), INTENT(in):: y(n),ynew(n),&1688 yerr(n), abstol(n),reltol(n)1689 LOGICAL, INTENT(in):: vectortol1561 REAL(kind=dp), INTENT(IN):: y(n), ynew(n), & 1562 yerr(n), abstol(n), reltol(n) 1563 LOGICAL, INTENT(IN):: vectortol 1690 1564 ! Local variables 1691 REAL(kind=dp):: err, scale,ymax1565 REAL(kind=dp):: err, scale, ymax 1692 1566 INTEGER :: i 1693 REAL(kind=dp), PARAMETER :: zero = 0.0_dp1567 REAL(kind=dp), PARAMETER :: zero = 0.0_dp 1694 1568 1695 1569 err = zero 1696 DO i=1, n1697 ymax = max(abs(y(i)), abs(ynew(i)))1570 DO i=1, n 1571 ymax = max(abs(y(i)), abs(ynew(i))) 1698 1572 IF (vectortol)THEN 1699 scale = abstol(i) + reltol(i)*ymax1573 scale = abstol(i) + reltol(i) * ymax 1700 1574 ELSE 1701 scale = abstol(1) + reltol(1)*ymax1575 scale = abstol(1) + reltol(1) * ymax 1702 1576 ENDIF 1703 err = err+ (yerr(i)/scale)**21577 err = err+ (yerr(i) /scale) ** 2 1704 1578 ENDDO 1705 1579 err = sqrt(err/n) 1706 1580 1707 ros_errornorm = max(err, 1.0d-10)1581 ros_errornorm = max(err, 1.0d-10) 1708 1582 1709 1583 END FUNCTION ros_errornorm … … 1711 1585 1712 1586 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1713 SUBROUTINE ros_funtimederivative(t, roundoff,y,&1714 fcn0, dfdt)1587 SUBROUTINE ros_funtimederivative(t, roundoff, y, & 1588 fcn0, dfdt) 1715 1589 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1716 1590 !~~~> the time partial derivative of the FUNCTION by finite differences … … 1718 1592 1719 1593 !~~~> input arguments 1720 REAL(kind=dp), INTENT(in):: t,roundoff,y(n),fcn0(n)1594 REAL(kind=dp), INTENT(IN):: t, roundoff, y(n), fcn0(n) 1721 1595 !~~~> output arguments 1722 REAL(kind=dp), INTENT(out):: dfdt(n)1596 REAL(kind=dp), INTENT(OUT):: dfdt(n) 1723 1597 !~~~> local variables 1724 1598 REAL(kind=dp):: delta 1725 REAL(kind=dp), PARAMETER :: one = 1.0_dp, deltamin = 1.0e-6_dp1726 1727 delta = sqrt(roundoff) *max(deltamin,abs(t))1728 CALL funtemplate(t+ delta, y,dfdt)1729 istatus(nfun) = istatus(nfun) + 11730 CALL waxpy(n, (- one),fcn0,1,dfdt,1)1731 CALL wscal(n, (one/delta),dfdt,1)1599 REAL(kind=dp), PARAMETER :: one = 1.0_dp, deltamin = 1.0e-6_dp 1600 1601 delta = sqrt(roundoff) * max(deltamin, abs(t)) 1602 CALL funtemplate(t+ delta, y, dfdt) 1603 istatus(nfun) = istatus(nfun) + 1 1604 CALL waxpy(n, (- one), fcn0, 1, dfdt, 1) 1605 CALL wscal(n, (one/delta), dfdt, 1) 1732 1606 1733 1607 END SUBROUTINE ros_funtimederivative … … 1735 1609 1736 1610 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1737 SUBROUTINE ros_preparematrix(h, direction,gam,&1738 jac0, ghimj,pivot,singular)1611 SUBROUTINE ros_preparematrix(h, direction, gam, & 1612 jac0, ghimj, pivot, singular) 1739 1613 ! --- --- --- --- --- --- --- --- --- --- --- --- --- 1740 1614 ! Prepares the LHS matrix for stage calculations … … 1748 1622 !~~~> input arguments 1749 1623 #ifdef full_algebra 1750 REAL(kind=dp), INTENT(in):: jac0(n,n)1624 REAL(kind=dp), INTENT(IN):: jac0(n, n) 1751 1625 #else 1752 REAL(kind=dp), INTENT(in):: jac0(lu_nonzero)1626 REAL(kind=dp), INTENT(IN):: jac0(lu_nonzero) 1753 1627 #endif 1754 REAL(kind=dp), INTENT(in):: gam1755 INTEGER, INTENT(in):: direction1628 REAL(kind=dp), INTENT(IN):: gam 1629 INTEGER, INTENT(IN):: direction 1756 1630 !~~~> output arguments 1757 1631 #ifdef full_algebra 1758 REAL(kind=dp), INTENT(out):: ghimj(n,n)1632 REAL(kind=dp), INTENT(OUT):: ghimj(n, n) 1759 1633 #else 1760 REAL(kind=dp), INTENT(out):: ghimj(lu_nonzero)1634 REAL(kind=dp), INTENT(OUT):: ghimj(lu_nonzero) 1761 1635 #endif 1762 LOGICAL, INTENT(out):: singular1763 INTEGER, INTENT(out):: pivot(n)1636 LOGICAL, INTENT(OUT):: singular 1637 INTEGER, INTENT(OUT):: pivot(n) 1764 1638 !~~~> inout arguments 1765 REAL(kind=dp), INTENT(inout):: h ! step size is decreased when lu fails1639 REAL(kind=dp), INTENT(INOUT):: h ! step size is decreased when lu fails 1766 1640 !~~~> local variables 1767 INTEGER :: i, ising,nconsecutive1641 INTEGER :: i, ising, nconsecutive 1768 1642 REAL(kind=dp):: ghinv 1769 REAL(kind=dp), PARAMETER :: one = 1.0_dp, half = 0.5_dp1643 REAL(kind=dp), PARAMETER :: one = 1.0_dp, half = 0.5_dp 1770 1644 1771 1645 nconsecutive = 0 1772 singular = . true.1646 singular = .TRUE. 1773 1647 1774 1648 DO WHILE (singular) 1775 1649 1776 !~~~> construct ghimj = 1/(h* gam)-jac01650 !~~~> construct ghimj = 1/(h* gam) - jac0 1777 1651 #ifdef full_algebra 1778 !slim: CALL wcopy(n* n,jac0,1,ghimj,1)1779 !slim: CALL wscal(n* n,(- one),ghimj,1)1652 !slim: CALL wcopy(n* n, jac0, 1, ghimj, 1) 1653 !slim: CALL wscal(n* n, (- one), ghimj, 1) 1780 1654 ghimj = - jac0 1781 ghinv = one/(direction* h*gam)1782 DO i=1, n1783 ghimj(i, i) = ghimj(i,i)+ ghinv1655 ghinv = one/(direction* h* gam) 1656 DO i=1, n 1657 ghimj(i, i) = ghimj(i, i) + ghinv 1784 1658 ENDDO 1785 1659 #else 1786 !slim: CALL wcopy(lu_nonzero, jac0,1,ghimj,1)1787 !slim: CALL wscal(lu_nonzero, (- one),ghimj,1)1660 !slim: CALL wcopy(lu_nonzero, jac0, 1, ghimj, 1) 1661 !slim: CALL wscal(lu_nonzero, (- one), ghimj, 1) 1788 1662 ghimj(1:lu_nonzero) = - jac0(1:lu_nonzero) 1789 ghinv = one/(direction* h*gam)1790 DO i=1, n1791 ghimj(lu_diag(i)) = ghimj(lu_diag(i)) + ghinv1663 ghinv = one/(direction* h* gam) 1664 DO i=1, n 1665 ghimj(lu_diag(i)) = ghimj(lu_diag(i)) + ghinv 1792 1666 ENDDO 1793 1667 #endif 1794 1668 !~~~> compute lu decomposition 1795 CALL ros_decomp( ghimj, pivot,ising)1669 CALL ros_decomp( ghimj, pivot, ising) 1796 1670 IF (ising == 0)THEN 1797 1671 !~~~> IF successful done 1798 singular = . false.1672 singular = .FALSE. 1799 1673 ELSE ! ising .ne. 0 1800 1674 !~~~> IF unsuccessful half the step size; IF 5 consecutive fails THEN RETURN 1801 istatus(nsng) = istatus(nsng) + 11675 istatus(nsng) = istatus(nsng) + 1 1802 1676 nconsecutive = nconsecutive+1 1803 singular = . true.1677 singular = .TRUE. 1804 1678 PRINT*,'Warning: LU Decomposition returned ISING = ',ISING 1805 1679 IF (nconsecutive <= 5)THEN ! less than 5 consecutive failed decompositions 1806 h = h* half1680 h = h* half 1807 1681 ELSE ! more than 5 consecutive failed decompositions 1808 1682 RETURN … … 1816 1690 1817 1691 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1818 SUBROUTINE ros_decomp( a, pivot,ising)1692 SUBROUTINE ros_decomp( a, pivot, ising) 1819 1693 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1820 1694 ! Template for the LU decomposition … … 1822 1696 !~~~> inout variables 1823 1697 #ifdef full_algebra 1824 REAL(kind=dp), INTENT(inout):: a(n,n)1698 REAL(kind=dp), INTENT(INOUT):: a(n, n) 1825 1699 #else 1826 REAL(kind=dp), INTENT(inout):: a(lu_nonzero)1700 REAL(kind=dp), INTENT(INOUT):: a(lu_nonzero) 1827 1701 #endif 1828 1702 !~~~> output variables 1829 INTEGER, INTENT(out):: pivot(n),ising1703 INTEGER, INTENT(OUT):: pivot(n), ising 1830 1704 1831 1705 #ifdef full_algebra 1832 CALL dgetrf( n, n,a,n,pivot,ising)1706 CALL dgetrf( n, n, a, n, pivot, ising) 1833 1707 #else 1834 CALL kppdecomp(a, ising)1708 CALL kppdecomp(a, ising) 1835 1709 pivot(1) = 1 1836 1710 #endif 1837 istatus(ndec) = istatus(ndec) + 11711 istatus(ndec) = istatus(ndec) + 1 1838 1712 1839 1713 END SUBROUTINE ros_decomp … … 1841 1715 1842 1716 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1843 SUBROUTINE ros_solve( a, pivot,b)1717 SUBROUTINE ros_solve( a, pivot, b) 1844 1718 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1845 1719 ! Template for the forward/backward substitution (using pre-computed LU decomposition) … … 1847 1721 !~~~> input variables 1848 1722 #ifdef full_algebra 1849 REAL(kind=dp), INTENT(in):: a(n,n)1723 REAL(kind=dp), INTENT(IN):: a(n, n) 1850 1724 INTEGER :: ising 1851 1725 #else 1852 REAL(kind=dp), INTENT(in):: a(lu_nonzero)1726 REAL(kind=dp), INTENT(IN):: a(lu_nonzero) 1853 1727 #endif 1854 INTEGER, INTENT(in):: pivot(n)1728 INTEGER, INTENT(IN):: pivot(n) 1855 1729 !~~~> inout variables 1856 REAL(kind=dp), INTENT(inout):: b(n)1730 REAL(kind=dp), INTENT(INOUT):: b(n) 1857 1731 1858 1732 #ifdef full_algebra 1859 1733 CALL DGETRS( 'N',N ,1,A,N,Pivot,b,N,ISING) 1860 IF (info < 0)THEN1861 print* ,"error in dgetrs. ising=",ising1734 IF (info < 0)THEN 1735 print* , "error in dgetrs. ising=", ising 1862 1736 ENDIF 1863 1737 #else 1864 CALL kppsolve( a, b)1738 CALL kppsolve( a, b) 1865 1739 #endif 1866 1740 1867 istatus(nsol) = istatus(nsol) + 11741 istatus(nsol) = istatus(nsol) + 1 1868 1742 1869 1743 END SUBROUTINE ros_solve … … 1893 1767 ! C(i,j) = ros_C( (i-1)*(i-2)/2 + j) 1894 1768 1895 ros_a(1) = (1.0_dp) /g1896 ros_c(1) = (- 2.0_dp) /g1769 ros_a(1) = (1.0_dp) /g 1770 ros_c(1) = (- 2.0_dp) /g 1897 1771 !~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true) 1898 1772 ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE) 1899 ros_newf(1) = . true.1900 ros_newf(2) = . true.1773 ros_newf(1) = .TRUE. 1774 ros_newf(2) = .TRUE. 1901 1775 !~~~> m_i = coefficients for new step solution 1902 ros_m(1) = (3.0_dp) /(2.0_dp*g)1903 ros_m(2) = (1.0_dp) /(2.0_dp*g)1776 ros_m(1) = (3.0_dp) /(2.0_dp* g) 1777 ros_m(2) = (1.0_dp) /(2.0_dp* g) 1904 1778 ! E_i = Coefficients for error estimator 1905 ros_e(1) = 1.0_dp/(2.0_dp* g)1906 ros_e(2) = 1.0_dp/(2.0_dp* g)1907 !~~~> ros_elo = estimator of local order - 1779 ros_e(1) = 1.0_dp/(2.0_dp* g) 1780 ros_e(2) = 1.0_dp/(2.0_dp* g) 1781 !~~~> ros_elo = estimator of local order - the minimum between the 1908 1782 ! main and the embedded scheme orders plus one 1909 1783 ros_elo = 2.0_dp 1910 !~~~> y_stage_i ~ y( t + h* alpha_i)1784 !~~~> y_stage_i ~ y( t + h* alpha_i) 1911 1785 ros_alpha(1) = 0.0_dp 1912 1786 ros_alpha(2) = 1.0_dp 1913 !~~~> gamma_i = \sum_j gamma_{i, j}1787 !~~~> gamma_i = \sum_j gamma_{i, j} 1914 1788 ros_gamma(1) = g 1915 ros_gamma(2) = -g1789 ros_gamma(2) = -g 1916 1790 1917 1791 END SUBROUTINE ros2 … … 1946 1820 !~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true) 1947 1821 ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE) 1948 ros_newf(1) = . true.1949 ros_newf(2) = . true.1950 ros_newf(3) = . false.1822 ros_newf(1) = .TRUE. 1823 ros_newf(2) = .TRUE. 1824 ros_newf(3) = .FALSE. 1951 1825 !~~~> m_i = coefficients for new step solution 1952 1826 ros_m(1) = 0.1e+01_dp … … 1957 1831 ros_e(2) = - 0.29079558716805469821718236208017e+01_dp 1958 1832 ros_e(3) = 0.22354069897811569627360909276199_dp 1959 !~~~> ros_elo = estimator of local order - 1833 !~~~> ros_elo = estimator of local order - the minimum between the 1960 1834 ! main and the embedded scheme orders plus 1 1961 1835 ros_elo = 3.0_dp 1962 !~~~> y_stage_i ~ y( t + h* alpha_i)1836 !~~~> y_stage_i ~ y( t + h* alpha_i) 1963 1837 ros_alpha(1) = 0.0_dp 1964 1838 ros_alpha(2) = 0.43586652150845899941601945119356_dp 1965 1839 ros_alpha(3) = 0.43586652150845899941601945119356_dp 1966 !~~~> gamma_i = \sum_j gamma_{i, j}1840 !~~~> gamma_i = \sum_j gamma_{i, j} 1967 1841 ros_gamma(1) = 0.43586652150845899941601945119356_dp 1968 1842 ros_gamma(2) = 0.24291996454816804366592249683314_dp … … 2007 1881 ros_a(6) = 0.0_dp 2008 1882 2009 ros_c(1) = -0.7137615036412310e+01_dp1883 ros_c(1) = -0.7137615036412310e+01_dp 2010 1884 ros_c(2) = 0.2580708087951457e+01_dp 2011 1885 ros_c(3) = 0.6515950076447975_dp 2012 ros_c(4) = -0.2137148994382534e+01_dp2013 ros_c(5) = -0.3214669691237626_dp2014 ros_c(6) = -0.6949742501781779_dp1886 ros_c(4) = -0.2137148994382534e+01_dp 1887 ros_c(5) = -0.3214669691237626_dp 1888 ros_c(6) = -0.6949742501781779_dp 2015 1889 !~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true) 2016 1890 ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE) 2017 ros_newf(1) = . true.2018 ros_newf(2) = . true.2019 ros_newf(3) = . true.2020 ros_newf(4) = . false.1891 ros_newf(1) = .TRUE. 1892 ros_newf(2) = .TRUE. 1893 ros_newf(3) = .TRUE. 1894 ros_newf(4) = .FALSE. 2021 1895 !~~~> m_i = coefficients for new step solution 2022 1896 ros_m(1) = 0.2255570073418735e+01_dp … … 2025 1899 ros_m(4) = 0.1093502252409163e+01_dp 2026 1900 !~~~> e_i = coefficients for error estimator 2027 ros_e(1) = -0.2815431932141155_dp2028 ros_e(2) = -0.7276199124938920e-01_dp2029 ros_e(3) = -0.1082196201495311_dp2030 ros_e(4) = -0.1093502252409163e+01_dp2031 !~~~> ros_elo = estimator of local order - 1901 ros_e(1) = -0.2815431932141155_dp 1902 ros_e(2) = -0.7276199124938920e-01_dp 1903 ros_e(3) = -0.1082196201495311_dp 1904 ros_e(4) = -0.1093502252409163e+01_dp 1905 !~~~> ros_elo = estimator of local order - the minimum between the 2032 1906 ! main and the embedded scheme orders plus 1 2033 1907 ros_elo = 4.0_dp 2034 !~~~> y_stage_i ~ y( t + h* alpha_i)1908 !~~~> y_stage_i ~ y( t + h* alpha_i) 2035 1909 ros_alpha(1) = 0.0_dp 2036 1910 ros_alpha(2) = 0.1145640000000000e+01_dp 2037 1911 ros_alpha(3) = 0.6552168638155900_dp 2038 1912 ros_alpha(4) = ros_alpha(3) 2039 !~~~> gamma_i = \sum_j gamma_{i, j}1913 !~~~> gamma_i = \sum_j gamma_{i, j} 2040 1914 ros_gamma(1) = 0.5728200000000000_dp 2041 ros_gamma(2) = -0.1769193891319233e+01_dp1915 ros_gamma(2) = -0.1769193891319233e+01_dp 2042 1916 ros_gamma(3) = 0.7592633437920482_dp 2043 ros_gamma(4) = -0.1049021087100450_dp1917 ros_gamma(4) = -0.1049021087100450_dp 2044 1918 2045 1919 END SUBROUTINE ros4 … … 2074 1948 ros_c(1) = 4.0_dp 2075 1949 ros_c(2) = 1.0_dp 2076 ros_c(3) = -1.0_dp1950 ros_c(3) = -1.0_dp 2077 1951 ros_c(4) = 1.0_dp 2078 ros_c(5) = -1.0_dp2079 ros_c(6) = -(8.0_dp/3.0_dp)1952 ros_c(5) = -1.0_dp 1953 ros_c(6) = -(8.0_dp/3.0_dp) 2080 1954 2081 1955 !~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true) 2082 1956 ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE) 2083 ros_newf(1) = . true.2084 ros_newf(2) = . false.2085 ros_newf(3) = . true.2086 ros_newf(4) = . true.1957 ros_newf(1) = .TRUE. 1958 ros_newf(2) = .FALSE. 1959 ros_newf(3) = .TRUE. 1960 ros_newf(4) = .TRUE. 2087 1961 !~~~> m_i = coefficients for new step solution 2088 1962 ros_m(1) = 2.0_dp … … 2095 1969 ros_e(3) = 0.0_dp 2096 1970 ros_e(4) = 1.0_dp 2097 !~~~> ros_elo = estimator of local order - 1971 !~~~> ros_elo = estimator of local order - the minimum between the 2098 1972 ! main and the embedded scheme orders plus 1 2099 1973 ros_elo = 3.0_dp 2100 !~~~> y_stage_i ~ y( t + h* alpha_i)1974 !~~~> y_stage_i ~ y( t + h* alpha_i) 2101 1975 ros_alpha(1) = 0.0_dp 2102 1976 ros_alpha(2) = 0.0_dp 2103 1977 ros_alpha(3) = 1.0_dp 2104 1978 ros_alpha(4) = 1.0_dp 2105 !~~~> gamma_i = \sum_j gamma_{i, j}1979 !~~~> gamma_i = \sum_j gamma_{i, j} 2106 1980 ros_gamma(1) = 0.5_dp 2107 1981 ros_gamma(2) = 1.5_dp … … 2129 2003 ros_s = 6 2130 2004 2131 !~~~> y_stage_i ~ y( t + h* alpha_i)2005 !~~~> y_stage_i ~ y( t + h* alpha_i) 2132 2006 ros_alpha(1) = 0.000_dp 2133 2007 ros_alpha(2) = 0.386_dp … … 2137 2011 ros_alpha(6) = 1.000_dp 2138 2012 2139 !~~~> gamma_i = \sum_j gamma_{i, j}2013 !~~~> gamma_i = \sum_j gamma_{i, j} 2140 2014 ros_gamma(1) = 0.2500000000000000_dp 2141 ros_gamma(2) = -0.1043000000000000_dp2015 ros_gamma(2) = -0.1043000000000000_dp 2142 2016 ros_gamma(3) = 0.1035000000000000_dp 2143 ros_gamma(4) = -0.3620000000000023e-01_dp2017 ros_gamma(4) = -0.3620000000000023e-01_dp 2144 2018 ros_gamma(5) = 0.0_dp 2145 2019 ros_gamma(6) = 0.0_dp … … 2160 2034 ros_a(8) = 0.6019134481288629e+01_dp 2161 2035 ros_a(9) = 0.1253708332932087e+02_dp 2162 ros_a(10) = -0.6878860361058950_dp2036 ros_a(10) = -0.6878860361058950_dp 2163 2037 ros_a(11) = ros_a(7) 2164 2038 ros_a(12) = ros_a(8) … … 2167 2041 ros_a(15) = 1.0_dp 2168 2042 2169 ros_c(1) = -0.5668800000000000e+01_dp2170 ros_c(2) = -0.2430093356833875e+01_dp2171 ros_c(3) = -0.2063599157091915_dp2172 ros_c(4) = -0.1073529058151375_dp2173 ros_c(5) = -0.9594562251023355e+01_dp2174 ros_c(6) = -0.2047028614809616e+02_dp2043 ros_c(1) = -0.5668800000000000e+01_dp 2044 ros_c(2) = -0.2430093356833875e+01_dp 2045 ros_c(3) = -0.2063599157091915_dp 2046 ros_c(4) = -0.1073529058151375_dp 2047 ros_c(5) = -0.9594562251023355e+01_dp 2048 ros_c(6) = -0.2047028614809616e+02_dp 2175 2049 ros_c(7) = 0.7496443313967647e+01_dp 2176 ros_c(8) = -0.1024680431464352e+02_dp2177 ros_c(9) = -0.3399990352819905e+02_dp2050 ros_c(8) = -0.1024680431464352e+02_dp 2051 ros_c(9) = -0.3399990352819905e+02_dp 2178 2052 ros_c(10) = 0.1170890893206160e+02_dp 2179 2053 ros_c(11) = 0.8083246795921522e+01_dp 2180 ros_c(12) = -0.7981132988064893e+01_dp2181 ros_c(13) = -0.3152159432874371e+02_dp2054 ros_c(12) = -0.7981132988064893e+01_dp 2055 ros_c(13) = -0.3152159432874371e+02_dp 2182 2056 ros_c(14) = 0.1631930543123136e+02_dp 2183 ros_c(15) = -0.6058818238834054e+01_dp2057 ros_c(15) = -0.6058818238834054e+01_dp 2184 2058 2185 2059 !~~~> m_i = coefficients for new step solution … … 2201 2075 !~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true) 2202 2076 ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE) 2203 ros_newf(1) = . true.2204 ros_newf(2) = . true.2205 ros_newf(3) = . true.2206 ros_newf(4) = . true.2207 ros_newf(5) = . true.2208 ros_newf(6) = . true.2209 2210 !~~~> ros_elo = estimator of local order - 2077 ros_newf(1) = .TRUE. 2078 ros_newf(2) = .TRUE. 2079 ros_newf(3) = .TRUE. 2080 ros_newf(4) = .TRUE. 2081 ros_newf(5) = .TRUE. 2082 ros_newf(6) = .TRUE. 2083 2084 !~~~> ros_elo = estimator of local order - the minimum between the 2211 2085 ! main and the embedded scheme orders plus 1 2212 2086 ros_elo = 4.0_dp … … 2271 2145 !~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true) 2272 2146 ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE) 2273 ros_newf(1) = . true.2274 ros_newf(2) = . true.2275 ros_newf(3) = . true.2276 ros_newf(4) = . true.2277 2278 !~~~> ros_elo = estimator of local order - 2147 ros_newf(1) = .TRUE. 2148 ros_newf(2) = .TRUE. 2149 ros_newf(3) = .TRUE. 2150 ros_newf(4) = .TRUE. 2151 2152 !~~~> ros_elo = estimator of local order - the minimum between the 2279 2153 ! main and the embedded scheme orders plus 1 2280 2154 ros_elo = 3.0_dp … … 2288 2162 END SUBROUTINE rosenbrock 2289 2163 2290 SUBROUTINE funtemplate( t, y,ydot)2164 SUBROUTINE funtemplate( t, y, ydot) 2291 2165 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 2292 2166 ! Template for the ODE function call. … … 2294 2168 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 2295 2169 !~~~> input variables 2296 REAL(kind=dp):: t, y(nvar)2170 REAL(kind=dp):: t, y(nvar) 2297 2171 !~~~> output variables 2298 2172 REAL(kind=dp):: ydot(nvar) … … 2302 2176 told = time 2303 2177 time = t 2304 CALL fun( y, fix,rconst,ydot)2178 CALL fun( y, fix, rconst, ydot) 2305 2179 time = told 2306 2180 2307 2181 END SUBROUTINE funtemplate 2308 2182 2309 SUBROUTINE jactemplate( t, y,jcb)2183 SUBROUTINE jactemplate( t, y, jcb) 2310 2184 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 2311 2185 ! Template for the ODE Jacobian call. … … 2313 2187 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 2314 2188 !~~~> input variables 2315 REAL(kind=dp):: t, y(nvar)2189 REAL(kind=dp):: t, y(nvar) 2316 2190 !~~~> output variables 2317 2191 #ifdef full_algebra 2318 REAL(kind=dp):: jv(lu_nonzero), jcb(nvar,nvar)2192 REAL(kind=dp):: jv(lu_nonzero), jcb(nvar, nvar) 2319 2193 #else 2320 2194 REAL(kind=dp):: jcb(lu_nonzero) … … 2323 2197 REAL(kind=dp):: told 2324 2198 #ifdef full_algebra 2325 INTEGER :: i, j2199 INTEGER :: i, j 2326 2200 #endif 2327 2201 … … 2329 2203 time = t 2330 2204 #ifdef full_algebra 2331 CALL jac_sp(y, fix,rconst,jv)2332 DO j=1, nvar2333 DO i=1, nvar2334 jcb(i, j) = 0.0_dp2205 CALL jac_sp(y, fix, rconst, jv) 2206 DO j=1, nvar 2207 DO i=1, nvar 2208 jcb(i, j) = 0.0_dp 2335 2209 ENDDO 2336 2210 ENDDO 2337 DO i=1, lu_nonzero2338 jcb(lu_irow(i), lu_icol(i)) = jv(i)2211 DO i=1, lu_nonzero 2212 jcb(lu_irow(i), lu_icol(i)) = jv(i) 2339 2213 ENDDO 2340 2214 #else 2341 CALL jac_sp( y, fix,rconst,jcb)2215 CALL jac_sp( y, fix, rconst, jcb) 2342 2216 #endif 2343 2217 time = told … … 2345 2219 END SUBROUTINE jactemplate 2346 2220 2347 SUBROUTINE chem_gasphase_integrate (time_step_len,conc,tempk,photo,ierrf,xnacc,xnrej,istatus,l_debug,pe) 2221 SUBROUTINE kppdecomp( jvs, ier) 2222 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 2223 ! sparse lu factorization 2224 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 2225 ! loop expansion generated by kp4 2226 2227 INTEGER :: ier 2228 REAL(kind=dp):: jvs(lu_nonzero), w(nvar), a 2229 INTEGER :: k, kk, j, jj 2230 2231 a = 0. 2232 ier = 0 2233 2234 ! i = 1 2235 ! i = 2 2236 RETURN 2237 2238 END SUBROUTINE kppdecomp 2239 2240 SUBROUTINE chem_gasphase_integrate (time_step_len, conc, tempi, qvapi, fakti, photo, ierrf, xnacc, xnrej, istatus, l_debug, pe, & 2241 icntrl_i, rcntrl_i) 2348 2242 2349 2243 IMPLICIT NONE 2350 2244 2351 REAL(dp), INTENT(IN) :: time_step_len 2352 REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: conc 2353 REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: photo 2354 REAL(dp), DIMENSION(:), INTENT(IN) :: tempk 2355 INTEGER, INTENT(OUT), OPTIONAL :: ierrf(:) 2356 INTEGER, INTENT(OUT), OPTIONAL :: xnacc(:) 2357 INTEGER, INTENT(OUT), OPTIONAL :: xnrej(:) 2358 INTEGER, INTENT(INOUT), OPTIONAL :: istatus(:) 2359 INTEGER, INTENT(IN), OPTIONAL :: pe 2360 LOGICAL, INTENT(IN), OPTIONAL :: l_debug 2245 REAL(dp), INTENT(IN) :: time_step_len 2246 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: conc 2247 REAL(dp), DIMENSION(:, :), INTENT(IN) :: photo 2248 REAL(dp), DIMENSION(:), INTENT(IN) :: tempi 2249 REAL(dp), DIMENSION(:), INTENT(IN) :: qvapi 2250 REAL(dp), DIMENSION(:), INTENT(IN) :: fakti 2251 INTEGER, INTENT(OUT), OPTIONAL :: ierrf(:) 2252 INTEGER, INTENT(OUT), OPTIONAL :: xnacc(:) 2253 INTEGER, INTENT(OUT), OPTIONAL :: xnrej(:) 2254 INTEGER, INTENT(INOUT), OPTIONAL :: istatus(:) 2255 INTEGER, INTENT(IN), OPTIONAL :: pe 2256 LOGICAL, INTENT(IN), OPTIONAL :: l_debug 2257 INTEGER, DIMENSION(nkppctrl), INTENT(IN), OPTIONAL :: icntrl_i 2258 REAL(dp), DIMENSION(nkppctrl), INTENT(IN), OPTIONAL :: rcntrl_i 2361 2259 2362 2260 INTEGER :: k ! loop variable 2363 REAL(dp) :: dt 2364 INTEGER, DIMENSION(20) :: istatus_u 2365 INTEGER :: ierr_u 2366 2367 2368 if (present (istatus)) istatus = 0 2369 2370 DO k=1,vl_glo,vl_dim 2261 REAL(dp) :: dt 2262 INTEGER, DIMENSION(20) :: istatus_u 2263 INTEGER :: ierr_u 2264 INTEGER :: istatf 2265 INTEGER :: vl_dim_lo 2266 2267 2268 IF (PRESENT (istatus)) istatus = 0 2269 IF (PRESENT (icntrl_i)) icntrl = icntrl_i 2270 IF (PRESENT (rcntrl_i)) rcntrl = rcntrl_i 2271 2272 vl_glo = size(tempi, 1) 2273 2274 vl_dim_lo = vl_dim 2275 DO k=1, vl_glo, vl_dim_lo 2371 2276 is = k 2372 ie = min(k+vl_dim-1,vl_glo) 2373 vl = ie-is+1 2374 2375 c(:) = conc(is,:) 2376 2377 temp = tempk(is) 2378 2379 phot(:) = photo(is,:) 2277 ie = min(k+ vl_dim_lo-1, vl_glo) 2278 vl = ie-is+ 1 2279 2280 c(:) = conc(is, :) 2281 2282 temp = tempi(is) 2283 2284 qvap = qvapi(is) 2285 2286 fakt = fakti(is) 2287 2288 CALL initialize 2289 2290 phot(:) = photo(is, :) 2380 2291 2381 2292 CALL update_rconst … … 2384 2295 2385 2296 ! integrate from t=0 to t=dt 2386 CALL integrate(0._dp, dt,icntrl,rcntrl,istatus_u = istatus_u,ierr_u=ierr_u) 2387 2388 2389 IF (PRESENT(l_debug) .AND. PRESENT(pe)) THEN 2390 IF (l_debug) CALL error_output(conc(is,:),ierr_u,pe) 2391 ENDIF 2392 conc(is,:) = c(:) 2393 2394 ! Return Diagnostic Information 2395 2396 if(PRESENT(ierrf)) ierrf(is) = ierr_u 2397 if(PRESENT(xnacc)) xnacc(is) = istatus_u(4) 2398 if(PRESENT(xnrej)) xnrej(is) = istatus_u(5) 2399 2400 if (PRESENT (istatus)) then 2401 istatus(1:8) = istatus(1:8) + istatus_u(1:8) 2297 CALL integrate(0._dp, dt, icntrl, rcntrl, istatus_u = istatus_u, ierr_u=ierr_u) 2298 2299 2300 IF (PRESENT(l_debug) .AND. PRESENT(pe)) THEN 2301 IF (l_debug) CALL error_output(conc(is, :), ierr_u, pe) 2302 ENDIF 2303 2304 conc(is, :) = c(:) 2305 2306 ! RETURN diagnostic information 2307 2308 IF (PRESENT(ierrf)) ierrf(is) = ierr_u 2309 IF (PRESENT(xnacc)) xnacc(is) = istatus_u(4) 2310 IF (PRESENT(xnrej)) xnrej(is) = istatus_u(5) 2311 2312 IF (PRESENT (istatus)) THEN 2313 istatus(1:8) = istatus(1:8) + istatus_u(1:8) 2402 2314 ENDIF 2403 2315 … … 2407 2319 ! Deallocate input arrays 2408 2320 2409 if (ALLOCATED(temp)) DEALLOCATE(temp) 2410 2411 data_loaded = .false. 2321 2322 data_loaded = .FALSE. 2412 2323 2413 2324 RETURN 2414 2325 END SUBROUTINE chem_gasphase_integrate 2415 2416 SUBROUTINE fill_temp(status,array)2417 2418 INTEGER, INTENT(OUT) :: status2419 REAL(dp), INTENT(IN), DIMENSION(:) :: array2420 2421 status = 02422 IF (.not. ALLOCATED(temp)) &2423 ALLOCATE(temp(size(array)))2424 2425 IF (data_loaded .AND. (vl_glo /= size(array,1))) THEN2426 status = 12427 RETURN2428 END IF2429 2430 vl_glo = size(array,1)2431 temp = array2432 data_loaded = .TRUE.2433 2434 RETURN2435 2436 END SUBROUTINE fill_temp2437 2326 2438 2327 END MODULE chem_gasphase_mod -
palm/trunk/SOURCE/data_output_2d.f90
r3294 r3298 25 25 ! ----------------- 26 26 ! $Id$ 27 ! minor formatting (kanani) 28 ! chem_data_output_2d subroutine added (basit) 29 ! 30 ! 3294 2018-10-01 02:37:10Z raasch 27 31 ! changes concerning modularization of ocean option 28 32 ! … … 252 256 ONLY: bulk_cloud_model, bcm_data_output_2d 253 257 258 USE chemistry_model_mod, & 259 ONLY: chem_data_output_2d 260 254 261 USE control_parameters, & 255 ONLY: data_output_2d_on_each_pe, data_output_xy,&262 ONLY: air_chemistry, data_output_2d_on_each_pe, data_output_xy, & 256 263 data_output_xz, data_output_yz, do2d, & 257 264 do2d_xy_last_time, do2d_xy_time_count, & … … 265 272 USE cpulog, & 266 273 ONLY: cpu_log, log_point 267 274 268 275 USE gust_mod, & 269 276 ONLY: gust_data_output_2d, gust_module_enabled … … 1326 1333 ENDIF 1327 1334 1335 IF ( .NOT. found .AND. air_chemistry ) THEN 1336 CALL chem_data_output_2d( av, do2d(av,if), found, grid, mode, & 1337 local_pf, two_d, nzb_do, nzt_do, fill_value ) 1338 ENDIF 1328 1339 ! 1329 1340 !-- User defined quantities -
palm/trunk/SOURCE/data_output_mask.f90
r3274 r3298 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Minor formatting (kanani) 28 ! Call for chem_data_output_mask added (basit) 29 ! 30 ! 3274 2018-09-24 15:42:55Z knoop 27 31 ! Modularization of all bulk cloud physics code components 28 32 ! … … 133 137 ONLY: e, nc, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, sa, & 134 138 tend, u, v, vpt, w, d_exner 135 139 136 140 USE averaging, & 137 141 ONLY: e_av, lpt_av, nc_av, nr_av, p_av, pc_av, pr_av, pt_av, q_av, & … … 141 145 USE basic_constants_and_equations_mod, & 142 146 ONLY: lv_d_cp 143 147 148 USE chemistry_model_mod, & 149 ONLY: chem_data_output_mask 150 144 151 USE control_parameters, & 145 ONLY: domask, domask_no, domask_time_count, mask_i,&152 ONLY: air_chemistry, domask, domask_no, domask_time_count, mask_i, & 146 153 mask_j, mask_k, mask_size, mask_size_l, mask_start_l, & 147 154 max_masks, message_string, mid, nz_do3d, simulated_time … … 151 158 USE indices, & 152 159 ONLY: nbgp, nxl, nxr, nyn, nys, nzb 153 160 154 161 USE kinds 155 162 … … 524 531 525 532 CASE DEFAULT 526 527 533 ! 528 534 !-- Radiation quantity … … 532 538 ENDIF 533 539 540 IF ( air_chemistry ) THEN 541 CALL chem_data_output_mask(av, domask(mid,av,ivar), found, & 542 local_pf ) 543 ENDIF 534 544 ! 535 545 !-- User defined quantity -
palm/trunk/SOURCE/date_and_time_mod.f90
r2718 r3298 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Minor formatting (kanani) 28 ! - Added Routines for DEFAULT mode of chemistry emissions (Russo) 29 ! - Added routine for reading-in real dates: format has to be in DDMMYYYY and 30 ! passed in the namelist parameter date_init (Russo) 31 ! - Added calculation of several time indices useful in several routines 32 ! of the model (Russo) 33 ! 34 ! 2718 2018-01-02 08:49:38Z maronga 27 35 ! Corrected "Former revisions" section 28 36 ! … … 46 54 !> This routine calculates all needed information on date and time used by 47 55 !> other modules 56 !> @todo Further testing and revision of routines for updating indices of 57 !> emissions in the default mode 58 !> @todo Add routine for recognizing leap years 59 !> @todo Add recognition of exact days of week (Monday, Tuesday, etc.) 60 !> @todo Reconsider whether to remove day_of_year_init from the namelist: we 61 !> already implemented changes for calculating it from date_init in 62 !> calc_date_and_time 48 63 !------------------------------------------------------------------------------! 49 64 MODULE date_and_time_mod 50 51 USE control_parameters, & 52 ONLY: time_since_reference_point 53 65 66 USE control_parameters, & 67 ONLY: coupling_start_time, days_since_reference_point, & 68 message_string, simulated_time, time_since_reference_point 69 54 70 USE kinds 55 71 72 56 73 IMPLICIT NONE 57 74 58 75 PRIVATE 59 60 PUBLIC calc_date_and_time, d_hours_day, d_seconds_hour, d_seconds_year, & 61 day_of_year, day_of_year_init, time_utc, time_utc_init 62 63 64 INTEGER(iwp) :: day_of_year !< day of the year (DOY) 65 INTEGER(iwp) :: day_of_year_init = 172 !< DOY at model start (default: 21 June) 66 67 REAL(wp) :: time_utc !< current model time in UTC 68 REAL(wp) :: time_utc_init = 43200.0_wp !< UTC time at model start 69 76 77 !-- Variables Declaration 78 79 INTEGER(iwp) :: day_of_year = 0 !< day of the year (DOY) 80 INTEGER(iwp) :: day_of_year_init = 0 !< DOY at model start (default: 0) 81 82 ! --- Most of these indices are updated by the routine calc_date_and_time according to the current date and time of the simulation 83 INTEGER(iwp) :: hour_of_year = 1 !< hour of the current year (1:8760(8784)) 84 INTEGER(iwp) :: hour_of_day=1 !< hour of the current day (1:24) 85 INTEGER(iwp) :: day_of_month=0 !< day of the current month (1:31) 86 INTEGER(iwp) :: month_of_year=0 !< month of the current year (1:12) 87 INTEGER(iwp) :: current_year=0 !< current year 88 INTEGER(iwp) :: hour_call_emis=0 !< index used to call the emissions just once every hour 89 90 INTEGER(iwp) :: ihour 91 INTEGER(iwp) :: index_mm !< index months of the default emission mode 92 INTEGER(iwp) :: index_dd !< index days of the default emission mode 93 INTEGER(iwp) :: index_hh !< index hours of the default emission mode 94 95 REAL(wp) :: time_utc !< current model time in UTC 96 REAL(wp) :: time_utc_emis !< current emission module time in UTC 97 REAL(wp) :: time_utc_init = 43200.0_wp !< UTC time at model start 98 REAL(wp) :: time_update !< used to calculate actual second of the simulation 99 70 100 REAL(wp), PARAMETER :: d_hours_day = 1.0_wp / 24.0_wp !< inverse of hours per day (1/24) 71 101 REAL(wp), PARAMETER :: d_seconds_hour = 1.0_wp / 3600.0_wp !< inverse of seconds per hour (1/3600) 72 102 REAL(wp), PARAMETER :: d_seconds_year = 1.0_wp / 31536000.0_wp !< inverse of the seconds per year (1/(365*86400)) 73 103 104 CHARACTER(len=8) :: date_init = "31122018" !< Starting date of simulation: We selected this because it was a monday 105 106 !> --- Parameters 107 INTEGER, PARAMETER, DIMENSION(12) :: days = (/31,28,31,30,31,30,31,31,30,31,30,31/) ! total number of days for each month (no leap year) 108 74 109 SAVE 75 110 111 !-- INTERFACES PART 112 !-- Read initial day and time of simulation 113 INTERFACE init_date_and_time 114 MODULE PROCEDURE init_date_and_time 115 END INTERFACE init_date_and_time 116 117 !-- Get hour index in the DEAFULT case of chemistry emissions : 118 INTERFACE time_default_indices 119 MODULE PROCEDURE time_mdh_indices 120 MODULE PROCEDURE time_hour_indices 121 END INTERFACE time_default_indices 122 123 !-- Calculate current date and time 124 INTERFACE calc_date_and_time 125 MODULE PROCEDURE calc_date_and_time 126 END INTERFACE 127 128 129 !-- Public Interfaces 130 PUBLIC calc_date_and_time, time_default_indices, init_date_and_time 131 132 !-- Public Variables 133 PUBLIC date_init, d_hours_day, d_seconds_hour, d_seconds_year, & 134 day_of_year, day_of_year_init, time_utc, time_utc_init, day_of_month, & 135 month_of_year, index_mm, index_dd, index_hh, hour_of_day, hour_of_year, & 136 hour_call_emis 137 76 138 CONTAINS 139 140 141 !------------------------------------------------------------------------------! 142 !> Reads starting date from namelist 143 !------------------------------------------------------------------------------! 144 145 SUBROUTINE init_date_and_time 146 147 IMPLICIT NONE 148 149 IF (day_of_year_init == 0) THEN 150 ! Day of the month at starting time 151 READ(UNIT=date_init(1:2),fmt=*)day_of_month 152 153 ! Month of the year at starting time 154 READ(UNIT=date_init(3:4),fmt=*)month_of_year 155 156 ! Year at starting time 157 READ(UNIT=date_init(5:8),fmt=*)current_year 158 159 ENDIF 160 161 END SUBROUTINE init_date_and_time 77 162 78 163 !------------------------------------------------------------------------------! 79 164 ! Description: 80 165 ! ------------ 81 !> Calculate current da y of the year and time in UTC166 !> Calculate current date and time of the simulation 82 167 !------------------------------------------------------------------------------! 83 168 … … 86 171 IMPLICIT NONE 87 172 88 ! 89 !-- Calculate current day of the year 90 day_of_year = day_of_year_init + INT(FLOOR( (time_utc_init + time_since_reference_point)& 91 / 86400.0_wp ), KIND=iwp) 92 173 !-- Variables Definition 174 INTEGER :: i_mon !< Index for going through the different months 175 176 !> Update simulation time in seconds 177 time_update = simulated_time-coupling_start_time 178 179 !-- Calculate current day of the simulated time 180 days_since_reference_point=INT(FLOOR( (time_utc_init + time_update) & 181 / 86400.0_wp ) ) 182 183 !-- Calculate actual UTC time 93 184 time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp) 94 185 186 !sB PRILIMINARY workaround for time_utc bug concerning time_since_reference_point: 187 time_utc_emis = MOD((time_utc_init + time_update), 86400.0_wp) 188 189 !-- Calculate initial day of the year: it is calculated only once. In fact, day_of_year_init is initialized to 0 and then a positive value is passed. This condition is also called only when day_of_year_init is not given in the namelist. 190 191 IF ( day_of_year_init == 0 ) THEN 192 193 !> Condition for printing an error when date_init is not provided when day_of_year_init is not given in the namelist or when the format of the date is not the one required by PALM. 194 IF ( day_of_month .GT. 0 .AND. day_of_month .LE. 31 .AND. month_of_year .GT. 0 .AND. month_of_year .LE. 12) THEN 195 196 IF ( month_of_year == 1 ) THEN !!month of year is read in input 197 198 day_of_year_init = day_of_month 199 200 ELSE 201 202 day_of_year_init= SUM(days( 1:(month_of_year-1) )) + day_of_month !day_of_month is read in input in this case 203 204 ENDIF 205 !kanani: Revise, we cannot force users to provide date_init, maybe set a default value? 206 ! ELSE 207 ! 208 ! message_string = 'date_init not provided in the namelist or' // & 209 ! ' given in the wrong format: MUST BE DDMMYYYY' 210 ! CALL message( 'calc_date_and_time', 'DT0100', 2, 2, 0, 6, 0 ) 211 212 ENDIF 213 214 ENDIF 215 216 !-- Calculate actual hour of the day: the first hour of the day is from 00:00:00 to 00:59:59. 217 218 hour_of_day = INT( FLOOR( time_utc_emis/3600.0_wp ) ) + 1 219 220 !-- Calculate current day of the year !TBD: considetr leap years 221 IF ( (day_of_year_init + days_since_reference_point) .GT. 365 ) THEN 222 223 day_of_year=INT(MOD((day_of_year_init + days_since_reference_point), 365.0_wp)) 224 225 ELSE 226 227 day_of_year = day_of_year_init + days_since_reference_point 228 229 ENDIF 230 231 ! 232 !-- Calculate current hour of the year 233 hour_of_year = ( (day_of_year-1) * 24 ) + hour_of_day !> actual hour of the year 234 235 236 ! 237 !-- UPDATE actual day of the month and month of the year 238 !> -------------------------------------------------------------------------------- 239 !> The first case is when date_init is not provided: we only know day_of_year_init 240 IF ( month_of_year == 0 .AND. day_of_month == 0) THEN 241 242 !> The first case is when date_init is not provided: we only know day_of_year_init 243 !DO i_mon=1,12 244 !IF (day_of_year .LE. SUM(days(1:i_mon))) THEN 245 IF ( day_of_year .LE. 31 ) THEN 246 247 month_of_year=1 248 day_of_month=day_of_year 249 250 ELSE 251 252 DO i_mon=2,12 !january is considered in the first case 253 IF ( day_of_year .LE. SUM(days(1:i_mon)) .AND. day_of_year .GT. SUM(days(1:(i_mon-1))) ) THEN 254 255 month_of_year=i_mon 256 257 day_of_month=INT(MOD(day_of_year, SUM(days(1:(i_mon-1))))) 258 259 GOTO 38 260 261 ENDIF 262 263 38 ENDDO 264 ENDIF 265 !> -------------------------------------------------------------------------------- 266 !> in the second condition both day of month and month_of_year are either given in input (passed to date_init) or we are in some day successive to the initial one, so that day_of_month has already be computed in previous step 267 !>TBD: something to calculate the current year is missing 268 ELSEIF ( day_of_month .GT. 0 .AND. day_of_month .LE. 31 .AND. month_of_year .GT. 0 .AND. month_of_year .LE. 12) THEN 269 270 !> calculate month_of_year. TBD: test the condition when day_of_year==31 271 272 IF (day_of_year==1) THEN !> this allows to turn from december to January when passing from a year to another 273 274 month_of_year = 1 275 276 ELSE IF (day_of_year .GT. 1 .AND. day_of_year .GT. SUM(days(1:month_of_year))) THEN 277 278 month_of_year = month_of_year + 1 279 280 ENDIF 281 282 !> calculate day_of_month 283 IF ( month_of_year == 1 ) THEN 284 285 day_of_month=day_of_year 286 287 ELSE 288 289 day_of_month=INT(MOD(day_of_year, SUM(days(1:(month_of_year-1))))) 290 291 ENDIF 292 293 294 295 296 ELSE 297 298 !> Condition when date_init is provided but it is given in the wrong format 299 message_string = 'date_init not provided in the namelist or' // & 300 ' given in the wrong format: MUST BE DDMMYYYY' 301 CALL message( 'calc_date_and_time', 'DT0101', 2, 2, 0, 6, 0 ) 302 303 ENDIF 304 305 END SUBROUTINE calc_date_and_time 306 307 !TBD: check the routines used for update of emission indices in the DEFAULT MODE 308 !------------------------------------------------------------------------------! 309 ! Description: 310 ! ------------ 311 !> This routine determines the time factor index in the mdh case of the DEFAULT chemistry emissions mode. 312 !------------------------------------------------------------------------------! 313 314 SUBROUTINE time_mdh_indices(daytype_mdh,mo, dd, hh, index_mm, index_dd, index_hh) 315 316 USE indices 317 318 IMPLICIT NONE 319 320 !> IN/OUTPUT 321 INTEGER, INTENT(INOUT) :: mo !> Month of year 322 INTEGER, INTENT(INOUT) :: dd !> Day of month 323 INTEGER, INTENT(INOUT) :: hh !> Hour of day 324 INTEGER, INTENT(INOUT) :: index_mm !> Index Month 325 INTEGER, INTENT(INOUT) :: index_dd !> Index Day 326 INTEGER, INTENT(INOUT) :: index_hh !> Index Hour 327 328 CHARACTER(len=80), INTENT(INOUT) :: daytype_mdh !> type of the day in mdh mode: one of 1-WORKDAY 329 ! 2-WEEKEND 330 ! 3-HOLIDAY 331 332 REAL(wp) :: frac_day=0 333 334 !> ------------------------------------------------------------------------ 335 336 INTEGER :: weekday 337 338 !> CONSTANTS 339 INTEGER, PARAMETER :: nmonth = 12 340 INTEGER, PARAMETER :: nday = 7 341 INTEGER, PARAMETER :: nhour = 24 342 343 frac_day= (dd-1)/nday !> indicates the week of the month, supposing the month starts on monday 344 345 ! 1:7 1:31 7 (0:30)/7 346 weekday = dd-( nday * (INT( CEILING( frac_day ) ) ) ) ! for now we let the year start on Monday. 347 348 !TBD: set weekday correct based on date 349 index_mm = mo 350 index_dd = nmonth + weekday !> Index of the days in the mdh mode (13:20) 351 352 SELECT CASE(TRIM(daytype_mdh)) 353 354 CASE ("workday") 355 356 index_hh = nmonth+ nday + hh 357 358 CASE ("weekend") 359 360 index_hh = nmonth+ nday + nhour + hh 361 362 CASE ("holiday") 363 364 index_hh = nmonth+ nday + 2*nhour + hh 365 END SELECT 366 367 368 END SUBROUTINE time_mdh_indices 369 370 !------------------------------------------------------------------------------! 371 ! Description: 372 ! ------------ 373 !> This routine determines the time factor index in the HOURLY case of the DEFAULT emissions mode. 374 !------------------------------------------------------------------------------! 375 376 SUBROUTINE time_hour_indices(mo,dd,hh,index_hh) 377 378 USE indices 379 380 IMPLICIT NONE 381 382 !> IN/OUTPUT 383 INTEGER, INTENT(INOUT) :: mo !> Month 384 INTEGER, INTENT(INOUT) :: hh !> Hour 385 INTEGER, INTENT(INOUT) :: dd !> Day 386 INTEGER, INTENT(INOUT) :: index_hh !> Index Hour 387 388 !> Additional Variables for calculateing indices 389 INTEGER :: index_mm !> Index Month 390 INTEGER :: index_dd !> Index Day 391 INTEGER :: i_mon !> Index for going through the different months 392 INTEGER :: sum_dd !> Sum days 393 394 !> CONSTANTS 395 INTEGER, PARAMETER :: nmonth=12 396 INTEGER, PARAMETER :: nday = 7 397 INTEGER, PARAMETER :: nhour = 24 398 INTEGER, PARAMETER, DIMENSION(12) :: days = (/31,28,31,30,31,30,31,31,30,31,30,31/) ! no leap year 399 400 401 index_mm = mo-1 402 index_dd = dd-1 403 sum_dd=0 404 405 IF (mo == 1) THEN 406 407 index_hh=(index_dd*nhour)+hh 408 409 ELSE 410 411 DO i_mon=1,index_mm 95 412 96 97 END SUBROUTINE calc_date_and_time 98 99 100 101 END MODULE date_and_time_mod 413 sum_dd=sum_dd+days(i_mon) 414 415 ENDDO 416 417 index_hh=(sum_dd*nhour)+(index_dd*nhour)+(hh) 418 419 ENDIF 420 421 422 END SUBROUTINE time_hour_indices 423 424 425 END MODULE date_and_time_mod -
palm/trunk/SOURCE/flow_statistics.f90
r3294 r3298 19 19 ! 20 20 ! Current revisions: 21 ! ----------------- 21 ! ------------------ 22 22 ! 23 23 ! … … 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Minor formatting (kanani) 28 ! - Added .AND. max_pr_cs > 0 before MPI_ALLREDUCE call (forkel) 29 ! - Data arrays, sums, sums_l, hom, hom_sum updated for chem species (basit) 30 ! - Directives of parallelism added for chemistry (basit) 31 ! - Call for chem_statistics added (basit) 32 ! 33 ! 3294 2018-10-01 02:37:10Z raasch 27 34 ! ocean renamed ocean_mode 28 35 ! … … 283 290 !------------------------------------------------------------------------------! 284 291 SUBROUTINE flow_statistics 285 292 286 293 287 294 USE arrays_3d, & … … 291 298 sa, u, ug, v, vg, vpt, w, w_subs, waterflux_output_conversion, & 292 299 zw, d_exner 293 300 294 301 USE basic_constants_and_equations_mod, & 295 ONLY: g, lv_d_cp 296 302 ONLY: g, lv_d_cp 303 304 305 USE chem_modules, & 306 ONLY: max_pr_cs 307 308 USE chemistry_model_mod, & 309 ONLY: chem_species, chem_statistics 310 297 311 USE control_parameters, & 298 ONLY: a verage_count_pr, cloud_droplets, do_sum,&312 ONLY: air_chemistry, average_count_pr, cloud_droplets, do_sum, & 299 313 dt_3d, humidity, initializing_actions, land_surface, & 300 314 large_scale_forcing, large_scale_subsidence, max_pr_user, & … … 303 317 use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, & 304 318 ws_scheme_mom, ws_scheme_sca 305 319 306 320 USE cpulog, & 307 321 ONLY: cpu_log, log_point 308 322 309 323 USE grid_variables, & 310 324 ONLY: ddx, ddy … … 1780 1794 ENDIF 1781 1795 ! 1796 !-- Calculate the chemistry module profiles 1797 IF ( air_chemistry ) THEN 1798 CALL chem_statistics( 'profiles', sr, tn ) 1799 ENDIF 1800 ! 1782 1801 !-- Calculate the user-defined profiles 1783 1802 CALL user_statistics( 'profiles', sr, tn ) … … 1797 1816 sums_l(:,pr_palm+1:pr_palm+max_pr_user,i) 1798 1817 ENDIF 1818 1819 IF ( air_chemistry ) THEN 1820 IF ( max_pr_cs > 0 ) THEN 1821 sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+ max_pr_cs,0) = & 1822 sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs,0) + & 1823 sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs,i) 1824 1825 ENDIF 1826 ENDIF 1799 1827 ENDDO 1800 1828 ENDIF … … 1811 1839 MPI_REAL, MPI_SUM, comm2d, ierr ) 1812 1840 ENDIF 1841 1842 IF ( air_chemistry .AND. max_pr_cs > 0 ) THEN 1843 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1844 CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+1,0), sums(nzb,pr_palm+1), & 1845 nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr ) 1846 ENDIF 1847 1813 1848 #else 1814 1849 sums = sums_l(:,:,0) … … 1875 1910 ngp_2dh_s_inner(k,sr) 1876 1911 ENDDO 1912 ENDIF 1913 1914 IF ( air_chemistry ) THEN 1915 IF ( max_pr_cs > 0 ) THEN 1916 DO k = nzb, nzt+1 1917 sums(k, pr_palm+1:pr_palm+max_pr_user+max_pr_cs) = & 1918 sums(k, pr_palm+1:pr_palm+max_pr_user+max_pr_cs) / & 1919 ngp_2dh_s_inner(k,sr) 1920 ENDDO 1921 ENDIF 1877 1922 ENDIF 1878 1923 … … 2018 2063 ENDIF 2019 2064 2065 IF ( air_chemistry ) THEN 2066 IF ( max_pr_cs > 0 ) THEN ! chem_spcs profiles 2067 hom(:, 1, pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs, sr) = & 2068 sums(:, pr_palm+max_pr_user+1:pr_palm+max_pr_user+max_pr_cs) 2069 ENDIF 2070 ENDIF 2020 2071 ! 2021 2072 !-- Determine the boundary layer height using two different schemes. -
palm/trunk/SOURCE/header.f90
r3294 r3298 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Minor formatting (kanani) 28 ! Add chemistry header (basit) 29 ! 30 ! 3294 2018-10-01 02:37:10Z raasch 27 31 ! changes concerning modularization of ocean option 28 32 ! … … 391 395 ONLY: bulk_cloud_model, bcm_header 392 396 397 USE chemistry_model_mod, & 398 ONLY: chem_header 399 393 400 USE control_parameters 394 401 395 402 USE cpulog, & 396 403 ONLY: log_point_s 397 404 398 405 USE date_and_time_mod, & 399 406 ONLY: day_of_year_init, time_utc_init … … 1123 1130 IF ( gust_module_enabled ) CALL gust_header ( io ) 1124 1131 1132 IF ( air_chemistry ) CALL chem_header ( io ) 1125 1133 ! 1126 1134 !-- Boundary conditions -
palm/trunk/SOURCE/init_3d_model.f90
- Property svn:mergeinfo changed
r3294 r3298 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Minor formatting (kanani) 28 ! Added CALL to the chem_emissions module (Russo) 29 ! 30 ! 3294 2018-10-01 02:37:10Z raasch 27 31 ! changes concerning modularization of ocean option 28 32 ! … … 501 505 !------------------------------------------------------------------------------! 502 506 SUBROUTINE init_3d_model 503 507 504 508 505 509 USE advec_ws … … 514 518 ONLY: bulk_cloud_model, bcm_init, bcm_init_arrays 515 519 516 USE chemistry_model_mod, & 517 ONLY: chem_emissions 520 USE chem_emissions_mod, & 521 ONLY: chem_emissions_init 522 523 USE chem_modules, & 524 ONLY: do_emis, max_pr_cs, nspec_out 518 525 519 526 USE control_parameters 520 527 521 528 USE flight_mod, & 522 529 ONLY: flight_init 523 530 524 531 USE grid_variables, & 525 532 ONLY: dx, dy, ddx2_mg, ddy2_mg … … 551 558 552 559 USE netcdf_data_input_mod, & 553 ONLY: init_3d, netcdf_data_input_init_3d 560 ONLY: chem_emis, chem_emis_att, init_3d, & 561 netcdf_data_input_init_3d, netcdf_data_input_interpolate 554 562 555 563 USE ocean_mod, & 556 564 ONLY: ocean_init, ocean_init_arrays 557 565 558 566 USE particle_attributes, & 559 567 ONLY: particle_advection 560 568 561 569 USE pegrid 562 570 563 571 USE plant_canopy_model_mod, & 564 572 ONLY: pcm_init … … 2261 2269 ! 2262 2270 !-- If required, set chemical emissions 2263 !-- (todo(FK): This should later on be CALLed time-dependently in init_3d_model) 2264 IF ( air_chemistry ) THEN 2265 CALL chem_emissions 2266 ENDIF 2267 2271 !-- Initialize values of cssws according to chemistry emission values 2272 IF ( air_chemistry .AND. do_emis ) THEN 2273 CALL chem_emissions_init( chem_emis_att, chem_emis, nspec_out ) 2274 ENDIF 2268 2275 ! 2269 2276 !-- Initialize radiation processes -
palm/trunk/SOURCE/modules.f90
r3294 r3298 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Minor formatting/clean-up (kanani) 28 ! Added some variables for time treatment (Russo) 29 ! 30 ! 3294 2018-10-01 02:37:10Z raasch 27 31 ! ocean renamed ocean_mode 28 32 ! … … 1401 1405 REAL(wp) :: cos_alpha_surface !< cosine of alpha_surface 1402 1406 REAL(wp) :: coupling_start_time = 0.0_wp !< namelist parameter 1407 REAL(wp) :: days_since_reference_point = 0.0_wp !< days after atmosphere-ocean coupling has been activated, 1408 !< or after spinup phase of LSM has been finished 1403 1409 REAL(wp) :: disturbance_amplitude = 0.25_wp !< namelist parameter 1404 1410 REAL(wp) :: disturbance_energy_limit = 0.01_wp !< namelist parameter -
palm/trunk/SOURCE/netcdf_data_input_mod.f90
r3257 r3298 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Modified EXPERT mode to PRE-PROCESSED mode (Russo) 28 ! Introduced Chemistry static netcdf file (Russo) 29 ! Added the routine for reading-in netcdf data for chemistry (Russo) 30 ! Added routines to get_variable interface specific for chemistry files (Russo) 31 ! 32 ! 3257 2018-09-17 17:11:46Z suehring 27 33 ! Adjust checks for building_type and building_id, which is necessary after 28 34 ! topography filtering (building_type and id can be modified by the filtering). … … 189 195 !> Modulue contains routines to input data according to Palm input data 190 196 !> standart using dynamic and static input files. 191 !> 197 !> @todo - Review Reading of netcdf files for chemistry 192 198 !> @todo - Order input alphabetically 193 199 !> @todo - Revise error messages and error numbers … … 341 347 342 348 END TYPE init_type 349 350 !-- Data type for the general information of chemistry emissions, do not dependent on the particular chemical species 351 TYPE chem_emis_att_type 352 353 !-DIMENSIONS 354 INTEGER(iwp) :: nspec !< number of chem species for which emission values are provided 355 INTEGER(iwp) :: ncat !< number of emission categories 356 INTEGER(iwp) :: nvoc !< number of VOCs components 357 INTEGER(iwp) :: npm !< number of PMs components 358 INTEGER(iwp) :: nnox=2 !< number of NOx components: NO and NO2 359 INTEGER(iwp) :: nsox=2 !< number of SOx components: SO and SO4 360 INTEGER(iwp) :: nhoursyear !< number of hours of a specific year in the HOURLY mode 361 ! of the default mode 362 INTEGER(iwp) :: nmonthdayhour !< number of month days and hours in the MDH mode 363 ! of the default mode 364 INTEGER(iwp) :: dt_emission !< Number of emissions timesteps for one year 365 ! in the pre-processed emissions case 366 !-- 1d emission input variables 367 CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: pm_name !< Names of PMs components 368 CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: cat_name !< Emission categories names 369 CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: species_name !< Names of emission chemical species 370 CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: voc_name !< Names of VOCs components 371 CHARACTER (LEN=25) :: units !< Units 372 373 INTEGER(iwp) :: i_hour !< indices for assigning the emission values at different timesteps 374 INTEGER(iwp),ALLOCATABLE, DIMENSION(:) :: cat_index !< Index of emission categories 375 INTEGER(iwp),ALLOCATABLE, DIMENSION(:) :: species_index !< Index of emission chem species 376 377 REAL(wp), DIMENSION(24) :: par_emis_time_factor !< time factors 378 ! for the parameterized mode: these are fixed for each hour 379 ! of a single day. 380 REAL(wp),ALLOCATABLE, DIMENSION(:) :: xm !< Molecular masses of emission chem species 381 382 !-- 2d emission input variables 383 REAL(wp),ALLOCATABLE, DIMENSION(:,:) :: hourly_emis_time_factor !< Time factors for HOURLY emissions (DEFAULT mode) 384 REAL(wp),ALLOCATABLE, DIMENSION(:,:) :: mdh_emis_time_factor !< Time factors for MDH emissions (DEFAULT mode) 385 REAL(wp),ALLOCATABLE, DIMENSION(:,:) :: nox_comp !< Composition of NO and NO2 386 REAL(wp),ALLOCATABLE, DIMENSION(:,:) :: sox_comp !< Composition of SO2 and SO4 387 REAL(wp),ALLOCATABLE, DIMENSION(:,:) :: voc_comp !< Composition of different VOC components (number not fixed) 388 389 !-- 3d emission input variables 390 REAL(wp),ALLOCATABLE, DIMENSION(:,:,:) :: pm_comp !< Composition of different PMs components (number not fixed) 391 392 END TYPE chem_emis_att_type 393 394 395 !-- Data type for the values of chemistry emissions ERUSSO 396 TYPE chem_emis_val_type 397 398 !REAL(wp),ALLOCATABLE, DIMENSION(:,:) :: stack_height !< stack height 399 400 !-- 3d emission input variables 401 REAL(wp),ALLOCATABLE, DIMENSION(:,:,:) :: default_emission_data !< Input Values emissions DEFAULT mode 402 403 !-- 4d emission input variables 404 REAL(wp),ALLOCATABLE, DIMENSION(:,:,:,:) :: preproc_emission_data !< Input Values emissions PRE-PROCESSED mode 405 406 END TYPE chem_emis_val_type 343 407 344 408 ! … … 546 610 TYPE(pars) :: water_pars_f !< input variable for water parameters 547 611 612 TYPE(chem_emis_att_type) :: chem_emis_att !< Input Information of Chemistry Emission Data from netcdf 613 TYPE(chem_emis_val_type), ALLOCATABLE, DIMENSION(:) :: chem_emis !< Input Chemistry Emission Data from netcdf 548 614 549 615 CHARACTER(LEN=3) :: char_lod = 'lod' !< name of level-of-detail attribute in NetCDF file … … 553 619 CHARACTER(LEN=100) :: input_file_static = 'PIDS_STATIC' !< Name of file which comprises static input data 554 620 CHARACTER(LEN=100) :: input_file_dynamic = 'PIDS_DYNAMIC' !< Name of file which comprises dynamic input data 621 CHARACTER(LEN=100) :: input_file_chem = 'PIDS_CHEM' !< Name of file which comprises chemistry input data 622 623 CHARACTER (LEN=25), ALLOCATABLE, DIMENSION(:) :: string_values !< output of string variables read from netcdf input files 624 625 INTEGER(iwp) :: id_emis !< NetCDF id of input file for chemistry emissions: TBD: It has to be removed 555 626 556 627 INTEGER(iwp) :: nc_stat !< return value of nf90 function call … … 558 629 LOGICAL :: input_pids_static = .FALSE. !< Flag indicating whether Palm-input-data-standard file containing static information exists 559 630 LOGICAL :: input_pids_dynamic = .FALSE. !< Flag indicating whether Palm-input-data-standard file containing dynamic information exists 631 LOGICAL :: input_pids_chem = .FALSE. !< Flag indicating whether Palm-input-data-standard file containing chemistry information exists 560 632 561 633 LOGICAL :: collective_read = .FALSE. !< Enable NetCDF collective read … … 581 653 MODULE PROCEDURE netcdf_data_input_check_static 582 654 END INTERFACE netcdf_data_input_check_static 655 656 INTERFACE netcdf_data_input_chemistry_data 657 MODULE PROCEDURE netcdf_data_input_chemistry_data 658 END INTERFACE netcdf_data_input_chemistry_data 583 659 584 660 INTERFACE netcdf_data_input_inquire_file … … 611 687 612 688 INTERFACE get_variable 689 MODULE PROCEDURE get_variable_string 613 690 MODULE PROCEDURE get_variable_1d_int 614 691 MODULE PROCEDURE get_variable_1d_real … … 619 696 MODULE PROCEDURE get_variable_3d_real 620 697 MODULE PROCEDURE get_variable_3d_real_dynamic 698 MODULE PROCEDURE get_variable_4d_to_3d_real 621 699 MODULE PROCEDURE get_variable_4d_real 700 MODULE PROCEDURE get_variable_5d_to_4d_real 622 701 END INTERFACE get_variable 623 702 … … 636 715 !-- Public variables 637 716 PUBLIC albedo_pars_f, albedo_type_f, basal_area_density_f, buildings_f, & 638 building_id_f, building_pars_f, building_type_f, init_3d, & 639 init_model, input_file_static, input_pids_static, & 717 building_id_f, building_pars_f, building_type_f, & 718 chem_emis, chem_emis_att, chem_emis_att_type, chem_emis_val_type, & 719 init_3d, init_model, input_file_static, input_pids_static, & 640 720 input_pids_dynamic, leaf_area_density_f, nest_offl, & 641 721 pavement_pars_f, pavement_subsurface_pars_f, pavement_type_f, & … … 648 728 !-- Public subroutines 649 729 PUBLIC netcdf_data_input_check_dynamic, netcdf_data_input_check_static, & 650 netcdf_data_input_ inquire_file, &730 netcdf_data_input_chemistry_data, netcdf_data_input_inquire_file, & 651 731 netcdf_data_input_init, netcdf_data_input_init_lsm, & 652 732 netcdf_data_input_init_3d, & … … 674 754 INQUIRE( FILE = TRIM( input_file_dynamic ) // TRIM( coupling_char ), & 675 755 EXIST = input_pids_dynamic ) 756 INQUIRE( FILE = TRIM( input_file_chem ) // TRIM( coupling_char ), & 757 EXIST = input_pids_chem ) 758 676 759 #endif 677 760 … … 756 839 757 840 END SUBROUTINE netcdf_data_input_init 841 842 !------------------------------------------------------------------------------! 843 ! Description: 844 ! ------------ 845 !> Reads Chemistry NETCDF Input data, such as emission values, emission species, etc. . 846 !------------------------------------------------------------------------------! 847 SUBROUTINE netcdf_data_input_chemistry_data(emt_att,emt) 848 849 USE chem_modules, & 850 ONLY: do_emis, mode_emis, time_fac_type, & 851 surface_csflux_name 852 853 USE control_parameters, & 854 ONLY: message_string 855 856 USE indices, & 857 ONLY: nz, nx, ny, nxl, nxr, nys, nyn, nzb, nzt 858 859 IMPLICIT NONE 860 861 TYPE(chem_emis_att_type), INTENT(INOUT) :: emt_att 862 TYPE(chem_emis_val_type), ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: emt 863 864 CHARACTER (LEN=80) :: units='' !< units of chemistry inputs 865 866 INTEGER(iwp) :: ispec !< index for number of emission species in input 867 868 INTEGER(iwp), ALLOCATABLE, DIMENSION(:) :: dum_var !< value of variable read from netcdf input 869 INTEGER(iwp) :: errno !< error number NF90_???? function 870 INTEGER(iwp) :: id_var !< variable id 871 ! INTEGER(iwp) :: id_emis !< NetCDF id of input file 872 INTEGER(iwp) :: num_vars !< number of variables in netcdf input file 873 INTEGER(iwp) :: len_dims,len_dims_2 !< Length of dimensions 874 875 INTEGER(iwp) :: max_string_length=25 !< Variable for the maximum length of a string 876 877 CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE :: var_names !< Name of Variables 878 879 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: dum_var_3d !< variable for storing temporary data of 3-dimensional 880 ! variables read from netcdf for chemistry emissions 881 882 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: dum_var_4d !< variable for storing temporary data of 5-dimensional 883 !< variables read from netcdf for chemistry emissions 884 !-- 885 !> Start the processing of the data 886 CALL location_message( 'starting allocation of chemistry emissions arrays', .FALSE. ) 887 888 !> Parameterized mode of the emissions 889 IF (TRIM(mode_emis)=="PARAMETERIZED" .OR. TRIM(mode_emis)=="parameterized") THEN 890 891 ispec=1 892 emt_att%nspec=0 893 894 !number of species 895 DO WHILE (TRIM( surface_csflux_name( ispec ) ) /= 'novalue' ) 896 897 emt_att%nspec=emt_att%nspec+1 898 ispec=ispec+1 899 900 ENDDO 901 902 !-- allocate emission values data type arrays 903 ALLOCATE(emt(emt_att%nspec)) 904 905 !-- Read EMISSION SPECIES NAMES 906 907 !Assign values 908 ALLOCATE(emt_att%species_name(emt_att%nspec)) 909 910 DO ispec=1,emt_att%nspec 911 emt_att%species_name(ispec) = TRIM(surface_csflux_name(ispec)) 912 ENDDO 913 914 !Assign Constant Values of time factors: 915 emt_att%par_emis_time_factor( : ) = (/ 0.01, 0.01, 0.01, 0.02, 0.03, 0.07, 0.09, 0.09, 0.05, 0.03, 0.03, 0.03, & 916 0.03, 0.03, 0.03, 0.04, 0.05, 0.09, 0.09, 0.09, 0.04, 0.02, 0.01, 0.01 /) 917 918 !> DEFAULT AND PRE-PROCESSED MODE 919 ELSE 920 921 #if defined ( __netcdf ) 922 IF ( .NOT. input_pids_chem ) RETURN 923 924 !-- Open file in read-only mode 925 CALL open_read_file( TRIM( input_file_chem ) // & 926 TRIM( coupling_char ), id_emis ) 927 !-- inquire number of variables 928 CALL inquire_num_variables( id_emis, num_vars ) 929 930 !-- Get General Dimension Lengths: only number of species and number of categories. 931 ! the other dimensions depend on the mode of the emissions or on the presence of specific components 932 !nspecies 933 CALL get_dimension_length( id_emis, emt_att%nspec, 'nspecies' ) 934 935 936 !-- Allocate emission values data type arrays 937 ALLOCATE(emt(1:emt_att%nspec)) 938 939 940 !-- Read EMISSION SPECIES NAMES 941 !Allocate Arrays 942 ALLOCATE(emt_att%species_name(emt_att%nspec)) 943 944 !Call get Variable 945 CALL get_variable( id_emis, 'emission_name', string_values, emt_att%nspec ) 946 emt_att%species_name=string_values 947 ! If allocated, Deallocate var_string, an array only used for reading-in strings 948 IF (ALLOCATED(string_values)) DEALLOCATE(string_values) 949 950 !-- Read EMISSION SPECIES INDEX 951 !Allocate Arrays 952 ALLOCATE(emt_att%species_index(emt_att%nspec)) 953 !Call get Variable 954 CALL get_variable( id_emis, 'emission_index', emt_att%species_index ) 955 956 957 !-- Now the routine has to distinguish between DEFAULT and PRE-PROCESSED chemistry emission modes 958 959 IF (TRIM(mode_emis)=="DEFAULT" .OR. TRIM(mode_emis)=="default") THEN 960 961 !number of categories 962 CALL get_dimension_length( id_emis, emt_att%ncat, 'ncat' ) 963 964 !-- Read EMISSION CATEGORIES INDEX 965 !Allocate Arrays 966 ALLOCATE(emt_att%cat_index(emt_att%ncat)) 967 !Call get Variable 968 CALL get_variable( id_emis, 'emission_cat_index', emt_att%cat_index ) 969 970 971 DO ispec=1,emt_att%nspec 972 !-- EMISSION_VOC_NAME (1-DIMENSIONAL) 973 IF (TRIM(emt_att%species_name(ispec))=="VOC" .OR. TRIM(emt_att%species_name(ispec))=="voc") THEN 974 !Allocate Array 975 CALL get_dimension_length( id_emis, emt_att%nvoc, 'nvoc' ) 976 ALLOCATE(emt_att%voc_name(1:emt_att%nvoc)) 977 !Read-in Variable 978 CALL get_variable( id_emis,"emission_voc_name",string_values, emt_att%nvoc) 979 emt_att%voc_name=string_values 980 IF (ALLOCATED(string_values)) DEALLOCATE(string_values) 981 982 !-- COMPOSITION VOC (2-DIMENSIONAL) 983 !Allocate Array 984 ALLOCATE(emt_att%voc_comp(1:emt_att%ncat,1:emt_att%nvoc)) 985 !Read-in Variable 986 ! CALL get_variable(id_emis,"composition_voc",emt%voc_comp,1,1,emt%ncat,emt%nvoc) 987 CALL get_variable(id_emis,"composition_voc",emt_att%voc_comp,1,emt_att%ncat,1,emt_att%nvoc) 988 ENDIF 989 990 !-- EMISSION_PM_NAME (1-DIMENSIONAL) 991 IF (TRIM(emt_att%species_name(ispec))=="PM" .OR. TRIM(emt_att%species_name(ispec))=="pm") THEN 992 CALL get_dimension_length( id_emis, emt_att%npm, 'npm' ) 993 ALLOCATE(emt_att%pm_name(1:emt_att%npm)) 994 !Read-in Variable 995 CALL get_variable( id_emis,"pm_name",string_values, emt_att%npm) 996 emt_att%pm_name=string_values 997 IF (ALLOCATED(string_values)) DEALLOCATE(string_values) 998 999 !-- COMPOSITION PM (3-DIMENSIONAL) 1000 !Allocate 1001 len_dims=3 !> number of PMs: PM1, PM2.5 and PM10 1002 ALLOCATE(emt_att%pm_comp(1:emt_att%ncat,1:emt_att%npm,1:len_dims)) 1003 !Read-in Variable 1004 CALL get_variable(id_emis,"composition_pm",emt_att%pm_comp,1,emt_att%ncat,1,emt_att%npm,1,len_dims) 1005 ENDIF 1006 1007 !-- COMPOSITION_NOX (2-DIMENSIONAL) 1008 IF (TRIM(emt_att%species_name(ispec))=="NOx" .OR. TRIM(emt_att%species_name(ispec))=="nox") THEN 1009 !Allocate array 1010 ALLOCATE(emt_att%nox_comp(1:emt_att%ncat,1:emt_att%nnox)) 1011 !Read-in Variable 1012 CALL get_variable(id_emis,"composition_nox",emt_att%nox_comp,1,emt_att%ncat,1,emt_att%nnox) 1013 ENDIF 1014 1015 !-- COMPOSITION-SOX (2-DIMENSIONAL) 1016 IF (TRIM(emt_att%species_name(ispec))=="SOx" .OR. TRIM(emt_att%species_name(ispec))=="sox") THEN 1017 ALLOCATE(emt_att%sox_comp(1:emt_att%ncat,1:emt_att%nsox)) 1018 !Read-in Variable 1019 CALL get_variable(id_emis,"composition_sox",emt_att%sox_comp,1,emt_att%ncat,1,emt_att%nsox) 1020 ENDIF 1021 ENDDO !>ispec 1022 1023 !-- For reading the emission time factors, the distinction between HOUR and MDH data is necessary 1024 1025 !-- EMISSION_TIME_SCALING_FACTORS 1026 !-- HOUR 1027 IF (TRIM(time_fac_type)=="HOUR" .OR. TRIM(time_fac_type)=="hour") THEN 1028 !-- Allocate Array 1029 CALL get_dimension_length( id_emis, emt_att%nhoursyear, 'nhoursyear' ) 1030 ALLOCATE(emt_att%hourly_emis_time_factor(1:emt_att%ncat,1:emt_att%nhoursyear)) 1031 !Read-in Variable 1032 CALL get_variable(id_emis,"emission_time_factors",emt_att%hourly_emis_time_factor,1,emt_att%ncat,1,emt_att%nhoursyear) 1033 1034 !-- MDH 1035 ELSE IF (TRIM(time_fac_type) == "MDH" .OR. TRIM(time_fac_type) == "mdh") THEN 1036 !-- Allocate Array 1037 CALL get_dimension_length( id_emis, emt_att%nmonthdayhour, 'nmonthdayhour' ) 1038 ALLOCATE(emt_att%mdh_emis_time_factor(1:emt_att%ncat,1:emt_att%nmonthdayhour)) 1039 !-- Read-in Variable 1040 CALL get_variable(id_emis,"emission_time_factors",emt_att%mdh_emis_time_factor,1,emt_att%ncat,1,emt_att%nmonthdayhour) 1041 1042 ELSE 1043 1044 message_string = 'We are in the DEFAULT chemistry emissions mode: ' // & 1045 ' !no time-factor type specified!' // & 1046 'Please, specify the value of time_fac_type:' // & 1047 ' either "MDH" or "HOUR"' 1048 CALL message( 'netcdf_data_input_chemistry_data', 'CM0200', 2, 2, 0, 6, 0 ) 1049 1050 1051 ENDIF 1052 1053 !-- Finally read-in the emission values and their units (DEFAULT mode) 1054 1055 DO ispec=1,emt_att%nspec 1056 1057 IF ( .NOT. ALLOCATED( emt(ispec)%default_emission_data ) ) ALLOCATE(emt(ispec)%default_emission_data(1:emt_att%ncat,1:ny+1,1:nx+1)) 1058 1059 ALLOCATE(dum_var_3d(1:emt_att%ncat,nys+1:nyn+1,nxl+1:nxr+1)) 1060 1061 CALL get_variable(id_emis,"emission_values",dum_var_3d,ispec,1,emt_att%ncat,nys,nyn,nxl,nxr) 1062 1063 emt(ispec)%default_emission_data(:,nys+1:nyn+1,nxl+1:nxr+1)=dum_var_3d(1:emt_att%ncat,nys+1:nyn+1,nxl+1:nxr+1) 1064 1065 DEALLOCATE (dum_var_3d) 1066 1067 ENDDO 1068 1069 !-- UNITS 1070 CALL get_attribute(id_emis,"units",emt_att%units,.FALSE.,"emission_values") 1071 1072 1073 !-- PRE-PROCESSED MODE -- 1074 1075 ELSE IF (TRIM(mode_emis)=="PRE-PROCESSED" .OR. TRIM(mode_emis)=="pre-processed") THEN 1076 !-- In the PRE-PROCESSED mode, only the VOC names, the VOC_composition, the emission values and their units remain to be read at this point 1077 1078 DO ispec=1,emt_att%nspec 1079 1080 !-- EMISSION_VOC_NAME (1-DIMENSIONAL) 1081 IF (TRIM(emt_att%species_name(ispec))=="VOC" .OR. TRIM(emt_att%species_name(ispec))=="voc") THEN 1082 !Allocate Array 1083 CALL get_dimension_length( id_emis, emt_att%nvoc, 'nvoc' ) 1084 ALLOCATE(emt_att%voc_name(1:emt_att%nvoc)) 1085 !Read-in Variable 1086 CALL get_variable( id_emis,"emission_voc_name",string_values, emt_att%nvoc) 1087 emt_att%voc_name=string_values 1088 IF (ALLOCATED(string_values)) DEALLOCATE(string_values) 1089 1090 !-- COMPOSITION VOC (2-DIMENSIONAL) 1091 !Allocate Array 1092 ALLOCATE(emt_att%voc_comp(1:emt_att%ncat,1:emt_att%nvoc)) 1093 !Read-in Variable 1094 CALL get_variable(id_emis,"composition_voc",emt_att%voc_comp,1,emt_att%ncat,1,emt_att%nvoc) 1095 ENDIF 1096 1097 ENDDO !> ispec 1098 1099 !-- EMISSION_VALUES (4-DIMENSIONAL) 1100 !Calculate temporal dimension length 1101 CALL get_dimension_length( id_emis, emt_att%dt_emission, 'time' ) 1102 1103 1104 DO ispec=1,emt_att%nspec 1105 1106 !Allocation for the entire domain has to be done only for the first processor between all the subdomains 1107 IF ( .NOT. ALLOCATED( emt(ispec)%preproc_emission_data ) ) ALLOCATE(emt(ispec)%preproc_emission_data(emt_att%dt_emission,1,1:ny+1,1:nx+1)) 1108 1109 !> allocate variable where to pass emission values read from netcdf 1110 ALLOCATE(dum_var_4d(1:emt_att%dt_emission,1,nys+1:nyn+1,nxl+1:nxr+1)) 1111 1112 !Read-in Variable 1113 CALL get_variable(id_emis,"emission_values",dum_var_4d,ispec,1,emt_att%dt_emission,1,1,nys,nyn,nxl,nxr) 1114 1115 1116 emt(ispec)%preproc_emission_data(:,1,nys+1:nyn+1,nxl+1:nxr+1)=dum_var_4d(1:emt_att%dt_emission,1,nys+1:nyn+1,nxl+1:nxr+1) 1117 1118 DEALLOCATE ( dum_var_4d ) 1119 1120 ENDDO 1121 1122 !-- UNITS 1123 CALL get_attribute(id_emis,"units",emt_att%units,.FALSE.,"emission_values") 1124 1125 ENDIF 1126 1127 CALL close_input_file( id_emis ) 1128 1129 #endif 1130 ENDIF 1131 1132 END SUBROUTINE netcdf_data_input_chemistry_data 758 1133 759 1134 !------------------------------------------------------------------------------! … … 4094 4469 ! Description: 4095 4470 ! ------------ 4471 !> Routine for reading-in a character string from the chem emissions netcdf input file. 4472 !------------------------------------------------------------------------------! 4473 4474 SUBROUTINE get_variable_string( id, variable_name, var_string, names_number) 4475 #if defined( __netcdf ) 4476 4477 USE indices 4478 USE pegrid 4479 4480 IMPLICIT NONE 4481 4482 CHARACTER (LEN=25), ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: var_string 4483 4484 CHARACTER(LEN=*) :: variable_name !> variable name 4485 4486 CHARACTER (LEN=1), ALLOCATABLE, DIMENSION(:,:) :: tmp_var_string !> variable to be read 4487 4488 4489 INTEGER(iwp), INTENT(IN) :: id !> file id 4490 4491 INTEGER(iwp), INTENT(IN) :: names_number !> number of names 4492 4493 INTEGER(iwp) :: id_var !> variable id 4494 4495 INTEGER(iwp) :: i,j !> index to go through the length of the dimensions 4496 4497 INTEGER(iwp) :: max_string_length=25 !> this is both the maximum length of a name, but also 4498 ! the number of the components of the first dimensions 4499 ! (rows) 4500 4501 4502 ALLOCATE(tmp_var_string(max_string_length,names_number)) 4503 4504 ALLOCATE(var_string(names_number)) 4505 4506 !-- Inquire variable id 4507 nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var ) 4508 4509 4510 !-- Get variable 4511 !-- Start cycle over the emission species 4512 DO i = 1, names_number 4513 !-- read the first letter of each component 4514 nc_stat = NF90_GET_VAR( id, id_var, var_string(i), start = (/ 1,i /), & 4515 count = (/ 1,1 /) ) 4516 CALL handle_error( 'get_variable_string', 701 ) 4517 4518 !-- Start cycle over charachters 4519 DO j = 1, max_string_length 4520 4521 !-- read the rest of the components of the name 4522 nc_stat = NF90_GET_VAR( id, id_var, tmp_var_string(j,i), start = (/ j,i /),& 4523 count = (/ 1,1 /) ) 4524 CALL handle_error( 'get_variable_string', 702 ) 4525 4526 IF ( iachar(tmp_var_string(j,i) ) == 0 ) THEN 4527 tmp_var_string(j,i)='' 4528 ENDIF 4529 4530 IF ( j>1 ) THEN 4531 !-- Concatenate first letter of the name and the others 4532 var_string(i)=TRIM(var_string(i)) // TRIM(tmp_var_string(j,i)) 4533 4534 ENDIF 4535 ENDDO 4536 ENDDO 4537 4538 #endif 4539 END SUBROUTINE get_variable_string 4540 4541 4542 !------------------------------------------------------------------------------! 4543 ! Description: 4544 ! ------------ 4096 4545 !> Reads a 1D integer variable from file. 4097 4546 !------------------------------------------------------------------------------! … … 4233 4682 #endif 4234 4683 ENDIF 4684 4685 4686 !Temporary solution for reading emission chemistry files: TBD: we should discuss whether remove it or not 4687 IF ( id==id_emis ) THEN 4688 4689 !-- Allocate temporary variable according to memory access on file. 4690 ALLOCATE( tmp(is:ie,js:je) ) 4691 4692 !-- Get variable 4693 nc_stat = NF90_GET_VAR( id, id_var, tmp, & 4694 start = (/ is, js /), & 4695 count = (/ ie-is+1 , je-js+1 /) ) 4696 4697 var=tmp 4698 4699 CALL handle_error( 'get_variable_2d_real', 530, variable_name ) !TBD: the error number shuld be changed, but since the solution is 4700 ! provisory, we give the same as below 4701 4702 DEALLOCATE( tmp ) 4703 4704 !> Original Subroutine part 4705 ELSE 4235 4706 ! 4236 4707 !-- Allocate temporary variable according to memory access on file. … … 4238 4709 ! 4239 4710 !-- Get variable 4240 nc_stat = NF90_GET_VAR( id, id_var, tmp, &4711 nc_stat = NF90_GET_VAR( id, id_var, tmp, & 4241 4712 start = (/ is+1, js+1 /), & 4242 4713 count = (/ ie-is + 1, je-js+1 /) ) 4243 4714 4244 CALL handle_error( 'get_variable_2d_real', 530, variable_name )4715 CALL handle_error( 'get_variable_2d_real', 530, variable_name ) 4245 4716 ! 4246 4717 !-- Resort data. Please note, dimension subscripts of var all start at 1. 4247 DO i = is, ie 4248 DO j = js, je 4249 var(j-js+1,i-is+1) = tmp(i,j) 4718 DO i = is, ie 4719 DO j = js, je 4720 var(j-js+1,i-is+1) = tmp(i,j) 4721 ENDDO 4250 4722 ENDDO 4251 ENDDO4252 4723 4253 DEALLOCATE( tmp ) 4254 4724 DEALLOCATE( tmp ) 4725 4726 ENDIF 4255 4727 #endif 4256 4728 END SUBROUTINE get_variable_2d_real … … 4647 5119 #endif 4648 5120 ENDIF 5121 5122 !Temporary solution for reading emission chemistry files: 5123 IF ( id==id_emis ) THEN 5124 5125 !-- Allocate temporary variable according to memory access on file. 5126 ALLOCATE( tmp(is:ie,js:je,k1s:k1e,k2s:k2e) ) 5127 5128 !-- Get variable 5129 nc_stat = NF90_GET_VAR( id, id_var, tmp, & 5130 start = (/ is, js, k1s+1, k2s+1 /), & 5131 count = (/ ie-is+1 , je-js+1, k1e-k1s+1, k2e-k2s+1 /) ) 5132 5133 var=tmp 5134 5135 CALL handle_error( 'get_variable_4d_real', 535, variable_name ) 5136 5137 DEALLOCATE( tmp ) 5138 5139 !> Original subroutine part 5140 ELSE 4649 5141 ! 4650 5142 !-- Allocate temporary variable according to memory access on file. … … 4652 5144 ! 4653 5145 !-- Get variable 4654 nc_stat = NF90_GET_VAR( id, id_var, tmp, &5146 nc_stat = NF90_GET_VAR( id, id_var, tmp, & 4655 5147 start = (/ is+1, js+1, k1s+1, k2s+1 /), & 4656 5148 count = (/ ie-is+1, je-js+1, & 4657 5149 k1e-k1s+1, k2e-k2s+1 /) ) 4658 5150 4659 CALL handle_error( 'get_variable_4d_real', 535, variable_name )5151 CALL handle_error( 'get_variable_4d_real', 535, variable_name ) 4660 5152 ! 4661 5153 !-- Resort data. Please note, dimension subscripts of var all start at 1. 4662 DO i = is, ie 4663 DO j = js, je 4664 DO k1 = k1s, k1e 4665 DO k2 = k2s, k2e 4666 var(k2-k2s+1,k1-k1s+1,j-js+1,i-is+1) = tmp(i,j,k1,k2) 5154 DO i = is, ie 5155 DO j = js, je 5156 DO k1 = k1s, k1e 5157 DO k2 = k2s, k2e 5158 var(k2-k2s+1,k1-k1s+1,j-js+1,i-is+1) = tmp(i,j,k1,k2) 5159 ENDDO 4667 5160 ENDDO 4668 5161 ENDDO 4669 5162 ENDDO 4670 ENDDO4671 5163 4672 DEALLOCATE( tmp ) 5164 DEALLOCATE( tmp ) 5165 ENDIF 4673 5166 #endif 4674 5167 END SUBROUTINE get_variable_4d_real 4675 5168 4676 5169 !------------------------------------------------------------------------------! 5170 ! Description: 5171 ! ------------ 5172 !> Reads a 4D float variable from file and store it to a 3-d variable. 5173 !------------------------------------------------------------------------------! 5174 SUBROUTINE get_variable_4d_to_3d_real( id, variable_name, var, ns, is, ie, js, je, & 5175 ks, ke ) 5176 5177 USE indices 5178 USE pegrid 5179 5180 IMPLICIT NONE 5181 5182 CHARACTER(LEN=*) :: variable_name !< variable name 5183 5184 INTEGER(iwp) :: ns !< start index for subdomain input along n dimension: ns coincides here with ne, since, we select only one value along the 1st dimension n 5185 5186 INTEGER(iwp) :: i !< index along x direction 5187 INTEGER(iwp) :: ie !< end index for subdomain input along x direction 5188 INTEGER(iwp) :: is !< start index for subdomain input along x direction 5189 INTEGER(iwp), INTENT(IN) :: id !< file id 5190 INTEGER(iwp) :: id_var !< variable id 5191 INTEGER(iwp) :: j !< index along y direction 5192 INTEGER(iwp) :: je !< end index for subdomain input along y direction 5193 INTEGER(iwp) :: js !< start index for subdomain input along y direction 5194 INTEGER(iwp) :: k !< index along any 4th dimension 5195 INTEGER(iwp) :: ke !< end index of 4th dimension 5196 INTEGER(iwp) :: ks !< start index of 4th dimension 5197 5198 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmp !< temporary variable to read data from file according 5199 !< to its reverse memory access 5200 5201 REAL(wp), DIMENSION(:,:,:), INTENT(INOUT) :: var !< variable where the read data have to be stored: one dimension is reduced in the process 5202 #if defined( __netcdf ) 5203 5204 ! 5205 !-- Inquire variable id 5206 nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var ) 5207 ! 5208 !-- Check for collective read-operation and set respective NetCDF flags if 5209 !-- required. 5210 IF ( collective_read ) THEN 5211 nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE) 5212 ENDIF 5213 5214 !Temporary solution for reading emission chemistry files: 5215 IF ( id==id_emis ) THEN 5216 5217 !-- Allocate temporary variable according to memory access on file. 5218 ALLOCATE( tmp(is:ie,js:je,ks:ke) ) 5219 5220 !-- Get variable 5221 nc_stat = NF90_GET_VAR( id, id_var, tmp(is:ie,js:je,ks:ke), & 5222 start = (/ ns, is, js+1, ks+1 /), & 5223 count = (/ 1, ie-is+1 , je-js+1, ke-ks+1 /) ) 5224 5225 var=tmp(:,:,:) 5226 5227 CALL handle_error( 'get_variable_4d_to_3d_real', 536, variable_name ) 5228 5229 DEALLOCATE( tmp ) 5230 5231 ELSE 5232 ! 5233 !-- Allocate temporary variable according to memory access on file. 5234 ALLOCATE( tmp(is:ie,js:je,ks:ke) ) 5235 ! 5236 !-- Get variable 5237 nc_stat = NF90_GET_VAR( id, id_var, tmp, & 5238 start = (/ ns+1, is+1, js+1, ks+1 /), & 5239 count = (/ 1, ie-is+1, je-js+1, ke-ks+1 /) ) 5240 5241 CALL handle_error( 'get_variable_4d_to_3d_real', 536, variable_name ) 5242 ! 5243 !-- Resort data. Please note, dimension subscripts of var all start at 1. 5244 DO i = is, ie 5245 DO j = js, je 5246 DO k = ks, ke 5247 var(k-ks+1,j-js+1,i-is+1) = tmp(i,j,k) 5248 ENDDO 5249 ENDDO 5250 ENDDO 5251 5252 DEALLOCATE( tmp ) 5253 5254 ENDIF 5255 #endif 5256 END SUBROUTINE get_variable_4d_to_3d_real 4677 5257 4678 5258 !------------------------------------------------------------------------------! … … 4755 5335 count = (/ count_1, count_2, count_3 /) ) 4756 5336 4757 CALL handle_error( 'get_variable_3d_real_dynamic', 53 6, variable_name )5337 CALL handle_error( 'get_variable_3d_real_dynamic', 537, variable_name ) 4758 5338 ! 4759 5339 !-- Resort data. Please note, dimension subscripts of var all start at 1. … … 4770 5350 END SUBROUTINE get_variable_3d_real_dynamic 4771 5351 5352 !------------------------------------------------------------------------------! 5353 ! Description: 5354 ! ------------ 5355 !> Reads a 5D float variable from file and store it to a 4-d variable. 5356 !------------------------------------------------------------------------------! 5357 SUBROUTINE get_variable_5d_to_4d_real( id, variable_name, var, & 5358 ns, ts, te, is, ie, js, je, ks, ke ) 5359 5360 USE indices 5361 USE pegrid 5362 5363 IMPLICIT NONE 5364 5365 CHARACTER(LEN=*) :: variable_name !< variable name 5366 5367 INTEGER(iwp) :: ns !< start index for subdomain input along n dimension: ns coincides here with ne, since, we select only one value along the 1st dimension n 5368 5369 INTEGER(iwp) :: t !< index along t direction 5370 INTEGER(iwp) :: te !< end index for subdomain input along t direction 5371 INTEGER(iwp) :: ts !< start index for subdomain input along t direction 5372 5373 INTEGER(iwp) :: i !< index along x direction 5374 INTEGER(iwp) :: ie !< end index for subdomain input along x direction 5375 INTEGER(iwp) :: is !< start index for subdomain input along x direction 5376 INTEGER(iwp), INTENT(IN) :: id !< file id 5377 INTEGER(iwp) :: id_var !< variable id 5378 INTEGER(iwp) :: j !< index along y direction 5379 INTEGER(iwp) :: je !< end index for subdomain input along y direction 5380 INTEGER(iwp) :: js !< start index for subdomain input along y direction 5381 INTEGER(iwp) :: k !< index along any 5th dimension 5382 INTEGER(iwp) :: ke !< end index of 5th dimension 5383 INTEGER(iwp) :: ks !< start index of 5th dimension 5384 5385 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tmp !< temporary variable to read data from file according 5386 ! to its reverse memory access 5387 REAL(wp), DIMENSION(:,:,:,:), INTENT(INOUT) :: var !< variable to be read 5388 #if defined( __netcdf ) 5389 ! 5390 !-- Inquire variable id 5391 nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var ) 5392 ! 5393 !-- Check for collective read-operation and set respective NetCDF flags if 5394 !-- required. 5395 IF ( collective_read ) THEN 5396 nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE) 5397 ENDIF 5398 5399 !Temporary solution for reading emission chemistry files: 5400 IF ( id==id_emis ) THEN 5401 5402 !-- Allocate temporary variable according to memory access on file. 5403 ALLOCATE( tmp(ts:te,1,js+1:je+1,ks+1:ke+1) ) 5404 5405 !-- Get variable 5406 nc_stat = NF90_GET_VAR( id, id_var, tmp(ts:te,1,js+1:je+1,ks+1:ke+1), & 5407 start = (/ ns, ts, 1, js+1, ks+1 /), & 5408 count = (/ 1, te-ts+1, 1, je-js+1, ke-ks+1 /) ) 5409 5410 var=tmp 5411 5412 CALL handle_error( 'get_variable_5d_to_4d_real', 538, variable_name ) 5413 5414 DEALLOCATE( tmp ) 5415 5416 !> Original Subroutine part 5417 ELSE 5418 ! 5419 !-- Allocate temporary variable according to memory access on file. 5420 ALLOCATE( tmp(ks:ke,js:je,is:is,ts:te) ) 5421 ! 5422 !-- Get variable 5423 nc_stat = NF90_GET_VAR( id, id_var, tmp, & 5424 start = (/ ks+1, js+1, is+1, ts+1, ns /), & 5425 count = (/ ke-ks+1, je-js+1, ie-is+1, te-ts+1, 1 /) ) 5426 5427 CALL handle_error( 'get_variable_5d_to_4d_real', 538, variable_name ) 5428 ! 5429 !-- Resort data. Please note, dimension subscripts of var all start at 1. 5430 5431 DO t = ts, te 5432 DO i = is, ie 5433 DO j = js, je 5434 DO k = ks, ke 5435 var(t-ts+1,i-is+1,j-js+1,k-ks+1) = tmp(k,j,i,t) 5436 ENDDO 5437 ENDDO 5438 ENDDO 5439 ENDDO 5440 5441 DEALLOCATE( tmp ) 5442 5443 ENDIF 5444 #endif 5445 END SUBROUTINE get_variable_5d_to_4d_real 4772 5446 4773 5447 … … 4789 5463 4790 5464 nc_stat = NF90_INQUIRE( id, NVARIABLES = num_vars ) 4791 CALL handle_error( 'inquire_num_variables', 53 7)5465 CALL handle_error( 'inquire_num_variables', 539 ) 4792 5466 4793 5467 #endif … … 4816 5490 ALLOCATE( varids(1:SIZE(var_names)) ) 4817 5491 nc_stat = NF90_INQ_VARIDS( id, NVARS = num_vars, VARIDS = varids ) 4818 CALL handle_error( 'inquire_variable_names', 5 38)5492 CALL handle_error( 'inquire_variable_names', 540 ) 4819 5493 4820 5494 DO i = 1, SIZE(var_names) 4821 5495 nc_stat = NF90_INQUIRE_VARIABLE( id, varids(i), NAME = var_names(i) ) 4822 CALL handle_error( 'inquire_variable_names', 5 38)5496 CALL handle_error( 'inquire_variable_names', 540 ) 4823 5497 ENDDO 4824 5498 -
palm/trunk/SOURCE/palm.f90
r3274 r3298 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Minor formatting (kanani) 28 ! - Added Call of date_and_time_init (Russo) 29 ! - Added Call of calc_date_and_time before call of init_3d where emissions 30 ! are initialized: 31 ! we have to know the time indices to initialize emission values (Russo) 32 ! - Added Call of netcdf_data_input_chemistry_data (Russo) 33 ! 34 ! 3274 2018-09-24 15:42:55Z knoop 27 35 ! Modularization of all bulk cloud physics code components 28 36 ! … … 241 249 USE arrays_3d 242 250 251 USE bulk_cloud_model_mod, & 252 ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert 253 254 USE chem_modules, & 255 ONLY: do_emis 256 243 257 USE chemistry_model_mod, & 244 ONLY: chem_ init245 246 USE chem_photolysis_mod, &247 ONLY: photolysis_init258 ONLY: chem_check_data_output_pr, chem_init 259 260 ! USE chem_photolysis_mod, & 261 ! ONLY: photolysis_init 248 262 249 263 USE control_parameters, & … … 259 273 ONLY: cpu_log, log_point, log_point_s, cpu_statistics 260 274 275 USE date_and_time_mod, & 276 ONLY: calc_date_and_time, init_date_and_time 277 261 278 USE indices, & 262 279 ONLY: nbgp … … 264 281 USE kinds 265 282 266 USE bulk_cloud_model_mod, &267 ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert268 269 283 USE multi_agent_system_mod, & 270 284 ONLY: agents_active, mas_last_actions 271 285 272 286 USE netcdf_data_input_mod, & 273 ONLY: netcdf_data_input_inquire_file, netcdf_data_input_init, & 287 ONLY: chem_emis, chem_emis_att, netcdf_data_input_chemistry_data, & 288 netcdf_data_input_inquire_file, netcdf_data_input_init, & 274 289 netcdf_data_input_surface_data, netcdf_data_input_topo 275 290 … … 410 425 !-- --> Needs to be moved!! What is the dependency about? 411 426 ! IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 412 IF ( air_chemistry ) THEN 427 IF ( air_chemistry .AND. do_emis ) THEN 428 429 IF ( do_emis ) THEN 430 CALL netcdf_data_input_chemistry_data(chem_emis_att,chem_emis) 431 ENDIF 432 413 433 CALL chem_init 414 CALL photolysis_init ! probably also required for restart 434 ! CALL photolysis_init ! probably also required for restart 435 436 CALL init_date_and_time !initialize the time of chemistry emissions 437 415 438 ENDIF 416 439 ! END IF … … 421 444 ! 422 445 !-- Initialize all necessary variables 446 CALL calc_date_and_time !this is required for chemistry emissions 447 423 448 CALL init_3d_model 424 449 -
palm/trunk/SOURCE/parin.f90
r3294 r3298 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Introduced reading of date_init to inipar.(Russo) 28 ! Add extra profiles for chemistry (basit) 29 ! 30 ! 3294 2018-10-01 02:37:10Z raasch 27 31 ! changes concerning modularization of ocean option 28 32 ! … … 421 425 ! ------------ 422 426 !> This subroutine reads variables controling the run from the NAMELIST files 427 !> 428 !> @todo: Revise max_pr_cs (profiles for chemistry) 423 429 !------------------------------------------------------------------------------! 424 430 SUBROUTINE parin … … 443 449 444 450 USE date_and_time_mod, & 445 ONLY: da y_of_year_init, time_utc_init451 ONLY: date_init, day_of_year_init, time_utc_init 446 452 447 453 USE dvrp_variables, & … … 548 554 cycle_mg, damp_level_1d, & 549 555 data_output_during_spinup, & 556 date_init, & 550 557 day_of_year_init, & 551 558 dissipation_1d, & … … 620 627 cycle_mg, damp_level_1d, & 621 628 data_output_during_spinup, & 629 date_init, & 622 630 day_of_year_init, & 623 631 dissipation_1d, & … … 1066 1074 ref_state(0:nz+1), sa_init(0:nz+1), ug(0:nz+1), & 1067 1075 u_init(0:nz+1), v_init(0:nz+1), vg(0:nz+1), & 1068 hom(0:nz+1,2,pr_palm+max_pr_user ,0:statistic_regions), &1069 hom_sum(0:nz+1,pr_palm+max_pr_user ,0:statistic_regions) )1076 hom(0:nz+1,2,pr_palm+max_pr_user+max_pr_cs,0:statistic_regions), & 1077 hom_sum(0:nz+1,pr_palm+max_pr_user+max_pr_cs,0:statistic_regions) ) 1070 1078 1071 1079 hom = 0.0_wp -
palm/trunk/SOURCE/prognostic_equations.f90
r3294 r3298 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Code added for decycling chemistry (basit) 28 ! 29 ! 3294 2018-10-01 02:37:10Z raasch 27 30 ! changes concerning modularization of ocean option 28 31 ! … … 335 338 336 339 USE chemistry_model_mod, & 337 ONLY: chem_integrate, chem_ prognostic_equations,&338 chem_species, nspec, nvar, spc_names340 ONLY: chem_integrate, chem_species, chem_prognostic_equations, & 341 nspec, nvar, spc_names, chem_boundary_conds 339 342 340 343 USE chem_modules, & … … 343 346 USE chem_photolysis_mod, & 344 347 ONLY: photolysis_control 348 349 USE chem_modules, & 350 ONLY: call_chem_at_all_substeps, chem_gasphase_on, cs_name 345 351 346 352 USE control_parameters, & … … 456 462 457 463 LOGICAL :: loop_start !< 458 INTEGER :: n, lsp !< lsp running index for chem spcs 464 INTEGER(iwp) :: n 465 INTEGER(iwp) :: lsp 466 INTEGER(iwp) :: lsp_usr !< lsp running index for chem spcs 459 467 460 468 … … 468 476 !-- concentrations of chemical species 469 477 IF ( air_chemistry ) THEN 478 lsp_usr = 1 470 479 ! 471 480 !-- If required, calculate photolysis frequencies - … … 492 501 !-- Loop over chemical species 493 502 CALL cpu_log( log_point_s(84), 'chemistry exch-horiz ', 'start' ) 494 DO n = 1, nspec 495 CALL exchange_horiz( chem_species(n)%conc, nbgp ) 503 DO lsp = 1, nspec 504 CALL exchange_horiz( chem_species(lsp)%conc, nbgp ) 505 lsp_usr = 1 506 DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' ) 507 IF ( TRIM(chem_species(lsp)%name) == TRIM(cs_name(lsp_usr)) ) THEN 508 509 CALL chem_boundary_conds( chem_species(lsp)%conc_p, & 510 chem_species(lsp)%conc_pr_init ) 511 512 ENDIF 513 lsp_usr = lsp_usr +1 514 ENDDO 515 496 516 ENDDO 497 517 CALL cpu_log( log_point_s(84), 'chemistry exch-horiz ', 'stop' ) -
palm/trunk/SOURCE/read_restart_data_mod.f90
r3294 r3298 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Added chemistry profiles for restart run (basit) 28 ! 29 ! 3294 2018-10-01 02:37:10Z raasch 27 30 ! changes concerning modularization of ocean option 28 31 ! … … 75 78 ! ------------ 76 79 !> Reads restart data from restart-file(s) (binary format). 80 !> 81 !> @todo: Revise max_pr_cs (profiles for chemistry) 77 82 !------------------------------------------------------------------------------! 78 83 MODULE read_restart_data_mod … … 80 85 81 86 USE control_parameters 87 88 USE chem_modules, & 89 ONLY: max_pr_cs 90 82 91 83 92 IMPLICIT NONE … … 266 275 v_init(0:nz+1), pt_init(0:nz+1), q_init(0:nz+1), & 267 276 ref_state(0:nz+1), s_init(0:nz+1), sa_init(0:nz+1), & 268 hom(0:nz+1,2,pr_palm+max_pr_user ,0:statistic_regions), &269 hom_sum(0:nz+1,pr_palm+max_pr_user ,0:statistic_regions) )277 hom(0:nz+1,2,pr_palm+max_pr_user+max_pr_cs,0:statistic_regions), & 278 hom_sum(0:nz+1,pr_palm+max_pr_user+max_pr_cs,0:statistic_regions) ) 270 279 ENDIF 271 280 -
palm/trunk/SOURCE/surface_mod.f90
r3294 r3298 26 26 ! ----------------- 27 27 ! $Id$ 28 ! - Minor formatting/clean-up (kanani) 29 ! - Removed initialization of cssws with surface_csflux values. Only the value 0 30 ! is now passed to all the chemical species of the selected mechanism. 31 ! The real inizialization is conducted in init_3d_model_mod.f90 (Russo) 32 ! - Added calls and variables for inizialization of chemistry emission 33 ! fluxes (Russo) 34 ! - Setting of cssws to zero moved (forkel) 35 ! - Set cssws to zero before setting values (forkel) 36 ! - Distinguish flux conversion factors for gases and PM (forkel) 37 ! - Some questions/todo's on conversion factors (forkel) 38 ! 39 ! 3294 2018-10-01 02:37:10Z raasch 28 40 ! changes concerning modularization of ocean option 29 41 ! … … 2071 2083 2072 2084 TYPE( surf_type ) :: surf !< respective surface type 2085 2073 2086 ! 2074 2087 !-- Store indices of respective surface element … … 2173 2186 2174 2187 DO lsp = 1, nvar 2175 IF ( air_chemistry ) surf%css(lsp,num_h) = 0.0_wp 2188 IF ( air_chemistry ) surf%css(lsp,num_h) = 0.0_wp 2189 ! 2190 !-- Ensure that fluxes of compounds which are not specified in 2191 !-- namelist are all zero --> kanani: revise 2192 IF ( air_chemistry ) surf%cssws(lsp,num_h) = 0.0_wp 2176 2193 ENDDO 2177 2194 ! … … 2220 2237 ELSE 2221 2238 surf%qsws(num_h) = wall_humidityflux(5) * & 2222 heatflux_input_conversion(k)2239 waterflux_input_conversion(k) 2223 2240 ENDIF 2224 2241 ENDIF … … 2227 2244 IF ( upward_facing ) THEN 2228 2245 IF ( constant_scalarflux ) THEN 2229 surf%ssws(num_h) = surface_scalarflux * rho_air_zw(k-1)2246 surf%ssws(num_h) = surface_scalarflux * rho_air_zw(k-1) 2230 2247 2231 2248 IF ( k-1 /= 0 ) & … … 3386 3403 nxr_on_file, nynf, nync, nyn_on_file, nysf, & 3387 3404 nysc, nys_on_file, found ) 3405 3388 3406 3389 3407 IMPLICIT NONE -
palm/trunk/SOURCE/time_integration.f90
r3294 r3298 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Formatting, clean-up, comments (kanani) 28 ! - Added CALL to chem_emissions_setup (Russo) 29 ! - Code added for decycling chemistry (basit) 30 ! 31 ! 3294 2018-10-01 02:37:10Z raasch 27 32 ! changes concerning modularization of ocean option 28 33 ! … … 380 385 ONLY: calc_mean_profile 381 386 387 USE chem_emissions_mod, & 388 ONLY: chem_emissions_setup 389 390 USE chem_modules, & 391 ONLY: bc_cs_t_val, call_chem_at_all_substeps, cs_name, & 392 constant_csflux, do_emis, nspec, nspec_out 393 382 394 USE chemistry_model_mod, & 383 ONLY: chem_emissions, chem_species 384 385 USE chem_modules, & 386 ONLY: nspec 395 ONLY: chem_boundary_conds, chem_species 387 396 388 397 USE control_parameters, & … … 426 435 ONLY: cpu_log, log_point, log_point_s 427 436 437 USE date_and_time_mod, & 438 ONLY: calc_date_and_time, hour_call_emis, hour_of_year 439 428 440 USE flight_mod, & 429 441 ONLY: flight_measurement … … 450 462 lsf_nesting_offline, lsf_nesting_offline_mass_conservation 451 463 452 USE netcdf_data_input_mod, &453 ONLY: nest_offl, netcdf_data_input_lsf454 455 464 USE multi_agent_system_mod, & 456 465 ONLY: agents_active, multi_agent_system 466 467 USE netcdf_data_input_mod, & 468 ONLY: chem_emis, chem_emis_att, nest_offl, netcdf_data_input_lsf 457 469 458 470 USE ocean_mod, & … … 525 537 IMPLICIT NONE 526 538 527 CHARACTER (LEN=9) :: time_to_string !< 528 529 INTEGER(iwp) :: n !< loop counter for chemistry species 539 CHARACTER (LEN=9) :: time_to_string !< 540 541 INTEGER(iwp) :: lsp !< 542 INTEGER(iwp) :: lsp_usr !< 543 INTEGER(iwp) :: n !< loop counter for chemistry species 530 544 531 545 REAL(wp) :: dt_3d_old !< temporary storage of timestep to be used for … … 611 625 bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1) 612 626 bc_q_t_val = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1) 627 IF ( air_chemistry ) THEN 628 DO lsp = 1, nspec 629 bc_cs_t_val = ( chem_species(lsp)%conc_pr_init(nzt+1) & 630 - chem_species(lsp)%conc_pr_init(nzt) ) & 631 / dzu(nzt+1) 632 ENDDO 633 ENDIF 613 634 ENDIF 614 635 ! … … 783 804 IF ( passive_scalar ) CALL exchange_horiz( s_p, nbgp ) 784 805 IF ( air_chemistry ) THEN 785 DO n = 1, nspec 786 CALL exchange_horiz( chem_species(n)%conc_p, nbgp ) 806 DO lsp = 1, nspec 807 CALL exchange_horiz( chem_species(lsp)%conc_p, nbgp ) 808 ! 809 !-- kanani: Push chem_boundary_conds after CALL boundary_conds 810 lsp_usr = 1 811 DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' ) 812 IF ( TRIM(chem_species(lsp)%name) == TRIM(cs_name(lsp_usr)) ) THEN 813 CALL chem_boundary_conds( chem_species(lsp)%conc_p, & 814 chem_species(lsp)%conc_pr_init ) 815 ENDIF 816 lsp_usr = lsp_usr + 1 817 ENDDO 787 818 ENDDO 788 819 ENDIF … … 1122 1153 ! 1123 1154 !-- If required, consider chemical emissions 1124 !-- (todo (FK): Implement hourly call of emissions, using time_utc from 1125 !-- data_and_time_mod.f90; 1126 !-- move the CALL to appropriate location) 1127 IF ( air_chemistry ) THEN 1128 CALL chem_emissions 1155 IF ( air_chemistry .AND. do_emis ) THEN 1156 ! 1157 !-- Update the time --> kanani: revise location of this CALL 1158 CALL calc_date_and_time 1159 ! 1160 !-- Call emission routine only once an hour 1161 IF (hour_of_year .GT. hour_call_emis ) THEN 1162 CALL chem_emissions_setup( chem_emis_att, chem_emis, nspec_out ) 1163 hour_call_emis = hour_of_year 1164 ENDIF 1129 1165 ENDIF 1130 1166 ! -
palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/src/lex.yy.c
r2696 r3298 10 10 #define YY_FLEX_MAJOR_VERSION 2 11 11 #define YY_FLEX_MINOR_VERSION 5 12 #define YY_FLEX_SUBMINOR_VERSION 3 712 #define YY_FLEX_SUBMINOR_VERSION 35 13 13 #if YY_FLEX_SUBMINOR_VERSION > 0 14 14 #define FLEX_BETA … … 143 143 /* Size of default input buffer. */ 144 144 #ifndef YY_BUF_SIZE 145 #ifdef __ia64__ 146 /* On IA-64, the buffer size is 16k, not 8k. 147 * Moreover, YY_BUF_SIZE is 2*YY_READ_BUF_SIZE in the general case. 148 * Ditto for the __ia64__ case accordingly. 149 */ 150 #define YY_BUF_SIZE 32768 151 #else 145 152 #define YY_BUF_SIZE 16384 153 #endif /* __ia64__ */ 146 154 #endif 147 155 … … 155 163 #endif 156 164 157 #ifndef YY_TYPEDEF_YY_SIZE_T 158 #define YY_TYPEDEF_YY_SIZE_T 159 typedef size_t yy_size_t; 160 #endif 161 162 extern yy_size_t yyleng; 165 extern int yyleng; 163 166 164 167 extern FILE *yyin, *yyout; … … 186 189 #define unput(c) yyunput( c, (yytext_ptr) ) 187 190 191 #ifndef YY_TYPEDEF_YY_SIZE_T 192 #define YY_TYPEDEF_YY_SIZE_T 193 typedef size_t yy_size_t; 194 #endif 195 188 196 #ifndef YY_STRUCT_YY_BUFFER_STATE 189 197 #define YY_STRUCT_YY_BUFFER_STATE … … 203 211 * characters. 204 212 */ 205 yy_size_t yy_n_chars;213 int yy_n_chars; 206 214 207 215 /* Whether we "own" the buffer - i.e., we know we created it, … … 273 281 /* yy_hold_char holds the character lost when yytext is formed. */ 274 282 static char yy_hold_char; 275 static yy_size_t yy_n_chars; /* number of characters read into yy_ch_buf */276 yy_size_t yyleng;283 static int yy_n_chars; /* number of characters read into yy_ch_buf */ 284 int yyleng; 277 285 278 286 /* Points to current character in buffer. */ … … 302 310 YY_BUFFER_STATE yy_scan_buffer (char *base,yy_size_t size ); 303 311 YY_BUFFER_STATE yy_scan_string (yyconst char *yy_str ); 304 YY_BUFFER_STATE yy_scan_bytes (yyconst char *bytes, yy_size_t len );312 YY_BUFFER_STATE yy_scan_bytes (yyconst char *bytes,int len ); 305 313 306 314 void *yyalloc (yy_size_t ); … … 801 809 return (x); \ 802 810 } 803 #line 8 04"lex.yy.c"811 #line 812 "lex.yy.c" 804 812 805 813 #define INITIAL 0 … … 864 872 void yyset_out (FILE * out_str ); 865 873 866 yy_size_t yyget_leng (void );874 int yyget_leng (void ); 867 875 868 876 char *yyget_text (void ); … … 906 914 /* Amount of stuff to slurp up with each read. */ 907 915 #ifndef YY_READ_BUF_SIZE 916 #ifdef __ia64__ 917 /* On IA-64, the buffer size is 16k, not 8k */ 918 #define YY_READ_BUF_SIZE 16384 919 #else 908 920 #define YY_READ_BUF_SIZE 8192 921 #endif /* __ia64__ */ 909 922 #endif 910 923 … … 1009 1022 #line 120 "scan.l" 1010 1023 1011 #line 10 12"lex.yy.c"1024 #line 1025 "lex.yy.c" 1012 1025 1013 1026 if ( !(yy_init) ) … … 1526 1539 ECHO; 1527 1540 YY_BREAK 1528 #line 15 29"lex.yy.c"1541 #line 1542 "lex.yy.c" 1529 1542 1530 1543 case YY_END_OF_BUFFER: … … 1710 1723 else 1711 1724 { 1712 yy_size_t num_to_read =1725 int num_to_read = 1713 1726 YY_CURRENT_BUFFER_LVALUE->yy_buf_size - number_to_move - 1; 1714 1727 … … 1717 1730 1718 1731 /* just a shorter name for the current buffer */ 1719 YY_BUFFER_STATE b = YY_CURRENT_BUFFER _LVALUE;1732 YY_BUFFER_STATE b = YY_CURRENT_BUFFER; 1720 1733 1721 1734 int yy_c_buf_p_offset = … … 1724 1737 if ( b->yy_is_our_buffer ) 1725 1738 { 1726 yy_size_t new_size = b->yy_buf_size * 2;1739 int new_size = b->yy_buf_size * 2; 1727 1740 1728 1741 if ( new_size <= 0 ) … … 1755 1768 /* Read in more data. */ 1756 1769 YY_INPUT( (&YY_CURRENT_BUFFER_LVALUE->yy_ch_buf[number_to_move]), 1757 (yy_n_chars), num_to_read );1770 (yy_n_chars), (size_t) num_to_read ); 1758 1771 1759 1772 YY_CURRENT_BUFFER_LVALUE->yy_n_chars = (yy_n_chars); … … 1850 1863 yy_is_jam = (yy_current_state == 198); 1851 1864 1852 1865 return yy_is_jam ? 0 : yy_current_state; 1853 1866 } 1854 1867 … … 1865 1878 { /* need to shift things up to make room */ 1866 1879 /* +2 for EOB chars. */ 1867 register yy_size_t number_to_move = (yy_n_chars) + 2;1880 register int number_to_move = (yy_n_chars) + 2; 1868 1881 register char *dest = &YY_CURRENT_BUFFER_LVALUE->yy_ch_buf[ 1869 1882 YY_CURRENT_BUFFER_LVALUE->yy_buf_size + 2]; … … 1914 1927 else 1915 1928 { /* need more input */ 1916 yy_size_t offset = (yy_c_buf_p) - (yytext_ptr);1929 int offset = (yy_c_buf_p) - (yytext_ptr); 1917 1930 ++(yy_c_buf_p); 1918 1931 … … 2074 2087 } 2075 2088 2089 #ifndef __cplusplus 2090 extern int isatty (int ); 2091 #endif /* __cplusplus */ 2092 2076 2093 /* Initializes or reinitializes a buffer. 2077 2094 * This function is sometimes called more than once on the same buffer, … … 2186 2203 static void yyensure_buffer_stack (void) 2187 2204 { 2188 yy_size_t num_to_alloc;2205 int num_to_alloc; 2189 2206 2190 2207 if (!(yy_buffer_stack)) { … … 2283 2300 * @return the newly allocated buffer state object. 2284 2301 */ 2285 YY_BUFFER_STATE yy_scan_bytes (yyconst char * yybytes, yy_size_t _yybytes_len )2302 YY_BUFFER_STATE yy_scan_bytes (yyconst char * yybytes, int _yybytes_len ) 2286 2303 { 2287 2304 YY_BUFFER_STATE b; … … 2370 2387 * 2371 2388 */ 2372 yy_size_t yyget_leng (void)2389 int yyget_leng (void) 2373 2390 { 2374 2391 return yyleng; -
palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/bin/kpp4palm.ksh
r2718 r3298 21 21 # Copyright 2017-2018 Leibniz Universitaet Hannover 22 22 #------------------------------------------------------------------------------# 23 # Nov. 2016: Initial Version of KPP chemistry convertor adapted for PALM 24 # by Klaus Ketelsen 25 # 26 # This code is a modified version of KP4 (Jöckel, P., Kerkweg, A., Pozzer, A., 27 # Sander, R., Tost, H., Riede, H., Baumgaertner, A., Gromov, S., and Kern, B., 28 # 2010: Development cycle 2 of the Modular Earth Submodel System (MESSy2), 29 # Geosci. Model Dev., 3, 717-752, https://doi.org/10.5194/gmd-3-717-2010). 30 # KP4 is part of the Modular Earth Submodel System (MESSy), which is is 31 # available under the GNU General Public License (GPL). 32 # 33 #------------------------------------------------------------------------------# 23 34 # 24 35 # Current revisions: … … 29 40 # ----------------- 30 41 # $Id$ 42 # forkel 25. September 2018: Added cat for $MECH to pass mechanism name to kpp4palm 43 # ketelsen 18. September 2018: Added cat for '#INLINE F90_GLOBAL' 44 # (moved here from mechanisms/def_MECH/chem_gasphase.kpp 45 # 46 # forkel: 14. September 2018: WCOPY removed 47 # ketelsen: July 2018: Adaptations for vektor mode 48 # forkel June 2018: re-established original case of subroutine names 49 # forkel May 2018: additional copying of chem_gasphase_mod.f90 into $DEFDIR 50 # forkel 20.04.2018: removed wlamch and wlamch_add from $KPP_SUBROUTINE_LIST 51 # (epsilon(one) is used now) 52 # forkel March 2017 53 # Re-introduced relative path for KPP_HOME 54 # Subroutine list adapted to lowercase subroutine names 55 # Added arr2, removed update_sun and k_3rd from subroutine list 56 # Renamed output file to chem_gasphase_mod 57 # Renamed this file from kp4/ksh to kpp4kpp.ksh 58 # changed location of def_mechanism directories to gasphase_preproc/mechanisms 59 # 60 # 61 # 2718 2018-01-02 08:49:38Z maronga 31 62 # Initial revision 32 63 # 33 64 # 65 ########################################################################## 34 66 # 35 67 # … … 63 95 # Default 64 96 97 MECH=smog 65 98 OUTDIR=`pwd`/../../../SOURCE 66 99 OUTFILE=chem_gasphase_mod … … 70 103 VLEN=1 71 104 KEEP="NO" 72 DE_INDEX="NO" 73 DE_INDEX_FAST="NO" 105 UPDT="NO" 106 DE_INDEX=0 107 DE_INDEX_FAST="YES" 74 108 75 109 export KPP_SOLVER=Rosenbrock … … 77 111 # get Command line option 78 112 79 echo xxxxxxxxxx 80 while getopts :d:ifkp:o:s:v:w: c # get options 113 while getopts :m:i:fkup:o:s:vl:w: c # get options 81 114 do case $c in 82 d) DEFDIR=$OPTARG;; # directory of definition files83 84 i) DE_INDEX= "YES";;# if set, deindexing115 m) MECH=$OPTARG;; # mechanism 116 117 i) DE_INDEX=$OPTARG;; # if set, deindexing 85 118 86 119 f) DE_INDEX_FAST="YES";; # if set, fast deindexing … … 92 125 p) PREFIX=$OPTARG;; # Name Prefix 93 126 94 s) KPP_SOLVER=$OPTARG;;