Changeset 3557 for palm/trunk/UTIL/inifor/src/inifor_transform.f90

Timestamp:

Nov 22, 2018 4:01:22 PM (5 years ago)

Author:

eckhard

Message:

inifor: Updated documentation

File:

• palm/trunk/UTIL/inifor/src/inifor_transform.f90 (modified) (47 diffs)

palm/trunk/UTIL/inifor/src/inifor_transform.f90

@@ -26 +26 @@
 ! -----------------
 ! $Id$
+! Updated documentation
+! 
+! 
+! 3537 2018-11-20 10:53:14Z eckhard
 ! bugfix: working precision added
 ! 
@@ -52 +56 @@
 ! Authors:
 ! --------
-! @author Eckhard Kadasch
+!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
 !
 ! Description:
@@ -110 +114 @@
        DO k = nz, LBOUND(out_arr, 3), -1
 
-          ! TODO: Remove IF clause and extrapolate based on a critical vertical
-          ! TODO: index marking the lower bound of COSMO-DE data coverage.
-          ! Check for negative interpolation weights indicating grid points 
-          ! below COSMO-DE domain and extrapolate from the top in such cells.
+!
+!--       TODO: Remove IF clause and extrapolate based on a critical vertical
+!--       TODO: index marking the lower bound of COSMO-DE data coverage.
+!--       Check for negative interpolation weights indicating grid points 
+!--       below COSMO-DE domain and extrapolate from the top in such cells.
           IF (outgrid % w_verti(i,j,k,1) < -1.0_dp .AND. k < nz)  THEN
              out_arr(i,j,k) = out_arr(i,j,k+1)
@@ -155 +160 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE interpolate_2d(invar, outvar, outgrid, ncvar)
-    ! I index 0-based for the indices of the outvar to be consistent with the
-    ! outgrid indices and interpolation weights.
+!
+!--    I index 0-based for the indices of the outvar to be consistent with the
+!--    outgrid indices and interpolation weights.
        TYPE(grid_definition), INTENT(IN)  ::  outgrid
        REAL(dp), INTENT(IN)               ::  invar(0:,0:,0:)
@@ -164 +170 @@
        INTEGER ::  i, j, k, l
 
-       ! TODO: check if input dimensions are consistent, i.e. ranges are correct
-       IF (UBOUND(outvar, 3) .GT. UBOUND(invar, 3))  THEN
+!
+!--    TODO: check if input dimensions are consistent, i.e. ranges are correct
+       IF ( UBOUND(outvar, 3) .GT. UBOUND(invar, 3) )  THEN
            message = "Output array for '" // TRIM(ncvar % name) // "' has ' more levels (" // &
               TRIM(str(UBOUND(outvar, 3))) // ") than input variable ("//&
@@ -190 +197 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Compute the horizontal average of the in_arr(:,:,:) and store it in
+!> out_arr(:)
+!------------------------------------------------------------------------------!
     SUBROUTINE average_2d(in_arr, out_arr, ii, jj)
        REAL(dp), INTENT(IN)              ::  in_arr(0:,0:,0:)
@@ -200 +213 @@
        IF (SIZE(ii) .NE. SIZE(jj))  THEN
           message = "Length of 'ii' and 'jj' index lists do not match." //     &
-             NEW_LINE(' ') // "ii has " // str(SIZE(ii)) // " elements, " //        &
+             NEW_LINE(' ') // "ii has " // str(SIZE(ii)) // " elements, " //   &
              NEW_LINE(' ') // "jj has " // str(SIZE(jj)) // "."
           CALL abort('average_2d', message)
@@ -211 +224 @@
              i = ii(l)
              j = jj(l)
-             out_arr(k) = out_arr(k) +&
-                         in_arr(i, j, k)
+             out_arr(k) = out_arr(k) + in_arr(i, j, k)
           END DO
 
@@ -223 +235 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Three-dimensional interpolation driver. Interpolates from the source_array to
+!> the given palm_grid.
+!>
+!> The routine separates horizontal and vertical
+!> interpolation. In the horizontal interpolation step, the source_array data is
+!> interpolated along COSMO grid levels onto the intermediate grid (vertically
+!> as coarse as COSMO, horizontally as fine as PALM).
+!------------------------------------------------------------------------------!
     SUBROUTINE interpolate_3d(source_array, palm_array, palm_intermediate, palm_grid)
        TYPE(grid_definition), INTENT(IN) ::  palm_intermediate, palm_grid
@@ -232 +255 @@
        nx = palm_intermediate % nx
        ny = palm_intermediate % ny
-       nlev = palm_intermediate % nz ! nlev
-
-       ! Interpolate from COSMO-DE to intermediate grid. Allocating with one
-       ! less point in the vertical, since scalars like T have 50 instead of 51
-       ! points in COSMO-DE.
+       nlev = palm_intermediate % nz
+
+!
+!--    Interpolate from COSMO to intermediate grid. Allocating with one
+!--    less point in the vertical, since scalars like T have 50 instead of 51
+!--    points in COSMO.
        ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nlev-1)) !
 
        CALL interpolate_2d(source_array, intermediate_array, palm_intermediate)
 
-       ! Interpolate from intermediate grid to palm_grid grid, includes
-       ! extrapolation for cells below COSMO-DE domain.
+!
+!--    Interpolate from intermediate grid to palm_grid grid, includes
+!--    extrapolation for cells below COSMO domain.
        CALL interpolate_1d(intermediate_array, palm_array, palm_grid)
 
@@ -250 +275 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Average data horizontally from the source_array over the region given by the
+!> averaging grid 'avg_grid' and store the result in 'profile_array'.
+!------------------------------------------------------------------------------!
     SUBROUTINE average_profile(source_array, profile_array, avg_grid)
        TYPE(grid_definition), INTENT(IN)          ::  avg_grid
@@ -265 +296 @@
           j_source = avg_grid % jjj(l)
 
-          DO k_profile = avg_grid % k_min, UBOUND(profile_array, 1) ! PALM levels
-
-             DO m = 1, 2 ! vertical interpolation neighbours
+!
+!--       Loop over PALM levels
+          DO k_profile = avg_grid % k_min, UBOUND(profile_array, 1)
+
+!
+!--          Loop over vertical interpolation neighbours
+             DO m = 1, 2
 
                 k_source = avg_grid % kkk(l, k_profile, m)
@@ -274 +309 @@
                    + avg_grid % w(l, k_profile, m)                             &
                    * source_array(i_source, j_source, k_source)
-
-             END DO ! m, vertical interpolation neighbours
-
-          END DO    ! k_profile, PALM levels
-
-       END DO       ! l, horizontal neighbours
+!
+!--          Loop over vertical interpolation neighbours m
+             END DO
+
+!
+!--       Loop over PALM levels k_profile
+          END DO
+
+!
+!--    Loop over horizontal neighbours l
+       END DO
 
        ni_columns = 1.0_dp / avg_grid % n_columns
        profile_array(:) = profile_array(:) * ni_columns
 
-       ! Extrapolate constant to the bottom
+!
+!--    Constant extrapolation to the bottom
        profile_array(1:avg_grid % k_min-1) = profile_array(avg_grid % k_min)
 
@@ -290 +331 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Extrapolates density linearly from the level 'k_min' downwards.
+!------------------------------------------------------------------------------!
     SUBROUTINE extrapolate_density(rho, avg_grid)
        REAL(dp), DIMENSION(:), INTENT(INOUT) ::  rho
@@ -308 +354 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Driver for extrapolating pressure from PALM level k_min downwards
+!------------------------------------------------------------------------------!
     SUBROUTINE extrapolate_pressure(p, rho, avg_grid)
        REAL(dp), DIMENSION(:), INTENT(IN)    ::  rho
@@ -405 +456 @@
        REAL(dp) :: ri
 
-       ! TODO check dimensions of lat/lon and y/x match
+!
+!--    TODO check dimensions of lat/lon and y/x match
 
        ri = 1.0_dp / r
@@ -438 +490 @@
        REAL(dp) :: ri
 
-       ! If this elemental function is called with a large array as xy, it is
-       ! computationally more efficient to precompute the inverse radius and
-       ! then muliply.
+!
+!--    If this elemental function is called with a large array as xy, it is
+!--    computationally more efficient to precompute the inverse radius and
+!--    then muliply.
        ri = 1.0_dp / r
 
