Home

Context Navigation

ocean_mod.f90 @ 4867

Last change on this file since 4867 was 4843, checked in by raasch, 4 years ago
local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates
Property svn:keywords set to `Id`
File size: 101.6 KB

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

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

Context Navigation

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

Download in other formats: