Changeset 3395 for palm/trunk/UTIL/inifor/src/grid.f90
- Timestamp:
- Oct 22, 2018 5:32:49 PM (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/UTIL/inifor/src/grid.f90
r3183 r3395 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Added computation of geostrophic winds form COSMO pressure fields 29 ! Introduced averaging grids and internal 'output' variables for computation of 30 ! geostrophic and large-scale forcings 31 ! Merged qv, uv and vg output variables and 'temperature' IO group into new 32 ! 'thermodynamics' IP group 33 ! Cleaned up terminal output, show some messages only with --debug option 34 ! Improved variable naming and readability 35 ! 36 ! 37 ! 3183 2018-07-27 14:25:55Z suehring 28 38 ! Introduced new PALM grid stretching 29 39 ! Updated variable names and metadata for PIDS v1.9 compatibility … … 61 71 ONLY: DATE, EARTH_RADIUS, TO_RADIANS, TO_DEGREES, PI, dp, hp, sp, & 62 72 SNAME, LNAME, PATH, FORCING_STEP, WATER_ID, FILL_ITERATIONS, & 63 BETA, P_SL, T_SL, BETA, RD, G, P_REF, RD_PALM, CP_PALM, RHO_L 73 BETA, P_SL, T_SL, BETA, RD, RV, G, P_REF, RD_PALM, CP_PALM, & 74 RHO_L, OMEGA, HECTO 64 75 USE io, & 65 76 ONLY: get_netcdf_variable, get_netcdf_attribute, & … … 68 79 ONLY: NF90_MAX_NAME, NF90_MAX_VAR_DIMS 69 80 USE transform, & 70 ONLY: rotate_to_cosmo, find_horizontal_neighbours,&81 ONLY: average_2d, rotate_to_cosmo, find_horizontal_neighbours,& 71 82 compute_horizontal_interp_weights, & 72 find_vertical_neighbours_and_weights, interpolate_2d, & 73 gamma_from_hemisphere, phic_to_phin, lamc_to_lamn, average_2d, & 74 project, centre_velocities, phi2phirot, rla2rlarot, uv2uvrot 83 find_vertical_neighbours_and_weights_interp, & 84 find_vertical_neighbours_and_weights_average, interpolate_2d, & 85 gamma_from_hemisphere, phic_to_phin, lamc_to_lamn, & 86 project, centre_velocities, phi2phirot, rla2rlarot, uv2uvrot, & 87 phirot2phi, rlarot2rla 75 88 USE types 76 89 USE util … … 80 93 SAVE 81 94 95 REAL(dp) :: averaging_angle = 0.0_dp !< latitudal and longitudal width of averaging regions [rad] 96 REAL(dp) :: averaging_width_ns = 0.0_dp !< longitudal width of averaging regions [m] 97 REAL(dp) :: averaging_width_ew = 0.0_dp !< latitudal width of averaging regions [m] 82 98 REAL(dp) :: phi_equat = 0.0_dp !< latitude of rotated equator of COSMO-DE grid [rad] 83 99 REAL(dp) :: phi_n = 0.0_dp !< latitude of rotated pole of COSMO-DE grid [rad] … … 87 103 REAL(dp) :: phi_cn = 0.0_dp !< latitude of the rotated pole relative to the COSMO-DE grid [rad] 88 104 REAL(dp) :: lambda_cn = 0.0_dp !< longitude of the rotated pole relative to the COSMO-DE grid [rad] 105 REAL(dp) :: lam_centre = 0.0_dp !< longitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad] 106 REAL(dp) :: phi_centre = 0.0_dp !< latitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad] 107 REAL(dp) :: lam_east = 0.0_dp !< longitude of the east central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad] 108 REAL(dp) :: lam_west = 0.0_dp !< longitude of the west central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad] 109 REAL(dp) :: phi_north = 0.0_dp !< latitude of the north central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad] 110 REAL(dp) :: phi_south = 0.0_dp !< latitude of the south central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad] 89 111 REAL(dp) :: gam = 0.0_dp !< angle for working around phirot2phi/rlarot2rla bug 90 112 REAL(dp) :: dx = 0.0_dp !< PALM-4U grid spacing in x direction [m] … … 100 122 REAL(dp) :: dyi = 0.0_dp !< inverse PALM-4U grid spacing in y direction [m^-1] 101 123 REAL(dp) :: dzi = 0.0_dp !< inverse PALM-4U grid spacing in z direction [m^-1] 124 REAL(dp) :: f3 = 0.0_dp !< Coriolis parameter 102 125 REAL(dp) :: lx = 0.0_dp !< PALM-4U domain size in x direction [m] 103 126 REAL(dp) :: ly = 0.0_dp !< PALM-4U domain size in y direction [m] 104 REAL(dp) :: ug = 0.0_dp !< geostrophic wind in x direction [m/s]105 REAL(dp) :: vg = 0.0_dp !< geostrophic wind in y direction [m/s]106 127 REAL(dp) :: p0 = 0.0_dp !< PALM-4U surface pressure, at z0 [Pa] 107 128 REAL(dp) :: x0 = 0.0_dp !< x coordinate of PALM-4U Earth tangent [m] … … 140 161 INTEGER :: dz_stretch_level_end_index(9) !< vertical grid level index until which the vertical grid spacing is stretched 141 162 INTEGER :: dz_stretch_level_start_index(9) !< vertical grid level index above which the vertical grid spacing is stretched 142 INTEGER :: i !< indexing variable143 INTEGER :: average_imin, average_imax, average_jmin, average_jmax !< index ranges for profile averaging144 INTEGER :: k !< indexing variable145 163 INTEGER :: nt !< number of output time steps 146 164 INTEGER :: nx !< number of PALM-4U grid points in x direction … … 164 182 LOGICAL :: profile_grids_required 165 183 166 INTEGER :: n_invar = 0 !< number of variables in the input variable table167 INTEGER :: n_outvar = 0 !< number of variables in the output variable table168 184 TYPE(nc_var), ALLOCATABLE, TARGET :: input_var_table(:) !< table of input variables 169 185 TYPE(nc_var), ALLOCATABLE, TARGET :: output_var_table(:) !< table of input variables … … 219 235 TYPE(grid_definition), TARGET :: w_south_intermediate !< 220 236 TYPE(grid_definition), TARGET :: w_top_intermediate !< 221 TYPE(grid_definition), TARGET :: scalar_profile_grid !< 222 TYPE(grid_definition), TARGET :: scalar_profile_intermediate !< 223 TYPE(grid_definition), TARGET :: w_profile_grid !< 224 TYPE(grid_definition), TARGET :: w_profile_intermediate !< 237 238 TYPE(grid_definition), TARGET :: north_averaged_scalar_profile 239 TYPE(grid_definition), TARGET :: south_averaged_scalar_profile 240 TYPE(grid_definition), TARGET :: west_averaged_scalar_profile 241 TYPE(grid_definition), TARGET :: east_averaged_scalar_profile 242 TYPE(grid_definition), TARGET :: averaged_scalar_profile 243 TYPE(grid_definition), TARGET :: averaged_w_profile 225 244 226 245 TYPE(io_group), ALLOCATABLE, TARGET :: io_group_list(:) !< List of I/O groups, which group together output variables that share the same input variable … … 237 256 CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: radiation_files 238 257 CHARACTER(LEN=SNAME) :: input_prefix !< Prefix of input files, e.g. 'laf' for COSMO-DE analyses 258 CHARACTER(LEN=SNAME) :: flow_prefix !< Prefix of flow input files, e.g. 'laf' for COSMO-DE analyses 259 CHARACTER(LEN=SNAME) :: soil_prefix !< Prefix of soil input files, e.g. 'laf' for COSMO-DE analyses 260 CHARACTER(LEN=SNAME) :: radiation_prefix !< Prefix of radiation input files, e.g. 'laf' for COSMO-DE analyses 261 CHARACTER(LEN=SNAME) :: soilmoisture_prefix !< Prefix of input files for soil moisture spin-up, e.g. 'laf' for COSMO-DE analyses 239 262 CHARACTER(LEN=SNAME) :: flow_suffix !< Suffix of flow input files, e.g. 'flow' 240 263 CHARACTER(LEN=SNAME) :: soil_suffix !< Suffix of soil input files, e.g. 'soil' … … 249 272 250 273 SUBROUTINE setup_parameters() 274 INTEGER :: k 251 275 252 276 ! … … 282 306 283 307 ! Default atmospheric parameters 284 ug = 0.0_dp285 vg = 0.0_dp308 cfg % ug = 0.0_dp 309 cfg % vg = 0.0_dp 286 310 cfg % p0 = P_SL 287 311 … … 296 320 input_prefix = 'laf' ! 'laf' for COSMO-DE analyses 297 321 cfg % flow_prefix = input_prefix 322 cfg % input_prefix = input_prefix 298 323 cfg % soil_prefix = input_prefix 299 324 cfg % radiation_prefix = input_prefix … … 304 329 radiation_suffix = '-rad' 305 330 soilmoisture_suffix = '-soilmoisture' 331 332 cfg % debug = .FALSE. 333 cfg % averaging_angle = 2.0_dp 306 334 ! 307 335 !------------------------------------------------------------------------------ … … 318 346 cfg % ic_mode = 'volume' 319 347 cfg % bc_mode = 'real' 320 321 ! Update default file names and domain centre 348 cfg % averaging_mode = 'level' 349 350 ! Overwrite defaults with user configuration 322 351 CALL parse_command_line_arguments( cfg ) 352 353 flow_prefix = TRIM(cfg % input_prefix) 354 radiation_prefix = TRIM(cfg % input_prefix) 355 soil_prefix = TRIM(cfg % input_prefix) 356 soilmoisture_prefix = TRIM(cfg % input_prefix) 357 IF (cfg % flow_prefix_is_set) flow_prefix = TRIM(cfg % flow_prefix) 358 IF (cfg % radiation_prefix_is_set) radiation_prefix = TRIM(cfg % radiation_prefix) 359 IF (cfg % soil_prefix_is_set) soil_prefix = TRIM(cfg % soil_prefix) 360 IF (cfg % soilmoisture_prefix_is_set) soilmoisture_prefix = TRIM(cfg % soilmoisture_prefix) 323 361 324 362 output_file % name = cfg % output_file … … 336 374 END IF 337 375 376 377 ! Set default file paths, if not specified by user. 338 378 CALL normalize_path(cfg % input_path) 339 379 IF (TRIM(cfg % hhl_file) == '') cfg % hhl_file = TRIM(cfg % input_path) // 'hhl.nc' … … 361 401 362 402 ! Generate input file lists 363 CALL get_input_file_list(cfg % start_date, start_hour_flow, end_hour, step_hour, & 364 cfg % input_path, cfg % flow_prefix, flow_suffix, flow_files) 365 CALL get_input_file_list(cfg % start_date, start_hour_soil, end_hour, step_hour, & 366 cfg % input_path, cfg % soil_prefix, soil_suffix, soil_files) 367 CALL get_input_file_list(cfg % start_date, start_hour_radiation, end_hour, step_hour, & 368 cfg % input_path, cfg % radiation_prefix, radiation_suffix, radiation_files) 369 CALL get_input_file_list(cfg % start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, & 370 cfg % input_path, cfg % soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files) 403 CALL get_input_file_list( & 404 cfg % start_date, start_hour_flow, end_hour, step_hour, & 405 cfg % input_path, flow_prefix, flow_suffix, flow_files) 406 CALL get_input_file_list( & 407 cfg % start_date, start_hour_soil, end_hour, step_hour, & 408 cfg % input_path, soil_prefix, soil_suffix, soil_files) 409 CALL get_input_file_list( & 410 cfg % start_date, start_hour_radiation, end_hour, step_hour, & 411 cfg % input_path, radiation_prefix, radiation_suffix, radiation_files) 412 CALL get_input_file_list( & 413 cfg % start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, & 414 cfg % input_path, soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files) 371 415 372 416 ! … … 499 543 500 544 CALL run_control('time', 'comp') 545 546 !------------------------------------------------------------------------------ 547 ! Section 4.3: INIFOR averaging domains 548 !------------------------------------------------------------------------------ 549 550 ! Compute coordiantes of the PALM centre in the source (COSMO) system 551 phi_centre = phirot2phi( & 552 phirot = project(0.5_dp*ly, y0, EARTH_RADIUS) * TO_DEGREES, & 553 rlarot = project(0.5_dp*lx, x0, EARTH_RADIUS) * TO_DEGREES, & 554 polphi = phi_cn * TO_DEGREES, pollam = lambda_cn * TO_DEGREES, & 555 polgam = gam * TO_DEGREES & 556 ) * TO_RADIANS 557 558 lam_centre = rlarot2rla( & 559 phirot = project(0.5_dp*ly, y0, EARTH_RADIUS) * TO_DEGREES, & 560 rlarot = project(0.5_dp*lx, x0, EARTH_RADIUS) * TO_DEGREES, & 561 polphi = phi_cn * TO_DEGREES, pollam = lambda_cn * TO_DEGREES, & 562 polgam = gam * TO_DEGREES & 563 ) * TO_RADIANS 564 565 message = "PALM-4U centre:" // NEW_LINE('') // & 566 " lon (lambda) = " // & 567 TRIM(real_to_str_f(lam_centre * TO_DEGREES)) // " deg" // NEW_LINE(' ') // & 568 " lat (phi ) = " // & 569 TRIM(real_to_str_f(phi_centre * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)" 570 CALL report( 'setup_parameters', message ) 571 572 573 ! Compute boundaries of the central averaging box 574 averaging_angle = cfg % averaging_angle * TO_RADIANS 575 lam_east = lam_centre + 0.5_dp * averaging_angle 576 lam_west = lam_centre - 0.5_dp * averaging_angle 577 phi_north = phi_centre + 0.5_dp * averaging_angle 578 phi_south = phi_centre - 0.5_dp * averaging_angle 579 averaging_width_ew = averaging_angle * COS(phi_centre) * EARTH_RADIUS 580 averaging_width_ns = averaging_angle * EARTH_RADIUS 581 582 ! Coriolis parameter 583 f3 = 2.0_dp * OMEGA * SIN( & 584 TO_RADIANS*phirot2phi( phi_centre * TO_DEGREES, lam_centre * TO_DEGREES,& 585 phi_n * TO_DEGREES, lambda_n * TO_DEGREES, & 586 gam * TO_DEGREES ) & 587 ) 501 588 502 589 END SUBROUTINE setup_parameters … … 603 690 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 604 691 x0=x0, y0=y0, z0 = z0, & 605 nx = nx, ny = ny, nz = nlev - 2)692 nx = nx, ny = ny, nz = nlev - 1) 606 693 607 694 IF (boundary_variables_required) THEN … … 824 911 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 825 912 x0=x0, y0=y0, z0 = z0, & 826 nx = 0, ny = ny, nz = nlev - 2)913 nx = 0, ny = ny, nz = nlev - 1) 827 914 828 915 CALL init_grid_definition('boundary intermediate', grid=w_west_intermediate, & … … 830 917 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 831 918 x0=x0, y0=y0, z0 = z0, & 832 nx = 0, ny = ny, nz = nlev - 2)919 nx = 0, ny = ny, nz = nlev - 1) 833 920 834 921 CALL init_grid_definition('boundary intermediate', grid=w_north_intermediate, & … … 836 923 ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & 837 924 x0=x0, y0=y0, z0 = z0, & 838 nx = nx, ny = 0, nz = nlev - 2)925 nx = nx, ny = 0, nz = nlev - 1) 839 926 840 927 CALL init_grid_definition('boundary intermediate', grid=w_south_intermediate, & … … 842 929 ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & 843 930 x0=x0, y0=y0, z0 = z0, & 844 nx = nx, ny = 0, nz = nlev - 2)931 nx = nx, ny = 0, nz = nlev - 1) 845 932 846 933 CALL init_grid_definition('boundary intermediate', grid=w_top_intermediate, & … … 848 935 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 849 936 x0=x0, y0=y0, z0 = z0, & 850 nx = nx, ny = ny, nz = nlev - 2)937 nx = nx, ny = ny, nz = nlev - 1) 851 938 END IF 852 939 … … 856 943 !------------------------------------------------------------------------------ 857 944 858 profile_grids_required = ( TRIM(cfg % ic_mode) == 'profile' .OR. & 945 CALL init_averaging_grid(averaged_scalar_profile, cosmo_grid, & 946 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, & 947 lonmin = lam_west, lonmax = lam_east, & 948 latmin = phi_south, latmax = phi_north, & 949 kind='scalar') 950 951 CALL init_averaging_grid(averaged_w_profile, cosmo_grid, & 952 x = 0.5_dp * lx, y = 0.5_dp * ly, z = zw, z0 = z0, & 953 lonmin = lam_west, lonmax = lam_east, & 954 latmin = phi_south, latmax = phi_north, & 955 kind='w') 956 957 CALL init_averaging_grid(south_averaged_scalar_profile, cosmo_grid, & 958 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, & 959 lonmin = lam_west, lonmax = lam_east, & 960 latmin = phi_centre - averaging_angle, & 961 latmax = phi_centre, & 962 kind='scalar') 963 964 CALL init_averaging_grid(north_averaged_scalar_profile, cosmo_grid, & 965 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, & 966 lonmin = lam_west, lonmax = lam_east, & 967 latmin = phi_centre, & 968 latmax = phi_centre + averaging_angle, & 969 kind='scalar') 970 971 CALL init_averaging_grid(west_averaged_scalar_profile, cosmo_grid, & 972 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, & 973 lonmin = lam_centre - averaging_angle, & 974 lonmax = lam_centre, & 975 latmin = phi_south, latmax = phi_north, & 976 kind='scalar') 977 978 CALL init_averaging_grid(east_averaged_scalar_profile, cosmo_grid, & 979 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, & 980 lonmin = lam_centre, & 981 lonmax = lam_centre + averaging_angle, & 982 latmin = phi_south, latmax = phi_north, & 983 kind='scalar') 984 985 profile_grids_required = ( TRIM(cfg % ic_mode) == 'profile' .OR. & 859 986 TRIM(cfg % bc_mode) == 'ideal' ) 860 987 861 IF (profile_grids_required) THEN 862 CALL init_grid_definition('boundary', grid=scalar_profile_grid, & 863 xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 864 ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 865 x0=x0, y0=y0, z0 = z0, & 866 nx = 0, ny = 0, nz = nz, z=z) 867 868 CALL init_grid_definition('boundary', grid=w_profile_grid, & 869 xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 870 ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 871 x0=x0, y0=y0, z0 = z0, & 872 nx = 0, ny = 0, nz = nz - 1, z=zw) 873 874 CALL init_grid_definition('boundary', grid=scalar_profile_intermediate,& 875 xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 876 ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 877 x0=x0, y0=y0, z0 = z0, & 878 nx = 0, ny = 0, nz = nlev - 2, z=z) 879 880 CALL init_grid_definition('boundary', grid=w_profile_intermediate, & 881 xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 882 ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 883 x0=x0, y0=y0, z0 = z0, & 884 nx = 0, ny = 0, nz = nlev - 2, z=zw) 885 END IF 988 !IF (profile_grids_required) THEN 989 ! ! Here, I use the 'boundary' kind, so the vertical coordinate uses the 990 ! ! PALM z coordinate. 991 ! CALL init_grid_definition('boundary', grid=scalar_profile_grid, & 992 ! xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 993 ! ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 994 ! x0=x0, y0=y0, z0 = z0, & 995 ! nx = 0, ny = 0, nz = nz, z=z) 996 997 ! CALL init_grid_definition('boundary', grid=w_profile_grid, & 998 ! xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 999 ! ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 1000 ! x0=x0, y0=y0, z0 = z0, & 1001 ! nx = 0, ny = 0, nz = nz - 1, z=zw) 1002 1003 ! CALL init_grid_definition('boundary', grid=scalar_profile_intermediate,& 1004 ! xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 1005 ! ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 1006 ! x0=x0, y0=y0, z0 = z0, & 1007 ! nx = 0, ny = 0, nz = nlev - 2, z=z) 1008 1009 ! CALL init_grid_definition('boundary', grid=w_profile_intermediate, & 1010 ! xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 1011 ! ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 1012 ! x0=x0, y0=y0, z0 = z0, & 1013 ! nx = 0, ny = 0, nz = nlev - 2, z=zw) 1014 !END IF 886 1015 887 1016 ! … … 930 1059 931 1060 IF (TRIM(cfg % ic_mode) == 'profile') THEN 932 CALL setup_averaging(cosmo_grid, palm_intermediate, & 933 average_imin, average_imax, average_jmin, average_jmax) 1061 !TODO: remove this conditional if not needed. 934 1062 END IF 935 1063 … … 938 1066 939 1067 940 SUBROUTINE setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, kind, ic_mode)1068 SUBROUTINE setup_interpolation(cosmo_grid, grid, intermediate_grid, kind, ic_mode) 941 1069 942 1070 TYPE(grid_definition), INTENT(IN), TARGET :: cosmo_grid 943 TYPE(grid_definition), INTENT(INOUT), TARGET :: palm_grid, palm_intermediate1071 TYPE(grid_definition), INTENT(INOUT), TARGET :: grid, intermediate_grid 944 1072 CHARACTER, INTENT(IN) :: kind 945 1073 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: ic_mode 946 1074 947 TYPE(grid_definition), POINTER :: grid 948 REAL(dp), DIMENSION(:), POINTER :: lat, lon 949 REAL(dp), DIMENSION(:,:,:), POINTER :: h 950 951 LOGICAL :: setup_vertical 952 953 setup_vertical = .TRUE. 954 IF (PRESENT(ic_mode)) THEN 955 IF (TRIM(ic_mode) == 'profile') setup_vertical = .FALSE. 956 END IF 1075 REAL(dp), DIMENSION(:), POINTER :: cosmo_lat, cosmo_lon 1076 REAL(dp), DIMENSION(:,:,:), POINTER :: cosmo_h 1077 1078 LOGICAL :: setup_volumetric 957 1079 958 1080 !------------------------------------------------------------------------------ … … 964 1086 CASE('s') ! scalars 965 1087 966 lat => cosmo_grid % lat967 lon => cosmo_grid % lon968 h => cosmo_grid % hfl1088 cosmo_lat => cosmo_grid % lat 1089 cosmo_lon => cosmo_grid % lon 1090 cosmo_h => cosmo_grid % hfl 969 1091 970 1092 CASE('w') ! vertical velocity 971 1093 972 lat => cosmo_grid % lat973 lon => cosmo_grid % lon974 h => cosmo_grid % hhl1094 cosmo_lat => cosmo_grid % lat 1095 cosmo_lon => cosmo_grid % lon 1096 cosmo_h => cosmo_grid % hhl 975 1097 976 1098 CASE('u') ! x velocity 977 1099 978 lat => cosmo_grid % lat979 lon => cosmo_grid % lonu980 h => cosmo_grid % hfl1100 cosmo_lat => cosmo_grid % lat 1101 cosmo_lon => cosmo_grid % lonu 1102 cosmo_h => cosmo_grid % hfl 981 1103 982 1104 CASE('v') ! y velocity 983 1105 984 lat => cosmo_grid % latv985 lon => cosmo_grid % lon986 h => cosmo_grid % hfl1106 cosmo_lat => cosmo_grid % latv 1107 cosmo_lon => cosmo_grid % lon 1108 cosmo_h => cosmo_grid % hfl 987 1109 988 1110 CASE DEFAULT … … 993 1115 END SELECT 994 1116 995 grid => palm_intermediate 996 997 CALL find_horizontal_neighbours(lat, lon, & 998 grid % clat, grid % clon, grid % ii, grid % jj) 999 1000 CALL compute_horizontal_interp_weights(lat, lon, & 1001 grid % clat, grid % clon, grid % ii, grid % jj, grid % w_horiz) 1117 CALL find_horizontal_neighbours(cosmo_lat, cosmo_lon, & 1118 intermediate_grid % clat, intermediate_grid % clon, & 1119 intermediate_grid % ii, intermediate_grid % jj) 1120 1121 CALL compute_horizontal_interp_weights(cosmo_lat, cosmo_lon, & 1122 intermediate_grid % clat, intermediate_grid % clon, & 1123 intermediate_grid % ii, intermediate_grid % jj, & 1124 intermediate_grid % w_horiz) 1002 1125 1003 1126 !------------------------------------------------------------------------------ … … 1005 1128 !------------------------------------------------------------------------------ 1006 1129 1007 IF (setup_vertical) THEN 1008 ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) ) 1009 grid % h(:,:,:) = - EARTH_RADIUS 1130 ! If profile initialization is chosen, we--somewhat counterintuitively-- 1131 ! don't need to compute vertical interpolation weights. At least, we 1132 ! don't need them on the intermediate grid, which fills the entire PALM 1133 ! domain volume. Instead we need vertical weights for the intermediate 1134 ! profile grids, which get computed in setup_averaging(). 1135 setup_volumetric = .TRUE. 1136 IF (PRESENT(ic_mode)) THEN 1137 IF (TRIM(ic_mode) == 'profile') setup_volumetric = .FALSE. 1138 END IF 1139 1140 IF (setup_volumetric) THEN 1141 ALLOCATE( intermediate_grid % h(0:intermediate_grid % nx, & 1142 0:intermediate_grid % ny, & 1143 0:intermediate_grid % nz) ) 1144 intermediate_grid % h(:,:,:) = - EARTH_RADIUS 1010 1145 1011 1146 ! For w points, use hhl, for scalars use hfl 1012 1147 ! compute the full heights for the intermediate grids 1013 CALL interpolate_2d(cosmo_ grid % hfl, grid % h,grid)1014 CALL find_vertical_neighbours_and_weights (palm_grid,grid)1148 CALL interpolate_2d(cosmo_h, intermediate_grid % h, intermediate_grid) 1149 CALL find_vertical_neighbours_and_weights_interp(grid, intermediate_grid) 1015 1150 END IF 1016 1151 1017 1152 END SUBROUTINE setup_interpolation 1018 1019 1020 SUBROUTINE setup_averaging(cosmo_grid, palm_intermediate, imin, imax, jmin, jmax)1021 1022 TYPE(grid_definition), INTENT(IN) :: cosmo_grid, palm_intermediate1023 INTEGER, INTENT(INOUT) :: imin, imax, jmin, jmax1024 1025 TYPE(grid_definition), POINTER :: grid1026 REAL(dp) :: lonmin_pos,lonmax_pos, latmin_pos, latmax_pos1027 REAL(dp) :: cosmo_dxi, cosmo_dyi1028 1029 cosmo_dxi = 1.0_dp / (cosmo_grid % lon(1) - cosmo_grid % lon(0))1030 cosmo_dyi = 1.0_dp / (cosmo_grid % lat(1) - cosmo_grid % lat(0))1031 1032 ! find horizontal index ranges for profile averaging1033 lonmin_pos = (MINVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi1034 lonmax_pos = (MAXVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi1035 latmin_pos = (MINVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi1036 latmax_pos = (MAXVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi1037 1038 imin = FLOOR(lonmin_pos)1039 imax = CEILING(lonmax_pos)1040 jmin = FLOOR(latmin_pos)1041 jmax = CEILING(latmax_pos)1042 1043 ! average heights for intermediate scalar and w profile grids1044 grid => scalar_profile_intermediate1045 ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) )1046 grid % h(:,:,:) = - EARTH_RADIUS1047 CALL average_2d(cosmo_grid % hfl, grid % h(0,0,:), imin, imax, jmin, jmax)1048 CALL find_vertical_neighbours_and_weights(scalar_profile_grid, grid)1049 1050 grid => w_profile_intermediate1051 ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) )1052 grid % h(:,:,:) = - EARTH_RADIUS1053 CALL average_2d(cosmo_grid % hhl, grid % h(0,0,:), imin, imax, jmin, jmax)1054 CALL find_vertical_neighbours_and_weights(w_profile_grid, grid)1055 1056 END SUBROUTINE setup_averaging1057 1058 1153 1059 1154 !------------------------------------------------------------------------------! … … 1291 1386 1292 1387 END SUBROUTINE init_grid_definition 1388 1389 1390 !------------------------------------------------------------------------------! 1391 ! Description: 1392 ! ------------ 1393 !> Initializes averagin_grid-type variables. 1394 !> 1395 !> Input parameters: 1396 !> ----------------- 1397 !> 1398 !> cosmo_grid : grid_definition-type varieble initialized as 'cosmo-de' grid 1399 !> providing longitudes and latitudes of the parent grid 1400 !> 1401 !> x, y : location of the profile to be averaged in the PALM system [m] 1402 !> 1403 !> z : vertical grid coordinates of the profile in the PALM system [m] 1404 !> 1405 !> <lat/lon><min/max>: boundaries of the averaging region in the parent system, 1406 !> i.e. the COSMO-DE rotated-pole system. [rad] 1407 !> 1408 !> kind : kind of quantity to be averaged using this averaging grid. 1409 !> Destinguishes COSMO-DE scalar and w-velocity levels. Note that finding the 1410 !> parent/COSMO columns for the region in get_cosmo_averaging_region() is 1411 !> independent of 'kind' b/ecause INIFOR uses column-centred u and v velocity 1412 !> components, which are computed in the preprocessing step. 1413 !> 1414 !> Output parameters: 1415 !> ------------------ 1416 !> avg_grid : averaging_grid variable to be initialized. 1417 !------------------------------------------------------------------------------! 1418 SUBROUTINE init_averaging_grid(avg_grid, cosmo_grid, x, y, z, z0, & 1419 lonmin, lonmax, latmin, latmax, kind) 1420 1421 TYPE(grid_definition), INTENT(INOUT) :: avg_grid 1422 TYPE(grid_definition), INTENT(IN) :: cosmo_grid 1423 REAL(dp), INTENT(IN) :: x, y, z0 1424 REAL(dp), INTENT(IN), TARGET :: z(:) 1425 REAL(dp), INTENT(IN) :: lonmin, lonmax, latmin, latmax 1426 1427 CHARACTER(LEN=*), INTENT(IN) :: kind 1428 1429 LOGICAL :: level_based_averaging 1430 1431 INTEGER :: i, j 1432 1433 ALLOCATE( avg_grid % x(1) ) 1434 ALLOCATE( avg_grid % y(1) ) 1435 avg_grid % x(1) = x 1436 avg_grid % y(1) = y 1437 avg_grid % z => z 1438 avg_grid % z0 = z0 1439 1440 avg_grid % nz = SIZE(z, 1) 1441 1442 ALLOCATE( avg_grid % lon(2) ) 1443 ALLOCATE( avg_grid % lat(2) ) 1444 avg_grid % lon(1:2) = (/lonmin, lonmax/) 1445 avg_grid % lat(1:2) = (/latmin, latmax/) 1446 1447 avg_grid % kind = TRIM(kind) 1448 1449 ! Find and store COSMO columns that fall into the coordinate range 1450 ! given by avg_grid % clon, % clat 1451 CALL get_cosmo_averaging_region(avg_grid, cosmo_grid) 1452 1453 ALLOCATE (avg_grid % kkk(avg_grid % n_columns, avg_grid % nz, 2) ) 1454 ALLOCATE (avg_grid % w(avg_grid % n_columns, avg_grid % nz, 2) ) 1455 ! Compute average COSMO levels in the averaging region 1456 SELECT CASE(avg_grid % kind) 1457 1458 CASE('scalar', 'u', 'v') 1459 avg_grid % cosmo_h => cosmo_grid % hfl 1460 1461 CASE('w') 1462 avg_grid % cosmo_h => cosmo_grid % hhl 1463 1464 CASE DEFAULT 1465 message = "Averaging grid kind '" // TRIM(avg_grid % kind) // & 1466 "' is not supported. Use 'scalar', 'u', or 'v'." 1467 CALL abort('get_cosmo_averaging_region', message) 1468 1469 END SELECT 1470 1471 ! For level-besed averaging, compute average heights 1472 level_based_averaging = .TRUE. 1473 IF (level_based_averaging) THEN 1474 ALLOCATE(avg_grid % h(1,1,SIZE(avg_grid % cosmo_h, 3)) ) 1475 1476 CALL average_2d(avg_grid % cosmo_h, avg_grid % h(1,1,:), & 1477 avg_grid % iii, avg_grid % jjj) 1478 END IF 1479 1480 ! Copy averaged grid to all COSMO columnts, leads to computing the same 1481 ! vertical interpolation weights for all columns and to COSMO grid level 1482 ! based averaging onto the averaged COSMO heights. 1483 IF ( TRIM(cfg % averaging_mode) == 'level' ) THEN 1484 FORALL( & 1485 i=1:SIZE(avg_grid % cosmo_h, 1), & 1486 j=1:SIZE(avg_grid % cosmo_h, 2) & 1487 ) avg_grid % cosmo_h(i,j,:) = avg_grid % h(1,1,:) 1488 END IF 1489 1490 ! Compute vertical weights and neighbours 1491 CALL find_vertical_neighbours_and_weights_average(avg_grid) 1492 1493 END SUBROUTINE init_averaging_grid 1494 1495 1496 SUBROUTINE get_cosmo_averaging_region(avg_grid, cosmo_grid) 1497 TYPE(grid_definition), INTENT(INOUT) :: avg_grid 1498 TYPE(grid_definition), TARGET, INTENT(IN) :: cosmo_grid 1499 1500 REAL(dp), DIMENSION(:), POINTER :: cosmo_lon, cosmo_lat 1501 REAL(dp) :: dlon, dlat 1502 1503 INTEGER :: i, j, imin, imax, jmin, jmax, l, nx, ny 1504 1505 1506 SELECT CASE( TRIM(avg_grid % kind) ) 1507 1508 CASE('scalar', 'w') 1509 cosmo_lon => cosmo_grid % lon 1510 cosmo_lat => cosmo_grid % lat 1511 1512 CASE('u') 1513 cosmo_lon => cosmo_grid % lonu 1514 cosmo_lat => cosmo_grid % lat 1515 1516 CASE('v') 1517 cosmo_lon => cosmo_grid % lon 1518 cosmo_lat => cosmo_grid % latv 1519 1520 CASE DEFAULT 1521 message = "Averaging grid kind '" // TRIM(avg_grid % kind) // & 1522 "' is not supported. Use 'scalar', 'u', or 'v'." 1523 CALL abort('get_cosmo_averaging_region', message) 1524 1525 END SELECT 1526 1527 1528 ! FIXME: make dlon, dlat parameters of the grid_defintion type 1529 dlon = cosmo_lon(1) - cosmo_lon(0) 1530 dlat = cosmo_lat(1) - cosmo_lat(0) 1531 1532 imin = CEILING( (avg_grid % lon(1) - cosmo_lon(0)) / dlon ) 1533 imax = FLOOR( (avg_grid % lon(2) - cosmo_lon(0)) / dlon ) 1534 1535 jmin = CEILING( (avg_grid % lat(1) - cosmo_lat(0)) / dlat ) 1536 jmax = FLOOR( (avg_grid % lat(2) - cosmo_lat(0)) / dlat ) 1537 1538 message = "Averaging over "// & 1539 TRIM(str(imin)) // " < i < " // TRIM(str(imax)) // & 1540 " and " // & 1541 TRIM(str(jmin)) // " < j < " // TRIM(str(jmax)) 1542 CALL report( 'get_cosmo_averaging_region', message ) 1543 1544 nx = imax - imin + 1 1545 ny = jmax - jmin + 1 1546 avg_grid % n_columns = nx * ny 1547 1548 ALLOCATE( avg_grid % iii(avg_grid % n_columns), & 1549 avg_grid % jjj(avg_grid % n_columns) ) 1550 1551 l = 0 1552 DO j = jmin, jmax 1553 DO i = imin, imax 1554 l = l + 1 1555 avg_grid % iii(l) = i 1556 avg_grid % jjj(l) = j 1557 END DO 1558 END DO 1559 1560 END SUBROUTINE get_cosmo_averaging_region 1293 1561 1294 1562 … … 1718 1986 ) 1719 1987 1720 !potential temperature, surface pressure, including nudging and subsidence 1988 !potential temperature, surface pressure, specific humidity including 1989 !nudging and subsidence, and geostrophic winds ug, vg 1721 1990 io_group_list(3) = init_io_group( & 1722 1991 in_files = flow_files, & 1723 out_vars = [output_var_table(3:8), output_var_table(42:42), & 1724 output_var_table(49:51)], & 1725 in_var_list = (/input_var_table(3), input_var_table(17)/), & 1726 kind = 'temperature' & 1727 ) 1728 1729 !specific humidity including nudging and subsidence 1730 io_group_list(4) = init_io_group( & 1731 in_files = flow_files, & 1732 out_vars = [output_var_table(9:14), output_var_table(52:54)], & 1733 in_var_list = input_var_table(4:4), & 1734 kind = 'scalar' & 1735 ) 1736 1737 !u and v velocity and geostrophic winds ug, vg 1992 out_vars = [output_var_table(56:64), & ! internal averaged density and pressure profiles 1993 output_var_table(3:8), output_var_table(9:14), & 1994 output_var_table(42:42), output_var_table(43:44), & 1995 output_var_table(49:51), output_var_table(52:54)], & 1996 in_var_list = (/input_var_table(3), input_var_table(17), & ! T, P, QV 1997 input_var_table(4) /), & 1998 kind = 'thermodynamics', & 1999 n_output_quantities = 4 & ! P, Theta, Rho, qv 2000 ) 2001 2002 !Moved to therodynamic io_group 2003 !io_group_list(4) = init_io_group( & 2004 ! in_files = flow_files, & 2005 ! out_vars = [output_var_table(9:14), output_var_table(52:54)], & 2006 ! in_var_list = input_var_table(4:4), & 2007 ! kind = 'scalar' & 2008 !) 2009 2010 !u and v velocity, ug anv vg moved to thermodynamic io group. 1738 2011 io_group_list(5) = init_io_group( & 1739 2012 in_files = flow_files, & 1740 out_vars = [output_var_table(15:26), output_var_table(4 3:46)], &2013 out_vars = [output_var_table(15:26), output_var_table(45:46)], & 1741 2014 !out_vars = output_var_table(15:20), & 1742 2015 in_var_list = input_var_table(5:6), & … … 1845 2118 1846 2119 1847 FUNCTION init_io_group(in_files, out_vars, in_var_list, kind) RESULT(group) 2120 FUNCTION init_io_group(in_files, out_vars, in_var_list, kind, & 2121 n_output_quantities) RESULT(group) 1848 2122 CHARACTER(LEN=PATH), INTENT(IN) :: in_files(:) 1849 2123 CHARACTER(LEN=*), INTENT(IN) :: kind 1850 2124 TYPE(nc_var), INTENT(IN) :: out_vars(:) 1851 2125 TYPE(nc_var), INTENT(IN) :: in_var_list(:) 2126 INTEGER, OPTIONAL :: n_output_quantities 1852 2127 1853 2128 TYPE(io_group) :: group … … 1855 2130 group % nt = SIZE(in_files) 1856 2131 group % nv = SIZE(out_vars) 2132 group % n_inputs = SIZE(in_var_list) 1857 2133 group % kind = TRIM(kind) 1858 2134 1859 ALLOCATE(group % in_var_list( SIZE(in_var_list) )) 2135 IF ( PRESENT(n_output_quantities) ) THEN 2136 group % n_output_quantities = n_output_quantities 2137 ELSE 2138 group % n_output_quantities = group % n_inputs 2139 END IF 2140 2141 ALLOCATE(group % in_var_list(group % n_inputs)) 1860 2142 ALLOCATE(group % in_files(group % nt)) 1861 2143 ALLOCATE(group % out_vars(group % nv)) … … 1876 2158 SUBROUTINE fini_grids() 1877 2159 1878 CALL report('fini_grids', 'Deallocating grids' )2160 CALL report('fini_grids', 'Deallocating grids', cfg % debug) 1879 2161 1880 2162 DEALLOCATE(x, y, z, xu, yv, zw, z_column, zw_column) … … 1904 2186 SUBROUTINE setup_variable_tables(ic_mode) 1905 2187 CHARACTER(LEN=*), INTENT(IN) :: ic_mode 2188 INTEGER :: n_invar = 0 !< number of variables in the input variable table 2189 INTEGER :: n_outvar = 0 !< number of variables in the output variable table 1906 2190 TYPE(nc_var), POINTER :: var 1907 2191 … … 1914 2198 1915 2199 n_invar = 17 1916 n_outvar = 552200 n_outvar = 64 1917 2201 ALLOCATE( input_var_table(n_invar) ) 1918 2202 ALLOCATE( output_var_table(n_outvar) ) … … 2052 2336 ) 2053 2337 IF (TRIM(ic_mode) == 'profile') THEN 2054 output_var_table(3) % grid => scalar_profile_grid 2055 output_var_table(3) % intermediate_grid => scalar_profile_intermediate 2338 output_var_table(3) % averaging_grid => averaged_scalar_profile 2056 2339 END IF 2057 2340 … … 2122 2405 units = "kg/kg", & 2123 2406 kind = "init scalar", & 2124 input_id = 1, &2407 input_id = 3, & 2125 2408 output_file = output_file, & 2126 2409 grid = palm_grid, & … … 2129 2412 ) 2130 2413 IF (TRIM(ic_mode) == 'profile') THEN 2131 output_var_table(9) % grid => scalar_profile_grid 2132 output_var_table(9) % intermediate_grid => scalar_profile_intermediate 2414 output_var_table(9) % averaging_grid => averaged_scalar_profile 2133 2415 END IF 2134 2416 … … 2139 2421 units = "kg/kg", & 2140 2422 kind = "left scalar", & 2141 input_id = 1, &2423 input_id = 3, & 2142 2424 output_file = output_file, & 2143 2425 grid = scalars_west_grid, & … … 2151 2433 units = "kg/kg", & 2152 2434 kind = "right scalar", & 2153 input_id = 1, &2435 input_id = 3, & 2154 2436 output_file = output_file, & 2155 2437 grid = scalars_east_grid, & … … 2163 2445 units = "kg/kg", & 2164 2446 kind = "north scalar", & 2165 input_id = 1, &2447 input_id = 3, & 2166 2448 output_file = output_file, & 2167 2449 grid = scalars_north_grid, & … … 2175 2457 units = "kg/kg", & 2176 2458 kind = "south scalar", & 2177 input_id = 1, &2459 input_id = 3, & 2178 2460 output_file = output_file, & 2179 2461 grid = scalars_south_grid, & … … 2187 2469 units = "kg/kg", & 2188 2470 kind = "top scalar", & 2189 input_id = 1, &2471 input_id = 3, & 2190 2472 output_file = output_file, & 2191 2473 grid = scalars_top_grid, & … … 2206 2488 ) 2207 2489 IF (TRIM(ic_mode) == 'profile') THEN 2208 output_var_table(15) % grid => scalar_profile_grid 2209 output_var_table(15) % intermediate_grid => scalar_profile_intermediate 2490 output_var_table(15) % averaging_grid => averaged_scalar_profile 2210 2491 END IF 2211 2492 … … 2283 2564 ) 2284 2565 IF (TRIM(ic_mode) == 'profile') THEN 2285 output_var_table(21) % grid => scalar_profile_grid 2286 output_var_table(21) % intermediate_grid => scalar_profile_intermediate 2566 output_var_table(21) % averaging_grid => averaged_scalar_profile 2287 2567 END IF 2288 2568 … … 2360 2640 ) 2361 2641 IF (TRIM(ic_mode) == 'profile') THEN 2362 output_var_table(27) % grid => w_profile_grid 2363 output_var_table(27) % intermediate_grid => w_profile_intermediate 2642 output_var_table(27) % averaging_grid => averaged_w_profile 2364 2643 END IF 2365 2644 … … 2546 2825 intermediate_grid = palm_intermediate & 2547 2826 ) 2827 output_var_table(42) % averaging_grid => averaged_scalar_profile 2548 2828 2549 2829 output_var_table(43) = init_nc_var( & … … 2552 2832 long_name = "geostrophic wind (u component)", & 2553 2833 units = "m/s", & 2554 kind = " constant scalar profile",&2834 kind = "geostrophic", & 2555 2835 input_id = 1, & 2556 2836 output_file = output_file, & 2557 grid = scalar_profile_grid,&2558 intermediate_grid = scalar_profile_intermediate&2837 grid = averaged_scalar_profile, & 2838 intermediate_grid = averaged_scalar_profile & 2559 2839 ) 2560 2840 … … 2564 2844 long_name = "geostrophic wind (v component)", & 2565 2845 units = "m/s", & 2566 kind = " constant scalar profile",&2846 kind = "geostrophic", & 2567 2847 input_id = 1, & 2568 2848 output_file = output_file, & 2569 grid = scalar_profile_grid,&2570 intermediate_grid = scalar_profile_intermediate&2849 grid = averaged_scalar_profile, & 2850 intermediate_grid = averaged_scalar_profile & 2571 2851 ) 2572 2852 … … 2576 2856 long_name = "wind component in x direction", & 2577 2857 units = "m/s", & 2578 kind = " large-scale scalar forcing",&2858 kind = "geostrophic", & 2579 2859 input_id = 1, & 2580 2860 output_file = output_file, & 2581 grid = scalar_profile_grid,&2582 intermediate_grid = scalar_profile_intermediate&2861 grid = averaged_scalar_profile, & 2862 intermediate_grid = averaged_scalar_profile & 2583 2863 ) 2584 2864 output_var_table(45) % to_be_processed = ls_forcing_variables_required … … 2592 2872 input_id = 1, & 2593 2873 output_file = output_file, & 2594 grid = scalar_profile_grid,&2595 intermediate_grid = scalar_profile_intermediate&2874 grid = averaged_scalar_profile, & 2875 intermediate_grid = averaged_scalar_profile & 2596 2876 ) 2597 2877 output_var_table(46) % to_be_processed = ls_forcing_variables_required … … 2605 2885 input_id = 1, & 2606 2886 output_file = output_file, & 2607 grid = w_profile_grid,&2608 intermediate_grid = w_profile_intermediate&2887 grid = averaged_scalar_profile, & 2888 intermediate_grid = averaged_scalar_profile & 2609 2889 ) 2610 2890 output_var_table(47) % to_be_processed = ls_forcing_variables_required … … 2618 2898 input_id = 1, & 2619 2899 output_file = output_file, & 2620 grid = w_profile_grid,&2621 intermediate_grid = w_profile_intermediate&2900 grid = averaged_w_profile, & 2901 intermediate_grid = averaged_w_profile & 2622 2902 ) 2623 2903 output_var_table(48) % to_be_processed = ls_forcing_variables_required … … 2632 2912 input_id = 1, & 2633 2913 output_file = output_file, & 2634 grid = scalar_profile_grid,&2635 intermediate_grid = scalar_profile_intermediate&2914 grid = averaged_scalar_profile, & 2915 intermediate_grid = averaged_scalar_profile & 2636 2916 ) 2637 2917 output_var_table(49) % to_be_processed = ls_forcing_variables_required … … 2645 2925 input_id = 1, & 2646 2926 output_file = output_file, & 2647 grid = scalar_profile_grid,&2648 intermediate_grid = scalar_profile_intermediate&2927 grid = averaged_scalar_profile, & 2928 intermediate_grid = averaged_scalar_profile & 2649 2929 ) 2650 2930 output_var_table(50) % to_be_processed = ls_forcing_variables_required … … 2658 2938 input_id = 1, & 2659 2939 output_file = output_file, & 2660 grid = scalar_profile_grid,&2661 intermediate_grid = scalar_profile_intermediate&2940 grid = averaged_scalar_profile, & 2941 intermediate_grid = averaged_scalar_profile & 2662 2942 ) 2663 2943 output_var_table(51) % to_be_processed = ls_forcing_variables_required … … 2669 2949 units = "kg/kg/s", & 2670 2950 kind = "large-scale scalar forcing", & 2671 input_id = 1, &2672 output_file = output_file, & 2673 grid = scalar_profile_grid,&2674 intermediate_grid = scalar_profile_intermediate&2951 input_id = 3, & 2952 output_file = output_file, & 2953 grid = averaged_scalar_profile, & 2954 intermediate_grid = averaged_scalar_profile & 2675 2955 ) 2676 2956 output_var_table(52) % to_be_processed = ls_forcing_variables_required … … 2683 2963 units = "kg/kg/s", & 2684 2964 kind = "large-scale scalar forcing", & 2685 input_id = 1, &2686 output_file = output_file, & 2687 grid = scalar_profile_grid,&2688 intermediate_grid = scalar_profile_intermediate&2965 input_id = 3, & 2966 output_file = output_file, & 2967 grid = averaged_scalar_profile, & 2968 intermediate_grid = averaged_scalar_profile & 2689 2969 ) 2690 2970 output_var_table(53) % to_be_processed = ls_forcing_variables_required … … 2696 2976 units = "kg/kg", & 2697 2977 kind = "large-scale scalar forcing", & 2698 input_id = 1, &2699 output_file = output_file, & 2700 grid = scalar_profile_grid,&2701 intermediate_grid = scalar_profile_intermediate&2978 input_id = 3, & 2979 output_file = output_file, & 2980 grid = averaged_scalar_profile, & 2981 intermediate_grid = averaged_scalar_profile & 2702 2982 ) 2703 2983 output_var_table(54) % to_be_processed = ls_forcing_variables_required … … 2711 2991 input_id = 1, & 2712 2992 output_file = output_file, & 2713 grid = scalar_profile_grid,&2714 intermediate_grid = scalar_profile_intermediate&2993 grid = averaged_scalar_profile, & 2994 intermediate_grid = averaged_scalar_profile & 2715 2995 ) 2716 2996 output_var_table(55) % to_be_processed = ls_forcing_variables_required 2717 2997 2998 2999 output_var_table(56) = init_nc_var( & 3000 name = 'internal_density_centre', & 3001 std_name = "", & 3002 long_name = "", & 3003 units = "", & 3004 kind = "internal profile", & 3005 input_id = 4, & 3006 output_file = output_file, & 3007 grid = averaged_scalar_profile, & 3008 intermediate_grid = averaged_scalar_profile & 3009 ) 3010 output_var_table(56) % averaging_grid => averaged_scalar_profile 3011 3012 3013 output_var_table(57) = init_nc_var( & 3014 name = 'internal_density_north', & 3015 std_name = "", & 3016 long_name = "", & 3017 units = "", & 3018 kind = "internal profile", & 3019 input_id = 4, & 3020 output_file = output_file, & 3021 grid = north_averaged_scalar_profile, & 3022 intermediate_grid = north_averaged_scalar_profile & 3023 ) 3024 output_var_table(57) % averaging_grid => north_averaged_scalar_profile 3025 output_var_table(57) % to_be_processed = .NOT. cfg % ug_is_set 3026 3027 3028 output_var_table(58) = init_nc_var( & 3029 name = 'internal_density_south', & 3030 std_name = "", & 3031 long_name = "", & 3032 units = "", & 3033 kind = "internal profile", & 3034 input_id = 4, & 3035 output_file = output_file, & 3036 grid = south_averaged_scalar_profile, & 3037 intermediate_grid = south_averaged_scalar_profile & 3038 ) 3039 output_var_table(58) % averaging_grid => south_averaged_scalar_profile 3040 output_var_table(58) % to_be_processed = .NOT. cfg % ug_is_set 3041 3042 3043 output_var_table(59) = init_nc_var( & 3044 name = 'internal_density_east', & 3045 std_name = "", & 3046 long_name = "", & 3047 units = "", & 3048 kind = "internal profile", & 3049 input_id = 4, & 3050 output_file = output_file, & 3051 grid = east_averaged_scalar_profile, & 3052 intermediate_grid = east_averaged_scalar_profile & 3053 ) 3054 output_var_table(59) % averaging_grid => east_averaged_scalar_profile 3055 output_var_table(59) % to_be_processed = .NOT. cfg % ug_is_set 3056 3057 3058 output_var_table(60) = init_nc_var( & 3059 name = 'internal_density_west', & 3060 std_name = "", & 3061 long_name = "", & 3062 units = "", & 3063 kind = "internal profile", & 3064 input_id = 4, & 3065 output_file = output_file, & 3066 grid = west_averaged_scalar_profile, & 3067 intermediate_grid = west_averaged_scalar_profile & 3068 ) 3069 output_var_table(60) % averaging_grid => west_averaged_scalar_profile 3070 output_var_table(60) % to_be_processed = .NOT. cfg % ug_is_set 3071 3072 output_var_table(61) = init_nc_var( & 3073 name = 'internal_pressure_north', & 3074 std_name = "", & 3075 long_name = "", & 3076 units = "", & 3077 kind = "internal profile", & 3078 input_id = 2, & 3079 output_file = output_file, & 3080 grid = north_averaged_scalar_profile, & 3081 intermediate_grid = north_averaged_scalar_profile & 3082 ) 3083 output_var_table(61) % averaging_grid => north_averaged_scalar_profile 3084 output_var_table(61) % to_be_processed = .NOT. cfg % ug_is_set 3085 3086 3087 output_var_table(62) = init_nc_var( & 3088 name = 'internal_pressure_south', & 3089 std_name = "", & 3090 long_name = "", & 3091 units = "", & 3092 kind = "internal profile", & 3093 input_id = 2, & 3094 output_file = output_file, & 3095 grid = south_averaged_scalar_profile, & 3096 intermediate_grid = south_averaged_scalar_profile & 3097 ) 3098 output_var_table(62) % averaging_grid => south_averaged_scalar_profile 3099 output_var_table(62) % to_be_processed = .NOT. cfg % ug_is_set 3100 3101 3102 output_var_table(63) = init_nc_var( & 3103 name = 'internal_pressure_east', & 3104 std_name = "", & 3105 long_name = "", & 3106 units = "", & 3107 kind = "internal profile", & 3108 input_id = 2, & 3109 output_file = output_file, & 3110 grid = east_averaged_scalar_profile, & 3111 intermediate_grid = east_averaged_scalar_profile & 3112 ) 3113 output_var_table(63) % averaging_grid => east_averaged_scalar_profile 3114 output_var_table(63) % to_be_processed = .NOT. cfg % ug_is_set 3115 3116 3117 output_var_table(64) = init_nc_var( & 3118 name = 'internal_pressure_west', & 3119 std_name = "", & 3120 long_name = "", & 3121 units = "", & 3122 kind = "internal profile", & 3123 input_id = 2, & 3124 output_file = output_file, & 3125 grid = west_averaged_scalar_profile, & 3126 intermediate_grid = west_averaged_scalar_profile & 3127 ) 3128 output_var_table(64) % averaging_grid => west_averaged_scalar_profile 3129 output_var_table(64) % to_be_processed = .NOT. cfg % ug_is_set 2718 3130 2719 3131 ! Attributes shared among all variables … … 2732 3144 !------------------------------------------------------------------------------! 2733 3145 FUNCTION init_nc_var(name, std_name, long_name, units, kind, input_id, & 2734 grid, intermediate_grid, output_file, is_profile 2735 )RESULT(var)3146 grid, intermediate_grid, output_file, is_profile) & 3147 RESULT(var) 2736 3148 2737 3149 CHARACTER(LEN=*), INTENT(IN) :: name, std_name, long_name, units, kind … … 2772 3184 var % dimvarids(1:3) = output_file % dimvarids_soil 2773 3185 var % to_be_processed = init_variables_required 3186 var % is_internal = .FALSE. 2774 3187 var % task = "interpolate_2d" 2775 3188 … … 2781 3194 var % dimvarids(1:3) = output_file % dimvarids_scl 2782 3195 var % to_be_processed = init_variables_required 3196 var % is_internal = .FALSE. 2783 3197 var % task = "interpolate_3d" 2784 3198 … … 2794 3208 var % dimvarids(3) = output_file % dimvarids_scl(3) 2795 3209 var % to_be_processed = init_variables_required 3210 var % is_internal = .FALSE. 2796 3211 var % task = "interpolate_3d" 2797 3212 … … 2807 3222 var % dimvarids(3) = output_file % dimvarids_scl(3) 2808 3223 var % to_be_processed = init_variables_required 3224 var % is_internal = .FALSE. 2809 3225 var % task = "interpolate_3d" 2810 3226 … … 2820 3236 var % dimvarids(3) = output_file % dimvarids_vel(3) 2821 3237 var % to_be_processed = init_variables_required 3238 var % is_internal = .FALSE. 2822 3239 var % task = "interpolate_3d" 2823 3240 … … 2829 3246 var % dimvarids(1) = output_file % dimvarids_scl(3) !z 2830 3247 var % to_be_processed = init_variables_required 3248 var % is_internal = .FALSE. 2831 3249 var % task = "average profile" 2832 3250 … … 2838 3256 var % dimvarids(1) = output_file % dimvarids_vel(3) !z 2839 3257 var % to_be_processed = init_variables_required 3258 var % is_internal = .FALSE. 2840 3259 var % task = "average profile" 2841 3260 … … 2848 3267 var % dimvarids(1:2) = output_file % dimvarids_soil(1:2) 2849 3268 var % to_be_processed = boundary_variables_required 3269 var % is_internal = .FALSE. 2850 3270 var % task = "interpolate_2d" 2851 3271 … … 2860 3280 var % dimvarids(2) = output_file % dimvarids_scl(3) 2861 3281 var % to_be_processed = boundary_variables_required 3282 var % is_internal = .FALSE. 2862 3283 var % task = "interpolate_3d" 2863 3284 … … 2872 3293 var % dimvarids(2) = output_file % dimvarids_scl(3) 2873 3294 var % to_be_processed = boundary_variables_required 3295 var % is_internal = .FALSE. 2874 3296 var % task = "interpolate_3d" 2875 3297 … … 2884 3306 var % dimvarids(2) = output_file % dimvarids_scl(2) 2885 3307 var % to_be_processed = boundary_variables_required 3308 var % is_internal = .FALSE. 2886 3309 var % task = "interpolate_3d" 2887 3310 … … 2896 3319 var % dimvarids(2) = output_file % dimvarids_scl(3) 2897 3320 var % to_be_processed = boundary_variables_required 3321 var % is_internal = .FALSE. 2898 3322 var % task = "interpolate_3d" 2899 3323 … … 2908 3332 var % dimvarids(2) = output_file % dimvarids_scl(3) 2909 3333 var % to_be_processed = boundary_variables_required 3334 var % is_internal = .FALSE. 2910 3335 var % task = "interpolate_3d" 2911 3336 … … 2920 3345 var % dimvarids(2) = output_file % dimvarids_scl(2) 2921 3346 var % to_be_processed = boundary_variables_required 3347 var % is_internal = .FALSE. 2922 3348 var % task = "interpolate_3d" 2923 3349 … … 2932 3358 var % dimvarids(2) = output_file % dimvarids_scl(3) 2933 3359 var % to_be_processed = boundary_variables_required 3360 var % is_internal = .FALSE. 2934 3361 var % task = "interpolate_3d" 2935 3362 … … 2944 3371 var % dimvarids(2) = output_file % dimvarids_scl(3) 2945 3372 var % to_be_processed = boundary_variables_required 3373 var % is_internal = .FALSE. 2946 3374 var % task = "interpolate_3d" 2947 3375 … … 2956 3384 var % dimvarids(2) = output_file % dimvarids_vel(2) 2957 3385 var % to_be_processed = boundary_variables_required 3386 var % is_internal = .FALSE. 2958 3387 var % task = "interpolate_3d" 2959 3388 … … 2968 3397 var % dimvarids(2) = output_file % dimvarids_vel(3) 2969 3398 var % to_be_processed = boundary_variables_required 3399 var % is_internal = .FALSE. 2970 3400 var % task = "interpolate_3d" 2971 3401 … … 2980 3410 var % dimvarids(2) = output_file % dimvarids_vel(3) 2981 3411 var % to_be_processed = boundary_variables_required 3412 var % is_internal = .FALSE. 2982 3413 var % task = "interpolate_3d" 2983 3414 … … 2988 3419 var % dimvarids(1) = output_file % dimvarid_time 2989 3420 var % to_be_processed = .TRUE. 2990 var % task = "average scalar" 3421 var % is_internal = .FALSE. 3422 var % task = "average profile" 2991 3423 2992 3424 CASE( 'constant scalar profile' ) … … 2998 3430 var % dimvarids(1) = output_file % dimvarids_scl(3) 2999 3431 var % to_be_processed = .TRUE. 3432 var % is_internal = .FALSE. 3000 3433 var % task = "set profile" 3001 3434 … … 3008 3441 var % dimvarids(1) = output_file % dimvarids_scl(3) 3009 3442 var % to_be_processed = ls_forcing_variables_required 3443 var % is_internal = .FALSE. 3010 3444 var % task = "average large-scale profile" 3445 3446 CASE( 'geostrophic' ) 3447 var % lod = -1 3448 var % ndim = 2 3449 var % dimids(2) = output_file % dimid_time !t 3450 var % dimids(1) = output_file % dimids_scl(3) !z 3451 var % dimvarids(2) = output_file % dimvarid_time 3452 var % dimvarids(1) = output_file % dimvarids_scl(3) 3453 var % to_be_processed = .TRUE. 3454 var % is_internal = .FALSE. 3455 var % task = "geostrophic winds" 3011 3456 3012 3457 CASE( 'large-scale w forcing' ) … … 3018 3463 var % dimvarids(1) = output_file % dimvarids_vel(3) 3019 3464 var % to_be_processed = ls_forcing_variables_required 3465 var % is_internal = .FALSE. 3020 3466 var % task = "average large-scale profile" 3467 3468 CASE( 'internal profile' ) 3469 var % lod = -1 3470 var % ndim = 2 3471 var % dimids(2) = output_file % dimid_time !t 3472 var % dimids(1) = output_file % dimids_scl(3) !z 3473 var % dimvarids(2) = output_file % dimvarid_time 3474 var % dimvarids(1) = output_file % dimvarids_scl(3) 3475 var % to_be_processed = .TRUE. 3476 var % is_internal = .TRUE. 3477 var % task = "internal profile" 3021 3478 3022 3479 CASE DEFAULT … … 3031 3488 SUBROUTINE fini_variables() 3032 3489 3033 CALL report('fini_variables', 'Deallocating variable table' )3490 CALL report('fini_variables', 'Deallocating variable table', cfg % debug) 3034 3491 DEALLOCATE( input_var_table ) 3035 3492 … … 3039 3496 SUBROUTINE fini_io_groups() 3040 3497 3041 CALL report('fini_io_groups', 'Deallocating IO groups' )3498 CALL report('fini_io_groups', 'Deallocating IO groups', cfg % debug) 3042 3499 DEALLOCATE( io_group_list ) 3043 3500 … … 3047 3504 SUBROUTINE fini_file_lists() 3048 3505 3049 CALL report('fini_file_lists', 'Deallocating file lists' )3506 CALL report('fini_file_lists', 'Deallocating file lists', cfg % debug) 3050 3507 DEALLOCATE( flow_files, soil_files, radiation_files, soil_moisture_files ) 3051 3508 … … 3134 3591 ! rotated velocities 3135 3592 DO k = 1, nz 3136 DO j = 2, ny3137 DO i = 2, nx3593 DO j = 1, ny 3594 DO i = 1, nx 3138 3595 CALL uv2uvrot( urot = preprocess_buffer(1) % array(i,j,k), & 3139 3596 vrot = preprocess_buffer(2) % array(i,j,k), & … … 3171 3628 CALL report('preprocess', message) 3172 3629 3173 CASE( 't emperature' ) ! P and T3630 CASE( 'thermodynamics' ) ! T, P, QV 3174 3631 nx = SIZE(input_buffer(1) % array, 1) 3175 3632 ny = SIZE(input_buffer(1) % array, 2) … … 3192 3649 3193 3650 input_buffer (2) % array(i,j,:) = & 3194 input_buffer (2) % array(i,j,:) + basic_state_pressure(:) 3651 HECTO * input_buffer (2) % array(i,j,:) + & 3652 basic_state_pressure(:) 3195 3653 3196 3654 END DO … … 3204 3662 3205 3663 END IF 3664 ! pressure 3665 input_buffer(2) % is_preprocessed = .TRUE. 3666 3667 ! Copy temperature to last input buffer array 3668 ALLOCATE( & 3669 input_buffer( group % n_output_quantities ) % array (nx, ny, nz) & 3670 ) 3671 input_buffer(group % n_output_quantities) % array(:,:,:) = & 3672 input_buffer(1) % array(:,:,:) 3206 3673 3207 3674 ! Convert absolute to potential temperature 3208 input_buffer(1) % array(:,:,:) = input_buffer(1) % array(:,:,:) * & 3209 EXP( RD_PALM/CP_PALM * LOG( P_REF / input_buffer(2) % array(:,:,:) )) 3210 3211 input_buffer(1:2) % is_preprocessed = .TRUE. 3675 CALL potential_temperature( & 3676 t = input_buffer(1) % array(:,:,:), & 3677 p = input_buffer(2) % array(:,:,:), & 3678 p_ref = P_REF, & 3679 r = RD_PALM, & 3680 cp = CP_PALM & 3681 ) 3682 3683 ! potential temperature 3684 input_buffer(1) % is_preprocessed = .TRUE. 3685 3686 ! Convert temperature copy to density 3687 CALL moist_density( & 3688 t_rho = input_buffer(group % n_output_quantities) % array(:,:,:), & 3689 p = input_buffer(2) % array(:,:,:), & 3690 qv = input_buffer(3) % array(:,:,:), & 3691 rd = RD, & 3692 rv = RV & 3693 ) 3694 3695 ! qv 3696 input_buffer(3) % is_preprocessed = .TRUE. 3697 3698 ! density 3699 input_buffer(group % n_output_quantities) % is_preprocessed = .TRUE. 3700 3701 3212 3702 message = "Input buffers for group '" // TRIM(group % kind) // "'"//& 3213 3703 " preprocessed sucessfully." … … 3336 3826 3337 3827 CASE DEFAULT 3338 message = " Bufferkind '" // TRIM(group % kind) // "' is not supported."3828 message = "IO group kind '" // TRIM(group % kind) // "' is not supported." 3339 3829 CALL abort('prerpocess', message) 3340 3830
Note: See TracChangeset
for help on using the changeset viewer.