Home

Context Navigation

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

Last change on this file since 3790 was 3767, checked in by raasch, 6 years ago
unused variables removed from rrd-subroutines parameter list
Property svn:keywords set to `Id`
File size: 82.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |