Home

Context Navigation

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

Last change on this file since 4810 was 4797, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 99.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |