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

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

inifor: Renamed source files for compatibilty with PALM build system

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