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

Last change on this file since 4495 was 4495, checked in by raasch, 5 years ago

restart data handling with MPI-IO added, first part

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