Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |