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

Last change on this file since 4004 was 3866, checked in by eckhard, 5 years ago

inifor: Use PALM's working precision; improved error handling, coding style, and comments

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