Changeset 3557 for palm/trunk/UTIL/inifor/src
- Timestamp:
- Nov 22, 2018 4:01:22 PM (6 years ago)
- Location:
- palm/trunk/UTIL/inifor/src
- Files:
-
- 8 edited
-
inifor.f90 (modified) (14 diffs)
-
inifor_control.f90 (modified) (8 diffs)
-
inifor_defs.f90 (modified) (4 diffs)
-
inifor_grid.f90 (modified) (89 diffs)
-
inifor_io.f90 (modified) (40 diffs)
-
inifor_transform.f90 (modified) (47 diffs)
-
inifor_types.f90 (modified) (7 diffs)
-
inifor_util.f90 (modified) (17 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/UTIL/inifor/src/inifor.f90
r3537 r3557 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Updated documentation 29 ! 30 ! 3537 2018-11-20 10:53:14Z eckhard 28 31 ! Print version number on program start 29 32 ! … … 54 57 ! Authors: 55 58 ! -------- 56 ! @author Eckhard Kadasch59 !> @author Eckhard Kadasch, (Deutscher Wetterdienst, Offenbach) 57 60 ! 58 61 ! Description: … … 84 87 IMPLICIT NONE 85 88 86 INTEGER :: igroup 87 INTEGER :: ivar 88 INTEGER :: iter 89 REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: output_arr 90 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_centre 91 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: ug_arr 92 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: vg_arr 93 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_north 94 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_south 95 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_east 96 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_west 97 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_north 98 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_south 99 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_east 100 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_west 101 REAL(dp), POINTER, DIMENSION(:) :: internal_arr 102 REAL(dp), POINTER, DIMENSION(:) :: ug_vg_arr 103 TYPE(nc_var), POINTER :: output_var 104 TYPE(io_group), POINTER :: group 105 TYPE(container), ALLOCATABLE :: input_buffer(:) 106 LOGICAL, SAVE :: ug_vg_have_been_computed = .FALSE. 107 LOGICAL, SAVE :: debugging_output = .TRUE. 89 INTEGER :: igroup !< loop index for IO groups 90 INTEGER :: ivar !< loop index for output variables 91 INTEGER :: iter !< loop index for time steps 92 93 REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: output_arr !< array buffer for interpolated quantities 94 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_centre !< density profile of the centre averaging domain 95 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: ug_arr !< geostrophic wind in x direction 96 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: vg_arr !< geostrophic wind in y direction 97 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_north !< density profile of the northern averaging domain 98 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_south !< density profile of the southern averaging domain 99 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_east !< density profile of the eastern averaging domain 100 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_west !< density profile of the western averaging domain 101 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_north !< pressure profile of the northern averaging domain 102 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_south !< pressure profile of the southern averaging domain 103 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_east !< pressure profile of the eastern averaging domain 104 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_west !< pressure profile of the western averaging domain 105 106 REAL(dp), POINTER, DIMENSION(:) :: internal_arr !< pointer to the currently processed internal array (density, pressure) 107 REAL(dp), POINTER, DIMENSION(:) :: ug_vg_arr !< pointer to the currently processed geostrophic wind component 108 109 TYPE(nc_var), POINTER :: output_var !< pointer to the currently processed output variable 110 TYPE(io_group), POINTER :: group !< pointer to the currently processed IO group 111 TYPE(container), ALLOCATABLE :: input_buffer(:) !< buffer of the current IO group's input arrays 112 113 LOGICAL, SAVE :: ug_vg_have_been_computed = .FALSE. !< flag for managing geostrophic wind allocation and computation 114 LOGICAL, SAVE :: debugging_output = .TRUE. !< flag controllging output of internal variables 108 115 109 116 !> \mainpage About INIFOR … … 116 123 CALL report('main_loop', 'Running INIFOR version ' // VERSION) 117 124 118 ! Initialize INIFOR's parameters from command-line interface and namelists 125 ! 126 !-- Initialize INIFOR's parameters from command-line interface and namelists 119 127 CALL setup_parameters() 120 128 121 ! Initialize all grids, including interpolation neighbours and weights 129 ! 130 !-- Initialize all grids, including interpolation neighbours and weights 122 131 CALL setup_grids() 123 132 CALL run_control('time', 'init') 124 133 125 ! Initialize the netCDF output file and define dimensions 134 ! 135 !-- Initialize the netCDF output file and define dimensions 126 136 CALL setup_netcdf_dimensions(output_file, palm_grid, cfg % start_date, & 127 137 origin_lon, origin_lat) 128 138 CALL run_control('time', 'write') 129 139 130 ! Set up the tables containing the input and output variables and set 131 ! the corresponding netCDF dimensions for each output variable 140 ! 141 !-- Set up the tables containing the input and output variables and set 142 !-- the corresponding netCDF dimensions for each output variable 132 143 CALL setup_variable_tables(cfg % ic_mode) 133 144 CALL run_control('time', 'write') 134 145 135 ! Add the output variables to the netCDF output file 146 ! 147 !-- Add the output variables to the netCDF output file 136 148 CALL setup_netcdf_variables(output_file % name, output_var_table, & 137 149 cfg % debug) … … 143 155 !- Section 2: Main loop 144 156 !------------------------------------------------------------------------------ 145 146 157 DO igroup = 1, SIZE(io_group_list) 147 158 … … 345 356 CALL run_control('time', 'comp') 346 357 347 348 !This case gets called twice, the first time for ug, the349 !second time for vg. We compute ug and vg at the first call350 !and keep vg (and ug for that matter) around for the second351 !call.358 ! 359 !-- This case gets called twice, the first time for ug, the 360 !-- second time for vg. We compute ug and vg at the first call 361 !-- and keep vg (and ug for that matter) around for the second 362 !-- call. 352 363 CASE ( 'geostrophic winds' ) 353 364 … … 360 371 vg_arr = cfg % vg 361 372 ELSE 362 CALL geostrophic_winds( p_north, p_south, p_east, & 363 p_west, rho_centre, f3, & 364 averaging_width_ew, & 365 averaging_width_ns, & 366 phi_n, lambda_n, phi_centre, lam_centre, & 373 CALL geostrophic_winds( p_north, p_south, p_east, & 374 p_west, rho_centre, f3, & 375 averaging_width_ew, & 376 averaging_width_ns, & 377 phi_n, lambda_n, & 378 phi_centre, lam_centre, & 367 379 ug_arr, vg_arr ) 368 380 END IF … … 372 384 END IF 373 385 374 ! Prepare output of geostrophic winds 386 ! 387 !-- Prepare output of geostrophic winds 375 388 SELECT CASE(TRIM(output_var % name)) 376 389 CASE ('ls_forcing_ug') … … 397 410 SELECT CASE (TRIM(output_var % name)) 398 411 399 !CASE('ls_forcing_ug')400 ! output_arr(1, 1, :) = ug401 402 !CASE('ls_forcing_vg')403 ! output_arr(1, 1, :) = vg404 405 412 CASE('nudging_tau') 406 413 output_arr(1, 1, :) = NUDGING_TAU … … 418 425 "has not been implemented, yet." 419 426 CALL abort('main loop', message) 420 !ALLOCATE( output_arr( 1, 1, 1:nz ) )421 427 422 428 CASE DEFAULT … … 444 450 END IF 445 451 446 END DO ! ouput variables 452 ! 453 !-- output variable loop 454 END DO 447 455 448 456 ug_vg_have_been_computed = .FALSE. … … 463 471 END IF 464 472 473 ! 474 !-- Keep input buffer around for averaged (radiation) and 475 !-- accumulated COSMO quantities (precipitation). 465 476 IF ( group % kind == 'running average' .OR. & 466 477 group % kind == 'accumulated' ) THEN 467 ! Keep input buffer around for averaged (radiation) and468 ! accumulated COSMO-DE quantities (precipitation).469 478 ELSE 470 479 CALL report('main loop', 'Deallocating input buffer', cfg % debug) … … 473 482 CALL run_control('time', 'alloc') 474 483 475 END DO ! time steps / input files 484 ! 485 !-- time steps / input files loop 486 END DO 476 487 477 488 IF (ALLOCATED(input_buffer)) THEN … … 491 502 CALL report('main loop', message, cfg % debug) 492 503 493 END IF ! IO group % to_be_processed 494 495 END DO ! IO groups 504 ! 505 !-- IO group % to_be_processed conditional 506 END IF 507 508 ! 509 !-- IO groups loop 510 END DO 496 511 497 512 !------------------------------------------------------------------------------ -
palm/trunk/UTIL/inifor/src/inifor_control.f90
r3447 r3557 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Updated documentation 29 ! 30 ! 31 ! 3447 2018-10-29 15:52:54Z eckhard 28 32 ! Renamed source files for compatibilty with PALM build system 29 33 ! … … 44 48 ! Authors: 45 49 ! -------- 46 ! @author Eckhard Kadasch50 !> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach) 47 51 ! 48 52 ! Description: … … 61 65 IMPLICIT NONE 62 66 63 CHARACTER (LEN=5000) :: message = '' 67 CHARACTER (LEN=5000) :: message = '' !< log message buffer 64 68 65 69 CONTAINS 66 70 71 !------------------------------------------------------------------------------! 72 ! Description: 73 ! ------------ 74 !> 75 !> report() is INIFOR's general logging routine. It prints the given 'message' 76 !> to the terminal and writes it to the INIFOR log file. 77 !> 78 !> You can use this routine to log events across INIFOR's code to log. For 79 !> warnings and abort messages, you may use the dedicated routines warn() and 80 !> abort() in this module. Both use report() and add specific behaviour to it. 81 !------------------------------------------------------------------------------! 67 82 SUBROUTINE report(routine, message, debug) 68 83 69 CHARACTER(LEN=*), INTENT(IN) :: routine 70 CHARACTER(LEN=*), INTENT(IN) :: message 71 LOGICAL, OPTIONAL, INTENT(IN) :: debug 72 INTEGER :: u 73 LOGICAL, SAVE :: is_first_run = .TRUE. 74 LOGICAL :: suppress_message 75 76 77 IF (is_first_run) THEN 84 CHARACTER(LEN=*), INTENT(IN) :: routine !< name of calling subroutine of function 85 CHARACTER(LEN=*), INTENT(IN) :: message !< log message 86 LOGICAL, OPTIONAL, INTENT(IN) :: debug !< flag the current message as debugging message 87 88 INTEGER :: u !< Fortran file unit for the log file 89 LOGICAL, SAVE :: is_first_run = .TRUE. !< control flag for file opening mode 90 LOGICAL :: suppress_message !< control falg for additional debugging log 91 92 93 IF ( is_first_run ) THEN 78 94 OPEN( NEWUNIT=u, FILE='inifor.log', STATUS='replace' ) 79 95 is_first_run = .FALSE. … … 84 100 85 101 suppress_message = .FALSE. 86 IF ( PRESENT(debug)) THEN87 IF ( .NOT. debug) suppress_message = .TRUE.102 IF ( PRESENT(debug) ) THEN 103 IF ( .NOT. debug ) suppress_message = .TRUE. 88 104 END IF 89 105 90 IF ( .NOT. suppress_message) THEN106 IF ( .NOT. suppress_message ) THEN 91 107 PRINT *, "inifor: " // TRIM(message) // " [ " // TRIM(routine) // " ]" 92 108 WRITE(u, *) TRIM(message) // " [ " // TRIM(routine) // " ]" … … 98 114 99 115 116 !------------------------------------------------------------------------------! 117 ! Description: 118 ! ------------ 119 !> 120 !> warn() prepends "WARNING:" the given 'message' and prints the result to the 121 !> terminal and writes it to the INIFOR logfile. 122 !> 123 !> You can use this routine for messaging issues, that still allow INIFOR to 124 !> continue. 125 !------------------------------------------------------------------------------! 100 126 SUBROUTINE warn(routine, message) 101 127 102 CHARACTER(LEN=*), INTENT(IN) :: routine 103 CHARACTER(LEN=*), INTENT(IN) :: message 128 CHARACTER(LEN=*), INTENT(IN) :: routine !< name of calling subroutine or function 129 CHARACTER(LEN=*), INTENT(IN) :: message !< log message 104 130 105 131 CALL report(routine, "WARNING: " // TRIM(message)) … … 108 134 109 135 136 !------------------------------------------------------------------------------! 137 ! Description: 138 ! ------------ 139 !> 140 !> abort() prepends "ERROR:" the given 'message' and prints the result to the 141 !> terminal, writes it to the INIFOR logfile, and exits INIFOR. 142 !> 143 !> You can use this routine for messaging issues, that are critical and prevent 144 !> INIFOR from continueing. 145 !------------------------------------------------------------------------------! 110 146 SUBROUTINE abort(routine, message) 111 147 112 CHARACTER(LEN=*), INTENT(IN) :: routine 113 CHARACTER(LEN=*), INTENT(IN) :: message 148 CHARACTER(LEN=*), INTENT(IN) :: routine !< name of calling subroutine or function 149 CHARACTER(LEN=*), INTENT(IN) :: message !< log message 114 150 115 151 CALL report(routine, "ERROR: " // TRIM(message) // " Stopping.") … … 119 155 120 156 157 !------------------------------------------------------------------------------! 158 ! Description: 159 ! ------------ 160 !> 161 !> print_version() prints the INIFOR version number and copyright notice. 162 !------------------------------------------------------------------------------! 121 163 SUBROUTINE print_version() 122 164 PRINT *, "INIFOR " // VERSION … … 125 167 126 168 169 !------------------------------------------------------------------------------! 170 ! Description: 171 ! ------------ 172 !> 173 !> run_control() measures the run times of various parts of INIFOR and 174 !> accumulates them in timing budgets. 175 !------------------------------------------------------------------------------! 127 176 SUBROUTINE run_control(mode, budget) 128 177 129 CHARACTER(LEN=*), INTENT(IN) :: mode, budget 130 REAL(dp), SAVE :: t0, t1 131 REAL(dp), SAVE :: t_comp=0.0_dp, & 132 t_alloc=0.0_dp, & 133 t_init=0.0_dp, & 134 t_read=0.0_dp, & 135 t_total=0.0_dp, & 136 t_write=0.0_dp 137 CHARACTER(LEN=*), PARAMETER :: fmt='(F6.2)' 178 CHARACTER(LEN=*), INTENT(IN) :: mode !< name of the calling mode 179 CHARACTER(LEN=*), INTENT(IN) :: budget !< name of the timing budget 180 181 REAL(dp), SAVE :: t0 !< begin of timing interval 182 REAL(dp), SAVE :: t1 !< end of timing interval 183 REAL(dp), SAVE :: t_comp = 0.0_dp !< computation timing budget 184 REAL(dp), SAVE :: t_alloc = 0.0_dp !< allocation timing budget 185 REAL(dp), SAVE :: t_init = 0.0_dp !< initialization timing budget 186 REAL(dp), SAVE :: t_read = 0.0_dp !< reading timing budget 187 REAL(dp), SAVE :: t_total = 0.0_dp !< total time 188 REAL(dp), SAVE :: t_write = 0.0_dp !< writing timing budget 189 190 CHARACTER(LEN=*), PARAMETER :: fmt='(F6.2)' !< floating-point output format 138 191 139 192 -
palm/trunk/UTIL/inifor/src/inifor_defs.f90
r3537 r3557 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Updated documentation 29 ! 30 ! 31 ! 3537 2018-11-20 10:53:14Z eckhard 28 32 ! Bumped version number 29 33 ! … … 59 63 ! Authors: 60 64 ! -------- 61 ! @author Eckhard Kadasch65 !> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach) 62 66 ! 63 67 ! Description: … … 68 72 69 73 IMPLICIT NONE 70 71 ! Parameters for type definitions 74 75 ! 76 !-- Parameters for type definitions 72 77 INTEGER, PARAMETER :: dp = 8 !< double precision (8 bytes = 64 bits) 73 78 INTEGER, PARAMETER :: sp = 4 !< single precision (4 bytes = 32 bits) … … 78 83 INTEGER, PARAMETER :: DATE = 10 !< length of date strings 79 84 80 ! Trigonomentry 85 ! 86 !-- Trigonomentry 81 87 REAL(dp), PARAMETER :: PI = 3.14159265358979323846264338_dp !< Ratio of a circle's circumference to its diamter [-] 82 88 REAL(dp), PARAMETER :: TO_RADIANS = PI / 180.0_dp !< Conversion factor from degrees to radiant [-] 83 89 REAL(dp), PARAMETER :: TO_DEGREES = 180.0_dp / PI !< Conversion factor from radians to degrees [-] 84 90 85 ! COSMO-DE parameters 86 INTEGER, PARAMETER :: WATER_ID = 9 !< Integer corresponding to the water soil type in COSMO-DE [-] 87 REAL(dp), PARAMETER :: EARTH_RADIUS = 6371229.0_dp !< Earth radius used in COSMO-DE [m] 88 REAL(dp), PARAMETER :: P_SL = 1e5_dp !< Reference pressure for computation of COSMO-DE's basic state pressure [Pa] 89 REAL(dp), PARAMETER :: T_SL = 288.15_dp !< Reference temperature for computation of COSMO-DE's basic state pressure [K] 90 REAL(dp), PARAMETER :: BETA = 42.0_dp !< logarithmic lapse rate, dT / d ln(p), for computation of COSMO-DE's basic state pressure [K] 91 REAL(dp), PARAMETER :: RD = 287.05_dp !< specific gas constant of dry air, used in computation of COSMO-DE's basic state [J/kg/K] 92 REAL(dp), PARAMETER :: RV = 461.51_dp !< specific gas constant of water vapor [J/kg/K] 93 REAL(dp), PARAMETER :: G = 9.80665_dp !< acceleration of Earth's gravity, used in computation of COSMO-DE's basic state [m/s/s] 94 REAL(dp), PARAMETER :: RHO_L = 1e3_dp !< density of liquid water, used to convert W_SO from [kg/m^2] to [m^3/m^3], in [kg/m^3] 95 REAL(dp), PARAMETER :: HECTO = 100_dp !< unit conversion factor from hPa to Pa 91 ! 92 !-- COSMO-DE parameters 93 INTEGER, PARAMETER :: WATER_ID = 9 !< Integer corresponding to the water soil type in COSMO-DE [-] 94 REAL(dp), PARAMETER :: EARTH_RADIUS = 6371229.0_dp !< Earth radius used in COSMO-DE [m] 95 REAL(dp), PARAMETER :: P_SL = 1e5_dp !< Reference pressure for computation of COSMO-DE's basic state pressure [Pa] 96 REAL(dp), PARAMETER :: T_SL = 288.15_dp !< Reference temperature for computation of COSMO-DE's basic state pressure [K] 97 REAL(dp), PARAMETER :: BETA = 42.0_dp !< logarithmic lapse rate, dT / d ln(p), for computation of COSMO-DE's basic 98 !< state pressure [K] 99 REAL(dp), PARAMETER :: RD = 287.05_dp !< specific gas constant of dry air, used in computation of COSMO-DE's basic 100 !< state [J/kg/K] 101 REAL(dp), PARAMETER :: RV = 461.51_dp !< specific gas constant of water vapor [J/kg/K] 102 REAL(dp), PARAMETER :: G = 9.80665_dp !< acceleration of Earth's gravity, used in computation of COSMO-DE's basic 103 !< state [m/s/s] 104 REAL(dp), PARAMETER :: RHO_L = 1e3_dp !< density of liquid water, used to convert W_SO from [kg/m^2] to [m^3/m^3], 105 !< in [kg/m^3] 106 REAL(dp), PARAMETER :: HECTO = 100_dp !< unit conversion factor from hPa to Pa 96 107 97 ! PALM-4U parameters 98 REAL(dp), PARAMETER :: OMEGA = 7.29e-5_dp !< angular velocity of Earth's rotation [s^-1] 99 REAL(dp), PARAMETER :: P_REF = 1e5_dp !< Reference pressure for potential temperature [Pa] 100 REAL(dp), PARAMETER :: RD_PALM = 287.0_dp !< specific gas constant of dry air, used in computation of PALM-4U's potential temperature [J/kg/K] 101 REAL(dp), PARAMETER :: CP_PALM = 1005.0_dp !< heat capacity of dry air at constant pressure, used in computation of PALM-4U's potential temperature [J/kg/K] 108 ! 109 !-- PALM-4U parameters 110 REAL(dp), PARAMETER :: OMEGA = 7.29e-5_dp !< angular velocity of Earth's rotation [s^-1] 111 REAL(dp), PARAMETER :: P_REF = 1e5_dp !< Reference pressure for potential temperature [Pa] 112 REAL(dp), PARAMETER :: RD_PALM = 287.0_dp !< specific gas constant of dry air, used in computation of PALM-4U's potential temperature [J/kg/K] 113 REAL(dp), PARAMETER :: CP_PALM = 1005.0_dp !< heat capacity of dry air at constant pressure, used in computation of PALM-4U's potential temperature [J/kg/K] 102 114 103 ! INIFOR parameters 104 INTEGER, PARAMETER :: FILL_ITERATIONS = 5 !< Number of iterations for extrapolating soil data into COSMO-DE water cells [-] 115 ! 116 !-- INIFOR parameters 117 INTEGER, PARAMETER :: FILL_ITERATIONS = 5 !< Number of iterations for extrapolating soil data into COSMO-DE 118 !< water cells [-] 105 119 INTEGER, PARAMETER :: FORCING_STEP = 1 !< Number of hours between forcing time steps [h] 106 120 REAL(dp), PARAMETER :: NUDGING_TAU = 21600.0_dp !< Nudging relaxation time scale [s] -
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 -
palm/trunk/UTIL/inifor/src/inifor_io.f90
r3537 r3557 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Updated documentation, removed unused subroutine write_netcdf_variable_2d() 29 ! 30 ! 31 ! 3537 2018-11-20 10:53:14Z eckhard 28 32 ! New routine get_netcdf_dim_vector() 29 33 ! … … 73 77 ! Authors: 74 78 ! -------- 75 ! @author Eckhard Kadasch79 !> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach) 76 80 ! 77 81 ! Description: … … 92 96 IMPLICIT NONE 93 97 98 !------------------------------------------------------------------------------! 99 ! Description: 100 ! ------------ 101 !> get_netcdf_variable() reads the netCDF data and metadate for the netCDF 102 !> variable 'in_var % name' from the file 'in_file'. The netCDF data array is 103 !> stored in the 'buffer' array and metadata added to the respective members of 104 !> the given 'in_var'. 105 !------------------------------------------------------------------------------! 94 106 INTERFACE get_netcdf_variable 95 107 MODULE PROCEDURE get_netcdf_variable_int … … 101 113 CONTAINS 102 114 115 !------------------------------------------------------------------------------! 116 ! Description: 117 ! ------------ 118 !> get_netcdf_variable_int() implements the integer variant for the 119 !> get_netcdf_variable interface. 120 !------------------------------------------------------------------------------! 103 121 SUBROUTINE get_netcdf_variable_int(in_file, in_var, buffer) 104 122 … … 139 157 140 158 159 !------------------------------------------------------------------------------! 160 ! Description: 161 ! ------------ 162 !> get_netcdf_variable_real() implements the real variant for the 163 !> get_netcdf_variable interface. 164 !------------------------------------------------------------------------------! 141 165 SUBROUTINE get_netcdf_variable_real(in_file, in_var, buffer) 142 166 … … 177 201 178 202 179 SUBROUTINE get_netcdf_dim_vector(filename, varname, array) 203 !------------------------------------------------------------------------------! 204 ! Description: 205 ! ------------ 206 !> get_netcdf_dim_vector() reads the coordinate array 'coordname' from the 207 !> netCDF file 'filename'. 208 !------------------------------------------------------------------------------! 209 SUBROUTINE get_netcdf_dim_vector(filename, coordname, coords) 180 210 181 211 CHARACTER(LEN=*), INTENT(IN) :: filename 182 CHARACTER(LEN=*), INTENT(IN) :: varname183 REAL(dp), ALLOCATABLE, INTENT(INOUT) :: array(:)212 CHARACTER(LEN=*), INTENT(IN) :: coordname 213 REAL(dp), ALLOCATABLE, INTENT(INOUT) :: coords(:) 184 214 185 215 INTEGER :: ncid, varid, dimlen … … 187 217 188 218 IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. & 189 nf90_inq_varid( ncid, varname, varid ) .EQ. NF90_NOERR ) THEN219 nf90_inq_varid( ncid, coordname, varid ) .EQ. NF90_NOERR ) THEN 190 220 191 221 CALL check(nf90_inquire_variable( ncid, varid, dimids = dimids )) 192 222 CALL check(nf90_inquire_dimension( ncid, dimids(1), len = dimlen )) 193 223 194 ALLOCATE( array(dimlen))195 CALL check(nf90_get_var( ncid, varid, array))224 ALLOCATE(coords(dimlen)) 225 CALL check(nf90_get_var( ncid, varid, coords)) 196 226 197 227 ELSE 198 228 199 message = "Failed to read '" // TRIM( varname) // &229 message = "Failed to read '" // TRIM(coordname) // & 200 230 "' from file '" // TRIM(filename) // "'." 201 231 CALL abort('get_netcdf_dim_vector', message) … … 206 236 207 237 238 !------------------------------------------------------------------------------! 239 ! Description: 240 ! ------------ 241 !> get_input_dimensions() reads dimensions metadata of the netCDF variable given 242 !> by 'in_var % name' into 'in_var' data structure. 243 !------------------------------------------------------------------------------! 208 244 SUBROUTINE get_input_dimensions(in_var, ncid) 209 245 … … 232 268 233 269 270 !------------------------------------------------------------------------------! 271 ! Description: 272 ! ------------ 273 !> get_netcdf_start_and_count() gets the start position and element counts for 274 !> the given netCDF file. This information is used in get_netcdf_variable_int() 275 !> and _real() for reading input variables.. 276 !------------------------------------------------------------------------------! 234 277 SUBROUTINE get_netcdf_start_and_count(in_var, start, count) 235 278 … … 251 294 start = (/ 1, 1, 1 /) 252 295 IF ( TRIM(in_var % name) .EQ. 'T_SO' ) THEN 253 ! Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8 296 ! 297 !-- Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8 254 298 in_var % dimlen(3) = in_var % dimlen(3) - 1 255 299 256 ! Start reading from second level, e.g. depth = 0.005 instead of 0.0 300 ! 301 !-- Start reading from second level, e.g. depth = 0.005 instead of 0.0 257 302 start(3) = 2 258 303 END IF … … 269 314 270 315 316 !------------------------------------------------------------------------------! 317 ! Description: 318 ! ------------ 319 !> Routine for defining netCDF variables in the dynamic driver, INIFOR's netCDF 320 !> output. 321 !------------------------------------------------------------------------------! 271 322 SUBROUTINE netcdf_define_variable(var, ncid) 272 323 … … 286 337 287 338 339 !------------------------------------------------------------------------------! 340 ! Description: 341 ! ------------ 342 !> netcdf_get_dimensions() reads in all dimensions and their lengths and stores 343 !> them in the given the 'var' data structure. This information is used later 344 !> for writing output variables in update_output(). 345 !------------------------------------------------------------------------------! 288 346 SUBROUTINE netcdf_get_dimensions(var, ncid) 289 347 … … 293 351 CHARACTER(SNAME) :: null 294 352 295 ! Remember dimension lenghts for NetCDF writing routine296 353 DO i = 1, var % ndim 297 354 CALL check(nf90_inquire_dimension(ncid, var % dimids(i), & … … 306 363 ! Description: 307 364 ! ------------ 308 !> This routine initializes Inifor. This includes parsing command-line 309 !> arguments, setting the names of the input and output file names as well as 310 !> the name of the input namelist and, subsequently, reading in and setting grid 311 !> parameters for the PALM-4U computational grid. 365 !> This routine parses and interpretes the command-line options and stores the 366 !> resulting settings in the 'cfg' data structure. 312 367 !------------------------------------------------------------------------------! 313 368 SUBROUTINE parse_command_line_arguments( cfg ) … … 330 385 IF (arg_count .GT. 0) THEN 331 386 332 ! Every option should have an argument.333 !IF ( MOD(arg_count, 2) .NE. 0 ) THEN334 ! message = "Syntax error in command line."335 ! CALL abort('parse_command_line_arguments', message)336 !END IF337 338 387 message = "The -clon and -clat command line options are depricated. " // & 339 388 "Please remove them form your inifor command and specify the " // & … … 360 409 cfg % debug = .TRUE. 361 410 362 ! Elevation of the PALM-4U domain above sea level363 411 CASE( '-z0', '-z', '--elevation' ) 364 412 CALL get_option_argument( i, arg ) 365 413 READ(arg, *) cfg % z0 366 414 367 ! surface pressure, at z0368 415 CASE( '-p0', '-r', '--surface-pressure' ) 369 416 cfg % p0_is_set = .TRUE. … … 371 418 READ(arg, *) cfg % p0 372 419 373 ! geostrophic wind in x direction374 420 CASE( '-ug', '-u', '--geostrophic-u' ) 375 421 cfg % ug_is_set = .TRUE. … … 377 423 READ(arg, *) cfg % ug 378 424 379 ! geostrophic wind in y direction380 425 CASE( '-vg', '-v', '--geostrophic-v' ) 381 426 cfg % vg_is_set = .TRUE. … … 383 428 READ(arg, *) cfg % vg 384 429 385 ! domain centre geographical longitude and latitude386 430 CASE( '-clon', '-clat' ) 387 431 CALL abort('parse_command_line_arguments', message) 388 !READ(arg, *) lambda_cg389 !lambda_cg = lambda_cg * TO_RADIANS390 !READ(arg, *) phi_cg391 !phi_cg = phi_cg * TO_RADIANS392 432 393 433 CASE( '-path', '-p', '--path' ) … … 444 484 cfg % namelist_file = TRIM(arg) 445 485 446 ! initial condition mode: 'profile' / 'volume'447 486 CASE( '-mode', '-i', '--init-mode' ) 448 487 CALL get_option_argument( i, arg ) 449 488 cfg % ic_mode = TRIM(arg) 450 489 451 ! boundary conditions / forcing mode: 'ideal' / 'real'452 490 CASE( '-f', '--forcing-mode' ) 453 491 CALL get_option_argument( i, arg ) … … 484 522 485 523 524 !------------------------------------------------------------------------------! 525 ! Description: 526 ! ------------ 527 !> Get the argument of the i'th command line option, which is at the location 528 !> i+1 of the argument list. 529 !------------------------------------------------------------------------------! 486 530 SUBROUTINE get_option_argument(i, arg) 487 531 CHARACTER(LEN=PATH), INTENT(INOUT) :: arg … … 494 538 495 539 540 !------------------------------------------------------------------------------! 541 ! Description: 542 ! ------------ 543 !> Checks the INIFOR configuration 'cfg' for plausibility. 544 !------------------------------------------------------------------------------! 496 545 SUBROUTINE validate_config(cfg) 497 546 TYPE(inifor_config), INTENT(IN) :: cfg … … 503 552 all_files_present = all_files_present .AND. file_present(cfg % soiltyp_file, 'SOILTYP') 504 553 505 ! Only check optional static driver file name, if it has been given. 554 ! 555 !-- Only check optional static driver file name, if it has been given. 506 556 IF (TRIM(cfg % static_driver_file) .NE. '') THEN 507 557 all_files_present = all_files_present .AND. file_present(cfg % static_driver_file, 'static driver') … … 557 607 558 608 609 !------------------------------------------------------------------------------! 610 ! Description: 611 ! ------------ 612 !> Check whether the given file is present on the filesystem. 613 !------------------------------------------------------------------------------! 559 614 LOGICAL FUNCTION file_present(filename, kind) 560 615 CHARACTER(LEN=PATH), INTENT(IN) :: filename … … 585 640 ! Description: 586 641 ! ------------ 587 !> This routine initializes the Inifor output file, i.e. the PALM-4U588 !> initializing and forcing data as a NetCDFfile.642 !> This routine initializes the dynamic driver file, i.e. INIFOR's netCDF output 643 !> file. 589 644 !> 590 645 !> Besides writing metadata, such as global attributes, coordinates, variables, 591 !> in the NetCDF file, various NetCDF IDs are saved for later, when Inifor646 !> in the netCDF file, various netCDF IDs are saved for later, when INIFOR 592 647 !> writes the actual data. 593 648 !------------------------------------------------------------------------------! … … 611 666 CALL report('setup_netcdf_dimensions', message) 612 667 613 ! Create the NetCDF file. NF90_CLOBBER selects overwrite mode. 668 ! 669 !-- Create the netCDF file as in netCDF-4/HDF5 format if __netcdf4 preprocessor flag is given 614 670 #if defined( __netcdf4 ) 615 671 CALL check(nf90_create(TRIM(output_file % name), OR(NF90_CLOBBER, NF90_HDF5), ncid)) … … 646 702 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat', TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)')))) 647 703 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon', TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)')))) 648 ! FIXME: This is the elevation relative to COSMO-DE/D2 sea level and does 649 ! FIXME: not necessarily comply with DHHN2016 (c.f. PALM Input Data 650 ! FIXME: Standard v1.9., origin_z) 704 ! 705 !-- FIXME: This is the elevation relative to COSMO-DE/D2 sea level and does 706 !-- FIXME: not necessarily comply with DHHN2016 (c.f. PALM Input Data 707 !-- FIXME: Standard v1.9., origin_z) 651 708 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_z', TRIM(real_to_str(z0, '(F18.13)')))) 652 709 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION))) … … 658 715 !- Section 2a: Define dimensions for cell centers (scalars in soil and atmosph.) 659 716 !------------------------------------------------------------------------------ 660 dimids = (/0, 0, 0/) ! reset dimids 661 CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) ) 662 CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) ) 663 CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) ) 664 output_file % dimids_scl = dimids ! save dimids for later 665 666 dimvarids = (/0, 0, 0/) ! reset dimvarids 667 CALL check(nf90_def_var(ncid, "x", NF90_FLOAT, dimids(1), dimvarids(1))) 668 CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell centers")) 669 CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m")) 670 671 CALL check(nf90_def_var(ncid, "y", NF90_FLOAT, dimids(2), dimvarids(2))) 672 CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell centers")) 673 CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m")) 674 675 CALL check(nf90_def_var(ncid, "z", NF90_FLOAT, dimids(3), dimvarids(3))) 676 CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell centers")) 677 CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m")) 678 output_file % dimvarids_scl = dimvarids ! save dimvarids for later 679 680 ! overwrite third dimid with the one of depth 717 ! 718 !-- reset dimids first 719 dimids = (/0, 0, 0/) 720 CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) ) 721 CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) ) 722 CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) ) 723 ! 724 !-- save dimids for later 725 output_file % dimids_scl = dimids 726 727 ! 728 !-- reset dimvarids first 729 dimvarids = (/0, 0, 0/) 730 CALL check(nf90_def_var(ncid, "x", NF90_FLOAT, dimids(1), dimvarids(1))) 731 CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell centers")) 732 CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m")) 733 734 CALL check(nf90_def_var(ncid, "y", NF90_FLOAT, dimids(2), dimvarids(2))) 735 CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell centers")) 736 CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m")) 737 738 CALL check(nf90_def_var(ncid, "z", NF90_FLOAT, dimids(3), dimvarids(3))) 739 CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell centers")) 740 CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m")) 741 ! 742 !-- save dimvarids for later 743 output_file % dimvarids_scl = dimvarids 744 745 ! 746 !-- overwrite third dimid with the one of depth 681 747 CALL check(nf90_def_dim(ncid, "zsoil", SIZE(palm_grid % depths), dimids(3)) ) 682 output_file % dimids_soil = dimids ! save dimids for later 683 684 ! overwrite third dimvarid with the one of depth 748 ! 749 !-- save dimids for later 750 output_file % dimids_soil = dimids 751 752 ! 753 !-- overwrite third dimvarid with the one of depth 685 754 CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file % dimids_soil(3), dimvarids(3))) 686 755 CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "depth_below_land")) 687 756 CALL check(nf90_put_att(ncid, dimvarids(3), "positive", "down")) 688 757 CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m")) 689 output_file % dimvarids_soil = dimvarids ! save dimvarids for later 758 ! 759 !-- save dimvarids for later 760 output_file % dimvarids_soil = dimvarids 690 761 ! 691 762 !------------------------------------------------------------------------------ 692 763 !- Section 2b: Define dimensions for cell faces/velocities 693 764 !------------------------------------------------------------------------------ 694 dimids = (/0, 0, 0/) ! reset dimids 695 CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) ) 696 CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) ) 697 CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) ) 698 output_file % dimids_vel = dimids ! save dimids for later 699 700 dimvarids = (/0, 0, 0/) ! reset dimvarids 701 CALL check(nf90_def_var(ncid, "xu", NF90_FLOAT, dimids(1), dimvarids(1))) 702 CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell faces")) 703 CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m")) 704 705 CALL check(nf90_def_var(ncid, "yv", NF90_FLOAT, dimids(2), dimvarids(2))) 706 CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell faces")) 707 CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m")) 708 709 CALL check(nf90_def_var(ncid, "zw", NF90_FLOAT, dimids(3), dimvarids(3))) 710 CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell faces")) 711 CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m")) 712 output_file % dimvarids_vel = dimvarids ! save dimvarids for later 765 ! 766 !-- reset dimids first 767 dimids = (/0, 0, 0/) 768 CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) ) 769 CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) ) 770 CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) ) 771 ! 772 !-- save dimids for later 773 output_file % dimids_vel = dimids 774 775 ! 776 !-- reset dimvarids first 777 dimvarids = (/0, 0, 0/) 778 CALL check(nf90_def_var(ncid, "xu", NF90_FLOAT, dimids(1), dimvarids(1))) 779 CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell faces")) 780 CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m")) 781 782 CALL check(nf90_def_var(ncid, "yv", NF90_FLOAT, dimids(2), dimvarids(2))) 783 CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell faces")) 784 CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m")) 785 786 CALL check(nf90_def_var(ncid, "zw", NF90_FLOAT, dimids(3), dimvarids(3))) 787 CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell faces")) 788 CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m")) 789 ! 790 !-- save dimvarids for later 791 output_file % dimvarids_vel = dimvarids 713 792 714 793 ! … … 739 818 CALL check(nf90_put_var(ncid, output_file % dimvarids_vel(3), palm_grid%zw)) 740 819 741 ! TODO Read in soil depths from input file before this. 820 ! 821 !-- TODO Read in soil depths from input file before this. 742 822 CALL check(nf90_put_var(ncid, output_file % dimvarids_soil(3), palm_grid%depths)) 743 823 744 ! Write time vector 824 ! 825 !-- Write time vector 745 826 CALL check(nf90_put_var(ncid, output_file % dimvarid_time, output_file % time)) 746 827 747 ! Close the file 828 ! 829 !-- Close the file 748 830 CALL check(nf90_close(ncid)) 749 831 … … 751 833 752 834 835 !------------------------------------------------------------------------------! 836 ! Description: 837 ! ------------ 838 !> Defines the netCDF variables to be written to the dynamic driver file 839 !------------------------------------------------------------------------------! 753 840 SUBROUTINE setup_netcdf_variables(filename, output_variable_table, debug) 754 841 … … 824 911 !- Section 1: Load input buffers for accumulated variables 825 912 !------------------------------------------------------------------------------ 913 ! 914 !-- radiation budgets, precipitation 826 915 IF (group % kind == 'running average' .OR. & 827 group % kind == 'accumulated') THEN ! radiation budgets, precipitation916 group % kind == 'accumulated') THEN 828 917 829 918 IF (SIZE(group % in_var_list) .GT. 1 ) THEN … … 835 924 END IF 836 925 837 ! use two buffer arrays 926 ! 927 !-- use two buffer arrays 838 928 nbuffers = 2 839 929 IF ( .NOT. ALLOCATED( buffer ) ) ALLOCATE( buffer(nbuffers) ) 840 930 841 ! chose correct buffer array 842 hour = iter - 1! hour of the day 931 ! 932 !-- hour of the day 933 hour = iter - 1 934 ! 935 !-- chose correct buffer array 843 936 buf_id = select_buffer(hour) 844 937 … … 859 952 ELSE 860 953 861 ! Allocate one input buffer per input_variable. If more quantities 862 ! have to be computed than input variables exist in this group, 863 ! allocate more buffers. For instance, in the thermodynamics group, 864 ! there are three input variabels (PP, T, Qv) and four quantities 865 ! necessart (P, Theta, Rho, qv) for the corresponding output fields 866 ! (p0, Theta, qv, ug, and vg) 954 ! 955 !-- Allocate one input buffer per input_variable. If more quantities 956 !-- have to be computed than input variables exist in this group, 957 !-- allocate more buffers. For instance, in the thermodynamics group, 958 !-- there are three input variabels (PP, T, Qv) and four quantities 959 !-- necessart (P, Theta, Rho, qv) for the corresponding output fields 960 !-- (p0, Theta, qv, ug, and vg) 867 961 nbuffers = MAX( group % n_inputs, group % n_output_quantities ) 868 962 ALLOCATE( buffer(nbuffers) ) 869 963 CALL run_control('time', 'alloc') 870 964 871 ! Read in all input variables, leave extra buffers-if any-untouched. 965 ! 966 !-- Read in all input variables, leave extra buffers-if any-untouched. 872 967 DO ivar = 1, group % n_inputs 873 968 874 969 input_var => group % in_var_list(ivar) 875 970 876 ! Check wheather P or PP is present in input file 971 ! 972 ! Check wheather P or PP is present in input file 877 973 IF (input_var % name == 'P') THEN 878 974 input_var % name = TRIM( get_pressure_varname(input_file) ) … … 891 987 892 988 989 !------------------------------------------------------------------------------! 990 ! Description: 991 ! ------------ 992 !> Select the appropriate buffer ID for accumulated COSMO input variables 993 !> depending on the current hour. 994 !------------------------------------------------------------------------------! 893 995 INTEGER FUNCTION select_buffer(hour) 894 996 INTEGER, INTENT(IN) :: hour … … 943 1045 944 1046 1047 !------------------------------------------------------------------------------! 1048 ! Description: 1049 ! ------------ 1050 !> Read the given global attribute form the given netCDF file. 1051 !------------------------------------------------------------------------------! 945 1052 FUNCTION get_netcdf_attribute(filename, attribute) RESULT(attribute_value) 946 1053 … … 948 1055 REAL(dp) :: attribute_value 949 1056 950 INTEGER :: ncid1057 INTEGER :: ncid 951 1058 952 1059 IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR ) THEN … … 966 1073 967 1074 1075 1076 !------------------------------------------------------------------------------! 1077 ! Description: 1078 ! ------------ 1079 !> Updates the dynamic driver with the interpolated field of the current 1080 !> variable at the current time step. 1081 !------------------------------------------------------------------------------! 968 1082 SUBROUTINE update_output(var, array, iter, output_file, cfg) 969 1083 TYPE(nc_var), INTENT(IN) :: var … … 980 1094 ) 981 1095 982 ! Skip time dimension for output 1096 ! 1097 !-- Skip time dimension for output 983 1098 ndim = var % ndim 984 1099 IF ( var_is_time_dependent ) ndim = var % ndim - 1 … … 990 1105 CALL check(nf90_open(output_file % name, NF90_WRITE, ncid)) 991 1106 992 ! Reduce dimension of output array according to variable kind 1107 ! 1108 !-- Reduce dimension of output array according to variable kind 993 1109 SELECT CASE (TRIM(var % kind)) 994 1110 … … 1086 1202 1087 1203 1088 SUBROUTINE write_netcdf_variable_2d(var, iter, output_file, buffer) 1089 TYPE(nc_var), INTENT(IN) :: var 1090 INTEGER, INTENT(IN) :: iter 1091 TYPE(nc_file), INTENT(IN) :: output_file 1092 REAL(dp), INTENT(IN) :: buffer(:,:,:) 1093 1094 INTEGER :: ncid, ndim_out, start(4), count(4) 1095 LOGICAL :: last_dimension_is_time 1096 1097 ndim_out = var % ndim 1098 1099 last_dimension_is_time = var % dimids( var % ndim ) == output_file % dimid_time 1100 IF ( last_dimension_is_time ) THEN 1101 ndim_out = ndim_out - 1 1102 END IF 1103 1104 start(:) = (/1,1,1,iter/) 1105 count(1:ndim_out) = var%dimlen(1:ndim_out) 1106 1107 CALL check(nf90_open(output_file % name, NF90_WRITE, ncid)) 1108 1109 IF (TRIM(var % kind) .EQ. 'left/right scalar') THEN 1110 1111 CALL check(nf90_put_var( ncid, var%varid, buffer(1,:,:) ) ) 1112 1113 ELSE IF (TRIM(var % kind) .EQ. 'north/south scalar') THEN 1114 1115 CALL check(nf90_put_var( ncid, var%varid, buffer(:,1,:) ) ) 1116 1117 ELSE IF (TRIM(var % kind) .EQ. 'top scalar') THEN 1118 1119 CALL check(nf90_put_var( ncid, var%varid, buffer(:,:,1) ) ) 1120 ELSE 1121 1122 CALL check(nf90_put_var( ncid, var%varid, buffer ) ) 1123 1124 END IF 1125 CALL check(nf90_close(ncid)) 1126 1127 END SUBROUTINE write_netcdf_variable_2d 1128 1129 1204 !------------------------------------------------------------------------------! 1205 ! Description: 1206 ! ------------ 1207 !> Checks the status of a netCDF API call and aborts if an error occured 1208 !------------------------------------------------------------------------------! 1130 1209 SUBROUTINE check(status) 1131 1210 -
palm/trunk/UTIL/inifor/src/inifor_transform.f90
r3537 r3557 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Updated documentation 29 ! 30 ! 31 ! 3537 2018-11-20 10:53:14Z eckhard 28 32 ! bugfix: working precision added 29 33 ! … … 52 56 ! Authors: 53 57 ! -------- 54 ! @author Eckhard Kadasch58 !> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach) 55 59 ! 56 60 ! Description: … … 110 114 DO k = nz, LBOUND(out_arr, 3), -1 111 115 112 ! TODO: Remove IF clause and extrapolate based on a critical vertical 113 ! TODO: index marking the lower bound of COSMO-DE data coverage. 114 ! Check for negative interpolation weights indicating grid points 115 ! below COSMO-DE domain and extrapolate from the top in such cells. 116 ! 117 !-- TODO: Remove IF clause and extrapolate based on a critical vertical 118 !-- TODO: index marking the lower bound of COSMO-DE data coverage. 119 !-- Check for negative interpolation weights indicating grid points 120 !-- below COSMO-DE domain and extrapolate from the top in such cells. 116 121 IF (outgrid % w_verti(i,j,k,1) < -1.0_dp .AND. k < nz) THEN 117 122 out_arr(i,j,k) = out_arr(i,j,k+1) … … 155 160 !------------------------------------------------------------------------------! 156 161 SUBROUTINE interpolate_2d(invar, outvar, outgrid, ncvar) 157 ! I index 0-based for the indices of the outvar to be consistent with the 158 ! outgrid indices and interpolation weights. 162 ! 163 !-- I index 0-based for the indices of the outvar to be consistent with the 164 !-- outgrid indices and interpolation weights. 159 165 TYPE(grid_definition), INTENT(IN) :: outgrid 160 166 REAL(dp), INTENT(IN) :: invar(0:,0:,0:) … … 164 170 INTEGER :: i, j, k, l 165 171 166 ! TODO: check if input dimensions are consistent, i.e. ranges are correct 167 IF (UBOUND(outvar, 3) .GT. UBOUND(invar, 3)) THEN 172 ! 173 !-- TODO: check if input dimensions are consistent, i.e. ranges are correct 174 IF ( UBOUND(outvar, 3) .GT. UBOUND(invar, 3) ) THEN 168 175 message = "Output array for '" // TRIM(ncvar % name) // "' has ' more levels (" // & 169 176 TRIM(str(UBOUND(outvar, 3))) // ") than input variable ("//& … … 190 197 191 198 199 !------------------------------------------------------------------------------! 200 ! Description: 201 ! ------------ 202 !> Compute the horizontal average of the in_arr(:,:,:) and store it in 203 !> out_arr(:) 204 !------------------------------------------------------------------------------! 192 205 SUBROUTINE average_2d(in_arr, out_arr, ii, jj) 193 206 REAL(dp), INTENT(IN) :: in_arr(0:,0:,0:) … … 200 213 IF (SIZE(ii) .NE. SIZE(jj)) THEN 201 214 message = "Length of 'ii' and 'jj' index lists do not match." // & 202 NEW_LINE(' ') // "ii has " // str(SIZE(ii)) // " elements, " // 215 NEW_LINE(' ') // "ii has " // str(SIZE(ii)) // " elements, " // & 203 216 NEW_LINE(' ') // "jj has " // str(SIZE(jj)) // "." 204 217 CALL abort('average_2d', message) … … 211 224 i = ii(l) 212 225 j = jj(l) 213 out_arr(k) = out_arr(k) +& 214 in_arr(i, j, k) 226 out_arr(k) = out_arr(k) + in_arr(i, j, k) 215 227 END DO 216 228 … … 223 235 224 236 237 !------------------------------------------------------------------------------! 238 ! Description: 239 ! ------------ 240 !> Three-dimensional interpolation driver. Interpolates from the source_array to 241 !> the given palm_grid. 242 !> 243 !> The routine separates horizontal and vertical 244 !> interpolation. In the horizontal interpolation step, the source_array data is 245 !> interpolated along COSMO grid levels onto the intermediate grid (vertically 246 !> as coarse as COSMO, horizontally as fine as PALM). 247 !------------------------------------------------------------------------------! 225 248 SUBROUTINE interpolate_3d(source_array, palm_array, palm_intermediate, palm_grid) 226 249 TYPE(grid_definition), INTENT(IN) :: palm_intermediate, palm_grid … … 232 255 nx = palm_intermediate % nx 233 256 ny = palm_intermediate % ny 234 nlev = palm_intermediate % nz ! nlev 235 236 ! Interpolate from COSMO-DE to intermediate grid. Allocating with one 237 ! less point in the vertical, since scalars like T have 50 instead of 51 238 ! points in COSMO-DE. 257 nlev = palm_intermediate % nz 258 259 ! 260 !-- Interpolate from COSMO to intermediate grid. Allocating with one 261 !-- less point in the vertical, since scalars like T have 50 instead of 51 262 !-- points in COSMO. 239 263 ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nlev-1)) ! 240 264 241 265 CALL interpolate_2d(source_array, intermediate_array, palm_intermediate) 242 266 243 ! Interpolate from intermediate grid to palm_grid grid, includes 244 ! extrapolation for cells below COSMO-DE domain. 267 ! 268 !-- Interpolate from intermediate grid to palm_grid grid, includes 269 !-- extrapolation for cells below COSMO domain. 245 270 CALL interpolate_1d(intermediate_array, palm_array, palm_grid) 246 271 … … 250 275 251 276 277 !------------------------------------------------------------------------------! 278 ! Description: 279 ! ------------ 280 !> Average data horizontally from the source_array over the region given by the 281 !> averaging grid 'avg_grid' and store the result in 'profile_array'. 282 !------------------------------------------------------------------------------! 252 283 SUBROUTINE average_profile(source_array, profile_array, avg_grid) 253 284 TYPE(grid_definition), INTENT(IN) :: avg_grid … … 265 296 j_source = avg_grid % jjj(l) 266 297 267 DO k_profile = avg_grid % k_min, UBOUND(profile_array, 1) ! PALM levels 268 269 DO m = 1, 2 ! vertical interpolation neighbours 298 ! 299 !-- Loop over PALM levels 300 DO k_profile = avg_grid % k_min, UBOUND(profile_array, 1) 301 302 ! 303 !-- Loop over vertical interpolation neighbours 304 DO m = 1, 2 270 305 271 306 k_source = avg_grid % kkk(l, k_profile, m) … … 274 309 + avg_grid % w(l, k_profile, m) & 275 310 * source_array(i_source, j_source, k_source) 276 277 END DO ! m, vertical interpolation neighbours 278 279 END DO ! k_profile, PALM levels 280 281 END DO ! l, horizontal neighbours 311 ! 312 !-- Loop over vertical interpolation neighbours m 313 END DO 314 315 ! 316 !-- Loop over PALM levels k_profile 317 END DO 318 319 ! 320 !-- Loop over horizontal neighbours l 321 END DO 282 322 283 323 ni_columns = 1.0_dp / avg_grid % n_columns 284 324 profile_array(:) = profile_array(:) * ni_columns 285 325 286 ! Extrapolate constant to the bottom 326 ! 327 !-- Constant extrapolation to the bottom 287 328 profile_array(1:avg_grid % k_min-1) = profile_array(avg_grid % k_min) 288 329 … … 290 331 291 332 333 !------------------------------------------------------------------------------! 334 ! Description: 335 ! ------------ 336 !> Extrapolates density linearly from the level 'k_min' downwards. 337 !------------------------------------------------------------------------------! 292 338 SUBROUTINE extrapolate_density(rho, avg_grid) 293 339 REAL(dp), DIMENSION(:), INTENT(INOUT) :: rho … … 308 354 309 355 356 !------------------------------------------------------------------------------! 357 ! Description: 358 ! ------------ 359 !> Driver for extrapolating pressure from PALM level k_min downwards 360 !------------------------------------------------------------------------------! 310 361 SUBROUTINE extrapolate_pressure(p, rho, avg_grid) 311 362 REAL(dp), DIMENSION(:), INTENT(IN) :: rho … … 405 456 REAL(dp) :: ri 406 457 407 ! TODO check dimensions of lat/lon and y/x match 458 ! 459 !-- TODO check dimensions of lat/lon and y/x match 408 460 409 461 ri = 1.0_dp / r … … 438 490 REAL(dp) :: ri 439 491 440 ! If this elemental function is called with a large array as xy, it is 441 ! computationally more efficient to precompute the inverse radius and 442 ! then muliply. 492 ! 493 !-- If this elemental function is called with a large array as xy, it is 494 !-- computationally more efficient to precompute the inverse radius and 495 !-- then muliply. 443 496 ri = 1.0_dp / r 444 497 … … 448 501 449 502 503 !------------------------------------------------------------------------------! 504 ! Description: 505 ! ------------ 506 !> For a rotated-pole system with the origin at geographical latitude 'phi_c', 507 !> compute the geographical latitude of its rotated north pole. 508 !------------------------------------------------------------------------------! 450 509 REAL(dp) FUNCTION phic_to_phin(phi_c) 451 510 REAL(dp), INTENT(IN) :: phi_c … … 456 515 457 516 517 !------------------------------------------------------------------------------! 518 ! Description: 519 ! ------------ 520 !> For a rotated-pole system with the origin at geographical latitude 'phi_c' 521 !> and longitude 'lam_c', compute the geographical longitude of its rotated 522 !> north pole. 523 !------------------------------------------------------------------------------! 458 524 REAL(dp) FUNCTION lamc_to_lamn(phi_c, lam_c) 459 525 REAL(dp), INTENT(IN) :: phi_c, lam_c … … 467 533 468 534 535 !------------------------------------------------------------------------------! 536 ! Description: 537 ! ------------ 538 !> Set gamma according to whether PALM domain is in the northern or southern 539 !> hemisphere of the COSMO rotated-pole system. Gamma assumes either the 540 !> value 0 or PI and is needed to work around around a bug in the 541 !> rotated-pole coordinate transformations. 542 !------------------------------------------------------------------------------! 469 543 REAL(dp) FUNCTION gamma_from_hemisphere(phi_cg, phi_ref) 470 REAL(dp), INTENT(IN) :: phi_cg, phi_ref 471 LOGICAL :: palm_centre_is_south_of_cosmo_origin 472 473 palm_centre_is_south_of_cosmo_origin = (phi_cg < phi_ref) 474 475 IF (palm_centre_is_south_of_cosmo_origin) THEN 544 REAL(dp), INTENT(IN) :: phi_cg 545 REAL(dp), INTENT(IN) :: phi_ref 546 547 LOGICAL :: palm_origin_is_south_of_cosmo_origin 548 549 palm_origin_is_south_of_cosmo_origin = (phi_cg < phi_ref) 550 551 IF (palm_origin_is_south_of_cosmo_origin) THEN 476 552 gamma_from_hemisphere = PI 477 553 ELSE … … 550 626 551 627 628 !------------------------------------------------------------------------------! 629 ! Description: 630 ! ------------ 631 !> Rotate the given vector field (x(:), y(:)) by the given 'angle'. 632 !------------------------------------------------------------------------------! 552 633 SUBROUTINE rotate_vector_field(x, y, angle) 553 634 REAL(dp), DIMENSION(:), INTENT(INOUT) :: x, y !< x and y coodrinate in arbitrary units … … 559 640 sine = SIN(angle * TO_RADIANS) 560 641 cosine = COS(angle * TO_RADIANS) 561 ! RESAHPE() fills columns first, so the rotation matrix becomes 562 ! 563 ! rotation = [ cosine -sine ] 564 ! [ sine cosine ] 642 ! 643 !-- RESAHPE() fills columns first, so the rotation matrix becomes 644 !-- 645 !-- rotation = [ cosine -sine ] 646 !-- [ sine cosine ] 565 647 rotation = RESHAPE( (/cosine, sine, -sine, cosine/), (/2, 2/) ) 566 648 … … 588 670 !> Datenbanken des DWD. 589 671 !> https://www.dwd.de/SharedDocs/downloads/DE/modelldokumentationen/nwv/cosmo_d2/cosmo_d2_dbbeschr_aktuell.pdf?__blob=publicationFile&v=2 590 ! >672 !------------------------------------------------------------------------------! 591 673 FUNCTION meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g) & 592 674 RESULT(delta) … … 664 746 DO j = 0, UBOUND(palm_clon, 2)!palm_grid % ny 665 747 DO i = 0, UBOUND(palm_clon, 1)!palm_grid % nx 666 ! Compute the floating point index corrseponding to PALM-4U grid point 667 ! location along target grid (COSMO-DE) axes. 748 ! 749 !-- Compute the floating point index corrseponding to PALM-4U grid point 750 !-- location along target grid (COSMO-DE) axes. 668 751 lonpos = (palm_clon(i,j) - lon0) * cosmo_dxi 669 752 latpos = (palm_clat(i,j) - lat0) * cosmo_dyi … … 693 776 694 777 778 !------------------------------------------------------------------------------! 779 ! Description: 780 ! ------------ 781 !> Computes linear vertical interpolation neighbour indices and weights for each 782 !> column of the given palm grid. 783 !------------------------------------------------------------------------------! 695 784 SUBROUTINE find_vertical_neighbours_and_weights_interp( palm_grid, & 696 785 palm_intermediate ) … … 709 798 nlev = palm_intermediate % nz 710 799 711 ! in each column of the fine grid, find vertical neighbours of every cell 800 ! 801 !-- in each column of the fine grid, find vertical neighbours of every cell 712 802 DO j = 0, ny 713 803 DO i = 0, nx … … 718 808 column_top = palm_intermediate % h(i,j,nlev) 719 809 720 ! scan through palm_grid column and set neighbour indices in 721 ! case current_height is either below column_base, in the current 722 ! cell, or above column_top. Keep increasing current cell index until 723 ! the current cell overlaps with the current_height. 810 ! 811 !-- scan through palm_grid column and set neighbour indices in 812 !-- case current_height is either below column_base, in the current 813 !-- cell, or above column_top. Keep increasing current cell index until 814 !-- the current cell overlaps with the current_height. 724 815 DO k = 1, nz 725 816 726 ! Memorize the top and bottom boundaries of the coarse cell and the 727 ! current height within it 817 ! 818 !-- Memorize the top and bottom boundaries of the coarse cell and the 819 !-- current height within it 728 820 current_height = palm_grid % z(k) + palm_grid % z0 729 821 h_top = palm_intermediate % h(i,j,k_intermediate+1) … … 738 830 ) 739 831 740 ! set default weights 832 ! 833 !-- set default weights 741 834 palm_grid % w_verti(i,j,k,1:2) = 0.0_dp 742 835 … … 755 848 756 849 ELSE 757 ! cycle through intermediate levels until current 758 ! intermediate-grid cell overlaps with current_height 850 ! 851 !-- cycle through intermediate levels until current 852 !-- intermediate-grid cell overlaps with current_height 759 853 DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1) 760 854 k_intermediate = k_intermediate + 1 … … 777 871 palm_grid % kk(i,j,k,2) = k_intermediate + 1 778 872 779 ! copmute vertical weights 873 ! 874 !-- compute vertical weights 780 875 weight = (h_top - current_height) / (h_top - h_bottom) 781 876 palm_grid % w_verti(i,j,k,1) = weight … … 791 886 792 887 888 !------------------------------------------------------------------------------! 889 ! Description: 890 ! ------------ 891 !> Computes linear vertical interpolation neighbour indices and weights for each 892 !> column of the given averaging grid. 893 !> 894 !> The difference to the _interp variant of this routine lies in how columns 895 !> are adressed. While the _interp variant loops over all PALM grid columns 896 !> given by combinations of all index indices (i,j), this routine loops over a 897 !> subset of COSMO columns, which are sequantlially stored in the index lists 898 !> iii(:) and jjj(:). 899 !------------------------------------------------------------------------------! 793 900 SUBROUTINE find_vertical_neighbours_and_weights_average( avg_grid ) 794 901 TYPE(grid_definition), INTENT(INOUT) :: avg_grid … … 805 912 nlev = SIZE(avg_grid % cosmo_h, 3) 806 913 807 ! in each column of the fine grid, find vertical neighbours of every cell 914 ! 915 !-- in each column of the fine grid, find vertical neighbours of every cell 808 916 DO l = 1, avg_grid % n_columns 809 917 … … 814 922 column_top = avg_grid % cosmo_h(i,j,nlev) 815 923 816 ! scan through avg_grid column until and set neighbour indices in 817 ! case current_height is either below column_base, in the current 818 ! cell, or above column_top. Keep increasing current cell index until 819 ! the current cell overlaps with the current_height. 924 ! 925 !-- scan through avg_grid column until and set neighbour indices in 926 !-- case current_height is either below column_base, in the current 927 !-- cell, or above column_top. Keep increasing current cell index until 928 !-- the current cell overlaps with the current_height. 820 929 k_intermediate = 1 !avg_grid % cosmo_h is indezed 1-based. 821 930 DO k_palm = 1, avg_grid % nz 822 931 823 ! Memorize the top and bottom boundaries of the coarse cell and the 824 ! current height within it 932 ! 933 !-- Memorize the top and bottom boundaries of the coarse cell and the 934 !-- current height within it 825 935 current_height = avg_grid % z(k_palm) + avg_grid % z0 826 936 h_top = avg_grid % cosmo_h(i,j,k_intermediate+1) 827 937 h_bottom = avg_grid % cosmo_h(i,j,k_intermediate) 828 938 829 point_is_above_grid = (current_height > column_top) !22000m, very unlikely 939 ! 940 !-- COSMO column top is located at 22000m, point_is_above_grid is very 941 !-- unlikely. 942 point_is_above_grid = (current_height > column_top) 830 943 point_is_below_grid = (current_height < column_base) 831 944 … … 835 948 ) 836 949 837 ! set default weights 950 ! 951 !-- set default weights 838 952 avg_grid % w(l,k_palm,1:2) = 0.0_dp 839 953 … … 852 966 avg_grid % k_min = MAX(k_palm + 1, avg_grid % k_min) 853 967 ELSE 854 ! cycle through intermediate levels until current 855 ! intermediate-grid cell overlaps with current_height 968 ! 969 !-- cycle through intermediate levels until current 970 !-- intermediate-grid cell overlaps with current_height 856 971 DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1) 857 972 k_intermediate = k_intermediate + 1 … … 865 980 END DO 866 981 867 ! k_intermediate = 48 indicates the last section (indices 48 and 49), i.e. 868 ! k_intermediate = 49 is not the beginning of a valid cell. 982 ! 983 !-- k_intermediate = 48 indicates the last section (indices 48 and 49), i.e. 984 !-- k_intermediate = 49 is not the beginning of a valid cell. 869 985 IF (k_intermediate > nlev-1) THEN 870 986 message = "Index " // TRIM(str(k_intermediate)) // & … … 876 992 avg_grid % kkk(l,k_palm,2) = k_intermediate + 1 877 993 878 ! copmute vertical weights 994 ! 995 !-- compute vertical weights 879 996 weight = (h_top - current_height) / (h_top - h_bottom) 880 997 avg_grid % w(l,k_palm,1) = weight … … 882 999 END IF 883 1000 884 END DO ! k, PALM levels 885 END DO ! l, averaging columns 1001 ! 1002 !-- Loop over PALM levels k 1003 END DO 1004 1005 ! 1006 !-- Loop over averaging columns l 1007 END DO 886 1008 887 1009 END SUBROUTINE find_vertical_neighbours_and_weights_average … … 931 1053 ! ii(i,j,1/2) ii(i,j,3/4) 932 1054 ! 1055 !------------------------------------------------------------------------------! 933 1056 SUBROUTINE compute_horizontal_interp_weights(cosmo_lat, cosmo_lon, & 934 1057 palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz) … … 950 1073 DO i = 0, UBOUND(palm_clon, 1) 951 1074 952 ! weight in lambda direction 1075 ! 1076 !-- weight in lambda direction 953 1077 wl = ( cosmo_lon(palm_ii(i,j,4)) - palm_clon(i,j) ) * cosmo_dxi 954 1078 955 ! weight in phi direction 1079 ! 1080 !-- weight in phi direction 956 1081 wp = ( cosmo_lat(palm_jj(i,j,2)) - palm_clat(i,j) ) * cosmo_dyi 957 1082 … … 988 1113 !> COSMO-DE the velocity points are staggared one half grid spaceing up-grid 989 1114 !> which means the first centre point has to be omitted and is set to zero. 1115 !------------------------------------------------------------------------------! 990 1116 SUBROUTINE centre_velocities(u_face, v_face, u_centre, v_centre) 991 1117 REAL(dp), DIMENSION(0:,0:,0:), INTENT(IN) :: u_face, v_face … … 1004 1130 1005 1131 1132 !------------------------------------------------------------------------------! 1133 ! Description: 1134 ! ------------ 1135 !> Compute the geographical latitude of a point given in rotated-pole cordinates 1136 !------------------------------------------------------------------------------! 1006 1137 FUNCTION phirot2phi (phirot, rlarot, polphi, pollam, polgam) 1007 1138 … … 1041 1172 1042 1173 1174 !------------------------------------------------------------------------------! 1175 ! Description: 1176 ! ------------ 1177 !> Compute the geographical latitude of a point given in rotated-pole cordinates 1178 !------------------------------------------------------------------------------! 1043 1179 FUNCTION phi2phirot (phi, rla, polphi, pollam) 1044 1180 … … 1072 1208 1073 1209 1210 !------------------------------------------------------------------------------! 1211 ! Description: 1212 ! ------------ 1213 !> Compute the geographical longitude of a point given in rotated-pole cordinates 1214 !------------------------------------------------------------------------------! 1074 1215 FUNCTION rlarot2rla(phirot, rlarot, polphi, pollam, polgam) 1075 1216 … … 1123 1264 1124 1265 1266 !------------------------------------------------------------------------------! 1267 ! Description: 1268 ! ------------ 1269 !> Compute the rotated-pole longitude of a point given in geographical cordinates 1270 !------------------------------------------------------------------------------! 1125 1271 FUNCTION rla2rlarot ( phi, rla, polphi, pollam, polgam ) 1126 1272 … … 1162 1308 1163 1309 1310 !------------------------------------------------------------------------------! 1311 ! Description: 1312 ! ------------ 1313 !> Rotate the given velocity vector (u,v) from the geographical to the 1314 !> rotated-pole system 1315 !------------------------------------------------------------------------------! 1164 1316 SUBROUTINE uv2uvrot(u, v, rlat, rlon, pollat, pollon, urot, vrot) 1165 1317 … … 1187 1339 1188 1340 1341 !------------------------------------------------------------------------------! 1342 ! Description: 1343 ! ------------ 1344 !> Rotate the given velocity vector (urot, vrot) from the rotated-pole to the 1345 !> geographical system 1346 !------------------------------------------------------------------------------! 1189 1347 SUBROUTINE uvrot2uv (urot, vrot, rlat, rlon, pollat, pollon, u, v) 1190 1348 -
palm/trunk/UTIL/inifor/src/inifor_types.f90
r3447 r3557 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Updated documentation 29 ! 30 ! 31 ! 3447 2018-10-29 15:52:54Z eckhard 28 32 ! Renamed source files for compatibilty with PALM build system 29 33 ! … … 51 55 ! Authors: 52 56 ! -------- 53 ! @author Eckhard Kadasch57 !> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach) 54 58 ! 55 59 ! Description: … … 66 70 IMPLICIT NONE 67 71 72 !------------------------------------------------------------------------------! 73 ! Description: 74 ! ------------ 75 !> Contaner for the INIFOR command-line configuration 76 !------------------------------------------------------------------------------! 68 77 TYPE inifor_config 69 78 CHARACTER(LEN=DATE) :: start_date !< String of the FORMAT YYYYMMDDHH indicating the start of the intended PALM-4U simulation … … 82 91 CHARACTER(LEN=SNAME) :: soilmoisture_prefix !< Prefix of input files for soil moisture spin-up, e.g 'laf' for COSMO-DE analyses 83 92 84 CHARACTER(LEN=SNAME) :: averaging_mode 85 CHARACTER(LEN=SNAME) :: bc_mode 86 CHARACTER(LEN=SNAME) :: ic_mode 87 CHARACTER(LEN=SNAME) :: rotation_method 88 89 REAL(dp) :: p0 90 REAL(dp) :: ug 91 REAL(dp) :: vg 92 REAL(dp) :: z0 !< Elevation of the PALM-4U domain above sea level [m]93 CHARACTER(LEN=SNAME) :: averaging_mode !< destinguishes between level-based and heigh-based averaging 94 CHARACTER(LEN=SNAME) :: bc_mode !< destinguishes realistic and idealistic forcing 95 CHARACTER(LEN=SNAME) :: ic_mode !< destinguishes volume and profile initialization 96 CHARACTER(LEN=SNAME) :: rotation_method !< selects method for velocity rotation 97 98 REAL(dp) :: p0 !< manually specified surface pressure [Pa] 99 REAL(dp) :: ug !< manually spefied geostrophic wind component in x direction [m/s] 100 REAL(dp) :: vg !< manually spefied geostrophic wind component in y direction [m/s] 101 REAL(dp) :: z0 !< elevation of the PALM-4U domain above sea level [m] 93 102 REAL(dp) :: averaging_angle !< latitudal and longitudal width of averaging regions [deg] 94 103 95 96 LOGICAL :: debug 97 LOGICAL :: p0_is_set 98 LOGICAL :: ug_is_set 99 LOGICAL :: vg_is_set 100 LOGICAL :: flow_prefix_is_set !< 101 LOGICAL :: input_prefix_is_set !< 102 LOGICAL :: radiation_prefix_is_set !< 103 LOGICAL :: soil_prefix_is_set !< 104 LOGICAL :: soilmoisture_prefix_is_set !< 104 LOGICAL :: debug !< indicates whether --debug option was given 105 LOGICAL :: p0_is_set !< indicates whether p0 was set manually 106 LOGICAL :: ug_is_set !< indicates whether ug was set manually 107 LOGICAL :: vg_is_set !< indicates whether vg was set manually 108 LOGICAL :: flow_prefix_is_set !< indicates whether the flow prefix was set manually 109 LOGICAL :: input_prefix_is_set !< indicates whether the input prefix was set manually 110 LOGICAL :: radiation_prefix_is_set !< indicates whether the radiation prefix was set manually 111 LOGICAL :: soil_prefix_is_set !< indicates whether the soil prefix was set manually 112 LOGICAL :: soilmoisture_prefix_is_set !< indicates whether the soilmoisture prefix was set manually 105 113 END TYPE inifor_config 106 114 115 116 !------------------------------------------------------------------------------! 117 ! Description: 118 ! ------------ 119 !> Container for grid data, in partucular coordinates, interpolation neighbours 120 !> and weights 121 !------------------------------------------------------------------------------! 107 122 TYPE grid_definition 108 123 CHARACTER(LEN=SNAME) :: name(3) !< names of the grid dimensions, e.g. (/'x', 'y', 'z'/) or (/'latitude', 'longitude', 'height'/) … … 152 167 153 168 169 !------------------------------------------------------------------------------! 170 ! Description: 171 ! ------------ 172 !> Container for name and dimensions of the netCDF output file 173 !------------------------------------------------------------------------------! 154 174 TYPE nc_file 155 CHARACTER(LEN=PATH) :: name !< file name156 INTEGER :: dimid_time !< NetCDF IDs of the time dimension157 INTEGER :: dimids_scl(3) !< NetCDF IDs of the grid dimensions for scalar points x, y, z158 INTEGER :: dimids_vel(3) !< NetCDF IDs of the grid dimensions for velocity points xu, yu, zu159 INTEGER :: dimids_soil(3) !< NetCDF IDs of the grid dimensions for soil points x, y, depth160 INTEGER :: dimvarid_time !< NetCDF IDs of the time variable161 INTEGER :: dimvarids_scl(3) !< NetCDF IDs of the grid coordinates of scalars x, y, z162 INTEGER :: dimvarids_vel(3) !< NetCDF IDs of the grid coordinates of velocities xu, yu, zu. Note that velocities are located at mix of both coordinates, e.g. u(xu, y, z).163 INTEGER :: dimvarids_soil(3) !< NetCDF IDs of the grid coordinates for soil points x, y, depth164 REAL(dp), POINTER :: time(:) !vector of output time steps175 CHARACTER(LEN=PATH) :: name !< file name 176 INTEGER :: dimid_time !< NetCDF IDs of the time dimension 177 INTEGER :: dimids_scl(3) !< NetCDF IDs of the grid dimensions for scalar points x, y, z 178 INTEGER :: dimids_vel(3) !< NetCDF IDs of the grid dimensions for velocity points xu, yu, zu 179 INTEGER :: dimids_soil(3) !< NetCDF IDs of the grid dimensions for soil points x, y, depth 180 INTEGER :: dimvarid_time !< NetCDF IDs of the time variable 181 INTEGER :: dimvarids_scl(3) !< NetCDF IDs of the grid coordinates of scalars x, y, z 182 INTEGER :: dimvarids_vel(3) !< NetCDF IDs of the grid coordinates of velocities xu, yu, zu. Note that velocities are located at mix of both coordinates, e.g. u(xu, y, z). 183 INTEGER :: dimvarids_soil(3) !< NetCDF IDs of the grid coordinates for soil points x, y, depth 184 REAL(dp), POINTER :: time(:) !< vector of output time steps 165 185 END TYPE nc_file 166 186 167 187 188 !------------------------------------------------------------------------------! 189 ! Description: 190 ! ------------ 191 !> Metadata container for netCDF variables 192 !------------------------------------------------------------------------------! 168 193 TYPE nc_var 169 194 INTEGER :: varid !< NetCDF ID of the variable … … 193 218 194 219 195 TYPE io_group !< Input/Output group, groups together output variabels that share their input variables. For instance, all boundary surfaces and initialization fields of the potential temperature are base on T and p. 220 !------------------------------------------------------------------------------! 221 ! Description: 222 ! ------------ 223 !> Input/Output group, groups together nc_var-type output variabels that share 224 !> input variables as well as lists of the netCDF files they are stored in. 225 !> For instance, all boundary surfaces and initialization fields of the 226 !> potential temperature are base on the input netCDF variables T and P. 227 !------------------------------------------------------------------------------! 228 TYPE io_group 196 229 INTEGER :: nt !< maximum number of output time steps across all output variables 197 230 INTEGER :: nv !< number of netCDF output variables … … 208 241 209 242 243 !------------------------------------------------------------------------------! 244 ! Description: 245 ! ------------ 246 !> Container for input data arrays. read_input_variables() allocates a 247 !> one-dimensional array of containers, to accomodate all inputs of the given 248 !> IO group in one variable. 249 !------------------------------------------------------------------------------! 210 250 TYPE container 211 REAL(dp), ALLOCATABLE :: array(:,:,:) 212 LOGICAL :: is_preprocessed = .FALSE. 251 REAL(dp), ALLOCATABLE :: array(:,:,:) !< generic data array 252 LOGICAL :: is_preprocessed = .FALSE. !< flag indicating whether input array has been preprocessed 213 253 END TYPE container 214 254 -
palm/trunk/UTIL/inifor/src/inifor_util.f90
r3447 r3557 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Updated documentation 29 ! 30 ! 31 ! 3447 2018-10-29 15:52:54Z eckhard 28 32 ! Renamed source files for compatibilty with PALM build system 29 33 ! … … 44 48 ! Authors: 45 49 ! -------- 46 ! @author Eckhard Kadasch50 !> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach) 47 51 ! 48 52 ! Description: … … 61 65 IMPLICIT NONE 62 66 67 !------------------------------------------------------------------------------! 68 ! Description: 69 ! ------------ 70 !> Fortran implementation of C's struct tm for representing points in time 71 !------------------------------------------------------------------------------! 63 72 TYPE, BIND(c) :: tm_struct 64 73 INTEGER(C_INT) :: tm_sec !< seconds after the minute [0, 61] … … 75 84 INTERFACE 76 85 86 !------------------------------------------------------------------------------! 87 ! Description: 88 ! ------------ 89 !> Interface to C's strptime function, which converts a string to a tm 90 !> structure. 91 !------------------------------------------------------------------------------! 77 92 FUNCTION strptime(string, format, timeinfo) BIND(c, NAME='strptime') 78 93 IMPORT :: C_CHAR, C_SIZE_T, tm_struct … … 87 102 88 103 104 !------------------------------------------------------------------------------! 105 ! Description: 106 ! ------------ 107 !> Interface to C's strftime function, which converts the given 'timeinfo' 108 !> structure to a string in the given 'format'. 109 !------------------------------------------------------------------------------! 89 110 FUNCTION strftime(string, string_len, format, timeinfo) BIND(c, NAME='strftime') 90 111 IMPORT :: C_CHAR, C_SIZE_T, tm_struct … … 101 122 102 123 124 !------------------------------------------------------------------------------! 125 ! Description: 126 ! ------------ 127 !> Interface to C's mktime function, which converts the given tm structure to 128 !> Unix time (number of seconds since 0 UTC, 1st January 1970). INIFOR uses the 129 !> side effect that a call to mktime noramlizes the given 'timeinfo' structure, 130 !> e.g. increments the date if hours overfow 24. 131 !------------------------------------------------------------------------------! 103 132 FUNCTION mktime(timeinfo) BIND(c, NAME='mktime') 104 133 IMPORT :: C_PTR, tm_struct … … 115 144 CONTAINS 116 145 146 !------------------------------------------------------------------------------! 147 ! Description: 148 ! ------------ 149 !> Takes a string of the form YYYYMMDDHH, adds the given number of hours to its 150 !> tm structure representation, and returns the corresponding string in the same 151 !> format 152 !------------------------------------------------------------------------------! 117 153 CHARACTER(LEN=DATE) FUNCTION add_hours_to(date_string, hours) 118 154 CHARACTER(LEN=DATE), INTENT(IN) :: date_string … … 143 179 144 180 181 !------------------------------------------------------------------------------! 182 ! Description: 183 ! ------------ 184 !> Print all members of the given tm structure 185 !------------------------------------------------------------------------------! 145 186 SUBROUTINE print_tm(timeinfo) 146 187 TYPE(tm_struct), INTENT(IN) :: timeinfo … … 158 199 159 200 201 !------------------------------------------------------------------------------! 202 ! Description: 203 ! ------------ 204 !> Initialize the given tm structure with zero values 205 !------------------------------------------------------------------------------! 160 206 SUBROUTINE init_tm(timeinfo) 161 207 TYPE(tm_struct), INTENT(INOUT) :: timeinfo … … 177 223 178 224 179 SUBROUTINE fake_output_3d(a) 180 181 REAL(dp), INTENT(INOUT) :: a(:,:,:) 182 REAL(dp) :: lxi, lyi 183 INTEGER :: i,j,k 184 185 lyi = 2.0_dp * PI / (SIZE(a, 2) - 1.0_dp) 186 lxi = 2.0_dp * PI / (SIZE(a, 1) - 1.0_dp) 187 188 DO k = 1, SIZE(a, 3) 189 DO j = 1, SIZE(a, 2) 190 DO i = 1, SIZE(a, 1) 191 a(i,j,k) = SIN(lxi * i) * COS(lyi * j) + k 192 END DO 193 END DO 194 END DO 195 196 END SUBROUTINE fake_output_3d 197 198 199 SUBROUTINE fake_output_2d(a, offset) 200 201 REAL(dp), INTENT(INOUT) :: a(:,:) 202 INTEGER, INTENT(IN) :: offset 203 REAL(dp) :: lxi, lyi 204 INTEGER :: i,j 205 206 lyi = 2.0_dp*PI / (SIZE(a, 2) - 1.0_dp) 207 lxi = 2.0_dp*PI / (SIZE(a, 1) - 1.0_dp) 208 209 a(:,:) = 1.0_dp 210 DO j = 1, SIZE(a, 2) 211 DO i = 1, SIZE(a, 1) 212 a(i,j) = SIN(lxi * i) * COS(lyi * j) + offset 213 END DO 214 END DO 215 216 END SUBROUTINE fake_output_2d 217 218 225 !------------------------------------------------------------------------------! 226 ! Description: 227 ! ------------ 228 !> Fill the given array with values equally spaced between and including start 229 !> and stop 230 !------------------------------------------------------------------------------! 219 231 SUBROUTINE linspace(start, stop, array) 220 232 … … 255 267 256 268 269 !------------------------------------------------------------------------------! 270 ! Description: 271 ! ------------ 272 !> 273 !------------------------------------------------------------------------------! 257 274 SUBROUTINE deaverage(avg_1, t1, avg_2, t2, avg_3, t3) 258 275 … … 270 287 271 288 289 !------------------------------------------------------------------------------! 290 ! Description: 291 ! ------------ 292 !> Compute the COSMO-DE/-D2 basic state pressure profile 293 !------------------------------------------------------------------------------! 272 294 SUBROUTINE get_basic_state(z, beta, p_sl, t_sl, rd, g, p0) 273 295 … … 278 300 REAL(dp), INTENT(IN) :: rd !< ideal gas constant of dry air [J/kg/K] 279 301 REAL(dp), INTENT(IN) :: g !< acceleration of Earth's gravity [m/s^2] 280 REAL(dp), INTENT(OUT) :: p0(1:) !< COSMO -DEbasic state pressure [Pa]302 REAL(dp), INTENT(OUT) :: p0(1:) !< COSMO basic state pressure [Pa] 281 303 REAL(dp) :: root_frac, factor !< precomputed factors 282 304 … … 285 307 286 308 p0(:) = p_sl * EXP( & 287 factor * ( 1.0_dp - SQRT( 1.0_dp - root_frac * z(:) ) ) &288 309 factor * ( 1.0_dp - SQRT( 1.0_dp - root_frac * z(:) ) ) & 310 ) 289 311 290 312 END SUBROUTINE get_basic_state … … 345 367 346 368 369 !------------------------------------------------------------------------------! 370 ! Description: 371 ! ------------ 372 !> Converts the given real value to a string 373 !------------------------------------------------------------------------------! 347 374 CHARACTER(LEN=16) FUNCTION real_to_str_f(val) 348 375 … … 355 382 356 383 384 !------------------------------------------------------------------------------! 385 ! Description: 386 ! ------------ 387 !> Converts the given integer value to a string 388 !------------------------------------------------------------------------------! 357 389 CHARACTER(LEN=10) FUNCTION str(val) 358 390 … … 365 397 366 398 367 CHARACTER(LEN=30) FUNCTION strs(vals) 368 369 INTEGER, INTENT(IN) :: vals(:) 370 INTEGER :: i 371 372 strs = '' 373 DO i = 1, SIZE(vals) 374 strs = TRIM(strs) // TRIM(str(vals(i))) 375 END DO 376 377 END FUNCTION strs 378 379 399 !------------------------------------------------------------------------------! 400 ! Description: 401 ! ------------ 402 !> If the given path is not conlcuded by a slash, add one. 403 !------------------------------------------------------------------------------! 380 404 SUBROUTINE normalize_path(path) 381 405
Note: See TracChangeset
for help on using the changeset viewer.
