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

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

Craik-Leibovich force and wave breaking effects added to ocean mode, header output for ocean mode

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