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

Last change on this file since 1697 was 1697, checked in by raasch, 6 years ago

FORTRAN an OpenMP errors removed
misplaced cpp-directive fixed
small E- and F-FORMAT changes to avoid informative compiler messages about insufficient field width

  • Property svn:keywords set to Id
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! FORTRAN and OpenMP errors removed
22!
23! Former revisions:
24! -----------------
25! $Id: surface_layer_fluxes.f90 1697 2015-10-28 17:14:10Z raasch $
26!
27! 1696 2015-10-27 10:03:34Z maronga
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 in 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 !< Bulk Richardson number
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
398!
399!--    Compute the absolute value of the horizontal velocity (relative to the
400!--    surface). This is required by all methods
401       !$OMP PARALLEL DO PRIVATE( k, z_mo )
402       !$acc kernels loop
403       DO  i = nxl, nxr
404          DO  j = nys, nyn
405
406             k   = nzb_s_inner(j,i)
407
408             uv_total(j,i) = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)      &
409                                         - u(k,j,i)   - u(k,j,i+1) ) )**2 +    &
410                              ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)           &
411                                         - v(k,j,i)   - v(k,j+1,i) ) )**2 )
412
413!
414!--          For too small values of the local wind, MOST does not work. A
415!--          threshold value is thus set if required
416!            uv_total(j,i) = MAX(0.01_wp,uv_total(j,i))
417
418          ENDDO
419       ENDDO
420
421!
422!--    Values of uv_total need to be exchanged at the ghost boundaries
423       !$acc update host( uv_total )
424       CALL exchange_horiz_2d( uv_total )
425       !$acc update device( uv_total )
426
427       IF ( TRIM( most_method ) /= 'circular' )  THEN
428     
429          !$OMP PARALLEL DO PRIVATE( k, z_mo )
430          !$acc kernels loop
431          DO  i = nxl, nxr
432             DO  j = nys, nyn
433
434                k   = nzb_s_inner(j,i)
435                z_mo = zu(k+1) - zw(k)
436
437!
438!--             Evaluate bulk Richardson number (calculation depends on
439!--             definition based on setting of boundary conditions
440                IF ( ibc_pt_b /= 1 )  THEN
441                   IF ( humidity )  THEN
442                      rib(j,i) = g * z_mo * ( vpt(k+1,j,i) - vpt(k,j,i) )      &
443                           / ( uv_total(j,i)**2 * vpt(k+1,j,i) )
444                   ELSE
445                      rib(j,i) = g * z_mo * ( pt(k+1,j,i) - pt(k,j,i) )        &
446                           / ( uv_total(j,i)**2 * pt(k+1,j,i) )
447                   ENDIF     
448                ELSE
449!
450!--                When using Neumann boundary conditions, the buoyancy flux
451!--                is required but cannot be calculated at the surface, as pt
452!--                and q are not known at the surface. Hence the values at
453!--                first grid level are used to estimate the buoyancy flux
454                   IF ( humidity )  THEN
455                      rib(j,i) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
456                                 * q(k+1,j,i) ) * shf(j,i) + 0.61_wp           &
457                                 * pt(k+1,j,i) * qsws(j,i) )   &
458                                 / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2 )
459                   ELSE
460                      rib(j,i) = - g * z_mo * shf(j,i)                         &
461                           / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2 )
462                   ENDIF
463                ENDIF 
464     
465             ENDDO
466          ENDDO
467
468       ENDIF
469
470!
471!--    Calculate the Obukhov length either using a Newton iteration
472!--    method, via a lookup table, or using the old circular way
473       IF ( TRIM( most_method ) == 'newton' )  THEN
474
475          !$OMP PARALLEL DO PRIVATE( k, z_mo )
476          !$acc kernels loop
477          DO  i = nxl, nxr
478             DO  j = nys, nyn
479
480                k   = nzb_s_inner(j,i)
481                z_mo = zu(k+1) - zw(k)
482
483!
484!--             Store current value in case the Newton iteration fails
485                ol_old = ol(j,i)
486
487!
488!--             Ensure that the bulk Richardson number and the Obukhov
489!--             lengtH have the same sign
490                IF ( rib(j,i) * ol(j,i) < 0.0_wp .OR.                        &
491                     ABS( ol(j,i) ) == ol_max )  THEN
492                   IF ( rib(j,i) > 0.0_wp ) ol(j,i) =  0.01_wp
493                   IF ( rib(j,i) < 0.0_wp ) ol(j,i) = -0.01_wp
494                ENDIF
495!
496!--             Iteration to find Obukhov length
497                iter = 0
498                DO
499                   iter = iter + 1
500!
501!--                In case of divergence, use the value of the previous time step
502                   IF ( iter > 1000 )  THEN
503                      ol(j,i) = ol_old
504                      EXIT
505                   ENDIF
506
507                   ol_m = ol(j,i)
508                   ol_l = ol_m - 0.001_wp * ol_m
509                   ol_u = ol_m + 0.001_wp * ol_m
510
511
512                   IF ( ibc_pt_b /= 1 )  THEN
513!
514!--                   Calculate f = Ri - [...]/[...]^2 = 0
515                      f = rib(j,i) - ( z_mo / ol_m ) * ( LOG( z_mo / z0h(j,i) )&
516                                                    - psi_h( z_mo / ol_m )     &
517                                                    + psi_h( z0h(j,i) / ol_m ) &
518                                                   )                           &
519                                                 / ( LOG( z_mo / z0(j,i) )     &
520                                                    - psi_m( z_mo / ol_m )     &
521                                                    + psi_m( z0(j,i) / ol_m )  &
522                                                    )**2
523
524!
525!--                    Calculate df/dL
526                       f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / z0h(j,i) ) &
527                                                   - psi_h( z_mo / ol_u )      &
528                                                   + psi_h( z0h(j,i) / ol_u )  &
529                                                 )                             &
530                                               / ( LOG( z_mo / z0(j,i) )       &
531                                                   - psi_m( z_mo / ol_u )      &
532                                                   + psi_m( z0(j,i) / ol_u )   &
533                                                 )**2                          &
534                              + ( z_mo / ol_l ) * ( LOG( z_mo / z0h(j,i) )     &
535                                                   - psi_h( z_mo / ol_l )      &
536                                                   + psi_h( z0h(j,i) / ol_l )  &
537                                                 )                             &
538                                               / ( LOG( z_mo / z0(j,i) )       &
539                                                   - psi_m( z_mo / ol_l )      &
540                                                   + psi_m( z0(j,i) / ol_l )   &
541                                                 )**2                          &
542                                ) / ( ol_u - ol_l )
543                   ELSE
544!
545!--                   Calculate f = Ri - 1 /[...]^3 = 0
546                      f = rib(j,i) - ( z_mo / ol_m ) / ( LOG( z_mo / z0(j,i) )&
547                                                    - psi_m( z_mo / ol_m )    &
548                                                    + psi_m( z0(j,i) / ol_m ) &
549                                                       )**3
550
551!
552!--                   Calculate df/dL
553                      f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / z0(j,i) )  &
554                                                   - psi_m( z_mo / ol_u )     &
555                                                   + psi_m( z0(j,i) / ol_u )  &
556                                                 )**3                         &
557                              + ( z_mo / ol_l ) / ( LOG( z_mo / z0(j,i) )     &
558                                                   - psi_m( z_mo / ol_l )     &
559                                                   + psi_m( z0(j,i) / ol_l )  &
560                                                 )**3                         &
561                                     ) / ( ol_u - ol_l )
562                   ENDIF
563!
564!--                Calculate new L
565                   ol(j,i) = ol_m - f / f_d_ol
566
567!
568!--                Ensure that the bulk Richardson number and the Obukhov
569!--                length have the same sign and ensure convergence.
570                   IF ( ol(j,i) * ol_m < 0.0_wp )  ol(j,i) = ol_m * 0.5_wp
571
572!
573!--                If unrealistic value occurs, set L to the maximum
574!--                value that is allowed
575                   IF ( ABS( ol(j,i) ) > ol_max )  THEN
576                      ol(j,i) = ol_max
577                      EXIT
578                   ENDIF
579!
580!--                Check for convergence
581                   IF ( ABS( ( ol(j,i) - ol_m ) / ol(j,i) ) < 1.0E-4_wp )  THEN
582                      EXIT
583                   ELSE
584                      CYCLE
585                   ENDIF
586
587                ENDDO
588                       
589             ENDDO
590          ENDDO
591
592       ELSEIF ( TRIM( most_method ) == 'lookup' )  THEN
593
594          !$OMP PARALLEL DO PRIVATE( k, z_mo )
595          !$acc kernels loop
596          DO  i = nxl, nxr
597             DO  j = nys, nyn
598
599                k   = nzb_s_inner(j,i)
600
601!
602!--             If the bulk Richardson number is outside the range of the lookup
603!--             table, set it to the exceeding threshold value
604                IF ( rib(j,i) < rib_min )  rib(j,i) = rib_min
605                IF ( rib(j,i) > rib_max )  rib(j,i) = rib_max
606
607!
608!--             Find the correct index bounds for linear interpolation. As the
609!--             Richardson number will not differ very much from time step to
610!--             time step , use the index from the last step and search in the
611!--             correct direction
612                l = l_bnd
613                IF ( rib_tab(l) - rib(j,i) > 0.0_wp )  THEN
614                   DO WHILE ( rib_tab(l-1) - rib(j,i) > 0.0_wp .AND. l > 0 )
615                      l = l-1
616                   ENDDO
617                ELSE
618                   DO WHILE ( rib_tab(l) - rib(j,i) < 0.0_wp                &
619                              .AND. l < num_steps-1 )
620                      l = l+1
621                   ENDDO
622                ENDIF
623                l_bnd = l
624
625!
626!--             Linear interpolation to find the correct value of z/L
627                ol(j,i) = ( ol_tab(l-1) + ( ol_tab(l) - ol_tab(l-1) )       &
628                            / (  rib_tab(l) - rib_tab(l-1) )                &
629                            * ( rib(j,i) - rib_tab(l-1) ) )
630
631             ENDDO
632          ENDDO
633
634       ELSEIF ( TRIM( most_method ) == 'circular' )  THEN
635
636          !$OMP PARALLEL DO PRIVATE( k, z_mo )
637          !$acc kernels loop
638          DO  i = nxl, nxr
639             DO  j = nys, nyn
640
641                k   = nzb_s_inner(j,i)
642                z_mo = zu(k+1) - zw(k)
643
644                IF ( .NOT. humidity )  THEN
645                   ol(j,i) =  ( pt(k+1,j,i) *  us(j,i)**2 ) / ( kappa * g      &
646                              * ts(j,i) + 1E-30_wp )
647                ELSEIF ( cloud_physics )  THEN
648
649                   ol(j,i) =  ( vpt(k+1,j,i) * us(j,i)**2 ) / ( kappa * g      &
650                              * ( ts(j,i) + 0.61_wp * pt1(j,i) * qs(j,i)       &
651                              + 0.61_wp * qv1(j,i) * ts(j,i) - ts(j,i)         &
652                              * ql(k+1,j,i) ) + 1E-30_wp )
653                ELSE
654                   ol(j,i) =  ( vpt(k+1,j,i) *  us(j,i)**2 ) / ( kappa * g     &
655                              * ( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i)    &
656                                  + 0.61_wp * q(k+1,j,i) * ts(j,i) ) + 1E-30_wp )
657                ENDIF
658!
659!--             Limit the value range of the Obukhov length.
660!--             This is necessary for very small velocities (u,v --> 0), because
661!--             the absolute value of ol can then become very small, which in
662!--             consequence would result in very large shear stresses and very
663!--             small momentum fluxes (both are generally unrealistic).
664                IF ( ( z_mo / ol(j,i) ) < zeta_min )  ol(j,i) = z_mo / zeta_min
665                IF ( ( z_mo / ol(j,i) ) > zeta_max )  ol(j,i) = z_mo / zeta_max
666
667             ENDDO
668          ENDDO
669
670       ENDIF
671
672!
673!--    Values of ol at ghost point locations are needed for the evaluation
674!--    of usws and vsws.
675       !$acc update host( ol )
676       CALL exchange_horiz_2d( ol )
677       !$acc update device( ol )
678
679    END SUBROUTINE calc_ol
680
681!
682!-- Calculate friction velocity u*
683    SUBROUTINE calc_us
684
685       IMPLICIT NONE
686
687       !$OMP PARALLEL DO PRIVATE( k, z_mo )
688       !$acc kernels loop
689       DO  i = nxlg, nxrg
690          DO  j = nysg, nyng
691
692             k   = nzb_s_inner(j,i)+1
693             z_mo = zu(k+1) - zw(k)
694
695!
696!--          Compute u* at the scalars' grid points
697             us(j,i) = kappa * uv_total(j,i) / ( LOG( z_mo / z0(j,i) )         &
698                                          - psi_m( z_mo / ol(j,i) )            &
699                                          + psi_m( z0(j,i) / ol(j,i) ) )
700          ENDDO
701       ENDDO
702
703    END SUBROUTINE calc_us
704
705!
706!-- Calculate potential temperature and specific humidity at first grid level
707    SUBROUTINE calc_pt_q
708
709       IMPLICIT NONE
710
711       DO  i = nxlg, nxrg
712          DO  j = nysg, nyng
713             k   = nzb_s_inner(j,i)+1
714             pt1(j,i) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
715             qv1(j,i) = q(k,j,i) - ql(k,j,i)
716          ENDDO
717       ENDDO
718
719    END SUBROUTINE calc_pt_q
720
721!
722!-- Calculate the other MOST scaling parameters theta*, q*, (qr*, nr*)
723    SUBROUTINE calc_scaling_parameters
724
725       IMPLICIT NONE
726
727!
728!--    Data information for accelerators
729       !$acc data present( e, nrsws, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt )  &
730       !$acc      present( q, qs, qsws, qrsws, shf, ts, u, us, usws, v )     &
731       !$acc      present( vpt, vsws, zu, zw, z0, z0h )
732!
733!--    Compute theta*
734       IF ( constant_heatflux )  THEN
735
736!
737!--       For a given heat flux in the surface layer:
738          !$OMP PARALLEL DO
739          !$acc kernels loop
740          DO  i = nxlg, nxrg
741             DO  j = nysg, nyng
742                ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30_wp )
743!
744!--             ts must be limited, because otherwise overflow may occur in case
745!--             of us=0 when computing ol further below
746                IF ( ts(j,i) < -1.05E5_wp )  ts(j,i) = -1.0E5_wp
747                IF ( ts(j,i) >   1.0E5_wp )  ts(j,i) =  1.0E5_wp
748             ENDDO
749          ENDDO
750
751       ELSE
752!
753!--       For a given surface temperature:
754          IF ( large_scale_forcing .AND. lsf_surf )  THEN
755             !$OMP PARALLEL DO
756             !$acc kernels loop
757             DO  i = nxlg, nxrg
758                DO  j = nysg, nyng
759                   k = nzb_s_inner(j,i)
760                   pt(k,j,i) = pt_surface
761                ENDDO
762             ENDDO
763          ENDIF
764
765          !$OMP PARALLEL DO PRIVATE( k, z_mo )
766          !$acc kernels loop
767          DO  i = nxlg, nxrg
768             DO  j = nysg, nyng
769
770                k   = nzb_s_inner(j,i)
771                z_mo = zu(k+1) - zw(k)
772
773                IF ( cloud_physics )  THEN
774                   ts(j,i) = kappa * ( pt1(j,i) - pt(k,j,i) )                  &
775                                     / ( LOG( z_mo / z0h(j,i) )                &
776                                         - psi_h( z_mo / ol(j,i) )             &
777                                         + psi_h( z0h(j,i) / ol(j,i) ) )
778                ELSE
779                   ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) )               &
780                                     / ( LOG( z_mo / z0h(j,i) )                &
781                                         - psi_h( z_mo / ol(j,i) )             &
782                                         + psi_h( z0h(j,i) / ol(j,i) ) )
783                ENDIF
784
785             ENDDO
786          ENDDO
787       ENDIF
788
789!
790!--    If required compute q*
791       IF ( humidity  .OR.  passive_scalar )  THEN
792          IF ( constant_waterflux )  THEN
793!
794!--          For a given water flux in the Prandtl layer:
795             !$OMP PARALLEL DO
796             !$acc kernels loop
797             DO  i = nxlg, nxrg
798                DO  j = nysg, nyng
799                   qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp )
800                ENDDO
801             ENDDO
802
803          ELSE
804             coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND.      &
805                             run_coupled )
806
807             IF ( large_scale_forcing .AND. lsf_surf )  THEN
808                !$OMP PARALLEL DO
809                !$acc kernels loop
810                DO  i = nxlg, nxrg
811                   DO  j = nysg, nyng
812                      k = nzb_s_inner(j,i)
813                      q(k,j,i) = q_surface
814                   ENDDO
815                ENDDO
816             ENDIF
817
818             !$OMP PARALLEL DO PRIVATE( e_s, k, z_mo )
819             !$acc kernels loop independent
820             DO  i = nxlg, nxrg
821                !$acc loop independent
822                DO  j = nysg, nyng
823
824                   k   = nzb_s_inner(j,i)
825                   z_mo = zu(k+1) - zw(k)
826
827!
828!--                Assume saturation for atmosphere coupled to ocean (but not
829!--                in case of precursor runs)
830                   IF ( coupled_run )  THEN
831                      e_s = 6.1_wp * &
832                              EXP( 0.07_wp * ( MIN(pt(k,j,i),pt(k+1,j,i))      &
833                                               - 273.15_wp ) )
834                      q(k,j,i) = 0.622_wp * e_s / ( surface_pressure - e_s )
835                   ENDIF
836
837                   IF ( cloud_physics )  THEN
838                      qs(j,i) = kappa * ( qv1(j,i) - q(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
843                   ELSE
844                      qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) )              &
845                                        / ( LOG( z_mo / z0h(j,i) )             &
846                                            - psi_h( z_mo / ol(j,i) )          &
847                                            + psi_h( z0h(j,i) / ol(j,i) ) )
848                   ENDIF
849
850                ENDDO
851             ENDDO
852          ENDIF
853       ENDIF
854
855
856!
857!--    If required compute qr* and nr*
858       IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation )  THEN
859
860          !$OMP PARALLEL DO PRIVATE( k, z_mo )
861          !$acc kernels loop independent
862          DO  i = nxlg, nxrg
863             !$acc loop independent
864             DO  j = nysg, nyng
865
866                k   = nzb_s_inner(j,i)
867                z_mo = zu(k+1) - zw(k)
868
869                   qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) )              &
870                                     / ( LOG( z_mo / z0h(j,i) )                &
871                                         - psi_h( z_mo / ol(j,i) )             &
872                                         + psi_h( z0h(j,i) / ol(j,i) ) )
873
874                   nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) )              &
875                                     / ( LOG( z_mo / z0h(j,i) )                &
876                                         - psi_h( z_mo / ol(j,i) )             &
877                                         + psi_h( z0h(j,i) / ol(j,i) ) )
878             ENDDO
879          ENDDO
880
881       ENDIF
882
883    END SUBROUTINE calc_scaling_parameters
884
885
886
887!
888!-- Calculate surface fluxes usws, vsws, shf, qsws, (qrsws, nrsws)
889    SUBROUTINE calc_surface_fluxes
890
891       IMPLICIT NONE
892
893       REAL(wp) :: ol_mid !< Grid-interpolated L
894
895!
896!--    Compute u'w' for the total model domain.
897!--    First compute the corresponding component of u* and square it.
898       !$OMP PARALLEL DO PRIVATE( k, ol_mid, z_mo )
899       !$acc kernels loop
900       DO  i = nxl, nxr
901          DO  j = nys, nyn
902
903             k   = nzb_u_inner(j,i)
904             z_mo = zu(k+1) - zw(k)
905!
906!--          Compute bulk Obukhov length for this point
907             ol_mid = 0.5_wp * ( ol(j,i-1) + ol(j,i) )
908
909             IF ( ol_mid == 0.0_wp )  THEN
910                ol_mid = MIN(ol(j,i-1), ol(j,i))
911             ENDIF
912
913             usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) )                     &
914                                 / ( LOG( z_mo / z0(j,i) )                     &
915                                     - psi_m( z_mo / ol_mid )                  &
916                                     + psi_m( z0(j,i) / ol_mid ) )
917
918             usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
919          ENDDO
920       ENDDO
921
922!
923!--    Compute v'w' for the total model domain.
924!--    First compute the corresponding component of u* and square it.
925       !$OMP PARALLEL DO PRIVATE( k, ol_mid, z_mo )
926       !$acc kernels loop
927       DO  i = nxl, nxr
928          DO  j = nys, nyn
929
930             k   = nzb_v_inner(j,i)
931             z_mo = zu(k+1) - zw(k)
932!
933!--          Compute bulk Obukhov length for this point
934             ol_mid = 0.5_wp * ( ol(j-1,i) + ol(j,i) )
935
936             IF ( ol_mid == 0.0_wp )  THEN
937                ol_mid = MIN(ol(j-1,i), ol(j-1,i))
938             ENDIF
939
940             vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) )                     &
941                                 / ( LOG( z_mo / z0(j,i) )                     &
942                                     - psi_m( z_mo / ol_mid )                  &
943                                     + psi_m( z0(j,i) / ol_mid ) )
944
945             vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
946
947          ENDDO
948       ENDDO
949
950!
951!--    Exchange the boundaries for the momentum fluxes (is this still required?)
952       !$acc update host( usws, vsws )
953       CALL exchange_horiz_2d( usws )
954       CALL exchange_horiz_2d( vsws )
955       !$acc update device( usws, vsws )
956
957!
958!--    Compute the vertical kinematic heat flux
959       IF ( .NOT. constant_heatflux .AND. ( simulated_time <=                &
960            skip_time_do_lsm .OR. .NOT. land_surface ) )  THEN
961          !$OMP PARALLEL DO
962          !$acc kernels loop independent
963          DO  i = nxlg, nxrg
964             !$acc loop independent
965             DO  j = nysg, nyng
966                shf(j,i) = -ts(j,i) * us(j,i)
967             ENDDO
968          ENDDO
969
970       ENDIF
971
972!
973!-- Compute the vertical water/scalar flux
974       IF ( .NOT. constant_waterflux .AND. ( humidity .OR. passive_scalar )    &
975            .AND. ( simulated_time <= skip_time_do_lsm .OR. .NOT.            &
976            land_surface ) )  THEN
977          !$OMP PARALLEL DO
978          !$acc kernels loop independent
979          DO  i = nxlg, nxrg
980             !$acc loop independent
981             DO  j = nysg, nyng
982                qsws(j,i) = -qs(j,i) * us(j,i)
983             ENDDO
984          ENDDO
985
986       ENDIF
987
988!
989!--    Compute (turbulent) fluxes of rain water content and rain drop conc.
990       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                    &
991            precipitation )  THEN
992          !$OMP PARALLEL DO
993          !$acc kernels loop independent
994          DO  i = nxlg, nxrg
995             !$acc loop independent
996             DO  j = nysg, nyng
997                qrsws(j,i) = -qrs(j,i) * us(j,i)
998                nrsws(j,i) = -nrs(j,i) * us(j,i)
999             ENDDO
1000          ENDDO
1001       ENDIF
1002
1003!
1004!--    Bottom boundary condition for the TKE
1005       IF ( ibc_e_b == 2 )  THEN
1006          !$OMP PARALLEL DO
1007          !$acc kernels loop independent
1008          DO  i = nxlg, nxrg
1009             !$acc loop independent
1010             DO  j = nysg, nyng
1011                e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1_wp )**2
1012!
1013!--             As a test: cm = 0.4
1014!               e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4_wp )**2
1015                e(nzb_s_inner(j,i),j,i)   = e(nzb_s_inner(j,i)+1,j,i)
1016             ENDDO
1017          ENDDO
1018       ENDIF
1019       !$acc end data
1020
1021    END SUBROUTINE calc_surface_fluxes
1022
1023
1024!
1025!-- Integrated stability function for momentum
1026    FUNCTION psi_m( zeta ) 
1027       
1028       USE kinds
1029
1030       IMPLICIT NONE
1031
1032       REAL(wp)            :: psi_m !< Integrated similarity function result
1033       REAL(wp)            :: zeta  !< Stability parameter z/L
1034       REAL(wp)            :: x     !< dummy variable
1035
1036       REAL(wp), PARAMETER :: a = 1.0_wp            !< constant
1037       REAL(wp), PARAMETER :: b = 0.66666666666_wp  !< constant
1038       REAL(wp), PARAMETER :: c = 5.0_wp            !< constant
1039       REAL(wp), PARAMETER :: d = 0.35_wp           !< constant
1040       REAL(wp), PARAMETER :: c_d_d = c / d         !< constant
1041       REAL(wp), PARAMETER :: bc_d_d = b * c / d    !< constant
1042
1043
1044       IF ( zeta < 0.0_wp )  THEN
1045          x = SQRT( SQRT(1.0_wp  - 16.0_wp * zeta ) )
1046          psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2    &
1047                  * ( 1.0_wp + x**2 ) * 0.125_wp )
1048       ELSE
1049
1050          psi_m = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - a * zeta         &
1051                   - bc_d_d
1052!
1053!--       Old version for stable conditions (only valid for z/L < 0.5)
1054!--       psi_m = - 5.0_wp * zeta
1055
1056       ENDIF
1057
1058    END FUNCTION psi_m
1059
1060
1061!
1062!-- Integrated stability function for heat and moisture
1063    FUNCTION psi_h( zeta ) 
1064       
1065       USE kinds
1066
1067       IMPLICIT NONE
1068
1069       REAL(wp)            :: psi_h !< Integrated similarity function result
1070       REAL(wp)            :: zeta  !< Stability parameter z/L
1071       REAL(wp)            :: x     !< dummy variable
1072
1073       REAL(wp), PARAMETER :: a = 1.0_wp            !< constant
1074       REAL(wp), PARAMETER :: b = 0.66666666666_wp  !< constant
1075       REAL(wp), PARAMETER :: c = 5.0_wp            !< constant
1076       REAL(wp), PARAMETER :: d = 0.35_wp           !< constant
1077       REAL(wp), PARAMETER :: c_d_d = c / d         !< constant
1078       REAL(wp), PARAMETER :: bc_d_d = b * c / d    !< constant
1079
1080
1081       IF ( zeta < 0.0_wp )  THEN
1082          x = SQRT(1.0_wp  - 16.0_wp * zeta )
1083          psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
1084       ELSE
1085          psi_h = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - (1.0_wp          &
1086                  + 0.66666666666_wp * a * zeta )**1.5_wp - bc_d_d             &
1087                  + 1.0_wp
1088!
1089!--       Old version for stable conditions (only valid for z/L < 0.5)
1090!--       psi_h = - 5.0_wp * zeta
1091       ENDIF
1092
1093    END FUNCTION psi_h
1094
1095 END MODULE surface_layer_fluxes_mod
Note: See TracBrowser for help on using the repository browser.