Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |