Changeset 3866 for palm/trunk/UTIL/inifor/src/inifor_transform.f90
- Timestamp:
- Apr 5, 2019 2:25:01 PM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/UTIL/inifor/src/inifor_transform.f90 (modified) (36 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/UTIL/inifor/src/inifor_transform.f90
r3785 r3866 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Use PALM's working precision 29 ! Improved coding style 30 ! 31 ! 32 ! 3785 2019-03-06 10:41:14Z eckhard 28 33 ! Remove basic state pressure before computing geostrophic wind 29 34 ! - Introduced new level-based profile averaging routine that does not rely on … … 106 111 USE inifor_control 107 112 USE inifor_defs, & 108 ONLY: BETA, dp, G, P_SL, PI, RD, T_SL, TO_DEGREES, TO_RADIANS113 ONLY: BETA, G, P_SL, PI, RD, T_SL, TO_DEGREES, TO_RADIANS, wp 109 114 USE inifor_types 110 115 USE inifor_util, & … … 116 121 117 122 118 119 120 REAL(dp), INTENT(IN) :: in_arr(:)121 REAL(dp), INTENT(OUT) :: out_arr(:)122 123 124 125 126 127 DOk = nz, LBOUND(out_arr, 1), -1128 129 ! 130 !-- 131 !-- 132 !-- 133 !-- 134 IF (outgrid % w(1,k,1) < -1.0_dp .AND. k < nz) THEN135 136 137 out_arr(k) = 0.0_dp138 DOl = 1, 2139 140 outgrid % w(1,k,l) * in_arr(outgrid %kkk(1,k,l) )141 142 143 144 145 123 SUBROUTINE interpolate_1d(in_arr, out_arr, outgrid) 124 TYPE(grid_definition), INTENT(IN) :: outgrid 125 REAL(wp), INTENT(IN) :: in_arr(:) 126 REAL(wp), INTENT(OUT) :: out_arr(:) 127 128 INTEGER :: k, l, nz 129 130 nz = UBOUND(out_arr, 1) 131 132 DO k = nz, LBOUND(out_arr, 1), -1 133 134 ! 135 !-- TODO: Remove IF clause and extrapolate based on a critical vertical 136 !-- TODO: index marking the lower bound of COSMO-DE data coverage. 137 !-- Check for negative interpolation weights indicating grid points 138 !-- below COSMO-DE domain and extrapolate from the top in such cells. 139 IF (outgrid%w(1,k,1) < -1.0_wp .AND. k < nz) THEN 140 out_arr(k) = out_arr(k+1) 141 ELSE 142 out_arr(k) = 0.0_wp 143 DO l = 1, 2 144 out_arr(k) = out_arr(k) + & 145 outgrid%w(1,k,l) * in_arr(outgrid%kkk(1,k,l) ) 146 ENDDO 147 ENDIF 148 ENDDO 149 150 END SUBROUTINE interpolate_1d 146 151 147 152 … … 158 163 !> invar : Array of source data 159 164 !> 160 !> outgrid %kk : Array of vertical neighbour indices. kk(i,j,k,:) contain the165 !> outgrid%kk : Array of vertical neighbour indices. kk(i,j,k,:) contain the 161 166 !> indices of the two vertical neighbors of PALM-4U point (i,j,k) on the 162 167 !> input grid corresponding to the source data invar. 163 168 !> 164 !> outgrid %w_verti : Array of weights for vertical linear interpolation169 !> outgrid%w_verti : Array of weights for vertical linear interpolation 165 170 !> corresponding to neighbour points indexed by kk. 166 171 !> … … 169 174 !> outvar : Array of interpolated data 170 175 !------------------------------------------------------------------------------! 171 172 173 REAL(dp), INTENT(IN) :: in_arr(0:,0:,0:)174 REAL(dp), INTENT(OUT) :: out_arr(0:,0:,:)175 176 177 178 179 180 DOj = LBOUND(out_arr, 2), UBOUND(out_arr, 2)181 DOi = LBOUND(out_arr, 1), UBOUND(out_arr, 1)182 DOk = nz, LBOUND(out_arr, 3), -1183 184 ! 185 !-- 186 !-- 187 !-- 188 !-- 189 IF (outgrid % w_verti(i,j,k,1) < -1.0_dp .AND. k < nz) THEN190 191 192 out_arr(i,j,k) = 0.0_dp193 DOl = 1, 2194 195 outgrid %w_verti(i,j,k,l) * &196 in_arr(i,j,outgrid %kk(i,j,k, l) )197 198 199 200 201 202 176 SUBROUTINE interpolate_1d_arr(in_arr, out_arr, outgrid) 177 TYPE(grid_definition), INTENT(IN) :: outgrid 178 REAL(wp), INTENT(IN) :: in_arr(0:,0:,0:) 179 REAL(wp), INTENT(OUT) :: out_arr(0:,0:,:) 180 181 INTEGER :: i, j, k, l, nz 182 183 nz = UBOUND(out_arr, 3) 184 185 DO j = LBOUND(out_arr, 2), UBOUND(out_arr, 2) 186 DO i = LBOUND(out_arr, 1), UBOUND(out_arr, 1) 187 DO k = nz, LBOUND(out_arr, 3), -1 188 189 ! 190 !-- TODO: Remove IF clause and extrapolate based on a critical vertical 191 !-- TODO: index marking the lower bound of COSMO-DE data coverage. 192 !-- Check for negative interpolation weights indicating grid points 193 !-- below COSMO-DE domain and extrapolate from the top in such cells. 194 IF (outgrid%w_verti(i,j,k,1) < -1.0_wp .AND. k < nz) THEN 195 out_arr(i,j,k) = out_arr(i,j,k+1) 196 ELSE 197 out_arr(i,j,k) = 0.0_wp 198 DO l = 1, 2 199 out_arr(i,j,k) = out_arr(i,j,k) + & 200 outgrid%w_verti(i,j,k,l) * & 201 in_arr(i,j,outgrid%kk(i,j,k, l) ) 202 ENDDO 203 ENDIF 204 ENDDO 205 ENDDO 206 ENDDO 207 END SUBROUTINE interpolate_1d_arr 203 208 204 209 … … 214 219 !> invar : Array of source data 215 220 !> 216 !> outgrid % ii, %jj : Array of neighbour indices in x and y direction.221 !> outgrid%ii,%jj : Array of neighbour indices in x and y direction. 217 222 !> ii(i,j,k,:), and jj(i,j,k,:) contain the four horizontal neighbour points 218 223 !> of PALM-4U point (i,j,k) on the input grid corresponding to the source … … 220 225 ! form of the interpolation weights.) 221 226 !> 222 !> outgrid %w_horiz: Array of weights for horizontal bi-linear interpolation227 !> outgrid%w_horiz: Array of weights for horizontal bi-linear interpolation 223 228 !> corresponding to neighbour points indexed by ii and jj. 224 229 !> … … 227 232 !> outvar : Array of interpolated data 228 233 !------------------------------------------------------------------------------! 229 SUBROUTINE interpolate_2d(invar, outvar, outgrid, ncvar) 230 ! 231 !-- I index 0-based for the indices of the outvar to be consistent with the 232 !-- outgrid indices and interpolation weights. 233 TYPE(grid_definition), INTENT(IN) :: outgrid 234 REAL(dp), INTENT(IN) :: invar(0:,0:,0:) 235 REAL(dp), INTENT(OUT) :: outvar(0:,0:,0:) 236 TYPE(nc_var), INTENT(IN), OPTIONAL :: ncvar 237 238 INTEGER :: i, j, k, l 239 240 ! 241 !-- TODO: check if input dimensions are consistent, i.e. ranges are correct 242 IF ( UBOUND(outvar, 3) .GT. UBOUND(invar, 3) ) THEN 243 message = "Output array for '" // TRIM(ncvar % name) // "' has ' more levels (" // & 244 TRIM(str(UBOUND(outvar, 3))) // ") than input variable ("//& 245 TRIM(str(UBOUND(invar, 3))) // ")." 246 CALL inifor_abort('interpolate_2d', message) 247 ENDIF 248 249 DO k = 0, UBOUND(outvar, 3) 250 DO j = 0, UBOUND(outvar, 2) 251 DO i = 0, UBOUND(outvar, 1) 252 outvar(i,j,k) = 0.0_dp 253 DO l = 1, 4 254 255 outvar(i,j,k) = outvar(i,j,k) + & 256 outgrid % w_horiz(i,j,l) * invar( outgrid % ii(i,j,l), & 257 outgrid % jj(i,j,l), & 258 k ) 259 ENDDO 234 SUBROUTINE interpolate_2d(invar, outvar, outgrid, ncvar) 235 ! 236 !-- I index 0-based for the indices of the outvar to be consistent with the 237 !-- outgrid indices and interpolation weights. 238 TYPE(grid_definition), INTENT(IN) :: outgrid 239 REAL(wp), INTENT(IN) :: invar(0:,0:,0:) 240 REAL(wp), INTENT(OUT) :: outvar(0:,0:,0:) 241 TYPE(nc_var), INTENT(IN), OPTIONAL :: ncvar 242 243 INTEGER :: i, j, k, l 244 245 ! 246 !-- TODO: check if input dimensions are consistent, i.e. ranges are correct 247 IF ( UBOUND(outvar, 3) .GT. UBOUND(invar, 3) ) THEN 248 message = "Output array for '" // TRIM(ncvar%name) // "' has ' more levels (" // & 249 TRIM(str(UBOUND(outvar, 3))) // ") than input variable ("//& 250 TRIM(str(UBOUND(invar, 3))) // ")." 251 CALL inifor_abort('interpolate_2d', message) 252 ENDIF 253 254 DO k = 0, UBOUND(outvar, 3) 255 DO j = 0, UBOUND(outvar, 2) 256 DO i = 0, UBOUND(outvar, 1) 257 outvar(i,j,k) = 0.0_wp 258 DO l = 1, 4 259 260 outvar(i,j,k) = outvar(i,j,k) + & 261 outgrid%w_horiz(i,j,l) * invar( outgrid%ii(i,j,l), & 262 outgrid%jj(i,j,l), & 263 k ) 260 264 ENDDO 261 ENDDO 262 ENDDO 265 ENDDO 266 ENDDO 267 ENDDO 263 268 264 269 END SUBROUTINE interpolate_2d 265 270 266 271 … … 271 276 !> out_arr(:) 272 277 !------------------------------------------------------------------------------! 273 SUBROUTINE average_2d(in_arr, out_arr, ii, jj) 274 REAL(dp), INTENT(IN) :: in_arr(0:,0:,0:) 275 REAL(dp), INTENT(OUT) :: out_arr(0:) 276 INTEGER, INTENT(IN), DIMENSION(:) :: ii, jj 277 278 INTEGER :: i, j, k, l 279 REAL(dp) :: ni 280 281 IF (SIZE(ii) /= SIZE(jj)) THEN 282 message = "Length of 'ii' and 'jj' index lists do not match." // & 283 NEW_LINE(' ') // "ii has " // str(SIZE(ii)) // " elements, " // & 284 NEW_LINE(' ') // "jj has " // str(SIZE(jj)) // "." 285 CALL inifor_abort('average_2d', message) 286 ENDIF 287 288 IF (SIZE(ii) == 0) THEN 289 message = "No columns to average over; " // & 290 "size of index lists 'ii' and 'jj' is zero." 291 CALL inifor_abort('average_2d', message) 292 ENDIF 293 294 DO k = 0, UBOUND(out_arr, 1) 295 296 out_arr(k) = 0.0_dp 297 DO l = 1, UBOUND(ii, 1) 298 i = ii(l) 299 j = jj(l) 300 out_arr(k) = out_arr(k) + in_arr(i, j, k) 301 ENDDO 302 278 SUBROUTINE average_2d(in_arr, out_arr, ii, jj) 279 REAL(wp), INTENT(IN) :: in_arr(0:,0:,0:) 280 REAL(wp), INTENT(OUT) :: out_arr(0:) 281 INTEGER, INTENT(IN), DIMENSION(:) :: ii, jj 282 283 INTEGER :: i, j, k, l 284 REAL(wp) :: ni 285 286 IF (SIZE(ii) /= SIZE(jj)) THEN 287 message = "Length of 'ii' and 'jj' index lists do not match." // & 288 NEW_LINE(' ') // "ii has " // str(SIZE(ii)) // " elements, " // & 289 NEW_LINE(' ') // "jj has " // str(SIZE(jj)) // "." 290 CALL inifor_abort('average_2d', message) 291 ENDIF 292 293 IF (SIZE(ii) == 0) THEN 294 message = "No columns to average over; " // & 295 "size of index lists 'ii' and 'jj' is zero." 296 CALL inifor_abort('average_2d', message) 297 ENDIF 298 299 DO k = 0, UBOUND(out_arr, 1) 300 301 out_arr(k) = 0.0_wp 302 DO l = 1, UBOUND(ii, 1) 303 i = ii(l) 304 j = jj(l) 305 out_arr(k) = out_arr(k) + in_arr(i, j, k) 303 306 ENDDO 304 307 305 ni = 1.0_dp / SIZE(ii) 306 out_arr(:) = out_arr(:) * ni 307 308 END SUBROUTINE average_2d 308 ENDDO 309 310 ni = 1.0_wp / SIZE(ii) 311 out_arr(:) = out_arr(:) * ni 312 313 END SUBROUTINE average_2d 309 314 310 315 … … 320 325 !> as coarse as COSMO, horizontally as fine as PALM). 321 326 !------------------------------------------------------------------------------! 322 323 324 REAL(dp), DIMENSION(:,:,:), INTENT(IN) :: source_array325 REAL(dp), DIMENSION(:,:,:), INTENT(OUT) :: palm_array326 REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: intermediate_array327 328 329 nx = palm_intermediate %nx330 ny = palm_intermediate %ny331 nlev = palm_intermediate %nz332 333 ! 334 !-- 335 !-- 336 !-- 337 338 339 340 341 ! 342 !-- 343 !-- 344 345 346 347 348 327 SUBROUTINE interpolate_3d(source_array, palm_array, palm_intermediate, palm_grid) 328 TYPE(grid_definition), INTENT(IN) :: palm_intermediate, palm_grid 329 REAL(wp), DIMENSION(:,:,:), INTENT(IN) :: source_array 330 REAL(wp), DIMENSION(:,:,:), INTENT(OUT) :: palm_array 331 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: intermediate_array 332 INTEGER :: nx, ny, nlev 333 334 nx = palm_intermediate%nx 335 ny = palm_intermediate%ny 336 nlev = palm_intermediate%nz 337 338 ! 339 !-- Interpolate from COSMO to intermediate grid. Allocating with one 340 !-- less point in the vertical, since scalars like T have 50 instead of 51 341 !-- points in COSMO. 342 ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nlev-1)) ! 343 344 CALL interpolate_2d(source_array, intermediate_array, palm_intermediate) 345 346 ! 347 !-- Interpolate from intermediate grid to palm_grid grid, includes 348 !-- extrapolation for cells below COSMO domain. 349 CALL interpolate_1d_arr(intermediate_array, palm_array, palm_grid) 350 351 DEALLOCATE(intermediate_array) 352 353 END SUBROUTINE interpolate_3d 349 354 350 355 … … 355 360 !> averaging grid 'avg_grid' and store the result in 'profile_array'. 356 361 !------------------------------------------------------------------------------! 357 SUBROUTINE interp_average_profile(source_array, profile_array, avg_grid) 358 TYPE(grid_definition), INTENT(IN) :: avg_grid 359 REAL(dp), DIMENSION(:,:,:), INTENT(IN) :: source_array 360 REAL(dp), DIMENSION(:), INTENT(OUT) :: profile_array 361 362 INTEGER :: i_source, j_source, k_profile, k_source, l, m 363 364 REAL :: ni_columns 365 366 profile_array(:) = 0.0_dp 367 368 DO l = 1, avg_grid % n_columns 369 i_source = avg_grid % iii(l) 370 j_source = avg_grid % jjj(l) 371 372 ! 373 !-- Loop over PALM levels 374 DO k_profile = avg_grid % k_min, UBOUND(profile_array, 1) 375 376 ! 377 !-- Loop over vertical interpolation neighbours 378 DO m = 1, 2 379 380 k_source = avg_grid % kkk(l, k_profile, m) 381 382 profile_array(k_profile) = profile_array(k_profile) & 383 + avg_grid % w(l, k_profile, m) & 384 * source_array(i_source, j_source, k_source) 385 ! 386 !-- Loop over vertical interpolation neighbours m 387 ENDDO 388 389 ! 390 !-- Loop over PALM levels k_profile 362 SUBROUTINE interp_average_profile(source_array, profile_array, avg_grid) 363 TYPE(grid_definition), INTENT(IN) :: avg_grid 364 REAL(wp), DIMENSION(:,:,:), INTENT(IN) :: source_array 365 REAL(wp), DIMENSION(:), INTENT(OUT) :: profile_array 366 367 INTEGER :: i_source, j_source, k_profile, k_source, l, m 368 369 REAL :: ni_columns 370 371 profile_array(:) = 0.0_wp 372 373 DO l = 1, avg_grid%n_columns 374 i_source = avg_grid%iii(l) 375 j_source = avg_grid%jjj(l) 376 377 ! 378 !-- Loop over PALM levels 379 DO k_profile = avg_grid%k_min, UBOUND(profile_array, 1) 380 381 ! 382 !-- Loop over vertical interpolation neighbours 383 DO m = 1, 2 384 385 k_source = avg_grid%kkk(l, k_profile, m) 386 387 profile_array(k_profile) = profile_array(k_profile) & 388 + avg_grid%w(l, k_profile, m) & 389 * source_array(i_source, j_source, k_source) 390 ! 391 !-- Loop over vertical interpolation neighbours m 391 392 ENDDO 392 393 393 394 ! 394 !-- Loop over horizontal neighbours l395 !-- Loop over PALM levels k_profile 395 396 ENDDO 396 397 397 ni_columns = 1.0_dp / avg_grid % n_columns 398 profile_array(:) = profile_array(:) * ni_columns 399 400 ! 401 !-- Constant extrapolation to the bottom 402 profile_array(1:avg_grid % k_min-1) = profile_array(avg_grid % k_min) 403 404 END SUBROUTINE interp_average_profile 398 ! 399 !-- Loop over horizontal neighbours l 400 ENDDO 401 402 ni_columns = 1.0_wp / avg_grid%n_columns 403 profile_array(:) = profile_array(:) * ni_columns 404 405 ! 406 !-- Constant extrapolation to the bottom 407 profile_array(1:avg_grid%k_min-1) = profile_array(avg_grid%k_min) 408 409 END SUBROUTINE interp_average_profile 405 410 406 411 … … 411 416 !> averaging grid 'avg_grid' and store the result in 'profile_array'. 412 417 !------------------------------------------------------------------------------! 413 414 415 416 REAL(dp), DIMENSION(:,:,:), INTENT(IN) :: source_array417 REAL(dp), DIMENSION(:), INTENT(OUT) :: profile_array418 419 420 421 REAL(dp) :: ni_columns422 423 424 425 426 427 428 429 430 431 432 433 profile_array(:) = 0.0_dp434 435 DO l = 1, avg_grid %n_columns436 437 i_source = avg_grid %iii(l)438 j_source = avg_grid %jjj(l)439 440 441 442 443 444 445 ni_columns = 1.0_dp / avg_grid %n_columns446 447 448 418 SUBROUTINE average_profile( source_array, profile_array, avg_grid ) 419 420 TYPE(grid_definition), INTENT(IN) :: avg_grid 421 REAL(wp), DIMENSION(:,:,:), INTENT(IN) :: source_array 422 REAL(wp), DIMENSION(:), INTENT(OUT) :: profile_array 423 424 INTEGER :: i_source, j_source, l, nz, nlev 425 426 REAL(wp) :: ni_columns 427 428 nlev = SIZE( source_array, 3 ) 429 nz = SIZE( profile_array, 1 ) 430 431 IF ( nlev /= nz ) THEN 432 message = "Lengths of input and output profiles do not match: " // & 433 "cosmo_pressure(" // TRIM( str( nlev ) ) // & 434 "), profile_array(" // TRIM( str( nz ) ) // ")." 435 CALL inifor_abort('average_pressure_perturbation', message) 436 ENDIF 437 438 profile_array(:) = 0.0_wp 439 440 DO l = 1, avg_grid%n_columns 441 442 i_source = avg_grid%iii(l) 443 j_source = avg_grid%jjj(l) 444 445 profile_array(:) = profile_array(:) & 446 + source_array(i_source, j_source, :) 447 448 ENDDO 449 450 ni_columns = 1.0_wp / avg_grid%n_columns 451 profile_array(:) = profile_array(:) * ni_columns 452 453 END SUBROUTINE average_profile 449 454 450 455 … … 456 461 !> averaging. 457 462 !------------------------------------------------------------------------------! 458 459 460 461 462 REAL(dp), DIMENSION(:,:,:), INTENT(IN) :: cosmo_pressure463 REAL(dp), DIMENSION(:), INTENT(OUT) :: profile_array464 465 466 467 REAL(dp) :: ni_columns468 REAL(dp), DIMENSION(:), ALLOCATABLE :: basic_state_pressure469 470 471 472 473 474 475 476 477 478 479 480 481 profile_array(:) = 0.0_dp482 483 DO l = 1, avg_grid %n_columns484 i_source = avg_grid %iii(l)485 j_source = avg_grid %jjj(l)486 487 ! 488 !-- 489 CALL get_basic_state( cosmo_grid %hfl(i_source,j_source,:), BETA, &490 491 492 493 494 495 496 ! 497 !-- 498 499 500 501 502 ni_columns = 1.0_dp / avg_grid %n_columns503 504 505 463 SUBROUTINE average_pressure_perturbation( cosmo_pressure, profile_array, & 464 cosmo_grid, avg_grid ) 465 466 TYPE(grid_definition), INTENT(IN) :: cosmo_grid, avg_grid 467 REAL(wp), DIMENSION(:,:,:), INTENT(IN) :: cosmo_pressure 468 REAL(wp), DIMENSION(:), INTENT(OUT) :: profile_array 469 470 INTEGER :: i_source, j_source, l, nz, nlev 471 472 REAL(wp) :: ni_columns 473 REAL(wp), DIMENSION(:), ALLOCATABLE :: basic_state_pressure 474 475 nlev = SIZE( cosmo_pressure, 3 ) 476 nz = SIZE( profile_array, 1 ) 477 478 IF ( nlev /= nz ) THEN 479 message = "Lengths of input and output profiles do not match: " // & 480 "cosmo_pressure(" // TRIM( str( nlev ) ) // & 481 "), profile_array(" // TRIM( str( nz ) ) // ")." 482 CALL inifor_abort('average_pressure_perturbation', message) 483 ENDIF 484 485 ALLOCATE( basic_state_pressure(nz) ) 486 profile_array(:) = 0.0_wp 487 488 DO l = 1, avg_grid%n_columns 489 i_source = avg_grid%iii(l) 490 j_source = avg_grid%jjj(l) 491 492 ! 493 !-- Compute pressure perturbation by removing COSMO basic state pressure 494 CALL get_basic_state( cosmo_grid%hfl(i_source,j_source,:), BETA, & 495 P_SL, T_SL, RD, G, basic_state_pressure ) 496 497 profile_array(:) = profile_array(:) & 498 + cosmo_pressure(i_source, j_source, :) & 499 - basic_state_pressure(:) 500 501 ! 502 !-- Loop over horizontal neighbours l 503 ENDDO 504 505 DEALLOCATE( basic_state_pressure ) 506 507 ni_columns = 1.0_wp / avg_grid%n_columns 508 profile_array(:) = profile_array(:) * ni_columns 509 510 END SUBROUTINE average_pressure_perturbation 506 511 507 512 … … 513 518 !> Extrapolates density linearly from the level 'k_min' downwards. 514 519 !------------------------------------------------------------------------------! 515 516 REAL(dp), DIMENSION(:), INTENT(INOUT) :: rho517 518 519 REAL(dp) :: drhodz, dz, zk, rhok520 521 522 k_min = avg_grid %k_min523 zk = avg_grid %z(k_min)524 525 dz = avg_grid % z(k_min + 1) - avg_grid %z(k_min)526 527 528 rho(1:k_min-1) = rhok + drhodz * (avg_grid %z(1:k_min-1) - zk)529 530 520 SUBROUTINE extrapolate_density(rho, avg_grid) 521 REAL(wp), DIMENSION(:), INTENT(INOUT) :: rho 522 TYPE(grid_definition), INTENT(IN) :: avg_grid 523 524 REAL(wp) :: drhodz, dz, zk, rhok 525 INTEGER :: k_min 526 527 k_min = avg_grid%k_min 528 zk = avg_grid%z(k_min) 529 rhok = rho(k_min) 530 dz = avg_grid%z(k_min + 1) - avg_grid%z(k_min) 531 drhodz = (rho(k_min + 1) - rho(k_min)) / dz 532 533 rho(1:k_min-1) = rhok + drhodz * (avg_grid%z(1:k_min-1) - zk) 534 535 END SUBROUTINE extrapolate_density 531 536 532 537 … … 536 541 !> Driver for extrapolating pressure from PALM level k_min downwards 537 542 !------------------------------------------------------------------------------! 538 539 REAL(dp), DIMENSION(:), INTENT(IN) :: rho540 REAL(dp), DIMENSION(:), INTENT(INOUT) :: p541 542 543 REAL(dp) :: drhodz, dz, zk, rhok544 545 546 k_min = avg_grid %k_min547 zk = avg_grid %z(k_min)548 549 dz = avg_grid % z(k_min + 1) - avg_grid %z(k_min)550 drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz551 552 DOk = 1, k_min-1553 554 avg_grid %z(k), G)555 556 557 543 SUBROUTINE extrapolate_pressure(p, rho, avg_grid) 544 REAL(wp), DIMENSION(:), INTENT(IN) :: rho 545 REAL(wp), DIMENSION(:), INTENT(INOUT) :: p 546 TYPE(grid_definition), INTENT(IN) :: avg_grid 547 548 REAL(wp) :: drhodz, dz, zk, rhok 549 INTEGER :: k, k_min 550 551 k_min = avg_grid%k_min 552 zk = avg_grid%z(k_min) 553 rhok = rho(k_min) 554 dz = avg_grid%z(k_min + 1) - avg_grid%z(k_min) 555 drhodz = 0.5_wp * (rho(k_min + 1) - rho(k_min)) / dz 556 557 DO k = 1, k_min-1 558 p(k) = constant_density_pressure(p(k_min), zk, rhok, drhodz, & 559 avg_grid%z(k), G) 560 ENDDO 561 562 END SUBROUTINE extrapolate_pressure 558 563 559 564 … … 564 569 !> extrapolated pressure at the surface. 565 570 !------------------------------------------------------------------------------! 566 567 REAL(dp), DIMENSION(:), INTENT(IN) :: rho568 REAL(dp), DIMENSION(:), INTENT(INOUT) :: p569 570 571 REAL(dp) :: drhodz, dz, zk, rhok572 573 574 k_min = avg_grid %k_min575 zk = avg_grid %z(k_min)576 577 dz = avg_grid % z(k_min + 1) - avg_grid %z(k_min)578 drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz579 580 581 0.0_dp, G)582 583 584 585 586 587 588 REAL(dp), INTENT(IN) :: pk, zk, rhok, drhodz, g, z589 REAL(dp) :: p590 591 592 593 571 SUBROUTINE get_surface_pressure(p, rho, avg_grid) 572 REAL(wp), DIMENSION(:), INTENT(IN) :: rho 573 REAL(wp), DIMENSION(:), INTENT(INOUT) :: p 574 TYPE(grid_definition), INTENT(IN) :: avg_grid 575 576 REAL(wp) :: drhodz, dz, zk, rhok 577 INTEGER :: k_min 578 579 k_min = avg_grid%k_min 580 zk = avg_grid%z(k_min) 581 rhok = rho(k_min) 582 dz = avg_grid%z(k_min + 1) - avg_grid%z(k_min) 583 drhodz = 0.5_wp * (rho(k_min + 1) - rho(k_min)) / dz 584 585 p(1) = constant_density_pressure(p(k_min), zk, rhok, drhodz, & 586 0.0_wp, G) 587 588 END SUBROUTINE get_surface_pressure 589 590 591 FUNCTION constant_density_pressure(pk, zk, rhok, drhodz, z, g) RESULT(p) 592 593 REAL(wp), INTENT(IN) :: pk, zk, rhok, drhodz, g, z 594 REAL(wp) :: p 595 596 p = pk + ( zk - z ) * g * ( rhok + 0.5*drhodz * (zk - z) ) 597 598 END FUNCTION constant_density_pressure 594 599 595 600 !-----------------------------------------------------------------------------! … … 599 604 !> vg. 600 605 !-----------------------------------------------------------------------------! 601 602 603 604 REAL( dp), DIMENSION(:), INTENT(IN) :: p_north, p_south, p_east, p_west, &606 SUBROUTINE geostrophic_winds(p_north, p_south, p_east, p_west, rho, f3, & 607 Lx, Ly, phi_n, lam_n, phi_g, lam_g, ug, vg) 608 609 REAL(wp), DIMENSION(:), INTENT(IN) :: p_north, p_south, p_east, p_west, & 605 610 rho 606 REAL( dp), INTENT(IN) :: f3, Lx, Ly, phi_n, lam_n, phi_g, lam_g607 REAL( dp), DIMENSION(:), INTENT(OUT) :: ug, vg608 REAL( dp) :: facx, facy609 610 facx = 1.0_ dp / (Lx * f3)611 facy = 1.0_ dp / (Ly * f3)611 REAL(wp), INTENT(IN) :: f3, Lx, Ly, phi_n, lam_n, phi_g, lam_g 612 REAL(wp), DIMENSION(:), INTENT(OUT) :: ug, vg 613 REAL(wp) :: facx, facy 614 615 facx = 1.0_wp / (Lx * f3) 616 facy = 1.0_wp / (Ly * f3) 612 617 ug(:) = - facy / rho(:) * (p_north(:) - p_south(:)) 613 618 vg(:) = facx / rho(:) * (p_east(:) - p_west(:)) … … 617 622 ) 618 623 619 624 END SUBROUTINE geostrophic_winds 620 625 621 626 … … 627 632 !> lngitude of a geographical system centered at x0 and y0. 628 633 !-----------------------------------------------------------------------------! 629 630 REAL(dp), INTENT(IN) :: x(:), y(:), x0, y0, r631 REAL(dp), INTENT(OUT) :: lat(:), lon(:)632 633 REAL(dp) :: ri634 635 ! 636 !-- 637 638 ri = 1.0_dp / r639 640 641 642 634 SUBROUTINE inv_plate_carree(x, y, x0, y0, r, lat, lon) 635 REAL(wp), INTENT(IN) :: x(:), y(:), x0, y0, r 636 REAL(wp), INTENT(OUT) :: lat(:), lon(:) 637 638 REAL(wp) :: ri 639 640 ! 641 !-- TODO check dimensions of lat/lon and y/x match 642 643 ri = 1.0_wp / r 644 645 lat(:) = (y(:) - y0) * ri 646 lon(:) = (x(:) - x0) * ri 647 END SUBROUTINE 643 648 644 649 … … 663 668 !> coordinate xy. 664 669 !------------------------------------------------------------------------------! 665 ELEMENTAL REAL(dp) FUNCTION project(xy, xy0, r)666 REAL(dp), INTENT(IN) :: xy, xy0, r667 REAL(dp) :: ri668 669 ! 670 !-- 671 !-- 672 !-- 673 ri = 1.0_dp / r674 675 676 677 670 ELEMENTAL REAL(wp) FUNCTION project(xy, xy0, r) 671 REAL(wp), INTENT(IN) :: xy, xy0, r 672 REAL(wp) :: ri 673 674 ! 675 !-- If this elemental function is called with a large array as xy, it is 676 !-- computationally more efficient to precompute the inverse radius and 677 !-- then muliply. 678 ri = 1.0_wp / r 679 680 project = (xy - xy0) * ri 681 682 END FUNCTION project 678 683 679 684 … … 684 689 !> compute the geographical latitude of its rotated north pole. 685 690 !------------------------------------------------------------------------------! 686 REAL(dp) FUNCTION phic_to_phin(phi_c)687 REAL(dp), INTENT(IN) :: phi_c688 689 phic_to_phin = 0.5_dp * PI - ABS(phi_c)690 691 691 REAL(wp) FUNCTION phic_to_phin(phi_c) 692 REAL(wp), INTENT(IN) :: phi_c 693 694 phic_to_phin = 0.5_wp * PI - ABS(phi_c) 695 696 END FUNCTION phic_to_phin 692 697 693 698 … … 699 704 !> north pole. 700 705 !------------------------------------------------------------------------------! 701 REAL(dp) FUNCTION lamc_to_lamn(phi_c, lam_c)702 REAL(dp), INTENT(IN) :: phi_c, lam_c703 704 705 IF (phi_c > 0.0_dp) THEN706 707 708 709 706 REAL(wp) FUNCTION lamc_to_lamn(phi_c, lam_c) 707 REAL(wp), INTENT(IN) :: phi_c, lam_c 708 709 lamc_to_lamn = lam_c 710 IF (phi_c > 0.0_wp) THEN 711 lamc_to_lamn = lam_c - SIGN(PI, lam_c) 712 ENDIF 713 714 END FUNCTION lamc_to_lamn 710 715 711 716 … … 718 723 !> rotated-pole coordinate transformations. 719 724 !------------------------------------------------------------------------------! 720 REAL(dp) FUNCTION gamma_from_hemisphere(phi_cg, phi_ref)721 REAL(dp), INTENT(IN) :: phi_cg722 REAL(dp), INTENT(IN) :: phi_ref723 724 725 726 727 728 729 730 731 gamma_from_hemisphere = 0.0_dp732 733 725 REAL(wp) FUNCTION gamma_from_hemisphere(phi_cg, phi_ref) 726 REAL(wp), INTENT(IN) :: phi_cg 727 REAL(wp), INTENT(IN) :: phi_ref 728 729 LOGICAL :: palm_origin_is_south_of_cosmo_origin 730 731 palm_origin_is_south_of_cosmo_origin = (phi_cg < phi_ref) 732 733 IF (palm_origin_is_south_of_cosmo_origin) THEN 734 gamma_from_hemisphere = PI 735 ELSE 736 gamma_from_hemisphere = 0.0_wp 737 ENDIF 738 END FUNCTION gamma_from_hemisphere 734 739 735 740 … … 760 765 !> phi(:,:), lam(:,:): geographical latitudes and logitudes 761 766 !------------------------------------------------------------------------------! 762 SUBROUTINE rotate_to_cosmo(phir, lamr, phip, lamp, phi, lam, gam) 763 REAL(dp), INTENT(IN) :: phir(0:), lamr(0:), phip, lamp, gam 764 REAL(dp), INTENT(OUT) :: phi(0:,0:), lam(0:,0:) 765 766 INTEGER :: i, j 767 SUBROUTINE rotate_to_cosmo(phir, lamr, phip, lamp, phi, lam, gam) 768 REAL(wp), INTENT(IN) :: phir(0:), lamr(0:), phip, lamp, gam 769 REAL(wp), INTENT(OUT) :: phi(0:,0:), lam(0:,0:) 770 771 INTEGER :: i, j 772 773 IF ( SIZE(phi, 1) .NE. SIZE(lam, 1) .OR. & 774 SIZE(phi, 2) .NE. SIZE(lam, 2) ) THEN 775 PRINT *, "inifor: rotate_to_cosmo: Dimensions of phi and lambda do not match. Dimensions are:" 776 PRINT *, "inifor: rotate_to_cosmo: phi: ", SIZE(phi, 1), SIZE(phi, 2) 777 PRINT *, "inifor: rotate_to_cosmo: lam: ", SIZE(lam, 1), SIZE(lam, 2) 778 STOP 779 ENDIF 780 781 IF ( SIZE(phir) .NE. SIZE(phi, 2) .OR. & 782 SIZE(lamr) .NE. SIZE(phi, 1) ) THEN 783 PRINT *, "inifor: rotate_to_cosmo: Dimensions of phir and lamr do not match. Dimensions are:" 784 PRINT *, "inifor: rotate_to_cosmo: phir: ", SIZE(phir), SIZE(phi, 2) 785 PRINT *, "inifor: rotate_to_cosmo: lamr: ", SIZE(lamr), SIZE(phi, 1) 786 STOP 787 ENDIF 788 789 DO j = 0, UBOUND(phir, 1) 790 DO i = 0, UBOUND(lamr, 1) 791 792 phi(i,j) = phirot2phi(phir(j) * TO_DEGREES, & 793 lamr(i) * TO_DEGREES, & 794 phip * TO_DEGREES, & 795 gam * TO_DEGREES) * TO_RADIANS 796 797 lam(i,j) = rlarot2rla(phir(j) * TO_DEGREES, & 798 lamr(i) * TO_DEGREES, & 799 phip * TO_DEGREES, & 800 lamp * TO_DEGREES, & 801 gam * TO_DEGREES) * TO_RADIANS 802 803 ENDDO 804 ENDDO 805 806 END SUBROUTINE rotate_to_cosmo 767 807 768 IF ( SIZE(phi, 1) .NE. SIZE(lam, 1) .OR. &769 SIZE(phi, 2) .NE. SIZE(lam, 2) ) THEN770 PRINT *, "inifor: rotate_to_cosmo: Dimensions of phi and lambda do not match. Dimensions are:"771 PRINT *, "inifor: rotate_to_cosmo: phi: ", SIZE(phi, 1), SIZE(phi, 2)772 PRINT *, "inifor: rotate_to_cosmo: lam: ", SIZE(lam, 1), SIZE(lam, 2)773 STOP774 ENDIF775 776 IF ( SIZE(phir) .NE. SIZE(phi, 2) .OR. &777 SIZE(lamr) .NE. SIZE(phi, 1) ) THEN778 PRINT *, "inifor: rotate_to_cosmo: Dimensions of phir and lamr do not match. Dimensions are:"779 PRINT *, "inifor: rotate_to_cosmo: phir: ", SIZE(phir), SIZE(phi, 2)780 PRINT *, "inifor: rotate_to_cosmo: lamr: ", SIZE(lamr), SIZE(phi, 1)781 STOP782 ENDIF783 784 DO j = 0, UBOUND(phir, 1)785 DO i = 0, UBOUND(lamr, 1)786 787 phi(i,j) = phirot2phi(phir(j) * TO_DEGREES, &788 lamr(i) * TO_DEGREES, &789 phip * TO_DEGREES, &790 gam * TO_DEGREES) * TO_RADIANS791 792 lam(i,j) = rlarot2rla(phir(j) * TO_DEGREES, &793 lamr(i) * TO_DEGREES, &794 phip * TO_DEGREES, &795 lamp * TO_DEGREES, &796 gam * TO_DEGREES) * TO_RADIANS797 798 ENDDO799 ENDDO800 801 END SUBROUTINE rotate_to_cosmo802 803 808 804 809 !------------------------------------------------------------------------------! … … 807 812 !> Rotate the given vector field (x(:), y(:)) by the given 'angle'. 808 813 !------------------------------------------------------------------------------! 809 810 REAL(dp), DIMENSION(:), INTENT(INOUT) :: x, y !< x and y coodrinate in arbitrary units811 REAL(dp), INTENT(IN) :: angle !< rotation angle [deg]812 813 814 REAL(dp) :: sine, cosine, v_rot(2), rotation(2,2)815 816 817 818 ! 819 !-- 820 !-- 821 !-- 822 !-- 823 824 825 DOi = LBOUND(x, 1), UBOUND(x, 1)826 827 828 829 830 831 832 833 834 814 SUBROUTINE rotate_vector_field(x, y, angle) 815 REAL(wp), DIMENSION(:), INTENT(INOUT) :: x, y !< x and y coodrinate in arbitrary units 816 REAL(wp), INTENT(IN) :: angle !< rotation angle [deg] 817 818 INTEGER :: i 819 REAL(wp) :: sine, cosine, v_rot(2), rotation(2,2) 820 821 sine = SIN(angle * TO_RADIANS) 822 cosine = COS(angle * TO_RADIANS) 823 ! 824 !-- RESAHPE() fills columns first, so the rotation matrix becomes 825 !-- 826 !-- rotation = [ cosine -sine ] 827 !-- [ sine cosine ] 828 rotation = RESHAPE( (/cosine, sine, -sine, cosine/), (/2, 2/) ) 829 830 DO i = LBOUND(x, 1), UBOUND(x, 1) 831 832 v_rot(:) = MATMUL(rotation, (/x(i), y(i)/)) 833 834 x(i) = v_rot(1) 835 y(i) = v_rot(2) 836 837 ENDDO 838 839 END SUBROUTINE rotate_vector_field 835 840 836 841 … … 847 852 !> https://www.dwd.de/SharedDocs/downloads/DE/modelldokumentationen/nwv/cosmo_d2/cosmo_d2_dbbeschr_aktuell.pdf?__blob=publicationFile&v=2 848 853 !------------------------------------------------------------------------------! 849 850 851 852 REAL(dp), INTENT(IN) :: phi_n, lam_n, phi_g, lam_g853 REAL(dp) :: delta854 855 856 857 858 859 854 FUNCTION meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g) & 855 RESULT(delta) 856 857 REAL(wp), INTENT(IN) :: phi_n, lam_n, phi_g, lam_g 858 REAL(wp) :: delta 859 860 delta = atan2( COS(phi_n) * SIN(lam_n - lam_g), & 861 COS(phi_g) * SIN(phi_n) - & 862 SIN(phi_g) * COS(phi_n) * COS(lam_n - lam_g) ) 863 864 END FUNCTION meridian_convergence_rotated 860 865 861 866 !------------------------------------------------------------------------------! … … 903 908 !> 904 909 !------------------------------------------------------------------------------! 905 906 907 908 909 REAL(dp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon910 REAL(dp), DIMENSION(0:,0:), INTENT(IN) :: palm_clat, palm_clon911 REAL(dp) :: cosmo_dxi, cosmo_dyi912 913 914 REAL(dp) :: lonpos, latpos, lon0, lat0915 916 917 918 919 cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))920 cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))921 922 DO j = 0, UBOUND(palm_clon, 2)!palm_grid %ny923 DO i = 0, UBOUND(palm_clon, 1)!palm_grid %nx924 ! 925 !-- 926 !-- 927 928 929 930 IF (lonpos < 0.0_dp .OR. latpos < 0.0_dp) THEN931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 910 SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon, & 911 palm_clat, palm_clon, & 912 palm_ii, palm_jj) 913 914 REAL(wp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon 915 REAL(wp), DIMENSION(0:,0:), INTENT(IN) :: palm_clat, palm_clon 916 REAL(wp) :: cosmo_dxi, cosmo_dyi 917 INTEGER, DIMENSION(0:,0:,1:), INTENT(OUT) :: palm_ii, palm_jj 918 919 REAL(wp) :: lonpos, latpos, lon0, lat0 920 INTEGER :: i, j 921 922 lon0 = cosmo_lon(0) 923 lat0 = cosmo_lat(0) 924 cosmo_dxi = 1.0_wp / (cosmo_lon(1) - cosmo_lon(0)) 925 cosmo_dyi = 1.0_wp / (cosmo_lat(1) - cosmo_lat(0)) 926 927 DO j = 0, UBOUND(palm_clon, 2)!palm_grid%ny 928 DO i = 0, UBOUND(palm_clon, 1)!palm_grid%nx 929 ! 930 !-- Compute the floating point index corrseponding to PALM-4U grid point 931 !-- location along target grid (COSMO-DE) axes. 932 lonpos = (palm_clon(i,j) - lon0) * cosmo_dxi 933 latpos = (palm_clat(i,j) - lat0) * cosmo_dyi 934 935 IF (lonpos < 0.0_wp .OR. latpos < 0.0_wp) THEN 936 message = "lonpos or latpos out of bounds " // & 937 "while finding interpolation neighbours!" // NEW_LINE(' ') // & 938 " (i,j) = (" // & 939 TRIM(str(i)) // ", " // TRIM(str(j)) // ")" // NEW_LINE(' ') //& 940 " lonpos " // TRIM(real_to_str(lonpos*TO_DEGREES)) // & 941 ", latpos " // TRIM(real_to_str(latpos*TO_DEGREES)) // NEW_LINE(' ') // & 942 " lon0 " // TRIM(real_to_str(lon0 *TO_DEGREES)) // & 943 ", lat0 " // TRIM(real_to_str(lat0*TO_DEGREES)) // NEW_LINE(' ') // & 944 " PALM lon " // TRIM(real_to_str(palm_clon(i,j)*TO_DEGREES)) // & 945 ", PALM lat " // TRIM(real_to_str(palm_clat(i,j)*TO_DEGREES)) 946 CALL inifor_abort('find_horizontal_neighbours', message) 947 ENDIF 948 949 palm_ii(i,j,1) = FLOOR(lonpos) 950 palm_ii(i,j,2) = FLOOR(lonpos) 951 palm_ii(i,j,3) = CEILING(lonpos) 952 palm_ii(i,j,4) = CEILING(lonpos) 953 954 palm_jj(i,j,1) = FLOOR(latpos) 955 palm_jj(i,j,2) = CEILING(latpos) 956 palm_jj(i,j,3) = CEILING(latpos) 957 palm_jj(i,j,4) = FLOOR(latpos) 958 ENDDO 959 ENDDO 960 961 END SUBROUTINE find_horizontal_neighbours 957 962 958 963 … … 963 968 !> column of the given palm grid. 964 969 !------------------------------------------------------------------------------! 965 970 SUBROUTINE find_vertical_neighbours_and_weights_interp( palm_grid, & 966 971 palm_intermediate ) 967 968 969 970 971 972 973 REAL(dp) :: current_height, column_base, column_top, h_top, h_bottom, &974 975 976 nx = palm_grid %nx977 ny = palm_grid %ny978 nz = palm_grid %nz979 nlev = palm_intermediate %nz980 981 ! 982 !-- 983 DOj = 0, ny984 DOi = 0, nx985 986 987 988 column_base = palm_intermediate %h(i,j,0)989 column_top = palm_intermediate %h(i,j,nlev)990 991 ! 992 !-- 993 !-- 994 !-- 995 !-- 996 DOk = 1, nz997 998 ! 999 !-- 1000 !-- 1001 current_height = palm_grid % z(k) + palm_grid %z01002 h_top = palm_intermediate %h(i,j,k_intermediate+1)1003 h_bottom = palm_intermediate %h(i,j,k_intermediate)1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 ! 1014 !-- 1015 palm_grid % w_verti(i,j,k,1:2) = 0.0_dp1016 1017 1018 1019 palm_grid %kk(i,j,k,1:2) = nlev1020 palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp1021 1022 1023 1024 1025 1026 1027 palm_grid %kk(i,j,k,1:2) = 01028 palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp1029 1030 1031 ! 1032 !-- 1033 !-- 1034 DOWHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)1035 1036 1037 h_top = palm_intermediate %h(i,j,k_intermediate+1)1038 h_bottom = palm_intermediate %h(i,j,k_intermediate)1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 972 TYPE(grid_definition), INTENT(INOUT) :: palm_grid 973 TYPE(grid_definition), INTENT(IN) :: palm_intermediate 974 975 INTEGER :: i, j, k, nx, ny, nz, nlev, k_intermediate 976 LOGICAL :: point_is_below_grid, point_is_above_grid, & 977 point_is_in_current_cell 978 REAL(wp) :: current_height, column_base, column_top, h_top, h_bottom, & 979 weight 980 981 nx = palm_grid%nx 982 ny = palm_grid%ny 983 nz = palm_grid%nz 984 nlev = palm_intermediate%nz 985 986 ! 987 !-- in each column of the fine grid, find vertical neighbours of every cell 988 DO j = 0, ny 989 DO i = 0, nx 990 991 k_intermediate = 0 992 993 column_base = palm_intermediate%h(i,j,0) 994 column_top = palm_intermediate%h(i,j,nlev) 995 996 ! 997 !-- scan through palm_grid column and set neighbour indices in 998 !-- case current_height is either below column_base, in the current 999 !-- cell, or above column_top. Keep increasing current cell index until 1000 !-- the current cell overlaps with the current_height. 1001 DO k = 1, nz 1002 1003 ! 1004 !-- Memorize the top and bottom boundaries of the coarse cell and the 1005 !-- current height within it 1006 current_height = palm_grid%z(k) + palm_grid%z0 1007 h_top = palm_intermediate%h(i,j,k_intermediate+1) 1008 h_bottom = palm_intermediate%h(i,j,k_intermediate) 1009 1010 point_is_above_grid = (current_height > column_top) !22000m, very unlikely 1011 point_is_below_grid = (current_height < column_base) 1012 1013 point_is_in_current_cell = ( & 1014 current_height >= h_bottom .AND. & 1015 current_height < h_top & 1016 ) 1017 1018 ! 1019 !-- set default weights 1020 palm_grid%w_verti(i,j,k,1:2) = 0.0_wp 1021 1022 IF (point_is_above_grid) THEN 1023 1024 palm_grid%kk(i,j,k,1:2) = nlev 1025 palm_grid%w_verti(i,j,k,1:2) = - 2.0_wp 1026 1027 message = "PALM-4U grid extends above COSMO-DE model top." 1028 CALL inifor_abort('find_vertical_neighbours_and_weights', message) 1029 1030 ELSE IF (point_is_below_grid) THEN 1031 1032 palm_grid%kk(i,j,k,1:2) = 0 1033 palm_grid%w_verti(i,j,k,1:2) = - 2.0_wp 1034 1035 ELSE 1036 ! 1037 !-- cycle through intermediate levels until current 1038 !-- intermediate-grid cell overlaps with current_height 1039 DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1) 1040 k_intermediate = k_intermediate + 1 1041 1042 h_top = palm_intermediate%h(i,j,k_intermediate+1) 1043 h_bottom = palm_intermediate%h(i,j,k_intermediate) 1044 point_is_in_current_cell = ( & 1045 current_height >= h_bottom .AND. & 1046 current_height < h_top & 1047 ) 1048 ENDDO 1049 1050 IF (k_intermediate > nlev-1) THEN 1051 message = "Index " // TRIM(str(k_intermediate)) // & 1052 " is above intermediate grid range." 1053 CALL inifor_abort('find_vertical_neighbours', message) 1054 ENDIF 1050 1055 1051 palm_grid % kk(i,j,k,1) = k_intermediate 1052 palm_grid % kk(i,j,k,2) = k_intermediate + 1 1053 1054 ! 1055 !-- compute vertical weights 1056 weight = (h_top - current_height) / (h_top - h_bottom) 1057 palm_grid % w_verti(i,j,k,1) = weight 1058 palm_grid % w_verti(i,j,k,2) = 1.0_dp - weight 1059 ENDIF 1060 1061 ENDDO 1056 palm_grid%kk(i,j,k,1) = k_intermediate 1057 palm_grid%kk(i,j,k,2) = k_intermediate + 1 1058 1059 ! 1060 !-- compute vertical weights 1061 weight = (h_top - current_height) / (h_top - h_bottom) 1062 palm_grid%w_verti(i,j,k,1) = weight 1063 palm_grid%w_verti(i,j,k,2) = 1.0_wp - weight 1064 ENDIF 1062 1065 1063 1066 ENDDO 1064 ENDDO 1065 1066 END SUBROUTINE find_vertical_neighbours_and_weights_interp 1067 1068 ENDDO 1069 ENDDO 1070 1071 END SUBROUTINE find_vertical_neighbours_and_weights_interp 1067 1072 1068 1073 … … 1079 1084 !> iii(:) and jjj(:). 1080 1085 !------------------------------------------------------------------------------! 1081 SUBROUTINE find_vertical_neighbours_and_weights_average( & 1082 avg_grid, level_based_averaging & 1083 ) 1084 1085 TYPE(grid_definition), INTENT(INOUT), TARGET :: avg_grid 1086 LOGICAL :: level_based_averaging 1087 1088 INTEGER :: i, j, k_palm, k_intermediate, l, nlev 1089 LOGICAL :: point_is_below_grid, point_is_above_grid, & 1090 point_is_in_current_cell 1091 REAL(dp) :: current_height, column_base, column_top, h_top, & 1092 h_bottom, weight 1093 REAL(dp), POINTER :: cosmo_h(:,:,:) 1094 1095 1096 avg_grid % k_min = LBOUND(avg_grid % z, 1) 1097 1098 nlev = SIZE(avg_grid % cosmo_h, 3) 1086 SUBROUTINE find_vertical_neighbours_and_weights_average( & 1087 avg_grid, level_based_averaging & 1088 ) 1089 1090 TYPE(grid_definition), INTENT(INOUT), TARGET :: avg_grid 1091 LOGICAL :: level_based_averaging 1092 1093 INTEGER :: i, j, k_palm, k_intermediate, l, nlev 1094 LOGICAL :: point_is_below_grid, point_is_above_grid, & 1095 point_is_in_current_cell 1096 REAL(wp) :: current_height, column_base, column_top, h_top, & 1097 h_bottom, weight 1098 REAL(wp), POINTER :: cosmo_h(:,:,:) 1099 1100 1101 avg_grid%k_min = LBOUND(avg_grid%z, 1) 1102 1103 nlev = SIZE(avg_grid%cosmo_h, 3) 1104 1105 IF (level_based_averaging) THEN 1106 cosmo_h => avg_grid%h 1107 ELSE 1108 cosmo_h => avg_grid%cosmo_h 1109 ENDIF 1110 1111 ! 1112 !-- in each column of the fine grid, find vertical neighbours of every cell 1113 DO l = 1, avg_grid%n_columns 1099 1114 1100 1115 IF (level_based_averaging) THEN 1101 cosmo_h => avg_grid % h 1116 i = 1 1117 j = 1 1102 1118 ELSE 1103 cosmo_h => avg_grid % cosmo_h 1119 i = avg_grid%iii(l) 1120 j = avg_grid%jjj(l) 1104 1121 ENDIF 1105 1122 1106 ! 1107 !-- in each column of the fine grid, find vertical neighbours of every cell 1108 DO l = 1, avg_grid % n_columns 1109 1110 IF (level_based_averaging) THEN 1111 i = 1 1112 j = 1 1123 column_base = cosmo_h(i,j,1) 1124 column_top = cosmo_h(i,j,nlev) 1125 1126 ! 1127 !-- scan through avg_grid column until and set neighbour indices in 1128 !-- case current_height is either below column_base, in the current 1129 !-- cell, or above column_top. Keep increasing current cell index until 1130 !-- the current cell overlaps with the current_height. 1131 k_intermediate = 1 !avg_grid%cosmo_h is indezed 1-based. 1132 DO k_palm = 1, avg_grid%nz 1133 1134 ! 1135 !-- Memorize the top and bottom boundaries of the coarse cell and the 1136 !-- current height within it 1137 current_height = avg_grid%z(k_palm) + avg_grid%z0 1138 h_top = cosmo_h(i,j,k_intermediate+1) 1139 h_bottom = cosmo_h(i,j,k_intermediate) 1140 1141 ! 1142 !-- COSMO column top is located at 22000m, point_is_above_grid is very 1143 !-- unlikely. 1144 point_is_above_grid = (current_height > column_top) 1145 point_is_below_grid = (current_height < column_base) 1146 1147 point_is_in_current_cell = ( & 1148 current_height >= h_bottom .AND. & 1149 current_height < h_top & 1150 ) 1151 1152 ! 1153 !-- set default weights 1154 avg_grid%w(l,k_palm,1:2) = 0.0_wp 1155 1156 IF (point_is_above_grid) THEN 1157 1158 avg_grid%kkk(l,k_palm,1:2) = nlev 1159 avg_grid%w(l,k_palm,1:2) = - 2.0_wp 1160 1161 message = "PALM-4U grid extends above COSMO-DE model top." 1162 CALL inifor_abort('find_vertical_neighbours_and_weights_average', message) 1163 1164 ELSE IF (point_is_below_grid) THEN 1165 1166 avg_grid%kkk(l,k_palm,1:2) = 0 1167 avg_grid%w(l,k_palm,1:2) = - 2.0_wp 1168 avg_grid%k_min = MAX(k_palm + 1, avg_grid%k_min) 1113 1169 ELSE 1114 i = avg_grid % iii(l) 1115 j = avg_grid % jjj(l) 1170 ! 1171 !-- cycle through intermediate levels until current 1172 !-- intermediate-grid cell overlaps with current_height 1173 DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1) 1174 k_intermediate = k_intermediate + 1 1175 1176 h_top = cosmo_h(i,j,k_intermediate+1) 1177 h_bottom = cosmo_h(i,j,k_intermediate) 1178 point_is_in_current_cell = ( & 1179 current_height >= h_bottom .AND. & 1180 current_height < h_top & 1181 ) 1182 ENDDO 1183 1184 ! 1185 !-- k_intermediate = 48 indicates the last section (indices 48 and 49), i.e. 1186 !-- k_intermediate = 49 is not the beginning of a valid cell. 1187 IF (k_intermediate > nlev-1) THEN 1188 message = "Index " // TRIM(str(k_intermediate)) // & 1189 " is above intermediate grid range." 1190 CALL inifor_abort('find_vertical_neighbours', message) 1191 ENDIF 1192 1193 avg_grid%kkk(l,k_palm,1) = k_intermediate 1194 avg_grid%kkk(l,k_palm,2) = k_intermediate + 1 1195 1196 ! 1197 !-- compute vertical weights 1198 weight = (h_top - current_height) / (h_top - h_bottom) 1199 avg_grid%w(l,k_palm,1) = weight 1200 avg_grid%w(l,k_palm,2) = 1.0_wp - weight 1116 1201 ENDIF 1117 1202 1118 column_base = cosmo_h(i,j,1) 1119 column_top = cosmo_h(i,j,nlev) 1120 1121 ! 1122 !-- scan through avg_grid column until and set neighbour indices in 1123 !-- case current_height is either below column_base, in the current 1124 !-- cell, or above column_top. Keep increasing current cell index until 1125 !-- the current cell overlaps with the current_height. 1126 k_intermediate = 1 !avg_grid % cosmo_h is indezed 1-based. 1127 DO k_palm = 1, avg_grid % nz 1128 1129 ! 1130 !-- Memorize the top and bottom boundaries of the coarse cell and the 1131 !-- current height within it 1132 current_height = avg_grid % z(k_palm) + avg_grid % z0 1133 h_top = cosmo_h(i,j,k_intermediate+1) 1134 h_bottom = cosmo_h(i,j,k_intermediate) 1135 1136 ! 1137 !-- COSMO column top is located at 22000m, point_is_above_grid is very 1138 !-- unlikely. 1139 point_is_above_grid = (current_height > column_top) 1140 point_is_below_grid = (current_height < column_base) 1141 1142 point_is_in_current_cell = ( & 1143 current_height >= h_bottom .AND. & 1144 current_height < h_top & 1145 ) 1146 1147 ! 1148 !-- set default weights 1149 avg_grid % w(l,k_palm,1:2) = 0.0_dp 1150 1151 IF (point_is_above_grid) THEN 1152 1153 avg_grid % kkk(l,k_palm,1:2) = nlev 1154 avg_grid % w(l,k_palm,1:2) = - 2.0_dp 1155 1156 message = "PALM-4U grid extends above COSMO-DE model top." 1157 CALL inifor_abort('find_vertical_neighbours_and_weights_average', message) 1158 1159 ELSE IF (point_is_below_grid) THEN 1160 1161 avg_grid % kkk(l,k_palm,1:2) = 0 1162 avg_grid % w(l,k_palm,1:2) = - 2.0_dp 1163 avg_grid % k_min = MAX(k_palm + 1, avg_grid % k_min) 1164 ELSE 1165 ! 1166 !-- cycle through intermediate levels until current 1167 !-- intermediate-grid cell overlaps with current_height 1168 DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1) 1169 k_intermediate = k_intermediate + 1 1170 1171 h_top = cosmo_h(i,j,k_intermediate+1) 1172 h_bottom = cosmo_h(i,j,k_intermediate) 1173 point_is_in_current_cell = ( & 1174 current_height >= h_bottom .AND. & 1175 current_height < h_top & 1176 ) 1177 ENDDO 1178 1179 ! 1180 !-- k_intermediate = 48 indicates the last section (indices 48 and 49), i.e. 1181 !-- k_intermediate = 49 is not the beginning of a valid cell. 1182 IF (k_intermediate > nlev-1) THEN 1183 message = "Index " // TRIM(str(k_intermediate)) // & 1184 " is above intermediate grid range." 1185 CALL inifor_abort('find_vertical_neighbours', message) 1186 ENDIF 1187 1188 avg_grid % kkk(l,k_palm,1) = k_intermediate 1189 avg_grid % kkk(l,k_palm,2) = k_intermediate + 1 1190 1191 ! 1192 !-- compute vertical weights 1193 weight = (h_top - current_height) / (h_top - h_bottom) 1194 avg_grid % w(l,k_palm,1) = weight 1195 avg_grid % w(l,k_palm,2) = 1.0_dp - weight 1196 ENDIF 1197 1198 ! 1199 !-- Loop over PALM levels k 1200 ENDDO 1201 1202 ! 1203 !-- Loop over averaging columns l 1203 ! 1204 !-- Loop over PALM levels k 1204 1205 ENDDO 1206 1207 ! 1208 !-- Loop over averaging columns l 1209 ENDDO 1205 1210 1206 1211 END SUBROUTINE find_vertical_neighbours_and_weights_average 1207 1212 1208 1213 !------------------------------------------------------------------------------! … … 1215 1220 !> Input parameters: 1216 1221 !> ----------------- 1217 !> palm_grid %clon : longitudes of PALM-4U scalars (cell centres) in COSMO-DE's rotated-pole grid [rad]1218 !> 1219 !> palm_grid %clat : latitudes of PALM-4U cell centres in COSMO-DE's rotated-pole grid [rad]1220 !> 1221 !> cosmo_grid %lon : rotated-pole longitudes of scalars (cell centres) of the COSMO-DE grid [rad]1222 !> 1223 !> cosmo_grid %lat : rotated-pole latitudes of scalars (cell centers) of the COSMO-DE grid [rad]1224 !> 1225 !> cosmo_grid %dxi : inverse grid spacing in the first dimension [m^-1]1226 !> 1227 !> cosmo_grid %dyi : inverse grid spacing in the second dimension [m^-1]1222 !> palm_grid%clon : longitudes of PALM-4U scalars (cell centres) in COSMO-DE's rotated-pole grid [rad] 1223 !> 1224 !> palm_grid%clat : latitudes of PALM-4U cell centres in COSMO-DE's rotated-pole grid [rad] 1225 !> 1226 !> cosmo_grid%lon : rotated-pole longitudes of scalars (cell centres) of the COSMO-DE grid [rad] 1227 !> 1228 !> cosmo_grid%lat : rotated-pole latitudes of scalars (cell centers) of the COSMO-DE grid [rad] 1229 !> 1230 !> cosmo_grid%dxi : inverse grid spacing in the first dimension [m^-1] 1231 !> 1232 !> cosmo_grid%dyi : inverse grid spacing in the second dimension [m^-1] 1228 1233 !> 1229 1234 !> Output parameters: 1230 1235 !> ------------------ 1231 !> palm_grid %w_horiz(:,:,1-4) : weights for bilinear horizontal interpolation1236 !> palm_grid%w_horiz(:,:,1-4) : weights for bilinear horizontal interpolation 1232 1237 ! 1233 1238 ! COSMO-DE grid … … 1251 1256 ! 1252 1257 !------------------------------------------------------------------------------! 1253 1254 1258 SUBROUTINE compute_horizontal_interp_weights(cosmo_lat, cosmo_lon, & 1259 palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz) 1255 1260 1256 REAL(dp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon1257 REAL(dp) :: cosmo_dxi, cosmo_dyi1258 REAL(dp), DIMENSION(0:,0:), INTENT(IN) :: palm_clat, palm_clon1259 1260 1261 REAL(dp), DIMENSION(0:,0:,1:), INTENT(OUT) :: palm_w_horiz1262 1263 REAL(dp) :: wl, wp1264 1265 1266 cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))1267 cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))1268 1269 DOj = 0, UBOUND(palm_clon, 2)1270 DOi = 0, UBOUND(palm_clon, 1)1271 1272 ! 1273 !-- 1274 wl= ( cosmo_lon(palm_ii(i,j,4)) - palm_clon(i,j) ) * cosmo_dxi1275 1276 ! 1277 !-- 1278 wp= ( cosmo_lat(palm_jj(i,j,2)) - palm_clat(i,j) ) * cosmo_dyi1279 1280 IF (wl > 1.0_dp .OR. wl < 0.0_dp) THEN1281 message = "Horizontal weight wl = " // TRIM(real_to_str(wl)) // &1282 1283 1284 1285 IF (wp > 1.0_dp .OR. wp < 0.0_dp) THEN1286 message = "Horizontal weight wp = " // TRIM(real_to_str(wp)) // &1287 1288 1289 1290 1291 palm_w_horiz(i,j,1) = wl * wp1292 palm_w_horiz(i,j,2) = wl * (1.0_dp - wp)1293 palm_w_horiz(i,j,3) = (1.0_dp - wl) * (1.0_dp - wp)1294 palm_w_horiz(i,j,4) = 1.0_dp - SUM( palm_w_horiz(i,j,1:3) )1295 1296 1297 1261 REAL(wp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon 1262 REAL(wp) :: cosmo_dxi, cosmo_dyi 1263 REAL(wp), DIMENSION(0:,0:), INTENT(IN) :: palm_clat, palm_clon 1264 INTEGER, DIMENSION(0:,0:,1:), INTENT(IN) :: palm_ii, palm_jj 1265 1266 REAL(wp), DIMENSION(0:,0:,1:), INTENT(OUT) :: palm_w_horiz 1267 1268 REAL(wp) :: wlambda, wphi 1269 INTEGER :: i, j 1270 1271 cosmo_dxi = 1.0_wp / (cosmo_lon(1) - cosmo_lon(0)) 1272 cosmo_dyi = 1.0_wp / (cosmo_lat(1) - cosmo_lat(0)) 1273 1274 DO j = 0, UBOUND(palm_clon, 2) 1275 DO i = 0, UBOUND(palm_clon, 1) 1276 1277 ! 1278 !-- weight in lambda direction 1279 wlambda = ( cosmo_lon(palm_ii(i,j,4)) - palm_clon(i,j) ) * cosmo_dxi 1280 1281 ! 1282 !-- weight in phi direction 1283 wphi = ( cosmo_lat(palm_jj(i,j,2)) - palm_clat(i,j) ) * cosmo_dyi 1284 1285 IF (wlambda > 1.0_wp .OR. wlambda < 0.0_wp) THEN 1286 message = "Horizontal weight wlambda = " // TRIM(real_to_str(wlambda)) // & 1287 " is out bounds." 1288 CALL inifor_abort('compute_horizontal_interp_weights', message) 1289 ENDIF 1290 IF (wphi > 1.0_wp .OR. wphi < 0.0_wp) THEN 1291 message = "Horizontal weight wphi = " // TRIM(real_to_str(wphi)) // & 1292 " is out bounds." 1293 CALL inifor_abort('compute_horizontal_interp_weights', message) 1294 ENDIF 1295 1296 palm_w_horiz(i,j,1) = wlambda * wphi 1297 palm_w_horiz(i,j,2) = wlambda * (1.0_wp - wphi) 1298 palm_w_horiz(i,j,3) = (1.0_wp - wlambda) * (1.0_wp - wphi) 1299 palm_w_horiz(i,j,4) = 1.0_wp - SUM( palm_w_horiz(i,j,1:3) ) 1300 1301 ENDDO 1302 ENDDO 1298 1303 1299 1304 END SUBROUTINE compute_horizontal_interp_weights 1300 1305 1301 1306 … … 1311 1316 !> which means the first centre point has to be omitted and is set to zero. 1312 1317 !------------------------------------------------------------------------------! 1313 1314 REAL(dp), DIMENSION(0:,0:,0:), INTENT(IN) :: u_face, v_face1315 REAL(dp), DIMENSION(0:,0:,0:), INTENT(OUT) :: u_centre, v_centre1316 1317 1318 1319 1320 1321 u_centre(0,:,:) = 0.0_dp1322 u_centre(1:,:,:) = 0.5_dp * ( u_face(0:nx-1,:,:) + u_face(1:,:,:) )1323 1324 v_centre(:,0,:) = 0.0_dp1325 v_centre(:,1:,:) = 0.5_dp * ( v_face(:,0:ny-1,:) + v_face(:,1:,:) )1326 1318 SUBROUTINE centre_velocities(u_face, v_face, u_centre, v_centre) 1319 REAL(wp), DIMENSION(0:,0:,0:), INTENT(IN) :: u_face, v_face 1320 REAL(wp), DIMENSION(0:,0:,0:), INTENT(OUT) :: u_centre, v_centre 1321 INTEGER :: nx, ny 1322 1323 nx = UBOUND(u_face, 1) 1324 ny = UBOUND(u_face, 2) 1325 1326 u_centre(0,:,:) = 0.0_wp 1327 u_centre(1:,:,:) = 0.5_wp * ( u_face(0:nx-1,:,:) + u_face(1:,:,:) ) 1328 1329 v_centre(:,0,:) = 0.0_wp 1330 v_centre(:,1:,:) = 0.5_wp * ( v_face(:,0:ny-1,:) + v_face(:,1:,:) ) 1331 END SUBROUTINE centre_velocities 1327 1332 1328 1333 … … 1332 1337 !> Compute the geographical latitude of a point given in rotated-pole cordinates 1333 1338 !------------------------------------------------------------------------------! 1334 FUNCTION phirot2phi (phirot, rlarot, polphi, polgam) 1335 1336 REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole 1337 REAL(dp), INTENT (IN) :: phirot !< latitude in the rotated system 1338 REAL(dp), INTENT (IN) :: rlarot !< longitude in the rotated system 1339 REAL(dp), INTENT (IN) :: polgam !< angle between the north poles of the systems 1340 1341 REAL(dp) :: phirot2phi !< latitude in the geographical system 1342 1343 REAL(dp) :: zsinpol, zcospol, zphis, zrlas, zarg, zgam 1344 1345 zsinpol = SIN(polphi * TO_RADIANS) 1346 zcospol = COS(polphi * TO_RADIANS) 1347 zphis = phirot * TO_RADIANS 1348 1349 IF (rlarot > 180.0_dp) THEN 1350 zrlas = rlarot - 360.0_dp 1351 ELSE 1352 zrlas = rlarot 1353 ENDIF 1354 zrlas = zrlas * TO_RADIANS 1339 FUNCTION phirot2phi (phirot, rlarot, polphi, polgam) 1340 1341 REAL(wp), INTENT (IN) :: polphi !< latitude of the rotated north pole 1342 REAL(wp), INTENT (IN) :: phirot !< latitude in the rotated system 1343 REAL(wp), INTENT (IN) :: rlarot !< longitude in the rotated system 1344 REAL(wp), INTENT (IN) :: polgam !< angle between the north poles of the systems 1345 1346 REAL(wp) :: phirot2phi !< latitude in the geographical system 1347 1348 REAL(wp) :: zsinpol, zcospol, zphis, zrlas, zarg, zgam 1349 1350 zsinpol = SIN(polphi * TO_RADIANS) 1351 zcospol = COS(polphi * TO_RADIANS) 1352 zphis = phirot * TO_RADIANS 1353 1354 IF (rlarot > 180.0_wp) THEN 1355 zrlas = rlarot - 360.0_wp 1356 ELSE 1357 zrlas = rlarot 1358 ENDIF 1359 zrlas = zrlas * TO_RADIANS 1360 1361 IF (polgam /= 0.0_wp) THEN 1362 zgam = polgam * TO_RADIANS 1363 zarg = zsinpol * SIN (zphis) + & 1364 zcospol * COS(zphis) * ( COS(zrlas) * COS(zgam) - & 1365 SIN(zgam) * SIN(zrlas) ) 1366 ELSE 1367 zarg = zcospol * COS (zphis) * COS (zrlas) + zsinpol * SIN (zphis) 1368 ENDIF 1369 1370 phirot2phi = ASIN (zarg) * TO_DEGREES 1371 1372 END FUNCTION phirot2phi 1373 1374 1375 !------------------------------------------------------------------------------! 1376 ! Description: 1377 ! ------------ 1378 !> Compute the geographical latitude of a point given in rotated-pole cordinates 1379 !------------------------------------------------------------------------------! 1380 FUNCTION phi2phirot (phi, rla, polphi, pollam) 1381 1382 REAL(wp), INTENT (IN) :: polphi !< latitude of the rotated north pole 1383 REAL(wp), INTENT (IN) :: pollam !< longitude of the rotated north pole 1384 REAL(wp), INTENT (IN) :: phi !< latitude in the geographical system 1385 REAL(wp), INTENT (IN) :: rla !< longitude in the geographical system 1386 1387 REAL(wp) :: phi2phirot !< longitude in the rotated system 1388 1389 REAL(wp) :: zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1 1390 1391 zsinpol = SIN(polphi * TO_RADIANS) 1392 zcospol = COS(polphi * TO_RADIANS) 1393 zlampol = pollam * TO_RADIANS 1394 zphi = phi * TO_RADIANS 1395 1396 IF (rla > 180.0_wp) THEN 1397 zrla1 = rla - 360.0_wp 1398 ELSE 1399 zrla1 = rla 1400 ENDIF 1401 zrla = zrla1 * TO_RADIANS 1402 1403 zarg1 = SIN(zphi) * zsinpol 1404 zarg2 = COS(zphi) * zcospol * COS(zrla - zlampol) 1405 1406 phi2phirot = ASIN(zarg1 + zarg2) * TO_DEGREES 1407 1408 END FUNCTION phi2phirot 1409 1410 1411 !------------------------------------------------------------------------------! 1412 ! Description: 1413 ! ------------ 1414 !> Compute the geographical longitude of a point given in rotated-pole cordinates 1415 !------------------------------------------------------------------------------! 1416 FUNCTION rlarot2rla(phirot, rlarot, polphi, pollam, polgam) 1417 1418 REAL(wp), INTENT (IN) :: polphi !< latitude of the rotated north pole 1419 REAL(wp), INTENT (IN) :: pollam !< longitude of the rotated north pole 1420 REAL(wp), INTENT (IN) :: phirot !< latitude in the rotated system 1421 REAL(wp), INTENT (IN) :: rlarot !< longitude in the rotated system 1422 REAL(wp), INTENT (IN) :: polgam !< angle between the north poles of the systems 1423 1424 REAL(wp) :: rlarot2rla !< latitude in the geographical system 1425 1426 REAL(wp) :: zsinpol, zcospol, zlampol, zphis, zrlas, zarg1, zarg2, zgam 1427 1428 zsinpol = SIN(TO_RADIANS * polphi) 1429 zcospol = COS(TO_RADIANS * polphi) 1430 zlampol = TO_RADIANS * pollam 1431 zphis = TO_RADIANS * phirot 1432 1433 IF (rlarot > 180.0_wp) THEN 1434 zrlas = rlarot - 360.0_wp 1435 ELSE 1436 zrlas = rlarot 1437 ENDIF 1438 zrlas = TO_RADIANS * zrlas 1439 1440 IF (polgam /= 0.0_wp) THEN 1441 zgam = TO_RADIANS * polgam 1442 zarg1 = SIN(zlampol) * (zcospol * SIN(zphis) - zsinpol*COS(zphis) * & 1443 (COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) - & 1444 COS(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) + & 1445 COS(zrlas) * SIN(zgam) ) 1446 1447 zarg2 = COS (zlampol) * (zcospol * SIN(zphis) - zsinpol*COS(zphis) * & 1448 (COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) + & 1449 SIN(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) + & 1450 COS(zrlas) * SIN(zgam) ) 1451 ELSE 1452 zarg1 = SIN (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis) + & 1453 zcospol * SIN(zphis)) - & 1454 COS (zlampol) * SIN(zrlas) * COS(zphis) 1455 zarg2 = COS (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis) + & 1456 zcospol * SIN(zphis)) + & 1457 SIN (zlampol) * SIN(zrlas) * COS(zphis) 1458 ENDIF 1459 1460 IF (zarg2 == 0.0_wp) zarg2 = 1.0E-20_wp 1461 1462 rlarot2rla = ATAN2(zarg1,zarg2) * TO_DEGREES 1355 1463 1356 IF (polgam /= 0.0_dp) THEN 1357 zgam = polgam * TO_RADIANS 1358 zarg = zsinpol * SIN (zphis) + & 1359 zcospol * COS(zphis) * ( COS(zrlas) * COS(zgam) - & 1360 SIN(zgam) * SIN(zrlas) ) 1361 ELSE 1362 zarg = zcospol * COS (zphis) * COS (zrlas) + zsinpol * SIN (zphis) 1363 ENDIF 1364 1365 phirot2phi = ASIN (zarg) * TO_DEGREES 1366 1367 END FUNCTION phirot2phi 1368 1369 1370 !------------------------------------------------------------------------------! 1371 ! Description: 1372 ! ------------ 1373 !> Compute the geographical latitude of a point given in rotated-pole cordinates 1374 !------------------------------------------------------------------------------! 1375 FUNCTION phi2phirot (phi, rla, polphi, pollam) 1376 1377 REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole 1378 REAL(dp), INTENT (IN) :: pollam !< longitude of the rotated north pole 1379 REAL(dp), INTENT (IN) :: phi !< latitude in the geographical system 1380 REAL(dp), INTENT (IN) :: rla !< longitude in the geographical system 1381 1382 REAL(dp) :: phi2phirot !< longitude in the rotated system 1383 1384 REAL(dp) :: zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1 1385 1386 zsinpol = SIN(polphi * TO_RADIANS) 1387 zcospol = COS(polphi * TO_RADIANS) 1388 zlampol = pollam * TO_RADIANS 1389 zphi = phi * TO_RADIANS 1390 1391 IF (rla > 180.0_dp) THEN 1392 zrla1 = rla - 360.0_dp 1393 ELSE 1394 zrla1 = rla 1395 ENDIF 1396 zrla = zrla1 * TO_RADIANS 1397 1398 zarg1 = SIN(zphi) * zsinpol 1399 zarg2 = COS(zphi) * zcospol * COS(zrla - zlampol) 1400 1401 phi2phirot = ASIN(zarg1 + zarg2) * TO_DEGREES 1402 1403 END FUNCTION phi2phirot 1404 1405 1406 !------------------------------------------------------------------------------! 1407 ! Description: 1408 ! ------------ 1409 !> Compute the geographical longitude of a point given in rotated-pole cordinates 1410 !------------------------------------------------------------------------------! 1411 FUNCTION rlarot2rla(phirot, rlarot, polphi, pollam, polgam) 1412 1413 REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole 1414 REAL(dp), INTENT (IN) :: pollam !< longitude of the rotated north pole 1415 REAL(dp), INTENT (IN) :: phirot !< latitude in the rotated system 1416 REAL(dp), INTENT (IN) :: rlarot !< longitude in the rotated system 1417 REAL(dp), INTENT (IN) :: polgam !< angle between the north poles of the systems 1418 1419 REAL(dp) :: rlarot2rla !< latitude in the geographical system 1420 1421 REAL(dp) :: zsinpol, zcospol, zlampol, zphis, zrlas, zarg1, zarg2, zgam 1422 1423 zsinpol = SIN(TO_RADIANS * polphi) 1424 zcospol = COS(TO_RADIANS * polphi) 1425 zlampol = TO_RADIANS * pollam 1426 zphis = TO_RADIANS * phirot 1427 1428 IF (rlarot > 180.0_dp) THEN 1429 zrlas = rlarot - 360.0_dp 1430 ELSE 1431 zrlas = rlarot 1432 ENDIF 1433 zrlas = TO_RADIANS * zrlas 1434 1435 IF (polgam /= 0.0_dp) THEN 1436 zgam = TO_RADIANS * polgam 1437 zarg1 = SIN(zlampol) * (zcospol * SIN(zphis) - zsinpol*COS(zphis) * & 1438 (COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) - & 1439 COS(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) + & 1440 COS(zrlas) * SIN(zgam) ) 1441 1442 zarg2 = COS (zlampol) * (zcospol * SIN(zphis) - zsinpol*COS(zphis) * & 1443 (COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) + & 1444 SIN(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) + & 1445 COS(zrlas) * SIN(zgam) ) 1446 ELSE 1447 zarg1 = SIN (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis) + & 1448 zcospol * SIN(zphis)) - & 1449 COS (zlampol) * SIN(zrlas) * COS(zphis) 1450 zarg2 = COS (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis) + & 1451 zcospol * SIN(zphis)) + & 1452 SIN (zlampol) * SIN(zrlas) * COS(zphis) 1453 ENDIF 1454 1455 IF (zarg2 == 0.0_dp) zarg2 = 1.0E-20_dp 1456 1457 rlarot2rla = ATAN2(zarg1,zarg2) * TO_DEGREES 1458 1459 END FUNCTION rlarot2rla 1464 END FUNCTION rlarot2rla 1460 1465 1461 1466 … … 1465 1470 !> Compute the rotated-pole longitude of a point given in geographical cordinates 1466 1471 !------------------------------------------------------------------------------! 1467 1468 1469 REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole1470 REAL(dp), INTENT (IN) :: pollam !< longitude of the rotated north pole1471 REAL(dp), INTENT (IN) :: phi !< latitude in geographical system1472 REAL(dp), INTENT (IN) :: rla !< longitude in geographical system1473 REAL(dp), INTENT (IN) :: polgam !< angle between the north poles of the systems1474 1475 REAL (KIND=dp) :: rla2rlarot !< latitude in the the rotated system1476 1477 REAL (KIND=dp) :: zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla11478 1479 1480 1481 1482 1483 1484 IF (rla > 180.0_dp) THEN1485 zrla1 = rla - 360.0_dp1486 1487 1488 1489 1490 1491 1492 1493 1494 IF (zarg2 == 0.0_dp) zarg2 = 1.0E-20_dp1495 1496 1497 1498 IF (polgam /= 0.0_dp ) THEN1499 1500 IF (rla2rlarot > 180._dp) rla2rlarot = rla2rlarot - 360.0_dp1501 1502 1503 1472 FUNCTION rla2rlarot ( phi, rla, polphi, pollam, polgam ) 1473 1474 REAL(wp), INTENT (IN) :: polphi !< latitude of the rotated north pole 1475 REAL(wp), INTENT (IN) :: pollam !< longitude of the rotated north pole 1476 REAL(wp), INTENT (IN) :: phi !< latitude in geographical system 1477 REAL(wp), INTENT (IN) :: rla !< longitude in geographical system 1478 REAL(wp), INTENT (IN) :: polgam !< angle between the north poles of the systems 1479 1480 REAL(wp) :: rla2rlarot !< latitude in the the rotated system 1481 1482 REAL(wp) :: zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1 1483 1484 zsinpol = SIN(polphi * TO_RADIANS) 1485 zcospol = COS(polphi * TO_RADIANS) 1486 zlampol = pollam * TO_RADIANS 1487 zphi = phi * TO_RADIANS 1488 1489 IF (rla > 180.0_wp) THEN 1490 zrla1 = rla - 360.0_wp 1491 ELSE 1492 zrla1 = rla 1493 ENDIF 1494 zrla = zrla1 * TO_RADIANS 1495 1496 zarg1 = - SIN (zrla-zlampol) * COS(zphi) 1497 zarg2 = - zsinpol * COS(zphi) * COS(zrla-zlampol) + zcospol * SIN(zphi) 1498 1499 IF (zarg2 == 0.0_wp) zarg2 = 1.0E-20_wp 1500 1501 rla2rlarot = ATAN2 (zarg1,zarg2) * TO_DEGREES 1502 1503 IF (polgam /= 0.0_wp ) THEN 1504 rla2rlarot = polgam + rla2rlarot 1505 IF (rla2rlarot > 180._wp) rla2rlarot = rla2rlarot - 360.0_wp 1506 ENDIF 1507 1508 END FUNCTION rla2rlarot 1504 1509 1505 1510 … … 1510 1515 !> rotated-pole system 1511 1516 !------------------------------------------------------------------------------! 1512 1513 1514 REAL(dp), INTENT (IN) :: u, v !< wind components in the true geographical system1515 REAL(dp), INTENT (IN) :: rlat, rlon !< coordinates in the true geographical system1516 REAL(dp), INTENT (IN) :: pollat, pollon !< latitude and longitude of the north pole of the rotated grid1517 1518 REAL(dp), INTENT (OUT) :: urot, vrot !< wind components in the rotated grid1519 1520 REAL (dp) :: zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm1521 1522 1523 1524 1525 1526 1527 1528 1529 znorm = 1.0_dp / SQRT(zarg1*zarg1 + zarg2*zarg2)1530 1531 1532 1533 1534 1517 SUBROUTINE uv2uvrot(u, v, rlat, rlon, pollat, pollon, urot, vrot) 1518 1519 REAL(wp), INTENT (IN) :: u, v !< wind components in the true geographical system 1520 REAL(wp), INTENT (IN) :: rlat, rlon !< coordinates in the true geographical system 1521 REAL(wp), INTENT (IN) :: pollat, pollon !< latitude and longitude of the north pole of the rotated grid 1522 1523 REAL(wp), INTENT (OUT) :: urot, vrot !< wind components in the rotated grid 1524 1525 REAL (wp) :: zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm 1526 1527 zsinpol = SIN(pollat * TO_RADIANS) 1528 zcospol = COS(pollat * TO_RADIANS) 1529 zlonp = (pollon-rlon) * TO_RADIANS 1530 zlat = rlat * TO_RADIANS 1531 1532 zarg1 = zcospol * SIN(zlonp) 1533 zarg2 = zsinpol * COS(zlat) - zcospol * SIN(zlat) * COS(zlonp) 1534 znorm = 1.0_wp / SQRT(zarg1*zarg1 + zarg2*zarg2) 1535 1536 urot = u * zarg2 * znorm - v * zarg1 * znorm 1537 vrot = u * zarg1 * znorm + v * zarg2 * znorm 1538 1539 END SUBROUTINE uv2uvrot 1535 1540 1536 1541 … … 1541 1546 !> geographical system 1542 1547 !------------------------------------------------------------------------------! 1543 1544 1545 REAL(dp), INTENT(IN) :: urot, vrot !< wind components in the rotated grid1546 REAL(dp), INTENT(IN) :: rlat, rlon !< latitude and longitude in the true geographical system1547 REAL(dp), INTENT(IN) :: pollat, pollon !< latitude and longitude of the north pole of the rotated grid1548 1549 REAL(dp), INTENT(OUT) :: u, v !< wind components in the true geographical system1550 1551 REAL(dp) :: zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm1552 1553 1554 1555 1556 1557 1558 1559 1560 znorm = 1.0_dp / SQRT(zarg1*zarg1 + zarg2*zarg2)1561 1562 1563 1564 1565 1548 SUBROUTINE uvrot2uv (urot, vrot, rlat, rlon, pollat, pollon, u, v) 1549 1550 REAL(wp), INTENT(IN) :: urot, vrot !< wind components in the rotated grid 1551 REAL(wp), INTENT(IN) :: rlat, rlon !< latitude and longitude in the true geographical system 1552 REAL(wp), INTENT(IN) :: pollat, pollon !< latitude and longitude of the north pole of the rotated grid 1553 1554 REAL(wp), INTENT(OUT) :: u, v !< wind components in the true geographical system 1555 1556 REAL(wp) :: zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm 1557 1558 zsinpol = SIN(pollat * TO_RADIANS) 1559 zcospol = COS(pollat * TO_RADIANS) 1560 zlonp = (pollon-rlon) * TO_RADIANS 1561 zlat = rlat * TO_RADIANS 1562 1563 zarg1 = zcospol * SIN(zlonp) 1564 zarg2 = zsinpol * COS(zlat) - zcospol * SIN(zlat) * COS(zlonp) 1565 znorm = 1.0_wp / SQRT(zarg1*zarg1 + zarg2*zarg2) 1566 1567 u = urot * zarg2 * znorm + vrot * zarg1 * znorm 1568 v = - urot * zarg1 * znorm + vrot * zarg2 * znorm 1569 1570 END SUBROUTINE uvrot2uv 1566 1571 1567 1572 END MODULE inifor_transform
Note: See TracChangeset
for help on using the changeset viewer.
