Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |