Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |