Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |