Changeset 1960

Timestamp:

Jul 12, 2016 4:34:24 PM (8 years ago)

Author:

suehring

Message:

Separate balance equations for humidity and passive_scalar

Location:

palm/trunk/SOURCE

Files:

• advec_s_bc.f90 (modified) (3 diffs)
• advec_ws.f90 (modified) (12 diffs)
• average_3d_data.f90 (modified) (2 diffs)
• boundary_conds.f90 (modified) (9 diffs)
• check_parameters.f90 (modified) (10 diffs)
• data_output_2d.f90 (modified) (24 diffs)
• data_output_3d.f90 (modified) (3 diffs)
• data_output_dvrp.f90 (modified) (5 diffs)
• data_output_mask.f90 (modified) (3 diffs)
• data_output_spectra.f90 (modified) (2 diffs)
• flow_statistics.f90 (modified) (36 diffs)
• header.f90 (modified) (10 diffs)
• inflow_turbulence.f90 (modified) (9 diffs)
• init_1d_model.f90 (modified) (5 diffs)
• init_3d_model.f90 (modified) (34 diffs)
• modules.f90 (modified) (13 diffs)
• netcdf_interface_mod.f90 (modified) (4 diffs)
• palm.f90 (modified) (2 diffs)
• parin.f90 (modified) (3 diffs)
• plant_canopy_model_mod.f90 (modified) (6 diffs)
• prognostic_equations.f90 (modified) (14 diffs)
• read_3d_binary.f90 (modified) (3 diffs)
• read_var_list.f90 (modified) (7 diffs)
• spectra_mod.f90 (modified) (3 diffs)
• sum_up_3d_data.f90 (modified) (4 diffs)
• surface_layer_fluxes_mod.f90 (modified) (7 diffs)
• swap_timelevel.f90 (modified) (8 diffs)
• time_integration.f90 (modified) (9 diffs)
• user_actions.f90 (modified) (3 diffs)
• user_spectra.f90 (modified) (3 diffs)
• virtual_flight_mod.f90 (modified) (7 diffs)
• write_3d_binary.f90 (modified) (3 diffs)
• write_var_list.f90 (modified) (2 diffs)

palm/trunk/SOURCE/advec_s_bc.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! New CASE statement to treat scalars and humidity separately
 ! 
 ! Former revisions:
@@ -120 +120 @@
 
        USE control_parameters,                                                    &
-           ONLY:  dt_3d, bc_pt_t_val, bc_q_t_val, ibc_pt_b, ibc_pt_t, ibc_q_t,    &
-                  message_string, pt_slope_offset, sloping_surface, u_gtrans,     &
-                  v_gtrans
+           ONLY:  dt_3d, bc_pt_t_val, bc_q_t_val, bc_s_t_val, ibc_pt_b, ibc_pt_t, &
+                  ibc_q_t, ibc_s_t, message_string, pt_slope_offset,              &
+                  sloping_surface, u_gtrans, v_gtrans
 
        USE cpulog,                                                                &
@@ -969 +969 @@
           ENDIF
 
+       ELSEIF ( sk_char == 's' )  THEN
+!
+!--       Specific scalar boundary condition at the bottom boundary.
+!--       Dirichlet (fixed surface humidity) or Neumann (i.e. zero gradient)
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
+                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
+                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
+             ENDDO
+          ENDDO
+
+!
+!--       Specific scalar boundary condition at the top boundary
+          IF ( ibc_s_t == 0 )  THEN
+!
+!--          Dirichlet
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
+                   sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
+                ENDDO
+             ENDDO
+
+          ELSE
+!
+!--          Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i) + bc_s_t_val * dzu(nzt+1)
+                   sk_p(nzt+3,j,i)   = sk_p(nzt+2,j,i) + bc_s_t_val * dzu(nzt+1)
+                ENDDO
+             ENDDO
+
+          ENDIF
+
        ELSEIF ( sk_char == 'qr' )  THEN
 

palm/trunk/SOURCE/advec_ws.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Separate humidity and passive scalar
 ! 
 ! Former revisions:
@@ -228 +228 @@
        USE arrays_3d,                                                          &
            ONLY:  diss_l_e, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr,         &
-                  diss_l_sa, diss_l_u, diss_l_v, diss_l_w,  flux_l_e,          &
-                  flux_l_nr, flux_l_pt, flux_l_q, flux_l_qr, flux_l_sa,        &
-                  flux_l_u, flux_l_v, flux_l_w, diss_s_e, diss_s_nr, diss_s_pt,&
-                  diss_s_q, diss_s_qr, diss_s_sa, diss_s_u, diss_s_v, diss_s_w,& 
-                  flux_s_e, flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr,         &
-                  flux_s_sa, flux_s_u, flux_s_v, flux_s_w
+                  diss_l_s, diss_l_sa, diss_l_u, diss_l_v, diss_l_w,  flux_l_e,&
+                  flux_l_nr, flux_l_pt, flux_l_q, flux_l_qr, flux_l_s,         &
+                  flux_l_sa, flux_l_u, flux_l_v, flux_l_w, diss_s_e, diss_s_nr,&
+                  diss_s_pt, diss_s_q, diss_s_qr, diss_s_s, diss_s_sa,         &
+                  diss_s_u, diss_s_v, diss_s_w, flux_s_e, flux_s_nr, flux_s_pt,&
+                  flux_s_q, flux_s_qr, flux_s_s, flux_s_sa, flux_s_u, flux_s_v,&
+                  flux_s_w
 
        USE constants,                                                          &
@@ -254 +255 @@
            ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsnrs_ws_l,&
                   sums_wspts_ws_l, sums_wsqrs_ws_l, sums_wsqs_ws_l,            &
-                  sums_wssas_ws_l,  sums_wsus_ws_l, sums_wsvs_ws_l  
+                  sums_wsss_ws_l, sums_wssas_ws_l,  sums_wsss_ws_l,               &
+                  sums_wsus_ws_l, sums_wsvs_ws_l  
 
 !
@@ -287 +289 @@
           sums_wspts_ws_l = 0.0_wp
 
-          IF ( humidity  .OR.  passive_scalar )  THEN
+          IF ( humidity  )  THEN
              ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
              sums_wsqs_ws_l = 0.0_wp
+          ENDIF
+          
+          IF ( passive_scalar )  THEN
+             ALLOCATE( sums_wsss_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+             sums_wsss_ws_l = 0.0_wp
           ENDIF
 
@@ -341 +348 @@
                        diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
 
-             IF ( humidity  .OR.  passive_scalar )  THEN
+             IF ( humidity )  THEN
                 ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1),            &
                           diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
@@ -347 +354 @@
                           diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
              ENDIF
+             
+             IF ( passive_scalar )  THEN
+                ALLOCATE( flux_s_s(nzb+1:nzt,0:threads_per_task-1),            &
+                          diss_s_s(nzb+1:nzt,0:threads_per_task-1) )
+                ALLOCATE( flux_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
+                          diss_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
+             ENDIF             
 
              IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -837 +851 @@
            ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsnrs_ws_l,&
                   sums_wspts_ws_l, sums_wsqrs_ws_l, sums_wsqs_ws_l,            &
-                  sums_wssas_ws_l,  sums_wsus_ws_l, sums_wsvs_ws_l  
+                  sums_wsss_ws_l, sums_wssas_ws_l,  sums_wsus_ws_l,            &
+                  sums_wsvs_ws_l  
 
        IMPLICIT NONE
@@ -854 +869 @@
        IF ( ws_scheme_sca )  THEN
           sums_wspts_ws_l = 0.0_wp
-          IF ( humidity  .OR.  passive_scalar )  sums_wsqs_ws_l = 0.0_wp
+          IF ( humidity       )  sums_wsqs_ws_l = 0.0_wp
+          IF ( passive_scalar )  sums_wsss_ws_l = 0.0_wp
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              sums_wsqrs_ws_l = 0.0_wp
@@ -898 +914 @@
        USE statistics,                                                         &
            ONLY:  sums_wsnrs_ws_l, sums_wspts_ws_l, sums_wsqrs_ws_l,           &
-                  sums_wsqs_ws_l, sums_wssas_ws_l, weight_substep
+                  sums_wsqs_ws_l, sums_wssas_ws_l, sums_wsss_ws_l,             &
+                  weight_substep
 
        IMPLICIT NONE
@@ -1665 +1682 @@
                                * weight_substep(intermediate_timestep_count)
              ENDDO
+             
+          CASE ( 's' )
+          
+             DO  k = nzb, nzt
+                sums_wsss_ws_l(k,tn)  = sums_wsss_ws_l(k,tn) +                &
+                                      ( flux_t(k) + diss_t(k) )               &
+                               * weight_substep(intermediate_timestep_count)
+             ENDDO
 
          END SELECT
@@ -3120 +3145 @@
        USE statistics,                                                        &
            ONLY:  sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,           &
-                  sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep
+                  sums_wsqrs_ws_l, sums_wsnrs_ws_l, sums_wsss_ws_l,           &
+                  weight_substep
 
        IMPLICIT NONE
@@ -3866 +3892 @@
                         *   weight_substep(intermediate_timestep_count)
                     ENDDO
+                 CASE ( 's' )
+                    DO  k = nzb, nzt
+                       sums_wsss_ws_l(k,tn) = sums_wsss_ws_l(k,tn)          &
+                        + ( flux_t(k) + diss_t(k) )                           &
+                        *   weight_substep(intermediate_timestep_count)
+                    ENDDO   
+                                   
 
               END SELECT

palm/trunk/SOURCE/average_3d_data.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Treat humidity and passive scalar separately
 ! 
 ! Former revisions:
@@ -498 +498 @@
              ENDDO
 
+          CASE ( 'ssws*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   ssws_av(j,i) = ssws_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
           CASE ( 't*' )
              DO  i = nxlg, nxrg

palm/trunk/SOURCE/boundary_conds.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Treat humidity and passive scalar separately
 ! 
 ! Former revisions:
@@ -135 +135 @@
     USE arrays_3d,                                                             &
         ONLY:  c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l,  &
-               dzu, e_p, nr_p, pt, pt_p, q, q_p, qr_p, sa, sa_p,               &
+               dzu, e_p, nr_p, pt, pt_p, q, q_p, qr_p, s, s_p, sa, sa_p,       &
                u, ug, u_init, u_m_l, u_m_n, u_m_r, u_m_s, u_p,                 &
                v, vg, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p,                 &
-               w, w_p, w_m_l, w_m_n, w_m_r, w_m_s,&
-               pt_init
+               w, w_p, w_m_l, w_m_n, w_m_r, w_m_s, pt_init
 
     USE control_parameters,                                                    &
-        ONLY:  bc_pt_t_val, bc_q_t_val, constant_diffusion,                    &
+        ONLY:  bc_pt_t_val, bc_q_t_val, bc_s_t_val, constant_diffusion,        &
                cloud_physics, dt_3d, humidity,                                 &
-               ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_sa_t, ibc_uv_b,       &
-               ibc_uv_t, inflow_l, inflow_n, inflow_r, inflow_s,               &
-               intermediate_timestep_count, large_scale_forcing,               &
+               ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, ibc_s_t,         &
+               ibc_sa_t, ibc_uv_b, ibc_uv_t, inflow_l, inflow_n, inflow_r,     &
+               inflow_s, intermediate_timestep_count, large_scale_forcing,     &
                microphysics_seifert, nest_domain, nest_bound_l, nest_bound_s,  &
                nudging, ocean, outflow_l, outflow_n, outflow_r, outflow_s,     &
@@ -302 +301 @@
 
 !
-!-- Boundary conditions for total water content or scalar,
+!-- Boundary conditions for total water content,
 !-- bottom and top boundary (see also temperature)
-    IF ( humidity  .OR.  passive_scalar )  THEN
+    IF ( humidity )  THEN
 !
 !--    Surface conditions for constant_humidity_flux
@@ -344 +343 @@
        ENDIF
     ENDIF
+!
+!-- Boundary conditions for scalar,
+!-- bottom and top boundary (see also temperature)
+    IF ( passive_scalar )  THEN
+!
+!--    Surface conditions for constant_humidity_flux
+       IF ( ibc_s_b == 0 ) THEN
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+                s_p(nzb_s_inner(j,i),j,i) = s(nzb_s_inner(j,i),j,i)
+             ENDDO
+          ENDDO
+       ELSE
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+                s_p(nzb_s_inner(j,i),j,i) = s_p(nzb_s_inner(j,i)+1,j,i)
+             ENDDO
+          ENDDO
+       ENDIF
+!
+!--    Top boundary
+       IF ( ibc_s_t == 0 ) THEN
+          s_p(nzt+1,:,:) = s(nzt+1,:,:)
+       ELSEIF ( ibc_s_t == 1 ) THEN
+          s_p(nzt+1,:,:) = s_p(nzt,:,:)   + bc_s_t_val * dzu(nzt+1)
+       ENDIF
+
+    ENDIF    
 !
 !-- In case of inflow or nest boundary at the south boundary the boundary for v
@@ -381 +408 @@
        pt_p(:,nys-1,:)     = pt_p(:,nys,:)
        IF ( .NOT. constant_diffusion     )  e_p(:,nys-1,:) = e_p(:,nys,:)
-       IF ( humidity  .OR.  passive_scalar )  THEN
+       IF ( humidity )  THEN
           q_p(:,nys-1,:) = q_p(:,nys,:)
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -388 +415 @@
           ENDIF
        ENDIF
+       IF ( passive_scalar )  s_p(:,nys-1,:) = s_p(:,nys,:)
     ELSEIF ( outflow_n )  THEN
        pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
        IF ( .NOT. constant_diffusion     )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
-       IF ( humidity  .OR.  passive_scalar )  THEN
+       IF ( humidity )  THEN
           q_p(:,nyn+1,:) = q_p(:,nyn,:)
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -398 +426 @@
           ENDIF
        ENDIF
+       IF ( passive_scalar )  s_p(:,nyn+1,:) = s_p(:,nyn,:)
     ELSEIF ( outflow_l )  THEN
        pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
        IF ( .NOT. constant_diffusion     )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
-       IF ( humidity  .OR.  passive_scalar )  THEN
+       IF ( humidity )  THEN
           q_p(:,:,nxl-1) = q_p(:,:,nxl)
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -408 +437 @@
           ENDIF
        ENDIF
+       IF ( passive_scalar )  s_p(:,:,nxl-1) = s_p(:,:,nxl)
     ELSEIF ( outflow_r )  THEN
        pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
        IF ( .NOT. constant_diffusion     )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
-       IF ( humidity .OR. passive_scalar )  THEN
+       IF ( humidity )  THEN
           q_p(:,:,nxr+1) = q_p(:,:,nxr)
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -418 +448 @@
           ENDIF
        ENDIF
+       IF ( passive_scalar )  s_p(:,:,nxr+1) = s_p(:,:,nxr)
     ENDIF
 

palm/trunk/SOURCE/check_parameters.f90

@@ -19 +19 @@
 ! 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 +1057 @@
     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 +1092 @@
 !--    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 +1278 @@
 !--    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 +1698 @@
 
 !
-!-- 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 +1725 @@
        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 +3656 @@
                         '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 +3662 @@
                         '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 +3676 @@
     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 +3755 @@
 
     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

palm/trunk/SOURCE/data_output_2d.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-!
+! Scalar surface flux added
+! Rename INTEGER variable s into s_ind, as s is already assigned to scalar
 !
 ! Former revisions:
@@ -138 +139 @@
     USE arrays_3d,                                                             &
         ONLY:  dzw, e, nr, ol, p, pt, precipitation_amount, precipitation_rate,&
-               prr,q, qc, ql, ql_c, ql_v, ql_vp, qr, qsws, rho, sa, shf, tend, &
-               ts, u, us, v, vpt, w, z0, z0h, z0q, zu, zw
+               prr,q, qc, ql, ql_c, ql_v, ql_vp, qr, qsws, rho, s, sa, shf,    &
+               ssws, tend, ts, u, us, v, vpt, w, z0, z0h, z0q, zu, zw
         
     USE averaging
@@ -224 +225 @@
     INTEGER(iwp) ::  nzt_do    !< upper limit of the data field (usually nzt+1)
     INTEGER(iwp) ::  psi       !< 
-    INTEGER(iwp) ::  s         !< 
+    INTEGER(iwp) ::  s_ind     !< 
     INTEGER(iwp) ::  sender    !< 
     INTEGER(iwp) ::  ind(4)    !< 
@@ -268 +269 @@
 
        CASE ( 'xy' )
-          s = 1
+          s_ind = 1
           ALLOCATE( level_z(nzb:nzt+1), local_2d(nxlg:nxrg,nysg:nyng) )
 
           IF ( netcdf_data_format > 4 )  THEN
              ns = 1
-             DO WHILE ( section(ns,s) /= -9999  .AND.  ns <= 100 )
+             DO WHILE ( section(ns,s_ind) /= -9999  .AND.  ns <= 100 )
                 ns = ns + 1
              ENDDO
@@ -298 +299 @@
 
        CASE ( 'xz' )
-          s = 2
+          s_ind = 2
           ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) )
 
           IF ( netcdf_data_format > 4 )  THEN
              ns = 1
-             DO WHILE ( section(ns,s) /= -9999  .AND.  ns <= 100 )
+             DO WHILE ( section(ns,s_ind) /= -9999  .AND.  ns <= 100 )
                 ns = ns + 1
              ENDDO
@@ -329 +330 @@
 
        CASE ( 'yz' )
-          s = 3
+          s_ind = 3
           ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) )
 
           IF ( netcdf_data_format > 4 )  THEN
              ns = 1
-             DO WHILE ( section(ns,s) /= -9999  .AND.  ns <= 100 )
+             DO WHILE ( section(ns,s_ind) /= -9999  .AND.  ns <= 100 )
                 ns = ns + 1
              ENDDO
@@ -1070 +1071 @@
              CASE ( 's_xy', 's_xz', 's_yz' )
                 IF ( av == 0 )  THEN
-                   to_be_resorted => q
+                   to_be_resorted => s
                 ELSE
                    to_be_resorted => s_av
@@ -1117 +1118 @@
                 two_d = .TRUE.
                 level_z(nzb+1) = zu(nzb+1)
+                
+             CASE ( 'ssws*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  ssws(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng 
+                         local_pf(i,j,nzb+1) =  ssws_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)                
 
              CASE ( 't*_xy' )        ! 2d-array
@@ -1303 +1322 @@
 !--       section mode chosen.
           is = 1
-   loop1: DO WHILE ( section(is,s) /= -9999  .OR.  two_d )
+   loop1: DO WHILE ( section(is,s_ind) /= -9999  .OR.  two_d )
 
              SELECT CASE ( mode )
@@ -1313 +1332 @@
                       layer_xy = nzb+1
                    ELSE
-                      layer_xy = section(is,s)
+                      layer_xy = section(is,s_ind)
                    ENDIF
 
@@ -1347 +1366 @@
 !
 !--                If required, carry out averaging along z
-                   IF ( section(is,s) == -1  .AND.  .NOT. two_d )  THEN
+                   IF ( section(is,s_ind) == -1  .AND.  .NOT. two_d )  THEN
 
                       local_2d = 0.0_wp
@@ -1537 +1556 @@
 !
 !--                If required, carry out averaging along y
-                   IF ( section(is,s) == -1 )  THEN
+                   IF ( section(is,s_ind) == -1 )  THEN
 
                       ALLOCATE( local_2d_l(nxlg:nxrg,nzb_do:nzt_do) )
@@ -1570 +1589 @@
 !--                   Just store the respective section on the local array
 !--                   (but only if it is available on this PE!)
-                      IF ( section(is,s) >= nys  .AND.  section(is,s) <= nyn ) &
+                      IF ( section(is,s_ind) >= nys  .AND.  section(is,s_ind) <= nyn ) &
                       THEN
-                         local_2d = local_pf(:,section(is,s),nzb_do:nzt_do)
+                         local_2d = local_pf(:,section(is,s_ind),nzb_do:nzt_do)
                       ENDIF
 
@@ -1584 +1603 @@
 !--                   sections reside. Cross sections averaged along y are
 !--                   output on the respective first PE along y (myidy=0).
-                      IF ( ( section(is,s) >= nys  .AND.                       &
-                             section(is,s) <= nyn )  .OR.                      &
-                           ( section(is,s) == -1  .AND.  myidy == 0 ) )  THEN
+                      IF ( ( section(is,s_ind) >= nys  .AND.                   &
+                             section(is,s_ind) <= nyn )  .OR.                  &
+                           ( section(is,s_ind) == -1  .AND.  myidy == 0 ) )  THEN
 #if defined( __netcdf )
 !
@@ -1615 +1634 @@
                          DO  i = 0, io_blocks-1
                             IF ( i == io_group )  THEN
-                               IF ( ( section(is,s) >= nys  .AND.              &
-                                      section(is,s) <= nyn )  .OR.             &
-                                    ( section(is,s) == -1  .AND.               &
+                               IF ( ( section(is,s_ind) >= nys  .AND.          &
+                                      section(is,s_ind) <= nyn )  .OR.         &
+                                    ( section(is,s_ind) == -1  .AND.           &
                                       nys-1 == -1 ) )                          &
                                THEN
@@ -1643 +1662 @@
 !
 !--                         Local array can be relocated directly.
-                            IF ( ( section(is,s) >= nys  .AND.                 &
-                                   section(is,s) <= nyn )  .OR.                &
-                                 ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
-                            THEN
+                            IF ( ( section(is,s_ind) >= nys  .AND.              &
+                                   section(is,s_ind) <= nyn )  .OR.             &
+                                 ( section(is,s_ind) == -1  .AND.               &
+                                   nys-1 == -1 ) )  THEN
                                total_2d(nxlg:nxrg,nzb_do:nzt_do) = local_2d
                             ENDIF
@@ -1691 +1710 @@
 !--                         If the cross section resides on the PE, send the
 !--                         local index limits, otherwise send -9999 to PE0.
-                            IF ( ( section(is,s) >= nys  .AND.                 &
-                                   section(is,s) <= nyn )  .OR.                &
-                                 ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
+                            IF ( ( section(is,s_ind) >= nys  .AND.              &
+                                   section(is,s_ind) <= nyn )  .OR.             &
+                                 ( section(is,s_ind) == -1  .AND.  nys-1 == -1 ) ) &
                             THEN
                                ind(1) = nxlg; ind(2) = nxrg
@@ -1756 +1775 @@
 !
 !--                If required, carry out averaging along x
-                   IF ( section(is,s) == -1 )  THEN
+                   IF ( section(is,s_ind) == -1 )  THEN
 
                       ALLOCATE( local_2d_l(nysg:nyng,nzb_do:nzt_do) )
@@ -1789 +1808 @@
 !--                   Just store the respective section on the local array
 !--                   (but only if it is available on this PE!)
-                      IF ( section(is,s) >= nxl  .AND.  section(is,s) <= nxr ) &
+                      IF ( section(is,s_ind) >= nxl  .AND.  section(is,s_ind) <= nxr ) &
                       THEN
-                         local_2d = local_pf(section(is,s),:,nzb_do:nzt_do)
+                         local_2d = local_pf(section(is,s_ind),:,nzb_do:nzt_do)
                       ENDIF
 
@@ -1803 +1822 @@
 !--                   sections reside. Cross sections averaged along x are
 !--                   output on the respective first PE along x (myidx=0).
-                      IF ( ( section(is,s) >= nxl  .AND.                       &
-                             section(is,s) <= nxr )  .OR.                      &
-                           ( section(is,s) == -1  .AND.  myidx == 0 ) )  THEN
+                      IF ( ( section(is,s_ind) >= nxl  .AND.                       &
+                             section(is,s_ind) <= nxr )  .OR.                      &
+                           ( section(is,s_ind) == -1  .AND.  myidx == 0 ) )  THEN
 #if defined( __netcdf )
 !
@@ -1834 +1853 @@
                          DO  i = 0, io_blocks-1
                             IF ( i == io_group )  THEN
-                               IF ( ( section(is,s) >= nxl  .AND.              &
-                                      section(is,s) <= nxr )  .OR.             &
-                                    ( section(is,s) == -1  .AND.               &
+                               IF ( ( section(is,s_ind) >= nxl  .AND.          &
+                                      section(is,s_ind) <= nxr )  .OR.         &
+                                    ( section(is,s_ind) == -1  .AND.           &
                                       nxl-1 == -1 ) )                          &
                                THEN
@@ -1862 +1881 @@
 !
 !--                         Local array can be relocated directly.
-                            IF ( ( section(is,s) >= nxl  .AND.                 &
-                                   section(is,s) <= nxr )   .OR.               &
-                                 ( section(is,s) == -1  .AND.  nxl-1 == -1 ) ) &
+                            IF ( ( section(is,s_ind) >= nxl  .AND.             &
+                                   section(is,s_ind) <= nxr )   .OR.           &
+                                 ( section(is,s_ind) == -1  .AND.  nxl-1 == -1 ) ) &
                             THEN
                                total_2d(nysg:nyng,nzb_do:nzt_do) = local_2d
@@ -1910 +1929 @@
 !--                         If the cross section resides on the PE, send the
 !--                         local index limits, otherwise send -9999 to PE0.
-                            IF ( ( section(is,s) >= nxl  .AND.                 &
-                                   section(is,s) <= nxr )  .OR.                &
-                                 ( section(is,s) == -1  .AND.  nxl-1 == -1 ) ) &
+                            IF ( ( section(is,s_ind) >= nxl  .AND.              &
+                                   section(is,s_ind) <= nxr )  .OR.             &
+                                 ( section(is,s_ind) == -1  .AND.  nxl-1 == -1 ) ) &
                             THEN
                                ind(1) = nysg; ind(2) = nyng
@@ -2148 +2167 @@
 !
 !-- Close plot output file.
-    file_id = 20 + s
+    file_id = 20 + s_ind
 
     IF ( data_output_2d_on_each_pe )  THEN

palm/trunk/SOURCE/data_output_3d.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-!
+! Scalar surface flux added
 !
 ! Former revisions:
@@ -119 +119 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  e, nr, p, pt, prr, q, qc, ql, ql_c, ql_v, qr, rho, sa, tend, u, &
-               v, vpt, w
+        ONLY:  e, nr, p, pt, prr, q, qc, ql, ql_c, ql_v, qr, rho, s, sa, tend, &
+               u, v, vpt, w
         
     USE averaging
@@ -566 +566 @@
           CASE ( 's' )
              IF ( av == 0 )  THEN
-                to_be_resorted => q
+                to_be_resorted => s
              ELSE
                 to_be_resorted => s_av

palm/trunk/SOURCE/data_output_dvrp.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Separate humidity and passive scalar
 ! 
 ! Former revisions:
@@ -112 +112 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  p, pt, q, ql, ts, u, us, v, w, zu
+        ONLY:  p, pt, q, ql, s, ts, u, us, v, w, zu
         
     USE cloud_parameters,                                                      &
@@ -336 +336 @@
 
                 CASE ( 'q', 'q_xy', 'q_xz', 'q_yz' )
-                   IF ( humidity  .OR.  passive_scalar )  THEN
+                   IF ( humidity )  THEN
                       DO  i = nxl_dvrp, nxr_dvrp+1
                          DO  j = nys_dvrp, nyn_dvrp+1
@@ -345 +345 @@
                       ENDDO            
                    ELSE                    
-                      message_string = 'if humidity/passive_scalar = '    //   & 
+                      message_string = 'if humidity = '    //                  & 
                             'FALSE output of ' // TRIM( output_variable ) //   &
                             'is not provided' 
@@ -366 +366 @@
                                   'is not provided' 
                       CALL message( 'data_output_dvrp', 'PA0184',&
+                                                                 0, 0, 0, 6, 0 )
+                   ENDIF
+
+                CASE ( 's', 's_xy', 's_xz', 's_yz' )
+                   IF ( passive_scalar )  THEN
+                      DO  i = nxl_dvrp, nxr_dvrp+1
+                         DO  j = nys_dvrp, nyn_dvrp+1
+                            DO  k = nzb, nz_do3d
+                               local_pf(i,j,k) = s(k,j,i)
+                            ENDDO
+                         ENDDO
+                      ENDDO            
+                   ELSE                    
+                      message_string = 'if passive_scalar = '    //            & 
+                            'FALSE output of ' // TRIM( output_variable ) //   &
+                            'is not provided' 
+                      CALL message( 'data_output_dvrp', 'PA0183',&
                                                                  0, 0, 0, 6, 0 )
                    ENDIF

palm/trunk/SOURCE/data_output_mask.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Separate humidity and passive scalar
 !
 ! Former revisions:
@@ -88 +88 @@
 #if defined( __netcdf )
     USE arrays_3d,                                                             &
-        ONLY:  e, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho, sa, tend, u,      &
+        ONLY:  e, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho, s, sa, tend, u,   &
                v, vpt, w
     
@@ -487 +487 @@
           CASE ( 's' )
              IF ( av == 0 )  THEN
-                to_be_resorted => q
+                to_be_resorted => s
              ELSE
                 to_be_resorted => s_av

palm/trunk/SOURCE/data_output_spectra.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Additional default spectra for passive scalar
 ! 
 ! Former revisions:
@@ -172 +172 @@
                 pr = 5
 
+             CASE ( 's' )
+                pr = 6
+
              CASE DEFAULT
 !

palm/trunk/SOURCE/flow_statistics.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Separate humidity and passive scalar
 !
 ! Former revisions:
@@ -197 +197 @@
     USE arrays_3d,                                                             &
         ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, ol, p, prho, prr, pt, q, qc, ql,&
-               qr, qs, qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt,  &
-               td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us,& 
-               usws, uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
+               qr, qs, qsws, qswst, rho, s, sa, ss, ssws, sswst, saswsb,       &
+               saswst, shf, td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q,        &
+               time_vert, ts, tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, &
+               vswst, w, w_subs, zw
         
     USE cloud_parameters,                                                      &
@@ -360 +361 @@
 
           DO  i = 0, threads_per_task-1
-             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
-             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
-             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
-                                                   sums_wsqs_ws_l(:,i) !w*q*
+             sums_l(:,17,i)                        = sums_wspts_ws_l(:,i) ! w*pt* 
+             IF ( ocean          ) sums_l(:,66,i)  = sums_wssas_ws_l(:,i) ! w*sa*
+             IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)  ! w*q*
+             IF ( passive_scalar ) sums_l(:,116,i) = sums_wsss_ws_l(:,i)  ! w*s*
           ENDDO
 
@@ -440 +441 @@
              DO  j =  nys, nyn
                 DO  k = nzb_s_inner(j,i), nzt+1
-                   sums_l(k,41,tn)  = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
+                   sums_l(k,117,tn)  = sums_l(k,117,tn) + s(k,j,i) * rmask(j,i,sr)
                 ENDDO
              ENDDO
@@ -465 +466 @@
              ENDIF
              IF ( passive_scalar )  THEN
-                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
+                sums_l(:,117,0) = sums_l(:,117,0) + sums_l(:,117,i)
              ENDIF
           ENDDO
@@ -506 +507 @@
        IF ( passive_scalar )  THEN
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
+          CALL MPI_ALLREDUCE( sums_l(nzb,117,0), sums(nzb,117), nzt+2-nzb,       &
                               MPI_REAL, MPI_SUM, comm2d, ierr )
        ENDIF
@@ -522 +523 @@
           ENDIF
        ENDIF
-       IF ( passive_scalar )  sums(:,41) = sums_l(:,41,0)
+       IF ( passive_scalar )  sums(:,117) = sums_l(:,117,0)
 #endif
 
@@ -560 +561 @@
 !
 !--    Passive scalar
-       IF ( passive_scalar )  hom(:,1,41,sr) = sums(:,41) /                    &
-            ngp_2dh_s_inner(:,sr)                    ! s (q)
+       IF ( passive_scalar )  hom(:,1,117,sr) = sums(:,117) /                  &
+            ngp_2dh_s_inner(:,sr)                    ! s 
 
 !
@@ -598 +599 @@
                                   ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
                 ENDIF
-
+                IF ( passive_scalar )  THEN
+                   sums_l(k,118,tn) = sums_l(k,118,tn) + &
+                                  ( s(k,j,i)-hom(k,1,117,sr) )**2 * rmask(j,i,sr)
+                ENDIF
 !
 !--             Higher moments
@@ -627 +631 @@
                 sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
                                             qs(j,i)   * rmask(j,i,sr)
+             ENDIF
+             IF ( passive_scalar )  THEN
+                sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
+                                            ss(j,i)   * rmask(j,i,sr)
              ENDIF
           ENDDO
@@ -738 +746 @@
 !--             Passive scalar flux
                 IF ( passive_scalar )  THEN
-                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
+                   sums_l(k,119,tn) = sums_l(k,119,tn)                         &
                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
-                                               * ( q(k+1,j,i) - q(k,j,i) )     &
+                                               * ( s(k+1,j,i) - s(k,j,i) )     &
                                                * ddzu(k+1) * rmask(j,i,sr)
                 ENDIF
@@ -784 +792 @@
                 ENDIF
                 IF ( passive_scalar )  THEN
-                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
-                                       qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
+                   sums_l(nzb,119,tn) = sums_l(nzb,119,tn) +                     &
+                                       ssws(j,i) * rmask(j,i,sr) ! w"s"
                 ENDIF
              ENDIF
@@ -863 +871 @@
                 ENDIF
                 IF ( passive_scalar )  THEN
-                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
-                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
+                   sums_l(nzt,119,tn) = sums_l(nzt,119,tn) + &
+                                       sswst(j,i) * rmask(j,i,sr) ! w"s"
                 ENDIF
              ENDIF
@@ -953 +961 @@
                 IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
                      .OR. sr /= 0 ) )  THEN
-                   pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +              &
-                                    q(k+1,j,i) - hom(k+1,1,41,sr) )
-                   sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *        &
+                   pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +             &
+                                    s(k+1,j,i) - hom(k+1,1,117,sr) )
+                   sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *      &
                                                        rmask(j,i,sr)
                 ENDIF
@@ -1011 +1019 @@
                                       q(k+1,j,i) - hom(k+1,1,41,sr) )
                      sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *      &
