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

Last change on this file since 3381 was 3381, checked in by raasch, 3 years ago

spin-up cooling for ocean surface implemented, see new parameter surface_cooling_spinup_time

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