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

Last change on this file since 4360 was 4360, checked in by suehring, 16 months ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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