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

Last change on this file since 3529 was 3421, checked in by gronemeier, 6 years ago

new surface-data output; renamed output variables (pt to theta, rho_air to rho, rho_ocean to rho_sea_water)

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