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

Last change on this file since 4677 was 4671, checked in by pavelkrc, 4 years ago

Radiative transfer model RTM version 4.1

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