Home

Context Navigation

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

Last change on this file since 3648 was 3636, checked in by raasch, 6 years ago
nopointer option removed
Property svn:keywords set to `Id`
File size: 82.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |