Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |