Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |