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

Last change on this file since 4180 was 4180, checked in by scharf, 5 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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