+                                       rmask(j,i,sr)
+                  ENDIF
+                  IF ( passive_scalar )  THEN
+                     pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +            &
+                                      s(k+1,j,i) - hom(k+1,1,117,sr) )
+                     sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *      &
                                        rmask(j,i,sr)
                   ENDIF
@@ -1279 +1293 @@
           sums(k,81:88)         = sums(k,81:88)         / ngp_2dh(sr)
           sums(k,89:114)        = sums(k,89:114)        / ngp_2dh(sr)
+          sums(k,116)           = sums(k,116)           / ngp_2dh(sr)
+          sums(k,119)           = sums(k,119)           / ngp_2dh(sr)
           IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
              sums(k,8:11)          = sums(k,8:11)          / ngp_2dh_s_inner(k,sr)
@@ -1287 +1303 @@
              sums(k,64)            = sums(k,64)            / ngp_2dh_s_inner(k,sr)
              sums(k,70:80)         = sums(k,70:80)         / ngp_2dh_s_inner(k,sr)
-             sums(k,115:pr_palm-2) = sums(k,115:pr_palm-2) / ngp_2dh_s_inner(k,sr)
+             sums(k,118)           = sums(k,118)           / ngp_2dh_s_inner(k,sr)
+             sums(k,120:pr_palm-2) = sums(k,120:pr_palm-2) / ngp_2dh_s_inner(k,sr)
           ENDIF
        ENDDO
@@ -1298 +1315 @@
                                     ngp_2dh(sr)
        sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
+                                    ngp_2dh(sr)
+       sums(nzb+13,pr_palm)       = sums(nzb+13,pr_palm)       / &    ! ss
                                     ngp_2dh(sr)
 !--    eges, e*
@@ -1443 +1462 @@
        hom(:,1,114,sr) = sums(:,114)            !: L
 
+       IF ( passive_scalar )  THEN
+          hom(:,1,119,sr) = sums(:,119)     ! w"s"
+          hom(:,1,116,sr) = sums(:,116)     ! w*s*
+          hom(:,1,120,sr) = sums(:,119) + sums(:,116)    ! ws 
+       ENDIF
+
        hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
                                        ! u*, w'u', w'v', t* (in last profile)
@@ -1589 +1614 @@
        ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
 
+       IF ( passive_scalar )  THEN
+          ts_value(24,sr) = hom(nzb+13,1,119,sr)       ! w"s" ( to do ! )
+          ts_value(25,sr) = hom(nzb+13,1,pr_palm,sr)   ! s*
+       ENDIF
+
 !
 !--    Collect land surface model timeseries
@@ -1660 +1690 @@
     USE arrays_3d,                                                             &
         ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
-               qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q,&
-               td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
-               uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
+               qsws, qswst, rho, s, sa, saswsb, saswst, shf, ss, ssws, sswst,  &
+               td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts,      &
+               tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, vswst, w,      &
+               w_subs, zw
                  
         
@@ -1833 +1864 @@
              sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
              IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
-             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
-                                                   sums_wsqs_ws_l(:,i) !w*q*
+             IF ( humidity       )  sums_l(:,49,i)  = sums_wsqs_ws_l(:,i) !w*q*
+             IF ( passive_scalar )  sums_l(:,116,i) = sums_wsss_ws_l(:,i) !w*s*
           ENDDO
 
@@ -1938 +1969 @@
        IF ( passive_scalar )  THEN
           !$OMP DO
-          !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l ) create( s1 )
+          !$acc parallel loop gang present( s, rflags_invers, rmask, sums_l ) create( s1 )
           DO  k = nzb, nzt+1
              s1 = 0
@@ -1944 +1975 @@
              DO  i = nxl, nxr
                 DO  j =  nys, nyn
-                   s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
-                ENDDO
-             ENDDO
-             sums_l(k,41,tn) = s1
+                   s1 = s1 + s(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                ENDDO
+             ENDDO
+             sums_l(k,117,tn) = s1
           ENDDO
           !$acc end parallel loop
@@ -1981 +2012 @@
              IF ( passive_scalar )  THEN
                 !$acc parallel present( sums_l )
-                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
+                sums_l(:,117,0) = sums_l(:,117,0) + sums_l(:,117,i)
                 !$acc end parallel
              ENDIF
@@ -2024 +2055 @@
        IF ( passive_scalar )  THEN
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
+          CALL MPI_ALLREDUCE( sums_l(nzb,117,0), sums(nzb,117), nzt+2-nzb,     &
                               MPI_REAL, MPI_SUM, comm2d, ierr )
        ENDIF
@@ -2053 +2084 @@
        IF ( passive_scalar )  THEN
           !$acc parallel present( sums, sums_l )
-          sums(:,41) = sums_l(:,41,0)
+          sums(:,117) = sums_l(:,117,0)
           !$acc end parallel
        ENDIF
@@ -2102 +2133 @@
        IF ( passive_scalar )  THEN
           !$acc parallel present( hom, ngp_2dh_s_inner, sums )
-          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
-          hom(:,1,41,sr) = sums(:,41)                ! s (q)
+          sums(:,117)     = sums(:,117) / ngp_2dh_s_inner(:,sr)
+          hom(:,1,117,sr) = sums(:,117)                ! s
           !$acc end parallel
        ENDIF
@@ -2230 +2261 @@
           ENDDO
           sums_l(nzb+12,pr_palm,tn) = s1
+          !$acc end parallel
+       ENDIF
+
+       IF ( passive_scalar )  THEN
+          !$acc parallel present( ss, rmask, sums_l ) create( s1 )
+          s1 = 0
+          !$acc loop vector collapse( 2 ) reduction( +: s1 )
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+                s1 = s1 + ss(j,i) * rmask(j,i,sr)
+             ENDDO
+          ENDDO
+          sums_l(nzb+13,pr_palm,tn) = s1
           !$acc end parallel
        ENDIF
@@ -2411 +2455 @@
        IF ( passive_scalar )  THEN
 
-          !$acc parallel loop gang present( ddzu, kh, q, rflags_invers, rmask, sums_l ) create( s1 )
+          !$acc parallel loop gang present( ddzu, kh, s, rflags_invers, rmask, sums_l ) create( s1 )
           DO  k = nzb, nzt_diff
              s1 = 0
@@ -2418 +2462 @@
                 DO  j = nys, nyn
                    s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
-                                    * ( q(k+1,j,i) - q(k,j,i) )                &
+                                    * ( s(k+1,j,i) - s(k,j,i) )                &
                                     * ddzu(k+1) * rmask(j,i,sr)                &
                                     * rflags_invers(j,i,k+1)
                 ENDDO
              ENDDO
-             sums_l(k,48,tn) = s1
+             sums_l(k,119,tn) = s1
           ENDDO
           !$acc end parallel loop
@@ -2532 +2576 @@
 
              !$OMP DO
-             !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
+             !$acc parallel present( ssws, rmask, sums_l ) create( s1 )
              s1 = 0
              !$acc loop vector collapse( 2 ) reduction( +: s1 )
              DO  i = nxl, nxr
                 DO  j =  nys, nyn
-                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
-                ENDDO
-             ENDDO
-             sums_l(nzb,48,tn) = s1
+                   s1 = s1 + ssws(j,i) * rmask(j,i,sr)  ! w"s"
+                ENDDO
+             ENDDO
+             sums_l(nzb,119,tn) = s1
              !$acc end parallel
 
@@ -2651 +2695 @@
 
              !$OMP DO
-             !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
+             !$acc parallel present( sswst, rmask, sums_l ) create( s1 )
              s1 = 0
              !$acc loop vector collapse( 2 ) reduction( +: s1 )
              DO  i = nxl, nxr
                 DO  j =  nys, nyn
-                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
-                ENDDO
-             ENDDO
-             sums_l(nzt,48,tn) = s1
+                   s1 = s1 + sswst(j,i) * rmask(j,i,sr) ! w"s"
+                ENDDO
+             ENDDO
+             sums_l(nzt,119,tn) = s1
              !$acc end parallel
 
@@ -2890 +2934 @@
        IF ( passive_scalar  .AND.  ( .NOT. ws_scheme_sca  .OR.  sr /= 0 ) )  THEN
 
-          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
+          !$acc parallel loop gang present( hom, s, rflags_invers, rmask, sums_l, w ) create( s1 )
           DO  k = nzb, nzt_diff
              s1 = 0
@@ -2896 +2940 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +          &
-                                        q(k+1,j,i) - hom(k+1,1,41,sr) )        &
+                   s1 = s1 + 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +          &
+                                        s(k+1,j,i) - hom(k+1,1,117,sr) )        &
                                     * w(k,j,i) * rmask(j,i,sr)                 &
                                     * rflags_invers(j,i,k+1)
@@ -2981 +3025 @@
                 ENDDO
                 sums_l(k,49,tn) = s1
+             ENDDO
+             !$acc end parallel loop
+
+          ENDIF
+
+          IF ( passive_scalar )  THEN
+
+             !$acc parallel loop gang present( hom, s, rflags_invers, rmask, sums_l, w ) create( s1 )
+             DO  k = nzb, nzt_diff
+                s1 = 0
+                !$acc loop vector collapse( 2 ) reduction( +: s1 )
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      s1 = s1 + 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +      &
+                                           s(k+1,j,i) - hom(k+1,1,117,sr) )    &
+                                       * w(k,j,i) * rmask(j,i,sr)              &
+                                       * rflags_invers(j,i,k+1)
+                   ENDDO
+                ENDDO
+                sums_l(k,116,tn) = s1
              ENDDO
              !$acc end parallel loop

palm/trunk/SOURCE/header.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Treat humidity and passive scalar separately.
+! Modify misleading information concerning humidity.
+! Bugfix, change unit for humidity flux.
 ! 
 ! Former revisions:
@@ -266 +268 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  pt_init, qsws, q_init, sa_init, shf, ug, vg, w_subs, zu, zw
+        ONLY:  pt_init, qsws, q_init, s_init, sa_init, shf, ug, vg, w_subs, zu,&
+               zw
         
     USE control_parameters
@@ -1022 +1025 @@
 
     IF ( passive_scalar )  THEN
-       IF ( ibc_q_b == 0 )  THEN
+       IF ( ibc_s_b == 0 )  THEN
           r_lower = 's(0)     = s_surface'
        ELSE
           r_lower = 's(0)     = s(1)'
        ENDIF
-       IF ( ibc_q_t == 0 )  THEN
+       IF ( ibc_s_t == 0 )  THEN
           r_upper =  's(nzt)   = s_top'
        ELSE
@@ -1052 +1055 @@
           ENDIF
        ENDIF
-       IF ( passive_scalar  .AND.  constant_waterflux )  THEN
-          WRITE ( io, 313 ) surface_waterflux
+       IF ( passive_scalar  .AND.  constant_scalarflux )  THEN
+          WRITE ( io, 313 ) surface_scalarflux
        ENDIF
     ENDIF
@@ -1070 +1073 @@
           WRITE ( io, 309 )  top_salinityflux
        ENDIF
-       IF ( humidity  .OR.  passive_scalar )  THEN
-          WRITE ( io, 315 )
-       ENDIF
+       IF ( humidity       )  WRITE ( io, 315 )
+       IF ( passive_scalar )  WRITE ( io, 315 )
     ENDIF
 
@@ -1083 +1085 @@
           WRITE ( io, 312 )
        ENDIF
-       IF ( passive_scalar  .AND.  .NOT. constant_waterflux )  THEN
+       IF ( passive_scalar  .AND.  .NOT. constant_scalarflux )  THEN
           WRITE ( io, 314 )
        ENDIF
@@ -1153 +1155 @@
 !-- Initial humidity profile
 !-- Building output strings, starting with surface humidity
-    IF ( humidity  .OR.  passive_scalar )  THEN
+    IF ( humidity )  THEN
        WRITE ( temperatures, '(E8.1)' )  q_surface
        gradients = '--------'
@@ -1181 +1183 @@
        ENDDO
 
-       IF ( humidity )  THEN
-          IF ( .NOT. nudging )  THEN
-             WRITE ( io, 421 )  TRIM( coordinates ), TRIM( temperatures ), &
-                                TRIM( gradients ), TRIM( slices )
-          ENDIF
-       ELSE
-          WRITE ( io, 422 )  TRIM( coordinates ), TRIM( temperatures ), &
+       IF ( .NOT. nudging )  THEN
+          WRITE ( io, 421 )  TRIM( coordinates ), TRIM( temperatures ),        &
                              TRIM( gradients ), TRIM( slices )
        ENDIF
     ENDIF
+!
+!-- Initial scalar profile
+!-- Building output strings, starting with surface humidity
+    IF ( passive_scalar )  THEN
+       WRITE ( temperatures, '(E8.1)' )  s_surface
+       gradients = '--------'
+       slices = '       0'
+       coordinates = '     0.0'
+       i = 1
+       DO  WHILE ( s_vertical_gradient_level_ind(i) /= -9999 )
+          
+          WRITE (coor_chr,'(E8.1,4X)')  s_init(q_vertical_gradient_level_ind(i))
+          temperatures = TRIM( temperatures ) // '  ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(E8.1,4X)')  s_vertical_gradient(i)
+          gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
+          
+          WRITE (coor_chr,'(I8,4X)')  s_vertical_gradient_level_ind(i)
+          slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
+          
+          WRITE (coor_chr,'(F8.1,4X)')  s_vertical_gradient_level(i)
+          coordinates = TRIM( coordinates ) // '  '  // TRIM( coor_chr )
+
+          IF ( i == 10 )  THEN
+             EXIT
+          ELSE
+             i = i + 1
+          ENDIF
+
+       ENDDO
+
+       WRITE ( io, 422 )  TRIM( coordinates ), TRIM( temperatures ),           &
+                          TRIM( gradients ), TRIM( slices )
+    ENDIF    
 
 !
@@ -1995 +2026 @@
 310 FORMAT (//'    1D-Model:'// &
              '       Rif value range:   ',F6.2,' <= rif <=',F6.2)
-311 FORMAT ('       Predefined constant humidity flux: ',E10.3,' m/s')
+311 FORMAT ('       Predefined constant humidity flux: ',E10.3,' kg/kg m/s')
 312 FORMAT ('       Predefined surface humidity')
 313 FORMAT ('       Predefined constant scalar flux: ',E10.3,' kg/(m**2 s)')
@@ -2158 +2189 @@
 430 FORMAT (//' Cloud physics quantities / methods:'/ &
               ' ----------------------------------'/)
-431 FORMAT ('    Humidity is treated as purely passive scalar (no condensati', &
-                 'on)')
+431 FORMAT ('    Humidity is considered, bu no condensation')
 432 FORMAT ('    Bulk scheme with liquid water potential temperature and'/ &
             '    total water content is used.'/ &

palm/trunk/SOURCE/inflow_turbulence.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Separate humidity and passive scalar
 ! 
 ! Former revisions:
@@ -71 +71 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  e, inflow_damping_factor, mean_inflow_profiles, pt, q, u, v, w
+        ONLY:  e, inflow_damping_factor, mean_inflow_profiles, pt, q, s, u, v, w
         
     USE control_parameters,                                                    &
@@ -98 +98 @@
     INTEGER(iwp) ::  prev     !< ID of sending PE for y-shift
 
-    REAL(wp), DIMENSION(nzb:nzt+1,6,nbgp)           ::                         &
+    REAL(wp), DIMENSION(nzb:nzt+1,7,nbgp)           ::                         &
        avpr               !< stores averaged profiles at recycling plane
-    REAL(wp), DIMENSION(nzb:nzt+1,6,nbgp)           ::                         &
+    REAL(wp), DIMENSION(nzb:nzt+1,7,nbgp)           ::                         &
        avpr_l             !< auxiliary variable to calculate avpr
-    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,6,nbgp) ::                         &
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,7,nbgp) ::                         &
        inflow_dist        !< turbulence signal of vars, added at inflow boundary
-    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,6,nbgp) ::                         &
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,7,nbgp) ::                         &
        local_inflow_dist  !< auxiliary variable for inflow_dist, used for yshift
 
@@ -112 +112 @@
 !-- Carry out spanwise averaging in the recycling plane
     avpr_l = 0.0_wp
-    ngp_pr = ( nzt - nzb + 2 ) * 6 * nbgp
+    ngp_pr = ( nzt - nzb + 2 ) * 7 * nbgp
     ngp_ifd = ngp_pr * ( nyn - nys + 1 + 2 * nbgp )
 
@@ -131 +131 @@
                 avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)
                 avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i)
-                IF ( humidity  .OR.  passive_scalar )                          &
+                IF ( humidity )                                                &
                    avpr_l(k,6,l) = avpr_l(k,6,l) + q(k,j,i)
+                IF ( passive_scalar )                                          &
+                   avpr_l(k,7,l) = avpr_l(k,7,l) + s(k,j,i)
 
              ENDDO
@@ -156 +158 @@
              avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)
              avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i)
-             IF ( humidity  .OR.  passive_scalar )                             &
+             IF ( humidity )                                                   &
                 avpr_l(k,6,l) = avpr_l(k,6,l) + q(k,j,i)
+             IF ( passive_scalar )                                             &
+                avpr_l(k,7,l) = avpr_l(k,7,l) + s(k,j,i)
 
           ENDDO
@@ -183 +187 @@
                 inflow_dist(k,j,4,l) = pt(k,j,i)  - avpr(k,4,l)
                 inflow_dist(k,j,5,l) = e(k,j,i)   - avpr(k,5,l)
-                IF ( humidity  .OR.  passive_scalar )                          &
+                IF ( humidity )                                                &
                    inflow_dist(k,j,6,l) = q(k,j,i) - avpr(k,6,l)
+                IF ( passive_scalar )                                          &
+                   inflow_dist(k,j,7,l) = s(k,j,i) - avpr(k,7,l)
             ENDDO
           ENDDO
@@ -201 +207 @@
              inflow_dist(k,j,4,l) = pt(k,j,i)  - avpr(k,4,l)
              inflow_dist(k,j,5,l) = e(k,j,i)   - avpr(k,5,l)
-             IF ( humidity  .OR.  passive_scalar )                             &
+             IF ( humidity )                                                   &
                 inflow_dist(k,j,6,l) = q(k,j,i) - avpr(k,6,l)
+             IF ( passive_scalar )                                             &
+                inflow_dist(k,j,7,l) = s(k,j,i) - avpr(k,7,l)
               
           ENDDO
@@ -276 +284 @@
              e(k,j,-nbgp:-1)  = MAX( e(k,j,-nbgp:-1), 0.0_wp )
 
-             IF ( humidity  .OR.  passive_scalar )                          &
+             IF ( humidity )                                                &
                 q(k,j,-nbgp:-1)  = mean_inflow_profiles(k,6) +              &
                         inflow_dist(k,j,6,1:nbgp) * inflow_damping_factor(k)
+             IF ( passive_scalar )                                          &
+                s(k,j,-nbgp:-1)  = mean_inflow_profiles(k,7) +              &
+                        inflow_dist(k,j,7,1:nbgp) * inflow_damping_factor(k)
 
           ENDDO

palm/trunk/SOURCE/init_1d_model.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Remove passive_scalar from IF-statements, as 1D-scalar profile is effectively
+! not used.
+! Formatting adjustment
 !
 ! Former revisions:
@@ -108 +110 @@
     USE control_parameters,                                                    &
         ONLY:  constant_diffusion, constant_flux_layer, f, humidity, kappa,    &
-               km_constant, mixing_length_1d, passive_scalar, prandtl_number,  &
+               km_constant, mixing_length_1d, prandtl_number,                  &
                roughness_length, simulated_time_chr, z0h_factor
 
@@ -194 +196 @@
     z01d  = roughness_length
     z0h1d = z0h_factor * z01d 
-    IF ( humidity .OR. passive_scalar )  qs1d = 0.0_wp
+    IF ( humidity )  qs1d = 0.0_wp
 
 !
@@ -232 +234 @@
                humidity, intermediate_timestep_count,                          &
                intermediate_timestep_count_max, f, g, ibc_e_b, kappa,          &  
-               mixing_length_1d, passive_scalar,                               &
+               mixing_length_1d,                                               &
                simulated_time_chr, timestep_scheme, tsc, zeta_max, zeta_min
                 
@@ -619 +621 @@
                 e1d(nzb) = e1d(nzb+1)
 
-                IF ( humidity .OR. passive_scalar ) THEN
+                IF ( humidity ) THEN
 !
 !--                Compute q* 
                    IF ( rif1d(1) >= 0.0_wp )  THEN
 !
-!--                Stable stratification
-                   qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /            &
+!--                   Stable stratification
+                      qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /         &
                           ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * rif1d(nzb+1) * &
                                           ( zu(nzb+1) - z0h1d ) / zu(nzb+1)    &
                           )
-                ELSE
-!
-!--                Unstable stratification
-                   a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) )
-                   b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /                 &
-                                      zu(nzb+1) * z0h1d )
-!
-!--                In the borderline case the formula for stable stratification 
-!--                must be applied, because otherwise a zero division would
-!--                occur in the argument of the logarithm.
-                   IF ( a == 1.0_wp  .OR.  b == 1.0_wp )  THEN
-                      qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /         &
-                             ( LOG( zu(nzb+1) / z0h1d ) +                      &
-                               5.0_wp * rif1d(nzb+1) *                         &
-                               ( zu(nzb+1) - z0h1d ) / zu(nzb+1)               &
-                             )
                    ELSE
-                      qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /         &
-                             LOG( (a-1.0_wp) / (a+1.0_wp) *                    &
-                                  (b+1.0_wp) / (b-1.0_wp) )
-                   ENDIF
-                ENDIF                
+!
+!--                   Unstable stratification
+                      a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) )
+                      b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /              &
+                                         zu(nzb+1) * z0h1d )
+!
+!--                   In the borderline case the formula for stable stratification 
+!--                   must be applied, because otherwise a zero division would
+!--                   occur in the argument of the logarithm.
+                      IF ( a == 1.0_wp  .OR.  b == 1.0_wp )  THEN
+                         qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /      &
+                                ( LOG( zu(nzb+1) / z0h1d ) +                   &
+                                  5.0_wp * rif1d(nzb+1) *                      &
+                                  ( zu(nzb+1) - z0h1d ) / zu(nzb+1)            &
+                                )
+                      ELSE
+                         qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /      &
+                                LOG( (a-1.0_wp) / (a+1.0_wp) *                 &
+                                     (b+1.0_wp) / (b-1.0_wp) )
+                      ENDIF
+                   ENDIF                
                 ELSE
                    qs1d = 0.0_wp

palm/trunk/SOURCE/init_3d_model.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Separate humidity and passive scalar
+! Increase dimension for mean_inflow_profiles
+! Remove inadvertent write-statement
+! Bugfix, large-scale forcing is still not implemented for passive scalars
 !
 ! Former revisions:
@@ -469 +472 @@
     ENDIF
 
-    IF ( humidity  .OR.  passive_scalar )  THEN
-!
-!--    2D-humidity/scalar arrays
+    IF ( humidity )  THEN
+!
+!--    2D-humidity
        ALLOCATE ( qs(nysg:nyng,nxlg:nxrg),                                     &
                   qsws(nysg:nyng,nxlg:nxrg),                                   &
@@ -477 +480 @@
 
 !
-!--    3D-humidity/scalar arrays
+!--    3D-humidity
 #if defined( __nopointer )
        ALLOCATE( q(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
@@ -489 +492 @@
 
 !
-!--    3D-arrays needed for humidity only
+!--    3D-arrays needed for humidity
        IF ( humidity )  THEN
 #if defined( __nopointer )
@@ -570 +573 @@
        ENDIF
 
+    ENDIF
+    
+    
+    IF ( passive_scalar )  THEN
+!
+!--    2D-scalar arrays
+       ALLOCATE ( ss(nysg:nyng,nxlg:nxrg),                                     &
+                  ssws(nysg:nyng,nxlg:nxrg),                                   &
+                  sswst(nysg:nyng,nxlg:nxrg) )
+
+!
+!--    3D scalar arrays
+#if defined( __nopointer )
+       ALLOCATE( s(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
+                 s_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
+                 ts_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+#else
+       ALLOCATE( s_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
+                 s_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
+                 s_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+#endif
     ENDIF
 
@@ -678 +702 @@
     w  => w_1;   w_p  => w_2;   tw_m  => w_3
 
-    IF ( humidity  .OR.  passive_scalar )  THEN
+    IF ( humidity )  THEN
        q => q_1;  q_p => q_2;  tq_m => q_3
        IF ( humidity )  THEN
@@ -696 +720 @@
        ENDIF
     ENDIF
+    
+    IF ( passive_scalar )  THEN
+       s => s_1;  s_p => s_2;  ts_m => s_3
+    ENDIF    
 
     IF ( ocean )  THEN
@@ -773 +801 @@
           ENDDO
 
-          IF ( humidity  .OR.  passive_scalar )  THEN
+          IF ( humidity )  THEN
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
@@ -788 +816 @@
 
              ENDIF
+          ENDIF
+          IF ( passive_scalar )  THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   s(:,j,i) = s_init
+                ENDDO
+             ENDDO   
           ENDIF
 
@@ -830 +865 @@
 !--       This could actually be computed more accurately in the 1D model.
 !--       Update when opportunity arises!
-          IF ( humidity  .OR.  passive_scalar )  THEN
+          IF ( humidity )  THEN
              qs = 0.0_wp
              IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -837 +872 @@
              ENDIF
           ENDIF
+!
+!--       Initialize scaling parameter for passive scalar
+          IF ( passive_scalar ) ss = 0.0_wp          
 
 !
@@ -881 +919 @@
              u_init  = unudge(:,1)
              v_init  = vnudge(:,1)
-             IF ( humidity  .OR.  passive_scalar )  THEN
+             IF ( humidity  )  THEN ! is passive_scalar correct???
                 q_init = qnudge(:,1)
              ENDIF
@@ -915 +953 @@
           ENDIF
 
-          IF ( humidity  .OR.  passive_scalar )  THEN
+          IF ( humidity )  THEN
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
@@ -931 +969 @@
 
              ENDIF
+          ENDIF
+          
+          IF ( passive_scalar )  THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   s(:,j,i) = s_init
+                ENDDO
+             ENDDO
           ENDIF
 
@@ -975 +1021 @@
           vsws  = 0.0_wp
           vswst = top_momentumflux_v
-          IF ( humidity  .OR.  passive_scalar ) qs = 0.0_wp
+          IF ( humidity )       qs = 0.0_wp
+          IF ( passive_scalar ) ss = 0.0_wp
 
 !
@@ -1015 +1062 @@
 !
 !--    Calculate virtual potential temperature
-       IF ( humidity ) vpt = pt * ( 1.0_wp + 0.61_wp * q )
+       IF ( humidity )  vpt = pt * ( 1.0_wp + 0.61_wp * q )
 
 !
@@ -1053 +1100 @@
 !
 !--       Store initial scalar profile
-          hom(:,1,26,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
+          hom(:,1,115,:) = SPREAD(  s(:,nys,nxl), 2, statistic_regions+1 )
        ENDIF
 
@@ -1116 +1163 @@
 !
 !--       Determine the near-surface water flux
-          IF ( humidity  .OR.  passive_scalar )  THEN
+          IF ( humidity )  THEN
              IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
                 qrsws = 0.0_wp
@@ -1138 +1185 @@
              ENDIF
           ENDIF
+!
+!--       Initialize the near-surface scalar flux
+          IF ( passive_scalar )  THEN
+             IF ( constant_scalarflux )  THEN
+                ssws   = surface_scalarflux
+!
+!--             Over topography surface_scalarflux is replaced by 
+!--             wall_scalarflux(0)
+                IF ( TRIM( topography ) /= 'flat' )  THEN
+                   wall_sflux = wall_scalarflux
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         IF ( nzb_s_inner(j,i) /= 0 )  ssws(j,i) = wall_sflux(0)
+                      ENDDO
+                   ENDDO
+                ENDIF
+             ENDIF
+          ENDIF          
 
        ENDIF
@@ -1152 +1217 @@
              tswst = top_heatflux
 
-             IF ( humidity  .OR.  passive_scalar )  THEN
+             IF ( humidity )  THEN
                 qswst = 0.0_wp
                 IF ( cloud_physics  .AND.  microphysics_seifert ) THEN
@@ -1159 +1224 @@
                 ENDIF
              ENDIF
+             IF ( passive_scalar )  sswst = 0.0_wp
 
              IF ( ocean )  THEN
@@ -1192 +1258 @@
           ENDIF
 
-          IF ( humidity  .OR.  passive_scalar )  THEN
+          IF ( humidity  )  THEN
              IF (  .NOT.  constant_waterflux )  qsws   = 0.0_wp
              IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -1199 +1265 @@
              ENDIF
           ENDIF
+          IF ( passive_scalar  .AND.  .NOT.  constant_scalarflux )  ssws = 0.0_wp
 
        ENDIF
@@ -1261 +1328 @@
 !--    If required, change the surface humidity/scalar at the start of the 3D
 !--    run
-       IF ( ( humidity .OR. passive_scalar ) .AND.                             &
-            q_surface_initial_change /= 0.0_wp )  THEN
+       IF ( humidity  .AND.  q_surface_initial_change /= 0.0_wp )              &
           q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
-       ENDIF
+          
+       IF ( passive_scalar .AND.  s_surface_initial_change /= 0.0_wp )         &
+          s(nzb,:,:) = s(nzb,:,:) + s_surface_initial_change
+        
 
 !
@@ -1271 +1340 @@
        e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
 
-       IF ( humidity  .OR.  passive_scalar )  THEN
+       IF ( humidity  )  THEN
           tq_m = 0.0_wp
           q_p = q
@@ -1281 +1350 @@
           ENDIF
        ENDIF
+       
+       IF ( passive_scalar )  THEN
+          ts_m = 0.0_wp
+          s_p  = s
+       ENDIF       
 
        IF ( ocean )  THEN
@@ -1323 +1397 @@
 #endif
        ENDDO
-       write(9,*) "EOF read binary"
-       flush(9)
 
 !
@@ -1335 +1407 @@
 !--       profiles from the prerun. Alternatively, prescribed profiles
 !--       for u,v-components can be used.
-          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,6) )
+          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,7) )
 
           IF ( use_prescribed_profile_data )  THEN
@@ -1346 +1418 @@
           mean_inflow_profiles(:,4) = hom_sum(:,4,0)       ! pt
           mean_inflow_profiles(:,5) = hom_sum(:,8,0)       ! e
-          mean_inflow_profiles(:,6) = hom_sum(:,41,0)      ! q
+          IF ( humidity )                                                      &
+             mean_inflow_profiles(:,6) = hom_sum(:,41,0)   ! q
+          IF ( passive_scalar )                                                &
+             mean_inflow_profiles(:,7) = hom_sum(:,115,0)   ! s
 
 !
@@ -1373 +1448 @@
                    pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
                    e(k,j,nxlg:-1)  = mean_inflow_profiles(k,5)
-                   IF ( humidity  .OR.  passive_scalar )                       &
+                   IF ( humidity )                                             &
                       q(k,j,nxlg:-1)  = mean_inflow_profiles(k,6)
+                   IF ( passive_scalar )                                       &
+                      s(k,j,nxlg:-1)  = mean_inflow_profiles(k,7)                      
                 ENDDO
              ENDDO
@@ -1459 +1536 @@
 !--    including ghost points)
        e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
-       IF ( humidity  .OR.  passive_scalar )  THEN
+       IF ( humidity )  THEN
           q_p = q
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -1466 +1543 @@
           ENDIF
        ENDIF
-       IF ( ocean )  sa_p = sa
+       IF ( passive_scalar )  s_p  = s
+       IF ( ocean          )  sa_p = sa
 
 !
@@ -1473 +1551 @@
 !--    there before they are set. 
        te_m = 0.0_wp; tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
-       IF ( humidity  .OR.  passive_scalar )  THEN
+       IF ( humidity )  THEN
           tq_m = 0.0_wp
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -1480 +1558 @@
           ENDIF
        ENDIF
-       IF ( ocean )  tsa_m = 0.0_wp
+       IF ( passive_scalar )  ts_m  = 0.0_wp
+       IF ( ocean          )  tsa_m = 0.0_wp
 
        CALL location_message( 'finished', .TRUE. )

palm/trunk/SOURCE/modules.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Separate humidity and passive scalar
+! +bc_s_t_val, diss_s_s, diss_l_s, flux_s_s, flux_l_s, s, sp, s1, s2, s3, ssws_av, 
+!  s_init, s_surf, sums_wsss_ws_l, ss, ssws, sswst, ts_m, wall_sflux
+! +constant_scalarflux, ibc_s_b, ibc_s_t, s_vertical_gradient_level_ind
+!
+! Unused variables removed
+! -gamma_x, gamma_y, gamma_z, var_x, var_y, var_z
+!
+! Change default values (in order to unify gradient calculation): 
+! pt_vertical_gradient_level, sa_vertical_gradient_level
 !
 ! Former revisions:
@@ -411 +420 @@
           dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor, l_grid,           &
           ptdf_x, ptdf_y, p_surf, pt_surf, pt_init, qsws_surf, q_init, q_surf, &
-          rdf, rdf_sc, ref_state, sa_init, shf_surf, timenudge, time_surf,     &
-          time_vert, tmp_tnudge, ug, u_init, u_nzb_p1_for_vfc, vg, v_init,     &
-          v_nzb_p1_for_vfc, w_subs, zu, zw
+          rdf, rdf_sc, ref_state, s_init, s_surf, sa_init, shf_surf,           &
+          timenudge, time_surf, time_vert, tmp_tnudge, ug, u_init,             &
+          u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, w_subs, zu, zw
 
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::                                   &
           c_u, c_v, c_w, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q,             &
-          diss_s_qr, diss_s_sa, diss_s_u, diss_s_v, diss_s_w, dzu_mg, dzw_mg,  &
-          flux_s_e, flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr, flux_s_sa,      &
-          flux_s_u, flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg,                   &
-          mean_inflow_profiles, nrs, nrsws, nrswst,                            &
+          diss_s_qr, diss_s_s, diss_s_sa, diss_s_u, diss_s_v, diss_s_w, dzu_mg,&
+          dzw_mg, flux_s_e, flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr,         &
+          flux_s_s, flux_s_sa, flux_s_u, flux_s_v, flux_s_w, f1_mg, f2_mg,     &
+          f3_mg, mean_inflow_profiles, nrs, nrsws, nrswst,                     &
           ol,                                                                  & !< Obukhov length
           precipitation_amount, precipitation_rate, ptnudge, pt_slope_ref,     &
           qnudge, qs, qsws, qswst, qswst_remote, qrs, qrsws, qrswst,           &
-          saswsb, saswst, shf, tnudge, td_lsa_lpt, td_lsa_q, td_sub_lpt,       &
+          saswsb, saswst, shf, ss, ssws, sswst, tnudge, td_lsa_lpt, td_lsa_q,  &
+          td_sub_lpt,                                                          &
           td_sub_q, total_2d_a, total_2d_o, ts, tswst, ug_vert, unudge, us,    &
           usws, uswst, vnudge, vg_vert, vsws, vswst, wnudge, wsubs_vert,       &
@@ -433 +443 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
           d, de_dx, de_dy, de_dz, diss, diss_l_e,                              &
-          diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr, diss_l_sa, diss_l_u,      &
-          diss_l_v, diss_l_w, flux_l_e, flux_l_nr, flux_l_pt, flux_l_q,        &
-          flux_l_qr, flux_l_sa, flux_l_u, flux_l_v, flux_l_w, kh, km,          &
-          l_wall, prr, p_loc, tend, tric,                                      &
+          diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr, diss_l_s, diss_l_sa,      &
+          diss_l_u, diss_l_v, diss_l_w, flux_l_e, flux_l_nr, flux_l_pt,        &
+          flux_l_q, flux_l_qr, flux_l_s, flux_l_sa, flux_l_u, flux_l_v,        &
+          flux_l_w, kh, km, l_wall, prr, p_loc, tend, tric,                    &
           u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l,       &
           w_m_n, w_m_r, w_m_s
@@ -443 +453 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
           e, e_p, nr, nr_p, p, prho, pt, pt_p, q, q_p, qc, ql, ql_c, ql_v,     &
-          ql_vp, qr, qr_p, rho, sa, sa_p, te_m, tnr_m, tpt_m, tq_m, tqr_m,     &
-          tsa_m, tu_m, tv_m, tw_m, u, u_p, v, v_p, vpt, w, w_p
+          ql_vp, qr, qr_p, rho, s, s_p, sa, sa_p, te_m, tnr_m, tpt_m, tq_m,    &
+          tqr_m, ts_m, tsa_m, tu_m, tv_m, tw_m, u, u_p, v, v_p, vpt, w, w_p
 #else
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
           e_1, e_2, e_3, p, prho_1, nr_1, nr_2, nr_3, pt_1, pt_2, pt_3, q_1,   &
           q_2, q_3, qc_1, ql_v, ql_vp, ql_1, ql_2, qr_1, qr_2, qr_3, rho_1,    &
-          sa_1, sa_2, sa_3, u_1, u_2, u_3, v_1, v_2, v_3, vpt_1, w_1, w_2, w_3
+          s_1, s_2, s_3, sa_1, sa_2, sa_3, u_1, u_2, u_3, v_1, v_2, v_3, vpt_1,&
+          w_1, w_2, w_3
 
     REAL(wp), DIMENSION(:,:,:), POINTER ::                                     &
           e, e_p, nr, nr_p, prho, pt, pt_p, q, q_p, qc, ql, ql_c, qr, qr_p,    &
-          rho, sa, sa_p, te_m, tnr_m, tpt_m, tq_m, tqr_m, tsa_m, tu_m, tv_m,   &
-          tw_m, u, u_p, v, v_p, vpt, w, w_p
+          rho, s, s_p, sa, sa_p, te_m, tnr_m, tpt_m, tq_m, tqr_m, ts_m,        &
+          tsa_m, tu_m, tv_m, tw_m, u, u_p, v, v_p, vpt, w, w_p
 #endif
 
     REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  rif_wall, tri
-
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: var_x, var_y, var_z, gamma_x,   &
-                                           gamma_y, gamma_z
 
 
@@ -481 +489 @@
                                               ol_av,                 & !< Avg. of Obukhov length
                                               qsws_av,               & !< Avg. of surface moisture flux
+                                              ssws_av,               & !< Avg. of surface scalar flux
                                               shf_av,                & !< Avg. of surface heat flux
                                               ts_av,                 & !< Avg. of characteristic temperature scale
@@ -631 +640 @@
                      dz_stretch_level_index, ensemble_member_nr = 0, gamma_mg, gathered_size, &
                      grid_level, ibc_e_b, ibc_p_b, ibc_p_t, &
-                     ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, &
+                     ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, ibc_s_t, &
                      ibc_sa_t, ibc_uv_b, ibc_uv_t, &
                      inflow_disturbance_begin = -1, inflow_disturbance_end = -1, &
@@ -657 +666 @@
                      pt_vertical_gradient_level_ind(10) = -9999, &
                      q_vertical_gradient_level_ind(10) = -9999, &
+                     s_vertical_gradient_level_ind(10) = -9999, &                     
                      sa_vertical_gradient_level_ind(10) = -9999, &
                      section(100,3), section_xy(100) = -9999, &
@@ -685 +695 @@
                 constant_top_momentumflux = .FALSE., &
                 constant_top_salinityflux = .TRUE., &
+                constant_scalarflux = .TRUE., &              
                 constant_waterflux = .TRUE., create_disturbances = .TRUE., &
                 data_output_2d_on_each_pe = .TRUE., & 
@@ -739 +750 @@
                  alpha_surface = 0.0_wp, atmos_ocean_sign = 1.0_wp, &
                  averaging_interval = 0.0_wp, averaging_interval_pr = 9999999.9_wp, &
-                 bc_pt_t_val, bc_q_t_val, bottom_salinityflux = 0.0_wp, &
+                 bc_pt_t_val, bc_q_t_val, bc_s_t_val, bottom_salinityflux = 0.0_wp, &
                  building_height = 50.0_wp, building_length_x = 50.0_wp, &
                  building_length_y = 50.0_wp, building_wall_left = 9999999.9_wp, &
@@ -813 +824 @@
                  mask_scale(3), &
                  pt_vertical_gradient(10) = 0.0_wp, &
-                 pt_vertical_gradient_level(10) = -9999999.9_wp, &
+                 pt_vertical_gradient_level(10) = -1.0_wp, &
                  q_vertical_gradient(10) = 0.0_wp, &
                  q_vertical_gradient_level(10) = -1.0_wp, &
@@ -819 +830 @@
                  s_vertical_gradient_level(10) = -1.0_wp, &
                  sa_vertical_gradient(10) = 0.0_wp, &
-                 sa_vertical_gradient_level(10) = -9999999.9_wp, &
+                 sa_vertical_gradient_level(10) = -1.0_wp, &
                  skip_time_domask(max_masks) = 9999999.9_wp, threshold(20) = 0.0_wp, &
                  time_domask(max_masks) = 0.0_wp, &
@@ -833 +844 @@
                  wall_humidityflux(0:4) = 0.0_wp, wall_nrflux(0:4) = 0.0_wp, &
                  wall_qflux(0:4) = 0.0_wp, wall_qrflux(0:4) = 0.0_wp, &
-                 wall_salinityflux(0:4) = 0.0_wp, wall_scalarflux(0:4) = 0.0_wp, &
+                 wall_salinityflux(0:4) = 0.0_wp, wall_sflux(0:4) = 0.0_wp, wall_scalarflux(0:4) = 0.0_wp, &
                  subs_vertical_gradient(10) = 0.0_wp, &
                  subs_vertical_gradient_level(10) = -9999999.9_wp
@@ -1255 +1266 @@
                                                   sums_wsnrs_ws_l,                &
                                                   sums_wspts_ws_l,                &
+                                                  sums_wsqs_ws_l,                 &   
+                                                  sums_wsqrs_ws_l,                &
                                                   sums_wssas_ws_l,                &
-                                                  sums_wsqs_ws_l,                 &
-                                                  sums_wsqrs_ws_l,                &
+                                                  sums_wsss_ws_l,                 &
                                                   sums_ls_l
                                               

palm/trunk/SOURCE/netcdf_interface_mod.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Additional labels and units for timeseries output of passive scalar-related 
+! quantities
 ! 
 ! Former revisions:
@@ -174 +175 @@
              'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'number2' /)
 
-    INTEGER(iwp) ::  dots_num  = 31  !< number of timeseries defined by default
-    INTEGER(iwp) ::  dots_soil = 24  !< starting index for soil-timeseries
-    INTEGER(iwp) ::  dots_rad  = 32  !< starting index for radiation-timeseries
+    INTEGER(iwp) ::  dots_num  = 33  !< number of timeseries defined by default
+    INTEGER(iwp) ::  dots_soil = 26  !< starting index for soil-timeseries
+    INTEGER(iwp) ::  dots_rad  = 34  !< starting index for radiation-timeseries
 
     CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label =                    &
@@ -186 +187 @@
              'wpt          ', 'pt(0)        ', 'pt(z_mo)     ',                &
              'w"u"0        ', 'w"v"0        ', 'w"q"0        ',                &
-             'ol           ', 'q*           ', 'ghf_eb       ',                &
-             'shf_eb       ', 'qsws_eb      ', 'qsws_liq_eb  ',                &
-             'qsws_soil_eb ', 'qsws_veg_eb  ', 'r_a          ',                &
-             'r_s          ', 'rad_net      ', 'rad_lw_in    ',                &
-             'rad_lw_out   ', 'rad_sw_in    ', 'rad_sw_out   ',                &
-             'rrtm_aldif   ', 'rrtm_aldir   ', 'rrtm_asdif   ',                &
-             'rrtm_asdir   ',                                                  &
-             ( 'unknown      ', i9 = 1, dots_max-40 ) /)
+             'ol           ', 'q*           ', 'w"s"         ',                &
+             's*           ', 'ghf_eb       ', 'shf_eb       ',                &
+             'qsws_eb      ', 'qsws_liq_eb  ', 'qsws_soil_eb ',                &
+             'qsws_veg_eb  ', 'r_a          ', 'r_s          ',                &
+             'rad_net      ', 'rad_lw_in    ', 'rad_lw_out   ',                &
+             'rad_sw_in    ', 'rad_sw_out   ', 'rrtm_aldif   ',                &
+             'rrtm_aldir   ', 'rrtm_asdif   ', 'rrtm_asdir   ',                &                                               
+             ( 'unknown      ', i9 = 1, dots_max-42 ) /)
 
     CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_unit =                     &
