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

Last change on this file since 4631 was 4535, checked in by raasch, 4 years ago

bugfix for restart data format query

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