Home

Context Navigation

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

Last change on this file since 4763 was 4731, checked in by schwenkel, 4 years ago
Move exchange_horiz from time_integration to modules
Property svn:keywords set to `Id`
File size: 94.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |