Home

Context Navigation

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

Last change on this file since 4235 was 4196, checked in by gronemeier, 5 years ago

Consider rotation of model domain for claculation of Coriolis force:

check_parameters: Overwrite rotation_angle from namelist by value from static driver;
coriolis: Consider rotation of model domain;
header: Write information about rotation angle;
modules: Added rotation_angle;
netcdf_interface_mod: replaced rotation angle from input-netCDF file by namelist parameter 'rotation_angle';
ocean_mod: Consider rotation of model domain for calculating the Stokes drift;
parin: added rotation_angle to initialization_parameters;

Property svn:keywords set to Id

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |