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

Last change on this file since 3352 was 3311, checked in by raasch, 6 years ago

Stokes drift is regarded in timestep calculation, check if ocean mode is used for invalid combinations

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