Changeset 2696 for palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90
- Timestamp:
- Dec 14, 2017 5:12:51 PM (6 years ago)
- 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 l sf_nudging_mod.f902 !------------------------------------------------------------------------------! 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. 4 4 ! 5 5 ! PALM is free software: you can redistribute it and/or modify it under the … … 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Forcing with larger-scale models implemented (MS) 28 ! 29 ! 2342 2017-08-08 11:00:43Z boeske 27 30 ! fixed check if surface forcing data is available until end of simulation 28 31 ! … … 46 49 47 50 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 51 54 52 55 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, & 54 59 humidity, intermediate_timestep_count, ibc_pt_b, ibc_q_b, & 55 60 large_scale_forcing, large_scale_subsidence, lsf_surf, lsf_vert,& 56 lsf_exception, message_string, n udging, passive_scalar,&61 lsf_exception, message_string, neutral, nudging, passive_scalar,& 57 62 pt_surface, ocean, q_surface, surface_pressure, topography, & 58 63 use_subsidence_tendencies 59 64 65 USE grid_variables 66 67 USE pegrid 60 68 61 69 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 63 72 64 73 USE kinds … … 70 79 ONLY: hom, statistic_regions, sums_ls_l, weight_substep 71 80 81 USE netcdf_data_input_mod, & 82 ONLY: force, netcdf_data_input_interpolate 72 83 73 84 INTEGER(iwp) :: nlsf = 1000 !< maximum number of profiles in LSF_DATA (large scale forcing) 74 85 INTEGER(iwp) :: ntnudge = 1000 !< maximum number of profiles in NUDGING_DATA (nudging) 86 87 REAL(wp) :: d_area_t 75 88 76 89 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ptnudge !< vertical profile of pot. temperature interpolated to vertical grid (nudging) … … 107 120 lsf_nudging_check_parameters, nudge_init, & 108 121 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 110 124 ! 111 125 !-- Public variables 112 126 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 114 129 115 130 INTERFACE ls_advec … … 125 140 CONTAINS 126 141 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 127 661 128 662 !------------------------------------------------------------------------------! … … 403 937 SUBROUTINE lsf_init 404 938 939 USE netcdf_data_input_mod, & 940 ONLY: netcdf_data_input_lsf 941 405 942 IMPLICIT NONE 406 943 … … 411 948 INTEGER(iwp) :: finput = 90 !< 412 949 INTEGER(iwp) :: k !< 413 INTEGER(iwp) :: nt !< 950 INTEGER(iwp) :: nt !< 951 INTEGER(iwp) :: t !< running index for time levels 414 952 415 953 REAL(wp) :: fac !< … … 432 970 REAL(wp) :: r_dummy !< 433 971 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 ) 486 1179 487 1180 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 535 1192 IF ( ierrn /= 0 ) THEN 536 1193 message_string = 'errors in file LSF_DATA' … … 538 1195 ENDIF 539 1196 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 ) 575 1247 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 ) 577 1259 ENDIF 578 1260 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 617 1339 618 1340 END SUBROUTINE lsf_init … … 1396 2118 END SUBROUTINE nudge_ref 1397 2119 2120 1398 2121 END MODULE lsf_nudging_mod
Note: See TracChangeset
for help on using the changeset viewer.
