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

inifor: Updated documentation

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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   
Note: See TracChangeset for help on using the changeset viewer.