source: palm/trunk/UTIL/inifor/src/inifor_transform.f90 @ 3557

Last change on this file since 3557 was 3557, checked in by eckhard, 3 years ago

inifor: Updated documentation

  • Property svn:keywords set to Id
File size: 52.6 KB
Line 
1!> @file src/inifor_transform.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2017-2018 Leibniz Universitaet Hannover
18! Copyright 2017-2018 Deutscher Wetterdienst Offenbach
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: inifor_transform.f90 3557 2018-11-22 16:01:22Z eckhard $
28! Updated documentation
29!
30!
31! 3537 2018-11-20 10:53:14Z eckhard
32! bugfix: working precision added
33!
34! 3447 2018-10-29 15:52:54Z eckhard
35! Renamed source files for compatibilty with PALM build system
36!
37!
38! 3395 2018-10-22 17:32:49Z eckhard
39! Switched addressing of averaging regions from index bounds to list of columns
40! Added routines for the computation of geostrophic winds including:
41!    - the downward extrapolation of density (linear) and
42!    - pressure (hydrostatic equation) and
43!    - rotation of geostrophic wind components to PALM frame of reference
44!
45! 3183 2018-07-27 14:25:55Z suehring
46! Introduced new PALM grid stretching
47! Removed unnecessary subroutine parameters
48! Renamed kcur to k_intermediate
49!
50!
51! 3182 2018-07-27 13:36:03Z suehring
52! Initial revision
53!
54!
55!
56! Authors:
57! --------
58!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
59!
60! Description:
61! ------------
62!> The transform module provides INIFOR's low-level transformation and
63!> interpolation routines. The rotated-pole transformation routines phirot2phi,
64!> phi2phirot, rlarot2rla, rla2rlarot, uv2uvrot, and uvrot2uv are adapted from
65!> int2lm's utility routines.
66!------------------------------------------------------------------------------!
67 MODULE transform
68
69    USE control
70    USE defs,                                                                  &
71        ONLY: G, TO_DEGREES, TO_RADIANS, PI, dp
72    USE types
73    USE util,                                                                  &       
74        ONLY: real_to_str, str
75
76    IMPLICIT NONE
77
78 CONTAINS
79
80!------------------------------------------------------------------------------!
81! Description:
82! ------------
83!> Interpolates linearly in the vertical direction in very column (i,j) of the
84!> output array outvar(i,j,:) using values of the source array invar. In cells
85!> that are outside the COSMO-DE domain, indicated by negative interpolation
86!> weights, extrapolate constantly from the cell above.
87!>
88!> Input parameters:
89!> -----------------
90!> invar : Array of source data
91!>
92!> outgrid % kk : Array of vertical neighbour indices. kk(i,j,k,:) contain the
93!>     indices of the two vertical neighbors of PALM-4U point (i,j,k) on the
94!>     input grid corresponding to the source data invar.
95!>
96!> outgrid % w_verti : Array of weights for vertical linear interpolation
97!>     corresponding to neighbour points indexed by kk.
98!>
99!> Output papameters:
100!> ------------------
101!> outvar : Array of interpolated data
102!------------------------------------------------------------------------------!
103    SUBROUTINE interpolate_1d(in_arr, out_arr, outgrid)
104       TYPE(grid_definition), INTENT(IN) ::  outgrid
105       REAL(dp), INTENT(IN)              ::  in_arr(0:,0:,0:)
106       REAL(dp), INTENT(OUT)             ::  out_arr(0:,0:,:)
107
108       INTEGER :: i, j, k, l, nz
109
110       nz = UBOUND(out_arr, 3)
111
112       DO j = LBOUND(out_arr, 2), UBOUND(out_arr, 2)
113       DO i = LBOUND(out_arr, 1), UBOUND(out_arr, 1)
114       DO k = nz, LBOUND(out_arr, 3), -1
115
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.
121          IF (outgrid % w_verti(i,j,k,1) < -1.0_dp .AND. k < nz)  THEN
122             out_arr(i,j,k) = out_arr(i,j,k+1)
123          ELSE
124             out_arr(i,j,k) = 0.0_dp
125             DO l = 1, 2
126                out_arr(i,j,k) = out_arr(i,j,k) +                              &
127                    outgrid % w_verti(i,j,k,l) *                               &
128                    in_arr(i,j,outgrid % kk(i,j,k, l) )
129             END DO
130          END IF
131       END DO
132       END DO
133       END DO
134    END SUBROUTINE interpolate_1d
135
136
137!------------------------------------------------------------------------------!
138! Description:
139! ------------
140!> Interpolates bi-linearly in horizontal planes on every k level of the output
141!> array outvar(:,:,k) using values of the source array invar(:,:,:). The source
142!> (invar) and interpolation array (outvar) need to have matching dimensions.
143!>
144!> Input parameters:
145!> -----------------
146!> invar : Array of source data
147!>
148!> outgrid % ii, % jj : Array of neighbour indices in x and y direction.
149!>     ii(i,j,k,:), and jj(i,j,k,:) contain the four horizontal neighbour points
150!>     of PALM-4U point (i,j,k) on the input grid corresponding to the source
151!>     data invar. (The outgrid carries the relationship with the ingrid in the
152!      form of the interpoaltion weights.)
153!>
154!> outgrid % w_horiz: Array of weights for horizontal bi-linear interpolation
155!>     corresponding to neighbour points indexed by ii and jj.
156!>
157!> Output papameters:
158!> ------------------
159!> outvar : Array of interpolated data
160!------------------------------------------------------------------------------!
161    SUBROUTINE interpolate_2d(invar, outvar, outgrid, ncvar)
162!
163!--    I index 0-based for the indices of the outvar to be consistent with the
164!--    outgrid indices and interpolation weights.
165       TYPE(grid_definition), INTENT(IN)  ::  outgrid
166       REAL(dp), INTENT(IN)               ::  invar(0:,0:,0:)
167       REAL(dp), INTENT(OUT)              ::  outvar(0:,0:,0:)
168       TYPE(nc_var), INTENT(IN), OPTIONAL ::  ncvar
169
170       INTEGER ::  i, j, k, l
171
172!
173!--    TODO: check if input dimensions are consistent, i.e. ranges are correct
174       IF ( UBOUND(outvar, 3) .GT. UBOUND(invar, 3) )  THEN
175           message = "Output array for '" // TRIM(ncvar % name) // "' has ' more levels (" // &
176              TRIM(str(UBOUND(outvar, 3))) // ") than input variable ("//&
177              TRIM(str(UBOUND(invar, 3))) // ")."
178           CALL abort('interpolate_2d', message) 
179       END IF
180
181       DO k = 0, UBOUND(outvar, 3)
182       DO j = 0, UBOUND(outvar, 2)
183       DO i = 0, UBOUND(outvar, 1)
184          outvar(i,j,k) = 0.0_dp
185          DO l = 1, 4
186             
187             outvar(i,j,k) = outvar(i,j,k) +                                   &
188                outgrid % w_horiz(i,j,l) * invar( outgrid % ii(i,j,l),         & 
189                                                  outgrid % jj(i,j,l),         &
190                                                  k )
191          END DO
192       END DO
193       END DO
194       END DO
195       
196    END SUBROUTINE interpolate_2d
197
198
199!------------------------------------------------------------------------------!
200! Description:
201! ------------
202!> Compute the horizontal average of the in_arr(:,:,:) and store it in
203!> out_arr(:)
204!------------------------------------------------------------------------------!
205    SUBROUTINE average_2d(in_arr, out_arr, ii, jj)
206       REAL(dp), INTENT(IN)              ::  in_arr(0:,0:,0:)
207       REAL(dp), INTENT(OUT)             ::  out_arr(0:)
208       INTEGER, INTENT(IN), DIMENSION(:) ::  ii, jj
209
210       INTEGER  ::  i, j, k, l
211       REAL(dp) ::  ni
212
213       IF (SIZE(ii) .NE. SIZE(jj))  THEN
214          message = "Length of 'ii' and 'jj' index lists do not match." //     &
215             NEW_LINE(' ') // "ii has " // str(SIZE(ii)) // " elements, " //   &
216             NEW_LINE(' ') // "jj has " // str(SIZE(jj)) // "."
217          CALL abort('average_2d', message)
218       END IF
219
220       DO k = 0, UBOUND(out_arr, 1)
221
222          out_arr(k) = 0.0_dp
223          DO l = 1, UBOUND(ii, 1)
224             i = ii(l)
225             j = jj(l)
226             out_arr(k) = out_arr(k) + in_arr(i, j, k)
227          END DO
228
229       END DO
230
231       ni = 1.0_dp / SIZE(ii)
232       out_arr(:) = out_arr(:) * ni
233
234    END SUBROUTINE average_2d
235
236
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!------------------------------------------------------------------------------!
248    SUBROUTINE interpolate_3d(source_array, palm_array, palm_intermediate, palm_grid)
249       TYPE(grid_definition), INTENT(IN) ::  palm_intermediate, palm_grid
250       REAL(dp), DIMENSION(:,:,:), INTENT(IN)  ::  source_array
251       REAL(dp), DIMENSION(:,:,:), INTENT(OUT) ::  palm_array
252       REAL(dp), DIMENSION(:,:,:), ALLOCATABLE ::  intermediate_array
253       INTEGER ::  nx, ny, nlev
254
255       nx = palm_intermediate % nx
256       ny = palm_intermediate % ny
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.
263       ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nlev-1)) !
264
265       CALL interpolate_2d(source_array, intermediate_array, palm_intermediate)
266
267!
268!--    Interpolate from intermediate grid to palm_grid grid, includes
269!--    extrapolation for cells below COSMO domain.
270       CALL interpolate_1d(intermediate_array, palm_array, palm_grid)
271
272       DEALLOCATE(intermediate_array)
273
274    END SUBROUTINE interpolate_3d
275
276
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!------------------------------------------------------------------------------!
283    SUBROUTINE average_profile(source_array, profile_array, avg_grid)
284       TYPE(grid_definition), INTENT(IN)          ::  avg_grid
285       REAL(dp), DIMENSION(:,:,:), INTENT(IN)     ::  source_array
286       REAL(dp), DIMENSION(:), INTENT(OUT)        ::  profile_array
287
288       INTEGER ::  i_source, j_source, k_profile, k_source, l, m
289
290       REAL ::  ni_columns
291
292       profile_array(:) = 0.0_dp
293
294       DO l = 1, avg_grid % n_columns
295          i_source = avg_grid % iii(l)
296          j_source = avg_grid % jjj(l)
297
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
305
306                k_source = avg_grid % kkk(l, k_profile, m)
307
308                profile_array(k_profile) = profile_array(k_profile)            &
309                   + avg_grid % w(l, k_profile, m)                             &
310                   * source_array(i_source, j_source, k_source)
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
322
323       ni_columns = 1.0_dp / avg_grid % n_columns
324       profile_array(:) = profile_array(:) * ni_columns
325
326!
327!--    Constant extrapolation to the bottom
328       profile_array(1:avg_grid % k_min-1) = profile_array(avg_grid % k_min)
329
330    END SUBROUTINE average_profile
331
332
333!------------------------------------------------------------------------------!
334! Description:
335! ------------
336!> Extrapolates density linearly from the level 'k_min' downwards.
337!------------------------------------------------------------------------------!
338    SUBROUTINE extrapolate_density(rho, avg_grid)
339       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  rho
340       TYPE(grid_definition), INTENT(IN)     ::  avg_grid
341
342       REAL(dp) ::  drhodz, dz, zk, rhok
343       INTEGER  ::  k_min
344
345       k_min  = avg_grid % k_min
346       zk     = avg_grid % z(k_min)
347       rhok   = rho(k_min)
348       dz     = avg_grid % z(k_min + 1) - avg_grid % z(k_min)
349       drhodz = (rho(k_min + 1) - rho(k_min)) / dz
350
351       rho(1:k_min-1) = rhok + drhodz * (avg_grid % z(1:k_min-1) - zk)
352
353    END SUBROUTINE extrapolate_density
354
355
356!------------------------------------------------------------------------------!
357! Description:
358! ------------
359!> Driver for extrapolating pressure from PALM level k_min downwards
360!------------------------------------------------------------------------------!
361    SUBROUTINE extrapolate_pressure(p, rho, avg_grid)
362       REAL(dp), DIMENSION(:), INTENT(IN)    ::  rho
363       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  p
364       TYPE(grid_definition), INTENT(IN)     ::  avg_grid
365
366       REAL(dp) ::  drhodz, dz, zk, rhok
367       INTEGER  ::  k, k_min
368
369       k_min = avg_grid % k_min
370       zk    = avg_grid % z(k_min)
371       rhok  = rho(k_min)
372       dz    = avg_grid % z(k_min + 1) - avg_grid % z(k_min)
373       drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz
374
375       DO k = 1, k_min-1
376          p(k) = constant_density_pressure(p(k_min), zk, rhok, drhodz,         &
377                                           avg_grid % z(k), G)
378       END DO
379
380    END SUBROUTINE extrapolate_pressure
381
382
383!------------------------------------------------------------------------------!
384! Description:
385! ------------
386!> Takes the averaged pressure profile <p> and sets the lowest entry to the
387!> extrapolated pressure at the surface.
388!------------------------------------------------------------------------------!
389    SUBROUTINE get_surface_pressure(p, rho, avg_grid)
390       REAL(dp), DIMENSION(:), INTENT(IN)    ::  rho
391       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  p
392       TYPE(grid_definition), INTENT(IN)     ::  avg_grid
393
394       REAL(dp) ::  drhodz, dz, zk, rhok
395       INTEGER  ::  k, k_min
396
397       k_min = avg_grid % k_min
398       zk    = avg_grid % z(k_min)
399       rhok  = rho(k_min)
400       dz    = avg_grid % z(k_min + 1) - avg_grid % z(k_min)
401       drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz
402
403       p(1) = constant_density_pressure(p(k_min), zk, rhok, drhodz,            &
404                                        0.0_dp, G)
405
406    END SUBROUTINE get_surface_pressure
407
408
409    FUNCTION constant_density_pressure(pk, zk, rhok, drhodz, z, g)  RESULT(p)
410
411       REAL(dp), INTENT(IN)  ::  pk, zk, rhok, drhodz, g, z
412       REAL(dp) ::  p
413
414       p = pk + ( zk - z ) * g * ( rhok + 0.5*drhodz * (zk - z) )
415
416    END FUNCTION constant_density_pressure
417
418!-----------------------------------------------------------------------------!
419! Description:
420! -----------
421!> This routine computes profiles of the two geostrophic wind components ug and
422!> vg.
423!-----------------------------------------------------------------------------!
424    SUBROUTINE geostrophic_winds(p_north, p_south, p_east, p_west, rho, f3,    &
425                                 Lx, Ly, phi_n, lam_n, phi_g, lam_g, ug, vg)
426
427    REAL(dp), DIMENSION(:), INTENT(IN)  ::  p_north, p_south, p_east, p_west,  &
428                                            rho
429    REAL(dp), INTENT(IN)                ::  f3, Lx, Ly, phi_n, lam_n, phi_g, lam_g
430    REAL(dp), DIMENSION(:), INTENT(OUT) ::  ug, vg
431    REAL(dp)                            ::  facx, facy
432
433    facx = 1.0_dp / (Lx * f3)
434    facy = 1.0_dp / (Ly * f3)
435    ug(:) = - facy / rho(:) * (p_north(:) - p_south(:))
436    vg(:) =   facx / rho(:) * (p_east(:) - p_west(:))
437
438    CALL rotate_vector_field(                                                  &
439       ug, vg, angle = meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g)&
440    )
441
442    END SUBROUTINE geostrophic_winds
443
444
445!-----------------------------------------------------------------------------!
446! Description:
447! -----------
448!> This routine computes the inverse Plate Carree projection, i.e. in projects
449!> Cartesian coordinates (x,y) onto a sphere. It returns the latitude and
450!> lngitude of a geographical system centered at x0 and y0.
451!-----------------------------------------------------------------------------!
452    SUBROUTINE inv_plate_carree(x, y, x0, y0, r, lat, lon)
453       REAL(dp), INTENT(IN)  ::  x(:), y(:), x0, y0, r
454       REAL(dp), INTENT(OUT) ::  lat(:), lon(:)
455       
456       REAL(dp) :: ri
457
458!
459!--    TODO check dimensions of lat/lon and y/x match
460
461       ri = 1.0_dp / r
462       
463       lat(:) = (y(:) - y0) * ri
464       lon(:) = (x(:) - x0) * ri
465    END SUBROUTINE 
466
467
468!-----------------------------------------------------------------------------!
469! Description:
470! ------------
471!> Computes the reverse Plate-Carree projection of a x or y position on a
472!> Cartesian grid.
473!>
474!> Input parameters:
475!> -----------------
476!> xy : x or y coordinate of the Cartasian grid point [m].
477!>
478!> xy0 : x or y coordinate that coincides with the origin of the underlying
479!>     sperical system (crossing point of the equator and prime meridian) [m].
480!>
481!> r : Radius of the of the underlying sphere, e.g. EARTH_RADIUS [m].
482!>
483!> Returns:
484!> --------
485!> project : Longitude (in case xy = x) or latitude (xy = y) of the given input
486!>     coordinate xy.
487!------------------------------------------------------------------------------!
488    ELEMENTAL REAL(dp) FUNCTION project(xy, xy0, r)
489       REAL(dp), INTENT(IN)  ::  xy, xy0, r
490       REAL(dp) :: ri
491
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.
496       ri = 1.0_dp / r
497
498       project = (xy - xy0) * ri
499
500    END FUNCTION project
501
502
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!------------------------------------------------------------------------------!
509    REAL(dp) FUNCTION phic_to_phin(phi_c)
510        REAL(dp), INTENT(IN) ::  phi_c
511
512        phic_to_phin = 0.5_dp * PI - ABS(phi_c)
513
514    END FUNCTION phic_to_phin
515
516   
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!------------------------------------------------------------------------------!
524    REAL(dp) FUNCTION lamc_to_lamn(phi_c, lam_c)
525       REAL(dp), INTENT(IN) ::  phi_c, lam_c
526       
527       lamc_to_lamn = lam_c
528       IF (phi_c > 0.0_dp)  THEN
529          lamc_to_lamn = lam_c - SIGN(PI, lam_c)
530       END IF
531
532    END FUNCTION lamc_to_lamn
533
534
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!------------------------------------------------------------------------------!
543    REAL(dp) FUNCTION gamma_from_hemisphere(phi_cg, phi_ref)
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
552           gamma_from_hemisphere = PI
553       ELSE
554           gamma_from_hemisphere = 0.0_dp
555       END IF
556    END FUNCTION gamma_from_hemisphere
557
558
559!------------------------------------------------------------------------------!
560! Description:
561! ------------
562!> Computes the geographical coordinates corresponding to the given rotated-pole
563!> coordinates.
564!>
565!> In INIFOR, this routine is used to convert coordinates between two
566!> rotated-pole systems: COSMO-DE's rotated-pole system, and one centred at the
567!> PALM-4U domain centre. In this case, the PALM-4U system is thought of as the
568!> rotated-pole system and the routine is used to rotate back to COSMO-DE's
569!> system which is thought of as the geographical one.
570!>
571!> Input parameters:
572!> -----------------
573!> phir(:), lamr(: ): latitudes and longitudes of the rotated-pole grid
574!>
575!> phip, lamp: latitude and longitude of the rotated north pole
576!>
577!> gam: "angle between the north poles. If [gam] is not present, the other
578!>       system is the real geographical system." (original phiro2rot
579!>       description)
580!>
581!> Output parameters:
582!> ------------------
583!> phi(:,:), lam(:,:): geographical latitudes and logitudes
584!------------------------------------------------------------------------------!
585    SUBROUTINE rotate_to_cosmo(phir, lamr, phip, lamp, phi, lam, gam)
586       REAL(dp), INTENT(IN)  ::  phir(0:), lamr(0:), phip, lamp, gam
587       REAL(dp), INTENT(OUT) ::  phi(0:,0:), lam(0:,0:)
588
589       INTEGER ::  i, j
590       
591       IF ( SIZE(phi, 1) .NE. SIZE(lam, 1) .OR. &
592            SIZE(phi, 2) .NE. SIZE(lam, 2) )  THEN
593          PRINT *, "inifor: rotate_to_cosmo: Dimensions of phi and lambda do not match. Dimensions are:"
594          PRINT *, "inifor: rotate_to_cosmo: phi: ", SIZE(phi, 1), SIZE(phi, 2)
595          PRINT *, "inifor: rotate_to_cosmo: lam: ", SIZE(lam, 1), SIZE(lam, 2)
596          STOP
597       END IF
598
599       IF ( SIZE(phir) .NE. SIZE(phi, 2) .OR. &
600            SIZE(lamr) .NE. SIZE(phi, 1) )  THEN
601          PRINT *, "inifor: rotate_to_cosmo: Dimensions of phir and lamr do not match. Dimensions are:"
602          PRINT *, "inifor: rotate_to_cosmo: phir: ", SIZE(phir), SIZE(phi, 2)
603          PRINT *, "inifor: rotate_to_cosmo: lamr: ", SIZE(lamr), SIZE(phi, 1)
604          STOP
605       END IF
606       
607       DO j = 0, UBOUND(phir, 1)
608          DO i = 0, UBOUND(lamr, 1)
609
610             phi(i,j) = phirot2phi(phir(j) * TO_DEGREES,                       &
611                                   lamr(i) * TO_DEGREES,                       &
612                                   phip * TO_DEGREES,                          &
613                                   lamp * TO_DEGREES,                          &
614                                   gam  * TO_DEGREES) * TO_RADIANS
615
616             lam(i,j) = rlarot2rla(phir(j) * TO_DEGREES,                       &
617                                   lamr(i) * TO_DEGREES,                       &
618                                   phip * TO_DEGREES,                          &
619                                   lamp * TO_DEGREES,                          &
620                                   gam  * TO_DEGREES) * TO_RADIANS
621
622          END DO
623       END DO
624
625    END SUBROUTINE rotate_to_cosmo
626       
627
628!------------------------------------------------------------------------------!
629! Description:
630! ------------
631!> Rotate the given vector field (x(:), y(:)) by the given 'angle'.
632!------------------------------------------------------------------------------!
633    SUBROUTINE rotate_vector_field(x, y, angle)
634       REAL(dp), DIMENSION(:), INTENT(INOUT) :: x, y  !< x and y coodrinate in arbitrary units
635       REAL(dp), INTENT(IN)                  :: angle !< rotation angle [deg]
636
637       INTEGER  :: i
638       REAL(dp) :: sine, cosine, v_rot(2), rotation(2,2)
639
640       sine = SIN(angle * TO_RADIANS)
641       cosine = COS(angle * TO_RADIANS)
642!
643!--    RESAHPE() fills columns first, so the rotation matrix becomes
644!--   
645!--    rotation = [ cosine   -sine  ]
646!--               [  sine    cosine ]
647       rotation = RESHAPE( (/cosine, sine, -sine, cosine/), (/2, 2/) )
648
649       DO i = LBOUND(x, 1), UBOUND(x, 1)
650
651          v_rot(:) = MATMUL(rotation, (/x(i), y(i)/))
652
653          x(i) = v_rot(1)
654          y(i) = v_rot(2)
655
656       END DO
657
658    END SUBROUTINE rotate_vector_field
659
660
661
662!------------------------------------------------------------------------------!
663! Description:
664! ------------
665!> This routine computes the local meridian convergence between a rotated-pole
666!> and a geographical system using the Eq. (6) given in the DWD manual
667!>
668!>    Baldauf et al. (2018), Beschreibung des operationelle Kürzestfrist-
669!>    vorhersagemodells COSMO-D2 und COSMO-D2-EPS und seiner Ausgabe in die
670!>    Datenbanken des DWD.
671!>    https://www.dwd.de/SharedDocs/downloads/DE/modelldokumentationen/nwv/cosmo_d2/cosmo_d2_dbbeschr_aktuell.pdf?__blob=publicationFile&v=2
672!------------------------------------------------------------------------------!
673    FUNCTION meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g)          &
674       RESULT(delta)
675
676       REAL(dp), INTENT(IN) ::  phi_n, lam_n, phi_g, lam_g
677       REAL(dp)             ::  delta
678
679       delta = atan2( COS(phi_n) * SIN(lam_n - lam_g),                         &
680                      COS(phi_g) * SIN(phi_n) -                                &
681                      SIN(phi_g) * COS(phi_n) * COS(lam_n - lam_g) )
682
683    END FUNCTION meridian_convergence_rotated
684
685!------------------------------------------------------------------------------!
686! Description:
687! ------------
688!> Compute indices of PALM-4U grid point neighbours in the target
689!> system (COSMO-DE) by rounding up and down. (i,j) are the indices of
690!> the PALM-4U grid and (ii(i,j,1-4), jj(i,j,1-4)) contain the indices
691!> of the its four neigbouring points in the COSMO-DE grid.
692!>
693!>
694!>                     COSMO-DE grid
695!>                     -------------
696!>           jj, lat
697!>              ^        j
698!>              |         \          i
699!>  jj(i,j,2/3) + ... 2 ---\--------/------ 3
700!>              |     | ^   \      /        |
701!>              |     | |wp  \    /         |
702!>              |     | v     \  /          |
703!>       latpos + ............ o/ (i,j)     |
704!>              |     |        :            |
705!>              |     |        :<----wl---->|
706!>  jj(i,j,1/4) + ... 1 -------:----------- 4
707!>              |     :        :            :
708!>              |     :        :            :
709!>              |     :      lonpos         :
710!>              L-----+--------+------------+------> ii, lon
711!>               ii(i,j,1/2)        ii(i,j,3/4)
712!>
713!>
714!> Input parameters:
715!> -----------------
716!> source_lat, source_lon : (rotated-pole) coordinates of the source grid (e.g.
717!>    COSMO-DE)
718!>
719!> source_dxi, source_dyi : inverse grid spacings of the source grid.
720!>
721!> target_lat, target_lon : (rotated-pole) coordinates of the target grid (e.g.
722!>    COSMO-DE)
723!>
724!> Output parameters:
725!> ------------------
726!> palm_ii, palm_jj : x and y index arrays of horizontal neighbour columns
727!>
728!------------------------------------------------------------------------------!
729    SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon,                &
730                                          palm_clat, palm_clon,                &
731                                          palm_ii, palm_jj)
732
733       REAL(dp), DIMENSION(0:), INTENT(IN)        ::  cosmo_lat, cosmo_lon
734       REAL(dp), DIMENSION(0:,0:), INTENT(IN)     ::  palm_clat, palm_clon
735       REAL(dp)                                   ::  cosmo_dxi, cosmo_dyi
736       INTEGER, DIMENSION(0:,0:,1:), INTENT(OUT)  ::  palm_ii, palm_jj
737
738       REAL(dp) ::  lonpos, latpos, lon0, lat0
739       INTEGER  ::  i, j
740
741       lon0 = cosmo_lon(0)
742       lat0 = cosmo_lat(0)
743       cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))
744       cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))
745
746       DO j = 0, UBOUND(palm_clon, 2)!palm_grid % ny
747       DO i = 0, UBOUND(palm_clon, 1)!palm_grid % nx
748!
749!--       Compute the floating point index corrseponding to PALM-4U grid point
750!--       location along target grid (COSMO-DE) axes.
751          lonpos = (palm_clon(i,j) - lon0) * cosmo_dxi
752          latpos = (palm_clat(i,j) - lat0) * cosmo_dyi
753
754          IF (lonpos < 0.0 .OR. latpos < 0.0)  THEN
755             PRINT *, " Error while finding neighbours: lonpos or latpos out of bounds!"
756             PRINT *, "     (i,j) = (", i, ",",j,")"
757             PRINT *, "      lonpos ", lonpos*TO_DEGREES, ", latpos ", latpos*TO_DEGREES
758             PRINT *, "        lon0 ", lon0  *TO_DEGREES,   ", lat0   ", lat0*TO_DEGREES
759             PRINT *, "    PALM lon ", palm_clon(i,j)*TO_DEGREES,   ", PALM lat ",palm_clat(i,j)*TO_DEGREES
760             STOP
761          END IF
762
763          palm_ii(i,j,1) = FLOOR(lonpos)
764          palm_ii(i,j,2) = FLOOR(lonpos)
765          palm_ii(i,j,3) = CEILING(lonpos)
766          palm_ii(i,j,4) = CEILING(lonpos)
767
768          palm_jj(i,j,1) = FLOOR(latpos)
769          palm_jj(i,j,2) = CEILING(latpos)
770          palm_jj(i,j,3) = CEILING(latpos)
771          palm_jj(i,j,4) = FLOOR(latpos)
772       END DO
773       END DO
774
775    END SUBROUTINE find_horizontal_neighbours
776
777   
778!------------------------------------------------------------------------------!
779! Description:
780! ------------
781!> Computes linear vertical interpolation neighbour indices and weights for each
782!> column of the given palm grid.
783!------------------------------------------------------------------------------!
784    SUBROUTINE find_vertical_neighbours_and_weights_interp( palm_grid,         &
785                                                            palm_intermediate )
786       TYPE(grid_definition), INTENT(INOUT) ::  palm_grid
787       TYPE(grid_definition), INTENT(IN)    ::  palm_intermediate
788
789       INTEGER  ::  i, j, k, nx, ny, nz, nlev, k_intermediate
790       LOGICAL  ::  point_is_below_grid, point_is_above_grid,                  &
791                    point_is_in_current_cell
792       REAL(dp) ::  current_height, column_base, column_top, h_top, h_bottom,  &
793                    weight
794
795       nx   = palm_grid % nx
796       ny   = palm_grid % ny
797       nz   = palm_grid % nz
798       nlev = palm_intermediate % nz
799
800!
801!--    in each column of the fine grid, find vertical neighbours of every cell
802       DO j = 0, ny
803       DO i = 0, nx
804
805          k_intermediate = 0
806
807          column_base = palm_intermediate % h(i,j,0)
808          column_top  = palm_intermediate % h(i,j,nlev)
809
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.
815          DO k = 1, nz
816
817!
818!--          Memorize the top and bottom boundaries of the coarse cell and the
819!--          current height within it
820             current_height = palm_grid % z(k) + palm_grid % z0
821             h_top    = palm_intermediate % h(i,j,k_intermediate+1)
822             h_bottom = palm_intermediate % h(i,j,k_intermediate)
823
824             point_is_above_grid = (current_height > column_top) !22000m, very unlikely
825             point_is_below_grid = (current_height < column_base)
826
827             point_is_in_current_cell = (                                      &
828                current_height >= h_bottom .AND.                               &
829                current_height <  h_top                                        &
830             )
831
832!
833!--          set default weights
834             palm_grid % w_verti(i,j,k,1:2) = 0.0_dp
835
836             IF (point_is_above_grid)  THEN
837
838                palm_grid % kk(i,j,k,1:2) = nlev
839                palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp
840
841                message = "PALM-4U grid extends above COSMO-DE model top."
842                CALL abort('find_vertical_neighbours_and_weights', message)
843
844             ELSE IF (point_is_below_grid)  THEN
845
846                palm_grid % kk(i,j,k,1:2) = 0
847                palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp
848
849             ELSE
850!
851!--             cycle through intermediate levels until current
852!--             intermediate-grid cell overlaps with current_height
853                DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
854                   k_intermediate = k_intermediate + 1
855
856                   h_top    = palm_intermediate % h(i,j,k_intermediate+1)
857                   h_bottom = palm_intermediate % h(i,j,k_intermediate)
858                   point_is_in_current_cell = (                                &
859                      current_height >= h_bottom .AND.                         &
860                      current_height <  h_top                                  &
861                   )
862                END DO
863
864                IF (k_intermediate > nlev-1)  THEN
865                   message = "Index " // TRIM(str(k_intermediate)) //          &
866                             " is above intermediate grid range."
867                   CALL abort('find_vertical_neighbours', message)
868                END IF
869   
870                palm_grid % kk(i,j,k,1) = k_intermediate
871                palm_grid % kk(i,j,k,2) = k_intermediate + 1
872
873!
874!--             compute vertical weights
875                weight = (h_top - current_height) / (h_top - h_bottom)
876                palm_grid % w_verti(i,j,k,1) = weight
877                palm_grid % w_verti(i,j,k,2) = 1.0_dp - weight
878             END IF
879
880          END DO
881
882       END DO
883       END DO
884
885    END SUBROUTINE find_vertical_neighbours_and_weights_interp
886
887
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!------------------------------------------------------------------------------!
900    SUBROUTINE find_vertical_neighbours_and_weights_average( avg_grid )
901       TYPE(grid_definition), INTENT(INOUT) ::  avg_grid
902
903       INTEGER  ::  i, j, k_palm, k_intermediate, l, nlev
904       LOGICAL  ::  point_is_below_grid, point_is_above_grid,                  &
905                    point_is_in_current_cell
906       REAL(dp) ::  current_height, column_base, column_top, h_top, h_bottom,  &
907                    weight
908
909
910       avg_grid % k_min = LBOUND(avg_grid % z, 1)
911
912       nlev = SIZE(avg_grid % cosmo_h, 3)
913
914!
915!--    in each column of the fine grid, find vertical neighbours of every cell
916       DO l = 1, avg_grid % n_columns
917
918          i = avg_grid % iii(l)
919          j = avg_grid % jjj(l)
920
921          column_base = avg_grid % cosmo_h(i,j,1)
922          column_top  = avg_grid % cosmo_h(i,j,nlev)
923
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.
929          k_intermediate = 1 !avg_grid % cosmo_h is indezed 1-based.
930          DO k_palm = 1, avg_grid % nz
931
932!
933!--          Memorize the top and bottom boundaries of the coarse cell and the
934!--          current height within it
935             current_height = avg_grid % z(k_palm) + avg_grid % z0
936             h_top    = avg_grid % cosmo_h(i,j,k_intermediate+1)
937             h_bottom = avg_grid % cosmo_h(i,j,k_intermediate)
938
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)
943             point_is_below_grid = (current_height < column_base)
944
945             point_is_in_current_cell = (                                      &
946                current_height >= h_bottom .AND.                               &
947                current_height <  h_top                                        &
948             )
949
950!
951!--          set default weights
952             avg_grid % w(l,k_palm,1:2) = 0.0_dp
953
954             IF (point_is_above_grid)  THEN
955
956                avg_grid % kkk(l,k_palm,1:2) = nlev
957                avg_grid % w(l,k_palm,1:2) = - 2.0_dp
958
959                message = "PALM-4U grid extends above COSMO-DE model top."
960                CALL abort('find_vertical_neighbours_and_weights_average', message)
961
962             ELSE IF (point_is_below_grid)  THEN
963
964                avg_grid % kkk(l,k_palm,1:2) = 0
965                avg_grid % w(l,k_palm,1:2) = - 2.0_dp
966                avg_grid % k_min = MAX(k_palm + 1, avg_grid % k_min)
967             ELSE
968!
969!--             cycle through intermediate levels until current
970!--             intermediate-grid cell overlaps with current_height
971                DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
972                   k_intermediate = k_intermediate + 1
973
974                   h_top    = avg_grid % cosmo_h(i,j,k_intermediate+1)
975                   h_bottom = avg_grid % cosmo_h(i,j,k_intermediate)
976                   point_is_in_current_cell = (                                &
977                      current_height >= h_bottom .AND.                         &
978                      current_height <  h_top                                  &
979                   )
980                END DO
981
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.
985                IF (k_intermediate > nlev-1)  THEN
986                   message = "Index " // TRIM(str(k_intermediate)) //          &
987                             " is above intermediate grid range."
988                   CALL abort('find_vertical_neighbours', message)
989                END IF
990   
991                avg_grid % kkk(l,k_palm,1) = k_intermediate
992                avg_grid % kkk(l,k_palm,2) = k_intermediate + 1
993
994!
995!--             compute vertical weights
996                weight = (h_top - current_height) / (h_top - h_bottom)
997                avg_grid % w(l,k_palm,1) = weight
998                avg_grid % w(l,k_palm,2) = 1.0_dp - weight
999             END IF
1000
1001!
1002!--       Loop over PALM levels k
1003          END DO
1004
1005!
1006!--       Loop over averaging columns l
1007       END DO
1008 
1009    END SUBROUTINE find_vertical_neighbours_and_weights_average
1010
1011!------------------------------------------------------------------------------!
1012! Description:
1013! ------------
1014!> Compute the four weights for horizontal bilinear interpolation given the
1015!> coordinates clon(i,j) clat(i,j) of the PALM-4U grid in the COSMO-DE
1016!> rotated-pole grid and the neightbour indices ii(i,j,1-4) and jj(i,j,1-4).
1017!>
1018!> Input parameters:
1019!> -----------------
1020!> palm_grid % clon : longitudes of PALM-4U scalars (cell centres) in COSMO-DE's rotated-pole grid [rad]
1021!>
1022!> palm_grid % clat : latitudes of PALM-4U cell centres in COSMO-DE's rotated-pole grid [rad]
1023!>
1024!> cosmo_grid % lon : rotated-pole longitudes of scalars (cell centres) of the COSMO-DE grid [rad]
1025!>
1026!> cosmo_grid % lat : rotated-pole latitudes of scalars (cell centers) of the COSMO-DE grid [rad]
1027!>
1028!> cosmo_grid % dxi : inverse grid spacing in the first dimension [m^-1]
1029!>
1030!> cosmo_grid % dyi : inverse grid spacing in the second dimension [m^-1]
1031!>
1032!> Output parameters:
1033!> ------------------
1034!> palm_grid % w_horiz(:,:,1-4) : weights for bilinear horizontal interpolation
1035!
1036!                               COSMO-DE grid
1037!                               -------------
1038!                     jj, lat
1039!                        ^        j
1040!                        |         \          i
1041!            jj(i,j,2/3) + ... 2 ---\--------/------ 3
1042!                        |     | ^   \      /        |
1043!                        |     | |wp  \    /         |
1044!                        |     | v     \  /          |
1045!                 latpos + ............ o/ (i,j)     |
1046!                        |     |        :            |
1047!                        |     |        :<----wl---->|
1048!            jj(i,j,1/4) + ... 1 -------:----------- 4
1049!                        |     :        :            :
1050!                        |     :        :            :
1051!                        |     :      lonpos         :
1052!                        L-----+--------+------------+------> ii, lon
1053!                         ii(i,j,1/2)        ii(i,j,3/4)
1054!         
1055!------------------------------------------------------------------------------!
1056    SUBROUTINE compute_horizontal_interp_weights(cosmo_lat, cosmo_lon,         &
1057       palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)
1058       
1059       REAL(dp), DIMENSION(0:), INTENT(IN)        ::  cosmo_lat, cosmo_lon
1060       REAL(dp)                                   ::  cosmo_dxi, cosmo_dyi
1061       REAL(dp), DIMENSION(0:,0:), INTENT(IN)     ::  palm_clat, palm_clon
1062       INTEGER, DIMENSION(0:,0:,1:), INTENT(IN)   ::  palm_ii, palm_jj
1063
1064       REAL(dp), DIMENSION(0:,0:,1:), INTENT(OUT) ::  palm_w_horiz
1065
1066       REAL(dp) ::  wl, wp
1067       INTEGER  ::  i, j
1068
1069       cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))
1070       cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))
1071
1072       DO j = 0, UBOUND(palm_clon, 2)
1073       DO i = 0, UBOUND(palm_clon, 1)
1074     
1075!
1076!--       weight in lambda direction
1077          wl = ( cosmo_lon(palm_ii(i,j,4)) - palm_clon(i,j) ) * cosmo_dxi
1078
1079!
1080!--       weight in phi direction
1081          wp = ( cosmo_lat(palm_jj(i,j,2)) - palm_clat(i,j) ) * cosmo_dyi
1082
1083          IF (wl > 1.0_dp .OR. wl < 0.0_dp)  THEN
1084              message = "Horizontal weight wl = " // TRIM(real_to_str(wl)) //   &
1085                        " is out bounds."
1086              CALL abort('compute_horizontal_interp_weights', message)
1087          END IF
1088          IF (wp > 1.0_dp .OR. wp < 0.0_dp)  THEN
1089              message = "Horizontal weight wp = " // TRIM(real_to_str(wp)) //   &
1090                        " is out bounds."
1091              CALL abort('compute_horizontal_interp_weights', message)
1092          END IF
1093
1094          palm_w_horiz(i,j,1) = wl * wp
1095          palm_w_horiz(i,j,2) = wl * (1.0_dp - wp)
1096          palm_w_horiz(i,j,3) = (1.0_dp - wl) * (1.0_dp - wp)
1097          palm_w_horiz(i,j,4) = 1.0_dp - SUM( palm_w_horiz(i,j,1:3) )
1098
1099       END DO
1100       END DO
1101       
1102    END SUBROUTINE compute_horizontal_interp_weights
1103
1104
1105!------------------------------------------------------------------------------!
1106! Description:
1107! ------------
1108!> Interpolates u and v components of velocities located at cell faces to the
1109!> cell centres by averaging neighbouring values.
1110!>
1111!> This routine is designed to be used with COSMO-DE arrays where there are the
1112!> same number of grid points for scalars (centres) and velocities (faces). In
1113!> COSMO-DE the velocity points are staggared one half grid spaceing up-grid
1114!> which means the first centre point has to be omitted and is set to zero.
1115!------------------------------------------------------------------------------!
1116    SUBROUTINE centre_velocities(u_face, v_face, u_centre, v_centre)
1117       REAL(dp), DIMENSION(0:,0:,0:), INTENT(IN)  ::  u_face, v_face
1118       REAL(dp), DIMENSION(0:,0:,0:), INTENT(OUT) ::  u_centre, v_centre
1119       INTEGER ::  nx, ny
1120
1121       nx = UBOUND(u_face, 1)
1122       ny = UBOUND(u_face, 2)
1123
1124       u_centre(0,:,:)  = 0.0_dp
1125       u_centre(1:,:,:) = 0.5_dp * ( u_face(0:nx-1,:,:) + u_face(1:,:,:) )
1126
1127       v_centre(:,0,:)  = 0.0_dp
1128       v_centre(:,1:,:) = 0.5_dp * ( v_face(:,0:ny-1,:) + v_face(:,1:,:) )
1129    END SUBROUTINE centre_velocities
1130
1131
1132!------------------------------------------------------------------------------!
1133! Description:
1134! ------------
1135!> Compute the geographical latitude of a point given in rotated-pole cordinates
1136!------------------------------------------------------------------------------!
1137    FUNCTION phirot2phi (phirot, rlarot, polphi, pollam, polgam)
1138   
1139       REAL(dp), INTENT (IN) ::  polphi      !< latitude of the rotated north pole
1140       REAL(dp), INTENT (IN) ::  pollam      !< longitude of the rotated north pole
1141       REAL(dp), INTENT (IN) ::  phirot      !< latitude in the rotated system
1142       REAL(dp), INTENT (IN) ::  rlarot      !< longitude in the rotated system
1143       REAL(dp), INTENT (IN) ::  polgam      !< angle between the north poles of the systems
1144
1145       REAL(dp)              ::  phirot2phi  !< latitude in the geographical system
1146       
1147       REAL(dp)              ::  zsinpol, zcospol, zphis, zrlas, zarg, zgam
1148   
1149       zsinpol = SIN(polphi * TO_RADIANS)
1150       zcospol = COS(polphi * TO_RADIANS)
1151       zphis   = phirot * TO_RADIANS
1152
1153       IF (rlarot > 180.0_dp)  THEN
1154          zrlas = rlarot - 360.0_dp
1155       ELSE
1156          zrlas = rlarot
1157       END IF
1158       zrlas = zrlas * TO_RADIANS
1159     
1160       IF (polgam /= 0.0_dp)  THEN
1161          zgam = polgam * TO_RADIANS
1162          zarg = zsinpol * SIN (zphis) +                                       &
1163                 zcospol * COS(zphis) * ( COS(zrlas) * COS(zgam) -             &
1164                                          SIN(zgam)  * SIN(zrlas) )
1165       ELSE
1166          zarg = zcospol * COS (zphis) * COS (zrlas) + zsinpol * SIN (zphis)
1167       END IF
1168     
1169       phirot2phi = ASIN (zarg) * TO_DEGREES
1170   
1171    END FUNCTION phirot2phi
1172
1173
1174!------------------------------------------------------------------------------!
1175! Description:
1176! ------------
1177!> Compute the geographical latitude of a point given in rotated-pole cordinates
1178!------------------------------------------------------------------------------!
1179    FUNCTION phi2phirot (phi, rla, polphi, pollam)
1180   
1181       REAL(dp), INTENT (IN) ::  polphi !< latitude of the rotated north pole
1182       REAL(dp), INTENT (IN) ::  pollam !< longitude of the rotated north pole
1183       REAL(dp), INTENT (IN) ::  phi    !< latitude in the geographical system
1184       REAL(dp), INTENT (IN) ::  rla    !< longitude in the geographical system
1185       
1186       REAL(dp) ::  phi2phirot          !< longitude in the rotated system
1187       
1188       REAL(dp) ::  zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1
1189       
1190       zsinpol = SIN(polphi * TO_RADIANS)
1191       zcospol = COS(polphi * TO_RADIANS)
1192       zlampol = pollam * TO_RADIANS
1193       zphi    = phi * TO_RADIANS
1194
1195       IF (rla > 180.0_dp)  THEN
1196          zrla1 = rla - 360.0_dp
1197       ELSE
1198          zrla1 = rla
1199       END IF
1200       zrla = zrla1 * TO_RADIANS
1201       
1202       zarg1 = SIN(zphi) * zsinpol
1203       zarg2 = COS(zphi) * zcospol * COS(zrla - zlampol)
1204       
1205       phi2phirot = ASIN(zarg1 + zarg2) * TO_DEGREES
1206   
1207    END FUNCTION phi2phirot
1208
1209
1210!------------------------------------------------------------------------------!
1211! Description:
1212! ------------
1213!> Compute the geographical longitude of a point given in rotated-pole cordinates
1214!------------------------------------------------------------------------------!
1215    FUNCTION rlarot2rla(phirot, rlarot, polphi, pollam, polgam)
1216   
1217       REAL(dp), INTENT (IN) ::  polphi !< latitude of the rotated north pole
1218       REAL(dp), INTENT (IN) ::  pollam !< longitude of the rotated north pole
1219       REAL(dp), INTENT (IN) ::  phirot !< latitude in the rotated system
1220       REAL(dp), INTENT (IN) ::  rlarot !< longitude in the rotated system
1221       REAL(dp), INTENT (IN) ::  polgam !< angle between the north poles of the systems
1222       
1223       REAL(dp) ::  rlarot2rla          !< latitude in the geographical system
1224       
1225       REAL(dp) ::  zsinpol, zcospol, zlampol, zphis, zrlas, zarg1, zarg2, zgam
1226       
1227       zsinpol = SIN(TO_RADIANS * polphi)
1228       zcospol = COS(TO_RADIANS * polphi)
1229       zlampol = TO_RADIANS * pollam
1230       zphis   = TO_RADIANS * phirot
1231
1232       IF (rlarot > 180.0_dp)  THEN
1233          zrlas = rlarot - 360.0_dp
1234       ELSE
1235          zrlas = rlarot
1236       END IF
1237       zrlas   = TO_RADIANS * zrlas
1238     
1239       IF (polgam /= 0.0_dp)  THEN
1240          zgam  = TO_RADIANS * polgam
1241          zarg1 = SIN(zlampol) * (zcospol * SIN(zphis) - zsinpol*COS(zphis) *  &
1242                  (COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) -        &
1243                  COS(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) +       &
1244                                                COS(zrlas) * SIN(zgam) )
1245       
1246          zarg2 = COS (zlampol) * (zcospol * SIN(zphis) - zsinpol*COS(zphis) * &
1247                  (COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) +        &
1248                  SIN(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) +       &
1249                                                COS(zrlas) * SIN(zgam) )
1250       ELSE
1251          zarg1   = SIN (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis)  +     &
1252                                      zcospol *              SIN(zphis)) -     &
1253                    COS (zlampol) *             SIN(zrlas) * COS(zphis)
1254          zarg2   = COS (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis)  +     &
1255                                      zcospol *              SIN(zphis)) +     &
1256                    SIN (zlampol) *             SIN(zrlas) * COS(zphis)
1257       END IF
1258     
1259       IF (zarg2 == 0.0_dp)  zarg2 = 1.0E-20_dp
1260     
1261       rlarot2rla = ATAN2(zarg1,zarg2) * TO_DEGREES
1262       
1263    END FUNCTION rlarot2rla
1264
1265
1266!------------------------------------------------------------------------------!
1267! Description:
1268! ------------
1269!> Compute the rotated-pole longitude of a point given in geographical cordinates
1270!------------------------------------------------------------------------------!
1271    FUNCTION rla2rlarot ( phi, rla, polphi, pollam, polgam )
1272
1273       REAL(dp), INTENT (IN) ::  polphi !< latitude of the rotated north pole
1274       REAL(dp), INTENT (IN) ::  pollam !< longitude of the rotated north pole
1275       REAL(dp), INTENT (IN) ::  phi    !< latitude in geographical system
1276       REAL(dp), INTENT (IN) ::  rla    !< longitude in geographical system
1277       REAL(dp), INTENT (IN) ::  polgam !< angle between the north poles of the systems
1278       
1279       REAL (KIND=dp) ::  rla2rlarot    !< latitude in the the rotated system
1280       
1281       REAL (KIND=dp) ::  zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1
1282       
1283       zsinpol = SIN(polphi * TO_RADIANS)
1284       zcospol = COS(polphi * TO_RADIANS)
1285       zlampol = pollam * TO_RADIANS
1286       zphi    = phi * TO_RADIANS
1287
1288       IF (rla > 180.0_dp)  THEN
1289          zrla1 = rla - 360.0_dp
1290       ELSE
1291          zrla1 = rla
1292       END IF
1293       zrla = zrla1 * TO_RADIANS
1294       
1295       zarg1 = - SIN (zrla-zlampol) * COS(zphi)
1296       zarg2 = - zsinpol * COS(zphi) * COS(zrla-zlampol) + zcospol * SIN(zphi)
1297       
1298       IF (zarg2 == 0.0_dp)  zarg2 = 1.0E-20_dp
1299       
1300       rla2rlarot = ATAN2 (zarg1,zarg2) * TO_DEGREES
1301       
1302       IF (polgam /= 0.0_dp )  THEN
1303          rla2rlarot = polgam + rla2rlarot
1304          IF (rla2rlarot > 180._dp)  rla2rlarot = rla2rlarot - 360.0_dp
1305       END IF
1306       
1307    END FUNCTION rla2rlarot
1308
1309
1310!------------------------------------------------------------------------------!
1311! Description:
1312! ------------
1313!> Rotate the given velocity vector (u,v) from the geographical to the
1314!> rotated-pole system
1315!------------------------------------------------------------------------------!
1316    SUBROUTINE uv2uvrot(u, v, rlat, rlon, pollat, pollon, urot, vrot)
1317   
1318       REAL(dp), INTENT (IN)  ::  u, v           !< wind components in the true geographical system
1319       REAL(dp), INTENT (IN)  ::  rlat, rlon     !< coordinates in the true geographical system
1320       REAL(dp), INTENT (IN)  ::  pollat, pollon !< latitude and longitude of the north pole of the rotated grid
1321       
1322       REAL(dp), INTENT (OUT) ::  urot, vrot     !< wind components in the rotated grid             
1323       
1324       REAL (dp) ::  zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm
1325       
1326       zsinpol = SIN(pollat * TO_RADIANS)
1327       zcospol = COS(pollat * TO_RADIANS)
1328       zlonp   = (pollon-rlon) * TO_RADIANS
1329       zlat    = rlat * TO_RADIANS
1330       
1331       zarg1 = zcospol * SIN(zlonp)
1332       zarg2 = zsinpol * COS(zlat) - zcospol * SIN(zlat) * COS(zlonp)
1333       znorm = 1.0_dp / SQRT(zarg1*zarg1 + zarg2*zarg2)
1334       
1335       urot = u * zarg2 * znorm - v * zarg1 * znorm
1336       vrot = u * zarg1 * znorm + v * zarg2 * znorm
1337   
1338    END SUBROUTINE uv2uvrot
1339
1340
1341!------------------------------------------------------------------------------!
1342! Description:
1343! ------------
1344!> Rotate the given velocity vector (urot, vrot) from the rotated-pole to the
1345!> geographical system
1346!------------------------------------------------------------------------------!
1347    SUBROUTINE uvrot2uv (urot, vrot, rlat, rlon, pollat, pollon, u, v)
1348   
1349       REAL(dp), INTENT(IN) ::  urot, vrot     !< wind components in the rotated grid
1350       REAL(dp), INTENT(IN) ::  rlat, rlon     !< latitude and longitude in the true geographical system
1351       REAL(dp), INTENT(IN) ::  pollat, pollon !< latitude and longitude of the north pole of the rotated grid
1352       
1353       REAL(dp), INTENT(OUT) ::  u, v          !< wind components in the true geographical system
1354       
1355       REAL(dp) ::  zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm
1356     
1357       zsinpol = SIN(pollat * TO_RADIANS)
1358       zcospol = COS(pollat * TO_RADIANS)
1359       zlonp   = (pollon-rlon) * TO_RADIANS
1360       zlat    = rlat * TO_RADIANS
1361     
1362       zarg1 = zcospol * SIN(zlonp)
1363       zarg2 = zsinpol * COS(zlat) - zcospol * SIN(zlat) * COS(zlonp)
1364       znorm = 1.0_dp / SQRT(zarg1*zarg1 + zarg2*zarg2)
1365     
1366       u =   urot * zarg2 * znorm + vrot * zarg1 * znorm
1367       v = - urot * zarg1 * znorm + vrot * zarg2 * znorm
1368   
1369    END SUBROUTINE uvrot2uv
1370
1371 END MODULE
1372
Note: See TracBrowser for help on using the repository browser.