Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 1932)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 1960)
@@ -19,5 +19,9 @@
 ! Current revisions:
 ! -----------------
-! 
+! Separate checks and calculations for humidity and passive scalar. Therefore,
+! a subroutine to calculate vertical gradients is introduced, in order  to reduce
+! redundance. 
+! Additional check - large-scale forcing is not implemented for passive scalars
+! so far.
 ! 
 ! Former revisions:
@@ -1053,10 +1057,4 @@
     ENDIF
 
-    IF ( passive_scalar  .AND.  humidity )  THEN
-       message_string = 'humidity = .TRUE. and passive_scalar = .TRUE. ' //    &
-                        'is not allowed simultaneously'
-       CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 )
-    ENDIF
-
 !
 !-- When plant canopy model is used, peform addtional checks
@@ -1094,9 +1092,7 @@
 !--    Initial profiles for 1D and 3D model, respectively (u,v further below)
        pt_init = pt_surface
-       IF ( humidity )  THEN
-          q_init  = q_surface
-       ENDIF
-       IF ( ocean )           sa_init = sa_surface
-       IF ( passive_scalar )  q_init  = s_surface
+       IF ( humidity       )  q_init  = q_surface
+       IF ( ocean          )  sa_init = sa_surface
+       IF ( passive_scalar )  s_init  = s_surface
 
 !
@@ -1282,191 +1278,33 @@
 !--    Compute initial temperature profile using the given temperature gradients
        IF (  .NOT.  neutral )  THEN
-
-          i = 1
-          gradient = 0.0_wp
-
-          IF (  .NOT.  ocean )  THEN
-
-             pt_vertical_gradient_level_ind(1) = 0
-             DO  k = 1, nzt+1
-                IF ( i < 11 )  THEN
-                   IF ( pt_vertical_gradient_level(i) < zu(k)  .AND.           &
-                        pt_vertical_gradient_level(i) >= 0.0_wp )  THEN
-                      gradient = pt_vertical_gradient(i) / 100.0_wp
-                      pt_vertical_gradient_level_ind(i) = k - 1
-                      i = i + 1
-                   ENDIF
-                ENDIF
-                IF ( gradient /= 0.0_wp )  THEN
-                   IF ( k /= 1 )  THEN
-                      pt_init(k) = pt_init(k-1) + dzu(k) * gradient
-                   ELSE
-                      pt_init(k) = pt_surface   + dzu(k) * gradient
-                   ENDIF
-                ELSE
-                   pt_init(k) = pt_init(k-1)
-                ENDIF
-             ENDDO
-
-          ELSE
-
-             pt_vertical_gradient_level_ind(1) = nzt+1
-             DO  k = nzt, 0, -1
-                IF ( i < 11 )  THEN
-                   IF ( pt_vertical_gradient_level(i) > zu(k)  .AND.           &
-                        pt_vertical_gradient_level(i) <= 0.0_wp )  THEN
-                      gradient = pt_vertical_gradient(i) / 100.0_wp
-                      pt_vertical_gradient_level_ind(i) = k + 1
-                      i = i + 1
-                   ENDIF
-                ENDIF
-                IF ( gradient /= 0.0_wp )  THEN
-                   IF ( k /= nzt )  THEN
-                      pt_init(k) = pt_init(k+1) - dzu(k+1) * gradient
-                   ELSE
-                      pt_init(k)   = pt_surface - 0.5_wp * dzu(k+1) * gradient
-                      pt_init(k+1) = pt_surface + 0.5_wp * dzu(k+1) * gradient
-                   ENDIF
-                ELSE
-                   pt_init(k) = pt_init(k+1)
-                ENDIF
-             ENDDO
-
-          ENDIF
-
-       ENDIF
-
-!
-!--    In case of no given temperature gradients, choose gradient of neutral
-!--    stratification
-       IF ( pt_vertical_gradient_level(1) == -9999999.9_wp )  THEN
-          pt_vertical_gradient_level(1) = 0.0_wp
-       ENDIF
-
-!
-!--    Store temperature gradient at the top boundary for possible Neumann
-!--    boundary condition
-       bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
-
-!
-!--    If required, compute initial humidity or scalar profile using the given
-!--    humidity/scalar gradient. In case of scalar transport, initially store
-!--    values of the scalar parameters on humidity parameters
+          CALL init_vertical_profiles( pt_vertical_gradient_level_ind,          &
+                                       pt_vertical_gradient_level,              &
+                                       pt_vertical_gradient, pt_init,           &
+                                       pt_surface, bc_pt_t_val )
+       ENDIF
+!
+!--    Compute initial humidity profile using the given humidity gradients
+       IF ( humidity )  THEN
+          CALL init_vertical_profiles( q_vertical_gradient_level_ind,          &
+                                       q_vertical_gradient_level,              &
+                                       q_vertical_gradient, q_init,            &
+                                       q_surface, bc_q_t_val )
+       ENDIF
+!
+!--    Compute initial scalar profile using the given scalar gradients
        IF ( passive_scalar )  THEN
