Changeset 1365 for palm/trunk/SOURCE/ls_forcing.f90
- Timestamp:
- Apr 22, 2014 3:03:56 PM (10 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/ls_forcing.f90 (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/ls_forcing.f90
r1354 r1365 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Usage of large scale forcing for pt and q enabled: 23 ! Added new subroutine ls_advec for horizontal large scale advection and large 24 ! scale subsidence, 25 ! error message in init_ls_forcing specified, 26 ! variable t renamed nt 23 27 ! 24 28 ! Former revisions: … … 61 65 62 66 PRIVATE 63 PUBLIC init_ls_forcing, ls_forcing_surf, ls_forcing_vert 67 PUBLIC init_ls_forcing, ls_forcing_surf, ls_forcing_vert, ls_advec 64 68 SAVE 65 69 70 INTERFACE ls_advec 71 MODULE PROCEDURE ls_advec 72 MODULE PROCEDURE ls_advec_ij 73 END INTERFACE ls_advec 66 74 67 75 CONTAINS … … 70 78 71 79 USE arrays_3d, & 72 ONLY: p_surf, pt_surf, q_surf, qsws_surf, shf_surf, time_surf, & 73 time_vert, ug_vert, vg_vert, wsubs_vert, zu 80 ONLY: p_surf, pt_lsa, pt_subs, pt_surf, q_lsa, q_subs, q_surf, & 81 qsws_surf, shf_surf, time_surf, time_vert, ug_vert, vg_vert, & 82 wsubs_vert, zu 74 83 75 84 USE control_parameters, & … … 77 86 78 87 USE indices, & 79 ONLY: n zb, nzt88 ONLY: ngp_sums_ls, nzb, nz, nzt 80 89 81 90 USE kinds 82 91 92 USE statistics, & 93 ONLY: sums_ls_l 94 83 95 84 96 IMPLICIT NONE 85 97 86 CHARACTER(100) :: chmess !: 87 CHARACTER(1) :: hash !: 88 89 INTEGER(iwp) :: ierrn !: 90 INTEGER(iwp) :: finput = 90 !: 91 INTEGER(iwp) :: k !: 92 INTEGER(iwp) :: t !: 93 94 REAL(wp) :: fac !: 95 REAL(wp) :: highheight !: 96 REAL(wp) :: highug_vert !: 97 REAL(wp) :: highvg_vert !: 98 REAL(wp) :: highwsubs_vert !: 99 REAL(wp) :: lowheight !: 100 REAL(wp) :: lowug_vert !: 101 REAL(wp) :: lowvg_vert !: 102 REAL(wp) :: lowwsubs_vert !: 103 REAL(wp) :: r_dummy !: 104 105 ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf), & 106 qsws_surf(0:nlsf), shf_surf(0:nlsf), time_vert(0:nlsf), & 107 time_surf(0:nlsf), ug_vert(nzb:nzt+1,0:nlsf), & 108 vg_vert(nzb:nzt+1,0:nlsf), wsubs_vert(nzb:nzt+1,0:nlsf) ) 109 110 p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp 111 qsws_surf = 0.0_wp; shf_surf = 0.0_wp; time_vert = 0.0_wp 112 time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp 113 wsubs_vert = 0.0_wp 114 98 CHARACTER(100) :: chmess !: 99 CHARACTER(1) :: hash !: 100 101 INTEGER(iwp) :: ierrn !: 102 INTEGER(iwp) :: finput = 90 !: 103 INTEGER(iwp) :: k !: 104 INTEGER(iwp) :: nt !: 105 106 REAL(wp) :: fac !: 107 REAL(wp) :: highheight !: 108 REAL(wp) :: highug_vert !: 109 REAL(wp) :: highvg_vert !: 110 REAL(wp) :: highwsubs_vert !: 111 REAL(wp) :: lowheight !: 112 REAL(wp) :: lowug_vert !: 113 REAL(wp) :: lowvg_vert !: 114 REAL(wp) :: lowwsubs_vert !: 115 REAL(wp) :: highpt_lsa !: 116 REAL(wp) :: lowpt_lsa !: 117 REAL(wp) :: highq_lsa !: 118 REAL(wp) :: lowq_lsa !: 119 REAL(wp) :: highpt_subs !: 120 REAL(wp) :: lowpt_subs !: 121 REAL(wp) :: highq_subs !: 122 REAL(wp) :: lowq_subs !: 123 REAL(wp) :: r_dummy !: 124 125 ALLOCATE( p_surf(0:nlsf), pt_lsa(nzb:nzt+1,0:nlsf), & 126 pt_subs(nzb:nzt+1,0:nlsf), pt_surf(0:nlsf), & 127 q_lsa(nzb:nzt+1,0:nlsf), q_subs(nzb:nzt+1,0:nlsf), & 128 q_surf(0:nlsf), qsws_surf(0:nlsf), shf_surf(0:nlsf), & 129 time_vert(0:nlsf), time_surf(0:nlsf), & 130 ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf), & 131 wsubs_vert(nzb:nzt+1,0:nlsf) ) 132 133 p_surf = 0.0_wp; pt_lsa = 0.0; pt_subs = 0.0; pt_surf = 0.0_wp 134 q_lsa = 0.0; q_subs = 0.0; q_surf = 0.0_wp; qsws_surf = 0.0_wp 135 shf_surf = 0.0_wp; time_vert = 0.0_wp; time_surf = 0.0_wp; 136 ug_vert = 0.0_wp; vg_vert = 0.0_wp; wsubs_vert = 0.0_wp 137 138 ! 139 !-- Array for storing large scale forcing and nudging tendencies at each 140 !-- timestep for data output 141 ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) ) 142 sums_ls_l = 0.0_wp 143 144 ngp_sums_ls = (nz+2)*6 115 145 116 146 OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', & … … 136 166 ! 137 167 !-- Surface values are read in 138 t = 0168 nt = 0 139 169 ierrn = 0 140 170 141 DO WHILE ( time_surf( t) < end_time )142 t =t + 1143 READ ( finput, *, IOSTAT = ierrn ) time_surf( t), shf_surf(t),&144 qsws_surf( t), pt_surf(t),&145 q_surf( t), p_surf(t)171 DO WHILE ( time_surf(nt) < end_time ) 172 nt = nt + 1 173 READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt), & 174 qsws_surf(nt), pt_surf(nt), & 175 q_surf(nt), p_surf(nt) 146 176 147 177 IF ( ierrn < 0 ) THEN … … 156 186 IF ( time_surf(1) > end_time ) THEN 157 187 WRITE ( message_string, * ) 'No time dependent surface variables in ',& 158 '&LSF_DATA for end of run found - ', &188 '&LSF_DATA for end of run found - ', & 159 189 'lsf_surf is set to FALSE' 160 190 CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 ) … … 171 201 !-- Profiles of ug, vg and w_subs are read in (large scale forcing) 172 202 173 t = 0174 DO WHILE ( time_vert( t) < end_time )175 t =t + 1203 nt = 0 204 DO WHILE ( time_vert(nt) < end_time ) 205 nt = nt + 1 176 206 hash = "#" 177 207 ierrn = 1 ! not zero … … 180 210 !-- from there onwards the profiles will be read 181 211 DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 182 READ ( finput, *, IOSTAT=ierrn ) hash, time_vert( t)212 READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt) 183 213 IF ( ierrn < 0 ) THEN 184 214 WRITE( message_string, * ) 'No time dependent vertical profiles',& … … 188 218 ENDDO 189 219 190 IF ( t == 1 .AND. time_vert(t) > end_time ) EXIT 191 192 READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert, & 193 lowwsubs_vert 220 IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT 221 222 READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert, & 223 lowwsubs_vert, lowpt_lsa, & 224 lowq_lsa, lowpt_subs, lowq_subs 194 225 IF ( ierrn /= 0 ) THEN 195 226 message_string = 'errors in file LSF_DATA' … … 197 228 ENDIF 198 229 199 READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, highvg_vert, & 200 highwsubs_vert 230 READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, & 231 highvg_vert, highwsubs_vert, & 232 highpt_lsa, highq_lsa, & 233 highpt_subs, highq_subs 201 234 202 235 IF ( ierrn /= 0 ) THEN … … 212 245 lowvg_vert = highvg_vert 213 246 lowwsubs_vert = highwsubs_vert 247 lowpt_lsa = highpt_lsa 248 lowq_lsa = highq_lsa 249 lowpt_subs = highpt_subs 250 lowq_subs = highq_subs 214 251 215 252 ierrn = 0 216 READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, & 217 highvg_vert, highwsubs_vert 253 READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, & 254 highvg_vert, highwsubs_vert, & 255 highpt_lsa, highq_lsa, & 256 highpt_subs, highq_subs 218 257 219 258 IF ( ierrn /= 0 ) THEN 220 message_string = 'errors in file LSF_DATA' 221 CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 ) 259 WRITE( message_string, * ) 'zu(nzt+1) = ', zu(nzt+1), 'm ', & 260 'is higher than the maximum height in LSF_DATA which ',& 261 'is ', lowheight, 'm. Interpolation on PALM ', & 262 'grid is not possible.' 263 CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 ) 222 264 ENDIF 223 265 … … 228 270 fac = (highheight-zu(k))/(highheight - lowheight) 229 271 230 ug_vert(k,t) = fac*lowug_vert + (1-fac)*highug_vert 231 vg_vert(k,t) = fac*lowvg_vert + (1-fac)*highvg_vert 232 wsubs_vert(k,t) = fac*lowwsubs_vert + (1-fac)*highwsubs_vert 272 ug_vert(k,nt) = fac*lowug_vert + (1.0_wp-fac)*highug_vert 273 vg_vert(k,nt) = fac*lowvg_vert + (1.0_wp-fac)*highvg_vert 274 wsubs_vert(k,nt) = fac*lowwsubs_vert + (1.0_wp-fac)*highwsubs_vert 275 276 pt_lsa(k,nt) = fac*lowpt_lsa + (1.0_wp-fac)*highpt_lsa 277 q_lsa(k,nt) = fac*lowq_lsa + (1.0_wp-fac)*highq_lsa 278 pt_subs(k,nt) = fac*lowpt_subs + (1.0_wp-fac)*highpt_subs 279 q_subs(k,nt) = fac*lowq_subs + (1.0_wp-fac)*highq_subs 233 280 234 281 ENDDO … … 241 288 242 289 IF ( time_vert(1) > end_time ) THEN 243 WRITE ( message_string, * ) 'Time dependent large scale profile ', &244 'forcing from&LSF_DATA sets in after end of ' , &290 WRITE ( message_string, * ) 'Time dependent large scale profile ', & 291 'forcing from&LSF_DATA sets in after end of ' , & 245 292 'simulation - lsf_vert is set to FALSE' 246 293 CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 ) … … 266 313 USE kinds 267 314 268 269 315 IMPLICIT NONE 270 316 271 INTEGER(iwp) :: t !:317 INTEGER(iwp) :: nt !: 272 318 273 319 REAL(wp) :: fac !: … … 276 322 ! 277 323 !-- Interpolation in time of LSF_DATA at the surface 278 t = 1279 DO WHILE ( time > time_surf( t) )280 t =t + 1324 nt = 1 325 DO WHILE ( time > time_surf(nt) ) 326 nt = nt + 1 281 327 ENDDO 282 IF ( time /= time_surf( t) ) THEN283 t =t - 1284 ENDIF 285 286 fac = ( time -time_surf( t) ) / ( time_surf(t+1) - time_surf(t) )328 IF ( time /= time_surf(nt) ) THEN 329 nt = nt - 1 330 ENDIF 331 332 fac = ( time -time_surf(nt) ) / ( time_surf(nt+1) - time_surf(nt) ) 287 333 288 334 IF ( ibc_pt_b == 0 ) THEN … … 290 336 !-- In case of Dirichlet boundary condition shf must not 291 337 !-- be set - it is calculated via MOST in prandtl_fluxes 292 pt_surface = pt_surf( t) + fac * ( pt_surf(t+1) - pt_surf(t) )338 pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) ) 293 339 294 340 ELSEIF ( ibc_pt_b == 1 ) THEN … … 296 342 !-- In case of Neumann boundary condition pt_surface is needed for 297 343 !-- calculation of reference density 298 shf = shf_surf( t) + fac * ( shf_surf(t+1) - shf_surf(t) )299 pt_surface = pt_surf( t) + fac * ( pt_surf(t+1) - pt_surf(t) )344 shf = shf_surf(nt) + fac * ( shf_surf(nt+1) - shf_surf(nt) ) 345 pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) ) 300 346 301 347 ENDIF … … 305 351 !-- In case of Dirichlet boundary condition qsws must not 306 352 !-- be set - it is calculated via MOST in prandtl_fluxes 307 q_surface = q_surf( t) + fac * ( q_surf(t+1) - q_surf(t) )353 q_surface = q_surf(nt) + fac * ( q_surf(nt+1) - q_surf(nt) ) 308 354 309 355 ELSEIF ( ibc_q_b == 1 ) THEN 310 356 311 qsws = qsws_surf( t) + fac * ( qsws_surf(t+1) - qsws_surf(t) )312 313 ENDIF 314 315 surface_pressure = p_surf( t) + fac * ( p_surf(t+1) - p_surf(t) )357 qsws = qsws_surf(nt) + fac * ( qsws_surf(nt+1) - qsws_surf(nt) ) 358 359 ENDIF 360 361 surface_pressure = p_surf(nt) + fac * ( p_surf(nt+1) - p_surf(nt) ) 316 362 317 363 END SUBROUTINE ls_forcing_surf … … 328 374 USE kinds 329 375 330 331 376 IMPLICIT NONE 332 377 333 INTEGER(iwp) :: t !:378 INTEGER(iwp) :: nt !: 334 379 335 380 REAL(wp) :: fac !: … … 338 383 ! 339 384 !-- Interpolation in time of LSF_DATA for ug, vg and w_subs 340 t = 1341 DO WHILE ( time > time_vert( t) )342 t =t + 1385 nt = 1 386 DO WHILE ( time > time_vert(nt) ) 387 nt = nt + 1 343 388 ENDDO 344 IF ( time /= time_vert( t) ) THEN345 t =t - 1346 ENDIF 347 348 fac = ( time-time_vert( t) ) / ( time_vert(t+1)-time_vert(t) )349 350 ug = ug_vert(:, t) + fac * ( ug_vert(:,t+1) - ug_vert(:,t) )351 vg = vg_vert(:, t) + fac * ( vg_vert(:,t+1) - vg_vert(:,t) )389 IF ( time /= time_vert(nt) ) THEN 390 nt = nt - 1 391 ENDIF 392 393 fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) ) 394 395 ug = ug_vert(:,nt) + fac * ( ug_vert(:,nt+1) - ug_vert(:,nt) ) 396 vg = vg_vert(:,nt) + fac * ( vg_vert(:,nt+1) - vg_vert(:,nt) ) 352 397 353 398 IF ( large_scale_subsidence ) THEN 354 w_subs = wsubs_vert(:,t) + fac * ( wsubs_vert(:,t+1) - wsubs_vert(:,t) ) 399 w_subs = wsubs_vert(:,nt) & 400 + fac * ( wsubs_vert(:,nt+1) - wsubs_vert(:,nt) ) 355 401 ENDIF 356 402 … … 358 404 359 405 406 !------------------------------------------------------------------------------! 407 ! Call for all grid points 408 !------------------------------------------------------------------------------! 409 SUBROUTINE ls_advec ( time, prog_var ) 410 411 USE arrays_3d, & 412 ONLY: pt_lsa, pt_subs, q_lsa, q_subs, tend, time_vert 413 414 USE control_parameters, & 415 ONLY: large_scale_subsidence, use_subsidence_tendencies 416 417 USE indices 418 419 USE kinds 420 421 IMPLICIT NONE 422 423 CHARACTER (LEN=*) :: prog_var !: 424 425 REAL(wp), INTENT(in) :: time !: 426 REAL(wp) :: fac !: 427 428 INTEGER(iwp) :: i !: 429 INTEGER(iwp) :: j !: 430 INTEGER(iwp) :: k !: 431 INTEGER(iwp) :: nt !: 432 433 ! 434 !-- Interpolation in time of LSF_DATA 435 nt = 1 436 DO WHILE ( time > time_vert(nt) ) 437 nt = nt + 1 438 ENDDO 439 IF ( time /= time_vert(nt) ) THEN 440 nt = nt - 1 441 ENDIF 442 443 fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) ) 444 445 ! 446 !-- Add horizontal large scale advection tendencies of pt and q 447 SELECT CASE ( prog_var ) 448 449 CASE ( 'pt' ) 450 451 DO i = nxl, nxr 452 DO j = nys, nyn 453 DO k = nzb_u_inner(j,i)+1, nzt 454 tend(k,j,i) = tend(k,j,i) + pt_lsa(k,nt) & 455 + fac * ( pt_lsa(k,nt+1) - pt_lsa(k,nt) ) 456 ENDDO 457 ENDDO 458 ENDDO 459 460 CASE ( 'q' ) 461 462 DO i = nxl, nxr 463 DO j = nys, nyn 464 DO k = nzb_u_inner(j,i)+1, nzt 465 tend(k,j,i) = tend(k,j,i) + q_lsa(k,nt) & 466 + fac * ( q_lsa(k,nt+1) - q_lsa(k,nt) ) 467 ENDDO 468 ENDDO 469 ENDDO 470 471 END SELECT 472 473 ! 474 !-- Subsidence of pt and q with prescribed subsidence tendencies 475 IF ( large_scale_subsidence .AND. use_subsidence_tendencies ) THEN 476 477 SELECT CASE ( prog_var ) 478 479 CASE ( 'pt' ) 480 481 DO i = nxl, nxr 482 DO j = nys, nyn 483 DO k = nzb_u_inner(j,i)+1, nzt 484 tend(k,j,i) = tend(k,j,i) + pt_subs(k,nt) + fac * & 485 ( pt_subs(k,nt+1) - pt_subs(k,nt) ) 486 ENDDO 487 ENDDO 488 ENDDO 489 490 CASE ( 'q' ) 491 492 DO i = nxl, nxr 493 DO j = nys, nyn 494 DO k = nzb_u_inner(j,i)+1, nzt 495 tend(k,j,i) = tend(k,j,i) + q_subs(k,nt) + fac * & 496 ( q_subs(k,nt+1) - q_subs(k,nt) ) 497 ENDDO 498 ENDDO 499 ENDDO 500 501 END SELECT 502 503 ENDIF 504 505 END SUBROUTINE ls_advec 506 507 508 !------------------------------------------------------------------------------! 509 ! Call for grid point i,j 510 !------------------------------------------------------------------------------! 511 SUBROUTINE ls_advec_ij ( i, j, time, prog_var ) 512 513 USE arrays_3d, & 514 ONLY: pt_lsa, pt_subs, q_lsa, q_subs, tend, time_vert 515 516 USE control_parameters, & 517 ONLY: large_scale_subsidence, use_subsidence_tendencies 518 519 USE indices 520 521 USE kinds 522 523 IMPLICIT NONE 524 525 CHARACTER (LEN=*) :: prog_var !: 526 527 REAL(wp), INTENT(in) :: time !: 528 REAL(wp) :: fac !: 529 530 INTEGER(iwp) :: i !: 531 INTEGER(iwp) :: j !: 532 INTEGER(iwp) :: k !: 533 INTEGER(iwp) :: nt !: 534 535 ! 536 !-- Interpolation in time of LSF_DATA 537 nt = 1 538 DO WHILE ( time > time_vert(nt) ) 539 nt = nt + 1 540 ENDDO 541 IF ( time /= time_vert(nt) ) THEN 542 nt = nt - 1 543 ENDIF 544 545 fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) ) 546 547 ! 548 !-- Add horizontal large scale advection tendencies of pt and q 549 SELECT CASE ( prog_var ) 550 551 CASE ( 'pt' ) 552 553 DO k = nzb_u_inner(j,i)+1, nzt 554 tend(k,j,i) = tend(k,j,i) + pt_lsa(k,nt) & 555 + fac * ( pt_lsa(k,nt+1) - pt_lsa(k,nt) ) 556 ENDDO 557 558 CASE ( 'q' ) 559 560 DO k = nzb_u_inner(j,i)+1, nzt 561 tend(k,j,i) = tend(k,j,i) + q_lsa(k,nt) & 562 + fac * ( q_lsa(k,nt+1) - q_lsa(k,nt) ) 563 ENDDO 564 565 END SELECT 566 567 ! 568 !-- Subsidence of pt and q with prescribed profiles 569 IF ( large_scale_subsidence .AND. use_subsidence_tendencies ) THEN 570 571 SELECT CASE ( prog_var ) 572 573 CASE ( 'pt' ) 574 575 DO k = nzb_u_inner(j,i)+1, nzt 576 tend(k,j,i) = tend(k,j,i) + pt_subs(k,nt) + fac * & 577 ( pt_subs(k,nt+1) - pt_subs(k,nt) ) 578 ENDDO 579 580 CASE ( 'q' ) 581 582 DO k = nzb_u_inner(j,i)+1, nzt 583 tend(k,j,i) = tend(k,j,i) + q_subs(k,nt) + fac * & 584 ( q_subs(k,nt+1) - q_subs(k,nt) ) 585 ENDDO 586 587 END SELECT 588 589 ENDIF 590 591 END SUBROUTINE ls_advec_ij 592 593 360 594 END MODULE ls_forcing_mod
Note: See TracChangeset
for help on using the changeset viewer.
