Index: /trunk/flexpart90.info
===================================================================
--- /trunk/flexpart90.info	(revision 20)
+++ /trunk/flexpart90.info	(revision 20)
@@ -0,0 +1,1 @@
+...is not flexpart77.
Index: /trunk/options/COMMAND
===================================================================
--- /trunk/options/COMMAND	(revision 30)
+++ /trunk/options/COMMAND	(revision 20)
@@ -1,198 +1,31 @@
-********************************************************************************
-*                                                                              *
-*      Input file for the Lagrangian particle dispersion model FLEXPART        *
-*                           Please select your options                         *
-*                                                                              *
-********************************************************************************
-
-1. __                3X, I2
-    1       
-   LDIRECT           1 FOR FORWARD SIMULATION, -1 FOR BACKWARD SIMULATION
-
-2. ________ ______   3X, I8, 1X, I6
-   20111210 000000
-   YYYYMMDD HHMISS   BEGINNING DATE OF SIMULATION
-
-3. ________ ______   3X, I8, 1X, I6
-   20111210 120000
-   YYYYMMDD HHMISS   ENDING DATE OF SIMULATION
-
-4. _____             3X, I5
-   10800 
-   SSSSS             OUTPUT EVERY SSSSS SECONDS
-
-5. _____             3X, I5
-   10800 
-   SSSSS             TIME AVERAGE OF OUTPUT (IN SSSSS SECONDS)
-
-6. _____             3X, I5
-     900 
-   SSSSS             SAMPLING RATE OF OUTPUT (IN SSSSS SECONDS)
-
-7. _________         3X, I9
-   999999999
-   SSSSSSSSS         TIME CONSTANT FOR PARTICLE SPLITTING (IN SECONDS)
-
-8. _____             3X, I5
-     900 
-   SSSSS             SYNCHRONISATION INTERVAL OF FLEXPART (IN SECONDS)
-
-9.  ---.--           4X, F6.4
-     -5.0
-    CTL              FACTOR, BY WHICH TIME STEP MUST BE SMALLER THAN TL
-
-10. ---              4X, I3
-      4
-    IFINE            DECREASE OF TIME STEP FOR VERTICAL MOTION BY FACTOR IFINE
-
-11. -                4X, I1
-    3  
-    IOUT             1 CONCENTRATION (RESIDENCE TIME FOR BACKWARD RUNS) OUTPUT, 2 MIXING RATIO OUTPUT, 3 BOTH,4 PLUME TRAJECT., 5=1+4
-
-12. -                4X, I1
-    0  
-    IPOUT            PARTICLE DUMP: 0 NO, 1 EVERY OUTPUT INTERVAL, 2 ONLY AT END
-
-13. _                4X, I1
-    1
-    LSUBGRID         SUBGRID TERRAIN EFFECT PARAMETERIZATION: 1 YES, 0 NO
-
-14. _                4X, I1
-    1
-    LCONVECTION      CONVECTION: 1 YES, 0 NO
-
-15. _                4X, I1
-    0
-    LAGESPECTRA      AGE SPECTRA: 1 YES, 0 NO
-
-16. _                4X, I1
-    0
-    IPIN             CONTINUE SIMULATION WITH DUMPED PARTICLE DATA: 1 YES, 0 NO
-
-17. _                
-    0                4X,I1 
-    IOFR             IOUTPUTFOREACHREL CREATE AN OUPUT FILE FOR EACH RELEASE LOCATION: 1 YES, 0 NO
-
-18. _                4X, I1
-    0
-    IFLUX            CALCULATE FLUXES: 1 YES, 0 NO
-
-19. _                4X, I1
-    0
-    MDOMAINFILL      DOMAIN-FILLING TRAJECTORY OPTION: 1 YES, 0 NO, 2 STRAT. O3 TRACER
-
-20. _                4X, I1
-    1
-    IND_SOURCE       1=MASS UNIT , 2=MASS MIXING RATIO UNIT 
-
-21. _                4X, I1
-    1
-    IND_RECEPTOR     1=MASS UNIT , 2=MASS MIXING RATIO UNIT 
-
-22. _                4X, I1
-    0
-    MQUASILAG        QUASILAGRANGIAN MODE TO TRACK INDIVIDUAL PARTICLES: 1 YES, 0 NO
-
-23. _                4X, I1
-    0
-    NESTED_OUTPUT    SHALL NESTED OUTPUT BE USED? 1 YES, 0 NO
-
-24. _                4X, I1
-    2
-    LINIT_COND       INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT
-
-25. _                4X, I1
-    0
-    SURF_ONLY        IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER
-
-
-1. Simulation direction, 1 for forward, -1 for backward in time 
-	(consult Seibert and Frank, 2004 for backward runs)
-
-2. Beginning date and time of simulation. Must be given in format
-   YYYYMMDD HHMISS, where YYYY is YEAR, MM is MONTH, DD is DAY, HH is HOUR,
-   MI is MINUTE and SS is SECOND. Current version utilizes UTC.
-
-3. Ending date and time of simulation. Same format as 3.
-
-4. Average concentrations are calculated every SSSSS seconds.
-
-5. The average concentrations are time averages of SSSSS seconds
-   duration. If SSSSS is 0, instantaneous concentrations are outputted.
-
-6. The concentrations are sampled every SSSSS seconds to calculate the time
-   average concentration. This period must be shorter than the averaging time.
-
-7. Time constant for particle splitting. Particles are split into two
-   after SSSSS seconds, 2xSSSSS seconds, 4xSSSSS seconds, and so on.
-
-8. All processes are synchronized with this time interval (lsynctime).
-   Therefore, all other time constants must be multiples of this value.
-   Output interval and time average of output must be at least twice lsynctime.
-
-9. CTL must be >1 for time steps shorter than the  Lagrangian time scale
-   If CTL<0, a purely random walk simulation is done
-
-10.IFINE=Reduction factor for time step used for vertical wind
-
-11.IOUT determines how the output shall be made: concentration
-   (ng/m3, Bq/m3), mixing ratio (pptv), or both, or plume trajectory mode,
-   or concentration + plume trajectory mode.
-   In plume trajectory mode, output is in the form of average trajectories.
-
-12.IPOUT determines whether particle positions are outputted (in addition
-   to the gridded concentrations or mixing ratios) or not.
-   0=no output, 1 output every output interval, 2 only at end of the
-   simulation
-
-13.Switch on/off subgridscale terrain parameterization (increase of
-   mixing heights due to subgridscale orographic variations)
-
-14.Switch on/off the convection parameterization
-
-15.Switch on/off the calculation of age spectra: if yes, the file AGECLASSES
-   must be available
-
-16. If IPIN=1, a file "partposit_end" from a previous run must be available in
-    the output directory. Particle positions are read in and previous simulation
-    is continued. If IPIN=0, no particles from a previous run are used
-
-17. IF IOUTPUTFOREACHRELEASE is set to 1, one output field for each location
-    in the RLEASE file is created. For backward calculation this should be
-    set to 1. For forward calculation both possibilities are applicable.
-
-18. If IFLUX is set to 1, fluxes of each species through each of the output
-    boxes are calculated. Six fluxes, corresponding to northward, southward,
-    eastward, westward, upward and downward are calculated for each grid cell of
-    the output grid. The control surfaces are placed in the middle of each
-    output grid cell. If IFLUX is set to 0, no fluxes are determined.
-
-19. If MDOMAINFILL is set to 1, the first box specified in file RELEASES is used
-    as the domain where domain-filling trajectory calculations are to be done.
-    Particles are initialized uniformly distributed (according to the air mass
-    distribution) in that domain at the beginning of the simulation, and are
-    created at the boundaries throughout the simulation period.
-
-20. IND_SOURCE switches between different units for concentrations at the source
-    NOTE that in backward simulations the release of computational particles 
-    takes place at the "receptor" and the sampling of particles at the "source".
-          1=mass units (for bwd-runs = concentration)
-          2=mass mixing ratio units 
-21. IND_RECEPTOR switches between different units for concentrations at the receptor
-          1=mass units (concentrations)
-          2=mass mixing ratio units 
-
-22. MQUASILAG indicates whether particles shall be numbered consecutively (1) or
-    with their release location number (0). The first option allows tracking of
-    individual particles using the partposit output files
-
-23. NESTED_OUTPUT decides whether model output shall be made also for a nested
-    output field (normally with higher resolution)
-
-24. LINIT_COND determines whether, for backward runs only, the sensitivity to initial
-    conditions shall be calculated and written to output files
-    0=no output, 1 or 2 determines in which units the initial conditions are provided.
-
-25. SURF_ONLY: When set to 1, concentration/emission sensitivity is written out only
-    for the surface layer; useful for instance when only footprint emission sensitivity is needed
-    but initial conditions are needed on a full 3-D grid
++++++++++++++ HEADER +++++++++++++++++
++++++++++++++ HEADER +++++++++++++++++
++++++++++++++ HEADER +++++++++++++++++
++++++++++++++ HEADER +++++++++++++++++
++++++++++++++ HEADER +++++++++++++++++
++++++++++++++ HEADER +++++++++++++++++
++++++++++++++ HEADER +++++++++++++++++
+-1
+20070930 180000
+20070931      0
+  3600
+  3600
+   900
+ 9999999
+900   SYNC
+-5.0  CTL
+4     IFINE
+5     IOUT
+0     IPOUT
+1     LSUBGRID
+1     LCONVECTION
+1     LAGESPECTRA
+0     IPIN
+1     IOFR
+0     IFLUX
+0     MDOMAINFILL
+1     IND_SOURCE
+2     IND_RECEPTOR
+0     MQUASILAG
+1     NESTED_OUTPUT
+2     LINIT_COND        INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT
Index: /trunk/options/COMMAND.alternative
===================================================================
--- /trunk/options/COMMAND.alternative	(revision 30)
+++ /trunk/options/COMMAND.alternative	(revision 20)
@@ -30,5 +30,4 @@
 0               NESTED_OUTPUT     SHALL NESTED OUTPUT BE USED? YES, 0 NO
 2               LINIT_COND        INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT
-0               SURF_ONLY         IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER
 
 
@@ -118,6 +117,2 @@
     conditions shall be calculated and written to output files
     0=no output, 1 or 2 determines in which units the initial conditions are provided.
-
-25. SURF_ONLY: When set to 1, concentration/emission sensitivity is written out only
-    for the surface layer; useful for instance when only footprint emission sensitivity is needed
-    but initial conditions are needed on a full 3-D grid
Index: /trunk/options/COMMAND.reference
===================================================================
--- /trunk/options/COMMAND.reference	(revision 30)
+++ /trunk/options/COMMAND.reference	(revision 20)
@@ -11,9 +11,9 @@
 
 2. ________ ______   3X, I8, 1X, I6
-   20110310 000000
+   20040720 000000
    YYYYMMDD HHMISS   BEGINNING DATE OF SIMULATION
 
 3. ________ ______   3X, I8, 1X, I6
-   20110310 120000
+   20040721 120000
    YYYYMMDD HHMISS   ENDING DATE OF SIMULATION
 
@@ -101,8 +101,4 @@
     2
     LINIT_COND       INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT
-
-25. _                4X, I1
-    0
-    SURF_ONLY        IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER
 
 
@@ -193,6 +189,2 @@
     conditions shall be calculated and written to output files
     0=no output, 1 or 2 determines in which units the initial conditions are provided.
-
-25. SURF_ONLY: When set to 1, concentration/emission sensitivity is written out only
-    for the surface layer; useful for instance when only footprint emission sensitivity is needed
-    but initial conditions are needed on a full 3-D grid
Index: /trunk/options/RELEASES
===================================================================
--- /trunk/options/RELEASES	(revision 30)
+++ /trunk/options/RELEASES	(revision 20)
@@ -1,54 +1,868 @@
-*************************************************************************
-*                                                                       *
-*                                                                       *
-*                                                                       *
-*   Input file for the Lagrangian particle dispersion model FLEXPART    *
-*                        Please select your options                     *
-*                                                                       *
-*                                                                       *
-*                                                                       *
-*************************************************************************
-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  1      
-___                        i3    Total number of species emitted
-
-  3
-___                        i3    Index of species in file SPECIES
-
-=========================================================================
-20111210 000001 
-________ ______            i8,1x,i6 Beginning date and time of release
-
-20111210 090000 
-________ ______            i8,1x,i6 Ending date and time of release
-
-   9.4048 
-____.____                  f9.4  Longitude [DEG] of lower left corner
-
-  48.5060
-____.____                  f9.4  Latitude [DEG] of lower left corner
-
-   9.5067
-____.____                  f9.4  Longitude [DEG] of upper right corner
-
-  48.5158
-____.____                  f9.4  Latitude [DEG] of upper right corner
-
-        2
-_________                  i9    1 for m above ground, 2 for m above sea level
-
-    0.00
-_____.___                  f10.3 Lower z-level (in m agl or m asl)
- 
-   10.00
-_____.___                  f10.3 Upper z-level (in m agl or m asl)
- 
-    20000                
-_________                  i9    Total number of particles to be released
-
-1.0000E00
-_.____E__                  e9.4  Total mass emitted
-
-RELEASE_TEST1
-________________________________________   character*40 comment
-+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++++++++++++++++++ HEADER ++++++++++++++++++++
++++++++++++++++++ HEADER ++++++++++++++++++++
++++++++++++++++++ HEADER ++++++++++++++++++++
++++++++++++++++++ HEADER ++++++++++++++++++++
++++++++++++++++++ HEADER ++++++++++++++++++++
++++++++++++++++++ HEADER ++++++++++++++++++++
++++++++++++++++++ HEADER ++++++++++++++++++++
++++++++++++++++++ HEADER ++++++++++++++++++++
++++++++++++++++++ HEADER ++++++++++++++++++++
++++++++++++++++++ HEADER ++++++++++++++++++++
++++++++++++++++++ HEADER ++++++++++++++++++++
+5
+31
+31
+31
+31
+31
+
+20040101 080000
+20040102 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040108 080000
+20040109 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040115 080000
+20040116 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040122 080000
+20040123 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040129 080000
+20040130 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040205 080000
+20040206 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040212 080000
+20040213 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040219 080000
+20040220 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040226 080000
+20040227 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040304 080000
+20040305 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040311 080000
+20040312 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040318 080000
+20040319 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040325 080000
+20040326 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040401 080000
+20040402 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040408 080000
+20040409 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040415 080000
+20040416 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040422 080000
+20040423 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040429 080000
+20040430 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040609 080000
+20040610 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040610 080000
+20040611 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040617 080000
+20040618 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040624 080000
+20040625 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040701 080000
+20040702 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040708 080000
+20040709 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040715 080000
+20040716 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040722 080000
+20040723 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040729 080000
+20040730 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040805 080000
+20040806 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040812 080000
+20040813 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040819 080000
+20040820 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040826 080000
+20040827 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040902 080000
+20040903 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040909 080000
+20040910 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040916 080000
+20040917 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040923 080000
+20040924 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20040930 080000
+20041001 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20041007 080000
+20041008 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20041014 080000
+20041015 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20041021 080000
+20041022 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20041028 080000
+20041029 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20041104 080000
+20041105 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20041111 080000
+20041112 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20041118 080000
+20041119 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20041124 080000
+20041125 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20041125 080000
+20041126 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20041208 080000
+20041209 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20041209 080000
+20041210 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20041216 080000
+20041217 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20041223 080000
+20041224 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
+20041230 080000
+20041231 080000
+   8.2500
+  58.3833
+   8.2500
+  58.3833
+1
+     0.000
+   100.000
+4000
+1
+1
+1
+1
+1
+rel_bj
+
Index: /trunk/options/SPECIES/SPECIES_001
===================================================================
--- /trunk/options/SPECIES/SPECIES_001	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_001	(revision 20)
@@ -7,10 +7,6 @@
 TRACER             Tracer name
 -999.9             Species half life
--9.9E-09           Below cloud scavenging - A
-                   Below cloud scavenging - B
--9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
--9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+-9.9E-09           Wet deposition - A
+                   Wet deposition - B
 -9.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_002
===================================================================
--- /trunk/options/SPECIES/SPECIES_002	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_002	(revision 20)
@@ -7,10 +7,6 @@
 O3                 Tracer name
 -999.9             Species half life
--9.9E-09           Below cloud scavenging - A
-                   Below cloud scavenging - B
--9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
--9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+-9.9E-09           Wet deposition - A
+                   Wet deposition - B
  1.5               Dry deposition (gases) - D
 1.0e-02            Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_003
===================================================================
--- /trunk/options/SPECIES/SPECIES_003	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_003	(revision 20)
@@ -7,10 +7,6 @@
 NO                 Tracer name
 -999.9             Species half life
- 8.0E-06           Below cloud scavenging - A
-0.62               Below cloud scavenging - B
-2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
-0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
-0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
-1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+ 8.0E-06           Wet deposition - A
+0.62               Wet deposition - B
  1.2               Dry deposition (gases) - D
 3.0e-03            Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_004
===================================================================
--- /trunk/options/SPECIES/SPECIES_004	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_004	(revision 20)
@@ -7,10 +7,6 @@
 NO2                Tracer name
 -999.9             Species half life
- 1.0E-05           Below cloud scavenging - A
-0.62               Below cloud scavenging - B
-2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
-0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
-0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
-1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+ 1.0E-05           Wet deposition - A
+0.62               Wet deposition - B
  1.6               Dry deposition (gases) - D
 1.0e-02            Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_005
===================================================================
--- /trunk/options/SPECIES/SPECIES_005	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_005	(revision 20)
@@ -7,10 +7,6 @@
 HNO3               Tracer name
 -999.9             Species half life
- 5.0E-05           Below cloud scavenging - A
-0.62               Below cloud scavenging - B
-2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
-0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
-0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
-1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+ 5.0E-05           Wet deposition - A
+0.62               Wet deposition - B
  1.9               Dry deposition (gases) - D
 1.0e+14            Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_006
===================================================================
--- /trunk/options/SPECIES/SPECIES_006	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_006	(revision 20)
@@ -7,10 +7,6 @@
 HNO2               Tracer name
 -999.9             Species half life
--9.9E-09           Below cloud scavenging - A
-                   Below cloud scavenging - B
--9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
--9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+-9.9E-09           Wet deposition - A
+                   Wet deposition - B
  1.6               Dry deposition (gases) - D
 1.0e+05            Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_007
===================================================================
--- /trunk/options/SPECIES/SPECIES_007	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_007	(revision 20)
@@ -7,10 +7,6 @@
 H2O2               Tracer name
 -999.9             Species half life
- 1.0E-04           Below cloud scavenging - A
-0.62               Below cloud scavenging - B
-2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
-0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
-0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
-1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+ 1.0E-04           Wet deposition - A
+0.62               Wet deposition - B
  1.4               Dry deposition (gases) - D
 1.0e+05            Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_008
===================================================================
--- /trunk/options/SPECIES/SPECIES_008	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_008	(revision 20)
@@ -7,10 +7,6 @@
 SO2                Tracer name
 -999.9             Species half life
--9.9E-09           Below cloud scavenging - A
-0.62               Below cloud scavenging - B
--9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
--9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+-9.9E-09           Wet deposition - A
+0.62               Wet deposition - B
  2.0               Dry deposition (gases) - D
 1.0e+05            Dry deposition (gases) - Henrys const.
@@ -43,5 +39,5 @@
 16       1.463       1.162
 17       1.426       1.141
-18       1.326       1.117
+18       1.326       1.117 
 19       1.225       1.109
 20       1.082       1.095
Index: /trunk/options/SPECIES/SPECIES_009
===================================================================
--- /trunk/options/SPECIES/SPECIES_009	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_009	(revision 20)
@@ -7,10 +7,6 @@
 HCHO               Tracer name
 -999.9             Species half life
--9.9E-09           Below cloud scavenging - A
-                   Below cloud scavenging - B
--9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
--9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+-9.9E-09           Wet deposition - A
+                   Wet deposition - B
  1.3               Dry deposition (gases) - D
 6.0e+03            Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_010
===================================================================
--- /trunk/options/SPECIES/SPECIES_010	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_010	(revision 20)
@@ -7,10 +7,6 @@
 PAN                Tracer name
 -999.9             Species half life
--9.9E-09           Below cloud scavenging - A
-                   Below cloud scavenging - B
--9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
--9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+-9.9E-09           Wet deposition - A
+                   Wet deposition - B
  2.6               Dry deposition (gases) - D
 3.6e+00            Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_011
===================================================================
--- /trunk/options/SPECIES/SPECIES_011	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_011	(revision 20)
@@ -7,10 +7,6 @@
 NH3                Tracer name
 -999.9             Species half life
- 9.9E-05           Below cloud scavenging - A
-0.62               Below cloud scavenging - B
-2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
-0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
-0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
-1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+ 9.9E-05           Wet deposition - A
+0.62               Wet deposition - B
  1.1               Dry deposition (gases) - D
 2.0e+14            Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_012
===================================================================
--- /trunk/options/SPECIES/SPECIES_012	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_012	(revision 20)
@@ -7,10 +7,6 @@
 SO4-aero           Tracer name
 -999.9             Species half life
- 5.0E-06           Below cloud scavenging - A
-0.62               Below cloud scavenging - B
-2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
-0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
-0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
-1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+ 5.0E-06           Wet deposition - A
+0.62               Wet deposition - B
 -9.9               Dry deposition (gases) - D
  1.0E-09           Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_013
===================================================================
--- /trunk/options/SPECIES/SPECIES_013	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_013	(revision 20)
@@ -7,10 +7,6 @@
 NO3-aero           Tracer name
 -999.9             Species half life
- 5.0E-06           Below cloud scavenging - A
-0.62               Below cloud scavenging - B
-2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
-0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
-0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
-1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+ 5.0E-06           Wet deposition - A
+0.62               Wet deposition - B
 -9.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_014
===================================================================
--- /trunk/options/SPECIES/SPECIES_014	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_014	(revision 20)
@@ -7,10 +7,6 @@
 I2-131             Tracer name
 691200.0           Species half life
- 8.0E-05           Below cloud scavenging - A
-0.62               Below cloud scavenging - B
-2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
-0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
-0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
-1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+ 8.0E-05           Wet deposition - A
+0.62               Wet deposition - B
  2.7               Dry deposition (gases) - D
 1.0e+05            Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_015
===================================================================
--- /trunk/options/SPECIES/SPECIES_015	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_015	(revision 20)
@@ -7,10 +7,6 @@
 I-131              Tracer name
 691200.0           Species half life
- 1.0E-04           Below cloud scavenging - A
-0.80               Below cloud scavenging - B
-2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
-0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
-0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
-1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+ 1.0E-04           Wet deposition - A
+0.80               Wet deposition - B
 -9.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_016
===================================================================
--- /trunk/options/SPECIES/SPECIES_016	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_016	(revision 20)
@@ -7,10 +7,6 @@
 Cs-137             Tracer name
 -999.9             Species half life
- 1.0E-04           Below cloud scavenging - A
-0.80               Below cloud scavenging - B
-2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
-0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
-0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
-1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+ 1.0E-04           Wet deposition - A
+0.80               Wet deposition - B
 -9.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_017
===================================================================
--- /trunk/options/SPECIES/SPECIES_017	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_017	(revision 20)
@@ -7,10 +7,6 @@
 Y-91               Tracer name
 5037120.0          Species half life
- 1.0E-04           Below cloud scavenging - A
-0.80               Below cloud scavenging - B
-2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
-0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
-0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
-1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+ 1.0E-04           Wet deposition - A
+0.80               Wet deposition - B
 -9.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_018
===================================================================
--- /trunk/options/SPECIES/SPECIES_018	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_018	(revision 20)
@@ -7,10 +7,6 @@
 Ru-106             Tracer name
 31536000.0         Species half life
- 1.0E-04           Below cloud scavenging - A
-0.80               Below cloud scavenging - B
-2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
-0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
-0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
-1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+ 1.0E-04           Wet deposition - A
+0.80               Wet deposition - B
 -9.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_019
===================================================================
--- /trunk/options/SPECIES/SPECIES_019	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_019	(revision 20)
@@ -7,10 +7,6 @@
 Kr-85              Tracer name
 -999.9             Species half life
--9.9E-09           Below cloud scavenging - A
-                   Below cloud scavenging - B
--9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
--9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+-9.9E-09           Wet deposition - A
+                   Wet deposition - B
 -9.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_020
===================================================================
--- /trunk/options/SPECIES/SPECIES_020	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_020	(revision 20)
@@ -7,10 +7,6 @@
 Sr-90              Tracer name
 -999.9             Species half life
- 1.0E-04           Below cloud scavenging - A
-0.80               Below cloud scavenging - B
-2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
-0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
-0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
-1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+ 1.0E-04           Wet deposition - A
+0.80               Wet deposition - B
 -9.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_021
===================================================================
--- /trunk/options/SPECIES/SPECIES_021	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_021	(revision 20)
@@ -6,11 +6,7 @@
 ****************************************************************************
 Xe-133             Tracer name
