Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |