Index: /palm/trunk/UTIL/inifor/README
===================================================================
--- /palm/trunk/UTIL/inifor/README	(revision 3800)
+++ /palm/trunk/UTIL/inifor/README	(revision 3801)
@@ -1,3 +1,3 @@
-# INIFOR - Mesoscale Interface for Initializing and Forcing PALM-4U (v1.4.7)
+# INIFOR - Mesoscale Interface for Initializing and Forcing PALM-4U
 
 INIFOR provides the meteorological fields required to initialize and drive the
Index: /palm/trunk/UTIL/inifor/src/inifor_defs.f90
===================================================================
--- /palm/trunk/UTIL/inifor/src/inifor_defs.f90	(revision 3800)
+++ /palm/trunk/UTIL/inifor/src/inifor_defs.f90	(revision 3801)
@@ -26,4 +26,9 @@
 ! -----------------
 ! $Id$
+! Defined netCDF variable names for COSMO grid
+! Bumped version number
+! 
+! 
+! 3779 2019-03-05 11:13:35Z eckhard
 ! Updated version number to 1.4.7, updated copyright note
 ! 
@@ -111,5 +116,5 @@
 
 ! 
-!-- COSMO-DE parameters
+!-- COSMO parameters
  INTEGER, PARAMETER  ::  WATER_ID = 9                !< Integer corresponding to the water soil type in COSMO-DE [-]
  REAL(dp), PARAMETER ::  EARTH_RADIUS = 6371229.0_dp !< Earth radius used in COSMO-DE [m]
@@ -134,4 +139,14 @@
  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]
 
+! 
+!-- COSMO netCDF names
+CHARACTER(SNAME) ::  NC_DEPTH_NAME = 'depth_2'
+CHARACTER(SNAME) ::  NC_HHL_NAME = 'HHL'
+CHARACTER(SNAME) ::  NC_RLAT_NAME = 'rlat'
+CHARACTER(SNAME) ::  NC_RLON_NAME = 'rlon'
+CHARACTER(SNAME) ::  NC_ROTATED_POLE_NAME = 'rotated_pole'
+CHARACTER(SNAME) ::  NC_POLE_LATITUDE_NAME = 'grid_north_pole_latitude'
+CHARACTER(SNAME) ::  NC_POLE_LONGITUDE_NAME = 'grid_north_pole_longitude'
+
 !
 !-- INIFOR parameters
@@ -140,5 +155,5 @@
  INTEGER, PARAMETER          ::  FORCING_STEP = 1             !< Number of hours between forcing time steps [h]
  REAL(dp), PARAMETER         ::  NUDGING_TAU = 21600.0_dp     !< Nudging relaxation time scale [s]
- CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.4.7'            !< INIFOR version number
+ CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.4.8'            !< INIFOR version number
  CHARACTER(LEN=*), PARAMETER ::  COPYRIGHT = 'Copyright 2017-2019 Leibniz Universitaet Hannover' // &
      ACHAR( 10 ) // ' Copyright 2017-2019 Deutscher Wetterdienst Offenbach' !< Copyright notice
Index: /palm/trunk/UTIL/inifor/src/inifor_grid.f90
===================================================================
--- /palm/trunk/UTIL/inifor/src/inifor_grid.f90	(revision 3800)
+++ /palm/trunk/UTIL/inifor/src/inifor_grid.f90	(revision 3801)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Read COSMO rotated pole from HHL file
+! Check for PALM and averaging domains extending outside COSMO domain
+! 
+! 3785 2019-03-06 10:41:14Z eckhard
 ! Assigned names to averaging grids
 ! Improved variable naming and minor clean-up
@@ -129,5 +133,5 @@
                RHO_L, OMEGA, HECTO
     USE inifor_io,                                                             &
-        ONLY:  get_netcdf_attribute, get_netcdf_dim_vector,                    &
+        ONLY:  get_cosmo_grid, get_netcdf_attribute, get_netcdf_dim_vector,    &
                get_netcdf_variable, parse_command_line_arguments,              &
                get_input_file_list, validate_config
@@ -195,4 +199,8 @@
     REAL(dp) ::  latmin_palm       = 0.0_dp       !< Minimunm latitude of PALM grid [COSMO rotated-pole rad]
     REAL(dp) ::  latmax_palm       = 0.0_dp       !< Maximum latitude of PALM grid [COSMO rotated-pole rad]
+    REAL(dp) ::  lonmin_tot        = 0.0_dp       !< Minimunm longitude of required COSMO data [COSMO rotated-pole rad]
+    REAL(dp) ::  lonmax_tot        = 0.0_dp       !< Maximum longitude of required COSMO data [COSMO rotated-pole rad]
+    REAL(dp) ::  latmin_tot        = 0.0_dp       !< Minimunm latitude of required COSMO data [COSMO rotated-pole rad]
+    REAL(dp) ::  latmax_tot        = 0.0_dp       !< Maximum latitude of required COSMO data [COSMO rotated-pole rad]
     REAL(dp) ::  latitude          = 0.0_dp       !< geographical latitude of the PALM-4U origin, from inipar namelist [deg]
     REAL(dp) ::  longitude         = 0.0_dp       !< geographical longitude of the PALM-4U origin, from inipar namelist [deg]
@@ -242,8 +250,9 @@
     LOGICAL ::  ls_forcing_variables_required !< flag controlling whether large-scale forcing variables are to be computed
     LOGICAL ::  surface_forcing_required      !< flag controlling whether surface forcing variables are to be computed
+    LOGICAL ::  palm_domain_outside_cosmo     !< indicates whether COSMO grid covers the PALM domain and the geostrophic averaging domains
 
     TYPE(nc_var), ALLOCATABLE, TARGET ::  input_var_table(:)  !< table of input variables
     TYPE(nc_var), ALLOCATABLE, TARGET ::  output_var_table(:) !< table of input variables
-    TYPE(nc_var)                      ::  cosmo_var           !< COSMO dummy variable, used for reading HHL, rlon, rlat
+    TYPE(nc_var) ::  cosmo_var                                !< COSMO dummy variable, used for reading HHL, rlon, rlat
 
     TYPE(grid_definition), TARGET ::  palm_grid                       !< PALM-4U grid in the target system (COSMO-DE rotated-pole)
@@ -356,10 +365,4 @@
 
 !
-!--    COSMO-DE and -D2 rotated pole position
-       phi_n     =   40.0_dp * TO_RADIANS
-       phi_equat =   50.0_dp * TO_RADIANS
-       lambda_n  = -170.0_dp * TO_RADIANS
-
-!
 !--    Defaultmain centre (_c) of the PALM-4U grid in the geographical system (_g)
        origin_lat = 52.325079_dp * TO_RADIANS ! south-west of Berlin, origin used for the Dec 2017 showcase simulation
@@ -531,42 +534,12 @@
        CALL report('setup_parameters', message)
 
-!
-!--    Read in COSMO heights of half layers (vertical cell faces)
-       cosmo_var % name = 'HHL'
-       CALL get_netcdf_variable(cfg % hhl_file, cosmo_var, hhl)
-       CALL get_netcdf_dim_vector(cfg % hhl_file, 'rlon', rlon)
-       CALL get_netcdf_dim_vector(cfg % hhl_file, 'rlat', rlat)
-       CALL get_netcdf_dim_vector(soil_files(1), 'depth_2', depths)
  CALL run_control('time', 'read')
 
-       CALL reverse(hhl)
-       nlon = SIZE(hhl, 1)
-       nlat = SIZE(hhl, 2)
-       nlev = SIZE(hhl, 3)
-       ndepths = SIZE(depths)
-
-       lonmin_cosmo = MINVAL(rlon) * TO_RADIANS
-       lonmax_cosmo = MAXVAL(rlon) * TO_RADIANS
-       latmin_cosmo = MINVAL(rlat) * TO_RADIANS
-       latmax_cosmo = MAXVAL(rlat) * TO_RADIANS
- CALL run_control('time', 'comp')
-
-!
-!--    Appoximate COSMO-DE heights of full layers (cell centres)
-       ALLOCATE( hfl(nlon, nlat, nlev-1) )
-       ALLOCATE( d_depth(ndepths), d_depth_rho_inv(ndepths) )
-
- CALL run_control('time', 'alloc')
-       CALL get_soil_layer_thickness( depths, d_depth )
-       d_depth_rho_inv = 1.0_dp / ( d_depth * RHO_L )
-
-!
-!--    Appoximate COSMO-DE heights of full layers (cell centres)
-       DO k = 1, nlev-1
-          hfl(:,:,k) = 0.5_dp * ( hhl(:,:,k) +                                 &
-                                  hhl(:,:,k+1) )
-       ENDDO
- CALL run_control('time', 'comp')
-
+       CALL get_cosmo_grid( cfg, soil_files(1), rlon, rlat, hhl, hfl, depths,  &
+                            d_depth, d_depth_rho_inv, phi_n, lambda_n,         &
+                            phi_equat,                                         &
+                            lonmin_cosmo, lonmax_cosmo,                        &
+                            latmin_cosmo, latmax_cosmo,                        &
+                            nlon, nlat, nlev, ndepths )
 
 
@@ -680,4 +653,21 @@
     averaging_width_ns = averaging_angle * EARTH_RADIUS
 
+    lonmin_tot = MIN(lam_centre - averaging_angle, lonmin_palm)
+    lonmax_tot = MAX(lam_centre + averaging_angle, lonmax_palm)
+    latmin_tot = MIN(phi_centre - averaging_angle, latmin_palm)
+    latmax_tot = MAX(phi_centre + averaging_angle, latmax_palm)
+
+    palm_domain_outside_cosmo = ANY(                                           &
+       (/ lonmin_tot,   -lonmax_tot,   latmin_tot,   -latmax_tot/) .LT.        &
+       (/ lonmin_cosmo, -lonmax_cosmo, latmin_cosmo, -latmax_cosmo/)           &
+    )
+
+    IF ( palm_domain_outside_cosmo )  THEN
+       message = 'PALM domain or geostrophic averaging domains extend ' //     &
+                 'outside COSMO domain.'
+       CALL inifor_abort( 'setup_parameters', message )
+    ENDIF
+
+
 !
 !-- Coriolis parameter
@@ -686,5 +676,5 @@
                               phi_n * TO_DEGREES,                              &
                               gam * TO_DEGREES )                               &
-       )
+    )
 
     END SUBROUTINE setup_parameters
@@ -694,5 +684,5 @@
 ! Description:
 ! ------------
-!> Defines the COSMO, PALM-4U, PALM-4U boundary grids, in partucular their
+!> Defines the COSMO, PALM-4U, PALM-4U boundary grids, in particular their
 !> coordinates and interpolation weights
 !------------------------------------------------------------------------------!
@@ -1104,4 +1094,5 @@
                latmin = phi_south, latmax = phi_north,                         &
                kind='scalar', name='east geostrophic scalar')
+
 
 !                                                                              
@@ -4079,42 +4070,4 @@
 
 
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Fills the thickness array of the COSMO soil layers. Since COSMO's (i.e. 
-!> TERRA_ML's [1]) soil layer boundaries follow the rule
-!>
-!>    depth(0) = 0.0, and
-!>    depth(k) = 0.01 * 3**(k-1), k in [1,2,3,...,7]
-!>
-!> and full levels are defined as the midpoints between two layer boundaries, 
-!> all except the first layer thicknesses equal the depth of the midpoint.
-!> 
-!> [1] A Description of the Nonhydrostatic Regional COSMO Model Part II :
-!>     Physical Parameterization*, Sect. 11 TERRA_ML.
-!>     http://www.cosmo-model.org/content/model/documentation/core/cosmoPhysParamtr.pdf)
-!> 
-!> Input parameters:
-!> -----------------
-!>
-!> depths: array of full soil layer depths (cell midpoints)
-!>
-!>
-!> Output parameters:
-!> ------------------
-!>
-!> d_depth: array of soil layer thicknesses
-!>
-!------------------------------------------------------------------------------!
-    SUBROUTINE get_soil_layer_thickness(depths, d_depth)
-       
-       REAL(dp), INTENT(IN)  ::  depths(:)
-       REAL(dp), INTENT(OUT) ::  d_depth(:)
-
-       d_depth(:) = depths(:)
-       d_depth(1) = 2.0_dp * depths(1)
-
-    END SUBROUTINE get_soil_layer_thickness
-
  END MODULE inifor_grid
 #endif
Index: /palm/trunk/UTIL/inifor/src/inifor_io.f90
===================================================================
--- /palm/trunk/UTIL/inifor/src/inifor_io.f90	(revision 3800)
+++ /palm/trunk/UTIL/inifor/src/inifor_io.f90	(revision 3801)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Added routine get_cosmo_grid() to read in COSMO rotated pole from COSMO domain
+! Moved get_soil_layer_thickness() here from inifor_grid
+! 
+! 3785 2019-03-06 10:41:14Z eckhard
 ! Temporariliy disabled height-based geostrophic wind averaging
 ! Improved variable naming
