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

Last change on this file since 4343 was 4329, checked in by motisi, 5 years ago

Renamed wall_flags_0 to wall_flags_static_0

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