Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |