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

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

all OpenACC directives and related parts removed from the code

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