Changeset 3557 for palm/trunk/UTIL


Ignore:
Timestamp:
Nov 22, 2018 4:01:22 PM (6 years ago)
Author:
eckhard
Message:

inifor: Updated documentation

Location:
palm/trunk/UTIL/inifor/src
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/UTIL/inifor/src/inifor.f90

    r3537 r3557  
    2626! -----------------
    2727! $Id$
     28! Updated documentation
     29!
     30! 3537 2018-11-20 10:53:14Z eckhard
    2831! Print version number on program start
    2932!
     
    5457! Authors:
    5558! --------
    56 ! @author Eckhard Kadasch
     59!> @author Eckhard Kadasch, (Deutscher Wetterdienst, Offenbach)
    5760!
    5861! Description:
     
    8487    IMPLICIT NONE
    8588   
    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
    108115   
    109116!> \mainpage About INIFOR
     
    116123    CALL report('main_loop', 'Running INIFOR version ' // VERSION)
    117124
    118     ! Initialize INIFOR's parameters from command-line interface and namelists
     125!
     126!-- Initialize INIFOR's parameters from command-line interface and namelists
    119127    CALL setup_parameters()
    120128
    121     ! Initialize all grids, including interpolation neighbours and weights
     129!
     130!-- Initialize all grids, including interpolation neighbours and weights
    122131    CALL setup_grids()
    123132 CALL run_control('time', 'init')
    124133
    125     ! Initialize the netCDF output file and define dimensions
     134!
     135!-- Initialize the netCDF output file and define dimensions
    126136    CALL setup_netcdf_dimensions(output_file, palm_grid, cfg % start_date,    &
    127137                                 origin_lon, origin_lat)
    128138 CALL run_control('time', 'write')
    129139
    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
    132143    CALL setup_variable_tables(cfg % ic_mode)
    133144 CALL run_control('time', 'write')
    134145
    135     ! Add the output variables to the netCDF output file
     146!
     147!-- Add the output variables to the netCDF output file
    136148    CALL setup_netcdf_variables(output_file % name, output_var_table,          &
    137149                                cfg % debug)
     
    143155!- Section 2: Main loop
    144156!------------------------------------------------------------------------------
    145 
    146157    DO igroup = 1, SIZE(io_group_list)
    147158
     
    345356 CALL run_control('time', 'comp')
    346357
    347 
    348                    ! This case gets called twice, the first time for ug, the
    349                    ! second time for vg. We compute ug and vg at the first call
    350                    ! and keep vg (and ug for that matter) around for the second
    351                    ! 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.
    352363                   CASE ( 'geostrophic winds' )
    353364
     
    360371                            vg_arr = cfg % vg
    361372                         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,    &
    367379                                                    ug_arr, vg_arr )
    368380                         END IF
     
    372384                      END IF
    373385
    374                       ! Prepare output of geostrophic winds
     386!
     387!--                   Prepare output of geostrophic winds
    375388                      SELECT CASE(TRIM(output_var % name))
    376389                      CASE ('ls_forcing_ug')
     
    397410                      SELECT CASE (TRIM(output_var % name))
    398411
    399                       !CASE('ls_forcing_ug')
    400                       !    output_arr(1, 1, :) = ug
    401 
    402                       !CASE('ls_forcing_vg')
    403                       !    output_arr(1, 1, :) = vg
    404 
    405412                      CASE('nudging_tau')
    406413                          output_arr(1, 1, :) = NUDGING_TAU
     
    418425                                "has not been implemented, yet."
    419426                      CALL abort('main loop', message)
    420                       !ALLOCATE( output_arr( 1, 1, 1:nz ) )
    421427
    422428                   CASE DEFAULT
     
    444450                END IF
    445451
    446              END DO ! ouput variables
     452!
     453!--          output variable loop
     454             END DO
    447455
    448456             ug_vg_have_been_computed = .FALSE.
     
    463471             END IF
    464472
     473!
     474!--          Keep input buffer around for averaged (radiation) and
     475!--          accumulated COSMO quantities (precipitation).
    465476             IF ( group % kind == 'running average' .OR. &
    466477                  group % kind == 'accumulated' )  THEN
    467                 ! Keep input buffer around for averaged (radiation) and
    468                 ! accumulated COSMO-DE quantities (precipitation).
    469478             ELSE
    470479                CALL report('main loop', 'Deallocating input buffer', cfg % debug)
     
    473482 CALL run_control('time', 'alloc')
    474483
    475           END DO ! time steps / input files
     484!
     485!--       time steps / input files loop
     486          END DO
    476487
    477488          IF (ALLOCATED(input_buffer))  THEN
     
    491502          CALL report('main loop', message, cfg % debug)
    492503
    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
    496511
    497512!------------------------------------------------------------------------------
  • palm/trunk/UTIL/inifor/src/inifor_control.f90

    r3447 r3557  
    2626! -----------------
    2727! $Id$
     28! Updated documentation
     29!
     30!
     31! 3447 2018-10-29 15:52:54Z eckhard
    2832! Renamed source files for compatibilty with PALM build system
    2933!
     
    4448! Authors:
    4549! --------
    46 ! @author Eckhard Kadasch
     50!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
    4751!
    4852! Description:
     
    6165    IMPLICIT NONE
    6266
    63     CHARACTER (LEN=5000) ::  message = ''
     67    CHARACTER (LEN=5000) ::  message = '' !< log message buffer
    6468
    6569 CONTAINS
    6670
     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!------------------------------------------------------------------------------!
    6782    SUBROUTINE report(routine, message, debug)
    6883
    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
    7894          OPEN( NEWUNIT=u, FILE='inifor.log', STATUS='replace' )
    7995          is_first_run = .FALSE.
     
    84100
    85101       suppress_message = .FALSE.
    86        IF (PRESENT(debug))  THEN
    87           IF (.NOT. debug)  suppress_message = .TRUE.
     102       IF ( PRESENT(debug) )  THEN
     103          IF ( .NOT. debug )  suppress_message = .TRUE.
    88104       END IF
    89105
    90        IF (.NOT. suppress_message)  THEN
     106       IF ( .NOT. suppress_message )  THEN
    91107          PRINT *, "inifor: " // TRIM(message) // "  [ " // TRIM(routine) // " ]"
    92108          WRITE(u, *)  TRIM(message) // "  [ " // TRIM(routine) // " ]"
     
    98114
    99115
     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!------------------------------------------------------------------------------!
    100126    SUBROUTINE warn(routine, message)
    101127
    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
    104130
    105131       CALL report(routine, "WARNING: " // TRIM(message))
     
    108134
    109135
     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!------------------------------------------------------------------------------!
    110146    SUBROUTINE abort(routine, message)
    111147
    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
    114150
    115151       CALL report(routine, "ERROR: " // TRIM(message) // " Stopping.")
     
    119155
    120156
     157!------------------------------------------------------------------------------!
     158! Description:
     159! ------------
     160!>
     161!> print_version() prints the INIFOR version number and copyright notice.
     162!------------------------------------------------------------------------------!
    121163    SUBROUTINE print_version()
    122164       PRINT *, "INIFOR " // VERSION
     
    125167
    126168
     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!------------------------------------------------------------------------------!
    127176    SUBROUTINE run_control(mode, budget)
    128177
    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
    138191
    139192
  • palm/trunk/UTIL/inifor/src/inifor_defs.f90

    r3537 r3557  
    2626! -----------------
    2727! $Id$
     28! Updated documentation
     29!
     30!
     31! 3537 2018-11-20 10:53:14Z eckhard
    2832! Bumped version number
    2933!
     
    5963! Authors:
    6064! --------
    61 ! @author Eckhard Kadasch
     65!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
    6266!
    6367! Description:
     
    6872 
    6973 IMPLICIT NONE
    70  
    71  ! Parameters for type definitions
     74
     75!
     76!-- Parameters for type definitions
    7277 INTEGER, PARAMETER  ::  dp    = 8   !< double precision (8 bytes = 64 bits)
    7378 INTEGER, PARAMETER  ::  sp    = 4   !< single precision (4 bytes = 32 bits)
     
    7883 INTEGER, PARAMETER  ::  DATE  = 10  !< length of date strings
    7984
    80  ! Trigonomentry
     85!
     86!-- Trigonomentry
    8187 REAL(dp), PARAMETER ::  PI = 3.14159265358979323846264338_dp !< Ratio of a circle's circumference to its diamter [-]
    8288 REAL(dp), PARAMETER ::  TO_RADIANS = PI / 180.0_dp           !< Conversion factor from degrees to radiant [-]
    8389 REAL(dp), PARAMETER ::  TO_DEGREES = 180.0_dp / PI           !< Conversion factor from radians to degrees [-]
    8490
    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
    96107
    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]
    102114
    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 [-]
    105119 INTEGER, PARAMETER          ::  FORCING_STEP = 1             !< Number of hours between forcing time steps [h]
    106120 REAL(dp), PARAMETER         ::  NUDGING_TAU = 21600.0_dp     !< Nudging relaxation time scale [s]
  • palm/trunk/UTIL/inifor/src/inifor_grid.f90

    r3537 r3557  
    2626! -----------------
    2727! $Id$
     28! Updated documentation
     29!
     30!
     31! 3537 2018-11-20 10:53:14Z eckhard
    2832! Read COSMO domain extents and soil depths from input files
    2933! Report averaging mode and debugging mode in log
     
    6973! Authors:
    7074! --------
    71 ! @author Eckhard Kadasch
     75!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
    7276!
    7377!------------------------------------------------------------------------------!
     
    193197    INTEGER ::  step_hour !< number of hours between forcing time steps
    194198
    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
    200203
    201204    TYPE(nc_var), ALLOCATABLE, TARGET ::  input_var_table(:)  !< table of input variables
    202205    TYPE(nc_var), ALLOCATABLE, TARGET ::  output_var_table(:) !< table of input variables
    203     TYPE(nc_var)                      ::  cosmo_var           !< COSMO-DE dummy variable, used for reading HHL, rlon, rlat
     206    TYPE(nc_var)                      ::  cosmo_var           !< COSMO dummy variable, used for reading HHL, rlon, rlat
    204207
    205208    TYPE(grid_definition), TARGET ::  palm_grid                       !< PALM-4U grid in the target system (COSMO-DE rotated-pole)
    206209    TYPE(grid_definition), TARGET ::  palm_intermediate               !< PALM-4U grid with coarse vertical grid wiht levels interpolated from COSMO-DE grid
    207210    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
    261263
    262264    TYPE(io_group), ALLOCATABLE, TARGET ::  io_group_list(:)  !< List of I/O groups, which group together output variables that share the same input variable
    263265 
    264266    NAMELIST /inipar/ nx, ny, nz, dx, dy, dz, longitude, latitude,             &
    265                       dz_max, dz_stretch_factor, dz_stretch_level,             & !< old grid stretching parameters
    266                       dz_stretch_level_start, dz_stretch_level_end               !< new grid stretching parameters
     267                      dz_max, dz_stretch_factor, dz_stretch_level,             &
     268                      dz_stretch_level_start, dz_stretch_level_end
    267269    NAMELIST /d3par/  end_time
    268270   
    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'
    283287                         
    284     TYPE(nc_file) ::  output_file
    285 
    286     TYPE(inifor_config) ::  cfg !< Container of INIFOR configuration
     288    TYPE(nc_file) ::  output_file !< metadata of the dynamic driver
     289
     290    TYPE(inifor_config) ::  cfg !< container of the INIFOR command-line configuration
    287291
    288292 CONTAINS
    289293   
     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!------------------------------------------------------------------------------!
    290302    SUBROUTINE setup_parameters()
    291        INTEGER ::  k
     303       INTEGER ::  k !< loop index
    292304
    293305!
     
    300312       start_hour_soilmoisture = - (4 * 7 * 24) - 2
    301313
    302        ! COSMO-DE default rotated pole
     314!
     315!--    COSMO-DE and -D2 rotated pole position
    303316       phi_n     =   40.0_dp * TO_RADIANS
    304317       phi_equat =   50.0_dp * TO_RADIANS
    305318       lambda_n  = -170.0_dp * TO_RADIANS
    306319
    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)
    308322       origin_lat = 52.325079_dp * TO_RADIANS ! south-west of Berlin, origin used for the Dec 2017 showcase simulation
    309323       origin_lon = 13.082744_dp * TO_RADIANS
    310324       cfg % z0 = 35.0_dp
    311325
    312        ! Default atmospheric parameters
     326!
     327!--    Default atmospheric parameters
    313328       cfg % ug = 0.0_dp
    314329       cfg % vg = 0.0_dp
    315330       cfg % p0 = P_SL
    316331
    317        ! Parameters for file names
     332!
     333!--    Parameters for file names
    318334       start_hour_flow = 0
    319335       start_hour_soil = 0
     
    323339       step_hour = FORCING_STEP
    324340
    325        input_prefix = 'laf'  ! 'laf' for COSMO-DE analyses
     341       input_prefix = 'laf'
    326342       cfg % flow_prefix = input_prefix
    327343       cfg % input_prefix = input_prefix
     
    342358!------------------------------------------------------------------------------
    343359
    344        ! Set default paths and modes
     360!
     361!--    Set default paths and modes
    345362       cfg % input_path         = './'
    346363       cfg % hhl_file           = ''
     
    353370       cfg % averaging_mode = 'level'
    354371
    355        ! Overwrite defaults with user configuration
     372!
     373!--    Overwrite defaults with user configuration
    356374       CALL parse_command_line_arguments( cfg )
    357375
     
    380398       END IF
    381399
    382 
    383        ! Set default file paths, if not specified by user.
     400!
     401!--    Set default file paths, if not specified by user.
    384402       CALL normalize_path(cfg % input_path)
    385403       IF (TRIM(cfg % hhl_file) == '')  cfg % hhl_file = TRIM(cfg % input_path) // 'hhl.nc'
     
    399417
    400418 CALL run_control('time', 'init')
    401        ! Read in namelist parameters
     419 !
     420 !--   Read in namelist parameters
    402421       OPEN(10, FILE=cfg % namelist_file)
    403422       READ(10, NML=inipar) ! nx, ny, nz, dx, dy, dz
     
    408427       end_hour = CEILING( end_time / 3600.0 * step_hour )
    409428
    410        ! Generate input file lists
     429!
     430!--    Generate input file lists
    411431       CALL get_input_file_list(                                               &
    412432          cfg % start_date, start_hour_flow, end_hour, step_hour,              &
     
    437457
    438458 CALL run_control('time', 'init')
    439        ! Read COSMO-DE soil type map
     459!
     460!--    Read COSMO soil type map
    440461       cosmo_var % name = 'SOILTYP'
    441462       CALL get_netcdf_variable(cfg % soiltyp_file, cosmo_var, soiltyp)
     
    465486       CALL report('setup_parameters', message)
    466487
    467 
    468        ! Read in COSMO-DE heights of half layers (vertical cell faces)
     488!
     489!--    Read in COSMO heights of half layers (vertical cell faces)
    469490       cosmo_var % name = 'HHL'
    470491       CALL get_netcdf_variable(cfg % hhl_file, cosmo_var, hhl)
     
    486507 CALL run_control('time', 'comp')
    487508
    488        ! Appoximate COSMO-DE heights of full layers (cell centres)
     509!
     510!--    Appoximate COSMO-DE heights of full layers (cell centres)
    489511       ALLOCATE( hfl(nlon, nlat, nlev-1) )
    490512       ALLOCATE( d_depth(ndepths), d_depth_rho_inv(ndepths) )
     
    494516       d_depth_rho_inv = 1.0_dp / ( d_depth * RHO_L )
    495517
    496        ! Appoximate COSMO-DE heights of full layers (cell centres)
     518!
     519!--    Appoximate COSMO-DE heights of full layers (cell centres)
    497520       DO k = 1, nlev-1
    498521          hfl(:,:,k) = 0.5_dp * ( hhl(:,:,k) +                                 &
     
    506529! Section 4.2: PALM-4U parameters
    507530!------------------------------------------------------------------------------
    508        ! PALM-4U domain extents
     531!
     532!--    PALM-4U domain extents
    509533       lx = (nx+1) * dx
    510534       ly = (ny+1) * dy
    511535       
    512        ! PALM-4U point of Earth tangency
     536!
     537!--    PALM-4U point of Earth tangency
    513538       x0 = 0.0_dp
    514539       y0 = 0.0_dp
    515540
    516        ! time vector
     541!
     542!--    time vector
    517543       nt = CEILING(end_time / (step_hour * 3600.0_dp)) + 1
    518544       ALLOCATE( time(nt) )
     
    521547 CALL run_control('time', 'init')
    522548
    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
    524551       phi_c    = TO_RADIANS *                                                 &
    525552                  phi2phirot( origin_lat * TO_DEGREES, origin_lon * TO_DEGREES,&
     
    530557                              0.0_dp )
    531558
    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.
    536564       gam = gamma_from_hemisphere(origin_lat, phi_equat)
    537565
    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.
    541570       phi_cn    = phic_to_phin(phi_c)
    542571       lambda_cn = lamc_to_lamn(phi_c, lambda_c)
     
    573602!------------------------------------------------------------------------------
    574603
    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
    576606    phi_centre = phirot2phi(                                                   &
    577607       phirot = project(0.5_dp*ly, y0, EARTH_RADIUS) * TO_DEGREES,             &
     
    595625    CALL report( 'setup_parameters', message )
    596626
    597 
    598     ! Compute boundaries of the central averaging box
     627!
     628!-- Compute boundaries of the central averaging box
    599629    averaging_angle = cfg % averaging_angle * TO_RADIANS
    600630    lam_east = lam_centre + 0.5_dp * averaging_angle
     
    605635    averaging_width_ns = averaging_angle * EARTH_RADIUS
    606636
    607     ! Coriolis parameter
     637!
     638!-- Coriolis parameter
    608639    f3 = 2.0_dp * OMEGA * SIN(                                                 &
    609640       TO_RADIANS*phirot2phi( phi_centre * TO_DEGREES, lam_centre * TO_DEGREES,&
     
    618649! Description:
    619650! ------------
    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
    621653!------------------------------------------------------------------------------!
    622     SUBROUTINE setup_grids() ! setup_grids(inifor_settings(with nx, ny, nz,...))
     654    SUBROUTINE setup_grids()
    623655       CHARACTER ::  interp_mode
    624656
     
    626658! Section 0: Define base PALM-4U grid coordinate vectors
    627659!------------------------------------------------------------------------------
    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.
    632664       ALLOCATE( x(0:nx), y(0:ny), z(1:nz), z_column(1:nz+1) )
    633665       CALL linspace(0.5_dp * dx, lx - 0.5_dp * dx, x)
     
    642674       z_top = z_column(nz+1)
    643675
    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.
    646679       ALLOCATE( xu(1:nx),  yv(1:ny), zw(1:nz-1), zw_column(1:nz))
    647680       CALL linspace(dx, lx - dx, xu)
     
    661694               nx=nx, ny=ny, nz=nz, z=z, zw=zw, ic_mode=cfg % ic_mode)
    662695
    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.
    664698       CALL init_grid_definition('cosmo-de', grid=cosmo_grid,                  &
    665699               xmin=lonmin, xmax=lonmax,                                       &
     
    668702               nx=nlon-1, ny=nlat-1, nz=nlev-1)
    669703
    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.
    674709       CALL init_grid_definition('palm intermediate', grid=palm_intermediate,  &
    675710               xmin=0.0_dp, xmax=lx,                                           &
     
    10071042               latmin = phi_south, latmax = phi_north,                         &
    10081043               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)  THEN
    1014        !   ! Here, I use the 'boundary' kind, so the vertical coordinate uses the
    1015        !   ! 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 IF
    10401044
    10411045!                                                                             
     
    10911095
    10921096
     1097!------------------------------------------------------------------------------!
     1098! Description:
     1099! ------------
     1100!> Driver for computing neighbour indices and weights for horizontal and
     1101!> vertical interpolation.
     1102!------------------------------------------------------------------------------!
    10931103    SUBROUTINE setup_interpolation(cosmo_grid, grid, intermediate_grid, kind, ic_mode)
    10941104
     
    11061116! Section 1: Horizontal interpolation                                       
    11071117!------------------------------------------------------------------------------
    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)
    11091120       SELECT CASE(kind)
    11101121
    1111        CASE('s') ! scalars
     1122!
     1123!--    scalars
     1124       CASE('s')
    11121125
    11131126          cosmo_lat => cosmo_grid % lat
    11141127          cosmo_lon => cosmo_grid % lon
    11151128          cosmo_h   => cosmo_grid % hfl
    1116 
    1117        CASE('w') ! vertical velocity
     1129!
     1130!--    vertical velocity
     1131       CASE('w')
    11181132
    11191133          cosmo_lat => cosmo_grid % lat
    11201134          cosmo_lon => cosmo_grid % lon
    11211135          cosmo_h   => cosmo_grid % hhl
    1122 
    1123        CASE('u') ! x velocity
     1136!
     1137!--    x velocity
     1138       CASE('u')
    11241139
    11251140          cosmo_lat => cosmo_grid % lat
     
    11271142          cosmo_h   => cosmo_grid % hfl
    11281143
    1129        CASE('v') ! y velocity
     1144!
     1145!--    y velocity
     1146       CASE('v')
    11301147
    11311148          cosmo_lat => cosmo_grid % latv
     
    11531170!------------------------------------------------------------------------------
    11541171
    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.
    11611179       IF (PRESENT(ic_mode))  THEN
    11621180          IF (TRIM(ic_mode) == 'profile')  setup_volumetric = .FALSE.
     
    11691187          intermediate_grid % h(:,:,:) = - EARTH_RADIUS
    11701188
    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
    11731192          CALL interpolate_2d(cosmo_h, intermediate_grid % h, intermediate_grid)
    11741193          CALL find_vertical_neighbours_and_weights_interp(grid, intermediate_grid)
     
    12461265           grid % z => z
    12471266
    1248            ! Allocate neighbour indices and weights
     1267!
     1268!--        Allocate neighbour indices and weights
    12491269           IF (TRIM(ic_mode) .NE. 'profile')  THEN
    12501270              ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) )
     
    12741294           )
    12751295
    1276            ! Allocate neighbour indices and weights
     1296!
     1297!--        Allocate neighbour indices and weights
    12771298           ALLOCATE( grid % ii(0:nx, 0:ny, 4),                                 &
    12781299                     grid % jj(0:nx, 0:ny, 4) )
     
    12831304           grid % w_horiz(:,:,:)   = 0.0_dp
    12841305       
    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.
    12871309        CASE('palm')
    12881310
     
    13011323           grid % name(3) = 'z'
    13021324
    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.
    13051328           ALLOCATE( grid % x(0:nx),   grid % y(0:ny) )
    13061329           ALLOCATE( grid % xu(1:nx),  grid % yv(1:ny) )
     
    13141337           grid % depths => depths
    13151338
    1316            ! Allocate neighbour indices and weights
     1339!
     1340!--        Allocate neighbour indices and weights
    13171341           IF (TRIM(ic_mode) .NE. 'profile')  THEN
    13181342              ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) )
     
    13291353           grid % name(3) = 'interpolated hhl or hfl'
    13301354
    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.
    13331358           ALLOCATE( grid % x(0:nx),   grid % y(0:ny) )
    13341359           ALLOCATE( grid % xu(1:nx),  grid % yv(1:ny) )
     
    13401365           grid % depths => depths
    13411366
    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
    13431369           ALLOCATE( grid % clon(0:nx, 0:ny),   grid % clat(0:nx, 0:ny)  )
    13441370           ALLOCATE( grid % clonu(1:nx, 0:ny),  grid % clatu(1:nx, 0:ny) )
    13451371           ALLOCATE( grid % clonv(0:nx, 1:ny),  grid % clatv(0:nx, 1:ny) )
    13461372
    1347            ! Compute rotated-pole coordinates of...
    1348            ! ... PALM-4U centres
     1373!
     1374!--        Compute rotated-pole coordinates of...
     1375!--        ... PALM-4U centres
    13491376           CALL rotate_to_cosmo(                                               &
    13501377              phir = project( grid % y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude
     
    13561383           )
    13571384
    1358            ! ... PALM-4U u winds
     1385!
     1386!--        ... PALM-4U u winds
    13591387           CALL rotate_to_cosmo(                                               &
    13601388              phir = project( grid % y,  y0, EARTH_RADIUS ), & ! = plate-carree latitude
     
    13661394           )
    13671395
    1368            ! ... PALM-4U v winds
     1396!
     1397!--        ... PALM-4U v winds
    13691398           CALL rotate_to_cosmo(                                               &
    13701399              phir = project( grid % yv, y0, EARTH_RADIUS ), & ! = plate-carree latitude
     
    13761405           )
    13771406
    1378            ! Allocate neighbour indices and weights
     1407!
     1408!--        Allocate neighbour indices and weights
    13791409           ALLOCATE( grid % ii(0:nx, 0:ny, 4),                                 &
    13801410                     grid % jj(0:nx, 0:ny, 4) )
     
    13981428           grid % latv(:) = grid % lat + 0.5_dp * (grid % ly / grid % ny)
    13991429
    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
    14021433           grid % hhl => hhl
    14031434           grid % hfl => hfl
     
    14721503       avg_grid % kind = TRIM(kind)
    14731504
    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
    14761508       CALL get_cosmo_averaging_region(avg_grid, cosmo_grid)
    14771509
    14781510       ALLOCATE (avg_grid % kkk(avg_grid % n_columns, avg_grid % nz, 2) )
    14791511       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
    14811514       SELECT CASE(avg_grid % kind)
    14821515
     
    14941527       END SELECT
    14951528
    1496        ! For level-besed averaging, compute average heights
     1529!
     1530!--    For level-besed averaging, compute average heights
    14971531       level_based_averaging = .TRUE.
    14981532       IF (level_based_averaging)  THEN
     
    15031537       END IF
    15041538
    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.
    15081543       IF ( TRIM(cfg % averaging_mode) == 'level' )  THEN
    15091544          FORALL(                                                              &
     
    15131548       END IF
    15141549
    1515        ! Compute vertical weights and neighbours
     1550!
     1551!--    Compute vertical weights and neighbours
    15161552       CALL find_vertical_neighbours_and_weights_average(avg_grid)
    15171553
     
    15501586       END SELECT
    15511587
    1552 
    1553        ! FIXME: make dlon, dlat parameters of the grid_defintion type
     1588!
     1589!--    FIXME: make dlon, dlat parameters of the grid_defintion type
    15541590       dlon = cosmo_lon(1) - cosmo_lon(0)
    15551591       dlat = cosmo_lat(1) - cosmo_lat(0)
     
    16101646       REAL(dp) ::  dz_level_end, dz_stretched
    16111647
    1612        INTEGER ::  dz_stretch_level_end_index(9)               !< vertical grid level index until which the vertical grid spacing is stretched
    1613        INTEGER ::  dz_stretch_level_start_index(9)             !< vertical grid level index above which the vertical grid spacing is stretched
     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
    16141650       INTEGER ::  dz_stretch_level_index = 0
    16151651       INTEGER ::  k, n, number_dz
     1652
    16161653!
    16171654!-- Compute height of u-levels from constant grid length and dz stretch factors
     
    18261863
    18271864
     1865!------------------------------------------------------------------------------!
    18281866! Description: [PALM subroutine]
    18291867! -----------------------------------------------------------------------------!
     
    19361974             
    19371975!
    1938 !--                stretch_factor_1 is taken to guarantee that the stretching
    1939 !--                procedure ends as close as possible to dz_stretch_level_end.
    1940 !--                stretch_factor_2 would guarantee that the stretched dz(n) is
    1941 !--                equal to dz(n+1) after l_rounded grid levels.
     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.
    19421980             IF (delta_total_new < delta_total_old) THEN
    19431981                dz_stretch_factor_array(n) = stretch_factor_1
     
    19742012 END SUBROUTINE calculate_stretching_factor
    19752013
     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!------------------------------------------------------------------------------!
    19762020    SUBROUTINE midpoints(z, zw)
    19772021
     
    19872031    END SUBROUTINE midpoints
    19882032
    1989 
     2033!------------------------------------------------------------------------------!
     2034! Description:
     2035! ------------
     2036!> Defines INFOR's IO groups.
     2037!------------------------------------------------------------------------------!
    19902038    SUBROUTINE setup_io_groups()
    19912039
     
    19942042       ngroups = 16
    19952043       ALLOCATE( io_group_list(ngroups) )
    1996        
    1997        !soil temp
     2044
     2045!
     2046!--    soil temp
    19982047       io_group_list(1) = init_io_group(                                       &
    19992048          in_files = soil_files,                                               &
     
    20032052       )
    20042053
    2005        !soil water
     2054!
     2055!--    soil water
    20062056       io_group_list(2) = init_io_group(                                       &
    20072057          in_files = soil_files,                                               &
     
    20112061       )
    20122062
    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
    20152066       io_group_list(3) = init_io_group(                                       &
    20162067          in_files = flow_files,                                               &
     
    20252076       )
    20262077
    2027        !Moved to therodynamic io_group
     2078!
     2079!--    Moved to therodynamic io_group
    20282080       !io_group_list(4) = init_io_group(                                       &
    20292081       !   in_files = flow_files,                                               &
     
    20332085       !)
    20342086
    2035        !u and v velocity, ug anv vg moved to thermodynamic io group.
     2087!
     2088!--    u and v velocity
    20362089       io_group_list(5) = init_io_group(                                       &
    20372090          in_files = flow_files,                                               &
     
    20432096       )
    20442097   
    2045        !!v velocity, deprecated!
     2098!
     2099!--    v velocity, deprecated!
    20462100       !io_group_list(6) = init_io_group(                                       &
    20472101       !   in_files = flow_files,                                               &
     
    20522106       !io_group_list(6) % to_be_processed = .FALSE.
    20532107   
    2054        !w velocity and subsidence and w nudging
     2108!
     2109!--    w velocity and subsidence and w nudging
    20552110       io_group_list(7) = init_io_group(                                       &
    20562111          in_files = flow_files,                                               &
     
    20592114          kind = 'scalar'                                                      &
    20602115       )
    2061 
    2062        !rain
     2116!
     2117!--    rain
    20632118       io_group_list(8) = init_io_group(                                       &
    20642119          in_files = soil_moisture_files,                                      &
     
    20682123       )
    20692124       io_group_list(8) % to_be_processed = .FALSE.
    2070 
    2071        !snow
     2125!
     2126!--    snow
    20722127       io_group_list(9) = init_io_group(                                       &
    20732128          in_files = soil_moisture_files,                                      &
     
    20772132       )
    20782133       io_group_list(9) % to_be_processed = .FALSE.
    2079 
    2080        !graupel
     2134!
     2135!--    graupel
    20812136       io_group_list(10) = init_io_group(                                      &
    20822137          in_files = soil_moisture_files,                                      &
     
    20862141       )
    20872142       io_group_list(10) % to_be_processed = .FALSE.
    2088 
    2089        !evapotranspiration
     2143!
     2144!--    evapotranspiration
    20902145       io_group_list(11) = init_io_group(                                      &
    20912146          in_files = soil_moisture_files,                                      &
     
    20952150       )
    20962151       io_group_list(11) % to_be_processed = .FALSE.
    2097 
    2098        !2m air temperature
     2152!
     2153!--    2m air temperature
    20992154       io_group_list(12) = init_io_group(                                      &
    21002155          in_files = soil_moisture_files,                                      &
     
    21032158          kind = 'surface'                                                     &
    21042159       )
    2105 
    2106        !incoming diffusive sw flux
     2160!
     2161!--    incoming diffusive sw flux
    21072162       io_group_list(13) = init_io_group(                                      &
    21082163          in_files = radiation_files,                                          &
     
    21122167       )
    21132168       io_group_list(13) % to_be_processed = .FALSE.
    2114 
    2115        !incoming direct sw flux
     2169!
     2170!--    incoming direct sw flux
    21162171       io_group_list(14) = init_io_group(                                      &
    21172172          in_files = radiation_files,                                          &
     
    21212176       )
    21222177       io_group_list(14) % to_be_processed = .FALSE.
    2123 
    2124        !sw radiation balance
     2178!
     2179!--    sw radiation balance
    21252180       io_group_list(15) = init_io_group(                                      &
    21262181          in_files = radiation_files,                                          &
     
    21302185       )
    21312186       io_group_list(15) % to_be_processed = .FALSE.
    2132 
    2133        !lw radiation balance
     2187!
     2188!--    lw radiation balance
    21342189       io_group_list(16) = init_io_group(                                      &
    21352190          in_files = radiation_files,                                          &
     
    21432198
    21442199
     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!------------------------------------------------------------------------------!
    21452211    FUNCTION init_io_group(in_files, out_vars, in_var_list, kind,              &
    21462212                           n_output_quantities) RESULT(group)
     
    21572223       group % n_inputs = SIZE(in_var_list)
    21582224       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.
    21602232       IF ( PRESENT(n_output_quantities) )  THEN
    21612233          group % n_output_quantities = n_output_quantities
     
    22072279! Description:
    22082280! ------------
    2209 !> Initializes the the variable list.
     2281!> Initializes the variable list.
    22102282!------------------------------------------------------------------------------!
    22112283    SUBROUTINE setup_variable_tables(ic_mode)
     
    27872859          intermediate_grid = palm_intermediate                                &
    27882860       )
    2789        ! Radiation fluxes and balances
     2861
    27902862       output_var_table(38) = init_nc_var(                                     &
    27912863          name              = 'rad_swd_dif_0',                                 &
     
    31543226       output_var_table(64) % to_be_processed = .NOT. cfg % ug_is_set
    31553227
    3156        ! Attributes shared among all variables
     3228!
     3229!--    Attributes shared among all variables
    31573230       output_var_table(:) % source = nc_source_text
    31583231
     
    31643237! Description:
    31653238! ------------
    3166 !> Initializes and nc_var varible with the given parameters. The 'kind'
     3239!> Initializes an nc_var varible with the given parameters. The 'kind'
    31673240!> parameter is used to infer the correct netCDF IDs and the level of detail,
    31683241!> 'lod', as defined by the PALM-4U input data standard.
     
    31993272       SELECT CASE( TRIM(out_var_kind) )
    32003273
    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.
    32043278       CASE( 'init soil' )
    32053279          var % nt              = 1
     
    35913665       
    35923666       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
    35943669          ALLOCATE( preprocess_buffer( SIZE(input_buffer) ) )
    35953670
    3596           ! Allocate u and v arrays with scalar dimensions
     3671!
     3672!--       Allocate u and v arrays with scalar dimensions
    35973673          nx = SIZE(input_buffer(1) % array, 1)
    35983674          ny = SIZE(input_buffer(1) % array, 2)
     
    36033679 CALL run_control('time', 'alloc')
    36043680
    3605           ! interpolate U and V to centres
     3681!
     3682!--       interpolate U and V to centres
    36063683          CALL centre_velocities( u_face = input_buffer(1) % array,            &
    36073684                                  v_face = input_buffer(2) % array,            &
     
    36133690
    36143691          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
    36173695             DO k = 1, nz
    36183696             DO j = 1, ny
     
    36373715          END SELECT
    36383716
    3639           ! set values
    36403717          input_buffer(1) % array(1,:,:) = 0.0_dp
    36413718          input_buffer(2) % array(1,:,:) = 0.0_dp
     
    36583735          nz = SIZE(input_buffer(1) % array, 3)
    36593736
    3660           ! Compute absolute pressure if presure perturbation has been read in.
     3737!
     3738!--       Compute absolute pressure if presure perturbation has been read in.
    36613739          IF ( TRIM(group % in_var_list(2) % name) == 'PP' )  THEN
    36623740             message = "Absolute pressure, P, not available, " //              &
     
    36733751                                     RD, G, basic_state_pressure)
    36743752
     3753!
     3754!--             Overwrite pressure perturbation with absolute pressure. HECTO
     3755!--             converts pressure perturbation from hPa to Pa.
    36753756                input_buffer (2) % array(i,j,:) =                              &
    36763757                   HECTO * input_buffer (2) % array(i,j,:) +                   &
     
    36873768
    36883769          END IF
    3689           ! pressure
     3770!
     3771!--       mark pressure as preprocessed
    36903772          input_buffer(2) % is_preprocessed = .TRUE.
    36913773
    3692           ! Copy temperature to last input buffer array
     3774!
     3775!--       Copy temperature to the last input buffer array
    36933776          ALLOCATE(                                                            &
    36943777              input_buffer( group % n_output_quantities ) % array (nx, ny, nz) &
     
    36973780              input_buffer(1) % array(:,:,:)
    36983781
    3699           ! Convert absolute to potential temperature
     3782!
     3783!--       Convert absolute in place to potential temperature
    37003784          CALL potential_temperature(                                          &
    37013785             t = input_buffer(1) % array(:,:,:),                               &
     
    37063790          )
    37073791
    3708           ! potential temperature
     3792!
     3793!--       mark potential temperature as preprocessed
    37093794          input_buffer(1) % is_preprocessed = .TRUE.
    37103795
    3711           ! Convert temperature copy to density
     3796!
     3797!--       Convert temperature copy to density
    37123798          CALL moist_density(                                                  &
    37133799             t_rho = input_buffer(group % n_output_quantities) % array(:,:,:), &
     
    37183804          )
    37193805
    3720           ! qv
     3806!
     3807!--       mark qv as preprocessed
    37213808          input_buffer(3) % is_preprocessed = .TRUE.
    37223809
    3723           ! density
     3810!
     3811!--       mark density as preprocessed
    37243812          input_buffer(group % n_output_quantities) % is_preprocessed = .TRUE.
    37253813
     
    37803868          SELECT CASE(dt)
    37813869
     3870!
     3871!--       input has been accumulated over one hour. Leave as is
     3872!--       input_buffer(1) % array(:,:,:) carrries one-hour integral
    37823873          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
    37863879          CASE(2)
    3787              !input has been accumulated over two hours. Subtract previous step
    3788              !input_buffer(1) % array(:,:,:) carrries one-hour integral
    3789              !input_buffer(2) % array(:,:,:) carrries two-hour integral
    37903880             CALL deaverage(                                                   &
    37913881                      avg_1 = input_buffer(1) % array(:,:,:), t1 = 1.0_dp,     &
    37923882                      avg_2 = input_buffer(2) % array(:,:,:), t2 = 1.0_dp,     &
    37933883                      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
    37963891          CASE(3)
    3797              !input has been accumulated over three hours. Subtract previous step
    3798              !input_buffer(1) % array(:,:,:) carrries three-hour integral
    3799              !input_buffer(2) % array(:,:,:) still carrries two-hour integral
    38003892             CALL deaverage(                                                   &
    38013893                     avg_1 = input_buffer(2) % array(:,:,:), t1 = 1.0_dp,      &
    38023894                     avg_2 = input_buffer(1) % array(:,:,:), t2 = 1.0_dp,      &
    38033895                     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
    38053898
    38063899          CASE DEFAULT
     
    38183911
    38193912          hour = iter - 1
    3820           dt = MODULO(hour, 3) + 1 ! averaging period
     3913!
     3914!--       averaging period
     3915          dt = MODULO(hour, 3) + 1
    38213916          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
    38233920          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
    38273926          CASE(2)
    3828              !input has been accumulated over two hours. Subtract previous step
    3829              !input_buffer(1) % array(:,:,:) carrries one-hour integral
    3830              !input_buffer(2) % array(:,:,:) carrries two-hour integral
    38313927             CALL deaverage( input_buffer(1) % array(:,:,:), 1.0_dp,           &
    38323928                             input_buffer(2) % array(:,:,:), 2.0_dp,           &
    38333929                             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
    38363937          CASE(3)
    3837              !input has been accumulated over three hours. Subtract previous step
    3838              !input_buffer(1) % array(:,:,:) carrries three-hour integral
    3839              !input_buffer(2) % array(:,:,:) still carrries two-hour integral
    38403938             CALL deaverage( input_buffer(2) % array(:,:,:), 2.0_dp,           &
    38413939                             input_buffer(1) % array(:,:,:), 3.0_dp,           &
    38423940                             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
    38443943
    38453944          CASE DEFAULT
     
    38603959
    38613960
     3961!------------------------------------------------------------------------------!
    38623962! Description:
    38633963! ------------
     
    39284028                END DO
    39294029
    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
    39314032                IF (n_cells > 0)  THEN
    39324033                   array(i,j,:) = column(:) / n_cells
  • palm/trunk/UTIL/inifor/src/inifor_io.f90

    r3537 r3557  
    2626! -----------------
    2727! $Id$
     28! Updated documentation, removed unused subroutine write_netcdf_variable_2d()
     29!
     30!
     31! 3537 2018-11-20 10:53:14Z eckhard
    2832! New routine get_netcdf_dim_vector()
    2933!
     
    7377! Authors:
    7478! --------
    75 ! @author Eckhard Kadasch
     79!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
    7680!
    7781! Description:
     
    9296    IMPLICIT NONE
    9397
     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!------------------------------------------------------------------------------!
    94106    INTERFACE get_netcdf_variable
    95107        MODULE PROCEDURE get_netcdf_variable_int
     
    101113 CONTAINS
    102114
     115!------------------------------------------------------------------------------!
     116! Description:
     117! ------------
     118!> get_netcdf_variable_int() implements the integer variant for the
     119!> get_netcdf_variable interface.
     120!------------------------------------------------------------------------------!
    103121    SUBROUTINE get_netcdf_variable_int(in_file, in_var, buffer)
    104122
     
    139157
    140158
     159!------------------------------------------------------------------------------!
     160! Description:
     161! ------------
     162!> get_netcdf_variable_real() implements the real variant for the
     163!> get_netcdf_variable interface.
     164!------------------------------------------------------------------------------!
    141165    SUBROUTINE get_netcdf_variable_real(in_file, in_var, buffer)
    142166
     
    177201
    178202
    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)
    180210
    181211       CHARACTER(LEN=*), INTENT(IN)         ::  filename
    182        CHARACTER(LEN=*), INTENT(IN)         ::  varname
    183        REAL(dp), ALLOCATABLE, INTENT(INOUT) ::  array(:)
     212       CHARACTER(LEN=*), INTENT(IN)         ::  coordname
     213       REAL(dp), ALLOCATABLE, INTENT(INOUT) ::  coords(:)
    184214
    185215       INTEGER ::  ncid, varid, dimlen
     
    187217
    188218       IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
    189             nf90_inq_varid( ncid, varname, varid ) .EQ. NF90_NOERR )  THEN
     219            nf90_inq_varid( ncid, coordname, varid ) .EQ. NF90_NOERR )  THEN
    190220
    191221          CALL check(nf90_inquire_variable( ncid, varid, dimids = dimids ))
    192222          CALL check(nf90_inquire_dimension( ncid, dimids(1), len = dimlen ))
    193223
    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))
    196226
    197227       ELSE
    198228
    199           message = "Failed to read '" // TRIM(varname) // &
     229          message = "Failed to read '" // TRIM(coordname) // &
    200230             "' from file '" // TRIM(filename) // "'."
    201231          CALL abort('get_netcdf_dim_vector', message)
     
    206236
    207237
     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!------------------------------------------------------------------------------!
    208244    SUBROUTINE get_input_dimensions(in_var, ncid)
    209245
     
    232268
    233269
     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!------------------------------------------------------------------------------!
    234277    SUBROUTINE get_netcdf_start_and_count(in_var, start, count)
    235278
     
    251294       start = (/ 1, 1, 1 /)
    252295       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
    254298          in_var % dimlen(3) = in_var % dimlen(3) - 1
    255299
    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
    257302          start(3) = 2
    258303       END IF
     
    269314
    270315
     316!------------------------------------------------------------------------------!
     317! Description:
     318! ------------
     319!> Routine for defining netCDF variables in the dynamic driver, INIFOR's netCDF
     320!> output.
     321!------------------------------------------------------------------------------!
    271322    SUBROUTINE netcdf_define_variable(var, ncid)
    272323
     
    286337   
    287338
     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!------------------------------------------------------------------------------!
    288346    SUBROUTINE netcdf_get_dimensions(var, ncid)
    289347
     
    293351        CHARACTER(SNAME)            ::  null
    294352
    295         ! Remember dimension lenghts for NetCDF writing routine
    296353        DO i = 1, var % ndim
    297354           CALL check(nf90_inquire_dimension(ncid, var % dimids(i), &
     
    306363! Description:
    307364! ------------
    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.
    312367!------------------------------------------------------------------------------!
    313368    SUBROUTINE parse_command_line_arguments( cfg )
     
    330385       IF (arg_count .GT. 0)  THEN
    331386
    332           ! Every option should have an argument.
    333           !IF ( MOD(arg_count, 2) .NE. 0 )  THEN
    334           !   message = "Syntax error in command line."
    335           !   CALL abort('parse_command_line_arguments', message)
    336           !END IF
    337          
    338387          message = "The -clon and -clat command line options are depricated. " // &
    339388             "Please remove them form your inifor command and specify the " // &
     
    360409                cfg % debug = .TRUE.
    361410
    362              ! Elevation of the PALM-4U domain above sea level
    363411             CASE( '-z0', '-z', '--elevation' )
    364412                CALL get_option_argument( i, arg )
    365413                READ(arg, *) cfg % z0
    366414
    367              ! surface pressure, at z0
    368415             CASE( '-p0', '-r', '--surface-pressure' )
    369416                cfg % p0_is_set = .TRUE.
     
    371418                READ(arg, *) cfg % p0
    372419
    373              ! geostrophic wind in x direction
    374420             CASE( '-ug', '-u', '--geostrophic-u' )
    375421                cfg % ug_is_set = .TRUE.
     
    377423                READ(arg, *) cfg % ug
    378424
    379              ! geostrophic wind in y direction
    380425             CASE( '-vg', '-v', '--geostrophic-v' )
    381426                cfg % vg_is_set = .TRUE.
     
    383428                READ(arg, *) cfg % vg
    384429
    385              ! domain centre geographical longitude and latitude
    386430             CASE( '-clon', '-clat' )
    387431                CALL abort('parse_command_line_arguments', message)         
    388                 !READ(arg, *) lambda_cg
    389                 !lambda_cg = lambda_cg * TO_RADIANS
    390                 !READ(arg, *) phi_cg
    391                 !phi_cg = phi_cg * TO_RADIANS
    392432
    393433             CASE( '-path', '-p', '--path' )
     
    444484                cfg % namelist_file = TRIM(arg)
    445485
    446              ! initial condition mode: 'profile' / 'volume'
    447486             CASE( '-mode', '-i', '--init-mode' )
    448487                CALL get_option_argument( i, arg )
    449488                cfg % ic_mode = TRIM(arg)
    450489
    451              ! boundary conditions / forcing mode: 'ideal' / 'real'
    452490             CASE( '-f', '--forcing-mode' )
    453491                CALL get_option_argument( i, arg )
     
    484522
    485523   
     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!------------------------------------------------------------------------------!
    486530   SUBROUTINE get_option_argument(i, arg)
    487531      CHARACTER(LEN=PATH), INTENT(INOUT) ::  arg
     
    494538
    495539
     540!------------------------------------------------------------------------------!
     541! Description:
     542! ------------
     543!> Checks the INIFOR configuration 'cfg' for plausibility.
     544!------------------------------------------------------------------------------!
    496545   SUBROUTINE validate_config(cfg)
    497546      TYPE(inifor_config), INTENT(IN) ::  cfg
     
    503552      all_files_present = all_files_present .AND. file_present(cfg % soiltyp_file, 'SOILTYP')
    504553
    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.
    506556      IF (TRIM(cfg % static_driver_file) .NE. '')  THEN
    507557         all_files_present = all_files_present .AND. file_present(cfg % static_driver_file, 'static driver')
     
    557607
    558608
     609!------------------------------------------------------------------------------!
     610! Description:
     611! ------------
     612!> Check whether the given file is present on the filesystem.
     613!------------------------------------------------------------------------------!
    559614   LOGICAL FUNCTION file_present(filename, kind)
    560615      CHARACTER(LEN=PATH), INTENT(IN) ::  filename
     
    585640! Description:
    586641! ------------
    587 !> This routine initializes the Inifor output file, i.e. the PALM-4U
    588 !> initializing and forcing data as a NetCDF file.
     642!> This routine initializes the dynamic driver file, i.e. INIFOR's netCDF output
     643!> file.
    589644!>
    590645!> Besides writing metadata, such as global attributes, coordinates, variables,
    591 !> in the NetCDF file, various NetCDF IDs are saved for later, when Inifor
     646!> in the netCDF file, various netCDF IDs are saved for later, when INIFOR
    592647!> writes the actual data.
    593648!------------------------------------------------------------------------------!
     
    611666       CALL report('setup_netcdf_dimensions', message)
    612667
    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
    614670#if defined( __netcdf4 )
    615671       CALL check(nf90_create(TRIM(output_file % name), OR(NF90_CLOBBER, NF90_HDF5), ncid))
     
    646702       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat',     TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)'))))
    647703       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)
    651708       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_z',       TRIM(real_to_str(z0, '(F18.13)'))))
    652709       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION)))
     
    658715!- Section 2a: Define dimensions for cell centers (scalars in soil and atmosph.)
    659716!------------------------------------------------------------------------------
    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
    681747       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
    685754       CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file % dimids_soil(3), dimvarids(3)))
    686755       CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "depth_below_land"))
    687756       CALL check(nf90_put_att(ncid, dimvarids(3), "positive", "down"))
    688757       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
    690761!
    691762!------------------------------------------------------------------------------
    692763!- Section 2b: Define dimensions for cell faces/velocities
    693764!------------------------------------------------------------------------------
    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
    713792
    714793!
     
    739818       CALL check(nf90_put_var(ncid, output_file % dimvarids_vel(3), palm_grid%zw))
    740819
    741        ! TODO Read in soil depths from input file before this.
     820!
     821!--    TODO Read in soil depths from input file before this.
    742822       CALL check(nf90_put_var(ncid, output_file % dimvarids_soil(3), palm_grid%depths))
    743823
    744        ! Write time vector
     824!
     825!--    Write time vector
    745826       CALL check(nf90_put_var(ncid, output_file % dimvarid_time, output_file % time))
    746827
    747        ! Close the file
     828!
     829!--    Close the file
    748830       CALL check(nf90_close(ncid))
    749831
     
    751833
    752834
     835!------------------------------------------------------------------------------!
     836! Description:
     837! ------------
     838!> Defines the netCDF variables to be written to the dynamic driver file
     839!------------------------------------------------------------------------------!
    753840    SUBROUTINE setup_netcdf_variables(filename, output_variable_table, debug)
    754841
     
    824911!- Section 1: Load input buffers for accumulated variables
    825912!------------------------------------------------------------------------------
     913!
     914!--    radiation budgets, precipitation
    826915       IF (group % kind == 'running average' .OR.                              &
    827            group % kind == 'accumulated')  THEN ! radiation budgets, precipitation
     916           group % kind == 'accumulated')  THEN
    828917
    829918          IF (SIZE(group % in_var_list) .GT. 1 ) THEN
     
    835924          END IF
    836925
    837           ! use two buffer arrays
     926!
     927!--       use two buffer arrays
    838928          nbuffers = 2
    839929          IF ( .NOT. ALLOCATED( buffer ) )  ALLOCATE( buffer(nbuffers) )
    840930
    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
    843936          buf_id = select_buffer(hour)
    844937
     
    859952       ELSE
    860953
    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)
    867961          nbuffers = MAX( group % n_inputs, group % n_output_quantities )
    868962          ALLOCATE( buffer(nbuffers) )
    869963 CALL run_control('time', 'alloc')
    870964         
    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.
    872967          DO ivar = 1, group % n_inputs
    873968
    874969             input_var => group % in_var_list(ivar)
    875970
    876              ! Check wheather P or PP is present in input file
     971!
     972!            Check wheather P or PP is present in input file
    877973             IF (input_var % name == 'P')  THEN
    878974                input_var % name = TRIM( get_pressure_varname(input_file) )
     
    891987
    892988
     989!------------------------------------------------------------------------------!
     990! Description:
     991! ------------
     992!> Select the appropriate buffer ID for accumulated COSMO input variables
     993!> depending on the current hour.
     994!------------------------------------------------------------------------------!
    893995    INTEGER FUNCTION select_buffer(hour)
    894996       INTEGER, INTENT(IN) ::  hour
     
    9431045
    9441046
     1047!------------------------------------------------------------------------------!
     1048! Description:
     1049! ------------
     1050!> Read the given global attribute form the given netCDF file.
     1051!------------------------------------------------------------------------------!
    9451052    FUNCTION get_netcdf_attribute(filename, attribute) RESULT(attribute_value)
    9461053
     
    9481055       REAL(dp)                     ::  attribute_value
    9491056
    950        INTEGER                      :: ncid
     1057       INTEGER                      ::  ncid
    9511058
    9521059       IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR )  THEN
     
    9661073
    9671074
     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!------------------------------------------------------------------------------!
    9681082    SUBROUTINE update_output(var, array, iter, output_file, cfg)
    9691083       TYPE(nc_var), INTENT(IN)  ::  var
     
    9801094       )
    9811095
    982        ! Skip time dimension for output
     1096!
     1097!--    Skip time dimension for output
    9831098       ndim = var % ndim
    9841099       IF ( var_is_time_dependent )  ndim = var % ndim - 1
     
    9901105       CALL check(nf90_open(output_file % name, NF90_WRITE, ncid))
    9911106
    992        ! Reduce dimension of output array according to variable kind
     1107!
     1108!--    Reduce dimension of output array according to variable kind
    9931109       SELECT CASE (TRIM(var % kind))
    9941110       
     
    10861202
    10871203
    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!------------------------------------------------------------------------------!
    11301209    SUBROUTINE check(status)
    11311210
  • palm/trunk/UTIL/inifor/src/inifor_transform.f90

    r3537 r3557  
    2626! -----------------
    2727! $Id$
     28! Updated documentation
     29!
     30!
     31! 3537 2018-11-20 10:53:14Z eckhard
    2832! bugfix: working precision added
    2933!
     
    5256! Authors:
    5357! --------
    54 ! @author Eckhard Kadasch
     58!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
    5559!
    5660! Description:
     
    110114       DO k = nz, LBOUND(out_arr, 3), -1
    111115
    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.
    116121          IF (outgrid % w_verti(i,j,k,1) < -1.0_dp .AND. k < nz)  THEN
    117122             out_arr(i,j,k) = out_arr(i,j,k+1)
     
    155160!------------------------------------------------------------------------------!
    156161    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.
    159165       TYPE(grid_definition), INTENT(IN)  ::  outgrid
    160166       REAL(dp), INTENT(IN)               ::  invar(0:,0:,0:)
     
    164170       INTEGER ::  i, j, k, l
    165171
    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
    168175           message = "Output array for '" // TRIM(ncvar % name) // "' has ' more levels (" // &
    169176              TRIM(str(UBOUND(outvar, 3))) // ") than input variable ("//&
     
    190197
    191198
     199!------------------------------------------------------------------------------!
     200! Description:
     201! ------------
     202!> Compute the horizontal average of the in_arr(:,:,:) and store it in
     203!> out_arr(:)
     204!------------------------------------------------------------------------------!
    192205    SUBROUTINE average_2d(in_arr, out_arr, ii, jj)
    193206       REAL(dp), INTENT(IN)              ::  in_arr(0:,0:,0:)
     
    200213       IF (SIZE(ii) .NE. SIZE(jj))  THEN
    201214          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, " //   &
    203216             NEW_LINE(' ') // "jj has " // str(SIZE(jj)) // "."
    204217          CALL abort('average_2d', message)
     
    211224             i = ii(l)
    212225             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)
    215227          END DO
    216228
     
    223235
    224236
     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!------------------------------------------------------------------------------!
    225248    SUBROUTINE interpolate_3d(source_array, palm_array, palm_intermediate, palm_grid)
    226249       TYPE(grid_definition), INTENT(IN) ::  palm_intermediate, palm_grid
     
    232255       nx = palm_intermediate % nx
    233256       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.
    239263       ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nlev-1)) !
    240264
    241265       CALL interpolate_2d(source_array, intermediate_array, palm_intermediate)
    242266
    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.
    245270       CALL interpolate_1d(intermediate_array, palm_array, palm_grid)
    246271
     
    250275
    251276
     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!------------------------------------------------------------------------------!
    252283    SUBROUTINE average_profile(source_array, profile_array, avg_grid)
    253284       TYPE(grid_definition), INTENT(IN)          ::  avg_grid
     
    265296          j_source = avg_grid % jjj(l)
    266297
    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
    270305
    271306                k_source = avg_grid % kkk(l, k_profile, m)
     
    274309                   + avg_grid % w(l, k_profile, m)                             &
    275310                   * 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
    282322
    283323       ni_columns = 1.0_dp / avg_grid % n_columns
    284324       profile_array(:) = profile_array(:) * ni_columns
    285325
    286        ! Extrapolate constant to the bottom
     326!
     327!--    Constant extrapolation to the bottom
    287328       profile_array(1:avg_grid % k_min-1) = profile_array(avg_grid % k_min)
    288329
     
    290331
    291332
     333!------------------------------------------------------------------------------!
     334! Description:
     335! ------------
     336!> Extrapolates density linearly from the level 'k_min' downwards.
     337!------------------------------------------------------------------------------!
    292338    SUBROUTINE extrapolate_density(rho, avg_grid)
    293339       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  rho
     
    308354
    309355
     356!------------------------------------------------------------------------------!
     357! Description:
     358! ------------
     359!> Driver for extrapolating pressure from PALM level k_min downwards
     360!------------------------------------------------------------------------------!
    310361    SUBROUTINE extrapolate_pressure(p, rho, avg_grid)
    311362       REAL(dp), DIMENSION(:), INTENT(IN)    ::  rho
     
    405456       REAL(dp) :: ri
    406457
    407        ! TODO check dimensions of lat/lon and y/x match
     458!
     459!--    TODO check dimensions of lat/lon and y/x match
    408460
    409461       ri = 1.0_dp / r
     
    438490       REAL(dp) :: ri
    439491
    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.
    443496       ri = 1.0_dp / r
    444497
     
    448501
    449502
     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!------------------------------------------------------------------------------!
    450509    REAL(dp) FUNCTION phic_to_phin(phi_c)
    451510        REAL(dp), INTENT(IN) ::  phi_c
     
    456515
    457516   
     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!------------------------------------------------------------------------------!
    458524    REAL(dp) FUNCTION lamc_to_lamn(phi_c, lam_c)
    459525       REAL(dp), INTENT(IN) ::  phi_c, lam_c
     
    467533
    468534
     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!------------------------------------------------------------------------------!
    469543    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
    476552           gamma_from_hemisphere = PI
    477553       ELSE
     
    550626       
    551627
     628!------------------------------------------------------------------------------!
     629! Description:
     630! ------------
     631!> Rotate the given vector field (x(:), y(:)) by the given 'angle'.
     632!------------------------------------------------------------------------------!
    552633    SUBROUTINE rotate_vector_field(x, y, angle)
    553634       REAL(dp), DIMENSION(:), INTENT(INOUT) :: x, y  !< x and y coodrinate in arbitrary units
     
    559640       sine = SIN(angle * TO_RADIANS)
    560641       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 ]
    565647       rotation = RESHAPE( (/cosine, sine, -sine, cosine/), (/2, 2/) )
    566648
     
    588670!>    Datenbanken des DWD.
    589671!>    https://www.dwd.de/SharedDocs/downloads/DE/modelldokumentationen/nwv/cosmo_d2/cosmo_d2_dbbeschr_aktuell.pdf?__blob=publicationFile&v=2
    590 !>
     672!------------------------------------------------------------------------------!
    591673    FUNCTION meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g)          &
    592674       RESULT(delta)
     
    664746       DO j = 0, UBOUND(palm_clon, 2)!palm_grid % ny
    665747       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.
    668751          lonpos = (palm_clon(i,j) - lon0) * cosmo_dxi
    669752          latpos = (palm_clat(i,j) - lat0) * cosmo_dyi
     
    693776
    694777   
     778!------------------------------------------------------------------------------!
     779! Description:
     780! ------------
     781!> Computes linear vertical interpolation neighbour indices and weights for each
     782!> column of the given palm grid.
     783!------------------------------------------------------------------------------!
    695784    SUBROUTINE find_vertical_neighbours_and_weights_interp( palm_grid,         &
    696785                                                            palm_intermediate )
     
    709798       nlev = palm_intermediate % nz
    710799
    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
    712802       DO j = 0, ny
    713803       DO i = 0, nx
     
    718808          column_top  = palm_intermediate % h(i,j,nlev)
    719809
    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.
    724815          DO k = 1, nz
    725816
    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
    728820             current_height = palm_grid % z(k) + palm_grid % z0
    729821             h_top    = palm_intermediate % h(i,j,k_intermediate+1)
     
    738830             )
    739831
    740              ! set default weights
     832!
     833!--          set default weights
    741834             palm_grid % w_verti(i,j,k,1:2) = 0.0_dp
    742835
     
    755848
    756849             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
    759853                DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
    760854                   k_intermediate = k_intermediate + 1
     
    777871                palm_grid % kk(i,j,k,2) = k_intermediate + 1
    778872
    779                 ! copmute vertical weights
     873!
     874!--             compute vertical weights
    780875                weight = (h_top - current_height) / (h_top - h_bottom)
    781876                palm_grid % w_verti(i,j,k,1) = weight
     
    791886
    792887
     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!------------------------------------------------------------------------------!
    793900    SUBROUTINE find_vertical_neighbours_and_weights_average( avg_grid )
    794901       TYPE(grid_definition), INTENT(INOUT) ::  avg_grid
     
    805912       nlev = SIZE(avg_grid % cosmo_h, 3)
    806913
    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
    808916       DO l = 1, avg_grid % n_columns
    809917
     
    814922          column_top  = avg_grid % cosmo_h(i,j,nlev)
    815923
    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.
    820929          k_intermediate = 1 !avg_grid % cosmo_h is indezed 1-based.
    821930          DO k_palm = 1, avg_grid % nz
    822931
    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
    825935             current_height = avg_grid % z(k_palm) + avg_grid % z0
    826936             h_top    = avg_grid % cosmo_h(i,j,k_intermediate+1)
    827937             h_bottom = avg_grid % cosmo_h(i,j,k_intermediate)
    828938
    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)
    830943             point_is_below_grid = (current_height < column_base)
    831944
     
    835948             )
    836949
    837              ! set default weights
     950!
     951!--          set default weights
    838952             avg_grid % w(l,k_palm,1:2) = 0.0_dp
    839953
     
    852966                avg_grid % k_min = MAX(k_palm + 1, avg_grid % k_min)
    853967             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
    856971                DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
    857972                   k_intermediate = k_intermediate + 1
     
    865980                END DO
    866981
    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.
    869985                IF (k_intermediate > nlev-1)  THEN
    870986                   message = "Index " // TRIM(str(k_intermediate)) //          &
     
    876992                avg_grid % kkk(l,k_palm,2) = k_intermediate + 1
    877993
    878                 ! copmute vertical weights
     994!
     995!--             compute vertical weights
    879996                weight = (h_top - current_height) / (h_top - h_bottom)
    880997                avg_grid % w(l,k_palm,1) = weight
     
    882999             END IF
    8831000
    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
    8861008 
    8871009    END SUBROUTINE find_vertical_neighbours_and_weights_average
     
    9311053!                         ii(i,j,1/2)        ii(i,j,3/4)
    9321054!         
     1055!------------------------------------------------------------------------------!
    9331056    SUBROUTINE compute_horizontal_interp_weights(cosmo_lat, cosmo_lon,         &
    9341057       palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)
     
    9501073       DO i = 0, UBOUND(palm_clon, 1)
    9511074     
    952           ! weight in lambda direction
     1075!
     1076!--       weight in lambda direction
    9531077          wl = ( cosmo_lon(palm_ii(i,j,4)) - palm_clon(i,j) ) * cosmo_dxi
    9541078
    955           ! weight in phi direction
     1079!
     1080!--       weight in phi direction
    9561081          wp = ( cosmo_lat(palm_jj(i,j,2)) - palm_clat(i,j) ) * cosmo_dyi
    9571082
     
    9881113!> COSMO-DE the velocity points are staggared one half grid spaceing up-grid
    9891114!> which means the first centre point has to be omitted and is set to zero.
     1115!------------------------------------------------------------------------------!
    9901116    SUBROUTINE centre_velocities(u_face, v_face, u_centre, v_centre)
    9911117       REAL(dp), DIMENSION(0:,0:,0:), INTENT(IN)  ::  u_face, v_face
     
    10041130
    10051131
     1132!------------------------------------------------------------------------------!
     1133! Description:
     1134! ------------
     1135!> Compute the geographical latitude of a point given in rotated-pole cordinates
     1136!------------------------------------------------------------------------------!
    10061137    FUNCTION phirot2phi (phirot, rlarot, polphi, pollam, polgam)
    10071138   
     
    10411172
    10421173
     1174!------------------------------------------------------------------------------!
     1175! Description:
     1176! ------------
     1177!> Compute the geographical latitude of a point given in rotated-pole cordinates
     1178!------------------------------------------------------------------------------!
    10431179    FUNCTION phi2phirot (phi, rla, polphi, pollam)
    10441180   
     
    10721208
    10731209
     1210!------------------------------------------------------------------------------!
     1211! Description:
     1212! ------------
     1213!> Compute the geographical longitude of a point given in rotated-pole cordinates
     1214!------------------------------------------------------------------------------!
    10741215    FUNCTION rlarot2rla(phirot, rlarot, polphi, pollam, polgam)
    10751216   
     
    11231264
    11241265
     1266!------------------------------------------------------------------------------!
     1267! Description:
     1268! ------------
     1269!> Compute the rotated-pole longitude of a point given in geographical cordinates
     1270!------------------------------------------------------------------------------!
    11251271    FUNCTION rla2rlarot ( phi, rla, polphi, pollam, polgam )
    11261272
     
    11621308
    11631309
     1310!------------------------------------------------------------------------------!
     1311! Description:
     1312! ------------
     1313!> Rotate the given velocity vector (u,v) from the geographical to the
     1314!> rotated-pole system
     1315!------------------------------------------------------------------------------!
    11641316    SUBROUTINE uv2uvrot(u, v, rlat, rlon, pollat, pollon, urot, vrot)
    11651317   
     
    11871339
    11881340
     1341!------------------------------------------------------------------------------!
     1342! Description:
     1343! ------------
     1344!> Rotate the given velocity vector (urot, vrot) from the rotated-pole to the
     1345!> geographical system
     1346!------------------------------------------------------------------------------!
    11891347    SUBROUTINE uvrot2uv (urot, vrot, rlat, rlon, pollat, pollon, u, v)
    11901348   
  • palm/trunk/UTIL/inifor/src/inifor_types.f90

    r3447 r3557  
    2626! -----------------
    2727! $Id$
     28! Updated documentation
     29!
     30!
     31! 3447 2018-10-29 15:52:54Z eckhard
    2832! Renamed source files for compatibilty with PALM build system
    2933!
     
    5155! Authors:
    5256! --------
    53 ! @author Eckhard Kadasch
     57!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
    5458!
    5559! Description:
     
    6670 IMPLICIT NONE
    6771
     72!------------------------------------------------------------------------------!
     73! Description:
     74! ------------
     75!> Contaner for the INIFOR command-line configuration
     76!------------------------------------------------------------------------------!
    6877 TYPE inifor_config
    6978    CHARACTER(LEN=DATE)  ::  start_date           !< String of the FORMAT YYYYMMDDHH indicating the start of the intended PALM-4U simulation
     
    8291    CHARACTER(LEN=SNAME) ::  soilmoisture_prefix  !< Prefix of input files for soil moisture spin-up, e.g 'laf' for COSMO-DE analyses
    8392
    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]
    93102    REAL(dp)             ::  averaging_angle      !< latitudal and longitudal width of averaging regions [deg]
    94103   
    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
    105113 END TYPE inifor_config
    106114
     115
     116!------------------------------------------------------------------------------!
     117! Description:
     118! ------------
     119!> Container for grid data, in partucular coordinates, interpolation neighbours
     120!> and weights
     121!------------------------------------------------------------------------------!
    107122 TYPE grid_definition
    108123    CHARACTER(LEN=SNAME)  ::  name(3)       !< names of the grid dimensions, e.g. (/'x', 'y', 'z'/) or (/'latitude', 'longitude', 'height'/)
     
    152167
    153168
     169!------------------------------------------------------------------------------!
     170! Description:
     171! ------------
     172!> Container for name and dimensions of the netCDF output file
     173!------------------------------------------------------------------------------!
    154174 TYPE nc_file
    155     CHARACTER(LEN=PATH)   ::  name          !< file name
    156     INTEGER               ::  dimid_time    !< NetCDF IDs of the time dimension
    157     INTEGER               ::  dimids_scl(3) !< NetCDF IDs of the grid dimensions for scalar points x, y, z
    158     INTEGER               ::  dimids_vel(3) !< NetCDF IDs of the grid dimensions for velocity points xu, yu, zu
    159     INTEGER               ::  dimids_soil(3)!< NetCDF IDs of the grid dimensions for soil points x, y, depth
    160     INTEGER               ::  dimvarid_time !< NetCDF IDs of the time variable
    161     INTEGER               ::  dimvarids_scl(3) !< NetCDF IDs of the grid coordinates of scalars x, y, z
    162     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, depth
    164     REAL(dp), POINTER     ::  time(:)       ! vector of output time steps
     175    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
    165185 END TYPE nc_file
    166186
    167187
     188!------------------------------------------------------------------------------!
     189! Description:
     190! ------------
     191!> Metadata container for netCDF variables
     192!------------------------------------------------------------------------------!
    168193 TYPE nc_var
    169194    INTEGER                               ::  varid     !< NetCDF ID of the variable
     
    193218
    194219
    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
    196229    INTEGER                          ::  nt             !< maximum number of output time steps across all output variables
    197230    INTEGER                          ::  nv             !< number of netCDF output variables
     
    208241
    209242
     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!------------------------------------------------------------------------------!
    210250 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
    213253 END TYPE container
    214254
  • palm/trunk/UTIL/inifor/src/inifor_util.f90

    r3447 r3557  
    2626! -----------------
    2727! $Id$
     28! Updated documentation
     29!
     30!
     31! 3447 2018-10-29 15:52:54Z eckhard
    2832! Renamed source files for compatibilty with PALM build system
    2933!
     
    4448! Authors:
    4549! --------
    46 ! @author Eckhard Kadasch
     50!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
    4751!
    4852! Description:
     
    6165    IMPLICIT NONE
    6266
     67!------------------------------------------------------------------------------!
     68! Description:
     69! ------------
     70!> Fortran implementation of C's struct tm for representing points in time
     71!------------------------------------------------------------------------------!
    6372    TYPE, BIND(c) :: tm_struct
    6473       INTEGER(C_INT) :: tm_sec     !< seconds after the minute [0, 61]
     
    7584    INTERFACE
    7685
     86!------------------------------------------------------------------------------!
     87! Description:
     88! ------------
     89!> Interface to C's strptime function, which converts a string to a tm
     90!> structure.
     91!------------------------------------------------------------------------------!
    7792       FUNCTION strptime(string, format, timeinfo) BIND(c, NAME='strptime')
    7893          IMPORT :: C_CHAR, C_SIZE_T, tm_struct
     
    87102
    88103
     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!------------------------------------------------------------------------------!
    89110       FUNCTION strftime(string, string_len, format, timeinfo) BIND(c, NAME='strftime')
    90111          IMPORT :: C_CHAR, C_SIZE_T, tm_struct
     
    101122
    102123
     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!------------------------------------------------------------------------------!
    103132       FUNCTION mktime(timeinfo) BIND(c, NAME='mktime')
    104133          IMPORT :: C_PTR, tm_struct
     
    115144 CONTAINS
    116145
     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!------------------------------------------------------------------------------!
    117153    CHARACTER(LEN=DATE) FUNCTION add_hours_to(date_string, hours)
    118154       CHARACTER(LEN=DATE), INTENT(IN)          ::  date_string
     
    143179
    144180
     181!------------------------------------------------------------------------------!
     182! Description:
     183! ------------
     184!> Print all members of the given tm structure
     185!------------------------------------------------------------------------------!
    145186    SUBROUTINE print_tm(timeinfo)
    146187       TYPE(tm_struct), INTENT(IN) :: timeinfo
     
    158199
    159200   
     201!------------------------------------------------------------------------------!
     202! Description:
     203! ------------
     204!> Initialize the given tm structure with zero values
     205!------------------------------------------------------------------------------!
    160206    SUBROUTINE init_tm(timeinfo)
    161207       TYPE(tm_struct), INTENT(INOUT) :: timeinfo
     
    177223
    178224
    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!------------------------------------------------------------------------------!
    219231    SUBROUTINE linspace(start, stop, array)
    220232
     
    255267
    256268
     269!------------------------------------------------------------------------------!
     270! Description:
     271! ------------
     272!>
     273!------------------------------------------------------------------------------!
    257274    SUBROUTINE deaverage(avg_1, t1, avg_2, t2, avg_3, t3)
    258275
     
    270287
    271288
     289!------------------------------------------------------------------------------!
     290! Description:
     291! ------------
     292!> Compute the COSMO-DE/-D2 basic state pressure profile
     293!------------------------------------------------------------------------------!
    272294    SUBROUTINE get_basic_state(z, beta, p_sl, t_sl, rd, g, p0)
    273295
     
    278300       REAL(dp), INTENT(IN)  ::  rd     !< ideal gas constant of dry air [J/kg/K]
    279301       REAL(dp), INTENT(IN)  ::  g      !< acceleration of Earth's gravity [m/s^2]
    280        REAL(dp), INTENT(OUT) ::  p0(1:) !< COSMO-DE basic state pressure [Pa]
     302       REAL(dp), INTENT(OUT) ::  p0(1:) !< COSMO basic state pressure [Pa]
    281303       REAL(dp) ::  root_frac, factor   !< precomputed factors
    282304
     
    285307
    286308       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       )
    289311
    290312    END SUBROUTINE get_basic_state
     
    345367
    346368
     369!------------------------------------------------------------------------------!
     370! Description:
     371! ------------
     372!> Converts the given real value to a string
     373!------------------------------------------------------------------------------!
    347374    CHARACTER(LEN=16) FUNCTION real_to_str_f(val)
    348375
     
    355382
    356383
     384!------------------------------------------------------------------------------!
     385! Description:
     386! ------------
     387!> Converts the given integer value to a string
     388!------------------------------------------------------------------------------!
    357389    CHARACTER(LEN=10) FUNCTION str(val)
    358390
     
    365397
    366398
    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!------------------------------------------------------------------------------!
    380404    SUBROUTINE normalize_path(path)
    381405       
Note: See TracChangeset for help on using the changeset viewer.