Changeset 3557 for palm/trunk/UTIL/inifor/src/inifor_grid.f90
- Timestamp:
- Nov 22, 2018 4:01:22 PM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/UTIL/inifor/src/inifor_grid.f90 (modified) (89 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/UTIL/inifor/src/inifor_grid.f90
r3537 r3557 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Updated documentation 29 ! 30 ! 31 ! 3537 2018-11-20 10:53:14Z eckhard 28 32 ! Read COSMO domain extents and soil depths from input files 29 33 ! Report averaging mode and debugging mode in log … … 69 73 ! Authors: 70 74 ! -------- 71 ! @author Eckhard Kadasch75 !> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach) 72 76 ! 73 77 !------------------------------------------------------------------------------! … … 193 197 INTEGER :: step_hour !< number of hours between forcing time steps 194 198 195 LOGICAL :: init_variables_required 196 LOGICAL :: boundary_variables_required 197 LOGICAL :: ls_forcing_variables_required 198 LOGICAL :: profile_grids_required 199 LOGICAL :: surface_forcing_required 199 LOGICAL :: init_variables_required !< flag controlling whether init variables are to be processed 200 LOGICAL :: boundary_variables_required !< flag controlling whether boundary grids are to be allocated and boundary variables are to be computed 201 LOGICAL :: ls_forcing_variables_required !< flag controlling whether large-scale forcing variables are to be computed 202 LOGICAL :: surface_forcing_required !< flag controlling whether surface forcing variables are to be computed 200 203 201 204 TYPE(nc_var), ALLOCATABLE, TARGET :: input_var_table(:) !< table of input variables 202 205 TYPE(nc_var), ALLOCATABLE, TARGET :: output_var_table(:) !< table of input variables 203 TYPE(nc_var) :: cosmo_var !< COSMO -DEdummy variable, used for reading HHL, rlon, rlat206 TYPE(nc_var) :: cosmo_var !< COSMO dummy variable, used for reading HHL, rlon, rlat 204 207 205 208 TYPE(grid_definition), TARGET :: palm_grid !< PALM-4U grid in the target system (COSMO-DE rotated-pole) 206 209 TYPE(grid_definition), TARGET :: palm_intermediate !< PALM-4U grid with coarse vertical grid wiht levels interpolated from COSMO-DE grid 207 210 TYPE(grid_definition), TARGET :: cosmo_grid !< target system (COSMO-DE rotated-pole) 208 TYPE(grid_definition), TARGET :: scalars_east_grid !< 209 TYPE(grid_definition), TARGET :: scalars_west_grid !< 210 TYPE(grid_definition), TARGET :: scalars_north_grid !< 211 TYPE(grid_definition), TARGET :: scalars_south_grid !< 212 TYPE(grid_definition), TARGET :: scalars_top_grid !< 213 TYPE(grid_definition), TARGET :: scalars_east_intermediate !< 214 TYPE(grid_definition), TARGET :: scalars_west_intermediate !< 215 TYPE(grid_definition), TARGET :: scalars_north_intermediate !< 216 TYPE(grid_definition), TARGET :: scalars_south_intermediate !< 217 TYPE(grid_definition), TARGET :: scalars_top_intermediate !< 218 TYPE(grid_definition), TARGET :: u_initial_grid !< 219 TYPE(grid_definition), TARGET :: u_east_grid !< 220 TYPE(grid_definition), TARGET :: u_west_grid !< 221 TYPE(grid_definition), TARGET :: u_north_grid !< 222 TYPE(grid_definition), TARGET :: u_south_grid !< 223 TYPE(grid_definition), TARGET :: u_top_grid !< 224 TYPE(grid_definition), TARGET :: u_initial_intermediate !< 225 TYPE(grid_definition), TARGET :: u_east_intermediate !< 226 TYPE(grid_definition), TARGET :: u_west_intermediate !< 227 TYPE(grid_definition), TARGET :: u_north_intermediate !< 228 TYPE(grid_definition), TARGET :: u_south_intermediate !< 229 TYPE(grid_definition), TARGET :: u_top_intermediate !< 230 TYPE(grid_definition), TARGET :: v_initial_grid !< 231 TYPE(grid_definition), TARGET :: v_east_grid !< 232 TYPE(grid_definition), TARGET :: v_west_grid !< 233 TYPE(grid_definition), TARGET :: v_north_grid !< 234 TYPE(grid_definition), TARGET :: v_south_grid !< 235 TYPE(grid_definition), TARGET :: v_top_grid !< 236 TYPE(grid_definition), TARGET :: v_initial_intermediate !< 237 TYPE(grid_definition), TARGET :: v_east_intermediate !< 238 TYPE(grid_definition), TARGET :: v_west_intermediate !< 239 TYPE(grid_definition), TARGET :: v_north_intermediate !< 240 TYPE(grid_definition), TARGET :: v_south_intermediate !< 241 TYPE(grid_definition), TARGET :: v_top_intermediate !< 242 TYPE(grid_definition), TARGET :: w_initial_grid !< 243 TYPE(grid_definition), TARGET :: w_east_grid !< 244 TYPE(grid_definition), TARGET :: w_west_grid !< 245 TYPE(grid_definition), TARGET :: w_north_grid !< 246 TYPE(grid_definition), TARGET :: w_south_grid !< 247 TYPE(grid_definition), TARGET :: w_top_grid !< 248 TYPE(grid_definition), TARGET :: w_initial_intermediate !< 249 TYPE(grid_definition), TARGET :: w_east_intermediate !< 250 TYPE(grid_definition), TARGET :: w_west_intermediate !< 251 TYPE(grid_definition), TARGET :: w_north_intermediate !< 252 TYPE(grid_definition), TARGET :: w_south_intermediate !< 253 TYPE(grid_definition), TARGET :: w_top_intermediate !< 254 255 TYPE(grid_definition), TARGET :: north_averaged_scalar_profile 256 TYPE(grid_definition), TARGET :: south_averaged_scalar_profile 257 TYPE(grid_definition), TARGET :: west_averaged_scalar_profile 258 TYPE(grid_definition), TARGET :: east_averaged_scalar_profile 259 TYPE(grid_definition), TARGET :: averaged_scalar_profile 260 TYPE(grid_definition), TARGET :: averaged_w_profile 211 TYPE(grid_definition), TARGET :: scalars_east_grid !< grid for eastern scalar boundary condition 212 TYPE(grid_definition), TARGET :: scalars_west_grid !< grid for western scalar boundary condition 213 TYPE(grid_definition), TARGET :: scalars_north_grid !< grid for northern scalar boundary condition 214 TYPE(grid_definition), TARGET :: scalars_south_grid !< grid for southern scalar boundary condition 215 TYPE(grid_definition), TARGET :: scalars_top_grid !< grid for top scalar boundary condition 216 TYPE(grid_definition), TARGET :: scalars_east_intermediate !< intermediate grid for eastern scalar boundary condition 217 TYPE(grid_definition), TARGET :: scalars_west_intermediate !< intermediate grid for western scalar boundary condition 218 TYPE(grid_definition), TARGET :: scalars_north_intermediate !< intermediate grid for northern scalar boundary condition 219 TYPE(grid_definition), TARGET :: scalars_south_intermediate !< intermediate grid for southern scalar boundary condition 220 TYPE(grid_definition), TARGET :: scalars_top_intermediate !< intermediate grid for top scalar boundary condition 221 TYPE(grid_definition), TARGET :: u_initial_grid !< grid for u initial condition 222 TYPE(grid_definition), TARGET :: u_east_grid !< grid for eastern u boundary condition 223 TYPE(grid_definition), TARGET :: u_west_grid !< grid for western u boundary condition 224 TYPE(grid_definition), TARGET :: u_north_grid !< grid for northern u boundary condition 225 TYPE(grid_definition), TARGET :: u_south_grid !< grid for southern u boundary condition 226 TYPE(grid_definition), TARGET :: u_top_grid !< grid for top u boundary condition 227 TYPE(grid_definition), TARGET :: u_initial_intermediate !< intermediate grid for u initial condition 228 TYPE(grid_definition), TARGET :: u_east_intermediate !< intermediate grid for eastern u boundary condition 229 TYPE(grid_definition), TARGET :: u_west_intermediate !< intermediate grid for western u boundary condition 230 TYPE(grid_definition), TARGET :: u_north_intermediate !< intermediate grid for northern u boundary condition 231 TYPE(grid_definition), TARGET :: u_south_intermediate !< intermediate grid for southern u boundary condition 232 TYPE(grid_definition), TARGET :: u_top_intermediate !< intermediate grid for top u boundary condition 233 TYPE(grid_definition), TARGET :: v_initial_grid !< grid for v initial condition 234 TYPE(grid_definition), TARGET :: v_east_grid !< grid for eastern v boundary condition 235 TYPE(grid_definition), TARGET :: v_west_grid !< grid for western v boundary condition 236 TYPE(grid_definition), TARGET :: v_north_grid !< grid for northern v boundary condition 237 TYPE(grid_definition), TARGET :: v_south_grid !< grid for southern v boundary condition 238 TYPE(grid_definition), TARGET :: v_top_grid !< grid for top v boundary condition 239 TYPE(grid_definition), TARGET :: v_initial_intermediate !< intermediate grid for v initial condition 240 TYPE(grid_definition), TARGET :: v_east_intermediate !< intermediate grid for eastern v boundary condition 241 TYPE(grid_definition), TARGET :: v_west_intermediate !< intermediate grid for western v boundary condition 242 TYPE(grid_definition), TARGET :: v_north_intermediate !< intermediate grid for northern v boundary condition 243 TYPE(grid_definition), TARGET :: v_south_intermediate !< intermediate grid for southern v boundary condition 244 TYPE(grid_definition), TARGET :: v_top_intermediate !< intermediate grid for top v boundary condition 245 TYPE(grid_definition), TARGET :: w_initial_grid !< grid for w initial condition 246 TYPE(grid_definition), TARGET :: w_east_grid !< grid for eastern w boundary condition 247 TYPE(grid_definition), TARGET :: w_west_grid !< grid for western w boundary condition 248 TYPE(grid_definition), TARGET :: w_north_grid !< grid for northern w boundary condition 249 TYPE(grid_definition), TARGET :: w_south_grid !< grid for southern w boundary condition 250 TYPE(grid_definition), TARGET :: w_top_grid !< grid for top w boundary condition 251 TYPE(grid_definition), TARGET :: w_initial_intermediate !< intermediate grid for w initial condition 252 TYPE(grid_definition), TARGET :: w_east_intermediate !< intermediate grid for eastern w boundary condition 253 TYPE(grid_definition), TARGET :: w_west_intermediate !< intermediate grid for western w boundary condition 254 TYPE(grid_definition), TARGET :: w_north_intermediate !< intermediate grid for northern w boundary condition 255 TYPE(grid_definition), TARGET :: w_south_intermediate !< intermediate grid for southern w boundary condition 256 TYPE(grid_definition), TARGET :: w_top_intermediate !< intermediate grid for top w boundary condition 257 TYPE(grid_definition), TARGET :: north_averaged_scalar_profile !< grid of the northern scalar averaging region 258 TYPE(grid_definition), TARGET :: south_averaged_scalar_profile !< grid of the southern scalar averaging region 259 TYPE(grid_definition), TARGET :: west_averaged_scalar_profile !< grid of the western scalar averaging region 260 TYPE(grid_definition), TARGET :: east_averaged_scalar_profile !< grid of the eastern scalar averaging region 261 TYPE(grid_definition), TARGET :: averaged_scalar_profile !< grid of the central scalar averaging region 262 TYPE(grid_definition), TARGET :: averaged_w_profile !< grid of the central w-velocity averaging region 261 263 262 264 TYPE(io_group), ALLOCATABLE, TARGET :: io_group_list(:) !< List of I/O groups, which group together output variables that share the same input variable 263 265 264 266 NAMELIST /inipar/ nx, ny, nz, dx, dy, dz, longitude, latitude, & 265 dz_max, dz_stretch_factor, dz_stretch_level, & !< old grid stretching parameters266 dz_stretch_level_start, dz_stretch_level_end !< new grid stretching parameters267 dz_max, dz_stretch_factor, dz_stretch_level, & 268 dz_stretch_level_start, dz_stretch_level_end 267 269 NAMELIST /d3par/ end_time 268 270 269 CHARACTER(LEN=LNAME) :: nc_source_text = '' !< Text describing the source of the output data, e.g. 'COSMO-DE analysis from ...' 270 CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: flow_files 271 CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: soil_moisture_files 272 CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: soil_files 273 CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: radiation_files 274 CHARACTER(LEN=SNAME) :: input_prefix !< Prefix of input files, e.g. 'laf' for COSMO-DE analyses 275 CHARACTER(LEN=SNAME) :: flow_prefix !< Prefix of flow input files, e.g. 'laf' for COSMO-DE analyses 276 CHARACTER(LEN=SNAME) :: soil_prefix !< Prefix of soil input files, e.g. 'laf' for COSMO-DE analyses 277 CHARACTER(LEN=SNAME) :: radiation_prefix !< Prefix of radiation input files, e.g. 'laf' for COSMO-DE analyses 278 CHARACTER(LEN=SNAME) :: soilmoisture_prefix !< Prefix of input files for soil moisture spin-up, e.g. 'laf' for COSMO-DE analyses 279 CHARACTER(LEN=SNAME) :: flow_suffix !< Suffix of flow input files, e.g. 'flow' 280 CHARACTER(LEN=SNAME) :: soil_suffix !< Suffix of soil input files, e.g. 'soil' 281 CHARACTER(LEN=SNAME) :: radiation_suffix !< Suffix of radiation input files, e.g. 'radiation' 282 CHARACTER(LEN=SNAME) :: soilmoisture_suffix !< Suffix of input files for soil moisture spin-up, e.g. 'soilmoisture' 271 CHARACTER(LEN=LNAME) :: nc_source_text = '' !< Text describing the source of the output data, e.g. 'COSMO-DE analysis from ...' 272 273 CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: flow_files !< list of atmospheric input files (<prefix>YYYYMMDDHH-flow.nc) 274 CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: soil_moisture_files !< list of precipitation input files (<prefix>YYYYMMDDHH-soilmoisture.nc) 275 CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: soil_files !< list of soil input files (temperature, moisture, <prefix>YYYYMMDDHH-soil.nc) 276 CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: radiation_files !< list of radiation input files (<prefix>YYYYMMDDHH-rad.nc) 277 278 CHARACTER(LEN=SNAME) :: input_prefix !< prefix of input files, e.g. 'laf' for COSMO-DE analyses 279 CHARACTER(LEN=SNAME) :: flow_prefix !< prefix of flow input files, e.g. 'laf' for COSMO-DE analyses 280 CHARACTER(LEN=SNAME) :: soil_prefix !< prefix of soil input files, e.g. 'laf' for COSMO-DE analyses 281 CHARACTER(LEN=SNAME) :: radiation_prefix !< prefix of radiation input files, e.g. 'laf' for COSMO-DE analyses 282 CHARACTER(LEN=SNAME) :: soilmoisture_prefix !< prefix of input files for soil moisture spin-up, e.g. 'laf' for COSMO-DE analyses 283 CHARACTER(LEN=SNAME) :: flow_suffix !< suffix of flow input files, e.g. 'flow' 284 CHARACTER(LEN=SNAME) :: soil_suffix !< suffix of soil input files, e.g. 'soil' 285 CHARACTER(LEN=SNAME) :: radiation_suffix !< suffix of radiation input files, e.g. 'radiation' 286 CHARACTER(LEN=SNAME) :: soilmoisture_suffix !< suffix of input files for soil moisture spin-up, e.g. 'soilmoisture' 283 287 284 TYPE(nc_file) :: output_file 285 286 TYPE(inifor_config) :: cfg !< Container of INIFORconfiguration288 TYPE(nc_file) :: output_file !< metadata of the dynamic driver 289 290 TYPE(inifor_config) :: cfg !< container of the INIFOR command-line configuration 287 291 288 292 CONTAINS 289 293 294 !------------------------------------------------------------------------------! 295 ! Description: 296 ! ------------ 297 !> This routine initializes INIFOR. This includes parsing command-line options, 298 !> setting the names of the input and output files, reading INIFOR's namelist 299 !> file as well as reading in and setting grid parameters for defining 300 !> interplation grids later in setup_grids(). 301 !------------------------------------------------------------------------------! 290 302 SUBROUTINE setup_parameters() 291 INTEGER :: k 303 INTEGER :: k !< loop index 292 304 293 305 ! … … 300 312 start_hour_soilmoisture = - (4 * 7 * 24) - 2 301 313 302 ! COSMO-DE default rotated pole 314 ! 315 !-- COSMO-DE and -D2 rotated pole position 303 316 phi_n = 40.0_dp * TO_RADIANS 304 317 phi_equat = 50.0_dp * TO_RADIANS 305 318 lambda_n = -170.0_dp * TO_RADIANS 306 319 307 ! Defaultmain centre (_c) of the PALM-4U grid in the geographical system (_g) 320 ! 321 !-- Defaultmain centre (_c) of the PALM-4U grid in the geographical system (_g) 308 322 origin_lat = 52.325079_dp * TO_RADIANS ! south-west of Berlin, origin used for the Dec 2017 showcase simulation 309 323 origin_lon = 13.082744_dp * TO_RADIANS 310 324 cfg % z0 = 35.0_dp 311 325 312 ! Default atmospheric parameters 326 ! 327 !-- Default atmospheric parameters 313 328 cfg % ug = 0.0_dp 314 329 cfg % vg = 0.0_dp 315 330 cfg % p0 = P_SL 316 331 317 ! Parameters for file names 332 ! 333 !-- Parameters for file names 318 334 start_hour_flow = 0 319 335 start_hour_soil = 0 … … 323 339 step_hour = FORCING_STEP 324 340 325 input_prefix = 'laf' ! 'laf' for COSMO-DE analyses341 input_prefix = 'laf' 326 342 cfg % flow_prefix = input_prefix 327 343 cfg % input_prefix = input_prefix … … 342 358 !------------------------------------------------------------------------------ 343 359 344 ! Set default paths and modes 360 ! 361 !-- Set default paths and modes 345 362 cfg % input_path = './' 346 363 cfg % hhl_file = '' … … 353 370 cfg % averaging_mode = 'level' 354 371 355 ! Overwrite defaults with user configuration 372 ! 373 !-- Overwrite defaults with user configuration 356 374 CALL parse_command_line_arguments( cfg ) 357 375 … … 380 398 END IF 381 399 382 383 !Set default file paths, if not specified by user.400 ! 401 !-- Set default file paths, if not specified by user. 384 402 CALL normalize_path(cfg % input_path) 385 403 IF (TRIM(cfg % hhl_file) == '') cfg % hhl_file = TRIM(cfg % input_path) // 'hhl.nc' … … 399 417 400 418 CALL run_control('time', 'init') 401 ! Read in namelist parameters 419 ! 420 !-- Read in namelist parameters 402 421 OPEN(10, FILE=cfg % namelist_file) 403 422 READ(10, NML=inipar) ! nx, ny, nz, dx, dy, dz … … 408 427 end_hour = CEILING( end_time / 3600.0 * step_hour ) 409 428 410 ! Generate input file lists 429 ! 430 !-- Generate input file lists 411 431 CALL get_input_file_list( & 412 432 cfg % start_date, start_hour_flow, end_hour, step_hour, & … … 437 457 438 458 CALL run_control('time', 'init') 439 ! Read COSMO-DE soil type map 459 ! 460 !-- Read COSMO soil type map 440 461 cosmo_var % name = 'SOILTYP' 441 462 CALL get_netcdf_variable(cfg % soiltyp_file, cosmo_var, soiltyp) … … 465 486 CALL report('setup_parameters', message) 466 487 467 468 ! Read in COSMO-DEheights of half layers (vertical cell faces)488 ! 489 !-- Read in COSMO heights of half layers (vertical cell faces) 469 490 cosmo_var % name = 'HHL' 470 491 CALL get_netcdf_variable(cfg % hhl_file, cosmo_var, hhl) … … 486 507 CALL run_control('time', 'comp') 487 508 488 ! Appoximate COSMO-DE heights of full layers (cell centres) 509 ! 510 !-- Appoximate COSMO-DE heights of full layers (cell centres) 489 511 ALLOCATE( hfl(nlon, nlat, nlev-1) ) 490 512 ALLOCATE( d_depth(ndepths), d_depth_rho_inv(ndepths) ) … … 494 516 d_depth_rho_inv = 1.0_dp / ( d_depth * RHO_L ) 495 517 496 ! Appoximate COSMO-DE heights of full layers (cell centres) 518 ! 519 !-- Appoximate COSMO-DE heights of full layers (cell centres) 497 520 DO k = 1, nlev-1 498 521 hfl(:,:,k) = 0.5_dp * ( hhl(:,:,k) + & … … 506 529 ! Section 4.2: PALM-4U parameters 507 530 !------------------------------------------------------------------------------ 508 ! PALM-4U domain extents 531 ! 532 !-- PALM-4U domain extents 509 533 lx = (nx+1) * dx 510 534 ly = (ny+1) * dy 511 535 512 ! PALM-4U point of Earth tangency 536 ! 537 !-- PALM-4U point of Earth tangency 513 538 x0 = 0.0_dp 514 539 y0 = 0.0_dp 515 540 516 ! time vector 541 ! 542 !-- time vector 517 543 nt = CEILING(end_time / (step_hour * 3600.0_dp)) + 1 518 544 ALLOCATE( time(nt) ) … … 521 547 CALL run_control('time', 'init') 522 548 523 ! Convert the PALM-4U origin coordinates to COSMO's rotated-pole grid 549 ! 550 !-- Convert the PALM-4U origin coordinates to COSMO's rotated-pole grid 524 551 phi_c = TO_RADIANS * & 525 552 phi2phirot( origin_lat * TO_DEGREES, origin_lon * TO_DEGREES,& … … 530 557 0.0_dp ) 531 558 532 ! Set gamma according to whether PALM domain is in the northern or southern 533 ! hemisphere of the COSMO-DE rotated-pole system. Gamma assumes either the 534 ! value 0 or PI and is needed to work around around a bug in the rotated-pole 535 ! coordinate transformations. 559 ! 560 !-- Set gamma according to whether PALM domain is in the northern or southern 561 !-- hemisphere of the COSMO rotated-pole system. Gamma assumes either the 562 !-- value 0 or PI and is needed to work around around a bug in the 563 !-- rotated-pole coordinate transformations. 536 564 gam = gamma_from_hemisphere(origin_lat, phi_equat) 537 565 538 ! Compute the north pole of the rotated-pole grid centred at the PALM-4U domain 539 ! centre. The resulting (phi_cn, lambda_cn) are coordinates in COSMO-DE's 540 ! rotated-pole grid. 566 ! 567 !-- Compute the north pole of the rotated-pole grid centred at the PALM-4U 568 !-- domain centre. The resulting (phi_cn, lambda_cn) are coordinates in 569 !-- COSMO-DE's rotated-pole grid. 541 570 phi_cn = phic_to_phin(phi_c) 542 571 lambda_cn = lamc_to_lamn(phi_c, lambda_c) … … 573 602 !------------------------------------------------------------------------------ 574 603 575 ! Compute coordiantes of the PALM centre in the source (COSMO) system 604 ! 605 !-- Compute coordiantes of the PALM centre in the source (COSMO) system 576 606 phi_centre = phirot2phi( & 577 607 phirot = project(0.5_dp*ly, y0, EARTH_RADIUS) * TO_DEGREES, & … … 595 625 CALL report( 'setup_parameters', message ) 596 626 597 598 !Compute boundaries of the central averaging box627 ! 628 !-- Compute boundaries of the central averaging box 599 629 averaging_angle = cfg % averaging_angle * TO_RADIANS 600 630 lam_east = lam_centre + 0.5_dp * averaging_angle … … 605 635 averaging_width_ns = averaging_angle * EARTH_RADIUS 606 636 607 ! Coriolis parameter 637 ! 638 !-- Coriolis parameter 608 639 f3 = 2.0_dp * OMEGA * SIN( & 609 640 TO_RADIANS*phirot2phi( phi_centre * TO_DEGREES, lam_centre * TO_DEGREES,& … … 618 649 ! Description: 619 650 ! ------------ 620 !> Initializes the COSMO-DE, PALM-4U, PALM-4U boundary grids. 651 !> Defines the COSMO, PALM-4U, PALM-4U boundary grids, in partucular their 652 !> coordinates and interpolation weights 621 653 !------------------------------------------------------------------------------! 622 SUBROUTINE setup_grids() ! setup_grids(inifor_settings(with nx, ny, nz,...))654 SUBROUTINE setup_grids() 623 655 CHARACTER :: interp_mode 624 656 … … 626 658 ! Section 0: Define base PALM-4U grid coordinate vectors 627 659 !------------------------------------------------------------------------------ 628 ! palm x y z, we allocate the column to nz+1 in order to include the top 629 ! scalar boundary. The interpolation grids will be associated with 630 ! a shorter column that omits the top element. 631 660 ! 661 !-- palm x y z, we allocate the column to nz+1 in order to include the top 662 !-- scalar boundary. The interpolation grids will be associated with 663 !-- a shorter column that omits the top element. 632 664 ALLOCATE( x(0:nx), y(0:ny), z(1:nz), z_column(1:nz+1) ) 633 665 CALL linspace(0.5_dp * dx, lx - 0.5_dp * dx, x) … … 642 674 z_top = z_column(nz+1) 643 675 644 ! palm xu yv zw, compared to the scalar grid, velocity coordinates 645 ! contain one element less. 676 ! 677 !-- palm xu yv zw, compared to the scalar grid, velocity coordinates 678 !-- contain one element less. 646 679 ALLOCATE( xu(1:nx), yv(1:ny), zw(1:nz-1), zw_column(1:nz)) 647 680 CALL linspace(dx, lx - dx, xu) … … 661 694 nx=nx, ny=ny, nz=nz, z=z, zw=zw, ic_mode=cfg % ic_mode) 662 695 663 ! Subtracting 1 because arrays will be allocated with nlon + 1 elements. 696 ! 697 !-- Subtracting 1 because arrays will be allocated with nlon + 1 elements. 664 698 CALL init_grid_definition('cosmo-de', grid=cosmo_grid, & 665 699 xmin=lonmin, xmax=lonmax, & … … 668 702 nx=nlon-1, ny=nlat-1, nz=nlev-1) 669 703 670 ! Define intermediate grid. This is the same as palm_grid except with a much 671 ! coarser vertical grid. The vertical levels are interpolated in each PALM-4U 672 ! column from COSMO-DE's secondary levels. The main levels are then computed as 673 ! the averages of the bounding secondary levels. 704 ! 705 !-- Define intermediate grid. This is the same as palm_grid except with a 706 !-- much coarser vertical grid. The vertical levels are interpolated in each 707 !-- PALM column from COSMO's secondary levels. The main levels are then 708 !-- computed as the averages of the bounding secondary levels. 674 709 CALL init_grid_definition('palm intermediate', grid=palm_intermediate, & 675 710 xmin=0.0_dp, xmax=lx, & … … 1007 1042 latmin = phi_south, latmax = phi_north, & 1008 1043 kind='scalar') 1009 1010 profile_grids_required = ( TRIM(cfg % ic_mode) == 'profile' .OR. &1011 TRIM(cfg % bc_mode) == 'ideal' )1012 1013 !IF (profile_grids_required) THEN1014 ! ! Here, I use the 'boundary' kind, so the vertical coordinate uses the1015 ! ! PALM z coordinate.1016 ! CALL init_grid_definition('boundary', grid=scalar_profile_grid, &1017 ! xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, &1018 ! ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, &1019 ! x0=x0, y0=y0, z0 = z0, &1020 ! nx = 0, ny = 0, nz = nz, z=z)1021 1022 ! CALL init_grid_definition('boundary', grid=w_profile_grid, &1023 ! xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, &1024 ! ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, &1025 ! x0=x0, y0=y0, z0 = z0, &1026 ! nx = 0, ny = 0, nz = nz - 1, z=zw)1027 1028 ! CALL init_grid_definition('boundary', grid=scalar_profile_intermediate,&1029 ! xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, &1030 ! ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, &1031 ! x0=x0, y0=y0, z0 = z0, &1032 ! nx = 0, ny = 0, nz = nlev - 2, z=z)1033 1034 ! CALL init_grid_definition('boundary', grid=w_profile_intermediate, &1035 ! xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, &1036 ! ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, &1037 ! x0=x0, y0=y0, z0 = z0, &1038 ! nx = 0, ny = 0, nz = nlev - 2, z=zw)1039 !END IF1040 1044 1041 1045 ! … … 1091 1095 1092 1096 1097 !------------------------------------------------------------------------------! 1098 ! Description: 1099 ! ------------ 1100 !> Driver for computing neighbour indices and weights for horizontal and 1101 !> vertical interpolation. 1102 !------------------------------------------------------------------------------! 1093 1103 SUBROUTINE setup_interpolation(cosmo_grid, grid, intermediate_grid, kind, ic_mode) 1094 1104 … … 1106 1116 ! Section 1: Horizontal interpolation 1107 1117 !------------------------------------------------------------------------------ 1108 ! Select horizontal coordinates according to kind of points (s/w, u, v) 1118 ! 1119 !-- Select horizontal coordinates according to kind of points (s/w, u, v) 1109 1120 SELECT CASE(kind) 1110 1121 1111 CASE('s') ! scalars 1122 ! 1123 !-- scalars 1124 CASE('s') 1112 1125 1113 1126 cosmo_lat => cosmo_grid % lat 1114 1127 cosmo_lon => cosmo_grid % lon 1115 1128 cosmo_h => cosmo_grid % hfl 1116 1117 CASE('w') ! vertical velocity 1129 ! 1130 !-- vertical velocity 1131 CASE('w') 1118 1132 1119 1133 cosmo_lat => cosmo_grid % lat 1120 1134 cosmo_lon => cosmo_grid % lon 1121 1135 cosmo_h => cosmo_grid % hhl 1122 1123 CASE('u') ! x velocity 1136 ! 1137 !-- x velocity 1138 CASE('u') 1124 1139 1125 1140 cosmo_lat => cosmo_grid % lat … … 1127 1142 cosmo_h => cosmo_grid % hfl 1128 1143 1129 CASE('v') ! y velocity 1144 ! 1145 !-- y velocity 1146 CASE('v') 1130 1147 1131 1148 cosmo_lat => cosmo_grid % latv … … 1153 1170 !------------------------------------------------------------------------------ 1154 1171 1155 ! If profile initialization is chosen, we--somewhat counterintuitively-- 1156 ! don't need to compute vertical interpolation weights. At least, we 1157 ! don't need them on the intermediate grid, which fills the entire PALM 1158 ! domain volume. Instead we need vertical weights for the intermediate 1159 ! profile grids, which get computed in setup_averaging(). 1160 setup_volumetric = .TRUE. 1172 ! 1173 !-- If profile initialization is chosen, we--somewhat counterintuitively-- 1174 !-- don't need to compute vertical interpolation weights. At least, we 1175 !-- don't need them on the intermediate grid, which fills the entire PALM 1176 !-- domain volume. Instead we need vertical weights for the intermediate 1177 !-- profile grids, which get computed in setup_averaging(). 1178 !-- setup_volumetric = .TRUE. 1161 1179 IF (PRESENT(ic_mode)) THEN 1162 1180 IF (TRIM(ic_mode) == 'profile') setup_volumetric = .FALSE. … … 1169 1187 intermediate_grid % h(:,:,:) = - EARTH_RADIUS 1170 1188 1171 ! For w points, use hhl, for scalars use hfl 1172 ! compute the full heights for the intermediate grids 1189 ! 1190 !-- For w points, use hhl, for scalars use hfl 1191 !-- compute the full heights for the intermediate grids 1173 1192 CALL interpolate_2d(cosmo_h, intermediate_grid % h, intermediate_grid) 1174 1193 CALL find_vertical_neighbours_and_weights_interp(grid, intermediate_grid) … … 1246 1265 grid % z => z 1247 1266 1248 ! Allocate neighbour indices and weights 1267 ! 1268 !-- Allocate neighbour indices and weights 1249 1269 IF (TRIM(ic_mode) .NE. 'profile') THEN 1250 1270 ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) ) … … 1274 1294 ) 1275 1295 1276 ! Allocate neighbour indices and weights 1296 ! 1297 !-- Allocate neighbour indices and weights 1277 1298 ALLOCATE( grid % ii(0:nx, 0:ny, 4), & 1278 1299 grid % jj(0:nx, 0:ny, 4) ) … … 1283 1304 grid % w_horiz(:,:,:) = 0.0_dp 1284 1305 1285 ! This mode initializes a Cartesian PALM-4U grid and adds the 1286 ! corresponding latitudes and longitudes of the rotated pole grid. 1306 ! 1307 !-- This mode initializes a Cartesian PALM-4U grid and adds the 1308 !-- corresponding latitudes and longitudes of the rotated pole grid. 1287 1309 CASE('palm') 1288 1310 … … 1301 1323 grid % name(3) = 'z' 1302 1324 1303 !TODO: Remove use of global dx, dy, dz variables. Consider 1304 !TODO: associating global x,y, and z arrays. 1325 ! 1326 !-- TODO: Remove use of global dx, dy, dz variables. Consider 1327 !-- TODO: associating global x,y, and z arrays. 1305 1328 ALLOCATE( grid % x(0:nx), grid % y(0:ny) ) 1306 1329 ALLOCATE( grid % xu(1:nx), grid % yv(1:ny) ) … … 1314 1337 grid % depths => depths 1315 1338 1316 ! Allocate neighbour indices and weights 1339 ! 1340 !-- Allocate neighbour indices and weights 1317 1341 IF (TRIM(ic_mode) .NE. 'profile') THEN 1318 1342 ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) ) … … 1329 1353 grid % name(3) = 'interpolated hhl or hfl' 1330 1354 1331 !TODO: Remove use of global dx, dy, dz variables. Consider 1332 !TODO: associating global x,y, and z arrays. 1355 ! 1356 !-- TODO: Remove use of global dx, dy, dz variables. Consider 1357 !-- TODO: associating global x,y, and z arrays. 1333 1358 ALLOCATE( grid % x(0:nx), grid % y(0:ny) ) 1334 1359 ALLOCATE( grid % xu(1:nx), grid % yv(1:ny) ) … … 1340 1365 grid % depths => depths 1341 1366 1342 ! Allocate rotated-pole coordinates, clon is for (c)osmo-de (lon)gitude 1367 ! 1368 !-- Allocate rotated-pole coordinates, clon is for (c)osmo-de (lon)gitude 1343 1369 ALLOCATE( grid % clon(0:nx, 0:ny), grid % clat(0:nx, 0:ny) ) 1344 1370 ALLOCATE( grid % clonu(1:nx, 0:ny), grid % clatu(1:nx, 0:ny) ) 1345 1371 ALLOCATE( grid % clonv(0:nx, 1:ny), grid % clatv(0:nx, 1:ny) ) 1346 1372 1347 ! Compute rotated-pole coordinates of... 1348 ! ... PALM-4U centres 1373 ! 1374 !-- Compute rotated-pole coordinates of... 1375 !-- ... PALM-4U centres 1349 1376 CALL rotate_to_cosmo( & 1350 1377 phir = project( grid % y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude … … 1356 1383 ) 1357 1384 1358 ! ... PALM-4U u winds 1385 ! 1386 !-- ... PALM-4U u winds 1359 1387 CALL rotate_to_cosmo( & 1360 1388 phir = project( grid % y, y0, EARTH_RADIUS ), & ! = plate-carree latitude … … 1366 1394 ) 1367 1395 1368 ! ... PALM-4U v winds 1396 ! 1397 !-- ... PALM-4U v winds 1369 1398 CALL rotate_to_cosmo( & 1370 1399 phir = project( grid % yv, y0, EARTH_RADIUS ), & ! = plate-carree latitude … … 1376 1405 ) 1377 1406 1378 ! Allocate neighbour indices and weights 1407 ! 1408 !-- Allocate neighbour indices and weights 1379 1409 ALLOCATE( grid % ii(0:nx, 0:ny, 4), & 1380 1410 grid % jj(0:nx, 0:ny, 4) ) … … 1398 1428 grid % latv(:) = grid % lat + 0.5_dp * (grid % ly / grid % ny) 1399 1429 1400 ! Point to heights of half levels (hhl) and compute heights of full 1401 ! levels (hfl) as arithmetic averages 1430 ! 1431 !-- Point to heights of half levels (hhl) and compute heights of full 1432 !-- levels (hfl) as arithmetic averages 1402 1433 grid % hhl => hhl 1403 1434 grid % hfl => hfl … … 1472 1503 avg_grid % kind = TRIM(kind) 1473 1504 1474 ! Find and store COSMO columns that fall into the coordinate range 1475 ! given by avg_grid % clon, % clat 1505 ! 1506 !-- Find and store COSMO columns that fall into the coordinate range 1507 !-- given by avg_grid % clon, % clat 1476 1508 CALL get_cosmo_averaging_region(avg_grid, cosmo_grid) 1477 1509 1478 1510 ALLOCATE (avg_grid % kkk(avg_grid % n_columns, avg_grid % nz, 2) ) 1479 1511 ALLOCATE (avg_grid % w(avg_grid % n_columns, avg_grid % nz, 2) ) 1480 ! Compute average COSMO levels in the averaging region 1512 ! 1513 !-- Compute average COSMO levels in the averaging region 1481 1514 SELECT CASE(avg_grid % kind) 1482 1515 … … 1494 1527 END SELECT 1495 1528 1496 ! For level-besed averaging, compute average heights 1529 ! 1530 !-- For level-besed averaging, compute average heights 1497 1531 level_based_averaging = .TRUE. 1498 1532 IF (level_based_averaging) THEN … … 1503 1537 END IF 1504 1538 1505 ! Copy averaged grid to all COSMO columns, leads to computing the same 1506 ! vertical interpolation weights for all columns and to COSMO grid level 1507 ! based averaging onto the averaged COSMO heights. 1539 ! 1540 !-- Copy averaged grid to all COSMO columns, leads to computing the same 1541 !-- vertical interpolation weights for all columns and to COSMO grid level 1542 !-- based averaging onto the averaged COSMO heights. 1508 1543 IF ( TRIM(cfg % averaging_mode) == 'level' ) THEN 1509 1544 FORALL( & … … 1513 1548 END IF 1514 1549 1515 ! Compute vertical weights and neighbours 1550 ! 1551 !-- Compute vertical weights and neighbours 1516 1552 CALL find_vertical_neighbours_and_weights_average(avg_grid) 1517 1553 … … 1550 1586 END SELECT 1551 1587 1552 1553 !FIXME: make dlon, dlat parameters of the grid_defintion type1588 ! 1589 !-- FIXME: make dlon, dlat parameters of the grid_defintion type 1554 1590 dlon = cosmo_lon(1) - cosmo_lon(0) 1555 1591 dlat = cosmo_lat(1) - cosmo_lat(0) … … 1610 1646 REAL(dp) :: dz_level_end, dz_stretched 1611 1647 1612 INTEGER :: dz_stretch_level_end_index(9) 1613 INTEGER :: dz_stretch_level_start_index(9) 1648 INTEGER :: dz_stretch_level_end_index(9) !< vertical grid level index until which the vertical grid spacing is stretched 1649 INTEGER :: dz_stretch_level_start_index(9) !< vertical grid level index above which the vertical grid spacing is stretched 1614 1650 INTEGER :: dz_stretch_level_index = 0 1615 1651 INTEGER :: k, n, number_dz 1652 1616 1653 ! 1617 1654 !-- Compute height of u-levels from constant grid length and dz stretch factors … … 1826 1863 1827 1864 1865 !------------------------------------------------------------------------------! 1828 1866 ! Description: [PALM subroutine] 1829 1867 ! -----------------------------------------------------------------------------! … … 1936 1974 1937 1975 ! 1938 !-- 1939 !-- 1940 !-- 1941 !-- 1976 !-- stretch_factor_1 is taken to guarantee that the stretching 1977 !-- procedure ends as close as possible to dz_stretch_level_end. 1978 !-- stretch_factor_2 would guarantee that the stretched dz(n) is 1979 !-- equal to dz(n+1) after l_rounded grid levels. 1942 1980 IF (delta_total_new < delta_total_old) THEN 1943 1981 dz_stretch_factor_array(n) = stretch_factor_1 … … 1974 2012 END SUBROUTINE calculate_stretching_factor 1975 2013 2014 !------------------------------------------------------------------------------! 2015 ! Description: 2016 ! ------------ 2017 !> Computes the midpoint between two neighbouring coordinates in the given 2018 !> coordinate vector 'z' and stores it in 'zw'. 2019 !------------------------------------------------------------------------------! 1976 2020 SUBROUTINE midpoints(z, zw) 1977 2021 … … 1987 2031 END SUBROUTINE midpoints 1988 2032 1989 2033 !------------------------------------------------------------------------------! 2034 ! Description: 2035 ! ------------ 2036 !> Defines INFOR's IO groups. 2037 !------------------------------------------------------------------------------! 1990 2038 SUBROUTINE setup_io_groups() 1991 2039 … … 1994 2042 ngroups = 16 1995 2043 ALLOCATE( io_group_list(ngroups) ) 1996 1997 !soil temp 2044 2045 ! 2046 !-- soil temp 1998 2047 io_group_list(1) = init_io_group( & 1999 2048 in_files = soil_files, & … … 2003 2052 ) 2004 2053 2005 !soil water 2054 ! 2055 !-- soil water 2006 2056 io_group_list(2) = init_io_group( & 2007 2057 in_files = soil_files, & … … 2011 2061 ) 2012 2062 2013 !potential temperature, surface pressure, specific humidity including 2014 !nudging and subsidence, and geostrophic winds ug, vg 2063 ! 2064 !-- potential temperature, surface pressure, specific humidity including 2065 !-- nudging and subsidence, and geostrophic winds ug, vg 2015 2066 io_group_list(3) = init_io_group( & 2016 2067 in_files = flow_files, & … … 2025 2076 ) 2026 2077 2027 !Moved to therodynamic io_group 2078 ! 2079 !-- Moved to therodynamic io_group 2028 2080 !io_group_list(4) = init_io_group( & 2029 2081 ! in_files = flow_files, & … … 2033 2085 !) 2034 2086 2035 !u and v velocity, ug anv vg moved to thermodynamic io group. 2087 ! 2088 !-- u and v velocity 2036 2089 io_group_list(5) = init_io_group( & 2037 2090 in_files = flow_files, & … … 2043 2096 ) 2044 2097 2045 !!v velocity, deprecated! 2098 ! 2099 !-- v velocity, deprecated! 2046 2100 !io_group_list(6) = init_io_group( & 2047 2101 ! in_files = flow_files, & … … 2052 2106 !io_group_list(6) % to_be_processed = .FALSE. 2053 2107 2054 !w velocity and subsidence and w nudging 2108 ! 2109 !-- w velocity and subsidence and w nudging 2055 2110 io_group_list(7) = init_io_group( & 2056 2111 in_files = flow_files, & … … 2059 2114 kind = 'scalar' & 2060 2115 ) 2061 2062 !rain2116 ! 2117 !-- rain 2063 2118 io_group_list(8) = init_io_group( & 2064 2119 in_files = soil_moisture_files, & … … 2068 2123 ) 2069 2124 io_group_list(8) % to_be_processed = .FALSE. 2070 2071 !snow2125 ! 2126 !-- snow 2072 2127 io_group_list(9) = init_io_group( & 2073 2128 in_files = soil_moisture_files, & … … 2077 2132 ) 2078 2133 io_group_list(9) % to_be_processed = .FALSE. 2079 2080 !graupel2134 ! 2135 !-- graupel 2081 2136 io_group_list(10) = init_io_group( & 2082 2137 in_files = soil_moisture_files, & … … 2086 2141 ) 2087 2142 io_group_list(10) % to_be_processed = .FALSE. 2088 2089 !evapotranspiration2143 ! 2144 !-- evapotranspiration 2090 2145 io_group_list(11) = init_io_group( & 2091 2146 in_files = soil_moisture_files, & … … 2095 2150 ) 2096 2151 io_group_list(11) % to_be_processed = .FALSE. 2097 2098 !2m air temperature2152 ! 2153 !-- 2m air temperature 2099 2154 io_group_list(12) = init_io_group( & 2100 2155 in_files = soil_moisture_files, & … … 2103 2158 kind = 'surface' & 2104 2159 ) 2105 2106 !incoming diffusive sw flux2160 ! 2161 !-- incoming diffusive sw flux 2107 2162 io_group_list(13) = init_io_group( & 2108 2163 in_files = radiation_files, & … … 2112 2167 ) 2113 2168 io_group_list(13) % to_be_processed = .FALSE. 2114 2115 !incoming direct sw flux2169 ! 2170 !-- incoming direct sw flux 2116 2171 io_group_list(14) = init_io_group( & 2117 2172 in_files = radiation_files, & … … 2121 2176 ) 2122 2177 io_group_list(14) % to_be_processed = .FALSE. 2123 2124 !sw radiation balance2178 ! 2179 !-- sw radiation balance 2125 2180 io_group_list(15) = init_io_group( & 2126 2181 in_files = radiation_files, & … … 2130 2185 ) 2131 2186 io_group_list(15) % to_be_processed = .FALSE. 2132 2133 !lw radiation balance2187 ! 2188 !-- lw radiation balance 2134 2189 io_group_list(16) = init_io_group( & 2135 2190 in_files = radiation_files, & … … 2143 2198 2144 2199 2200 !------------------------------------------------------------------------------! 2201 ! Description: 2202 ! ------------ 2203 !> Initializes an IO group with the list of input files, the lists of 2204 !> input and output variables, and sets the number of output quantities. 2205 !> 2206 !> In this context, output quantities refers to the physical quantities required 2207 !> to interpolate netCDF output variables, not the output variables itsleft. For 2208 !> instance, the output variables ls_forcing_left/_right/_north/_south all rely 2209 !> on the output quantity Theta. 2210 !------------------------------------------------------------------------------! 2145 2211 FUNCTION init_io_group(in_files, out_vars, in_var_list, kind, & 2146 2212 n_output_quantities) RESULT(group) … … 2157 2223 group % n_inputs = SIZE(in_var_list) 2158 2224 group % kind = TRIM(kind) 2159 2225 ! 2226 !-- For the 'thermodynamics' IO group, one quantity more than input variables 2227 !-- is needed to compute all output variables of the IO group. Concretely, in 2228 !-- preprocess() the density is computed from T,P or PP,QV in adddition to 2229 !-- the variables Theta, p, qv. In read_input_variables(), 2230 !-- n_output_quantities is used to allocate the correct number of input 2231 !-- buffers. 2160 2232 IF ( PRESENT(n_output_quantities) ) THEN 2161 2233 group % n_output_quantities = n_output_quantities … … 2207 2279 ! Description: 2208 2280 ! ------------ 2209 !> Initializes the thevariable list.2281 !> Initializes the variable list. 2210 2282 !------------------------------------------------------------------------------! 2211 2283 SUBROUTINE setup_variable_tables(ic_mode) … … 2787 2859 intermediate_grid = palm_intermediate & 2788 2860 ) 2789 ! Radiation fluxes and balances 2861 2790 2862 output_var_table(38) = init_nc_var( & 2791 2863 name = 'rad_swd_dif_0', & … … 3154 3226 output_var_table(64) % to_be_processed = .NOT. cfg % ug_is_set 3155 3227 3156 ! Attributes shared among all variables 3228 ! 3229 !-- Attributes shared among all variables 3157 3230 output_var_table(:) % source = nc_source_text 3158 3231 … … 3164 3237 ! Description: 3165 3238 ! ------------ 3166 !> Initializes an dnc_var varible with the given parameters. The 'kind'3239 !> Initializes an nc_var varible with the given parameters. The 'kind' 3167 3240 !> parameter is used to infer the correct netCDF IDs and the level of detail, 3168 3241 !> 'lod', as defined by the PALM-4U input data standard. … … 3199 3272 SELECT CASE( TRIM(out_var_kind) ) 3200 3273 3201 !TODO: Using global module variables 'init_variables_required' and 3202 !TODO: 'boundary_variables_required'. Encapsulate in settings type 3203 !TODO: and pass into init_nc_var. 3274 ! 3275 !-- TODO: Using global module variables 'init_variables_required' and 3276 !-- TODO: 'boundary_variables_required'. Encapsulate in settings type 3277 !-- TODO: and pass into init_nc_var. 3204 3278 CASE( 'init soil' ) 3205 3279 var % nt = 1 … … 3591 3665 3592 3666 CASE( 'velocities' ) 3593 ! Allocate a compute buffer with the same number of arrays as the input 3667 ! 3668 !-- Allocate a compute buffer with the same number of arrays as the input 3594 3669 ALLOCATE( preprocess_buffer( SIZE(input_buffer) ) ) 3595 3670 3596 ! Allocate u and v arrays with scalar dimensions 3671 ! 3672 !-- Allocate u and v arrays with scalar dimensions 3597 3673 nx = SIZE(input_buffer(1) % array, 1) 3598 3674 ny = SIZE(input_buffer(1) % array, 2) … … 3603 3679 CALL run_control('time', 'alloc') 3604 3680 3605 ! interpolate U and V to centres 3681 ! 3682 !-- interpolate U and V to centres 3606 3683 CALL centre_velocities( u_face = input_buffer(1) % array, & 3607 3684 v_face = input_buffer(2) % array, & … … 3613 3690 3614 3691 CASE('rotated-pole') 3615 ! rotate U and V to PALM-4U orientation and overwrite U and V with 3616 ! rotated velocities 3692 ! 3693 !-- rotate U and V to PALM-4U orientation and overwrite U and V with 3694 !-- rotated velocities 3617 3695 DO k = 1, nz 3618 3696 DO j = 1, ny … … 3637 3715 END SELECT 3638 3716 3639 ! set values3640 3717 input_buffer(1) % array(1,:,:) = 0.0_dp 3641 3718 input_buffer(2) % array(1,:,:) = 0.0_dp … … 3658 3735 nz = SIZE(input_buffer(1) % array, 3) 3659 3736 3660 ! Compute absolute pressure if presure perturbation has been read in. 3737 ! 3738 !-- Compute absolute pressure if presure perturbation has been read in. 3661 3739 IF ( TRIM(group % in_var_list(2) % name) == 'PP' ) THEN 3662 3740 message = "Absolute pressure, P, not available, " // & … … 3673 3751 RD, G, basic_state_pressure) 3674 3752 3753 ! 3754 !-- Overwrite pressure perturbation with absolute pressure. HECTO 3755 !-- converts pressure perturbation from hPa to Pa. 3675 3756 input_buffer (2) % array(i,j,:) = & 3676 3757 HECTO * input_buffer (2) % array(i,j,:) + & … … 3687 3768 3688 3769 END IF 3689 ! pressure 3770 ! 3771 !-- mark pressure as preprocessed 3690 3772 input_buffer(2) % is_preprocessed = .TRUE. 3691 3773 3692 ! Copy temperature to last input buffer array 3774 ! 3775 !-- Copy temperature to the last input buffer array 3693 3776 ALLOCATE( & 3694 3777 input_buffer( group % n_output_quantities ) % array (nx, ny, nz) & … … 3697 3780 input_buffer(1) % array(:,:,:) 3698 3781 3699 ! Convert absolute to potential temperature 3782 ! 3783 !-- Convert absolute in place to potential temperature 3700 3784 CALL potential_temperature( & 3701 3785 t = input_buffer(1) % array(:,:,:), & … … 3706 3790 ) 3707 3791 3708 ! potential temperature 3792 ! 3793 !-- mark potential temperature as preprocessed 3709 3794 input_buffer(1) % is_preprocessed = .TRUE. 3710 3795 3711 ! Convert temperature copy to density 3796 ! 3797 !-- Convert temperature copy to density 3712 3798 CALL moist_density( & 3713 3799 t_rho = input_buffer(group % n_output_quantities) % array(:,:,:), & … … 3718 3804 ) 3719 3805 3720 ! qv 3806 ! 3807 !-- mark qv as preprocessed 3721 3808 input_buffer(3) % is_preprocessed = .TRUE. 3722 3809 3723 ! density 3810 ! 3811 !-- mark density as preprocessed 3724 3812 input_buffer(group % n_output_quantities) % is_preprocessed = .TRUE. 3725 3813 … … 3780 3868 SELECT CASE(dt) 3781 3869 3870 ! 3871 !-- input has been accumulated over one hour. Leave as is 3872 !-- input_buffer(1) % array(:,:,:) carrries one-hour integral 3782 3873 CASE(1) 3783 !input has been accumulated over one hour. Leave as is 3784 !input_buffer(1) % array(:,:,:) carrries one-hour integral 3785 3874 3875 ! 3876 !-- input has been accumulated over two hours. Subtract previous step 3877 !-- input_buffer(1) % array(:,:,:) carrries one-hour integral 3878 !-- input_buffer(2) % array(:,:,:) carrries two-hour integral 3786 3879 CASE(2) 3787 !input has been accumulated over two hours. Subtract previous step3788 !input_buffer(1) % array(:,:,:) carrries one-hour integral3789 !input_buffer(2) % array(:,:,:) carrries two-hour integral3790 3880 CALL deaverage( & 3791 3881 avg_1 = input_buffer(1) % array(:,:,:), t1 = 1.0_dp, & 3792 3882 avg_2 = input_buffer(2) % array(:,:,:), t2 = 1.0_dp, & 3793 3883 avg_3 = input_buffer(1) % array(:,:,:), t3 = 1.0_dp ) 3794 !input_buffer(1) % array(:,:,:) carrries one-hour integral of second hour 3795 3884 ! 3885 !-- input_buffer(1) % array(:,:,:) carrries one-hour integral of second hour 3886 3887 ! 3888 !-- input has been accumulated over three hours. Subtract previous step 3889 !-- input_buffer(1) % array(:,:,:) carrries three-hour integral 3890 !-- input_buffer(2) % array(:,:,:) still carrries two-hour integral 3796 3891 CASE(3) 3797 !input has been accumulated over three hours. Subtract previous step3798 !input_buffer(1) % array(:,:,:) carrries three-hour integral3799 !input_buffer(2) % array(:,:,:) still carrries two-hour integral3800 3892 CALL deaverage( & 3801 3893 avg_1 = input_buffer(2) % array(:,:,:), t1 = 1.0_dp, & 3802 3894 avg_2 = input_buffer(1) % array(:,:,:), t2 = 1.0_dp, & 3803 3895 avg_3 = input_buffer(1) % array(:,:,:), t3 = 1.0_dp ) 3804 !input_buffer(1) % array(:,:,:) carrries one-hour integral of third hourA 3896 ! 3897 !-- input_buffer(1) % array(:,:,:) carrries one-hour integral of third hourA 3805 3898 3806 3899 CASE DEFAULT … … 3818 3911 3819 3912 hour = iter - 1 3820 dt = MODULO(hour, 3) + 1 ! averaging period 3913 ! 3914 !-- averaging period 3915 dt = MODULO(hour, 3) + 1 3821 3916 SELECT CASE(dt) 3822 3917 ! 3918 !-- input has been accumulated over one hour. Leave as is 3919 !-- input_buffer(1) % array(:,:,:) carrries one-hour integral 3823 3920 CASE(1) 3824 !input has been accumulated over one hour. Leave as is 3825 !input_buffer(1) % array(:,:,:) carrries one-hour integral 3826 3921 3922 ! 3923 !-- input has been accumulated over two hours. Subtract previous step 3924 !-- input_buffer(1) % array(:,:,:) carrries one-hour integral 3925 !-- input_buffer(2) % array(:,:,:) carrries two-hour integral 3827 3926 CASE(2) 3828 !input has been accumulated over two hours. Subtract previous step3829 !input_buffer(1) % array(:,:,:) carrries one-hour integral3830 !input_buffer(2) % array(:,:,:) carrries two-hour integral3831 3927 CALL deaverage( input_buffer(1) % array(:,:,:), 1.0_dp, & 3832 3928 input_buffer(2) % array(:,:,:), 2.0_dp, & 3833 3929 input_buffer(1) % array(:,:,:), 1.0_dp) 3834 !input_buffer(1) % array(:,:,:) carrries one-hour integral of second hour 3835 3930 ! 3931 !-- input_buffer(1) % array(:,:,:) carrries one-hour integral of second hour 3932 3933 ! 3934 !-- input has been accumulated over three hours. Subtract previous step 3935 !-- input_buffer(1) % array(:,:,:) carrries three-hour integral 3936 !-- input_buffer(2) % array(:,:,:) still carrries two-hour integral 3836 3937 CASE(3) 3837 !input has been accumulated over three hours. Subtract previous step3838 !input_buffer(1) % array(:,:,:) carrries three-hour integral3839 !input_buffer(2) % array(:,:,:) still carrries two-hour integral3840 3938 CALL deaverage( input_buffer(2) % array(:,:,:), 2.0_dp, & 3841 3939 input_buffer(1) % array(:,:,:), 3.0_dp, & 3842 3940 input_buffer(1) % array(:,:,:), 1.0_dp) 3843 !input_buffer(1) % array(:,:,:) carrries one-hour integral of third hourA 3941 ! 3942 !-- input_buffer(1) % array(:,:,:) carrries one-hour integral of third hourA 3844 3943 3845 3944 CASE DEFAULT … … 3860 3959 3861 3960 3961 !------------------------------------------------------------------------------! 3862 3962 ! Description: 3863 3963 ! ------------ … … 3928 4028 END DO 3929 4029 3930 ! Overwrite if at least one non-water neighbour cell is available 4030 ! 4031 !-- Overwrite if at least one non-water neighbour cell is available 3931 4032 IF (n_cells > 0) THEN 3932 4033 array(i,j,:) = column(:) / n_cells
Note: See TracChangeset
for help on using the changeset viewer.
