Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |