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

Last change on this file since 4631 was 4594, checked in by suehring, 4 years ago

Bugfix, add acc directives for scalar-roughness length; Include k index in OMP PRIVATE statements

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