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

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

several bugfixes related to the new surface layer routine and land-surface-radiation interaction

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