Ignore:
Timestamp:
Dec 14, 2017 5:12:51 PM (4 years ago)
Author:
kanani
Message:

Merge of branch palm4u into trunk

Location:
palm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90

    r2342 r2696  
    1 !> @file lsf_nudging_mod.f90
    2 !------------------------------------------------------------------------------!
    3 ! This file is part of PALM.
     1!> @file large_scale_forcing_nudging_mod.f90
     2!------------------------------------------------------------------------------!
     3! This file is part of the PALM model system.
    44!
    55! PALM is free software: you can redistribute it and/or modify it under the
     
    2525! -----------------
    2626! $Id$
     27! Forcing with larger-scale models implemented (MS)
     28!
     29! 2342 2017-08-08 11:00:43Z boeske
    2730! fixed check if surface forcing data is available until end of simulation
    2831!
     
    4649
    4750    USE arrays_3d,                                                             &
    48         ONLY:  heatflux_input_conversion, pt, pt_init, q, q_init, tend,        &
    49                u, u_init, ug, v, v_init, vg, w_subs,                           &
    50                waterflux_input_conversion, zu                 
     51        ONLY:  dzw, e, heatflux_input_conversion, pt, pt_init, q, q_init, s,   &
     52               tend, u, u_init, ug, v, v_init, vg, w, w_subs,                  &
     53               waterflux_input_conversion, zu, zw                 
    5154
    5255    USE control_parameters,                                                    &
    53         ONLY:  bc_lr, bc_ns, bc_pt_b, bc_q_b, data_output_pr, dt_3d, end_time, & 
     56        ONLY:  bc_lr, bc_ns, bc_pt_b, bc_q_b, constant_diffusion,              &
     57               data_output_pr, dt_3d, end_time, forcing,                       &
     58               force_bound_l, force_bound_n, force_bound_r, force_bound_s,     & 
    5459               humidity, intermediate_timestep_count, ibc_pt_b, ibc_q_b,       &
    5560               large_scale_forcing, large_scale_subsidence, lsf_surf, lsf_vert,&
    56                lsf_exception, message_string, nudging, passive_scalar,         &
     61               lsf_exception, message_string, neutral, nudging, passive_scalar,&
    5762               pt_surface, ocean, q_surface, surface_pressure, topography,     &
    5863               use_subsidence_tendencies
    5964
     65    USE grid_variables
     66
     67    USE pegrid
    6068
    6169    USE indices,                                                               &
    62         ONLY:  ngp_sums_ls, nxl, nxr, nys, nyn, nzb, nz, nzt, wall_flags_0
     70        ONLY:  nbgp, ngp_sums_ls, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys,     &
     71               nysv, nysg, nyn, nyng, nzb, nz, nzt, wall_flags_0
    6372
    6473    USE kinds
     
    7079        ONLY:  hom, statistic_regions, sums_ls_l, weight_substep
    7180
     81    USE netcdf_data_input_mod,                                                 &
     82        ONLY:  force, netcdf_data_input_interpolate
    7283
    7384    INTEGER(iwp) ::  nlsf = 1000                       !< maximum number of profiles in LSF_DATA (large scale forcing)
    7485    INTEGER(iwp) ::  ntnudge = 1000                    !< maximum number of profiles in NUDGING_DATA (nudging)
     86
     87    REAL(wp) ::  d_area_t
    7588
    7689    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ptnudge     !< vertical profile of pot. temperature interpolated to vertical grid (nudging)
     
    107120           lsf_nudging_check_parameters, nudge_init,                           &
    108121           lsf_nudging_check_data_output_pr, lsf_nudging_header,               &
    109            calc_tnudge, nudge, nudge_ref
     122           calc_tnudge, nudge, nudge_ref, forcing_bc_mass_conservation,        &
     123           forcing_bc
    110124!
    111125!-- Public variables
    112126    PUBLIC qsws_surf, shf_surf, td_lsa_lpt, td_lsa_q, td_sub_lpt,              &
    113            td_sub_q, time_vert
     127           td_sub_q, time_vert, force
     128
    114129
    115130    INTERFACE ls_advec
     
    125140 CONTAINS
    126141
     142
     143!------------------------------------------------------------------------------!
     144! Description:
     145! ------------
     146!> @todo Missing subroutine description.
     147!------------------------------------------------------------------------------!
     148    SUBROUTINE forcing_bc_mass_conservation
     149
     150       USE control_parameters,                                                 &
     151           ONLY:  volume_flow
     152
     153       IMPLICIT NONE
     154
     155       INTEGER(iwp) ::  i !<
     156       INTEGER(iwp) ::  j !<
     157       INTEGER(iwp) ::  k !<
     158
     159       REAL(wp) ::  w_correct !<
     160       REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !<
     161
     162       volume_flow   = 0.0_wp
     163       volume_flow_l = 0.0_wp
     164
     165       d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
     166
     167       IF ( force_bound_l )  THEN
     168          i = nxl
     169          DO  j = nys, nyn
     170             DO  k = nzb+1, nzt
     171                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy   &
     172                                     * MERGE( 1.0_wp, 0.0_wp,                  &
     173                                              BTEST( wall_flags_0(k,j,i), 1 ) )
     174             ENDDO
     175          ENDDO
     176       ENDIF
     177       IF ( force_bound_r )  THEN
     178          i = nxr+1
     179          DO  j = nys, nyn
     180             DO  k = nzb+1, nzt
     181                volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy   &
     182                                     * MERGE( 1.0_wp, 0.0_wp,                  &
     183                                              BTEST( wall_flags_0(k,j,i), 1 ) )
     184             ENDDO
     185          ENDDO
     186       ENDIF
     187       IF ( force_bound_s )  THEN
     188          j = nys
     189          DO  i = nxl, nxr
     190             DO  k = nzb+1, nzt
     191                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx   &
     192                                     * MERGE( 1.0_wp, 0.0_wp,                  &
     193                                              BTEST( wall_flags_0(k,j,i), 2 ) )
     194             ENDDO
     195          ENDDO
     196       ENDIF
     197       IF ( force_bound_n )  THEN
     198          j = nyn+1
     199          DO  i = nxl, nxr
     200             DO  k = nzb+1, nzt
     201                volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx   &
     202                                     * MERGE( 1.0_wp, 0.0_wp,                  &
     203                                              BTEST( wall_flags_0(k,j,i), 2 ) )
     204             ENDDO
     205          ENDDO
     206       ENDIF
     207!
     208!--       Top boundary
     209       k = nzt
     210       DO  i = nxl, nxr
     211          DO  j = nys, nyn
     212             volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
     213          ENDDO
     214       ENDDO
     215
     216#if defined( __parallel )
     217       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     218       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM,      &
     219                           comm2d, ierr )
     220#else
     221       volume_flow = volume_flow_l
     222#endif
     223
     224       w_correct = SUM( volume_flow ) * d_area_t
     225
     226       DO  i = nxl, nxr
     227          DO  j = nys, nyn
     228             DO  k = nzt, nzt + 1
     229                w(k,j,i) = w(k,j,i) + w_correct
     230             ENDDO
     231          ENDDO
     232       ENDDO
     233
     234    END SUBROUTINE forcing_bc_mass_conservation
     235
     236
     237!------------------------------------------------------------------------------!
     238! Description:
     239! ------------
     240!> @todo Missing subroutine description.
     241!------------------------------------------------------------------------------!
     242    SUBROUTINE forcing_bc
     243
     244       USE control_parameters,                                                 &
     245           ONLY:  force_bound_l, force_bound_n, force_bound_r, force_bound_s,  &
     246                  humidity, neutral, passive_scalar, simulated_time
     247
     248       USE netcdf_data_input_mod,                                              &
     249           ONLY:  force                             
     250
     251       IMPLICIT NONE
     252
     253       INTEGER(iwp) ::  i !< running index x-direction
     254       INTEGER(iwp) ::  j !< running index y-direction
     255       INTEGER(iwp) ::  k !< running index z-direction
     256       INTEGER(iwp) ::  t !< running index for time levels
     257
     258       REAL(wp) ::  ddt_lsf !< inverse value of time resolution of forcing data
     259       REAL(wp) ::  t_ref   !< time past since last reference step
     260
     261!
     262!--    If required, interpolate and/or extrapolate data vertically. This is
     263!--    required as Inifor outputs only equidistant vertical data.
     264       IF ( ANY( zu(1:nzt+1) /= force%zu_atmos(1:force%nzu) ) )  THEN
     265          IF ( .NOT. force%interpolated )  THEN
     266
     267             DO  t = 0, 1
     268                IF ( force_bound_l )  THEN
     269                   CALL netcdf_data_input_interpolate( force%u_left(t,:,:),    &
     270                                                       zu(nzb+1:nzt+1),        &
     271                                                       force%zu_atmos )
     272                   CALL netcdf_data_input_interpolate( force%v_left(t,:,:),    &
     273                                                       zu(nzb+1:nzt+1),        &
     274                                                       force%zu_atmos )
     275                   CALL netcdf_data_input_interpolate( force%w_left(t,:,:),    &
     276                                                       zw(nzb+1:nzt+1),        &
     277                                                       force%zw_atmos )
     278                   IF ( .NOT. neutral )                                        &
     279                      CALL netcdf_data_input_interpolate( force%pt_left(t,:,:),&
     280                                                          zu(nzb+1:nzt+1),     &
     281                                                          force%zu_atmos )
     282                   IF ( humidity )                                             &
     283                      CALL netcdf_data_input_interpolate( force%q_left(t,:,:), &
     284                                                          zu(nzb+1:nzt+1),     &
     285                                                          force%zu_atmos )
     286                ENDIF
     287                IF ( force_bound_r )  THEN
     288                   CALL netcdf_data_input_interpolate( force%u_right(t,:,:),   &
     289                                                       zu(nzb+1:nzt+1),        &
     290                                                       force%zu_atmos )
     291                   CALL netcdf_data_input_interpolate( force%v_right(t,:,:),   &
     292                                                       zu(nzb+1:nzt+1),        &
     293                                                       force%zu_atmos )
     294                   CALL netcdf_data_input_interpolate( force%w_right(t,:,:),   &
     295                                                       zw(nzb+1:nzt+1),        &
     296                                                       force%zw_atmos )
     297                   IF ( .NOT. neutral )                                        &
     298                      CALL netcdf_data_input_interpolate( force%pt_right(t,:,:),&
     299                                                          zu(nzb+1:nzt+1),     &
     300                                                          force%zu_atmos )
     301                   IF ( humidity )                                             &
     302                      CALL netcdf_data_input_interpolate( force%q_right(t,:,:),&
     303                                                          zu(nzb+1:nzt+1),     &
     304                                                          force%zu_atmos )
     305                ENDIF
     306                IF ( force_bound_n )  THEN
     307                   CALL netcdf_data_input_interpolate( force%u_north(t,:,:),   &
     308                                                       zu(nzb+1:nzt+1),        &
     309                                                       force%zu_atmos )
     310                   CALL netcdf_data_input_interpolate( force%v_north(t,:,:),   &
     311                                                       zu(nzb+1:nzt+1),        &
     312                                                       force%zu_atmos )
     313                   CALL netcdf_data_input_interpolate( force%w_north(t,:,:),   &
     314                                                       zw(nzb+1:nzt+1),        &
     315                                                       force%zw_atmos )
     316                   IF ( .NOT. neutral )                                        &
     317                      CALL netcdf_data_input_interpolate( force%pt_north(t,:,:),&
     318                                                          zu(nzb+1:nzt+1),     &
     319                                                          force%zu_atmos )
     320                   IF ( humidity )                                             &
     321                      CALL netcdf_data_input_interpolate( force%q_north(t,:,:),&
     322                                                          zu(nzb+1:nzt+1),     &
     323                                                          force%zu_atmos )
     324                ENDIF
     325                IF ( force_bound_s )  THEN
     326                   CALL netcdf_data_input_interpolate( force%u_south(t,:,:),   &
     327                                                       zu(nzb+1:nzt+1),        &
     328                                                       force%zu_atmos )
     329                   CALL netcdf_data_input_interpolate( force%v_south(t,:,:),   &
     330                                                       zu(nzb+1:nzt+1),        &
     331                                                       force%zu_atmos )
     332                   CALL netcdf_data_input_interpolate( force%w_south(t,:,:),   &
     333                                                       zw(nzb+1:nzt+1),        &
     334                                                       force%zw_atmos )
     335                   IF ( .NOT. neutral )                                        &
     336                      CALL netcdf_data_input_interpolate( force%pt_south(t,:,:),&
     337                                                          zu(nzb+1:nzt+1),     &
     338                                                          force%zu_atmos )
     339                   IF ( humidity )                                             &
     340                      CALL netcdf_data_input_interpolate( force%q_south(t,:,:),&
     341                                                          zu(nzb+1:nzt+1),     &
     342                                                          force%zu_atmos )
     343                ENDIF
     344             ENDDO
     345!
     346!--          Note, no interpolation of top boundary. Just use initial value.
     347!--          No physical meaningful extrapolation possible if only one layer is
     348!--          given.
     349
     350             force%interpolated = .TRUE.
     351          ENDIF
     352       ENDIF
     353     
     354!
     355!--    Calculate time interval of forcing data       
     356       ddt_lsf = 1.0_wp / ( force%time(force%tind_p) - force%time(force%tind) )
     357!
     358!--    Calculate reziproke time past since last reference step. Please note,
     359!--    as simulated time is still not updated, the actual time here is
     360!--    simulated time + dt_3d
     361       t_ref = simulated_time + dt_3d - force%time(force%tind)
     362
     363       IF ( force_bound_l )  THEN
     364
     365          DO  j = nys, nyn
     366             DO  k = nzb+1, nzt+1
     367                u(k,j,nxlg:nxl) = force%u_left(0,k,j) + ddt_lsf * t_ref *      &
     368                         ( force%u_left(1,k,j) - force%u_left(0,k,j) ) *       &
     369                           MERGE( 1.0_wp, 0.0_wp,                              &
     370                                  BTEST( wall_flags_0(k,j,nxlg:nxl), 1 ) )
     371             ENDDO
     372          ENDDO
     373
     374          DO  j = nys, nyn
     375             DO  k = nzb+1, nzt
     376                w(k,j,nxlg:nxl-1) = force%w_left(0,k,j) + ddt_lsf * t_ref *    &
     377                         ( force%w_left(1,k,j) - force%w_left(0,k,j) ) *       &
     378                           MERGE( 1.0_wp, 0.0_wp,                              &
     379                                  BTEST( wall_flags_0(k,j,nxlg:nxl-1), 3 ) )
     380             ENDDO
     381          ENDDO
     382
     383          DO  j = nysv, nyn
     384             DO  k = nzb+1, nzt+1
     385                v(k,j,nxlg:nxl-1) = force%v_left(0,k,j) + ddt_lsf * t_ref *    &
     386                         ( force%v_left(1,k,j) - force%v_left(0,k,j) ) *       &
     387                           MERGE( 1.0_wp, 0.0_wp,                              &
     388                                  BTEST( wall_flags_0(k,j,nxlg:nxl-1), 2 ) )
     389             ENDDO
     390          ENDDO
     391
     392          IF ( .NOT. neutral )  THEN
     393             DO  j = nys, nyn
     394                DO  k = nzb+1, nzt+1
     395                   pt(k,j,nxlg:nxl-1) = force%pt_left(0,k,j) + ddt_lsf *       &
     396                                                                  t_ref   *    &
     397                       ( force%pt_left(1,k,j) - force%pt_left(0,k,j) )
     398 
     399                ENDDO
     400             ENDDO
     401          ENDIF
     402
     403          IF ( humidity )  THEN
     404             DO  j = nys, nyn
     405                DO  k = nzb+1, nzt+1
     406                   q(k,j,nxlg:nxl-1) = force%q_left(0,k,j) + ddt_lsf   *       &
     407                                                                  t_ref   *    &
     408                       ( force%q_left(1,k,j) - force%q_left(0,k,j) )
     409 
     410                ENDDO
     411             ENDDO
     412          ENDIF
     413
     414       ENDIF
     415
     416       IF ( force_bound_r )  THEN
     417
     418          DO  j = nys, nyn
     419             DO  k = nzb+1, nzt+1
     420                u(k,j,nxr+1:nxrg) = force%u_right(0,k,j) + ddt_lsf * t_ref *   &
     421                        ( force%u_right(1,k,j) - force%u_right(0,k,j) ) *      &
     422                           MERGE( 1.0_wp, 0.0_wp,                              &
     423                                  BTEST( wall_flags_0(k,j,nxr+1:nxrg), 1 ) )
     424
     425             ENDDO
     426          ENDDO
     427          DO  j = nys, nyn
     428             DO  k = nzb+1, nzt
     429                w(k,j,nxr+1:nxrg) = force%w_right(0,k,j) + ddt_lsf * t_ref *   &
     430                        ( force%w_right(1,k,j) - force%w_right(0,k,j) ) *      &
     431                           MERGE( 1.0_wp, 0.0_wp,                              &
     432                                  BTEST( wall_flags_0(k,j,nxr+1:nxrg), 3 ) )
     433             ENDDO
     434          ENDDO
     435
     436          DO  j = nysv, nyn
     437             DO  k = nzb+1, nzt+1
     438                v(k,j,nxr+1:nxrg) = force%v_right(0,k,j) + ddt_lsf * t_ref *   &
     439                        ( force%v_right(1,k,j) - force%v_right(0,k,j) ) *      &
     440                           MERGE( 1.0_wp, 0.0_wp,                              &
     441                                  BTEST( wall_flags_0(k,j,nxr+1:nxrg), 2 ) )
     442             ENDDO
     443          ENDDO
     444
     445          IF ( .NOT. neutral )  THEN
     446             DO  j = nys, nyn
     447                DO  k = nzb+1, nzt+1
     448                   pt(k,j,nxr+1:nxrg) = force%pt_right(0,k,j) + ddt_lsf *      &
     449                                                                  t_ref   *    &
     450                     ( force%pt_right(1,k,j) - force%pt_right(0,k,j) )
     451 
     452                ENDDO
     453             ENDDO
     454          ENDIF
     455
     456          IF ( humidity )  THEN
     457             DO  j = nys, nyn
     458                DO  k = nzb+1, nzt+1
     459                   q(k,j,nxr+1:nxrg) = force%q_right(0,k,j) + ddt_lsf   *      &
     460                                                                  t_ref   *    &
     461                       ( force%q_right(1,k,j) - force%q_right(0,k,j) )
     462 
     463                ENDDO
     464             ENDDO
     465          ENDIF
     466
     467       ENDIF
     468
     469       IF ( force_bound_s )  THEN
     470
     471          DO  i = nxl, nxr
     472             DO  k = nzb+1, nzt+1
     473                v(k,nysg:nys,i)   = force%v_south(0,k,i) + ddt_lsf * t_ref *   &
     474                        ( force%v_south(1,k,i) - force%v_south(0,k,i) ) *      &
     475                           MERGE( 1.0_wp, 0.0_wp,                              &
     476                                  BTEST( wall_flags_0(k,nysg:nys,i), 2 ) )
     477             ENDDO
     478          ENDDO
     479
     480          DO  i = nxl, nxr
     481             DO  k = nzb+1, nzt
     482                w(k,nysg:nys-1,i) = force%w_south(0,k,i) + ddt_lsf * t_ref *   &
     483                        ( force%w_south(1,k,i) - force%w_south(0,k,i) ) *      &
     484                           MERGE( 1.0_wp, 0.0_wp,                              &
     485                                  BTEST( wall_flags_0(k,nysg:nys-1,i), 3 ) )
     486             ENDDO
     487          ENDDO
     488
     489          DO  i = nxlu, nxr
     490             DO  k = nzb+1, nzt+1
     491                u(k,nysg:nys-1,i) = force%u_south(0,k,i) + ddt_lsf * t_ref *   &
     492                        ( force%u_south(1,k,i) - force%u_south(0,k,i) ) *      &
     493                           MERGE( 1.0_wp, 0.0_wp,                              &
     494                                  BTEST( wall_flags_0(k,nysg:nys-1,i), 1 ) )
     495             ENDDO
     496          ENDDO
     497
     498          IF ( .NOT. neutral )  THEN
     499             DO  i = nxl, nxr
     500                DO  k = nzb+1, nzt+1
     501                   pt(k,nysg:nys-1,i) = force%pt_south(0,k,i) + ddt_lsf *      &
     502                                                                  t_ref   *    &
     503                     ( force%pt_south(1,k,i) - force%pt_south(0,k,i) )
     504 
     505                ENDDO
     506             ENDDO
     507          ENDIF
     508
     509          IF ( humidity )  THEN
     510             DO  i = nxl, nxr
     511                DO  k = nzb+1, nzt+1
     512                   q(k,nysg:nys-1,i) = force%q_south(0,k,i) + ddt_lsf   *      &
     513                                                                  t_ref   *    &
     514                       ( force%q_south(1,k,i) - force%q_south(0,k,i) )
     515 
     516                ENDDO
     517             ENDDO
     518          ENDIF
     519
     520       ENDIF
     521
     522       IF ( force_bound_n )  THEN
     523
     524          DO  i = nxl, nxr
     525             DO  k = nzb+1, nzt+1
     526                v(k,nyn+1:nyng,i)   = force%v_north(0,k,i) + ddt_lsf * t_ref * &
     527                        ( force%v_north(1,k,i) - force%v_north(0,k,i) ) *      &
     528                           MERGE( 1.0_wp, 0.0_wp,                              &
     529                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 2 ) )
     530             ENDDO
     531          ENDDO
     532          DO  i = nxl, nxr
     533             DO  k = nzb+1, nzt
     534                w(k,nyn+1:nyng,i) = force%w_north(0,k,i) + ddt_lsf * t_ref *   &
     535                        ( force%w_north(1,k,i) - force%w_north(0,k,i) ) *      &
     536                           MERGE( 1.0_wp, 0.0_wp,                              &
     537                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 3 ) )
     538             ENDDO
     539          ENDDO
     540
     541          DO  i = nxlu, nxr
     542             DO  k = nzb+1, nzt+1
     543                u(k,nyn+1:nyng,i) = force%u_north(0,k,i) + ddt_lsf * t_ref *   &
     544                        ( force%u_north(1,k,i) - force%u_north(0,k,i) ) *      &
     545                           MERGE( 1.0_wp, 0.0_wp,                              &
     546                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 1 ) )
     547
     548             ENDDO
     549          ENDDO
     550
     551          IF ( .NOT. neutral )  THEN
     552             DO  i = nxl, nxr
     553                DO  k = nzb+1, nzt+1
     554                   pt(k,nyn+1:nyng,i) = force%pt_north(0,k,i) + ddt_lsf *      &
     555                                                                  t_ref   *    &
     556                     ( force%pt_north(1,k,i) - force%pt_north(0,k,i) )
     557 
     558                ENDDO
     559             ENDDO
     560          ENDIF
     561
     562          IF ( humidity )  THEN
     563             DO  i = nxl, nxr
     564                DO  k = nzb+1, nzt+1
     565                   q(k,nyn+1:nyng,i) = force%q_north(0,k,i) + ddt_lsf   *      &
     566                                                                  t_ref   *    &
     567                       ( force%q_north(1,k,i) - force%q_north(0,k,i) )
     568 
     569                ENDDO
     570             ENDDO
     571          ENDIF
     572
     573       ENDIF
     574!
     575!--    Top boundary.
     576!--    Please note, only map Inifor data on model top in case the numeric is
     577!--    identical to the Inifor grid. At the top boundary an extrapolation is
     578!--    not possible.
     579       IF ( ANY( zu(1:nzt+1) /= force%zu_atmos(1:force%nzu) ) )  THEN
     580          DO  i = nxlu, nxr
     581             DO  j = nys, nyn
     582                u(nzt+1,j,i) = force%u_top(0,j,i) + ddt_lsf * t_ref *          &
     583                           ( force%u_top(1,j,i) - force%u_top(0,j,i) ) *       &
     584                              MERGE( 1.0_wp, 0.0_wp,                           &
     585                                     BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
     586             ENDDO
     587          ENDDO
     588
     589          DO  i = nxl, nxr
     590             DO  j = nysv, nyn
     591                v(nzt+1,j,i) = force%v_top(0,j,i) + ddt_lsf * t_ref *          &
     592                           ( force%v_top(1,j,i) - force%v_top(0,j,i) ) *       &
     593                              MERGE( 1.0_wp, 0.0_wp,                           &
     594                                     BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
     595             ENDDO
     596          ENDDO
     597
     598          DO  i = nxl, nxr
     599             DO  j = nys, nyn
     600                w(nzt:nzt+1,j,i) = force%w_top(0,j,i) + ddt_lsf * t_ref *      &
     601                           ( force%w_top(1,j,i) - force%w_top(0,j,i) ) *       &
     602                              MERGE( 1.0_wp, 0.0_wp,                           &
     603                                     BTEST( wall_flags_0(nzt:nzt+1,j,i), 3 ) )
     604             ENDDO
     605          ENDDO
     606
     607
     608          IF ( .NOT. neutral )  THEN
     609             DO  i = nxl, nxr
     610                DO  j = nys, nyn
     611                   pt(nzt+1,j,i) = force%pt_top(0,j,i) + ddt_lsf * t_ref *     &
     612                           ( force%pt_top(1,j,i) - force%pt_top(0,j,i) )
     613                ENDDO
     614             ENDDO
     615          ENDIF
     616
     617          IF ( humidity )  THEN
     618             DO  i = nxl, nxr
     619                DO  j = nys, nyn
     620                   q(nzt+1,j,i) = force%q_top(0,j,i) + ddt_lsf * t_ref *       &
     621                           ( force%q_top(1,j,i) - force%q_top(0,j,i) )
     622                ENDDO
     623             ENDDO
     624          ENDIF
     625       ENDIF
     626!
     627!--    Moreover, set Neumann boundary condition for subgrid-scale TKE and
     628!--    passive scalar
     629       IF ( .NOT. constant_diffusion )  THEN
     630          IF (  force_bound_l )  e(:,:,nxl-1) = e(:,:,nxl)
     631          IF (  force_bound_r )  e(:,:,nxr+1) = e(:,:,nxr)
     632          IF (  force_bound_s )  e(:,nys-1,:) = e(:,nys,:)
     633          IF (  force_bound_n )  e(:,nyn+1,:) = e(:,nyn,:)
     634          e(nzt+1,:,:) = e(nzt,:,:)
     635       ENDIF
     636       IF ( passive_scalar )  THEN
     637          IF (  force_bound_l )  s(:,:,nxl-1) = s(:,:,nxl)
     638          IF (  force_bound_r )  s(:,:,nxr+1) = s(:,:,nxr)
     639          IF (  force_bound_s )  s(:,nys-1,:) = s(:,nys,:)
     640          IF (  force_bound_n )  s(:,nyn+1,:) = s(:,nyn,:)
     641       ENDIF
     642
     643
     644       CALL exchange_horiz( u, nbgp )
     645       CALL exchange_horiz( v, nbgp )
     646       CALL exchange_horiz( w, nbgp )
     647       IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
     648       IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
     649
     650!
     651!--    Set surface pressure. Please note, time-dependent surface
     652!--    pressure would require changes in anelastic approximation and
     653!--    treatment of fluxes.
     654!--    For the moment, comment this out!
     655!      surface_pressure = force%surface_pressure(force%tind) +                 &
     656!                                                      ddt_lsf * t_ref *       &
     657!                                    ( force%surface_pressure(force%tind_p)    &
     658!                                    - force%surface_pressure(force%tind) )
     659
     660    END SUBROUTINE forcing_bc
    127661
    128662!------------------------------------------------------------------------------!
     
    403937    SUBROUTINE lsf_init
    404938
     939       USE netcdf_data_input_mod,                                              &
     940           ONLY:  netcdf_data_input_lsf 
     941
    405942       IMPLICIT NONE
    406943
     
    411948       INTEGER(iwp) ::  finput = 90   !<
    412949       INTEGER(iwp) ::  k             !<
    413        INTEGER(iwp) ::  nt             !<
     950       INTEGER(iwp) ::  nt            !<
     951       INTEGER(iwp) ::  t             !< running index for time levels
    414952
    415953       REAL(wp) ::  fac               !<
     
    432970       REAL(wp) ::  r_dummy           !<
    433971
    434        ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),              &
    435                  qsws_surf(0:nlsf), shf_surf(0:nlsf),                          &
    436                  td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf),     &
    437                  td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf),     &
    438                  time_vert(0:nlsf), time_surf(0:nlsf),                         &
    439                  ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),         &
    440                  wsubs_vert(nzb:nzt+1,0:nlsf) )
    441 
    442        p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
    443        shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
    444        td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
    445        time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
    446        wsubs_vert = 0.0_wp
    447 
    448 !
    449 !--    Array for storing large scale forcing and nudging tendencies at each
    450 !--    timestep for data output
    451        ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
    452        sums_ls_l = 0.0_wp
    453 
    454        ngp_sums_ls = (nz+2)*6
    455 
    456        OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
    457               FORM='FORMATTED', IOSTAT=ierrn )
    458 
    459        IF ( ierrn /= 0 )  THEN
    460           message_string = 'file LSF_DATA does not exist'
    461           CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
    462        ENDIF
    463 
    464        ierrn = 0
    465 !
    466 !--    First three lines of LSF_DATA contain header
    467        READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
    468        READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
    469        READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
    470 
    471        IF ( ierrn /= 0 )  THEN
    472           message_string = 'errors in file LSF_DATA'
    473           CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
    474        ENDIF
    475 
    476 !
    477 !--    Surface values are read in
    478        nt     = 0
    479        ierrn = 0
    480 
    481        DO WHILE ( time_surf(nt) < end_time )
    482           nt = nt + 1
    483           READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),      &
    484                                              qsws_surf(nt), pt_surf(nt),       &
    485                                              q_surf(nt), p_surf(nt)
     972       IF ( forcing )  THEN
     973!
     974!--       Allocate arrays for reading boundary values. Arrays will incorporate 2
     975!--       time levels in order to interpolate in between.
     976          IF ( force_bound_l )  THEN
     977             ALLOCATE( force%u_left(0:1,nzb+1:nzt+1,nys:nyn)  )
     978             ALLOCATE( force%v_left(0:1,nzb+1:nzt+1,nysv:nyn) )
     979             ALLOCATE( force%w_left(0:1,nzb+1:nzt,nys:nyn)    )
     980             IF ( humidity )      ALLOCATE( force%q_left(0:1,nzb+1:nzt+1,nys:nyn)  )
     981             IF ( .NOT. neutral ) ALLOCATE( force%pt_left(0:1,nzb+1:nzt+1,nys:nyn) )
     982          ENDIF
     983          IF ( force_bound_r )  THEN
     984             ALLOCATE( force%u_right(0:1,nzb+1:nzt+1,nys:nyn)  )
     985             ALLOCATE( force%v_right(0:1,nzb+1:nzt+1,nysv:nyn) )
     986             ALLOCATE( force%w_right(0:1,nzb+1:nzt,nys:nyn)    )
     987             IF ( humidity )      ALLOCATE( force%q_right(0:1,nzb+1:nzt+1,nys:nyn)  )
     988             IF ( .NOT. neutral ) ALLOCATE( force%pt_right(0:1,nzb+1:nzt+1,nys:nyn) )
     989          ENDIF
     990          IF ( force_bound_n )  THEN
     991             ALLOCATE( force%u_north(0:1,nzb+1:nzt+1,nxlu:nxr) )
     992             ALLOCATE( force%v_north(0:1,nzb+1:nzt+1,nxl:nxr)  )
     993             ALLOCATE( force%w_north(0:1,nzb+1:nzt,nxl:nxr)    )
     994             IF ( humidity )      ALLOCATE( force%q_north(0:1,nzb+1:nzt+1,nxl:nxr)  )
     995             IF ( .NOT. neutral ) ALLOCATE( force%pt_north(0:1,nzb+1:nzt+1,nxl:nxr) )
     996          ENDIF
     997          IF ( force_bound_s )  THEN
     998             ALLOCATE( force%u_south(0:1,nzb+1:nzt+1,nxlu:nxr) )
     999             ALLOCATE( force%v_south(0:1,nzb+1:nzt+1,nxl:nxr)  )
     1000             ALLOCATE( force%w_south(0:1,nzb+1:nzt,nxl:nxr)    )
     1001             IF ( humidity )      ALLOCATE( force%q_south(0:1,nzb+1:nzt+1,nxl:nxr)  )
     1002             IF ( .NOT. neutral ) ALLOCATE( force%pt_south(0:1,nzb+1:nzt+1,nxl:nxr) )
     1003          ENDIF
     1004         
     1005          ALLOCATE( force%u_top(0:1,nys:nyn,nxlu:nxr) )
     1006          ALLOCATE( force%v_top(0:1,nysv:nyn,nxl:nxr) )
     1007          ALLOCATE( force%w_top(0:1,nys:nyn,nxl:nxr)  )
     1008          IF ( humidity )      ALLOCATE( force%q_top(0:1,nys:nyn,nxl:nxr)  )
     1009          IF ( .NOT. neutral ) ALLOCATE( force%pt_top(0:1,nys:nyn,nxl:nxr) )
     1010
     1011!
     1012!--       Initial call of input. Time array, initial 3D data of u, v, w,
     1013!--       potential temperature, as well as mixing ratio, will be read.
     1014!--       Moreover, data at lateral and top boundary will be read. 
     1015          CALL netcdf_data_input_lsf
     1016!
     1017!--       Please note, at the moment INIFOR assumes only an equidistant vertical
     1018!--       grid. In case of vertical grid stretching, input of inital 3D data
     1019!--       need to be inter- and/or extrapolated.
     1020!--       Therefore, check if zw grid on file is identical to numeric zw grid.   
     1021          IF ( ANY( zu(1:nzt+1) /= force%zu_atmos(1:force%nzu) ) )  THEN
     1022!
     1023!--          Also data at the boundaries need to be inter/extrapolated at both
     1024!--          time levels
     1025             DO  t = 0, 1
     1026                IF ( force_bound_l )  THEN
     1027                   CALL netcdf_data_input_interpolate( force%u_left(t,:,:),    &
     1028                                                       zu(1:nzt+1),            &
     1029                                                       force%zu_atmos )
     1030                   CALL netcdf_data_input_interpolate( force%v_left(t,:,:),    &
     1031                                                       zu(1:nzt+1),            &
     1032                                                       force%zu_atmos )
     1033                   CALL netcdf_data_input_interpolate( force%w_left(t,:,:),    &
     1034                                                       zw(1:nzt+1),            &
     1035                                                       force%zw_atmos )
     1036                   IF ( .NOT. neutral )                                        &
     1037                      CALL netcdf_data_input_interpolate( force%pt_left(t,:,:),&
     1038                                                          zu(1:nzt+1),         &
     1039                                                          force%zu_atmos )
     1040                   IF ( humidity )                                             &
     1041                      CALL netcdf_data_input_interpolate( force%q_left(t,:,:), &
     1042                                                          zu(1:nzt+1),         &
     1043                                                          force%zu_atmos )
     1044                ENDIF
     1045                IF ( force_bound_r )  THEN
     1046                   CALL netcdf_data_input_interpolate( force%u_right(t,:,:),   &
     1047                                                       zu(1:nzt+1),            &
     1048                                                       force%zu_atmos )
     1049                   CALL netcdf_data_input_interpolate( force%v_right(t,:,:),   &
     1050                                                       zu(1:nzt+1),            &
     1051                                                       force%zu_atmos )
     1052                   CALL netcdf_data_input_interpolate( force%w_right(t,:,:),   &
     1053                                                       zw(1:nzt+1),            &
     1054                                                       force%zw_atmos )
     1055                   IF ( .NOT. neutral )                                        &
     1056                      CALL netcdf_data_input_interpolate( force%pt_right(t,:,:),&
     1057                                                          zu(1:nzt+1),          &
     1058                                                          force%zu_atmos )
     1059                   IF ( humidity )                                             &
     1060                      CALL netcdf_data_input_interpolate( force%q_right(t,:,:),&
     1061                                                          zu(1:nzt+1),         &
     1062                                                          force%zu_atmos )
     1063                ENDIF
     1064                IF ( force_bound_n )  THEN
     1065                   CALL netcdf_data_input_interpolate( force%u_north(t,:,:),   &
     1066                                                       zu(1:nzt+1),            &
     1067                                                       force%zu_atmos )
     1068                   CALL netcdf_data_input_interpolate( force%v_north(t,:,:),   &
     1069                                                       zu(1:nzt+1),            &
     1070                                                       force%zu_atmos )
     1071                   CALL netcdf_data_input_interpolate( force%w_north(t,:,:),   &
     1072                                                       zw(1:nzt+1),            &
     1073                                                       force%zw_atmos )
     1074                   IF ( .NOT. neutral )                                        &
     1075                      CALL netcdf_data_input_interpolate( force%pt_north(t,:,:),&
     1076                                                          zu(1:nzt+1),          &
     1077                                                          force%zu_atmos )
     1078                   IF ( humidity )                                             &
     1079                      CALL netcdf_data_input_interpolate( force%q_north(t,:,:),&
     1080                                                          zu(1:nzt+1),         &
     1081                                                          force%zu_atmos )
     1082                ENDIF
     1083                IF ( force_bound_s )  THEN
     1084                   CALL netcdf_data_input_interpolate( force%u_south(t,:,:),   &
     1085                                                       zu(1:nzt+1),            &
     1086                                                       force%zu_atmos )
     1087                   CALL netcdf_data_input_interpolate( force%v_south(t,:,:),   &
     1088                                                       zu(1:nzt+1),            &
     1089                                                       force%zu_atmos )
     1090                   CALL netcdf_data_input_interpolate( force%w_south(t,:,:),   &
     1091                                                       zw(1:nzt+1),            &
     1092                                                       force%zw_atmos )
     1093                   IF ( .NOT. neutral )                                        &
     1094                      CALL netcdf_data_input_interpolate( force%pt_south(t,:,:),&
     1095                                                          zu(1:nzt+1),          &
     1096                                                          force%zu_atmos )
     1097                   IF ( humidity )                                             &
     1098                      CALL netcdf_data_input_interpolate( force%q_south(t,:,:),&
     1099                                                          zu(1:nzt+1),         &
     1100                                                          force%zu_atmos )
     1101                ENDIF
     1102             ENDDO
     1103          ENDIF
     1104
     1105!
     1106!--       Exchange ghost points
     1107          CALL exchange_horiz( u, nbgp )
     1108          CALL exchange_horiz( v, nbgp )
     1109          CALL exchange_horiz( w, nbgp )
     1110          IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
     1111          IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
     1112!
     1113!--       At lateral boundaries, set also initial boundary conditions
     1114          IF ( force_bound_l )  THEN
     1115             u(:,:,nxl)   = u(:,:,nxlu)
     1116             v(:,:,nxl-1) = v(:,:,nxl)
     1117             w(:,:,nxl-1) = w(:,:,nxl)
     1118             IF ( .NOT. neutral )  pt(:,:,nxl-1) = pt(:,:,nxl)
     1119             IF ( humidity      )  q(:,:,nxl-1)  = q(:,:,nxl)
     1120          ENDIF
     1121          IF ( force_bound_r )  THEN
     1122             u(:,:,nxr+1) = u(:,:,nxr)
     1123             v(:,:,nxr+1) = v(:,:,nxr)
     1124             w(:,:,nxr+1) = w(:,:,nxr)
     1125             IF ( .NOT. neutral )  pt(:,:,nxr+1) = pt(:,:,nxr)
     1126             IF ( humidity      )  q(:,:,nxr+1)  = q(:,:,nxr)
     1127          ENDIF
     1128          IF ( force_bound_s )  THEN
     1129             u(:,nys-1,:) = u(:,nys,:)
     1130             v(:,nys,:)   = v(:,nysv,:)
     1131             w(:,nys-1,:) = w(:,nys,:)
     1132             IF ( .NOT. neutral )  pt(:,nys-1,:) = pt(:,nys,:)
     1133             IF ( humidity      )  q(:,nys-1,:)  = q(:,nys,:)
     1134          ENDIF
     1135          IF ( force_bound_n )  THEN
     1136             u(:,nyn+1,:) = u(:,nyn,:)
     1137             v(:,nyn+1,:) = v(:,nyn,:)
     1138             w(:,nyn+1,:) = w(:,nyn,:)
     1139             IF ( .NOT. neutral )  pt(:,nyn+1,:) = pt(:,nyn,:)
     1140             IF ( humidity      )  q(:,nyn+1,:)  = q(:,nyn,:)
     1141          ENDIF
     1142
     1143!
     1144!--       After 3D data is initialized, ensure mass conservation
     1145          CALL forcing_bc_mass_conservation
     1146!
     1147!--       Initialize surface pressure. Please note, time-dependent surface
     1148!--       pressure would require changes in anelastic approximation and
     1149!--       treatment of fluxes.
     1150!--       For the moment, comment this out!
     1151!         surface_pressure = force%surface_pressure(0)
     1152
     1153       ELSE
     1154
     1155          ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),           &
     1156                    qsws_surf(0:nlsf), shf_surf(0:nlsf),                       &
     1157                    td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf),  &
     1158                    td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf),  &
     1159                    time_vert(0:nlsf), time_surf(0:nlsf),                      &
     1160                    ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),      &
     1161                    wsubs_vert(nzb:nzt+1,0:nlsf) )
     1162
     1163          p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
     1164          shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
     1165          td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
     1166          time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
     1167          wsubs_vert = 0.0_wp
     1168
     1169!
     1170!--       Array for storing large scale forcing and nudging tendencies at each
     1171!--       timestep for data output
     1172          ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
     1173          sums_ls_l = 0.0_wp
     1174
     1175          ngp_sums_ls = (nz+2)*6
     1176
     1177          OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
     1178                 FORM='FORMATTED', IOSTAT=ierrn )
    4861179
    4871180          IF ( ierrn /= 0 )  THEN
    488             WRITE ( message_string, * ) 'No time dependent surface variables ',&
    489                               'in&LSF_DATA for end of run found'
    490 
    491              CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
    492           ENDIF
    493        ENDDO
    494 
    495        IF ( time_surf(1) > end_time )  THEN
    496           WRITE ( message_string, * ) 'Time dependent surface variables in ',  &
    497                                       '&LSF_DATA set in after end of ' ,      &
    498                                       'simulation - lsf_surf is set to FALSE'
    499           CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
    500           lsf_surf = .FALSE.
    501        ENDIF
    502 
    503 !
    504 !--    Go to the end of the list with surface variables
    505        DO WHILE ( ierrn == 0 )
    506           READ ( finput, *, IOSTAT = ierrn ) r_dummy
    507        ENDDO
    508 
    509 !
    510 !--    Profiles of ug, vg and w_subs are read in (large scale forcing)
    511 
    512        nt = 0
    513        DO WHILE ( time_vert(nt) < end_time )
    514           nt = nt + 1
    515           hash = "#"
    516           ierrn = 1 ! not zero
    517 !
    518 !--       Search for the next line consisting of "# time",
    519 !--       from there onwards the profiles will be read
    520           DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
    521              READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
    522              IF ( ierrn < 0 )  THEN
    523                 WRITE( message_string, * ) 'No time dependent vertical profiles',&
    524                                  ' in&LSF_DATA for end of run found'
    525                 CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
    526              ENDIF
    527           ENDDO
    528 
    529           IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
    530 
    531           READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,  &
    532                                            lowwsubs_vert, low_td_lsa_lpt,      &
    533                                            low_td_lsa_q, low_td_sub_lpt,       &
    534                                            low_td_sub_q
     1181             message_string = 'file LSF_DATA does not exist'
     1182             CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
     1183          ENDIF
     1184
     1185          ierrn = 0
     1186!
     1187!--       First three lines of LSF_DATA contain header
     1188          READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
     1189          READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
     1190          READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
     1191
    5351192          IF ( ierrn /= 0 )  THEN
    5361193             message_string = 'errors in file LSF_DATA'
     
    5381195          ENDIF
    5391196
    540           READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,            &
    541                                            highvg_vert, highwsubs_vert,        &
    542                                            high_td_lsa_lpt, high_td_lsa_q,     &
    543                                            high_td_sub_lpt, high_td_sub_q
    544      
    545           IF ( ierrn /= 0 )  THEN
    546              message_string = 'errors in file LSF_DATA'
    547              CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
    548           ENDIF
    549 
    550 
    551           DO  k = nzb, nzt+1
    552              IF ( highheight < zu(k) )  THEN
    553                 lowheight      = highheight
    554                 lowug_vert     = highug_vert
    555                 lowvg_vert     = highvg_vert
    556                 lowwsubs_vert  = highwsubs_vert
    557                 low_td_lsa_lpt = high_td_lsa_lpt
    558                 low_td_lsa_q   = high_td_lsa_q
    559                 low_td_sub_lpt = high_td_sub_lpt
    560                 low_td_sub_q   = high_td_sub_q
    561 
    562                 ierrn = 0
    563                 READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,      &
    564                                                  highvg_vert, highwsubs_vert,  &
    565                                                  high_td_lsa_lpt,              &
    566                                                  high_td_lsa_q,                &
    567                                                  high_td_sub_lpt, high_td_sub_q
    568 
    569                 IF ( ierrn /= 0 )  THEN
    570                    WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ', &
    571                         'is higher than the maximum height in LSF_DATA which ',&
    572                         'is ', lowheight, 'm. Interpolation on PALM ',         &
    573                         'grid is not possible.'
    574                    CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
     1197!
     1198!--       Surface values are read in
     1199          nt     = 0
     1200          ierrn = 0
     1201
     1202          DO WHILE ( time_surf(nt) < end_time )
     1203             nt = nt + 1
     1204             READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),   &
     1205                                                qsws_surf(nt), pt_surf(nt),    &
     1206                                                q_surf(nt), p_surf(nt)
     1207
     1208             IF ( ierrn /= 0 )  THEN
     1209               WRITE ( message_string, * ) 'No time dependent surface variables ' //&
     1210                                 'in&LSF_DATA for end of run found'
     1211
     1212                CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
     1213             ENDIF
     1214          ENDDO
     1215
     1216          IF ( time_surf(1) > end_time )  THEN
     1217             WRITE ( message_string, * ) 'Time dependent surface variables in ' // &
     1218                                         '&LSF_DATA set in after end of ' ,        &
     1219                                         'simulation - lsf_surf is set to FALSE'
     1220             CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
     1221             lsf_surf = .FALSE.
     1222          ENDIF
     1223
     1224!
     1225!--       Go to the end of the list with surface variables
     1226          DO WHILE ( ierrn == 0 )
     1227             READ ( finput, *, IOSTAT = ierrn ) r_dummy
     1228          ENDDO
     1229
     1230!
     1231!--       Profiles of ug, vg and w_subs are read in (large scale forcing)
     1232
     1233          nt = 0
     1234          DO WHILE ( time_vert(nt) < end_time )
     1235             nt = nt + 1
     1236             hash = "#"
     1237             ierrn = 1 ! not zero
     1238!
     1239!--          Search for the next line consisting of "# time",
     1240!--          from there onwards the profiles will be read
     1241             DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
     1242                READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
     1243                IF ( ierrn < 0 )  THEN
     1244                   WRITE( message_string, * ) 'No time dependent vertical profiles',&
     1245                                    ' in&LSF_DATA for end of run found'
     1246                   CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
    5751247                ENDIF
    576 
     1248             ENDDO
     1249
     1250             IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
     1251
     1252             READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,&
     1253                                              lowwsubs_vert, low_td_lsa_lpt,    &
     1254                                              low_td_lsa_q, low_td_sub_lpt,     &
     1255                                              low_td_sub_q
     1256             IF ( ierrn /= 0 )  THEN
     1257                message_string = 'errors in file LSF_DATA'
     1258                CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
    5771259             ENDIF
    5781260
    579 !
    580 !--          Interpolation of prescribed profiles in space
    581              fac = (highheight-zu(k))/(highheight - lowheight)
    582 
    583              ug_vert(k,nt)    = fac * lowug_vert                               &
    584                                 + ( 1.0_wp - fac ) * highug_vert
    585              vg_vert(k,nt)    = fac * lowvg_vert                               &
    586                                 + ( 1.0_wp - fac ) * highvg_vert
    587              wsubs_vert(k,nt) = fac * lowwsubs_vert                            &
    588                                 + ( 1.0_wp - fac ) * highwsubs_vert
    589 
    590              td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt                           &
    591                                 + ( 1.0_wp - fac ) * high_td_lsa_lpt
    592              td_lsa_q(k,nt)   = fac * low_td_lsa_q                             &
    593                                 + ( 1.0_wp - fac ) * high_td_lsa_q
    594              td_sub_lpt(k,nt) = fac * low_td_sub_lpt                           &
    595                                 + ( 1.0_wp - fac ) * high_td_sub_lpt
    596              td_sub_q(k,nt)   = fac * low_td_sub_q                             &
    597                                 + ( 1.0_wp - fac ) * high_td_sub_q
    598 
    599           ENDDO
    600 
    601        ENDDO
    602 
    603 !
    604 !--    Large scale vertical velocity has to be zero at the surface
    605        wsubs_vert(nzb,:) = 0.0_wp
    606  
    607        IF ( time_vert(1) > end_time )  THEN
    608           WRITE ( message_string, * ) 'Time dependent large scale profile ',   &
    609                              'forcing from&LSF_DATA sets in after end of ' ,   &
    610                              'simulation - lsf_vert is set to FALSE'
    611           CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
    612           lsf_vert = .FALSE.
    613        ENDIF
    614 
    615        CLOSE( finput )
    616 
     1261             READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,         &
     1262                                              highvg_vert, highwsubs_vert,     &
     1263                                              high_td_lsa_lpt, high_td_lsa_q,  &
     1264                                              high_td_sub_lpt, high_td_sub_q
     1265         
     1266             IF ( ierrn /= 0 )  THEN
     1267                message_string = 'errors in file LSF_DATA'
     1268                CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
     1269             ENDIF
     1270
     1271
     1272             DO  k = nzb, nzt+1
     1273                IF ( highheight < zu(k) )  THEN
     1274                   lowheight      = highheight
     1275                   lowug_vert     = highug_vert
     1276                   lowvg_vert     = highvg_vert
     1277                   lowwsubs_vert  = highwsubs_vert
     1278                   low_td_lsa_lpt = high_td_lsa_lpt
     1279                   low_td_lsa_q   = high_td_lsa_q
     1280                   low_td_sub_lpt = high_td_sub_lpt
     1281                   low_td_sub_q   = high_td_sub_q
     1282
     1283                   ierrn = 0
     1284                   READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,    &
     1285                                                    highvg_vert, highwsubs_vert,&
     1286                                                    high_td_lsa_lpt,            &
     1287                                                    high_td_lsa_q,              &
     1288                                                    high_td_sub_lpt, high_td_sub_q
     1289
     1290                   IF ( ierrn /= 0 )  THEN
     1291                      WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ', &
     1292                           'is higher than the maximum height in LSF_DATA which ',&
     1293                           'is ', lowheight, 'm. Interpolation on PALM ',         &
     1294                           'grid is not possible.'
     1295                      CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
     1296                   ENDIF
     1297
     1298                ENDIF
     1299
     1300!
     1301!--             Interpolation of prescribed profiles in space
     1302                fac = (highheight-zu(k))/(highheight - lowheight)
     1303
     1304                ug_vert(k,nt)    = fac * lowug_vert                            &
     1305                                   + ( 1.0_wp - fac ) * highug_vert
     1306                vg_vert(k,nt)    = fac * lowvg_vert                            &
     1307                                   + ( 1.0_wp - fac ) * highvg_vert
     1308                wsubs_vert(k,nt) = fac * lowwsubs_vert                         &
     1309                                   + ( 1.0_wp - fac ) * highwsubs_vert
     1310
     1311                td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt                        &
     1312                                   + ( 1.0_wp - fac ) * high_td_lsa_lpt
     1313                td_lsa_q(k,nt)   = fac * low_td_lsa_q                          &
     1314                                   + ( 1.0_wp - fac ) * high_td_lsa_q
     1315                td_sub_lpt(k,nt) = fac * low_td_sub_lpt                        &
     1316                                   + ( 1.0_wp - fac ) * high_td_sub_lpt
     1317                td_sub_q(k,nt)   = fac * low_td_sub_q                          &
     1318                                   + ( 1.0_wp - fac ) * high_td_sub_q
     1319
     1320             ENDDO
     1321
     1322          ENDDO
     1323
     1324!
     1325!--       Large scale vertical velocity has to be zero at the surface
     1326          wsubs_vert(nzb,:) = 0.0_wp
     1327   
     1328          IF ( time_vert(1) > end_time )  THEN
     1329             WRITE ( message_string, * ) 'Time dependent large scale profile ',&
     1330                                'forcing from&LSF_DATA sets in after end of ' ,&
     1331                                'simulation - lsf_vert is set to FALSE'
     1332             CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
     1333             lsf_vert = .FALSE.
     1334          ENDIF
     1335
     1336          CLOSE( finput )
     1337
     1338       ENDIF
    6171339
    6181340    END SUBROUTINE lsf_init
     
    13962118    END SUBROUTINE nudge_ref
    13972119
     2120
    13982121 END MODULE lsf_nudging_mod
Note: See TracChangeset for help on using the changeset viewer.