@@ -203 +204 @@
              'K m/s        ', 'K            ', 'K            ',                &
              'm2/s2        ', 'm2/s2        ', 'kg m/s       ',                &
-             'm            ', 'kg/kg        ', '             ',                &
+             'm            ', 'kg/kg        ', 'kg m/(kg s)  ',                &
+             'kg/kg        ', '             ', '             ',                &
              '             ', '             ', '             ',                &
-             '             ', 'W/m2         ', 's/m          ',                &
-             '             ', 'W/m2         ', 'W/m2         ',                &
+             'W/m2         ', 's/m          ', '             ',                &
              'W/m2         ', 'W/m2         ', 'W/m2         ',                &
+             'W/m2         ', 'W/m2         ', '             ',                &
              '             ', '             ', '             ',                &
-             '             ',                                                  &
-             ( 'unknown      ', i9 = 1, dots_max-40 ) /)
+             ( 'unknown      ', i9 = 1, dots_max-42 ) /)
 
     CHARACTER (LEN=9), DIMENSION(300) ::  dopr_unit = 'unknown'

palm/trunk/SOURCE/palm.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Separate humidity and passive scalar
 ! 
 ! Former revisions:
@@ -361 +361 @@
           ENDIF
           IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp )
-          IF (humidity  .OR.  passive_scalar)  THEN
-             CALL exchange_horiz( q, nbgp )
-          ENDIF
+          IF ( humidity       )  CALL exchange_horiz( q, nbgp )
+          IF ( passive_scalar )  CALL exchange_horiz( s, nbgp )
        ENDIF
 

palm/trunk/SOURCE/parin.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Allocation of s_init
 ! 
 ! Former revisions:
@@ -218 +218 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  pt_init, q_init, ref_state, sa_init, ug, u_init, v_init,        &
+        ONLY:  pt_init, q_init, ref_state, s_init, sa_init, ug, u_init, v_init,&
                vg
 
@@ -597 +597 @@
 !--          ATTENTION: in case of changes to the following statement please
 !--                  also check the allocate statement in routine read_var_list
-             ALLOCATE( pt_init(0:nz+1), q_init(0:nz+1),                       &
+             ALLOCATE( pt_init(0:nz+1), q_init(0:nz+1), s_init(0:nz+1),        &
                        ref_state(0:nz+1), sa_init(0:nz+1), ug(0:nz+1),         &
                        u_init(0:nz+1), v_init(0:nz+1), vg(0:nz+1),             &

palm/trunk/SOURCE/plant_canopy_model_mod.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Separate humidity and passive scalar
 ! 
 ! Former revisions:
@@ -86 +86 @@
  
     USE arrays_3d,                                                             &
-        ONLY:  dzu, dzw, e, q, shf, tend, u, v, w, zu, zw 
+        ONLY:  dzu, dzw, e, q, s, shf, tend, u, v, w, zu, zw 
 
     USE indices,                                                               &
@@ -818 +818 @@
 
 !
-!--       scalar concentration
+!--       humidity
           CASE ( 5 )
              DO  i = nxl, nxr
@@ -866 +866 @@
                 ENDDO
              ENDDO 
+!
+!--       scalar concentration
+          CASE ( 7 )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
+                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
+                      tend(k,j,i) = tend(k,j,i) -                              &
+                                       lsec *                                  &
+                                       lad_s(kk,j,i) *                         &
+                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
+                                                          u(k,j,i+1) )         &
+                                             )**2 +                            &
+                                             ( 0.5_wp * ( v(k,j,i) +           &
+                                                          v(k,j+1,i) )         &
+                                             )**2 +                            &
+                                             ( 0.5_wp * ( w(k-1,j,i) +         & 
+                                                          w(k,j,i) )           &
+                                             )**2                              &
+                                           ) *                                 &
+                                       ( s(k,j,i) - lsc )
+                   ENDDO
+                ENDDO
+             ENDDO    
+
 
 
@@ -1146 +1171 @@
 
 !
-!--       scalar concentration 
+!--       humidity
           CASE ( 5 )
              DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
@@ -1186 +1211 @@
                                  e(k,j,i)
              ENDDO
+             
+!
+!--       scalar concentration 
+          CASE ( 7 )
+             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
+                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
+                tend(k,j,i) = tend(k,j,i) -                                    &
+                                 lsec *                                        &
+                                 lad_s(kk,j,i) *                               &
+                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
+                                                    u(k,j,i+1) )               &
+                                       )**2  +                                 &
+                                       ( 0.5_wp * ( v(k,j,i) +                 &
+                                                    v(k,j+1,i) )               &
+                                       )**2 +                                  &
+                                       ( 0.5_wp * ( w(k-1,j,i) +               &
+                                                    w(k,j,i) )                 &
+                                       )**2                                    &
+                                     ) *                                       &
+                                 ( s(k,j,i) - lsc )
+             ENDDO                
 
        CASE DEFAULT

palm/trunk/SOURCE/prognostic_equations.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Separate humidity and passive scalar
 ! 
 ! Former revisions:
@@ -192 +192 @@
     USE arrays_3d,                                                             &
         ONLY:  diss_l_e, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr,            &
-               diss_l_sa, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q,            &
-               diss_s_qr, diss_s_sa, e, e_p, flux_s_e, flux_s_nr, flux_s_pt,   &
-               flux_s_q, flux_s_qr, flux_s_sa, flux_l_e, flux_l_nr,            &
-               flux_l_pt, flux_l_q, flux_l_qr, flux_l_sa, nr, nr_p, nrsws,     &
-               nrswst, pt, ptdf_x, ptdf_y, pt_init, pt_p, prho, q, q_init,     &
-               q_p, qsws, qswst, qr, qr_p, qrsws, qrswst, rdf, rdf_sc,         &
-               ref_state, rho, sa, sa_init, sa_p, saswsb, saswst, shf, tend,   &
-               te_m, tnr_m, tpt_m, tq_m, tqr_m, tsa_m, tswst, tu_m, tv_m,      &
+               diss_l_s, diss_l_sa, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q,  &
+               diss_s_qr, diss_s_s, diss_s_sa, e, e_p, flux_s_e, flux_s_nr,    &
+               flux_s_pt, flux_s_q, flux_s_qr, flux_s_s, flux_s_sa, flux_l_e,  &
+               flux_l_nr, flux_l_pt, flux_l_q, flux_l_qr, flux_l_s, flux_l_sa, &
+               nr, nr_p, nrsws, nrswst, pt, ptdf_x, ptdf_y, pt_init, pt_p,     &
+               prho, q, q_init, q_p, qsws, qswst, qr, qr_p, qrsws, qrswst, rdf,&
+               rdf_sc, ref_state, rho, s, s_init, s_p, sa, sa_init, sa_p,      &
+               saswsb, saswst, shf, ssws, sswst, tend,                         &
+               te_m, tnr_m, tpt_m, tq_m, tqr_m, ts_m, tsa_m, tswst, tu_m, tv_m,&
                tw_m, u, ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
         
@@ -216 +217 @@
                use_upstream_for_tke, wall_heatflux,                            &
                wall_nrflux, wall_qflux, wall_qflux, wall_qflux, wall_qrflux,   &
-               wall_salinityflux, ws_scheme_mom, ws_scheme_sca
+               wall_salinityflux, wall_sflux, ws_scheme_mom, ws_scheme_sca
 
     USE cpulog,                                                                &
@@ -724 +725 @@
 
 !
-!--       If required, compute prognostic equation for total water content /
-!--       scalar
-          IF ( humidity  .OR.  passive_scalar )  THEN
+!--       If required, compute prognostic equation for total water content
+          IF ( humidity )  THEN
 
 !
@@ -745 +745 @@
 
 !
-!--          Sink or source of scalar concentration due to canopy elements
+!--          Sink or source of humidity due to canopy elements
              IF ( plant_canopy )  CALL pcm_tendency( i, j, 5 )
 
@@ -881 +881 @@
 
           ENDIF
-
+          
+!
+!--       If required, compute prognostic equation for scalar
+          IF ( passive_scalar )  THEN
+!
+!--          Tendency-terms for total water content / scalar
+             tend(:,j,i) = 0.0_wp
+             IF ( timestep_scheme(1:5) == 'runge' ) &
+             THEN
+                IF ( ws_scheme_sca )  THEN
+                   CALL advec_s_ws( i, j, s, 's', flux_s_s, & 
+                                diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
+                ELSE 
+                   CALL advec_s_pw( i, j, s )
+                ENDIF
+             ELSE
+                CALL advec_s_up( i, j, s )
+             ENDIF
+             CALL diffusion_s( i, j, s, ssws, sswst, wall_sflux )
+
+!
+!--          Sink or source of scalar concentration due to canopy elements
+             IF ( plant_canopy )  CALL pcm_tendency( i, j, 7 )
+
+!
+!--          Large scale advection, still need to be extended for scalars
+!              IF ( large_scale_forcing )  THEN
+!                 CALL ls_advec( i, j, simulated_time, 's' )
+!              ENDIF
+
+!
+!--          Nudging, still need to be extended for scalars
+!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' ) 
+
+!
+!--          If required compute influence of large-scale subsidence/ascent.
+!--          Note, the last argument is of no meaning in this case, as it is 
+!--          only used in conjunction with large_scale_forcing, which is to 
+!--          date not implemented for scalars.
+             IF ( large_scale_subsidence  .AND.                                &
+                  .NOT. use_subsidence_tendencies  .AND.                       &
+                  .NOT. large_scale_forcing )  THEN
+                CALL subsidence( i, j, tend, s, s_init, 3 )
+             ENDIF
+
+             CALL user_actions( i, j, 's-tendency' )
+
+!
+!--          Prognostic equation for scalar
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                s_p(k,j,i) = s(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +       &
+                                                  tsc(3) * ts_m(k,j,i) )       &
+                                      - tsc(5) * rdf_sc(k) *                   &
+                                        ( s(k,j,i) - s_init(k) )
+                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
+             ENDDO
+
+!
+!--          Calculate tendencies for the next Runge-Kutta step
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( intermediate_timestep_count == 1 )  THEN
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      ts_m(k,j,i) = tend(k,j,i)
+                   ENDDO
+                ELSEIF ( intermediate_timestep_count < &
+                         intermediate_timestep_count_max )  THEN
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
+                                     5.3125_wp * ts_m(k,j,i)
+                   ENDDO
+                ENDIF
+             ENDIF
+
+          ENDIF          
 !
 !--       If required, compute prognostic equation for turbulent kinetic 
@@ -1432 +1505 @@
 
 !
-!-- If required, compute prognostic equation for total water content / scalar
-    IF ( humidity  .OR.  passive_scalar )  THEN
-
-       CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
+!-- If required, compute prognostic equation for total water content
+    IF ( humidity )  THEN
+
+       CALL cpu_log( log_point(29), 'q-equation', 'start' )
 
 !
@@ -1470 +1543 @@
        
 !
-!--    Sink or source of scalar concentration due to canopy elements
+!--    Sink or source of humidity due to canopy elements
        IF ( plant_canopy ) CALL pcm_tendency( 5 )
 
@@ -1493 +1566 @@
 
 !
-!--    Prognostic equation for total water content / scalar
+!--    Prognostic equation for total water content
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -1529 +1602 @@
        ENDIF
 
-       CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
+       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
 
 !
@@ -1686 +1759 @@
 
     ENDIF
-
+!
+!-- If required, compute prognostic equation for scalar
+    IF ( passive_scalar )  THEN
+
+       CALL cpu_log( log_point(66), 's-equation', 'start' )
+
+!
+!--    Scalar/q-tendency terms with communication
+       sbt = tsc(2)
+       IF ( scalar_advec == 'bc-scheme' )  THEN
+
+          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--          Bott-Chlond scheme always uses Euler time step. Thus:
+             sbt = 1.0_wp
+          ENDIF
+          tend = 0.0_wp
+          CALL advec_s_bc( s, 's' )
+
+       ENDIF
+
+!
+!--    Scalar/q-tendency terms with no communication
+       IF ( scalar_advec /= 'bc-scheme' )  THEN
+          tend = 0.0_wp
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( ws_scheme_sca )  THEN
+                CALL advec_s_ws( s, 's' )
+             ELSE
+                CALL advec_s_pw( s )
+             ENDIF
+          ELSE
+             CALL advec_s_up( s )
+          ENDIF
+       ENDIF
+
+       CALL diffusion_s( s, ssws, sswst, wall_sflux )
+       
+!
+!--    Sink or source of humidity due to canopy elements
+       IF ( plant_canopy ) CALL pcm_tendency( 7 )
+
+!
+!--    Large scale advection. Not implemented for scalars so far.
+!        IF ( large_scale_forcing )  THEN
+!           CALL ls_advec( simulated_time, 'q' )
+!        ENDIF
+
+!
+!--    Nudging. Not implemented for scalars so far.
+!        IF ( nudging )  CALL nudge( simulated_time, 'q' ) 
+
+!
+!--    If required compute influence of large-scale subsidence/ascent.
+!--    Not implemented for scalars so far.
+       IF ( large_scale_subsidence  .AND.                                      &
+            .NOT. use_subsidence_tendencies  .AND.                             &
+            .NOT. large_scale_forcing )  THEN
+         CALL subsidence( tend, s, s_init, 3 )
+       ENDIF
+
+       CALL user_actions( 's-tendency' )
+
+!
+!--    Prognostic equation for total water content
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                s_p(k,j,i) = s(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +          &
+                                                  tsc(3) * ts_m(k,j,i) )       &
+                                      - tsc(5) * rdf_sc(k) *                   &
+                                        ( s(k,j,i) - s_init(k) )
+                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
+             ENDDO
+          ENDDO
+       ENDDO
+
+!
+!--    Calculate tendencies for the next Runge-Kutta step
+       IF ( timestep_scheme(1:5) == 'runge' )  THEN
+          IF ( intermediate_timestep_count == 1 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      ts_m(k,j,i) = tend(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSEIF ( intermediate_timestep_count < &
+                   intermediate_timestep_count_max )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * ts_m(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+       ENDIF
+
+       CALL cpu_log( log_point(66), 's-equation', 'stop' )
+
+    ENDIF
 !
 !-- If required, compute prognostic equation for turbulent kinetic 
@@ -2226 +2401 @@
 
 !
-!-- If required, compute prognostic equation for total water content / scalar
-    IF ( humidity  .OR.  passive_scalar )  THEN
-
-       CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
+!-- If required, compute prognostic equation for total water content
+    IF ( humidity )  THEN
+
+       CALL cpu_log( log_point(29), 'q-equation', 'start' )
 
 !
@@ -2307 +2482 @@
        ENDDO
 
-       CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
+       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
 
 !
@@ -2430 +2605 @@
     ENDIF
 
+!
+!-- If required, compute prognostic equation for scalar
+    IF ( passive_scalar )  THEN
+
+       CALL cpu_log( log_point(66), 's-equation', 'start' )
+
+!
+!--    Scalar/q-tendency terms with communication
+       sbt = tsc(2)
+       IF ( scalar_advec == 'bc-scheme' )  THEN
+
+          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--          Bott-Chlond scheme always uses Euler time step. Thus:
+             sbt = 1.0_wp
+          ENDIF
+          tend = 0.0_wp
+          CALL advec_s_bc( s, 's' )
+
+       ENDIF
+
+!
+!--    Scalar/q-tendency terms with no communication
+       IF ( scalar_advec /= 'bc-scheme' )  THEN
+          tend = 0.0_wp
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( ws_scheme_sca )  THEN
+                CALL advec_s_ws( s, 's' )
+             ELSE
+                CALL advec_s_pw( s )
+             ENDIF
+          ELSE
+             CALL advec_s_up( s )
+          ENDIF
+       ENDIF
+
+       CALL diffusion_s( s, ssws, sswst, wall_sflux )
+
+!
+!--    Sink or source of scalar concentration due to canopy elements
+       IF ( plant_canopy ) CALL pcm_tendency( 7 )
+
+!
+!--    Large scale advection. Not implemented so far.
+!        IF ( large_scale_forcing )  THEN
+!           CALL ls_advec( simulated_time, 's' )
+!        ENDIF
+
+!
+!--    Nudging. Not implemented so far.
+!        IF ( nudging )  CALL nudge( simulated_time, 's' ) 
+
+!
+!--    If required compute influence of large-scale subsidence/ascent.
+!--    Not implemented so far.
+       IF ( large_scale_subsidence  .AND.                                      &
+            .NOT. use_subsidence_tendencies  .AND.                             &
+            .NOT. large_scale_forcing )  THEN
+         CALL subsidence( tend, s, s_init, 3 )
+       ENDIF
+
+       CALL user_actions( 's-tendency' )
+
+!
+!--    Prognostic equation for total water content / scalar
+       DO  i = i_left, i_right
+          DO  j = j_south, j_north
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                s_p(k,j,i) = s(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +          &
+                                                  tsc(3) * ts_m(k,j,i) )       &
+                                      - tsc(5) * rdf_sc(k) *                   &
+                                        ( s(k,j,i) - s_init(k) )
+                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
+!
+!--             Tendencies for the next Runge-Kutta step
+                IF ( runge_step == 1 )  THEN
+                   ts_m(k,j,i) = tend(k,j,i)
+                ELSEIF ( runge_step == 2 )  THEN
+                   ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * ts_m(k,j,i)
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
+
+       CALL cpu_log( log_point(66), 's-equation', 'stop' )
+
+    ENDIF
 !
 !-- If required, compute prognostic equation for turbulent kinetic

palm/trunk/SOURCE/read_3d_binary.f90

@@ -107 +107 @@
         ONLY:  e, kh, km, ol, p, pt, q, ql, qc, nr, nrs, nrsws, nrswst,        &
                prr, precipitation_amount, qr,                                  &
-               qrs, qrsws, qrswst, qs, qsws, qswst, sa, saswsb, saswst,        &
-               rif_wall, shf, ts, tswst, u, u_m_l, u_m_n, u_m_r, u_m_s, us,    &
-               usws, uswst, v, v_m_l, v_m_n, v_m_r, v_m_s, vpt, vsws, vswst,   &
-               w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h, z0q
+               qrs, qrsws, qrswst, qs, qsws, qswst, s, sa, saswsb, saswst,     &
+               ss, ssws, sswst, rif_wall, shf, ss, ts, tswst, u, u_m_l, u_m_n, &
+               u_m_r, u_m_s, us, usws, uswst, v, v_m_l, v_m_n, v_m_r, v_m_s,   &
+               vpt, vsws, vswst, w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h, z0q
 
     USE averaging
@@ -985 +985 @@
                    rif_wall(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp,:) = &
                             tmp_4d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp,:)
+                            
+                CASE ( 's' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   s(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 's_av' )
@@ -1067 +1072 @@
                       ENDIF
                    ENDIF
-
+                   
+                CASE ( 'ss' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   ss(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'ssws' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   ssws(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'ssws_av' )
+                   IF ( .NOT. ALLOCATED( ssws_av ) )  THEN
+                      ALLOCATE( ssws_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   ssws_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'sswst' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   sswst(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 
+                 
                 CASE ( 'ts' )
                    IF ( k == 1 )  READ ( 13 )  tmp_2d

palm/trunk/SOURCE/read_var_list.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Separate humidity and passive scalar
+! Remove unused variables from ONLY list
 !
 ! Former revisions:
@@ -175 +176 @@
     USE arrays_3d,                                                             &
         ONLY:  inflow_damping_factor, mean_inflow_profiles, pt_init,           &
-               q_init, ref_state, sa_init, u_init, ug, v_init, vg
+               q_init, ref_state, s_init, sa_init, u_init, ug, v_init, vg
 
     USE control_parameters
@@ -292 +293 @@
        ALLOCATE( ug(0:nz+1), u_init(0:nz+1), vg(0:nz+1),                       &
                  v_init(0:nz+1), pt_init(0:nz+1), q_init(0:nz+1),              &
-                 ref_state(0:nz+1), sa_init(0:nz+1),                           &
+                 ref_state(0:nz+1), s_init(0:nz+1), sa_init(0:nz+1),           &
                  hom(0:nz+1,2,pr_palm+max_pr_user,0:statistic_regions),        &
                  hom_sum(0:nz+1,pr_palm+max_pr_user,0:statistic_regions) )
@@ -486 +487 @@
           CASE ( 'mean_inflow_profiles' )
              IF ( .NOT. ALLOCATED( mean_inflow_profiles ) )  THEN
-                ALLOCATE( mean_inflow_profiles(0:nz+1,6) )
+                ALLOCATE( mean_inflow_profiles(0:nz+1,7) )
              ENDIF
              READ ( 13 )  mean_inflow_profiles
@@ -597 +598 @@
           CASE ( 'run_coupled' )
              READ ( 13 )  run_coupled
+          CASE ( 's_init' )
+             READ ( 13 )  s_init
+          CASE ( 's_surface' )
+             READ ( 13 )  s_surface
+          CASE ( 's_surface_initial_change' )
+             READ ( 13 )  s_surface_initial_change
+          CASE ( 's_vertical_gradient' )
+             READ ( 13 )  s_vertical_gradient
+          CASE ( 's_vertical_gradient_level' )
+             READ ( 13 )  s_vertical_gradient_level
+          CASE ( 's_vertical_gradient_level_ind' )
+             READ ( 13 )  s_vertical_gradient_level_ind 
           CASE ( 'sa_init' )
              READ ( 13 )  sa_init
@@ -617 +630 @@
           CASE ( 'surface_waterflux' )
              READ ( 13 )  surface_waterflux      
-          CASE ( 's_surface' )
-             READ ( 13 )  s_surface
-          CASE ( 's_surface_initial_change' )
-             READ ( 13 )  s_surface_initial_change
-          CASE ( 's_vertical_gradient' )
-             READ ( 13 )  s_vertical_gradient
-          CASE ( 's_vertical_gradient_level' )
-             READ ( 13 )  s_vertical_gradient_level
           CASE ( 'time_coupling' )
              READ ( 13 )  time_coupling
@@ -789 +794 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  inflow_damping_factor, mean_inflow_profiles, pt_init,           &
-               q_init, ref_state, sa_init, u_init, ug, v_init, vg
+        ONLY:  inflow_damping_factor, mean_inflow_profiles, ref_state, ug, vg
 
     USE control_parameters

palm/trunk/SOURCE/spectra_mod.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Additional default spectra for passive scalar
 ! 
 ! Former revisions:
@@ -480 +480 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  d, pt, q, u, v, w
+           ONLY:  d, pt, q, s, u, v, w
 
        USE indices,                                                            &
@@ -532 +532 @@
           pr = 41
           d(nzb+1:nzt,nys:nyn,nxl:nxr) = q(nzb+1:nzt,nys:nyn,nxl:nxr)
+          
+       CASE ( 's' )
+          pr = 117
+          d(nzb+1:nzt,nys:nyn,nxl:nxr) = s(nzb+1:nzt,nys:nyn,nxl:nxr)
        
        CASE DEFAULT

palm/trunk/SOURCE/sum_up_3d_data.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Scalar surface flux added
 ! 
 ! Former revisions:
@@ -101 +101 @@
     USE arrays_3d,                                                             &
         ONLY:  dzw, e, nr, ol, p, pt, precipitation_rate, q, qc, ql, ql_c,     &
-               ql_v, qr, qsws, rho, sa, shf, ts, u, us, v, vpt, w, z0, z0h, z0q
+               ql_v, qr, qsws, rho, sa, shf, ssws, ts, u, us, v, vpt, w, z0,   &
+               z0h, z0q
 
     USE averaging,                                                             &
@@ -107 +108 @@
                precipitation_rate_av, pt_av, q_av, qc_av, ql_av, ql_c_av,      &
                ql_v_av, ql_vp_av, qr_av, qsws_av, qv_av, rho_av, s_av, sa_av,  &
-               shf_av, ts_av, u_av, us_av, v_av, vpt_av, w_av, z0_av, z0h_av,  &
-               z0q_av
+               shf_av, ssws_av, ts_av, u_av, us_av, v_av, vpt_av, w_av, z0_av, &
+               z0h_av, z0q_av
 
     USE cloud_parameters,                                                      &
@@ -949 +950 @@
              ENDDO
 
+          CASE ( 'ssws*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   ssws_av(j,i) = ssws_av(j,i) + ssws(j,i)
+                ENDDO
+             ENDDO
+
           CASE ( 't*' )
              DO  i = nxlg, nxrg

palm/trunk/SOURCE/surface_layer_fluxes_mod.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-!  
+! Treat humidity and passive scalar separately
 ! 
 ! Former revisions:
@@ -141 +141 @@
     USE arrays_3d,                                                             &
         ONLY:  e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, &
-               shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h, z0q
+               s, shf, ss, ssws, ts, u, us, usws, v, vpt, vsws, zu, zw, z0,    &
+               z0h, z0q
 
     USE cloud_parameters,                                                      &
@@ -152 +153 @@
 
     USE control_parameters,                                                    &
-        ONLY:  cloud_physics, constant_heatflux, constant_waterflux,           &
-               coupling_mode, g, humidity, ibc_e_b, ibc_pt_b,                  &
-               initializing_actions, kappa, intermediate_timestep_count,       &
+        ONLY:  cloud_physics, constant_heatflux, constant_scalarflux,          &
+               constant_waterflux, coupling_mode, g, humidity, ibc_e_b,        &
+               ibc_pt_b, initializing_actions, kappa,                          &
+               intermediate_timestep_count,                                    &
                intermediate_timestep_count_max, large_scale_forcing, lsf_surf, &
                message_string, microphysics_seifert, most_method, neutral,     &
@@ -859 +861 @@
 !
 !--    If required compute q*
-       IF ( humidity  .OR.  passive_scalar )  THEN
+       IF ( humidity )  THEN
           IF ( constant_waterflux )  THEN
 !
@@ -922 +924 @@
           ENDIF
        ENDIF
+       
+!
+!--    If required compute s*
+       IF ( passive_scalar )  THEN
+          IF ( constant_scalarflux )  THEN
+!
+!--          For a given water flux in the surface layer
+             !$OMP PARALLEL DO
+             !$acc kernels loop private( j )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   ss(j,i) = -ssws(j,i) / ( us(j,i) + 1E-30_wp )
+                ENDDO
+             ENDDO
+          ENDIF
+       ENDIF
 
 
@@ -1043 +1061 @@
 
 !
-!-- Compute the vertical water/scalar flux
-       IF (  .NOT.  constant_waterflux  .AND.  ( humidity  .OR.                &
-             passive_scalar )  .AND.  ( simulated_time <= skip_time_do_lsm     &
+!--    Compute the vertical water flux
+       IF (  .NOT.  constant_waterflux  .AND.  humidity  .AND.                 &
+             ( simulated_time <= skip_time_do_lsm                              &
             .OR.  .NOT.  land_surface ) )  THEN
           !$OMP PARALLEL DO
@@ -1057 +1075 @@
 
        ENDIF
+       
+!
+!--    Compute the vertical scalar flux
+       IF (  .NOT.  constant_scalarflux  .AND.  passive_scalar  .AND.          &
+             ( simulated_time <= skip_time_do_lsm                              &
+            .OR.  .NOT.  land_surface ) )  THEN
+          !$OMP PARALLEL DO
+          !$acc kernels loop independent present( qs, qsws, us )
+          DO  i = nxlg, nxrg
+             !$acc loop independent
+             DO  j = nysg, nyng
+                ssws(j,i) = -ss(j,i) * us(j,i)
+             ENDDO
+          ENDDO
+
+       ENDIF       
 
 !

palm/trunk/SOURCE/swap_timelevel.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-!
+! Separate humidity and passive scalar
 !
 ! Former revisions:
@@ -83 +83 @@
 #if defined( __nopointer )
     USE arrays_3d,                                                             &
-        ONLY:  e, e_p, nr, nr_p, pt, pt_p, q, q_p, qr, qr_p, sa, sa_p, u, u_p, &
-               v, v_p, w, w_p
+        ONLY:  e, e_p, nr, nr_p, pt, pt_p, q, q_p, qr, qr_p, s, s_p, sa, sa_p, &
+               u, u_p, v, v_p, w, w_p
 #else
     USE arrays_3d,                                                             &
         ONLY:  e, e_1, e_2, e_p, nr, nr_1, nr_2, nr_p, pt, pt_1, pt_2, pt_p, q,&
-               q_1, q_2, q_p, qr, qr_1, qr_2, qr_p, sa, sa_1, sa_2, sa_p, u,   &
-               u_1, u_2, u_p, v, v_1, v_2, v_p, w, w_1, w_2, w_p
+               q_1, q_2, q_p, qr, qr_1, qr_2, qr_p, s, s_1, s_2, s_p, sa, sa_1,&
+               sa_2, sa_p, u, u_1, u_2, u_p, v, v_1, v_2, v_p, w, w_1, w_2, w_p
 
 #endif
@@ -162 +162 @@
        sa = sa_p
     ENDIF
-    IF ( humidity  .OR.  passive_scalar )  THEN
+    IF ( humidity )  THEN
        q = q_p             
        IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -169 +169 @@
        ENDIF
     ENDIF
+    IF ( passive_scalar )  s = s_p             
 
     IF ( land_surface )  THEN
@@ -195 +196 @@
              sa => sa_1;  sa_p => sa_2
           ENDIF
-          IF ( humidity  .OR.  passive_scalar )  THEN
+          IF ( humidity )  THEN
              q => q_1;    q_p => q_2
              IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -202 +203 @@
              ENDIF
           ENDIF
+          IF ( passive_scalar )  THEN
+             s => s_1;    s_p => s_2
+          ENDIF
 
           IF ( land_surface )  THEN
@@ -223 +227 @@
              sa => sa_2;  sa_p => sa_1
           ENDIF
-          IF ( humidity  .OR.  passive_scalar )  THEN
+          IF ( humidity )  THEN
              q => q_2;    q_p => q_1
              IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -230 +234 @@
              ENDIF
           ENDIF
+          IF ( passive_scalar )  THEN
+             s => s_2;    s_p => s_1
+          ENDIF
 
           IF ( land_surface )  THEN

palm/trunk/SOURCE/time_integration.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Separate humidity and passive scalar
 ! 
 ! Former revisions:
@@ -225 +225 @@
     USE arrays_3d,                                                             &
         ONLY:  diss, dzu, e, e_p, nr_p, prho, pt, pt_p, pt_init, q_init, q,    &
-               ql, ql_c, ql_v, ql_vp, qr_p, q_p, ref_state, rho, sa_p, tend,   &
-               u, u_p, v, vpt, v_p, w, w_p
+               ql, ql_c, ql_v, ql_vp, qr_p, q_p, ref_state, rho, s, s_p, sa_p, &
+               tend, u, u_p, v, vpt, v_p, w, w_p
 
     USE calc_mean_profile_mod,                                                 &
@@ -532 +532 @@
                    CALL exchange_horiz( sa_p, nbgp )
                    CALL exchange_horiz( rho, nbgp )
-                  CALL exchange_horiz( prho, nbgp )
-                ENDIF
-                IF (humidity  .OR.  passive_scalar)  THEN
+                   CALL exchange_horiz( prho, nbgp )
+                ENDIF
+                IF ( humidity )  THEN
                    CALL exchange_horiz( q_p, nbgp )
                    IF ( cloud_physics .AND. microphysics_seifert )  THEN
@@ -551 +551 @@
                    CALL exchange_horiz( diss, nbgp )
                 ENDIF
+                IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )
 
                 IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
@@ -601 +602 @@
                   CALL exchange_horiz( prho, nbgp )
                 ENDIF
-                IF (humidity  .OR.  passive_scalar)  THEN
+                IF ( humidity )  THEN
                    CALL exchange_horiz( q_p, nbgp )
                    IF ( cloud_physics .AND. microphysics_seifert )  THEN
@@ -618 +619 @@
                    CALL exchange_horiz( diss, nbgp )
                 ENDIF
+                IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )
 
                 IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
@@ -690 +692 @@
                 CALL exchange_horiz( prho, nbgp )
              ENDIF
-             IF (humidity  .OR.  passive_scalar)  THEN
+             IF ( humidity )  THEN
                 CALL exchange_horiz( q_p, nbgp )
                 IF ( cloud_physics .AND. microphysics_seifert )  THEN
@@ -707 +709 @@
                 CALL exchange_horiz( diss, nbgp )
              ENDIF
+             IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )      
 
              IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
@@ -746 +749 @@
                 CALL exchange_horiz( v, nbgp )
                 CALL exchange_horiz( w, nbgp )
-                IF ( .NOT. neutral )  THEN
-                   CALL exchange_horiz( pt, nbgp )
-                ENDIF
-                IF ( humidity  .OR.  passive_scalar )  THEN
-                   CALL exchange_horiz( q, nbgp )
-                ENDIF
-                IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp )
+                IF ( .NOT. neutral            )  CALL exchange_horiz( pt, nbgp )
+                IF ( humidity                 )  CALL exchange_horiz( q, nbgp  )
+                IF ( passive_scalar           )  CALL exchange_horiz( s, nbgp  )
+                IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp  )
              ENDIF
 !

palm/trunk/SOURCE/user_actions.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! New CASE statement for scalar tendency
 ! 
 ! Former revisions:
@@ -170 +170 @@
 
           CASE ( 'q-tendency' )
+          
+          
+          CASE ( 's-tendency' )          
 
 
@@ -228 +231 @@
 
           CASE ( 'q-tendency' )
+          
+          
+          CASE ( 's-tendency' )
 
 

palm/trunk/SOURCE/user_spectra.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Consider additional spectra for passive scalar
 ! 
 ! Former revisions:
@@ -89 +89 @@
        SELECT CASE ( TRIM( data_output_sp(m) ) )
           
-          CASE ( 'u', 'v', 'w', 'pt', 'q' )
+          CASE ( 'u', 'v', 'w', 'pt', 'q', 's' )
 !--          Not allowed here since these are the standard quantities used in 
 !--          preprocess_spectra.
@@ -108 +108 @@
        SELECT CASE ( TRIM( data_output_sp(m) ) )
 
-          CASE ( 'u', 'v', 'w', 'pt', 'q' )
+          CASE ( 'u', 'v', 'w', 'pt', 'q', 's' )
 !--          Not allowed here since these are the standard quantities used in 
 !--          data_output_spectra.

palm/trunk/SOURCE/virtual_flight_mod.f90

@@ -19 +19 @@
 ! Current revisions:
 ! ------------------
-! 
+! Separate humidity and passive scalar.
+! Bugfix concerning labeling of timeseries.
 ! 
 ! Former revisions:
@@ -326 +327 @@
        
        INTEGER(iwp) ::  i               !< loop variable
+       INTEGER(iwp) ::  id_pt           !< identifyer for labeling
+       INTEGER(iwp) ::  id_q            !< identifyer for labeling
+       INTEGER(iwp) ::  id_ql           !< identifyer for labeling
+       INTEGER(iwp) ::  id_s            !< identifyer for labeling        
+       INTEGER(iwp) ::  id_u = 1        !< identifyer for labeling  
+       INTEGER(iwp) ::  id_v = 2        !< identifyer for labeling
+       INTEGER(iwp) ::  id_w = 3        !< identifyer for labeling
        INTEGER(iwp) ::  k               !< index variable
        
@@ -332 +340 @@
 !--    Define output quanities, at least three variables are measured (u,v,w)
        num_var_fl = 3
-       IF ( .NOT. neutral )                      num_var_fl = num_var_fl + 1
-       IF ( humidity .OR. passive_scalar )       num_var_fl = num_var_fl + 1
-       IF ( cloud_physics .OR. cloud_droplets )  num_var_fl = num_var_fl + 1  
+       IF ( .NOT. neutral                     )  THEN
+          num_var_fl = num_var_fl + 1
+          id_pt      = num_var_fl
+       ENDIF
+       IF ( humidity                          )  THEN
+          num_var_fl = num_var_fl + 1
+          id_q       = num_var_fl
+       ENDIF
+       IF ( cloud_physics .OR. cloud_droplets )  THEN
+          num_var_fl = num_var_fl + 1  
+          id_ql      = num_var_fl
+       ENDIF
+       IF ( passive_scalar                    )  THEN
+          num_var_fl = num_var_fl + 1
+          id_s       = num_var_fl
+       ENDIF
 !
 !--    Write output strings for dimensions x, y, z
@@ -364 +385 @@
           DO i=1, num_var_fl
 
-             SELECT CASE ( i )
-                      
-                CASE( 1 )
-                   dofl_label(k) = TRIM( label_leg ) // '_u'
-                   dofl_unit(k)  = 'm/s'
-                   k             = k + 1
-
-                CASE( 2 )
-                   dofl_label(k) = TRIM( label_leg ) // '_v'
-                   dofl_unit(k)  = 'm/s'
-                   k             = k + 1
-
-                CASE( 3 )
-                   dofl_label(k) = TRIM( label_leg ) // '_w'
-                   dofl_unit(k)  = 'm/s'
-                   k             = k + 1
-
-                CASE( 4 )
-                   dofl_label(k) = TRIM( label_leg ) // '_pt'
-                   dofl_unit(k)  = 'K'
-                   k             = k + 1
-
-                CASE( 5 )
-                   dofl_label(k) = TRIM( label_leg ) // '_q'
-                   dofl_unit(k)  = 'kg/kg'
-                   k             = k + 1
-
-                CASE( 6 )
-                   dofl_label(k) = TRIM( label_leg ) // '_ql'
-                   dofl_unit(k)  = 'kg/kg'
-                   k             = k + 1
-
-             END SELECT
-
+             IF ( i == id_u      )  THEN         
+                dofl_label(k) = TRIM( label_leg ) // '_u'
+                dofl_unit(k)  = 'm/s'
+                k             = k + 1
+             ELSEIF ( i == id_v  )  THEN       
+                dofl_label(k) = TRIM( label_leg ) // '_v'
+                dofl_unit(k)  = 'm/s'
+                k             = k + 1
+             ELSEIF ( i == id_w  )  THEN          
+                dofl_label(k) = TRIM( label_leg ) // '_w'
+                dofl_unit(k)  = 'm/s'
+                k             = k + 1
+             ELSEIF ( i == id_pt )  THEN       
+                dofl_label(k) = TRIM( label_leg ) // '_pt'
+                dofl_unit(k)  = 'K'
+                k             = k + 1
+             ELSEIF ( i == id_q  )  THEN       
+                dofl_label(k) = TRIM( label_leg ) // '_q'
+                dofl_unit(k)  = 'kg/kg'
+                k             = k + 1
+             ELSEIF ( i == id_ql )  THEN       
+                dofl_label(k) = TRIM( label_leg ) // '_ql'
+                dofl_unit(k)  = 'kg/kg'
+                k             = k + 1
+             ELSEIF ( i == id_s  )  THEN                          
+                dofl_label(k) = TRIM( label_leg ) // '_s'
+                dofl_unit(k)  = 'kg/kg'
+                k             = k + 1
+             ENDIF
           ENDDO
           
@@ -421 +437 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, pt, q, ql, u, v, w, zu, zw
+           ONLY:  ddzu, ddzw, pt, q, ql, s, u, v, w, zu, zw
 
        USE control_parameters,                                                 &
@@ -596 +612 @@
                 ENDIF
 !
-!--             Humidity or passive scalar
-                IF ( humidity .OR. passive_scalar )  THEN
+!--             Humidity
+                IF ( humidity )  THEN
                    CALL interpolate_xyz( q, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
                    var_index = var_index + 1
@@ -605 +621 @@
                 IF ( cloud_physics .OR. cloud_droplets )  THEN
                    CALL interpolate_xyz( ql, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
+                   var_index = var_index + 1
+                ENDIF
+!
+!--             Passive scalar
+                IF ( passive_scalar )  THEN
+                   CALL interpolate_xyz( s, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
                    var_index = var_index + 1
                 ENDIF

palm/trunk/SOURCE/write_3d_binary.f90

@@ -100 +100 @@
         ONLY:  e, kh, km, ol, p, pt, q, ql, qc, nr, nrs, nrsws, nrswst,        &
                prr, precipitation_amount, qr,                                  &
-               qrs, qrsws, qrswst, qs, qsws, qswst, sa, saswsb, saswst,        &
-               rif_wall, shf, ts, tswst, u, u_m_l, u_m_n, u_m_r, u_m_s, us,    &
-               usws, uswst, v, v_m_l, v_m_n, v_m_r, v_m_s, vpt, vsws, vswst,   &
-               w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h, z0q
+               qrs, qrsws, qrswst, qs, qsws, qswst, s, sa, ssws, sswst,        &
+               saswsb, saswst, rif_wall, shf, ss, ts, tswst, u, u_m_l, u_m_n,  &
+               u_m_r, u_m_s, us, usws, uswst, v, v_m_l, v_m_n, v_m_r, v_m_s,   &
+               vpt, vsws, vswst, w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h, z0q
         
     USE averaging
@@ -237 +237 @@
        WRITE ( 14 )  'pt_av               ';  WRITE ( 14 )  pt_av
     ENDIF
-    IF ( humidity  .OR. passive_scalar )  THEN
+    IF ( humidity )  THEN
        WRITE ( 14 )  'q                   ';  WRITE ( 14 )  q 
        IF ( ALLOCATED( q_av ) )  THEN
@@ -277 +277 @@
        WRITE ( 14 )  'qswst               ';  WRITE ( 14 ) qswst
     ENDIF
+    IF ( passive_scalar )  THEN
+       WRITE ( 14 )  's                   ';  WRITE ( 14 )  s 
+       IF ( ALLOCATED( s_av ) )  THEN
+          WRITE ( 14 )  's_av                ';  WRITE ( 14 )  s_av
+       ENDIF
+       WRITE ( 14 )  'ss                  ';  WRITE ( 14 )  ss
+       WRITE ( 14 )  'ssws                ';  WRITE ( 14 )  ssws
+       IF ( ALLOCATED( ssws_av ) )  THEN
+          WRITE ( 14 )  'ssws_av             ';  WRITE ( 14 )  ssws_av
+       ENDIF
+       WRITE ( 14 )  'sswst               ';  WRITE ( 14 ) sswst
+    ENDIF    
     IF ( land_surface )  THEN
        IF ( ALLOCATED( qsws_eb_av ) )  THEN

palm/trunk/SOURCE/write_var_list.f90

@@ -154 +154 @@
     USE arrays_3d,                                                             &
         ONLY:  inflow_damping_factor, mean_inflow_profiles, pt_init,           &
-               q_init, ref_state, sa_init, u_init, ug, v_init, vg
+               q_init, ref_state, s_init, sa_init, u_init, ug, v_init, vg
 
     USE control_parameters
@@ -504 +504 @@
     WRITE ( 14 )  'run_coupled                   '
     WRITE ( 14 )  run_coupled
+    WRITE ( 14 )  's_init                        '
+    WRITE ( 14 )  s_init
+    WRITE ( 14 )  's_surface                     '
+    WRITE ( 14 )  s_surface
+    WRITE ( 14 )  's_surface_initial_change      '
+    WRITE ( 14 )  s_surface_initial_change
+    WRITE ( 14 )  's_vertical_gradient           '
+    WRITE ( 14 )  s_vertical_gradient
+    WRITE ( 14 )  's_vertical_gradient_level     '
+    WRITE ( 14 )  s_vertical_gradient_level
+    WRITE ( 14 )  's_vertical_gradient_level_ind '
+    WRITE ( 14 )  s_vertical_gradient_level_ind 
     WRITE ( 14 )  'sa_init                       '
     WRITE ( 14 )  sa_init