@@ -41,6 +45,4 @@
 ! Moved get_input_file_list() here from grid module, added check for presence of
 !    input files
-!
-! 
 ! 
 ! 
@@ -116,5 +118,8 @@
     USE inifor_control
     USE inifor_defs,                                                           &
-        ONLY:  DATE, SNAME, PATH, PI, dp, hp, TO_RADIANS, TO_DEGREES, VERSION
+        ONLY:  DATE, SNAME, PATH, PI, dp, hp, TO_RADIANS, TO_DEGREES, VERSION, &
+               NC_DEPTH_NAME, NC_HHL_NAME, NC_RLAT_NAME, NC_RLON_NAME,         &
+               NC_ROTATED_POLE_NAME, NC_POLE_LATITUDE_NAME,                    &
+               NC_POLE_LONGITUDE_NAME, RHO_L
     USE inifor_types
     USE inifor_util,                                                           &
@@ -742,4 +747,122 @@
 
 
+   SUBROUTINE get_cosmo_grid( cfg, soil_file, rlon, rlat, hhl, hfl, depths, &
+                              d_depth, d_depth_rho_inv, phi_n, lambda_n,       &
+                              phi_equat,                                       &
+                              lonmin_cosmo, lonmax_cosmo,                      &
+                              latmin_cosmo, latmax_cosmo,                      &
+                              nlon, nlat, nlev, ndepths )
+
+      TYPE(inifor_config), INTENT(IN)                      ::  cfg
+      CHARACTER(LEN=PATH), INTENT(IN)                      ::  soil_file !< list of soil input files (temperature, moisture, <prefix>YYYYMMDDHH-soil.nc)
+      REAL(dp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  rlon      !< longitudes of COSMO-DE's rotated-pole grid
+      REAL(dp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  rlat      !< latitudes of COSMO-DE's rotated-pole grid
+      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(OUT) ::  hhl       !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file
+      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(OUT) ::  hfl       !< heights of full layers (cell centres) above sea level in COSMO-DE, computed as arithmetic average of hhl
+      REAL(dp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  depths    !< COSMO-DE's TERRA-ML soil layer depths
+      REAL(dp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  d_depth
+      REAL(dp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  d_depth_rho_inv
+      REAL(dp), INTENT(OUT)                                ::  phi_n
+      REAL(dp), INTENT(OUT)                                ::  phi_equat
+      REAL(dp), INTENT(OUT)                                ::  lambda_n
+      REAL(dp), INTENT(OUT)                                ::  lonmin_cosmo !< Minimunm longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
+      REAL(dp), INTENT(OUT)                                ::  lonmax_cosmo !< Maximum longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
+      REAL(dp), INTENT(OUT)                                ::  latmin_cosmo !< Minimunm latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
+      REAL(dp), INTENT(OUT)                                ::  latmax_cosmo !< Maximum latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
+      INTEGER, INTENt(OUT)                                 ::  nlon, nlat, nlev, ndepths
+
+      TYPE(nc_var) ::  cosmo_var !< COSMO dummy variable, used for reading HHL, rlon, rlat
+      INTEGER      ::  k
+
+!
+!--   Read in COSMO's heights of half layers (vertical cell faces)
+      cosmo_var % name = NC_HHL_NAME
+      CALL get_netcdf_variable( cfg % hhl_file, cosmo_var, hhl )
+      CALL get_netcdf_dim_vector( cfg % hhl_file, NC_RLON_NAME, rlon )
+      CALL get_netcdf_dim_vector( cfg % hhl_file, NC_RLAT_NAME, rlat )
+      CALL get_netcdf_dim_vector( soil_file, NC_DEPTH_NAME, depths)
+ CALL run_control( 'time', 'read' )
+
+      CALL reverse( hhl )
+      nlon = SIZE( hhl, 1 )
+      nlat = SIZE( hhl, 2 )
+      nlev = SIZE( hhl, 3 )
+      ndepths = SIZE( depths )
+
+ CALL run_control( 'time', 'comp' )
+
+      ALLOCATE( hfl( nlon, nlat, nlev-1 ) )
+      ALLOCATE( d_depth( ndepths ), d_depth_rho_inv( ndepths ) )
+ CALL run_control('time', 'alloc')
+
+      CALL get_soil_layer_thickness( depths, d_depth )
+      d_depth_rho_inv = 1.0_dp / ( d_depth * RHO_L )
+
+!
+!--   Compute COSMO's heights of full layers (cell centres)
+      DO k = 1, nlev-1
+         hfl(:,:,k) = 0.5_dp * ( hhl(:,:,k) +                                  &
+                                 hhl(:,:,k+1) )
+      ENDDO
+!
+!--   COSMO rotated pole coordinates
+      phi_n = TO_RADIANS                                                       &
+            * get_netcdf_variable_attribute( cfg % hhl_file,                   &
+                                             NC_ROTATED_POLE_NAME,             &
+                                             NC_POLE_LATITUDE_NAME )
+
+      lambda_n = TO_RADIANS                                                    &
+               * get_netcdf_variable_attribute( cfg % hhl_file,                &
+                                                NC_ROTATED_POLE_NAME,          &
+                                                NC_POLE_LONGITUDE_NAME )
+
+      phi_equat = 90.0_dp * TO_RADIANS - phi_n
+
+      lonmin_cosmo = MINVAL( rlon ) * TO_RADIANS
+      lonmax_cosmo = MAXVAL( rlon ) * TO_RADIANS
+      latmin_cosmo = MINVAL( rlat ) * TO_RADIANS
+      latmax_cosmo = MAXVAL( rlat ) * TO_RADIANS
+ CALL run_control('time', 'comp')
+
+   END SUBROUTINE get_cosmo_grid
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Fills the thickness array of the COSMO soil layers. Since COSMO's (i.e. 
+!> TERRA_ML's [1]) soil layer boundaries follow the rule
+!>
+!>    depth(0) = 0.0, and
+!>    depth(k) = 0.01 * 3**(k-1), k in [1,2,3,...,7]
+!>
+!> and full levels are defined as the midpoints between two layer boundaries, 
+!> all except the first layer thicknesses equal the depth of the midpoint.
+!> 
+!> [1] A Description of the Nonhydrostatic Regional COSMO Model Part II :
+!>     Physical Parameterization*, Sect. 11 TERRA_ML.
+!>     http://www.cosmo-model.org/content/model/documentation/core/cosmoPhysParamtr.pdf)
+!> 
+!> Input parameters:
+!> -----------------
+!>
+!> depths: array of full soil layer depths (cell midpoints)
+!>
+!>
+!> Output parameters:
+!> ------------------
+!>
+!> d_depth: array of soil layer thicknesses
+!>
+!------------------------------------------------------------------------------!
+    SUBROUTINE get_soil_layer_thickness(depths, d_depth)
+
+       REAL(dp), INTENT(IN)  ::  depths(:)
+       REAL(dp), INTENT(OUT) ::  d_depth(:)
+
+       d_depth(:) = depths(:)
+       d_depth(1) = 2.0_dp * depths(1)
+
+    END SUBROUTINE get_soil_layer_thickness
 !------------------------------------------------------------------------------!
 ! Description:
@@ -1191,4 +1314,33 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Read the attribute of the given variable form the given netCDF file.
+!------------------------------------------------------------------------------!
+    FUNCTION get_netcdf_variable_attribute(filename, varname, attribute)       &
+       RESULT(attribute_value)
+
+       CHARACTER(LEN=*), INTENT(IN) ::  filename, varname, attribute
+       REAL(dp)                     ::  attribute_value
+
+       INTEGER                      ::  ncid, varid
+
+       IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR )  THEN
+
+          CALL check( nf90_inq_varid( ncid, TRIM( varname ), varid ) )
+          CALL check( nf90_get_att( ncid, varid, TRIM( attribute ),            &
+                      attribute_value ) )
+          CALL check( nf90_close( ncid ) )
+
+       ELSE
+
+          message = "Failed to read '" // TRIM( varname ) // ":" //            &
+                    TRIM( attribute ) // "' from file '" // TRIM(filename) // "'."
+          CALL inifor_abort('get_netcdf_variable_attribute', message)
+
+       ENDIF
+
+    END FUNCTION get_netcdf_variable_attribute
 
 !------------------------------------------------------------------------------!
