Home

Context Navigation

ocean_mod.f90 @ 4867

Last change on this file since 4867 was 4843, checked in by raasch, 4 years ago
local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates
Property svn:keywords set to `Id`
File size: 101.6 KB

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

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

Context Navigation

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

Download in other formats: