Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |