Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |