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

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

several bugfixes in particle model and serial mode

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