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

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

some changes in land surface model, radiation model, nudging and some minor updates

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