Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |