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

Last change on this file since 3303 was 3303, checked in by raasch, 3 years ago

salinity allowed to be switched off, bugfix for swapping in case of ocean mode

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