!> @file src/transform.f90
!------------------------------------------------------------------------------!
! This file is part of the PALM model system.
!
! PALM is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software
! Foundation, either version 3 of the License, or (at your option) any later
! version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with
! PALM. If not, see
and sets the lowest entry to the !> extrapolated pressure at the surface. !------------------------------------------------------------------------------! SUBROUTINE get_surface_pressure(p, rho, avg_grid) REAL(dp), DIMENSION(:), INTENT(IN) :: rho REAL(dp), DIMENSION(:), INTENT(INOUT) :: p TYPE(grid_definition), INTENT(IN) :: avg_grid REAL(dp) :: drhodz, dz, zk, rhok INTEGER :: k, k_min k_min = avg_grid % k_min zk = avg_grid % z(k_min) rhok = rho(k_min) dz = avg_grid % z(k_min + 1) - avg_grid % z(k_min) drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz p(1) = constant_density_pressure(p(k_min), zk, rhok, drhodz, & 0.0, G) END SUBROUTINE get_surface_pressure FUNCTION constant_density_pressure(pk, zk, rhok, drhodz, z, g) RESULT(p) REAL(dp), INTENT(IN) :: pk, zk, rhok, drhodz, g, z REAL(dp) :: p p = pk + ( zk - z ) * g * ( rhok + 0.5*drhodz * (zk - z) ) END FUNCTION constant_density_pressure !-----------------------------------------------------------------------------! ! Description: ! ----------- !> This routine computes profiles of the two geostrophic wind components ug and !> vg. !-----------------------------------------------------------------------------! SUBROUTINE geostrophic_winds(p_north, p_south, p_east, p_west, rho, f3, & Lx, Ly, phi_n, lam_n, phi_g, lam_g, ug, vg) REAL(dp), DIMENSION(:), INTENT(IN) :: p_north, p_south, p_east, p_west, & rho REAL(dp), INTENT(IN) :: f3, Lx, Ly, phi_n, lam_n, phi_g, lam_g REAL(dp), DIMENSION(:), INTENT(OUT) :: ug, vg REAL(dp) :: facx, facy facx = 1.0_dp / (Lx * f3) facy = 1.0_dp / (Ly * f3) ug(:) = - facy / rho(:) * (p_north(:) - p_south(:)) vg(:) = facx / rho(:) * (p_east(:) - p_west(:)) CALL rotate_vector_field( & ug, vg, angle = meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g)& ) END SUBROUTINE geostrophic_winds !-----------------------------------------------------------------------------! ! Description: ! ----------- !> This routine computes the inverse Plate Carree projection, i.e. in projects !> Cartesian coordinates (x,y) onto a sphere. It returns the latitude and !> lngitude of a geographical system centered at x0 and y0. !-----------------------------------------------------------------------------! SUBROUTINE inv_plate_carree(x, y, x0, y0, r, lat, lon) REAL(dp), INTENT(IN) :: x(:), y(:), x0, y0, r REAL(dp), INTENT(OUT) :: lat(:), lon(:) REAL(dp) :: ri ! TODO check dimensions of lat/lon and y/x match ri = 1.0_dp / r lat(:) = (y(:) - y0) * ri lon(:) = (x(:) - x0) * ri END SUBROUTINE !-----------------------------------------------------------------------------! ! Description: ! ------------ !> Computes the reverse Plate-Carree projection of a x or y position on a !> Cartesian grid. !> !> Input parameters: !> ----------------- !> xy : x or y coordinate of the Cartasian grid point [m]. !> !> xy0 : x or y coordinate that coincides with the origin of the underlying !> sperical system (crossing point of the equator and prime meridian) [m]. !> !> r : Radius of the of the underlying sphere, e.g. EARTH_RADIUS [m]. !> !> Returns: !> -------- !> project : Longitude (in case xy = x) or latitude (xy = y) of the given input !> coordinate xy. !------------------------------------------------------------------------------! ELEMENTAL REAL(dp) FUNCTION project(xy, xy0, r) REAL(dp), INTENT(IN) :: xy, xy0, r 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. ri = 1.0_dp / r project = (xy - xy0) * ri END FUNCTION project REAL(dp) FUNCTION phic_to_phin(phi_c) REAL(dp), INTENT(IN) :: phi_c phic_to_phin = 0.5_dp * PI - ABS(phi_c) END FUNCTION phic_to_phin REAL(dp) FUNCTION lamc_to_lamn(phi_c, lam_c) REAL(dp), INTENT(IN) :: phi_c, lam_c lamc_to_lamn = lam_c IF (phi_c > 0.0_dp) THEN lamc_to_lamn = lam_c - SIGN(PI, lam_c) END IF END FUNCTION lamc_to_lamn 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 gamma_from_hemisphere = PI ELSE gamma_from_hemisphere = 0.0_dp END IF END FUNCTION gamma_from_hemisphere !------------------------------------------------------------------------------! ! Description: ! ------------ !> Computes the geographical coordinates corresponding to the given rotated-pole !> coordinates. !> !> In INIFOR, this routine is used to convert coordinates between two !> rotated-pole systems: COSMO-DE's rotated-pole system, and one centred at the !> PALM-4U domain centre. In this case, the PALM-4U system is thought of as the !> rotated-pole system and the routine is used to rotate back to COSMO-DE's !> system which is thought of as the geographical one. !> !> Input parameters: !> ----------------- !> phir(:), lamr(: ): latitudes and longitudes of the rotated-pole grid !> !> phip, lamp: latitude and longitude of the rotated north pole !> !> gam: "angle between the north poles. If [gam] is not present, the other !> system is the real geographical system." (original phiro2rot !> description) !> !> Output parameters: !> ------------------ !> phi(:,:), lam(:,:): geographical latitudes and logitudes !------------------------------------------------------------------------------! SUBROUTINE rotate_to_cosmo(phir, lamr, phip, lamp, phi, lam, gam) REAL(dp), INTENT(IN) :: phir(0:), lamr(0:), phip, lamp, gam REAL(dp), INTENT(OUT) :: phi(0:,0:), lam(0:,0:) INTEGER :: i, j IF ( SIZE(phi, 1) .NE. SIZE(lam, 1) .OR. & SIZE(phi, 2) .NE. SIZE(lam, 2) ) THEN PRINT *, "inifor: rotate_to_cosmo: Dimensions of phi and lambda do not match. Dimensions are:" PRINT *, "inifor: rotate_to_cosmo: phi: ", SIZE(phi, 1), SIZE(phi, 2) PRINT *, "inifor: rotate_to_cosmo: lam: ", SIZE(lam, 1), SIZE(lam, 2) STOP END IF IF ( SIZE(phir) .NE. SIZE(phi, 2) .OR. & SIZE(lamr) .NE. SIZE(phi, 1) ) THEN PRINT *, "inifor: rotate_to_cosmo: Dimensions of phir and lamr do not match. Dimensions are:" PRINT *, "inifor: rotate_to_cosmo: phir: ", SIZE(phir), SIZE(phi, 2) PRINT *, "inifor: rotate_to_cosmo: lamr: ", SIZE(lamr), SIZE(phi, 1) STOP END IF DO j = 0, UBOUND(phir, 1) DO i = 0, UBOUND(lamr, 1) phi(i,j) = phirot2phi(phir(j) * TO_DEGREES, & lamr(i) * TO_DEGREES, & phip * TO_DEGREES, & lamp * TO_DEGREES, & gam * TO_DEGREES) * TO_RADIANS lam(i,j) = rlarot2rla(phir(j) * TO_DEGREES, & lamr(i) * TO_DEGREES, & phip * TO_DEGREES, & lamp * TO_DEGREES, & gam * TO_DEGREES) * TO_RADIANS END DO END DO END SUBROUTINE rotate_to_cosmo SUBROUTINE rotate_vector_field(x, y, angle) REAL(dp), DIMENSION(:), INTENT(INOUT) :: x, y !< x and y coodrinate in arbitrary units REAL(dp), INTENT(IN) :: angle !< rotation angle [deg] INTEGER :: i REAL(dp) :: sine, cosine, v_rot(2), rotation(2,2) sine = SIN(angle * TO_RADIANS) cosine = COS(angle * TO_RADIANS) ! RESAHPE() fills columns first, so the rotation matrix becomes ! ! rotation = [ cosine -sine ] ! [ sine cosine ] rotation = RESHAPE( (/cosine, sine, -sine, cosine/), (/2, 2/) ) DO i = LBOUND(x, 1), UBOUND(x, 1) v_rot(:) = MATMUL(rotation, (/x(i), y(i)/)) x(i) = v_rot(1) y(i) = v_rot(2) END DO END SUBROUTINE rotate_vector_field !------------------------------------------------------------------------------! ! Description: ! ------------ !> This routine computes the local meridian convergence between a rotated-pole !> and a geographical system using the Eq. (6) given in the DWD manual !> !> Baldauf et al. (2018), Beschreibung des operationelle Kürzestfrist- !> vorhersagemodells COSMO-D2 und COSMO-D2-EPS und seiner Ausgabe in die !> 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) REAL(dp), INTENT(IN) :: phi_n, lam_n, phi_g, lam_g REAL(dp) :: delta delta = atan2( COS(phi_n) * SIN(lam_n - lam_g), & COS(phi_g) * SIN(phi_n) - & SIN(phi_g) * COS(phi_n) * COS(lam_n - lam_g) ) END FUNCTION meridian_convergence_rotated !------------------------------------------------------------------------------! ! Description: ! ------------ !> Compute indices of PALM-4U grid point neighbours in the target !> system (COSMO-DE) by rounding up and down. (i,j) are the indices of !> the PALM-4U grid and (ii(i,j,1-4), jj(i,j,1-4)) contain the indices !> of the its four neigbouring points in the COSMO-DE grid. !> !> !> COSMO-DE grid !> ------------- !> jj, lat !> ^ j !> | \ i !> jj(i,j,2/3) + ... 2 ---\--------/------ 3 !> | | ^ \ / | !> | | |wp \ / | !> | | v \ / | !> latpos + ............ o/ (i,j) | !> | | : | !> | | :<----wl---->| !> jj(i,j,1/4) + ... 1 -------:----------- 4 !> | : : : !> | : : : !> | : lonpos : !> L-----+--------+------------+------> ii, lon !> ii(i,j,1/2) ii(i,j,3/4) !> !> !> Input parameters: !> ----------------- !> source_lat, source_lon : (rotated-pole) coordinates of the source grid (e.g. !> COSMO-DE) !> !> source_dxi, source_dyi : inverse grid spacings of the source grid. !> !> target_lat, target_lon : (rotated-pole) coordinates of the target grid (e.g. !> COSMO-DE) !> !> Output parameters: !> ------------------ !> palm_ii, palm_jj : x and y index arrays of horizontal neighbour columns !> !------------------------------------------------------------------------------! SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon, & palm_clat, palm_clon, & palm_ii, palm_jj) REAL(dp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon REAL(dp), DIMENSION(0:,0:), INTENT(IN) :: palm_clat, palm_clon REAL(dp) :: cosmo_dxi, cosmo_dyi INTEGER, DIMENSION(0:,0:,1:), INTENT(OUT) :: palm_ii, palm_jj REAL(dp) :: lonpos, latpos, lon0, lat0 INTEGER :: i, j lon0 = cosmo_lon(0) lat0 = cosmo_lat(0) cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0)) cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0)) 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. lonpos = (palm_clon(i,j) - lon0) * cosmo_dxi latpos = (palm_clat(i,j) - lat0) * cosmo_dyi IF (lonpos < 0.0 .OR. latpos < 0.0) THEN PRINT *, " Error while finding neighbours: lonpos or latpos out of bounds!" PRINT *, " (i,j) = (", i, ",",j,")" PRINT *, " lonpos ", lonpos*TO_DEGREES, ", latpos ", latpos*TO_DEGREES PRINT *, " lon0 ", lon0 *TO_DEGREES, ", lat0 ", lat0*TO_DEGREES PRINT *, " PALM lon ", palm_clon(i,j)*TO_DEGREES, ", PALM lat ",palm_clat(i,j)*TO_DEGREES STOP END IF palm_ii(i,j,1) = FLOOR(lonpos) palm_ii(i,j,2) = FLOOR(lonpos) palm_ii(i,j,3) = CEILING(lonpos) palm_ii(i,j,4) = CEILING(lonpos) palm_jj(i,j,1) = FLOOR(latpos) palm_jj(i,j,2) = CEILING(latpos) palm_jj(i,j,3) = CEILING(latpos) palm_jj(i,j,4) = FLOOR(latpos) END DO END DO END SUBROUTINE find_horizontal_neighbours SUBROUTINE find_vertical_neighbours_and_weights_interp( palm_grid, & palm_intermediate ) TYPE(grid_definition), INTENT(INOUT) :: palm_grid TYPE(grid_definition), INTENT(IN) :: palm_intermediate INTEGER :: i, j, k, nx, ny, nz, nlev, k_intermediate LOGICAL :: point_is_below_grid, point_is_above_grid, & point_is_in_current_cell REAL(dp) :: current_height, column_base, column_top, h_top, h_bottom, & weight nx = palm_grid % nx ny = palm_grid % ny nz = palm_grid % nz nlev = palm_intermediate % nz ! in each column of the fine grid, find vertical neighbours of every cell DO j = 0, ny DO i = 0, nx k_intermediate = 0 column_base = palm_intermediate % h(i,j,0) 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. DO k = 1, nz ! 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) h_bottom = palm_intermediate % h(i,j,k_intermediate) point_is_above_grid = (current_height > column_top) !22000m, very unlikely point_is_below_grid = (current_height < column_base) point_is_in_current_cell = ( & current_height >= h_bottom .AND. & current_height < h_top & ) ! set default weights palm_grid % w_verti(i,j,k,1:2) = 0.0_dp IF (point_is_above_grid) THEN palm_grid % kk(i,j,k,1:2) = nlev palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp message = "PALM-4U grid extends above COSMO-DE model top." CALL abort('find_vertical_neighbours_and_weights', message) ELSE IF (point_is_below_grid) THEN palm_grid % kk(i,j,k,1:2) = 0 palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp ELSE ! 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 h_top = palm_intermediate % h(i,j,k_intermediate+1) h_bottom = palm_intermediate % h(i,j,k_intermediate) point_is_in_current_cell = ( & current_height >= h_bottom .AND. & current_height < h_top & ) END DO IF (k_intermediate > nlev-1) THEN message = "Index " // TRIM(str(k_intermediate)) // & " is above intermediate grid range." CALL abort('find_vertical_neighbours', message) END IF palm_grid % kk(i,j,k,1) = k_intermediate palm_grid % kk(i,j,k,2) = k_intermediate + 1 ! copmute vertical weights weight = (h_top - current_height) / (h_top - h_bottom) palm_grid % w_verti(i,j,k,1) = weight palm_grid % w_verti(i,j,k,2) = 1.0_dp - weight END IF END DO END DO END DO END SUBROUTINE find_vertical_neighbours_and_weights_interp SUBROUTINE find_vertical_neighbours_and_weights_average( avg_grid ) TYPE(grid_definition), INTENT(INOUT) :: avg_grid INTEGER :: i, j, k_palm, k_intermediate, l, nlev LOGICAL :: point_is_below_grid, point_is_above_grid, & point_is_in_current_cell REAL(dp) :: current_height, column_base, column_top, h_top, h_bottom, & weight avg_grid % k_min = LBOUND(avg_grid % z, 1) nlev = SIZE(avg_grid % cosmo_h, 3) ! in each column of the fine grid, find vertical neighbours of every cell DO l = 1, avg_grid % n_columns i = avg_grid % iii(l) j = avg_grid % jjj(l) column_base = avg_grid % cosmo_h(i,j,1) 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. 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 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 point_is_below_grid = (current_height < column_base) point_is_in_current_cell = ( & current_height >= h_bottom .AND. & current_height < h_top & ) ! set default weights avg_grid % w(l,k_palm,1:2) = 0.0_dp IF (point_is_above_grid) THEN avg_grid % kkk(l,k_palm,1:2) = nlev avg_grid % w(l,k_palm,1:2) = - 2.0_dp message = "PALM-4U grid extends above COSMO-DE model top." CALL abort('find_vertical_neighbours_and_weights_average', message) ELSE IF (point_is_below_grid) THEN avg_grid % kkk(l,k_palm,1:2) = 0 avg_grid % w(l,k_palm,1:2) = - 2.0_dp 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 DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1) k_intermediate = k_intermediate + 1 h_top = avg_grid % cosmo_h(i,j,k_intermediate+1) h_bottom = avg_grid % cosmo_h(i,j,k_intermediate) point_is_in_current_cell = ( & current_height >= h_bottom .AND. & current_height < h_top & ) 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. IF (k_intermediate > nlev-1) THEN message = "Index " // TRIM(str(k_intermediate)) // & " is above intermediate grid range." CALL abort('find_vertical_neighbours', message) END IF avg_grid % kkk(l,k_palm,1) = k_intermediate avg_grid % kkk(l,k_palm,2) = k_intermediate + 1 ! copmute vertical weights weight = (h_top - current_height) / (h_top - h_bottom) avg_grid % w(l,k_palm,1) = weight avg_grid % w(l,k_palm,2) = 1.0_dp - weight END IF END DO ! k, PALM levels END DO ! l, averaging columns END SUBROUTINE find_vertical_neighbours_and_weights_average !------------------------------------------------------------------------------! ! Description: ! ------------ !> Compute the four weights for horizontal bilinear interpolation given the !> coordinates clon(i,j) clat(i,j) of the PALM-4U grid in the COSMO-DE !> rotated-pole grid and the neightbour indices ii(i,j,1-4) and jj(i,j,1-4). !> !> Input parameters: !> ----------------- !> palm_grid % clon : longitudes of PALM-4U scalars (cell centres) in COSMO-DE's rotated-pole grid [rad] !> !> palm_grid % clat : latitudes of PALM-4U cell centres in COSMO-DE's rotated-pole grid [rad] !> !> cosmo_grid % lon : rotated-pole longitudes of scalars (cell centres) of the COSMO-DE grid [rad] !> !> cosmo_grid % lat : rotated-pole latitudes of scalars (cell centers) of the COSMO-DE grid [rad] !> !> cosmo_grid % dxi : inverse grid spacing in the first dimension [m^-1] !> !> cosmo_grid % dyi : inverse grid spacing in the second dimension [m^-1] !> !> Output parameters: !> ------------------ !> palm_grid % w_horiz(:,:,1-4) : weights for bilinear horizontal interpolation ! ! COSMO-DE grid ! ------------- ! jj, lat ! ^ j ! | \ i ! jj(i,j,2/3) + ... 2 ---\--------/------ 3 ! | | ^ \ / | ! | | |wp \ / | ! | | v \ / | ! latpos + ............ o/ (i,j) | ! | | : | ! | | :<----wl---->| ! jj(i,j,1/4) + ... 1 -------:----------- 4 ! | : : : ! | : : : ! | : lonpos : ! L-----+--------+------------+------> ii, lon ! 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) REAL(dp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon REAL(dp) :: cosmo_dxi, cosmo_dyi REAL(dp), DIMENSION(0:,0:), INTENT(IN) :: palm_clat, palm_clon INTEGER, DIMENSION(0:,0:,1:), INTENT(IN) :: palm_ii, palm_jj REAL(dp), DIMENSION(0:,0:,1:), INTENT(OUT) :: palm_w_horiz REAL(dp) :: wl, wp INTEGER :: i, j cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0)) cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0)) DO j = 0, UBOUND(palm_clon, 2) DO i = 0, UBOUND(palm_clon, 1) ! weight in lambda direction wl = ( cosmo_lon(palm_ii(i,j,4)) - palm_clon(i,j) ) * cosmo_dxi ! weight in phi direction wp = ( cosmo_lat(palm_jj(i,j,2)) - palm_clat(i,j) ) * cosmo_dyi IF (wl > 1.0_dp .OR. wl < 0.0_dp) THEN message = "Horizontal weight wl = " // TRIM(real_to_str(wl)) // & " is out bounds." CALL abort('compute_horizontal_interp_weights', message) END IF IF (wp > 1.0_dp .OR. wp < 0.0_dp) THEN message = "Horizontal weight wp = " // TRIM(real_to_str(wp)) // & " is out bounds." CALL abort('compute_horizontal_interp_weights', message) END IF palm_w_horiz(i,j,1) = wl * wp palm_w_horiz(i,j,2) = wl * (1.0_dp - wp) palm_w_horiz(i,j,3) = (1.0_dp - wl) * (1.0_dp - wp) palm_w_horiz(i,j,4) = 1.0_dp - SUM( palm_w_horiz(i,j,1:3) ) END DO END DO END SUBROUTINE compute_horizontal_interp_weights !------------------------------------------------------------------------------! ! Description: ! ------------ !> Interpolates u and v components of velocities located at cell faces to the !> cell centres by averaging neighbouring values. !> !> This routine is designed to be used with COSMO-DE arrays where there are the !> same number of grid points for scalars (centres) and velocities (faces). In !> 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 REAL(dp), DIMENSION(0:,0:,0:), INTENT(OUT) :: u_centre, v_centre INTEGER :: nx, ny nx = UBOUND(u_face, 1) ny = UBOUND(u_face, 2) u_centre(0,:,:) = 0.0_dp u_centre(1:,:,:) = 0.5_dp * ( u_face(0:nx-1,:,:) + u_face(1:,:,:) ) v_centre(:,0,:) = 0.0_dp v_centre(:,1:,:) = 0.5_dp * ( v_face(:,0:ny-1,:) + v_face(:,1:,:) ) END SUBROUTINE centre_velocities FUNCTION phirot2phi (phirot, rlarot, polphi, pollam, polgam) REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole REAL(dp), INTENT (IN) :: pollam !< longitude of the rotated north pole REAL(dp), INTENT (IN) :: phirot !< latitude in the rotated system REAL(dp), INTENT (IN) :: rlarot !< longitude in the rotated system REAL(dp), INTENT (IN) :: polgam !< angle between the north poles of the systems REAL(dp) :: phirot2phi !< latitude in the geographical system REAL(dp) :: zsinpol, zcospol, zphis, zrlas, zarg, zgam zsinpol = SIN(polphi * TO_RADIANS) zcospol = COS(polphi * TO_RADIANS) zphis = phirot * TO_RADIANS IF (rlarot > 180.0_dp) THEN zrlas = rlarot - 360.0_dp ELSE zrlas = rlarot END IF zrlas = zrlas * TO_RADIANS IF (polgam /= 0.0_dp) THEN zgam = polgam * TO_RADIANS zarg = zsinpol * SIN (zphis) + & zcospol * COS(zphis) * ( COS(zrlas) * COS(zgam) - & SIN(zgam) * SIN(zrlas) ) ELSE zarg = zcospol * COS (zphis) * COS (zrlas) + zsinpol * SIN (zphis) END IF phirot2phi = ASIN (zarg) * TO_DEGREES END FUNCTION phirot2phi FUNCTION phi2phirot (phi, rla, polphi, pollam) REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole REAL(dp), INTENT (IN) :: pollam !< longitude of the rotated north pole REAL(dp), INTENT (IN) :: phi !< latitude in the geographical system REAL(dp), INTENT (IN) :: rla !< longitude in the geographical system REAL(dp) :: phi2phirot !< longitude in the rotated system REAL(dp) :: zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1 zsinpol = SIN(polphi * TO_RADIANS) zcospol = COS(polphi * TO_RADIANS) zlampol = pollam * TO_RADIANS zphi = phi * TO_RADIANS IF (rla > 180.0_dp) THEN zrla1 = rla - 360.0_dp ELSE zrla1 = rla END IF zrla = zrla1 * TO_RADIANS zarg1 = SIN(zphi) * zsinpol zarg2 = COS(zphi) * zcospol * COS(zrla - zlampol) phi2phirot = ASIN(zarg1 + zarg2) * TO_DEGREES END FUNCTION phi2phirot FUNCTION rlarot2rla(phirot, rlarot, polphi, pollam, polgam) REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole REAL(dp), INTENT (IN) :: pollam !< longitude of the rotated north pole REAL(dp), INTENT (IN) :: phirot !< latitude in the rotated system REAL(dp), INTENT (IN) :: rlarot !< longitude in the rotated system REAL(dp), INTENT (IN) :: polgam !< angle between the north poles of the systems REAL(dp) :: rlarot2rla !< latitude in the geographical system REAL(dp) :: zsinpol, zcospol, zlampol, zphis, zrlas, zarg1, zarg2, zgam zsinpol = SIN(TO_RADIANS * polphi) zcospol = COS(TO_RADIANS * polphi) zlampol = TO_RADIANS * pollam zphis = TO_RADIANS * phirot IF (rlarot > 180.0_dp) THEN zrlas = rlarot - 360.0_dp ELSE zrlas = rlarot END IF zrlas = TO_RADIANS * zrlas IF (polgam /= 0.0_dp) THEN zgam = TO_RADIANS * polgam zarg1 = SIN(zlampol) * (zcospol * SIN(zphis) - zsinpol*COS(zphis) * & (COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) - & COS(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) + & COS(zrlas) * SIN(zgam) ) zarg2 = COS (zlampol) * (zcospol * SIN(zphis) - zsinpol*COS(zphis) * & (COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) + & SIN(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) + & COS(zrlas) * SIN(zgam) ) ELSE zarg1 = SIN (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis) + & zcospol * SIN(zphis)) - & COS (zlampol) * SIN(zrlas) * COS(zphis) zarg2 = COS (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis) + & zcospol * SIN(zphis)) + & SIN (zlampol) * SIN(zrlas) * COS(zphis) END IF IF (zarg2 == 0.0_dp) zarg2 = 1.0E-20_dp rlarot2rla = ATAN2(zarg1,zarg2) * TO_DEGREES END FUNCTION rlarot2rla FUNCTION rla2rlarot ( phi, rla, polphi, pollam, polgam ) REAL(dp), INTENT (IN) :: polphi !< latitude of the rotated north pole REAL(dp), INTENT (IN) :: pollam !< longitude of the rotated north pole REAL(dp), INTENT (IN) :: phi !< latitude in geographical system REAL(dp), INTENT (IN) :: rla !< longitude in geographical system REAL(dp), INTENT (IN) :: polgam !< angle between the north poles of the systems REAL (KIND=dp) :: rla2rlarot !< latitude in the the rotated system REAL (KIND=dp) :: zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1 zsinpol = SIN(polphi * TO_RADIANS) zcospol = COS(polphi * TO_RADIANS) zlampol = pollam * TO_RADIANS zphi = phi * TO_RADIANS IF (rla > 180.0_dp) THEN zrla1 = rla - 360.0_dp ELSE zrla1 = rla END IF zrla = zrla1 * TO_RADIANS zarg1 = - SIN (zrla-zlampol) * COS(zphi) zarg2 = - zsinpol * COS(zphi) * COS(zrla-zlampol) + zcospol * SIN(zphi) IF (zarg2 == 0.0_dp) zarg2 = 1.0E-20_dp rla2rlarot = ATAN2 (zarg1,zarg2) * TO_DEGREES IF (polgam /= 0.0_dp ) THEN rla2rlarot = polgam + rla2rlarot IF (rla2rlarot > 180._dp) rla2rlarot = rla2rlarot - 360.0_dp END IF END FUNCTION rla2rlarot SUBROUTINE uv2uvrot(u, v, rlat, rlon, pollat, pollon, urot, vrot) REAL(dp), INTENT (IN) :: u, v !< wind components in the true geographical system REAL(dp), INTENT (IN) :: rlat, rlon !< coordinates in the true geographical system REAL(dp), INTENT (IN) :: pollat, pollon !< latitude and longitude of the north pole of the rotated grid REAL(dp), INTENT (OUT) :: urot, vrot !< wind components in the rotated grid REAL (dp) :: zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm zsinpol = SIN(pollat * TO_RADIANS) zcospol = COS(pollat * TO_RADIANS) zlonp = (pollon-rlon) * TO_RADIANS zlat = rlat * TO_RADIANS zarg1 = zcospol * SIN(zlonp) zarg2 = zsinpol * COS(zlat) - zcospol * SIN(zlat) * COS(zlonp) znorm = 1.0_dp / SQRT(zarg1*zarg1 + zarg2*zarg2) urot = u * zarg2 * znorm - v * zarg1 * znorm vrot = u * zarg1 * znorm + v * zarg2 * znorm END SUBROUTINE uv2uvrot SUBROUTINE uvrot2uv (urot, vrot, rlat, rlon, pollat, pollon, u, v) REAL(dp), INTENT(IN) :: urot, vrot !< wind components in the rotated grid REAL(dp), INTENT(IN) :: rlat, rlon !< latitude and longitude in the true geographical system REAL(dp), INTENT(IN) :: pollat, pollon !< latitude and longitude of the north pole of the rotated grid REAL(dp), INTENT(OUT) :: u, v !< wind components in the true geographical system REAL(dp) :: zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm zsinpol = SIN(pollat * TO_RADIANS) zcospol = COS(pollat * TO_RADIANS) zlonp = (pollon-rlon) * TO_RADIANS zlat = rlat * TO_RADIANS zarg1 = zcospol * SIN(zlonp) zarg2 = zsinpol * COS(zlat) - zcospol * SIN(zlat) * COS(zlonp) znorm = 1.0_dp / SQRT(zarg1*zarg1 + zarg2*zarg2) u = urot * zarg2 * znorm + vrot * zarg1 * znorm v = - urot * zarg1 * znorm + vrot * zarg2 * znorm END SUBROUTINE uvrot2uv END MODULE