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

Last change on this file since 4366 was 4366, checked in by raasch, 3 years ago

code vectorization for NEC Aurora: vectorized version of Temperton FFT, vectorization of Newtor iteration for calculating the Obukhov length

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