Home

Context Navigation

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

Last change on this file since 4499 was 4495, checked in by raasch, 5 years ago
restart data handling with MPI-IO added, first part
Property svn:keywords set to `Id`
File size: 92.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |