Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |