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

Last change on this file since 4281 was 4258, checked in by suehring, 20 months ago

Land-surface model: Revise limitation for soil moisture in case it exceeds its saturation value; Revise initialization of soil moisture and temperature in a nested run in case dynamic input information is available. This case, the soil within the child domains can be initialized separately; As part of this revision, migrate the netcdf input of soil temperature / moisture to this module, as well as the routine to inter/extrapolate soil profiles between different grids.; Plant-canopy: Check if any LAD is prescribed when plant-canopy model is applied, in order to avoid crashes in the radiation transfer model; Surface-layer fluxes: Initialization of Obukhov length also at vertical surfaces (if allocated); Urban-surface model: Add checks to ensure that relative fractions of walls, windowns and green surfaces sum-u to one; Revise message calls dealing with local checks

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