Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |