source: palm/trunk/SOURCE/ocean_mod.f90

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

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

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