source: palm/trunk/SOURCE/surface_mod.f90 @ 4595

Last change on this file since 4595 was 4595, checked in by suehring, 10 months ago

Fix accidently commented subroutine

  • Property svn:keywords set to Id
File size: 270.4 KB
RevLine 
[3833]1!> @file surface_mod.f90
[4559]2!--------------------------------------------------------------------------------------------------!
[3833]3! This file is part of the PALM model system.
4!
[4559]5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
[3833]8!
[4559]9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
[3833]12!
[4559]13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
[3833]15!
[4360]16! Copyright 1997-2020 Leibniz Universitaet Hannover
[4559]17!--------------------------------------------------------------------------------------------------!
[3833]18!
19!
20! Current revisions:
[4559]21! -----------------
[4586]22!
23!
[3833]24! Former revisions:
25! -----------------
26! $Id: surface_mod.f90 4595 2020-07-09 17:18:21Z suehring $
[4595]27! Fix accidently commented subroutine
28!
29! 4594 2020-07-09 15:01:00Z suehring
[4594]30! Bugfix, add acc directives for scalar-roughness length
31!
32! 4593 2020-07-09 12:48:18Z suehring
[4593]33! Add arrays for pre-calculated ln(z/z0)
34!
35! 4586 2020-07-01 16:16:43Z gronemeier
[4586]36! renamed Richardson flux number into gradient Richardson number (1D model)
37!
38! 4559 2020-06-11 08:51:48Z raasch
[4559]39! File re-formatted to follow the PALM coding standard
40!
41! 4535 2020-05-15 12:07:23Z raasch
42! Bugfix for restart data format query
43!
[4535]44! 4521 2020-05-06 11:39:49Z schwenkel
[4521]45! Rename variable
[4559]46!
[4521]47! 4517 2020-05-03 14:29:30Z raasch
[4559]48! Added restart with MPI-IO for reading local arrays
49!
[4517]50! 4502 2020-04-17 16:14:16Z schwenkel
[4502]51! Implementation of ice microphysics
[4559]52!
[4502]53! 4495 2020-04-13 20:11:20Z raasch
[4559]54! Restart data handling with MPI-IO added
55!
[4495]56! 4366 2020-01-09 08:12:43Z raasch
[4559]57! Workaround implemented to avoid vectorization bug on NEC Aurora
58!
[4366]59! 4360 2020-01-07 11:25:50Z suehring
[4355]60! Fix also remaining message calls.
[4559]61!
[4355]62! 4354 2019-12-19 16:10:18Z suehring
[4354]63! Bugfix in message call and specify error number
[4559]64!
[4354]65! 4346 2019-12-18 11:55:56Z motisi
[4559]66! Introduction of wall_flags_total_0, which currently sets bits based on static topography
67! information used in wall_flags_static_0
68!
[4346]69! 4331 2019-12-10 18:25:02Z suehring
[4331]70! -pt_2m - array is moved to diagnostic_output_quantities
[4559]71!
[4331]72! 4329 2019-12-10 15:46:36Z motisi
[4329]73! Renamed wall_flags_0 to wall_flags_static_0
[4559]74!
[4329]75! 4245 2019-09-30 08:40:37Z pavelkrc
[4182]76! Corrected "Former revisions" section
[4559]77!
[4182]78! 4168 2019-08-16 13:50:17Z suehring
[4559]79! Remove functions get_topography_top_index. These are now replaced by precalculated arrays because
80! of too much CPU-time consumption
81!
[4168]82! 4159 2019-08-15 13:31:35Z suehring
[4159]83! Surface classification revised and adjusted to changes in init_grid
[4559]84!
[4159]85! 4156 2019-08-14 09:18:14Z schwenkel
[4156]86! Bugfix in case of cloud microphysics morrison
[4559]87!
[4156]88! 4150 2019-08-08 20:00:47Z suehring
[4150]89! Generic routine to initialize single surface properties added
[4559]90!
[4150]91! 4104 2019-07-17 17:08:20Z suehring
[4559]92! Bugfix, initialization of index space for boundary data structure accidantly run over ghost
93! points, causing a segmentation fault.
94!
[4104]95! 3943 2019-05-02 09:50:41Z maronga
[4102]96! - Revise initialization of the boundary data structure
97! - Add new data structure to set boundary conditions at vertical walls
[4559]98!
[4102]99! 3943 2019-05-02 09:50:41Z maronga
[3943]100! Removed qsws_eb as it is no longer needed.
[4559]101!
[3943]102! 3933 2019-04-25 12:33:20Z kanani
[3933]103! Add (de)allocation of pt_2m,
[4559]104! Bugfix: initialize pt_2m
105!
[3933]106! 3833 2019-03-28 15:04:04Z forkel
[4559]107! Added USE chem_gasphase_mod (chem_modules will not transport nvar and nspec anymore)
108!
[3833]109! 3772 2019-02-28 15:51:57Z suehring
[4559]110! Small change in the todo's
111!
[3772]112! 3767 2019-02-27 08:18:02Z raasch
[4559]113! Unused variables removed from rrd-subroutine parameter list
114!
[3833]115! 3761 2019-02-25 15:31:42Z raasch
[4559]116! OpenACC directives added to avoid compiler warnings about unused variables, unused variable
117! removed
118!
[3833]119! 3745 2019-02-15 18:57:56Z suehring
120! +waste_heat
[4559]121!
[3833]122! 3744 2019-02-15 18:38:58Z suehring
123! OpenACC port for SPEC
[4559]124!
[4182]125! 2233 2017-05-30 18:08:54Z suehring
126! Initial revision
[3833]127!
[4182]128!
[3833]129! Description:
130! ------------
[4559]131!> Surface module defines derived data structures to treat surface-bounded grid cells. Three
132!> different types of surfaces are defined: default surfaces, natural surfaces, and urban surfaces.
133!> The module encompasses the allocation and initialization of surface arrays, and handles reading
134!> and writing restart data. In addition, a further derived data structure is defined, in order to
135!> set boundary conditions at surfaces.
136!> @todo For the moment, downward-facing surfaces are only classified as default type
137!> @todo Clean up urban-surface variables (some of them are not used any more)
138!> @todo Revise initialization of surface fluxes (especially for chemistry)
[3833]139!> @todo Get rid-off deallocation routines in restarts
[4559]140!--------------------------------------------------------------------------------------------------!
[3833]141 MODULE surface_mod
142
[4559]143    USE arrays_3d,                                                                                 &
144        ONLY:  heatflux_input_conversion,                                                          &
145               momentumflux_input_conversion,                                                      &
146               rho_air,                                                                            &
147               rho_air_zw,                                                                         &
148               zu,                                                                                 &
149               zw,                                                                                 &
150               waterflux_input_conversion
[3833]151
[4559]152    USE chem_gasphase_mod,                                                                         &
153        ONLY:  nvar,                                                                               &
154               spc_names
[3833]155
156    USE chem_modules
157
158    USE control_parameters
159
[4559]160    USE indices,                                                                                   &
161        ONLY:  nxl,                                                                                &
162               nxlg,                                                                               &
163               nxr,                                                                                &
164               nxrg,                                                                               &
165               nys,                                                                                &
166               nysg,                                                                               &
167               nyn,                                                                                &
168               nyng,                                                                               &
169               nzb,                                                                                &
170               nzt,                                                                                &
[4346]171               wall_flags_total_0
[3833]172
[4559]173    USE grid_variables,                                                                            &
174        ONLY:  dx,                                                                                 &
175               dy
[3833]176
177    USE kinds
178
[4559]179    USE model_1d_mod,                                                                              &
[4586]180        ONLY:  ri1d,                                                                               &
[4559]181               us1d,                                                                               &
182               usws1d,                                                                             &
183               vsws1d
[3833]184
[4495]185    USE restart_data_mpi_io_mod,                                                                   &
[4559]186        ONLY:  rd_mpi_io_surface_filetypes,                                                        &
187               rrd_mpi_io,                                                                         &
188               rrd_mpi_io_global_array,                                                            &
189               rrd_mpi_io_surface,                                                                 &
190               total_number_of_surface_values,                                                     &
191               wrd_mpi_io,                                                                         &
192               wrd_mpi_io_global_array,                                                            &
193               wrd_mpi_io_surface
[3833]194
195    IMPLICIT NONE
196
197!
[4559]198!-- Data type used to identify grid-points where horizontal boundary conditions are applied
[3833]199    TYPE bc_type
[4559]200       INTEGER(iwp) ::  ioff  !< offset value in x-direction, used to determine index of surface element
201       INTEGER(iwp) ::  joff  !< offset value in y-direction, used to determine index of surface element
202       INTEGER(iwp) ::  koff  !< offset value in z-direction, used to determine index of surface element
203       INTEGER(iwp) ::  ns    !< number of surface elements on the PE
[3833]204
[4559]205       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i  !< x-index linking to the PALM 3D-grid
206       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j  !< y-index linking to the PALM 3D-grid
207       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k  !< z-index linking to the PALM 3D-grid
[3833]208
[4559]209       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  end_index    !< end index within surface data type for given (j,i)
210       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  start_index  !< start index within surface data type for given (j,i)
[3833]211
212    END TYPE bc_type
213!
[4502]214!-- Data type used to identify and treat surface-bounded grid points
[3833]215    TYPE surf_type
216
[4559]217       INTEGER(iwp) ::  ioff  !< offset value in x-direction, used to determine index of surface element
218       INTEGER(iwp) ::  joff  !< offset value in y-direction, used to determine index of surface element
219       INTEGER(iwp) ::  koff  !< offset value in z-direction, used to determine index of surface element
220       INTEGER(iwp) ::  ns    !< number of surface elements on the PE
[4502]221
[4559]222       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i  !< x-index linking to the PALM 3D-grid
223       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j  !< y-index linking to the PALM 3D-grid
224       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k  !< z-index linking to the PALM 3D-grid
[3833]225
[4502]226       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  facing  !< Bit indicating surface orientation
[3833]227
[4559]228       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  start_index  !< Start index within surface data type for given (j,i)
229       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  end_index    !< End index within surface data type for given (j,i)
[4502]230
[4559]231       LOGICAL ::  albedo_from_ascii = .FALSE. !< flag indicating that albedo for urban surfaces is input via ASCII format
232                                               !< (just for a workaround)
[3833]233
[4559]234       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z_mo     !< surface-layer height
235       REAL(wp), DIMENSION(:), ALLOCATABLE ::  uvw_abs  !< absolute surface-parallel velocity
236       REAL(wp), DIMENSION(:), ALLOCATABLE ::  us       !< friction velocity
237       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ts       !< scaling parameter temerature
238       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qs       !< scaling parameter humidity
239       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ss       !< scaling parameter passive scalar
240       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qcs      !< scaling parameter qc
241       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ncs      !< scaling parameter nc
242       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qis      !< scaling parameter qi
243       REAL(wp), DIMENSION(:), ALLOCATABLE ::  nis      !< scaling parameter ni
244       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qrs      !< scaling parameter qr
245       REAL(wp), DIMENSION(:), ALLOCATABLE ::  nrs      !< scaling parameter nr
[3833]246
[4559]247       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ol       !< Obukhov length
248       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rib      !< Richardson bulk number
[3833]249
[4559]250       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z0       !< roughness length for momentum
251       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z0h      !< roughness length for heat
252       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z0q      !< roughness length for humidity
[4502]253
[4559]254       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt1      !< potential temperature at first grid level
255       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qv1      !< mixing ratio at first grid level
256       REAL(wp), DIMENSION(:), ALLOCATABLE ::  vpt1     !< virtual potential temperature at first grid level
257
258       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  css    !< scaling parameter chemical species
[3833]259!
[4593]260!--    Pre-defined arrays for ln(z/z0)
261       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ln_z_z0  !< ln(z/z0)
262       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ln_z_z0h !< ln(z/z0h)
263       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ln_z_z0q !< ln(z/z0q)
264!
[3833]265!--    Define arrays for surface fluxes
[4559]266       REAL(wp), DIMENSION(:), ALLOCATABLE ::  usws     !< vertical momentum flux for u-component at horizontal surfaces
267       REAL(wp), DIMENSION(:), ALLOCATABLE ::  vsws     !< vertical momentum flux for v-component at horizontal surfaces
[3833]268
[4559]269       REAL(wp), DIMENSION(:), ALLOCATABLE ::  shf      !< surface flux sensible heat
270       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws     !< surface flux latent heat
271       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ssws     !< surface flux passive scalar
272       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qcsws    !< surface flux qc
273       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ncsws    !< surface flux nc
274       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qisws    !< surface flux qi
275       REAL(wp), DIMENSION(:), ALLOCATABLE ::  nisws    !< surface flux ni
276       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qrsws    !< surface flux qr
277       REAL(wp), DIMENSION(:), ALLOCATABLE ::  nrsws    !< surface flux nr
278       REAL(wp), DIMENSION(:), ALLOCATABLE ::  sasws    !< surface flux salinity
[3833]279!--    Added for SALSA:
[4559]280       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  answs  !< surface flux aerosol number: dim 1: flux, dim 2: bin
281       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  amsws  !< surface flux aerosol mass:   dim 1: flux, dim 2: bin
282       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gtsws  !< surface flux gesous tracers: dim 1: flux, dim 2: gas
[4502]283
[4559]284       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  cssws  !< surface flux chemical species
[3833]285!
286!--    Required for horizontal walls in production_e
[4559]287       REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_0  !< virtual velocity component (see production_e_init for further explanation)
288       REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_0  !< virtual velocity component (see production_e_init for further explanation)
[3833]289
[4559]290       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  mom_flux_uv   !< momentum flux usvs and vsus at vertical surfaces
291                                                               !< (used in diffusion_u and diffusion_v)
292       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  mom_flux_w    !< momentum flux wsus and wsvs at vertical surfaces
293                                                               !< (used in diffusion_w)
294       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  mom_flux_tke  !< momentum flux usvs, vsus, wsus, wsvs at vertical surfaces at grid
295                                                               !< center (used in production_e)
[3833]296!
297!--    Variables required for LSM as well as for USM
298       CHARACTER(LEN=40), DIMENSION(:), ALLOCATABLE ::  building_type_name    !< building type name at surface element
299       CHARACTER(LEN=40), DIMENSION(:), ALLOCATABLE ::  pavement_type_name    !< pavement type name at surface element
300       CHARACTER(LEN=40), DIMENSION(:), ALLOCATABLE ::  vegetation_type_name  !< water type at name surface element
301       CHARACTER(LEN=40), DIMENSION(:), ALLOCATABLE ::  water_type_name       !< water type at name surface element
[4502]302
[3833]303       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nzt_pavement     !< top index for pavement in soil
304       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  building_type    !< building type at surface element
305       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  pavement_type    !< pavement type at surface element
306       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  vegetation_type  !< vegetation type at surface element
307       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  water_type       !< water type at surface element
[4502]308
[4559]309       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  albedo_type  !< albedo type, for each fraction
310                                                                  !< (wall,green,window or vegetation,pavement water)
[3833]311
[4559]312       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  building_surface  !< flag parameter indicating that the surface element is covered
313                                                                 !< by buildings (no LSM actions, not implemented yet)
314       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  building_covered  !< flag indicating that buildings are on top of orography,
315                                                                 !< only used for vertical surfaces in LSM
[3833]316       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  pavement_surface    !< flag parameter for pavements
317       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  water_surface       !< flag parameter for water surfaces
318       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  vegetation_surface  !< flag parameter for natural land surfaces
[4502]319
[4559]320       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  albedo  !< broadband albedo for each surface fraction
321                                                         !< (LSM: vegetation, water, pavement; USM: wall, green, window)
322       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  emissivity  !< emissivity of the surface, for each fraction
323                                                             !< (LSM: vegetation, water, pavement; USM: wall, green, window)
324       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  frac  !< relative surface fraction
325                                                       !< (LSM: vegetation, water, pavement; USM: wall, green, window)
[3833]326
[4559]327       REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  aldif       !< albedo for longwave diffusive radiation, solar angle of 60 degrees
328       REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  aldir       !< albedo for longwave direct radiation, solar angle of 60 degrees
329       REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  asdif       !< albedo for shortwave diffusive radiation, solar angle of 60 deg.
330       REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  asdir       !< albedo for shortwave direct radiation, solar angle of 60 degrees
331       REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  rrtm_aldif  !< albedo for longwave diffusive radiation, solar angle of 60 degrees
332       REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  rrtm_aldir  !< albedo for longwave direct radiation, solar angle of 60 degrees
333       REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  rrtm_asdif  !< albedo for shortwave diffusive radiation, solar angle of 60 deg.
334       REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  rrtm_asdir  !< albedo for shortwave direct radiation, solar angle of 60 degrees
[3833]335
[4559]336       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  q_surface        !< skin-surface mixing ratio
337       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  pt_surface       !< skin-surface temperature
338       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  vpt_surface      !< skin-surface virtual temperature
339       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rad_net          !< net radiation
340       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rad_net_l        !< net radiation, used in USM
341       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h         !< heat conductivity of soil/ wall (W/m/K)
342       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_green   !< heat conductivity of green soil (W/m/K)
343       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_window  !< heat conductivity of windows (W/m/K)
344       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_def     !< default heat conductivity of soil (W/m/K)
[3833]345
[4559]346       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_lw_in   !< incoming longwave radiation
347       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_lw_out  !< emitted longwave radiation
348       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_lw_dif  !< incoming longwave radiation from sky
349       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_lw_ref  !< incoming longwave radiation from reflection
350       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_lw_res  !< resedual longwave radiation in surface after last reflection step
351       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_sw_in   !< incoming shortwave radiation
352       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_sw_out  !< emitted shortwave radiation
353       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_sw_dir  !< direct incoming shortwave radiation
354       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_sw_dif  !< diffuse incoming shortwave radiation
355       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_sw_ref  !< incoming shortwave radiation from reflection
356       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_sw_res  !< resedual shortwave radiation in surface after last reflection step
[3833]357
[4559]358       REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_liq             !< liquid water coverage (of vegetated area)
359       REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_veg             !< vegetation coverage
360       REAL(wp), DIMENSION(:), ALLOCATABLE ::  f_sw_in           !< fraction of absorbed shortwave radiation by the surface layer
361                                                                 !< (not implemented yet)
362       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ghf               !< ground heat flux
363       REAL(wp), DIMENSION(:), ALLOCATABLE ::  g_d               !< coefficient for dependence of r_canopy
364                                                                 !< on water vapour pressure deficit
365       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lai               !< leaf area index
366       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surface_u  !< coupling between surface and soil (depends on vegetation type)
367                                                                 !< (W/m2/K)
368       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surface_s  !< coupling between surface and soil (depends on vegetation type)
369                                                                 !< (W/m2/K)
370       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_liq          !< surface flux of latent heat (liquid water portion)
371       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_soil         !< surface flux of latent heat (soil portion)
372       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_veg          !< surface flux of latent heat (vegetation portion)
[3833]373
[4559]374       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_a           !< aerodynamic resistance
375       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_a_green     !< aerodynamic resistance at green fraction
376       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_a_window    !< aerodynamic resistance at window fraction
377       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_canopy      !< canopy resistance
378       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_soil        !< soil resistance
379       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_soil_min    !< minimum soil resistance
380       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_s           !< total surface resistance (combination of r_soil and r_canopy)
381       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_canopy_min  !< minimum canopy (stomatal) resistance
[4502]382
[4559]383       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_10cm  !< near surface air potential temperature at distance 10 cm from
384                                                        !< the surface (K)
[4502]385
[4559]386       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  alpha_vg         !< coef. of Van Genuchten
387       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_w         !< hydraulic diffusivity of soil (?)
388       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gamma_w          !< hydraulic conductivity of soil (W/m/K)
389       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gamma_w_sat      !< hydraulic conductivity at saturation
390       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  l_vg             !< coef. of Van Genuchten
391       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_fc             !< soil moisture at field capacity (m3/m3)
392       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_res            !< residual soil moisture
393       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_sat            !< saturation soil moisture (m3/m3)
394       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_wilt           !< soil moisture at permanent wilting point (m3/m3)
395       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_vg             !< coef. Van Genuchten
396       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_total_def  !< default volumetric heat capacity of the (soil) layer (J/m3/K)
397       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_total      !< volumetric heat capacity of the actual soil matrix (J/m3/K)
398       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  root_fr          !< root fraction within the soil layers
[4502]399
[3833]400!--    Indoor model variables
[4559]401       REAL(wp), DIMENSION(:), ALLOCATABLE ::  waste_heat  !< waste heat
[3833]402!
403!--    Urban surface variables
[4559]404       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  surface_types  !< array of types of wall parameters
[3833]405
[4559]406       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  isroof_surf   !< flag indicating roof surfaces
407       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  ground_level  !< flag indicating ground floor level surfaces
[3833]408
409       REAL(wp), DIMENSION(:), ALLOCATABLE ::  target_temp_summer  !< indoor target temperature summer
[4502]410       REAL(wp), DIMENSION(:), ALLOCATABLE ::  target_temp_winter  !< indoor target temperature summer
[3833]411
412       REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_surface           !< heat capacity of the wall surface skin (J/m2/K)
413       REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_surface_green     !< heat capacity of the green surface skin (J/m2/K)
414       REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_surface_window    !< heat capacity of the window surface skin (J/m2/K)
415       REAL(wp), DIMENSION(:), ALLOCATABLE ::  green_type_roof     !< type of the green roof
416       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surf         !< heat conductivity between air and surface (W/m2/K)
417       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surf_green   !< heat conductivity between air and green surface (W/m2/K)
418       REAL(wp), DIMENSION(:), ALLOCATABLE ::  lambda_surf_window  !< heat conductivity between air and window surface (W/m2/K)
419       REAL(wp), DIMENSION(:), ALLOCATABLE ::  thickness_wall      !< thickness of the wall, roof and soil layers
420       REAL(wp), DIMENSION(:), ALLOCATABLE ::  thickness_green     !< thickness of the green wall, roof and soil layers
421       REAL(wp), DIMENSION(:), ALLOCATABLE ::  thickness_window    !< thickness of the window wall, roof and soil layers
422       REAL(wp), DIMENSION(:), ALLOCATABLE ::  transmissivity      !< transmissivity of windows
423
[4559]424       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutsl  !< reflected shortwave radiation for local surface in i-th reflection
425       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutll  !< reflected + emitted longwave radiation for local surface
426                                                          !< in i-th reflection
427       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfhf     !< total radiation flux incoming to minus outgoing from local surface
[3833]428
[4559]429       REAL(wp), DIMENSION(:), ALLOCATABLE ::  tt_surface_wall_m    !< surface temperature tendency (K)
430       REAL(wp), DIMENSION(:), ALLOCATABLE ::  tt_surface_window_m  !< window surface temperature tendency (K)
431       REAL(wp), DIMENSION(:), ALLOCATABLE ::  tt_surface_green_m   !< green surface temperature tendency (K)
432       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wshf                 !< kinematic wall heat flux of sensible heat
433                                                                    !< (actually no longer needed)
434       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wshf_eb              !< wall heat flux of sensible heat in wall normal direction
[3833]435
[4559]436       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb          !< wall ground heat flux
437       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb_window   !< window ground heat flux
438       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb_green    !< green ground heat flux
439       REAL(wp), DIMENSION(:), ALLOCATABLE ::  iwghf_eb         !< indoor wall ground heat flux
440       REAL(wp), DIMENSION(:), ALLOCATABLE ::  iwghf_eb_window  !< indoor window ground heat flux
[3833]441
[4559]442       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_lw_out_change_0  !<
[3833]443
[4559]444       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinsw   !< shortwave radiation falling to local surface including radiation
445                                                          !< from reflections
446       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutsw  !< total shortwave radiation outgoing from nonvirtual surfaces surfaces
447                                                          !< after all reflection
448       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlw   !< longwave radiation falling to local surface including radiation from
449                                                          !< reflections
450       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutlw  !< total longwave radiation outgoing from nonvirtual surfaces surfaces
451                                                          !< after all reflection
[4502]452
[4559]453       REAL(wp), DIMENSION(:), ALLOCATABLE ::  n_vg_green      !< vangenuchten parameters
454       REAL(wp), DIMENSION(:), ALLOCATABLE ::  alpha_vg_green  !< vangenuchten parameters
455       REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_vg_green      !< vangenuchten parameters
[3833]456
457
[4559]458       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_wall         !< volumetric heat capacity of the material ( J m-3 K-1 )
459                                                                    !< (= 2.19E6)
460       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_wall            !< wall grid spacing (center-center)
461       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_wall           !< 1/dz_wall
462       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_wall_stag       !< wall grid spacing (edge-edge)
463       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_wall_stag      !< 1/dz_wall_stag
464       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_wall_m          !< t_wall prognostic array
465       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw                 !< wall layer depths (m)
466       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_window       !< volumetric heat capacity of the window material ( J m-3 K-1 )
467                                                                    !< (= 2.19E6)
468       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_window          !< window grid spacing (center-center)
469       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_window         !< 1/dz_window
470       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_window_stag     !< window grid spacing (edge-edge)
471       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_window_stag    !< 1/dz_window_stag
472       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_window_m        !< t_window prognostic array
473       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw_window          !< window layer depths (m)
474       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_green        !< volumetric heat capacity of the green material ( J m-3 K-1 )
475                                                                    !< (= 2.19E6)
476       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_total_green  !< volumetric heat capacity of the moist green material
477                                                                    !< ( J m-3 K-1 ) (= 2.19E6)
478       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_green           !< green grid spacing (center-center)
479       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_green          !< 1/dz_green
480       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_green_stag      !< green grid spacing (edge-edge)
481       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_green_stag     !< 1/dz_green_stag
482       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_green_m         !< t_green prognostic array
483       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw_green           !< green layer depths (m)
[3833]484
[4559]485       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gamma_w_green_sat  !< hydraulic conductivity
486       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_w_green     !<
487       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gamma_w_green      !< hydraulic conductivity
488       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tswc_h_m           !<
[3833]489
490
[4559]491!-- Arrays for time averages
492       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_net_av          !< average of rad_net_l
493       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinsw_av         !< average of sw radiation falling to local surface including
494                                                                   !< radiation from reflections
495       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlw_av         !< average of lw radiation falling to local surface including
496                                                                   !< radiation from reflections
497       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinswdir_av      !< average of direct sw radiation falling to local surface
498       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinswdif_av      !< average of diffuse sw radiation from sky and model boundary
499                                                                   !< falling to local surface
500       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlwdif_av      !< average of diffuse lw radiation from sky and model boundary
501                                                                   !< falling to local surface
502       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinswref_av      !< average of sw radiation falling to surface from reflections
503       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlwref_av      !< average of lw radiation falling to surface from reflections
504       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutsw_av        !< average of total sw radiation outgoing from nonvirtual
505                                                                   !< surfaces surfaces after all reflection
506       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutlw_av        !< average of total lw radiation outgoing from nonvirtual
507                                                                   !< surfaces after all reflection
508       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfins_av          !< average of array of residua of sw radiation absorbed in
509                                                                   !< surface after last reflection
510       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinl_av          !< average of array of residua of lw radiation absorbed in
511                                                                   !< surface after last reflection
512       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfhf_av           !< average of total radiation flux incoming to minus outgoing
513                                                                   !< from local surface
514       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb_av          !< average of wghf_eb
515       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb_window_av   !< average of wghf_eb window
516       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb_green_av    !< average of wghf_eb window
517       REAL(wp), DIMENSION(:), ALLOCATABLE ::  iwghf_eb_av         !< indoor average of wghf_eb
518       REAL(wp), DIMENSION(:), ALLOCATABLE ::  iwghf_eb_window_av  !< indoor average of wghf_eb window
519       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wshf_eb_av          !< average of wshf_eb
520       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_av             !< average of qsws
521       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_veg_av         !< average of qsws_veg_eb
522       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_liq_av         !< average of qsws_liq_eb
523       REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_surf_wall_av      !< average of wall surface temperature (K)
524       REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_surf_av           !< average of wall surface temperature (K)
525       REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_surf_window_av    !< average of window surface temperature (K)
526       REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_surf_green_av     !< average of green wall surface temperature (K)
[3833]527
[4559]528       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_10cm_av  !< average of theta_10cm (K)
[4502]529
[4559]530       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  t_wall_av    !< Average of t_wall
531       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  t_window_av  !< Average of t_window
532       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  t_green_av   !< Average of t_green
533       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  swc_av       !< Average of swc
[3833]534
535    END TYPE surf_type
536
[4559]537    TYPE (bc_type), DIMENSION(0:1)  ::  bc_h  !< boundary condition data type, horizontal upward- and downward facing surfaces
538    TYPE (bc_type), DIMENSION(0:3)  ::  bc_v  !< boundary condition data type, vertical surfaces
[3833]539
540    TYPE (surf_type), DIMENSION(0:2), TARGET ::  surf_def_h  !< horizontal default surfaces (Up, Down, and Top)
541    TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_def_v  !< vertical default surfaces (North, South, East, West)
542    TYPE (surf_type)                , TARGET ::  surf_lsm_h  !< horizontal natural land surfaces, so far only upward-facing
543    TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_lsm_v  !< vertical land surfaces (North, South, East, West)
544    TYPE (surf_type)                , TARGET ::  surf_usm_h  !< horizontal urban surfaces, so far only upward-facing
545    TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_usm_v  !< vertical urban surfaces (North, South, East, West)
546
[4559]547    INTEGER(iwp), PARAMETER ::  ind_veg_wall  = 0  !< index for vegetation / wall-surface fraction, used for access of albedo,
548                                                   !< emissivity, etc., for each surface type
549    INTEGER(iwp), PARAMETER ::  ind_pav_green = 1  !< index for pavement / green-wall surface fraction, used for access of albedo,
550                                                   !< emissivity, etc., for each surface type
551    INTEGER(iwp), PARAMETER ::  ind_wat_win   = 2  !< index for water / window-surface fraction, used for access of albedo,
552                                                   !< emissivity, etc., for each surface type
[3833]553
[4559]554    INTEGER(iwp) ::  ns_h_on_file(0:2)  !< total number of horizontal surfaces with the same facing, required for writing
555                                        !< restart data
556    INTEGER(iwp) ::  ns_v_on_file(0:3)  !< total number of vertical surfaces with the same facing, required for writing restart data
[3833]557
[4559]558    LOGICAL ::  vertical_surfaces_exist     = .FALSE.  !< flag indicating that there are vertical urban/land surfaces
559                                                       !< in the domain (required to activiate RTM)
[3833]560
[4559]561    LOGICAL ::  surf_bulk_cloud_model       = .FALSE.  !< use cloud microphysics
562    LOGICAL ::  surf_microphysics_morrison  = .FALSE.  !< use 2-moment Morrison (add. prog. eq. for nc and qc)
563    LOGICAL ::  surf_microphysics_seifert   = .FALSE.  !< use 2-moment Seifert and Beheng scheme
564    LOGICAL ::  surf_microphysics_ice_phase = .FALSE.  !< use 2-moment Seifert and Beheng scheme
[3833]565
566
567    SAVE
568
569    PRIVATE
[4502]570
[3833]571    INTERFACE init_bc
572       MODULE PROCEDURE init_bc
573    END INTERFACE init_bc
574
[4150]575    INTERFACE init_single_surface_properties
576       MODULE PROCEDURE init_single_surface_properties
577    END INTERFACE init_single_surface_properties
[4502]578
[3833]579    INTERFACE init_surfaces
580       MODULE PROCEDURE init_surfaces
581    END INTERFACE init_surfaces
582
583    INTERFACE init_surface_arrays
584       MODULE PROCEDURE init_surface_arrays
585    END INTERFACE init_surface_arrays
586
587    INTERFACE surface_rrd_local
[4517]588       MODULE PROCEDURE surface_rrd_local_ftn
589       MODULE PROCEDURE surface_rrd_local_mpi
[3833]590    END INTERFACE surface_rrd_local
591
592    INTERFACE surface_wrd_local
593       MODULE PROCEDURE surface_wrd_local
594    END INTERFACE surface_wrd_local
595
596    INTERFACE surface_last_actions
597       MODULE PROCEDURE surface_last_actions
598    END INTERFACE surface_last_actions
599
600    INTERFACE surface_restore_elements
601       MODULE PROCEDURE surface_restore_elements_1d
602       MODULE PROCEDURE surface_restore_elements_2d
603    END INTERFACE surface_restore_elements
604
605#if defined( _OPENACC )
606    INTERFACE enter_surface_arrays
607       MODULE PROCEDURE enter_surface_arrays
608    END INTERFACE
609
610    INTERFACE exit_surface_arrays
611       MODULE PROCEDURE exit_surface_arrays
612    END INTERFACE
613#endif
614
615!
616!-- Public variables
[4559]617    PUBLIC bc_h,                                                                                   &
618           bc_v,                                                                                   &
619           ind_pav_green,                                                                          &
620           ind_veg_wall,                                                                           &
621           ind_wat_win,                                                                            &
622           ns_h_on_file,                                                                           &
623           ns_v_on_file,                                                                           &
624           surf_def_h,                                                                             &
625           surf_def_v,                                                                             &
626           surf_lsm_h,                                                                             &
627           surf_lsm_v,                                                                             &
628           surf_usm_h,                                                                             &
629           surf_usm_v,                                                                             &
630           surf_type,                                                                              &
631           vertical_surfaces_exist,                                                                &
632           surf_bulk_cloud_model,                                                                  &
633           surf_microphysics_morrison,                                                             &
634           surf_microphysics_seifert,                                                              &
635           surf_microphysics_ice_phase
[3833]636!
637!-- Public subroutines and functions
[4168]638    PUBLIC init_bc,                                                                                &
[4150]639           init_single_surface_properties,                                                         &
640           init_surfaces,                                                                          &
641           init_surface_arrays,                                                                    &
642           surface_last_actions,                                                                   &
643           surface_rrd_local,                                                                      &
644           surface_restore_elements,                                                               &
645           surface_wrd_local
[3833]646
647#if defined( _OPENACC )
[4168]648    PUBLIC enter_surface_arrays,                                                                   &
649           exit_surface_arrays
[3833]650#endif
651
652 CONTAINS
653
[4559]654!--------------------------------------------------------------------------------------------------!
[3833]655! Description:
656! ------------
[4559]657!> Initialize data type for setting boundary conditions at horizontal and vertical surfaces.
658!--------------------------------------------------------------------------------------------------!
659 SUBROUTINE init_bc
[3833]660
[4559]661    IMPLICIT NONE
[3833]662
[4559]663    INTEGER(iwp) ::  i  !< loop index along x-direction
664    INTEGER(iwp) ::  j  !< loop index along y-direction
665    INTEGER(iwp) ::  k  !< loop index along y-direction
666    INTEGER(iwp) ::  l  !< running index for differently aligned surfaces
[3833]667
[4559]668    INTEGER(iwp), DIMENSION(0:1) ::  num_h          !< number of horizontal surfaces on subdomain
669    INTEGER(iwp), DIMENSION(0:1) ::  num_h_kji      !< number of horizontal surfaces at (j,i)-grid point
670    INTEGER(iwp), DIMENSION(0:1) ::  start_index_h  !< local start index of horizontal surface elements
[4502]671
[4559]672    INTEGER(iwp), DIMENSION(0:3) ::  num_v          !< number of vertical surfaces on subdomain
673    INTEGER(iwp), DIMENSION(0:3) ::  num_v_kji      !< number of vertical surfaces at (j,i)-grid point
674    INTEGER(iwp), DIMENSION(0:3) ::  start_index_v  !< local start index of vertical surface elements
[3833]675!
[4559]676!-- Set offset indices, i.e. index difference between surface element and surface-bounded grid point.
677!-- Horizontal surfaces - no horizontal offsets
678    bc_h(:)%ioff = 0
679    bc_h(:)%joff = 0
[3833]680!
[4559]681!-- Horizontal surfaces, upward facing (0) and downward facing (1)
682    bc_h(0)%koff = -1
683    bc_h(1)%koff = 1
[3833]684!
[4559]685!-- Vertical surfaces - no vertical offset
686    bc_v(0:3)%koff = 0
[3833]687!
[4559]688!-- North- and southward facing - no offset in x
689    bc_v(0:1)%ioff = 0
[3833]690!
[4559]691!-- Northward facing offset in y
692    bc_v(0)%joff = -1
[3833]693!
[4559]694!-- Southward facing offset in y
695    bc_v(1)%joff = 1
[3833]696!
[4559]697!-- East- and westward facing - no offset in y
698    bc_v(2:3)%joff = 0
[3833]699!
[4559]700!-- Eastward facing offset in x
701    bc_v(2)%ioff = -1
[3833]702!
[4559]703!-- Westward facing offset in y
704    bc_v(3)%ioff = 1
[3833]705!
[4559]706!-- Initialize data structure for horizontal surfaces, i.e. count the number of surface elements,
707!-- allocate and initialize the respective index arrays, and set the respective start and end
708!-- indices at each (j,i)-location. The index space is defined also over the ghost points, so that
709!-- e.g. boundary conditions for diagnostic quanitities can be set on ghost points so that no
710!-- exchange is required any more.
711    DO  l = 0, 1
[4102]712!
[4559]713!--    Count the number of upward- and downward-facing surfaces on subdomain
714       num_h(l) = 0
715       DO  i = nxlg, nxrg
716          DO  j = nysg, nyng
717             DO  k = nzb+1, nzt
[4502]718!
[4559]719!--             Check if current gridpoint belongs to the atmosphere
720                IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
721                   IF ( .NOT. BTEST( wall_flags_total_0(k+bc_h(l)%koff, j+bc_h(l)%joff,            &
722                                     i+bc_h(l)%ioff), 0 ) )  num_h(l) = num_h(l) + 1
723                ENDIF
[4102]724             ENDDO
725          ENDDO
[4559]726       ENDDO
[4502]727!
[4559]728!--    Save the number of horizontal surface elements
729       bc_h(l)%ns = num_h(l)
[3833]730!
[4559]731!--    ALLOCATE arrays for horizontal surfaces
732       ALLOCATE( bc_h(l)%i(1:bc_h(l)%ns) )
733       ALLOCATE( bc_h(l)%j(1:bc_h(l)%ns) )
734       ALLOCATE( bc_h(l)%k(1:bc_h(l)%ns) )
735       ALLOCATE( bc_h(l)%start_index(nysg:nyng,nxlg:nxrg) )
736       ALLOCATE( bc_h(l)%end_index(nysg:nyng,nxlg:nxrg) )
737       bc_h(l)%start_index = 1
738       bc_h(l)%end_index   = 0
[4502]739
[4559]740       num_h(l)         = 1
741       start_index_h(l) = 1
742       DO  i = nxlg, nxrg
743          DO  j = nysg, nyng
[4502]744
[4559]745             num_h_kji(l) = 0
746             DO  k = nzb+1, nzt
[4502]747!
[4559]748!--             Check if current gridpoint belongs to the atmosphere
749                IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
[4502]750!
[4559]751!--                Upward-facing
752                   IF ( .NOT. BTEST( wall_flags_total_0(k+bc_h(l)%koff, j+bc_h(l)%joff,            &
753                                     i+bc_h(l)%ioff), 0 ) )  THEN
754                      bc_h(l)%i(num_h(l)) = i
755                      bc_h(l)%j(num_h(l)) = j
756                      bc_h(l)%k(num_h(l)) = k
757                      num_h_kji(l)        = num_h_kji(l) + 1
758                      num_h(l)            = num_h(l) + 1
[3833]759                   ENDIF
[4559]760                ENDIF
[3833]761             ENDDO
[4559]762             bc_h(l)%start_index(j,i) = start_index_h(l)
763             bc_h(l)%end_index(j,i)   = bc_h(l)%start_index(j,i) + num_h_kji(l) - 1
764             start_index_h(l)         = bc_h(l)%end_index(j,i) + 1
[4102]765          ENDDO
766       ENDDO
[4559]767    ENDDO
[3833]768
[4102]769!
[4559]770!-- Initialize data structure for vertical surfaces, i.e. count the number of surface elements,
771!-- allocate and initialize the respective index arrays, and set the respective start and end
772!-- indices at each (j,i)-location.
773    DO  l = 0, 3
[4102]774!
[4559]775!--    Count the number of upward- and downward-facing surfaces on subdomain
776       num_v(l) = 0
777       DO  i = nxl, nxr
778          DO  j = nys, nyn
779             DO  k = nzb+1, nzt
[4502]780!
[4559]781!--             Check if current gridpoint belongs to the atmosphere
782                IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
783                   IF ( .NOT. BTEST( wall_flags_total_0(k+bc_v(l)%koff, j+bc_v(l)%joff,            &
784                                     i+bc_v(l)%ioff), 0 ) )  num_v(l) = num_v(l) + 1
785                ENDIF
[4102]786             ENDDO
[3833]787          ENDDO
[4559]788       ENDDO
[4502]789!
[4559]790!--    Save the number of horizontal surface elements
791       bc_v(l)%ns = num_v(l)
[4102]792!
[4559]793!--    ALLOCATE arrays for horizontal surfaces. In contrast to the horizontal surfaces, the index
794!--    space is not defined over the ghost points.
795       ALLOCATE( bc_v(l)%i(1:bc_v(l)%ns) )
796       ALLOCATE( bc_v(l)%j(1:bc_v(l)%ns) )
797       ALLOCATE( bc_v(l)%k(1:bc_v(l)%ns) )
798       ALLOCATE( bc_v(l)%start_index(nys:nyn,nxl:nxr) )
799       ALLOCATE( bc_v(l)%end_index(nys:nyn,nxl:nxr) )
800       bc_v(l)%start_index = 1
801       bc_v(l)%end_index   = 0
[4502]802
[4559]803       num_v(l)         = 1
804       start_index_v(l) = 1
805       DO  i = nxl, nxr
806          DO  j = nys, nyn
[4502]807
[4559]808             num_v_kji(l) = 0
809             DO  k = nzb+1, nzt
[4502]810!
[4559]811!--             Check if current gridpoint belongs to the atmosphere
812                IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
[4502]813!
[4559]814!--                Upward-facing
815                   IF ( .NOT. BTEST( wall_flags_total_0(k+bc_v(l)%koff, j+bc_v(l)%joff,            &
816                                     i+bc_v(l)%ioff), 0 ) )  THEN
817                      bc_v(l)%i(num_v(l)) = i
818                      bc_v(l)%j(num_v(l)) = j
819                      bc_v(l)%k(num_v(l)) = k
820                      num_v_kji(l)        = num_v_kji(l) + 1
821                      num_v(l)            = num_v(l) + 1
[4102]822                   ENDIF
[4559]823                ENDIF
[4102]824             ENDDO
[4559]825             bc_v(l)%start_index(j,i) = start_index_v(l)
826             bc_v(l)%end_index(j,i)   = bc_v(l)%start_index(j,i) + num_v_kji(l) - 1
827             start_index_v(l)         = bc_v(l)%end_index(j,i) + 1
[4102]828          ENDDO
[3833]829       ENDDO
[4559]830    ENDDO
[3833]831
832
[4559]833 END SUBROUTINE init_bc
[3833]834
835
[4559]836!--------------------------------------------------------------------------------------------------!
[3833]837! Description:
838! ------------
[4559]839!> Initialize horizontal and vertical surfaces. Counts the number of default-, natural and urban
840!> surfaces and allocates memory, respectively.
841!--------------------------------------------------------------------------------------------------!
842 SUBROUTINE init_surface_arrays
[3833]843
844
[4559]845    USE pegrid
[3833]846
847
[4559]848    IMPLICIT NONE
[3833]849
[4559]850    INTEGER(iwp) ::  i         !< running index x-direction
851    INTEGER(iwp) ::  j         !< running index y-direction
852    INTEGER(iwp) ::  k         !< running index z-direction
853    INTEGER(iwp) ::  l         !< index variable for surface facing
854    INTEGER(iwp) ::  num_lsm_h !< number of horizontally-aligned natural surfaces
855    INTEGER(iwp) ::  num_usm_h !< number of horizontally-aligned urban surfaces
[3833]856
[4559]857    INTEGER(iwp), DIMENSION(0:2) ::  num_def_h !< number of horizontally-aligned default surfaces
858    INTEGER(iwp), DIMENSION(0:3) ::  num_def_v !< number of vertically-aligned default surfaces
859    INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v !< number of vertically-aligned natural surfaces
860    INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v !< number of vertically-aligned urban surfaces
[3833]861
[4559]862    INTEGER(iwp) ::  num_surf_v_l !< number of vertically-aligned local urban/land surfaces
863    INTEGER(iwp) ::  num_surf_v   !< number of vertically-aligned total urban/land surfaces
[3833]864
[4559]865    LOGICAL ::  building             !< flag indicating building grid point
866    LOGICAL ::  terrain              !< flag indicating natural terrain grid point
867    LOGICAL ::  unresolved_building  !< flag indicating a grid point where actually a building is
868                                     !< defined but not resolved by the vertical grid
[3833]869
[4559]870    num_def_h = 0
871    num_def_v = 0
872    num_lsm_h = 0
873    num_lsm_v = 0
874    num_usm_h = 0
875    num_usm_v = 0
[3833]876!
[4559]877!-- Surfaces are classified according to the input data read from static input file. If no input
878!-- file is present, all surfaces are classified either as natural, urban, or default, depending on
879!-- the setting of land_surface and urban_surface. To control this, use the control flag
880!-- topo_no_distinct
[3833]881!
[4559]882!-- Count number of horizontal surfaces on local domain
883    DO  i = nxl, nxr
884       DO  j = nys, nyn
885          DO  k = nzb+1, nzt
[3833]886!
[4559]887!--          Check if current gridpoint belongs to the atmosphere
888             IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
[3833]889!
[4559]890!--             Check if grid point adjoins to any upward-facing horizontal surface, e.g. the Earth
891!--             surface, plane roofs, or ceilings.
[3833]892
[4559]893                IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i), 0 ) )  THEN
[3833]894!
[4559]895!--                Determine flags indicating a terrain surface, a building surface,
896                   terrain  = BTEST( wall_flags_total_0(k-1,j,i), 5 )  .OR.  topo_no_distinct
897                   building = BTEST( wall_flags_total_0(k-1,j,i), 6 )  .OR.  topo_no_distinct
[3833]898!
[4559]899!--                Unresolved_building indicates a surface with equal height as terrain but with a
900!--                non-grid resolved building on top. These surfaces will be flagged as urban
901!--                surfaces.
902                   unresolved_building = BTEST( wall_flags_total_0(k-1,j,i), 5 )  .AND.            &
903                                         BTEST( wall_flags_total_0(k-1,j,i), 6 )
[4159]904!
[4559]905!--                Land-surface type
906                   IF ( land_surface  .AND.  terrain  .AND.  .NOT. unresolved_building )  THEN
907                      num_lsm_h    = num_lsm_h    + 1
[3833]908!
[4559]909!--                Urban surface tpye
910                   ELSEIF ( urban_surface  .AND.  building )  THEN
911                      num_usm_h = num_usm_h + 1
[3833]912!
[4559]913!--                Default-surface type
914                   ELSEIF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
915                      num_def_h(0) = num_def_h(0) + 1
[3833]916!
[4559]917!--                Unclassifified surface-grid point. Give error message.
918                   ELSE
919                      WRITE( message_string, * ) 'Unclassified upward-facing surface element '//   &
920                                                 'at grid point (k,j,i) = ', k, j, i
921                      CALL message( 'surface_mod', 'PA0698', 1, 2, myid, 6, 0 )
922                   ENDIF
[3833]923
[4559]924                ENDIF
[3833]925!
[4559]926!--             Check for top-fluxes
927                IF ( k == nzt  .AND.  use_top_fluxes )  THEN
928                   num_def_h(2) = num_def_h(2) + 1
[3833]929!
[4559]930!--             Check for any other downward-facing surface. So far only for default surface type.
931                ELSEIF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i), 0 ) )  THEN
932                   num_def_h(1) = num_def_h(1) + 1
933                ENDIF
[3833]934
[4559]935             ENDIF
[3833]936          ENDDO
937       ENDDO
[4559]938    ENDDO
[3833]939!
[4559]940!-- Count number of vertical surfaces on local domain
941    DO  i = nxl, nxr
942       DO  j = nys, nyn
943          DO  k = nzb+1, nzt
944             IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
[3833]945!
[4559]946!--             Northward-facing
947                IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i), 0 ) )  THEN
[3833]948!
[4559]949!--                Determine flags indicating terrain or building
[3833]950
[4559]951                   terrain  = BTEST( wall_flags_total_0(k,j-1,i), 5 )  .OR.  topo_no_distinct
952                   building = BTEST( wall_flags_total_0(k,j-1,i), 6 )  .OR.  topo_no_distinct
[4159]953
[4559]954                   unresolved_building = BTEST( wall_flags_total_0(k,j-1,i), 5 )  .AND.            &
955                                         BTEST( wall_flags_total_0(k,j-1,i), 6 )
[4502]956
[4559]957                   IF (  land_surface  .AND.  terrain  .AND.  .NOT. unresolved_building )  THEN
958                      num_lsm_v(0) = num_lsm_v(0) + 1
959                   ELSEIF ( urban_surface  .AND.  building )  THEN
960                      num_usm_v(0) = num_usm_v(0) + 1
[3833]961!
[4559]962!--                Default-surface type
963                   ELSEIF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
964                      num_def_v(0) = num_def_v(0) + 1
[3833]965!
[4559]966!--                Unclassifified surface-grid point. Give error message.
967                   ELSE
968                      WRITE( message_string, * ) 'Unclassified northward-facing surface ' //       &
969                                                 'element at grid point (k,j,i) = ', k, j, i
970                      CALL message( 'surface_mod', 'PA0698', 1, 2, myid, 6, 0 )
[3833]971
972                   ENDIF
[4559]973                ENDIF
[3833]974!
[4559]975!--             Southward-facing
976                IF ( .NOT. BTEST( wall_flags_total_0(k,j+1,i), 0 ) )  THEN
[3833]977!
[4559]978!--                Determine flags indicating terrain or building
979                   terrain  = BTEST( wall_flags_total_0(k,j+1,i), 5 )  .OR.  topo_no_distinct
980                   building = BTEST( wall_flags_total_0(k,j+1,i), 6 )  .OR.  topo_no_distinct
[4502]981
[4559]982                   unresolved_building = BTEST( wall_flags_total_0(k,j+1,i), 5 )  .AND.            &
983                                         BTEST( wall_flags_total_0(k,j+1,i), 6 )
[4502]984
[4559]985                   IF (  land_surface  .AND.  terrain  .AND.  .NOT. unresolved_building )  THEN
986                      num_lsm_v(1) = num_lsm_v(1) + 1
987                   ELSEIF ( urban_surface  .AND.  building )  THEN
988                      num_usm_v(1) = num_usm_v(1) + 1
[3833]989!
[4559]990!--                Default-surface type
991                   ELSEIF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
992                      num_def_v(1) = num_def_v(1) + 1
[3833]993!
[4559]994!--                Unclassifified surface-grid point. Give error message.
995                   ELSE
996                      WRITE( message_string, * ) 'Unclassified southward-facing surface ' //       &
997                                                 'element at grid point (k,j,i) = ', k, j, i
998                      CALL message( 'surface_mod', 'PA0698', 1, 2, myid, 6, 0 )
[3833]999
1000                   ENDIF
[4559]1001                ENDIF
[3833]1002!
[4559]1003!--             Eastward-facing
1004                IF ( .NOT. BTEST( wall_flags_total_0(k,j,i-1), 0 ) )  THEN
[3833]1005!
[4559]1006!--                Determine flags indicating terrain or building
1007                   terrain  = BTEST( wall_flags_total_0(k,j,i-1), 5 )  .OR.  topo_no_distinct
1008                   building = BTEST( wall_flags_total_0(k,j,i-1), 6 )  .OR.  topo_no_distinct
[4502]1009
[4559]1010                   unresolved_building = BTEST( wall_flags_total_0(k,j,i-1), 5 )  .AND.            &
1011                                         BTEST( wall_flags_total_0(k,j,i-1), 6 )
[4502]1012
[4559]1013                   IF ( land_surface  .AND.  terrain  .AND.  .NOT. unresolved_building )  THEN
1014                      num_lsm_v(2) = num_lsm_v(2) + 1
1015                   ELSEIF ( urban_surface  .AND.  building )  THEN
1016                      num_usm_v(2) = num_usm_v(2) + 1
[3833]1017!
[4559]1018!--                Default-surface type
1019                   ELSEIF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
1020                      num_def_v(2) = num_def_v(2) + 1
[3833]1021!
[4559]1022!--                Unclassifified surface-grid point. Give error message.
1023                   ELSE
1024                      WRITE( message_string, * ) 'Unclassified eastward-facing surface ' //        &
1025                                                 'element at grid point (k,j,i) = ', k, j, i
1026                      CALL message( 'surface_mod', 'PA0698', 1, 2, myid, 6, 0 )
[3833]1027
1028                   ENDIF
[4559]1029                ENDIF
[3833]1030!
[4559]1031!--             Westward-facing
1032                IF ( .NOT. BTEST( wall_flags_total_0(k,j,i+1), 0 ) )  THEN
[3833]1033!
[4559]1034!--                Determine flags indicating terrain or building
1035                   terrain  = BTEST( wall_flags_total_0(k,j,i+1), 5 )  .OR.  topo_no_distinct
1036                   building = BTEST( wall_flags_total_0(k,j,i+1), 6 )  .OR.  topo_no_distinct
[4502]1037
[4559]1038                   unresolved_building = BTEST( wall_flags_total_0(k,j,i+1), 5 )  .AND.            &
1039                                         BTEST( wall_flags_total_0(k,j,i+1), 6 )
[4502]1040
[4559]1041                   IF ( land_surface  .AND.  terrain  .AND.  .NOT. unresolved_building )  THEN
1042                      num_lsm_v(3) = num_lsm_v(3) + 1
1043                   ELSEIF ( urban_surface  .AND.  building )  THEN
1044                      num_usm_v(3) = num_usm_v(3) + 1
[3833]1045!
[4559]1046!--                Default-surface type
1047                   ELSEIF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
1048                      num_def_v(3) = num_def_v(3) + 1
[3833]1049!
[4559]1050!--                Unclassifified surface-grid point. Give error message.
1051                   ELSE
1052                      WRITE( message_string, * ) 'Unclassified westward-facing surface ' //        &
1053                                                 'element at grid point (k,j,i) = ', k, j, i
1054                      CALL message( 'surface_mod', 'PA0698', 1, 2, myid, 6, 0 )
[3833]1055
1056                   ENDIF
1057                ENDIF
[4559]1058             ENDIF
[3833]1059          ENDDO
1060       ENDDO
[4559]1061    ENDDO
[3833]1062
1063!
[4559]1064!-- Store number of surfaces per core.
1065!-- Horizontal surface, default type, upward facing
1066    surf_def_h(0)%ns = num_def_h(0)
[3833]1067!
[4559]1068!-- Horizontal surface, default type, downward facing
1069    surf_def_h(1)%ns = num_def_h(1)
[3833]1070!
[4559]1071!-- Horizontal surface, default type, top downward facing
1072    surf_def_h(2)%ns = num_def_h(2)
[3833]1073!
[4559]1074!-- Horizontal surface, natural type, so far only upward-facing
1075    surf_lsm_h%ns    = num_lsm_h
[3833]1076!
[4559]1077!-- Horizontal surface, urban type, so far only upward-facing
1078    surf_usm_h%ns    = num_usm_h
[3833]1079!
[4559]1080!-- Vertical surface, default type, northward facing
1081    surf_def_v(0)%ns = num_def_v(0)
[3833]1082!
[4559]1083!-- Vertical surface, default type, southward facing
1084    surf_def_v(1)%ns = num_def_v(1)
[3833]1085!
[4559]1086!-- Vertical surface, default type, eastward facing
1087    surf_def_v(2)%ns = num_def_v(2)
[3833]1088!
[4559]1089!-- Vertical surface, default type, westward facing
1090    surf_def_v(3)%ns = num_def_v(3)
[3833]1091!
[4559]1092!-- Vertical surface, natural type, northward facing
1093    surf_lsm_v(0)%ns = num_lsm_v(0)
[3833]1094!
[4559]1095!-- Vertical surface, natural type, southward facing
1096    surf_lsm_v(1)%ns = num_lsm_v(1)
[3833]1097!
[4559]1098!-- Vertical surface, natural type, eastward facing
1099    surf_lsm_v(2)%ns = num_lsm_v(2)
[3833]1100!
[4559]1101!-- Vertical surface, natural type, westward facing
1102    surf_lsm_v(3)%ns = num_lsm_v(3)
[3833]1103!
[4559]1104!-- Vertical surface, urban type, northward facing
1105    surf_usm_v(0)%ns = num_usm_v(0)
[3833]1106!
[4559]1107!-- Vertical surface, urban type, southward facing
1108    surf_usm_v(1)%ns = num_usm_v(1)
[3833]1109!
[4559]1110!-- Vertical surface, urban type, eastward facing
1111    surf_usm_v(2)%ns = num_usm_v(2)
[3833]1112!
[4559]1113!-- Vertical surface, urban type, westward facing
1114    surf_usm_v(3)%ns = num_usm_v(3)
[3833]1115!
[4559]1116!-- Allocate required attributes for horizontal surfaces - default type.
1117!-- Upward-facing (l=0) and downward-facing (l=1).
1118    DO  l = 0, 1
1119       CALL allocate_surface_attributes_h ( surf_def_h(l), nys, nyn, nxl, nxr )
1120    ENDDO
[3833]1121!
[4559]1122!-- Allocate required attributes for model top
1123    CALL allocate_surface_attributes_h_top ( surf_def_h(2), nys, nyn, nxl, nxr )
[3833]1124!
[4559]1125!-- Allocate required attributes for horizontal surfaces - natural type.
1126    CALL allocate_surface_attributes_h ( surf_lsm_h, nys, nyn, nxl, nxr )
[3833]1127!
[4559]1128!-- Allocate required attributes for horizontal surfaces - urban type.
1129    CALL allocate_surface_attributes_h ( surf_usm_h, nys, nyn, nxl, nxr )
[3833]1130
1131!
[4559]1132!-- Allocate required attributes for vertical surfaces.
1133!-- Northward-facing (l=0), southward-facing (l=1), eastward-facing (l=2) and westward-facing (l=3).
1134!-- Default type.
1135    DO  l = 0, 3
1136       CALL allocate_surface_attributes_v ( surf_def_v(l), nys, nyn, nxl, nxr )
1137    ENDDO
[3833]1138!
[4559]1139!-- Natural type
1140    DO  l = 0, 3
1141       CALL allocate_surface_attributes_v ( surf_lsm_v(l), nys, nyn, nxl, nxr )
1142    ENDDO
[3833]1143!
[4559]1144!-- Urban type
1145    DO  l = 0, 3
1146       CALL allocate_surface_attributes_v ( surf_usm_v(l), nys, nyn, nxl, nxr )
1147    ENDDO
[3833]1148!
[4559]1149!-- Set the flag for the existence of vertical urban/land surfaces
1150    num_surf_v_l = 0
1151    DO  l = 0, 3
1152       num_surf_v_l = num_surf_v_l + surf_usm_v(l)%ns + surf_lsm_v(l)%ns
1153    ENDDO
[3833]1154
1155#if defined( __parallel )
[4559]1156    CALL MPI_ALLREDUCE( num_surf_v_l, num_surf_v, 1, MPI_INTEGER, MPI_SUM, comm2d, ierr)
[3833]1157#else
[4559]1158    num_surf_v = num_surf_v_l
[3833]1159#endif
[4559]1160    IF ( num_surf_v > 0 )  vertical_surfaces_exist = .TRUE.
[3833]1161
[4502]1162
[4559]1163 END SUBROUTINE init_surface_arrays
[3833]1164
1165
[4559]1166!--------------------------------------------------------------------------------------------------!
[3833]1167! Description:
1168! ------------
1169!> Enter horizontal and vertical surfaces.
[4559]1170!--------------------------------------------------------------------------------------------------!
[3833]1171#if defined( _OPENACC )
[4559]1172 SUBROUTINE enter_surface_arrays
[3833]1173
[4559]1174    IMPLICIT NONE
[3833]1175
[4559]1176    INTEGER(iwp) ::  l  !<
[4502]1177
[4559]1178    !$ACC ENTER DATA &
1179    !$ACC COPYIN(surf_def_h(0:2)) &
1180    !$ACC COPYIN(surf_def_v(0:3)) &
1181    !$ACC COPYIN(surf_lsm_h) &
1182    !$ACC COPYIN(surf_lsm_v(0:3)) &
1183    !$ACC COPYIN(surf_usm_h) &
1184    !$ACC COPYIN(surf_usm_v(0:3))
1185!
1186!-- Copy data in surf_def_h(0:2)
1187    DO  l = 0, 1
1188       CALL enter_surface_attributes_h( surf_def_h(l) )
1189    ENDDO
1190    CALL enter_surface_attributes_h_top( surf_def_h(2) )
1191!
1192!-- Copy data in surf_def_v(0:3)
1193    DO  l = 0, 3
1194       CALL enter_surface_attributes_v( surf_def_v(l) )
1195    ENDDO
1196!
1197!-- Copy data in surf_lsm_h
1198    CALL enter_surface_attributes_h( surf_lsm_h )
1199!
1200!-- Copy data in surf_lsm_v(0:3)
1201    DO  l = 0, 3
1202       CALL enter_surface_attributes_v( surf_lsm_v(l) )
1203    ENDDO
1204!
1205!-- Copy data in surf_usm_h
1206    CALL enter_surface_attributes_h( surf_usm_h )
1207!
1208!-- Copy data in surf_usm_v(0:3)
1209    DO  l = 0, 3
1210       CALL enter_surface_attributes_v( surf_usm_v(l) )
1211    ENDDO
[3833]1212
[4559]1213 END SUBROUTINE enter_surface_arrays
[3833]1214#endif
1215
[4559]1216!--------------------------------------------------------------------------------------------------!
[3833]1217! Description:
1218! ------------
1219!> Exit horizontal and vertical surfaces.
[4559]1220!--------------------------------------------------------------------------------------------------!
[3833]1221#if defined( _OPENACC )
[4559]1222 SUBROUTINE exit_surface_arrays
[3833]1223
[4559]1224    IMPLICIT NONE
[3833]1225
[4559]1226    INTEGER(iwp) ::  l  !<
1227!
1228!-- Delete data in surf_def_h(0:2)
1229    DO  l = 0, 1
1230       CALL exit_surface_attributes_h( surf_def_h(l) )
1231    ENDDO
1232    CALL exit_surface_attributes_h( surf_def_h(2) )
1233!
1234!-- Delete data in surf_def_v(0:3)
1235    DO  l = 0, 3
1236       CALL exit_surface_attributes_v( surf_def_v(l) )
1237    ENDDO
1238!
1239!-- Delete data in surf_lsm_h
1240    CALL exit_surface_attributes_h( surf_lsm_h )
1241!
1242!-- Delete data in surf_lsm_v(0:3)
1243    DO  l = 0, 3
1244       CALL exit_surface_attributes_v( surf_lsm_v(l) )
1245    ENDDO
1246!
1247!-- Delete data in surf_usm_h
1248    CALL exit_surface_attributes_h( surf_usm_h )
1249!
1250!-- Delete data in surf_usm_v(0:3)
1251    DO  l = 0, 3
1252       CALL exit_surface_attributes_v( surf_usm_v(l) )
1253    ENDDO
[4502]1254
[4559]1255    !$ACC EXIT DATA &
1256    !$ACC DELETE(surf_def_h(0:2)) &
1257    !$ACC DELETE(surf_def_v(0:3)) &
1258    !$ACC DELETE(surf_lsm_h) &
1259    !$ACC DELETE(surf_lsm_v(0:3)) &
1260    !$ACC DELETE(surf_usm_h) &
1261    !$ACC DELETE(surf_usm_v(0:3))
[3833]1262
[4559]1263 END SUBROUTINE exit_surface_arrays
[3833]1264#endif
1265
[4559]1266!--------------------------------------------------------------------------------------------------!
[3833]1267! Description:
1268! ------------
[4559]1269!> Deallocating memory for upward and downward-facing horizontal surface types, except for top
1270!> fluxes.
1271!--------------------------------------------------------------------------------------------------!
1272 SUBROUTINE deallocate_surface_attributes_h( surfaces )
[3833]1273
[4559]1274    IMPLICIT NONE
[3833]1275
1276
[4559]1277    TYPE(surf_type) ::  surfaces  !< respective surface type
[3833]1278
1279
[4559]1280    DEALLOCATE ( surfaces%start_index )
1281    DEALLOCATE ( surfaces%end_index )
[3833]1282!
[4559]1283!-- Indices to locate surface element
1284    DEALLOCATE ( surfaces%i )
1285    DEALLOCATE ( surfaces%j )
1286    DEALLOCATE ( surfaces%k )
[3833]1287!
[4559]1288!-- Surface-layer height
1289    DEALLOCATE ( surfaces%z_mo )
[3833]1290!
[4559]1291!-- Surface orientation
1292    DEALLOCATE ( surfaces%facing )
[3833]1293!
[4559]1294!-- Surface-parallel wind velocity
1295    DEALLOCATE ( surfaces%uvw_abs )
[3833]1296!
[4593]1297!-- Pre-calculated ln(z/z0)
1298    DEALLOCATE ( surfaces%ln_z_z0  )
1299    DEALLOCATE ( surfaces%ln_z_z0h )
1300    DEALLOCATE ( surfaces%ln_z_z0q )
1301!
[4559]1302!-- Roughness
1303    DEALLOCATE ( surfaces%z0 )
1304    DEALLOCATE ( surfaces%z0h )
1305    DEALLOCATE ( surfaces%z0q )
[3833]1306!
[4559]1307!-- Friction velocity
1308    DEALLOCATE ( surfaces%us )
[3833]1309!
[4559]1310!-- Stability parameter
1311    DEALLOCATE ( surfaces%ol )
[3833]1312!
[4559]1313!-- Bulk Richardson number
1314    DEALLOCATE ( surfaces%rib )
[3833]1315!
[4559]1316!-- Vertical momentum fluxes of u and v
1317    DEALLOCATE ( surfaces%usws )
1318    DEALLOCATE ( surfaces%vsws )
[3833]1319!
[4559]1320!-- Required in production_e
1321    IF ( .NOT. constant_diffusion )  THEN
1322       DEALLOCATE ( surfaces%u_0 )
1323       DEALLOCATE ( surfaces%v_0 )
1324    ENDIF
[3833]1325!
[4559]1326!-- Characteristic temperature and surface flux of sensible heat
1327    DEALLOCATE ( surfaces%ts )
1328    DEALLOCATE ( surfaces%shf )
[3833]1329!
[4559]1330!-- Surface temperature
1331    DEALLOCATE ( surfaces%pt_surface )
[3833]1332!
[4559]1333!-- Characteristic humidity and surface flux of latent heat
1334    IF ( humidity )  THEN
1335       DEALLOCATE ( surfaces%qs )
1336       DEALLOCATE ( surfaces%qsws )
1337       DEALLOCATE ( surfaces%q_surface   )
1338       DEALLOCATE ( surfaces%vpt_surface )
1339    ENDIF
[3833]1340!
[4559]1341!-- Characteristic scalar and surface flux of scalar
1342    IF ( passive_scalar )  THEN
1343       DEALLOCATE ( surfaces%ss )
1344       DEALLOCATE ( surfaces%ssws )
1345    ENDIF
[3833]1346!
[4559]1347!-- Scaling parameter (cs*) and surface flux of chemical species
1348    IF ( air_chemistry )  THEN
1349       DEALLOCATE ( surfaces%css )
1350       DEALLOCATE ( surfaces%cssws )
1351    ENDIF
[3833]1352!
[4559]1353!-- Arrays for storing potential temperature and mixing ratio at first grid level
1354    DEALLOCATE ( surfaces%pt1 )
1355    DEALLOCATE ( surfaces%qv1 )
1356    DEALLOCATE ( surfaces%vpt1 )
[4502]1357
[3833]1358!
[4502]1359!--
[4559]1360    IF ( surf_bulk_cloud_model .AND. surf_microphysics_morrison)  THEN
1361       DEALLOCATE ( surfaces%qcs )
1362       DEALLOCATE ( surfaces%ncs )
1363       DEALLOCATE ( surfaces%qcsws )
1364       DEALLOCATE ( surfaces%ncsws )
1365    ENDIF
[3833]1366!
[4502]1367!--
[4559]1368    IF ( surf_bulk_cloud_model .AND. surf_microphysics_seifert)  THEN
1369       DEALLOCATE ( surfaces%qrs )
1370       DEALLOCATE ( surfaces%nrs )
1371       DEALLOCATE ( surfaces%qrsws )
1372       DEALLOCATE ( surfaces%nrsws )
1373    ENDIF
[3833]1374!
[4502]1375!--
[4559]1376    IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_phase)  THEN
1377       DEALLOCATE ( surfaces%qis )
1378       DEALLOCATE ( surfaces%nis )
1379       DEALLOCATE ( surfaces%qisws )
1380       DEALLOCATE ( surfaces%nisws )
1381    ENDIF
[4502]1382!
[4559]1383!-- Salinity surface flux
1384    IF ( ocean_mode )  DEALLOCATE ( surfaces%sasws )
[3833]1385
[4559]1386 END SUBROUTINE deallocate_surface_attributes_h
[3833]1387
1388
[4559]1389!--------------------------------------------------------------------------------------------------!
[3833]1390! Description:
1391! ------------
[4559]1392!> Allocating memory for upward and downward-facing horizontal surface types, except for top fluxes.
1393!--------------------------------------------------------------------------------------------------!
1394 SUBROUTINE allocate_surface_attributes_h( surfaces, nys_l, nyn_l, nxl_l, nxr_l )
[3833]1395
[4559]1396    IMPLICIT NONE
[3833]1397
[4559]1398    INTEGER(iwp) ::  nyn_l  !< north bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1399    INTEGER(iwp) ::  nys_l  !< south bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1400    INTEGER(iwp) ::  nxl_l  !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1401    INTEGER(iwp) ::  nxr_l  !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
[3833]1402
[4559]1403    TYPE(surf_type) ::  surfaces  !< respective surface type
[3833]1404
1405!
[4559]1406!-- Allocate arrays for start and end index of horizontal surface type for each (j,i)-grid point.
1407!-- This is required e.g. in diffion_x, which is called for each (j,i). In order to find the
1408!-- location where the respective flux is store within the surface-type, start- and end-index are
1409!-- stored for each (j,i). For example, each (j,i) can have several entries where fluxes for
1410!-- horizontal surfaces might be stored, e.g. for overhanging structures where several upward-facing
1411!-- surfaces might exist for given (j,i). If no surface of respective type exist at current (j,i),
1412!-- set indicies such that loop in diffusion routines will not be entered.
1413    ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
1414    ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l) )
1415    surfaces%start_index = 0
1416    surfaces%end_index   = -1
[3833]1417!
[4559]1418!-- Indices to locate surface element
1419    ALLOCATE ( surfaces%i(1:surfaces%ns) )
1420    ALLOCATE ( surfaces%j(1:surfaces%ns) )
1421    ALLOCATE ( surfaces%k(1:surfaces%ns) )
[3833]1422!
[4559]1423!-- Surface-layer height
1424    ALLOCATE ( surfaces%z_mo(1:surfaces%ns) )
[3833]1425!
[4559]1426!-- Surface orientation
1427    ALLOCATE ( surfaces%facing(1:surfaces%ns) )
[3833]1428!
[4559]1429!-- Surface-parallel wind velocity
1430    ALLOCATE ( surfaces%uvw_abs(1:surfaces%ns) )
[3833]1431!
[4593]1432!-- Precalculated ln(z/z0)
1433    ALLOCATE( surfaces%ln_z_z0(1:surfaces%ns)  )
1434    ALLOCATE( surfaces%ln_z_z0h(1:surfaces%ns) )
1435    ALLOCATE( surfaces%ln_z_z0q(1:surfaces%ns) )
1436!
[4559]1437!-- Roughness
1438    ALLOCATE ( surfaces%z0(1:surfaces%ns)  )
1439    ALLOCATE ( surfaces%z0h(1:surfaces%ns) )
1440    ALLOCATE ( surfaces%z0q(1:surfaces%ns) )
[3833]1441!
[4559]1442!-- Friction velocity
1443    ALLOCATE ( surfaces%us(1:surfaces%ns) )
[3833]1444!
[4559]1445!-- Stability parameter
1446    ALLOCATE ( surfaces%ol(1:surfaces%ns) )
[3833]1447!
[4559]1448!-- Bulk Richardson number
1449    ALLOCATE ( surfaces%rib(1:surfaces%ns) )
[3833]1450!
[4559]1451!-- Vertical momentum fluxes of u and v
1452    ALLOCATE ( surfaces%usws(1:surfaces%ns) )
1453    ALLOCATE ( surfaces%vsws(1:surfaces%ns) )
[3833]1454!
[4559]1455!-- Required in production_e
1456    IF ( .NOT. constant_diffusion )  THEN
1457       ALLOCATE ( surfaces%u_0(1:surfaces%ns) )
1458       ALLOCATE ( surfaces%v_0(1:surfaces%ns) )
1459    ENDIF
[3833]1460!
[4559]1461!-- Characteristic temperature and surface flux of sensible heat
1462    ALLOCATE ( surfaces%ts(1:surfaces%ns) )
1463    ALLOCATE ( surfaces%shf(1:surfaces%ns) )
[3833]1464!
[4559]1465!-- Surface temperature
1466    ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) )
[3833]1467!
[4559]1468!-- Characteristic humidity, surface flux of latent heat, and surface virtual potential temperature
1469    IF ( humidity )  THEN
1470       ALLOCATE ( surfaces%qs(1:surfaces%ns) )
1471       ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
1472       ALLOCATE ( surfaces%q_surface(1:surfaces%ns) )
1473       ALLOCATE ( surfaces%vpt_surface(1:surfaces%ns) )
1474    ENDIF
[3833]1475
1476!
[4559]1477!-- Characteristic scalar and surface flux of scalar
1478    IF ( passive_scalar )  THEN
1479       ALLOCATE ( surfaces%ss(1:surfaces%ns) )
1480       ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
1481    ENDIF
[3833]1482!
[4559]1483!-- Scaling parameter (cs*) and surface flux of chemical species
1484    IF ( air_chemistry )  THEN
1485       ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns) )
1486       ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
1487    ENDIF
[3833]1488!
[4559]1489!-- Arrays for storing potential temperature and mixing ratio at first grid level
1490    ALLOCATE ( surfaces%pt1(1:surfaces%ns) )
1491    ALLOCATE ( surfaces%qv1(1:surfaces%ns) )
1492    ALLOCATE ( surfaces%vpt1(1:surfaces%ns) )
[3833]1493!
[4502]1494!--
[4559]1495    IF ( surf_bulk_cloud_model .AND. surf_microphysics_morrison)  THEN
1496       ALLOCATE ( surfaces%qcs(1:surfaces%ns) )
1497       ALLOCATE ( surfaces%ncs(1:surfaces%ns) )
1498       ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
1499       ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
1500    ENDIF
[3833]1501!
[4502]1502!--
[4559]1503    IF ( surf_bulk_cloud_model .AND. surf_microphysics_seifert)  THEN
1504       ALLOCATE ( surfaces%qrs(1:surfaces%ns) )
1505       ALLOCATE ( surfaces%nrs(1:surfaces%ns) )
1506       ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
1507       ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
1508    ENDIF
[4502]1509
[3833]1510!
[4502]1511!--
[4559]1512    IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_phase)  THEN
1513       ALLOCATE ( surfaces%qis(1:surfaces%ns) )
1514       ALLOCATE ( surfaces%nis(1:surfaces%ns) )
1515       ALLOCATE ( surfaces%qisws(1:surfaces%ns) )
1516       ALLOCATE ( surfaces%nisws(1:surfaces%ns) )
1517    ENDIF
[4502]1518
1519!
[4559]1520!-- Salinity surface flux
1521    IF ( ocean_mode )  ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
[3833]1522
[4559]1523 END SUBROUTINE allocate_surface_attributes_h
[3833]1524
1525
[4559]1526!--------------------------------------------------------------------------------------------------!
[3833]1527! Description:
1528! ------------
[4559]1529!> Exit memory for upward and downward-facing horizontal surface types, except for top fluxes.
1530!--------------------------------------------------------------------------------------------------!
[3833]1531#if defined( _OPENACC )
[4559]1532 SUBROUTINE exit_surface_attributes_h( surfaces )
[3833]1533
[4559]1534    IMPLICIT NONE
[4502]1535
[4559]1536    TYPE(surf_type) ::  surfaces  !< respective surface type
[4502]1537
[4559]1538    !$ACC EXIT DATA &
1539    !$ACC DELETE(surfaces%start_index(nys:nyn,nxl:nxr)) &
1540    !$ACC DELETE(surfaces%end_index(nys:nyn,nxl:nxr)) &
1541    !$ACC DELETE(surfaces%i(1:surfaces%ns)) &
1542    !$ACC DELETE(surfaces%j(1:surfaces%ns)) &
1543    !$ACC DELETE(surfaces%k(1:surfaces%ns)) &
1544    !$ACC DELETE(surfaces%z_mo(1:surfaces%ns)) &
1545    !$ACC DELETE(surfaces%uvw_abs(1:surfaces%ns)) &
[4593]1546    !$ACC DELETE(surfaces%ln_z_z0(1:surfaces%ns)) &
1547    !$ACC DELETE(surfaces%ln_z_z0h(1:surfaces%ns)) &
1548    !$ACC DELETE(surfaces%ln_z_z0q(1:surfaces%ns)) &
[4559]1549    !$ACC DELETE(surfaces%z0(1:surfaces%ns)) &
[4594]1550    !$ACC DELETE(surfaces%z0h(1:surfaces%ns)) &
1551    !$ACC DELETE(surfaces%z0q(1:surfaces%ns)) &
[4559]1552    !$ACC COPYOUT(surfaces%us(1:surfaces%ns)) &
1553    !$ACC COPYOUT(surfaces%ol(1:surfaces%ns)) &
1554    !$ACC DELETE(surfaces%rib(1:surfaces%ns)) &
1555    !$ACC COPYOUT(surfaces%usws(1:surfaces%ns)) &
1556    !$ACC COPYOUT(surfaces%vsws(1:surfaces%ns)) &
1557    !$ACC COPYOUT(surfaces%ts(1:surfaces%ns)) &
1558    !$ACC COPYOUT(surfaces%shf(1:surfaces%ns)) &
1559    !$ACC DELETE(surfaces%pt_surface(1:surfaces%ns)) &
1560    !$ACC DELETE(surfaces%pt1(1:surfaces%ns)) &
1561    !$ACC DELETE(surfaces%qv1(1:surfaces%ns))
1562
1563    IF ( .NOT. constant_diffusion )  THEN
[3833]1564       !$ACC EXIT DATA &
[4559]1565       !$ACC DELETE(surfaces%u_0(1:surfaces%ns)) &
1566       !$ACC DELETE(surfaces%v_0(1:surfaces%ns))
1567    ENDIF
[3833]1568
[4559]1569 END SUBROUTINE exit_surface_attributes_h
[3833]1570#endif
1571
[4559]1572!--------------------------------------------------------------------------------------------------!
[3833]1573! Description:
1574! ------------
[4559]1575!> Enter memory for upward and downward-facing horizontal surface types, except for top fluxes.
1576!--------------------------------------------------------------------------------------------------!
[3833]1577#if defined( _OPENACC )
[4559]1578 SUBROUTINE enter_surface_attributes_h( surfaces )
[3833]1579
[4559]1580    IMPLICIT NONE
[3833]1581
[4559]1582    TYPE(surf_type) ::  surfaces  !< respective surface type
[3833]1583
[4559]1584    !$ACC ENTER DATA &
1585    !$ACC COPYIN(surfaces%start_index(nys:nyn,nxl:nxr)) &
1586    !$ACC COPYIN(surfaces%end_index(nys:nyn,nxl:nxr)) &
1587    !$ACC COPYIN(surfaces%i(1:surfaces%ns)) &
1588    !$ACC COPYIN(surfaces%j(1:surfaces%ns)) &
1589    !$ACC COPYIN(surfaces%k(1:surfaces%ns)) &
1590    !$ACC COPYIN(surfaces%z_mo(1:surfaces%ns)) &
1591    !$ACC COPYIN(surfaces%uvw_abs(1:surfaces%ns)) &
[4593]1592    !$ACC COPYIN(surfaces%ln_z_z0(1:surfaces%ns)) &
1593    !$ACC COPYIN(surfaces%ln_z_z0h(1:surfaces%ns)) &
1594    !$ACC COPYIN(surfaces%ln_z_z0q(1:surfaces%ns)) &
[4559]1595    !$ACC COPYIN(surfaces%z0(1:surfaces%ns)) &
[4594]1596    !$ACC COPYIN(surfaces%z0h(1:surfaces%ns)) &
1597    !$ACC COPYIN(surfaces%z0q(1:surfaces%ns)) &
[4559]1598    !$ACC COPYIN(surfaces%us(1:surfaces%ns)) &
1599    !$ACC COPYIN(surfaces%ol(1:surfaces%ns)) &
1600    !$ACC COPYIN(surfaces%rib(1:surfaces%ns)) &
1601    !$ACC COPYIN(surfaces%usws(1:surfaces%ns)) &
1602    !$ACC COPYIN(surfaces%vsws(1:surfaces%ns)) &
1603    !$ACC COPYIN(surfaces%ts(1:surfaces%ns)) &
1604    !$ACC COPYIN(surfaces%shf(1:surfaces%ns)) &
1605    !$ACC COPYIN(surfaces%pt1(1:surfaces%ns)) &
1606    !$ACC COPYIN(surfaces%qv1(1:surfaces%ns)) &
1607    !$ACC COPYIN(surfaces%pt_surface(1:surfaces%ns))
1608
1609    IF ( .NOT. constant_diffusion )  THEN
[3833]1610       !$ACC ENTER DATA &
[4559]1611       !$ACC COPYIN(surfaces%u_0(1:surfaces%ns)) &
1612       !$ACC COPYIN(surfaces%v_0(1:surfaces%ns))
1613    ENDIF
[3833]1614
[4559]1615 END SUBROUTINE enter_surface_attributes_h
[3833]1616#endif
1617
[4559]1618!--------------------------------------------------------------------------------------------------!
[3833]1619! Description:
1620! ------------
[4502]1621!> Deallocating memory for model-top fluxes
[4559]1622!--------------------------------------------------------------------------------------------------!
1623 SUBROUTINE deallocate_surface_attributes_h_top( surfaces )
[3833]1624
[4559]1625    IMPLICIT NONE
[3833]1626
1627
[4559]1628    TYPE(surf_type) ::  surfaces !< respective surface type
[3833]1629
[4559]1630    DEALLOCATE ( surfaces%start_index )
1631    DEALLOCATE ( surfaces%end_index )
[3833]1632!
[4559]1633!-- Indices to locate surface (model-top) element
1634    DEALLOCATE ( surfaces%i )
1635    DEALLOCATE ( surfaces%j )
1636    DEALLOCATE ( surfaces%k )
[3833]1637
[4559]1638    IF ( .NOT. constant_diffusion )  THEN
1639       DEALLOCATE ( surfaces%u_0 )
1640       DEALLOCATE ( surfaces%v_0 )
1641    ENDIF
[3833]1642!
[4559]1643!-- Vertical momentum fluxes of u and v
1644    DEALLOCATE ( surfaces%usws )
1645    DEALLOCATE ( surfaces%vsws )
[3833]1646!
[4559]1647!-- Sensible heat flux
1648    DEALLOCATE ( surfaces%shf )
[3833]1649!
[4559]1650!-- Latent heat flux
1651    IF ( humidity .OR. coupling_mode == 'ocean_to_atmosphere')  THEN
1652       DEALLOCATE ( surfaces%qsws )
1653    ENDIF
[3833]1654!
[4559]1655!-- Scalar flux
1656    IF ( passive_scalar )  THEN
1657       DEALLOCATE ( surfaces%ssws )
1658    ENDIF
[3833]1659!
[4559]1660!-- Chemical species flux
1661    IF ( air_chemistry )  THEN
1662       DEALLOCATE ( surfaces%cssws )
1663    ENDIF
[3833]1664!
[4502]1665!--
[4559]1666    IF ( surf_bulk_cloud_model .AND. surf_microphysics_morrison)  THEN
1667       DEALLOCATE ( surfaces%qcsws )
1668       DEALLOCATE ( surfaces%ncsws )
1669    ENDIF
[3833]1670!
[4502]1671!--
[4559]1672    IF ( surf_bulk_cloud_model .AND. surf_microphysics_seifert)  THEN
1673       DEALLOCATE ( surfaces%qrsws )
1674       DEALLOCATE ( surfaces%nrsws )
1675    ENDIF
[4502]1676
[3833]1677!
[4502]1678!--
[4559]1679    IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_phase)  THEN
1680       DEALLOCATE ( surfaces%qisws )
1681       DEALLOCATE ( surfaces%nisws )
1682    ENDIF
[4502]1683!
[4559]1684!-- Salinity flux
1685    IF ( ocean_mode )  DEALLOCATE ( surfaces%sasws )
[3833]1686
[4559]1687 END SUBROUTINE deallocate_surface_attributes_h_top
[3833]1688
1689
[4559]1690!--------------------------------------------------------------------------------------------------!
[3833]1691! Description:
1692! ------------
[4502]1693!> Allocating memory for model-top fluxes
[4559]1694!--------------------------------------------------------------------------------------------------!
1695 SUBROUTINE allocate_surface_attributes_h_top( surfaces, nys_l, nyn_l, nxl_l, nxr_l )
[3833]1696
[4559]1697    IMPLICIT NONE
[3833]1698
[4559]1699    INTEGER(iwp) ::  nyn_l  !< north bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1700    INTEGER(iwp) ::  nys_l  !< south bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1701    INTEGER(iwp) ::  nxl_l  !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1702    INTEGER(iwp) ::  nxr_l  !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
[3833]1703
[4559]1704    TYPE(surf_type) ::  surfaces !< respective surface type
[3833]1705
[4559]1706    ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
1707    ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l) )
1708    surfaces%start_index = 0
1709    surfaces%end_index   = -1
[3833]1710!
[4559]1711!-- Indices to locate surface (model-top) element
1712    ALLOCATE ( surfaces%i(1:surfaces%ns) )
1713    ALLOCATE ( surfaces%j(1:surfaces%ns) )
1714    ALLOCATE ( surfaces%k(1:surfaces%ns) )
[3833]1715
[4559]1716    IF ( .NOT. constant_diffusion )  THEN
1717       ALLOCATE ( surfaces%u_0(1:surfaces%ns) )
1718       ALLOCATE ( surfaces%v_0(1:surfaces%ns) )
1719    ENDIF
[3833]1720!
[4559]1721!-- Vertical momentum fluxes of u and v
1722    ALLOCATE ( surfaces%usws(1:surfaces%ns) )
1723    ALLOCATE ( surfaces%vsws(1:surfaces%ns) )
[3833]1724!
[4559]1725!-- Sensible heat flux
1726    ALLOCATE ( surfaces%shf(1:surfaces%ns) )
[3833]1727!
[4559]1728!-- Latent heat flux
1729    IF ( humidity .OR. coupling_mode == 'ocean_to_atmosphere')  THEN
1730       ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
1731    ENDIF
[3833]1732!
[4559]1733!-- Scalar flux
1734    IF ( passive_scalar )  THEN
1735       ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
1736    ENDIF
[3833]1737!
[4559]1738!-- Chemical species flux
1739    IF ( air_chemistry )  THEN
1740       ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
1741    ENDIF
[3833]1742!
[4502]1743!--
[4559]1744    IF ( surf_bulk_cloud_model .AND. surf_microphysics_morrison )  THEN
1745       ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
1746       ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
1747    ENDIF
[3833]1748!
[4502]1749!--
[4559]1750    IF ( surf_bulk_cloud_model .AND. surf_microphysics_seifert )  THEN
1751       ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
1752       ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
1753    ENDIF
[4502]1754
[3833]1755!
[4502]1756!--
[4559]1757    IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_phase )  THEN
1758       ALLOCATE ( surfaces%qisws(1:surfaces%ns) )
1759       ALLOCATE ( surfaces%nisws(1:surfaces%ns) )
1760    ENDIF
[4502]1761!
[4559]1762!-- Salinity flux
1763    IF ( ocean_mode )  ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
[3833]1764
[4559]1765 END SUBROUTINE allocate_surface_attributes_h_top
[3833]1766
1767
[4559]1768!--------------------------------------------------------------------------------------------------!
[3833]1769! Description:
1770! ------------
1771!> Exit memory for model-top fluxes.
[4559]1772!--------------------------------------------------------------------------------------------------!
[3833]1773#if defined( _OPENACC )
[4559]1774 SUBROUTINE exit_surface_attributes_h_top( surfaces )
[3833]1775
[4559]1776    IMPLICIT NONE
[4502]1777
[4559]1778    TYPE(surf_type) ::  surfaces  !< respective surface type
[4502]1779
[4559]1780    !$ACC EXIT DATA &
1781    !$ACC DELETE(surfaces%start_index(nys:nyn,nxl:nxr)) &
1782    !$ACC DELETE(surfaces%end_index(nys:nyn,nxl:nxr)) &
1783    !$ACC DELETE(surfaces%i(1:surfaces%ns)) &
1784    !$ACC DELETE(surfaces%j(1:surfaces%ns)) &
1785    !$ACC DELETE(surfaces%k(1:surfaces%ns)) &
1786    !$ACC DELETE(surfaces%usws(1:surfaces%ns)) &
1787    !$ACC DELETE(surfaces%vsws(1:surfaces%ns)) &
1788    !$ACC DELETE(surfaces%shf(1:surfaces%ns))
1789
1790    IF ( .NOT. constant_diffusion )  THEN
[3833]1791       !$ACC EXIT DATA &
[4559]1792       !$ACC DELETE(surfaces%u_0(1:surfaces%ns)) &
1793       !$ACC DELETE(surfaces%v_0(1:surfaces%ns))
1794    ENDIF
[3833]1795
[4559]1796 END SUBROUTINE exit_surface_attributes_h_top
[3833]1797#endif
1798
[4559]1799!--------------------------------------------------------------------------------------------------!
[3833]1800! Description:
1801! ------------
1802!> Enter memory for model-top fluxes.
[4559]1803!--------------------------------------------------------------------------------------------------!
[3833]1804#if defined( _OPENACC )
[4559]1805 SUBROUTINE enter_surface_attributes_h_top( surfaces )
[3833]1806
[4559]1807    IMPLICIT NONE
[3833]1808
[4559]1809    TYPE(surf_type) ::  surfaces  !< respective surface type
[3833]1810
[4559]1811    !$ACC ENTER DATA &
1812    !$ACC COPYIN(surfaces%start_index(nys:nyn,nxl:nxr)) &
1813    !$ACC COPYIN(surfaces%end_index(nys:nyn,nxl:nxr)) &
1814    !$ACC COPYIN(surfaces%i(1:surfaces%ns)) &
1815    !$ACC COPYIN(surfaces%j(1:surfaces%ns)) &
1816    !$ACC COPYIN(surfaces%k(1:surfaces%ns)) &
1817    !$ACC COPYIN(surfaces%usws(1:surfaces%ns)) &
1818    !$ACC COPYIN(surfaces%vsws(1:surfaces%ns)) &
1819    !$ACC COPYIN(surfaces%shf(1:surfaces%ns))
1820
1821    IF ( .NOT. constant_diffusion )  THEN
[3833]1822       !$ACC ENTER DATA &
[4559]1823       !$ACC COPYIN(surfaces%u_0(1:surfaces%ns)) &
1824       !$ACC COPYIN(surfaces%v_0(1:surfaces%ns))
1825    ENDIF
[3833]1826
[4559]1827 END SUBROUTINE enter_surface_attributes_h_top
[3833]1828#endif
1829
[4559]1830!--------------------------------------------------------------------------------------------------!
[3833]1831! Description:
1832! ------------
[4502]1833!> Deallocating memory for vertical surface types.
[4559]1834!--------------------------------------------------------------------------------------------------!
1835 SUBROUTINE deallocate_surface_attributes_v( surfaces )
[3833]1836
[4559]1837    IMPLICIT NONE
[3833]1838
1839
[4559]1840    TYPE(surf_type) ::  surfaces !< respective surface type
[3833]1841
1842!
[4559]1843!-- Allocate arrays for start and end index of vertical surface type for each (j,i)-grid point. This
1844!-- is required in diffion_x, which is called for each (j,i). In order to find the location where
1845!-- the respective flux is store within the surface-type, start- and end-index are stored for each
1846!-- (j,i). For example, each (j,i) can have several entries where fluxes for vertical surfaces might
1847!-- be stored. In the flat case, where no vertical walls exit, set indicies such that loop in
1848!-- diffusion routines will not be entered.
1849    DEALLOCATE ( surfaces%start_index )
1850    DEALLOCATE ( surfaces%end_index )
[3833]1851!
[4559]1852!-- Indices to locate surface element.
1853    DEALLOCATE ( surfaces%i )
1854    DEALLOCATE ( surfaces%j )
1855    DEALLOCATE ( surfaces%k )
[3833]1856!
[4559]1857!-- Surface-layer height
1858    DEALLOCATE ( surfaces%z_mo )
[3833]1859!
[4559]1860!-- Surface orientation
1861    DEALLOCATE ( surfaces%facing )
[3833]1862!
[4559]1863!-- Surface parallel wind velocity
1864    DEALLOCATE ( surfaces%uvw_abs )
[3833]1865!
[4593]1866!-- Precalculated ln(z/z0)
1867    DEALLOCATE ( surfaces%ln_z_z0  )
1868    DEALLOCATE ( surfaces%ln_z_z0h )
1869    DEALLOCATE ( surfaces%ln_z_z0q )
1870!
[4559]1871!-- Roughness
1872    DEALLOCATE ( surfaces%z0 )
1873    DEALLOCATE ( surfaces%z0h )
1874    DEALLOCATE ( surfaces%z0q )
[3833]1875!
[4559]1876!-- Friction velocity
1877    DEALLOCATE ( surfaces%us )
[3833]1878!
[4559]1879!-- Allocate Obukhov length and bulk Richardson number. Actually, at vertical surfaces these are
1880!-- only required for natural surfaces.
1881!-- For natural land surfaces
1882    DEALLOCATE( surfaces%ol )
1883    DEALLOCATE( surfaces%rib )
[3833]1884!
[4559]1885!-- Allocate arrays for surface momentum fluxes for u and v. For u at north- and south-facing
1886!-- surfaces, for v at east- and west-facing surfaces.
1887    DEALLOCATE ( surfaces%mom_flux_uv )
[3833]1888!
[4559]1889!-- Allocate array for surface momentum flux for w - wsus and wsvs
1890    DEALLOCATE ( surfaces%mom_flux_w )
[3833]1891!
[4559]1892!-- Allocate array for surface momentum flux for subgrid-scale tke wsus and wsvs; first index usvs
1893!-- or vsws, second index for wsus or wsvs, depending on surface.
1894    DEALLOCATE ( surfaces%mom_flux_tke )
[3833]1895!
[4559]1896!-- Characteristic temperature and surface flux of sensible heat
1897    DEALLOCATE ( surfaces%ts )
1898    DEALLOCATE ( surfaces%shf )
[3833]1899!
[4559]1900!-- Surface temperature
1901    DEALLOCATE ( surfaces%pt_surface )
[3833]1902!
[4559]1903!-- Characteristic humidity and surface flux of latent heat
1904    IF ( humidity )  THEN
1905       DEALLOCATE ( surfaces%qs )
1906       DEALLOCATE ( surfaces%qsws )
1907       DEALLOCATE ( surfaces%q_surface )
1908       DEALLOCATE ( surfaces%vpt_surface )
1909    ENDIF
[3833]1910!
[4559]1911!-- Characteristic scalar and surface flux of scalar
1912    IF ( passive_scalar )  THEN
1913       DEALLOCATE ( surfaces%ss )
1914       DEALLOCATE ( surfaces%ssws )
1915    ENDIF
[3833]1916!
[4559]1917!-- Scaling parameter (cs*) and surface flux of chemical species
1918    IF ( air_chemistry )  THEN
[4593]1919       DEALLOCATE ( surfaces%css )
1920       DEALLOCATE ( surfaces%cssws )
[4559]1921    ENDIF
[3833]1922!
[4559]1923!-- Arrays for storing potential temperature and mixing ratio at first grid level
1924    DEALLOCATE ( surfaces%pt1 )
1925    DEALLOCATE ( surfaces%qv1 )
1926    DEALLOCATE ( surfaces%vpt1 )
[3833]1927
[4559]1928    IF ( surf_bulk_cloud_model .AND. surf_microphysics_morrison)  THEN
1929       DEALLOCATE ( surfaces%qcs )
1930       DEALLOCATE ( surfaces%ncs )
1931       DEALLOCATE ( surfaces%qcsws )
1932       DEALLOCATE ( surfaces%ncsws )
1933    ENDIF
[3833]1934
[4559]1935    IF ( surf_bulk_cloud_model .AND. surf_microphysics_seifert)  THEN
1936       DEALLOCATE ( surfaces%qrs )
1937       DEALLOCATE ( surfaces%nrs )
1938       DEALLOCATE ( surfaces%qrsws )
1939       DEALLOCATE ( surfaces%nrsws )
1940    ENDIF
[4502]1941
[4559]1942    IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_phase)  THEN
1943       DEALLOCATE ( surfaces%qis )
1944       DEALLOCATE ( surfaces%nis )
1945       DEALLOCATE ( surfaces%qisws )
1946       DEALLOCATE ( surfaces%nisws )
1947    ENDIF
[4502]1948
[3833]1949!
[4559]1950!-- Salinity surface flux
1951    IF ( ocean_mode )  DEALLOCATE ( surfaces%sasws )
[3833]1952
[4559]1953 END SUBROUTINE deallocate_surface_attributes_v
[3833]1954
1955
[4559]1956!--------------------------------------------------------------------------------------------------!
[3833]1957! Description:
1958! ------------
[4502]1959!> Allocating memory for vertical surface types.
[4559]1960!--------------------------------------------------------------------------------------------------!
1961 SUBROUTINE allocate_surface_attributes_v( surfaces, nys_l, nyn_l, nxl_l, nxr_l )
[3833]1962
[4559]1963    IMPLICIT NONE
[3833]1964
[4559]1965    INTEGER(iwp) ::  nyn_l  !< north bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1966    INTEGER(iwp) ::  nys_l  !< south bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1967    INTEGER(iwp) ::  nxl_l  !< west bound of local 2d array start/end_index, is equal to nyn, except for restart-array
1968    INTEGER(iwp) ::  nxr_l  !< east bound of local 2d array start/end_index, is equal to nyn, except for restart-array
[3833]1969
[4559]1970    TYPE(surf_type) ::  surfaces !< respective surface type
[3833]1971
1972!
[4559]1973!-- Allocate arrays for start and end index of vertical surface type for each (j,i)-grid point. This
1974!-- is required in diffion_x, which is called for each (j,i). In order to find the location where
1975!-- the respective flux is store within the surface-type, start- and end-index are stored for each
1976!-- (j,i). For example, each (j,i) can have several entries where fluxes for vertical surfaces might
1977!-- be stored. In the flat case, where no vertical walls exit, set indicies such that loop in
1978!-- diffusion routines will not be entered.
1979    ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
1980    ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l) )
1981    surfaces%start_index = 0
1982    surfaces%end_index   = -1
[3833]1983!
[4559]1984!-- Indices to locate surface element.
1985    ALLOCATE ( surfaces%i(1:surfaces%ns) )
1986    ALLOCATE ( surfaces%j(1:surfaces%ns) )
1987    ALLOCATE ( surfaces%k(1:surfaces%ns) )
[3833]1988!
[4559]1989!-- Surface-layer height
1990    ALLOCATE ( surfaces%z_mo(1:surfaces%ns) )
[3833]1991!
[4559]1992!-- Surface orientation
1993    ALLOCATE ( surfaces%facing(1:surfaces%ns) )
[3833]1994!
[4559]1995!-- Surface parallel wind velocity
1996    ALLOCATE ( surfaces%uvw_abs(1:surfaces%ns) )
[3833]1997!
[4593]1998!-- Precalculated ln(z/z0)
1999    ALLOCATE( surfaces%ln_z_z0(1:surfaces%ns) )
2000    ALLOCATE( surfaces%ln_z_z0h(1:surfaces%ns) )
2001    ALLOCATE( surfaces%ln_z_z0q(1:surfaces%ns) )
2002!
[4559]2003!-- Roughness
2004    ALLOCATE ( surfaces%z0(1:surfaces%ns)  )
2005    ALLOCATE ( surfaces%z0h(1:surfaces%ns) )
2006    ALLOCATE ( surfaces%z0q(1:surfaces%ns) )
[3833]2007
2008!
[4559]2009!-- Friction velocity
2010    ALLOCATE ( surfaces%us(1:surfaces%ns) )
[3833]2011!
[4559]2012!-- Allocate Obukhov length and bulk Richardson number. Actually, at vertical surfaces these are
2013!-- only required for natural surfaces.
2014!-- For natural land surfaces
2015    ALLOCATE( surfaces%ol(1:surfaces%ns) )
2016    ALLOCATE( surfaces%rib(1:surfaces%ns) )
[3833]2017!
[4559]2018!-- Allocate arrays for surface momentum fluxes for u and v. For u at north- and south-facing
2019!-- surfaces, for v at east- and west-facing surfaces.
2020    ALLOCATE ( surfaces%mom_flux_uv(1:surfaces%ns) )
[3833]2021!
[4559]2022!-- Allocate array for surface momentum flux for w - wsus and wsvs
2023    ALLOCATE ( surfaces%mom_flux_w(1:surfaces%ns) )
[3833]2024!
[4559]2025!-- Allocate array for surface momentum flux for subgrid-scale tke wsus and wsvs; first index usvs
2026!-- or vsws, second index for wsus or wsvs, depending on surface.
2027    ALLOCATE ( surfaces%mom_flux_tke(0:1,1:surfaces%ns) )
[3833]2028!
[4559]2029!-- Characteristic temperature and surface flux of sensible heat
2030    ALLOCATE ( surfaces%ts(1:surfaces%ns) )
2031    ALLOCATE ( surfaces%shf(1:surfaces%ns) )
[3833]2032!
[4559]2033!-- Surface temperature
2034    ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) )
[3833]2035!
[4559]2036!-- Characteristic humidity and surface flux of latent heat
2037    IF ( humidity )  THEN
2038       ALLOCATE ( surfaces%qs(1:surfaces%ns) )
2039       ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
2040       ALLOCATE ( surfaces%q_surface(1:surfaces%ns) )
2041       ALLOCATE ( surfaces%vpt_surface(1:surfaces%ns) )
2042    ENDIF
[3833]2043!
[4559]2044!-- Characteristic scalar and surface flux of scalar
2045    IF ( passive_scalar )  THEN
2046       ALLOCATE ( surfaces%ss(1:surfaces%ns) )
2047       ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
2048    ENDIF
[3833]2049!
[4559]2050!-- Scaling parameter (cs*) and surface flux of chemical species
2051    IF ( air_chemistry )  THEN
[4593]2052       ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns) )
2053       ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
[4559]2054    ENDIF
[3833]2055!
[4559]2056!-- Arrays for storing potential temperature and mixing ratio at first grid level
2057    ALLOCATE ( surfaces%pt1(1:surfaces%ns) )
2058    ALLOCATE ( surfaces%qv1(1:surfaces%ns) )
2059    ALLOCATE ( surfaces%vpt1(1:surfaces%ns) )
[3833]2060
[4559]2061    IF ( surf_bulk_cloud_model .AND. surf_microphysics_morrison)  THEN
2062       ALLOCATE ( surfaces%qcs(1:surfaces%ns) )
2063       ALLOCATE ( surfaces%ncs(1:surfaces%ns) )
2064       ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
2065       ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
2066    ENDIF
[3833]2067
[4559]2068    IF ( surf_bulk_cloud_model .AND. surf_microphysics_seifert)  THEN
2069       ALLOCATE ( surfaces%qrs(1:surfaces%ns) )
2070       ALLOCATE ( surfaces%nrs(1:surfaces%ns) )
2071       ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
2072       ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
2073    ENDIF
[4502]2074
[4559]2075    IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_phase)  THEN
2076       ALLOCATE ( surfaces%qis(1:surfaces%ns) )
2077       ALLOCATE ( surfaces%nis(1:surfaces%ns) )
2078       ALLOCATE ( surfaces%qisws(1:surfaces%ns) )
2079       ALLOCATE ( surfaces%nisws(1:surfaces%ns) )
2080    ENDIF
[3833]2081!
[4559]2082!-- Salinity surface flux
2083    IF ( ocean_mode )  ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
[3833]2084
[4559]2085 END SUBROUTINE allocate_surface_attributes_v
[3833]2086
2087
[4559]2088!--------------------------------------------------------------------------------------------------!
[3833]2089! Description:
2090! ------------
[4502]2091!> Exit memory for vertical surface types.
[4559]2092!--------------------------------------------------------------------------------------------------!
[3833]2093#if defined( _OPENACC )
[4559]2094 SUBROUTINE exit_surface_attributes_v( surfaces )
[3833]2095
[4559]2096    IMPLICIT NONE
[3833]2097
[4559]2098    TYPE(surf_type) ::  surfaces  !< respective surface type
[3833]2099
[4559]2100    !$ACC EXIT DATA &
2101    !$ACC DELETE(surfaces%start_index(nys:nyn,nxl:nxr)) &
2102    !$ACC DELETE(surfaces%end_index(nys:nyn,nxl:nxr)) &
2103    !$ACC DELETE(surfaces%i(1:surfaces%ns)) &
2104    !$ACC DELETE(surfaces%j(1:surfaces%ns)) &
2105    !$ACC DELETE(surfaces%k(1:surfaces%ns)) &
2106    !$ACC DELETE(surfaces%uvw_abs(1:surfaces%ns)) &
[4593]2107    !$ACC DELETE(surfaces%ln_z_z0(1:surfaces%ns) ) &
2108    !$ACC DELETE(surfaces%ln_z_z0h(1:surfaces%ns) ) &
2109    !$ACC DELETE(surfaces%ln_z_z0q(1:surfaces%ns) ) &
[4559]2110    !$ACC DELETE(surfaces%z0(1:surfaces%ns)) &
[4594]2111    !$ACC DELETE(surfaces%z0h(1:surfaces%ns)) &
2112    !$ACC DELETE(surfaces%z0q(1:surfaces%ns)) &
[4559]2113    !$ACC DELETE(surfaces%rib(1:surfaces%ns)) &
2114    !$ACC DELETE(surfaces%mom_flux_uv(1:surfaces%ns)) &
2115    !$ACC DELETE(surfaces%mom_flux_w(1:surfaces%ns)) &
2116    !$ACC DELETE(surfaces%mom_flux_tke(0:1,1:surfaces%ns)) &
2117    !$ACC DELETE(surfaces%ts(1:surfaces%ns)) &
2118    !$ACC DELETE(surfaces%shf(1:surfaces%ns)) &
2119    !$ACC DELETE(surfaces%pt1(1:surfaces%ns)) &
2120    !$ACC DELETE(surfaces%qv1(1:surfaces%ns))
[3833]2121
[4559]2122 END SUBROUTINE exit_surface_attributes_v
[3833]2123#endif
2124
[4559]2125!--------------------------------------------------------------------------------------------------!
[3833]2126! Description:
2127! ------------
[4502]2128!> Enter memory for vertical surface types.
[4559]2129!--------------------------------------------------------------------------------------------------!
[3833]2130#if defined( _OPENACC )
[4559]2131 SUBROUTINE enter_surface_attributes_v( surfaces )
[4502]2132
[4559]2133    IMPLICIT NONE
[4502]2134
[4559]2135    TYPE(surf_type) ::  surfaces  !< respective surface type
[4502]2136
[4559]2137    !$ACC ENTER DATA &
2138    !$ACC COPYIN(surfaces%start_index(nys:nyn,nxl:nxr)) &
2139    !$ACC COPYIN(surfaces%end_index(nys:nyn,nxl:nxr)) &
2140    !$ACC COPYIN(surfaces%i(1:surfaces%ns)) &
2141    !$ACC COPYIN(surfaces%j(1:surfaces%ns)) &
2142    !$ACC COPYIN(surfaces%k(1:surfaces%ns)) &
2143    !$ACC COPYIN(surfaces%uvw_abs(1:surfaces%ns)) &
[4593]2144    !$ACC COPYIN(surfaces%ln_z_z0(1:surfaces%ns) ) &
2145    !$ACC COPYIN(surfaces%ln_z_z0h(1:surfaces%ns) ) &
2146    !$ACC COPYIN(surfaces%ln_z_z0q(1:surfaces%ns) ) &
[4559]2147    !$ACC COPYIN(surfaces%z0(1:surfaces%ns)) &
[4594]2148    !$ACC COPYIN(surfaces%z0h(1:surfaces%ns)) &
2149    !$ACC COPYIN(surfaces%z0q(1:surfaces%ns)) &
[4559]2150    !$ACC COPYIN(surfaces%rib(1:surfaces%ns)) &
2151    !$ACC COPYIN(surfaces%mom_flux_uv(1:surfaces%ns)) &
2152    !$ACC COPYIN(surfaces%mom_flux_w(1:surfaces%ns)) &
2153    !$ACC COPYIN(surfaces%mom_flux_tke(0:1,1:surfaces%ns)) &
2154    !$ACC COPYIN(surfaces%ts(1:surfaces%ns)) &
2155    !$ACC COPYIN(surfaces%shf(1:surfaces%ns)) &
2156    !$ACC COPYIN(surfaces%pt1(1:surfaces%ns)) &
2157    !$ACC COPYIN(surfaces%qv1(1:surfaces%ns))
[4502]2158
[4559]2159 END SUBROUTINE enter_surface_attributes_v
[3833]2160#endif
2161
[4559]2162!--------------------------------------------------------------------------------------------------!
[3833]2163! Description:
2164! ------------
[4559]2165!> Initialize surface elements, i.e. set initial values for surface fluxes, friction velocity,
2166!> calcuation of start/end indices, etc. Please note, further initialization concerning special
2167!> surface characteristics, e.g. soil- and vegatation type, building type, etc.,
2168!> is done in the land-surface and urban-surface module, respectively.
2169!--------------------------------------------------------------------------------------------------!
2170 SUBROUTINE init_surfaces
[3833]2171
[4559]2172    IMPLICIT NONE
[3833]2173
[4559]2174    INTEGER(iwp) ::  i  !< running index x-direction
2175    INTEGER(iwp) ::  j  !< running index y-direction
2176    INTEGER(iwp) ::  k  !< running index z-direction
[3833]2177
[4559]2178    INTEGER(iwp)  ::  start_index_lsm_h  !< dummy to determing local start index in surface type for given (j,i), for horizontal
2179                                         !< natural surfaces
2180    INTEGER(iwp)  ::  start_index_usm_h  !< dummy to determing local start index in surface type for given (j,i), for horizontal
2181                                         !< urban surfaces
[3833]2182
[4559]2183    INTEGER(iwp)  ::  num_lsm_h      !< current number of horizontal surface element, natural type
2184    INTEGER(iwp)  ::  num_lsm_h_kji  !< dummy to determing local end index in surface type for given (j,i), for for horizonal
2185                                     !< natural surfaces
2186    INTEGER(iwp)  ::  num_usm_h      !< current number of horizontal surface element, urban type
2187    INTEGER(iwp)  ::  num_usm_h_kji  !< dummy to determing local end index in surface type for given (j,i), for for horizonal urban
2188                                     !< surfaces
[3833]2189
[4559]2190    INTEGER(iwp), DIMENSION(0:2) ::  num_def_h          !< current number of horizontal surface element, default type
2191    INTEGER(iwp), DIMENSION(0:2) ::  num_def_h_kji      !< dummy to determing local end index in surface type for given (j,i),
2192                                                        !< for horizonal default surfaces
2193    INTEGER(iwp), DIMENSION(0:2) ::  start_index_def_h  !< dummy to determing local start index in surface type for given (j,i),
2194                                                        !< for horizontal default surfaces
[4502]2195
[4559]2196    INTEGER(iwp), DIMENSION(0:3) ::  num_def_v      !< current number of vertical surface element, default type
2197    INTEGER(iwp), DIMENSION(0:3) ::  num_def_v_kji  !< dummy to determing local end index in surface type for given (j,i),
2198                                                    !< for vertical default surfaces
2199    INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v      !< current number of vertical surface element, natural type
2200    INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v_kji  !< dummy to determing local end index in surface type for given (j,i),
2201                                                    !< for vertical natural surfaces
2202    INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v      !< current number of vertical surface element, urban type
2203    INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v_kji  !< dummy to determing local end index in surface type for given (j,i),
2204                                                    !< for vertical urban surfaces
[3833]2205
[4559]2206    INTEGER(iwp), DIMENSION(0:3) ::  start_index_def_v  !< dummy to determing local start index in surface type for given (j,i),
2207                                                        !< for vertical default surfaces
2208    INTEGER(iwp), DIMENSION(0:3) ::  start_index_lsm_v  !< dummy to determing local start index in surface type for given (j,i),
2209                                                        !< for vertical natural surfaces
2210    INTEGER(iwp), DIMENSION(0:3) ::  start_index_usm_v  !< dummy to determing local start index in surface type for given (j,i),
2211                                                        !< for vertical urban surfaces
[3833]2212
[4559]2213    LOGICAL ::  building             !< flag indicating building grid point
2214    LOGICAL ::  terrain              !< flag indicating natural terrain grid point
2215    LOGICAL ::  unresolved_building  !< flag indicating a grid point where actually a building is defined but not resolved by the
2216                                     !< vertical grid
[3833]2217!
[4559]2218!-- Set offset indices, i.e. index difference between surface element and surface-bounded grid point.
2219!-- Upward facing - no horizontal offsets
2220    surf_def_h(0:2)%ioff = 0
2221    surf_def_h(0:2)%joff = 0
[3833]2222
[4559]2223    surf_lsm_h%ioff = 0
2224    surf_lsm_h%joff = 0
[3833]2225
[4559]2226    surf_usm_h%ioff = 0
2227    surf_usm_h%joff = 0
[3833]2228!
[4559]2229!-- Upward facing vertical offsets
2230    surf_def_h(0)%koff = -1
2231    surf_lsm_h%koff    = -1
2232    surf_usm_h%koff    = -1
[3833]2233!
[4559]2234!-- Downward facing vertical offset
2235    surf_def_h(1:2)%koff = 1
[3833]2236!
[4559]2237!-- Vertical surfaces - no vertical offset
2238    surf_def_v(0:3)%koff = 0
2239    surf_lsm_v(0:3)%koff = 0
2240    surf_usm_v(0:3)%koff = 0
[3833]2241!
[4559]2242!-- North- and southward facing - no offset in x
2243    surf_def_v(0:1)%ioff = 0
2244    surf_lsm_v(0:1)%ioff = 0
2245    surf_usm_v(0:1)%ioff = 0
[3833]2246!
[4559]2247!-- Northward facing offset in y
2248    surf_def_v(0)%joff = -1
2249    surf_lsm_v(0)%joff = -1
2250    surf_usm_v(0)%joff = -1
[3833]2251!
[4559]2252!-- Southward facing offset in y
2253    surf_def_v(1)%joff = 1
2254    surf_lsm_v(1)%joff = 1
2255    surf_usm_v(1)%joff = 1
[3833]2256
2257!
[4559]2258!-- East- and westward facing - no offset in y
2259    surf_def_v(2:3)%joff = 0
2260    surf_lsm_v(2:3)%joff = 0
2261    surf_usm_v(2:3)%joff = 0
[3833]2262!
[4559]2263!-- Eastward facing offset in x
2264    surf_def_v(2)%ioff = -1
2265    surf_lsm_v(2)%ioff = -1
2266    surf_usm_v(2)%ioff = -1
[3833]2267!
[4559]2268!-- Westward facing offset in y
2269    surf_def_v(3)%ioff = 1
2270    surf_lsm_v(3)%ioff = 1
2271    surf_usm_v(3)%ioff = 1
[3833]2272
2273!
[4559]2274!-- Initialize surface attributes, store indicies, surfaces orientation, etc.,
2275    num_def_h(0:2) = 1
2276    num_def_v(0:3) = 1
[3833]2277
[4559]2278    num_lsm_h      = 1
2279    num_lsm_v(0:3) = 1
[3833]2280
[4559]2281    num_usm_h      = 1
2282    num_usm_v(0:3) = 1
[3833]2283
[4559]2284    start_index_def_h(0:2) = 1
2285    start_index_def_v(0:3) = 1
[3833]2286
[4559]2287    start_index_lsm_h      = 1
2288    start_index_lsm_v(0:3) = 1
[3833]2289
[4559]2290    start_index_usm_h      = 1
2291    start_index_usm_v(0:3) = 1
[3833]2292
[4559]2293    DO  i = nxl, nxr
2294       DO  j = nys, nyn
2295          num_def_h_kji = 0
2296          num_def_v_kji = 0
2297          num_lsm_h_kji = 0
2298          num_lsm_v_kji = 0
2299          num_usm_h_kji = 0
2300          num_usm_v_kji = 0
[3833]2301
[4559]2302          DO  k = nzb+1, nzt
[3833]2303!
[4559]2304!--          Check if current gridpoint belongs to the atmosphere
2305             IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
[3833]2306!
[4559]2307!--             Upward-facing surface. Distinguish between differet surface types.
2308!--             To do, think about method to flag natural and non-natural surfaces.
2309                IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i), 0 ) )  THEN
[3833]2310!
[4559]2311!--                Determine flags indicating terrain or building
2312                   terrain  = BTEST( wall_flags_total_0(k-1,j,i), 5 )  .OR.  topo_no_distinct
2313                   building = BTEST( wall_flags_total_0(k-1,j,i), 6 )  .OR.  topo_no_distinct
[4502]2314
[3833]2315!
[4559]2316!--                Unresolved_building indicates a surface with equal height as terrain but with a
2317!--                non-grid resolved building on top. These surfaces will be flagged as urban
2318!--                surfaces.
2319                   unresolved_building = BTEST( wall_flags_total_0(k-1,j,i), 5 ) .AND.             &
2320                                         BTEST( wall_flags_total_0(k-1,j,i), 6 )
[4159]2321!
[4559]2322!--                Natural surface type
2323                   IF ( land_surface  .AND.  terrain  .AND.  .NOT. unresolved_building )  THEN
2324                      CALL initialize_horizontal_surfaces( k, j, i, surf_lsm_h, num_lsm_h,         &
2325                                                           num_lsm_h_kji, .TRUE., .FALSE. )
[3833]2326!
[4559]2327!--                Urban surface tpye
2328                   ELSEIF ( urban_surface  .AND.  building )  THEN
2329                      CALL initialize_horizontal_surfaces( k, j, i, surf_usm_h, num_usm_h,         &
2330                                                           num_usm_h_kji, .TRUE., .FALSE. )
[3833]2331!
[4559]2332!--                Default surface type
2333                   ELSE
2334                      CALL initialize_horizontal_surfaces( k, j, i, surf_def_h(0), num_def_h(0),   &
2335                                                           num_def_h_kji(0), .TRUE., .FALSE. )
[4502]2336                   ENDIF
[4559]2337                ENDIF
[3833]2338!
[4559]2339!--             Downward-facing surface, first, model top. Please note, for the moment,
2340!--             downward-facing surfaces are always of default type
2341                IF ( k == nzt  .AND.  use_top_fluxes )  THEN
2342                   CALL initialize_top( k, j, i, surf_def_h(2), num_def_h(2), num_def_h_kji(2) )
[3833]2343!
[4559]2344!--             Check for any other downward-facing surface. So far only for default surface type.
2345                ELSEIF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i), 0 ) )  THEN
2346                   CALL initialize_horizontal_surfaces( k, j, i, surf_def_h(1), num_def_h(1),      &
2347                                                        num_def_h_kji(1), .FALSE., .TRUE. )
2348                ENDIF
[3833]2349!
[4559]2350!--             Check for vertical walls and, if required, initialize it.
2351!               Start with northward-facing surface.
2352                IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i), 0 ) )  THEN
[3833]2353!
[4559]2354!--                Determine flags indicating terrain or building
2355                   terrain  = BTEST( wall_flags_total_0(k,j-1,i), 5 )  .OR.  topo_no_distinct
2356                   building = BTEST( wall_flags_total_0(k,j-1,i), 6 )  .OR.  topo_no_distinct
[4159]2357
[4559]2358                   unresolved_building = BTEST( wall_flags_total_0(k,j-1,i), 5 )  .AND.            &
2359                                         BTEST( wall_flags_total_0(k,j-1,i), 6 )
[4502]2360
[4559]2361                   IF ( land_surface  .AND.  terrain  .AND.  .NOT. unresolved_building )  THEN
2362                      CALL initialize_vertical_surfaces( k, j, i, surf_lsm_v(0), num_lsm_v(0),     &
2363                                                         num_lsm_v_kji(0), .FALSE., .FALSE.,       &
2364                                                         .FALSE., .TRUE. )
2365
2366                   ELSEIF ( urban_surface  .AND.  building )  THEN
2367                      CALL initialize_vertical_surfaces( k, j, i, surf_usm_v(0), num_usm_v(0),     &
2368                                                         num_usm_v_kji(0), .FALSE., .FALSE.,       &
2369                                                         .FALSE., .TRUE. )
2370                   ELSE
2371                      CALL initialize_vertical_surfaces( k, j, i, surf_def_v(0), num_def_v(0),     &
2372                                                         num_def_v_kji(0), .FALSE., .FALSE.,       &
2373                                                         .FALSE., .TRUE. )
[3833]2374                   ENDIF
[4559]2375                ENDIF
[3833]2376!
[4559]2377!--             Southward-facing surface
2378                IF ( .NOT. BTEST( wall_flags_total_0(k,j+1,i), 0 ) )  THEN
[3833]2379!
[4559]2380!--                Determine flags indicating terrain or building
2381                   terrain  = BTEST( wall_flags_total_0(k,j+1,i), 5 )  .OR.  topo_no_distinct
2382                   building = BTEST( wall_flags_total_0(k,j+1,i), 6 )  .OR.  topo_no_distinct
[4502]2383
[4559]2384                   unresolved_building = BTEST( wall_flags_total_0(k,j+1,i), 5 )  .AND.            &
2385                                         BTEST( wall_flags_total_0(k,j+1,i), 6 )
[4502]2386
[4559]2387                   IF ( land_surface  .AND.  terrain  .AND.  .NOT. unresolved_building )  THEN
2388                      CALL initialize_vertical_surfaces( k, j, i, surf_lsm_v(1), num_lsm_v(1),     &
2389                                                         num_lsm_v_kji(1), .FALSE., .FALSE.,       &
2390                                                         .TRUE., .FALSE. )
2391
2392                   ELSEIF ( urban_surface  .AND.  building )  THEN
2393                      CALL initialize_vertical_surfaces( k, j, i, surf_usm_v(1), num_usm_v(1),     &
2394                                                         num_usm_v_kji(1), .FALSE., .FALSE.,       &
2395                                                         .TRUE., .FALSE. )
2396                   ELSE
2397                      CALL initialize_vertical_surfaces( k, j, i, surf_def_v(1), num_def_v(1),     &
2398                                                         num_def_v_kji(1), .FALSE., .FALSE.,       &
2399                                                         .TRUE., .FALSE. )
[3833]2400                   ENDIF
[4559]2401                ENDIF
[3833]2402!
[4559]2403!--             Eastward-facing surface
2404                IF ( .NOT. BTEST( wall_flags_total_0(k,j,i-1), 0 ) )  THEN
[3833]2405!
[4559]2406!--                Determine flags indicating terrain or building
2407                   terrain  = BTEST( wall_flags_total_0(k,j,i-1), 5 )  .OR.  topo_no_distinct
2408                   building = BTEST( wall_flags_total_0(k,j,i-1), 6 )  .OR.  topo_no_distinct
[4502]2409
[4559]2410                   unresolved_building = BTEST( wall_flags_total_0(k,j,i-1), 5 )  .AND.            &
2411                                         BTEST( wall_flags_total_0(k,j,i-1), 6 )
[4502]2412
[4559]2413                   IF ( land_surface  .AND.  terrain  .AND.  .NOT. unresolved_building )  THEN
2414                      CALL initialize_vertical_surfaces( k, j, i, surf_lsm_v(2), num_lsm_v(2),     &
2415                                                         num_lsm_v_kji(2), .TRUE., .FALSE.,        &
2416                                                         .FALSE., .FALSE. )
2417
2418                   ELSEIF ( urban_surface  .AND.  building )  THEN
2419                      CALL initialize_vertical_surfaces( k, j, i, surf_usm_v(2), num_usm_v(2),     &
2420                                                         num_usm_v_kji(2), .TRUE., .FALSE.,        &
2421                                                         .FALSE., .FALSE. )
2422                   ELSE
2423                      CALL initialize_vertical_surfaces( k, j, i, surf_def_v(2), num_def_v(2),     &
2424                                                         num_def_v_kji(2), .TRUE., .FALSE.,        &
2425                                                         .FALSE., .FALSE. )
[4502]2426                   ENDIF
[4559]2427                ENDIF
[4502]2428!
[4559]2429!--             Westward-facing surface
2430                IF ( .NOT. BTEST( wall_flags_total_0(k,j,i+1), 0 ) )  THEN
[3833]2431!
[4559]2432!--                Determine flags indicating terrain or building
2433                   terrain  = BTEST( wall_flags_total_0(k,j,i+1), 5 )  .OR.  topo_no_distinct
2434                   building = BTEST( wall_flags_total_0(k,j,i+1), 6 )  .OR.  topo_no_distinct
[4502]2435
[4559]2436                   unresolved_building = BTEST( wall_flags_total_0(k,j,i+1), 5 )  .AND.            &
2437                                         BTEST( wall_flags_total_0(k,j,i+1), 6 )
[4502]2438
[4559]2439                   IF ( land_surface  .AND.  terrain  .AND.  .NOT. unresolved_building )  THEN
2440                      CALL initialize_vertical_surfaces( k, j, i, surf_lsm_v(3), num_lsm_v(3),     &
2441                                                         num_lsm_v_kji(3), .FALSE., .TRUE.,        &
2442                                                        .FALSE., .FALSE. )
2443
2444                   ELSEIF ( urban_surface  .AND.  building )  THEN
2445                      CALL initialize_vertical_surfaces( k, j, i, surf_usm_v(3), num_usm_v(3),     &
2446                                                         num_usm_v_kji(3), .FALSE., .TRUE.,        &
2447                                                        .FALSE., .FALSE. )
2448                   ELSE
2449                      CALL initialize_vertical_surfaces( k, j, i, surf_def_v(3), num_def_v(3),     &
2450                                                         num_def_v_kji(3), .FALSE., .TRUE.,        &
2451                                                        .FALSE., .FALSE. )
[3833]2452                   ENDIF
2453                ENDIF
[4559]2454             ENDIF
[3833]2455
[4502]2456
[4559]2457          ENDDO
[3833]2458!
[4559]2459!--       Determine start- and end-index at grid point (j,i). Also, for horizontal surfaces more
2460!--       than 1 horizontal surface element can exist at grid point (j,i) if overhanging structures
2461!--       are present.
2462!--       Upward-facing surfaces
2463          surf_def_h(0)%start_index(j,i) = start_index_def_h(0)
2464          surf_def_h(0)%end_index(j,i)   = surf_def_h(0)%start_index(j,i) + num_def_h_kji(0) - 1
2465          start_index_def_h(0)           = surf_def_h(0)%end_index(j,i) + 1
[3833]2466!
[4559]2467!--       ATTENTION:
2468!--       Workaround to prevent vectorization bug on NEC Aurora
2469          IF ( start_index_def_h(0) < -99999 )  THEN
2470             PRINT*, 'i=', i, ' j=',j, ' s=',surf_def_h(0)%start_index(j,i),                       &
2471                    ' e=', surf_def_h(0)%end_index(j,i)
2472          ENDIF
[4366]2473!
[4559]2474!--       Downward-facing surfaces, except model top
2475          surf_def_h(1)%start_index(j,i) = start_index_def_h(1)
2476          surf_def_h(1)%end_index(j,i)   = surf_def_h(1)%start_index(j,i) + num_def_h_kji(1) - 1
2477          start_index_def_h(1)           = surf_def_h(1)%end_index(j,i) + 1
[3833]2478!
[4559]2479!--       Downward-facing surfaces -- model top fluxes
2480          surf_def_h(2)%start_index(j,i) = start_index_def_h(2)
2481          surf_def_h(2)%end_index(j,i)   = surf_def_h(2)%start_index(j,i) + num_def_h_kji(2) - 1
2482          start_index_def_h(2)           = surf_def_h(2)%end_index(j,i) + 1
[3833]2483!
[4559]2484!--       Horizontal natural land surfaces
2485          surf_lsm_h%start_index(j,i)    = start_index_lsm_h
2486          surf_lsm_h%end_index(j,i)      = surf_lsm_h%start_index(j,i) + num_lsm_h_kji - 1
2487          start_index_lsm_h              = surf_lsm_h%end_index(j,i) + 1
[3833]2488!
[4559]2489!--       Horizontal urban surfaces
2490          surf_usm_h%start_index(j,i)    = start_index_usm_h
2491          surf_usm_h%end_index(j,i)      = surf_usm_h%start_index(j,i) + num_usm_h_kji - 1
2492          start_index_usm_h              = surf_usm_h%end_index(j,i) + 1
[3833]2493
2494!
[4559]2495!--       Vertical surfaces - Default type
2496          surf_def_v(0)%start_index(j,i) = start_index_def_v(0)
2497          surf_def_v(1)%start_index(j,i) = start_index_def_v(1)
2498          surf_def_v(2)%start_index(j,i) = start_index_def_v(2)
2499          surf_def_v(3)%start_index(j,i) = start_index_def_v(3)
2500          surf_def_v(0)%end_index(j,i)   = start_index_def_v(0) + num_def_v_kji(0) - 1
2501          surf_def_v(1)%end_index(j,i)   = start_index_def_v(1) + num_def_v_kji(1) - 1
2502          surf_def_v(2)%end_index(j,i)   = start_index_def_v(2) + num_def_v_kji(2) - 1
2503          surf_def_v(3)%end_index(j,i)   = start_index_def_v(3) + num_def_v_kji(3) - 1
2504          start_index_def_v(0)           = surf_def_v(0)%end_index(j,i) + 1
2505          start_index_def_v(1)           = surf_def_v(1)%end_index(j,i) + 1
2506          start_index_def_v(2)           = surf_def_v(2)%end_index(j,i) + 1
2507          start_index_def_v(3)           = surf_def_v(3)%end_index(j,i) + 1
[3833]2508!
[4559]2509!--       Natural type
2510          surf_lsm_v(0)%start_index(j,i) = start_index_lsm_v(0)
2511          surf_lsm_v(1)%start_index(j,i) = start_index_lsm_v(1)
2512          surf_lsm_v(2)%start_index(j,i) = start_index_lsm_v(2)
2513          surf_lsm_v(3)%start_index(j,i) = start_index_lsm_v(3)
2514          surf_lsm_v(0)%end_index(j,i)   = start_index_lsm_v(0) + num_lsm_v_kji(0) - 1
2515          surf_lsm_v(1)%end_index(j,i)   = start_index_lsm_v(1) + num_lsm_v_kji(1) - 1
2516          surf_lsm_v(2)%end_index(j,i)   = start_index_lsm_v(2) + num_lsm_v_kji(2) - 1
2517          surf_lsm_v(3)%end_index(j,i)   = start_index_lsm_v(3) + num_lsm_v_kji(3) - 1
2518          start_index_lsm_v(0)           = surf_lsm_v(0)%end_index(j,i) + 1
2519          start_index_lsm_v(1)           = surf_lsm_v(1)%end_index(j,i) + 1
2520          start_index_lsm_v(2)           = surf_lsm_v(2)%end_index(j,i) + 1
2521          start_index_lsm_v(3)           = surf_lsm_v(3)%end_index(j,i) + 1
[3833]2522!
[4559]2523!--       Urban type
2524          surf_usm_v(0)%start_index(j,i) = start_index_usm_v(0)
2525          surf_usm_v(1)%start_index(j,i) = start_index_usm_v(1)
2526          surf_usm_v(2)%start_index(j,i) = start_index_usm_v(2)
2527          surf_usm_v(3)%start_index(j,i) = start_index_usm_v(3)
2528          surf_usm_v(0)%end_index(j,i)   = start_index_usm_v(0) + num_usm_v_kji(0) - 1
2529          surf_usm_v(1)%end_index(j,i)   = start_index_usm_v(1) + num_usm_v_kji(1) - 1
2530          surf_usm_v(2)%end_index(j,i)   = start_index_usm_v(2) + num_usm_v_kji(2) - 1
2531          surf_usm_v(3)%end_index(j,i)   = start_index_usm_v(3) + num_usm_v_kji(3) - 1
2532          start_index_usm_v(0)           = surf_usm_v(0)%end_index(j,i) + 1
2533          start_index_usm_v(1)           = surf_usm_v(1)%end_index(j,i) + 1
2534          start_index_usm_v(2)           = surf_usm_v(2)%end_index(j,i) + 1
2535          start_index_usm_v(3)           = surf_usm_v(3)%end_index(j,i) + 1
[3833]2536
2537
2538       ENDDO
[4559]2539    ENDDO
[3833]2540
[4559]2541 CONTAINS
[3833]2542
[4559]2543!--------------------------------------------------------------------------------------------------!
[3833]2544! Description:
2545! ------------
[4559]2546!> Initialize horizontal surface elements, upward- and downward-facing. Note, horizontal surface
2547!> type also comprises model-top fluxes, which are, initialized in a different routine.
2548!--------------------------------------------------------------------------------------------------!
2549 SUBROUTINE initialize_horizontal_surfaces( k, j, i, surf, num_h, num_h_kji, upward_facing,        &
2550                                            downward_facing )
[3833]2551
[4559]2552    IMPLICIT NONE
[3833]2553
[4559]2554    INTEGER(iwp)  ::  i          !< running index x-direction
2555    INTEGER(iwp)  ::  j          !< running index y-direction
2556    INTEGER(iwp)  ::  k          !< running index z-direction
2557    INTEGER(iwp)  ::  num_h      !< current number of surface element
2558    INTEGER(iwp)  ::  num_h_kji  !< dummy increment
2559    INTEGER(iwp)  ::  lsp        !< running index chemical species
2560    INTEGER(iwp)  ::  lsp_pr     !< running index chemical species??
[3833]2561
[4559]2562    LOGICAL       ::  upward_facing    !< flag indicating upward-facing surface
2563    LOGICAL       ::  downward_facing  !< flag indicating downward-facing surface
[3833]2564
[4559]2565    TYPE(surf_type) ::  surf  !< respective surface type
[3833]2566
2567!
[4559]2568!-- Store indices of respective surface element
2569    surf%i(num_h) = i
2570    surf%j(num_h) = j
2571    surf%k(num_h) = k
[3833]2572!
[4559]2573!-- Surface orientation, bit 0 is set to 1 for upward-facing surfaces, bit 1 is for downward-facing
2574!-- surfaces.
2575    IF ( upward_facing   )  surf%facing(num_h) = IBSET( surf%facing(num_h), 0 )
2576    IF ( downward_facing )  surf%facing(num_h) = IBSET( surf%facing(num_h), 1 )
[3833]2577!
[4559]2578!-- Initialize surface-layer height
2579    IF ( upward_facing )  THEN
2580       surf%z_mo(num_h) = zu(k) - zw(k-1)
2581    ELSE
2582       surf%z_mo(num_h) = zw(k) - zu(k)
2583    ENDIF
[4502]2584
[4559]2585    surf%z0(num_h)  = roughness_length
2586    surf%z0h(num_h) = z0h_factor * roughness_length
2587    surf%z0q(num_h) = z0h_factor * roughness_length
[3833]2588!
[4559]2589!-- Initialization in case of 1D pre-cursor run
2590    IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
2591       IF ( .NOT. constant_diffusion )  THEN
2592          IF ( constant_flux_layer )  THEN
[4586]2593             surf%ol(num_h)   = surf%z_mo(num_h) / ( ri1d(nzb+1) + 1.0E-20_wp )
[4559]2594             surf%us(num_h)   = us1d
2595             surf%usws(num_h) = usws1d
2596             surf%vsws(num_h) = vsws1d
2597          ELSE
2598             surf%ol(num_h)   = surf%z_mo(num_h) / zeta_min
2599             surf%us(num_h)   = 0.0_wp
2600             surf%usws(num_h) = 0.0_wp
2601             surf%vsws(num_h) = 0.0_wp
2602          ENDIF
2603       ELSE
2604          surf%ol(num_h)   = surf%z_mo(num_h) / zeta_min
2605          surf%us(num_h)   = 0.0_wp
2606          surf%usws(num_h) = 0.0_wp
2607          surf%vsws(num_h) = 0.0_wp
2608       ENDIF
[3833]2609!
[4559]2610!-- Initialization in all other cases
2611    ELSE
[3833]2612
[4559]2613       surf%ol(num_h) = surf%z_mo(num_h) / zeta_min
[3833]2614!
[4559]2615!--    Very small number is required for calculation of Obukhov length at first timestep
2616       surf%us(num_h)   = 1E-30_wp
2617       surf%usws(num_h) = 0.0_wp
2618       surf%vsws(num_h) = 0.0_wp
[4502]2619
[4559]2620    ENDIF
[3833]2621
[4559]2622    surf%rib(num_h)     = 0.0_wp
2623    surf%uvw_abs(num_h) = 0.0_wp
[4593]2624!
2625!-- Initialize ln(z/z0)
2626    surf%ln_z_z0(num_h)  = LOG( surf%z_mo(num_h) / surf%z0(num_h) )
2627    surf%ln_z_z0h(num_h) = LOG( surf%z_mo(num_h) / surf%z0h(num_h) )
2628    surf%ln_z_z0q(num_h) = LOG( surf%z_mo(num_h) / surf%z0q(num_h) )
[3833]2629
[4559]2630    IF ( .NOT. constant_diffusion )  THEN
2631       surf%u_0(num_h) = 0.0_wp
2632       surf%v_0(num_h) = 0.0_wp
2633    ENDIF
[3833]2634
[4559]2635    surf%ts(num_h) = 0.0_wp
[3833]2636!
[4559]2637!-- Set initial value for surface temperature
2638    surf%pt_surface(num_h) = pt_surface
[4502]2639
[4559]2640    IF ( humidity )  THEN
2641       surf%qs(num_h)   = 0.0_wp
2642       IF ( surf_bulk_cloud_model .AND. surf_microphysics_morrison)  THEN
2643          surf%qcs(num_h) = 0.0_wp
2644          surf%ncs(num_h) = 0.0_wp
[4502]2645
[4559]2646          surf%qcsws(num_h) = 0.0_wp
2647          surf%ncsws(num_h) = 0.0_wp
[3833]2648
[4559]2649       ENDIF
2650       IF ( surf_bulk_cloud_model .AND. surf_microphysics_seifert)  THEN
2651          surf%qrs(num_h) = 0.0_wp
2652          surf%nrs(num_h) = 0.0_wp
[4502]2653
[4559]2654          surf%qrsws(num_h) = 0.0_wp
2655          surf%nrsws(num_h) = 0.0_wp
[3833]2656
[4559]2657          surf%pt1(num_h)  = 0.0_wp
2658          surf%qv1(num_h)  = 0.0_wp
2659          surf%vpt1(num_h) = 0.0_wp
[4502]2660
[4559]2661       ENDIF
[4502]2662
[4559]2663       IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_phase)  THEN
2664          surf%qis(num_h) = 0.0_wp
2665          surf%nis(num_h) = 0.0_wp
[4502]2666
[4559]2667          surf%qisws(num_h) = 0.0_wp
2668          surf%nisws(num_h) = 0.0_wp
2669       ENDIF
[4502]2670
2671
[4559]2672       surf%q_surface(num_h)   = q_surface
2673       surf%vpt_surface(num_h) = surf%pt_surface(num_h) *                                          &
2674                                ( 1.0_wp + 0.61_wp * surf%q_surface(num_h) )
2675    ENDIF
[3833]2676
[4559]2677    IF ( passive_scalar )  surf%ss(num_h) = 0.0_wp
[3833]2678
[4559]2679    DO  lsp = 1, nvar
2680       IF ( air_chemistry )  surf%css(lsp,num_h)   = 0.0_wp
[3833]2681!
[4559]2682!--    Ensure that fluxes of compounds which are not specified in namelist are all zero
2683!--    --> kanani: revise
2684       IF ( air_chemistry )  surf%cssws(lsp,num_h) = 0.0_wp
2685    ENDDO
[3833]2686!
[4559]2687!-- Inititalize surface fluxes of sensible and latent heat, as well as passive scalar
2688    IF ( use_surface_fluxes )  THEN
[3833]2689
[4559]2690       IF ( upward_facing )  THEN
2691          IF ( constant_heatflux )  THEN
[4502]2692!
[4559]2693!--          Initialize surface heatflux. However, skip this for now if random_heatflux is set.
2694!--          This case, shf is initialized later.
2695             IF ( .NOT. random_heatflux )  THEN
2696                surf%shf(num_h) = surface_heatflux * heatflux_input_conversion(k-1)
[3833]2697!
[4559]2698!--             Check if surface heat flux might be replaced by prescribed wall heatflux
2699                IF ( k-1 /= 0 )  THEN
2700                   surf%shf(num_h) = wall_heatflux(0) * heatflux_input_conversion(k-1)
2701                ENDIF
2702             ENDIF
2703          ELSE
2704             surf%shf(num_h) = 0.0_wp
2705          ENDIF
[3833]2706!
[4559]2707!--    Set heat-flux at downward-facing surfaces
2708       ELSE
2709          surf%shf(num_h) = wall_heatflux(5) * heatflux_input_conversion(k)
2710       ENDIF
[3833]2711
[4559]2712       IF ( humidity )  THEN
2713          IF ( upward_facing )  THEN
2714             IF ( constant_waterflux )  THEN
2715                surf%qsws(num_h) = surface_waterflux * waterflux_input_conversion(k-1)
2716                IF ( k-1 /= 0 )  THEN
2717                   surf%qsws(num_h) = wall_humidityflux(0) * waterflux_input_conversion(k-1)
[3833]2718                ENDIF
[4559]2719             ELSE
2720                surf%qsws(num_h) = 0.0_wp
2721             ENDIF
2722          ELSE
2723             surf%qsws(num_h) = wall_humidityflux(5) * waterflux_input_conversion(k)
2724          ENDIF
2725       ENDIF
[3833]2726
[4559]2727       IF ( passive_scalar )  THEN
2728          IF ( upward_facing )  THEN
2729             IF ( constant_scalarflux )  THEN
2730                surf%ssws(num_h) = surface_scalarflux  * rho_air_zw(k-1)
2731
2732                IF ( k-1 /= 0 )  surf%ssws(num_h) = wall_scalarflux(0) * rho_air_zw(k-1)
2733             ELSE
2734                surf%ssws(num_h) = 0.0_wp
2735             ENDIF
2736          ELSE
2737             surf%ssws(num_h) = wall_scalarflux(5) * rho_air_zw(k)
2738          ENDIF
2739       ENDIF
2740
2741       IF ( air_chemistry )  THEN
2742          lsp_pr = 1
2743          DO  WHILE ( TRIM( surface_csflux_name( lsp_pr ) ) /= 'novalue' ) !<'novalue' is the default
2744             DO  lsp = 1, nvar
2745!
2746!--             Assign surface flux for each variable species
2747                IF ( TRIM( spc_names(lsp) ) == TRIM( surface_csflux_name(lsp_pr) ) )  THEN
[3833]2748                   IF ( upward_facing )  THEN
[4559]2749                      IF ( constant_csflux(lsp_pr) )  THEN
2750                         surf%cssws(lsp,num_h) = surface_csflux(lsp_pr) * rho_air_zw(k-1)
[3833]2751
[4559]2752                         IF ( k-1 /= 0 )  surf%cssws(lsp,num_h) = wall_csflux(lsp,0) *             &
2753                                                                  rho_air_zw(k-1)
[3833]2754                      ELSE
[4559]2755                         surf%cssws(lsp,num_h) = 0.0_wp
[3833]2756                      ENDIF
2757                   ELSE
[4559]2758                      surf%cssws(lsp,num_h) = wall_csflux(lsp,5) * rho_air_zw(k)
[3833]2759                   ENDIF
2760                ENDIF
[4559]2761             ENDDO
2762             lsp_pr = lsp_pr + 1
2763          ENDDO
2764       ENDIF
[3833]2765
[4559]2766       IF ( ocean_mode )  THEN
2767          IF ( upward_facing )  THEN
2768             surf%sasws(num_h) = bottom_salinityflux * rho_air_zw(k-1)
2769          ELSE
2770             surf%sasws(num_h) = 0.0_wp
2771          ENDIF
2772       ENDIF
2773    ENDIF
[3833]2774!
[4559]2775!-- Increment surface indices
2776    num_h     = num_h + 1
2777    num_h_kji = num_h_kji + 1
[3833]2778
2779
[4559]2780 END SUBROUTINE initialize_horizontal_surfaces
[3833]2781
2782
[4559]2783!--------------------------------------------------------------------------------------------------!
[3833]2784! Description:
2785! ------------
[4559]2786!> Initialize model-top fluxes. Currently, only the heatflux and salinity flux can be prescribed,
2787!> latent flux is zero in this case!
2788!--------------------------------------------------------------------------------------------------!
2789 SUBROUTINE initialize_top( k, j, i, surf, num_h, num_h_kji )
[3833]2790
[4559]2791    IMPLICIT NONE
[3833]2792
[4559]2793    INTEGER(iwp)  ::  i          !< running index x-direction
2794    INTEGER(iwp)  ::  j          !< running index y-direction
2795    INTEGER(iwp)  ::  k          !< running index z-direction
2796    INTEGER(iwp)  ::  num_h      !< current number of surface element
2797    INTEGER(iwp)  ::  num_h_kji  !< dummy increment
2798    INTEGER(iwp)  ::  lsp        !< running index for chemical species
[3833]2799
[4559]2800    TYPE( surf_type ) ::  surf   !< respective surface type
[3833]2801!
[4559]2802!-- Store indices of respective surface element
2803    surf%i(num_h) = i
2804    surf%j(num_h) = j
2805    surf%k(num_h) = k
[3833]2806!
[4559]2807!-- Initialize top heat flux
2808    IF ( constant_top_heatflux )  surf%shf(num_h) = top_heatflux * heatflux_input_conversion(nzt+1)
[3833]2809!
[4559]2810!-- Initialization in case of a coupled model run
2811    IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
2812       surf%shf(num_h) = 0.0_wp
2813       surf%qsws(num_h) = 0.0_wp
2814    ENDIF
[3833]2815!
[4559]2816!-- Prescribe latent heat flux at the top
2817    IF ( humidity )  THEN
2818       surf%qsws(num_h) = 0.0_wp
2819       IF ( surf_bulk_cloud_model  .AND.  surf_microphysics_morrison )  THEN
2820          surf%ncsws(num_h) = 0.0_wp
2821          surf%qcsws(num_h) = 0.0_wp
2822       ENDIF
2823       IF ( surf_bulk_cloud_model  .AND.  surf_microphysics_seifert )  THEN
2824          surf%nrsws(num_h) = 0.0_wp
2825          surf%qrsws(num_h) = 0.0_wp
2826       ENDIF
2827       IF ( surf_bulk_cloud_model  .AND.  surf_microphysics_ice_phase )  THEN
2828          surf%nisws(num_h) = 0.0_wp
2829          surf%qisws(num_h) = 0.0_wp
2830       ENDIF
2831    ENDIF
[3833]2832!
[4559]2833!-- Prescribe top scalar flux
2834    IF ( passive_scalar .AND. constant_top_scalarflux )  surf%ssws(num_h) = top_scalarflux *       &
2835                                                                            rho_air_zw(nzt+1)
[3833]2836!
[4559]2837!-- Prescribe top chemical species' flux
2838    DO  lsp = 1, nvar
2839       IF ( air_chemistry  .AND.  constant_top_csflux(lsp) )  THEN
2840          surf%cssws(lsp,num_h) = top_csflux(lsp) * rho_air_zw(nzt+1)
2841       ENDIF
2842    ENDDO
[3833]2843!
[4559]2844!-- Prescribe top salinity flux
2845    IF ( ocean_mode .AND. constant_top_salinityflux)  surf%sasws(num_h) = top_salinityflux *       &
2846                                                                          rho_air_zw(nzt+1)
[3833]2847!
[4559]2848!-- Top momentum fluxes
2849    IF ( constant_top_momentumflux )  THEN
2850       surf%usws(num_h) = top_momentumflux_u * momentumflux_input_conversion(nzt+1)
2851       surf%vsws(num_h) = top_momentumflux_v * momentumflux_input_conversion(nzt+1)
2852    ENDIF
[3833]2853!
[4559]2854!-- Increment surface indices
2855    num_h     = num_h + 1
2856    num_h_kji = num_h_kji + 1
[3833]2857
2858
[4559]2859 END SUBROUTINE initialize_top
[3833]2860
2861
[4559]2862!--------------------------------------------------------------------------------------------------!
[3833]2863! Description:
2864! ------------
[4502]2865!> Initialize vertical surface elements.
[4559]2866!--------------------------------------------------------------------------------------------------!
2867 SUBROUTINE initialize_vertical_surfaces( k, j, i, surf, num_v, num_v_kji, east_facing,            &
2868                                          west_facing, south_facing, north_facing )
[3833]2869
[4559]2870    IMPLICIT NONE
[3833]2871
[4559]2872    INTEGER(iwp) ::  component  !< index of wall_fluxes_ array for respective orientation
2873    INTEGER(iwp) ::  i          !< running index x-direction
2874    INTEGER(iwp) ::  j          !< running index x-direction
2875    INTEGER(iwp) ::  k          !< running index x-direction
2876    INTEGER(iwp) ::  num_v      !< current number of surface element
2877    INTEGER(iwp) ::  num_v_kji  !< current number of surface element at (j,i)
2878    INTEGER(iwp) ::  lsp        !< running index for chemical species
[3833]2879
2880
[4559]2881    LOGICAL ::  east_facing   !< flag indicating east-facing surfaces
2882    LOGICAL ::  north_facing  !< flag indicating north-facing surfaces
2883    LOGICAL ::  south_facing  !< flag indicating south-facing surfaces
2884    LOGICAL ::  west_facing   !< flag indicating west-facing surfaces
[3833]2885
[4559]2886    TYPE( surf_type ) ::  surf  !< respective surface type
[3833]2887
2888!
[4559]2889!-- Store indices of respective wall element
2890    surf%i(num_v) = i
2891    surf%j(num_v) = j
2892    surf%k(num_v) = k
[3833]2893!
[4559]2894!-- Initialize surface-layer height, or more precisely, distance to surface
2895    IF ( north_facing  .OR.  south_facing )  THEN
2896       surf%z_mo(num_v) = 0.5_wp * dy
2897    ELSE
2898       surf%z_mo(num_v) = 0.5_wp * dx
2899    ENDIF
[3833]2900
[4559]2901    surf%facing(num_v) = 0
[3833]2902!
[4559]2903!-- Surface orientation. Moreover, set component id to map wall_heatflux, etc., on surface type
2904!-- (further below)
2905    IF ( north_facing )  THEN
2906       surf%facing(num_v) = 5 !IBSET( surf%facing(num_v), 0 )
2907       component          = 4
2908    ENDIF
[3833]2909
[4559]2910    IF ( south_facing )  THEN
2911       surf%facing(num_v) = 6 !IBSET( surf%facing(num_v), 1 )
2912       component          = 3
2913    ENDIF
[3833]2914
[4559]2915    IF ( east_facing )  THEN
2916       surf%facing(num_v) = 7 !IBSET( surf%facing(num_v), 2 )
2917       component          = 2
2918    ENDIF
[3833]2919
[4559]2920    IF ( west_facing )  THEN
2921       surf%facing(num_v) = 8 !IBSET( surf%facing(num_v), 3 )
2922       component          = 1
2923    ENDIF
[3833]2924
[4502]2925
[4559]2926    surf%z0(num_v)  = roughness_length
2927    surf%z0h(num_v) = z0h_factor * roughness_length
2928    surf%z0q(num_v) = z0h_factor * roughness_length
[3833]2929