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

Last change on this file since 4477 was 4370, checked in by raasch, 5 years ago

bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added

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