Home

Context Navigation

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

Last change on this file since 3382 was 3381, checked in by raasch, 6 years ago
spin-up cooling for ocean surface implemented, see new parameter surface_cooling_spinup_time
Property svn:keywords set to `Id`
File size: 82.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |