Changeset 24


Ignore:
Timestamp:
May 23, 2014, 11:48:41 AM (10 years ago)
Author:
igpis
Message:

version 9.2 beta. Changes from HH, AST, MC, NIK, IP. Changes in vert transform. New SPECIES input includes scavenging coefficients

Location:
trunk
Files:
1 deleted
47 edited

Legend:

Unmodified
Added
Removed

  • trunk/options/COMMAND

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

  • trunk/options/COMMAND.alternative

    r4 r24  
    30300               NESTED_OUTPUT     SHALL NESTED OUTPUT BE USED? YES, 0 NO
    31312               LINIT_COND        INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT
     320               SURF_ONLY         IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER
    3233
    3334
    … …  
    117118    conditions shall be calculated and written to output files
    118119    0=no output, 1 or 2 determines in which units the initial conditions are provided.
     120
     12125. SURF_ONLY: When set to 1, concentration/emission sensitivity is written out only
     122    for the surface layer; useful for instance when only footprint emission sensitivity is needed
     123    but initial conditions are needed on a full 3-D grid

  • trunk/options/COMMAND.reference

    r4 r24  
    1111
    12122. ________ ______   3X, I8, 1X, I6
    13    20040720 000000
     13   20110310 000000
    1414   YYYYMMDD HHMISS   BEGINNING DATE OF SIMULATION
    1515
    16163. ________ ______   3X, I8, 1X, I6
    17    20040721 120000
     17   20110310 120000
    1818   YYYYMMDD HHMISS   ENDING DATE OF SIMULATION
    1919
    … …  
    101101    2
    102102    LINIT_COND       INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT
     103
     10425. _                4X, I1
     105    0
     106    SURF_ONLY        IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER
    103107
    104108
    … …  
    189193    conditions shall be calculated and written to output files
    190194    0=no output, 1 or 2 determines in which units the initial conditions are provided.
     195
     19625. SURF_ONLY: When set to 1, concentration/emission sensitivity is written out only
     197    for the surface layer; useful for instance when only footprint emission sensitivity is needed
     198    but initial conditions are needed on a full 3-D grid

  • trunk/options/RELEASES

    r4 r24  
    1 +++++++++++++++++ HEADER ++++++++++++++++++++
    2 +++++++++++++++++ HEADER ++++++++++++++++++++
    3 +++++++++++++++++ HEADER ++++++++++++++++++++
    4 +++++++++++++++++ HEADER ++++++++++++++++++++
    5 +++++++++++++++++ HEADER ++++++++++++++++++++
    6 +++++++++++++++++ HEADER ++++++++++++++++++++
    7 +++++++++++++++++ HEADER ++++++++++++++++++++
    8 +++++++++++++++++ HEADER ++++++++++++++++++++
    9 +++++++++++++++++ HEADER ++++++++++++++++++++
    10 +++++++++++++++++ HEADER ++++++++++++++++++++
    11 +++++++++++++++++ HEADER ++++++++++++++++++++
    12 5
    13 31
    14 31
    15 31
    16 31
    17 31
     1*************************************************************************
     2*                                                                       *
     3*                                                                       *
     4*                                                                       *
     5*   Input file for the Lagrangian particle dispersion model FLEXPART    *
     6*                        Please select your options                     *
     7*                                                                       *
     8*                                                                       *
     9*                                                                       *
     10*************************************************************************
     11+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     12  1     
     13___                        i3    Total number of species emitted
    1814
    19 20040101 080000
    20 20040102 080000
    21    8.2500
    22   58.3833
    23    8.2500
    24   58.3833
    25 1
    26      0.000
    27    100.000
    28 4000
    29 1
    30 1
    31 1
    32 1
    33 1
    34 rel_bj
     15  3
     16___                        i3    Index of species in file SPECIES
    3517
    36 20040108 080000
    37 20040109 080000
    38    8.2500
    39   58.3833
    40    8.2500
    41   58.3833
    42 1
    43      0.000
    44    100.000
    45 4000
    46 1
    47 1
    48 1
    49 1
    50 1
    51 rel_bj
     18=========================================================================
     1920111210 000001
     20________ ______            i8,1x,i6 Beginning date and time of release
    5221
    53 20040115 080000
    54 20040116 080000
    55    8.2500
    56   58.3833
    57    8.2500
    58   58.3833
    59 1
    60      0.000
    61    100.000
    62 4000
    63 1
    64 1
    65 1
    66 1
    67 1
    68 rel_bj
     2220111210 090000
     23________ ______            i8,1x,i6 Ending date and time of release
    6924
    70 20040122 080000
    71 20040123 080000
    72    8.2500
    73   58.3833
    74    8.2500
    75   58.3833
    76 1
    77      0.000
    78    100.000
    79 4000
    80 1
    81 1
    82 1
    83 1
    84 1
    85 rel_bj
     25   9.4048
     26____.____                  f9.4  Longitude [DEG] of lower left corner
    8627
    87 20040129 080000
    88 20040130 080000
    89    8.2500
    90   58.3833
    91    8.2500
    92   58.3833
    93 1
    94      0.000
    95    100.000
    96 4000
    97 1
    98 1
    99 1
    100 1
    101 1
    102 rel_bj
     28  48.5060
     29____.____                  f9.4  Latitude [DEG] of lower left corner
    10330
    104 20040205 080000
    105 20040206 080000
    106    8.2500
    107   58.3833
    108    8.2500
    109   58.3833
    110 1
    111      0.000
    112    100.000
    113 4000
    114 1
    115 1
    116 1
    117 1
    118 1
    119 rel_bj
     31   9.5067
     32____.____                  f9.4  Longitude [DEG] of upper right corner
    12033
    121 20040212 080000
    122 20040213 080000
    123    8.2500
    124   58.3833
    125    8.2500
    126   58.3833
    127 1
    128      0.000
    129    100.000
    130 4000
    131 1
    132 1
    133 1
    134 1
    135 1
    136 rel_bj
     34  48.5158
     35____.____                  f9.4  Latitude [DEG] of upper right corner
    13736
    138 20040219 080000
    139 20040220 080000
    140    8.2500
    141   58.3833
    142    8.2500
    143   58.3833
    144 1
    145      0.000
    146    100.000
    147 4000
    148 1
    149 1
    150 1
    151 1
    152 1
    153 rel_bj
     37        2
     38_________                  i9    1 for m above ground, 2 for m above sea level
    15439
    155 20040226 080000
    156 20040227 080000
    157    8.2500
    158   58.3833
    159    8.2500
    160   58.3833
    161 1
    162      0.000
    163    100.000
    164 4000
    165 1
    166 1
    167 1
    168 1
    169 1
    170 rel_bj
     40    0.00
     41_____.___                  f10.3 Lower z-level (in m agl or m asl)
     42 
     43   10.00
     44_____.___                  f10.3 Upper z-level (in m agl or m asl)
     45 
     46    20000               
     47_________                  i9    Total number of particles to be released
    17148
    172 20040304 080000
    173 20040305 080000
    174    8.2500
    175   58.3833
    176    8.2500
    177   58.3833
    178 1
    179      0.000
    180    100.000
    181 4000
    182 1
    183 1
    184 1
    185 1
    186 1
    187 rel_bj
     491.0000E00
     50_.____E__                  e9.4  Total mass emitted
    18851
    189 20040311 080000
    190 20040312 080000
    191    8.2500
    192   58.3833
    193    8.2500
    194   58.3833
    195 1
    196      0.000
    197    100.000
    198 4000
    199 1
    200 1
    201 1
    202 1
    203 1
    204 rel_bj
    205 
    206 20040318 080000
    207 20040319 080000
    208    8.2500
    209   58.3833
    210    8.2500
    211   58.3833
    212 1
    213      0.000
    214    100.000
    215 4000
    216 1
    217 1
    218 1
    219 1
    220 1
    221 rel_bj
    222 
    223 20040325 080000
    224 20040326 080000
    225    8.2500
    226   58.3833
    227    8.2500
    228   58.3833
    229 1
    230      0.000
    231    100.000
    232 4000
    233 1
    234 1
    235 1
    236 1
    237 1
    238 rel_bj
    239 
    240 20040401 080000
    241 20040402 080000
    242    8.2500
    243   58.3833
    244    8.2500
    245   58.3833
    246 1
    247      0.000
    248    100.000
    249 4000
    250 1
    251 1
    252 1
    253 1
    254 1
    255 rel_bj
    256 
    257 20040408 080000
    258 20040409 080000
    259    8.2500
    260   58.3833
    261    8.2500
    262   58.3833
    263 1
    264      0.000
    265    100.000
    266 4000
    267 1
    268 1
    269 1
    270 1
    271 1
    272 rel_bj
    273 
    274 20040415 080000
    275 20040416 080000
    276    8.2500
    277   58.3833
    278    8.2500
    279   58.3833
    280 1
    281      0.000
    282    100.000
    283 4000
    284 1
    285 1
    286 1
    287 1
    288 1
    289 rel_bj
    290 
    291 20040422 080000
    292 20040423 080000
    293    8.2500
    294   58.3833
    295    8.2500
    296   58.3833
    297 1
    298      0.000
    299    100.000
    300 4000
    301 1
    302 1
    303 1
    304 1
    305 1
    306 rel_bj
    307 
    308 20040429 080000
    309 20040430 080000
    310    8.2500
    311   58.3833
    312    8.2500
    313   58.3833
    314 1
    315      0.000
    316    100.000
    317 4000
    318 1
    319 1
    320 1
    321 1
    322 1
    323 rel_bj
    324 
    325 20040609 080000
    326 20040610 080000
    327    8.2500
    328   58.3833
    329    8.2500
    330   58.3833
    331 1
    332      0.000
    333    100.000
    334 4000
    335 1
    336 1
    337 1
    338 1
    339 1
    340 rel_bj
    341 
    342 20040610 080000
    343 20040611 080000
    344    8.2500
    345   58.3833
    346    8.2500
    347   58.3833
    348 1
    349      0.000
    350    100.000
    351 4000
    352 1
    353 1
    354 1
    355 1
    356 1
    357 rel_bj
    358 
    359 20040617 080000
    360 20040618 080000
    361    8.2500
    362   58.3833
    363    8.2500
    364   58.3833
    365 1
    366      0.000
    367    100.000
    368 4000
    369 1
    370 1
    371 1
    372 1
    373 1
    374 rel_bj
    375 
    376 20040624 080000
    377 20040625 080000
    378    8.2500
    379   58.3833
    380    8.2500
    381   58.3833
    382 1
    383      0.000
    384    100.000
    385 4000
    386 1
    387 1
    388 1
    389 1
    390 1
    391 rel_bj
    392 
    393 20040701 080000
    394 20040702 080000
    395    8.2500
    396   58.3833
    397    8.2500
    398   58.3833
    399 1
    400      0.000
    401    100.000
    402 4000
    403 1
    404 1
    405 1
    406 1
    407 1
    408 rel_bj
    409 
    410 20040708 080000
    411 20040709 080000
    412    8.2500
    413   58.3833
    414    8.2500
    415   58.3833
    416 1
    417      0.000
    418    100.000
    419 4000
    420 1
    421 1
    422 1
    423 1
    424 1
    425 rel_bj
    426 
    427 20040715 080000
    428 20040716 080000
    429    8.2500
    430   58.3833
    431    8.2500
    432   58.3833
    433 1
    434      0.000
    435    100.000
    436 4000
    437 1
    438 1
    439 1
    440 1
    441 1
    442 rel_bj
    443 
    444 20040722 080000
    445 20040723 080000
    446    8.2500
    447   58.3833
    448    8.2500
    449   58.3833
    450 1
    451      0.000
    452    100.000
    453 4000
    454 1
    455 1
    456 1
    457 1
    458 1
    459 rel_bj
    460 
    461 20040729 080000
    462 20040730 080000
    463    8.2500
    464   58.3833
    465    8.2500
    466   58.3833
    467 1
    468      0.000
    469    100.000
    470 4000
    471 1
    472 1
    473 1
    474 1
    475 1
    476 rel_bj
    477 
    478 20040805 080000
    479 20040806 080000
    480    8.2500
    481   58.3833
    482    8.2500
    483   58.3833
    484 1
    485      0.000
    486    100.000
    487 4000
    488 1
    489 1
    490 1
    491 1
    492 1
    493 rel_bj
    494 
    495 20040812 080000
    496 20040813 080000
    497    8.2500
    498   58.3833
    499    8.2500
    500   58.3833
    501 1
    502      0.000
    503    100.000
    504 4000
    505 1
    506 1
    507 1
    508 1
    509 1
    510 rel_bj
    511 
    512 20040819 080000
    513 20040820 080000
    514    8.2500
    515   58.3833
    516    8.2500
    517   58.3833
    518 1
    519      0.000
    520    100.000
    521 4000
    522 1
    523 1
    524 1
    525 1
    526 1
    527 rel_bj
    528 
    529 20040826 080000
    530 20040827 080000
    531    8.2500
    532   58.3833
    533    8.2500
    534   58.3833
    535 1
    536      0.000
    537    100.000
    538 4000
    539 1
    540 1
    541 1
    542 1
    543 1
    544 rel_bj
    545 
    546 20040902 080000
    547 20040903 080000
    548    8.2500
    549   58.3833
    550    8.2500
    551   58.3833
    552 1
    553      0.000
    554    100.000
    555 4000
    556 1
    557 1
    558 1
    559 1
    560 1
    561 rel_bj
    562 
    563 20040909 080000
    564 20040910 080000
    565    8.2500
    566   58.3833
    567    8.2500
    568   58.3833
    569 1
    570      0.000
    571    100.000
    572 4000
    573 1
    574 1
    575 1
    576 1
    577 1
    578 rel_bj
    579 
    580 20040916 080000
    581 20040917 080000
    582    8.2500
    583   58.3833
    584    8.2500
    585   58.3833
    586 1
    587      0.000
    588    100.000
    589 4000
    590 1
    591 1
    592 1
    593 1
    594 1
    595 rel_bj
    596 
    597 20040923 080000
    598 20040924 080000
    599    8.2500
    600   58.3833
    601    8.2500
    602   58.3833
    603 1
    604      0.000
    605    100.000
    606 4000
    607 1
    608 1
    609 1
    610 1
    611 1
    612 rel_bj
    613 
    614 20040930 080000
    615 20041001 080000
    616    8.2500
    617   58.3833
    618    8.2500
    619   58.3833
    620 1
    621      0.000
    622    100.000
    623 4000
    624 1
    625 1
    626 1
    627 1
    628 1
    629 rel_bj
    630 
    631 20041007 080000
    632 20041008 080000
    633    8.2500
    634   58.3833
    635    8.2500
    636   58.3833
    637 1
    638      0.000
    639    100.000
    640 4000
    641 1
    642 1
    643 1
    644 1
    645 1
    646 rel_bj
    647 
    648 20041014 080000
    649 20041015 080000
    650    8.2500
    651   58.3833
    652    8.2500
    653   58.3833
    654 1
    655      0.000
    656    100.000
    657 4000
    658 1
    659 1
    660 1
    661 1
    662 1
    663 rel_bj
    664 
    665 20041021 080000
    666 20041022 080000
    667    8.2500
    668   58.3833
    669    8.2500
    670   58.3833
    671 1
    672      0.000
    673    100.000
    674 4000
    675 1
    676 1
    677 1
    678 1
    679 1
    680 rel_bj
    681 
    682 20041028 080000
    683 20041029 080000
    684    8.2500
    685   58.3833
    686    8.2500
    687   58.3833
    688 1
    689      0.000
    690    100.000
    691 4000
    692 1
    693 1
    694 1
    695 1
    696 1
    697 rel_bj
    698 
    699 20041104 080000
    700 20041105 080000
    701    8.2500
    702   58.3833
    703    8.2500
    704   58.3833
    705 1
    706      0.000
    707    100.000
    708 4000
    709 1
    710 1
    711 1
    712 1
    713 1
    714 rel_bj
    715 
    716 20041111 080000
    717 20041112 080000
    718    8.2500
    719   58.3833
    720    8.2500
    721   58.3833
    722 1
    723      0.000
    724    100.000
    725 4000
    726 1
    727 1
    728 1
    729 1
    730 1
    731 rel_bj
    732 
    733 20041118 080000
    734 20041119 080000
    735    8.2500
    736   58.3833
    737    8.2500
    738   58.3833
    739 1
    740      0.000
    741    100.000
    742 4000
    743 1
    744 1
    745 1
    746 1
    747 1
    748 rel_bj
    749 
    750 20041124 080000
    751 20041125 080000
    752    8.2500
    753   58.3833
    754    8.2500
    755   58.3833
    756 1
    757      0.000
    758    100.000
    759 4000
    760 1
    761 1
    762 1
    763 1
    764 1
    765 rel_bj
    766 
    767 20041125 080000
    768 20041126 080000
    769    8.2500
    770   58.3833
    771    8.2500
    772   58.3833
    773 1
    774      0.000
    775    100.000
    776 4000
    777 1
    778 1
    779 1
    780 1
    781 1
    782 rel_bj
    783 
    784 20041208 080000
    785 20041209 080000
    786    8.2500
    787   58.3833
    788    8.2500
    789   58.3833
    790 1
    791      0.000
    792    100.000
    793 4000
    794 1
    795 1
    796 1
    797 1
    798 1
    799 rel_bj
    800 
    801 20041209 080000
    802 20041210 080000
    803    8.2500
    804   58.3833
    805    8.2500
    806   58.3833
    807 1
    808      0.000
    809    100.000
    810 4000
    811 1
    812 1
    813 1
    814 1
    815 1
    816 rel_bj
    817 
    818 20041216 080000
    819 20041217 080000
    820    8.2500
    821   58.3833
    822    8.2500
    823   58.3833
    824 1
    825      0.000
    826    100.000
    827 4000
    828 1
    829 1
    830 1
    831 1
    832 1
    833 rel_bj
    834 
    835 20041223 080000
    836 20041224 080000
    837    8.2500
    838   58.3833
    839    8.2500
    840   58.3833
    841 1
    842      0.000
    843    100.000
    844 4000
    845 1
    846 1
    847 1
    848 1
    849 1
    850 rel_bj
    851 
    852 20041230 080000
    853 20041231 080000
    854    8.2500
    855   58.3833
    856    8.2500
    857   58.3833
    858 1
    859      0.000
    860    100.000
    861 4000
    862 1
    863 1
    864 1
    865 1
    866 1
    867 rel_bj
    868 
     52RELEASE_TEST1
     53________________________________________   character*40 comment
     54+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  • trunk/options/SPECIES/SPECIES_001

    r4 r24  
    99-9.9E-09           Wet deposition - A
    1010                   Wet deposition - B
     11-9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     12-9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     13-9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
     14-9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_002

    r4 r24  
    99-9.9E-09           Wet deposition - A
    1010                   Wet deposition - B
     11-9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     12-9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     13-9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
     14-9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 1.5               Dry deposition (gases) - D
    12161.0e-02            Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_003

    r4 r24  
    99 8.0E-06           Wet deposition - A
    10100.62               Wet deposition - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.1                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 1.2               Dry deposition (gases) - D
    12163.0e-03            Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_004

    r4 r24  
    99 1.0E-05           Wet deposition - A
    10100.62               Wet deposition - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.1                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 1.6               Dry deposition (gases) - D
    12161.0e-02            Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_005

    r4 r24  
    99 5.0E-05           Wet deposition - A
    10100.62               Wet deposition - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.1                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 1.9               Dry deposition (gases) - D
    12161.0e+14            Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_006

    r4 r24  
    99-9.9E-09           Wet deposition - A
    1010                   Wet deposition - B
     11-9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     12-9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     13-9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
     14-9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 1.6               Dry deposition (gases) - D
    12161.0e+05            Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_007

    r4 r24  
    99 1.0E-04           Wet deposition - A
    10100.62               Wet deposition - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.1                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 1.4               Dry deposition (gases) - D
    12161.0e+05            Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_008

    r4 r24  
    99-9.9E-09           Wet deposition - A
    10100.62               Wet deposition - B
     11-9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     12-9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     13-9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
     14-9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 2.0               Dry deposition (gases) - D
    12161.0e+05            Dry deposition (gases) - Henrys const.
    … …  
    394316       1.463       1.162
    404417       1.426       1.141
    41 18       1.326       1.117 
     4518       1.326       1.117
    424619       1.225       1.109
    434720       1.082       1.095

  • trunk/options/SPECIES/SPECIES_009

    r4 r24  
    99-9.9E-09           Wet deposition - A
    1010                   Wet deposition - B
     11-9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     12-9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     13-9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
     14-9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 1.3               Dry deposition (gases) - D
    12166.0e+03            Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_010

    r4 r24  
    99-9.9E-09           Wet deposition - A
    1010                   Wet deposition - B
     11-9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     12-9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     13-9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
     14-9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 2.6               Dry deposition (gases) - D
    12163.6e+00            Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_011

    r4 r24  
    99 9.9E-05           Wet deposition - A
    10100.62               Wet deposition - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.1                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 1.1               Dry deposition (gases) - D
    12162.0e+14            Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_012

    r4 r24  
    99 5.0E-06           Wet deposition - A
    10100.62               Wet deposition - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.1                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216 1.0E-09           Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_013

    r4 r24  
    99 5.0E-06           Wet deposition - A
    10100.62               Wet deposition - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.1                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_014

    r4 r24  
    99 8.0E-05           Wet deposition - A
    10100.62               Wet deposition - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.1                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115 2.7               Dry deposition (gases) - D
    12161.0e+05            Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_015

    r4 r24  
    99 1.0E-04           Wet deposition - A
    10100.80               Wet deposition - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.1                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_016

    r4 r24  
    99 1.0E-04           Wet deposition - A
    10100.80               Wet deposition - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.1                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_017

    r4 r24  
    99 1.0E-04           Wet deposition - A
    10100.80               Wet deposition - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.1                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_018

    r4 r24  
    99 1.0E-04           Wet deposition - A
    10100.80               Wet deposition - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.1                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_019

    r4 r24  
    99-9.9E-09           Wet deposition - A
    1010                   Wet deposition - B
     11-9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     12-9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     13-9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
     14-9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_020

    r4 r24  
    99 1.0E-04           Wet deposition - A
    10100.80               Wet deposition - B
     112.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     120.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     130.1                In-cloud scavenging - Ci (S_i=Ci/cl)
     141.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_021

    r4 r24  
    66****************************************************************************
    77Xe-133             Tracer name
    8 198720.0           Species half life
     8453168.0           Species half life
    99-9.9E-09           Wet deposition - A
    1010                   Wet deposition - B
     11-9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     12-9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     13-9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
     14-9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_022

    r4 r24  
    99-9.9E-09           Wet deposition - A
    1010                   Wet deposition - B
     11-9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     12-9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     13-9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
     14-9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_024

    r4 r24  
    99-9.9E-09           Wet deposition - A
    1010                   Wet deposition - B
     11-9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
     12-9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
     13-9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
     14-9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
    1115-9.9               Dry deposition (gases) - D
    1216                   Dry deposition (gases) - Henrys const.

  • trunk/options/SPECIES/SPECIES_025

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

  • trunk/src/FLEXPART.f90

    r20 r24  
    6161
    6262  !
    63   flexversion='Version 9.1.8  (2013-12-08)'
     63  flexversion='Version 9.2 beta (2014-05-23)'
    6464  !verbosity=0
    6565  ! Read the pathnames where input/output files are stored
    … …  
    111111  ! Print the GPL License statement
    112112  !*******************************************************
    113   print*,'Welcome to FLEXPART', trim(flexversion)
     113  print*,'Welcome to FLEXPART ', trim(flexversion)
    114114  print*,'FLEXPART is free software released under the GNU Genera'// &
    115115       'l Public License.'

  • trunk/src/calcpv.f90

    r4 r24  
    2121
    2222subroutine calcpv(n,uuh,vvh,pvh)
    23   !                  i  i   i   o
     23  !               i  i   i   o
    2424  !*****************************************************************************
    2525  !                                                                            *
    … …  
    4949  integer :: nlck
    5050  real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
    51   real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
    52   real :: pvavr,ppml(nuvzmax)
     51  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt
     52  real :: pvavr,ppml(0:nxmax-1,0:nymax-1,nuvzmax),ppmk(0:nxmax-1,0:nymax-1,nuvzmax)
    5353  real :: thup,thdn
    5454  real,parameter :: eps=1.e-5, p0=101325
    … …  
    6262  ! Loop over entire grid
    6363  !**********************
     64  do kl=1,nuvz
     65    do jy=0,nymin1
     66      do ix=0,nxmin1
     67         ppml(ix,jy,kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
     68      enddo
     69    enddo
     70  enddo
     71  ppmk=(100000./ppml)**kappa
     72
    6473  do jy=0,nymin1
    6574    if (sglobal.and.jy.eq.0) goto 10
    … …  
    106115      end if
    107116      jux=jumpx
    108   ! Precalculate pressure values for efficiency
    109       do kl=1,nuvz
    110         ppml(kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
    111       end do
    112117  !
    113118  ! Loop over the vertical
    … …  
    115120
    116121      do kl=1,nuvz
    117         ppmk=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
    118         theta=tth(ix,jy,kl,n)*(100000./ppmk)**kappa
     122        theta=tth(ix,jy,kl,n)*ppmk(ix,jy,kl)
    119123        klvrp=kl+1
    120124        klvrm=kl-1
    … …  
    125129        if (klvrp.gt.nuvz) klvrp=nuvz
    126130        if (klvrm.lt.1) klvrm=1
    127         ppmk=akz(klvrp)+bkz(klvrp)*ps(ix,jy,1,n)
    128         thetap=tth(ix,jy,klvrp,n)*(100000./ppmk)**kappa
    129         ppmk=akz(klvrm)+bkz(klvrm)*ps(ix,jy,1,n)
    130         thetam=tth(ix,jy,klvrm,n)*(100000./ppmk)**kappa
    131         dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm))
     131        thetap=tth(ix,jy,klvrp,n)*ppmk(ix,jy,klvrp)
     132        thetam=tth(ix,jy,klvrm,n)*ppmk(ix,jy,klvrm)
     133        dthetadp=(thetap-thetam)/(ppml(ix,jy,klvrp)-ppml(ix,jy,klvrm))
    132134
    133135  ! Compute vertical position at pot. temperature surface on subgrid
    … …  
    156158          kch=kch+1
    157159          k=kup
    158           ppmk=akz(k)+bkz(k)*ps(ivr,jy,1,n)
    159           thdn=tth(ivr,jy,k,n)*(100000./ppmk)**kappa
    160           ppmk=akz(k+1)+bkz(k+1)*ps(ivr,jy,1,n)
    161           thup=tth(ivr,jy,k+1,n)*(100000./ppmk)**kappa
    162 
    163 
    164       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    165            ((thdn.le.theta).and.(thup.ge.theta))) then
    166               dt1=abs(theta-thdn)
    167               dt2=abs(theta-thup)
    168               dt=dt1+dt2
    169               if (dt.lt.eps) then   ! Avoid division by zero error
    170                 dt1=0.5             ! G.W., 10.4.1996
    171                 dt2=0.5
    172                 dt=1.0
    173               endif
    174           vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
    175               goto 20
     160          thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k)
     161          thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1)
     162
     163
     164          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
     165          ((thdn.le.theta).and.(thup.ge.theta))) then
     166            dt1=abs(theta-thdn)
     167            dt2=abs(theta-thup)
     168            dt=dt1+dt2
     169            if (dt.lt.eps) then   ! Avoid division by zero error
     170              dt1=0.5             ! G.W., 10.4.1996
     171              dt2=0.5
     172              dt=1.0
    176173            endif
     174            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
     175            goto 20
     176          endif
    17717741        continue
    178178  ! Downward branch
    … …  
    181181          kch=kch+1
    182182          k=kdn
    183           ppmk=akz(k)+bkz(k)*ps(ivr,jy,1,n)
    184           thdn=tth(ivr,jy,k,n)*(100000./ppmk)**kappa
    185           ppmk=akz(k+1)+bkz(k+1)*ps(ivr,jy,1,n)
    186           thup=tth(ivr,jy,k+1,n)*(100000./ppmk)**kappa
    187 
    188       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    189            ((thdn.le.theta).and.(thup.ge.theta))) then
    190               dt1=abs(theta-thdn)
    191               dt2=abs(theta-thup)
    192               dt=dt1+dt2
    193               if (dt.lt.eps) then   ! Avoid division by zero error
    194                 dt1=0.5             ! G.W., 10.4.1996
    195                 dt2=0.5
    196                 dt=1.0
    197               endif
    198           vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
    199               goto 20
     183          thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k)
     184          thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1)
     185
     186          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
     187          ((thdn.le.theta).and.(thup.ge.theta))) then
     188            dt1=abs(theta-thdn)
     189            dt2=abs(theta-thup)
     190            dt=dt1+dt2
     191            if (dt.lt.eps) then   ! Avoid division by zero error
     192              dt1=0.5             ! G.W., 10.4.1996
     193              dt2=0.5
     194              dt=1.0
    200195            endif
    201             goto 40
     196            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
     197            goto 20
     198          endif
     199          goto 40
    202200  ! This section used when no values were found
    20320121      continue
    … …  
    236234          kch=kch+1
    237235          k=kup
    238           ppmk=akz(k)+bkz(k)*ps(ix,j,1,n)
    239           thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa
    240           ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n)
    241           thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa
    242       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    243            ((thdn.le.theta).and.(thup.ge.theta))) then
    244               dt1=abs(theta-thdn)
    245               dt2=abs(theta-thup)
    246               dt=dt1+dt2
    247               if (dt.lt.eps) then   ! Avoid division by zero error
    248                 dt1=0.5             ! G.W., 10.4.1996
    249                 dt2=0.5
    250                 dt=1.0
    251               endif
    252               uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
    253               goto 50
     236          thdn=tth(ix,j,k,n)*ppmk(ix,j,k)
     237          thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1)
     238          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
     239          ((thdn.le.theta).and.(thup.ge.theta))) then
     240            dt1=abs(theta-thdn)
     241            dt2=abs(theta-thup)
     242            dt=dt1+dt2
     243            if (dt.lt.eps) then   ! Avoid division by zero error
     244              dt1=0.5             ! G.W., 10.4.1996
     245              dt2=0.5
     246              dt=1.0
    254247            endif
     248            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
     249            goto 50
     250          endif
    25525171        continue
    256252  ! Downward branch
    … …  
    259255          kch=kch+1
    260256          k=kdn
    261           ppmk=akz(k)+bkz(k)*ps(ix,j,1,n)
    262           thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa
    263           ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n)
    264           thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa
    265       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    266            ((thdn.le.theta).and.(thup.ge.theta))) then
    267               dt1=abs(theta-thdn)
    268               dt2=abs(theta-thup)
    269               dt=dt1+dt2
    270               if (dt.lt.eps) then   ! Avoid division by zero error
    271                 dt1=0.5             ! G.W., 10.4.1996
    272                 dt2=0.5
    273                 dt=1.0
    274               endif
    275               uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
    276               goto 50
     257          thdn=tth(ix,j,k,n)*ppmk(ix,j,k)
     258          thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1)
     259          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
     260          ((thdn.le.theta).and.(thup.ge.theta))) then
     261            dt1=abs(theta-thdn)
     262            dt2=abs(theta-thup)
     263            dt=dt1+dt2
     264            if (dt.lt.eps) then   ! Avoid division by zero error
     265              dt1=0.5             ! G.W., 10.4.1996
     266              dt2=0.5
     267              dt=1.0
    277268            endif
    278             goto 70
     269            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
     270            goto 50
     271          endif
     272          goto 70
    279273  ! This section used when no values were found
    28027451      continue
    281275  ! Must use uu at current level and lat. juy becomes smaller by 1
    282         uy(jj)=uuh(ix,jy,kl)
    283         juy=juy-1
     276          uy(jj)=uuh(ix,jy,kl)
     277          juy=juy-1
    284278  ! Otherwise OK
    28527950        continue

  • trunk/src/calcpv_nests.f90

    r4 r24  
    2121
    2222subroutine calcpv_nests(l,n,uuhn,vvhn,pvhn)
    23   !                        i i  i    i    o
     23  !                     i i  i    i    o
    2424  !*****************************************************************************
    2525  !                                                                            *
    … …  
    4848  integer :: nlck,l
    4949  real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
    50   real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
    51   real :: ppml(nuvzmax)
     50  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt
     51  real :: ppml(0:nxmaxn-1,0:nymaxn-1,nuvzmax),ppmk(0:nxmaxn-1,0:nymaxn-1,nuvzmax)
    5252  real :: thup,thdn
    5353  real,parameter :: eps=1.e-5,p0=101325
    … …  
    6161  ! Loop over entire grid
    6262  !**********************
     63  do kl=1,nuvz
     64    do jy=0,nyn(l)-1
     65      do ix=0,nxn(l)-1
     66         ppml(ix,jy,kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
     67      enddo
     68    enddo
     69  enddo
     70  ppmk=(100000./ppml)**kappa
     71
    6372  do jy=0,nyn(l)-1
    6473    phi = (ylat0n(l) + jy * dyn(l)) * pi / 180.
    … …  
    8897      if (ix.eq.0.or.ix.eq.nxn(l)-1) jumpx=1
    8998      jux=jumpx
    90   ! Precalculate pressure values for efficiency
    91       do kl=1,nuvz
    92         ppml(kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
    93       end do
    9499  !
    95100  ! Loop over the vertical
    … …  
    97102
    98103      do kl=1,nuvz
    99         ppmk=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
    100         theta=tthn(ix,jy,kl,n,l)*(100000./ppmk)**kappa
     104        theta=tthn(ix,jy,kl,n,l)*ppmk(ix,jy,kl)
    101105        klvrp=kl+1
    102106        klvrm=kl-1
    … …  
    107111        if (klvrp.gt.nuvz) klvrp=nuvz
    108112        if (klvrm.lt.1) klvrm=1
    109         ppmk=akz(klvrp)+bkz(klvrp)*psn(ix,jy,1,n,l)
    110         thetap=tthn(ix,jy,klvrp,n,l)*(100000./ppmk)**kappa
    111         ppmk=akz(klvrm)+bkz(klvrm)*psn(ix,jy,1,n,l)
    112         thetam=tthn(ix,jy,klvrm,n,l)*(100000./ppmk)**kappa
    113         dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm))
     113        thetap=tthn(ix,jy,klvrp,n,l)*ppmk(ix,jy,klvrp)
     114        thetam=tthn(ix,jy,klvrm,n,l)*ppmk(ix,jy,klvrm)
     115        dthetadp=(thetap-thetam)/(ppml(ix,jy,klvrp)-ppml(ix,jy,klvrm))
    114116
    115117  ! Compute vertical position at pot. temperature surface on subgrid
    … …  
    134136          kch=kch+1
    135137          k=kup
    136           ppmk=akz(k)+bkz(k)*psn(ivr,jy,1,n,l)
    137           thdn=tthn(ivr,jy,k,n,l)*(100000./ppmk)**kappa
    138           ppmk=akz(k+1)+bkz(k+1)*psn(ivr,jy,1,n,l)
    139           thup=tthn(ivr,jy,k+1,n,l)*(100000./ppmk)**kappa
     138          thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k)
     139          thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1)
    140140
    141141      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    … …  
    158158          kch=kch+1
    159159          k=kdn
    160           ppmk=akz(k)+bkz(k)*psn(ivr,jy,1,n,l)
    161           thdn=tthn(ivr,jy,k,n,l)*(100000./ppmk)**kappa
    162           ppmk=akz(k+1)+bkz(k+1)*psn(ivr,jy,1,n,l)
    163           thup=tthn(ivr,jy,k+1,n,l)*(100000./ppmk)**kappa
     160          thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k)
     161          thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1)
    164162      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    165163           ((thdn.le.theta).and.(thup.ge.theta))) then
    … …  
    212210          kch=kch+1
    213211          k=kup
    214           ppmk=akz(k)+bkz(k)*psn(ix,j,1,n,l)
    215           thdn=tthn(ix,j,k,n,l)*(100000./ppmk)**kappa
    216           ppmk=akz(k+1)+bkz(k+1)*psn(ix,j,1,n,l)
    217           thup=tthn(ix,j,k+1,n,l)*(100000./ppmk)**kappa
     212          thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k)
     213          thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1)
    218214      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    219215           ((thdn.le.theta).and.(thup.ge.theta))) then
    … …  
    235231          kch=kch+1
    236232          k=kdn
    237           ppmk=akz(k)+bkz(k)*psn(ix,j,1,n,l)
    238           thdn=tthn(ix,j,k,n,l)*(100000./ppmk)**kappa
    239           ppmk=akz(k+1)+bkz(k+1)*psn(ix,j,1,n,l)
    240           thup=tthn(ix,j,k+1,n,l)*(100000./ppmk)**kappa
     233          thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k)
     234          thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1)
    241235      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    242236           ((thdn.le.theta).and.(thup.ge.theta))) then

  • trunk/src/com_mod.f90

    r20 r24  
    340340  real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2)
    341341  real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2)
    342   !scavenging NIK, PS
    343342  integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2)
    344343  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)
    347344
    348345
    … …  
    359356  !      rainout  conv/lsp dominated  2/3
    360357  !      washout  conv/lsp dominated  4/5
    361 ! PS 2013
    362 !c icloudbot (m)        cloud bottom height
    363 !c icloudthck (m)       cloud thickness     
    364358
    365359  ! pplev for the GFS version

  • trunk/src/gridcheck.f90

    r20 r24  
    8484  real(kind=4) :: xaux1,xaux2,yaux1,yaux2
    8585  real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
    86   integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
     86  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parID
    8787  !HSO  end
    8888  integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
    89   real :: sizesouth,sizenorth,xauxa,pint
     89  real :: sizesouth,sizenorth,xauxa,pint,conversion_factor
    9090
    9191  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
    … …  
    9898
    9999  integer :: isec1(56),isec2(22+nxmax+nymax)
    100   !real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
    101   real(kind=4) :: zsec2(184),zsec4(jpunp)
     100  real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
    102101  character(len=1) :: opt
    103102
    … …  
    172171  call grib_check(iret,gribFunction,gribErrorMsg)
    173172  call grib_get_int(igrib,'level',valSurf,iret)
     173  call grib_check(iret,gribFunction,gribErrorMsg)
     174  call grib_get_int(igrib,'paramId',parId,iret)
    174175  call grib_check(iret,gribFunction,gribErrorMsg)
    175176
    … …  
    193194  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
    194195    isec1(6)=135         ! indicatorOfParameter
     196  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
     197    isec1(6)=135         ! indicatorOfParameter
    195198  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    196199    isec1(6)=151         ! indicatorOfParameter
    … …  
    205208  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    206209    isec1(6)=141         ! indicatorOfParameter
    207   elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
     210  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
    208211    isec1(6)=164         ! indicatorOfParameter
    209   elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
     212  elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
    210213    isec1(6)=142         ! indicatorOfParameter
    211214  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    … …  
    215218  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    216219    isec1(6)=176         ! indicatorOfParameter
    217   elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
     220  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
    218221    isec1(6)=180         ! indicatorOfParameter
    219   elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
     222  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
    220223    isec1(6)=181         ! indicatorOfParameter
    221224  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    222225    isec1(6)=129         ! indicatorOfParameter
    223   elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
     226  elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
    224227    isec1(6)=160         ! indicatorOfParameter
    225228  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    … …  
    230233         parCat,parNum,typSurf
    231234  endif
     235  if(parId .ne. isec1(6) .and. parId .ne. 77) then
     236    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
     237!    stop
     238  endif
    232239
    233240  endif
    … …  
    265272
    266273  !HSO  get the second part of the grid dimensions only from GRiB1 messages
    267   if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
     274  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
    268275    call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
    269276         xaux2in,iret)
    … …  
    291298    dyconst=180./(dy*r_earth*pi)
    292299    gotGrid=1
    293 
    294300  ! Check whether fields are global
    295301  ! If they contain the poles, specify polar stereographic map
    … …  
    443449  write(*,*)
    444450  write(*,'(a)') 'Mother domain:'
    445   write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
     451  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
    446452       xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
    447   write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
     453  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range: ', &
    448454       ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
    449455  write(*,*)
    … …  
    467473    k=nlev_ec+1+numskip+i
    468474    akm(nwz-i+1)=zsec2(j)
     475  !   write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)
    469476    bkm(nwz-i+1)=zsec2(k)
    470477  end do
    … …  
    553560
    554561end subroutine gridcheck
     562

  • trunk/src/gridcheck_nests.f90

    r4 r24  
    4848  integer :: igrib
    4949  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
     50  integer :: parID !added by mc for making it consistent with new gridcheck.f90
    5051  integer :: gotGrib
    5152  !HSO  end
    … …  
    5556  real(kind=4) :: xaux1,xaux2,yaux1,yaux2
    5657  real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
     58  real :: conversion_factor !added by mc to make it consistent with new gridchek.f90
    5759
    5860  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
    … …  
    150152  call grib_get_int(igrib,'level',valSurf,iret)
    151153  call grib_check(iret,gribFunction,gribErrorMsg)
     154  call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consisitent with new grid_check.f90
     155  call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consisitent with new  grid_check.f90
    152156
    153157  !print*,discipl,parCat,parNum,typSurf,valSurf
    … …  
    170174  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
    171175    isec1(6)=135         ! indicatorOfParameter
     176  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added bymc to make it consistent with new gridchek.f90
     177    isec1(6)=135         ! indicatorOfParameter    !                     ! " " " " " " " " " " " " " " " " " " " " " " " "  " " "
    172178  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    173179    isec1(6)=151         ! indicatorOfParameter
    … …  
    182188  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    183189    isec1(6)=141         ! indicatorOfParameter
    184   elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
     190  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consistent with new gridchek.f90
    185191    isec1(6)=164         ! indicatorOfParameter
    186   elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
     192 elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consistent with new gridchek.f90
    187193    isec1(6)=142         ! indicatorOfParameter
    188194  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    … …  
    192198  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    193199    isec1(6)=176         ! indicatorOfParameter
    194   elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
     200  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consistent with new gridchek.f90
    195201    isec1(6)=180         ! indicatorOfParameter
    196   elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
     202  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consistent with new gridchek.f90
    197203    isec1(6)=181         ! indicatorOfParameter
    198204  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    199205    isec1(6)=129         ! indicatorOfParameter
    200   elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
     206  elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consistent with new gridchek.f90
    201207    isec1(6)=160         ! indicatorOfParameter
    202208  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    … …  
    207213         parCat,parNum,typSurf
    208214  endif
     215  if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consistent with new gridchek.f90
     216    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
     217!    stop
     218  endif
    209219
    210220  endif
    … …  
    237247
    238248  !HSO  get the second part of the grid dimensions only from GRiB1 messages
    239   if ((gribVer.eq.1).and.(gotGrib.eq.0)) then
    240     call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
     249  if (isec1(6) .eq. 167 .and. (gotGrib.eq.0)) then !added by mc to make it consistent with new gridchek.f90 note that gotGrid must be changed in gotGrib!!
     250    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', & !comment by mc: note that this was in the (if (ifield.eq.1) ..end above in gridchek.f90 see line 257
    241251         xaux1in,iret)
    242252    call grib_check(iret,gribFunction,gribErrorMsg)
    … …  
    254264    yaux1=yaux1in
    255265    yaux2=yaux2in
    256     if(xaux1.gt.180) xaux1=xaux1-360.0
    257     if(xaux2.gt.180) xaux2=xaux2-360.0
    258     if(xaux1.lt.-180) xaux1=xaux1+360.0
    259     if(xaux2.lt.-180) xaux2=xaux2+360.0
    260     if (xaux2.lt.xaux1) xaux2=xaux2+360.
     266    if(xaux1.gt.180.) xaux1=xaux1-360.0
     267    if(xaux2.gt.180.) xaux2=xaux2-360.0
     268    if(xaux1.lt.-180.) xaux1=xaux1+360.0
     269    if(xaux2.lt.-180.) xaux2=xaux2+360.0
     270    if (xaux2.lt.xaux1) xaux2=xaux2+360.0
    261271    xlon0n(l)=xaux1
    262272    ylat0n(l)=yaux1
    263273    dxn(l)=(xaux2-xaux1)/real(nxn(l)-1)
    264274    dyn(l)=(yaux2-yaux1)/real(nyn(l)-1)
    265     gotGrib=1
     275    gotGrib=1 !commetn by mc note tahthere gotGRIB is used instead of gotGrid!!!
    266276  endif ! ifield.eq.1
    267277
    … …  
    341351
    342352  write(*,'(a,i2)') 'Nested domain #: ',l
    343   write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
     353  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
    344354       xlon0n(l),' to ',xlon0n(l)+(nxn(l)-1)*dxn(l), &
    345355       '   Grid distance: ',dxn(l)
    346   write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
     356  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range: ', &
    347357       ylat0n(l),' to ',ylat0n(l)+(nyn(l)-1)*dyn(l), &
    348358       '   Grid distance: ',dyn(l)

  • trunk/src/interpol_rain.f90

    r20 r24  
    2020!**********************************************************************
    2121
    22 !subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
    23 !       ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
    24 !  !                          i   i   i    i    i     i   i
    25 !  !i    i    i  i    i     i      i      i     o     o     o
    26 
    27       subroutine interpol_rain(yy1,yy2,yy3,iy1,iy2,nxmax,nymax,nzmax,nx, &
    28      ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, &
    29      intiy1,intiy2,icmv)
    30 !                               i   i   i    i    i     i   i
    31 !     i    i    i  i    i     i      i      i     o     o     o
    32 
     22subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
     23       ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
     24  !                          i   i   i    i    i     i   i
     25  !i    i    i  i    i     i      i      i     o     o     o
    3326  !****************************************************************************
    3427  !                                                                           *
    … …  
    4740  !                                                                           *
    4841  !     30 August 1996                                                        *
    49   !
    50   !*  Petra Seibert, 2011/2012:
    51   !*  Add interpolation of cloud bottom and cloud thickness
    52   !*  for fix to SE's new wet scavenging scheme  *
     42  !                                                                           *
    5343  !****************************************************************************
    5444  !                                                                           *
    … …  
    8070
    8171  integer :: nx,ny,nxmax,nymax,nzmax,memind(2),m,ix,jy,ixp,jyp
    82   !integer :: itime,itime1,itime2,level,indexh
    83   integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4
    84   integer :: intiy1,intiy2,ipsum,icmv 
     72  integer :: itime,itime1,itime2,level,indexh
    8573  real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2)
    8674  real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2)
    8775  real :: yy3(0:nxmax-1,0:nymax-1,nzmax,2)
    88   integer iy1(0:nxmax-1,0:nymax-1,2),iy2(0:nxmax-1,0:nymax-1,2)
    89   real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),yi1(2),yi2(2)
    90   !real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
    91   real :: xt,yt,yint1,yint2,yint3,yint4,p1,p2,p3,p4 
    92   !real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
     76  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
     77  real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
    9378
    9479
    … …  
    142127         + p3*yy3(ix ,jyp,level,indexh) &
    143128         + p4*yy3(ixp,jyp,level,indexh)
    144 
    145 !CPS clouds:
    146         ip1=1
    147         ip2=1
    148         ip3=1
    149         ip4=1
    150         if (iy1(ix ,jy ,indexh) .eq. icmv) ip1=0
    151         if (iy1(ixp,jy ,indexh) .eq. icmv) ip2=0
    152         if (iy1(ix ,jyp,indexh) .eq. icmv) ip3=0
    153         if (iy1(ixp,jyp,indexh) .eq. icmv) ip4=0
    154         ipsum= ip1+ip2+ip3+ip4
    155         if (ipsum .eq. 0) then
    156           yi1(m)=icmv
    157         else
    158           yi1(m)=(ip1*p1*iy1(ix ,jy ,indexh)    &
    159              + ip2*p2*iy1(ixp,jy ,indexh)       &
    160              + ip3*p3*iy1(ix ,jyp,indexh)       &
    161              + ip4*p4*iy1(ixp,jyp,indexh))/ipsum
    162         endif
    163        
    164         ip1=1
    165         ip2=1
    166         ip3=1
    167         ip4=1
    168         if (iy2(ix ,jy ,indexh) .eq. icmv) ip1=0
    169         if (iy2(ixp,jy ,indexh) .eq. icmv) ip2=0
    170         if (iy2(ix ,jyp,indexh) .eq. icmv) ip3=0
    171         if (iy2(ixp,jyp,indexh) .eq. icmv) ip4=0
    172         ipsum= ip1+ip2+ip3+ip4
    173         if (ipsum .eq. 0) then
    174           yi2(m)=icmv
    175         else
    176           yi2(m)=(ip1*p1*iy2(ix ,jy ,indexh)  &
    177              + ip2*p2*iy2(ixp,jy ,indexh)     &
    178              + ip3*p3*iy2(ix ,jyp,indexh)     &
    179              + ip4*p4*iy2(ixp,jyp,indexh))/ipsum
    180         endif
    181 !CPS end clouds
    182 
    183129  end do
    184130
    … …  
    187133  ! 2.) Temporal interpolation (linear)
    188134  !************************************
    189 
    190       if (abs(itime) .lt. abs(itime1)) then
    191         print*,'interpol_rain.f90'
    192         print*,itime,itime1,itime2
    193         stop 'ITIME PROBLEM'
    194       endif
    195 
    196135
    197136  dt1=real(itime-itime1)
    … …  
    204143
    205144
    206 !PS clouds:
    207       intiy1=(yi1(1)*dt2 + yi1(2)*dt1)/dt
    208       if (yi1(1) .eq. float(icmv)) intiy1=yi1(2)
    209       if (yi1(2) .eq. float(icmv)) intiy1=yi1(1)
    210 
    211       intiy2=(yi2(1)*dt2 + yi2(2)*dt1)/dt
    212       if (yi2(1) .eq. float(icmv)) intiy2=yi2(2) 
    213       if (yi2(2) .eq. float(icmv)) intiy2=yi2(1)
    214      
    215       if (intiy1 .ne. icmv .and. intiy2 .ne. icmv) then
    216         intiy2 = intiy1 + intiy2 ! convert cloud thickness to cloud top
    217       else
    218         intiy1=icmv
    219         intiy2=icmv
    220       endif
    221 !PS end clouds
    222 
    223145end subroutine interpol_rain

  • trunk/src/makefile

    r20 r24  
    11SHELL = /bin/bash
    2 TARGET = local
    3 WINDS=ecmwf
    4 #WINDS=gfs
    5 #WINDS=fnl
     2MAIN = FP_ecmwf_gfortran
    63
    74FC       = gfortran
     5INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
     6LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
     7LIBPATH2 =   /usr/lib/x86_64-linux-gnu/
    88FFLAGS   =   -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
     9#FFLAGS   =   -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
     10LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
    911
    10 ifeq ($(TARGET),dmz)   
    11   # options for ganglia
    12   INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
    13   LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
    14   LIBPATH2 = /usr/lib/x86_64-linux-gnu/
    15   MAIN = FLEXPART_dmz_
    16 endif
    17 ifeq ($(TARGET),local)
    18   # local options       
    19   #libs_dir=/.../flexpart/libs/
    20   libs_dir=/Users/ignacio/flexpart/libs/
    21   INCPATH  = $(libs_dir)/grib_api-1.9.9_dir/include
    22   LIBPATH1 = $(libs_dir)/grib_api-1.9.9_dir/lib
    23   LIBPATH2 = $(libs_dir)/jasper_dir/lib
    24   MAIN = FLEXPART_local_
    25 endif
    26 
    27 LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
    2812
    2913MODOBJS = \
    … …  
    3721OBJECTS = \
    3822writeheader.o  writeheader_txt.o   writeheader_surf.o       assignland.o\
    39                part0.o \
     23calcpar.o               part0.o \
    4024caldate.o               partdep.o \
    4125coordtrafo.o            psih.o \
    … …  
    5337getrc.o                 readreleases.o \
    5438getvdep.o               readspecies.o \
    55 interpol_misslev.o      \
    56 conccalc.o              \
     39interpol_misslev.o      readwind.o \
     40conccalc.o              richardson.o \
    5741concoutput.o  concoutput_surf.o          scalev.o \
    5842pbl_profile.o           readOHfield.o\
    5943juldate.o               timemanager.o \
    6044interpol_vdep.o         interpol_rain.o \
    61 partoutput.o \
     45verttransform.o         partoutput.o \
    6246hanna.o                 wetdepokernel.o \
    6347mean.o                  wetdepo.o \
    6448hanna_short.o           windalign.o \
     49obukhov.o               gridcheck.o \
    6550hanna1.o                initialize.o \
    66                            calcpar_nests.o \
     51                        gridcheck_nests.o \
     52readwind_nests.o        calcpar_nests.o \
    6753verttransform_nests.o   interpol_all_nests.o \
    6854interpol_wind_nests.o   interpol_misslev_nests.o \
    6955interpol_vdep_nests.o   interpol_rain_nests.o \
    70 getvdep_nests.o   gridcheck_nests.o \
    71 readwind_nests.o  \
     56getvdep_nests.o \
    7257readageclasses.o        readpartpositions.o \
    7358calcfluxes.o            fluxoutput.o \
    7459qvsat.o                 skplin.o \
     60convmix.o               calcmatrix.o \
    7561convect43c.o               redist.o \
    7662sort2.o                 distance.o \
    … …  
    9177
    9278
    93 ifeq ($(WINDS),ecmwf)
    94   OBJECTS_WINDS = \
    95   calcpar.o          readwind.o \
    96   richardson.o       verttransform.o \
    97   obukhov.o          gridcheck.o  \
    98   convmix.o          calcmatrix.o
    99 endif
    100 
    101 ifeq ($(WINDS),gfs)
    102   OBJECTS_WINDS = \
    103   calcpar_gfs.o          readwind_gfs.o \
    104   richardson_gfs.o       verttransform_gfs.o \
    105   obukhov_gfs.o          gridcheck_gfs.o  \
    106   convmix_gfs.o          calcmatrix_gfs.o
    107 endif
    108 
    109 $(MAIN): $(MODOBJS) $(OBJECTS)  $(OBJECTS_WINDS)
    110         $(FC) *.o -o $(MAIN)_$(WINDS) $(LDFLAGS)
     79$(MAIN): $(MODOBJS) $(OBJECTS)
     80        $(FC) *.o -o $(MAIN) $(LDFLAGS)
    11181
    11282$(OBJECTS): $(MODOBJS)
    … …  
    11888        rm *.o *.mod
    11989
     90cleanall:
     91        rm *.o *.mod $(MAIN)

  • trunk/src/par_mod.f90

    r20 r24  
    121121  !*********************************************
    122122 
    123   integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !FNL
     123  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !FNL XF
     124  integer,parameter :: nxmax=361,nymax=181,nuvzmax=152,nwzmax=152,nzmax=152 !ECMWF new
    124125  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF
    125126  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26
    … …  
    154155
    155156  !integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
    156   !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151 !ECMWF
    157   integer,parameter :: maxnests=1, nxmaxn=201, nymaxn=161 ! FNL XF
     157  integer,parameter :: maxnests=1,nxmaxn=351,nymaxn=351 !ECMWF
     158  !integer,parameter :: maxnests=1, nxmaxn=201, nymaxn=161 ! FNL XF
    158159  ! maxnests                maximum number of nested grids
    159160  ! nxmaxn,nymaxn           maximum dimension of nested wind fields in
    … …  
    197198  !**************************************************
    198199
    199   !integer,parameter :: maxpart=4000000
    200   integer,parameter :: maxpart=2000
    201   integer,parameter :: maxspec=6
     200  integer,parameter :: maxpart=15000000
     201  integer,parameter :: maxspec=4
    202202
    203203

  • trunk/src/readcommand.f90

    r20 r24  
    8080  character(len=50) :: line
    8181  logical :: old
    82   logical :: nmlout=.true. !.false.
     82  logical :: nml_COMMAND=.true. , nmlout=.true. !.false.
    8383  integer :: readerror
    8484
    8585  namelist /command/ &
    86     ldirect, &
    87     ibdate,ibtime, &
    88     iedate,ietime, &
    89     loutstep, &
    90     loutaver, &
    91     loutsample, &
    92     itsplit, &
    93     lsynctime, &
    94     ctl, &
    95     ifine, &
    96     iout, &
    97     ipout, &
    98     lsubgrid, &
    99     lconvection, &
    100     lagespectra, &
    101     ipin, &
    102     ioutputforeachrelease, &
    103     iflux, &
    104     mdomainfill, &
    105     ind_source, &
    106     ind_receptor, &
    107     mquasilag, &
    108     nested_output, &
    109     linit_cond, &
    110     surf_only   
     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   
    111111
    112112  ! Presetting namelist command
    … …  
    144144  open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
    145145         form='formatted',iostat=readerror)
    146   ! If fail, check if file does not exist
    147   if (readerror.ne.0) then
    148    
    149     print*,'***ERROR: file COMMAND not found in '
    150     print*, path(1)(1:length(1))//'COMMAND'
    151     print*, 'Check your pathnames file.'
    152     stop
    153 
    154   endif
    155 
    156   read(unitcommand,command,iostat=readerror)
     146   ! If fail, check if file does not exist
     147   if (readerror.ne.0) then
     148     print*,'***ERROR: file COMMAND not found in '
     149     print*, path(1)(1:length(1))//'COMMAND'
     150     print*, 'Check your pathnames file.'
     151     stop
     152   endif
     153
     154   ! print error code
     155   !write(*,*) 'readcommand > readerror open=' , readerror
     156   !probe first line 
     157   read (unitcommand,901) line
     158   !write(*,*) 'index(line,COMMAND) =', index(line,'COMMAND')
     159
     160 
     161   !default is namelist input
     162   ! distinguish namelist from fixed text input
     163   if (index(line,'COMMAND') .eq. 0) then
     164   nml_COMMAND = .false.
     165   !write(*,*) 'COMMAND file does not contain the string COMMAND in the first line'     
     166   endif
     167   !write(*,*) 'readcommand > read as namelist? ' , nml_COMMAND
     168   rewind(unitcommand)
     169   read(unitcommand,command,iostat=readerror)
     170
    157171  close(unitcommand)
    158172
     173  !write(*,*) 'readcommand > readerror read=' , readerror
    159174  ! If error in namelist format, try to open with old input code
    160   if (readerror.ne.0) then
     175  ! if (readerror.ne.0) then
     176  ! IP 21/5/2014 the previous line cause the old long format
     177  ! to be confused with namelist input
     178 
     179  ! use text input
     180  if (nml_COMMAND .eqv. .false.) then
    161181
    162182    open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
    … …  
    173193    if (index(line,'LDIRECT') .eq. 0) then
    174194      old = .false.
     195    !write(*,*) 'readcommand old short'
    175196    else
    176197      old = .true.
     198    !write(*,*) 'readcommand old long'
    177199    endif
    178200    rewind(unitcommand)
    179201
     202
    180203    ! Read parameters
    181204    !****************
    … …  
    231254    if (old) call skplin(3,unitcommand)
    232255    read(unitcommand,*) linit_cond
     256    if (old) call skplin(3,unitcommand)
     257    read(unitcommand,*) surf_only
    233258    close(unitcommand)
    234259
    … …  
    237262  ! write command file in namelist format to output directory if requested
    238263  if (nmlout.eqv..true.) then
    239     !open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist.out',status='new',err=1000)
    240264    open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist',err=1000)
    241265    write(unitcommand,nml=command)

  • trunk/src/readreleases.f90

    r20 r24  
    5555  ! xpoint2,ypoint2     geograf. coordinates of upper right corner of release  *
    5656  !                     area                                                   *
    57   ! weta, wetb          parameters to determine the wet scavenging coefficient *
     57  ! weta, wetb          parameters to determine the below-cloud scavenging     *
     58  ! weta_in, wetb_in    parameters to determine the in-cloud scavenging        *
     59  ! wetc_in, wetd_in    parameters to determine the in-cloud scavenging        *
    5860  ! zpoint1,zpoint2     height range, over which release takes place           *
    5961  ! num_min_discrete    if less, release cannot be randomized and happens at   *
    … …  
    272274        vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
    273275      end do
    274       write(*,*) 'Average setting velocity: ',i,vsetaver(i)
     276      write(*,*) 'Average settling velocity: ',i,vsetaver(i)
    275277    endif
    276278
    … …  
    284286    if (weta(i).gt.0.)  then
    285287      WETDEP=.true.
    286       write (*,*) 'Wetdeposition switched on: ',weta(i),i
    287     endif
     288      write (*,*) 'Below-cloud scavenging is switched on'
     289      write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
     290    else
     291      write (*,*) 'Below-cloud scavenging is switched OFF'
     292    endif
     293   
     294! NIK 31.01.2013 + 10.12.2013
     295    if (weta_in(i).gt.0.)  then
     296      WETDEP=.true.
     297      write (*,*) 'In-cloud scavenging is switched on'
     298      write (*,*) 'In-cloud scavenging coefficients: ',weta_in(i),wetb_in(i), wetc_in(i), wetd_in(i),i
     299    else
     300      write (*,*) 'In-cloud scavenging is switched OFF'
     301    endif
     302
    288303    if (ohreact(i).gt.0) then
    289304      OHREA=.true.

  • trunk/src/readspecies.f90

    r4 r24  
    3131  !   11 July 1996                                                             *
    3232  !                                                                            *
     33  !   Changes:                                                                 *
     34  !   N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging *
     35  !                                                                            *
    3336  !*****************************************************************************
    3437  !                                                                            *
    … …  
    3639  ! decaytime(maxtable)  half time for radiological decay                      *
    3740  ! specname(maxtable)   names of chemical species, radionuclides              *
    38   ! wetscava, wetscavb   Parameters for determining scavenging coefficient     *
     41  ! weta, wetb           Parameters for determining below-cloud scavenging     *
     42  ! weta_in              Parameter for determining in-cloud scavenging         *
     43  ! wetb_in              Parameter for determining in-cloud scavenging         *
     44  ! wetc_in              Parameter for determining in-cloud scavenging         *
     45  ! wetd_in              Parameter for determining in-cloud scavenging         *
    3946  ! ohreact              OH reaction rate                                      *
    4047  ! id_spec              SPECIES number as referenced in RELEASE file          *
    … …  
    7885    read(unitspecies,'(f18.2)',end=22) wetb(pos_spec)
    7986  !  write(*,*) wetb(pos_spec)
     87
     88!*** NIK 31.01.2013: including in-cloud scavening parameters
     89   read(unitspecies,'(e18.1)',end=22) weta_in(pos_spec)
     90  !  write(*,*) weta_in(pos_spec)
     91   read(unitspecies,'(f18.2)',end=22) wetb_in(pos_spec)
     92  !  write(*,*) wetb_in(pos_spec)
     93   read(unitspecies,'(f18.2)',end=22) wetc_in(pos_spec)
     94  !  write(*,*) wetc_in(pos_spec)
     95   read(unitspecies,'(f18.2)',end=22) wetd_in(pos_spec)
     96  !  write(*,*) wetd_in(pos_spec)
     97
    8098    read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
    8199  !  write(*,*) reldiff(pos_spec)

  • trunk/src/readwind.f90

    r20 r24  
    7171  integer :: iret
    7272  integer :: igrib
    73   integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
     73  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
    7474  integer :: gotGrid
    7575  !HSO  end
    … …  
    9090  integer :: isec1(56),isec2(22+nxmax+nymax)
    9191  real(kind=4) :: zsec4(jpunp)
    92   real(kind=4) :: xaux,yaux,xaux0,yaux0
     92  real(kind=4) :: xaux,yaux
    9393  real(kind=8) :: xauxin,yauxin
    9494  real,parameter :: eps=1.e-4
    9595  real(kind=4) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
    96   real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
     96  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
    9797
    9898  logical :: hflswitch,strswitch
    … …  
    101101  character(len=24) :: gribErrorMsg = 'Error reading grib file'
    102102  character(len=20) :: gribFunction = 'readwind'
    103 
    104   !HSO conversion of ECMWF etadot to etadot*dp/deta
    105   logical :: etacon=.false.
    106   real,parameter :: p00=101325.
    107   real :: dak,dbk
    108103
    109104  hflswitch=.false.
    … …  
    155150  endif
    156151
     152  conversion_factor=1.
     153
    157154  else
    158155
    … …  
    168165  call grib_check(iret,gribFunction,gribErrorMsg)
    169166  call grib_get_int(igrib,'level',valSurf,iret)
     167  call grib_check(iret,gribFunction,gribErrorMsg)
     168  call grib_get_int(igrib,'paramId',parId,iret)
    170169  call grib_check(iret,gribFunction,gribErrorMsg)
    171170
    … …  
    177176  isec1(8)=-1
    178177  isec1(8)=valSurf     ! level
     178  conversion_factor=1.
    179179  if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
    180180    isec1(6)=130         ! indicatorOfParameter
    … …  
    188188    isec1(6)=134         ! indicatorOfParameter
    189189  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
     190    isec1(6)=135         ! indicatorOfParameter
     191  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
    190192    isec1(6)=135         ! indicatorOfParameter
    191193  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    … …  
    201203  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    202204    isec1(6)=141         ! indicatorOfParameter
    203   elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
     205    conversion_factor=1000.
     206  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
    204207    isec1(6)=164         ! indicatorOfParameter
    205   elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
     208  elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
    206209    isec1(6)=142         ! indicatorOfParameter
    207210  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    208211    isec1(6)=143         ! indicatorOfParameter
     212    conversion_factor=1000.
    209213  elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
    210214    isec1(6)=146         ! indicatorOfParameter
    211215  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    212216    isec1(6)=176         ! indicatorOfParameter
    213   elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
     217  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
    214218    isec1(6)=180         ! indicatorOfParameter
    215   elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
     219  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
    216220    isec1(6)=181         ! indicatorOfParameter
    217221  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    218222    isec1(6)=129         ! indicatorOfParameter
    219   elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
     223  elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
    220224    isec1(6)=160         ! indicatorOfParameter
    221225  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    … …  
    223227    isec1(6)=172         ! indicatorOfParameter
    224228  else
    225     print*,'***ERROR: undefined GRiB2 message found!',discipl, &
     229    print*,'***WARNING: undefined GRiB2 message found!',discipl, &
    226230         parCat,parNum,typSurf
     231  endif
     232  if(parId .ne. isec1(6) .and. parId .ne. 77) then
     233    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
     234!    stop
    227235  endif
    228236
    … …  
    237245  !HSO  get the required fields from section 2 in a gribex compatible manner
    238246  if (ifield.eq.1) then
    239   call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
    240        isec2(2),iret)
    241   call grib_check(iret,gribFunction,gribErrorMsg)
    242   call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
    243        isec2(3),iret)
    244   call grib_check(iret,gribFunction,gribErrorMsg)
    245   call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
    246        isec2(12))
     247  call grib_get_int(igrib,'numberOfPointsAlongAParallel',isec2(2),iret)
     248  call grib_check(iret,gribFunction,gribErrorMsg)
     249  call grib_get_int(igrib,'numberOfPointsAlongAMeridian',isec2(3),iret)
     250  call grib_check(iret,gribFunction,gribErrorMsg)
     251  call grib_get_int(igrib,'numberOfVerticalCoordinateValues',isec2(12))
    247252  call grib_check(iret,gribFunction,gribErrorMsg)
    248253  ! CHECK GRID SPECIFICATIONS
    … …  
    250255  if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
    251256  if(isec2(12)/2-1.ne.nlev_ec) &
    252        stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
     257  stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
    253258  endif ! ifield
    254259
    255260  !HSO  get the second part of the grid dimensions only from GRiB1 messages
    256   if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
     261  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
    257262    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
    258263         xauxin,iret)
    … …  
    261266         yauxin,iret)
    262267    call grib_check(iret,gribFunction,gribErrorMsg)
     268    if (xauxin.gt.180.) xauxin=xauxin-360.0
     269    if (xauxin.lt.-180.) xauxin=xauxin+360.0
     270
    263271    xaux=xauxin+real(nxshift)*dx
    264272    yaux=yauxin
    265     xaux0=xlon0
    266     yaux0=ylat0
    267     if(xaux.lt.0.) xaux=xaux+360.
    268     if(yaux.lt.0.) yaux=yaux+360.
    269     if(xaux0.lt.0.) xaux0=xaux0+360.
    270     if(yaux0.lt.0.) yaux0=yaux0+360.
    271     if(abs(xaux-xaux0).gt.eps) &
    272          stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
    273     if(abs(yaux-yaux0).gt.eps) &
    274          stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
     273    if (xaux.gt.180.) xaux=xaux-360.0
     274    if(abs(xaux-xlon0).gt.eps) &
     275    stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
     276    if(abs(yaux-ylat0).gt.eps) &
     277    stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
    275278    gotGrid=1
    276279  endif ! gotGrid
    … …  
    298301           zsec4(nxfield*(ny-j-1)+i+1)
    299302      if(isec1(6).eq.141) sd(i,j,1,n)= &!! SNOW DEPTH
    300            zsec4(nxfield*(ny-j-1)+i+1)
     303           zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
    301304      if(isec1(6).eq.151) msl(i,j,1,n)= &!! SEA LEVEL PRESS.
    302305           zsec4(nxfield*(ny-j-1)+i+1)
    … …  
    316319      endif
    317320      if(isec1(6).eq.143) then                      !! CONVECTIVE PREC.
    318         convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
     321        convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
    319322        if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
    320323      endif
    … …  
    366369      do j=0,nymin1
    367370        wwh(i,j,nlev_ec+1)=0.
    368       end do
    369     end do
    370   endif
    371 
    372   ! convert from ECMWF etadot to etadot*dp/deta as needed by FLEXPART
    373   if(etacon.eqv..true.) then
    374     do k=1,nwzmax
    375       dak=akm(k+1)-akm(k)
    376       dbk=bkm(k+1)-bkm(k)
    377       do i=0,nxmin1
    378         do j=0,nymin1
    379           wwh(i,j,k)=2*wwh(i,j,k)*ps(i,j,1,n)*(dak/ps(i,j,1,n)+dbk)/(dak/p00+dbk)
    380           if (k.gt.1) then
    381             wwh(i,j,k)=wwh(i,j,k)-wwh(i,j,k-1)
    382           endif
    383         end do
    384371      end do
    385372    end do
    … …  
    479466
    480467end subroutine readwind
     468

  • trunk/src/readwind_nests.f90

    r4 r24  
    5151  integer :: igrib
    5252  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
     53  integer :: parId !!added by mc for making it consistent with new readwind.f90
    5354  integer :: gotGrid
    5455  !HSO  end
    … …  
    6970  integer :: isec1(56),isec2(22+nxmaxn+nymaxn)
    7071  real(kind=4) :: zsec4(jpunp)
     72  real(kind=4) :: xaux,yaux
    7173  real(kind=8) :: xauxin,yauxin
    72   real(kind=4) :: xaux,yaux,xaux0,yaux0
     74  real,parameter :: eps=1.e-4
    7375  real :: ewss(0:nxmaxn-1,0:nymaxn-1),nsss(0:nxmaxn-1,0:nymaxn-1)
    7476  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
     77  real :: conversion_factor !added by mc to make it consistent with new gridchek.f90
    7578
    7679  logical :: hflswitch,strswitch
    … …  
    134137  endif
    135138
     139  conversion_factor=1.
     140
     141
    136142  else
    137143
    … …  
    148154  call grib_get_int(igrib,'level',valSurf,iret)
    149155  call grib_check(iret,gribFunction,gribErrorMsg)
     156  call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consisitent with new readwind.f90
     157  call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consisitent with new readwind.f90
    150158
    151159  !print*,discipl,parCat,parNum,typSurf,valSurf
    … …  
    156164  isec1(8)=-1
    157165  isec1(8)=valSurf     ! level
     166   conversion_factor=1.
    158167  if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
    159168    isec1(6)=130         ! indicatorOfParameter
    … …  
    166175  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
    167176    isec1(6)=134         ! indicatorOfParameter
    168   elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
     177  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot !
    169178    isec1(6)=135         ! indicatorOfParameter
     179  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added by mc to make it consisitent with new readwind.f90
     180    isec1(6)=135         ! indicatorOfParameter    !added by mc to make it consisitent with new readwind.f90
    170181  elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
    171182    isec1(6)=151         ! indicatorOfParameter
    … …  
    180191  elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
    181192    isec1(6)=141         ! indicatorOfParameter
    182   elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
     193    conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
     194  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consisitent with new readwind.f90
    183195    isec1(6)=164         ! indicatorOfParameter
    184   elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
     196  elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consisitent with new readwind.f90
    185197    isec1(6)=142         ! indicatorOfParameter
    186198  elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
    187199    isec1(6)=143         ! indicatorOfParameter
     200    conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
    188201  elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
    189202    isec1(6)=146         ! indicatorOfParameter
    190203  elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    191204    isec1(6)=176         ! indicatorOfParameter
    192   elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
     205  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consisitent with new readwind.f90
    193206    isec1(6)=180         ! indicatorOfParameter
    194   elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
     207  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consisitent with new readwind.f90
    195208    isec1(6)=181         ! indicatorOfParameter
    196209  elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
    197210    isec1(6)=129         ! indicatorOfParameter
    198   elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
     211   elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consisitent with new readwind.f90
    199212    isec1(6)=160         ! indicatorOfParameter
    200213  elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
    … …  
    202215    isec1(6)=172         ! indicatorOfParameter
    203216  else
    204     print*,'***ERROR: undefined GRiB2 message found!',discipl, &
     217    print*,'***WARNING: undefined GRiB2 message found!',discipl, &
    205218         parCat,parNum,typSurf
     219  endif
     220  if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consisitent with new readwind.f90
     221    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)  !
     222!    stop
    206223  endif
    207224
    … …  
    227244  ! CHECK GRID SPECIFICATIONS
    228245  if(isec2(2).ne.nxn(l)) stop &
    229        'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
     246  'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
    230247  if(isec2(3).ne.nyn(l)) stop &
    231        'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
     248  'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
    232249  if(isec2(12)/2-1.ne.nlev_ec) stop 'READWIND: VERTICAL DISCRET&
    233        &IZATION NOT CONSISTENT FOR A NESTING LEVEL'
     250  IZATION NOT CONSISTENT FOR A NESTING LEVEL'
    234251  endif ! ifield
    235252
    236253  !HSO  get the second part of the grid dimensions only from GRiB1 messages
    237   if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
     254 if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consisitent with new readwind.f90
    238255    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
    239256         xauxin,iret)
    … …  
    242259         yauxin,iret)
    243260    call grib_check(iret,gribFunction,gribErrorMsg)
     261    if (xauxin.gt.180.) xauxin=xauxin-360.0
     262    if (xauxin.lt.-180.) xauxin=xauxin+360.0
     263
    244264    xaux=xauxin
    245265    yaux=yauxin
    246     xaux0=xlon0n(l)
    247     yaux0=ylat0n(l)
    248     if(xaux.lt.0.) xaux=xaux+360.
    249     if(yaux.lt.0.) yaux=yaux+360.
    250     if(xaux0.lt.0.) xaux0=xaux0+360.
    251     if(yaux0.lt.0.) yaux0=yaux0+360.
    252     if(xaux.ne.xaux0) &
    253          stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT FOR A NES&
    254          &TING LEVEL'
    255     if(yaux.ne.yaux0) &
    256          stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT FOR A NEST&
    257          &ING LEVEL'
     266    if (abs(xaux-xlon0n(l)).gt.eps) &
     267    stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT FOR A NESTING LEVEL'
     268    if (abs(yaux-ylat0n(l)).gt.eps) &
     269    stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT FOR A NESTING LEVEL'
    258270    gotGrid=1
    259271  endif
    … …  
    280292             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
    281293        if(isec1(6).eq.141) sdn(i,j,1,n,l)= &!! SNOW DEPTH
    282              zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
     294             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90!
    283295        if(isec1(6).eq.151) msln(i,j,1,n,l)= &!! SEA LEVEL PRESS.
    284296             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
    … …  
    298310        endif
    299311        if(isec1(6).eq.143) then                         !! CONVECTIVE PREC.
    300           convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
     312          convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90
    301313          if (convprecn(i,j,1,n,l).lt.0.) convprecn(i,j,1,n,l)=0.
    302314        endif

  • trunk/src/verttransform.f90

    r20 r24  
    2121
    2222subroutine verttransform(n,uuh,vvh,wwh,pvh)
    23   !                         i  i   i   i   i
     23  !                      i  i   i   i   i
    2424  !*****************************************************************************
    2525  !                                                                            *
    … …  
    4949  ! Sabine Eckhardt, March 2007
    5050  ! added the variable cloud for use with scavenging - descr. in com_mod
    51   ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
    52   ! note that also other subroutines are affected by the fix
    5351  !*****************************************************************************
    5452  !                                                                            *
    … …  
    7270
    7371  integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
    74  integer :: rain_cloud_above,kz_inv !SE
    75   integer icloudtop !PS
     72  integer :: rain_cloud_above,kz_inv
    7673  real :: f_qvsat,pressure
    77  !real :: rh,lsp,convp
    78   real :: rh,lsp,convp,prec,rhmin 
    79   real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
     74  real :: rh,lsp,convp
     75  real :: rhoh(nuvzmax),pinmconv(nzmax)
    8076  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
    8177  real :: xlon,ylat,xlonr,dzdx,dzdy
    82   real :: dzdx1,dzdx2,dzdy1,dzdy2, precmin
     78  real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
    8379  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
    8480  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
    … …  
    8783  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
    8884  real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
    89   logical lconvectprec
    9085  real,parameter :: const=r_air/ga
    91   parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
    9286
    9387  logical :: init = .true.
    … …  
    9791  ! If verttransform is called the first time, initialize heights of the   *
    9892  ! z levels in meter. The heights are the heights of model levels, where  *
    99   ! u,v,T and qv are given, and of the interfaces, where w is given. So,   *
    100   ! the vertical resolution in the z system is doubled. As reference point,*
    101   ! the lower left corner of the grid is used.                             *
    102   ! Unlike in the eta system, no difference between heights for u,v and    *
    103   ! heights for w exists.                                                  *
     93  ! u,v,T and qv are given.                                                *
    10494  !*************************************************************************
    105 
    106 
    107   !      do 897 kz=1,nuvz
    108   !       write (*,*) 'akz: ',akz(kz),'bkz',bkz(kz)
    109   !897     continue
    11095
    11196  if (init) then
    … …  
    11398  ! Search for a point with high surface pressure (i.e. not above significant topography)
    11499  ! Then, use this point to construct a reference z profile, to be used at all times
    115   !*****************************************************************************
     100  !**************************************************************************************
    116101
    117102    do jy=0,nymin1
    … …  
    128113
    129114    tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
    130          ps(ixm,jym,1,n))
     115    ps(ixm,jym,1,n))
    131116    pold=ps(ixm,jym,1,n)
    132117    height(1)=0.
    … …  
    136121      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
    137122
    138 
    139   ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
    140   ! upon the transformation to z levels. In order to save computer memory, this is
    141   ! not done anymore in the standard version. However, this option can still be
    142   ! switched on by replacing the following lines with those below, that are
    143   ! currently commented out.
    144   ! Note that two more changes are necessary in this subroutine below.
    145   ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
    146   !*****************************************************************************
    147 
    148123      if (abs(tv-tvold).gt.0.2) then
    149         height(kz)= &
    150              height(kz-1)+const*log(pold/pint)* &
    151              (tv-tvold)/log(tv/tvold)
     124        height(kz)=height(kz-1)+const*log(pold/pint)* &
     125        (tv-tvold)/log(tv/tvold)
    152126      else
    153         height(kz)=height(kz-1)+ &
    154              const*log(pold/pint)*tv
     127        height(kz)=height(kz-1)+const*log(pold/pint)*tv
    155128      endif
    156 
    157   ! Switch on following lines to use doubled vertical resolution
    158   !*************************************************************
    159   !    if (abs(tv-tvold).gt.0.2) then
    160   !      height((kz-1)*2)=
    161   !    +      height(max((kz-2)*2,1))+const*log(pold/pint)*
    162   !    +      (tv-tvold)/log(tv/tvold)
    163   !    else
    164   !      height((kz-1)*2)=height(max((kz-2)*2,1))+
    165   !    +      const*log(pold/pint)*tv
    166   !    endif
    167   ! End doubled vertical resolution
    168129
    169130      tvold=tv
    170131      pold=pint
    171132    end do
    172 
    173 
    174   ! Switch on following lines to use doubled vertical resolution
    175   !*************************************************************
    176   !  do 7 kz=3,nz-1,2
    177   !    height(kz)=0.5*(height(kz-1)+height(kz+1))
    178   !  height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
    179   ! End doubled vertical resolution
    180133
    181134
    … …  
    204157  do jy=0,nymin1
    205158    do ix=0,nxmin1
    206       tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
    207            ps(ix,jy,1,n))
     159      tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ps(ix,jy,1,n))
    208160      pold=ps(ix,jy,1,n)
    209       uvzlev(1)=0.
     161      uvwzlev(ix,jy,1)=0.
    210162      wzlev(1)=0.
    211163      rhoh(1)=pold/(r_air*tvold)
    … …  
    221173
    222174        if (abs(tv-tvold).gt.0.2) then
    223           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
    224                (tv-tvold)/log(tv/tvold)
     175          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
     176          (tv-tvold)/log(tv/tvold)
    225177        else
    226           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
     178          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
    227179        endif
    228180
    … …  
    233185
    234186      do kz=2,nwz-1
    235         wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
    236       end do
    237       wzlev(nwz)=wzlev(nwz-1)+ &
    238            uvzlev(nuvz)-uvzlev(nuvz-1)
    239 
    240       uvwzlev(ix,jy,1)=0.0
    241       do kz=2,nuvz
    242         uvwzlev(ix,jy,kz)=uvzlev(kz)
    243       end do
    244 
    245   ! Switch on following lines to use doubled vertical resolution
    246   ! Switch off the three lines above.
    247   !*************************************************************
    248   !22          uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
    249   !     do 23 kz=2,nwz
    250   !23          uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
    251   ! End doubled vertical resolution
     187        wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2.
     188      end do
     189      wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1)
    252190
    253191  ! pinmconv=(h2-h1)/(p2-p1)
    254192
    255193      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
    256            ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
    257            (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
     194      ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
     195      (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
    258196      do kz=2,nz-1
    259197        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
    260              ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
    261              (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
     198        ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
     199        (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
    262200      end do
    263201      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
    264            ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
    265            (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
     202      ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
     203      (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
    266204
    267205  ! Levels, where u,v,t and q are given
    … …  
    283221      do iz=2,nz-1
    284222        do kz=kmin,nuvz
    285           if(height(iz).gt.uvzlev(nuvz)) then
     223          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
    286224            uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
    287225            vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
    … …  
    292230            goto 30
    293231          endif
    294           if ((height(iz).gt.uvzlev(kz-1)).and. &
    295                (height(iz).le.uvzlev(kz))) then
    296            dz1=height(iz)-uvzlev(kz-1)
    297            dz2=uvzlev(kz)-height(iz)
     232          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
     233          (height(iz).le.uvwzlev(ix,jy,kz))) then
     234           dz1=height(iz)-uvwzlev(ix,jy,kz-1)
     235           dz2=uvwzlev(ix,jy,kz)-height(iz)
    298236           dz=dz1+dz2
    299237           uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
    … …  
    339277
    340278      drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
    341            (height(2)-height(1))
     279      (height(2)-height(1))
    342280      do kz=2,nz-1
    343281        drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
    344              (height(kz+1)-height(kz-1))
     282        (height(kz+1)-height(kz-1))
    345283      end do
    346284      drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
    … …  
    356294
    357295  do jy=1,ny-2
     296    cosf=cos((real(jy)*dy+ylat0)*pi180)
    358297    do ix=1,nx-2
    359298
    … …  
    361300      do iz=2,nz-1
    362301
    363         ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
     302        ui=uu(ix,jy,iz,n)*dxconst/cosf
    364303        vi=vv(ix,jy,iz,n)*dyconst
    365304
    366305        do kz=kmin,nz
    367306          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
    368                (height(iz).le.uvwzlev(ix,jy,kz))) then
     307          (height(iz).le.uvwzlev(ix,jy,kz))) then
    369308            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
    370309            dz2=uvwzlev(ix,jy,kz)-height(iz)
    … …  
    409348        do iz=1,nz
    410349          call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
    411                vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
    412                vvpol(ix,jy,iz,n))
     350          vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
    413351        end do
    414352      end do
    … …  
    424362      xlon=xlon0+real(nx/2-1)*dx
    425363      xlonr=xlon*pi/180.
    426       ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
    427            vv(nx/2-1,nymin1,iz,n)**2)
     364      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2)
    428365      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
    429         ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
    430              vv(nx/2-1,nymin1,iz,n))-xlonr
     366        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
    431367      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
    432368        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
    433              vv(nx/2-1,nymin1,iz,n))-xlonr
     369        vv(nx/2-1,nymin1,iz,n))-xlonr
    434370      else
    435371        ddpol=pi/2-xlonr
    … …  
    444380      uuaux=-ffpol*sin(xlonr+ddpol)
    445381      vvaux=-ffpol*cos(xlonr+ddpol)
    446       call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
    447            vvpolaux)
    448 
     382      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
    449383      jy=nymin1
    450384      do ix=0,nxmin1
    … …  
    485419        do iz=1,nz
    486420          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
    487                vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
    488                vvpol(ix,jy,iz,n))
     421          vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
    489422        end do
    490423      end do
    … …  
    499432      xlon=xlon0+real(nx/2-1)*dx
    500433      xlonr=xlon*pi/180.
    501       ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
    502            vv(nx/2-1,0,iz,n)**2)
     434      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2)
    503435      if (vv(nx/2-1,0,iz,n).lt.0.) then
    504         ddpol=atan(uu(nx/2-1,0,iz,n)/ &
    505              vv(nx/2-1,0,iz,n))+xlonr
     436        ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
    506437      else if (vv(nx/2-1,0,iz,n).gt.0.) then
    507         ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
    508              vv(nx/2-1,0,iz,n))+xlonr
     438        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
    509439      else
    510440        ddpol=pi/2-xlonr
    … …  
    519449      uuaux=+ffpol*sin(xlonr-ddpol)
    520450      vvaux=-ffpol*cos(xlonr-ddpol)
    521       call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
    522            vvpolaux)
     451      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
    523452
    524453      jy=0
    … …  
    551480  !   create a cloud and rainout/washout field, clouds occur where rh>80%
    552481  !   total cloudheight is stored at level 0
    553 
    554 
    555 
    556482  do jy=0,nymin1
    557483    do ix=0,nxmin1
    558 
    559 
    560 
    561   !    rain_cloud_above=0
    562   !    lsp=lsprec(ix,jy,1,n)
    563   !    convp=convprec(ix,jy,1,n)
    564   !    cloudsh(ix,jy,n)=0
    565   !    do kz_inv=1,nz-1
    566   !       kz=nz-kz_inv+1
    567   !       pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
    568   !       rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
    569   !       clouds(ix,jy,kz,n)=0
    570   !       if (rh.gt.0.8) then ! in cloud
    571   !          if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
    572   !             rain_cloud_above=1
    573   !             cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
    574   !                  height(kz)-height(kz-1)
    575   !             if (lsp.ge.convp) then
    576   !                clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
    577   !             else
    578   !                clouds(ix,jy,kz,n)=2 ! convp dominated rainout
    579   !             endif
    580   !          else ! no precipitation
    581   !                clouds(ix,jy,kz,n)=1 ! cloud
    582   !          endif
    583   !       else ! no cloud
    584   !          if (rain_cloud_above.eq.1) then ! scavenging
    585   !             if (lsp.ge.convp) then
    586   !                clouds(ix,jy,kz,n)=5 ! lsp dominated washout
    587   !             else
    588   !                clouds(ix,jy,kz,n)=4 ! convp dominated washout
    589   !             endif
    590   !          endif
    591   !       endif
    592   !    end do
    593 
    594 
    595    ! PS 3012
    596 
    597              lsp=lsprec(ix,jy,1,n)
    598           convp=convprec(ix,jy,1,n)
    599           prec=lsp+convp
    600           if (lsp.gt.convp) then !  prectype='lsp'
    601             lconvectprec = .false.
    602           else ! prectype='cp '
    603             lconvectprec = .true.
    604           endif
    605           rhmin = 0.90 ! standard condition for presence of clouds
    606 !PS       note that original by Sabine Eckhart was 80%
    607 !PS       however, for T<-20 C we consider saturation over ice
    608 !PS       so I think 90% should be enough         
    609           icloudbot(ix,jy,n)=icmv
    610           icloudtop=icmv ! this is just a local variable
    611 98        do kz=1,nz
    612             pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
    613             rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
    614 !ps            if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
    615             if (rh .gt. rhmin) then
    616               if (icloudbot(ix,jy,n) .eq. icmv) then
    617                 icloudbot(ix,jy,n)=nint(height(kz))
    618               endif
    619               icloudtop=nint(height(kz)) ! use int to save memory
     484      rain_cloud_above=0
     485      lsp=lsprec(ix,jy,1,n)
     486      convp=convprec(ix,jy,1,n)
     487      cloudsh(ix,jy,n)=0
     488      do kz_inv=1,nz-1
     489         kz=nz-kz_inv+1
     490         pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
     491         rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
     492         clouds(ix,jy,kz,n)=0
     493         if (rh.gt.0.8) then ! in cloud
     494            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
     495               rain_cloud_above=1
     496               cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
     497               height(kz)-height(kz-1)
     498               if (lsp.ge.convp) then
     499                  clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
     500               else
     501                  clouds(ix,jy,kz,n)=2 ! convp dominated rainout
     502               endif
     503            else ! no precipitation
     504                  clouds(ix,jy,kz,n)=1 ! cloud
    620505            endif
    621           enddo
    622 
    623 !PS try to get a cloud thicker than 50 m
    624 !PS if there is at least .01 mm/h  - changed to 0.002 and put into
    625 !PS parameter precpmin       
    626           if ((icloudbot(ix,jy,n) .eq. icmv .or. &
    627               icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
    628               prec .gt. precmin) then
    629             rhmin = rhmin - 0.05
    630             if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
    631           endif
    632 !PS implement a rough fix for badly represented convection
    633 !PS is based on looking at a limited set of comparison data
    634           if (lconvectprec .and. icloudtop .lt. 6000 .and. &
    635              prec .gt. precmin) then
    636 
    637             if (convp .lt. 0.1) then
    638               icloudbot(ix,jy,n) = 500
    639               icloudtop =         8000
    640             else
    641               icloudbot(ix,jy,n) = 0
    642               icloudtop =      10000
     506         else ! no cloud
     507            if (rain_cloud_above.eq.1) then ! scavenging
     508               if (lsp.ge.convp) then
     509                  clouds(ix,jy,kz,n)=5 ! lsp dominated washout
     510               else
     511                  clouds(ix,jy,kz,n)=4 ! convp dominated washout
     512               endif
    643513            endif
    644           endif
    645           if (icloudtop .ne. icmv) then
    646             icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
    647           else
    648             icloudthck(ix,jy,n) = icmv
    649           endif
    650 !PS  get rid of too thin clouds     
    651           if (icloudthck(ix,jy,n) .lt. 50) then
    652             icloudbot(ix,jy,n)=icmv
    653             icloudthck(ix,jy,n)=icmv
    654           endif
    655 
    656 
     514         endif
     515      end do
    657516    end do
    658517  end do
    659518
    660   !do 102 kz=1,nuvz
    661   !write(an,'(i02)') kz+10
    662   !write(*,*) nuvz,nymin1,nxmin1,'--',an,'--'
    663   !open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')
    664   !do 101 jy=0,nymin1
    665   !    write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)
    666   !101   continue
    667   ! close(4)
    668   !102   continue
    669 
    670   ! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')
    671   ! do 103 jy=0,nymin1
    672   !     write (4,*)
    673   !+       (height(kz),kz=1,nuvz)
    674   !103   continue
    675   ! close(4)
    676 
    677   !open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')
    678   ! do 104 jy=0,nymin1
    679   !     write (4,*)
    680   !+       (r_air*tt(ix,jy,1,n)*rho(ix,jy,1,n),ix=1,nxmin1)
    681   !104   continue
    682   ! close(4)
    683 
    684 
    685519end subroutine verttransform

  • trunk/src/verttransform_gfs.f90

    r4 r24  
    2121
    2222subroutine verttransform(n,uuh,vvh,wwh,pvh)
    23   !                         i  i   i   i   i
     23  !                      i  i   i   i   i
    2424  !*****************************************************************************
    2525  !                                                                            *
    … …  
    7171  real :: f_qvsat,pressure
    7272  real :: rh,lsp,convp
    73   real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
     73  real :: rhoh(nuvzmax),pinmconv(nzmax)
    7474  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
    7575  real :: xlon,ylat,xlonr,dzdx,dzdy
    76   real :: dzdx1,dzdx2,dzdy1,dzdy2
     76  real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
    7777  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
    7878  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
    … …  
    9292  ! If verttransform is called the first time, initialize heights of the   *
    9393  ! z levels in meter. The heights are the heights of model levels, where  *
    94   ! u,v,T and qv are given, and of the interfaces, where w is given. So,   *
    95   ! the vertical resolution in the z system is doubled. As reference point,*
    96   ! the lower left corner of the grid is used.                             *
    97   ! Unlike in the eta system, no difference between heights for u,v and    *
    98   ! heights for w exists.                                                  *
     94  ! u,v,T and qv are given.                                                *
    9995  !*************************************************************************
    10096
    … …  
    118114
    119115    tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
    120          ps(ixm,jym,1,n))
     116    ps(ixm,jym,1,n))
    121117    pold=ps(ixm,jym,1,n)
    122118    height(1)=0.
    … …  
    126122      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
    127123
    128 
    129   ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
    130   ! upon the transformation to z levels. In order to save computer memory, this is
    131   ! not done anymore in the standard version. However, this option can still be
    132   ! switched on by replacing the following lines with those below, that are
    133   ! currently commented out.
    134   ! Note that two more changes are necessary in this subroutine below.
    135   ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
    136   !*****************************************************************************
    137 
    138124      if (abs(tv-tvold).gt.0.2) then
    139         height(kz)= &
    140              height(kz-1)+const*log(pold/pint)* &
    141              (tv-tvold)/log(tv/tvold)
     125        height(kz)=height(kz-1)+const*log(pold/pint)* &
     126        (tv-tvold)/log(tv/tvold)
    142127      else
    143         height(kz)=height(kz-1)+ &
    144              const*log(pold/pint)*tv
     128        height(kz)=height(kz-1)+const*log(pold/pint)*tv
    145129      endif
    146 
    147   ! Switch on following lines to use doubled vertical resolution
    148   !*************************************************************
    149   !    if (abs(tv-tvold).gt.0.2) then
    150   !      height((kz-1)*2)=
    151   !    +      height(max((kz-2)*2,1))+const*log(pold/pint)*
    152   !    +      (tv-tvold)/log(tv/tvold)
    153   !    else
    154   !      height((kz-1)*2)=height(max((kz-2)*2,1))+
    155   !    +      const*log(pold/pint)*tv
    156   !    endif
    157   ! End doubled vertical resolution
    158130
    159131      tvold=tv
    160132      pold=pint
    161133    end do
    162 
    163 
    164   ! Switch on following lines to use doubled vertical resolution
    165   !*************************************************************
    166   !  do 7 kz=3,nz-1,2
    167   !    height(kz)=0.5*(height(kz-1)+height(kz+1))
    168   !  height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
    169   ! End doubled vertical resolution
    170134
    171135
    … …  
    198162      llev = 0
    199163      do i=1,nuvz
    200        if (ps(ix,jy,1,n).lt.akz(i)) llev=i
     164        if (ps(ix,jy,1,n).lt.akz(i)) llev=i
    201165      end do
    202166       llev = llev+1
    … …  
    212176      tvold=tth(ix,jy,llev,n)*(1.+0.608*qvh(ix,jy,llev,n))
    213177      pold=akz(llev)
    214       uvzlev(llev)=0.
    215178      wzlev(llev)=0.
    216179      uvwzlev(ix,jy,llev)=0.
    … …  
    223186
    224187        if (abs(tv-tvold).gt.0.2) then
    225           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
    226                (tv-tvold)/log(tv/tvold)
     188          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
     189          (tv-tvold)/log(tv/tvold)
    227190        else
    228           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
     191          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
    229192        endif
    230         wzlev(kz)=uvzlev(kz)
    231         uvwzlev(ix,jy,kz)=uvzlev(kz)
     193        wzlev(kz)=uvwzlev(ix,jy,kz)
    232194
    233195        tvold=tv
    234196        pold=pint
    235197      end do
    236 
    237 
    238   ! Switch on following lines to use doubled vertical resolution
    239   ! Switch off the three lines above.
    240   !*************************************************************
    241   !22          uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
    242   !     do 23 kz=2,nwz
    243   !23          uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
    244   ! End doubled vertical resolution
    245198
    246199  ! pinmconv=(h2-h1)/(p2-p1)
    … …  
    279232      do iz=2,nz-1
    280233        do kz=kmin,nuvz
    281           if(height(iz).gt.uvzlev(nuvz)) then
     234          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
    282235            uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
    283236            vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
    … …  
    289242            goto 30
    290243          endif
    291           if ((height(iz).gt.uvzlev(kz-1)).and. &
    292                (height(iz).le.uvzlev(kz))) then
    293            dz1=height(iz)-uvzlev(kz-1)
    294            dz2=uvzlev(kz)-height(iz)
    295            dz=dz1+dz2
    296            uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
    297            vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
    298            tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
    299                 +tth(ix,jy,kz,n)*dz1)/dz
    300            qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
    301                 +qvh(ix,jy,kz,n)*dz1)/dz
    302            pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
    303            rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
    304            pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz
     244          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
     245          (height(iz).le.uvwzlev(ix,jy,kz))) then
     246            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
     247            dz2=uvwzlev(ix,jy,kz)-height(iz)
     248            dz=dz1+dz2
     249            uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
     250            vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
     251            tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
     252            +tth(ix,jy,kz,n)*dz1)/dz
     253            qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
     254            +qvh(ix,jy,kz,n)*dz1)/dz
     255            pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
     256            rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
     257            pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz
    305258          endif
    306259        end do
    … …  
    318271        do kz=kmin,nwz
    319272          if ((height(iz).gt.wzlev(kz-1)).and. &
    320                (height(iz).le.wzlev(kz))) then
    321            dz1=height(iz)-wzlev(kz-1)
    322            dz2=wzlev(kz)-height(iz)
    323            dz=dz1+dz2
    324            ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
    325                 +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
    326 
     273          (height(iz).le.wzlev(kz))) then
     274            dz1=height(iz)-wzlev(kz-1)
     275            dz2=wzlev(kz)-height(iz)
     276            dz=dz1+dz2
     277            ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
     278            +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
    327279          endif
    328280        end do
    … …  
    337289      do kz=2,nz-1
    338290        drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
    339              (height(kz+1)-height(kz-1))
     291        (height(kz+1)-height(kz-1))
    340292      end do
    341293      drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
    … …  
    351303
    352304  do jy=1,ny-2
     305    cosf=cos((real(jy)*dy+ylat0)*pi180)
    353306    do ix=1,nx-2
    354307
    … …  
    367320      do iz=2,nz-1
    368321
    369         ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
     322        ui=uu(ix,jy,iz,n)*dxconst/cosf
    370323        vi=vv(ix,jy,iz,n)*dyconst
    371324
    372325        do kz=kmin,nz
    373326          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
    374                (height(iz).le.uvwzlev(ix,jy,kz))) then
     327          (height(iz).le.uvwzlev(ix,jy,kz))) then
    375328            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
    376329            dz2=uvwzlev(ix,jy,kz)-height(iz)
    … …  
    426379      xlon=xlon0+real(nx/2-1)*dx
    427380      xlonr=xlon*pi/180.
    428       ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
    429            vv(nx/2-1,nymin1,iz,n)**2)
    430       if(vv(nx/2-1,nymin1,iz,n).lt.0.) then
    431         ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
    432              vv(nx/2-1,nymin1,iz,n))-xlonr
     381      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2)
     382      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
     383        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
    433384      elseif (vv(nx/2-1,nymin1,iz,n).gt.0.) then
    434385        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
    435              vv(nx/2-1,nymin1,iz,n))-xlonr
     386        vv(nx/2-1,nymin1,iz,n))-xlonr
    436387      else
    437388        ddpol=pi/2-xlonr
    … …  
    446397      uuaux=-ffpol*sin(xlonr+ddpol)
    447398      vvaux=-ffpol*cos(xlonr+ddpol)
    448       call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
    449            vvpolaux)
    450 
     399      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
    451400      jy=nymin1
    452401      do ix=0,nxmin1
    … …  
    460409  ! ward parallel of latitude
    461410
    462   do iz=1,nz
     411    do iz=1,nz
    463412      wdummy=0.
    464413      jy=ny-2
    … …  
    471420        ww(ix,jy,iz,n)=wdummy
    472421      end do
    473   end do
     422    end do
    474423
    475424  endif
    … …  
    487436        do iz=1,nz
    488437          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
    489                vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
    490                vvpol(ix,jy,iz,n))
     438          vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
    491439        end do
    492440      end do
    … …  
    498446      xlon=xlon0+real(nx/2-1)*dx
    499447      xlonr=xlon*pi/180.
    500       ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
    501            vv(nx/2-1,0,iz,n)**2)
     448      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2)
    502449      if(vv(nx/2-1,0,iz,n).lt.0.) then
    503         ddpol=atan(uu(nx/2-1,0,iz,n)/ &
    504              vv(nx/2-1,0,iz,n))+xlonr
     450        ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
    505451      elseif (vv(nx/2-1,0,iz,n).gt.0.) then
    506         ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
    507              vv(nx/2-1,0,iz,n))-xlonr
     452        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))-xlonr
    508453      else
    509454        ddpol=pi/2-xlonr
    … …  
    518463      uuaux=+ffpol*sin(xlonr-ddpol)
    519464      vvaux=-ffpol*cos(xlonr-ddpol)
    520       call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
    521            vvpolaux)
     465      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
    522466
    523467      jy=0
    … …  
    562506         clouds(ix,jy,kz,n)=0
    563507         if (rh.gt.0.8) then ! in cloud
    564             if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
    565                rain_cloud_above=1
    566                cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
    567                     height(kz)-height(kz-1)
    568                if (lsp.ge.convp) then
    569                   clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
    570                else
    571                   clouds(ix,jy,kz,n)=2 ! convp dominated rainout
    572                endif
    573             else ! no precipitation
    574                   clouds(ix,jy,kz,n)=1 ! cloud
    575             endif
     508           if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
     509              rain_cloud_above=1
     510              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
     511              if (lsp.ge.convp) then
     512                 clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
     513              else
     514                 clouds(ix,jy,kz,n)=2 ! convp dominated rainout
     515              endif
     516           else ! no precipitation
     517             clouds(ix,jy,kz,n)=1 ! cloud
     518           endif
    576519         else ! no cloud
    577             if (rain_cloud_above.eq.1) then ! scavenging
    578                if (lsp.ge.convp) then
    579                   clouds(ix,jy,kz,n)=5 ! lsp dominated washout
    580                else
    581                   clouds(ix,jy,kz,n)=4 ! convp dominated washout
    582                endif
    583             endif
     520           if (rain_cloud_above.eq.1) then ! scavenging
     521             if (lsp.ge.convp) then
     522               clouds(ix,jy,kz,n)=5 ! lsp dominated washout
     523             else
     524               clouds(ix,jy,kz,n)=4 ! convp dominated washout
     525             endif
     526           endif
    584527         endif
    585528      end do

  • trunk/src/verttransform_nests.f90

    r4 r24  
    2121
    2222subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn)
    23   !                               i   i    i    i   i
     23  !                            i   i    i    i   i
    2424  !*****************************************************************************
    2525  !                                                                            *
    … …  
    7272  integer :: rain_cloud_above,kz_inv
    7373  real :: f_qvsat,pressure,rh,lsp,convp
    74   real :: uvzlev(nuvzmax),wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax)
     74  real :: wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax)
    7575  real :: uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax)
    76   real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
     76  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi,cosf
    7777  real :: dzdx,dzdy
    7878  real :: dzdx1,dzdx2,dzdy1,dzdy2
    … …  
    9696
    9797      tvold=tt2n(ix,jy,1,n,l)*(1.+0.378*ew(td2n(ix,jy,1,n,l))/ &
    98            psn(ix,jy,1,n,l))
     98      psn(ix,jy,1,n,l))
    9999      pold=psn(ix,jy,1,n,l)
    100       uvzlev(1)=0.
     100      uvwzlev(ix,jy,1)=0.
    101101      wzlev(1)=0.
    102102      rhoh(1)=pold/(r_air*tvold)
    … …  
    112112
    113113        if (abs(tv-tvold).gt.0.2) then
    114           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
    115                (tv-tvold)/log(tv/tvold)
     114          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
     115          (tv-tvold)/log(tv/tvold)
    116116        else
    117           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
     117          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
    118118        endif
    119119
    … …  
    124124
    125125      do kz=2,nwz-1
    126         wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
    127       end do
    128       wzlev(nwz)=wzlev(nwz-1)+ &
    129            uvzlev(nuvz)-uvzlev(nuvz-1)
    130 
    131   ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
    132   ! upon the transformation to z levels. In order to save computer memory, this is
    133   ! not done anymore in the standard version. However, this option can still be
    134   ! switched on by replacing the following lines with those below, that are
    135   ! currently commented out.
    136   ! Note that one change is also necessary in gridcheck.f,
    137   ! and three changes in verttransform.f
    138   !*****************************************************************************
    139       uvwzlev(ix,jy,1)=0.0
    140       do kz=2,nuvz
    141         uvwzlev(ix,jy,kz)=uvzlev(kz)
    142       end do
    143 
    144   ! Switch on following lines to use doubled vertical resolution
    145   ! Switch off the three lines above.
    146   !*************************************************************
    147   !22          uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
    148   !     do 23 kz=2,nwz
    149   !23          uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
    150   ! End doubled vertical resolution
     126        wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2.
     127      end do
     128      wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1)
    151129
    152130  ! pinmconv=(h2-h1)/(p2-p1)
    153131
    154132      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
    155            ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
    156            (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
     133      ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
     134      (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
    157135      do kz=2,nz-1
    158136        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
    159              ((aknew(kz+1)+bknew(kz+1)*psn(ix,jy,1,n,l))- &
    160              (aknew(kz-1)+bknew(kz-1)*psn(ix,jy,1,n,l)))
     137        ((aknew(kz+1)+bknew(kz+1)*psn(ix,jy,1,n,l))- &
     138        (aknew(kz-1)+bknew(kz-1)*psn(ix,jy,1,n,l)))
    161139      end do
    162140      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
    163            ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
    164            (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
     141      ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
     142      (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
    165143
    166144
    … …  
    183161      do iz=2,nz-1
    184162        do kz=kmin,nuvz
    185           if(height(iz).gt.uvzlev(nuvz)) then
     163          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
    186164            uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
    187165            vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
    … …  
    192170            goto 30
    193171          endif
    194           if ((height(iz).gt.uvzlev(kz-1)).and. &
    195                (height(iz).le.uvzlev(kz))) then
    196            dz1=height(iz)-uvzlev(kz-1)
    197            dz2=uvzlev(kz)-height(iz)
     172          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
     173          (height(iz).le.uvwzlev(ix,jy,kz))) then
     174           dz1=height(iz)-uvwzlev(ix,jy,kz-1)
     175           dz2=uvwzlev(ix,jy,kz)-height(iz)
    198176           dz=dz1+dz2
    199177           uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
    200                 uuhn(ix,jy,kz,l)*dz1)/dz
     178           uuhn(ix,jy,kz,l)*dz1)/dz
    201179           vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
    202                 vvhn(ix,jy,kz,l)*dz1)/dz
     180           vvhn(ix,jy,kz,l)*dz1)/dz
    203181           ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
    204                 tthn(ix,jy,kz,n,l)*dz1)/dz
     182           tthn(ix,jy,kz,n,l)*dz1)/dz
    205183           qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
    206                 qvhn(ix,jy,kz,n,l)*dz1)/dz
     184           qvhn(ix,jy,kz,n,l)*dz1)/dz
    207185           pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
    208                 pvhn(ix,jy,kz,l)*dz1)/dz
     186           pvhn(ix,jy,kz,l)*dz1)/dz
    209187           rhon(ix,jy,iz,n,l)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
    210188           kmin=kz
    … …  
    225203        do kz=kmin,nwz
    226204          if ((height(iz).gt.wzlev(kz-1)).and. &
    227                (height(iz).le.wzlev(kz))) then
    228            dz1=height(iz)-wzlev(kz-1)
    229            dz2=wzlev(kz)-height(iz)
    230            dz=dz1+dz2
    231            wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
    232                 +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
    233            kmin=kz
    234            goto 40
     205          (height(iz).le.wzlev(kz))) then
     206            dz1=height(iz)-wzlev(kz-1)
     207            dz2=wzlev(kz)-height(iz)
     208            dz=dz1+dz2
     209            wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
     210            +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
     211            kmin=kz
     212            goto 40
    235213          endif
    236214        end do
    … …  
    242220
    243221      drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
    244            (height(2)-height(1))
     222      (height(2)-height(1))
    245223      do kz=2,nz-1
    246224        drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
    247              rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
     225        rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
    248226      end do
    249227      drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
    … …  
    259237
    260238  do jy=1,nyn(l)-2
     239    cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
    261240    do ix=1,nxn(l)-2
    262241
    … …  
    264243      do iz=2,nz-1
    265244
    266         ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/ &
    267              cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
     245        ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf
    268246        vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)
    269247
    … …  
    319297               rain_cloud_above=1
    320298               cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ &
    321                     height(kz)-height(kz-1)
     299               height(kz)-height(kz-1)
    322300               if (lsp.ge.convp) then
    323301                  cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout

  • trunk/src/wetdepo.f90

    r20 r24  
    2121
    2222subroutine wetdepo(itime,ltsample,loutnext)
    23   !                     i      i        i
     23  !                  i      i        i
    2424  !*****************************************************************************
    2525  !                                                                            *
    … …  
    3939  ! Code may not be correct for decay of deposition!                           *
    4040  !                                                                            *
    41   ! Modification by Sabine Eckhart to introduce a new in-/below-cloud
    42   ! scheme, not dated
    43   ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
    4441  !*****************************************************************************
    4542  !                                                                            *
    … …  
    7471
    7572  integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
    76   integer :: ngrid,itage,nage,hz,il,interp_time, n , clouds_v !NIK scheme
    77   integer :: kz !PS scheme
    78   integer :: ks, kp, n1,n2, icbot,ictop, indcloud
    79   integer :: scheme_number ! NIK==1, PS ==2
     73  integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
     74  integer :: ks, kp
    8075  real :: S_i, act_temp, cl, cle ! in cloud scavenging
    8176  real :: clouds_h ! cloud height for the specific grid point
    82   real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav, precsub,f
     77  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
    8378  real :: wetdeposit(maxspec),restmass
    8479  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
    … …  
    10398
    10499  do jpart=1,numpart
     100
    105101    if (itra1(jpart).eq.-999999999) goto 20
    106102    if(ldirect.eq.1)then
    … …  
    114110      if (itage.lt.lage(nage)) goto 33
    115111    end do
    116 33   continue
     11233  continue
    117113
    118114
    … …  
    123119    do j=numbnests,1,-1
    124120      if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
    125            (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
     121      (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
    126122        ngrid=j
    127123        goto 23
    128124      endif
    129125    end do
    130 23   continue
     12623  continue
    131127
    132128
    … …  
    151147    interp_time=nint(itime-0.5*ltsample)
    152148
    153 !    PS nest case still needs to be implemented!!   
    154 !    if (ngrid.eq.0) then
    155 !      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
    156 !           1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
    157 !           memtime(1),memtime(2),interp_time,lsp,convp,cc)
    158            call interpol_rain(lsprec,convprec,tcc,     &
    159             icloudbot,icloudthck,nxmax,nymax,1,nx,ny,         &
    160             memind,sngl(xtra1(jpart)),sngl(ytra1(jpart)),1,memtime(1), &
    161             memtime(2),interp_time,lsp,convp,cc,icbot,ictop,icmv) 
    162 !    else
    163 !      call interpol_rain_nests(lsprecn,convprecn,tccn, &
    164 !           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
    165 !           memtime(1),memtime(2),interp_time,lsp,convp,cc)
    166 !    endif
    167 
    168 
    169  !   if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
    170  !PS 2012: subtract a small value, eg 0.01 mm/h, to remove spurious precip
    171         prec = lsp+convp
    172         precsub = 0.01
    173         if (prec .lt. precsub) then
    174           goto 20
    175         else
    176           f = (prec-precsub)/prec
    177           lsp = f*lsp
    178           convp = f*convp
    179         endif
    180 
    181 
     149    if (ngrid.eq.0) then
     150      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
     151      1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
     152      memtime(1),memtime(2),interp_time,lsp,convp,cc)
     153    else
     154      call interpol_rain_nests(lsprecn,convprecn,tccn, &
     155      nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
     156      memtime(1),memtime(2),interp_time,lsp,convp,cc)
     157    endif
     158
     159    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
    182160  ! get the level were the actual particle is in
    183       do il=2,nz
    184         if (height(il).gt.ztra1(jpart)) then
    185           !hz=il-1
    186           kz=il-1
    187           goto 26
    188         endif
    189       end do
    190 26     continue
    191 
    192   n=memind(2)
    193   if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
    194        n=memind(1)
     161    do il=2,nz
     162      if (height(il).gt.ztra1(jpart)) then
     163        hz=il-1
     164        goto 26
     165      endif
     166    end do
     16726  continue
     168
     169    n=memind(2)
     170    if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
     171    n=memind(1)
    195172
    196173  ! if there is no precipitation or the particle is above the clouds no
    197174  ! scavenging is done
    198175
    199 !old scheme
    200 !  if (ngrid.eq.0) then
    201 !     clouds_v=clouds(ix,jy,hz,n)
    202 !     clouds_h=cloudsh(ix,jy,n)
    203 !  else
    204 !     clouds_v=cloudsn(ix,jy,hz,n,ngrid)
    205 !     clouds_h=cloudsnh(ix,jy,n,ngrid)
    206 !  endif
    207 !  !write(*,*) 'there is
    208 !  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
    209 !  if (clouds_v.le.1) goto 20
    210 !  !write (*,*) 'there is scavenging'
    211 
    212   ! PS: part of 2011/2012 fix
    213 
    214         if (ztra1(jpart) .le. float(ictop)) then
    215           if (ztra1(jpart) .gt. float(icbot)) then
    216             indcloud = 2 ! in-cloud
    217           else
    218             indcloud = 1 ! below-cloud
    219           endif
    220         elseif (ictop .eq. icmv) then
    221           indcloud = 0 ! no cloud found, use old scheme
    222         else
    223           goto 20 ! above cloud
    224         endif
    225 
     176    if (ngrid.eq.0) then
     177      clouds_v=clouds(ix,jy,hz,n)
     178      clouds_h=cloudsh(ix,jy,n)
     179    else
     180      clouds_v=cloudsn(ix,jy,hz,n,ngrid)
     181      clouds_h=cloudsnh(ix,jy,n,ngrid)
     182    endif
     183  !write(*,*) 'there is
     184  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
     185    if (clouds_v.le.1) goto 20
     186  !write (*,*) 'there is scavenging'
    226187
    227188  ! 1) Parameterization of the the area fraction of the grid cell where the
    … …  
    267228  !**********************************************************
    268229
    269     do ks=1,nspec            ! loop over species
    270       wetdeposit(ks)=0.
    271 
    272    
    273     !conflicting changes to the same routine: 1=NIK 2 =PS
    274     scheme_number=2   
    275     if (scheme_number.eq.1) then !NIK
    276 
    277       if (weta(ks).gt.0.) then
    278         if (clouds_v.ge.4) then
    279   !          BELOW CLOUD SCAVENGING
    280   !          for aerosols and not highliy soluble substances weta=5E-6
    281           wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
    282   !        write(*,*) 'bel. wetscav: ',wetscav
    283 
    284         else ! below_cloud clouds_v is lt 4 and gt 1 -> in cloud scavenging
    285   !        IN CLOUD SCAVENGING
    286   ! BUGFIX tt for nested fields should be ttn
    287   ! sec may 2008
    288           if (ngrid.gt.0) then
    289              act_temp=ttn(ix,jy,hz,n,ngrid)
    290           else
    291              act_temp=tt(ix,jy,hz,n)
    292           endif
     230    do ks=1,nspec      ! loop over species
     231      wetdeposit(ks)=0.
     232      wetscav=0.   
     233
     234  ! NIK 09.12.13: allowed to turn off either below cloud or in-cloud by including TWO separate blocks of if tests
     235
     236
     237  !****i*******************
     238  ! BELOW CLOUD SCAVENGING
     239  !*********************** 
     240      if (weta(ks).gt.0. .and. clouds_v.ge.4) then !if positive below-cloud coefficient given in SPECIES file, and if in below cloud regime
     241 !       for aerosols and not highliy soluble substances weta=5E-6
     242        wetscav=weta(ks)*prec**wetb(ks)  ! scavenging coefficient
     243      endif !end below-cloud
     244
     245
     246  !********************
     247  ! IN CLOUD SCAVENGING
     248  !********************
     249      if (weta_in(ks).gt.0. .and. clouds_v.lt.4) then !if positive in-cloud coefficient given in SPECIES file, and if in in-cloud regime
     250
     251        if (ngrid.gt.0) then
     252          act_temp=ttn(ix,jy,hz,n,ngrid)
     253        else
     254          act_temp=tt(ix,jy,hz,n)
     255        endif
    293256
    294257! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening
    … …  
    297260! wetc_in=0.9 (default)
    298261! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling)
    299           cl=weta_in(ks)*prec**wetb_in(ks)
    300           if (dquer(ks).gt.0) then ! is particle
    301             S_i=wetc_in(ks)/cl
    302            else ! is gas
    303             cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
    304             S_i=1/cle
    305            endif
    306            wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
    307   !         write(*,*) 'in. wetscav:'
    308   !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h
     262        cl=weta_in(ks)*prec**wetb_in(ks)
     263        if (dquer(ks).gt.0) then ! is particle
     264          S_i=wetc_in(ks)/cl
     265        else ! is gas
     266          cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
     267          S_i=1/cle
    309268        endif
    310 
    311 
    312   !      if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
    313   !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
     269        wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
     270  !     write(*,*) 'in. wetscav:'
     271  !     +  ,wetscav,cle,cl,act_temp,prec,clouds_h
     272       
     273      endif ! end in-cloud
     274
     275
     276     
     277  !**************************************************
     278  ! CALCULATE DEPOSITION
     279  !**************************************************
     280  !     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
     281  !     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
     282
     283      if (wetscav.gt.0.) then
    314284        wetdeposit(ks)=xmass1(jpart,ks)* &
    315              (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! wet deposition
    316   !      new particle mass:
    317   !      if (wetdeposit(ks).gt.0) then
    318   !         write(*,*) 'wetdepo: ',wetdeposit(ks),ks
    319   !      endif
    320         restmass = xmass1(jpart,ks)-wetdeposit(ks)
    321          if (ioutputforeachrelease.eq.1) then
    322            kp=npoint(jpart)
    323          else
    324            kp=1
    325          endif
    326         if (restmass .gt. smallnum) then
    327           xmass1(jpart,ks)=restmass
    328   !ccccccccccccccc depostatistic
    329   !       wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
    330   !ccccccccccccccc depostatistic
    331         else
    332           xmass1(jpart,ks)=0.
    333         endif
    334   ! Correct deposited mass to the last time step when radioactive decay of
    335   ! gridded deposited mass was calculated
    336        if (decay(ks).gt.0.) then
    337           wetdeposit(ks)=wetdeposit(ks) &
    338                *exp(abs(ldeltat)*decay(ks))
    339        endif
    340       else  ! weta(k)
    341          wetdeposit(ks)=0.
    342       endif ! weta(k)
    343 
    344      elseif (scheme_number.eq.2) then ! PS
    345 
    346 !PS          indcloud=0 ! Use this for FOR TESTING,
    347 !PS          will skip the new in/below cloud method !!!             
    348 
    349           if (weta(ks).gt.0.) then
    350             if (indcloud .eq. 1) then ! BELOW CLOUD SCAVENGING
    351 !C               for aerosols and not highliy soluble substances weta=5E-6
    352               wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
    353 !c             write(*,*) 'bel. wetscav: ',wetscav
    354             elseif (indcloud .eq. 2) then !  IN CLOUD SCAVENGING
    355               if (ngrid.gt.0) then
    356                  act_temp=ttn(ix,jy,kz,n,ngrid)
    357               else
    358                  act_temp=tt(ix,jy,kz,n)
    359               endif
    360 
    361 ! from NIK             
    362 ! weta_in=2.0E-07 (default)
    363 ! wetb_in=0.36 (default)
    364 ! wetc_in=0.9 (default)
    365 
    366 
    367               cl=2E-7*prec**0.36
    368               if (dquer(ks).gt.0) then ! is particle
    369                 S_i=0.9/cl
    370               else ! is gas
    371                 cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
    372                 S_i=1/cle
    373               endif
    374               wetscav=S_i*prec/3.6E6/(ictop-icbot) ! 3.6e6 converts mm/h to m/s
    375             else ! PS: no cloud diagnosed, old scheme,
    376 !CPS          using with fixed a,b for simplicity, one may wish to change!!
    377               wetscav = 1.e-4*prec**0.62
    378             endif
    379 
    380 
    381             wetdeposit(ks)=xmass1(jpart,ks)*    &
    382             ! (1.-exp(-wetscav*abs(ltsample)))*fraction  ! wet deposition
    383              (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! fraction = grfraction (PS)
    384             restmass = xmass1(jpart,ks)-wetdeposit(ks)
    385             if (ioutputforeachrelease.eq.1) then
    386               kp=npoint(jpart)
    387             else
    388               kp=1
    389             endif
    390             if (restmass .gt. smallnum) then
    391               xmass1(jpart,ks)=restmass
    392 !cccccccccccccccc depostatistic
    393 !c            wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
    394 !cccccccccccccccc depostatistic
    395             else
    396               xmass1(jpart,ks)=0.
    397             endif
    398 !C Correct deposited mass to the last time step when radioactive decay of
    399 !C gridded deposited mass was calculated
    400             if (decay(ks).gt.0.) then
    401               wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
    402             endif
    403           else  ! weta(k)<0
    404              wetdeposit(ks)=0.
    405           endif
    406 
    407      endif !on scheme
    408 
    409 
    410 
    411     end do
    412 
    413   ! Sabine Eckhard, June 2008 create deposition runs only for forward runs
     285        (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! wet deposition
     286      else ! if no scavenging
     287        wetdeposit(ks)=0.
     288      endif
     289
     290
     291      restmass = xmass1(jpart,ks)-wetdeposit(ks)
     292      if (ioutputforeachrelease.eq.1) then
     293        kp=npoint(jpart)
     294      else
     295        kp=1
     296      endif
     297      if (restmass .gt. smallnum) then
     298        xmass1(jpart,ks)=restmass
     299  !   depostatistic
     300  !   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
     301  !   depostatistic
     302      else
     303        xmass1(jpart,ks)=0.
     304      endif
     305  !   Correct deposited mass to the last time step when radioactive decay of
     306  !   gridded deposited mass was calculated
     307      if (decay(ks).gt.0.) then
     308        wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
     309      endif
     310
     311
     312
     313    end do !all species
     314
     315  ! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
    414316  ! Add the wet deposition to accumulated amount on output grid and nested output grid
    415317  !*****************************************************************************
    416318
    417 !    if (ldirect.eq.1) then
    418 !    call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
    419 !         real(ytra1(jpart)),nage,kp)
    420 !    if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
    421 !         wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), &
    422 !         nage,kp)
    423 !    endif
    424 
    425    !PS
    426         if (ldirect.eq.1) then
    427           call wetdepokernel(nclass(jpart),wetdeposit,     &
    428              sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp)
    429           if (nested_output.eq.1)                          &
    430            call wetdepokernel_nest(nclass(jpart),wetdeposit, &
    431              sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp) 
    432         endif
    433 
     319    if (ldirect.eq.1) then
     320      call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
     321      real(ytra1(jpart)),nage,kp)
     322      if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
     323        wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
     324    endif
    434325
    43532620  continue
    436   end do
     327  end do ! all particles
    437328
    438329end subroutine wetdepo

  • trunk/src/writeheader_nest.f90

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