source: palm/trunk/UTIL/inifor/src/transform.f90 @ 3395

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

inifor: Added computation of geostrophic winds from COSMO input

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