Index: palm/trunk/INSTALL/ex_lsm_clearsky_rc
===================================================================
--- palm/trunk/INSTALL/ex_lsm_clearsky_rc	(revision 2607)
+++ palm/trunk/INSTALL/ex_lsm_clearsky_rc	(revision 2608)
@@ -1,10 +1,10 @@
 
  ******************************    --------------------------------------------
- * PALM 4.0  Rev: 2532        *    atmosphere - run without 1D - prerun
+ * PALM 4.0  Rev: 2603:26     *    atmosphere - run without 1D - prerun
  ******************************    --------------------------------------------
 
- Date:                 18-10-17    Run:       ex_lsm_clearsky     
- Time:                 17:42:56    Run-No.:   00
- Run on host:              imuk
+ Date:                 13-11-17    Run:       ex_lsm_clearsky     
+ Time:                 14:31:43    Run-No.:   00
+ Run on host:             lcmuk
  Number of PEs:               4    Processor grid (x,y): (   2,   2) calculated
  ------------------------------------------------------------------------------
@@ -88,7 +88,4 @@
 
     Albedo is set for land surface type: mixed farming, tall grassland
-    Geograph. longitude            :   lambda =  0.0 degr
-    Day of the year at model start :   day_init = 172
-    UTC time at model start        :   time_utc_init = 43200.0 s
     Timestep: dt_radiation =  60.00  s
 
@@ -237,8 +234,11 @@
  -------------------
 
-    Geograph. latitude  :   phi    = 55.0 degr
+    Geograph. latitude  :   latitude  = 55.0 degr
+    Geograph. longitude :   longitude =  0.0 degr
     Angular velocity    :   omega  = 0.729E-04 rad/s
     Coriolis parameter  :   f      =  0.000119 1/s
                             f*     =  0.000084 1/s
+    Day of the year at model start :   day_init = 172
+    UTC time at model start        :   time_utc_init = 43200.0 s
 
     Gravity             :   g      =  9.8 m/s**2
@@ -284,15 +284,15 @@
   0      0 00:00:00.00   8.6700A   5.1898D -0.2288D  0.1182  0.100  0.00   0.200E-11    800.  0.119E+02  0.103E-02  0.611E-03  0.888E-13    9   9  23   11   7  19    5  26   7    0.000    0.000       0
   0      1 00:00:08.67   8.6900A   5.1806   0.2159   0.1162  0.001  0.00   0.000E+00    550.  0.119E+02  0.867E-03  0.287E-05  0.440E-15    9   9  24    4   5  26    5  26   8    0.000    0.000       0
-  0      2 00:00:17.35   8.6900A   5.1762   0.1915   0.1151  0.021  1.11  -0.474E+02    100.  0.119E+02  0.793E-03  0.759E-05  0.407E-15    5  16  32    4   5  27    5  26   9    0.000    0.000       0
-  0      3 00:00:26.04   8.7000A   5.1723   0.1786   0.1113  0.021  0.72  -0.545E+01    100.  0.119E+02  0.751E-03  0.253E-05  0.386E-15    5  16  33    4   5  28    5  26  10    0.000    0.000       0
-  0      4 00:00:34.75   8.7000A   5.1732  -0.1675   0.1061  0.028  0.70  -0.372E+01    100.  0.119E+02  0.721E-03  0.408E-05  0.358E-15    5  16  34    8  29  32    5  26  11    0.000    0.000       0
-  0      5 00:00:43.45   8.7000A   5.1719  -0.1671   0.1028  0.034  0.94  -0.297E+01    250.  0.119E+02  0.698E-03  0.487E-05  0.337E-15    5  16  35    8  29  33    5  26  11    0.000    0.000       0
-  0      6 00:00:52.15   8.7000A   5.1701  -0.1621   0.1039  0.040  0.94  -0.250E+01    250.  0.119E+02  0.680E-03  0.567E-05  0.321E-15    5  16  36    8  29  34    5  26  12    0.000    0.000       0
-  0      7 00:01:00.85   8.7100A   5.1675  -0.1565   0.1043  0.046  1.22  -0.218E+01    550.  0.119E+02  0.665E-03  0.649E-05  0.306E-15    5  16  37    8  29  35    5  26  13    0.000    0.000       0
-  0      8 00:01:09.56   8.7100A   5.1643   0.1515   0.1040  0.051  1.23  -0.195E+01    550.  0.119E+02  0.652E-03  0.733E-05  0.295E-15    5  16  38   11   3  22    5  26  14    0.000    0.000       0
-  0      9 00:01:18.27   8.7100A   5.1636  -0.1506   0.1034  0.057  1.21  -0.172E+01    550.  0.119E+02  0.640E-03  0.812E-05  0.282E-15    5  16  38    8  29  36    5  26  15    0.000    0.000       0
-  0     10 00:01:26.98   8.7100A   5.1636  -0.1516   0.1024  0.062  1.22  -0.158E+01    550.  0.119E+02  0.630E-03  0.898E-05  0.270E-15    5  16  39    8  29  37    5  26  16    0.000    0.000       0
-  0     11 00:01:35.69   8.7200A   5.1632  -0.1516   0.1013  0.067  1.22  -0.148E+01    550.  0.119E+02  0.620E-03  0.980E-05  0.260E-15    5  16   0    8  29  38    5  26  17    0.000    0.000       0
-  0     12 00:01:44.41   8.7200A   5.1625  -0.1506   0.0999  0.072  1.23  -0.139E+01    550.  0.119E+02  0.612E-03  0.106E-04  0.250E-15    5  16   1    8  29  39    5  26  18    0.000    0.000       0
-  0     13 00:01:53.13   8.7200A   5.1614  -0.1488   0.0984  0.077  1.23  -0.132E+01    550.  0.119E+02  0.604E-03  0.115E-04  0.239E-15    5  16   2    8  29   0    5  26  19    0.000    0.000       0
-  0     14 00:02:01.85   8.7200A   5.1600  -0.1462   0.0966  0.082  1.24  -0.126E+01    550.  0.119E+02  0.596E-03  0.123E-04  0.231E-15    5  16   3    8  29   1    5  26  20    0.000    0.000       0
+  0      2 00:00:17.35   8.6900A   5.1762   0.1915   0.1151  0.021  1.06  -0.495E+02    100.  0.119E+02  0.793E-03  0.744E-05  0.407E-15    5  16  32    4   5  27    5  26   9    0.000    0.000       0
+  0      3 00:00:26.04   8.7000A   5.1723   0.1786   0.1113  0.021  0.70  -0.516E+01    100.  0.119E+02  0.751E-03  0.247E-05  0.385E-15    5  16  33    4   5  28    5  26  10    0.000    0.000       0
+  0      4 00:00:34.75   8.7000A   5.1732  -0.1675   0.1061  0.027  0.69  -0.360E+01    100.  0.119E+02  0.721E-03  0.396E-05  0.358E-15    5  16  34    8  29  32    5  26  11    0.000    0.000       0
+  0      5 00:00:43.45   8.7000A   5.1719  -0.1671   0.1028  0.034  0.93  -0.287E+01    250.  0.119E+02  0.698E-03  0.471E-05  0.337E-15    5  16  35    8  29  33    5  26  11    0.000    0.000       0
+  0      6 00:00:52.15   8.7000A   5.1701  -0.1621   0.1039  0.039  0.93  -0.242E+01    250.  0.119E+02  0.680E-03  0.548E-05  0.321E-15    5  16  36    8  29  34    5  26  12    0.000    0.000       0
+  0      7 00:01:00.85   8.7100A   5.1675  -0.1565   0.1043  0.045  1.21  -0.211E+01    550.  0.119E+02  0.665E-03  0.627E-05  0.306E-15    5  16  37    8  29  35    5  26  13    0.000    0.000       0
+  0      8 00:01:09.56   8.7100A   5.1643   0.1515   0.1040  0.051  1.21  -0.189E+01    550.  0.119E+02  0.652E-03  0.709E-05  0.294E-15    5  16  38   11   3  22    5  26  14    0.000    0.000       0
+  0      9 00:01:18.27   8.7100A   5.1636  -0.1506   0.1034  0.056  1.20  -0.166E+01    550.  0.119E+02  0.640E-03  0.785E-05  0.282E-15    5  16  38    8  29  36    5  26  15    0.000    0.000       0
+  0     10 00:01:26.98   8.7100A   5.1635  -0.1516   0.1024  0.061  1.20  -0.153E+01    550.  0.119E+02  0.630E-03  0.868E-05  0.269E-15    5  16  39    8  29  37    5  26  16    0.000    0.000       0
+  0     11 00:01:35.69   8.7200A   5.1632  -0.1516   0.1013  0.067  1.20  -0.143E+01    550.  0.119E+02  0.620E-03  0.949E-05  0.260E-15    5  16   0    8  29  38    5  26  17    0.000    0.000       0
+  0     12 00:01:44.41   8.7200A   5.1624  -0.1506   0.0999  0.072  1.21  -0.135E+01    550.  0.119E+02  0.612E-03  0.103E-04  0.249E-15    5  16   1    8  29  39    5  26  18    0.000    0.000       0
+  0     13 00:01:53.13   8.7200A   5.1614  -0.1488   0.0983  0.077  1.22  -0.128E+01    550.  0.119E+02  0.604E-03  0.111E-04  0.239E-15    5  16   2    8  29   0    5  26  19    0.000    0.000       0
+  0     14 00:02:01.85   8.7200A   5.1600  -0.1462   0.0966  0.081  1.22  -0.122E+01    550.  0.119E+02  0.596E-03  0.119E-04  0.230E-15    5  16   3    8  29   1    5  26  20    0.000    0.000       0
Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 2607)
+++ palm/trunk/SOURCE/Makefile	(revision 2608)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# Added diagnostic_quantities_mod 
+# 
+# 2600 2017-11-01 14:11:20Z raasch
 # comment line concerning bound checks removed
 # 