@@ -448 +501 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> For a rotated-pole system with the origin at geographical latitude 'phi_c',
+!> compute the geographical latitude of its rotated north pole.
+!------------------------------------------------------------------------------!
     REAL(dp) FUNCTION phic_to_phin(phi_c)
         REAL(dp), INTENT(IN) ::  phi_c
@@ -456 +515 @@
 
     
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> For a rotated-pole system with the origin at geographical latitude 'phi_c'
+!> and longitude 'lam_c', compute the geographical longitude of its rotated
+!> north pole.
+!------------------------------------------------------------------------------!
     REAL(dp) FUNCTION lamc_to_lamn(phi_c, lam_c)
        REAL(dp), INTENT(IN) ::  phi_c, lam_c
@@ -467 +533 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Set gamma according to whether PALM domain is in the northern or southern
+!> hemisphere of the COSMO rotated-pole system. Gamma assumes either the
+!> value 0 or PI and is needed to work around around a bug in the
+!> rotated-pole coordinate transformations.
+!------------------------------------------------------------------------------!
     REAL(dp) FUNCTION gamma_from_hemisphere(phi_cg, phi_ref)
-       REAL(dp), INTENT(IN) ::  phi_cg, phi_ref
-       LOGICAL ::  palm_centre_is_south_of_cosmo_origin
-       
-       palm_centre_is_south_of_cosmo_origin = (phi_cg < phi_ref)
-
-       IF (palm_centre_is_south_of_cosmo_origin)  THEN
+       REAL(dp), INTENT(IN) ::  phi_cg
+       REAL(dp), INTENT(IN) ::  phi_ref
+
+       LOGICAL ::  palm_origin_is_south_of_cosmo_origin
+       
+       palm_origin_is_south_of_cosmo_origin = (phi_cg < phi_ref)
+
+       IF (palm_origin_is_south_of_cosmo_origin)  THEN
            gamma_from_hemisphere = PI
        ELSE
@@ -550 +626 @@
        
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Rotate the given vector field (x(:), y(:)) by the given 'angle'.
+!------------------------------------------------------------------------------!
     SUBROUTINE rotate_vector_field(x, y, angle)
        REAL(dp), DIMENSION(:), INTENT(INOUT) :: x, y  !< x and y coodrinate in arbitrary units
@@ -559 +640 @@
        sine = SIN(angle * TO_RADIANS)
        cosine = COS(angle * TO_RADIANS)
-       ! RESAHPE() fills columns first, so the rotation matrix becomes
-       !
-       ! rotation = [ cosine   -sine  ]
-       !            [  sine    cosine ]
+!
+!--    RESAHPE() fills columns first, so the rotation matrix becomes
+!--    
+!--    rotation = [ cosine   -sine  ]
+!--               [  sine    cosine ]
        rotation = RESHAPE( (/cosine, sine, -sine, cosine/), (/2, 2/) )
 
@@ -588 +670 @@
 !>    Datenbanken des DWD.
 !>    https://www.dwd.de/SharedDocs/downloads/DE/modelldokumentationen/nwv/cosmo_d2/cosmo_d2_dbbeschr_aktuell.pdf?__blob=publicationFile&v=2
-!>
+!------------------------------------------------------------------------------!
     FUNCTION meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g)          &
        RESULT(delta)
@@ -664 +746 @@
        DO j = 0, UBOUND(palm_clon, 2)!palm_grid % ny
        DO i = 0, UBOUND(palm_clon, 1)!palm_grid % nx
-          ! Compute the floating point index corrseponding to PALM-4U grid point
-          ! location along target grid (COSMO-DE) axes.
+!
+!--       Compute the floating point index corrseponding to PALM-4U grid point
+!--       location along target grid (COSMO-DE) axes.
           lonpos = (palm_clon(i,j) - lon0) * cosmo_dxi
           latpos = (palm_clat(i,j) - lat0) * cosmo_dyi
@@ -693 +776 @@
 
     
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Computes linear vertical interpolation neighbour indices and weights for each
+!> column of the given palm grid.
+!------------------------------------------------------------------------------!
     SUBROUTINE find_vertical_neighbours_and_weights_interp( palm_grid,         &
                                                             palm_intermediate )
@@ -709 +798 @@
        nlev = palm_intermediate % nz
 
-       ! in each column of the fine grid, find vertical neighbours of every cell
+!
+!--    in each column of the fine grid, find vertical neighbours of every cell
        DO j = 0, ny
        DO i = 0, nx
