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

Last change on this file since 4313 was 4298, checked in by suehring, 4 years ago

Bugfixes: Calculation of 2-m temperature adjusted to the case the 2-m level is above the first grid level; salsa: close netcdf input files after reading; open netcdf input files with read-only attribute instead of write attribute

  • Property svn:keywords set to Id
File size: 72.6 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-2019 Leibniz Universitaet Hannover
18!
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! ------------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: surface_layer_fluxes_mod.f90 4298 2019-11-21 15:59:16Z suehring $
28! Calculation of 2-m temperature adjusted to the case the 2-m level is above
29! the first grid point.
30!
31! 4258 2019-10-07 13:29:08Z suehring
32! Initialization of Obukhov lenght also at vertical surfaces (if allocated).
33!
34! 4237 2019-09-25 11:33:42Z knoop
35! Added missing OpenMP directives
36!
37! 4186 2019-08-23 16:06:14Z suehring
38! - To enable limitation of Obukhov length, move it before exit-cycle construct.
39!   Further, change the limit to 10E-5 in order to get rid-off unrealistic
40!   peaks in the heat fluxes during nighttime
41! - Unused variable removed
42!
43! 4182 2019-08-22 15:20:23Z scharf
44! Corrected "Former revisions" section
45!
46! 3987 2019-05-22 09:52:13Z kanani
47! Introduce alternative switch for debug output during timestepping
48!
49! 3885 2019-04-11 11:29:34Z kanani
50! Changes related to global restructuring of location messages and introduction
51! of additional debug messages
52!
53! 3881 2019-04-10 09:31:22Z suehring
54! Assure that Obukhov length does not become zero
55!
56! 3834 2019-03-28 15:40:15Z forkel
57! added USE chem_gasphase_mod
58!
59! 3787 2019-03-07 08:43:54Z raasch
60! unused variables removed
61!
62! 3745 2019-02-15 18:57:56Z suehring
63! Bugfix, missing calculation of 10cm temperature at vertical building walls,
64! required for indoor model
65!
66! 3744 2019-02-15 18:38:58Z suehring
67! Some interface calls moved to module_interface + cleanup
68!
69! 3668 2019-01-14 12:49:24Z maronga
70! Removed methods "circular" and "lookup"
71!
72! 3655 2019-01-07 16:51:22Z knoop
73! OpenACC port for SPEC
74!
75! Revision 1.1  1998/01/23 10:06:06  raasch
76! Initial revision
77!
78!
79! Description:
80! ------------
81!> Diagnostic computation of vertical fluxes in the constant flux layer from the
82!> values of the variables at grid point k=1 based on Newton iteration
83!>
84!> @todo (re)move large_scale_forcing actions
85!> @todo check/optimize OpenMP directives
86!> @todo simplify if conditions (which flux need to be computed in which case)
87!------------------------------------------------------------------------------!
88 MODULE surface_layer_fluxes_mod
89
90    USE arrays_3d,                                                             &
91        ONLY:  e, kh, nc, nr, pt, q, ql, qc, qr, s, u, v, vpt, w, zu, zw,      &
92               drho_air_zw, rho_air_zw, d_exner
93
94    USE basic_constants_and_equations_mod,                                     &
95        ONLY:  g, kappa, lv_d_cp, pi, rd_d_rv
96
97    USE chem_gasphase_mod,                                                     &
98        ONLY:  nvar
99
100    USE chem_modules,                                                          &
101        ONLY:  constant_csflux
102
103    USE cpulog
104
105    USE control_parameters,                                                    &
106        ONLY:  air_chemistry, cloud_droplets,                                  &
107               constant_heatflux, constant_scalarflux,                         &
108               constant_waterflux, coupling_mode,                              &
109               debug_output_timestep,                                          &
110               do_output_at_2m, humidity,                                      &
111               ibc_e_b, ibc_pt_b, indoor_model,                                &
112               land_surface, large_scale_forcing, lsf_surf, message_string,    &
113               neutral, passive_scalar, pt_surface, q_surface,                 &
114               run_coupled, surface_pressure, simulated_time,                  &
115               time_since_reference_point, urban_surface,                      &
116               use_free_convection_scaling, zeta_max, zeta_min
117
118    USE grid_variables,                                                        &
119        ONLY:  dx, dy 
120
121    USE indices,                                                               &
122        ONLY:  nzt
123
124    USE kinds
125
126    USE bulk_cloud_model_mod,                                                  &
127        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
128
129    USE pegrid
130
131    USE land_surface_model_mod,                                                &
132        ONLY:  aero_resist_kray, skip_time_do_lsm
133
134    USE surface_mod,                                                           &
135        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_type,     &
136                surf_usm_h, surf_usm_v
137       
138
139    IMPLICIT NONE
140
141    INTEGER(iwp) ::  i              !< loop index x direction
142    INTEGER(iwp) ::  j              !< loop index y direction
143    INTEGER(iwp) ::  k              !< loop index z direction
144    INTEGER(iwp) ::  l              !< loop index for surf type
145
146    LOGICAL      ::  coupled_run       !< Flag for coupled atmosphere-ocean runs
147    LOGICAL      ::  downward = .FALSE.!< Flag indicating downward-facing horizontal surface
148    LOGICAL      ::  mom_uv  = .FALSE. !< Flag indicating calculation of usvs and vsus at vertical surfaces
149    LOGICAL      ::  mom_w   = .FALSE. !< Flag indicating calculation of wsus and wsvs at vertical surfaces
150    LOGICAL      ::  mom_tke = .FALSE. !< Flag indicating calculation of momentum fluxes at vertical surfaces used for TKE production
151    LOGICAL      ::  surf_vertical     !< Flag indicating vertical surfaces
152
153    REAL(wp)     ::  e_s,               & !< Saturation water vapor pressure
154                     ol_max = 1.0E6_wp, & !< Maximum Obukhov length
155                     z_mo                 !< Height of the constant flux layer where MOST is assumed
156
157    TYPE(surf_type), POINTER ::  surf     !< surf-type array, used to generalize subroutines
158
159
160    SAVE
161
162    PRIVATE
163
164    PUBLIC init_surface_layer_fluxes, phi_m, surface_layer_fluxes
165
166    INTERFACE init_surface_layer_fluxes
167       MODULE PROCEDURE init_surface_layer_fluxes
168    END INTERFACE init_surface_layer_fluxes
169
170    INTERFACE phi_m
171       MODULE PROCEDURE phi_m
172    END INTERFACE phi_m
173
174    INTERFACE surface_layer_fluxes
175       MODULE PROCEDURE surface_layer_fluxes
176    END INTERFACE surface_layer_fluxes
177
178
179 CONTAINS
180
181
182!------------------------------------------------------------------------------!
183! Description:
184! ------------
185!> Main routine to compute the surface fluxes
186!------------------------------------------------------------------------------!
187    SUBROUTINE surface_layer_fluxes
188
189       IMPLICIT NONE
190
191
192       IF ( debug_output_timestep )  CALL debug_message( 'surface_layer_fluxes', 'start' )
193
194       surf_vertical = .FALSE. !< flag indicating vertically orientated surface elements
195       downward      = .FALSE. !< flag indicating downward-facing surface elements
196!
197!--    Derive potential temperature and specific humidity at first grid level
198!--    from the fields pt and q
199!
200!--    First call for horizontal default-type surfaces (l=0 - upward facing,
201!--    l=1 - downward facing)
202       DO  l = 0, 1
203          IF ( surf_def_h(l)%ns >= 1  )  THEN
204             surf => surf_def_h(l)
205             CALL calc_pt_q
206             IF ( .NOT. neutral )  THEN
207                CALL calc_pt_surface
208                IF ( humidity )  THEN
209                   CALL calc_q_surface
210                   CALL calc_vpt_surface
211                ENDIF
212             ENDIF
213          ENDIF
214       ENDDO
215!
216!--    Call for natural-type horizontal surfaces
217       IF ( surf_lsm_h%ns >= 1  )  THEN
218          surf => surf_lsm_h
219          CALL calc_pt_q
220       ENDIF
221
222!
223!--    Call for urban-type horizontal surfaces
224       IF ( surf_usm_h%ns >= 1  )  THEN
225          surf => surf_usm_h
226          CALL calc_pt_q
227       ENDIF
228
229!
230!--    Call for natural-type vertical surfaces
231       DO  l = 0, 3
232          IF ( surf_lsm_v(l)%ns >= 1  )  THEN
233             surf => surf_lsm_v(l)
234             CALL calc_pt_q
235          ENDIF
236
237!--       Call for urban-type vertical surfaces
238          IF ( surf_usm_v(l)%ns >= 1  )  THEN
239             surf => surf_usm_v(l)
240             CALL calc_pt_q
241          ENDIF
242       ENDDO
243
244!
245!--    First, calculate the new Obukhov length, then new friction velocity,
246!--    followed by the new scaling parameters (th*, q*, etc.), and the new
247!--    surface fluxes if required. Note, each routine is called for different surface types.
248!--    First call for default-type horizontal surfaces, for natural- and
249!--    urban-type surfaces. Note, at this place only upward-facing horizontal
250!--    surfaces are treated.
251
252!
253!--    Default-type upward-facing horizontal surfaces
254       IF ( surf_def_h(0)%ns >= 1  )  THEN
255          surf => surf_def_h(0)
256          CALL calc_uvw_abs
257          IF ( .NOT. neutral )  CALL calc_ol
258          CALL calc_us
259          CALL calc_scaling_parameters
260          CALL calc_surface_fluxes
261          IF ( do_output_at_2m )  THEN
262             CALL calc_pt_near_surface ( '2m' )
263          ENDIF
264       ENDIF
265!
266!--    Natural-type horizontal surfaces
267       IF ( surf_lsm_h%ns >= 1  )  THEN
268          surf => surf_lsm_h
269          CALL calc_uvw_abs
270          IF ( .NOT. neutral )  CALL calc_ol
271          CALL calc_us
272          CALL calc_scaling_parameters
273          CALL calc_surface_fluxes
274          IF ( do_output_at_2m )  THEN
275             CALL calc_pt_near_surface ( '2m' )
276          ENDIF
277       ENDIF
278!
279!--    Urban-type horizontal surfaces
280       IF ( surf_usm_h%ns >= 1  )  THEN
281          surf => surf_usm_h
282          CALL calc_uvw_abs
283          IF ( .NOT. neutral )  CALL calc_ol
284          CALL calc_us
285          CALL calc_scaling_parameters
286          CALL calc_surface_fluxes
287          IF ( do_output_at_2m )  THEN
288             CALL calc_pt_near_surface ( '2m' )
289          ENDIF
290!
291!--       Calculate 10cm temperature, required in indoor model
292          IF ( indoor_model )  CALL calc_pt_near_surface ( '10cm' )
293       ENDIF
294
295!
296!--    Treat downward-facing horizontal surfaces. Note, so far, these are
297!--    always default type. Stratification is not considered
298!--    in this case, hence, no further distinction between different
299!--    most_method is required. 
300       IF ( surf_def_h(1)%ns >= 1  )  THEN
301          downward = .TRUE.
302          surf => surf_def_h(1)
303          CALL calc_uvw_abs
304          CALL calc_us
305          CALL calc_surface_fluxes
306          downward = .FALSE.
307       ENDIF
308!
309!--    Calculate surfaces fluxes at vertical surfaces for momentum
310!--    and subgrid-scale TKE.
311!--    No stability is considered. Therefore, scaling parameters and Obukhov-
312!--    length do not need to be calculated and no distinction in 'circular',
313!--    'Newton' or 'lookup' is necessary so far.
314!--    Note, this will change if stability is once considered.
315       surf_vertical = .TRUE.
316!
317!--    Calculate horizontal momentum fluxes at north- and south-facing
318!--    surfaces(usvs).
319!--    For default-type surfaces
320       mom_uv = .TRUE. 
321       DO  l = 0, 1
322          IF ( surf_def_v(l)%ns >= 1  )  THEN
323             surf => surf_def_v(l)
324!
325!--          Compute surface-parallel velocity
326             CALL calc_uvw_abs_v_ugrid
327!
328!--          Compute respective friction velocity on staggered grid
329             CALL calc_us
330!
331!--          Compute respective surface fluxes for momentum and TKE
332             CALL calc_surface_fluxes
333          ENDIF
334       ENDDO
335!
336!--    For natural-type surfaces. Please note, even though stability is not
337!--    considered for the calculation of momentum fluxes at vertical surfaces,
338!--    scaling parameters and Obukov length are calculated nevertheless in this
339!--    case. This is due to the requirement of ts in parameterization of heat
340!--    flux in land-surface model in case of aero_resist_kray is not true.
341       IF ( .NOT. aero_resist_kray )  THEN
342          DO  l = 0, 1
343             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
344                surf => surf_lsm_v(l)
345!
346!--             Compute surface-parallel velocity
347                CALL calc_uvw_abs_v_ugrid
348!
349!--             Compute Obukhov length
350                IF ( .NOT. neutral )  CALL calc_ol
351!
352!--             Compute respective friction velocity on staggered grid
353                CALL calc_us
354!
355!--             Compute scaling parameters
356                CALL calc_scaling_parameters
357!
358!--             Compute respective surface fluxes for momentum and TKE
359                CALL calc_surface_fluxes
360             ENDIF
361          ENDDO
362!
363!--    No ts is required, so scaling parameters and Obukhov length do not need
364!--    to be computed.
365       ELSE
366          DO  l = 0, 1
367             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
368                surf => surf_lsm_v(l)
369!   
370!--             Compute surface-parallel velocity
371                CALL calc_uvw_abs_v_ugrid
372!
373!--             Compute respective friction velocity on staggered grid
374                CALL calc_us
375!
376!--             Compute respective surface fluxes for momentum and TKE
377                CALL calc_surface_fluxes
378             ENDIF
379          ENDDO
380       ENDIF
381!
382!--    For urban-type surfaces
383       DO  l = 0, 1
384          IF ( surf_usm_v(l)%ns >= 1  )  THEN
385             surf => surf_usm_v(l)
386!
387!--          Compute surface-parallel velocity
388             CALL calc_uvw_abs_v_ugrid
389!
390!--          Compute respective friction velocity on staggered grid
391             CALL calc_us
392!
393!--          Compute respective surface fluxes for momentum and TKE
394             CALL calc_surface_fluxes
395!
396!--          Calculate 10cm temperature, required in indoor model
397             IF ( indoor_model )  CALL calc_pt_near_surface ( '10cm' )
398          ENDIF
399       ENDDO
400!
401!--    Calculate horizontal momentum fluxes at east- and west-facing
402!--    surfaces (vsus).
403!--    For default-type surfaces
404       DO  l = 2, 3
405          IF ( surf_def_v(l)%ns >= 1  )  THEN
406             surf => surf_def_v(l)
407!
408!--          Compute surface-parallel velocity
409             CALL calc_uvw_abs_v_vgrid
410!
411!--          Compute respective friction velocity on staggered grid
412             CALL calc_us
413!
414!--          Compute respective surface fluxes for momentum and TKE
415             CALL calc_surface_fluxes
416             
417          ENDIF
418       ENDDO
419!
420!--    For natural-type surfaces. Please note, even though stability is not
421!--    considered for the calculation of momentum fluxes at vertical surfaces,
422!--    scaling parameters and Obukov length are calculated nevertheless in this
423!--    case. This is due to the requirement of ts in parameterization of heat
424!--    flux in land-surface model in case of aero_resist_kray is not true.
425       IF ( .NOT. aero_resist_kray )  THEN
426          DO  l = 2, 3
427             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
428                surf => surf_lsm_v(l)
429!
430!--             Compute surface-parallel velocity
431                CALL calc_uvw_abs_v_vgrid
432!
433!--             Compute Obukhov length
434                IF ( .NOT. neutral )  CALL calc_ol
435!
436!--             Compute respective friction velocity on staggered grid
437                CALL calc_us
438!
439!--             Compute scaling parameters
440                CALL calc_scaling_parameters
441!
442!--             Compute respective surface fluxes for momentum and TKE
443                CALL calc_surface_fluxes
444             ENDIF
445          ENDDO
446       ELSE
447          DO  l = 2, 3
448             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
449                surf => surf_lsm_v(l)
450!   
451!--             Compute surface-parallel velocity
452                CALL calc_uvw_abs_v_vgrid
453!
454!--             Compute respective friction velocity on staggered grid
455                CALL calc_us
456!
457!--             Compute respective surface fluxes for momentum and TKE
458                CALL calc_surface_fluxes
459             ENDIF
460          ENDDO
461       ENDIF
462!
463!--    For urban-type surfaces
464       DO  l = 2, 3
465          IF ( surf_usm_v(l)%ns >= 1  )  THEN
466             surf => surf_usm_v(l)
467!
468!--          Compute surface-parallel velocity
469             CALL calc_uvw_abs_v_vgrid
470!
471!--          Compute respective friction velocity on staggered grid
472             CALL calc_us
473!
474!--          Compute respective surface fluxes for momentum and TKE
475             CALL calc_surface_fluxes
476!
477!--          Calculate 10cm temperature, required in indoor model             
478             IF ( indoor_model )  CALL calc_pt_near_surface ( '10cm' )
479          ENDIF
480       ENDDO
481       mom_uv = .FALSE.
482!
483!--    Calculate horizontal momentum fluxes of w (wsus and wsvs) at vertial
484!--    surfaces.
485       mom_w = .TRUE.
486!
487!--    Default-type surfaces
488       DO  l = 0, 3
489          IF ( surf_def_v(l)%ns >= 1  )  THEN
490             surf => surf_def_v(l)
491             CALL calc_uvw_abs_v_wgrid
492             CALL calc_us
493             CALL calc_surface_fluxes
494          ENDIF
495       ENDDO 
496!
497!--    Natural-type surfaces
498       DO  l = 0, 3
499          IF ( surf_lsm_v(l)%ns >= 1  )  THEN
500             surf => surf_lsm_v(l)
501             CALL calc_uvw_abs_v_wgrid
502             CALL calc_us
503             CALL calc_surface_fluxes
504          ENDIF
505       ENDDO 
506!
507!--    Urban-type surfaces
508       DO  l = 0, 3
509          IF ( surf_usm_v(l)%ns >= 1  )  THEN
510             surf => surf_usm_v(l)
511             CALL calc_uvw_abs_v_wgrid
512             CALL calc_us
513             CALL calc_surface_fluxes
514          ENDIF
515       ENDDO 
516       mom_w = .FALSE.   
517!
518!--    Calculate momentum fluxes usvs, vsus, wsus and wsvs at vertical
519!--    surfaces for TKE production. Note, here, momentum fluxes are defined
520!--    at grid center and are not staggered as before.
521       mom_tke = .TRUE.
522!
523!--    Default-type surfaces
524       DO  l = 0, 3
525          IF ( surf_def_v(l)%ns >= 1  )  THEN
526             surf => surf_def_v(l)
527             CALL calc_uvw_abs_v_sgrid
528             CALL calc_us
529             CALL calc_surface_fluxes
530          ENDIF
531       ENDDO 
532!
533!--    Natural-type surfaces
534       DO  l = 0, 3
535          IF ( surf_lsm_v(l)%ns >= 1  )  THEN
536             surf => surf_lsm_v(l)
537             CALL calc_uvw_abs_v_sgrid
538             CALL calc_us
539             CALL calc_surface_fluxes
540          ENDIF
541       ENDDO 
542!
543!--    Urban-type surfaces
544       DO  l = 0, 3
545          IF ( surf_usm_v(l)%ns >= 1  )  THEN
546             surf => surf_usm_v(l)
547             CALL calc_uvw_abs_v_sgrid
548             CALL calc_us
549             CALL calc_surface_fluxes
550          ENDIF
551       ENDDO 
552       mom_tke = .FALSE.
553
554       IF ( debug_output_timestep )  CALL debug_message( 'surface_layer_fluxes', 'end' )
555
556    END SUBROUTINE surface_layer_fluxes
557
558
559!------------------------------------------------------------------------------!
560! Description:
561! ------------
562!> Initializing actions for the surface layer routine.
563!------------------------------------------------------------------------------!
564    SUBROUTINE init_surface_layer_fluxes
565
566       IMPLICIT NONE
567
568       INTEGER(iwp) ::  l  !< running index for vertical surface orientation
569
570       CALL location_message( 'initializing surface layer', 'start' )
571
572!
573!--    In case of runs with neutral statification, set Obukhov length to a
574!--    large value
575       IF ( neutral )  THEN
576          IF ( surf_def_h(0)%ns >= 1 )  surf_def_h(0)%ol = 1.0E10_wp
577          IF ( surf_lsm_h%ns    >= 1 )  surf_lsm_h%ol    = 1.0E10_wp
578          IF ( surf_usm_h%ns    >= 1 )  surf_usm_h%ol    = 1.0E10_wp
579         
580          DO  l = 0, 3
581             IF ( surf_def_v(l)%ns >= 1  .AND.                                 &
582                  ALLOCATED( surf_def_v(l)%ol ) )  surf_def_v(l)%ol = 1.0E10_wp
583             IF ( surf_lsm_v(l)%ns >= 1  .AND.                                 &
584                  ALLOCATED( surf_lsm_v(l)%ol ) )  surf_lsm_v(l)%ol = 1.0E10_wp
585             IF ( surf_usm_v(l)%ns >= 1  .AND.                                 &
586                  ALLOCATED( surf_usm_v(l)%ol ) )  surf_usm_v(l)%ol = 1.0E10_wp 
587          ENDDO
588         
589       ENDIF
590
591       CALL location_message( 'initializing surface layer', 'finished' )
592
593    END SUBROUTINE init_surface_layer_fluxes
594
595
596!------------------------------------------------------------------------------!
597! Description:
598! ------------
599!> Compute the absolute value of the horizontal velocity (relative to the
600!> surface) for horizontal surface elements. This is required by all methods.
601!------------------------------------------------------------------------------!
602    SUBROUTINE calc_uvw_abs
603   
604       IMPLICIT NONE
605
606       INTEGER(iwp) ::  i             !< running index x direction
607       INTEGER(iwp) ::  ibit          !< flag to mask computation of relative velocity in case of downward-facing surfaces
608       INTEGER(iwp) ::  j             !< running index y direction
609       INTEGER(iwp) ::  k             !< running index z direction
610       INTEGER(iwp) ::  m             !< running index surface elements
611
612       REAL(wp)     :: w_lfc          !< local free convection velocity scale
613!
614!--    ibit is 1 for upward-facing surfaces, zero for downward-facing surfaces.
615       ibit = MERGE( 1, 0, .NOT. downward )
616
617       !$OMP PARALLEL DO PRIVATE(i, j, k, w_lfc)
618       !$ACC PARALLEL LOOP PRIVATE(i, j, k, w_lfc) &
619       !$ACC PRESENT(surf, u, v)
620       DO  m = 1, surf%ns
621
622          i   = surf%i(m)           
623          j   = surf%j(m)
624          k   = surf%k(m)
625       
626!
627!--       Calculate free convection velocity scale w_lfc is
628!--       use_free_convection_scaling = .T.. This will maintain a horizontal
629!--       velocity even for very weak wind convective conditions. SIGN is used
630!--       to set w_lfc to zero under stable conditions.
631          IF ( use_free_convection_scaling )  THEN
632             w_lfc = ABS(g / surf%pt1(m) * surf%z_mo(m) * surf%shf(m))
633             w_lfc = ( 0.5_wp * ( w_lfc + SIGN(w_lfc,surf%shf(m)) ) )**(0.33333_wp)
634          ELSE
635             w_lfc = 0.0_wp
636          ENDIF
637
638!
639!--       Compute the absolute value of the horizontal velocity.
640!--       (relative to the surface in case the lower surface is the ocean).
641!--       Please note, in new surface modelling concept the index values changed,
642!--       i.e. the reference grid point is not the surface-grid point itself but
643!--       the first grid point outside of the topography.
644!--       Note, in case of coupled ocean-atmosphere simulations relative velocity
645!--       with respect to the ocean surface is used, hence, (k-1,j,i) values
646!--       are used to calculate the absolute velocity.
647!--       However, this do not apply for downward-facing walls. To mask this,
648!--       use ibit, which checks for upward/downward-facing surfaces.
649          surf%uvw_abs(m) = SQRT(                                              &
650                              ( 0.5_wp * (   u(k,j,i)   + u(k,j,i+1)           &
651                                        -  ( u(k-1,j,i) + u(k-1,j,i+1)         &
652                                           ) * ibit                            &
653                                         )                                     &
654                              )**2 +                                           &
655                              ( 0.5_wp * (   v(k,j,i)   + v(k,j+1,i)           &
656                                        -  ( v(k-1,j,i) + v(k-1,j+1,i)         &
657                                           ) * ibit                            &
658                                         )                                     &
659                              )**2  + w_lfc**2                                 &
660                                )
661
662                               
663
664       ENDDO
665
666    END SUBROUTINE calc_uvw_abs
667
668
669!------------------------------------------------------------------------------!
670! Description:
671! ------------
672!> Compute the absolute value of the horizontal velocity (relative to the
673!> surface) for horizontal surface elements. This is required by all methods.
674!------------------------------------------------------------------------------!
675    SUBROUTINE calc_uvw_abs_v_ugrid
676
677       IMPLICIT NONE
678
679       INTEGER(iwp) ::  i   !< running index x direction
680       INTEGER(iwp) ::  j   !< running index y direction
681       INTEGER(iwp) ::  k   !< running index z direction
682       INTEGER(iwp) ::  m   !< running index surface elements
683
684       REAL(wp)     ::  u_i !< u-component on xu-grid
685       REAL(wp)     ::  w_i !< w-component on xu-grid
686
687
688       DO  m = 1, surf%ns
689          i   = surf%i(m)           
690          j   = surf%j(m)
691          k   = surf%k(m)
692!
693!--       Compute the absolute value of the surface parallel velocity on u-grid.
694          u_i  = u(k,j,i)
695          w_i  = 0.25_wp * ( w(k-1,j,i-1) + w(k-1,j,i) +                       &
696                             w(k,j,i-1)   + w(k,j,i) ) 
697
698          surf%uvw_abs(m) = SQRT( u_i**2 + w_i**2 ) 
699
700       ENDDO
701
702    END SUBROUTINE calc_uvw_abs_v_ugrid
703
704!------------------------------------------------------------------------------!
705! Description:
706! ------------
707!> Compute the absolute value of the horizontal velocity (relative to the
708!> surface) for horizontal surface elements. This is required by all methods.
709!------------------------------------------------------------------------------!
710    SUBROUTINE calc_uvw_abs_v_vgrid
711
712       IMPLICIT NONE
713
714       INTEGER(iwp) ::  i   !< running index x direction
715       INTEGER(iwp) ::  j   !< running index y direction
716       INTEGER(iwp) ::  k   !< running index z direction
717       INTEGER(iwp) ::  m   !< running index surface elements
718
719       REAL(wp)     ::  v_i !< v-component on yv-grid
720       REAL(wp)     ::  w_i !< w-component on yv-grid
721
722
723       DO  m = 1, surf%ns
724          i   = surf%i(m)           
725          j   = surf%j(m)
726          k   = surf%k(m)
727
728          v_i  = u(k,j,i)
729          w_i  = 0.25_wp * ( w(k-1,j-1,i) + w(k-1,j,i) +                       &
730                             w(k,j-1,i)   + w(k,j,i) ) 
731
732          surf%uvw_abs(m) = SQRT( v_i**2 + w_i**2 ) 
733
734       ENDDO
735
736    END SUBROUTINE calc_uvw_abs_v_vgrid
737
738!------------------------------------------------------------------------------!
739! Description:
740! ------------
741!> Compute the absolute value of the horizontal velocity (relative to the
742!> surface) for horizontal surface elements. This is required by all methods.
743!------------------------------------------------------------------------------!
744    SUBROUTINE calc_uvw_abs_v_wgrid
745
746       IMPLICIT NONE
747
748       INTEGER(iwp) ::  i   !< running index x direction
749       INTEGER(iwp) ::  j   !< running index y direction
750       INTEGER(iwp) ::  k   !< running index z direction
751       INTEGER(iwp) ::  m   !< running index surface elements
752
753       REAL(wp)     ::  u_i !< u-component on x-zw-grid
754       REAL(wp)     ::  v_i !< v-component on y-zw-grid
755       REAL(wp)     ::  w_i !< w-component on zw-grid
756!
757!--    North- (l=0) and south-facing (l=1) surfaces
758       IF ( l == 0  .OR.  l == 1 )  THEN
759          DO  m = 1, surf%ns
760             i   = surf%i(m)           
761             j   = surf%j(m)
762             k   = surf%k(m)
763
764             u_i  = 0.25_wp * ( u(k+1,j,i+1) + u(k+1,j,i) +                    &
765                                u(k,j,i+1)   + u(k,j,i) )
766             v_i  = 0.0_wp
767             w_i  = w(k,j,i)
768
769             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
770          ENDDO
771!
772!--    East- (l=2) and west-facing (l=3) surfaces
773       ELSE
774          DO  m = 1, surf%ns
775             i   = surf%i(m)           
776             j   = surf%j(m)
777             k   = surf%k(m)
778
779             u_i  = 0.0_wp
780             v_i  = 0.25_wp * ( v(k+1,j+1,i) + v(k+1,j,i) +                    &
781                                v(k,j+1,i)   + v(k,j,i) )
782             w_i  = w(k,j,i)
783
784             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
785          ENDDO
786       ENDIF           
787
788    END SUBROUTINE calc_uvw_abs_v_wgrid
789
790!------------------------------------------------------------------------------!
791! Description:
792! ------------
793!> Compute the absolute value of the horizontal velocity (relative to the
794!> surface) for horizontal surface elements. This is required by all methods.
795!------------------------------------------------------------------------------!
796    SUBROUTINE calc_uvw_abs_v_sgrid
797
798       IMPLICIT NONE
799
800       INTEGER(iwp) ::  i   !< running index x direction
801       INTEGER(iwp) ::  j   !< running index y direction
802       INTEGER(iwp) ::  k   !< running index z direction
803       INTEGER(iwp) ::  m   !< running index surface elements
804
805       REAL(wp)     ::  u_i !< u-component on scalar grid
806       REAL(wp)     ::  v_i !< v-component on scalar grid
807       REAL(wp)     ::  w_i !< w-component on scalar grid
808
809!
810!--    North- (l=0) and south-facing (l=1) walls
811       IF ( l == 0  .OR.  l == 1 )  THEN
812          DO  m = 1, surf%ns
813             i   = surf%i(m)           
814             j   = surf%j(m)
815             k   = surf%k(m)
816
817             u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
818             v_i = 0.0_wp
819             w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
820
821             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
822          ENDDO
823!
824!--    East- (l=2) and west-facing (l=3) walls
825       ELSE
826          DO  m = 1, surf%ns
827             i   = surf%i(m)           
828             j   = surf%j(m)
829             k   = surf%k(m)
830
831             u_i = 0.0_wp
832             v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
833             w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
834
835             surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 
836          ENDDO
837       ENDIF 
838
839    END SUBROUTINE calc_uvw_abs_v_sgrid
840
841
842!------------------------------------------------------------------------------!
843! Description:
844! ------------
845!> Calculate the Obukhov length (L) and Richardson flux number (z/L)
846!------------------------------------------------------------------------------!
847    SUBROUTINE calc_ol
848
849       IMPLICIT NONE
850
851       INTEGER(iwp) ::  iter    !< Newton iteration step
852       INTEGER(iwp) ::  m       !< loop variable over all horizontal wall elements
853
854       REAL(wp)     :: f,      & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
855                       f_d_ol, & !< Derivative of f
856                       ol_l,   & !< Lower bound of L for Newton iteration
857                       ol_m,   & !< Previous value of L for Newton iteration
858                       ol_old, & !< Previous time step value of L
859                       ol_u      !< Upper bound of L for Newton iteration
860
861!
862!--    Evaluate bulk Richardson number (calculation depends on
863!--    definition based on setting of boundary conditions
864       IF ( ibc_pt_b /= 1 )  THEN
865          IF ( humidity )  THEN
866             !$OMP PARALLEL DO PRIVATE( z_mo )
867             DO  m = 1, surf%ns
868
869                z_mo = surf%z_mo(m)
870
871                surf%rib(m) = g * z_mo                                         &
872                              * ( surf%vpt1(m) - surf%vpt_surface(m) )         &
873                              / ( surf%uvw_abs(m)**2 * surf%vpt1(m)            &
874                              + 1.0E-20_wp )
875             ENDDO
876          ELSE
877             !$OMP PARALLEL DO PRIVATE( z_mo )
878             DO  m = 1, surf%ns
879
880                z_mo = surf%z_mo(m)
881
882                surf%rib(m) = g * z_mo                                         &
883                              * ( surf%pt1(m) - surf%pt_surface(m) )           &
884                              / ( surf%uvw_abs(m)**2 * surf%pt1(m) + 1.0E-20_wp )
885             ENDDO
886          ENDIF
887       ELSE
888          IF ( humidity )  THEN
889             !$OMP PARALLEL DO PRIVATE( k, z_mo )
890             DO  m = 1, surf%ns
891
892                k   = surf%k(m)
893
894                z_mo = surf%z_mo(m)
895
896                surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp                &
897                           * surf%qv1(m) ) * surf%shf(m) + 0.61_wp             &
898                           * surf%pt1(m) * surf%qsws(m) ) *                    &
899                             drho_air_zw(k-1)                /                 &
900                         ( surf%uvw_abs(m)**3 * surf%vpt1(m) * kappa**2        &
901                           + 1.0E-20_wp )
902             ENDDO
903          ELSE
904             !$OMP PARALLEL DO PRIVATE( k, z_mo )
905             !$ACC PARALLEL LOOP PRIVATE(k, z_mo) &
906             !$ACC PRESENT(surf, drho_air_zw)
907             DO  m = 1, surf%ns
908
909                k   = surf%k(m)
910
911                z_mo = surf%z_mo(m)
912
913                surf%rib(m) = - g * z_mo * surf%shf(m) *                    &
914                                     drho_air_zw(k-1)            /          &
915                        ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2       &
916                        + 1.0E-20_wp )
917             ENDDO
918          ENDIF
919       ENDIF
920
921!
922!--    Calculate the Obukhov length using Newton iteration
923       !$OMP PARALLEL DO PRIVATE(i, j, z_mo) &
924       !$OMP PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol)
925       !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
926       !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
927       !$ACC PRESENT(surf)
928       DO  m = 1, surf%ns
929
930          i   = surf%i(m)           
931          j   = surf%j(m)
932
933          z_mo = surf%z_mo(m)
934
935!
936!--       Store current value in case the Newton iteration fails
937          ol_old = surf%ol(m)
938
939!
940!--       Ensure that the bulk Richardson number and the Obukhov
941!--       length have the same sign
942          IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.                      &
943               ABS( surf%ol(m) ) == ol_max )  THEN
944             IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
945             IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
946          ENDIF
947!
948!--       Iteration to find Obukhov length
949          iter = 0
950          DO
951             iter = iter + 1
952!
953!--          In case of divergence, use the value of the previous time step
954             IF ( iter > 1000 )  THEN
955                surf%ol(m) = ol_old
956                EXIT
957             ENDIF
958
959             ol_m = surf%ol(m)
960             ol_l = ol_m - 0.001_wp * ol_m
961             ol_u = ol_m + 0.001_wp * ol_m
962
963
964             IF ( ibc_pt_b /= 1 )  THEN
965!
966!--             Calculate f = Ri - [...]/[...]^2 = 0
967                f = surf%rib(m) - ( z_mo / ol_m ) * (                          &
968                                              LOG( z_mo / surf%z0h(m) )        &
969                                              - psi_h( z_mo / ol_m )           &
970                                              + psi_h( surf%z0h(m) /           &
971                                                       ol_m )                  & 
972                                                     )                         &
973                                           / ( LOG( z_mo / surf%z0(m) )        &
974                                              - psi_m( z_mo / ol_m )           &
975                                              + psi_m( surf%z0(m) /  ol_m )    &
976                                                 )**2
977
978!
979!--              Calculate df/dL
980                 f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo /                  &
981                                                          surf%z0h(m) )        &
982                                         - psi_h( z_mo / ol_u )                &
983                                         + psi_h( surf%z0h(m) / ol_u )         &
984                                           )                                   &
985                                         / ( LOG( z_mo / surf%z0(m) )          &
986                                         - psi_m( z_mo / ol_u )                &
987                                         + psi_m( surf%z0(m) / ol_u )          &
988                                           )**2                                &
989                        + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) )        &
990                                         - psi_h( z_mo / ol_l )                &
991                                         + psi_h( surf%z0h(m) / ol_l )         &
992                                           )                                   &
993                                         / ( LOG( z_mo / surf%z0(m) )          &
994                                         - psi_m( z_mo / ol_l )                &
995                                         + psi_m( surf%z0(m) / ol_l )          &
996                                           )**2                                &
997                          ) / ( ol_u - ol_l )
998             ELSE
999!
1000!--             Calculate f = Ri - 1 /[...]^3 = 0
1001                f = surf%rib(m) - ( z_mo / ol_m ) /                            &
1002                                             ( LOG( z_mo / surf%z0(m) )        &
1003                                         - psi_m( z_mo / ol_m )                &
1004                                         + psi_m( surf%z0(m) / ol_m )          &
1005                                             )**3
1006
1007!
1008!--             Calculate df/dL
1009                f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / surf%z0(m) )      &
1010                                         - psi_m( z_mo / ol_u )                &
1011                                         + psi_m( surf%z0(m) / ol_u )          &
1012                                                  )**3                         &
1013                           + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )      &
1014                                         - psi_m( z_mo / ol_l )                &
1015                                         + psi_m( surf%z0(m) / ol_l )          &
1016                                            )**3                               &
1017                          ) / ( ol_u - ol_l )
1018             ENDIF
1019!
1020!--          Calculate new L
1021             surf%ol(m) = ol_m - f / f_d_ol
1022
1023!
1024!--          Ensure that the bulk Richardson number and the Obukhov
1025!--          length have the same sign and ensure convergence.
1026             IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
1027!
1028!--          If unrealistic value occurs, set L to the maximum
1029!--          value that is allowed
1030             IF ( ABS( surf%ol(m) ) > ol_max )  THEN
1031                surf%ol(m) = ol_max
1032                EXIT
1033             ENDIF
1034!
1035!--          Assure that Obukhov length does not become zero. If the limit is
1036!--          reached, exit the loop.
1037             IF ( ABS( surf%ol(m) ) < 1E-5_wp )  THEN
1038                surf%ol(m) = SIGN( 1E-5_wp, surf%ol(m) )
1039                EXIT
1040             ENDIF
1041!
1042!--          Check for convergence
1043             IF ( ABS( ( surf%ol(m) - ol_m ) /  surf%ol(m) ) < 1.0E-4_wp )  THEN
1044                EXIT
1045             ELSE
1046                CYCLE
1047             ENDIF             
1048
1049          ENDDO
1050       ENDDO
1051
1052    END SUBROUTINE calc_ol
1053
1054!
1055!-- Calculate friction velocity u*
1056    SUBROUTINE calc_us
1057
1058       IMPLICIT NONE
1059
1060       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1061
1062!
1063!--    Compute u* at horizontal surfaces at the scalars' grid points
1064       IF ( .NOT. surf_vertical )  THEN
1065!
1066!--       Compute u* at upward-facing surfaces
1067          IF ( .NOT. downward )  THEN
1068             !$OMP PARALLEL  DO PRIVATE( z_mo )
1069             !$ACC PARALLEL LOOP PRIVATE(z_mo) &
1070             !$ACC PRESENT(surf)
1071             DO  m = 1, surf%ns
1072
1073                z_mo = surf%z_mo(m)
1074!
1075!--             Compute u* at the scalars' grid points
1076                surf%us(m) = kappa * surf%uvw_abs(m) /                         &
1077                            ( LOG( z_mo / surf%z0(m) )                         &
1078                           - psi_m( z_mo / surf%ol(m) )                        &
1079                           + psi_m( surf%z0(m) / surf%ol(m) ) )
1080   
1081             ENDDO
1082!
1083!--       Compute u* at downward-facing surfaces. This case, do not consider
1084!--       any stability.
1085          ELSE
1086             !$OMP PARALLEL  DO PRIVATE( z_mo )
1087             !$ACC PARALLEL LOOP PRIVATE(z_mo) &
1088             !$ACC PRESENT(surf)
1089             DO  m = 1, surf%ns
1090
1091                z_mo = surf%z_mo(m)
1092!
1093!--             Compute u* at the scalars' grid points
1094                surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
1095   
1096             ENDDO
1097          ENDIF
1098!
1099!--    Compute u* at vertical surfaces at the u/v/v grid, respectively.
1100!--    No stability is considered in this case.
1101       ELSE
1102          !$OMP PARALLEL DO PRIVATE( z_mo )
1103          !$ACC PARALLEL LOOP PRIVATE(z_mo) &
1104          !$ACC PRESENT(surf)
1105          DO  m = 1, surf%ns
1106             z_mo = surf%z_mo(m)
1107
1108             surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
1109          ENDDO
1110       ENDIF
1111
1112    END SUBROUTINE calc_us
1113
1114!
1115!-- Calculate potential temperature, specific humidity, and virtual potential
1116!-- temperature at first grid level
1117    SUBROUTINE calc_pt_q
1118
1119       IMPLICIT NONE
1120
1121       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1122
1123       !$OMP PARALLEL DO PRIVATE( i, j, k )
1124       !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
1125       !$ACC PRESENT(surf, pt)
1126       DO  m = 1, surf%ns 
1127
1128          i   = surf%i(m)           
1129          j   = surf%j(m)
1130          k   = surf%k(m)
1131
1132#ifndef _OPENACC
1133          IF ( bulk_cloud_model ) THEN
1134             surf%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
1135             surf%qv1(m) = q(k,j,i) - ql(k,j,i)
1136          ELSEIF( cloud_droplets ) THEN
1137             surf%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
1138             surf%qv1(m) = q(k,j,i) 
1139          ELSE
1140#endif
1141             surf%pt1(m) = pt(k,j,i)
1142#ifndef _OPENACC
1143             IF ( humidity )  THEN
1144                surf%qv1(m) = q(k,j,i)
1145             ELSE
1146#endif
1147                surf%qv1(m) = 0.0_wp
1148#ifndef _OPENACC
1149             ENDIF
1150          ENDIF
1151
1152          IF ( humidity )  THEN
1153             surf%vpt1(m) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
1154          ENDIF
1155#endif
1156         
1157       ENDDO
1158
1159    END SUBROUTINE calc_pt_q
1160
1161
1162!
1163!-- Set potential temperature at surface grid level.
1164!-- ( only for upward-facing surfs )
1165    SUBROUTINE calc_pt_surface
1166
1167       IMPLICIT NONE
1168
1169       INTEGER(iwp) ::  k_off    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
1170       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1171       
1172       k_off = surf%koff
1173       !$OMP PARALLEL DO PRIVATE( i, j, k )
1174       !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
1175       !$ACC PRESENT(surf, pt)
1176       DO  m = 1, surf%ns 
1177
1178          i   = surf%i(m)           
1179          j   = surf%j(m)
1180          k   = surf%k(m)
1181
1182          surf%pt_surface(m) = pt(k+k_off,j,i)
1183
1184       ENDDO
1185
1186    END SUBROUTINE calc_pt_surface
1187
1188!
1189!-- Set mixing ratio at surface grid level. ( Only for upward-facing surfs. )
1190    SUBROUTINE calc_q_surface
1191
1192       IMPLICIT NONE
1193
1194       INTEGER(iwp) ::  k_off   !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
1195       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1196       
1197       k_off = surf%koff
1198       !$OMP PARALLEL DO PRIVATE( i, j, k )
1199       DO  m = 1, surf%ns 
1200
1201          i   = surf%i(m)           
1202          j   = surf%j(m)
1203          k   = surf%k(m)
1204
1205          surf%q_surface(m) = q(k+k_off,j,i)
1206
1207       ENDDO
1208
1209    END SUBROUTINE calc_q_surface
1210   
1211!
1212!-- Set virtual potential temperature at surface grid level.
1213!-- ( only for upward-facing surfs )
1214    SUBROUTINE calc_vpt_surface
1215
1216       IMPLICIT NONE
1217
1218       INTEGER(iwp) ::  k_off    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
1219       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
1220       
1221       k_off = surf%koff
1222       !$OMP PARALLEL DO PRIVATE( i, j, k )
1223       DO  m = 1, surf%ns 
1224
1225          i   = surf%i(m)           
1226          j   = surf%j(m)
1227          k   = surf%k(m)
1228
1229          surf%vpt_surface(m) = vpt(k+k_off,j,i)
1230
1231       ENDDO
1232
1233    END SUBROUTINE calc_vpt_surface
1234   
1235!
1236!-- Calculate the other MOST scaling parameters theta*, q*, (qc*, qr*, nc*, nr*)
1237    SUBROUTINE calc_scaling_parameters
1238
1239       IMPLICIT NONE
1240
1241
1242       INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements
1243       INTEGER(iwp)  ::  lsp     !< running index for chemical species
1244!
1245!--    Compute theta* at horizontal surfaces
1246       IF ( constant_heatflux  .AND.  .NOT. surf_vertical )  THEN
1247!
1248!--       For a given heat flux in the surface layer:
1249
1250          !$OMP PARALLEL DO PRIVATE( i, j, k )
1251          !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
1252          !$ACC PRESENT(surf, drho_air_zw)
1253          DO  m = 1, surf%ns 
1254
1255             i   = surf%i(m)           
1256             j   = surf%j(m)
1257             k   = surf%k(m)
1258
1259             surf%ts(m) = -surf%shf(m) * drho_air_zw(k-1) /                    &
1260                                  ( surf%us(m) + 1E-30_wp )
1261
1262!
1263!--          ts must be limited, because otherwise overflow may occur in case
1264!--          of us=0 when computing ol further below
1265             IF ( surf%ts(m) < -1.05E5_wp )  surf%ts(m) = -1.0E5_wp
1266             IF ( surf%ts(m) >  1.0E5_wp  )  surf%ts(m) =  1.0E5_wp
1267
1268          ENDDO
1269
1270       ELSEIF ( .NOT. surf_vertical ) THEN
1271!
1272!--       For a given surface temperature:
1273          IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
1274
1275             !$OMP PARALLEL DO PRIVATE( i, j, k )
1276             DO  m = 1, surf%ns 
1277                i   = surf%i(m)           
1278                j   = surf%j(m)
1279                k   = surf%k(m)
1280
1281                pt(k-1,j,i) = pt_surface
1282             ENDDO
1283          ENDIF
1284
1285          !$OMP PARALLEL DO PRIVATE( z_mo )
1286          DO  m = 1, surf%ns   
1287
1288             z_mo = surf%z_mo(m)
1289
1290             surf%ts(m) = kappa * ( surf%pt1(m) - surf%pt_surface(m) )      &
1291                                  / ( LOG( z_mo / surf%z0h(m) )             &
1292                                      - psi_h( z_mo / surf%ol(m) )          &
1293                                      + psi_h( surf%z0h(m) / surf%ol(m) ) )
1294
1295          ENDDO
1296
1297       ENDIF
1298!
1299!--    Compute theta* at vertical surfaces. This is only required in case of
1300!--    land-surface model, in order to compute aerodynamical resistance.
1301       IF ( surf_vertical )  THEN
1302          !$OMP PARALLEL DO PRIVATE( i, j )
1303          DO  m = 1, surf%ns 
1304
1305             i          =  surf%i(m)           
1306             j          =  surf%j(m)
1307             surf%ts(m) = -surf%shf(m) / ( surf%us(m) + 1E-30_wp )
1308!
1309!--          ts must be limited, because otherwise overflow may occur in case
1310!--          of us=0 when computing ol further below
1311             IF ( surf%ts(m) < -1.05E5_wp )  surf%ts(m) = -1.0E5_wp
1312             IF ( surf%ts(m) >  1.0E5_wp  )  surf%ts(m) =  1.0E5_wp
1313
1314          ENDDO
1315       ENDIF
1316
1317!
1318!--    If required compute q* at horizontal surfaces
1319       IF ( humidity )  THEN
1320          IF ( constant_waterflux  .AND.  .NOT. surf_vertical )  THEN
1321!
1322!--          For a given water flux in the surface layer
1323             !$OMP PARALLEL DO PRIVATE( i, j, k )
1324             DO  m = 1, surf%ns 
1325
1326                i   = surf%i(m)           
1327                j   = surf%j(m)
1328                k   = surf%k(m)
1329                surf%qs(m) = -surf%qsws(m) * drho_air_zw(k-1) /                &
1330                                               ( surf%us(m) + 1E-30_wp )
1331
1332             ENDDO
1333
1334          ELSEIF ( .NOT. surf_vertical )  THEN
1335             coupled_run = ( coupling_mode == 'atmosphere_to_ocean'  .AND.     &
1336                             run_coupled )
1337
1338             IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
1339                !$OMP PARALLEL DO PRIVATE( i, j, k )
1340                DO  m = 1, surf%ns 
1341
1342                   i   = surf%i(m)           
1343                   j   = surf%j(m)
1344                   k   = surf%k(m)
1345                   q(k-1,j,i) = q_surface
1346                   
1347                ENDDO
1348             ENDIF
1349
1350!
1351!--          Assume saturation for atmosphere coupled to ocean (but not
1352!--          in case of precursor runs)
1353             IF ( coupled_run )  THEN
1354                !$OMP PARALLEL DO PRIVATE( i, j, k, e_s )
1355                DO  m = 1, surf%ns   
1356                   i   = surf%i(m)           
1357                   j   = surf%j(m)
1358                   k   = surf%k(m)
1359                   e_s = 6.1_wp * &
1360                              EXP( 0.07_wp * ( MIN(pt(k-1,j,i),pt(k,j,i))      &
1361                                               - 273.15_wp ) )
1362                   q(k-1,j,i) = rd_d_rv * e_s / ( surface_pressure - e_s )
1363                ENDDO
1364             ENDIF
1365
1366             IF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
1367               !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1368                DO  m = 1, surf%ns   
1369
1370                   i   = surf%i(m)           
1371                   j   = surf%j(m)
1372                   k   = surf%k(m)
1373   
1374                   z_mo = surf%z_mo(m)
1375
1376                   surf%qs(m) = kappa * ( surf%qv1(m) - surf%q_surface(m) )    &
1377                                        / ( LOG( z_mo / surf%z0q(m) )          &
1378                                            - psi_h( z_mo / surf%ol(m) )       &
1379                                            + psi_h( surf%z0q(m) /             &
1380                                                     surf%ol(m) ) )
1381                ENDDO
1382             ELSE
1383                !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1384                DO  m = 1, surf%ns   
1385
1386                   i   = surf%i(m)           
1387                   j   = surf%j(m)
1388                   k   = surf%k(m)
1389   
1390                   z_mo = surf%z_mo(m)
1391
1392                   surf%qs(m) = kappa * ( q(k,j,i) - q(k-1,j,i) )              &
1393                                        / ( LOG( z_mo / surf%z0q(m) )          &
1394                                            - psi_h( z_mo / surf%ol(m) )       &
1395                                            + psi_h( surf%z0q(m) /             &
1396                                                     surf%ol(m) ) )
1397                ENDDO
1398             ENDIF
1399          ENDIF
1400!
1401!--       Compute q* at vertical surfaces
1402          IF ( surf_vertical )  THEN
1403             !$OMP PARALLEL DO PRIVATE( i, j )
1404             DO  m = 1, surf%ns 
1405
1406                i          =  surf%i(m)           
1407                j          =  surf%j(m)
1408                surf%qs(m) = -surf%qsws(m) / ( surf%us(m) + 1E-30_wp )
1409
1410             ENDDO
1411          ENDIF
1412       ENDIF
1413       
1414!
1415!--    If required compute s*
1416       IF ( passive_scalar )  THEN
1417!
1418!--       At horizontal surfaces
1419          IF ( constant_scalarflux  .AND.  .NOT. surf_vertical )  THEN
1420!
1421!--          For a given scalar flux in the surface layer
1422             !$OMP PARALLEL DO PRIVATE( i, j )
1423             DO  m = 1, surf%ns 
1424                i   = surf%i(m)           
1425                j   = surf%j(m)
1426                surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
1427             ENDDO
1428          ENDIF
1429!
1430!--       At vertical surfaces
1431          IF ( surf_vertical )  THEN
1432             !$OMP PARALLEL DO PRIVATE( i, j )
1433             DO  m = 1, surf%ns 
1434                i          =  surf%i(m)           
1435                j          =  surf%j(m)
1436                surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
1437             ENDDO
1438          ENDIF
1439       ENDIF
1440
1441!
1442!--    If required compute cs* (chemical species)
1443       IF ( air_chemistry  )  THEN 
1444!
1445!--       At horizontal surfaces                             
1446          DO  lsp = 1, nvar
1447             IF ( constant_csflux(lsp)  .AND.  .NOT.  surf_vertical )  THEN
1448!--             For a given chemical species' flux in the surface layer
1449                !$OMP PARALLEL DO PRIVATE( i, j )
1450                DO  m = 1, surf%ns 
1451                   i   = surf%i(m)           
1452                   j   = surf%j(m)
1453                   surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
1454                ENDDO
1455             ENDIF
1456          ENDDO
1457!
1458!--       At vertical surfaces
1459          IF ( surf_vertical )  THEN
1460             DO  lsp = 1, nvar
1461                !$OMP PARALLEL DO PRIVATE( i, j )
1462                DO  m = 1, surf%ns 
1463                   i   = surf%i(m)           
1464                   j   = surf%j(m)
1465                   surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
1466                ENDDO
1467             ENDDO
1468          ENDIF
1469       ENDIF
1470
1471!
1472!--    If required compute qc* and nc*
1473       IF ( bulk_cloud_model  .AND.  microphysics_morrison  .AND.              &
1474            .NOT. surf_vertical )  THEN
1475          !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1476          DO  m = 1, surf%ns   
1477
1478             i    = surf%i(m)           
1479             j    = surf%j(m)
1480             k    = surf%k(m)
1481
1482             z_mo = surf%z_mo(m)
1483
1484             surf%qcs(m) = kappa * ( qc(k,j,i) - qc(k-1,j,i) )                 &
1485                                 / ( LOG( z_mo / surf%z0q(m) )                 &
1486                                 - psi_h( z_mo / surf%ol(m) )                  &
1487                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
1488
1489             surf%ncs(m) = kappa * ( nc(k,j,i) - nc(k-1,j,i) )                 &
1490                                 / ( LOG( z_mo / surf%z0q(m) )                 &
1491                                 - psi_h( z_mo / surf%ol(m) )                  &
1492                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
1493          ENDDO
1494
1495       ENDIF
1496
1497!
1498!--    If required compute qr* and nr*
1499       IF ( bulk_cloud_model  .AND.  microphysics_seifert  .AND.               &
1500            .NOT. surf_vertical )  THEN
1501          !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1502          DO  m = 1, surf%ns   
1503
1504             i    = surf%i(m)           
1505             j    = surf%j(m)
1506             k    = surf%k(m)
1507
1508             z_mo = surf%z_mo(m)
1509
1510             surf%qrs(m) = kappa * ( qr(k,j,i) - qr(k-1,j,i) )                 &
1511                                 / ( LOG( z_mo / surf%z0q(m) )                 &
1512                                 - psi_h( z_mo / surf%ol(m) )                  &
1513                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
1514
1515             surf%nrs(m) = kappa * ( nr(k,j,i) - nr(k-1,j,i) )                 &
1516                                 / ( LOG( z_mo / surf%z0q(m) )                 &
1517                                 - psi_h( z_mo / surf%ol(m) )                  &
1518                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
1519          ENDDO
1520
1521       ENDIF
1522
1523    END SUBROUTINE calc_scaling_parameters
1524
1525
1526
1527!
1528!-- Calculate surface fluxes usws, vsws, shf, qsws, (qcsws, qrsws, ncsws, nrsws)
1529    SUBROUTINE calc_surface_fluxes
1530
1531       IMPLICIT NONE
1532
1533       INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements
1534       INTEGER(iwp)  ::  lsp     !< running index for chemical species
1535
1536       REAL(wp)                            ::  dum     !< dummy to precalculate logarithm
1537       REAL(wp)                            ::  flag_u  !< flag indicating u-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
1538       REAL(wp)                            ::  flag_v  !< flag indicating v-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
1539       REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_i     !< u-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
1540       REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_i     !< v-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
1541       REAL(wp), DIMENSION(:), ALLOCATABLE ::  w_i     !< w-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
1542
1543!
1544!--    Calcuate surface fluxes at horizontal walls
1545       IF ( .NOT. surf_vertical )  THEN
1546!
1547!--       Compute u'w' for the total model domain at upward-facing surfaces.
1548!--       First compute the corresponding component of u* and square it.
1549          IF ( .NOT. downward )  THEN
1550             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1551             !$ACC PARALLEL LOOP PRIVATE(i, j, k, z_mo) &
1552             !$ACC PRESENT(surf, u, rho_air_zw)
1553             DO  m = 1, surf%ns 
1554   
1555                i = surf%i(m)           
1556                j = surf%j(m)
1557                k = surf%k(m)
1558
1559                z_mo = surf%z_mo(m)
1560
1561                surf%usws(m) = kappa * ( u(k,j,i) - u(k-1,j,i) )               &
1562                              / ( LOG( z_mo / surf%z0(m) )                     &
1563                                  - psi_m( z_mo / surf%ol(m) )                 &
1564                                  + psi_m( surf%z0(m) / surf%ol(m) ) )
1565!
1566!--             Please note, the computation of usws is not fully accurate. Actually
1567!--             a further interpolation of us onto the u-grid, where usws is defined,
1568!--             is required. However, this is not done as this would require several
1569!--             data transfers between 2D-grid and the surf-type.
1570!--             The impact of the missing interpolation is negligible as several
1571!--             tests had shown.
1572!--             Same also for ol. 
1573                surf%usws(m) = -surf%usws(m) * surf%us(m) * rho_air_zw(k-1)
1574
1575             ENDDO
1576!
1577!--       At downward-facing surfaces
1578          ELSE
1579             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1580             DO  m = 1, surf%ns 
1581   
1582                i = surf%i(m)           
1583                j = surf%j(m)
1584                k = surf%k(m)
1585
1586                z_mo = surf%z_mo(m)
1587
1588                surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) )
1589                surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k)
1590
1591             ENDDO     
1592          ENDIF
1593
1594!
1595!--       Compute v'w' for the total model domain.
1596!--       First compute the corresponding component of u* and square it.
1597!--       Upward-facing surfaces
1598          IF ( .NOT. downward )  THEN
1599             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1600             !$ACC PARALLEL LOOP PRIVATE(i, j, k, z_mo) &
1601             !$ACC PRESENT(surf, v, rho_air_zw)
1602             DO  m = 1, surf%ns 
1603                i = surf%i(m)           
1604                j = surf%j(m)
1605                k = surf%k(m)
1606
1607                z_mo = surf%z_mo(m)
1608
1609                surf%vsws(m) = kappa * ( v(k,j,i) - v(k-1,j,i) )               &
1610                           / ( LOG( z_mo / surf%z0(m) )                        &
1611                               - psi_m( z_mo / surf%ol(m) )                    &
1612                               + psi_m( surf%z0(m) / surf%ol(m) ) )
1613!
1614!--             Please note, the computation of vsws is not fully accurate. Actually
1615!--             a further interpolation of us onto the v-grid, where vsws is defined,
1616!--             is required. However, this is not done as this would require several
1617!--             data transfers between 2D-grid and the surf-type.
1618!--             The impact of the missing interpolation is negligible as several
1619!--             tests had shown.
1620!--             Same also for ol. 
1621                surf%vsws(m) = -surf%vsws(m) * surf%us(m) * rho_air_zw(k-1)
1622             ENDDO
1623!
1624!--       Downward-facing surfaces
1625          ELSE
1626             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1627             DO  m = 1, surf%ns 
1628                i = surf%i(m)           
1629                j = surf%j(m)
1630                k = surf%k(m)
1631
1632                z_mo = surf%z_mo(m)
1633
1634                surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) )
1635                surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k)
1636             ENDDO
1637          ENDIF
1638!
1639!--       Compute the vertical kinematic heat flux
1640          IF (  .NOT.  constant_heatflux  .AND.  ( ( time_since_reference_point&
1641               <=  skip_time_do_lsm  .AND. simulated_time > 0.0_wp ) .OR.      &
1642               .NOT.  land_surface )  .AND.  .NOT. urban_surface  .AND.        &
1643               .NOT. downward )  THEN
1644             !$OMP PARALLEL DO PRIVATE( i, j, k )
1645             DO  m = 1, surf%ns 
1646                i    = surf%i(m)           
1647                j    = surf%j(m)
1648                k    = surf%k(m)
1649                surf%shf(m) = -surf%ts(m) * surf%us(m) * rho_air_zw(k-1)
1650             ENDDO
1651          ENDIF
1652!
1653!--       Compute the vertical water flux
1654          IF (  .NOT.  constant_waterflux  .AND.  humidity  .AND.              &
1655               ( ( time_since_reference_point <= skip_time_do_lsm  .AND.       &
1656               simulated_time > 0.0_wp ) .OR.  .NOT.  land_surface  )  .AND.   &
1657               .NOT.  urban_surface  .AND.  .NOT. downward )  THEN
1658             !$OMP PARALLEL DO PRIVATE( i, j, k )
1659             DO  m = 1, surf%ns 
1660                i    = surf%i(m)           
1661                j    = surf%j(m)
1662                k    = surf%k(m)
1663                surf%qsws(m) = -surf%qs(m) * surf%us(m) * rho_air_zw(k-1)
1664             ENDDO
1665          ENDIF
1666!
1667!--       Compute the vertical scalar flux
1668          IF (  .NOT.  constant_scalarflux  .AND.  passive_scalar  .AND.       &
1669                .NOT.  downward )  THEN
1670             !$OMP PARALLEL DO PRIVATE( i, j )
1671             DO  m = 1, surf%ns   
1672
1673                i    = surf%i(m)           
1674                j    = surf%j(m)
1675                surf%ssws(m) = -surf%ss(m) * surf%us(m)
1676
1677             ENDDO
1678          ENDIF   
1679!
1680!--       Compute the vertical chemical species' flux
1681          DO  lsp = 1, nvar
1682             IF (  .NOT.  constant_csflux(lsp)  .AND.  air_chemistry  .AND.    &
1683                   .NOT.  downward )  THEN
1684                !$OMP PARALLEL DO PRIVATE( i, j )
1685                DO  m = 1, surf%ns   
1686                   i     = surf%i(m)           
1687                   j     = surf%j(m)
1688                   surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m)
1689                ENDDO
1690             ENDIF   
1691          ENDDO
1692
1693!
1694!--       Compute (turbulent) fluxes of cloud water content and cloud drop conc.
1695          IF ( bulk_cloud_model  .AND.  microphysics_morrison  .AND.           &
1696               .NOT. downward)  THEN
1697             !$OMP PARALLEL DO PRIVATE( i, j )
1698             DO  m = 1, surf%ns   
1699
1700                i    = surf%i(m)           
1701                j    = surf%j(m)
1702
1703                surf%qcsws(m) = -surf%qcs(m) * surf%us(m)
1704                surf%ncsws(m) = -surf%ncs(m) * surf%us(m)
1705             ENDDO
1706          ENDIF   
1707!
1708!--       Compute (turbulent) fluxes of rain water content and rain drop conc.
1709          IF ( bulk_cloud_model  .AND.  microphysics_seifert  .AND.            &
1710               .NOT. downward)  THEN
1711             !$OMP PARALLEL DO PRIVATE( i, j )
1712             DO  m = 1, surf%ns   
1713
1714                i    = surf%i(m)           
1715                j    = surf%j(m)
1716
1717                surf%qrsws(m) = -surf%qrs(m) * surf%us(m)
1718                surf%nrsws(m) = -surf%nrs(m) * surf%us(m)
1719             ENDDO
1720          ENDIF
1721
1722!
1723!--       Bottom boundary condition for the TKE.
1724          IF ( ibc_e_b == 2 )  THEN
1725             !$OMP PARALLEL DO PRIVATE( i, j, k )
1726             DO  m = 1, surf%ns   
1727
1728                i    = surf%i(m)           
1729                j    = surf%j(m)
1730                k    = surf%k(m)
1731
1732                e(k,j,i) = ( surf%us(m) / 0.1_wp )**2
1733!
1734!--             As a test: cm = 0.4
1735!               e(k,j,i) = ( us(j,i) / 0.4_wp )**2
1736                e(k-1,j,i)   = e(k,j,i)
1737
1738             ENDDO
1739          ENDIF
1740!
1741!--    Calcuate surface fluxes at vertical surfaces. No stability is considered.
1742       ELSE
1743!
1744!--       Compute usvs l={0,1} and vsus l={2,3}
1745          IF ( mom_uv )  THEN
1746!
1747!--          Generalize computation by introducing flags. At north- and south-
1748!--          facing surfaces u-component is used, at east- and west-facing
1749!--          surfaces v-component is used.
1750             flag_u = MERGE( 1.0_wp, 0.0_wp, l == 0  .OR.  l == 1 )   
1751             flag_v = MERGE( 1.0_wp, 0.0_wp, l == 2  .OR.  l == 3 )   
1752             !$OMP PARALLEL  DO PRIVATE( i, j, k, z_mo )
1753             DO  m = 1, surf%ns 
1754                i = surf%i(m)           
1755                j = surf%j(m)
1756                k = surf%k(m)
1757
1758                z_mo = surf%z_mo(m)
1759
1760                surf%mom_flux_uv(m) = kappa *                                  &
1761                                ( flag_u * u(k,j,i) + flag_v * v(k,j,i) )  /   &
1762                                                        LOG( z_mo / surf%z0(m) )
1763
1764               surf%mom_flux_uv(m) =                                           &
1765                                    - surf%mom_flux_uv(m) * surf%us(m)
1766             ENDDO
1767          ENDIF
1768!
1769!--       Compute wsus l={0,1} and wsvs l={2,3}
1770          IF ( mom_w )  THEN
1771             !$OMP PARALLEL  DO PRIVATE( i, j, k, z_mo )
1772             DO  m = 1, surf%ns 
1773                i = surf%i(m)           
1774                j = surf%j(m)
1775                k = surf%k(m)
1776
1777                z_mo = surf%z_mo(m)
1778
1779                surf%mom_flux_w(m) = kappa * w(k,j,i) / LOG( z_mo / surf%z0(m) )
1780
1781                surf%mom_flux_w(m) =                                           &
1782                                     - surf%mom_flux_w(m) * surf%us(m)
1783             ENDDO
1784          ENDIF
1785!
1786!--       Compute momentum fluxes used for subgrid-scale TKE production at
1787!--       vertical surfaces. In constrast to the calculated momentum fluxes at
1788!--       vertical surfaces before, which are defined on the u/v/w-grid,
1789!--       respectively), the TKE fluxes are defined at the scalar grid.
1790!--       
1791          IF ( mom_tke )  THEN
1792!
1793!--          Precalculate velocity components at scalar grid point.
1794             ALLOCATE( u_i(1:surf%ns) )
1795             ALLOCATE( v_i(1:surf%ns) )
1796             ALLOCATE( w_i(1:surf%ns) )
1797
1798             IF ( l == 0  .OR.  l == 1 )  THEN
1799                !$OMP PARALLEL DO PRIVATE( i, j, k )
1800                DO  m = 1, surf%ns 
1801                   i = surf%i(m)           
1802                   j = surf%j(m)
1803                   k = surf%k(m)
1804
1805                   u_i(m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
1806                   v_i(m) = 0.0_wp
1807                   w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
1808                ENDDO
1809             ELSE
1810                !$OMP PARALLEL DO PRIVATE( i, j, k )
1811                DO  m = 1, surf%ns 
1812                   i = surf%i(m)           
1813                   j = surf%j(m)
1814                   k = surf%k(m)
1815
1816                   u_i(m) = 0.0_wp
1817                   v_i(m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
1818                   w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
1819                ENDDO
1820             ENDIF
1821
1822             !$OMP PARALLEL DO PRIVATE( i, j, dum, z_mo )
1823             DO  m = 1, surf%ns 
1824                i = surf%i(m)           
1825                j = surf%j(m)
1826
1827                z_mo = surf%z_mo(m)
1828
1829                dum = kappa / LOG( z_mo / surf%z0(m) )
1830!
1831!--             usvs (l=0,1) and vsus (l=2,3)
1832                surf%mom_flux_tke(0,m) = dum * ( u_i(m) + v_i(m) )
1833!
1834!--             wsvs (l=0,1) and wsus (l=2,3)
1835                surf%mom_flux_tke(1,m) = dum * w_i(m)
1836
1837                surf%mom_flux_tke(0:1,m) =                                     &
1838                               - surf%mom_flux_tke(0:1,m) * surf%us(m)
1839             ENDDO
1840!
1841!--          Deallocate temporary arrays
1842             DEALLOCATE( u_i )             
1843             DEALLOCATE( v_i ) 
1844             DEALLOCATE( w_i ) 
1845          ENDIF
1846       ENDIF
1847
1848    END SUBROUTINE calc_surface_fluxes
1849
1850   
1851!------------------------------------------------------------------------------!
1852! Description:
1853! ------------
1854!> Calculates temperature near surface (10 cm) for indoor model or 2 m
1855!> temperature for output
1856!------------------------------------------------------------------------------!
1857    SUBROUTINE calc_pt_near_surface ( z_char )
1858
1859       IMPLICIT NONE
1860
1861       CHARACTER (LEN = *), INTENT(IN) :: z_char !< string identifier to identify z level
1862       INTEGER(iwp)                    :: i      !< grid index x-dimension
1863       INTEGER(iwp)                    :: j      !< grid index y-dimension
1864       INTEGER(iwp)                    :: k      !< grid index z-dimension
1865       INTEGER(iwp)                    :: kk     !< running index along the z-dimension
1866       INTEGER(iwp)                    :: m      !< running index for surface elements
1867
1868       
1869       SELECT CASE ( z_char)
1870           
1871       
1872          CASE ( '10cm' )
1873
1874             DO  m = 1, surf%ns
1875
1876                i = surf%i(m)
1877                j = surf%j(m)
1878                k = surf%k(m)
1879
1880                surf%pt_10cm(m) = surf%pt_surface(m) + surf%ts(m) / kappa      &
1881                                   * ( LOG( 0.1_wp /  surf%z0h(m) )            &
1882                                     - psi_h( 0.1_wp / surf%ol(m) )            &
1883                                     + psi_h( surf%z0h(m) / surf%ol(m) ) )
1884                                   
1885             ENDDO
1886
1887!
1888!--       2-m temperature. Note, this is only calculated for output reasons at
1889!--       horizontal upward-facing surfaces.
1890          CASE ( '2m' )
1891     
1892             DO  m = 1, surf%ns
1893
1894                i = surf%i(m)
1895                j = surf%j(m)
1896                k = surf%k(m)
1897!
1898!--             If 2-m level is below the first grid level, MOST is
1899!--             used for calculation of 2-m temperature.
1900                IF ( surf%z_mo(m) > 2.0_wp )  THEN
1901                   surf%pt_2m(m) = surf%pt_surface(m) + surf%ts(m) / kappa     &
1902                                   * ( LOG( 2.0_wp /  surf%z0h(m) )            &
1903                                     - psi_h( 2.0_wp / surf%ol(m) )            &
1904                                     + psi_h( surf%z0h(m) / surf%ol(m) ) )
1905!
1906!--             If 2-m level is above the first grid level, 2-m temperature
1907!--             is linearly interpolated between the two nearest vertical grid
1908!--             levels. Note, since 2-m temperature is only computed for
1909!--             horizontal upward-facing surfaces, only a vertical
1910!--             interpolation is necessary.
1911                ELSE
1912!
1913!--                zw(k-1) defines the height of the surface.
1914                   kk = k
1915                   DO WHILE ( zu(kk) - zw(k-1) < 2.0_wp  .AND.  kk <= nzt )
1916                      kk = kk + 1 
1917                   ENDDO               
1918!
1919!--                kk defines the index of the first grid level >= 2m.
1920                   surf%pt_2m(m) = pt(kk-1,j,i) +                              &
1921                                 ( zw(k-1) + 2.0_wp - zu(kk-1)     ) *         &
1922                                 ( pt(kk,j,i)       - pt(kk-1,j,i) ) /         &
1923                                 ( zu(kk)           - zu(kk-1)     )
1924                ENDIF
1925
1926             ENDDO
1927         
1928       
1929       END SELECT
1930
1931    END SUBROUTINE calc_pt_near_surface
1932   
1933
1934!
1935!-- Integrated stability function for momentum
1936    FUNCTION psi_m( zeta ) 
1937       !$ACC ROUTINE SEQ
1938       
1939       USE kinds
1940
1941       IMPLICIT NONE
1942
1943       REAL(wp)            :: psi_m !< Integrated similarity function result
1944       REAL(wp)            :: zeta  !< Stability parameter z/L
1945       REAL(wp)            :: x     !< dummy variable
1946
1947       REAL(wp), PARAMETER :: a = 1.0_wp            !< constant
1948       REAL(wp), PARAMETER :: b = 0.66666666666_wp  !< constant
1949       REAL(wp), PARAMETER :: c = 5.0_wp            !< constant
1950       REAL(wp), PARAMETER :: d = 0.35_wp           !< constant
1951       REAL(wp), PARAMETER :: c_d_d = c / d         !< constant
1952       REAL(wp), PARAMETER :: bc_d_d = b * c / d    !< constant
1953
1954
1955       IF ( zeta < 0.0_wp )  THEN
1956          x = SQRT( SQRT( 1.0_wp  - 16.0_wp * zeta ) )
1957          psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2    &
1958                  * ( 1.0_wp + x**2 ) * 0.125_wp )
1959       ELSE
1960
1961          psi_m = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - a * zeta         &
1962                   - bc_d_d
1963!
1964!--       Old version for stable conditions (only valid for z/L < 0.5)
1965!--       psi_m = - 5.0_wp * zeta
1966
1967       ENDIF
1968
1969    END FUNCTION psi_m
1970
1971
1972!
1973!-- Integrated stability function for heat and moisture
1974    FUNCTION psi_h( zeta ) 
1975       !$ACC ROUTINE SEQ
1976       
1977       USE kinds
1978
1979       IMPLICIT NONE
1980
1981       REAL(wp)            :: psi_h !< Integrated similarity function result
1982       REAL(wp)            :: zeta  !< Stability parameter z/L
1983       REAL(wp)            :: x     !< dummy variable
1984
1985       REAL(wp), PARAMETER :: a = 1.0_wp            !< constant
1986       REAL(wp), PARAMETER :: b = 0.66666666666_wp  !< constant
1987       REAL(wp), PARAMETER :: c = 5.0_wp            !< constant
1988       REAL(wp), PARAMETER :: d = 0.35_wp           !< constant
1989       REAL(wp), PARAMETER :: c_d_d = c / d         !< constant
1990       REAL(wp), PARAMETER :: bc_d_d = b * c / d    !< constant
1991
1992
1993       IF ( zeta < 0.0_wp )  THEN
1994          x = SQRT( 1.0_wp  - 16.0_wp * zeta )
1995          psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
1996       ELSE
1997          psi_h = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - (1.0_wp          &
1998                  + 0.66666666666_wp * a * zeta )**1.5_wp - bc_d_d             &
1999                  + 1.0_wp
2000!
2001!--       Old version for stable conditions (only valid for z/L < 0.5)
2002!--       psi_h = - 5.0_wp * zeta
2003       ENDIF
2004
2005    END FUNCTION psi_h
2006
2007
2008!------------------------------------------------------------------------------!
2009! Description:
2010! ------------
2011!> Calculates stability function for momentum
2012!>
2013!> @author Hauke Wurps
2014!------------------------------------------------------------------------------!
2015    FUNCTION phi_m( zeta ) 
2016       !$ACC ROUTINE SEQ
2017   
2018       IMPLICIT NONE
2019   
2020       REAL(wp)            :: phi_m         !< Value of the function
2021       REAL(wp)            :: zeta          !< Stability parameter z/L
2022   
2023       REAL(wp), PARAMETER :: a = 16.0_wp   !< constant
2024       REAL(wp), PARAMETER :: c = 5.0_wp    !< constant
2025   
2026       IF ( zeta < 0.0_wp )  THEN
2027          phi_m = 1.0_wp / SQRT( SQRT( 1.0_wp - a * zeta ) )
2028       ELSE
2029          phi_m = 1.0_wp + c * zeta
2030       ENDIF
2031   
2032    END FUNCTION phi_m
2033
2034 END MODULE surface_layer_fluxes_mod
Note: See TracBrowser for help on using the repository browser.