@@ -384,5 +387,5 @@
 	data_output_ptseries.f90 data_output_spectra.f90 data_output_flight.f90\
 	data_output_tseries.f90 data_output_2d.f90 data_output_3d.f90 \
-	date_and_time_mod.f90 \
+	date_and_time_mod.f90  diagnostic_quantities_mod.f90 \
 	diffusion_e.f90 diffusion_s.f90 diffusion_u.f90 diffusion_v.f90 \
 	diffusion_w.f90 diffusivities.f90 disturb_field.f90 \
@@ -482,5 +485,5 @@
 buoyancy.o: modules.o mod_kinds.o
 calc_mean_profile.o: modules.o mod_kinds.o
-calc_liquid_water_content.o: modules.o mod_kinds.o
+calc_liquid_water_content.o: modules.o mod_kinds.o diagnostic_quantities_mod.o
 calc_radiation.o: modules.o mod_kinds.o
 check_for_restart.o: modules.o mod_kinds.o pmc_interface_mod.o
@@ -513,4 +516,5 @@
    netcdf_interface_mod.o land_surface_model_mod.o urban_surface_mod.o
 date_and_time_mod.o: modules.o mod_kinds.o 
+diagnostic_quantities_mod.o: modules.o mod_kinds.o
 diffusion_e.o: modules.o mod_kinds.o microphysics_mod.o \
    mod_particle_attributes.o surface_mod.o
@@ -552,5 +556,5 @@
 init_slope.o: modules.o mod_kinds.o
 interaction_droplets_ptq.o: modules.o mod_kinds.o
-land_surface_model_mod.o: modules.o mod_kinds.o radiation_model_mod.o surface_mod.o
+land_surface_model_mod.o: modules.o mod_kinds.o radiation_model_mod.o surface_mod.o diagnostic_quantities_mod.o
 large_scale_forcing_nudging_mod.o: modules.o cpulog_mod.o mod_kinds.o surface_mod.o
 local_stop.o: modules.o mod_kinds.o pmc_interface_mod.o
