Changes in / [30:20]


Ignore:
Files:
1 added
5 deleted
57 edited

Legend:

Unmodified
Added
Removed

  • /trunk/options/COMMAND

    r30 r20  
    1 ********************************************************************************
    2 *                                                                              *
    3 *      Input file for the Lagrangian particle dispersion model FLEXPART        *
    4 *                           Please select your options                         *
    5 *                                                                              *
    6 ********************************************************************************
    7 
    8 1. __                3X, I2
    9     1       
    10    LDIRECT           1 FOR FORWARD SIMULATION, -1 FOR BACKWARD SIMULATION
    11 
    12 2. ________ ______   3X, I8, 1X, I6
    13    20111210 000000
    14    YYYYMMDD HHMISS   BEGINNING DATE OF SIMULATION
    15 
    16 3. ________ ______   3X, I8, 1X, I6
    17    20111210 120000
    18    YYYYMMDD HHMISS   ENDING DATE OF SIMULATION
    19 
    20 4. _____             3X, I5
    21    10800
    22    SSSSS             OUTPUT EVERY SSSSS SECONDS
    23 
    24 5. _____             3X, I5
    25    10800
    26    SSSSS             TIME AVERAGE OF OUTPUT (IN SSSSS SECONDS)
    27 
    28 6. _____             3X, I5
    29      900
    30    SSSSS             SAMPLING RATE OF OUTPUT (IN SSSSS SECONDS)
    31 
    32 7. _________         3X, I9
    33    999999999
    34    SSSSSSSSS         TIME CONSTANT FOR PARTICLE SPLITTING (IN SECONDS)
    35 
    36 8. _____             3X, I5
    37      900
    38    SSSSS             SYNCHRONISATION INTERVAL OF FLEXPART (IN SECONDS)
    39 
    40 9.  ---.--           4X, F6.4
    41      -5.0
    42     CTL              FACTOR, BY WHICH TIME STEP MUST BE SMALLER THAN TL
    43 
    44 10. ---              4X, I3
    45       4
    46     IFINE            DECREASE OF TIME STEP FOR VERTICAL MOTION BY FACTOR IFINE
    47 
    48 11. -                4X, I1
    49     3 
    50     IOUT             1 CONCENTRATION (RESIDENCE TIME FOR BACKWARD RUNS) OUTPUT, 2 MIXING RATIO OUTPUT, 3 BOTH,4 PLUME TRAJECT., 5=1+4
    51 
    52 12. -                4X, I1
    53     0 
    54     IPOUT            PARTICLE DUMP: 0 NO, 1 EVERY OUTPUT INTERVAL, 2 ONLY AT END
    55 
    56 13. _                4X, I1
    57     1
    58     LSUBGRID         SUBGRID TERRAIN EFFECT PARAMETERIZATION: 1 YES, 0 NO
    59 
    60 14. _                4X, I1
    61     1
    62     LCONVECTION      CONVECTION: 1 YES, 0 NO
    63 
    64 15. _                4X, I1
    65     0
    66     LAGESPECTRA      AGE SPECTRA: 1 YES, 0 NO
    67 
    68 16. _                4X, I1
    69     0
    70     IPIN             CONTINUE SIMULATION WITH DUMPED PARTICLE DATA: 1 YES, 0 NO
    71 
    72 17. _               
    73     0                4X,I1
    74     IOFR             IOUTPUTFOREACHREL CREATE AN OUPUT FILE FOR EACH RELEASE LOCATION: 1 YES, 0 NO
    75 
    76 18. _                4X, I1
    77     0
    78     IFLUX            CALCULATE FLUXES: 1 YES, 0 NO
    79 
    80 19. _                4X, I1
    81     0
    82     MDOMAINFILL      DOMAIN-FILLING TRAJECTORY OPTION: 1 YES, 0 NO, 2 STRAT. O3 TRACER
    83 
    84 20. _                4X, I1
    85     1
    86     IND_SOURCE       1=MASS UNIT , 2=MASS MIXING RATIO UNIT
    87 
    88 21. _                4X, I1
    89     1
    90     IND_RECEPTOR     1=MASS UNIT , 2=MASS MIXING RATIO UNIT
    91 
    92 22. _                4X, I1
    93     0
    94     MQUASILAG        QUASILAGRANGIAN MODE TO TRACK INDIVIDUAL PARTICLES: 1 YES, 0 NO
    95 
    96 23. _                4X, I1
    97     0
    98     NESTED_OUTPUT    SHALL NESTED OUTPUT BE USED? 1 YES, 0 NO
    99 
    100 24. _                4X, I1
    101     2
    102     LINIT_COND       INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT
    103 
    104 25. _                4X, I1
    105     0
    106     SURF_ONLY        IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER
    107 
    108 
    109 1. Simulation direction, 1 for forward, -1 for backward in time
    110         (consult Seibert and Frank, 2004 for backward runs)
    111 
    112 2. Beginning date and time of simulation. Must be given in format
    113    YYYYMMDD HHMISS, where YYYY is YEAR, MM is MONTH, DD is DAY, HH is HOUR,
    114    MI is MINUTE and SS is SECOND. Current version utilizes UTC.
    115 
    116 3. Ending date and time of simulation. Same format as 3.
    117 
    118 4. Average concentrations are calculated every SSSSS seconds.
    119 
    120 5. The average concentrations are time averages of SSSSS seconds
    121    duration. If SSSSS is 0, instantaneous concentrations are outputted.
    122 
    123 6. The concentrations are sampled every SSSSS seconds to calculate the time
    124    average concentration. This period must be shorter than the averaging time.
    125 
    126 7. Time constant for particle splitting. Particles are split into two
    127    after SSSSS seconds, 2xSSSSS seconds, 4xSSSSS seconds, and so on.
    128 
    129 8. All processes are synchronized with this time interval (lsynctime).
    130    Therefore, all other time constants must be multiples of this value.
    131    Output interval and time average of output must be at least twice lsynctime.
    132 
    133 9. CTL must be >1 for time steps shorter than the  Lagrangian time scale
    134    If CTL<0, a purely random walk simulation is done
    135 
    136 10.IFINE=Reduction factor for time step used for vertical wind
    137 
    138 11.IOUT determines how the output shall be made: concentration
    139    (ng/m3, Bq/m3), mixing ratio (pptv), or both, or plume trajectory mode,
    140    or concentration + plume trajectory mode.
    141    In plume trajectory mode, output is in the form of average trajectories.
    142 
    143 12.IPOUT determines whether particle positions are outputted (in addition
    144    to the gridded concentrations or mixing ratios) or not.
    145    0=no output, 1 output every output interval, 2 only at end of the
    146    simulation
    147 
    148 13.Switch on/off subgridscale terrain parameterization (increase of
    149    mixing heights due to subgridscale orographic variations)
    150 
    151 14.Switch on/off the convection parameterization
    152 
    153 15.Switch on/off the calculation of age spectra: if yes, the file AGECLASSES
    154    must be available
    155 
    156 16. If IPIN=1, a file "partposit_end" from a previous run must be available in
    157     the output directory. Particle positions are read in and previous simulation
    158     is continued. If IPIN=0, no particles from a previous run are used
    159 
    160 17. IF IOUTPUTFOREACHRELEASE is set to 1, one output field for each location
    161     in the RLEASE file is created. For backward calculation this should be
    162     set to 1. For forward calculation both possibilities are applicable.
    163 
    164 18. If IFLUX is set to 1, fluxes of each species through each of the output
    165     boxes are calculated. Six fluxes, corresponding to northward, southward,
    166     eastward, westward, upward and downward are calculated for each grid cell of
    167     the output grid. The control surfaces are placed in the middle of each
    168     output grid cell. If IFLUX is set to 0, no fluxes are determined.
    169 
    170 19. If MDOMAINFILL is set to 1, the first box specified in file RELEASES is used
    171     as the domain where domain-filling trajectory calculations are to be done.
    172     Particles are initialized uniformly distributed (according to the air mass
    173     distribution) in that domain at the beginning of the simulation, and are
    174     created at the boundaries throughout the simulation period.
    175 
    176 20. IND_SOURCE switches between different units for concentrations at the source
    177     NOTE that in backward simulations the release of computational particles
    178     takes place at the "receptor" and the sampling of particles at the "source".
    179           1=mass units (for bwd-runs = concentration)
    180           2=mass mixing ratio units
    181 21. IND_RECEPTOR switches between different units for concentrations at the receptor
    182           1=mass units (concentrations)
    183           2=mass mixing ratio units
    184 
    185 22. MQUASILAG indicates whether particles shall be numbered consecutively (1) or
    186     with their release location number (0). The first option allows tracking of
    187     individual particles using the partposit output files
    188 
    189 23. NESTED_OUTPUT decides whether model output shall be made also for a nested
    190     output field (normally with higher resolution)
    191 
    192 24. LINIT_COND determines whether, for backward runs only, the sensitivity to initial
    193     conditions shall be calculated and written to output files
    194     0=no output, 1 or 2 determines in which units the initial conditions are provided.
    195 
    196 25. SURF_ONLY: When set to 1, concentration/emission sensitivity is written out only
    197     for the surface layer; useful for instance when only footprint emission sensitivity is needed
    198     but initial conditions are needed on a full 3-D grid
     1+++++++++++++ HEADER +++++++++++++++++
     2+++++++++++++ HEADER +++++++++++++++++
     3+++++++++++++ HEADER +++++++++++++++++
     4+++++++++++++ HEADER +++++++++++++++++
     5+++++++++++++ HEADER +++++++++++++++++
     6+++++++++++++ HEADER +++++++++++++++++
     7+++++++++++++ HEADER +++++++++++++++++
     8-1
     920070930 180000
     1020070931      0
     11  3600
     12  3600
     13   900
     14 9999999
     15900   SYNC
     16-5.0  CTL
     174     IFINE
     185     IOUT
     190     IPOUT
     201     LSUBGRID
     211     LCONVECTION
     221     LAGESPECTRA
     230     IPIN
     241     IOFR
     250     IFLUX
     260     MDOMAINFILL
     271     IND_SOURCE
     282     IND_RECEPTOR
     290     MQUASILAG
     301     NESTED_OUTPUT
     312     LINIT_COND        INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT

  • /trunk/options/COMMAND.alternative

    r30 r20  
    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
    32 0               SURF_ONLY         IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER
    3332
    3433
    … …  
    118117    conditions shall be calculated and written to output files
    119118    0=no output, 1 or 2 determines in which units the initial conditions are provided.
    120 
    121 25. SURF_ONLY: When set to 1, concentration/emission sensitivity is written out only
    122     for the surface layer; useful for instance when only footprint emission sensitivity is needed
    123     but initial conditions are needed on a full 3-D grid

  • /trunk/options/COMMAND.reference

    r30 r20  
    1111
    12122. ________ ______   3X, I8, 1X, I6
    13    20110310 000000
     13   20040720 000000
    1414   YYYYMMDD HHMISS   BEGINNING DATE OF SIMULATION
    1515
    16163. ________ ______   3X, I8, 1X, I6
    17    20110310 120000
     17   20040721 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 
    104 25. _                4X, I1
    105     0
    106     SURF_ONLY        IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER
    107103
    108104
    … …  
    193189    conditions shall be calculated and written to output files
    194190    0=no output, 1 or 2 determines in which units the initial conditions are provided.
    195 
    196 25. SURF_ONLY: When set to 1, concentration/emission sensitivity is written out only
    197     for the surface layer; useful for instance when only footprint emission sensitivity is needed
    198     but initial conditions are needed on a full 3-D grid

  • /trunk/options/RELEASES

    r30 r20  
    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
    14 
    15   3
    16 ___                        i3    Index of species in file SPECIES
    17 
    18 =========================================================================
    19 20111210 000001
    20 ________ ______            i8,1x,i6 Beginning date and time of release
    21 
    22 20111210 090000
    23 ________ ______            i8,1x,i6 Ending date and time of release
    24 
    25    9.4048
    26 ____.____                  f9.4  Longitude [DEG] of lower left corner
    27 
    28   48.5060
    29 ____.____                  f9.4  Latitude [DEG] of lower left corner
    30 
    31    9.5067
    32 ____.____                  f9.4  Longitude [DEG] of upper right corner
    33 
    34   48.5158
    35 ____.____                  f9.4  Latitude [DEG] of upper right corner
    36 
    37         2
    38 _________                  i9    1 for m above ground, 2 for m above sea level
    39 
    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
    48 
    49 1.0000E00
    50 _.____E__                  e9.4  Total mass emitted
    51 
    52 RELEASE_TEST1
    53 ________________________________________   character*40 comment
    54 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     1+++++++++++++++++ HEADER ++++++++++++++++++++
     2+++++++++++++++++ HEADER ++++++++++++++++++++
     3+++++++++++++++++ HEADER ++++++++++++++++++++
     4+++++++++++++++++ HEADER ++++++++++++++++++++
     5+++++++++++++++++ HEADER ++++++++++++++++++++
     6+++++++++++++++++ HEADER ++++++++++++++++++++
     7+++++++++++++++++ HEADER ++++++++++++++++++++
     8+++++++++++++++++ HEADER ++++++++++++++++++++
     9+++++++++++++++++ HEADER ++++++++++++++++++++
     10+++++++++++++++++ HEADER ++++++++++++++++++++
     11+++++++++++++++++ HEADER ++++++++++++++++++++
     125
     1331
     1431
     1531
     1631
     1731
     18
     1920040101 080000
     2020040102 080000
     21   8.2500
     22  58.3833
     23   8.2500
     24  58.3833
     251
     26     0.000
     27   100.000
     284000
     291
     301
     311
     321
     331
     34rel_bj
     35
     3620040108 080000
     3720040109 080000
     38   8.2500
     39  58.3833
     40   8.2500
     41  58.3833
     421
     43     0.000
     44   100.000
     454000
     461
     471
     481
     491
     501
     51rel_bj
     52
     5320040115 080000
     5420040116 080000
     55   8.2500
     56  58.3833
     57   8.2500
     58  58.3833
     591
     60     0.000
     61   100.000
     624000
     631
     641
     651
     661
     671
     68rel_bj
     69
     7020040122 080000
     7120040123 080000
     72   8.2500
     73  58.3833
     74   8.2500
     75  58.3833
     761
     77     0.000
     78   100.000
     794000
     801
     811
     821
     831
     841
     85rel_bj
     86
     8720040129 080000
     8820040130 080000
     89   8.2500
     90  58.3833
     91   8.2500
     92  58.3833
     931
     94     0.000
     95   100.000
     964000
     971
     981
     991
     1001
     1011
     102rel_bj
     103
     10420040205 080000
     10520040206 080000
     106   8.2500
     107  58.3833
     108   8.2500
     109  58.3833
     1101
     111     0.000
     112   100.000
     1134000
     1141
     1151
     1161
     1171
     1181
     119rel_bj
     120
     12120040212 080000
     12220040213 080000
     123   8.2500
     124  58.3833
     125   8.2500
     126  58.3833
     1271
     128     0.000
     129   100.000
     1304000
     1311
     1321
     1331
     1341
     1351
     136rel_bj
     137
     13820040219 080000
     13920040220 080000
     140   8.2500
     141  58.3833
     142   8.2500
     143  58.3833
     1441
     145     0.000
     146   100.000
     1474000
     1481
     1491
     1501
     1511
     1521
     153rel_bj
     154
     15520040226 080000
     15620040227 080000
     157   8.2500
     158  58.3833
     159   8.2500
     160  58.3833
     1611
     162     0.000
     163   100.000
     1644000
     1651
     1661
     1671
     1681
     1691
     170rel_bj
     171
     17220040304 080000
     17320040305 080000
     174   8.2500
     175  58.3833
     176   8.2500
     177  58.3833
     1781
     179     0.000
     180   100.000
     1814000
     1821
     1831
     1841
     1851
     1861
     187rel_bj
     188
     18920040311 080000
     19020040312 080000
     191   8.2500
     192  58.3833
     193   8.2500
     194  58.3833
     1951
     196     0.000
     197   100.000
     1984000
     1991
     2001
     2011
     2021
     2031
     204rel_bj
     205
     20620040318 080000
     20720040319 080000
     208   8.2500
     209  58.3833
     210   8.2500
     211  58.3833
     2121
     213     0.000
     214   100.000
     2154000
     2161
     2171
     2181
     2191
     2201
     221rel_bj
     222
     22320040325 080000
     22420040326 080000
     225   8.2500
     226  58.3833
     227   8.2500
     228  58.3833
     2291
     230     0.000
     231   100.000
     2324000
     2331
     2341
     2351
     2361
     2371
     238rel_bj
     239
     24020040401 080000
     24120040402 080000
     242   8.2500
     243  58.3833
     244   8.2500
     245  58.3833
     2461
     247     0.000
     248   100.000
     2494000
     2501
     2511
     2521
     2531
     2541
     255rel_bj
     256
     25720040408 080000
     25820040409 080000
     259   8.2500
     260  58.3833
     261   8.2500
     262  58.3833
     2631
     264     0.000
     265   100.000
     2664000
     2671
     2681
     2691
     2701
     2711
     272rel_bj
     273
     27420040415 080000
     27520040416 080000
     276   8.2500
     277  58.3833
     278   8.2500
     279  58.3833
     2801
     281     0.000
     282   100.000
     2834000
     2841
     2851
     2861
     2871
     2881
     289rel_bj
     290
     29120040422 080000
     29220040423 080000
     293   8.2500
     294  58.3833
     295   8.2500
     296  58.3833
     2971
     298     0.000
     299   100.000
     3004000
     3011
     3021
     3031
     3041
     3051
     306rel_bj
     307
     30820040429 080000
     30920040430 080000
     310   8.2500
     311  58.3833
     312   8.2500
     313  58.3833
     3141
     315     0.000
     316   100.000
     3174000
     3181
     3191
     3201
     3211
     3221
     323rel_bj
     324
     32520040609 080000
     32620040610 080000
     327   8.2500
     328  58.3833
     329   8.2500
     330  58.3833
     3311
     332     0.000
     333   100.000
     3344000
     3351
     3361
     3371
     3381
     3391
     340rel_bj
     341
     34220040610 080000
     34320040611 080000
     344   8.2500
     345  58.3833
     346   8.2500
     347  58.3833
     3481
     349     0.000
     350   100.000
     3514000
     3521
     3531
     3541
     3551
     3561
     357rel_bj
     358
     35920040617 080000
     36020040618 080000
     361   8.2500
     362  58.3833
     363   8.2500
     364  58.3833
     3651
     366     0.000
     367   100.000
     3684000
     3691
     3701
     3711
     3721
     3731
     374rel_bj
     375
     37620040624 080000
     37720040625 080000
     378   8.2500
     379  58.3833
     380   8.2500
     381  58.3833
     3821
     383     0.000
     384   100.000
     3854000
     3861
     3871
     3881
     3891
     3901
     391rel_bj
     392
     39320040701 080000
     39420040702 080000
     395   8.2500
     396  58.3833
     397   8.2500
     398  58.3833
     3991
     400     0.000
     401   100.000
     4024000
     4031
     4041
     4051
     4061
     4071
     408rel_bj
     409
     41020040708 080000
     41120040709 080000
     412   8.2500
     413  58.3833
     414   8.2500
     415  58.3833
     4161
     417     0.000
     418   100.000
     4194000
     4201
     4211
     4221
     4231
     4241
     425rel_bj
     426
     42720040715 080000
     42820040716 080000
     429   8.2500
     430  58.3833
     431   8.2500
     432  58.3833
     4331
     434     0.000
     435   100.000
     4364000
     4371
     4381
     4391
     4401
     4411
     442rel_bj
     443
     44420040722 080000
     44520040723 080000
     446   8.2500
     447  58.3833
     448   8.2500
     449  58.3833
     4501
     451     0.000
     452   100.000
     4534000
     4541
     4551
     4561
     4571
     4581
     459rel_bj
     460
     46120040729 080000
     46220040730 080000
     463   8.2500
     464  58.3833
     465   8.2500
     466  58.3833
     4671
     468     0.000
     469   100.000
     4704000
     4711
     4721
     4731
     4741
     4751
     476rel_bj
     477
     47820040805 080000
     47920040806 080000
     480   8.2500
     481  58.3833
     482   8.2500
     483  58.3833
     4841
     485     0.000
     486   100.000
     4874000
     4881
     4891
     4901
     4911
     4921
     493rel_bj
     494
     49520040812 080000
     49620040813 080000
     497   8.2500
     498  58.3833
     499   8.2500
     500  58.3833
     5011
     502     0.000
     503   100.000
     5044000
     5051
     5061
     5071
     5081
     5091
     510rel_bj
     511
     51220040819 080000
     51320040820 080000
     514   8.2500
     515  58.3833
     516   8.2500
     517  58.3833
     5181
     519     0.000
     520   100.000
     5214000
     5221
     5231
     5241
     5251
     5261
     527rel_bj
     528
     52920040826 080000
     53020040827 080000
     531   8.2500
     532  58.3833
     533   8.2500
     534  58.3833
     5351
     536     0.000
     537   100.000
     5384000
     5391
     5401
     5411
     5421
     5431
     544rel_bj
     545
     54620040902 080000
     54720040903 080000
     548   8.2500
     549  58.3833
     550   8.2500
     551  58.3833
     5521
     553     0.000
     554   100.000
     5554000
     5561
     5571
     5581
     5591
     5601
     561rel_bj
     562
     56320040909 080000
     56420040910 080000
     565   8.2500
     566  58.3833
     567   8.2500
     568  58.3833
     5691
     570     0.000
     571   100.000
     5724000
     5731
     5741
     5751
     5761
     5771
     578rel_bj
     579
     58020040916 080000
     58120040917 080000
     582   8.2500
     583  58.3833
     584   8.2500
     585  58.3833
     5861
     587     0.000
     588   100.000
     5894000
     5901
     5911
     5921
     5931
     5941
     595rel_bj
     596
     59720040923 080000
     59820040924 080000
     599   8.2500
     600  58.3833
     601   8.2500
     602  58.3833
     6031
     604     0.000
     605   100.000
     6064000
     6071
     6081
     6091
     6101
     6111
     612rel_bj
     613
     61420040930 080000
     61520041001 080000
     616   8.2500
     617  58.3833
     618   8.2500
     619  58.3833
     6201
     621     0.000
     622   100.000
     6234000
     6241
     6251
     6261
     6271
     6281
     629rel_bj
     630
     63120041007 080000
     63220041008 080000
     633   8.2500
     634  58.3833
     635   8.2500
     636  58.3833
     6371
     638     0.000
     639   100.000
     6404000
     6411
     6421
     6431
     6441
     6451
     646rel_bj
     647
     64820041014 080000
     64920041015 080000
     650   8.2500
     651  58.3833
     652   8.2500
     653  58.3833
     6541
     655     0.000
     656   100.000
     6574000
     6581
     6591
     6601
     6611
     6621
     663rel_bj
     664
     66520041021 080000
     66620041022 080000
     667   8.2500
     668  58.3833
     669   8.2500
     670  58.3833
     6711
     672     0.000
     673   100.000
     6744000
     6751
     6761
     6771
     6781
     6791
     680rel_bj
     681
     68220041028 080000
     68320041029 080000
     684   8.2500
     685  58.3833
     686   8.2500
     687  58.3833
     6881
     689     0.000
     690   100.000
     6914000
     6921
     6931
     6941
     6951
     6961
     697rel_bj
     698
     69920041104 080000
     70020041105 080000
     701   8.2500
     702  58.3833
     703   8.2500
     704  58.3833
     7051
     706     0.000
     707   100.000
     7084000
     7091
     7101
     7111
     7121
     7131
     714rel_bj
     715
     71620041111 080000
     71720041112 080000
     718   8.2500
     719  58.3833
     720   8.2500
     721  58.3833
     7221
     723     0.000
     724   100.000
     7254000
     7261
     7271
     7281
     7291
     7301
     731rel_bj
     732
     73320041118 080000
     73420041119 080000
     735   8.2500
     736  58.3833
     737   8.2500
     738  58.3833
     7391
     740     0.000
     741   100.000
     7424000
     7431
     7441
     7451
     7461
     7471
     748rel_bj
     749
     75020041124 080000
     75120041125 080000
     752   8.2500
     753  58.3833
     754   8.2500
     755  58.3833
     7561
     757     0.000
     758   100.000
     7594000
     7601
     7611
     7621
     7631
     7641
     765rel_bj
     766
     76720041125 080000
     76820041126 080000
     769   8.2500
     770  58.3833
     771   8.2500
     772  58.3833
     7731
     774     0.000
     775   100.000
     7764000
     7771
     7781
     7791
     7801
     7811
     782rel_bj
     783
     78420041208 080000
     78520041209 080000
     786   8.2500
     787  58.3833
     788   8.2500
     789  58.3833
     7901
     791     0.000
     792   100.000
     7934000
     7941
     7951
     7961
     7971
     7981
     799rel_bj
     800
     80120041209 080000
     80220041210 080000
     803   8.2500
     804  58.3833
     805   8.2500
     806  58.3833
     8071
     808     0.000
     809   100.000
     8104000
     8111
     8121
     8131
     8141
     8151
     816rel_bj
     817
     81820041216 080000
     81920041217 080000
     820   8.2500
     821  58.3833
     822   8.2500
     823  58.3833
     8241
     825     0.000
     826   100.000
     8274000
     8281
     8291
     8301
     8311
     8321
     833rel_bj
     834
     83520041223 080000
     83620041224 080000
     837   8.2500
     838  58.3833
     839   8.2500
     840  58.3833
     8411
     842     0.000
     843   100.000
     8444000
     8451
     8461
     8471
     8481
     8491
     850rel_bj
     851
     85220041230 080000
     85320041231 080000
     854   8.2500
     855  58.3833
     856   8.2500
     857  58.3833
     8581
     859     0.000
     860   100.000
     8614000
     8621
     8631
     8641
     8651
     8661
     867rel_bj
     868

  • /trunk/options/SPECIES/SPECIES_001

    r30 r20  
    77TRACER             Tracer name
    88-999.9             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)
     9-9.9E-09           Wet deposition - A
     10                   Wet deposition - B
    1511-9.9               Dry deposition (gases) - D
    1612                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_002

    r30 r20  
    77O3                 Tracer name
    88-999.9             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)
     9-9.9E-09           Wet deposition - A
     10                   Wet deposition - B
    1511 1.5               Dry deposition (gases) - D
    16121.0e-02            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_003

    r30 r20  
    77NO                 Tracer name
    88-999.9             Species half life
    9  8.0E-06           Below cloud scavenging - A
    10 0.62               Below cloud scavenging - B
    11 2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
    12 0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
    13 0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
    14 1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
     9 8.0E-06           Wet deposition - A
     100.62               Wet deposition - B
    1511 1.2               Dry deposition (gases) - D
    16123.0e-03            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_004

    r30 r20  
    77NO2                Tracer name
    88-999.9             Species half life
    9  1.0E-05           Below cloud scavenging - A
    10 0.62               Below cloud scavenging - B
    11 2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
    12 0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
    13 0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
    14 1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
     9 1.0E-05           Wet deposition - A
     100.62               Wet deposition - B
    1511 1.6               Dry deposition (gases) - D
    16121.0e-02            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_005

    r30 r20  
    77HNO3               Tracer name
    88-999.9             Species half life
    9  5.0E-05           Below cloud scavenging - A
    10 0.62               Below cloud scavenging - B
    11 2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
    12 0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
    13 0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
    14 1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
     9 5.0E-05           Wet deposition - A
     100.62               Wet deposition - B
    1511 1.9               Dry deposition (gases) - D
    16121.0e+14            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_006

    r30 r20  
    77HNO2               Tracer name
    88-999.9             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)
     9-9.9E-09           Wet deposition - A
     10                   Wet deposition - B
    1511 1.6               Dry deposition (gases) - D
    16121.0e+05            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_007

    r30 r20  
    77H2O2               Tracer name
    88-999.9             Species half life
    9  1.0E-04           Below cloud scavenging - A
    10 0.62               Below cloud scavenging - B
    11 2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
    12 0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
    13 0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
    14 1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
     9 1.0E-04           Wet deposition - A
     100.62               Wet deposition - B
    1511 1.4               Dry deposition (gases) - D
    16121.0e+05            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_008

    r30 r20  
    77SO2                Tracer name
    88-999.9             Species half life
    9 -9.9E-09           Below cloud scavenging - A
    10 0.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)
     9-9.9E-09           Wet deposition - A
     100.62               Wet deposition - B
    1511 2.0               Dry deposition (gases) - D
    16121.0e+05            Dry deposition (gases) - Henrys const.
    … …  
    433916       1.463       1.162
    444017       1.426       1.141
    45 18       1.326       1.117
     4118       1.326       1.117 
    464219       1.225       1.109
    474320       1.082       1.095

  • /trunk/options/SPECIES/SPECIES_009

    r30 r20  
    77HCHO               Tracer name
    88-999.9             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)
     9-9.9E-09           Wet deposition - A
     10                   Wet deposition - B
    1511 1.3               Dry deposition (gases) - D
    16126.0e+03            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_010

    r30 r20  
    77PAN                Tracer name
    88-999.9             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)
     9-9.9E-09           Wet deposition - A
     10                   Wet deposition - B
    1511 2.6               Dry deposition (gases) - D
    16123.6e+00            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_011

    r30 r20  
    77NH3                Tracer name
    88-999.9             Species half life
    9  9.9E-05           Below cloud scavenging - A
    10 0.62               Below cloud scavenging - B
    11 2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
    12 0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
    13 0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
    14 1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
     9 9.9E-05           Wet deposition - A
     100.62               Wet deposition - B
    1511 1.1               Dry deposition (gases) - D
    16122.0e+14            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_012

    r30 r20  
    77SO4-aero           Tracer name
    88-999.9             Species half life
    9  5.0E-06           Below cloud scavenging - A
    10 0.62               Below cloud scavenging - B
    11 2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
    12 0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
    13 0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
    14 1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
     9 5.0E-06           Wet deposition - A
     100.62               Wet deposition - B
    1511-9.9               Dry deposition (gases) - D
    1612 1.0E-09           Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_013

    r30 r20  
    77NO3-aero           Tracer name
    88-999.9             Species half life
    9  5.0E-06           Below cloud scavenging - A
    10 0.62               Below cloud scavenging - B
    11 2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
    12 0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
    13 0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
    14 1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
     9 5.0E-06           Wet deposition - A
     100.62               Wet deposition - B
    1511-9.9               Dry deposition (gases) - D
    1612                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_014

    r30 r20  
    77I2-131             Tracer name
    88691200.0           Species half life
    9  8.0E-05           Below cloud scavenging - A
    10 0.62               Below cloud scavenging - B
    11 2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
    12 0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
    13 0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
    14 1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
     9 8.0E-05           Wet deposition - A
     100.62               Wet deposition - B
    1511 2.7               Dry deposition (gases) - D
    16121.0e+05            Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_015

    r30 r20  
    77I-131              Tracer name
    88691200.0           Species half life
    9  1.0E-04           Below cloud scavenging - A
    10 0.80               Below cloud scavenging - B
    11 2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
    12 0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
    13 0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
    14 1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
     9 1.0E-04           Wet deposition - A
     100.80               Wet deposition - B
    1511-9.9               Dry deposition (gases) - D
    1612                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_016

    r30 r20  
    77Cs-137             Tracer name
    88-999.9             Species half life
    9  1.0E-04           Below cloud scavenging - A
    10 0.80               Below cloud scavenging - B
    11 2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
    12 0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
    13 0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
    14 1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
     9 1.0E-04           Wet deposition - A
     100.80               Wet deposition - B
    1511-9.9               Dry deposition (gases) - D
    1612                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_017

    r30 r20  
    77Y-91               Tracer name
    885037120.0          Species half life
    9  1.0E-04           Below cloud scavenging - A
    10 0.80               Below cloud scavenging - B
    11 2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
    12 0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
    13 0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
    14 1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
     9 1.0E-04           Wet deposition - A
     100.80               Wet deposition - B
    1511-9.9               Dry deposition (gases) - D
    1612                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_018

    r30 r20  
    77Ru-106             Tracer name
    8831536000.0         Species half life
    9  1.0E-04           Below cloud scavenging - A
    10 0.80               Below cloud scavenging - B
    11 2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
    12 0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
    13 0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
    14 1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
     9 1.0E-04           Wet deposition - A
     100.80               Wet deposition - B
    1511-9.9               Dry deposition (gases) - D
    1612                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_019

    r30 r20  
    77Kr-85              Tracer name
    88-999.9             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)
     9-9.9E-09           Wet deposition - A
     10                   Wet deposition - B
    1511-9.9               Dry deposition (gases) - D
    1612                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_020

    r30 r20  
    77Sr-90              Tracer name
    88-999.9             Species half life
    9  1.0E-04           Below cloud scavenging - A
    10 0.80               Below cloud scavenging - B
    11 2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
    12 0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
    13 0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
    14 1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
     9 1.0E-04           Wet deposition - A
     100.80               Wet deposition - B
    1511-9.9               Dry deposition (gases) - D
    1612                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_021

    r30 r20  
    66****************************************************************************
    77Xe-133             Tracer name
    8 453168.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)
     8198720.0           Species half life
     9-9.9E-09           Wet deposition - A
     10                   Wet deposition - B
    1511-9.9               Dry deposition (gases) - D
    1612                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_022

    r30 r20  
    77CO                 Tracer name
    88-999.9             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)
     9-9.9E-09           Wet deposition - A
     10                   Wet deposition - B
    1511-9.9               Dry deposition (gases) - D
    1612                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_024

    r30 r20  
    77AIRTRACER             Tracer name
    88-999.9             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)
     9-9.9E-09           Wet deposition - A
     10                   Wet deposition - B
    1511-9.9               Dry deposition (gases) - D
    1612                   Dry deposition (gases) - Henrys const.

  • /trunk/options/SPECIES/SPECIES_025

    r30 r20  
    992.0E-05            Wet deposition - A
    10100.8                Wet deposition - B
    11 2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
    12 0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
    13 0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
    14 1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1511-9.9               Dry deposition (gases) - D
    1612                   Dry deposition (gases) - Henrys const.

  • /trunk/src/FLEXPART.f90

    r30 r20  
    4545  use conv_mod
    4646
    47 #ifdef NETCDF_OUTPUT
    48   use netcdf_output_mod, only: writeheader_netcdf
    49 #endif
    50 
    51 
    5247  implicit none
    5348
    … …  
    6560  call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
    6661
    67   ! FLEXPART version string
    68   flexversion='Version 9.2 beta (2014-07-01)'
    69   verbosity=0
    70 
     62  !
     63  flexversion='Version 9.1.8  (2013-12-08)'
     64  !verbosity=0
    7165  ! Read the pathnames where input/output files are stored
    7266  !*******************************************************
    … …  
    8276    call getarg(1,arg1)
    8377    pathfile=arg1
     78    verbosity=0
    8479    if (arg1(1:1).eq.'-') then
    85       write(pathfile,'(a11)') './pathnames'
    86       inline_options=arg1
     80        write(pathfile,'(a11)') './pathnames'
     81        inline_options=arg1
    8782    endif
    8883  case (0)
    8984    write(pathfile,'(a11)') './pathnames'
     85    verbosity=0
    9086  end select
    9187 
     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
    92111  ! Print the GPL License statement
    93112  !*******************************************************
    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'
     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'
    120119  endif
    121120  call readpaths(pathfile)
     121   
    122122 
    123123  if (verbosity.gt.1) then !show clock info
    … …  
    131131  endif
    132132
     133
    133134  ! Read the user specifications for the current model run
    134135  !*******************************************************
    135136
    136137  if (verbosity.gt.0) then
    137     write(*,*) 'call readcommand'
     138        WRITE(*,*) 'call readcommand'
    138139  endif
    139140  call readcommand
    140141  if (verbosity.gt.0) then
    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
     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
    149151
    150152  ! Read the age classes to be used
    151153  !********************************
    152154  if (verbosity.gt.0) then
    153     write(*,*) 'call readageclasses'
     155        WRITE(*,*) 'call readageclasses'
    154156  endif
    155157  call readageclasses
    156158
    157159  if (verbosity.gt.1) then   
    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
     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
    160162  endif     
     163 
     164
    161165
    162166  ! Read, which wind fields are available within the modelling period
    … …  
    164168
    165169  if (verbosity.gt.0) then
    166     write(*,*) 'call readavailable'
     170        WRITE(*,*) 'call readavailable'
    167171  endif 
    168172  call readavailable
    … …  
    173177 
    174178  if (verbosity.gt.0) then
    175      write(*,*) 'call gridcheck'
     179     WRITE(*,*) 'call gridcheck'
    176180  endif
    177181
    … …  
    179183
    180184  if (verbosity.gt.1) then   
    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
     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
    183187  endif     
    184 
    185   if (verbosity.gt.0) then
    186     write(*,*) 'call gridcheck_nests'
     188 
     189
     190  if (verbosity.gt.0) then
     191        WRITE(*,*) 'call gridcheck_nests'
    187192  endif 
    188193  call gridcheck_nests
    189194
     195
    190196  ! Read the output grid specifications
    191197  !************************************
    192198
    193199  if (verbosity.gt.0) then
    194     write(*,*) 'call readoutgrid'
     200        WRITE(*,*) 'call readoutgrid'
    195201  endif
    196202
    … …  
    198204
    199205  if (nested_output.eq.1) then
    200     call readoutgrid_nest
     206          call readoutgrid_nest
    201207    if (verbosity.gt.0) then
    202       write(*,*) '# readoutgrid_nest'
     208        WRITE(*,*) '# readoutgrid_nest'
    203209    endif
    204210  endif
    … …  
    226232  call readlanduse
    227233
     234
    228235  ! Assign fractional cover of landuse classes to each ECMWF grid point
    229236  !********************************************************************
    … …  
    234241  call assignland
    235242
     243
     244
    236245  ! Read the coordinates of the release locations
    237246  !**********************************************
    … …  
    242251  call readreleases
    243252
     253
    244254  ! Read and compute surface resistances to dry deposition of gases
    245255  !****************************************************************
    … …  
    257267    print*,'call coordtrafo'
    258268  endif
     269
    259270
    260271  ! Initialize all particles to non-existent
    … …  
    284295  endif
    285296
     297
    286298  ! Calculate volume, surface area, etc., of all output grid cells
    287299  ! Allocate fluxes and OHfield if necessary
    288300  !***************************************************************
    289301
     302
    290303  if (verbosity.gt.0) then
    291304    print*,'call outgrid_init'
    … …  
    294307  if (nested_output.eq.1) call outgrid_init_nest
    295308
     309
    296310  ! Read the OH field
    297311  !******************
    … …  
    299313  if (OHREA.eqv..TRUE.) then
    300314    if (verbosity.gt.0) then
    301       print*,'call readOHfield'
     315       print*,'call readOHfield'
    302316    endif
    303     call readOHfield
     317       call readOHfield
    304318  endif
    305319
    … …  
    308322  !******************************************************************
    309323
    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
    327324
    328325  if (verbosity.gt.0) then
    … …  
    335332  !if (nested_output.eq.1) call writeheader_nest
    336333  if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
     334
    337335  if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
    338336  if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
    339337
    340   if (lnetcdfout.ne.1) then
    341     open(unitdates,file=path(2)(1:length(2))//'dates')
    342   end if
     338
     339
     340  !open(unitdates,file=path(2)(1:length(2))//'dates')
    343341
    344342  if (verbosity.gt.0) then
    … …  
    348346  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj
    349347
     348
    350349  ! Releases can only start and end at discrete times (multiples of lsynctime)
    351350  !***************************************************************************
    … …  
    355354  endif
    356355  do i=1,numpoint
    357     ireleasestart(i)=nint(real(ireleasestart(i))/real(lsynctime))*lsynctime
    358     ireleaseend(i)=nint(real(ireleaseend(i))/real(lsynctime))*lsynctime
     356    ireleasestart(i)=nint(real(ireleasestart(i))/ &
     357         real(lsynctime))*lsynctime
     358    ireleaseend(i)=nint(real(ireleaseend(i))/ &
     359         real(lsynctime))*lsynctime
    359360  end do
     361
    360362
    361363  ! Initialize cloud-base mass fluxes for the convection scheme
    … …  
    379381  end do
    380382
     383
    381384  ! Calculate particle trajectories
    382385  !********************************
    … …  
    384387  if (verbosity.gt.0) then
    385388     if (verbosity.gt.1) then   
    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
     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
    388391     endif
    389392     if (info_flag.eq.1) then
    390        print*, 'Info only mode (stop)'   
    391        stop
     393         print*, 'info only mode (stop)'   
     394         stop
    392395     endif
    393396     print*,'call timemanager'
    … …  
    396399  call timemanager
    397400
    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
     401
     402  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
     403       &XPART MODEL RUN!'
    406404
    407405end program flexpart

  • /trunk/src/advance.f90

    r30 r20  
    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
    469465
    470466      if (zt.gt.h) then
    … …  
    703699  endif
    704700
    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
    715701
    716702  ! Check position: If trajectory outside model domain, terminate it
    … …  
    718704
    719705  if ((xt.lt.0.).or.(xt.ge.real(nxmin1)).or.(yt.lt.0.).or. &
    720        (yt.gt.real(nymin1))) then
     706       (yt.ge.real(nymin1))) then
    721707    nstop=3
    722708    return
    … …  
    873859  endif
    874860
    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 
    886861  ! Check position: If trajectory outside model domain, terminate it
    887862  !*****************************************************************
    888863
    889864  if ((xt.lt.0.).or.(xt.ge.real(nxmin1)).or.(yt.lt.0.).or. &
    890        (yt.gt.real(nymin1))) then
     865       (yt.ge.real(nymin1))) then
    891866    nstop=3
    892867    return

  • /trunk/src/calcpv.f90

    r30 r20  
    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
    52   real :: pvavr,ppml(0:nxmax-1,0:nymax-1,nuvzmax),ppmk(0:nxmax-1,0:nymax-1,nuvzmax)
     51  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
     52  real :: pvavr,ppml(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 
    7364  do jy=0,nymin1
    7465    if (sglobal.and.jy.eq.0) goto 10
    … …  
    115106      end if
    116107      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
    117112  !
    118113  ! Loop over the vertical
    … …  
    120115
    121116      do kl=1,nuvz
    122         theta=tth(ix,jy,kl,n)*ppmk(ix,jy,kl)
     117        ppmk=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
     118        theta=tth(ix,jy,kl,n)*(100000./ppmk)**kappa
    123119        klvrp=kl+1
    124120        klvrm=kl-1
    … …  
    129125        if (klvrp.gt.nuvz) klvrp=nuvz
    130126        if (klvrm.lt.1) klvrm=1
    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))
     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))
    134132
    135133  ! Compute vertical position at pot. temperature surface on subgrid
    … …  
    158156          kch=kch+1
    159157          k=kup
    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
     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
    173176            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           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
     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
    195200            endif
    196             vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
    197             goto 20
    198           endif
    199           goto 40
     201            goto 40
    200202  ! This section used when no values were found
    20120321      continue
    … …  
    234236          kch=kch+1
    235237          k=kup
    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
     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
    247254            endif
    248             uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
    249             goto 50
    250           endif
    25125571        continue
    252256  ! Downward branch
    … …  
    255259          kch=kch+1
    256260          k=kdn
    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
     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
    268277            endif
    269             uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
    270             goto 50
    271           endif
    272           goto 70
     278            goto 70
    273279  ! This section used when no values were found
    27428051      continue
    275281  ! Must use uu at current level and lat. juy becomes smaller by 1
    276           uy(jj)=uuh(ix,jy,kl)
    277           juy=juy-1
     282        uy(jj)=uuh(ix,jy,kl)
     283        juy=juy-1
    278284  ! Otherwise OK
    27928550        continue

  • /trunk/src/calcpv_nests.f90

    r30 r20  
    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
    51   real :: ppml(0:nxmaxn-1,0:nymaxn-1,nuvzmax),ppmk(0:nxmaxn-1,0:nymaxn-1,nuvzmax)
     50  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
     51  real :: ppml(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 
    7263  do jy=0,nyn(l)-1
    7364    phi = (ylat0n(l) + jy * dyn(l)) * pi / 180.
    … …  
    9788      if (ix.eq.0.or.ix.eq.nxn(l)-1) jumpx=1
    9889      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
    9994  !
    10095  ! Loop over the vertical
    … …  
    10297
    10398      do kl=1,nuvz
    104         theta=tthn(ix,jy,kl,n,l)*ppmk(ix,jy,kl)
     99        ppmk=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
     100        theta=tthn(ix,jy,kl,n,l)*(100000./ppmk)**kappa
    105101        klvrp=kl+1
    106102        klvrm=kl-1
    … …  
    111107        if (klvrp.gt.nuvz) klvrp=nuvz
    112108        if (klvrm.lt.1) klvrm=1
    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))
     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))
    116114
    117115  ! Compute vertical position at pot. temperature surface on subgrid
    … …  
    136134          kch=kch+1
    137135          k=kup
    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)
     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
    140140
    141141      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    … …  
    158158          kch=kch+1
    159159          k=kdn
    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)
     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
    162164      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    163165           ((thdn.le.theta).and.(thup.ge.theta))) then
    … …  
    210212          kch=kch+1
    211213          k=kup
    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)
     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
    214218      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    215219           ((thdn.le.theta).and.(thup.ge.theta))) then
    … …  
    231235          kch=kch+1
    232236          k=kdn
    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)
     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
    235241      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    236242           ((thdn.le.theta).and.(thup.ge.theta))) then

  • /trunk/src/com_mod.f90

    r30 r20  
    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)
    117115
    118116  integer :: nageclass,lage(maxageclass)
    … …  
    342340  real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2)
    343341  real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2)
     342  !scavenging NIK, PS
    344343  integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2)
    345344  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)
    346347
    347348
    … …  
    358359  !      rainout  conv/lsp dominated  2/3
    359360  !      washout  conv/lsp dominated  4/5
     361! PS 2013
     362!c icloudbot (m)        cloud bottom height
     363!c icloudthck (m)       cloud thickness     
    360364
    361365  ! pplev for the GFS version
    … …  
    684688
    685689  !********************
    686   ! Verbosity, testing flags, namelist I/O
     690  ! Verbosity, testing flags
    687691  !********************   
    688692  integer :: verbosity=0
    689693  integer :: info_flag=0
    690   integer :: count_clock, count_clock0,  count_rate, count_max
    691   real    :: tins
    692   logical :: nmlout=.true.
     694  INTEGER :: count_clock, count_clock0,  count_rate, count_max
    693695   
    694696

  • /trunk/src/gridcheck.f90

    r30 r20  
    8484  real(kind=4) :: xaux1,xaux2,yaux1,yaux2
    8585  real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
    86   integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parID
     86  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
    8787  !HSO  end
    8888  integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
    89   real :: sizesouth,sizenorth,xauxa,pint,conversion_factor
     89  real :: sizesouth,sizenorth,xauxa,pint
    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)
     100  !real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
     101  real(kind=4) :: zsec2(184),zsec4(jpunp)
    101102  character(len=1) :: opt
    102103
    … …  
    171172  call grib_check(iret,gribFunction,gribErrorMsg)
    172173  call grib_get_int(igrib,'level',valSurf,iret)
    173   call grib_check(iret,gribFunction,gribErrorMsg)
    174   call grib_get_int(igrib,'paramId',parId,iret)
    175174  call grib_check(iret,gribFunction,gribErrorMsg)
    176175
    … …  
    194193  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
    195194    isec1(6)=135         ! indicatorOfParameter
    196   elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
    197     isec1(6)=135         ! indicatorOfParameter
    198195  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    199196    isec1(6)=151         ! indicatorOfParameter
    … …  
    208205  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    209206    isec1(6)=141         ! indicatorOfParameter
    210   elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
     207  elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
    211208    isec1(6)=164         ! indicatorOfParameter
    212   elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
     209  elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
    213210    isec1(6)=142         ! indicatorOfParameter
    214211  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    … …  
    218215  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    219216    isec1(6)=176         ! indicatorOfParameter
    220   elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
     217  elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
    221218    isec1(6)=180         ! indicatorOfParameter
    222   elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
     219  elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
    223220    isec1(6)=181         ! indicatorOfParameter
    224221  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    225222    isec1(6)=129         ! indicatorOfParameter
    226   elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
     223  elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
    227224    isec1(6)=160         ! indicatorOfParameter
    228225  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    … …  
    233230         parCat,parNum,typSurf
    234231  endif
    235   if(parId .ne. isec1(6) .and. parId .ne. 77) then
    236     write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
    237 !    stop
    238   endif
    239232
    240233  endif
    … …  
    272265
    273266  !HSO  get the second part of the grid dimensions only from GRiB1 messages
    274   if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
     267  if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
    275268    call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
    276269         xaux2in,iret)
    … …  
    298291    dyconst=180./(dy*r_earth*pi)
    299292    gotGrid=1
     293
    300294  ! Check whether fields are global
    301295  ! If they contain the poles, specify polar stereographic map
    … …  
    443437  !********************
    444438
    445   write(*,'(a,2i7)') ' Vertical levels in ECMWF data: ', &
     439  write(*,*)
     440  write(*,*)
     441  write(*,'(a,2i7)') '# of vertical levels in ECMWF data: ', &
    446442       nuvz+1,nwz
    447443  write(*,*)
    448   write(*,'(a)') ' Mother domain:'
    449   write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
     444  write(*,'(a)') 'Mother domain:'
     445  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
    450446       xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
    451   write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range : ', &
     447  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
    452448       ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
    453449  write(*,*)
    … …  
    471467    k=nlev_ec+1+numskip+i
    472468    akm(nwz-i+1)=zsec2(j)
    473   !   write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)
    474469    bkm(nwz-i+1)=zsec2(k)
    475470  end do
    … …  
    558553
    559554end subroutine gridcheck
    560 

  • /trunk/src/gridcheck_gfs.f90

    r30 r20  
    423423  write(*,*)
    424424  write(*,*)
    425   write(*,'(a,2i7)') 'Vertical levels in NCEP data: ', &
     425  write(*,'(a,2i7)') '# of 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

    r30 r20  
    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
    5150  integer :: gotGrib
    5251  !HSO  end
    … …  
    5655  real(kind=4) :: xaux1,xaux2,yaux1,yaux2
    5756  real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
    58   real :: conversion_factor !added by mc to make it consistent with new gridchek.f90
    5957
    6058  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
    … …  
    152150  call grib_get_int(igrib,'level',valSurf,iret)
    153151  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
    156152
    157153  !print*,discipl,parCat,parNum,typSurf,valSurf
    … …  
    174170  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
    175171    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    !                     ! " " " " " " " " " " " " " " " " " " " " " " " "  " " "
    178172  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    179173    isec1(6)=151         ! indicatorOfParameter
    … …  
    188182  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    189183    isec1(6)=141         ! indicatorOfParameter
    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
     184  elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
    191185    isec1(6)=164         ! indicatorOfParameter
    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
     186  elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
    193187    isec1(6)=142         ! indicatorOfParameter
    194188  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    … …  
    198192  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    199193    isec1(6)=176         ! indicatorOfParameter
    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
     194  elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
    201195    isec1(6)=180         ! indicatorOfParameter
    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
     196  elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
    203197    isec1(6)=181         ! indicatorOfParameter
    204198  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    205199    isec1(6)=129         ! indicatorOfParameter
    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
     200  elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
    207201    isec1(6)=160         ! indicatorOfParameter
    208202  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    … …  
    213207         parCat,parNum,typSurf
    214208  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
    219209
    220210  endif
    … …  
    247237
    248238  !HSO  get the second part of the grid dimensions only from GRiB1 messages
    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
     239  if ((gribVer.eq.1).and.(gotGrib.eq.0)) then
     240    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
    251241         xaux1in,iret)
    252242    call grib_check(iret,gribFunction,gribErrorMsg)
    … …  
    264254    yaux1=yaux1in
    265255    yaux2=yaux2in
    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
     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.
    271261    xlon0n(l)=xaux1
    272262    ylat0n(l)=yaux1
    273263    dxn(l)=(xaux2-xaux1)/real(nxn(l)-1)
    274264    dyn(l)=(yaux2-yaux1)/real(nyn(l)-1)
    275     gotGrib=1 !commetn by mc note tahthere gotGRIB is used instead of gotGrid!!!
     265    gotGrib=1
    276266  endif ! ifield.eq.1
    277267
    … …  
    350340  !********************
    351341
    352   write(*,'(a,i2,a)') ' Nested domain ',l,':'
    353   write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
     342  write(*,'(a,i2)') 'Nested domain #: ',l
     343  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
    354344       xlon0n(l),' to ',xlon0n(l)+(nxn(l)-1)*dxn(l), &
    355345       '   Grid distance: ',dxn(l)
    356   write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range : ', &
     346  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
    357347       ylat0n(l),' to ',ylat0n(l)+(nyn(l)-1)*dyn(l), &
    358348       '   Grid distance: ',dyn(l)

  • /trunk/src/interpol_rain.f90

    r30 r20  
    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
     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
    2633  !****************************************************************************
    2734  !                                                                           *
    … …  
    4047  !                                                                           *
    4148  !     30 August 1996                                                        *
    42   !                                                                           *
     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  *
    4353  !****************************************************************************
    4454  !                                                                           *
    … …  
    7080
    7181  integer :: nx,ny,nxmax,nymax,nzmax,memind(2),m,ix,jy,ixp,jyp
    72   integer :: itime,itime1,itime2,level,indexh
     82  !integer :: itime,itime1,itime2,level,indexh
     83  integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4
     84  integer :: intiy1,intiy2,ipsum,icmv 
    7385  real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2)
    7486  real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2)
    7587  real :: yy3(0:nxmax-1,0:nymax-1,nzmax,2)
    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
     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
    7893
    7994
    … …  
    127142         + p3*yy3(ix ,jyp,level,indexh) &
    128143         + 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
    129183  end do
    130184
    … …  
    134188  !************************************
    135189
     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
     196
    136197  dt1=real(itime-itime1)
    137198  dt2=real(itime2-itime)
    … …  
    143204
    144205
     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
    145223end subroutine interpol_rain

  • /trunk/src/makefile

    r30 r20  
    11SHELL = /bin/bash
    2 MAIN = FP_ecmwf_gfortran
     2TARGET = local
     3WINDS=ecmwf
     4#WINDS=gfs
     5#WINDS=fnl
    36
    47FC       = gfortran
    5 INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
    6 LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
    7 LIBPATH2 =   /usr/lib/x86_64-linux-gnu/
    88FFLAGS   =   -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
    9 #FFLAGS   =   -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
     9
     10ifeq ($(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_
     16endif
     17ifeq ($(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_
     25endif
     26
    1027LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
    1128
    … …  
    2037OBJECTS = \
    2138writeheader.o  writeheader_txt.o   writeheader_surf.o       assignland.o\
    22 calcpar.o               part0.o \
     39               part0.o \
    2340caldate.o               partdep.o \
    2441coordtrafo.o            psih.o \
    … …  
    3653getrc.o                 readreleases.o \
    3754getvdep.o               readspecies.o \
    38 interpol_misslev.o      readwind.o \
    39 conccalc.o              richardson.o \
     55interpol_misslev.o      \
     56conccalc.o              \
    4057concoutput.o  concoutput_surf.o          scalev.o \
    4158pbl_profile.o           readOHfield.o\
    4259juldate.o               timemanager.o \
    4360interpol_vdep.o         interpol_rain.o \
    44 verttransform.o         partoutput.o \
     61partoutput.o \
    4562hanna.o                 wetdepokernel.o \
    4663mean.o                  wetdepo.o \
    4764hanna_short.o           windalign.o \
    48 obukhov.o               gridcheck.o \
    4965hanna1.o                initialize.o \
    50                         gridcheck_nests.o \
    51 readwind_nests.o        calcpar_nests.o \
     66                           calcpar_nests.o \
    5267verttransform_nests.o   interpol_all_nests.o \
    5368interpol_wind_nests.o   interpol_misslev_nests.o \
    5469interpol_vdep_nests.o   interpol_rain_nests.o \
    55 getvdep_nests.o \
     70getvdep_nests.o   gridcheck_nests.o \
     71readwind_nests.o  \
    5672readageclasses.o        readpartpositions.o \
    5773calcfluxes.o            fluxoutput.o \
    5874qvsat.o                 skplin.o \
    59 convmix.o               calcmatrix.o \
    6075convect43c.o               redist.o \
    6176sort2.o                 distance.o \
    … …  
    7691
    7792
    78 $(MAIN): $(MODOBJS) $(OBJECTS)
    79         $(FC) *.o -o $(MAIN) $(LDFLAGS)
     93ifeq ($(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
     99endif
     100
     101ifeq ($(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
     107endif
     108
     109$(MAIN): $(MODOBJS) $(OBJECTS)  $(OBJECTS_WINDS)
     110        $(FC) *.o -o $(MAIN)_$(WINDS) $(LDFLAGS)
    80111
    81112$(OBJECTS): $(MODOBJS)
    … …  
    87118        rm *.o *.mod
    88119
    89 cleanall:
    90         rm *.o *.mod $(MAIN)

  • /trunk/src/outgrid_init.f90

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

  • /trunk/src/par_mod.f90

    r30 r20  
    121121  !*********************************************
    122122 
    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
     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
    126125  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26
    127126  !integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
    … …  
    155154
    156155  !integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
    157   integer,parameter :: maxnests=1,nxmaxn=351,nymaxn=351 !ECMWF
    158   !integer,parameter :: maxnests=1, nxmaxn=201, nymaxn=161 ! FNL XF
     156  !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151 !ECMWF
     157  integer,parameter :: maxnests=1, nxmaxn=201, nymaxn=161 ! FNL XF
    159158  ! maxnests                maximum number of nested grids
    160159  ! nxmaxn,nymaxn           maximum dimension of nested wind fields in
    … …  
    198197  !**************************************************
    199198
    200   integer,parameter :: maxpart=150000
    201   integer,parameter :: maxspec=4
     199  !integer,parameter :: maxpart=4000000
     200  integer,parameter :: maxpart=2000
     201  integer,parameter :: maxspec=6
    202202
    203203

  • /trunk/src/readageclasses.f90

    r30 r20  
    2828  !                                                                            *
    2929  !     Author: A. Stohl                                                       *
     30  !                                                                            *
    3031  !     20 March 2000                                                          *
    31   !     HSO, 1 July 2014                                                       *
    32   !     Added optional namelist input                                          *
    3332  !                                                                            *
    3433  !*****************************************************************************
    … …  
    4746  integer :: i
    4847
    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
    5848
    5949  ! If age spectra calculation is switched off, set number of age classes
    … …  
    6757  endif
    6858
     59
    6960  ! If age spectra claculation is switched on,
    7061  ! open the AGECLASSSES file and read user options
    7162  !************************************************
    7263
    73   open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES',form='formatted',status='old',err=999)
     64  open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES', &
     65       status='old',err=999)
    7466
    75   ! try to read in as a namelist
    76   read(unitageclasses,ageclass,iostat=readerror)
    77   close(unitageclasses)
     67  do i=1,13
     68    read(unitageclasses,*)
     69  end do
     70  read(unitageclasses,*) nageclass
    7871
    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
    9872
    9973  if (nageclass.gt.maxageclass) then
    … …  
    10680  endif
    10781
     82  read(unitageclasses,*) lage(1)
    10883  if (lage(1).le.0) then
    10984    write(*,*) ' #### FLEXPART MODEL ERROR! AGE OF FIRST      #### '
    … …  
    11489
    11590  do i=2,nageclass
     91    read(unitageclasses,*) lage(i)
    11692    if (lage(i).le.lage(i-1)) then
    11793      write(*,*) ' #### FLEXPART MODEL ERROR! AGE CLASSES     #### '
    … …  
    129105  stop
    130106
    131 1000  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 
    137107end subroutine readageclasses

  • /trunk/src/readavailable.f90

    r30 r20  
    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

    r30 r20  
    2929  !                                                                            *
    3030  !     18 May 1996                                                            *
    31   !     HSO, 1 July 2014                                                       *
    32   !     Added optional namelist input                                          *
    3331  !                                                                            *
    3432  !*****************************************************************************
    … …  
    8280  character(len=50) :: line
    8381  logical :: old
     82  logical :: nmlout=.true. !.false.
    8483  integer :: readerror
    8584
    8685  namelist /command/ &
    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   
     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   
    113111
    114112  ! Presetting namelist command
    115   ldirect=0
     113  ldirect=1
    116114  ibdate=20000101
    117115  ibtime=0
    … …  
    139137  nested_output=0
    140138  linit_cond=0
    141   lnetcdfout=0
    142139  surf_only=0
    143140
    … …  
    145142  ! Namelist input first: try to read as namelist file
    146143  !**************************************************************************
    147   open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old',form='formatted',err=999)
    148 
    149   ! try namelist input (default)
     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
    150156  read(unitcommand,command,iostat=readerror)
    151157  close(unitcommand)
    152158
    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)
     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)
    157164
    158165    ! Check the format of the COMMAND file (either in free format,
    … …  
    166173    if (index(line,'LDIRECT') .eq. 0) then
    167174      old = .false.
    168       write(*,*) 'COMMAND in old short format, please update to namelist format'
    169175    else
    170176      old = .true.
    171       write(*,*) 'COMMAND in old long format, please update to namelist format'
    172177    endif
    173178    rewind(unitcommand)
    174179
    175 
    176180    ! Read parameters
    177181    !****************
    … …  
    179183    call skplin(7,unitcommand)
    180184    if (old) call skplin(1,unitcommand)
     185
    181186    read(unitcommand,*) ldirect
    182187    if (old) call skplin(3,unitcommand)
    … …  
    226231    if (old) call skplin(3,unitcommand)
    227232    read(unitcommand,*) linit_cond
    228     if (old) call skplin(3,unitcommand)
    229     read(unitcommand,*) surf_only
    230233    close(unitcommand)
    231234
    … …  
    234237  ! write command file in namelist format to output directory if requested
    235238  if (nmlout.eqv..true.) then
     239    !open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist.out',status='new',err=1000)
    236240    open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist',err=1000)
    237241    write(unitcommand,nml=command)
    238242    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)
    239246  endif
    240247
    … …  
    346353  endif
    347354
    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 
    359355  ! Check whether a valid option for gridded model output has been chosen
    360356  !**********************************************************************
    … …  
    362358  if ((iout.lt.1).or.(iout.gt.5)) then
    363359    write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
    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                       #### '
     360    write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4, OR 5!          #### '
    367361    stop
    368362  endif
    … …  
    376370    stop
    377371  endif
     372
    378373
    379374
    … …  
    597592
    5985931000   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "COMMAND"    #### '
    599   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    600   write(*,'(a)') path(2)(1:length(1))
    601   stop
     594       write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     595        write(*,'(a)') path(2)(1:length(1))
     596        stop
    602597end subroutine readcommand

  • /trunk/src/readoutgrid.f90

    r30 r20  
    2929  !                                                                            *
    3030  !     4 June 1996                                                            *
    31   !     HSO, 1 July 2014
    32   !     Added optional namelist input
    3331  !                                                                            *
    3432  !*****************************************************************************
    … …  
    5553  real,parameter :: eps=1.e-4
    5654
    57   ! namelist variables
    58   integer, parameter :: maxoutlev=500
    59   integer :: readerror
    60   real,allocatable, dimension (:) :: outheights
    6155
    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
    7656
    7757  ! Open the OUTGRID file and read output grid specifications
    7858  !**********************************************************
    7959
    80   open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old',form='formatted',err=999)
     60  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old', &
     61       err=999)
    8162
    82   ! try namelist input
    83   read(unitoutgrid,outgrid,iostat=readerror)
    84   close(unitoutgrid)
    8563
    86   if ((dxout.le.0).or.(readerror.ne.0)) then
     64  call skplin(5,unitoutgrid)
    8765
    88     readerror=1
    8966
    90     open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old',err=999)
     67  ! 1.  Read horizontal grid specifications
     68  !****************************************
    9169
    92     call skplin(5,unitoutgrid)
     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
    9382
    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
    11183
    11284  ! Check validity of output grid (shall be within model domain)
    … …  
    11991  if ((outlon0+eps.lt.xlon0).or.(outlat0+eps.lt.ylat0) &
    12092       .or.(xr.gt.xr1+eps).or.(yr.gt.yr1+eps)) then
     93    write(*,*) 'outlon0,outlat0:'
    12194    write(*,*) outlon0,outlat0
     95    write(*,*) 'xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout:'
    12296    write(*,*) xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout
    12397    write(*,*) ' #### FLEXPART MODEL ERROR! PART OF OUTPUT    ####'
    … …  
    130104  ! 2. Count Vertical levels of output grid
    131105  !****************************************
    132 
    133   if (readerror.ne.0) then
    134     j=0
    135 100 j=j+1
     106  j=0
     107100   j=j+1
    136108    do i=1,3
    137109      read(unitoutgrid,*,end=99)
    … …  
    140112    if (outhelp.eq.0.) goto 99
    141113    goto 100
    142 99  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
     11499   numzgrid=j-1
    149115
    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'
     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)
    154126
    155127  ! 2. Vertical levels of output grid
    156128  !**********************************
    157129
    158   if (readerror.ne.0) then
     130  j=0
     1311000   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
     139990   numzgrid=j-1
    159140
    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
    196141
    197142  ! Check whether vertical levels are specified in ascending order
    … …  
    215160  end do
    216161
     162
    217163  xoutshift=xlon0-outlon0
    218164  youtshift=ylat0-outlat0
     165  close(unitoutgrid)
    219166
    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'
     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'
    230182  return
    231183
    232184
    233 999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
     185999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
    234186  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    235   write(*,'(a)') path(1)(1:length(1))
    236   stop
    237 
    238 1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
    239   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    240   write(*,'(a)') path(2)(1:length(2))
     187  write(*,*) ' #### xxx/flexpart/options                    #### '
    241188  stop
    242189

  • /trunk/src/readoutgrid_nest.f90

    r30 r20  
    5353  real,parameter :: eps=1.e-4
    5454
    55   integer :: readerror
    5655
    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
    6556
    6657  ! Open the OUTGRID file and read output grid specifications
    6758  !**********************************************************
    6859
    69   open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',form='formatted',status='old',err=999)
     60  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST', &
     61       status='old',err=999)
    7062
    71   ! try namelist input
    72   read(unitoutgrid,outgridn,iostat=readerror)
    73   close(unitoutgrid)
    7463
    75   if ((dxoutn.le.0).or.(readerror.ne.0)) then
     64  call skplin(5,unitoutgrid)
    7665
    77     open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',status='old',err=999)
    78     call skplin(5,unitoutgrid)
    7966
    80     ! 1.  Read horizontal grid specifications
    81     !****************************************
     67  ! 1.  Read horizontal grid specifications
     68  !****************************************
    8269
    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
     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
    9582
    96     close(unitoutgrid)
    97   endif
    9883
    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'
     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'
    11293
    11394  ! Check validity of output grid (shall be within model domain)
    … …  
    129110  xoutshiftn=xlon0-outlon0n
    130111  youtshiftn=ylat0-outlat0n
     112  close(unitoutgrid)
    131113  return
    132114
    133 999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
     115
     116999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE OUTGRID_NEST #### '
    134117  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    135   write(*,'(a)') path(1)(1:length(1))
    136   stop
    137 
    138 1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
    139   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    140   write(*,'(a)') path(2)(1:length(2))
     118  write(*,*) ' #### xxx/flexpart/options                    #### '
    141119  stop
    142120

  • /trunk/src/readpaths.f90

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

  • /trunk/src/readreceptors.f90

    r30 r20  
    2727  !                                                                            *
    2828  !     Author: A. Stohl                                                       *
     29  !                                                                            *
    2930  !     1 August 1996                                                          *
    30   !     HSO, 14 August 2013
    31   !     Added optional namelist input
    3231  !                                                                            *
    3332  !*****************************************************************************
    … …  
    5251  character(len=16) :: receptor
    5352
    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
    6353
    6454  ! For backward runs, do not allow receptor output. Thus, set number of receptors to zero
    … …  
    7464  !************************************************************
    7565
    76   open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS',form='formatted',status='old',err=999)
     66  open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS', &
     67       status='old',err=999)
    7768
    78   ! try namelist input
    79   read(unitreceptor,receptors,iostat=readerror)
     69  call skplin(5,unitreceptor)
    8070
    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
    8571
    86   if ((lon.lt.-900).or.(readerror.ne.0)) then
     72  ! Read the names and coordinates of the receptors
     73  !************************************************
    8774
    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
    96 100 j=j+1
     75  j=0
     76100   j=j+1
    9777    read(unitreceptor,*,end=99)
    9878    read(unitreceptor,*,end=99)
    … …  
    10989    endif
    11090    if (j.gt.maxreceptor) then
    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   #### '
     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   #### '
    11595    endif
    11696    receptorname(j)=receptor
    … …  
    120100    ym=r_earth*dy/180.*pi
    121101    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 
    130102    goto 100
    131103
    132 99  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 
     10499   numreceptor=j-1
     105  close(unitreceptor)
    173106  return
    174107
    175108
    176 999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
     109999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
    177110  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    178111  write(*,'(a)') path(1)(1:length(1))
    179112  stop
    180113
    181 1000 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 
    186114end subroutine readreceptors

  • /trunk/src/readreleases.f90

    r30 r20  
    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
    3735  !                                                                            *
    3836  !*****************************************************************************
    … …  
    5755  ! xpoint2,ypoint2     geograf. coordinates of upper right corner of release  *
    5856  !                     area                                                   *
    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        *
     57  ! weta, wetb          parameters to determine the wet scavenging coefficient *
    6258  ! zpoint1,zpoint2     height range, over which release takes place           *
    6359  ! num_min_discrete    if less, release cannot be randomized and happens at   *
    … …  
    7369  implicit none
    7470
    75   integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat
     71  integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
    7672  integer,parameter :: num_min_discrete=100
    7773  real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
    … …  
    8076  logical :: old
    8177
    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
     78  !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
     90901   format (a)
     91  if (index(line,'Total') .eq. 0) then
     92    old = .false.
     93  else
     94    old = .true.
     95  endif
     96  rewind(unitreleases)
     97
     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)
     110  end do
    107111
    108112  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 
    120   !sec, read release to find how many releasepoints should be allocated
    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
    146 901 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
    168 100 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)
     113100   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
    174136    read(unitreleases,*,err=998) xdum
    175137    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 
    200 25  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 
    217   rewind(unitreleases)
    218 
    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.
    264138  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
     14525   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
    265193
    266194  !now save the information
    … …  
    273201  end do
    274202
    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
     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)
     218  do i=1,nspec
     219    read(unitreleases,*,err=998) specnum_rel
    286220    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 
    294   do i=1,nspec
    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
     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
    318239      write(*,*) 'For mixing ratio output, valid molar weight'
    319240      write(*,*) 'must be specified for all simulated species.'
    … …  
    323244    endif
    324245
    325     ! Radioactive decay
    326     !******************
     246
     247  ! Radioactive decay
     248  !******************
    327249
    328250    decay(i)=0.693147/decay(i) !conversion half life to decay constant
    329251
    330252
    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     !****************************
     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  !****************************
    338261
    339262    vsetaver(i)=0.
    340263    cunningham(i)=0.
    341264    dquer(i)=dquer(i)*1000000.         ! Conversion m to um
    342     if (density(i).gt.0.) then         ! Additional parameters
    343       call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
     265    if (density(i).gt.0.) then                  ! Additional parameters
     266     call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
    344267      do j=1,ni
    345268        fract(i,j)=fracth(j)
    … …  
    349272        vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
    350273      end do
    351       write(*,*) 'Average settling velocity: ',i,vsetaver(i)
    352     endif
    353 
    354     ! Dry deposition for constant deposition velocity
    355     !************************************************
     274      write(*,*) 'Average setting velocity: ',i,vsetaver(i)
     275    endif
     276
     277  ! Dry deposition for constant deposition velocity
     278  !************************************************
    356279
    357280    dryvel(i)=dryvel(i)*0.01         ! conversion to m/s
    358281
    359     ! Check if wet deposition or OH reaction shall be calculated
    360     !***********************************************************
     282  ! Check if wet deposition or OH reaction shall be calculated
     283  !***********************************************************
    361284    if (weta(i).gt.0.)  then
    362285      WETDEP=.true.
    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 
     286      write (*,*) 'Wetdeposition switched on: ',weta(i),i
     287    endif
    378288    if (ohreact(i).gt.0) then
    379289      OHREA=.true.
    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
     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
    384296      DRYDEP=.true.
    385297      DRYDEPSPEC(i)=.true.
    … …  
    388300  end do
    389301
    390   if (WETDEP.or.DRYDEP) DEP=.true.
     302    if (WETDEP.or.DRYDEP) DEP=.true.
    391303
    392304  ! Read specifications for each release point
    393305  !*******************************************
    394   numpoints=numpoint
     306
    395307  numpoint=0
    396308  numpartmax=0
    397309  releaserate=0.
    398 101   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 
     3101000   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)
    430339  else
    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
     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
    495352
    496353  ! If a release point contains no particles, stop and issue error message
    … …  
    563420  endif
    564421
     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
    565427  ! Determine the release rate (particles per second) and total number
    566428  ! of particles released during the simulation
    … …  
    574436  endif
    575437  numpartmax=numpartmax+npart(numpoint)
    576   goto 101
     438  goto 1000
     439
    577440
    578441250   close(unitreleases)
    579442
    580   if (nmlout.eqv..true.) then
    581     close(unitreleasesout)
    582   endif
    583 
    584   write (*,*) 'Particles allocated (maxpart)  : ',maxpart
    585   write (*,*) 'Particles released (numpartmax): ',numpartmax
     443  write (*,*) ' Particles allocated for this run: ',maxpart, ', released in simulation: ',  numpartmax
    586444  numpoint=numpoint-1
    587445
    … …  
    591449    maxpointspec_act=1
    592450  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   !************************************************************************
    597451
    598452  if (releaserate.gt. &
    … …  
    652506  stop
    653507
    654 1000 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 
    659508end subroutine readreleases

  • /trunk/src/readspecies.f90

    r30 r20  
    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
    3832  !                                                                            *
    3933  !*****************************************************************************
    … …  
    4236  ! decaytime(maxtable)  half time for radiological decay                      *
    4337  ! specname(maxtable)   names of chemical species, radionuclides              *
    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         *
     38  ! wetscava, wetscavb   Parameters for determining scavenging coefficient     *
    4939  ! ohreact              OH reaction rate                                      *
    5040  ! id_spec              SPECIES number as referenced in RELEASE file          *
    … …  
    6555  logical :: spec_found
    6656
    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 
    10057  ! Open the SPECIES file and read species names and properties
    10158  !************************************************************
    10259  specnum(pos_spec)=id_spec
    10360  write(aspecnumb,'(i3.3)') specnum(pos_spec)
    104   open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998)
     61  open(unitspecies,file= &
     62       path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old', &
     63       err=998)
    10564  !write(*,*) 'reading SPECIES',specnum(pos_spec)
    10665
    10766  ASSSPEC=.FALSE.
    10867
    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
     68  do i=1,6
     69    read(unitspecies,*)
     70  end do
    12271
    12372    read(unitspecies,'(a10)',end=22) species(pos_spec)
    … …  
    12978    read(unitspecies,'(f18.2)',end=22) wetb(pos_spec)
    13079  !  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 
    14280    read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
    14381  !  write(*,*) reldiff(pos_spec)
    … …  
    162100    read(unitspecies,'(f18.2)',end=22) kao(pos_spec)
    163101  !       write(*,*) kao(pos_spec)
     102    i=pos_spec
    164103
    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)
     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
     106    endif
    184107
    185   else
     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
    186121
    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
     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
    206124
    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
     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
    226131    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
    23913220   continue
    240133
    … …  
    242135  ! Read in daily and day-of-week variation of emissions, if available
    243136  !*******************************************************************
    244   ! HSO: This is not yet implemented as namelist parameters
    245   ! default values set to 1
    246137
    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
     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
    257146
    258147    read(unitspecies,*,end=22)
    … …  
    265154    end do
    266155
    267   endif
     15622   close(unitspecies)
    268157
    269 22 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
     158   return
    279159
    280160996   write(*,*) '#####################################################'
    … …  
    305185  stop
    306186
    307 1000 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
    311187
    312188end subroutine readspecies

  • /trunk/src/readwind.f90

    r30 r20  
    7171  integer :: iret
    7272  integer :: igrib
    73   integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
     73  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
    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
     92  real(kind=4) :: xaux,yaux,xaux0,yaux0
    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,conversion_factor
     96  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
    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
    103108
    104109  hflswitch=.false.
    … …  
    150155  endif
    151156
    152   conversion_factor=1.
    153 
    154157  else
    155158
    … …  
    165168  call grib_check(iret,gribFunction,gribErrorMsg)
    166169  call grib_get_int(igrib,'level',valSurf,iret)
    167   call grib_check(iret,gribFunction,gribErrorMsg)
    168   call grib_get_int(igrib,'paramId',parId,iret)
    169170  call grib_check(iret,gribFunction,gribErrorMsg)
    170171
    … …  
    176177  isec1(8)=-1
    177178  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
    192190    isec1(6)=135         ! indicatorOfParameter
    193191  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    … …  
    203201  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    204202    isec1(6)=141         ! indicatorOfParameter
    205     conversion_factor=1000.
    206   elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
     203  elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
    207204    isec1(6)=164         ! indicatorOfParameter
    208   elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
     205  elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
    209206    isec1(6)=142         ! indicatorOfParameter
    210207  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    211208    isec1(6)=143         ! indicatorOfParameter
    212     conversion_factor=1000.
    213209  elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
    214210    isec1(6)=146         ! indicatorOfParameter
    215211  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    216212    isec1(6)=176         ! indicatorOfParameter
    217   elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
     213  elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
    218214    isec1(6)=180         ! indicatorOfParameter
    219   elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
     215  elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
    220216    isec1(6)=181         ! indicatorOfParameter
    221217  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    222218    isec1(6)=129         ! indicatorOfParameter
    223   elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
     219  elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
    224220    isec1(6)=160         ! indicatorOfParameter
    225221  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    … …  
    227223    isec1(6)=172         ! indicatorOfParameter
    228224  else
    229     print*,'***WARNING: undefined GRiB2 message found!',discipl, &
     225    print*,'***ERROR: undefined GRiB2 message found!',discipl, &
    230226         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
    235227  endif
    236228
    … …  
    245237  !HSO  get the required fields from section 2 in a gribex compatible manner
    246238  if (ifield.eq.1) then
    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))
     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))
    252247  call grib_check(iret,gribFunction,gribErrorMsg)
    253248  ! CHECK GRID SPECIFICATIONS
    … …  
    255250  if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
    256251  if(isec2(12)/2-1.ne.nlev_ec) &
    257   stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
     252       stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
    258253  endif ! ifield
    259254
    260255  !HSO  get the second part of the grid dimensions only from GRiB1 messages
    261   if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
     256  if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
    262257    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
    263258         xauxin,iret)
    … …  
    266261         yauxin,iret)
    267262    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 
    271263    xaux=xauxin+real(nxshift)*dx
    272264    yaux=yauxin
    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'
     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'
    278275    gotGrid=1
    279276  endif ! gotGrid
    … …  
    301298           zsec4(nxfield*(ny-j-1)+i+1)
    302299      if(isec1(6).eq.141) sd(i,j,1,n)= &!! SNOW DEPTH
    303            zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
     300           zsec4(nxfield*(ny-j-1)+i+1)
    304301      if(isec1(6).eq.151) msl(i,j,1,n)= &!! SEA LEVEL PRESS.
    305302           zsec4(nxfield*(ny-j-1)+i+1)
    … …  
    319316      endif
    320317      if(isec1(6).eq.143) then                      !! CONVECTIVE PREC.
    321         convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
     318        convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
    322319        if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
    323320      endif
    … …  
    369366      do j=0,nymin1
    370367        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
    371384      end do
    372385    end do
    … …  
    466479
    467480end subroutine readwind
    468 

  • /trunk/src/readwind_nests.f90

    r30 r20  
    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
    5453  integer :: gotGrid
    5554  !HSO  end
    … …  
    7069  integer :: isec1(56),isec2(22+nxmaxn+nymaxn)
    7170  real(kind=4) :: zsec4(jpunp)
    72   real(kind=4) :: xaux,yaux
    7371  real(kind=8) :: xauxin,yauxin
    74   real,parameter :: eps=1.e-4
     72  real(kind=4) :: xaux,yaux,xaux0,yaux0
    7573  real :: ewss(0:nxmaxn-1,0:nymaxn-1),nsss(0:nxmaxn-1,0:nymaxn-1)
    7674  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
    77   real :: conversion_factor !added by mc to make it consistent with new gridchek.f90
    7875
    7976  logical :: hflswitch,strswitch
    … …  
    137134  endif
    138135
    139   conversion_factor=1.
    140 
    141 
    142136  else
    143137
    … …  
    154148  call grib_get_int(igrib,'level',valSurf,iret)
    155149  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
    158150
    159151  !print*,discipl,parCat,parNum,typSurf,valSurf
    … …  
    164156  isec1(8)=-1
    165157  isec1(8)=valSurf     ! level
    166    conversion_factor=1.
    167158  if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
    168159    isec1(6)=130         ! indicatorOfParameter
    … …  
    175166  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
    176167    isec1(6)=134         ! indicatorOfParameter
    177   elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot !
     168  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
    178169    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
    181170  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    182171    isec1(6)=151         ! indicatorOfParameter
    … …  
    191180  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    192181    isec1(6)=141         ! indicatorOfParameter
    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
     182  elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
    195183    isec1(6)=164         ! indicatorOfParameter
    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
     184  elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
    197185    isec1(6)=142         ! indicatorOfParameter
    198186  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    199187    isec1(6)=143         ! indicatorOfParameter
    200     conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
    201188  elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
    202189    isec1(6)=146         ! indicatorOfParameter
    203190  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    204191    isec1(6)=176         ! indicatorOfParameter
    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
     192  elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
    206193    isec1(6)=180         ! indicatorOfParameter
    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
     194  elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
    208195    isec1(6)=181         ! indicatorOfParameter
    209196  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    210197    isec1(6)=129         ! indicatorOfParameter
    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
     198  elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
    212199    isec1(6)=160         ! indicatorOfParameter
    213200  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    … …  
    215202    isec1(6)=172         ! indicatorOfParameter
    216203  else
    217     print*,'***WARNING: undefined GRiB2 message found!',discipl, &
     204    print*,'***ERROR: undefined GRiB2 message found!',discipl, &
    218205         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
    223206  endif
    224207
    … …  
    244227  ! CHECK GRID SPECIFICATIONS
    245228  if(isec2(2).ne.nxn(l)) stop &
    246   'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
     229       'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
    247230  if(isec2(3).ne.nyn(l)) stop &
    248   'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
     231       'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
    249232  if(isec2(12)/2-1.ne.nlev_ec) stop 'READWIND: VERTICAL DISCRET&
    250   IZATION NOT CONSISTENT FOR A NESTING LEVEL'
     233       &IZATION NOT CONSISTENT FOR A NESTING LEVEL'
    251234  endif ! ifield
    252235
    253236  !HSO  get the second part of the grid dimensions only from GRiB1 messages
    254  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consisitent with new readwind.f90
     237  if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
    255238    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
    256239         xauxin,iret)
    … …  
    259242         yauxin,iret)
    260243    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 
    264244    xaux=xauxin
    265245    yaux=yauxin
    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'
     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'
    270258    gotGrid=1
    271259  endif
    … …  
    292280             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
    293281        if(isec1(6).eq.141) sdn(i,j,1,n,l)= &!! SNOW DEPTH
    294              zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90!
     282             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
    295283        if(isec1(6).eq.151) msln(i,j,1,n,l)= &!! SEA LEVEL PRESS.
    296284             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
    … …  
    310298        endif
    311299        if(isec1(6).eq.143) then                         !! CONVECTIVE PREC.
    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
     300          convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
    313301          if (convprecn(i,j,1,n,l).lt.0.) convprecn(i,j,1,n,l)=0.
    314302        endif

  • /trunk/src/timemanager.f90

    r30 r20  
    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                         *
    47   !  Changes Stefan Henne, Harald Sodemann, 2013-2014                          *
    48   !   added netcdf output code                                                 *
     45  !new interface between flexpart and convection scheme                        *
     46  !Emanuel's latest subroutine convect43c.f is used                            *
    4947  !*****************************************************************************
    5048  !                                                                            *
    5149  ! Variables:                                                                 *
    52   ! dep                .true. if either wet or dry deposition is switched on   *
     50  ! DEP                .true. if either wet or dry deposition is switched on   *
    5351  ! decay(maxspec) [1/s] decay constant for radioactive decay                  *
    54   ! drydep             .true. if dry deposition is switched on                 *
     52  ! DRYDEP             .true. if dry deposition is switched on                 *
    5553  ! ideltas [s]        modelling period                                        *
    5654  ! itime [s]          actual temporal position of calculation                 *
    … …  
    7270  ! prob               probability of absorption at ground due to dry          *
    7371  !                    deposition                                              *
    74   ! wetdep             .true. if wet deposition is switched on                 *
     72  ! WETDEP             .true. if wet deposition is switched on                 *
    7573  ! weight             weight for each concentration sample (1/2 or 1)         *
    7674  ! uap(maxpart),ucp(maxpart),uzp(maxpart) = random velocities due to          *
    … …  
    9492  use par_mod
    9593  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
    9994
    10095  implicit none
    … …  
    134129
    135130
    136   itime=0
    137   !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
    138   write(*,46) float(itime)/3600,itime,numpart
    139131  if (verbosity.gt.0) then
    140132    write (*,*) 'timemanager> starting simulation'
    … …  
    146138
    147139  do itime=0,ideltas,lsynctime
     140
    148141
    149142  ! Computation of wet deposition, OH reaction and mass transfer
    … …  
    157150  !********************************************************************
    158151
    159     if (wetdep .and. itime .ne. 0 .and. numpart .gt. 0) then
     152    if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
    160153        if (verbosity.gt.0) then
    161154           write (*,*) 'timemanager> call wetdepo'
    … …  
    164157    endif
    165158
    166     if (ohrea .and. itime .ne. 0 .and. numpart .gt. 0) &
     159    if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) &
    167160         call ohreaction(itime,lsynctime,loutnext)
    168161
    169     if (assspec .and. itime .ne. 0 .and. numpart .gt. 0) then
     162    if (ASSSPEC .and. itime .ne. 0 .and. numpart .gt. 0) then
    170163       stop 'associated species not yet implemented!'
    171164  !     call transferspec(itime,lsynctime,loutnext)
    … …  
    245238  !***********************************************************************
    246239
    247     if (dep.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
     240    if (DEP.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
    248241      do ks=1,nspec
    249242      do kp=1,maxpointspec_act
    … …  
    355348        if ((iout.le.3.).or.(iout.eq.5)) then
    356349          if (surf_only.ne.1) then
    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
     350          call concoutput(itime,outnum,gridtotalunc, &
     351               wetgridtotalunc,drygridtotalunc)
    364352          else 
    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
     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
    382365          endif
    383366
    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
     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)
    401369          outnum=0.
    402370        endif
    403371        if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
    404372        if (iflux.eq.1) call fluxoutput(itime)
    405         !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
    406         write(*,46) float(itime)/3600,itime,numpart
    407 45      format(i9,' SECONDS SIMULATED: ',i8, ' PARTICLES:    Uncertainty: ',3f7.3)
    408 46      format(' Simulated ',f7.1,' hours (',i9,' s), ',i8, ' particles')
     373        write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc, &
     374             drygridtotalunc
     37545      format(i9,' SECONDS SIMULATED: ',i8, &
     376             ' PARTICLES:    Uncertainty: ',3f7.3)
    409377        if (ipout.ge.1) call partoutput(itime)    ! dump particle positions
    410378        loutnext=loutnext+loutstep
    … …  
    504472  !****************************
    505473
    506         xold=real(xtra1(j))
    507         yold=real(ytra1(j))
     474        xold=xtra1(j)
     475        yold=ytra1(j)
    508476        zold=ztra1(j)
    509477
    … …  
    545513            endif
    546514
    547             if (drydepspec(ks)) then        ! dry deposition
     515            if (DRYDEPSPEC(ks)) then        ! dry deposition
    548516              drydeposit(ks)=xmass1(j,ks)*prob(ks)*decfact
    549517              xmass1(j,ks)=xmass1(j,ks)*(1.-prob(ks))*decfact
    … …  
    556524            endif
    557525
     526
    558527            if (mdomainfill.eq.0) then
    559528              if (xmass(npoint(j),ks).gt.0.) &
    … …  
    567536          if (xmassfract.lt.0.0001) then   ! terminate all particles carrying less mass
    568537            itra1(j)=-999999999
    569             if (verbosity.gt.0) then
    570               print*,'terminated particle ',j,' for small mass'
    571             endif
    572538          endif
    573539
    574540  !        Sabine Eckhardt, June 2008
    575541  !        don't create depofield for backward runs
    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
     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)
    581548          endif
    582549
    … …  
    585552
    586553          if (abs(itra1(j)-itramem(j)).ge.lage(nageclass)) then
    587             if (linit_cond.ge.1) call initial_cond_calc(itime+lsynctime,j)
     554            if (linit_cond.ge.1) &
     555                 call initial_cond_calc(itime+lsynctime,j)
    588556            itra1(j)=-999999999
    589             if (verbosity.gt.0) then
    590               print*,'terminated particle ',j,' for age'
    591             endif
    592557          endif
    593558        endif
    … …  
    617582
    618583  if (iflux.eq.1) then
    619     deallocate(flux)
     584      deallocate(flux)
    620585  endif
    621   if (ohrea.eqv..TRUE.) then
    622     deallocate(OH_field,OH_field_height)
     586  if (OHREA.eqv..TRUE.) then
     587      deallocate(OH_field,OH_field_height)
    623588  endif
    624589  if (ldirect.gt.0) then
    625     deallocate(drygridunc,wetgridunc)
     590  deallocate(drygridunc,wetgridunc)
    626591  endif
    627592  deallocate(gridunc)
    … …  
    630595  deallocate(xmasssave)
    631596  if (nested_output.eq.1) then
    632     deallocate(orooutn, arean, volumen)
    633     if (ldirect.gt.0) then
    634       deallocate(griduncn,drygriduncn,wetgriduncn)
    635     endif
     597     deallocate(orooutn, arean, volumen)
     598     if (ldirect.gt.0) then
     599     deallocate(griduncn,drygriduncn,wetgriduncn)
     600     endif
    636601  endif
    637602  deallocate(outheight,outheighthalf)
    638   deallocate(oroout,area,volume)
     603  deallocate(oroout, area, volume)
    639604
    640605end subroutine timemanager

  • /trunk/src/verttransform.f90

    r30 r20  
    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
    5153  !*****************************************************************************
    5254  !                                                                            *
    … …  
    7072
    7173  integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
    72   integer :: rain_cloud_above,kz_inv
     74 integer :: rain_cloud_above,kz_inv !SE
     75  integer icloudtop !PS
    7376  real :: f_qvsat,pressure
    74   real :: rh,lsp,convp
    75   real :: rhoh(nuvzmax),pinmconv(nzmax)
     77 !real :: rh,lsp,convp
     78  real :: rh,lsp,convp,prec,rhmin 
     79  real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
    7680  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
    7781  real :: xlon,ylat,xlonr,dzdx,dzdy
    78   real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
     82  real :: dzdx1,dzdx2,dzdy1,dzdy2, precmin
    7983  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
    8084  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
    … …  
    8387  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
    8488  real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
     89  logical lconvectprec
    8590  real,parameter :: const=r_air/ga
     91  parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
    8692
    8793  logical :: init = .true.
    … …  
    9197  ! If verttransform is called the first time, initialize heights of the   *
    9298  ! z levels in meter. The heights are the heights of model levels, where  *
    93   ! u,v,T and qv are given.                                                *
     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.                                                  *
    94104  !*************************************************************************
     105
     106
     107  !      do 897 kz=1,nuvz
     108  !       write (*,*) 'akz: ',akz(kz),'bkz',bkz(kz)
     109  !897     continue
    95110
    96111  if (init) then
    … …  
    98113  ! Search for a point with high surface pressure (i.e. not above significant topography)
    99114  ! Then, use this point to construct a reference z profile, to be used at all times
    100   !**************************************************************************************
     115  !*****************************************************************************
    101116
    102117    do jy=0,nymin1
    … …  
    113128
    114129    tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
    115     ps(ixm,jym,1,n))
     130         ps(ixm,jym,1,n))
    116131    pold=ps(ixm,jym,1,n)
    117132    height(1)=0.
    … …  
    121136      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
    122137
     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
    123148      if (abs(tv-tvold).gt.0.2) then
    124         height(kz)=height(kz-1)+const*log(pold/pint)* &
    125         (tv-tvold)/log(tv/tvold)
     149        height(kz)= &
     150             height(kz-1)+const*log(pold/pint)* &
     151             (tv-tvold)/log(tv/tvold)
    126152      else
    127         height(kz)=height(kz-1)+const*log(pold/pint)*tv
     153        height(kz)=height(kz-1)+ &
     154             const*log(pold/pint)*tv
    128155      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
    129168
    130169      tvold=tv
    131170      pold=pint
    132171    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
    133180
    134181
    … …  
    157204  do jy=0,nymin1
    158205    do ix=0,nxmin1
    159       tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ps(ix,jy,1,n))
     206      tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
     207           ps(ix,jy,1,n))
    160208      pold=ps(ix,jy,1,n)
    161       uvwzlev(ix,jy,1)=0.
     209      uvzlev(1)=0.
    162210      wzlev(1)=0.
    163211      rhoh(1)=pold/(r_air*tvold)
    … …  
    173221
    174222        if (abs(tv-tvold).gt.0.2) then
    175           uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
    176           (tv-tvold)/log(tv/tvold)
     223          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
     224               (tv-tvold)/log(tv/tvold)
    177225        else
    178           uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
     226          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
    179227        endif
    180228
    … …  
    185233
    186234      do kz=2,nwz-1
    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)
     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
    190252
    191253  ! pinmconv=(h2-h1)/(p2-p1)
    192254
    193255      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
    194       ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
    195       (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
     256           ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
     257           (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
    196258      do kz=2,nz-1
    197259        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
    198         ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
    199         (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
     260             ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
     261             (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
    200262      end do
    201263      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
    202       ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
    203       (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
     264           ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
     265           (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
    204266
    205267  ! Levels, where u,v,t and q are given
    … …  
    221283      do iz=2,nz-1
    222284        do kz=kmin,nuvz
    223           if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
     285          if(height(iz).gt.uvzlev(nuvz)) then
    224286            uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
    225287            vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
    … …  
    230292            goto 30
    231293          endif
    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)
     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)
    236298           dz=dz1+dz2
    237299           uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
    … …  
    277339
    278340      drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
    279       (height(2)-height(1))
     341           (height(2)-height(1))
    280342      do kz=2,nz-1
    281343        drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
    282         (height(kz+1)-height(kz-1))
     344             (height(kz+1)-height(kz-1))
    283345      end do
    284346      drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
    … …  
    294356
    295357  do jy=1,ny-2
    296     cosf=cos((real(jy)*dy+ylat0)*pi180)
    297358    do ix=1,nx-2
    298359
    … …  
    300361      do iz=2,nz-1
    301362
    302         ui=uu(ix,jy,iz,n)*dxconst/cosf
     363        ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
    303364        vi=vv(ix,jy,iz,n)*dyconst
    304365
    305366        do kz=kmin,nz
    306367          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
    307           (height(iz).le.uvwzlev(ix,jy,kz))) then
     368               (height(iz).le.uvwzlev(ix,jy,kz))) then
    308369            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
    309370            dz2=uvwzlev(ix,jy,kz)-height(iz)
    … …  
    348409        do iz=1,nz
    349410          call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
    350           vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
     411               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
     412               vvpol(ix,jy,iz,n))
    351413        end do
    352414      end do
    … …  
    362424      xlon=xlon0+real(nx/2-1)*dx
    363425      xlonr=xlon*pi/180.
    364       ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2)
     426      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
     427           vv(nx/2-1,nymin1,iz,n)**2)
    365428      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
    366         ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
     429        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
     430             vv(nx/2-1,nymin1,iz,n))-xlonr
    367431      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
    368432        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
    369         vv(nx/2-1,nymin1,iz,n))-xlonr
     433             vv(nx/2-1,nymin1,iz,n))-xlonr
    370434      else
    371435        ddpol=pi/2-xlonr
    … …  
    380444      uuaux=-ffpol*sin(xlonr+ddpol)
    381445      vvaux=-ffpol*cos(xlonr+ddpol)
    382       call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
     446      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
     447           vvpolaux)
     448
    383449      jy=nymin1
    384450      do ix=0,nxmin1
    … …  
    419485        do iz=1,nz
    420486          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
    421           vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
     487               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
     488               vvpol(ix,jy,iz,n))
    422489        end do
    423490      end do
    … …  
    432499      xlon=xlon0+real(nx/2-1)*dx
    433500      xlonr=xlon*pi/180.
    434       ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2)
     501      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
     502           vv(nx/2-1,0,iz,n)**2)
    435503      if (vv(nx/2-1,0,iz,n).lt.0.) then
    436         ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
     504        ddpol=atan(uu(nx/2-1,0,iz,n)/ &
     505             vv(nx/2-1,0,iz,n))+xlonr
    437506      else if (vv(nx/2-1,0,iz,n).gt.0.) then
    438         ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
     507        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
     508             vv(nx/2-1,0,iz,n))+xlonr
    439509      else
    440510        ddpol=pi/2-xlonr
    … …  
    449519      uuaux=+ffpol*sin(xlonr-ddpol)
    450520      vvaux=-ffpol*cos(xlonr-ddpol)
    451       call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
     521      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
     522           vvpolaux)
    452523
    453524      jy=0
    … …  
    480551  !   create a cloud and rainout/washout field, clouds occur where rh>80%
    481552  !   total cloudheight is stored at level 0
     553
     554
     555
    482556  do jy=0,nymin1
    483557    do ix=0,nxmin1
    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
     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
     61198        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
    505620            endif
    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
     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
    513643            endif
    514          endif
    515       end do
     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
    516657    end do
    517658  end do
    518659
     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
    519685end subroutine verttransform

  • /trunk/src/verttransform_gfs.f90

    r30 r20  
    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 :: rhoh(nuvzmax),pinmconv(nzmax)
     73  real :: uvzlev(nuvzmax),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,cosf
     76  real :: dzdx1,dzdx2,dzdy1,dzdy2
    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.                                                *
     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.                                                  *
    9599  !*************************************************************************
    96100
    … …  
    114118
    115119    tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
    116     ps(ixm,jym,1,n))
     120         ps(ixm,jym,1,n))
    117121    pold=ps(ixm,jym,1,n)
    118122    height(1)=0.
    … …  
    122126      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
    123127
     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
    124138      if (abs(tv-tvold).gt.0.2) then
    125         height(kz)=height(kz-1)+const*log(pold/pint)* &
    126         (tv-tvold)/log(tv/tvold)
     139        height(kz)= &
     140             height(kz-1)+const*log(pold/pint)* &
     141             (tv-tvold)/log(tv/tvold)
    127142      else
    128         height(kz)=height(kz-1)+const*log(pold/pint)*tv
     143        height(kz)=height(kz-1)+ &
     144             const*log(pold/pint)*tv
    129145      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
    130158
    131159      tvold=tv
    132160      pold=pint
    133161    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
    134170
    135171
    … …  
    162198      llev = 0
    163199      do i=1,nuvz
    164         if (ps(ix,jy,1,n).lt.akz(i)) llev=i
     200       if (ps(ix,jy,1,n).lt.akz(i)) llev=i
    165201      end do
    166202       llev = llev+1
    … …  
    176212      tvold=tth(ix,jy,llev,n)*(1.+0.608*qvh(ix,jy,llev,n))
    177213      pold=akz(llev)
     214      uvzlev(llev)=0.
    178215      wzlev(llev)=0.
    179216      uvwzlev(ix,jy,llev)=0.
    … …  
    186223
    187224        if (abs(tv-tvold).gt.0.2) then
    188           uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
    189           (tv-tvold)/log(tv/tvold)
     225          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
     226               (tv-tvold)/log(tv/tvold)
    190227        else
    191           uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
     228          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
    192229        endif
    193         wzlev(kz)=uvwzlev(ix,jy,kz)
     230        wzlev(kz)=uvzlev(kz)
     231        uvwzlev(ix,jy,kz)=uvzlev(kz)
    194232
    195233        tvold=tv
    196234        pold=pint
    197235      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
    198245
    199246  ! pinmconv=(h2-h1)/(p2-p1)
    … …  
    232279      do iz=2,nz-1
    233280        do kz=kmin,nuvz
    234           if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
     281          if(height(iz).gt.uvzlev(nuvz)) then
    235282            uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
    236283            vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
    … …  
    242289            goto 30
    243290          endif
    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
     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
    258305          endif
    259306        end do
    … …  
    271318        do kz=kmin,nwz
    272319          if ((height(iz).gt.wzlev(kz-1)).and. &
    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
     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
    279327          endif
    280328        end do
    … …  
    289337      do kz=2,nz-1
    290338        drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
    291         (height(kz+1)-height(kz-1))
     339             (height(kz+1)-height(kz-1))
    292340      end do
    293341      drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
    … …  
    303351
    304352  do jy=1,ny-2
    305     cosf=cos((real(jy)*dy+ylat0)*pi180)
    306353    do ix=1,nx-2
    307354
    … …  
    320367      do iz=2,nz-1
    321368
    322         ui=uu(ix,jy,iz,n)*dxconst/cosf
     369        ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
    323370        vi=vv(ix,jy,iz,n)*dyconst
    324371
    325372        do kz=kmin,nz
    326373          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
    327           (height(iz).le.uvwzlev(ix,jy,kz))) then
     374               (height(iz).le.uvwzlev(ix,jy,kz))) then
    328375            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
    329376            dz2=uvwzlev(ix,jy,kz)-height(iz)
    … …  
    379426      xlon=xlon0+real(nx/2-1)*dx
    380427      xlonr=xlon*pi/180.
    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
     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
    384433      elseif (vv(nx/2-1,nymin1,iz,n).gt.0.) then
    385434        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
    386         vv(nx/2-1,nymin1,iz,n))-xlonr
     435             vv(nx/2-1,nymin1,iz,n))-xlonr
    387436      else
    388437        ddpol=pi/2-xlonr
    … …  
    397446      uuaux=-ffpol*sin(xlonr+ddpol)
    398447      vvaux=-ffpol*cos(xlonr+ddpol)
    399       call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
     448      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
     449           vvpolaux)
     450
    400451      jy=nymin1
    401452      do ix=0,nxmin1
    … …  
    409460  ! ward parallel of latitude
    410461
    411     do iz=1,nz
     462  do iz=1,nz
    412463      wdummy=0.
    413464      jy=ny-2
    … …  
    420471        ww(ix,jy,iz,n)=wdummy
    421472      end do
    422     end do
     473  end do
    423474
    424475  endif
    … …  
    436487        do iz=1,nz
    437488          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
    438           vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
     489               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
     490               vvpol(ix,jy,iz,n))
    439491        end do
    440492      end do
    … …  
    446498      xlon=xlon0+real(nx/2-1)*dx
    447499      xlonr=xlon*pi/180.
    448       ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2)
     500      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
     501           vv(nx/2-1,0,iz,n)**2)
    449502      if(vv(nx/2-1,0,iz,n).lt.0.) then
    450         ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
     503        ddpol=atan(uu(nx/2-1,0,iz,n)/ &
     504             vv(nx/2-1,0,iz,n))+xlonr
    451505      elseif (vv(nx/2-1,0,iz,n).gt.0.) then
    452         ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))-xlonr
     506        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
     507             vv(nx/2-1,0,iz,n))-xlonr
    453508      else
    454509        ddpol=pi/2-xlonr
    … …  
    463518      uuaux=+ffpol*sin(xlonr-ddpol)
    464519      vvaux=-ffpol*cos(xlonr-ddpol)
    465       call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
     520      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
     521           vvpolaux)
    466522
    467523      jy=0
    … …  
    506562         clouds(ix,jy,kz,n)=0
    507563         if (rh.gt.0.8) then ! in cloud
    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
     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
    519576         else ! no cloud
    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
     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
    527584         endif
    528585      end do

  • /trunk/src/verttransform_nests.f90

    r30 r20  
    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 :: wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax)
     74  real :: uvzlev(nuvzmax),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,cosf
     76  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
    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       uvwzlev(ix,jy,1)=0.
     100      uvzlev(1)=0.
    101101      wzlev(1)=0.
    102102      rhoh(1)=pold/(r_air*tvold)
    … …  
    112112
    113113        if (abs(tv-tvold).gt.0.2) then
    114           uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
    115           (tv-tvold)/log(tv/tvold)
     114          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
     115               (tv-tvold)/log(tv/tvold)
    116116        else
    117           uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
     117          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
    118118        endif
    119119
    … …  
    124124
    125125      do kz=2,nwz-1
    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)
     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
    129151
    130152  ! pinmconv=(h2-h1)/(p2-p1)
    131153
    132154      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
    133       ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
    134       (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
     155           ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
     156           (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
    135157      do kz=2,nz-1
    136158        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
    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)))
     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)))
    139161      end do
    140162      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
    141       ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
    142       (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
     163           ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
     164           (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
    143165
    144166
    … …  
    161183      do iz=2,nz-1
    162184        do kz=kmin,nuvz
    163           if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
     185          if(height(iz).gt.uvzlev(nuvz)) then
    164186            uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
    165187            vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
    … …  
    170192            goto 30
    171193          endif
    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)
     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)
    176198           dz=dz1+dz2
    177199           uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
    178            uuhn(ix,jy,kz,l)*dz1)/dz
     200                uuhn(ix,jy,kz,l)*dz1)/dz
    179201           vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
    180            vvhn(ix,jy,kz,l)*dz1)/dz
     202                vvhn(ix,jy,kz,l)*dz1)/dz
    181203           ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
    182            tthn(ix,jy,kz,n,l)*dz1)/dz
     204                tthn(ix,jy,kz,n,l)*dz1)/dz
    183205           qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
    184            qvhn(ix,jy,kz,n,l)*dz1)/dz
     206                qvhn(ix,jy,kz,n,l)*dz1)/dz
    185207           pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
    186            pvhn(ix,jy,kz,l)*dz1)/dz
     208                pvhn(ix,jy,kz,l)*dz1)/dz
    187209           rhon(ix,jy,iz,n,l)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
    188210           kmin=kz
    … …  
    203225        do kz=kmin,nwz
    204226          if ((height(iz).gt.wzlev(kz-1)).and. &
    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
     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
    213235          endif
    214236        end do
    … …  
    220242
    221243      drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
    222       (height(2)-height(1))
     244           (height(2)-height(1))
    223245      do kz=2,nz-1
    224246        drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
    225         rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
     247             rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
    226248      end do
    227249      drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
    … …  
    237259
    238260  do jy=1,nyn(l)-2
    239     cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
    240261    do ix=1,nxn(l)-2
    241262
    … …  
    243264      do iz=2,nz-1
    244265
    245         ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf
     266        ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/ &
     267             cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
    246268        vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)
    247269
    … …  
    297319               rain_cloud_above=1
    298320               cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ &
    299                height(kz)-height(kz-1)
     321                    height(kz)-height(kz-1)
    300322               if (lsp.ge.convp) then
    301323                  cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout

  • /trunk/src/wetdepo.f90

    r30 r20  
    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
    4144  !*****************************************************************************
    4245  !                                                                            *
    … …  
    7174
    7275  integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
    73   integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
    74   integer :: ks, kp
     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
    7580  real :: S_i, act_temp, cl, cle ! in cloud scavenging
    7681  real :: clouds_h ! cloud height for the specific grid point
    77   real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
     82  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav, precsub,f
    7883  real :: wetdeposit(maxspec),restmass
    7984  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
    … …  
    98103
    99104  do jpart=1,numpart
    100 
    101105    if (itra1(jpart).eq.-999999999) goto 20
    102106    if(ldirect.eq.1)then
    … …  
    110114      if (itage.lt.lage(nage)) goto 33
    111115    end do
    112 33  continue
     11633   continue
    113117
    114118
    … …  
    119123    do j=numbnests,1,-1
    120124      if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
    121       (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
     125           (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
    122126        ngrid=j
    123127        goto 23
    124128      endif
    125129    end do
    126 23  continue
     13023   continue
    127131
    128132
    … …  
    147151    interp_time=nint(itime-0.5*ltsample)
    148152
    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
     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
    160182  ! get the level were the actual particle is in
    161     do il=2,nz
    162       if (height(il).gt.ztra1(jpart)) then
    163         hz=il-1
    164         goto 26
    165       endif
    166     end do
    167 26  continue
    168 
    169     n=memind(2)
    170     if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
    171     n=memind(1)
     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
     19026     continue
     191
     192  n=memind(2)
     193  if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
     194       n=memind(1)
    172195
    173196  ! if there is no precipitation or the particle is above the clouds no
    174197  ! scavenging is done
    175198
    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'
     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
    187226
    188227  ! 1) Parameterization of the the area fraction of the grid cell where the
    … …  
    228267  !**********************************************************
    229268
    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
     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
    256293
    257294! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening
    … …  
    260297! wetc_in=0.9 (default)
    261298! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling)
    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
    268         endif
    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
     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
     309        endif
     310
     311
     312  !      if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
     313  !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
    284314        wetdeposit(ks)=xmass1(jpart,ks)* &
    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
     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
    316414  ! Add the wet deposition to accumulated amount on output grid and nested output grid
    317415  !*****************************************************************************
    318416
    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
     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
    325434
    32643520  continue
    327   end do ! all particles
     436  end do
    328437
    329438end subroutine wetdepo

  • /trunk/src/writeheader_nest.f90

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