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

Last change on this file since 3156 was 3152, checked in by suehring, 6 years ago

Further adjustments for surface structure

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