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

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

added support for water and paved surfaced in land surface model / minor changes

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