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

Last change on this file since 3721 was 3719, checked in by kanani, 6 years ago

Correct and clean-up cpu_logs, some overlapping counts (chemistry_model_mod, disturb_heatflux, large_scale_forcing_nudging_mod, ocean_mod, palm, prognostic_equations, synthetic_turbulence_generator_mod, time_integration, time_integration_spinup, turbulence_closure_mod)

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