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

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

last commit documented

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