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

Last change on this file since 2119 was 2119, checked in by raasch, 7 years ago

last commit documented

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