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

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

inifor: Removed unused variables, improved coding style

  • Property svn:keywords set to Id
File size: 59.6 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 3785 2019-03-06 10:41:14Z 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 :: 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             ENDDO
142          ENDIF
143       ENDDO
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             ENDDO
198          ENDIF
199       ENDDO
200       ENDDO
201       ENDDO
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       ENDIF
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          ENDDO
260       ENDDO
261       ENDDO
262       ENDDO
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       ENDIF
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       ENDIF
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          ENDDO
302
303       ENDDO
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             ENDDO
388
389!
390!--       Loop over PALM levels k_profile
391          ENDDO
392
393!
394!--    Loop over horizontal neighbours l
395       ENDDO
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, avg_grid )
414
415       TYPE(grid_definition), INTENT(IN)          ::  avg_grid
416       REAL(dp), DIMENSION(:,:,:), INTENT(IN)     ::  source_array
417       REAL(dp), DIMENSION(:), INTENT(OUT)        ::  profile_array
418
419       INTEGER ::  i_source, j_source, l, nz, nlev
420
421       REAL(dp) ::  ni_columns
422
423       nlev = SIZE( source_array, 3 )
424       nz   = SIZE( profile_array, 1 )
425
426       IF ( nlev /= nz )  THEN
427          message = "Lengths of input and output profiles do not match: " //   &
428                    "cosmo_pressure(" // TRIM( str( nlev ) ) //                &
429                    "), profile_array(" // TRIM( str( nz ) )  // ")."
430          CALL inifor_abort('average_pressure_perturbation', message)
431       ENDIF
432     
433       profile_array(:) = 0.0_dp
434
435       DO l = 1, avg_grid % n_columns
436
437          i_source = avg_grid % iii(l)
438          j_source = avg_grid % jjj(l)
439
440          profile_array(:) = profile_array(:)                                  &
441                           + source_array(i_source, j_source, :)
442
443       ENDDO
444
445       ni_columns = 1.0_dp / avg_grid % n_columns
446       profile_array(:) = profile_array(:) * ni_columns
447
448    END SUBROUTINE average_profile
449
450
451!------------------------------------------------------------------------------!
452! Description:
453! ------------
454!> This is a sister routine to average_profile() and differes from it in that
455!> it removes the COSMO basic state pressure from the input array before
456!> averaging.
457!------------------------------------------------------------------------------!
458    SUBROUTINE average_pressure_perturbation( cosmo_pressure, profile_array,   &
459                                              cosmo_grid, avg_grid )
460
461       TYPE(grid_definition), INTENT(IN)          ::  cosmo_grid, avg_grid
462       REAL(dp), DIMENSION(:,:,:), INTENT(IN)     ::  cosmo_pressure
463       REAL(dp), DIMENSION(:), INTENT(OUT)        ::  profile_array
464
465       INTEGER ::  i_source, j_source, l, nz, nlev
466
467       REAL(dp)                            ::  ni_columns
468       REAL(dp), DIMENSION(:), ALLOCATABLE ::  basic_state_pressure
469
470       nlev = SIZE( cosmo_pressure, 3 )
471       nz   = SIZE( profile_array, 1 )
472
473       IF ( nlev /= nz )  THEN
474          message = "Lengths of input and output profiles do not match: " //   &
475                    "cosmo_pressure(" // TRIM( str( nlev ) ) //                &
476                    "), profile_array(" // TRIM( str( nz ) )  // ")."
477          CALL inifor_abort('average_pressure_perturbation', message)
478       ENDIF
479
480       ALLOCATE( basic_state_pressure(nz) )
481       profile_array(:) = 0.0_dp
482
483       DO l = 1, avg_grid % n_columns
484          i_source = avg_grid % iii(l)
485          j_source = avg_grid % jjj(l)
486
487!
488!--       Compute pressure perturbation by removing COSMO basic state pressure
489          CALL get_basic_state( cosmo_grid % hfl(i_source,j_source,:), BETA,   &
490                                P_SL, T_SL, RD, G, basic_state_pressure )
491
492          profile_array(:) = profile_array(:)                                  &
493                           + cosmo_pressure(i_source, j_source, :)             &
494                           - basic_state_pressure(:)
495
496!
497!--    Loop over horizontal neighbours l
498       ENDDO
499
500       DEALLOCATE( basic_state_pressure )
501
502       ni_columns = 1.0_dp / avg_grid % n_columns
503       profile_array(:) = profile_array(:) * ni_columns
504
505    END SUBROUTINE average_pressure_perturbation
506
507
508
509
510!------------------------------------------------------------------------------!
511! Description:
512! ------------
513!> Extrapolates density linearly from the level 'k_min' downwards.
514!------------------------------------------------------------------------------!
515    SUBROUTINE extrapolate_density(rho, avg_grid)
516       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  rho
517       TYPE(grid_definition), INTENT(IN)     ::  avg_grid
518
519       REAL(dp) ::  drhodz, dz, zk, rhok
520       INTEGER  ::  k_min
521
522       k_min  = avg_grid % k_min
523       zk     = avg_grid % z(k_min)
524       rhok   = rho(k_min)
525       dz     = avg_grid % z(k_min + 1) - avg_grid % z(k_min)
526       drhodz = (rho(k_min + 1) - rho(k_min)) / dz
527
528       rho(1:k_min-1) = rhok + drhodz * (avg_grid % z(1:k_min-1) - zk)
529
530    END SUBROUTINE extrapolate_density
531
532
533!------------------------------------------------------------------------------!
534! Description:
535! ------------
536!> Driver for extrapolating pressure from PALM level k_min downwards
537!------------------------------------------------------------------------------!
538    SUBROUTINE extrapolate_pressure(p, rho, avg_grid)
539       REAL(dp), DIMENSION(:), INTENT(IN)    ::  rho
540       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  p
541       TYPE(grid_definition), INTENT(IN)     ::  avg_grid
542
543       REAL(dp) ::  drhodz, dz, zk, rhok
544       INTEGER  ::  k, k_min
545
546       k_min = avg_grid % k_min
547       zk    = avg_grid % z(k_min)
548       rhok  = rho(k_min)
549       dz    = avg_grid % z(k_min + 1) - avg_grid % z(k_min)
550       drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz
551
552       DO k = 1, k_min-1
553          p(k) = constant_density_pressure(p(k_min), zk, rhok, drhodz,         &
554                                           avg_grid % z(k), G)
555       ENDDO
556
557    END SUBROUTINE extrapolate_pressure
558
559
560!------------------------------------------------------------------------------!
561! Description:
562! ------------
563!> Takes the averaged pressure profile <p> and sets the lowest entry to the
564!> extrapolated pressure at the surface.
565!------------------------------------------------------------------------------!
566    SUBROUTINE get_surface_pressure(p, rho, avg_grid)
567       REAL(dp), DIMENSION(:), INTENT(IN)    ::  rho
568       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  p
569       TYPE(grid_definition), INTENT(IN)     ::  avg_grid
570
571       REAL(dp) ::  drhodz, dz, zk, rhok
572       INTEGER  ::  k_min
573
574       k_min = avg_grid % k_min
575       zk    = avg_grid % z(k_min)
576       rhok  = rho(k_min)
577       dz    = avg_grid % z(k_min + 1) - avg_grid % z(k_min)
578       drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz
579
580       p(1) = constant_density_pressure(p(k_min), zk, rhok, drhodz,            &
581                                        0.0_dp, G)
582
583    END SUBROUTINE get_surface_pressure
584
585
586    FUNCTION constant_density_pressure(pk, zk, rhok, drhodz, z, g)  RESULT(p)
587
588       REAL(dp), INTENT(IN)  ::  pk, zk, rhok, drhodz, g, z
589       REAL(dp) ::  p
590
591       p = pk + ( zk - z ) * g * ( rhok + 0.5*drhodz * (zk - z) )
592
593    END FUNCTION constant_density_pressure
594
595!-----------------------------------------------------------------------------!
596! Description:
597! -----------
598!> This routine computes profiles of the two geostrophic wind components ug and
599!> vg.
600!-----------------------------------------------------------------------------!
601    SUBROUTINE geostrophic_winds(p_north, p_south, p_east, p_west, rho, f3,    &
602                                 Lx, Ly, phi_n, lam_n, phi_g, lam_g, ug, vg)
603
604    REAL(dp), DIMENSION(:), INTENT(IN)  ::  p_north, p_south, p_east, p_west,  &
605                                            rho
606    REAL(dp), INTENT(IN)                ::  f3, Lx, Ly, phi_n, lam_n, phi_g, lam_g
607    REAL(dp), DIMENSION(:), INTENT(OUT) ::  ug, vg
608    REAL(dp)                            ::  facx, facy
609
610    facx = 1.0_dp / (Lx * f3)
611    facy = 1.0_dp / (Ly * f3)
612    ug(:) = - facy / rho(:) * (p_north(:) - p_south(:))
613    vg(:) =   facx / rho(:) * (p_east(:) - p_west(:))
614
615    CALL rotate_vector_field(                                                  &
616       ug, vg, angle = meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g)&
617    )
618
619    END SUBROUTINE geostrophic_winds
620
621
622!-----------------------------------------------------------------------------!
623! Description:
624! -----------
625!> This routine computes the inverse Plate Carree projection, i.e. in projects
626!> Cartesian coordinates (x,y) onto a sphere. It returns the latitude and
627!> lngitude of a geographical system centered at x0 and y0.
628!-----------------------------------------------------------------------------!
629    SUBROUTINE inv_plate_carree(x, y, x0, y0, r, lat, lon)
630       REAL(dp), INTENT(IN)  ::  x(:), y(:), x0, y0, r
631       REAL(dp), INTENT(OUT) ::  lat(:), lon(:)
632       
633       REAL(dp) :: ri
634
635!
636!--    TODO check dimensions of lat/lon and y/x match
637
638       ri = 1.0_dp / r
639       
640       lat(:) = (y(:) - y0) * ri
641       lon(:) = (x(:) - x0) * ri
642    END SUBROUTINE 
643
644
645!-----------------------------------------------------------------------------!
646! Description:
647! ------------
648!> Computes the reverse Plate-Carree projection of a x or y position on a
649!> Cartesian grid.
650!>
651!> Input parameters:
652!> -----------------
653!> xy : x or y coordinate of the Cartasian grid point [m].
654!>
655!> xy0 : x or y coordinate that coincides with the origin of the underlying
656!>     sperical system (crossing point of the equator and prime meridian) [m].
657!>
658!> r : Radius of the of the underlying sphere, e.g. EARTH_RADIUS [m].
659!>
660!> Returns:
661!> --------
662!> project : Longitude (in case xy = x) or latitude (xy = y) of the given input
663!>     coordinate xy.
664!------------------------------------------------------------------------------!
665    ELEMENTAL REAL(dp) FUNCTION project(xy, xy0, r)
666       REAL(dp), INTENT(IN)  ::  xy, xy0, r
667       REAL(dp) :: ri
668
669!
670!--    If this elemental function is called with a large array as xy, it is
671!--    computationally more efficient to precompute the inverse radius and
672!--    then muliply.
673       ri = 1.0_dp / r
674
675       project = (xy - xy0) * ri
676
677    END FUNCTION project
678
679
680!------------------------------------------------------------------------------!
681! Description:
682! ------------
683!> For a rotated-pole system with the origin at geographical latitude 'phi_c',
684!> compute the geographical latitude of its rotated north pole.
685!------------------------------------------------------------------------------!
686    REAL(dp) FUNCTION phic_to_phin(phi_c)
687        REAL(dp), INTENT(IN) ::  phi_c
688
689        phic_to_phin = 0.5_dp * PI - ABS(phi_c)
690
691    END FUNCTION phic_to_phin
692
693   
694!------------------------------------------------------------------------------!
695! Description:
696! ------------
697!> For a rotated-pole system with the origin at geographical latitude 'phi_c'
698!> and longitude 'lam_c', compute the geographical longitude of its rotated
699!> north pole.
700!------------------------------------------------------------------------------!
701    REAL(dp) FUNCTION lamc_to_lamn(phi_c, lam_c)
702       REAL(dp), INTENT(IN) ::  phi_c, lam_c
703       
704       lamc_to_lamn = lam_c
705       IF (phi_c > 0.0_dp)  THEN
706          lamc_to_lamn = lam_c - SIGN(PI, lam_c)
707       ENDIF
708
709    END FUNCTION lamc_to_lamn
710
711
712!------------------------------------------------------------------------------!
713! Description:
714! ------------
715!> Set gamma according to whether PALM domain is in the northern or southern
716!> hemisphere of the COSMO rotated-pole system. Gamma assumes either the
717!> value 0 or PI and is needed to work around around a bug in the
718!> rotated-pole coordinate transformations.
719!------------------------------------------------------------------------------!
720    REAL(dp) FUNCTION gamma_from_hemisphere(phi_cg, phi_ref)
721       REAL(dp), INTENT(IN) ::  phi_cg
722       REAL(dp), INTENT(IN) ::  phi_ref
723
724       LOGICAL ::  palm_origin_is_south_of_cosmo_origin
725       
726       palm_origin_is_south_of_cosmo_origin = (phi_cg < phi_ref)
727
728       IF (palm_origin_is_south_of_cosmo_origin)  THEN
729           gamma_from_hemisphere = PI
730       ELSE
731           gamma_from_hemisphere = 0.0_dp
732       ENDIF
733    END FUNCTION gamma_from_hemisphere
734
735
736!------------------------------------------------------------------------------!
737! Description:
738! ------------
739!> Computes the geographical coordinates corresponding to the given rotated-pole
740!> coordinates.
741!>
742!> In INIFOR, this routine is used to convert coordinates between two
743!> rotated-pole systems: COSMO-DE's rotated-pole system, and one centred at the
744!> PALM-4U domain centre. In this case, the PALM-4U system is thought of as the
745!> rotated-pole system and the routine is used to rotate back to COSMO-DE's
746!> system which is thought of as the geographical one.
747!>
748!> Input parameters:
749!> -----------------
750!> phir(:), lamr(: ): latitudes and longitudes of the rotated-pole grid
751!>
752!> phip, lamp: latitude and longitude of the rotated north pole
753!>
754!> gam: "angle between the north poles. If [gam] is not present, the other
755!>       system is the real geographical system." (original phiro2rot
756!>       description)
757!>
758!> Output parameters:
759!> ------------------
760!> phi(:,:), lam(:,:): geographical latitudes and logitudes
761!------------------------------------------------------------------------------!
762    SUBROUTINE rotate_to_cosmo(phir, lamr, phip, lamp, phi, lam, gam)
763       REAL(dp), INTENT(IN)  ::  phir(0:), lamr(0:), phip, lamp, gam
764       REAL(dp), INTENT(OUT) ::  phi(0:,0:), lam(0:,0:)
765
766       INTEGER ::  i, j
767       
768       IF ( SIZE(phi, 1) .NE. SIZE(lam, 1) .OR. &
769            SIZE(phi, 2) .NE. SIZE(lam, 2) )  THEN
770          PRINT *, "inifor: rotate_to_cosmo: Dimensions of phi and lambda do not match. Dimensions are:"
771          PRINT *, "inifor: rotate_to_cosmo: phi: ", SIZE(phi, 1), SIZE(phi, 2)
772          PRINT *, "inifor: rotate_to_cosmo: lam: ", SIZE(lam, 1), SIZE(lam, 2)
773          STOP
774       ENDIF
775
776       IF ( SIZE(phir) .NE. SIZE(phi, 2) .OR. &
777            SIZE(lamr) .NE. SIZE(phi, 1) )  THEN
778          PRINT *, "inifor: rotate_to_cosmo: Dimensions of phir and lamr do not match. Dimensions are:"
779          PRINT *, "inifor: rotate_to_cosmo: phir: ", SIZE(phir), SIZE(phi, 2)
780          PRINT *, "inifor: rotate_to_cosmo: lamr: ", SIZE(lamr), SIZE(phi, 1)
781          STOP
782       ENDIF
783       
784       DO j = 0, UBOUND(phir, 1)
785          DO i = 0, UBOUND(lamr, 1)
786
787             phi(i,j) = phirot2phi(phir(j) * TO_DEGREES,                       &
788                                   lamr(i) * TO_DEGREES,                       &
789                                   phip * TO_DEGREES,                          &
790                                   gam  * TO_DEGREES) * TO_RADIANS
791
792             lam(i,j) = rlarot2rla(phir(j) * TO_DEGREES,                       &
793                                   lamr(i) * TO_DEGREES,                       &
794                                   phip * TO_DEGREES,                          &
795                                   lamp * TO_DEGREES,                          &
796                                   gam  * TO_DEGREES) * TO_RADIANS
797
798          ENDDO
799       ENDDO
800
801    END SUBROUTINE rotate_to_cosmo
802       
803
804!------------------------------------------------------------------------------!
805! Description:
806! ------------
807!> Rotate the given vector field (x(:), y(:)) by the given 'angle'.
808!------------------------------------------------------------------------------!
809    SUBROUTINE rotate_vector_field(x, y, angle)
810       REAL(dp), DIMENSION(:), INTENT(INOUT) :: x, y  !< x and y coodrinate in arbitrary units
811       REAL(dp), INTENT(IN)                  :: angle !< rotation angle [deg]
812
813       INTEGER  :: i
814       REAL(dp) :: sine, cosine, v_rot(2), rotation(2,2)
815
816       sine = SIN(angle * TO_RADIANS)
817       cosine = COS(angle * TO_RADIANS)
818!
819!--    RESAHPE() fills columns first, so the rotation matrix becomes
820!--   
821!--    rotation = [ cosine   -sine  ]
822!--               [  sine    cosine ]
823       rotation = RESHAPE( (/cosine, sine, -sine, cosine/), (/2, 2/) )
824
825       DO i = LBOUND(x, 1), UBOUND(x, 1)
826
827          v_rot(:) = MATMUL(rotation, (/x(i), y(i)/))
828
829          x(i) = v_rot(1)
830          y(i) = v_rot(2)
831
832       ENDDO
833
834    END SUBROUTINE rotate_vector_field
835
836
837
838!------------------------------------------------------------------------------!
839! Description:
840! ------------
841!> This routine computes the local meridian convergence between a rotated-pole
842!> and a geographical system using the Eq. (6) given in the DWD manual
843!>
844!>    Baldauf et al. (2018), Beschreibung des operationelle Kürzestfrist-
845!>    vorhersagemodells COSMO-D2 und COSMO-D2-EPS und seiner Ausgabe in die
846!>    Datenbanken des DWD.
847!>    https://www.dwd.de/SharedDocs/downloads/DE/modelldokumentationen/nwv/cosmo_d2/cosmo_d2_dbbeschr_aktuell.pdf?__blob=publicationFile&v=2
848!------------------------------------------------------------------------------!
849    FUNCTION meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g)          &
850       RESULT(delta)
851
852       REAL(dp), INTENT(IN) ::  phi_n, lam_n, phi_g, lam_g
853       REAL(dp)             ::  delta
854
855       delta = atan2( COS(phi_n) * SIN(lam_n - lam_g),                         &
856                      COS(phi_g) * SIN(phi_n) -                                &
857                      SIN(phi_g) * COS(phi_n) * COS(lam_n - lam_g) )
858
859    END FUNCTION meridian_convergence_rotated
860
861!------------------------------------------------------------------------------!
862! Description:
863! ------------
864!> Compute indices of PALM-4U grid point neighbours in the target
865!> system (COSMO-DE) by rounding up and down. (i,j) are the indices of
866!> the PALM-4U grid and (ii(i,j,1-4), jj(i,j,1-4)) contain the indices
867!> of the its four neigbouring points in the COSMO-DE grid.
868!>
869!>
870!>                     COSMO-DE grid
871!>                     -------------
872!>           jj, lat
873!>              ^        j
874!>              |         \          i
875!>  jj(i,j,2/3) + ... 2 ---\--------/------ 3
876!>              |     | ^   \      /        |
877!>              |     | |wp  \    /         |
878!>              |     | v     \  /          |
879!>       latpos + ............ o/ (i,j)     |
880!>              |     |        :            |
881!>              |     |        :<----wl---->|
882!>  jj(i,j,1/4) + ... 1 -------:----------- 4
883!>              |     :        :            :
884!>              |     :        :            :
885!>              |     :      lonpos         :
886!>              L-----+--------+------------+------> ii, lon
887!>               ii(i,j,1/2)        ii(i,j,3/4)
888!>
889!>
890!> Input parameters:
891!> -----------------
892!> source_lat, source_lon : (rotated-pole) coordinates of the source grid (e.g.
893!>    COSMO-DE)
894!>
895!> source_dxi, source_dyi : inverse grid spacings of the source grid.
896!>
897!> target_lat, target_lon : (rotated-pole) coordinates of the target grid (e.g.
898!>    COSMO-DE)
899!>
900!> Output parameters:
901!> ------------------
902!> palm_ii, palm_jj : x and y index arrays of horizontal neighbour columns
903!>
904!------------------------------------------------------------------------------!
905    SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon,                &
906                                          palm_clat, palm_clon,                &
907                                          palm_ii, palm_jj)
908
909       REAL(dp), DIMENSION(0:), INTENT(IN)        ::  cosmo_lat, cosmo_lon
910       REAL(dp), DIMENSION(0:,0:), INTENT(IN)     ::  palm_clat, palm_clon
911       REAL(dp)                                   ::  cosmo_dxi, cosmo_dyi
912       INTEGER, DIMENSION(0:,0:,1:), INTENT(OUT)  ::  palm_ii, palm_jj
913
914       REAL(dp) ::  lonpos, latpos, lon0, lat0
915       INTEGER  ::  i, j
916
917       lon0 = cosmo_lon(0)
918       lat0 = cosmo_lat(0)
919       cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))
920       cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))
921
922       DO j = 0, UBOUND(palm_clon, 2)!palm_grid % ny
923       DO i = 0, UBOUND(palm_clon, 1)!palm_grid % nx
924!
925!--       Compute the floating point index corrseponding to PALM-4U grid point
926!--       location along target grid (COSMO-DE) axes.
927          lonpos = (palm_clon(i,j) - lon0) * cosmo_dxi
928          latpos = (palm_clat(i,j) - lat0) * cosmo_dyi
929
930          IF (lonpos < 0.0_dp .OR. latpos < 0.0_dp)  THEN
931             message = "lonpos or latpos out of bounds " //                    &
932                "while finding interpolation neighbours!" // NEW_LINE(' ') //  &
933                "          (i,j) = (" //                                       &
934                TRIM(str(i)) // ", " // TRIM(str(j)) // ")" // NEW_LINE(' ') //&
935                "          lonpos " // TRIM(real_to_str(lonpos*TO_DEGREES)) // &
936                ", latpos " // TRIM(real_to_str(latpos*TO_DEGREES)) // NEW_LINE(' ') // &
937                "          lon0 " // TRIM(real_to_str(lon0  *TO_DEGREES)) //   &
938                ", lat0   " // TRIM(real_to_str(lat0*TO_DEGREES)) // NEW_LINE(' ') // &
939                "          PALM lon " // TRIM(real_to_str(palm_clon(i,j)*TO_DEGREES)) // &
940                ", PALM lat " // TRIM(real_to_str(palm_clat(i,j)*TO_DEGREES))
941             CALL inifor_abort('find_horizontal_neighbours', message)
942          ENDIF
943
944          palm_ii(i,j,1) = FLOOR(lonpos)
945          palm_ii(i,j,2) = FLOOR(lonpos)
946          palm_ii(i,j,3) = CEILING(lonpos)
947          palm_ii(i,j,4) = CEILING(lonpos)
948
949          palm_jj(i,j,1) = FLOOR(latpos)
950          palm_jj(i,j,2) = CEILING(latpos)
951          palm_jj(i,j,3) = CEILING(latpos)
952          palm_jj(i,j,4) = FLOOR(latpos)
953       ENDDO
954       ENDDO
955
956    END SUBROUTINE find_horizontal_neighbours
957
958   
959!------------------------------------------------------------------------------!
960! Description:
961! ------------
962!> Computes linear vertical interpolation neighbour indices and weights for each
963!> column of the given palm grid.
964!------------------------------------------------------------------------------!
965    SUBROUTINE find_vertical_neighbours_and_weights_interp( palm_grid,         &
966                                                            palm_intermediate )
967       TYPE(grid_definition), INTENT(INOUT) ::  palm_grid
968       TYPE(grid_definition), INTENT(IN)    ::  palm_intermediate
969
970       INTEGER  ::  i, j, k, nx, ny, nz, nlev, k_intermediate
971       LOGICAL  ::  point_is_below_grid, point_is_above_grid,                  &
972                    point_is_in_current_cell
973       REAL(dp) ::  current_height, column_base, column_top, h_top, h_bottom,  &
974                    weight
975
976       nx   = palm_grid % nx
977       ny   = palm_grid % ny
978       nz   = palm_grid % nz
979       nlev = palm_intermediate % nz
980
981!
982!--    in each column of the fine grid, find vertical neighbours of every cell
983       DO j = 0, ny
984       DO i = 0, nx
985
986          k_intermediate = 0
987
988          column_base = palm_intermediate % h(i,j,0)
989          column_top  = palm_intermediate % h(i,j,nlev)
990
991!
992!--       scan through palm_grid column and set neighbour indices in
993!--       case current_height is either below column_base, in the current
994!--       cell, or above column_top. Keep increasing current cell index until
995!--       the current cell overlaps with the current_height.
996          DO k = 1, nz
997
998!
999!--          Memorize the top and bottom boundaries of the coarse cell and the
1000!--          current height within it
1001             current_height = palm_grid % z(k) + palm_grid % z0
1002             h_top    = palm_intermediate % h(i,j,k_intermediate+1)
1003             h_bottom = palm_intermediate % h(i,j,k_intermediate)
1004
1005             point_is_above_grid = (current_height > column_top) !22000m, very unlikely
1006             point_is_below_grid = (current_height < column_base)
1007
1008             point_is_in_current_cell = (                                      &
1009                current_height >= h_bottom .AND.                               &
1010                current_height <  h_top                                        &
1011             )
1012
1013!
1014!--          set default weights
1015             palm_grid % w_verti(i,j,k,1:2) = 0.0_dp
1016
1017             IF (point_is_above_grid)  THEN
1018
1019                palm_grid % kk(i,j,k,1:2) = nlev
1020                palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp
1021
1022                message = "PALM-4U grid extends above COSMO-DE model top."
1023                CALL inifor_abort('find_vertical_neighbours_and_weights', message)
1024
1025             ELSE IF (point_is_below_grid)  THEN
1026
1027                palm_grid % kk(i,j,k,1:2) = 0
1028                palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp
1029
1030             ELSE
1031!
1032!--             cycle through intermediate levels until current
1033!--             intermediate-grid cell overlaps with current_height
1034                DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
1035                   k_intermediate = k_intermediate + 1
1036
1037                   h_top    = palm_intermediate % h(i,j,k_intermediate+1)
1038                   h_bottom = palm_intermediate % h(i,j,k_intermediate)
1039                   point_is_in_current_cell = (                                &
1040                      current_height >= h_bottom .AND.                         &
1041                      current_height <  h_top                                  &
1042                   )
1043                ENDDO
1044
1045                IF (k_intermediate > nlev-1)  THEN
1046                   message = "Index " // TRIM(str(k_intermediate)) //          &
1047                             " is above intermediate grid range."
1048                   CALL inifor_abort('find_vertical_neighbours', message)
1049                ENDIF
1050   
1051                palm_grid % kk(i,j,k,1) = k_intermediate
1052                palm_grid % kk(i,j,k,2) = k_intermediate + 1
1053
1054!
1055!--             compute vertical weights
1056                weight = (h_top - current_height) / (h_top - h_bottom)
1057                palm_grid % w_verti(i,j,k,1) = weight
1058                palm_grid % w_verti(i,j,k,2) = 1.0_dp - weight
1059             ENDIF
1060
1061          ENDDO
1062
1063       ENDDO
1064       ENDDO
1065
1066    END SUBROUTINE find_vertical_neighbours_and_weights_interp
1067
1068
1069!------------------------------------------------------------------------------!
1070! Description:
1071! ------------
1072!> Computes linear vertical interpolation neighbour indices and weights for each
1073!> column of the given averaging grid.
1074!>
1075!> The difference to the _interp variant of this routine lies in how columns
1076!> are adressed. While the _interp variant loops over all PALM grid columns
1077!> given by combinations of all index indices (i,j), this routine loops over a
1078!> subset of COSMO columns, which are sequantlially stored in the index lists
1079!> iii(:) and jjj(:).
1080!------------------------------------------------------------------------------!
1081    SUBROUTINE find_vertical_neighbours_and_weights_average(                   &
1082       avg_grid, level_based_averaging                                         &
1083    )
1084
1085       TYPE(grid_definition), INTENT(INOUT), TARGET ::  avg_grid
1086       LOGICAL                                      ::  level_based_averaging
1087
1088       INTEGER           ::  i, j, k_palm, k_intermediate, l, nlev
1089       LOGICAL           ::  point_is_below_grid, point_is_above_grid,         &
1090                             point_is_in_current_cell
1091       REAL(dp)          ::  current_height, column_base, column_top, h_top,   &
1092                             h_bottom, weight
1093       REAL(dp), POINTER ::  cosmo_h(:,:,:)
1094
1095
1096       avg_grid % k_min = LBOUND(avg_grid % z, 1)
1097
1098       nlev = SIZE(avg_grid % cosmo_h, 3)
1099
1100       IF (level_based_averaging)  THEN
1101          cosmo_h => avg_grid % h
1102       ELSE
1103          cosmo_h => avg_grid % cosmo_h
1104       ENDIF
1105
1106!
1107!--    in each column of the fine grid, find vertical neighbours of every cell
1108       DO l = 1, avg_grid % n_columns
1109
1110          IF (level_based_averaging)  THEN
1111             i = 1
1112             j = 1
1113          ELSE
1114             i = avg_grid % iii(l)
1115             j = avg_grid % jjj(l)
1116          ENDIF
1117
1118          column_base = cosmo_h(i,j,1)
1119          column_top  = cosmo_h(i,j,nlev)
1120
1121!
1122!--       scan through avg_grid column until and set neighbour indices in
1123!--       case current_height is either below column_base, in the current
1124!--       cell, or above column_top. Keep increasing current cell index until
1125!--       the current cell overlaps with the current_height.
1126          k_intermediate = 1 !avg_grid % cosmo_h is indezed 1-based.
1127          DO k_palm = 1, avg_grid % nz
1128
1129!
1130!--          Memorize the top and bottom boundaries of the coarse cell and the
1131!--          current height within it
1132             current_height = avg_grid % z(k_palm) + avg_grid % z0
1133             h_top    = cosmo_h(i,j,k_intermediate+1)
1134             h_bottom = cosmo_h(i,j,k_intermediate)
1135
1136!
1137!--          COSMO column top is located at 22000m, point_is_above_grid is very
1138!--          unlikely.
1139             point_is_above_grid = (current_height > column_top)
1140             point_is_below_grid = (current_height < column_base)
1141
1142             point_is_in_current_cell = (                                      &
1143                current_height >= h_bottom .AND.                               &
1144                current_height <  h_top                                        &
1145             )
1146
1147!
1148!--          set default weights
1149             avg_grid % w(l,k_palm,1:2) = 0.0_dp
1150
1151             IF (point_is_above_grid)  THEN
1152
1153                avg_grid % kkk(l,k_palm,1:2) = nlev
1154                avg_grid % w(l,k_palm,1:2) = - 2.0_dp
1155
1156                message = "PALM-4U grid extends above COSMO-DE model top."
1157                CALL inifor_abort('find_vertical_neighbours_and_weights_average', message)
1158
1159             ELSE IF (point_is_below_grid)  THEN
1160
1161                avg_grid % kkk(l,k_palm,1:2) = 0
1162                avg_grid % w(l,k_palm,1:2) = - 2.0_dp
1163                avg_grid % k_min = MAX(k_palm + 1, avg_grid % k_min)
1164             ELSE
1165!
1166!--             cycle through intermediate levels until current
1167!--             intermediate-grid cell overlaps with current_height
1168                DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
1169                   k_intermediate = k_intermediate + 1
1170
1171                   h_top    = cosmo_h(i,j,k_intermediate+1)
1172                   h_bottom = cosmo_h(i,j,k_intermediate)
1173                   point_is_in_current_cell = (                                &
1174                      current_height >= h_bottom .AND.                         &
1175                      current_height <  h_top                                  &
1176                   )
1177                ENDDO
1178
1179!
1180!--             k_intermediate = 48 indicates the last section (indices 48 and 49), i.e.
1181!--             k_intermediate = 49 is not the beginning of a valid cell.
1182                IF (k_intermediate > nlev-1)  THEN
1183                   message = "Index " // TRIM(str(k_intermediate)) //          &
1184                             " is above intermediate grid range."
1185                   CALL inifor_abort('find_vertical_neighbours', message)
1186                ENDIF
1187   
1188                avg_grid % kkk(l,k_palm,1) = k_intermediate
1189                avg_grid % kkk(l,k_palm,2) = k_intermediate + 1
1190
1191!
1192!--             compute vertical weights
1193                weight = (h_top - current_height) / (h_top - h_bottom)
1194                avg_grid % w(l,k_palm,1) = weight
1195                avg_grid % w(l,k_palm,2) = 1.0_dp - weight
1196             ENDIF
1197
1198!
1199!--       Loop over PALM levels k
1200          ENDDO
1201
1202!
1203!--       Loop over averaging columns l
1204       ENDDO
1205 
1206    END SUBROUTINE find_vertical_neighbours_and_weights_average
1207
1208!------------------------------------------------------------------------------!
1209! Description:
1210! ------------
1211!> Compute the four weights for horizontal bilinear interpolation given the
1212!> coordinates clon(i,j) clat(i,j) of the PALM-4U grid in the COSMO-DE
1213!> rotated-pole grid and the neightbour indices ii(i,j,1-4) and jj(i,j,1-4).
1214!>
1215!> Input parameters:
1216!> -----------------
1217!> palm_grid % clon : longitudes of PALM-4U scalars (cell centres) in COSMO-DE's rotated-pole grid [rad]
1218!>
1219!> palm_grid % clat : latitudes of PALM-4U cell centres in COSMO-DE's rotated-pole grid [rad]
1220!>
1221!> cosmo_grid % lon : rotated-pole longitudes of scalars (cell centres) of the COSMO-DE grid [rad]
1222!>
1223!> cosmo_grid % lat : rotated-pole latitudes of scalars (cell centers) of the COSMO-DE grid [rad]
1224!>
1225!> cosmo_grid % dxi : inverse grid spacing in the first dimension [m^-1]
1226!>
1227!> cosmo_grid % dyi : inverse grid spacing in the second dimension [m^-1]
1228!>
1229!> Output parameters:
1230!> ------------------
1231!> palm_grid % w_horiz(:,:,1-4) : weights for bilinear horizontal interpolation
1232!
1233!                               COSMO-DE grid
1234!                               -------------
1235!                     jj, lat
1236!                        ^        j
1237!                        |         \          i
1238!            jj(i,j,2/3) + ... 2 ---\--------/------ 3
1239!                        |     | ^   \      /        |
1240!                        |     | |wp  \    /         |
1241!                        |     | v     \  /          |
1242!                 latpos + ............ o/ (i,j)     |
1243!                        |     |        :            |
1244!                        |     |        :<----wl---->|
1245!            jj(i,j,1/4) + ... 1 -------:----------- 4
1246!                        |     :        :            :
1247!                        |     :        :            :
1248!                        |     :      lonpos         :
1249!                        L-----+--------+------------+------> ii, lon
1250!                         ii(i,j,1/2)        ii(i,j,3/4)
1251!         
1252!------------------------------------------------------------------------------!
1253    SUBROUTINE compute_horizontal_interp_weights(cosmo_lat, cosmo_lon,         &
1254       palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)
1255       
1256       REAL(dp), DIMENSION(0:), INTENT(IN)        ::  cosmo_lat, cosmo_lon
1257       REAL(dp)                                   ::  cosmo_dxi, cosmo_dyi
1258       REAL(dp), DIMENSION(0:,0:), INTENT(IN)     ::  palm_clat, palm_clon
1259       INTEGER, DIMENSION(0:,0:,1:), INTENT(IN)   ::  palm_ii, palm_jj
1260
1261       REAL(dp), DIMENSION(0:,0:,1:), INTENT(OUT) ::  palm_w_horiz
1262
1263       REAL(dp) ::  wl, wp
1264       INTEGER  ::  i, j
1265
1266       cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0))
1267       cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0))
1268
1269       DO j = 0, UBOUND(palm_clon, 2)
1270       DO i = 0, UBOUND(palm_clon, 1)
1271     
1272!
1273!--       weight in lambda direction
1274          wl = ( cosmo_lon(palm_ii(i,j,4)) - palm_clon(i,j) ) * cosmo_dxi
1275
1276!
1277!--       weight in phi direction
1278          wp = ( cosmo_lat(palm_jj(i,j,2)) - palm_clat(i,j) ) * cosmo_dyi
1279
1280          IF (wl > 1.0_dp .OR. wl < 0.0_dp)  THEN
1281              message = "Horizontal weight wl = " // TRIM(real_to_str(wl)) //   &
1282                        " is out bounds."
1283              CALL inifor_abort('compute_horizontal_interp_weights', message)
1284          ENDIF
1285          IF (wp > 1.0_dp .OR. wp < 0.0_dp)  THEN
1286              message = "Horizontal weight wp = " // TRIM(real_to_str(wp)) //   &
1287                        " is out bounds."
1288              CALL inifor_abort('compute_horizontal_interp_weights', message)
1289          ENDIF
1290
1291          palm_w_horiz(i,j,1) = wl * wp
1292          palm_w_horiz(i,j,2) = wl * (1.0_dp - wp)
1293          palm_w_horiz(i,j,3) = (1.0_dp - wl) * (1.0_dp - wp)
1294          palm_w_horiz(i,j,4) = 1.0_dp - SUM( palm_w_horiz(i,j,1:3) )
1295
1296       ENDDO
1297       ENDDO
1298       
1299    END SUBROUTINE compute_horizontal_interp_weights
1300
1301
1302!------------------------------------------------------------------------------!
1303! Description:
1304! ------------
1305!> Interpolates u and v components of velocities located at cell faces to the
1306!> cell centres by averaging neighbouring values.
1307!>
1308!> This routine is designed to be used with COSMO-DE arrays where there are the
1309!> same number of grid points for scalars (centres) and velocities (faces). In
1310!> COSMO-DE the velocity points are staggared one half grid spaceing up-grid
1311!> which means the first centre point has to be omitted and is set to zero.
1312!------------------------------------------------------------------------------!
1313    SUBROUTINE centre_velocities(u_face, v_face, u_centre, v_centre)
1314       REAL(dp), DIMENSION(0:,0:,0:), INTENT(IN)  ::  u_face, v_face
1315       REAL(dp), DIMENSION(0:,0:,0:), INTENT(OUT) ::  u_centre, v_centre
1316       INTEGER ::  nx, ny
1317
1318       nx = UBOUND(u_face, 1)
1319       ny = UBOUND(u_face, 2)
1320
1321       u_centre(0,:,:)  = 0.0_dp
1322       u_centre(1:,:,:) = 0.5_dp * ( u_face(0:nx-1,:,:) + u_face(1:,:,:) )
1323
1324       v_centre(:,0,:)  = 0.0_dp
1325       v_centre(:,1:,:) = 0.5_dp * ( v_face(:,0:ny-1,:) + v_face(:,1:,:) )
1326    END SUBROUTINE centre_velocities
1327
1328
1329!------------------------------------------------------------------------------!
1330! Description:
1331! ------------
1332!> Compute the geographical latitude of a point given in rotated-pole cordinates
1333!------------------------------------------------------------------------------!
1334    FUNCTION phirot2phi (phirot, rlarot, polphi, polgam)
1335   
1336       REAL(dp), INTENT (IN) ::  polphi      !< latitude of the rotated north pole
1337       REAL(dp), INTENT (IN) ::  phirot      !< latitude in the rotated system
1338       REAL(dp), INTENT (IN) ::  rlarot      !< longitude in the rotated system
1339       REAL(dp), INTENT (IN) ::  polgam      !< angle between the north poles of the systems
1340
1341       REAL(dp)              ::  phirot2phi  !< latitude in the geographical system
1342       
1343       REAL(dp)              ::  zsinpol, zcospol, zphis, zrlas, zarg, zgam
1344   
1345       zsinpol = SIN(polphi * TO_RADIANS)
1346       zcospol = COS(polphi * TO_RADIANS)
1347       zphis   = phirot * TO_RADIANS
1348
1349       IF (rlarot > 180.0_dp)  THEN
1350          zrlas = rlarot - 360.0_dp
1351       ELSE
1352          zrlas = rlarot
1353       ENDIF
1354       zrlas = zrlas * TO_RADIANS
1355     
1356       IF (polgam /= 0.0_dp)  THEN
1357          zgam = polgam * TO_RADIANS
1358          zarg = zsinpol * SIN (zphis) +                                       &
1359                 zcospol * COS(zphis) * ( COS(zrlas) * COS(zgam) -             &
1360                                          SIN(zgam)  * SIN(zrlas) )
1361       ELSE
1362          zarg = zcospol * COS (zphis) * COS (zrlas) + zsinpol * SIN (zphis)
1363       ENDIF
1364     
1365       phirot2phi = ASIN (zarg) * TO_DEGREES
1366   
1367    END FUNCTION phirot2phi
1368
1369
1370!------------------------------------------------------------------------------!
1371! Description:
1372! ------------
1373!> Compute the geographical latitude of a point given in rotated-pole cordinates
1374!------------------------------------------------------------------------------!
1375    FUNCTION phi2phirot (phi, rla, polphi, pollam)
1376   
1377       REAL(dp), INTENT (IN) ::  polphi !< latitude of the rotated north pole
1378       REAL(dp), INTENT (IN) ::  pollam !< longitude of the rotated north pole
1379       REAL(dp), INTENT (IN) ::  phi    !< latitude in the geographical system
1380       REAL(dp), INTENT (IN) ::  rla    !< longitude in the geographical system
1381       
1382       REAL(dp) ::  phi2phirot          !< longitude in the rotated system
1383       
1384       REAL(dp) ::  zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1
1385       
1386       zsinpol = SIN(polphi * TO_RADIANS)
1387       zcospol = COS(polphi * TO_RADIANS)
1388       zlampol = pollam * TO_RADIANS
1389       zphi    = phi * TO_RADIANS
1390
1391       IF (rla > 180.0_dp)  THEN
1392          zrla1 = rla - 360.0_dp
1393       ELSE
1394          zrla1 = rla
1395       ENDIF
1396       zrla = zrla1 * TO_RADIANS
1397       
1398       zarg1 = SIN(zphi) * zsinpol
1399       zarg2 = COS(zphi) * zcospol * COS(zrla - zlampol)
1400       
1401       phi2phirot = ASIN(zarg1 + zarg2) * TO_DEGREES
1402   
1403    END FUNCTION phi2phirot
1404
1405
1406!------------------------------------------------------------------------------!
1407! Description:
1408! ------------
1409!> Compute the geographical longitude of a point given in rotated-pole cordinates
1410!------------------------------------------------------------------------------!
1411    FUNCTION rlarot2rla(phirot, rlarot, polphi, pollam, polgam)
1412   
1413       REAL(dp), INTENT (IN) ::  polphi !< latitude of the rotated north pole
1414       REAL(dp), INTENT (IN) ::  pollam !< longitude of the rotated north pole
1415       REAL(dp), INTENT (IN) ::  phirot !< latitude in the rotated system
1416       REAL(dp), INTENT (IN) ::  rlarot !< longitude in the rotated system
1417       REAL(dp), INTENT (IN) ::  polgam !< angle between the north poles of the systems
1418       
1419       REAL(dp) ::  rlarot2rla          !< latitude in the geographical system
1420       
1421       REAL(dp) ::  zsinpol, zcospol, zlampol, zphis, zrlas, zarg1, zarg2, zgam
1422       
1423       zsinpol = SIN(TO_RADIANS * polphi)
1424       zcospol = COS(TO_RADIANS * polphi)
1425       zlampol = TO_RADIANS * pollam
1426       zphis   = TO_RADIANS * phirot
1427
1428       IF (rlarot > 180.0_dp)  THEN
1429          zrlas = rlarot - 360.0_dp
1430       ELSE
1431          zrlas = rlarot
1432       ENDIF
1433       zrlas   = TO_RADIANS * zrlas
1434     
1435       IF (polgam /= 0.0_dp)  THEN
1436          zgam  = TO_RADIANS * polgam
1437          zarg1 = SIN(zlampol) * (zcospol * SIN(zphis) - zsinpol*COS(zphis) *  &
1438                  (COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) -        &
1439                  COS(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) +       &
1440                                                COS(zrlas) * SIN(zgam) )
1441       
1442          zarg2 = COS (zlampol) * (zcospol * SIN(zphis) - zsinpol*COS(zphis) * &
1443                  (COS(zrlas) * COS(zgam) - SIN(zrlas) * SIN(zgam)) ) +        &
1444                  SIN(zlampol) * COS(zphis) * ( SIN(zrlas) * COS(zgam) +       &
1445                                                COS(zrlas) * SIN(zgam) )
1446       ELSE
1447          zarg1   = SIN (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis)  +     &
1448                                      zcospol *              SIN(zphis)) -     &
1449                    COS (zlampol) *             SIN(zrlas) * COS(zphis)
1450          zarg2   = COS (zlampol) * (-zsinpol * COS(zrlas) * COS(zphis)  +     &
1451                                      zcospol *              SIN(zphis)) +     &
1452                    SIN (zlampol) *             SIN(zrlas) * COS(zphis)
1453       ENDIF
1454     
1455       IF (zarg2 == 0.0_dp)  zarg2 = 1.0E-20_dp
1456     
1457       rlarot2rla = ATAN2(zarg1,zarg2) * TO_DEGREES
1458       
1459    END FUNCTION rlarot2rla
1460
1461
1462!------------------------------------------------------------------------------!
1463! Description:
1464! ------------
1465!> Compute the rotated-pole longitude of a point given in geographical cordinates
1466!------------------------------------------------------------------------------!
1467    FUNCTION rla2rlarot ( phi, rla, polphi, pollam, polgam )
1468
1469       REAL(dp), INTENT (IN) ::  polphi !< latitude of the rotated north pole
1470       REAL(dp), INTENT (IN) ::  pollam !< longitude of the rotated north pole
1471       REAL(dp), INTENT (IN) ::  phi    !< latitude in geographical system
1472       REAL(dp), INTENT (IN) ::  rla    !< longitude in geographical system
1473       REAL(dp), INTENT (IN) ::  polgam !< angle between the north poles of the systems
1474       
1475       REAL (KIND=dp) ::  rla2rlarot    !< latitude in the the rotated system
1476       
1477       REAL (KIND=dp) ::  zsinpol, zcospol, zlampol, zphi, zrla, zarg1, zarg2, zrla1
1478       
1479       zsinpol = SIN(polphi * TO_RADIANS)
1480       zcospol = COS(polphi * TO_RADIANS)
1481       zlampol = pollam * TO_RADIANS
1482       zphi    = phi * TO_RADIANS
1483
1484       IF (rla > 180.0_dp)  THEN
1485          zrla1 = rla - 360.0_dp
1486       ELSE
1487          zrla1 = rla
1488       ENDIF
1489       zrla = zrla1 * TO_RADIANS
1490       
1491       zarg1 = - SIN (zrla-zlampol) * COS(zphi)
1492       zarg2 = - zsinpol * COS(zphi) * COS(zrla-zlampol) + zcospol * SIN(zphi)
1493       
1494       IF (zarg2 == 0.0_dp)  zarg2 = 1.0E-20_dp
1495       
1496       rla2rlarot = ATAN2 (zarg1,zarg2) * TO_DEGREES
1497       
1498       IF (polgam /= 0.0_dp )  THEN
1499          rla2rlarot = polgam + rla2rlarot
1500          IF (rla2rlarot > 180._dp)  rla2rlarot = rla2rlarot - 360.0_dp
1501       ENDIF
1502       
1503    END FUNCTION rla2rlarot
1504
1505
1506!------------------------------------------------------------------------------!
1507! Description:
1508! ------------
1509!> Rotate the given velocity vector (u,v) from the geographical to the
1510!> rotated-pole system
1511!------------------------------------------------------------------------------!
1512    SUBROUTINE uv2uvrot(u, v, rlat, rlon, pollat, pollon, urot, vrot)
1513   
1514       REAL(dp), INTENT (IN)  ::  u, v           !< wind components in the true geographical system
1515       REAL(dp), INTENT (IN)  ::  rlat, rlon     !< coordinates in the true geographical system
1516       REAL(dp), INTENT (IN)  ::  pollat, pollon !< latitude and longitude of the north pole of the rotated grid
1517       
1518       REAL(dp), INTENT (OUT) ::  urot, vrot     !< wind components in the rotated grid             
1519       
1520       REAL (dp) ::  zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm
1521       
1522       zsinpol = SIN(pollat * TO_RADIANS)
1523       zcospol = COS(pollat * TO_RADIANS)
1524       zlonp   = (pollon-rlon) * TO_RADIANS
1525       zlat    = rlat * TO_RADIANS
1526       
1527       zarg1 = zcospol * SIN(zlonp)
1528       zarg2 = zsinpol * COS(zlat) - zcospol * SIN(zlat) * COS(zlonp)
1529       znorm = 1.0_dp / SQRT(zarg1*zarg1 + zarg2*zarg2)
1530       
1531       urot = u * zarg2 * znorm - v * zarg1 * znorm
1532       vrot = u * zarg1 * znorm + v * zarg2 * znorm
1533   
1534    END SUBROUTINE uv2uvrot
1535
1536
1537!------------------------------------------------------------------------------!
1538! Description:
1539! ------------
1540!> Rotate the given velocity vector (urot, vrot) from the rotated-pole to the
1541!> geographical system
1542!------------------------------------------------------------------------------!
1543    SUBROUTINE uvrot2uv (urot, vrot, rlat, rlon, pollat, pollon, u, v)
1544   
1545       REAL(dp), INTENT(IN) ::  urot, vrot     !< wind components in the rotated grid
1546       REAL(dp), INTENT(IN) ::  rlat, rlon     !< latitude and longitude in the true geographical system
1547       REAL(dp), INTENT(IN) ::  pollat, pollon !< latitude and longitude of the north pole of the rotated grid
1548       
1549       REAL(dp), INTENT(OUT) ::  u, v          !< wind components in the true geographical system
1550       
1551       REAL(dp) ::  zsinpol, zcospol, zlonp, zlat, zarg1, zarg2, znorm
1552     
1553       zsinpol = SIN(pollat * TO_RADIANS)
1554       zcospol = COS(pollat * TO_RADIANS)
1555       zlonp   = (pollon-rlon) * TO_RADIANS
1556       zlat    = rlat * TO_RADIANS
1557     
1558       zarg1 = zcospol * SIN(zlonp)
1559       zarg2 = zsinpol * COS(zlat) - zcospol * SIN(zlat) * COS(zlonp)
1560       znorm = 1.0_dp / SQRT(zarg1*zarg1 + zarg2*zarg2)
1561     
1562       u =   urot * zarg2 * znorm + vrot * zarg1 * znorm
1563       v = - urot * zarg1 * znorm + vrot * zarg2 * znorm
1564   
1565    END SUBROUTINE uvrot2uv
1566
1567 END MODULE inifor_transform
1568#endif
1569
Note: See TracBrowser for help on using the repository browser.