source: palm/trunk/SOURCE/surface_layer_fluxes_mod.f90 @ 1916

Last change on this file since 1916 was 1916, checked in by suehring, 5 years ago

last commit documented

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