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

Last change on this file since 3768 was 3767, checked in by raasch, 2 years ago

unused variables removed from rrd-subroutines parameter list

  • Property svn:keywords set to Id
File size: 82.0 KB
Line 
1!> @file ocean_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2017-2018 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: ocean_mod.f90 3767 2019-02-27 08:18:02Z raasch $
27! unused variable for file index and tmp_2d removed from rrd-subroutine parameter
28! list
29!
30! 3719 2019-02-06 13:10:18Z kanani
31! Changed log_point to log_point_s, otherwise this overlaps with
32! 'all progn.equations' cpu measurement.
33!
34! 3684 2019-01-20 20:20:58Z knoop
35! nopointer option removed
36!
37! 3614 2018-12-10 07:05:46Z raasch
38! unused variables removed
39!
40! 3568 2018-11-27 16:07:59Z raasch
41! bugifx: calculate equation of state for seawater even if salinity is switched
42!         off
43!
44! 3421 2018-10-24 18:39:32Z gronemeier
45! Renamed output variables
46!
47! 3381 2018-10-19 13:09:06Z raasch
48! spin-up cooling for ocean surface implemented, see new parameter
49! surface_cooling_spinup_time
50!
51! 3311 2018-10-05 12:34:56Z raasch
52! check if ocean mode is used for invalid combinations
53!
54! 3303 2018-10-03 12:04:15Z raasch
55! salinity allowed to be switched off
56!
57! 3302 2018-10-03 02:39:40Z raasch
58! Craik Leibovich force (Stokes drift) + wave breaking effect added
59!
60! 3294 2018-10-01 02:37:10Z raasch
61! initial revision
62!
63!
64! Authors:
65! --------
66! @author Siegfried Raasch
67!
68! Description:
69! ------------
70!> This module contains relevant code for PALM's ocean mode, e.g. the
71!> prognostic equation for salinity, the equation of state for seawater,
72!> the Craik Leibovich force (Stokes force), and wave breaking effects
73!------------------------------------------------------------------------------!
74 MODULE ocean_mod
75 
76
77    USE arrays_3d,                                                             &
78        ONLY:  prho, prho_1, rho_ocean, rho_1, sa, sa_init, sa_1, sa_2, sa_3,  &
79               sa_p, tsa_m
80
81    USE control_parameters,                                                    &
82        ONLY:  atmos_ocean_sign, bottom_salinityflux,                          &
83               constant_top_salinityflux, ocean_mode, top_salinityflux,        &
84               wall_salinityflux
85
86    USE kinds
87
88
89    IMPLICIT NONE
90
91    CHARACTER (LEN=20) ::  bc_sa_t = 'neumann'  !< namelist parameter
92
93    INTEGER(iwp) ::  ibc_sa_t   !< integer flag for bc_sa_t
94    INTEGER(iwp) ::  iran_ocean = -1234567  !< random number used for wave breaking effects
95
96    INTEGER(iwp) ::  sa_vertical_gradient_level_ind(10) = -9999  !< grid index values of sa_vertical_gradient_level(s)
97
98    LOGICAL ::  salinity = .TRUE.             !< switch for using salinity
99    LOGICAL ::  stokes_force = .FALSE.        !< switch to switch on the Stokes force
100    LOGICAL ::  wave_breaking = .FALSE.       !< switch to switch on wave breaking effects
101    LOGICAL ::  surface_cooling_switched_off = .FALSE.  !< variable to check if surface heat flux has been switched off
102
103    REAL(wp) ::  alpha_wave_breaking = 3.0_wp !< coefficient for wave breaking generated turbulence from Noh et al. (2004), JPO
104    REAL(wp) ::  prho_reference               !< reference state of potential density at ocean surface
105    REAL(wp) ::  sa_surface = 35.0_wp         !< surface salinity, namelist parameter
106    REAL(wp) ::  sa_vertical_gradient(10) = 0.0_wp               !< namelist parameter
107    REAL(wp) ::  sa_vertical_gradient_level(10) = -999999.9_wp   !< namelist parameter
108    REAL(wp) ::  stokes_waveheight = 0.0_wp  !< wave height assumed for Stokes drift velocity
109    REAL(wp) ::  stokes_wavelength = 0.0_wp  !< wavelength assumed for Stokes drift velocity
110    REAL(wp) ::  surface_cooling_spinup_time = 999999.9_wp  !< time after which surface heat flux is switched off
111    REAL(wp) ::  timescale_wave_breaking     !< time scale of random forcing
112    REAL(wp) ::  u_star_wave_breaking        !< to store the absolute value of friction velocity at the ocean surface
113
114    REAL(wp), DIMENSION(12), PARAMETER ::  nom =                               &
115                          (/ 9.9984085444849347D2,   7.3471625860981584D0,     &
116                            -5.3211231792841769D-2,  3.6492439109814549D-4,    &
117                             2.5880571023991390D0,  -6.7168282786692354D-3,    &
118                             1.9203202055760151D-3,  1.1798263740430364D-2,    &
119                             9.8920219266399117D-8,  4.6996642771754730D-6,    &
120                            -2.5862187075154352D-8, -3.2921414007960662D-12 /)
121                          !< constants used in equation of state for seawater
122
123    REAL(wp), DIMENSION(13), PARAMETER ::  den =                               &
124                          (/ 1.0D0,                  7.2815210113327091D-3,    &
125                            -4.4787265461983921D-5,  3.3851002965802430D-7,    &
126                             1.3651202389758572D-10, 1.7632126669040377D-3,    &
127                            -8.8066583251206474D-6, -1.8832689434804897D-10,   &
128                             5.7463776745432097D-6,  1.4716275472242334D-9,    &
129                             6.7103246285651894D-6, -2.4461698007024582D-17,   &
130                            -9.1534417604289062D-18 /)
131                          !< constants used in equation of state for seawater
132
133    SAVE
134
135    PUBLIC ::  bc_sa_t, ibc_sa_t, prho_reference, sa_surface,                  &
136               sa_vertical_gradient, sa_vertical_gradient_level,               &
137               sa_vertical_gradient_level_ind, stokes_force, wave_breaking
138
139
140    INTERFACE eqn_state_seawater
141       MODULE PROCEDURE eqn_state_seawater
142       MODULE PROCEDURE eqn_state_seawater_ij
143    END INTERFACE eqn_state_seawater
144
145    INTERFACE eqn_state_seawater_func
146       MODULE PROCEDURE eqn_state_seawater_func
147    END INTERFACE eqn_state_seawater_func
148
149    INTERFACE ocean_check_parameters
150       MODULE PROCEDURE ocean_check_parameters
151    END INTERFACE ocean_check_parameters
152
153    INTERFACE ocean_check_data_output
154       MODULE PROCEDURE ocean_check_data_output
155    END INTERFACE ocean_check_data_output
156
157    INTERFACE ocean_check_data_output_pr
158       MODULE PROCEDURE ocean_check_data_output_pr
159    END INTERFACE ocean_check_data_output_pr
160
161    INTERFACE ocean_define_netcdf_grid
162       MODULE PROCEDURE ocean_define_netcdf_grid
163    END INTERFACE ocean_define_netcdf_grid
164
165    INTERFACE ocean_data_output_2d
166       MODULE PROCEDURE ocean_data_output_2d
167    END INTERFACE ocean_data_output_2d
168
169    INTERFACE ocean_data_output_3d
170       MODULE PROCEDURE ocean_data_output_3d
171    END INTERFACE ocean_data_output_3d
172
173    INTERFACE ocean_header
174       MODULE PROCEDURE ocean_header
175    END INTERFACE ocean_header
176
177    INTERFACE ocean_init
178       MODULE PROCEDURE ocean_init
179    END INTERFACE ocean_init
180
181    INTERFACE ocean_init_arrays
182       MODULE PROCEDURE ocean_init_arrays
183    END INTERFACE ocean_init_arrays
184
185    INTERFACE ocean_parin
186       MODULE PROCEDURE ocean_parin
187    END INTERFACE ocean_parin
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_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, nyng, 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    prho => prho_1
1193    rho_ocean  => rho_1  ! routines calc_mean_profile and diffusion_e require
1194                         ! density to be a pointer
1195
1196!
1197!-- Initial assignment of pointers
1198    IF ( salinity )  THEN
1199       sa => sa_1;  sa_p => sa_2;  tsa_m => sa_3
1200    ELSE
1201       sa => sa_1;  sa_p => sa_1;  tsa_m => sa_3
1202    ENDIF
1203
1204 END SUBROUTINE ocean_init_arrays
1205
1206
1207!------------------------------------------------------------------------------!
1208! Description:
1209! ------------
1210!> Initialization of quantities needed for the ocean mode
1211!------------------------------------------------------------------------------!
1212 SUBROUTINE ocean_init
1213
1214
1215    USE arrays_3d,                                                             &
1216        ONLY:  dzu, dzw, hyp, pt_init, ref_state, u_stokes_zu, u_stokes_zw,    &
1217               v_stokes_zu, v_stokes_zw, zu, zw
1218
1219    USE basic_constants_and_equations_mod,                                     &
1220        ONLY:  g
1221
1222    USE basic_constants_and_equations_mod,                                     &
1223        ONLY:  pi
1224
1225    USE control_parameters,                                                    &
1226        ONLY:  initializing_actions, molecular_viscosity, rho_surface,         &
1227               rho_reference, surface_pressure, top_momentumflux_u,            &
1228               top_momentumflux_v, use_single_reference_value
1229
1230    USE indices,                                                               &
1231        ONLY:  nxl, nxlg, nxrg, nyng, nys, nysg, nzb, nzt
1232
1233    USE kinds
1234
1235    USE pegrid,                                                                &
1236        ONLY:  myid
1237
1238    USE statistics,                                                            &
1239        ONLY:  hom, statistic_regions
1240
1241    IMPLICIT NONE
1242
1243    INTEGER(iwp) ::  i  !< loop index
1244    INTEGER(iwp) ::  j  !< loop index
1245    INTEGER(iwp) ::  k  !< loop index
1246    INTEGER(iwp) ::  n  !< loop index
1247
1248    REAL(wp) ::  alpha !< angle of surface stress
1249    REAL(wp) ::  dum   !< dummy argument
1250    REAL(wp) ::  pt_l  !< local scalar for pt used in equation of state function
1251    REAL(wp) ::  sa_l  !< local scalar for sa used in equation of state function
1252    REAL(wp) ::  velocity_amplitude  !< local scalar for amplitude of Stokes drift velocity
1253    REAL(wp) ::  x     !< temporary variable to store surface stress along x
1254    REAL(wp) ::  y     !< temporary variable to store surface stress along y
1255
1256    REAL(wp), DIMENSION(nzb:nzt+1) ::  rho_ocean_init  !< local array for initial density
1257
1258    ALLOCATE( hyp(nzb:nzt+1) )
1259
1260
1261!
1262!-- In case of no restart run, calculate the inital salinity profilevcusing the
1263!-- given salinity gradients
1264    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1265
1266       sa_init = sa_surface
1267!
1268!--    Last arguments gives back the gradient at top level to be used as
1269!--    possible Neumann boundary condition. This is not realized for the ocean
1270!--    mode, therefore a dummy argument is used.
1271       IF ( salinity )  THEN
1272          CALL init_vertical_profiles( sa_vertical_gradient_level_ind,          &
1273                                       sa_vertical_gradient_level,              &
1274                                       sa_vertical_gradient, sa_init,           &
1275                                       sa_surface, dum )
1276       ENDIF
1277    ENDIF
1278
1279!
1280!-- Initialize required 3d-arrays
1281    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
1282         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1283!
1284!--    Initialization via computed 1D-model profiles
1285       IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 )  THEN
1286
1287          DO  i = nxlg, nxrg
1288             DO  j = nysg, nyng
1289                sa(:,j,i) = sa_init
1290             ENDDO
1291          ENDDO
1292
1293       ENDIF
1294!
1295!--    Store initial profiles for output purposes etc.
1296!--    Store initial salinity profile
1297       hom(:,1,26,:)  = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
1298!
1299!--    Initialize old and new time levels.
1300       tsa_m = 0.0_wp
1301       sa_p  = sa
1302
1303    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
1304
1305!
1306!--    Initialize new time levels (only done in order to set boundary values
1307!--    including ghost points)
1308       sa_p = sa
1309!
1310!--    Allthough tendency arrays are set in prognostic_equations, they have
1311!--    have to be predefined here because they are used (but multiplied with 0)
1312!--    there before they are set.
1313       tsa_m = 0.0_wp
1314
1315    ENDIF
1316
1317!
1318!-- Set water density near the ocean surface
1319    rho_surface = 1027.62_wp
1320
1321!
1322!-- Set kinematic viscosity to sea water at 20C.
1323!-- This changes the default value that is given for air!
1324    molecular_viscosity = 1.05E-6_wp
1325
1326!
1327!-- Change sign of buoyancy/stability terms because density gradient is used
1328!-- instead of the potential temperature gradient to calculate the buoyancy
1329    atmos_ocean_sign = -1.0_wp
1330
1331!
1332!-- Calculate initial vertical profile of hydrostatic pressure (in Pa)
1333!-- and the reference density (used later in buoyancy term)
1334!-- First step: Calculate pressure using reference density
1335    hyp(nzt+1) = surface_pressure * 100.0_wp
1336    hyp(nzt)   = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1)
1337    rho_ocean_init(nzt)   = rho_surface
1338    rho_ocean_init(nzt+1) = rho_surface  ! only required for output
1339
1340    DO  k = nzt-1, 1, -1
1341       hyp(k) = hyp(k+1) + rho_surface * g * dzu(k)
1342    ENDDO
1343    hyp(0) = hyp(1) + rho_surface * g * dzu(1)
1344
1345!
1346!-- Second step: Iteratively calculate in situ density (based on presssure)
1347!-- and pressure (based on in situ density)
1348    DO  n = 1, 5
1349
1350       rho_reference = rho_surface * 0.5_wp * dzu(nzt+1)
1351
1352       DO  k = nzt, 0, -1
1353
1354          sa_l = 0.5_wp * ( sa_init(k) + sa_init(k+1) )
1355          pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) )
1356
1357          rho_ocean_init(k) = eqn_state_seawater_func( hyp(k), pt_l, sa_l )
1358
1359          rho_reference = rho_reference + rho_ocean_init(k) * dzu(k+1)
1360
1361       ENDDO
1362
1363       rho_reference = rho_reference / ( zw(nzt) - zu(nzb) )
1364
1365       hyp(nzt) = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1)
1366       DO  k = nzt-1, 0, -1
1367          hyp(k) = hyp(k+1) + g * 0.5_wp * ( rho_ocean_init(k)                 &
1368                                           + rho_ocean_init(k+1) ) * dzu(k+1)
1369       ENDDO
1370
1371    ENDDO
1372
1373!
1374!-- Calculate the reference potential density
1375    prho_reference = 0.0_wp
1376    DO  k = 0, nzt
1377
1378       sa_l = 0.5_wp * ( sa_init(k) + sa_init(k+1) )
1379       pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) )
1380
1381       prho_reference = prho_reference + dzu(k+1) * &
1382                        eqn_state_seawater_func( 0.0_wp, pt_l, sa_l )
1383
1384    ENDDO
1385
1386    prho_reference = prho_reference / ( zu(nzt) - zu(nzb) )
1387
1388!
1389!-- Calculate the 3d array of initial in situ and potential density,
1390!-- based on the initial temperature and salinity profile
1391    CALL eqn_state_seawater
1392
1393!
1394!-- Store initial density profile
1395    hom(:,1,77,:)  = SPREAD( rho_ocean_init(:), 2, statistic_regions+1 )
1396
1397!
1398!-- Set the reference state to be used in the buoyancy terms
1399    IF ( use_single_reference_value )  THEN
1400       ref_state(:) = rho_reference
1401    ELSE
1402       ref_state(:) = rho_ocean_init(:)
1403    ENDIF
1404
1405!
1406!-- Calculate the Stokes drift velocity profile
1407    IF ( stokes_force )  THEN
1408
1409!
1410!--    First, calculate angle of surface stress
1411       x = -top_momentumflux_u
1412       y = -top_momentumflux_v
1413       IF ( x == 0.0_wp )  THEN
1414          IF ( y > 0.0_wp )  THEN
1415             alpha = pi / 2.0_wp
1416          ELSEIF ( y < 0.0_wp )  THEN
1417             alpha = 3.0_wp * pi / 2.0_wp
1418          ENDIF
1419       ELSE
1420          IF ( x < 0.0_wp )  THEN
1421             alpha = ATAN( y / x ) + pi
1422          ELSE
1423             IF ( y < 0.0_wp )  THEN
1424                alpha = ATAN( y / x ) + 2.0_wp * pi
1425             ELSE
1426                alpha = ATAN( y / x )
1427             ENDIF
1428          ENDIF
1429       ENDIF
1430
1431       velocity_amplitude = ( pi * stokes_waveheight / stokes_wavelength )**2 *&
1432                            SQRT( g * stokes_wavelength / ( 2.0_wp * pi ) )
1433
1434       DO  k = nzb, nzt
1435          u_stokes_zu(k) = velocity_amplitude * COS( alpha ) *                 &
1436                           EXP( 4.0_wp * pi * zu(k) / stokes_wavelength )
1437          u_stokes_zw(k) = velocity_amplitude * COS( alpha ) *                 &
1438                           EXP( 4.0_wp * pi * zw(k) / stokes_wavelength )
1439          v_stokes_zu(k) = velocity_amplitude * SIN( alpha ) *                 &
1440                           EXP( 4.0_wp * pi * zu(k) / stokes_wavelength )
1441          v_stokes_zw(k) = velocity_amplitude * SIN( alpha ) *                 &
1442                           EXP( 4.0_wp * pi * zw(k) / stokes_wavelength )
1443       ENDDO
1444       u_stokes_zu(nzt+1) = u_stokes_zw(nzt) ! because zu(nzt+1) changes the sign
1445       u_stokes_zw(nzt+1) = u_stokes_zw(nzt) ! because zw(nzt+1) changes the sign
1446       v_stokes_zu(nzt+1) = v_stokes_zw(nzt) ! because zu(nzt+1) changes the sign
1447       v_stokes_zw(nzt+1) = v_stokes_zw(nzt) ! because zw(nzt+1) changes the sign
1448
1449    ENDIF
1450
1451!
1452!-- Wave breaking effects
1453    IF ( wave_breaking )  THEN
1454!
1455!--    Calculate friction velocity at ocean surface
1456       u_star_wave_breaking = SQRT( SQRT( top_momentumflux_u**2 +              &
1457                                          top_momentumflux_v**2 ) )
1458!
1459!--    Set the time scale of random forcing. The vertical grid spacing at the
1460!--    ocean surface is assumed as the length scale of turbulence.
1461!--    Formula follows Noh et al. (2004), JPO
1462       timescale_wave_breaking = 0.1_wp * dzw(nzt) / alpha_wave_breaking /     &
1463                                 u_star_wave_breaking
1464!
1465!--    Set random number seeds differently on the processor cores in order to
1466!--    create different random number sequences
1467       iran_ocean = iran_ocean + myid
1468    ENDIF
1469
1470 END SUBROUTINE ocean_init
1471
1472
1473!------------------------------------------------------------------------------!
1474! Description:
1475! ------------
1476!> Prognostic equation for salinity.
1477!> Vector-optimized version
1478!------------------------------------------------------------------------------!
1479 SUBROUTINE ocean_prognostic_equations
1480
1481    USE advec_s_bc_mod,                                                        &
1482        ONLY:  advec_s_bc
1483
1484    USE advec_s_pw_mod,                                                        &
1485        ONLY:  advec_s_pw
1486
1487    USE advec_s_up_mod,                                                        &
1488        ONLY:  advec_s_up
1489
1490    USE advec_ws,                                                              &
1491        ONLY:  advec_s_ws
1492
1493    USE arrays_3d,                                                             &
1494        ONLY:  rdf_sc, tend, tsa_m
1495
1496    USE control_parameters,                                                    &
1497        ONLY:  dt_3d, intermediate_timestep_count,                             &
1498               intermediate_timestep_count_max, scalar_advec, simulated_time,  &
1499               timestep_scheme, tsc, ws_scheme_sca
1500
1501    USE cpulog,                                                                &
1502        ONLY:  cpu_log, log_point_s
1503
1504    USE diffusion_s_mod,                                                       &
1505        ONLY:  diffusion_s
1506
1507    USE indices,                                                               &
1508        ONLY:  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
1509
1510    USE surface_mod,                                                           &
1511        ONLY:  surf_def_v, surf_def_h, surf_lsm_h, surf_lsm_v, surf_usm_h,     &
1512               surf_usm_v
1513
1514    IMPLICIT NONE
1515
1516    INTEGER(iwp) ::  i       !< loop index
1517    INTEGER(iwp) ::  j       !< loop index
1518    INTEGER(iwp) ::  k       !< loop index
1519
1520    REAL(wp)     ::  sbt     !< weighting factor for sub-time step
1521
1522!
1523!-- Switch of the surface heat flux, if requested
1524    IF ( surface_cooling_spinup_time /= 999999.9_wp )  THEN
1525       IF ( .NOT. surface_cooling_switched_off  .AND.                          &
1526            simulated_time >= surface_cooling_spinup_time )  THEN
1527
1528          surf_def_h(2)%shf = 0.0_wp
1529          surface_cooling_switched_off = .TRUE.
1530
1531       ENDIF
1532    ENDIF
1533
1534!
1535!-- Compute prognostic equations for the ocean mode
1536!-- First, start with salinity
1537    IF ( salinity )  THEN
1538
1539       CALL cpu_log( log_point_s(20), 'sa-equation', 'start' )
1540
1541!
1542!--    sa-tendency terms with communication
1543       sbt = tsc(2)
1544       IF ( scalar_advec == 'bc-scheme' )  THEN
1545
1546          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1547!
1548!--          Bott-Chlond scheme always uses Euler time step. Thus:
1549             sbt = 1.0_wp
1550          ENDIF
1551          tend = 0.0_wp
1552          CALL advec_s_bc( sa, 'sa' )
1553
1554       ENDIF
1555
1556!
1557!--    sa-tendency terms with no communication
1558       IF ( scalar_advec /= 'bc-scheme' )  THEN
1559          tend = 0.0_wp
1560          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1561             IF ( ws_scheme_sca )  THEN
1562                CALL advec_s_ws( sa, 'sa' )
1563             ELSE
1564                CALL advec_s_pw( sa )
1565             ENDIF
1566          ELSE
1567             CALL advec_s_up( sa )
1568          ENDIF
1569       ENDIF
1570
1571       CALL diffusion_s( sa,                                                   &
1572                         surf_def_h(0)%sasws, surf_def_h(1)%sasws,             &
1573                         surf_def_h(2)%sasws,                                  &
1574                         surf_lsm_h%sasws,    surf_usm_h%sasws,                &
1575                         surf_def_v(0)%sasws, surf_def_v(1)%sasws,             &
1576                         surf_def_v(2)%sasws, surf_def_v(3)%sasws,             &
1577                         surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws,             &
1578                         surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws,             &
1579                         surf_usm_v(0)%sasws, surf_usm_v(1)%sasws,             &
1580                         surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
1581
1582!       CALL user_actions( 'sa-tendency' ) ToDo: find general solution for dependency between modules
1583
1584!
1585!--    Prognostic equation for salinity
1586       DO  i = nxl, nxr
1587          DO  j = nys, nyn
1588             DO  k = nzb+1, nzt
1589                sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +      &
1590                                                      tsc(3) * tsa_m(k,j,i) )  &
1591                                                  - tsc(5) * rdf_sc(k) *       &
1592                                                    ( sa(k,j,i) - sa_init(k) ) &
1593                                          )                                    &
1594                                            * MERGE( 1.0_wp, 0.0_wp,           &
1595                                               BTEST( wall_flags_0(k,j,i), 0 ) &
1596                                                   )
1597                IF ( sa_p(k,j,i) < 0.0_wp )  sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
1598             ENDDO
1599          ENDDO
1600       ENDDO
1601
1602!
1603!--    Calculate tendencies for the next Runge-Kutta step
1604       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1605          IF ( intermediate_timestep_count == 1 )  THEN
1606             DO  i = nxl, nxr
1607                DO  j = nys, nyn
1608                   DO  k = nzb+1, nzt
1609                      tsa_m(k,j,i) = tend(k,j,i)
1610                   ENDDO
1611                ENDDO
1612             ENDDO
1613          ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) &
1614          THEN
1615             DO  i = nxl, nxr
1616                DO  j = nys, nyn
1617                   DO  k = nzb+1, nzt
1618                      tsa_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
1619                                        5.3125_wp * tsa_m(k,j,i)
1620                   ENDDO
1621                ENDDO
1622             ENDDO
1623          ENDIF
1624       ENDIF
1625
1626       CALL cpu_log( log_point_s(20), 'sa-equation', 'stop' )
1627
1628    ENDIF
1629
1630!
1631!-- Calculate density by the equation of state for seawater
1632    CALL cpu_log( log_point_s(21), 'eqns-seawater', 'start' )
1633    CALL eqn_state_seawater
1634    CALL cpu_log( log_point_s(21), 'eqns-seawater', 'stop' )
1635
1636 END SUBROUTINE ocean_prognostic_equations
1637
1638
1639!------------------------------------------------------------------------------!
1640! Description:
1641! ------------
1642!> Prognostic equations for ocean mode (so far, salinity only)
1643!> Cache-optimized version
1644!------------------------------------------------------------------------------!
1645 SUBROUTINE ocean_prognostic_equations_ij( i, j, i_omp_start, tn )
1646
1647    USE advec_s_pw_mod,                                                        &
1648        ONLY:  advec_s_pw
1649
1650    USE advec_s_up_mod,                                                        &
1651        ONLY:  advec_s_up
1652
1653    USE advec_ws,                                                              &
1654        ONLY:  advec_s_ws
1655
1656    USE arrays_3d,                                                             &
1657        ONLY:  diss_l_sa, diss_s_sa, flux_l_sa, flux_s_sa, rdf_sc, tend, tsa_m
1658
1659    USE control_parameters,                                                    &
1660        ONLY:  dt_3d, intermediate_timestep_count,                             &
1661               intermediate_timestep_count_max, simulated_time,                &
1662               timestep_scheme, tsc, ws_scheme_sca
1663
1664    USE diffusion_s_mod,                                                       &
1665        ONLY:  diffusion_s
1666
1667    USE indices,                                                               &
1668        ONLY:  nzb, nzt, wall_flags_0
1669
1670    USE surface_mod,                                                           &
1671        ONLY:  surf_def_v, surf_def_h, surf_lsm_h, surf_lsm_v, surf_usm_h,     &
1672               surf_usm_v
1673
1674    IMPLICIT NONE
1675
1676    INTEGER(iwp) ::  i             !< loop index x direction
1677    INTEGER(iwp) ::  i_omp_start   !< first loop index of i-loop in calling    &
1678                                   !< routine prognostic_equations
1679    INTEGER(iwp) ::  j             !< loop index y direction
1680    INTEGER(iwp) ::  k             !< loop index z direction
1681    INTEGER(iwp) ::  tn            !< task number of openmp task
1682
1683
1684!
1685!-- Switch of the surface heat flux, if requested
1686    IF ( surface_cooling_spinup_time /= 999999.9_wp )  THEN
1687       IF ( .NOT. surface_cooling_switched_off  .AND.                          &
1688            simulated_time >= surface_cooling_spinup_time )  THEN
1689
1690          surf_def_h(2)%shf = 0.0_wp
1691          surface_cooling_switched_off = .TRUE.
1692
1693       ENDIF
1694    ENDIF
1695
1696!
1697!-- Compute prognostic equations for the ocean mode
1698!-- First, start with tendency-terms for salinity
1699    IF ( salinity )  THEN
1700
1701       tend(:,j,i) = 0.0_wp
1702       IF ( timestep_scheme(1:5) == 'runge' ) &
1703       THEN
1704          IF ( ws_scheme_sca )  THEN
1705             CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa,  diss_s_sa, flux_l_sa,&
1706                              diss_l_sa, i_omp_start, tn  )
1707          ELSE
1708             CALL advec_s_pw( i, j, sa )
1709          ENDIF
1710       ELSE
1711          CALL advec_s_up( i, j, sa )
1712       ENDIF
1713       CALL diffusion_s( i, j, sa,                                             &
1714                         surf_def_h(0)%sasws, surf_def_h(1)%sasws,             &
1715                         surf_def_h(2)%sasws,                                  &
1716                         surf_lsm_h%sasws,    surf_usm_h%sasws,                &
1717                         surf_def_v(0)%sasws, surf_def_v(1)%sasws,             &
1718                         surf_def_v(2)%sasws, surf_def_v(3)%sasws,             &
1719                         surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws,             &
1720                         surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws,             &
1721                         surf_usm_v(0)%sasws, surf_usm_v(1)%sasws,             &
1722                         surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
1723
1724!       CALL user_actions( i, j, 'sa-tendency' ) ToDo: find general solution for dependency between modules
1725
1726!
1727!--    Prognostic equation for salinity
1728       DO  k = nzb+1, nzt
1729
1730          sa_p(k,j,i) = sa(k,j,i) + ( dt_3d *                                  &
1731                                              ( tsc(2) * tend(k,j,i) +         &
1732                                                tsc(3) * tsa_m(k,j,i) )        &
1733                                    - tsc(5) * rdf_sc(k)                       &
1734                                             * ( sa(k,j,i) - sa_init(k) )      &
1735                                    ) * MERGE( 1.0_wp, 0.0_wp,                 &
1736                                               BTEST( wall_flags_0(k,j,i), 0 ) )
1737
1738          IF ( sa_p(k,j,i) < 0.0_wp )  sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
1739
1740       ENDDO
1741
1742!
1743!--    Calculate tendencies for the next Runge-Kutta step
1744       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1745          IF ( intermediate_timestep_count == 1 )  THEN
1746             DO  k = nzb+1, nzt
1747                tsa_m(k,j,i) = tend(k,j,i)
1748             ENDDO
1749          ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) &
1750          THEN
1751             DO  k = nzb+1, nzt
1752                tsa_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +                    &
1753                                  5.3125_wp * tsa_m(k,j,i)
1754             ENDDO
1755          ENDIF
1756       ENDIF
1757
1758    ENDIF
1759
1760!
1761!-- Calculate density by the equation of state for seawater
1762    CALL eqn_state_seawater( i, j )
1763
1764 END SUBROUTINE ocean_prognostic_equations_ij
1765
1766
1767!------------------------------------------------------------------------------!
1768! Description:
1769! ------------
1770!> Swapping of timelevels.
1771!------------------------------------------------------------------------------!
1772 SUBROUTINE ocean_swap_timelevel( mod_count )
1773
1774    IMPLICIT NONE
1775
1776    INTEGER, INTENT(IN) ::  mod_count  !< flag defining where pointers point to
1777
1778
1779    SELECT CASE ( mod_count )
1780
1781       CASE ( 0 )
1782          IF ( salinity )  THEN
1783             sa => sa_1;    sa_p => sa_2
1784          ENDIF
1785
1786       CASE ( 1 )
1787          IF ( salinity )  THEN
1788             sa => sa_2;    sa_p => sa_1
1789          ENDIF
1790
1791    END SELECT
1792
1793 END SUBROUTINE ocean_swap_timelevel
1794
1795
1796!------------------------------------------------------------------------------!
1797! Description:
1798! ------------
1799!> This routine reads the respective restart data for the ocean module.
1800!------------------------------------------------------------------------------!
1801 SUBROUTINE ocean_rrd_global( found )
1802
1803
1804    USE control_parameters,                                                    &
1805        ONLY: length, restart_string
1806
1807
1808    IMPLICIT NONE
1809
1810    LOGICAL, INTENT(OUT)  ::  found
1811
1812
1813    found = .TRUE.
1814
1815    SELECT CASE ( restart_string(1:length) )
1816
1817       CASE ( 'bc_sa_t' )
1818          READ ( 13 )  bc_sa_t
1819
1820       CASE ( 'bottom_salinityflux' )
1821          READ ( 13 )  bottom_salinityflux
1822
1823       CASE ( 'salinity' )
1824          READ ( 13 )  salinity
1825
1826       CASE ( 'sa_init' )
1827          READ ( 13 )  sa_init
1828
1829       CASE ( 'sa_surface' )
1830          READ ( 13 )  sa_surface
1831
1832       CASE ( 'sa_vertical_gradient' )
1833          READ ( 13 )  sa_vertical_gradient
1834
1835       CASE ( 'sa_vertical_gradient_level' )
1836          READ ( 13 )  sa_vertical_gradient_level
1837
1838       CASE ( 'stokes_waveheight' )
1839          READ ( 13 )  stokes_waveheight
1840
1841       CASE ( 'stokes_wavelength' )
1842          READ ( 13 )  stokes_wavelength
1843
1844       CASE ( 'surface_cooling_spinup_time' )
1845          READ ( 13 )  surface_cooling_spinup_time
1846
1847       CASE ( 'top_salinityflux' )
1848          READ ( 13 )  top_salinityflux
1849
1850       CASE ( 'wall_salinityflux' )
1851          READ ( 13 )  wall_salinityflux
1852
1853       CASE ( 'wave_breaking' )
1854          READ ( 13 )  wave_breaking
1855
1856       CASE DEFAULT
1857
1858          found = .FALSE.
1859
1860    END SELECT
1861
1862 END SUBROUTINE ocean_rrd_global
1863
1864
1865!------------------------------------------------------------------------------!
1866! Description:
1867! ------------
1868!> This routine reads the respective restart data for the ocean module.
1869!------------------------------------------------------------------------------!
1870 SUBROUTINE ocean_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,           &
1871                             nxr_on_file, nynf, nync, nyn_on_file, nysf,       &
1872                             nysc, nys_on_file, tmp_3d, found )
1873
1874    USE averaging,                                                             &
1875        ONLY:  rho_ocean_av, sa_av
1876
1877    USE control_parameters,                                                    &
1878        ONLY:  length, restart_string
1879
1880    USE indices,                                                               &
1881        ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzt
1882
1883    USE pegrid
1884
1885
1886    IMPLICIT NONE
1887
1888    INTEGER(iwp) ::  k               !<
1889    INTEGER(iwp) ::  nxlc            !<
1890    INTEGER(iwp) ::  nxlf            !<
1891    INTEGER(iwp) ::  nxl_on_file     !<
1892    INTEGER(iwp) ::  nxrc            !<
1893    INTEGER(iwp) ::  nxrf            !<
1894    INTEGER(iwp) ::  nxr_on_file     !<
1895    INTEGER(iwp) ::  nync            !<
1896    INTEGER(iwp) ::  nynf            !<
1897    INTEGER(iwp) ::  nyn_on_file     !<
1898    INTEGER(iwp) ::  nysc            !<
1899    INTEGER(iwp) ::  nysf            !<
1900    INTEGER(iwp) ::  nys_on_file     !<
1901
1902    LOGICAL, INTENT(OUT)  ::  found
1903
1904    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
1905
1906
1907    found = .TRUE.
1908
1909    SELECT CASE ( restart_string(1:length) )
1910
1911       CASE ( 'rho_ocean_av' )
1912          IF ( .NOT. ALLOCATED( rho_ocean_av ) )  THEN
1913             ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1914          ENDIF
1915          IF ( k == 1 )  READ ( 13 )  tmp_3d
1916          rho_ocean_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =            &
1917                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1918
1919       CASE ( 'sa' )
1920          IF ( k == 1 )  READ ( 13 )  tmp_3d
1921          sa(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                      &
1922                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1923
1924       CASE ( 'sa_av' )
1925          IF ( .NOT. ALLOCATED( sa_av ) )  THEN
1926              ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1927          ENDIF
1928          IF ( k == 1 )  READ ( 13 )  tmp_3d
1929          sa_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                   &
1930                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1931
1932       CASE DEFAULT
1933          found = .FALSE.
1934
1935    END SELECT
1936
1937 END SUBROUTINE ocean_rrd_local
1938
1939
1940!------------------------------------------------------------------------------!
1941! Description:
1942! ------------
1943!> This routine writes the respective restart data for the ocean module.
1944!------------------------------------------------------------------------------!
1945 SUBROUTINE ocean_wrd_global
1946
1947
1948    IMPLICIT NONE
1949
1950    CALL wrd_write_string( 'bc_sa_t' )
1951    WRITE ( 14 )  bc_sa_t
1952
1953    CALL wrd_write_string( 'bottom_salinityflux' )
1954    WRITE ( 14 )  bottom_salinityflux
1955
1956    CALL wrd_write_string( 'salinity' )
1957    WRITE ( 14 )  salinity
1958
1959    CALL wrd_write_string( 'sa_init' )
1960    WRITE ( 14 )  sa_init
1961
1962    CALL wrd_write_string( 'sa_surface' )
1963    WRITE ( 14 )  sa_surface
1964
1965    CALL wrd_write_string( 'sa_vertical_gradient' )
1966    WRITE ( 14 )  sa_vertical_gradient
1967
1968    CALL wrd_write_string( 'sa_vertical_gradient_level' )
1969    WRITE ( 14 )  sa_vertical_gradient_level
1970
1971    CALL wrd_write_string( 'stokes_waveheight' )
1972    WRITE ( 14 )  stokes_waveheight
1973
1974    CALL wrd_write_string( 'stokes_wavelength' )
1975    WRITE ( 14 )  stokes_wavelength
1976
1977    CALL wrd_write_string( 'surface_cooling_spinup_time' )
1978    WRITE ( 14 )  surface_cooling_spinup_time
1979
1980    CALL wrd_write_string( 'top_salinityflux' )
1981    WRITE ( 14 )  top_salinityflux
1982
1983    CALL wrd_write_string( 'wall_salinityflux' )
1984    WRITE ( 14 )  wall_salinityflux
1985
1986    CALL wrd_write_string( 'wave_breaking' )
1987    WRITE ( 14 )  wave_breaking
1988
1989 END SUBROUTINE ocean_wrd_global
1990
1991
1992!------------------------------------------------------------------------------!
1993! Description:
1994! ------------
1995!> This routine writes the respective restart data for the ocean module.
1996!------------------------------------------------------------------------------!
1997 SUBROUTINE ocean_wrd_local
1998
1999    USE averaging,                                                             &
2000        ONLY:  rho_ocean_av, sa_av
2001
2002    IMPLICIT NONE
2003
2004    IF ( ALLOCATED( rho_ocean_av ) )  THEN
2005       CALL wrd_write_string( 'rho_ocean_av' )
2006       WRITE ( 14 )  rho_ocean_av
2007    ENDIF
2008
2009    CALL wrd_write_string( 'sa' )
2010    WRITE ( 14 )  sa
2011
2012    IF ( ALLOCATED( sa_av ) )  THEN
2013       CALL wrd_write_string( 'sa_av' )
2014       WRITE ( 14 )  sa_av
2015    ENDIF
2016
2017 END SUBROUTINE ocean_wrd_local
2018
2019
2020!------------------------------------------------------------------------------!
2021! Description:
2022! ------------
2023!> This routine calculates the Craik Leibovich vortex force and the additional
2024!> effect of the Stokes drift on the Coriolis force
2025!> Call for all gridpoints.
2026!------------------------------------------------------------------------------!
2027 SUBROUTINE stokes_drift_terms( component )
2028
2029    USE arrays_3d,                                                             &
2030        ONLY:  ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu,              &
2031               v_stokes_zw, w, tend
2032
2033    USE control_parameters,                                                    &
2034        ONLY:  f, fs, message_string
2035
2036    USE grid_variables,                                                        &
2037        ONLY:  ddx, ddy
2038
2039    USE indices,                                                               &
2040        ONLY:  nxl, nxr, nys, nysv, nyn, nzb, nzt
2041
2042    IMPLICIT NONE
2043
2044    INTEGER(iwp) ::  component  !< component of momentum equation
2045    INTEGER(iwp) ::  i          !< loop index along x
2046    INTEGER(iwp) ::  j          !< loop index along y
2047    INTEGER(iwp) ::  k          !< loop index along z
2048
2049
2050!
2051!-- Compute Stokes terms for the respective velocity components
2052    SELECT CASE ( component )
2053
2054!
2055!--    u-component
2056       CASE ( 1 )
2057          DO  i = nxl, nxr
2058             DO  j = nysv, nyn
2059                DO  k = nzb+1, nzt
2060                   tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * (              &
2061                                   0.5 * ( v(k,j+1,i) - v(k,j+1,i-1)           &
2062                                         + v(k,j,i)   - v(k,j,i-1)   ) * ddx   &
2063                                 - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) )   * ddy   &
2064                                                                )              &
2065                                 + f * v_stokes_zu(k)
2066                ENDDO
2067             ENDDO
2068          ENDDO
2069
2070!
2071!--    v-component
2072       CASE ( 2 )
2073          DO  i = nxl, nxr
2074             DO  j = nysv, nyn
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             ENDDO
2084          ENDDO
2085
2086!
2087!--    w-component
2088       CASE ( 3 )
2089          DO  i = nxl, nxr
2090             DO  j = nys, nyn
2091                DO  k = nzb+1, nzt
2092                   tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * (              &
2093                                             0.5 * ( u(k+1,j,i) - u(k,j,i)     &
2094                                                   + u(k+1,j,i+1) - u(k,j,i+1) &
2095                                                   ) * ddzu(k+1)               &
2096                                           - 0.5 * ( w(k,j,i+1) - w(k,j,i-1)   &
2097                                                   ) * ddx      )              &
2098                                             - v_stokes_zw(k) * (              &
2099                                             0.5 * ( w(k,j+1,i) - w(k,j-1,i)   &
2100                                                   ) * ddy                     &
2101                                           - 0.5 * ( v(k+1,j,i) - v(k,j,i)     &
2102                                                   + v(k+1,j+1,i) - v(k,j+1,i) &
2103                                                   ) * ddzu(k)  )              &
2104                                           + fs * u_stokes_zw(k)
2105                ENDDO
2106             ENDDO
2107          ENDDO
2108
2109       CASE DEFAULT
2110          WRITE( message_string, * ) 'wrong component of Stokes force: ',      &
2111                                     component
2112          CALL message( 'stokes_drift_terms', 'PA0091', 1, 2, 0, 6, 0 )
2113
2114    END SELECT
2115
2116 END SUBROUTINE stokes_drift_terms
2117
2118
2119!------------------------------------------------------------------------------!
2120! Description:
2121! ------------
2122!> This routine calculates the Craik Leibovich vortex force and the additional
2123!> effect of the Stokes drift on the Coriolis force
2124!> Call for gridpoints i,j.
2125!------------------------------------------------------------------------------!
2126
2127 SUBROUTINE stokes_drift_terms_ij( i, j, component )
2128
2129    USE arrays_3d,                                                             &
2130        ONLY:  ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu,              &
2131               v_stokes_zw, w, tend
2132
2133    USE control_parameters,                                                    &
2134        ONLY:  f, fs, message_string
2135
2136    USE grid_variables,                                                        &
2137        ONLY:  ddx, ddy
2138
2139    USE indices,                                                               &
2140        ONLY:  nzb, nzt
2141
2142    IMPLICIT NONE
2143
2144    INTEGER(iwp) ::  component  !< component of momentum equation
2145    INTEGER(iwp) ::  i          !< loop index along x
2146    INTEGER(iwp) ::  j          !< loop index along y
2147    INTEGER(iwp) ::  k          !< loop incex along z
2148
2149
2150!
2151!-- Compute Stokes terms for the respective velocity components
2152    SELECT CASE ( component )
2153
2154!
2155!--    u-component
2156       CASE ( 1 )
2157          DO  k = nzb+1, nzt
2158             tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * (                    &
2159                                     0.5 * ( v(k,j+1,i) - v(k,j+1,i-1)         &
2160                                           + v(k,j,i)   - v(k,j,i-1)   ) * ddx &
2161                                   - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) )   * ddy &
2162                                                          )                    &
2163                                       + f * v_stokes_zu(k)
2164          ENDDO
2165!
2166!--    v-component
2167       CASE ( 2 )
2168          DO  k = nzb+1, nzt
2169             tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * (                    &
2170                                     0.5 * ( v(k,j,i+1) - v(k,j,i-1) )   * ddx &
2171                                   - 0.5 * ( u(k,j,i) - u(k,j-1,i)             &
2172                                           + u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
2173                                                          )                    &
2174                                       - f * u_stokes_zu(k)
2175          ENDDO
2176
2177!
2178!--    w-component
2179       CASE ( 3 )
2180          DO  k = nzb+1, nzt
2181             tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * (              &
2182                                     0.5 * ( u(k+1,j,i) - u(k,j,i)     &
2183                                                   + u(k+1,j,i+1) - u(k,j,i+1) &
2184                                                   ) * ddzu(k+1)               &
2185                                           - 0.5 * ( w(k,j,i+1) - w(k,j,i-1)   &
2186                                                   ) * ddx )                   &
2187                                       - v_stokes_zw(k) * (                    &
2188                                             0.5 * ( w(k,j+1,i) - w(k,j-1,i)   &
2189                                                   ) * ddy                     &
2190                                           - 0.5 * ( v(k+1,j,i) - v(k,j,i)     &
2191                                                   + v(k+1,j+1,i) - v(k,j+1,i) &
2192                                                   ) * ddzu(k)  )              &
2193                                       + fs * u_stokes_zw(k)
2194          ENDDO
2195
2196       CASE DEFAULT
2197          WRITE( message_string, * ) ' wrong component: ', component
2198          CALL message( 'stokes_drift_terms', 'PA0091', 1, 2, 0, 6, 0 )
2199
2200    END SELECT
2201
2202 END SUBROUTINE stokes_drift_terms_ij
2203
2204
2205!------------------------------------------------------------------------------!
2206! Description:
2207! ------------
2208!> This routine calculates turbulence generated by wave breaking near the ocean
2209!> surface, following a parameterization given in Noh et al. (2004), JPO
2210!> Call for all gridpoints.
2211!> TODO: so far, this routine only works if the model time step has about the
2212!>       same value as the time scale of wave breaking!
2213!------------------------------------------------------------------------------!
2214 SUBROUTINE wave_breaking_term( component )
2215
2216    USE arrays_3d,                                                             &
2217        ONLY:  u_p, v_p
2218
2219    USE control_parameters,                                                    &
2220        ONLY:  dt_3d, message_string
2221
2222    USE indices,                                                               &
2223        ONLY:  nxl, nxlu, nxr, nys, nysv, nyn, nzt
2224
2225    IMPLICIT NONE
2226
2227    INTEGER(iwp) ::  component  !< component of momentum equation
2228    INTEGER(iwp) ::  i          !< loop index along x
2229    INTEGER(iwp) ::  j          !< loop index along y
2230
2231    REAL(wp) ::  random_gauss  !< function that creates a random number with a
2232                               !< Gaussian distribution
2233
2234
2235!
2236!-- Compute wave breaking terms for the respective velocity components.
2237!-- Velocities are directly manipulated, since this is not a real force
2238    SELECT CASE ( component )
2239
2240!
2241!--    u-component
2242       CASE ( 1 )
2243          DO  i = nxlu, nxr
2244             DO  j = nys, nyn
2245                u_p(nzt,j,i) = u_p(nzt,j,i) +                                  &
2246                               ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &
2247                               * alpha_wave_breaking * u_star_wave_breaking    &
2248                               / timescale_wave_breaking * dt_3d
2249             ENDDO
2250          ENDDO
2251!
2252!--    v-component
2253       CASE ( 2 )
2254          DO  i = nxl, nxr
2255             DO  j = nysv, nyn
2256                v_p(nzt,j,i) = v_p(nzt,j,i) +                                  &
2257                               ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &
2258                               * alpha_wave_breaking * u_star_wave_breaking    &
2259                               / timescale_wave_breaking * dt_3d
2260             ENDDO
2261          ENDDO
2262
2263       CASE DEFAULT
2264          WRITE( message_string, * ) 'wrong component of wave breaking: ',     &
2265                                     component
2266          CALL message( 'stokes_drift_terms', 'PA0466', 1, 2, 0, 6, 0 )
2267
2268    END SELECT
2269
2270 END SUBROUTINE wave_breaking_term
2271
2272
2273!------------------------------------------------------------------------------!
2274! Description:
2275! ------------
2276!> This routine calculates turbulence generated by wave breaking near the ocean
2277!> surface, following a parameterization given in Noh et al. (2004), JPO
2278!> Call for gridpoint i,j.
2279!> TODO: so far, this routine only works if the model time step has about the
2280!>       same value as the time scale of wave breaking!
2281!------------------------------------------------------------------------------!
2282 SUBROUTINE wave_breaking_term_ij( i, j, component )
2283
2284    USE arrays_3d,                                                             &
2285        ONLY:  u_p, v_p
2286
2287    USE control_parameters,                                                    &
2288        ONLY:  dt_3d, message_string
2289
2290    USE indices,                                                               &
2291        ONLY:  nzt
2292
2293    IMPLICIT NONE
2294
2295    INTEGER(iwp) ::  component  !< component of momentum equation
2296    INTEGER(iwp) ::  i          !< loop index along x
2297    INTEGER(iwp) ::  j          !< loop index along y
2298
2299    REAL(wp) ::  random_gauss  !< function that creates a random number with a
2300                               !< Gaussian distribution
2301
2302!
2303!-- Compute wave breaking terms for the respective velocity components
2304    SELECT CASE ( component )
2305
2306!
2307!--    u-/v-component
2308       CASE ( 1 )
2309          u_p(nzt,j,i) = u_p(nzt,j,i) +                                        &
2310                         ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp )       &
2311                         * alpha_wave_breaking * u_star_wave_breaking          &
2312                         / timescale_wave_breaking * dt_3d
2313
2314       CASE ( 2 )
2315          v_p(nzt,j,i) = v_p(nzt,j,i) +                                        &
2316                         ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp )       &
2317                         * alpha_wave_breaking * u_star_wave_breaking          &
2318                         / timescale_wave_breaking * dt_3d
2319
2320       CASE DEFAULT
2321          WRITE( message_string, * ) 'wrong component of wave breaking: ',     &
2322                                     component
2323          CALL message( 'stokes_drift_terms', 'PA0466', 1, 2, 0, 6, 0 )
2324
2325    END SELECT
2326
2327 END SUBROUTINE wave_breaking_term_ij
2328
2329
2330 END MODULE ocean_mod
Note: See TracBrowser for help on using the repository browser.