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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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.