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

Last change on this file since 3137 was 3130, checked in by gronemeier, 6 years ago

calculate km according to MOST within the surface layer

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