Home

Context Navigation

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

Last change on this file since 4427 was 4370, checked in by raasch, 5 years ago
bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added
Property svn:keywords set to `Id`
File size: 89.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |