Home

Context Navigation

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

Last change on this file since 4734 was 4731, checked in by schwenkel, 4 years ago
Move exchange_horiz from time_integration to modules
Property svn:keywords set to `Id`
File size: 94.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |