Index: trunk/flexpart90.info
===================================================================
--- trunk/flexpart90.info	(revision 20)
+++ 	(revision )
@@ -1,1 +1,0 @@
-...is not flexpart77.
Index: trunk/options/COMMAND
===================================================================
--- trunk/options/COMMAND	(revision 20)
+++ trunk/options/COMMAND	(revision 24)
@@ -1,31 +1,198 @@
-+++++++++++++ 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
+********************************************************************************
+*                                                                              *
+*      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
Index: trunk/options/COMMAND.alternative
===================================================================
--- trunk/options/COMMAND.alternative	(revision 20)
+++ trunk/options/COMMAND.alternative	(revision 24)
@@ -30,4 +30,5 @@
 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
 
 
@@ -117,2 +118,6 @@
     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 20)
+++ trunk/options/COMMAND.reference	(revision 24)
@@ -11,9 +11,9 @@
 
 2. ________ ______   3X, I8, 1X, I6
-   20040720 000000
+   20110310 000000
    YYYYMMDD HHMISS   BEGINNING DATE OF SIMULATION
 
 3. ________ ______   3X, I8, 1X, I6
-   20040721 120000
+   20110310 120000
    YYYYMMDD HHMISS   ENDING DATE OF SIMULATION
 
@@ -101,4 +101,8 @@
     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
 
 
@@ -189,2 +193,6 @@
     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 20)
+++ trunk/options/RELEASES	(revision 24)
@@ -1,868 +1,54 @@
-+++++++++++++++++ HEADER ++++++++++++++++++++
-+++++++++++++++++ HEADER ++++++++++++++++++++
-+++++++++++++++++ HEADER ++++++++++++++++++++
-+++++++++++++++++ HEADER ++++++++++++++++++++
-+++++++++++++++++ HEADER ++++++++++++++++++++
-+++++++++++++++++ HEADER ++++++++++++++++++++
-+++++++++++++++++ HEADER ++++++++++++++++++++
-+++++++++++++++++ HEADER ++++++++++++++++++++
-+++++++++++++++++ HEADER ++++++++++++++++++++
-+++++++++++++++++ HEADER ++++++++++++++++++++
-+++++++++++++++++ HEADER ++++++++++++++++++++
-5
-31
-31
-31
-31
-31
+*************************************************************************
+*                                                                       *
+*                                                                       *
+*                                                                       *
+*   Input file for the Lagrangian particle dispersion model FLEXPART    *
+*                        Please select your options                     *
+*                                                                       *
+*                                                                       *
+*                                                                       *
+*************************************************************************
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  1      
+___                        i3    Total number of species emitted
 
-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
+  3
+___                        i3    Index of species in file SPECIES
 
-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
+=========================================================================
+20111210 000001 
+________ ______            i8,1x,i6 Beginning date and time of release
 
-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
+20111210 090000 
+________ ______            i8,1x,i6 Ending date and time of release
 
-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
+   9.4048 
+____.____                  f9.4  Longitude [DEG] of lower left corner
 
-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
+  48.5060
+____.____                  f9.4  Latitude [DEG] of lower left corner
 
-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
+   9.5067
+____.____                  f9.4  Longitude [DEG] of upper right corner
 
-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
+  48.5158
+____.____                  f9.4  Latitude [DEG] of upper right corner
 
-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
+        2
+_________                  i9    1 for m above ground, 2 for m above sea level
 
-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
+    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
 
-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
+1.0000E00
+_.____E__                  e9.4  Total mass emitted
 
-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
-
+RELEASE_TEST1
+________________________________________   character*40 comment
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Index: trunk/options/SPECIES/SPECIES_001
===================================================================
--- trunk/options/SPECIES/SPECIES_001	(revision 20)
+++ trunk/options/SPECIES/SPECIES_001	(revision 24)
@@ -9,4 +9,8 @@
 -9.9E-09           Wet deposition - A
                    Wet deposition - 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.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: trunk/options/SPECIES/SPECIES_002
===================================================================
--- trunk/options/SPECIES/SPECIES_002	(revision 20)
+++ trunk/options/SPECIES/SPECIES_002	(revision 24)
@@ -9,4 +9,8 @@
 -9.9E-09           Wet deposition - A
                    Wet deposition - 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)
  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 20)
+++ trunk/options/SPECIES/SPECIES_003	(revision 24)
@@ -9,4 +9,8 @@
  8.0E-06           Wet deposition - A
 0.62               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.1                In-cloud scavenging - Ci (S_i=Ci/cl)
+1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
  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 20)
+++ trunk/options/SPECIES/SPECIES_004	(revision 24)
@@ -9,4 +9,8 @@
  1.0E-05           Wet deposition - A
 0.62               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.1                In-cloud scavenging - Ci (S_i=Ci/cl)
+1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
  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 20)
+++ trunk/options/SPECIES/SPECIES_005	(revision 24)
@@ -9,4 +9,8 @@
  5.0E-05           Wet deposition - A
 0.62               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.1                In-cloud scavenging - Ci (S_i=Ci/cl)
+1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
  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 20)
+++ trunk/options/SPECIES/SPECIES_006	(revision 24)
@@ -9,4 +9,8 @@
 -9.9E-09           Wet deposition - A
                    Wet deposition - 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)
  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 20)
+++ trunk/options/SPECIES/SPECIES_007	(revision 24)
@@ -9,4 +9,8 @@
  1.0E-04           Wet deposition - A
 0.62               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.1                In-cloud scavenging - Ci (S_i=Ci/cl)
+1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
  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 20)
+++ trunk/options/SPECIES/SPECIES_008	(revision 24)
@@ -9,4 +9,8 @@
 -9.9E-09           Wet deposition - A
 0.62               Wet deposition - 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)
  2.0               Dry deposition (gases) - D
 1.0e+05            Dry deposition (gases) - Henrys const.
@@ -39,5 +43,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 20)
+++ trunk/options/SPECIES/SPECIES_009	(revision 24)
@@ -9,4 +9,8 @@
 -9.9E-09           Wet deposition - A
                    Wet deposition - 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)
  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 20)
+++ trunk/options/SPECIES/SPECIES_010	(revision 24)
@@ -9,4 +9,8 @@
 -9.9E-09           Wet deposition - A
                    Wet deposition - 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)
  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 20)
+++ trunk/options/SPECIES/SPECIES_011	(revision 24)
@@ -9,4 +9,8 @@
  9.9E-05           Wet deposition - A
 0.62               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.1                In-cloud scavenging - Ci (S_i=Ci/cl)
+1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
  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 20)
+++ trunk/options/SPECIES/SPECIES_012	(revision 24)
@@ -9,4 +9,8 @@
  5.0E-06           Wet deposition - A
 0.62               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.1                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
  1.0E-09           Dry deposition (gases) - Henrys const.
Index: trunk/options/SPECIES/SPECIES_013
===================================================================
--- trunk/options/SPECIES/SPECIES_013	(revision 20)
+++ trunk/options/SPECIES/SPECIES_013	(revision 24)
@@ -9,4 +9,8 @@
  5.0E-06           Wet deposition - A
 0.62               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.1                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/options/SPECIES/SPECIES_014
===================================================================
--- trunk/options/SPECIES/SPECIES_014	(revision 20)
+++ trunk/options/SPECIES/SPECIES_014	(revision 24)
@@ -9,4 +9,8 @@
  8.0E-05           Wet deposition - A
 0.62               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.1                In-cloud scavenging - Ci (S_i=Ci/cl)
+1.                 In-cloud scavenging - Di (wetscav=S_i*prec/3.6E6/clouds_h/Di)
  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 20)
+++ trunk/options/SPECIES/SPECIES_015	(revision 24)
@@ -9,4 +9,8 @@
  1.0E-04           Wet deposition - A
 0.80               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.1                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/options/SPECIES/SPECIES_016
===================================================================
--- trunk/options/SPECIES/SPECIES_016	(revision 20)
+++ trunk/options/SPECIES/SPECIES_016	(revision 24)
@@ -9,4 +9,8 @@
  1.0E-04           Wet deposition - A
 0.80               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.1                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/options/SPECIES/SPECIES_017
===================================================================
--- trunk/options/SPECIES/SPECIES_017	(revision 20)
+++ trunk/options/SPECIES/SPECIES_017	(revision 24)
@@ -9,4 +9,8 @@
  1.0E-04           Wet deposition - A
 0.80               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.1                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/options/SPECIES/SPECIES_018
===================================================================
--- trunk/options/SPECIES/SPECIES_018	(revision 20)
+++ trunk/options/SPECIES/SPECIES_018	(revision 24)
@@ -9,4 +9,8 @@
  1.0E-04           Wet deposition - A
 0.80               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.1                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/options/SPECIES/SPECIES_019
===================================================================
--- trunk/options/SPECIES/SPECIES_019	(revision 20)
+++ trunk/options/SPECIES/SPECIES_019	(revision 24)
@@ -9,4 +9,8 @@
 -9.9E-09           Wet deposition - A
                    Wet deposition - 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.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: trunk/options/SPECIES/SPECIES_020
===================================================================
--- trunk/options/SPECIES/SPECIES_020	(revision 20)
+++ trunk/options/SPECIES/SPECIES_020	(revision 24)
@@ -9,4 +9,8 @@
  1.0E-04           Wet deposition - A
 0.80               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.1                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/options/SPECIES/SPECIES_021
===================================================================
--- trunk/options/SPECIES/SPECIES_021	(revision 20)
+++ trunk/options/SPECIES/SPECIES_021	(revision 24)
@@ -6,7 +6,11 @@
 ****************************************************************************
 Xe-133             Tracer name
-198720.0           Species half life
+453168.0           Species half life
 -9.9E-09           Wet deposition - A
                    Wet deposition - 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.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: trunk/options/SPECIES/SPECIES_022
===================================================================
--- trunk/options/SPECIES/SPECIES_022	(revision 20)
+++ trunk/options/SPECIES/SPECIES_022	(revision 24)
@@ -9,4 +9,8 @@
 -9.9E-09           Wet deposition - A
                    Wet deposition - 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.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: trunk/options/SPECIES/SPECIES_024
===================================================================
--- trunk/options/SPECIES/SPECIES_024	(revision 20)
+++ trunk/options/SPECIES/SPECIES_024	(revision 24)
@@ -9,4 +9,8 @@
 -9.9E-09           Wet deposition - A
                    Wet deposition - 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.9               Dry deposition (gases) - D
                    Dry deposition (gases) - Henrys const.
Index: trunk/options/SPECIES/SPECIES_025
===================================================================
--- trunk/options/SPECIES/SPECIES_025	(revision 20)
+++ trunk/options/SPECIES/SPECIES_025	(revision 24)
@@ -9,4 +9,8 @@
 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.1                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 20)
+++ trunk/src/FLEXPART.f90	(revision 24)
@@ -61,5 +61,5 @@
 
   ! 
-  flexversion='Version 9.1.8  (2013-12-08)'
+  flexversion='Version 9.2 beta (2014-05-23)'
   !verbosity=0
   ! Read the pathnames where input/output files are stored
@@ -111,5 +111,5 @@
   ! Print the GPL License statement
   !*******************************************************
-  print*,'Welcome to FLEXPART', trim(flexversion)
+  print*,'Welcome to FLEXPART ', trim(flexversion)
   print*,'FLEXPART is free software released under the GNU Genera'// &
        'l Public License.'
Index: trunk/src/calcpv.f90
===================================================================
--- trunk/src/calcpv.f90	(revision 20)
+++ trunk/src/calcpv.f90	(revision 24)
@@ -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,ppmk
-  real :: pvavr,ppml(nuvzmax)
+  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 :: thup,thdn
   real,parameter :: eps=1.e-5, p0=101325
@@ -62,4 +62,13 @@
   ! 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
@@ -106,8 +115,4 @@
       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
@@ -115,6 +120,5 @@
 
       do kl=1,nuvz
-        ppmk=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
-        theta=tth(ix,jy,kl,n)*(100000./ppmk)**kappa
+        theta=tth(ix,jy,kl,n)*ppmk(ix,jy,kl)
         klvrp=kl+1
         klvrm=kl-1
@@ -125,9 +129,7 @@
         if (klvrp.gt.nuvz) klvrp=nuvz
         if (klvrm.lt.1) klvrm=1
-        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))
+        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))
 
   ! Compute vertical position at pot. temperature surface on subgrid
@@ -156,23 +158,21 @@
           kch=kch+1
           k=kup
-          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
+          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
             endif
+            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
+            goto 20
+          endif
 41        continue
   ! Downward branch
@@ -181,23 +181,21 @@
           kch=kch+1
           k=kdn
-          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
+          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
             endif
-            goto 40
+            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
+            goto 20
+          endif
+          goto 40
   ! This section used when no values were found
 21      continue
@@ -236,21 +234,19 @@
           kch=kch+1
           k=kup
-          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
+          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
             endif
+            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
+            goto 50
+          endif
 71        continue
   ! Downward branch
@@ -259,27 +255,25 @@
           kch=kch+1
           k=kdn
-          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
+          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
             endif
-            goto 70
+            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
+            goto 50
+          endif
+          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 20)
+++ trunk/src/calcpv_nests.f90	(revision 24)
@@ -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,ppmk
-  real :: ppml(nuvzmax)
+  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 :: thup,thdn
   real,parameter :: eps=1.e-5,p0=101325
@@ -61,4 +61,13 @@
   ! 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.
@@ -88,8 +97,4 @@
       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
@@ -97,6 +102,5 @@
 
       do kl=1,nuvz
-        ppmk=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
-        theta=tthn(ix,jy,kl,n,l)*(100000./ppmk)**kappa
+        theta=tthn(ix,jy,kl,n,l)*ppmk(ix,jy,kl)
         klvrp=kl+1
         klvrm=kl-1
@@ -107,9 +111,7 @@
         if (klvrp.gt.nuvz) klvrp=nuvz
         if (klvrm.lt.1) klvrm=1
-        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))
+        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))
 
   ! Compute vertical position at pot. temperature surface on subgrid
@@ -134,8 +136,6 @@
           kch=kch+1
           k=kup
-          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
+          thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k)
+          thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1)
 
       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
@@ -158,8 +158,6 @@
           kch=kch+1
           k=kdn
-          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
+          thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k)
+          thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1)
       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
            ((thdn.le.theta).and.(thup.ge.theta))) then
@@ -212,8 +210,6 @@
           kch=kch+1
           k=kup
-          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
+          thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k)
+          thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1)
       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
            ((thdn.le.theta).and.(thup.ge.theta))) then
@@ -235,8 +231,6 @@
           kch=kch+1
           k=kdn
-          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
+          thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k)
+          thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1)
       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 20)
+++ trunk/src/com_mod.f90	(revision 24)
@@ -340,9 +340,6 @@
   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)
 
 
@@ -359,7 +356,4 @@
   !      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
Index: trunk/src/gridcheck.f90
===================================================================
--- trunk/src/gridcheck.f90	(revision 20)
+++ trunk/src/gridcheck.f90	(revision 24)
@@ -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
+  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parID
   !HSO  end
   integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
-  real :: sizesouth,sizenorth,xauxa,pint
+  real :: sizesouth,sizenorth,xauxa,pint,conversion_factor
 
   ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
@@ -98,6 +98,5 @@
 
   integer :: isec1(56),isec2(22+nxmax+nymax)
-  !real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
-  real(kind=4) :: zsec2(184),zsec4(jpunp)
+  real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
   character(len=1) :: opt
 
@@ -172,4 +171,6 @@
   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)
 
@@ -193,4 +194,6 @@
   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
@@ -205,7 +208,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)) then ! CC
+  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
     isec1(6)=164         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
+  elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
     isec1(6)=142         ! indicatorOfParameter
   elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
@@ -215,11 +218,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)) then ! EWSS
+  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
     isec1(6)=180         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
+  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) 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)) then ! SDO
+  elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
     isec1(6)=160         ! indicatorOfParameter
   elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
@@ -230,4 +233,8 @@
          parCat,parNum,typSurf
   endif
+  if(parId .ne. isec1(6) .and. parId .ne. 77) then
+    write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
+!    stop
+  endif
 
   endif
@@ -265,5 +272,5 @@
 
   !HSO  get the second part of the grid dimensions only from GRiB1 messages
-  if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
+  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
     call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
          xaux2in,iret)
@@ -291,5 +298,4 @@
     dyconst=180./(dy*r_earth*pi)
     gotGrid=1
-
   ! Check whether fields are global
   ! If they contain the poles, specify polar stereographic map
@@ -443,7 +449,7 @@
   write(*,*)
   write(*,'(a)') 'Mother domain:'
-  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
+  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  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.5,a,f10.5,a,f10.5)') '  Latitude range: ', &
        ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
   write(*,*)
@@ -467,4 +473,5 @@
     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
@@ -553,2 +560,3 @@
 
 end subroutine gridcheck
+
Index: trunk/src/gridcheck_nests.f90
===================================================================
--- trunk/src/gridcheck_nests.f90	(revision 20)
+++ trunk/src/gridcheck_nests.f90	(revision 24)
@@ -48,4 +48,5 @@
   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
@@ -55,4 +56,5 @@
   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
@@ -150,4 +152,6 @@
   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
@@ -170,4 +174,6 @@
   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
@@ -182,7 +188,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)) then ! CC
+  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
     isec1(6)=164         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
+ 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
     isec1(6)=142         ! indicatorOfParameter
   elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
@@ -192,11 +198,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)) then ! EWSS
+  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
     isec1(6)=180         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
+  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
     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)) then ! SDO
+  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
     isec1(6)=160         ! indicatorOfParameter
   elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
@@ -207,4 +213,8 @@
          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
@@ -237,6 +247,6 @@
 
   !HSO  get the second part of the grid dimensions only from GRiB1 messages
-  if ((gribVer.eq.1).and.(gotGrib.eq.0)) then
-    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
+  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
          xaux1in,iret)
     call grib_check(iret,gribFunction,gribErrorMsg)
@@ -254,14 +264,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.
+    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
     xlon0n(l)=xaux1
     ylat0n(l)=yaux1
     dxn(l)=(xaux2-xaux1)/real(nxn(l)-1)
     dyn(l)=(yaux2-yaux1)/real(nyn(l)-1)
-    gotGrib=1
+    gotGrib=1 !commetn by mc note tahthere gotGRIB is used instead of gotGrid!!!
   endif ! ifield.eq.1
 
@@ -341,8 +351,8 @@
 
   write(*,'(a,i2)') 'Nested domain #: ',l
-  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
+  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
        xlon0n(l),' to ',xlon0n(l)+(nxn(l)-1)*dxn(l), &
        '   Grid distance: ',dxn(l)
-  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
+  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  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 20)
+++ trunk/src/interpol_rain.f90	(revision 24)
@@ -20,15 +20,8 @@
 !**********************************************************************
 
-!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
-
+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
   !****************************************************************************
   !                                                                           *
@@ -47,8 +40,5 @@
   !                                                                           *
   !     30 August 1996                                                        *
-  !
-  !*  Petra Seibert, 2011/2012:
-  !*  Add interpolation of cloud bottom and cloud thickness
-  !*  for fix to SE's new wet scavenging scheme  *
+  !                                                                           *
   !****************************************************************************
   !                                                                           *
@@ -80,15 +70,10 @@
 
   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,ip1,ip2,ip3,ip4 
-  integer :: intiy1,intiy2,ipsum,icmv  
+  integer :: itime,itime1,itime2,level,indexh
   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)
-  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
+  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
+  real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
 
 
@@ -142,43 +127,4 @@
          + 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
 
@@ -187,11 +133,4 @@
   ! 2.) Temporal interpolation (linear)
   !************************************
-
-      if (abs(itime) .lt. abs(itime1)) then
-        print*,'interpol_rain.f90'
-        print*,itime,itime1,itime2
-        stop 'ITIME PROBLEM'
-      endif
-
 
   dt1=real(itime-itime1)
@@ -204,20 +143,3 @@
 
 
-!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 20)
+++ trunk/src/makefile	(revision 24)
@@ -1,29 +1,13 @@
 SHELL = /bin/bash
-TARGET = local
-WINDS=ecmwf
-#WINDS=gfs
-#WINDS=fnl
+MAIN = FP_ecmwf_gfortran
 
 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)
+LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
 
-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
 
 MODOBJS = \
@@ -37,5 +21,5 @@
 OBJECTS = \
 writeheader.o  writeheader_txt.o   writeheader_surf.o       assignland.o\
-               part0.o \
+calcpar.o               part0.o \
 caldate.o               partdep.o \
 coordtrafo.o            psih.o \
@@ -53,24 +37,26 @@
 getrc.o                 readreleases.o \
 getvdep.o               readspecies.o \
-interpol_misslev.o      \
-conccalc.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 \
-partoutput.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 \
-                           calcpar_nests.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   gridcheck_nests.o \
-readwind_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 \
@@ -91,22 +77,6 @@
 
 
-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)
+$(MAIN): $(MODOBJS) $(OBJECTS)
+	$(FC) *.o -o $(MAIN) $(LDFLAGS)
 
 $(OBJECTS): $(MODOBJS)
@@ -118,2 +88,4 @@
 	rm *.o *.mod
 
+cleanall:
+	rm *.o *.mod $(MAIN)
Index: trunk/src/par_mod.f90
===================================================================
--- trunk/src/par_mod.f90	(revision 20)
+++ trunk/src/par_mod.f90	(revision 24)
@@ -121,5 +121,6 @@
   !*********************************************
   
-  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 !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=26,nwzmax=26,nzmax=26
@@ -154,6 +155,6 @@
 
   !integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
-  !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151 !ECMWF
-  integer,parameter :: maxnests=1, nxmaxn=201, nymaxn=161 ! FNL XF
+  integer,parameter :: maxnests=1,nxmaxn=351,nymaxn=351 !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
@@ -197,7 +198,6 @@
   !**************************************************
 
-  !integer,parameter :: maxpart=4000000
-  integer,parameter :: maxpart=2000
-  integer,parameter :: maxspec=6
+  integer,parameter :: maxpart=15000000
+  integer,parameter :: maxspec=4
 
 
Index: trunk/src/readcommand.f90
===================================================================
--- trunk/src/readcommand.f90	(revision 20)
+++ trunk/src/readcommand.f90	(revision 24)
@@ -80,33 +80,33 @@
   character(len=50) :: line
   logical :: old
-  logical :: nmlout=.true. !.false.
+  logical :: nml_COMMAND=.true. , 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, &
-    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
@@ -144,19 +144,39 @@
   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)
+   ! 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
+
+   ! print error code
+   !write(*,*) 'readcommand > readerror open=' , readerror
+   !probe first line  
+   read (unitcommand,901) line
+   !write(*,*) 'index(line,COMMAND) =', index(line,'COMMAND') 
+
+ 
+   !default is namelist input 
+   ! distinguish namelist from fixed text input
+   if (index(line,'COMMAND') .eq. 0) then
+   nml_COMMAND = .false. 
+   !write(*,*) 'COMMAND file does not contain the string COMMAND in the first line'      
+   endif
+   !write(*,*) 'readcommand > read as namelist? ' , nml_COMMAND
+   rewind(unitcommand)
+   read(unitcommand,command,iostat=readerror)
+
   close(unitcommand)
 
+  !write(*,*) 'readcommand > readerror read=' , readerror
   ! If error in namelist format, try to open with old input code
-  if (readerror.ne.0) then
+  ! if (readerror.ne.0) then
+  ! IP 21/5/2014 the previous line cause the old long format 
+  ! to be confused with namelist input
+  
+  ! use text input 
+  if (nml_COMMAND .eqv. .false.) then 
 
     open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
@@ -173,9 +193,12 @@
     if (index(line,'LDIRECT') .eq. 0) then
       old = .false.
+    !write(*,*) 'readcommand old short'
     else
       old = .true.
+    !write(*,*) 'readcommand old long'
     endif
     rewind(unitcommand)
 
+
     ! Read parameters
     !****************
@@ -231,4 +254,6 @@
     if (old) call skplin(3,unitcommand)
     read(unitcommand,*) linit_cond
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) surf_only
     close(unitcommand)
 
@@ -237,5 +262,4 @@
   ! 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)
Index: trunk/src/readreleases.f90
===================================================================
--- trunk/src/readreleases.f90	(revision 20)
+++ trunk/src/readreleases.f90	(revision 24)
@@ -55,5 +55,7 @@
   ! xpoint2,ypoint2     geograf. coordinates of upper right corner of release  *
   !                     area                                                   *
-  ! weta, wetb          parameters to determine the wet scavenging coefficient *
+  ! 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        *
   ! zpoint1,zpoint2     height range, over which release takes place           *
   ! num_min_discrete    if less, release cannot be randomized and happens at   *
@@ -272,5 +274,5 @@
         vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
       end do
-      write(*,*) 'Average setting velocity: ',i,vsetaver(i)
+      write(*,*) 'Average settling velocity: ',i,vsetaver(i)
     endif
 
@@ -284,6 +286,19 @@
     if (weta(i).gt.0.)  then
       WETDEP=.true.
-      write (*,*) 'Wetdeposition switched on: ',weta(i),i
-    endif
+      write (*,*) 'Below-cloud scavenging is switched on'
+      write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
+    else
+      write (*,*) 'Below-cloud scavenging is switched OFF'
+    endif
+    
+! NIK 31.01.2013 + 10.12.2013
+    if (weta_in(i).gt.0.)  then
+      WETDEP=.true.
+      write (*,*) 'In-cloud scavenging is switched on'
+      write (*,*) 'In-cloud scavenging coefficients: ',weta_in(i),wetb_in(i), wetc_in(i), wetd_in(i),i
+    else
+      write (*,*) 'In-cloud scavenging is switched OFF' 
+    endif
+
     if (ohreact(i).gt.0) then
       OHREA=.true.
Index: trunk/src/readspecies.f90
===================================================================
--- trunk/src/readspecies.f90	(revision 20)
+++ trunk/src/readspecies.f90	(revision 24)
@@ -31,4 +31,7 @@
   !   11 July 1996                                                             *
   !                                                                            *
+  !   Changes:                                                                 *
+  !   N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging *
+  !                                                                            *
   !*****************************************************************************
   !                                                                            *
@@ -36,5 +39,9 @@
   ! decaytime(maxtable)  half time for radiological decay                      *
   ! specname(maxtable)   names of chemical species, radionuclides              *
-  ! wetscava, wetscavb   Parameters for determining scavenging coefficient     *
+  ! 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         *
   ! ohreact              OH reaction rate                                      *
   ! id_spec              SPECIES number as referenced in RELEASE file          *
@@ -78,4 +85,15 @@
     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)
Index: trunk/src/readwind.f90
===================================================================
--- trunk/src/readwind.f90	(revision 20)
+++ trunk/src/readwind.f90	(revision 24)
@@ -71,5 +71,5 @@
   integer :: iret
   integer :: igrib
-  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
+  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
   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,xaux0,yaux0
+  real(kind=4) :: xaux,yaux
   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
+  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
 
   logical :: hflswitch,strswitch
@@ -101,9 +101,4 @@
   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.
@@ -155,4 +150,6 @@
   endif
 
+  conversion_factor=1.
+
   else
 
@@ -168,4 +165,6 @@
   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)
 
@@ -177,4 +176,5 @@
   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,4 +188,6 @@
     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
@@ -201,21 +203,23 @@
   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)) then ! CC
+    conversion_factor=1000.
+  elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
     isec1(6)=164         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
+  elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) 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)) then ! EWSS
+  elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
     isec1(6)=180         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
+  elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) 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)) then ! SDO
+  elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
     isec1(6)=160         ! indicatorOfParameter
   elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
@@ -223,6 +227,10 @@
     isec1(6)=172         ! indicatorOfParameter
   else
-    print*,'***ERROR: undefined GRiB2 message found!',discipl, &
+    print*,'***WARNING: 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
 
@@ -237,12 +245,9 @@
   !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
@@ -250,9 +255,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 ((gribVer.eq.1).and.(gotGrid.eq.0)) then
+  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
     call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
          xauxin,iret)
@@ -261,16 +266,14 @@
          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
-    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'
+    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'
     gotGrid=1
   endif ! gotGrid
@@ -298,5 +301,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)
+           zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
       if(isec1(6).eq.151) msl(i,j,1,n)= &!! SEA LEVEL PRESS.
            zsec4(nxfield*(ny-j-1)+i+1)
@@ -316,5 +319,5 @@
       endif
       if(isec1(6).eq.143) then                      !! CONVECTIVE PREC.
-        convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
+        convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
         if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
       endif
@@ -366,20 +369,4 @@
       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
@@ -479,2 +466,3 @@
 
 end subroutine readwind
+
Index: trunk/src/readwind_nests.f90
===================================================================
--- trunk/src/readwind_nests.f90	(revision 20)
+++ trunk/src/readwind_nests.f90	(revision 24)
@@ -51,4 +51,5 @@
   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
@@ -69,8 +70,10 @@
   integer :: isec1(56),isec2(22+nxmaxn+nymaxn)
   real(kind=4) :: zsec4(jpunp)
+  real(kind=4) :: xaux,yaux
   real(kind=8) :: xauxin,yauxin
-  real(kind=4) :: xaux,yaux,xaux0,yaux0
+  real,parameter :: eps=1.e-4
   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
@@ -134,4 +137,7 @@
   endif
 
+  conversion_factor=1.
+
+
   else
 
@@ -148,4 +154,6 @@
   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
@@ -156,4 +164,5 @@
   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
@@ -166,6 +175,8 @@
   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
@@ -180,21 +191,23 @@
   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)) then ! CC
+    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
     isec1(6)=164         ! indicatorOfParameter
-  elseif ((parCat.eq.1).and.(parNum.eq.9)) then ! LSP
+  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
     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)) then ! EWSS
+  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
     isec1(6)=180         ! indicatorOfParameter
-  elseif ((parCat.eq.2).and.(parNum.eq.18)) then ! NSSS
+  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
     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)) then ! SDO
+   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
     isec1(6)=160         ! indicatorOfParameter
   elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
@@ -202,6 +215,10 @@
     isec1(6)=172         ! indicatorOfParameter
   else
-    print*,'***ERROR: undefined GRiB2 message found!',discipl, &
+    print*,'***WARNING: 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
 
@@ -227,13 +244,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 ((gribVer.eq.1).and.(gotGrid.eq.0)) then
+ if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consisitent with new readwind.f90
     call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
          xauxin,iret)
@@ -242,18 +259,13 @@
          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
-    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'
+    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'
     gotGrid=1
   endif
@@ -280,5 +292,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)
+             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90!
         if(isec1(6).eq.151) msln(i,j,1,n,l)= &!! SEA LEVEL PRESS.
              zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
@@ -298,5 +310,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)
+          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
           if (convprecn(i,j,1,n,l).lt.0.) convprecn(i,j,1,n,l)=0.
         endif
Index: trunk/src/verttransform.f90
===================================================================
--- trunk/src/verttransform.f90	(revision 20)
+++ trunk/src/verttransform.f90	(revision 24)
@@ -21,5 +21,5 @@
 
 subroutine verttransform(n,uuh,vvh,wwh,pvh)
-  !                         i  i   i   i   i
+  !                      i  i   i   i   i
   !*****************************************************************************
   !                                                                            *
@@ -49,6 +49,4 @@
   ! 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
   !*****************************************************************************
   !                                                                            *
@@ -72,13 +70,11 @@
 
   integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
- integer :: rain_cloud_above,kz_inv !SE
-  integer icloudtop !PS
+  integer :: rain_cloud_above,kz_inv
   real :: f_qvsat,pressure
- !real :: rh,lsp,convp
-  real :: rh,lsp,convp,prec,rhmin  
-  real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
+  real :: rh,lsp,convp
+  real :: 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, precmin
+  real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
   real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
   real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
@@ -87,7 +83,5 @@
   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.
@@ -97,15 +91,6 @@
   ! 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, 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.                                                  *
+  ! u,v,T and qv are given.                                                *
   !*************************************************************************
-
-
-  !      do 897 kz=1,nuvz
-  !       write (*,*) 'akz: ',akz(kz),'bkz',bkz(kz)
-  !897     continue
 
   if (init) then
@@ -113,5 +98,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
@@ -128,5 +113,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.
@@ -136,46 +121,14 @@
       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
 
 
@@ -204,8 +157,7 @@
   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)
-      uvzlev(1)=0.
+      uvwzlev(ix,jy,1)=0.
       wzlev(1)=0.
       rhoh(1)=pold/(r_air*tvold)
@@ -221,8 +173,8 @@
 
         if (abs(tv-tvold).gt.0.2) then
-          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
-               (tv-tvold)/log(tv/tvold)
+          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
+          (tv-tvold)/log(tv/tvold)
         else
-          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
+          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
         endif
 
@@ -233,35 +185,21 @@
 
       do kz=2,nwz-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
+        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)
 
   ! 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
@@ -283,5 +221,5 @@
       do iz=2,nz-1
         do kz=kmin,nuvz
-          if(height(iz).gt.uvzlev(nuvz)) then
+          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
             uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
             vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
@@ -292,8 +230,8 @@
             goto 30
           endif
-          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)
+          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
@@ -339,8 +277,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)
@@ -356,4 +294,5 @@
 
   do jy=1,ny-2
+    cosf=cos((real(jy)*dy+ylat0)*pi180)
     do ix=1,nx-2
 
@@ -361,10 +300,10 @@
       do iz=2,nz-1
 
-        ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
+        ui=uu(ix,jy,iz,n)*dxconst/cosf
         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)
@@ -409,6 +348,5 @@
         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
@@ -424,12 +362,10 @@
       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
@@ -444,7 +380,5 @@
       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
@@ -485,6 +419,5 @@
         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
@@ -499,12 +432,9 @@
       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
@@ -519,6 +449,5 @@
       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
@@ -551,135 +480,40 @@
   !   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
-  !          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
+      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
-          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
+         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
-          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
-
-
+         endif
+      end do
     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 20)
+++ trunk/src/verttransform_gfs.f90	(revision 24)
@@ -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 :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
+  real :: 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
+  real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
   real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
   real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
@@ -92,9 +92,5 @@
   ! 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, 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.                                                  *
+  ! u,v,T and qv are given.                                                *
   !*************************************************************************
 
@@ -118,5 +114,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.
@@ -126,46 +122,14 @@
       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
 
 
@@ -198,5 +162,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
@@ -212,5 +176,4 @@
       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.
@@ -223,24 +186,14 @@
 
         if (abs(tv-tvold).gt.0.2) then
-          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
-               (tv-tvold)/log(tv/tvold)
+          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
+          (tv-tvold)/log(tv/tvold)
         else
-          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
+          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
         endif
-        wzlev(kz)=uvzlev(kz)
-        uvwzlev(ix,jy,kz)=uvzlev(kz)
+        wzlev(kz)=uvwzlev(ix,jy,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)
@@ -279,5 +232,5 @@
       do iz=2,nz-1
         do kz=kmin,nuvz
-          if(height(iz).gt.uvzlev(nuvz)) then
+          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
             uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
             vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
@@ -289,18 +242,18 @@
             goto 30
           endif
-          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
+          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
           endif
         end do
@@ -318,11 +271,10 @@
         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
@@ -337,5 +289,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)
@@ -351,4 +303,5 @@
 
   do jy=1,ny-2
+    cosf=cos((real(jy)*dy+ylat0)*pi180)
     do ix=1,nx-2
 
@@ -367,10 +320,10 @@
       do iz=2,nz-1
 
-        ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
+        ui=uu(ix,jy,iz,n)*dxconst/cosf
         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)
@@ -426,12 +379,10 @@
       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
@@ -446,7 +397,5 @@
       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
@@ -460,5 +409,5 @@
   ! ward parallel of latitude
 
-  do iz=1,nz
+    do iz=1,nz
       wdummy=0.
       jy=ny-2
@@ -471,5 +420,5 @@
         ww(ix,jy,iz,n)=wdummy
       end do
-  end do
+    end do
 
   endif
@@ -487,6 +436,5 @@
         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
@@ -498,12 +446,9 @@
       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
@@ -518,6 +463,5 @@
       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
@@ -562,24 +506,23 @@
          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 20)
+++ trunk/src/verttransform_nests.f90	(revision 24)
@@ -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 :: uvzlev(nuvzmax),wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax)
+  real :: 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
+  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi,cosf
   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)
-      uvzlev(1)=0.
+      uvwzlev(ix,jy,1)=0.
       wzlev(1)=0.
       rhoh(1)=pold/(r_air*tvold)
@@ -112,8 +112,8 @@
 
         if (abs(tv-tvold).gt.0.2) then
-          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
-               (tv-tvold)/log(tv/tvold)
+          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
+          (tv-tvold)/log(tv/tvold)
         else
-          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
+          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
         endif
 
@@ -124,43 +124,21 @@
 
       do kz=2,nwz-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
+        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)
 
   ! 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)))
 
 
@@ -183,5 +161,5 @@
       do iz=2,nz-1
         do kz=kmin,nuvz
-          if(height(iz).gt.uvzlev(nuvz)) then
+          if(height(iz).gt.uvwzlev(ix,jy,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)
@@ -192,19 +170,19 @@
             goto 30
           endif
-          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)
+          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
            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
@@ -225,12 +203,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
@@ -242,8 +220,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)
@@ -259,4 +237,5 @@
 
   do jy=1,nyn(l)-2
+    cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
     do ix=1,nxn(l)-2
 
@@ -264,6 +243,5 @@
       do iz=2,nz-1
 
-        ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/ &
-             cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
+        ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf
         vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)
 
@@ -319,5 +297,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 20)
+++ trunk/src/wetdepo.f90	(revision 24)
@@ -21,5 +21,5 @@
 
 subroutine wetdepo(itime,ltsample,loutnext)
-  !                     i      i        i
+  !                  i      i        i
   !*****************************************************************************
   !                                                                            *
@@ -39,7 +39,4 @@
   ! 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
   !*****************************************************************************
   !                                                                            *
@@ -74,11 +71,9 @@
 
   integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
-  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
+  integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
+  integer :: ks, kp
   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, precsub,f
+  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
   real :: wetdeposit(maxspec),restmass
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
@@ -103,4 +98,5 @@
 
   do jpart=1,numpart
+
     if (itra1(jpart).eq.-999999999) goto 20
     if(ldirect.eq.1)then
@@ -114,5 +110,5 @@
       if (itage.lt.lage(nage)) goto 33
     end do
-33   continue
+33  continue
 
 
@@ -123,10 +119,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
 
 
@@ -151,77 +147,42 @@
     interp_time=nint(itime-0.5*ltsample)
 
-!    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
-
-
+    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
   ! get the level were the actual particle is in
-      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)
+    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)
 
   ! if there is no precipitation or the particle is above the clouds no
   ! scavenging is done
 
-!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
-
+    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'
 
   ! 1) Parameterization of the the area fraction of the grid cell where the
@@ -267,28 +228,30 @@
   !**********************************************************
 
-    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
+    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
 
 ! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening 
@@ -297,142 +260,70 @@
 ! 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
+        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
-
-
-  !      if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
-  !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
+        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
         wetdeposit(ks)=xmass1(jpart,ks)* &
-             (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
+        (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
   ! 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
-
-   !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
-
+    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
 
 20  continue
-  end do
+  end do ! all particles
 
 end subroutine wetdepo
Index: trunk/src/writeheader_nest.f90
===================================================================
--- trunk/src/writeheader_nest.f90	(revision 20)
+++ trunk/src/writeheader_nest.f90	(revision 24)
@@ -67,7 +67,7 @@
 
   if (ldirect.eq.1) then
-    write(unitheader) ibdate,ibtime,'FLEXPART V8.2'
+    write(unitheader) ibdate,ibtime, trim(flexversion)
   else
-    write(unitheader) iedate,ietime,'FLEXPART V8.2'
+    write(unitheader) iedate,ietime, trim(flexversion)
   endif
 
