Home

Context Navigation

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

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