Changeset 24
- Timestamp:
- May 23, 2014, 11:48:41 AM (10 years ago)
- Location:
- trunk
- Files:
-
- 1 deleted
- 47 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/options/COMMAND
r4 r24 1 +++++++++++++ HEADER +++++++++++++++++ 2 +++++++++++++ HEADER +++++++++++++++++ 3 +++++++++++++ HEADER +++++++++++++++++ 4 +++++++++++++ HEADER +++++++++++++++++ 5 +++++++++++++ HEADER +++++++++++++++++ 6 +++++++++++++ HEADER +++++++++++++++++ 7 +++++++++++++ HEADER +++++++++++++++++ 8 -1 9 20070930 180000 10 20070931 0 11 3600 12 3600 13 900 14 9999999 15 900 SYNC 16 -5.0 CTL 17 4 IFINE 18 5 IOUT 19 0 IPOUT 20 1 LSUBGRID 21 1 LCONVECTION 22 1 LAGESPECTRA 23 0 IPIN 24 1 IOFR 25 0 IFLUX 26 0 MDOMAINFILL 27 1 IND_SOURCE 28 2 IND_RECEPTOR 29 0 MQUASILAG 30 1 NESTED_OUTPUT 31 2 LINIT_COND INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT 1 ******************************************************************************** 2 * * 3 * Input file for the Lagrangian particle dispersion model FLEXPART * 4 * Please select your options * 5 * * 6 ******************************************************************************** 7 8 1. __ 3X, I2 9 1 10 LDIRECT 1 FOR FORWARD SIMULATION, -1 FOR BACKWARD SIMULATION 11 12 2. ________ ______ 3X, I8, 1X, I6 13 20111210 000000 14 YYYYMMDD HHMISS BEGINNING DATE OF SIMULATION 15 16 3. ________ ______ 3X, I8, 1X, I6 17 20111210 120000 18 YYYYMMDD HHMISS ENDING DATE OF SIMULATION 19 20 4. _____ 3X, I5 21 10800 22 SSSSS OUTPUT EVERY SSSSS SECONDS 23 24 5. _____ 3X, I5 25 10800 26 SSSSS TIME AVERAGE OF OUTPUT (IN SSSSS SECONDS) 27 28 6. _____ 3X, I5 29 900 30 SSSSS SAMPLING RATE OF OUTPUT (IN SSSSS SECONDS) 31 32 7. _________ 3X, I9 33 999999999 34 SSSSSSSSS TIME CONSTANT FOR PARTICLE SPLITTING (IN SECONDS) 35 36 8. _____ 3X, I5 37 900 38 SSSSS SYNCHRONISATION INTERVAL OF FLEXPART (IN SECONDS) 39 40 9. ---.-- 4X, F6.4 41 -5.0 42 CTL FACTOR, BY WHICH TIME STEP MUST BE SMALLER THAN TL 43 44 10. --- 4X, I3 45 4 46 IFINE DECREASE OF TIME STEP FOR VERTICAL MOTION BY FACTOR IFINE 47 48 11. - 4X, I1 49 3 50 IOUT 1 CONCENTRATION (RESIDENCE TIME FOR BACKWARD RUNS) OUTPUT, 2 MIXING RATIO OUTPUT, 3 BOTH,4 PLUME TRAJECT., 5=1+4 51 52 12. - 4X, I1 53 0 54 IPOUT PARTICLE DUMP: 0 NO, 1 EVERY OUTPUT INTERVAL, 2 ONLY AT END 55 56 13. _ 4X, I1 57 1 58 LSUBGRID SUBGRID TERRAIN EFFECT PARAMETERIZATION: 1 YES, 0 NO 59 60 14. _ 4X, I1 61 1 62 LCONVECTION CONVECTION: 1 YES, 0 NO 63 64 15. _ 4X, I1 65 0 66 LAGESPECTRA AGE SPECTRA: 1 YES, 0 NO 67 68 16. _ 4X, I1 69 0 70 IPIN CONTINUE SIMULATION WITH DUMPED PARTICLE DATA: 1 YES, 0 NO 71 72 17. _ 73 0 4X,I1 74 IOFR IOUTPUTFOREACHREL CREATE AN OUPUT FILE FOR EACH RELEASE LOCATION: 1 YES, 0 NO 75 76 18. _ 4X, I1 77 0 78 IFLUX CALCULATE FLUXES: 1 YES, 0 NO 79 80 19. _ 4X, I1 81 0 82 MDOMAINFILL DOMAIN-FILLING TRAJECTORY OPTION: 1 YES, 0 NO, 2 STRAT. O3 TRACER 83 84 20. _ 4X, I1 85 1 86 IND_SOURCE 1=MASS UNIT , 2=MASS MIXING RATIO UNIT 87 88 21. _ 4X, I1 89 1 90 IND_RECEPTOR 1=MASS UNIT , 2=MASS MIXING RATIO UNIT 91 92 22. _ 4X, I1 93 0 94 MQUASILAG QUASILAGRANGIAN MODE TO TRACK INDIVIDUAL PARTICLES: 1 YES, 0 NO 95 96 23. _ 4X, I1 97 0 98 NESTED_OUTPUT SHALL NESTED OUTPUT BE USED? 1 YES, 0 NO 99 100 24. _ 4X, I1 101 2 102 LINIT_COND INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT 103 104 25. _ 4X, I1 105 0 106 SURF_ONLY IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER 107 108 109 1. Simulation direction, 1 for forward, -1 for backward in time 110 (consult Seibert and Frank, 2004 for backward runs) 111 112 2. Beginning date and time of simulation. Must be given in format 113 YYYYMMDD HHMISS, where YYYY is YEAR, MM is MONTH, DD is DAY, HH is HOUR, 114 MI is MINUTE and SS is SECOND. Current version utilizes UTC. 115 116 3. Ending date and time of simulation. Same format as 3. 117 118 4. Average concentrations are calculated every SSSSS seconds. 119 120 5. The average concentrations are time averages of SSSSS seconds 121 duration. If SSSSS is 0, instantaneous concentrations are outputted. 122 123 6. The concentrations are sampled every SSSSS seconds to calculate the time 124 average concentration. This period must be shorter than the averaging time. 125 126 7. Time constant for particle splitting. Particles are split into two 127 after SSSSS seconds, 2xSSSSS seconds, 4xSSSSS seconds, and so on. 128 129 8. All processes are synchronized with this time interval (lsynctime). 130 Therefore, all other time constants must be multiples of this value. 131 Output interval and time average of output must be at least twice lsynctime. 132 133 9. CTL must be >1 for time steps shorter than the Lagrangian time scale 134 If CTL<0, a purely random walk simulation is done 135 136 10.IFINE=Reduction factor for time step used for vertical wind 137 138 11.IOUT determines how the output shall be made: concentration 139 (ng/m3, Bq/m3), mixing ratio (pptv), or both, or plume trajectory mode, 140 or concentration + plume trajectory mode. 141 In plume trajectory mode, output is in the form of average trajectories. 142 143 12.IPOUT determines whether particle positions are outputted (in addition 144 to the gridded concentrations or mixing ratios) or not. 145 0=no output, 1 output every output interval, 2 only at end of the 146 simulation 147 148 13.Switch on/off subgridscale terrain parameterization (increase of 149 mixing heights due to subgridscale orographic variations) 150 151 14.Switch on/off the convection parameterization 152 153 15.Switch on/off the calculation of age spectra: if yes, the file AGECLASSES 154 must be available 155 156 16. If IPIN=1, a file "partposit_end" from a previous run must be available in 157 the output directory. Particle positions are read in and previous simulation 158 is continued. If IPIN=0, no particles from a previous run are used 159 160 17. IF IOUTPUTFOREACHRELEASE is set to 1, one output field for each location 161 in the RLEASE file is created. For backward calculation this should be 162 set to 1. For forward calculation both possibilities are applicable. 163 164 18. If IFLUX is set to 1, fluxes of each species through each of the output 165 boxes are calculated. Six fluxes, corresponding to northward, southward, 166 eastward, westward, upward and downward are calculated for each grid cell of 167 the output grid. The control surfaces are placed in the middle of each 168 output grid cell. If IFLUX is set to 0, no fluxes are determined. 169 170 19. If MDOMAINFILL is set to 1, the first box specified in file RELEASES is used 171 as the domain where domain-filling trajectory calculations are to be done. 172 Particles are initialized uniformly distributed (according to the air mass 173 distribution) in that domain at the beginning of the simulation, and are 174 created at the boundaries throughout the simulation period. 175 176 20. IND_SOURCE switches between different units for concentrations at the source 177 NOTE that in backward simulations the release of computational particles 178 takes place at the "receptor" and the sampling of particles at the "source". 179 1=mass units (for bwd-runs = concentration) 180 2=mass mixing ratio units 181 21. IND_RECEPTOR switches between different units for concentrations at the receptor 182 1=mass units (concentrations) 183 2=mass mixing ratio units 184 185 22. MQUASILAG indicates whether particles shall be numbered consecutively (1) or 186 with their release location number (0). The first option allows tracking of 187 individual particles using the partposit output files 188 189 23. NESTED_OUTPUT decides whether model output shall be made also for a nested 190 output field (normally with higher resolution) 191 192 24. LINIT_COND determines whether, for backward runs only, the sensitivity to initial 193 conditions shall be calculated and written to output files 194 0=no output, 1 or 2 determines in which units the initial conditions are provided. 195 196 25. SURF_ONLY: When set to 1, concentration/emission sensitivity is written out only 197 for the surface layer; useful for instance when only footprint emission sensitivity is needed 198 but initial conditions are needed on a full 3-D grid -
trunk/options/COMMAND.alternative
r4 r24 30 30 0 NESTED_OUTPUT SHALL NESTED OUTPUT BE USED? YES, 0 NO 31 31 2 LINIT_COND INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT 32 0 SURF_ONLY IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER 32 33 33 34 … … 117 118 conditions shall be calculated and written to output files 118 119 0=no output, 1 or 2 determines in which units the initial conditions are provided. 120 121 25. SURF_ONLY: When set to 1, concentration/emission sensitivity is written out only 122 for the surface layer; useful for instance when only footprint emission sensitivity is needed 123 but initial conditions are needed on a full 3-D grid -
trunk/options/COMMAND.reference
r4 r24 11 11 12 12 2. ________ ______ 3X, I8, 1X, I6 13 20 040720 00000013 20110310 000000 14 14 YYYYMMDD HHMISS BEGINNING DATE OF SIMULATION 15 15 16 16 3. ________ ______ 3X, I8, 1X, I6 17 20 04072112000017 20110310 120000 18 18 YYYYMMDD HHMISS ENDING DATE OF SIMULATION 19 19 … … 101 101 2 102 102 LINIT_COND INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT 103 104 25. _ 4X, I1 105 0 106 SURF_ONLY IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER 103 107 104 108 … … 189 193 conditions shall be calculated and written to output files 190 194 0=no output, 1 or 2 determines in which units the initial conditions are provided. 195 196 25. SURF_ONLY: When set to 1, concentration/emission sensitivity is written out only 197 for the surface layer; useful for instance when only footprint emission sensitivity is needed 198 but initial conditions are needed on a full 3-D grid -
trunk/options/RELEASES
r4 r24 1 +++++++++++++++++ HEADER ++++++++++++++++++++ 2 +++++++++++++++++ HEADER ++++++++++++++++++++ 3 +++++++++++++++++ HEADER ++++++++++++++++++++ 4 +++++++++++++++++ HEADER ++++++++++++++++++++ 5 +++++++++++++++++ HEADER ++++++++++++++++++++ 6 +++++++++++++++++ HEADER ++++++++++++++++++++ 7 +++++++++++++++++ HEADER ++++++++++++++++++++ 8 +++++++++++++++++ HEADER ++++++++++++++++++++ 9 +++++++++++++++++ HEADER ++++++++++++++++++++ 10 +++++++++++++++++ HEADER ++++++++++++++++++++ 11 +++++++++++++++++ HEADER ++++++++++++++++++++ 12 5 13 31 14 31 15 31 16 31 17 31 1 ************************************************************************* 2 * * 3 * * 4 * * 5 * Input file for the Lagrangian particle dispersion model FLEXPART * 6 * Please select your options * 7 * * 8 * * 9 * * 10 ************************************************************************* 11 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 12 1 13 ___ i3 Total number of species emitted 18 14 19 20040101 080000 20 20040102 080000 21 8.2500 22 58.3833 23 8.2500 24 58.3833 25 1 26 0.000 27 100.000 28 4000 29 1 30 1 31 1 32 1 33 1 34 rel_bj 15 3 16 ___ i3 Index of species in file SPECIES 35 17 36 20040108 080000 37 20040109 080000 38 8.2500 39 58.3833 40 8.2500 41 58.3833 42 1 43 0.000 44 100.000 45 4000 46 1 47 1 48 1 49 1 50 1 51 rel_bj 18 ========================================================================= 19 20111210 000001 20 ________ ______ i8,1x,i6 Beginning date and time of release 52 21 53 20040115 080000 54 20040116 080000 55 8.2500 56 58.3833 57 8.2500 58 58.3833 59 1 60 0.000 61 100.000 62 4000 63 1 64 1 65 1 66 1 67 1 68 rel_bj 22 20111210 090000 23 ________ ______ i8,1x,i6 Ending date and time of release 69 24 70 20040122 080000 71 20040123 080000 72 8.2500 73 58.3833 74 8.2500 75 58.3833 76 1 77 0.000 78 100.000 79 4000 80 1 81 1 82 1 83 1 84 1 85 rel_bj 25 9.4048 26 ____.____ f9.4 Longitude [DEG] of lower left corner 86 27 87 20040129 080000 88 20040130 080000 89 8.2500 90 58.3833 91 8.2500 92 58.3833 93 1 94 0.000 95 100.000 96 4000 97 1 98 1 99 1 100 1 101 1 102 rel_bj 28 48.5060 29 ____.____ f9.4 Latitude [DEG] of lower left corner 103 30 104 20040205 080000 105 20040206 080000 106 8.2500 107 58.3833 108 8.2500 109 58.3833 110 1 111 0.000 112 100.000 113 4000 114 1 115 1 116 1 117 1 118 1 119 rel_bj 31 9.5067 32 ____.____ f9.4 Longitude [DEG] of upper right corner 120 33 121 20040212 080000 122 20040213 080000 123 8.2500 124 58.3833 125 8.2500 126 58.3833 127 1 128 0.000 129 100.000 130 4000 131 1 132 1 133 1 134 1 135 1 136 rel_bj 34 48.5158 35 ____.____ f9.4 Latitude [DEG] of upper right corner 137 36 138 20040219 080000 139 20040220 080000 140 8.2500 141 58.3833 142 8.2500 143 58.3833 144 1 145 0.000 146 100.000 147 4000 148 1 149 1 150 1 151 1 152 1 153 rel_bj 37 2 38 _________ i9 1 for m above ground, 2 for m above sea level 154 39 155 20040226 080000 156 20040227 080000 157 8.2500 158 58.3833 159 8.2500 160 58.3833 161 1 162 0.000 163 100.000 164 4000 165 1 166 1 167 1 168 1 169 1 170 rel_bj 40 0.00 41 _____.___ f10.3 Lower z-level (in m agl or m asl) 42 43 10.00 44 _____.___ f10.3 Upper z-level (in m agl or m asl) 45 46 20000 47 _________ i9 Total number of particles to be released 171 48 172 20040304 080000 173 20040305 080000 174 8.2500 175 58.3833 176 8.2500 177 58.3833 178 1 179 0.000 180 100.000 181 4000 182 1 183 1 184 1 185 1 186 1 187 rel_bj 49 1.0000E00 50 _.____E__ e9.4 Total mass emitted 188 51 189 20040311 080000 190 20040312 080000 191 8.2500 192 58.3833 193 8.2500 194 58.3833 195 1 196 0.000 197 100.000 198 4000 199 1 200 1 201 1 202 1 203 1 204 rel_bj 205 206 20040318 080000 207 20040319 080000 208 8.2500 209 58.3833 210 8.2500 211 58.3833 212 1 213 0.000 214 100.000 215 4000 216 1 217 1 218 1 219 1 220 1 221 rel_bj 222 223 20040325 080000 224 20040326 080000 225 8.2500 226 58.3833 227 8.2500 228 58.3833 229 1 230 0.000 231 100.000 232 4000 233 1 234 1 235 1 236 1 237 1 238 rel_bj 239 240 20040401 080000 241 20040402 080000 242 8.2500 243 58.3833 244 8.2500 245 58.3833 246 1 247 0.000 248 100.000 249 4000 250 1 251 1 252 1 253 1 254 1 255 rel_bj 256 257 20040408 080000 258 20040409 080000 259 8.2500 260 58.3833 261 8.2500 262 58.3833 263 1 264 0.000 265 100.000 266 4000 267 1 268 1 269 1 270 1 271 1 272 rel_bj 273 274 20040415 080000 275 20040416 080000 276 8.2500 277 58.3833 278 8.2500 279 58.3833 280 1 281 0.000 282 100.000 283 4000 284 1 285 1 286 1 287 1 288 1 289 rel_bj 290 291 20040422 080000 292 20040423 080000 293 8.2500 294 58.3833 295 8.2500 296 58.3833 297 1 298 0.000 299 100.000 300 4000 301 1 302 1 303 1 304 1 305 1 306 rel_bj 307 308 20040429 080000 309 20040430 080000 310 8.2500 311 58.3833 312 8.2500 313 58.3833 314 1 315 0.000 316 100.000 317 4000 318 1 319 1 320 1 321 1 322 1 323 rel_bj 324 325 20040609 080000 326 20040610 080000 327 8.2500 328 58.3833 329 8.2500 330 58.3833 331 1 332 0.000 333 100.000 334 4000 335 1 336 1 337 1 338 1 339 1 340 rel_bj 341 342 20040610 080000 343 20040611 080000 344 8.2500 345 58.3833 346 8.2500 347 58.3833 348 1 349 0.000 350 100.000 351 4000 352 1 353 1 354 1 355 1 356 1 357 rel_bj 358 359 20040617 080000 360 20040618 080000 361 8.2500 362 58.3833 363 8.2500 364 58.3833 365 1 366 0.000 367 100.000 368 4000 369 1 370 1 371 1 372 1 373 1 374 rel_bj 375 376 20040624 080000 377 20040625 080000 378 8.2500 379 58.3833 380 8.2500 381 58.3833 382 1 383 0.000 384 100.000 385 4000 386 1 387 1 388 1 389 1 390 1 391 rel_bj 392 393 20040701 080000 394 20040702 080000 395 8.2500 396 58.3833 397 8.2500 398 58.3833 399 1 400 0.000 401 100.000 402 4000 403 1 404 1 405 1 406 1 407 1 408 rel_bj 409 410 20040708 080000 411 20040709 080000 412 8.2500 413 58.3833 414 8.2500 415 58.3833 416 1 417 0.000 418 100.000 419 4000 420 1 421 1 422 1 423 1 424 1 425 rel_bj 426 427 20040715 080000 428 20040716 080000 429 8.2500 430 58.3833 431 8.2500 432 58.3833 433 1 434 0.000 435 100.000 436 4000 437 1 438 1 439 1 440 1 441 1 442 rel_bj 443 444 20040722 080000 445 20040723 080000 446 8.2500 447 58.3833 448 8.2500 449 58.3833 450 1 451 0.000 452 100.000 453 4000 454 1 455 1 456 1 457 1 458 1 459 rel_bj 460 461 20040729 080000 462 20040730 080000 463 8.2500 464 58.3833 465 8.2500 466 58.3833 467 1 468 0.000 469 100.000 470 4000 471 1 472 1 473 1 474 1 475 1 476 rel_bj 477 478 20040805 080000 479 20040806 080000 480 8.2500 481 58.3833 482 8.2500 483 58.3833 484 1 485 0.000 486 100.000 487 4000 488 1 489 1 490 1 491 1 492 1 493 rel_bj 494 495 20040812 080000 496 20040813 080000 497 8.2500 498 58.3833 499 8.2500 500 58.3833 501 1 502 0.000 503 100.000 504 4000 505 1 506 1 507 1 508 1 509 1 510 rel_bj 511 512 20040819 080000 513 20040820 080000 514 8.2500 515 58.3833 516 8.2500 517 58.3833 518 1 519 0.000 520 100.000 521 4000 522 1 523 1 524 1 525 1 526 1 527 rel_bj 528 529 20040826 080000 530 20040827 080000 531 8.2500 532 58.3833 533 8.2500 534 58.3833 535 1 536 0.000 537 100.000 538 4000 539 1 540 1 541 1 542 1 543 1 544 rel_bj 545 546 20040902 080000 547 20040903 080000 548 8.2500 549 58.3833 550 8.2500 551 58.3833 552 1 553 0.000 554 100.000 555 4000 556 1 557 1 558 1 559 1 560 1 561 rel_bj 562 563 20040909 080000 564 20040910 080000 565 8.2500 566 58.3833 567 8.2500 568 58.3833 569 1 570 0.000 571 100.000 572 4000 573 1 574 1 575 1 576 1 577 1 578 rel_bj 579 580 20040916 080000 581 20040917 080000 582 8.2500 583 58.3833 584 8.2500 585 58.3833 586 1 587 0.000 588 100.000 589 4000 590 1 591 1 592 1 593 1 594 1 595 rel_bj 596 597 20040923 080000 598 20040924 080000 599 8.2500 600 58.3833 601 8.2500 602 58.3833 603 1 604 0.000 605 100.000 606 4000 607 1 608 1 609 1 610 1 611 1 612 rel_bj 613 614 20040930 080000 615 20041001 080000 616 8.2500 617 58.3833 618 8.2500 619 58.3833 620 1 621 0.000 622 100.000 623 4000 624 1 625 1 626 1 627 1 628 1 629 rel_bj 630 631 20041007 080000 632 20041008 080000 633 8.2500 634 58.3833 635 8.2500 636 58.3833 637 1 638 0.000 639 100.000 640 4000 641 1 642 1 643 1 644 1 645 1 646 rel_bj 647 648 20041014 080000 649 20041015 080000 650 8.2500 651 58.3833 652 8.2500 653 58.3833 654 1 655 0.000 656 100.000 657 4000 658 1 659 1 660 1 661 1 662 1 663 rel_bj 664 665 20041021 080000 666 20041022 080000 667 8.2500 668 58.3833 669 8.2500 670 58.3833 671 1 672 0.000 673 100.000 674 4000 675 1 676 1 677 1 678 1 679 1 680 rel_bj 681 682 20041028 080000 683 20041029 080000 684 8.2500 685 58.3833 686 8.2500 687 58.3833 688 1 689 0.000 690 100.000 691 4000 692 1 693 1 694 1 695 1 696 1 697 rel_bj 698 699 20041104 080000 700 20041105 080000 701 8.2500 702 58.3833 703 8.2500 704 58.3833 705 1 706 0.000 707 100.000 708 4000 709 1 710 1 711 1 712 1 713 1 714 rel_bj 715 716 20041111 080000 717 20041112 080000 718 8.2500 719 58.3833 720 8.2500 721 58.3833 722 1 723 0.000 724 100.000 725 4000 726 1 727 1 728 1 729 1 730 1 731 rel_bj 732 733 20041118 080000 734 20041119 080000 735 8.2500 736 58.3833 737 8.2500 738 58.3833 739 1 740 0.000 741 100.000 742 4000 743 1 744 1 745 1 746 1 747 1 748 rel_bj 749 750 20041124 080000 751 20041125 080000 752 8.2500 753 58.3833 754 8.2500 755 58.3833 756 1 757 0.000 758 100.000 759 4000 760 1 761 1 762 1 763 1 764 1 765 rel_bj 766 767 20041125 080000 768 20041126 080000 769 8.2500 770 58.3833 771 8.2500 772 58.3833 773 1 774 0.000 775 100.000 776 4000 777 1 778 1 779 1 780 1 781 1 782 rel_bj 783 784 20041208 080000 785 20041209 080000 786 8.2500 787 58.3833 788 8.2500 789 58.3833 790 1 791 0.000 792 100.000 793 4000 794 1 795 1 796 1 797 1 798 1 799 rel_bj 800 801 20041209 080000 802 20041210 080000 803 8.2500 804 58.3833 805 8.2500 806 58.3833 807 1 808 0.000 809 100.000 810 4000 811 1 812 1 813 1 814 1 815 1 816 rel_bj 817 818 20041216 080000 819 20041217 080000 820 8.2500 821 58.3833 822 8.2500 823 58.3833 824 1 825 0.000 826 100.000 827 4000 828 1 829 1 830 1 831 1 832 1 833 rel_bj 834 835 20041223 080000 836 20041224 080000 837 8.2500 838 58.3833 839 8.2500 840 58.3833 841 1 842 0.000 843 100.000 844 4000 845 1 846 1 847 1 848 1 849 1 850 rel_bj 851 852 20041230 080000 853 20041231 080000 854 8.2500 855 58.3833 856 8.2500 857 58.3833 858 1 859 0.000 860 100.000 861 4000 862 1 863 1 864 1 865 1 866 1 867 rel_bj 868 52 RELEASE_TEST1 53 ________________________________________ character*40 comment 54 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -
trunk/options/SPECIES/SPECIES_001
r4 r24 9 9 -9.9E-09 Wet deposition - A 10 10 Wet deposition - B 11 -9.9E-09 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 -9.9 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 -9.9 In-cloud scavenging - Ci (S_i=Ci/cl) 14 -9.9 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 -9.9 Dry deposition (gases) - D 12 16 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_002
r4 r24 9 9 -9.9E-09 Wet deposition - A 10 10 Wet deposition - B 11 -9.9E-09 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 -9.9 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 -9.9 In-cloud scavenging - Ci (S_i=Ci/cl) 14 -9.9 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 1.5 Dry deposition (gases) - D 12 16 1.0e-02 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_003
r4 r24 9 9 8.0E-06 Wet deposition - A 10 10 0.62 Wet deposition - B 11 2.0E-7 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 0.36 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 0.1 In-cloud scavenging - Ci (S_i=Ci/cl) 14 1. In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 1.2 Dry deposition (gases) - D 12 16 3.0e-03 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_004
r4 r24 9 9 1.0E-05 Wet deposition - A 10 10 0.62 Wet deposition - B 11 2.0E-7 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 0.36 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 0.1 In-cloud scavenging - Ci (S_i=Ci/cl) 14 1. In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 1.6 Dry deposition (gases) - D 12 16 1.0e-02 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_005
r4 r24 9 9 5.0E-05 Wet deposition - A 10 10 0.62 Wet deposition - B 11 2.0E-7 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 0.36 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 0.1 In-cloud scavenging - Ci (S_i=Ci/cl) 14 1. In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 1.9 Dry deposition (gases) - D 12 16 1.0e+14 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_006
r4 r24 9 9 -9.9E-09 Wet deposition - A 10 10 Wet deposition - B 11 -9.9E-09 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 -9.9 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 -9.9 In-cloud scavenging - Ci (S_i=Ci/cl) 14 -9.9 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 1.6 Dry deposition (gases) - D 12 16 1.0e+05 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_007
r4 r24 9 9 1.0E-04 Wet deposition - A 10 10 0.62 Wet deposition - B 11 2.0E-7 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 0.36 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 0.1 In-cloud scavenging - Ci (S_i=Ci/cl) 14 1. In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 1.4 Dry deposition (gases) - D 12 16 1.0e+05 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_008
r4 r24 9 9 -9.9E-09 Wet deposition - A 10 10 0.62 Wet deposition - B 11 -9.9E-09 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 -9.9 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 -9.9 In-cloud scavenging - Ci (S_i=Ci/cl) 14 -9.9 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 2.0 Dry deposition (gases) - D 12 16 1.0e+05 Dry deposition (gases) - Henrys const. … … 39 43 16 1.463 1.162 40 44 17 1.426 1.141 41 18 1.326 1.117 45 18 1.326 1.117 42 46 19 1.225 1.109 43 47 20 1.082 1.095 -
trunk/options/SPECIES/SPECIES_009
r4 r24 9 9 -9.9E-09 Wet deposition - A 10 10 Wet deposition - B 11 -9.9E-09 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 -9.9 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 -9.9 In-cloud scavenging - Ci (S_i=Ci/cl) 14 -9.9 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 1.3 Dry deposition (gases) - D 12 16 6.0e+03 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_010
r4 r24 9 9 -9.9E-09 Wet deposition - A 10 10 Wet deposition - B 11 -9.9E-09 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 -9.9 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 -9.9 In-cloud scavenging - Ci (S_i=Ci/cl) 14 -9.9 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 2.6 Dry deposition (gases) - D 12 16 3.6e+00 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_011
r4 r24 9 9 9.9E-05 Wet deposition - A 10 10 0.62 Wet deposition - B 11 2.0E-7 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 0.36 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 0.1 In-cloud scavenging - Ci (S_i=Ci/cl) 14 1. In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 1.1 Dry deposition (gases) - D 12 16 2.0e+14 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_012
r4 r24 9 9 5.0E-06 Wet deposition - A 10 10 0.62 Wet deposition - B 11 2.0E-7 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 0.36 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 0.1 In-cloud scavenging - Ci (S_i=Ci/cl) 14 1. In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 -9.9 Dry deposition (gases) - D 12 16 1.0E-09 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_013
r4 r24 9 9 5.0E-06 Wet deposition - A 10 10 0.62 Wet deposition - B 11 2.0E-7 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 0.36 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 0.1 In-cloud scavenging - Ci (S_i=Ci/cl) 14 1. In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 -9.9 Dry deposition (gases) - D 12 16 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_014
r4 r24 9 9 8.0E-05 Wet deposition - A 10 10 0.62 Wet deposition - B 11 2.0E-7 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 0.36 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 0.1 In-cloud scavenging - Ci (S_i=Ci/cl) 14 1. In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 2.7 Dry deposition (gases) - D 12 16 1.0e+05 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_015
r4 r24 9 9 1.0E-04 Wet deposition - A 10 10 0.80 Wet deposition - B 11 2.0E-7 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 0.36 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 0.1 In-cloud scavenging - Ci (S_i=Ci/cl) 14 1. In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 -9.9 Dry deposition (gases) - D 12 16 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_016
r4 r24 9 9 1.0E-04 Wet deposition - A 10 10 0.80 Wet deposition - B 11 2.0E-7 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 0.36 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 0.1 In-cloud scavenging - Ci (S_i=Ci/cl) 14 1. In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 -9.9 Dry deposition (gases) - D 12 16 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_017
r4 r24 9 9 1.0E-04 Wet deposition - A 10 10 0.80 Wet deposition - B 11 2.0E-7 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 0.36 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 0.1 In-cloud scavenging - Ci (S_i=Ci/cl) 14 1. In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 -9.9 Dry deposition (gases) - D 12 16 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_018
r4 r24 9 9 1.0E-04 Wet deposition - A 10 10 0.80 Wet deposition - B 11 2.0E-7 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 0.36 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 0.1 In-cloud scavenging - Ci (S_i=Ci/cl) 14 1. In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 -9.9 Dry deposition (gases) - D 12 16 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_019
r4 r24 9 9 -9.9E-09 Wet deposition - A 10 10 Wet deposition - B 11 -9.9E-09 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 -9.9 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 -9.9 In-cloud scavenging - Ci (S_i=Ci/cl) 14 -9.9 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 -9.9 Dry deposition (gases) - D 12 16 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_020
r4 r24 9 9 1.0E-04 Wet deposition - A 10 10 0.80 Wet deposition - B 11 2.0E-7 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 0.36 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 0.1 In-cloud scavenging - Ci (S_i=Ci/cl) 14 1. In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 -9.9 Dry deposition (gases) - D 12 16 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_021
r4 r24 6 6 **************************************************************************** 7 7 Xe-133 Tracer name 8 198720.0 Species half life8 453168.0 Species half life 9 9 -9.9E-09 Wet deposition - A 10 10 Wet deposition - B 11 -9.9E-09 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 -9.9 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 -9.9 In-cloud scavenging - Ci (S_i=Ci/cl) 14 -9.9 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 -9.9 Dry deposition (gases) - D 12 16 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_022
r4 r24 9 9 -9.9E-09 Wet deposition - A 10 10 Wet deposition - B 11 -9.9E-09 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 -9.9 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 -9.9 In-cloud scavenging - Ci (S_i=Ci/cl) 14 -9.9 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 -9.9 Dry deposition (gases) - D 12 16 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_024
r4 r24 9 9 -9.9E-09 Wet deposition - A 10 10 Wet deposition - B 11 -9.9E-09 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 -9.9 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 -9.9 In-cloud scavenging - Ci (S_i=Ci/cl) 14 -9.9 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 -9.9 Dry deposition (gases) - D 12 16 Dry deposition (gases) - Henrys const. -
trunk/options/SPECIES/SPECIES_025
r4 r24 9 9 2.0E-05 Wet deposition - A 10 10 0.8 Wet deposition - B 11 2.0E-7 In-cloud scavenging - Ai (cl=Ai*prec**Bi) 12 0.36 In-cloud scavenging - Bi (cl=Ai*prec**Bi) 13 0.1 In-cloud scavenging - Ci (S_i=Ci/cl) 14 1. In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di) 11 15 -9.9 Dry deposition (gases) - D 12 16 Dry deposition (gases) - Henrys const. -
trunk/src/FLEXPART.f90
r20 r24 61 61 62 62 ! 63 flexversion='Version 9. 1.8 (2013-12-08)'63 flexversion='Version 9.2 beta (2014-05-23)' 64 64 !verbosity=0 65 65 ! Read the pathnames where input/output files are stored … … 111 111 ! Print the GPL License statement 112 112 !******************************************************* 113 print*,'Welcome to FLEXPART ', trim(flexversion)113 print*,'Welcome to FLEXPART ', trim(flexversion) 114 114 print*,'FLEXPART is free software released under the GNU Genera'// & 115 115 'l Public License.' -
trunk/src/calcpv.f90
r4 r24 21 21 22 22 subroutine calcpv(n,uuh,vvh,pvh) 23 ! 23 ! i i i o 24 24 !***************************************************************************** 25 25 ! * … … 49 49 integer :: nlck 50 50 real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f 51 real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt ,ppmk52 real :: pvavr,ppml( nuvzmax)51 real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt 52 real :: pvavr,ppml(0:nxmax-1,0:nymax-1,nuvzmax),ppmk(0:nxmax-1,0:nymax-1,nuvzmax) 53 53 real :: thup,thdn 54 54 real,parameter :: eps=1.e-5, p0=101325 … … 62 62 ! Loop over entire grid 63 63 !********************** 64 do kl=1,nuvz 65 do jy=0,nymin1 66 do ix=0,nxmin1 67 ppml(ix,jy,kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n) 68 enddo 69 enddo 70 enddo 71 ppmk=(100000./ppml)**kappa 72 64 73 do jy=0,nymin1 65 74 if (sglobal.and.jy.eq.0) goto 10 … … 106 115 end if 107 116 jux=jumpx 108 ! Precalculate pressure values for efficiency109 do kl=1,nuvz110 ppml(kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n)111 end do112 117 ! 113 118 ! Loop over the vertical … … 115 120 116 121 do kl=1,nuvz 117 ppmk=akz(kl)+bkz(kl)*ps(ix,jy,1,n) 118 theta=tth(ix,jy,kl,n)*(100000./ppmk)**kappa 122 theta=tth(ix,jy,kl,n)*ppmk(ix,jy,kl) 119 123 klvrp=kl+1 120 124 klvrm=kl-1 … … 125 129 if (klvrp.gt.nuvz) klvrp=nuvz 126 130 if (klvrm.lt.1) klvrm=1 127 ppmk=akz(klvrp)+bkz(klvrp)*ps(ix,jy,1,n) 128 thetap=tth(ix,jy,klvrp,n)*(100000./ppmk)**kappa 129 ppmk=akz(klvrm)+bkz(klvrm)*ps(ix,jy,1,n) 130 thetam=tth(ix,jy,klvrm,n)*(100000./ppmk)**kappa 131 dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm)) 131 thetap=tth(ix,jy,klvrp,n)*ppmk(ix,jy,klvrp) 132 thetam=tth(ix,jy,klvrm,n)*ppmk(ix,jy,klvrm) 133 dthetadp=(thetap-thetam)/(ppml(ix,jy,klvrp)-ppml(ix,jy,klvrm)) 132 134 133 135 ! Compute vertical position at pot. temperature surface on subgrid … … 156 158 kch=kch+1 157 159 k=kup 158 ppmk=akz(k)+bkz(k)*ps(ivr,jy,1,n) 159 thdn=tth(ivr,jy,k,n)*(100000./ppmk)**kappa 160 ppmk=akz(k+1)+bkz(k+1)*ps(ivr,jy,1,n) 161 thup=tth(ivr,jy,k+1,n)*(100000./ppmk)**kappa 162 163 164 if (((thdn.ge.theta).and.(thup.le.theta)).or. & 165 ((thdn.le.theta).and.(thup.ge.theta))) then 166 dt1=abs(theta-thdn) 167 dt2=abs(theta-thup) 168 dt=dt1+dt2 169 if (dt.lt.eps) then ! Avoid division by zero error 170 dt1=0.5 ! G.W., 10.4.1996 171 dt2=0.5 172 dt=1.0 173 endif 174 vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt 175 goto 20 160 thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k) 161 thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1) 162 163 164 if (((thdn.ge.theta).and.(thup.le.theta)).or. & 165 ((thdn.le.theta).and.(thup.ge.theta))) then 166 dt1=abs(theta-thdn) 167 dt2=abs(theta-thup) 168 dt=dt1+dt2 169 if (dt.lt.eps) then ! Avoid division by zero error 170 dt1=0.5 ! G.W., 10.4.1996 171 dt2=0.5 172 dt=1.0 176 173 endif 174 vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt 175 goto 20 176 endif 177 177 41 continue 178 178 ! Downward branch … … 181 181 kch=kch+1 182 182 k=kdn 183 ppmk=akz(k)+bkz(k)*ps(ivr,jy,1,n) 184 thdn=tth(ivr,jy,k,n)*(100000./ppmk)**kappa 185 ppmk=akz(k+1)+bkz(k+1)*ps(ivr,jy,1,n) 186 thup=tth(ivr,jy,k+1,n)*(100000./ppmk)**kappa 187 188 if (((thdn.ge.theta).and.(thup.le.theta)).or. & 189 ((thdn.le.theta).and.(thup.ge.theta))) then 190 dt1=abs(theta-thdn) 191 dt2=abs(theta-thup) 192 dt=dt1+dt2 193 if (dt.lt.eps) then ! Avoid division by zero error 194 dt1=0.5 ! G.W., 10.4.1996 195 dt2=0.5 196 dt=1.0 197 endif 198 vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt 199 goto 20 183 thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k) 184 thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1) 185 186 if (((thdn.ge.theta).and.(thup.le.theta)).or. & 187 ((thdn.le.theta).and.(thup.ge.theta))) then 188 dt1=abs(theta-thdn) 189 dt2=abs(theta-thup) 190 dt=dt1+dt2 191 if (dt.lt.eps) then ! Avoid division by zero error 192 dt1=0.5 ! G.W., 10.4.1996 193 dt2=0.5 194 dt=1.0 200 195 endif 201 goto 40 196 vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt 197 goto 20 198 endif 199 goto 40 202 200 ! This section used when no values were found 203 201 21 continue … … 236 234 kch=kch+1 237 235 k=kup 238 ppmk=akz(k)+bkz(k)*ps(ix,j,1,n) 239 thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa 240 ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n) 241 thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa 242 if (((thdn.ge.theta).and.(thup.le.theta)).or. & 243 ((thdn.le.theta).and.(thup.ge.theta))) then 244 dt1=abs(theta-thdn) 245 dt2=abs(theta-thup) 246 dt=dt1+dt2 247 if (dt.lt.eps) then ! Avoid division by zero error 248 dt1=0.5 ! G.W., 10.4.1996 249 dt2=0.5 250 dt=1.0 251 endif 252 uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt 253 goto 50 236 thdn=tth(ix,j,k,n)*ppmk(ix,j,k) 237 thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1) 238 if (((thdn.ge.theta).and.(thup.le.theta)).or. & 239 ((thdn.le.theta).and.(thup.ge.theta))) then 240 dt1=abs(theta-thdn) 241 dt2=abs(theta-thup) 242 dt=dt1+dt2 243 if (dt.lt.eps) then ! Avoid division by zero error 244 dt1=0.5 ! G.W., 10.4.1996 245 dt2=0.5 246 dt=1.0 254 247 endif 248 uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt 249 goto 50 250 endif 255 251 71 continue 256 252 ! Downward branch … … 259 255 kch=kch+1 260 256 k=kdn 261 ppmk=akz(k)+bkz(k)*ps(ix,j,1,n) 262 thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa 263 ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n) 264 thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa 265 if (((thdn.ge.theta).and.(thup.le.theta)).or. & 266 ((thdn.le.theta).and.(thup.ge.theta))) then 267 dt1=abs(theta-thdn) 268 dt2=abs(theta-thup) 269 dt=dt1+dt2 270 if (dt.lt.eps) then ! Avoid division by zero error 271 dt1=0.5 ! G.W., 10.4.1996 272 dt2=0.5 273 dt=1.0 274 endif 275 uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt 276 goto 50 257 thdn=tth(ix,j,k,n)*ppmk(ix,j,k) 258 thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1) 259 if (((thdn.ge.theta).and.(thup.le.theta)).or. & 260 ((thdn.le.theta).and.(thup.ge.theta))) then 261 dt1=abs(theta-thdn) 262 dt2=abs(theta-thup) 263 dt=dt1+dt2 264 if (dt.lt.eps) then ! Avoid division by zero error 265 dt1=0.5 ! G.W., 10.4.1996 266 dt2=0.5 267 dt=1.0 277 268 endif 278 goto 70 269 uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt 270 goto 50 271 endif 272 goto 70 279 273 ! This section used when no values were found 280 274 51 continue 281 275 ! Must use uu at current level and lat. juy becomes smaller by 1 282 uy(jj)=uuh(ix,jy,kl)283 juy=juy-1276 uy(jj)=uuh(ix,jy,kl) 277 juy=juy-1 284 278 ! Otherwise OK 285 279 50 continue -
trunk/src/calcpv_nests.f90
r4 r24 21 21 22 22 subroutine calcpv_nests(l,n,uuhn,vvhn,pvhn) 23 ! 23 ! i i i i o 24 24 !***************************************************************************** 25 25 ! * … … 48 48 integer :: nlck,l 49 49 real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f 50 real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt ,ppmk51 real :: ppml( nuvzmax)50 real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt 51 real :: ppml(0:nxmaxn-1,0:nymaxn-1,nuvzmax),ppmk(0:nxmaxn-1,0:nymaxn-1,nuvzmax) 52 52 real :: thup,thdn 53 53 real,parameter :: eps=1.e-5,p0=101325 … … 61 61 ! Loop over entire grid 62 62 !********************** 63 do kl=1,nuvz 64 do jy=0,nyn(l)-1 65 do ix=0,nxn(l)-1 66 ppml(ix,jy,kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l) 67 enddo 68 enddo 69 enddo 70 ppmk=(100000./ppml)**kappa 71 63 72 do jy=0,nyn(l)-1 64 73 phi = (ylat0n(l) + jy * dyn(l)) * pi / 180. … … 88 97 if (ix.eq.0.or.ix.eq.nxn(l)-1) jumpx=1 89 98 jux=jumpx 90 ! Precalculate pressure values for efficiency91 do kl=1,nuvz92 ppml(kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)93 end do94 99 ! 95 100 ! Loop over the vertical … … 97 102 98 103 do kl=1,nuvz 99 ppmk=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l) 100 theta=tthn(ix,jy,kl,n,l)*(100000./ppmk)**kappa 104 theta=tthn(ix,jy,kl,n,l)*ppmk(ix,jy,kl) 101 105 klvrp=kl+1 102 106 klvrm=kl-1 … … 107 111 if (klvrp.gt.nuvz) klvrp=nuvz 108 112 if (klvrm.lt.1) klvrm=1 109 ppmk=akz(klvrp)+bkz(klvrp)*psn(ix,jy,1,n,l) 110 thetap=tthn(ix,jy,klvrp,n,l)*(100000./ppmk)**kappa 111 ppmk=akz(klvrm)+bkz(klvrm)*psn(ix,jy,1,n,l) 112 thetam=tthn(ix,jy,klvrm,n,l)*(100000./ppmk)**kappa 113 dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm)) 113 thetap=tthn(ix,jy,klvrp,n,l)*ppmk(ix,jy,klvrp) 114 thetam=tthn(ix,jy,klvrm,n,l)*ppmk(ix,jy,klvrm) 115 dthetadp=(thetap-thetam)/(ppml(ix,jy,klvrp)-ppml(ix,jy,klvrm)) 114 116 115 117 ! Compute vertical position at pot. temperature surface on subgrid … … 134 136 kch=kch+1 135 137 k=kup 136 ppmk=akz(k)+bkz(k)*psn(ivr,jy,1,n,l) 137 thdn=tthn(ivr,jy,k,n,l)*(100000./ppmk)**kappa 138 ppmk=akz(k+1)+bkz(k+1)*psn(ivr,jy,1,n,l) 139 thup=tthn(ivr,jy,k+1,n,l)*(100000./ppmk)**kappa 138 thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k) 139 thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1) 140 140 141 141 if (((thdn.ge.theta).and.(thup.le.theta)).or. & … … 158 158 kch=kch+1 159 159 k=kdn 160 ppmk=akz(k)+bkz(k)*psn(ivr,jy,1,n,l) 161 thdn=tthn(ivr,jy,k,n,l)*(100000./ppmk)**kappa 162 ppmk=akz(k+1)+bkz(k+1)*psn(ivr,jy,1,n,l) 163 thup=tthn(ivr,jy,k+1,n,l)*(100000./ppmk)**kappa 160 thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k) 161 thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1) 164 162 if (((thdn.ge.theta).and.(thup.le.theta)).or. & 165 163 ((thdn.le.theta).and.(thup.ge.theta))) then … … 212 210 kch=kch+1 213 211 k=kup 214 ppmk=akz(k)+bkz(k)*psn(ix,j,1,n,l) 215 thdn=tthn(ix,j,k,n,l)*(100000./ppmk)**kappa 216 ppmk=akz(k+1)+bkz(k+1)*psn(ix,j,1,n,l) 217 thup=tthn(ix,j,k+1,n,l)*(100000./ppmk)**kappa 212 thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k) 213 thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1) 218 214 if (((thdn.ge.theta).and.(thup.le.theta)).or. & 219 215 ((thdn.le.theta).and.(thup.ge.theta))) then … … 235 231 kch=kch+1 236 232 k=kdn 237 ppmk=akz(k)+bkz(k)*psn(ix,j,1,n,l) 238 thdn=tthn(ix,j,k,n,l)*(100000./ppmk)**kappa 239 ppmk=akz(k+1)+bkz(k+1)*psn(ix,j,1,n,l) 240 thup=tthn(ix,j,k+1,n,l)*(100000./ppmk)**kappa 233 thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k) 234 thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1) 241 235 if (((thdn.ge.theta).and.(thup.le.theta)).or. & 242 236 ((thdn.le.theta).and.(thup.ge.theta))) then -
trunk/src/com_mod.f90
r20 r24 340 340 real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2) 341 341 real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2) 342 !scavenging NIK, PS343 342 integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2) 344 343 integer :: cloudsh(0:nxmax-1,0:nymax-1,2) 345 integer icloudbot(0:nxmax-1,0:nymax-1,2)346 integer icloudthck(0:nxmax-1,0:nymax-1,2)347 344 348 345 … … 359 356 ! rainout conv/lsp dominated 2/3 360 357 ! washout conv/lsp dominated 4/5 361 ! PS 2013362 !c icloudbot (m) cloud bottom height363 !c icloudthck (m) cloud thickness364 358 365 359 ! pplev for the GFS version -
trunk/src/gridcheck.f90
r20 r24 84 84 real(kind=4) :: xaux1,xaux2,yaux1,yaux2 85 85 real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in 86 integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl 86 integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parID 87 87 !HSO end 88 88 integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip 89 real :: sizesouth,sizenorth,xauxa,pint 89 real :: sizesouth,sizenorth,xauxa,pint,conversion_factor 90 90 91 91 ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING … … 98 98 99 99 integer :: isec1(56),isec2(22+nxmax+nymax) 100 !real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp) 101 real(kind=4) :: zsec2(184),zsec4(jpunp) 100 real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp) 102 101 character(len=1) :: opt 103 102 … … 172 171 call grib_check(iret,gribFunction,gribErrorMsg) 173 172 call grib_get_int(igrib,'level',valSurf,iret) 173 call grib_check(iret,gribFunction,gribErrorMsg) 174 call grib_get_int(igrib,'paramId',parId,iret) 174 175 call grib_check(iret,gribFunction,gribErrorMsg) 175 176 … … 193 194 elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot 194 195 isec1(6)=135 ! indicatorOfParameter 196 elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot 197 isec1(6)=135 ! indicatorOfParameter 195 198 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP 196 199 isec1(6)=151 ! indicatorOfParameter … … 205 208 elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD 206 209 isec1(6)=141 ! indicatorOfParameter 207 elseif ((parCat.eq.6).and.(parNum.eq.1) ) then ! CC210 elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC 208 211 isec1(6)=164 ! indicatorOfParameter 209 elseif ((parCat.eq.1).and.(parNum.eq.9) ) then ! LSP212 elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP 210 213 isec1(6)=142 ! indicatorOfParameter 211 214 elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP … … 215 218 elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR 216 219 isec1(6)=176 ! indicatorOfParameter 217 elseif ((parCat.eq.2).and.(parNum.eq.17) ) then ! EWSS220 elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS 218 221 isec1(6)=180 ! indicatorOfParameter 219 elseif ((parCat.eq.2).and.(parNum.eq.18) ) then ! NSSS222 elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS 220 223 isec1(6)=181 ! indicatorOfParameter 221 224 elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO 222 225 isec1(6)=129 ! indicatorOfParameter 223 elseif ((parCat.eq.3).and.(parNum.eq.7) ) then ! SDO226 elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO 224 227 isec1(6)=160 ! indicatorOfParameter 225 228 elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & … … 230 233 parCat,parNum,typSurf 231 234 endif 235 if(parId .ne. isec1(6) .and. parId .ne. 77) then 236 write(*,*) 'parId',parId, 'isec1(6)',isec1(6) 237 ! stop 238 endif 232 239 233 240 endif … … 265 272 266 273 !HSO get the second part of the grid dimensions only from GRiB1 messages 267 if ( (gribVer.eq.1).and.(gotGrid.eq.0)) then274 if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then 268 275 call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', & 269 276 xaux2in,iret) … … 291 298 dyconst=180./(dy*r_earth*pi) 292 299 gotGrid=1 293 294 300 ! Check whether fields are global 295 301 ! If they contain the poles, specify polar stereographic map … … 443 449 write(*,*) 444 450 write(*,'(a)') 'Mother domain:' 445 write(*,'(a,f10. 2,a1,f10.2,a,f10.2)') ' Longitude range: ', &451 write(*,'(a,f10.5,a,f10.5,a,f10.5)') ' Longitude range: ', & 446 452 xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx 447 write(*,'(a,f10. 2,a1,f10.2,a,f10.2)') ' Latitude range: ', &453 write(*,'(a,f10.5,a,f10.5,a,f10.5)') ' Latitude range: ', & 448 454 ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy 449 455 write(*,*) … … 467 473 k=nlev_ec+1+numskip+i 468 474 akm(nwz-i+1)=zsec2(j) 475 ! write (*,*) 'ifield:',ifield,k,j,zsec2(10+j) 469 476 bkm(nwz-i+1)=zsec2(k) 470 477 end do … … 553 560 554 561 end subroutine gridcheck 562 -
trunk/src/gridcheck_nests.f90
r4 r24 48 48 integer :: igrib 49 49 integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl 50 integer :: parID !added by mc for making it consistent with new gridcheck.f90 50 51 integer :: gotGrib 51 52 !HSO end … … 55 56 real(kind=4) :: xaux1,xaux2,yaux1,yaux2 56 57 real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in 58 real :: conversion_factor !added by mc to make it consistent with new gridchek.f90 57 59 58 60 ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING … … 150 152 call grib_get_int(igrib,'level',valSurf,iret) 151 153 call grib_check(iret,gribFunction,gribErrorMsg) 154 call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consisitent with new grid_check.f90 155 call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consisitent with new grid_check.f90 152 156 153 157 !print*,discipl,parCat,parNum,typSurf,valSurf … … 170 174 elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot 171 175 isec1(6)=135 ! indicatorOfParameter 176 elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added bymc to make it consistent with new gridchek.f90 177 isec1(6)=135 ! indicatorOfParameter ! ! " " " " " " " " " " " " " " " " " " " " " " " " " " " 172 178 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP 173 179 isec1(6)=151 ! indicatorOfParameter … … 182 188 elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD 183 189 isec1(6)=141 ! indicatorOfParameter 184 elseif ((parCat.eq.6).and.(parNum.eq.1) ) then ! CC190 elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consistent with new gridchek.f90 185 191 isec1(6)=164 ! indicatorOfParameter 186 elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP192 elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consistent with new gridchek.f90 187 193 isec1(6)=142 ! indicatorOfParameter 188 194 elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP … … 192 198 elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR 193 199 isec1(6)=176 ! indicatorOfParameter 194 elseif ((parCat.eq.2).and.(parNum.eq.17) ) then ! EWSS200 elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consistent with new gridchek.f90 195 201 isec1(6)=180 ! indicatorOfParameter 196 elseif ((parCat.eq.2).and.(parNum.eq.18) ) then ! NSSS202 elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consistent with new gridchek.f90 197 203 isec1(6)=181 ! indicatorOfParameter 198 204 elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO 199 205 isec1(6)=129 ! indicatorOfParameter 200 elseif ((parCat.eq.3).and.(parNum.eq.7) ) then ! SDO206 elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consistent with new gridchek.f90 201 207 isec1(6)=160 ! indicatorOfParameter 202 208 elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & … … 207 213 parCat,parNum,typSurf 208 214 endif 215 if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consistent with new gridchek.f90 216 write(*,*) 'parId',parId, 'isec1(6)',isec1(6) 217 ! stop 218 endif 209 219 210 220 endif … … 237 247 238 248 !HSO get the second part of the grid dimensions only from GRiB1 messages 239 if ( (gribVer.eq.1).and.(gotGrib.eq.0)) then240 call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', & 249 if (isec1(6) .eq. 167 .and. (gotGrib.eq.0)) then !added by mc to make it consistent with new gridchek.f90 note that gotGrid must be changed in gotGrib!! 250 call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', & !comment by mc: note that this was in the (if (ifield.eq.1) ..end above in gridchek.f90 see line 257 241 251 xaux1in,iret) 242 252 call grib_check(iret,gribFunction,gribErrorMsg) … … 254 264 yaux1=yaux1in 255 265 yaux2=yaux2in 256 if(xaux1.gt.180 ) xaux1=xaux1-360.0257 if(xaux2.gt.180 ) xaux2=xaux2-360.0258 if(xaux1.lt.-180 ) xaux1=xaux1+360.0259 if(xaux2.lt.-180 ) xaux2=xaux2+360.0260 if (xaux2.lt.xaux1) xaux2=xaux2+360. 266 if(xaux1.gt.180.) xaux1=xaux1-360.0 267 if(xaux2.gt.180.) xaux2=xaux2-360.0 268 if(xaux1.lt.-180.) xaux1=xaux1+360.0 269 if(xaux2.lt.-180.) xaux2=xaux2+360.0 270 if (xaux2.lt.xaux1) xaux2=xaux2+360.0 261 271 xlon0n(l)=xaux1 262 272 ylat0n(l)=yaux1 263 273 dxn(l)=(xaux2-xaux1)/real(nxn(l)-1) 264 274 dyn(l)=(yaux2-yaux1)/real(nyn(l)-1) 265 gotGrib=1 275 gotGrib=1 !commetn by mc note tahthere gotGRIB is used instead of gotGrid!!! 266 276 endif ! ifield.eq.1 267 277 … … 341 351 342 352 write(*,'(a,i2)') 'Nested domain #: ',l 343 write(*,'(a,f10. 2,a1,f10.2,a,f10.2)') ' Longitude range: ', &353 write(*,'(a,f10.5,a,f10.5,a,f10.5)') ' Longitude range: ', & 344 354 xlon0n(l),' to ',xlon0n(l)+(nxn(l)-1)*dxn(l), & 345 355 ' Grid distance: ',dxn(l) 346 write(*,'(a,f10. 2,a1,f10.2,a,f10.2)') ' Latitude range: ', &356 write(*,'(a,f10.5,a,f10.5,a,f10.5)') ' Latitude range: ', & 347 357 ylat0n(l),' to ',ylat0n(l)+(nyn(l)-1)*dyn(l), & 348 358 ' Grid distance: ',dyn(l) -
trunk/src/interpol_rain.f90
r20 r24 20 20 !********************************************************************** 21 21 22 !subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, & 23 ! ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3) 24 ! ! i i i i i i i 25 ! !i i i i i i i i o o o 26 27 subroutine interpol_rain(yy1,yy2,yy3,iy1,iy2,nxmax,nymax,nzmax,nx, & 28 ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, & 29 intiy1,intiy2,icmv) 30 ! i i i i i i i 31 ! i i i i i i i i o o o 32 22 subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, & 23 ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3) 24 ! i i i i i i i 25 !i i i i i i i i o o o 33 26 !**************************************************************************** 34 27 ! * … … 47 40 ! * 48 41 ! 30 August 1996 * 49 ! 50 !* Petra Seibert, 2011/2012: 51 !* Add interpolation of cloud bottom and cloud thickness 52 !* for fix to SE's new wet scavenging scheme * 42 ! * 53 43 !**************************************************************************** 54 44 ! * … … 80 70 81 71 integer :: nx,ny,nxmax,nymax,nzmax,memind(2),m,ix,jy,ixp,jyp 82 !integer :: itime,itime1,itime2,level,indexh 83 integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4 84 integer :: intiy1,intiy2,ipsum,icmv 72 integer :: itime,itime1,itime2,level,indexh 85 73 real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2) 86 74 real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2) 87 75 real :: yy3(0:nxmax-1,0:nymax-1,nzmax,2) 88 integer iy1(0:nxmax-1,0:nymax-1,2),iy2(0:nxmax-1,0:nymax-1,2) 89 real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),yi1(2),yi2(2) 90 !real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2) 91 real :: xt,yt,yint1,yint2,yint3,yint4,p1,p2,p3,p4 92 !real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4 76 real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2) 77 real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4 93 78 94 79 … … 142 127 + p3*yy3(ix ,jyp,level,indexh) & 143 128 + p4*yy3(ixp,jyp,level,indexh) 144 145 !CPS clouds:146 ip1=1147 ip2=1148 ip3=1149 ip4=1150 if (iy1(ix ,jy ,indexh) .eq. icmv) ip1=0151 if (iy1(ixp,jy ,indexh) .eq. icmv) ip2=0152 if (iy1(ix ,jyp,indexh) .eq. icmv) ip3=0153 if (iy1(ixp,jyp,indexh) .eq. icmv) ip4=0154 ipsum= ip1+ip2+ip3+ip4155 if (ipsum .eq. 0) then156 yi1(m)=icmv157 else158 yi1(m)=(ip1*p1*iy1(ix ,jy ,indexh) &159 + ip2*p2*iy1(ixp,jy ,indexh) &160 + ip3*p3*iy1(ix ,jyp,indexh) &161 + ip4*p4*iy1(ixp,jyp,indexh))/ipsum162 endif163 164 ip1=1165 ip2=1166 ip3=1167 ip4=1168 if (iy2(ix ,jy ,indexh) .eq. icmv) ip1=0169 if (iy2(ixp,jy ,indexh) .eq. icmv) ip2=0170 if (iy2(ix ,jyp,indexh) .eq. icmv) ip3=0171 if (iy2(ixp,jyp,indexh) .eq. icmv) ip4=0172 ipsum= ip1+ip2+ip3+ip4173 if (ipsum .eq. 0) then174 yi2(m)=icmv175 else176 yi2(m)=(ip1*p1*iy2(ix ,jy ,indexh) &177 + ip2*p2*iy2(ixp,jy ,indexh) &178 + ip3*p3*iy2(ix ,jyp,indexh) &179 + ip4*p4*iy2(ixp,jyp,indexh))/ipsum180 endif181 !CPS end clouds182 183 129 end do 184 130 … … 187 133 ! 2.) Temporal interpolation (linear) 188 134 !************************************ 189 190 if (abs(itime) .lt. abs(itime1)) then191 print*,'interpol_rain.f90'192 print*,itime,itime1,itime2193 stop 'ITIME PROBLEM'194 endif195 196 135 197 136 dt1=real(itime-itime1) … … 204 143 205 144 206 !PS clouds:207 intiy1=(yi1(1)*dt2 + yi1(2)*dt1)/dt208 if (yi1(1) .eq. float(icmv)) intiy1=yi1(2)209 if (yi1(2) .eq. float(icmv)) intiy1=yi1(1)210 211 intiy2=(yi2(1)*dt2 + yi2(2)*dt1)/dt212 if (yi2(1) .eq. float(icmv)) intiy2=yi2(2)213 if (yi2(2) .eq. float(icmv)) intiy2=yi2(1)214 215 if (intiy1 .ne. icmv .and. intiy2 .ne. icmv) then216 intiy2 = intiy1 + intiy2 ! convert cloud thickness to cloud top217 else218 intiy1=icmv219 intiy2=icmv220 endif221 !PS end clouds222 223 145 end subroutine interpol_rain -
trunk/src/makefile
r20 r24 1 1 SHELL = /bin/bash 2 TARGET = local 3 WINDS=ecmwf 4 #WINDS=gfs 5 #WINDS=fnl 2 MAIN = FP_ecmwf_gfortran 6 3 7 4 FC = gfortran 5 INCPATH = /xnilu_wrk/flex_wrk/bin64/grib_api/include 6 LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib 7 LIBPATH2 = /usr/lib/x86_64-linux-gnu/ 8 8 FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH) 9 #FFLAGS = -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH) 10 LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper 9 11 10 ifeq ($(TARGET),dmz)11 # options for ganglia12 INCPATH = /xnilu_wrk/flex_wrk/bin64/grib_api/include13 LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib14 LIBPATH2 = /usr/lib/x86_64-linux-gnu/15 MAIN = FLEXPART_dmz_16 endif17 ifeq ($(TARGET),local)18 # local options19 #libs_dir=/.../flexpart/libs/20 libs_dir=/Users/ignacio/flexpart/libs/21 INCPATH = $(libs_dir)/grib_api-1.9.9_dir/include22 LIBPATH1 = $(libs_dir)/grib_api-1.9.9_dir/lib23 LIBPATH2 = $(libs_dir)/jasper_dir/lib24 MAIN = FLEXPART_local_25 endif26 27 LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper28 12 29 13 MODOBJS = \ … … 37 21 OBJECTS = \ 38 22 writeheader.o writeheader_txt.o writeheader_surf.o assignland.o\ 39 part0.o \23 calcpar.o part0.o \ 40 24 caldate.o partdep.o \ 41 25 coordtrafo.o psih.o \ … … 53 37 getrc.o readreleases.o \ 54 38 getvdep.o readspecies.o \ 55 interpol_misslev.o \56 conccalc.o \39 interpol_misslev.o readwind.o \ 40 conccalc.o richardson.o \ 57 41 concoutput.o concoutput_surf.o scalev.o \ 58 42 pbl_profile.o readOHfield.o\ 59 43 juldate.o timemanager.o \ 60 44 interpol_vdep.o interpol_rain.o \ 61 partoutput.o \45 verttransform.o partoutput.o \ 62 46 hanna.o wetdepokernel.o \ 63 47 mean.o wetdepo.o \ 64 48 hanna_short.o windalign.o \ 49 obukhov.o gridcheck.o \ 65 50 hanna1.o initialize.o \ 66 calcpar_nests.o \ 51 gridcheck_nests.o \ 52 readwind_nests.o calcpar_nests.o \ 67 53 verttransform_nests.o interpol_all_nests.o \ 68 54 interpol_wind_nests.o interpol_misslev_nests.o \ 69 55 interpol_vdep_nests.o interpol_rain_nests.o \ 70 getvdep_nests.o gridcheck_nests.o \ 71 readwind_nests.o \ 56 getvdep_nests.o \ 72 57 readageclasses.o readpartpositions.o \ 73 58 calcfluxes.o fluxoutput.o \ 74 59 qvsat.o skplin.o \ 60 convmix.o calcmatrix.o \ 75 61 convect43c.o redist.o \ 76 62 sort2.o distance.o \ … … 91 77 92 78 93 ifeq ($(WINDS),ecmwf) 94 OBJECTS_WINDS = \ 95 calcpar.o readwind.o \ 96 richardson.o verttransform.o \ 97 obukhov.o gridcheck.o \ 98 convmix.o calcmatrix.o 99 endif 100 101 ifeq ($(WINDS),gfs) 102 OBJECTS_WINDS = \ 103 calcpar_gfs.o readwind_gfs.o \ 104 richardson_gfs.o verttransform_gfs.o \ 105 obukhov_gfs.o gridcheck_gfs.o \ 106 convmix_gfs.o calcmatrix_gfs.o 107 endif 108 109 $(MAIN): $(MODOBJS) $(OBJECTS) $(OBJECTS_WINDS) 110 $(FC) *.o -o $(MAIN)_$(WINDS) $(LDFLAGS) 79 $(MAIN): $(MODOBJS) $(OBJECTS) 80 $(FC) *.o -o $(MAIN) $(LDFLAGS) 111 81 112 82 $(OBJECTS): $(MODOBJS) … … 118 88 rm *.o *.mod 119 89 90 cleanall: 91 rm *.o *.mod $(MAIN) -
trunk/src/par_mod.f90
r20 r24 121 121 !********************************************* 122 122 123 integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !FNL 123 !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !FNL XF 124 integer,parameter :: nxmax=361,nymax=181,nuvzmax=152,nwzmax=152,nzmax=152 !ECMWF new 124 125 !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF 125 126 !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26 … … 154 155 155 156 !integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0 156 !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151 !ECMWF157 integer,parameter :: maxnests=1, nxmaxn=201, nymaxn=161 ! FNL XF157 integer,parameter :: maxnests=1,nxmaxn=351,nymaxn=351 !ECMWF 158 !integer,parameter :: maxnests=1, nxmaxn=201, nymaxn=161 ! FNL XF 158 159 ! maxnests maximum number of nested grids 159 160 ! nxmaxn,nymaxn maximum dimension of nested wind fields in … … 197 198 !************************************************** 198 199 199 !integer,parameter :: maxpart=4000000 200 integer,parameter :: maxpart=2000 201 integer,parameter :: maxspec=6 200 integer,parameter :: maxpart=15000000 201 integer,parameter :: maxspec=4 202 202 203 203 -
trunk/src/readcommand.f90
r20 r24 80 80 character(len=50) :: line 81 81 logical :: old 82 logical :: nml out=.true.!.false.82 logical :: nml_COMMAND=.true. , nmlout=.true. !.false. 83 83 integer :: readerror 84 84 85 85 namelist /command/ & 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 86 ldirect, & 87 ibdate,ibtime, & 88 iedate,ietime, & 89 loutstep, & 90 loutaver, & 91 loutsample, & 92 itsplit, & 93 lsynctime, & 94 ctl, & 95 ifine, & 96 iout, & 97 ipout, & 98 lsubgrid, & 99 lconvection, & 100 lagespectra, & 101 ipin, & 102 ioutputforeachrelease, & 103 iflux, & 104 mdomainfill, & 105 ind_source, & 106 ind_receptor, & 107 mquasilag, & 108 nested_output, & 109 linit_cond, & 110 surf_only 111 111 112 112 ! Presetting namelist command … … 144 144 open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', & 145 145 form='formatted',iostat=readerror) 146 ! If fail, check if file does not exist 147 if (readerror.ne.0) then 148 149 print*,'***ERROR: file COMMAND not found in ' 150 print*, path(1)(1:length(1))//'COMMAND' 151 print*, 'Check your pathnames file.' 152 stop 153 154 endif 155 156 read(unitcommand,command,iostat=readerror) 146 ! If fail, check if file does not exist 147 if (readerror.ne.0) then 148 print*,'***ERROR: file COMMAND not found in ' 149 print*, path(1)(1:length(1))//'COMMAND' 150 print*, 'Check your pathnames file.' 151 stop 152 endif 153 154 ! print error code 155 !write(*,*) 'readcommand > readerror open=' , readerror 156 !probe first line 157 read (unitcommand,901) line 158 !write(*,*) 'index(line,COMMAND) =', index(line,'COMMAND') 159 160 161 !default is namelist input 162 ! distinguish namelist from fixed text input 163 if (index(line,'COMMAND') .eq. 0) then 164 nml_COMMAND = .false. 165 !write(*,*) 'COMMAND file does not contain the string COMMAND in the first line' 166 endif 167 !write(*,*) 'readcommand > read as namelist? ' , nml_COMMAND 168 rewind(unitcommand) 169 read(unitcommand,command,iostat=readerror) 170 157 171 close(unitcommand) 158 172 173 !write(*,*) 'readcommand > readerror read=' , readerror 159 174 ! If error in namelist format, try to open with old input code 160 if (readerror.ne.0) then 175 ! if (readerror.ne.0) then 176 ! IP 21/5/2014 the previous line cause the old long format 177 ! to be confused with namelist input 178 179 ! use text input 180 if (nml_COMMAND .eqv. .false.) then 161 181 162 182 open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', & … … 173 193 if (index(line,'LDIRECT') .eq. 0) then 174 194 old = .false. 195 !write(*,*) 'readcommand old short' 175 196 else 176 197 old = .true. 198 !write(*,*) 'readcommand old long' 177 199 endif 178 200 rewind(unitcommand) 179 201 202 180 203 ! Read parameters 181 204 !**************** … … 231 254 if (old) call skplin(3,unitcommand) 232 255 read(unitcommand,*) linit_cond 256 if (old) call skplin(3,unitcommand) 257 read(unitcommand,*) surf_only 233 258 close(unitcommand) 234 259 … … 237 262 ! write command file in namelist format to output directory if requested 238 263 if (nmlout.eqv..true.) then 239 !open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist.out',status='new',err=1000)240 264 open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist',err=1000) 241 265 write(unitcommand,nml=command) -
trunk/src/readreleases.f90
r20 r24 55 55 ! xpoint2,ypoint2 geograf. coordinates of upper right corner of release * 56 56 ! area * 57 ! weta, wetb parameters to determine the wet scavenging coefficient * 57 ! weta, wetb parameters to determine the below-cloud scavenging * 58 ! weta_in, wetb_in parameters to determine the in-cloud scavenging * 59 ! wetc_in, wetd_in parameters to determine the in-cloud scavenging * 58 60 ! zpoint1,zpoint2 height range, over which release takes place * 59 61 ! num_min_discrete if less, release cannot be randomized and happens at * … … 272 274 vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j) 273 275 end do 274 write(*,*) 'Average sett ing velocity: ',i,vsetaver(i)276 write(*,*) 'Average settling velocity: ',i,vsetaver(i) 275 277 endif 276 278 … … 284 286 if (weta(i).gt.0.) then 285 287 WETDEP=.true. 286 write (*,*) 'Wetdeposition switched on: ',weta(i),i 287 endif 288 write (*,*) 'Below-cloud scavenging is switched on' 289 write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i 290 else 291 write (*,*) 'Below-cloud scavenging is switched OFF' 292 endif 293 294 ! NIK 31.01.2013 + 10.12.2013 295 if (weta_in(i).gt.0.) then 296 WETDEP=.true. 297 write (*,*) 'In-cloud scavenging is switched on' 298 write (*,*) 'In-cloud scavenging coefficients: ',weta_in(i),wetb_in(i), wetc_in(i), wetd_in(i),i 299 else 300 write (*,*) 'In-cloud scavenging is switched OFF' 301 endif 302 288 303 if (ohreact(i).gt.0) then 289 304 OHREA=.true. -
trunk/src/readspecies.f90
r4 r24 31 31 ! 11 July 1996 * 32 32 ! * 33 ! Changes: * 34 ! N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging * 35 ! * 33 36 !***************************************************************************** 34 37 ! * … … 36 39 ! decaytime(maxtable) half time for radiological decay * 37 40 ! specname(maxtable) names of chemical species, radionuclides * 38 ! wetscava, wetscavb Parameters for determining scavenging coefficient * 41 ! weta, wetb Parameters for determining below-cloud scavenging * 42 ! weta_in Parameter for determining in-cloud scavenging * 43 ! wetb_in Parameter for determining in-cloud scavenging * 44 ! wetc_in Parameter for determining in-cloud scavenging * 45 ! wetd_in Parameter for determining in-cloud scavenging * 39 46 ! ohreact OH reaction rate * 40 47 ! id_spec SPECIES number as referenced in RELEASE file * … … 78 85 read(unitspecies,'(f18.2)',end=22) wetb(pos_spec) 79 86 ! write(*,*) wetb(pos_spec) 87 88 !*** NIK 31.01.2013: including in-cloud scavening parameters 89 read(unitspecies,'(e18.1)',end=22) weta_in(pos_spec) 90 ! write(*,*) weta_in(pos_spec) 91 read(unitspecies,'(f18.2)',end=22) wetb_in(pos_spec) 92 ! write(*,*) wetb_in(pos_spec) 93 read(unitspecies,'(f18.2)',end=22) wetc_in(pos_spec) 94 ! write(*,*) wetc_in(pos_spec) 95 read(unitspecies,'(f18.2)',end=22) wetd_in(pos_spec) 96 ! write(*,*) wetd_in(pos_spec) 97 80 98 read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec) 81 99 ! write(*,*) reldiff(pos_spec) -
trunk/src/readwind.f90
r20 r24 71 71 integer :: iret 72 72 integer :: igrib 73 integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl 73 integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId 74 74 integer :: gotGrid 75 75 !HSO end … … 90 90 integer :: isec1(56),isec2(22+nxmax+nymax) 91 91 real(kind=4) :: zsec4(jpunp) 92 real(kind=4) :: xaux,yaux ,xaux0,yaux092 real(kind=4) :: xaux,yaux 93 93 real(kind=8) :: xauxin,yauxin 94 94 real,parameter :: eps=1.e-4 95 95 real(kind=4) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1) 96 real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1 96 real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor 97 97 98 98 logical :: hflswitch,strswitch … … 101 101 character(len=24) :: gribErrorMsg = 'Error reading grib file' 102 102 character(len=20) :: gribFunction = 'readwind' 103 104 !HSO conversion of ECMWF etadot to etadot*dp/deta105 logical :: etacon=.false.106 real,parameter :: p00=101325.107 real :: dak,dbk108 103 109 104 hflswitch=.false. … … 155 150 endif 156 151 152 conversion_factor=1. 153 157 154 else 158 155 … … 168 165 call grib_check(iret,gribFunction,gribErrorMsg) 169 166 call grib_get_int(igrib,'level',valSurf,iret) 167 call grib_check(iret,gribFunction,gribErrorMsg) 168 call grib_get_int(igrib,'paramId',parId,iret) 170 169 call grib_check(iret,gribFunction,gribErrorMsg) 171 170 … … 177 176 isec1(8)=-1 178 177 isec1(8)=valSurf ! level 178 conversion_factor=1. 179 179 if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T 180 180 isec1(6)=130 ! indicatorOfParameter … … 188 188 isec1(6)=134 ! indicatorOfParameter 189 189 elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot 190 isec1(6)=135 ! indicatorOfParameter 191 elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot 190 192 isec1(6)=135 ! indicatorOfParameter 191 193 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP … … 201 203 elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD 202 204 isec1(6)=141 ! indicatorOfParameter 203 elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC 205 conversion_factor=1000. 206 elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC 204 207 isec1(6)=164 ! indicatorOfParameter 205 elseif ((parCat.eq.1).and.(parNum.eq.9) ) then ! LSP208 elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP 206 209 isec1(6)=142 ! indicatorOfParameter 207 210 elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP 208 211 isec1(6)=143 ! indicatorOfParameter 212 conversion_factor=1000. 209 213 elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF 210 214 isec1(6)=146 ! indicatorOfParameter 211 215 elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR 212 216 isec1(6)=176 ! indicatorOfParameter 213 elseif ((parCat.eq.2).and.(parNum.eq.17) ) then ! EWSS217 elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS 214 218 isec1(6)=180 ! indicatorOfParameter 215 elseif ((parCat.eq.2).and.(parNum.eq.18) ) then ! NSSS219 elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS 216 220 isec1(6)=181 ! indicatorOfParameter 217 221 elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO 218 222 isec1(6)=129 ! indicatorOfParameter 219 elseif ((parCat.eq.3).and.(parNum.eq.7) ) then ! SDO223 elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO 220 224 isec1(6)=160 ! indicatorOfParameter 221 225 elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & … … 223 227 isec1(6)=172 ! indicatorOfParameter 224 228 else 225 print*,'*** ERROR: undefined GRiB2 message found!',discipl, &229 print*,'***WARNING: undefined GRiB2 message found!',discipl, & 226 230 parCat,parNum,typSurf 231 endif 232 if(parId .ne. isec1(6) .and. parId .ne. 77) then 233 write(*,*) 'parId',parId, 'isec1(6)',isec1(6) 234 ! stop 227 235 endif 228 236 … … 237 245 !HSO get the required fields from section 2 in a gribex compatible manner 238 246 if (ifield.eq.1) then 239 call grib_get_int(igrib,'numberOfPointsAlongAParallel', & 240 isec2(2),iret) 241 call grib_check(iret,gribFunction,gribErrorMsg) 242 call grib_get_int(igrib,'numberOfPointsAlongAMeridian', & 243 isec2(3),iret) 244 call grib_check(iret,gribFunction,gribErrorMsg) 245 call grib_get_int(igrib,'numberOfVerticalCoordinateValues', & 246 isec2(12)) 247 call grib_get_int(igrib,'numberOfPointsAlongAParallel',isec2(2),iret) 248 call grib_check(iret,gribFunction,gribErrorMsg) 249 call grib_get_int(igrib,'numberOfPointsAlongAMeridian',isec2(3),iret) 250 call grib_check(iret,gribFunction,gribErrorMsg) 251 call grib_get_int(igrib,'numberOfVerticalCoordinateValues',isec2(12)) 247 252 call grib_check(iret,gribFunction,gribErrorMsg) 248 253 ! CHECK GRID SPECIFICATIONS … … 250 255 if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT' 251 256 if(isec2(12)/2-1.ne.nlev_ec) & 252 257 stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT' 253 258 endif ! ifield 254 259 255 260 !HSO get the second part of the grid dimensions only from GRiB1 messages 256 if ( (gribVer.eq.1).and.(gotGrid.eq.0)) then261 if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then 257 262 call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', & 258 263 xauxin,iret) … … 261 266 yauxin,iret) 262 267 call grib_check(iret,gribFunction,gribErrorMsg) 268 if (xauxin.gt.180.) xauxin=xauxin-360.0 269 if (xauxin.lt.-180.) xauxin=xauxin+360.0 270 263 271 xaux=xauxin+real(nxshift)*dx 264 272 yaux=yauxin 265 xaux0=xlon0 266 yaux0=ylat0 267 if(xaux.lt.0.) xaux=xaux+360. 268 if(yaux.lt.0.) yaux=yaux+360. 269 if(xaux0.lt.0.) xaux0=xaux0+360. 270 if(yaux0.lt.0.) yaux0=yaux0+360. 271 if(abs(xaux-xaux0).gt.eps) & 272 stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT' 273 if(abs(yaux-yaux0).gt.eps) & 274 stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT' 273 if (xaux.gt.180.) xaux=xaux-360.0 274 if(abs(xaux-xlon0).gt.eps) & 275 stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT' 276 if(abs(yaux-ylat0).gt.eps) & 277 stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT' 275 278 gotGrid=1 276 279 endif ! gotGrid … … 298 301 zsec4(nxfield*(ny-j-1)+i+1) 299 302 if(isec1(6).eq.141) sd(i,j,1,n)= &!! SNOW DEPTH 300 zsec4(nxfield*(ny-j-1)+i+1) 303 zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor 301 304 if(isec1(6).eq.151) msl(i,j,1,n)= &!! SEA LEVEL PRESS. 302 305 zsec4(nxfield*(ny-j-1)+i+1) … … 316 319 endif 317 320 if(isec1(6).eq.143) then !! CONVECTIVE PREC. 318 convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1) 321 convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor 319 322 if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0. 320 323 endif … … 366 369 do j=0,nymin1 367 370 wwh(i,j,nlev_ec+1)=0. 368 end do369 end do370 endif371 372 ! convert from ECMWF etadot to etadot*dp/deta as needed by FLEXPART373 if(etacon.eqv..true.) then374 do k=1,nwzmax375 dak=akm(k+1)-akm(k)376 dbk=bkm(k+1)-bkm(k)377 do i=0,nxmin1378 do j=0,nymin1379 wwh(i,j,k)=2*wwh(i,j,k)*ps(i,j,1,n)*(dak/ps(i,j,1,n)+dbk)/(dak/p00+dbk)380 if (k.gt.1) then381 wwh(i,j,k)=wwh(i,j,k)-wwh(i,j,k-1)382 endif383 end do384 371 end do 385 372 end do … … 479 466 480 467 end subroutine readwind 468 -
trunk/src/readwind_nests.f90
r4 r24 51 51 integer :: igrib 52 52 integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl 53 integer :: parId !!added by mc for making it consistent with new readwind.f90 53 54 integer :: gotGrid 54 55 !HSO end … … 69 70 integer :: isec1(56),isec2(22+nxmaxn+nymaxn) 70 71 real(kind=4) :: zsec4(jpunp) 72 real(kind=4) :: xaux,yaux 71 73 real(kind=8) :: xauxin,yauxin 72 real (kind=4) :: xaux,yaux,xaux0,yaux074 real,parameter :: eps=1.e-4 73 75 real :: ewss(0:nxmaxn-1,0:nymaxn-1),nsss(0:nxmaxn-1,0:nymaxn-1) 74 76 real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1 77 real :: conversion_factor !added by mc to make it consistent with new gridchek.f90 75 78 76 79 logical :: hflswitch,strswitch … … 134 137 endif 135 138 139 conversion_factor=1. 140 141 136 142 else 137 143 … … 148 154 call grib_get_int(igrib,'level',valSurf,iret) 149 155 call grib_check(iret,gribFunction,gribErrorMsg) 156 call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consisitent with new readwind.f90 157 call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consisitent with new readwind.f90 150 158 151 159 !print*,discipl,parCat,parNum,typSurf,valSurf … … 156 164 isec1(8)=-1 157 165 isec1(8)=valSurf ! level 166 conversion_factor=1. 158 167 if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T 159 168 isec1(6)=130 ! indicatorOfParameter … … 166 175 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP 167 176 isec1(6)=134 ! indicatorOfParameter 168 elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot 177 elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot ! 169 178 isec1(6)=135 ! indicatorOfParameter 179 elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added by mc to make it consisitent with new readwind.f90 180 isec1(6)=135 ! indicatorOfParameter !added by mc to make it consisitent with new readwind.f90 170 181 elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP 171 182 isec1(6)=151 ! indicatorOfParameter … … 180 191 elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD 181 192 isec1(6)=141 ! indicatorOfParameter 182 elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC 193 conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90 194 elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consisitent with new readwind.f90 183 195 isec1(6)=164 ! indicatorOfParameter 184 elseif ((parCat.eq.1).and.(parNum.eq.9) ) then ! LSP196 elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consisitent with new readwind.f90 185 197 isec1(6)=142 ! indicatorOfParameter 186 198 elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP 187 199 isec1(6)=143 ! indicatorOfParameter 200 conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90 188 201 elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF 189 202 isec1(6)=146 ! indicatorOfParameter 190 203 elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR 191 204 isec1(6)=176 ! indicatorOfParameter 192 elseif ((parCat.eq.2).and.(parNum.eq.17) ) then ! EWSS205 elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consisitent with new readwind.f90 193 206 isec1(6)=180 ! indicatorOfParameter 194 elseif ((parCat.eq.2).and.(parNum.eq.18) ) then ! NSSS207 elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consisitent with new readwind.f90 195 208 isec1(6)=181 ! indicatorOfParameter 196 209 elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO 197 210 isec1(6)=129 ! indicatorOfParameter 198 elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO211 elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consisitent with new readwind.f90 199 212 isec1(6)=160 ! indicatorOfParameter 200 213 elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. & … … 202 215 isec1(6)=172 ! indicatorOfParameter 203 216 else 204 print*,'*** ERROR: undefined GRiB2 message found!',discipl, &217 print*,'***WARNING: undefined GRiB2 message found!',discipl, & 205 218 parCat,parNum,typSurf 219 endif 220 if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consisitent with new readwind.f90 221 write(*,*) 'parId',parId, 'isec1(6)',isec1(6) ! 222 ! stop 206 223 endif 207 224 … … 227 244 ! CHECK GRID SPECIFICATIONS 228 245 if(isec2(2).ne.nxn(l)) stop & 229 246 'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL' 230 247 if(isec2(3).ne.nyn(l)) stop & 231 248 'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL' 232 249 if(isec2(12)/2-1.ne.nlev_ec) stop 'READWIND: VERTICAL DISCRET& 233 &IZATION NOT CONSISTENT FOR A NESTING LEVEL'250 IZATION NOT CONSISTENT FOR A NESTING LEVEL' 234 251 endif ! ifield 235 252 236 253 !HSO get the second part of the grid dimensions only from GRiB1 messages 237 if ((gribVer.eq.1).and.(gotGrid.eq.0)) then254 if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consisitent with new readwind.f90 238 255 call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', & 239 256 xauxin,iret) … … 242 259 yauxin,iret) 243 260 call grib_check(iret,gribFunction,gribErrorMsg) 261 if (xauxin.gt.180.) xauxin=xauxin-360.0 262 if (xauxin.lt.-180.) xauxin=xauxin+360.0 263 244 264 xaux=xauxin 245 265 yaux=yauxin 246 xaux0=xlon0n(l) 247 yaux0=ylat0n(l) 248 if(xaux.lt.0.) xaux=xaux+360. 249 if(yaux.lt.0.) yaux=yaux+360. 250 if(xaux0.lt.0.) xaux0=xaux0+360. 251 if(yaux0.lt.0.) yaux0=yaux0+360. 252 if(xaux.ne.xaux0) & 253 stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT FOR A NES& 254 &TING LEVEL' 255 if(yaux.ne.yaux0) & 256 stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT FOR A NEST& 257 &ING LEVEL' 266 if (abs(xaux-xlon0n(l)).gt.eps) & 267 stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT FOR A NESTING LEVEL' 268 if (abs(yaux-ylat0n(l)).gt.eps) & 269 stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT FOR A NESTING LEVEL' 258 270 gotGrid=1 259 271 endif … … 280 292 zsec4(nxn(l)*(nyn(l)-j-1)+i+1) 281 293 if(isec1(6).eq.141) sdn(i,j,1,n,l)= &!! SNOW DEPTH 282 zsec4(nxn(l)*(nyn(l)-j-1)+i+1) 294 zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90! 283 295 if(isec1(6).eq.151) msln(i,j,1,n,l)= &!! SEA LEVEL PRESS. 284 296 zsec4(nxn(l)*(nyn(l)-j-1)+i+1) … … 298 310 endif 299 311 if(isec1(6).eq.143) then !! CONVECTIVE PREC. 300 convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1) 312 convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90 301 313 if (convprecn(i,j,1,n,l).lt.0.) convprecn(i,j,1,n,l)=0. 302 314 endif -
trunk/src/verttransform.f90
r20 r24 21 21 22 22 subroutine verttransform(n,uuh,vvh,wwh,pvh) 23 ! 23 ! i i i i i 24 24 !***************************************************************************** 25 25 ! * … … 49 49 ! Sabine Eckhardt, March 2007 50 50 ! added the variable cloud for use with scavenging - descr. in com_mod 51 ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification52 ! note that also other subroutines are affected by the fix53 51 !***************************************************************************** 54 52 ! * … … 72 70 73 71 integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym 74 integer :: rain_cloud_above,kz_inv !SE 75 integer icloudtop !PS 72 integer :: rain_cloud_above,kz_inv 76 73 real :: f_qvsat,pressure 77 !real :: rh,lsp,convp 78 real :: rh,lsp,convp,prec,rhmin 79 real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax) 74 real :: rh,lsp,convp 75 real :: rhoh(nuvzmax),pinmconv(nzmax) 80 76 real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi 81 77 real :: xlon,ylat,xlonr,dzdx,dzdy 82 real :: dzdx1,dzdx2,dzdy1,dzdy2, precmin78 real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf 83 79 real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy 84 80 real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax) … … 87 83 real :: wwh(0:nxmax-1,0:nymax-1,nwzmax) 88 84 real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax) 89 logical lconvectprec90 85 real,parameter :: const=r_air/ga 91 parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics92 86 93 87 logical :: init = .true. … … 97 91 ! If verttransform is called the first time, initialize heights of the * 98 92 ! z levels in meter. The heights are the heights of model levels, where * 99 ! u,v,T and qv are given, and of the interfaces, where w is given. So, * 100 ! the vertical resolution in the z system is doubled. As reference point,* 101 ! the lower left corner of the grid is used. * 102 ! Unlike in the eta system, no difference between heights for u,v and * 103 ! heights for w exists. * 93 ! u,v,T and qv are given. * 104 94 !************************************************************************* 105 106 107 ! do 897 kz=1,nuvz108 ! write (*,*) 'akz: ',akz(kz),'bkz',bkz(kz)109 !897 continue110 95 111 96 if (init) then … … 113 98 ! Search for a point with high surface pressure (i.e. not above significant topography) 114 99 ! Then, use this point to construct a reference z profile, to be used at all times 115 !***************************************************************************** 100 !************************************************************************************** 116 101 117 102 do jy=0,nymin1 … … 128 113 129 114 tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ & 130 115 ps(ixm,jym,1,n)) 131 116 pold=ps(ixm,jym,1,n) 132 117 height(1)=0. … … 136 121 tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n)) 137 122 138 139 ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled140 ! upon the transformation to z levels. In order to save computer memory, this is141 ! not done anymore in the standard version. However, this option can still be142 ! switched on by replacing the following lines with those below, that are143 ! currently commented out.144 ! Note that two more changes are necessary in this subroutine below.145 ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.146 !*****************************************************************************147 148 123 if (abs(tv-tvold).gt.0.2) then 149 height(kz)= & 150 height(kz-1)+const*log(pold/pint)* & 151 (tv-tvold)/log(tv/tvold) 124 height(kz)=height(kz-1)+const*log(pold/pint)* & 125 (tv-tvold)/log(tv/tvold) 152 126 else 153 height(kz)=height(kz-1)+ & 154 const*log(pold/pint)*tv 127 height(kz)=height(kz-1)+const*log(pold/pint)*tv 155 128 endif 156 157 ! Switch on following lines to use doubled vertical resolution158 !*************************************************************159 ! if (abs(tv-tvold).gt.0.2) then160 ! height((kz-1)*2)=161 ! + height(max((kz-2)*2,1))+const*log(pold/pint)*162 ! + (tv-tvold)/log(tv/tvold)163 ! else164 ! height((kz-1)*2)=height(max((kz-2)*2,1))+165 ! + const*log(pold/pint)*tv166 ! endif167 ! End doubled vertical resolution168 129 169 130 tvold=tv 170 131 pold=pint 171 132 end do 172 173 174 ! Switch on following lines to use doubled vertical resolution175 !*************************************************************176 ! do 7 kz=3,nz-1,2177 ! height(kz)=0.5*(height(kz-1)+height(kz+1))178 ! height(nz)=height(nz-1)+height(nz-1)-height(nz-2)179 ! End doubled vertical resolution180 133 181 134 … … 204 157 do jy=0,nymin1 205 158 do ix=0,nxmin1 206 tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ & 207 ps(ix,jy,1,n)) 159 tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ps(ix,jy,1,n)) 208 160 pold=ps(ix,jy,1,n) 209 uv zlev(1)=0.161 uvwzlev(ix,jy,1)=0. 210 162 wzlev(1)=0. 211 163 rhoh(1)=pold/(r_air*tvold) … … 221 173 222 174 if (abs(tv-tvold).gt.0.2) then 223 uv zlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &224 175 uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* & 176 (tv-tvold)/log(tv/tvold) 225 177 else 226 uv zlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv178 uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv 227 179 endif 228 180 … … 233 185 234 186 do kz=2,nwz-1 235 wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2. 236 end do 237 wzlev(nwz)=wzlev(nwz-1)+ & 238 uvzlev(nuvz)-uvzlev(nuvz-1) 239 240 uvwzlev(ix,jy,1)=0.0 241 do kz=2,nuvz 242 uvwzlev(ix,jy,kz)=uvzlev(kz) 243 end do 244 245 ! Switch on following lines to use doubled vertical resolution 246 ! Switch off the three lines above. 247 !************************************************************* 248 !22 uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz) 249 ! do 23 kz=2,nwz 250 !23 uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz) 251 ! End doubled vertical resolution 187 wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2. 188 end do 189 wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1) 252 190 253 191 ! pinmconv=(h2-h1)/(p2-p1) 254 192 255 193 pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ & 256 257 194 ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- & 195 (aknew(1)+bknew(1)*ps(ix,jy,1,n))) 258 196 do kz=2,nz-1 259 197 pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ & 260 261 198 ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- & 199 (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n))) 262 200 end do 263 201 pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ & 264 265 202 ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- & 203 (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n))) 266 204 267 205 ! Levels, where u,v,t and q are given … … 283 221 do iz=2,nz-1 284 222 do kz=kmin,nuvz 285 if(height(iz).gt.uv zlev(nuvz)) then223 if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then 286 224 uu(ix,jy,iz,n)=uu(ix,jy,nz,n) 287 225 vv(ix,jy,iz,n)=vv(ix,jy,nz,n) … … 292 230 goto 30 293 231 endif 294 if ((height(iz).gt.uv zlev(kz-1)).and. &295 (height(iz).le.uvzlev(kz))) then296 dz1=height(iz)-uv zlev(kz-1)297 dz2=uv zlev(kz)-height(iz)232 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. & 233 (height(iz).le.uvwzlev(ix,jy,kz))) then 234 dz1=height(iz)-uvwzlev(ix,jy,kz-1) 235 dz2=uvwzlev(ix,jy,kz)-height(iz) 298 236 dz=dz1+dz2 299 237 uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz … … 339 277 340 278 drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ & 341 279 (height(2)-height(1)) 342 280 do kz=2,nz-1 343 281 drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ & 344 282 (height(kz+1)-height(kz-1)) 345 283 end do 346 284 drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n) … … 356 294 357 295 do jy=1,ny-2 296 cosf=cos((real(jy)*dy+ylat0)*pi180) 358 297 do ix=1,nx-2 359 298 … … 361 300 do iz=2,nz-1 362 301 363 ui=uu(ix,jy,iz,n)*dxconst/cos ((real(jy)*dy+ylat0)*pi180)302 ui=uu(ix,jy,iz,n)*dxconst/cosf 364 303 vi=vv(ix,jy,iz,n)*dyconst 365 304 366 305 do kz=kmin,nz 367 306 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. & 368 307 (height(iz).le.uvwzlev(ix,jy,kz))) then 369 308 dz1=height(iz)-uvwzlev(ix,jy,kz-1) 370 309 dz2=uvwzlev(ix,jy,kz)-height(iz) … … 409 348 do iz=1,nz 410 349 call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), & 411 vv(ix,jy,iz,n),uupol(ix,jy,iz,n), & 412 vvpol(ix,jy,iz,n)) 350 vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n)) 413 351 end do 414 352 end do … … 424 362 xlon=xlon0+real(nx/2-1)*dx 425 363 xlonr=xlon*pi/180. 426 ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ & 427 vv(nx/2-1,nymin1,iz,n)**2) 364 ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2) 428 365 if (vv(nx/2-1,nymin1,iz,n).lt.0.) then 429 ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ & 430 vv(nx/2-1,nymin1,iz,n))-xlonr 366 ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr 431 367 else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then 432 368 ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ & 433 369 vv(nx/2-1,nymin1,iz,n))-xlonr 434 370 else 435 371 ddpol=pi/2-xlonr … … 444 380 uuaux=-ffpol*sin(xlonr+ddpol) 445 381 vvaux=-ffpol*cos(xlonr+ddpol) 446 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, & 447 vvpolaux) 448 382 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux) 449 383 jy=nymin1 450 384 do ix=0,nxmin1 … … 485 419 do iz=1,nz 486 420 call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), & 487 vv(ix,jy,iz,n),uupol(ix,jy,iz,n), & 488 vvpol(ix,jy,iz,n)) 421 vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n)) 489 422 end do 490 423 end do … … 499 432 xlon=xlon0+real(nx/2-1)*dx 500 433 xlonr=xlon*pi/180. 501 ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ & 502 vv(nx/2-1,0,iz,n)**2) 434 ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2) 503 435 if (vv(nx/2-1,0,iz,n).lt.0.) then 504 ddpol=atan(uu(nx/2-1,0,iz,n)/ & 505 vv(nx/2-1,0,iz,n))+xlonr 436 ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr 506 437 else if (vv(nx/2-1,0,iz,n).gt.0.) then 507 ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ & 508 vv(nx/2-1,0,iz,n))+xlonr 438 ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr 509 439 else 510 440 ddpol=pi/2-xlonr … … 519 449 uuaux=+ffpol*sin(xlonr-ddpol) 520 450 vvaux=-ffpol*cos(xlonr-ddpol) 521 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, & 522 vvpolaux) 451 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux) 523 452 524 453 jy=0 … … 551 480 ! create a cloud and rainout/washout field, clouds occur where rh>80% 552 481 ! total cloudheight is stored at level 0 553 554 555 556 482 do jy=0,nymin1 557 483 do ix=0,nxmin1 558 559 560 561 ! rain_cloud_above=0 562 ! lsp=lsprec(ix,jy,1,n) 563 ! convp=convprec(ix,jy,1,n) 564 ! cloudsh(ix,jy,n)=0 565 ! do kz_inv=1,nz-1 566 ! kz=nz-kz_inv+1 567 ! pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n) 568 ! rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n)) 569 ! clouds(ix,jy,kz,n)=0 570 ! if (rh.gt.0.8) then ! in cloud 571 ! if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation 572 ! rain_cloud_above=1 573 ! cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ & 574 ! height(kz)-height(kz-1) 575 ! if (lsp.ge.convp) then 576 ! clouds(ix,jy,kz,n)=3 ! lsp dominated rainout 577 ! else 578 ! clouds(ix,jy,kz,n)=2 ! convp dominated rainout 579 ! endif 580 ! else ! no precipitation 581 ! clouds(ix,jy,kz,n)=1 ! cloud 582 ! endif 583 ! else ! no cloud 584 ! if (rain_cloud_above.eq.1) then ! scavenging 585 ! if (lsp.ge.convp) then 586 ! clouds(ix,jy,kz,n)=5 ! lsp dominated washout 587 ! else 588 ! clouds(ix,jy,kz,n)=4 ! convp dominated washout 589 ! endif 590 ! endif 591 ! endif 592 ! end do 593 594 595 ! PS 3012 596 597 lsp=lsprec(ix,jy,1,n) 598 convp=convprec(ix,jy,1,n) 599 prec=lsp+convp 600 if (lsp.gt.convp) then ! prectype='lsp' 601 lconvectprec = .false. 602 else ! prectype='cp ' 603 lconvectprec = .true. 604 endif 605 rhmin = 0.90 ! standard condition for presence of clouds 606 !PS note that original by Sabine Eckhart was 80% 607 !PS however, for T<-20 C we consider saturation over ice 608 !PS so I think 90% should be enough 609 icloudbot(ix,jy,n)=icmv 610 icloudtop=icmv ! this is just a local variable 611 98 do kz=1,nz 612 pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n) 613 rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n)) 614 !ps if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz) 615 if (rh .gt. rhmin) then 616 if (icloudbot(ix,jy,n) .eq. icmv) then 617 icloudbot(ix,jy,n)=nint(height(kz)) 618 endif 619 icloudtop=nint(height(kz)) ! use int to save memory 484 rain_cloud_above=0 485 lsp=lsprec(ix,jy,1,n) 486 convp=convprec(ix,jy,1,n) 487 cloudsh(ix,jy,n)=0 488 do kz_inv=1,nz-1 489 kz=nz-kz_inv+1 490 pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n) 491 rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n)) 492 clouds(ix,jy,kz,n)=0 493 if (rh.gt.0.8) then ! in cloud 494 if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation 495 rain_cloud_above=1 496 cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ & 497 height(kz)-height(kz-1) 498 if (lsp.ge.convp) then 499 clouds(ix,jy,kz,n)=3 ! lsp dominated rainout 500 else 501 clouds(ix,jy,kz,n)=2 ! convp dominated rainout 502 endif 503 else ! no precipitation 504 clouds(ix,jy,kz,n)=1 ! cloud 620 505 endif 621 enddo 622 623 !PS try to get a cloud thicker than 50 m 624 !PS if there is at least .01 mm/h - changed to 0.002 and put into 625 !PS parameter precpmin 626 if ((icloudbot(ix,jy,n) .eq. icmv .or. & 627 icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. & 628 prec .gt. precmin) then 629 rhmin = rhmin - 0.05 630 if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum. 631 endif 632 !PS implement a rough fix for badly represented convection 633 !PS is based on looking at a limited set of comparison data 634 if (lconvectprec .and. icloudtop .lt. 6000 .and. & 635 prec .gt. precmin) then 636 637 if (convp .lt. 0.1) then 638 icloudbot(ix,jy,n) = 500 639 icloudtop = 8000 640 else 641 icloudbot(ix,jy,n) = 0 642 icloudtop = 10000 506 else ! no cloud 507 if (rain_cloud_above.eq.1) then ! scavenging 508 if (lsp.ge.convp) then 509 clouds(ix,jy,kz,n)=5 ! lsp dominated washout 510 else 511 clouds(ix,jy,kz,n)=4 ! convp dominated washout 512 endif 643 513 endif 644 endif 645 if (icloudtop .ne. icmv) then 646 icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n) 647 else 648 icloudthck(ix,jy,n) = icmv 649 endif 650 !PS get rid of too thin clouds 651 if (icloudthck(ix,jy,n) .lt. 50) then 652 icloudbot(ix,jy,n)=icmv 653 icloudthck(ix,jy,n)=icmv 654 endif 655 656 514 endif 515 end do 657 516 end do 658 517 end do 659 518 660 !do 102 kz=1,nuvz661 !write(an,'(i02)') kz+10662 !write(*,*) nuvz,nymin1,nxmin1,'--',an,'--'663 !open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')664 !do 101 jy=0,nymin1665 ! write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)666 !101 continue667 ! close(4)668 !102 continue669 670 ! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')671 ! do 103 jy=0,nymin1672 ! write (4,*)673 !+ (height(kz),kz=1,nuvz)674 !103 continue675 ! close(4)676 677 !open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')678 ! do 104 jy=0,nymin1679 ! write (4,*)680 !+ (r_air*tt(ix,jy,1,n)*rho(ix,jy,1,n),ix=1,nxmin1)681 !104 continue682 ! close(4)683 684 685 519 end subroutine verttransform -
trunk/src/verttransform_gfs.f90
r4 r24 21 21 22 22 subroutine verttransform(n,uuh,vvh,wwh,pvh) 23 ! 23 ! i i i i i 24 24 !***************************************************************************** 25 25 ! * … … 71 71 real :: f_qvsat,pressure 72 72 real :: rh,lsp,convp 73 real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)73 real :: rhoh(nuvzmax),pinmconv(nzmax) 74 74 real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi 75 75 real :: xlon,ylat,xlonr,dzdx,dzdy 76 real :: dzdx1,dzdx2,dzdy1,dzdy2 76 real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf 77 77 real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy 78 78 real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax) … … 92 92 ! If verttransform is called the first time, initialize heights of the * 93 93 ! z levels in meter. The heights are the heights of model levels, where * 94 ! u,v,T and qv are given, and of the interfaces, where w is given. So, * 95 ! the vertical resolution in the z system is doubled. As reference point,* 96 ! the lower left corner of the grid is used. * 97 ! Unlike in the eta system, no difference between heights for u,v and * 98 ! heights for w exists. * 94 ! u,v,T and qv are given. * 99 95 !************************************************************************* 100 96 … … 118 114 119 115 tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ & 120 116 ps(ixm,jym,1,n)) 121 117 pold=ps(ixm,jym,1,n) 122 118 height(1)=0. … … 126 122 tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n)) 127 123 128 129 ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled130 ! upon the transformation to z levels. In order to save computer memory, this is131 ! not done anymore in the standard version. However, this option can still be132 ! switched on by replacing the following lines with those below, that are133 ! currently commented out.134 ! Note that two more changes are necessary in this subroutine below.135 ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.136 !*****************************************************************************137 138 124 if (abs(tv-tvold).gt.0.2) then 139 height(kz)= & 140 height(kz-1)+const*log(pold/pint)* & 141 (tv-tvold)/log(tv/tvold) 125 height(kz)=height(kz-1)+const*log(pold/pint)* & 126 (tv-tvold)/log(tv/tvold) 142 127 else 143 height(kz)=height(kz-1)+ & 144 const*log(pold/pint)*tv 128 height(kz)=height(kz-1)+const*log(pold/pint)*tv 145 129 endif 146 147 ! Switch on following lines to use doubled vertical resolution148 !*************************************************************149 ! if (abs(tv-tvold).gt.0.2) then150 ! height((kz-1)*2)=151 ! + height(max((kz-2)*2,1))+const*log(pold/pint)*152 ! + (tv-tvold)/log(tv/tvold)153 ! else154 ! height((kz-1)*2)=height(max((kz-2)*2,1))+155 ! + const*log(pold/pint)*tv156 ! endif157 ! End doubled vertical resolution158 130 159 131 tvold=tv 160 132 pold=pint 161 133 end do 162 163 164 ! Switch on following lines to use doubled vertical resolution165 !*************************************************************166 ! do 7 kz=3,nz-1,2167 ! height(kz)=0.5*(height(kz-1)+height(kz+1))168 ! height(nz)=height(nz-1)+height(nz-1)-height(nz-2)169 ! End doubled vertical resolution170 134 171 135 … … 198 162 llev = 0 199 163 do i=1,nuvz 200 if (ps(ix,jy,1,n).lt.akz(i)) llev=i164 if (ps(ix,jy,1,n).lt.akz(i)) llev=i 201 165 end do 202 166 llev = llev+1 … … 212 176 tvold=tth(ix,jy,llev,n)*(1.+0.608*qvh(ix,jy,llev,n)) 213 177 pold=akz(llev) 214 uvzlev(llev)=0.215 178 wzlev(llev)=0. 216 179 uvwzlev(ix,jy,llev)=0. … … 223 186 224 187 if (abs(tv-tvold).gt.0.2) then 225 uv zlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &226 188 uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* & 189 (tv-tvold)/log(tv/tvold) 227 190 else 228 uv zlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv191 uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv 229 192 endif 230 wzlev(kz)=uvzlev(kz) 231 uvwzlev(ix,jy,kz)=uvzlev(kz) 193 wzlev(kz)=uvwzlev(ix,jy,kz) 232 194 233 195 tvold=tv 234 196 pold=pint 235 197 end do 236 237 238 ! Switch on following lines to use doubled vertical resolution239 ! Switch off the three lines above.240 !*************************************************************241 !22 uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)242 ! do 23 kz=2,nwz243 !23 uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)244 ! End doubled vertical resolution245 198 246 199 ! pinmconv=(h2-h1)/(p2-p1) … … 279 232 do iz=2,nz-1 280 233 do kz=kmin,nuvz 281 if(height(iz).gt.uv zlev(nuvz)) then234 if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then 282 235 uu(ix,jy,iz,n)=uu(ix,jy,nz,n) 283 236 vv(ix,jy,iz,n)=vv(ix,jy,nz,n) … … 289 242 goto 30 290 243 endif 291 if ((height(iz).gt.uv zlev(kz-1)).and. &292 (height(iz).le.uvzlev(kz))) then293 dz1=height(iz)-uvzlev(kz-1)294 dz2=uvzlev(kz)-height(iz)295 dz=dz1+dz2296 uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz297 vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz298 tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &299 300 qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &301 302 pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz303 rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz304 pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz244 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. & 245 (height(iz).le.uvwzlev(ix,jy,kz))) then 246 dz1=height(iz)-uvwzlev(ix,jy,kz-1) 247 dz2=uvwzlev(ix,jy,kz)-height(iz) 248 dz=dz1+dz2 249 uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz 250 vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz 251 tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 & 252 +tth(ix,jy,kz,n)*dz1)/dz 253 qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 & 254 +qvh(ix,jy,kz,n)*dz1)/dz 255 pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz 256 rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz 257 pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz 305 258 endif 306 259 end do … … 318 271 do kz=kmin,nwz 319 272 if ((height(iz).gt.wzlev(kz-1)).and. & 320 (height(iz).le.wzlev(kz))) then 321 dz1=height(iz)-wzlev(kz-1) 322 dz2=wzlev(kz)-height(iz) 323 dz=dz1+dz2 324 ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 & 325 +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz 326 273 (height(iz).le.wzlev(kz))) then 274 dz1=height(iz)-wzlev(kz-1) 275 dz2=wzlev(kz)-height(iz) 276 dz=dz1+dz2 277 ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 & 278 +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz 327 279 endif 328 280 end do … … 337 289 do kz=2,nz-1 338 290 drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ & 339 291 (height(kz+1)-height(kz-1)) 340 292 end do 341 293 drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n) … … 351 303 352 304 do jy=1,ny-2 305 cosf=cos((real(jy)*dy+ylat0)*pi180) 353 306 do ix=1,nx-2 354 307 … … 367 320 do iz=2,nz-1 368 321 369 ui=uu(ix,jy,iz,n)*dxconst/cos ((real(jy)*dy+ylat0)*pi180)322 ui=uu(ix,jy,iz,n)*dxconst/cosf 370 323 vi=vv(ix,jy,iz,n)*dyconst 371 324 372 325 do kz=kmin,nz 373 326 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. & 374 327 (height(iz).le.uvwzlev(ix,jy,kz))) then 375 328 dz1=height(iz)-uvwzlev(ix,jy,kz-1) 376 329 dz2=uvwzlev(ix,jy,kz)-height(iz) … … 426 379 xlon=xlon0+real(nx/2-1)*dx 427 380 xlonr=xlon*pi/180. 428 ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ & 429 vv(nx/2-1,nymin1,iz,n)**2) 430 if(vv(nx/2-1,nymin1,iz,n).lt.0.) then 431 ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ & 432 vv(nx/2-1,nymin1,iz,n))-xlonr 381 ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2) 382 if (vv(nx/2-1,nymin1,iz,n).lt.0.) then 383 ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr 433 384 elseif (vv(nx/2-1,nymin1,iz,n).gt.0.) then 434 385 ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ & 435 386 vv(nx/2-1,nymin1,iz,n))-xlonr 436 387 else 437 388 ddpol=pi/2-xlonr … … 446 397 uuaux=-ffpol*sin(xlonr+ddpol) 447 398 vvaux=-ffpol*cos(xlonr+ddpol) 448 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, & 449 vvpolaux) 450 399 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux) 451 400 jy=nymin1 452 401 do ix=0,nxmin1 … … 460 409 ! ward parallel of latitude 461 410 462 do iz=1,nz411 do iz=1,nz 463 412 wdummy=0. 464 413 jy=ny-2 … … 471 420 ww(ix,jy,iz,n)=wdummy 472 421 end do 473 end do422 end do 474 423 475 424 endif … … 487 436 do iz=1,nz 488 437 call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), & 489 vv(ix,jy,iz,n),uupol(ix,jy,iz,n), & 490 vvpol(ix,jy,iz,n)) 438 vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n)) 491 439 end do 492 440 end do … … 498 446 xlon=xlon0+real(nx/2-1)*dx 499 447 xlonr=xlon*pi/180. 500 ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ & 501 vv(nx/2-1,0,iz,n)**2) 448 ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2) 502 449 if(vv(nx/2-1,0,iz,n).lt.0.) then 503 ddpol=atan(uu(nx/2-1,0,iz,n)/ & 504 vv(nx/2-1,0,iz,n))+xlonr 450 ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr 505 451 elseif (vv(nx/2-1,0,iz,n).gt.0.) then 506 ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ & 507 vv(nx/2-1,0,iz,n))-xlonr 452 ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))-xlonr 508 453 else 509 454 ddpol=pi/2-xlonr … … 518 463 uuaux=+ffpol*sin(xlonr-ddpol) 519 464 vvaux=-ffpol*cos(xlonr-ddpol) 520 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, & 521 vvpolaux) 465 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux) 522 466 523 467 jy=0 … … 562 506 clouds(ix,jy,kz,n)=0 563 507 if (rh.gt.0.8) then ! in cloud 564 if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation 565 rain_cloud_above=1 566 cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ & 567 height(kz)-height(kz-1) 568 if (lsp.ge.convp) then 569 clouds(ix,jy,kz,n)=3 ! lsp dominated rainout 570 else 571 clouds(ix,jy,kz,n)=2 ! convp dominated rainout 572 endif 573 else ! no precipitation 574 clouds(ix,jy,kz,n)=1 ! cloud 575 endif 508 if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation 509 rain_cloud_above=1 510 cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1) 511 if (lsp.ge.convp) then 512 clouds(ix,jy,kz,n)=3 ! lsp dominated rainout 513 else 514 clouds(ix,jy,kz,n)=2 ! convp dominated rainout 515 endif 516 else ! no precipitation 517 clouds(ix,jy,kz,n)=1 ! cloud 518 endif 576 519 else ! no cloud 577 578 579 580 581 582 583 520 if (rain_cloud_above.eq.1) then ! scavenging 521 if (lsp.ge.convp) then 522 clouds(ix,jy,kz,n)=5 ! lsp dominated washout 523 else 524 clouds(ix,jy,kz,n)=4 ! convp dominated washout 525 endif 526 endif 584 527 endif 585 528 end do -
trunk/src/verttransform_nests.f90
r4 r24 21 21 22 22 subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn) 23 ! 23 ! i i i i i 24 24 !***************************************************************************** 25 25 ! * … … 72 72 integer :: rain_cloud_above,kz_inv 73 73 real :: f_qvsat,pressure,rh,lsp,convp 74 real :: uvzlev(nuvzmax),wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax)74 real :: wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax) 75 75 real :: uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax) 76 real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi 76 real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi,cosf 77 77 real :: dzdx,dzdy 78 78 real :: dzdx1,dzdx2,dzdy1,dzdy2 … … 96 96 97 97 tvold=tt2n(ix,jy,1,n,l)*(1.+0.378*ew(td2n(ix,jy,1,n,l))/ & 98 98 psn(ix,jy,1,n,l)) 99 99 pold=psn(ix,jy,1,n,l) 100 uv zlev(1)=0.100 uvwzlev(ix,jy,1)=0. 101 101 wzlev(1)=0. 102 102 rhoh(1)=pold/(r_air*tvold) … … 112 112 113 113 if (abs(tv-tvold).gt.0.2) then 114 uv zlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &115 114 uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* & 115 (tv-tvold)/log(tv/tvold) 116 116 else 117 uv zlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv117 uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv 118 118 endif 119 119 … … 124 124 125 125 do kz=2,nwz-1 126 wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2. 127 end do 128 wzlev(nwz)=wzlev(nwz-1)+ & 129 uvzlev(nuvz)-uvzlev(nuvz-1) 130 131 ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled 132 ! upon the transformation to z levels. In order to save computer memory, this is 133 ! not done anymore in the standard version. However, this option can still be 134 ! switched on by replacing the following lines with those below, that are 135 ! currently commented out. 136 ! Note that one change is also necessary in gridcheck.f, 137 ! and three changes in verttransform.f 138 !***************************************************************************** 139 uvwzlev(ix,jy,1)=0.0 140 do kz=2,nuvz 141 uvwzlev(ix,jy,kz)=uvzlev(kz) 142 end do 143 144 ! Switch on following lines to use doubled vertical resolution 145 ! Switch off the three lines above. 146 !************************************************************* 147 !22 uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz) 148 ! do 23 kz=2,nwz 149 !23 uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz) 150 ! End doubled vertical resolution 126 wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2. 127 end do 128 wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1) 151 129 152 130 ! pinmconv=(h2-h1)/(p2-p1) 153 131 154 132 pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ & 155 156 133 ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- & 134 (aknew(1)+bknew(1)*psn(ix,jy,1,n,l))) 157 135 do kz=2,nz-1 158 136 pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ & 159 160 137 ((aknew(kz+1)+bknew(kz+1)*psn(ix,jy,1,n,l))- & 138 (aknew(kz-1)+bknew(kz-1)*psn(ix,jy,1,n,l))) 161 139 end do 162 140 pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ & 163 164 141 ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- & 142 (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l))) 165 143 166 144 … … 183 161 do iz=2,nz-1 184 162 do kz=kmin,nuvz 185 if(height(iz).gt.uv zlev(nuvz)) then163 if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then 186 164 uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l) 187 165 vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l) … … 192 170 goto 30 193 171 endif 194 if ((height(iz).gt.uv zlev(kz-1)).and. &195 (height(iz).le.uvzlev(kz))) then196 dz1=height(iz)-uv zlev(kz-1)197 dz2=uv zlev(kz)-height(iz)172 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. & 173 (height(iz).le.uvwzlev(ix,jy,kz))) then 174 dz1=height(iz)-uvwzlev(ix,jy,kz-1) 175 dz2=uvwzlev(ix,jy,kz)-height(iz) 198 176 dz=dz1+dz2 199 177 uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ & 200 178 uuhn(ix,jy,kz,l)*dz1)/dz 201 179 vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ & 202 180 vvhn(ix,jy,kz,l)*dz1)/dz 203 181 ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ & 204 182 tthn(ix,jy,kz,n,l)*dz1)/dz 205 183 qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ & 206 184 qvhn(ix,jy,kz,n,l)*dz1)/dz 207 185 pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ & 208 186 pvhn(ix,jy,kz,l)*dz1)/dz 209 187 rhon(ix,jy,iz,n,l)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz 210 188 kmin=kz … … 225 203 do kz=kmin,nwz 226 204 if ((height(iz).gt.wzlev(kz-1)).and. & 227 228 dz1=height(iz)-wzlev(kz-1)229 dz2=wzlev(kz)-height(iz)230 dz=dz1+dz2231 wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &232 233 kmin=kz234 goto 40205 (height(iz).le.wzlev(kz))) then 206 dz1=height(iz)-wzlev(kz-1) 207 dz2=wzlev(kz)-height(iz) 208 dz=dz1+dz2 209 wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 & 210 +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz 211 kmin=kz 212 goto 40 235 213 endif 236 214 end do … … 242 220 243 221 drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ & 244 222 (height(2)-height(1)) 245 223 do kz=2,nz-1 246 224 drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- & 247 225 rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1)) 248 226 end do 249 227 drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l) … … 259 237 260 238 do jy=1,nyn(l)-2 239 cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180) 261 240 do ix=1,nxn(l)-2 262 241 … … 264 243 do iz=2,nz-1 265 244 266 ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/ & 267 cos((real(jy)*dyn(l)+ylat0n(l))*pi180) 245 ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf 268 246 vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l) 269 247 … … 319 297 rain_cloud_above=1 320 298 cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ & 321 299 height(kz)-height(kz-1) 322 300 if (lsp.ge.convp) then 323 301 cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout -
trunk/src/wetdepo.f90
r20 r24 21 21 22 22 subroutine wetdepo(itime,ltsample,loutnext) 23 ! 23 ! i i i 24 24 !***************************************************************************** 25 25 ! * … … 39 39 ! Code may not be correct for decay of deposition! * 40 40 ! * 41 ! Modification by Sabine Eckhart to introduce a new in-/below-cloud42 ! scheme, not dated43 ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification44 41 !***************************************************************************** 45 42 ! * … … 74 71 75 72 integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy 76 integer :: ngrid,itage,nage,hz,il,interp_time, n , clouds_v !NIK scheme 77 integer :: kz !PS scheme 78 integer :: ks, kp, n1,n2, icbot,ictop, indcloud 79 integer :: scheme_number ! NIK==1, PS ==2 73 integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v 74 integer :: ks, kp 80 75 real :: S_i, act_temp, cl, cle ! in cloud scavenging 81 76 real :: clouds_h ! cloud height for the specific grid point 82 real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav , precsub,f77 real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav 83 78 real :: wetdeposit(maxspec),restmass 84 79 real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled … … 103 98 104 99 do jpart=1,numpart 100 105 101 if (itra1(jpart).eq.-999999999) goto 20 106 102 if(ldirect.eq.1)then … … 114 110 if (itage.lt.lage(nage)) goto 33 115 111 end do 116 33 112 33 continue 117 113 118 114 … … 123 119 do j=numbnests,1,-1 124 120 if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. & 125 121 (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then 126 122 ngrid=j 127 123 goto 23 128 124 endif 129 125 end do 130 23 126 23 continue 131 127 132 128 … … 151 147 interp_time=nint(itime-0.5*ltsample) 152 148 153 ! PS nest case still needs to be implemented!! 154 ! if (ngrid.eq.0) then 155 ! call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, & 156 ! 1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, & 157 ! memtime(1),memtime(2),interp_time,lsp,convp,cc) 158 call interpol_rain(lsprec,convprec,tcc, & 159 icloudbot,icloudthck,nxmax,nymax,1,nx,ny, & 160 memind,sngl(xtra1(jpart)),sngl(ytra1(jpart)),1,memtime(1), & 161 memtime(2),interp_time,lsp,convp,cc,icbot,ictop,icmv) 162 ! else 163 ! call interpol_rain_nests(lsprecn,convprecn,tccn, & 164 ! nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, & 165 ! memtime(1),memtime(2),interp_time,lsp,convp,cc) 166 ! endif 167 168 169 ! if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20 170 !PS 2012: subtract a small value, eg 0.01 mm/h, to remove spurious precip 171 prec = lsp+convp 172 precsub = 0.01 173 if (prec .lt. precsub) then 174 goto 20 175 else 176 f = (prec-precsub)/prec 177 lsp = f*lsp 178 convp = f*convp 179 endif 180 181 149 if (ngrid.eq.0) then 150 call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, & 151 1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, & 152 memtime(1),memtime(2),interp_time,lsp,convp,cc) 153 else 154 call interpol_rain_nests(lsprecn,convprecn,tccn, & 155 nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, & 156 memtime(1),memtime(2),interp_time,lsp,convp,cc) 157 endif 158 159 if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20 182 160 ! get the level were the actual particle is in 183 do il=2,nz 184 if (height(il).gt.ztra1(jpart)) then 185 !hz=il-1 186 kz=il-1 187 goto 26 188 endif 189 end do 190 26 continue 191 192 n=memind(2) 193 if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) & 194 n=memind(1) 161 do il=2,nz 162 if (height(il).gt.ztra1(jpart)) then 163 hz=il-1 164 goto 26 165 endif 166 end do 167 26 continue 168 169 n=memind(2) 170 if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) & 171 n=memind(1) 195 172 196 173 ! if there is no precipitation or the particle is above the clouds no 197 174 ! scavenging is done 198 175 199 !old scheme 200 ! if (ngrid.eq.0) then 201 ! clouds_v=clouds(ix,jy,hz,n) 202 ! clouds_h=cloudsh(ix,jy,n) 203 ! else 204 ! clouds_v=cloudsn(ix,jy,hz,n,ngrid) 205 ! clouds_h=cloudsnh(ix,jy,n,ngrid) 206 ! endif 207 ! !write(*,*) 'there is 208 ! ! + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz 209 ! if (clouds_v.le.1) goto 20 210 ! !write (*,*) 'there is scavenging' 211 212 ! PS: part of 2011/2012 fix 213 214 if (ztra1(jpart) .le. float(ictop)) then 215 if (ztra1(jpart) .gt. float(icbot)) then 216 indcloud = 2 ! in-cloud 217 else 218 indcloud = 1 ! below-cloud 219 endif 220 elseif (ictop .eq. icmv) then 221 indcloud = 0 ! no cloud found, use old scheme 222 else 223 goto 20 ! above cloud 224 endif 225 176 if (ngrid.eq.0) then 177 clouds_v=clouds(ix,jy,hz,n) 178 clouds_h=cloudsh(ix,jy,n) 179 else 180 clouds_v=cloudsn(ix,jy,hz,n,ngrid) 181 clouds_h=cloudsnh(ix,jy,n,ngrid) 182 endif 183 !write(*,*) 'there is 184 ! + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz 185 if (clouds_v.le.1) goto 20 186 !write (*,*) 'there is scavenging' 226 187 227 188 ! 1) Parameterization of the the area fraction of the grid cell where the … … 267 228 !********************************************************** 268 229 269 do ks=1,nspec ! loop over species 270 wetdeposit(ks)=0. 271 272 273 !conflicting changes to the same routine: 1=NIK 2 =PS 274 scheme_number=2 275 if (scheme_number.eq.1) then !NIK 276 277 if (weta(ks).gt.0.) then 278 if (clouds_v.ge.4) then 279 ! BELOW CLOUD SCAVENGING 280 ! for aerosols and not highliy soluble substances weta=5E-6 281 wetscav=weta(ks)*prec**wetb(ks) ! scavenging coeff. 282 ! write(*,*) 'bel. wetscav: ',wetscav 283 284 else ! below_cloud clouds_v is lt 4 and gt 1 -> in cloud scavenging 285 ! IN CLOUD SCAVENGING 286 ! BUGFIX tt for nested fields should be ttn 287 ! sec may 2008 288 if (ngrid.gt.0) then 289 act_temp=ttn(ix,jy,hz,n,ngrid) 290 else 291 act_temp=tt(ix,jy,hz,n) 292 endif 230 do ks=1,nspec ! loop over species 231 wetdeposit(ks)=0. 232 wetscav=0. 233 234 ! NIK 09.12.13: allowed to turn off either below cloud or in-cloud by including TWO separate blocks of if tests 235 236 237 !****i******************* 238 ! BELOW CLOUD SCAVENGING 239 !*********************** 240 if (weta(ks).gt.0. .and. clouds_v.ge.4) then !if positive below-cloud coefficient given in SPECIES file, and if in below cloud regime 241 ! for aerosols and not highliy soluble substances weta=5E-6 242 wetscav=weta(ks)*prec**wetb(ks) ! scavenging coefficient 243 endif !end below-cloud 244 245 246 !******************** 247 ! IN CLOUD SCAVENGING 248 !******************** 249 if (weta_in(ks).gt.0. .and. clouds_v.lt.4) then !if positive in-cloud coefficient given in SPECIES file, and if in in-cloud regime 250 251 if (ngrid.gt.0) then 252 act_temp=ttn(ix,jy,hz,n,ngrid) 253 else 254 act_temp=tt(ix,jy,hz,n) 255 endif 293 256 294 257 ! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening … … 297 260 ! wetc_in=0.9 (default) 298 261 ! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling) 299 cl=weta_in(ks)*prec**wetb_in(ks) 300 if (dquer(ks).gt.0) then ! is particle 301 S_i=wetc_in(ks)/cl 302 else ! is gas 303 cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl 304 S_i=1/cle 305 endif 306 wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks) 307 ! write(*,*) 'in. wetscav:' 308 ! + ,wetscav,cle,cl,act_temp,prec,clouds_h 262 cl=weta_in(ks)*prec**wetb_in(ks) 263 if (dquer(ks).gt.0) then ! is particle 264 S_i=wetc_in(ks)/cl 265 else ! is gas 266 cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl 267 S_i=1/cle 309 268 endif 310 311 312 ! if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!' 313 ! + ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v 269 wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks) 270 ! write(*,*) 'in. wetscav:' 271 ! + ,wetscav,cle,cl,act_temp,prec,clouds_h 272 273 endif ! end in-cloud 274 275 276 277 !************************************************** 278 ! CALCULATE DEPOSITION 279 !************************************************** 280 ! if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!' 281 ! + ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v 282 283 if (wetscav.gt.0.) then 314 284 wetdeposit(ks)=xmass1(jpart,ks)* & 315 (1.-exp(-wetscav*abs(ltsample)))*grfraction ! wet deposition 316 ! new particle mass: 317 ! if (wetdeposit(ks).gt.0) then 318 ! write(*,*) 'wetdepo: ',wetdeposit(ks),ks 319 ! endif 320 restmass = xmass1(jpart,ks)-wetdeposit(ks) 321 if (ioutputforeachrelease.eq.1) then 322 kp=npoint(jpart) 323 else 324 kp=1 325 endif 326 if (restmass .gt. smallnum) then 327 xmass1(jpart,ks)=restmass 328 !ccccccccccccccc depostatistic 329 ! wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks) 330 !ccccccccccccccc depostatistic 331 else 332 xmass1(jpart,ks)=0. 333 endif 334 ! Correct deposited mass to the last time step when radioactive decay of 335 ! gridded deposited mass was calculated 336 if (decay(ks).gt.0.) then 337 wetdeposit(ks)=wetdeposit(ks) & 338 *exp(abs(ldeltat)*decay(ks)) 339 endif 340 else ! weta(k) 341 wetdeposit(ks)=0. 342 endif ! weta(k) 343 344 elseif (scheme_number.eq.2) then ! PS 345 346 !PS indcloud=0 ! Use this for FOR TESTING, 347 !PS will skip the new in/below cloud method !!! 348 349 if (weta(ks).gt.0.) then 350 if (indcloud .eq. 1) then ! BELOW CLOUD SCAVENGING 351 !C for aerosols and not highliy soluble substances weta=5E-6 352 wetscav=weta(ks)*prec**wetb(ks) ! scavenging coeff. 353 !c write(*,*) 'bel. wetscav: ',wetscav 354 elseif (indcloud .eq. 2) then ! IN CLOUD SCAVENGING 355 if (ngrid.gt.0) then 356 act_temp=ttn(ix,jy,kz,n,ngrid) 357 else 358 act_temp=tt(ix,jy,kz,n) 359 endif 360 361 ! from NIK 362 ! weta_in=2.0E-07 (default) 363 ! wetb_in=0.36 (default) 364 ! wetc_in=0.9 (default) 365 366 367 cl=2E-7*prec**0.36 368 if (dquer(ks).gt.0) then ! is particle 369 S_i=0.9/cl 370 else ! is gas 371 cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl 372 S_i=1/cle 373 endif 374 wetscav=S_i*prec/3.6E6/(ictop-icbot) ! 3.6e6 converts mm/h to m/s 375 else ! PS: no cloud diagnosed, old scheme, 376 !CPS using with fixed a,b for simplicity, one may wish to change!! 377 wetscav = 1.e-4*prec**0.62 378 endif 379 380 381 wetdeposit(ks)=xmass1(jpart,ks)* & 382 ! (1.-exp(-wetscav*abs(ltsample)))*fraction ! wet deposition 383 (1.-exp(-wetscav*abs(ltsample)))*grfraction ! fraction = grfraction (PS) 384 restmass = xmass1(jpart,ks)-wetdeposit(ks) 385 if (ioutputforeachrelease.eq.1) then 386 kp=npoint(jpart) 387 else 388 kp=1 389 endif 390 if (restmass .gt. smallnum) then 391 xmass1(jpart,ks)=restmass 392 !cccccccccccccccc depostatistic 393 !c wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks) 394 !cccccccccccccccc depostatistic 395 else 396 xmass1(jpart,ks)=0. 397 endif 398 !C Correct deposited mass to the last time step when radioactive decay of 399 !C gridded deposited mass was calculated 400 if (decay(ks).gt.0.) then 401 wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks)) 402 endif 403 else ! weta(k)<0 404 wetdeposit(ks)=0. 405 endif 406 407 endif !on scheme 408 409 410 411 end do 412 413 ! Sabine Eckhard, June 2008 create deposition runs only for forward runs 285 (1.-exp(-wetscav*abs(ltsample)))*grfraction ! wet deposition 286 else ! if no scavenging 287 wetdeposit(ks)=0. 288 endif 289 290 291 restmass = xmass1(jpart,ks)-wetdeposit(ks) 292 if (ioutputforeachrelease.eq.1) then 293 kp=npoint(jpart) 294 else 295 kp=1 296 endif 297 if (restmass .gt. smallnum) then 298 xmass1(jpart,ks)=restmass 299 ! depostatistic 300 ! wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks) 301 ! depostatistic 302 else 303 xmass1(jpart,ks)=0. 304 endif 305 ! Correct deposited mass to the last time step when radioactive decay of 306 ! gridded deposited mass was calculated 307 if (decay(ks).gt.0.) then 308 wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks)) 309 endif 310 311 312 313 end do !all species 314 315 ! Sabine Eckhardt, June 2008 create deposition runs only for forward runs 414 316 ! Add the wet deposition to accumulated amount on output grid and nested output grid 415 317 !***************************************************************************** 416 318 417 ! if (ldirect.eq.1) then 418 ! call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), & 419 ! real(ytra1(jpart)),nage,kp) 420 ! if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), & 421 ! wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), & 422 ! nage,kp) 423 ! endif 424 425 !PS 426 if (ldirect.eq.1) then 427 call wetdepokernel(nclass(jpart),wetdeposit, & 428 sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp) 429 if (nested_output.eq.1) & 430 call wetdepokernel_nest(nclass(jpart),wetdeposit, & 431 sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp) 432 endif 433 319 if (ldirect.eq.1) then 320 call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), & 321 real(ytra1(jpart)),nage,kp) 322 if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), & 323 wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp) 324 endif 434 325 435 326 20 continue 436 end do 327 end do ! all particles 437 328 438 329 end subroutine wetdepo -
trunk/src/writeheader_nest.f90
r4 r24 67 67 68 68 if (ldirect.eq.1) then 69 write(unitheader) ibdate,ibtime, 'FLEXPART V8.2'69 write(unitheader) ibdate,ibtime, trim(flexversion) 70 70 else 71 write(unitheader) iedate,ietime, 'FLEXPART V8.2'71 write(unitheader) iedate,ietime, trim(flexversion) 72 72 endif 73 73
Note: See TracChangeset
for help on using the changeset viewer.