@@ -570,5 +574,5 @@
 	mod_kinds.o mod_particle_attributes.o
 lpm_droplet_condensation.o: modules.o cpulog_mod.o lpm_collision_kernels.o \
-	mod_kinds.o mod_particle_attributes.o
+	mod_kinds.o mod_particle_attributes.o diagnostic_quantities_mod.o
 lpm_exchange_horiz.o: modules.o cpulog_mod.o lpm_pack_arrays.o mod_kinds.o \
 	mod_particle_attributes.o netcdf_interface_mod.o
@@ -576,5 +580,5 @@
 	netcdf_interface_mod.o random_function_mod.o mod_particle_attributes.o \
 	lpm_exchange_horiz.o lpm_pack_arrays.o random_generator_parallel_mod.o \
-        surface_mod.o
+        surface_mod.o diagnostic_quantities_mod.o
 lpm_init_sgs_tke.o: modules.o mod_kinds.o mod_particle_attributes.o surface_mod.o
 lpm_merging.o: modules.o cpulog_mod.o mod_kinds.o mod_particle_attributes.o
@@ -588,5 +592,5 @@
 lpm_write_restart_file.o: modules.o mod_kinds.o mod_particle_attributes.o
 message.o: modules.o mod_kinds.o pmc_interface_mod.o
-microphysics_mod.o: modules.o cpulog_mod.o mod_kinds.o surface_mod.o
+microphysics_mod.o: modules.o cpulog_mod.o mod_kinds.o surface_mod.o diagnostic_quantities_mod.o
 model_1d_mod.o: modules.o mod_kinds.o
 modules.o: modules.f90 mod_kinds.o
Index: palm/trunk/SOURCE/calc_liquid_water_content.f90
===================================================================
--- palm/trunk/SOURCE/calc_liquid_water_content.f90	(revision 2607)
+++ palm/trunk/SOURCE/calc_liquid_water_content.f90	(revision 2608)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Calculation of supersaturation in external module (diagnostic_quantities_mod).
+! Change: correct calculation of saturation specific humidity to saturation
+! mixing ratio (the factor of 0.378 vanishes).
+! 
+! 2292 2017-06-20 09:51:42Z schwenkel
 ! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
 ! includes two more prognostic equations for cloud drop concentration (nc)  
@@ -82,11 +87,11 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  hyp, pt, q, qc, ql, qr
-
-    USE cloud_parameters,                                                      &
-        ONLY:  l_d_cp, l_d_r, t_d_pt
+        ONLY:  q, qc, ql, qr
 
     USE control_parameters,                                                    &
         ONLY:  microphysics_morrison, microphysics_seifert
+
+    USE diagnostic_quantities_mod,                                             &
+        ONLY:  e_s, magnus, q_s, sat, supersaturation, t_l
 
     USE indices,                                                               &
@@ -104,8 +109,4 @@
     INTEGER(iwp) ::  k !<
 
-    REAL(wp) ::  alpha !<
-    REAL(wp) ::  e_s   !<
-    REAL(wp) ::  q_s   !<
-    REAL(wp) ::  t_l   !<
 
     DO  i = nxlg, nxrg
@@ -114,24 +115,7 @@
 
 !
-!--          Compute the liquid water temperature
-             t_l = t_d_pt(k) * pt(k,j,i)
-
-!
-!--          Compute saturation water vapor pressure at t_l
-             e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /          &
-                                                ( t_l - 35.86_wp ) )
-
-!
-!--          Compute approximation of saturation humidity
-             q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
-
-!
-!--          Correction factor
-             alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
-
-!
-!--          Correction of the approximated value
-!--          (see: Cuijpers + Duynkerke, 1993, JAS, 23)
-             q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / ( 1.0_wp + alpha * q_s )
+!--          Call calculation of supersaturation located 
+!--          in diagnostic_quantities_mod
+             CALL supersaturation ( i, j, k )
 
 !
Index: palm/trunk/SOURCE/diagnostic_quantities_mod.f90
===================================================================
--- palm/trunk/SOURCE/diagnostic_quantities_mod.f90	(revision 2608)
+++ palm/trunk/SOURCE/diagnostic_quantities_mod.f90	(revision 2608)
@@ -0,0 +1,126 @@
+!> @file diagnostic_quantities_mod.f90
+!------------------------------------------------------------------------------!
+! This file is part of PALM.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 1997-2017 Leibniz Universitaet Hannover
+!------------------------------------------------------------------------------!
+!
+! Current revisions:
+! -----------------
+! 
+! 
+! Former revisions:
+! -----------------! 
+! $Id$
+! Initial revision
+! 
+! 
+! Description:
+! ------------
+!> This module contains subroutines and functions for the calculation of 
+!> diagnostic quantities. Especially moisture quantities such as the saturation
+!> mixining ratio is calculated
+!------------------------------------------------------------------------------!
+MODULE diagnostic_quantities_mod
+ 
+
+   USE kinds
+   
+   IMPLICIT NONE
+
+   REAL(wp) ::  alpha   !< correction factor 
+   REAL(wp) ::  e_s     !< saturation water vapor pressure
+   REAL(wp) ::  q_s     !< saturation mixing ratio
+   REAL(wp) ::  sat     !< supersaturation
+   REAL(wp) ::  t_l     !< actual temperature
+
+   PRIVATE
+   PUBLIC  e_s, magnus, q_s, sat, supersaturation, t_l
+
+    INTERFACE supersaturation
+       MODULE PROCEDURE supersaturation
+    END INTERFACE supersaturation
+
+ CONTAINS
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Computation of the diagnostic supersaturation sat, actual temperature t_l
+!< and saturation water vapor mixing ratio q_
+!------------------------------------------------------------------------------!
+    SUBROUTINE supersaturation ( i,j,k )
+
+       USE arrays_3d,                                                          &
+           ONLY:  hyp, pt, q, qc, qr
+
+       USE cloud_parameters,                                                   &
+           ONLY:  l_d_cp, l_d_r, t_d_pt                                              
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !<
+       INTEGER(iwp) ::  j                 !<
+       INTEGER(iwp) ::  k                 !<
+!
+!--    Actual liquid water temperature:
+       t_l = t_d_pt(k) * pt(k,j,i)
+!
+!--    Calculate water vapor saturation pressure
+       e_s = magnus( t_l )
+!
+!--    Computation of saturation humidity:
+       q_s   = 0.622_wp * e_s / ( hyp(k) - e_s )
+!
+!--    Correction factor
+       alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
+!
+!--    Correction of the approximated value
+!--    (see: Cuijpers + Duynkerke, 1993, JAS, 23)
+       q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / ( 1.0_wp + alpha * q_s )
+!
+!--    Supersaturation:
+       sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
+
+    END SUBROUTINE supersaturation
+
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This function computes the magnus function (Press et al., 1992).
+!> The magnus formula is needed to calculate the saturation vapor pressure
+!------------------------------------------------------------------------------!
+
+    FUNCTION magnus( temperature )
+
+       IMPLICIT NONE
+
+       REAL(wp)     ::  magnus            !<
+       REAL(wp)     ::  temperature       !<
+
+!
+!--    Saturation vapor pressure at t_l:
+       magnus =  611.2_wp * EXP( 17.62_wp * ( temperature - 273.15_wp ) /      & 
+                                            ( temperature - 29.65_wp  ) )
+
+!        magnus = 610.78_wp * EXP( 17.269_wp * ( temperature - 273.16_wp ) /     &
+!                                              ( temperature - 35.86_wp )        &
+!                                )
+
+    END FUNCTION magnus
+
+
+END MODULE diagnostic_quantities_mod
Index: palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2607)
+++ palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2608)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Calculation of magnus equation in external module (diagnostic_quantities_mod).
+! Adjust calculation of vapor pressure and saturation mixing ratio that it is
+! consistent with formulations in other parts of PALM.
+! 
+! 2575 2017-10-24 09:57:58Z maronga
 ! Pavement parameterization revised
 ! 
@@ -1410,5 +1415,7 @@
  SUBROUTINE lsm_energy_balance( horizontal, l )
 
- 
+    USE diagnostic_quantities_mod,                                             &
+        ONLY:  magnus 
+
     USE pegrid
 
@@ -1648,8 +1655,9 @@
           f2 = 1.0E-20_wp
        ENDIF
-!
-!--    Calculate water vapour pressure at saturation 
-       e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( surf_t_surface%var_1d(m) &
-                     - 273.16_wp ) / ( surf_t_surface%var_1d(m) - 35.86_wp ) )
+
+!
+!--    Calculate water vapour pressure at saturation and convert to hPa
+       e_s = 0.01_wp * magnus( surf_t_surface%var_1d(m) )
+
 !
 !--    f3: correction for vapour pressure deficit
@@ -1657,5 +1665,5 @@
 !
 !--       Calculate vapour pressure
-          e  = qv1 * surface_pressure / 0.622_wp
+          e  = qv1 * surface_pressure / ( qv1 + 0.622_wp )
           f3 = EXP ( - surf%g_d(m) * (e_s - e) )
        ELSE
@@ -1698,5 +1706,5 @@
 !
 !--    Calculate saturation specific humidity
-       q_s = 0.622_wp * e_s / surface_pressure
+       q_s = 0.622_wp * e_s / ( surface_pressure - e_s )
 !
 !--    In case of dewfall, set evapotranspiration to zero
@@ -1733,9 +1741,9 @@
 !
 !--    Calculate derivative of q_s for Taylor series expansion
-       e_s_dt = e_s * ( 17.269_wp / ( surf_t_surface%var_1d(m) - 35.86_wp) -   &
-                        17.269_wp*( surf_t_surface%var_1d(m) - 273.16_wp)      &
-                       / ( surf_t_surface%var_1d(m) - 35.86_wp)**2 )
-
-       dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
+       e_s_dt = e_s * ( 17.62_wp / ( surf_t_surface%var_1d(m) - 29.65_wp) -   &
+                        17.62_wp*( surf_t_surface%var_1d(m) - 273.15_wp)      &
+                       / ( surf_t_surface%var_1d(m) - 29.65_wp)**2 )
+
+       dq_s_dt = 0.622_wp * e_s_dt / ( surface_pressure - e_s_dt )
 !
 !--    Add LW up so that it can be removed in prognostic equation
@@ -1988,4 +1996,6 @@
     SUBROUTINE calc_q_surface
 
+       USE diagnostic_quantities_mod
+
        IMPLICIT NONE
 
@@ -1997,15 +2007,11 @@
           j   = surf%j(m)
           k   = surf%k(m)
-
-!
-!--       Calculate water vapour pressure at saturation
-          e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp *                         &
-                                 ( surf_t_surface_p%var_1d(m) - 273.16_wp ) /  &
-                                 ( surf_t_surface_p%var_1d(m) - 35.86_wp  )    &
-                                         )
+!
+!--       Calculate water vapour pressure at saturation and convert to hPa
+          e_s = 0.01_wp * magnus( surf_t_surface_p%var_1d(m) )                   
 
 !
 !--       Calculate specific humidity at saturation
-          q_s = 0.622_wp * e_s / surface_pressure 
+          q_s = 0.622_wp * e_s / ( surface_pressure - e_s )
 
           resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) )
Index: palm/trunk/SOURCE/lpm_droplet_condensation.f90
===================================================================
--- palm/trunk/SOURCE/lpm_droplet_condensation.f90	(revision 2607)
+++ palm/trunk/SOURCE/lpm_droplet_condensation.f90	(revision 2608)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Calculation of magnus equation in external module (diagnostic_quantities_mod).
+! 
+! 2375 2017-08-29 14:10:28Z schwenkel
 ! Changed ONLY-dependencies
 !
@@ -129,4 +132,7 @@
     USE cpulog,                                                                &
         ONLY:  cpu_log, log_point_s
+
+    USE diagnostic_quantities_mod,                                             &
+        ONLY:  magnus
 
     USE grid_variables,                                                        &
