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

Last change on this file since 1930 was 1930, checked in by suehring, 5 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 43.2 KB
Line 
1!> @file surface_layer_fluxes_mod.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!
18!--------------------------------------------------------------------------------!
19! Current revisions:
20! ------------------
21
22!
23! Former revisions:
24! -----------------
25! $Id: surface_layer_fluxes_mod.f90 1930 2016-06-09 16:32:12Z suehring $
26!
27! 1929 2016-06-09 16:25:25Z suehring
28! Bugfix: avoid segmentation fault in case one grid point is horizontally
29! completely surrounded by topography
30!
31! 1920 2016-05-30 10:50:15Z suehring
32! Avoid segmentation fault (see change in 1915) by different initialization of
33! us instead of adding a very small number in the denominator
34!
35! 1915 2016-05-27 11:05:02Z suehring
36! Bugfix: avoid segmentation fault in case of most_method = 'circular' at first
37! timestep
38!
39! 1850 2016-04-08 13:29:27Z maronga
40! Module renamed
41!
42!
43! 1822 2016-04-07 07:49:42Z hoffmann
44! icloud_scheme replaced by microphysics_*
45!
46! 1788 2016-03-10 11:01:04Z maronga
47! Added parameter z0q which replaces z0h in the similarity functions for
48! humidity.
49! Syntax layout improved.
50!
51! 1757 2016-02-22 15:49:32Z maronga
52! Minor fixes.
53!
54! 1749 2016-02-09 12:19:56Z raasch
55! further OpenACC adjustments
56!
57! 1747 2016-02-08 12:25:53Z raasch
58! adjustments for OpenACC usage
59!
60! 1709 2015-11-04 14:47:01Z maronga
61! Bugfix: division by zero could occur when calculating rib at low wind speeds
62! Bugfix: calculation of uv_total for neutral = .T., initial value for ol for
63! neutral = .T.
64!
65! 1705 2015-11-02 14:28:56Z maronga
66! Typo removed
67!
68! 1697 2015-10-28 17:14:10Z raasch
69! FORTRAN and OpenMP errors removed
70!
71! 1696 2015-10-27 10:03:34Z maronga
72! Modularized and completely re-written version of prandtl_fluxes.f90. In the
73! course of the re-writing two additional methods have been implemented. See
74! updated description.
75!
76! 1551 2015-03-03 14:18:16Z maronga
77! Removed land surface model part. The surface fluxes are now always calculated
78! within prandtl_fluxes, based on the given surface temperature/humidity (which
79! is either provided by the land surface model, by large scale forcing data, or
80! directly prescribed by the user.
81!
82! 1496 2014-12-02 17:25:50Z maronga
83! Adapted for land surface model
84!
85! 1494 2014-11-21 17:14:03Z maronga
86! Bugfixes: qs is now calculated before calculation of Rif. calculation of
87! buoyancy flux in Rif corrected (added missing humidity term), allow use of
88! topography for coupled runs (not tested)
89!
90! 1361 2014-04-16 15:17:48Z hoffmann
91! Bugfix: calculation of turbulent fluxes of rain water content (qrsws) and rain
92! drop concentration (nrsws) added
93!
94! 1340 2014-03-25 19:45:13Z kanani
95! REAL constants defined as wp-kind
96!
97! 1320 2014-03-20 08:40:49Z raasch
98! ONLY-attribute added to USE-statements,
99! kind-parameters added to all INTEGER and REAL declaration statements,
100! kinds are defined in new module kinds,
101! old module precision_kind is removed,
102! revision history before 2012 removed,
103! comment fields (!:) to be used for variable explanations added to
104! all variable declaration statements
105!
106! 1276 2014-01-15 13:40:41Z heinze
107! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
108!
109! 1257 2013-11-08 15:18:40Z raasch
110! openACC "kernels do" replaced by "kernels loop", "loop independent" added
111!
112! 1036 2012-10-22 13:43:42Z raasch
113! code put under GPL (PALM 3.9)
114!
115! 1015 2012-09-27 09:23:24Z raasch
116! OpenACC statements added
117!
118! 978 2012-08-09 08:28:32Z fricke
119! roughness length for scalar quantities z0h added
120!
121! Revision 1.1  1998/01/23 10:06:06  raasch
122! Initial revision
123!
124!
125! Description:
126! ------------
127!> Diagnostic computation of vertical fluxes in the constant flux layer from the
128!> values of the variables at grid point k=1. Three different methods are
129!> available:
130!> 1) the "old" version (most_method = 'circular') which is fast, but inaccurate
131!> 2) a Newton iteration method (most_method = 'newton'), which is accurate, but
132!>    slower
133!> 3) a method using a lookup table which is fast and accurate. Note, however,
134!>    that this method cannot be used in case of roughness heterogeneity
135!>
136!> @todo (re)move large_scale_forcing actions
137!> @todo check/optimize OpenMP and OpenACC directives
138!------------------------------------------------------------------------------!
139 MODULE surface_layer_fluxes_mod
140
141    USE arrays_3d,                                                             &
142        ONLY:  e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, &
143               shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h, z0q
144
145    USE cloud_parameters,                                                      &
146        ONLY:  l_d_cp, pt_d_t
147
148    USE constants,                                                             &
149        ONLY:  pi
150
151    USE cpulog
152
153    USE control_parameters,                                                    &
154        ONLY:  cloud_physics, constant_heatflux, constant_waterflux,           &
155               coupling_mode, g, humidity, ibc_e_b, ibc_pt_b,                  &
156               initializing_actions, kappa, intermediate_timestep_count,       &
157               intermediate_timestep_count_max, large_scale_forcing, lsf_surf, &
158               message_string, microphysics_seifert, most_method, neutral,     &
159               passive_scalar, pt_surface, q_surface, run_coupled,             &
160               surface_pressure, simulated_time, terminate_run, zeta_max,      &
161               zeta_min
162
163    USE indices,                                                               &
164        ONLY:  nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_s_inner,        &
165               nzb_u_inner, nzb_v_inner
166
167    USE kinds
168
169    USE pegrid
170
171    USE land_surface_model_mod,                                                &
172        ONLY:  land_surface, skip_time_do_lsm
173
174    IMPLICIT NONE
175
176    INTEGER(iwp) ::  i            !< loop index x direction
177    INTEGER(iwp) ::  j            !< loop index y direction
178    INTEGER(iwp) ::  k            !< loop index z direction
179
180    INTEGER(iwp), PARAMETER     :: num_steps = 15000  !< number of steps in the lookup table
181
182    LOGICAL      ::  coupled_run  !< Flag for coupled atmosphere-ocean runs
183
184    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt1,      & !< Potential temperature at first grid level (required for cloud_physics = .T.)
185                                             qv1,      & !< Specific humidity at first grid level (required for cloud_physics = .T.)
186                                             uv_total    !< Total velocity at first grid level
187
188    REAL(wp), DIMENSION(0:num_steps-1) :: rib_tab,  & !< Lookup table bulk Richardson number
189                                          ol_tab      !< Lookup table values of L
190
191    REAL(wp)     ::  e_s,               & !< Saturation water vapor pressure
192                     l_bnd  = 7500,     & !< Lookup table index of the last time step
193                     ol_max = 1.0E6_wp, & !< Maximum Obukhov length
194                     rib_max,           & !< Maximum Richardson number in lookup table
195                     rib_min,           & !< Minimum Richardson number in lookup table
196                     z_mo                 !< Height of the constant flux layer where MOST is assumed
197
198
199    SAVE
200
201    PRIVATE
202
203    PUBLIC init_surface_layer_fluxes, pt1, qv1, surface_layer_fluxes, uv_total
204
205    INTERFACE init_surface_layer_fluxes
206       MODULE PROCEDURE init_surface_layer_fluxes
207    END INTERFACE init_surface_layer_fluxes
208
209    INTERFACE surface_layer_fluxes
210       MODULE PROCEDURE surface_layer_fluxes
211    END INTERFACE surface_layer_fluxes
212
213
214 CONTAINS
215
216
217!------------------------------------------------------------------------------!
218! Description:
219! ------------
220!> Main routine to compute the surface fluxes
221!------------------------------------------------------------------------------!
222    SUBROUTINE surface_layer_fluxes
223
224       IMPLICIT NONE
225
226!
227!--    In case cloud physics is used, it is required to derive potential
228!--    temperature and specific humidity at first grid level from the fields pt
229!--    and q
230       IF ( cloud_physics )  THEN
231          CALL calc_pt_q
232       ENDIF
233
234!
235!--    First, calculate the new Obukhov length, then new friction velocity,
236!--    followed by the new scaling parameters (th*, q*, etc.), and the new
237!--    surface fluxes if required. The old routine ("circular") requires a
238!--    different order of calls as the scaling parameters from the previous time
239!--    steps are used to calculate the Obukhov length
240
241!
242!--    Depending on setting of most_method use the "old" routine
243       IF ( most_method == 'circular' )  THEN
244
245          CALL calc_scaling_parameters
246
247          CALL calc_uv_total
248
249          IF ( .NOT. neutral )  THEN
250             CALL calc_ol
251          ENDIF
252
253          CALL calc_us
254
255          CALL calc_surface_fluxes
256
257!
258!--    Use either Newton iteration or a lookup table for the bulk Richardson
259!--    number to calculate the Obukhov length
260       ELSEIF ( most_method == 'newton'  .OR.  most_method == 'lookup' )  THEN
261
262          CALL calc_uv_total
263
264          IF ( .NOT. neutral )  THEN
265             CALL calc_ol
266          ENDIF
267
268          CALL calc_us
269
270          CALL calc_scaling_parameters
271
272          CALL calc_surface_fluxes
273
274       ENDIF
275
276    END SUBROUTINE surface_layer_fluxes
277
278
279!------------------------------------------------------------------------------!
280! Description:
281! ------------
282!> Initializing actions for the surface layer routine. Basically, this involves
283!> the preparation of a lookup table for the the bulk Richardson number vs
284!> Obukhov length L when using the lookup table method.
285!------------------------------------------------------------------------------!
286    SUBROUTINE init_surface_layer_fluxes
287
288       IMPLICIT NONE
289
290       INTEGER(iwp) :: l,          & !< Index for loop to create lookup table
291                       num_steps_n   !< Number of non-stretched zeta steps
292
293       LOGICAL :: terminate_run_l = .FALSE.    !< Flag to terminate run (global)
294
295       REAL(wp), PARAMETER ::  zeta_stretch = -10.0_wp !< Start of stretching in the free convection limit
296                               
297       REAL(wp), DIMENSION(:), ALLOCATABLE :: zeta_tmp
298
299
300       REAL(wp) :: zeta_step,            & !< Increment of zeta
301                   regr      = 1.01_wp,  & !< Stretching factor of zeta_step in the free convection limit
302                   regr_old  = 1.0E9_wp, & !< Stretching factor of last iteration step
303                   z0h_min   = 0.0_wp,   & !< Minimum value of z0h to create table
304                   z0_min    = 0.0_wp      !< Minimum value of z0 to create table
305!
306!--    When cloud physics is used, arrays for storing potential temperature and
307!--    specific humidity at first grid level are required
308       IF ( cloud_physics )  THEN
309          ALLOCATE ( pt1(nysg:nyng,nxlg:nxrg) )
310          ALLOCATE ( qv1(nysg:nyng,nxlg:nxrg) )
311       ENDIF
312
313!
314!--    Allocate field for storing the horizontal velocity
315       ALLOCATE ( uv_total(nysg:nyng,nxlg:nxrg) )
316
317
318!
319!--    In case of runs with neutral statification, set Obukhov length to a
320!--    large value
321       IF ( neutral ) ol = 1.0E10_wp
322
323       IF ( most_method == 'lookup' )  THEN
324
325!
326!--       Check for roughness heterogeneity. In that case terminate run and
327!--       inform user
328          IF ( MINVAL( z0h ) /= MAXVAL( z0h )  .OR.                            &
329               MINVAL( z0  ) /= MAXVAL( z0  ) )  THEN
330             terminate_run_l = .TRUE.
331          ENDIF
332
333#if defined( __parallel )
334!
335!--       Make a logical OR for all processes. Force termiation of model if result
336!--       is TRUE
337          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
338          CALL MPI_ALLREDUCE( terminate_run_l, terminate_run, 1, MPI_LOGICAL,  &
339                              MPI_LOR, comm2d, ierr )
340#else
341          terminate_run = terminate_run_l
342#endif
343
344          IF ( terminate_run )  THEN
345             message_string = 'most_method = "lookup" cannot be used in ' //   &
346                              'combination with a prescribed roughness '  //   &
347                              'heterogeneity'
348             CALL message( 'surface_layer_fluxes', 'PA0417', 1, 2, 0, 6, 0 )
349          ENDIF
350
351          ALLOCATE(  zeta_tmp(0:num_steps-1) )
352
353!
354!--       Use the lowest possible value for z_mo
355          k    = MINVAL(nzb_s_inner)
356          z_mo = zu(k+1) - zw(k)
357
358!
359!--       Calculate z/L range from zeta_stretch to zeta_max using 90% of the
360!--       available steps (num_steps). The calculation is done with negative
361!--       values of zeta in order to simplify the stretching in the free
362!--       convection limit for the remaining 10% of steps.
363          zeta_tmp(0) = - zeta_max
364          num_steps_n = ( num_steps * 9 / 10 ) - 1
365          zeta_step   = (zeta_max - zeta_stretch) / REAL(num_steps_n)
366
367          DO l = 1, num_steps_n
368             zeta_tmp(l) = zeta_tmp(l-1) + zeta_step
369          ENDDO
370
371!
372!--       Calculate stretching factor for the free convection range
373          DO  WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
374             regr_old = regr
375             regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr )   &
376                    )**( 10.0_wp / REAL(num_steps) )
377          ENDDO
378
379!
380!--       Calculate z/L range from zeta_min to zeta_stretch
381          DO l = num_steps_n+1, num_steps-1
382             zeta_tmp(l) = zeta_tmp(l-1) + zeta_step
383             zeta_step = zeta_step * regr
384          ENDDO
385
386!
387!--       Save roughness lengths to temporary variables
388          z0h_min = z0h(nys,nxl)
389          z0_min  = z0(nys,nxl)
390         
391!
392!--       Calculate lookup table for the Richardson number versus Obukhov length
393!--       The Richardson number (rib) is defined depending on the choice of
394!--       boundary conditions for temperature
395          IF ( ibc_pt_b == 1 )  THEN
396             DO l = 0, num_steps-1
397                ol_tab(l)  = - z_mo / zeta_tmp(num_steps-1-l)
398                rib_tab(l) = z_mo / ol_tab(l)  / ( LOG( z_mo / z0_min )        &
399                                                - psi_m( z_mo / ol_tab(l) )    &
400                                                + psi_m( z0_min / ol_tab(l) )  &
401                                                  )**3
402             ENDDO 
403          ELSE
404             DO l = 0, num_steps-1
405                ol_tab(l)  = - z_mo / zeta_tmp(num_steps-1-l)
406                rib_tab(l) = z_mo / ol_tab(l)  * ( LOG( z_mo / z0h_min )       &
407                                              - psi_h( z_mo / ol_tab(l) )      &
408                                              + psi_h( z0h_min / ol_tab(l) )   &
409                                            )                                  &
410                                          / ( LOG( z_mo / z0_min )             &
411                                              - psi_m( z_mo / ol_tab(l) )      &
412                                              + psi_m( z0_min / ol_tab(l) )    &
413                                            )**2
414             ENDDO
415          ENDIF
416
417!
418!--       Determine minimum values of rib in the lookup table. Set upper limit
419!--       to critical Richardson number (0.25)
420          rib_min  = MINVAL(rib_tab)
421          rib_max  = 0.25 !MAXVAL(rib_tab)
422
423          DEALLOCATE( zeta_tmp )
424       ENDIF
425
426    END SUBROUTINE init_surface_layer_fluxes
427
428
429!------------------------------------------------------------------------------!
430! Description:
431! ------------
432!> Compute the absolute value of the horizontal velocity (relative to the
433!> surface). This is required by all methods
434!------------------------------------------------------------------------------!
435    SUBROUTINE calc_uv_total
436
437       IMPLICIT NONE
438
439
440       !$OMP PARALLEL DO PRIVATE( k )
441       !$acc kernels loop present( nzb_s_inner, u, uv_total, v ) private( j, k )
442       DO  i = nxl, nxr
443          DO  j = nys, nyn
444
445             k   = nzb_s_inner(j,i)
446             uv_total(j,i) = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)      &
447                                         - u(k,j,i)   - u(k,j,i+1) ) )**2 +    &
448                              ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)           &
449                                         - v(k,j,i)   - v(k,j+1,i) ) )**2 )
450
451!
452!--          For too small values of the local wind, MOST does not work. A
453!--          threshold value is thus set if required
454!            uv_total(j,i) = MAX(0.01_wp,uv_total(j,i))
455
456          ENDDO
457       ENDDO
458
459!
460!--    Values of uv_total need to be exchanged at the ghost boundaries
461       !$acc update host( uv_total )
462       CALL exchange_horiz_2d( uv_total )
463       !$acc update device( uv_total )
464
465    END SUBROUTINE calc_uv_total
466
467
468!------------------------------------------------------------------------------!
469! Description:
470! ------------
471!> Calculate the Obukhov length (L) and Richardson flux number (z/L)
472!------------------------------------------------------------------------------!
473    SUBROUTINE calc_ol
474
475       IMPLICIT NONE
476
477       INTEGER(iwp) :: iter,  & !< Newton iteration step
478                       l        !< look index
479
480       REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) :: rib !< Bulk Richardson number
481
482       REAL(wp)     :: f,      & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
483                       f_d_ol, & !< Derivative of f
484                       ol_l,   & !< Lower bound of L for Newton iteration
485                       ol_m,   & !< Previous value of L for Newton iteration
486                       ol_old, & !< Previous time step value of L
487                       ol_u      !< Upper bound of L for Newton iteration
488
489       IF ( TRIM( most_method ) /= 'circular' )  THEN
490     
491          !$acc data present( nzb_s_inner, pt, q, qsws, rib, shf, uv_total, vpt, zu, zw )
492
493          !$OMP PARALLEL DO PRIVATE( k, z_mo )
494          !$acc kernels loop private( j, k, z_mo )
495          DO  i = nxl, nxr
496             DO  j = nys, nyn
497
498                k   = nzb_s_inner(j,i)
499                z_mo = zu(k+1) - zw(k)
500
501!
502!--             Evaluate bulk Richardson number (calculation depends on
503!--             definition based on setting of boundary conditions
504                IF ( ibc_pt_b /= 1 )  THEN
505                   IF ( humidity )  THEN
506                      rib(j,i) = g * z_mo * ( vpt(k+1,j,i) - vpt(k,j,i) )      &
507                           / ( uv_total(j,i)**2 * vpt(k+1,j,i) + 1.0E-20_wp )
508                   ELSE
509                      rib(j,i) = g * z_mo * ( pt(k+1,j,i) - pt(k,j,i) )        &
510                           / ( uv_total(j,i)**2 * pt(k+1,j,i)  + 1.0E-20_wp )
511                   ENDIF     
512                ELSE
513!
514!--                When using Neumann boundary conditions, the buoyancy flux
515!--                is required but cannot be calculated at the surface, as pt
516!--                and q are not known at the surface. Hence the values at
517!--                first grid level are used to estimate the buoyancy flux
518                   IF ( humidity )  THEN
519                      rib(j,i) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
520                                 * q(k+1,j,i) ) * shf(j,i) + 0.61_wp           &
521                                 * pt(k+1,j,i) * qsws(j,i) )   &
522                                 / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2&
523                                 + 1.0E-20_wp)
524                   ELSE
525                      rib(j,i) = - g * z_mo * shf(j,i)                         &
526                           / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2       &
527                           + 1.0E-20_wp )
528                   ENDIF
529                ENDIF 
530     
531             ENDDO
532          ENDDO
533          !$acc end data
534
535       ENDIF
536
537!
538!--    Calculate the Obukhov length either using a Newton iteration
539!--    method, via a lookup table, or using the old circular way
540       IF ( TRIM( most_method ) == 'newton' )  THEN
541
542          !$OMP PARALLEL DO PRIVATE( k, z_mo )
543          !# WARNING: does not work on GPU so far because of DO-loop with
544          !#          undetermined iterations
545          !!!!!!$acc kernels loop
546          DO  i = nxl, nxr
547             DO  j = nys, nyn
548
549                k   = nzb_s_inner(j,i)
550                z_mo = zu(k+1) - zw(k)
551
552!
553!--             Store current value in case the Newton iteration fails
554                ol_old = ol(j,i)
555
556!
557!--             Ensure that the bulk Richardson number and the Obukhov
558!--             lengtH have the same sign
559                IF ( rib(j,i) * ol(j,i) < 0.0_wp  .OR.                         &
560                     ABS( ol(j,i) ) == ol_max )  THEN
561                   IF ( rib(j,i) > 0.0_wp ) ol(j,i) =  0.01_wp
562                   IF ( rib(j,i) < 0.0_wp ) ol(j,i) = -0.01_wp
563                ENDIF
564!
565!--             Iteration to find Obukhov length
566                iter = 0
567                DO
568                   iter = iter + 1
569!
570!--                In case of divergence, use the value of the previous time step
571                   IF ( iter > 1000 )  THEN
572                      ol(j,i) = ol_old
573                      EXIT
574                   ENDIF
575
576                   ol_m = ol(j,i)
577                   ol_l = ol_m - 0.001_wp * ol_m
578                   ol_u = ol_m + 0.001_wp * ol_m
579
580
581                   IF ( ibc_pt_b /= 1 )  THEN
582!
583!--                   Calculate f = Ri - [...]/[...]^2 = 0
584                      f = rib(j,i) - ( z_mo / ol_m ) * ( LOG( z_mo / z0h(j,i) )&
585                                                    - psi_h( z_mo / ol_m )     &
586                                                    + psi_h( z0h(j,i) / ol_m ) &
587                                                   )                           &
588                                                 / ( LOG( z_mo / z0(j,i) )     &
589                                                    - psi_m( z_mo / ol_m )     &
590                                                    + psi_m( z0(j,i) / ol_m )  &
591                                                    )**2
592
593!
594!--                    Calculate df/dL
595                       f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / z0h(j,i) ) &
596                                                   - psi_h( z_mo / ol_u )      &
597                                                   + psi_h( z0h(j,i) / ol_u )  &
598                                                 )                             &
599                                               / ( LOG( z_mo / z0(j,i) )       &
600                                                   - psi_m( z_mo / ol_u )      &
601                                                   + psi_m( z0(j,i) / ol_u )   &
602                                                 )**2                          &
603                              + ( z_mo / ol_l ) * ( LOG( z_mo / z0h(j,i) )     &
604                                                   - psi_h( z_mo / ol_l )      &
605                                                   + psi_h( z0h(j,i) / ol_l )  &
606                                                 )                             &
607                                               / ( LOG( z_mo / z0(j,i) )       &
608                                                   - psi_m( z_mo / ol_l )      &
609                                                   + psi_m( z0(j,i) / ol_l )   &
610                                                 )**2                          &
611                                ) / ( ol_u - ol_l )
612                   ELSE
613!
614!--                   Calculate f = Ri - 1 /[...]^3 = 0
615                      f = rib(j,i) - ( z_mo / ol_m ) / ( LOG( z_mo / z0(j,i) )&
616                                                    - psi_m( z_mo / ol_m )    &
617                                                    + psi_m( z0(j,i) / ol_m ) &
618                                                       )**3
619
620!
621!--                   Calculate df/dL
622                      f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / z0(j,i) )  &
623                                                   - psi_m( z_mo / ol_u )     &
624                                                   + psi_m( z0(j,i) / ol_u )  &
625                                                 )**3                         &
626                              + ( z_mo / ol_l ) / ( LOG( z_mo / z0(j,i) )     &
627                                                   - psi_m( z_mo / ol_l )     &
628                                                   + psi_m( z0(j,i) / ol_l )  &
629                                                 )**3                         &
630                                     ) / ( ol_u - ol_l )
631                   ENDIF
632!
633!--                Calculate new L
634                   ol(j,i) = ol_m - f / f_d_ol
635
636!
637!--                Ensure that the bulk Richardson number and the Obukhov
638!--                length have the same sign and ensure convergence.
639                   IF ( ol(j,i) * ol_m < 0.0_wp )  ol(j,i) = ol_m * 0.5_wp
640
641!
642!--                If unrealistic value occurs, set L to the maximum
643!--                value that is allowed
644                   IF ( ABS( ol(j,i) ) > ol_max )  THEN
645                      ol(j,i) = ol_max
646                      EXIT
647                   ENDIF
648!
649!--                Check for convergence
650                   IF ( ABS( ( ol(j,i) - ol_m ) / ol(j,i) ) < 1.0E-4_wp )  THEN
651                      EXIT
652                   ELSE
653                      CYCLE
654                   ENDIF
655
656                ENDDO
657                       
658             ENDDO
659          ENDDO
660
661       ELSEIF ( TRIM( most_method ) == 'lookup' )  THEN
662
663          !$OMP PARALLEL DO PRIVATE( k, z_mo )
664          !# WARNING: does not work on GPU so far because of DO  WHILE construct
665          !!!!!!$acc kernels loop
666          DO  i = nxl, nxr
667             DO  j = nys, nyn
668
669!
670!--             If the bulk Richardson number is outside the range of the lookup
671!--             table, set it to the exceeding threshold value
672                IF ( rib(j,i) < rib_min )  rib(j,i) = rib_min
673                IF ( rib(j,i) > rib_max )  rib(j,i) = rib_max
674
675!
676!--             Find the correct index bounds for linear interpolation. As the
677!--             Richardson number will not differ very much from time step to
678!--             time step , use the index from the last step and search in the
679!--             correct direction
680                l = l_bnd
681                IF ( rib_tab(l) - rib(j,i) > 0.0_wp )  THEN
682                   DO  WHILE ( rib_tab(l-1) - rib(j,i) > 0.0_wp  .AND.  l > 0 )
683                      l = l-1
684                   ENDDO
685                ELSE
686                   DO  WHILE ( rib_tab(l) - rib(j,i) < 0.0_wp                  &
687                              .AND.  l < num_steps-1 )
688                      l = l+1
689                   ENDDO
690                ENDIF
691                l_bnd = l
692
693!
694!--             Linear interpolation to find the correct value of z/L
695                ol(j,i) = ( ol_tab(l-1) + ( ol_tab(l) - ol_tab(l-1) )          &
696                            / (  rib_tab(l) - rib_tab(l-1) )                   &
697                            * ( rib(j,i) - rib_tab(l-1) ) )
698
699             ENDDO
700          ENDDO
701
702       ELSEIF ( TRIM( most_method ) == 'circular' )  THEN
703
704          !$OMP PARALLEL DO PRIVATE( k, z_mo )
705          !$acc kernels loop present( nzb_s_inner, ol, pt, pt1, q, ql, qs, qv1, ts, us, vpt, zu, zw ) private( j, k, z_mo )
706          DO  i = nxl, nxr
707             DO  j = nys, nyn
708
709                k   = nzb_s_inner(j,i)
710                z_mo = zu(k+1) - zw(k)
711
712                IF ( .NOT. humidity )  THEN
713                   ol(j,i) =  ( pt(k+1,j,i) *  us(j,i)**2 ) / ( kappa * g      &
714                              * ts(j,i) + 1E-30_wp )
715                ELSEIF ( cloud_physics )  THEN
716
717                   ol(j,i) =  ( vpt(k+1,j,i) * us(j,i)**2 ) / ( kappa * g      &
718                              * ( ts(j,i) + 0.61_wp * pt1(j,i) * qs(j,i)       &
719                              + 0.61_wp * qv1(j,i) * ts(j,i) - ts(j,i)         &
720                              * ql(k+1,j,i) ) + 1E-30_wp )
721                ELSE
722                   ol(j,i) =  ( vpt(k+1,j,i) *  us(j,i)**2 ) / ( kappa * g     &
723                              * ( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i)    &
724                                  + 0.61_wp * q(k+1,j,i) * ts(j,i) ) + 1E-30_wp )
725                ENDIF
726!
727!--             Limit the value range of the Obukhov length.
728!--             This is necessary for very small velocities (u,v --> 0), because
729!--             the absolute value of ol can then become very small, which in
730!--             consequence would result in very large shear stresses and very
731!--             small momentum fluxes (both are generally unrealistic).
732                IF ( ( z_mo / ( ol(j,i) + 1E-30_wp ) ) < zeta_min )            &
733                   ol(j,i) = z_mo / zeta_min
734                IF ( ( z_mo / ( ol(j,i) + 1E-30_wp ) ) > zeta_max )            &
735                   ol(j,i) = z_mo / zeta_max
736             ENDDO
737          ENDDO
738
739       ENDIF
740
741!
742!--    Values of ol at ghost point locations are needed for the evaluation
743!--    of usws and vsws.
744       !$acc update host( ol )
745       CALL exchange_horiz_2d( ol )
746       !$acc update device( ol )
747
748    END SUBROUTINE calc_ol
749
750!
751!-- Calculate friction velocity u*
752    SUBROUTINE calc_us
753
754       IMPLICIT NONE
755
756       !$OMP PARALLEL DO PRIVATE( k, z_mo )
757       !$acc kernels loop present( nzb_s_inner, ol, us, uv_total, zu, zw, z0 ) private( j, k, z_mo )
758       DO  i = nxlg, nxrg
759          DO  j = nysg, nyng
760
761             k   = nzb_s_inner(j,i)+1
762             z_mo = zu(k+1) - zw(k)
763
764!
765!--          Compute u* at the scalars' grid points
766             us(j,i) = kappa * uv_total(j,i) / ( LOG( z_mo / z0(j,i) )         &
767                                          - psi_m( z_mo / ol(j,i) )            &
768                                          + psi_m( z0(j,i) / ol(j,i) ) )
769          ENDDO
770       ENDDO
771
772    END SUBROUTINE calc_us
773
774!
775!-- Calculate potential temperature and specific humidity at first grid level
776    SUBROUTINE calc_pt_q
777
778       IMPLICIT NONE
779
780       !$acc kernels loop present( nzb_s_inner, pt, pt1, pt_d_t, q, ql, qv1 ) private( j, k )
781       DO  i = nxlg, nxrg
782          DO  j = nysg, nyng
783             k   = nzb_s_inner(j,i)+1
784             pt1(j,i) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
785             qv1(j,i) = q(k,j,i) - ql(k,j,i)
786          ENDDO
787       ENDDO
788
789    END SUBROUTINE calc_pt_q
790
791!
792!-- Calculate the other MOST scaling parameters theta*, q*, (qr*, nr*)
793    SUBROUTINE calc_scaling_parameters
794
795       IMPLICIT NONE
796
797!
798!--    Data information for accelerators
799       !$acc data present( e, nrsws, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt )  &
800       !$acc      present( q, qs, qsws, qrsws, shf, ts, u, us, usws, v )     &
801       !$acc      present( vpt, vsws, zu, zw, z0, z0h )
802!
803!--    Compute theta*
804       IF ( constant_heatflux )  THEN
805
806!
807!--       For a given heat flux in the surface layer:
808          !$OMP PARALLEL DO
809          !$acc kernels loop private( j )
810          DO  i = nxlg, nxrg
811             DO  j = nysg, nyng
812                ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30_wp )
813!
814!--             ts must be limited, because otherwise overflow may occur in case
815!--             of us=0 when computing ol further below
816                IF ( ts(j,i) < -1.05E5_wp )  ts(j,i) = -1.0E5_wp
817                IF ( ts(j,i) >   1.0E5_wp )  ts(j,i) =  1.0E5_wp
818             ENDDO
819          ENDDO
820
821       ELSE
822!
823!--       For a given surface temperature:
824          IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
825             !$OMP PARALLEL DO
826             !$acc kernels loop private( j, k )
827             DO  i = nxlg, nxrg
828                DO  j = nysg, nyng
829                   k = nzb_s_inner(j,i)
830                   pt(k,j,i) = pt_surface
831                ENDDO
832             ENDDO
833          ENDIF
834
835          !$OMP PARALLEL DO PRIVATE( k, z_mo )
836          !$acc kernels loop present( nzb_s_inner, ol, pt, pt1, ts, zu, zw, z0h ) private( j, k, z_mo )
837          DO  i = nxlg, nxrg
838             DO  j = nysg, nyng
839
840                k   = nzb_s_inner(j,i)
841                z_mo = zu(k+1) - zw(k)
842
843                IF ( cloud_physics )  THEN
844                   ts(j,i) = kappa * ( pt1(j,i) - pt(k,j,i) )                  &
845                                     / ( LOG( z_mo / z0h(j,i) )                &
846                                         - psi_h( z_mo / ol(j,i) )             &
847                                         + psi_h( z0h(j,i) / ol(j,i) ) )
848                ELSE
849                   ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) )               &
850                                     / ( LOG( z_mo / z0h(j,i) )                &
851                                         - psi_h( z_mo / ol(j,i) )             &
852                                         + psi_h( z0h(j,i) / ol(j,i) ) )
853                ENDIF
854
855             ENDDO
856          ENDDO
857       ENDIF
858
859!
860!--    If required compute q*
861       IF ( humidity  .OR.  passive_scalar )  THEN
862          IF ( constant_waterflux )  THEN
863!
864!--          For a given water flux in the surface layer
865             !$OMP PARALLEL DO
866             !$acc kernels loop private( j )
867             DO  i = nxlg, nxrg
868                DO  j = nysg, nyng
869                   qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp )
870                ENDDO
871             ENDDO
872
873          ELSE
874             coupled_run = ( coupling_mode == 'atmosphere_to_ocean'  .AND.     &
875                             run_coupled )
876
877             IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
878                !$OMP PARALLEL DO
879                !$acc kernels loop private( j, k )
880                DO  i = nxlg, nxrg
881                   DO  j = nysg, nyng
882                      k = nzb_s_inner(j,i)
883                      q(k,j,i) = q_surface
884                   ENDDO
885                ENDDO
886             ENDIF
887
888             !$OMP PARALLEL DO PRIVATE( e_s, k, z_mo )
889             !$acc kernels loop independent present( nzb_s_inner, ol, pt, q, qs, qv1, zu, zw, z0q ) private( e_s, j, k, z_mo )
890             DO  i = nxlg, nxrg
891                !$acc loop independent
892                DO  j = nysg, nyng
893
894                   k   = nzb_s_inner(j,i)
895                   z_mo = zu(k+1) - zw(k)
896
897!
898!--                Assume saturation for atmosphere coupled to ocean (but not
899!--                in case of precursor runs)
900                   IF ( coupled_run )  THEN
901                      e_s = 6.1_wp * &
902                              EXP( 0.07_wp * ( MIN(pt(k,j,i),pt(k+1,j,i))      &
903                                               - 273.15_wp ) )
904                      q(k,j,i) = 0.622_wp * e_s / ( surface_pressure - e_s )
905                   ENDIF
906
907                   IF ( cloud_physics )  THEN
908                      qs(j,i) = kappa * ( qv1(j,i) - q(k,j,i) )                &
909                                        / ( LOG( z_mo / z0q(j,i) )             &
910                                            - psi_h( z_mo / ol(j,i) )          &
911                                            + psi_h( z0q(j,i) / ol(j,i) ) )
912
913                   ELSE
914                      qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) )              &
915                                        / ( LOG( z_mo / z0q(j,i) )             &
916                                            - psi_h( z_mo / ol(j,i) )          &
917                                            + psi_h( z0q(j,i) / ol(j,i) ) )
918                   ENDIF
919
920                ENDDO
921             ENDDO
922          ENDIF
923       ENDIF
924
925
926!
927!--    If required compute qr* and nr*
928       IF ( cloud_physics  .AND.  microphysics_seifert )   &
929       THEN
930
931          !$OMP PARALLEL DO PRIVATE( k, z_mo )
932          !$acc kernels loop independent present( nr, nrs, nzb_s_inner, ol, qr, qrs, zu, zw, z0q ) private( j, k, z_mo )
933          DO  i = nxlg, nxrg
934             !$acc loop independent
935             DO  j = nysg, nyng
936
937                k   = nzb_s_inner(j,i)
938                z_mo = zu(k+1) - zw(k)
939
940                qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) )              &
941                                 / ( LOG( z_mo / z0q(j,i) )                 &
942                                     - psi_h( z_mo / ol(j,i) )              &
943                                     + psi_h( z0q(j,i) / ol(j,i) ) )
944
945                nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) )              &
946                                 / ( LOG( z_mo / z0q(j,i) )                 &
947                                     - psi_h( z_mo / ol(j,i) )              &
948                                     + psi_h( z0q(j,i) / ol(j,i) ) )
949             ENDDO
950          ENDDO
951
952       ENDIF
953       !$acc end data
954
955    END SUBROUTINE calc_scaling_parameters
956
957
958
959!
960!-- Calculate surface fluxes usws, vsws, shf, qsws, (qrsws, nrsws)
961    SUBROUTINE calc_surface_fluxes
962
963       IMPLICIT NONE
964
965       REAL(wp) :: ol_mid !< Grid-interpolated L
966
967!
968!--    Compute u'w' for the total model domain.
969!--    First compute the corresponding component of u* and square it.
970       !$OMP PARALLEL DO PRIVATE( k, ol_mid, z_mo )
971       !$acc kernels loop present( nzb_u_inner, ol, u, us, usws, zu, zw, z0 ) private( j, k, z_mo )
972       DO  i = nxl, nxr
973          DO  j = nys, nyn
974
975             k   = nzb_u_inner(j,i)
976             z_mo = zu(k+1) - zw(k)
977!
978!--          Compute bulk Obukhov length for this point
979             ol_mid = 0.5_wp * ( ol(j,i-1) + ol(j,i) )
980
981             IF ( ol_mid == 0.0_wp )  THEN
982                ol_mid = MIN(ol(j,i-1), ol(j,i))
983             ENDIF
984
985             usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) )                     &
986                                 / ( LOG( z_mo / z0(j,i) )                     &
987                                     - psi_m( z_mo / ol_mid )                  &
988                                     + psi_m( z0(j,i) / ol_mid ) )
989
990             usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
991          ENDDO
992       ENDDO
993
994!
995!--    Compute v'w' for the total model domain.
996!--    First compute the corresponding component of u* and square it.
997       !$OMP PARALLEL DO PRIVATE( k, ol_mid, z_mo )
998       !$acc kernels loop present( nzb_v_inner, ol, v, us, vsws, zu, zw, z0 ) private( j, k, ol_mid, z_mo )
999       DO  i = nxl, nxr
1000          DO  j = nys, nyn
1001
1002             k   = nzb_v_inner(j,i)
1003             z_mo = zu(k+1) - zw(k)
1004!
1005!--          Compute bulk Obukhov length for this point
1006             ol_mid = 0.5_wp * ( ol(j-1,i) + ol(j,i) )
1007
1008             IF ( ol_mid == 0.0_wp )  THEN
1009                ol_mid = MIN(ol(j-1,i), ol(j-1,i))
1010             ENDIF
1011
1012             vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) )                     &
1013                                 / ( LOG( z_mo / z0(j,i) )                     &
1014                                     - psi_m( z_mo / ol_mid )                  &
1015                                     + psi_m( z0(j,i) / ol_mid ) )
1016
1017             vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
1018
1019          ENDDO
1020       ENDDO
1021
1022!
1023!--    Exchange the boundaries for the momentum fluxes (is this still required?)
1024       !$acc update host( usws, vsws )
1025       CALL exchange_horiz_2d( usws )
1026       CALL exchange_horiz_2d( vsws )
1027       !$acc update device( usws, vsws )
1028
1029!
1030!--    Compute the vertical kinematic heat flux
1031       IF (  .NOT.  constant_heatflux .AND.  ( simulated_time <=               &
1032            skip_time_do_lsm  .OR.  .NOT.  land_surface ) )  THEN
1033          !$OMP PARALLEL DO
1034          !$acc kernels loop independent present( shf, ts, us )
1035          DO  i = nxlg, nxrg
1036             !$acc loop independent
1037             DO  j = nysg, nyng
1038                shf(j,i) = -ts(j,i) * us(j,i)
1039             ENDDO
1040          ENDDO
1041
1042       ENDIF
1043
1044!
1045!-- Compute the vertical water/scalar flux
1046       IF (  .NOT.  constant_waterflux  .AND.  ( humidity  .OR.                &
1047             passive_scalar )  .AND.  ( simulated_time <= skip_time_do_lsm     &
1048            .OR.  .NOT.  land_surface ) )  THEN
1049          !$OMP PARALLEL DO
1050          !$acc kernels loop independent present( qs, qsws, us )
1051          DO  i = nxlg, nxrg
1052             !$acc loop independent
1053             DO  j = nysg, nyng
1054                qsws(j,i) = -qs(j,i) * us(j,i)
1055             ENDDO
1056          ENDDO
1057
1058       ENDIF
1059
1060!
1061!--    Compute (turbulent) fluxes of rain water content and rain drop conc.
1062       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
1063          !$OMP PARALLEL DO
1064          !$acc kernels loop independent present( nrs, nrsws, qrs, qrsws, us )
1065          DO  i = nxlg, nxrg
1066             !$acc loop independent
1067             DO  j = nysg, nyng
1068                qrsws(j,i) = -qrs(j,i) * us(j,i)
1069                nrsws(j,i) = -nrs(j,i) * us(j,i)
1070             ENDDO
1071          ENDDO
1072       ENDIF
1073
1074!
1075!--    Bottom boundary condition for the TKE
1076       IF ( ibc_e_b == 2 )  THEN
1077          !$OMP PARALLEL DO
1078          !$acc kernels loop independent present( e, nzb_s_inner, us )
1079          DO  i = nxlg, nxrg
1080             !$acc loop independent
1081             DO  j = nysg, nyng
1082                k = nzb_s_inner(j,i)
1083                e(k+1,j,i) = ( us(j,i) / 0.1_wp )**2
1084!
1085!--             As a test: cm = 0.4
1086!               e(k+1,j,i) = ( us(j,i) / 0.4_wp )**2
1087                e(k,j,i)   = e(k+1,j,i)
1088             ENDDO
1089          ENDDO
1090       ENDIF
1091
1092    END SUBROUTINE calc_surface_fluxes
1093
1094
1095!
1096!-- Integrated stability function for momentum
1097    FUNCTION psi_m( zeta ) 
1098       
1099       USE kinds
1100
1101       IMPLICIT NONE
1102
1103       REAL(wp)            :: psi_m !< Integrated similarity function result
1104       REAL(wp)            :: zeta  !< Stability parameter z/L
1105       REAL(wp)            :: x     !< dummy variable
1106
1107       REAL(wp), PARAMETER :: a = 1.0_wp            !< constant
1108       REAL(wp), PARAMETER :: b = 0.66666666666_wp  !< constant
1109       REAL(wp), PARAMETER :: c = 5.0_wp            !< constant
1110       REAL(wp), PARAMETER :: d = 0.35_wp           !< constant
1111       REAL(wp), PARAMETER :: c_d_d = c / d         !< constant
1112       REAL(wp), PARAMETER :: bc_d_d = b * c / d    !< constant
1113
1114
1115       IF ( zeta < 0.0_wp )  THEN
1116          x = SQRT( SQRT( 1.0_wp  - 16.0_wp * zeta ) )
1117          psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2    &
1118                  * ( 1.0_wp + x**2 ) * 0.125_wp )
1119       ELSE
1120
1121          psi_m = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - a * zeta         &
1122                   - bc_d_d
1123!
1124!--       Old version for stable conditions (only valid for z/L < 0.5)
1125!--       psi_m = - 5.0_wp * zeta
1126
1127       ENDIF
1128
1129    END FUNCTION psi_m
1130
1131
1132!
1133!-- Integrated stability function for heat and moisture
1134    FUNCTION psi_h( zeta ) 
1135       
1136       USE kinds
1137
1138       IMPLICIT NONE
1139
1140       REAL(wp)            :: psi_h !< Integrated similarity function result
1141       REAL(wp)            :: zeta  !< Stability parameter z/L
1142       REAL(wp)            :: x     !< dummy variable
1143
1144       REAL(wp), PARAMETER :: a = 1.0_wp            !< constant
1145       REAL(wp), PARAMETER :: b = 0.66666666666_wp  !< constant
1146       REAL(wp), PARAMETER :: c = 5.0_wp            !< constant
1147       REAL(wp), PARAMETER :: d = 0.35_wp           !< constant
1148       REAL(wp), PARAMETER :: c_d_d = c / d         !< constant
1149       REAL(wp), PARAMETER :: bc_d_d = b * c / d    !< constant
1150
1151
1152       IF ( zeta < 0.0_wp )  THEN
1153          x = SQRT( 1.0_wp  - 16.0_wp * zeta )
1154          psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
1155       ELSE
1156          psi_h = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - (1.0_wp          &
1157                  + 0.66666666666_wp * a * zeta )**1.5_wp - bc_d_d             &
1158                  + 1.0_wp
1159!
1160!--       Old version for stable conditions (only valid for z/L < 0.5)
1161!--       psi_h = - 5.0_wp * zeta
1162       ENDIF
1163
1164    END FUNCTION psi_h
1165
1166 END MODULE surface_layer_fluxes_mod
Note: See TracBrowser for help on using the repository browser.