Changes in / [20:30]


Ignore:
Files:
424 added
1 deleted
57 edited

Legend:

Unmodified
Added
Removed

  • /trunk/options/COMMAND

    r20 r30  
    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
     81. __                3X, I2
     9    1       
     10   LDIRECT           1 FOR FORWARD SIMULATION, -1 FOR BACKWARD SIMULATION
     11
     122. ________ ______   3X, I8, 1X, I6
     13   20111210 000000
     14   YYYYMMDD HHMISS   BEGINNING DATE OF SIMULATION
     15
     163. ________ ______   3X, I8, 1X, I6
     17   20111210 120000
     18   YYYYMMDD HHMISS   ENDING DATE OF SIMULATION
     19
     204. _____             3X, I5
     21   10800
     22   SSSSS             OUTPUT EVERY SSSSS SECONDS
     23
     245. _____             3X, I5
     25   10800
     26   SSSSS             TIME AVERAGE OF OUTPUT (IN SSSSS SECONDS)
     27
     286. _____             3X, I5
     29     900
     30   SSSSS             SAMPLING RATE OF OUTPUT (IN SSSSS SECONDS)
     31
     327. _________         3X, I9
     33   999999999
     34   SSSSSSSSS         TIME CONSTANT FOR PARTICLE SPLITTING (IN SECONDS)
     35
     368. _____             3X, I5
     37     900
     38   SSSSS             SYNCHRONISATION INTERVAL OF FLEXPART (IN SECONDS)
     39
     409.  ---.--           4X, F6.4
     41     -5.0
     42    CTL              FACTOR, BY WHICH TIME STEP MUST BE SMALLER THAN TL
     43
     4410. ---              4X, I3
     45      4
     46    IFINE            DECREASE OF TIME STEP FOR VERTICAL MOTION BY FACTOR IFINE
     47
     4811. -                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
     5212. -                4X, I1
     53    0 
     54    IPOUT            PARTICLE DUMP: 0 NO, 1 EVERY OUTPUT INTERVAL, 2 ONLY AT END
     55
     5613. _                4X, I1
     57    1
     58    LSUBGRID         SUBGRID TERRAIN EFFECT PARAMETERIZATION: 1 YES, 0 NO
     59
     6014. _                4X, I1
     61    1
     62    LCONVECTION      CONVECTION: 1 YES, 0 NO
     63
     6415. _                4X, I1
     65    0
     66    LAGESPECTRA      AGE SPECTRA: 1 YES, 0 NO
     67
     6816. _                4X, I1
     69    0
     70    IPIN             CONTINUE SIMULATION WITH DUMPED PARTICLE DATA: 1 YES, 0 NO
     71
     7217. _               
     73    0                4X,I1
     74    IOFR             IOUTPUTFOREACHREL CREATE AN OUPUT FILE FOR EACH RELEASE LOCATION: 1 YES, 0 NO
     75
     7618. _                4X, I1
     77    0
     78    IFLUX            CALCULATE FLUXES: 1 YES, 0 NO
     79
     8019. _                4X, I1
     81    0
     82    MDOMAINFILL      DOMAIN-FILLING TRAJECTORY OPTION: 1 YES, 0 NO, 2 STRAT. O3 TRACER
     83
     8420. _                4X, I1
     85    1
     86    IND_SOURCE       1=MASS UNIT , 2=MASS MIXING RATIO UNIT
     87
     8821. _                4X, I1
     89    1
     90    IND_RECEPTOR     1=MASS UNIT , 2=MASS MIXING RATIO UNIT
     91
     9222. _                4X, I1
     93    0
     94    MQUASILAG        QUASILAGRANGIAN MODE TO TRACK INDIVIDUAL PARTICLES: 1 YES, 0 NO
     95
     9623. _                4X, I1
     97    0
     98    NESTED_OUTPUT    SHALL NESTED OUTPUT BE USED? 1 YES, 0 NO
     99
     10024. _                4X, I1
     101    2
     102    LINIT_COND       INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT
     103
     10425. _                4X, I1
     105    0
     106    SURF_ONLY        IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER
     107
     108
     1091. Simulation direction, 1 for forward, -1 for backward in time
     110        (consult Seibert and Frank, 2004 for backward runs)
     111
     1122. 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
     1163. Ending date and time of simulation. Same format as 3.
     117
     1184. Average concentrations are calculated every SSSSS seconds.
     119
     1205. The average concentrations are time averages of SSSSS seconds
     121   duration. If SSSSS is 0, instantaneous concentrations are outputted.
     122
     1236. The concentrations are sampled every SSSSS seconds to calculate the time
     124   average concentration. This period must be shorter than the averaging time.
     125
     1267. Time constant for particle splitting. Particles are split into two
     127   after SSSSS seconds, 2xSSSSS seconds, 4xSSSSS seconds, and so on.
     128
     1298. 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
     1339. 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
     13610.IFINE=Reduction factor for time step used for vertical wind
     137
     13811.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
     14312.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
     14813.Switch on/off subgridscale terrain parameterization (increase of
     149   mixing heights due to subgridscale orographic variations)
     150
     15114.Switch on/off the convection parameterization
     152
     15315.Switch on/off the calculation of age spectra: if yes, the file AGECLASSES
     154   must be available
     155
     15616. 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
     16017. 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
     16418. 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
     17019. 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
     17620. 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
     18121. IND_RECEPTOR switches between different units for concentrations at the receptor
     182          1=mass units (concentrations)
     183          2=mass mixing ratio units
     184
     18522. 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
     18923. NESTED_OUTPUT decides whether model output shall be made also for a nested
     190    output field (normally with higher resolution)
     191
     19224. 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
     19625. 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

    r20 r30  
    30300               NESTED_OUTPUT     SHALL NESTED OUTPUT BE USED? YES, 0 NO
    31312               LINIT_COND        INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT
     320               SURF_ONLY         IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER
    3233
    3334
    … …  
    117118    conditions shall be calculated and written to output files
    118119    0=no output, 1 or 2 determines in which units the initial conditions are provided.
     120
     12125. 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

    r20 r30  
    1111
    12122. ________ ______   3X, I8, 1X, I6
    13    20040720 000000
     13   20110310 000000
    1414   YYYYMMDD HHMISS   BEGINNING DATE OF SIMULATION
    1515
    16163. ________ ______   3X, I8, 1X, I6
    17    20040721 120000
     17   20110310 120000
    1818   YYYYMMDD HHMISS   ENDING DATE OF SIMULATION
    1919
    … …  
    101101    2
    102102    LINIT_COND       INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT
     103
     10425. _                4X, I1
     105    0
     106    SURF_ONLY        IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER
    103107
    104108
    … …  
    189193    conditions shall be calculated and written to output files
    190194    0=no output, 1 or 2 determines in which units the initial conditions are provided.
     195
     19625. 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

    r20 r30  
    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
    1814
    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
    3517
    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=========================================================================
     1920111210 000001
     20________ ______            i8,1x,i6 Beginning date and time of release
    5221
    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
     2220111210 090000
     23________ ______            i8,1x,i6 Ending date and time of release
    6924
    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
    8627
    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
    10330
    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
    12033
    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
    13736
    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
    15439
    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
    17148
    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
     491.0000E00
     50_.____E__                  e9.4  Total mass emitted
    18851
    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 
     52RELEASE_TEST1
     53________________________________________   character*40 comment
     54+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  • /trunk/options/SPECIES/SPECIES_001

    r20 r30  
    77TRACER             Tracer name
    88-999.9             Species half life
    9 -9.9E-09           Wet deposition - A
    10                    Wet deposition - B
     9-9.9E-09           Below cloud scavenging - A
     10                   Below cloud scavenging - 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)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_002

    r20 r30  
    77O3                 Tracer name
    88-999.9             Species half life
    9 -9.9E-09           Wet deposition - A
    10                    Wet deposition - B
     9-9.9E-09           Below cloud scavenging - A
     10                   Below cloud scavenging - 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)
    1115 1.5               Dry deposition (gases) - D
    12161.0e-02            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_003

    r20 r30  
    77NO                 Tracer name
    88-999.9             Species half life
    9  8.0E-06           Wet deposition - A
    10 0.62               Wet deposition - B
     9 8.0E-06           Below cloud scavenging - A
     100.62               Below cloud scavenging - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.9                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 1.2               Dry deposition (gases) - D
    12163.0e-03            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_004

    r20 r30  
    77NO2                Tracer name
    88-999.9             Species half life
    9  1.0E-05           Wet deposition - A
    10 0.62               Wet deposition - B
     9 1.0E-05           Below cloud scavenging - A
     100.62               Below cloud scavenging - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.9                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 1.6               Dry deposition (gases) - D
    12161.0e-02            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_005

    r20 r30  
    77HNO3               Tracer name
    88-999.9             Species half life
    9  5.0E-05           Wet deposition - A
    10 0.62               Wet deposition - B
     9 5.0E-05           Below cloud scavenging - A
     100.62               Below cloud scavenging - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.9                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 1.9               Dry deposition (gases) - D
    12161.0e+14            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_006

    r20 r30  
    77HNO2               Tracer name
    88-999.9             Species half life
    9 -9.9E-09           Wet deposition - A
    10                    Wet deposition - B
     9-9.9E-09           Below cloud scavenging - A
     10                   Below cloud scavenging - 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)
    1115 1.6               Dry deposition (gases) - D
    12161.0e+05            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_007

    r20 r30  
    77H2O2               Tracer name
    88-999.9             Species half life
    9  1.0E-04           Wet deposition - A
    10 0.62               Wet deposition - B
     9 1.0E-04           Below cloud scavenging - A
     100.62               Below cloud scavenging - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.9                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 1.4               Dry deposition (gases) - D
    12161.0e+05            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_008

    r20 r30  
    77SO2                Tracer name
    88-999.9             Species half life
    9 -9.9E-09           Wet deposition - A
    10 0.62               Wet deposition - B
     9-9.9E-09           Below cloud scavenging - A
     100.62               Below cloud scavenging - 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)
    1115 2.0               Dry deposition (gases) - D
    12161.0e+05            Dry deposition (gases) - Henrys const.
    … …  
    394316       1.463       1.162
    404417       1.426       1.141
    41 18       1.326       1.117 
     4518       1.326       1.117
    424619       1.225       1.109
    434720       1.082       1.095

  • /trunk/options/SPECIES/SPECIES_009

    r20 r30  
    77HCHO               Tracer name
    88-999.9             Species half life
    9 -9.9E-09           Wet deposition - A
    10                    Wet deposition - B
     9-9.9E-09           Below cloud scavenging - A
     10                   Below cloud scavenging - 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)
    1115 1.3               Dry deposition (gases) - D
    12166.0e+03            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_010

    r20 r30  
    77PAN                Tracer name
    88-999.9             Species half life
    9 -9.9E-09           Wet deposition - A
    10                    Wet deposition - B
     9-9.9E-09           Below cloud scavenging - A
     10                   Below cloud scavenging - 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)
    1115 2.6               Dry deposition (gases) - D
    12163.6e+00            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_011

    r20 r30  
    77NH3                Tracer name
    88-999.9             Species half life
    9  9.9E-05           Wet deposition - A
    10 0.62               Wet deposition - B
     9 9.9E-05           Below cloud scavenging - A
     100.62               Below cloud scavenging - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.9                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 1.1               Dry deposition (gases) - D
    12162.0e+14            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_012

    r20 r30  
    77SO4-aero           Tracer name
    88-999.9             Species half life
    9  5.0E-06           Wet deposition - A
    10 0.62               Wet deposition - B
     9 5.0E-06           Below cloud scavenging - A
     100.62               Below cloud scavenging - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.9                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216 1.0E-09           Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_013

    r20 r30  
    77NO3-aero           Tracer name
    88-999.9             Species half life
    9  5.0E-06           Wet deposition - A
    10 0.62               Wet deposition - B
     9 5.0E-06           Below cloud scavenging - A
     100.62               Below cloud scavenging - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.9                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_014

    r20 r30  
    77I2-131             Tracer name
    88691200.0           Species half life
    9  8.0E-05           Wet deposition - A
    10 0.62               Wet deposition - B
     9 8.0E-05           Below cloud scavenging - A
     100.62               Below cloud scavenging - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.9                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 2.7               Dry deposition (gases) - D
    12161.0e+05            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_015

    r20 r30  
    77I-131              Tracer name
    88691200.0           Species half life
    9  1.0E-04           Wet deposition - A
    10 0.80               Wet deposition - B
     9 1.0E-04           Below cloud scavenging - A
     100.80               Below cloud scavenging - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.9                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_016

    r20 r30  
    77Cs-137             Tracer name
    88-999.9             Species half life
    9  1.0E-04           Wet deposition - A
    10 0.80               Wet deposition - B
     9 1.0E-04           Below cloud scavenging - A
     100.80               Below cloud scavenging - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.9                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_017

    r20 r30  
    77Y-91               Tracer name
    885037120.0          Species half life
    9  1.0E-04           Wet deposition - A
    10 0.80               Wet deposition - B
     9 1.0E-04           Below cloud scavenging - A
     100.80               Below cloud scavenging - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.9                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_018

    r20 r30  
    77Ru-106             Tracer name
    8831536000.0         Species half life
    9  1.0E-04           Wet deposition - A
    10 0.80               Wet deposition - B
     9 1.0E-04           Below cloud scavenging - A
     100.80               Below cloud scavenging - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.9                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_019

    r20 r30  
    77Kr-85              Tracer name
    88-999.9             Species half life
    9 -9.9E-09           Wet deposition - A
    10                    Wet deposition - B
     9-9.9E-09           Below cloud scavenging - A
     10                   Below cloud scavenging - 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)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_020

    r20 r30  
    77Sr-90              Tracer name
    88-999.9             Species half life
    9  1.0E-04           Wet deposition - A
    10 0.80               Wet deposition - B
     9 1.0E-04           Below cloud scavenging - A
     100.80               Below cloud scavenging - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.9                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_021

    r20 r30  
    66****************************************************************************
    77Xe-133             Tracer name
    8 198720.0           Species half life
    9 -9.9E-09           Wet deposition - A
    10                    Wet deposition - B
     8453168.0           Species half life
     9-9.9E-09           Below cloud scavenging - A
     10                   Below cloud scavenging - 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)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_022

    r20 r30  
    77CO                 Tracer name
    88-999.9             Species half life
    9 -9.9E-09           Wet deposition - A
    10                    Wet deposition - B
     9-9.9E-09           Below cloud scavenging - A
     10                   Below cloud scavenging - 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)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_024

    r20 r30  
    77AIRTRACER             Tracer name
    88-999.9             Species half life
    9 -9.9E-09           Wet deposition - A
    10                    Wet deposition - B
     9-9.9E-09           Below cloud scavenging - A
     10                   Below cloud scavenging - 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)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_025

    r20 r30  
    992.0E-05            Wet deposition - A
    10100.8                Wet deposition - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.9                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • /trunk/src/FLEXPART.f90

    r20 r30  
    4545  use conv_mod
    4646
     47#ifdef NETCDF_OUTPUT
     48  use netcdf_output_mod, only: writeheader_netcdf
     49#endif
     50
     51
    4752  implicit none
    4853
    … …  
    6065  call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
    6166
    62   !
    63   flexversion='Version 9.1.8  (2013-12-08)'
    64   !verbosity=0
     67  ! FLEXPART version string
     68  flexversion='Version 9.2 beta (2014-07-01)'
     69  verbosity=0
     70
    6571  ! Read the pathnames where input/output files are stored
    6672  !*******************************************************
    … …  
    7682    call getarg(1,arg1)
    7783    pathfile=arg1
    78     verbosity=0
    7984    if (arg1(1:1).eq.'-') then
    80         write(pathfile,'(a11)') './pathnames'
    81         inline_options=arg1
     85      write(pathfile,'(a11)') './pathnames'
     86      inline_options=arg1
    8287    endif
    8388  case (0)
    8489    write(pathfile,'(a11)') './pathnames'
    85     verbosity=0
    8690  end select
    8791 
    88     if (inline_options(1:1).eq.'-') then
    89       print*, 'inline options=', inline_options       
    90       if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
    91          print*, 'verbose mode 1: additional information will be displayed'
    92          verbosity=1
    93       endif
    94       if (trim(inline_options).eq.'-v2') then
    95          print*, 'verbose mode 2: additional information will be displayed'
    96          verbosity=2
    97       endif
    98       if (trim(inline_options).eq.'-i') then
    99          print*, 'info mode: will provide run specific information and stop'
    100          verbosity=1
    101          info_flag=1
    102       endif
    103       if (trim(inline_options).eq.'-i2') then
    104          print*, 'info mode: will provide run specific information and stop'
    105          verbosity=2
    106          info_flag=1
    107       endif
    108     endif
    109 
    110 
    11192  ! Print the GPL License statement
    11293  !*******************************************************
    113   print*,'Welcome to FLEXPART', trim(flexversion)
    114   print*,'FLEXPART is free software released under the GNU Genera'// &
    115        'l Public License.'
    116            
    117   if (verbosity.gt.0) then
    118         WRITE(*,*) 'call readpaths'
     94  print*,'Welcome to FLEXPART ', trim(flexversion)
     95  print*,'FLEXPART is free software released under the GNU General Public License.'
     96 
     97  if (inline_options(1:1).eq.'-') then
     98    if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
     99       print*, 'Verbose mode 1: display detailed information during run'
     100       verbosity=1
     101    endif
     102    if (trim(inline_options).eq.'-v2') then
     103       print*, 'Verbose mode 2: display more detailed information during run'
     104       verbosity=2
     105    endif
     106    if (trim(inline_options).eq.'-i') then
     107       print*, 'Info mode: provide detailed run specific information and stop'
     108       verbosity=1
     109       info_flag=1
     110    endif
     111    if (trim(inline_options).eq.'-i2') then
     112       print*, 'Info mode: provide more detailed run specific information and stop'
     113       verbosity=2
     114       info_flag=1
     115    endif
     116  endif
     117           
     118  if (verbosity.gt.0) then
     119    write(*,*) 'call readpaths'
    119120  endif
    120121  call readpaths(pathfile)
    121    
    122122 
    123123  if (verbosity.gt.1) then !show clock info
    … …  
    131131  endif
    132132
    133 
    134133  ! Read the user specifications for the current model run
    135134  !*******************************************************
    136135
    137136  if (verbosity.gt.0) then
    138         WRITE(*,*) 'call readcommand'
     137    write(*,*) 'call readcommand'
    139138  endif
    140139  call readcommand
    141140  if (verbosity.gt.0) then
    142         WRITE(*,*) '    ldirect=', ldirect
    143         WRITE(*,*) '    ibdate,ibtime=',ibdate,ibtime
    144         WRITE(*,*) '    iedate,ietime=', iedate,ietime
    145         if (verbosity.gt.1) then   
    146                 CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    147                 WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
    148         endif     
    149   endif
    150 
     141    write(*,*) '    ldirect=', ldirect
     142    write(*,*) '    ibdate,ibtime=',ibdate,ibtime
     143    write(*,*) '    iedate,ietime=', iedate,ietime
     144    if (verbosity.gt.1) then   
     145      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
     146      write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
     147    endif     
     148  endif
    151149
    152150  ! Read the age classes to be used
    153151  !********************************
    154152  if (verbosity.gt.0) then
    155         WRITE(*,*) 'call readageclasses'
     153    write(*,*) 'call readageclasses'
    156154  endif
    157155  call readageclasses
    158156
    159157  if (verbosity.gt.1) then   
    160                 CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    161                 WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
     158    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
     159    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
    162160  endif     
    163  
    164 
    165161
    166162  ! Read, which wind fields are available within the modelling period
    … …  
    168164
    169165  if (verbosity.gt.0) then
    170         WRITE(*,*) 'call readavailable'
     166    write(*,*) 'call readavailable'
    171167  endif 
    172168  call readavailable
    … …  
    177173 
    178174  if (verbosity.gt.0) then
    179      WRITE(*,*) 'call gridcheck'
     175     write(*,*) 'call gridcheck'
    180176  endif
    181177
    … …  
    183179
    184180  if (verbosity.gt.1) then   
    185      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    186      WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
     181    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
     182    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
    187183  endif     
    188  
    189 
    190   if (verbosity.gt.0) then
    191         WRITE(*,*) 'call gridcheck_nests'
     184
     185  if (verbosity.gt.0) then
     186    write(*,*) 'call gridcheck_nests'
    192187  endif 
    193188  call gridcheck_nests
    194189
    195 
    196190  ! Read the output grid specifications
    197191  !************************************
    198192
    199193  if (verbosity.gt.0) then
    200         WRITE(*,*) 'call readoutgrid'
     194    write(*,*) 'call readoutgrid'
    201195  endif
    202196
    … …  
    204198
    205199  if (nested_output.eq.1) then
    206           call readoutgrid_nest
     200    call readoutgrid_nest
    207201    if (verbosity.gt.0) then
    208         WRITE(*,*) '# readoutgrid_nest'
     202      write(*,*) '# readoutgrid_nest'
    209203    endif
    210204  endif
    … …  
    232226  call readlanduse
    233227
    234 
    235228  ! Assign fractional cover of landuse classes to each ECMWF grid point
    236229  !********************************************************************
    … …  
    241234  call assignland
    242235
    243 
    244 
    245236  ! Read the coordinates of the release locations
    246237  !**********************************************
    … …  
    251242  call readreleases
    252243
    253 
    254244  ! Read and compute surface resistances to dry deposition of gases
    255245  !****************************************************************
    … …  
    267257    print*,'call coordtrafo'
    268258  endif
    269 
    270259
    271260  ! Initialize all particles to non-existent
    … …  
    295284  endif
    296285
    297 
    298286  ! Calculate volume, surface area, etc., of all output grid cells
    299287  ! Allocate fluxes and OHfield if necessary
    300288  !***************************************************************
    301289
    302 
    303290  if (verbosity.gt.0) then
    304291    print*,'call outgrid_init'
    … …  
    307294  if (nested_output.eq.1) call outgrid_init_nest
    308295
    309 
    310296  ! Read the OH field
    311297  !******************
    … …  
    313299  if (OHREA.eqv..TRUE.) then
    314300    if (verbosity.gt.0) then
    315        print*,'call readOHfield'
    316     endif
    317        call readOHfield
     301      print*,'call readOHfield'
     302    endif
     303    call readOHfield
    318304  endif
    319305
    … …  
    322308  !******************************************************************
    323309
     310  if (lnetcdfout.eq.1) then
     311#ifdef NETCDF_OUTPUT
     312    call writeheader_netcdf(lnest = .false.)
     313#endif
     314  else
     315    call writeheader
     316  end if
     317
     318  if (nested_output.eq.1) then
     319    if (lnetcdfout.eq.1) then
     320#ifdef NETCDF_OUTPUT
     321      call writeheader_netcdf(lnest = .true.)
     322#endif
     323    else
     324      call writeheader_nest
     325    endif
     326  endif
    324327
    325328  if (verbosity.gt.0) then
    … …  
    332335  !if (nested_output.eq.1) call writeheader_nest
    333336  if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
    334 
    335337  if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
    336338  if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
    337339
    338 
    339 
    340   !open(unitdates,file=path(2)(1:length(2))//'dates')
     340  if (lnetcdfout.ne.1) then
     341    open(unitdates,file=path(2)(1:length(2))//'dates')
     342  end if
    341343
    342344  if (verbosity.gt.0) then
    … …  
    346348  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj
    347349
    348 
    349350  ! Releases can only start and end at discrete times (multiples of lsynctime)
    350351  !***************************************************************************
    … …  
    354355  endif
    355356  do i=1,numpoint
    356     ireleasestart(i)=nint(real(ireleasestart(i))/ &
    357          real(lsynctime))*lsynctime
    358     ireleaseend(i)=nint(real(ireleaseend(i))/ &
    359          real(lsynctime))*lsynctime
     357    ireleasestart(i)=nint(real(ireleasestart(i))/real(lsynctime))*lsynctime
     358    ireleaseend(i)=nint(real(ireleaseend(i))/real(lsynctime))*lsynctime
    360359  end do
    361 
    362360
    363361  ! Initialize cloud-base mass fluxes for the convection scheme
    … …  
    381379  end do
    382380
    383 
    384381  ! Calculate particle trajectories
    385382  !********************************
    … …  
    387384  if (verbosity.gt.0) then
    388385     if (verbosity.gt.1) then   
    389        CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    390        WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
     386       call system_clock(count_clock, count_rate, count_max)
     387       write(*,*) 'System clock',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
    391388     endif
    392389     if (info_flag.eq.1) then
    393          print*, 'info only mode (stop)'   
    394          stop
     390       print*, 'Info only mode (stop)'   
     391       stop
    395392     endif
    396393     print*,'call timemanager'
    … …  
    399396  call timemanager
    400397
    401 
    402   write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
    403        &XPART MODEL RUN!'
     398  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLEXPART MODEL RUN!'
     399
     400  ! output wall time
     401  if (verbosity .gt. 0) then
     402    call system_clock(count_clock,count_rate)
     403    tins=(count_clock - count_clock0)/real(count_rate)
     404    print*,'Wall time ',tins,'s, ',tins/60,'min, ',tins/3600,'h.'
     405  endif
    404406
    405407end program flexpart

  • /trunk/src/advance.f90

    r20 r30  
    463463      dcwsave=dcwsave+vp*dt
    464464      zt=zt+w*dt*real(ldirect)
     465
     466      ! HSO/AL: Particle managed to go over highest level -> interpolation error in goto 700
     467      !          alias interpol_wind (division by zero)
     468      if (zt.ge.height(nz)) zt=height(nz)-100.*eps
    465469
    466470      if (zt.gt.h) then
    … …  
    699703  endif
    700704
     705  ! HSO/AL: Prevent particles from disappearing at the pole
     706  !******************************************************************
     707
     708  if ( yt.lt.0. ) then
     709    xt=mod(xt+180.,360.)
     710    yt=-yt
     711  else if ( yt.gt.real(nymin1) ) then
     712    xt=mod(xt+180.,360.)
     713    yt=2*real(nymin1)-yt
     714  endif
    701715
    702716  ! Check position: If trajectory outside model domain, terminate it
    … …  
    704718
    705719  if ((xt.lt.0.).or.(xt.ge.real(nxmin1)).or.(yt.lt.0.).or. &
    706        (yt.ge.real(nymin1))) then
     720       (yt.gt.real(nymin1))) then
    707721    nstop=3
    708722    return
    … …  
    859873  endif
    860874
     875  ! HSO/AL: Prevent particles from disappearing at the pole
     876  !******************************************************************
     877
     878  if ( yt.lt.0. ) then
     879    xt=mod(xt+180.,360.)
     880    yt=-yt
     881  else if ( yt.gt.real(nymin1) ) then
     882    xt=mod(xt+180.,360.)
     883    yt=2*real(nymin1)-yt
     884  endif
     885
    861886  ! Check position: If trajectory outside model domain, terminate it
    862887  !*****************************************************************
    863888
    864889  if ((xt.lt.0.).or.(xt.ge.real(nxmin1)).or.(yt.lt.0.).or. &
    865        (yt.ge.real(nymin1))) then
     890       (yt.gt.real(nymin1))) then
    866891    nstop=3
    867892    return

  • /trunk/src/calcpv.f90

    r20 r30  
    2121
    2222subroutine calcpv(n,uuh,vvh,pvh)
    23   !                  i  i   i   o
     23  !               i  i   i   o
    2424  !*****************************************************************************
    2525  !                                                                            *
    … …  
    4949  integer :: nlck
    5050  real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
    51   real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
    52   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)
    5353  real :: thup,thdn
    5454  real,parameter :: eps=1.e-5, p0=101325
    … …  
    6262  ! Loop over entire grid
    6363  !**********************
     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
    6473  do jy=0,nymin1
    6574    if (sglobal.and.jy.eq.0) goto 10
    … …  
    106115      end if
    107116      jux=jumpx
    108   ! Precalculate pressure values for efficiency
    109       do kl=1,nuvz
    110         ppml(kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
    111       end do
    112117  !
    113118  ! Loop over the vertical
    … …  
    115120
    116121      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)
    119123        klvrp=kl+1
    120124        klvrm=kl-1
    … …  
    125129        if (klvrp.gt.nuvz) klvrp=nuvz
    126130        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))
    132134
    133135  ! Compute vertical position at pot. temperature surface on subgrid
    … …  
    156158          kch=kch+1
    157159          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
    176173            endif
     174            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
     175            goto 20
     176          endif
    17717741        continue
    178178  ! Downward branch
    … …  
    181181          kch=kch+1
    182182          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
    200195            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
    202200  ! This section used when no values were found
    20320121      continue
    … …  
    236234          kch=kch+1
    237235          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
    254247            endif
     248            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
     249            goto 50
     250          endif
    25525171        continue
    256252  ! Downward branch
    … …  
    259255          kch=kch+1
    260256          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
    277268            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
    279273  ! This section used when no values were found
    28027451      continue
    281275  ! Must use uu at current level and lat. juy becomes smaller by 1
    282         uy(jj)=uuh(ix,jy,kl)
    283         juy=juy-1
     276          uy(jj)=uuh(ix,jy,kl)
     277          juy=juy-1
    284278  ! Otherwise OK
    28527950        continue

  • /trunk/src/calcpv_nests.f90

    r20 r30  
    2121
    2222subroutine calcpv_nests(l,n,uuhn,vvhn,pvhn)
    23   !                        i i  i    i    o
     23  !                     i i  i    i    o
    2424  !*****************************************************************************
    2525  !                                                                            *
    … …  
    4848  integer :: nlck,l
    4949  real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
    50   real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
    51   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)
    5252  real :: thup,thdn
    5353  real,parameter :: eps=1.e-5,p0=101325
    … …  
    6161  ! Loop over entire grid
    6262  !**********************
     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
    6372  do jy=0,nyn(l)-1
    6473    phi = (ylat0n(l) + jy * dyn(l)) * pi / 180.
    … …  
    8897      if (ix.eq.0.or.ix.eq.nxn(l)-1) jumpx=1
    8998      jux=jumpx
    90   ! Precalculate pressure values for efficiency
    91       do kl=1,nuvz
    92         ppml(kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
    93       end do
    9499  !
    95100  ! Loop over the vertical
    … …  
    97102
    98103      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)
    101105        klvrp=kl+1
    102106        klvrm=kl-1
    … …  
    107111        if (klvrp.gt.nuvz) klvrp=nuvz
    108112        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))
    114116
    115117  ! Compute vertical position at pot. temperature surface on subgrid
    … …  
    134136          kch=kch+1
    135137          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)
    140140
    141141      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    … …  
    158158          kch=kch+1
    159159          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)
    164162      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    165163           ((thdn.le.theta).and.(thup.ge.theta))) then
    … …  
    212210          kch=kch+1
    213211          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)
    218214      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    219215           ((thdn.le.theta).and.(thup.ge.theta))) then
    … …  
    235231          kch=kch+1
    236232          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)
    241235      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    242236           ((thdn.le.theta).and.(thup.ge.theta))) then

  • /trunk/src/com_mod.f90

    r20 r30  
    113113  ! lagespectra  1 if age spectra calculation switched on, 2 if not
    114114
     115  integer :: lnetcdfout
     116  ! lnetcdfout   1 for netcdf grid output, 0 if not. Set in COMMAND (namelist input)
    115117
    116118  integer :: nageclass,lage(maxageclass)
    … …  
    340342  real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2)
    341343  real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2)
    342   !scavenging NIK, PS
    343344  integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2)
    344345  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)
    347346
    348347
    … …  
    359358  !      rainout  conv/lsp dominated  2/3
    360359  !      washout  conv/lsp dominated  4/5
    361 ! PS 2013
    362 !c icloudbot (m)        cloud bottom height
    363 !c icloudthck (m)       cloud thickness     
    364360
    365361  ! pplev for the GFS version
    … …  
    688684
    689685  !********************
    690   ! Verbosity, testing flags
     686  ! Verbosity, testing flags, namelist I/O
    691687  !********************   
    692688  integer :: verbosity=0
    693689  integer :: info_flag=0
    694   INTEGER :: count_clock, count_clock0,  count_rate, count_max
     690  integer :: count_clock, count_clock0,  count_rate, count_max
     691  real    :: tins
     692  logical :: nmlout=.true.
    695693   
    696694

  • /trunk/src/gridcheck.f90

    r20 r30  
    8484  real(kind=4) :: xaux1,xaux2,yaux1,yaux2
    8585  real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
    86   integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
     86  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parID
    8787  !HSO  end
    8888  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
    9090
    9191  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
    … …  
    9898
    9999  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)
    102101  character(len=1) :: opt
    103102
    … …  
    172171  call grib_check(iret,gribFunction,gribErrorMsg)
    173172  call grib_get_int(igrib,'level',valSurf,iret)
     173  call grib_check(iret,gribFunction,gribErrorMsg)
     174  call grib_get_int(igrib,'paramId',parId,iret)
    174175  call grib_check(iret,gribFunction,gribErrorMsg)
    175176
    … …  
    193194  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
    194195    isec1(6)=135         ! indicatorOfParameter
     196  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
     197    isec1(6)=135         ! indicatorOfParameter
    195198  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    196199    isec1(6)=151         ! indicatorOfParameter
    … …  
    205208  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    206209    isec1(6)=141         ! indicatorOfParameter
    207   elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
     210  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
    208211    isec1(6)=164         ! indicatorOfParameter
    209   elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
     212  elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
    210213    isec1(6)=142         ! indicatorOfParameter
    211214  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    … …  
    215218  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    216219    isec1(6)=176         ! indicatorOfParameter
    217   elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
     220  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
    218221    isec1(6)=180         ! indicatorOfParameter
    219   elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
     222  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
    220223    isec1(6)=181         ! indicatorOfParameter
    221224  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    222225    isec1(6)=129         ! indicatorOfParameter
    223   elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
     226  elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
    224227    isec1(6)=160         ! indicatorOfParameter
    225228  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    … …  
    230233         parCat,parNum,typSurf
    231234  endif
     235  if(parId .ne. isec1(6) .and. parId .ne. 77) then
     236    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
     237!    stop
     238  endif
    232239
    233240  endif
    … …  
    265272
    266273  !HSO  get the second part of the grid dimensions only from GRiB1 messages
    267   if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
     274  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
    268275    call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
    269276         xaux2in,iret)
    … …  
    291298    dyconst=180./(dy*r_earth*pi)
    292299    gotGrid=1
    293 
    294300  ! Check whether fields are global
    295301  ! If they contain the poles, specify polar stereographic map
    … …  
    437443  !********************
    438444
    439   write(*,*)
    440   write(*,*)
    441   write(*,'(a,2i7)') '# of vertical levels in ECMWF data: ', &
     445  write(*,'(a,2i7)') ' Vertical levels in ECMWF data: ', &
    442446       nuvz+1,nwz
    443447  write(*,*)
    444   write(*,'(a)') 'Mother domain:'
    445   write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
     448  write(*,'(a)') ' Mother domain:'
     449  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
    446450       xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
    447   write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
     451  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range : ', &
    448452       ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
    449453  write(*,*)
    … …  
    467471    k=nlev_ec+1+numskip+i
    468472    akm(nwz-i+1)=zsec2(j)
     473  !   write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)
    469474    bkm(nwz-i+1)=zsec2(k)
    470475  end do
    … …  
    553558
    554559end subroutine gridcheck
     560

  • /trunk/src/gridcheck_gfs.f90

    r20 r30  
    423423  write(*,*)
    424424  write(*,*)
    425   write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', &
     425  write(*,'(a,2i7)') 'Vertical levels in NCEP data: ', &
    426426       nuvz,nwz
    427427  write(*,*)
    … …  
    429429  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
    430430       xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
    431   write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
     431  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range : ', &
    432432       ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
    433433  write(*,*)

  • /trunk/src/gridcheck_nests.f90

    r20 r30  
    4848  integer :: igrib
    4949  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
     50  integer :: parID !added by mc for making it consistent with new gridcheck.f90
    5051  integer :: gotGrib
    5152  !HSO  end
    … …  
    5556  real(kind=4) :: xaux1,xaux2,yaux1,yaux2
    5657  real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
     58  real :: conversion_factor !added by mc to make it consistent with new gridchek.f90
    5759
    5860  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
    … …  
    150152  call grib_get_int(igrib,'level',valSurf,iret)
    151153  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
    152156
    153157  !print*,discipl,parCat,parNum,typSurf,valSurf
    … …  
    170174  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
    171175    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    !                     ! " " " " " " " " " " " " " " " " " " " " " " " "  " " "
    172178  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    173179    isec1(6)=151         ! indicatorOfParameter
    … …  
    182188  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    183189    isec1(6)=141         ! indicatorOfParameter
    184   elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
     190  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
    185191    isec1(6)=164         ! indicatorOfParameter
    186   elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
     192 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
    187193    isec1(6)=142         ! indicatorOfParameter
    188194  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    … …  
    192198  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    193199    isec1(6)=176         ! indicatorOfParameter
    194   elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
     200  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
    195201    isec1(6)=180         ! indicatorOfParameter
    196   elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
     202  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
    197203    isec1(6)=181         ! indicatorOfParameter
    198204  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    199205    isec1(6)=129         ! indicatorOfParameter
    200   elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
     206  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
    201207    isec1(6)=160         ! indicatorOfParameter
    202208  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    … …  
    207213         parCat,parNum,typSurf
    208214  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
    209219
    210220  endif
    … …  
    237247
    238248  !HSO  get the second part of the grid dimensions only from GRiB1 messages
    239   if ((gribVer.eq.1).and.(gotGrib.eq.0)) then
    240     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
    241251         xaux1in,iret)
    242252    call grib_check(iret,gribFunction,gribErrorMsg)
    … …  
    254264    yaux1=yaux1in
    255265    yaux2=yaux2in
    256     if(xaux1.gt.180) xaux1=xaux1-360.0
    257     if(xaux2.gt.180) xaux2=xaux2-360.0
    258     if(xaux1.lt.-180) xaux1=xaux1+360.0
    259     if(xaux2.lt.-180) xaux2=xaux2+360.0
    260     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
    261271    xlon0n(l)=xaux1
    262272    ylat0n(l)=yaux1
    263273    dxn(l)=(xaux2-xaux1)/real(nxn(l)-1)
    264274    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!!!
    266276  endif ! ifield.eq.1
    267277
    … …  
    340350  !********************
    341351
    342   write(*,'(a,i2)') 'Nested domain #: ',l
    343   write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
     352  write(*,'(a,i2,a)') ' Nested domain ',l,':'
     353  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
    344354       xlon0n(l),' to ',xlon0n(l)+(nxn(l)-1)*dxn(l), &
    345355       '   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 : ', &
    347357       ylat0n(l),' to ',ylat0n(l)+(nyn(l)-1)*dyn(l), &
    348358       '   Grid distance: ',dyn(l)

  • /trunk/src/interpol_rain.f90

    r20 r30  
    2020!**********************************************************************
    2121
    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 
     22subroutine 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
    3326  !****************************************************************************
    3427  !                                                                           *
    … …  
    4740  !                                                                           *
    4841  !     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  !                                                                           *
    5343  !****************************************************************************
    5444  !                                                                           *
    … …  
    8070
    8171  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
    8573  real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2)
    8674  real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2)
    8775  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
    9378
    9479
    … …  
    142127         + p3*yy3(ix ,jyp,level,indexh) &
    143128         + p4*yy3(ixp,jyp,level,indexh)
    144 
    145 !CPS clouds:
    146         ip1=1
    147         ip2=1
    148         ip3=1
    149         ip4=1
    150         if (iy1(ix ,jy ,indexh) .eq. icmv) ip1=0
    151         if (iy1(ixp,jy ,indexh) .eq. icmv) ip2=0
    152         if (iy1(ix ,jyp,indexh) .eq. icmv) ip3=0
    153         if (iy1(ixp,jyp,indexh) .eq. icmv) ip4=0
    154         ipsum= ip1+ip2+ip3+ip4
    155         if (ipsum .eq. 0) then
    156           yi1(m)=icmv
    157         else
    158           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))/ipsum
    162         endif
    163        
    164         ip1=1
    165         ip2=1
    166         ip3=1
    167         ip4=1
    168         if (iy2(ix ,jy ,indexh) .eq. icmv) ip1=0
    169         if (iy2(ixp,jy ,indexh) .eq. icmv) ip2=0
    170         if (iy2(ix ,jyp,indexh) .eq. icmv) ip3=0
    171         if (iy2(ixp,jyp,indexh) .eq. icmv) ip4=0
    172         ipsum= ip1+ip2+ip3+ip4
    173         if (ipsum .eq. 0) then
    174           yi2(m)=icmv
    175         else
    176           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))/ipsum
    180         endif
    181 !CPS end clouds
    182 
    183129  end do
    184130
    … …  
    187133  ! 2.) Temporal interpolation (linear)
    188134  !************************************
    189 
    190       if (abs(itime) .lt. abs(itime1)) then
    191         print*,'interpol_rain.f90'
    192         print*,itime,itime1,itime2
    193         stop 'ITIME PROBLEM'
    194       endif
    195 
    196135
    197136  dt1=real(itime-itime1)
    … …  
    204143
    205144
    206 !PS clouds:
    207       intiy1=(yi1(1)*dt2 + yi1(2)*dt1)/dt
    208       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)/dt
    212       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) then
    216         intiy2 = intiy1 + intiy2 ! convert cloud thickness to cloud top
    217       else
    218         intiy1=icmv
    219         intiy2=icmv
    220       endif
    221 !PS end clouds
    222 
    223145end subroutine interpol_rain

  • /trunk/src/makefile

    r20 r30  
    11SHELL = /bin/bash
    2 TARGET = local
    3 WINDS=ecmwf
    4 #WINDS=gfs
    5 #WINDS=fnl
     2MAIN = FP_ecmwf_gfortran
    63
    74FC       = gfortran
     5INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
     6LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
     7LIBPATH2 =   /usr/lib/x86_64-linux-gnu/
    88FFLAGS   =   -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
    9 
    10 ifeq ($(TARGET),dmz)   
    11   # options for ganglia
    12   INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
    13   LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
    14   LIBPATH2 = /usr/lib/x86_64-linux-gnu/
    15   MAIN = FLEXPART_dmz_
    16 endif
    17 ifeq ($(TARGET),local)
    18   # local options       
    19   #libs_dir=/.../flexpart/libs/
    20   libs_dir=/Users/ignacio/flexpart/libs/
    21   INCPATH  = $(libs_dir)/grib_api-1.9.9_dir/include
    22   LIBPATH1 = $(libs_dir)/grib_api-1.9.9_dir/lib
    23   LIBPATH2 = $(libs_dir)/jasper_dir/lib
    24   MAIN = FLEXPART_local_
    25 endif
    26 
     9#FFLAGS   =   -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
    2710LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
    2811
    … …  
    3720OBJECTS = \
    3821writeheader.o  writeheader_txt.o   writeheader_surf.o       assignland.o\
    39                part0.o \
     22calcpar.o               part0.o \
    4023caldate.o               partdep.o \
    4124coordtrafo.o            psih.o \
    … …  
    5336getrc.o                 readreleases.o \
    5437getvdep.o               readspecies.o \
    55 interpol_misslev.o      \
    56 conccalc.o              \
     38interpol_misslev.o      readwind.o \
     39conccalc.o              richardson.o \
    5740concoutput.o  concoutput_surf.o          scalev.o \
    5841pbl_profile.o           readOHfield.o\
    5942juldate.o               timemanager.o \
    6043interpol_vdep.o         interpol_rain.o \
    61 partoutput.o \
     44verttransform.o         partoutput.o \
    6245hanna.o                 wetdepokernel.o \
    6346mean.o                  wetdepo.o \
    6447hanna_short.o           windalign.o \
     48obukhov.o               gridcheck.o \
    6549hanna1.o                initialize.o \
    66                            calcpar_nests.o \
     50                        gridcheck_nests.o \
     51readwind_nests.o        calcpar_nests.o \
    6752verttransform_nests.o   interpol_all_nests.o \
    6853interpol_wind_nests.o   interpol_misslev_nests.o \
    6954interpol_vdep_nests.o   interpol_rain_nests.o \
    70 getvdep_nests.o   gridcheck_nests.o \
    71 readwind_nests.o  \
     55getvdep_nests.o \
    7256readageclasses.o        readpartpositions.o \
    7357calcfluxes.o            fluxoutput.o \
    7458qvsat.o                 skplin.o \
     59convmix.o               calcmatrix.o \
    7560convect43c.o               redist.o \
    7661sort2.o                 distance.o \
    … …  
    9176
    9277
    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)
     78$(MAIN): $(MODOBJS) $(OBJECTS)
     79        $(FC) *.o -o $(MAIN) $(LDFLAGS)
    11180
    11281$(OBJECTS): $(MODOBJS)
    … …  
    11887        rm *.o *.mod
    11988
     89cleanall:
     90        rm *.o *.mod $(MAIN)

  • /trunk/src/outgrid_init.f90

    r20 r30  
    225225  !     maxspec,maxpointspec_act,nclassunc,maxageclass
    226226
    227   write (*,*) ' Allocating fields for nested and global output (x,y): ', &
    228        max(numxgrid,numxgridn),max(numygrid,numygridn)
     227  write (*,*) 'Allocating fields for global output (x,y): ', numxgrid,numygrid
     228  write (*,*) 'Allocating fields for nested output (x,y): ', numxgridn,numygridn
    229229
    230230  ! allocate fields for concoutput with maximum dimension of outgrid

  • /trunk/src/par_mod.f90

    r20 r30  
    121121  !*********************************************
    122122 
    123   integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !FNL
    124   !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF
     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
     125  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF
    125126  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26
    126127  !integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
    … …  
    154155
    155156  !integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
    156   !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151 !ECMWF
    157   integer,parameter :: maxnests=1, nxmaxn=201, nymaxn=161 ! FNL XF
     157  integer,parameter :: maxnests=1,nxmaxn=351,nymaxn=351 !ECMWF
     158  !integer,parameter :: maxnests=1, nxmaxn=201, nymaxn=161 ! FNL XF
    158159  ! maxnests                maximum number of nested grids
    159160  ! nxmaxn,nymaxn           maximum dimension of nested wind fields in
    … …  
    197198  !**************************************************
    198199
    199   !integer,parameter :: maxpart=4000000
    200   integer,parameter :: maxpart=2000
    201   integer,parameter :: maxspec=6
     200  integer,parameter :: maxpart=150000
     201  integer,parameter :: maxspec=4
    202202
    203203

  • /trunk/src/readageclasses.f90

    r20 r30  
    2828  !                                                                            *
    2929  !     Author: A. Stohl                                                       *
    30   !                                                                            *
    3130  !     20 March 2000                                                          *
     31  !     HSO, 1 July 2014                                                       *
     32  !     Added optional namelist input                                          *
    3233  !                                                                            *
    3334  !*****************************************************************************
    … …  
    4647  integer :: i
    4748
     49  ! namelist help variables
     50  integer :: readerror
     51
     52  ! namelist declaration
     53  namelist /ageclass/ &
     54    nageclass, &
     55    lage
     56
     57  nageclass=-1 ! preset to negative value to identify failed namelist input
    4858
    4959  ! If age spectra calculation is switched off, set number of age classes
    … …  
    5767  endif
    5868
    59 
    6069  ! If age spectra claculation is switched on,
    6170  ! open the AGECLASSSES file and read user options
    6271  !************************************************
    6372
    64   open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES', &
    65        status='old',err=999)
     73  open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES',form='formatted',status='old',err=999)
    6674
    67   do i=1,13
    68     read(unitageclasses,*)
    69   end do
    70   read(unitageclasses,*) nageclass
     75  ! try to read in as a namelist
     76  read(unitageclasses,ageclass,iostat=readerror)
     77  close(unitageclasses)
    7178
     79  if ((nageclass.lt.0).or.(readerror.ne.0)) then
     80    open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES',status='old',err=999)
     81    do i=1,13
     82      read(unitageclasses,*)
     83    end do
     84    read(unitageclasses,*) nageclass
     85    read(unitageclasses,*) lage(1)
     86    do i=2,nageclass
     87      read(unitageclasses,*) lage(i)
     88    end do
     89    close(unitageclasses)
     90  endif
     91
     92  ! write ageclasses file in namelist format to output directory if requested
     93  if (nmlout.eqv..true.) then
     94    open(unitageclasses,file=path(2)(1:length(2))//'AGECLASSES.namelist',err=1000)
     95    write(unitageclasses,nml=ageclass)
     96    close(unitageclasses)
     97  endif
    7298
    7399  if (nageclass.gt.maxageclass) then
    … …  
    80106  endif
    81107
    82   read(unitageclasses,*) lage(1)
    83108  if (lage(1).le.0) then
    84109    write(*,*) ' #### FLEXPART MODEL ERROR! AGE OF FIRST      #### '
    … …  
    89114
    90115  do i=2,nageclass
    91     read(unitageclasses,*) lage(i)
    92116    if (lage(i).le.lage(i-1)) then
    93117      write(*,*) ' #### FLEXPART MODEL ERROR! AGE CLASSES     #### '
    … …  
    105129  stop
    106130
     1311000  write(*,*) ' #### FLEXPART MODEL ERROR! FILE "AGECLASSES" #### '
     132  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     133  write(*,'(a)') path(2)(1:length(2))
     134  stop
     135
     136
    107137end subroutine readageclasses

  • /trunk/src/readavailable.f90

    r20 r30  
    123123
    124124  do k=1,numbnests
    125   print*,length(numpath+2*(k-1)+1),length(numpath+2*(k-1)+2),length(4),length(3)
    126   print*,path(numpath+2*(k-1)+2)(1:length(numpath+2*(k-1)+2))
     125  !print*,length(numpath+2*(k-1)+1),length(numpath+2*(k-1)+2),length(4),length(3)
     126  !print*,path(numpath+2*(k-1)+2)(1:length(numpath+2*(k-1)+2))
    127127    open(unitavailab,file=path(numpath+2*(k-1)+2) &
    128128         (1:length(numpath+2*(k-1)+2)),status='old',err=998)

  • /trunk/src/readcommand.f90

    r20 r30  
    2929  !                                                                            *
    3030  !     18 May 1996                                                            *
     31  !     HSO, 1 July 2014                                                       *
     32  !     Added optional namelist input                                          *
    3133  !                                                                            *
    3234  !*****************************************************************************
    … …  
    8082  character(len=50) :: line
    8183  logical :: old
    82   logical :: nmlout=.true. !.false.
    8384  integer :: readerror
    8485
    8586  namelist /command/ &
    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   
     87  ldirect, &
     88  ibdate,ibtime, &
     89  iedate,ietime, &
     90  loutstep, &
     91  loutaver, &
     92  loutsample, &
     93  itsplit, &
     94  lsynctime, &
     95  ctl, &
     96  ifine, &
     97  iout, &
     98  ipout, &
     99  lsubgrid, &
     100  lconvection, &
     101  lagespectra, &
     102  ipin, &
     103  ioutputforeachrelease, &
     104  iflux, &
     105  mdomainfill, &
     106  ind_source, &
     107  ind_receptor, &
     108  mquasilag, &
     109  nested_output, &
     110  linit_cond, &
     111  lnetcdfout, &
     112  surf_only   
    111113
    112114  ! Presetting namelist command
    113   ldirect=1
     115  ldirect=0
    114116  ibdate=20000101
    115117  ibtime=0
    … …  
    137139  nested_output=0
    138140  linit_cond=0
     141  lnetcdfout=0
    139142  surf_only=0
    140143
    … …  
    142145  ! Namelist input first: try to read as namelist file
    143146  !**************************************************************************
    144   open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
    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 
     147  open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old',form='formatted',err=999)
     148
     149  ! try namelist input (default)
    156150  read(unitcommand,command,iostat=readerror)
    157151  close(unitcommand)
    158152
    159   ! If error in namelist format, try to open with old input code
    160   if (readerror.ne.0) then
    161 
    162     open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
    163          err=999)
     153  ! distinguish namelist from fixed text input
     154  if ((readerror.ne.0).or.(ldirect.eq.0)) then ! parse as text file format
     155 
     156    open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', err=999)
    164157
    165158    ! Check the format of the COMMAND file (either in free format,
    … …  
    173166    if (index(line,'LDIRECT') .eq. 0) then
    174167      old = .false.
     168      write(*,*) 'COMMAND in old short format, please update to namelist format'
    175169    else
    176170      old = .true.
     171      write(*,*) 'COMMAND in old long format, please update to namelist format'
    177172    endif
    178173    rewind(unitcommand)
    179174
     175
    180176    ! Read parameters
    181177    !****************
    … …  
    183179    call skplin(7,unitcommand)
    184180    if (old) call skplin(1,unitcommand)
    185 
    186181    read(unitcommand,*) ldirect
    187182    if (old) call skplin(3,unitcommand)
    … …  
    231226    if (old) call skplin(3,unitcommand)
    232227    read(unitcommand,*) linit_cond
     228    if (old) call skplin(3,unitcommand)
     229    read(unitcommand,*) surf_only
    233230    close(unitcommand)
    234231
    … …  
    237234  ! write command file in namelist format to output directory if requested
    238235  if (nmlout.eqv..true.) then
    239     !open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist.out',status='new',err=1000)
    240236    open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist',err=1000)
    241237    write(unitcommand,nml=command)
    242238    close(unitcommand)
    243      ! open(unitheader,file=path(2)(1:length(2))//'header_nml',status='new',err=999)
    244      ! write(unitheader,NML=COMMAND)
    245      !close(unitheader)
    246239  endif
    247240
    … …  
    353346  endif
    354347
     348!  check for netcdf output switch (use for non-namelist input only!)
     349  if (iout.ge.8) then
     350     lnetcdfout = 1
     351     iout = iout - 8
     352#ifndef NETCDF_OUTPUT
     353     print*,'ERROR: netcdf output not activated during compile time but used in COMMAND file!'
     354     print*,'Please recompile with netcdf library or use standard output format.'
     355     stop
     356#endif
     357  endif
     358
    355359  ! Check whether a valid option for gridded model output has been chosen
    356360  !**********************************************************************
    … …  
    358362  if ((iout.lt.1).or.(iout.gt.5)) then
    359363    write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
    360     write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4, OR 5!          #### '
     364    write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4, OR 5 FOR       #### '
     365    write(*,*) ' #### STANDARD FLEXPART OUTPUT OR  9 - 13     #### '
     366    write(*,*) ' #### FOR NETCDF OUTPUT                       #### '
    361367    stop
    362368  endif
    … …  
    370376    stop
    371377  endif
    372 
    373378
    374379
    … …  
    592597
    5935981000   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "COMMAND"    #### '
    594        write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    595         write(*,'(a)') path(2)(1:length(1))
    596         stop
     599  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     600  write(*,'(a)') path(2)(1:length(1))
     601  stop
    597602end subroutine readcommand

  • /trunk/src/readoutgrid.f90

    r20 r30  
    2929  !                                                                            *
    3030  !     4 June 1996                                                            *
     31  !     HSO, 1 July 2014
     32  !     Added optional namelist input
    3133  !                                                                            *
    3234  !*****************************************************************************
    … …  
    5355  real,parameter :: eps=1.e-4
    5456
    55 
     57  ! namelist variables
     58  integer, parameter :: maxoutlev=500
     59  integer :: readerror
     60  real,allocatable, dimension (:) :: outheights
     61
     62  ! declare namelist
     63  namelist /outgrid/ &
     64    outlon0,outlat0, &
     65    numxgrid,numygrid, &
     66    dxout,dyout, &
     67    outheights
     68
     69  ! allocate large array for reading input
     70  allocate(outheights(maxoutlev),stat=stat)
     71  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheights'
     72
     73  ! helps identifying failed namelist input
     74  dxout=-1.0
     75  outheights=-1.0
    5676
    5777  ! Open the OUTGRID file and read output grid specifications
    5878  !**********************************************************
    5979
    60   open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old', &
    61        err=999)
    62 
    63 
    64   call skplin(5,unitoutgrid)
    65 
    66 
    67   ! 1.  Read horizontal grid specifications
    68   !****************************************
    69 
    70   call skplin(3,unitoutgrid)
    71   read(unitoutgrid,'(4x,f11.4)') outlon0
    72   call skplin(3,unitoutgrid)
    73   read(unitoutgrid,'(4x,f11.4)') outlat0
    74   call skplin(3,unitoutgrid)
    75   read(unitoutgrid,'(4x,i5)') numxgrid
    76   call skplin(3,unitoutgrid)
    77   read(unitoutgrid,'(4x,i5)') numygrid
    78   call skplin(3,unitoutgrid)
    79   read(unitoutgrid,'(4x,f12.5)') dxout
    80   call skplin(3,unitoutgrid)
    81   read(unitoutgrid,'(4x,f12.5)') dyout
    82 
     80  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old',form='formatted',err=999)
     81
     82  ! try namelist input
     83  read(unitoutgrid,outgrid,iostat=readerror)
     84  close(unitoutgrid)
     85
     86  if ((dxout.le.0).or.(readerror.ne.0)) then
     87
     88    readerror=1
     89
     90    open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old',err=999)
     91
     92    call skplin(5,unitoutgrid)
     93
     94    ! 1.  Read horizontal grid specifications
     95    !****************************************
     96
     97    call skplin(3,unitoutgrid)
     98    read(unitoutgrid,'(4x,f11.4)') outlon0
     99    call skplin(3,unitoutgrid)
     100    read(unitoutgrid,'(4x,f11.4)') outlat0
     101    call skplin(3,unitoutgrid)
     102    read(unitoutgrid,'(4x,i5)') numxgrid
     103    call skplin(3,unitoutgrid)
     104    read(unitoutgrid,'(4x,i5)') numygrid
     105    call skplin(3,unitoutgrid)
     106    read(unitoutgrid,'(4x,f12.5)') dxout
     107    call skplin(3,unitoutgrid)
     108    read(unitoutgrid,'(4x,f12.5)') dyout
     109
     110  endif
    83111
    84112  ! Check validity of output grid (shall be within model domain)
    … …  
    91119  if ((outlon0+eps.lt.xlon0).or.(outlat0+eps.lt.ylat0) &
    92120       .or.(xr.gt.xr1+eps).or.(yr.gt.yr1+eps)) then
    93     write(*,*) 'outlon0,outlat0:'
    94121    write(*,*) outlon0,outlat0
    95     write(*,*) 'xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout:'
    96122    write(*,*) xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout
    97123    write(*,*) ' #### FLEXPART MODEL ERROR! PART OF OUTPUT    ####'
    … …  
    104130  ! 2. Count Vertical levels of output grid
    105131  !****************************************
    106   j=0
    107 100   j=j+1
     132
     133  if (readerror.ne.0) then
     134    j=0
     135100 j=j+1
    108136    do i=1,3
    109137      read(unitoutgrid,*,end=99)
    … …  
    112140    if (outhelp.eq.0.) goto 99
    113141    goto 100
    114 99   numzgrid=j-1
    115 
    116     allocate(outheight(numzgrid) &
    117          ,stat=stat)
    118     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    119     allocate(outheighthalf(numzgrid) &
    120          ,stat=stat)
    121     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    122 
    123 
    124   rewind(unitoutgrid)
    125   call skplin(29,unitoutgrid)
     14299  numzgrid=j-1
     143  else
     144    do i=1,maxoutlev
     145      if (outheights(i).lt.0) exit
     146    end do
     147    numzgrid=i-1
     148  end if
     149
     150  allocate(outheight(numzgrid),stat=stat)
     151  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheight'
     152  allocate(outheighthalf(numzgrid),stat=stat)
     153  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheighthalf'
    126154
    127155  ! 2. Vertical levels of output grid
    128156  !**********************************
    129157
    130   j=0
    131 1000   j=j+1
    132     do i=1,3
    133       read(unitoutgrid,*,end=990)
    134     end do
    135     read(unitoutgrid,'(4x,f7.1)',end=990) outhelp
    136     if (outhelp.eq.0.) goto 99
    137     outheight(j)=outhelp
    138     goto 1000
    139 990   numzgrid=j-1
    140 
     158  if (readerror.ne.0) then
     159
     160    rewind(unitoutgrid)
     161    call skplin(29,unitoutgrid)
     162
     163    do j=1,numzgrid
     164      do i=1,3
     165        read(unitoutgrid,*)
     166      end do
     167      read(unitoutgrid,'(4x,f7.1)') outhelp
     168      outheight(j)=outhelp
     169      outheights(j)=outhelp
     170    end do
     171    close(unitoutgrid)
     172
     173  else
     174
     175    do j=1,numzgrid
     176      outheight(j)=outheights(j)
     177    end do
     178
     179  endif
     180
     181  ! write outgrid file in namelist format to output directory if requested
     182  if (nmlout.eqv..true.) then
     183    ! reallocate outheights with actually required dimension for namelist writing
     184    deallocate(outheights)
     185    allocate(outheights(numzgrid),stat=stat)
     186    if (stat.ne.0) write(*,*)'ERROR: could not allocate outheights'
     187
     188    do j=1,numzgrid
     189      outheights(j)=outheight(j)
     190    end do
     191
     192    open(unitoutgrid,file=path(2)(1:length(2))//'OUTGRID.namelist',err=1000)
     193    write(unitoutgrid,nml=outgrid)
     194    close(unitoutgrid)
     195  endif
    141196
    142197  ! Check whether vertical levels are specified in ascending order
    … …  
    160215  end do
    161216
    162 
    163217  xoutshift=xlon0-outlon0
    164218  youtshift=ylat0-outlat0
    165   close(unitoutgrid)
    166 
    167     allocate(oroout(0:numxgrid-1,0:numygrid-1) &
    168          ,stat=stat)
    169     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    170     allocate(area(0:numxgrid-1,0:numygrid-1) &
    171          ,stat=stat)
    172     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    173     allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid) &
    174          ,stat=stat)
    175     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    176     allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid) &
    177          ,stat=stat)
    178     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    179     allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid) &
    180          ,stat=stat)
    181     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
     219
     220  allocate(oroout(0:numxgrid-1,0:numygrid-1),stat=stat)
     221  if (stat.ne.0) write(*,*)'ERROR: could not allocate oroout'
     222  allocate(area(0:numxgrid-1,0:numygrid-1),stat=stat)
     223  if (stat.ne.0) write(*,*)'ERROR: could not allocate area'
     224  allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
     225  if (stat.ne.0) write(*,*)'ERROR: could not allocate volume'
     226  allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
     227  if (stat.ne.0) write(*,*)'ERROR: could not allocate areaeast'
     228  allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
     229  if (stat.ne.0) write(*,*)'ERROR: could not allocate areanorth'
    182230  return
    183231
    184232
    185 999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
     233999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
    186234  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    187   write(*,*) ' #### xxx/flexpart/options                    #### '
     235  write(*,'(a)') path(1)(1:length(1))
    188236  stop
    189237
     2381000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
     239  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     240  write(*,'(a)') path(2)(1:length(2))
     241  stop
     242
    190243end subroutine readoutgrid

  • /trunk/src/readoutgrid_nest.f90

    r20 r30  
    5353  real,parameter :: eps=1.e-4
    5454
     55  integer :: readerror
    5556
     57  ! declare namelist
     58  namelist /outgridn/ &
     59    outlon0n,outlat0n, &
     60    numxgridn,numygridn, &
     61    dxoutn,dyoutn
     62
     63  ! helps identifying failed namelist input
     64  dxoutn=-1.0
    5665
    5766  ! Open the OUTGRID file and read output grid specifications
    5867  !**********************************************************
    5968
    60   open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST', &
    61        status='old',err=999)
     69  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',form='formatted',status='old',err=999)
    6270
     71  ! try namelist input
     72  read(unitoutgrid,outgridn,iostat=readerror)
     73  close(unitoutgrid)
    6374
    64   call skplin(5,unitoutgrid)
     75  if ((dxoutn.le.0).or.(readerror.ne.0)) then
    6576
     77    open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',status='old',err=999)
     78    call skplin(5,unitoutgrid)
    6679
    67   ! 1.  Read horizontal grid specifications
    68   !****************************************
     80    ! 1.  Read horizontal grid specifications
     81    !****************************************
    6982
    70   call skplin(3,unitoutgrid)
    71   read(unitoutgrid,'(4x,f11.4)') outlon0n
    72   call skplin(3,unitoutgrid)
    73   read(unitoutgrid,'(4x,f11.4)') outlat0n
    74   call skplin(3,unitoutgrid)
    75   read(unitoutgrid,'(4x,i5)') numxgridn
    76   call skplin(3,unitoutgrid)
    77   read(unitoutgrid,'(4x,i5)') numygridn
    78   call skplin(3,unitoutgrid)
    79   read(unitoutgrid,'(4x,f12.5)') dxoutn
    80   call skplin(3,unitoutgrid)
    81   read(unitoutgrid,'(4x,f12.5)') dyoutn
     83    call skplin(3,unitoutgrid)
     84    read(unitoutgrid,'(4x,f11.4)') outlon0n
     85    call skplin(3,unitoutgrid)
     86    read(unitoutgrid,'(4x,f11.4)') outlat0n
     87    call skplin(3,unitoutgrid)
     88    read(unitoutgrid,'(4x,i5)') numxgridn
     89    call skplin(3,unitoutgrid)
     90    read(unitoutgrid,'(4x,i5)') numygridn
     91    call skplin(3,unitoutgrid)
     92    read(unitoutgrid,'(4x,f12.5)') dxoutn
     93    call skplin(3,unitoutgrid)
     94    read(unitoutgrid,'(4x,f12.5)') dyoutn
    8295
     96    close(unitoutgrid)
     97  endif
    8398
    84     allocate(orooutn(0:numxgridn-1,0:numygridn-1) &
    85          ,stat=stat)
    86     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    87     allocate(arean(0:numxgridn-1,0:numygridn-1) &
    88          ,stat=stat)
    89     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    90     allocate(volumen(0:numxgridn-1,0:numygridn-1,numzgrid) &
    91          ,stat=stat)
    92     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
     99  ! write outgrid_nest file in namelist format to output directory if requested
     100  if (nmlout.eqv..true.) then
     101    open(unitoutgrid,file=path(2)(1:length(2))//'OUTGRID_NEST.namelist',err=1000)
     102    write(unitoutgrid,nml=outgridn)
     103    close(unitoutgrid)
     104  endif
     105
     106  allocate(orooutn(0:numxgridn-1,0:numygridn-1),stat=stat)
     107  if (stat.ne.0) write(*,*)'ERROR: could not allocate orooutn'
     108  allocate(arean(0:numxgridn-1,0:numygridn-1),stat=stat)
     109  if (stat.ne.0) write(*,*)'ERROR: could not allocate arean'
     110  allocate(volumen(0:numxgridn-1,0:numygridn-1,numzgrid),stat=stat)
     111  if (stat.ne.0) write(*,*)'ERROR: could not allocate volumen'
    93112
    94113  ! Check validity of output grid (shall be within model domain)
    … …  
    110129  xoutshiftn=xlon0-outlon0n
    111130  youtshiftn=ylat0-outlat0n
    112   close(unitoutgrid)
    113131  return
    114132
     133999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
     134  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     135  write(*,'(a)') path(1)(1:length(1))
     136  stop
    115137
    116 999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE OUTGRID_NEST #### '
     1381000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
    117139  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    118   write(*,*) ' #### xxx/flexpart/options                    #### '
     140  write(*,'(a)') path(2)(1:length(2))
    119141  stop
    120142

  • /trunk/src/readpaths.f90

    r20 r30  
    8888    length(numpath+2*(i-1)+2)=index(path(numpath+2*(i-1)+2),' ')-1
    8989  end do
    90   print*,length(5),length(6)
    9190
    9291

  • /trunk/src/readreceptors.f90

    r20 r30  
    2727  !                                                                            *
    2828  !     Author: A. Stohl                                                       *
    29   !                                                                            *
    3029  !     1 August 1996                                                          *
     30  !     HSO, 14 August 2013
     31  !     Added optional namelist input
    3132  !                                                                            *
    3233  !*****************************************************************************
    … …  
    5152  character(len=16) :: receptor
    5253
     54  integer :: readerror
     55  real :: lon,lat   ! for namelist input, lon/lat are used instead of x,y
     56  integer,parameter :: unitreceptorout=2
     57
     58  ! declare namelist
     59  namelist /receptors/ &
     60    receptor, lon, lat
     61
     62  lon=-999.9
    5363
    5464  ! For backward runs, do not allow receptor output. Thus, set number of receptors to zero
    … …  
    6474  !************************************************************
    6575
    66   open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS', &
    67        status='old',err=999)
     76  open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS',form='formatted',status='old',err=999)
    6877
    69   call skplin(5,unitreceptor)
     78  ! try namelist input
     79  read(unitreceptor,receptors,iostat=readerror)
    7080
     81  ! prepare namelist output if requested
     82  if (nmlout.eqv..true.) then
     83    open(unitreceptorout,file=path(2)(1:length(2))//'RECEPTORS.namelist',access='append',status='new',err=1000)
     84  endif
    7185
    72   ! Read the names and coordinates of the receptors
    73   !************************************************
     86  if ((lon.lt.-900).or.(readerror.ne.0)) then
    7487
    75   j=0
    76 100   j=j+1
     88    close(unitreceptor)
     89    open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS',status='old',err=999)
     90    call skplin(5,unitreceptor)
     91
     92    ! Read the names and coordinates of the receptors
     93    !************************************************
     94
     95    j=0
     96100 j=j+1
    7797    read(unitreceptor,*,end=99)
    7898    read(unitreceptor,*,end=99)
    … …  
    89109    endif
    90110    if (j.gt.maxreceptor) then
    91     write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
    92     write(*,*) ' #### POINTS ARE GIVEN.                       #### '
    93     write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
    94     write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
     111      write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
     112      write(*,*) ' #### POINTS ARE GIVEN.                       #### '
     113      write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
     114      write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
    95115    endif
    96116    receptorname(j)=receptor
    … …  
    100120    ym=r_earth*dy/180.*pi
    101121    receptorarea(j)=xm*ym
     122
     123    ! write receptors file in namelist format to output directory if requested
     124    if (nmlout.eqv..true.) then
     125      lon=x
     126      lat=y
     127      write(unitreceptorout,nml=receptors)
     128    endif
     129
    102130    goto 100
    103131
    104 99   numreceptor=j-1
    105   close(unitreceptor)
     13299  numreceptor=j-1
     133
     134  else ! continue with namelist input
     135
     136    j=0
     137    do while (readerror.eq.0)
     138      j=j+1
     139      lon=-999.9
     140      read(unitreceptor,receptors,iostat=readerror)
     141      if ((lon.lt.-900).or.(readerror.ne.0)) then
     142        readerror=1
     143      else
     144        if (j.gt.maxreceptor) then
     145          write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
     146          write(*,*) ' #### POINTS ARE GIVEN.                       #### '
     147          write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
     148          write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
     149        endif
     150        receptorname(j)=receptor
     151        xreceptor(j)=(lon-xlon0)/dx       ! transform to grid coordinates
     152        yreceptor(j)=(lat-ylat0)/dy
     153        xm=r_earth*cos(lat*pi/180.)*dx/180.*pi
     154        ym=r_earth*dy/180.*pi
     155        receptorarea(j)=xm*ym
     156      endif
     157
     158      ! write receptors file in namelist format to output directory if requested
     159      if (nmlout.eqv..true.) then
     160        write(unitreceptorout,nml=receptors)
     161      endif
     162
     163    end do
     164    numreceptor=j-1
     165    close(unitreceptor)
     166
     167  endif
     168
     169  if (nmlout.eqv..true.) then
     170    close(unitreceptorout)
     171  endif
     172
    106173  return
    107174
    108175
    109 999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
     176999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
    110177  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    111178  write(*,'(a)') path(1)(1:length(1))
    112179  stop
    113180
     1811000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"    #### '
     182  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     183  write(*,'(a)') path(2)(1:length(2))
     184  stop
     185
    114186end subroutine readreceptors

  • /trunk/src/readreleases.f90

    r20 r30  
    3333  !     Update: 29 January 2001                                                *
    3434  !     Release altitude can be either in magl or masl                         *
     35  !     HSO, 12 August 2013
     36  !     Added optional namelist input
    3537  !                                                                            *
    3638  !*****************************************************************************
    … …  
    5557  ! xpoint2,ypoint2     geograf. coordinates of upper right corner of release  *
    5658  !                     area                                                   *
    57   ! weta, wetb          parameters to determine the wet scavenging coefficient *
     59  ! weta, wetb          parameters to determine the below-cloud scavenging     *
     60  ! weta_in, wetb_in    parameters to determine the in-cloud scavenging        *
     61  ! wetc_in, wetd_in    parameters to determine the in-cloud scavenging        *
    5862  ! zpoint1,zpoint2     height range, over which release takes place           *
    5963  ! num_min_discrete    if less, release cannot be randomized and happens at   *
    … …  
    6973  implicit none
    7074
    71   integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
     75  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat
    7276  integer,parameter :: num_min_discrete=100
    7377  real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
    … …  
    7680  logical :: old
    7781
     82  ! help variables for namelist reading
     83  integer :: numpoints, parts, readerror
     84  integer*2 :: zkind
     85  integer :: idate1, itime1, idate2, itime2
     86  real :: lon1,lon2,lat1,lat2,z1,z2
     87  character*40 :: comment
     88  integer,parameter :: unitreleasesout=2
     89  real,allocatable, dimension (:) :: mass
     90  integer,allocatable, dimension (:) :: specnum_rel,specnum_rel2
     91
     92  ! declare namelists
     93  namelist /releases_ctrl/ &
     94    nspec, &
     95    specnum_rel
     96
     97  namelist /release/ &
     98    idate1, itime1, &
     99    idate2, itime2, &
     100    lon1, lon2, &
     101    lat1, lat2, &
     102    z1, z2, &
     103    zkind, &
     104    mass, &
     105    parts, &
     106    comment
     107
     108  numpoint=0
     109
     110  ! allocate with maxspec for first input loop
     111  allocate(mass(maxspec),stat=stat)
     112  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
     113  allocate(specnum_rel(maxspec),stat=stat)
     114  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
     115
     116  ! presetting namelist releases_ctrl
     117  nspec = -1  ! use negative value to determine failed namelist input
     118  specnum_rel = 0
     119
    78120  !sec, read release to find how many releasepoints should be allocated
    79 
    80   open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old', &
    81        err=999)
    82 
    83   ! Check the format of the RELEASES file (either in free format,
    84   ! or using a formatted mask)
    85   ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
    86   !**************************************************************************
    87 
    88   call skplin(12,unitreleases)
    89   read (unitreleases,901) line
    90 901   format (a)
    91   if (index(line,'Total') .eq. 0) then
    92     old = .false.
    93   else
    94     old = .true.
    95   endif
     121  open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',form='formatted',err=999)
     122
     123  ! check if namelist input provided
     124  read(unitreleases,releases_ctrl,iostat=readerror)
     125
     126  ! prepare namelist output if requested
     127  if (nmlout.eqv..true.) then
     128    open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',access='append',status='new',err=1000)
     129  endif
     130
     131  if ((readerror.ne.0).or.(nspec.lt.0)) then
     132
     133    ! no namelist format, close file and allow reopening in old format
     134    close(unitreleases)
     135    open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',err=999)
     136
     137    readerror=1 ! indicates old format
     138
     139    ! Check the format of the RELEASES file (either in free format,
     140    ! or using a formatted mask)
     141    ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
     142    !**************************************************************************
     143
     144    call skplin(12,unitreleases)
     145    read (unitreleases,901) line
     146901 format (a)
     147    if (index(line,'Total') .eq. 0) then
     148      old = .false.
     149    else
     150      old = .true.
     151    endif
     152    rewind(unitreleases)
     153
     154
     155    ! Skip first 11 lines (file header)
     156    !**********************************
     157
     158    call skplin(11,unitreleases)
     159
     160    read(unitreleases,*,err=998) nspec
     161    if (old) call skplin(2,unitreleases)
     162    do i=1,nspec
     163      read(unitreleases,*,err=998) specnum_rel(i)
     164      if (old) call skplin(2,unitreleases)
     165    end do
     166
     167    numpoint=0
     168100 numpoint=numpoint+1
     169    read(unitreleases,*,end=25)
     170    read(unitreleases,*,err=998,end=25) idum,idum
     171    if (old) call skplin(2,unitreleases)
     172    read(unitreleases,*,err=998) idum,idum
     173    if (old) call skplin(2,unitreleases)
     174    read(unitreleases,*,err=998) xdum
     175    if (old) call skplin(2,unitreleases)
     176    read(unitreleases,*,err=998) xdum
     177    if (old) call skplin(2,unitreleases)
     178    read(unitreleases,*,err=998) xdum
     179    if (old) call skplin(2,unitreleases)
     180    read(unitreleases,*,err=998) xdum
     181    if (old) call skplin(2,unitreleases)
     182    read(unitreleases,*,err=998) idum
     183    if (old) call skplin(2,unitreleases)
     184    read(unitreleases,*,err=998) xdum
     185    if (old) call skplin(2,unitreleases)
     186    read(unitreleases,*,err=998) xdum
     187    if (old) call skplin(2,unitreleases)
     188    read(unitreleases,*,err=998) idum
     189    if (old) call skplin(2,unitreleases)
     190    do i=1,nspec
     191      read(unitreleases,*,err=998) xdum
     192      if (old) call skplin(2,unitreleases)
     193    end do
     194    !save compoint only for the first 1000 release points
     195    read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
     196    if (old) call skplin(1,unitreleases)
     197
     198    goto 100
     199
     20025  numpoint=numpoint-1
     201
     202  else
     203
     204    readerror=0
     205    do while (readerror.eq.0)
     206      idate1=-1
     207      read(unitreleases,release,iostat=readerror)
     208      if ((idate1.lt.0).or.(readerror.ne.0)) then
     209        readerror=1
     210      else
     211        numpoint=numpoint+1
     212      endif
     213    end do
     214    readerror=0
     215  endif ! if namelist input
     216
    96217  rewind(unitreleases)
    97218
    98 
    99   ! Skip first 11 lines (file header)
    100   !**********************************
    101 
    102   call skplin(11,unitreleases)
    103 
    104 
    105   read(unitreleases,*,err=998) nspec
    106   if (old) call skplin(2,unitreleases)
    107   do i=1,nspec
    108     read(unitreleases,*,err=998) specnum_rel
    109     if (old) call skplin(2,unitreleases)
     219  ! allocate arrays of matching size for number of species (namelist output)
     220  deallocate(mass)
     221  allocate(mass(nspec),stat=stat)
     222  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
     223  allocate(specnum_rel2(nspec),stat=stat)
     224  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel2'
     225  specnum_rel2=specnum_rel(1:nspec)
     226  deallocate(specnum_rel)
     227  allocate(specnum_rel(nspec),stat=stat)
     228  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
     229  specnum_rel=specnum_rel2
     230  deallocate(specnum_rel2)
     231
     232  !allocate memory for numpoint releaspoints
     233  allocate(ireleasestart(numpoint),stat=stat)
     234  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleasestart'
     235  allocate(ireleaseend(numpoint),stat=stat)
     236  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleaseend'
     237  allocate(xpoint1(numpoint),stat=stat)
     238  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint1'
     239  allocate(xpoint2(numpoint),stat=stat)
     240  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint2'
     241  allocate(ypoint1(numpoint),stat=stat)
     242  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint1'
     243  allocate(ypoint2(numpoint),stat=stat)
     244  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint2'
     245  allocate(zpoint1(numpoint),stat=stat)
     246  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint1'
     247  allocate(zpoint2(numpoint),stat=stat)
     248  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint2'
     249  allocate(kindz(numpoint),stat=stat)
     250  if (stat.ne.0) write(*,*)'ERROR: could not allocate kindz'
     251  allocate(xmass(numpoint,maxspec),stat=stat)
     252  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmass'
     253  allocate(rho_rel(numpoint),stat=stat)
     254  if (stat.ne.0) write(*,*)'ERROR: could not allocate rho_rel'
     255  allocate(npart(numpoint),stat=stat)
     256  if (stat.ne.0) write(*,*)'ERROR: could not allocate npart'
     257  allocate(xmasssave(numpoint),stat=stat)
     258  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmasssave'
     259
     260  write (*,*) 'Releasepoints : ', numpoint
     261
     262  do i=1,numpoint
     263    xmasssave(i)=0.
    110264  end do
    111 
    112   numpoint=0
    113 100   numpoint=numpoint+1
    114   read(unitreleases,*,end=25)
    115   read(unitreleases,*,err=998,end=25) idum,idum
    116   if (old) call skplin(2,unitreleases)
    117   read(unitreleases,*,err=998) idum,idum
    118   if (old) call skplin(2,unitreleases)
    119   read(unitreleases,*,err=998) xdum
    120   if (old) call skplin(2,unitreleases)
    121   read(unitreleases,*,err=998) xdum
    122   if (old) call skplin(2,unitreleases)
    123   read(unitreleases,*,err=998) xdum
    124   if (old) call skplin(2,unitreleases)
    125   read(unitreleases,*,err=998) xdum
    126   if (old) call skplin(2,unitreleases)
    127   read(unitreleases,*,err=998) idum
    128   if (old) call skplin(2,unitreleases)
    129   read(unitreleases,*,err=998) xdum
    130   if (old) call skplin(2,unitreleases)
    131   read(unitreleases,*,err=998) xdum
    132   if (old) call skplin(2,unitreleases)
    133   read(unitreleases,*,err=998) idum
    134   if (old) call skplin(2,unitreleases)
    135   do i=1,nspec
    136     read(unitreleases,*,err=998) xdum
    137     if (old) call skplin(2,unitreleases)
    138   end do
    139   !save compoint only for the first 1000 release points
    140   read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
    141   if (old) call skplin(1,unitreleases)
    142 
    143   goto 100
    144 
    145 25   numpoint=numpoint-1
    146 
    147   !allocate memory for numpoint releaspoint
    148     allocate(ireleasestart(numpoint) &
    149          ,stat=stat)
    150     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    151     allocate(ireleaseend(numpoint) &
    152          ,stat=stat)
    153     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    154     allocate(xpoint1(numpoint) &
    155          ,stat=stat)
    156     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    157     allocate(xpoint2(numpoint) &
    158          ,stat=stat)
    159     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    160     allocate(ypoint1(numpoint) &
    161          ,stat=stat)
    162     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    163     allocate(ypoint2(numpoint) &
    164          ,stat=stat)
    165     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    166     allocate(zpoint1(numpoint) &
    167          ,stat=stat)
    168     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    169     allocate(zpoint2(numpoint) &
    170          ,stat=stat)
    171     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    172     allocate(kindz(numpoint) &
    173          ,stat=stat)
    174     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    175     allocate(xmass(numpoint,maxspec) &
    176          ,stat=stat)
    177     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    178     allocate(rho_rel(numpoint) &
    179          ,stat=stat)
    180     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    181     allocate(npart(numpoint) &
    182          ,stat=stat)
    183     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    184     allocate(xmasssave(numpoint) &
    185          ,stat=stat)
    186     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    187 
    188    write (*,*) ' Releasepoints allocated: ', numpoint
    189 
    190    do i=1,numpoint
    191      xmasssave(i)=0.
    192    end do
    193265
    194266  !now save the information
    … …  
    201273  end do
    202274
    203   rewind(unitreleases)
    204 
    205 
    206   ! Skip first 11 lines (file header)
    207   !**********************************
    208 
    209   call skplin(11,unitreleases)
    210 
    211 
    212   ! Assign species-specific parameters needed for physical processes
    213   !*************************************************************************
    214 
    215   read(unitreleases,*,err=998) nspec
    216   if (nspec.gt.maxspec) goto 994
    217   if (old) call skplin(2,unitreleases)
     275  if (readerror.ne.0) then
     276    ! Skip first 11 lines (file header)
     277    !**********************************
     278
     279    call skplin(11,unitreleases)
     280
     281    ! Assign species-specific parameters needed for physical processes
     282    !*************************************************************************
     283
     284    read(unitreleases,*,err=998) nspec
     285    if (nspec.gt.maxspec) goto 994
     286    if (old) call skplin(2,unitreleases)
     287  endif
     288
     289  ! namelist output
     290  if (nmlout.eqv..true.) then
     291    write(unitreleasesout,nml=releases_ctrl)
     292  endif
     293
    218294  do i=1,nspec
    219     read(unitreleases,*,err=998) specnum_rel
    220     if (old) call skplin(2,unitreleases)
    221     call readspecies(specnum_rel,i)
    222 
    223   ! For backward runs, only 1 species is allowed
    224   !*********************************************
    225 
    226   !if ((ldirect.lt.0).and.(nspec.gt.1)) then
    227   !write(*,*) '#####################################################'
    228   !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
    229   !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
    230   !write(*,*) '#####################################################'
    231   !  stop
    232   !endif
    233 
    234   ! Molecular weight
    235   !*****************
    236 
    237     if (((iout.eq.2).or.(iout.eq.3)).and. &
    238          (weightmolar(i).lt.0.)) then
     295    if (readerror.ne.0) then
     296      read(unitreleases,*,err=998) specnum_rel(i)
     297      if (old) call skplin(2,unitreleases)
     298      call readspecies(specnum_rel(i),i)
     299    else
     300      call readspecies(specnum_rel(i),i)
     301    endif
     302
     303    ! For backward runs, only 1 species is allowed
     304    !*********************************************
     305
     306    !if ((ldirect.lt.0).and.(nspec.gt.1)) then
     307    !write(*,*) '#####################################################'
     308    !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
     309    !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
     310    !write(*,*) '#####################################################'
     311    !  stop
     312    !endif
     313
     314    ! Molecular weight
     315    !*****************
     316
     317    if (((iout.eq.2).or.(iout.eq.3)).and.(weightmolar(i).lt.0.)) then
    239318      write(*,*) 'For mixing ratio output, valid molar weight'
    240319      write(*,*) 'must be specified for all simulated species.'
    … …  
    244323    endif
    245324
    246 
    247   ! Radioactive decay
    248   !******************
     325    ! Radioactive decay
     326    !******************
    249327
    250328    decay(i)=0.693147/decay(i) !conversion half life to decay constant
    251329
    252330
    253   ! Dry deposition of gases
    254   !************************
    255 
    256     if (reldiff(i).gt.0.) &
    257          rm(i)=1./(henry(i)/3000.+100.*f0(i))    ! mesophyll resistance
    258 
    259   ! Dry deposition of particles
    260   !****************************
     331    ! Dry deposition of gases
     332    !************************
     333
     334    if (reldiff(i).gt.0.) rm(i)=1./(henry(i)/3000.+100.*f0(i))    ! mesophyll resistance
     335
     336    ! Dry deposition of particles
     337    !****************************
    261338
    262339    vsetaver(i)=0.
    263340    cunningham(i)=0.
    264341    dquer(i)=dquer(i)*1000000.         ! Conversion m to um
    265     if (density(i).gt.0.) then                  ! Additional parameters
    266      call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
     342    if (density(i).gt.0.) then         ! Additional parameters
     343      call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
    267344      do j=1,ni
    268345        fract(i,j)=fracth(j)
    … …  
    272349        vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
    273350      end do
    274       write(*,*) 'Average setting velocity: ',i,vsetaver(i)
    275     endif
    276 
    277   ! Dry deposition for constant deposition velocity
    278   !************************************************
     351      write(*,*) 'Average settling velocity: ',i,vsetaver(i)
     352    endif
     353
     354    ! Dry deposition for constant deposition velocity
     355    !************************************************
    279356
    280357    dryvel(i)=dryvel(i)*0.01         ! conversion to m/s
    281358
    282   ! Check if wet deposition or OH reaction shall be calculated
    283   !***********************************************************
     359    ! Check if wet deposition or OH reaction shall be calculated
     360    !***********************************************************
    284361    if (weta(i).gt.0.)  then
    285362      WETDEP=.true.
    286       write (*,*) 'Wetdeposition switched on: ',weta(i),i
    287     endif
     363      write (*,*) 'Below-cloud scavenging: ON'
     364      write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
     365    else
     366      write (*,*) 'Below-cloud scavenging: OFF'
     367    endif
     368   
     369    ! NIK 31.01.2013 + 10.12.2013
     370    if (weta_in(i).gt.0.)  then
     371      WETDEP=.true.
     372      write (*,*) 'In-cloud scavenging: ON'
     373      write (*,*) 'In-cloud scavenging coefficients: ',weta_in(i),wetb_in(i), wetc_in(i), wetd_in(i),i
     374    else
     375      write (*,*) 'In-cloud scavenging: OFF'
     376    endif
     377
    288378    if (ohreact(i).gt.0) then
    289379      OHREA=.true.
    290       write (*,*) 'OHreaction switched on: ',ohreact(i),i
    291     endif
    292 
    293 
    294     if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or. &
    295          (dryvel(i).gt.0.)) then
     380      write (*,*) 'OHreaction: ON (',ohreact(i),i,')'
     381    endif
     382
     383    if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or.(dryvel(i).gt.0.)) then
    296384      DRYDEP=.true.
    297385      DRYDEPSPEC(i)=.true.
    … …  
    300388  end do
    301389
    302     if (WETDEP.or.DRYDEP) DEP=.true.
     390  if (WETDEP.or.DRYDEP) DEP=.true.
    303391
    304392  ! Read specifications for each release point
    305393  !*******************************************
    306 
     394  numpoints=numpoint
    307395  numpoint=0
    308396  numpartmax=0
    309397  releaserate=0.
    310 1000   numpoint=numpoint+1
    311   read(unitreleases,*,end=250)
    312   read(unitreleases,*,err=998,end=250) id1,it1
    313   if (old) call skplin(2,unitreleases)
    314   read(unitreleases,*,err=998) id2,it2
    315   if (old) call skplin(2,unitreleases)
    316   read(unitreleases,*,err=998) xpoint1(numpoint)
    317   if (old) call skplin(2,unitreleases)
    318   read(unitreleases,*,err=998) ypoint1(numpoint)
    319   if (old) call skplin(2,unitreleases)
    320   read(unitreleases,*,err=998) xpoint2(numpoint)
    321   if (old) call skplin(2,unitreleases)
    322   read(unitreleases,*,err=998) ypoint2(numpoint)
    323   if (old) call skplin(2,unitreleases)
    324   read(unitreleases,*,err=998) kindz(numpoint)
    325   if (old) call skplin(2,unitreleases)
    326   read(unitreleases,*,err=998) zpoint1(numpoint)
    327   if (old) call skplin(2,unitreleases)
    328   read(unitreleases,*,err=998) zpoint2(numpoint)
    329   if (old) call skplin(2,unitreleases)
    330   read(unitreleases,*,err=998) npart(numpoint)
    331   if (old) call skplin(2,unitreleases)
    332   do i=1,nspec
    333     read(unitreleases,*,err=998) xmass(numpoint,i)
    334     if (old) call skplin(2,unitreleases)
    335   end do
    336   !save compoint only for the first 1000 release points
    337   if (numpoint.le.1000) then
    338     read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
     398101   numpoint=numpoint+1
     399
     400  if (readerror.lt.1) then ! reading namelist format
     401
     402    if (numpoint.gt.numpoints) goto 250
     403    zkind = 1
     404    mass = 0
     405    parts = 0
     406    comment = ' '
     407    read(unitreleases,release,iostat=readerror)
     408    id1=idate1
     409    it1=itime1
     410    id2=idate2
     411    it2=itime2
     412    xpoint1(numpoint)=lon1
     413    xpoint2(numpoint)=lon2
     414    ypoint1(numpoint)=lat1
     415    ypoint2(numpoint)=lat2
     416    zpoint1(numpoint)=z1
     417    zpoint2(numpoint)=z2
     418    kindz(numpoint)=zkind
     419    do i=1,nspec
     420      xmass(numpoint,i)=mass(i)
     421    end do
     422    npart(numpoint)=parts
     423    compoint(min(1001,numpoint))=comment
     424
     425    ! namelist output
     426    if (nmlout.eqv..true.) then
     427      write(unitreleasesout,nml=release)
     428    endif
     429
    339430  else
    340     read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
    341   endif
    342   if (old) call skplin(1,unitreleases)
    343   if (numpoint.le.1000) then
    344     if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
    345          (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
    346          (compoint(numpoint)(1:8).eq.'        ')) goto 250
    347   else
    348     if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
    349          (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
    350   endif
    351 
     431
     432    read(unitreleases,*,end=250)
     433    read(unitreleases,*,err=998,end=250) id1,it1
     434    if (old) call skplin(2,unitreleases)
     435    read(unitreleases,*,err=998) id2,it2
     436    if (old) call skplin(2,unitreleases)
     437    read(unitreleases,*,err=998) xpoint1(numpoint)
     438    if (old) call skplin(2,unitreleases)
     439    read(unitreleases,*,err=998) ypoint1(numpoint)
     440    if (old) call skplin(2,unitreleases)
     441    read(unitreleases,*,err=998) xpoint2(numpoint)
     442    if (old) call skplin(2,unitreleases)
     443    read(unitreleases,*,err=998) ypoint2(numpoint)
     444    if (old) call skplin(2,unitreleases)
     445    read(unitreleases,*,err=998) kindz(numpoint)
     446    if (old) call skplin(2,unitreleases)
     447    read(unitreleases,*,err=998) zpoint1(numpoint)
     448    if (old) call skplin(2,unitreleases)
     449    read(unitreleases,*,err=998) zpoint2(numpoint)
     450    if (old) call skplin(2,unitreleases)
     451    read(unitreleases,*,err=998) npart(numpoint)
     452    if (old) call skplin(2,unitreleases)
     453    do i=1,nspec
     454      read(unitreleases,*,err=998) xmass(numpoint,i)
     455      if (old) call skplin(2,unitreleases)
     456      mass(i)=xmass(numpoint,i)
     457    end do
     458    !save compoint only for the first 1000 release points
     459    if (numpoint.le.1000) then
     460      read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
     461      comment=compoint(numpoint)(1:40)
     462    else
     463      read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
     464      comment=compoint(1001)(1:40)
     465    endif
     466    if (old) call skplin(1,unitreleases)
     467
     468    ! namelist output
     469    if (nmlout.eqv..true.) then
     470      idate1=id1
     471      itime1=it1
     472      idate2=id2
     473      itime2=it2
     474      lon1=xpoint1(numpoint)
     475      lon2=xpoint2(numpoint)
     476      lat1=ypoint1(numpoint)
     477      lat2=ypoint2(numpoint)
     478      z1=zpoint1(numpoint)
     479      z2=zpoint2(numpoint)
     480      zkind=kindz(numpoint)
     481      parts=npart(numpoint)
     482      write(unitreleasesout,nml=release)
     483    endif
     484
     485    if (numpoint.le.1000) then
     486      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
     487           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
     488           (compoint(numpoint)(1:8).eq.'        ')) goto 250
     489    else
     490      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
     491           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
     492    endif
     493
     494  endif ! if namelist format
    352495
    353496  ! If a release point contains no particles, stop and issue error message
    … …  
    420563  endif
    421564
    422 
    423   ! Check, whether the total number of particles may exceed totally allowed
    424   ! number of particles at some time during the simulation
    425   !************************************************************************
    426 
    427565  ! Determine the release rate (particles per second) and total number
    428566  ! of particles released during the simulation
    … …  
    436574  endif
    437575  numpartmax=numpartmax+npart(numpoint)
    438   goto 1000
    439 
     576  goto 101
    440577
    441578250   close(unitreleases)
    442579
    443   write (*,*) ' Particles allocated for this run: ',maxpart, ', released in simulation: ',  numpartmax
     580  if (nmlout.eqv..true.) then
     581    close(unitreleasesout)
     582  endif
     583
     584  write (*,*) 'Particles allocated (maxpart)  : ',maxpart
     585  write (*,*) 'Particles released (numpartmax): ',numpartmax
    444586  numpoint=numpoint-1
    445587
    … …  
    449591    maxpointspec_act=1
    450592  endif
     593
     594  ! Check, whether the total number of particles may exceed totally allowed
     595  ! number of particles at some time during the simulation
     596  !************************************************************************
    451597
    452598  if (releaserate.gt. &
    … …  
    506652  stop
    507653
     6541000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RELEASES"    #### '
     655  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     656  write(*,'(a)') path(2)(1:length(2))
     657  stop
     658
    508659end subroutine readreleases

  • /trunk/src/readspecies.f90

    r20 r30  
    3030  !                                                                            *
    3131  !   11 July 1996                                                             *
     32  !
     33  !   Changes:                                                                 *
     34  !   N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging *
     35  !                                                                            *
     36  !   HSO, 13 August 2013
     37  !   added optional namelist input
    3238  !                                                                            *
    3339  !*****************************************************************************
    … …  
    3642  ! decaytime(maxtable)  half time for radiological decay                      *
    3743  ! specname(maxtable)   names of chemical species, radionuclides              *
    38   ! wetscava, wetscavb   Parameters for determining scavenging coefficient     *
     44  ! weta, wetb           Parameters for determining below-cloud scavenging     *
     45  ! weta_in              Parameter for determining in-cloud scavenging         *
     46  ! wetb_in              Parameter for determining in-cloud scavenging         *
     47  ! wetc_in              Parameter for determining in-cloud scavenging         *
     48  ! wetd_in              Parameter for determining in-cloud scavenging         *
    3949  ! ohreact              OH reaction rate                                      *
    4050  ! id_spec              SPECIES number as referenced in RELEASE file          *
    … …  
    5565  logical :: spec_found
    5666
     67  character(len=16) :: pspecies
     68  real :: pdecay, pweta, pwetb, preldiff, phenry, pf0, pdensity, pdquer
     69  real :: pdsigma, pdryvel, pweightmolar, pohreact, pspec_ass, pkao
     70  real :: pweta_in, pwetb_in, pwetc_in, pwetd_in
     71  integer :: readerror
     72
     73  ! declare namelist
     74  namelist /species_params/ &
     75   pspecies, pdecay, pweta, pwetb, &
     76   pweta_in, pwetb_in, pwetc_in, pwetd_in, &
     77   preldiff, phenry, pf0, pdensity, pdquer, &
     78   pdsigma, pdryvel, pweightmolar, pohreact, pspec_ass, pkao
     79
     80  pspecies=" "
     81  pdecay=-999.9
     82  pweta=-9.9E-09
     83  pwetb=0.0
     84  pweta_in=-9.9E-09
     85  pwetb_in=-9.9E-09
     86  pwetc_in=-9.9E-09
     87  pwetd_in=-9.9E-09
     88  preldiff=-9.9
     89  phenry=0.0
     90  pf0=0.0
     91  pdensity=-9.9E09
     92  pdquer=0.0
     93  pdsigma=0.0
     94  pdryvel=-9.99
     95  pohreact=-9.9E-09
     96  pspec_ass=-9
     97  pkao=-99.99
     98  pweightmolar=-789.0 ! read failure indicator value
     99
    57100  ! Open the SPECIES file and read species names and properties
    58101  !************************************************************
    59102  specnum(pos_spec)=id_spec
    60103  write(aspecnumb,'(i3.3)') specnum(pos_spec)
    61   open(unitspecies,file= &
    62        path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old', &
    63        err=998)
     104  open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998)
    64105  !write(*,*) 'reading SPECIES',specnum(pos_spec)
    65106
    66107  ASSSPEC=.FALSE.
    67108
    68   do i=1,6
    69     read(unitspecies,*)
    70   end do
     109  ! try namelist input
     110  read(unitspecies,species_params,iostat=readerror)
     111  close(unitspecies)
     112   
     113  if ((pweightmolar.eq.-789.0).or.(readerror.ne.0)) then ! no namelist found
     114
     115    readerror=1
     116
     117    open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',err=998)
     118
     119    do i=1,6
     120      read(unitspecies,*)
     121    end do
    71122
    72123    read(unitspecies,'(a10)',end=22) species(pos_spec)
    … …  
    78129    read(unitspecies,'(f18.2)',end=22) wetb(pos_spec)
    79130  !  write(*,*) wetb(pos_spec)
     131
     132  !*** NIK 31.01.2013: including in-cloud scavening parameters
     133   read(unitspecies,'(e18.1)',end=22) weta_in(pos_spec)
     134  !  write(*,*) weta_in(pos_spec)
     135   read(unitspecies,'(f18.2)',end=22) wetb_in(pos_spec)
     136  !  write(*,*) wetb_in(pos_spec)
     137   read(unitspecies,'(f18.2)',end=22) wetc_in(pos_spec)
     138  !  write(*,*) wetc_in(pos_spec)
     139   read(unitspecies,'(f18.2)',end=22) wetd_in(pos_spec)
     140  !  write(*,*) wetd_in(pos_spec)
     141
    80142    read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
    81143  !  write(*,*) reldiff(pos_spec)
    … …  
    100162    read(unitspecies,'(f18.2)',end=22) kao(pos_spec)
    101163  !       write(*,*) kao(pos_spec)
    102     i=pos_spec
    103 
    104     if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
    105        if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
     164
     165    pspecies=species(pos_spec)
     166    pdecay=decay(pos_spec)
     167    pweta=weta(pos_spec)
     168    pwetb=wetb(pos_spec)
     169    pweta_in=weta_in(pos_spec)
     170    pwetb_in=wetb_in(pos_spec)
     171    pwetc_in=wetc_in(pos_spec)
     172    pwetd_in=wetd_in(pos_spec)
     173    preldiff=reldiff(pos_spec)
     174    phenry=henry(pos_spec)
     175    pf0=f0(pos_spec)
     176    pdensity=density(pos_spec)
     177    pdquer=dquer(pos_spec)
     178    pdsigma=dsigma(pos_spec)
     179    pdryvel=dryvel(pos_spec)
     180    pweightmolar=weightmolar(pos_spec)
     181    pohreact=ohreact(pos_spec)
     182    pspec_ass=spec_ass(pos_spec)
     183    pkao=kao(pos_spec)
     184
     185  else
     186
     187    species(pos_spec)=pspecies
     188    decay(pos_spec)=pdecay
     189    weta(pos_spec)=pweta
     190    wetb(pos_spec)=pwetb
     191    weta_in(pos_spec)=pweta_in
     192    wetb_in(pos_spec)=pwetb_in
     193    wetc_in(pos_spec)=pwetc_in
     194    wetd_in(pos_spec)=pwetd_in
     195    reldiff(pos_spec)=preldiff
     196    henry(pos_spec)=phenry
     197    f0(pos_spec)=pf0
     198    density(pos_spec)=pdensity
     199    dquer(pos_spec)=pdquer
     200    dsigma(pos_spec)=pdsigma
     201    dryvel(pos_spec)=pdryvel
     202    weightmolar(pos_spec)=pweightmolar
     203    ohreact(pos_spec)=pohreact
     204    spec_ass(pos_spec)=pspec_ass
     205    kao(pos_spec)=pkao
     206
     207  endif
     208
     209  i=pos_spec
     210
     211  if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
     212   if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
     213  endif
     214
     215  if (spec_ass(pos_spec).gt.0) then
     216    spec_found=.FALSE.
     217    do j=1,pos_spec-1
     218      if (spec_ass(pos_spec).eq.specnum(j)) then
     219        spec_ass(pos_spec)=j
     220        spec_found=.TRUE.
     221        ASSSPEC=.TRUE.
     222      endif
     223    end do
     224    if (spec_found.eqv..FALSE.) then
     225      goto 997
    106226    endif
    107 
    108     if (spec_ass(pos_spec).gt.0) then
    109        spec_found=.FALSE.
    110        do j=1,pos_spec-1
    111           if (spec_ass(pos_spec).eq.specnum(j)) then
    112              spec_ass(pos_spec)=j
    113              spec_found=.TRUE.
    114              ASSSPEC=.TRUE.
    115           endif
    116        end do
    117        if (spec_found.eqv..FALSE.) then
    118           goto 997
    119        endif
    120     endif
    121 
    122     if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
    123     if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
    124 
    125     if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
    126       write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
    127       write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
    128       write(*,*) '#### PARTICLE AND GAS.                       ####'
    129       write(*,*) '#### SPECIES NUMBER',aspecnumb
    130       stop
    131     endif
     227  endif
     228
     229  if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
     230  if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
     231
     232  if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
     233    write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
     234    write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
     235    write(*,*) '#### PARTICLE AND GAS.                       ####'
     236    write(*,*) '#### SPECIES NUMBER',aspecnumb
     237    stop
     238  endif
    13223920   continue
    133240
    … …  
    135242  ! Read in daily and day-of-week variation of emissions, if available
    136243  !*******************************************************************
    137 
    138     do j=1,24           ! initialize everything to no variation
    139       area_hour(i,j)=1.
    140       point_hour(i,j)=1.
    141     end do
    142     do j=1,7
    143       area_dow(i,j)=1.
    144       point_dow(i,j)=1.
    145     end do
     244  ! HSO: This is not yet implemented as namelist parameters
     245  ! default values set to 1
     246
     247  do j=1,24           ! initialize everything to no variation
     248    area_hour(i,j)=1.
     249    point_hour(i,j)=1.
     250  end do
     251  do j=1,7
     252    area_dow(i,j)=1.
     253    point_dow(i,j)=1.
     254  end do
     255
     256  if (readerror.ne.0) then ! text format input
    146257
    147258    read(unitspecies,*,end=22)
    … …  
    154265    end do
    155266
    156 22   close(unitspecies)
    157 
    158    return
     267  endif
     268
     26922 close(unitspecies)
     270
     271  ! namelist output if requested
     272  if (nmlout.eqv..true.) then
     273    open(unitspecies,file=path(2)(1:length(2))//'SPECIES_'//aspecnumb//'.namelist',access='append',status='new',err=1000)
     274    write(unitspecies,nml=species_params)
     275    close(unitspecies)
     276  endif
     277
     278  return
    159279
    160280996   write(*,*) '#####################################################'
    … …  
    185305  stop
    186306
     3071000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "SPECIES_',aspecnumb,'.namelist'
     308  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     309  write(*,'(a)') path(2)(1:length(2))
     310  stop
    187311
    188312end subroutine readspecies

  • /trunk/src/readwind.f90

    r20 r30  
    7171  integer :: iret
    7272  integer :: igrib
    73   integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
     73  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
    7474  integer :: gotGrid
    7575  !HSO  end
    … …  
    9090  integer :: isec1(56),isec2(22+nxmax+nymax)
    9191  real(kind=4) :: zsec4(jpunp)
    92   real(kind=4) :: xaux,yaux,xaux0,yaux0
     92  real(kind=4) :: xaux,yaux
    9393  real(kind=8) :: xauxin,yauxin
    9494  real,parameter :: eps=1.e-4
    9595  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
    9797
    9898  logical :: hflswitch,strswitch
    … …  
    101101  character(len=24) :: gribErrorMsg = 'Error reading grib file'
    102102  character(len=20) :: gribFunction = 'readwind'
    103 
    104   !HSO conversion of ECMWF etadot to etadot*dp/deta
    105   logical :: etacon=.false.
    106   real,parameter :: p00=101325.
    107   real :: dak,dbk
    108103
    109104  hflswitch=.false.
    … …  
    155150  endif
    156151
     152  conversion_factor=1.
     153
    157154  else
    158155
    … …  
    168165  call grib_check(iret,gribFunction,gribErrorMsg)
    169166  call grib_get_int(igrib,'level',valSurf,iret)
     167  call grib_check(iret,gribFunction,gribErrorMsg)
     168  call grib_get_int(igrib,'paramId',parId,iret)
    170169  call grib_check(iret,gribFunction,gribErrorMsg)
    171170
    … …  
    177176  isec1(8)=-1
    178177  isec1(8)=valSurf     ! level
     178  conversion_factor=1.
    179179  if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
    180180    isec1(6)=130         ! indicatorOfParameter
    … …  
    188188    isec1(6)=134         ! indicatorOfParameter
    189189  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
    190192    isec1(6)=135         ! indicatorOfParameter
    191193  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    … …  
    201203  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    202204    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
    204207    isec1(6)=164         ! indicatorOfParameter
    205   elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
     208  elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
    206209    isec1(6)=142         ! indicatorOfParameter
    207210  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    208211    isec1(6)=143         ! indicatorOfParameter
     212    conversion_factor=1000.
    209213  elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
    210214    isec1(6)=146         ! indicatorOfParameter
    211215  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    212216    isec1(6)=176         ! indicatorOfParameter
    213   elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
     217  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
    214218    isec1(6)=180         ! indicatorOfParameter
    215   elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
     219  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
    216220    isec1(6)=181         ! indicatorOfParameter
    217221  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    218222    isec1(6)=129         ! indicatorOfParameter
    219   elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
     223  elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
    220224    isec1(6)=160         ! indicatorOfParameter
    221225  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    … …  
    223227    isec1(6)=172         ! indicatorOfParameter
    224228  else
    225     print*,'***ERROR: undefined GRiB2 message found!',discipl, &
     229    print*,'***WARNING: undefined GRiB2 message found!',discipl, &
    226230         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
    227235  endif
    228236
    … …  
    237245  !HSO  get the required fields from section 2 in a gribex compatible manner
    238246  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))
    247252  call grib_check(iret,gribFunction,gribErrorMsg)
    248253  ! CHECK GRID SPECIFICATIONS
    … …  
    250255  if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
    251256  if(isec2(12)/2-1.ne.nlev_ec) &
    252        stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
     257  stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
    253258  endif ! ifield
    254259
    255260  !HSO  get the second part of the grid dimensions only from GRiB1 messages
    256   if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
     261  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
    257262    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
    258263         xauxin,iret)
    … …  
    261266         yauxin,iret)
    262267    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
    263271    xaux=xauxin+real(nxshift)*dx
    264272    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'
    275278    gotGrid=1
    276279  endif ! gotGrid
    … …  
    298301           zsec4(nxfield*(ny-j-1)+i+1)
    299302      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
    301304      if(isec1(6).eq.151) msl(i,j,1,n)= &!! SEA LEVEL PRESS.
    302305           zsec4(nxfield*(ny-j-1)+i+1)
    … …  
    316319      endif
    317320      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
    319322        if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
    320323      endif
    … …  
    366369      do j=0,nymin1
    367370        wwh(i,j,nlev_ec+1)=0.
    368       end do
    369     end do
    370   endif
    371 
    372   ! convert from ECMWF etadot to etadot*dp/deta as needed by FLEXPART
    373   if(etacon.eqv..true.) then
    374     do k=1,nwzmax
    375       dak=akm(k+1)-akm(k)
    376       dbk=bkm(k+1)-bkm(k)
    377       do i=0,nxmin1
    378         do j=0,nymin1
    379           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) then
    381             wwh(i,j,k)=wwh(i,j,k)-wwh(i,j,k-1)
    382           endif
    383         end do
    384371      end do
    385372    end do
    … …  
    479466
    480467end subroutine readwind
     468

  • /trunk/src/readwind_nests.f90

    r20 r30  
    5151  integer :: igrib
    5252  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
     53  integer :: parId !!added by mc for making it consistent with new readwind.f90
    5354  integer :: gotGrid
    5455  !HSO  end
    … …  
    6970  integer :: isec1(56),isec2(22+nxmaxn+nymaxn)
    7071  real(kind=4) :: zsec4(jpunp)
     72  real(kind=4) :: xaux,yaux
    7173  real(kind=8) :: xauxin,yauxin
    72   real(kind=4) :: xaux,yaux,xaux0,yaux0
     74  real,parameter :: eps=1.e-4
    7375  real :: ewss(0:nxmaxn-1,0:nymaxn-1),nsss(0:nxmaxn-1,0:nymaxn-1)
    7476  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
     77  real :: conversion_factor !added by mc to make it consistent with new gridchek.f90
    7578
    7679  logical :: hflswitch,strswitch
    … …  
    134137  endif
    135138
     139  conversion_factor=1.
     140
     141
    136142  else
    137143
    … …  
    148154  call grib_get_int(igrib,'level',valSurf,iret)
    149155  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
    150158
    151159  !print*,discipl,parCat,parNum,typSurf,valSurf
    … …  
    156164  isec1(8)=-1
    157165  isec1(8)=valSurf     ! level
     166   conversion_factor=1.
    158167  if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
    159168    isec1(6)=130         ! indicatorOfParameter
    … …  
    166175  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
    167176    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 !
    169178    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
    170181  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    171182    isec1(6)=151         ! indicatorOfParameter
    … …  
    180191  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    181192    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
    183195    isec1(6)=164         ! indicatorOfParameter
    184   elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
     196  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
    185197    isec1(6)=142         ! indicatorOfParameter
    186198  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    187199    isec1(6)=143         ! indicatorOfParameter
     200    conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
    188201  elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
    189202    isec1(6)=146         ! indicatorOfParameter
    190203  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    191204    isec1(6)=176         ! indicatorOfParameter
    192   elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
     205  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
    193206    isec1(6)=180         ! indicatorOfParameter
    194   elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
     207  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
    195208    isec1(6)=181         ! indicatorOfParameter
    196209  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    197210    isec1(6)=129         ! indicatorOfParameter
    198   elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
     211   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
    199212    isec1(6)=160         ! indicatorOfParameter
    200213  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    … …  
    202215    isec1(6)=172         ! indicatorOfParameter
    203216  else
    204     print*,'***ERROR: undefined GRiB2 message found!',discipl, &
     217    print*,'***WARNING: undefined GRiB2 message found!',discipl, &
    205218         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
    206223  endif
    207224
    … …  
    227244  ! CHECK GRID SPECIFICATIONS
    228245  if(isec2(2).ne.nxn(l)) stop &
    229        'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
     246  'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
    230247  if(isec2(3).ne.nyn(l)) stop &
    231        'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
     248  'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
    232249  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'
    234251  endif ! ifield
    235252
    236253  !HSO  get the second part of the grid dimensions only from GRiB1 messages
    237   if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
     254 if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consisitent with new readwind.f90
    238255    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
    239256         xauxin,iret)
    … …  
    242259         yauxin,iret)
    243260    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
    244264    xaux=xauxin
    245265    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'
    258270    gotGrid=1
    259271  endif
    … …  
    280292             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
    281293        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!
    283295        if(isec1(6).eq.151) msln(i,j,1,n,l)= &!! SEA LEVEL PRESS.
    284296             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
    … …  
    298310        endif
    299311        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
    301313          if (convprecn(i,j,1,n,l).lt.0.) convprecn(i,j,1,n,l)=0.
    302314        endif

  • /trunk/src/timemanager.f90

    r20 r30  
    4343  !     call convection BEFORE new fields are read in BWD mode                 *
    4444  !  Changes Caroline Forster, Feb 2005                                        *
    45   !new interface between flexpart and convection scheme                        *
    46   !Emanuel's latest subroutine convect43c.f is used                            *
     45  !   new interface between flexpart and convection scheme                     *
     46  !   Emanuel's latest subroutine convect43c.f is used                         *
     47  !  Changes Stefan Henne, Harald Sodemann, 2013-2014                          *
     48  !   added netcdf output code                                                 *
    4749  !*****************************************************************************
    4850  !                                                                            *
    4951  ! Variables:                                                                 *
    50   ! DEP                .true. if either wet or dry deposition is switched on   *
     52  ! dep                .true. if either wet or dry deposition is switched on   *
    5153  ! decay(maxspec) [1/s] decay constant for radioactive decay                  *
    52   ! DRYDEP             .true. if dry deposition is switched on                 *
     54  ! drydep             .true. if dry deposition is switched on                 *
    5355  ! ideltas [s]        modelling period                                        *
    5456  ! itime [s]          actual temporal position of calculation                 *
    … …  
    7072  ! prob               probability of absorption at ground due to dry          *
    7173  !                    deposition                                              *
    72   ! WETDEP             .true. if wet deposition is switched on                 *
     74  ! wetdep             .true. if wet deposition is switched on                 *
    7375  ! weight             weight for each concentration sample (1/2 or 1)         *
    7476  ! uap(maxpart),ucp(maxpart),uzp(maxpart) = random velocities due to          *
    … …  
    9294  use par_mod
    9395  use com_mod
     96#ifdef NETCDF_OUTPUT
     97  use netcdf_output_mod, only: concoutput_netcdf, concoutput_nest_netcdf,concoutput_surf_netcdf, concoutput_surf_nest_netcdf
     98#endif
    9499
    95100  implicit none
    … …  
    129134
    130135
     136  itime=0
     137  !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
     138  write(*,46) float(itime)/3600,itime,numpart
    131139  if (verbosity.gt.0) then
    132140    write (*,*) 'timemanager> starting simulation'
    … …  
    138146
    139147  do itime=0,ideltas,lsynctime
    140 
    141148
    142149  ! Computation of wet deposition, OH reaction and mass transfer
    … …  
    150157  !********************************************************************
    151158
    152     if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
     159    if (wetdep .and. itime .ne. 0 .and. numpart .gt. 0) then
    153160        if (verbosity.gt.0) then
    154161           write (*,*) 'timemanager> call wetdepo'
    … …  
    157164    endif
    158165
    159     if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) &
     166    if (ohrea .and. itime .ne. 0 .and. numpart .gt. 0) &
    160167         call ohreaction(itime,lsynctime,loutnext)
    161168
    162     if (ASSSPEC .and. itime .ne. 0 .and. numpart .gt. 0) then
     169    if (assspec .and. itime .ne. 0 .and. numpart .gt. 0) then
    163170       stop 'associated species not yet implemented!'
    164171  !     call transferspec(itime,lsynctime,loutnext)
    … …  
    238245  !***********************************************************************
    239246
    240     if (DEP.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
     247    if (dep.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
    241248      do ks=1,nspec
    242249      do kp=1,maxpointspec_act
    … …  
    348355        if ((iout.le.3.).or.(iout.eq.5)) then
    349356          if (surf_only.ne.1) then
    350           call concoutput(itime,outnum,gridtotalunc, &
    351                wetgridtotalunc,drygridtotalunc)
     357            if (lnetcdfout.eq.1) then
     358#ifdef NETCDF_OUTPUT
     359              call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
     360#endif
     361            else
     362              call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
     363            endif
    352364          else 
    353   if (verbosity.eq.1) then
    354      print*,'call concoutput_surf '
    355      CALL SYSTEM_CLOCK(count_clock)
    356      WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
    357   endif
    358           call concoutput_surf(itime,outnum,gridtotalunc, &
    359                wetgridtotalunc,drygridtotalunc)
    360   if (verbosity.eq.1) then
    361      print*,'called concoutput_surf '
    362      CALL SYSTEM_CLOCK(count_clock)
    363      WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
    364   endif
     365            if (verbosity.eq.1) then
     366             print*,'call concoutput_surf '
     367             call system_clock(count_clock)
     368             write(*,*) 'system clock',count_clock - count_clock0   
     369            endif
     370            if (lnetcdfout.eq.1) then
     371#ifdef NETCDF_OUTPUT
     372              call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
     373#endif
     374            else
     375              call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
     376              if (verbosity.eq.1) then
     377                print*,'called concoutput_surf '
     378                call system_clock(count_clock)
     379                write(*,*) 'system clock',count_clock - count_clock0   
     380              endif
     381            endif
    365382          endif
    366383
    367           if ((nested_output.eq.1).and.(surf_only.ne.1)) call concoutput_nest(itime,outnum)
    368           if ((nested_output.eq.1).and.(surf_only.eq.1)) call concoutput_surf_nest(itime,outnum)
     384          if (nested_output .eq. 1) then
     385            if (lnetcdfout.eq.0) then
     386              if (surf_only.ne.1) then
     387                call concoutput_nest(itime,outnum)
     388              else
     389                call concoutput_surf_nest(itime,outnum)
     390              endif
     391            else
     392#ifdef NETCDF_OUTPUT
     393              if (surf_only.ne.1) then
     394                call concoutput_nest_netcdf(itime,outnum)
     395              else
     396                call concoutput_surf_nest_netcdf(itime,outnum)
     397              endif
     398#endif
     399            endif
     400          endif
    369401          outnum=0.
    370402        endif
    371403        if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
    372404        if (iflux.eq.1) call fluxoutput(itime)
    373         write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc, &
    374              drygridtotalunc
    375 45      format(i9,' SECONDS SIMULATED: ',i8, &
    376              ' PARTICLES:    Uncertainty: ',3f7.3)
     405        !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
     406        write(*,46) float(itime)/3600,itime,numpart
     40745      format(i9,' SECONDS SIMULATED: ',i8, ' PARTICLES:    Uncertainty: ',3f7.3)
     40846      format(' Simulated ',f7.1,' hours (',i9,' s), ',i8, ' particles')
    377409        if (ipout.ge.1) call partoutput(itime)    ! dump particle positions
    378410        loutnext=loutnext+loutstep
    … …  
    472504  !****************************
    473505
    474         xold=xtra1(j)
    475         yold=ytra1(j)
     506        xold=real(xtra1(j))
     507        yold=real(ytra1(j))
    476508        zold=ztra1(j)
    477509
    … …  
    513545            endif
    514546
    515             if (DRYDEPSPEC(ks)) then        ! dry deposition
     547            if (drydepspec(ks)) then        ! dry deposition
    516548              drydeposit(ks)=xmass1(j,ks)*prob(ks)*decfact
    517549              xmass1(j,ks)=xmass1(j,ks)*(1.-prob(ks))*decfact
    … …  
    524556            endif
    525557
    526 
    527558            if (mdomainfill.eq.0) then
    528559              if (xmass(npoint(j),ks).gt.0.) &
    … …  
    536567          if (xmassfract.lt.0.0001) then   ! terminate all particles carrying less mass
    537568            itra1(j)=-999999999
     569            if (verbosity.gt.0) then
     570              print*,'terminated particle ',j,' for small mass'
     571            endif
    538572          endif
    539573
    540574  !        Sabine Eckhardt, June 2008
    541575  !        don't create depofield for backward runs
    542           if (DRYDEP.AND.(ldirect.eq.1)) then
    543             call drydepokernel(nclass(j),drydeposit,real(xtra1(j)), &
    544                  real(ytra1(j)),nage,kp)
    545             if (nested_output.eq.1) call drydepokernel_nest( &
    546                  nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)), &
    547                  nage,kp)
     576          if (drydep.AND.(ldirect.eq.1)) then
     577            call drydepokernel(nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)),nage,kp)
     578            if (nested_output.eq.1) then
     579              call drydepokernel_nest(nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)),nage,kp)
     580            endif
    548581          endif
    549582
    … …  
    552585
    553586          if (abs(itra1(j)-itramem(j)).ge.lage(nageclass)) then
    554             if (linit_cond.ge.1) &
    555                  call initial_cond_calc(itime+lsynctime,j)
     587            if (linit_cond.ge.1) call initial_cond_calc(itime+lsynctime,j)
    556588            itra1(j)=-999999999
     589            if (verbosity.gt.0) then
     590              print*,'terminated particle ',j,' for age'
     591            endif
    557592          endif
    558593        endif
    … …  
    582617
    583618  if (iflux.eq.1) then
    584       deallocate(flux)
     619    deallocate(flux)
    585620  endif
    586   if (OHREA.eqv..TRUE.) then
    587       deallocate(OH_field,OH_field_height)
     621  if (ohrea.eqv..TRUE.) then
     622    deallocate(OH_field,OH_field_height)
    588623  endif
    589624  if (ldirect.gt.0) then
    590   deallocate(drygridunc,wetgridunc)
     625    deallocate(drygridunc,wetgridunc)
    591626  endif
    592627  deallocate(gridunc)
    … …  
    595630  deallocate(xmasssave)
    596631  if (nested_output.eq.1) then
    597      deallocate(orooutn, arean, volumen)
    598      if (ldirect.gt.0) then
    599      deallocate(griduncn,drygriduncn,wetgriduncn)
    600      endif
     632    deallocate(orooutn, arean, volumen)
     633    if (ldirect.gt.0) then
     634      deallocate(griduncn,drygriduncn,wetgriduncn)
     635    endif
    601636  endif
    602637  deallocate(outheight,outheighthalf)
    603   deallocate(oroout, area, volume)
     638  deallocate(oroout,area,volume)
    604639
    605640end subroutine timemanager

  • /trunk/src/verttransform.f90

    r20 r30  
    2121
    2222subroutine verttransform(n,uuh,vvh,wwh,pvh)
    23   !                         i  i   i   i   i
     23  !                      i  i   i   i   i
    2424  !*****************************************************************************
    2525  !                                                                            *
    … …  
    4949  ! Sabine Eckhardt, March 2007
    5050  ! added the variable cloud for use with scavenging - descr. in com_mod
    51   ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
    52   ! note that also other subroutines are affected by the fix
    5351  !*****************************************************************************
    5452  !                                                                            *
    … …  
    7270
    7371  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
    7673  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)
    8076  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
    8177  real :: xlon,ylat,xlonr,dzdx,dzdy
    82   real :: dzdx1,dzdx2,dzdy1,dzdy2, precmin
     78  real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
    8379  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
    8480  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
    … …  
    8783  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
    8884  real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
    89   logical lconvectprec
    9085  real,parameter :: const=r_air/ga
    91   parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
    9286
    9387  logical :: init = .true.
    … …  
    9791  ! If verttransform is called the first time, initialize heights of the   *
    9892  ! 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.                                                *
    10494  !*************************************************************************
    105 
    106 
    107   !      do 897 kz=1,nuvz
    108   !       write (*,*) 'akz: ',akz(kz),'bkz',bkz(kz)
    109   !897     continue
    11095
    11196  if (init) then
    … …  
    11398  ! Search for a point with high surface pressure (i.e. not above significant topography)
    11499  ! Then, use this point to construct a reference z profile, to be used at all times
    115   !*****************************************************************************
     100  !**************************************************************************************
    116101
    117102    do jy=0,nymin1
    … …  
    128113
    129114    tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
    130          ps(ixm,jym,1,n))
     115    ps(ixm,jym,1,n))
    131116    pold=ps(ixm,jym,1,n)
    132117    height(1)=0.
    … …  
    136121      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
    137122
    138 
    139   ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
    140   ! upon the transformation to z levels. In order to save computer memory, this is
    141   ! not done anymore in the standard version. However, this option can still be
    142   ! switched on by replacing the following lines with those below, that are
    143   ! 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 
    148123      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)
    152126      else
    153         height(kz)=height(kz-1)+ &
    154              const*log(pold/pint)*tv
     127        height(kz)=height(kz-1)+const*log(pold/pint)*tv
    155128      endif
    156 
    157   ! Switch on following lines to use doubled vertical resolution
    158   !*************************************************************
    159   !    if (abs(tv-tvold).gt.0.2) then
    160   !      height((kz-1)*2)=
    161   !    +      height(max((kz-2)*2,1))+const*log(pold/pint)*
    162   !    +      (tv-tvold)/log(tv/tvold)
    163   !    else
    164   !      height((kz-1)*2)=height(max((kz-2)*2,1))+
    165   !    +      const*log(pold/pint)*tv
    166   !    endif
    167   ! End doubled vertical resolution
    168129
    169130      tvold=tv
    170131      pold=pint
    171132    end do
    172 
    173 
    174   ! Switch on following lines to use doubled vertical resolution
    175   !*************************************************************
    176   !  do 7 kz=3,nz-1,2
    177   !    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 resolution
    180133
    181134
    … …  
    204157  do jy=0,nymin1
    205158    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))
    208160      pold=ps(ix,jy,1,n)
    209       uvzlev(1)=0.
     161      uvwzlev(ix,jy,1)=0.
    210162      wzlev(1)=0.
    211163      rhoh(1)=pold/(r_air*tvold)
    … …  
    221173
    222174        if (abs(tv-tvold).gt.0.2) then
    223           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
    224                (tv-tvold)/log(tv/tvold)
     175          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
     176          (tv-tvold)/log(tv/tvold)
    225177        else
    226           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
     178          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
    227179        endif
    228180
    … …  
    233185
    234186      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)
    252190
    253191  ! pinmconv=(h2-h1)/(p2-p1)
    254192
    255193      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
    256            ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
    257            (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
     194      ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
     195      (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
    258196      do kz=2,nz-1
    259197        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
    260              ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
    261              (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
     198        ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
     199        (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
    262200      end do
    263201      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
    264            ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
    265            (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
     202      ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
     203      (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
    266204
    267205  ! Levels, where u,v,t and q are given
    … …  
    283221      do iz=2,nz-1
    284222        do kz=kmin,nuvz
    285           if(height(iz).gt.uvzlev(nuvz)) then
     223          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
    286224            uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
    287225            vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
    … …  
    292230            goto 30
    293231          endif
    294           if ((height(iz).gt.uvzlev(kz-1)).and. &
    295                (height(iz).le.uvzlev(kz))) then
    296            dz1=height(iz)-uvzlev(kz-1)
    297            dz2=uvzlev(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)
    298236           dz=dz1+dz2
    299237           uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
    … …  
    339277
    340278      drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
    341            (height(2)-height(1))
     279      (height(2)-height(1))
    342280      do kz=2,nz-1
    343281        drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
    344              (height(kz+1)-height(kz-1))
     282        (height(kz+1)-height(kz-1))
    345283      end do
    346284      drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
    … …  
    356294
    357295  do jy=1,ny-2
     296    cosf=cos((real(jy)*dy+ylat0)*pi180)
    358297    do ix=1,nx-2
    359298
    … …  
    361300      do iz=2,nz-1
    362301
    363         ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
     302        ui=uu(ix,jy,iz,n)*dxconst/cosf
    364303        vi=vv(ix,jy,iz,n)*dyconst
    365304
    366305        do kz=kmin,nz
    367306          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
    368                (height(iz).le.uvwzlev(ix,jy,kz))) then
     307          (height(iz).le.uvwzlev(ix,jy,kz))) then
    369308            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
    370309            dz2=uvwzlev(ix,jy,kz)-height(iz)
    … …  
    409348        do iz=1,nz
    410349          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))
    413351        end do
    414352      end do
    … …  
    424362      xlon=xlon0+real(nx/2-1)*dx
    425363      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)
    428365      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
    431367      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
    432368        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
    433              vv(nx/2-1,nymin1,iz,n))-xlonr
     369        vv(nx/2-1,nymin1,iz,n))-xlonr
    434370      else
    435371        ddpol=pi/2-xlonr
    … …  
    444380      uuaux=-ffpol*sin(xlonr+ddpol)
    445381      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)
    449383      jy=nymin1
    450384      do ix=0,nxmin1
    … …  
    485419        do iz=1,nz
    486420          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))
    489422        end do
    490423      end do
    … …  
    499432      xlon=xlon0+real(nx/2-1)*dx
    500433      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)
    503435      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
    506437      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
    509439      else
    510440        ddpol=pi/2-xlonr
    … …  
    519449      uuaux=+ffpol*sin(xlonr-ddpol)
    520450      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)
    523452
    524453      jy=0
    … …  
    551480  !   create a cloud and rainout/washout field, clouds occur where rh>80%
    552481  !   total cloudheight is stored at level 0
    553 
    554 
    555 
    556482  do jy=0,nymin1
    557483    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
    620505            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
    643513            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
    657516    end do
    658517  end do
    659518
    660   !do 102 kz=1,nuvz
    661   !write(an,'(i02)') kz+10
    662   !write(*,*) nuvz,nymin1,nxmin1,'--',an,'--'
    663   !open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')
    664   !do 101 jy=0,nymin1
    665   !    write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)
    666   !101   continue
    667   ! close(4)
    668   !102   continue
    669 
    670   ! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')
    671   ! do 103 jy=0,nymin1
    672   !     write (4,*)
    673   !+       (height(kz),kz=1,nuvz)
    674   !103   continue
    675   ! close(4)
    676 
    677   !open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')
    678   ! do 104 jy=0,nymin1
    679   !     write (4,*)
    680   !+       (r_air*tt(ix,jy,1,n)*rho(ix,jy,1,n),ix=1,nxmin1)
    681   !104   continue
    682   ! close(4)
    683 
    684 
    685519end subroutine verttransform

  • /trunk/src/verttransform_gfs.f90

    r20 r30  
    2121
    2222subroutine verttransform(n,uuh,vvh,wwh,pvh)
    23   !                         i  i   i   i   i
     23  !                      i  i   i   i   i
    2424  !*****************************************************************************
    2525  !                                                                            *
    … …  
    7171  real :: f_qvsat,pressure
    7272  real :: rh,lsp,convp
    73   real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
     73  real :: rhoh(nuvzmax),pinmconv(nzmax)
    7474  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
    7575  real :: xlon,ylat,xlonr,dzdx,dzdy
    76   real :: dzdx1,dzdx2,dzdy1,dzdy2
     76  real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
    7777  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
    7878  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
    … …  
    9292  ! If verttransform is called the first time, initialize heights of the   *
    9393  ! 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.                                                *
    9995  !*************************************************************************
    10096
    … …  
    118114
    119115    tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
    120          ps(ixm,jym,1,n))
     116    ps(ixm,jym,1,n))
    121117    pold=ps(ixm,jym,1,n)
    122118    height(1)=0.
    … …  
    126122      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
    127123
    128 
    129   ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
    130   ! upon the transformation to z levels. In order to save computer memory, this is
    131   ! not done anymore in the standard version. However, this option can still be
    132   ! switched on by replacing the following lines with those below, that are
    133   ! 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 
    138124      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)
    142127      else
    143         height(kz)=height(kz-1)+ &
    144              const*log(pold/pint)*tv
     128        height(kz)=height(kz-1)+const*log(pold/pint)*tv
    145129      endif
    146 
    147   ! Switch on following lines to use doubled vertical resolution
    148   !*************************************************************
    149   !    if (abs(tv-tvold).gt.0.2) then
    150   !      height((kz-1)*2)=
    151   !    +      height(max((kz-2)*2,1))+const*log(pold/pint)*
    152   !    +      (tv-tvold)/log(tv/tvold)
    153   !    else
    154   !      height((kz-1)*2)=height(max((kz-2)*2,1))+
    155   !    +      const*log(pold/pint)*tv
    156   !    endif
    157   ! End doubled vertical resolution
    158130
    159131      tvold=tv
    160132      pold=pint
    161133    end do
    162 
    163 
    164   ! Switch on following lines to use doubled vertical resolution
    165   !*************************************************************
    166   !  do 7 kz=3,nz-1,2
    167   !    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 resolution
    170134
    171135
    … …  
    198162      llev = 0
    199163      do i=1,nuvz
    200        if (ps(ix,jy,1,n).lt.akz(i)) llev=i
     164        if (ps(ix,jy,1,n).lt.akz(i)) llev=i
    201165      end do
    202166       llev = llev+1
    … …  
    212176      tvold=tth(ix,jy,llev,n)*(1.+0.608*qvh(ix,jy,llev,n))
    213177      pold=akz(llev)
    214       uvzlev(llev)=0.
    215178      wzlev(llev)=0.
    216179      uvwzlev(ix,jy,llev)=0.
    … …  
    223186
    224187        if (abs(tv-tvold).gt.0.2) then
    225           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
    226                (tv-tvold)/log(tv/tvold)
     188          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
     189          (tv-tvold)/log(tv/tvold)
    227190        else
    228           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
     191          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
    229192        endif
    230         wzlev(kz)=uvzlev(kz)
    231         uvwzlev(ix,jy,kz)=uvzlev(kz)
     193        wzlev(kz)=uvwzlev(ix,jy,kz)
    232194
    233195        tvold=tv
    234196        pold=pint
    235197      end do
    236 
    237 
    238   ! Switch on following lines to use doubled vertical resolution
    239   ! Switch off the three lines above.
    240   !*************************************************************
    241   !22          uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
    242   !     do 23 kz=2,nwz
    243   !23          uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
    244   ! End doubled vertical resolution
    245198
    246199  ! pinmconv=(h2-h1)/(p2-p1)
    … …  
    279232      do iz=2,nz-1
    280233        do kz=kmin,nuvz
    281           if(height(iz).gt.uvzlev(nuvz)) then
     234          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
    282235            uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
    283236            vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
    … …  
    289242            goto 30
    290243          endif
    291           if ((height(iz).gt.uvzlev(kz-1)).and. &
    292                (height(iz).le.uvzlev(kz))) then
    293            dz1=height(iz)-uvzlev(kz-1)
    294            dz2=uvzlev(kz)-height(iz)
    295            dz=dz1+dz2
    296            uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
    297            vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
    298            tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
    299                 +tth(ix,jy,kz,n)*dz1)/dz
    300            qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
    301                 +qvh(ix,jy,kz,n)*dz1)/dz
    302            pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
    303            rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
    304            pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz
     244          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
    305258          endif
    306259        end do
    … …  
    318271        do kz=kmin,nwz
    319272          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
    327279          endif
    328280        end do
    … …  
    337289      do kz=2,nz-1
    338290        drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
    339              (height(kz+1)-height(kz-1))
     291        (height(kz+1)-height(kz-1))
    340292      end do
    341293      drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
    … …  
    351303
    352304  do jy=1,ny-2
     305    cosf=cos((real(jy)*dy+ylat0)*pi180)
    353306    do ix=1,nx-2
    354307
    … …  
    367320      do iz=2,nz-1
    368321
    369         ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
     322        ui=uu(ix,jy,iz,n)*dxconst/cosf
    370323        vi=vv(ix,jy,iz,n)*dyconst
    371324
    372325        do kz=kmin,nz
    373326          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
    374                (height(iz).le.uvwzlev(ix,jy,kz))) then
     327          (height(iz).le.uvwzlev(ix,jy,kz))) then
    375328            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
    376329            dz2=uvwzlev(ix,jy,kz)-height(iz)
    … …  
    426379      xlon=xlon0+real(nx/2-1)*dx
    427380      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
    433384      elseif (vv(nx/2-1,nymin1,iz,n).gt.0.) then
    434385        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
    435              vv(nx/2-1,nymin1,iz,n))-xlonr
     386        vv(nx/2-1,nymin1,iz,n))-xlonr
    436387      else
    437388        ddpol=pi/2-xlonr
    … …  
    446397      uuaux=-ffpol*sin(xlonr+ddpol)
    447398      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)
    451400      jy=nymin1
    452401      do ix=0,nxmin1
    … …  
    460409  ! ward parallel of latitude
    461410
    462   do iz=1,nz
     411    do iz=1,nz
    463412      wdummy=0.
    464413      jy=ny-2
    … …  
    471420        ww(ix,jy,iz,n)=wdummy
    472421      end do
    473   end do
     422    end do
    474423
    475424  endif
    … …  
    487436        do iz=1,nz
    488437          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))
    491439        end do
    492440      end do
    … …  
    498446      xlon=xlon0+real(nx/2-1)*dx
    499447      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)
    502449      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
    505451      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
    508453      else
    509454        ddpol=pi/2-xlonr
    … …  
    518463      uuaux=+ffpol*sin(xlonr-ddpol)
    519464      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)
    522466
    523467      jy=0
    … …  
    562506         clouds(ix,jy,kz,n)=0
    563507         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
    576519         else ! no cloud
    577             if (rain_cloud_above.eq.1) then ! scavenging
    578                if (lsp.ge.convp) then
    579                   clouds(ix,jy,kz,n)=5 ! lsp dominated washout
    580                else
    581                   clouds(ix,jy,kz,n)=4 ! convp dominated washout
    582                endif
    583             endif
     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
    584527         endif
    585528      end do

  • /trunk/src/verttransform_nests.f90

    r20 r30  
    2121
    2222subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn)
    23   !                               i   i    i    i   i
     23  !                            i   i    i    i   i
    2424  !*****************************************************************************
    2525  !                                                                            *
    … …  
    7272  integer :: rain_cloud_above,kz_inv
    7373  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)
    7575  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
    7777  real :: dzdx,dzdy
    7878  real :: dzdx1,dzdx2,dzdy1,dzdy2
    … …  
    9696
    9797      tvold=tt2n(ix,jy,1,n,l)*(1.+0.378*ew(td2n(ix,jy,1,n,l))/ &
    98            psn(ix,jy,1,n,l))
     98      psn(ix,jy,1,n,l))
    9999      pold=psn(ix,jy,1,n,l)
    100       uvzlev(1)=0.
     100      uvwzlev(ix,jy,1)=0.
    101101      wzlev(1)=0.
    102102      rhoh(1)=pold/(r_air*tvold)
    … …  
    112112
    113113        if (abs(tv-tvold).gt.0.2) then
    114           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
    115                (tv-tvold)/log(tv/tvold)
     114          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
     115          (tv-tvold)/log(tv/tvold)
    116116        else
    117           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
     117          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
    118118        endif
    119119
    … …  
    124124
    125125      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)
    151129
    152130  ! pinmconv=(h2-h1)/(p2-p1)
    153131
    154132      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
    155            ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
    156            (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
     133      ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
     134      (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
    157135      do kz=2,nz-1
    158136        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
    159              ((aknew(kz+1)+bknew(kz+1)*psn(ix,jy,1,n,l))- &
    160              (aknew(kz-1)+bknew(kz-1)*psn(ix,jy,1,n,l)))
     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)))
    161139      end do
    162140      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
    163            ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
    164            (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
     141      ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
     142      (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
    165143
    166144
    … …  
    183161      do iz=2,nz-1
    184162        do kz=kmin,nuvz
    185           if(height(iz).gt.uvzlev(nuvz)) then
     163          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
    186164            uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
    187165            vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
    … …  
    192170            goto 30
    193171          endif
    194           if ((height(iz).gt.uvzlev(kz-1)).and. &
    195                (height(iz).le.uvzlev(kz))) then
    196            dz1=height(iz)-uvzlev(kz-1)
    197            dz2=uvzlev(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)
    198176           dz=dz1+dz2
    199177           uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
    200                 uuhn(ix,jy,kz,l)*dz1)/dz
     178           uuhn(ix,jy,kz,l)*dz1)/dz
    201179           vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
    202                 vvhn(ix,jy,kz,l)*dz1)/dz
     180           vvhn(ix,jy,kz,l)*dz1)/dz
    203181           ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
    204                 tthn(ix,jy,kz,n,l)*dz1)/dz
     182           tthn(ix,jy,kz,n,l)*dz1)/dz
    205183           qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
    206                 qvhn(ix,jy,kz,n,l)*dz1)/dz
     184           qvhn(ix,jy,kz,n,l)*dz1)/dz
    207185           pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
    208                 pvhn(ix,jy,kz,l)*dz1)/dz
     186           pvhn(ix,jy,kz,l)*dz1)/dz
    209187           rhon(ix,jy,iz,n,l)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
    210188           kmin=kz
    … …  
    225203        do kz=kmin,nwz
    226204          if ((height(iz).gt.wzlev(kz-1)).and. &
    227                (height(iz).le.wzlev(kz))) then
    228            dz1=height(iz)-wzlev(kz-1)
    229            dz2=wzlev(kz)-height(iz)
    230            dz=dz1+dz2
    231            wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
    232                 +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
    233            kmin=kz
    234            goto 40
     205          (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
    235213          endif
    236214        end do
    … …  
    242220
    243221      drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
    244            (height(2)-height(1))
     222      (height(2)-height(1))
    245223      do kz=2,nz-1
    246224        drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
    247              rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
     225        rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
    248226      end do
    249227      drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
    … …  
    259237
    260238  do jy=1,nyn(l)-2
     239    cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
    261240    do ix=1,nxn(l)-2
    262241
    … …  
    264243      do iz=2,nz-1
    265244
    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
    268246        vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)
    269247
    … …  
    319297               rain_cloud_above=1
    320298               cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ &
    321                     height(kz)-height(kz-1)
     299               height(kz)-height(kz-1)
    322300               if (lsp.ge.convp) then
    323301                  cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout

  • /trunk/src/wetdepo.f90

    r20 r30  
    2121
    2222subroutine wetdepo(itime,ltsample,loutnext)
    23   !                     i      i        i
     23  !                  i      i        i
    2424  !*****************************************************************************
    2525  !                                                                            *
    … …  
    3939  ! Code may not be correct for decay of deposition!                           *
    4040  !                                                                            *
    41   ! Modification by Sabine Eckhart to introduce a new in-/below-cloud
    42   ! scheme, not dated
    43   ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
    4441  !*****************************************************************************
    4542  !                                                                            *
    … …  
    7471
    7572  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
    8075  real :: S_i, act_temp, cl, cle ! in cloud scavenging
    8176  real :: clouds_h ! cloud height for the specific grid point
    82   real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav, precsub,f
     77  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
    8378  real :: wetdeposit(maxspec),restmass
    8479  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
    … …  
    10398
    10499  do jpart=1,numpart
     100
    105101    if (itra1(jpart).eq.-999999999) goto 20
    106102    if(ldirect.eq.1)then
    … …  
    114110      if (itage.lt.lage(nage)) goto 33
    115111    end do
    116 33   continue
     11233  continue
    117113
    118114
    … …  
    123119    do j=numbnests,1,-1
    124120      if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
    125            (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
     121      (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
    126122        ngrid=j
    127123        goto 23
    128124      endif
    129125    end do
    130 23   continue
     12623  continue
    131127
    132128
    … …  
    151147    interp_time=nint(itime-0.5*ltsample)
    152148
    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
    182160  ! 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
     16726  continue
     168
     169    n=memind(2)
     170    if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
     171    n=memind(1)
    195172
    196173  ! if there is no precipitation or the particle is above the clouds no
    197174  ! scavenging is done
    198175
    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'
    226187
    227188  ! 1) Parameterization of the the area fraction of the grid cell where the
    … …  
    267228  !**********************************************************
    268229
    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
    293256
    294257! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening
    … …  
    297260! wetc_in=0.9 (default)
    298261! 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
    309268        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
    314284        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
    414316  ! Add the wet deposition to accumulated amount on output grid and nested output grid
    415317  !*****************************************************************************
    416318
    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
    434325
    43532620  continue
    436   end do
     327  end do ! all particles
    437328
    438329end subroutine wetdepo

  • /trunk/src/writeheader_nest.f90

    r20 r30  
    6767
    6868  if (ldirect.eq.1) then
    69     write(unitheader) ibdate,ibtime,'FLEXPART V8.2'
     69    write(unitheader) ibdate,ibtime, trim(flexversion)
    7070  else
    71     write(unitheader) iedate,ietime,'FLEXPART V8.2'
     71    write(unitheader) iedate,ietime, trim(flexversion)
    7272  endif
    7373
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG