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

Last change on this file since 1921 was 1921, checked in by suehring, 8 years ago

last commit documented

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