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

Last change on this file since 4583 was 4562, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

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