@@ -718 +808 @@
           column_top  = palm_intermediate % h(i,j,nlev)
 
-          ! scan through palm_grid column and set neighbour indices in
-          ! case current_height is either below column_base, in the current
-          ! cell, or above column_top. Keep increasing current cell index until
-          ! the current cell overlaps with the current_height.
+!
+!--       scan through palm_grid column and set neighbour indices in
+!--       case current_height is either below column_base, in the current
+!--       cell, or above column_top. Keep increasing current cell index until
+!--       the current cell overlaps with the current_height.
           DO k = 1, nz
 
-             ! Memorize the top and bottom boundaries of the coarse cell and the
-             ! current height within it
+!
+!--          Memorize the top and bottom boundaries of the coarse cell and the
+!--          current height within it
              current_height = palm_grid % z(k) + palm_grid % z0
              h_top    = palm_intermediate % h(i,j,k_intermediate+1)
@@ -738 +830 @@
              )
 
-             ! set default weights
+!
+!--          set default weights
              palm_grid % w_verti(i,j,k,1:2) = 0.0_dp
 
@@ -755 +848 @@
 
              ELSE
-                ! cycle through intermediate levels until current
-                ! intermediate-grid cell overlaps with current_height
+!
+!--             cycle through intermediate levels until current
+!--             intermediate-grid cell overlaps with current_height
                 DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
                    k_intermediate = k_intermediate + 1
@@ -777 +871 @@
                 palm_grid % kk(i,j,k,2) = k_intermediate + 1
 
-                ! copmute vertical weights
+!
+!--             compute vertical weights
                 weight = (h_top - current_height) / (h_top - h_bottom)
                 palm_grid % w_verti(i,j,k,1) = weight
@@ -791 +886 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Computes linear vertical interpolation neighbour indices and weights for each
+!> column of the given averaging grid.
+!>
+!> The difference to the _interp variant of this routine lies in how columns
+!> are adressed. While the _interp variant loops over all PALM grid columns
+!> given by combinations of all index indices (i,j), this routine loops over a
+!> subset of COSMO columns, which are sequantlially stored in the index lists
+!> iii(:) and jjj(:).
+!------------------------------------------------------------------------------!
     SUBROUTINE find_vertical_neighbours_and_weights_average( avg_grid )
        TYPE(grid_definition), INTENT(INOUT) ::  avg_grid
@@ -805 +912 @@
        nlev = SIZE(avg_grid % cosmo_h, 3)
 
-       ! in each column of the fine grid, find vertical neighbours of every cell
+!
+!--    in each column of the fine grid, find vertical neighbours of every cell
        DO l = 1, avg_grid % n_columns
 
@@ -814 +922 @@
           column_top  = avg_grid % cosmo_h(i,j,nlev)
 
-          ! scan through avg_grid column until and set neighbour indices in
-          ! case current_height is either below column_base, in the current
-          ! cell, or above column_top. Keep increasing current cell index until
-          ! the current cell overlaps with the current_height.
+!
+!--       scan through avg_grid column until and set neighbour indices in
+!--       case current_height is either below column_base, in the current
+!--       cell, or above column_top. Keep increasing current cell index until
+!--       the current cell overlaps with the current_height.
           k_intermediate = 1 !avg_grid % cosmo_h is indezed 1-based.
           DO k_palm = 1, avg_grid % nz
 
-             ! Memorize the top and bottom boundaries of the coarse cell and the
-             ! current height within it
+!
+!--          Memorize the top and bottom boundaries of the coarse cell and the
+!--          current height within it
              current_height = avg_grid % z(k_palm) + avg_grid % z0
              h_top    = avg_grid % cosmo_h(i,j,k_intermediate+1)
              h_bottom = avg_grid % cosmo_h(i,j,k_intermediate)
 
-             point_is_above_grid = (current_height > column_top) !22000m, very unlikely
+!
+!--          COSMO column top is located at 22000m, point_is_above_grid is very
+!--          unlikely.
+             point_is_above_grid = (current_height > column_top)
              point_is_below_grid = (current_height < column_base)
 
@@ -835 +948 @@
              )
 
-             ! set default weights
+!
+!--          set default weights
              avg_grid % w(l,k_palm,1:2) = 0.0_dp
 
@@ -852 +966 @@
                 avg_grid % k_min = MAX(k_palm + 1, avg_grid % k_min)
              ELSE
-                ! cycle through intermediate levels until current
-                ! intermediate-grid cell overlaps with current_height
+!
+!--             cycle through intermediate levels until current
+!--             intermediate-grid cell overlaps with current_height
                 DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
                    k_intermediate = k_intermediate + 1
@@ -865 +980 @@
                 END DO
 
-                ! k_intermediate = 48 indicates the last section (indices 48 and 49), i.e.
-                ! k_intermediate = 49 is not the beginning of a valid cell.
+!
+!--             k_intermediate = 48 indicates the last section (indices 48 and 49), i.e.
+!--             k_intermediate = 49 is not the beginning of a valid cell.
                 IF (k_intermediate > nlev-1)  THEN
                    message = "Index " // TRIM(str(k_intermediate)) //          &
@@ -876 +992 @@
                 avg_grid % kkk(l,k_palm,2) = k_intermediate + 1
 
-                ! copmute vertical weights
+!
+!--             compute vertical weights
                 weight = (h_top - current_height) / (h_top - h_bottom)
                 avg_grid % w(l,k_palm,1) = weight
@@ -882 +999 @@
              END IF
 
-          END DO ! k, PALM levels
-       END DO ! l, averaging columns
+!
+!--       Loop over PALM levels k
+          END DO
+
+!
+!--       Loop over averaging columns l
+       END DO
  
     END SUBROUTINE find_vertical_neighbours_and_weights_average
@@ -931 +1053 @@
 !                         ii(i,j,1/2)        ii(i,j,3/4)
 !          
+!------------------------------------------------------------------------------!
     SUBROUTINE compute_horizontal_interp_weights(cosmo_lat, cosmo_lon,         &
        palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)
@@ -950 +1073 @@
        DO i = 0, UBOUND(palm_clon, 1)
       
-          ! weight in lambda direction
+!
+!--       weight in lambda direction
           wl = ( cosmo_lon(palm_ii(i,j,4)) - palm_clon(i,j) ) * cosmo_dxi
 
-          ! weight in phi direction
+!
+!--       weight in phi direction
           wp = ( cosmo_lat(palm_jj(i,j,2)) - palm_clat(i,j) ) * cosmo_dyi
 
@@ -988 +1113 @@
 !> COSMO-DE the velocity points are staggared one half grid spaceing up-grid
 !> which means the first centre point has to be omitted and is set to zero.
+!------------------------------------------------------------------------------!
     SUBROUTINE centre_velocities(u_face, v_face, u_centre, v_centre)
        REAL(dp), DIMENSION(0:,0:,0:), INTENT(IN)  ::  u_face, v_face
@@ -1004 +1130 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Compute the geographical latitude of a point given in rotated-pole cordinates
+!------------------------------------------------------------------------------!
     FUNCTION phirot2phi (phirot, rlarot, polphi, pollam, polgam)
     
@@ -1041 +1172 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Compute the geographical latitude of a point given in rotated-pole cordinates
+!------------------------------------------------------------------------------!
     FUNCTION phi2phirot (phi, rla, polphi, pollam)
     
@@ -1072 +1208 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Compute the geographical longitude of a point given in rotated-pole cordinates
+!------------------------------------------------------------------------------!
     FUNCTION rlarot2rla(phirot, rlarot, polphi, pollam, polgam)
     
@@ -1123 +1264 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Compute the rotated-pole longitude of a point given in geographical cordinates
+!------------------------------------------------------------------------------!
     FUNCTION rla2rlarot ( phi, rla, polphi, pollam, polgam )
 
@@ -1162 +1308 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Rotate the given velocity vector (u,v) from the geographical to the
+!> rotated-pole system
+!------------------------------------------------------------------------------!
     SUBROUTINE uv2uvrot(u, v, rlat, rlon, pollat, pollon, urot, vrot)
     
@@ -1187 +1339 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Rotate the given velocity vector (urot, vrot) from the rotated-pole to the
+!> geographical system
+!------------------------------------------------------------------------------!
     SUBROUTINE uvrot2uv (urot, vrot, rlat, rlon, pollat, pollon, u, v)