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

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

inifor: COSMO-D2 support

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