-          bc_q_b                    = bc_s_b
-          bc_q_t                    = bc_s_t
-          q_surface                 = s_surface
-          q_surface_initial_change  = s_surface_initial_change
-          q_vertical_gradient       = s_vertical_gradient
-          q_vertical_gradient_level = s_vertical_gradient_level
-          surface_waterflux         = surface_scalarflux
-          wall_humidityflux         = wall_scalarflux
-       ENDIF
-
-       IF ( humidity  .OR.  passive_scalar )  THEN
-
-          i = 1
-          gradient = 0.0_wp
-
-          IF (  .NOT.  ocean )  THEN
-
-             q_vertical_gradient_level_ind(1) = 0
-             DO  k = 1, nzt+1
-                IF ( i < 11 )  THEN
-                   IF ( q_vertical_gradient_level(i) < zu(k)  .AND.            &
-                        q_vertical_gradient_level(i) >= 0.0_wp )  THEN
-                      gradient = q_vertical_gradient(i) / 100.0_wp
-                      q_vertical_gradient_level_ind(i) = k - 1
-                      i = i + 1
-                   ENDIF
-                ENDIF
-                IF ( gradient /= 0.0_wp )  THEN
-                   IF ( k /= 1 )  THEN
-                      q_init(k) = q_init(k-1) + dzu(k) * gradient
-                   ELSE
-                      q_init(k) = q_init(k-1) + dzu(k) * gradient
-                   ENDIF
-                ELSE
-                   q_init(k) = q_init(k-1)
-                ENDIF
-   !
-   !--          Avoid negative humidities
-                IF ( q_init(k) < 0.0_wp )  THEN
-                   q_init(k) = 0.0_wp
-                ENDIF
-             ENDDO
-
-          ELSE
-
-             q_vertical_gradient_level_ind(1) = nzt+1
-             DO  k = nzt, 0, -1
-                IF ( i < 11 )  THEN
-                   IF ( q_vertical_gradient_level(i) > zu(k)  .AND.            &
-                        q_vertical_gradient_level(i) <= 0.0_wp )  THEN
-                      gradient = q_vertical_gradient(i) / 100.0_wp
-                      q_vertical_gradient_level_ind(i) = k + 1
-                      i = i + 1
-                   ENDIF
-                ENDIF
-                IF ( gradient /= 0.0_wp )  THEN
-                   IF ( k /= nzt )  THEN
-                      q_init(k) = q_init(k+1) - dzu(k+1) * gradient
-                   ELSE
-                      q_init(k)   = q_surface - 0.5_wp * dzu(k+1) * gradient
-                      q_init(k+1) = q_surface + 0.5_wp * dzu(k+1) * gradient
-                   ENDIF
-                ELSE
-                   q_init(k) = q_init(k+1)
-                ENDIF
-   !
-   !--          Avoid negative humidities
-                IF ( q_init(k) < 0.0_wp )  THEN
-                   q_init(k) = 0.0_wp
-                ENDIF
-             ENDDO
-
-          ENDIF
-!
-!--       In case of no given humidity gradients, choose zero gradient
-!--       conditions
-          IF ( q_vertical_gradient_level(1) == -1.0_wp )  THEN
-             q_vertical_gradient_level(1) = 0.0_wp
-          ENDIF
-!
-!--       Store humidity, rain water content and rain drop concentration 
-!--       gradient at the top boundary for possile Neumann boundary condition
-          bc_q_t_val  = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
-       ENDIF
-
+          CALL init_vertical_profiles( s_vertical_gradient_level_ind,          &
+                                       s_vertical_gradient_level,              &
+                                       s_vertical_gradient, s_init,            &
+                                       s_surface, bc_s_t_val )
+       ENDIF
 !
 !--    If required, compute initial salinity profile using the given salinity
 !--    gradients
        IF ( ocean )  THEN
-
-          i = 1
-          gradient = 0.0_wp
-
-          sa_vertical_gradient_level_ind(1) = nzt+1
-          DO  k = nzt, 0, -1
-             IF ( i < 11 )  THEN
-                IF ( sa_vertical_gradient_level(i) > zu(k)  .AND.              &
-                     sa_vertical_gradient_level(i) <= 0.0_wp )  THEN
-                   gradient = sa_vertical_gradient(i) / 100.0_wp
-                   sa_vertical_gradient_level_ind(i) = k + 1
-                   i = i + 1
-                ENDIF
-             ENDIF
-             IF ( gradient /= 0.0_wp )  THEN
-                IF ( k /= nzt )  THEN
-                   sa_init(k) = sa_init(k+1) - dzu(k+1) * gradient
-                ELSE
-                   sa_init(k)   = sa_surface - 0.5_wp * dzu(k+1) * gradient
-                   sa_init(k+1) = sa_surface + 0.5_wp * dzu(k+1) * gradient
-                ENDIF
-             ELSE
-                sa_init(k) = sa_init(k+1)
-             ENDIF
-          ENDDO
-
+          CALL init_vertical_profiles( sa_vertical_gradient_level_ind,          &
+                                       sa_vertical_gradient_level,              &
+                                       sa_vertical_gradient, s_init,            &
+                                       sa_surface, -999.0_wp )
        ENDIF
 
@@ -1860,32 +1698,10 @@
 
 !
-!-- In case of humidity or passive scalar, set boundary conditions for total 
-!-- water content / scalar
-    IF ( humidity  .OR.  passive_scalar ) THEN
-       IF ( humidity )  THEN
-          sq = 'q'
-       ELSE
-          sq = 's'
-       ENDIF
-       IF ( bc_q_b == 'dirichlet' )  THEN
-          ibc_q_b = 0
-       ELSEIF ( bc_q_b == 'neumann' )  THEN
-          ibc_q_b = 1
-       ELSE
-          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
-                           '_b ="' // TRIM( bc_q_b ) // '"'
-          CALL message( 'check_parameters', 'PA0071', 1, 2, 0, 6, 0 )
-       ENDIF
-       IF ( bc_q_t == 'dirichlet' )  THEN
-          ibc_q_t = 0
-       ELSEIF ( bc_q_t == 'neumann' )  THEN
-          ibc_q_t = 1
-       ELSEIF ( bc_q_t == 'nested' )  THEN
-          ibc_q_t = 3
-       ELSE
-          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
-                           '_t ="' // TRIM( bc_q_t ) // '"'
-          CALL message( 'check_parameters', 'PA0072', 1, 2, 0, 6, 0 )
-       ENDIF
+!-- Set boundary conditions for total water content
+    IF ( humidity )  THEN
+       CALL check_bc_scalars( 'q', bc_q_b, bc_q_t, ibc_q_b, ibc_q_t,           &
+                              'PA0071', 'PA0072', 'PA0073', 'PA0074',          &
+                              surface_waterflux, constant_waterflux,           &
+                              q_surface_initial_change )
 
        IF ( surface_waterflux == 9999999.9_wp  )  THEN
@@ -1909,21 +1725,11 @@
        ENDIF
 
-!
-!--    A given surface humidity implies Dirichlet boundary condition for
-!--    humidity. In this case specification of a constant water flux is
-!--    forbidden.
-       IF ( ibc_q_b == 0  .AND.  constant_waterflux )  THEN
-          message_string = 'boundary condition: bc_' // TRIM( sq ) // '_b ' // &
-                           '= "' // TRIM( bc_q_b ) // '" is not allowed wi' // &
-                           'th prescribed surface flux'
-          CALL message( 'check_parameters', 'PA0073', 1, 2, 0, 6, 0 )
-       ENDIF
-       IF ( constant_waterflux  .AND.  q_surface_initial_change /= 0.0_wp )  THEN
-          WRITE( message_string, * )  'a prescribed surface flux is not allo', &
-                 'wed with ', sq, '_surface_initial_change (/=0) = ',          &
-                 q_surface_initial_change
-          CALL message( 'check_parameters', 'PA0074', 1, 2, 0, 6, 0 )
-       ENDIF
-
+    ENDIF
+    
+    IF ( passive_scalar )  THEN
+       CALL check_bc_scalars( 's', bc_s_b, bc_s_t, ibc_s_b, ibc_s_t,           &
+                              'PA0071', 'PA0072', 'PA0073', 'PA0074',          &
+                              surface_scalarflux, constant_scalarflux,         &
+                              s_surface_initial_change )
     ENDIF
 !
@@ -3850,5 +3656,5 @@
                         'the usage of large scale forcing from external file.'
        CALL message( 'check_parameters', 'PA0375', 1, 2, 0, 6, 0 )
-     ENDIF
+    ENDIF
 
     IF ( large_scale_forcing  .AND.  (  .NOT.  humidity ) )  THEN
@@ -3856,5 +3662,11 @@
                         'file LSF_DATA requires humidity = .T..'
        CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 )
-     ENDIF
+    ENDIF
+
+    IF ( large_scale_forcing  .AND.  passive_scalar )  THEN
+       message_string = 'The usage of large scale forcing from external &'//   & 
+                        'file LSF_DATA is not implemented for psassive scalars'
+       CALL message( 'check_parameters', 'PA0440', 1, 2, 0, 6, 0 )
+    ENDIF
 
     IF ( large_scale_forcing  .AND.  topography /= 'flat' )  THEN
@@ -3864,5 +3676,5 @@
     ENDIF
 
-    IF ( large_scale_forcing  .AND.  ocean  )  THEN
+    IF ( large_scale_forcing  .AND.  ocean )  THEN
        message_string = 'The usage of large scale forcing from external &'//   & 
                         'file LSF_DATA is not implemented for ocean runs'
@@ -3943,4 +3755,182 @@
 
     END SUBROUTINE check_dt_do
+    
+    
+    
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Inititalizes the vertical profiles of scalar quantities.
+!------------------------------------------------------------------------------!
+
+    SUBROUTINE init_vertical_profiles( vertical_gradient_level_ind,            &
+                                       vertical_gradient_level,                &
+                                       vertical_gradient,                      &
+                                       pr_init, surface_value, bc_t_val )
+
+
+       IMPLICIT NONE    
+    
+       INTEGER(iwp) ::  i     !< counter
+       INTEGER(iwp), DIMENSION(1:10) ::  vertical_gradient_level_ind !< vertical grid indices for gradient levels
+       
+       REAL(wp)     ::  bc_t_val      !< model top gradient       
+       REAL(wp)     ::  gradient      !< vertica gradient of the respective quantity 
+       REAL(wp)     ::  surface_value !< surface value of the respecitve quantity
+
+       
+       REAL(wp), DIMENSION(0:nz+1) ::  pr_init                 !< initialisation profile
+       REAL(wp), DIMENSION(1:10)   ::  vertical_gradient       !< given vertical gradient
+       REAL(wp), DIMENSION(1:10)   ::  vertical_gradient_level !< given vertical gradient level
+       
+       i = 1
+       gradient = 0.0_wp
+
+       IF (  .NOT.  ocean )  THEN
+
+          vertical_gradient_level_ind(1) = 0
+          DO  k = 1, nzt+1
+             IF ( i < 11 )  THEN
+                IF ( vertical_gradient_level(i) < zu(k)  .AND.            &
+                     vertical_gradient_level(i) >= 0.0_wp )  THEN
+                   gradient = vertical_gradient(i) / 100.0_wp
+                   vertical_gradient_level_ind(i) = k - 1
+                   i = i + 1
+                ENDIF
+             ENDIF
+             IF ( gradient /= 0.0_wp )  THEN
+                IF ( k /= 1 )  THEN
+                   pr_init(k) = pr_init(k-1) + dzu(k) * gradient
+                ELSE
+                   pr_init(k) = pr_init(k-1) + dzu(k) * gradient
+                ENDIF
+             ELSE
+                pr_init(k) = pr_init(k-1)
+             ENDIF
+   !
+   !--       Avoid negative values
+             IF ( pr_init(k) < 0.0_wp )  THEN
+                pr_init(k) = 0.0_wp
+             ENDIF
+          ENDDO
+
+       ELSE
+
+          vertical_gradient_level_ind(1) = nzt+1
+          DO  k = nzt, 0, -1
+             IF ( i < 11 )  THEN
+                IF ( vertical_gradient_level(i) > zu(k)  .AND.            &
+                     vertical_gradient_level(i) <= 0.0_wp )  THEN
+                   gradient = vertical_gradient(i) / 100.0_wp
+                   vertical_gradient_level_ind(i) = k + 1
+                   i = i + 1
+                ENDIF
+             ENDIF
+             IF ( gradient /= 0.0_wp )  THEN
+                IF ( k /= nzt )  THEN
+                   pr_init(k) = pr_init(k+1) - dzu(k+1) * gradient
+                ELSE
+                   pr_init(k)   = surface_value - 0.5_wp * dzu(k+1) * gradient
+                   pr_init(k+1) = surface_value + 0.5_wp * dzu(k+1) * gradient
+                ENDIF
+             ELSE
+                pr_init(k) = pr_init(k+1)
+             ENDIF
+   !
+   !--       Avoid negative humidities
+             IF ( pr_init(k) < 0.0_wp )  THEN
+                pr_init(k) = 0.0_wp
+             ENDIF
+          ENDDO
+
+       ENDIF
+       
+!
+!--    In case of no given gradients, choose zero gradient conditions
+       IF ( vertical_gradient_level(1) == -1.0_wp )  THEN
+          vertical_gradient_level(1) = 0.0_wp
+       ENDIF
+!
+!--    Store gradient at the top boundary for possible Neumann boundary condition
+       bc_t_val  = ( pr_init(nzt+1) - pr_init(nzt) ) / dzu(nzt+1)
+    
+    END SUBROUTINE init_vertical_profiles
+    
+    
+     
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check the bottom and top boundary conditions for humidity and scalars.
+!------------------------------------------------------------------------------!
+
+    SUBROUTINE check_bc_scalars( sq, bc_b, bc_t, ibc_b, ibc_t,                 &
+                                 err_nr_b, err_nr_t, err_nr_3, err_nr_4,       &
+                                 surface_flux, constant_flux, surface_initial_change )
+
+
+       IMPLICIT NONE    
+    
+       CHARACTER (LEN=1)   ::  sq                       !< 
+       CHARACTER (LEN=*)   ::  bc_b
+       CHARACTER (LEN=*)   ::  bc_t
+       CHARACTER (LEN=*)   ::  err_nr_b
+       CHARACTER (LEN=*)   ::  err_nr_t
+       CHARACTER (LEN=*)   ::  err_nr_3
+       CHARACTER (LEN=*)   ::  err_nr_4
+       
+       INTEGER(iwp)        ::  ibc_b
+       INTEGER(iwp)        ::  ibc_t
+       
+       LOGICAL             ::  constant_flux
+       
+       REAL(wp)            ::  surface_flux
+       REAL(wp)            ::  surface_initial_change
+       
+
+       IF ( bc_b == 'dirichlet' )  THEN
+          ibc_b = 0
+       ELSEIF ( bc_b == 'neumann' )  THEN
+          ibc_b = 1
+       ELSE
+          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
+                           '_b ="' // TRIM( bc_b ) // '"'
+          CALL message( 'check_parameters', err_nr_b, 1, 2, 0, 6, 0 )
+       ENDIF
+       
+       IF ( bc_t == 'dirichlet' )  THEN
+          ibc_t = 0
+       ELSEIF ( bc_t == 'neumann' )  THEN
+          ibc_t = 1
+       ELSEIF ( bc_t == 'nested' )  THEN
+          ibc_t = 3
+       ELSE
+          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
+                           '_t ="' // TRIM( bc_t ) // '"'
+          CALL message( 'check_parameters', err_nr_t, 1, 2, 0, 6, 0 )
+       ENDIF
+  
+!
+!--    A given surface value implies Dirichlet boundary condition for
+!--    the respective quantity. In this case specification of a constant flux is
+!--    forbidden.
+       IF ( ibc_b == 0  .AND.  constant_flux )  THEN
+          message_string = 'boundary condition: bc_' // TRIM( sq ) // '_b ' // &
+                           '= "' // TRIM( bc_b ) // '" is not allowed wi' // &
+                           'th prescribed surface flux'
+          CALL message( 'check_parameters', err_nr_3, 1, 2, 0, 6, 0 )
+       ENDIF
+       IF ( constant_waterflux  .AND.  surface_initial_change /= 0.0_wp )  THEN
+          WRITE( message_string, * )  'a prescribed surface flux is not allo', &
+                 'wed with ', sq, '_surface_initial_change (/=0) = ',          &
+                 surface_initial_change
+          CALL message( 'check_parameters', err_nr_4, 1, 2, 0, 6, 0 )
+       ENDIF  
+       
+    
+    END SUBROUTINE check_bc_scalars   
+    
+    
 
  END SUBROUTINE check_parameters
