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

Last change on this file since 1691 was 1691, checked in by maronga, 8 years ago

various bugfixes and modifications of the atmosphere-land-surface-radiation interaction. Completely re-written routine to calculate surface fluxes (surface_layer_fluxes.f90) that replaces prandtl_fluxes. Minor formatting corrections and renamings

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