Home

Context Navigation

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

Last change on this file since 1754 was 1750, checked in by raasch, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 42.4 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 1750 2016-02-09 12:24:00Z maronga $
26	!
27	! 1749 2016-02-09 12:19:56Z raasch
28	! further OpenACC adjustments
29	!
30	! 1747 2016-02-08 12:25:53Z raasch
31	! adjustments for OpenACC usage
32	!
33	! 1709 2015-11-04 14:47:01Z maronga
34	! Bugfix: division by zero could occur when calculating rib at low wind speeds
35	! Bugfix: calculation of uv_total for neutral = .T., initial value for ol for
36	! neutral = .T.
37	!
38	! 1705 2015-11-02 14:28:56Z maronga
39	! Typo removed
40	!
41	! 1697 2015-10-28 17:14:10Z raasch
42	! FORTRAN and OpenMP errors removed
43	!
44	! 1696 2015-10-27 10:03:34Z maronga
45	! Modularized and completely re-written version of prandtl_fluxes.f90. In the
46	! course of the re-writing two additional methods have been implemented. See
47	! updated description.
48	!
49	! 1551 2015-03-03 14:18:16Z maronga
50	! Removed land surface model part. The surface fluxes are now always calculated
51	! within prandtl_fluxes, based on the given surface temperature/humidity (which
52	! is either provided by the land surface model, by large scale forcing data, or
53	! directly prescribed by the user.
54	!
55	! 1496 2014-12-02 17:25:50Z maronga
56	! Adapted for land surface model
57	!
58	! 1494 2014-11-21 17:14:03Z maronga
59	! Bugfixes: qs is now calculated before calculation of Rif. calculation of
60	! buoyancy flux in Rif corrected (added missing humidity term), allow use of
61	! topography for coupled runs (not tested)
62	!
63	! 1361 2014-04-16 15:17:48Z hoffmann
64	! Bugfix: calculation of turbulent fluxes of rain water content (qrsws) and rain
65	! drop concentration (nrsws) added
66	!
67	! 1340 2014-03-25 19:45:13Z kanani
68	! REAL constants defined as wp-kind
69	!
70	! 1320 2014-03-20 08:40:49Z raasch
71	! ONLY-attribute added to USE-statements,
72	! kind-parameters added to all INTEGER and REAL declaration statements,
73	! kinds are defined in new module kinds,
74	! old module precision_kind is removed,
75	! revision history before 2012 removed,
76	! comment fields (!:) to be used for variable explanations added to
77	! all variable declaration statements
78	!
79	! 1276 2014-01-15 13:40:41Z heinze
80	! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
81	!
82	! 1257 2013-11-08 15:18:40Z raasch
83	! openACC "kernels do" replaced by "kernels loop", "loop independent" added
84	!
85	! 1036 2012-10-22 13:43:42Z raasch
86	! code put under GPL (PALM 3.9)
87	!
88	! 1015 2012-09-27 09:23:24Z raasch
89	! OpenACC statements added
90	!
91	! 978 2012-08-09 08:28:32Z fricke
92	! roughness length for scalar quantities z0h added
93	!
94	! Revision 1.1 1998/01/23 10:06:06 raasch
95	! Initial revision
96	!
97	!
98	! Description:
99	! ------------
100	!> Diagnostic computation of vertical fluxes in the constant flux layer from the
101	!> values of the variables at grid point k=1. Three different methods are
102	!> available:
103	!> 1) the "old" version (most_method = 'circular') which is fast, but inaccurate
104	!> 2) a Newton iteration method (most_method = 'newton'), which is accurate, but
105	!> slower
106	!> 3) a method using a lookup table which is fast and accurate. Note, however,
107	!> that this method cannot be used in case of roughness heterogeneity
108	!>
109	!> @todo (re)move large_scale_forcing actions
110	!> @todo check/optimize OpenMP and OpenACC directives
111	!------------------------------------------------------------------------------!
112	MODULE surface_layer_fluxes_mod
113
114	USE arrays_3d, &
115	ONLY: e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, &
116	shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h
117
118	USE cloud_parameters, &
119	ONLY: l_d_cp, pt_d_t
120
121	USE constants, &
122	ONLY: pi
123
124	USE cpulog
125
126	USE control_parameters, &
127	ONLY: cloud_physics, constant_heatflux, constant_waterflux, &
128	coupling_mode, g, humidity, ibc_e_b, ibc_pt_b, icloud_scheme, &
129	initializing_actions, kappa, intermediate_timestep_count, &
130	intermediate_timestep_count_max, large_scale_forcing, lsf_surf, &
131	message_string, most_method, neutral, passive_scalar, &
132	precipitation, pt_surface, q_surface, run_coupled, &
133	surface_pressure, simulated_time, terminate_run, zeta_max, &
134	zeta_min
135
136	USE indices, &
137	ONLY: nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_s_inner, &
138	nzb_u_inner, nzb_v_inner
139
140	USE kinds
141
142	USE pegrid
143
144	USE land_surface_model_mod, &
145	ONLY: land_surface, skip_time_do_lsm
146
147	IMPLICIT NONE
148
149	INTEGER(iwp) :: i !< loop index x direction
150	INTEGER(iwp) :: j !< loop index y direction
151	INTEGER(iwp) :: k !< loop index z direction
152
153	INTEGER(iwp), PARAMETER :: num_steps = 15000 !< number of steps in the lookup table
154
155	LOGICAL :: coupled_run !< Flag for coupled atmosphere-ocean runs
156
157	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt1, & !< Potential temperature at first grid level (required for cloud_physics = .T.)
158	qv1, & !< Specific humidity at first grid level (required for cloud_physics = .T.)
159	uv_total !< Total velocity at first grid level
160
161	REAL(wp), DIMENSION(0:num_steps-1) :: rib_tab, & !< Lookup table bulk Richardson number
162	ol_tab !< Lookup table values of L
163
164	REAL(wp) :: e_s, & !< Saturation water vapor pressure
165	l_bnd = 7500, & !< Lookup table index of the last time step
166	ol_max = 1.0E6_wp, & !< Maximum Obukhov length
167	rib_max, & !< Maximum Richardson number in lookup table
168	rib_min, & !< Minimum Richardson number in lookup table
169	z_mo !< Height of the constant flux layer where MOST is assumed
170
171
172	SAVE
173
174	PRIVATE
175
176	PUBLIC init_surface_layer_fluxes, pt1, qv1, surface_layer_fluxes, uv_total
177
178	INTERFACE init_surface_layer_fluxes
179	MODULE PROCEDURE init_surface_layer_fluxes
180	END INTERFACE init_surface_layer_fluxes
181
182	INTERFACE surface_layer_fluxes
183	MODULE PROCEDURE surface_layer_fluxes
184	END INTERFACE surface_layer_fluxes
185
186
187	CONTAINS
188
189
190	!------------------------------------------------------------------------------!
191	! Description:
192	! ------------
193	!> Main routine to compute the surface fluxes
194	!------------------------------------------------------------------------------!
195	SUBROUTINE surface_layer_fluxes
196
197	IMPLICIT NONE
198
199	!
200	!-- In case cloud physics is used, it is required to derive potential
201	!-- temperature and specific humidity at first grid level from the fields pt
202	!-- and q
203	IF ( cloud_physics ) THEN
204	CALL calc_pt_q
205	ENDIF
206
207	!
208	!-- First, calculate the new Obukhov length, then new friction velocity,
209	!-- followed by the new scaling parameters (th, q, etc.), and the new
210	!-- surface fluxes if required. The old routine ("circular") requires a
211	!-- different order of calls as the scaling parameters from the previous time
212	!-- steps are used to calculate the Obukhov length
213
214	!
215	!-- Depending on setting of most_method use the "old" routine
216	IF ( most_method == 'circular' ) THEN
217
218	CALL calc_scaling_parameters
219
220	CALL calc_uv_total
221
222	IF ( .NOT. neutral ) THEN
223	CALL calc_ol
224	ENDIF
225
226	CALL calc_us
227
228	CALL calc_surface_fluxes
229
230	!
231	!-- Use either Newton iteration or a lookup table for the bulk Richardson
232	!-- number to calculate the Obukhov length
233	ELSEIF ( most_method == 'newton' .OR. most_method == 'lookup' ) THEN
234
235	CALL calc_uv_total
236
237	IF ( .NOT. neutral ) THEN
238	CALL calc_ol
239	ENDIF
240
241	CALL calc_us
242
243	CALL calc_scaling_parameters
244
245	CALL calc_surface_fluxes
246
247	ENDIF
248
249	END SUBROUTINE surface_layer_fluxes
250
251
252	!------------------------------------------------------------------------------!
253	! Description:
254	! ------------
255	!> Initializing actions for the surface layer routine. Basically, this involves
256	!> the preparation of a lookup table for the the bulk Richardson number vs
257	!> Obukhov length L when using the lookup table method.
258	!------------------------------------------------------------------------------!
259	SUBROUTINE init_surface_layer_fluxes
260
261	IMPLICIT NONE
262
263	INTEGER(iwp) :: l, & !< Index for loop to create lookup table
264	num_steps_n !< Number of non-stretched zeta steps
265
266	LOGICAL :: terminate_run_l = .FALSE. !< Flag to terminate run (global)
267
268	REAL(wp), PARAMETER :: zeta_stretch = -10.0_wp !< Start of stretching in the free convection limit
269
270	REAL(wp), DIMENSION(:), ALLOCATABLE :: zeta_tmp
271
272
273	REAL(wp) :: zeta_step, & !< Increment of zeta
274	regr = 1.01_wp, & !< Stretching factor of zeta_step in the free convection limit
275	regr_old = 1.0E9_wp, & !< Stretching factor of last iteration step
276	z0h_min = 0.0_wp, & !< Minimum value of z0h to create table
277	z0_min = 0.0_wp !< Minimum value of z0 to create table
278	!
279	!-- When cloud physics is used, arrays for storing potential temperature and
280	!-- specific humidity at first grid level are required
281	IF ( cloud_physics ) THEN
282	ALLOCATE ( pt1(nysg:nyng,nxlg:nxrg) )
283	ALLOCATE ( qv1(nysg:nyng,nxlg:nxrg) )
284	ENDIF
285
286	!
287	!-- Allocate field for storing the horizontal velocity
288	ALLOCATE ( uv_total(nysg:nyng,nxlg:nxrg) )
289
290
291	!
292	!-- In case of runs with neutral statification, set Obukhov length to a
293	!-- large value
294	IF ( neutral ) ol = 1.0E10_wp
295
296	IF ( most_method == 'lookup' ) THEN
297
298	!
299	!-- Check for roughness heterogeneity. In that case terminate run and
300	!-- inform user
301	IF ( MINVAL( z0h ) /= MAXVAL( z0h ) .OR. &
302	MINVAL( z0 ) /= MAXVAL( z0 ) ) THEN
303	terminate_run_l = .TRUE.
304	ENDIF
305
306	#if defined( __parallel )
307	!
308	!-- Make a logical OR for all processes. Force termiation of model if result
309	!-- is TRUE
310	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
311	CALL MPI_ALLREDUCE( terminate_run_l, terminate_run, 1, MPI_LOGICAL, &
312	MPI_LOR, comm2d, ierr )
313	#else
314	terminate_run = terminate_run_l
315	#endif
316
317	IF ( terminate_run ) THEN
318	message_string = 'most_method = "lookup" cannot be used in ' // &
319	'combination with a prescribed roughness ' // &
320	'heterogeneity'
321	CALL message( 'surface_layer_fluxes', 'PA0417', 1, 2, 0, 6, 0 )
322	ENDIF
323
324	ALLOCATE( zeta_tmp(0:num_steps-1) )
325
326	!
327	!-- Use the lowest possible value for z_mo
328	k = MINVAL(nzb_s_inner)
329	z_mo = zu(k+1) - zw(k)
330
331	!
332	!-- Calculate z/L range from zeta_stretch to zeta_max using 90% of the
333	!-- available steps (num_steps). The calculation is done with negative
334	!-- values of zeta in order to simplify the stretching in the free
335	!-- convection limit for the remaining 10% of steps.
336	zeta_tmp(0) = - zeta_max
337	num_steps_n = ( num_steps * 9 / 10 ) - 1
338	zeta_step = (zeta_max - zeta_stretch) / REAL(num_steps_n)
339
340	DO l = 1, num_steps_n
341	zeta_tmp(l) = zeta_tmp(l-1) + zeta_step
342	ENDDO
343
344	!
345	!-- Calculate stretching factor for the free convection range
346	DO WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
347	regr_old = regr
348	regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr ) &
349	)**( 10.0_wp / REAL(num_steps) )
350	ENDDO
351
352	!
353	!-- Calculate z/L range from zeta_min to zeta_stretch
354	DO l = num_steps_n+1, num_steps-1
355	zeta_tmp(l) = zeta_tmp(l-1) + zeta_step
356	zeta_step = zeta_step * regr
357	ENDDO
358
359	!
360	!-- Invert array and switch sign, then calculate Obukhov length and bulk
361	!-- Richardson number
362	z0h_min = MINVAL(z0h)
363	z0_min = MINVAL(z0)
364
365	!
366	!-- Calculate lookup table for the Richardson number versus Obukhov length
367	!-- The Richardson number (rib) is defined depending on the choice of
368	!-- boundary conditions for temperature
369	IF ( ibc_pt_b == 1 ) THEN
370	DO l = 0, num_steps-1
371	ol_tab(l) = - z_mo / zeta_tmp(num_steps-1-l)
372	rib_tab(l) = z_mo / ol_tab(l) / ( LOG( z_mo / z0_min ) &
373	- psi_m( z_mo / ol_tab(l) ) &
374	+ psi_m( z0_min / ol_tab(l) ) &
375	)**3
376	ENDDO
377	ELSE
378	DO l = 0, num_steps-1
379	ol_tab(l) = - z_mo / zeta_tmp(num_steps-1-l)
380	rib_tab(l) = z_mo / ol_tab(l) * ( LOG( z_mo / z0h_min ) &
381	- psi_h( z_mo / ol_tab(l) ) &
382	+ psi_h( z0h_min / ol_tab(l) ) &
383	) &
384	/ ( LOG( z_mo / z0_min ) &
385	- psi_m( z_mo / ol_tab(l) ) &
386	+ psi_m( z0_min / ol_tab(l) ) &
387	)**2
388	ENDDO
389	ENDIF
390
391	!
392	!-- Determine minimum values of rib in the lookup table. Set upper limit
393	!-- to critical Richardson number (0.25)
394	rib_min = MINVAL(rib_tab)
395	rib_max = 0.25 !MAXVAL(rib_tab)
396
397	DEALLOCATE( zeta_tmp )
398	ENDIF
399
400	END SUBROUTINE init_surface_layer_fluxes
401
402
403	!------------------------------------------------------------------------------!
404	! Description:
405	! ------------
406	!> Compute the absolute value of the horizontal velocity (relative to the
407	!> surface). This is required by all methods
408	!------------------------------------------------------------------------------!
409	SUBROUTINE calc_uv_total
410
411	IMPLICIT NONE
412
413
414	!$OMP PARALLEL DO PRIVATE( k )
415	!$acc kernels loop present( nzb_s_inner, u, uv_total, v ) private( j, k )
416	DO i = nxl, nxr
417	DO j = nys, nyn
418
419	k = nzb_s_inner(j,i)
420	uv_total(j,i) = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) &
421	- u(k,j,i) - u(k,j,i+1) ) )**2 + &
422	( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) &
423	- v(k,j,i) - v(k,j+1,i) ) )**2 )
424
425	!
426	!-- For too small values of the local wind, MOST does not work. A
427	!-- threshold value is thus set if required
428	! uv_total(j,i) = MAX(0.01_wp,uv_total(j,i))
429
430	ENDDO
431	ENDDO
432
433	!
434	!-- Values of uv_total need to be exchanged at the ghost boundaries
435	!$acc update host( uv_total )
436	CALL exchange_horiz_2d( uv_total )
437	!$acc update device( uv_total )
438
439	END SUBROUTINE calc_uv_total
440
441
442	!------------------------------------------------------------------------------!
443	! Description:
444	! ------------
445	!> Calculate the Obukhov length (L) and Richardson flux number (z/L)
446	!------------------------------------------------------------------------------!
447	SUBROUTINE calc_ol
448
449	IMPLICIT NONE
450
451	INTEGER(iwp) :: iter, & !< Newton iteration step
452	l !< look index
453
454	REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) :: rib !< Bulk Richardson number
455
456	REAL(wp) :: f, & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
457	f_d_ol, & !< Derivative of f
458	ol_l, & !< Lower bound of L for Newton iteration
459	ol_m, & !< Previous value of L for Newton iteration
460	ol_old, & !< Previous time step value of L
461	ol_u !< Upper bound of L for Newton iteration
462
463	IF ( TRIM( most_method ) /= 'circular' ) THEN
464
465	!$acc data present( nzb_s_inner, pt, q, qsws, rib, shf, uv_total, vpt, zu, zw )
466
467	!$OMP PARALLEL DO PRIVATE( k, z_mo )
468	!$acc kernels loop private( j, k, z_mo )
469	DO i = nxl, nxr
470	DO j = nys, nyn
471
472	k = nzb_s_inner(j,i)
473	z_mo = zu(k+1) - zw(k)
474
475	!
476	!-- Evaluate bulk Richardson number (calculation depends on
477	!-- definition based on setting of boundary conditions
478	IF ( ibc_pt_b /= 1 ) THEN
479	IF ( humidity ) THEN
480	rib(j,i) = g * z_mo * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
481	/ ( uv_total(j,i)*2 vpt(k+1,j,i) + 1.0E-20_wp )
482	ELSE
483	rib(j,i) = g * z_mo * ( pt(k+1,j,i) - pt(k,j,i) ) &
484	/ ( uv_total(j,i)*2 pt(k+1,j,i) + 1.0E-20_wp )
485	ENDIF
486	ELSE
487	!
488	!-- When using Neumann boundary conditions, the buoyancy flux
489	!-- is required but cannot be calculated at the surface, as pt
490	!-- and q are not known at the surface. Hence the values at
491	!-- first grid level are used to estimate the buoyancy flux
492	IF ( humidity ) THEN
493	rib(j,i) = - g * z_mo * ( ( 1.0_wp + 0.61_wp &
494	* q(k+1,j,i) ) * shf(j,i) + 0.61_wp &
495	* pt(k+1,j,i) * qsws(j,i) ) &
496	/ ( uv_total(j,i)*3 vpt(k+1,j,i) * kappa**2&
497	+ 1.0E-20_wp)
498	ELSE
499	rib(j,i) = - g * z_mo * shf(j,i) &
500	/ ( uv_total(j,i)*3 pt(k+1,j,i) * kappa**2 &
501	+ 1.0E-20_wp )
502	ENDIF
503	ENDIF
504
505	ENDDO
506	ENDDO
507	!$acc end data
508
509	ENDIF
510
511	!
512	!-- Calculate the Obukhov length either using a Newton iteration
513	!-- method, via a lookup table, or using the old circular way
514	IF ( TRIM( most_method ) == 'newton' ) THEN
515
516	!$OMP PARALLEL DO PRIVATE( k, z_mo )
517	!# WARNING: does not work on GPU so far because of DO-loop with
518	!# undetermined iterations
519	!!!!!!$acc kernels loop
520	DO i = nxl, nxr
521	DO j = nys, nyn
522
523	k = nzb_s_inner(j,i)
524	z_mo = zu(k+1) - zw(k)
525
526	!
527	!-- Store current value in case the Newton iteration fails
528	ol_old = ol(j,i)
529
530	!
531	!-- Ensure that the bulk Richardson number and the Obukhov
532	!-- lengtH have the same sign
533	IF ( rib(j,i) * ol(j,i) < 0.0_wp .OR. &
534	ABS( ol(j,i) ) == ol_max ) THEN
535	IF ( rib(j,i) > 0.0_wp ) ol(j,i) = 0.01_wp
536	IF ( rib(j,i) < 0.0_wp ) ol(j,i) = -0.01_wp
537	ENDIF
538	!
539	!-- Iteration to find Obukhov length
540	iter = 0
541	DO
542	iter = iter + 1
543	!
544	!-- In case of divergence, use the value of the previous time step
545	IF ( iter > 1000 ) THEN
546	ol(j,i) = ol_old
547	EXIT
548	ENDIF
549
550	ol_m = ol(j,i)
551	ol_l = ol_m - 0.001_wp * ol_m
552	ol_u = ol_m + 0.001_wp * ol_m
553
554
555	IF ( ibc_pt_b /= 1 ) THEN
556	!
557	!-- Calculate f = Ri - [...]/[...]^2 = 0
558	f = rib(j,i) - ( z_mo / ol_m ) * ( LOG( z_mo / z0h(j,i) )&
559	- psi_h( z_mo / ol_m ) &
560	+ psi_h( z0h(j,i) / ol_m ) &
561	) &
562	/ ( LOG( z_mo / z0(j,i) ) &
563	- psi_m( z_mo / ol_m ) &
564	+ psi_m( z0(j,i) / ol_m ) &
565	)**2
566
567	!
568	!-- Calculate df/dL
569	f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / z0h(j,i) ) &
570	- psi_h( z_mo / ol_u ) &
571	+ psi_h( z0h(j,i) / ol_u ) &
572	) &
573	/ ( LOG( z_mo / z0(j,i) ) &
574	- psi_m( z_mo / ol_u ) &
575	+ psi_m( z0(j,i) / ol_u ) &
576	)**2 &
577	+ ( z_mo / ol_l ) * ( LOG( z_mo / z0h(j,i) ) &
578	- psi_h( z_mo / ol_l ) &
579	+ psi_h( z0h(j,i) / ol_l ) &
580	) &
581	/ ( LOG( z_mo / z0(j,i) ) &
582	- psi_m( z_mo / ol_l ) &
583	+ psi_m( z0(j,i) / ol_l ) &
584	)**2 &
585	) / ( ol_u - ol_l )
586	ELSE
587	!
588	!-- Calculate f = Ri - 1 /[...]^3 = 0
589	f = rib(j,i) - ( z_mo / ol_m ) / ( LOG( z_mo / z0(j,i) )&
590	- psi_m( z_mo / ol_m ) &
591	+ psi_m( z0(j,i) / ol_m ) &
592	)**3
593
594	!
595	!-- Calculate df/dL
596	f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / z0(j,i) ) &
597	- psi_m( z_mo / ol_u ) &
598	+ psi_m( z0(j,i) / ol_u ) &
599	)**3 &
600	+ ( z_mo / ol_l ) / ( LOG( z_mo / z0(j,i) ) &
601	- psi_m( z_mo / ol_l ) &
602	+ psi_m( z0(j,i) / ol_l ) &
603	)**3 &
604	) / ( ol_u - ol_l )
605	ENDIF
606	!
607	!-- Calculate new L
608	ol(j,i) = ol_m - f / f_d_ol
609
610	!
611	!-- Ensure that the bulk Richardson number and the Obukhov
612	!-- length have the same sign and ensure convergence.
613	IF ( ol(j,i) * ol_m < 0.0_wp ) ol(j,i) = ol_m * 0.5_wp
614
615	!
616	!-- If unrealistic value occurs, set L to the maximum
617	!-- value that is allowed
618	IF ( ABS( ol(j,i) ) > ol_max ) THEN
619	ol(j,i) = ol_max
620	EXIT
621	ENDIF
622	!
623	!-- Check for convergence
624	IF ( ABS( ( ol(j,i) - ol_m ) / ol(j,i) ) < 1.0E-4_wp ) THEN
625	EXIT
626	ELSE
627	CYCLE
628	ENDIF
629
630	ENDDO
631
632	ENDDO
633	ENDDO
634
635	ELSEIF ( TRIM( most_method ) == 'lookup' ) THEN
636
637	!$OMP PARALLEL DO PRIVATE( k, z_mo )
638	!# WARNING: does not work on GPU so far because of DO WHILE construct
639	!!!!!!$acc kernels loop
640	DO i = nxl, nxr
641	DO j = nys, nyn
642
643	k = nzb_s_inner(j,i)
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	e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1_wp )**2
1058	!
1059	!-- As a test: cm = 0.4
1060	! e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4_wp )**2
1061	e(nzb_s_inner(j,i),j,i) = e(nzb_s_inner(j,i)+1,j,i)
1062	ENDDO
1063	ENDDO
1064	ENDIF
1065
1066	END SUBROUTINE calc_surface_fluxes
1067
1068
1069	!
1070	!-- Integrated stability function for momentum
1071	FUNCTION psi_m( zeta )
1072
1073	USE kinds
1074
1075	IMPLICIT NONE
1076
1077	REAL(wp) :: psi_m !< Integrated similarity function result
1078	REAL(wp) :: zeta !< Stability parameter z/L
1079	REAL(wp) :: x !< dummy variable
1080
1081	REAL(wp), PARAMETER :: a = 1.0_wp !< constant
1082	REAL(wp), PARAMETER :: b = 0.66666666666_wp !< constant
1083	REAL(wp), PARAMETER :: c = 5.0_wp !< constant
1084	REAL(wp), PARAMETER :: d = 0.35_wp !< constant
1085	REAL(wp), PARAMETER :: c_d_d = c / d !< constant
1086	REAL(wp), PARAMETER :: bc_d_d = b * c / d !< constant
1087
1088
1089	IF ( zeta < 0.0_wp ) THEN
1090	x = SQRT( SQRT(1.0_wp - 16.0_wp * zeta ) )
1091	psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2 &
1092	* ( 1.0_wp + x*2 ) 0.125_wp )
1093	ELSE
1094
1095	psi_m = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - a * zeta &
1096	- bc_d_d
1097	!
1098	!-- Old version for stable conditions (only valid for z/L < 0.5)
1099	!-- psi_m = - 5.0_wp * zeta
1100
1101	ENDIF
1102
1103	END FUNCTION psi_m
1104
1105
1106	!
1107	!-- Integrated stability function for heat and moisture
1108	FUNCTION psi_h( zeta )
1109
1110	USE kinds
1111
1112	IMPLICIT NONE
1113
1114	REAL(wp) :: psi_h !< Integrated similarity function result
1115	REAL(wp) :: zeta !< Stability parameter z/L
1116	REAL(wp) :: x !< dummy variable
1117
1118	REAL(wp), PARAMETER :: a = 1.0_wp !< constant
1119	REAL(wp), PARAMETER :: b = 0.66666666666_wp !< constant
1120	REAL(wp), PARAMETER :: c = 5.0_wp !< constant
1121	REAL(wp), PARAMETER :: d = 0.35_wp !< constant
1122	REAL(wp), PARAMETER :: c_d_d = c / d !< constant
1123	REAL(wp), PARAMETER :: bc_d_d = b * c / d !< constant
1124
1125
1126	IF ( zeta < 0.0_wp ) THEN
1127	x = SQRT(1.0_wp - 16.0_wp * zeta )
1128	psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
1129	ELSE
1130	psi_h = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - (1.0_wp &
1131	+ 0.66666666666_wp * a * zeta )**1.5_wp - bc_d_d &
1132	+ 1.0_wp
1133	!
1134	!-- Old version for stable conditions (only valid for z/L < 0.5)
1135	!-- psi_h = - 5.0_wp * zeta
1136	ENDIF
1137
1138	END FUNCTION psi_h
1139
1140	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 |