Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |