source: palm/trunk/SOURCE/surface_layer_fluxes_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 8 years ago

added _mod string to several filenames to meet the naming convection for modules

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