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

Last change on this file since 1851 was 1851, checked in by maronga, 5 years ago

last commit documented

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