Changeset 3395 for palm/trunk/UTIL/inifor/src/transform.f90
- Timestamp:
- Oct 22, 2018 5:32:49 PM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/UTIL/inifor/src/transform.f90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/UTIL/inifor/src/transform.f90
r3183 r3395 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Switched addressing of averaging regions from index bounds to list of columns 29 ! Added routines for the computation of geostrophic winds including: 30 ! - the downward extrapolation of density (linear) and 31 ! - pressure (hydrostatic equation) and 32 ! - rotation of geostrophic wind components to PALM frame of reference 33 ! 34 ! 3183 2018-07-27 14:25:55Z suehring 28 35 ! Introduced new PALM grid stretching 29 36 ! Removed unnecessary subroutine parameters … … 51 58 USE control 52 59 USE defs, & 53 ONLY: TO_DEGREES, TO_RADIANS, PI, dp60 ONLY: G, TO_DEGREES, TO_RADIANS, PI, dp 54 61 USE types 55 62 USE util, & … … 176 183 177 184 178 SUBROUTINE average_2d(in_arr, out_arr, i min, imax, jmin, jmax)179 REAL(dp), INTENT(IN) :: in_arr(0:,0:,0:)180 REAL(dp), INTENT(OUT) :: out_arr(0:)181 INTEGER, INTENT(IN) :: imin, imax, jmin, jmax182 183 INTEGER :: i, j, k 185 SUBROUTINE average_2d(in_arr, out_arr, ii, jj) 186 REAL(dp), INTENT(IN) :: in_arr(0:,0:,0:) 187 REAL(dp), INTENT(OUT) :: out_arr(0:) 188 INTEGER, INTENT(IN), DIMENSION(:) :: ii, jj 189 190 INTEGER :: i, j, k, l 184 191 REAL(dp) :: ni 185 186 IF (imin < 0) CALL abort('average_2d', "imin < 0.") 187 IF (jmin < 0) CALL abort('average_2d', "jmin < 0.") 188 IF (imax > UBOUND(in_arr, 1)) CALL abort('average_2d', "imax out of i bound.") 189 IF (jmax > UBOUND(in_arr, 2)) CALL abort('average_2d', "jmax out of j bound.") 192 193 IF (SIZE(ii) .NE. SIZE(jj)) THEN 194 message = "Length of 'ii' and 'jj' index lists do not match." // & 195 NEW_LINE(' ') // "ii has " // str(SIZE(ii)) // " elements, " // & 196 NEW_LINE(' ') // "jj has " // str(SIZE(jj)) // "." 197 CALL abort('average_2d', message) 198 END IF 190 199 191 200 DO k = 0, UBOUND(out_arr, 1) 192 201 193 202 out_arr(k) = 0.0_dp 194 DO j = jmin, jmax 195 DO i = imin, imax 196 out_arr(k) = out_arr(k) + in_arr(i, j, k) 203 DO l = 1, UBOUND(ii, 1) 204 i = ii(l) 205 j = jj(l) 206 out_arr(k) = out_arr(k) +& 207 in_arr(i, j, k) 197 208 END DO 198 END DO 199 200 END DO 201 202 ! devide by number of grid points 203 ni = 1.0_dp / ( (imax - imin + 1) * (jmax - jmin + 1) ) 209 210 END DO 211 212 ni = 1.0_dp / SIZE(ii) 204 213 out_arr(:) = out_arr(:) * ni 205 214 … … 212 221 REAL(dp), DIMENSION(:,:,:), INTENT(OUT) :: palm_array 213 222 REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: intermediate_array 214 INTEGER :: nx, ny, n z223 INTEGER :: nx, ny, nlev 215 224 216 225 nx = palm_intermediate % nx 217 226 ny = palm_intermediate % ny 218 n z= palm_intermediate % nz ! nlev227 nlev = palm_intermediate % nz ! nlev 219 228 220 229 ! Interpolate from COSMO-DE to intermediate grid. Allocating with one 221 230 ! less point in the vertical, since scalars like T have 50 instead of 51 222 231 ! points in COSMO-DE. 223 ALLOCATE(intermediate_array(0:nx, 0:ny, 0:n z-1)) !232 ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nlev-1)) ! 224 233 225 234 CALL interpolate_2d(source_array, intermediate_array, palm_intermediate) … … 234 243 235 244 236 SUBROUTINE average_profile(source_array, profile_array, imin, imax, jmin, jmax,& 237 palm_intermediate, palm_grid) 238 TYPE(grid_definition), INTENT(IN) :: palm_intermediate, palm_grid 245 SUBROUTINE average_profile(source_array, profile_array, avg_grid) 246 TYPE(grid_definition), INTENT(IN) :: avg_grid 239 247 REAL(dp), DIMENSION(:,:,:), INTENT(IN) :: source_array 240 INTEGER, INTENT(IN) :: imin, imax, jmin, jmax 241 REAL(dp), DIMENSION(0:,0:,0:), INTENT(OUT) :: profile_array 242 REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: intermediate_array 243 INTEGER :: nx, ny, nz 244 245 nx = palm_intermediate % nx 246 ny = palm_intermediate % ny 247 nz = palm_intermediate % nz 248 ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nz-1)) 249 intermediate_array(:,:,:) = 0.0_dp 250 251 ! average input array to intermediate profile 252 CALL average_2d(source_array, intermediate_array(0,0,:), imin, imax, jmin, jmax) 253 254 ! vertically interpolate to ouput array 255 CALL interpolate_1d(intermediate_array, profile_array, palm_grid) 256 257 DEALLOCATE(intermediate_array) 248 REAL(dp), DIMENSION(:), INTENT(OUT) :: profile_array 249 250 INTEGER :: i_source, j_source, k_profile, k_source, l, m 251 252 REAL :: ni_columns 253 254 profile_array(:) = 0.0_dp 255 256 DO l = 1, avg_grid % n_columns 257 i_source = avg_grid % iii(l) 258 j_source = avg_grid % jjj(l) 259 260 DO k_profile = avg_grid % k_min, UBOUND(profile_array, 1) ! PALM levels 261 262 DO m = 1, 2 ! vertical interpolation neighbours 263 264 k_source = avg_grid % kkk(l, k_profile, m) 265 266 profile_array(k_profile) = profile_array(k_profile) & 267 + avg_grid % w(l, k_profile, m) & 268 * source_array(i_source, j_source, k_source) 269 270 END DO ! m, vertical interpolation neighbours 271 272 END DO ! k_profile, PALM levels 273 274 END DO ! l, horizontal neighbours 275 276 ni_columns = 1.0_dp / avg_grid % n_columns 277 profile_array(:) = profile_array(:) * ni_columns 278 279 ! Extrapolate constant to the bottom 280 profile_array(1:avg_grid % k_min-1) = profile_array(avg_grid % k_min) 258 281 259 282 END SUBROUTINE average_profile 260 283 284 285 SUBROUTINE extrapolate_density(rho, avg_grid) 286 REAL(dp), DIMENSION(:), INTENT(INOUT) :: rho 287 TYPE(grid_definition), INTENT(IN) :: avg_grid 288 289 REAL(dp) :: drhodz, dz, zk, rhok 290 INTEGER :: k_min 291 292 k_min = avg_grid % k_min 293 zk = avg_grid % z(k_min) 294 rhok = rho(k_min) 295 dz = avg_grid % z(k_min + 1) - avg_grid % z(k_min) 296 drhodz = (rho(k_min + 1) - rho(k_min)) / dz 297 298 rho(1:k_min-1) = rhok + drhodz * (avg_grid % z(1:k_min-1) - zk) 299 300 END SUBROUTINE extrapolate_density 301 302 303 SUBROUTINE extrapolate_pressure(p, rho, avg_grid) 304 REAL(dp), DIMENSION(:), INTENT(IN) :: rho 305 REAL(dp), DIMENSION(:), INTENT(INOUT) :: p 306 TYPE(grid_definition), INTENT(IN) :: avg_grid 307 308 REAL(dp) :: drhodz, dz, zk, rhok 309 INTEGER :: k, k_min 310 311 k_min = avg_grid % k_min 312 zk = avg_grid % z(k_min) 313 rhok = rho(k_min) 314 dz = avg_grid % z(k_min + 1) - avg_grid % z(k_min) 315 drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz 316 317 DO k = 1, k_min-1 318 p(k) = constant_density_pressure(p(k_min), zk, rhok, drhodz, & 319 avg_grid % z(k), G) 320 END DO 321 322 END SUBROUTINE extrapolate_pressure 323 324 325 !------------------------------------------------------------------------------! 326 ! Description: 327 ! ------------ 328 !> Takes the averaged pressure profile <p> and sets the lowest entry to the 329 !> extrapolated pressure at the surface. 330 !------------------------------------------------------------------------------! 331 SUBROUTINE get_surface_pressure(p, rho, avg_grid) 332 REAL(dp), DIMENSION(:), INTENT(IN) :: rho 333 REAL(dp), DIMENSION(:), INTENT(INOUT) :: p 334 TYPE(grid_definition), INTENT(IN) :: avg_grid 335 336 REAL(dp) :: drhodz, dz, zk, rhok 337 INTEGER :: k, k_min 338 339 k_min = avg_grid % k_min 340 zk = avg_grid % z(k_min) 341 rhok = rho(k_min) 342 dz = avg_grid % z(k_min + 1) - avg_grid % z(k_min) 343 drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz 344 345 p(1) = constant_density_pressure(p(k_min), zk, rhok, drhodz, & 346 0.0, G) 347 348 END SUBROUTINE get_surface_pressure 349 350 351 FUNCTION constant_density_pressure(pk, zk, rhok, drhodz, z, g) RESULT(p) 352 353 REAL(dp), INTENT(IN) :: pk, zk, rhok, drhodz, g, z 354 REAL(dp) :: p 355 356 p = pk + ( zk - z ) * g * ( rhok + 0.5*drhodz * (zk - z) ) 357 358 END FUNCTION constant_density_pressure 359 360 !-----------------------------------------------------------------------------! 361 ! Description: 362 ! ----------- 363 !> This routine computes profiles of the two geostrophic wind components ug and 364 !> vg. 365 !-----------------------------------------------------------------------------! 366 SUBROUTINE geostrophic_winds(p_north, p_south, p_east, p_west, rho, f3, & 367 Lx, Ly, phi_n, lam_n, phi_g, lam_g, ug, vg) 368 369 REAL(dp), DIMENSION(:), INTENT(IN) :: p_north, p_south, p_east, p_west, & 370 rho 371 REAL(dp), INTENT(IN) :: f3, Lx, Ly, phi_n, lam_n, phi_g, lam_g 372 REAL(dp), DIMENSION(:), INTENT(OUT) :: ug, vg 373 REAL(dp) :: facx, facy 374 375 facx = 1.0_dp / (Lx * f3) 376 facy = 1.0_dp / (Ly * f3) 377 ug(:) = - facy / rho(:) * (p_north(:) - p_south(:)) 378 vg(:) = facx / rho(:) * (p_east(:) - p_west(:)) 379 380 CALL rotate_vector_field( & 381 ug, vg, angle = meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g)& 382 ) 383 384 END SUBROUTINE geostrophic_winds 261 385 262 386 … … 419 543 420 544 545 SUBROUTINE rotate_vector_field(x, y, angle) 546 REAL(dp), DIMENSION(:), INTENT(INOUT) :: x, y !< x and y coodrinate in arbitrary units 547 REAL(dp), INTENT(IN) :: angle !< rotation angle [deg] 548 549 INTEGER :: i 550 REAL(dp) :: sine, cosine, v_rot(2), rotation(2,2) 551 552 sine = SIN(angle * TO_RADIANS) 553 cosine = COS(angle * TO_RADIANS) 554 ! RESAHPE() fills columns first, so the rotation matrix becomes 555 ! 556 ! rotation = [ cosine -sine ] 557 ! [ sine cosine ] 558 rotation = RESHAPE( (/cosine, sine, -sine, cosine/), (/2, 2/) ) 559 560 DO i = LBOUND(x, 1), UBOUND(x, 1) 561 562 v_rot(:) = MATMUL(rotation, (/x(i), y(i)/)) 563 564 x(i) = v_rot(1) 565 y(i) = v_rot(2) 566 567 END DO 568 569 END SUBROUTINE rotate_vector_field 570 571 572 573 !------------------------------------------------------------------------------! 574 ! Description: 575 ! ------------ 576 !> This routine computes the local meridian convergence between a rotated-pole 577 !> and a geographical system using the Eq. (6) given in the DWD manual 578 !> 579 !> Baldauf et al. (2018), Beschreibung des operationelle KuÌrzestfrist- 580 !> vorhersagemodells COSMO-D2 und COSMO-D2-EPS und seiner Ausgabe in die 581 !> Datenbanken des DWD. 582 !> https://www.dwd.de/SharedDocs/downloads/DE/modelldokumentationen/nwv/cosmo_d2/cosmo_d2_dbbeschr_aktuell.pdf?__blob=publicationFile&v=2 583 !> 584 FUNCTION meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g) & 585 RESULT(delta) 586 587 REAL(dp), INTENT(IN) :: phi_n, lam_n, phi_g, lam_g 588 REAL(dp) :: delta 589 590 delta = atan2( COS(phi_n) * SIN(lam_n - lam_g), & 591 COS(phi_g) * SIN(phi_n) - & 592 SIN(phi_g) * COS(phi_n) * COS(lam_n - lam_g) ) 593 594 END FUNCTION meridian_convergence_rotated 421 595 422 596 !------------------------------------------------------------------------------! … … 512 686 513 687 514 SUBROUTINE find_vertical_neighbours_and_weights(palm_grid, palm_intermediate) 688 SUBROUTINE find_vertical_neighbours_and_weights_interp( palm_grid, & 689 palm_intermediate ) 515 690 TYPE(grid_definition), INTENT(INOUT) :: palm_grid 516 691 TYPE(grid_definition), INTENT(IN) :: palm_intermediate … … 528 703 529 704 ! in each column of the fine grid, find vertical neighbours of every cell 705 DO j = 0, ny 530 706 DO i = 0, nx 531 DO j = 0, ny532 707 533 708 k_intermediate = 0 … … 536 711 column_top = palm_intermediate % h(i,j,nlev) 537 712 538 ! scan through palm_grid column untiland set neighbour indices in713 ! scan through palm_grid column and set neighbour indices in 539 714 ! case current_height is either below column_base, in the current 540 715 ! cell, or above column_top. Keep increasing current cell index until … … 580 755 h_top = palm_intermediate % h(i,j,k_intermediate+1) 581 756 h_bottom = palm_intermediate % h(i,j,k_intermediate) 757 point_is_in_current_cell = ( & 758 current_height >= h_bottom .AND. & 759 current_height < h_top & 760 ) 761 END DO 762 763 IF (k_intermediate > nlev-1) THEN 764 message = "Index " // TRIM(str(k_intermediate)) // & 765 " is above intermediate grid range." 766 CALL abort('find_vertical_neighbours', message) 767 END IF 768 769 palm_grid % kk(i,j,k,1) = k_intermediate 770 palm_grid % kk(i,j,k,2) = k_intermediate + 1 771 772 ! copmute vertical weights 773 weight = (h_top - current_height) / (h_top - h_bottom) 774 palm_grid % w_verti(i,j,k,1) = weight 775 palm_grid % w_verti(i,j,k,2) = 1.0_dp - weight 776 END IF 777 778 END DO 779 780 END DO 781 END DO 782 783 END SUBROUTINE find_vertical_neighbours_and_weights_interp 784 785 786 SUBROUTINE find_vertical_neighbours_and_weights_average( avg_grid ) 787 TYPE(grid_definition), INTENT(INOUT) :: avg_grid 788 789 INTEGER :: i, j, k_palm, k_intermediate, l, nlev 790 LOGICAL :: point_is_below_grid, point_is_above_grid, & 791 point_is_in_current_cell 792 REAL(dp) :: current_height, column_base, column_top, h_top, h_bottom, & 793 weight 794 795 796 avg_grid % k_min = LBOUND(avg_grid % z, 1) 797 798 nlev = SIZE(avg_grid % cosmo_h, 3) 799 800 ! in each column of the fine grid, find vertical neighbours of every cell 801 DO l = 1, avg_grid % n_columns 802 803 i = avg_grid % iii(l) 804 j = avg_grid % jjj(l) 805 806 column_base = avg_grid % cosmo_h(i,j,1) 807 column_top = avg_grid % cosmo_h(i,j,nlev) 808 809 ! scan through avg_grid column until and set neighbour indices in 810 ! case current_height is either below column_base, in the current 811 ! cell, or above column_top. Keep increasing current cell index until 812 ! the current cell overlaps with the current_height. 813 k_intermediate = 1 !avg_grid % cosmo_h is indezed 1-based. 814 DO k_palm = 1, avg_grid % nz 815 816 ! Memorize the top and bottom boundaries of the coarse cell and the 817 ! current height within it 818 current_height = avg_grid % z(k_palm) + avg_grid % z0 819 h_top = avg_grid % cosmo_h(i,j,k_intermediate+1) 820 h_bottom = avg_grid % cosmo_h(i,j,k_intermediate) 821 822 point_is_above_grid = (current_height > column_top) !22000m, very unlikely 823 point_is_below_grid = (current_height < column_base) 824 825 point_is_in_current_cell = ( & 826 current_height >= h_bottom .AND. & 827 current_height < h_top & 828 ) 829 830 ! set default weights 831 avg_grid % w(l,k_palm,1:2) = 0.0_dp 832 833 IF (point_is_above_grid) THEN 834 835 avg_grid % kkk(l,k_palm,1:2) = nlev 836 avg_grid % w(l,k_palm,1:2) = - 2.0_dp 837 838 message = "PALM-4U grid extends above COSMO-DE model top." 839 CALL abort('find_vertical_neighbours_and_weights_average', message) 840 841 ELSE IF (point_is_below_grid) THEN 842 843 avg_grid % kkk(l,k_palm,1:2) = 0 844 avg_grid % w(l,k_palm,1:2) = - 2.0_dp 845 avg_grid % k_min = MAX(k_palm + 1, avg_grid % k_min) 846 ELSE 847 ! cycle through intermediate levels until current 848 ! intermediate-grid cell overlaps with current_height 849 DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1) 850 k_intermediate = k_intermediate + 1 851 852 h_top = avg_grid % cosmo_h(i,j,k_intermediate+1) 853 h_bottom = avg_grid % cosmo_h(i,j,k_intermediate) 582 854 point_is_in_current_cell = ( & 583 855 current_height >= h_bottom .AND. & … … 594 866 END IF 595 867 596 palm_grid % kk(i,j,k,1) = k_intermediate597 palm_grid % kk(i,j,k,2) = k_intermediate + 1868 avg_grid % kkk(l,k_palm,1) = k_intermediate 869 avg_grid % kkk(l,k_palm,2) = k_intermediate + 1 598 870 599 871 ! copmute vertical weights 600 872 weight = (h_top - current_height) / (h_top - h_bottom) 601 palm_grid % w_verti(i,j,k,1) = weight602 palm_grid % w_verti(i,j,k,2) = 1.0_dp - weight873 avg_grid % w(l,k_palm,1) = weight 874 avg_grid % w(l,k_palm,2) = 1.0_dp - weight 603 875 END IF 604 876 605 END DO 606 607 END DO 608 END DO 609 610 END SUBROUTINE find_vertical_neighbours_and_weights 877 END DO ! k, PALM levels 878 END DO ! l, averaging columns 879 880 END SUBROUTINE find_vertical_neighbours_and_weights_average 611 881 612 882 !------------------------------------------------------------------------------!
Note: See TracChangeset
for help on using the changeset viewer.
