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

Last change on this file since 4361 was 4346, checked in by motisi, 4 years ago

Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0

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