Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |