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

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