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

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

Small error handling improvements

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