source: palm/trunk/SOURCE/surface_layer_fluxes_mod.f90 @ 3597

Last change on this file since 3597 was 3597, checked in by maronga, 3 years ago

revised calculation of near surface air potential temperature

  • Property svn:keywords set to Id
File size: 89.8 KB
Line 
1!> @file surface_layer_fluxes_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2018 Leibniz Universitaet Hannover
18!
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! ------------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: surface_layer_fluxes_mod.f90 3597 2018-12-04 08:40:18Z maronga $
28! Added routine for calculating near surface air potential temperature (moved
29! from urban_surface_mod)
30!
31! 3547 2018-11-21 13:21:24Z suehring
32! variables documented
33!
34! 3361 2018-10-16 20:39:37Z knoop
35! Modularization of all bulk cloud physics code components
36!
37! 3271 2018-09-24 08:20:34Z suehring
38! Comment revised
39!
40! 3157 2018-07-19 21:08:49Z maronga
41! Added local free convection velocity scale w_lfc in calculation of uvw_abs.
42! This can be switche on/off by the user via the flag namelist parameter
43! use_free_convection_scaling
44!
45! 3152 2018-07-19 13:26:52Z suehring
46! q_surface is now part of surface structure
47!
48! 3149 2018-07-19 05:48:09Z maronga
49! Correct accidental last commit
50!
51! 3148 2018-07-19 05:45:25Z maronga
52! Major bugfix in calculation of Obukhov length
53!
54! 3146 2018-07-18 22:36:19Z maronga
55! move phi_m from turbulence_closure_mod
56!
57! 3045 2018-05-28 07:55:41Z Giersch
58! Error message revised
59!
60! 2766 2018-01-22 17:17:47Z kanani
61! Removed preprocessor directive __chem
62!
63! 2718 2018-01-02 08:49:38Z maronga
64! Corrected "Former revisions" section
65!
66! 2696 2017-12-14 17:12:51Z kanani
67! - Change in file header (GPL part)
68! - Implementation of chemistry module (FK)
69! - Added calculation of pt1 and qv1 for all surface types. Added calculation of
70!   pt_surface for default-type surfaces (BM)
71! - Add flag to disable computation of qsws in case of urban surface (MS)
72!
73! 2547 2017-10-16 12:41:56Z schwenkel
74! extended by cloud_droplets option
75!
76! 2321 2017-07-24 15:57:07Z schwenkel
77! Bugfix: Correct index in lookup table for Obukhov length
78!
79! 2299 2017-06-29 10:14:38Z suehring
80! Adjusted for allow separate spinups of LSM and atmosphere code
81!
82! 2292 2017-06-20 09:51:42Z schwenkel
83! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
84! includes two more prognostic equations for cloud drop concentration (nc) 
85! and cloud water content (qc).
86!
87! 2281 2017-06-13 11:34:50Z suehring
88! Clean-up unnecessary index access to surface type
89!
90! 2233 2017-05-30 18:08:54Z suehring
91!
92! 2232 2017-05-30 17:47:52Z suehring
93! Adjustments to new surface concept
94! OpenMP bugfix
95!
96! 2118 2017-01-17 16:38:49Z raasch
97! OpenACC directives and related code removed
98!
99! 2091 2016-12-21 16:38:18Z suehring
100! Bugfix in calculation of vsws ( incorrect linear interpolation of us )
101!
102! 2076 2016-12-02 13:54:20Z raasch
103! further openmp bugfix for lookup method
104!
105! 2073 2016-11-30 14:34:05Z raasch
106! openmp bugfix for lookup method
107!
108! 2037 2016-10-26 11:15:40Z knoop
109! Anelastic approximation implemented
110!
111! 2011 2016-09-19 17:29:57Z kanani
112! Flag urban_surface is now defined in module control_parameters.
113!
114! 2007 2016-08-24 15:47:17Z kanani
115! Account for urban surface model in computation of vertical kinematic heatflux
116!
117! 2000 2016-08-20 18:09:15Z knoop
118! Forced header and separation lines into 80 columns
119!
120! 1992 2016-08-12 15:14:59Z suehring
121! Minor bug, declaration of look-up index as INTEGER
122!
123! 1960 2016-07-12 16:34:24Z suehring
124! Treat humidity and passive scalar separately
125!
126! 1929 2016-06-09 16:25:25Z suehring
127! Bugfix: avoid segmentation fault in case one grid point is horizontally
128! completely surrounded by topography
129!
130! 1920 2016-05-30 10:50:15Z suehring
131! Avoid segmentation fault (see change in 1915) by different initialization of
132! us instead of adding a very small number in the denominator
133!
134! 1915 2016-05-27 11:05:02Z suehring
135! Bugfix: avoid segmentation fault in case of most_method = 'circular' at first
136! timestep
137!
138! 1850 2016-04-08 13:29:27Z maronga
139! Module renamed
140!
141!
142! 1822 2016-04-07 07:49:42Z hoffmann
143! icloud_scheme replaced by microphysics_*
144!
145! 1788 2016-03-10 11:01:04Z maronga
146! Added parameter z0q which replaces z0h in the similarity functions for
147! humidity.
148! Syntax layout improved.
149!
150! 1757 2016-02-22 15:49:32Z maronga
151! Minor fixes.
152!
153! 1749 2016-02-09 12:19:56Z raasch
154! further OpenACC adjustments
155!
156! 1747 2016-02-08 12:25:53Z raasch
157! adjustments for OpenACC usage
158!
159! 1709 2015-11-04 14:47:01Z maronga
160! Bugfix: division by zero could occur when calculating rib at low wind speeds
161! Bugfix: calculation of uv_total for neutral = .T., initial value for ol for
162! neutral = .T.
163!
164! 1705 2015-11-02 14:28:56Z maronga
165! Typo removed
166!
167! 1697 2015-10-28 17:14:10Z raasch
168! FORTRAN and OpenMP errors removed
169!
170! 1696 2015-10-27 10:03:34Z maronga
171! Modularized and completely re-written version of prandtl_fluxes.f90. In the
172! course of the re-writing two additional methods have been implemented. See
173! updated description.
174!
175! 1551 2015-03-03 14:18:16Z maronga
176! Removed land surface model part. The surface fluxes are now always calculated
177! within prandtl_fluxes, based on the given surface temperature/humidity (which
178! is either provided by the land surface model, by large scale forcing data, or
179! directly prescribed by the user.
180!
181! 1496 2014-12-02 17:25:50Z maronga
182! Adapted for land surface model
183!
184! 1494 2014-11-21 17:14:03Z maronga
185! Bugfixes: qs is now calculated before calculation of Rif. calculation of
186! buoyancy flux in Rif corrected (added missing humidity term), allow use of
187! topography for coupled runs (not tested)
188!
189! 1361 2014-04-16 15:17:48Z hoffmann
190! Bugfix: calculation of turbulent fluxes of rain water content (qrsws) and rain
191! drop concentration (nrsws) added
192!
193! 1340 2014-03-25 19:45:13Z kanani
194! REAL constants defined as wp-kind
195!
196! 1320 2014-03-20 08:40:49Z raasch
197! ONLY-attribute added to USE-statements,
198! kind-parameters added to all INTEGER and REAL declaration statements,
199! kinds are defined in new module kinds,
200! old module precision_kind is removed,
201! revision history before 2012 removed,
202! comment fields (!:) to be used for variable explanations added to
203! all variable declaration statements
204!
205! 1276 2014-01-15 13:40:41Z heinze
206! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
207!
208! 1257 2013-11-08 15:18:40Z raasch
209! openACC "kernels do" replaced by "kernels loop", "loop independent" added
210!
211! 1036 2012-10-22 13:43:42Z raasch
212! code put under GPL (PALM 3.9)
213!
214! 1015 2012-09-27 09:23:24Z raasch
215! OpenACC statements added
216!
217! 978 2012-08-09 08:28:32Z fricke
218! roughness length for scalar quantities z0h added
219!
220! Revision 1.1  1998/01/23 10:06:06  raasch
221! Initial revision
222!
223!
224! Description:
225! ------------
226!> Diagnostic computation of vertical fluxes in the constant flux layer from the
227!> values of the variables at grid point k=1. Three different methods are
228!> available:
229!> 1) the "old" version (most_method = 'circular') which is fast, but inaccurate
230!> 2) a Newton iteration method (most_method = 'newton'), which is accurate, but
231!>    slower
232!> 3) a method using a lookup table which is fast and accurate. Note, however,
233!>    that this method cannot be used in case of roughness heterogeneity
234!>
235!> @todo (re)move large_scale_forcing actions
236!> @todo check/optimize OpenMP directives
237!> @todo simplify if conditions (which flux need to be computed in which case)
238!------------------------------------------------------------------------------!
239 MODULE surface_layer_fluxes_mod
240
241    USE arrays_3d,                                                             &
242        ONLY:  e, kh, nc, nr, pt, q, ql, qc, qr, s, u, v, vpt, w, zu, zw,      &
243               drho_air_zw, rho_air_zw, d_exner
244
245    USE basic_constants_and_equations_mod,                                     &
246        ONLY:  g, kappa, lv_d_cp, pi, rd_d_rv
247
248    USE chem_modules,                                                          &
249        ONLY:  constant_csflux, nvar
250
251    USE cpulog
252
253    USE control_parameters,                                                    &
254        ONLY:  air_chemistry, cloud_droplets,                                  &
255               constant_heatflux, constant_scalarflux,                         &
256               constant_waterflux, coupling_mode, do_output_at_2m, humidity,   &
257               ibc_e_b, ibc_pt_b, indoor_model, initializing_actions,          &
258               intermediate_timestep_count, intermediate_timestep_count_max,   &
259               land_surface, large_scale_forcing, lsf_surf, message_string,    &
260               most_method, neutral, passive_scalar, pt_surface, q_surface,    &
261               run_coupled, surface_pressure, simulated_time, terminate_run,   &
262               time_since_reference_point, urban_surface,                      &
263               use_free_convection_scaling, zeta_max, zeta_min
264
265    USE grid_variables,                                                        &
266        ONLY:  dx, dy 
267
268    USE indices,                                                               &
269        ONLY:  nxl, nxr, nys, nyn, nzb
270
271    USE kinds
272
273    USE bulk_cloud_model_mod,                                                  &
274        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
275
276    USE pegrid
277
278    USE land_surface_model_mod,                                                &
279        ONLY:  aero_resist_kray, skip_time_do_lsm
280
281    USE surface_mod,                                                           &
282        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_type,     &
283                surf_usm_h, surf_usm_v
284       
285
286    IMPLICIT NONE
287
288    INTEGER(iwp) ::  i              !< loop index x direction
289    INTEGER(iwp) ::  j              !< loop index y direction
290    INTEGER(iwp) ::  k              !< loop index z direction
291    INTEGER(iwp) ::  l              !< loop index for surf type
292    INTEGER(iwp) ::  li_bnd  = 7500 !< Lookup table index of the last time step
293
294    INTEGER(iwp), PARAMETER ::  num_steps = 15000  !< number of steps in the lookup table
295
296    LOGICAL      ::  coupled_run       !< Flag for coupled atmosphere-ocean runs
297    LOGICAL      ::  downward = .FALSE.!< Flag indicating downward-facing horizontal surface
298    LOGICAL      ::  mom_uv  = .FALSE. !< Flag indicating calculation of usvs and vsus at vertical surfaces
299    LOGICAL      ::  mom_w   = .FALSE. !< Flag indicating calculation of wsus and wsvs at vertical surfaces
300    LOGICAL      ::  mom_tke = .FALSE. !< Flag indicating calculation of momentum fluxes at vertical surfaces used for TKE production
301    LOGICAL      ::  surf_vertical     !< Flag indicating vertical surfaces
302
303    REAL(wp), DIMENSION(0:num_steps-1) :: rib_tab,  & !< Lookup table bulk Richardson number
304                                          ol_tab      !< Lookup table values of L
305
306    REAL(wp)     ::  e_s,               & !< Saturation water vapor pressure
307                     ol_max = 1.0E6_wp, & !< Maximum Obukhov length
308                     rib_max,           & !< Maximum Richardson number in lookup table
309                     rib_min,           & !< Minimum Richardson number in lookup table
310                     z_mo                 !< Height of the constant flux layer where MOST is assumed
311
312    TYPE(surf_type), POINTER ::  surf     !< surf-type array, used to generalize subroutines
313
314
315    SAVE
316
317    PRIVATE
318
319    PUBLIC init_surface_layer_fluxes, phi_m, surface_layer_fluxes
320
321    INTERFACE init_surface_layer_fluxes
322       MODULE PROCEDURE init_surface_layer_fluxes
323    END INTERFACE init_surface_layer_fluxes
324
325    INTERFACE phi_m
326       MODULE PROCEDURE phi_m
327    END INTERFACE phi_m
328
329    INTERFACE surface_layer_fluxes
330       MODULE PROCEDURE surface_layer_fluxes
331    END INTERFACE surface_layer_fluxes
332
333
334 CONTAINS
335
336
337!------------------------------------------------------------------------------!
338! Description:
339! ------------
340!> Main routine to compute the surface fluxes
341!------------------------------------------------------------------------------!
342    SUBROUTINE surface_layer_fluxes
343
344       IMPLICIT NONE
345
346       surf_vertical = .FALSE. !< flag indicating vertically orientated surface elements
347       downward      = .FALSE. !< flag indicating downward-facing surface elements
348!
349!--    Derive potential temperature and specific humidity at first grid level
350!--    from the fields pt and q
351!
352!--    First call for horizontal default-type surfaces (l=0 - upward facing,
353!--    l=1 - downward facing)
354       DO  l = 0, 1
355          IF ( surf_def_h(l)%ns >= 1  )  THEN
356             surf => surf_def_h(l)
357             CALL calc_pt_q
358             IF ( .NOT. neutral )  THEN
359                CALL calc_pt_surface
360                IF ( humidity )  THEN
361                   CALL calc_q_surface
362                   CALL calc_vpt_surface
363                ENDIF
364             ENDIF
365          ENDIF
366       ENDDO
367!
368!--    Call for natural-type horizontal surfaces
369       IF ( surf_lsm_h%ns >= 1  )  THEN
370          surf => surf_lsm_h
371          CALL calc_pt_q
372       ENDIF
373
374!
375!--    Call for urban-type horizontal surfaces
376       IF ( surf_usm_h%ns >= 1  )  THEN
377          surf => surf_usm_h
378          CALL calc_pt_q
379       ENDIF
380
381!
382!--    Call for natural-type vertical surfaces
383       DO  l = 0, 3
384          IF ( surf_lsm_v(l)%ns >= 1  )  THEN
385             surf => surf_lsm_v(l)
386             CALL calc_pt_q
387          ENDIF
388
389!--       Call for urban-type vertical surfaces
390          IF ( surf_usm_v(l)%ns >= 1  )  THEN
391             surf => surf_usm_v(l)
392             CALL calc_pt_q
393          ENDIF
394       ENDDO
395
396!
397!--    First, calculate the new Obukhov length, then new friction velocity,
398!--    followed by the new scaling parameters (th*, q*, etc.), and the new
399!--    surface fluxes if required. The old routine ("circular") requires a
400!--    different order of calls as the scaling parameters from the previous time
401!--    steps are used to calculate the Obukhov length
402
403!
404!--    Depending on setting of most_method use the "old" routine
405!--    Note, each routine is called for different surface types.
406!--    First call for default-type horizontal surfaces, for natural- and
407!--    urban-type surfaces. Note, at this place only upward-facing horizontal
408!--    surfaces are treted.
409       IF ( most_method == 'circular' )  THEN
410!
411!--       Default-type upward-facing horizontal surfaces
412          IF ( surf_def_h(0)%ns >= 1  )  THEN
413             surf => surf_def_h(0)
414             CALL calc_scaling_parameters
415             CALL calc_uvw_abs
416             IF ( .NOT. neutral )  CALL calc_ol
417             CALL calc_us
418             CALL calc_surface_fluxes
419             
420             IF ( do_output_at_2m )  THEN
421                CALL calc_pt_near_surface ( '2m' )
422             ENDIF
423          ENDIF
424!
425!--       Natural-type horizontal surfaces
426          IF ( surf_lsm_h%ns >= 1  )  THEN
427             surf => surf_lsm_h
428             CALL calc_scaling_parameters
429             CALL calc_uvw_abs
430             IF ( .NOT. neutral )  CALL calc_ol
431             CALL calc_us
432             CALL calc_surface_fluxes
433             IF ( do_output_at_2m )  THEN
434                CALL calc_pt_near_surface ( '2m' )
435             ENDIF
436          ENDIF
437!
438!--       Urban-type horizontal surfaces
439          IF ( surf_usm_h%ns >= 1  )  THEN
440             surf => surf_usm_h
441             CALL calc_scaling_parameters
442             CALL calc_uvw_abs
443             IF ( .NOT. neutral )  CALL calc_ol
444             CALL calc_us
445             CALL calc_surface_fluxes
446             IF ( do_output_at_2m )  THEN
447                CALL calc_pt_near_surface ( '2m' )
448             ENDIF
449             IF ( indoor_model )  THEN
450                CALL calc_pt_near_surface ( '10cm' )
451             ENDIF
452          ENDIF
453!
454!--    Use either Newton iteration or a lookup table for the bulk Richardson
455!--    number to calculate the Obukhov length
456       ELSEIF ( most_method == 'newton'  .OR.  most_method == 'lookup' )  THEN
457!
458!--       Default-type upward-facing horizontal surfaces
459          IF ( surf_def_h(0)%ns >= 1  )  THEN
460             surf => surf_def_h(0)
461             CALL calc_uvw_abs
462             IF ( .NOT. neutral )  CALL calc_ol
463             CALL calc_us
464             CALL calc_scaling_parameters
465             CALL calc_surface_fluxes
466             IF ( do_output_at_2m )  THEN
467                CALL calc_pt_near_surface ( '2m' )
468             ENDIF
469          ENDIF
470!
471!--       Natural-type horizontal surfaces
472          IF ( surf_lsm_h%ns >= 1  )  THEN
473             surf => surf_lsm_h
474             CALL calc_uvw_abs
475             IF ( .NOT. neutral )  CALL calc_ol
476             CALL calc_us
477             CALL calc_scaling_parameters
478             CALL calc_surface_fluxes
479             IF ( do_output_at_2m )  THEN
480                CALL calc_pt_near_surface ( '2m' )
481             ENDIF
482          ENDIF
483!
484!--       Urban-type horizontal surfaces
485          IF ( surf_usm_h%ns >= 1  )  THEN
486             surf => surf_usm_h
487             CALL calc_uvw_abs
488             IF ( .NOT. neutral )  CALL calc_ol
489             CALL calc_us
490             CALL calc_scaling_parameters
491             CALL calc_surface_fluxes
492             IF ( do_output_at_2m )  THEN
493                CALL calc_pt_near_surface ( '2m' )
494             ENDIF
495             IF ( indoor_model )  THEN
496                CALL calc_pt_near_surface ( '10cm' )
497             ENDIF
498          ENDIF
499
500       ENDIF
501!
502!--    Treat downward-facing horizontal surfaces. Note, so far, these are
503!--    always default type. Stratification is not considered
504!--    in this case, hence, no further distinction between different
505!--    most_method is required. 
506       IF ( surf_def_h(1)%ns >= 1  )  THEN
507          downward = .TRUE.
508          surf => surf_def_h(1)
509          CALL calc_uvw_abs
510          CALL calc_us
511          CALL calc_surface_fluxes
512          downward = .FALSE.
513       ENDIF
514!
515!--    Calculate surfaces fluxes at vertical surfaces for momentum
516!--    and subgrid-scale TKE.
517!--    No stability is considered. Therefore, scaling parameters and Obukhov-
518!--    length do not need to be calculated and no distinction in 'circular',
519!--    'Newton' or 'lookup' is necessary so far.
520!--    Note, this will change if stability is once considered.
521       surf_vertical = .TRUE.
522!
523!--    Calculate horizontal momentum fluxes at north- and south-facing
524!--    surfaces(usvs).
525!--    For default-type surfaces
526       mom_uv = .TRUE. 
527       DO  l = 0, 1
528          IF ( surf_def_v(l)%ns >= 1  )  THEN
529             surf => surf_def_v(l)
530!
531!--          Compute surface-parallel velocity
532             CALL calc_uvw_abs_v_ugrid
533!
534!--          Compute respective friction velocity on staggered grid
535             CALL calc_us
536!
537!--          Compute respective surface fluxes for momentum and TKE
538             CALL calc_surface_fluxes
539          ENDIF
540       ENDDO
541!
542!--    For natural-type surfaces. Please note, even though stability is not
543!--    considered for the calculation of momentum fluxes at vertical surfaces,
544!--    scaling parameters and Obukov length are calculated nevertheless in this
545!--    case. This is due to the requirement of ts in parameterization of heat
546!--    flux in land-surface model in case of aero_resist_kray is not true.
547       IF ( .NOT. aero_resist_kray )  THEN
548          IF ( most_method == 'circular' )  THEN
549             DO  l = 0, 1
550                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
551                   surf => surf_lsm_v(l)
552!
553!--                Compute scaling parameters
554                   CALL calc_scaling_parameters
555!
556!--                Compute surface-parallel velocity
557                   CALL calc_uvw_abs_v_ugrid
558!
559!--                Compute Obukhov length
560                   IF ( .NOT. neutral )  CALL calc_ol
561!
562!--                Compute respective friction velocity on staggered grid
563                   CALL calc_us
564!
565!--                Compute respective surface fluxes for momentum and TKE
566                   CALL calc_surface_fluxes
567                ENDIF
568             ENDDO
569          ELSE
570             DO  l = 0, 1
571                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
572                   surf => surf_lsm_v(l)
573!
574!--                Compute surface-parallel velocity
575                   CALL calc_uvw_abs_v_ugrid
576!
577!--                Compute Obukhov length
578                   IF ( .NOT. neutral )  CALL calc_ol
579!
580!--                Compute respective friction velocity on staggered grid
581                   CALL calc_us
582!
583!--                Compute scaling parameters
584                   CALL calc_scaling_parameters
585!
586!--                Compute respective surface fluxes for momentum and TKE
587                   CALL calc_surface_fluxes
588                ENDIF
589             ENDDO
590          ENDIF
591!
592!--    No ts is required, so scaling parameters and Obukhov length do not need
593!--    to be computed.
594       ELSE
595          DO  l = 0, 1
596             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
597                surf => surf_lsm_v(l)
598!   
599!--             Compute surface-parallel velocity
600                CALL calc_uvw_abs_v_ugrid
601!
602!--             Compute respective friction velocity on staggered grid
603                CALL calc_us
604!
605!--             Compute respective surface fluxes for momentum and TKE
606                CALL calc_surface_fluxes
607             ENDIF
608          ENDDO
609       ENDIF
610!
611!--    For urban-type surfaces
612       DO  l = 0, 1
613          IF ( surf_usm_v(l)%ns >= 1  )  THEN
614             surf => surf_usm_v(l)
615!
616!--          Compute surface-parallel velocity
617             CALL calc_uvw_abs_v_ugrid
618!
619!--          Compute respective friction velocity on staggered grid
620             CALL calc_us
621!
622!--          Compute respective surface fluxes for momentum and TKE
623             CALL calc_surface_fluxes
624          ENDIF
625       ENDDO
626!
627!--    Calculate horizontal momentum fluxes at east- and west-facing
628!--    surfaces (vsus).
629!--    For default-type surfaces
630       DO  l = 2, 3
631          IF ( surf_def_v(l)%ns >= 1  )  THEN
632             surf => surf_def_v(l)
633!
634!--          Compute surface-parallel velocity
635             CALL calc_uvw_abs_v_vgrid
636!
637!--          Compute respective friction velocity on staggered grid
638             CALL calc_us
639!
640!--          Compute respective surface fluxes for momentum and TKE
641             CALL calc_surface_fluxes
642          ENDIF
643       ENDDO
644!
645!--    For natural-type surfaces. Please note, even though stability is not
646!--    considered for the calculation of momentum fluxes at vertical surfaces,
647!--    scaling parameters and Obukov length are calculated nevertheless in this
648!--    case. This is due to the requirement of ts in parameterization of heat
649!--    flux in land-surface model in case of aero_resist_kray is not true.
650       IF ( .NOT. aero_resist_kray )  THEN
651          IF ( most_method == 'circular' )  THEN
652             DO  l = 2, 3
653                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
654                   surf => surf_lsm_v(l)
655!
656!--                Compute scaling parameters
657                   CALL calc_scaling_parameters
658!
659!--                Compute surface-parallel velocity
660                   CALL calc_uvw_abs_v_vgrid
661!
662!--                Compute Obukhov length
663                   IF ( .NOT. neutral )  CALL calc_ol
664!
665!--                Compute respective friction velocity on staggered grid
666                   CALL calc_us
667!
668!--                Compute respective surface fluxes for momentum and TKE
669                   CALL calc_surface_fluxes
670                ENDIF
671             ENDDO
672          ELSE
673             DO  l = 2, 3
674                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
675                   surf => surf_lsm_v(l)
676!
677!--                Compute surface-parallel velocity
678                   CALL calc_uvw_abs_v_vgrid
679!
680!--                Compute Obukhov length
681                   IF ( .NOT. neutral )  CALL calc_ol
682!
683!--                Compute respective friction velocity on staggered grid
684                   CALL calc_us
685!
686!--                Compute scaling parameters
687                   CALL calc_scaling_parameters
688!
689!--                Compute respective surface fluxes for momentum and TKE
690                   CALL calc_surface_fluxes
691                ENDIF
692             ENDDO
693          ENDIF
694       ELSE
695          DO  l = 2, 3
696             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
697                surf => surf_lsm_v(l)
698!   
699!--             Compute surface-parallel velocity
700                CALL calc_uvw_abs_v_vgrid
701!
702!--             Compute respective friction velocity on staggered grid
703                CALL calc_us
704!
705!--             Compute respective surface fluxes for momentum and TKE
706                CALL calc_surface_fluxes
707             ENDIF
708          ENDDO
709       ENDIF
710!
711!--    For urban-type surfaces
712       DO  l = 2, 3
713          IF ( surf_usm_v(l)%ns >= 1  )  THEN
714             surf => surf_usm_v(l)
715!
716!--          Compute surface-parallel velocity
717             CALL calc_uvw_abs_v_vgrid
718!
719!--          Compute respective friction velocity on staggered grid
720             CALL calc_us
721!
722!--          Compute respective surface fluxes for momentum and TKE
723             CALL calc_surface_fluxes
724          ENDIF
725       ENDDO
726       mom_uv = .FALSE.
727!
728!--    Calculate horizontal momentum fluxes of w (wsus and wsvs) at vertial
729!--    surfaces.
730       mom_w = .TRUE.
731!
732!--    Default-type surfaces
733       DO  l = 0, 3
734          IF ( surf_def_v(l)%ns >= 1  )  THEN
735             surf => surf_def_v(l)
736             CALL calc_uvw_abs_v_wgrid
737             CALL calc_us
738             CALL calc_surface_fluxes
739          ENDIF
740       ENDDO
741!
742!--    Natural-type surfaces
743       DO  l = 0, 3
744          IF ( surf_lsm_v(l)%ns >= 1  )  THEN
745             surf => surf_lsm_v(l)
746             CALL calc_uvw_abs_v_wgrid
747             CALL calc_us
748             CALL calc_surface_fluxes
749          ENDIF
750       ENDDO
751!
752!--    Urban-type surfaces
753       DO  l = 0, 3
754          IF ( surf_usm_v(l)%ns >= 1  )  THEN
755             surf => surf_usm_v(l)
756             CALL calc_uvw_abs_v_wgrid
757             CALL calc_us
758             CALL calc_surface_fluxes
759          ENDIF
760       ENDDO
761       mom_w = .FALSE.   
762!
763!--    Calculate momentum fluxes usvs, vsus, wsus and wsvs at vertical
764!--    surfaces for TKE production. Note, here, momentum fluxes are defined
765!--    at grid center and are not staggered as before.
766       mom_tke = .TRUE.
767!
768!--    Default-type surfaces
769       DO  l = 0, 3
770          IF ( surf_def_v(l)%ns >= 1  )  THEN
771             surf => surf_def_v(l)
772             CALL calc_uvw_abs_v_sgrid
773             CALL calc_us
774             CALL calc_surface_fluxes
775          ENDIF
776       ENDDO
777!
778!--    Natural-type surfaces
779       DO  l = 0, 3
780          IF ( surf_lsm_v(l)%ns >= 1  )  THEN
781             surf => surf_lsm_v(l)
782             CALL calc_uvw_abs_v_sgrid
783             CALL calc_us
784             CALL calc_surface_fluxes
785          ENDIF
786       ENDDO
787!
788!--    Urban-type surfaces
789       DO  l = 0, 3
790          IF ( surf_usm_v(l)%ns >= 1  )  THEN
791             surf => surf_usm_v(l)
792             CALL calc_uvw_abs_v_sgrid
793             CALL calc_us
794             CALL calc_surface_fluxes
795          ENDIF
796       ENDDO
797       mom_tke = .FALSE.
798 
799
800    END SUBROUTINE surface_layer_fluxes
801
802
803!------------------------------------------------------------------------------!
804! Description:
805! ------------
806!> Initializing actions for the surface layer routine. Basically, this involves
807!> the preparation of a lookup table for the the bulk Richardson number vs
808!> Obukhov length L when using the lookup table method.
809!------------------------------------------------------------------------------!
810    SUBROUTINE init_surface_layer_fluxes
811
812       IMPLICIT NONE
813
814       INTEGER(iwp) :: li,         & !< Index for loop to create lookup table
815                       num_steps_n   !< Number of non-stretched zeta steps
816
817       LOGICAL :: terminate_run_l = .FALSE.    !< Flag to terminate run (global)
818
819       REAL(wp), PARAMETER ::  zeta_stretch = -10.0_wp !< Start of stretching in the free convection limit
820                               
821       REAL(wp), DIMENSION(:), ALLOCATABLE :: zeta_tmp
822
823
824       REAL(wp) :: zeta_step,            & !< Increment of zeta
825                   regr      = 1.01_wp,  & !< Stretching factor of zeta_step in the free convection limit
826                   regr_old  = 1.0E9_wp, & !< Stretching factor of last iteration step
827                   z0h_min   = 0.0_wp,   & !< Minimum value of z0h to create table
828                   z0_min    = 0.0_wp      !< Minimum value of z0 to create table
829
830
831
832
833!
834!--    In case of runs with neutral statification, set Obukhov length to a
835!--    large value
836       IF ( neutral )  THEN
837          IF ( surf_def_h(0)%ns >= 1 )  surf_def_h(0)%ol = 1.0E10_wp
838          IF ( surf_lsm_h%ns    >= 1 )  surf_lsm_h%ol    = 1.0E10_wp
839          IF ( surf_usm_h%ns    >= 1 )  surf_usm_h%ol    = 1.0E10_wp
840       ENDIF
841
842       IF ( most_method == 'lookup' )  THEN
843
844!
845!--       Check for roughness heterogeneity. In that case terminate run and
846!--       inform user. Check for both, natural and non-natural walls.
847          IF ( surf_def_h(0)%ns >= 1 )  THEN
848             IF ( MINVAL( surf_def_h(0)%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR. &
849                  MINVAL( surf_def_h(0)%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
850                terminate_run_l = .TRUE.
851             ENDIF
852          ENDIF
853          IF ( surf_lsm_h%ns >= 1 )  THEN
854             IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h )  .OR.       &
855                  MINVAL( surf_lsm_h%z0  ) /= MAXVAL( surf_lsm_h%z0  ) )  THEN
856                terminate_run_l = .TRUE.
857             ENDIF
858          ENDIF
859          IF ( surf_usm_h%ns >= 1 )  THEN
860             IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_usm_h%z0h )  .OR.       &
861                  MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_usm_h%z0  ) )  THEN
862                terminate_run_l = .TRUE.
863             ENDIF
864          ENDIF
865!
866!--       Check roughness homogeneity between differt surface types.
867          IF ( surf_lsm_h%ns >= 1  .AND.  surf_def_h(0)%ns >= 1 )  THEN
868             IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR.    &
869                  MINVAL( surf_lsm_h%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
870                terminate_run_l = .TRUE.
871             ENDIF
872          ENDIF
873          IF ( surf_usm_h%ns >= 1  .AND.  surf_def_h(0)%ns >= 1 )  THEN
874             IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR.    &
875                  MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
876                terminate_run_l = .TRUE.
877             ENDIF
878          ENDIF
879          IF ( surf_usm_h%ns >= 1  .AND.  surf_lsm_h%ns >= 1 )  THEN
880             IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h )  .OR.       &
881                  MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_lsm_h%z0  ) )  THEN
882                terminate_run_l = .TRUE.
883             ENDIF
884          ENDIF
885
886
887#if defined( __parallel )
888!
889!--       Make a logical OR for all processes. Force termiation of model if result
890!--       is TRUE
891          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
892          CALL MPI_ALLREDUCE( terminate_run_l, terminate_run, 1, MPI_LOGICAL,  &
893                              MPI_LOR, comm2d, ierr )
894#else
895          terminate_run = terminate_run_l
896#endif
897
898          IF ( terminate_run )  THEN
899             message_string = 'most_method = "lookup" cannot be used in ' //   &
900                              'combination with a prescribed roughness '  //   &
901                              'heterogeneity'
902             CALL message( 'surface_layer_fluxes', 'PA0116', 1, 2, 0, 6, 0 )
903          ENDIF
904
905          ALLOCATE(  zeta_tmp(0:num_steps-1) )
906
907!
908!--       Use the lowest possible value for z_mo
909          k    = nzb
910          z_mo = zu(k+1) - zw(k)
911
912!
913!--       Calculate z/L range from zeta_stretch to zeta_max using 90% of the
914!--       available steps (num_steps). The calculation is done with negative
915!--       values of zeta in order to simplify the stretching in the free
916!--       convection limit for the remaining 10% of steps.
917          zeta_tmp(0) = - zeta_max
918          num_steps_n = ( num_steps * 9 / 10 ) - 1
919          zeta_step   = (zeta_max - zeta_stretch) / REAL(num_steps_n)
920
921          DO li = 1, num_steps_n
922             zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
923          ENDDO
924
925!
926!--       Calculate stretching factor for the free convection range
927          DO  WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
928             regr_old = regr
929             regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr )   &
930                    )**( 10.0_wp / REAL(num_steps) )
931          ENDDO
932
933!
934!--       Calculate z/L range from zeta_min to zeta_stretch
935          DO li = num_steps_n+1, num_steps-1
936             zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
937             zeta_step = zeta_step * regr
938          ENDDO
939
940!
941!--       Save roughness lengths to temporary variables
942          IF ( surf_def_h(0)%ns >= 1  )  THEN
943             z0h_min = surf_def_h(0)%z0h(1)
944             z0_min  = surf_def_h(0)%z0(1)
945          ELSEIF ( surf_lsm_h%ns >= 1 )  THEN
946             z0h_min = surf_lsm_h%z0h(1)
947             z0_min  = surf_lsm_h%z0(1)
948          ELSEIF ( surf_usm_h%ns >= 1 )  THEN
949             z0h_min = surf_usm_h%z0h(1)
950             z0_min  = surf_usm_h%z0(1)
951          ENDIF         
952!
953!--       Calculate lookup table for the Richardson number versus Obukhov length
954!--       The Richardson number (rib) is defined depending on the choice of
955!--       boundary conditions for temperature
956          IF ( ibc_pt_b == 1 )  THEN
957             DO li = 0, num_steps-1
958                ol_tab(li)  = - z_mo / zeta_tmp(num_steps-1-li)
959                rib_tab(li) = z_mo / ol_tab(li)  / ( LOG( z_mo / z0_min )       &
960                                                - psi_m( z_mo / ol_tab(li) )    &
961                                                + psi_m( z0_min / ol_tab(li) )  &
962                                                  )**3
963             ENDDO 
964          ELSE
965             DO li = 0, num_steps-1
966                ol_tab(li)  = - z_mo / zeta_tmp(num_steps-1-li)
967                rib_tab(li) = z_mo / ol_tab(li)  * ( LOG( z_mo / z0h_min )     &
968                                              - psi_h( z_mo / ol_tab(li) )     &
969                                              + psi_h( z0h_min / ol_tab(li) )  &
970                                            )                                  &
971                                          / ( LOG( z_mo / z0_min )             &
972                                              - psi_m( z_mo / ol_tab(li) )     &
973                                              + psi_m( z0_min / ol_tab(li) )   &
974                                            )**2
975             ENDDO
976          ENDIF
977
978!
979!--       Determine minimum values of rib in the lookup table. Set upper limit
980!--       to critical Richardson number (0.25)
981          rib_min  = MINVAL(rib_tab)
982          rib_max  = 0.25 !MAXVAL(rib_tab)
983
984          DEALLOCATE( zeta_tmp )
985       ENDIF
986
987    END SUBROUTINE init_surface_layer_fluxes
988
989
990!------------------------------------------------------------------------------!
991! Description:
992! ------------
993!> Compute the absolute value of the horizontal velocity (relative to the
994!> surface) for horizontal surface elements. This is required by all methods.
995!------------------------------------------------------------------------------!
996    SUBROUTINE calc_uvw_abs
997   
998       IMPLICIT NONE
999
1000       INTEGER(iwp) ::  i             !< running index x direction
1001       INTEGER(iwp) ::  ibit          !< flag to mask computation of relative velocity in case of downward-facing surfaces
1002       INTEGER(iwp) ::  j             !< running index y direction
1003       INTEGER(iwp) ::  k             !< running index z direction
1004       INTEGER(iwp) ::  m             !< running index surface elements
1005
1006       REAL(wp)     :: w_lfc          !< local free convection velocity scale
1007!
1008!--    ibit is 1 for upward-facing surfaces, zero for downward-facing surfaces.
1009       ibit = MERGE( 1, 0, .NOT. downward )
1010
1011       DO  m = 1, surf%ns
1012
1013          i   = surf%i(m)           
1014          j   = surf%j(m)
1015          k   = surf%k(m)
1016       
1017!
1018!--       Calculate free convection velocity scale w_lfc is
1019!--       use_free_convection_scaling = .T.. This will maintain a horizontal
1020!--       velocity even for very weak wind convective conditions. SIGN is used
1021!--       to set w_lfc to zero under stable conditions.
1022          IF ( use_free_convection_scaling )  THEN
1023             w_lfc = ABS(g / surf%pt1(m) * surf%z_mo(m) * surf%shf(m))
1024             w_lfc = ( 0.5_wp * ( w_lfc + SIGN(w_lfc,surf%shf(m)) ) )**(0.33333_wp)
1025          ELSE
1026             w_lfc = 0.0_wp
1027          ENDIF
1028
1029!
1030!--       Compute the absolute value of the horizontal velocity.
1031!--       (relative to the surface in case the lower surface is the ocean).
1032!--       Please note, in new surface modelling concept the index values changed,
1033!--       i.e. the reference grid point is not the surface-grid point itself but
1034!--       the first grid point outside of the topography.
1035!--       Note, in case of coupled ocean-atmosphere simulations relative velocity
1036!--       with respect to the ocean surface is used, hence, (k-1,j,i) values
1037!--       are used to calculate the absolute velocity.
1038!--       However, this do not apply for downward-facing walls. To mask this,
1039!--       use ibit, which checks for upward/downward-facing surfaces.
1040          surf%uvw_abs(m) = SQRT(                                              &
1041                              ( 0.5_wp * (   u(k,j,i)   + u(k,j,i+1)           &
1042                                        -  ( u(k-1,j,i) + u(k-1,j,i+1)         &
1043                                           ) * ibit                            &
1044                                         )                                     &
1045                              )**2 +                                           &
1046                              ( 0.5_wp * (   v(k,j,i)   + v(k,j+1,i)           &
1047                                        -  ( v(k-1,j,i) + v(k-1,j+1,i)         &
1048                                           ) * ibit                            &
1049                                         )                                     &
1050                              )**2  + w_lfc**2                                 &
1051                                )
1052
1053                               
1054
1055       ENDDO
1056
1057    END SUBROUTINE calc_uvw_abs
1058
1059
1060!------------------------------------------------------------------------------!
1061! Description:
1062! ------------
1063!> Compute the absolute value of the horizontal velocity (relative to the
1064!> surface) for horizontal surface elements. This is required by all methods.
1065!------------------------------------------------------------------------------!
1066    SUBROUTINE calc_uvw_abs_v_ugrid
1067
1068       IMPLICIT NONE
1069
1070       INTEGER(iwp) ::  i   !< running index x direction
1071       INTEGER(iwp) ::  j   !< running index y direction
1072       INTEGER(iwp) ::  k   !< running index z direction
1073       INTEGER(iwp) ::  m   !< running index surface elements
1074
1075       REAL(wp)     ::  u_i !< u-component on xu-grid
1076       REAL(wp)     ::  w_i !< w-component on xu-grid
1077
1078
1079       DO  m = 1, surf%ns
1080          i   = surf%i(m)           
1081          j   = surf%j(m)
1082          k   = surf%k(m)
1083!
1084!--       Compute the absolute value of the surface parallel velocity on u-grid.
1085          u_i  = u(k,j,i)
1086          w_i  = 0.25_wp * ( w(k-1,j,i-1) + w(k-1,j,i) +                       &
1087                             w(k,j,i-1)   + w(k,j,i) ) 
1088
1089          surf%uvw_abs(m) = SQRT( u_i**2 + w_i**2 ) 
1090
1091       ENDDO
1092
1093    END SUBROUTINE calc_uvw_abs_v_ugrid
1094
1095!------------------------------------------------------------------------------!
1096! Description:
1097! ------------
1098!> Compute the absolute value of the horizontal velocity (relative to the
1099!> surface) for horizontal surface elements. This is required by all methods.
1100!------------------------------------------------------------------------------!
1101    SUBROUTINE calc_uvw_abs_v_vgrid
1102
1103       IMPLICIT NONE
1104
1105       INTEGER(iwp) ::  i   !< running index x direction
1106       INTEGER(iwp) ::  j   !< running index y direction
1107       INTEGER(iwp) ::  k   !< running index z direction
1108       INTEGER(iwp) ::  m   !< running index surface elements
1109
1110       REAL(wp)     ::  v_i !< v-component on yv-grid
1111       REAL(wp)     ::  w_i !< w-component on yv-grid
1112
1113
1114       DO  m = 1, surf%ns
1115          i   = surf%i(m)           
1116          j   = surf%j(m)
1117          k   = surf%k(m)
1118
1119          v_i  = u(k,j,i)
1120          w_i  = 0.25_wp * ( w(k-1,j-1,i) + w(k-1,j,i) +                       &
1121                             w(k,j-1,i)   + w(k,j,i) ) 
1122
1123          surf%uvw_abs(m) = SQRT( v_i**2 + w_i**2 ) 
1124
1125       ENDDO
1126
1127    END SUBROUTINE calc_uvw_abs_v_vgrid
1128
1129!------------------------------------------------------------------------------!
1130! Description:
1131! ------------
1132!> Compute the absolute value of the horizontal velocity (relative to the
1133!> surface) for horizontal surface elements. This is required by all methods.
1134!------------------------------------------------------------------------------!
1135    SUBROUTINE calc_uvw_abs_v_wgrid
1136
1137       IMPLICIT NONE
1138
1139       INTEGER(iwp) ::  i   !< running index x direction
1140       INTEGER(iwp) ::  j   !< running index y direction
1141       INTEGER(iwp) ::  k   !< running index z direction
1142       INTEGER(iwp) ::  m   !< running index surface elements
1143
1144       REAL(wp)     ::  u_i !< u-component on x-zw-grid
1145       REAL(wp)     ::  v_i !< v-component on y-zw-grid
1146       REAL(wp)     ::  w_i !< w-component on zw-grid
1147!
1148!--    North- (l=0) and south-facing (l=1) surfaces
1149       IF ( l == 0  .OR.  l == 1 )  THEN
1150          DO  m = 1, surf%ns
1151             i   = surf%i(m)           
1152             j   = surf%j(m)
1153             k   = surf%k(m)
1154
1155             u_i  = 0.25_wp * ( u(k+1,j,i+1) + u(k+1,j,i) +                    &
1156                                u(k,j,i+1)   + u(k,j,i) )
1157             v_i  = 0.0_wp
1158             w_i  = w(k,j,i)
1159
1160             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
1161          ENDDO
1162!
1163!--    East- (l=2) and west-facing (l=3) surfaces
1164       ELSE
1165          DO  m = 1, surf%ns
1166             i   = surf%i(m)           
1167             j   = surf%j(m)
1168             k   = surf%k(m)
1169
1170             u_i  = 0.0_wp
1171             v_i  = 0.25_wp * ( v(k+1,j+1,i) + v(k+1,j,i) +                    &
1172                                v(k,j+1,i)   + v(k,j,i) )
1173             w_i  = w(k,j,i)
1174
1175             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
1176          ENDDO
1177       ENDIF           
1178
1179    END SUBROUTINE calc_uvw_abs_v_wgrid
1180
1181!------------------------------------------------------------------------------!
1182! Description:
1183! ------------
1184!> Compute the absolute value of the horizontal velocity (relative to the
1185!> surface) for horizontal surface elements. This is required by all methods.
1186!------------------------------------------------------------------------------!
1187    SUBROUTINE calc_uvw_abs_v_sgrid
1188
1189       IMPLICIT NONE
1190
1191       INTEGER(iwp) ::  i   !< running index x direction
1192       INTEGER(iwp) ::  j   !< running index y direction
1193       INTEGER(iwp) ::  k   !< running index z direction
1194       INTEGER(iwp) ::  m   !< running index surface elements
1195
1196       REAL(wp)     ::  u_i !< u-component on scalar grid
1197       REAL(wp)     ::  v_i !< v-component on scalar grid
1198       REAL(wp)     ::  w_i !< w-component on scalar grid
1199
1200!
1201!--    North- (l=0) and south-facing (l=1) walls
1202       IF ( l == 0  .OR.  l == 1 )  THEN
1203          DO  m = 1, surf%ns
1204             i   = surf%i(m)           
1205             j   = surf%j(m)
1206             k   = surf%k(m)
1207
1208             u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
1209             v_i = 0.0_wp
1210             w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
1211
1212             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
1213          ENDDO
1214!
1215!--    East- (l=2) and west-facing (l=3) walls
1216       ELSE
1217          DO  m = 1, surf%ns
1218             i   = surf%i(m)           
1219             j   = surf%j(m)
1220             k   = surf%k(m)
1221
1222             u_i = 0.0_wp
1223             v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
1224             w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
1225
1226             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
1227          ENDDO
1228       ENDIF 
1229
1230    END SUBROUTINE calc_uvw_abs_v_sgrid
1231
1232
1233!------------------------------------------------------------------------------!
1234! Description:
1235! ------------
1236!> Calculate the Obukhov length (L) and Richardson flux number (z/L)
1237!------------------------------------------------------------------------------!
1238    SUBROUTINE calc_ol
1239
1240       IMPLICIT NONE
1241
1242       INTEGER(iwp) ::  iter    !< Newton iteration step
1243       INTEGER(iwp) ::  li      !< look index
1244       INTEGER(iwp) ::  m       !< loop variable over all horizontal wall elements
1245
1246       REAL(wp)     :: f,      & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
1247                       f_d_ol, & !< Derivative of f
1248                       ol_l,   & !< Lower bound of L for Newton iteration
1249                       ol_m,   & !< Previous value of L for Newton iteration
1250                       ol_old, & !< Previous time step value of L
1251                       ol_u      !< Upper bound of L for Newton iteration
1252
1253       IF ( TRIM( most_method ) /= 'circular' )  THEN
1254!
1255!--       Evaluate bulk Richardson number (calculation depends on
1256!--       definition based on setting of boundary conditions
1257          IF ( ibc_pt_b /= 1 )  THEN
1258             IF ( humidity )  THEN
1259                !$OMP PARALLEL DO PRIVATE( z_mo )
1260                DO  m = 1, surf%ns
1261
1262                   z_mo = surf%z_mo(m)
1263
1264                   surf%rib(m) = g * z_mo                                      &
1265                                 * ( surf%vpt1(m) - surf%vpt_surface(m) )      &
1266                                 / ( surf%uvw_abs(m)**2 * surf%vpt1(m)         &
1267                                 + 1.0E-20_wp )
1268                ENDDO
1269             ELSE
1270                !$OMP PARALLEL DO PRIVATE( z_mo )
1271                DO  m = 1, surf%ns
1272
1273                   z_mo = surf%z_mo(m)
1274
1275                   surf%rib(m) = g * z_mo                                      &
1276                                 * ( surf%pt1(m) - surf%pt_surface(m) )        &
1277                                 / ( surf%uvw_abs(m)**2 * surf%pt1(m)          &
1278                                 + 1.0E-20_wp )
1279                ENDDO
1280             ENDIF
1281          ELSE
1282             IF ( humidity )  THEN
1283                !$OMP PARALLEL DO PRIVATE( k, z_mo )
1284                DO  m = 1, surf%ns
1285
1286                   k   = surf%k(m)
1287
1288                   z_mo = surf%z_mo(m)
1289
1290                   surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
1291                           * surf%qv1(m) ) * surf%shf(m) + 0.61_wp             &
1292                           * surf%pt1(m) * surf%qsws(m) ) *                    &
1293                             drho_air_zw(k-1)                /                 &
1294                         ( surf%uvw_abs(m)**3 * surf%vpt1(m) * kappa**2        &
1295                           + 1.0E-20_wp )
1296                ENDDO
1297             ELSE
1298                !$OMP PARALLEL DO PRIVATE( k, z_mo )
1299                DO  m = 1, surf%ns
1300
1301                   k   = surf%k(m)
1302
1303                   z_mo = surf%z_mo(m)
1304
1305                   surf%rib(m) = - g * z_mo * surf%shf(m) *                    &
1306                                        drho_air_zw(k-1)            /          &
1307                        ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2          &
1308                           + 1.0E-20_wp )
1309                ENDDO
1310             ENDIF
1311          ENDIF
1312
1313       ENDIF
1314
1315
1316!
1317!--    Calculate the Obukhov length either using a Newton iteration
1318!--    method, via a lookup table, or using the old circular way
1319       IF ( TRIM( most_method ) == 'newton' )  THEN
1320
1321          DO  m = 1, surf%ns
1322
1323             i   = surf%i(m)           
1324             j   = surf%j(m)
1325
1326             z_mo = surf%z_mo(m)
1327
1328!
1329!--          Store current value in case the Newton iteration fails
1330             ol_old = surf%ol(m)
1331
1332!
1333!--          Ensure that the bulk Richardson number and the Obukhov
1334!--          length have the same sign
1335             IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.                      &
1336                  ABS( surf%ol(m) ) == ol_max )  THEN
1337                IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
1338                IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
1339             ENDIF
1340!
1341!--          Iteration to find Obukhov length
1342             iter = 0
1343             DO
1344                iter = iter + 1
1345!
1346!--             In case of divergence, use the value of the previous time step
1347                IF ( iter > 1000 )  THEN
1348                   surf%ol(m) = ol_old
1349                   EXIT
1350                ENDIF
1351
1352                ol_m = surf%ol(m)
1353                ol_l = ol_m - 0.001_wp * ol_m
1354                ol_u = ol_m + 0.001_wp * ol_m
1355
1356
1357                IF ( ibc_pt_b /= 1 )  THEN
1358!
1359!--                Calculate f = Ri - [...]/[...]^2 = 0
1360                   f = surf%rib(m) - ( z_mo / ol_m ) * (                       &
1361                                                 LOG( z_mo / surf%z0h(m) )     &
1362                                                 - psi_h( z_mo / ol_m )        &
1363                                                 + psi_h( surf%z0h(m) /        &
1364                                                          ol_m )               &
1365                                                               )               &
1366                                              / ( LOG( z_mo / surf%z0(m) )     &
1367                                                 - psi_m( z_mo / ol_m )        &
1368                                                 + psi_m( surf%z0(m) /         &
1369                                                          ol_m )               &
1370                                                 )**2
1371
1372!
1373!--                 Calculate df/dL
1374                    f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo /               &
1375                                                             surf%z0h(m) )     &
1376                                            - psi_h( z_mo / ol_u )             &
1377                                            + psi_h( surf%z0h(m) / ol_u )      &
1378                                              )                                &
1379                                            / ( LOG( z_mo / surf%z0(m) )       &
1380                                            - psi_m( z_mo / ol_u )             &
1381                                            + psi_m( surf%z0(m) / ol_u )       &
1382                                              )**2                             &
1383                           + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) )     &
1384                                            - psi_h( z_mo / ol_l )             &
1385                                            + psi_h( surf%z0h(m) / ol_l )      &
1386                                              )                                &
1387                                            / ( LOG( z_mo / surf%z0(m) )       &
1388                                            - psi_m( z_mo / ol_l )             &
1389                                            + psi_m( surf%z0(m) / ol_l )       &
1390                                              )**2                             &
1391                             ) / ( ol_u - ol_l )
1392                ELSE
1393!
1394!--                Calculate f = Ri - 1 /[...]^3 = 0
1395                   f = surf%rib(m) - ( z_mo / ol_m ) /                         &
1396                                                ( LOG( z_mo / surf%z0(m) )     &
1397                                            - psi_m( z_mo / ol_m )             &
1398                                            + psi_m( surf%z0(m) / ol_m )       &
1399                                                )**3
1400
1401!
1402!--                Calculate df/dL
1403                   f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo /                &
1404                                                            surf%z0(m) )       &
1405                                            - psi_m( z_mo / ol_u )             &
1406                                            + psi_m( surf%z0(m) / ol_u )       &
1407                                                     )**3                      &
1408                           + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )      &
1409                                            - psi_m( z_mo / ol_l )             &
1410                                            + psi_m( surf%z0(m) / ol_l )       &
1411                                               )**3                            &
1412                                  ) / ( ol_u - ol_l )
1413                ENDIF
1414!
1415!--             Calculate new L
1416                surf%ol(m) = ol_m - f / f_d_ol
1417
1418!
1419!--             Ensure that the bulk Richardson number and the Obukhov
1420!--             length have the same sign and ensure convergence.
1421                IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
1422!
1423!--             If unrealistic value occurs, set L to the maximum
1424!--             value that is allowed
1425                IF ( ABS( surf%ol(m) ) > ol_max )  THEN
1426                   surf%ol(m) = ol_max
1427                   EXIT
1428                ENDIF
1429!
1430!--             Check for convergence
1431                IF ( ABS( ( surf%ol(m) - ol_m ) /                              &
1432                     surf%ol(m) ) < 1.0E-4_wp )  THEN
1433                   EXIT
1434                ELSE
1435                   CYCLE
1436                ENDIF
1437
1438             ENDDO
1439          ENDDO
1440
1441       ELSEIF ( TRIM( most_method ) == 'lookup' )  THEN
1442
1443          !$OMP PARALLEL DO PRIVATE( i, j, z_mo, li ) FIRSTPRIVATE( li_bnd ) LASTPRIVATE( li_bnd )
1444          DO  m = 1, surf%ns
1445
1446             i   = surf%i(m)           
1447             j   = surf%j(m)
1448!
1449!--          If the bulk Richardson number is outside the range of the lookup
1450!--          table, set it to the exceeding threshold value
1451             IF ( surf%rib(m) < rib_min )  surf%rib(m) = rib_min
1452             IF ( surf%rib(m) > rib_max )  surf%rib(m) = rib_max
1453
1454!
1455!--          Find the correct index bounds for linear interpolation. As the
1456!--          Richardson number will not differ very much from time step to
1457!--          time step , use the index from the last step and search in the
1458!--          correct direction
1459             li = li_bnd
1460             IF ( rib_tab(li) - surf%rib(m) > 0.0_wp )  THEN
1461                DO  WHILE ( rib_tab(li-1) - surf%rib(m) > 0.0_wp  .AND.  li > 0 )
1462                   li = li-1
1463                ENDDO
1464             ELSE
1465                DO  WHILE ( rib_tab(li) - surf%rib(m) < 0.0_wp                 &
1466                           .AND.  li < num_steps-1 )
1467                   li = li+1
1468                ENDDO
1469             ENDIF
1470             li_bnd = li
1471
1472!
1473!--          Linear interpolation to find the correct value of z/L
1474             surf%ol(m) = ( ol_tab(li-1) + ( ol_tab(li) - ol_tab(li-1) )       &
1475                         / (  rib_tab(li) - rib_tab(li-1) )                    &
1476                         * ( surf%rib(m) - rib_tab(li-1) ) )
1477
1478          ENDDO
1479
1480       ELSEIF ( TRIM( most_method ) == 'circular' )  THEN
1481
1482          IF ( .NOT. humidity )  THEN
1483             !$OMP PARALLEL DO PRIVATE( z_mo )
1484             DO  m = 1, surf%ns
1485
1486                z_mo = surf%z_mo(m)
1487
1488                surf%ol(m) =  ( surf%pt1(m) *  surf%us(m)**2 ) /               &
1489                                       ( kappa * g *                           &
1490                                         surf%ts(m) + 1E-30_wp )
1491!
1492!--             Limit the value range of the Obukhov length.
1493!--             This is necessary for very small velocities (u,v --> 1), because
1494!--             the absolute value of ol can then become very small, which in
1495!--             consequence would result in very large shear stresses and very
1496!--             small momentum fluxes (both are generally unrealistic).
1497                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
1498                   surf%ol(m) = z_mo / zeta_min
1499                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
1500                   surf%ol(m) = z_mo / zeta_max
1501
1502             ENDDO
1503          ELSEIF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
1504             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1505             DO  m = 1, surf%ns
1506
1507                i   = surf%i(m)           
1508                j   = surf%j(m)
1509                k   = surf%k(m)
1510
1511                z_mo = surf%z_mo(m)
1512
1513
1514                surf%ol(m) =  ( surf%vpt1(m) * surf%us(m)**2 ) /               &
1515                    ( kappa * g * ( surf%ts(m) +                               &
1516                         0.61_wp * surf%pt1(m) * surf%us(m)                    &
1517                       + 0.61_wp * surf%qv1(m) * surf%ts(m) -                  &
1518                                   surf%ts(m)  * ql(k,j,i) ) + 1E-30_wp )
1519!
1520!--             Limit the value range of the Obukhov length.
1521!--             This is necessary for very small velocities (u,v --> 1), because
1522!--             the absolute value of ol can then become very small, which in
1523!--             consequence would result in very large shear stresses and very
1524!--             small momentum fluxes (both are generally unrealistic).
1525                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
1526                   surf%ol(m) = z_mo / zeta_min
1527                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
1528                   surf%ol(m) = z_mo / zeta_max
1529
1530             ENDDO
1531          ELSE
1532
1533             !$OMP PARALLEL DO PRIVATE( z_mo )
1534             DO  m = 1, surf%ns
1535
1536                z_mo = surf%z_mo(m)
1537
1538                surf%ol(m) =  ( surf%vpt1(m) *  surf%us(m)**2 ) /              &
1539                    ( kappa * g * ( surf%ts(m) + 0.61_wp * surf%pt1(m) *       &
1540                                    surf%qs(m) + 0.61_wp * surf%qv1(m) *       &
1541                                    surf%ts(m) ) + 1E-30_wp )
1542
1543!
1544!--             Limit the value range of the Obukhov length.
1545!--             This is necessary for very small velocities (u,v --> 1), because
1546!--             the absolute value of ol can then become very small, which in
1547!--             consequence would result in very large shear stresses and very
1548!--             small momentum fluxes (both are generally unrealistic).
1549                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
1550                   surf%ol(m) = z_mo / zeta_min
1551                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
1552                   surf%ol(m) = z_mo / zeta_max
1553
1554             ENDDO
1555
1556          ENDIF
1557
1558       ENDIF
1559
1560    END SUBROUTINE calc_ol
1561
1562!
1563!-- Calculate friction velocity u*
1564    SUBROUTINE calc_us
1565
1566       IMPLICIT NONE
1567
1568       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1569
1570!
1571!--    Compute u* at horizontal surfaces at the scalars' grid points
1572       IF ( .NOT. surf_vertical )  THEN
1573!
1574!--       Compute u* at upward-facing surfaces
1575          IF ( .NOT. downward )  THEN
1576             !$OMP PARALLEL  DO PRIVATE( z_mo )
1577             DO  m = 1, surf%ns
1578
1579                z_mo = surf%z_mo(m)
1580!
1581!--             Compute u* at the scalars' grid points
1582                surf%us(m) = kappa * surf%uvw_abs(m) /                         &
1583                            ( LOG( z_mo / surf%z0(m) )                         &
1584                           - psi_m( z_mo / surf%ol(m) )                        &
1585                           + psi_m( surf%z0(m) / surf%ol(m) ) )
1586   
1587             ENDDO
1588!
1589!--       Compute u* at downward-facing surfaces. This case, do not consider
1590!--       any stability.
1591          ELSE
1592             !$OMP PARALLEL  DO PRIVATE( z_mo )
1593             DO  m = 1, surf%ns
1594
1595                z_mo = surf%z_mo(m)
1596!
1597!--             Compute u* at the scalars' grid points
1598                surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
1599   
1600             ENDDO
1601          ENDIF
1602!
1603!--    Compute u* at vertical surfaces at the u/v/v grid, respectively.
1604!--    No stability is considered in this case.
1605       ELSE
1606          !$OMP PARALLEL DO PRIVATE( z_mo )
1607          DO  m = 1, surf%ns
1608             z_mo = surf%z_mo(m)
1609
1610             surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
1611          ENDDO
1612       ENDIF
1613
1614    END SUBROUTINE calc_us
1615
1616!
1617!-- Calculate potential temperature, specific humidity, and virtual potential
1618!-- temperature at first grid level
1619    SUBROUTINE calc_pt_q
1620
1621       IMPLICIT NONE
1622
1623       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1624
1625       !$OMP PARALLEL DO PRIVATE( i, j, k )
1626       DO  m = 1, surf%ns
1627
1628          i   = surf%i(m)           
1629          j   = surf%j(m)
1630          k   = surf%k(m)
1631
1632          IF ( bulk_cloud_model ) THEN
1633             surf%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
1634             surf%qv1(m) = q(k,j,i) - ql(k,j,i)
1635          ELSEIF( cloud_droplets ) THEN
1636             surf%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
1637             surf%qv1(m) = q(k,j,i) 
1638          ELSE
1639             surf%pt1(m) = pt(k,j,i)
1640             IF ( humidity )  THEN
1641                surf%qv1(m) = q(k,j,i)
1642             ELSE
1643                surf%qv1(m) = 0.0_wp
1644             ENDIF
1645          ENDIF
1646
1647          IF ( humidity )  THEN
1648             surf%vpt1(m) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
1649          ENDIF
1650         
1651       ENDDO
1652
1653    END SUBROUTINE calc_pt_q
1654
1655
1656!
1657!-- Set potential temperature at surface grid level.
1658!-- ( only for upward-facing surfs )
1659    SUBROUTINE calc_pt_surface
1660
1661       IMPLICIT NONE
1662
1663       INTEGER(iwp) ::  k_off    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
1664       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1665       
1666       k_off = surf%koff
1667       !$OMP PARALLEL DO PRIVATE( i, j, k )
1668       DO  m = 1, surf%ns
1669
1670          i   = surf%i(m)           
1671          j   = surf%j(m)
1672          k   = surf%k(m)
1673
1674          surf%pt_surface(m) = pt(k+k_off,j,i)
1675
1676       ENDDO
1677
1678    END SUBROUTINE calc_pt_surface
1679
1680!
1681!-- Set mixing ratio at surface grid level. ( Only for upward-facing surfs. )
1682    SUBROUTINE calc_q_surface
1683
1684       IMPLICIT NONE
1685
1686       INTEGER(iwp) ::  k_off   !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
1687       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1688       
1689       k_off = surf%koff
1690       !$OMP PARALLEL DO PRIVATE( i, j, k )
1691       DO  m = 1, surf%ns
1692
1693          i   = surf%i(m)           
1694          j   = surf%j(m)
1695          k   = surf%k(m)
1696
1697          surf%q_surface(m) = q(k+k_off,j,i)
1698
1699       ENDDO
1700
1701    END SUBROUTINE calc_q_surface
1702   
1703!
1704!-- Set virtual potential temperature at surface grid level.
1705!-- ( only for upward-facing surfs )
1706    SUBROUTINE calc_vpt_surface
1707
1708       IMPLICIT NONE
1709
1710       INTEGER(iwp) ::  k_off    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
1711       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1712       
1713       k_off = surf%koff
1714       !$OMP PARALLEL DO PRIVATE( i, j, k )
1715       DO  m = 1, surf%ns
1716
1717          i   = surf%i(m)           
1718          j   = surf%j(m)
1719          k   = surf%k(m)
1720
1721          surf%vpt_surface(m) = vpt(k+k_off,j,i)
1722
1723       ENDDO
1724
1725    END SUBROUTINE calc_vpt_surface
1726   
1727!
1728!-- Calculate the other MOST scaling parameters theta*, q*, (qc*, qr*, nc*, nr*)
1729    SUBROUTINE calc_scaling_parameters
1730
1731       IMPLICIT NONE
1732
1733
1734       INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements
1735       INTEGER(iwp)  ::  lsp     !< running index for chemical species
1736!
1737!--    Compute theta* at horizontal surfaces
1738       IF ( constant_heatflux  .AND.  .NOT. surf_vertical )  THEN
1739!
1740!--       For a given heat flux in the surface layer:
1741
1742          !$OMP PARALLEL DO PRIVATE( i, j, k )
1743          DO  m = 1, surf%ns
1744
1745             i   = surf%i(m)           
1746             j   = surf%j(m)
1747             k   = surf%k(m)
1748
1749             surf%ts(m) = -surf%shf(m) * drho_air_zw(k-1) /                    &
1750                                  ( surf%us(m) + 1E-30_wp )
1751
1752!
1753!--          ts must be limited, because otherwise overflow may occur in case
1754!--          of us=0 when computing ol further below
1755             IF ( surf%ts(m) < -1.05E5_wp )  surf%ts(m) = -1.0E5_wp
1756             IF ( surf%ts(m) >  1.0E5_wp  )  surf%ts(m) =  1.0E5_wp
1757
1758          ENDDO
1759
1760       ELSEIF ( .NOT. surf_vertical ) THEN
1761!
1762!--       For a given surface temperature:
1763          IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
1764
1765             !$OMP PARALLEL DO PRIVATE( i, j, k )
1766             DO  m = 1, surf%ns
1767                i   = surf%i(m)           
1768                j   = surf%j(m)
1769                k   = surf%k(m)
1770
1771                pt(k-1,j,i) = pt_surface
1772             ENDDO
1773          ENDIF
1774
1775          !$OMP PARALLEL DO PRIVATE( z_mo )
1776          DO  m = 1, surf%ns   
1777
1778             z_mo = surf%z_mo(m)
1779
1780             surf%ts(m) = kappa * ( surf%pt1(m) - surf%pt_surface(m) )      &
1781                                  / ( LOG( z_mo / surf%z0h(m) )             &
1782                                      - psi_h( z_mo / surf%ol(m) )          &
1783                                      + psi_h( surf%z0h(m) / surf%ol(m) ) )
1784
1785          ENDDO
1786
1787       ENDIF
1788!
1789!--    Compute theta* at vertical surfaces. This is only required in case of
1790!--    land-surface model, in order to compute aerodynamical resistance.
1791       IF ( surf_vertical )  THEN
1792          !$OMP PARALLEL DO PRIVATE( i, j )
1793          DO  m = 1, surf%ns
1794
1795             i          =  surf%i(m)           
1796             j          =  surf%j(m)
1797             surf%ts(m) = -surf%shf(m) / ( surf%us(m) + 1E-30_wp )
1798!
1799!--          ts must be limited, because otherwise overflow may occur in case
1800!--          of us=0 when computing ol further below
1801             IF ( surf%ts(m) < -1.05E5_wp )  surf%ts(m) = -1.0E5_wp
1802             IF ( surf%ts(m) >  1.0E5_wp  )  surf%ts(m) =  1.0E5_wp
1803
1804          ENDDO
1805       ENDIF
1806
1807!
1808!--    If required compute q* at horizontal surfaces
1809       IF ( humidity )  THEN
1810          IF ( constant_waterflux  .AND.  .NOT. surf_vertical )  THEN
1811!
1812!--          For a given water flux in the surface layer
1813             !$OMP PARALLEL DO PRIVATE( i, j, k )
1814             DO  m = 1, surf%ns 
1815
1816                i   = surf%i(m)           
1817                j   = surf%j(m)
1818                k   = surf%k(m)
1819                surf%qs(m) = -surf%qsws(m) * drho_air_zw(k-1) /                &
1820                                               ( surf%us(m) + 1E-30_wp )
1821
1822             ENDDO
1823
1824          ELSEIF ( .NOT. surf_vertical )  THEN
1825             coupled_run = ( coupling_mode == 'atmosphere_to_ocean'  .AND.     &
1826                             run_coupled )
1827
1828             IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
1829                !$OMP PARALLEL DO PRIVATE( i, j, k )
1830                DO  m = 1, surf%ns 
1831
1832                   i   = surf%i(m)           
1833                   j   = surf%j(m)
1834                   k   = surf%k(m)
1835                   q(k-1,j,i) = q_surface
1836                   
1837                ENDDO
1838             ENDIF
1839
1840!
1841!--          Assume saturation for atmosphere coupled to ocean (but not
1842!--          in case of precursor runs)
1843             IF ( coupled_run )  THEN
1844                !$OMP PARALLEL DO PRIVATE( i, j, k, e_s )
1845                DO  m = 1, surf%ns   
1846                   i   = surf%i(m)           
1847                   j   = surf%j(m)
1848                   k   = surf%k(m)
1849                   e_s = 6.1_wp * &
1850                              EXP( 0.07_wp * ( MIN(pt(k-1,j,i),pt(k,j,i))      &
1851                                               - 273.15_wp ) )
1852                   q(k-1,j,i) = rd_d_rv * e_s / ( surface_pressure - e_s )
1853                ENDDO
1854             ENDIF
1855
1856             IF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
1857               !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1858                DO  m = 1, surf%ns   
1859
1860                   i   = surf%i(m)           
1861                   j   = surf%j(m)
1862                   k   = surf%k(m)
1863   
1864                   z_mo = surf%z_mo(m)
1865
1866                   surf%qs(m) = kappa * ( surf%qv1(m) - surf%q_surface(m) )    &
1867                                        / ( LOG( z_mo / surf%z0q(m) )          &
1868                                            - psi_h( z_mo / surf%ol(m) )       &
1869                                            + psi_h( surf%z0q(m) /             &
1870                                                     surf%ol(m) ) )
1871                ENDDO
1872             ELSE
1873                !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1874                DO  m = 1, surf%ns   
1875
1876                   i   = surf%i(m)           
1877                   j   = surf%j(m)
1878                   k   = surf%k(m)
1879   
1880                   z_mo = surf%z_mo(m)
1881
1882                   surf%qs(m) = kappa * ( q(k,j,i) - q(k-1,j,i) )              &
1883                                        / ( LOG( z_mo / surf%z0q(m) )          &
1884                                            - psi_h( z_mo / surf%ol(m) )       &
1885                                            + psi_h( surf%z0q(m) /             &
1886                                                     surf%ol(m) ) )
1887                ENDDO
1888             ENDIF
1889          ENDIF
1890!
1891!--       Compute q* at vertical surfaces
1892          IF ( surf_vertical )  THEN
1893             !$OMP PARALLEL DO PRIVATE( i, j )
1894             DO  m = 1, surf%ns 
1895
1896                i          =  surf%i(m)           
1897                j          =  surf%j(m)
1898                surf%qs(m) = -surf%qsws(m) / ( surf%us(m) + 1E-30_wp )
1899
1900             ENDDO
1901          ENDIF
1902       ENDIF
1903       
1904!
1905!--    If required compute s*
1906       IF ( passive_scalar )  THEN
1907!
1908!--       At horizontal surfaces
1909          IF ( constant_scalarflux  .AND.  .NOT. surf_vertical )  THEN
1910!
1911!--          For a given scalar flux in the surface layer
1912             !$OMP PARALLEL DO PRIVATE( i, j )
1913             DO  m = 1, surf%ns 
1914                i   = surf%i(m)           
1915                j   = surf%j(m)
1916                surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
1917             ENDDO
1918          ENDIF
1919!
1920!--       At vertical surfaces
1921          IF ( surf_vertical )  THEN
1922             !$OMP PARALLEL DO PRIVATE( i, j )
1923             DO  m = 1, surf%ns 
1924                i          =  surf%i(m)           
1925                j          =  surf%j(m)
1926                surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
1927             ENDDO
1928          ENDIF
1929       ENDIF
1930
1931!
1932!--    If required compute cs* (chemical species)
1933       IF ( air_chemistry  )  THEN 
1934!
1935!--       At horizontal surfaces                             
1936          DO  lsp = 1, nvar
1937             IF ( constant_csflux(lsp)  .AND.  .NOT.  surf_vertical )  THEN
1938!--             For a given chemical species' flux in the surface layer
1939                !$OMP PARALLEL DO PRIVATE( i, j )
1940                DO  m = 1, surf%ns 
1941                   i   = surf%i(m)           
1942                   j   = surf%j(m)
1943                   surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
1944                ENDDO
1945             ENDIF
1946          ENDDO
1947!
1948!--       At vertical surfaces
1949          IF ( surf_vertical )  THEN
1950             DO  lsp = 1, nvar
1951                !$OMP PARALLEL DO PRIVATE( i, j )
1952                DO  m = 1, surf%ns 
1953                   i   = surf%i(m)           
1954                   j   = surf%j(m)
1955                   surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
1956                ENDDO
1957             ENDDO
1958          ENDIF
1959       ENDIF
1960
1961!
1962!--    If required compute qc* and nc*
1963       IF ( bulk_cloud_model  .AND.  microphysics_morrison  .AND.              &
1964            .NOT. surf_vertical )  THEN
1965          !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1966          DO  m = 1, surf%ns   
1967
1968             i    = surf%i(m)           
1969             j    = surf%j(m)
1970             k    = surf%k(m)
1971
1972             z_mo = surf%z_mo(m)
1973
1974             surf%qcs(m) = kappa * ( qc(k,j,i) - qc(k-1,j,i) )                 &
1975                                 / ( LOG( z_mo / surf%z0q(m) )                 &
1976                                 - psi_h( z_mo / surf%ol(m) )                  &
1977                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
1978
1979             surf%ncs(m) = kappa * ( nc(k,j,i) - nc(k-1,j,i) )                 &
1980                                 / ( LOG( z_mo / surf%z0q(m) )                 &
1981                                 - psi_h( z_mo / surf%ol(m) )                  &
1982                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
1983          ENDDO
1984
1985       ENDIF
1986
1987!
1988!--    If required compute qr* and nr*
1989       IF ( bulk_cloud_model  .AND.  microphysics_seifert  .AND.               &
1990            .NOT. surf_vertical )  THEN
1991          !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1992          DO  m = 1, surf%ns   
1993
1994             i    = surf%i(m)           
1995             j    = surf%j(m)
1996             k    = surf%k(m)
1997
1998             z_mo = surf%z_mo(m)
1999
2000             surf%qrs(m) = kappa * ( qr(k,j,i) - qr(k-1,j,i) )                 &
2001                                 / ( LOG( z_mo / surf%z0q(m) )                 &
2002                                 - psi_h( z_mo / surf%ol(m) )                  &
2003                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
2004
2005             surf%nrs(m) = kappa * ( nr(k,j,i) - nr(k-1,j,i) )                 &
2006                                 / ( LOG( z_mo / surf%z0q(m) )                 &
2007                                 - psi_h( z_mo / surf%ol(m) )                  &
2008                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
2009          ENDDO
2010
2011       ENDIF
2012
2013    END SUBROUTINE calc_scaling_parameters
2014
2015
2016
2017!
2018!-- Calculate surface fluxes usws, vsws, shf, qsws, (qcsws, qrsws, ncsws, nrsws)
2019    SUBROUTINE calc_surface_fluxes
2020
2021       IMPLICIT NONE
2022
2023       INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements
2024       INTEGER(iwp)  ::  lsp     !< running index for chemical species
2025
2026       REAL(wp)                            ::  dum     !< dummy to precalculate logarithm
2027       REAL(wp)                            ::  flag_u  !< flag indicating u-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
2028       REAL(wp)                            ::  flag_v  !< flag indicating v-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
2029       REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_i     !< u-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
2030       REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_i     !< v-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
2031       REAL(wp), DIMENSION(:), ALLOCATABLE ::  w_i     !< w-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
2032
2033!
2034!--    Calcuate surface fluxes at horizontal walls
2035       IF ( .NOT. surf_vertical )  THEN
2036!
2037!--       Compute u'w' for the total model domain at upward-facing surfaces.
2038!--       First compute the corresponding component of u* and square it.
2039          IF ( .NOT. downward )  THEN
2040             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
2041             DO  m = 1, surf%ns 
2042   
2043                i = surf%i(m)           
2044                j = surf%j(m)
2045                k = surf%k(m)
2046
2047                z_mo = surf%z_mo(m)
2048
2049                surf%usws(m) = kappa * ( u(k,j,i) - u(k-1,j,i) )               &
2050                              / ( LOG( z_mo / surf%z0(m) )                     &
2051                                  - psi_m( z_mo / surf%ol(m) )                 &
2052                                  + psi_m( surf%z0(m) / surf%ol(m) ) )
2053!
2054!--             Please note, the computation of usws is not fully accurate. Actually
2055!--             a further interpolation of us onto the u-grid, where usws is defined,
2056!--             is required. However, this is not done as this would require several
2057!--             data transfers between 2D-grid and the surf-type.
2058!--             The impact of the missing interpolation is negligible as several
2059!--             tests had shown.
2060!--             Same also for ol. 
2061                surf%usws(m) = -surf%usws(m) * surf%us(m) * rho_air_zw(k-1)
2062
2063             ENDDO
2064!
2065!--       At downward-facing surfaces
2066          ELSE
2067             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
2068             DO  m = 1, surf%ns 
2069   
2070                i = surf%i(m)           
2071                j = surf%j(m)
2072                k = surf%k(m)
2073
2074                z_mo = surf%z_mo(m)
2075
2076                surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) )
2077                surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k)
2078
2079             ENDDO     
2080          ENDIF
2081
2082!
2083!--       Compute v'w' for the total model domain.
2084!--       First compute the corresponding component of u* and square it.
2085!--       Upward-facing surfaces
2086          IF ( .NOT. downward )  THEN
2087             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
2088             DO  m = 1, surf%ns 
2089                i = surf%i(m)           
2090                j = surf%j(m)
2091                k = surf%k(m)
2092
2093                z_mo = surf%z_mo(m)
2094
2095                surf%vsws(m) = kappa * ( v(k,j,i) - v(k-1,j,i) )               &
2096                           / ( LOG( z_mo / surf%z0(m) )                        &
2097                               - psi_m( z_mo / surf%ol(m) )                    &
2098                               + psi_m( surf%z0(m) / surf%ol(m) ) )
2099!
2100!--             Please note, the computation of vsws is not fully accurate. Actually
2101!--             a further interpolation of us onto the v-grid, where vsws is defined,
2102!--             is required. However, this is not done as this would require several
2103!--             data transfers between 2D-grid and the surf-type.
2104!--             The impact of the missing interpolation is negligible as several
2105!--             tests had shown.
2106!--             Same also for ol. 
2107                surf%vsws(m) = -surf%vsws(m) * surf%us(m) * rho_air_zw(k-1)
2108             ENDDO
2109!
2110!--       Downward-facing surfaces
2111          ELSE
2112             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
2113             DO  m = 1, surf%ns 
2114                i = surf%i(m)           
2115                j = surf%j(m)
2116                k = surf%k(m)
2117
2118                z_mo = surf%z_mo(m)
2119
2120                surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) )
2121                surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k)
2122             ENDDO
2123          ENDIF
2124!
2125!--       Compute the vertical kinematic heat flux
2126          IF (  .NOT.  constant_heatflux  .AND.  ( ( time_since_reference_point&
2127               <=  skip_time_do_lsm  .AND. simulated_time > 0.0_wp ) .OR.      &
2128               .NOT.  land_surface )  .AND.  .NOT. urban_surface  .AND.        &
2129               .NOT. downward )  THEN
2130             !$OMP PARALLEL DO PRIVATE( i, j, k )
2131             DO  m = 1, surf%ns
2132                i    = surf%i(m)           
2133                j    = surf%j(m)
2134                k    = surf%k(m)
2135                surf%shf(m) = -surf%ts(m) * surf%us(m) * rho_air_zw(k-1)
2136             ENDDO
2137          ENDIF
2138!
2139!--       Compute the vertical water flux
2140          IF (  .NOT.  constant_waterflux  .AND.  humidity  .AND.              &
2141               ( ( time_since_reference_point <= skip_time_do_lsm  .AND.       &
2142               simulated_time > 0.0_wp ) .OR.  .NOT.  land_surface  )  .AND.   &
2143               .NOT.  urban_surface  .AND.  .NOT. downward )  THEN
2144             !$OMP PARALLEL DO PRIVATE( i, j, k )
2145             DO  m = 1, surf%ns
2146                i    = surf%i(m)           
2147                j    = surf%j(m)
2148                k    = surf%k(m)
2149                surf%qsws(m) = -surf%qs(m) * surf%us(m) * rho_air_zw(k-1)
2150             ENDDO
2151          ENDIF
2152!
2153!--       Compute the vertical scalar flux
2154          IF (  .NOT.  constant_scalarflux  .AND.  passive_scalar  .AND.       &
2155                .NOT.  downward )  THEN
2156             !$OMP PARALLEL DO PRIVATE( i, j )
2157             DO  m = 1, surf%ns   
2158
2159                i    = surf%i(m)           
2160                j    = surf%j(m)
2161                surf%ssws(m) = -surf%ss(m) * surf%us(m)
2162
2163             ENDDO
2164          ENDIF   
2165!
2166!--       Compute the vertical chemical species' flux
2167          DO  lsp = 1, nvar
2168             IF (  .NOT.  constant_csflux(lsp)  .AND.  air_chemistry  .AND.    &
2169                   .NOT.  downward )  THEN
2170                !$OMP PARALLEL DO PRIVATE( i, j )
2171                DO  m = 1, surf%ns   
2172                   i     = surf%i(m)           
2173                   j     = surf%j(m)
2174                   surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m)
2175                ENDDO
2176             ENDIF   
2177          ENDDO
2178
2179!
2180!--       Compute (turbulent) fluxes of cloud water content and cloud drop conc.
2181          IF ( bulk_cloud_model  .AND.  microphysics_morrison  .AND.           &
2182               .NOT. downward)  THEN
2183             !$OMP PARALLEL DO PRIVATE( i, j )
2184             DO  m = 1, surf%ns   
2185
2186                i    = surf%i(m)           
2187                j    = surf%j(m)
2188
2189                surf%qcsws(m) = -surf%qcs(m) * surf%us(m)
2190                surf%ncsws(m) = -surf%ncs(m) * surf%us(m)
2191             ENDDO
2192          ENDIF   
2193!
2194!--       Compute (turbulent) fluxes of rain water content and rain drop conc.
2195          IF ( bulk_cloud_model  .AND.  microphysics_seifert  .AND.            &
2196               .NOT. downward)  THEN
2197             !$OMP PARALLEL DO PRIVATE( i, j )
2198             DO  m = 1, surf%ns   
2199
2200                i    = surf%i(m)           
2201                j    = surf%j(m)
2202
2203                surf%qrsws(m) = -surf%qrs(m) * surf%us(m)
2204                surf%nrsws(m) = -surf%nrs(m) * surf%us(m)
2205             ENDDO
2206          ENDIF
2207
2208!
2209!--       Bottom boundary condition for the TKE.
2210          IF ( ibc_e_b == 2 )  THEN
2211             !$OMP PARALLEL DO PRIVATE( i, j, k )
2212             DO  m = 1, surf%ns   
2213
2214                i    = surf%i(m)           
2215                j    = surf%j(m)
2216                k    = surf%k(m)
2217
2218                e(k,j,i) = ( surf%us(m) / 0.1_wp )**2
2219!
2220!--             As a test: cm = 0.4
2221!               e(k,j,i) = ( us(j,i) / 0.4_wp )**2
2222                e(k-1,j,i)   = e(k,j,i)
2223
2224             ENDDO
2225          ENDIF
2226!
2227!--    Calcuate surface fluxes at vertical surfaces. No stability is considered.
2228       ELSE
2229!
2230!--       Compute usvs l={0,1} and vsus l={2,3}
2231          IF ( mom_uv )  THEN
2232!
2233!--          Generalize computation by introducing flags. At north- and south-
2234!--          facing surfaces u-component is used, at east- and west-facing
2235!--          surfaces v-component is used.
2236             flag_u = MERGE( 1.0_wp, 0.0_wp, l == 0  .OR.  l == 1 )   
2237             flag_v = MERGE( 1.0_wp, 0.0_wp, l == 2  .OR.  l == 3 )   
2238             !$OMP PARALLEL  DO PRIVATE( i, j, k, z_mo )
2239             DO  m = 1, surf%ns 
2240                i = surf%i(m)           
2241                j = surf%j(m)
2242                k = surf%k(m)
2243
2244                z_mo = surf%z_mo(m)
2245
2246                surf%mom_flux_uv(m) = kappa *                                  &
2247                                ( flag_u * u(k,j,i) + flag_v * v(k,j,i) )  /   &
2248                                                        LOG( z_mo / surf%z0(m) )
2249
2250               surf%mom_flux_uv(m) =                                           &
2251                                    - surf%mom_flux_uv(m) * surf%us(m)
2252             ENDDO
2253          ENDIF
2254!
2255!--       Compute wsus l={0,1} and wsvs l={2,3}
2256          IF ( mom_w )  THEN
2257             !$OMP PARALLEL  DO PRIVATE( i, j, k, z_mo )
2258             DO  m = 1, surf%ns 
2259                i = surf%i(m)           
2260                j = surf%j(m)
2261                k = surf%k(m)
2262
2263                z_mo = surf%z_mo(m)
2264
2265                surf%mom_flux_w(m) = kappa * w(k,j,i) / LOG( z_mo / surf%z0(m) )
2266
2267                surf%mom_flux_w(m) =                                           &
2268                                     - surf%mom_flux_w(m) * surf%us(m)
2269             ENDDO
2270          ENDIF
2271!
2272!--       Compute momentum fluxes used for subgrid-scale TKE production at
2273!--       vertical surfaces. In constrast to the calculated momentum fluxes at
2274!--       vertical surfaces before, which are defined on the u/v/w-grid,
2275!--       respectively), the TKE fluxes are defined at the scalar grid.
2276!--       
2277          IF ( mom_tke )  THEN
2278!
2279!--          Precalculate velocity components at scalar grid point.
2280             ALLOCATE( u_i(1:surf%ns) )
2281             ALLOCATE( v_i(1:surf%ns) )
2282             ALLOCATE( w_i(1:surf%ns) )
2283
2284             IF ( l == 0  .OR.  l == 1 )  THEN
2285                !$OMP PARALLEL DO PRIVATE( i, j, k )
2286                DO  m = 1, surf%ns 
2287                   i = surf%i(m)           
2288                   j = surf%j(m)
2289                   k = surf%k(m)
2290
2291                   u_i(m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
2292                   v_i(m) = 0.0_wp
2293                   w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
2294                ENDDO
2295             ELSE
2296                !$OMP PARALLEL DO PRIVATE( i, j, k )
2297                DO  m = 1, surf%ns 
2298                   i = surf%i(m)           
2299                   j = surf%j(m)
2300                   k = surf%k(m)
2301
2302                   u_i(m) = 0.0_wp
2303                   v_i(m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
2304                   w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
2305                ENDDO
2306             ENDIF
2307
2308             !$OMP PARALLEL DO PRIVATE( i, j, dum, z_mo )
2309             DO  m = 1, surf%ns 
2310                i = surf%i(m)           
2311                j = surf%j(m)
2312
2313                z_mo = surf%z_mo(m)
2314
2315                dum = kappa / LOG( z_mo / surf%z0(m) )
2316!
2317!--             usvs (l=0,1) and vsus (l=2,3)
2318                surf%mom_flux_tke(0,m) = dum * ( u_i(m) + v_i(m) )
2319!
2320!--             wsvs (l=0,1) and wsus (l=2,3)
2321                surf%mom_flux_tke(1,m) = dum * w_i(m)
2322
2323                surf%mom_flux_tke(0:1,m) =                                     &
2324                               - surf%mom_flux_tke(0:1,m) * surf%us(m)
2325             ENDDO
2326!
2327!--          Deallocate temporary arrays
2328             DEALLOCATE( u_i )             
2329             DEALLOCATE( v_i ) 
2330             DEALLOCATE( w_i ) 
2331          ENDIF
2332       ENDIF
2333
2334    END SUBROUTINE calc_surface_fluxes
2335
2336   
2337!------------------------------------------------------------------------------!
2338! Description:
2339! ------------
2340!> Calculates temperature near surface (10 cm) for indoor model or 2 m
2341!> temperature for output
2342!------------------------------------------------------------------------------!
2343    SUBROUTINE calc_pt_near_surface ( z_char )
2344
2345       IMPLICIT NONE
2346
2347       CHARACTER (LEN = *), INTENT(IN)       :: z_char          !< string identifier to identify z level
2348       INTEGER(iwp)                          :: i, j, k, l, m   !< running indices
2349
2350       
2351       SELECT CASE ( z_char)
2352           
2353       
2354          CASE ( '10cm' )
2355
2356             DO  m = 1, surf%ns
2357
2358                i = surf%i(m)
2359                j = surf%j(m)
2360                k = surf%k(m)
2361
2362                surf%pt_10cm(m) = surf%pt_surface(m) + surf%ts(m) / kappa  &
2363                                   * ( log( 0.1_wp /  surf%z0h(m) )                    &
2364                                  - psi_h( 0.1_wp / surf%ol(m) )                    &
2365                                     + psi_h( surf%z0h(m) / surf%ol(m) ) )
2366
2367             ENDDO
2368
2369
2370          CASE ( '2m' )
2371     
2372             DO  m = 1, surf%ns
2373
2374                i = surf%i(m)
2375                j = surf%j(m)
2376                k = surf%k(m)
2377
2378                surf%pt_2m(m) = surf%pt_surface(m) + surf%ts(m) / kappa  &
2379                                   * ( log( 2.0_wp /  surf%z0h(m) )                    &
2380                                     - psi_h( 2.0_wp / surf%ol(m) )                    &
2381                                     + psi_h( surf%z0h(m) / surf%ol(m) ) )
2382
2383             ENDDO
2384         
2385       
2386       END SELECT
2387
2388    END SUBROUTINE calc_pt_near_surface
2389   
2390
2391!
2392!-- Integrated stability function for momentum
2393    FUNCTION psi_m( zeta ) 
2394       
2395       USE kinds
2396
2397       IMPLICIT NONE
2398
2399       REAL(wp)            :: psi_m !< Integrated similarity function result
2400       REAL(wp)            :: zeta  !< Stability parameter z/L
2401       REAL(wp)            :: x     !< dummy variable
2402
2403       REAL(wp), PARAMETER :: a = 1.0_wp            !< constant
2404       REAL(wp), PARAMETER :: b = 0.66666666666_wp  !< constant
2405       REAL(wp), PARAMETER :: c = 5.0_wp            !< constant
2406       REAL(wp), PARAMETER :: d = 0.35_wp           !< constant
2407       REAL(wp), PARAMETER :: c_d_d = c / d         !< constant
2408       REAL(wp), PARAMETER :: bc_d_d = b * c / d    !< constant
2409
2410
2411       IF ( zeta < 0.0_wp )  THEN
2412          x = SQRT( SQRT( 1.0_wp  - 16.0_wp * zeta ) )
2413          psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2    &
2414                  * ( 1.0_wp + x**2 ) * 0.125_wp )
2415       ELSE
2416
2417          psi_m = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - a * zeta         &
2418                   - bc_d_d
2419!
2420!--       Old version for stable conditions (only valid for z/L < 0.5)
2421!--       psi_m = - 5.0_wp * zeta
2422
2423       ENDIF
2424
2425    END FUNCTION psi_m
2426
2427
2428!
2429!-- Integrated stability function for heat and moisture
2430    FUNCTION psi_h( zeta ) 
2431       
2432       USE kinds
2433
2434       IMPLICIT NONE
2435
2436       REAL(wp)            :: psi_h !< Integrated similarity function result
2437       REAL(wp)            :: zeta  !< Stability parameter z/L
2438       REAL(wp)            :: x     !< dummy variable
2439
2440       REAL(wp), PARAMETER :: a = 1.0_wp            !< constant
2441       REAL(wp), PARAMETER :: b = 0.66666666666_wp  !< constant
2442       REAL(wp), PARAMETER :: c = 5.0_wp            !< constant
2443       REAL(wp), PARAMETER :: d = 0.35_wp           !< constant
2444       REAL(wp), PARAMETER :: c_d_d = c / d         !< constant
2445       REAL(wp), PARAMETER :: bc_d_d = b * c / d    !< constant
2446
2447
2448       IF ( zeta < 0.0_wp )  THEN
2449          x = SQRT( 1.0_wp  - 16.0_wp * zeta )
2450          psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
2451       ELSE
2452          psi_h = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - (1.0_wp          &
2453                  + 0.66666666666_wp * a * zeta )**1.5_wp - bc_d_d             &
2454                  + 1.0_wp
2455!
2456!--       Old version for stable conditions (only valid for z/L < 0.5)
2457!--       psi_h = - 5.0_wp * zeta
2458       ENDIF
2459
2460    END FUNCTION psi_h
2461
2462
2463!------------------------------------------------------------------------------!
2464! Description:
2465! ------------
2466!> Calculates stability function for momentum
2467!>
2468!> @author Hauke Wurps
2469!------------------------------------------------------------------------------!
2470    FUNCTION phi_m( zeta ) 
2471   
2472       IMPLICIT NONE
2473   
2474       REAL(wp)            :: phi_m         !< Value of the function
2475       REAL(wp)            :: zeta          !< Stability parameter z/L
2476   
2477       REAL(wp), PARAMETER :: a = 16.0_wp   !< constant
2478       REAL(wp), PARAMETER :: c = 5.0_wp    !< constant
2479   
2480       IF ( zeta < 0.0_wp )  THEN
2481          phi_m = 1.0_wp / SQRT( SQRT( 1.0_wp - a * zeta ) )
2482       ELSE
2483          phi_m = 1.0_wp + c * zeta
2484       ENDIF
2485   
2486    END FUNCTION phi_m
2487
2488 END MODULE surface_layer_fluxes_mod
Note: See TracBrowser for help on using the repository browser.