Changeset 1960 for palm


Ignore:
Timestamp:
Jul 12, 2016 4:34:24 PM (5 years ago)
Author:
suehring
Message:

Separate balance equations for humidity and passive_scalar

Location:
palm/trunk/SOURCE
Files:
33 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/advec_s_bc.f90

    r1874 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! New CASE statement to treat scalars and humidity separately
    2222!
    2323! Former revisions:
     
    120120
    121121       USE control_parameters,                                                    &
    122            ONLY:  dt_3d, bc_pt_t_val, bc_q_t_val, ibc_pt_b, ibc_pt_t, ibc_q_t,    &
    123                   message_string, pt_slope_offset, sloping_surface, u_gtrans,     &
    124                   v_gtrans
     122           ONLY:  dt_3d, bc_pt_t_val, bc_q_t_val, bc_s_t_val, ibc_pt_b, ibc_pt_t, &
     123                  ibc_q_t, ibc_s_t, message_string, pt_slope_offset,              &
     124                  sloping_surface, u_gtrans, v_gtrans
    125125
    126126       USE cpulog,                                                                &
     
    969969          ENDIF
    970970
     971       ELSEIF ( sk_char == 's' )  THEN
     972!
     973!--       Specific scalar boundary condition at the bottom boundary.
     974!--       Dirichlet (fixed surface humidity) or Neumann (i.e. zero gradient)
     975          DO  i = nxl, nxr
     976             DO  j = nys, nyn
     977                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
     978                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     979                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     980             ENDDO
     981          ENDDO
     982
     983!
     984!--       Specific scalar boundary condition at the top boundary
     985          IF ( ibc_s_t == 0 )  THEN
     986!
     987!--          Dirichlet
     988             DO  i = nxl, nxr
     989                DO  j = nys, nyn
     990                   sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
     991                   sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
     992                ENDDO
     993             ENDDO
     994
     995          ELSE
     996!
     997!--          Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
     998             DO  i = nxl, nxr
     999                DO  j = nys, nyn
     1000                   sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i) + bc_s_t_val * dzu(nzt+1)
     1001                   sk_p(nzt+3,j,i)   = sk_p(nzt+2,j,i) + bc_s_t_val * dzu(nzt+1)
     1002                ENDDO
     1003             ENDDO
     1004
     1005          ENDIF
     1006
    9711007       ELSEIF ( sk_char == 'qr' )  THEN
    9721008
  • palm/trunk/SOURCE/advec_ws.f90

    r1943 r1960  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Separate humidity and passive scalar
    2222!
    2323! Former revisions:
     
    228228       USE arrays_3d,                                                          &
    229229           ONLY:  diss_l_e, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr,         &
    230                   diss_l_sa, diss_l_u, diss_l_v, diss_l_w,  flux_l_e,          &
    231                   flux_l_nr, flux_l_pt, flux_l_q, flux_l_qr, flux_l_sa,        &
    232                   flux_l_u, flux_l_v, flux_l_w, diss_s_e, diss_s_nr, diss_s_pt,&
    233                   diss_s_q, diss_s_qr, diss_s_sa, diss_s_u, diss_s_v, diss_s_w,&
    234                   flux_s_e, flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr,         &
    235                   flux_s_sa, flux_s_u, flux_s_v, flux_s_w
     230                  diss_l_s, diss_l_sa, diss_l_u, diss_l_v, diss_l_w,  flux_l_e,&
     231                  flux_l_nr, flux_l_pt, flux_l_q, flux_l_qr, flux_l_s,         &
     232                  flux_l_sa, flux_l_u, flux_l_v, flux_l_w, diss_s_e, diss_s_nr,&
     233                  diss_s_pt, diss_s_q, diss_s_qr, diss_s_s, diss_s_sa,         &
     234                  diss_s_u, diss_s_v, diss_s_w, flux_s_e, flux_s_nr, flux_s_pt,&
     235                  flux_s_q, flux_s_qr, flux_s_s, flux_s_sa, flux_s_u, flux_s_v,&
     236                  flux_s_w
    236237
    237238       USE constants,                                                          &
     
    254255           ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsnrs_ws_l,&
    255256                  sums_wspts_ws_l, sums_wsqrs_ws_l, sums_wsqs_ws_l,            &
    256                   sums_wssas_ws_l,  sums_wsus_ws_l, sums_wsvs_ws_l 
     257                  sums_wsss_ws_l, sums_wssas_ws_l,  sums_wsss_ws_l,               &
     258                  sums_wsus_ws_l, sums_wsvs_ws_l 
    257259
    258260!
     
    287289          sums_wspts_ws_l = 0.0_wp
    288290
    289           IF ( humidity  .OR.  passive_scalar )  THEN
     291          IF ( humidity  )  THEN
    290292             ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    291293             sums_wsqs_ws_l = 0.0_wp
     294          ENDIF
     295         
     296          IF ( passive_scalar )  THEN
     297             ALLOCATE( sums_wsss_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     298             sums_wsss_ws_l = 0.0_wp
    292299          ENDIF
    293300
     
    341348                       diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    342349
    343              IF ( humidity  .OR.  passive_scalar )  THEN
     350             IF ( humidity )  THEN
    344351                ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1),            &
    345352                          diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
     
    347354                          diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    348355             ENDIF
     356             
     357             IF ( passive_scalar )  THEN
     358                ALLOCATE( flux_s_s(nzb+1:nzt,0:threads_per_task-1),            &
     359                          diss_s_s(nzb+1:nzt,0:threads_per_task-1) )
     360                ALLOCATE( flux_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
     361                          diss_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
     362             ENDIF             
    349363
    350364             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     
    837851           ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsnrs_ws_l,&
    838852                  sums_wspts_ws_l, sums_wsqrs_ws_l, sums_wsqs_ws_l,            &
    839                   sums_wssas_ws_l,  sums_wsus_ws_l, sums_wsvs_ws_l 
     853                  sums_wsss_ws_l, sums_wssas_ws_l,  sums_wsus_ws_l,            &
     854                  sums_wsvs_ws_l 
    840855
    841856       IMPLICIT NONE
     
    854869       IF ( ws_scheme_sca )  THEN
    855870          sums_wspts_ws_l = 0.0_wp
    856           IF ( humidity  .OR.  passive_scalar )  sums_wsqs_ws_l = 0.0_wp
     871          IF ( humidity       )  sums_wsqs_ws_l = 0.0_wp
     872          IF ( passive_scalar )  sums_wsss_ws_l = 0.0_wp
    857873          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    858874             sums_wsqrs_ws_l = 0.0_wp
     
    898914       USE statistics,                                                         &
    899915           ONLY:  sums_wsnrs_ws_l, sums_wspts_ws_l, sums_wsqrs_ws_l,           &
    900                   sums_wsqs_ws_l, sums_wssas_ws_l, weight_substep
     916                  sums_wsqs_ws_l, sums_wssas_ws_l, sums_wsss_ws_l,             &
     917                  weight_substep
    901918
    902919       IMPLICIT NONE
     
    16651682                               * weight_substep(intermediate_timestep_count)
    16661683             ENDDO
     1684             
     1685          CASE ( 's' )
     1686         
     1687             DO  k = nzb, nzt
     1688                sums_wsss_ws_l(k,tn)  = sums_wsss_ws_l(k,tn) +                &
     1689                                      ( flux_t(k) + diss_t(k) )               &
     1690                               * weight_substep(intermediate_timestep_count)
     1691             ENDDO
    16671692
    16681693         END SELECT
     
    31203145       USE statistics,                                                        &
    31213146           ONLY:  sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,           &
    3122                   sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep
     3147                  sums_wsqrs_ws_l, sums_wsnrs_ws_l, sums_wsss_ws_l,           &
     3148                  weight_substep
    31233149
    31243150       IMPLICIT NONE
     
    38663892                        *   weight_substep(intermediate_timestep_count)
    38673893                    ENDDO
     3894                 CASE ( 's' )
     3895                    DO  k = nzb, nzt
     3896                       sums_wsss_ws_l(k,tn) = sums_wsss_ws_l(k,tn)          &
     3897                        + ( flux_t(k) + diss_t(k) )                           &
     3898                        *   weight_substep(intermediate_timestep_count)
     3899                    ENDDO   
     3900                                   
    38683901
    38693902              END SELECT
  • palm/trunk/SOURCE/average_3d_data.f90

    r1818 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Treat humidity and passive scalar separately
    2222!
    2323! Former revisions:
     
    498498             ENDDO
    499499
     500          CASE ( 'ssws*' )
     501             DO  i = nxlg, nxrg
     502                DO  j = nysg, nyng
     503                   ssws_av(j,i) = ssws_av(j,i) / REAL( average_count_3d, KIND=wp )
     504                ENDDO
     505             ENDDO
     506
    500507          CASE ( 't*' )
    501508             DO  i = nxlg, nxrg
  • palm/trunk/SOURCE/boundary_conds.f90

    r1933 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Treat humidity and passive scalar separately
    2222!
    2323! Former revisions:
     
    135135    USE arrays_3d,                                                             &
    136136        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,  &
    137                dzu, e_p, nr_p, pt, pt_p, q, q_p, qr_p, sa, sa_p,               &
     137               dzu, e_p, nr_p, pt, pt_p, q, q_p, qr_p, s, s_p, sa, sa_p,       &
    138138               u, ug, u_init, u_m_l, u_m_n, u_m_r, u_m_s, u_p,                 &
    139139               v, vg, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p,                 &
    140                w, w_p, w_m_l, w_m_n, w_m_r, w_m_s,&
    141                pt_init
     140               w, w_p, w_m_l, w_m_n, w_m_r, w_m_s, pt_init
    142141
    143142    USE control_parameters,                                                    &
    144         ONLY:  bc_pt_t_val, bc_q_t_val, constant_diffusion,                    &
     143        ONLY:  bc_pt_t_val, bc_q_t_val, bc_s_t_val, constant_diffusion,        &
    145144               cloud_physics, dt_3d, humidity,                                 &
    146                ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_sa_t, ibc_uv_b,       &
    147                ibc_uv_t, inflow_l, inflow_n, inflow_r, inflow_s,               &
    148                intermediate_timestep_count, large_scale_forcing,               &
     145               ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, ibc_s_t,         &
     146               ibc_sa_t, ibc_uv_b, ibc_uv_t, inflow_l, inflow_n, inflow_r,     &
     147               inflow_s, intermediate_timestep_count, large_scale_forcing,     &
    149148               microphysics_seifert, nest_domain, nest_bound_l, nest_bound_s,  &
    150149               nudging, ocean, outflow_l, outflow_n, outflow_r, outflow_s,     &
     
    302301
    303302!
    304 !-- Boundary conditions for total water content or scalar,
     303!-- Boundary conditions for total water content,
    305304!-- bottom and top boundary (see also temperature)
    306     IF ( humidity  .OR.  passive_scalar )  THEN
     305    IF ( humidity )  THEN
    307306!
    308307!--    Surface conditions for constant_humidity_flux
     
    344343       ENDIF
    345344    ENDIF
     345!
     346!-- Boundary conditions for scalar,
     347!-- bottom and top boundary (see also temperature)
     348    IF ( passive_scalar )  THEN
     349!
     350!--    Surface conditions for constant_humidity_flux
     351       IF ( ibc_s_b == 0 ) THEN
     352          DO  i = nxlg, nxrg
     353             DO  j = nysg, nyng
     354                s_p(nzb_s_inner(j,i),j,i) = s(nzb_s_inner(j,i),j,i)
     355             ENDDO
     356          ENDDO
     357       ELSE
     358          DO  i = nxlg, nxrg
     359             DO  j = nysg, nyng
     360                s_p(nzb_s_inner(j,i),j,i) = s_p(nzb_s_inner(j,i)+1,j,i)
     361             ENDDO
     362          ENDDO
     363       ENDIF
     364!
     365!--    Top boundary
     366       IF ( ibc_s_t == 0 ) THEN
     367          s_p(nzt+1,:,:) = s(nzt+1,:,:)
     368       ELSEIF ( ibc_s_t == 1 ) THEN
     369          s_p(nzt+1,:,:) = s_p(nzt,:,:)   + bc_s_t_val * dzu(nzt+1)
     370       ENDIF
     371
     372    ENDIF   
    346373!
    347374!-- In case of inflow or nest boundary at the south boundary the boundary for v
     
    381408       pt_p(:,nys-1,:)     = pt_p(:,nys,:)
    382409       IF ( .NOT. constant_diffusion     )  e_p(:,nys-1,:) = e_p(:,nys,:)
    383        IF ( humidity  .OR.  passive_scalar )  THEN
     410       IF ( humidity )  THEN
    384411          q_p(:,nys-1,:) = q_p(:,nys,:)
    385412          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     
    388415          ENDIF
    389416       ENDIF
     417       IF ( passive_scalar )  s_p(:,nys-1,:) = s_p(:,nys,:)
    390418    ELSEIF ( outflow_n )  THEN
    391419       pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
    392420       IF ( .NOT. constant_diffusion     )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
    393        IF ( humidity  .OR.  passive_scalar )  THEN
     421       IF ( humidity )  THEN
    394422          q_p(:,nyn+1,:) = q_p(:,nyn,:)
    395423          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     
    398426          ENDIF
    399427       ENDIF
     428       IF ( passive_scalar )  s_p(:,nyn+1,:) = s_p(:,nyn,:)
    400429    ELSEIF ( outflow_l )  THEN
    401430       pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
    402431       IF ( .NOT. constant_diffusion     )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
    403        IF ( humidity  .OR.  passive_scalar )  THEN
     432       IF ( humidity )  THEN
    404433          q_p(:,:,nxl-1) = q_p(:,:,nxl)
    405434          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     
    408437          ENDIF
    409438       ENDIF
     439       IF ( passive_scalar )  s_p(:,:,nxl-1) = s_p(:,:,nxl)
    410440    ELSEIF ( outflow_r )  THEN
    411441       pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
    412442       IF ( .NOT. constant_diffusion     )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
    413        IF ( humidity .OR. passive_scalar )  THEN
     443       IF ( humidity )  THEN
    414444          q_p(:,:,nxr+1) = q_p(:,:,nxr)
    415445          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     
    418448          ENDIF
    419449       ENDIF
     450       IF ( passive_scalar )  s_p(:,:,nxr+1) = s_p(:,:,nxr)
    420451    ENDIF
    421452
  • palm/trunk/SOURCE/check_parameters.f90

    r1932 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Separate checks and calculations for humidity and passive scalar. Therefore,
     22! a subroutine to calculate vertical gradients is introduced, in order  to reduce
     23! redundance.
     24! Additional check - large-scale forcing is not implemented for passive scalars
     25! so far.
    2226!
    2327! Former revisions:
     
    10531057    ENDIF
    10541058
    1055     IF ( passive_scalar  .AND.  humidity )  THEN
    1056        message_string = 'humidity = .TRUE. and passive_scalar = .TRUE. ' //    &
    1057                         'is not allowed simultaneously'
    1058        CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 )
    1059     ENDIF
    1060 
    10611059!
    10621060!-- When plant canopy model is used, peform addtional checks
     
    10941092!--    Initial profiles for 1D and 3D model, respectively (u,v further below)
    10951093       pt_init = pt_surface
    1096        IF ( humidity )  THEN
    1097           q_init  = q_surface
    1098        ENDIF
    1099        IF ( ocean )           sa_init = sa_surface
    1100        IF ( passive_scalar )  q_init  = s_surface
     1094       IF ( humidity       )  q_init  = q_surface
     1095       IF ( ocean          )  sa_init = sa_surface
     1096       IF ( passive_scalar )  s_init  = s_surface
    11011097
    11021098!
     
    12821278!--    Compute initial temperature profile using the given temperature gradients
    12831279       IF (  .NOT.  neutral )  THEN
    1284 
    1285           i = 1
    1286           gradient = 0.0_wp
    1287 
    1288           IF (  .NOT.  ocean )  THEN
    1289 
    1290              pt_vertical_gradient_level_ind(1) = 0
    1291              DO  k = 1, nzt+1
    1292                 IF ( i < 11 )  THEN
    1293                    IF ( pt_vertical_gradient_level(i) < zu(k)  .AND.           &
    1294                         pt_vertical_gradient_level(i) >= 0.0_wp )  THEN
    1295                       gradient = pt_vertical_gradient(i) / 100.0_wp
    1296                       pt_vertical_gradient_level_ind(i) = k - 1
    1297                       i = i + 1
    1298                    ENDIF
    1299                 ENDIF
    1300                 IF ( gradient /= 0.0_wp )  THEN
    1301                    IF ( k /= 1 )  THEN
    1302                       pt_init(k) = pt_init(k-1) + dzu(k) * gradient
    1303                    ELSE
    1304                       pt_init(k) = pt_surface   + dzu(k) * gradient
    1305                    ENDIF
    1306                 ELSE
    1307                    pt_init(k) = pt_init(k-1)
    1308                 ENDIF
    1309              ENDDO
    1310 
    1311           ELSE
    1312 
    1313              pt_vertical_gradient_level_ind(1) = nzt+1
    1314              DO  k = nzt, 0, -1
    1315                 IF ( i < 11 )  THEN
    1316                    IF ( pt_vertical_gradient_level(i) > zu(k)  .AND.           &
    1317                         pt_vertical_gradient_level(i) <= 0.0_wp )  THEN
    1318                       gradient = pt_vertical_gradient(i) / 100.0_wp
    1319                       pt_vertical_gradient_level_ind(i) = k + 1
    1320                       i = i + 1
    1321                    ENDIF
    1322                 ENDIF
    1323                 IF ( gradient /= 0.0_wp )  THEN
    1324                    IF ( k /= nzt )  THEN
    1325                       pt_init(k) = pt_init(k+1) - dzu(k+1) * gradient
    1326                    ELSE
    1327                       pt_init(k)   = pt_surface - 0.5_wp * dzu(k+1) * gradient
    1328                       pt_init(k+1) = pt_surface + 0.5_wp * dzu(k+1) * gradient
    1329                    ENDIF
    1330                 ELSE
    1331                    pt_init(k) = pt_init(k+1)
    1332                 ENDIF
    1333              ENDDO
    1334 
    1335           ENDIF
    1336 
    1337        ENDIF
    1338 
    1339 !
    1340 !--    In case of no given temperature gradients, choose gradient of neutral
    1341 !--    stratification
    1342        IF ( pt_vertical_gradient_level(1) == -9999999.9_wp )  THEN
    1343           pt_vertical_gradient_level(1) = 0.0_wp
    1344        ENDIF
    1345 
    1346 !
    1347 !--    Store temperature gradient at the top boundary for possible Neumann
    1348 !--    boundary condition
    1349        bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
    1350 
    1351 !
    1352 !--    If required, compute initial humidity or scalar profile using the given
    1353 !--    humidity/scalar gradient. In case of scalar transport, initially store
    1354 !--    values of the scalar parameters on humidity parameters
     1280          CALL init_vertical_profiles( pt_vertical_gradient_level_ind,          &
     1281                                       pt_vertical_gradient_level,              &
     1282                                       pt_vertical_gradient, pt_init,           &
     1283                                       pt_surface, bc_pt_t_val )
     1284       ENDIF
     1285!
     1286!--    Compute initial humidity profile using the given humidity gradients
     1287       IF ( humidity )  THEN
     1288          CALL init_vertical_profiles( q_vertical_gradient_level_ind,          &
     1289                                       q_vertical_gradient_level,              &
     1290                                       q_vertical_gradient, q_init,            &
     1291                                       q_surface, bc_q_t_val )
     1292       ENDIF
     1293!
     1294!--    Compute initial scalar profile using the given scalar gradients
    13551295       IF ( passive_scalar )  THEN
    1356           bc_q_b                    = bc_s_b
    1357           bc_q_t                    = bc_s_t
    1358           q_surface                 = s_surface
    1359           q_surface_initial_change  = s_surface_initial_change
    1360           q_vertical_gradient       = s_vertical_gradient
    1361           q_vertical_gradient_level = s_vertical_gradient_level
    1362           surface_waterflux         = surface_scalarflux
    1363           wall_humidityflux         = wall_scalarflux
    1364        ENDIF
    1365 
    1366        IF ( humidity  .OR.  passive_scalar )  THEN
    1367 
    1368           i = 1
    1369           gradient = 0.0_wp
    1370 
    1371           IF (  .NOT.  ocean )  THEN
    1372 
    1373              q_vertical_gradient_level_ind(1) = 0
    1374              DO  k = 1, nzt+1
    1375                 IF ( i < 11 )  THEN
    1376                    IF ( q_vertical_gradient_level(i) < zu(k)  .AND.            &
    1377                         q_vertical_gradient_level(i) >= 0.0_wp )  THEN
    1378                       gradient = q_vertical_gradient(i) / 100.0_wp
    1379                       q_vertical_gradient_level_ind(i) = k - 1
    1380                       i = i + 1
    1381                    ENDIF
    1382                 ENDIF
    1383                 IF ( gradient /= 0.0_wp )  THEN
    1384                    IF ( k /= 1 )  THEN
    1385                       q_init(k) = q_init(k-1) + dzu(k) * gradient
    1386                    ELSE
    1387                       q_init(k) = q_init(k-1) + dzu(k) * gradient
    1388                    ENDIF
    1389                 ELSE
    1390                    q_init(k) = q_init(k-1)
    1391                 ENDIF
    1392    !
    1393    !--          Avoid negative humidities
    1394                 IF ( q_init(k) < 0.0_wp )  THEN
    1395                    q_init(k) = 0.0_wp
    1396                 ENDIF
    1397              ENDDO
    1398 
    1399           ELSE
    1400 
    1401              q_vertical_gradient_level_ind(1) = nzt+1
    1402              DO  k = nzt, 0, -1
    1403                 IF ( i < 11 )  THEN
    1404                    IF ( q_vertical_gradient_level(i) > zu(k)  .AND.            &
    1405                         q_vertical_gradient_level(i) <= 0.0_wp )  THEN
    1406                       gradient = q_vertical_gradient(i) / 100.0_wp
    1407                       q_vertical_gradient_level_ind(i) = k + 1
    1408                       i = i + 1
    1409                    ENDIF
    1410                 ENDIF
    1411                 IF ( gradient /= 0.0_wp )  THEN
    1412                    IF ( k /= nzt )  THEN
    1413                       q_init(k) = q_init(k+1) - dzu(k+1) * gradient
    1414                    ELSE
    1415                       q_init(k)   = q_surface - 0.5_wp * dzu(k+1) * gradient
    1416                       q_init(k+1) = q_surface + 0.5_wp * dzu(k+1) * gradient
    1417                    ENDIF
    1418                 ELSE
    1419                    q_init(k) = q_init(k+1)
    1420                 ENDIF
    1421    !
    1422    !--          Avoid negative humidities
    1423                 IF ( q_init(k) < 0.0_wp )  THEN
    1424                    q_init(k) = 0.0_wp
    1425                 ENDIF
    1426              ENDDO
    1427 
    1428           ENDIF
    1429 !
    1430 !--       In case of no given humidity gradients, choose zero gradient
    1431 !--       conditions
    1432           IF ( q_vertical_gradient_level(1) == -1.0_wp )  THEN
    1433              q_vertical_gradient_level(1) = 0.0_wp
    1434           ENDIF
    1435 !
    1436 !--       Store humidity, rain water content and rain drop concentration
    1437 !--       gradient at the top boundary for possile Neumann boundary condition
    1438           bc_q_t_val  = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
    1439        ENDIF
    1440 
     1296          CALL init_vertical_profiles( s_vertical_gradient_level_ind,          &
     1297                                       s_vertical_gradient_level,              &
     1298                                       s_vertical_gradient, s_init,            &
     1299                                       s_surface, bc_s_t_val )
     1300       ENDIF
    14411301!
    14421302!--    If required, compute initial salinity profile using the given salinity
    14431303!--    gradients
    14441304       IF ( ocean )  THEN
    1445 
    1446           i = 1
    1447           gradient = 0.0_wp
    1448 
    1449           sa_vertical_gradient_level_ind(1) = nzt+1
    1450           DO  k = nzt, 0, -1
    1451              IF ( i < 11 )  THEN
    1452                 IF ( sa_vertical_gradient_level(i) > zu(k)  .AND.              &
    1453                      sa_vertical_gradient_level(i) <= 0.0_wp )  THEN
    1454                    gradient = sa_vertical_gradient(i) / 100.0_wp
    1455                    sa_vertical_gradient_level_ind(i) = k + 1
    1456                    i = i + 1
    1457                 ENDIF
    1458              ENDIF
    1459              IF ( gradient /= 0.0_wp )  THEN
    1460                 IF ( k /= nzt )  THEN
    1461                    sa_init(k) = sa_init(k+1) - dzu(k+1) * gradient
    1462                 ELSE
    1463                    sa_init(k)   = sa_surface - 0.5_wp * dzu(k+1) * gradient
    1464                    sa_init(k+1) = sa_surface + 0.5_wp * dzu(k+1) * gradient
    1465                 ENDIF
    1466              ELSE
    1467                 sa_init(k) = sa_init(k+1)
    1468              ENDIF
    1469           ENDDO
    1470 
     1305          CALL init_vertical_profiles( sa_vertical_gradient_level_ind,          &
     1306                                       sa_vertical_gradient_level,              &
     1307                                       sa_vertical_gradient, s_init,            &
     1308                                       sa_surface, -999.0_wp )
    14711309       ENDIF
    14721310
     
    18601698
    18611699!
    1862 !-- In case of humidity or passive scalar, set boundary conditions for total
    1863 !-- water content / scalar
    1864     IF ( humidity  .OR.  passive_scalar ) THEN
    1865        IF ( humidity )  THEN
    1866           sq = 'q'
    1867        ELSE
    1868           sq = 's'
    1869        ENDIF
    1870        IF ( bc_q_b == 'dirichlet' )  THEN
    1871           ibc_q_b = 0
    1872        ELSEIF ( bc_q_b == 'neumann' )  THEN
    1873           ibc_q_b = 1
    1874        ELSE
    1875           message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
    1876                            '_b ="' // TRIM( bc_q_b ) // '"'
    1877           CALL message( 'check_parameters', 'PA0071', 1, 2, 0, 6, 0 )
    1878        ENDIF
    1879        IF ( bc_q_t == 'dirichlet' )  THEN
    1880           ibc_q_t = 0
    1881        ELSEIF ( bc_q_t == 'neumann' )  THEN
    1882           ibc_q_t = 1
    1883        ELSEIF ( bc_q_t == 'nested' )  THEN
    1884           ibc_q_t = 3
    1885        ELSE
    1886           message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
    1887                            '_t ="' // TRIM( bc_q_t ) // '"'
    1888           CALL message( 'check_parameters', 'PA0072', 1, 2, 0, 6, 0 )
    1889        ENDIF
     1700!-- Set boundary conditions for total water content
     1701    IF ( humidity )  THEN
     1702       CALL check_bc_scalars( 'q', bc_q_b, bc_q_t, ibc_q_b, ibc_q_t,           &
     1703                              'PA0071', 'PA0072', 'PA0073', 'PA0074',          &
     1704                              surface_waterflux, constant_waterflux,           &
     1705                              q_surface_initial_change )
    18901706
    18911707       IF ( surface_waterflux == 9999999.9_wp  )  THEN
     
    19091725       ENDIF
    19101726
    1911 !
    1912 !--    A given surface humidity implies Dirichlet boundary condition for
    1913 !--    humidity. In this case specification of a constant water flux is
    1914 !--    forbidden.
    1915        IF ( ibc_q_b == 0  .AND.  constant_waterflux )  THEN
    1916           message_string = 'boundary condition: bc_' // TRIM( sq ) // '_b ' // &
    1917                            '= "' // TRIM( bc_q_b ) // '" is not allowed wi' // &
    1918                            'th prescribed surface flux'
    1919           CALL message( 'check_parameters', 'PA0073', 1, 2, 0, 6, 0 )
    1920        ENDIF
    1921        IF ( constant_waterflux  .AND.  q_surface_initial_change /= 0.0_wp )  THEN
    1922           WRITE( message_string, * )  'a prescribed surface flux is not allo', &
    1923                  'wed with ', sq, '_surface_initial_change (/=0) = ',          &
    1924                  q_surface_initial_change
    1925           CALL message( 'check_parameters', 'PA0074', 1, 2, 0, 6, 0 )
    1926        ENDIF
    1927 
     1727    ENDIF
     1728   
     1729    IF ( passive_scalar )  THEN
     1730       CALL check_bc_scalars( 's', bc_s_b, bc_s_t, ibc_s_b, ibc_s_t,           &
     1731                              'PA0071', 'PA0072', 'PA0073', 'PA0074',          &
     1732                              surface_scalarflux, constant_scalarflux,         &
     1733                              s_surface_initial_change )
    19281734    ENDIF
    19291735!
     
    38503656                        'the usage of large scale forcing from external file.'
    38513657       CALL message( 'check_parameters', 'PA0375', 1, 2, 0, 6, 0 )
    3852      ENDIF
     3658    ENDIF
    38533659
    38543660    IF ( large_scale_forcing  .AND.  (  .NOT.  humidity ) )  THEN
     
    38563662                        'file LSF_DATA requires humidity = .T..'
    38573663       CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 )
    3858      ENDIF
     3664    ENDIF
     3665
     3666    IF ( large_scale_forcing  .AND.  passive_scalar )  THEN
     3667       message_string = 'The usage of large scale forcing from external &'//   &
     3668                        'file LSF_DATA is not implemented for psassive scalars'
     3669       CALL message( 'check_parameters', 'PA0440', 1, 2, 0, 6, 0 )
     3670    ENDIF
    38593671
    38603672    IF ( large_scale_forcing  .AND.  topography /= 'flat' )  THEN
     
    38643676    ENDIF
    38653677
    3866     IF ( large_scale_forcing  .AND.  ocean  )  THEN
     3678    IF ( large_scale_forcing  .AND.  ocean )  THEN
    38673679       message_string = 'The usage of large scale forcing from external &'//   &
    38683680                        'file LSF_DATA is not implemented for ocean runs'
     
    39433755
    39443756    END SUBROUTINE check_dt_do
     3757   
     3758   
     3759   
     3760   
     3761!------------------------------------------------------------------------------!
     3762! Description:
     3763! ------------
     3764!> Inititalizes the vertical profiles of scalar quantities.
     3765!------------------------------------------------------------------------------!
     3766
     3767    SUBROUTINE init_vertical_profiles( vertical_gradient_level_ind,            &
     3768                                       vertical_gradient_level,                &
     3769                                       vertical_gradient,                      &
     3770                                       pr_init, surface_value, bc_t_val )
     3771
     3772
     3773       IMPLICIT NONE   
     3774   
     3775       INTEGER(iwp) ::  i     !< counter
     3776       INTEGER(iwp), DIMENSION(1:10) ::  vertical_gradient_level_ind !< vertical grid indices for gradient levels
     3777       
     3778       REAL(wp)     ::  bc_t_val      !< model top gradient       
     3779       REAL(wp)     ::  gradient      !< vertica gradient of the respective quantity
     3780       REAL(wp)     ::  surface_value !< surface value of the respecitve quantity
     3781
     3782       
     3783       REAL(wp), DIMENSION(0:nz+1) ::  pr_init                 !< initialisation profile
     3784       REAL(wp), DIMENSION(1:10)   ::  vertical_gradient       !< given vertical gradient
     3785       REAL(wp), DIMENSION(1:10)   ::  vertical_gradient_level !< given vertical gradient level
     3786       
     3787       i = 1
     3788       gradient = 0.0_wp
     3789
     3790       IF (  .NOT.  ocean )  THEN
     3791
     3792          vertical_gradient_level_ind(1) = 0
     3793          DO  k = 1, nzt+1
     3794             IF ( i < 11 )  THEN
     3795                IF ( vertical_gradient_level(i) < zu(k)  .AND.            &
     3796                     vertical_gradient_level(i) >= 0.0_wp )  THEN
     3797                   gradient = vertical_gradient(i) / 100.0_wp
     3798                   vertical_gradient_level_ind(i) = k - 1
     3799                   i = i + 1
     3800                ENDIF
     3801             ENDIF
     3802             IF ( gradient /= 0.0_wp )  THEN
     3803                IF ( k /= 1 )  THEN
     3804                   pr_init(k) = pr_init(k-1) + dzu(k) * gradient
     3805                ELSE
     3806                   pr_init(k) = pr_init(k-1) + dzu(k) * gradient
     3807                ENDIF
     3808             ELSE
     3809                pr_init(k) = pr_init(k-1)
     3810             ENDIF
     3811   !
     3812   !--       Avoid negative values
     3813             IF ( pr_init(k) < 0.0_wp )  THEN
     3814                pr_init(k) = 0.0_wp
     3815             ENDIF
     3816          ENDDO
     3817
     3818       ELSE
     3819
     3820          vertical_gradient_level_ind(1) = nzt+1
     3821          DO  k = nzt, 0, -1
     3822             IF ( i < 11 )  THEN
     3823                IF ( vertical_gradient_level(i) > zu(k)  .AND.            &
     3824                     vertical_gradient_level(i) <= 0.0_wp )  THEN
     3825                   gradient = vertical_gradient(i) / 100.0_wp
     3826                   vertical_gradient_level_ind(i) = k + 1
     3827                   i = i + 1
     3828                ENDIF
     3829             ENDIF
     3830             IF ( gradient /= 0.0_wp )  THEN
     3831                IF ( k /= nzt )  THEN
     3832                   pr_init(k) = pr_init(k+1) - dzu(k+1) * gradient
     3833                ELSE
     3834                   pr_init(k)   = surface_value - 0.5_wp * dzu(k+1) * gradient
     3835                   pr_init(k+1) = surface_value + 0.5_wp * dzu(k+1) * gradient
     3836                ENDIF
     3837             ELSE
     3838                pr_init(k) = pr_init(k+1)
     3839             ENDIF
     3840   !
     3841   !--       Avoid negative humidities
     3842             IF ( pr_init(k) < 0.0_wp )  THEN
     3843                pr_init(k) = 0.0_wp
     3844             ENDIF
     3845          ENDDO
     3846
     3847       ENDIF
     3848       
     3849!
     3850!--    In case of no given gradients, choose zero gradient conditions
     3851       IF ( vertical_gradient_level(1) == -1.0_wp )  THEN
     3852          vertical_gradient_level(1) = 0.0_wp
     3853       ENDIF
     3854!
     3855!--    Store gradient at the top boundary for possible Neumann boundary condition
     3856       bc_t_val  = ( pr_init(nzt+1) - pr_init(nzt) ) / dzu(nzt+1)
     3857   
     3858    END SUBROUTINE init_vertical_profiles
     3859   
     3860   
     3861     
     3862!------------------------------------------------------------------------------!
     3863! Description:
     3864! ------------
     3865!> Check the bottom and top boundary conditions for humidity and scalars.
     3866!------------------------------------------------------------------------------!
     3867
     3868    SUBROUTINE check_bc_scalars( sq, bc_b, bc_t, ibc_b, ibc_t,                 &
     3869                                 err_nr_b, err_nr_t, err_nr_3, err_nr_4,       &
     3870                                 surface_flux, constant_flux, surface_initial_change )
     3871
     3872
     3873       IMPLICIT NONE   
     3874   
     3875       CHARACTER (LEN=1)   ::  sq                       !<
     3876       CHARACTER (LEN=*)   ::  bc_b
     3877       CHARACTER (LEN=*)   ::  bc_t
     3878       CHARACTER (LEN=*)   ::  err_nr_b
     3879       CHARACTER (LEN=*)   ::  err_nr_t
     3880       CHARACTER (LEN=*)   ::  err_nr_3
     3881       CHARACTER (LEN=*)   ::  err_nr_4
     3882       
     3883       INTEGER(iwp)        ::  ibc_b
     3884       INTEGER(iwp)        ::  ibc_t
     3885       
     3886       LOGICAL             ::  constant_flux
     3887       
     3888       REAL(wp)            ::  surface_flux
     3889       REAL(wp)            ::  surface_initial_change
     3890       
     3891
     3892       IF ( bc_b == 'dirichlet' )  THEN
     3893          ibc_b = 0
     3894       ELSEIF ( bc_b == 'neumann' )  THEN
     3895          ibc_b = 1
     3896       ELSE
     3897          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
     3898                           '_b ="' // TRIM( bc_b ) // '"'
     3899          CALL message( 'check_parameters', err_nr_b, 1, 2, 0, 6, 0 )
     3900       ENDIF
     3901       
     3902       IF ( bc_t == 'dirichlet' )  THEN
     3903          ibc_t = 0
     3904       ELSEIF ( bc_t == 'neumann' )  THEN
     3905          ibc_t = 1
     3906       ELSEIF ( bc_t == 'nested' )  THEN
     3907          ibc_t = 3
     3908       ELSE
     3909          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
     3910                           '_t ="' // TRIM( bc_t ) // '"'
     3911          CALL message( 'check_parameters', err_nr_t, 1, 2, 0, 6, 0 )
     3912       ENDIF
     3913 
     3914!
     3915!--    A given surface value implies Dirichlet boundary condition for
     3916!--    the respective quantity. In this case specification of a constant flux is
     3917!--    forbidden.
     3918       IF ( ibc_b == 0  .AND.  constant_flux )  THEN
     3919          message_string = 'boundary condition: bc_' // TRIM( sq ) // '_b ' // &
     3920                           '= "' // TRIM( bc_b ) // '" is not allowed wi' // &
     3921                           'th prescribed surface flux'
     3922          CALL message( 'check_parameters', err_nr_3, 1, 2, 0, 6, 0 )
     3923       ENDIF
     3924       IF ( constant_waterflux  .AND.  surface_initial_change /= 0.0_wp )  THEN
     3925          WRITE( message_string, * )  'a prescribed surface flux is not allo', &
     3926                 'wed with ', sq, '_surface_initial_change (/=0) = ',          &
     3927                 surface_initial_change
     3928          CALL message( 'check_parameters', err_nr_4, 1, 2, 0, 6, 0 )
     3929       ENDIF 
     3930       
     3931   
     3932    END SUBROUTINE check_bc_scalars   
     3933   
     3934   
    39453935
    39463936 END SUBROUTINE check_parameters
  • palm/trunk/SOURCE/data_output_2d.f90

    r1852 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Scalar surface flux added
     22! Rename INTEGER variable s into s_ind, as s is already assigned to scalar
    2223!
    2324! Former revisions:
     
    138139    USE arrays_3d,                                                             &
    139140        ONLY:  dzw, e, nr, ol, p, pt, precipitation_amount, precipitation_rate,&
    140                prr,q, qc, ql, ql_c, ql_v, ql_vp, qr, qsws, rho, sa, shf, tend, &
    141                ts, u, us, v, vpt, w, z0, z0h, z0q, zu, zw
     141               prr,q, qc, ql, ql_c, ql_v, ql_vp, qr, qsws, rho, s, sa, shf,    &
     142               ssws, tend, ts, u, us, v, vpt, w, z0, z0h, z0q, zu, zw
    142143       
    143144    USE averaging
     
    224225    INTEGER(iwp) ::  nzt_do    !< upper limit of the data field (usually nzt+1)
    225226    INTEGER(iwp) ::  psi       !<
    226     INTEGER(iwp) ::  s         !<
     227    INTEGER(iwp) ::  s_ind     !<
    227228    INTEGER(iwp) ::  sender    !<
    228229    INTEGER(iwp) ::  ind(4)    !<
     
    268269
    269270       CASE ( 'xy' )
    270           s = 1
     271          s_ind = 1
    271272          ALLOCATE( level_z(nzb:nzt+1), local_2d(nxlg:nxrg,nysg:nyng) )
    272273
    273274          IF ( netcdf_data_format > 4 )  THEN
    274275             ns = 1
    275              DO WHILE ( section(ns,s) /= -9999  .AND.  ns <= 100 )
     276             DO WHILE ( section(ns,s_ind) /= -9999  .AND.  ns <= 100 )
    276277                ns = ns + 1
    277278             ENDDO
     
    298299
    299300       CASE ( 'xz' )
    300           s = 2
     301          s_ind = 2
    301302          ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) )
    302303
    303304          IF ( netcdf_data_format > 4 )  THEN
    304305             ns = 1
    305              DO WHILE ( section(ns,s) /= -9999  .AND.  ns <= 100 )
     306             DO WHILE ( section(ns,s_ind) /= -9999  .AND.  ns <= 100 )
    306307                ns = ns + 1
    307308             ENDDO
     
    329330
    330331       CASE ( 'yz' )
    331           s = 3
     332          s_ind = 3
    332333          ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) )
    333334
    334335          IF ( netcdf_data_format > 4 )  THEN
    335336             ns = 1
    336              DO WHILE ( section(ns,s) /= -9999  .AND.  ns <= 100 )
     337             DO WHILE ( section(ns,s_ind) /= -9999  .AND.  ns <= 100 )
    337338                ns = ns + 1
    338339             ENDDO
     
    10701071             CASE ( 's_xy', 's_xz', 's_yz' )
    10711072                IF ( av == 0 )  THEN
    1072                    to_be_resorted => q
     1073                   to_be_resorted => s
    10731074                ELSE
    10741075                   to_be_resorted => s_av
     
    11171118                two_d = .TRUE.
    11181119                level_z(nzb+1) = zu(nzb+1)
     1120               
     1121             CASE ( 'ssws*_xy' )        ! 2d-array
     1122                IF ( av == 0 ) THEN
     1123                   DO  i = nxlg, nxrg
     1124                      DO  j = nysg, nyng
     1125                         local_pf(i,j,nzb+1) =  ssws(j,i)
     1126                      ENDDO
     1127                   ENDDO
     1128                ELSE
     1129                   DO  i = nxlg, nxrg
     1130                      DO  j = nysg, nyng
     1131                         local_pf(i,j,nzb+1) =  ssws_av(j,i)
     1132                      ENDDO
     1133                   ENDDO
     1134                ENDIF
     1135                resorted = .TRUE.
     1136                two_d = .TRUE.
     1137                level_z(nzb+1) = zu(nzb+1)               
    11191138
    11201139             CASE ( 't*_xy' )        ! 2d-array
     
    13031322!--       section mode chosen.
    13041323          is = 1
    1305    loop1: DO WHILE ( section(is,s) /= -9999  .OR.  two_d )
     1324   loop1: DO WHILE ( section(is,s_ind) /= -9999  .OR.  two_d )
    13061325
    13071326             SELECT CASE ( mode )
     
    13131332                      layer_xy = nzb+1
    13141333                   ELSE
    1315                       layer_xy = section(is,s)
     1334                      layer_xy = section(is,s_ind)
    13161335                   ENDIF
    13171336
     
    13471366!
    13481367!--                If required, carry out averaging along z
    1349                    IF ( section(is,s) == -1  .AND.  .NOT. two_d )  THEN
     1368                   IF ( section(is,s_ind) == -1  .AND.  .NOT. two_d )  THEN
    13501369
    13511370                      local_2d = 0.0_wp
     
    15371556!
    15381557!--                If required, carry out averaging along y
    1539                    IF ( section(is,s) == -1 )  THEN
     1558                   IF ( section(is,s_ind) == -1 )  THEN
    15401559
    15411560                      ALLOCATE( local_2d_l(nxlg:nxrg,nzb_do:nzt_do) )
     
    15701589!--                   Just store the respective section on the local array
    15711590!--                   (but only if it is available on this PE!)
    1572                       IF ( section(is,s) >= nys  .AND.  section(is,s) <= nyn ) &
     1591                      IF ( section(is,s_ind) >= nys  .AND.  section(is,s_ind) <= nyn ) &
    15731592                      THEN
    1574                          local_2d = local_pf(:,section(is,s),nzb_do:nzt_do)
     1593                         local_2d = local_pf(:,section(is,s_ind),nzb_do:nzt_do)
    15751594                      ENDIF
    15761595
     
    15841603!--                   sections reside. Cross sections averaged along y are
    15851604!--                   output on the respective first PE along y (myidy=0).
    1586                       IF ( ( section(is,s) >= nys  .AND.                       &
    1587                              section(is,s) <= nyn )  .OR.                      &
    1588                            ( section(is,s) == -1  .AND.  myidy == 0 ) )  THEN
     1605                      IF ( ( section(is,s_ind) >= nys  .AND.                   &
     1606                             section(is,s_ind) <= nyn )  .OR.                  &
     1607                           ( section(is,s_ind) == -1  .AND.  myidy == 0 ) )  THEN
    15891608#if defined( __netcdf )
    15901609!
     
    16151634                         DO  i = 0, io_blocks-1
    16161635                            IF ( i == io_group )  THEN
    1617                                IF ( ( section(is,s) >= nys  .AND.              &
    1618                                       section(is,s) <= nyn )  .OR.             &
    1619                                     ( section(is,s) == -1  .AND.               &
     1636                               IF ( ( section(is,s_ind) >= nys  .AND.          &
     1637                                      section(is,s_ind) <= nyn )  .OR.         &
     1638                                    ( section(is,s_ind) == -1  .AND.           &
    16201639                                      nys-1 == -1 ) )                          &
    16211640                               THEN
     
    16431662!
    16441663!--                         Local array can be relocated directly.
    1645                             IF ( ( section(is,s) >= nys  .AND.                 &
    1646                                    section(is,s) <= nyn )  .OR.                &
    1647                                  ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
    1648                             THEN
     1664                            IF ( ( section(is,s_ind) >= nys  .AND.              &
     1665                                   section(is,s_ind) <= nyn )  .OR.             &
     1666                                 ( section(is,s_ind) == -1  .AND.              &
     1667                                   nys-1 == -1 ) )  THEN
    16491668                               total_2d(nxlg:nxrg,nzb_do:nzt_do) = local_2d
    16501669                            ENDIF
     
    16911710!--                         If the cross section resides on the PE, send the
    16921711!--                         local index limits, otherwise send -9999 to PE0.
    1693                             IF ( ( section(is,s) >= nys  .AND.                 &
    1694                                    section(is,s) <= nyn )  .OR.                &
    1695                                  ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
     1712                            IF ( ( section(is,s_ind) >= nys  .AND.              &
     1713                                   section(is,s_ind) <= nyn )  .OR.             &
     1714                                 ( section(is,s_ind) == -1  .AND.  nys-1 == -1 ) ) &
    16961715                            THEN
    16971716                               ind(1) = nxlg; ind(2) = nxrg
     
    17561775!
    17571776!--                If required, carry out averaging along x
    1758                    IF ( section(is,s) == -1 )  THEN
     1777                   IF ( section(is,s_ind) == -1 )  THEN
    17591778
    17601779                      ALLOCATE( local_2d_l(nysg:nyng,nzb_do:nzt_do) )
     
    17891808!--                   Just store the respective section on the local array
    17901809!--                   (but only if it is available on this PE!)
    1791                       IF ( section(is,s) >= nxl  .AND.  section(is,s) <= nxr ) &
     1810                      IF ( section(is,s_ind) >= nxl  .AND.  section(is,s_ind) <= nxr ) &
    17921811                      THEN
    1793                          local_2d = local_pf(section(is,s),:,nzb_do:nzt_do)
     1812                         local_2d = local_pf(section(is,s_ind),:,nzb_do:nzt_do)
    17941813                      ENDIF
    17951814
     
    18031822!--                   sections reside. Cross sections averaged along x are
    18041823!--                   output on the respective first PE along x (myidx=0).
    1805                       IF ( ( section(is,s) >= nxl  .AND.                       &
    1806                              section(is,s) <= nxr )  .OR.                      &
    1807                            ( section(is,s) == -1  .AND.  myidx == 0 ) )  THEN
     1824                      IF ( ( section(is,s_ind) >= nxl  .AND.                       &
     1825                             section(is,s_ind) <= nxr )  .OR.                      &
     1826                           ( section(is,s_ind) == -1  .AND.  myidx == 0 ) )  THEN
    18081827#if defined( __netcdf )
    18091828!
     
    18341853                         DO  i = 0, io_blocks-1
    18351854                            IF ( i == io_group )  THEN
    1836                                IF ( ( section(is,s) >= nxl  .AND.              &
    1837                                       section(is,s) <= nxr )  .OR.             &
    1838                                     ( section(is,s) == -1  .AND.               &
     1855                               IF ( ( section(is,s_ind) >= nxl  .AND.          &
     1856                                      section(is,s_ind) <= nxr )  .OR.         &
     1857                                    ( section(is,s_ind) == -1  .AND.           &
    18391858                                      nxl-1 == -1 ) )                          &
    18401859                               THEN
     
    18621881!
    18631882!--                         Local array can be relocated directly.
    1864                             IF ( ( section(is,s) >= nxl  .AND.                 &
    1865                                    section(is,s) <= nxr )   .OR.               &
    1866                                  ( section(is,s) == -1  .AND.  nxl-1 == -1 ) ) &
     1883                            IF ( ( section(is,s_ind) >= nxl  .AND.             &
     1884                                   section(is,s_ind) <= nxr )   .OR.           &
     1885                                 ( section(is,s_ind) == -1  .AND.  nxl-1 == -1 ) ) &
    18671886                            THEN
    18681887                               total_2d(nysg:nyng,nzb_do:nzt_do) = local_2d
     
    19101929!--                         If the cross section resides on the PE, send the
    19111930!--                         local index limits, otherwise send -9999 to PE0.
    1912                             IF ( ( section(is,s) >= nxl  .AND.                 &
    1913                                    section(is,s) <= nxr )  .OR.                &
    1914                                  ( section(is,s) == -1  .AND.  nxl-1 == -1 ) ) &
     1931                            IF ( ( section(is,s_ind) >= nxl  .AND.              &
     1932                                   section(is,s_ind) <= nxr )  .OR.             &
     1933                                 ( section(is,s_ind) == -1  .AND.  nxl-1 == -1 ) ) &
    19151934                            THEN
    19161935                               ind(1) = nysg; ind(2) = nyng
     
    21482167!
    21492168!-- Close plot output file.
    2150     file_id = 20 + s
     2169    file_id = 20 + s_ind
    21512170
    21522171    IF ( data_output_2d_on_each_pe )  THEN
  • palm/trunk/SOURCE/data_output_3d.f90

    r1852 r1960  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Scalar surface flux added
    2222!
    2323! Former revisions:
     
    119119
    120120    USE arrays_3d,                                                             &
    121         ONLY:  e, nr, p, pt, prr, q, qc, ql, ql_c, ql_v, qr, rho, sa, tend, u, &
    122                v, vpt, w
     121        ONLY:  e, nr, p, pt, prr, q, qc, ql, ql_c, ql_v, qr, rho, s, sa, tend, &
     122               u, v, vpt, w
    123123       
    124124    USE averaging
     
    566566          CASE ( 's' )
    567567             IF ( av == 0 )  THEN
    568                 to_be_resorted => q
     568                to_be_resorted => s
    569569             ELSE
    570570                to_be_resorted => s_av
  • palm/trunk/SOURCE/data_output_dvrp.f90

    r1874 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Separate humidity and passive scalar
    2222!
    2323! Former revisions:
     
    112112
    113113    USE arrays_3d,                                                             &
    114         ONLY:  p, pt, q, ql, ts, u, us, v, w, zu
     114        ONLY:  p, pt, q, ql, s, ts, u, us, v, w, zu
    115115       
    116116    USE cloud_parameters,                                                      &
     
    336336
    337337                CASE ( 'q', 'q_xy', 'q_xz', 'q_yz' )
    338                    IF ( humidity  .OR.  passive_scalar )  THEN
     338                   IF ( humidity )  THEN
    339339                      DO  i = nxl_dvrp, nxr_dvrp+1
    340340                         DO  j = nys_dvrp, nyn_dvrp+1
     
    345345                      ENDDO           
    346346                   ELSE                   
    347                       message_string = 'if humidity/passive_scalar = '    //   &
     347                      message_string = 'if humidity = '    //                  &
    348348                            'FALSE output of ' // TRIM( output_variable ) //   &
    349349                            'is not provided'
     
    366366                                  'is not provided'
    367367                      CALL message( 'data_output_dvrp', 'PA0184',&
     368                                                                 0, 0, 0, 6, 0 )
     369                   ENDIF
     370
     371                CASE ( 's', 's_xy', 's_xz', 's_yz' )
     372                   IF ( passive_scalar )  THEN
     373                      DO  i = nxl_dvrp, nxr_dvrp+1
     374                         DO  j = nys_dvrp, nyn_dvrp+1
     375                            DO  k = nzb, nz_do3d
     376                               local_pf(i,j,k) = s(k,j,i)
     377                            ENDDO
     378                         ENDDO
     379                      ENDDO           
     380                   ELSE                   
     381                      message_string = 'if passive_scalar = '    //            &
     382                            'FALSE output of ' // TRIM( output_variable ) //   &
     383                            'is not provided'
     384                      CALL message( 'data_output_dvrp', 'PA0183',&
    368385                                                                 0, 0, 0, 6, 0 )
    369386                   ENDIF
  • palm/trunk/SOURCE/data_output_mask.f90

    r1818 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Separate humidity and passive scalar
    2222!
    2323! Former revisions:
     
    8888#if defined( __netcdf )
    8989    USE arrays_3d,                                                             &
    90         ONLY:  e, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho, sa, tend, u,      &
     90        ONLY:  e, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho, s, sa, tend, u,   &
    9191               v, vpt, w
    9292   
     
    487487          CASE ( 's' )
    488488             IF ( av == 0 )  THEN
    489                 to_be_resorted => q
     489                to_be_resorted => s
    490490             ELSE
    491491                to_be_resorted => s_av
  • palm/trunk/SOURCE/data_output_spectra.f90

    r1834 r1960  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Additional default spectra for passive scalar
    2222!
    2323! Former revisions:
     
    172172                pr = 5
    173173
     174             CASE ( 's' )
     175                pr = 6
     176
    174177             CASE DEFAULT
    175178!
  • palm/trunk/SOURCE/flow_statistics.f90

    r1919 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Separate humidity and passive scalar
    2222!
    2323! Former revisions:
     
    197197    USE arrays_3d,                                                             &
    198198        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, ol, p, prho, prr, pt, q, qc, ql,&
    199                qr, qs, qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt,  &
    200                td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us,&
    201                usws, uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
     199               qr, qs, qsws, qswst, rho, s, sa, ss, ssws, sswst, saswsb,       &
     200               saswst, shf, td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q,        &
     201               time_vert, ts, tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, &
     202               vswst, w, w_subs, zw
    202203       
    203204    USE cloud_parameters,                                                      &
     
    360361
    361362          DO  i = 0, threads_per_task-1
    362              sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
    363              IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
    364              IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
    365                                                    sums_wsqs_ws_l(:,i) !w*q*
     363             sums_l(:,17,i)                        = sums_wspts_ws_l(:,i) ! w*pt*
     364             IF ( ocean          ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
     365             IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)  ! w*q*
     366             IF ( passive_scalar ) sums_l(:,116,i) = sums_wsss_ws_l(:,i)  ! w*s*
    366367          ENDDO
    367368
     
    440441             DO  j =  nys, nyn
    441442                DO  k = nzb_s_inner(j,i), nzt+1
    442                    sums_l(k,41,tn)  = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
     443                   sums_l(k,117,tn)  = sums_l(k,117,tn) + s(k,j,i) * rmask(j,i,sr)
    443444                ENDDO
    444445             ENDDO
     
    465466             ENDIF
    466467             IF ( passive_scalar )  THEN
    467                 sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
     468                sums_l(:,117,0) = sums_l(:,117,0) + sums_l(:,117,i)
    468469             ENDIF
    469470          ENDDO
     
    506507       IF ( passive_scalar )  THEN
    507508          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    508           CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
     509          CALL MPI_ALLREDUCE( sums_l(nzb,117,0), sums(nzb,117), nzt+2-nzb,       &
    509510                              MPI_REAL, MPI_SUM, comm2d, ierr )
    510511       ENDIF
     
    522523          ENDIF
    523524       ENDIF
    524        IF ( passive_scalar )  sums(:,41) = sums_l(:,41,0)
     525       IF ( passive_scalar )  sums(:,117) = sums_l(:,117,0)
    525526#endif
    526527
     
    560561!
    561562!--    Passive scalar
    562        IF ( passive_scalar )  hom(:,1,41,sr) = sums(:,41) /                    &
    563             ngp_2dh_s_inner(:,sr)                    ! s (q)
     563       IF ( passive_scalar )  hom(:,1,117,sr) = sums(:,117) /                  &
     564            ngp_2dh_s_inner(:,sr)                    ! s
    564565
    565566!
     
    598599                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
    599600                ENDIF
    600 
     601                IF ( passive_scalar )  THEN
     602                   sums_l(k,118,tn) = sums_l(k,118,tn) + &
     603                                  ( s(k,j,i)-hom(k,1,117,sr) )**2 * rmask(j,i,sr)
     604                ENDIF
    601605!
    602606!--             Higher moments
     
    627631                sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
    628632                                            qs(j,i)   * rmask(j,i,sr)
     633             ENDIF
     634             IF ( passive_scalar )  THEN
     635                sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
     636                                            ss(j,i)   * rmask(j,i,sr)
    629637             ENDIF
    630638          ENDDO
     
    738746!--             Passive scalar flux
    739747                IF ( passive_scalar )  THEN
    740                    sums_l(k,48,tn) = sums_l(k,48,tn)                           &
     748                   sums_l(k,119,tn) = sums_l(k,119,tn)                         &
    741749                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    742                                                * ( q(k+1,j,i) - q(k,j,i) )     &
     750                                               * ( s(k+1,j,i) - s(k,j,i) )     &
    743751                                               * ddzu(k+1) * rmask(j,i,sr)
    744752                ENDIF
     
    784792                ENDIF
    785793                IF ( passive_scalar )  THEN
    786                    sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
    787                                        qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
     794                   sums_l(nzb,119,tn) = sums_l(nzb,119,tn) +                     &
     795                                       ssws(j,i) * rmask(j,i,sr) ! w"s"
    788796                ENDIF
    789797             ENDIF
     
    863871                ENDIF
    864872                IF ( passive_scalar )  THEN
    865                    sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
    866                                        qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
     873                   sums_l(nzt,119,tn) = sums_l(nzt,119,tn) + &
     874                                       sswst(j,i) * rmask(j,i,sr) ! w"s"
    867875                ENDIF
    868876             ENDIF
     
    953961                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
    954962                     .OR. sr /= 0 ) )  THEN
    955                    pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +              &
    956                                     q(k+1,j,i) - hom(k+1,1,41,sr) )
    957                    sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *        &
     963                   pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +             &
     964                                    s(k+1,j,i) - hom(k+1,1,117,sr) )
     965                   sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *      &
    958966                                                       rmask(j,i,sr)
    959967                ENDIF
     
    10111019                                      q(k+1,j,i) - hom(k+1,1,41,sr) )
    10121020                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *      &
     1021                                       rmask(j,i,sr)
     1022                  ENDIF
     1023                  IF ( passive_scalar )  THEN
     1024                     pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +            &
     1025                                      s(k+1,j,i) - hom(k+1,1,117,sr) )
     1026                     sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *      &
    10131027                                       rmask(j,i,sr)
    10141028                  ENDIF
     
    12791293          sums(k,81:88)         = sums(k,81:88)         / ngp_2dh(sr)
    12801294          sums(k,89:114)        = sums(k,89:114)        / ngp_2dh(sr)
     1295          sums(k,116)           = sums(k,116)           / ngp_2dh(sr)
     1296          sums(k,119)           = sums(k,119)           / ngp_2dh(sr)
    12811297          IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
    12821298             sums(k,8:11)          = sums(k,8:11)          / ngp_2dh_s_inner(k,sr)
     
    12871303             sums(k,64)            = sums(k,64)            / ngp_2dh_s_inner(k,sr)
    12881304             sums(k,70:80)         = sums(k,70:80)         / ngp_2dh_s_inner(k,sr)
    1289              sums(k,115:pr_palm-2) = sums(k,115:pr_palm-2) / ngp_2dh_s_inner(k,sr)
     1305             sums(k,118)           = sums(k,118)           / ngp_2dh_s_inner(k,sr)
     1306             sums(k,120:pr_palm-2) = sums(k,120:pr_palm-2) / ngp_2dh_s_inner(k,sr)
    12901307          ENDIF
    12911308       ENDDO
     
    12981315                                    ngp_2dh(sr)
    12991316       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
     1317                                    ngp_2dh(sr)
     1318       sums(nzb+13,pr_palm)       = sums(nzb+13,pr_palm)       / &    ! ss
    13001319                                    ngp_2dh(sr)
    13011320!--    eges, e*
     
    14431462       hom(:,1,114,sr) = sums(:,114)            !: L
    14441463
     1464       IF ( passive_scalar )  THEN
     1465          hom(:,1,119,sr) = sums(:,119)     ! w"s"
     1466          hom(:,1,116,sr) = sums(:,116)     ! w*s*
     1467          hom(:,1,120,sr) = sums(:,119) + sums(:,116)    ! ws
     1468       ENDIF
     1469
    14451470       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
    14461471                                       ! u*, w'u', w'v', t* (in last profile)
     
    15891614       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
    15901615
     1616       IF ( passive_scalar )  THEN
     1617          ts_value(24,sr) = hom(nzb+13,1,119,sr)       ! w"s" ( to do ! )
     1618          ts_value(25,sr) = hom(nzb+13,1,pr_palm,sr)   ! s*
     1619       ENDIF
     1620
    15911621!
    15921622!--    Collect land surface model timeseries
     
    16601690    USE arrays_3d,                                                             &
    16611691        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
    1662                qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q,&
    1663                td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
    1664                uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
     1692               qsws, qswst, rho, s, sa, saswsb, saswst, shf, ss, ssws, sswst,  &
     1693               td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts,      &
     1694               tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, vswst, w,      &
     1695               w_subs, zw
    16651696                 
    16661697       
     
    18331864             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
    18341865             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
    1835              IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
    1836                                                    sums_wsqs_ws_l(:,i) !w*q*
     1866             IF ( humidity       )  sums_l(:,49,i)  = sums_wsqs_ws_l(:,i) !w*q*
     1867             IF ( passive_scalar )  sums_l(:,116,i) = sums_wsss_ws_l(:,i) !w*s*
    18371868          ENDDO
    18381869
     
    19381969       IF ( passive_scalar )  THEN
    19391970          !$OMP DO
    1940           !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l ) create( s1 )
     1971          !$acc parallel loop gang present( s, rflags_invers, rmask, sums_l ) create( s1 )
    19411972          DO  k = nzb, nzt+1
    19421973             s1 = 0
     
    19441975             DO  i = nxl, nxr
    19451976                DO  j =  nys, nyn
    1946                    s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    1947                 ENDDO
    1948              ENDDO
    1949              sums_l(k,41,tn) = s1
     1977                   s1 = s1 + s(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     1978                ENDDO
     1979             ENDDO
     1980             sums_l(k,117,tn) = s1
    19501981          ENDDO
    19511982          !$acc end parallel loop
     
    19812012             IF ( passive_scalar )  THEN
    19822013                !$acc parallel present( sums_l )
    1983                 sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
     2014                sums_l(:,117,0) = sums_l(:,117,0) + sums_l(:,117,i)
    19842015                !$acc end parallel
    19852016             ENDIF
     
    20242055       IF ( passive_scalar )  THEN
    20252056          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2026           CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
     2057          CALL MPI_ALLREDUCE( sums_l(nzb,117,0), sums(nzb,117), nzt+2-nzb,     &
    20272058                              MPI_REAL, MPI_SUM, comm2d, ierr )
    20282059       ENDIF
     
    20532084       IF ( passive_scalar )  THEN
    20542085          !$acc parallel present( sums, sums_l )
    2055           sums(:,41) = sums_l(:,41,0)
     2086          sums(:,117) = sums_l(:,117,0)
    20562087          !$acc end parallel
    20572088       ENDIF
     
    21022133       IF ( passive_scalar )  THEN
    21032134          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
    2104           sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
    2105           hom(:,1,41,sr) = sums(:,41)                ! s (q)
     2135          sums(:,117)     = sums(:,117) / ngp_2dh_s_inner(:,sr)
     2136          hom(:,1,117,sr) = sums(:,117)                ! s
    21062137          !$acc end parallel
    21072138       ENDIF
     
    22302261          ENDDO
    22312262          sums_l(nzb+12,pr_palm,tn) = s1
     2263          !$acc end parallel
     2264       ENDIF
     2265
     2266       IF ( passive_scalar )  THEN
     2267          !$acc parallel present( ss, rmask, sums_l ) create( s1 )
     2268          s1 = 0
     2269          !$acc loop vector collapse( 2 ) reduction( +: s1 )
     2270          DO  i = nxl, nxr
     2271             DO  j =  nys, nyn
     2272                s1 = s1 + ss(j,i) * rmask(j,i,sr)
     2273             ENDDO
     2274          ENDDO
     2275          sums_l(nzb+13,pr_palm,tn) = s1
    22322276          !$acc end parallel
    22332277       ENDIF
     
    24112455       IF ( passive_scalar )  THEN
    24122456
    2413           !$acc parallel loop gang present( ddzu, kh, q, rflags_invers, rmask, sums_l ) create( s1 )
     2457          !$acc parallel loop gang present( ddzu, kh, s, rflags_invers, rmask, sums_l ) create( s1 )
    24142458          DO  k = nzb, nzt_diff
    24152459             s1 = 0
     
    24182462                DO  j = nys, nyn
    24192463                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
    2420                                     * ( q(k+1,j,i) - q(k,j,i) )                &
     2464                                    * ( s(k+1,j,i) - s(k,j,i) )                &
    24212465                                    * ddzu(k+1) * rmask(j,i,sr)                &
    24222466                                    * rflags_invers(j,i,k+1)
    24232467                ENDDO
    24242468             ENDDO
    2425              sums_l(k,48,tn) = s1
     2469             sums_l(k,119,tn) = s1
    24262470          ENDDO
    24272471          !$acc end parallel loop
     
    25322576
    25332577             !$OMP DO
    2534              !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
     2578             !$acc parallel present( ssws, rmask, sums_l ) create( s1 )
    25352579             s1 = 0
    25362580             !$acc loop vector collapse( 2 ) reduction( +: s1 )
    25372581             DO  i = nxl, nxr
    25382582                DO  j =  nys, nyn
    2539                    s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
    2540                 ENDDO
    2541              ENDDO
    2542              sums_l(nzb,48,tn) = s1
     2583                   s1 = s1 + ssws(j,i) * rmask(j,i,sr)  ! w"s"
     2584                ENDDO
     2585             ENDDO
     2586             sums_l(nzb,119,tn) = s1
    25432587             !$acc end parallel
    25442588
     
    26512695
    26522696             !$OMP DO
    2653              !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
     2697             !$acc parallel present( sswst, rmask, sums_l ) create( s1 )
    26542698             s1 = 0
    26552699             !$acc loop vector collapse( 2 ) reduction( +: s1 )
    26562700             DO  i = nxl, nxr
    26572701                DO  j =  nys, nyn
    2658                    s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
    2659                 ENDDO
    2660              ENDDO
    2661              sums_l(nzt,48,tn) = s1
     2702                   s1 = s1 + sswst(j,i) * rmask(j,i,sr) ! w"s"
     2703                ENDDO
     2704             ENDDO
     2705             sums_l(nzt,119,tn) = s1
    26622706             !$acc end parallel
    26632707
     
    28902934       IF ( passive_scalar  .AND.  ( .NOT. ws_scheme_sca  .OR.  sr /= 0 ) )  THEN
    28912935
    2892           !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
     2936          !$acc parallel loop gang present( hom, s, rflags_invers, rmask, sums_l, w ) create( s1 )
    28932937          DO  k = nzb, nzt_diff
    28942938             s1 = 0
     
    28962940             DO  i = nxl, nxr
    28972941                DO  j = nys, nyn
    2898                    s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +          &
    2899                                         q(k+1,j,i) - hom(k+1,1,41,sr) )        &
     2942                   s1 = s1 + 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +          &
     2943                                        s(k+1,j,i) - hom(k+1,1,117,sr) )        &
    29002944                                    * w(k,j,i) * rmask(j,i,sr)                 &
    29012945                                    * rflags_invers(j,i,k+1)
     
    29813025                ENDDO
    29823026                sums_l(k,49,tn) = s1
     3027             ENDDO
     3028             !$acc end parallel loop
     3029
     3030          ENDIF
     3031
     3032          IF ( passive_scalar )  THEN
     3033
     3034             !$acc parallel loop gang present( hom, s, rflags_invers, rmask, sums_l, w ) create( s1 )
     3035             DO  k = nzb, nzt_diff
     3036                s1 = 0
     3037                !$acc loop vector collapse( 2 ) reduction( +: s1 )
     3038                DO  i = nxl, nxr
     3039                   DO  j = nys, nyn
     3040                      s1 = s1 + 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +      &
     3041                                           s(k+1,j,i) - hom(k+1,1,117,sr) )    &
     3042                                       * w(k,j,i) * rmask(j,i,sr)              &
     3043                                       * rflags_invers(j,i,k+1)
     3044                   ENDDO
     3045                ENDDO
     3046                sums_l(k,116,tn) = s1
    29833047             ENDDO
    29843048             !$acc end parallel loop
  • palm/trunk/SOURCE/header.f90

    r1958 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Treat humidity and passive scalar separately.
     22! Modify misleading information concerning humidity.
     23! Bugfix, change unit for humidity flux.
    2224!
    2325! Former revisions:
     
    266268
    267269    USE arrays_3d,                                                             &
    268         ONLY:  pt_init, qsws, q_init, sa_init, shf, ug, vg, w_subs, zu, zw
     270        ONLY:  pt_init, qsws, q_init, s_init, sa_init, shf, ug, vg, w_subs, zu,&
     271               zw
    269272       
    270273    USE control_parameters
     
    10221025
    10231026    IF ( passive_scalar )  THEN
    1024        IF ( ibc_q_b == 0 )  THEN
     1027       IF ( ibc_s_b == 0 )  THEN
    10251028          r_lower = 's(0)     = s_surface'
    10261029       ELSE
    10271030          r_lower = 's(0)     = s(1)'
    10281031       ENDIF
    1029        IF ( ibc_q_t == 0 )  THEN
     1032       IF ( ibc_s_t == 0 )  THEN
    10301033          r_upper =  's(nzt)   = s_top'
    10311034       ELSE
     
    10521055          ENDIF
    10531056       ENDIF
    1054        IF ( passive_scalar  .AND.  constant_waterflux )  THEN
    1055           WRITE ( io, 313 ) surface_waterflux
     1057       IF ( passive_scalar  .AND.  constant_scalarflux )  THEN
     1058          WRITE ( io, 313 ) surface_scalarflux
    10561059       ENDIF
    10571060    ENDIF
     
    10701073          WRITE ( io, 309 )  top_salinityflux
    10711074       ENDIF
    1072        IF ( humidity  .OR.  passive_scalar )  THEN
    1073           WRITE ( io, 315 )
    1074        ENDIF
     1075       IF ( humidity       )  WRITE ( io, 315 )
     1076       IF ( passive_scalar )  WRITE ( io, 315 )
    10751077    ENDIF
    10761078
     
    10831085          WRITE ( io, 312 )
    10841086       ENDIF
    1085        IF ( passive_scalar  .AND.  .NOT. constant_waterflux )  THEN
     1087       IF ( passive_scalar  .AND.  .NOT. constant_scalarflux )  THEN
    10861088          WRITE ( io, 314 )
    10871089       ENDIF
     
    11531155!-- Initial humidity profile
    11541156!-- Building output strings, starting with surface humidity
    1155     IF ( humidity  .OR.  passive_scalar )  THEN
     1157    IF ( humidity )  THEN
    11561158       WRITE ( temperatures, '(E8.1)' )  q_surface
    11571159       gradients = '--------'
     
    11811183       ENDDO
    11821184
    1183        IF ( humidity )  THEN
    1184           IF ( .NOT. nudging )  THEN
    1185              WRITE ( io, 421 )  TRIM( coordinates ), TRIM( temperatures ), &
    1186                                 TRIM( gradients ), TRIM( slices )
    1187           ENDIF
    1188        ELSE
    1189           WRITE ( io, 422 )  TRIM( coordinates ), TRIM( temperatures ), &
     1185       IF ( .NOT. nudging )  THEN
     1186          WRITE ( io, 421 )  TRIM( coordinates ), TRIM( temperatures ),        &
    11901187                             TRIM( gradients ), TRIM( slices )
    11911188       ENDIF
    11921189    ENDIF
     1190!
     1191!-- Initial scalar profile
     1192!-- Building output strings, starting with surface humidity
     1193    IF ( passive_scalar )  THEN
     1194       WRITE ( temperatures, '(E8.1)' )  s_surface
     1195       gradients = '--------'
     1196       slices = '       0'
     1197       coordinates = '     0.0'
     1198       i = 1
     1199       DO  WHILE ( s_vertical_gradient_level_ind(i) /= -9999 )
     1200         
     1201          WRITE (coor_chr,'(E8.1,4X)')  s_init(q_vertical_gradient_level_ind(i))
     1202          temperatures = TRIM( temperatures ) // '  ' // TRIM( coor_chr )
     1203
     1204          WRITE (coor_chr,'(E8.1,4X)')  s_vertical_gradient(i)
     1205          gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
     1206         
     1207          WRITE (coor_chr,'(I8,4X)')  s_vertical_gradient_level_ind(i)
     1208          slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
     1209         
     1210          WRITE (coor_chr,'(F8.1,4X)')  s_vertical_gradient_level(i)
     1211          coordinates = TRIM( coordinates ) // '  '  // TRIM( coor_chr )
     1212
     1213          IF ( i == 10 )  THEN
     1214             EXIT
     1215          ELSE
     1216             i = i + 1
     1217          ENDIF
     1218
     1219       ENDDO
     1220
     1221       WRITE ( io, 422 )  TRIM( coordinates ), TRIM( temperatures ),           &
     1222                          TRIM( gradients ), TRIM( slices )
     1223    ENDIF   
    11931224
    11941225!
     
    19952026310 FORMAT (//'    1D-Model:'// &
    19962027             '       Rif value range:   ',F6.2,' <= rif <=',F6.2)
    1997 311 FORMAT ('       Predefined constant humidity flux: ',E10.3,' m/s')
     2028311 FORMAT ('       Predefined constant humidity flux: ',E10.3,' kg/kg m/s')
    19982029312 FORMAT ('       Predefined surface humidity')
    19992030313 FORMAT ('       Predefined constant scalar flux: ',E10.3,' kg/(m**2 s)')
     
    21582189430 FORMAT (//' Cloud physics quantities / methods:'/ &
    21592190              ' ----------------------------------'/)
    2160 431 FORMAT ('    Humidity is treated as purely passive scalar (no condensati', &
    2161                  'on)')
     2191431 FORMAT ('    Humidity is considered, bu no condensation')
    21622192432 FORMAT ('    Bulk scheme with liquid water potential temperature and'/ &
    21632193            '    total water content is used.'/ &
  • palm/trunk/SOURCE/inflow_turbulence.f90

    r1818 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Separate humidity and passive scalar
    2222!
    2323! Former revisions:
     
    7171
    7272    USE arrays_3d,                                                             &
    73         ONLY:  e, inflow_damping_factor, mean_inflow_profiles, pt, q, u, v, w
     73        ONLY:  e, inflow_damping_factor, mean_inflow_profiles, pt, q, s, u, v, w
    7474       
    7575    USE control_parameters,                                                    &
     
    9898    INTEGER(iwp) ::  prev     !< ID of sending PE for y-shift
    9999
    100     REAL(wp), DIMENSION(nzb:nzt+1,6,nbgp)           ::                         &
     100    REAL(wp), DIMENSION(nzb:nzt+1,7,nbgp)           ::                         &
    101101       avpr               !< stores averaged profiles at recycling plane
    102     REAL(wp), DIMENSION(nzb:nzt+1,6,nbgp)           ::                         &
     102    REAL(wp), DIMENSION(nzb:nzt+1,7,nbgp)           ::                         &
    103103       avpr_l             !< auxiliary variable to calculate avpr
    104     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,6,nbgp) ::                         &
     104    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,7,nbgp) ::                         &
    105105       inflow_dist        !< turbulence signal of vars, added at inflow boundary
    106     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,6,nbgp) ::                         &
     106    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,7,nbgp) ::                         &
    107107       local_inflow_dist  !< auxiliary variable for inflow_dist, used for yshift
    108108
     
    112112!-- Carry out spanwise averaging in the recycling plane
    113113    avpr_l = 0.0_wp
    114     ngp_pr = ( nzt - nzb + 2 ) * 6 * nbgp
     114    ngp_pr = ( nzt - nzb + 2 ) * 7 * nbgp
    115115    ngp_ifd = ngp_pr * ( nyn - nys + 1 + 2 * nbgp )
    116116
     
    131131                avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)
    132132                avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i)
    133                 IF ( humidity  .OR.  passive_scalar )                          &
     133                IF ( humidity )                                                &
    134134                   avpr_l(k,6,l) = avpr_l(k,6,l) + q(k,j,i)
     135                IF ( passive_scalar )                                          &
     136                   avpr_l(k,7,l) = avpr_l(k,7,l) + s(k,j,i)
    135137
    136138             ENDDO
     
    156158             avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)
    157159             avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i)
    158              IF ( humidity  .OR.  passive_scalar )                             &
     160             IF ( humidity )                                                   &
    159161                avpr_l(k,6,l) = avpr_l(k,6,l) + q(k,j,i)
     162             IF ( passive_scalar )                                             &
     163                avpr_l(k,7,l) = avpr_l(k,7,l) + s(k,j,i)
    160164
    161165          ENDDO
     
    183187                inflow_dist(k,j,4,l) = pt(k,j,i)  - avpr(k,4,l)
    184188                inflow_dist(k,j,5,l) = e(k,j,i)   - avpr(k,5,l)
    185                 IF ( humidity  .OR.  passive_scalar )                          &
     189                IF ( humidity )                                                &
    186190                   inflow_dist(k,j,6,l) = q(k,j,i) - avpr(k,6,l)
     191                IF ( passive_scalar )                                          &
     192                   inflow_dist(k,j,7,l) = s(k,j,i) - avpr(k,7,l)
    187193            ENDDO
    188194          ENDDO
     
    201207             inflow_dist(k,j,4,l) = pt(k,j,i)  - avpr(k,4,l)
    202208             inflow_dist(k,j,5,l) = e(k,j,i)   - avpr(k,5,l)
    203              IF ( humidity  .OR.  passive_scalar )                             &
     209             IF ( humidity )                                                   &
    204210                inflow_dist(k,j,6,l) = q(k,j,i) - avpr(k,6,l)
     211             IF ( passive_scalar )                                             &
     212                inflow_dist(k,j,7,l) = s(k,j,i) - avpr(k,7,l)
    205213             
    206214          ENDDO
     
    276284             e(k,j,-nbgp:-1)  = MAX( e(k,j,-nbgp:-1), 0.0_wp )
    277285
    278              IF ( humidity  .OR.  passive_scalar )                          &
     286             IF ( humidity )                                                &
    279287                q(k,j,-nbgp:-1)  = mean_inflow_profiles(k,6) +              &
    280288                        inflow_dist(k,j,6,1:nbgp) * inflow_damping_factor(k)
     289             IF ( passive_scalar )                                          &
     290                s(k,j,-nbgp:-1)  = mean_inflow_profiles(k,7) +              &
     291                        inflow_dist(k,j,7,1:nbgp) * inflow_damping_factor(k)
    281292
    282293          ENDDO
  • palm/trunk/SOURCE/init_1d_model.f90

    r1818 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Remove passive_scalar from IF-statements, as 1D-scalar profile is effectively
     22! not used.
     23! Formatting adjustment
    2224!
    2325! Former revisions:
     
    108110    USE control_parameters,                                                    &
    109111        ONLY:  constant_diffusion, constant_flux_layer, f, humidity, kappa,    &
    110                km_constant, mixing_length_1d, passive_scalar, prandtl_number,  &
     112               km_constant, mixing_length_1d, prandtl_number,                  &
    111113               roughness_length, simulated_time_chr, z0h_factor
    112114
     
    194196    z01d  = roughness_length
    195197    z0h1d = z0h_factor * z01d
    196     IF ( humidity .OR. passive_scalar )  qs1d = 0.0_wp
     198    IF ( humidity )  qs1d = 0.0_wp
    197199
    198200!
     
    232234               humidity, intermediate_timestep_count,                          &
    233235               intermediate_timestep_count_max, f, g, ibc_e_b, kappa,          & 
    234                mixing_length_1d, passive_scalar,                               &
     236               mixing_length_1d,                                               &
    235237               simulated_time_chr, timestep_scheme, tsc, zeta_max, zeta_min
    236238               
     
    619621                e1d(nzb) = e1d(nzb+1)
    620622
    621                 IF ( humidity .OR. passive_scalar ) THEN
     623                IF ( humidity ) THEN
    622624!
    623625!--                Compute q*
    624626                   IF ( rif1d(1) >= 0.0_wp )  THEN
    625627!
    626 !--                Stable stratification
    627                    qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /            &
     628!--                   Stable stratification
     629                      qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /         &
    628630                          ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * rif1d(nzb+1) * &
    629631                                          ( zu(nzb+1) - z0h1d ) / zu(nzb+1)    &
    630632                          )
    631                 ELSE
    632 !
    633 !--                Unstable stratification
    634                    a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) )
    635                    b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /                 &
    636                                       zu(nzb+1) * z0h1d )
    637 !
    638 !--                In the borderline case the formula for stable stratification
    639 !--                must be applied, because otherwise a zero division would
    640 !--                occur in the argument of the logarithm.
    641                    IF ( a == 1.0_wp  .OR.  b == 1.0_wp )  THEN
    642                       qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /         &
    643                              ( LOG( zu(nzb+1) / z0h1d ) +                      &
    644                                5.0_wp * rif1d(nzb+1) *                         &
    645                                ( zu(nzb+1) - z0h1d ) / zu(nzb+1)               &
    646                              )
    647633                   ELSE
    648                       qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /         &
    649                              LOG( (a-1.0_wp) / (a+1.0_wp) *                    &
    650                                   (b+1.0_wp) / (b-1.0_wp) )
    651                    ENDIF
    652                 ENDIF               
     634!
     635!--                   Unstable stratification
     636                      a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) )
     637                      b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /              &
     638                                         zu(nzb+1) * z0h1d )
     639!
     640!--                   In the borderline case the formula for stable stratification
     641!--                   must be applied, because otherwise a zero division would
     642!--                   occur in the argument of the logarithm.
     643                      IF ( a == 1.0_wp  .OR.  b == 1.0_wp )  THEN
     644                         qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /      &
     645                                ( LOG( zu(nzb+1) / z0h1d ) +                   &
     646                                  5.0_wp * rif1d(nzb+1) *                      &
     647                                  ( zu(nzb+1) - z0h1d ) / zu(nzb+1)            &
     648                                )
     649                      ELSE
     650                         qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /      &
     651                                LOG( (a-1.0_wp) / (a+1.0_wp) *                 &
     652                                     (b+1.0_wp) / (b-1.0_wp) )
     653                      ENDIF
     654                   ENDIF               
    653655                ELSE
    654656                   qs1d = 0.0_wp
  • palm/trunk/SOURCE/init_3d_model.f90

    r1958 r1960  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Separate humidity and passive scalar
     22! Increase dimension for mean_inflow_profiles
     23! Remove inadvertent write-statement
     24! Bugfix, large-scale forcing is still not implemented for passive scalars
    2225!
    2326! Former revisions:
     
    469472    ENDIF
    470473
    471     IF ( humidity  .OR.  passive_scalar )  THEN
    472 !
    473 !--    2D-humidity/scalar arrays
     474    IF ( humidity )  THEN
     475!
     476!--    2D-humidity
    474477       ALLOCATE ( qs(nysg:nyng,nxlg:nxrg),                                     &
    475478                  qsws(nysg:nyng,nxlg:nxrg),                                   &
     
    477480
    478481!
    479 !--    3D-humidity/scalar arrays
     482!--    3D-humidity
    480483#if defined( __nopointer )
    481484       ALLOCATE( q(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
     
    489492
    490493!
    491 !--    3D-arrays needed for humidity only
     494!--    3D-arrays needed for humidity
    492495       IF ( humidity )  THEN
    493496#if defined( __nopointer )
     
    570573       ENDIF
    571574
     575    ENDIF
     576   
     577   
     578    IF ( passive_scalar )  THEN
     579!
     580!--    2D-scalar arrays
     581       ALLOCATE ( ss(nysg:nyng,nxlg:nxrg),                                     &
     582                  ssws(nysg:nyng,nxlg:nxrg),                                   &
     583                  sswst(nysg:nyng,nxlg:nxrg) )
     584
     585!
     586!--    3D scalar arrays
     587#if defined( __nopointer )
     588       ALLOCATE( s(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
     589                 s_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
     590                 ts_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     591#else
     592       ALLOCATE( s_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
     593                 s_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
     594                 s_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     595#endif
    572596    ENDIF
    573597
     
    678702    w  => w_1;   w_p  => w_2;   tw_m  => w_3
    679703
    680     IF ( humidity  .OR.  passive_scalar )  THEN
     704    IF ( humidity )  THEN
    681705       q => q_1;  q_p => q_2;  tq_m => q_3
    682706       IF ( humidity )  THEN
     
    696720       ENDIF
    697721    ENDIF
     722   
     723    IF ( passive_scalar )  THEN
     724       s => s_1;  s_p => s_2;  ts_m => s_3
     725    ENDIF   
    698726
    699727    IF ( ocean )  THEN
     
    773801          ENDDO
    774802
    775           IF ( humidity  .OR.  passive_scalar )  THEN
     803          IF ( humidity )  THEN
    776804             DO  i = nxlg, nxrg
    777805                DO  j = nysg, nyng
     
    788816
    789817             ENDIF
     818          ENDIF
     819          IF ( passive_scalar )  THEN
     820             DO  i = nxlg, nxrg
     821                DO  j = nysg, nyng
     822                   s(:,j,i) = s_init
     823                ENDDO
     824             ENDDO   
    790825          ENDIF
    791826
     
    830865!--       This could actually be computed more accurately in the 1D model.
    831866!--       Update when opportunity arises!
    832           IF ( humidity  .OR.  passive_scalar )  THEN
     867          IF ( humidity )  THEN
    833868             qs = 0.0_wp
    834869             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     
    837872             ENDIF
    838873          ENDIF
     874!
     875!--       Initialize scaling parameter for passive scalar
     876          IF ( passive_scalar ) ss = 0.0_wp         
    839877
    840878!
     
    881919             u_init  = unudge(:,1)
    882920             v_init  = vnudge(:,1)
    883              IF ( humidity  .OR.  passive_scalar )  THEN
     921             IF ( humidity  )  THEN ! is passive_scalar correct???
    884922                q_init = qnudge(:,1)
    885923             ENDIF
     
    915953          ENDIF
    916954
    917           IF ( humidity  .OR.  passive_scalar )  THEN
     955          IF ( humidity )  THEN
    918956             DO  i = nxlg, nxrg
    919957                DO  j = nysg, nyng
     
    931969
    932970             ENDIF
     971          ENDIF
     972         
     973          IF ( passive_scalar )  THEN
     974             DO  i = nxlg, nxrg
     975                DO  j = nysg, nyng
     976                   s(:,j,i) = s_init
     977                ENDDO
     978             ENDDO
    933979          ENDIF
    934980
     
    9751021          vsws  = 0.0_wp
    9761022          vswst = top_momentumflux_v
    977           IF ( humidity  .OR.  passive_scalar ) qs = 0.0_wp
     1023          IF ( humidity )       qs = 0.0_wp
     1024          IF ( passive_scalar ) ss = 0.0_wp
    9781025
    9791026!
     
    10151062!
    10161063!--    Calculate virtual potential temperature
    1017        IF ( humidity ) vpt = pt * ( 1.0_wp + 0.61_wp * q )
     1064       IF ( humidity )  vpt = pt * ( 1.0_wp + 0.61_wp * q )
    10181065
    10191066!
     
    10531100!
    10541101!--       Store initial scalar profile
    1055           hom(:,1,26,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
     1102          hom(:,1,115,:) = SPREAD(  s(:,nys,nxl), 2, statistic_regions+1 )
    10561103       ENDIF
    10571104
     
    11161163!
    11171164!--       Determine the near-surface water flux
    1118           IF ( humidity  .OR.  passive_scalar )  THEN
     1165          IF ( humidity )  THEN
    11191166             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    11201167                qrsws = 0.0_wp
     
    11381185             ENDIF
    11391186          ENDIF
     1187!
     1188!--       Initialize the near-surface scalar flux
     1189          IF ( passive_scalar )  THEN
     1190             IF ( constant_scalarflux )  THEN
     1191                ssws   = surface_scalarflux
     1192!
     1193!--             Over topography surface_scalarflux is replaced by
     1194!--             wall_scalarflux(0)
     1195                IF ( TRIM( topography ) /= 'flat' )  THEN
     1196                   wall_sflux = wall_scalarflux
     1197                   DO  i = nxlg, nxrg
     1198                      DO  j = nysg, nyng
     1199                         IF ( nzb_s_inner(j,i) /= 0 )  ssws(j,i) = wall_sflux(0)
     1200                      ENDDO
     1201                   ENDDO
     1202                ENDIF
     1203             ENDIF
     1204          ENDIF         
    11401205
    11411206       ENDIF
     
    11521217             tswst = top_heatflux
    11531218
    1154              IF ( humidity  .OR.  passive_scalar )  THEN
     1219             IF ( humidity )  THEN
    11551220                qswst = 0.0_wp
    11561221                IF ( cloud_physics  .AND.  microphysics_seifert ) THEN
     
    11591224                ENDIF
    11601225             ENDIF
     1226             IF ( passive_scalar )  sswst = 0.0_wp
    11611227
    11621228             IF ( ocean )  THEN
     
    11921258          ENDIF
    11931259
    1194           IF ( humidity  .OR.  passive_scalar )  THEN
     1260          IF ( humidity  )  THEN
    11951261             IF (  .NOT.  constant_waterflux )  qsws   = 0.0_wp
    11961262             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     
    11991265             ENDIF
    12001266          ENDIF
     1267          IF ( passive_scalar  .AND.  .NOT.  constant_scalarflux )  ssws = 0.0_wp
    12011268
    12021269       ENDIF
     
    12611328!--    If required, change the surface humidity/scalar at the start of the 3D
    12621329!--    run
    1263        IF ( ( humidity .OR. passive_scalar ) .AND.                             &
    1264             q_surface_initial_change /= 0.0_wp )  THEN
     1330       IF ( humidity  .AND.  q_surface_initial_change /= 0.0_wp )              &
    12651331          q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
    1266        ENDIF
     1332         
     1333       IF ( passive_scalar .AND.  s_surface_initial_change /= 0.0_wp )         &
     1334          s(nzb,:,:) = s(nzb,:,:) + s_surface_initial_change
     1335       
    12671336
    12681337!
     
    12711340       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
    12721341
    1273        IF ( humidity  .OR.  passive_scalar )  THEN
     1342       IF ( humidity  )  THEN
    12741343          tq_m = 0.0_wp
    12751344          q_p = q
     
    12811350          ENDIF
    12821351       ENDIF
     1352       
     1353       IF ( passive_scalar )  THEN
     1354          ts_m = 0.0_wp
     1355          s_p  = s
     1356       ENDIF       
    12831357
    12841358       IF ( ocean )  THEN
     
    13231397#endif
    13241398       ENDDO
    1325        write(9,*) "EOF read binary"
    1326        flush(9)
    13271399
    13281400!
     
    13351407!--       profiles from the prerun. Alternatively, prescribed profiles
    13361408!--       for u,v-components can be used.
    1337           ALLOCATE( mean_inflow_profiles(nzb:nzt+1,6) )
     1409          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,7) )
    13381410
    13391411          IF ( use_prescribed_profile_data )  THEN
     
    13461418          mean_inflow_profiles(:,4) = hom_sum(:,4,0)       ! pt
    13471419          mean_inflow_profiles(:,5) = hom_sum(:,8,0)       ! e
    1348           mean_inflow_profiles(:,6) = hom_sum(:,41,0)      ! q
     1420          IF ( humidity )                                                      &
     1421             mean_inflow_profiles(:,6) = hom_sum(:,41,0)   ! q
     1422          IF ( passive_scalar )                                                &
     1423             mean_inflow_profiles(:,7) = hom_sum(:,115,0)   ! s
    13491424
    13501425!
     
    13731448                   pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
    13741449                   e(k,j,nxlg:-1)  = mean_inflow_profiles(k,5)
    1375                    IF ( humidity  .OR.  passive_scalar )                       &
     1450                   IF ( humidity )                                             &
    13761451                      q(k,j,nxlg:-1)  = mean_inflow_profiles(k,6)
     1452                   IF ( passive_scalar )                                       &
     1453                      s(k,j,nxlg:-1)  = mean_inflow_profiles(k,7)                     
    13771454                ENDDO
    13781455             ENDDO
     
    14591536!--    including ghost points)
    14601537       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
    1461        IF ( humidity  .OR.  passive_scalar )  THEN
     1538       IF ( humidity )  THEN
    14621539          q_p = q
    14631540          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     
    14661543          ENDIF
    14671544       ENDIF
    1468        IF ( ocean )  sa_p = sa
     1545       IF ( passive_scalar )  s_p  = s
     1546       IF ( ocean          )  sa_p = sa
    14691547
    14701548!
     
    14731551!--    there before they are set.
    14741552       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
    1475        IF ( humidity  .OR.  passive_scalar )  THEN
     1553       IF ( humidity )  THEN
    14761554          tq_m = 0.0_wp
    14771555          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     
    14801558          ENDIF
    14811559       ENDIF
    1482        IF ( ocean )  tsa_m = 0.0_wp
     1560       IF ( passive_scalar )  ts_m  = 0.0_wp
     1561       IF ( ocean          )  tsa_m = 0.0_wp
    14831562
    14841563       CALL location_message( 'finished', .TRUE. )
  • palm/trunk/SOURCE/modules.f90

    r1958 r1960  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Separate humidity and passive scalar
     22! +bc_s_t_val, diss_s_s, diss_l_s, flux_s_s, flux_l_s, s, sp, s1, s2, s3, ssws_av,
     23!  s_init, s_surf, sums_wsss_ws_l, ss, ssws, sswst, ts_m, wall_sflux
     24! +constant_scalarflux, ibc_s_b, ibc_s_t, s_vertical_gradient_level_ind
     25!
     26! Unused variables removed
     27! -gamma_x, gamma_y, gamma_z, var_x, var_y, var_z
     28!
     29! Change default values (in order to unify gradient calculation):
     30! pt_vertical_gradient_level, sa_vertical_gradient_level
    2231!
    2332! Former revisions:
     
    411420          dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor, l_grid,           &
    412421          ptdf_x, ptdf_y, p_surf, pt_surf, pt_init, qsws_surf, q_init, q_surf, &
    413           rdf, rdf_sc, ref_state, sa_init, shf_surf, timenudge, time_surf,     &
    414           time_vert, tmp_tnudge, ug, u_init, u_nzb_p1_for_vfc, vg, v_init,     &
    415           v_nzb_p1_for_vfc, w_subs, zu, zw
     422          rdf, rdf_sc, ref_state, s_init, s_surf, sa_init, shf_surf,           &
     423          timenudge, time_surf, time_vert, tmp_tnudge, ug, u_init,             &
     424          u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, w_subs, zu, zw
    416425
    417426    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::                                   &
    418427          c_u, c_v, c_w, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q,             &
    419           diss_s_qr, diss_s_sa, diss_s_u, diss_s_v, diss_s_w, dzu_mg, dzw_mg,  &
    420           flux_s_e, flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr, flux_s_sa,      &
    421           flux_s_u, flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg,                   &
    422           mean_inflow_profiles, nrs, nrsws, nrswst,                            &
     428          diss_s_qr, diss_s_s, diss_s_sa, diss_s_u, diss_s_v, diss_s_w, dzu_mg,&
     429          dzw_mg, flux_s_e, flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr,         &
     430          flux_s_s, flux_s_sa, flux_s_u, flux_s_v, flux_s_w, f1_mg, f2_mg,     &
     431          f3_mg, mean_inflow_profiles, nrs, nrsws, nrswst,                     &
    423432          ol,                                                                  & !< Obukhov length
    424433          precipitation_amount, precipitation_rate, ptnudge, pt_slope_ref,     &
    425434          qnudge, qs, qsws, qswst, qswst_remote, qrs, qrsws, qrswst,           &
    426           saswsb, saswst, shf, tnudge, td_lsa_lpt, td_lsa_q, td_sub_lpt,       &
     435          saswsb, saswst, shf, ss, ssws, sswst, tnudge, td_lsa_lpt, td_lsa_q,  &
     436          td_sub_lpt,                                                          &
    427437          td_sub_q, total_2d_a, total_2d_o, ts, tswst, ug_vert, unudge, us,    &
    428438          usws, uswst, vnudge, vg_vert, vsws, vswst, wnudge, wsubs_vert,       &
     
    433443    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
    434444          d, de_dx, de_dy, de_dz, diss, diss_l_e,                              &
    435           diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr, diss_l_sa, diss_l_u,      &
    436           diss_l_v, diss_l_w, flux_l_e, flux_l_nr, flux_l_pt, flux_l_q,        &
    437           flux_l_qr, flux_l_sa, flux_l_u, flux_l_v, flux_l_w, kh, km,          &
    438           l_wall, prr, p_loc, tend, tric,                                      &
     445          diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr, diss_l_s, diss_l_sa,      &
     446          diss_l_u, diss_l_v, diss_l_w, flux_l_e, flux_l_nr, flux_l_pt,        &
     447          flux_l_q, flux_l_qr, flux_l_s, flux_l_sa, flux_l_u, flux_l_v,        &
     448          flux_l_w, kh, km, l_wall, prr, p_loc, tend, tric,                    &
    439449          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,       &
    440450          w_m_n, w_m_r, w_m_s
     
    443453    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
    444454          e, e_p, nr, nr_p, p, prho, pt, pt_p, q, q_p, qc, ql, ql_c, ql_v,     &
    445           ql_vp, qr, qr_p, rho, sa, sa_p, te_m, tnr_m, tpt_m, tq_m, tqr_m,     &
    446           tsa_m, tu_m, tv_m, tw_m, u, u_p, v, v_p, vpt, w, w_p
     455          ql_vp, qr, qr_p, rho, s, s_p, sa, sa_p, te_m, tnr_m, tpt_m, tq_m,    &
     456          tqr_m, ts_m, tsa_m, tu_m, tv_m, tw_m, u, u_p, v, v_p, vpt, w, w_p
    447457#else
    448458    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
    449459          e_1, e_2, e_3, p, prho_1, nr_1, nr_2, nr_3, pt_1, pt_2, pt_3, q_1,   &
    450460          q_2, q_3, qc_1, ql_v, ql_vp, ql_1, ql_2, qr_1, qr_2, qr_3, rho_1,    &
    451           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
     461          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,&
     462          w_1, w_2, w_3
    452463
    453464    REAL(wp), DIMENSION(:,:,:), POINTER ::                                     &
    454465          e, e_p, nr, nr_p, prho, pt, pt_p, q, q_p, qc, ql, ql_c, qr, qr_p,    &
    455           rho, sa, sa_p, te_m, tnr_m, tpt_m, tq_m, tqr_m, tsa_m, tu_m, tv_m,   &
    456           tw_m, u, u_p, v, v_p, vpt, w, w_p
     466          rho, s, s_p, sa, sa_p, te_m, tnr_m, tpt_m, tq_m, tqr_m, ts_m,        &
     467          tsa_m, tu_m, tv_m, tw_m, u, u_p, v, v_p, vpt, w, w_p
    457468#endif
    458469
    459470    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  rif_wall, tri
    460 
    461     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: var_x, var_y, var_z, gamma_x,   &
    462                                            gamma_y, gamma_z
    463471
    464472
     
    481489                                              ol_av,                 & !< Avg. of Obukhov length
    482490                                              qsws_av,               & !< Avg. of surface moisture flux
     491                                              ssws_av,               & !< Avg. of surface scalar flux
    483492                                              shf_av,                & !< Avg. of surface heat flux
    484493                                              ts_av,                 & !< Avg. of characteristic temperature scale
     
    631640                     dz_stretch_level_index, ensemble_member_nr = 0, gamma_mg, gathered_size, &
    632641                     grid_level, ibc_e_b, ibc_p_b, ibc_p_t, &
    633                      ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, &
     642                     ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, ibc_s_t, &
    634643                     ibc_sa_t, ibc_uv_b, ibc_uv_t, &
    635644                     inflow_disturbance_begin = -1, inflow_disturbance_end = -1, &
     
    657666                     pt_vertical_gradient_level_ind(10) = -9999, &
    658667                     q_vertical_gradient_level_ind(10) = -9999, &
     668                     s_vertical_gradient_level_ind(10) = -9999, &                     
    659669                     sa_vertical_gradient_level_ind(10) = -9999, &
    660670                     section(100,3), section_xy(100) = -9999, &
     
    685695                constant_top_momentumflux = .FALSE., &
    686696                constant_top_salinityflux = .TRUE., &
     697                constant_scalarflux = .TRUE., &             
    687698                constant_waterflux = .TRUE., create_disturbances = .TRUE., &
    688699                data_output_2d_on_each_pe = .TRUE., &
     
    739750                 alpha_surface = 0.0_wp, atmos_ocean_sign = 1.0_wp, &
    740751                 averaging_interval = 0.0_wp, averaging_interval_pr = 9999999.9_wp, &
    741                  bc_pt_t_val, bc_q_t_val, bottom_salinityflux = 0.0_wp, &
     752                 bc_pt_t_val, bc_q_t_val, bc_s_t_val, bottom_salinityflux = 0.0_wp, &
    742753                 building_height = 50.0_wp, building_length_x = 50.0_wp, &
    743754                 building_length_y = 50.0_wp, building_wall_left = 9999999.9_wp, &
     
    813824                 mask_scale(3), &
    814825                 pt_vertical_gradient(10) = 0.0_wp, &
    815                  pt_vertical_gradient_level(10) = -9999999.9_wp, &
     826                 pt_vertical_gradient_level(10) = -1.0_wp, &
    816827                 q_vertical_gradient(10) = 0.0_wp, &
    817828                 q_vertical_gradient_level(10) = -1.0_wp, &
     
    819830                 s_vertical_gradient_level(10) = -1.0_wp, &
    820831                 sa_vertical_gradient(10) = 0.0_wp, &
    821                  sa_vertical_gradient_level(10) = -9999999.9_wp, &
     832                 sa_vertical_gradient_level(10) = -1.0_wp, &
    822833                 skip_time_domask(max_masks) = 9999999.9_wp, threshold(20) = 0.0_wp, &
    823834                 time_domask(max_masks) = 0.0_wp, &
     
    833844                 wall_humidityflux(0:4) = 0.0_wp, wall_nrflux(0:4) = 0.0_wp, &
    834845                 wall_qflux(0:4) = 0.0_wp, wall_qrflux(0:4) = 0.0_wp, &
    835                  wall_salinityflux(0:4) = 0.0_wp, wall_scalarflux(0:4) = 0.0_wp, &
     846                 wall_salinityflux(0:4) = 0.0_wp, wall_sflux(0:4) = 0.0_wp, wall_scalarflux(0:4) = 0.0_wp, &
    836847                 subs_vertical_gradient(10) = 0.0_wp, &
    837848                 subs_vertical_gradient_level(10) = -9999999.9_wp
     
    12551266                                                  sums_wsnrs_ws_l,                &
    12561267                                                  sums_wspts_ws_l,                &
     1268                                                  sums_wsqs_ws_l,                 &   
     1269                                                  sums_wsqrs_ws_l,                &
    12571270                                                  sums_wssas_ws_l,                &
    1258                                                   sums_wsqs_ws_l,                 &
    1259                                                   sums_wsqrs_ws_l,                &
     1271                                                  sums_wsss_ws_l,                 &
    12601272                                                  sums_ls_l
    12611273                                             
  • palm/trunk/SOURCE/netcdf_interface_mod.f90

    r1958 r1960  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Additional labels and units for timeseries output of passive scalar-related
     22! quantities
    2223!
    2324! Former revisions:
     
    174175             'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'number2' /)
    175176
    176     INTEGER(iwp) ::  dots_num  = 31  !< number of timeseries defined by default
    177     INTEGER(iwp) ::  dots_soil = 24  !< starting index for soil-timeseries
    178     INTEGER(iwp) ::  dots_rad  = 32  !< starting index for radiation-timeseries
     177    INTEGER(iwp) ::  dots_num  = 33  !< number of timeseries defined by default
     178    INTEGER(iwp) ::  dots_soil = 26  !< starting index for soil-timeseries
     179    INTEGER(iwp) ::  dots_rad  = 34  !< starting index for radiation-timeseries
    179180
    180181    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label =                    &
     
    186187             'wpt          ', 'pt(0)        ', 'pt(z_mo)     ',                &
    187188             'w"u"0        ', 'w"v"0        ', 'w"q"0        ',                &
    188              'ol           ', 'q*           ', 'ghf_eb       ',                &
    189              'shf_eb       ', 'qsws_eb      ', 'qsws_liq_eb  ',                &
    190              'qsws_soil_eb ', 'qsws_veg_eb  ', 'r_a          ',                &
    191              'r_s          ', 'rad_net      ', 'rad_lw_in    ',                &
    192              'rad_lw_out   ', 'rad_sw_in    ', 'rad_sw_out   ',                &
    193              'rrtm_aldif   ', 'rrtm_aldir   ', 'rrtm_asdif   ',                &
    194              'rrtm_asdir   ',                                                  &
    195              ( 'unknown      ', i9 = 1, dots_max-40 ) /)
     189             'ol           ', 'q*           ', 'w"s"         ',                &
     190             's*           ', 'ghf_eb       ', 'shf_eb       ',                &
     191             'qsws_eb      ', 'qsws_liq_eb  ', 'qsws_soil_eb ',                &
     192             'qsws_veg_eb  ', 'r_a          ', 'r_s          ',                &
     193             'rad_net      ', 'rad_lw_in    ', 'rad_lw_out   ',                &
     194             'rad_sw_in    ', 'rad_sw_out   ', 'rrtm_aldif   ',                &
     195             'rrtm_aldir   ', 'rrtm_asdif   ', 'rrtm_asdir   ',                &                                               
     196             ( 'unknown      ', i9 = 1, dots_max-42 ) /)
    196197
    197198    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_unit =                     &
     
    203204             'K m/s        ', 'K            ', 'K            ',                &
    204205             'm2/s2        ', 'm2/s2        ', 'kg m/s       ',                &
    205              'm            ', 'kg/kg        ', '             ',                &
     206             'm            ', 'kg/kg        ', 'kg m/(kg s)  ',                &
     207             'kg/kg        ', '             ', '             ',                &
    206208             '             ', '             ', '             ',                &
    207              '             ', 'W/m2         ', 's/m          ',                &
    208              '             ', 'W/m2         ', 'W/m2         ',                &
     209             'W/m2         ', 's/m          ', '             ',                &
    209210             'W/m2         ', 'W/m2         ', 'W/m2         ',                &
     211             'W/m2         ', 'W/m2         ', '             ',                &
    210212             '             ', '             ', '             ',                &
    211              '             ',                                                  &
    212              ( 'unknown      ', i9 = 1, dots_max-40 ) /)
     213             ( 'unknown      ', i9 = 1, dots_max-42 ) /)
    213214
    214215    CHARACTER (LEN=9), DIMENSION(300) ::  dopr_unit = 'unknown'
  • palm/trunk/SOURCE/palm.f90

    r1933 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Separate humidity and passive scalar
    2222!
    2323! Former revisions:
     
    361361          ENDIF
    362362          IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp )
    363           IF (humidity  .OR.  passive_scalar)  THEN
    364              CALL exchange_horiz( q, nbgp )
    365           ENDIF
     363          IF ( humidity       )  CALL exchange_horiz( q, nbgp )
     364          IF ( passive_scalar )  CALL exchange_horiz( s, nbgp )
    366365       ENDIF
    367366
  • palm/trunk/SOURCE/parin.f90

    r1958 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Allocation of s_init
    2222!
    2323! Former revisions:
     
    218218
    219219    USE arrays_3d,                                                             &
    220         ONLY:  pt_init, q_init, ref_state, sa_init, ug, u_init, v_init,        &
     220        ONLY:  pt_init, q_init, ref_state, s_init, sa_init, ug, u_init, v_init,&
    221221               vg
    222222
     
    597597!--          ATTENTION: in case of changes to the following statement please
    598598!--                  also check the allocate statement in routine read_var_list
    599              ALLOCATE( pt_init(0:nz+1), q_init(0:nz+1),                       &
     599             ALLOCATE( pt_init(0:nz+1), q_init(0:nz+1), s_init(0:nz+1),        &
    600600                       ref_state(0:nz+1), sa_init(0:nz+1), ug(0:nz+1),         &
    601601                       u_init(0:nz+1), v_init(0:nz+1), vg(0:nz+1),             &
  • palm/trunk/SOURCE/plant_canopy_model_mod.f90

    r1954 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Separate humidity and passive scalar
    2222!
    2323! Former revisions:
     
    8686 
    8787    USE arrays_3d,                                                             &
    88         ONLY:  dzu, dzw, e, q, shf, tend, u, v, w, zu, zw
     88        ONLY:  dzu, dzw, e, q, s, shf, tend, u, v, w, zu, zw
    8989
    9090    USE indices,                                                               &
     
    818818
    819819!
    820 !--       scalar concentration
     820!--       humidity
    821821          CASE ( 5 )
    822822             DO  i = nxl, nxr
     
    866866                ENDDO
    867867             ENDDO
     868!
     869!--       scalar concentration
     870          CASE ( 7 )
     871             DO  i = nxl, nxr
     872                DO  j = nys, nyn
     873                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
     874                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
     875                      tend(k,j,i) = tend(k,j,i) -                              &
     876                                       lsec *                                  &
     877                                       lad_s(kk,j,i) *                         &
     878                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
     879                                                          u(k,j,i+1) )         &
     880                                             )**2 +                            &
     881                                             ( 0.5_wp * ( v(k,j,i) +           &
     882                                                          v(k,j+1,i) )         &
     883                                             )**2 +                            &
     884                                             ( 0.5_wp * ( w(k-1,j,i) +         &
     885                                                          w(k,j,i) )           &
     886                                             )**2                              &
     887                                           ) *                                 &
     888                                       ( s(k,j,i) - lsc )
     889                   ENDDO
     890                ENDDO
     891             ENDDO   
     892
    868893
    869894
     
    11461171
    11471172!
    1148 !--       scalar concentration
     1173!--       humidity
    11491174          CASE ( 5 )
    11501175             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
     
    11861211                                 e(k,j,i)
    11871212             ENDDO
     1213             
     1214!
     1215!--       scalar concentration
     1216          CASE ( 7 )
     1217             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
     1218                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
     1219                tend(k,j,i) = tend(k,j,i) -                                    &
     1220                                 lsec *                                        &
     1221                                 lad_s(kk,j,i) *                               &
     1222                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
     1223                                                    u(k,j,i+1) )               &
     1224                                       )**2  +                                 &
     1225                                       ( 0.5_wp * ( v(k,j,i) +                 &
     1226                                                    v(k,j+1,i) )               &
     1227                                       )**2 +                                  &
     1228                                       ( 0.5_wp * ( w(k-1,j,i) +               &
     1229                                                    w(k,j,i) )                 &
     1230                                       )**2                                    &
     1231                                     ) *                                       &
     1232                                 ( s(k,j,i) - lsc )
     1233             ENDDO               
    11881234
    11891235       CASE DEFAULT
  • palm/trunk/SOURCE/prognostic_equations.f90

    r1917 r1960  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Separate humidity and passive scalar
    2222!
    2323! Former revisions:
     
    192192    USE arrays_3d,                                                             &
    193193        ONLY:  diss_l_e, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr,            &
    194                diss_l_sa, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q,            &
    195                diss_s_qr, diss_s_sa, e, e_p, flux_s_e, flux_s_nr, flux_s_pt,   &
    196                flux_s_q, flux_s_qr, flux_s_sa, flux_l_e, flux_l_nr,            &
    197                flux_l_pt, flux_l_q, flux_l_qr, flux_l_sa, nr, nr_p, nrsws,     &
    198                nrswst, pt, ptdf_x, ptdf_y, pt_init, pt_p, prho, q, q_init,     &
    199                q_p, qsws, qswst, qr, qr_p, qrsws, qrswst, rdf, rdf_sc,         &
    200                ref_state, rho, sa, sa_init, sa_p, saswsb, saswst, shf, tend,   &
    201                te_m, tnr_m, tpt_m, tq_m, tqr_m, tsa_m, tswst, tu_m, tv_m,      &
     194               diss_l_s, diss_l_sa, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q,  &
     195               diss_s_qr, diss_s_s, diss_s_sa, e, e_p, flux_s_e, flux_s_nr,    &
     196               flux_s_pt, flux_s_q, flux_s_qr, flux_s_s, flux_s_sa, flux_l_e,  &
     197               flux_l_nr, flux_l_pt, flux_l_q, flux_l_qr, flux_l_s, flux_l_sa, &
     198               nr, nr_p, nrsws, nrswst, pt, ptdf_x, ptdf_y, pt_init, pt_p,     &
     199               prho, q, q_init, q_p, qsws, qswst, qr, qr_p, qrsws, qrswst, rdf,&
     200               rdf_sc, ref_state, rho, s, s_init, s_p, sa, sa_init, sa_p,      &
     201               saswsb, saswst, shf, ssws, sswst, tend,                         &
     202               te_m, tnr_m, tpt_m, tq_m, tqr_m, ts_m, tsa_m, tswst, tu_m, tv_m,&
    202203               tw_m, u, ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
    203204       
     
    216217               use_upstream_for_tke, wall_heatflux,                            &
    217218               wall_nrflux, wall_qflux, wall_qflux, wall_qflux, wall_qrflux,   &
    218                wall_salinityflux, ws_scheme_mom, ws_scheme_sca
     219               wall_salinityflux, wall_sflux, ws_scheme_mom, ws_scheme_sca
    219220
    220221    USE cpulog,                                                                &
     
    724725
    725726!
    726 !--       If required, compute prognostic equation for total water content /
    727 !--       scalar
    728           IF ( humidity  .OR.  passive_scalar )  THEN
     727!--       If required, compute prognostic equation for total water content
     728          IF ( humidity )  THEN
    729729
    730730!
     
    745745
    746746!
    747 !--          Sink or source of scalar concentration due to canopy elements
     747!--          Sink or source of humidity due to canopy elements
    748748             IF ( plant_canopy )  CALL pcm_tendency( i, j, 5 )
    749749
     
    881881
    882882          ENDIF
    883 
     883         
     884!
     885!--       If required, compute prognostic equation for scalar
     886          IF ( passive_scalar )  THEN
     887!
     888!--          Tendency-terms for total water content / scalar
     889             tend(:,j,i) = 0.0_wp
     890             IF ( timestep_scheme(1:5) == 'runge' ) &
     891             THEN
     892                IF ( ws_scheme_sca )  THEN
     893                   CALL advec_s_ws( i, j, s, 's', flux_s_s, &
     894                                diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
     895                ELSE
     896                   CALL advec_s_pw( i, j, s )
     897                ENDIF
     898             ELSE
     899                CALL advec_s_up( i, j, s )
     900             ENDIF
     901             CALL diffusion_s( i, j, s, ssws, sswst, wall_sflux )
     902
     903!
     904!--          Sink or source of scalar concentration due to canopy elements
     905             IF ( plant_canopy )  CALL pcm_tendency( i, j, 7 )
     906
     907!
     908!--          Large scale advection, still need to be extended for scalars
     909!              IF ( large_scale_forcing )  THEN
     910!                 CALL ls_advec( i, j, simulated_time, 's' )
     911!              ENDIF
     912
     913!
     914!--          Nudging, still need to be extended for scalars
     915!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' )
     916
     917!
     918!--          If required compute influence of large-scale subsidence/ascent.
     919!--          Note, the last argument is of no meaning in this case, as it is
     920!--          only used in conjunction with large_scale_forcing, which is to
     921!--          date not implemented for scalars.
     922             IF ( large_scale_subsidence  .AND.                                &
     923                  .NOT. use_subsidence_tendencies  .AND.                       &
     924                  .NOT. large_scale_forcing )  THEN
     925                CALL subsidence( i, j, tend, s, s_init, 3 )
     926             ENDIF
     927
     928             CALL user_actions( i, j, 's-tendency' )
     929
     930!
     931!--          Prognostic equation for scalar
     932             DO  k = nzb_s_inner(j,i)+1, nzt
     933                s_p(k,j,i) = s(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +       &
     934                                                  tsc(3) * ts_m(k,j,i) )       &
     935                                      - tsc(5) * rdf_sc(k) *                   &
     936                                        ( s(k,j,i) - s_init(k) )
     937                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
     938             ENDDO
     939
     940!
     941!--          Calculate tendencies for the next Runge-Kutta step
     942             IF ( timestep_scheme(1:5) == 'runge' )  THEN
     943                IF ( intermediate_timestep_count == 1 )  THEN
     944                   DO  k = nzb_s_inner(j,i)+1, nzt
     945                      ts_m(k,j,i) = tend(k,j,i)
     946                   ENDDO
     947                ELSEIF ( intermediate_timestep_count < &
     948                         intermediate_timestep_count_max )  THEN
     949                   DO  k = nzb_s_inner(j,i)+1, nzt
     950                      ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
     951                                     5.3125_wp * ts_m(k,j,i)
     952                   ENDDO
     953                ENDIF
     954             ENDIF
     955
     956          ENDIF         
    884957!
    885958!--       If required, compute prognostic equation for turbulent kinetic
     
    14321505
    14331506!
    1434 !-- If required, compute prognostic equation for total water content / scalar
    1435     IF ( humidity  .OR.  passive_scalar )  THEN
    1436 
    1437        CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
     1507!-- If required, compute prognostic equation for total water content
     1508    IF ( humidity )  THEN
     1509
     1510       CALL cpu_log( log_point(29), 'q-equation', 'start' )
    14381511
    14391512!
     
    14701543       
    14711544!
    1472 !--    Sink or source of scalar concentration due to canopy elements
     1545!--    Sink or source of humidity due to canopy elements
    14731546       IF ( plant_canopy ) CALL pcm_tendency( 5 )
    14741547
     
    14931566
    14941567!
    1495 !--    Prognostic equation for total water content / scalar
     1568!--    Prognostic equation for total water content
    14961569       DO  i = nxl, nxr
    14971570          DO  j = nys, nyn
     
    15291602       ENDIF
    15301603
    1531        CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
     1604       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
    15321605
    15331606!
     
    16861759
    16871760    ENDIF
    1688 
     1761!
     1762!-- If required, compute prognostic equation for scalar
     1763    IF ( passive_scalar )  THEN
     1764
     1765       CALL cpu_log( log_point(66), 's-equation', 'start' )
     1766
     1767!
     1768!--    Scalar/q-tendency terms with communication
     1769       sbt = tsc(2)
     1770       IF ( scalar_advec == 'bc-scheme' )  THEN
     1771
     1772          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     1773!
     1774!--          Bott-Chlond scheme always uses Euler time step. Thus:
     1775             sbt = 1.0_wp
     1776          ENDIF
     1777          tend = 0.0_wp
     1778          CALL advec_s_bc( s, 's' )
     1779
     1780       ENDIF
     1781
     1782!
     1783!--    Scalar/q-tendency terms with no communication
     1784       IF ( scalar_advec /= 'bc-scheme' )  THEN
     1785          tend = 0.0_wp
     1786          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1787             IF ( ws_scheme_sca )  THEN
     1788                CALL advec_s_ws( s, 's' )
     1789             ELSE
     1790                CALL advec_s_pw( s )
     1791             ENDIF
     1792          ELSE
     1793             CALL advec_s_up( s )
     1794          ENDIF
     1795       ENDIF
     1796
     1797       CALL diffusion_s( s, ssws, sswst, wall_sflux )
     1798       
     1799!
     1800!--    Sink or source of humidity due to canopy elements
     1801       IF ( plant_canopy ) CALL pcm_tendency( 7 )
     1802
     1803!
     1804!--    Large scale advection. Not implemented for scalars so far.
     1805!        IF ( large_scale_forcing )  THEN
     1806!           CALL ls_advec( simulated_time, 'q' )
     1807!        ENDIF
     1808
     1809!
     1810!--    Nudging. Not implemented for scalars so far.
     1811!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
     1812
     1813!
     1814!--    If required compute influence of large-scale subsidence/ascent.
     1815!--    Not implemented for scalars so far.
     1816       IF ( large_scale_subsidence  .AND.                                      &
     1817            .NOT. use_subsidence_tendencies  .AND.                             &
     1818            .NOT. large_scale_forcing )  THEN
     1819         CALL subsidence( tend, s, s_init, 3 )
     1820       ENDIF
     1821
     1822       CALL user_actions( 's-tendency' )
     1823
     1824!
     1825!--    Prognostic equation for total water content
     1826       DO  i = nxl, nxr
     1827          DO  j = nys, nyn
     1828             DO  k = nzb_s_inner(j,i)+1, nzt
     1829                s_p(k,j,i) = s(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +          &
     1830                                                  tsc(3) * ts_m(k,j,i) )       &
     1831                                      - tsc(5) * rdf_sc(k) *                   &
     1832                                        ( s(k,j,i) - s_init(k) )
     1833                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
     1834             ENDDO
     1835          ENDDO
     1836       ENDDO
     1837
     1838!
     1839!--    Calculate tendencies for the next Runge-Kutta step
     1840       IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1841          IF ( intermediate_timestep_count == 1 )  THEN
     1842             DO  i = nxl, nxr
     1843                DO  j = nys, nyn
     1844                   DO  k = nzb_s_inner(j,i)+1, nzt
     1845                      ts_m(k,j,i) = tend(k,j,i)
     1846                   ENDDO
     1847                ENDDO
     1848             ENDDO
     1849          ELSEIF ( intermediate_timestep_count < &
     1850                   intermediate_timestep_count_max )  THEN
     1851             DO  i = nxl, nxr
     1852                DO  j = nys, nyn
     1853                   DO  k = nzb_s_inner(j,i)+1, nzt
     1854                      ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * ts_m(k,j,i)
     1855                   ENDDO
     1856                ENDDO
     1857             ENDDO
     1858          ENDIF
     1859       ENDIF
     1860
     1861       CALL cpu_log( log_point(66), 's-equation', 'stop' )
     1862
     1863    ENDIF
    16891864!
    16901865!-- If required, compute prognostic equation for turbulent kinetic
     
    22262401
    22272402!
    2228 !-- If required, compute prognostic equation for total water content / scalar
    2229     IF ( humidity  .OR.  passive_scalar )  THEN
    2230 
    2231        CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
     2403!-- If required, compute prognostic equation for total water content
     2404    IF ( humidity )  THEN
     2405
     2406       CALL cpu_log( log_point(29), 'q-equation', 'start' )
    22322407
    22332408!
     
    23072482       ENDDO
    23082483
    2309        CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
     2484       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
    23102485
    23112486!
     
    24302605    ENDIF
    24312606
     2607!
     2608!-- If required, compute prognostic equation for scalar
     2609    IF ( passive_scalar )  THEN
     2610
     2611       CALL cpu_log( log_point(66), 's-equation', 'start' )
     2612
     2613!
     2614!--    Scalar/q-tendency terms with communication
     2615       sbt = tsc(2)
     2616       IF ( scalar_advec == 'bc-scheme' )  THEN
     2617
     2618          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     2619!
     2620!--          Bott-Chlond scheme always uses Euler time step. Thus:
     2621             sbt = 1.0_wp
     2622          ENDIF
     2623          tend = 0.0_wp
     2624          CALL advec_s_bc( s, 's' )
     2625
     2626       ENDIF
     2627
     2628!
     2629!--    Scalar/q-tendency terms with no communication
     2630       IF ( scalar_advec /= 'bc-scheme' )  THEN
     2631          tend = 0.0_wp
     2632          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     2633             IF ( ws_scheme_sca )  THEN
     2634                CALL advec_s_ws( s, 's' )
     2635             ELSE
     2636                CALL advec_s_pw( s )
     2637             ENDIF
     2638          ELSE
     2639             CALL advec_s_up( s )
     2640          ENDIF
     2641       ENDIF
     2642
     2643       CALL diffusion_s( s, ssws, sswst, wall_sflux )
     2644
     2645!
     2646!--    Sink or source of scalar concentration due to canopy elements
     2647       IF ( plant_canopy ) CALL pcm_tendency( 7 )
     2648
     2649!
     2650!--    Large scale advection. Not implemented so far.
     2651!        IF ( large_scale_forcing )  THEN
     2652!           CALL ls_advec( simulated_time, 's' )
     2653!        ENDIF
     2654
     2655!
     2656!--    Nudging. Not implemented so far.
     2657!        IF ( nudging )  CALL nudge( simulated_time, 's' )
     2658
     2659!
     2660!--    If required compute influence of large-scale subsidence/ascent.
     2661!--    Not implemented so far.
     2662       IF ( large_scale_subsidence  .AND.                                      &
     2663            .NOT. use_subsidence_tendencies  .AND.                             &
     2664            .NOT. large_scale_forcing )  THEN
     2665         CALL subsidence( tend, s, s_init, 3 )
     2666       ENDIF
     2667
     2668       CALL user_actions( 's-tendency' )
     2669
     2670!
     2671!--    Prognostic equation for total water content / scalar
     2672       DO  i = i_left, i_right
     2673          DO  j = j_south, j_north
     2674             DO  k = nzb_s_inner(j,i)+1, nzt
     2675                s_p(k,j,i) = s(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +          &
     2676                                                  tsc(3) * ts_m(k,j,i) )       &
     2677                                      - tsc(5) * rdf_sc(k) *                   &
     2678                                        ( s(k,j,i) - s_init(k) )
     2679                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
     2680!
     2681!--             Tendencies for the next Runge-Kutta step
     2682                IF ( runge_step == 1 )  THEN
     2683                   ts_m(k,j,i) = tend(k,j,i)
     2684                ELSEIF ( runge_step == 2 )  THEN
     2685                   ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * ts_m(k,j,i)
     2686                ENDIF
     2687             ENDDO
     2688          ENDDO
     2689       ENDDO
     2690
     2691       CALL cpu_log( log_point(66), 's-equation', 'stop' )
     2692
     2693    ENDIF
    24322694!
    24332695!-- If required, compute prognostic equation for turbulent kinetic
  • palm/trunk/SOURCE/read_3d_binary.f90

    r1852 r1960  
    107107        ONLY:  e, kh, km, ol, p, pt, q, ql, qc, nr, nrs, nrsws, nrswst,        &
    108108               prr, precipitation_amount, qr,                                  &
    109                qrs, qrsws, qrswst, qs, qsws, qswst, sa, saswsb, saswst,        &
    110                rif_wall, shf, ts, tswst, u, u_m_l, u_m_n, u_m_r, u_m_s, us,    &
    111                usws, uswst, v, v_m_l, v_m_n, v_m_r, v_m_s, vpt, vsws, vswst,   &
    112                w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h, z0q
     109               qrs, qrsws, qrswst, qs, qsws, qswst, s, sa, saswsb, saswst,     &
     110               ss, ssws, sswst, rif_wall, shf, ss, ts, tswst, u, u_m_l, u_m_n, &
     111               u_m_r, u_m_s, us, usws, uswst, v, v_m_l, v_m_n, v_m_r, v_m_s,   &
     112               vpt, vsws, vswst, w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h, z0q
    113113
    114114    USE averaging
     
    985985                   rif_wall(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp,:) = &
    986986                            tmp_4d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp,:)
     987                           
     988                CASE ( 's' )
     989                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     990                   s(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     991                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    987992
    988993                CASE ( 's_av' )
     
    10671072                      ENDIF
    10681073                   ENDIF
    1069 
     1074                   
     1075                CASE ( 'ss' )
     1076                   IF ( k == 1 )  READ ( 13 )  tmp_2d
     1077                   ss(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     1078                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     1079
     1080                CASE ( 'ssws' )
     1081                   IF ( k == 1 )  READ ( 13 )  tmp_2d
     1082                   ssws(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     1083                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     1084
     1085                CASE ( 'ssws_av' )
     1086                   IF ( .NOT. ALLOCATED( ssws_av ) )  THEN
     1087                      ALLOCATE( ssws_av(nysg:nyng,nxlg:nxrg) )
     1088                   ENDIF 
     1089                   IF ( k == 1 )  READ ( 13 )  tmp_2d
     1090                   ssws_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     1091                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     1092
     1093                CASE ( 'sswst' )
     1094                   IF ( k == 1 )  READ ( 13 )  tmp_2d
     1095                   sswst(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     1096                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     1097                 
    10701098                CASE ( 'ts' )
    10711099                   IF ( k == 1 )  READ ( 13 )  tmp_2d
  • palm/trunk/SOURCE/read_var_list.f90

    r1958 r1960  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Separate humidity and passive scalar
     22! Remove unused variables from ONLY list
    2223!
    2324! Former revisions:
     
    175176    USE arrays_3d,                                                             &
    176177        ONLY:  inflow_damping_factor, mean_inflow_profiles, pt_init,           &
    177                q_init, ref_state, sa_init, u_init, ug, v_init, vg
     178               q_init, ref_state, s_init, sa_init, u_init, ug, v_init, vg
    178179
    179180    USE control_parameters
     
    292293       ALLOCATE( ug(0:nz+1), u_init(0:nz+1), vg(0:nz+1),                       &
    293294                 v_init(0:nz+1), pt_init(0:nz+1), q_init(0:nz+1),              &
    294                  ref_state(0:nz+1), sa_init(0:nz+1),                           &
     295                 ref_state(0:nz+1), s_init(0:nz+1), sa_init(0:nz+1),           &
    295296                 hom(0:nz+1,2,pr_palm+max_pr_user,0:statistic_regions),        &
    296297                 hom_sum(0:nz+1,pr_palm+max_pr_user,0:statistic_regions) )
     
    486487          CASE ( 'mean_inflow_profiles' )
    487488             IF ( .NOT. ALLOCATED( mean_inflow_profiles ) )  THEN
    488                 ALLOCATE( mean_inflow_profiles(0:nz+1,6) )
     489                ALLOCATE( mean_inflow_profiles(0:nz+1,7) )
    489490             ENDIF
    490491             READ ( 13 )  mean_inflow_profiles
     
    597598          CASE ( 'run_coupled' )
    598599             READ ( 13 )  run_coupled
     600          CASE ( 's_init' )
     601             READ ( 13 )  s_init
     602          CASE ( 's_surface' )
     603             READ ( 13 )  s_surface
     604          CASE ( 's_surface_initial_change' )
     605             READ ( 13 )  s_surface_initial_change
     606          CASE ( 's_vertical_gradient' )
     607             READ ( 13 )  s_vertical_gradient
     608          CASE ( 's_vertical_gradient_level' )
     609             READ ( 13 )  s_vertical_gradient_level
     610          CASE ( 's_vertical_gradient_level_ind' )
     611             READ ( 13 )  s_vertical_gradient_level_ind
    599612          CASE ( 'sa_init' )
    600613             READ ( 13 )  sa_init
     
    617630          CASE ( 'surface_waterflux' )
    618631             READ ( 13 )  surface_waterflux     
    619           CASE ( 's_surface' )
    620              READ ( 13 )  s_surface
    621           CASE ( 's_surface_initial_change' )
    622              READ ( 13 )  s_surface_initial_change
    623           CASE ( 's_vertical_gradient' )
    624              READ ( 13 )  s_vertical_gradient
    625           CASE ( 's_vertical_gradient_level' )
    626              READ ( 13 )  s_vertical_gradient_level
    627632          CASE ( 'time_coupling' )
    628633             READ ( 13 )  time_coupling
     
    789794
    790795    USE arrays_3d,                                                             &
    791         ONLY:  inflow_damping_factor, mean_inflow_profiles, pt_init,           &
    792                q_init, ref_state, sa_init, u_init, ug, v_init, vg
     796        ONLY:  inflow_damping_factor, mean_inflow_profiles, ref_state, ug, vg
    793797
    794798    USE control_parameters
  • palm/trunk/SOURCE/spectra_mod.f90

    r1834 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Additional default spectra for passive scalar
    2222!
    2323! Former revisions:
     
    480480
    481481       USE arrays_3d,                                                          &
    482            ONLY:  d, pt, q, u, v, w
     482           ONLY:  d, pt, q, s, u, v, w
    483483
    484484       USE indices,                                                            &
     
    532532          pr = 41
    533533          d(nzb+1:nzt,nys:nyn,nxl:nxr) = q(nzb+1:nzt,nys:nyn,nxl:nxr)
     534         
     535       CASE ( 's' )
     536          pr = 117
     537          d(nzb+1:nzt,nys:nyn,nxl:nxr) = s(nzb+1:nzt,nys:nyn,nxl:nxr)
    534538       
    535539       CASE DEFAULT
  • palm/trunk/SOURCE/sum_up_3d_data.f90

    r1950 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Scalar surface flux added
    2222!
    2323! Former revisions:
     
    101101    USE arrays_3d,                                                             &
    102102        ONLY:  dzw, e, nr, ol, p, pt, precipitation_rate, q, qc, ql, ql_c,     &
    103                ql_v, qr, qsws, rho, sa, shf, ts, u, us, v, vpt, w, z0, z0h, z0q
     103               ql_v, qr, qsws, rho, sa, shf, ssws, ts, u, us, v, vpt, w, z0,   &
     104               z0h, z0q
    104105
    105106    USE averaging,                                                             &
     
    107108               precipitation_rate_av, pt_av, q_av, qc_av, ql_av, ql_c_av,      &
    108109               ql_v_av, ql_vp_av, qr_av, qsws_av, qv_av, rho_av, s_av, sa_av,  &
    109                shf_av, ts_av, u_av, us_av, v_av, vpt_av, w_av, z0_av, z0h_av, &
    110                z0q_av
     110               shf_av, ssws_av, ts_av, u_av, us_av, v_av, vpt_av, w_av, z0_av, &
     111               z0h_av, z0q_av
    111112
    112113    USE cloud_parameters,                                                      &
     
    949950             ENDDO
    950951
     952          CASE ( 'ssws*' )
     953             DO  i = nxlg, nxrg
     954                DO  j = nysg, nyng
     955                   ssws_av(j,i) = ssws_av(j,i) + ssws(j,i)
     956                ENDDO
     957             ENDDO
     958
    951959          CASE ( 't*' )
    952960             DO  i = nxlg, nxrg
  • palm/trunk/SOURCE/surface_layer_fluxes_mod.f90

    r1930 r1960  
    1919! Current revisions:
    2020! ------------------
    21 !  
     21! Treat humidity and passive scalar separately
    2222!
    2323! Former revisions:
     
    141141    USE arrays_3d,                                                             &
    142142        ONLY:  e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, &
    143                shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h, z0q
     143               s, shf, ss, ssws, ts, u, us, usws, v, vpt, vsws, zu, zw, z0,    &
     144               z0h, z0q
    144145
    145146    USE cloud_parameters,                                                      &
     
    152153
    153154    USE control_parameters,                                                    &
    154         ONLY:  cloud_physics, constant_heatflux, constant_waterflux,           &
    155                coupling_mode, g, humidity, ibc_e_b, ibc_pt_b,                  &
    156                initializing_actions, kappa, intermediate_timestep_count,       &
     155        ONLY:  cloud_physics, constant_heatflux, constant_scalarflux,          &
     156               constant_waterflux, coupling_mode, g, humidity, ibc_e_b,        &
     157               ibc_pt_b, initializing_actions, kappa,                          &
     158               intermediate_timestep_count,                                    &
    157159               intermediate_timestep_count_max, large_scale_forcing, lsf_surf, &
    158160               message_string, microphysics_seifert, most_method, neutral,     &
     
    859861!
    860862!--    If required compute q*
    861        IF ( humidity  .OR.  passive_scalar )  THEN
     863       IF ( humidity )  THEN
    862864          IF ( constant_waterflux )  THEN
    863865!
     
    922924          ENDIF
    923925       ENDIF
     926       
     927!
     928!--    If required compute s*
     929       IF ( passive_scalar )  THEN
     930          IF ( constant_scalarflux )  THEN
     931!
     932!--          For a given water flux in the surface layer
     933             !$OMP PARALLEL DO
     934             !$acc kernels loop private( j )
     935             DO  i = nxlg, nxrg
     936                DO  j = nysg, nyng
     937                   ss(j,i) = -ssws(j,i) / ( us(j,i) + 1E-30_wp )
     938                ENDDO
     939             ENDDO
     940          ENDIF
     941       ENDIF
    924942
    925943
     
    10431061
    10441062!
    1045 !-- Compute the vertical water/scalar flux
    1046        IF (  .NOT.  constant_waterflux  .AND.  ( humidity  .OR.                &
    1047              passive_scalar )  .AND.  ( simulated_time <= skip_time_do_lsm     &
     1063!--    Compute the vertical water flux
     1064       IF (  .NOT.  constant_waterflux  .AND.  humidity  .AND.                 &
     1065             ( simulated_time <= skip_time_do_lsm                              &
    10481066            .OR.  .NOT.  land_surface ) )  THEN
    10491067          !$OMP PARALLEL DO
     
    10571075
    10581076       ENDIF
     1077       
     1078!
     1079!--    Compute the vertical scalar flux
     1080       IF (  .NOT.  constant_scalarflux  .AND.  passive_scalar  .AND.          &
     1081             ( simulated_time <= skip_time_do_lsm                              &
     1082            .OR.  .NOT.  land_surface ) )  THEN
     1083          !$OMP PARALLEL DO
     1084          !$acc kernels loop independent present( qs, qsws, us )
     1085          DO  i = nxlg, nxrg
     1086             !$acc loop independent
     1087             DO  j = nysg, nyng
     1088                ssws(j,i) = -ss(j,i) * us(j,i)
     1089             ENDDO
     1090          ENDDO
     1091
     1092       ENDIF       
    10591093
    10601094!
  • palm/trunk/SOURCE/swap_timelevel.f90

    r1823 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Separate humidity and passive scalar
    2222!
    2323! Former revisions:
     
    8383#if defined( __nopointer )
    8484    USE arrays_3d,                                                             &
    85         ONLY:  e, e_p, nr, nr_p, pt, pt_p, q, q_p, qr, qr_p, sa, sa_p, u, u_p, &
    86                v, v_p, w, w_p
     85        ONLY:  e, e_p, nr, nr_p, pt, pt_p, q, q_p, qr, qr_p, s, s_p, sa, sa_p, &
     86               u, u_p, v, v_p, w, w_p
    8787#else
    8888    USE arrays_3d,                                                             &
    8989        ONLY:  e, e_1, e_2, e_p, nr, nr_1, nr_2, nr_p, pt, pt_1, pt_2, pt_p, q,&
    90                q_1, q_2, q_p, qr, qr_1, qr_2, qr_p, sa, sa_1, sa_2, sa_p, u,   &
    91                u_1, u_2, u_p, v, v_1, v_2, v_p, w, w_1, w_2, w_p
     90               q_1, q_2, q_p, qr, qr_1, qr_2, qr_p, s, s_1, s_2, s_p, sa, sa_1,&
     91               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
    9292
    9393#endif
     
    162162       sa = sa_p
    163163    ENDIF
    164     IF ( humidity  .OR.  passive_scalar )  THEN
     164    IF ( humidity )  THEN
    165165       q = q_p             
    166166       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     
    169169       ENDIF
    170170    ENDIF
     171    IF ( passive_scalar )  s = s_p             
    171172
    172173    IF ( land_surface )  THEN
     
    195196             sa => sa_1;  sa_p => sa_2
    196197          ENDIF
    197           IF ( humidity  .OR.  passive_scalar )  THEN
     198          IF ( humidity )  THEN
    198199             q => q_1;    q_p => q_2
    199200             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     
    202203             ENDIF
    203204          ENDIF
     205          IF ( passive_scalar )  THEN
     206             s => s_1;    s_p => s_2
     207          ENDIF
    204208
    205209          IF ( land_surface )  THEN
     
    223227             sa => sa_2;  sa_p => sa_1
    224228          ENDIF
    225           IF ( humidity  .OR.  passive_scalar )  THEN
     229          IF ( humidity )  THEN
    226230             q => q_2;    q_p => q_1
    227231             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     
    230234             ENDIF
    231235          ENDIF
     236          IF ( passive_scalar )  THEN
     237             s => s_2;    s_p => s_1
     238          ENDIF
    232239
    233240          IF ( land_surface )  THEN
  • palm/trunk/SOURCE/time_integration.f90

    r1958 r1960  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Separate humidity and passive scalar
    2222!
    2323! Former revisions:
     
    225225    USE arrays_3d,                                                             &
    226226        ONLY:  diss, dzu, e, e_p, nr_p, prho, pt, pt_p, pt_init, q_init, q,    &
    227                ql, ql_c, ql_v, ql_vp, qr_p, q_p, ref_state, rho, sa_p, tend,  &
    228                u, u_p, v, vpt, v_p, w, w_p
     227               ql, ql_c, ql_v, ql_vp, qr_p, q_p, ref_state, rho, s, s_p, sa_p, &
     228               tend, u, u_p, v, vpt, v_p, w, w_p
    229229
    230230    USE calc_mean_profile_mod,                                                 &
     
    532532                   CALL exchange_horiz( sa_p, nbgp )
    533533                   CALL exchange_horiz( rho, nbgp )
    534                   CALL exchange_horiz( prho, nbgp )
    535                 ENDIF
    536                 IF (humidity  .OR.  passive_scalar)  THEN
     534                   CALL exchange_horiz( prho, nbgp )
     535                ENDIF
     536                IF ( humidity )  THEN
    537537                   CALL exchange_horiz( q_p, nbgp )
    538538                   IF ( cloud_physics .AND. microphysics_seifert )  THEN
     
    551551                   CALL exchange_horiz( diss, nbgp )
    552552                ENDIF
     553                IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )
    553554
    554555                IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
     
    601602                  CALL exchange_horiz( prho, nbgp )
    602603                ENDIF
    603                 IF (humidity  .OR.  passive_scalar)  THEN
     604                IF ( humidity )  THEN
    604605                   CALL exchange_horiz( q_p, nbgp )
    605606                   IF ( cloud_physics .AND. microphysics_seifert )  THEN
     
    618619                   CALL exchange_horiz( diss, nbgp )
    619620                ENDIF
     621                IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )
    620622
    621623                IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
     
    690692                CALL exchange_horiz( prho, nbgp )
    691693             ENDIF
    692              IF (humidity  .OR.  passive_scalar)  THEN
     694             IF ( humidity )  THEN
    693695                CALL exchange_horiz( q_p, nbgp )
    694696                IF ( cloud_physics .AND. microphysics_seifert )  THEN
     
    707709                CALL exchange_horiz( diss, nbgp )
    708710             ENDIF
     711             IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )     
    709712
    710713             IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
     
    746749                CALL exchange_horiz( v, nbgp )
    747750                CALL exchange_horiz( w, nbgp )
    748                 IF ( .NOT. neutral )  THEN
    749                    CALL exchange_horiz( pt, nbgp )
    750                 ENDIF
    751                 IF ( humidity  .OR.  passive_scalar )  THEN
    752                    CALL exchange_horiz( q, nbgp )
    753                 ENDIF
    754                 IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp )
     751                IF ( .NOT. neutral            )  CALL exchange_horiz( pt, nbgp )
     752                IF ( humidity                 )  CALL exchange_horiz( q, nbgp  )
     753                IF ( passive_scalar           )  CALL exchange_horiz( s, nbgp  )
     754                IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp  )
    755755             ENDIF
    756756!
  • palm/trunk/SOURCE/user_actions.f90

    r1874 r1960  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! New CASE statement for scalar tendency
    2222!
    2323! Former revisions:
     
    170170
    171171          CASE ( 'q-tendency' )
     172         
     173         
     174          CASE ( 's-tendency' )         
    172175
    173176
     
    228231
    229232          CASE ( 'q-tendency' )
     233         
     234         
     235          CASE ( 's-tendency' )
    230236
    231237
  • palm/trunk/SOURCE/user_spectra.f90

    r1834 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Consider additional spectra for passive scalar
    2222!
    2323! Former revisions:
     
    8989       SELECT CASE ( TRIM( data_output_sp(m) ) )
    9090         
    91           CASE ( 'u', 'v', 'w', 'pt', 'q' )
     91          CASE ( 'u', 'v', 'w', 'pt', 'q', 's' )
    9292!--          Not allowed here since these are the standard quantities used in
    9393!--          preprocess_spectra.
     
    108108       SELECT CASE ( TRIM( data_output_sp(m) ) )
    109109
    110           CASE ( 'u', 'v', 'w', 'pt', 'q' )
     110          CASE ( 'u', 'v', 'w', 'pt', 'q', 's' )
    111111!--          Not allowed here since these are the standard quantities used in
    112112!--          data_output_spectra.
  • palm/trunk/SOURCE/virtual_flight_mod.f90

    r1958 r1960  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Separate humidity and passive scalar.
     22! Bugfix concerning labeling of timeseries.
    2223!
    2324! Former revisions:
     
    326327       
    327328       INTEGER(iwp) ::  i               !< loop variable
     329       INTEGER(iwp) ::  id_pt           !< identifyer for labeling
     330       INTEGER(iwp) ::  id_q            !< identifyer for labeling
     331       INTEGER(iwp) ::  id_ql           !< identifyer for labeling
     332       INTEGER(iwp) ::  id_s            !< identifyer for labeling       
     333       INTEGER(iwp) ::  id_u = 1        !< identifyer for labeling 
     334       INTEGER(iwp) ::  id_v = 2        !< identifyer for labeling
     335       INTEGER(iwp) ::  id_w = 3        !< identifyer for labeling
    328336       INTEGER(iwp) ::  k               !< index variable
    329337       
     
    332340!--    Define output quanities, at least three variables are measured (u,v,w)
    333341       num_var_fl = 3
    334        IF ( .NOT. neutral )                      num_var_fl = num_var_fl + 1
    335        IF ( humidity .OR. passive_scalar )       num_var_fl = num_var_fl + 1
    336        IF ( cloud_physics .OR. cloud_droplets )  num_var_fl = num_var_fl + 1 
     342       IF ( .NOT. neutral                     )  THEN
     343          num_var_fl = num_var_fl + 1
     344          id_pt      = num_var_fl
     345       ENDIF
     346       IF ( humidity                          )  THEN
     347          num_var_fl = num_var_fl + 1
     348          id_q       = num_var_fl
     349       ENDIF
     350       IF ( cloud_physics .OR. cloud_droplets )  THEN
     351          num_var_fl = num_var_fl + 1 
     352          id_ql      = num_var_fl
     353       ENDIF
     354       IF ( passive_scalar                    )  THEN
     355          num_var_fl = num_var_fl + 1
     356          id_s       = num_var_fl
     357       ENDIF
    337358!
    338359!--    Write output strings for dimensions x, y, z
     
    364385          DO i=1, num_var_fl
    365386
    366              SELECT CASE ( i )
    367                      
    368                 CASE( 1 )
    369                    dofl_label(k) = TRIM( label_leg ) // '_u'
    370                    dofl_unit(k)  = 'm/s'
    371                    k             = k + 1
    372 
    373                 CASE( 2 )
    374                    dofl_label(k) = TRIM( label_leg ) // '_v'
    375                    dofl_unit(k)  = 'm/s'
    376                    k             = k + 1
    377 
    378                 CASE( 3 )
    379                    dofl_label(k) = TRIM( label_leg ) // '_w'
    380                    dofl_unit(k)  = 'm/s'
    381                    k             = k + 1
    382 
    383                 CASE( 4 )
    384                    dofl_label(k) = TRIM( label_leg ) // '_pt'
    385                    dofl_unit(k)  = 'K'
    386                    k             = k + 1
    387 
    388                 CASE( 5 )
    389                    dofl_label(k) = TRIM( label_leg ) // '_q'
    390                    dofl_unit(k)  = 'kg/kg'
    391                    k             = k + 1
    392 
    393                 CASE( 6 )
    394                    dofl_label(k) = TRIM( label_leg ) // '_ql'
    395                    dofl_unit(k)  = 'kg/kg'
    396                    k             = k + 1
    397 
    398              END SELECT
    399 
     387             IF ( i == id_u      )  THEN         
     388                dofl_label(k) = TRIM( label_leg ) // '_u'
     389                dofl_unit(k)  = 'm/s'
     390                k             = k + 1
     391             ELSEIF ( i == id_v  )  THEN       
     392                dofl_label(k) = TRIM( label_leg ) // '_v'
     393                dofl_unit(k)  = 'm/s'
     394                k             = k + 1
     395             ELSEIF ( i == id_w  )  THEN         
     396                dofl_label(k) = TRIM( label_leg ) // '_w'
     397                dofl_unit(k)  = 'm/s'
     398                k             = k + 1
     399             ELSEIF ( i == id_pt )  THEN       
     400                dofl_label(k) = TRIM( label_leg ) // '_pt'
     401                dofl_unit(k)  = 'K'
     402                k             = k + 1
     403             ELSEIF ( i == id_q  )  THEN       
     404                dofl_label(k) = TRIM( label_leg ) // '_q'
     405                dofl_unit(k)  = 'kg/kg'
     406                k             = k + 1
     407             ELSEIF ( i == id_ql )  THEN       
     408                dofl_label(k) = TRIM( label_leg ) // '_ql'
     409                dofl_unit(k)  = 'kg/kg'
     410                k             = k + 1
     411             ELSEIF ( i == id_s  )  THEN                         
     412                dofl_label(k) = TRIM( label_leg ) // '_s'
     413                dofl_unit(k)  = 'kg/kg'
     414                k             = k + 1
     415             ENDIF
    400416          ENDDO
    401417         
     
    421437
    422438       USE arrays_3d,                                                          &
    423            ONLY:  ddzu, ddzw, pt, q, ql, u, v, w, zu, zw
     439           ONLY:  ddzu, ddzw, pt, q, ql, s, u, v, w, zu, zw
    424440
    425441       USE control_parameters,                                                 &
     
    596612                ENDIF
    597613!
    598 !--             Humidity or passive scalar
    599                 IF ( humidity .OR. passive_scalar )  THEN
     614!--             Humidity
     615                IF ( humidity )  THEN
    600616                   CALL interpolate_xyz( q, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
    601617                   var_index = var_index + 1
     
    605621                IF ( cloud_physics .OR. cloud_droplets )  THEN
    606622                   CALL interpolate_xyz( ql, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
     623                   var_index = var_index + 1
     624                ENDIF
     625!
     626!--             Passive scalar
     627                IF ( passive_scalar )  THEN
     628                   CALL interpolate_xyz( s, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
    607629                   var_index = var_index + 1
    608630                ENDIF
  • palm/trunk/SOURCE/write_3d_binary.f90

    r1852 r1960  
    100100        ONLY:  e, kh, km, ol, p, pt, q, ql, qc, nr, nrs, nrsws, nrswst,        &
    101101               prr, precipitation_amount, qr,                                  &
    102                qrs, qrsws, qrswst, qs, qsws, qswst, sa, saswsb, saswst,        &
    103                rif_wall, shf, ts, tswst, u, u_m_l, u_m_n, u_m_r, u_m_s, us,    &
    104                usws, uswst, v, v_m_l, v_m_n, v_m_r, v_m_s, vpt, vsws, vswst,   &
    105                w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h, z0q
     102               qrs, qrsws, qrswst, qs, qsws, qswst, s, sa, ssws, sswst,        &
     103               saswsb, saswst, rif_wall, shf, ss, ts, tswst, u, u_m_l, u_m_n,  &
     104               u_m_r, u_m_s, us, usws, uswst, v, v_m_l, v_m_n, v_m_r, v_m_s,   &
     105               vpt, vsws, vswst, w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h, z0q
    106106       
    107107    USE averaging
     
    237237       WRITE ( 14 )  'pt_av               ';  WRITE ( 14 )  pt_av
    238238    ENDIF
    239     IF ( humidity  .OR. passive_scalar )  THEN
     239    IF ( humidity )  THEN
    240240       WRITE ( 14 )  'q                   ';  WRITE ( 14 )  q
    241241       IF ( ALLOCATED( q_av ) )  THEN
     
    277277       WRITE ( 14 )  'qswst               ';  WRITE ( 14 ) qswst
    278278    ENDIF
     279    IF ( passive_scalar )  THEN
     280       WRITE ( 14 )  's                   ';  WRITE ( 14 )  s
     281       IF ( ALLOCATED( s_av ) )  THEN
     282          WRITE ( 14 )  's_av                ';  WRITE ( 14 )  s_av
     283       ENDIF
     284       WRITE ( 14 )  'ss                  ';  WRITE ( 14 )  ss
     285       WRITE ( 14 )  'ssws                ';  WRITE ( 14 )  ssws
     286       IF ( ALLOCATED( ssws_av ) )  THEN
     287          WRITE ( 14 )  'ssws_av             ';  WRITE ( 14 )  ssws_av
     288       ENDIF
     289       WRITE ( 14 )  'sswst               ';  WRITE ( 14 ) sswst
     290    ENDIF   
    279291    IF ( land_surface )  THEN
    280292       IF ( ALLOCATED( qsws_eb_av ) )  THEN
  • palm/trunk/SOURCE/write_var_list.f90

    r1958 r1960  
    154154    USE arrays_3d,                                                             &
    155155        ONLY:  inflow_damping_factor, mean_inflow_profiles, pt_init,           &
    156                q_init, ref_state, sa_init, u_init, ug, v_init, vg
     156               q_init, ref_state, s_init, sa_init, u_init, ug, v_init, vg
    157157
    158158    USE control_parameters
     
    504504    WRITE ( 14 )  'run_coupled                   '
    505505    WRITE ( 14 )  run_coupled
     506    WRITE ( 14 )  's_init                        '
     507    WRITE ( 14 )  s_init
     508    WRITE ( 14 )  's_surface                     '
     509    WRITE ( 14 )  s_surface
     510    WRITE ( 14 )  's_surface_initial_change      '
     511    WRITE ( 14 )  s_surface_initial_change
     512    WRITE ( 14 )  's_vertical_gradient           '
     513    WRITE ( 14 )  s_vertical_gradient
     514    WRITE ( 14 )  's_vertical_gradient_level     '
     515    WRITE ( 14 )  s_vertical_gradient_level
     516    WRITE ( 14 )  's_vertical_gradient_level_ind '
     517    WRITE ( 14 )  s_vertical_gradient_level_ind
    506518    WRITE ( 14 )  'sa_init                       '
    507519    WRITE ( 14 )  sa_init
Note: See TracChangeset for help on using the changeset viewer.