@@ -199,5 +205,5 @@
 !
 !-- Saturation vapor pressure (Eq. 10 in Bolton, 1980)
-    e_s = 611.2_wp * EXP( 17.62_wp * ( t_int - 273.15_wp ) / ( t_int - 29.65_wp ) )
+    e_s = magnus( t_int )
 !
 !-- Current vapor pressure
Index: palm/trunk/SOURCE/lpm_init.f90
===================================================================
--- palm/trunk/SOURCE/lpm_init.f90	(revision 2607)
+++ palm/trunk/SOURCE/lpm_init.f90	(revision 2608)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Calculation of magnus equation in external module (diagnostic_quantities_mod).
+! 
+! 2606 2017-11-10 10:36:31Z schwenkel
 ! Changed particle box locations: center of particle box now coincides 
 ! with scalar grid point of same index.
@@ -986,4 +989,8 @@
     USE constants,                                                             &
         ONLY: pi
+
+    USE diagnostic_quantities_mod,                                             &
+        ONLY:  magnus
+
 
     USE kinds
@@ -1124,6 +1131,6 @@
              t_int  = pt(kp,jp,ip) * ( hyp(kp) / 100000.0_wp )**0.286_wp
 
-             e_s = 611.0_wp * EXP( l_d_rv * ( 3.6609E-3_wp - 1.0_wp / t_int ) )
-             e_a = q(kp,jp,ip) * hyp(kp) / ( 0.378_wp * q(kp,jp,ip) + 0.622_wp )
+             e_s = magnus( t_int )
+             e_a = q(kp,jp,ip) * hyp(kp) / ( q(kp,jp,ip) + 0.622_wp )
 
              sigma   = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
Index: palm/trunk/SOURCE/microphysics_mod.f90
===================================================================
--- palm/trunk/SOURCE/microphysics_mod.f90	(revision 2607)
+++ palm/trunk/SOURCE/microphysics_mod.f90	(revision 2608)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Calculation of supersaturation in external module (diagnostic_quantities_mod).
+! Change: correct calculation of saturation specific humidity to saturation
+! mixing ratio (the factor of 0.378 vanishes).
+! 
+! 2522 2017-10-05 14:20:37Z schwenkel
 ! Minor bugfix
 ! 
@@ -548,4 +553,7 @@
            ONLY:  cpu_log, log_point_s
 
+       USE diagnostic_quantities_mod,                                          &
+           ONLY: e_s, magnus, q_s, sat, supersaturation, t_l
+
        USE indices,                                                            &
            ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
@@ -564,19 +572,13 @@
        REAL(wp)     ::  activ             !<
        REAL(wp)     ::  afactor           !<
-       REAL(wp)     ::  alpha             !<
        REAL(wp)     ::  beta_act          !<
        REAL(wp)     ::  bfactor           !<
-       REAL(wp)     ::  e_s               !<
        REAL(wp)     ::  k_act             !<
        REAL(wp)     ::  n_act             !<
        REAL(wp)     ::  n_ccn             !<
-       REAL(wp)     ::  q_s               !<
        REAL(wp)     ::  s_0               !<
-       REAL(wp)     ::  sat               !<
        REAL(wp)     ::  sat_max           !<
        REAL(wp)     ::  sigma             !<
        REAL(wp)     ::  sigma_act         !<
-       REAL(wp)     ::  t_int             !<
-       REAL(wp)     ::  t_l               !<
 
        CALL cpu_log( log_point_s(65), 'activation', 'start' )
@@ -587,25 +589,7 @@
 
 !
-!--             Actual liquid water temperature:
-                t_l = t_d_pt(k) * pt(k,j,i)
-
-!
-!--             Calculate actual temperature
-                t_int = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
-!
-!--             Saturation vapor pressure at t_l:
-                e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
-                                       ( t_l - 35.86_wp )                   &
-                                     )
-!
-!--             Computation of saturation humidity:
-                q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
-                alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
-                q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /               &
-                        ( 1.0_wp + alpha * q_s )
-
-!--             Supersaturation:
-                sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
-
+!--             Call calculation of supersaturation located 
+!--             in diagnostic_quantities_mod
+                CALL supersaturation ( i, j, k )
 !
 !--             Prescribe parameters for activation
