Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |