Home

Context Navigation

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

Last change on this file since 3574 was 3568, checked in by raasch, 6 years ago
bugifx: calculate equation of state for seawater even if salinity is switched off
Property svn:keywords set to `Id`
File size: 82.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |