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

Last change on this file since 1823 was 1823, checked in by hoffmann, 8 years ago

last commit documented

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