source: palm/trunk/SOURCE/surface_layer_fluxes_mod.f90

Last change on this file was 4828, checked in by Giersch, 3 years ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

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