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

Last change on this file since 3418 was 3418, checked in by kanani, 3 years ago

Add green facades, update building data base, fix for thin walls in spinup

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