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

Last change on this file since 4808 was 4797, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

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