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

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

modularization of the ocean code

  • Property svn:keywords set to Id
File size: 60.4 KB
Line 
1!> @file ocean_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2017-2018 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: ocean_mod.f90 3294 2018-10-01 02:37:10Z raasch $
27! initial revision
28!
29!
30!
31!
32! Authors:
33! --------
34! @author Siegfried Raasch
35!
36! Description:
37! ------------
38!> This module contains relevant code for PALM's ocean option, e.g. the
39!> prognostic equation for salinity, the equation of state for seawater,
40!> and the Craik Leibovich force (Stokes force)
41!>
42!------------------------------------------------------------------------------!
43 MODULE ocean_mod
44 
45
46#if defined( __nopointer )
47    USE arrays_3d,                                                             &
48        ONLY:  prho, rho_ocean, sa, sa_init, sa_p, tsa_m
49#else
50    USE arrays_3d,                                                             &
51        ONLY:  prho, prho_1, rho_ocean, rho_1, sa, sa_init, sa_1, sa_2, sa_3,  &
52               sa_p, tsa_m
53#endif
54
55    USE control_parameters,                                                    &
56        ONLY:  atmos_ocean_sign, bottom_salinityflux,                          &
57               constant_top_salinityflux, ocean_mode, top_salinityflux,        &
58               wall_salinityflux
59
60    USE kinds
61
62
63    IMPLICIT NONE
64
65    CHARACTER (LEN=20) ::  bc_sa_t = 'neumann'  !< namelist parameter
66
67    INTEGER(iwp) ::  ibc_sa_t   !< integer flag for bc_sa_t
68
69    INTEGER(iwp) ::  sa_vertical_gradient_level_ind(10) = -9999  !< grid index values of sa_vertical_gradient_level(s)
70
71    LOGICAL ::  stokes_force = .FALSE.              !< switch to switch on the Stokes force
72
73    REAL(wp) ::  langmuir_number = 0.34_wp    !< turbulent Langmuir number, default from Li et al., 2005
74    REAL(wp) ::  prho_reference               !< reference state of potential density at ocean surface
75    REAL(wp) ::  sa_surface = 35.0_wp         !< surface salinity, namelist parameter
76    REAL(wp) ::  sa_vertical_gradient(10) = 0.0_wp               !< namelist parameter
77    REAL(wp) ::  sa_vertical_gradient_level(10) = -999999.9_wp   !< namelist parameter
78    REAL(wp) ::  stokes_waveheight =  1.0_wp  !< typical wave height assumed for Stokes velocity
79    REAL(wp) ::  stokes_wavelength = 40.0_wp  !< typical wavelength assumed for Stokes velocity
80
81    REAL(wp), DIMENSION(12), PARAMETER ::  nom =                               &
82                          (/ 9.9984085444849347D2,   7.3471625860981584D0,     &
83                            -5.3211231792841769D-2,  3.6492439109814549D-4,    &
84                             2.5880571023991390D0,  -6.7168282786692354D-3,    &
85                             1.9203202055760151D-3,  1.1798263740430364D-2,    &
86                             9.8920219266399117D-8,  4.6996642771754730D-6,    &
87                            -2.5862187075154352D-8, -3.2921414007960662D-12 /)
88                          !< constants used in equation of state for seawater
89
90    REAL(wp), DIMENSION(13), PARAMETER ::  den =                               &
91                          (/ 1.0D0,                  7.2815210113327091D-3,    &
92                            -4.4787265461983921D-5,  3.3851002965802430D-7,    &
93                             1.3651202389758572D-10, 1.7632126669040377D-3,    &
94                            -8.8066583251206474D-6, -1.8832689434804897D-10,   &
95                             5.7463776745432097D-6,  1.4716275472242334D-9,    &
96                             6.7103246285651894D-6, -2.4461698007024582D-17,   &
97                            -9.1534417604289062D-18 /)
98                          !< constants used in equation of state for seawater
99
100    SAVE
101
102    PUBLIC ::  bc_sa_t, ibc_sa_t, langmuir_number, prho_reference, sa_surface, &
103               sa_vertical_gradient, sa_vertical_gradient_level,               &
104               sa_vertical_gradient_level_ind, stokes_force,                   &
105               stokes_waveheight, stokes_wavelength
106
107
108    INTERFACE eqn_state_seawater
109       MODULE PROCEDURE eqn_state_seawater
110       MODULE PROCEDURE eqn_state_seawater_ij
111    END INTERFACE eqn_state_seawater
112
113    INTERFACE eqn_state_seawater_func
114       MODULE PROCEDURE eqn_state_seawater_func
115    END INTERFACE eqn_state_seawater_func
116
117    INTERFACE ocean_check_parameters
118       MODULE PROCEDURE ocean_check_parameters
119    END INTERFACE ocean_check_parameters
120
121    INTERFACE ocean_check_data_output
122       MODULE PROCEDURE ocean_check_data_output
123    END INTERFACE ocean_check_data_output
124
125    INTERFACE ocean_check_data_output_pr
126       MODULE PROCEDURE ocean_check_data_output_pr
127    END INTERFACE ocean_check_data_output_pr
128
129    INTERFACE ocean_define_netcdf_grid
130       MODULE PROCEDURE ocean_define_netcdf_grid
131    END INTERFACE ocean_define_netcdf_grid
132
133    INTERFACE ocean_data_output_2d
134       MODULE PROCEDURE ocean_data_output_2d
135    END INTERFACE ocean_data_output_2d
136
137    INTERFACE ocean_data_output_3d
138       MODULE PROCEDURE ocean_data_output_3d
139    END INTERFACE ocean_data_output_3d
140
141    INTERFACE ocean_init
142       MODULE PROCEDURE ocean_init
143    END INTERFACE ocean_init
144
145    INTERFACE ocean_init_arrays
146       MODULE PROCEDURE ocean_init_arrays
147    END INTERFACE ocean_init_arrays
148
149    INTERFACE ocean_parin
150       MODULE PROCEDURE ocean_parin
151    END INTERFACE ocean_parin
152
153    INTERFACE ocean_prognostic_equations
154       MODULE PROCEDURE ocean_prognostic_equations
155       MODULE PROCEDURE ocean_prognostic_equations_ij
156    END INTERFACE ocean_prognostic_equations
157
158    INTERFACE ocean_swap_timelevel
159       MODULE PROCEDURE ocean_swap_timelevel
160    END INTERFACE ocean_swap_timelevel
161
162    INTERFACE ocean_rrd_global
163       MODULE PROCEDURE ocean_rrd_global
164    END INTERFACE ocean_rrd_global
165
166    INTERFACE ocean_rrd_local
167       MODULE PROCEDURE ocean_rrd_local
168    END INTERFACE ocean_rrd_local
169
170    INTERFACE ocean_wrd_global
171       MODULE PROCEDURE ocean_wrd_global
172    END INTERFACE ocean_wrd_global
173
174    INTERFACE ocean_wrd_local
175       MODULE PROCEDURE ocean_wrd_local
176    END INTERFACE ocean_wrd_local
177
178    INTERFACE ocean_3d_data_averaging
179       MODULE PROCEDURE ocean_3d_data_averaging
180    END INTERFACE ocean_3d_data_averaging
181
182    SAVE
183
184    PRIVATE
185!
186!-- Add INTERFACES that must be available to other modules (alphabetical order)
187    PUBLIC eqn_state_seawater, ocean_check_data_output,                        &
188           ocean_check_data_output_pr, ocean_check_parameters,                 &
189           ocean_data_output_2d, ocean_data_output_3d,                         &
190           ocean_define_netcdf_grid, ocean_init, ocean_init_arrays,            &
191           ocean_parin, ocean_prognostic_equations, ocean_swap_timelevel,      &
192           ocean_rrd_global, ocean_rrd_local, ocean_wrd_global,                &
193           ocean_wrd_local, ocean_3d_data_averaging
194
195
196 CONTAINS
197
198!------------------------------------------------------------------------------!
199! Description:
200! ------------
201!> Equation of state for seawater as a function of potential temperature,
202!> salinity, and pressure.
203!> For coefficients see Jackett et al., 2006: J. Atm. Ocean Tech.
204!> eqn_state_seawater calculates the potential density referred at hyp(0).
205!> eqn_state_seawater_func calculates density.
206!> TODO: so far, routine is not adjusted to use topography
207!------------------------------------------------------------------------------!
208 SUBROUTINE eqn_state_seawater
209
210    USE arrays_3d,                                                             &
211        ONLY:  hyp, prho, pt_p, rho_ocean, sa_p
212    USE indices,                                                               &
213        ONLY:  nxl, nxr, nyn, nys, nzb, nzt
214
215    USE surface_mod,                                                           &
216        ONLY :  bc_h
217
218    IMPLICIT NONE
219
220    INTEGER(iwp) ::  i       !< running index x direction
221    INTEGER(iwp) ::  j       !< running index y direction
222    INTEGER(iwp) ::  k       !< running index z direction
223    INTEGER(iwp) ::  m       !< running index surface elements
224
225    REAL(wp) ::  pden   !< temporary scalar user for calculating density
226    REAL(wp) ::  pnom   !< temporary scalar user for calculating density
227    REAL(wp) ::  p1     !< temporary scalar user for calculating density
228    REAL(wp) ::  p2     !< temporary scalar user for calculating density
229    REAL(wp) ::  p3     !< temporary scalar user for calculating density
230    REAL(wp) ::  pt1    !< temporary scalar user for calculating density
231    REAL(wp) ::  pt2    !< temporary scalar user for calculating density
232    REAL(wp) ::  pt3    !< temporary scalar user for calculating density
233    REAL(wp) ::  pt4    !< temporary scalar user for calculating density
234    REAL(wp) ::  sa1    !< temporary scalar user for calculating density
235    REAL(wp) ::  sa15   !< temporary scalar user for calculating density
236    REAL(wp) ::  sa2    !< temporary scalar user for calculating density
237
238
239    DO  i = nxl, nxr
240       DO  j = nys, nyn
241          DO  k = nzb+1, nzt
242!
243!--          Pressure is needed in dbar
244             p1 = hyp(k) * 1E-4_wp
245             p2 = p1 * p1
246             p3 = p2 * p1
247
248!
249!--          Temperature needed in degree Celsius
250             pt1 = pt_p(k,j,i) - 273.15_wp
251             pt2 = pt1 * pt1
252             pt3 = pt1 * pt2
253             pt4 = pt2 * pt2
254
255             sa1  = sa_p(k,j,i)
256             sa15 = sa1 * SQRT( sa1 )
257             sa2  = sa1 * sa1
258
259             pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     +       &
260                    nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 +       &
261                    nom(7)*sa2
262
263             pden = den(1)           + den(2)*pt1     + den(3)*pt2     +       &
264                    den(4)*pt3       + den(5)*pt4     + den(6)*sa1     +       &
265                    den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    +       &
266                    den(10)*sa15*pt2
267!
268!--          Potential density (without pressure terms)
269             prho(k,j,i) = pnom / pden
270
271             pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  +       &
272                    nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
273
274             pden = pden +             den(11)*p1     + den(12)*p2*pt3 +       &
275                    den(13)*p3*pt1
276
277!
278!--          In-situ density
279             rho_ocean(k,j,i) = pnom / pden
280
281          ENDDO
282
283!
284!--       Neumann conditions are assumed at top boundary and currently also at
285!--       bottom boundary (see comment lines below)
286          prho(nzt+1,j,i)      = prho(nzt,j,i)
287          rho_ocean(nzt+1,j,i) = rho_ocean(nzt,j,i)
288
289       ENDDO
290    ENDDO
291!
292!-- Neumann conditions at up/downward-facing surfaces
293    !$OMP PARALLEL DO PRIVATE( i, j, k )
294    DO  m = 1, bc_h(0)%ns
295       i = bc_h(0)%i(m)
296       j = bc_h(0)%j(m)
297       k = bc_h(0)%k(m)
298       prho(k-1,j,i)      = prho(k,j,i)
299       rho_ocean(k-1,j,i) = rho_ocean(k,j,i)
300    ENDDO
301!
302!-- Downward facing surfaces
303    !$OMP PARALLEL DO PRIVATE( i, j, k )
304    DO  m = 1, bc_h(1)%ns
305       i = bc_h(1)%i(m)
306       j = bc_h(1)%j(m)
307       k = bc_h(1)%k(m)
308       prho(k+1,j,i)      = prho(k,j,i)
309       rho_ocean(k+1,j,i) = rho_ocean(k,j,i)
310    ENDDO
311
312 END SUBROUTINE eqn_state_seawater
313
314
315!------------------------------------------------------------------------------!
316! Description:
317! ------------
318!> Same as above, but for grid point i,j
319!------------------------------------------------------------------------------!
320 SUBROUTINE eqn_state_seawater_ij( i, j )
321
322    USE arrays_3d,                                                             &
323        ONLY:  hyp, prho, pt_p, rho_ocean, sa_p
324
325    USE indices,                                                               &
326        ONLY:  nzb, nzt
327
328    USE surface_mod,                                                           &
329       ONLY :  bc_h
330
331    IMPLICIT NONE
332
333    INTEGER(iwp) ::  i       !< running index x direction
334    INTEGER(iwp) ::  j       !< running index y direction
335    INTEGER(iwp) ::  k       !< running index z direction
336    INTEGER(iwp) ::  m       !< running index surface elements
337    INTEGER(iwp) ::  surf_e  !< end index of surface elements at (j,i)-gridpoint
338    INTEGER(iwp) ::  surf_s  !< start index of surface elements at (j,i)-gridpoint
339
340    REAL(wp) ::  pden   !< temporary scalar user for calculating density
341    REAL(wp) ::  pnom   !< temporary scalar user for calculating density
342    REAL(wp) ::  p1     !< temporary scalar user for calculating density
343    REAL(wp) ::  p2     !< temporary scalar user for calculating density
344    REAL(wp) ::  p3     !< temporary scalar user for calculating density
345    REAL(wp) ::  pt1    !< temporary scalar user for calculating density
346    REAL(wp) ::  pt2    !< temporary scalar user for calculating density
347    REAL(wp) ::  pt3    !< temporary scalar user for calculating density
348    REAL(wp) ::  pt4    !< temporary scalar user for calculating density
349    REAL(wp) ::  sa1    !< temporary scalar user for calculating density
350    REAL(wp) ::  sa15   !< temporary scalar user for calculating density
351    REAL(wp) ::  sa2    !< temporary scalar user for calculating density
352
353    DO  k = nzb+1, nzt
354!
355!--    Pressure is needed in dbar
356       p1 = hyp(k) * 1E-4_wp
357       p2 = p1 * p1
358       p3 = p2 * p1
359!
360!--    Temperature needed in degree Celsius
361       pt1 = pt_p(k,j,i) - 273.15_wp
362       pt2 = pt1 * pt1
363       pt3 = pt1 * pt2
364       pt4 = pt2 * pt2
365
366       sa1  = sa_p(k,j,i)
367       sa15 = sa1 * SQRT( sa1 )
368       sa2  = sa1 * sa1
369
370       pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     +             &
371              nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + nom(7)*sa2
372
373       pden = den(1)           + den(2)*pt1     + den(3)*pt2     +             &
374              den(4)*pt3       + den(5)*pt4     + den(6)*sa1     +             &
375              den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    +             &
376              den(10)*sa15*pt2
377!
378!--    Potential density (without pressure terms)
379       prho(k,j,i) = pnom / pden
380
381       pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  +             &
382              nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
383       pden = pden +             den(11)*p1     + den(12)*p2*pt3 +             &
384              den(13)*p3*pt1
385
386!
387!--    In-situ density
388       rho_ocean(k,j,i) = pnom / pden
389
390    ENDDO
391!
392!-- Neumann conditions at up/downward-facing walls
393    surf_s = bc_h(0)%start_index(j,i)
394    surf_e = bc_h(0)%end_index(j,i)
395    DO  m = surf_s, surf_e
396       k = bc_h(0)%k(m)
397       prho(k-1,j,i)      = prho(k,j,i)
398       rho_ocean(k-1,j,i) = rho_ocean(k,j,i)
399    ENDDO
400!
401!-- Downward facing surfaces
402    surf_s = bc_h(1)%start_index(j,i)
403    surf_e = bc_h(1)%end_index(j,i)
404    DO  m = surf_s, surf_e
405       k = bc_h(1)%k(m)
406       prho(k+1,j,i)      = prho(k,j,i)
407       rho_ocean(k+1,j,i) = rho_ocean(k,j,i)
408    ENDDO
409!
410!-- Neumann condition are assumed at top boundary
411    prho(nzt+1,j,i)      = prho(nzt,j,i)
412    rho_ocean(nzt+1,j,i) = rho_ocean(nzt,j,i)
413
414 END SUBROUTINE eqn_state_seawater_ij
415
416
417!------------------------------------------------------------------------------!
418! Description:
419! ------------
420!> Equation of state for seawater as a function
421!------------------------------------------------------------------------------!
422 REAL(wp) FUNCTION eqn_state_seawater_func( p, pt, sa )
423
424    IMPLICIT NONE
425
426    REAL(wp) ::  p      !< temporary scalar user for calculating density
427    REAL(wp) ::  p1     !< temporary scalar user for calculating density
428    REAL(wp) ::  p2     !< temporary scalar user for calculating density
429    REAL(wp) ::  p3     !< temporary scalar user for calculating density
430    REAL(wp) ::  pt     !< temporary scalar user for calculating density
431    REAL(wp) ::  pt1    !< temporary scalar user for calculating density
432    REAL(wp) ::  pt2    !< temporary scalar user for calculating density
433    REAL(wp) ::  pt3    !< temporary scalar user for calculating density
434    REAL(wp) ::  pt4    !< temporary scalar user for calculating density
435    REAL(wp) ::  sa     !< temporary scalar user for calculating density
436    REAL(wp) ::  sa15   !< temporary scalar user for calculating density
437    REAL(wp) ::  sa2    !< temporary scalar user for calculating density
438
439!
440!-- Pressure is needed in dbar
441    p1 = p  * 1.0E-4_wp
442    p2 = p1 * p1
443    p3 = p2 * p1
444
445!
446!-- Temperature needed in degree Celsius
447    pt1 = pt - 273.15_wp
448    pt2 = pt1 * pt1
449    pt3 = pt1 * pt2
450    pt4 = pt2 * pt2
451
452    sa15 = sa * SQRT( sa )
453    sa2  = sa * sa
454
455
456    eqn_state_seawater_func =                                                  &
457         ( nom(1)        + nom(2)*pt1       + nom(3)*pt2    + nom(4)*pt3     + &
458           nom(5)*sa     + nom(6)*sa*pt1    + nom(7)*sa2    + nom(8)*p1      + &
459           nom(9)*p1*pt2 + nom(10)*p1*sa    + nom(11)*p2    + nom(12)*p2*pt2   &
460         ) /                                                                   &
461         ( den(1)        + den(2)*pt1       + den(3)*pt2    + den(4)*pt3     + &
462           den(5)*pt4    + den(6)*sa        + den(7)*sa*pt1 + den(8)*sa*pt3  + &
463           den(9)*sa15   + den(10)*sa15*pt2 + den(11)*p1    + den(12)*p2*pt3 + &
464           den(13)*p3*pt1                                                      &
465         )
466
467
468 END FUNCTION eqn_state_seawater_func
469
470
471!------------------------------------------------------------------------------!
472! Description:
473! ------------
474!> Reads the ocean parameters namelist
475!------------------------------------------------------------------------------!
476 SUBROUTINE ocean_parin
477
478    USE control_parameters,                                                    &
479        ONLY:  message_string
480
481    IMPLICIT NONE
482
483    CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
484
485
486    NAMELIST /ocean_parameters/  bc_sa_t, bottom_salinityflux,                 &
487             langmuir_number, sa_surface, sa_vertical_gradient,                &
488             sa_vertical_gradient_level, stokes_force, stokes_waveheight,      &
489             stokes_wavelength, top_salinityflux, wall_salinityflux
490
491!
492!-- Try to find the namelist
493    REWIND ( 11 )
494    line = ' '
495    DO WHILE ( INDEX( line, '&ocean_parameters' ) == 0 )
496       READ ( 11, '(A)', END=12 )  line
497    ENDDO
498    BACKSPACE ( 11 )
499
500!
501!-- Read namelist
502    READ ( 11, ocean_parameters, ERR = 10 )
503!
504!-- Set switch that enables PALM's ocean mode
505    ocean_mode = .TRUE.
506
507    GOTO 12
508
509 10 BACKSPACE( 11 )
510    READ( 11 , '(A)') line
511    CALL parin_fail_message( 'ocean_parameters', line )
512
513 12 CONTINUE
514
515 END SUBROUTINE ocean_parin
516
517!------------------------------------------------------------------------------!
518! Description:
519! ------------
520!> Check parameters routine for the ocean mode
521!------------------------------------------------------------------------------!
522 SUBROUTINE ocean_check_parameters
523
524    USE control_parameters,                                                    &
525        ONLY:  coupling_char, coupling_mode, message_string, use_top_fluxes
526
527    IMPLICIT NONE
528
529
530!
531!-- Check ocean setting
532    IF ( TRIM( coupling_mode ) == 'uncoupled'  .AND.                           &
533         TRIM( coupling_char ) == '_O' .AND.                                   &
534         .NOT. ocean_mode )  THEN
535
536!
537!--    Check whether an (uncoupled) atmospheric run has been declared as an
538!--    ocean run (this setting is done via palmrun-option -y)
539       message_string = 'ocean = .F. does not allow coupling_char = "' //      &
540                        TRIM( coupling_char ) // '" set by palmrun-option "-y"'
541       CALL message( 'check_parameters', 'PA0317', 1, 2, 0, 6, 0 )
542
543    ENDIF
544
545!
546!-- Ocean version must use flux boundary conditions at the top
547    IF ( .NOT. use_top_fluxes )  THEN
548       message_string = 'use_top_fluxes must be .TRUE. in ocean mode'
549       CALL message( 'ocean_check_parameters', 'PA0042', 1, 2, 0, 6, 0 )
550    ENDIF
551
552!
553!-- Boundary conditions for salinity
554    IF ( bc_sa_t == 'dirichlet' )  THEN
555       ibc_sa_t = 0
556    ELSEIF ( bc_sa_t == 'neumann' )  THEN
557       ibc_sa_t = 1
558    ELSE
559       message_string = 'unknown boundary condition: bc_sa_t = "' //           &
560                        TRIM( bc_sa_t ) // '"'
561       CALL message( 'ocean_check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
562    ENDIF
563
564    IF ( top_salinityflux == 9999999.9_wp )  constant_top_salinityflux = .FALSE.
565
566    IF ( ibc_sa_t == 1  .AND.  top_salinityflux == 9999999.9_wp )  THEN
567       message_string = 'boundary condition: bc_sa_t = "' //                   &
568                        TRIM( bc_sa_t ) // '" requires to set top_salinityflux'
569       CALL message( 'ocean_check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
570    ENDIF
571
572!
573!-- A fixed salinity at the top implies Dirichlet boundary condition for
574!-- salinity. In this case specification of a constant salinity flux is
575!-- forbidden.
576    IF ( ibc_sa_t == 0  .AND.  constant_top_salinityflux  .AND.                &
577         top_salinityflux /= 0.0_wp )  THEN
578       message_string = 'boundary condition: bc_sa_t = "' //                   &
579                        TRIM( bc_sa_t ) // '" is not allowed with ' //         &
580                        'top_salinityflux /= 0.0'
581       CALL message( 'ocean_check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
582    ENDIF
583
584 END SUBROUTINE ocean_check_parameters
585
586
587!------------------------------------------------------------------------------!
588! Description:
589! ------------
590!> Check data output.
591!------------------------------------------------------------------------------!
592 SUBROUTINE ocean_check_data_output( var, unit )
593 
594    IMPLICIT NONE
595
596    CHARACTER (LEN=*) ::  unit     !< unit of output variable
597    CHARACTER (LEN=*) ::  var      !< name of output variable
598
599
600    SELECT CASE ( TRIM( var ) )
601
602       CASE ( 'rho_ocean' )
603          unit = 'kg/m3'
604
605       CASE ( 'sa' )
606          unit = 'psu'
607
608       CASE DEFAULT
609          unit = 'illegal'
610
611    END SELECT
612
613 END SUBROUTINE ocean_check_data_output
614
615
616!------------------------------------------------------------------------------!
617! Description:
618! ------------
619!> Check data output of profiles
620!------------------------------------------------------------------------------!
621 SUBROUTINE ocean_check_data_output_pr( variable, var_count, unit, dopr_unit )
622
623    USE arrays_3d,                                                             &
624        ONLY:  zu, zw
625
626    USE control_parameters,                                                    &
627        ONLY:  data_output_pr, message_string
628
629    USE indices
630
631    USE profil_parameter
632
633    USE statistics
634
635    IMPLICIT NONE
636
637    CHARACTER (LEN=*) ::  unit      !<
638    CHARACTER (LEN=*) ::  variable  !<
639    CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
640
641    INTEGER(iwp) ::  var_count     !<
642
643    SELECT CASE ( TRIM( variable ) )
644
645       CASE ( 'prho' )
646          dopr_index(var_count) = 71
647          dopr_unit             = 'kg/m3'
648          hom(:,2,71,:) = SPREAD( zu, 2, statistic_regions+1 )
649          unit = dopr_unit
650
651       CASE ( 'rho_ocean' )
652          dopr_index(var_count) = 64
653          dopr_unit             = 'kg/m3'
654          hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
655          IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
656             dopr_initial_index(var_count) = 77
657             hom(:,2,77,:)             = SPREAD( zu, 2, statistic_regions+1 )
658             hom(nzb,2,77,:)           = 0.0_wp    ! because zu(nzb) is negative
659             data_output_pr(var_count) = data_output_pr(var_count)(2:)
660          ENDIF
661          unit = dopr_unit
662
663       CASE ( 'sa', '#sa' )
664          dopr_index(var_count) = 23
665          dopr_unit             = 'psu'
666          hom(:,2,23,:) = SPREAD( zu, 2, statistic_regions+1 )
667          IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
668             dopr_initial_index(var_count) = 26
669             hom(:,2,26,:)             = SPREAD( zu, 2, statistic_regions+1 )
670             hom(nzb,2,26,:)           = 0.0_wp    ! because zu(nzb) is negative
671             data_output_pr(var_count) = data_output_pr(var_count)(2:)
672          ENDIF
673          unit = dopr_unit
674
675       CASE ( 'w"sa"' )
676          dopr_index(var_count) = 65
677          dopr_unit             = 'psu m/s'
678          hom(:,2,65,:) = SPREAD( zw, 2, statistic_regions+1 )
679          unit = dopr_unit
680
681       CASE ( 'w*sa*' )
682          dopr_index(var_count) = 66
683          dopr_unit             = 'psu m/s'
684          hom(:,2,66,:) = SPREAD( zw, 2, statistic_regions+1 )
685          unit = dopr_unit
686
687       CASE ( 'wsa' )
688          dopr_index(var_count) = 67
689          dopr_unit             = 'psu m/s'
690          hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
691          unit = dopr_unit
692
693       CASE DEFAULT
694          unit = 'illegal'
695
696    END SELECT
697
698
699 END SUBROUTINE ocean_check_data_output_pr
700
701
702!------------------------------------------------------------------------------!
703! Description:
704! ------------
705!> Define appropriate grid for netcdf variables.
706!> It is called out from subroutine netcdf.
707!------------------------------------------------------------------------------!
708 SUBROUTINE ocean_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
709   
710    IMPLICIT NONE
711
712    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< x grid of output variable
713    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !< y grid of output variable
714    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !< z grid of output variable
715    CHARACTER (LEN=*), INTENT(IN)  ::  var      !< name of output variable
716
717    LOGICAL, INTENT(OUT) ::  found   !< flag if output variable is found
718
719    found  = .TRUE.
720
721!
722!-- Check for the grid
723    SELECT CASE ( TRIM( var ) )
724
725       CASE ( 'rho_ocean', 'sa' )
726          grid_x = 'x'
727          grid_y = 'y'
728          grid_z = 'zu'
729
730       CASE DEFAULT
731          found  = .FALSE.
732          grid_x = 'none'
733          grid_y = 'none'
734          grid_z = 'none'
735
736    END SELECT
737
738 END SUBROUTINE ocean_define_netcdf_grid
739
740
741!------------------------------------------------------------------------------!
742! Description:
743! ------------
744!> Average 3D data.
745!------------------------------------------------------------------------------!
746 SUBROUTINE ocean_3d_data_averaging( mode, variable )
747 
748
749    USE arrays_3d,                                                             &
750        ONLY:  rho_ocean, sa
751
752    USE averaging,                                                             &
753        ONLY:  rho_ocean_av, sa_av
754
755    USE control_parameters,                                                    &
756        ONLY:  average_count_3d
757
758    USE indices,                                                               &
759        ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt
760
761    IMPLICIT NONE
762
763    CHARACTER (LEN=*) ::  mode       !< flag defining mode 'allocate', 'sum' or 'average'
764    CHARACTER (LEN=*) ::  variable   !< name of variable
765
766    INTEGER(iwp) ::  i   !< loop index
767    INTEGER(iwp) ::  j   !< loop index
768    INTEGER(iwp) ::  k   !< loop index
769
770    IF ( mode == 'allocate' )  THEN
771
772       SELECT CASE ( TRIM( variable ) )
773
774          CASE ( 'rho_ocean' )
775             IF ( .NOT. ALLOCATED( rho_ocean_av ) )  THEN
776                ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
777             ENDIF
778             rho_ocean_av = 0.0_wp
779
780          CASE ( 'sa' )
781             IF ( .NOT. ALLOCATED( sa_av ) )  THEN
782                ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
783             ENDIF
784             sa_av = 0.0_wp
785
786          CASE DEFAULT
787             CONTINUE
788
789       END SELECT
790
791    ELSEIF ( mode == 'sum' )  THEN
792
793       SELECT CASE ( TRIM( variable ) )
794
795          CASE ( 'rho_ocean' )
796             IF ( ALLOCATED( rho_ocean_av ) )  THEN
797                DO  i = nxlg, nxrg
798                   DO  j = nysg, nyng
799                      DO  k = nzb, nzt+1
800                         rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) +           &
801                                               rho_ocean(k,j,i)
802                      ENDDO
803                   ENDDO
804                ENDDO
805             ENDIF
806
807          CASE ( 'sa' )
808             IF ( ALLOCATED( sa_av ) )  THEN
809                DO  i = nxlg, nxrg
810                   DO  j = nysg, nyng
811                      DO  k = nzb, nzt+1
812                         sa_av(k,j,i) = sa_av(k,j,i) + sa(k,j,i)
813                      ENDDO
814                   ENDDO
815                ENDDO
816             ENDIF
817
818          CASE DEFAULT
819             CONTINUE
820
821       END SELECT
822
823    ELSEIF ( mode == 'average' )  THEN
824
825       SELECT CASE ( TRIM( variable ) )
826
827          CASE ( 'rho_ocean' )
828             IF ( ALLOCATED( rho_ocean_av ) ) THEN
829                DO  i = nxlg, nxrg
830                   DO  j = nysg, nyng
831                      DO  k = nzb, nzt+1
832                         rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) /           &
833                                               REAL( average_count_3d, KIND=wp )
834                      ENDDO
835                   ENDDO
836                ENDDO
837             ENDIF
838
839          CASE ( 'sa' )
840             IF ( ALLOCATED( sa_av ) ) THEN
841                DO  i = nxlg, nxrg
842                   DO  j = nysg, nyng
843                      DO  k = nzb, nzt+1
844                         sa_av(k,j,i) = sa_av(k,j,i) /                         &
845                                        REAL( average_count_3d, KIND=wp )
846                      ENDDO
847                   ENDDO
848                ENDDO
849             ENDIF
850
851       END SELECT
852
853    ENDIF
854
855 END SUBROUTINE ocean_3d_data_averaging
856
857
858!------------------------------------------------------------------------------!
859! Description:
860! ------------
861!> Define 2D output variables.
862!------------------------------------------------------------------------------!
863 SUBROUTINE ocean_data_output_2d( av, variable, found, grid, mode, local_pf,   &
864                                  nzb_do, nzt_do )
865 
866    USE arrays_3d,                                                             &
867        ONLY:  rho_ocean, sa
868
869    USE averaging,                                                             &
870        ONLY:  rho_ocean_av, sa_av
871
872    USE indices,                                                               &
873        ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
874
875    IMPLICIT NONE
876
877    CHARACTER (LEN=*) ::  grid       !< name of vertical grid
878    CHARACTER (LEN=*) ::  mode       !< either 'xy', 'xz' or 'yz'
879    CHARACTER (LEN=*) ::  variable   !< name of variable
880
881    INTEGER(iwp) ::  av        !< flag for (non-)average output
882    INTEGER(iwp) ::  flag_nr   !< number of masking flag
883    INTEGER(iwp) ::  i         !< loop index
884    INTEGER(iwp) ::  j         !< loop index
885    INTEGER(iwp) ::  k         !< loop index
886    INTEGER(iwp) ::  nzb_do    !< vertical output index (bottom)
887    INTEGER(iwp) ::  nzt_do    !< vertical output index (top)
888
889    LOGICAL ::  found   !< flag if output variable is found
890    LOGICAL ::  resorted  !< flag if output is already resorted
891
892    REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
893
894    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< local
895       !< array to which output data is resorted to
896
897    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
898   
899    found = .TRUE.
900    resorted = .FALSE.
901!
902!-- Set masking flag for topography for not resorted arrays
903    flag_nr = 0
904
905    SELECT CASE ( TRIM( variable ) )
906
907       CASE ( 'rho_ocean_xy', 'rho_ocean_xz', 'rho_ocean_yz' )
908          IF ( av == 0 )  THEN
909             to_be_resorted => rho_ocean
910          ELSE
911             IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
912                ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
913                rho_ocean_av = REAL( fill_value, KIND = wp )
914             ENDIF
915             to_be_resorted => rho_ocean_av
916          ENDIF
917
918       CASE ( 'sa_xy', 'sa_xz', 'sa_yz' )
919          IF ( av == 0 )  THEN
920             to_be_resorted => sa
921          ELSE
922             IF ( .NOT. ALLOCATED( sa_av ) ) THEN
923                ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
924                sa_av = REAL( fill_value, KIND = wp )
925             ENDIF
926             to_be_resorted => sa_av
927          ENDIF
928          IF ( mode == 'xy' ) grid = 'zu'
929
930       CASE DEFAULT
931          found = .FALSE.
932          grid  = 'none'
933
934    END SELECT
935
936    IF ( found .AND. .NOT. resorted )  THEN
937!       DO  i = nxl, nxr
938!          DO  j = nys, nyn
939!             DO  k = nzb_do, nzt_do
940!                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
941!                                         REAL( fill_value, KIND = wp ),        &
942!                                         BTEST( wall_flags_0(k,j,i), flag_nr ) )
943!             ENDDO
944!          ENDDO
945!       ENDDO
946    ENDIF
947 
948 END SUBROUTINE ocean_data_output_2d
949
950 
951!------------------------------------------------------------------------------!
952! Description:
953! ------------
954!> Define 3D output variables.
955!------------------------------------------------------------------------------!
956 SUBROUTINE ocean_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
957 
958
959    USE arrays_3d,                                                             &
960        ONLY:  rho_ocean, sa
961
962    USE averaging,                                                             &
963        ONLY:  rho_ocean_av, sa_av
964
965    USE indices,                                                               &
966        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,           &
967               wall_flags_0
968
969    IMPLICIT NONE
970
971    CHARACTER (LEN=*) ::  variable   !< name of variable
972
973    INTEGER(iwp) ::  av        !< flag for (non-)average output
974    INTEGER(iwp) ::  flag_nr   !< number of masking flag
975    INTEGER(iwp) ::  i         !< loop index
976    INTEGER(iwp) ::  j         !< loop index
977    INTEGER(iwp) ::  k         !< loop index
978    INTEGER(iwp) ::  nzb_do    !< lower limit of the data output (usually 0)
979    INTEGER(iwp) ::  nzt_do    !< vertical upper limit of the data output (usually nz_do3d)
980
981    LOGICAL ::  found     !< flag if output variable is found
982    LOGICAL ::  resorted  !< flag if output is already resorted
983
984    REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
985
986    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< local
987                                  !< array to which output data is resorted to
988
989    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
990
991    found = .TRUE.
992    resorted = .FALSE.
993!
994!-- Set masking flag for topography for not resorted arrays
995    flag_nr = 0
996
997    SELECT CASE ( TRIM( variable ) )
998
999       CASE ( 'rho_ocean' )
1000          IF ( av == 0 )  THEN
1001             to_be_resorted => rho_ocean
1002          ELSE
1003             IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
1004                ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1005                rho_ocean_av = REAL( fill_value, KIND = wp )
1006             ENDIF
1007             to_be_resorted => rho_ocean_av
1008          ENDIF
1009
1010       CASE ( 'sa' )
1011          IF ( av == 0 )  THEN
1012             to_be_resorted => sa
1013          ELSE
1014             IF ( .NOT. ALLOCATED( sa_av ) ) THEN
1015                ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1016                sa_av = REAL( fill_value, KIND = wp )
1017             ENDIF
1018             to_be_resorted => sa_av
1019          ENDIF
1020
1021       CASE DEFAULT
1022          found = .FALSE.
1023
1024    END SELECT
1025
1026
1027    IF ( found  .AND.  .NOT. resorted )  THEN
1028       DO  i = nxl, nxr
1029          DO  j = nys, nyn
1030             DO  k = nzb_do, nzt_do
1031                local_pf(i,j,k) = MERGE(                                       &
1032                                   to_be_resorted(k,j,i),                      &
1033                                   REAL( fill_value, KIND = wp ),              &
1034                                   BTEST( wall_flags_0(k,j,i), flag_nr ) )
1035             ENDDO
1036          ENDDO
1037       ENDDO
1038       resorted = .TRUE.
1039    ENDIF
1040
1041 END SUBROUTINE ocean_data_output_3d
1042
1043
1044!------------------------------------------------------------------------------!
1045! Description:
1046! ------------
1047!> Allocate arrays and assign pointers.
1048!------------------------------------------------------------------------------!
1049 SUBROUTINE ocean_init_arrays
1050
1051    USE indices,                                                               &
1052        ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt
1053
1054    USE pmc_interface,                                                         &
1055        ONLY:  nested_run
1056
1057    IMPLICIT NONE
1058
1059#if defined( __nopointer )
1060    ALLOCATE( prho(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
1061              rho_ocean(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
1062              sa(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                            &
1063              sa_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
1064              tsa_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1065#else
1066    ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                        &
1067              rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                         &
1068              sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
1069              sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
1070              sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1071
1072    prho => prho_1
1073    rho_ocean  => rho_1  ! routines calc_mean_profile and diffusion_e require
1074                         ! density to be a pointer
1075#endif
1076
1077#if ! defined( __nopointer )
1078!
1079!-- Initial assignment of pointers
1080    sa  => sa_1;   sa_p  => sa_2;   tsa_m  => sa_3
1081#endif
1082
1083 END SUBROUTINE ocean_init_arrays
1084
1085
1086!------------------------------------------------------------------------------!
1087! Description:
1088! ------------
1089!> Initialization of quantities needed for the ocean mode
1090!------------------------------------------------------------------------------!
1091 SUBROUTINE ocean_init
1092
1093
1094    USE arrays_3d,                                                             &
1095        ONLY:  dzu, hyp, pt_init, ref_state, zu, zw
1096
1097    USE basic_constants_and_equations_mod,                                     &
1098        ONLY:  g
1099
1100    USE control_parameters,                                                    &
1101        ONLY:  initializing_actions, molecular_viscosity, rho_surface,         &
1102               rho_reference, surface_pressure, use_single_reference_value
1103
1104    USE indices,                                                               &
1105        ONLY:  nxl, nxlg, nxrg, nyng, nys, nysg, nzb, nzt
1106
1107    USE kinds
1108
1109    USE pegrid
1110
1111    USE statistics,                                                            &
1112        ONLY:  hom, statistic_regions
1113
1114    IMPLICIT NONE
1115
1116    INTEGER(iwp) ::  i  !< loop index
1117    INTEGER(iwp) ::  j  !< loop index
1118    INTEGER(iwp) ::  k  !< loop index
1119    INTEGER(iwp) ::  n  !< loop index
1120
1121    REAL(wp)     ::  dum   !< dummy argument
1122    REAL(wp)     ::  pt_l  !< local scalar for pt used in equation of state function
1123    REAL(wp)     ::  sa_l  !< local scalar for sa used in equation of state function
1124
1125    REAL(wp), DIMENSION(nzb:nzt+1) ::  rho_ocean_init  !< local array for initial density
1126
1127    ALLOCATE( hyp(nzb:nzt+1) )
1128
1129
1130!
1131!-- In case of no restart run, calculate the inital salinity profilevcusing the
1132!-- given salinity gradients
1133    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1134
1135       sa_init = sa_surface
1136!
1137!--    Last arguments gives back the gradient at top level to be used as
1138!--    possible Neumann boundary condition. This is not realized for the ocean
1139!--    mode, therefore a dummy argument is used.
1140       CALL init_vertical_profiles( sa_vertical_gradient_level_ind,          &
1141                                    sa_vertical_gradient_level,              &
1142                                    sa_vertical_gradient, sa_init,           &
1143                                    sa_surface, dum )
1144    ENDIF
1145
1146!
1147!-- Initialize required 3d-arrays
1148    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
1149         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1150!
1151!--    Initialization via computed 1D-model profiles
1152       IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 )  THEN
1153
1154          DO  i = nxlg, nxrg
1155             DO  j = nysg, nyng
1156                sa(:,j,i) = sa_init
1157             ENDDO
1158          ENDDO
1159
1160       ENDIF
1161!
1162!--    Store initial profiles for output purposes etc.
1163!--    Store initial salinity profile
1164       hom(:,1,26,:)  = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
1165!
1166!--    Initialize old and new time levels.
1167       tsa_m = 0.0_wp
1168       sa_p  = sa
1169
1170    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
1171
1172!
1173!--    Initialize new time levels (only done in order to set boundary values
1174!--    including ghost points)
1175       sa_p = sa
1176!
1177!--    Allthough tendency arrays are set in prognostic_equations, they have
1178!--    have to be predefined here because they are used (but multiplied with 0)
1179!--    there before they are set.
1180       tsa_m = 0.0_wp
1181
1182    ENDIF
1183
1184!
1185!-- Set water density near the ocean surface
1186    rho_surface = 1027.62_wp
1187
1188!
1189!-- Set kinematic viscosity to sea water at 20C.
1190!-- This changes the default value that is given for air!
1191    molecular_viscosity = 1.05E-6_wp
1192
1193!
1194!-- Change sign of buoyancy/stability terms because density gradient is used
1195!-- instead of the potential temperature gradient to calculate the buoyancy
1196    atmos_ocean_sign = -1.0_wp
1197
1198!
1199!-- Calculate initial vertical profile of hydrostatic pressure (in Pa)
1200!-- and the reference density (used later in buoyancy term)
1201!-- First step: Calculate pressure using reference density
1202    hyp(nzt+1) = surface_pressure * 100.0_wp
1203    hyp(nzt)   = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1)
1204    rho_ocean_init(nzt)   = rho_surface
1205    rho_ocean_init(nzt+1) = rho_surface  ! only required for output
1206
1207    DO  k = nzt-1, 1, -1
1208       hyp(k) = hyp(k+1) + rho_surface * g * dzu(k)
1209    ENDDO
1210    hyp(0) = hyp(1) + rho_surface * g * dzu(1)
1211
1212!
1213!-- Second step: Iteratively calculate in situ density (based on presssure)
1214!-- and pressure (based on in situ density)
1215    DO  n = 1, 5
1216
1217       rho_reference = rho_surface * 0.5_wp * dzu(nzt+1)
1218
1219       DO  k = nzt, 0, -1
1220
1221          sa_l = 0.5_wp * ( sa_init(k) + sa_init(k+1) )
1222          pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) )
1223
1224          rho_ocean_init(k) = eqn_state_seawater_func( hyp(k), pt_l, sa_l )
1225
1226          rho_reference = rho_reference + rho_ocean_init(k) * dzu(k+1)
1227
1228       ENDDO
1229
1230       rho_reference = rho_reference / ( zw(nzt) - zu(nzb) )
1231
1232       hyp(nzt) = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1)
1233       DO  k = nzt-1, 0, -1
1234          hyp(k) = hyp(k+1) + g * 0.5_wp * ( rho_ocean_init(k)                 &
1235                                           + rho_ocean_init(k+1) ) * dzu(k+1)
1236       ENDDO
1237
1238    ENDDO
1239
1240!
1241!-- Calculate the reference potential density
1242    prho_reference = 0.0_wp
1243    DO  k = 0, nzt
1244
1245       sa_l = 0.5_wp * ( sa_init(k) + sa_init(k+1) )
1246       pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) )
1247
1248       prho_reference = prho_reference + dzu(k+1) * &
1249                        eqn_state_seawater_func( 0.0_wp, pt_l, sa_l )
1250
1251    ENDDO
1252
1253    prho_reference = prho_reference / ( zu(nzt) - zu(nzb) )
1254
1255!
1256!-- Calculate the 3d array of initial in situ and potential density,
1257!-- based on the initial temperature and salinity profile
1258    CALL eqn_state_seawater
1259
1260!
1261!-- Store initial density profile
1262    hom(:,1,77,:)  = SPREAD( rho_ocean_init(:), 2, statistic_regions+1 )
1263
1264!
1265!-- Set the reference state to be used in the buoyancy terms
1266    IF ( use_single_reference_value )  THEN
1267       ref_state(:) = rho_reference
1268    ELSE
1269       ref_state(:) = rho_ocean_init(:)
1270    ENDIF
1271
1272
1273 END SUBROUTINE ocean_init
1274
1275
1276!------------------------------------------------------------------------------!
1277! Description:
1278! ------------
1279!> Prognostic equation for salinity.
1280!> Vector-optimized version
1281!------------------------------------------------------------------------------!
1282 SUBROUTINE ocean_prognostic_equations
1283
1284    USE advec_s_bc_mod,                                                        &
1285        ONLY:  advec_s_bc
1286
1287    USE advec_s_pw_mod,                                                        &
1288        ONLY:  advec_s_pw
1289
1290    USE advec_s_up_mod,                                                        &
1291        ONLY:  advec_s_up
1292
1293    USE advec_ws,                                                              &
1294        ONLY:  advec_s_ws
1295
1296    USE arrays_3d,                                                             &
1297        ONLY:  rdf_sc, tend, tsa_m
1298
1299    USE control_parameters,                                                    &
1300        ONLY:  dt_3d, intermediate_timestep_count, intermediate_timestep_count_max,   &
1301               scalar_advec, timestep_scheme, tsc, ws_scheme_sca
1302
1303    USE cpulog,                                                                &
1304        ONLY:  cpu_log, log_point
1305
1306    USE diffusion_s_mod,                                                       &
1307        ONLY:  diffusion_s
1308
1309    USE indices,                                                               &
1310        ONLY:  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
1311
1312    USE surface_mod,                                                           &
1313        ONLY:  surf_def_v, surf_def_h, surf_lsm_h, surf_lsm_v, surf_usm_h,     &
1314               surf_usm_v
1315
1316    USE user_actions_mod,                                                      &
1317        ONLY:  user_actions
1318
1319    IMPLICIT NONE
1320
1321    INTEGER(iwp) ::  i       !< loop index
1322    INTEGER(iwp) ::  j       !< loop index
1323    INTEGER(iwp) ::  k       !< loop index
1324
1325    REAL(wp)     ::  sbt     !< weighting factor for sub-time step
1326
1327!
1328!-- Compute prognostic equations for the ocean mode
1329!-- First, start with salinity
1330    CALL cpu_log( log_point(37), 'sa-equation', 'start' )
1331
1332!
1333!-- sa-tendency terms with communication
1334    sbt = tsc(2)
1335    IF ( scalar_advec == 'bc-scheme' )  THEN
1336
1337       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1338!
1339!--       Bott-Chlond scheme always uses Euler time step. Thus:
1340          sbt = 1.0_wp
1341       ENDIF
1342       tend = 0.0_wp
1343       CALL advec_s_bc( sa, 'sa' )
1344
1345    ENDIF
1346
1347!
1348!-- sa-tendency terms with no communication
1349    IF ( scalar_advec /= 'bc-scheme' )  THEN
1350       tend = 0.0_wp
1351       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1352          IF ( ws_scheme_sca )  THEN
1353             CALL advec_s_ws( sa, 'sa' )
1354          ELSE
1355             CALL advec_s_pw( sa )
1356          ENDIF
1357       ELSE
1358          CALL advec_s_up( sa )
1359       ENDIF
1360    ENDIF
1361
1362    CALL diffusion_s( sa,                                                      &
1363                      surf_def_h(0)%sasws, surf_def_h(1)%sasws,                &
1364                      surf_def_h(2)%sasws,                                     &
1365                      surf_lsm_h%sasws,    surf_usm_h%sasws,                   &
1366                      surf_def_v(0)%sasws, surf_def_v(1)%sasws,                &
1367                      surf_def_v(2)%sasws, surf_def_v(3)%sasws,                &
1368                      surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws,                &
1369                      surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws,                &
1370                      surf_usm_v(0)%sasws, surf_usm_v(1)%sasws,                &
1371                      surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
1372
1373    CALL user_actions( 'sa-tendency' )
1374
1375!
1376!-- Prognostic equation for salinity
1377    DO  i = nxl, nxr
1378       DO  j = nys, nyn
1379          DO  k = nzb+1, nzt
1380             sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +         &
1381                                                   tsc(3) * tsa_m(k,j,i) )     &
1382                                               - tsc(5) * rdf_sc(k) *          &
1383                                                 ( sa(k,j,i) - sa_init(k) )    &
1384                                       )                                       &
1385                                         * MERGE( 1.0_wp, 0.0_wp,              &
1386                                               BTEST( wall_flags_0(k,j,i), 0 ) &
1387                                                )
1388             IF ( sa_p(k,j,i) < 0.0_wp )  sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
1389          ENDDO
1390       ENDDO
1391    ENDDO
1392
1393!
1394!-- Calculate tendencies for the next Runge-Kutta step
1395    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1396       IF ( intermediate_timestep_count == 1 )  THEN
1397          DO  i = nxl, nxr
1398             DO  j = nys, nyn
1399                DO  k = nzb+1, nzt
1400                   tsa_m(k,j,i) = tend(k,j,i)
1401                ENDDO
1402             ENDDO
1403          ENDDO
1404       ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )&
1405       THEN
1406          DO  i = nxl, nxr
1407             DO  j = nys, nyn
1408                DO  k = nzb+1, nzt
1409                   tsa_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
1410                                     5.3125_wp * tsa_m(k,j,i)
1411                ENDDO
1412             ENDDO
1413         ENDDO
1414       ENDIF
1415    ENDIF
1416
1417    CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
1418
1419!
1420!-- Calculate density by the equation of state for seawater
1421    CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
1422    CALL eqn_state_seawater
1423    CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
1424
1425 END SUBROUTINE ocean_prognostic_equations
1426
1427
1428!------------------------------------------------------------------------------!
1429! Description:
1430! ------------
1431!> Prognostic equations for ocean mode (so far, salinity only)
1432!> Cache-optimized version
1433!------------------------------------------------------------------------------!
1434 SUBROUTINE ocean_prognostic_equations_ij( i, j, i_omp_start, tn )
1435
1436    USE advec_s_pw_mod,                                                        &
1437        ONLY:  advec_s_pw
1438
1439    USE advec_s_up_mod,                                                        &
1440        ONLY:  advec_s_up
1441
1442    USE advec_ws,                                                              &
1443        ONLY:  advec_s_ws
1444
1445    USE arrays_3d,                                                             &
1446        ONLY:  diss_l_sa, diss_s_sa, flux_l_sa, flux_s_sa, rdf_sc, tend, tsa_m
1447
1448    USE control_parameters,                                                    &
1449        ONLY:  dt_3d, intermediate_timestep_count,                             &
1450               intermediate_timestep_count_max, timestep_scheme, tsc,          &
1451               ws_scheme_sca
1452
1453    USE diffusion_s_mod,                                                       &
1454        ONLY:  diffusion_s
1455
1456    USE indices,                                                               &
1457        ONLY:  nzb, nzt, wall_flags_0
1458
1459    USE surface_mod,                                                           &
1460        ONLY:  surf_def_v, surf_def_h, surf_lsm_h, surf_lsm_v, surf_usm_h,     &
1461               surf_usm_v
1462
1463    USE user_actions_mod,                                                      &
1464        ONLY:  user_actions
1465
1466    IMPLICIT NONE
1467
1468    INTEGER(iwp) ::  i             !< loop index x direction
1469    INTEGER(iwp) ::  i_omp_start   !< first loop index of i-loop in calling    &
1470                                   !< routine prognostic_equations
1471    INTEGER(iwp) ::  j             !< loop index y direction
1472    INTEGER(iwp) ::  k             !< loop index z direction
1473    INTEGER(iwp) ::  tn            !< task number of openmp task
1474
1475!
1476!-- Compute prognostic equations for the ocean mode
1477!-- First, start with tendency-terms for salinity
1478    tend(:,j,i) = 0.0_wp
1479    IF ( timestep_scheme(1:5) == 'runge' ) &
1480    THEN
1481       IF ( ws_scheme_sca )  THEN
1482          CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa,  diss_s_sa, flux_l_sa,   &
1483                           diss_l_sa, i_omp_start, tn  )
1484       ELSE
1485          CALL advec_s_pw( i, j, sa )
1486       ENDIF
1487    ELSE
1488       CALL advec_s_up( i, j, sa )
1489    ENDIF
1490    CALL diffusion_s( i, j, sa,                                       &
1491                      surf_def_h(0)%sasws, surf_def_h(1)%sasws,       &
1492                      surf_def_h(2)%sasws,                            &
1493                      surf_lsm_h%sasws,    surf_usm_h%sasws,          &
1494                      surf_def_v(0)%sasws, surf_def_v(1)%sasws,       &
1495                      surf_def_v(2)%sasws, surf_def_v(3)%sasws,       &
1496                      surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws,       &
1497                      surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws,       &
1498                      surf_usm_v(0)%sasws, surf_usm_v(1)%sasws,       &
1499                      surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
1500
1501    CALL user_actions( i, j, 'sa-tendency' )
1502
1503!
1504!-- Prognostic equation for salinity
1505    DO  k = nzb+1, nzt
1506
1507       sa_p(k,j,i) = sa(k,j,i) + ( dt_3d *                                     &
1508                                           ( tsc(2) * tend(k,j,i) +            &
1509                                             tsc(3) * tsa_m(k,j,i) )           &
1510                                 - tsc(5) * rdf_sc(k)                          &
1511                                          * ( sa(k,j,i) - sa_init(k) )         &
1512                                 ) * MERGE( 1.0_wp, 0.0_wp,                    &
1513                                            BTEST( wall_flags_0(k,j,i), 0 ) )
1514
1515       IF ( sa_p(k,j,i) < 0.0_wp )  sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
1516
1517    ENDDO
1518
1519!
1520!-- Calculate tendencies for the next Runge-Kutta step
1521    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1522       IF ( intermediate_timestep_count == 1 )  THEN
1523          DO  k = nzb+1, nzt
1524             tsa_m(k,j,i) = tend(k,j,i)
1525          ENDDO
1526       ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )&
1527       THEN
1528          DO  k = nzb+1, nzt
1529             tsa_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
1530                               5.3125_wp * tsa_m(k,j,i)
1531          ENDDO
1532       ENDIF
1533    ENDIF
1534
1535!
1536!-- Calculate density by the equation of state for seawater
1537    CALL eqn_state_seawater( i, j )
1538
1539 END SUBROUTINE ocean_prognostic_equations_ij
1540
1541
1542!------------------------------------------------------------------------------!
1543! Description:
1544! ------------
1545!> Swapping of timelevels.
1546!------------------------------------------------------------------------------!
1547 SUBROUTINE ocean_swap_timelevel( mod_count )
1548
1549    IMPLICIT NONE
1550
1551    INTEGER, INTENT(IN) ::  mod_count  !< flag defining where pointers point to
1552
1553#if defined( __nopointer )
1554
1555    sa = sa_p
1556
1557#else
1558   
1559    SELECT CASE ( mod_count )
1560
1561       CASE ( 0 )
1562          sa => sa_1;    sa_p => sa_2
1563
1564       CASE ( 1 )
1565          sa => sa_2;    sa_p => sa_1
1566
1567    END SELECT
1568
1569#endif
1570
1571 END SUBROUTINE ocean_swap_timelevel
1572
1573
1574!------------------------------------------------------------------------------!
1575! Description:
1576! ------------
1577!> This routine reads the respective restart data for the ocean module.
1578!------------------------------------------------------------------------------!
1579 SUBROUTINE ocean_rrd_global( found )
1580
1581
1582    USE control_parameters,                                                    &
1583        ONLY: length, restart_string
1584
1585
1586    IMPLICIT NONE
1587
1588    LOGICAL, INTENT(OUT)  ::  found
1589
1590
1591    found = .TRUE.
1592
1593    SELECT CASE ( restart_string(1:length) )
1594
1595       CASE ( 'bc_sa_t' )
1596          READ ( 13 )  bc_sa_t
1597
1598       CASE ( 'bottom_salinityflux' )
1599          READ ( 13 )  bottom_salinityflux
1600
1601       CASE ( 'langmuir_number' )
1602          READ ( 13 )  langmuir_number
1603
1604       CASE ( 'sa_init' )
1605          READ ( 13 )  sa_init
1606
1607       CASE ( 'sa_surface' )
1608          READ ( 13 )  sa_surface
1609
1610       CASE ( 'sa_vertical_gradient' )
1611          READ ( 13 )  sa_vertical_gradient
1612
1613       CASE ( 'sa_vertical_gradient_level' )
1614          READ ( 13 )  sa_vertical_gradient_level
1615
1616       CASE ( 'stokes_force' )
1617          READ ( 13 )  stokes_force
1618
1619       CASE ( 'stokes_waveheight' )
1620          READ ( 13 )  stokes_waveheight
1621
1622       CASE ( 'stokes_wavelength' )
1623          READ ( 13 )  stokes_wavelength
1624
1625       CASE ( 'top_salinityflux' )
1626          READ ( 13 )  top_salinityflux
1627
1628       CASE ( 'wall_salinityflux' )
1629          READ ( 13 )  wall_salinityflux
1630
1631       CASE DEFAULT
1632
1633          found = .FALSE.
1634
1635    END SELECT
1636
1637 END SUBROUTINE ocean_rrd_global
1638
1639
1640!------------------------------------------------------------------------------!
1641! Description:
1642! ------------
1643!> This routine reads the respective restart data for the ocean module.
1644!------------------------------------------------------------------------------!
1645 SUBROUTINE ocean_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,        &
1646                             nxr_on_file, nynf, nync, nyn_on_file, nysf,       &
1647                             nysc, nys_on_file, tmp_2d, tmp_3d, found )
1648
1649    USE averaging,                                                             &
1650        ONLY:  rho_ocean_av, sa_av
1651
1652    USE control_parameters,                                                    &
1653        ONLY:  length, restart_string
1654
1655    USE indices,                                                               &
1656        ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzt
1657
1658    USE pegrid
1659
1660
1661    IMPLICIT NONE
1662
1663    INTEGER(iwp) ::  i               !<
1664    INTEGER(iwp) ::  k               !<
1665    INTEGER(iwp) ::  nxlc            !<
1666    INTEGER(iwp) ::  nxlf            !<
1667    INTEGER(iwp) ::  nxl_on_file     !<
1668    INTEGER(iwp) ::  nxrc            !<
1669    INTEGER(iwp) ::  nxrf            !<
1670    INTEGER(iwp) ::  nxr_on_file     !<
1671    INTEGER(iwp) ::  nync            !<
1672    INTEGER(iwp) ::  nynf            !<
1673    INTEGER(iwp) ::  nyn_on_file     !<
1674    INTEGER(iwp) ::  nysc            !<
1675    INTEGER(iwp) ::  nysf            !<
1676    INTEGER(iwp) ::  nys_on_file     !<
1677
1678    LOGICAL, INTENT(OUT)  ::  found
1679
1680    REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d   !<
1681    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
1682
1683
1684    found = .TRUE.
1685
1686    SELECT CASE ( restart_string(1:length) )
1687
1688       CASE ( 'rho_ocean_av' )
1689          IF ( .NOT. ALLOCATED( rho_ocean_av ) )  THEN
1690             ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1691          ENDIF
1692          IF ( k == 1 )  READ ( 13 )  tmp_3d
1693          rho_ocean_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =            &
1694                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1695
1696       CASE ( 'sa' )
1697          IF ( k == 1 )  READ ( 13 )  tmp_3d
1698          sa(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                      &
1699                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1700
1701       CASE ( 'sa_av' )
1702          IF ( .NOT. ALLOCATED( sa_av ) )  THEN
1703              ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1704          ENDIF
1705          IF ( k == 1 )  READ ( 13 )  tmp_3d
1706          sa_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                   &
1707                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1708
1709       CASE DEFAULT
1710          found = .FALSE.
1711
1712    END SELECT
1713
1714 END SUBROUTINE ocean_rrd_local
1715
1716
1717!------------------------------------------------------------------------------!
1718! Description:
1719! ------------
1720!> This routine writes the respective restart data for the ocean module.
1721!------------------------------------------------------------------------------!
1722 SUBROUTINE ocean_wrd_global
1723
1724
1725    IMPLICIT NONE
1726
1727    CALL wrd_write_string( 'bc_sa_t' )
1728    WRITE ( 14 )  bc_sa_t
1729
1730    CALL wrd_write_string( 'bottom_salinityflux' )
1731    WRITE ( 14 )  bottom_salinityflux
1732
1733    CALL wrd_write_string( 'langmuir_number' )
1734    WRITE ( 14 )  langmuir_number
1735
1736    CALL wrd_write_string( 'sa_init' )
1737    WRITE ( 14 )  sa_init
1738
1739    CALL wrd_write_string( 'sa_surface' )
1740    WRITE ( 14 )  sa_surface
1741
1742    CALL wrd_write_string( 'sa_vertical_gradient' )
1743    WRITE ( 14 )  sa_vertical_gradient
1744
1745    CALL wrd_write_string( 'sa_vertical_gradient_level' )
1746    WRITE ( 14 )  sa_vertical_gradient_level
1747
1748    CALL wrd_write_string( 'stokes_force' )
1749    WRITE ( 14 )  stokes_force
1750
1751    CALL wrd_write_string( 'stokes_waveheight' )
1752    WRITE ( 14 )  stokes_waveheight
1753
1754    CALL wrd_write_string( 'stokes_wavelength' )
1755    WRITE ( 14 )  stokes_wavelength
1756
1757    CALL wrd_write_string( 'top_salinityflux' )
1758    WRITE ( 14 )  top_salinityflux
1759
1760    CALL wrd_write_string( 'wall_salinityflux' )
1761    WRITE ( 14 )  wall_salinityflux
1762
1763 END SUBROUTINE ocean_wrd_global
1764
1765
1766!------------------------------------------------------------------------------!
1767! Description:
1768! ------------
1769!> This routine writes the respective restart data for the ocean module.
1770!------------------------------------------------------------------------------!
1771 SUBROUTINE ocean_wrd_local
1772
1773    USE averaging,                                                             &
1774        ONLY:  rho_ocean_av, sa_av
1775
1776    IMPLICIT NONE
1777
1778    IF ( ALLOCATED( rho_ocean_av ) )  THEN
1779       CALL wrd_write_string( 'rho_ocean_av' )
1780       WRITE ( 14 )  rho_ocean_av
1781    ENDIF
1782
1783    CALL wrd_write_string( 'sa' )
1784    WRITE ( 14 )  sa
1785
1786    IF ( ALLOCATED( sa_av ) )  THEN
1787       CALL wrd_write_string( 'sa_av' )
1788       WRITE ( 14 )  sa_av
1789    ENDIF
1790
1791 END SUBROUTINE ocean_wrd_local
1792
1793
1794 END MODULE ocean_mod
Note: See TracBrowser for help on using the repository browser.