Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |