Home

Context Navigation

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

Last change on this file since 4529 was 4517, checked in by raasch, 5 years ago
added restart with MPI-IO for reading local arrays
Property svn:keywords set to `Id`
File size: 93.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |