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

Last change on this file since 4411 was 4370, checked in by raasch, 4 years ago

bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added

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