Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |