Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |