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

Last change on this file since 3780 was 3745, checked in by suehring, 6 years ago

document last commit

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