-453168.0           Species half life
--9.9E-09           Below cloud scavenging - A
-                   Below cloud scavenging - B
--9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
--9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+198720.0           Species half life
+-9.9E-09           Wet deposition - A
+                   Wet deposition - B
 -9.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_022
===================================================================
--- /trunk/options/SPECIES/SPECIES_022	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_022	(revision 20)
@@ -7,10 +7,6 @@
 CO                 Tracer name
 -999.9             Species half life
--9.9E-09           Below cloud scavenging - A
-                   Below cloud scavenging - B
--9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
--9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+-9.9E-09           Wet deposition - A
+                   Wet deposition - B
 -9.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_024
===================================================================
--- /trunk/options/SPECIES/SPECIES_024	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_024	(revision 20)
@@ -7,10 +7,6 @@
 AIRTRACER             Tracer name
 -999.9             Species half life
--9.9E-09           Below cloud scavenging - A
-                   Below cloud scavenging - B
--9.9E-09           In-cloud scavenging - Ai (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
--9.9               In-cloud scavenging - Ci (S_i=Ci/cl)
--9.9               In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
+-9.9E-09           Wet deposition - A
+                   Wet deposition - B
 -9.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: /trunk/options/SPECIES/SPECIES_025
===================================================================
--- /trunk/options/SPECIES/SPECIES_025	(revision 30)
+++ /trunk/options/SPECIES/SPECIES_025	(revision 20)
@@ -9,8 +9,4 @@
 2.0E-05            Wet deposition - A
 0.8                Wet deposition - B
-2.0E-7             In-cloud scavenging - Ai (cl=Ai*prec**Bi)
-0.36               In-cloud scavenging - Bi (cl=Ai*prec**Bi)
-0.9                In-cloud scavenging - Ci (S_i=Ci/cl)
-1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
 -9.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: /trunk/src/FLEXPART.f90
===================================================================
--- /trunk/src/FLEXPART.f90	(revision 30)
+++ /trunk/src/FLEXPART.f90	(revision 20)
@@ -45,9 +45,4 @@
   use conv_mod
 
-#ifdef NETCDF_OUTPUT
-  use netcdf_output_mod, only: writeheader_netcdf
-#endif
-
-
   implicit none
 
@@ -65,8 +60,7 @@
   call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
 
-  ! FLEXPART version string
-  flexversion='Version 9.2 beta (2014-07-01)'
-  verbosity=0
-
+  ! 
+  flexversion='Version 9.1.8  (2013-12-08)'
+  !verbosity=0
   ! Read the pathnames where input/output files are stored
   !*******************************************************
@@ -82,42 +76,48 @@
     call getarg(1,arg1)
     pathfile=arg1
+    verbosity=0
     if (arg1(1:1).eq.'-') then
-      write(pathfile,'(a11)') './pathnames'
-      inline_options=arg1 
+        write(pathfile,'(a11)') './pathnames'
+        inline_options=arg1 
     endif
   case (0)
     write(pathfile,'(a11)') './pathnames'
+    verbosity=0
   end select
   
+    if (inline_options(1:1).eq.'-') then
+      print*, 'inline options=', inline_options       
+      if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
+         print*, 'verbose mode 1: additional information will be displayed'
+         verbosity=1
+      endif
+      if (trim(inline_options).eq.'-v2') then
+         print*, 'verbose mode 2: additional information will be displayed'
+         verbosity=2
+      endif
+      if (trim(inline_options).eq.'-i') then
+         print*, 'info mode: will provide run specific information and stop'
+         verbosity=1
+         info_flag=1
+      endif
+      if (trim(inline_options).eq.'-i2') then
+         print*, 'info mode: will provide run specific information and stop'
+         verbosity=2
+         info_flag=1
+      endif
+    endif
+
+
   ! Print the GPL License statement
   !*******************************************************
-  print*,'Welcome to FLEXPART ', trim(flexversion)
-  print*,'FLEXPART is free software released under the GNU General Public License.'
- 
-  if (inline_options(1:1).eq.'-') then
-    if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
-       print*, 'Verbose mode 1: display detailed information during run'
-       verbosity=1
-    endif
-    if (trim(inline_options).eq.'-v2') then
-       print*, 'Verbose mode 2: display more detailed information during run'
-       verbosity=2
-    endif
-    if (trim(inline_options).eq.'-i') then
-       print*, 'Info mode: provide detailed run specific information and stop'
-       verbosity=1
-       info_flag=1
-    endif
-    if (trim(inline_options).eq.'-i2') then
-       print*, 'Info mode: provide more detailed run specific information and stop'
-       verbosity=2
-       info_flag=1
-    endif
-  endif
-           
-  if (verbosity.gt.0) then
-    write(*,*) 'call readpaths'
+  print*,'Welcome to FLEXPART', trim(flexversion)
+  print*,'FLEXPART is free software released under the GNU Genera'// &
+       'l Public License.'
+            
+  if (verbosity.gt.0) then
+        WRITE(*,*) 'call readpaths'
   endif 
   call readpaths(pathfile)
+   
  
   if (verbosity.gt.1) then !show clock info 
@@ -131,32 +131,36 @@
   endif
 
+
   ! Read the user specifications for the current model run
   !*******************************************************
 
   if (verbosity.gt.0) then
-    write(*,*) 'call readcommand'
+        WRITE(*,*) 'call readcommand'
   endif
   call readcommand
   if (verbosity.gt.0) then
-    write(*,*) '    ldirect=', ldirect
-    write(*,*) '    ibdate,ibtime=',ibdate,ibtime
-    write(*,*) '    iedate,ietime=', iedate,ietime
-    if (verbosity.gt.1) then   
-      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
-      write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
-    endif     
-  endif
+        WRITE(*,*) '    ldirect=', ldirect
+        WRITE(*,*) '    ibdate,ibtime=',ibdate,ibtime
+        WRITE(*,*) '    iedate,ietime=', iedate,ietime
+        if (verbosity.gt.1) then   
+                CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+                WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+        endif     
+  endif
+
 
   ! Read the age classes to be used
   !********************************
   if (verbosity.gt.0) then
-    write(*,*) 'call readageclasses'
+        WRITE(*,*) 'call readageclasses'
   endif
   call readageclasses
 
   if (verbosity.gt.1) then   
-    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
-    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+                CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+                WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
   endif     
+ 
+
 
   ! Read, which wind fields are available within the modelling period
@@ -164,5 +168,5 @@
 
   if (verbosity.gt.0) then
-    write(*,*) 'call readavailable'
+        WRITE(*,*) 'call readavailable'
   endif  
   call readavailable
@@ -173,5 +177,5 @@
  
   if (verbosity.gt.0) then
-     write(*,*) 'call gridcheck'
+     WRITE(*,*) 'call gridcheck'
   endif
 
@@ -179,18 +183,20 @@
 
   if (verbosity.gt.1) then   
-    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
-    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+     CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+     WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
   endif      
-
-  if (verbosity.gt.0) then
-    write(*,*) 'call gridcheck_nests'
+  
+
+  if (verbosity.gt.0) then
+        WRITE(*,*) 'call gridcheck_nests'
   endif  
   call gridcheck_nests
 
+
   ! Read the output grid specifications
   !************************************
 
   if (verbosity.gt.0) then
-    write(*,*) 'call readoutgrid'
+        WRITE(*,*) 'call readoutgrid'
   endif
 
@@ -198,7 +204,7 @@
 
   if (nested_output.eq.1) then
-    call readoutgrid_nest
+          call readoutgrid_nest
     if (verbosity.gt.0) then
-      write(*,*) '# readoutgrid_nest'
+        WRITE(*,*) '# readoutgrid_nest'
     endif
   endif
@@ -226,4 +232,5 @@
   call readlanduse
 
+
   ! Assign fractional cover of landuse classes to each ECMWF grid point
   !********************************************************************
@@ -234,4 +241,6 @@
   call assignland
 
+
+
   ! Read the coordinates of the release locations
   !**********************************************
@@ -242,4 +251,5 @@
   call readreleases
 
+
   ! Read and compute surface resistances to dry deposition of gases
   !****************************************************************
@@ -257,4 +267,5 @@
     print*,'call coordtrafo'
   endif
+
 
   ! Initialize all particles to non-existent
@@ -284,8 +295,10 @@
   endif
 
+
   ! Calculate volume, surface area, etc., of all output grid cells
   ! Allocate fluxes and OHfield if necessary
   !***************************************************************
 
+
   if (verbosity.gt.0) then
     print*,'call outgrid_init'
@@ -294,4 +307,5 @@
   if (nested_output.eq.1) call outgrid_init_nest
 
+
   ! Read the OH field
   !******************
@@ -299,7 +313,7 @@
   if (OHREA.eqv..TRUE.) then
     if (verbosity.gt.0) then
-      print*,'call readOHfield'
+       print*,'call readOHfield'
     endif
-    call readOHfield
+       call readOHfield
   endif
 
@@ -308,21 +322,4 @@
   !******************************************************************
 
-  if (lnetcdfout.eq.1) then 
-#ifdef NETCDF_OUTPUT
-    call writeheader_netcdf(lnest = .false.)
-#endif
-  else 
-    call writeheader
-  end if
-
-  if (nested_output.eq.1) then
-    if (lnetcdfout.eq.1) then
-#ifdef NETCDF_OUTPUT
-      call writeheader_netcdf(lnest = .true.)
-#endif
-    else
-      call writeheader_nest
-    endif
-  endif
 
   if (verbosity.gt.0) then
@@ -335,10 +332,11 @@
   !if (nested_output.eq.1) call writeheader_nest
   if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
+
   if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
   if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
 
-  if (lnetcdfout.ne.1) then
-    open(unitdates,file=path(2)(1:length(2))//'dates')
-  end if
+
+
+  !open(unitdates,file=path(2)(1:length(2))//'dates')
 
   if (verbosity.gt.0) then
@@ -348,4 +346,5 @@
   if ((iout.eq.4).or.(iout.eq.5)) call openouttraj
 
+
   ! Releases can only start and end at discrete times (multiples of lsynctime)
   !***************************************************************************
@@ -355,7 +354,10 @@
   endif
   do i=1,numpoint
-    ireleasestart(i)=nint(real(ireleasestart(i))/real(lsynctime))*lsynctime
-    ireleaseend(i)=nint(real(ireleaseend(i))/real(lsynctime))*lsynctime
+    ireleasestart(i)=nint(real(ireleasestart(i))/ &
+         real(lsynctime))*lsynctime
+    ireleaseend(i)=nint(real(ireleaseend(i))/ &
+         real(lsynctime))*lsynctime
   end do
+
 
   ! Initialize cloud-base mass fluxes for the convection scheme
@@ -379,4 +381,5 @@
   end do
 
+
   ! Calculate particle trajectories
   !********************************
@@ -384,10 +387,10 @@
   if (verbosity.gt.0) then
      if (verbosity.gt.1) then   
-       call system_clock(count_clock, count_rate, count_max)
-       write(*,*) 'System clock',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+       CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+       WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
      endif
      if (info_flag.eq.1) then
-       print*, 'Info only mode (stop)'    
-       stop
+         print*, 'info only mode (stop)'    
+         stop
      endif
      print*,'call timemanager'
@@ -396,12 +399,7 @@
   call timemanager
 
-  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLEXPART MODEL RUN!'
-
-  ! output wall time
-  if (verbosity .gt. 0) then
-    call system_clock(count_clock,count_rate)
-    tins=(count_clock - count_clock0)/real(count_rate)
-    print*,'Wall time ',tins,'s, ',tins/60,'min, ',tins/3600,'h.'
-  endif
+
+  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
+       &XPART MODEL RUN!'
 
 end program flexpart
Index: /trunk/src/advance.f90
===================================================================
--- /trunk/src/advance.f90	(revision 30)
+++ /trunk/src/advance.f90	(revision 20)
@@ -463,8 +463,4 @@
       dcwsave=dcwsave+vp*dt
       zt=zt+w*dt*real(ldirect)
-
-      ! HSO/AL: Particle managed to go over highest level -> interpolation error in goto 700
-      !          alias interpol_wind (division by zero)
-      if (zt.ge.height(nz)) zt=height(nz)-100.*eps
 
       if (zt.gt.h) then
@@ -703,14 +699,4 @@
   endif
 
-  ! HSO/AL: Prevent particles from disappearing at the pole
-  !******************************************************************
-
-  if ( yt.lt.0. ) then
-    xt=mod(xt+180.,360.)
-    yt=-yt
-  else if ( yt.gt.real(nymin1) ) then
-    xt=mod(xt+180.,360.)
-    yt=2*real(nymin1)-yt
-  endif
 
   ! Check position: If trajectory outside model domain, terminate it
@@ -718,5 +704,5 @@
 
   if ((xt.lt.0.).or.(xt.ge.real(nxmin1)).or.(yt.lt.0.).or. &
-       (yt.gt.real(nymin1))) then
+       (yt.ge.real(nymin1))) then
     nstop=3
     return
@@ -873,20 +859,9 @@
   endif
 
-  ! HSO/AL: Prevent particles from disappearing at the pole
-  !******************************************************************
-
-  if ( yt.lt.0. ) then
-    xt=mod(xt+180.,360.)
-    yt=-yt
-  else if ( yt.gt.real(nymin1) ) then
-    xt=mod(xt+180.,360.)
-    yt=2*real(nymin1)-yt
-  endif
-
   ! Check position: If trajectory outside model domain, terminate it
   !*****************************************************************
 
   if ((xt.lt.0.).or.(xt.ge.real(nxmin1)).or.(yt.lt.0.).or. &
-       (yt.gt.real(nymin1))) then
+       (yt.ge.real(nymin1))) then
     nstop=3
     return
Index: /trunk/src/calcpv.f90
===================================================================
--- /trunk/src/calcpv.f90	(revision 30)
+++ /trunk/src/calcpv.f90	(revision 20)
@@ -21,5 +21,5 @@
 
 subroutine calcpv(n,uuh,vvh,pvh)
-  !               i  i   i   o
+  !                  i  i   i   o
   !*****************************************************************************
   !                                                                            *
@@ -49,6 +49,6 @@
   integer :: nlck
   real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
-  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt
-  real :: pvavr,ppml(0:nxmax-1,0:nymax-1,nuvzmax),ppmk(0:nxmax-1,0:nymax-1,nuvzmax)
+  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
+  real :: pvavr,ppml(nuvzmax)
   real :: thup,thdn
   real,parameter :: eps=1.e-5, p0=101325
@@ -62,13 +62,4 @@
   ! Loop over entire grid
   !**********************
-  do kl=1,nuvz
-    do jy=0,nymin1
-      do ix=0,nxmin1
-         ppml(ix,jy,kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
-      enddo
-    enddo
-  enddo
-  ppmk=(100000./ppml)**kappa
-
   do jy=0,nymin1
     if (sglobal.and.jy.eq.0) goto 10
@@ -115,4 +106,8 @@
       end if
       jux=jumpx
+  ! Precalculate pressure values for efficiency
+      do kl=1,nuvz
+        ppml(kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
+      end do
   !
   ! Loop over the vertical
@@ -120,5 +115,6 @@
 
       do kl=1,nuvz
-        theta=tth(ix,jy,kl,n)*ppmk(ix,jy,kl)
+        ppmk=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
+        theta=tth(ix,jy,kl,n)*(100000./ppmk)**kappa
         klvrp=kl+1
         klvrm=kl-1
@@ -129,7 +125,9 @@
         if (klvrp.gt.nuvz) klvrp=nuvz
         if (klvrm.lt.1) klvrm=1
-        thetap=tth(ix,jy,klvrp,n)*ppmk(ix,jy,klvrp)
-        thetam=tth(ix,jy,klvrm,n)*ppmk(ix,jy,klvrm)
-        dthetadp=(thetap-thetam)/(ppml(ix,jy,klvrp)-ppml(ix,jy,klvrm))
+        ppmk=akz(klvrp)+bkz(klvrp)*ps(ix,jy,1,n)
+        thetap=tth(ix,jy,klvrp,n)*(100000./ppmk)**kappa
+        ppmk=akz(klvrm)+bkz(klvrm)*ps(ix,jy,1,n)
+        thetam=tth(ix,jy,klvrm,n)*(100000./ppmk)**kappa
+        dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm))
 
   ! Compute vertical position at pot. temperature surface on subgrid
@@ -158,21 +156,23 @@
           kch=kch+1
           k=kup
-          thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k)
-          thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1)
-
-
-          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
-          ((thdn.le.theta).and.(thup.ge.theta))) then
-            dt1=abs(theta-thdn)
-            dt2=abs(theta-thup)
-            dt=dt1+dt2
-            if (dt.lt.eps) then   ! Avoid division by zero error
-              dt1=0.5             ! G.W., 10.4.1996
-              dt2=0.5
-              dt=1.0
+          ppmk=akz(k)+bkz(k)*ps(ivr,jy,1,n)
+          thdn=tth(ivr,jy,k,n)*(100000./ppmk)**kappa
+          ppmk=akz(k+1)+bkz(k+1)*ps(ivr,jy,1,n)
+          thup=tth(ivr,jy,k+1,n)*(100000./ppmk)**kappa
+
+
+      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
+           ((thdn.le.theta).and.(thup.ge.theta))) then
+              dt1=abs(theta-thdn)
+              dt2=abs(theta-thup)
+              dt=dt1+dt2
+              if (dt.lt.eps) then   ! Avoid division by zero error
+                dt1=0.5             ! G.W., 10.4.1996
+                dt2=0.5
+                dt=1.0
+              endif
+          vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
+              goto 20
             endif
-            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
-            goto 20
-          endif
 41        continue
   ! Downward branch
@@ -181,21 +181,23 @@
           kch=kch+1
           k=kdn
-          thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k)
-          thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1)
-
-          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
-          ((thdn.le.theta).and.(thup.ge.theta))) then
-            dt1=abs(theta-thdn)
-            dt2=abs(theta-thup)
-            dt=dt1+dt2
-            if (dt.lt.eps) then   ! Avoid division by zero error
-              dt1=0.5             ! G.W., 10.4.1996
-              dt2=0.5
-              dt=1.0
+          ppmk=akz(k)+bkz(k)*ps(ivr,jy,1,n)
+          thdn=tth(ivr,jy,k,n)*(100000./ppmk)**kappa
+          ppmk=akz(k+1)+bkz(k+1)*ps(ivr,jy,1,n)
+          thup=tth(ivr,jy,k+1,n)*(100000./ppmk)**kappa
+
+      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
+           ((thdn.le.theta).and.(thup.ge.theta))) then
+              dt1=abs(theta-thdn)
+              dt2=abs(theta-thup)
+              dt=dt1+dt2
+              if (dt.lt.eps) then   ! Avoid division by zero error
+                dt1=0.5             ! G.W., 10.4.1996
+                dt2=0.5
+                dt=1.0
+              endif
+          vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
+              goto 20
             endif
-            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
-            goto 20
-          endif
-          goto 40
+            goto 40
   ! This section used when no values were found
 21      continue
@@ -234,19 +236,21 @@
           kch=kch+1
           k=kup
-          thdn=tth(ix,j,k,n)*ppmk(ix,j,k)
-          thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1)
-          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
-          ((thdn.le.theta).and.(thup.ge.theta))) then
-            dt1=abs(theta-thdn)
-            dt2=abs(theta-thup)
-            dt=dt1+dt2
-            if (dt.lt.eps) then   ! Avoid division by zero error
-              dt1=0.5             ! G.W., 10.4.1996
-              dt2=0.5
-              dt=1.0
+          ppmk=akz(k)+bkz(k)*ps(ix,j,1,n)
+          thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa
+          ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n)
+          thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa
+      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
+           ((thdn.le.theta).and.(thup.ge.theta))) then
+              dt1=abs(theta-thdn)
+              dt2=abs(theta-thup)
+              dt=dt1+dt2
+              if (dt.lt.eps) then   ! Avoid division by zero error
+                dt1=0.5             ! G.W., 10.4.1996
+                dt2=0.5
+                dt=1.0
+              endif
+              uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
+              goto 50
             endif
-            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
-            goto 50
-          endif
 71        continue
   ! Downward branch
@@ -255,25 +259,27 @@
           kch=kch+1
           k=kdn
-          thdn=tth(ix,j,k,n)*ppmk(ix,j,k)
-          thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1)
-          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
-          ((thdn.le.theta).and.(thup.ge.theta))) then
-            dt1=abs(theta-thdn)
-            dt2=abs(theta-thup)
-            dt=dt1+dt2
-            if (dt.lt.eps) then   ! Avoid division by zero error
-              dt1=0.5             ! G.W., 10.4.1996
-              dt2=0.5
-              dt=1.0
+          ppmk=akz(k)+bkz(k)*ps(ix,j,1,n)
+          thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa
+          ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n)
+          thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa
+      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
+           ((thdn.le.theta).and.(thup.ge.theta))) then
+              dt1=abs(theta-thdn)
+              dt2=abs(theta-thup)
+              dt=dt1+dt2
+              if (dt.lt.eps) then   ! Avoid division by zero error
+                dt1=0.5             ! G.W., 10.4.1996
+                dt2=0.5
+                dt=1.0
+              endif
+              uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
+              goto 50
             endif
-            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
-            goto 50
-          endif
-          goto 70
+            goto 70
   ! This section used when no values were found
 51      continue
   ! Must use uu at current level and lat. juy becomes smaller by 1
-          uy(jj)=uuh(ix,jy,kl)
-          juy=juy-1
+        uy(jj)=uuh(ix,jy,kl)
+        juy=juy-1
   ! Otherwise OK
 50        continue
Index: /trunk/src/calcpv_nests.f90
===================================================================
--- /trunk/src/calcpv_nests.f90	(revision 30)
+++ /trunk/src/calcpv_nests.f90	(revision 20)
@@ -21,5 +21,5 @@
 
 subroutine calcpv_nests(l,n,uuhn,vvhn,pvhn)
-  !                     i i  i    i    o
+  !                        i i  i    i    o
   !*****************************************************************************
   !                                                                            *
@@ -48,6 +48,6 @@
   integer :: nlck,l
   real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
-  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt
-  real :: ppml(0:nxmaxn-1,0:nymaxn-1,nuvzmax),ppmk(0:nxmaxn-1,0:nymaxn-1,nuvzmax)
+  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
+  real :: ppml(nuvzmax)
   real :: thup,thdn
   real,parameter :: eps=1.e-5,p0=101325
@@ -61,13 +61,4 @@
   ! Loop over entire grid
   !**********************
-  do kl=1,nuvz
-    do jy=0,nyn(l)-1
-      do ix=0,nxn(l)-1
-         ppml(ix,jy,kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
-      enddo
-    enddo
-  enddo
-  ppmk=(100000./ppml)**kappa
-
   do jy=0,nyn(l)-1
     phi = (ylat0n(l) + jy * dyn(l)) * pi / 180.
@@ -97,4 +88,8 @@
       if (ix.eq.0.or.ix.eq.nxn(l)-1) jumpx=1
       jux=jumpx
+  ! Precalculate pressure values for efficiency
+      do kl=1,nuvz
+        ppml(kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
+      end do
   !
   ! Loop over the vertical
@@ -102,5 +97,6 @@
 
       do kl=1,nuvz
-        theta=tthn(ix,jy,kl,n,l)*ppmk(ix,jy,kl)
+        ppmk=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
+        theta=tthn(ix,jy,kl,n,l)*(100000./ppmk)**kappa
         klvrp=kl+1
         klvrm=kl-1
@@ -111,7 +107,9 @@
         if (klvrp.gt.nuvz) klvrp=nuvz
         if (klvrm.lt.1) klvrm=1
-        thetap=tthn(ix,jy,klvrp,n,l)*ppmk(ix,jy,klvrp)
-        thetam=tthn(ix,jy,klvrm,n,l)*ppmk(ix,jy,klvrm)
-        dthetadp=(thetap-thetam)/(ppml(ix,jy,klvrp)-ppml(ix,jy,klvrm))
+        ppmk=akz(klvrp)+bkz(klvrp)*psn(ix,jy,1,n,l)
+        thetap=tthn(ix,jy,klvrp,n,l)*(100000./ppmk)**kappa
+        ppmk=akz(klvrm)+bkz(klvrm)*psn(ix,jy,1,n,l)
+        thetam=tthn(ix,jy,klvrm,n,l)*(100000./ppmk)**kappa
+        dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm))
 
   ! Compute vertical position at pot. temperature surface on subgrid
@@ -136,6 +134,8 @@
           kch=kch+1
           k=kup
-          thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k)
-          thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1)
+          ppmk=akz(k)+bkz(k)*psn(ivr,jy,1,n,l)
+          thdn=tthn(ivr,jy,k,n,l)*(100000./ppmk)**kappa
+          ppmk=akz(k+1)+bkz(k+1)*psn(ivr,jy,1,n,l)
+          thup=tthn(ivr,jy,k+1,n,l)*(100000./ppmk)**kappa
 
       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
@@ -158,6 +158,8 @@
           kch=kch+1
           k=kdn
-          thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k)
-          thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1)
+          ppmk=akz(k)+bkz(k)*psn(ivr,jy,1,n,l)
+          thdn=tthn(ivr,jy,k,n,l)*(100000./ppmk)**kappa
+          ppmk=akz(k+1)+bkz(k+1)*psn(ivr,jy,1,n,l)
+          thup=tthn(ivr,jy,k+1,n,l)*(100000./ppmk)**kappa
       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
            ((thdn.le.theta).and.(thup.ge.theta))) then
@@ -210,6 +212,8 @@
           kch=kch+1
           k=kup
-          thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k)
-          thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1)
+          ppmk=akz(k)+bkz(k)*psn(ix,j,1,n,l)
+          thdn=tthn(ix,j,k,n,l)*(100000./ppmk)**kappa
+          ppmk=akz(k+1)+bkz(k+1)*psn(ix,j,1,n,l)
+          thup=tthn(ix,j,k+1,n,l)*(100000./ppmk)**kappa
       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
            ((thdn.le.theta).and.(thup.ge.theta))) then
@@ -231,6 +235,8 @@
           kch=kch+1
           k=kdn
-          thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k)
-          thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1)
+          ppmk=akz(k)+bkz(k)*psn(ix,j,1,n,l)
+          thdn=tthn(ix,j,k,n,l)*(100000./ppmk)**kappa
+          ppmk=akz(k+1)+bkz(k+1)*psn(ix,j,1,n,l)
+          thup=tthn(ix,j,k+1,n,l)*(100000./ppmk)**kappa
       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
            ((thdn.le.theta).and.(thup.ge.theta))) then
Index: /trunk/src/com_mod.f90
===================================================================
--- /trunk/src/com_mod.f90	(revision 30)
+++ /trunk/src/com_mod.f90	(revision 20)
@@ -113,6 +113,4 @@
   ! lagespectra  1 if age spectra calculation switched on, 2 if not
 
-  integer :: lnetcdfout
-  ! lnetcdfout   1 for netcdf grid output, 0 if not. Set in COMMAND (namelist input)
 
   integer :: nageclass,lage(maxageclass)
@@ -342,6 +340,9 @@
   real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2)
   real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2)
+  !scavenging NIK, PS
   integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2)
   integer :: cloudsh(0:nxmax-1,0:nymax-1,2)
+      integer icloudbot(0:nxmax-1,0:nymax-1,2)
+      integer icloudthck(0:nxmax-1,0:nymax-1,2)
 
 
@@ -358,4 +359,7 @@
   !      rainout  conv/lsp dominated  2/3
   !      washout  conv/lsp dominated  4/5
+! PS 2013
+!c icloudbot (m)        cloud bottom height
+!c icloudthck (m)       cloud thickness     
 
   ! pplev for the GFS version
@@ -684,11 +688,9 @@
 
   !********************
-  ! Verbosity, testing flags, namelist I/O
+  ! Verbosity, testing flags
   !********************   
   integer :: verbosity=0
   integer :: info_flag=0
-  integer :: count_clock, count_clock0,  count_rate, count_max
-  real    :: tins
-  logical :: nmlout=.true.
+  INTEGER :: count_clock, count_clock0,  count_rate, count_max
    
 
Index: /trunk/src/gridcheck.f90
===================================================================
--- /trunk/src/gridcheck.f90	(revision 30)
+++ /trunk/src/gridcheck.f90	(revision 20)
@@ -84,8 +84,8 @@
   real(kind=4) :: xaux1,xaux2,yaux1,yaux2
   real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
-  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parID
+  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
   !HSO  end
   integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
-  real :: sizesouth,sizenorth,xauxa,pint,conversion_factor
+  real :: sizesouth,sizenorth,xauxa,pint
 
   ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
@@ -98,5 +98,6 @@
 
   integer :: isec1(56),isec2(22+nxmax+nymax)
-  real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
+  !real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
+  real(kind=4) :: zsec2(184),zsec4(jpunp)
   character(len=1) :: opt
 
@@ -171,6 +172,4 @@
   call grib_check(iret,gribFunction,gribErrorMsg)
   call grib_get_int(igrib,'level',valSurf,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'paramId',parId,iret)
   call grib_check(iret,gribFunction,gribErrorMsg)
 
@@ -194,6 +193,4 @@
   elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
     isec1(6)=135         ! indicatorOfParameter
-  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
-    isec1(6)=135         ! indicatorOfParameter
   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
     isec1(6)=151         ! indicatorOfParameter
@@ -208,7 +205,7 @@
   elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
     isec1(6)=141         ! indicatorOfParameter
-  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
+  elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
     isec1(6)=164         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
+  elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
     isec1(6)=142         ! indicatorOfParameter
   elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
@@ -218,11 +215,11 @@
   elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
     isec1(6)=176         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
+  elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
     isec1(6)=180         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
+  elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
     isec1(6)=181         ! indicatorOfParameter
   elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
     isec1(6)=129         ! indicatorOfParameter
-  elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
+  elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
     isec1(6)=160         ! indicatorOfParameter
   elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
@@ -233,8 +230,4 @@
          parCat,parNum,typSurf
   endif
-  if(parId .ne. isec1(6) .and. parId .ne. 77) then
-    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
-!    stop
-  endif
 
   endif
@@ -272,5 +265,5 @@
 
   !HSO  get the second part of the grid dimensions only from GRiB1 messages
-  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
+  if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
     call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
          xaux2in,iret)
@@ -298,4 +291,5 @@
     dyconst=180./(dy*r_earth*pi)
     gotGrid=1
+
   ! Check whether fields are global
   ! If they contain the poles, specify polar stereographic map
@@ -443,11 +437,13 @@
   !********************
 
-  write(*,'(a,2i7)') ' Vertical levels in ECMWF data: ', &
+  write(*,*)
+  write(*,*)
+  write(*,'(a,2i7)') '# of vertical levels in ECMWF data: ', &
        nuvz+1,nwz
   write(*,*)
-  write(*,'(a)') ' Mother domain:'
-  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
+  write(*,'(a)') 'Mother domain:'
+  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
        xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
-  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range : ', &
+  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
        ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
   write(*,*)
@@ -471,5 +467,4 @@
     k=nlev_ec+1+numskip+i
     akm(nwz-i+1)=zsec2(j)
-  !   write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)
     bkm(nwz-i+1)=zsec2(k)
   end do
@@ -558,3 +553,2 @@
 
 end subroutine gridcheck
-
Index: /trunk/src/gridcheck_gfs.f90
===================================================================
--- /trunk/src/gridcheck_gfs.f90	(revision 30)
+++ /trunk/src/gridcheck_gfs.f90	(revision 20)
@@ -423,5 +423,5 @@
   write(*,*)
   write(*,*)
-  write(*,'(a,2i7)') 'Vertical levels in NCEP data: ', &
+  write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', &
        nuvz,nwz
   write(*,*)
@@ -429,5 +429,5 @@
   write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
        xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
-  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range : ', &
+  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
        ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
   write(*,*)
Index: /trunk/src/gridcheck_nests.f90
===================================================================
--- /trunk/src/gridcheck_nests.f90	(revision 30)
+++ /trunk/src/gridcheck_nests.f90	(revision 20)
@@ -48,5 +48,4 @@
   integer :: igrib
   integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
-  integer :: parID !added by mc for making it consistent with new gridcheck.f90
   integer :: gotGrib
   !HSO  end
@@ -56,5 +55,4 @@
   real(kind=4) :: xaux1,xaux2,yaux1,yaux2
   real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
-  real :: conversion_factor !added by mc to make it consistent with new gridchek.f90
 
   ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
@@ -152,6 +150,4 @@
   call grib_get_int(igrib,'level',valSurf,iret)
   call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consisitent with new grid_check.f90
-  call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consisitent with new  grid_check.f90
 
   !print*,discipl,parCat,parNum,typSurf,valSurf
@@ -174,6 +170,4 @@
   elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
     isec1(6)=135         ! indicatorOfParameter
-  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added bymc to make it consistent with new gridchek.f90
-    isec1(6)=135         ! indicatorOfParameter    !                     ! " " " " " " " " " " " " " " " " " " " " " " " "  " " " 
   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
     isec1(6)=151         ! indicatorOfParameter
@@ -188,7 +182,7 @@
   elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
     isec1(6)=141         ! indicatorOfParameter
-  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
+  elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
     isec1(6)=164         ! indicatorOfParameter
- 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
+  elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
     isec1(6)=142         ! indicatorOfParameter
   elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
@@ -198,11 +192,11 @@
   elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
     isec1(6)=176         ! indicatorOfParameter
-  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
+  elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
     isec1(6)=180         ! indicatorOfParameter
-  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
+  elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
     isec1(6)=181         ! indicatorOfParameter
   elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
     isec1(6)=129         ! indicatorOfParameter
-  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
+  elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
     isec1(6)=160         ! indicatorOfParameter
   elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
@@ -213,8 +207,4 @@
          parCat,parNum,typSurf
   endif
-  if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consistent with new gridchek.f90
-    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
-!    stop
-  endif
 
   endif
@@ -247,6 +237,6 @@
 
   !HSO  get the second part of the grid dimensions only from GRiB1 messages
-  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!!
-    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
+  if ((gribVer.eq.1).and.(gotGrib.eq.0)) then
+    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
          xaux1in,iret)
     call grib_check(iret,gribFunction,gribErrorMsg)
@@ -264,14 +254,14 @@
     yaux1=yaux1in
     yaux2=yaux2in
-    if(xaux1.gt.180.) xaux1=xaux1-360.0
-    if(xaux2.gt.180.) xaux2=xaux2-360.0
-    if(xaux1.lt.-180.) xaux1=xaux1+360.0
-    if(xaux2.lt.-180.) xaux2=xaux2+360.0
-    if (xaux2.lt.xaux1) xaux2=xaux2+360.0
+    if(xaux1.gt.180) xaux1=xaux1-360.0
+    if(xaux2.gt.180) xaux2=xaux2-360.0
+    if(xaux1.lt.-180) xaux1=xaux1+360.0
+    if(xaux2.lt.-180) xaux2=xaux2+360.0
+    if (xaux2.lt.xaux1) xaux2=xaux2+360.
     xlon0n(l)=xaux1
     ylat0n(l)=yaux1
     dxn(l)=(xaux2-xaux1)/real(nxn(l)-1)
     dyn(l)=(yaux2-yaux1)/real(nyn(l)-1)
-    gotGrib=1 !commetn by mc note tahthere gotGRIB is used instead of gotGrid!!!
+    gotGrib=1
   endif ! ifield.eq.1
 
@@ -350,9 +340,9 @@
   !********************
 
-  write(*,'(a,i2,a)') ' Nested domain ',l,':'
-  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
+  write(*,'(a,i2)') 'Nested domain #: ',l
+  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
        xlon0n(l),' to ',xlon0n(l)+(nxn(l)-1)*dxn(l), &
        '   Grid distance: ',dxn(l)
-  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range : ', &
+  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
        ylat0n(l),' to ',ylat0n(l)+(nyn(l)-1)*dyn(l), &
        '   Grid distance: ',dyn(l)
Index: /trunk/src/interpol_rain.f90
===================================================================
--- /trunk/src/interpol_rain.f90	(revision 30)
+++ /trunk/src/interpol_rain.f90	(revision 20)
@@ -20,8 +20,15 @@
 !**********************************************************************
 
-subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
-       ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
-  !                          i   i   i    i    i     i   i
-  !i    i    i  i    i     i      i      i     o     o     o
+!subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
+!       ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
+!  !                          i   i   i    i    i     i   i
+!  !i    i    i  i    i     i      i      i     o     o     o
+
+      subroutine interpol_rain(yy1,yy2,yy3,iy1,iy2,nxmax,nymax,nzmax,nx, &
+     ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, &
+     intiy1,intiy2,icmv)
+!                               i   i   i    i    i     i   i 
+!     i    i    i  i    i     i      i      i     o     o     o
+
   !****************************************************************************
   !                                                                           *
@@ -40,5 +47,8 @@
   !                                                                           *
   !     30 August 1996                                                        *
-  !                                                                           *
+  !
+  !*  Petra Seibert, 2011/2012:
+  !*  Add interpolation of cloud bottom and cloud thickness
+  !*  for fix to SE's new wet scavenging scheme  *
   !****************************************************************************
   !                                                                           *
@@ -70,10 +80,15 @@
 
   integer :: nx,ny,nxmax,nymax,nzmax,memind(2),m,ix,jy,ixp,jyp
-  integer :: itime,itime1,itime2,level,indexh
+  !integer :: itime,itime1,itime2,level,indexh
+  integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4 
+  integer :: intiy1,intiy2,ipsum,icmv  
   real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2)
   real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2)
   real :: yy3(0:nxmax-1,0:nymax-1,nzmax,2)
-  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
-  real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
+  integer iy1(0:nxmax-1,0:nymax-1,2),iy2(0:nxmax-1,0:nymax-1,2)
+  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),yi1(2),yi2(2)
+  !real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
+  real :: xt,yt,yint1,yint2,yint3,yint4,p1,p2,p3,p4  
+  !real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
 
 
@@ -127,4 +142,43 @@
          + p3*yy3(ix ,jyp,level,indexh) &
          + p4*yy3(ixp,jyp,level,indexh)
+
+!CPS clouds:
+        ip1=1
+        ip2=1
+        ip3=1
+        ip4=1
+        if (iy1(ix ,jy ,indexh) .eq. icmv) ip1=0
+        if (iy1(ixp,jy ,indexh) .eq. icmv) ip2=0
+        if (iy1(ix ,jyp,indexh) .eq. icmv) ip3=0
+        if (iy1(ixp,jyp,indexh) .eq. icmv) ip4=0
+        ipsum= ip1+ip2+ip3+ip4
+        if (ipsum .eq. 0) then
+          yi1(m)=icmv
+        else
+          yi1(m)=(ip1*p1*iy1(ix ,jy ,indexh)    &
+             + ip2*p2*iy1(ixp,jy ,indexh)       &
+             + ip3*p3*iy1(ix ,jyp,indexh)       &
+             + ip4*p4*iy1(ixp,jyp,indexh))/ipsum
+        endif
+        
+        ip1=1
+        ip2=1
+        ip3=1
+        ip4=1
+        if (iy2(ix ,jy ,indexh) .eq. icmv) ip1=0
+        if (iy2(ixp,jy ,indexh) .eq. icmv) ip2=0
+        if (iy2(ix ,jyp,indexh) .eq. icmv) ip3=0
+        if (iy2(ixp,jyp,indexh) .eq. icmv) ip4=0
+        ipsum= ip1+ip2+ip3+ip4
+        if (ipsum .eq. 0) then
+          yi2(m)=icmv
+        else
+          yi2(m)=(ip1*p1*iy2(ix ,jy ,indexh)  &
+             + ip2*p2*iy2(ixp,jy ,indexh)     &
+             + ip3*p3*iy2(ix ,jyp,indexh)     &
+             + ip4*p4*iy2(ixp,jyp,indexh))/ipsum
+        endif
+!CPS end clouds
+
   end do
 
@@ -134,4 +188,11 @@
   !************************************
 
+      if (abs(itime) .lt. abs(itime1)) then
+        print*,'interpol_rain.f90'
+        print*,itime,itime1,itime2
+        stop 'ITIME PROBLEM'
+      endif
+
+
   dt1=real(itime-itime1)
   dt2=real(itime2-itime)
@@ -143,3 +204,20 @@
 
 
+!PS clouds:
+      intiy1=(yi1(1)*dt2 + yi1(2)*dt1)/dt
+      if (yi1(1) .eq. float(icmv)) intiy1=yi1(2) 
+      if (yi1(2) .eq. float(icmv)) intiy1=yi1(1) 
+
+      intiy2=(yi2(1)*dt2 + yi2(2)*dt1)/dt
+      if (yi2(1) .eq. float(icmv)) intiy2=yi2(2)  
+      if (yi2(2) .eq. float(icmv)) intiy2=yi2(1) 
+      
+      if (intiy1 .ne. icmv .and. intiy2 .ne. icmv) then
+        intiy2 = intiy1 + intiy2 ! convert cloud thickness to cloud top
+      else
+        intiy1=icmv
+        intiy2=icmv
+      endif
+!PS end clouds
+
 end subroutine interpol_rain
Index: /trunk/src/makefile
===================================================================
--- /trunk/src/makefile	(revision 30)
+++ /trunk/src/makefile	(revision 20)
@@ -1,11 +1,28 @@
 SHELL = /bin/bash
-MAIN = FP_ecmwf_gfortran
+TARGET = local
+WINDS=ecmwf
+#WINDS=gfs
+#WINDS=fnl
 
 FC       = gfortran
-INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
-LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
-LIBPATH2 =   /usr/lib/x86_64-linux-gnu/
 FFLAGS   =   -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-#FFLAGS   =   -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
+
+ifeq ($(TARGET),dmz)	
+  # options for ganglia
+  INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
+  LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
+  LIBPATH2 = /usr/lib/x86_64-linux-gnu/
+  MAIN = FLEXPART_dmz_
+endif
+ifeq ($(TARGET),local)
+  # local options 	
+  #libs_dir=/.../flexpart/libs/
+  libs_dir=/Users/ignacio/flexpart/libs/
+  INCPATH  = $(libs_dir)/grib_api-1.9.9_dir/include
+  LIBPATH1 = $(libs_dir)/grib_api-1.9.9_dir/lib
+  LIBPATH2 = $(libs_dir)/jasper_dir/lib
+  MAIN = FLEXPART_local_
+endif
+
 LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
 
@@ -20,5 +37,5 @@
 OBJECTS = \
 writeheader.o  writeheader_txt.o   writeheader_surf.o       assignland.o\
-calcpar.o               part0.o \
+               part0.o \
 caldate.o               partdep.o \
 coordtrafo.o            psih.o \
@@ -36,26 +53,24 @@
 getrc.o                 readreleases.o \
 getvdep.o               readspecies.o \
-interpol_misslev.o      readwind.o \
-conccalc.o              richardson.o \
+interpol_misslev.o      \
+conccalc.o              \
 concoutput.o  concoutput_surf.o          scalev.o \
 pbl_profile.o           readOHfield.o\
 juldate.o               timemanager.o \
 interpol_vdep.o         interpol_rain.o \
-verttransform.o         partoutput.o \
+partoutput.o \
 hanna.o                 wetdepokernel.o \
 mean.o                  wetdepo.o \
 hanna_short.o           windalign.o \
-obukhov.o               gridcheck.o \
 hanna1.o                initialize.o \
-                        gridcheck_nests.o \
-readwind_nests.o        calcpar_nests.o \
+                           calcpar_nests.o \
 verttransform_nests.o   interpol_all_nests.o \
 interpol_wind_nests.o   interpol_misslev_nests.o \
 interpol_vdep_nests.o   interpol_rain_nests.o \
-getvdep_nests.o \
+getvdep_nests.o   gridcheck_nests.o \
+readwind_nests.o  \
 readageclasses.o        readpartpositions.o \
 calcfluxes.o            fluxoutput.o \
 qvsat.o                 skplin.o \
-convmix.o               calcmatrix.o \
 convect43c.o               redist.o \
 sort2.o                 distance.o \
@@ -76,6 +91,22 @@
 
 
-$(MAIN): $(MODOBJS) $(OBJECTS)
-	$(FC) *.o -o $(MAIN) $(LDFLAGS)
+ifeq ($(WINDS),ecmwf)
+  OBJECTS_WINDS = \
+  calcpar.o          readwind.o \
+  richardson.o       verttransform.o \
+  obukhov.o          gridcheck.o  \
+  convmix.o          calcmatrix.o 
+endif
+
+ifeq ($(WINDS),gfs)
+  OBJECTS_WINDS = \
+  calcpar_gfs.o          readwind_gfs.o \
+  richardson_gfs.o       verttransform_gfs.o \
+  obukhov_gfs.o          gridcheck_gfs.o  \
+  convmix_gfs.o          calcmatrix_gfs.o
+endif
+
+$(MAIN): $(MODOBJS) $(OBJECTS)  $(OBJECTS_WINDS)
+	$(FC) *.o -o $(MAIN)_$(WINDS) $(LDFLAGS)
 
 $(OBJECTS): $(MODOBJS)
@@ -87,4 +118,2 @@
 	rm *.o *.mod
 
-cleanall:
-	rm *.o *.mod $(MAIN)
Index: unk/src/makefile.gfs
===================================================================
--- /trunk/src/makefile.gfs	(revision 30)
+++ 	(revision )
@@ -1,98 +1,0 @@
-SHELL = /bin/bash
-MAIN = FLEXPART_GFS_GFORTRAN
-#
-#TARGET = other
-
-FC       = gfortran
-INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
-LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
-LIBPATH2 = /flex_wrk/flexpart/lib64/gfortran/lib/
-FFLAGS   =   -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-#
-
-#ifeq ($(TARGET),other)	
-#  libs_dir=/Users/ignacio/flexpart/libs/
-#  INCPATH  = $(libs_dir)/grib_api-1.9.9_dir/include
-#  LIBPATH1 = $(libs_dir)/grib_api-1.9.9_dir/lib
-#  LIBPATH2 = $(libs_dir)/jasper_dir/lib
-#  #MAIN = FLEXPART_laptop
-#endif
-
-MODOBJS = \
-par_mod.o               com_mod.o \
-conv_mod.o              hanna_mod.o \
-interpol_mod.o          cmapf_mod.o \
-unc_mod.o               oh_mod.o \
-xmass_mod.o             flux_mod.o \
-point_mod.o             outg_mod.o
-
-OBJECTS = \
-writeheader.o    writeheader_txt.o   writeheader_surf.o       assignland.o\
-calcpar_gfs.o           part0.o \
-caldate.o               partdep.o \
-coordtrafo.o            psih.o \
-raerod.o \
-drydepokernel.o         random.o \
-erf.o                   readavailable.o \
-ew.o                    readcommand.o \
-advance.o               readdepo.o \
-releaseparticles.o      psim.o     \
-FLEXPART.o              readlanduse.o \
-getfields.o             init_domainfill.o\
-interpol_wind.o         readoutgrid.o \
-interpol_all.o          readpaths.o \
-getrb.o                 readreceptors.o \
-getrc.o                 readreleases.o \
-getvdep.o               readspecies.o \
-interpol_misslev.o      readwind_gfs.o \
-conccalc.o              richardson_gfs.o \
-concoutput.o    concoutput_surf.o         scalev.o \
-pbl_profile.o           readOHfield.o\
-juldate.o               timemanager.o \
-interpol_vdep.o         interpol_rain.o \
-verttransform_gfs.o         partoutput.o \
-hanna.o                 wetdepokernel.o \
-mean.o                  wetdepo.o \
-hanna_short.o           windalign.o \
-obukhov_gfs.o           gridcheck_gfs.o \
-hanna1.o                initialize.o \
-                        gridcheck_nests.o \
-readwind_nests.o        calcpar_nests.o \
-verttransform_nests.o   interpol_all_nests.o \
-interpol_wind_nests.o   interpol_misslev_nests.o \
-interpol_vdep_nests.o   interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o        readpartpositions.o \
-calcfluxes.o            fluxoutput.o \
-qvsat.o                 skplin.o \
-convmix_gfs.o               calcmatrix_gfs.o \
-convect43c.o            redist.o \
-sort2.o                 distance.o \
-centerofmass.o          plumetraj.o \
-openouttraj.o           calcpv.o \
-calcpv_nests.o          distance2.o \
-clustering.o            interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o           outgrid_init.o \
-openreceptors.o         boundcond_domainfill.o\
-partoutput_short.o      readoutgrid_nest.o \
-outgrid_init_nest.o     writeheader_nest.o writeheader_nest_surf.o\
-concoutput_nest.o  concoutput_surf_nest.o     wetdepokernel_nest.o \
-drydepokernel_nest.o    zenithangle.o \
-ohreaction.o            getvdep_nests.o \
-initial_cond_calc.o     initial_cond_output.o \
-dynamic_viscosity.o     get_settling.o
-
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
-	$(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-%.o: %.f90
-	$(FC) -c $(FFLAGS) $<
-
-clean:
-	rm *.o *.mod
-
Index: unk/src/makefile.netcdf
===================================================================
--- /trunk/src/makefile.netcdf	(revision 30)
+++ 	(revision )
@@ -1,95 +1,0 @@
-SHELL = /bin/bash
-MAIN = FP_ecmwf_gfortran_netcdf
-
-FC       = gfortran
-INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
-LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
-LIBPATH2 =   /usr/lib/x86_64-linux-gnu/
-NCLIBPATH = /usr/local/lib
-FFLAGS   =   -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-# netcdf output is included by using the complier switch -DNETCDF_OUTPUT and the C preprocessor -cpp
-FFLAGS   =   -O2 -cpp -DNETCDF_OUTPUT -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-#FFLAGS   =   -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-#LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -L$(NCLIBPATH) -lnetcdf -lnetcdff -lhdf5 -lgrib_api_f90 -lgrib_api -lm -ljasper
-
-MODOBJS = \
-par_mod.o    com_mod.o \
-conv_mod.o              hanna_mod.o \
-interpol_mod.o          cmapf_mod.o \
-unc_mod.o               oh_mod.o \
-xmass_mod.o             flux_mod.o \
-point_mod.o             outg_mod.o \
-netcdf_output_mod.o
-
-OBJECTS = \
-writeheader.o  writeheader_txt.o   writeheader_surf.o       assignland.o\
-calcpar.o               part0.o \
-caldate.o               partdep.o \
-coordtrafo.o            psih.o \
-raerod.o \
-drydepokernel.o         random.o \
-erf.o                   readavailable.o \
-ew.o                    readcommand.o \
-advance.o               readdepo.o \
-releaseparticles.o      psim.o     \
-FLEXPART.o              readlanduse.o \
-getfields.o             init_domainfill.o\
-interpol_wind.o         readoutgrid.o \
-interpol_all.o          readpaths.o \
-getrb.o                 readreceptors.o \
-getrc.o                 readreleases.o \
-getvdep.o               readspecies.o \
-interpol_misslev.o      readwind.o \
-conccalc.o              richardson.o \
-concoutput.o  concoutput_surf.o          scalev.o \
-pbl_profile.o           readOHfield.o\
-juldate.o               timemanager.o \
-interpol_vdep.o         interpol_rain.o \
-verttransform.o         partoutput.o \
-hanna.o                 wetdepokernel.o \
-mean.o                  wetdepo.o \
-hanna_short.o           windalign.o \
-obukhov.o               gridcheck.o \
-hanna1.o                initialize.o \
-                        gridcheck_nests.o \
-readwind_nests.o        calcpar_nests.o \
-verttransform_nests.o   interpol_all_nests.o \
-interpol_wind_nests.o   interpol_misslev_nests.o \
-interpol_vdep_nests.o   interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o        readpartpositions.o \
-calcfluxes.o            fluxoutput.o \
-qvsat.o                 skplin.o \
-convmix.o               calcmatrix.o \
-convect43c.o               redist.o \
-sort2.o                 distance.o \
-centerofmass.o          plumetraj.o \
-openouttraj.o           calcpv.o \
-calcpv_nests.o          distance2.o \
-clustering.o            interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o           outgrid_init.o \
-openreceptors.o         boundcond_domainfill.o\
-partoutput_short.o      readoutgrid_nest.o \
-outgrid_init_nest.o     writeheader_nest.o writeheader_nest_surf.o \
-concoutput_nest.o concoutput_surf_nest.o      wetdepokernel_nest.o \
-drydepokernel_nest.o    zenithangle.o \
-ohreaction.o            getvdep_nests.o \
-initial_cond_calc.o     initial_cond_output.o \
-dynamic_viscosity.o     get_settling.o
-
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
-	$(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-%.o: %.f90
-	$(FC) -c $(FFLAGS) $<
-
-clean:
-	rm *.o *.mod
-
-cleanall:
-	rm *.o *.mod $(MAIN)
Index: /trunk/src/outgrid_init.f90
===================================================================
--- /trunk/src/outgrid_init.f90	(revision 30)
+++ /trunk/src/outgrid_init.f90	(revision 20)
@@ -225,6 +225,6 @@
   !     maxspec,maxpointspec_act,nclassunc,maxageclass
 
-  write (*,*) 'Allocating fields for global output (x,y): ', numxgrid,numygrid
-  write (*,*) 'Allocating fields for nested output (x,y): ', numxgridn,numygridn 
+  write (*,*) ' Allocating fields for nested and global output (x,y): ', &
+       max(numxgrid,numxgridn),max(numygrid,numygridn) 
 
   ! allocate fields for concoutput with maximum dimension of outgrid
Index: /trunk/src/par_mod.f90
===================================================================
--- /trunk/src/par_mod.f90	(revision 30)
+++ /trunk/src/par_mod.f90	(revision 20)
@@ -121,7 +121,6 @@
   !*********************************************
   
-  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !FNL XF
-  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=152,nwzmax=152,nzmax=152 !ECMWF new 
-  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF
+  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !FNL
+  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF
   !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26
   !integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
@@ -155,6 +154,6 @@
 
   !integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
-  integer,parameter :: maxnests=1,nxmaxn=351,nymaxn=351 !ECMWF
-  !integer,parameter :: maxnests=1, nxmaxn=201, nymaxn=161 ! FNL XF
+  !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151 !ECMWF
+  integer,parameter :: maxnests=1, nxmaxn=201, nymaxn=161 ! FNL XF
   ! maxnests                maximum number of nested grids
   ! nxmaxn,nymaxn           maximum dimension of nested wind fields in
@@ -198,6 +197,7 @@
   !**************************************************
 
-  integer,parameter :: maxpart=150000
-  integer,parameter :: maxspec=4
+  !integer,parameter :: maxpart=4000000
+  integer,parameter :: maxpart=2000
+  integer,parameter :: maxspec=6
 
 
Index: unk/src/pathnames
===================================================================
--- /trunk/src/pathnames	(revision 30)
+++ 	(revision )
@@ -1,6 +1,0 @@
-../options/
-../output/
-/
-/flex_wrk/WIND_FIELDS/AVAILABLE_ECMWF_OPER_fields_global
-============================================
-
Index: /trunk/src/readageclasses.f90
===================================================================
--- /trunk/src/readageclasses.f90	(revision 30)
+++ /trunk/src/readageclasses.f90	(revision 20)
@@ -28,7 +28,6 @@
   !                                                                            *
   !     Author: A. Stohl                                                       *
+  !                                                                            *
   !     20 March 2000                                                          *
-  !     HSO, 1 July 2014                                                       *
-  !     Added optional namelist input                                          *
   !                                                                            *
   !*****************************************************************************
@@ -47,13 +46,4 @@
   integer :: i
 
-  ! namelist help variables
-  integer :: readerror
-
-  ! namelist declaration
-  namelist /ageclass/ &
-    nageclass, &
-    lage
-
-  nageclass=-1 ! preset to negative value to identify failed namelist input
 
   ! If age spectra calculation is switched off, set number of age classes
@@ -67,33 +57,17 @@
   endif
 
+
   ! If age spectra claculation is switched on,
   ! open the AGECLASSSES file and read user options
   !************************************************
 
-  open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES',form='formatted',status='old',err=999)
+  open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES', &
+       status='old',err=999)
 
-  ! try to read in as a namelist
-  read(unitageclasses,ageclass,iostat=readerror)
-  close(unitageclasses)
+  do i=1,13
+    read(unitageclasses,*)
+  end do
+  read(unitageclasses,*) nageclass
 
-  if ((nageclass.lt.0).or.(readerror.ne.0)) then
-    open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES',status='old',err=999)
-    do i=1,13
-      read(unitageclasses,*)
-    end do
-    read(unitageclasses,*) nageclass
-    read(unitageclasses,*) lage(1)
-    do i=2,nageclass
-      read(unitageclasses,*) lage(i)
-    end do
-    close(unitageclasses)
-  endif
-
-  ! write ageclasses file in namelist format to output directory if requested
-  if (nmlout.eqv..true.) then
-    open(unitageclasses,file=path(2)(1:length(2))//'AGECLASSES.namelist',err=1000)
-    write(unitageclasses,nml=ageclass)
-    close(unitageclasses)
-  endif
 
   if (nageclass.gt.maxageclass) then
@@ -106,4 +80,5 @@
   endif
 
+  read(unitageclasses,*) lage(1)
   if (lage(1).le.0) then
     write(*,*) ' #### FLEXPART MODEL ERROR! AGE OF FIRST      #### '
@@ -114,4 +89,5 @@
 
   do i=2,nageclass
+    read(unitageclasses,*) lage(i)
     if (lage(i).le.lage(i-1)) then
       write(*,*) ' #### FLEXPART MODEL ERROR! AGE CLASSES     #### '
@@ -129,9 +105,3 @@
   stop
 
-1000  write(*,*) ' #### FLEXPART MODEL ERROR! FILE "AGECLASSES" #### '
-  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(2)(1:length(2))
-  stop
-
-
 end subroutine readageclasses
Index: /trunk/src/readavailable.f90
===================================================================
--- /trunk/src/readavailable.f90	(revision 30)
+++ /trunk/src/readavailable.f90	(revision 20)
@@ -123,6 +123,6 @@
 
   do k=1,numbnests
-  !print*,length(numpath+2*(k-1)+1),length(numpath+2*(k-1)+2),length(4),length(3)
-  !print*,path(numpath+2*(k-1)+2)(1:length(numpath+2*(k-1)+2))
+  print*,length(numpath+2*(k-1)+1),length(numpath+2*(k-1)+2),length(4),length(3)
+  print*,path(numpath+2*(k-1)+2)(1:length(numpath+2*(k-1)+2))
     open(unitavailab,file=path(numpath+2*(k-1)+2) &
          (1:length(numpath+2*(k-1)+2)),status='old',err=998)
Index: /trunk/src/readcommand.f90
===================================================================
--- /trunk/src/readcommand.f90	(revision 30)
+++ /trunk/src/readcommand.f90	(revision 20)
@@ -29,6 +29,4 @@
   !                                                                            *
   !     18 May 1996                                                            *
-  !     HSO, 1 July 2014                                                       *
-  !     Added optional namelist input                                          *
   !                                                                            *
   !*****************************************************************************
@@ -82,36 +80,36 @@
   character(len=50) :: line
   logical :: old
+  logical :: nmlout=.true. !.false.
   integer :: readerror
 
   namelist /command/ &
-  ldirect, &
-  ibdate,ibtime, &
-  iedate,ietime, &
-  loutstep, &
-  loutaver, &
-  loutsample, &
-  itsplit, &
-  lsynctime, &
-  ctl, &
-  ifine, &
-  iout, &
-  ipout, &
-  lsubgrid, &
-  lconvection, &
-  lagespectra, &
-  ipin, &
-  ioutputforeachrelease, &
-  iflux, &
-  mdomainfill, &
-  ind_source, &
-  ind_receptor, &
-  mquasilag, &
-  nested_output, &
-  linit_cond, &
-  lnetcdfout, &
-  surf_only    
+    ldirect, &
+    ibdate,ibtime, &
+    iedate,ietime, &
+    loutstep, &
+    loutaver, &
+    loutsample, &
+    itsplit, &
+    lsynctime, &
+    ctl, &
+    ifine, &
+    iout, &
+    ipout, &
+    lsubgrid, &
+    lconvection, &
+    lagespectra, &
+    ipin, &
+    ioutputforeachrelease, &
+    iflux, &
+    mdomainfill, &
+    ind_source, &
+    ind_receptor, &
+    mquasilag, &
+    nested_output, &
+    linit_cond, &
+    surf_only    
 
   ! Presetting namelist command
-  ldirect=0
+  ldirect=1
   ibdate=20000101
   ibtime=0
@@ -139,5 +137,4 @@
   nested_output=0
   linit_cond=0
-  lnetcdfout=0
   surf_only=0 
 
@@ -145,14 +142,24 @@
   ! Namelist input first: try to read as namelist file
   !**************************************************************************
-  open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old',form='formatted',err=999)
-
-  ! try namelist input (default)
+  open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
+         form='formatted',iostat=readerror)
+  ! If fail, check if file does not exist
+  if (readerror.ne.0) then
+    
+    print*,'***ERROR: file COMMAND not found in ' 
+    print*, path(1)(1:length(1))//'COMMAND'
+    print*, 'Check your pathnames file.'
+    stop
+
+  endif
+
   read(unitcommand,command,iostat=readerror)
   close(unitcommand)
 
-  ! distinguish namelist from fixed text input
-  if ((readerror.ne.0).or.(ldirect.eq.0)) then ! parse as text file format
- 
-    open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', err=999)
+  ! If error in namelist format, try to open with old input code
+  if (readerror.ne.0) then
+
+    open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
+         err=999)
 
     ! Check the format of the COMMAND file (either in free format,
@@ -166,12 +173,9 @@
     if (index(line,'LDIRECT') .eq. 0) then
       old = .false.
-      write(*,*) 'COMMAND in old short format, please update to namelist format'
     else
       old = .true.
-      write(*,*) 'COMMAND in old long format, please update to namelist format'
     endif
     rewind(unitcommand)
 
-
     ! Read parameters
     !****************
@@ -179,4 +183,5 @@
     call skplin(7,unitcommand)
     if (old) call skplin(1,unitcommand)
+
     read(unitcommand,*) ldirect
     if (old) call skplin(3,unitcommand)
@@ -226,6 +231,4 @@
     if (old) call skplin(3,unitcommand)
     read(unitcommand,*) linit_cond
-    if (old) call skplin(3,unitcommand)
-    read(unitcommand,*) surf_only
     close(unitcommand)
 
@@ -234,7 +237,11 @@
   ! write command file in namelist format to output directory if requested
   if (nmlout.eqv..true.) then
+    !open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist.out',status='new',err=1000)
     open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist',err=1000)
     write(unitcommand,nml=command)
     close(unitcommand)
+     ! open(unitheader,file=path(2)(1:length(2))//'header_nml',status='new',err=999)
+     ! write(unitheader,NML=COMMAND)
+     !close(unitheader) 
   endif
 
@@ -346,15 +353,4 @@
   endif
 
-!  check for netcdf output switch (use for non-namelist input only!)
-  if (iout.ge.8) then
-     lnetcdfout = 1
-     iout = iout - 8
-#ifndef NETCDF_OUTPUT
-     print*,'ERROR: netcdf output not activated during compile time but used in COMMAND file!'
-     print*,'Please recompile with netcdf library or use standard output format.'
-     stop
-#endif
-  endif
-
   ! Check whether a valid option for gridded model output has been chosen
   !**********************************************************************
@@ -362,7 +358,5 @@
   if ((iout.lt.1).or.(iout.gt.5)) then
     write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
-    write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4, OR 5 FOR       #### '
-    write(*,*) ' #### STANDARD FLEXPART OUTPUT OR  9 - 13     #### '
-    write(*,*) ' #### FOR NETCDF OUTPUT                       #### '
+    write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4, OR 5!          #### '
     stop
   endif
@@ -376,4 +370,5 @@
     stop
   endif
+
 
 
@@ -597,6 +592,6 @@
 
 1000   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "COMMAND"    #### '
-  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(2)(1:length(1))
-  stop
+       write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+        write(*,'(a)') path(2)(1:length(1))
+        stop
 end subroutine readcommand
Index: /trunk/src/readoutgrid.f90
===================================================================
--- /trunk/src/readoutgrid.f90	(revision 30)
+++ /trunk/src/readoutgrid.f90	(revision 20)
@@ -29,6 +29,4 @@
   !                                                                            *
   !     4 June 1996                                                            *
-  !     HSO, 1 July 2014
-  !     Added optional namelist input
   !                                                                            *
   !*****************************************************************************
@@ -55,58 +53,32 @@
   real,parameter :: eps=1.e-4
 
-  ! namelist variables
-  integer, parameter :: maxoutlev=500
-  integer :: readerror
-  real,allocatable, dimension (:) :: outheights
 
-  ! declare namelist
-  namelist /outgrid/ &
-    outlon0,outlat0, &
-    numxgrid,numygrid, &
-    dxout,dyout, &
-    outheights
-
-  ! allocate large array for reading input
-  allocate(outheights(maxoutlev),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheights'
-
-  ! helps identifying failed namelist input
-  dxout=-1.0
-  outheights=-1.0
 
   ! Open the OUTGRID file and read output grid specifications
   !**********************************************************
 
-  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old',form='formatted',err=999)
+  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old', &
+       err=999)
 
-  ! try namelist input
-  read(unitoutgrid,outgrid,iostat=readerror)
-  close(unitoutgrid)
 
-  if ((dxout.le.0).or.(readerror.ne.0)) then
+  call skplin(5,unitoutgrid)
 
-    readerror=1
 
-    open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old',err=999)
+  ! 1.  Read horizontal grid specifications
+  !****************************************
 
-    call skplin(5,unitoutgrid)
+  call skplin(3,unitoutgrid)
+  read(unitoutgrid,'(4x,f11.4)') outlon0
+  call skplin(3,unitoutgrid)
+  read(unitoutgrid,'(4x,f11.4)') outlat0
+  call skplin(3,unitoutgrid)
+  read(unitoutgrid,'(4x,i5)') numxgrid
+  call skplin(3,unitoutgrid)
+  read(unitoutgrid,'(4x,i5)') numygrid
+  call skplin(3,unitoutgrid)
+  read(unitoutgrid,'(4x,f12.5)') dxout
+  call skplin(3,unitoutgrid)
+  read(unitoutgrid,'(4x,f12.5)') dyout
 
-    ! 1.  Read horizontal grid specifications
-    !****************************************
-
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f11.4)') outlon0
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f11.4)') outlat0
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,i5)') numxgrid
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,i5)') numygrid
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f12.5)') dxout
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f12.5)') dyout
-
-  endif
 
   ! Check validity of output grid (shall be within model domain)
@@ -119,5 +91,7 @@
   if ((outlon0+eps.lt.xlon0).or.(outlat0+eps.lt.ylat0) &
        .or.(xr.gt.xr1+eps).or.(yr.gt.yr1+eps)) then
+    write(*,*) 'outlon0,outlat0:'
     write(*,*) outlon0,outlat0
+    write(*,*) 'xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout:'
     write(*,*) xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout
     write(*,*) ' #### FLEXPART MODEL ERROR! PART OF OUTPUT    ####'
@@ -130,8 +104,6 @@
   ! 2. Count Vertical levels of output grid
   !****************************************
-
-  if (readerror.ne.0) then
-    j=0
-100 j=j+1
+  j=0
+100   j=j+1
     do i=1,3
       read(unitoutgrid,*,end=99)
@@ -140,58 +112,31 @@
     if (outhelp.eq.0.) goto 99
     goto 100
-99  numzgrid=j-1
-  else
-    do i=1,maxoutlev
-      if (outheights(i).lt.0) exit
-    end do
-    numzgrid=i-1
-  end if
+99   numzgrid=j-1
 
-  allocate(outheight(numzgrid),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheight'
-  allocate(outheighthalf(numzgrid),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheighthalf'
+    allocate(outheight(numzgrid) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
+    allocate(outheighthalf(numzgrid) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
+
+
+  rewind(unitoutgrid)
+  call skplin(29,unitoutgrid)
 
   ! 2. Vertical levels of output grid
   !**********************************
 
-  if (readerror.ne.0) then
+  j=0
+1000   j=j+1
+    do i=1,3
+      read(unitoutgrid,*,end=990)
+    end do
+    read(unitoutgrid,'(4x,f7.1)',end=990) outhelp
+    if (outhelp.eq.0.) goto 99
+    outheight(j)=outhelp
+    goto 1000
+990   numzgrid=j-1
 
-    rewind(unitoutgrid)
-    call skplin(29,unitoutgrid)
-
-    do j=1,numzgrid
-      do i=1,3
-        read(unitoutgrid,*)
-      end do
-      read(unitoutgrid,'(4x,f7.1)') outhelp
-      outheight(j)=outhelp
-      outheights(j)=outhelp
-    end do
-    close(unitoutgrid)
-
-  else
-
-    do j=1,numzgrid
-      outheight(j)=outheights(j)
-    end do
-
-  endif
-
-  ! write outgrid file in namelist format to output directory if requested
-  if (nmlout.eqv..true.) then
-    ! reallocate outheights with actually required dimension for namelist writing
-    deallocate(outheights)
-    allocate(outheights(numzgrid),stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outheights'
-
-    do j=1,numzgrid
-      outheights(j)=outheight(j)
-    end do
-
-    open(unitoutgrid,file=path(2)(1:length(2))//'OUTGRID.namelist',err=1000)
-    write(unitoutgrid,nml=outgrid)
-    close(unitoutgrid)
-  endif
 
   ! Check whether vertical levels are specified in ascending order
@@ -215,28 +160,30 @@
   end do
 
+
   xoutshift=xlon0-outlon0
   youtshift=ylat0-outlat0
+  close(unitoutgrid)
 
-  allocate(oroout(0:numxgrid-1,0:numygrid-1),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate oroout'
-  allocate(area(0:numxgrid-1,0:numygrid-1),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate area'
-  allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate volume'
-  allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate areaeast'
-  allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate areanorth'
+    allocate(oroout(0:numxgrid-1,0:numygrid-1) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
+    allocate(area(0:numxgrid-1,0:numygrid-1) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
+    allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
+    allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
+    allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
   return
 
 
-999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
+999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(1)(1:length(1))
-  stop
-
-1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
-  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(2)(1:length(2))
+  write(*,*) ' #### xxx/flexpart/options                    #### '
   stop
 
Index: /trunk/src/readoutgrid_nest.f90
===================================================================
--- /trunk/src/readoutgrid_nest.f90	(revision 30)
+++ /trunk/src/readoutgrid_nest.f90	(revision 20)
@@ -53,61 +53,42 @@
   real,parameter :: eps=1.e-4
 
-  integer :: readerror
 
-  ! declare namelist
-  namelist /outgridn/ &
-    outlon0n,outlat0n, &
-    numxgridn,numygridn, &
-    dxoutn,dyoutn
-
-  ! helps identifying failed namelist input
-  dxoutn=-1.0
 
   ! Open the OUTGRID file and read output grid specifications
   !**********************************************************
 
-  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',form='formatted',status='old',err=999)
+  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST', &
+       status='old',err=999)
 
-  ! try namelist input
-  read(unitoutgrid,outgridn,iostat=readerror)
-  close(unitoutgrid)
 
-  if ((dxoutn.le.0).or.(readerror.ne.0)) then
+  call skplin(5,unitoutgrid)
 
-    open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',status='old',err=999)
-    call skplin(5,unitoutgrid)
 
-    ! 1.  Read horizontal grid specifications
-    !****************************************
+  ! 1.  Read horizontal grid specifications
+  !****************************************
 
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f11.4)') outlon0n
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f11.4)') outlat0n
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,i5)') numxgridn
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,i5)') numygridn
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f12.5)') dxoutn
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f12.5)') dyoutn
+  call skplin(3,unitoutgrid)
+  read(unitoutgrid,'(4x,f11.4)') outlon0n
+  call skplin(3,unitoutgrid)
+  read(unitoutgrid,'(4x,f11.4)') outlat0n
+  call skplin(3,unitoutgrid)
+  read(unitoutgrid,'(4x,i5)') numxgridn
+  call skplin(3,unitoutgrid)
+  read(unitoutgrid,'(4x,i5)') numygridn
+  call skplin(3,unitoutgrid)
+  read(unitoutgrid,'(4x,f12.5)') dxoutn
+  call skplin(3,unitoutgrid)
+  read(unitoutgrid,'(4x,f12.5)') dyoutn
 
-    close(unitoutgrid)
-  endif
 
-  ! write outgrid_nest file in namelist format to output directory if requested
-  if (nmlout.eqv..true.) then
-    open(unitoutgrid,file=path(2)(1:length(2))//'OUTGRID_NEST.namelist',err=1000)
-    write(unitoutgrid,nml=outgridn)
-    close(unitoutgrid)
-  endif
-
-  allocate(orooutn(0:numxgridn-1,0:numygridn-1),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate orooutn'
-  allocate(arean(0:numxgridn-1,0:numygridn-1),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate arean'
-  allocate(volumen(0:numxgridn-1,0:numygridn-1,numzgrid),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate volumen'
+    allocate(orooutn(0:numxgridn-1,0:numygridn-1) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
+    allocate(arean(0:numxgridn-1,0:numygridn-1) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
+    allocate(volumen(0:numxgridn-1,0:numygridn-1,numzgrid) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
 
   ! Check validity of output grid (shall be within model domain)
@@ -129,14 +110,11 @@
   xoutshiftn=xlon0-outlon0n
   youtshiftn=ylat0-outlat0n
+  close(unitoutgrid)
   return
 
-999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
+
+999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE OUTGRID_NEST #### '
   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(1)(1:length(1))
-  stop
-
-1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
-  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(2)(1:length(2))
+  write(*,*) ' #### xxx/flexpart/options                    #### '
   stop
 
Index: /trunk/src/readpaths.f90
===================================================================
--- /trunk/src/readpaths.f90	(revision 30)
+++ /trunk/src/readpaths.f90	(revision 20)
@@ -88,4 +88,5 @@
     length(numpath+2*(i-1)+2)=index(path(numpath+2*(i-1)+2),' ')-1
   end do
+  print*,length(5),length(6)
 
 
Index: /trunk/src/readreceptors.f90
===================================================================
--- /trunk/src/readreceptors.f90	(revision 30)
+++ /trunk/src/readreceptors.f90	(revision 20)
@@ -27,7 +27,6 @@
   !                                                                            *
   !     Author: A. Stohl                                                       *
+  !                                                                            *
   !     1 August 1996                                                          *
-  !     HSO, 14 August 2013
-  !     Added optional namelist input
   !                                                                            *
   !*****************************************************************************
@@ -52,13 +51,4 @@
   character(len=16) :: receptor
 
-  integer :: readerror
-  real :: lon,lat   ! for namelist input, lon/lat are used instead of x,y
-  integer,parameter :: unitreceptorout=2
-
-  ! declare namelist
-  namelist /receptors/ &
-    receptor, lon, lat
-
-  lon=-999.9
 
   ! For backward runs, do not allow receptor output. Thus, set number of receptors to zero
@@ -74,25 +64,15 @@
   !************************************************************
 
-  open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS',form='formatted',status='old',err=999)
+  open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS', &
+       status='old',err=999)
 
-  ! try namelist input
-  read(unitreceptor,receptors,iostat=readerror)
+  call skplin(5,unitreceptor)
 
-  ! prepare namelist output if requested
-  if (nmlout.eqv..true.) then
-    open(unitreceptorout,file=path(2)(1:length(2))//'RECEPTORS.namelist',access='append',status='new',err=1000)
-  endif
 
-  if ((lon.lt.-900).or.(readerror.ne.0)) then
+  ! Read the names and coordinates of the receptors
+  !************************************************
 
-    close(unitreceptor)
-    open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS',status='old',err=999)
-    call skplin(5,unitreceptor)
-
-    ! Read the names and coordinates of the receptors
-    !************************************************
-
-    j=0
-100 j=j+1
+  j=0
+100   j=j+1
     read(unitreceptor,*,end=99)
     read(unitreceptor,*,end=99)
@@ -109,8 +89,8 @@
     endif
     if (j.gt.maxreceptor) then
-      write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
-      write(*,*) ' #### POINTS ARE GIVEN.                       #### '
-      write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
-      write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
+    write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
+    write(*,*) ' #### POINTS ARE GIVEN.                       #### '
+    write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
+    write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
     endif
     receptorname(j)=receptor
@@ -120,67 +100,15 @@
     ym=r_earth*dy/180.*pi
     receptorarea(j)=xm*ym
-
-    ! write receptors file in namelist format to output directory if requested
-    if (nmlout.eqv..true.) then
-      lon=x
-      lat=y
-      write(unitreceptorout,nml=receptors)
-    endif
-
     goto 100
 
-99  numreceptor=j-1
-
-  else ! continue with namelist input
-
-    j=0
-    do while (readerror.eq.0)
-      j=j+1
-      lon=-999.9
-      read(unitreceptor,receptors,iostat=readerror)
-      if ((lon.lt.-900).or.(readerror.ne.0)) then
-        readerror=1
-      else
-        if (j.gt.maxreceptor) then
-          write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
-          write(*,*) ' #### POINTS ARE GIVEN.                       #### '
-          write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
-          write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
-        endif
-        receptorname(j)=receptor
-        xreceptor(j)=(lon-xlon0)/dx       ! transform to grid coordinates
-        yreceptor(j)=(lat-ylat0)/dy
-        xm=r_earth*cos(lat*pi/180.)*dx/180.*pi
-        ym=r_earth*dy/180.*pi
-        receptorarea(j)=xm*ym
-      endif
-
-      ! write receptors file in namelist format to output directory if requested
-      if (nmlout.eqv..true.) then
-        write(unitreceptorout,nml=receptors)
-      endif
-
-    end do
-    numreceptor=j-1
-    close(unitreceptor)
-
-  endif
-
-  if (nmlout.eqv..true.) then
-    close(unitreceptorout)
-  endif
-
+99   numreceptor=j-1
+  close(unitreceptor)
   return
 
 
-999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
+999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
   write(*,'(a)') path(1)(1:length(1))
   stop
 
-1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"    #### '
-  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(2)(1:length(2))
-  stop
-
 end subroutine readreceptors
Index: /trunk/src/readreleases.f90
===================================================================
--- /trunk/src/readreleases.f90	(revision 30)
+++ /trunk/src/readreleases.f90	(revision 20)
@@ -33,6 +33,4 @@
   !     Update: 29 January 2001                                                *
   !     Release altitude can be either in magl or masl                         *
-  !     HSO, 12 August 2013
-  !     Added optional namelist input
   !                                                                            *
   !*****************************************************************************
@@ -57,7 +55,5 @@
   ! xpoint2,ypoint2     geograf. coordinates of upper right corner of release  *
   !                     area                                                   *
-  ! weta, wetb          parameters to determine the below-cloud scavenging     *
-  ! weta_in, wetb_in    parameters to determine the in-cloud scavenging        *
-  ! wetc_in, wetd_in    parameters to determine the in-cloud scavenging        *
+  ! weta, wetb          parameters to determine the wet scavenging coefficient *
   ! zpoint1,zpoint2     height range, over which release takes place           *
   ! num_min_discrete    if less, release cannot be randomized and happens at   *
@@ -73,5 +69,5 @@
   implicit none
 
-  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat
+  integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
   integer,parameter :: num_min_discrete=100
   real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
@@ -80,187 +76,119 @@
   logical :: old
 
-  ! help variables for namelist reading
-  integer :: numpoints, parts, readerror
-  integer*2 :: zkind
-  integer :: idate1, itime1, idate2, itime2
-  real :: lon1,lon2,lat1,lat2,z1,z2
-  character*40 :: comment
-  integer,parameter :: unitreleasesout=2
-  real,allocatable, dimension (:) :: mass
-  integer,allocatable, dimension (:) :: specnum_rel,specnum_rel2
-
-  ! declare namelists
-  namelist /releases_ctrl/ &
-    nspec, &
-    specnum_rel
-
-  namelist /release/ &
-    idate1, itime1, &
-    idate2, itime2, &
-    lon1, lon2, &
-    lat1, lat2, &
-    z1, z2, &
-    zkind, &
-    mass, &
-    parts, &
-    comment
+  !sec, read release to find how many releasepoints should be allocated
+
+  open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old', &
+       err=999)
+
+  ! Check the format of the RELEASES file (either in free format,
+  ! or using a formatted mask)
+  ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
+  !**************************************************************************
+
+  call skplin(12,unitreleases)
+  read (unitreleases,901) line
+901   format (a)
+  if (index(line,'Total') .eq. 0) then
+    old = .false.
+  else
+    old = .true.
+  endif
+  rewind(unitreleases)
+
+
+  ! Skip first 11 lines (file header)
+  !**********************************
+
+  call skplin(11,unitreleases)
+
+
+  read(unitreleases,*,err=998) nspec
+  if (old) call skplin(2,unitreleases)
+  do i=1,nspec
+    read(unitreleases,*,err=998) specnum_rel
+    if (old) call skplin(2,unitreleases)
+  end do
 
   numpoint=0
-
-  ! allocate with maxspec for first input loop
-  allocate(mass(maxspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
-  allocate(specnum_rel(maxspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
-
-  ! presetting namelist releases_ctrl
-  nspec = -1  ! use negative value to determine failed namelist input
-  specnum_rel = 0
-
-  !sec, read release to find how many releasepoints should be allocated
-  open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',form='formatted',err=999)
-
-  ! check if namelist input provided
-  read(unitreleases,releases_ctrl,iostat=readerror)
-
-  ! prepare namelist output if requested
-  if (nmlout.eqv..true.) then
-    open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',access='append',status='new',err=1000)
-  endif
-
-  if ((readerror.ne.0).or.(nspec.lt.0)) then
-
-    ! no namelist format, close file and allow reopening in old format
-    close(unitreleases)
-    open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',err=999)
-
-    readerror=1 ! indicates old format
-
-    ! Check the format of the RELEASES file (either in free format,
-    ! or using a formatted mask)
-    ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
-    !**************************************************************************
-
-    call skplin(12,unitreleases)
-    read (unitreleases,901) line
-901 format (a)
-    if (index(line,'Total') .eq. 0) then
-      old = .false.
-    else
-      old = .true.
-    endif
-    rewind(unitreleases)
-
-
-    ! Skip first 11 lines (file header)
-    !**********************************
-
-    call skplin(11,unitreleases)
-
-    read(unitreleases,*,err=998) nspec
-    if (old) call skplin(2,unitreleases)
-    do i=1,nspec
-      read(unitreleases,*,err=998) specnum_rel(i)
-      if (old) call skplin(2,unitreleases)
-    end do
-
-    numpoint=0
-100 numpoint=numpoint+1
-    read(unitreleases,*,end=25)
-    read(unitreleases,*,err=998,end=25) idum,idum
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) idum,idum
-    if (old) call skplin(2,unitreleases)
+100   numpoint=numpoint+1
+  read(unitreleases,*,end=25)
+  read(unitreleases,*,err=998,end=25) idum,idum
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) idum,idum
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) xdum
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) xdum
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) xdum
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) xdum
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) idum
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) xdum
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) xdum
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) idum
+  if (old) call skplin(2,unitreleases)
+  do i=1,nspec
     read(unitreleases,*,err=998) xdum
     if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) xdum
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) xdum
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) xdum
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) idum
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) xdum
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) xdum
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) idum
-    if (old) call skplin(2,unitreleases)
-    do i=1,nspec
-      read(unitreleases,*,err=998) xdum
-      if (old) call skplin(2,unitreleases)
-    end do
-    !save compoint only for the first 1000 release points
-    read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
-    if (old) call skplin(1,unitreleases)
-
-    goto 100
-
-25  numpoint=numpoint-1
-
-  else 
-
-    readerror=0
-    do while (readerror.eq.0) 
-      idate1=-1
-      read(unitreleases,release,iostat=readerror)
-      if ((idate1.lt.0).or.(readerror.ne.0)) then
-        readerror=1
-      else
-        numpoint=numpoint+1
-      endif
-    end do
-    readerror=0
-  endif ! if namelist input
-
-  rewind(unitreleases)
-
-  ! allocate arrays of matching size for number of species (namelist output)
-  deallocate(mass)
-  allocate(mass(nspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
-  allocate(specnum_rel2(nspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel2'
-  specnum_rel2=specnum_rel(1:nspec)
-  deallocate(specnum_rel)
-  allocate(specnum_rel(nspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
-  specnum_rel=specnum_rel2
-  deallocate(specnum_rel2)
-
-  !allocate memory for numpoint releaspoints
-  allocate(ireleasestart(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleasestart'
-  allocate(ireleaseend(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleaseend'
-  allocate(xpoint1(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint1'
-  allocate(xpoint2(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint2'
-  allocate(ypoint1(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint1'
-  allocate(ypoint2(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint2'
-  allocate(zpoint1(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint1'
-  allocate(zpoint2(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint2'
-  allocate(kindz(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate kindz'
-  allocate(xmass(numpoint,maxspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmass'
-  allocate(rho_rel(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate rho_rel'
-  allocate(npart(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate npart'
-  allocate(xmasssave(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmasssave'
-
-  write (*,*) 'Releasepoints : ', numpoint
-
-  do i=1,numpoint
-    xmasssave(i)=0.
   end do
+  !save compoint only for the first 1000 release points
+  read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
+  if (old) call skplin(1,unitreleases)
+
+  goto 100
+
+25   numpoint=numpoint-1
+
+  !allocate memory for numpoint releaspoint
+    allocate(ireleasestart(numpoint) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
+    allocate(ireleaseend(numpoint) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
+    allocate(xpoint1(numpoint) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
+    allocate(xpoint2(numpoint) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
+    allocate(ypoint1(numpoint) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
+    allocate(ypoint2(numpoint) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
+    allocate(zpoint1(numpoint) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
+    allocate(zpoint2(numpoint) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
+    allocate(kindz(numpoint) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
+    allocate(xmass(numpoint,maxspec) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
+    allocate(rho_rel(numpoint) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
+    allocate(npart(numpoint) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
+    allocate(xmasssave(numpoint) &
+         ,stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
+
+   write (*,*) ' Releasepoints allocated: ', numpoint
+
+   do i=1,numpoint
+     xmasssave(i)=0.
+   end do
 
   !now save the information
@@ -273,47 +201,40 @@
   end do
 
-  if (readerror.ne.0) then
-    ! Skip first 11 lines (file header)
-    !**********************************
-
-    call skplin(11,unitreleases)
-
-    ! Assign species-specific parameters needed for physical processes
-    !*************************************************************************
-
-    read(unitreleases,*,err=998) nspec
-    if (nspec.gt.maxspec) goto 994
+  rewind(unitreleases)
+
+
+  ! Skip first 11 lines (file header)
+  !**********************************
+
+  call skplin(11,unitreleases)
+
+
+  ! Assign species-specific parameters needed for physical processes
+  !*************************************************************************
+
+  read(unitreleases,*,err=998) nspec
+  if (nspec.gt.maxspec) goto 994
+  if (old) call skplin(2,unitreleases)
+  do i=1,nspec
+    read(unitreleases,*,err=998) specnum_rel
     if (old) call skplin(2,unitreleases)
-  endif
-
-  ! namelist output
-  if (nmlout.eqv..true.) then
-    write(unitreleasesout,nml=releases_ctrl)
-  endif
-
-  do i=1,nspec
-    if (readerror.ne.0) then
-      read(unitreleases,*,err=998) specnum_rel(i)
-      if (old) call skplin(2,unitreleases)
-      call readspecies(specnum_rel(i),i)
-    else
-      call readspecies(specnum_rel(i),i)
-    endif
-
-    ! For backward runs, only 1 species is allowed
-    !*********************************************
-
-    !if ((ldirect.lt.0).and.(nspec.gt.1)) then
-    !write(*,*) '#####################################################'
-    !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
-    !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
-    !write(*,*) '#####################################################'
-    !  stop
-    !endif
-
-    ! Molecular weight
-    !*****************
-
-    if (((iout.eq.2).or.(iout.eq.3)).and.(weightmolar(i).lt.0.)) then
+    call readspecies(specnum_rel,i)
+
+  ! For backward runs, only 1 species is allowed
+  !*********************************************
+
+  !if ((ldirect.lt.0).and.(nspec.gt.1)) then
+  !write(*,*) '#####################################################'
+  !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
+  !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
+  !write(*,*) '#####################################################'
+  !  stop
+  !endif
+
+  ! Molecular weight
+  !*****************
+
+    if (((iout.eq.2).or.(iout.eq.3)).and. &
+         (weightmolar(i).lt.0.)) then
       write(*,*) 'For mixing ratio output, valid molar weight'
       write(*,*) 'must be specified for all simulated species.'
@@ -323,23 +244,25 @@
     endif
 
-    ! Radioactive decay
-    !******************
+
+  ! Radioactive decay
+  !******************
 
     decay(i)=0.693147/decay(i) !conversion half life to decay constant
 
 
-    ! Dry deposition of gases
-    !************************
-
-    if (reldiff(i).gt.0.) rm(i)=1./(henry(i)/3000.+100.*f0(i))    ! mesophyll resistance
-
-    ! Dry deposition of particles
-    !****************************
+  ! Dry deposition of gases
+  !************************
+
+    if (reldiff(i).gt.0.) &
+         rm(i)=1./(henry(i)/3000.+100.*f0(i))    ! mesophyll resistance
+
+  ! Dry deposition of particles
+  !****************************
 
     vsetaver(i)=0.
     cunningham(i)=0.
     dquer(i)=dquer(i)*1000000.         ! Conversion m to um
-    if (density(i).gt.0.) then         ! Additional parameters
-      call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
+    if (density(i).gt.0.) then                  ! Additional parameters
+     call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
       do j=1,ni
         fract(i,j)=fracth(j)
@@ -349,37 +272,26 @@
         vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
       end do
-      write(*,*) 'Average settling velocity: ',i,vsetaver(i)
-    endif
-
-    ! Dry deposition for constant deposition velocity
-    !************************************************
+      write(*,*) 'Average setting velocity: ',i,vsetaver(i)
+    endif
+
+  ! Dry deposition for constant deposition velocity
+  !************************************************
 
     dryvel(i)=dryvel(i)*0.01         ! conversion to m/s
 
-    ! Check if wet deposition or OH reaction shall be calculated
-    !***********************************************************
+  ! Check if wet deposition or OH reaction shall be calculated
+  !***********************************************************
     if (weta(i).gt.0.)  then
       WETDEP=.true.
-      write (*,*) 'Below-cloud scavenging: ON'
-      write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
-    else
-      write (*,*) 'Below-cloud scavenging: OFF'
-    endif
-    
-    ! NIK 31.01.2013 + 10.12.2013
-    if (weta_in(i).gt.0.)  then
-      WETDEP=.true.
-      write (*,*) 'In-cloud scavenging: ON'
-      write (*,*) 'In-cloud scavenging coefficients: ',weta_in(i),wetb_in(i), wetc_in(i), wetd_in(i),i
-    else
-      write (*,*) 'In-cloud scavenging: OFF' 
-    endif
-
+      write (*,*) 'Wetdeposition switched on: ',weta(i),i
+    endif
     if (ohreact(i).gt.0) then
       OHREA=.true.
-      write (*,*) 'OHreaction: ON (',ohreact(i),i,')'
-    endif
-
-    if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or.(dryvel(i).gt.0.)) then
+      write (*,*) 'OHreaction switched on: ',ohreact(i),i
+    endif
+
+
+    if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or. &
+         (dryvel(i).gt.0.)) then
       DRYDEP=.true.
       DRYDEPSPEC(i)=.true.
@@ -388,109 +300,54 @@
   end do
 
-  if (WETDEP.or.DRYDEP) DEP=.true.
+    if (WETDEP.or.DRYDEP) DEP=.true.
 
   ! Read specifications for each release point
   !*******************************************
-  numpoints=numpoint
+
   numpoint=0
   numpartmax=0
   releaserate=0.
-101   numpoint=numpoint+1
-
-  if (readerror.lt.1) then ! reading namelist format
-
-    if (numpoint.gt.numpoints) goto 250
-    zkind = 1
-    mass = 0
-    parts = 0
-    comment = ' '
-    read(unitreleases,release,iostat=readerror)
-    id1=idate1
-    it1=itime1
-    id2=idate2
-    it2=itime2
-    xpoint1(numpoint)=lon1
-    xpoint2(numpoint)=lon2
-    ypoint1(numpoint)=lat1
-    ypoint2(numpoint)=lat2
-    zpoint1(numpoint)=z1
-    zpoint2(numpoint)=z2
-    kindz(numpoint)=zkind
-    do i=1,nspec
-      xmass(numpoint,i)=mass(i)
-    end do
-    npart(numpoint)=parts
-    compoint(min(1001,numpoint))=comment
-
-    ! namelist output
-    if (nmlout.eqv..true.) then
-      write(unitreleasesout,nml=release)
-    endif
-
+1000   numpoint=numpoint+1
+  read(unitreleases,*,end=250)
+  read(unitreleases,*,err=998,end=250) id1,it1
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) id2,it2
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) xpoint1(numpoint)
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) ypoint1(numpoint)
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) xpoint2(numpoint)
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) ypoint2(numpoint)
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) kindz(numpoint)
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) zpoint1(numpoint)
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) zpoint2(numpoint)
+  if (old) call skplin(2,unitreleases)
+  read(unitreleases,*,err=998) npart(numpoint)
+  if (old) call skplin(2,unitreleases)
+  do i=1,nspec
+    read(unitreleases,*,err=998) xmass(numpoint,i)
+    if (old) call skplin(2,unitreleases)
+  end do
+  !save compoint only for the first 1000 release points
+  if (numpoint.le.1000) then
+    read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
   else
-
-    read(unitreleases,*,end=250)
-    read(unitreleases,*,err=998,end=250) id1,it1
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) id2,it2
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) xpoint1(numpoint)
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) ypoint1(numpoint)
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) xpoint2(numpoint)
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) ypoint2(numpoint)
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) kindz(numpoint)
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) zpoint1(numpoint)
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) zpoint2(numpoint)
-    if (old) call skplin(2,unitreleases)
-    read(unitreleases,*,err=998) npart(numpoint)
-    if (old) call skplin(2,unitreleases)
-    do i=1,nspec
-      read(unitreleases,*,err=998) xmass(numpoint,i)
-      if (old) call skplin(2,unitreleases)
-      mass(i)=xmass(numpoint,i)
-    end do
-    !save compoint only for the first 1000 release points
-    if (numpoint.le.1000) then
-      read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
-      comment=compoint(numpoint)(1:40)
-    else
-      read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
-      comment=compoint(1001)(1:40)
-    endif
-    if (old) call skplin(1,unitreleases)
-
-    ! namelist output
-    if (nmlout.eqv..true.) then
-      idate1=id1
-      itime1=it1
-      idate2=id2
-      itime2=it2
-      lon1=xpoint1(numpoint)
-      lon2=xpoint2(numpoint)
-      lat1=ypoint1(numpoint)
-      lat2=ypoint2(numpoint)
-      z1=zpoint1(numpoint)
-      z2=zpoint2(numpoint)
-      zkind=kindz(numpoint)
-      parts=npart(numpoint)
-      write(unitreleasesout,nml=release)
-    endif
-
-    if (numpoint.le.1000) then
-      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
-           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
-           (compoint(numpoint)(1:8).eq.'        ')) goto 250
-    else
-      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
-           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
-    endif
-
-  endif ! if namelist format
+    read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
+  endif
+  if (old) call skplin(1,unitreleases)
+  if (numpoint.le.1000) then
+    if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
+         (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
+         (compoint(numpoint)(1:8).eq.'        ')) goto 250
+  else
+    if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
+         (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
+  endif
+
 
   ! If a release point contains no particles, stop and issue error message
@@ -563,4 +420,9 @@
   endif
 
+
+  ! Check, whether the total number of particles may exceed totally allowed
+  ! number of particles at some time during the simulation
+  !************************************************************************
+
   ! Determine the release rate (particles per second) and total number
   ! of particles released during the simulation
@@ -574,14 +436,10 @@
   endif
   numpartmax=numpartmax+npart(numpoint)
-  goto 101
+  goto 1000
+
 
 250   close(unitreleases)
 
-  if (nmlout.eqv..true.) then
-    close(unitreleasesout)
-  endif
-
-  write (*,*) 'Particles allocated (maxpart)  : ',maxpart
-  write (*,*) 'Particles released (numpartmax): ',numpartmax
+  write (*,*) ' Particles allocated for this run: ',maxpart, ', released in simulation: ',  numpartmax
   numpoint=numpoint-1
 
@@ -591,8 +449,4 @@
     maxpointspec_act=1
   endif
-
-  ! Check, whether the total number of particles may exceed totally allowed
-  ! number of particles at some time during the simulation
-  !************************************************************************
 
   if (releaserate.gt. &
@@ -652,8 +506,3 @@
   stop
 
-1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RELEASES"    #### '
-  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(2)(1:length(2))
-  stop
-
 end subroutine readreleases
Index: /trunk/src/readspecies.f90
===================================================================
--- /trunk/src/readspecies.f90	(revision 30)
+++ /trunk/src/readspecies.f90	(revision 20)
@@ -30,10 +30,4 @@
   !                                                                            *
   !   11 July 1996                                                             *
-  !
-  !   Changes:                                                                 *
-  !   N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging *
-  !                                                                            *
-  !   HSO, 13 August 2013
-  !   added optional namelist input
   !                                                                            *
   !*****************************************************************************
@@ -42,9 +36,5 @@
   ! decaytime(maxtable)  half time for radiological decay                      *
   ! specname(maxtable)   names of chemical species, radionuclides              *
-  ! weta, wetb           Parameters for determining below-cloud scavenging     *
-  ! weta_in              Parameter for determining in-cloud scavenging         *
-  ! wetb_in              Parameter for determining in-cloud scavenging         *
-  ! wetc_in              Parameter for determining in-cloud scavenging         *
-  ! wetd_in              Parameter for determining in-cloud scavenging         *
+  ! wetscava, wetscavb   Parameters for determining scavenging coefficient     *
   ! ohreact              OH reaction rate                                      *
   ! id_spec              SPECIES number as referenced in RELEASE file          *
@@ -65,59 +55,18 @@
   logical :: spec_found
 
-  character(len=16) :: pspecies
-  real :: pdecay, pweta, pwetb, preldiff, phenry, pf0, pdensity, pdquer
-  real :: pdsigma, pdryvel, pweightmolar, pohreact, pspec_ass, pkao
-  real :: pweta_in, pwetb_in, pwetc_in, pwetd_in
-  integer :: readerror
-
-  ! declare namelist
-  namelist /species_params/ &
-   pspecies, pdecay, pweta, pwetb, &
-   pweta_in, pwetb_in, pwetc_in, pwetd_in, &
-   preldiff, phenry, pf0, pdensity, pdquer, &
-   pdsigma, pdryvel, pweightmolar, pohreact, pspec_ass, pkao
-
-  pspecies=" "
-  pdecay=-999.9
-  pweta=-9.9E-09
-  pwetb=0.0
-  pweta_in=-9.9E-09
-  pwetb_in=-9.9E-09
-  pwetc_in=-9.9E-09
-  pwetd_in=-9.9E-09
-  preldiff=-9.9
-  phenry=0.0
-  pf0=0.0
-  pdensity=-9.9E09
-  pdquer=0.0
-  pdsigma=0.0
-  pdryvel=-9.99
-  pohreact=-9.9E-09
-  pspec_ass=-9
-  pkao=-99.99
-  pweightmolar=-789.0 ! read failure indicator value
-
   ! Open the SPECIES file and read species names and properties
   !************************************************************
   specnum(pos_spec)=id_spec
   write(aspecnumb,'(i3.3)') specnum(pos_spec)
-  open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998)
+  open(unitspecies,file= &
+       path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old', &
+       err=998)
   !write(*,*) 'reading SPECIES',specnum(pos_spec)
 
   ASSSPEC=.FALSE.
 
-  ! try namelist input
-  read(unitspecies,species_params,iostat=readerror)
-  close(unitspecies)
-   
-  if ((pweightmolar.eq.-789.0).or.(readerror.ne.0)) then ! no namelist found
-
-    readerror=1
-
-    open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',err=998)
-
-    do i=1,6
-      read(unitspecies,*)
-    end do
+  do i=1,6
+    read(unitspecies,*)
+  end do
 
     read(unitspecies,'(a10)',end=22) species(pos_spec)
@@ -129,15 +78,4 @@
     read(unitspecies,'(f18.2)',end=22) wetb(pos_spec)
   !  write(*,*) wetb(pos_spec)
-
-  !*** NIK 31.01.2013: including in-cloud scavening parameters
-   read(unitspecies,'(e18.1)',end=22) weta_in(pos_spec)
-  !  write(*,*) weta_in(pos_spec)
-   read(unitspecies,'(f18.2)',end=22) wetb_in(pos_spec)
-  !  write(*,*) wetb_in(pos_spec)
-   read(unitspecies,'(f18.2)',end=22) wetc_in(pos_spec)
-  !  write(*,*) wetc_in(pos_spec)
-   read(unitspecies,'(f18.2)',end=22) wetd_in(pos_spec)
-  !  write(*,*) wetd_in(pos_spec)
-
     read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
   !  write(*,*) reldiff(pos_spec)
@@ -162,79 +100,34 @@
     read(unitspecies,'(f18.2)',end=22) kao(pos_spec)
   !       write(*,*) kao(pos_spec)
+    i=pos_spec
 
-    pspecies=species(pos_spec)
-    pdecay=decay(pos_spec)
-    pweta=weta(pos_spec)
-    pwetb=wetb(pos_spec)
-    pweta_in=weta_in(pos_spec)
-    pwetb_in=wetb_in(pos_spec)
-    pwetc_in=wetc_in(pos_spec)
-    pwetd_in=wetd_in(pos_spec)
-    preldiff=reldiff(pos_spec)
-    phenry=henry(pos_spec)
-    pf0=f0(pos_spec)
-    pdensity=density(pos_spec)
-    pdquer=dquer(pos_spec)
-    pdsigma=dsigma(pos_spec)
-    pdryvel=dryvel(pos_spec)
-    pweightmolar=weightmolar(pos_spec)
-    pohreact=ohreact(pos_spec)
-    pspec_ass=spec_ass(pos_spec)
-    pkao=kao(pos_spec)
+    if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
+       if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
+    endif
 
-  else
+    if (spec_ass(pos_spec).gt.0) then
+       spec_found=.FALSE.
+       do j=1,pos_spec-1
+          if (spec_ass(pos_spec).eq.specnum(j)) then
+             spec_ass(pos_spec)=j
+             spec_found=.TRUE.
+             ASSSPEC=.TRUE.
+          endif
+       end do
+       if (spec_found.eqv..FALSE.) then
+          goto 997
+       endif
+    endif
 
-    species(pos_spec)=pspecies
-    decay(pos_spec)=pdecay
-    weta(pos_spec)=pweta
-    wetb(pos_spec)=pwetb
-    weta_in(pos_spec)=pweta_in
-    wetb_in(pos_spec)=pwetb_in
-    wetc_in(pos_spec)=pwetc_in
-    wetd_in(pos_spec)=pwetd_in
-    reldiff(pos_spec)=preldiff
-    henry(pos_spec)=phenry
-    f0(pos_spec)=pf0
-    density(pos_spec)=pdensity
-    dquer(pos_spec)=pdquer
-    dsigma(pos_spec)=pdsigma
-    dryvel(pos_spec)=pdryvel
-    weightmolar(pos_spec)=pweightmolar
-    ohreact(pos_spec)=pohreact
-    spec_ass(pos_spec)=pspec_ass
-    kao(pos_spec)=pkao
+    if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
+    if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
 
-  endif
-
-  i=pos_spec
-
-  if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
-   if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
-  endif
-
-  if (spec_ass(pos_spec).gt.0) then
-    spec_found=.FALSE.
-    do j=1,pos_spec-1
-      if (spec_ass(pos_spec).eq.specnum(j)) then
-        spec_ass(pos_spec)=j
-        spec_found=.TRUE.
-        ASSSPEC=.TRUE.
-      endif
-    end do
-    if (spec_found.eqv..FALSE.) then
-      goto 997
+    if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
+      write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
+      write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
+      write(*,*) '#### PARTICLE AND GAS.                       ####'
+      write(*,*) '#### SPECIES NUMBER',aspecnumb
+      stop
     endif
-  endif
-
-  if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
-  if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
-
-  if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
-    write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
-    write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
-    write(*,*) '#### PARTICLE AND GAS.                       ####'
-    write(*,*) '#### SPECIES NUMBER',aspecnumb
-    stop
-  endif
 20   continue
 
@@ -242,17 +135,13 @@
   ! Read in daily and day-of-week variation of emissions, if available
   !*******************************************************************
-  ! HSO: This is not yet implemented as namelist parameters
-  ! default values set to 1
 
-  do j=1,24           ! initialize everything to no variation
-    area_hour(i,j)=1.
-    point_hour(i,j)=1.
-  end do
-  do j=1,7
-    area_dow(i,j)=1.
-    point_dow(i,j)=1.
-  end do
-
-  if (readerror.ne.0) then ! text format input
+    do j=1,24           ! initialize everything to no variation
+      area_hour(i,j)=1.
+      point_hour(i,j)=1.
+    end do
+    do j=1,7
+      area_dow(i,j)=1.
+      point_dow(i,j)=1.
+    end do
 
     read(unitspecies,*,end=22)
@@ -265,16 +154,7 @@
     end do
 
-  endif
+22   close(unitspecies)
 
-22 close(unitspecies)
-
-  ! namelist output if requested
-  if (nmlout.eqv..true.) then
-    open(unitspecies,file=path(2)(1:length(2))//'SPECIES_'//aspecnumb//'.namelist',access='append',status='new',err=1000)
-    write(unitspecies,nml=species_params)
-    close(unitspecies)
-  endif
-
-  return
+   return
 
 996   write(*,*) '#####################################################'
@@ -305,8 +185,4 @@
   stop
 
-1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "SPECIES_',aspecnumb,'.namelist'
-  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(2)(1:length(2))
-  stop
 
 end subroutine readspecies
Index: /trunk/src/readwind.f90
===================================================================
--- /trunk/src/readwind.f90	(revision 30)
+++ /trunk/src/readwind.f90	(revision 20)
@@ -71,5 +71,5 @@
   integer :: iret
   integer :: igrib
-  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
+  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
   integer :: gotGrid
   !HSO  end
@@ -90,9 +90,9 @@
   integer :: isec1(56),isec2(22+nxmax+nymax)
   real(kind=4) :: zsec4(jpunp)
-  real(kind=4) :: xaux,yaux
+  real(kind=4) :: xaux,yaux,xaux0,yaux0
   real(kind=8) :: xauxin,yauxin
   real,parameter :: eps=1.e-4
   real(kind=4) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
-  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
+  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
 
   logical :: hflswitch,strswitch
@@ -101,4 +101,9 @@
   character(len=24) :: gribErrorMsg = 'Error reading grib file'
   character(len=20) :: gribFunction = 'readwind'
+
+  !HSO conversion of ECMWF etadot to etadot*dp/deta
+  logical :: etacon=.false.
+  real,parameter :: p00=101325.
+  real :: dak,dbk
 
   hflswitch=.false.
@@ -150,6 +155,4 @@
   endif
 
-  conversion_factor=1.
-
   else
 
@@ -165,6 +168,4 @@
   call grib_check(iret,gribFunction,gribErrorMsg)
   call grib_get_int(igrib,'level',valSurf,iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'paramId',parId,iret)
   call grib_check(iret,gribFunction,gribErrorMsg)
 
@@ -176,5 +177,4 @@
   isec1(8)=-1
   isec1(8)=valSurf     ! level
-  conversion_factor=1.
   if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
     isec1(6)=130         ! indicatorOfParameter
@@ -188,6 +188,4 @@
     isec1(6)=134         ! indicatorOfParameter
   elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
-    isec1(6)=135         ! indicatorOfParameter
-  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
     isec1(6)=135         ! indicatorOfParameter
   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
@@ -203,23 +201,21 @@
   elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
     isec1(6)=141         ! indicatorOfParameter
-    conversion_factor=1000.
-  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
+  elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
     isec1(6)=164         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
+  elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
     isec1(6)=142         ! indicatorOfParameter
   elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
     isec1(6)=143         ! indicatorOfParameter
-    conversion_factor=1000.
   elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
     isec1(6)=146         ! indicatorOfParameter
   elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
     isec1(6)=176         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
+  elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
     isec1(6)=180         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
+  elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
     isec1(6)=181         ! indicatorOfParameter
   elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
     isec1(6)=129         ! indicatorOfParameter
-  elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
+  elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
     isec1(6)=160         ! indicatorOfParameter
   elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
@@ -227,10 +223,6 @@
     isec1(6)=172         ! indicatorOfParameter
   else
-    print*,'***WARNING: undefined GRiB2 message found!',discipl, &
+    print*,'***ERROR: undefined GRiB2 message found!',discipl, &
          parCat,parNum,typSurf
-  endif
-  if(parId .ne. isec1(6) .and. parId .ne. 77) then
-    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
-!    stop
   endif
 
@@ -245,9 +237,12 @@
   !HSO  get the required fields from section 2 in a gribex compatible manner
   if (ifield.eq.1) then
-  call grib_get_int(igrib,'numberOfPointsAlongAParallel',isec2(2),iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'numberOfPointsAlongAMeridian',isec2(3),iret)
-  call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'numberOfVerticalCoordinateValues',isec2(12))
+  call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
+       isec2(2),iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+  call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
+       isec2(3),iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+  call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
+       isec2(12))
   call grib_check(iret,gribFunction,gribErrorMsg)
   ! CHECK GRID SPECIFICATIONS
@@ -255,9 +250,9 @@
   if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
   if(isec2(12)/2-1.ne.nlev_ec) &
-  stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
+       stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
   endif ! ifield
 
   !HSO  get the second part of the grid dimensions only from GRiB1 messages
-  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
+  if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
     call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
          xauxin,iret)
@@ -266,14 +261,16 @@
          yauxin,iret)
     call grib_check(iret,gribFunction,gribErrorMsg)
-    if (xauxin.gt.180.) xauxin=xauxin-360.0
-    if (xauxin.lt.-180.) xauxin=xauxin+360.0
-
     xaux=xauxin+real(nxshift)*dx
     yaux=yauxin
-    if (xaux.gt.180.) xaux=xaux-360.0
-    if(abs(xaux-xlon0).gt.eps) &
-    stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
-    if(abs(yaux-ylat0).gt.eps) &
-    stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
+    xaux0=xlon0
+    yaux0=ylat0
+    if(xaux.lt.0.) xaux=xaux+360.
+    if(yaux.lt.0.) yaux=yaux+360.
+    if(xaux0.lt.0.) xaux0=xaux0+360.
+    if(yaux0.lt.0.) yaux0=yaux0+360.
+    if(abs(xaux-xaux0).gt.eps) &
+         stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
+    if(abs(yaux-yaux0).gt.eps) &
+         stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
     gotGrid=1
   endif ! gotGrid
@@ -301,5 +298,5 @@
            zsec4(nxfield*(ny-j-1)+i+1)
       if(isec1(6).eq.141) sd(i,j,1,n)= &!! SNOW DEPTH
-           zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
+           zsec4(nxfield*(ny-j-1)+i+1)
       if(isec1(6).eq.151) msl(i,j,1,n)= &!! SEA LEVEL PRESS.
            zsec4(nxfield*(ny-j-1)+i+1)
@@ -319,5 +316,5 @@
       endif
       if(isec1(6).eq.143) then                      !! CONVECTIVE PREC.
-        convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
+        convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
         if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
       endif
@@ -369,4 +366,20 @@
       do j=0,nymin1
         wwh(i,j,nlev_ec+1)=0.
+      end do
+    end do
+  endif
+
+  ! convert from ECMWF etadot to etadot*dp/deta as needed by FLEXPART
+  if(etacon.eqv..true.) then
+    do k=1,nwzmax
+      dak=akm(k+1)-akm(k)
+      dbk=bkm(k+1)-bkm(k)
+      do i=0,nxmin1
+        do j=0,nymin1 
+          wwh(i,j,k)=2*wwh(i,j,k)*ps(i,j,1,n)*(dak/ps(i,j,1,n)+dbk)/(dak/p00+dbk)
+          if (k.gt.1) then 
+            wwh(i,j,k)=wwh(i,j,k)-wwh(i,j,k-1)
+          endif
+        end do
       end do
     end do
@@ -466,3 +479,2 @@
 
 end subroutine readwind
-
Index: /trunk/src/readwind_nests.f90
===================================================================
--- /trunk/src/readwind_nests.f90	(revision 30)
+++ /trunk/src/readwind_nests.f90	(revision 20)
@@ -51,5 +51,4 @@
   integer :: igrib
   integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
-  integer :: parId !!added by mc for making it consistent with new readwind.f90
   integer :: gotGrid
   !HSO  end
@@ -70,10 +69,8 @@
   integer :: isec1(56),isec2(22+nxmaxn+nymaxn)
   real(kind=4) :: zsec4(jpunp)
-  real(kind=4) :: xaux,yaux
   real(kind=8) :: xauxin,yauxin
-  real,parameter :: eps=1.e-4
+  real(kind=4) :: xaux,yaux,xaux0,yaux0
   real :: ewss(0:nxmaxn-1,0:nymaxn-1),nsss(0:nxmaxn-1,0:nymaxn-1)
   real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
-  real :: conversion_factor !added by mc to make it consistent with new gridchek.f90
 
   logical :: hflswitch,strswitch
@@ -137,7 +134,4 @@
   endif
 
-  conversion_factor=1.
-
-
   else
 
@@ -154,6 +148,4 @@
   call grib_get_int(igrib,'level',valSurf,iret)
   call grib_check(iret,gribFunction,gribErrorMsg)
-  call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consisitent with new readwind.f90
-  call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consisitent with new readwind.f90
 
   !print*,discipl,parCat,parNum,typSurf,valSurf
@@ -164,5 +156,4 @@
   isec1(8)=-1
   isec1(8)=valSurf     ! level
-   conversion_factor=1.
   if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
     isec1(6)=130         ! indicatorOfParameter
@@ -175,8 +166,6 @@
   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
     isec1(6)=134         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot !
+  elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
     isec1(6)=135         ! indicatorOfParameter
-  elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added by mc to make it consisitent with new readwind.f90
-    isec1(6)=135         ! indicatorOfParameter    !added by mc to make it consisitent with new readwind.f90
   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
     isec1(6)=151         ! indicatorOfParameter
@@ -191,23 +180,21 @@
   elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
     isec1(6)=141         ! indicatorOfParameter
-    conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
-  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
+  elseif ((parCat.eq.6).and.(parNum.eq.1)) then ! CC
     isec1(6)=164         ! indicatorOfParameter
-  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
+  elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
     isec1(6)=142         ! indicatorOfParameter
   elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
     isec1(6)=143         ! indicatorOfParameter
-    conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
   elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
     isec1(6)=146         ! indicatorOfParameter
   elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
     isec1(6)=176         ! indicatorOfParameter
-  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
+  elseif ((parCat.eq.2).and.(parNum.eq.17)) then ! EWSS
     isec1(6)=180         ! indicatorOfParameter
-  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
+  elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
     isec1(6)=181         ! indicatorOfParameter
   elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
     isec1(6)=129         ! indicatorOfParameter
-   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
+  elseif ((parCat.eq.3).and.(parNum.eq.7)) then ! SDO
     isec1(6)=160         ! indicatorOfParameter
   elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
@@ -215,10 +202,6 @@
     isec1(6)=172         ! indicatorOfParameter
   else
-    print*,'***WARNING: undefined GRiB2 message found!',discipl, &
+    print*,'***ERROR: undefined GRiB2 message found!',discipl, &
          parCat,parNum,typSurf
-  endif
-  if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consisitent with new readwind.f90
-    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)  !
-!    stop
   endif
 
@@ -244,13 +227,13 @@
   ! CHECK GRID SPECIFICATIONS
   if(isec2(2).ne.nxn(l)) stop &
-  'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
+       'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
   if(isec2(3).ne.nyn(l)) stop &
-  'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
+       'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
   if(isec2(12)/2-1.ne.nlev_ec) stop 'READWIND: VERTICAL DISCRET&
-  IZATION NOT CONSISTENT FOR A NESTING LEVEL'
+       &IZATION NOT CONSISTENT FOR A NESTING LEVEL'
   endif ! ifield
 
   !HSO  get the second part of the grid dimensions only from GRiB1 messages
- if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consisitent with new readwind.f90
+  if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
     call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
          xauxin,iret)
@@ -259,13 +242,18 @@
          yauxin,iret)
     call grib_check(iret,gribFunction,gribErrorMsg)
-    if (xauxin.gt.180.) xauxin=xauxin-360.0
-    if (xauxin.lt.-180.) xauxin=xauxin+360.0
-
     xaux=xauxin
     yaux=yauxin
-    if (abs(xaux-xlon0n(l)).gt.eps) &
-    stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT FOR A NESTING LEVEL'
-    if (abs(yaux-ylat0n(l)).gt.eps) &
-    stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT FOR A NESTING LEVEL'
+    xaux0=xlon0n(l)
+    yaux0=ylat0n(l)
+    if(xaux.lt.0.) xaux=xaux+360.
+    if(yaux.lt.0.) yaux=yaux+360.
+    if(xaux0.lt.0.) xaux0=xaux0+360.
+    if(yaux0.lt.0.) yaux0=yaux0+360.
+    if(xaux.ne.xaux0) &
+         stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT FOR A NES&
+         &TING LEVEL'
+    if(yaux.ne.yaux0) &
+         stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT FOR A NEST&
+         &ING LEVEL'
     gotGrid=1
   endif
@@ -292,5 +280,5 @@
              zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
         if(isec1(6).eq.141) sdn(i,j,1,n,l)= &!! SNOW DEPTH
-             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90!
+             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
         if(isec1(6).eq.151) msln(i,j,1,n,l)= &!! SEA LEVEL PRESS.
              zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
@@ -310,5 +298,5 @@
         endif
         if(isec1(6).eq.143) then                         !! CONVECTIVE PREC.
-          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
+          convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
           if (convprecn(i,j,1,n,l).lt.0.) convprecn(i,j,1,n,l)=0.
         endif
Index: /trunk/src/timemanager.f90
===================================================================
--- /trunk/src/timemanager.f90	(revision 30)
+++ /trunk/src/timemanager.f90	(revision 20)
@@ -43,14 +43,12 @@
   !     call convection BEFORE new fields are read in BWD mode                 *
   !  Changes Caroline Forster, Feb 2005                                        *
-  !   new interface between flexpart and convection scheme                     *
-  !   Emanuel's latest subroutine convect43c.f is used                         *
-  !  Changes Stefan Henne, Harald Sodemann, 2013-2014                          *
-  !   added netcdf output code                                                 *
+  !new interface between flexpart and convection scheme                        *
+  !Emanuel's latest subroutine convect43c.f is used                            *
   !*****************************************************************************
   !                                                                            *
   ! Variables:                                                                 *
-  ! dep                .true. if either wet or dry deposition is switched on   *
+  ! DEP                .true. if either wet or dry deposition is switched on   *
   ! decay(maxspec) [1/s] decay constant for radioactive decay                  *
-  ! drydep             .true. if dry deposition is switched on                 *
+  ! DRYDEP             .true. if dry deposition is switched on                 *
   ! ideltas [s]        modelling period                                        *
   ! itime [s]          actual temporal position of calculation                 *
@@ -72,5 +70,5 @@
   ! prob               probability of absorption at ground due to dry          *
   !                    deposition                                              *
-  ! wetdep             .true. if wet deposition is switched on                 *
+  ! WETDEP             .true. if wet deposition is switched on                 *
   ! weight             weight for each concentration sample (1/2 or 1)         *
   ! uap(maxpart),ucp(maxpart),uzp(maxpart) = random velocities due to          *
@@ -94,7 +92,4 @@
   use par_mod
   use com_mod
-#ifdef NETCDF_OUTPUT
-  use netcdf_output_mod, only: concoutput_netcdf, concoutput_nest_netcdf,concoutput_surf_netcdf, concoutput_surf_nest_netcdf
-#endif 
 
   implicit none
@@ -134,7 +129,4 @@
 
 
-  itime=0
-  !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
-  write(*,46) float(itime)/3600,itime,numpart
   if (verbosity.gt.0) then
     write (*,*) 'timemanager> starting simulation'
@@ -146,4 +138,5 @@
 
   do itime=0,ideltas,lsynctime
+
 
   ! Computation of wet deposition, OH reaction and mass transfer
@@ -157,5 +150,5 @@
   !********************************************************************
 
-    if (wetdep .and. itime .ne. 0 .and. numpart .gt. 0) then
+    if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
         if (verbosity.gt.0) then
            write (*,*) 'timemanager> call wetdepo'
@@ -164,8 +157,8 @@
     endif
 
-    if (ohrea .and. itime .ne. 0 .and. numpart .gt. 0) &
+    if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) &
          call ohreaction(itime,lsynctime,loutnext)
 
-    if (assspec .and. itime .ne. 0 .and. numpart .gt. 0) then
+    if (ASSSPEC .and. itime .ne. 0 .and. numpart .gt. 0) then
        stop 'associated species not yet implemented!'
   !     call transferspec(itime,lsynctime,loutnext)
@@ -245,5 +238,5 @@
   !***********************************************************************
 
-    if (dep.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
+    if (DEP.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
       do ks=1,nspec
       do kp=1,maxpointspec_act
@@ -355,56 +348,31 @@
         if ((iout.le.3.).or.(iout.eq.5)) then
           if (surf_only.ne.1) then 
-            if (lnetcdfout.eq.1) then 
-#ifdef NETCDF_OUTPUT
-              call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
-#endif 
-            else 
-              call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
-            endif
+          call concoutput(itime,outnum,gridtotalunc, &
+               wetgridtotalunc,drygridtotalunc)
           else  
-            if (verbosity.eq.1) then
-             print*,'call concoutput_surf '
-             call system_clock(count_clock)
-             write(*,*) 'system clock',count_clock - count_clock0   
-            endif
-            if (lnetcdfout.eq.1) then 
-#ifdef NETCDF_OUTPUT
-              call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
-#endif 
-            else
-              call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
-              if (verbosity.eq.1) then
-                print*,'called concoutput_surf '
-                call system_clock(count_clock)
-                write(*,*) 'system clock',count_clock - count_clock0   
-              endif
-            endif
+  if (verbosity.eq.1) then
+     print*,'call concoutput_surf '
+     CALL SYSTEM_CLOCK(count_clock)
+     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
+  endif
+          call concoutput_surf(itime,outnum,gridtotalunc, &
+               wetgridtotalunc,drygridtotalunc)
+  if (verbosity.eq.1) then
+     print*,'called concoutput_surf '
+     CALL SYSTEM_CLOCK(count_clock)
+     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
+  endif
           endif
 
-          if (nested_output .eq. 1) then
-            if (lnetcdfout.eq.0) then
-              if (surf_only.ne.1) then
-                call concoutput_nest(itime,outnum)
-              else 
-                call concoutput_surf_nest(itime,outnum)
-              endif
-            else
-#ifdef NETCDF_OUTPUT
-              if (surf_only.ne.1) then
-                call concoutput_nest_netcdf(itime,outnum)
-              else 
-                call concoutput_surf_nest_netcdf(itime,outnum)
-              endif
-#endif
-            endif
-          endif
+          if ((nested_output.eq.1).and.(surf_only.ne.1)) call concoutput_nest(itime,outnum)
+          if ((nested_output.eq.1).and.(surf_only.eq.1)) call concoutput_surf_nest(itime,outnum)
           outnum=0.
         endif
         if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
         if (iflux.eq.1) call fluxoutput(itime)
-        !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
-        write(*,46) float(itime)/3600,itime,numpart
-45      format(i9,' SECONDS SIMULATED: ',i8, ' PARTICLES:    Uncertainty: ',3f7.3)
-46      format(' Simulated ',f7.1,' hours (',i9,' s), ',i8, ' particles')
+        write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc, &
+             drygridtotalunc
+45      format(i9,' SECONDS SIMULATED: ',i8, &
+             ' PARTICLES:    Uncertainty: ',3f7.3)
         if (ipout.ge.1) call partoutput(itime)    ! dump particle positions
         loutnext=loutnext+loutstep
@@ -504,6 +472,6 @@
   !****************************
 
-        xold=real(xtra1(j))
-        yold=real(ytra1(j))
+        xold=xtra1(j)
+        yold=ytra1(j)
         zold=ztra1(j)
 
@@ -545,5 +513,5 @@
             endif
 
-            if (drydepspec(ks)) then        ! dry deposition
+            if (DRYDEPSPEC(ks)) then        ! dry deposition
               drydeposit(ks)=xmass1(j,ks)*prob(ks)*decfact
               xmass1(j,ks)=xmass1(j,ks)*(1.-prob(ks))*decfact
@@ -556,4 +524,5 @@
             endif
 
+
             if (mdomainfill.eq.0) then
               if (xmass(npoint(j),ks).gt.0.) &
@@ -567,16 +536,14 @@
           if (xmassfract.lt.0.0001) then   ! terminate all particles carrying less mass
             itra1(j)=-999999999
-            if (verbosity.gt.0) then
-              print*,'terminated particle ',j,' for small mass'
-            endif
           endif
 
   !        Sabine Eckhardt, June 2008
   !        don't create depofield for backward runs
-          if (drydep.AND.(ldirect.eq.1)) then
-            call drydepokernel(nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)),nage,kp)
-            if (nested_output.eq.1) then
-              call drydepokernel_nest(nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)),nage,kp)
-            endif
+          if (DRYDEP.AND.(ldirect.eq.1)) then
+            call drydepokernel(nclass(j),drydeposit,real(xtra1(j)), &
+                 real(ytra1(j)),nage,kp)
+            if (nested_output.eq.1) call drydepokernel_nest( &
+                 nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)), &
+                 nage,kp)
           endif
 
@@ -585,9 +552,7 @@
 
           if (abs(itra1(j)-itramem(j)).ge.lage(nageclass)) then
-            if (linit_cond.ge.1) call initial_cond_calc(itime+lsynctime,j)
+            if (linit_cond.ge.1) &
+                 call initial_cond_calc(itime+lsynctime,j)
             itra1(j)=-999999999
-            if (verbosity.gt.0) then
-              print*,'terminated particle ',j,' for age'
-            endif
           endif
         endif
@@ -617,11 +582,11 @@
 
   if (iflux.eq.1) then
-    deallocate(flux)
+      deallocate(flux)
   endif
-  if (ohrea.eqv..TRUE.) then
-    deallocate(OH_field,OH_field_height)
+  if (OHREA.eqv..TRUE.) then
+      deallocate(OH_field,OH_field_height)
   endif
   if (ldirect.gt.0) then
-    deallocate(drygridunc,wetgridunc)
+  deallocate(drygridunc,wetgridunc)
   endif
   deallocate(gridunc)
@@ -630,11 +595,11 @@
   deallocate(xmasssave)
   if (nested_output.eq.1) then
-    deallocate(orooutn, arean, volumen)
-    if (ldirect.gt.0) then
-      deallocate(griduncn,drygriduncn,wetgriduncn)
-    endif
+     deallocate(orooutn, arean, volumen)
+     if (ldirect.gt.0) then
+     deallocate(griduncn,drygriduncn,wetgriduncn)
+     endif
   endif
   deallocate(outheight,outheighthalf)
-  deallocate(oroout,area,volume)
+  deallocate(oroout, area, volume)
 
 end subroutine timemanager
Index: /trunk/src/verttransform.f90
===================================================================
--- /trunk/src/verttransform.f90	(revision 30)
+++ /trunk/src/verttransform.f90	(revision 20)
@@ -21,5 +21,5 @@
 
 subroutine verttransform(n,uuh,vvh,wwh,pvh)
-  !                      i  i   i   i   i
+  !                         i  i   i   i   i
   !*****************************************************************************
   !                                                                            *
@@ -49,4 +49,6 @@
   ! Sabine Eckhardt, March 2007
   ! added the variable cloud for use with scavenging - descr. in com_mod
+  ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
+  ! note that also other subroutines are affected by the fix
   !*****************************************************************************
   !                                                                            *
@@ -70,11 +72,13 @@
 
   integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
-  integer :: rain_cloud_above,kz_inv
+ integer :: rain_cloud_above,kz_inv !SE
+  integer icloudtop !PS
   real :: f_qvsat,pressure
-  real :: rh,lsp,convp
-  real :: rhoh(nuvzmax),pinmconv(nzmax)
+ !real :: rh,lsp,convp
+  real :: rh,lsp,convp,prec,rhmin  
+  real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
   real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
   real :: xlon,ylat,xlonr,dzdx,dzdy
-  real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
+  real :: dzdx1,dzdx2,dzdy1,dzdy2, precmin
   real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
   real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
@@ -83,5 +87,7 @@
   real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
   real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
+  logical lconvectprec
   real,parameter :: const=r_air/ga
+  parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
 
   logical :: init = .true.
@@ -91,6 +97,15 @@
   ! If verttransform is called the first time, initialize heights of the   *
   ! z levels in meter. The heights are the heights of model levels, where  *
-  ! u,v,T and qv are given.                                                *
+  ! u,v,T and qv are given, and of the interfaces, where w is given. So,   *
+  ! the vertical resolution in the z system is doubled. As reference point,*
+  ! the lower left corner of the grid is used.                             *
+  ! Unlike in the eta system, no difference between heights for u,v and    *
+  ! heights for w exists.                                                  *
   !*************************************************************************
+
+
+  !      do 897 kz=1,nuvz
+  !       write (*,*) 'akz: ',akz(kz),'bkz',bkz(kz)
+  !897     continue
 
   if (init) then
@@ -98,5 +113,5 @@
   ! Search for a point with high surface pressure (i.e. not above significant topography)
   ! Then, use this point to construct a reference z profile, to be used at all times
-  !**************************************************************************************
+  !*****************************************************************************
 
     do jy=0,nymin1
@@ -113,5 +128,5 @@
 
     tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
-    ps(ixm,jym,1,n))
+         ps(ixm,jym,1,n))
     pold=ps(ixm,jym,1,n)
     height(1)=0.
@@ -121,14 +136,46 @@
       tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
 
+
+  ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
+  ! upon the transformation to z levels. In order to save computer memory, this is
+  ! not done anymore in the standard version. However, this option can still be
+  ! switched on by replacing the following lines with those below, that are
+  ! currently commented out.
+  ! Note that two more changes are necessary in this subroutine below.
+  ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
+  !*****************************************************************************
+
       if (abs(tv-tvold).gt.0.2) then
-        height(kz)=height(kz-1)+const*log(pold/pint)* &
-        (tv-tvold)/log(tv/tvold)
+        height(kz)= &
+             height(kz-1)+const*log(pold/pint)* &
+             (tv-tvold)/log(tv/tvold)
       else
-        height(kz)=height(kz-1)+const*log(pold/pint)*tv
+        height(kz)=height(kz-1)+ &
+             const*log(pold/pint)*tv
       endif
+
+  ! Switch on following lines to use doubled vertical resolution
+  !*************************************************************
+  !    if (abs(tv-tvold).gt.0.2) then
+  !      height((kz-1)*2)=
+  !    +      height(max((kz-2)*2,1))+const*log(pold/pint)*
+  !    +      (tv-tvold)/log(tv/tvold)
+  !    else
+  !      height((kz-1)*2)=height(max((kz-2)*2,1))+
+  !    +      const*log(pold/pint)*tv
+  !    endif
+  ! End doubled vertical resolution
 
       tvold=tv
       pold=pint
     end do
+
+
+  ! Switch on following lines to use doubled vertical resolution
+  !*************************************************************
+  !  do 7 kz=3,nz-1,2
+  !    height(kz)=0.5*(height(kz-1)+height(kz+1))
+  !  height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
+  ! End doubled vertical resolution
 
 
@@ -157,7 +204,8 @@
   do jy=0,nymin1
     do ix=0,nxmin1
-      tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ps(ix,jy,1,n))
+      tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
+           ps(ix,jy,1,n))
       pold=ps(ix,jy,1,n)
-      uvwzlev(ix,jy,1)=0.
+      uvzlev(1)=0.
       wzlev(1)=0.
       rhoh(1)=pold/(r_air*tvold)
@@ -173,8 +221,8 @@
 
         if (abs(tv-tvold).gt.0.2) then
-          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
-          (tv-tvold)/log(tv/tvold)
+          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
+               (tv-tvold)/log(tv/tvold)
         else
-          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
+          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
         endif
 
@@ -185,21 +233,35 @@
 
       do kz=2,nwz-1
-        wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2.
-      end do
-      wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1)
+        wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
+      end do
+      wzlev(nwz)=wzlev(nwz-1)+ &
+           uvzlev(nuvz)-uvzlev(nuvz-1)
+
+      uvwzlev(ix,jy,1)=0.0
+      do kz=2,nuvz
+        uvwzlev(ix,jy,kz)=uvzlev(kz)
+      end do
+
+  ! Switch on following lines to use doubled vertical resolution
+  ! Switch off the three lines above.
+  !*************************************************************
+  !22          uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
+  !     do 23 kz=2,nwz
+  !23          uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
+  ! End doubled vertical resolution
 
   ! pinmconv=(h2-h1)/(p2-p1)
 
       pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
-      ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
-      (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
+           ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
+           (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
       do kz=2,nz-1
         pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
-        ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
-        (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
+             ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
+             (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
       end do
       pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
-      ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
-      (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
+           ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
+           (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
 
   ! Levels, where u,v,t and q are given
@@ -221,5 +283,5 @@
       do iz=2,nz-1
         do kz=kmin,nuvz
-          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
+          if(height(iz).gt.uvzlev(nuvz)) then
             uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
             vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
@@ -230,8 +292,8 @@
             goto 30
           endif
-          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
-          (height(iz).le.uvwzlev(ix,jy,kz))) then
-           dz1=height(iz)-uvwzlev(ix,jy,kz-1)
-           dz2=uvwzlev(ix,jy,kz)-height(iz)
+          if ((height(iz).gt.uvzlev(kz-1)).and. &
+               (height(iz).le.uvzlev(kz))) then
+           dz1=height(iz)-uvzlev(kz-1)
+           dz2=uvzlev(kz)-height(iz)
            dz=dz1+dz2
            uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
@@ -277,8 +339,8 @@
 
       drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
-      (height(2)-height(1))
+           (height(2)-height(1))
       do kz=2,nz-1
         drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
-        (height(kz+1)-height(kz-1))
+             (height(kz+1)-height(kz-1))
       end do
       drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
@@ -294,5 +356,4 @@
 
   do jy=1,ny-2
-    cosf=cos((real(jy)*dy+ylat0)*pi180)
     do ix=1,nx-2
 
@@ -300,10 +361,10 @@
       do iz=2,nz-1
 
-        ui=uu(ix,jy,iz,n)*dxconst/cosf
+        ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
         vi=vv(ix,jy,iz,n)*dyconst
 
         do kz=kmin,nz
           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
-          (height(iz).le.uvwzlev(ix,jy,kz))) then
+               (height(iz).le.uvwzlev(ix,jy,kz))) then
             dz1=height(iz)-uvwzlev(ix,jy,kz-1)
             dz2=uvwzlev(ix,jy,kz)-height(iz)
@@ -348,5 +409,6 @@
         do iz=1,nz
           call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
-          vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
+               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+               vvpol(ix,jy,iz,n))
         end do
       end do
@@ -362,10 +424,12 @@
       xlon=xlon0+real(nx/2-1)*dx
       xlonr=xlon*pi/180.
-      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2)
+      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
+           vv(nx/2-1,nymin1,iz,n)**2)
       if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
-        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
+        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
+             vv(nx/2-1,nymin1,iz,n))-xlonr
       else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
         ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
-        vv(nx/2-1,nymin1,iz,n))-xlonr
+             vv(nx/2-1,nymin1,iz,n))-xlonr
       else
         ddpol=pi/2-xlonr
@@ -380,5 +444,7 @@
       uuaux=-ffpol*sin(xlonr+ddpol)
       vvaux=-ffpol*cos(xlonr+ddpol)
-      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+           vvpolaux)
+
       jy=nymin1
       do ix=0,nxmin1
@@ -419,5 +485,6 @@
         do iz=1,nz
           call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
-          vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
+               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+               vvpol(ix,jy,iz,n))
         end do
       end do
@@ -432,9 +499,12 @@
       xlon=xlon0+real(nx/2-1)*dx
       xlonr=xlon*pi/180.
-      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2)
+      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
+           vv(nx/2-1,0,iz,n)**2)
       if (vv(nx/2-1,0,iz,n).lt.0.) then
-        ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
+        ddpol=atan(uu(nx/2-1,0,iz,n)/ &
+             vv(nx/2-1,0,iz,n))+xlonr
       else if (vv(nx/2-1,0,iz,n).gt.0.) then
-        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
+        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
+             vv(nx/2-1,0,iz,n))+xlonr
       else
         ddpol=pi/2-xlonr
@@ -449,5 +519,6 @@
       uuaux=+ffpol*sin(xlonr-ddpol)
       vvaux=-ffpol*cos(xlonr-ddpol)
-      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+           vvpolaux)
 
       jy=0
@@ -480,40 +551,135 @@
   !   create a cloud and rainout/washout field, clouds occur where rh>80%
   !   total cloudheight is stored at level 0
+
+
+
   do jy=0,nymin1
     do ix=0,nxmin1
-      rain_cloud_above=0
-      lsp=lsprec(ix,jy,1,n)
-      convp=convprec(ix,jy,1,n)
-      cloudsh(ix,jy,n)=0
-      do kz_inv=1,nz-1
-         kz=nz-kz_inv+1
-         pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
-         rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
-         clouds(ix,jy,kz,n)=0
-         if (rh.gt.0.8) then ! in cloud
-            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
-               rain_cloud_above=1
-               cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
-               height(kz)-height(kz-1)
-               if (lsp.ge.convp) then
-                  clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
-               else
-                  clouds(ix,jy,kz,n)=2 ! convp dominated rainout
-               endif
-            else ! no precipitation
-                  clouds(ix,jy,kz,n)=1 ! cloud
+
+
+
+  !    rain_cloud_above=0
+  !    lsp=lsprec(ix,jy,1,n)
+  !    convp=convprec(ix,jy,1,n)
+  !    cloudsh(ix,jy,n)=0
+  !    do kz_inv=1,nz-1
+  !       kz=nz-kz_inv+1
+  !       pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
+  !       rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
+  !       clouds(ix,jy,kz,n)=0
+  !       if (rh.gt.0.8) then ! in cloud
+  !          if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+  !             rain_cloud_above=1
+  !             cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
+  !                  height(kz)-height(kz-1)
+  !             if (lsp.ge.convp) then
+  !                clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
+  !             else
+  !                clouds(ix,jy,kz,n)=2 ! convp dominated rainout
+  !             endif
+  !          else ! no precipitation
+  !                clouds(ix,jy,kz,n)=1 ! cloud
+  !          endif
+  !       else ! no cloud
+  !          if (rain_cloud_above.eq.1) then ! scavenging
+  !             if (lsp.ge.convp) then
+  !                clouds(ix,jy,kz,n)=5 ! lsp dominated washout
+  !             else
+  !                clouds(ix,jy,kz,n)=4 ! convp dominated washout
+  !             endif
+  !          endif
+  !       endif
+  !    end do
+
+
+   ! PS 3012
+
+             lsp=lsprec(ix,jy,1,n)
+          convp=convprec(ix,jy,1,n)
+          prec=lsp+convp
+          if (lsp.gt.convp) then !  prectype='lsp'
+            lconvectprec = .false.
+          else ! prectype='cp '
+            lconvectprec = .true.
+          endif
+          rhmin = 0.90 ! standard condition for presence of clouds
+!PS       note that original by Sabine Eckhart was 80%
+!PS       however, for T<-20 C we consider saturation over ice
+!PS       so I think 90% should be enough          
+          icloudbot(ix,jy,n)=icmv
+          icloudtop=icmv ! this is just a local variable
+98        do kz=1,nz
+            pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
+            rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
+!ps            if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
+            if (rh .gt. rhmin) then
+              if (icloudbot(ix,jy,n) .eq. icmv) then
+                icloudbot(ix,jy,n)=nint(height(kz))
+              endif
+              icloudtop=nint(height(kz)) ! use int to save memory
             endif
-         else ! no cloud
-            if (rain_cloud_above.eq.1) then ! scavenging
-               if (lsp.ge.convp) then
-                  clouds(ix,jy,kz,n)=5 ! lsp dominated washout
-               else
-                  clouds(ix,jy,kz,n)=4 ! convp dominated washout
-               endif
+          enddo
+
+!PS try to get a cloud thicker than 50 m 
+!PS if there is at least .01 mm/h  - changed to 0.002 and put into
+!PS parameter precpmin        
+          if ((icloudbot(ix,jy,n) .eq. icmv .or. &
+              icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
+              prec .gt. precmin) then
+            rhmin = rhmin - 0.05
+            if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
+          endif
+!PS implement a rough fix for badly represented convection
+!PS is based on looking at a limited set of comparison data
+          if (lconvectprec .and. icloudtop .lt. 6000 .and. &
+             prec .gt. precmin) then 
+
+            if (convp .lt. 0.1) then
+              icloudbot(ix,jy,n) = 500
+              icloudtop =         8000
+            else
+              icloudbot(ix,jy,n) = 0
+              icloudtop =      10000
             endif
-         endif
-      end do
+          endif
+          if (icloudtop .ne. icmv) then
+            icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
+          else
+            icloudthck(ix,jy,n) = icmv
+          endif
+!PS  get rid of too thin clouds      
+          if (icloudthck(ix,jy,n) .lt. 50) then
+            icloudbot(ix,jy,n)=icmv
+            icloudthck(ix,jy,n)=icmv
+          endif
+
+
     end do
   end do
 
+  !do 102 kz=1,nuvz
+  !write(an,'(i02)') kz+10
+  !write(*,*) nuvz,nymin1,nxmin1,'--',an,'--'
+  !open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')
+  !do 101 jy=0,nymin1
+  !    write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)
+  !101   continue
+  ! close(4)
+  !102   continue
+
+  ! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')
+  ! do 103 jy=0,nymin1
+  !     write (4,*)
+  !+       (height(kz),kz=1,nuvz)
+  !103   continue
+  ! close(4)
+
+  !open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')
+  ! do 104 jy=0,nymin1
+  !     write (4,*)
+  !+       (r_air*tt(ix,jy,1,n)*rho(ix,jy,1,n),ix=1,nxmin1)
+  !104   continue
+  ! close(4)
+
+
 end subroutine verttransform
Index: /trunk/src/verttransform_gfs.f90
===================================================================
--- /trunk/src/verttransform_gfs.f90	(revision 30)
+++ /trunk/src/verttransform_gfs.f90	(revision 20)
@@ -21,5 +21,5 @@
 
 subroutine verttransform(n,uuh,vvh,wwh,pvh)
-  !                      i  i   i   i   i
+  !                         i  i   i   i   i
   !*****************************************************************************
   !                                                                            *
@@ -71,8 +71,8 @@
   real :: f_qvsat,pressure
   real :: rh,lsp,convp
-  real :: rhoh(nuvzmax),pinmconv(nzmax)
+  real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
   real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
   real :: xlon,ylat,xlonr,dzdx,dzdy
-  real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
+  real :: dzdx1,dzdx2,dzdy1,dzdy2
   real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
   real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
@@ -92,5 +92,9 @@
   ! If verttransform is called the first time, initialize heights of the   *
   ! z levels in meter. The heights are the heights of model levels, where  *
-  ! u,v,T and qv are given.                                                *
+  ! u,v,T and qv are given, and of the interfaces, where w is given. So,   *
+  ! the vertical resolution in the z system is doubled. As reference point,*
+  ! the lower left corner of the grid is used.                             *
+  ! Unlike in the eta system, no difference between heights for u,v and    *
+  ! heights for w exists.                                                  *
   !*************************************************************************
 
@@ -114,5 +118,5 @@
 
     tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
-    ps(ixm,jym,1,n))
+         ps(ixm,jym,1,n))
     pold=ps(ixm,jym,1,n)
     height(1)=0.
@@ -122,14 +126,46 @@
       tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
 
+
+  ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
+  ! upon the transformation to z levels. In order to save computer memory, this is
+  ! not done anymore in the standard version. However, this option can still be
+  ! switched on by replacing the following lines with those below, that are
+  ! currently commented out.
+  ! Note that two more changes are necessary in this subroutine below.
+  ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
+  !*****************************************************************************
+
       if (abs(tv-tvold).gt.0.2) then
-        height(kz)=height(kz-1)+const*log(pold/pint)* &
-        (tv-tvold)/log(tv/tvold)
+        height(kz)= &
+             height(kz-1)+const*log(pold/pint)* &
+             (tv-tvold)/log(tv/tvold)
       else
-        height(kz)=height(kz-1)+const*log(pold/pint)*tv
+        height(kz)=height(kz-1)+ &
+             const*log(pold/pint)*tv
       endif
+
+  ! Switch on following lines to use doubled vertical resolution
+  !*************************************************************
+  !    if (abs(tv-tvold).gt.0.2) then
+  !      height((kz-1)*2)=
+  !    +      height(max((kz-2)*2,1))+const*log(pold/pint)*
+  !    +      (tv-tvold)/log(tv/tvold)
+  !    else
+  !      height((kz-1)*2)=height(max((kz-2)*2,1))+
+  !    +      const*log(pold/pint)*tv
+  !    endif
+  ! End doubled vertical resolution
 
       tvold=tv
       pold=pint
     end do
+
+
+  ! Switch on following lines to use doubled vertical resolution
+  !*************************************************************
+  !  do 7 kz=3,nz-1,2
+  !    height(kz)=0.5*(height(kz-1)+height(kz+1))
+  !  height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
+  ! End doubled vertical resolution
 
 
@@ -162,5 +198,5 @@
       llev = 0
       do i=1,nuvz
-        if (ps(ix,jy,1,n).lt.akz(i)) llev=i
+       if (ps(ix,jy,1,n).lt.akz(i)) llev=i
       end do
        llev = llev+1
@@ -176,4 +212,5 @@
       tvold=tth(ix,jy,llev,n)*(1.+0.608*qvh(ix,jy,llev,n))
       pold=akz(llev)
+      uvzlev(llev)=0.
       wzlev(llev)=0.
       uvwzlev(ix,jy,llev)=0.
@@ -186,14 +223,24 @@
 
         if (abs(tv-tvold).gt.0.2) then
-          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
-          (tv-tvold)/log(tv/tvold)
+          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
+               (tv-tvold)/log(tv/tvold)
         else
-          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
+          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
         endif
-        wzlev(kz)=uvwzlev(ix,jy,kz)
+        wzlev(kz)=uvzlev(kz)
+        uvwzlev(ix,jy,kz)=uvzlev(kz)
 
         tvold=tv
         pold=pint
       end do
+
+
+  ! Switch on following lines to use doubled vertical resolution
+  ! Switch off the three lines above.
+  !*************************************************************
+  !22          uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
+  !     do 23 kz=2,nwz
+  !23          uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
+  ! End doubled vertical resolution
 
   ! pinmconv=(h2-h1)/(p2-p1)
@@ -232,5 +279,5 @@
       do iz=2,nz-1
         do kz=kmin,nuvz
-          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
+          if(height(iz).gt.uvzlev(nuvz)) then
             uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
             vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
@@ -242,18 +289,18 @@
             goto 30
           endif
-          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
-          (height(iz).le.uvwzlev(ix,jy,kz))) then
-            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
-            dz2=uvwzlev(ix,jy,kz)-height(iz)
-            dz=dz1+dz2
-            uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
-            vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
-            tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
-            +tth(ix,jy,kz,n)*dz1)/dz
-            qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
-            +qvh(ix,jy,kz,n)*dz1)/dz
-            pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
-            rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
-            pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz
+          if ((height(iz).gt.uvzlev(kz-1)).and. &
+               (height(iz).le.uvzlev(kz))) then
+           dz1=height(iz)-uvzlev(kz-1)
+           dz2=uvzlev(kz)-height(iz)
+           dz=dz1+dz2
+           uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
+           vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
+           tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
+                +tth(ix,jy,kz,n)*dz1)/dz
+           qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
+                +qvh(ix,jy,kz,n)*dz1)/dz
+           pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
+           rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
+           pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz
           endif
         end do
@@ -271,10 +318,11 @@
         do kz=kmin,nwz
           if ((height(iz).gt.wzlev(kz-1)).and. &
-          (height(iz).le.wzlev(kz))) then
-            dz1=height(iz)-wzlev(kz-1)
-            dz2=wzlev(kz)-height(iz)
-            dz=dz1+dz2
-            ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
-            +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
+               (height(iz).le.wzlev(kz))) then
+           dz1=height(iz)-wzlev(kz-1)
+           dz2=wzlev(kz)-height(iz)
+           dz=dz1+dz2
+           ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
+                +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
+
           endif
         end do
@@ -289,5 +337,5 @@
       do kz=2,nz-1
         drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
-        (height(kz+1)-height(kz-1))
+             (height(kz+1)-height(kz-1))
       end do
       drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
@@ -303,5 +351,4 @@
 
   do jy=1,ny-2
-    cosf=cos((real(jy)*dy+ylat0)*pi180)
     do ix=1,nx-2
 
@@ -320,10 +367,10 @@
       do iz=2,nz-1
 
-        ui=uu(ix,jy,iz,n)*dxconst/cosf
+        ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
         vi=vv(ix,jy,iz,n)*dyconst
 
         do kz=kmin,nz
           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
-          (height(iz).le.uvwzlev(ix,jy,kz))) then
+               (height(iz).le.uvwzlev(ix,jy,kz))) then
             dz1=height(iz)-uvwzlev(ix,jy,kz-1)
             dz2=uvwzlev(ix,jy,kz)-height(iz)
@@ -379,10 +426,12 @@
       xlon=xlon0+real(nx/2-1)*dx
       xlonr=xlon*pi/180.
-      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2)
-      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
-        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
+      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
+           vv(nx/2-1,nymin1,iz,n)**2)
+      if(vv(nx/2-1,nymin1,iz,n).lt.0.) then
+        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
+             vv(nx/2-1,nymin1,iz,n))-xlonr
       elseif (vv(nx/2-1,nymin1,iz,n).gt.0.) then
         ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
-        vv(nx/2-1,nymin1,iz,n))-xlonr
+             vv(nx/2-1,nymin1,iz,n))-xlonr
       else
         ddpol=pi/2-xlonr
@@ -397,5 +446,7 @@
       uuaux=-ffpol*sin(xlonr+ddpol)
       vvaux=-ffpol*cos(xlonr+ddpol)
-      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+           vvpolaux)
+
       jy=nymin1
       do ix=0,nxmin1
@@ -409,5 +460,5 @@
   ! ward parallel of latitude
 
-    do iz=1,nz
+  do iz=1,nz
       wdummy=0.
       jy=ny-2
@@ -420,5 +471,5 @@
         ww(ix,jy,iz,n)=wdummy
       end do
-    end do
+  end do
 
   endif
@@ -436,5 +487,6 @@
         do iz=1,nz
           call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
-          vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
+               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+               vvpol(ix,jy,iz,n))
         end do
       end do
@@ -446,9 +498,12 @@
       xlon=xlon0+real(nx/2-1)*dx
       xlonr=xlon*pi/180.
-      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2)
+      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
+           vv(nx/2-1,0,iz,n)**2)
       if(vv(nx/2-1,0,iz,n).lt.0.) then
-        ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
+        ddpol=atan(uu(nx/2-1,0,iz,n)/ &
+             vv(nx/2-1,0,iz,n))+xlonr
       elseif (vv(nx/2-1,0,iz,n).gt.0.) then
-        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))-xlonr
+        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
+             vv(nx/2-1,0,iz,n))-xlonr
       else
         ddpol=pi/2-xlonr
@@ -463,5 +518,6 @@
       uuaux=+ffpol*sin(xlonr-ddpol)
       vvaux=-ffpol*cos(xlonr-ddpol)
-      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+           vvpolaux)
 
       jy=0
@@ -506,23 +562,24 @@
          clouds(ix,jy,kz,n)=0
          if (rh.gt.0.8) then ! in cloud
-           if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
-              rain_cloud_above=1
-              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
-              if (lsp.ge.convp) then
-                 clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
-              else
-                 clouds(ix,jy,kz,n)=2 ! convp dominated rainout
-              endif
-           else ! no precipitation
-             clouds(ix,jy,kz,n)=1 ! cloud
-           endif
+            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+               rain_cloud_above=1
+               cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
+                    height(kz)-height(kz-1)
+               if (lsp.ge.convp) then
+                  clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
+               else
+                  clouds(ix,jy,kz,n)=2 ! convp dominated rainout
+               endif
+            else ! no precipitation
+                  clouds(ix,jy,kz,n)=1 ! cloud
+            endif
          else ! no cloud
-           if (rain_cloud_above.eq.1) then ! scavenging
-             if (lsp.ge.convp) then
-               clouds(ix,jy,kz,n)=5 ! lsp dominated washout
-             else
-               clouds(ix,jy,kz,n)=4 ! convp dominated washout
-             endif
-           endif
+            if (rain_cloud_above.eq.1) then ! scavenging
+               if (lsp.ge.convp) then
+                  clouds(ix,jy,kz,n)=5 ! lsp dominated washout
+               else
+                  clouds(ix,jy,kz,n)=4 ! convp dominated washout
+               endif
+            endif
          endif
       end do
Index: /trunk/src/verttransform_nests.f90
===================================================================
--- /trunk/src/verttransform_nests.f90	(revision 30)
+++ /trunk/src/verttransform_nests.f90	(revision 20)
@@ -21,5 +21,5 @@
 
 subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn)
-  !                            i   i    i    i   i
+  !                               i   i    i    i   i
   !*****************************************************************************
   !                                                                            *
@@ -72,7 +72,7 @@
   integer :: rain_cloud_above,kz_inv
   real :: f_qvsat,pressure,rh,lsp,convp
-  real :: wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax)
+  real :: uvzlev(nuvzmax),wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax)
   real :: uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax)
-  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi,cosf
+  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
   real :: dzdx,dzdy
   real :: dzdx1,dzdx2,dzdy1,dzdy2
@@ -96,7 +96,7 @@
 
       tvold=tt2n(ix,jy,1,n,l)*(1.+0.378*ew(td2n(ix,jy,1,n,l))/ &
-      psn(ix,jy,1,n,l))
+           psn(ix,jy,1,n,l))
       pold=psn(ix,jy,1,n,l)
-      uvwzlev(ix,jy,1)=0.
+      uvzlev(1)=0.
       wzlev(1)=0.
       rhoh(1)=pold/(r_air*tvold)
@@ -112,8 +112,8 @@
 
         if (abs(tv-tvold).gt.0.2) then
-          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
-          (tv-tvold)/log(tv/tvold)
+          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
+               (tv-tvold)/log(tv/tvold)
         else
-          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
+          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
         endif
 
@@ -124,21 +124,43 @@
 
       do kz=2,nwz-1
-        wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2.
-      end do
-      wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1)
+        wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
+      end do
+      wzlev(nwz)=wzlev(nwz-1)+ &
+           uvzlev(nuvz)-uvzlev(nuvz-1)
+
+  ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
+  ! upon the transformation to z levels. In order to save computer memory, this is
+  ! not done anymore in the standard version. However, this option can still be
+  ! switched on by replacing the following lines with those below, that are
+  ! currently commented out.
+  ! Note that one change is also necessary in gridcheck.f,
+  ! and three changes in verttransform.f
+  !*****************************************************************************
+      uvwzlev(ix,jy,1)=0.0
+      do kz=2,nuvz
+        uvwzlev(ix,jy,kz)=uvzlev(kz)
+      end do
+
+  ! Switch on following lines to use doubled vertical resolution
+  ! Switch off the three lines above.
+  !*************************************************************
+  !22          uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
+  !     do 23 kz=2,nwz
+  !23          uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
+  ! End doubled vertical resolution
 
   ! pinmconv=(h2-h1)/(p2-p1)
 
       pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
-      ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
-      (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
+           ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
+           (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
       do kz=2,nz-1
         pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
-        ((aknew(kz+1)+bknew(kz+1)*psn(ix,jy,1,n,l))- &
-        (aknew(kz-1)+bknew(kz-1)*psn(ix,jy,1,n,l)))
+             ((aknew(kz+1)+bknew(kz+1)*psn(ix,jy,1,n,l))- &
+             (aknew(kz-1)+bknew(kz-1)*psn(ix,jy,1,n,l)))
       end do
       pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
-      ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
-      (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
+           ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
+           (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
 
 
@@ -161,5 +183,5 @@
       do iz=2,nz-1
         do kz=kmin,nuvz
-          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
+          if(height(iz).gt.uvzlev(nuvz)) then
             uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
             vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
@@ -170,19 +192,19 @@
             goto 30
           endif
-          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
-          (height(iz).le.uvwzlev(ix,jy,kz))) then
-           dz1=height(iz)-uvwzlev(ix,jy,kz-1)
-           dz2=uvwzlev(ix,jy,kz)-height(iz)
+          if ((height(iz).gt.uvzlev(kz-1)).and. &
+               (height(iz).le.uvzlev(kz))) then
+           dz1=height(iz)-uvzlev(kz-1)
+           dz2=uvzlev(kz)-height(iz)
            dz=dz1+dz2
            uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
-           uuhn(ix,jy,kz,l)*dz1)/dz
+                uuhn(ix,jy,kz,l)*dz1)/dz
            vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
-           vvhn(ix,jy,kz,l)*dz1)/dz
+                vvhn(ix,jy,kz,l)*dz1)/dz
            ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
-           tthn(ix,jy,kz,n,l)*dz1)/dz
+                tthn(ix,jy,kz,n,l)*dz1)/dz
            qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
-           qvhn(ix,jy,kz,n,l)*dz1)/dz
+                qvhn(ix,jy,kz,n,l)*dz1)/dz
            pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
-           pvhn(ix,jy,kz,l)*dz1)/dz
+                pvhn(ix,jy,kz,l)*dz1)/dz
            rhon(ix,jy,iz,n,l)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
            kmin=kz
@@ -203,12 +225,12 @@
         do kz=kmin,nwz
           if ((height(iz).gt.wzlev(kz-1)).and. &
-          (height(iz).le.wzlev(kz))) then
-            dz1=height(iz)-wzlev(kz-1)
-            dz2=wzlev(kz)-height(iz)
-            dz=dz1+dz2
-            wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
-            +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
-            kmin=kz
-            goto 40
+               (height(iz).le.wzlev(kz))) then
+           dz1=height(iz)-wzlev(kz-1)
+           dz2=wzlev(kz)-height(iz)
+           dz=dz1+dz2
+           wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
+                +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
+           kmin=kz
+           goto 40
           endif
         end do
@@ -220,8 +242,8 @@
 
       drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
-      (height(2)-height(1))
+           (height(2)-height(1))
       do kz=2,nz-1
         drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
-        rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
+             rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
       end do
       drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
@@ -237,5 +259,4 @@
 
   do jy=1,nyn(l)-2
-    cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
     do ix=1,nxn(l)-2
 
@@ -243,5 +264,6 @@
       do iz=2,nz-1
 
-        ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf
+        ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/ &
+             cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
         vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)
 
@@ -297,5 +319,5 @@
                rain_cloud_above=1
                cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ &
-               height(kz)-height(kz-1)
+                    height(kz)-height(kz-1)
                if (lsp.ge.convp) then
                   cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout
Index: /trunk/src/wetdepo.f90
===================================================================
--- /trunk/src/wetdepo.f90	(revision 30)
+++ /trunk/src/wetdepo.f90	(revision 20)
@@ -21,5 +21,5 @@
 
 subroutine wetdepo(itime,ltsample,loutnext)
-  !                  i      i        i
+  !                     i      i        i
   !*****************************************************************************
   !                                                                            *
@@ -39,4 +39,7 @@
   ! Code may not be correct for decay of deposition!                           *
   !                                                                            *
+  ! Modification by Sabine Eckhart to introduce a new in-/below-cloud
+  ! scheme, not dated
+  ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
   !*****************************************************************************
   !                                                                            *
@@ -71,9 +74,11 @@
 
   integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
-  integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
-  integer :: ks, kp
+  integer :: ngrid,itage,nage,hz,il,interp_time, n , clouds_v !NIK scheme
+  integer :: kz !PS scheme
+  integer :: ks, kp, n1,n2, icbot,ictop, indcloud
+  integer :: scheme_number ! NIK==1, PS ==2
   real :: S_i, act_temp, cl, cle ! in cloud scavenging
   real :: clouds_h ! cloud height for the specific grid point
-  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
+  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav, precsub,f
   real :: wetdeposit(maxspec),restmass
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
@@ -98,5 +103,4 @@
 
   do jpart=1,numpart
-
     if (itra1(jpart).eq.-999999999) goto 20
     if(ldirect.eq.1)then
@@ -110,5 +114,5 @@
       if (itage.lt.lage(nage)) goto 33
     end do
-33  continue
+33   continue
 
 
@@ -119,10 +123,10 @@
     do j=numbnests,1,-1
       if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
-      (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
+           (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
         ngrid=j
         goto 23
       endif
     end do
-23  continue
+23   continue
 
 
@@ -147,42 +151,77 @@
     interp_time=nint(itime-0.5*ltsample)
 
-    if (ngrid.eq.0) then
-      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
-      1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
-      memtime(1),memtime(2),interp_time,lsp,convp,cc)
-    else
-      call interpol_rain_nests(lsprecn,convprecn,tccn, &
-      nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
-      memtime(1),memtime(2),interp_time,lsp,convp,cc)
-    endif
-
-    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
+!    PS nest case still needs to be implemented!!    
+!    if (ngrid.eq.0) then
+!      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
+!           1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
+!           memtime(1),memtime(2),interp_time,lsp,convp,cc)
+           call interpol_rain(lsprec,convprec,tcc,     &
+            icloudbot,icloudthck,nxmax,nymax,1,nx,ny,         &
+            memind,sngl(xtra1(jpart)),sngl(ytra1(jpart)),1,memtime(1), &
+            memtime(2),interp_time,lsp,convp,cc,icbot,ictop,icmv)  
+!    else
+!      call interpol_rain_nests(lsprecn,convprecn,tccn, &
+!           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
+!           memtime(1),memtime(2),interp_time,lsp,convp,cc)
+!    endif
+
+
+ !   if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
+ !PS 2012: subtract a small value, eg 0.01 mm/h, to remove spurious precip
+        prec = lsp+convp
+        precsub = 0.01
+        if (prec .lt. precsub) then
+          goto 20
+        else
+          f = (prec-precsub)/prec
+          lsp = f*lsp
+          convp = f*convp
+        endif
+
+
   ! get the level were the actual particle is in
-    do il=2,nz
-      if (height(il).gt.ztra1(jpart)) then
-        hz=il-1
-        goto 26
-      endif
-    end do
-26  continue
-
-    n=memind(2)
-    if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
-    n=memind(1)
+      do il=2,nz
+        if (height(il).gt.ztra1(jpart)) then
+          !hz=il-1
+          kz=il-1
+          goto 26
+        endif
+      end do
+26     continue
+
+  n=memind(2)
+  if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
+       n=memind(1)
 
   ! if there is no precipitation or the particle is above the clouds no
   ! scavenging is done
 
-    if (ngrid.eq.0) then
-      clouds_v=clouds(ix,jy,hz,n)
-      clouds_h=cloudsh(ix,jy,n)
-    else
-      clouds_v=cloudsn(ix,jy,hz,n,ngrid)
-      clouds_h=cloudsnh(ix,jy,n,ngrid)
-    endif
-  !write(*,*) 'there is
-  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
-    if (clouds_v.le.1) goto 20
-  !write (*,*) 'there is scavenging'
+!old scheme
+!  if (ngrid.eq.0) then
+!     clouds_v=clouds(ix,jy,hz,n)
+!     clouds_h=cloudsh(ix,jy,n)
+!  else
+!     clouds_v=cloudsn(ix,jy,hz,n,ngrid)
+!     clouds_h=cloudsnh(ix,jy,n,ngrid)
+!  endif
+!  !write(*,*) 'there is
+!  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
+!  if (clouds_v.le.1) goto 20
+!  !write (*,*) 'there is scavenging'
+
+  ! PS: part of 2011/2012 fix 
+
+        if (ztra1(jpart) .le. float(ictop)) then
+          if (ztra1(jpart) .gt. float(icbot)) then
+            indcloud = 2 ! in-cloud
+          else
+            indcloud = 1 ! below-cloud
+          endif
+        elseif (ictop .eq. icmv) then
+          indcloud = 0 ! no cloud found, use old scheme
+        else
+          goto 20 ! above cloud
+        endif
+
 
   ! 1) Parameterization of the the area fraction of the grid cell where the
@@ -228,30 +267,28 @@
   !**********************************************************
 
-    do ks=1,nspec      ! loop over species
-      wetdeposit(ks)=0. 
-      wetscav=0.   
-
-  ! NIK 09.12.13: allowed to turn off either below cloud or in-cloud by including TWO separate blocks of if tests
-
-
-  !****i*******************
-  ! BELOW CLOUD SCAVENGING
-  !***********************  
-      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
- !       for aerosols and not highliy soluble substances weta=5E-6
-        wetscav=weta(ks)*prec**wetb(ks)  ! scavenging coefficient 
-      endif !end below-cloud
-
-
-  !********************
-  ! IN CLOUD SCAVENGING
-  !********************
-      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
-
-        if (ngrid.gt.0) then
-          act_temp=ttn(ix,jy,hz,n,ngrid)
-        else
-          act_temp=tt(ix,jy,hz,n)
-        endif
+    do ks=1,nspec            ! loop over species
+      wetdeposit(ks)=0.
+
+    
+    !conflicting changes to the same routine: 1=NIK 2 =PS
+    scheme_number=2   
+    if (scheme_number.eq.1) then !NIK 
+
+      if (weta(ks).gt.0.) then
+        if (clouds_v.ge.4) then
+  !          BELOW CLOUD SCAVENGING
+  !          for aerosols and not highliy soluble substances weta=5E-6
+          wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
+  !        write(*,*) 'bel. wetscav: ',wetscav
+
+        else ! below_cloud clouds_v is lt 4 and gt 1 -> in cloud scavenging
+  !        IN CLOUD SCAVENGING
+  ! BUGFIX tt for nested fields should be ttn
+  ! sec may 2008
+          if (ngrid.gt.0) then
+             act_temp=ttn(ix,jy,hz,n,ngrid)
+          else
+             act_temp=tt(ix,jy,hz,n)
+          endif
 
 ! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening 
@@ -260,70 +297,142 @@
 ! wetc_in=0.9 (default)
 ! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling)
-        cl=weta_in(ks)*prec**wetb_in(ks)
-        if (dquer(ks).gt.0) then ! is particle
-          S_i=wetc_in(ks)/cl
-        else ! is gas
-          cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
-          S_i=1/cle
-        endif
-        wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
-  !     write(*,*) 'in. wetscav:'
-  !     +  ,wetscav,cle,cl,act_temp,prec,clouds_h
-       
-      endif ! end in-cloud
-
-
-      
-  !**************************************************
-  ! CALCULATE DEPOSITION 
-  !**************************************************
-  !     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
-  !     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
-
-      if (wetscav.gt.0.) then
+          cl=weta_in(ks)*prec**wetb_in(ks)
+          if (dquer(ks).gt.0) then ! is particle
+            S_i=wetc_in(ks)/cl
+           else ! is gas
+            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
+            S_i=1/cle
+           endif
+           wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
+  !         write(*,*) 'in. wetscav:'
+  !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h
+        endif
+
+
+  !      if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
+  !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
         wetdeposit(ks)=xmass1(jpart,ks)* &
-        (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! wet deposition
-      else ! if no scavenging
-        wetdeposit(ks)=0.
-      endif
-
-
-      restmass = xmass1(jpart,ks)-wetdeposit(ks)
-      if (ioutputforeachrelease.eq.1) then
-        kp=npoint(jpart)
-      else
-        kp=1
-      endif
-      if (restmass .gt. smallnum) then
-        xmass1(jpart,ks)=restmass
-  !   depostatistic
-  !   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
-  !   depostatistic
-      else
-        xmass1(jpart,ks)=0.
-      endif
-  !   Correct deposited mass to the last time step when radioactive decay of
-  !   gridded deposited mass was calculated
-      if (decay(ks).gt.0.) then
-        wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
-      endif
-
-
-
-    end do !all species
-
-  ! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
+             (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! wet deposition
+  !      new particle mass:
+  !      if (wetdeposit(ks).gt.0) then
+  !         write(*,*) 'wetdepo: ',wetdeposit(ks),ks
+  !      endif
+        restmass = xmass1(jpart,ks)-wetdeposit(ks)
+         if (ioutputforeachrelease.eq.1) then
+           kp=npoint(jpart)
+         else
+           kp=1
+         endif
+        if (restmass .gt. smallnum) then
+          xmass1(jpart,ks)=restmass
+  !ccccccccccccccc depostatistic
+  !       wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
+  !ccccccccccccccc depostatistic
+        else
+          xmass1(jpart,ks)=0.
+        endif
+  ! Correct deposited mass to the last time step when radioactive decay of
+  ! gridded deposited mass was calculated
+       if (decay(ks).gt.0.) then
+          wetdeposit(ks)=wetdeposit(ks) &
+               *exp(abs(ldeltat)*decay(ks))
+       endif
+      else  ! weta(k)
+         wetdeposit(ks)=0.
+      endif ! weta(k)
+
+     elseif (scheme_number.eq.2) then ! PS
+
+!PS          indcloud=0 ! Use this for FOR TESTING, 
+!PS          will skip the new in/below cloud method !!!             
+
+          if (weta(ks).gt.0.) then
+            if (indcloud .eq. 1) then ! BELOW CLOUD SCAVENGING
+!C               for aerosols and not highliy soluble substances weta=5E-6
+              wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
+!c             write(*,*) 'bel. wetscav: ',wetscav
+            elseif (indcloud .eq. 2) then !  IN CLOUD SCAVENGING
+              if (ngrid.gt.0) then
+                 act_temp=ttn(ix,jy,kz,n,ngrid)
+              else
+                 act_temp=tt(ix,jy,kz,n)
+              endif
+
+! from NIK              
+! weta_in=2.0E-07 (default)
+! wetb_in=0.36 (default)
+! wetc_in=0.9 (default)
+
+
+              cl=2E-7*prec**0.36
+              if (dquer(ks).gt.0) then ! is particle
+                S_i=0.9/cl
+              else ! is gas
+                cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
+                S_i=1/cle
+              endif
+              wetscav=S_i*prec/3.6E6/(ictop-icbot) ! 3.6e6 converts mm/h to m/s
+            else ! PS: no cloud diagnosed, old scheme, 
+!CPS          using with fixed a,b for simplicity, one may wish to change!!
+              wetscav = 1.e-4*prec**0.62
+            endif
+
+
+            wetdeposit(ks)=xmass1(jpart,ks)*    &
+            ! (1.-exp(-wetscav*abs(ltsample)))*fraction  ! wet deposition
+             (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! fraction = grfraction (PS) 
+            restmass = xmass1(jpart,ks)-wetdeposit(ks)
+            if (ioutputforeachrelease.eq.1) then
+              kp=npoint(jpart)
+            else
+              kp=1
+            endif
+            if (restmass .gt. smallnum) then
+              xmass1(jpart,ks)=restmass
+!cccccccccccccccc depostatistic
+!c            wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
+!cccccccccccccccc depostatistic
+            else
+              xmass1(jpart,ks)=0.
+            endif
+!C Correct deposited mass to the last time step when radioactive decay of 
+!C gridded deposited mass was calculated
+            if (decay(ks).gt.0.) then
+              wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
+            endif
+          else  ! weta(k)<0
+             wetdeposit(ks)=0.
+          endif 
+
+     endif !on scheme
+
+
+
+    end do
+
+  ! Sabine Eckhard, June 2008 create deposition runs only for forward runs
   ! Add the wet deposition to accumulated amount on output grid and nested output grid
   !*****************************************************************************
 
-    if (ldirect.eq.1) then
-      call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
-      real(ytra1(jpart)),nage,kp)
-      if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
-        wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
-    endif
+!    if (ldirect.eq.1) then
+!    call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
+!         real(ytra1(jpart)),nage,kp)
+!    if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
+!         wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), &
+!         nage,kp)
+!    endif
+
+   !PS
+        if (ldirect.eq.1) then
+          call wetdepokernel(nclass(jpart),wetdeposit,     &
+             sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp)
+          if (nested_output.eq.1)                          &
+           call wetdepokernel_nest(nclass(jpart),wetdeposit, &
+             sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp)  
+        endif
+
 
 20  continue
-  end do ! all particles
+  end do
 
 end subroutine wetdepo
Index: /trunk/src/writeheader_nest.f90
===================================================================
--- /trunk/src/writeheader_nest.f90	(revision 30)
+++ /trunk/src/writeheader_nest.f90	(revision 20)
@@ -67,7 +67,7 @@
 
   if (ldirect.eq.1) then
-    write(unitheader) ibdate,ibtime, trim(flexversion)
+    write(unitheader) ibdate,ibtime,'FLEXPART V8.2'
   else
-    write(unitheader) iedate,ietime, trim(flexversion)
+    write(unitheader) iedate,ietime,'FLEXPART V8.2'
   endif
 
