Home

Context Navigation

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

Last change on this file since 4791 was 4768, checked in by suehring, 4 years ago
Enable 3D data output also with 64-bit precision
Property svn:keywords set to `Id`
File size: 94.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |