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

Last change on this file since 4289 was 4272, checked in by schwenkel, 5 years ago

further modularization of boundary conditions: moving boundary conditions to their modules

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