Home

Context Navigation

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

Last change on this file since 3888 was 3873, checked in by knoop, 6 years ago
Moved ocean_mode specific code from advec_ws to ocean_mod + implemented ocean_actions
Property svn:keywords set to `Id`
File size: 84.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |