Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |