Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |