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

Last change on this file since 1750 was 1750, checked in by raasch, 8 years ago

last commit documented

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