Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |