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

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

last commit documented / copyright update

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