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

Last change on this file since 3787 was 3787, checked in by raasch, 2 years ago

unused variables removed

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