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

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

inifor: bugfix: Fixes issue #815 with geostrophic wind profiles

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