Home

Context Navigation

source: palm/trunk/SOURCE/ocean_mod.f90 @ 4109

Last change on this file since 4109 was 4109, checked in by suehring, 5 years ago
Control discretization of advection term: separate initialization of WS advection flags for momentum and scalars. In this context, resort the bits and do some minor formatting; Make initialization of scalar-advection flags more flexible, i.e. introduce an arguemnt list to indicate non-cyclic boundaries (required for decycled scalars such as chemical species or aerosols); Introduce extended 'degradation zones', where horizontal advection of passive scalars is discretized by first-order scheme at all grid points that in the vicinity of buildings (<= 3 grid points). Even though no building is within the numerical stencil, first-order scheme is used. At fourth and fifth grid point the order of the horizontal advection scheme is successively upgraded. These extended degradation zones are used to avoid stationary numerical oscillations, which are responsible for high concentration maxima that may appear under shear-free stable conditions. Therefore, an additional 3D interger array used to store control flags is introduced; Change interface for scalar advection routine; Bugfix, avoid uninitialized value sk_num in vector version of WS scalar advection; Chemistry: Decycling boundary conditions are only set at the ghost points not on the prognostic grid points; Land-surface model: Relax checks for non-consistent initialization in case static or dynamic input is provided. For example, soil_temperature or deep_soil_temperature is not mandatory any more if dynamic input is available. Also, improper settings of x_type in namelist are only checked if no static file is available.
Property svn:keywords set to `Id`
File size: 86.5 KB

Line
1	!> @file ocean_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 2017-2018 Leibniz Universitaet Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	! Pass integer flag array as well as boundary flags to WS scalar advection
23	! routine
24	!
25	! Former revisions:
26	! -----------------
27	! $Id: ocean_mod.f90 4109 2019-07-22 17:00:34Z suehring $
28	! implemented ocean_actions
29	!
30	! 3767 2019-02-27 08:18:02Z raasch
31	! unused variable for file index and tmp_2d removed from rrd-subroutine parameter
32	! list
33	!
34	! 3719 2019-02-06 13:10:18Z kanani
35	! Changed log_point to log_point_s, otherwise this overlaps with
36	! 'all progn.equations' cpu measurement.
37	!
38	! 3684 2019-01-20 20:20:58Z knoop
39	! nopointer option removed
40	!
41	! 3614 2018-12-10 07:05:46Z raasch
42	! unused variables removed
43	!
44	! 3568 2018-11-27 16:07:59Z raasch
45	! bugifx: calculate equation of state for seawater even if salinity is switched
46	! off
47	!
48	! 3421 2018-10-24 18:39:32Z gronemeier
49	! Renamed output variables
50	!
51	! 3381 2018-10-19 13:09:06Z raasch
52	! spin-up cooling for ocean surface implemented, see new parameter
53	! surface_cooling_spinup_time
54	!
55	! 3311 2018-10-05 12:34:56Z raasch
56	! check if ocean mode is used for invalid combinations
57	!
58	! 3303 2018-10-03 12:04:15Z raasch
59	! salinity allowed to be switched off
60	!
61	! 3302 2018-10-03 02:39:40Z raasch
62	! Craik Leibovich force (Stokes drift) + wave breaking effect added
63	!
64	! 3294 2018-10-01 02:37:10Z raasch
65	! initial revision
66	!
67	!
68	! Authors:
69	! --------
70	! @author Siegfried Raasch
71	!
72	! Description:
73	! ------------
74	!> This module contains relevant code for PALM's ocean mode, e.g. the
75	!> prognostic equation for salinity, the equation of state for seawater,
76	!> the Craik Leibovich force (Stokes force), and wave breaking effects
77	!------------------------------------------------------------------------------!
78	MODULE ocean_mod
79
80
81	USE arrays_3d, &
82	ONLY: prho, prho_1, rho_ocean, rho_1, sa, sa_init, sa_1, sa_2, sa_3, &
83	sa_p, tsa_m, flux_l_sa, flux_s_sa, diss_l_sa, diss_s_sa
84
85	USE control_parameters, &
86	ONLY: atmos_ocean_sign, bottom_salinityflux, &
87	constant_top_salinityflux, ocean_mode, top_salinityflux, &
88	wall_salinityflux, loop_optimization, ws_scheme_sca
89
90	USE kinds
91
92	USE pegrid, &
93	ONLY: threads_per_task
94
95	USE statistics, &
96	ONLY: sums_wssas_ws_l
97
98
99	IMPLICIT NONE
100
101	CHARACTER (LEN=20) :: bc_sa_t = 'neumann' !< namelist parameter
102
103	INTEGER(iwp) :: ibc_sa_t !< integer flag for bc_sa_t
104	INTEGER(iwp) :: iran_ocean = -1234567 !< random number used for wave breaking effects
105
106	INTEGER(iwp) :: sa_vertical_gradient_level_ind(10) = -9999 !< grid index values of sa_vertical_gradient_level(s)
107
108	LOGICAL :: salinity = .TRUE. !< switch for using salinity
109	LOGICAL :: stokes_force = .FALSE. !< switch to switch on the Stokes force
110	LOGICAL :: wave_breaking = .FALSE. !< switch to switch on wave breaking effects
111	LOGICAL :: surface_cooling_switched_off = .FALSE. !< variable to check if surface heat flux has been switched off
112
113	REAL(wp) :: alpha_wave_breaking = 3.0_wp !< coefficient for wave breaking generated turbulence from Noh et al. (2004), JPO
114	REAL(wp) :: prho_reference !< reference state of potential density at ocean surface
115	REAL(wp) :: sa_surface = 35.0_wp !< surface salinity, namelist parameter
116	REAL(wp) :: sa_vertical_gradient(10) = 0.0_wp !< namelist parameter
117	REAL(wp) :: sa_vertical_gradient_level(10) = -999999.9_wp !< namelist parameter
118	REAL(wp) :: stokes_waveheight = 0.0_wp !< wave height assumed for Stokes drift velocity
119	REAL(wp) :: stokes_wavelength = 0.0_wp !< wavelength assumed for Stokes drift velocity
120	REAL(wp) :: surface_cooling_spinup_time = 999999.9_wp !< time after which surface heat flux is switched off
121	REAL(wp) :: timescale_wave_breaking !< time scale of random forcing
122	REAL(wp) :: u_star_wave_breaking !< to store the absolute value of friction velocity at the ocean surface
123
124	REAL(wp), DIMENSION(12), PARAMETER :: nom = &
125	(/ 9.9984085444849347D2, 7.3471625860981584D0, &
126	-5.3211231792841769D-2, 3.6492439109814549D-4, &
127	2.5880571023991390D0, -6.7168282786692354D-3, &
128	1.9203202055760151D-3, 1.1798263740430364D-2, &
129	9.8920219266399117D-8, 4.6996642771754730D-6, &
130	-2.5862187075154352D-8, -3.2921414007960662D-12 /)
131	!< constants used in equation of state for seawater
132
133	REAL(wp), DIMENSION(13), PARAMETER :: den = &
134	(/ 1.0D0, 7.2815210113327091D-3, &
135	-4.4787265461983921D-5, 3.3851002965802430D-7, &
136	1.3651202389758572D-10, 1.7632126669040377D-3, &
137	-8.8066583251206474D-6, -1.8832689434804897D-10, &
138	5.7463776745432097D-6, 1.4716275472242334D-9, &
139	6.7103246285651894D-6, -2.4461698007024582D-17, &
140	-9.1534417604289062D-18 /)
141	!< constants used in equation of state for seawater
142
143	SAVE
144
145	PUBLIC :: bc_sa_t, ibc_sa_t, prho_reference, sa_surface, &
146	sa_vertical_gradient, sa_vertical_gradient_level, &
147	sa_vertical_gradient_level_ind, stokes_force, wave_breaking
148
149
150	INTERFACE eqn_state_seawater
151	MODULE PROCEDURE eqn_state_seawater
152	MODULE PROCEDURE eqn_state_seawater_ij
153	END INTERFACE eqn_state_seawater
154
155	INTERFACE eqn_state_seawater_func
156	MODULE PROCEDURE eqn_state_seawater_func
157	END INTERFACE eqn_state_seawater_func
158
159	INTERFACE ocean_check_parameters
160	MODULE PROCEDURE ocean_check_parameters
161	END INTERFACE ocean_check_parameters
162
163	INTERFACE ocean_check_data_output
164	MODULE PROCEDURE ocean_check_data_output
165	END INTERFACE ocean_check_data_output
166
167	INTERFACE ocean_check_data_output_pr
168	MODULE PROCEDURE ocean_check_data_output_pr
169	END INTERFACE ocean_check_data_output_pr
170
171	INTERFACE ocean_define_netcdf_grid
172	MODULE PROCEDURE ocean_define_netcdf_grid
173	END INTERFACE ocean_define_netcdf_grid
174
175	INTERFACE ocean_data_output_2d
176	MODULE PROCEDURE ocean_data_output_2d
177	END INTERFACE ocean_data_output_2d
178
179	INTERFACE ocean_data_output_3d
180	MODULE PROCEDURE ocean_data_output_3d
181	END INTERFACE ocean_data_output_3d
182
183	INTERFACE ocean_header
184	MODULE PROCEDURE ocean_header
185	END INTERFACE ocean_header
186
187	INTERFACE ocean_init
188	MODULE PROCEDURE ocean_init
189	END INTERFACE ocean_init
190
191	INTERFACE ocean_init_arrays
192	MODULE PROCEDURE ocean_init_arrays
193	END INTERFACE ocean_init_arrays
194
195	INTERFACE ocean_parin
196	MODULE PROCEDURE ocean_parin
197	END INTERFACE ocean_parin
198
199	INTERFACE ocean_actions
200	MODULE PROCEDURE ocean_actions
201	MODULE PROCEDURE ocean_actions_ij
202	END INTERFACE ocean_actions
203
204	INTERFACE ocean_prognostic_equations
205	MODULE PROCEDURE ocean_prognostic_equations
206	MODULE PROCEDURE ocean_prognostic_equations_ij
207	END INTERFACE ocean_prognostic_equations
208
209	INTERFACE ocean_swap_timelevel
210	MODULE PROCEDURE ocean_swap_timelevel
211	END INTERFACE ocean_swap_timelevel
212
213	INTERFACE ocean_rrd_global
214	MODULE PROCEDURE ocean_rrd_global
215	END INTERFACE ocean_rrd_global
216
217	INTERFACE ocean_rrd_local
218	MODULE PROCEDURE ocean_rrd_local
219	END INTERFACE ocean_rrd_local
220
221	INTERFACE ocean_wrd_global
222	MODULE PROCEDURE ocean_wrd_global
223	END INTERFACE ocean_wrd_global
224
225	INTERFACE ocean_wrd_local
226	MODULE PROCEDURE ocean_wrd_local
227	END INTERFACE ocean_wrd_local
228
229	INTERFACE ocean_3d_data_averaging
230	MODULE PROCEDURE ocean_3d_data_averaging
231	END INTERFACE ocean_3d_data_averaging
232
233	INTERFACE stokes_drift_terms
234	MODULE PROCEDURE stokes_drift_terms
235	MODULE PROCEDURE stokes_drift_terms_ij
236	END INTERFACE stokes_drift_terms
237
238	INTERFACE wave_breaking_term
239	MODULE PROCEDURE wave_breaking_term
240	MODULE PROCEDURE wave_breaking_term_ij
241	END INTERFACE wave_breaking_term
242
243	PRIVATE
244	!
245	!-- Add INTERFACES that must be available to other modules (alphabetical order)
246	PUBLIC eqn_state_seawater, ocean_actions, ocean_check_data_output, &
247	ocean_check_data_output_pr, ocean_check_parameters, &
248	ocean_data_output_2d, ocean_data_output_3d, &
249	ocean_define_netcdf_grid, ocean_header, ocean_init, &
250	ocean_init_arrays, ocean_parin, ocean_prognostic_equations, &
251	ocean_swap_timelevel, ocean_rrd_global, ocean_rrd_local, &
252	ocean_wrd_global, ocean_wrd_local, ocean_3d_data_averaging, &
253	stokes_drift_terms, wave_breaking_term
254
255
256	CONTAINS
257
258	!------------------------------------------------------------------------------!
259	! Description:
260	! ------------
261	!> Equation of state for seawater as a function of potential temperature,
262	!> salinity, and pressure.
263	!> For coefficients see Jackett et al., 2006: J. Atm. Ocean Tech.
264	!> eqn_state_seawater calculates the potential density referred at hyp(0).
265	!> eqn_state_seawater_func calculates density.
266	!> TODO: so far, routine is not adjusted to use topography
267	!------------------------------------------------------------------------------!
268	SUBROUTINE eqn_state_seawater
269
270	USE arrays_3d, &
271	ONLY: hyp, prho, pt_p, rho_ocean, sa_p
272	USE indices, &
273	ONLY: nxl, nxr, nyn, nys, nzb, nzt
274
275	USE surface_mod, &
276	ONLY : bc_h
277
278	IMPLICIT NONE
279
280	INTEGER(iwp) :: i !< running index x direction
281	INTEGER(iwp) :: j !< running index y direction
282	INTEGER(iwp) :: k !< running index z direction
283	INTEGER(iwp) :: m !< running index surface elements
284
285	REAL(wp) :: pden !< temporary scalar user for calculating density
286	REAL(wp) :: pnom !< temporary scalar user for calculating density
287	REAL(wp) :: p1 !< temporary scalar user for calculating density
288	REAL(wp) :: p2 !< temporary scalar user for calculating density
289	REAL(wp) :: p3 !< temporary scalar user for calculating density
290	REAL(wp) :: pt1 !< temporary scalar user for calculating density
291	REAL(wp) :: pt2 !< temporary scalar user for calculating density
292	REAL(wp) :: pt3 !< temporary scalar user for calculating density
293	REAL(wp) :: pt4 !< temporary scalar user for calculating density
294	REAL(wp) :: sa1 !< temporary scalar user for calculating density
295	REAL(wp) :: sa15 !< temporary scalar user for calculating density
296	REAL(wp) :: sa2 !< temporary scalar user for calculating density
297
298
299	DO i = nxl, nxr
300	DO j = nys, nyn
301	DO k = nzb+1, nzt
302	!
303	!-- Pressure is needed in dbar
304	p1 = hyp(k) * 1E-4_wp
305	p2 = p1 * p1
306	p3 = p2 * p1
307
308	!
309	!-- Temperature needed in degree Celsius
310	pt1 = pt_p(k,j,i) - 273.15_wp
311	pt2 = pt1 * pt1
312	pt3 = pt1 * pt2
313	pt4 = pt2 * pt2
314
315	sa1 = sa_p(k,j,i)
316	sa15 = sa1 * SQRT( sa1 )
317	sa2 = sa1 * sa1
318
319	pnom = nom(1) + nom(2)pt1 + nom(3)pt2 + &
320	nom(4)pt3 + nom(5)sa1 + nom(6)sa1pt1 + &
321	nom(7)*sa2
322
323	pden = den(1) + den(2)pt1 + den(3)pt2 + &
324	den(4)pt3 + den(5)pt4 + den(6)*sa1 + &
325	den(7)sa1pt1 + den(8)sa1pt3 + den(9)*sa15 + &
326	den(10)sa15pt2
327	!
328	!-- Potential density (without pressure terms)
329	prho(k,j,i) = pnom / pden
330
331	pnom = pnom + nom(8)p1 + nom(9)p1*pt2 + &
332	nom(10)p1sa1 + nom(11)p2 + nom(12)p2*pt2
333
334	pden = pden + den(11)p1 + den(12)p2*pt3 + &
335	den(13)p3pt1
336
337	!
338	!-- In-situ density
339	rho_ocean(k,j,i) = pnom / pden
340
341	ENDDO
342
343	!
344	!-- Neumann conditions are assumed at top boundary and currently also at
345	!-- bottom boundary (see comment lines below)
346	prho(nzt+1,j,i) = prho(nzt,j,i)
347	rho_ocean(nzt+1,j,i) = rho_ocean(nzt,j,i)
348
349	ENDDO
350	ENDDO
351	!
352	!-- Neumann conditions at up/downward-facing surfaces
353	!$OMP PARALLEL DO PRIVATE( i, j, k )
354	DO m = 1, bc_h(0)%ns
355	i = bc_h(0)%i(m)
356	j = bc_h(0)%j(m)
357	k = bc_h(0)%k(m)
358	prho(k-1,j,i) = prho(k,j,i)
359	rho_ocean(k-1,j,i) = rho_ocean(k,j,i)
360	ENDDO
361	!
362	!-- Downward facing surfaces
363	!$OMP PARALLEL DO PRIVATE( i, j, k )
364	DO m = 1, bc_h(1)%ns
365	i = bc_h(1)%i(m)
366	j = bc_h(1)%j(m)
367	k = bc_h(1)%k(m)
368	prho(k+1,j,i) = prho(k,j,i)
369	rho_ocean(k+1,j,i) = rho_ocean(k,j,i)
370	ENDDO
371
372	END SUBROUTINE eqn_state_seawater
373
374
375	!------------------------------------------------------------------------------!
376	! Description:
377	! ------------
378	!> Same as above, but for grid point i,j
379	!------------------------------------------------------------------------------!
380	SUBROUTINE eqn_state_seawater_ij( i, j )
381
382	USE arrays_3d, &
383	ONLY: hyp, prho, pt_p, rho_ocean, sa_p
384
385	USE indices, &
386	ONLY: nzb, nzt
387
388	USE surface_mod, &
389	ONLY : bc_h
390
391	IMPLICIT NONE
392
393	INTEGER(iwp) :: i !< running index x direction
394	INTEGER(iwp) :: j !< running index y direction
395	INTEGER(iwp) :: k !< running index z direction
396	INTEGER(iwp) :: m !< running index surface elements
397	INTEGER(iwp) :: surf_e !< end index of surface elements at (j,i)-gridpoint
398	INTEGER(iwp) :: surf_s !< start index of surface elements at (j,i)-gridpoint
399
400	REAL(wp) :: pden !< temporary scalar user for calculating density
401	REAL(wp) :: pnom !< temporary scalar user for calculating density
402	REAL(wp) :: p1 !< temporary scalar user for calculating density
403	REAL(wp) :: p2 !< temporary scalar user for calculating density
404	REAL(wp) :: p3 !< temporary scalar user for calculating density
405	REAL(wp) :: pt1 !< temporary scalar user for calculating density
406	REAL(wp) :: pt2 !< temporary scalar user for calculating density
407	REAL(wp) :: pt3 !< temporary scalar user for calculating density
408	REAL(wp) :: pt4 !< temporary scalar user for calculating density
409	REAL(wp) :: sa1 !< temporary scalar user for calculating density
410	REAL(wp) :: sa15 !< temporary scalar user for calculating density
411	REAL(wp) :: sa2 !< temporary scalar user for calculating density
412
413	DO k = nzb+1, nzt
414	!
415	!-- Pressure is needed in dbar
416	p1 = hyp(k) * 1E-4_wp
417	p2 = p1 * p1
418	p3 = p2 * p1
419	!
420	!-- Temperature needed in degree Celsius
421	pt1 = pt_p(k,j,i) - 273.15_wp
422	pt2 = pt1 * pt1
423	pt3 = pt1 * pt2
424	pt4 = pt2 * pt2
425
426	sa1 = sa_p(k,j,i)
427	sa15 = sa1 * SQRT( sa1 )
428	sa2 = sa1 * sa1
429
430	pnom = nom(1) + nom(2)pt1 + nom(3)pt2 + &
431	nom(4)pt3 + nom(5)sa1 + nom(6)sa1pt1 + nom(7)*sa2
432
433	pden = den(1) + den(2)pt1 + den(3)pt2 + &
434	den(4)pt3 + den(5)pt4 + den(6)*sa1 + &
435	den(7)sa1pt1 + den(8)sa1pt3 + den(9)*sa15 + &
436	den(10)sa15pt2
437	!
438	!-- Potential density (without pressure terms)
439	prho(k,j,i) = pnom / pden
440
441	pnom = pnom + nom(8)p1 + nom(9)p1*pt2 + &
442	nom(10)p1sa1 + nom(11)p2 + nom(12)p2*pt2
443	pden = pden + den(11)p1 + den(12)p2*pt3 + &
444	den(13)p3pt1
445
446	!
447	!-- In-situ density
448	rho_ocean(k,j,i) = pnom / pden
449
450	ENDDO
451	!
452	!-- Neumann conditions at up/downward-facing walls
453	surf_s = bc_h(0)%start_index(j,i)
454	surf_e = bc_h(0)%end_index(j,i)
455	DO m = surf_s, surf_e
456	k = bc_h(0)%k(m)
457	prho(k-1,j,i) = prho(k,j,i)
458	rho_ocean(k-1,j,i) = rho_ocean(k,j,i)
459	ENDDO
460	!
461	!-- Downward facing surfaces
462	surf_s = bc_h(1)%start_index(j,i)
463	surf_e = bc_h(1)%end_index(j,i)
464	DO m = surf_s, surf_e
465	k = bc_h(1)%k(m)
466	prho(k+1,j,i) = prho(k,j,i)
467	rho_ocean(k+1,j,i) = rho_ocean(k,j,i)
468	ENDDO
469	!
470	!-- Neumann condition are assumed at top boundary
471	prho(nzt+1,j,i) = prho(nzt,j,i)
472	rho_ocean(nzt+1,j,i) = rho_ocean(nzt,j,i)
473
474	END SUBROUTINE eqn_state_seawater_ij
475
476
477	!------------------------------------------------------------------------------!
478	! Description:
479	! ------------
480	!> Equation of state for seawater as a function
481	!------------------------------------------------------------------------------!
482	REAL(wp) FUNCTION eqn_state_seawater_func( p, pt, sa )
483
484	IMPLICIT NONE
485
486	REAL(wp) :: p !< temporary scalar user for calculating density
487	REAL(wp) :: p1 !< temporary scalar user for calculating density
488	REAL(wp) :: p2 !< temporary scalar user for calculating density
489	REAL(wp) :: p3 !< temporary scalar user for calculating density
490	REAL(wp) :: pt !< temporary scalar user for calculating density
491	REAL(wp) :: pt1 !< temporary scalar user for calculating density
492	REAL(wp) :: pt2 !< temporary scalar user for calculating density
493	REAL(wp) :: pt3 !< temporary scalar user for calculating density
494	REAL(wp) :: pt4 !< temporary scalar user for calculating density
495	REAL(wp) :: sa !< temporary scalar user for calculating density
496	REAL(wp) :: sa15 !< temporary scalar user for calculating density
497	REAL(wp) :: sa2 !< temporary scalar user for calculating density
498
499	!
500	!-- Pressure is needed in dbar
501	p1 = p * 1.0E-4_wp
502	p2 = p1 * p1
503	p3 = p2 * p1
504
505	!
506	!-- Temperature needed in degree Celsius
507	pt1 = pt - 273.15_wp
508	pt2 = pt1 * pt1
509	pt3 = pt1 * pt2
510	pt4 = pt2 * pt2
511
512	sa15 = sa * SQRT( sa )
513	sa2 = sa * sa
514
515
516	eqn_state_seawater_func = &
517	( nom(1) + nom(2)pt1 + nom(3)pt2 + nom(4)*pt3 + &
518	nom(5)sa + nom(6)sapt1 + nom(7)sa2 + nom(8)*p1 + &
519	nom(9)p1pt2 + nom(10)p1sa + nom(11)p2 + nom(12)p2*pt2 &
520	) / &
521	( den(1) + den(2)pt1 + den(3)pt2 + den(4)*pt3 + &
522	den(5)pt4 + den(6)sa + den(7)sapt1 + den(8)sapt3 + &
523	den(9)sa15 + den(10)sa15pt2 + den(11)p1 + den(12)p2pt3 + &
524	den(13)p3pt1 &
525	)
526
527
528	END FUNCTION eqn_state_seawater_func
529
530
531	!------------------------------------------------------------------------------!
532	! Description:
533	! ------------
534	!> Reads the ocean parameters namelist
535	!------------------------------------------------------------------------------!
536	SUBROUTINE ocean_parin
537
538	IMPLICIT NONE
539
540	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
541
542
543	NAMELIST /ocean_parameters/ bc_sa_t, bottom_salinityflux, salinity, &
544	sa_surface, sa_vertical_gradient, sa_vertical_gradient_level, &
545	stokes_waveheight, stokes_wavelength, surface_cooling_spinup_time,&
546	top_salinityflux, wall_salinityflux, wave_breaking
547
548	!
549	!-- Try to find the namelist
550	REWIND ( 11 )
551	line = ' '
552	DO WHILE ( INDEX( line, '&ocean_parameters' ) == 0 )
553	READ ( 11, '(A)', END=12 ) line
554	ENDDO
555	BACKSPACE ( 11 )
556
557	!
558	!-- Read namelist
559	READ ( 11, ocean_parameters, ERR = 10 )
560	!
561	!-- Set switch that enables PALM's ocean mode
562	ocean_mode = .TRUE.
563
564	GOTO 12
565
566	10 BACKSPACE( 11 )
567	READ( 11 , '(A)') line
568	CALL parin_fail_message( 'ocean_parameters', line )
569
570	12 CONTINUE
571
572	END SUBROUTINE ocean_parin
573
574	!------------------------------------------------------------------------------!
575	! Description:
576	! ------------
577	!> Check parameters routine for the ocean mode
578	!------------------------------------------------------------------------------!
579	SUBROUTINE ocean_check_parameters
580
581	USE control_parameters, &
582	ONLY: coupling_char, coupling_mode, initializing_actions, &
583	message_string, use_top_fluxes
584
585	USE pmc_interface, &
586	ONLY: nested_run
587
588	IMPLICIT NONE
589
590
591	!
592	!-- Check for invalid combinations
593	IF ( nested_run ) THEN
594	message_string = 'ocean mode not allowed for nesting'
595	CALL message( 'ocean_check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
596	ENDIF
597
598	IF ( TRIM( initializing_actions ) == 'cyclic_fill' ) THEN
599	message_string = 'ocean mode does not allow cyclic-fill initialization'
600	CALL message( 'ocean_check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
601	ENDIF
602
603	!
604	!-- Check ocean setting
605	IF ( TRIM( coupling_mode ) == 'uncoupled' .AND. &
606	TRIM( coupling_char ) == '_O' .AND. &
607	.NOT. ocean_mode ) THEN
608
609	!
610	!-- Check whether an (uncoupled) atmospheric run has been declared as an
611	!-- ocean run (this setting is done via palmrun-option -y)
612	message_string = 'ocean mode does not allow coupling_char = "' // &
613	TRIM( coupling_char ) // '" set by palmrun-option "-y"'
614	CALL message( 'ocean_check_parameters', 'PA0317', 1, 2, 0, 6, 0 )
615
616	ENDIF
617
618	!
619	!-- Ocean version must use flux boundary conditions at the top
620	IF ( .NOT. use_top_fluxes ) THEN
621	message_string = 'use_top_fluxes must be .TRUE. in ocean mode'
622	CALL message( 'ocean_check_parameters', 'PA0042', 1, 2, 0, 6, 0 )
623	ENDIF
624
625	!
626	!-- Boundary conditions for salinity
627	IF ( bc_sa_t == 'dirichlet' ) THEN
628	ibc_sa_t = 0
629	ELSEIF ( bc_sa_t == 'neumann' ) THEN
630	ibc_sa_t = 1
631	ELSE
632	message_string = 'unknown boundary condition: bc_sa_t = "' // &
633	TRIM( bc_sa_t ) // '"'
634	CALL message( 'ocean_check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
635	ENDIF
636
637	IF ( top_salinityflux == 9999999.9_wp ) constant_top_salinityflux = .FALSE.
638
639	IF ( .NOT. salinity ) THEN
640	IF ( ( bottom_salinityflux /= 0.0_wp .AND. &
641	bottom_salinityflux /= 9999999.9_wp ) .OR. &
642	( top_salinityflux /= 0.0_wp .AND. &
643	top_salinityflux /= 9999999.9_wp ) ) &
644	THEN
645	message_string = 'salinityflux must not be set for ocean run ' // &
646	'without salinity'
647	CALL message( 'ocean_check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
648	ENDIF
649	ENDIF
650
651	IF ( ibc_sa_t == 1 .AND. top_salinityflux == 9999999.9_wp ) THEN
652	message_string = 'boundary condition: bc_sa_t = "' // &
653	TRIM( bc_sa_t ) // '" requires to set top_salinityflux'
654	CALL message( 'ocean_check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
655	ENDIF
656
657	!
658	!-- A fixed salinity at the top implies Dirichlet boundary condition for
659	!-- salinity. In this case specification of a constant salinity flux is
660	!-- forbidden.
661	IF ( ibc_sa_t == 0 .AND. constant_top_salinityflux .AND. &
662	top_salinityflux /= 0.0_wp ) THEN
663	message_string = 'boundary condition: bc_sa_t = "' // &
664	TRIM( bc_sa_t ) // '" is not allowed with ' // &
665	'top_salinityflux /= 0.0'
666	CALL message( 'ocean_check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
667	ENDIF
668
669	!
670	!-- Check if Stokes force is to be used
671	IF ( stokes_waveheight /= 0.0_wp .AND. stokes_wavelength /= 0.0_wp ) THEN
672	stokes_force = .TRUE.
673	ELSE
674	IF ( ( stokes_waveheight <= 0.0_wp .AND. stokes_wavelength > 0.0_wp ) &
675	.OR. &
676	( stokes_waveheight > 0.0_wp .AND. stokes_wavelength <= 0.0_wp ) &
677	.OR. &
678	( stokes_waveheight < 0.0_wp .AND. stokes_wavelength < 0.0_wp ) ) &
679	THEN
680	message_string = 'wrong settings for stokes_wavelength and/or ' // &
681	'stokes_waveheight'
682	CALL message( 'ocean_check_parameters', 'PA0460', 1, 2, 0, 6, 0 )
683	ENDIF
684	ENDIF
685
686	END SUBROUTINE ocean_check_parameters
687
688
689	!------------------------------------------------------------------------------!
690	! Description:
691	! ------------
692	!> Check data output.
693	!------------------------------------------------------------------------------!
694	SUBROUTINE ocean_check_data_output( var, unit )
695
696	IMPLICIT NONE
697
698	CHARACTER (LEN=*) :: unit !< unit of output variable
699	CHARACTER (LEN=*) :: var !< name of output variable
700
701
702	SELECT CASE ( TRIM( var ) )
703
704	CASE ( 'rho_sea_water' )
705	unit = 'kg/m3'
706
707	CASE ( 'sa' )
708	unit = 'psu'
709
710	CASE DEFAULT
711	unit = 'illegal'
712
713	END SELECT
714
715	END SUBROUTINE ocean_check_data_output
716
717
718	!------------------------------------------------------------------------------!
719	! Description:
720	! ------------
721	!> Check data output of profiles
722	!------------------------------------------------------------------------------!
723	SUBROUTINE ocean_check_data_output_pr( variable, var_count, unit, dopr_unit )
724
725	USE arrays_3d, &
726	ONLY: zu, zw
727
728	USE control_parameters, &
729	ONLY: data_output_pr
730
731	USE indices
732
733	USE profil_parameter
734
735	USE statistics
736
737	IMPLICIT NONE
738
739	CHARACTER (LEN=*) :: unit !<
740	CHARACTER (LEN=*) :: variable !<
741	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
742
743	INTEGER(iwp) :: var_count !<
744
745	SELECT CASE ( TRIM( variable ) )
746
747	CASE ( 'prho' )
748	dopr_index(var_count) = 71
749	dopr_unit = 'kg/m3'
750	hom(:,2,71,:) = SPREAD( zu, 2, statistic_regions+1 )
751	unit = dopr_unit
752
753	CASE ( 'rho_sea_water' )
754	dopr_index(var_count) = 64
755	dopr_unit = 'kg/m3'
756	hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
757	IF ( data_output_pr(var_count)(1:1) == '#' ) THEN
758	dopr_initial_index(var_count) = 77
759	hom(:,2,77,:) = SPREAD( zu, 2, statistic_regions+1 )
760	hom(nzb,2,77,:) = 0.0_wp ! because zu(nzb) is negative
761	data_output_pr(var_count) = data_output_pr(var_count)(2:)
762	ENDIF
763	unit = dopr_unit
764
765	CASE ( 'sa', '#sa' )
766	dopr_index(var_count) = 23
767	dopr_unit = 'psu'
768	hom(:,2,23,:) = SPREAD( zu, 2, statistic_regions+1 )
769	IF ( data_output_pr(var_count)(1:1) == '#' ) THEN
770	dopr_initial_index(var_count) = 26
771	hom(:,2,26,:) = SPREAD( zu, 2, statistic_regions+1 )
772	hom(nzb,2,26,:) = 0.0_wp ! because zu(nzb) is negative
773	data_output_pr(var_count) = data_output_pr(var_count)(2:)
774	ENDIF
775	unit = dopr_unit
776
777	CASE ( 'w"sa"' )
778	dopr_index(var_count) = 65
779	dopr_unit = 'psu m/s'
780	hom(:,2,65,:) = SPREAD( zw, 2, statistic_regions+1 )
781	unit = dopr_unit
782
783	CASE ( 'wsa' )
784	dopr_index(var_count) = 66
785	dopr_unit = 'psu m/s'
786	hom(:,2,66,:) = SPREAD( zw, 2, statistic_regions+1 )
787	unit = dopr_unit
788
789	CASE ( 'wsa' )
790	dopr_index(var_count) = 67
791	dopr_unit = 'psu m/s'
792	hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
793	unit = dopr_unit
794
795	CASE DEFAULT
796	unit = 'illegal'
797
798	END SELECT
799
800
801	END SUBROUTINE ocean_check_data_output_pr
802
803
804	!------------------------------------------------------------------------------!
805	! Description:
806	! ------------
807	!> Define appropriate grid for netcdf variables.
808	!> It is called out from subroutine netcdf.
809	!------------------------------------------------------------------------------!
810	SUBROUTINE ocean_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
811
812	IMPLICIT NONE
813
814	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< x grid of output variable
815	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< y grid of output variable
816	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< z grid of output variable
817	CHARACTER (LEN=*), INTENT(IN) :: var !< name of output variable
818
819	LOGICAL, INTENT(OUT) :: found !< flag if output variable is found
820
821	found = .TRUE.
822
823	!
824	!-- Check for the grid
825	SELECT CASE ( TRIM( var ) )
826
827	CASE ( 'rho_sea_water', 'rho_sea_water_xy', 'rho_sea_water_xz', &
828	'rho_sea_water_yz', 'sa', 'sa_xy', 'sa_xz', 'sa_yz' )
829	grid_x = 'x'
830	grid_y = 'y'
831	grid_z = 'zu'
832
833	CASE DEFAULT
834	found = .FALSE.
835	grid_x = 'none'
836	grid_y = 'none'
837	grid_z = 'none'
838
839	END SELECT
840
841	END SUBROUTINE ocean_define_netcdf_grid
842
843
844	!------------------------------------------------------------------------------!
845	! Description:
846	! ------------
847	!> Average 3D data.
848	!------------------------------------------------------------------------------!
849	SUBROUTINE ocean_3d_data_averaging( mode, variable )
850
851
852	USE arrays_3d, &
853	ONLY: rho_ocean, sa
854
855	USE averaging, &
856	ONLY: rho_ocean_av, sa_av
857
858	USE control_parameters, &
859	ONLY: average_count_3d
860
861	USE indices, &
862	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt
863
864	IMPLICIT NONE
865
866	CHARACTER (LEN=*) :: mode !< flag defining mode 'allocate', 'sum' or 'average'
867	CHARACTER (LEN=*) :: variable !< name of variable
868
869	INTEGER(iwp) :: i !< loop index
870	INTEGER(iwp) :: j !< loop index
871	INTEGER(iwp) :: k !< loop index
872
873	IF ( mode == 'allocate' ) THEN
874
875	SELECT CASE ( TRIM( variable ) )
876
877	CASE ( 'rho_sea_water' )
878	IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
879	ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
880	ENDIF
881	rho_ocean_av = 0.0_wp
882
883	CASE ( 'sa' )
884	IF ( .NOT. ALLOCATED( sa_av ) ) THEN
885	ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
886	ENDIF
887	sa_av = 0.0_wp
888
889	CASE DEFAULT
890	CONTINUE
891
892	END SELECT
893
894	ELSEIF ( mode == 'sum' ) THEN
895
896	SELECT CASE ( TRIM( variable ) )
897
898	CASE ( 'rho_sea_water' )
899	IF ( ALLOCATED( rho_ocean_av ) ) THEN
900	DO i = nxlg, nxrg
901	DO j = nysg, nyng
902	DO k = nzb, nzt+1
903	rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) + &
904	rho_ocean(k,j,i)
905	ENDDO
906	ENDDO
907	ENDDO
908	ENDIF
909
910	CASE ( 'sa' )
911	IF ( ALLOCATED( sa_av ) ) THEN
912	DO i = nxlg, nxrg
913	DO j = nysg, nyng
914	DO k = nzb, nzt+1
915	sa_av(k,j,i) = sa_av(k,j,i) + sa(k,j,i)
916	ENDDO
917	ENDDO
918	ENDDO
919	ENDIF
920
921	CASE DEFAULT
922	CONTINUE
923
924	END SELECT
925
926	ELSEIF ( mode == 'average' ) THEN
927
928	SELECT CASE ( TRIM( variable ) )
929
930	CASE ( 'rho_sea_water' )
931	IF ( ALLOCATED( rho_ocean_av ) ) THEN
932	DO i = nxlg, nxrg
933	DO j = nysg, nyng
934	DO k = nzb, nzt+1
935	rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) / &
936	REAL( average_count_3d, KIND=wp )
937	ENDDO
938	ENDDO
939	ENDDO
940	ENDIF
941
942	CASE ( 'sa' )
943	IF ( ALLOCATED( sa_av ) ) THEN
944	DO i = nxlg, nxrg
945	DO j = nysg, nyng
946	DO k = nzb, nzt+1
947	sa_av(k,j,i) = sa_av(k,j,i) / &
948	REAL( average_count_3d, KIND=wp )
949	ENDDO
950	ENDDO
951	ENDDO
952	ENDIF
953
954	END SELECT
955
956	ENDIF
957
958	END SUBROUTINE ocean_3d_data_averaging
959
960
961	!------------------------------------------------------------------------------!
962	! Description:
963	! ------------
964	!> Define 2D output variables.
965	!------------------------------------------------------------------------------!
966	SUBROUTINE ocean_data_output_2d( av, variable, found, grid, mode, local_pf, &
967	nzb_do, nzt_do )
968
969	USE arrays_3d, &
970	ONLY: rho_ocean, sa
971
972	USE averaging, &
973	ONLY: rho_ocean_av, sa_av
974
975	USE indices, &
976	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, &
977	wall_flags_0
978
979	IMPLICIT NONE
980
981	CHARACTER (LEN=*) :: grid !< name of vertical grid
982	CHARACTER (LEN=*) :: mode !< either 'xy', 'xz' or 'yz'
983	CHARACTER (LEN=*) :: variable !< name of variable
984
985	INTEGER(iwp) :: av !< flag for (non-)average output
986	INTEGER(iwp) :: flag_nr !< number of masking flag
987	INTEGER(iwp) :: i !< loop index
988	INTEGER(iwp) :: j !< loop index
989	INTEGER(iwp) :: k !< loop index
990	INTEGER(iwp) :: nzb_do !< vertical output index (bottom)
991	INTEGER(iwp) :: nzt_do !< vertical output index (top)
992
993	LOGICAL :: found !< flag if output variable is found
994	LOGICAL :: resorted !< flag if output is already resorted
995
996	REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute
997
998	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local
999	!< array to which output data is resorted to
1000
1001	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable
1002
1003	found = .TRUE.
1004	resorted = .FALSE.
1005	!
1006	!-- Set masking flag for topography for not resorted arrays
1007	flag_nr = 0
1008
1009	SELECT CASE ( TRIM( variable ) )
1010
1011	CASE ( 'rho_sea_water_xy', 'rho_sea_water_xz', 'rho_sea_water_yz' )
1012	IF ( av == 0 ) THEN
1013	to_be_resorted => rho_ocean
1014	ELSE
1015	IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
1016	ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1017	rho_ocean_av = REAL( fill_value, KIND = wp )
1018	ENDIF
1019	to_be_resorted => rho_ocean_av
1020	ENDIF
1021
1022	CASE ( 'sa_xy', 'sa_xz', 'sa_yz' )
1023	IF ( av == 0 ) THEN
1024	to_be_resorted => sa
1025	ELSE
1026	IF ( .NOT. ALLOCATED( sa_av ) ) THEN
1027	ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1028	sa_av = REAL( fill_value, KIND = wp )
1029	ENDIF
1030	to_be_resorted => sa_av
1031	ENDIF
1032	IF ( mode == 'xy' ) grid = 'zu'
1033
1034	CASE DEFAULT
1035	found = .FALSE.
1036	grid = 'none'
1037
1038	END SELECT
1039
1040	IF ( found .AND. .NOT. resorted ) THEN
1041	DO i = nxl, nxr
1042	DO j = nys, nyn
1043	DO k = nzb_do, nzt_do
1044	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
1045	REAL( fill_value, KIND = wp ), &
1046	BTEST( wall_flags_0(k,j,i), flag_nr ) )
1047	ENDDO
1048	ENDDO
1049	ENDDO
1050	resorted = .TRUE.
1051	ENDIF
1052
1053	END SUBROUTINE ocean_data_output_2d
1054
1055
1056	!------------------------------------------------------------------------------!
1057	! Description:
1058	! ------------
1059	!> Define 3D output variables.
1060	!------------------------------------------------------------------------------!
1061	SUBROUTINE ocean_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
1062
1063
1064	USE arrays_3d, &
1065	ONLY: rho_ocean, sa
1066
1067	USE averaging, &
1068	ONLY: rho_ocean_av, sa_av
1069
1070	USE indices, &
1071	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, &
1072	wall_flags_0
1073
1074	IMPLICIT NONE
1075
1076	CHARACTER (LEN=*) :: variable !< name of variable
1077
1078	INTEGER(iwp) :: av !< flag for (non-)average output
1079	INTEGER(iwp) :: flag_nr !< number of masking flag
1080	INTEGER(iwp) :: i !< loop index
1081	INTEGER(iwp) :: j !< loop index
1082	INTEGER(iwp) :: k !< loop index
1083	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
1084	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
1085
1086	LOGICAL :: found !< flag if output variable is found
1087	LOGICAL :: resorted !< flag if output is already resorted
1088
1089	REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute
1090
1091	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local
1092	!< array to which output data is resorted to
1093
1094	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable
1095
1096	found = .TRUE.
1097	resorted = .FALSE.
1098	!
1099	!-- Set masking flag for topography for not resorted arrays
1100	flag_nr = 0
1101
1102	SELECT CASE ( TRIM( variable ) )
1103
1104	CASE ( 'rho_sea_water' )
1105	IF ( av == 0 ) THEN
1106	to_be_resorted => rho_ocean
1107	ELSE
1108	IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
1109	ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1110	rho_ocean_av = REAL( fill_value, KIND = wp )
1111	ENDIF
1112	to_be_resorted => rho_ocean_av
1113	ENDIF
1114
1115	CASE ( 'sa' )
1116	IF ( av == 0 ) THEN
1117	to_be_resorted => sa
1118	ELSE
1119	IF ( .NOT. ALLOCATED( sa_av ) ) THEN
1120	ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1121	sa_av = REAL( fill_value, KIND = wp )
1122	ENDIF
1123	to_be_resorted => sa_av
1124	ENDIF
1125
1126	CASE DEFAULT
1127	found = .FALSE.
1128
1129	END SELECT
1130
1131
1132	IF ( found .AND. .NOT. resorted ) THEN
1133	DO i = nxl, nxr
1134	DO j = nys, nyn
1135	DO k = nzb_do, nzt_do
1136	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
1137	REAL( fill_value, KIND = wp ), &
1138	BTEST( wall_flags_0(k,j,i), flag_nr ) )
1139	ENDDO
1140	ENDDO
1141	ENDDO
1142	resorted = .TRUE.
1143	ENDIF
1144
1145	END SUBROUTINE ocean_data_output_3d
1146
1147
1148	!------------------------------------------------------------------------------!
1149	! Description:
1150	! ------------
1151	!> Header output for ocean parameters
1152	!------------------------------------------------------------------------------!
1153	SUBROUTINE ocean_header( io )
1154
1155
1156	IMPLICIT NONE
1157
1158	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
1159
1160	!
1161	!-- Write ocean header
1162	WRITE( io, 1 )
1163	IF ( stokes_force ) WRITE( io, 2 ) stokes_waveheight, stokes_wavelength
1164	IF ( wave_breaking ) THEN
1165	WRITE( io, 3 ) alpha_wave_breaking, timescale_wave_breaking
1166	ENDIF
1167	IF ( .NOT. salinity ) WRITE( io, 4 )
1168	IF ( surface_cooling_spinup_time /= 999999.9_wp ) THEN
1169	WRITE( io, 5 ) surface_cooling_spinup_time
1170	ENDIF
1171
1172	1 FORMAT (//' Ocean settings:'/ &
1173	' ------------------------------------------'/)
1174	2 FORMAT (' --> Craik-Leibovich vortex force and Stokes drift switched', &
1175	' on'/ &
1176	' waveheight: ',F4.1,' m wavelength: ',F6.1,' m')
1177	3 FORMAT (' --> wave breaking generated turbulence switched on'/ &
1178	' alpha: ',F4.1/ &
1179	' timescale:',F5.1,' s')
1180	4 FORMAT (' --> prognostic salinity equation is switched off' )
1181	5 FORMAT (' --> surface heat flux is switched off after ',F8.1,' s')
1182
1183	END SUBROUTINE ocean_header
1184
1185
1186	!------------------------------------------------------------------------------!
1187	! Description:
1188	! ------------
1189	!> Allocate arrays and assign pointers.
1190	!------------------------------------------------------------------------------!
1191	SUBROUTINE ocean_init_arrays
1192
1193	USE indices, &
1194	ONLY: nxlg, nxrg, nyn, nyng, nys, nysg, nzb, nzt
1195
1196	IMPLICIT NONE
1197
1198	ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
1199	rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
1200	sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
1201	sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1202
1203	IF ( salinity ) THEN
1204	ALLOCATE( sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1205	ENDIF
1206
1207	IF ( ws_scheme_sca ) THEN
1208	ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) )
1209	sums_wssas_ws_l = 0.0_wp
1210	ENDIF
1211
1212	IF ( loop_optimization /= 'vector' ) THEN
1213
1214	IF ( ws_scheme_sca ) THEN
1215	ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1) )
1216	ALLOCATE( diss_s_sa(nzb+1:nzt,0:threads_per_task-1) )
1217	ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
1218	ALLOCATE( diss_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
1219	ENDIF
1220
1221	ENDIF
1222
1223	prho => prho_1
1224	rho_ocean => rho_1 ! routines calc_mean_profile and diffusion_e require
1225	! density to be a pointer
1226
1227	!
1228	!-- Initial assignment of pointers
1229	IF ( salinity ) THEN
1230	sa => sa_1; sa_p => sa_2; tsa_m => sa_3
1231	ELSE
1232	sa => sa_1; sa_p => sa_1; tsa_m => sa_3
1233	ENDIF
1234
1235	END SUBROUTINE ocean_init_arrays
1236
1237
1238	!------------------------------------------------------------------------------!
1239	! Description:
1240	! ------------
1241	!> Initialization of quantities needed for the ocean mode
1242	!------------------------------------------------------------------------------!
1243	SUBROUTINE ocean_init
1244
1245
1246	USE arrays_3d, &
1247	ONLY: dzu, dzw, hyp, pt_init, ref_state, u_stokes_zu, u_stokes_zw, &
1248	v_stokes_zu, v_stokes_zw, zu, zw
1249
1250	USE basic_constants_and_equations_mod, &
1251	ONLY: g
1252
1253	USE basic_constants_and_equations_mod, &
1254	ONLY: pi
1255
1256	USE control_parameters, &
1257	ONLY: initializing_actions, molecular_viscosity, rho_surface, &
1258	rho_reference, surface_pressure, top_momentumflux_u, &
1259	top_momentumflux_v, use_single_reference_value
1260
1261	USE indices, &
1262	ONLY: nxl, nxlg, nxrg, nyng, nys, nysg, nzb, nzt
1263
1264	USE kinds
1265
1266	USE pegrid, &
1267	ONLY: myid
1268
1269	USE statistics, &
1270	ONLY: hom, statistic_regions
1271
1272	IMPLICIT NONE
1273
1274	INTEGER(iwp) :: i !< loop index
1275	INTEGER(iwp) :: j !< loop index
1276	INTEGER(iwp) :: k !< loop index
1277	INTEGER(iwp) :: n !< loop index
1278
1279	REAL(wp) :: alpha !< angle of surface stress
1280	REAL(wp) :: dum !< dummy argument
1281	REAL(wp) :: pt_l !< local scalar for pt used in equation of state function
1282	REAL(wp) :: sa_l !< local scalar for sa used in equation of state function
1283	REAL(wp) :: velocity_amplitude !< local scalar for amplitude of Stokes drift velocity
1284	REAL(wp) :: x !< temporary variable to store surface stress along x
1285	REAL(wp) :: y !< temporary variable to store surface stress along y
1286
1287	REAL(wp), DIMENSION(nzb:nzt+1) :: rho_ocean_init !< local array for initial density
1288
1289	ALLOCATE( hyp(nzb:nzt+1) )
1290
1291
1292	!
1293	!-- In case of no restart run, calculate the inital salinity profilevcusing the
1294	!-- given salinity gradients
1295	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
1296
1297	sa_init = sa_surface
1298	!
1299	!-- Last arguments gives back the gradient at top level to be used as
1300	!-- possible Neumann boundary condition. This is not realized for the ocean
1301	!-- mode, therefore a dummy argument is used.
1302	IF ( salinity ) THEN
1303	CALL init_vertical_profiles( sa_vertical_gradient_level_ind, &
1304	sa_vertical_gradient_level, &
1305	sa_vertical_gradient, sa_init, &
1306	sa_surface, dum )
1307	ENDIF
1308	ENDIF
1309
1310	!
1311	!-- Initialize required 3d-arrays
1312	IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
1313	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
1314	!
1315	!-- Initialization via computed 1D-model profiles
1316	IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 ) THEN
1317
1318	DO i = nxlg, nxrg
1319	DO j = nysg, nyng
1320	sa(:,j,i) = sa_init
1321	ENDDO
1322	ENDDO
1323
1324	ENDIF
1325	!
1326	!-- Store initial profiles for output purposes etc.
1327	!-- Store initial salinity profile
1328	hom(:,1,26,:) = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
1329	!
1330	!-- Initialize old and new time levels.
1331	tsa_m = 0.0_wp
1332	sa_p = sa
1333
1334	ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN
1335
1336	!
1337	!-- Initialize new time levels (only done in order to set boundary values
1338	!-- including ghost points)
1339	sa_p = sa
1340	!
1341	!-- Allthough tendency arrays are set in prognostic_equations, they have
1342	!-- have to be predefined here because they are used (but multiplied with 0)
1343	!-- there before they are set.
1344	tsa_m = 0.0_wp
1345
1346	ENDIF
1347
1348	!
1349	!-- Set water density near the ocean surface
1350	rho_surface = 1027.62_wp
1351
1352	!
1353	!-- Set kinematic viscosity to sea water at 20C.
1354	!-- This changes the default value that is given for air!
1355	molecular_viscosity = 1.05E-6_wp
1356
1357	!
1358	!-- Change sign of buoyancy/stability terms because density gradient is used
1359	!-- instead of the potential temperature gradient to calculate the buoyancy
1360	atmos_ocean_sign = -1.0_wp
1361
1362	!
1363	!-- Calculate initial vertical profile of hydrostatic pressure (in Pa)
1364	!-- and the reference density (used later in buoyancy term)
1365	!-- First step: Calculate pressure using reference density
1366	hyp(nzt+1) = surface_pressure * 100.0_wp
1367	hyp(nzt) = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1)
1368	rho_ocean_init(nzt) = rho_surface
1369	rho_ocean_init(nzt+1) = rho_surface ! only required for output
1370
1371	DO k = nzt-1, 1, -1
1372	hyp(k) = hyp(k+1) + rho_surface * g * dzu(k)
1373	ENDDO
1374	hyp(0) = hyp(1) + rho_surface * g * dzu(1)
1375
1376	!
1377	!-- Second step: Iteratively calculate in situ density (based on presssure)
1378	!-- and pressure (based on in situ density)
1379	DO n = 1, 5
1380
1381	rho_reference = rho_surface * 0.5_wp * dzu(nzt+1)
1382
1383	DO k = nzt, 0, -1
1384
1385	sa_l = 0.5_wp * ( sa_init(k) + sa_init(k+1) )
1386	pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) )
1387
1388	rho_ocean_init(k) = eqn_state_seawater_func( hyp(k), pt_l, sa_l )
1389
1390	rho_reference = rho_reference + rho_ocean_init(k) * dzu(k+1)
1391
1392	ENDDO
1393
1394	rho_reference = rho_reference / ( zw(nzt) - zu(nzb) )
1395
1396	hyp(nzt) = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1)
1397	DO k = nzt-1, 0, -1
1398	hyp(k) = hyp(k+1) + g * 0.5_wp * ( rho_ocean_init(k) &
1399	+ rho_ocean_init(k+1) ) * dzu(k+1)
1400	ENDDO
1401
1402	ENDDO
1403
1404	!
1405	!-- Calculate the reference potential density
1406	prho_reference = 0.0_wp
1407	DO k = 0, nzt
1408
1409	sa_l = 0.5_wp * ( sa_init(k) + sa_init(k+1) )
1410	pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) )
1411
1412	prho_reference = prho_reference + dzu(k+1) * &
1413	eqn_state_seawater_func( 0.0_wp, pt_l, sa_l )
1414
1415	ENDDO
1416
1417	prho_reference = prho_reference / ( zu(nzt) - zu(nzb) )
1418
1419	!
1420	!-- Calculate the 3d array of initial in situ and potential density,
1421	!-- based on the initial temperature and salinity profile
1422	CALL eqn_state_seawater
1423
1424	!
1425	!-- Store initial density profile
1426	hom(:,1,77,:) = SPREAD( rho_ocean_init(:), 2, statistic_regions+1 )
1427
1428	!
1429	!-- Set the reference state to be used in the buoyancy terms
1430	IF ( use_single_reference_value ) THEN
1431	ref_state(:) = rho_reference
1432	ELSE
1433	ref_state(:) = rho_ocean_init(:)
1434	ENDIF
1435
1436	!
1437	!-- Calculate the Stokes drift velocity profile
1438	IF ( stokes_force ) THEN
1439
1440	!
1441	!-- First, calculate angle of surface stress
1442	x = -top_momentumflux_u
1443	y = -top_momentumflux_v
1444	IF ( x == 0.0_wp ) THEN
1445	IF ( y > 0.0_wp ) THEN
1446	alpha = pi / 2.0_wp
1447	ELSEIF ( y < 0.0_wp ) THEN
1448	alpha = 3.0_wp * pi / 2.0_wp
1449	ENDIF
1450	ELSE
1451	IF ( x < 0.0_wp ) THEN
1452	alpha = ATAN( y / x ) + pi
1453	ELSE
1454	IF ( y < 0.0_wp ) THEN
1455	alpha = ATAN( y / x ) + 2.0_wp * pi
1456	ELSE
1457	alpha = ATAN( y / x )
1458	ENDIF
1459	ENDIF
1460	ENDIF
1461
1462	velocity_amplitude = ( pi * stokes_waveheight / stokes_wavelength )*2 &
1463	SQRT( g * stokes_wavelength / ( 2.0_wp * pi ) )
1464
1465	DO k = nzb, nzt
1466	u_stokes_zu(k) = velocity_amplitude * COS( alpha ) * &
1467	EXP( 4.0_wp * pi * zu(k) / stokes_wavelength )
1468	u_stokes_zw(k) = velocity_amplitude * COS( alpha ) * &
1469	EXP( 4.0_wp * pi * zw(k) / stokes_wavelength )
1470	v_stokes_zu(k) = velocity_amplitude * SIN( alpha ) * &
1471	EXP( 4.0_wp * pi * zu(k) / stokes_wavelength )
1472	v_stokes_zw(k) = velocity_amplitude * SIN( alpha ) * &
1473	EXP( 4.0_wp * pi * zw(k) / stokes_wavelength )
1474	ENDDO
1475	u_stokes_zu(nzt+1) = u_stokes_zw(nzt) ! because zu(nzt+1) changes the sign
1476	u_stokes_zw(nzt+1) = u_stokes_zw(nzt) ! because zw(nzt+1) changes the sign
1477	v_stokes_zu(nzt+1) = v_stokes_zw(nzt) ! because zu(nzt+1) changes the sign
1478	v_stokes_zw(nzt+1) = v_stokes_zw(nzt) ! because zw(nzt+1) changes the sign
1479
1480	ENDIF
1481
1482	!
1483	!-- Wave breaking effects
1484	IF ( wave_breaking ) THEN
1485	!
1486	!-- Calculate friction velocity at ocean surface
1487	u_star_wave_breaking = SQRT( SQRT( top_momentumflux_u**2 + &
1488	top_momentumflux_v**2 ) )
1489	!
1490	!-- Set the time scale of random forcing. The vertical grid spacing at the
1491	!-- ocean surface is assumed as the length scale of turbulence.
1492	!-- Formula follows Noh et al. (2004), JPO
1493	timescale_wave_breaking = 0.1_wp * dzw(nzt) / alpha_wave_breaking / &
1494	u_star_wave_breaking
1495	!
1496	!-- Set random number seeds differently on the processor cores in order to
1497	!-- create different random number sequences
1498	iran_ocean = iran_ocean + myid
1499	ENDIF
1500
1501	END SUBROUTINE ocean_init
1502
1503
1504	!------------------------------------------------------------------------------!
1505	! Description:
1506	! ------------
1507	!> Call for all grid points
1508	!------------------------------------------------------------------------------!
1509	SUBROUTINE ocean_actions( location )
1510
1511
1512	CHARACTER (LEN=*), INTENT(IN) :: location !< call location string
1513
1514	SELECT CASE ( location )
1515
1516	CASE ( 'before_timestep' )
1517
1518	IF ( ws_scheme_sca ) THEN
1519	sums_wssas_ws_l = 0.0_wp
1520	ENDIF
1521
1522	CASE DEFAULT
1523	CONTINUE
1524
1525	END SELECT
1526
1527	END SUBROUTINE ocean_actions
1528
1529
1530	!------------------------------------------------------------------------------!
1531	! Description:
1532	! ------------
1533	!> Call for grid points i,j
1534	!------------------------------------------------------------------------------!
1535	SUBROUTINE ocean_actions_ij( i, j, location )
1536
1537
1538	INTEGER(iwp), INTENT(IN) :: i !< grid index in x-direction
1539	INTEGER(iwp), INTENT(IN) :: j !< grid index in y-direction
1540	CHARACTER (LEN=*), INTENT(IN) :: location !< call location string
1541	INTEGER(iwp) :: dummy !< call location string
1542
1543	IF ( ocean_mode ) dummy = i + j
1544
1545	SELECT CASE ( location )
1546
1547	CASE ( 'before_timestep' )
1548
1549	IF ( ws_scheme_sca ) THEN
1550	sums_wssas_ws_l = 0.0_wp
1551	ENDIF
1552
1553	CASE DEFAULT
1554	CONTINUE
1555
1556	END SELECT
1557
1558	END SUBROUTINE ocean_actions_ij
1559
1560
1561	!------------------------------------------------------------------------------!
1562	! Description:
1563	! ------------
1564	!> Prognostic equation for salinity.
1565	!> Vector-optimized version
1566	!------------------------------------------------------------------------------!
1567	SUBROUTINE ocean_prognostic_equations
1568
1569	USE advec_s_bc_mod, &
1570	ONLY: advec_s_bc
1571
1572	USE advec_s_pw_mod, &
1573	ONLY: advec_s_pw
1574
1575	USE advec_s_up_mod, &
1576	ONLY: advec_s_up
1577
1578	USE advec_ws, &
1579	ONLY: advec_s_ws
1580
1581	USE arrays_3d, &
1582	ONLY: rdf_sc, tend, tsa_m
1583
1584	USE control_parameters, &
1585	ONLY: bc_dirichlet_l, &
1586	bc_dirichlet_n, &
1587	bc_dirichlet_r, &
1588	bc_dirichlet_s, &
1589	bc_radiation_l, &
1590	bc_radiation_n, &
1591	bc_radiation_r, &
1592	bc_radiation_s, &
1593	dt_3d, intermediate_timestep_count, &
1594	intermediate_timestep_count_max, scalar_advec, simulated_time, &
1595	timestep_scheme, tsc, ws_scheme_sca
1596
1597	USE cpulog, &
1598	ONLY: cpu_log, log_point_s
1599
1600	USE diffusion_s_mod, &
1601	ONLY: diffusion_s
1602
1603	USE indices, &
1604	ONLY: advc_flags_s, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
1605
1606	USE surface_mod, &
1607	ONLY: surf_def_v, surf_def_h, surf_lsm_h, surf_lsm_v, surf_usm_h, &
1608	surf_usm_v
1609
1610	IMPLICIT NONE
1611
1612	INTEGER(iwp) :: i !< loop index
1613	INTEGER(iwp) :: j !< loop index
1614	INTEGER(iwp) :: k !< loop index
1615
1616	REAL(wp) :: sbt !< weighting factor for sub-time step
1617
1618	!
1619	!-- Switch of the surface heat flux, if requested
1620	IF ( surface_cooling_spinup_time /= 999999.9_wp ) THEN
1621	IF ( .NOT. surface_cooling_switched_off .AND. &
1622	simulated_time >= surface_cooling_spinup_time ) THEN
1623
1624	surf_def_h(2)%shf = 0.0_wp
1625	surface_cooling_switched_off = .TRUE.
1626
1627	ENDIF
1628	ENDIF
1629
1630	!
1631	!-- Compute prognostic equations for the ocean mode
1632	!-- First, start with salinity
1633	IF ( salinity ) THEN
1634
1635	CALL cpu_log( log_point_s(20), 'sa-equation', 'start' )
1636
1637	!
1638	!-- sa-tendency terms with communication
1639	sbt = tsc(2)
1640	IF ( scalar_advec == 'bc-scheme' ) THEN
1641
1642	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1643	!
1644	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1645	sbt = 1.0_wp
1646	ENDIF
1647	tend = 0.0_wp
1648	CALL advec_s_bc( sa, 'sa' )
1649
1650	ENDIF
1651
1652	!
1653	!-- sa-tendency terms with no communication
1654	IF ( scalar_advec /= 'bc-scheme' ) THEN
1655	tend = 0.0_wp
1656	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1657	IF ( ws_scheme_sca ) THEN
1658	CALL advec_s_ws( advc_flags_s, sa, 'sa', &
1659	bc_dirichlet_l .OR. bc_radiation_l, &
1660	bc_dirichlet_n .OR. bc_radiation_n, &
1661	bc_dirichlet_r .OR. bc_radiation_r, &
1662	bc_dirichlet_s .OR. bc_radiation_s )
1663	ELSE
1664	CALL advec_s_pw( sa )
1665	ENDIF
1666	ELSE
1667	CALL advec_s_up( sa )
1668	ENDIF
1669	ENDIF
1670
1671	CALL diffusion_s( sa, &
1672	surf_def_h(0)%sasws, surf_def_h(1)%sasws, &
1673	surf_def_h(2)%sasws, &
1674	surf_lsm_h%sasws, surf_usm_h%sasws, &
1675	surf_def_v(0)%sasws, surf_def_v(1)%sasws, &
1676	surf_def_v(2)%sasws, surf_def_v(3)%sasws, &
1677	surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws, &
1678	surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws, &
1679	surf_usm_v(0)%sasws, surf_usm_v(1)%sasws, &
1680	surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
1681
1682	! CALL user_actions( 'sa-tendency' ) ToDo: find general solution for dependency between modules
1683
1684	!
1685	!-- Prognostic equation for salinity
1686	DO i = nxl, nxr
1687	DO j = nys, nyn
1688	DO k = nzb+1, nzt
1689	sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
1690	tsc(3) * tsa_m(k,j,i) ) &
1691	- tsc(5) * rdf_sc(k) * &
1692	( sa(k,j,i) - sa_init(k) ) &
1693	) &
1694	* MERGE( 1.0_wp, 0.0_wp, &
1695	BTEST( wall_flags_0(k,j,i), 0 ) &
1696	)
1697	IF ( sa_p(k,j,i) < 0.0_wp ) sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
1698	ENDDO
1699	ENDDO
1700	ENDDO
1701
1702	!
1703	!-- Calculate tendencies for the next Runge-Kutta step
1704	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1705	IF ( intermediate_timestep_count == 1 ) THEN
1706	DO i = nxl, nxr
1707	DO j = nys, nyn
1708	DO k = nzb+1, nzt
1709	tsa_m(k,j,i) = tend(k,j,i)
1710	ENDDO
1711	ENDDO
1712	ENDDO
1713	ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) &
1714	THEN
1715	DO i = nxl, nxr
1716	DO j = nys, nyn
1717	DO k = nzb+1, nzt
1718	tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1719	5.3125_wp * tsa_m(k,j,i)
1720	ENDDO
1721	ENDDO
1722	ENDDO
1723	ENDIF
1724	ENDIF
1725
1726	CALL cpu_log( log_point_s(20), 'sa-equation', 'stop' )
1727
1728	ENDIF
1729
1730	!
1731	!-- Calculate density by the equation of state for seawater
1732	CALL cpu_log( log_point_s(21), 'eqns-seawater', 'start' )
1733	CALL eqn_state_seawater
1734	CALL cpu_log( log_point_s(21), 'eqns-seawater', 'stop' )
1735
1736	END SUBROUTINE ocean_prognostic_equations
1737
1738
1739	!------------------------------------------------------------------------------!
1740	! Description:
1741	! ------------
1742	!> Prognostic equations for ocean mode (so far, salinity only)
1743	!> Cache-optimized version
1744	!------------------------------------------------------------------------------!
1745	SUBROUTINE ocean_prognostic_equations_ij( i, j, i_omp_start, tn )
1746
1747	USE advec_s_pw_mod, &
1748	ONLY: advec_s_pw
1749
1750	USE advec_s_up_mod, &
1751	ONLY: advec_s_up
1752
1753	USE advec_ws, &
1754	ONLY: advec_s_ws
1755
1756	USE arrays_3d, &
1757	ONLY: diss_l_sa, diss_s_sa, flux_l_sa, flux_s_sa, rdf_sc, tend, tsa_m
1758
1759	USE control_parameters, &
1760	ONLY: bc_dirichlet_l, &
1761	bc_dirichlet_n, &
1762	bc_dirichlet_r, &
1763	bc_dirichlet_s, &
1764	bc_radiation_l, &
1765	bc_radiation_n, &
1766	bc_radiation_r, &
1767	bc_radiation_s, &
1768	dt_3d, intermediate_timestep_count, &
1769	intermediate_timestep_count_max, simulated_time, &
1770	timestep_scheme, tsc, ws_scheme_sca
1771
1772	USE diffusion_s_mod, &
1773	ONLY: diffusion_s
1774
1775	USE indices, &
1776	ONLY: advc_flags_s, nzb, nzt, wall_flags_0
1777
1778	USE surface_mod, &
1779	ONLY: surf_def_v, surf_def_h, surf_lsm_h, surf_lsm_v, surf_usm_h, &
1780	surf_usm_v
1781
1782	IMPLICIT NONE
1783
1784	INTEGER(iwp) :: i !< loop index x direction
1785	INTEGER(iwp) :: i_omp_start !< first loop index of i-loop in calling &
1786	!< routine prognostic_equations
1787	INTEGER(iwp) :: j !< loop index y direction
1788	INTEGER(iwp) :: k !< loop index z direction
1789	INTEGER(iwp) :: tn !< task number of openmp task
1790
1791
1792	!
1793	!-- Switch of the surface heat flux, if requested
1794	IF ( surface_cooling_spinup_time /= 999999.9_wp ) THEN
1795	IF ( .NOT. surface_cooling_switched_off .AND. &
1796	simulated_time >= surface_cooling_spinup_time ) THEN
1797
1798	surf_def_h(2)%shf = 0.0_wp
1799	surface_cooling_switched_off = .TRUE.
1800
1801	ENDIF
1802	ENDIF
1803
1804	!
1805	!-- Compute prognostic equations for the ocean mode
1806	!-- First, start with tendency-terms for salinity
1807	IF ( salinity ) THEN
1808
1809	tend(:,j,i) = 0.0_wp
1810	IF ( timestep_scheme(1:5) == 'runge' ) &
1811	THEN
1812	IF ( ws_scheme_sca ) THEN
1813	CALL advec_s_ws( advc_flags_s, &
1814	i, j, sa, 'sa', flux_s_sa, diss_s_sa, flux_l_sa,&
1815	diss_l_sa, i_omp_start, tn, &
1816	bc_dirichlet_l .OR. bc_radiation_l, &
1817	bc_dirichlet_n .OR. bc_radiation_n, &
1818	bc_dirichlet_r .OR. bc_radiation_r, &
1819	bc_dirichlet_s .OR. bc_radiation_s )
1820	ELSE
1821	CALL advec_s_pw( i, j, sa )
1822	ENDIF
1823	ELSE
1824	CALL advec_s_up( i, j, sa )
1825	ENDIF
1826	CALL diffusion_s( i, j, sa, &
1827	surf_def_h(0)%sasws, surf_def_h(1)%sasws, &
1828	surf_def_h(2)%sasws, &
1829	surf_lsm_h%sasws, surf_usm_h%sasws, &
1830	surf_def_v(0)%sasws, surf_def_v(1)%sasws, &
1831	surf_def_v(2)%sasws, surf_def_v(3)%sasws, &
1832	surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws, &
1833	surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws, &
1834	surf_usm_v(0)%sasws, surf_usm_v(1)%sasws, &
1835	surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
1836
1837	! CALL user_actions( i, j, 'sa-tendency' ) ToDo: find general solution for dependency between modules
1838
1839	!
1840	!-- Prognostic equation for salinity
1841	DO k = nzb+1, nzt
1842
1843	sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * &
1844	( tsc(2) * tend(k,j,i) + &
1845	tsc(3) * tsa_m(k,j,i) ) &
1846	- tsc(5) * rdf_sc(k) &
1847	* ( sa(k,j,i) - sa_init(k) ) &
1848	) * MERGE( 1.0_wp, 0.0_wp, &
1849	BTEST( wall_flags_0(k,j,i), 0 ) )
1850
1851	IF ( sa_p(k,j,i) < 0.0_wp ) sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
1852
1853	ENDDO
1854
1855	!
1856	!-- Calculate tendencies for the next Runge-Kutta step
1857	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1858	IF ( intermediate_timestep_count == 1 ) THEN
1859	DO k = nzb+1, nzt
1860	tsa_m(k,j,i) = tend(k,j,i)
1861	ENDDO
1862	ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) &
1863	THEN
1864	DO k = nzb+1, nzt
1865	tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1866	5.3125_wp * tsa_m(k,j,i)
1867	ENDDO
1868	ENDIF
1869	ENDIF
1870
1871	ENDIF
1872
1873	!
1874	!-- Calculate density by the equation of state for seawater
1875	CALL eqn_state_seawater( i, j )
1876
1877	END SUBROUTINE ocean_prognostic_equations_ij
1878
1879
1880	!------------------------------------------------------------------------------!
1881	! Description:
1882	! ------------
1883	!> Swapping of timelevels.
1884	!------------------------------------------------------------------------------!
1885	SUBROUTINE ocean_swap_timelevel( mod_count )
1886
1887	IMPLICIT NONE
1888
1889	INTEGER, INTENT(IN) :: mod_count !< flag defining where pointers point to
1890
1891
1892	SELECT CASE ( mod_count )
1893
1894	CASE ( 0 )
1895	IF ( salinity ) THEN
1896	sa => sa_1; sa_p => sa_2
1897	ENDIF
1898
1899	CASE ( 1 )
1900	IF ( salinity ) THEN
1901	sa => sa_2; sa_p => sa_1
1902	ENDIF
1903
1904	END SELECT
1905
1906	END SUBROUTINE ocean_swap_timelevel
1907
1908
1909	!------------------------------------------------------------------------------!
1910	! Description:
1911	! ------------
1912	!> This routine reads the respective restart data for the ocean module.
1913	!------------------------------------------------------------------------------!
1914	SUBROUTINE ocean_rrd_global( found )
1915
1916
1917	USE control_parameters, &
1918	ONLY: length, restart_string
1919
1920
1921	IMPLICIT NONE
1922
1923	LOGICAL, INTENT(OUT) :: found
1924
1925
1926	found = .TRUE.
1927
1928	SELECT CASE ( restart_string(1:length) )
1929
1930	CASE ( 'bc_sa_t' )
1931	READ ( 13 ) bc_sa_t
1932
1933	CASE ( 'bottom_salinityflux' )
1934	READ ( 13 ) bottom_salinityflux
1935
1936	CASE ( 'salinity' )
1937	READ ( 13 ) salinity
1938
1939	CASE ( 'sa_init' )
1940	READ ( 13 ) sa_init
1941
1942	CASE ( 'sa_surface' )
1943	READ ( 13 ) sa_surface
1944
1945	CASE ( 'sa_vertical_gradient' )
1946	READ ( 13 ) sa_vertical_gradient
1947
1948	CASE ( 'sa_vertical_gradient_level' )
1949	READ ( 13 ) sa_vertical_gradient_level
1950
1951	CASE ( 'stokes_waveheight' )
1952	READ ( 13 ) stokes_waveheight
1953
1954	CASE ( 'stokes_wavelength' )
1955	READ ( 13 ) stokes_wavelength
1956
1957	CASE ( 'surface_cooling_spinup_time' )
1958	READ ( 13 ) surface_cooling_spinup_time
1959
1960	CASE ( 'top_salinityflux' )
1961	READ ( 13 ) top_salinityflux
1962
1963	CASE ( 'wall_salinityflux' )
1964	READ ( 13 ) wall_salinityflux
1965
1966	CASE ( 'wave_breaking' )
1967	READ ( 13 ) wave_breaking
1968
1969	CASE DEFAULT
1970
1971	found = .FALSE.
1972
1973	END SELECT
1974
1975	END SUBROUTINE ocean_rrd_global
1976
1977
1978	!------------------------------------------------------------------------------!
1979	! Description:
1980	! ------------
1981	!> This routine reads the respective restart data for the ocean module.
1982	!------------------------------------------------------------------------------!
1983	SUBROUTINE ocean_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, &
1984	nxr_on_file, nynf, nync, nyn_on_file, nysf, &
1985	nysc, nys_on_file, tmp_3d, found )
1986
1987	USE averaging, &
1988	ONLY: rho_ocean_av, sa_av
1989
1990	USE control_parameters, &
1991	ONLY: length, restart_string
1992
1993	USE indices, &
1994	ONLY: nbgp, nxlg, nxrg, nyng, nysg, nzb, nzt
1995
1996	USE pegrid
1997
1998
1999	IMPLICIT NONE
2000
2001	INTEGER(iwp) :: k !<
2002	INTEGER(iwp) :: nxlc !<
2003	INTEGER(iwp) :: nxlf !<
2004	INTEGER(iwp) :: nxl_on_file !<
2005	INTEGER(iwp) :: nxrc !<
2006	INTEGER(iwp) :: nxrf !<
2007	INTEGER(iwp) :: nxr_on_file !<
2008	INTEGER(iwp) :: nync !<
2009	INTEGER(iwp) :: nynf !<
2010	INTEGER(iwp) :: nyn_on_file !<
2011	INTEGER(iwp) :: nysc !<
2012	INTEGER(iwp) :: nysf !<
2013	INTEGER(iwp) :: nys_on_file !<
2014
2015	LOGICAL, INTENT(OUT) :: found
2016
2017	REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !<
2018
2019
2020	found = .TRUE.
2021
2022	SELECT CASE ( restart_string(1:length) )
2023
2024	CASE ( 'rho_ocean_av' )
2025	IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
2026	ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2027	ENDIF
2028	IF ( k == 1 ) READ ( 13 ) tmp_3d
2029	rho_ocean_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2030	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2031
2032	CASE ( 'sa' )
2033	IF ( k == 1 ) READ ( 13 ) tmp_3d
2034	sa(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2035	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2036
2037	CASE ( 'sa_av' )
2038	IF ( .NOT. ALLOCATED( sa_av ) ) THEN
2039	ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2040	ENDIF
2041	IF ( k == 1 ) READ ( 13 ) tmp_3d
2042	sa_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2043	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2044
2045	CASE DEFAULT
2046	found = .FALSE.
2047
2048	END SELECT
2049
2050	END SUBROUTINE ocean_rrd_local
2051
2052
2053	!------------------------------------------------------------------------------!
2054	! Description:
2055	! ------------
2056	!> This routine writes the respective restart data for the ocean module.
2057	!------------------------------------------------------------------------------!
2058	SUBROUTINE ocean_wrd_global
2059
2060
2061	IMPLICIT NONE
2062
2063	CALL wrd_write_string( 'bc_sa_t' )
2064	WRITE ( 14 ) bc_sa_t
2065
2066	CALL wrd_write_string( 'bottom_salinityflux' )
2067	WRITE ( 14 ) bottom_salinityflux
2068
2069	CALL wrd_write_string( 'salinity' )
2070	WRITE ( 14 ) salinity
2071
2072	CALL wrd_write_string( 'sa_init' )
2073	WRITE ( 14 ) sa_init
2074
2075	CALL wrd_write_string( 'sa_surface' )
2076	WRITE ( 14 ) sa_surface
2077
2078	CALL wrd_write_string( 'sa_vertical_gradient' )
2079	WRITE ( 14 ) sa_vertical_gradient
2080
2081	CALL wrd_write_string( 'sa_vertical_gradient_level' )
2082	WRITE ( 14 ) sa_vertical_gradient_level
2083
2084	CALL wrd_write_string( 'stokes_waveheight' )
2085	WRITE ( 14 ) stokes_waveheight
2086
2087	CALL wrd_write_string( 'stokes_wavelength' )
2088	WRITE ( 14 ) stokes_wavelength
2089
2090	CALL wrd_write_string( 'surface_cooling_spinup_time' )
2091	WRITE ( 14 ) surface_cooling_spinup_time
2092
2093	CALL wrd_write_string( 'top_salinityflux' )
2094	WRITE ( 14 ) top_salinityflux
2095
2096	CALL wrd_write_string( 'wall_salinityflux' )
2097	WRITE ( 14 ) wall_salinityflux
2098
2099	CALL wrd_write_string( 'wave_breaking' )
2100	WRITE ( 14 ) wave_breaking
2101
2102	END SUBROUTINE ocean_wrd_global
2103
2104
2105	!------------------------------------------------------------------------------!
2106	! Description:
2107	! ------------
2108	!> This routine writes the respective restart data for the ocean module.
2109	!------------------------------------------------------------------------------!
2110	SUBROUTINE ocean_wrd_local
2111
2112	USE averaging, &
2113	ONLY: rho_ocean_av, sa_av
2114
2115	IMPLICIT NONE
2116
2117	IF ( ALLOCATED( rho_ocean_av ) ) THEN
2118	CALL wrd_write_string( 'rho_ocean_av' )
2119	WRITE ( 14 ) rho_ocean_av
2120	ENDIF
2121
2122	CALL wrd_write_string( 'sa' )
2123	WRITE ( 14 ) sa
2124
2125	IF ( ALLOCATED( sa_av ) ) THEN
2126	CALL wrd_write_string( 'sa_av' )
2127	WRITE ( 14 ) sa_av
2128	ENDIF
2129
2130	END SUBROUTINE ocean_wrd_local
2131
2132
2133	!------------------------------------------------------------------------------!
2134	! Description:
2135	! ------------
2136	!> This routine calculates the Craik Leibovich vortex force and the additional
2137	!> effect of the Stokes drift on the Coriolis force
2138	!> Call for all gridpoints.
2139	!------------------------------------------------------------------------------!
2140	SUBROUTINE stokes_drift_terms( component )
2141
2142	USE arrays_3d, &
2143	ONLY: ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu, &
2144	v_stokes_zw, w, tend
2145
2146	USE control_parameters, &
2147	ONLY: f, fs, message_string
2148
2149	USE grid_variables, &
2150	ONLY: ddx, ddy
2151
2152	USE indices, &
2153	ONLY: nxl, nxr, nys, nysv, nyn, nzb, nzt
2154
2155	IMPLICIT NONE
2156
2157	INTEGER(iwp) :: component !< component of momentum equation
2158	INTEGER(iwp) :: i !< loop index along x
2159	INTEGER(iwp) :: j !< loop index along y
2160	INTEGER(iwp) :: k !< loop index along z
2161
2162
2163	!
2164	!-- Compute Stokes terms for the respective velocity components
2165	SELECT CASE ( component )
2166
2167	!
2168	!-- u-component
2169	CASE ( 1 )
2170	DO i = nxl, nxr
2171	DO j = nysv, nyn
2172	DO k = nzb+1, nzt
2173	tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * ( &
2174	0.5 * ( v(k,j+1,i) - v(k,j+1,i-1) &
2175	+ v(k,j,i) - v(k,j,i-1) ) * ddx &
2176	- 0.5 * ( u(k,j+1,i) - u(k,j-1,i) ) * ddy &
2177	) &
2178	+ f * v_stokes_zu(k)
2179	ENDDO
2180	ENDDO
2181	ENDDO
2182
2183	!
2184	!-- v-component
2185	CASE ( 2 )
2186	DO i = nxl, nxr
2187	DO j = nysv, nyn
2188	DO k = nzb+1, nzt
2189	tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * ( &
2190	0.5 * ( v(k,j,i+1) - v(k,j,i-1) ) * ddx &
2191	- 0.5 * ( u(k,j,i) - u(k,j-1,i) &
2192	+ u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
2193	) &
2194	- f * u_stokes_zu(k)
2195	ENDDO
2196	ENDDO
2197	ENDDO
2198
2199	!
2200	!-- w-component
2201	CASE ( 3 )
2202	DO i = nxl, nxr
2203	DO j = nys, nyn
2204	DO k = nzb+1, nzt
2205	tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * ( &
2206	0.5 * ( u(k+1,j,i) - u(k,j,i) &
2207	+ u(k+1,j,i+1) - u(k,j,i+1) &
2208	) * ddzu(k+1) &
2209	- 0.5 * ( w(k,j,i+1) - w(k,j,i-1) &
2210	) * ddx ) &
2211	- v_stokes_zw(k) * ( &
2212	0.5 * ( w(k,j+1,i) - w(k,j-1,i) &
2213	) * ddy &
2214	- 0.5 * ( v(k+1,j,i) - v(k,j,i) &
2215	+ v(k+1,j+1,i) - v(k,j+1,i) &
2216	) * ddzu(k) ) &
2217	+ fs * u_stokes_zw(k)
2218	ENDDO
2219	ENDDO
2220	ENDDO
2221
2222	CASE DEFAULT
2223	WRITE( message_string, * ) 'wrong component of Stokes force: ', &
2224	component
2225	CALL message( 'stokes_drift_terms', 'PA0091', 1, 2, 0, 6, 0 )
2226
2227	END SELECT
2228
2229	END SUBROUTINE stokes_drift_terms
2230
2231
2232	!------------------------------------------------------------------------------!
2233	! Description:
2234	! ------------
2235	!> This routine calculates the Craik Leibovich vortex force and the additional
2236	!> effect of the Stokes drift on the Coriolis force
2237	!> Call for gridpoints i,j.
2238	!------------------------------------------------------------------------------!
2239
2240	SUBROUTINE stokes_drift_terms_ij( i, j, component )
2241
2242	USE arrays_3d, &
2243	ONLY: ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu, &
2244	v_stokes_zw, w, tend
2245
2246	USE control_parameters, &
2247	ONLY: f, fs, message_string
2248
2249	USE grid_variables, &
2250	ONLY: ddx, ddy
2251
2252	USE indices, &
2253	ONLY: nzb, nzt
2254
2255	IMPLICIT NONE
2256
2257	INTEGER(iwp) :: component !< component of momentum equation
2258	INTEGER(iwp) :: i !< loop index along x
2259	INTEGER(iwp) :: j !< loop index along y
2260	INTEGER(iwp) :: k !< loop incex along z
2261
2262
2263	!
2264	!-- Compute Stokes terms for the respective velocity components
2265	SELECT CASE ( component )
2266
2267	!
2268	!-- u-component
2269	CASE ( 1 )
2270	DO k = nzb+1, nzt
2271	tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * ( &
2272	0.5 * ( v(k,j+1,i) - v(k,j+1,i-1) &
2273	+ v(k,j,i) - v(k,j,i-1) ) * ddx &
2274	- 0.5 * ( u(k,j+1,i) - u(k,j-1,i) ) * ddy &
2275	) &
2276	+ f * v_stokes_zu(k)
2277	ENDDO
2278	!
2279	!-- v-component
2280	CASE ( 2 )
2281	DO k = nzb+1, nzt
2282	tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * ( &
2283	0.5 * ( v(k,j,i+1) - v(k,j,i-1) ) * ddx &
2284	- 0.5 * ( u(k,j,i) - u(k,j-1,i) &
2285	+ u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
2286	) &
2287	- f * u_stokes_zu(k)
2288	ENDDO
2289
2290	!
2291	!-- w-component
2292	CASE ( 3 )
2293	DO k = nzb+1, nzt
2294	tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * ( &
2295	0.5 * ( u(k+1,j,i) - u(k,j,i) &
2296	+ u(k+1,j,i+1) - u(k,j,i+1) &
2297	) * ddzu(k+1) &
2298	- 0.5 * ( w(k,j,i+1) - w(k,j,i-1) &
2299	) * ddx ) &
2300	- v_stokes_zw(k) * ( &
2301	0.5 * ( w(k,j+1,i) - w(k,j-1,i) &
2302	) * ddy &
2303	- 0.5 * ( v(k+1,j,i) - v(k,j,i) &
2304	+ v(k+1,j+1,i) - v(k,j+1,i) &
2305	) * ddzu(k) ) &
2306	+ fs * u_stokes_zw(k)
2307	ENDDO
2308
2309	CASE DEFAULT
2310	WRITE( message_string, * ) ' wrong component: ', component
2311	CALL message( 'stokes_drift_terms', 'PA0091', 1, 2, 0, 6, 0 )
2312
2313	END SELECT
2314
2315	END SUBROUTINE stokes_drift_terms_ij
2316
2317
2318	!------------------------------------------------------------------------------!
2319	! Description:
2320	! ------------
2321	!> This routine calculates turbulence generated by wave breaking near the ocean
2322	!> surface, following a parameterization given in Noh et al. (2004), JPO
2323	!> Call for all gridpoints.
2324	!> TODO: so far, this routine only works if the model time step has about the
2325	!> same value as the time scale of wave breaking!
2326	!------------------------------------------------------------------------------!
2327	SUBROUTINE wave_breaking_term( component )
2328
2329	USE arrays_3d, &
2330	ONLY: u_p, v_p
2331
2332	USE control_parameters, &
2333	ONLY: dt_3d, message_string
2334
2335	USE indices, &
2336	ONLY: nxl, nxlu, nxr, nys, nysv, nyn, nzt
2337
2338	IMPLICIT NONE
2339
2340	INTEGER(iwp) :: component !< component of momentum equation
2341	INTEGER(iwp) :: i !< loop index along x
2342	INTEGER(iwp) :: j !< loop index along y
2343
2344	REAL(wp) :: random_gauss !< function that creates a random number with a
2345	!< Gaussian distribution
2346
2347
2348	!
2349	!-- Compute wave breaking terms for the respective velocity components.
2350	!-- Velocities are directly manipulated, since this is not a real force
2351	SELECT CASE ( component )
2352
2353	!
2354	!-- u-component
2355	CASE ( 1 )
2356	DO i = nxlu, nxr
2357	DO j = nys, nyn
2358	u_p(nzt,j,i) = u_p(nzt,j,i) + &
2359	( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &
2360	* alpha_wave_breaking * u_star_wave_breaking &
2361	/ timescale_wave_breaking * dt_3d
2362	ENDDO
2363	ENDDO
2364	!
2365	!-- v-component
2366	CASE ( 2 )
2367	DO i = nxl, nxr
2368	DO j = nysv, nyn
2369	v_p(nzt,j,i) = v_p(nzt,j,i) + &
2370	( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &
2371	* alpha_wave_breaking * u_star_wave_breaking &
2372	/ timescale_wave_breaking * dt_3d
2373	ENDDO
2374	ENDDO
2375
2376	CASE DEFAULT
2377	WRITE( message_string, * ) 'wrong component of wave breaking: ', &
2378	component
2379	CALL message( 'stokes_drift_terms', 'PA0466', 1, 2, 0, 6, 0 )
2380
2381	END SELECT
2382
2383	END SUBROUTINE wave_breaking_term
2384
2385
2386	!------------------------------------------------------------------------------!
2387	! Description:
2388	! ------------
2389	!> This routine calculates turbulence generated by wave breaking near the ocean
2390	!> surface, following a parameterization given in Noh et al. (2004), JPO
2391	!> Call for gridpoint i,j.
2392	!> TODO: so far, this routine only works if the model time step has about the
2393	!> same value as the time scale of wave breaking!
2394	!------------------------------------------------------------------------------!
2395	SUBROUTINE wave_breaking_term_ij( i, j, component )
2396
2397	USE arrays_3d, &
2398	ONLY: u_p, v_p
2399
2400	USE control_parameters, &
2401	ONLY: dt_3d, message_string
2402
2403	USE indices, &
2404	ONLY: nzt
2405
2406	IMPLICIT NONE
2407
2408	INTEGER(iwp) :: component !< component of momentum equation
2409	INTEGER(iwp) :: i !< loop index along x
2410	INTEGER(iwp) :: j !< loop index along y
2411
2412	REAL(wp) :: random_gauss !< function that creates a random number with a
2413	!< Gaussian distribution
2414
2415	!
2416	!-- Compute wave breaking terms for the respective velocity components
2417	SELECT CASE ( component )
2418
2419	!
2420	!-- u-/v-component
2421	CASE ( 1 )
2422	u_p(nzt,j,i) = u_p(nzt,j,i) + &
2423	( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &
2424	* alpha_wave_breaking * u_star_wave_breaking &
2425	/ timescale_wave_breaking * dt_3d
2426
2427	CASE ( 2 )
2428	v_p(nzt,j,i) = v_p(nzt,j,i) + &
2429	( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &
2430	* alpha_wave_breaking * u_star_wave_breaking &
2431	/ timescale_wave_breaking * dt_3d
2432
2433	CASE DEFAULT
2434	WRITE( message_string, * ) 'wrong component of wave breaking: ', &
2435	component
2436	CALL message( 'stokes_drift_terms', 'PA0466', 1, 2, 0, 6, 0 )
2437
2438	END SELECT
2439
2440	END SUBROUTINE wave_breaking_term_ij
2441
2442
2443	END MODULE ocean_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |