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

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