Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |