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

Last change on this file since 3444 was 3444, checked in by schwenkel, 3 years ago

Bugfix in surface_wrd_local for cloud and rain water flux

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