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

Last change on this file since 4692 was 4691, checked in by suehring, 4 years ago

Reference model topography to the lowest grid point also in ASCII input case

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