@@ -636,6 +620,6 @@
 !--                Curvature effect (afactor) with surface tension
 !--                parameterization by Straka (2009)
-                   sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
-                   afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
+                   sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
+                   afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
 !
 !--                Solute effect (bfactor)
@@ -695,4 +679,7 @@
            ONLY:  cpu_log, log_point_s
 
+       USE diagnostic_quantities_mod,                                          &
+           ONLY: e_s, magnus, q_s, sat, supersaturation, t_l
+
        USE indices,                                                            &
            ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
@@ -709,16 +696,11 @@
        INTEGER(iwp) ::  k                 !<
 
-       REAL(wp)     ::  alpha             !<
        REAL(wp)     ::  cond              !<
        REAL(wp)     ::  cond_max          !<
        REAL(wp)     ::  dc                !<
-       REAL(wp)     ::  e_s               !<
        REAL(wp)     ::  evap              !<
        REAL(wp)     ::  evap_nc           !<
        REAL(wp)     ::  g_fac             !<
        REAL(wp)     ::  nc_0              !<
-       REAL(wp)     ::  q_s               !<
-       REAL(wp)     ::  sat               !<
-       REAL(wp)     ::  t_l               !<
        REAL(wp)     ::  temp              !<
        REAL(wp)     ::  xc                !<
@@ -730,21 +712,7 @@
              DO  k = nzb+1, nzt
 !
-!--             Actual liquid water temperature:
-                t_l = t_d_pt(k) * pt(k,j,i)
-!
-!--             Saturation vapor pressure at t_l:
-                e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /       &
-                                       ( t_l - 35.86_wp )                      &
-                                     )
-!
-!--             Computation of saturation humidity:
-                q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
-                alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
-                q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /                  &
-                        ( 1.0_wp + alpha * q_s )
-
-!--             Supersaturation:
-                sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
-
+!--             Call calculation of supersaturation located 
+!--             in diagnostic_quantities_mod
+                CALL supersaturation ( i, j, k )
 !
 !--             Actual temperature:
@@ -1200,4 +1168,7 @@
            ONLY:  cpu_log, log_point_s
 
+       USE diagnostic_quantities_mod,                                          &
+           ONLY: e_s, magnus, q_s, sat, supersaturation, t_l
+
        USE indices,                                                            &
            ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
@@ -1211,7 +1182,5 @@
        INTEGER(iwp) ::  k                 !<
 
-       REAL(wp)     ::  alpha             !<
        REAL(wp)     ::  dr                !<
-       REAL(wp)     ::  e_s               !<
        REAL(wp)     ::  evap              !<
        REAL(wp)     ::  evap_nr           !<
@@ -1224,7 +1193,4 @@
        REAL(wp)     ::  mu_r_5d2          !<
        REAL(wp)     ::  nr_0              !<
-       REAL(wp)     ::  q_s               !<
-       REAL(wp)     ::  sat               !<
-       REAL(wp)     ::  t_l               !<
        REAL(wp)     ::  temp              !<
        REAL(wp)     ::  xr                !<
@@ -1240,21 +1206,9 @@
 
                 IF ( qr(k,j,i) > eps_sb )  THEN
-!
-!--                Actual liquid water temperature:
-                   t_l = t_d_pt(k) * pt(k,j,i)
-!
-!--                Saturation vapor pressure at t_l:
-                   e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
-                                          ( t_l - 35.86_wp )                   &
-                                        )
-!
-!--                Computation of saturation humidity:
-                   q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
-                   alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
-                   q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /               &
-                           ( 1.0_wp + alpha * q_s )
-!
-!--                Supersaturation:
-                   sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
+
+!
+!--                Call calculation of supersaturation located 
+!--                in diagnostic_quantities_mod
+                   CALL supersaturation ( i, j, k )
 !
 !--                Evaporation needs only to be calculated in subsaturated regions
