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

Last change on this file since 4265 was 4196, checked in by gronemeier, 5 years ago

Consider rotation of model domain for claculation of Coriolis force:

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