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

Last change on this file since 4109 was 4109, checked in by suehring, 2 years ago

Control discretization of advection term: separate initialization of WS advection flags for momentum and scalars. In this context, resort the bits and do some minor formatting; Make initialization of scalar-advection flags more flexible, i.e. introduce an arguemnt list to indicate non-cyclic boundaries (required for decycled scalars such as chemical species or aerosols); Introduce extended 'degradation zones', where horizontal advection of passive scalars is discretized by first-order scheme at all grid points that in the vicinity of buildings (<= 3 grid points). Even though no building is within the numerical stencil, first-order scheme is used. At fourth and fifth grid point the order of the horizontal advection scheme is successively upgraded. These extended degradation zones are used to avoid stationary numerical oscillations, which are responsible for high concentration maxima that may appear under shear-free stable conditions. Therefore, an additional 3D interger array used to store control flags is introduced; Change interface for scalar advection routine; Bugfix, avoid uninitialized value sk_num in vector version of WS scalar advection; Chemistry: Decycling boundary conditions are only set at the ghost points not on the prognostic grid points; Land-surface model: Relax checks for non-consistent initialization in case static or dynamic input is provided. For example, soil_temperature or deep_soil_temperature is not mandatory any more if dynamic input is available. Also, improper settings of x_type in namelist are only checked if no static file is available.

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