Home

Context Navigation

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

Last change on this file since 3301 was 3294, checked in by raasch, 6 years ago
modularization of the ocean code
Property svn:keywords set to `Id`
File size: 60.4 KB

Line
1	!> @file ocean_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 2017-2018 Leibniz Universitaet Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: ocean_mod.f90 3294 2018-10-01 02:37:10Z gronemeier $
27	! initial revision
28	!
29	!
30	!
31	!
32	! Authors:
33	! --------
34	! @author Siegfried Raasch
35	!
36	! Description:
37	! ------------
38	!> This module contains relevant code for PALM's ocean option, e.g. the
39	!> prognostic equation for salinity, the equation of state for seawater,
40	!> and the Craik Leibovich force (Stokes force)
41	!>
42	!------------------------------------------------------------------------------!
43	MODULE ocean_mod
44
45
46	#if defined( __nopointer )
47	USE arrays_3d, &
48	ONLY: prho, rho_ocean, sa, sa_init, sa_p, tsa_m
49	#else
50	USE arrays_3d, &
51	ONLY: prho, prho_1, rho_ocean, rho_1, sa, sa_init, sa_1, sa_2, sa_3, &
52	sa_p, tsa_m
53	#endif
54
55	USE control_parameters, &
56	ONLY: atmos_ocean_sign, bottom_salinityflux, &
57	constant_top_salinityflux, ocean_mode, top_salinityflux, &
58	wall_salinityflux
59
60	USE kinds
61
62
63	IMPLICIT NONE
64
65	CHARACTER (LEN=20) :: bc_sa_t = 'neumann' !< namelist parameter
66
67	INTEGER(iwp) :: ibc_sa_t !< integer flag for bc_sa_t
68
69	INTEGER(iwp) :: sa_vertical_gradient_level_ind(10) = -9999 !< grid index values of sa_vertical_gradient_level(s)
70
71	LOGICAL :: stokes_force = .FALSE. !< switch to switch on the Stokes force
72
73	REAL(wp) :: langmuir_number = 0.34_wp !< turbulent Langmuir number, default from Li et al., 2005
74	REAL(wp) :: prho_reference !< reference state of potential density at ocean surface
75	REAL(wp) :: sa_surface = 35.0_wp !< surface salinity, namelist parameter
76	REAL(wp) :: sa_vertical_gradient(10) = 0.0_wp !< namelist parameter
77	REAL(wp) :: sa_vertical_gradient_level(10) = -999999.9_wp !< namelist parameter
78	REAL(wp) :: stokes_waveheight = 1.0_wp !< typical wave height assumed for Stokes velocity
79	REAL(wp) :: stokes_wavelength = 40.0_wp !< typical wavelength assumed for Stokes velocity
80
81	REAL(wp), DIMENSION(12), PARAMETER :: nom = &
82	(/ 9.9984085444849347D2, 7.3471625860981584D0, &
83	-5.3211231792841769D-2, 3.6492439109814549D-4, &
84	2.5880571023991390D0, -6.7168282786692354D-3, &
85	1.9203202055760151D-3, 1.1798263740430364D-2, &
86	9.8920219266399117D-8, 4.6996642771754730D-6, &
87	-2.5862187075154352D-8, -3.2921414007960662D-12 /)
88	!< constants used in equation of state for seawater
89
90	REAL(wp), DIMENSION(13), PARAMETER :: den = &
91	(/ 1.0D0, 7.2815210113327091D-3, &
92	-4.4787265461983921D-5, 3.3851002965802430D-7, &
93	1.3651202389758572D-10, 1.7632126669040377D-3, &
94	-8.8066583251206474D-6, -1.8832689434804897D-10, &
95	5.7463776745432097D-6, 1.4716275472242334D-9, &
96	6.7103246285651894D-6, -2.4461698007024582D-17, &
97	-9.1534417604289062D-18 /)
98	!< constants used in equation of state for seawater
99
100	SAVE
101
102	PUBLIC :: bc_sa_t, ibc_sa_t, langmuir_number, prho_reference, sa_surface, &
103	sa_vertical_gradient, sa_vertical_gradient_level, &
104	sa_vertical_gradient_level_ind, stokes_force, &
105	stokes_waveheight, stokes_wavelength
106
107
108	INTERFACE eqn_state_seawater
109	MODULE PROCEDURE eqn_state_seawater
110	MODULE PROCEDURE eqn_state_seawater_ij
111	END INTERFACE eqn_state_seawater
112
113	INTERFACE eqn_state_seawater_func
114	MODULE PROCEDURE eqn_state_seawater_func
115	END INTERFACE eqn_state_seawater_func
116
117	INTERFACE ocean_check_parameters
118	MODULE PROCEDURE ocean_check_parameters
119	END INTERFACE ocean_check_parameters
120
121	INTERFACE ocean_check_data_output
122	MODULE PROCEDURE ocean_check_data_output
123	END INTERFACE ocean_check_data_output
124
125	INTERFACE ocean_check_data_output_pr
126	MODULE PROCEDURE ocean_check_data_output_pr
127	END INTERFACE ocean_check_data_output_pr
128
129	INTERFACE ocean_define_netcdf_grid
130	MODULE PROCEDURE ocean_define_netcdf_grid
131	END INTERFACE ocean_define_netcdf_grid
132
133	INTERFACE ocean_data_output_2d
134	MODULE PROCEDURE ocean_data_output_2d
135	END INTERFACE ocean_data_output_2d
136
137	INTERFACE ocean_data_output_3d
138	MODULE PROCEDURE ocean_data_output_3d
139	END INTERFACE ocean_data_output_3d
140
141	INTERFACE ocean_init
142	MODULE PROCEDURE ocean_init
143	END INTERFACE ocean_init
144
145	INTERFACE ocean_init_arrays
146	MODULE PROCEDURE ocean_init_arrays
147	END INTERFACE ocean_init_arrays
148
149	INTERFACE ocean_parin
150	MODULE PROCEDURE ocean_parin
151	END INTERFACE ocean_parin
152
153	INTERFACE ocean_prognostic_equations
154	MODULE PROCEDURE ocean_prognostic_equations
155	MODULE PROCEDURE ocean_prognostic_equations_ij
156	END INTERFACE ocean_prognostic_equations
157
158	INTERFACE ocean_swap_timelevel
159	MODULE PROCEDURE ocean_swap_timelevel
160	END INTERFACE ocean_swap_timelevel
161
162	INTERFACE ocean_rrd_global
163	MODULE PROCEDURE ocean_rrd_global
164	END INTERFACE ocean_rrd_global
165
166	INTERFACE ocean_rrd_local
167	MODULE PROCEDURE ocean_rrd_local
168	END INTERFACE ocean_rrd_local
169
170	INTERFACE ocean_wrd_global
171	MODULE PROCEDURE ocean_wrd_global
172	END INTERFACE ocean_wrd_global
173
174	INTERFACE ocean_wrd_local
175	MODULE PROCEDURE ocean_wrd_local
176	END INTERFACE ocean_wrd_local
177
178	INTERFACE ocean_3d_data_averaging
179	MODULE PROCEDURE ocean_3d_data_averaging
180	END INTERFACE ocean_3d_data_averaging
181
182	SAVE
183
184	PRIVATE
185	!
186	!-- Add INTERFACES that must be available to other modules (alphabetical order)
187	PUBLIC eqn_state_seawater, ocean_check_data_output, &
188	ocean_check_data_output_pr, ocean_check_parameters, &
189	ocean_data_output_2d, ocean_data_output_3d, &
190	ocean_define_netcdf_grid, ocean_init, ocean_init_arrays, &
191	ocean_parin, ocean_prognostic_equations, ocean_swap_timelevel, &
192	ocean_rrd_global, ocean_rrd_local, ocean_wrd_global, &
193	ocean_wrd_local, ocean_3d_data_averaging
194
195
196	CONTAINS
197
198	!------------------------------------------------------------------------------!
199	! Description:
200	! ------------
201	!> Equation of state for seawater as a function of potential temperature,
202	!> salinity, and pressure.
203	!> For coefficients see Jackett et al., 2006: J. Atm. Ocean Tech.
204	!> eqn_state_seawater calculates the potential density referred at hyp(0).
205	!> eqn_state_seawater_func calculates density.
206	!> TODO: so far, routine is not adjusted to use topography
207	!------------------------------------------------------------------------------!
208	SUBROUTINE eqn_state_seawater
209
210	USE arrays_3d, &
211	ONLY: hyp, prho, pt_p, rho_ocean, sa_p
212	USE indices, &
213	ONLY: nxl, nxr, nyn, nys, nzb, nzt
214
215	USE surface_mod, &
216	ONLY : bc_h
217
218	IMPLICIT NONE
219
220	INTEGER(iwp) :: i !< running index x direction
221	INTEGER(iwp) :: j !< running index y direction
222	INTEGER(iwp) :: k !< running index z direction
223	INTEGER(iwp) :: m !< running index surface elements
224
225	REAL(wp) :: pden !< temporary scalar user for calculating density
226	REAL(wp) :: pnom !< temporary scalar user for calculating density
227	REAL(wp) :: p1 !< temporary scalar user for calculating density
228	REAL(wp) :: p2 !< temporary scalar user for calculating density
229	REAL(wp) :: p3 !< temporary scalar user for calculating density
230	REAL(wp) :: pt1 !< temporary scalar user for calculating density
231	REAL(wp) :: pt2 !< temporary scalar user for calculating density
232	REAL(wp) :: pt3 !< temporary scalar user for calculating density
233	REAL(wp) :: pt4 !< temporary scalar user for calculating density
234	REAL(wp) :: sa1 !< temporary scalar user for calculating density
235	REAL(wp) :: sa15 !< temporary scalar user for calculating density
236	REAL(wp) :: sa2 !< temporary scalar user for calculating density
237
238
239	DO i = nxl, nxr
240	DO j = nys, nyn
241	DO k = nzb+1, nzt
242	!
243	!-- Pressure is needed in dbar
244	p1 = hyp(k) * 1E-4_wp
245	p2 = p1 * p1
246	p3 = p2 * p1
247
248	!
249	!-- Temperature needed in degree Celsius
250	pt1 = pt_p(k,j,i) - 273.15_wp
251	pt2 = pt1 * pt1
252	pt3 = pt1 * pt2
253	pt4 = pt2 * pt2
254
255	sa1 = sa_p(k,j,i)
256	sa15 = sa1 * SQRT( sa1 )
257	sa2 = sa1 * sa1
258
259	pnom = nom(1) + nom(2)pt1 + nom(3)pt2 + &
260	nom(4)pt3 + nom(5)sa1 + nom(6)sa1pt1 + &
261	nom(7)*sa2
262
263	pden = den(1) + den(2)pt1 + den(3)pt2 + &
264	den(4)pt3 + den(5)pt4 + den(6)*sa1 + &
265	den(7)sa1pt1 + den(8)sa1pt3 + den(9)*sa15 + &
266	den(10)sa15pt2
267	!
268	!-- Potential density (without pressure terms)
269	prho(k,j,i) = pnom / pden
270
271	pnom = pnom + nom(8)p1 + nom(9)p1*pt2 + &
272	nom(10)p1sa1 + nom(11)p2 + nom(12)p2*pt2
273
274	pden = pden + den(11)p1 + den(12)p2*pt3 + &
275	den(13)p3pt1
276
277	!
278	!-- In-situ density
279	rho_ocean(k,j,i) = pnom / pden
280
281	ENDDO
282
283	!
284	!-- Neumann conditions are assumed at top boundary and currently also at
285	!-- bottom boundary (see comment lines below)
286	prho(nzt+1,j,i) = prho(nzt,j,i)
287	rho_ocean(nzt+1,j,i) = rho_ocean(nzt,j,i)
288
289	ENDDO
290	ENDDO
291	!
292	!-- Neumann conditions at up/downward-facing surfaces
293	!$OMP PARALLEL DO PRIVATE( i, j, k )
294	DO m = 1, bc_h(0)%ns
295	i = bc_h(0)%i(m)
296	j = bc_h(0)%j(m)
297	k = bc_h(0)%k(m)
298	prho(k-1,j,i) = prho(k,j,i)
299	rho_ocean(k-1,j,i) = rho_ocean(k,j,i)
300	ENDDO
301	!
302	!-- Downward facing surfaces
303	!$OMP PARALLEL DO PRIVATE( i, j, k )
304	DO m = 1, bc_h(1)%ns
305	i = bc_h(1)%i(m)
306	j = bc_h(1)%j(m)
307	k = bc_h(1)%k(m)
308	prho(k+1,j,i) = prho(k,j,i)
309	rho_ocean(k+1,j,i) = rho_ocean(k,j,i)
310	ENDDO
311
312	END SUBROUTINE eqn_state_seawater
313
314
315	!------------------------------------------------------------------------------!
316	! Description:
317	! ------------
318	!> Same as above, but for grid point i,j
319	!------------------------------------------------------------------------------!
320	SUBROUTINE eqn_state_seawater_ij( i, j )
321
322	USE arrays_3d, &
323	ONLY: hyp, prho, pt_p, rho_ocean, sa_p
324
325	USE indices, &
326	ONLY: nzb, nzt
327
328	USE surface_mod, &
329	ONLY : bc_h
330
331	IMPLICIT NONE
332
333	INTEGER(iwp) :: i !< running index x direction
334	INTEGER(iwp) :: j !< running index y direction
335	INTEGER(iwp) :: k !< running index z direction
336	INTEGER(iwp) :: m !< running index surface elements
337	INTEGER(iwp) :: surf_e !< end index of surface elements at (j,i)-gridpoint
338	INTEGER(iwp) :: surf_s !< start index of surface elements at (j,i)-gridpoint
339
340	REAL(wp) :: pden !< temporary scalar user for calculating density
341	REAL(wp) :: pnom !< temporary scalar user for calculating density
342	REAL(wp) :: p1 !< temporary scalar user for calculating density
343	REAL(wp) :: p2 !< temporary scalar user for calculating density
344	REAL(wp) :: p3 !< temporary scalar user for calculating density
345	REAL(wp) :: pt1 !< temporary scalar user for calculating density
346	REAL(wp) :: pt2 !< temporary scalar user for calculating density
347	REAL(wp) :: pt3 !< temporary scalar user for calculating density
348	REAL(wp) :: pt4 !< temporary scalar user for calculating density
349	REAL(wp) :: sa1 !< temporary scalar user for calculating density
350	REAL(wp) :: sa15 !< temporary scalar user for calculating density
351	REAL(wp) :: sa2 !< temporary scalar user for calculating density
352
353	DO k = nzb+1, nzt
354	!
355	!-- Pressure is needed in dbar
356	p1 = hyp(k) * 1E-4_wp
357	p2 = p1 * p1
358	p3 = p2 * p1
359	!
360	!-- Temperature needed in degree Celsius
361	pt1 = pt_p(k,j,i) - 273.15_wp
362	pt2 = pt1 * pt1
363	pt3 = pt1 * pt2
364	pt4 = pt2 * pt2
365
366	sa1 = sa_p(k,j,i)
367	sa15 = sa1 * SQRT( sa1 )
368	sa2 = sa1 * sa1
369
370	pnom = nom(1) + nom(2)pt1 + nom(3)pt2 + &
371	nom(4)pt3 + nom(5)sa1 + nom(6)sa1pt1 + nom(7)*sa2
372
373	pden = den(1) + den(2)pt1 + den(3)pt2 + &
374	den(4)pt3 + den(5)pt4 + den(6)*sa1 + &
375	den(7)sa1pt1 + den(8)sa1pt3 + den(9)*sa15 + &
376	den(10)sa15pt2
377	!
378	!-- Potential density (without pressure terms)
379	prho(k,j,i) = pnom / pden
380
381	pnom = pnom + nom(8)p1 + nom(9)p1*pt2 + &
382	nom(10)p1sa1 + nom(11)p2 + nom(12)p2*pt2
383	pden = pden + den(11)p1 + den(12)p2*pt3 + &
384	den(13)p3pt1
385
386	!
387	!-- In-situ density
388	rho_ocean(k,j,i) = pnom / pden
389
390	ENDDO
391	!
392	!-- Neumann conditions at up/downward-facing walls
393	surf_s = bc_h(0)%start_index(j,i)
394	surf_e = bc_h(0)%end_index(j,i)
395	DO m = surf_s, surf_e
396	k = bc_h(0)%k(m)
397	prho(k-1,j,i) = prho(k,j,i)
398	rho_ocean(k-1,j,i) = rho_ocean(k,j,i)
399	ENDDO
400	!
401	!-- Downward facing surfaces
402	surf_s = bc_h(1)%start_index(j,i)
403	surf_e = bc_h(1)%end_index(j,i)
404	DO m = surf_s, surf_e
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	!-- Neumann condition are assumed at top boundary
411	prho(nzt+1,j,i) = prho(nzt,j,i)
412	rho_ocean(nzt+1,j,i) = rho_ocean(nzt,j,i)
413
414	END SUBROUTINE eqn_state_seawater_ij
415
416
417	!------------------------------------------------------------------------------!
418	! Description:
419	! ------------
420	!> Equation of state for seawater as a function
421	!------------------------------------------------------------------------------!
422	REAL(wp) FUNCTION eqn_state_seawater_func( p, pt, sa )
423
424	IMPLICIT NONE
425
426	REAL(wp) :: p !< temporary scalar user for calculating density
427	REAL(wp) :: p1 !< temporary scalar user for calculating density
428	REAL(wp) :: p2 !< temporary scalar user for calculating density
429	REAL(wp) :: p3 !< temporary scalar user for calculating density
430	REAL(wp) :: pt !< temporary scalar user for calculating density
431	REAL(wp) :: pt1 !< temporary scalar user for calculating density
432	REAL(wp) :: pt2 !< temporary scalar user for calculating density
433	REAL(wp) :: pt3 !< temporary scalar user for calculating density
434	REAL(wp) :: pt4 !< temporary scalar user for calculating density
435	REAL(wp) :: sa !< temporary scalar user for calculating density
436	REAL(wp) :: sa15 !< temporary scalar user for calculating density
437	REAL(wp) :: sa2 !< temporary scalar user for calculating density
438
439	!
440	!-- Pressure is needed in dbar
441	p1 = p * 1.0E-4_wp
442	p2 = p1 * p1
443	p3 = p2 * p1
444
445	!
446	!-- Temperature needed in degree Celsius
447	pt1 = pt - 273.15_wp
448	pt2 = pt1 * pt1
449	pt3 = pt1 * pt2
450	pt4 = pt2 * pt2
451
452	sa15 = sa * SQRT( sa )
453	sa2 = sa * sa
454
455
456	eqn_state_seawater_func = &
457	( nom(1) + nom(2)pt1 + nom(3)pt2 + nom(4)*pt3 + &
458	nom(5)sa + nom(6)sapt1 + nom(7)sa2 + nom(8)*p1 + &
459	nom(9)p1pt2 + nom(10)p1sa + nom(11)p2 + nom(12)p2*pt2 &
460	) / &
461	( den(1) + den(2)pt1 + den(3)pt2 + den(4)*pt3 + &
462	den(5)pt4 + den(6)sa + den(7)sapt1 + den(8)sapt3 + &
463	den(9)sa15 + den(10)sa15pt2 + den(11)p1 + den(12)p2pt3 + &
464	den(13)p3pt1 &
465	)
466
467
468	END FUNCTION eqn_state_seawater_func
469
470
471	!------------------------------------------------------------------------------!
472	! Description:
473	! ------------
474	!> Reads the ocean parameters namelist
475	!------------------------------------------------------------------------------!
476	SUBROUTINE ocean_parin
477
478	USE control_parameters, &
479	ONLY: message_string
480
481	IMPLICIT NONE
482
483	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
484
485
486	NAMELIST /ocean_parameters/ bc_sa_t, bottom_salinityflux, &
487	langmuir_number, sa_surface, sa_vertical_gradient, &
488	sa_vertical_gradient_level, stokes_force, stokes_waveheight, &
489	stokes_wavelength, top_salinityflux, wall_salinityflux
490
491	!
492	!-- Try to find the namelist
493	REWIND ( 11 )
494	line = ' '
495	DO WHILE ( INDEX( line, '&ocean_parameters' ) == 0 )
496	READ ( 11, '(A)', END=12 ) line
497	ENDDO
498	BACKSPACE ( 11 )
499
500	!
501	!-- Read namelist
502	READ ( 11, ocean_parameters, ERR = 10 )
503	!
504	!-- Set switch that enables PALM's ocean mode
505	ocean_mode = .TRUE.
506
507	GOTO 12
508
509	10 BACKSPACE( 11 )
510	READ( 11 , '(A)') line
511	CALL parin_fail_message( 'ocean_parameters', line )
512
513	12 CONTINUE
514
515	END SUBROUTINE ocean_parin
516
517	!------------------------------------------------------------------------------!
518	! Description:
519	! ------------
520	!> Check parameters routine for the ocean mode
521	!------------------------------------------------------------------------------!
522	SUBROUTINE ocean_check_parameters
523
524	USE control_parameters, &
525	ONLY: coupling_char, coupling_mode, message_string, use_top_fluxes
526
527	IMPLICIT NONE
528
529
530	!
531	!-- Check ocean setting
532	IF ( TRIM( coupling_mode ) == 'uncoupled' .AND. &
533	TRIM( coupling_char ) == '_O' .AND. &
534	.NOT. ocean_mode ) THEN
535
536	!
537	!-- Check whether an (uncoupled) atmospheric run has been declared as an
538	!-- ocean run (this setting is done via palmrun-option -y)
539	message_string = 'ocean = .F. does not allow coupling_char = "' // &
540	TRIM( coupling_char ) // '" set by palmrun-option "-y"'
541	CALL message( 'check_parameters', 'PA0317', 1, 2, 0, 6, 0 )
542
543	ENDIF
544
545	!
546	!-- Ocean version must use flux boundary conditions at the top
547	IF ( .NOT. use_top_fluxes ) THEN
548	message_string = 'use_top_fluxes must be .TRUE. in ocean mode'
549	CALL message( 'ocean_check_parameters', 'PA0042', 1, 2, 0, 6, 0 )
550	ENDIF
551
552	!
553	!-- Boundary conditions for salinity
554	IF ( bc_sa_t == 'dirichlet' ) THEN
555	ibc_sa_t = 0
556	ELSEIF ( bc_sa_t == 'neumann' ) THEN
557	ibc_sa_t = 1
558	ELSE
559	message_string = 'unknown boundary condition: bc_sa_t = "' // &
560	TRIM( bc_sa_t ) // '"'
561	CALL message( 'ocean_check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
562	ENDIF
563
564	IF ( top_salinityflux == 9999999.9_wp ) constant_top_salinityflux = .FALSE.
565
566	IF ( ibc_sa_t == 1 .AND. top_salinityflux == 9999999.9_wp ) THEN
567	message_string = 'boundary condition: bc_sa_t = "' // &
568	TRIM( bc_sa_t ) // '" requires to set top_salinityflux'
569	CALL message( 'ocean_check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
570	ENDIF
571
572	!
573	!-- A fixed salinity at the top implies Dirichlet boundary condition for
574	!-- salinity. In this case specification of a constant salinity flux is
575	!-- forbidden.
576	IF ( ibc_sa_t == 0 .AND. constant_top_salinityflux .AND. &
577	top_salinityflux /= 0.0_wp ) THEN
578	message_string = 'boundary condition: bc_sa_t = "' // &
579	TRIM( bc_sa_t ) // '" is not allowed with ' // &
580	'top_salinityflux /= 0.0'
581	CALL message( 'ocean_check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
582	ENDIF
583
584	END SUBROUTINE ocean_check_parameters
585
586
587	!------------------------------------------------------------------------------!
588	! Description:
589	! ------------
590	!> Check data output.
591	!------------------------------------------------------------------------------!
592	SUBROUTINE ocean_check_data_output( var, unit )
593
594	IMPLICIT NONE
595
596	CHARACTER (LEN=*) :: unit !< unit of output variable
597	CHARACTER (LEN=*) :: var !< name of output variable
598
599
600	SELECT CASE ( TRIM( var ) )
601
602	CASE ( 'rho_ocean' )
603	unit = 'kg/m3'
604
605	CASE ( 'sa' )
606	unit = 'psu'
607
608	CASE DEFAULT
609	unit = 'illegal'
610
611	END SELECT
612
613	END SUBROUTINE ocean_check_data_output
614
615
616	!------------------------------------------------------------------------------!
617	! Description:
618	! ------------
619	!> Check data output of profiles
620	!------------------------------------------------------------------------------!
621	SUBROUTINE ocean_check_data_output_pr( variable, var_count, unit, dopr_unit )
622
623	USE arrays_3d, &
624	ONLY: zu, zw
625
626	USE control_parameters, &
627	ONLY: data_output_pr, message_string
628
629	USE indices
630
631	USE profil_parameter
632
633	USE statistics
634
635	IMPLICIT NONE
636
637	CHARACTER (LEN=*) :: unit !<
638	CHARACTER (LEN=*) :: variable !<
639	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
640
641	INTEGER(iwp) :: var_count !<
642
643	SELECT CASE ( TRIM( variable ) )
644
645	CASE ( 'prho' )
646	dopr_index(var_count) = 71
647	dopr_unit = 'kg/m3'
648	hom(:,2,71,:) = SPREAD( zu, 2, statistic_regions+1 )
649	unit = dopr_unit
650
651	CASE ( 'rho_ocean' )
652	dopr_index(var_count) = 64
653	dopr_unit = 'kg/m3'
654	hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
655	IF ( data_output_pr(var_count)(1:1) == '#' ) THEN
656	dopr_initial_index(var_count) = 77
657	hom(:,2,77,:) = SPREAD( zu, 2, statistic_regions+1 )
658	hom(nzb,2,77,:) = 0.0_wp ! because zu(nzb) is negative
659	data_output_pr(var_count) = data_output_pr(var_count)(2:)
660	ENDIF
661	unit = dopr_unit
662
663	CASE ( 'sa', '#sa' )
664	dopr_index(var_count) = 23
665	dopr_unit = 'psu'
666	hom(:,2,23,:) = SPREAD( zu, 2, statistic_regions+1 )
667	IF ( data_output_pr(var_count)(1:1) == '#' ) THEN
668	dopr_initial_index(var_count) = 26
669	hom(:,2,26,:) = SPREAD( zu, 2, statistic_regions+1 )
670	hom(nzb,2,26,:) = 0.0_wp ! because zu(nzb) is negative
671	data_output_pr(var_count) = data_output_pr(var_count)(2:)
672	ENDIF
673	unit = dopr_unit
674
675	CASE ( 'w"sa"' )
676	dopr_index(var_count) = 65
677	dopr_unit = 'psu m/s'
678	hom(:,2,65,:) = SPREAD( zw, 2, statistic_regions+1 )
679	unit = dopr_unit
680
681	CASE ( 'wsa' )
682	dopr_index(var_count) = 66
683	dopr_unit = 'psu m/s'
684	hom(:,2,66,:) = SPREAD( zw, 2, statistic_regions+1 )
685	unit = dopr_unit
686
687	CASE ( 'wsa' )
688	dopr_index(var_count) = 67
689	dopr_unit = 'psu m/s'
690	hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
691	unit = dopr_unit
692
693	CASE DEFAULT
694	unit = 'illegal'
695
696	END SELECT
697
698
699	END SUBROUTINE ocean_check_data_output_pr
700
701
702	!------------------------------------------------------------------------------!
703	! Description:
704	! ------------
705	!> Define appropriate grid for netcdf variables.
706	!> It is called out from subroutine netcdf.
707	!------------------------------------------------------------------------------!
708	SUBROUTINE ocean_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
709
710	IMPLICIT NONE
711
712	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< x grid of output variable
713	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< y grid of output variable
714	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< z grid of output variable
715	CHARACTER (LEN=*), INTENT(IN) :: var !< name of output variable
716
717	LOGICAL, INTENT(OUT) :: found !< flag if output variable is found
718
719	found = .TRUE.
720
721	!
722	!-- Check for the grid
723	SELECT CASE ( TRIM( var ) )
724
725	CASE ( 'rho_ocean', 'sa' )
726	grid_x = 'x'
727	grid_y = 'y'
728	grid_z = 'zu'
729
730	CASE DEFAULT
731	found = .FALSE.
732	grid_x = 'none'
733	grid_y = 'none'
734	grid_z = 'none'
735
736	END SELECT
737
738	END SUBROUTINE ocean_define_netcdf_grid
739
740
741	!------------------------------------------------------------------------------!
742	! Description:
743	! ------------
744	!> Average 3D data.
745	!------------------------------------------------------------------------------!
746	SUBROUTINE ocean_3d_data_averaging( mode, variable )
747
748
749	USE arrays_3d, &
750	ONLY: rho_ocean, sa
751
752	USE averaging, &
753	ONLY: rho_ocean_av, sa_av
754
755	USE control_parameters, &
756	ONLY: average_count_3d
757
758	USE indices, &
759	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt
760
761	IMPLICIT NONE
762
763	CHARACTER (LEN=*) :: mode !< flag defining mode 'allocate', 'sum' or 'average'
764	CHARACTER (LEN=*) :: variable !< name of variable
765
766	INTEGER(iwp) :: i !< loop index
767	INTEGER(iwp) :: j !< loop index
768	INTEGER(iwp) :: k !< loop index
769
770	IF ( mode == 'allocate' ) THEN
771
772	SELECT CASE ( TRIM( variable ) )
773
774	CASE ( 'rho_ocean' )
775	IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
776	ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
777	ENDIF
778	rho_ocean_av = 0.0_wp
779
780	CASE ( 'sa' )
781	IF ( .NOT. ALLOCATED( sa_av ) ) THEN
782	ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
783	ENDIF
784	sa_av = 0.0_wp
785
786	CASE DEFAULT
787	CONTINUE
788
789	END SELECT
790
791	ELSEIF ( mode == 'sum' ) THEN
792
793	SELECT CASE ( TRIM( variable ) )
794
795	CASE ( 'rho_ocean' )
796	IF ( ALLOCATED( rho_ocean_av ) ) THEN
797	DO i = nxlg, nxrg
798	DO j = nysg, nyng
799	DO k = nzb, nzt+1
800	rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) + &
801	rho_ocean(k,j,i)
802	ENDDO
803	ENDDO
804	ENDDO
805	ENDIF
806
807	CASE ( 'sa' )
808	IF ( ALLOCATED( sa_av ) ) THEN
809	DO i = nxlg, nxrg
810	DO j = nysg, nyng
811	DO k = nzb, nzt+1
812	sa_av(k,j,i) = sa_av(k,j,i) + sa(k,j,i)
813	ENDDO
814	ENDDO
815	ENDDO
816	ENDIF
817
818	CASE DEFAULT
819	CONTINUE
820
821	END SELECT
822
823	ELSEIF ( mode == 'average' ) THEN
824
825	SELECT CASE ( TRIM( variable ) )
826
827	CASE ( 'rho_ocean' )
828	IF ( ALLOCATED( rho_ocean_av ) ) THEN
829	DO i = nxlg, nxrg
830	DO j = nysg, nyng
831	DO k = nzb, nzt+1
832	rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) / &
833	REAL( average_count_3d, KIND=wp )
834	ENDDO
835	ENDDO
836	ENDDO
837	ENDIF
838
839	CASE ( 'sa' )
840	IF ( ALLOCATED( sa_av ) ) THEN
841	DO i = nxlg, nxrg
842	DO j = nysg, nyng
843	DO k = nzb, nzt+1
844	sa_av(k,j,i) = sa_av(k,j,i) / &
845	REAL( average_count_3d, KIND=wp )
846	ENDDO
847	ENDDO
848	ENDDO
849	ENDIF
850
851	END SELECT
852
853	ENDIF
854
855	END SUBROUTINE ocean_3d_data_averaging
856
857
858	!------------------------------------------------------------------------------!
859	! Description:
860	! ------------
861	!> Define 2D output variables.
862	!------------------------------------------------------------------------------!
863	SUBROUTINE ocean_data_output_2d( av, variable, found, grid, mode, local_pf, &
864	nzb_do, nzt_do )
865
866	USE arrays_3d, &
867	ONLY: rho_ocean, sa
868
869	USE averaging, &
870	ONLY: rho_ocean_av, sa_av
871
872	USE indices, &
873	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
874
875	IMPLICIT NONE
876
877	CHARACTER (LEN=*) :: grid !< name of vertical grid
878	CHARACTER (LEN=*) :: mode !< either 'xy', 'xz' or 'yz'
879	CHARACTER (LEN=*) :: variable !< name of variable
880
881	INTEGER(iwp) :: av !< flag for (non-)average output
882	INTEGER(iwp) :: flag_nr !< number of masking flag
883	INTEGER(iwp) :: i !< loop index
884	INTEGER(iwp) :: j !< loop index
885	INTEGER(iwp) :: k !< loop index
886	INTEGER(iwp) :: nzb_do !< vertical output index (bottom)
887	INTEGER(iwp) :: nzt_do !< vertical output index (top)
888
889	LOGICAL :: found !< flag if output variable is found
890	LOGICAL :: resorted !< flag if output is already resorted
891
892	REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute
893
894	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local
895	!< array to which output data is resorted to
896
897	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable
898
899	found = .TRUE.
900	resorted = .FALSE.
901	!
902	!-- Set masking flag for topography for not resorted arrays
903	flag_nr = 0
904
905	SELECT CASE ( TRIM( variable ) )
906
907	CASE ( 'rho_ocean_xy', 'rho_ocean_xz', 'rho_ocean_yz' )
908	IF ( av == 0 ) THEN
909	to_be_resorted => rho_ocean
910	ELSE
911	IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
912	ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
913	rho_ocean_av = REAL( fill_value, KIND = wp )
914	ENDIF
915	to_be_resorted => rho_ocean_av
916	ENDIF
917
918	CASE ( 'sa_xy', 'sa_xz', 'sa_yz' )
919	IF ( av == 0 ) THEN
920	to_be_resorted => sa
921	ELSE
922	IF ( .NOT. ALLOCATED( sa_av ) ) THEN
923	ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
924	sa_av = REAL( fill_value, KIND = wp )
925	ENDIF
926	to_be_resorted => sa_av
927	ENDIF
928	IF ( mode == 'xy' ) grid = 'zu'
929
930	CASE DEFAULT
931	found = .FALSE.
932	grid = 'none'
933
934	END SELECT
935
936	IF ( found .AND. .NOT. resorted ) THEN
937	! DO i = nxl, nxr
938	! DO j = nys, nyn
939	! DO k = nzb_do, nzt_do
940	! local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
941	! REAL( fill_value, KIND = wp ), &
942	! BTEST( wall_flags_0(k,j,i), flag_nr ) )
943	! ENDDO
944	! ENDDO
945	! ENDDO
946	ENDIF
947
948	END SUBROUTINE ocean_data_output_2d
949
950
951	!------------------------------------------------------------------------------!
952	! Description:
953	! ------------
954	!> Define 3D output variables.
955	!------------------------------------------------------------------------------!
956	SUBROUTINE ocean_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
957
958
959	USE arrays_3d, &
960	ONLY: rho_ocean, sa
961
962	USE averaging, &
963	ONLY: rho_ocean_av, sa_av
964
965	USE indices, &
966	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, &
967	wall_flags_0
968
969	IMPLICIT NONE
970
971	CHARACTER (LEN=*) :: variable !< name of variable
972
973	INTEGER(iwp) :: av !< flag for (non-)average output
974	INTEGER(iwp) :: flag_nr !< number of masking flag
975	INTEGER(iwp) :: i !< loop index
976	INTEGER(iwp) :: j !< loop index
977	INTEGER(iwp) :: k !< loop index
978	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
979	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
980
981	LOGICAL :: found !< flag if output variable is found
982	LOGICAL :: resorted !< flag if output is already resorted
983
984	REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute
985
986	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local
987	!< array to which output data is resorted to
988
989	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable
990
991	found = .TRUE.
992	resorted = .FALSE.
993	!
994	!-- Set masking flag for topography for not resorted arrays
995	flag_nr = 0
996
997	SELECT CASE ( TRIM( variable ) )
998
999	CASE ( 'rho_ocean' )
1000	IF ( av == 0 ) THEN
1001	to_be_resorted => rho_ocean
1002	ELSE
1003	IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
1004	ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1005	rho_ocean_av = REAL( fill_value, KIND = wp )
1006	ENDIF
1007	to_be_resorted => rho_ocean_av
1008	ENDIF
1009
1010	CASE ( 'sa' )
1011	IF ( av == 0 ) THEN
1012	to_be_resorted => sa
1013	ELSE
1014	IF ( .NOT. ALLOCATED( sa_av ) ) THEN
1015	ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1016	sa_av = REAL( fill_value, KIND = wp )
1017	ENDIF
1018	to_be_resorted => sa_av
1019	ENDIF
1020
1021	CASE DEFAULT
1022	found = .FALSE.
1023
1024	END SELECT
1025
1026
1027	IF ( found .AND. .NOT. resorted ) THEN
1028	DO i = nxl, nxr
1029	DO j = nys, nyn
1030	DO k = nzb_do, nzt_do
1031	local_pf(i,j,k) = MERGE( &
1032	to_be_resorted(k,j,i), &
1033	REAL( fill_value, KIND = wp ), &
1034	BTEST( wall_flags_0(k,j,i), flag_nr ) )
1035	ENDDO
1036	ENDDO
1037	ENDDO
1038	resorted = .TRUE.
1039	ENDIF
1040
1041	END SUBROUTINE ocean_data_output_3d
1042
1043
1044	!------------------------------------------------------------------------------!
1045	! Description:
1046	! ------------
1047	!> Allocate arrays and assign pointers.
1048	!------------------------------------------------------------------------------!
1049	SUBROUTINE ocean_init_arrays
1050
1051	USE indices, &
1052	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt
1053
1054	USE pmc_interface, &
1055	ONLY: nested_run
1056
1057	IMPLICIT NONE
1058
1059	#if defined( __nopointer )
1060	ALLOCATE( prho(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
1061	rho_ocean(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
1062	sa(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
1063	sa_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
1064	tsa_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1065	#else
1066	ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
1067	rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
1068	sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
1069	sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
1070	sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1071
1072	prho => prho_1
1073	rho_ocean => rho_1 ! routines calc_mean_profile and diffusion_e require
1074	! density to be a pointer
1075	#endif
1076
1077	#if ! defined( __nopointer )
1078	!
1079	!-- Initial assignment of pointers
1080	sa => sa_1; sa_p => sa_2; tsa_m => sa_3
1081	#endif
1082
1083	END SUBROUTINE ocean_init_arrays
1084
1085
1086	!------------------------------------------------------------------------------!
1087	! Description:
1088	! ------------
1089	!> Initialization of quantities needed for the ocean mode
1090	!------------------------------------------------------------------------------!
1091	SUBROUTINE ocean_init
1092
1093
1094	USE arrays_3d, &
1095	ONLY: dzu, hyp, pt_init, ref_state, zu, zw
1096
1097	USE basic_constants_and_equations_mod, &
1098	ONLY: g
1099
1100	USE control_parameters, &
1101	ONLY: initializing_actions, molecular_viscosity, rho_surface, &
1102	rho_reference, surface_pressure, use_single_reference_value
1103
1104	USE indices, &
1105	ONLY: nxl, nxlg, nxrg, nyng, nys, nysg, nzb, nzt
1106
1107	USE kinds
1108
1109	USE pegrid
1110
1111	USE statistics, &
1112	ONLY: hom, statistic_regions
1113
1114	IMPLICIT NONE
1115
1116	INTEGER(iwp) :: i !< loop index
1117	INTEGER(iwp) :: j !< loop index
1118	INTEGER(iwp) :: k !< loop index
1119	INTEGER(iwp) :: n !< loop index
1120
1121	REAL(wp) :: dum !< dummy argument
1122	REAL(wp) :: pt_l !< local scalar for pt used in equation of state function
1123	REAL(wp) :: sa_l !< local scalar for sa used in equation of state function
1124
1125	REAL(wp), DIMENSION(nzb:nzt+1) :: rho_ocean_init !< local array for initial density
1126
1127	ALLOCATE( hyp(nzb:nzt+1) )
1128
1129
1130	!
1131	!-- In case of no restart run, calculate the inital salinity profilevcusing the
1132	!-- given salinity gradients
1133	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
1134
1135	sa_init = sa_surface
1136	!
1137	!-- Last arguments gives back the gradient at top level to be used as
1138	!-- possible Neumann boundary condition. This is not realized for the ocean
1139	!-- mode, therefore a dummy argument is used.
1140	CALL init_vertical_profiles( sa_vertical_gradient_level_ind, &
1141	sa_vertical_gradient_level, &
1142	sa_vertical_gradient, sa_init, &
1143	sa_surface, dum )
1144	ENDIF
1145
1146	!
1147	!-- Initialize required 3d-arrays
1148	IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
1149	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
1150	!
1151	!-- Initialization via computed 1D-model profiles
1152	IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 ) THEN
1153
1154	DO i = nxlg, nxrg
1155	DO j = nysg, nyng
1156	sa(:,j,i) = sa_init
1157	ENDDO
1158	ENDDO
1159
1160	ENDIF
1161	!
1162	!-- Store initial profiles for output purposes etc.
1163	!-- Store initial salinity profile
1164	hom(:,1,26,:) = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
1165	!
1166	!-- Initialize old and new time levels.
1167	tsa_m = 0.0_wp
1168	sa_p = sa
1169
1170	ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN
1171
1172	!
1173	!-- Initialize new time levels (only done in order to set boundary values
1174	!-- including ghost points)
1175	sa_p = sa
1176	!
1177	!-- Allthough tendency arrays are set in prognostic_equations, they have
1178	!-- have to be predefined here because they are used (but multiplied with 0)
1179	!-- there before they are set.
1180	tsa_m = 0.0_wp
1181
1182	ENDIF
1183
1184	!
1185	!-- Set water density near the ocean surface
1186	rho_surface = 1027.62_wp
1187
1188	!
1189	!-- Set kinematic viscosity to sea water at 20C.
1190	!-- This changes the default value that is given for air!
1191	molecular_viscosity = 1.05E-6_wp
1192
1193	!
1194	!-- Change sign of buoyancy/stability terms because density gradient is used
1195	!-- instead of the potential temperature gradient to calculate the buoyancy
1196	atmos_ocean_sign = -1.0_wp
1197
1198	!
1199	!-- Calculate initial vertical profile of hydrostatic pressure (in Pa)
1200	!-- and the reference density (used later in buoyancy term)
1201	!-- First step: Calculate pressure using reference density
1202	hyp(nzt+1) = surface_pressure * 100.0_wp
1203	hyp(nzt) = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1)
1204	rho_ocean_init(nzt) = rho_surface
1205	rho_ocean_init(nzt+1) = rho_surface ! only required for output
1206
1207	DO k = nzt-1, 1, -1
1208	hyp(k) = hyp(k+1) + rho_surface * g * dzu(k)
1209	ENDDO
1210	hyp(0) = hyp(1) + rho_surface * g * dzu(1)
1211
1212	!
1213	!-- Second step: Iteratively calculate in situ density (based on presssure)
1214	!-- and pressure (based on in situ density)
1215	DO n = 1, 5
1216
1217	rho_reference = rho_surface * 0.5_wp * dzu(nzt+1)
1218
1219	DO k = nzt, 0, -1
1220
1221	sa_l = 0.5_wp * ( sa_init(k) + sa_init(k+1) )
1222	pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) )
1223
1224	rho_ocean_init(k) = eqn_state_seawater_func( hyp(k), pt_l, sa_l )
1225
1226	rho_reference = rho_reference + rho_ocean_init(k) * dzu(k+1)
1227
1228	ENDDO
1229
1230	rho_reference = rho_reference / ( zw(nzt) - zu(nzb) )
1231
1232	hyp(nzt) = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1)
1233	DO k = nzt-1, 0, -1
1234	hyp(k) = hyp(k+1) + g * 0.5_wp * ( rho_ocean_init(k) &
1235	+ rho_ocean_init(k+1) ) * dzu(k+1)
1236	ENDDO
1237
1238	ENDDO
1239
1240	!
1241	!-- Calculate the reference potential density
1242	prho_reference = 0.0_wp
1243	DO k = 0, nzt
1244
1245	sa_l = 0.5_wp * ( sa_init(k) + sa_init(k+1) )
1246	pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) )
1247
1248	prho_reference = prho_reference + dzu(k+1) * &
1249	eqn_state_seawater_func( 0.0_wp, pt_l, sa_l )
1250
1251	ENDDO
1252
1253	prho_reference = prho_reference / ( zu(nzt) - zu(nzb) )
1254
1255	!
1256	!-- Calculate the 3d array of initial in situ and potential density,
1257	!-- based on the initial temperature and salinity profile
1258	CALL eqn_state_seawater
1259
1260	!
1261	!-- Store initial density profile
1262	hom(:,1,77,:) = SPREAD( rho_ocean_init(:), 2, statistic_regions+1 )
1263
1264	!
1265	!-- Set the reference state to be used in the buoyancy terms
1266	IF ( use_single_reference_value ) THEN
1267	ref_state(:) = rho_reference
1268	ELSE
1269	ref_state(:) = rho_ocean_init(:)
1270	ENDIF
1271
1272
1273	END SUBROUTINE ocean_init
1274
1275
1276	!------------------------------------------------------------------------------!
1277	! Description:
1278	! ------------
1279	!> Prognostic equation for salinity.
1280	!> Vector-optimized version
1281	!------------------------------------------------------------------------------!
1282	SUBROUTINE ocean_prognostic_equations
1283
1284	USE advec_s_bc_mod, &
1285	ONLY: advec_s_bc
1286
1287	USE advec_s_pw_mod, &
1288	ONLY: advec_s_pw
1289
1290	USE advec_s_up_mod, &
1291	ONLY: advec_s_up
1292
1293	USE advec_ws, &
1294	ONLY: advec_s_ws
1295
1296	USE arrays_3d, &
1297	ONLY: rdf_sc, tend, tsa_m
1298
1299	USE control_parameters, &
1300	ONLY: dt_3d, intermediate_timestep_count, intermediate_timestep_count_max, &
1301	scalar_advec, timestep_scheme, tsc, ws_scheme_sca
1302
1303	USE cpulog, &
1304	ONLY: cpu_log, log_point
1305
1306	USE diffusion_s_mod, &
1307	ONLY: diffusion_s
1308
1309	USE indices, &
1310	ONLY: nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
1311
1312	USE surface_mod, &
1313	ONLY: surf_def_v, surf_def_h, surf_lsm_h, surf_lsm_v, surf_usm_h, &
1314	surf_usm_v
1315
1316	USE user_actions_mod, &
1317	ONLY: user_actions
1318
1319	IMPLICIT NONE
1320
1321	INTEGER(iwp) :: i !< loop index
1322	INTEGER(iwp) :: j !< loop index
1323	INTEGER(iwp) :: k !< loop index
1324
1325	REAL(wp) :: sbt !< weighting factor for sub-time step
1326
1327	!
1328	!-- Compute prognostic equations for the ocean mode
1329	!-- First, start with salinity
1330	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
1331
1332	!
1333	!-- sa-tendency terms with communication
1334	sbt = tsc(2)
1335	IF ( scalar_advec == 'bc-scheme' ) THEN
1336
1337	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1338	!
1339	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1340	sbt = 1.0_wp
1341	ENDIF
1342	tend = 0.0_wp
1343	CALL advec_s_bc( sa, 'sa' )
1344
1345	ENDIF
1346
1347	!
1348	!-- sa-tendency terms with no communication
1349	IF ( scalar_advec /= 'bc-scheme' ) THEN
1350	tend = 0.0_wp
1351	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1352	IF ( ws_scheme_sca ) THEN
1353	CALL advec_s_ws( sa, 'sa' )
1354	ELSE
1355	CALL advec_s_pw( sa )
1356	ENDIF
1357	ELSE
1358	CALL advec_s_up( sa )
1359	ENDIF
1360	ENDIF
1361
1362	CALL diffusion_s( sa, &
1363	surf_def_h(0)%sasws, surf_def_h(1)%sasws, &
1364	surf_def_h(2)%sasws, &
1365	surf_lsm_h%sasws, surf_usm_h%sasws, &
1366	surf_def_v(0)%sasws, surf_def_v(1)%sasws, &
1367	surf_def_v(2)%sasws, surf_def_v(3)%sasws, &
1368	surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws, &
1369	surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws, &
1370	surf_usm_v(0)%sasws, surf_usm_v(1)%sasws, &
1371	surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
1372
1373	CALL user_actions( 'sa-tendency' )
1374
1375	!
1376	!-- Prognostic equation for salinity
1377	DO i = nxl, nxr
1378	DO j = nys, nyn
1379	DO k = nzb+1, nzt
1380	sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
1381	tsc(3) * tsa_m(k,j,i) ) &
1382	- tsc(5) * rdf_sc(k) * &
1383	( sa(k,j,i) - sa_init(k) ) &
1384	) &
1385	* MERGE( 1.0_wp, 0.0_wp, &
1386	BTEST( wall_flags_0(k,j,i), 0 ) &
1387	)
1388	IF ( sa_p(k,j,i) < 0.0_wp ) sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
1389	ENDDO
1390	ENDDO
1391	ENDDO
1392
1393	!
1394	!-- Calculate tendencies for the next Runge-Kutta step
1395	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1396	IF ( intermediate_timestep_count == 1 ) THEN
1397	DO i = nxl, nxr
1398	DO j = nys, nyn
1399	DO k = nzb+1, nzt
1400	tsa_m(k,j,i) = tend(k,j,i)
1401	ENDDO
1402	ENDDO
1403	ENDDO
1404	ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )&
1405	THEN
1406	DO i = nxl, nxr
1407	DO j = nys, nyn
1408	DO k = nzb+1, nzt
1409	tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1410	5.3125_wp * tsa_m(k,j,i)
1411	ENDDO
1412	ENDDO
1413	ENDDO
1414	ENDIF
1415	ENDIF
1416
1417	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
1418
1419	!
1420	!-- Calculate density by the equation of state for seawater
1421	CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
1422	CALL eqn_state_seawater
1423	CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
1424
1425	END SUBROUTINE ocean_prognostic_equations
1426
1427
1428	!------------------------------------------------------------------------------!
1429	! Description:
1430	! ------------
1431	!> Prognostic equations for ocean mode (so far, salinity only)
1432	!> Cache-optimized version
1433	!------------------------------------------------------------------------------!
1434	SUBROUTINE ocean_prognostic_equations_ij( i, j, i_omp_start, tn )
1435
1436	USE advec_s_pw_mod, &
1437	ONLY: advec_s_pw
1438
1439	USE advec_s_up_mod, &
1440	ONLY: advec_s_up
1441
1442	USE advec_ws, &
1443	ONLY: advec_s_ws
1444
1445	USE arrays_3d, &
1446	ONLY: diss_l_sa, diss_s_sa, flux_l_sa, flux_s_sa, rdf_sc, tend, tsa_m
1447
1448	USE control_parameters, &
1449	ONLY: dt_3d, intermediate_timestep_count, &
1450	intermediate_timestep_count_max, timestep_scheme, tsc, &
1451	ws_scheme_sca
1452
1453	USE diffusion_s_mod, &
1454	ONLY: diffusion_s
1455
1456	USE indices, &
1457	ONLY: nzb, nzt, wall_flags_0
1458
1459	USE surface_mod, &
1460	ONLY: surf_def_v, surf_def_h, surf_lsm_h, surf_lsm_v, surf_usm_h, &
1461	surf_usm_v
1462
1463	USE user_actions_mod, &
1464	ONLY: user_actions
1465
1466	IMPLICIT NONE
1467
1468	INTEGER(iwp) :: i !< loop index x direction
1469	INTEGER(iwp) :: i_omp_start !< first loop index of i-loop in calling &
1470	!< routine prognostic_equations
1471	INTEGER(iwp) :: j !< loop index y direction
1472	INTEGER(iwp) :: k !< loop index z direction
1473	INTEGER(iwp) :: tn !< task number of openmp task
1474
1475	!
1476	!-- Compute prognostic equations for the ocean mode
1477	!-- First, start with tendency-terms for salinity
1478	tend(:,j,i) = 0.0_wp
1479	IF ( timestep_scheme(1:5) == 'runge' ) &
1480	THEN
1481	IF ( ws_scheme_sca ) THEN
1482	CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa, diss_s_sa, flux_l_sa, &
1483	diss_l_sa, i_omp_start, tn )
1484	ELSE
1485	CALL advec_s_pw( i, j, sa )
1486	ENDIF
1487	ELSE
1488	CALL advec_s_up( i, j, sa )
1489	ENDIF
1490	CALL diffusion_s( i, j, sa, &
1491	surf_def_h(0)%sasws, surf_def_h(1)%sasws, &
1492	surf_def_h(2)%sasws, &
1493	surf_lsm_h%sasws, surf_usm_h%sasws, &
1494	surf_def_v(0)%sasws, surf_def_v(1)%sasws, &
1495	surf_def_v(2)%sasws, surf_def_v(3)%sasws, &
1496	surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws, &
1497	surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws, &
1498	surf_usm_v(0)%sasws, surf_usm_v(1)%sasws, &
1499	surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
1500
1501	CALL user_actions( i, j, 'sa-tendency' )
1502
1503	!
1504	!-- Prognostic equation for salinity
1505	DO k = nzb+1, nzt
1506
1507	sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * &
1508	( tsc(2) * tend(k,j,i) + &
1509	tsc(3) * tsa_m(k,j,i) ) &
1510	- tsc(5) * rdf_sc(k) &
1511	* ( sa(k,j,i) - sa_init(k) ) &
1512	) * MERGE( 1.0_wp, 0.0_wp, &
1513	BTEST( wall_flags_0(k,j,i), 0 ) )
1514
1515	IF ( sa_p(k,j,i) < 0.0_wp ) sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
1516
1517	ENDDO
1518
1519	!
1520	!-- Calculate tendencies for the next Runge-Kutta step
1521	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1522	IF ( intermediate_timestep_count == 1 ) THEN
1523	DO k = nzb+1, nzt
1524	tsa_m(k,j,i) = tend(k,j,i)
1525	ENDDO
1526	ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )&
1527	THEN
1528	DO k = nzb+1, nzt
1529	tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1530	5.3125_wp * tsa_m(k,j,i)
1531	ENDDO
1532	ENDIF
1533	ENDIF
1534
1535	!
1536	!-- Calculate density by the equation of state for seawater
1537	CALL eqn_state_seawater( i, j )
1538
1539	END SUBROUTINE ocean_prognostic_equations_ij
1540
1541
1542	!------------------------------------------------------------------------------!
1543	! Description:
1544	! ------------
1545	!> Swapping of timelevels.
1546	!------------------------------------------------------------------------------!
1547	SUBROUTINE ocean_swap_timelevel( mod_count )
1548
1549	IMPLICIT NONE
1550
1551	INTEGER, INTENT(IN) :: mod_count !< flag defining where pointers point to
1552
1553	#if defined( __nopointer )
1554
1555	sa = sa_p
1556
1557	#else
1558
1559	SELECT CASE ( mod_count )
1560
1561	CASE ( 0 )
1562	sa => sa_1; sa_p => sa_2
1563
1564	CASE ( 1 )
1565	sa => sa_2; sa_p => sa_1
1566
1567	END SELECT
1568
1569	#endif
1570
1571	END SUBROUTINE ocean_swap_timelevel
1572
1573
1574	!------------------------------------------------------------------------------!
1575	! Description:
1576	! ------------
1577	!> This routine reads the respective restart data for the ocean module.
1578	!------------------------------------------------------------------------------!
1579	SUBROUTINE ocean_rrd_global( found )
1580
1581
1582	USE control_parameters, &
1583	ONLY: length, restart_string
1584
1585
1586	IMPLICIT NONE
1587
1588	LOGICAL, INTENT(OUT) :: found
1589
1590
1591	found = .TRUE.
1592
1593	SELECT CASE ( restart_string(1:length) )
1594
1595	CASE ( 'bc_sa_t' )
1596	READ ( 13 ) bc_sa_t
1597
1598	CASE ( 'bottom_salinityflux' )
1599	READ ( 13 ) bottom_salinityflux
1600
1601	CASE ( 'langmuir_number' )
1602	READ ( 13 ) langmuir_number
1603
1604	CASE ( 'sa_init' )
1605	READ ( 13 ) sa_init
1606
1607	CASE ( 'sa_surface' )
1608	READ ( 13 ) sa_surface
1609
1610	CASE ( 'sa_vertical_gradient' )
1611	READ ( 13 ) sa_vertical_gradient
1612
1613	CASE ( 'sa_vertical_gradient_level' )
1614	READ ( 13 ) sa_vertical_gradient_level
1615
1616	CASE ( 'stokes_force' )
1617	READ ( 13 ) stokes_force
1618
1619	CASE ( 'stokes_waveheight' )
1620	READ ( 13 ) stokes_waveheight
1621
1622	CASE ( 'stokes_wavelength' )
1623	READ ( 13 ) stokes_wavelength
1624
1625	CASE ( 'top_salinityflux' )
1626	READ ( 13 ) top_salinityflux
1627
1628	CASE ( 'wall_salinityflux' )
1629	READ ( 13 ) wall_salinityflux
1630
1631	CASE DEFAULT
1632
1633	found = .FALSE.
1634
1635	END SELECT
1636
1637	END SUBROUTINE ocean_rrd_global
1638
1639
1640	!------------------------------------------------------------------------------!
1641	! Description:
1642	! ------------
1643	!> This routine reads the respective restart data for the ocean module.
1644	!------------------------------------------------------------------------------!
1645	SUBROUTINE ocean_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, &
1646	nxr_on_file, nynf, nync, nyn_on_file, nysf, &
1647	nysc, nys_on_file, tmp_2d, tmp_3d, found )
1648
1649	USE averaging, &
1650	ONLY: rho_ocean_av, sa_av
1651
1652	USE control_parameters, &
1653	ONLY: length, restart_string
1654
1655	USE indices, &
1656	ONLY: nbgp, nxlg, nxrg, nyng, nysg, nzb, nzt
1657
1658	USE pegrid
1659
1660
1661	IMPLICIT NONE
1662
1663	INTEGER(iwp) :: i !<
1664	INTEGER(iwp) :: k !<
1665	INTEGER(iwp) :: nxlc !<
1666	INTEGER(iwp) :: nxlf !<
1667	INTEGER(iwp) :: nxl_on_file !<
1668	INTEGER(iwp) :: nxrc !<
1669	INTEGER(iwp) :: nxrf !<
1670	INTEGER(iwp) :: nxr_on_file !<
1671	INTEGER(iwp) :: nync !<
1672	INTEGER(iwp) :: nynf !<
1673	INTEGER(iwp) :: nyn_on_file !<
1674	INTEGER(iwp) :: nysc !<
1675	INTEGER(iwp) :: nysf !<
1676	INTEGER(iwp) :: nys_on_file !<
1677
1678	LOGICAL, INTENT(OUT) :: found
1679
1680	REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d !<
1681	REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !<
1682
1683
1684	found = .TRUE.
1685
1686	SELECT CASE ( restart_string(1:length) )
1687
1688	CASE ( 'rho_ocean_av' )
1689	IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
1690	ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1691	ENDIF
1692	IF ( k == 1 ) READ ( 13 ) tmp_3d
1693	rho_ocean_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1694	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1695
1696	CASE ( 'sa' )
1697	IF ( k == 1 ) READ ( 13 ) tmp_3d
1698	sa(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1699	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1700
1701	CASE ( 'sa_av' )
1702	IF ( .NOT. ALLOCATED( sa_av ) ) THEN
1703	ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1704	ENDIF
1705	IF ( k == 1 ) READ ( 13 ) tmp_3d
1706	sa_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
1707	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1708
1709	CASE DEFAULT
1710	found = .FALSE.
1711
1712	END SELECT
1713
1714	END SUBROUTINE ocean_rrd_local
1715
1716
1717	!------------------------------------------------------------------------------!
1718	! Description:
1719	! ------------
1720	!> This routine writes the respective restart data for the ocean module.
1721	!------------------------------------------------------------------------------!
1722	SUBROUTINE ocean_wrd_global
1723
1724
1725	IMPLICIT NONE
1726
1727	CALL wrd_write_string( 'bc_sa_t' )
1728	WRITE ( 14 ) bc_sa_t
1729
1730	CALL wrd_write_string( 'bottom_salinityflux' )
1731	WRITE ( 14 ) bottom_salinityflux
1732
1733	CALL wrd_write_string( 'langmuir_number' )
1734	WRITE ( 14 ) langmuir_number
1735
1736	CALL wrd_write_string( 'sa_init' )
1737	WRITE ( 14 ) sa_init
1738
1739	CALL wrd_write_string( 'sa_surface' )
1740	WRITE ( 14 ) sa_surface
1741
1742	CALL wrd_write_string( 'sa_vertical_gradient' )
1743	WRITE ( 14 ) sa_vertical_gradient
1744
1745	CALL wrd_write_string( 'sa_vertical_gradient_level' )
1746	WRITE ( 14 ) sa_vertical_gradient_level
1747
1748	CALL wrd_write_string( 'stokes_force' )
1749	WRITE ( 14 ) stokes_force
1750
1751	CALL wrd_write_string( 'stokes_waveheight' )
1752	WRITE ( 14 ) stokes_waveheight
1753
1754	CALL wrd_write_string( 'stokes_wavelength' )
1755	WRITE ( 14 ) stokes_wavelength
1756
1757	CALL wrd_write_string( 'top_salinityflux' )
1758	WRITE ( 14 ) top_salinityflux
1759
1760	CALL wrd_write_string( 'wall_salinityflux' )
1761	WRITE ( 14 ) wall_salinityflux
1762
1763	END SUBROUTINE ocean_wrd_global
1764
1765
1766	!------------------------------------------------------------------------------!
1767	! Description:
1768	! ------------
1769	!> This routine writes the respective restart data for the ocean module.
1770	!------------------------------------------------------------------------------!
1771	SUBROUTINE ocean_wrd_local
1772
1773	USE averaging, &
1774	ONLY: rho_ocean_av, sa_av
1775
1776	IMPLICIT NONE
1777
1778	IF ( ALLOCATED( rho_ocean_av ) ) THEN
1779	CALL wrd_write_string( 'rho_ocean_av' )
1780	WRITE ( 14 ) rho_ocean_av
1781	ENDIF
1782
1783	CALL wrd_write_string( 'sa' )
1784	WRITE ( 14 ) sa
1785
1786	IF ( ALLOCATED( sa_av ) ) THEN
1787	CALL wrd_write_string( 'sa_av' )
1788	WRITE ( 14 ) sa_av
1789	ENDIF
1790
1791	END SUBROUTINE ocean_wrd_local
1792
1793
1794	END MODULE ocean_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |