Home

Context Navigation

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

Last change on this file since 3738 was 3719, checked in by kanani, 6 years ago
Correct and clean-up cpu_logs, some overlapping counts (chemistry_model_mod, disturb_heatflux, large_scale_forcing_nudging_mod, ocean_mod, palm, prognostic_equations, synthetic_turbulence_generator_mod, time_integration, time_integration_spinup, turbulence_closure_mod)
Property svn:keywords set to `Id`
File size: 82.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |