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

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

bugfix in write_var_list (restart runs not possible)

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