Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |