source: palm/trunk/SOURCE/surface_layer_fluxes.f90 @ 1692

Last change on this file since 1692 was 1692, checked in by maronga, 6 years ago

last commit documented

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