Home

Context Navigation

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

Last change on this file since 2169 was 2119, checked in by raasch, 8 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 suehring $
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.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |