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

Last change on this file since 4796 was 4756, checked in by eckhard, 4 years ago

inifor: Fixed an error in surface pressure extrapolation

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