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

Last change on this file since 3656 was 3636, checked in by raasch, 6 years ago

nopointer option removed

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