Changeset 3182 for palm/trunk/UTIL/inifor/src/transform.f90
- Timestamp:
- Jul 27, 2018 1:36:03 PM (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/UTIL/inifor/src/transform.f90
r2718 r3182 21 21 ! Current revisions: 22 22 ! ----------------- 23 ! Introduced new PALM grid stretching 24 ! Removed unnecessary subroutine parameters 25 ! Renamed kcur to k_intermediate 23 26 ! 24 27 ! … … 80 83 TYPE(grid_definition), INTENT(IN) :: outgrid 81 84 REAL(dp), INTENT(IN) :: in_arr(0:,0:,0:) 82 REAL(dp), INTENT(OUT) :: out_arr(0:,0:,0:) 83 84 INTEGER :: i, j, k, l, nx, ny, nz 85 86 nx = UBOUND(out_arr, 1) 87 ny = UBOUND(out_arr, 2) 85 REAL(dp), INTENT(OUT) :: out_arr(0:,0:,:) 86 87 INTEGER :: i, j, k, l, nz 88 88 89 nz = UBOUND(out_arr, 3) 89 90 90 DO j = 0, ny91 DO i = 0, nx92 DO k = nz, 0, -191 DO j = LBOUND(out_arr, 2), UBOUND(out_arr, 2) 92 DO i = LBOUND(out_arr, 1), UBOUND(out_arr, 1) 93 DO k = nz, LBOUND(out_arr, 3), -1 93 94 94 95 ! TODO: Remove IF clause and extrapolate based on a critical vertical … … 101 102 out_arr(i,j,k) = 0.0_dp 102 103 DO l = 1, 2 103 out_arr(i,j,k) = out_arr(i,j,k) + 104 outgrid % w_verti(i,j,k,l) * 104 out_arr(i,j,k) = out_arr(i,j,k) + & 105 outgrid % w_verti(i,j,k,l) * & 105 106 in_arr(i,j,outgrid % kk(i,j,k, l) ) 106 107 END DO … … 139 140 ! I index 0-based for the indices of the outvar to be consistent with the 140 141 ! outgrid indices and interpolation weights. 141 TYPE(grid_definition), INTENT(IN) :: outgrid142 REAL(dp), INTENT(IN) :: invar(0:,0:,0:)143 REAL(dp), INTENT(OUT) :: outvar(0:,0:,0:)142 TYPE(grid_definition), INTENT(IN) :: outgrid 143 REAL(dp), INTENT(IN) :: invar(0:,0:,0:) 144 REAL(dp), INTENT(OUT) :: outvar(0:,0:,0:) 144 145 TYPE(nc_var), INTENT(IN), OPTIONAL :: ncvar 145 146 … … 413 414 414 415 END SUBROUTINE rotate_to_cosmo 416 415 417 416 418 … … 427 429 !> ------------- 428 430 !> jj, lat 429 !> ^ j430 !> | \ i431 !> ^ j 432 !> | \ i 431 433 !> jj(i,j,2/3) + ... 2 ---\--------/------ 3 432 434 !> | | ^ \ / | … … 459 461 !> 460 462 !------------------------------------------------------------------------------! 461 SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon, cosmo_dxi, & 462 cosmo_dyi, palm_clat, palm_clon, palm_ii, palm_jj) 463 SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon, & 464 palm_clat, palm_clon, & 465 palm_ii, palm_jj) 463 466 464 467 REAL(dp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon 465 468 REAL(dp), DIMENSION(0:,0:), INTENT(IN) :: palm_clat, palm_clon 466 REAL(dp) , INTENT(IN):: cosmo_dxi, cosmo_dyi469 REAL(dp) :: cosmo_dxi, cosmo_dyi 467 470 INTEGER, DIMENSION(0:,0:,1:), INTENT(OUT) :: palm_ii, palm_jj 468 471 … … 472 475 lon0 = cosmo_lon(0) 473 476 lat0 = cosmo_lat(0) 477 cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0)) 478 cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0)) 474 479 475 480 DO j = 0, UBOUND(palm_clon, 2)!palm_grid % ny … … 508 513 TYPE(grid_definition), INTENT(IN) :: palm_intermediate 509 514 510 INTEGER :: i, j, k, nx, ny, nz, nlev, k cur515 INTEGER :: i, j, k, nx, ny, nz, nlev, k_intermediate 511 516 LOGICAL :: point_is_below_grid, point_is_above_grid, & 512 517 point_is_in_current_cell … … 523 528 DO j = 0, ny 524 529 525 k cur= 0530 k_intermediate = 0 526 531 527 532 column_base = palm_intermediate % h(i,j,0) … … 532 537 ! cell, or above column_top. Keep increasing current cell index until 533 538 ! the current cell overlaps with the current_height. 534 DO k = 0, nz539 DO k = 1, nz 535 540 536 541 ! Memorize the top and bottom boundaries of the coarse cell and the 537 542 ! current height within it 538 543 current_height = palm_grid % z(k) + palm_grid % z0 539 h_top = palm_intermediate % h(i,j,k cur+1)540 h_bottom = palm_intermediate % h(i,j,k cur)544 h_top = palm_intermediate % h(i,j,k_intermediate+1) 545 h_bottom = palm_intermediate % h(i,j,k_intermediate) 541 546 542 547 point_is_above_grid = (current_height > column_top) !22000m, very unlikely … … 556 561 palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp 557 562 563 message = "PALM-4U grid extends above COSMO-DE model top." 564 CALL abort('find_vertical_neighbours_and_weights', message) 565 558 566 ELSE IF (point_is_below_grid) THEN 559 567 … … 564 572 ! cycle through intermediate levels until current 565 573 ! intermediate-grid cell overlaps with current_height 566 DO WHILE (.NOT. point_is_in_current_cell .AND. k cur<= nlev-1)567 k cur = kcur+ 1568 569 h_top = palm_intermediate % h(i,j,k cur+1)570 h_bottom = palm_intermediate % h(i,j,k cur)574 DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1) 575 k_intermediate = k_intermediate + 1 576 577 h_top = palm_intermediate % h(i,j,k_intermediate+1) 578 h_bottom = palm_intermediate % h(i,j,k_intermediate) 571 579 point_is_in_current_cell = ( & 572 580 current_height >= h_bottom .AND. & … … 575 583 END DO 576 584 577 ! kcur = 48 indicates the last section (indices 48 and 49), i.e. 578 ! kcur = 49 is not the beginning of a valid cell. 579 IF (kcur > nlev-1) THEN 580 message = "Index " // TRIM(str(kcur)) // " is above intermediate grid range." 585 ! k_intermediate = 48 indicates the last section (indices 48 and 49), i.e. 586 ! k_intermediate = 49 is not the beginning of a valid cell. 587 IF (k_intermediate > nlev-1) THEN 588 message = "Index " // TRIM(str(k_intermediate)) // & 589 " is above intermediate grid range." 581 590 CALL abort('find_vertical_neighbours', message) 582 591 END IF 583 592 584 palm_grid % kk(i,j,k,1) = k cur585 palm_grid % kk(i,j,k,2) = k cur+ 1593 palm_grid % kk(i,j,k,1) = k_intermediate 594 palm_grid % kk(i,j,k,2) = k_intermediate + 1 586 595 587 596 ! copmute vertical weights … … 643 652 ! 644 653 SUBROUTINE compute_horizontal_interp_weights(cosmo_lat, cosmo_lon, & 645 cosmo_dxi, cosmo_dyi,palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)654 palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz) 646 655 647 656 REAL(dp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon 648 REAL(dp) , INTENT(IN):: cosmo_dxi, cosmo_dyi657 REAL(dp) :: cosmo_dxi, cosmo_dyi 649 658 REAL(dp), DIMENSION(0:,0:), INTENT(IN) :: palm_clat, palm_clon 650 659 INTEGER, DIMENSION(0:,0:,1:), INTENT(IN) :: palm_ii, palm_jj … … 654 663 REAL(dp) :: wl, wp 655 664 INTEGER :: i, j 665 666 cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0)) 667 cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0)) 656 668 657 669 DO j = 0, UBOUND(palm_clon, 2)
Note: See TracChangeset
for help on using the changeset viewer.