Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |