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

Last change on this file since 4517 was 4517, checked in by raasch, 14 months ago

added restart with MPI-IO for reading local arrays

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