Changeset 3801


Ignore:
Timestamp:
Mar 15, 2019 5:14:25 PM (3 years ago)
Author:
eckhard
Message:

inifor: Read COSMO rotated pole from dataset, check if PALM domain exceeds COSMO domain

Location:
palm/trunk/UTIL/inifor
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/UTIL/inifor/README

    r3779 r3801  
    1 # INIFOR - Mesoscale Interface for Initializing and Forcing PALM-4U (v1.4.7)
     1# INIFOR - Mesoscale Interface for Initializing and Forcing PALM-4U
    22
    33INIFOR provides the meteorological fields required to initialize and drive the
  • palm/trunk/UTIL/inifor/src/inifor_defs.f90

    r3779 r3801  
    2626! -----------------
    2727! $Id$
     28! Defined netCDF variable names for COSMO grid
     29! Bumped version number
     30!
     31!
     32! 3779 2019-03-05 11:13:35Z eckhard
    2833! Updated version number to 1.4.7, updated copyright note
    2934!
     
    111116
    112117!
    113 !-- COSMO-DE parameters
     118!-- COSMO parameters
    114119 INTEGER, PARAMETER  ::  WATER_ID = 9                !< Integer corresponding to the water soil type in COSMO-DE [-]
    115120 REAL(dp), PARAMETER ::  EARTH_RADIUS = 6371229.0_dp !< Earth radius used in COSMO-DE [m]
     
    134139 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]
    135140
     141!
     142!-- COSMO netCDF names
     143CHARACTER(SNAME) ::  NC_DEPTH_NAME = 'depth_2'
     144CHARACTER(SNAME) ::  NC_HHL_NAME = 'HHL'
     145CHARACTER(SNAME) ::  NC_RLAT_NAME = 'rlat'
     146CHARACTER(SNAME) ::  NC_RLON_NAME = 'rlon'
     147CHARACTER(SNAME) ::  NC_ROTATED_POLE_NAME = 'rotated_pole'
     148CHARACTER(SNAME) ::  NC_POLE_LATITUDE_NAME = 'grid_north_pole_latitude'
     149CHARACTER(SNAME) ::  NC_POLE_LONGITUDE_NAME = 'grid_north_pole_longitude'
     150
    136151!
    137152!-- INIFOR parameters
     
    140155 INTEGER, PARAMETER          ::  FORCING_STEP = 1             !< Number of hours between forcing time steps [h]
    141156 REAL(dp), PARAMETER         ::  NUDGING_TAU = 21600.0_dp     !< Nudging relaxation time scale [s]
    142  CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.4.7'            !< INIFOR version number
     157 CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.4.8'            !< INIFOR version number
    143158 CHARACTER(LEN=*), PARAMETER ::  COPYRIGHT = 'Copyright 2017-2019 Leibniz Universitaet Hannover' // &
    144159     ACHAR( 10 ) // ' Copyright 2017-2019 Deutscher Wetterdienst Offenbach' !< Copyright notice
  • palm/trunk/UTIL/inifor/src/inifor_grid.f90

    r3785 r3801  
    2626! -----------------
    2727! $Id$
     28! Read COSMO rotated pole from HHL file
     29! Check for PALM and averaging domains extending outside COSMO domain
     30!
     31! 3785 2019-03-06 10:41:14Z eckhard
    2832! Assigned names to averaging grids
    2933! Improved variable naming and minor clean-up
     
    129133               RHO_L, OMEGA, HECTO
    130134    USE inifor_io,                                                             &
    131         ONLY:  get_netcdf_attribute, get_netcdf_dim_vector,                    &
     135        ONLY:  get_cosmo_grid, get_netcdf_attribute, get_netcdf_dim_vector,    &
    132136               get_netcdf_variable, parse_command_line_arguments,              &
    133137               get_input_file_list, validate_config
     
    195199    REAL(dp) ::  latmin_palm       = 0.0_dp       !< Minimunm latitude of PALM grid [COSMO rotated-pole rad]
    196200    REAL(dp) ::  latmax_palm       = 0.0_dp       !< Maximum latitude of PALM grid [COSMO rotated-pole rad]
     201    REAL(dp) ::  lonmin_tot        = 0.0_dp       !< Minimunm longitude of required COSMO data [COSMO rotated-pole rad]
     202    REAL(dp) ::  lonmax_tot        = 0.0_dp       !< Maximum longitude of required COSMO data [COSMO rotated-pole rad]
     203    REAL(dp) ::  latmin_tot        = 0.0_dp       !< Minimunm latitude of required COSMO data [COSMO rotated-pole rad]
     204    REAL(dp) ::  latmax_tot        = 0.0_dp       !< Maximum latitude of required COSMO data [COSMO rotated-pole rad]
    197205    REAL(dp) ::  latitude          = 0.0_dp       !< geographical latitude of the PALM-4U origin, from inipar namelist [deg]
    198206    REAL(dp) ::  longitude         = 0.0_dp       !< geographical longitude of the PALM-4U origin, from inipar namelist [deg]
     
    242250    LOGICAL ::  ls_forcing_variables_required !< flag controlling whether large-scale forcing variables are to be computed
    243251    LOGICAL ::  surface_forcing_required      !< flag controlling whether surface forcing variables are to be computed
     252    LOGICAL ::  palm_domain_outside_cosmo     !< indicates whether COSMO grid covers the PALM domain and the geostrophic averaging domains
    244253
    245254    TYPE(nc_var), ALLOCATABLE, TARGET ::  input_var_table(:)  !< table of input variables
    246255    TYPE(nc_var), ALLOCATABLE, TARGET ::  output_var_table(:) !< table of input variables
    247     TYPE(nc_var)                      ::  cosmo_var           !< COSMO dummy variable, used for reading HHL, rlon, rlat
     256    TYPE(nc_var) ::  cosmo_var                                !< COSMO dummy variable, used for reading HHL, rlon, rlat
    248257
    249258    TYPE(grid_definition), TARGET ::  palm_grid                       !< PALM-4U grid in the target system (COSMO-DE rotated-pole)
     
    356365
    357366!
    358 !--    COSMO-DE and -D2 rotated pole position
    359        phi_n     =   40.0_dp * TO_RADIANS
    360        phi_equat =   50.0_dp * TO_RADIANS
    361        lambda_n  = -170.0_dp * TO_RADIANS
    362 
    363 !
    364367!--    Defaultmain centre (_c) of the PALM-4U grid in the geographical system (_g)
    365368       origin_lat = 52.325079_dp * TO_RADIANS ! south-west of Berlin, origin used for the Dec 2017 showcase simulation
     
    531534       CALL report('setup_parameters', message)
    532535
    533 !
    534 !--    Read in COSMO heights of half layers (vertical cell faces)
    535        cosmo_var % name = 'HHL'
    536        CALL get_netcdf_variable(cfg % hhl_file, cosmo_var, hhl)
    537        CALL get_netcdf_dim_vector(cfg % hhl_file, 'rlon', rlon)
    538        CALL get_netcdf_dim_vector(cfg % hhl_file, 'rlat', rlat)
    539        CALL get_netcdf_dim_vector(soil_files(1), 'depth_2', depths)
    540536 CALL run_control('time', 'read')
    541537
    542        CALL reverse(hhl)
    543        nlon = SIZE(hhl, 1)
    544        nlat = SIZE(hhl, 2)
    545        nlev = SIZE(hhl, 3)
    546        ndepths = SIZE(depths)
    547 
    548        lonmin_cosmo = MINVAL(rlon) * TO_RADIANS
    549        lonmax_cosmo = MAXVAL(rlon) * TO_RADIANS
    550        latmin_cosmo = MINVAL(rlat) * TO_RADIANS
    551        latmax_cosmo = MAXVAL(rlat) * TO_RADIANS
    552  CALL run_control('time', 'comp')
    553 
    554 !
    555 !--    Appoximate COSMO-DE heights of full layers (cell centres)
    556        ALLOCATE( hfl(nlon, nlat, nlev-1) )
    557        ALLOCATE( d_depth(ndepths), d_depth_rho_inv(ndepths) )
    558 
    559  CALL run_control('time', 'alloc')
    560        CALL get_soil_layer_thickness( depths, d_depth )
    561        d_depth_rho_inv = 1.0_dp / ( d_depth * RHO_L )
    562 
    563 !
    564 !--    Appoximate COSMO-DE heights of full layers (cell centres)
    565        DO k = 1, nlev-1
    566           hfl(:,:,k) = 0.5_dp * ( hhl(:,:,k) +                                 &
    567                                   hhl(:,:,k+1) )
    568        ENDDO
    569  CALL run_control('time', 'comp')
    570 
     538       CALL get_cosmo_grid( cfg, soil_files(1), rlon, rlat, hhl, hfl, depths,  &
     539                            d_depth, d_depth_rho_inv, phi_n, lambda_n,         &
     540                            phi_equat,                                         &
     541                            lonmin_cosmo, lonmax_cosmo,                        &
     542                            latmin_cosmo, latmax_cosmo,                        &
     543                            nlon, nlat, nlev, ndepths )
    571544
    572545
     
    680653    averaging_width_ns = averaging_angle * EARTH_RADIUS
    681654
     655    lonmin_tot = MIN(lam_centre - averaging_angle, lonmin_palm)
     656    lonmax_tot = MAX(lam_centre + averaging_angle, lonmax_palm)
     657    latmin_tot = MIN(phi_centre - averaging_angle, latmin_palm)
     658    latmax_tot = MAX(phi_centre + averaging_angle, latmax_palm)
     659
     660    palm_domain_outside_cosmo = ANY(                                           &
     661       (/ lonmin_tot,   -lonmax_tot,   latmin_tot,   -latmax_tot/) .LT.        &
     662       (/ lonmin_cosmo, -lonmax_cosmo, latmin_cosmo, -latmax_cosmo/)           &
     663    )
     664
     665    IF ( palm_domain_outside_cosmo )  THEN
     666       message = 'PALM domain or geostrophic averaging domains extend ' //     &
     667                 'outside COSMO domain.'
     668       CALL inifor_abort( 'setup_parameters', message )
     669    ENDIF
     670
     671
    682672!
    683673!-- Coriolis parameter
     
    686676                              phi_n * TO_DEGREES,                              &
    687677                              gam * TO_DEGREES )                               &
    688        )
     678    )
    689679
    690680    END SUBROUTINE setup_parameters
     
    694684! Description:
    695685! ------------
    696 !> Defines the COSMO, PALM-4U, PALM-4U boundary grids, in partucular their
     686!> Defines the COSMO, PALM-4U, PALM-4U boundary grids, in particular their
    697687!> coordinates and interpolation weights
    698688!------------------------------------------------------------------------------!
     
    11041094               latmin = phi_south, latmax = phi_north,                         &
    11051095               kind='scalar', name='east geostrophic scalar')
     1096
    11061097
    11071098!                                                                             
     
    40794070
    40804071
    4081 !------------------------------------------------------------------------------!
    4082 ! Description:
    4083 ! ------------
    4084 !> Fills the thickness array of the COSMO soil layers. Since COSMO's (i.e.
    4085 !> TERRA_ML's [1]) soil layer boundaries follow the rule
    4086 !>
    4087 !>    depth(0) = 0.0, and
    4088 !>    depth(k) = 0.01 * 3**(k-1), k in [1,2,3,...,7]
    4089 !>
    4090 !> and full levels are defined as the midpoints between two layer boundaries,
    4091 !> all except the first layer thicknesses equal the depth of the midpoint.
    4092 !>
    4093 !> [1] A Description of the Nonhydrostatic Regional COSMO Model Part II :
    4094 !>     Physical Parameterization*, Sect. 11 TERRA_ML.
    4095 !>     http://www.cosmo-model.org/content/model/documentation/core/cosmoPhysParamtr.pdf)
    4096 !>
    4097 !> Input parameters:
    4098 !> -----------------
    4099 !>
    4100 !> depths: array of full soil layer depths (cell midpoints)
    4101 !>
    4102 !>
    4103 !> Output parameters:
    4104 !> ------------------
    4105 !>
    4106 !> d_depth: array of soil layer thicknesses
    4107 !>
    4108 !------------------------------------------------------------------------------!
    4109     SUBROUTINE get_soil_layer_thickness(depths, d_depth)
    4110        
    4111        REAL(dp), INTENT(IN)  ::  depths(:)
    4112        REAL(dp), INTENT(OUT) ::  d_depth(:)
    4113 
    4114        d_depth(:) = depths(:)
    4115        d_depth(1) = 2.0_dp * depths(1)
    4116 
    4117     END SUBROUTINE get_soil_layer_thickness
    4118 
    41194072 END MODULE inifor_grid
    41204073#endif
  • palm/trunk/UTIL/inifor/src/inifor_io.f90

    r3785 r3801  
    2626! -----------------
    2727! $Id$
     28! Added routine get_cosmo_grid() to read in COSMO rotated pole from COSMO domain
     29! Moved get_soil_layer_thickness() here from inifor_grid
     30!
     31! 3785 2019-03-06 10:41:14Z eckhard
    2832! Temporariliy disabled height-based geostrophic wind averaging
    2933! Improved variable naming
     
    4145! Moved get_input_file_list() here from grid module, added check for presence of
    4246!    input files
    43 !
    44 !
    4547!
    4648!
     
    116118    USE inifor_control
    117119    USE inifor_defs,                                                           &
    118         ONLY:  DATE, SNAME, PATH, PI, dp, hp, TO_RADIANS, TO_DEGREES, VERSION
     120        ONLY:  DATE, SNAME, PATH, PI, dp, hp, TO_RADIANS, TO_DEGREES, VERSION, &
     121               NC_DEPTH_NAME, NC_HHL_NAME, NC_RLAT_NAME, NC_RLON_NAME,         &
     122               NC_ROTATED_POLE_NAME, NC_POLE_LATITUDE_NAME,                    &
     123               NC_POLE_LONGITUDE_NAME, RHO_L
    119124    USE inifor_types
    120125    USE inifor_util,                                                           &
     
    742747
    743748
     749   SUBROUTINE get_cosmo_grid( cfg, soil_file, rlon, rlat, hhl, hfl, depths, &
     750                              d_depth, d_depth_rho_inv, phi_n, lambda_n,       &
     751                              phi_equat,                                       &
     752                              lonmin_cosmo, lonmax_cosmo,                      &
     753                              latmin_cosmo, latmax_cosmo,                      &
     754                              nlon, nlat, nlev, ndepths )
     755
     756      TYPE(inifor_config), INTENT(IN)                      ::  cfg
     757      CHARACTER(LEN=PATH), INTENT(IN)                      ::  soil_file !< list of soil input files (temperature, moisture, <prefix>YYYYMMDDHH-soil.nc)
     758      REAL(dp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  rlon      !< longitudes of COSMO-DE's rotated-pole grid
     759      REAL(dp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  rlat      !< latitudes of COSMO-DE's rotated-pole grid
     760      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(OUT) ::  hhl       !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file
     761      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
     762      REAL(dp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  depths    !< COSMO-DE's TERRA-ML soil layer depths
     763      REAL(dp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  d_depth
     764      REAL(dp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  d_depth_rho_inv
     765      REAL(dp), INTENT(OUT)                                ::  phi_n
     766      REAL(dp), INTENT(OUT)                                ::  phi_equat
     767      REAL(dp), INTENT(OUT)                                ::  lambda_n
     768      REAL(dp), INTENT(OUT)                                ::  lonmin_cosmo !< Minimunm longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
     769      REAL(dp), INTENT(OUT)                                ::  lonmax_cosmo !< Maximum longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
     770      REAL(dp), INTENT(OUT)                                ::  latmin_cosmo !< Minimunm latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
     771      REAL(dp), INTENT(OUT)                                ::  latmax_cosmo !< Maximum latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
     772      INTEGER, INTENt(OUT)                                 ::  nlon, nlat, nlev, ndepths
     773
     774      TYPE(nc_var) ::  cosmo_var !< COSMO dummy variable, used for reading HHL, rlon, rlat
     775      INTEGER      ::  k
     776
     777!
     778!--   Read in COSMO's heights of half layers (vertical cell faces)
     779      cosmo_var % name = NC_HHL_NAME
     780      CALL get_netcdf_variable( cfg % hhl_file, cosmo_var, hhl )
     781      CALL get_netcdf_dim_vector( cfg % hhl_file, NC_RLON_NAME, rlon )
     782      CALL get_netcdf_dim_vector( cfg % hhl_file, NC_RLAT_NAME, rlat )
     783      CALL get_netcdf_dim_vector( soil_file, NC_DEPTH_NAME, depths)
     784 CALL run_control( 'time', 'read' )
     785
     786      CALL reverse( hhl )
     787      nlon = SIZE( hhl, 1 )
     788      nlat = SIZE( hhl, 2 )
     789      nlev = SIZE( hhl, 3 )
     790      ndepths = SIZE( depths )
     791
     792 CALL run_control( 'time', 'comp' )
     793
     794      ALLOCATE( hfl( nlon, nlat, nlev-1 ) )
     795      ALLOCATE( d_depth( ndepths ), d_depth_rho_inv( ndepths ) )
     796 CALL run_control('time', 'alloc')
     797
     798      CALL get_soil_layer_thickness( depths, d_depth )
     799      d_depth_rho_inv = 1.0_dp / ( d_depth * RHO_L )
     800
     801!
     802!--   Compute COSMO's heights of full layers (cell centres)
     803      DO k = 1, nlev-1
     804         hfl(:,:,k) = 0.5_dp * ( hhl(:,:,k) +                                  &
     805                                 hhl(:,:,k+1) )
     806      ENDDO
     807!
     808!--   COSMO rotated pole coordinates
     809      phi_n = TO_RADIANS                                                       &
     810            * get_netcdf_variable_attribute( cfg % hhl_file,                   &
     811                                             NC_ROTATED_POLE_NAME,             &
     812                                             NC_POLE_LATITUDE_NAME )
     813
     814      lambda_n = TO_RADIANS                                                    &
     815               * get_netcdf_variable_attribute( cfg % hhl_file,                &
     816                                                NC_ROTATED_POLE_NAME,          &
     817                                                NC_POLE_LONGITUDE_NAME )
     818
     819      phi_equat = 90.0_dp * TO_RADIANS - phi_n
     820
     821      lonmin_cosmo = MINVAL( rlon ) * TO_RADIANS
     822      lonmax_cosmo = MAXVAL( rlon ) * TO_RADIANS
     823      latmin_cosmo = MINVAL( rlat ) * TO_RADIANS
     824      latmax_cosmo = MAXVAL( rlat ) * TO_RADIANS
     825 CALL run_control('time', 'comp')
     826
     827   END SUBROUTINE get_cosmo_grid
     828
     829
     830!------------------------------------------------------------------------------!
     831! Description:
     832! ------------
     833!> Fills the thickness array of the COSMO soil layers. Since COSMO's (i.e.
     834!> TERRA_ML's [1]) soil layer boundaries follow the rule
     835!>
     836!>    depth(0) = 0.0, and
     837!>    depth(k) = 0.01 * 3**(k-1), k in [1,2,3,...,7]
     838!>
     839!> and full levels are defined as the midpoints between two layer boundaries,
     840!> all except the first layer thicknesses equal the depth of the midpoint.
     841!>
     842!> [1] A Description of the Nonhydrostatic Regional COSMO Model Part II :
     843!>     Physical Parameterization*, Sect. 11 TERRA_ML.
     844!>     http://www.cosmo-model.org/content/model/documentation/core/cosmoPhysParamtr.pdf)
     845!>
     846!> Input parameters:
     847!> -----------------
     848!>
     849!> depths: array of full soil layer depths (cell midpoints)
     850!>
     851!>
     852!> Output parameters:
     853!> ------------------
     854!>
     855!> d_depth: array of soil layer thicknesses
     856!>
     857!------------------------------------------------------------------------------!
     858    SUBROUTINE get_soil_layer_thickness(depths, d_depth)
     859
     860       REAL(dp), INTENT(IN)  ::  depths(:)
     861       REAL(dp), INTENT(OUT) ::  d_depth(:)
     862
     863       d_depth(:) = depths(:)
     864       d_depth(1) = 2.0_dp * depths(1)
     865
     866    END SUBROUTINE get_soil_layer_thickness
    744867!------------------------------------------------------------------------------!
    745868! Description:
     
    11911314
    11921315
     1316!------------------------------------------------------------------------------!
     1317! Description:
     1318! ------------
     1319!> Read the attribute of the given variable form the given netCDF file.
     1320!------------------------------------------------------------------------------!
     1321    FUNCTION get_netcdf_variable_attribute(filename, varname, attribute)       &
     1322       RESULT(attribute_value)
     1323
     1324       CHARACTER(LEN=*), INTENT(IN) ::  filename, varname, attribute
     1325       REAL(dp)                     ::  attribute_value
     1326
     1327       INTEGER                      ::  ncid, varid
     1328
     1329       IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR )  THEN
     1330
     1331          CALL check( nf90_inq_varid( ncid, TRIM( varname ), varid ) )
     1332          CALL check( nf90_get_att( ncid, varid, TRIM( attribute ),            &
     1333                      attribute_value ) )
     1334          CALL check( nf90_close( ncid ) )
     1335
     1336       ELSE
     1337
     1338          message = "Failed to read '" // TRIM( varname ) // ":" //            &
     1339                    TRIM( attribute ) // "' from file '" // TRIM(filename) // "'."
     1340          CALL inifor_abort('get_netcdf_variable_attribute', message)
     1341
     1342       ENDIF
     1343
     1344    END FUNCTION get_netcdf_variable_attribute
    11931345
    11941346!------------------------------------------------------------------------------!
Note: See TracChangeset for help on using the changeset viewer.