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

Last change on this file since 3873 was 3873, checked in by knoop, 2 years ago

Moved ocean_mode specific code from advec_ws to ocean_mod + implemented ocean_actions

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