Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |