Home

Context Navigation

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

Last change on this file since 4245 was 4196, checked in by gronemeier, 5 years ago

Consider rotation of model domain for claculation of Coriolis force:

check_parameters: Overwrite rotation_angle from namelist by value from static driver;
coriolis: Consider rotation of model domain;
header: Write information about rotation angle;
modules: Added rotation_angle;
netcdf_interface_mod: replaced rotation angle from input-netCDF file by namelist parameter 'rotation_angle';
ocean_mod: Consider rotation of model domain for calculating the Stokes drift;
parin: added rotation_angle to initialization_parameters;

Property svn:keywords set to Id

File size: 87.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |