Home

Context Navigation

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

Last change on this file since 1777 was 1758, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 42.3 KB

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