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

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

last commit documented

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