Ignore:
Timestamp:
Oct 15, 2018 2:21:08 PM (3 years ago)
Author:
suehring
Message:

Offline nesting revised and separated from large_scale_forcing_mod; Time-dependent synthetic turbulence generator; bugfixes in USM/LSM radiation model and init_pegrid

File:
1 edited

Legend:

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

    r3294 r3347  
    2525! -----------------
    2626! $Id$
     27! Offline nesting parts are moved to an independent module nesting_offl_mod
     28!
     29! 3341 2018-10-15 10:31:27Z suehring
    2730! ocean renamed ocean_mode
    2831!
     
    8790
    8891    USE control_parameters,                                                    &
    89         ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
    90                bc_lr, bc_ns, bc_pt_b, bc_q_b, constant_diffusion,              &
     92        ONLY:  bc_lr, bc_ns, bc_pt_b, bc_q_b, constant_diffusion,              &
    9193               constant_heatflux, constant_waterflux,                          &
    9294               data_output_pr, dt_3d, end_time,                                &
     
    9496               ibc_pt_b, ibc_q_b,                                              &
    9597               large_scale_forcing, large_scale_subsidence, lsf_surf, lsf_vert,&
    96                lsf_exception, message_string, nesting_offline, neutral,        &
     98               lsf_exception, message_string, neutral,                         &
    9799               nudging, passive_scalar, pt_surface, ocean_mode, q_surface,     &
    98100               surface_heatflux, surface_pressure, surface_waterflux,          &
    99101               topography, use_subsidence_tendencies
     102               
     103    USE cpulog,                                                                &
     104        ONLY:  cpu_log, log_point
    100105
    101106    USE grid_variables
     
    107112    USE kinds
    108113
    109     USE netcdf_data_input_mod,                                                 &
    110         ONLY:  nest_offl
    111 
    112114    USE pegrid
    113115
     
    120122    INTEGER(iwp) ::  nlsf = 1000                       !< maximum number of profiles in LSF_DATA (large scale forcing)
    121123    INTEGER(iwp) ::  ntnudge = 1000                    !< maximum number of profiles in NUDGING_DATA (nudging)
    122 
    123     REAL(wp) ::  d_area_t
    124124
    125125    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ptnudge     !< vertical profile of pot. temperature interpolated to vertical grid (nudging)
     
    156156           lsf_nudging_check_parameters, nudge_init,                           &
    157157           lsf_nudging_check_data_output_pr, lsf_nudging_header,               &
    158            lsf_nesting_offline, lsf_nesting_offline_mass_conservation,         &
    159158           nudge, nudge_ref
    160159           
     
    177176 CONTAINS
    178177
    179 
    180 !------------------------------------------------------------------------------!
    181 ! Description:
    182 ! ------------
    183 !> In this subroutine a constant mass within the model domain is guaranteed.
    184 !> Larger-scale models may be based on a compressible equation system, which is
    185 !> not consistent with PALMs incompressible equation system. In order to avoid
    186 !> a decrease or increase of mass during the simulation, non-divergent flow
    187 !> through the lateral and top boundaries is compensated by the vertical wind
    188 !> component at the top boundary.
    189 !------------------------------------------------------------------------------!
    190     SUBROUTINE lsf_nesting_offline_mass_conservation
    191 
    192        USE control_parameters,                                                 &
    193            ONLY:  volume_flow
    194 
    195        IMPLICIT NONE
    196 
    197        INTEGER(iwp) ::  i !< grid index in x-direction
    198        INTEGER(iwp) ::  j !< grid index in y-direction
    199        INTEGER(iwp) ::  k !< grid index in z-direction
    200 
    201        REAL(wp) ::  w_correct                       !< vertical velocity increment required to compensate non-divergent flow through the boundaries
    202        REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !< local volume flow
    203 
    204        volume_flow   = 0.0_wp
    205        volume_flow_l = 0.0_wp
    206 
    207        d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
    208 
    209        IF ( bc_dirichlet_l )  THEN
    210           i = nxl
    211           DO  j = nys, nyn
    212              DO  k = nzb+1, nzt
    213                 volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy   &
    214                                      * MERGE( 1.0_wp, 0.0_wp,                  &
    215                                               BTEST( wall_flags_0(k,j,i), 1 ) )
    216              ENDDO
    217           ENDDO
    218        ENDIF
    219        IF ( bc_dirichlet_r )  THEN
    220           i = nxr+1
    221           DO  j = nys, nyn
    222              DO  k = nzb+1, nzt
    223                 volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy   &
    224                                      * MERGE( 1.0_wp, 0.0_wp,                  &
    225                                               BTEST( wall_flags_0(k,j,i), 1 ) )
    226              ENDDO
    227           ENDDO
    228        ENDIF
    229        IF ( bc_dirichlet_s )  THEN
    230           j = nys
    231           DO  i = nxl, nxr
    232              DO  k = nzb+1, nzt
    233                 volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx   &
    234                                      * MERGE( 1.0_wp, 0.0_wp,                  &
    235                                               BTEST( wall_flags_0(k,j,i), 2 ) )
    236              ENDDO
    237           ENDDO
    238        ENDIF
    239        IF ( bc_dirichlet_n )  THEN
    240           j = nyn+1
    241           DO  i = nxl, nxr
    242              DO  k = nzb+1, nzt
    243                 volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx   &
    244                                      * MERGE( 1.0_wp, 0.0_wp,                  &
    245                                               BTEST( wall_flags_0(k,j,i), 2 ) )
    246              ENDDO
    247           ENDDO
    248        ENDIF
    249 !
    250 !--    Top boundary
    251        k = nzt
    252        DO  i = nxl, nxr
    253           DO  j = nys, nyn
    254              volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
    255           ENDDO
    256        ENDDO
    257 
    258 #if defined( __parallel )
    259        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    260        CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM,      &
    261                            comm2d, ierr )
    262 #else
    263        volume_flow = volume_flow_l
    264 #endif
    265 
    266        w_correct = SUM( volume_flow ) * d_area_t
    267 
    268        DO  i = nxl, nxr
    269           DO  j = nys, nyn
    270              DO  k = nzt, nzt + 1
    271                 w(k,j,i) = w(k,j,i) + w_correct
    272              ENDDO
    273           ENDDO
    274        ENDDO
    275 
    276     END SUBROUTINE lsf_nesting_offline_mass_conservation
    277 
    278 
    279 !------------------------------------------------------------------------------!
    280 ! Description:
    281 ! ------------
    282 !> Set the lateral and top boundary conditions in case the PALM domain is
    283 !> nested offline in a mesoscale model.
    284 !------------------------------------------------------------------------------!
    285     SUBROUTINE lsf_nesting_offline
    286 
    287        USE control_parameters,                                                 &
    288            ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,              &
    289                   bc_dirichlet_s, humidity, neutral, passive_scalar, rans_mode,&
    290                   rans_tke_e, time_since_reference_point                     
    291 
    292        IMPLICIT NONE
    293 
    294        INTEGER(iwp) ::  i !< running index x-direction
    295        INTEGER(iwp) ::  j !< running index y-direction
    296        INTEGER(iwp) ::  k !< running index z-direction
    297 
    298        REAL(wp) ::  ddt_lsf !< inverse value of time resolution of forcing data
    299        REAL(wp) ::  t_ref   !< time past since last reference step
    300      
    301 !
    302 !--    Calculate time interval of forcing data       
    303        ddt_lsf = 1.0_wp / ( nest_offl%time(nest_offl%tind_p) -                 &
    304                             nest_offl%time(nest_offl%tind) )
    305 !
    306 !--    Calculate reziproke time past since last reference step. Please note,
    307 !--    the time coordinate is still not updated, so that the actual time need
    308 !--    to be incremented by dt_3d. Moreover, note that the simulation time
    309 !--    passed since simulation start is time_since_reference_point, not
    310 !--    simulated_time!
    311        t_ref = time_since_reference_point + dt_3d -                            &
    312                                             nest_offl%time(nest_offl%tind)
    313                                            
    314        IF ( bc_dirichlet_l )  THEN
    315 
    316           DO  j = nys, nyn
    317              DO  k = nzb+1, nzt
    318                 u(k,j,nxlg:nxl) = nest_offl%u_left(0,k,j) + ddt_lsf * t_ref  * &
    319                        ( nest_offl%u_left(1,k,j) - nest_offl%u_left(0,k,j) ) * &
    320                          MERGE( 1.0_wp, 0.0_wp,                                &
    321                                 BTEST( wall_flags_0(k,j,nxlg:nxl), 1 ) )
    322              ENDDO
    323           ENDDO
    324 
    325           DO  j = nys, nyn
    326              DO  k = nzb+1, nzt-1
    327                 w(k,j,nxlg:nxl-1) = nest_offl%w_left(0,k,j) + ddt_lsf * t_ref *&
    328                        ( nest_offl%w_left(1,k,j) - nest_offl%w_left(0,k,j) )  *&
    329                          MERGE( 1.0_wp, 0.0_wp,                                &
    330                                 BTEST( wall_flags_0(k,j,nxlg:nxl-1), 3 ) )
    331              ENDDO
    332           ENDDO
    333 
    334           DO  j = nysv, nyn
    335              DO  k = nzb+1, nzt
    336                 v(k,j,nxlg:nxl-1) = nest_offl%v_left(0,k,j) + ddt_lsf * t_ref *&
    337                        ( nest_offl%v_left(1,k,j) - nest_offl%v_left(0,k,j) )  *&
    338                          MERGE( 1.0_wp, 0.0_wp,                                &
    339                                 BTEST( wall_flags_0(k,j,nxlg:nxl-1), 2 ) )
    340              ENDDO
    341           ENDDO
    342 
    343           IF ( .NOT. neutral )  THEN
    344              DO  j = nys, nyn
    345                 DO  k = nzb+1, nzt
    346                    pt(k,j,nxlg:nxl-1) = nest_offl%pt_left(0,k,j) + ddt_lsf *   &
    347                                                                    t_ref   *   &
    348                        ( nest_offl%pt_left(1,k,j) - nest_offl%pt_left(0,k,j) )
    349  
    350                 ENDDO
    351              ENDDO
    352           ENDIF
    353 
    354           IF ( humidity )  THEN
    355              DO  j = nys, nyn
    356                 DO  k = nzb+1, nzt
    357                    q(k,j,nxlg:nxl-1) = nest_offl%q_left(0,k,j) + ddt_lsf   *   &
    358                                                                  t_ref     *   &
    359                        ( nest_offl%q_left(1,k,j) - nest_offl%q_left(0,k,j) )
    360  
    361                 ENDDO
    362              ENDDO
    363           ENDIF
    364 
    365        ENDIF
    366 
    367        IF ( bc_dirichlet_r )  THEN
    368 
    369           DO  j = nys, nyn
    370              DO  k = nzb+1, nzt
    371                 u(k,j,nxr+1:nxrg) = nest_offl%u_right(0,k,j) + ddt_lsf * t_ref *&
    372                       ( nest_offl%u_right(1,k,j) - nest_offl%u_right(0,k,j) )  *&
    373                          MERGE( 1.0_wp, 0.0_wp,                                 &
    374                                 BTEST( wall_flags_0(k,j,nxr+1:nxrg), 1 ) )
    375 
    376              ENDDO
    377           ENDDO
    378           DO  j = nys, nyn
    379              DO  k = nzb+1, nzt-1
    380                 w(k,j,nxr+1:nxrg) = nest_offl%w_right(0,k,j) + ddt_lsf * t_ref *&
    381                       ( nest_offl%w_right(1,k,j) - nest_offl%w_right(0,k,j) )  *&
    382                          MERGE( 1.0_wp, 0.0_wp,                                 &
    383                                 BTEST( wall_flags_0(k,j,nxr+1:nxrg), 3 ) )
    384              ENDDO
    385           ENDDO
    386 
    387           DO  j = nysv, nyn
    388              DO  k = nzb+1, nzt
    389                 v(k,j,nxr+1:nxrg) = nest_offl%v_right(0,k,j) + ddt_lsf * t_ref *&
    390                       ( nest_offl%v_right(1,k,j) - nest_offl%v_right(0,k,j) )  *&
    391                          MERGE( 1.0_wp, 0.0_wp,                                 &
    392                                 BTEST( wall_flags_0(k,j,nxr+1:nxrg), 2 ) )
    393              ENDDO
    394           ENDDO
    395 
    396           IF ( .NOT. neutral )  THEN
    397              DO  j = nys, nyn
    398                 DO  k = nzb+1, nzt
    399                    pt(k,j,nxr+1:nxrg) = nest_offl%pt_right(0,k,j) + ddt_lsf *  &
    400                                                                     t_ref   *  &
    401                      ( nest_offl%pt_right(1,k,j) - nest_offl%pt_right(0,k,j) )
    402  
    403                 ENDDO
    404              ENDDO
    405           ENDIF
    406 
    407           IF ( humidity )  THEN
    408              DO  j = nys, nyn
    409                 DO  k = nzb+1, nzt
    410                    q(k,j,nxr+1:nxrg) = nest_offl%q_right(0,k,j) + ddt_lsf   *  &
    411                                                                     t_ref   *  &
    412                        ( nest_offl%q_right(1,k,j) - nest_offl%q_right(0,k,j) )
    413  
    414                 ENDDO
    415              ENDDO
    416           ENDIF
    417 
    418        ENDIF
    419 
    420        IF ( bc_dirichlet_s )  THEN
    421 
    422           DO  i = nxl, nxr
    423              DO  k = nzb+1, nzt
    424                 v(k,nysg:nys,i)   = nest_offl%v_south(0,k,i) + ddt_lsf * t_ref *&
    425                       ( nest_offl%v_south(1,k,i) - nest_offl%v_south(0,k,i) )  *&
    426                          MERGE( 1.0_wp, 0.0_wp,                                 &
    427                                 BTEST( wall_flags_0(k,nysg:nys,i), 2 ) )
    428              ENDDO
    429           ENDDO
    430 
    431           DO  i = nxl, nxr
    432              DO  k = nzb+1, nzt-1
    433                 w(k,nysg:nys-1,i) = nest_offl%w_south(0,k,i) + ddt_lsf * t_ref  *&
    434                         ( nest_offl%w_south(1,k,i) - nest_offl%w_south(0,k,i) ) *&
    435                            MERGE( 1.0_wp, 0.0_wp,                                &
    436                                   BTEST( wall_flags_0(k,nysg:nys-1,i), 3 ) )
    437              ENDDO
    438           ENDDO
    439 
    440           DO  i = nxlu, nxr
    441              DO  k = nzb+1, nzt
    442                 u(k,nysg:nys-1,i) = nest_offl%u_south(0,k,i) + ddt_lsf * t_ref  *&
    443                         ( nest_offl%u_south(1,k,i) - nest_offl%u_south(0,k,i) ) *&
    444                            MERGE( 1.0_wp, 0.0_wp,                                &
    445                                   BTEST( wall_flags_0(k,nysg:nys-1,i), 1 ) )
    446              ENDDO
    447           ENDDO
    448 
    449           IF ( .NOT. neutral )  THEN
    450              DO  i = nxl, nxr
    451                 DO  k = nzb+1, nzt
    452                    pt(k,nysg:nys-1,i) = nest_offl%pt_south(0,k,i) + ddt_lsf *  &
    453                                                                     t_ref   *  &
    454                      ( nest_offl%pt_south(1,k,i) - nest_offl%pt_south(0,k,i) )
    455  
    456                 ENDDO
    457              ENDDO
    458           ENDIF
    459 
    460           IF ( humidity )  THEN
    461              DO  i = nxl, nxr
    462                 DO  k = nzb+1, nzt
    463                    q(k,nysg:nys-1,i) = nest_offl%q_south(0,k,i) + ddt_lsf   *  &
    464                                                                     t_ref   *  &
    465                        ( nest_offl%q_south(1,k,i) - nest_offl%q_south(0,k,i) )
    466  
    467                 ENDDO
    468              ENDDO
    469           ENDIF
    470 
    471        ENDIF
    472 
    473        IF ( bc_dirichlet_n )  THEN
    474 
    475           DO  i = nxl, nxr
    476              DO  k = nzb+1, nzt
    477                 v(k,nyn+1:nyng,i)   = nest_offl%v_north(0,k,i) + ddt_lsf * t_ref *&
    478                         ( nest_offl%v_north(1,k,i) - nest_offl%v_north(0,k,i) )  *&
    479                            MERGE( 1.0_wp, 0.0_wp,                                 &
    480                                   BTEST( wall_flags_0(k,nyn+1:nyng,i), 2 ) )
    481              ENDDO
    482           ENDDO
    483           DO  i = nxl, nxr
    484              DO  k = nzb+1, nzt-1
    485                 w(k,nyn+1:nyng,i) = nest_offl%w_north(0,k,i) + ddt_lsf * t_ref  *&
    486                         ( nest_offl%w_north(1,k,i) - nest_offl%w_north(0,k,i) ) *&
    487                            MERGE( 1.0_wp, 0.0_wp,                                &
    488                                   BTEST( wall_flags_0(k,nyn+1:nyng,i), 3 ) )
    489              ENDDO
    490           ENDDO
    491 
    492           DO  i = nxlu, nxr
    493              DO  k = nzb+1, nzt
    494                 u(k,nyn+1:nyng,i) = nest_offl%u_north(0,k,i) + ddt_lsf * t_ref  *&
    495                         ( nest_offl%u_north(1,k,i) - nest_offl%u_north(0,k,i) ) *&
    496                            MERGE( 1.0_wp, 0.0_wp,                                &
    497                                   BTEST( wall_flags_0(k,nyn+1:nyng,i), 1 ) )
    498 
    499              ENDDO
    500           ENDDO
    501 
    502           IF ( .NOT. neutral )  THEN
    503              DO  i = nxl, nxr
    504                 DO  k = nzb+1, nzt
    505                    pt(k,nyn+1:nyng,i) = nest_offl%pt_north(0,k,i) + ddt_lsf *  &
    506                                                                     t_ref   *  &
    507                      ( nest_offl%pt_north(1,k,i) - nest_offl%pt_north(0,k,i) )
    508  
    509                 ENDDO
    510              ENDDO
    511           ENDIF
    512 
    513           IF ( humidity )  THEN
    514              DO  i = nxl, nxr
    515                 DO  k = nzb+1, nzt
    516                    q(k,nyn+1:nyng,i) = nest_offl%q_north(0,k,i) + ddt_lsf   *  &
    517                                                                     t_ref   *  &
    518                        ( nest_offl%q_north(1,k,i) - nest_offl%q_north(0,k,i) )
    519  
    520                 ENDDO
    521              ENDDO
    522           ENDIF
    523 
    524        ENDIF
    525 !
    526 !--    Top boundary.
    527        DO  i = nxlu, nxr
    528           DO  j = nys, nyn
    529              u(nzt+1,j,i) = nest_offl%u_top(0,j,i) + ddt_lsf * t_ref *         &
    530                         ( nest_offl%u_top(1,j,i) - nest_offl%u_top(0,j,i) ) *  &
    531                            MERGE( 1.0_wp, 0.0_wp,                              &
    532                                   BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
    533           ENDDO
    534        ENDDO
    535 
    536        DO  i = nxl, nxr
    537           DO  j = nysv, nyn
    538              v(nzt+1,j,i) = nest_offl%v_top(0,j,i) + ddt_lsf * t_ref *         &
    539                         ( nest_offl%v_top(1,j,i) - nest_offl%v_top(0,j,i) ) *  &
    540                            MERGE( 1.0_wp, 0.0_wp,                              &
    541                                   BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
    542           ENDDO
    543        ENDDO
    544 
    545        DO  i = nxl, nxr
    546           DO  j = nys, nyn
    547              w(nzt:nzt+1,j,i) = nest_offl%w_top(0,j,i) + ddt_lsf * t_ref *     &
    548                         ( nest_offl%w_top(1,j,i) - nest_offl%w_top(0,j,i) ) *  &
    549                            MERGE( 1.0_wp, 0.0_wp,                              &
    550                                   BTEST( wall_flags_0(nzt:nzt+1,j,i), 3 ) )
    551           ENDDO
    552        ENDDO
    553 
    554 
    555        IF ( .NOT. neutral )  THEN
    556           DO  i = nxl, nxr
    557              DO  j = nys, nyn
    558                 pt(nzt+1,j,i) = nest_offl%pt_top(0,j,i) + ddt_lsf * t_ref *    &
    559                         ( nest_offl%pt_top(1,j,i) - nest_offl%pt_top(0,j,i) )
    560              ENDDO
    561           ENDDO
    562        ENDIF
    563 
    564        IF ( humidity )  THEN
    565           DO  i = nxl, nxr
    566              DO  j = nys, nyn
    567                 q(nzt+1,j,i) = nest_offl%q_top(0,j,i) + ddt_lsf * t_ref *      &
    568                         ( nest_offl%q_top(1,j,i) - nest_offl%q_top(0,j,i) )
    569              ENDDO
    570           ENDDO
    571        ENDIF
    572 !
    573 !--    At the edges( left-south, left-north, right-south and right-north) set
    574 !--    data on ghost points.
    575        IF ( bc_dirichlet_l  .AND.  bc_dirichlet_s )  THEN
    576           DO  i = 1, nbgp
    577              u(:,nys-i,nxlg:nxl)   = u(:,nys,nxlg:nxl)
    578              w(:,nys-i,nxlg:nxl-1) = w(:,nys,nxlg:nxl-1)
    579              IF ( .NOT. neutral )  pt(:,nys-i,nxlg:nxl-1) = pt(:,nys,nxlg:nxl-1)
    580              IF ( humidity      )  q(:,nys-i,nxlg:nxl-1)  = q(:,nys,nxlg:nxl-1)
    581           ENDDO
    582           DO  i = 1, nbgp+1
    583              v(:,nysv-i,nxlg:nxl-1) = v(:,nysv,nxlg:nxl-1)
    584           ENDDO
    585        ENDIF
    586        IF ( bc_dirichlet_l  .AND.  bc_dirichlet_n )  THEN
    587           DO  i = 1, nbgp
    588              u(:,nyn+i,nxlg:nxl)   = u(:,nyn,nxlg:nxl)
    589              v(:,nyn+i,nxlg:nxl-1) = v(:,nyn,nxlg:nxl-1)
    590              w(:,nyn+i,nxlg:nxl-1) = w(:,nyn,nxlg:nxl-1)
    591              IF ( .NOT. neutral )  pt(:,nyn+i,nxlg:nxl-1) = pt(:,nyn,nxlg:nxl-1)
    592              IF ( humidity      )  q(:,nyn+i,nxlg:nxl-1)  = q(:,nyn,nxlg:nxl-1)
    593           ENDDO
    594        ENDIF
    595        IF ( bc_dirichlet_r  .AND.  bc_dirichlet_s )  THEN
    596           DO  i = 1, nbgp
    597              u(:,nys-i,nxr+1:nxrg) = u(:,nys,nxr+1:nxrg)
    598              w(:,nys-i,nxr+1:nxrg) = w(:,nys,nxr+1:nxrg)
    599              IF ( .NOT. neutral )  pt(:,nys-i,nxr+1:nxrg) = pt(:,nys,nxr+1:nxrg)
    600              IF ( humidity      )  q(:,nys-i,nxr+1:nxrg)  = q(:,nys,nxr+1:nxrg)
    601           ENDDO
    602           DO  i = 1, nbgp+1
    603              v(:,nysv-i,nxr+1:nxrg) = v(:,nysv,nxr+1:nxrg)
    604           ENDDO
    605        ENDIF
    606        IF ( bc_dirichlet_r  .AND.  bc_dirichlet_n )  THEN
    607           DO  i = 1, nbgp
    608              u(:,nyn+i,nxr+1:nxrg) = u(:,nyn,nxr+1:nxrg)
    609              v(:,nyn+i,nxr+1:nxrg) = v(:,nyn,nxr+1:nxrg)
    610              w(:,nyn+i,nxr+1:nxrg) = w(:,nyn,nxr+1:nxrg)
    611              IF ( .NOT. neutral )  pt(:,nyn+i,nxr+1:nxrg) = pt(:,nyn,nxr+1:nxrg)
    612              IF ( humidity      )  q(:,nyn+i,nxr+1:nxrg)  = q(:,nyn,nxr+1:nxrg)
    613           ENDDO
    614        ENDIF
    615 !
    616 !--    Moreover, set Neumann boundary condition for subgrid-scale TKE,
    617 !--    passive scalar, dissipation, and chemical species if required
    618        IF ( rans_mode  .AND.  rans_tke_e )  THEN
    619           IF (  bc_dirichlet_l )  diss(:,:,nxl-1) = diss(:,:,nxl)
    620           IF (  bc_dirichlet_r )  diss(:,:,nxr+1) = diss(:,:,nxr)
    621           IF (  bc_dirichlet_s )  diss(:,nys-1,:) = diss(:,nys,:)
    622           IF (  bc_dirichlet_n )  diss(:,nyn+1,:) = diss(:,nyn,:)
    623        ENDIF
    624        IF ( .NOT. constant_diffusion )  THEN
    625           IF (  bc_dirichlet_l )  e(:,:,nxl-1) = e(:,:,nxl)
    626           IF (  bc_dirichlet_r )  e(:,:,nxr+1) = e(:,:,nxr)
    627           IF (  bc_dirichlet_s )  e(:,nys-1,:) = e(:,nys,:)
    628           IF (  bc_dirichlet_n )  e(:,nyn+1,:) = e(:,nyn,:)
    629           e(nzt+1,:,:) = e(nzt,:,:)
    630        ENDIF
    631        IF ( passive_scalar )  THEN
    632           IF (  bc_dirichlet_l )  s(:,:,nxl-1) = s(:,:,nxl)
    633           IF (  bc_dirichlet_r )  s(:,:,nxr+1) = s(:,:,nxr)
    634           IF (  bc_dirichlet_s )  s(:,nys-1,:) = s(:,nys,:)
    635           IF (  bc_dirichlet_n )  s(:,nyn+1,:) = s(:,nyn,:)
    636        ENDIF
    637 
    638 
    639 
    640        CALL exchange_horiz( u, nbgp )
    641        CALL exchange_horiz( v, nbgp )
    642        CALL exchange_horiz( w, nbgp )
    643        IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
    644        IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
    645 
    646 !
    647 !--    Set surface pressure. Please note, time-dependent surface
    648 !--    pressure would require changes in anelastic approximation and
    649 !--    treatment of fluxes.
    650 !--    For the moment, comment this out!
    651 !      surface_pressure = nest_offl%surface_pressure(nest_offl%tind) +                 &
    652 !                                                      ddt_lsf * t_ref *       &
    653 !                                    ( nest_offl%surface_pressure(nest_offl%tind_p)    &
    654 !                                    - nest_offl%surface_pressure(nest_offl%tind) )
    655 
    656     END SUBROUTINE lsf_nesting_offline
    657178
    658179!------------------------------------------------------------------------------!
     
    935456           ONLY:  bc_lr_cyc, bc_ns_cyc
    936457
    937        USE netcdf_data_input_mod,                                              &
    938            ONLY:  netcdf_data_input_lsf 
    939 
    940458       IMPLICIT NONE
    941459
     
    967485       REAL(wp) ::  r_dummy           !<
    968486
    969        IF ( nesting_offline )  THEN
    970 !
    971 !--       Allocate arrays for geostrophic wind components. Arrays will
    972 !--       incorporate 2 time levels in order to interpolate in between. Please
    973 !--       note, forcing using geostrophic wind components is only required in
    974 !--       case of cyclic boundary conditions.
    975           IF ( bc_lr_cyc  .AND.  bc_ns_cyc )  THEN
    976              ALLOCATE( nest_offl%ug(0:1,nzb:nzt+1) )
    977              ALLOCATE( nest_offl%vg(0:1,nzb:nzt+1) )
     487       ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),              &
     488                 qsws_surf(0:nlsf), shf_surf(0:nlsf),                          &
     489                 td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf),     &
     490                 td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf),     &
     491                 time_vert(0:nlsf), time_surf(0:nlsf),                         &
     492                 ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),         &
     493                 wsubs_vert(nzb:nzt+1,0:nlsf) )
     494
     495       p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
     496       shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
     497       td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
     498       time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
     499       wsubs_vert = 0.0_wp
     500
     501!
     502!--    Array for storing large scale forcing and nudging tendencies at each
     503!--    timestep for data output
     504       ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
     505       sums_ls_l = 0.0_wp
     506
     507       ngp_sums_ls = (nz+2)*6
     508
     509       OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
     510              FORM='FORMATTED', IOSTAT=ierrn )
     511
     512       IF ( ierrn /= 0 )  THEN
     513          message_string = 'file LSF_DATA does not exist'
     514          CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
     515       ENDIF
     516
     517       ierrn = 0
     518!
     519!--    First three lines of LSF_DATA contain header
     520       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
     521       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
     522       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
     523
     524       IF ( ierrn /= 0 )  THEN
     525          message_string = 'errors in file LSF_DATA'
     526          CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
     527       ENDIF
     528
     529!
     530!--    Surface values are read in
     531       nt     = 0
     532       ierrn = 0
     533
     534       DO WHILE ( time_surf(nt) < end_time )
     535          nt = nt + 1
     536          READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),      &
     537                                             qsws_surf(nt), pt_surf(nt),       &
     538                                             q_surf(nt), p_surf(nt)
     539
     540          IF ( ierrn /= 0 )  THEN
     541            WRITE ( message_string, * ) 'No time dependent surface ' //        &
     542                              'variables in & LSF_DATA for end of run found'
     543
     544             CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
    978545          ENDIF
    979 !
    980 !--       Allocate arrays for reading boundary values. Arrays will incorporate 2
    981 !--       time levels in order to interpolate in between.
    982           IF ( bc_dirichlet_l )  THEN
    983              ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn)  )
    984              ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
    985              ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
    986              IF ( humidity )      ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn)  )
    987              IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
    988           ENDIF
    989           IF ( bc_dirichlet_r )  THEN
    990              ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn)  )
    991              ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
    992              ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
    993              IF ( humidity )      ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn)  )
    994              IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
    995           ENDIF
    996           IF ( bc_dirichlet_n )  THEN
    997              ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
    998              ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr)  )
    999              ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
    1000              IF ( humidity )      ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr)  )
    1001              IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
    1002           ENDIF
    1003           IF ( bc_dirichlet_s )  THEN
    1004              ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
    1005              ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr)  )
    1006              ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr)    )
    1007              IF ( humidity )      ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr)  )
    1008              IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
    1009           ENDIF
    1010          
    1011           ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
    1012           ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
    1013           ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr)  )
    1014           IF ( humidity )      ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr)  )
    1015           IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
    1016 
    1017 !
    1018 !--       Read COSMO data at lateral and top boundaries
    1019           CALL netcdf_data_input_lsf
    1020 !
    1021 !--       Write COSMO data at lateral and top boundaries
    1022           CALL lsf_nesting_offline
    1023 !
    1024 !--       After 3D data is initialized, ensure mass conservation
    1025           CALL lsf_nesting_offline_mass_conservation
    1026 !
    1027 !--       Initialize surface pressure. Please note, time-dependent surface
    1028 !--       pressure would require changes in anelastic approximation and
    1029 !--       treatment of fluxes.
    1030 !--       For the moment, comment this out!
    1031 !         surface_pressure = nest_offl%surface_pressure(0)
    1032 
    1033        ELSE
    1034 
    1035           ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),           &
    1036                     qsws_surf(0:nlsf), shf_surf(0:nlsf),                       &
    1037                     td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf),  &
    1038                     td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf),  &
    1039                     time_vert(0:nlsf), time_surf(0:nlsf),                      &
    1040                     ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),      &
    1041                     wsubs_vert(nzb:nzt+1,0:nlsf) )
    1042 
    1043           p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
    1044           shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
    1045           td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
    1046           time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
    1047           wsubs_vert = 0.0_wp
    1048 
    1049 !
    1050 !--       Array for storing large scale forcing and nudging tendencies at each
    1051 !--       timestep for data output
    1052           ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
    1053           sums_ls_l = 0.0_wp
    1054 
    1055           ngp_sums_ls = (nz+2)*6
    1056 
    1057           OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
    1058                  FORM='FORMATTED', IOSTAT=ierrn )
    1059 
    1060           IF ( ierrn /= 0 )  THEN
    1061              message_string = 'file LSF_DATA does not exist'
    1062              CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
    1063           ENDIF
    1064 
    1065           ierrn = 0
    1066 !
    1067 !--       First three lines of LSF_DATA contain header
    1068           READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
    1069           READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
    1070           READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
    1071 
     546       ENDDO
     547
     548       IF ( time_surf(1) > end_time )  THEN
     549          WRITE ( message_string, * ) 'Time dependent surface variables in ' //&
     550                                      '&LSF_DATA set in after end of ' ,       &
     551                                      'simulation - lsf_surf is set to FALSE'
     552          CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
     553          lsf_surf = .FALSE.
     554       ENDIF
     555
     556!
     557!--    Go to the end of the list with surface variables
     558       DO WHILE ( ierrn == 0 )
     559          READ ( finput, *, IOSTAT = ierrn ) r_dummy
     560       ENDDO
     561
     562!
     563!--    Profiles of ug, vg and w_subs are read in (large scale forcing)
     564
     565       nt = 0
     566       DO WHILE ( time_vert(nt) < end_time )
     567          nt = nt + 1
     568          hash = "#"
     569          ierrn = 1 ! not zero
     570!
     571!--       Search for the next line consisting of "# time",
     572!--       from there onwards the profiles will be read
     573          DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
     574             READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
     575             IF ( ierrn < 0 )  THEN
     576                WRITE( message_string, * ) 'No time dependent vertical profiles',&
     577                                 ' in & LSF_DATA for end of run found'
     578                CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
     579             ENDIF
     580          ENDDO
     581
     582          IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
     583
     584          READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,  &
     585                                           lowwsubs_vert, low_td_lsa_lpt,      &
     586                                           low_td_lsa_q, low_td_sub_lpt,       &
     587                                           low_td_sub_q
    1072588          IF ( ierrn /= 0 )  THEN
    1073589             message_string = 'errors in file LSF_DATA'
     
    1075591          ENDIF
    1076592
    1077 !
    1078 !--       Surface values are read in
    1079           nt     = 0
    1080           ierrn = 0
    1081 
    1082           DO WHILE ( time_surf(nt) < end_time )
    1083              nt = nt + 1
    1084              READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),   &
    1085                                                 qsws_surf(nt), pt_surf(nt),    &
    1086                                                 q_surf(nt), p_surf(nt)
    1087 
    1088              IF ( ierrn /= 0 )  THEN
    1089                WRITE ( message_string, * ) 'No time dependent surface ' //     &
    1090                                  'variables in & LSF_DATA for end of run found'
    1091 
    1092                 CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
     593          READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,            &
     594                                           highvg_vert, highwsubs_vert,        &
     595                                           high_td_lsa_lpt, high_td_lsa_q,     &
     596                                           high_td_sub_lpt, high_td_sub_q
     597       
     598          IF ( ierrn /= 0 )  THEN
     599             message_string = 'errors in file LSF_DATA'
     600             CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
     601          ENDIF
     602
     603
     604          DO  k = nzb, nzt+1
     605             IF ( highheight < zu(k) )  THEN
     606                lowheight      = highheight
     607                lowug_vert     = highug_vert
     608                lowvg_vert     = highvg_vert
     609                lowwsubs_vert  = highwsubs_vert
     610                low_td_lsa_lpt = high_td_lsa_lpt
     611                low_td_lsa_q   = high_td_lsa_q
     612                low_td_sub_lpt = high_td_sub_lpt
     613                low_td_sub_q   = high_td_sub_q
     614
     615                ierrn = 0
     616                READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,      &
     617                                                 highvg_vert, highwsubs_vert,  &
     618                                                 high_td_lsa_lpt,              &
     619                                                 high_td_lsa_q,                &
     620                                                 high_td_sub_lpt, high_td_sub_q
     621
     622                IF ( ierrn /= 0 )  THEN
     623                   WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ',     &
     624                        'is higher than the maximum height in LSF_DATA ',      &
     625                        'which is ', lowheight, 'm. Interpolation on PALM ',   &
     626                        'grid is not possible.'
     627                   CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
     628                ENDIF
     629
    1093630             ENDIF
     631
     632!
     633!--          Interpolation of prescribed profiles in space
     634             fac = (highheight-zu(k))/(highheight - lowheight)
     635
     636             ug_vert(k,nt)    = fac * lowug_vert                               &
     637                                + ( 1.0_wp - fac ) * highug_vert
     638             vg_vert(k,nt)    = fac * lowvg_vert                               &
     639                                + ( 1.0_wp - fac ) * highvg_vert
     640             wsubs_vert(k,nt) = fac * lowwsubs_vert                            &
     641                                + ( 1.0_wp - fac ) * highwsubs_vert
     642
     643             td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt                           &
     644                                + ( 1.0_wp - fac ) * high_td_lsa_lpt
     645             td_lsa_q(k,nt)   = fac * low_td_lsa_q                             &
     646                                + ( 1.0_wp - fac ) * high_td_lsa_q
     647             td_sub_lpt(k,nt) = fac * low_td_sub_lpt                           &
     648                                + ( 1.0_wp - fac ) * high_td_sub_lpt
     649             td_sub_q(k,nt)   = fac * low_td_sub_q                             &
     650                                + ( 1.0_wp - fac ) * high_td_sub_q
     651
    1094652          ENDDO
    1095653
    1096           IF ( time_surf(1) > end_time )  THEN
    1097              WRITE ( message_string, * ) 'Time dependent surface variables in ' // &
    1098                                          '&LSF_DATA set in after end of ' ,        &
    1099                                          'simulation - lsf_surf is set to FALSE'
    1100              CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
    1101              lsf_surf = .FALSE.
    1102           ENDIF
    1103 
    1104 !
    1105 !--       Go to the end of the list with surface variables
    1106           DO WHILE ( ierrn == 0 )
    1107              READ ( finput, *, IOSTAT = ierrn ) r_dummy
    1108           ENDDO
    1109 
    1110 !
    1111 !--       Profiles of ug, vg and w_subs are read in (large scale forcing)
    1112 
    1113           nt = 0
    1114           DO WHILE ( time_vert(nt) < end_time )
    1115              nt = nt + 1
    1116              hash = "#"
    1117              ierrn = 1 ! not zero
    1118 !
    1119 !--          Search for the next line consisting of "# time",
    1120 !--          from there onwards the profiles will be read
    1121              DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
    1122                 READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
    1123                 IF ( ierrn < 0 )  THEN
    1124                    WRITE( message_string, * ) 'No time dependent vertical profiles',&
    1125                                     ' in & LSF_DATA for end of run found'
    1126                    CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
    1127                 ENDIF
    1128              ENDDO
    1129 
    1130              IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
    1131 
    1132              READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,&
    1133                                               lowwsubs_vert, low_td_lsa_lpt,    &
    1134                                               low_td_lsa_q, low_td_sub_lpt,     &
    1135                                               low_td_sub_q
    1136              IF ( ierrn /= 0 )  THEN
    1137                 message_string = 'errors in file LSF_DATA'
    1138                 CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
    1139              ENDIF
    1140 
    1141              READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,         &
    1142                                               highvg_vert, highwsubs_vert,     &
    1143                                               high_td_lsa_lpt, high_td_lsa_q,  &
    1144                                               high_td_sub_lpt, high_td_sub_q
    1145          
    1146              IF ( ierrn /= 0 )  THEN
    1147                 message_string = 'errors in file LSF_DATA'
    1148                 CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
    1149              ENDIF
    1150 
    1151 
    1152              DO  k = nzb, nzt+1
    1153                 IF ( highheight < zu(k) )  THEN
    1154                    lowheight      = highheight
    1155                    lowug_vert     = highug_vert
    1156                    lowvg_vert     = highvg_vert
    1157                    lowwsubs_vert  = highwsubs_vert
    1158                    low_td_lsa_lpt = high_td_lsa_lpt
    1159                    low_td_lsa_q   = high_td_lsa_q
    1160                    low_td_sub_lpt = high_td_sub_lpt
    1161                    low_td_sub_q   = high_td_sub_q
    1162 
    1163                    ierrn = 0
    1164                    READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,    &
    1165                                                     highvg_vert, highwsubs_vert,&
    1166                                                     high_td_lsa_lpt,            &
    1167                                                     high_td_lsa_q,              &
    1168                                                     high_td_sub_lpt, high_td_sub_q
    1169 
    1170                    IF ( ierrn /= 0 )  THEN
    1171                       WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ',  &
    1172                            'is higher than the maximum height in LSF_DATA ',   &
    1173                            'which is ', lowheight, 'm. Interpolation on PALM ',&
    1174                            'grid is not possible.'
    1175                       CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
    1176                    ENDIF
    1177 
    1178                 ENDIF
    1179 
    1180 !
    1181 !--             Interpolation of prescribed profiles in space
    1182                 fac = (highheight-zu(k))/(highheight - lowheight)
    1183 
    1184                 ug_vert(k,nt)    = fac * lowug_vert                            &
    1185                                    + ( 1.0_wp - fac ) * highug_vert
    1186                 vg_vert(k,nt)    = fac * lowvg_vert                            &
    1187                                    + ( 1.0_wp - fac ) * highvg_vert
    1188                 wsubs_vert(k,nt) = fac * lowwsubs_vert                         &
    1189                                    + ( 1.0_wp - fac ) * highwsubs_vert
    1190 
    1191                 td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt                        &
    1192                                    + ( 1.0_wp - fac ) * high_td_lsa_lpt
    1193                 td_lsa_q(k,nt)   = fac * low_td_lsa_q                          &
    1194                                    + ( 1.0_wp - fac ) * high_td_lsa_q
    1195                 td_sub_lpt(k,nt) = fac * low_td_sub_lpt                        &
    1196                                    + ( 1.0_wp - fac ) * high_td_sub_lpt
    1197                 td_sub_q(k,nt)   = fac * low_td_sub_q                          &
    1198                                    + ( 1.0_wp - fac ) * high_td_sub_q
    1199 
    1200              ENDDO
    1201 
    1202           ENDDO
    1203 
    1204 !
    1205 !--       Large scale vertical velocity has to be zero at the surface
    1206           wsubs_vert(nzb,:) = 0.0_wp
     654       ENDDO
     655
     656!
     657!--    Large scale vertical velocity has to be zero at the surface
     658       wsubs_vert(nzb,:) = 0.0_wp
    1207659   
    1208           IF ( time_vert(1) > end_time )  THEN
    1209              WRITE ( message_string, * ) 'Time dependent large scale profile ',&
    1210                                 'forcing from&LSF_DATA sets in after end of ' ,&
    1211                                 'simulation - lsf_vert is set to FALSE'
    1212              CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
    1213              lsf_vert = .FALSE.
    1214           ENDIF
    1215 
    1216           CLOSE( finput )
    1217 
    1218        ENDIF
     660       IF ( time_vert(1) > end_time )  THEN
     661          WRITE ( message_string, * ) 'Time dependent large scale profile ',   &
     662                             'forcing from&LSF_DATA sets in after end of ' ,   &
     663                             'simulation - lsf_vert is set to FALSE'
     664          CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
     665          lsf_vert = .FALSE.
     666       ENDIF
     667
     668       CLOSE( finput )
    1219669
    1220670    END SUBROUTINE lsf_init
Note: See TracChangeset for help on using the changeset viewer.