Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |