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

Last change on this file since 3698 was 3680, checked in by knoop, 6 years ago

Added cpp-option netcdf to inifor

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