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

Last change on this file since 3613 was 3613, checked in by eckhard, 6 years ago

inifor: Average initial profiles over the PALM, not the geostrophic, region

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