Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |