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

Last change on this file since 3593 was 3568, checked in by raasch, 6 years ago

bugifx: calculate equation of state for seawater even if salinity is switched off

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